]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity (Ruben)
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCEventInspector.cxx
CommitLineData
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//====================================================================
39AliFMDMCEventInspector::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),
f7cfc454 49 fHCentVsBin(0),
50 fProduction("")
e1f47419 51{
52 //
53 // Constructor
54 //
55}
56
57//____________________________________________________________________
58AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
59 : AliFMDEventInspector("fmdEventInspector"),
60 fHVertex(0),
61 fHPhiR(0),
e308a636 62 fHB(0),
63 fHBvsPart(0),
64 fHBvsBin(0),
65 fHBvsCent(0),
66 fHVzComp(0),
67 fHCentVsPart(0),
f7cfc454 68 fHCentVsBin(0),
69 fProduction("")
e1f47419 70{
71 //
72 // Constructor
73 //
74 // Parameters:
75 // name Name of object
76 //
77}
78
79//____________________________________________________________________
80AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
81 : AliFMDEventInspector(o),
82 fHVertex(0),
83 fHPhiR(0),
e308a636 84 fHB(0),
85 fHBvsPart(0),
86 fHBvsBin(0),
87 fHBvsCent(0),
88 fHVzComp(0),
89 fHCentVsPart(0),
f7cfc454 90 fHCentVsBin(0),
91 fProduction("")
e1f47419 92{
93 //
94 // Copy constructor
95 //
96 // Parameters:
97 // o Object to copy from
98 //
99}
100
101//____________________________________________________________________
102AliFMDMCEventInspector::~AliFMDMCEventInspector()
103{
104 //
105 // Destructor
106 //
107}
108//____________________________________________________________________
109AliFMDMCEventInspector&
110AliFMDMCEventInspector::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//____________________________________________________________________
126void
127AliFMDMCEventInspector::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
e308a636 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);
e1f47419 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
e308a636 160 fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
e1f47419 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);
e308a636 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);
e1f47419 183
e308a636 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);
e1f47419 221}
222
f7cfc454 223//____________________________________________________________________
224void
225AliFMDMCEventInspector::StoreInformation()
226{
227 // Store information about running conditions in the output list
228 if (!fList) return;
229 AliFMDEventInspector::StoreInformation();
230 TNamed* mc = new TNamed("mc", fProduction.Data());
231 mc->SetUniqueID(1);
232 fList->Add(mc);
233}
234
235namespace
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//____________________________________________________________________
247void
248AliFMDMCEventInspector::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
e1f47419 280//____________________________________________________________________
281UInt_t
282AliFMDMCEventInspector::ProcessMC(AliMCEvent* event,
283 UInt_t& triggers,
284 UShort_t& ivz,
285 Double_t& vz,
286 Double_t& b,
e308a636 287 Int_t& npart,
288 Int_t& nbin,
e1f47419 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();
e308a636 316 AliCollisionGeometry* colGeometry =
317 dynamic_cast<AliCollisionGeometry*>(genHeader);
e1f47419 318 AliGenPythiaEventHeader* pythiaHeader =
319 dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
320 AliGenDPMjetEventHeader* dpmHeader =
321 dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
322 AliGenGeVSimEventHeader* gevHeader =
323 dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
e1f47419 324 AliGenHerwigEventHeader* herwigHeader =
325 dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
e308a636 326 // AliGenHijingEventHeader* hijingHeader =
327 // dynamic_cast<AliGenHijingEventHeader*>(genHeader);
e1f47419 328 // AliGenHydjetEventHeader* hydjetHeader =
329 // dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
e308a636 330 // AliGenEposEventHeader* eposHeader =
331 // dynamic_cast<AliGenEposEventHeader*>(genHeader);
e1f47419 332
333 // Check if this is a single diffractive event
e308a636 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 }
e1f47419 346 if(pythiaHeader) {
347 Int_t pythiaType = pythiaHeader->ProcessType();
1f480471 348 // 92 and 93 are SD
349 // 94 is DD
e1f47419 350 if (pythiaType==92 || pythiaType==93) sd = kTRUE;
e308a636 351 b = pythiaHeader->GetImpactParameter();
352 npart = 2; // Always 2 protons
353 nbin = 1; // Always 1 binary collision
e1f47419 354 }
e308a636 355 if(dpmHeader) { // Also an AliCollisionGeometry
e1f47419 356 Int_t processType = dpmHeader->ProcessType();
1f480471 357 // 1 & 4 are ND
358 // 5 & 6 are SD
359 // 7 is DD
e1f47419 360 if (processType == 5 || processType == 6) sd = kTRUE;
e1f47419 361 }
362 if (gevHeader) {
363 phi = gevHeader->GetEventPlane();
364 }
e1f47419 365 if (herwigHeader) {
366 Int_t processType = herwigHeader->ProcessType();
367 // This is a guess
368 if (processType == 5 || processType == 6) sd = kTRUE;
e308a636 369 npart = 2; // Always 2 protons
370 nbin = 1; // Always 1 binary collision
e1f47419 371 }
e308a636 372 // if (hijingHeader) {
373 // b = hijingHeader->ImpactParameter();
374 // phi = hijingHeader->ReactionPlaneAngle();
375 // }
e1f47419 376 // if (hydjetHeader) {
377 // b = hydjetHeader->ImpactParameter();
378 // phi = hydjetHeader->ReactionPlaneAngle();
379 // }
e308a636 380 // if (eposHeader) {
381 // b = eposHeader->ImpactParameter();
382 // phi = eposHeader->ReactionPlaneAngle();
383 // }
e1f47419 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
1f480471 395 // Do a check on particles
396 sd = IsSingleDiffractive(event->Stack());
397
e1f47419 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);
e308a636 412 fHBvsPart->Fill(b, npart);
413 fHBvsBin->Fill(b, nbin);
e1f47419 414
415 // Check for the vertex bin
d2286e26 416 ivz = fVtxAxis.FindBin(vz);
e1f47419 417 if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
418 if (fDebug > 3) {
419 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
d2286e26 420 vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
e1f47419 421 ivz = 0;
422 return kBadVertex;
423 }
424
425
426 return kOk;
427}
1f480471 428
429//____________________________________________________________________
430namespace {
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//____________________________________________________________________
442Bool_t
443AliFMDMCEventInspector::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
5bb5d1f6 518 // Rapidity shift
1f480471 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}
e1f47419 527//____________________________________________________________________
528Bool_t
e308a636 529AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd,
530 Double_t& cent,
531 UShort_t& qual) const
e1f47419 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 //
e308a636 543 cent = -1;
544 qual = 0;
e1f47419 545 AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
e308a636 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");
e1f47419 553
554 return true;
555}
556
e308a636 557//____________________________________________________________________
558Bool_t
559AliFMDMCEventInspector::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}
e1f47419 570//
571// EOF
572//
573