[u/mrichter/AliRoot.git] / PWG4 / totEt / macros / hadEt / CorrEfficiency.C

//Christine Nattrass, University of Tennessee at Knoxville
//This macro is for calculating the single track efficiency for the hadronic transverse energy measurement
//Uses the output of AliAnalysisTaskHadEt
//This is not actually what gets used in the correction class AliAnalysisHadEtCorrections - that is done in the macro GetCorrections.C - but this is useful for making plots and playing around with different options

//this function calculates the efficiency propagating errors properly
TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name) 
{
    if(!numerator){
      cerr<<"Error:  numerator does not exist!"<<endl;
      return NULL;
    }
    if(!denominator){
      cerr<<"Error:  denominator does not exist!"<<endl;
      return NULL;
    }
    TH1D* result = (TH1D*)numerator->Clone(name);
    Int_t nbins = numerator->GetNbinsX();
    for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
      Double_t numeratorVal = numerator->GetBinContent(ibin);
      Double_t denominatorVal = denominator->GetBinContent(ibin);
      // Check if the errors are right or the thing is scaled
      Double_t numeratorValErr = numerator->GetBinError(ibin);
      if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
	numeratorVal /= rescale;
      }
      Double_t denominatorValErr = denominator->GetBinError(ibin);
      if (!(denominatorValErr==0. || denominatorVal==0. )) {
	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
	denominatorVal /= rescale;
      }
      Double_t quotient = 0.;
      if (denominatorVal!=0.) {
	quotient = numeratorVal/denominatorVal;
      }
      Double_t quotientErr=0;
      quotientErr = TMath::Sqrt(
				(numeratorVal+1.0)/(denominatorVal+2.0)*
				((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
      result->SetBinContent(ibin,quotient);
      result->SetBinError(ibin,quotientErr);
      //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
    }
    return result;
}


//This is a somewhat messy function that gets the efficiency for different particles
TH1D *GetHisto(float cut = 0.12, char *name, int mycase, bool eta, int color, int marker,bool TPC, bool ITS){
  //TFile *file = new TFile("Et.ESD.new.sim.merged.root");
  TFile *file = new TFile("Et.ESD.new.sim.LHC10d4.pp.merged.root");
  TList *list = file->FindObject("out2");
  char *myname = "ITS";
  if(TPC&&!ITS) myname = "TPC";
  if(TPC&&ITS) myname = "TPCITS";
  cout<<"Using tracks from "<<myname<<" for efficiency"<<endl;
  TH2F *numeratorParent; 
  switch(mycase){
  case 0:
    numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoHadron");
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")));
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")));
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
    break;
  case 1://pion
    numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoPion");
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
    break;
  case 2://kaon
    numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")))->Clone("RecoKaon");
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
    break;
  case 3://proton
    numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")))->Clone("RecoProton");
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
    break;
  case 4://electron
    numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EPlus")))->Clone("RecoElectron");
    numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EMinus")));
    break;
  }
  TH2F *denominatorParent; 
  switch(mycase){
  case 0:
    denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedChargedHadron"))->Clone("RecoHadron");
//     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoHadron");
//     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
//     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
//     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKPlus"));
//     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedProton"));
//     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
    break;
  case 1://pion
    denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoPion");
    denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
    break;
  case 2://kaon
    denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedKPlus"))->Clone("RecoKaon");
    denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
    break;
  case 3://proton
    denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedProton"))->Clone("RecoProton");
    denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
    break;
  case 4://electron
    denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedEPlus"))->Clone("RecoElectron");
    denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedEMinus"));
    break;
  }
  numeratorParent->Sumw2();
  denominatorParent->Sumw2();
  TH1D *denominator;
  TH1D *numerator;
  if(eta){
    int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
    int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
    cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
    denominator = denominatorParent->ProjectionX(Form("garbage%s",name),lowbin,highbin);
    numerator = numeratorParent->ProjectionX(name,lowbin,highbin);
  }
  else{
    int lowbin = denominatorParent->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
    int highbin = denominatorParent->GetXaxis()->GetNbins();
    cout<<"Here Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
    numerator = numeratorParent->ProjectionY(name,lowbin,highbin);
    denominator = denominatorParent->ProjectionY(Form("denominator%s",name),lowbin,highbin);
  }
  delete numeratorParent;
  delete denominatorParent;
  //numerator->Divide(denominator);
  TH1D *result = bayneseffdiv((TH1D*) numerator,(TH1D*)denominator,name);
  //result->Rebin(2);
  //result->Scale(0.5);
  result->SetYTitle("Efficiency");
  result->GetYaxis()->SetTitleOffset(0.8);
  result->GetXaxis()->SetTitleOffset(0.8);
  result->GetYaxis()->SetLabelSize(0.05);
  result->GetXaxis()->SetLabelSize(0.05);
  result->GetYaxis()->SetTitleSize(0.05);
  result->GetXaxis()->SetTitleSize(0.05);
  result->SetMarkerColor(color);
  result->SetLineColor(color);
  result->SetMarkerStyle(marker);
  //result->Draw("e");
  return result;

}


//this is a method that makes pretty plots
void CorrEfficiency(char *prodname= "LHC10d4", char *shortprodname = "LHC10d4 PYTHIA D6T 7 TeV p+p", bool TPC = true,bool ITS = true, bool eta = false){

  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  TCanvas *c = new TCanvas("c","c",600,400);
  c->SetTopMargin(0.02);
  c->SetRightMargin(0.02);
  c->SetBorderSize(0);
  c->SetFillColor(0);
  c->SetFillColor(0);
  c->SetBorderMode(0);
  c->SetFrameFillColor(0);
  c->SetFrameBorderMode(0);
  //c->SetLogx();

  int colortotal = 1;
  int colorpi = 2;
  int colork = 3;
  int colorp = 4;
  int phosmarker = 20;
  int emcalmarker = 24;
  float ptcut1 = 0.05;
  float ptcut2 = 0.1;
  float phoscut = 0.12;
  float emcalcut = 0.7;
  if(!eta){
    phoscut = 0.1;
    emcalcut = 0.15;
  }
  TH1D *PHOStotal = GetHisto(phoscut,"PHOStotal",0,eta,colortotal,phosmarker,TPC,ITS);
  TH1D *PHOSpi = GetHisto(phoscut,"PHOSpi",1,eta,colorpi,phosmarker,TPC,ITS);
  TH1D *PHOSp = GetHisto(phoscut,"PHOSp",2,eta,colork,phosmarker,TPC,ITS);
  TH1D *PHOSk = GetHisto(phoscut,"PHOSk",3,eta,colorp,phosmarker,TPC,ITS);
  if(eta) PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));
