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