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