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                                   Double_t&         c,
 290                                   Int_t&            npart,
 291                                   Int_t&            nbin,
 292                                   Double_t&         phiR)
 293 {
 294   //
 295   // Process the event
 296   //
 297   // Parameters:
 298   //   event     Input event
 299   //   triggers  On return, the triggers fired
 300   //   ivz       On return, the found vertex bin (1-based).  A zero
 301   //                  means outside of the defined vertex range
 302   //   vz        On return, the z position of the interaction
 303   //   b         On return, the impact parameter - < 0 if not available
 304   //   c         On return, centrality estimate - < 0 if not available
 305   //   phiR      On return, the event plane angle - < 0 if not available
 306   //
 307   // Return:
 308   //    0 (or kOk) on success, otherwise a bit mask of error codes
 309   //
 310
 311   // Check that we have an event
 312   if (!event) {
 313     AliWarning("No MC event found for input event");
 314     return kNoEvent;
 315   }
 316
 317   //Assign MC only triggers : True NSD etc.
 318   AliHeader*               header          = event->Header();
 319   AliGenEventHeader*       genHeader       = header->GenEventHeader();
 320   AliCollisionGeometry*    colGeometry     =
 321     dynamic_cast<AliCollisionGeometry*>(genHeader);
 322   AliGenPythiaEventHeader* pythiaHeader    =
 323     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 324   AliGenDPMjetEventHeader* dpmHeader       =
 325     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 326   AliGenGeVSimEventHeader* gevHeader       =
 327     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
 328   AliGenHerwigEventHeader* herwigHeader    =
 329     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
 330   // AliGenHijingEventHeader* hijingHeader    =
 331   //   dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 332   // AliGenHydjetEventHeader* hydjetHeader    =
 333   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
 334   // AliGenEposEventHeader*   eposHeader      =
 335   //   dynamic_cast<AliGenEposEventHeader*>(genHeader);
 336
 337   // Check if this is a single diffractive event
 338   Bool_t   sd    = kFALSE;
 339   Double_t phi   = -1111;
 340   npart          = 0;
 341   nbin           = 0;
 342   b              = -1;
 343   c              = -1;
 344   if (colGeometry) {
 345     b     = colGeometry->ImpactParameter();
 346     phi   = colGeometry->ReactionPlaneAngle();
 347     npart = (colGeometry->ProjectileParticipants() +
 348              colGeometry->TargetParticipants());
 349     nbin  = colGeometry->NN();
 350   }
 351   if (fDebug && !colGeometry) {
 352     AliWarningF("Collision header of class %s is not a CollisionHeader",
 353                 genHeader->ClassName());
 354   }
 355
 356   if(pythiaHeader) {
 357     Int_t pythiaType = pythiaHeader->ProcessType();
 358     // 92 and 93 are SD
 359     // 94 is DD
 360     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
 361     b     = pythiaHeader->GetImpactParameter();
 362     npart = 2; // Always 2 protons
 363     nbin  = 1; // Always 1 binary collision
 364   }
 365   if (b >= 0) {
 366     if      (b <  3.5)  c = 2.5; //0-5%
 367     else if (b <  4.95) c = 7.5; //5-10%
 368     else if (b <  6.98) c = 15; //10-20%
 369     else if (b <  8.55) c = 25; //20-30%
 370     else if (b <  9.88) c = 35; //30-40%
 371     else if (b < 11.04) c = 45; //40-50%
 372     else                c = 55; //50-60%
 373   }
 374   if(dpmHeader) { // Also an AliCollisionGeometry
 375     Int_t processType = dpmHeader->ProcessType();
 376     // 1 & 4 are ND
 377     // 5 & 6 are SD
 378     // 7 is DD
 379     if (processType == 5 || processType == 6)  sd = kTRUE;
 380   }
 381   if (gevHeader) {
 382     phi  = gevHeader->GetEventPlane();
 383   }
 384   if (herwigHeader) {
 385     Int_t processType = herwigHeader->ProcessType();
 386     // This is a guess
 387     if (processType == 5 || processType == 6)  sd = kTRUE;
 388     npart = 2; // Always 2 protons
 389     nbin  = 1; // Always 1 binary collision
 390   }
 391   // if (hijingHeader) {
 392   // b    = hijingHeader->ImpactParameter();
 393   // phi  = hijingHeader->ReactionPlaneAngle();
 394   // }
 395   // if (hydjetHeader) {
 396   //   b    = hydjetHeader->ImpactParameter();
 397   //   phi  = hydjetHeader->ReactionPlaneAngle();
 398   // }
 399   // if (eposHeader) {
 400   //   b    = eposHeader->ImpactParameter();
 401   // phi  = eposHeader->ReactionPlaneAngle();
 402   // }
 403
 404   // Normalize event plane angle to [0,2pi]
 405   if (phi <= -1111) phiR = -1;
 406   else {
 407     while (true) {
 408       if      (phi < 0)             phi += 2*TMath::Pi();
 409       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
 410       else                          break;
 411     }
 412   }
 413
 414   // Do a check on particles
 415   sd = IsSingleDiffractive(event->Stack());
 416
 417   // Set NSD flag
 418   if(!sd) {
 419     triggers |= AliAODForwardMult::kMCNSD;
 420     fHTriggers->Fill(kMCNSD+0.5);
 421   }
 422
 423   // Get the primary vertex from EG
 424   TArrayF vtx;
 425   genHeader->PrimaryVertex(vtx);
 426   vz = vtx[2];
 427
 428   if (fDebug) {
 429     AliInfoF("vz=%f, phiR=%f, b=%f, npart=%d, nbin=%d",
 430              vz, phiR, b, npart, nbin);
 431   }
 432   fHVertex->Fill(vz);
 433   fHPhiR->Fill(phiR);
 434   fHB->Fill(b);
 435   fHBvsPart->Fill(b, npart);
 436   fHBvsBin->Fill(b, nbin);
 437
 438   if(fUseDisplacedVertices) {
 439     // Put the vertex at fixed locations
 440     Double_t zvtx  = vz;
 441     Double_t ratio = zvtx/37.5;
 442     if(ratio > 0) ratio = ratio + 0.5;
 443     if(ratio < 0) ratio = ratio - 0.5;
 444     Int_t ratioInt = Int_t(ratio);
 445     zvtx = 37.5*((Double_t)ratioInt);
 446     if(TMath::Abs(zvtx) > 999)
 447       return kBadVertex;
 448   }
 449
 450   // Check for the vertex bin
 451   ivz = fVtxAxis.FindBin(vz);
 452   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
 453     if (fDebug > 3) {
 454       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 455                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
 456     ivz = 0;
 457     return kBadVertex;
 458   }
 459
 460
 461   return kOk;
 462 }
 463
 464 //____________________________________________________________________
 465 namespace {
 466   Double_t rapidity(TParticle* p, Double_t mass)
 467   {
 468     Double_t pt = p->Pt();
 469     Double_t pz = p->Pz();
 470     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
 471     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
 472     return .5 * TMath::Log((ee + pz) / (ee - pz));
 473   }
 474 }
 475
 476 //____________________________________________________________________
 477 Bool_t
 478 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
 479                                             Double_t xiMin,
 480                                             Double_t xiMax) const
 481 {
 482   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
 483   //
 484   // This is re-implemented here to be indendent of the PWG0 library.
 485   TParticle* p1 = 0;     // Particle with least y
 486   TParticle* p2 = 0;     // Particle with largest y
 487   Double_t   y1 =  1e10; // y of p1
 488   Double_t   y2 = -1e10; // y of p2
 489
 490   // Loop over primaries
 491   for (Int_t i = 0; i < stack->GetNprimary(); i++) {
 492     TParticle* p = stack->Particle(i);
 493     if (!p) continue;
 494
 495     Int_t pdg = TMath::Abs(p->GetPdgCode());
 496     Int_t c1  = p->GetFirstDaughter();
 497     Int_t s   = p->GetStatusCode();
 498     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
 499
 500     // Select final state charm and beauty
 501     if (c1 > -1 || s != 1) mfl = 0;
 502
 503     // Check if this is a primary, pi0, Sigma0, ???, or most
 504     // significant digit is larger than or equal to 4
 505     if (!(stack->IsPhysicalPrimary(i) ||
 506           pdg == 111  ||
 507           pdg == 3212 ||
 508           pdg == 3124 ||
 509           mfl >= 4)) continue;
 510
 511     Int_t m1 = p->GetFirstMother();
 512     if (m1 > 0) {
 513       TParticle* pm1  = stack->Particle(m1);
 514       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
 515       // Check if mother is a p0, Simga0, or ???
 516       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
 517     }
 518
 519     // Calculate the rapidity of the particle
 520     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
 521     Double_t yy = rapidity(p, mm);
 522
 523     // Check if the rapidity of this particle is further out than any
 524     // of the preceding particles
 525     if (yy < y1) {
 526       y1 = yy;
 527       p1 = p;
 528     }
 529     if (yy > y2) {
 530       y2 = yy;
 531       p2 = p;
 532     }
 533   }
 534   if (!p1 || !p2) return false;
 535
 536   // Calculate rapidities assuming protons
 537   y1            = TMath::Abs(rapidity(p1, 0.938));
 538   y2            = TMath::Abs(rapidity(p2, 0.938));
 539
 540   // Check if both or just one is a proton
 541   Int_t    pdg1 = p1->GetPdgCode();
 542   Int_t    pdg2 = p2->GetPdgCode();
 543   Int_t    arm  = -99999;
 544   if (pdg1 == 2212 && pdg2 == 2212)
 545     arm = (y1 > y2 ? 0 : 1);
 546   else if (pdg1 == 2212)
 547     arm = 0;
 548   else if (pdg2 == 2212)
 549     arm = 1;
 550   else
 551     return false;
 552
 553   // Rapidity shift
 554   Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0);
 555   Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
 556
 557   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
 558   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
 559
 560   return false;
 561 }
 562
 563 //____________________________________________________________________
 564 Bool_t
 565 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz,
 566                                        Double_t cent,  Double_t b,
 567                                        Int_t    npart, Int_t    nbin)
 568 {
 569   fHVzComp->Fill(trueVz, vz);
 570   fHBvsCent->Fill(b, cent);
 571   fHCentVsPart->Fill(npart, cent);
 572   fHCentVsBin->Fill(nbin, cent);
 573
 574   return true;
 575 }
 576 //
 577 // EOF
 578 //
 579