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