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

//********************************************************************
//     Example (very naive for the moment) of using the ESD classes
//        It demonstrates the idea 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 <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 <AliESD.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;

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

   ::Info("AliESDComparison.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.,300.);
   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 gAlice->GetRunLoader();
      delete gAlice;
      gAlice=0;
   }
   sprintf(fname,"%s/galice.root",dir);
   AliRunLoader *rl = AliRunLoader::Open(fname);
   if (rl == 0x0) {
      ::Error("AliESDComparison.C","Can not open session !");
      return 1;
   }
   if (rl->LoadgAlice()) {
      ::Error("AliESDComparison.C","LoadgAlice returned error !");
      delete rl;
      return 1;
   }
   if (rl->LoadHeader()) {
      ::Error("AliESDComparison.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("AliESDComparison.C","Can't AliESDs.root !"); 
      delete rl;
      return 1;
   }
   AliESD* event = new AliESD;
   TTree* tree = (TTree*) ef->Get("esdTree");
   if (!tree) {
      ::Error("AliESDComparison.C", "no ESD tree found");
      delete rl;
      return 1;
   }
   tree->SetBranchAddress("ESD", &event);

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

   //******* 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);

         UInt_t status=AliESDtrack::kESDpid;
         status|=AliESDtrack::kITSpid; 
         status|=AliESDtrack::kTPCpid; 
         status|=AliESDtrack::kTOFpid; 

        if ((t->GetStatus()&status) == status) {
           nsel++;

           Double_t p=t->GetP();
	   Double_t dedx=t->GetTPCsignal();
           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);

           Double_t rcc=0.;
           Int_t i;
           for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
           if (rcc==0.) continue;

	   //Here we apply Bayes' formula
           Double_t w[10];
           for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;

           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[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
	      pisel++;
	      piG->Fill(p);
	      if (TMath::Abs(code)!=211) 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;
   ef->Close();

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

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

   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->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->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->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("AliESDComparison.root","RECREATE");
   c1->Write();
   fc.Close();

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

   delete rl;
   return 0;
}
Commit	Line	Data
0717295b	1	//********************************************************************
0717295b	2	// Example (very naive for the moment) of using the ESD classes
42345b19	3	// It demonstrates the idea of the "combined PID".
42345b19	4	// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
0717295b	5	//********************************************************************
	6
	7	#if !defined( __CINT__) \|\| defined(__MAKECINT__)
4fa1e629	8	#include <TMath.h>
4fa1e629	9	#include <TError.h>
0717295b	10	#include <Riostream.h>
4fa1e629	11	#include <TH1F.h>
	12	#include <TH2F.h>
	13	#include <TParticle.h>
	14	#include <TCanvas.h>
	15	#include <TBenchmark.h>
8b462fd8	16	#include <TFile.h>
8b462fd8	17	#include <TTree.h>
4fa1e629	18	#include <TROOT.h>
	19
	20	#include <AliStack.h>
	21	#include <AliRunLoader.h>
	22	#include <AliRun.h>
	23	#include <AliESD.h>
0717295b	24	#endif
	25
	26	extern AliRun *gAlice;
4fa1e629	27	extern TBenchmark *gBenchmark;
42345b19	28	extern TROOT *gROOT;
0717295b	29
4fa1e629	30	static Int_t allpisel=0;
	31	static Int_t allkasel=0;
	32	static Int_t allprsel=0;
	33	static Int_t allnsel=0;
	34
3f7a17bb	35	Int_t AliESDComparison(const Char_t *dir=".") {
4fa1e629	36	gBenchmark->Start("AliESDComparison");
	37
	38	::Info("AliESDComparison.C","Doing comparison...");
	39
	40	Double_t pi=0.2,pa=3;
f2eae663	41
42345b19	42	TH2F tpcHist=(TH2F)gROOT->FindObject("tpcHist");
42345b19	43	if (!tpcHist)
f2eae663	44	tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,300.);
0717295b	45	tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
42345b19	46	tpcHist->SetMarkerStyle(8);
42345b19	47	tpcHist->SetMarkerSize(0.3);
0717295b	48
42345b19	49	const Char_t *hname[]={
	50	"piG","piR","piF","piGood","piFake",
	51	"kaG","kaR","kaF","kaGood","kaFake",
	52	"prG","prR","prF","prGood","prFake"
	53	};
	54	Int_t nh=sizeof(hname)/sizeof(const Char_t *);
	55	TH1F *hprt=new TH1F[nh];
	56
	57	for (Int_t i=0; i<nh; i++) {
	58	hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
f2eae663	59	if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
42345b19	60	}
	61	TH1F *piG=hprt[0];
	62	TH1F *piR=hprt[1];
	63	TH1F *piF=hprt[2];
	64	TH1F *piGood=hprt[3];
	65	piGood->SetTitle("Combined PID for pions");
	66	piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
	67	TH1F *piFake=hprt[4];
	68	piFake->SetLineColor(2);
	69	TH1F *kaG=hprt[5];
	70	TH1F *kaR=hprt[6];
	71	TH1F *kaF=hprt[7];
	72	TH1F *kaGood=hprt[8];
	73	kaGood->SetTitle("Combined PID for kaons");
	74	kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
	75	TH1F *kaFake=hprt[9];
	76	kaFake->SetLineColor(2);
	77	TH1F *prG=hprt[10];
	78	TH1F *prR=hprt[11];
	79	TH1F *prF=hprt[12];
	80	TH1F *prGood=hprt[13];
	81	prGood->SetTitle("Combined PID for protons");
	82	prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
	83	TH1F *prFake=hprt[14];
	84	prFake->SetLineColor(2);
	85
	86	delete[] hprt;
	87
	88	Char_t fname[100];
0717295b	89
	90	if (gAlice) {
	91	delete gAlice->GetRunLoader();
	92	delete gAlice;
	93	gAlice=0;
	94	}
42345b19	95	sprintf(fname,"%s/galice.root",dir);
42345b19	96	AliRunLoader *rl = AliRunLoader::Open(fname);
0717295b	97	if (rl == 0x0) {
4fa1e629	98	::Error("AliESDComparison.C","Can not open session !");
0717295b	99	return 1;
	100	}
	101	if (rl->LoadgAlice()) {
4fa1e629	102	::Error("AliESDComparison.C","LoadgAlice returned error !");
0717295b	103	delete rl;
	104	return 1;
	105	}
	106	if (rl->LoadHeader()) {
4fa1e629	107	::Error("AliESDComparison.C","LoadHeader returned error !");
0717295b	108	delete rl;
	109	return 1;
	110	}
	111	rl->LoadKinematics();
0717295b	112
42345b19	113	sprintf(fname,"%s/AliESDs.root",dir);
42345b19	114	TFile *ef=TFile::Open(fname);
4fa1e629	115	if (!ef \|\| !ef->IsOpen()) {
	116	::Error("AliESDComparison.C","Can't AliESDs.root !");
	117	delete rl;
	118	return 1;
	119	}
8b462fd8	120	AliESD* event = new AliESD;
	121	TTree* tree = (TTree*) ef->Get("esdTree");
	122	if (!tree) {
	123	::Error("AliESDComparison.C", "no ESD tree found");
	124	delete rl;
	125	return 1;
	126	}
	127	tree->SetBranchAddress("ESD", &event);
0717295b	128
	129	//****** Tentative particle type "concentrations"
	130	Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
	131
	132	//******* The loop over events
4fa1e629	133	Int_t e=0;
8b462fd8	134	while (tree->GetEvent(e)) {
4fa1e629	135	cout<<endl<<endl<<"******* Processing event number: "<<e<<"*****\n";
0717295b	136
8b462fd8	137	rl->GetEvent(e);
4fa1e629	138
4fa1e629	139	e++;
0717295b	140
4fa1e629	141	Int_t ntrk=event->GetNumberOfTracks();
4fa1e629	142	cerr<<"Number of ESD tracks : "<<ntrk<<endl;
0717295b	143
4fa1e629	144	Int_t pisel=0,kasel=0,prsel=0,nsel=0;
0717295b	145
4fa1e629	146	AliStack *stack = rl->Stack();
0717295b	147
4fa1e629	148	while (ntrk--) {
4fa1e629	149	AliESDtrack *t=event->GetTrack(ntrk);
0717295b	150
4fa1e629	151	UInt_t status=AliESDtrack::kESDpid;
	152	status\|=AliESDtrack::kITSpid;
	153	status\|=AliESDtrack::kTPCpid;
	154	status\|=AliESDtrack::kTOFpid;
0717295b	155
4fa1e629	156	if ((t->GetStatus()&status) == status) {
4fa1e629	157	nsel++;
f2eae663	158
4fa1e629	159	Double_t p=t->GetP();
	160	Double_t dedx=t->GetTPCsignal();
	161	tpcHist->Fill(p,dedx,1);
0717295b	162
4fa1e629	163	Int_t lab=TMath::Abs(t->GetLabel());
	164	TParticle *part=stack->Particle(lab);
	165	Int_t code=part->GetPdgCode();
0717295b	166
4fa1e629	167	Double_t r[10]; t->GetESDpid(r);
0717295b	168
4fa1e629	169	Double_t rcc=0.;
	170	Int_t i;
	171	for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
	172	if (rcc==0.) continue;
0717295b	173
4fa1e629	174	//Here we apply Bayes' formula
	175	Double_t w[10];
	176	for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
0717295b	177
4fa1e629	178	if (TMath::Abs(code)==2212) prR->Fill(p);
	179	if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
	180	prsel++;
	181	prG->Fill(p);
	182	if (TMath::Abs(code)!=2212) prF->Fill(p);
	183	}
0717295b	184
4fa1e629	185	if (TMath::Abs(code)==321) kaR->Fill(p);
	186	if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
	187	kasel++;
	188	kaG->Fill(p);
	189	if (TMath::Abs(code)!=321) kaF->Fill(p);
	190	}
0717295b	191
4fa1e629	192	if (TMath::Abs(code)==211) piR->Fill(p);
	193	if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
	194	pisel++;
	195	piG->Fill(p);
	196	if (TMath::Abs(code)!=211) piF->Fill(p);
	197	}
	198	}
	199	}
