bf2dd7d2c7e648619da495c95bcb3c6ebb0e0792
[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 "TH2F.h"
15   #include "TCanvas.h"
16   #include "TStopwatch.h"
17   #include "TParticle.h"
18
19   #include "AliRun.h"
20
21   #include "AliESD.h"
22
23 #endif
24
25 extern AliRun *gAlice;
26
27 Int_t AliESDv0Analysis(Int_t nev=1) { 
28    TH1F *hm=new TH1F("hm","Effective Mass",40,1.065,1.165);
29    hm->SetXTitle("Mass (GeV/c**2)");
30
31    TFile *ef=TFile::Open("AliESDs.root");
32    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
33    AliESD* event = new AliESD;
34    TTree* tree = (TTree*) ef->Get("esdTree");
35    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
36    tree->SetBranchAddress("ESD", &event);
37
38    TStopwatch timer;
39    Int_t rc=0,n=0;
40
41    //****** Tentative particle type "concentrations"
42    Double_t c[5]={0.0, 0.0, 0.1, 0.1, 0.1};
43
44    //******* The loop over events
45    while (tree->GetEvent(n)) {
46
47      cerr<<"Processing event number : "<<n++<<endl;
48
49      Int_t nv0=event->GetNumberOfV0s();
50      cerr<<"Number of ESD v0s : "<<nv0<<endl; 
51
52      while (nv0--) {
53        AliESDv0 *v0=event->GetV0(nv0);
54        Int_t pi=v0->GetPindex();
55        AliESDtrack *t=event->GetTrack(pi);
56        Int_t isProton=1;
57        if ((t->GetStatus()&AliESDtrack::kESDpid)!=0) {
58          Double_t r[10]; t->GetESDpid(r);
59          Double_t rcc=0.;
60          Int_t i;
61          for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
62          if (rcc==0.) continue;
63          //Here we apply Bayes' formula
64          Double_t w[10];
65          for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
66
67          if (w[4]<w[3]) isProton=0;
68          if (w[4]<w[2]) isProton=0;
69          if (w[4]<w[1]) isProton=0;
70          if (w[4]<w[0]) isProton=0;
71        }
72        if (!isProton) continue;
73        v0->ChangeMassHypothesis(3122);
74        Double_t mass=v0->GetEffMass();
75        hm->Fill(mass);
76      } 
77    }
78
79    delete event;
80    ef->Close();
81
82    timer.Stop(); timer.Print();
83
84    hm->Draw();
85
86    ef->Close();
87
88
89    return rc;
90 }