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