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