New geometry (support bars and frames) to be draw
[u/mrichter/AliRoot.git] / STEER / AliESDanalysis.C
CommitLineData
ae982df3 1//********************************************************************
8c6a71ab 2// Example (very naive for the moment) of the data analysis
ae982df3 3// using the ESD classes
8c6a71ab 4// It demonstrates the idea of the "combined PID".
ae982df3 5// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
6//********************************************************************
7
c630aafd 8#if !defined( __CINT__) || defined(__MAKECINT__)
ae982df3 9 #include <Riostream.h>
10 #include "TKey.h"
11 #include "TFile.h"
c630aafd 12 #include "TH1F.h"
ae982df3 13 #include "TH2F.h"
14 #include "TCanvas.h"
15 #include "TStopwatch.h"
c630aafd 16 #include "TParticle.h"
17
18 #include "AliRun.h"
489609b5 19 #include "AliStack.h"
c630aafd 20 #include "AliRunLoader.h"
21 #include "AliLoader.h"
ae982df3 22
23 #include "AliESD.h"
24#endif
25
c630aafd 26extern AliRun *gAlice;
27
ae982df3 28Int_t AliESDanalysis(Int_t nev=1) {
8c6a71ab 29 TH2F *tpcHist=
bb2ceb1f 30 new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
c630aafd 31 tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
8c6a71ab 32 tpcHist->SetMarkerStyle(8);
33 tpcHist->SetMarkerSize(0.3);
34
bb2ceb1f 35 TH1F *piG=new TH1F("piG","",20,0.,4.);
36 TH1F *piR=new TH1F("piR","",20,0.,4.);
37 TH1F *piF=new TH1F("piF","",20,0.,4.);
38 TH1F *piGood=new TH1F("piGood","Combined PID for pions",20,0.,4.);
c630aafd 39 piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
bb2ceb1f 40 TH1F *piFake=new TH1F("piFake","",20,0.,4.); piFake->SetLineColor(2);
c630aafd 41
bb2ceb1f 42 TH1F *kaG=new TH1F("kaG","",20,0.,4.);
43 TH1F *kaR=new TH1F("kaR","",20,0.,4.);
44 TH1F *kaF=new TH1F("kaF","",20,0.,4.);
45 TH1F *kaGood=new TH1F("kaGood","Combined PID for kaons",20,0.,4.);
c630aafd 46 kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
bb2ceb1f 47 TH1F *kaFake=new TH1F("kaFake","",20,0.,4.); kaFake->SetLineColor(2);
c630aafd 48
bb2ceb1f 49 TH1F *prG=new TH1F("prG","",20,0.,4.);
50 TH1F *prR=new TH1F("prR","",20,0.,4.);
51 TH1F *prF=new TH1F("prF","",20,0.,4.);
52 TH1F *prGood=new TH1F("prGood","Combined PID for protons",20,0.,4.);
c630aafd 53 prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
bb2ceb1f 54 TH1F *prFake=new TH1F("prFake","",20,0.,4.); prFake->SetLineColor(2);
c630aafd 55
56 if (gAlice) {
57 delete gAlice->GetRunLoader();
58 delete gAlice;
59 gAlice=0;
60 }
61 AliRunLoader *rl = AliRunLoader::Open("galice.root");
62 if (rl == 0x0) {
63 cerr<<"Can not open session"<<endl;
64 return 1;
65 }
66 if (rl->LoadgAlice()) {
67 cerr<<"LoadgAlice returned error"<<endl;
68 delete rl;
69 return 1;
70 }
71 if (rl->LoadHeader()) {
72 cerr<<"LoadHeader returned error"<<endl;
73 delete rl;
74 return 1;
75 }
76 rl->LoadKinematics();
489609b5 77 AliStack* stack = rl->Stack();
c630aafd 78
ae982df3 79 TFile *ef=TFile::Open("AliESDs.root");
80 if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
81
82 TStopwatch timer;
83 Int_t rc=0,n=0;
84 TKey *key=0;
85 TIter next(ef->GetListOfKeys());
86
8c6a71ab 87 //****** Tentative particle type "concentrations"
c630aafd 88 Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
8c6a71ab 89
ae982df3 90 //******* The loop over events
91 while ((key=(TKey*)next())!=0) {
c630aafd 92 rl->GetEvent(n);
93
ae982df3 94 cerr<<"Processing event number : "<<n++<<endl;
c630aafd 95
ae982df3 96 AliESD *event=(AliESD*)key->ReadObj();
97
98 Int_t ntrk=event->GetNumberOfTracks();
99 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
100 //****** The loop over tracks
c630aafd 101
102 Int_t pisel=0,kasel=0,prsel=0,nsel=0;
103
ae982df3 104 while (ntrk--) {
105 AliESDtrack *t=event->GetTrack(ntrk);
8c6a71ab 106
ae982df3 107 Double_t p=t->GetP();
8c6a71ab 108
ae982df3 109 if (t->GetStatus()&AliESDtrack::kTPCin) {
110 Double_t dedx=t->GetTPCsignal();
111 tpcHist->Fill(p,dedx,1);
112 }
8c6a71ab 113
c630aafd 114 UInt_t status=AliESDtrack::kESDpid;
115 status|=AliESDtrack::kTPCpid;
116 status|=AliESDtrack::kTOFpid;
117
118 if ((t->GetStatus()&status) == status) {
119 nsel++;
120
121 Int_t lab=TMath::Abs(t->GetLabel());
489609b5 122 TParticle *part=stack->Particle(lab);
c630aafd 123 Int_t code=part->GetPdgCode();
124
125 Double_t r[10]; t->GetESDpid(r);
126
127 Double_t rcc=0.;
8c6a71ab 128 Int_t i;
c630aafd 129 for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
130 if (rcc==0.) continue;
8c6a71ab 131
132 //Here we apply Bayes' formula
133 Double_t w[10];
c630aafd 134 for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
135
136 if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
137 prsel++;
138 prG->Fill(p);
139 if (TMath::Abs(code)==2212) prR->Fill(p);
140 else prF->Fill(p);
141 }
142
143 if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
144 kasel++;
145 kaG->Fill(p);
146 if (TMath::Abs(code)==321) kaR->Fill(p);
147 else kaF->Fill(p);
148 }
149
150 if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
151 pisel++;
152 piG->Fill(p);
153 if (TMath::Abs(code)==211) piR->Fill(p);
154 else piF->Fill(p);
155 }
8c6a71ab 156
ae982df3 157 }
8c6a71ab 158
ae982df3 159 }
160 delete event;
c630aafd 161 cerr<<"Number of selected ESD tracks : "<<nsel<<endl;
162 cerr<<"Number of selected pion ESD tracks : "<<pisel<<endl;
163 cerr<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
164 cerr<<"Number of selected proton ESD tracks : "<<prsel<<endl;
ae982df3 165 }
c630aafd 166
ae982df3 167 timer.Stop(); timer.Print();
168
169 TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
c630aafd 170 c1->Divide(1,4);
ae982df3 171
172 c1->cd(1);
173 tpcHist->Draw();
c630aafd 174
ae982df3 175 c1->cd(2);
c630aafd 176 piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
177 piGood->Divide(piR,piG,1,1,"b");
178 piFake->Divide(piF,piG,1,1,"b");
179 piGood->Draw("hist");
180 piFake->Draw("same");
181
182 c1->cd(3);
183 kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
184 kaGood->Divide(kaR,kaG,1,1,"b");
185 kaFake->Divide(kaF,kaG,1,1,"b");
186 kaGood->Draw("hist");
187 kaFake->Draw("same");
188
189 c1->cd(4);
190 prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
191 prGood->Divide(prR,prG,1,1,"b");
192 prFake->Divide(prF,prG,1,1,"b");
193 prGood->Draw("hist");
194 prFake->Draw("same");
ae982df3 195
196 ef->Close();
197
c630aafd 198 delete rl;
199
ae982df3 200 return rc;
201}