o) added selector for multiplicity distribution
[u/mrichter/AliRoot.git] / PWG0 / AliPWG0Helper.cxx
1 /* $Id$ */
2
3 #include <AliPWG0Helper.h>
4
5 #include <TParticle.h>
6 #include <TParticlePDG.h>
7 #include <TH1.h>
8 #include <TH3.h>
9
10 #include <AliLog.h>
11 #include <AliESD.h>
12 #include <AliESDVertex.h>
13
14 //____________________________________________________________________
15 ClassImp(AliPWG0Helper)
16
17 //____________________________________________________________________
18 Bool_t AliPWG0Helper::IsEventTriggered(AliESD* aEsd)
19 {
20   // check if the event was triggered
21   //
22   // MB should be
23   // ITS_SPD_GFO_L0  : 32
24   // VZERO_OR_LEFT   : 1
25   // VZERO_OR_RIGHT  : 2
26
27   ULong64_t triggerMask = aEsd->GetTriggerMask();
28
29   if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2)))
30     return kTRUE;
31
32   return kFALSE;
33 }
34
35 //____________________________________________________________________
36 Bool_t AliPWG0Helper::IsVertexReconstructed(AliESD* aEsd)
37 {
38   // checks if the vertex is reasonable
39
40   const AliESDVertex* vtxESD = aEsd->GetVertex();
41
42   // the vertex should be reconstructed
43   if (strcmp(vtxESD->GetName(), "default")==0)
44     return kFALSE;
45
46   Double_t vtx_res[3];
47   vtx_res[0] = vtxESD->GetXRes();
48   vtx_res[1] = vtxESD->GetYRes();
49   vtx_res[2] = vtxESD->GetZRes();
50
51   if (vtx_res[2]==0 || vtx_res[2]>0.1)
52     return kFALSE;
53
54   return kTRUE;
55 }
56
57 //____________________________________________________________________
58 Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t debug)
59 {
60   //
61   // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
62   // shall be counted as a primary particle
63   //
64   // This function or a equivalent should be available in some common place of AliRoot
65   //
66
67   // if the particle has a daughter primary, we do not want to count it
68   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
69   {
70     if (debug)
71       printf("Dropping particle because it has a daughter among the primaries.\n");
72     return kFALSE;
73   }
74
75   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
76
77   // skip quarks and gluon
78   if (pdgCode <= 10 || pdgCode == 21)
79   {
80     if (debug)
81       printf("Dropping particle because it is a quark or gluon.\n");
82     return kFALSE;
83   }
84
85   if (strcmp(aParticle->GetName(),"XXX") == 0)
86   {
87     if (debug)
88       printf("WARNING: There is a particle named XXX.\n");
89     return kFALSE;
90   }
91
92   TParticlePDG* pdgPart = aParticle->GetPDG();
93
94   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
95   {
96     if (debug)
97       printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
98     return kFALSE;
99   }
100
101   if (pdgPart->Charge() == 0)
102   {
103     if (debug)
104       printf("Dropping particle because it is not charged.\n");
105     return kFALSE;
106   }
107
108   return kTRUE;
109 }
110
111 //____________________________________________________________________
112 void AliPWG0Helper::CreateProjections(TH3* hist)
113 {
114   // create projections of 3d hists to all 2d combinations
115   // the histograms are not returned, just use them from memory or use this to create them in a file
116
117   TH1* proj = hist->Project3D("yx");
118   proj->SetXTitle(hist->GetXaxis()->GetTitle());
119   proj->SetYTitle(hist->GetYaxis()->GetTitle());
120
121   proj = hist->Project3D("zx");
122   proj->SetXTitle(hist->GetXaxis()->GetTitle());
123   proj->SetYTitle(hist->GetZaxis()->GetTitle());
124
125   proj = hist->Project3D("zy");
126   proj->SetXTitle(hist->GetYaxis()->GetTitle());
127   proj->SetYTitle(hist->GetZaxis()->GetTitle());
128 }
129
130 //____________________________________________________________________
131 void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors)
132 {
133   // create projections of the 3d hists divides them
134   // axis decides to which plane, if axis is 0 to all planes
135   // the histograms are not returned, just use them from memory or use this to create them in a file
136
137   if (axis == 0)
138   {
139     CreateDividedProjections(hist, hist2, "yx", putErrors);
140     CreateDividedProjections(hist, hist2, "zx", putErrors);
141     CreateDividedProjections(hist, hist2, "zy", putErrors);
142
143     return;
144   }
145
146   TH1* proj = hist->Project3D(axis);
147
148   if (strlen(axis) == 2)
149   {
150     proj->SetYTitle(GetAxisTitle(hist, axis[0]));
151     proj->SetXTitle(GetAxisTitle(hist, axis[1]));
152   }
153   else if (strlen(axis) == 1)
154     proj->SetXTitle(GetAxisTitle(hist, axis[0]));
155
156   TH1* proj2 = hist2->Project3D(axis);
157   if (strlen(axis) == 2)
158   {
159     proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
160     proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
161   }
162   else if (strlen(axis) == 1)
163     proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
164
165   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
166   division->Divide(proj2);
167
168   if (putErrors)
169   {
170     division->Sumw2();
171     if (division->GetDimension() == 1)
172     {
173       Int_t nBins = division->GetNbinsX();
174       for (Int_t i = 0; i <= nBins; ++i)
175         if (proj2->GetBinContent(i) != 0)
176           division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
177     }
178     else if (division->GetDimension() == 2)
179     {
180       Int_t nBinsX = division->GetNbinsX();
181       Int_t nBinsY = division->GetNbinsY();
182       for (Int_t i = 0; i <= nBinsX; ++i)
183         for (Int_t j = 0; j <= nBinsY; ++j)
184           if (proj2->GetBinContent(i, j) != 0)
185             division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
186     }
187   }
188 }
189
190 //____________________________________________________________________
191 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
192 {
193   // returns the title of the axis given in axis (x, y, z)
194
195   if (axis == 'x')
196     return hist->GetXaxis()->GetTitle();
197   else if (axis == 'y')
198     return hist->GetYaxis()->GetTitle();
199   else if (axis == 'z')
200     return hist->GetZaxis()->GetTitle();
201
202   return 0;
203 }