updated drawplots macro
[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)
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");
140     CreateDividedProjections(hist, hist2, "zx");
141     CreateDividedProjections(hist, hist2, "zy");
142
143     return;
144   }
145
146   TH1* proj = hist->Project3D(axis);
147   proj->SetYTitle(GetAxisTitle(hist, axis[0]));
148   proj->SetXTitle(GetAxisTitle(hist, axis[1]));
149
150   TH1* proj2 = hist2->Project3D(axis);
151   proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
152   proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
153
154   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
155   division->Divide(proj2);
156 }
157
158 //____________________________________________________________________
159 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
160 {
161   // returns the title of the axis given in axis (x, y, z)
162
163   if (axis == 'x')
164     return hist->GetXaxis()->GetTitle();
165   else if (axis == 'y')
166     return hist->GetYaxis()->GetTitle();
167   else if (axis == 'z')
168     return hist->GetZaxis()->GetTitle();
169
170   return 0;
171 }