adding flag for non-field data
[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 #include <AliVertexerTracks.h>
24 #include <AliMultiplicity.h>
25
26 #include <AliGenEventHeader.h>
27 #include <AliGenPythiaEventHeader.h>
28 #include <AliGenCocktailEventHeader.h>
29 #include <AliGenDPMjetEventHeader.h>
30 #include <AliESDVZERO.h>
31
32 //____________________________________________________________________
33 ClassImp(AliPWG0Helper)
34
35 Int_t AliPWG0Helper::fgLastProcessType = -1;
36
37 //____________________________________________________________________
38 Bool_t AliPWG0Helper::IsEventTriggered(const AliESDEvent* aEsd, Trigger trigger)
39 {
40   // checks if an event has been triggered
41   // this function implements the "offline" methods that use the ESD, other trigger requests are passed to the function prototype with ULong_t
42
43   Int_t firedChips = 0;
44   Bool_t v0A = kFALSE;
45   Bool_t v0C = kFALSE;
46
47   // offline triggers have to be dealt with here, because we need the esd pointer
48   if (trigger == kOfflineFASTOR || trigger == kOfflineMB1 || trigger == kOfflineMB2 || trigger == kOfflineMB3)
49   {
50       const AliMultiplicity* mult = aEsd->GetMultiplicity();
51       if (!mult)
52       {
53         Printf("AliPWG0Helper::IsEventTriggered: ERROR: AliMultiplicity not available");
54         return kFALSE;
55       }
56       firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
57   }
58   if (trigger == kOfflineMB1 || trigger == kOfflineMB2 || trigger == kOfflineMB3)
59   {
60     AliESDVZERO* v0Data = aEsd->GetVZEROData();
61     if (!v0Data)
62     {
63       Printf("AliPWG0Helper::IsEventTriggered: ERROR: AliESDVZERO not available");
64       return kFALSE;
65     }
66     for (Int_t i=0; i<32; i++)
67     {
68       if (v0Data->BBTriggerV0A(i))
69         v0A = kTRUE;
70       if (v0Data->BBTriggerV0C(i))
71         v0C = kTRUE;
72     }
73   }
74       
75   switch (trigger)
76   {
77     case kOfflineFASTOR:
78     {
79       if (firedChips > 0)
80         return kTRUE;
81       break;
82     }
83     case kOfflineMB1:
84     {
85       if ((firedChips > 0) || v0A || v0C)
86         return kTRUE;
87       break;
88     }
89     case kOfflineMB2:
90     {
91       if ((firedChips > 0) && (v0A || v0C))
92         return kTRUE;
93       break;
94     }
95     case kOfflineMB3:
96     {
97       if ((firedChips > 0) && v0A && v0C)
98         return kTRUE;
99       break;
100     }
101     default:
102     {
103       return IsEventTriggered(aEsd->GetTriggerMask(), trigger);
104       break;
105     }
106   }
107   
108   return kFALSE;
109 }
110
111 //____________________________________________________________________
112 Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
113 {
114   // check if the event was triggered
115   //
116   // this function needs the branch fTriggerMask
117   
118   // definitions from p-p.cfg
119   ULong64_t spdFO = (1 << 14);
120   ULong64_t v0left = (1 << 10);
121   ULong64_t v0right = (1 << 11);
122
123   switch (trigger)
124   {
125     case kAcceptAll:
126     {
127       return kTRUE;
128       break;
129     }
130     case kMB1:
131     {
132       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
133         return kTRUE;
134       break;
135     }
136     case kMB2:
137     {
138       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
139         return kTRUE;
140       break;
141     }
142     case kMB3:
143     {
144       if (triggerMask & spdFO && (triggerMask & v0left) && (triggerMask & v0right))
145         return kTRUE;
146       break;
147     }
148     case kSPDFASTOR:
149     {
150       if (triggerMask & spdFO)
151         return kTRUE;
152       break;
153     }
154     default:
155       Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
156       break;
157   }
158
159   return kFALSE;
160 }
161
162 //____________________________________________________________________
163 Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
164 {
165     // Checks if a vertex meets the needed quality criteria
166
167   Float_t requiredZResolution = -1;
168   if (analysisMode & kSPD || analysisMode & kTPCITS)
169   {
170     requiredZResolution = 0.1;
171   }
172   else if (analysisMode & kTPC)
173     requiredZResolution = 10.;
174
175   // check resolution
176   Double_t zRes = vertex->GetZRes();
177
178   if (zRes > requiredZResolution) {
179     if (debug)
180       Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
181     return kFALSE;
182   }
183
184   return kTRUE;
185 }
186
187 //____________________________________________________________________
188 const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC)
189 {
190   // Get the vertex from the ESD and returns it if the vertex is valid
191   //
192   // Second argument decides which vertex is used (this selects
193   // also the quality criteria that are applied)
194
195   const AliESDVertex* vertex = 0;
196   if (analysisMode & kSPD || analysisMode & kTPCITS)
197   {
198     vertex = aEsd->GetPrimaryVertexSPD();
199     if (debug)
200       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
201   }
202   else if (analysisMode & kTPC)
203   {
204     if(bRedoTPC){
205       if (debug)
206         Printf("AliPWG0Helper::GetVertex: Redoing vertex");
207       Double_t kBz = aEsd->GetMagneticField();
208       AliVertexerTracks vertexer(kBz);
209       vertexer.SetTPCMode();
210       AliESDVertex *vTPC = vertexer.FindPrimaryVertex(aEsd);
211       aEsd->SetPrimaryVertexTPC(vTPC);
212       for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
213         AliESDtrack *t = aEsd->GetTrack(i);
214         t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
215       }
216       delete vTPC;
217     }
218
219     vertex = aEsd->GetPrimaryVertexTPC();
220     if (debug)
221       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
222   }
223   else
224     Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
225
226   if (!vertex) {
227     if (debug)
228       Printf("AliPWG0Helper::GetVertex: No vertex found in ESD");
229     return 0;
230   }
231
232   // check Ncontributors
233   if (vertex->GetNContributors() <= 0) {
234     if (debug){
235       Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
236       Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices());
237       vertex->Print();
238     }
239     return 0;
240   }
241
242   // check resolution
243   Double_t zRes = vertex->GetZRes();
244   if (zRes == 0) {
245     Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0.");
246     return 0;
247   }
248
249   if (debug)
250   {
251     Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
252     vertex->Print();
253   }
254
255   return vertex;
256 }
257
258 //____________________________________________________________________
259 Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
260 {
261   //
262   // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
263   // shall be counted as a primary particle
264   //
265   // This function or a equivalent should be available in some common place of AliRoot
266   //
267   // WARNING: Call this function only for particles that are among the particles from the event generator!
268   // --> stack->Particle(id) with id < stack->GetNprimary()
269
270   // if the particle has a daughter primary, we do not want to count it
271   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
272   {
273     if (adebug)
274       printf("Dropping particle because it has a daughter among the primaries.\n");
275     return kFALSE;
276   }
277
278   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
279   
280
281   // skip quarks and gluon
282   if (pdgCode <= 10 || pdgCode == 21)
283   {
284     if (adebug)
285       printf("Dropping particle because it is a quark or gluon.\n");
286     return kFALSE;
287   }
288
289   Int_t status = aParticle->GetStatusCode();
290   // skip non final state particles..
291   if(status!=1){
292     if (adebug)
293       printf("Dropping particle because it is not a final state particle.\n");
294     return kFALSE;
295   }
296
297   if (strcmp(aParticle->GetName(),"XXX") == 0)
298   {
299     Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
300     return kFALSE;
301   }
302
303   TParticlePDG* pdgPart = aParticle->GetPDG();
304
305   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
306   {
307     Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
308     return kFALSE;
309   }
310
311   if (pdgPart->Charge() == 0)
312   {
313     if (adebug)
314       printf("Dropping particle because it is not charged.\n");
315     return kFALSE;
316   }
317
318   return kTRUE;
319 }
320
321 //____________________________________________________________________
322 void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
323 {
324   // create projections of 3d hists to all 2d combinations
325   // the histograms are not returned, just use them from memory or use this to create them in a file
326
327   TH1* proj = hist->Project3D("yx");
328   proj->SetXTitle(hist->GetXaxis()->GetTitle());
329   proj->SetYTitle(hist->GetYaxis()->GetTitle());
330   if (save)
331     proj->Write();
332
333   proj = hist->Project3D("zx");
334   proj->SetXTitle(hist->GetXaxis()->GetTitle());
335   proj->SetYTitle(hist->GetZaxis()->GetTitle());
336   if (save)
337     proj->Write();
338
339   proj = hist->Project3D("zy");
340   proj->SetXTitle(hist->GetYaxis()->GetTitle());
341   proj->SetYTitle(hist->GetZaxis()->GetTitle());
342   if (save)
343     proj->Write();
344 }
345
346 //____________________________________________________________________
347 void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
348 {
349   // create projections of the 3d hists divides them
350   // axis decides to which plane, if axis is 0 to all planes
351   // the histograms are not returned, just use them from memory or use this to create them in a file
352
353   if (axis == 0)
354   {
355     CreateDividedProjections(hist, hist2, "yx", putErrors, save);
356     CreateDividedProjections(hist, hist2, "zx", putErrors, save);
357     CreateDividedProjections(hist, hist2, "zy", putErrors, save);
358
359     return;
360   }
361
362   TH1* proj = hist->Project3D(axis);
363
364   if (strlen(axis) == 2)
365   {
366     proj->SetYTitle(GetAxisTitle(hist, axis[0]));
367     proj->SetXTitle(GetAxisTitle(hist, axis[1]));
368   }
369   else if (strlen(axis) == 1)
370     proj->SetXTitle(GetAxisTitle(hist, axis[0]));
371
372   TH1* proj2 = hist2->Project3D(axis);
373   if (strlen(axis) == 2)
374   {
375     proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
376     proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
377   }
378   else if (strlen(axis) == 1)
379     proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
380
381   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
382   //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()));
383   //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()));
384   division->Divide(proj, proj2, 1, 1, "B");
385   division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
386
387   if (putErrors)
388   {
389     division->Sumw2();
390     if (division->GetDimension() == 1)
391     {
392       Int_t nBins = division->GetNbinsX();
393       for (Int_t i = 1; i <= nBins; ++i)
394         if (proj2->GetBinContent(i) != 0)
395           division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
396     }
397     else if (division->GetDimension() == 2)
398     {
399       Int_t nBinsX = division->GetNbinsX();
400       Int_t nBinsY = division->GetNbinsY();
401       for (Int_t i = 1; i <= nBinsX; ++i)
402         for (Int_t j = 1; j <= nBinsY; ++j)
403           if (proj2->GetBinContent(i, j) != 0)
404             division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
405     }
406   }
407
408   if (save)
409   {
410     proj->Write();
411     proj2->Write();
412     division->Write();
413   }
414 }
415
416 //____________________________________________________________________
417 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
418 {
419   // returns the title of the axis given in axis (x, y, z)
420
421   if (axis == 'x')
422     return hist->GetXaxis()->GetTitle();
423   else if (axis == 'y')
424     return hist->GetYaxis()->GetTitle();
425   else if (axis == 'z')
426     return hist->GetZaxis()->GetTitle();
427
428   return 0;
429 }
430
431
432 AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
433
434   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
435
436   if (!pythiaGenHeader) {
437     printf("AliPWG0Helper::GetProcessType : Unknown gen Header type). \n");
438     return kInvalidProcess;
439   }
440
441
442   Int_t pythiaType = pythiaGenHeader->ProcessType();
443   fgLastProcessType = pythiaType;
444   MCProcessType globalType = kInvalidProcess;  
445
446
447   if (adebug) {
448     printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaType);
449   }
450
451
452   if(pythiaType==92||pythiaType==93){
453     globalType = kSD;
454   }
455   else if(pythiaType==94){
456     globalType = kDD;
457   }
458   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
459   else {
460     globalType = kND;
461   }
462   return globalType;
463 }
464
465
466 AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
467   //
468   // get the process type of the event.
469   //
470
471   // can only read pythia headers, either directly or from cocktalil header
472   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
473
474   if (!dpmJetGenHeader) {
475     printf("AliPWG0Helper::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
476     return kInvalidProcess;
477   }
478
479   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
480   fgLastProcessType = dpmJetType;
481   MCProcessType globalType = kInvalidProcess;  
482
483
484   if (adebug) {
485     printf("AliPWG0Helper::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
486   }
487
488
489   if(dpmJetType == 1){ // this is explicitly inelastic
490     globalType = kND;
491   }  
492   else if(dpmJetType==5||dpmJetType==6){
493     globalType = kSD;
494   }
495   else if(dpmJetType==7||dpmJetType==4){// DD and double pomeron
496     globalType = kDD;
497   }
498   return globalType;
499 }
500
501
502 AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
503   //
504   // get the process type of the event.
505   //
506
507
508   // Check for simple headers first
509
510   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
511   if (pythiaGenHeader) {
512     return GetPythiaEventProcessType(pythiaGenHeader,adebug);
513   }
514
515   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
516   if (dpmJetGenHeader) {
517     return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
518   }
519   
520
521   // check for cocktail
522
523   AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
524   if (!genCocktailHeader) {
525     printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
526     return kInvalidProcess;
527   }
528
529   TList* headerList = genCocktailHeader->GetHeaders();
530   if (!headerList) {
531     return kInvalidProcess;
532   }
533
534   for (Int_t i=0; i<headerList->GetEntries(); i++) {
535
536     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
537     if (pythiaGenHeader) {
538       return GetPythiaEventProcessType(pythiaGenHeader,adebug);
539     }
540
541     dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
542     if (dpmJetGenHeader) {
543       return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
544     }
545   }
546   return kInvalidProcess;
547 }
548
549
550
551 //____________________________________________________________________
552 TParticle* AliPWG0Helper::FindPrimaryMother(AliStack* stack, Int_t label)
553 {
554   //
555   // Finds the first mother among the primary particles of the particle identified by <label>,
556   // i.e. the primary that "caused" this particle
557   //
558
559   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
560   if (motherLabel < 0)
561     return 0;
562
563   return stack->Particle(motherLabel);
564 }
565
566 //____________________________________________________________________
567 Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
568 {
569   //
570   // Finds the first mother among the primary particles of the particle identified by <label>,
571   // i.e. the primary that "caused" this particle
572   //
573   // returns its label
574   //
575
576   Int_t nPrim  = stack->GetNprimary();
577
578   while (label >= nPrim)
579   {
580     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
581
582     TParticle* particle = stack->Particle(label);
583     if (!particle)
584     {
585       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
586       return -1;
587     }
588  
589     // find mother
590     if (particle->GetMother(0) < 0)
591     {
592       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
593       return -1;
594     }
595
596     label = particle->GetMother(0);
597   }
598
599   return label;
600 }
601
602 //____________________________________________________________________
603 void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
604 {
605   //
606   // normalizes a 1-d histogram to its bin width
607   //
608
609   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
610   {
611     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
612     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
613   }
614 }
615
616 //____________________________________________________________________
617 void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
618 {
619   //
620   // normalizes a 2-d histogram to its bin width (x width * y width)
621   //
622
623   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
624     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
625     {
626       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
627       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
628       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
629     }
630 }
631
632 //____________________________________________________________________
633 void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
634 {
635   //
636   // Prints the given configuration
637   //
638
639   TString str(">>>> Running with >");
640
641   if (analysisMode & kSPD)
642     str += "SPD-only";
643     
644   if (analysisMode & kTPC)
645      str += "TPC-only";
646     
647   if (analysisMode & kTPCITS)
648      str += "Global tracking";
649
650   if (analysisMode & kFieldOn)
651   {
652      str += " (with field)";
653   }
654   else
655      str += " (WITHOUT field)";
656   
657   str += "< and trigger >";
658
659   switch (trigger)
660   {
661     case kAcceptAll : str += "ACCEPT ALL (bypass!)"; break;
662     case kMB1 : str += "MB1"; break;
663     case kMB2 : str += "MB2"; break;
664     case kMB3 : str += "MB3"; break;
665     case kSPDFASTOR : str += "SPD FASTOR"; break;
666     case kOfflineMB1 : str += "Offline MB1"; break;
667     case kOfflineMB2 : str += "Offline MB2"; break;
668     case kOfflineMB3 : str += "Offline MB3"; break;
669     case kOfflineFASTOR : str += "Offline SPD FASTOR"; break;
670   }
671
672   str += "< <<<<";
673
674   Printf("%s", str.Data());
675 }
676