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::kTOFpid;
 155
 156         if ((t->GetStatus()&status) == status) {
 157            nsel++;
 158
 159            Double_t p=t->GetP();
 160            Double_t dedx=t->GetTPCsignal();
 161            tpcHist->Fill(p,dedx,1);
 162
 163            Int_t lab=TMath::Abs(t->GetLabel());
 164            TParticle *part=stack->Particle(lab);
 165            Int_t code=part->GetPdgCode();
 166
 167            Double_t r[10]; t->GetESDpid(r);
 168
 169            Double_t rcc=0.;
 170            Int_t i;
 171            for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
 172            if (rcc==0.) continue;
 173
 174            //Here we apply Bayes' formula
 175            Double_t w[10];
 176            for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
 177
 178            if (TMath::Abs(code)==2212) prR->Fill(p);
 179            if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
 180               prsel++;
 181               prG->Fill(p);
 182               if (TMath::Abs(code)!=2212) prF->Fill(p);
 183            }
 184
 185            if (TMath::Abs(code)==321) kaR->Fill(p);
 186            if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
 187               kasel++;
 188               kaG->Fill(p);
 189               if (TMath::Abs(code)!=321) kaF->Fill(p);
 190            }
 191
 192            if (TMath::Abs(code)==211) piR->Fill(p);
 193            if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
 194               pisel++;
 195               piG->Fill(p);
 196               if (TMath::Abs(code)!=211) piF->Fill(p);
 197            }
 198         }
 199       }
 200       cout<<"Number of selected ESD tracks : "<<nsel<<endl;
 201       cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
 202       cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
 203       cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
 204
 205       allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
 206
 207    } // ***** End of the loop over events
 208
 209    delete event;
 210    ef->Close();
 211
 212    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
 213    if (!c1) {
 214       c1=new TCanvas("c1","",0,0,600,1200);
 215       c1->Divide(1,4);
 216    }
 217
 218    c1->cd(1);
 219    tpcHist->Draw();
 220
 221    c1->cd(2);
 222    //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
 223    piGood->Add(piG,piF,1,-1);
 224    piGood->Divide(piGood,piR,1,1,"b");
 225    piFake->Divide(piF,piG,1,1,"b");
 226    piGood->Draw("hist");
 227    piFake->Draw("same");
 228
 229    c1->cd(3);
 230    //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
 231    kaGood->Add(kaG,kaF,1,-1);
 232    kaGood->Divide(kaGood,kaR,1,1,"b");
 233    kaFake->Divide(kaF,kaG,1,1,"b");
 234    kaGood->Draw("hist");
 235    kaFake->Draw("same");
 236
 237    c1->cd(4);
 238    //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
 239    prGood->Add(prG,prF,1,-1);
 240    prGood->Divide(prGood,prR,1,1,"b");
 241    prFake->Divide(prF,prG,1,1,"b");
 242    prGood->Draw("hist");
 243    prFake->Draw("same");
 244
 245    c1->Update();
 246
 247    cout<<endl;
 248    cout<<endl;
 249    e=(Int_t)piR->GetEntries();
 250    Int_t o=(Int_t)piG->GetEntries();
 251    if (e*o) cout<<"Efficiency (contamination) for pions : "<<
 252       piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
 253    e=(Int_t)kaR->GetEntries();
 254    o=(Int_t)kaG->GetEntries();
 255    if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
 256       kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
 257    e=(Int_t)prR->GetEntries();
 258    o=(Int_t)prG->GetEntries();
 259    if (e*o) cout<<"Efficiency (contamination) for protons : "<<
 260       prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
 261    cout<<endl;
 262
 263    cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
 264    cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
 265    cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
 266    cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
 267
 268    ef->Close();
 269    TFile fc("AliESDComparison.root","RECREATE");
 270    c1->Write();
 271    fc.Close();
 272
 273    gBenchmark->Stop("AliESDComparison");
 274    gBenchmark->Show("AliESDComparison");
 275
 276    delete rl;
 277    return 0;
 278 }