//if(ITS&&!TPC){PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));}
//else{PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.0),PHOStotal->GetXaxis()->FindBin(3.0));}
  PHOStotal->SetMinimum(0.0);
  PHOStotal->SetMaximum(1.0);
  //parameters[centbin][0]*exp(-pow(parameters[centbin][1]/pt,parameters[centbin][2]))
  TF1 *func = new TF1("func","[0]*exp(-pow([1]/x,[2]))");
  func->SetParameter(0,.9);
  func->SetParameter(1,.05);
  func->SetParLimits(1,1e-3,1);
  func->SetParameter(2,.1);
  //PHOStotal->Fit(func);
  PHOStotal->Draw();
  PHOSpi->Draw("same");
  PHOSp->Draw("same");
  PHOSk->Draw("same");
  TH1D *EMCALtotal = GetHisto(emcalcut,"EMCALtotal",0,eta,colortotal,emcalmarker,TPC,ITS);
  TH1D *EMCALpi = GetHisto(emcalcut,"EMCALpi",1,eta,colorpi,emcalmarker,TPC,ITS);
  TH1D *EMCALp = GetHisto(emcalcut,"EMCALp",2,eta,colork,emcalmarker,TPC,ITS);
  TH1D *EMCALk = GetHisto(emcalcut,"EMCALk",3,eta,colorp,emcalmarker,TPC,ITS);
  EMCALtotal->Draw("same");
  EMCALpi->Draw("same");
  EMCALp->Draw("same");
  EMCALk->Draw("same");


  TLegend *leg = new  TLegend(0.22651,0.247312,0.370805,0.438172);
  leg->AddEntry(PHOStotal,"#pi,K,p");
  leg->AddEntry(PHOSpi,"#pi^{#pm}");
  leg->AddEntry(PHOSk,"K^{#pm}");
  leg->AddEntry(PHOSp,"p,#bar{p}");
  leg->SetFillStyle(0);
  leg->SetFillColor(0);
  leg->SetBorderSize(0);
  leg->SetTextSize(0.06);
 leg->Draw();

  TLine *line = new TLine(0.150,0.0,0.150,1.0);
  if(eta)line->Draw();
  line->SetLineWidth(3.0);
  //line->SetLineColor(TColor::kYellow);
  line->SetLineStyle(2);
  TLatex *tex = new TLatex(0.497269,0.0513196,prodname);
  tex->SetTextSize(0.0537634);
  tex->Draw();
  TLatex *tex3 = new TLatex(1.16186,0.28348,"Closed symbols |#eta|<0.12 (PHOS)");
  tex3->SetTextSize(0.0537634);
  tex3->Draw();
  TLatex *tex4 = new TLatex(1.16186,0.213221,"Open symbols |#eta|<0.70 (EMCal)");
  tex4->SetTextSize(0.0537634);
  tex4->Draw();
  TLatex *tex2 = new TLatex(0.164016,0.860826,"TPC cut-off 150 MeV/c");
  tex2->SetTextSize(0.0537634);
  if(eta) tex2->Draw();


  TLine *line2 = new TLine(0.10,0.0,0.10,1.0);
  line2->SetLineWidth(3.0);
  TLatex *tex5 = new TLatex(0.10817,0.924976,"ITS cut-off 100 MeV/c");
  tex5->SetTextSize(0.0537634);
  line2->SetLineStyle(2);
  tex5->SetTextColor(4);
  line2->SetLineColor(4);
  if(!TPC && eta){
    line2->Draw();
    tex5->Draw();
  }
  if(!TPC){
    c->SaveAs("pics/CorrEfficiency.eps");
    c->SaveAs("pics/CorrEfficiency.png");
  }
  else{
    c->SaveAs("pics/CorrEfficiencyTPCITS.eps");
    c->SaveAs("pics/CorrEfficiencyTPCITS.png");
  }
}
Commit	Line	Data
b610fb37	1	//Christine Nattrass, University of Tennessee at Knoxville
	2	//This macro is for calculating the single track efficiency for the hadronic transverse energy measurement
	3	//Uses the output of AliAnalysisTaskHadEt
	4	//This is not actually what gets used in the correction class AliAnalysisHadEtCorrections - that is done in the macro GetCorrections.C - but this is useful for making plots and playing around with different options
	5
	6	//this function calculates the efficiency propagating errors properly
