PWGLF/FORWARD/analysis2/AliFMDMCEventInspector.cxx

   1 //
   2 // This class inspects the event
   3 //
   4 // Input:
   5 //   - AliESDFMD object possibly corrected for sharing
   6 //
   7 // Output:
   8 //   - A histogram of v_z of events with triggers.
   9 //   - A histogram of v_z of events with vertex and triggers
  10 //   - A histogram of trigger counters
  11 //
  12 // Note, that these are added to the master output list
  13 //
  14 // Corrections used:
  15 //   - None
  16 //
  17 #include "AliFMDMCEventInspector.h"
  18 #include "AliLog.h"
  19 #include "AliAODForwardMult.h"
  20 #include "AliForwardUtil.h"
  21 #include "AliCentrality.h"
  22 #include "AliESDEvent.h"
  23 #include "AliMCEvent.h"
  24 #include "AliStack.h"
  25 #include "AliGenPythiaEventHeader.h"
  26 #include "AliGenDPMjetEventHeader.h"
  27 #include "AliGenHijingEventHeader.h"
  28 // #include "AliGenHydjetEventHeader.h"
  29 #include "AliGenEposEventHeader.h"
  30 #include "AliGenHerwigEventHeader.h"
  31 #include "AliGenGeVSimEventHeader.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliMCEventHandler.h"
  34 #include "AliHeader.h"
  35 #include <TList.h>
  36 #include <TH2F.h>
  37 #include <TParticle.h>
  38 #include <TMath.h>
  39 #include <TParameter.h>
  40
  41 //====================================================================
  42 AliFMDMCEventInspector::AliFMDMCEventInspector()
  43   : AliFMDEventInspector(),
  44     fHVertex(0),
  45     fHPhiR(0),
  46     fHB(0),
  47     fHMcC(0),
  48     fHBvsPart(0),
  49     fHBvsBin(0),
  50     fHBvsCent(0),
  51     fHVzComp(0),
  52     fHCentVsPart(0),
  53     fHCentVsBin(0),
  54     fHCentVsMcC(0),
  55     fProduction("")
  56 {
  57   //
  58   // Constructor
  59   //
  60 }
  61
  62 //____________________________________________________________________
  63 AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
  64   : AliFMDEventInspector("fmdEventInspector"),
  65     fHVertex(0),
  66     fHPhiR(0),
  67     fHB(0),
  68     fHMcC(0),
  69     fHBvsPart(0),
  70     fHBvsBin(0),
  71     fHBvsCent(0),
  72     fHVzComp(0),
  73     fHCentVsPart(0),
  74     fHCentVsBin(0),
  75     fHCentVsMcC(0),
  76     fProduction("")
  77 {
  78   //
  79   // Constructor
  80   //
  81   // Parameters:
  82   //   name Name of object
  83   //
  84 }
  85
  86 //____________________________________________________________________
  87 AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
  88   : AliFMDEventInspector(o),
  89     fHVertex(0),
  90     fHPhiR(0),
  91     fHB(0),
  92     fHMcC(0),
  93     fHBvsPart(0),
  94     fHBvsBin(0),
  95     fHBvsCent(0),
  96     fHVzComp(0),
  97     fHCentVsPart(0),
  98     fHCentVsBin(0),
  99     fHCentVsMcC(0),
 100     fProduction("")
 101 {
 102   //
 103   // Copy constructor
 104   //
 105   // Parameters:
 106   //   o Object to copy from
 107   //
 108 }
 109
 110 //____________________________________________________________________
 111 AliFMDMCEventInspector::~AliFMDMCEventInspector()
 112 {
 113   //
 114   // Destructor
 115   //
 116 }
 117 //____________________________________________________________________
 118 AliFMDMCEventInspector&
 119 AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
 120 {
 121   //
 122   // Assignement operator
 123   //
 124   // Parameters:
 125   //   o Object to assign from
 126   //
 127   // Return:
 128   //    Reference to this object
 129   //
 130   AliFMDEventInspector::operator=(o);
 131   return *this;
 132 }
 133
 134 //____________________________________________________________________
 135 void
 136 AliFMDMCEventInspector::SetupForData(const TAxis& vtxAxis)
 137 {
 138   //
 139   // Initialize the object
 140   //
 141   // Parameters:
 142   //   vtxAxis Vertex axis in use
 143   //
 144
 145   // We temporary disable the displaced vertices so we can initialize
 146   // the routine ourselves.
 147   Bool_t saveDisplaced  = fUseDisplacedVertices;
 148   fUseDisplacedVertices = false;
 149   AliFMDEventInspector::SetupForData(vtxAxis);
 150   fUseDisplacedVertices = saveDisplaced;
 151
 152   Int_t    maxPart = 450;
 153   Int_t    maxBin  = 225;
 154   Int_t    maxB    = 25;
 155   Int_t    nVtx = vtxAxis.GetNbins();
 156   Double_t lVtx = vtxAxis.GetXmin();
 157   Double_t hVtx = vtxAxis.GetXmax();
 158   fHVertex = new TH1F("vertex", "True vertex distribution", nVtx, lVtx, hVtx);
 159   fHVertex->SetXTitle("v_{z} [cm]");
 160   fHVertex->SetYTitle("# of events");
 161   fHVertex->SetFillColor(kGreen+1);
 162   fHVertex->SetFillStyle(3001);
 163   fHVertex->SetDirectory(0);
 164   // fHVertex->Sumw2();
 165   fList->Add(fHVertex);
 166
 167   fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
 168   fHPhiR->SetXTitle("#Phi_{R} [radians]");
 169   fHPhiR->SetYTitle("# of events");
 170   fHPhiR->SetFillColor(kGreen+1);
 171   fHPhiR->SetFillStyle(3001);
 172   fHPhiR->SetDirectory(0);
 173   fList->Add(fHPhiR);
 174
 175   fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
 176   fHB->SetXTitle("b [fm]");
 177   fHB->SetYTitle("# of events");
 178   fHB->SetFillColor(kGreen+1);
 179   fHB->SetFillStyle(3001);
 180   fHB->SetDirectory(0);
 181   fList->Add(fHB);
 182
 183   fHMcC = static_cast<TH1F*>(fHCent->Clone("mcC"));
 184   fHMcC->SetFillColor(kCyan+2);
 185   fHMcC->SetDirectory(0);
 186   fList->Add(fHMcC);
 187
 188   fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
 189                        5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
 190   fHBvsPart->SetXTitle("b [fm]");
 191   fHBvsPart->SetYTitle("# of participants");
 192   fHBvsPart->SetZTitle("Events");
 193   fHBvsPart->SetDirectory(0);
 194   fList->Add(fHBvsPart);
 195
 196   fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
 197                        5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
 198   fHBvsBin->SetXTitle("b [fm]");
 199   fHBvsBin->SetYTitle("# of binary collisions");
 200   fHBvsBin->SetZTitle("Events");
 201   fHBvsBin->SetDirectory(0);
 202   fList->Add(fHBvsBin);
 203
 204   fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
 205                        5*maxB, 0, maxB, fCentAxis->GetNbins(),
 206                        fCentAxis->GetXbins()->GetArray());
 207   fHBvsCent->SetXTitle("b [fm]");
 208   fHBvsCent->SetYTitle("Centrality [%]");
 209   fHBvsCent->SetZTitle("Event");
 210   fHBvsCent->SetDirectory(0);
 211   fList->Add(fHBvsCent);
 212
 213
 214   fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
 215                       10*nVtx, lVtx, hVtx, 10*nVtx, lVtx, hVtx);
 216   fHVzComp->SetXTitle("True v_{z} [cm]");
 217   fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
 218   fHVzComp->SetZTitle("Events");
 219   fHVzComp->SetDirectory(0);
 220   fList->Add(fHVzComp);
 221
 222   fHCentVsPart = new TH2F("centralityVsParticipans",
 223                           "# of participants vs Centrality",
 224                           maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(),
 225                           fCentAxis->GetXbins()->GetArray());
 226   fHCentVsPart->SetXTitle("Participants");
 227   fHCentVsPart->SetYTitle("Centrality [%]");
 228   fHCentVsPart->SetZTitle("Event");
 229   fHCentVsPart->SetDirectory(0);
 230   fList->Add(fHCentVsPart);
 231
 232   fHCentVsBin = new TH2F("centralityVsBinary",
 233                          "# of binary collisions vs Centrality",
 234                          maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(),
 235                          fCentAxis->GetXbins()->GetArray());
 236   fHCentVsBin->SetXTitle("Binary collisions");
 237   fHCentVsBin->SetYTitle("Centrality [%]");
 238   fHCentVsBin->SetZTitle("Event");
 239   fHCentVsBin->SetDirectory(0);
 240   fList->Add(fHCentVsBin);
 241
 242   Int_t    nC = fHCent->GetNbinsX();
 243   Double_t cL = fHCent->GetXaxis()->GetXmin();
 244   Double_t cH = fHCent->GetXaxis()->GetXmax();
 245   fHCentVsMcC = new TH2F("centralityRecoVsMC",
 246                          "Centrality from reconstruction vs MC derived",
 247                          nC, cL, cH, nC, cL, cH);
 248   fHCentVsMcC->SetDirectory(0);
 249   fHCentVsMcC->SetStats(0);
 250   fHCentVsMcC->SetXTitle("Centralty from Reco [%]");
 251   fHCentVsMcC->SetYTitle("Centralty derived from Impact Par. [%]");
 252   fHCentVsMcC->SetZTitle("Events");
 253   fList->Add(fHCentVsMcC);
 254
 255   if (fUseDisplacedVertices) fDisplacedVertex.SetupForData(fList, "", true);
 256 }
 257
 258 //____________________________________________________________________
 259 void AliFMDMCEventInspector::StoreInformation()
 260 {
 261   // Store information about running conditions in the output list
 262   if (!fList) return;
 263   AliFMDEventInspector::StoreInformation();
 264   Bool_t mc = true;
 265   fList->Add(AliForwardUtil::MakeParameter("mc", mc));
 266   // , fProduction.Data());
 267 }
 268
 269 namespace
 270 {
 271   TString& AppendField(TString& s, bool test, const char* f)
 272   {
 273     if (!test) return s;
 274     if (!s.IsNull()) s.Append(",");
 275     s.Append(f);
 276     return s;
 277   }
 278 }
 279
 280 //____________________________________________________________________
 281 void
 282 AliFMDMCEventInspector::ReadProductionDetails(AliMCEvent* event)
 283 {
 284   //Assign MC only triggers : True NSD etc.
 285   AliHeader*               header          = event->Header();
 286   AliGenEventHeader*       genHeader       = header->GenEventHeader();
 287   AliCollisionGeometry*    colGeometry     =
 288     dynamic_cast<AliCollisionGeometry*>(genHeader);
 289   AliGenPythiaEventHeader* pythiaHeader    =
 290     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 291   AliGenDPMjetEventHeader* dpmHeader       =
 292     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 293   AliGenGeVSimEventHeader* gevHeader       =
 294     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
 295   AliGenHerwigEventHeader* herwigHeader    =
 296     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
 297   AliGenHijingEventHeader* hijingHeader    =
 298     dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 299   // AliGenHydjetEventHeader* hydjetHeader    =
 300   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
 301   AliGenEposEventHeader*   eposHeader      =
 302     dynamic_cast<AliGenEposEventHeader*>(genHeader);
 303
 304   AppendField(fProduction, colGeometry,  "Geometry");
 305   AppendField(fProduction, pythiaHeader, "Pythia");
 306   AppendField(fProduction, dpmHeader,    "DPM");
 307   AppendField(fProduction, gevHeader,    "GevSim");
 308   AppendField(fProduction, herwigHeader, "Herwig");
 309   AppendField(fProduction, hijingHeader, "Hijing");
 310   // AppendField(fProduction, hydjetHeader, "Hydjet");
 311   AppendField(fProduction, eposHeader,   "EPOS");
 312 }
 313
 314 //____________________________________________________________________
 315 UInt_t
 316 AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
 317                                   UInt_t&           triggers,
 318                                   UShort_t&         ivz,
 319                                   Double_t&         vz,
 320                                   Double_t&         b,
 321                                   Double_t&         c,
 322                                   Int_t&            npart,
 323                                   Int_t&            nbin,
 324                                   Double_t&         phiR)
 325 {
 326   //
 327   // Process the event
 328   //
 329   // Parameters:
 330   //   event     Input event
 331   //   triggers  On return, the triggers fired
 332   //   ivz       On return, the found vertex bin (1-based).  A zero
 333   //                  means outside of the defined vertex range
 334   //   vz        On return, the z position of the interaction
 335   //   b         On return, the impact parameter - < 0 if not available
 336   //   c         On return, centrality estimate - < 0 if not available
 337   //   phiR      On return, the event plane angle - < 0 if not available
 338   //
 339   // Return:
 340   //    0 (or kOk) on success, otherwise a bit mask of error codes
 341   //
 342
 343   // Check that we have an event
 344   if (!event) {
 345     AliWarning("No MC event found for input event");
 346     return kNoEvent;
 347   }
 348
 349   //Assign MC only triggers : True NSD etc.
 350   AliHeader*               header          = event->Header();
 351   AliGenEventHeader*       genHeader       = header->GenEventHeader();
 352   AliCollisionGeometry*    colGeometry     =
 353     dynamic_cast<AliCollisionGeometry*>(genHeader);
 354   AliGenPythiaEventHeader* pythiaHeader    =
 355     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 356   AliGenDPMjetEventHeader* dpmHeader       =
 357     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 358   AliGenGeVSimEventHeader* gevHeader       =
 359     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
 360   AliGenHerwigEventHeader* herwigHeader    =
 361     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
 362   // AliGenHijingEventHeader* hijingHeader    =
 363   //   dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 364   // AliGenHydjetEventHeader* hydjetHeader    =
 365   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
 366   // AliGenEposEventHeader*   eposHeader      =
 367   //   dynamic_cast<AliGenEposEventHeader*>(genHeader);
 368
 369   // Check if this is a single diffractive event
 370   Bool_t   sd    = kFALSE;
 371   Double_t phi   = -1111;
 372   npart          = 0;
 373   nbin           = 0;
 374   b              = -1;
 375   c              = -1;
 376   phiR           = -1111;
 377   if (colGeometry) {
 378     b     = colGeometry->ImpactParameter();
 379     phi   = colGeometry->ReactionPlaneAngle();
 380     npart = (colGeometry->ProjectileParticipants() +
 381              colGeometry->TargetParticipants());
 382     nbin  = colGeometry->NN();
 383   }
 384   if (fDebug && !colGeometry) {
 385     AliWarningF("Collision header of class %s is not a CollisionHeader",
 386                 genHeader->ClassName());
 387   }
 388
 389   if(pythiaHeader) {
 390     Int_t pythiaType = pythiaHeader->ProcessType();
 391     // 92 and 93 are SD
 392     // 94 is DD
 393     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
 394     b     = pythiaHeader->GetImpactParameter();
 395     npart = 2; // Always 2 protons
 396     nbin  = 1; // Always 1 binary collision
 397   }
 398   if (b >= 0) {
 399     if      (b <  3.5)  c = 2.5; //0-5%
 400     else if (b <  4.95) c = 7.5; //5-10%
 401     else if (b <  6.98) c = 15; //10-20%
 402     else if (b <  8.55) c = 25; //20-30%
 403     else if (b <  9.88) c = 35; //30-40%
 404     else if (b < 11.04) c = 45; //40-50%
 405     else                c = 55; //50-60%
 406   }
 407   if(dpmHeader) { // Also an AliCollisionGeometry
 408     Int_t processType = dpmHeader->ProcessType();
 409     // 1 & 4 are ND
 410     // 5 & 6 are SD
 411     // 7 is DD
 412     if (processType == 5 || processType == 6)  sd = kTRUE;
 413   }
 414   if (gevHeader) {
 415     phi  = gevHeader->GetEventPlane();
 416   }
 417   if (herwigHeader) {
 418     Int_t processType = herwigHeader->ProcessType();
 419     // This is a guess
 420     if (processType == 5 || processType == 6)  sd = kTRUE;
 421     npart = 2; // Always 2 protons
 422     nbin  = 1; // Always 1 binary collision
 423   }
 424   // if (hijingHeader) {
 425   // b    = hijingHeader->ImpactParameter();
 426   // phi  = hijingHeader->ReactionPlaneAngle();
 427   // }
 428   // if (hydjetHeader) {
 429   //   b    = hydjetHeader->ImpactParameter();
 430   //   phi  = hydjetHeader->ReactionPlaneAngle();
 431   // }
 432   // if (eposHeader) {
 433   //   b    = eposHeader->ImpactParameter();
 434   // phi  = eposHeader->ReactionPlaneAngle();
 435   // }
 436
 437   // Normalize event plane angle to [0,2pi]
 438   if (phi <= -1111) phiR = -1;
 439   else {
 440     while (true) {
 441       if      (phi < 0)             phi += 2*TMath::Pi();
 442       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
 443       else                          break;
 444     }
 445     phiR = phi;
 446   }
 447
 448   // Do a check on particles
 449   sd = IsSingleDiffractive(event->Stack());
 450
 451   // Set NSD flag
 452   if(!sd) {
 453     triggers |= AliAODForwardMult::kMCNSD;
 454     fHTriggers->Fill(kMCNSD+0.5);
 455   }
 456
 457   // Get the primary vertex from EG
 458   TArrayF vtx;
 459   genHeader->PrimaryVertex(vtx);
 460   vz = vtx[2];
 461
 462   DMSG(fDebug, 2, "vz=%f, phiR=%f, b=%f, npart=%d, nbin=%d",
 463        vz, phiR, b, npart, nbin);
 464
 465   fHVertex->Fill(vz);
 466   fHPhiR->Fill(phiR);
 467   fHB->Fill(b);
 468   fHMcC->Fill(c);
 469   fHBvsPart->Fill(b, npart);
 470   fHBvsBin->Fill(b, nbin);
 471
 472   if(fUseDisplacedVertices) {
 473 #if 0
 474     // Put the vertex at fixed locations
 475     Double_t zvtx  = vz;
 476     Double_t ratio = zvtx/37.5;
 477     if(ratio > 0) ratio = ratio + 0.5;
 478     if(ratio < 0) ratio = ratio - 0.5;
 479     Int_t ratioInt = Int_t(ratio);
 480     zvtx = 37.5*((Double_t)ratioInt);
 481     if(TMath::Abs(zvtx) > 999)
 482       return kBadVertex;
 483 #endif
 484     if (!fDisplacedVertex.ProcessMC(event))
 485       return kBadVertex;
 486     if (fDisplacedVertex.IsSatellite())
 487       vz = fDisplacedVertex.GetVertexZ();
 488   }
 489
 490   // Check for the vertex bin
 491   ivz = fVtxAxis.FindBin(vz);
 492   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
 493     if (fDebug > 3) {
 494       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 495                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
 496     ivz = 0;
 497     return kBadVertex;
 498   }
 499
 500
 501   return kOk;
 502 }
 503 //____________________________________________________________________
 504 Bool_t
 505 AliFMDMCEventInspector::ReadCentrality(const AliESDEvent& esd,
 506                                        Double_t& cent,
 507                                        UShort_t& qual) const
 508 {
 509   //
 510   // Read centrality from event
 511   //
 512   // Parameters:
 513   //    esd  Event
 514   //    cent On return, the centrality or negative if not found
 515   //
 516   // Return:
 517   //    False on error, true otherwise
 518   //
 519   Bool_t ret = AliFMDEventInspector::ReadCentrality(esd, cent, qual);
 520   if (qual != 0) {
 521     AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
 522     if (!centObj)  return ret;
 523
 524     // For MC, we allow `bad' centrality selections
 525     cent = centObj->GetCentralityPercentileUnchecked(fCentMethod);
 526   }
 527   return ret;
 528 }
 529
 530 //____________________________________________________________________
 531 namespace {
 532   Double_t rapidity(TParticle* p, Double_t mass)
 533   {
 534     Double_t pt = p->Pt();
 535     Double_t pz = p->Pz();
 536     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
 537     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
 538     return .5 * TMath::Log((ee + pz) / (ee - pz));
 539   }
 540 }
 541
 542 //____________________________________________________________________
 543 Bool_t
 544 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
 545                                             Double_t xiMin,
 546                                             Double_t xiMax) const
 547 {
 548   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
 549   //
 550   // This is re-implemented here to be indendent of the PWG0 library.
 551   TParticle* p1 = 0;     // Particle with least y
 552   TParticle* p2 = 0;     // Particle with largest y
 553   Double_t   y1 =  1e10; // y of p1
 554   Double_t   y2 = -1e10; // y of p2
 555
 556   // Loop over primaries
 557   for (Int_t i = 0; i < stack->GetNprimary(); i++) {
 558     TParticle* p = stack->Particle(i);
 559     if (!p) continue;
 560
 561     Int_t pdg = TMath::Abs(p->GetPdgCode());
 562     Int_t c1  = p->GetFirstDaughter();
 563     Int_t s   = p->GetStatusCode();
 564     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
 565
 566     // Select final state charm and beauty
 567     if (c1 > -1 || s != 1) mfl = 0;
 568
 569     // Check if this is a primary, pi0, Sigma0, ???, or most
 570     // significant digit is larger than or equal to 4
 571     if (!(stack->IsPhysicalPrimary(i) ||
 572           pdg == 111  ||
 573           pdg == 3212 ||
 574           pdg == 3124 ||
 575           mfl >= 4)) continue;
 576
 577     Int_t m1 = p->GetFirstMother();
 578     if (m1 > 0) {
 579       TParticle* pm1  = stack->Particle(m1);
 580       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
 581       // Check if mother is a p0, Simga0, or ???
 582       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
 583     }
 584
 585     // Calculate the rapidity of the particle
 586     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
 587     Double_t yy = rapidity(p, mm);
 588
 589     // Check if the rapidity of this particle is further out than any
 590     // of the preceding particles
 591     if (yy < y1) {
 592       y1 = yy;
 593       p1 = p;
 594     }
 595     if (yy > y2) {
 596       y2 = yy;
 597       p2 = p;
 598     }
 599   }
 600   if (!p1 || !p2) return false;
 601
 602   // Calculate rapidities assuming protons
 603   y1            = TMath::Abs(rapidity(p1, 0.938));
 604   y2            = TMath::Abs(rapidity(p2, 0.938));
 605
 606   // Check if both or just one is a proton
 607   Int_t    pdg1 = p1->GetPdgCode();
 608   Int_t    pdg2 = p2->GetPdgCode();
 609   Int_t    arm  = -99999;
 610   if (pdg1 == 2212 && pdg2 == 2212)
 611     arm = (y1 > y2 ? 0 : 1);
 612   else if (pdg1 == 2212)
 613     arm = 0;
 614   else if (pdg2 == 2212)
 615     arm = 1;
 616   else
 617     return false;
 618
 619   // Rapidity shift
 620   Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0);
 621   Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
 622
 623   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
 624   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
 625
 626   return false;
 627 }
 628
 629 //____________________________________________________________________
 630 Bool_t
 631 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz,
 632                                        Double_t cent,  Double_t mcC,
 633                                        Double_t b,
 634                                        Int_t    npart, Int_t    nbin)
 635 {
 636   fHVzComp->Fill(trueVz, vz);
 637   fHBvsCent->Fill(b, cent);
 638   fHCentVsPart->Fill(npart, cent);
 639   fHCentVsBin->Fill(nbin, cent);
 640   fHCentVsMcC->Fill(cent, mcC);
 641
 642   return true;
 643 }
 644
 645
 646 //
 647 // EOF
 648 //
 649