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