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