102b0586b407ba47950616de3a8b39f715ff1d72
[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 <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 }