Introducing hyperons in ESD (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>
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
25extern AliRun *gAlice;
26
27Int_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}