// determine TOF efficiency as a function of x-y position at the location
// of the TOF detector. the input root file is hists_runXXXXXX.root


void caleff(int R, float MCut, int MCdata){

  int DEBUG = 0;

  float StatCut = 100;

  float TOFMatchCut = MCut; // +/- TOFMatchCUT to find tof point for track

  char rfile[128];
  sprintf(rfile,"EFF/hists_run%06d.root",R);

  if (MCdata){
    sprintf(rfile,"EFF/hists_MC_run%06d.root",R);
  }

  TFile *RF = new TFile(rfile, "READ");

  TH2F *xy = (TH2F*)RF->Get("xy");

  char hnam[128];

  float EfficienciesX[82][82];
  float EfficienciesY[82][82];
  float BinPosX[82][82];
  float BinPosY[82][82];


  char s[2][10] = { "x", "y" };
  char outf[128];
  sprintf(outf,"EFF/position_efficiency_%02d_run%06d.dat",(int)(TOFMatchCut*10.), R);
  if (MCdata){
    sprintf(outf,"EFF/position_efficiency_MC_%02d_run%06d.dat",(int)(TOFMatchCut*10.), R);
  }
  ofstream OutF(outf);

  for (int j=0; j<82; j++){

    float x = j*3.-123. + 1.5;

    for (int k=0; k<82; k++){

      float y = k*3.-123 + 1.5;

      EfficienciesX[j][k] = 0.;
      EfficienciesY[j][k] = 0.;
      BinPosX[j][k] = x;
      BinPosY[j][k] = y;


      //cout<<j+1<<"  "<<k+1<<endl;
      float Ntot = xy->GetBinContent(j+1, k+1);

      for (int iloc = 0; iloc<2; iloc++) {

	sprintf(hnam, "xy_p_vs_d%s%02d%02d", s[iloc], j, k);
	
	TH2F *h = (TH2F*)RF->Get(hnam);

	if (Ntot<100.){
	  OutF<<setw(2)<<j<<" "<<setw(2)<<k;
	  OutF<<"   "<<x<<"  "<<y<<"     "<<iloc;
	  OutF<<"   "<<"0 "<<"   "<<"0 "<<"     "<<"0"<<"     "<<"0"<<endl;
	  continue;
	}
	
	int nbins = h->GetYaxis()->GetNbins();
	TH1D *hproj = h->ProjectionX("hproj", 3,nbins+1); // bin 3 is >0.5GeV
	int nbinsp = hproj->GetXaxis()->GetNbins();
	
	float Integral = hproj->Integral(1,nbinsp);
	float AllHits = hproj->Integral(0,nbinsp+1);
	
	if (Integral<StatCut){
	  OutF<<setw(2)<<j<<" "<<setw(2)<<k;
	  OutF<<"   "<<x<<"  "<<y<<"     "<<iloc;
	  OutF<<"   "<<"0 "<<"   "<<"0 "<<"     "<<"0"<<"     "<<"0"<<endl;
	  continue;
	}
	
	char ptit[128];
	sprintf(ptit, "(%6.1f / %6.1f) projection xbin%02d ybin%02d", x, y, j, k);
	hproj->SetTitle(ptit);
	sprintf(ptit,"#Delta %s [cm]", s[iloc]);
	hproj->GetXaxis()->SetTitle(ptit);
	
	int bp = hproj->GetMaximumBin();
	float p = hproj->GetBinCenter(bp);
	
	if (DEBUG) {
	  h->Draw("colz");
	  gPad->SetGrid();
	  gPad->Update();
	  getchar();
	}
	hproj->Draw(); 
	hproj->Fit("gaus", "Q", "R", p-4., p+4.);
	TF1 *f1 = hproj->GetFunction("gaus");
	float pos = f1->GetParameter(1);
	float sig = f1->GetParameter(2);
	
	float low = pos - 4.;
	float hig = pos + 4.;
	int binlow = hproj->FindBin(low);
	int binhig = hproj->FindBin(hig);

	binlow = hproj->FindBin(-TOFMatchCut);
	binhig = hproj->FindBin(TOFMatchCut);
	
	float Integral3sig = hproj->Integral(binlow, binhig);
	float MatchEfficiency = Integral3sig/Ntot;
	float dMatchEfficiency = TMath::Sqrt(( 
					      (Ntot-Integral3sig) + (Integral3sig*Integral3sig)/
					      Ntot )/Ntot/Ntot);

	if (iloc == 0){
	  EfficienciesX[j][k] = MatchEfficiency;
	} else{
	  EfficienciesY[j][k] = MatchEfficiency;
	}
	BinPosX[j][k] = x;
	BinPosY[j][k] = y;
	
	//cout<<setw(2)<<j<<" / "<<setw(2)<<k;
	//cout<<"      (x/y) = "<<x<<"/"<<y;
	//cout<<"     "<<s[iloc]<<" "<<MatchEfficiency<<"      "<<pos<<" "<<sig<<endl;
	
	OutF<<setw(2)<<j<<" "<<setw(2)<<k;
	OutF<<"   "<<x<<"  "<<y<<"     "<<iloc;
	OutF<<"   "<<pos<<"   "<<sig<<"     "<<MatchEfficiency<<"     "<<dMatchEfficiency<<endl;
	
	if (DEBUG) {
	  gPad->SetGrid();
	  gPad->Update();
	  getchar();
	}
      }
    }
  }

  OutF.close();
  

}