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