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