[u/mrichter/AliRoot.git] / STEER / AliPIDComparison.C

//*************************************************************************
//                Example of using the PID classes
//  which calculates the efficiency and contamination of the combined PID.
//
//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
//*************************************************************************

#if !defined( __CINT__) || defined(__MAKECINT__)
  #include <TMath.h>
  #include <TError.h>
  #include <Riostream.h>
  #include <TF1.h>
  #include <TH1F.h>
  #include <TH2F.h>
  #include <TParticle.h>
  #include <TCanvas.h>
  #include <TBenchmark.h>
  #include <TFile.h>
  #include <TTree.h>
  #include <TROOT.h>

  #include <AliStack.h>
  #include <AliRunLoader.h>
  #include <AliRun.h>
  #include <AliTPCpidESD.h>
  #include <AliESDEvent.h>
#endif

extern AliRun *gAlice;
extern TBenchmark *gBenchmark;
extern TROOT *gROOT;

static Int_t allpisel=0;
static Int_t allkasel=0;
static Int_t allprsel=0;
static Int_t allnsel=0;


Double_t tpcBethe(Double_t *x, Double_t *par) {
  //
  // This is the TPC Bethe-Bloch function normalised to 1 at the minimum
  //
  Double_t p=x[0];   //momentum
  Double_t m=par[0]; //mass
  Double_t bg=p/m;
  return AliTPCpidESD::Bethe(bg);  
}


