STEER/AliESDComparison.C

   1 //********************************************************************
   2 //     Example (very naive for the moment) of using the ESD classes
   3 //        It demonstrates the idea of the "combined PID".
   4 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
   5 //********************************************************************
   6
   7 #if !defined( __CINT__) || defined(__MAKECINT__)
   8   #include <TMath.h>
   9   #include <TError.h>
  10   #include <Riostream.h>
  11   #include <TH1F.h>
  12   #include <TH2F.h>
  13   #include <TParticle.h>
  14   #include <TCanvas.h>
  15   #include <TBenchmark.h>
  16   #include <TFile.h>
  17   #include <TTree.h>
  18   #include <TROOT.h>
  19
  20   #include <AliStack.h>
  21   #include <AliRunLoader.h>
  22   #include <AliRun.h>
  23   #include <AliESD.h>
  24 #endif
  25
  26 extern AliRun *gAlice;
  27 extern TBenchmark *gBenchmark;
  28 extern TROOT *gROOT;
  29
  30 static Int_t allpisel=0;
  31 static Int_t allkasel=0;
  32 static Int_t allprsel=0;
  33 static Int_t allnsel=0;
  34
  35 Int_t AliESDComparison(const Char_t *dir=".") {
  36    gBenchmark->Start("AliESDComparison");
  37
  38    ::Info("AliESDComparison.C","Doing comparison...");
  39
  40    Double_t pi=0.2,pa=3;
  41
  42    TH2F *tpcHist=(TH2F*)gROOT->FindObject("tpcHist");
  43    if (!tpcHist)
  44      tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,300.);
  45    tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
  46    tpcHist->SetMarkerStyle(8);
  47    tpcHist->SetMarkerSize(0.3);
  48
  49    const Char_t *hname[]={
  50     "piG","piR","piF","piGood","piFake",
  51     "kaG","kaR","kaF","kaGood","kaFake",
  52     "prG","prR","prF","prGood","prFake"
  53    };
  54    Int_t nh=sizeof(hname)/sizeof(const Char_t *);
  55    TH1F **hprt=new TH1F*[nh];
  56
  57    for (Int_t i=0; i<nh; i++) {
  58      hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
  59      if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
  60    }
  61    TH1F *piG=hprt[0];
  62    TH1F *piR=hprt[1];
  63    TH1F *piF=hprt[2];
  64    TH1F *piGood=hprt[3];
  65         piGood->SetTitle("Combined PID for pions");
  66         piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
  67    TH1F *piFake=hprt[4];
  68         piFake->SetLineColor(2);
  69    TH1F *kaG=hprt[5];
  70    TH1F *kaR=hprt[6];
  71    TH1F *kaF=hprt[7];
  72    TH1F *kaGood=hprt[8];
  73         kaGood->SetTitle("Combined PID for kaons");
  74         kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
  75    TH1F *kaFake=hprt[9];
  76         kaFake->SetLineColor(2);
  77    TH1F *prG=hprt[10];
  78    TH1F *prR=hprt[11];
  79    TH1F *prF=hprt[12];
  80    TH1F *prGood=hprt[13];
  81         prGood->SetTitle("Combined PID for protons");
  82         prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
  83    TH1F *prFake=hprt[14];
  84         prFake->SetLineColor(2);
  85
  86    delete[] hprt;
  87
  88    Char_t fname[100];
  89
  90    if (gAlice) {
  91       delete gAlice->GetRunLoader();
  92       delete gAlice;
  93       gAlice=0;
  94    }
  95    sprintf(fname,"%s/galice.root",dir);
  96    AliRunLoader *rl = AliRunLoader::Open(fname);
  97    if (rl == 0x0) {
  98       ::Error("AliESDComparison.C","Can not open session !");
  99       return 1;
 100    }
 101    if (rl->LoadgAlice()) {
 102       ::Error("AliESDComparison.C","LoadgAlice returned error !");
 103       delete rl;
 104       return 1;
 105    }
 106    if (rl->LoadHeader()) {
 107       ::Error("AliESDComparison.C","LoadHeader returned error !");
 108       delete rl;
 109       return 1;
 110    }
 111    rl->LoadKinematics();
 112
 113    sprintf(fname,"%s/AliESDs.root",dir);
 114    TFile *ef=TFile::Open(fname);
 115    if (!ef || !ef->IsOpen()) {
 116       ::Error("AliESDComparison.C","Can't AliESDs.root !");
 117       delete rl;
 118       return 1;
 119    }
 120    AliESD* event = new AliESD;
 121    TTree* tree = (TTree*) ef->Get("esdTree");
 122    if (!tree) {
 123       ::Error("AliESDComparison.C", "no ESD tree found");
 124       delete rl;
 125       return 1;
 126    }
 127    tree->SetBranchAddress("ESD", &event);
 128
 129    //****** Tentative particle type "concentrations"
 130    Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
 131
 132    //******* The loop over events
 133    Int_t e=0;
 134    while (tree->GetEvent(e)) {
 135       cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
 136
 137       rl->GetEvent(e);
 138
 139       e++;
 140
 141       Int_t ntrk=event->GetNumberOfTracks();
 142       cerr<<"Number of ESD tracks : "<<ntrk<<endl;
 143
 144       Int_t pisel=0,kasel=0,prsel=0,nsel=0;
 145
 146       AliStack *stack = rl->Stack();
 147
 148       while (ntrk--) {
 149          AliESDtrack *t=event->GetTrack(ntrk);
 150
 151          UInt_t status=AliESDtrack::kESDpid;
 152          status|=AliESDtrack::kITSpid;
 153          status|=AliESDtrack::kTPCpid;
 154          status|=AliESDtrack::kTRDpid;
 155          status|=AliESDtrack::kTOFpid;
 156
 157         if ((t->GetStatus()&status) == status) {
 158            nsel++;
 159
 160            Double_t p=t->GetP();
 161            Double_t dedx=t->GetTPCsignal();
 162            tpcHist->Fill(p,dedx,1);
 163
 164            Int_t lab=TMath::Abs(t->GetLabel());
 165            TParticle *part=stack->Particle(lab);
 166            Int_t code=part->GetPdgCode();
 167
 168            Double_t r[10]; t->GetESDpid(r);
 169            //t->GetTRDpid(r);
 170
 171            Double_t rcc=0.;
 172            Int_t i;
 173            for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
 174            if (rcc==0.) continue;
 175
 176            //Here we apply Bayes' formula
 177            Double_t w[10];
 178            for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
 179
 180            if (TMath::Abs(code)==2212) prR->Fill(p);
 181            if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
 182               prsel++;
 183               prG->Fill(p);
 184               if (TMath::Abs(code)!=2212) prF->Fill(p);
 185            }
 186
 187            if (TMath::Abs(code)==321) kaR->Fill(p);
 188            if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
 189               kasel++;
 190               kaG->Fill(p);
 191               if (TMath::Abs(code)!=321) kaF->Fill(p);
 192            }
 193
 194            if (TMath::Abs(code)==211) piR->Fill(p);
 195            if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
 196               pisel++;
 197               piG->Fill(p);
 198               if (TMath::Abs(code)!=211) piF->Fill(p);
 199            }
 200         }
 201       }
 202       cout<<"Number of selected ESD tracks : "<<nsel<<endl;
 203       cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
 204       cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
 205       cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
 206
 207       allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
 208
 209    } // ***** End of the loop over events
 210
 211    delete event;
 212    ef->Close();
 213
 214    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
 215    if (!c1) {
 216       c1=new TCanvas("c1","",0,0,600,1200);
 217       c1->Divide(1,4);
 218    }
 219
 220    c1->cd(1);
 221    tpcHist->Draw();
 222
 223    c1->cd(2);
 224    //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
 225    piGood->Add(piG,piF,1,-1);
 226    piGood->Divide(piGood,piR,1,1,"b");
 227    piFake->Divide(piF,piG,1,1,"b");
 228    piGood->Draw("hist");
 229    piFake->Draw("same");
 230
 231    c1->cd(3);
 232    //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
 233    kaGood->Add(kaG,kaF,1,-1);
 234    kaGood->Divide(kaGood,kaR,1,1,"b");
 235    kaFake->Divide(kaF,kaG,1,1,"b");
 236    kaGood->Draw("hist");
 237    kaFake->Draw("same");
 238
 239    c1->cd(4);
 240    //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
 241    prGood->Add(prG,prF,1,-1);
 242    prGood->Divide(prGood,prR,1,1,"b");
 243    prFake->Divide(prF,prG,1,1,"b");
 244    prGood->Draw("hist");
 245    prFake->Draw("same");
 246
 247    c1->Update();
 248
 249    cout<<endl;
 250    cout<<endl;
 251    e=(Int_t)piR->GetEntries();
 252    Int_t o=(Int_t)piG->GetEntries();
 253    if (e*o) cout<<"Efficiency (contamination) for pions : "<<
 254       piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
 255    e=(Int_t)kaR->GetEntries();
 256    o=(Int_t)kaG->GetEntries();
 257    if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
 258       kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
 259    e=(Int_t)prR->GetEntries();
 260    o=(Int_t)prG->GetEntries();
 261    if (e*o) cout<<"Efficiency (contamination) for protons : "<<
 262       prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
 263    cout<<endl;
 264
 265    cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
 266    cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
 267    cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
 268    cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
 269
 270    ef->Close();
 271    TFile fc("AliESDComparison.root","RECREATE");
 272    c1->Write();
 273    fc.Close();
 274
 275    gBenchmark->Stop("AliESDComparison");
 276    gBenchmark->Show("AliESDComparison");
 277
 278    delete rl;
 279    return 0;
 280 }