]>
Commit | Line | Data |
---|---|---|
04a7657f | 1 | /* $Id$ */ |
2 | ||
3 | #include <AliPWG0Helper.h> | |
4 | ||
5 | #include <TParticle.h> | |
6 | #include <TParticlePDG.h> | |
847489f7 | 7 | #include <TH1.h> |
dd367a14 | 8 | #include <TH2.h> |
25db2d85 | 9 | #include <TH3.h> |
116083b1 | 10 | #include <TList.h> |
9cd015f7 | 11 | #include <TTree.h> |
12 | #include <TBranch.h> | |
13 | #include <TLeaf.h> | |
116083b1 | 14 | |
15 | #include <AliHeader.h> | |
16 | #include <AliStack.h> | |
17 | #include <AliLog.h> | |
04a7657f | 18 | |
19 | #include <AliLog.h> | |
20 | #include <AliESD.h> | |
770a1f1d | 21 | #include <AliESDEvent.h> |
04a7657f | 22 | #include <AliESDVertex.h> |
23 | ||
116083b1 | 24 | #include <AliGenEventHeader.h> |
25 | #include <AliGenPythiaEventHeader.h> | |
26 | #include <AliGenCocktailEventHeader.h> | |
27 | ||
04a7657f | 28 | //____________________________________________________________________ |
29 | ClassImp(AliPWG0Helper) | |
30 | ||
31 | //____________________________________________________________________ | |
116083b1 | 32 | Bool_t AliPWG0Helper::IsEventTriggered(const AliESD* aEsd, Trigger trigger) |
33 | { | |
34 | // see function with ULong64_t argument | |
35 | ||
36 | ULong64_t triggerMask = aEsd->GetTriggerMask(); | |
37 | return IsEventTriggered(triggerMask, trigger); | |
38 | } | |
39 | ||
40 | //____________________________________________________________________ | |
41 | Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger) | |
04a7657f | 42 | { |
43 | // check if the event was triggered | |
44 | // | |
5c495d37 | 45 | // this function needs the branch fTriggerMask |
dd367a14 | 46 | |
47 | ||
48 | // definitions from p-p.cfg | |
49 | ULong64_t spdFO = (1 << 14); | |
50 | ULong64_t v0left = (1 << 11); | |
51 | ULong64_t v0right = (1 << 12); | |
04a7657f | 52 | |
e9c3977b | 53 | switch (trigger) |
54 | { | |
55 | case kMB1: | |
56 | { | |
dd367a14 | 57 | if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))) |
e9c3977b | 58 | return kTRUE; |
59 | break; | |
60 | } | |
61 | case kMB2: | |
62 | { | |
dd367a14 | 63 | if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right))) |
e9c3977b | 64 | return kTRUE; |
65 | break; | |
66 | } | |
67 | } | |
04a7657f | 68 | |
69 | return kFALSE; | |
70 | } | |
71 | ||
72 | //____________________________________________________________________ | |
116083b1 | 73 | Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESD* aEsd) |
74 | { | |
75 | // see function with AliESDVertex argument | |
76 | ||
77 | const AliESDVertex* vtxESD = aEsd->GetVertex(); | |
78 | return IsVertexReconstructed(vtxESD); | |
79 | } | |
80 | ||
81 | //____________________________________________________________________ | |
82 | Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD) | |
04a7657f | 83 | { |
84 | // checks if the vertex is reasonable | |
5c495d37 | 85 | // |
86 | // this function needs the branches fSPDVertex* | |
87 | ||
116083b1 | 88 | if (!vtxESD) |
89 | return kFALSE; | |
04a7657f | 90 | |
91 | // the vertex should be reconstructed | |
92 | if (strcmp(vtxESD->GetName(), "default")==0) | |
93 | return kFALSE; | |
94 | ||
dd367a14 | 95 | // check Ncontributors |
96 | if (vtxESD->GetNContributors() <= 0) | |
97 | return kFALSE; | |
98 | ||
04a7657f | 99 | Double_t vtx_res[3]; |
100 | vtx_res[0] = vtxESD->GetXRes(); | |
101 | vtx_res[1] = vtxESD->GetYRes(); | |
102 | vtx_res[2] = vtxESD->GetZRes(); | |
103 | ||
104 | if (vtx_res[2]==0 || vtx_res[2]>0.1) | |
105 | return kFALSE; | |
106 | ||
107 | return kTRUE; | |
108 | } | |
109 | ||
770a1f1d | 110 | //____________________________________________________________________ |
ebf31fda | 111 | const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug) |
770a1f1d | 112 | { |
113 | // Get the vertex from the ESD and returns it if the vertex is valid | |
114 | // | |
115 | // Second argument decides which vertex is used (this selects | |
116 | // also the quality criteria that are applied) | |
117 | ||
118 | const AliESDVertex* vertex = 0; | |
119 | Float_t requiredZResolution = -1; | |
120 | if (analysisMode == kSPD || analysisMode == kTPCITS) | |
121 | { | |
122 | vertex = aEsd->GetVertex(); | |
123 | requiredZResolution = 0.1; | |
ebf31fda | 124 | if (debug) |
125 | Printf("AliPWG0Helper::GetVertex: Returning SPD vertex"); | |
770a1f1d | 126 | } |
127 | else if (analysisMode == kTPC) | |
128 | { | |
129 | vertex = aEsd->GetPrimaryVertex(); | |
130 | requiredZResolution = 0.6; | |
ebf31fda | 131 | if (debug) |
132 | Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); | |
770a1f1d | 133 | } |
134 | else | |
135 | Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode); | |
136 | ||
ebf31fda | 137 | if (!vertex) { |
138 | if (debug) | |
139 | Printf("AliPWG0Helper::GetVertex: No vertex found in ESD"); | |
770a1f1d | 140 | return 0; |
ebf31fda | 141 | } |
770a1f1d | 142 | |
143 | // check Ncontributors | |
ebf31fda | 144 | if (vertex->GetNContributors() <= 0) { |
145 | if (debug) | |
146 | Printf("AliPWG0Helper::GetVertex: NContributors() <= 0"); | |
770a1f1d | 147 | return 0; |
ebf31fda | 148 | } |
770a1f1d | 149 | |
150 | // check resolution | |
151 | Double_t zRes = vertex->GetZRes(); | |
152 | ||
ebf31fda | 153 | if (zRes == 0 || zRes > requiredZResolution) { |
154 | if (debug) | |
155 | Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution); | |
770a1f1d | 156 | return 0; |
ebf31fda | 157 | } |
158 | ||
159 | if (debug) | |
160 | Printf("AliPWG0Helper::GetVertex: Returning valid vertex"); | |
770a1f1d | 161 | |
162 | return vertex; | |
163 | } | |
164 | ||
04a7657f | 165 | //____________________________________________________________________ |
7584d357 | 166 | Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug) |
04a7657f | 167 | { |
168 | // | |
25db2d85 | 169 | // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack) |
170 | // shall be counted as a primary particle | |
171 | // | |
04a7657f | 172 | // This function or a equivalent should be available in some common place of AliRoot |
173 | // | |
dd367a14 | 174 | // WARNING: Call this function only for particles that are among the particles from the event generator! |
175 | // --> stack->Particle(id) with id < stack->GetNprimary() | |
04a7657f | 176 | |
177 | // if the particle has a daughter primary, we do not want to count it | |
178 | if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries) | |
179 | { | |
7584d357 | 180 | if (adebug) |
25db2d85 | 181 | printf("Dropping particle because it has a daughter among the primaries.\n"); |
04a7657f | 182 | return kFALSE; |
183 | } | |
184 | ||
185 | Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode()); | |
186 | ||
187 | // skip quarks and gluon | |
188 | if (pdgCode <= 10 || pdgCode == 21) | |
189 | { | |
7584d357 | 190 | if (adebug) |
25db2d85 | 191 | printf("Dropping particle because it is a quark or gluon.\n"); |
04a7657f | 192 | return kFALSE; |
193 | } | |
194 | ||
195 | if (strcmp(aParticle->GetName(),"XXX") == 0) | |
196 | { | |
7584d357 | 197 | if (adebug) |
25db2d85 | 198 | printf("WARNING: There is a particle named XXX.\n"); |
04a7657f | 199 | return kFALSE; |
200 | } | |
201 | ||
202 | TParticlePDG* pdgPart = aParticle->GetPDG(); | |
203 | ||
204 | if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0) | |
205 | { | |
7584d357 | 206 | if (adebug) |
25db2d85 | 207 | printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode); |
04a7657f | 208 | return kFALSE; |
209 | } | |
210 | ||
211 | if (pdgPart->Charge() == 0) | |
212 | { | |
7584d357 | 213 | if (adebug) |
25db2d85 | 214 | printf("Dropping particle because it is not charged.\n"); |
04a7657f | 215 | return kFALSE; |
216 | } | |
217 | ||
218 | return kTRUE; | |
219 | } | |
847489f7 | 220 | |
221 | //____________________________________________________________________ | |
29771dc8 | 222 | void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save) |
847489f7 | 223 | { |
224 | // create projections of 3d hists to all 2d combinations | |
225 | // the histograms are not returned, just use them from memory or use this to create them in a file | |
226 | ||
227 | TH1* proj = hist->Project3D("yx"); | |
228 | proj->SetXTitle(hist->GetXaxis()->GetTitle()); | |
229 | proj->SetYTitle(hist->GetYaxis()->GetTitle()); | |
29771dc8 | 230 | if (save) |
231 | proj->Write(); | |
847489f7 | 232 | |
233 | proj = hist->Project3D("zx"); | |
234 | proj->SetXTitle(hist->GetXaxis()->GetTitle()); | |
235 | proj->SetYTitle(hist->GetZaxis()->GetTitle()); | |
29771dc8 | 236 | if (save) |
237 | proj->Write(); | |
847489f7 | 238 | |
239 | proj = hist->Project3D("zy"); | |
240 | proj->SetXTitle(hist->GetYaxis()->GetTitle()); | |
241 | proj->SetYTitle(hist->GetZaxis()->GetTitle()); | |
29771dc8 | 242 | if (save) |
243 | proj->Write(); | |
847489f7 | 244 | } |
1afae8ff | 245 | |
246 | //____________________________________________________________________ | |
29771dc8 | 247 | void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save) |
1afae8ff | 248 | { |
249 | // create projections of the 3d hists divides them | |
250 | // axis decides to which plane, if axis is 0 to all planes | |
251 | // the histograms are not returned, just use them from memory or use this to create them in a file | |
252 | ||
253 | if (axis == 0) | |
254 | { | |
29771dc8 | 255 | CreateDividedProjections(hist, hist2, "yx", putErrors, save); |
256 | CreateDividedProjections(hist, hist2, "zx", putErrors, save); | |
257 | CreateDividedProjections(hist, hist2, "zy", putErrors, save); | |
1afae8ff | 258 | |
259 | return; | |
260 | } | |
261 | ||
262 | TH1* proj = hist->Project3D(axis); | |
0ab29cfa | 263 | |
264 | if (strlen(axis) == 2) | |
265 | { | |
266 | proj->SetYTitle(GetAxisTitle(hist, axis[0])); | |
267 | proj->SetXTitle(GetAxisTitle(hist, axis[1])); | |
268 | } | |
269 | else if (strlen(axis) == 1) | |
270 | proj->SetXTitle(GetAxisTitle(hist, axis[0])); | |
1afae8ff | 271 | |
272 | TH1* proj2 = hist2->Project3D(axis); | |
0ab29cfa | 273 | if (strlen(axis) == 2) |
274 | { | |
275 | proj2->SetYTitle(GetAxisTitle(hist2, axis[0])); | |
276 | proj2->SetXTitle(GetAxisTitle(hist2, axis[1])); | |
277 | } | |
278 | else if (strlen(axis) == 1) | |
279 | proj2->SetXTitle(GetAxisTitle(hist2, axis[0])); | |
1afae8ff | 280 | |
281 | TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName()))); | |
29771dc8 | 282 | //printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX())); |
283 | //printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY())); | |
9e952c39 | 284 | division->Divide(proj, proj2, 1, 1, "B"); |
29771dc8 | 285 | division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle())); |
0ab29cfa | 286 | |
287 | if (putErrors) | |
288 | { | |
289 | division->Sumw2(); | |
290 | if (division->GetDimension() == 1) | |
291 | { | |
292 | Int_t nBins = division->GetNbinsX(); | |
29771dc8 | 293 | for (Int_t i = 1; i <= nBins; ++i) |
0ab29cfa | 294 | if (proj2->GetBinContent(i) != 0) |
295 | division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i)); | |
296 | } | |
297 | else if (division->GetDimension() == 2) | |
298 | { | |
299 | Int_t nBinsX = division->GetNbinsX(); | |
300 | Int_t nBinsY = division->GetNbinsY(); | |
29771dc8 | 301 | for (Int_t i = 1; i <= nBinsX; ++i) |
302 | for (Int_t j = 1; j <= nBinsY; ++j) | |
0ab29cfa | 303 | if (proj2->GetBinContent(i, j) != 0) |
304 | division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j)); | |
305 | } | |
306 | } | |
29771dc8 | 307 | |
308 | if (save) | |
309 | { | |
310 | proj->Write(); | |
311 | proj2->Write(); | |
312 | division->Write(); | |
313 | } | |
1afae8ff | 314 | } |
315 | ||
92d2d8ad | 316 | //____________________________________________________________________ |
25db2d85 | 317 | const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis) |
92d2d8ad | 318 | { |
319 | // returns the title of the axis given in axis (x, y, z) | |
320 | ||
321 | if (axis == 'x') | |
322 | return hist->GetXaxis()->GetTitle(); | |
323 | else if (axis == 'y') | |
324 | return hist->GetYaxis()->GetTitle(); | |
325 | else if (axis == 'z') | |
326 | return hist->GetZaxis()->GetTitle(); | |
327 | ||
328 | return 0; | |
329 | } | |
116083b1 | 330 | |
331 | //____________________________________________________________________ | |
332 | Int_t AliPWG0Helper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) { | |
333 | // | |
334 | // get the process type of the event. | |
335 | // | |
336 | ||
337 | // can only read pythia headers, either directly or from cocktalil header | |
338 | AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader()); | |
339 | ||
340 | if (!pythiaGenHeader) { | |
341 | ||
342 | AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader()); | |
343 | if (!genCocktailHeader) { | |
344 | printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n"); | |
345 | return -1; | |
346 | } | |
347 | ||
348 | TList* headerList = genCocktailHeader->GetHeaders(); | |
349 | if (!headerList) { | |
350 | return -1; | |
351 | } | |
352 | ||
353 | for (Int_t i=0; i<headerList->GetEntries(); i++) { | |
354 | pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i)); | |
355 | if (pythiaGenHeader) | |
356 | break; | |
357 | } | |
358 | ||
359 | if (!pythiaGenHeader) { | |
360 | printf("AliPWG0Helper::GetProcessType : Could not find Pythia header. \n"); | |
361 | return -1; | |
362 | } | |
363 | } | |
364 | ||
365 | if (adebug) { | |
366 | printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType()); | |
367 | } | |
368 | ||
369 | return pythiaGenHeader->ProcessType(); | |
370 | } | |
371 | ||
372 | //____________________________________________________________________ | |
373 | TParticle* AliPWG0Helper::FindPrimaryMother(AliStack* stack, Int_t label) | |
374 | { | |
375 | // | |
376 | // Finds the first mother among the primary particles of the particle identified by <label>, | |
377 | // i.e. the primary that "caused" this particle | |
378 | // | |
379 | ||
380 | Int_t motherLabel = FindPrimaryMotherLabel(stack, label); | |
381 | if (motherLabel < 0) | |
382 | return 0; | |
383 | ||
384 | return stack->Particle(motherLabel); | |
385 | } | |
386 | ||
387 | //____________________________________________________________________ | |
388 | Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label) | |
389 | { | |
390 | // | |
391 | // Finds the first mother among the primary particles of the particle identified by <label>, | |
392 | // i.e. the primary that "caused" this particle | |
393 | // | |
394 | // returns its label | |
395 | // | |
396 | ||
397 | Int_t nPrim = stack->GetNprimary(); | |
398 | ||
399 | while (label >= nPrim) | |
400 | { | |
401 | //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0)); | |
402 | ||
403 | TParticle* particle = stack->Particle(label); | |
404 | if (!particle) | |
405 | { | |
406 | AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label)); | |
407 | return -1; | |
408 | } | |
409 | ||
410 | // find mother | |
411 | if (particle->GetMother(0) < 0) | |
412 | { | |
413 | AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label)); | |
414 | return -1; | |
415 | } | |
416 | ||
417 | label = particle->GetMother(0); | |
418 | } | |
419 | ||
420 | return label; | |
421 | } | |
9cd015f7 | 422 | |
423 | void AliPWG0Helper::SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug) | |
424 | { | |
425 | // Function to switch on/off all data members of a top level branch | |
426 | // this is needed for branches without a trailing dot ".", for those | |
427 | // the root functionality with regular expressions does not work. | |
428 | // Usage e.g. | |
429 | // chain->SetBranchStatus("*", 0); | |
430 | // SetBranchStatusRecursive(chain,"SPDVertex",1); | |
431 | // You need to give the full name of the top level branch zou want to access | |
432 | //========================================================================== | |
433 | // Author Christian.Klein-Boesing@cern.ch | |
434 | ||
435 | if (!tree) | |
436 | return; | |
437 | ||
438 | TBranch *br = tree->GetBranch(bname); | |
439 | if(!br) { | |
440 | Printf("AliPWG0Helper::SetBranchStatusRecursive: Branch %s not found", bname); | |
441 | } | |
442 | ||
443 | TObjArray *leaves = tree->GetListOfLeaves(); | |
444 | Int_t nleaves = leaves->GetEntries(); | |
445 | TLeaf *leaf = 0; | |
446 | TBranch *branch = 0; | |
447 | TBranch *mother = 0; | |
448 | for (Int_t i=0;i<nleaves;i++) { | |
449 | // the matched entry e.g. SPDVertex is its own Mother | |
450 | leaf = (TLeaf*)leaves->UncheckedAt(i); | |
451 | branch = (TBranch*)leaf->GetBranch(); | |
452 | mother = branch->GetMother(); | |
453 | if (mother==br){ | |
454 | if (debug) | |
455 | Printf(">>>> AliPWG0Helper::SetBranchStatusRecursive: Setting status %s %s to %d", mother->GetName(), leaf->GetName(), status); | |
456 | if (status) branch->ResetBit(kDoNotProcess); | |
457 | else branch->SetBit(kDoNotProcess); | |
458 | } | |
459 | } | |
460 | } | |
dd367a14 | 461 | |
462 | //____________________________________________________________________ | |
463 | void AliPWG0Helper::NormalizeToBinWidth(TH1* hist) | |
464 | { | |
465 | // | |
466 | // normalizes a 1-d histogram to its bin width | |
467 | // | |
468 | ||
469 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
470 | { | |
471 | hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); | |
472 | hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); | |
473 | } | |
474 | } | |
475 | ||
476 | //____________________________________________________________________ | |
477 | void AliPWG0Helper::NormalizeToBinWidth(TH2* hist) | |
478 | { | |
479 | // | |
480 | // normalizes a 2-d histogram to its bin width (x width * y width) | |
481 | // | |
482 | ||
483 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
484 | for (Int_t j=1; j<=hist->GetNbinsY(); ++j) | |
485 | { | |
486 | Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); | |
487 | hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); | |
488 | hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); | |
489 | } | |
490 | } |