Access to the stack using AliRunLoader (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDanalysis.C
index 7593cb6958706928fc49f022353e17c4cda28196..9a990c736f9c195b71b204bfff2252cb8f8d34d7 100644 (file)
@@ -5,45 +5,77 @@
 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 //********************************************************************
 
-#ifndef __CINT__
+#if !defined( __CINT__) || defined(__MAKECINT__)
   #include <Riostream.h>
   #include "TKey.h"
   #include "TFile.h"
+  #include "TH1F.h"
   #include "TH2F.h"
   #include "TCanvas.h"
   #include "TStopwatch.h"
+  #include "TParticle.h"
+
+  #include "AliRun.h"
+  #include "AliStack.h"
+  #include "AliRunLoader.h"
+  #include "AliLoader.h"
 
   #include "AliESD.h"
 #endif
 
+extern AliRun *gAlice;
+
 Int_t AliESDanalysis(Int_t nev=1) { 
   TH2F *tpcHist=
      new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
-  tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
+   tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
      tpcHist->SetMarkerStyle(8); 
      tpcHist->SetMarkerSize(0.3);
  
-  TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
-     elHist->SetMarkerStyle(8); 
-     elHist->SetMarkerSize(0.3);
-
-  TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
-     piHist->SetMarkerColor(2); 
-     piHist->SetMarkerStyle(8); 
-     piHist->SetMarkerSize(0.3);
-
-  TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
-     kaHist->SetMarkerColor(4); 
-     kaHist->SetMarkerStyle(8); 
-     kaHist->SetMarkerSize(0.3);
-
-  TH2F *prHist=new 
-   TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
-   prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
-     prHist->SetMarkerColor(6); 
-     prHist->SetMarkerStyle(8); 
-     prHist->SetMarkerSize(0.3);
-   
+  TH1F *piG=new TH1F("piG","",20,0.,4.);
+  TH1F *piR=new TH1F("piR","",20,0.,4.);
+  TH1F *piF=new TH1F("piF","",20,0.,4.);
+  TH1F *piGood=new TH1F("piGood","Combined PID for pions",20,0.,4.); 
+  piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
+  TH1F *piFake=new TH1F("piFake","",20,0.,4.); piFake->SetLineColor(2);
+
+  TH1F *kaG=new TH1F("kaG","",20,0.,4.);
+  TH1F *kaR=new TH1F("kaR","",20,0.,4.);
+  TH1F *kaF=new TH1F("kaF","",20,0.,4.);
+  TH1F *kaGood=new TH1F("kaGood","Combined PID for kaons",20,0.,4.); 
+  kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
+  TH1F *kaFake=new TH1F("kaFake","",20,0.,4.); kaFake->SetLineColor(2);
+
+  TH1F *prG=new TH1F("prG","",20,0.,4.);
+  TH1F *prR=new TH1F("prR","",20,0.,4.);
+  TH1F *prF=new TH1F("prF","",20,0.,4.);
+  TH1F *prGood=new TH1F("prGood","Combined PID for protons",20,0.,4.); 
+  prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
+  TH1F *prFake=new TH1F("prFake","",20,0.,4.); prFake->SetLineColor(2);
+
+   if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice;
+      gAlice=0;
+   }
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (rl == 0x0) {
+      cerr<<"Can not open session"<<endl;
+      return 1;
+   }
+   if (rl->LoadgAlice()) {
+      cerr<<"LoadgAlice returned error"<<endl;
+      delete rl;
+      return 1;
+   }
+   if (rl->LoadHeader()) {
+      cerr<<"LoadHeader returned error"<<endl;
+      delete rl;
+      return 1;
+   }
+   rl->LoadKinematics();
+   AliStack* stack = rl->Stack();
+
    TFile *ef=TFile::Open("AliESDs.root");
    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
 
@@ -53,16 +85,22 @@ Int_t AliESDanalysis(Int_t nev=1) {
    TIter next(ef->GetListOfKeys());
 
    //****** Tentative particle type "concentrations"
-   Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
+   Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
 
    //******* The loop over events
    while ((key=(TKey*)next())!=0) {
+     rl->GetEvent(n);
+
      cerr<<"Processing event number : "<<n++<<endl;
+
      AliESD *event=(AliESD*)key->ReadObj();
 
      Int_t ntrk=event->GetNumberOfTracks();
      cerr<<"Number of ESD tracks : "<<ntrk<<endl; 
      //****** The loop over tracks
+
+     Int_t pisel=0,kasel=0,prsel=0,nsel=0;
+
      while (ntrk--) {
        AliESDtrack *t=event->GetTrack(ntrk);
 
@@ -73,41 +111,91 @@ Int_t AliESDanalysis(Int_t nev=1) {
          tpcHist->Fill(p,dedx,1);
        }
 
-       if (t->GetStatus()&AliESDtrack::kESDpid) {
-         Double_t dedx=t->GetTPCsignal();
-        Double_t r[10]; t->GetESDpid(r);
-         Double_t rc=0.;
+       UInt_t status=AliESDtrack::kESDpid;
+       status|=AliESDtrack::kTPCpid; 
+       status|=AliESDtrack::kTOFpid; 
+
+       if ((t->GetStatus()&status) == status) {
+         nsel++;
+
+         Int_t lab=TMath::Abs(t->GetLabel());
+         TParticle *part=stack->Particle(lab);
+         Int_t code=part->GetPdgCode();
+
+         Double_t r[10]; t->GetESDpid(r);
+
+         Double_t rcc=0.;
          Int_t i;
-         for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
-         if (rc==0.) continue;
+         for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
+         if (rcc==0.) continue;
 
         //Here we apply Bayes' formula
          Double_t w[10];
-         for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
+         for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
+
+         if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
+           prsel++;
+           prG->Fill(p);
+            if (TMath::Abs(code)==2212) prR->Fill(p);
+            else prF->Fill(p);
+         }
+
+         if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
+           kasel++;
+           kaG->Fill(p);
+            if (TMath::Abs(code)==321) kaR->Fill(p);
+            else kaF->Fill(p);
+         }
+
+        if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
+           pisel++;
+           piG->Fill(p);
+           if (TMath::Abs(code)==211) piR->Fill(p);
+            else piF->Fill(p);
+         }
 
-         if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
-         if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
-         if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
-         if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
        }
 
      } 
      delete event;
+     cerr<<"Number of selected ESD tracks : "<<nsel<<endl;
+     cerr<<"Number of selected pion ESD tracks : "<<pisel<<endl;
+     cerr<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
+     cerr<<"Number of selected proton ESD tracks : "<<prsel<<endl;
    }
+
    timer.Stop(); timer.Print();
 
    TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
-   c1->Divide(1,2);
+   c1->Divide(1,4);
 
    c1->cd(1);
    tpcHist->Draw();
+
    c1->cd(2);
-   prHist->Draw();
-   kaHist->Draw("same");
-   piHist->Draw("same");
-   elHist->Draw("same");
+   piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
+   piGood->Divide(piR,piG,1,1,"b");
+   piFake->Divide(piF,piG,1,1,"b");
+   piGood->Draw("hist");
+   piFake->Draw("same");
+
+   c1->cd(3);
+   kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
+   kaGood->Divide(kaR,kaG,1,1,"b");
+   kaFake->Divide(kaF,kaG,1,1,"b");
+   kaGood->Draw("hist");
+   kaFake->Draw("same");
+
+   c1->cd(4);
+   prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
+   prGood->Divide(prR,prG,1,1,"b");
+   prFake->Divide(prF,prG,1,1,"b");
+   prGood->Draw("hist");
+   prFake->Draw("same");
 
    ef->Close();
 
+   delete rl;
+
    return rc;
 }