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