Bug corrected.
[u/mrichter/AliRoot.git] / 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
60        Int_t protonIdx=v0->GetPindex();
61        Int_t pionIdx  =v0->GetNindex();
62       
63        v0->ChangeMassHypothesis(3122);
64        Double_t mass=v0->GetEffMass();
65        if (mass>1.17) {  //check also the LambdaBar hypothesis
66           v0->ChangeMassHypothesis(-3122);
67           mass=v0->GetEffMass();
68           if (mass>1.17) continue;
69           Int_t tmp=protonIdx; protonIdx=pionIdx; pionIdx=tmp;
70        } 
71
72        AliESDtrack *protonTrk=event->GetTrack(protonIdx);
73        AliESDtrack *pionTrk  =event->GetTrack(pionIdx);
74
75        if (protonTrk->GetP()<0.5) continue;
76
77        // Check if the "proton track" is a proton
78        if ((protonTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
79          Double_t r[10]; protonTrk->GetESDpid(r);
80          pid.SetProbabilities(r);
81          Double_t pp=pid.GetProbability(AliPID::kProton);
82          if (pp < pid.GetProbability(AliPID::kElectron)) continue;
83          if (pp < pid.GetProbability(AliPID::kMuon)) continue;
84          if (pp < pid.GetProbability(AliPID::kPion)) continue;
85          if (pp < pid.GetProbability(AliPID::kKaon)) continue;
86        }
87  
88        //Check if the "pion track" is a pion
89        if ((pionTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
90          Double_t r[10]; pionTrk->GetESDpid(r);
91          pid.SetProbabilities(r);
92          Double_t ppi=pid.GetProbability(AliPID::kPion);
93          if (ppi < pid.GetProbability(AliPID::kElectron)) continue;
94          if (ppi < pid.GetProbability(AliPID::kMuon)) continue;
95          if (ppi < pid.GetProbability(AliPID::kKaon)) continue;
96          if (ppi < pid.GetProbability(AliPID::kProton)) continue;
97         }
98
99        hm->Fill(mass);
100      } 
101
102    }
103
104    delete event;
105    delete tree;
106    ef->Close();
107
108    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
109    if (!c1) {
110       c1=new TCanvas();
111    }
112
113    if (hm->GetEntries()>100) hm->Fit("gaus","","",1.11,1.12);
114    else hm->Draw();
115
116    c1->Update();
117
118    return rc;
119 }