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