Coverity fixes
[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 <AliESDRun.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::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
39 {
40     // Checks if a vertex meets the needed quality criteria
41
42   Float_t requiredZResolution = -1;
43   if (analysisMode & kSPD || analysisMode & kTPCITS || analysisMode & kTPCSPD)
44   {
45     // disable cut on resolution
46     requiredZResolution = 1000;
47   }
48   else if (analysisMode & kTPC)
49     requiredZResolution = 10.;
50
51   // check resolution
52   Double_t zRes = vertex->GetZRes();
53
54   if (zRes > requiredZResolution) {
55     if (debug)
56       Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
57     return kFALSE;
58   }
59   
60   if (vertex->IsFromVertexerZ())
61   {
62     if (vertex->GetDispersion() > 0.02) 
63     {
64       if (debug)
65         Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02);
66       return kFALSE;
67     }
68   }
69
70   return kTRUE;
71 }
72
73 //____________________________________________________________________
74 const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
75 {
76   // Get the vertex from the ESD and returns it if the vertex is valid
77   //
78   // Second argument decides which vertex is used (this selects
79   // also the quality criteria that are applied)
80
81   const AliESDVertex* vertex = 0;
82   if (analysisMode & kSPD)
83   {
84     vertex = aEsd->GetPrimaryVertexSPD();
85     if (debug)
86       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
87   }
88   else if (analysisMode & kTPCITS || analysisMode & kTPCSPD)
89   {
90     vertex = aEsd->GetPrimaryVertexTracks();
91     if (debug)
92       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
93     if (!vertex || vertex->GetNContributors() <= 0)
94     {
95       if (debug)
96         Printf("AliPWG0Helper::GetVertex: Vertex from tracks has no contributors. Falling back to SPD vertex.");
97       vertex = aEsd->GetPrimaryVertexSPD();
98     }
99   }
100   else if (analysisMode & kTPC)
101   {
102     vertex = aEsd->GetPrimaryVertexTPC();
103     if (debug)
104       Printf("AliPWG0Helper::GetVertex: Returning vertex from TPC-only tracks");
105   }
106   else
107     Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
108
109   if (!vertex) {
110     if (debug)
111       Printf("AliPWG0Helper::GetVertex: No vertex found in ESD");
112     return 0;
113   }
114
115   // check Ncontributors
116   if (vertex->GetNContributors() <= 0) {
117     if (debug){
118       Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
119       Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices());
120       vertex->Print();
121     }
122     return 0;
123   }
124
125   // check resolution
126   Double_t zRes = vertex->GetZRes();
127   if (zRes == 0) {
128     Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0.");
129     return 0;
130   }
131
132   if (debug)
133   {
134     Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
135     vertex->Print();
136   }
137
138   return vertex;
139 }
140
141 //____________________________________________________________________
142 Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
143 {
144   //
145   // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
146   // shall be counted as a primary particle
147   //
148   // This function or a equivalent should be available in some common place of AliRoot
149   //
150   // WARNING: Call this function only for particles that are among the particles from the event generator!
151   // --> stack->Particle(id) with id < stack->GetNprimary()
152
153   // if the particle has a daughter primary, we do not want to count it
154   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
155   {
156     if (adebug)
157       printf("Dropping particle because it has a daughter among the primaries.\n");
158     return kFALSE;
159   }
160
161   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
162   
163
164   // skip quarks and gluon
165   if (pdgCode <= 10 || pdgCode == 21)
166   {
167     if (adebug)
168       printf("Dropping particle because it is a quark or gluon.\n");
169     return kFALSE;
170   }
171
172   Int_t status = aParticle->GetStatusCode();
173   // skip non final state particles..
174   if(status!=1){
175     if (adebug)
176       printf("Dropping particle because it is not a final state particle.\n");
177     return kFALSE;
178   }
179
180   if (strcmp(aParticle->GetName(),"XXX") == 0)
181   {
182     Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
183     return kFALSE;
184   }
185
186   TParticlePDG* pdgPart = aParticle->GetPDG();
187
188   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
189   {
190     Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
191     return kFALSE;
192   }
193
194   if (pdgPart->Charge() == 0)
195   {
196     if (adebug)
197       printf("Dropping particle because it is not charged.\n");
198     return kFALSE;
199   }
200
201   return kTRUE;
202 }
203
204 //____________________________________________________________________
205 void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
206 {
207   // create projections of 3d hists to all 2d combinations
208   // the histograms are not returned, just use them from memory or use this to create them in a file
209
210   TH1* proj = hist->Project3D("yx");
211   proj->SetXTitle(hist->GetXaxis()->GetTitle());
212   proj->SetYTitle(hist->GetYaxis()->GetTitle());
213   if (save)
214     proj->Write();
215
216   proj = hist->Project3D("zx");
217   proj->SetXTitle(hist->GetXaxis()->GetTitle());
218   proj->SetYTitle(hist->GetZaxis()->GetTitle());
219   if (save)
220     proj->Write();
221
222   proj = hist->Project3D("zy");
223   proj->SetXTitle(hist->GetYaxis()->GetTitle());
224   proj->SetYTitle(hist->GetZaxis()->GetTitle());
225   if (save)
226     proj->Write();
227 }
228
229 //____________________________________________________________________
230 void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
231 {
232   // create projections of the 3d hists divides them
233   // axis decides to which plane, if axis is 0 to all planes
234   // the histograms are not returned, just use them from memory or use this to create them in a file
235
236   if (axis == 0)
237   {
238     CreateDividedProjections(hist, hist2, "yx", putErrors, save);
239     CreateDividedProjections(hist, hist2, "zx", putErrors, save);
240     CreateDividedProjections(hist, hist2, "zy", putErrors, save);
241
242     return;
243   }
244
245   TH1* proj = hist->Project3D(axis);
246
247   if (strlen(axis) == 2)
248   {
249     proj->SetYTitle(GetAxisTitle(hist, axis[0]));
250     proj->SetXTitle(GetAxisTitle(hist, axis[1]));
251   }
252   else if (strlen(axis) == 1)
253     proj->SetXTitle(GetAxisTitle(hist, axis[0]));
254
255   TH1* proj2 = hist2->Project3D(axis);
256   if (strlen(axis) == 2)
257   {
258     proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
259     proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
260   }
261   else if (strlen(axis) == 1)
262     proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
263
264   TH1* division = static_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
265   //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()));
266   //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()));
267   division->Divide(proj, proj2, 1, 1, "B");
268   division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
269
270   if (putErrors)
271   {
272     division->Sumw2();
273     if (division->GetDimension() == 1)
274     {
275       Int_t nBins = division->GetNbinsX();
276       for (Int_t i = 1; i <= nBins; ++i)
277         if (proj2->GetBinContent(i) != 0)
278           division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
279     }
280     else if (division->GetDimension() == 2)
281     {
282       Int_t nBinsX = division->GetNbinsX();
283       Int_t nBinsY = division->GetNbinsY();
284       for (Int_t i = 1; i <= nBinsX; ++i)
285         for (Int_t j = 1; j <= nBinsY; ++j)
286           if (proj2->GetBinContent(i, j) != 0)
287             division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
288     }
289   }
290
291   if (save)
292   {
293     proj->Write();
294     proj2->Write();
295     division->Write();
296   }
297 }
298
299 //____________________________________________________________________
300 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
301 {
302   // returns the title of the axis given in axis (x, y, z)
303
304   if (axis == 'x')
305     return hist->GetXaxis()->GetTitle();
306   else if (axis == 'y')
307     return hist->GetYaxis()->GetTitle();
308   else if (axis == 'z')
309     return hist->GetZaxis()->GetTitle();
310
311   return 0;
312 }
313
314
315 AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
316
317   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
318
319   if (!pythiaGenHeader) {
320     printf("AliPWG0Helper::GetProcessType : Unknown gen Header type). \n");
321     return kInvalidProcess;
322   }
323
324
325   Int_t pythiaType = pythiaGenHeader->ProcessType();
326   fgLastProcessType = pythiaType;
327   MCProcessType globalType = kInvalidProcess;  
328
329
330   if (adebug) {
331     printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaType);
332   }
333
334
335   if(pythiaType==92||pythiaType==93){
336     globalType = kSD;
337   }
338   else if(pythiaType==94){
339     globalType = kDD;
340   }
341   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
342   else {
343     globalType = kND;
344   }
345   return globalType;
346 }
347
348
349 AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
350   //
351   // get the process type of the event.
352   //
353
354   // can only read pythia headers, either directly or from cocktalil header
355   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
356
357   if (!dpmJetGenHeader) {
358     printf("AliPWG0Helper::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
359     return kInvalidProcess;
360   }
361
362   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
363   fgLastProcessType = dpmJetType;
364   MCProcessType globalType = kInvalidProcess;  
365
366
367   if (adebug) {
368     printf("AliPWG0Helper::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
369   }
370
371
372   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
373     globalType = kND;
374   }  
375   else if (dpmJetType==5 || dpmJetType==6) {
376     globalType = kSD;
377   }
378   else if (dpmJetType==7) {
379     globalType = kDD;
380   }
381   return globalType;
382 }
383
384
385 AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
386   //
387   // get the process type of the event.
388   //
389
390
391   // Check for simple headers first
392
393   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
394   if (pythiaGenHeader) {
395     return GetPythiaEventProcessType(pythiaGenHeader,adebug);
396   }
397
398   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
399   if (dpmJetGenHeader) {
400     return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
401   }
402   
403
404   // check for cocktail
405
406   AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
407   if (!genCocktailHeader) {
408     printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
409     return kInvalidProcess;
410   }
411
412   TList* headerList = genCocktailHeader->GetHeaders();
413   if (!headerList) {
414     return kInvalidProcess;
415   }
416
417   for (Int_t i=0; i<headerList->GetEntries(); i++) {
418
419     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
420     if (pythiaGenHeader) {
421       return GetPythiaEventProcessType(pythiaGenHeader,adebug);
422     }
423
424     dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
425     if (dpmJetGenHeader) {
426       return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
427     }
428   }
429   return kInvalidProcess;
430 }
431
432
433
434 //____________________________________________________________________
435 TParticle* AliPWG0Helper::FindPrimaryMother(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
442   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
443   if (motherLabel < 0)
444     return 0;
445
446   return stack->Particle(motherLabel);
447 }
448
449 //____________________________________________________________________
450 Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
451 {
452   //
453   // Finds the first mother among the primary particles of the particle identified by <label>,
454   // i.e. the primary that "caused" this particle
455   //
456   // returns its label
457   //
458
459   Int_t nPrim  = stack->GetNprimary();
460
461   while (label >= nPrim)
462   {
463     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
464
465     TParticle* particle = stack->Particle(label);
466     if (!particle)
467     {
468       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
469       return -1;
470     }
471  
472     // find mother
473     if (particle->GetMother(0) < 0)
474     {
475       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
476       return -1;
477     }
478
479     label = particle->GetMother(0);
480   }
481
482   return label;
483 }
484
485 //____________________________________________________________________
486 void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
487 {
488   //
489   // normalizes a 1-d histogram to its bin width
490   //
491
492   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
493   {
494     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
495     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
496   }
497 }
498
499 //____________________________________________________________________
500 void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
501 {
502   //
503   // normalizes a 2-d histogram to its bin width (x width * y width)
504   //
505
506   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
507     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
508     {
509       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
510       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
511       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
512     }
513 }
514
515 //____________________________________________________________________
516 void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment)
517 {
518   //
519   // Prints the given configuration
520   //
521
522   TString str(">>>> Running with >");
523
524   if (analysisMode & kSPD)
525     str += "SPD-only";
526     
527   if (analysisMode & kSPDOnlyL0)
528     str += " (only L0 clusters)";
529   
530   if (analysisMode & kTPC)
531      str += "TPC-only";
532     
533   if (analysisMode & kTPCITS)
534      str += "Global tracking";
535   
536   if (analysisMode & kTPCSPD) 
537     str += "Tracks and tracklets";
538
539   if (analysisMode & kFieldOn)
540   {
541      str += " (with field)";
542   }
543   else
544      str += " (WITHOUT field)";
545   
546   str += "< and trigger >";
547   str += AliTriggerAnalysis::GetTriggerName(trigger);
548   str += "< and diffractive treatment >";
549   
550   switch (diffTreatment)
551   {
552     case kMCFlags:
553       str += "MC flags";
554       break;
555     
556     case kUA5Cuts:
557       str += "UA5 cuts";
558       break;
559        
560     case kE710Cuts:
561       str += "E710 cuts";
562       break;
563     
564     case kALICEHadronLevel:
565       str += "ALICE hadron level";
566       break;
567   }
568   
569   str += "< <<<<";
570
571   Printf("%s", str.Data());
572 }
573
574 //____________________________________________________________________
575 Double_t AliPWG0Helper::Rapidity(Double_t pt, Double_t pz, Double_t m)
576 {
577   //
578   // calculates rapidity keeping the sign in case E == pz
579   //
580
581   Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m);
582   if (energy != TMath::Abs(pz))
583     return 0.5*TMath::Log((energy+pz)/(energy-pz));
584
585   Printf("W- mt=0");
586   return TMath::Sign(1.e30,pz);
587 }
588
589 //____________________________________________________________________
590 Bool_t AliPWG0Helper::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
591 {
592   //
593   // return if process is single diffractive on hadron level
594   // 
595   // xiMax and xiMin cut on M^2/s
596   //
597   // Based on code from Martin Poghoysan
598   //
599   
600   TParticle* part1 = 0;
601   TParticle* part2 = 0;
602   
603   Double_t smallestY = 1e10;
604   Double_t largestY  = -1e10;
605   
606   for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
607   {
608     TParticle* part = stack->Particle(iParticle);
609     if (!part)
610       continue;
611
612     Int_t pdg = TMath::Abs(part->GetPdgCode());
613
614     Int_t child1 = part->GetFirstDaughter();
615     Int_t ist    = part->GetStatusCode();
616
617     Int_t mfl  = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg))));
618     if (child1 > -1 || ist != 1)
619       mfl = 0; // select final state charm and beauty
620     if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4)) 
621       continue;
622     Int_t imother = part->GetFirstMother();
623     if (imother>0)
624     {
625       TParticle *partM = stack->Particle(imother);
626       Int_t pdgM=TMath::Abs(partM->GetPdgCode());
627       if (pdgM==111 || pdgM==3124 || pdgM==3212) 
628         continue;
629     }
630     
631     Double_t y = 0;
632
633     // fix for problem with getting mass of particle 3124
634     if (pdg != 3124)
635       y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
636     else 
637       y = Rapidity(part->Pt(), part->Pz(), 1.5195);
638       
639     if (y < smallestY)
640     {
641       smallestY = y;
642       part1 = part;
643     }
644     
645     if (y > largestY)
646     {
647       largestY = y;
648       part2 = part;
649     }
650   }
651   
652   if (part1 == 0 || part2 == 0)
653     return kFALSE;
654
655   Int_t pdg1 = part1->GetPdgCode();
656   Int_t pdg2 = part2->GetPdgCode();
657
658   Double_t pt1 = part1->Pt();
659   Double_t pt2 = part2->Pt();
660   Double_t pz1 = part1->Pz();
661   Double_t pz2 = part2->Pz();
662   
663   Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
664   Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
665   
666   Int_t arm = -99999;
667   if (pdg1 == 2212 && pdg2 == 2212)
668   {
669     if (y1 > y2) 
670       arm = 0;
671     else
672       arm = 1;
673   }
674   else if (pdg1 == 2212) 
675     arm = 0;
676   else if (pdg2 == 2212) 
677     arm = 1;
678
679   Double_t M02s = 1. - 2 * part1->Energy() / cms;
680   Double_t M12s = 1. - 2 * part2->Energy() / cms;
681
682   if (arm == 0 && M02s > xiMin && M02s < xiMax)
683     return kTRUE;
684   else if (arm == 1 && M12s > xiMin && M12s < xiMax)
685     return kTRUE;
686
687   return kFALSE;
688 }
689
690 //____________________________________________________________________
691 AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, AliPWG0Helper::DiffTreatment diffTreatment)
692 {
693   // 
694   // get process type
695   //
696   // diffTreatment:
697   //   kMCFlags: use MC flags
698   //   kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
699   //   kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND
700   //   kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
701   //
702
703   MCProcessType mcProcessType = GetEventProcessType(header);
704   
705   if (diffTreatment == kMCFlags)
706     return mcProcessType;
707     
708   if (!esd)
709   {
710     Printf("ERROR: AliPWG0Helper::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
711     return kInvalidProcess;
712   }
713     
714   Float_t cms = esd->GetESDRun()->GetBeamEnergy();
715   if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
716     cms *= 2;
717   //Printf("cms = %f", cms);
718
719   if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
720   {
721     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
722       return kSD;
723   }
724   else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
725   {
726     if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
727       return kSD;
728   }
729   else if (diffTreatment == kALICEHadronLevel)
730   {
731     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
732       return kSD;
733   }
734   
735   if (mcProcessType == kSD)
736     return kND;
737     
738   return mcProcessType;
739 }