First version of combined PID (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 #ifndef __CINT__
9   #include <Riostream.h>
10   #include "TKey.h"
11   #include "TFile.h"
12   #include "TH2F.h"
13   #include "TCanvas.h"
14   #include "TStopwatch.h"
15
16   #include "AliESD.h"
17 #endif
18
19 Int_t AliESDanalysis(Int_t nev=1) { 
20   TH2F *tpcHist=
21      new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
22   tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
23      tpcHist->SetMarkerStyle(8); 
24      tpcHist->SetMarkerSize(0.3);
25  
26   TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
27      elHist->SetMarkerStyle(8); 
28      elHist->SetMarkerSize(0.3);
29
30   TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
31      piHist->SetMarkerColor(2); 
32      piHist->SetMarkerStyle(8); 
33      piHist->SetMarkerSize(0.3);
34
35   TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
36      kaHist->SetMarkerColor(4); 
37      kaHist->SetMarkerStyle(8); 
38      kaHist->SetMarkerSize(0.3);
39
40   TH2F *prHist=new 
41    TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
42    prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
43      prHist->SetMarkerColor(6); 
44      prHist->SetMarkerStyle(8); 
45      prHist->SetMarkerSize(0.3);
46    
47    TFile *ef=TFile::Open("AliESDs.root");
48    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
49
50    TStopwatch timer;
51    Int_t rc=0,n=0;
52    TKey *key=0;
53    TIter next(ef->GetListOfKeys());
54
55    //****** Tentative particle type "concentrations"
56    Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
57
58    //******* The loop over events
59    while ((key=(TKey*)next())!=0) {
60      cerr<<"Processing event number : "<<n++<<endl;
61      AliESD *event=(AliESD*)key->ReadObj();
62
63      Int_t ntrk=event->GetNumberOfTracks();
64      cerr<<"Number of ESD tracks : "<<ntrk<<endl; 
65      //****** The loop over tracks
66      while (ntrk--) {
67        AliESDtrack *t=event->GetTrack(ntrk);
68
69        Double_t p=t->GetP();
70
71        if (t->GetStatus()&AliESDtrack::kTPCin) {
72          Double_t dedx=t->GetTPCsignal();
73          tpcHist->Fill(p,dedx,1);
74        }
75
76        if (t->GetStatus()&AliESDtrack::kESDpid) {
77          Double_t dedx=t->GetTPCsignal();
78          Double_t r[10]; t->GetESDpid(r);
79          Double_t rc=0.;
80          Int_t i;
81          for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
82          if (rc==0.) continue;
83
84          //Here we apply Bayes' formula
85          Double_t w[10];
86          for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
87
88          if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
89          if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
90          if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
91          if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
92        }
93
94      } 
95      delete event;
96    }
97    timer.Stop(); timer.Print();
98
99    TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
100    c1->Divide(1,2);
101
102    c1->cd(1);
103    tpcHist->Draw();
104    c1->cd(2);
105    prHist->Draw();
106    kaHist->Draw("same");
107    piHist->Draw("same");
108    elHist->Draw("same");
109
110    ef->Close();
111
112    return rc;
113 }