a lot of work on the analysis
[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 <TH3F.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)
59 {
60   //
61   // Returns if the given particle is a primary particle
62   // This function or a equivalent should be available in some common place of AliRoot
63   //
64
65   // if the particle has a daughter primary, we do not want to count it
66   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
67   {
68     //AliDebug(AliLog::kDebug+1, "Dropping particle because it has a daughter among the primaries.");
69     return kFALSE;
70   }
71
72   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
73
74   // skip quarks and gluon
75   if (pdgCode <= 10 || pdgCode == 21)
76   {
77     //AliDebug(AliLog::kDebug+1, "Dropping particle because it is a quark or gluon.");
78     return kFALSE;
79   }
80
81   if (strcmp(aParticle->GetName(),"XXX") == 0)
82   {
83     //AliDebug(AliLog::kDebug, Form("WARNING: There is a particle named XXX."));
84     return kFALSE;
85   }
86
87   TParticlePDG* pdgPart = aParticle->GetPDG();
88
89   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
90   {
91     //AliDebug(AliLog::kDebug, Form("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode));
92     return kFALSE;
93   }
94
95   if (pdgPart->Charge() == 0)
96   {
97     //AliDebug(AliLog::kDebug+1, "Dropping particle because it is not charged.");
98     return kFALSE;
99   }
100
101   return kTRUE;
102 }
103
104 //____________________________________________________________________
105 void AliPWG0Helper::CreateProjections(TH3F* hist)
106 {
107   // create projections of 3d hists to all 2d combinations
108   // the histograms are not returned, just use them from memory or use this to create them in a file
109
110   TH1* proj = hist->Project3D("yx");
111   proj->SetXTitle(hist->GetXaxis()->GetTitle());
112   proj->SetYTitle(hist->GetYaxis()->GetTitle());
113
114   proj = hist->Project3D("zx");
115   proj->SetXTitle(hist->GetXaxis()->GetTitle());
116   proj->SetYTitle(hist->GetZaxis()->GetTitle());
117
118   proj = hist->Project3D("zy");
119   proj->SetXTitle(hist->GetYaxis()->GetTitle());
120   proj->SetYTitle(hist->GetZaxis()->GetTitle());
121 }
122
123 //____________________________________________________________________
124 void AliPWG0Helper::CreateDividedProjections(TH3F* hist, TH3F* hist2, const char* axis)
125 {
126   // create projections of the 3d hists divides them
127   // axis decides to which plane, if axis is 0 to all planes
128   // the histograms are not returned, just use them from memory or use this to create them in a file
129
130   if (axis == 0)
131   {
132     CreateDividedProjections(hist, hist2, "yx");
133     CreateDividedProjections(hist, hist2, "zx");
134     CreateDividedProjections(hist, hist2, "zy");
135
136     return;
137   }
138
139   TH1* proj = hist->Project3D(axis);
140   proj->SetXTitle(hist->GetXaxis()->GetTitle());
141   proj->SetYTitle(hist->GetYaxis()->GetTitle());
142
143   TH1* proj2 = hist2->Project3D(axis);
144   proj2->SetXTitle(hist2->GetXaxis()->GetTitle());
145   proj2->SetYTitle(hist2->GetYaxis()->GetTitle());
146
147   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
148   division->Divide(proj2);
149 }
150