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 <AliESDEvent.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    AliESDEvent* event = new AliESDEvent();
 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    event->ReadFromTree(tree);
 128
 129    //****** Tentative particle type "concentrations"
 130    Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
 131    //Double_t c[5]={0.2, 0.2, 0.2, 0.2, 0.2};
 132    AliPID::SetPriors(c);
 133
 134
 135    //******* The loop over events
 136    Int_t e=0;
 137    while (tree->GetEvent(e)) {
 138       cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
 139
 140       rl->GetEvent(e);
 141
 142       e++;
 143
 144       Int_t ntrk=event->GetNumberOfTracks();
 145       cerr<<"Number of ESD tracks : "<<ntrk<<endl;
 146
 147       Int_t pisel=0,kasel=0,prsel=0,nsel=0;
 148
 149       AliStack *stack = rl->Stack();
 150
 151       while (ntrk--) {
 152          AliESDtrack *t=event->GetTrack(ntrk);
 153
 154          UInt_t status=AliESDtrack::kESDpid;
 155          status|=AliESDtrack::kITSpid;
 156          status|=AliESDtrack::kTPCpid;
 157          status|=AliESDtrack::kTRDpid;
 158          status|=AliESDtrack::kTOFpid;
 159
 160         if ((t->GetStatus()&status) == status) {
 161            nsel++;
 162
 163            Double_t p=t->GetP();
 164            Double_t dedx=t->GetTPCsignal();
 165            tpcHist->Fill(p,dedx,1);
 166
 167            Int_t lab=TMath::Abs(t->GetLabel());
 168            TParticle *part=stack->Particle(lab);
 169            Int_t code=part->GetPdgCode();
 170
 171            Double_t r[10]; t->GetESDpid(r);
 172            //t->GetTRDpid(r);
 173            //t->GetTPCpid(r);
 174
 175            AliPID pid(r);
 176
 177            Double_t w[10];
 178            w[0]=pid.GetProbability(AliPID::kElectron);
 179            w[1]=pid.GetProbability(AliPID::kMuon);
 180            w[2]=pid.GetProbability(AliPID::kPion);
 181            w[3]=pid.GetProbability(AliPID::kKaon);
 182            w[4]=pid.GetProbability(AliPID::kProton);
 183
 184
 185            if (TMath::Abs(code)==2212) prR->Fill(p);
 186            if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
 187               prsel++;
 188               prG->Fill(p);
 189               if (TMath::Abs(code)!=2212) prF->Fill(p);
 190            }
 191
 192            if (TMath::Abs(code)==321) kaR->Fill(p);
 193            if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
 194               kasel++;
 195               kaG->Fill(p);
 196               if (TMath::Abs(code)!=321) kaF->Fill(p);
 197            }
 198
 199            if (TMath::Abs(code)==211) piR->Fill(p);
 200            if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion
 201               pisel++;
 202               piG->Fill(p);
 203               if (TMath::Abs(code)!=211) piF->Fill(p);
 204            }
 205         }
 206       }
 207       cout<<"Number of selected ESD tracks : "<<nsel<<endl;
 208       cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
 209       cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
 210       cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
 211
 212       allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
 213
 214    } // ***** End of the loop over events
 215
 216    delete event;
 217    delete tree;
 218    ef->Close();
 219
 220    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
 221    if (!c1) {
 222       c1=new TCanvas("c1","",0,0,600,1200);
 223       c1->Divide(1,4);
 224    }
 225
 226    c1->cd(1);
 227    tpcHist->Draw();
 228
 229    c1->cd(2);
 230    //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
 231    piGood->Add(piG,piF,1,-1);
 232    piGood->Divide(piGood,piR,1,1,"b");
 233    piFake->Divide(piF,piG,1,1,"b");
 234    piGood->Draw("hist");
 235    piFake->Draw("same");
 236
 237    c1->cd(3);
 238    //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
 239    kaGood->Add(kaG,kaF,1,-1);
 240    kaGood->Divide(kaGood,kaR,1,1,"b");
 241    kaFake->Divide(kaF,kaG,1,1,"b");
 242    kaGood->Draw("hist");
 243    kaFake->Draw("same");
 244
 245    c1->cd(4);
 246    //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
 247    prGood->Add(prG,prF,1,-1);
 248    prGood->Divide(prGood,prR,1,1,"b");
 249    prFake->Divide(prF,prG,1,1,"b");
 250    prGood->Draw("hist");
 251    prFake->Draw("same");
 252
 253    c1->Update();
 254
 255    cout<<endl;
 256    cout<<endl;
 257    e=(Int_t)piR->GetEntries();
 258    Int_t o=(Int_t)piG->GetEntries();
 259    if (e*o) cout<<"Efficiency (contamination) for pions : "<<
 260       piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
 261    e=(Int_t)kaR->GetEntries();
 262    o=(Int_t)kaG->GetEntries();
 263    if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
 264       kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
 265    e=(Int_t)prR->GetEntries();
 266    o=(Int_t)prG->GetEntries();
 267    if (e*o) cout<<"Efficiency (contamination) for protons : "<<
 268       prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
 269    cout<<endl;
 270
 271    cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
 272    cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
 273    cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
 274    cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
 275
 276    ef->Close();
 277    TFile fc("AliESDComparison.root","RECREATE");
 278    c1->Write();
 279    fc.Close();
 280
 281    gBenchmark->Stop("AliESDComparison");
 282    gBenchmark->Show("AliESDComparison");
 283
 284    delete rl;
 285    return 0;
 286 }