4fa1e629	200	cout<<"Number of selected ESD tracks : "<<nsel<<endl;
	201	cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
	202	cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
	203	cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
0717295b	204
4fa1e629	205	allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
f2eae663	206
4fa1e629	207	} // ***** End of the loop over events
0717295b	208
8b462fd8	209	delete event;
4fa1e629	210	ef->Close();
0717295b	211
42345b19	212	TCanvas c1=(TCanvas)gROOT->FindObject("c1");
4fa1e629	213	if (!c1) {
	214	c1=new TCanvas("c1","",0,0,600,1200);
	215	c1->Divide(1,4);
	216	}
0717295b	217
	218	c1->cd(1);
	219	tpcHist->Draw();
	220
	221	c1->cd(2);
42345b19	222	//piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
f2eae663	223	piGood->Add(piG,piF,1,-1);
f2eae663	224	piGood->Divide(piGood,piR,1,1,"b");
0717295b	225	piFake->Divide(piF,piG,1,1,"b");
	226	piGood->Draw("hist");
	227	piFake->Draw("same");
	228
	229	c1->cd(3);
42345b19	230	//kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
f2eae663	231	kaGood->Add(kaG,kaF,1,-1);
f2eae663	232	kaGood->Divide(kaGood,kaR,1,1,"b");
0717295b	233	kaFake->Divide(kaF,kaG,1,1,"b");
	234	kaGood->Draw("hist");
	235	kaFake->Draw("same");
	236
	237	c1->cd(4);
