Improved version of V0 analysis (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDv0Analysis.C
CommitLineData
e23730c7 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>
8b462fd8 11 #include <TTree.h>
9d0d65be 12 #include <TFile.h>
13 #include <TH1F.h>
14 #include <TCanvas.h>
e23730c7 15
16 #include "AliESD.h"
17
18#endif
19
9d0d65be 20Int_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";
e23730c7 31
8b462fd8 32 AliESD* event = new AliESD;
9d0d65be 33
34
8b462fd8 35 TTree* tree = (TTree*) ef->Get("esdTree");
36 if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
37 tree->SetBranchAddress("ESD", &event);
e23730c7 38
e23730c7 39 Int_t rc=0,n=0;
e23730c7 40
41 //****** Tentative particle type "concentrations"
9d0d65be 42 Double_t c[5]={0.0, 0.0, 1, 0, 1};
e23730c7 43
44 //******* The loop over events
9d0d65be 45 while (tree->GetEvent(n))
46 {
e23730c7 47
48 cerr<<"Processing event number : "<<n++<<endl;
49
e23730c7 50 Int_t nv0=event->GetNumberOfV0s();
51 cerr<<"Number of ESD v0s : "<<nv0<<endl;
52
53 while (nv0--) {
54 AliESDv0 *v0=event->GetV0(nv0);
9d0d65be 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);
e23730c7 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
9d0d65be 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;
e23730c7 88 }
9d0d65be 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
e23730c7 107 hm->Fill(mass);
108 }
9d0d65be 109
e23730c7 110 }
111
8b462fd8 112 delete event;
113 ef->Close();
114
9d0d65be 115 TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
116 if (!c1) {
117 c1=new TCanvas();
118 }
e23730c7 119
9d0d65be 120 if (hm->GetEntries()>100) hm->Fit("gaus","","",1.11,1.12);
121 else hm->Draw();
e23730c7 122
9d0d65be 123 c1->Update();
e23730c7 124
125 return rc;
126}