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