42345b19	238	//prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
f2eae663	239	prGood->Add(prG,prF,1,-1);
f2eae663	240	prGood->Divide(prGood,prR,1,1,"b");
0717295b	241	prFake->Divide(prF,prG,1,1,"b");
	242	prGood->Draw("hist");
	243	prFake->Draw("same");
	244
4fa1e629	245	c1->Update();
	246
	247	cout<<endl;
	248	cout<<endl;
	249	e=(Int_t)piR->GetEntries();
	250	Int_t o=(Int_t)piG->GetEntries();
	251	if (e*o) cout<<"Efficiency (contamination) for pions : "<<
	252	piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
	253	e=(Int_t)kaR->GetEntries();
	254	o=(Int_t)kaG->GetEntries();
	255	if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
	256	kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
	257	e=(Int_t)prR->GetEntries();
	258	o=(Int_t)prG->GetEntries();
	259	if (e*o) cout<<"Efficiency (contamination) for protons : "<<
	260	prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
	261	cout<<endl;
	262
	263	cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
	264	cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
	265	cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
	266	cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
	267
0717295b	268	ef->Close();
3f7a17bb	269	TFile fc("AliESDComparison.root","RECREATE");
	270	c1->Write();
	271	fc.Close();
0717295b	272
4fa1e629	273	gBenchmark->Stop("AliESDComparison");
4fa1e629	274	gBenchmark->Show("AliESDComparison");
0717295b	275
4fa1e629	276	delete rl;
4fa1e629	277	return 0;
0717295b	278	}