Adding track references at decay points (M.Ivanov)
[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 "TKey.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
34    TStopwatch timer;
35    Int_t rc=0,n=0;
36    TKey *key=0;
37    TIter next(ef->GetListOfKeys());
38
39    //****** Tentative particle type "concentrations"
40    Double_t c[5]={0.0, 0.0, 0.1, 0.1, 0.1};
41
42    //******* The loop over events
43    while ((key=(TKey*)next())!=0) {
44
45      cerr<<"Processing event number : "<<n++<<endl;
46
47      AliESD *event=(AliESD*)key->ReadObj();
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      delete event;
78    }
79
80    timer.Stop(); timer.Print();
81
82    hm->Draw();
83
84    ef->Close();
85
86
87    return rc;
88 }