possibility to read process type from phojet events (Christian)
[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 <TH2.h>
9 #include <TH3.h>
10 #include <TList.h>
11 #include <TTree.h>
12 #include <TBranch.h>
13 #include <TLeaf.h>
14
15 #include <AliHeader.h>
16 #include <AliStack.h>
17 #include <AliLog.h>
18
19 #include <AliLog.h>
20 #include <AliESD.h>
21 #include <AliESDEvent.h>
22 #include <AliESDVertex.h>
23
24 #include <AliGenEventHeader.h>
25 #include <AliGenPythiaEventHeader.h>
26 #include <AliGenCocktailEventHeader.h>
27 #include <AliGenDPMjetEventHeader.h>
28
29 //____________________________________________________________________
30 ClassImp(AliPWG0Helper)
31
32 //____________________________________________________________________
33 Bool_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 //____________________________________________________________________
42 Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
43 {
44   // check if the event was triggered
45   //
46   // this function needs the branch fTriggerMask
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);
53
54   switch (trigger)
55   {
56     case kMB1:
57     {
58       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
59         return kTRUE;
60       break;
61     }
62     case kMB2:
63     {
64       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
65         return kTRUE;
66       break;
67     }
68   }
69
70   return kFALSE;
71 }
72
73 //____________________________________________________________________
74 Bool_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 //____________________________________________________________________
83 Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
84 {
85   // checks if the vertex is reasonable
86   //
87   // this function needs the branches fSPDVertex*
88
89   if (!vtxESD)
90     return kFALSE;
91
92   // the vertex should be reconstructed
93   if (strcmp(vtxESD->GetName(), "default")==0)
94     return kFALSE;
95
96   // check Ncontributors
97   if (vtxESD->GetNContributors() <= 0)
98     return kFALSE;
99
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 //____________________________________________________________________
112 const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
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   {
123     vertex = aEsd->GetPrimaryVertexSPD();
124     requiredZResolution = 0.1;
125     if (debug)
126       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
127   }
128   else if (analysisMode == kTPC) 
129   {
130     vertex = aEsd->GetPrimaryVertexTPC();
131     requiredZResolution = 0.6;
132     if (debug)
133       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
134   }
135   else
136     Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
137
138   if (!vertex) {
139     if (debug)
140       Printf("AliPWG0Helper::GetVertex: No vertex found in ESD");
141     return 0;
142   }
143
144   // check Ncontributors
145   if (vertex->GetNContributors() <= 0) {
146     if (debug)
147       Printf("AliPWG0Helper::GetVertex: NContributors() <= 0");
148     return 0;
149   }
150
151   // check resolution
152   Double_t zRes = vertex->GetZRes();
153
154   if (zRes == 0 || zRes > requiredZResolution) {
155     if (debug)
156       Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
157     return 0;
158   }
159
160   if (debug)
161     Printf("AliPWG0Helper::GetVertex: Returning valid vertex");
162
163   return vertex;
164 }
165
166 //____________________________________________________________________
167 Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
168 {
169   //
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   //
173   // This function or a equivalent should be available in some common place of AliRoot
174   //
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()
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   {
181     if (adebug)
182       printf("Dropping particle because it has a daughter among the primaries.\n");
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   {
191     if (adebug)
192       printf("Dropping particle because it is a quark or gluon.\n");
193     return kFALSE;
194   }
195
196   if (strcmp(aParticle->GetName(),"XXX") == 0)
197   {
198     if (adebug)
199       printf("WARNING: There is a particle named XXX.\n");
200     return kFALSE;
201   }
202
203   TParticlePDG* pdgPart = aParticle->GetPDG();
204
205   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
206   {
207     if (adebug)
208       printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
209     return kFALSE;
210   }
211
212   if (pdgPart->Charge() == 0)
213   {
214     if (adebug)
215       printf("Dropping particle because it is not charged.\n");
216     return kFALSE;
217   }
218
219   return kTRUE;
220 }
221
222 //____________________________________________________________________
223 void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
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());
231   if (save)
232     proj->Write();
233
234   proj = hist->Project3D("zx");
235   proj->SetXTitle(hist->GetXaxis()->GetTitle());
236   proj->SetYTitle(hist->GetZaxis()->GetTitle());
237   if (save)
238     proj->Write();
239
240   proj = hist->Project3D("zy");
241   proj->SetXTitle(hist->GetYaxis()->GetTitle());
242   proj->SetYTitle(hist->GetZaxis()->GetTitle());
243   if (save)
244     proj->Write();
245 }
246
247 //____________________________________________________________________
248 void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
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   {
256     CreateDividedProjections(hist, hist2, "yx", putErrors, save);
257     CreateDividedProjections(hist, hist2, "zx", putErrors, save);
258     CreateDividedProjections(hist, hist2, "zy", putErrors, save);
259
260     return;
261   }
262
263   TH1* proj = hist->Project3D(axis);
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]));
272
273   TH1* proj2 = hist2->Project3D(axis);
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]));
281
282   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
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()));
285   division->Divide(proj, proj2, 1, 1, "B");
286   division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
287
288   if (putErrors)
289   {
290     division->Sumw2();
291     if (division->GetDimension() == 1)
292     {
293       Int_t nBins = division->GetNbinsX();
294       for (Int_t i = 1; i <= nBins; ++i)
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();
302       for (Int_t i = 1; i <= nBinsX; ++i)
303         for (Int_t j = 1; j <= nBinsY; ++j)
304           if (proj2->GetBinContent(i, j) != 0)
305             division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
306     }
307   }
308
309   if (save)
310   {
311     proj->Write();
312     proj2->Write();
313     division->Write();
314   }
315 }
316
317 //____________________________________________________________________
318 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
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 }
331
332
333 Int_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
366 Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
367   //
368   // get the process type of the event.
369   //
370
371   // can only read pythia headers, either directly or from cocktalil header
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
401 Int_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
409   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
410   if (pythiaGenHeader) {
411     return GetPythiaEventProcessType(pythiaGenHeader,adebug);
412   }
413
414   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
415   if (dpmJetGenHeader) {
416     return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
417   }
418   
419
420   // check for cocktail
421
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   }
427
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);
438     }
439
440     dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
441     if (dpmJetGenHeader) {
442       return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
443     }
444   }
445   return kInvalidProcess;
446 }
447
448
449
450 //____________________________________________________________________
451 TParticle* 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 //____________________________________________________________________
466 Int_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 }
500
501 void 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 }
539
540 //____________________________________________________________________
541 void 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 //____________________________________________________________________
555 void 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 }
569
570 //____________________________________________________________________
571 void 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   {
581     case kInvalid: str += "invalid setting"; break;
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