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