2 // This class inspects the event
5 // - AliESDFMD object possibly corrected for sharing
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
12 // Note, that these are added to the master output list
17 #include "AliFMDMCEventInspector.h"
19 #include "AliAODForwardMult.h"
20 #include "AliForwardUtil.h"
21 #include "AliCentrality.h"
22 #include "AliESDEvent.h"
23 #include "AliMCEvent.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"
35 #include <TParticle.h>
38 //====================================================================
39 AliFMDMCEventInspector::AliFMDMCEventInspector()
40 : AliFMDEventInspector(),
57 //____________________________________________________________________
58 AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
59 : AliFMDEventInspector("fmdEventInspector"),
75 // name Name of object
79 //____________________________________________________________________
80 AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
81 : AliFMDEventInspector(o),
97 // o Object to copy from
101 //____________________________________________________________________
102 AliFMDMCEventInspector::~AliFMDMCEventInspector()
108 //____________________________________________________________________
109 AliFMDMCEventInspector&
110 AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
113 // Assignement operator
116 // o Object to assign from
119 // Reference to this object
121 AliFMDEventInspector::operator=(o);
125 //____________________________________________________________________
127 AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
130 // Initialize the object
133 // vtxAxis Vertex axis in use
135 AliFMDEventInspector::Init(vtxAxis);
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);
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);
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);
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);
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);
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);
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);
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);
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);
223 //____________________________________________________________________
225 AliFMDMCEventInspector::StoreInformation(Int_t runNo)
227 // Store information about running conditions in the output list
229 AliFMDEventInspector::StoreInformation(runNo);
230 TNamed* mc = new TNamed("mc", fProduction.Data());
237 TString& AppendField(TString& s, bool test, const char* f)
240 if (!s.IsNull()) s.Append(",");
246 //____________________________________________________________________
248 AliFMDMCEventInspector::ReadProductionDetails(AliMCEvent* event)
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);
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");
280 //____________________________________________________________________
282 AliFMDMCEventInspector::ProcessMC(AliMCEvent* 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
304 // 0 (or kOk) on success, otherwise a bit mask of error codes
307 // Check that we have an event
309 AliWarning("No MC event found for input event");
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);
333 // Check if this is a single diffractive event
335 Double_t phi = -1111;
340 b = colGeometry->ImpactParameter();
341 phi = colGeometry->ReactionPlaneAngle();
342 npart = (colGeometry->ProjectileParticipants() +
343 colGeometry->TargetParticipants());
344 nbin = colGeometry->NN();
347 Int_t pythiaType = pythiaHeader->ProcessType();
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
355 if(dpmHeader) { // Also an AliCollisionGeometry
356 Int_t processType = dpmHeader->ProcessType();
360 if (processType == 5 || processType == 6) sd = kTRUE;
363 phi = gevHeader->GetEventPlane();
366 Int_t processType = herwigHeader->ProcessType();
368 if (processType == 5 || processType == 6) sd = kTRUE;
369 npart = 2; // Always 2 protons
370 nbin = 1; // Always 1 binary collision
372 // if (hijingHeader) {
373 // b = hijingHeader->ImpactParameter();
374 // phi = hijingHeader->ReactionPlaneAngle();
376 // if (hydjetHeader) {
377 // b = hydjetHeader->ImpactParameter();
378 // phi = hydjetHeader->ReactionPlaneAngle();
381 // b = eposHeader->ImpactParameter();
382 // phi = eposHeader->ReactionPlaneAngle();
385 // Normalize event plane angle to [0,2pi]
386 if (phi <= -1111) phiR = -1;
389 if (phi < 0) phi += 2*TMath::Pi();
390 else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
395 // Do a check on particles
396 sd = IsSingleDiffractive(event->Stack());
400 triggers |= AliAODForwardMult::kMCNSD;
401 fHTriggers->Fill(kMCNSD+0.5);
404 // Get the primary vertex from EG
406 genHeader->PrimaryVertex(vtx);
412 fHBvsPart->Fill(b, npart);
413 fHBvsBin->Fill(b, nbin);
415 // Check for the vertex bin
416 ivz = fVtxAxis.FindBin(vz);
417 if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
419 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
420 vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
429 //____________________________________________________________________
431 Double_t rapidity(TParticle* p, Double_t mass)
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));
441 //____________________________________________________________________
443 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
445 Double_t xiMax) const
447 // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
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
455 // Loop over primaries
456 for (Int_t i = 0; i < stack->GetNprimary(); i++) {
457 TParticle* p = stack->Particle(i);
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))));
465 // Select final state charm and beauty
466 if (c1 > -1 || s != 1) mfl = 0;
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) ||
476 Int_t m1 = p->GetFirstMother();
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;
484 // Calculate the rapidity of the particle
485 Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
486 Double_t yy = rapidity(p, mm);
488 // Check if the rapidity of this particle is further out than any
489 // of the preceding particles
499 if (!p1 || !p2) return false;
501 // Calculate rapidities assuming protons
502 y1 = TMath::Abs(rapidity(p1, 0.938));
503 y2 = TMath::Abs(rapidity(p2, 0.938));
505 // Check if both or just one is a proton
506 Int_t pdg1 = p1->GetPdgCode();
507 Int_t pdg2 = p2->GetPdgCode();
509 if (pdg1 == 2212 && pdg2 == 2212)
510 arm = (y1 > y2 ? 0 : 1);
511 else if (pdg1 == 2212)
513 else if (pdg2 == 2212)
519 Double_t m02s = 1 - 2 * p1->Energy() / fEnergy;
520 Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
522 if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
523 if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
527 //____________________________________________________________________
529 AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd,
531 UShort_t& qual) const
534 // Read centrality from event
538 // cent On return, the centrality or negative if not found
541 // False on error, true otherwise
545 AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
546 if (!centObj) return true;
548 qual = centObj->GetQuality();
549 if (qual == 0x8) // Ignore ZDC outliers
550 cent = centObj->GetCentralityPercentileUnchecked("V0M");
552 cent = centObj->GetCentralityPercentile("V0M");
557 //____________________________________________________________________
559 AliFMDMCEventInspector::CompareResults(Double_t vz, Double_t trueVz,
560 Double_t cent, Double_t b,
561 Int_t npart, Int_t nbin)
563 fHVzComp->Fill(trueVz, vz);
564 fHBvsCent->Fill(b, cent);
565 fHCentVsPart->Fill(npart, cent);
566 fHCentVsBin->Fill(nbin, cent);