STEER/AliESDv0Analysis.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 //            applied to the Lambda0 reconstruction.
   6 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
   7 //********************************************************************
   8
   9 #if !defined( __CINT__) || defined(__MAKECINT__)
  10   #include <Riostream.h>
  11   #include <TTree.h>
  12   #include "TFile.h"
  13   #include "TH1F.h"
  14   #include "TH2F.h"
  15   #include "TCanvas.h"
  16   #include "TStopwatch.h"
  17   #include "TParticle.h"
  18
  19   #include "AliRun.h"
  20
  21   #include "AliESD.h"
  22
  23 #endif
  24
  25 extern AliRun *gAlice;
  26
  27 Int_t AliESDv0Analysis(Int_t nev=1) {
  28    TH1F *hm=new TH1F("hm","Effective Mass",40,1.065,1.165);
  29    hm->SetXTitle("Mass (GeV/c**2)");
  30
  31    TFile *ef=TFile::Open("AliESDs.root");
  32    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
  33    AliESD* event = new AliESD;
  34    TTree* tree = (TTree*) ef->Get("esdTree");
  35    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
  36    tree->SetBranchAddress("ESD", &event);
  37
  38    TStopwatch timer;
  39    Int_t rc=0,n=0;
  40
  41    //****** Tentative particle type "concentrations"
  42    Double_t c[5]={0.0, 0.0, 0.1, 0.1, 0.1};
  43
  44    //******* The loop over events
  45    while (tree->GetEvent(n)) {
  46
  47      cerr<<"Processing event number : "<<n++<<endl;
  48
  49      Int_t nv0=event->GetNumberOfV0s();
  50      cerr<<"Number of ESD v0s : "<<nv0<<endl;
  51
  52      while (nv0--) {
  53        AliESDv0 *v0=event->GetV0(nv0);
  54        Int_t pi=v0->GetPindex();
  55        AliESDtrack *t=event->GetTrack(pi);
  56        Int_t isProton=1;
  57        if ((t->GetStatus()&AliESDtrack::kESDpid)!=0) {
  58          Double_t r[10]; t->GetESDpid(r);
  59          Double_t rcc=0.;
  60          Int_t i;
  61          for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
  62          if (rcc==0.) continue;
  63          //Here we apply Bayes' formula
  64          Double_t w[10];
  65          for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
  66
  67          if (w[4]<w[3]) isProton=0;
  68          if (w[4]<w[2]) isProton=0;
  69          if (w[4]<w[1]) isProton=0;
  70          if (w[4]<w[0]) isProton=0;
  71        }
  72        if (!isProton) continue;
  73        v0->ChangeMassHypothesis(3122);
  74        Double_t mass=v0->GetEffMass();
  75        hm->Fill(mass);
  76      }
  77    }
  78
  79    delete event;
  80    ef->Close();
  81
  82    timer.Stop(); timer.Print();
  83
  84    hm->Draw();
  85
  86    ef->Close();
  87
  88
  89    return rc;
  90 }