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 #ifndef __CINT__
   9   #include <Riostream.h>
  10   #include "TKey.h"
  11   #include "TFile.h"
  12   #include "TH2F.h"
  13   #include "TCanvas.h"
  14   #include "TStopwatch.h"
  15
  16   #include "AliESD.h"
  17 #endif
  18
  19 Int_t AliESDanalysis(Int_t nev=1) {
  20   TH2F *tpcHist=
  21      new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
  22   tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
  23      tpcHist->SetMarkerStyle(8);
  24      tpcHist->SetMarkerSize(0.3);
  25
  26   TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
  27      elHist->SetMarkerStyle(8);
  28      elHist->SetMarkerSize(0.3);
  29
  30   TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
  31      piHist->SetMarkerColor(2);
  32      piHist->SetMarkerStyle(8);
  33      piHist->SetMarkerSize(0.3);
  34
  35   TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
  36      kaHist->SetMarkerColor(4);
  37      kaHist->SetMarkerStyle(8);
  38      kaHist->SetMarkerSize(0.3);
  39
  40   TH2F *prHist=new
  41    TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
  42    prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
  43      prHist->SetMarkerColor(6);
  44      prHist->SetMarkerStyle(8);
  45      prHist->SetMarkerSize(0.3);
  46
  47    TFile *ef=TFile::Open("AliESDs.root");
  48    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
  49
  50    TStopwatch timer;
  51    Int_t rc=0,n=0;
  52    TKey *key=0;
  53    TIter next(ef->GetListOfKeys());
  54
  55    //****** Tentative particle type "concentrations"
  56    Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
  57
  58    //******* The loop over events
  59    while ((key=(TKey*)next())!=0) {
  60      cerr<<"Processing event number : "<<n++<<endl;
  61      AliESD *event=(AliESD*)key->ReadObj();
  62
  63      Int_t ntrk=event->GetNumberOfTracks();
  64      cerr<<"Number of ESD tracks : "<<ntrk<<endl;
  65      //****** The loop over tracks
  66      while (ntrk--) {
  67        AliESDtrack *t=event->GetTrack(ntrk);
  68
  69        Double_t p=t->GetP();
  70
  71        if (t->GetStatus()&AliESDtrack::kTPCin) {
  72          Double_t dedx=t->GetTPCsignal();
  73          tpcHist->Fill(p,dedx,1);
  74        }
  75
  76        if (t->GetStatus()&AliESDtrack::kESDpid) {
  77          Double_t dedx=t->GetTPCsignal();
  78          Double_t r[10]; t->GetESDpid(r);
  79          Double_t rc=0.;
  80          Int_t i;
  81          for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
  82          if (rc==0.) continue;
  83
  84          //Here we apply Bayes' formula
  85          Double_t w[10];
  86          for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
  87
  88          if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
  89          if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
  90          if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
  91          if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
  92        }
  93
  94      }
  95      delete event;
  96    }
  97    timer.Stop(); timer.Print();
  98
  99    TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
 100    c1->Divide(1,2);
 101
 102    c1->cd(1);
 103    tpcHist->Draw();
 104    c1->cd(2);
 105    prHist->Draw();
 106    kaHist->Draw("same");
 107    piHist->Draw("same");
 108    elHist->Draw("same");
 109
 110    ef->Close();
 111
 112    return rc;
 113 }