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 }