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 |
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"; |
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 | } |