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 <TROOT.h>
  12   #include <TTree.h>
  13   #include <TFile.h>
  14   #include <TH1F.h>
  15   #include <TCanvas.h>
  16
  17   #include "AliESDEvent.h"
  18   #include "AliESDv0.h"
  19 #endif
  20
  21 extern TROOT *gROOT;
  22
  23 Int_t AliESDv0Analysis(const Char_t *dir=".") {
  24    TH1F *hm=(TH1F*)gROOT->FindObject("hm");
  25    if (!hm) {
  26       hm=new TH1F("hm","Lambda+LambdaBar Effective Mass",60,1.065,1.165);
  27       hm->SetXTitle("Mass (GeV/c**2)");
  28    }
  29    Char_t fname[100];
  30    sprintf(fname,"%s/AliESDs.root",dir);
  31    TFile *ef=TFile::Open(fname);
  32    if (!ef||!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
  33    cerr<<"\n****** "<<fname<<" ******\n";
  34
  35    AliESDEvent* event = new AliESDEvent();
  36
  37
  38    TTree* tree = (TTree*) ef->Get("esdTree");
  39    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
  40    event->ReadFromTree(tree);
  41
  42    Int_t rc=0,n=0;
  43
  44    //****** Tentative particle type "concentrations"
  45    Double_t c[5]={0.0, 0.0, 1, 0, 1};
  46    AliPID pid;
  47    pid.SetPriors(c);
  48
  49    //******* The loop over events
  50     while (tree->GetEvent(n)) {
  51
  52      cerr<<"Processing event number : "<<n++<<endl;
  53
  54      Int_t nv0=event->GetNumberOfV0s();
  55      cerr<<"Number of ESD v0s : "<<nv0<<endl;
  56
  57      while (nv0--) {
  58        AliESDv0 *v0=event->GetV0(nv0);
  59        if (v0->GetOnFlyStatus()) continue;
  60
  61        Int_t protonIdx=v0->GetPindex();
  62        Int_t pionIdx  =v0->GetNindex();
  63
  64        v0->ChangeMassHypothesis(3122);
  65        Double_t mass=v0->GetEffMass();
  66        if (mass>1.17) {  //check also the LambdaBar hypothesis
  67           v0->ChangeMassHypothesis(-3122);
  68           mass=v0->GetEffMass();
  69           if (mass>1.17) continue;
  70           Int_t tmp=protonIdx; protonIdx=pionIdx; pionIdx=tmp;
  71        }
  72
  73        AliESDtrack *protonTrk=event->GetTrack(protonIdx);
  74        AliESDtrack *pionTrk  =event->GetTrack(pionIdx);
  75
  76        if (protonTrk->GetP()<0.5) continue;
  77
  78        // Check if the "proton track" is a proton
  79        if ((protonTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
  80          Double_t r[10]; protonTrk->GetESDpid(r);
  81          pid.SetProbabilities(r);
  82          Double_t pp=pid.GetProbability(AliPID::kProton);
  83          if (pp < pid.GetProbability(AliPID::kElectron)) continue;
  84          if (pp < pid.GetProbability(AliPID::kMuon)) continue;
  85          if (pp < pid.GetProbability(AliPID::kPion)) continue;
  86          if (pp < pid.GetProbability(AliPID::kKaon)) continue;
  87        }
  88
  89        //Check if the "pion track" is a pion
  90        if ((pionTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
  91          Double_t r[10]; pionTrk->GetESDpid(r);
  92          pid.SetProbabilities(r);
  93          Double_t ppi=pid.GetProbability(AliPID::kPion);
  94          if (ppi < pid.GetProbability(AliPID::kElectron)) continue;
  95          if (ppi < pid.GetProbability(AliPID::kMuon)) continue;
  96          if (ppi < pid.GetProbability(AliPID::kKaon)) continue;
  97          if (ppi < pid.GetProbability(AliPID::kProton)) continue;
  98         }
  99
 100        hm->Fill(mass);
 101      }
 102
 103    }
 104
 105    delete event;
 106    delete tree;
 107    ef->Close();
 108
 109    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
 110    if (!c1) {
 111       c1=new TCanvas();
 112    }
 113
 114    if (hm->GetEntries()>100) hm->Fit("gaus","","",1.11,1.12);
 115    else hm->Draw();
 116
 117    c1->Update();
 118
 119    return rc;
 120 }