Fixes for coverity.
[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        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 }