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