]>
Commit | Line | Data |
---|---|---|
e1f47419 | 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" | |
1f480471 | 24 | #include "AliStack.h" |
e1f47419 | 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> | |
e308a636 | 34 | #include <TH2F.h> |
1f480471 | 35 | #include <TParticle.h> |
36 | #include <TMath.h> | |
e308a636 | 37 | |
e1f47419 | 38 | //==================================================================== |
39 | AliFMDMCEventInspector::AliFMDMCEventInspector() | |
40 | : AliFMDEventInspector(), | |
41 | fHVertex(0), | |
42 | fHPhiR(0), | |
e308a636 | 43 | fHB(0), |
44 | fHBvsPart(0), | |
45 | fHBvsBin(0), | |
46 | fHBvsCent(0), | |
47 | fHVzComp(0), | |
48 | fHCentVsPart(0), | |
49 | fHCentVsBin(0) | |
e1f47419 | 50 | { |
51 | // | |
52 | // Constructor | |
53 | // | |
54 | } | |
55 | ||
56 | //____________________________________________________________________ | |
57 | AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */) | |
58 | : AliFMDEventInspector("fmdEventInspector"), | |
59 | fHVertex(0), | |
60 | fHPhiR(0), | |
e308a636 | 61 | fHB(0), |
62 | fHBvsPart(0), | |
63 | fHBvsBin(0), | |
64 | fHBvsCent(0), | |
65 | fHVzComp(0), | |
66 | fHCentVsPart(0), | |
67 | fHCentVsBin(0) | |
e1f47419 | 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), | |
e308a636 | 82 | fHB(0), |
83 | fHBvsPart(0), | |
84 | fHBvsBin(0), | |
85 | fHBvsCent(0), | |
86 | fHVzComp(0), | |
87 | fHCentVsPart(0), | |
88 | fHCentVsBin(0) | |
e1f47419 | 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 | ||
e308a636 | 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); | |
e1f47419 | 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 | ||
e308a636 | 157 | fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB); |
e1f47419 | 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); | |
e308a636 | 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); | |
e1f47419 | 180 | |
e308a636 | 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); | |
e1f47419 | 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, | |
e308a636 | 227 | Int_t& npart, |
228 | Int_t& nbin, | |
e1f47419 | 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(); | |
e308a636 | 256 | AliCollisionGeometry* colGeometry = |
257 | dynamic_cast<AliCollisionGeometry*>(genHeader); | |
e1f47419 | 258 | AliGenPythiaEventHeader* pythiaHeader = |
259 | dynamic_cast<AliGenPythiaEventHeader*>(genHeader); | |
260 | AliGenDPMjetEventHeader* dpmHeader = | |
261 | dynamic_cast<AliGenDPMjetEventHeader*>(genHeader); | |
262 | AliGenGeVSimEventHeader* gevHeader = | |
263 | dynamic_cast<AliGenGeVSimEventHeader*>(genHeader); | |
e1f47419 | 264 | AliGenHerwigEventHeader* herwigHeader = |
265 | dynamic_cast<AliGenHerwigEventHeader*>(genHeader); | |
e308a636 | 266 | // AliGenHijingEventHeader* hijingHeader = |
267 | // dynamic_cast<AliGenHijingEventHeader*>(genHeader); | |
e1f47419 | 268 | // AliGenHydjetEventHeader* hydjetHeader = |
269 | // dynamic_cast<AliGenHydjetEventHeader*>(genHeader); | |
e308a636 | 270 | // AliGenEposEventHeader* eposHeader = |
271 | // dynamic_cast<AliGenEposEventHeader*>(genHeader); | |
e1f47419 | 272 | |
273 | // Check if this is a single diffractive event | |
e308a636 | 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 | } | |
e1f47419 | 286 | if(pythiaHeader) { |
287 | Int_t pythiaType = pythiaHeader->ProcessType(); | |
1f480471 | 288 | // 92 and 93 are SD |
289 | // 94 is DD | |
e1f47419 | 290 | if (pythiaType==92 || pythiaType==93) sd = kTRUE; |
e308a636 | 291 | b = pythiaHeader->GetImpactParameter(); |
292 | npart = 2; // Always 2 protons | |
293 | nbin = 1; // Always 1 binary collision | |
e1f47419 | 294 | } |
e308a636 | 295 | if(dpmHeader) { // Also an AliCollisionGeometry |
e1f47419 | 296 | Int_t processType = dpmHeader->ProcessType(); |
1f480471 | 297 | // 1 & 4 are ND |
298 | // 5 & 6 are SD | |
299 | // 7 is DD | |
e1f47419 | 300 | if (processType == 5 || processType == 6) sd = kTRUE; |
e1f47419 | 301 | } |
302 | if (gevHeader) { | |
303 | phi = gevHeader->GetEventPlane(); | |
304 | } | |
e1f47419 | 305 | if (herwigHeader) { |
306 | Int_t processType = herwigHeader->ProcessType(); | |
307 | // This is a guess | |
308 | if (processType == 5 || processType == 6) sd = kTRUE; | |
e308a636 | 309 | npart = 2; // Always 2 protons |
310 | nbin = 1; // Always 1 binary collision | |
e1f47419 | 311 | } |
e308a636 | 312 | // if (hijingHeader) { |
313 | // b = hijingHeader->ImpactParameter(); | |
314 | // phi = hijingHeader->ReactionPlaneAngle(); | |
315 | // } | |
e1f47419 | 316 | // if (hydjetHeader) { |
317 | // b = hydjetHeader->ImpactParameter(); | |
318 | // phi = hydjetHeader->ReactionPlaneAngle(); | |
319 | // } | |
e308a636 | 320 | // if (eposHeader) { |
321 | // b = eposHeader->ImpactParameter(); | |
322 | // phi = eposHeader->ReactionPlaneAngle(); | |
323 | // } | |
e1f47419 | 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 | ||
1f480471 | 335 | // Do a check on particles |
336 | sd = IsSingleDiffractive(event->Stack()); | |
337 | ||
e1f47419 | 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); | |
e308a636 | 352 | fHBvsPart->Fill(b, npart); |
353 | fHBvsBin->Fill(b, nbin); | |
e1f47419 | 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 | } | |
1f480471 | 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 | ||
5bb5d1f6 | 459 | // Rapidity shift |
1f480471 | 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 | } | |
e1f47419 | 468 | //____________________________________________________________________ |
469 | Bool_t | |
e308a636 | 470 | AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, |
471 | Double_t& cent, | |
472 | UShort_t& qual) const | |
e1f47419 | 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 | // | |
e308a636 | 484 | cent = -1; |
485 | qual = 0; | |
e1f47419 | 486 | AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality(); |
e308a636 | 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"); | |
e1f47419 | 494 | |
495 | return true; | |
496 | } | |
497 | ||
e308a636 | 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 | } | |
e1f47419 | 511 | // |
512 | // EOF | |
513 | // | |
514 |