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.,3.,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.,3.);
  35   TH1F *piR=new TH1F("piR","",20,0.,3.);
  36   TH1F *piF=new TH1F("piF","",20,0.,3.);
  37   TH1F *piGood=new TH1F("piGood","Combined PID for pions",20,0.,3.);
  38   piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
  39   TH1F *piFake=new TH1F("piFake","",20,0.,3.); piFake->SetLineColor(2);
  40
  41   TH1F *kaG=new TH1F("kaG","",20,0.,3.);
  42   TH1F *kaR=new TH1F("kaR","",20,0.,3.);
  43   TH1F *kaF=new TH1F("kaF","",20,0.,3.);
  44   TH1F *kaGood=new TH1F("kaGood","Combined PID for kaons",20,0.,3.);
  45   kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
  46   TH1F *kaFake=new TH1F("kaFake","",20,0.,3.); kaFake->SetLineColor(2);
  47
  48   TH1F *prG=new TH1F("prG","",20,0.,3.);
  49   TH1F *prR=new TH1F("prR","",20,0.,3.);
  50   TH1F *prF=new TH1F("prF","",20,0.,3.);
  51   TH1F *prGood=new TH1F("prGood","Combined PID for protons",20,0.,3.);
  52   prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
  53   TH1F *prFake=new TH1F("prFake","",20,0.,3.); 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 }