86d9d0dc	7	TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name)
	8	{
	9	if(!numerator){
	10	cerr<<"Error: numerator does not exist!"<<endl;
	11	return NULL;
	12	}
	13	if(!denominator){
	14	cerr<<"Error: denominator does not exist!"<<endl;
	15	return NULL;
	16	}
	17	TH1D* result = (TH1D*)numerator->Clone(name);
	18	Int_t nbins = numerator->GetNbinsX();
	19	for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
	20	Double_t numeratorVal = numerator->GetBinContent(ibin);
	21	Double_t denominatorVal = denominator->GetBinContent(ibin);
	22	// Check if the errors are right or the thing is scaled
	23	Double_t numeratorValErr = numerator->GetBinError(ibin);
	24	if (!(numeratorValErr==0. \|\| numeratorVal ==0.) ) {
	25	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
	26	numeratorVal /= rescale;
	27	}
	28	Double_t denominatorValErr = denominator->GetBinError(ibin);
	29	if (!(denominatorValErr==0. \|\| denominatorVal==0. )) {
	30	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
	31	denominatorVal /= rescale;
	32	}
	33	Double_t quotient = 0.;
	34	if (denominatorVal!=0.) {
	35	quotient = numeratorVal/denominatorVal;
	36	}
	37	Double_t quotientErr=0;
	38	quotientErr = TMath::Sqrt(
	39	(numeratorVal+1.0)/(denominatorVal+2.0)*
	40	((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
	41	result->SetBinContent(ibin,quotient);
	42	result->SetBinError(ibin,quotientErr);
	43	//cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
	44	}
	45	return result;
	46	}
	47
	48
b610fb37	49	//This is a somewhat messy function that gets the efficiency for different particles
86d9d0dc	50	TH1D GetHisto(float cut = 0.12, char name, int mycase, bool eta, int color, int marker,bool TPC, bool ITS){
	51	//TFile *file = new TFile("Et.ESD.new.sim.merged.root");
	52	TFile *file = new TFile("Et.ESD.new.sim.LHC10d4.pp.merged.root");
	53	TList *list = file->FindObject("out2");
	54	char *myname = "ITS";
	55	if(TPC&&!ITS) myname = "TPC";
	56	if(TPC&&ITS) myname = "TPCITS";
	57	cout<<"Using tracks from "<<myname<<" for efficiency"<<endl;
	58	TH2F *numeratorParent;
	59	switch(mycase){
	60	case 0:
	61	numeratorParent = (TH2F)((TH2F) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoHadron");
	62	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
	63	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
	64	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")));
	65	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")));
	66	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
	67	break;
	68	case 1://pion
	69	numeratorParent = (TH2F)((TH2F) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoPion");
	70	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
	71	break;
	72	case 2://kaon
	73	numeratorParent = (TH2F)((TH2F) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")))->Clone("RecoKaon");
	74	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
	75	break;
	76	case 3://proton
	77	numeratorParent = (TH2F)((TH2F) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")))->Clone("RecoProton");
	78	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
	79	break;
	80	case 4://electron
	81	numeratorParent = (TH2F)((TH2F) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EPlus")))->Clone("RecoElectron");
	82	numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EMinus")));
	83	break;
	84	}
	85	TH2F *denominatorParent;
	86	switch(mycase){
	87	case 0:
	88	denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedChargedHadron"))->Clone("RecoHadron");
	89	// denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoHadron");
	90	// denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
	91	// denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
	92	// denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKPlus"));
	93	// denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedProton"));
	94	// denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
	95	break;
	96	case 1://pion
	97	denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoPion");
	98	denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
	99	break;
	100	case 2://kaon
	101	denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedKPlus"))->Clone("RecoKaon");
	102	denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
	103	break;
	104	case 3://proton
	105	denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedProton"))->Clone("RecoProton");
	106	denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
	107	break;
	108	case 4://electron
	109	denominatorParent = (TH2F)((TH2F) out2->FindObject("EtNSimulatedEPlus"))->Clone("RecoElectron");
	110	denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedEMinus"));
	111	break;
	112	}
	113	numeratorParent->Sumw2();
114	denominatorParent->Sumw2();
115	TH1D *denominator;
116	TH1D *numerator;
117	if(eta){
118	int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
119	int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
120	cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
121	denominator = denominatorParent->ProjectionX(Form("garbage%s",name),lowbin,highbin);
122	numerator = numeratorParent->ProjectionX(name,lowbin,highbin);
123	}
124	else{
125	int lowbin = denominatorParent->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
126	int highbin = denominatorParent->GetXaxis()->GetNbins();
b610fb37	127	cout<<"Here Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
86d9d0dc	128	numerator = numeratorParent->ProjectionY(name,lowbin,highbin);
	129	denominator = denominatorParent->ProjectionY(Form("denominator%s",name),lowbin,highbin);
	130	}
	131	delete numeratorParent;
	132	delete denominatorParent;
	133	//numerator->Divide(denominator);
	134	TH1D result = bayneseffdiv((TH1D) numerator,(TH1D*)denominator,name);
	135	//result->Rebin(2);
	136	//result->Scale(0.5);
	137	result->SetYTitle("Efficiency");
	138	result->GetYaxis()->SetTitleOffset(0.8);
	139	result->GetXaxis()->SetTitleOffset(0.8);
	140	result->GetYaxis()->SetLabelSize(0.05);
	141	result->GetXaxis()->SetLabelSize(0.05);
	142	result->GetYaxis()->SetTitleSize(0.05);
	143	result->GetXaxis()->SetTitleSize(0.05);
	144	result->SetMarkerColor(color);
	145	result->SetLineColor(color);
	146	result->SetMarkerStyle(marker);
	147	//result->Draw("e");
	148	return result;
	149
	150	}
	151
b610fb37	152
	153	//this is a method that makes pretty plots
	154	void CorrEfficiency(char prodname= "LHC10d4", char shortprodname = "LHC10d4 PYTHIA D6T 7 TeV p+p", bool TPC = true,bool ITS = true, bool eta = false){
86d9d0dc	155
	156	gStyle->SetOptTitle(0);
	157	gStyle->SetOptStat(0);
	158	gStyle->SetOptFit(0);
	159	TCanvas *c = new TCanvas("c","c",600,400);
	160	c->SetTopMargin(0.02);
	161	c->SetRightMargin(0.02);
	162	c->SetBorderSize(0);
	163	c->SetFillColor(0);
	164	c->SetFillColor(0);
	165	c->SetBorderMode(0);
	166	c->SetFrameFillColor(0);
	167	c->SetFrameBorderMode(0);
	168	//c->SetLogx();
	169
	170	int colortotal = 1;
	171	int colorpi = 2;
	172	int colork = 3;
	173	int colorp = 4;
	174	int phosmarker = 20;
	175	int emcalmarker = 24;
	176	float ptcut1 = 0.05;
	177	float ptcut2 = 0.1;
b610fb37	178	float phoscut = 0.12;
	179	float emcalcut = 0.7;
	180	if(!eta){
	181	phoscut = 0.1;
	182	emcalcut = 0.15;
	183	}
	184	TH1D *PHOStotal = GetHisto(phoscut,"PHOStotal",0,eta,colortotal,phosmarker,TPC,ITS);
	185	TH1D *PHOSpi = GetHisto(phoscut,"PHOSpi",1,eta,colorpi,phosmarker,TPC,ITS);
	186	TH1D *PHOSp = GetHisto(phoscut,"PHOSp",2,eta,colork,phosmarker,TPC,ITS);
	187	TH1D *PHOSk = GetHisto(phoscut,"PHOSk",3,eta,colorp,phosmarker,TPC,ITS);
	188	if(eta) PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));
	189	//if(ITS&&!TPC){PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));}
	190	//else{PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.0),PHOStotal->GetXaxis()->FindBin(3.0));}
86d9d0dc	191	PHOStotal->SetMinimum(0.0);
	192	PHOStotal->SetMaximum(1.0);
	193	//parameters[centbin][0]*exp(-pow(parameters[centbin][1]/pt,parameters[centbin][2]))
	194	TF1 func = new TF1("func","[0]exp(-pow([1]/x,[2]))");
	195	func->SetParameter(0,.9);
	196	func->SetParameter(1,.05);
	197	func->SetParLimits(1,1e-3,1);
	198	func->SetParameter(2,.1);
	199	//PHOStotal->Fit(func);
	200	PHOStotal->Draw();
	201	PHOSpi->Draw("same");
	202	PHOSp->Draw("same");
	203	PHOSk->Draw("same");
b610fb37	204	TH1D *EMCALtotal = GetHisto(emcalcut,"EMCALtotal",0,eta,colortotal,emcalmarker,TPC,ITS);
	205	TH1D *EMCALpi = GetHisto(emcalcut,"EMCALpi",1,eta,colorpi,emcalmarker,TPC,ITS);
	206	TH1D *EMCALp = GetHisto(emcalcut,"EMCALp",2,eta,colork,emcalmarker,TPC,ITS);
	207	TH1D *EMCALk = GetHisto(emcalcut,"EMCALk",3,eta,colorp,emcalmarker,TPC,ITS);
86d9d0dc	208	EMCALtotal->Draw("same");
	209	EMCALpi->Draw("same");
	210	EMCALp->Draw("same");
	211	EMCALk->Draw("same");
	212
	213
	214	TLegend *leg = new TLegend(0.22651,0.247312,0.370805,0.438172);
	215	leg->AddEntry(PHOStotal,"#pi,K,p");
	216	leg->AddEntry(PHOSpi,"#pi^{#pm}");
	217	leg->AddEntry(PHOSk,"K^{#pm}");
	218	leg->AddEntry(PHOSp,"p,#bar{p}");
	219	leg->SetFillStyle(0);
	220	leg->SetFillColor(0);
	221	leg->SetBorderSize(0);
	222	leg->SetTextSize(0.06);
	223	leg->Draw();
	224
b610fb37	225	TLine *line = new TLine(0.150,0.0,0.150,1.0);
b610fb37	226	if(eta)line->Draw();
86d9d0dc	227	line->SetLineWidth(3.0);
	228	//line->SetLineColor(TColor::kYellow);
	229	line->SetLineStyle(2);
	230	TLatex *tex = new TLatex(0.497269,0.0513196,prodname);
	231	tex->SetTextSize(0.0537634);
	232	tex->Draw();
	233	TLatex *tex3 = new TLatex(1.16186,0.28348,"Closed symbols \|#eta\|<0.12 (PHOS)");
	234	tex3->SetTextSize(0.0537634);
	235	tex3->Draw();
	236	TLatex *tex4 = new TLatex(1.16186,0.213221,"Open symbols \|#eta\|<0.70 (EMCal)");
	237	tex4->SetTextSize(0.0537634);
	238	tex4->Draw();
b610fb37	239	TLatex *tex2 = new TLatex(0.164016,0.860826,"TPC cut-off 150 MeV/c");
86d9d0dc	240	tex2->SetTextSize(0.0537634);
b610fb37	241	if(eta) tex2->Draw();
	242
	243
	244	TLine *line2 = new TLine(0.10,0.0,0.10,1.0);
	245	line2->SetLineWidth(3.0);
	246	TLatex *tex5 = new TLatex(0.10817,0.924976,"ITS cut-off 100 MeV/c");
	247	tex5->SetTextSize(0.0537634);
	248	line2->SetLineStyle(2);
	249	tex5->SetTextColor(4);
	250	line2->SetLineColor(4);
	251	if(!TPC && eta){
	252	line2->Draw();
	253	tex5->Draw();
	254	}
	255	if(!TPC){
86d9d0dc	256	c->SaveAs("pics/CorrEfficiency.eps");
	257	c->SaveAs("pics/CorrEfficiency.png");
	258	}
	259	else{
	260	c->SaveAs("pics/CorrEfficiencyTPCITS.eps");
	261	c->SaveAs("pics/CorrEfficiencyTPCITS.png");
	262	}
	263	}