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)
  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)
  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 = dynamic_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 & kFieldOn)
 537   {
 538      str += " (with field)";
 539   }
 540   else
 541      str += " (WITHOUT field)";
 542
 543   str += "< and trigger >";
 544   str += AliTriggerAnalysis::GetTriggerName(trigger);
 545   str += "< and diffractive treatment >";
 546
 547   switch (diffTreatment)
 548   {
 549     case kMCFlags:
 550       str += "MC flags";
 551       break;
 552
 553     case kUA5Cuts:
 554       str += "UA5 cuts";
 555       break;
 556
 557     case kE710Cuts:
 558       str += "E710 cuts";
 559       break;
 560
 561     case kALICEHadronLevel:
 562       str += "ALICE hadron level";
 563       break;
 564   }
 565
 566   str += "< <<<<";
 567
 568   Printf("%s", str.Data());
 569 }
 570
 571 //____________________________________________________________________
 572 Double_t AliPWG0Helper::Rapidity(Double_t pt, Double_t pz, Double_t m)
 573 {
 574   //
 575   // calculates rapidity keeping the sign in case E == pz
 576   //
 577
 578   Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m);
 579   if (energy != TMath::Abs(pz))
 580     return 0.5*TMath::Log((energy+pz)/(energy-pz));
 581
 582   Printf("W- mt=0");
 583   return TMath::Sign(1.e30,pz);
 584 }
 585
 586 //____________________________________________________________________
 587 Bool_t AliPWG0Helper::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
 588 {
 589   //
 590   // return if process is single diffractive on hadron level
 591   //
 592   // xiMax and xiMin cut on M^2/s
 593   //
 594   // Based on code from Martin Poghoysan
 595   //
 596
 597   TParticle* part1 = 0;
 598   TParticle* part2 = 0;
 599
 600   Double_t smallestY = 1e10;
 601   Double_t largestY  = -1e10;
 602
 603   for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
 604   {
 605     TParticle* part = stack->Particle(iParticle);
 606     if (!part)
 607       continue;
 608
 609     Int_t pdg = TMath::Abs(part->GetPdgCode());
 610
 611     Int_t child1 = part->GetFirstDaughter();
 612     Int_t ist    = part->GetStatusCode();
 613
 614     Int_t mfl  = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg))));
 615     if (child1 > -1 || ist != 1)
 616       mfl = 0; // select final state charm and beauty
 617     if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4))
 618       continue;
 619     Int_t imother = part->GetFirstMother();
 620     if (imother>0)
 621     {
 622       TParticle *partM = stack->Particle(imother);
 623       Int_t pdgM=TMath::Abs(partM->GetPdgCode());
 624       if (pdgM==111 || pdgM==3124 || pdgM==3212)
 625         continue;
 626     }
 627
 628     Double_t y = 0;
 629
 630     // fix for problem with getting mass of particle 3124
 631     if (pdg != 3124)
 632       y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
 633     else
 634       y = Rapidity(part->Pt(), part->Pz(), 1.5195);
 635
 636     if (y < smallestY)
 637     {
 638       smallestY = y;
 639       part1 = part;
 640     }
 641
 642     if (y > largestY)
 643     {
 644       largestY = y;
 645       part2 = part;
 646     }
 647   }
 648
 649   if (part1 == 0 || part2 == 0)
 650     return kFALSE;
 651
 652   Int_t pdg1 = part1->GetPdgCode();
 653   Int_t pdg2 = part2->GetPdgCode();
 654
 655   Double_t pt1 = part1->Pt();
 656   Double_t pt2 = part2->Pt();
 657   Double_t pz1 = part1->Pz();
 658   Double_t pz2 = part2->Pz();
 659
 660   Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
 661   Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
 662
 663   Int_t arm = -99999;
 664   if (pdg1 == 2212 && pdg2 == 2212)
 665   {
 666     if (y1 > y2)
 667       arm = 0;
 668     else
 669       arm = 1;
 670   }
 671   else if (pdg1 == 2212)
 672     arm = 0;
 673   else if (pdg2 == 2212)
 674     arm = 1;
 675
 676   Double_t M02s = 1. - 2 * part1->Energy() / cms;
 677   Double_t M12s = 1. - 2 * part2->Energy() / cms;
 678
 679   if (arm == 0 && M02s > xiMin && M02s < xiMax)
 680     return kTRUE;
 681   else if (arm == 1 && M12s > xiMin && M12s < xiMax)
 682     return kTRUE;
 683
 684   return kFALSE;
 685 }
 686
 687 //____________________________________________________________________
 688 AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, AliPWG0Helper::DiffTreatment diffTreatment)
 689 {
 690   //
 691   // get process type
 692   //
 693   // diffTreatment:
 694   //   kMCFlags: use MC flags
 695   //   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
 696   //   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
 697   //   kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
 698   //
 699
 700   MCProcessType mcProcessType = GetEventProcessType(header);
 701
 702   if (diffTreatment == kMCFlags)
 703     return mcProcessType;
 704
 705   Float_t cms = esd->GetESDRun()->GetBeamEnergy();
 706   if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
 707     cms *= 2;
 708   //Printf("cms = %f", cms);
 709
 710   if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
 711   {
 712     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
 713       return kSD;
 714   }
 715   else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
 716   {
 717     if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
 718       return kSD;
 719   }
 720   else if (diffTreatment == kALICEHadronLevel)
 721   {
 722     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
 723       return kSD;
 724   }
 725
 726   if (mcProcessType == kSD)
 727     return kND;
 728
 729   return mcProcessType;
 730 }