Int_t AliPIDComparison(const Char_t *dir=".") { 
   gBenchmark->Start("AliPIDComparison");

   ::Info("AliPIDComparison.C","Doing comparison...");

   Double_t pi=0.2,pa=3;

   TH2F *tpcHist=(TH2F*)gROOT->FindObject("tpcHist");
   if (!tpcHist)
     tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,6.);
   tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
   tpcHist->SetMarkerStyle(8); 
   tpcHist->SetMarkerSize(0.3);
 
   const Char_t *hname[]={
    "piG","piR","piF","piGood","piFake",
    "kaG","kaR","kaF","kaGood","kaFake",
    "prG","prR","prF","prGood","prFake"
   };
   Int_t nh=sizeof(hname)/sizeof(const Char_t *);
   TH1F **hprt=new TH1F*[nh]; 

   for (Int_t i=0; i<nh; i++) {
     hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
     if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
   }
   TH1F *piG=hprt[0];
   TH1F *piR=hprt[1];
   TH1F *piF=hprt[2];
   TH1F *piGood=hprt[3]; 
        piGood->SetTitle("Combined PID for pions");
        piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
   TH1F *piFake=hprt[4]; 
        piFake->SetLineColor(2);
   TH1F *kaG=hprt[5];
   TH1F *kaR=hprt[6];
   TH1F *kaF=hprt[7];
   TH1F *kaGood=hprt[8];
        kaGood->SetTitle("Combined PID for kaons"); 
        kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
   TH1F *kaFake=hprt[9]; 
        kaFake->SetLineColor(2);
   TH1F *prG=hprt[10];
   TH1F *prR=hprt[11];
   TH1F *prF=hprt[12];
   TH1F *prGood=hprt[13];
        prGood->SetTitle("Combined PID for protons"); 
        prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
   TH1F *prFake=hprt[14]; 
        prFake->SetLineColor(2);

   delete[] hprt;

   Char_t fname[100];

   if (gAlice) {
      delete AliRunLoader::Instance();
      delete gAlice;
      gAlice=0;
   }
   sprintf(fname,"%s/galice.root",dir);
   AliRunLoader *rl = AliRunLoader::Open(fname);
   if (rl == 0x0) {
      ::Error("AliPIDComparison.C","Can not open session !");
      return 1;
   }
   if (rl->LoadgAlice()) {
      ::Error("AliPIDComparison.C","LoadgAlice returned error !");
      delete rl;
      return 1;
   }
   if (rl->LoadHeader()) {
      ::Error("AliPIDComparison.C","LoadHeader returned error !");
      delete rl;
      return 1;
   }
   rl->LoadKinematics();

   sprintf(fname,"%s/AliESDs.root",dir);
   TFile *ef=TFile::Open(fname);
   if (!ef || !ef->IsOpen()) {
      ::Error("AliPIDComparison.C","Can't AliESDs.root !"); 
      delete rl;
      return 1;
   }
   AliESDEvent* event = new AliESDEvent();
   TTree* tree = (TTree*) ef->Get("esdTree");
   if (!tree) {
      ::Error("AliPIDComparison.C", "no ESD tree found");
      delete rl;
      return 1;
   }
   event->ReadFromTree(tree);


   //****** Tentative particle type "concentrations"
   Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
   AliPID::SetPriors(c);


   //******* The loop over events
   Int_t e=0;
   while (tree->GetEvent(e)) {
      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";

      rl->GetEvent(e);
 
      e++;

      Int_t ntrk=event->GetNumberOfTracks();
      cerr<<"Number of ESD tracks : "<<ntrk<<endl; 

      Int_t pisel=0,kasel=0,prsel=0,nsel=0;

      AliStack *stack = rl->Stack();

      while (ntrk--) {
        AliESDtrack *t=event->GetTrack(ntrk);

	//*** Some track quality cuts ****
        if (t->GetITSclusters(0) < 6 ) continue;
        if (t->GetTPCclusters(0) < 60) continue;
        //if (t->GetTRDclusters(0) < 60) continue;

        if (!t->IsOn(AliESDtrack::kITSpid)) continue;
        if (!t->IsOn(AliESDtrack::kTPCpid)) continue;
        //if (!t->IsOn(AliESDtrack::kTRDpid)) continue;
        if (!t->IsOn(AliESDtrack::kTOFpid)) continue;
        {
           nsel++;

           Double_t p=t->GetInnerParam()->GetP();
	   Double_t dedx=t->GetTPCsignal()/47.;
	   tpcHist->Fill(p,dedx,1);

           Int_t lab=TMath::Abs(t->GetLabel());
           TParticle *part=stack->Particle(lab);
           Int_t code=part->GetPdgCode();

           Double_t r[10]; t->GetESDpid(r);

           AliPID pid(r);

           Double_t w[10];
           w[0]=pid.GetProbability(AliPID::kElectron);
           w[1]=pid.GetProbability(AliPID::kMuon);
           w[2]=pid.GetProbability(AliPID::kPion);
           w[3]=pid.GetProbability(AliPID::kKaon);
           w[4]=pid.GetProbability(AliPID::kProton);

           if (TMath::Abs(code)==2212) prR->Fill(p);
           if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
	      prsel++;
	      prG->Fill(p);
              if (TMath::Abs(code)!=2212) prF->Fill(p);
           }

           if (TMath::Abs(code)==321) kaR->Fill(p);
           if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
	      kasel++;
	      kaG->Fill(p);
              if (TMath::Abs(code)!=321) kaF->Fill(p);
           }

	   if (TMath::Abs(code)==211) piR->Fill(p);
	   if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion
	      pisel++;
	      piG->Fill(p);
	      if (TMath::Abs(code)!=211) piF->Fill(p);
           }
	   /*
           if (TMath::Abs(code)==11) piR->Fill(p);
           if (w[0]>w[4] && w[0]>w[3] && w[0]>w[3] && w[0]>w[1]) {//electron
	      pisel++;
	      piG->Fill(p);
              if (TMath::Abs(code)!=11) piF->Fill(p);
           }
	   */
	}
      }
      cout<<"Number of selected ESD tracks : "<<nsel<<endl;
      cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
      cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
      cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;

      allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;

   } // ***** End of the loop over events

   delete event;
   delete tree;
   ef->Close();

   TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
   if (!c1) {
      c1=new TCanvas("c1","",0,0,600,1000);
      c1->Divide(1,4);
   }

   c1->cd(1);
   tpcHist->Draw();

   TF1 *fpi=(TF1*)gROOT->FindObject("piBethe");
   if (!fpi) {
      fpi=new TF1("piBethe",tpcBethe,pi,pa,1);
      fpi->SetLineColor(2);
      fpi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));
   }
   fpi->Draw("same");

   TF1 *fka=(TF1*)gROOT->FindObject("kaBethe");
   if (!fka) {
      fka=new TF1("kaBethe",tpcBethe,pi,pa,1);
      fka->SetLineColor(3);
      fka->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
   }
   fka->Draw("same");

   TF1 *fpr=(TF1*)gROOT->FindObject("prBethe");
   if (!fpr) {
      fpr=new TF1("prBethe",tpcBethe,pi,pa,1);
      fpr->SetLineColor(4);
      fpr->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
   }
   fpr->Draw("same");


   c1->cd(2);
   //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
   piGood->Add(piG,piF,1,-1);
   piGood->Divide(piGood,piR,1,1,"b");
   piFake->Divide(piF,piG,1,1,"b");
   piGood->SetMinimum(0);
   piGood->Draw("hist");
   piFake->Draw("same");

   c1->cd(3);
   //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
   kaGood->Add(kaG,kaF,1,-1);
   kaGood->Divide(kaGood,kaR,1,1,"b");
   kaFake->Divide(kaF,kaG,1,1,"b");
   kaGood->SetMinimum(0);
   kaGood->Draw("hist");
   kaFake->Draw("same");

   c1->cd(4);
   //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
   prGood->Add(prG,prF,1,-1);
   prGood->Divide(prGood,prR,1,1,"b");
   prFake->Divide(prF,prG,1,1,"b");
   prGood->SetMinimum(0);
   prGood->Draw("hist");
   prFake->Draw("same");

   c1->Update();

   cout<<endl;
   cout<<endl;
   e=(Int_t)piR->GetEntries();
   Int_t o=(Int_t)piG->GetEntries();
   if (e*o) cout<<"Efficiency (contamination) for pions : "<<
      piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl; 
   e=(Int_t)kaR->GetEntries();
   o=(Int_t)kaG->GetEntries();
   if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
      kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl; 
   e=(Int_t)prR->GetEntries();
   o=(Int_t)prG->GetEntries();
   if (e*o) cout<<"Efficiency (contamination) for protons : "<<
      prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
   cout<<endl; 

   cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
   cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
   cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
   cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;

   ef->Close();
   TFile fc("AliPIDComparison.root","RECREATE");
   c1->Write();
   fc.Close();

   gBenchmark->Stop("AliPIDComparison");
   gBenchmark->Show("AliPIDComparison");

   delete rl;
   return 0;
}
Commit	Line	Data
8781a4f1	1	//*************************************************************************
	2	// Example of using the PID classes
	3	// which calculates the efficiency and contamination of the combined PID.
	4	//
	5	// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
	6	//*************************************************************************
	7
	8	#if !defined( __CINT__) \|\| defined(__MAKECINT__)
	9	#include <TMath.h>
	10	#include <TError.h>
	11	#include <Riostream.h>
	12	#include <TF1.h>
	13	#include <TH1F.h>
	14	#include <TH2F.h>
	15	#include <TParticle.h>
	16	#include <TCanvas.h>
	17	#include <TBenchmark.h>
	18	#include <TFile.h>
	19	#include <TTree.h>
	20	#include <TROOT.h>
	21
	22	#include <AliStack.h>
	23	#include <AliRunLoader.h>
	24	#include <AliRun.h>
	25	#include <AliTPCpidESD.h>
	26	#include <AliESDEvent.h>
	27	#endif
	28
	29	extern AliRun *gAlice;
	30	extern TBenchmark *gBenchmark;
	31	extern TROOT *gROOT;
	32
	33	static Int_t allpisel=0;
	34	static Int_t allkasel=0;
	35	static Int_t allprsel=0;
	36	static Int_t allnsel=0;
	37
	38
	39	Double_t tpcBethe(Double_t x, Double_t par) {
	40	//
	41	// This is the TPC Bethe-Bloch function normalised to 1 at the minimum
	42	//
	43	Double_t p=x[0]; //momentum
	44	Double_t m=par[0]; //mass
	45	Double_t bg=p/m;
	46	return AliTPCpidESD::Bethe(bg);
	47	}
	48
	49
	50	Int_t AliPIDComparison(const Char_t *dir=".") {
	51	gBenchmark->Start("AliPIDComparison");
	52
	53	::Info("AliPIDComparison.C","Doing comparison...");
	54
	55	Double_t pi=0.2,pa=3;
	56
	57	TH2F tpcHist=(TH2F)gROOT->FindObject("tpcHist");
	58	if (!tpcHist)
	59	tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,6.);
	60	tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
	61	tpcHist->SetMarkerStyle(8);
	62	tpcHist->SetMarkerSize(0.3);
	63
	64	const Char_t *hname[]={
65	"piG","piR","piF","piGood","piFake",
66	"kaG","kaR","kaF","kaGood","kaFake",
67	"prG","prR","prF","prGood","prFake"
68	};
69	Int_t nh=sizeof(hname)/sizeof(const Char_t *);
70	TH1F *hprt=new TH1F[nh];
71
72	for (Int_t i=0; i<nh; i++) {
73	hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
74	if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
75	}
76	TH1F *piG=hprt[0];
77	TH1F *piR=hprt[1];
78	TH1F *piF=hprt[2];
79	TH1F *piGood=hprt[3];
80	piGood->SetTitle("Combined PID for pions");
81	piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
82	TH1F *piFake=hprt[4];
83	piFake->SetLineColor(2);
84	TH1F *kaG=hprt[5];
85	TH1F *kaR=hprt[6];
86	TH1F *kaF=hprt[7];
87	TH1F *kaGood=hprt[8];
88	kaGood->SetTitle("Combined PID for kaons");
89	kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
90	TH1F *kaFake=hprt[9];
91	kaFake->SetLineColor(2);
92	TH1F *prG=hprt[10];
93	TH1F *prR=hprt[11];
94	TH1F *prF=hprt[12];
95	TH1F *prGood=hprt[13];
96	prGood->SetTitle("Combined PID for protons");
97	prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
98	TH1F *prFake=hprt[14];
99	prFake->SetLineColor(2);
100
101	delete[] hprt;
102
103	Char_t fname[100];
104
105	if (gAlice) {
33c3c91a	106	delete AliRunLoader::Instance();
8781a4f1	107	delete gAlice;
	108	gAlice=0;
	109	}
	110	sprintf(fname,"%s/galice.root",dir);
	111	AliRunLoader *rl = AliRunLoader::Open(fname);
	112	if (rl == 0x0) {
	113	::Error("AliPIDComparison.C","Can not open session !");
	114	return 1;
	115	}
	116	if (rl->LoadgAlice()) {
	117	::Error("AliPIDComparison.C","LoadgAlice returned error !");
	118	delete rl;
	119	return 1;
	120	}
	121	if (rl->LoadHeader()) {
	122	::Error("AliPIDComparison.C","LoadHeader returned error !");
	123	delete rl;
	124	return 1;
	125	}
	126	rl->LoadKinematics();
	127
	128	sprintf(fname,"%s/AliESDs.root",dir);
	129	TFile *ef=TFile::Open(fname);
	130	if (!ef \|\| !ef->IsOpen()) {
	131	::Error("AliPIDComparison.C","Can't AliESDs.root !");
	132	delete rl;
	133	return 1;
	134	}
	135	AliESDEvent* event = new AliESDEvent();
	136	TTree* tree = (TTree*) ef->Get("esdTree");
	137	if (!tree) {
	138	::Error("AliPIDComparison.C", "no ESD tree found");
	139	delete rl;
	140	return 1;
	141	}
	142	event->ReadFromTree(tree);
	143
	144
	145	//****** Tentative particle type "concentrations"
	146	Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
	147	AliPID::SetPriors(c);
	148
	149
	150	//******* The loop over events
	151	Int_t e=0;
	152	while (tree->GetEvent(e)) {
	153	cout<<endl<<endl<<"******* Processing event number: "<<e<<"*****\n";
	154
	155	rl->GetEvent(e);
	156
	157	e++;
	158
	159	Int_t ntrk=event->GetNumberOfTracks();
	160	cerr<<"Number of ESD tracks : "<<ntrk<<endl;
	161
	162	Int_t pisel=0,kasel=0,prsel=0,nsel=0;
	163
	164	AliStack *stack = rl->Stack();
	165
	166	while (ntrk--) {
	167	AliESDtrack *t=event->GetTrack(ntrk);
	168
	169	//* Some track quality cuts **
	170	if (t->GetITSclusters(0) < 6 ) continue;
171	if (t->GetTPCclusters(0) < 60) continue;
172	//if (t->GetTRDclusters(0) < 60) continue;
173
174	if (!t->IsOn(AliESDtrack::kITSpid)) continue;
175	if (!t->IsOn(AliESDtrack::kTPCpid)) continue;
176	//if (!t->IsOn(AliESDtrack::kTRDpid)) continue;
177	if (!t->IsOn(AliESDtrack::kTOFpid)) continue;
178	{
179	nsel++;
180
181	Double_t p=t->GetInnerParam()->GetP();
182	Double_t dedx=t->GetTPCsignal()/47.;
183	tpcHist->Fill(p,dedx,1);
184
185	Int_t lab=TMath::Abs(t->GetLabel());
186	TParticle *part=stack->Particle(lab);
187	Int_t code=part->GetPdgCode();
188
189	Double_t r[10]; t->GetESDpid(r);
190
191	AliPID pid(r);
192
193	Double_t w[10];
194	w[0]=pid.GetProbability(AliPID::kElectron);
195	w[1]=pid.GetProbability(AliPID::kMuon);
196	w[2]=pid.GetProbability(AliPID::kPion);
197	w[3]=pid.GetProbability(AliPID::kKaon);
198	w[4]=pid.GetProbability(AliPID::kProton);
199
200	if (TMath::Abs(code)==2212) prR->Fill(p);
201	if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
202	prsel++;
203	prG->Fill(p);
204	if (TMath::Abs(code)!=2212) prF->Fill(p);
205	}
206
207	if (TMath::Abs(code)==321) kaR->Fill(p);
208	if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
209	kasel++;
210	kaG->Fill(p);
211	if (TMath::Abs(code)!=321) kaF->Fill(p);
212	}
213
214	if (TMath::Abs(code)==211) piR->Fill(p);
215	if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion
216	pisel++;
217	piG->Fill(p);
218	if (TMath::Abs(code)!=211) piF->Fill(p);
219	}
220	/*
221	if (TMath::Abs(code)==11) piR->Fill(p);
222	if (w[0]>w[4] && w[0]>w[3] && w[0]>w[3] && w[0]>w[1]) {//electron
223	pisel++;
224	piG->Fill(p);
225	if (TMath::Abs(code)!=11) piF->Fill(p);
226	}
227	*/
228	}
229	}
230	cout<<"Number of selected ESD tracks : "<<nsel<<endl;
231	cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
232	cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
233	cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
234
235	allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
236
237	} // ***** End of the loop over events
238
239	delete event;
240	delete tree;
241	ef->Close();
242
243	TCanvas c1=(TCanvas)gROOT->FindObject("c1");
244	if (!c1) {
245	c1=new TCanvas("c1","",0,0,600,1000);
246	c1->Divide(1,4);
247	}
248
249	c1->cd(1);
250	tpcHist->Draw();
251
252	TF1 fpi=(TF1)gROOT->FindObject("piBethe");
253	if (!fpi) {
254	fpi=new TF1("piBethe",tpcBethe,pi,pa,1);
255	fpi->SetLineColor(2);
256	fpi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));
257	}
258	fpi->Draw("same");
259
260	TF1 fka=(TF1)gROOT->FindObject("kaBethe");
261	if (!fka) {
262	fka=new TF1("kaBethe",tpcBethe,pi,pa,1);
263	fka->SetLineColor(3);
264	fka->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
265	}
266	fka->Draw("same");
267
268	TF1 fpr=(TF1)gROOT->FindObject("prBethe");
269	if (!fpr) {
270	fpr=new TF1("prBethe",tpcBethe,pi,pa,1);
271	fpr->SetLineColor(4);
272	fpr->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
273	}
274	fpr->Draw("same");
275
276
277	c1->cd(2);
278	//piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
279	piGood->Add(piG,piF,1,-1);
280	piGood->Divide(piGood,piR,1,1,"b");
281	piFake->Divide(piF,piG,1,1,"b");
282	piGood->SetMinimum(0);
283	piGood->Draw("hist");
284	piFake->Draw("same");
285
286	c1->cd(3);
287	//kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
288	kaGood->Add(kaG,kaF,1,-1);
289	kaGood->Divide(kaGood,kaR,1,1,"b");
290	kaFake->Divide(kaF,kaG,1,1,"b");
291	kaGood->SetMinimum(0);
292	kaGood->Draw("hist");
293	kaFake->Draw("same");
294
295	c1->cd(4);
296	//prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
297	prGood->Add(prG,prF,1,-1);
298	prGood->Divide(prGood,prR,1,1,"b");
299	prFake->Divide(prF,prG,1,1,"b");
300	prGood->SetMinimum(0);
301	prGood->Draw("hist");
302	prFake->Draw("same");
303
304	c1->Update();
305
306	cout<<endl;
307	cout<<endl;
308	e=(Int_t)piR->GetEntries();
309	Int_t o=(Int_t)piG->GetEntries();
310	if (e*o) cout<<"Efficiency (contamination) for pions : "<<
311	piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
312	e=(Int_t)kaR->GetEntries();
313	o=(Int_t)kaG->GetEntries();
314	if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
315	kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
316	e=(Int_t)prR->GetEntries();
317	o=(Int_t)prG->GetEntries();
318	if (e*o) cout<<"Efficiency (contamination) for protons : "<<
319	prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
320	cout<<endl;
321
322	cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
323	cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
324	cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
325	cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
326
327	ef->Close();
328	TFile fc("AliPIDComparison.root","RECREATE");
329	c1->Write();
330	fc.Close();
331
332	gBenchmark->Stop("AliPIDComparison");
333	gBenchmark->Show("AliPIDComparison");
334
335	delete rl;
336	return 0;
337	}