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  = AllowDisplaced();
 149   if (saveDisplaced) SetVertexMethod(kNormal);
 150   AliFMDEventInspector::SetupForData(vtxAxis);
 151   if (saveDisplaced) SetVertexMethod(kDisplaced);
 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 (AllowDisplaced()) 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    = false;
 361   Double_t phi   = -1111;
 362   Bool_t   egSD  = false;
 363   npart          = 0;
 364   nbin           = 0;
 365   b              = -1;
 366   c              = -1;
 367   phiR           = -1111;
 368   if (colGeometry) {
 369     b     = colGeometry->ImpactParameter();
 370     phi   = colGeometry->ReactionPlaneAngle();
 371     npart = (colGeometry->ProjectileParticipants() +
 372              colGeometry->TargetParticipants());
 373     nbin  = colGeometry->NN();
 374   }
 375   if (fDebug && !colGeometry) {
 376     AliWarningF("Collision header of class %s is not a CollisionHeader",
 377                 genHeader->ClassName());
 378   }
 379
 380   if(pythiaHeader) {
 381     Int_t pythiaType = pythiaHeader->ProcessType();
 382     egSD = true; // We have SD flag in EG
 383     if (pythiaType <= 100) {
 384       // Pythia6
 385       // 92 and 93 are SD
 386       // 94 is DD
 387       if (pythiaType==92 || pythiaType==93) sd = true;
 388     }
 389     else {
 390       // Pythia8
 391       if (pythiaType==103 || pythiaType==104) sd = true;
 392     }
 393     b     = pythiaHeader->GetImpactParameter();
 394     npart = 2; // Always 2 protons
 395     nbin  = 1; // Always 1 binary collision
 396   }
 397   if (b >= 0) {
 398 #if 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 #else
 407     // Updated 4th of November 2014 from
 408     // cern.ch/twiki/bin/view/ALICE/CentStudies#Tables_with_centrality_bins_AN1
 409     Float_t np=0;
 410     Float_t nc=0;
 411     if      (0.00 >= b  && b < 1.57)  { c=0.5;  np=403.8; nc=1861; }
 412     else if (1.57 >= b  && b < 2.22)  { c=1.5;  np=393.6; nc=1766; }
 413     else if (2.22 >= b  && b < 2.71)  { c=2.5;  np=382.9; nc=1678; }
 414     else if (2.71 >= b  && b < 3.13)  { c=3.5;  np=372;   nc=1597; }
 415     else if (3.13 >= b  && b < 3.50)  { c=4.5;  np=361.1; nc=1520; }
 416     else if (3.50 >= b  && b < 4.94)  { c=7.5;  np=329.4; nc=1316; }
 417     else if (4.94 >= b  && b < 6.05)  { c=12.5; np=281.2; nc=1032; }
 418     else if (6.05 >= b  && b < 6.98)  { c=17.5; np=239;   nc=809.8; }
 419     else if (6.98 >= b  && b < 7.81)  { c=22.5; np=202.1; nc=629.6; }
 420     else if (7.81 >= b  && b < 8.55)  { c=27.5; np=169.5; nc=483.7; }
 421     else if (8.55 >= b  && b < 9.23)  { c=32.5; np=141;   nc=366.7; }
 422     else if (9.23 >= b  && b < 9.88)  { c=37.5; np=116;   nc=273.4; }
 423     else if (9.88 >= b  && b < 10.47) { c=42.5; np=94.11; nc=199.4; }
 424     else if (10.47 >= b && b < 11.04) { c=47.5; np=75.3;  nc=143.1; }
 425     else if (11.04 >= b && b < 11.58) { c=52.5; np=59.24; nc=100.1; }
 426     else if (11.58 >= b && b < 12.09) { c=57.5; np=45.58; nc=68.46; }
 427     else if (12.09 >= b && b < 12.58) { c=62.5; np=34.33; nc=45.79; }
 428     else if (12.58 >= b && b < 13.05) { c=67.5; np=25.21; nc=29.92; }
 429     else if (13.05 >= b && b < 13.52) { c=72.5; np=17.96; nc=19.08; }
 430     else if (13.52 >= b && b < 13.97) { c=77.5; np=12.58; nc=12.07; }
 431     else if (13.97 >= b && b < 14.43) { c=82.5; np=8.812; nc=7.682; }
 432     else if (14.43 >= b && b < 14.96) { c=87.5; np=6.158; nc=4.904; }
 433     else if (14.96 >= b && b < 15.67) { c=92.5; np=4.376; nc=3.181; }
 434     else if (15.67 >= b && b < 20.00) { c=97.5; np=3.064; nc=1.994; }
 435     // Be careful to round off
 436     if (npart <= 0) npart = Int_t(np+.5);
 437     if (nbin  <= 0) nbin  = Int_t(nc+.5)/2;
 438 #endif
 439   }
 440   if(dpmHeader) { // Also an AliCollisionGeometry
 441     Int_t processType = dpmHeader->ProcessType();
 442     egSD = true; // We have SD flag in EG
 443     // 1 & 4 are ND
 444     // 5 & 6 are SD
 445     // 7 is DD
 446     if (processType == 5 || processType == 6)  sd = kTRUE;
 447 #if 0
 448       // The below - or rather a different implementation with some
 449       // errors - was proposed by Cvetan - I don't think it's right
 450       // though.  See also
 451       //
 452       //   https://cern.ch/twiki/pub/ALICE/PAPaperCentrality/normalization.pdf
 453       //   https://cern.ch/twiki/bin/view/ALICE/PAMCProductionStudies
 454       //
 455       Int_t nsd1=0, nsd2=0, ndd=0;
 456       Int_t npP = dpm->ProjectileParticipants();
 457       Int_t npT = dpm->TargetParticipants();
 458       // Get the numbeer of single and double diffractive participants
 459       dpm->GetNDiffractive(nsd1,nsd2,ndd);
 460       // Check if all partipants are single/double diffractive
 461       if      ((ndd == 0) && ((npP == nsd1) || (npT == nsd2)))  sd = true;
 462       else if (ndd == (npP + npT))                              dd = true;
 463       // Printf("Projectile: %3d (%3d) Target: %3d (%3d) DD: %3d Process: %d",
 464       //             npP, nsd1, npT, nsd2, ndd, type);
 465 #endif
 466   }
 467   if (gevHeader) {
 468     phi  = gevHeader->GetEventPlane();
 469   }
 470   if (herwigHeader) {
 471     Int_t processType = herwigHeader->ProcessType();
 472     egSD = true; // We have SD flag in EG
 473     // This is a guess
 474     if (processType == 5 || processType == 6)  sd = kTRUE;
 475     npart = 2; // Always 2 protons
 476     nbin  = 1; // Always 1 binary collision
 477   }
 478   // if (hijingHeader) {
 479   // b    = hijingHeader->ImpactParameter();
 480   // phi  = hijingHeader->ReactionPlaneAngle();
 481   // }
 482   // if (hydjetHeader) {
 483   //   b    = hydjetHeader->ImpactParameter();
 484   //   phi  = hydjetHeader->ReactionPlaneAngle();
 485   // }
 486   // if (eposHeader) {
 487   //   b    = eposHeader->ImpactParameter();
 488   // phi  = eposHeader->ReactionPlaneAngle();
 489   // }
 490
 491   // Normalize event plane angle to [0,2pi]
 492   if (phi <= -1111) phiR = -1;
 493   else {
 494     while (true) {
 495       if      (phi < 0)             phi += 2*TMath::Pi();
 496       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
 497       else                          break;
 498     }
 499     phiR = phi;
 500   }
 501
 502   // Do a check on particles - but only if the EG does not give it or
 503   // the impact parameter is large enough
 504   if (!egSD && (b < 0 || b > 15)) sd = IsSingleDiffractive(event->Stack());
 505
 506   // Set NSD flag
 507   if(!sd) {
 508     triggers |= AliAODForwardMult::kMCNSD;
 509     fHTriggers->Fill(kMCNSD+0.5);
 510   }
 511
 512   // Get the primary vertex from EG
 513   TArrayF vtx;
 514   genHeader->PrimaryVertex(vtx);
 515   vz = vtx[2];
 516
 517   DMSG(fDebug, 2, "vz=%f, phiR=%f, b=%f, npart=%d, nbin=%d",
 518        vz, phiR, b, npart, nbin);
 519
 520   fHVertex->Fill(vz);
 521   fHPhiR->Fill(phiR);
 522   fHB->Fill(b);
 523   fHMcC->Fill(c);
 524   fHBvsPart->Fill(b, npart);
 525   fHBvsBin->Fill(b, nbin);
 526
 527   if(AllowDisplaced()) {
 528 #if 0
 529     // Put the vertex at fixed locations
 530     Double_t zvtx  = vz;
 531     Double_t ratio = zvtx/37.5;
 532     if(ratio > 0) ratio = ratio + 0.5;
 533     if(ratio < 0) ratio = ratio - 0.5;
 534     Int_t ratioInt = Int_t(ratio);
 535     zvtx = 37.5*((Double_t)ratioInt);
 536     if(TMath::Abs(zvtx) > 999)
 537       return kBadVertex;
 538 #endif
 539     if (!fDisplacedVertex.ProcessMC(event))
 540       return kBadVertex;
 541     if (fDisplacedVertex.IsSatellite())
 542       vz = fDisplacedVertex.GetVertexZ();
 543   }
 544
 545   // Check for the vertex bin
 546   ivz = fVtxAxis.FindBin(vz);
 547   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
 548     if (fDebug > 3) {
 549       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 550                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
 551     ivz = 0;
 552     return kBadVertex;
 553   }
 554
 555
 556   return kOk;
 557 }
 558 //____________________________________________________________________
 559 Bool_t
 560 AliFMDMCEventInspector::ReadCentrality(const AliESDEvent& esd,
 561                                        Double_t& cent,
 562                                        UShort_t& qual) const
 563 {
 564   //
 565   // Read centrality from event
 566   //
 567   // Parameters:
 568   //    esd  Event
 569   //    cent On return, the centrality or negative if not found
 570   //
 571   // Return:
 572   //    False on error, true otherwise
 573   //
 574   Bool_t ret = AliFMDEventInspector::ReadCentrality(esd, cent, qual);
 575   if (qual != 0) {
 576     AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
 577     if (!centObj)  return ret;
 578
 579     // For MC, we allow `bad' centrality selections
 580     cent = centObj->GetCentralityPercentileUnchecked(fCentMethod);
 581   }
 582   return ret;
 583 }
 584
 585 //____________________________________________________________________
 586 namespace {
 587   Double_t rapidity(TParticle* p, Double_t mass)
 588   {
 589     Double_t pt = p->Pt();
 590     Double_t pz = p->Pz();
 591     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
 592     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
 593     return .5 * TMath::Log((ee + pz) / (ee - pz));
 594   }
 595 }
 596
 597 //____________________________________________________________________
 598 Bool_t
 599 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
 600                                             Double_t xiMin,
 601                                             Double_t xiMax) const
 602 {
 603   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
 604   //
 605   // This is re-implemented here to be indendent of the PWG0 library.
 606   TParticle* p1 = 0;     // Particle with least y
 607   TParticle* p2 = 0;     // Particle with largest y
 608   Double_t   y1 =  1e10; // y of p1
 609   Double_t   y2 = -1e10; // y of p2
 610
 611   // Loop over primaries
 612   for (Int_t i = 0; i < stack->GetNprimary(); i++) {
 613     TParticle* p = stack->Particle(i);
 614     if (!p) continue;
 615
 616     Int_t pdg = TMath::Abs(p->GetPdgCode());
 617     Int_t c1  = p->GetFirstDaughter();
 618     Int_t s   = p->GetStatusCode();
 619     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
 620
 621     // Select final state charm and beauty
 622     if (c1 > -1 || s != 1) mfl = 0;
 623
 624     // Check if this is a primary, pi0, Sigma0, ???, or most
 625     // significant digit is larger than or equal to 4
 626     if (!(stack->IsPhysicalPrimary(i) ||
 627           pdg == 111  ||
 628           pdg == 3212 ||
 629           pdg == 3124 ||
 630           mfl >= 4)) continue;
 631
 632     Int_t m1 = p->GetFirstMother();
 633     if (m1 > 0) {
 634       TParticle* pm1  = stack->Particle(m1);
 635       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
 636       // Check if mother is a p0, Simga0, or ???
 637       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
 638     }
 639
 640     // Calculate the rapidity of the particle
 641     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
 642     Double_t yy = rapidity(p, mm);
 643
 644     // Check if the rapidity of this particle is further out than any
 645     // of the preceding particles
 646     if (yy < y1) {
 647       y1 = yy;
 648       p1 = p;
 649     }
 650     if (yy > y2) {
 651       y2 = yy;
 652       p2 = p;
 653     }
 654   }
 655   if (!p1 || !p2) return false;
 656
 657   // Calculate rapidities assuming protons
 658   y1            = TMath::Abs(rapidity(p1, 0.938));
 659   y2            = TMath::Abs(rapidity(p2, 0.938));
 660
 661   // Check if both or just one is a proton
 662   Int_t    pdg1 = p1->GetPdgCode();
 663   Int_t    pdg2 = p2->GetPdgCode();
 664   Int_t    arm  = -99999;
 665   if (pdg1 == 2212 && pdg2 == 2212)
 666     arm = (y1 > y2 ? 0 : 1);
 667   else if (pdg1 == 2212)
 668     arm = 0;
 669   else if (pdg2 == 2212)
 670     arm = 1;
 671   else
 672     return false;
 673
 674   // Rapidity shift
 675   Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0);
 676   Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
 677
 678   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
 679   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
 680
 681   return false;
 682 }
 683
 684 //____________________________________________________________________
 685 Bool_t
 686 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz,
 687                                        Double_t cent,  Double_t mcC,
 688                                        Double_t b,
 689                                        Int_t    npart, Int_t    nbin)
 690 {
 691   fHVzComp->Fill(trueVz, vz);
 692   fHBvsCent->Fill(b, cent);
 693   fHCentVsPart->Fill(npart, cent);
 694   fHCentVsBin->Fill(nbin, cent);
 695   fHCentVsMcC->Fill(cent, mcC);
 696
 697   return true;
 698 }
 699
 700
 701 //
 702 // EOF
 703 //
 704