]>
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" | |
489609b5 | 19 | #include "AliStack.h" |
c630aafd | 20 | #include "AliRunLoader.h" |
21 | #include "AliLoader.h" | |
ae982df3 | 22 | |
23 | #include "AliESD.h" | |
24 | #endif | |
25 | ||
c630aafd | 26 | extern AliRun *gAlice; |
27 | ||
ae982df3 | 28 | Int_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 | } |