Access to the stack using AliRunLoader (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDanalysis.C
1 //********************************************************************
2 //     Example (very naive for the moment) of the data analysis 
3 //                    using the ESD classes
4 //       It demonstrates the idea of the "combined PID".
5 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
6 //********************************************************************
7
8 #if !defined( __CINT__) || defined(__MAKECINT__)
9   #include <Riostream.h>
10   #include "TKey.h"
11   #include "TFile.h"
12   #include "TH1F.h"
13   #include "TH2F.h"
14   #include "TCanvas.h"
15   #include "TStopwatch.h"
16   #include "TParticle.h"
17
18   #include "AliRun.h"
19   #include "AliStack.h"
20   #include "AliRunLoader.h"
21   #include "AliLoader.h"
22
23   #include "AliESD.h"
24 #endif
25
26 extern AliRun *gAlice;
27
28 Int_t AliESDanalysis(Int_t nev=1) { 
29   TH2F *tpcHist=
30      new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
31    tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
32      tpcHist->SetMarkerStyle(8); 
33      tpcHist->SetMarkerSize(0.3);
34  
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.); 
39   piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
40   TH1F *piFake=new TH1F("piFake","",20,0.,4.); piFake->SetLineColor(2);
41
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.); 
46   kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
47   TH1F *kaFake=new TH1F("kaFake","",20,0.,4.); kaFake->SetLineColor(2);
48
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.); 
53   prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
54   TH1F *prFake=new TH1F("prFake","",20,0.,4.); prFake->SetLineColor(2);
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();
77    AliStack* stack = rl->Stack();
78
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
87    //****** Tentative particle type "concentrations"
88    Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
89
90    //******* The loop over events
91    while ((key=(TKey*)next())!=0) {
92      rl->GetEvent(n);
93
94      cerr<<"Processing event number : "<<n++<<endl;
95
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
101
102      Int_t pisel=0,kasel=0,prsel=0,nsel=0;
103
104      while (ntrk--) {
105        AliESDtrack *t=event->GetTrack(ntrk);
106
107        Double_t p=t->GetP();
108
109        if (t->GetStatus()&AliESDtrack::kTPCin) {
110          Double_t dedx=t->GetTPCsignal();
111          tpcHist->Fill(p,dedx,1);
112        }
113
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());
122          TParticle *part=stack->Particle(lab);
123          Int_t code=part->GetPdgCode();
124
125          Double_t r[10]; t->GetESDpid(r);
126
127          Double_t rcc=0.;
128          Int_t i;
129          for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
130          if (rcc==0.) continue;
131
132          //Here we apply Bayes' formula
133          Double_t w[10];
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          }
156
157        }
158
159      } 
160      delete event;
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;
165    }
166
167    timer.Stop(); timer.Print();
168
169    TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
170    c1->Divide(1,4);
171
172    c1->cd(1);
173    tpcHist->Draw();
174
175    c1->cd(2);
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");
195
196    ef->Close();
197
198    delete rl;
199
200    return rc;
201 }