TRD included in the ESD chain (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 "AliRunLoader.h"
20   #include "AliLoader.h"
21
22   #include "AliESD.h"
23 #endif
24
25 extern AliRun *gAlice;
26
27 Int_t AliESDanalysis(Int_t nev=1) { 
28   TH2F *tpcHist=
29      new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
30    tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
31      tpcHist->SetMarkerStyle(8); 
32      tpcHist->SetMarkerSize(0.3);
33  
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.); 
38   piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
39   TH1F *piFake=new TH1F("piFake","",20,0.,4.); piFake->SetLineColor(2);
40
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.); 
45   kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
46   TH1F *kaFake=new TH1F("kaFake","",20,0.,4.); kaFake->SetLineColor(2);
47
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.); 
52   prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
53   TH1F *prFake=new TH1F("prFake","",20,0.,4.); prFake->SetLineColor(2);
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
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
85    //****** Tentative particle type "concentrations"
86    Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
87
88    //******* The loop over events
89    while ((key=(TKey*)next())!=0) {
90      rl->GetEvent(n);
91
92      cerr<<"Processing event number : "<<n++<<endl;
93
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
99
100      Int_t pisel=0,kasel=0,prsel=0,nsel=0;
101
102      while (ntrk--) {
103        AliESDtrack *t=event->GetTrack(ntrk);
104
105        Double_t p=t->GetP();
106
107        if (t->GetStatus()&AliESDtrack::kTPCin) {
108          Double_t dedx=t->GetTPCsignal();
109          tpcHist->Fill(p,dedx,1);
110        }
111
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.;
126          Int_t i;
127          for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
128          if (rcc==0.) continue;
129
130          //Here we apply Bayes' formula
131          Double_t w[10];
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          }
154
155        }
156
157      } 
158      delete event;
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;
163    }
164
165    timer.Stop(); timer.Print();
166
167    TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
168    c1->Divide(1,4);
169
170    c1->cd(1);
171    tpcHist->Draw();
172
173    c1->cd(2);
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");
193
194    ef->Close();
195
196    delete rl;
197
198    return rc;
199 }