]> 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
Mior fixes
[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),
49 fHCentVsBin(0)
e1f47419 50{
51 //
52 // Constructor
53 //
54}
55
56//____________________________________________________________________
57AliFMDMCEventInspector::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//____________________________________________________________________
78AliFMDMCEventInspector::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//____________________________________________________________________
99AliFMDMCEventInspector::~AliFMDMCEventInspector()
100{
101 //
102 // Destructor
103 //
104}
105//____________________________________________________________________
106AliFMDMCEventInspector&
107AliFMDMCEventInspector::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//____________________________________________________________________
123void
124AliFMDMCEventInspector::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//____________________________________________________________________
221UInt_t
222AliFMDMCEventInspector::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
d2286e26 356 ivz = fVtxAxis.FindBin(vz);
e1f47419 357 if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
358 if (fDebug > 3) {
359 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
d2286e26 360 vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
e1f47419 361 ivz = 0;
362 return kBadVertex;
363 }
364
365
366 return kOk;
367}
1f480471 368
369//____________________________________________________________________
370namespace {
371 Double_t rapidity(TParticle* p, Double_t mass)
372 {
373 Double_t pt = p->Pt();
374 Double_t pz = p->Pz();
375 Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
376 if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
377 return .5 * TMath::Log((ee + pz) / (ee - pz));
378 }
379}
380
381//____________________________________________________________________
382Bool_t
383AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
384 Double_t xiMin,
385 Double_t xiMax) const
386{
387 // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
388 //
389 // This is re-implemented here to be indendent of the PWG0 library.
390 TParticle* p1 = 0; // Particle with least y
391 TParticle* p2 = 0; // Particle with largest y
392 Double_t y1 = 1e10; // y of p1
393 Double_t y2 = -1e10; // y of p2
394
395 // Loop over primaries
396 for (Int_t i = 0; i < stack->GetNprimary(); i++) {
397 TParticle* p = stack->Particle(i);
398 if (!p) continue;
399
400 Int_t pdg = TMath::Abs(p->GetPdgCode());
401 Int_t c1 = p->GetFirstDaughter();
402 Int_t s = p->GetStatusCode();
403 Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
404
405 // Select final state charm and beauty
406 if (c1 > -1 || s != 1) mfl = 0;
407
408 // Check if this is a primary, pi0, Sigma0, ???, or most
409 // significant digit is larger than or equal to 4
410 if (!(stack->IsPhysicalPrimary(i) ||
411 pdg == 111 ||
412 pdg == 3212 ||
413 pdg == 3124 ||
414 mfl >= 4)) continue;
415
416 Int_t m1 = p->GetFirstMother();
417 if (m1 > 0) {
418 TParticle* pm1 = stack->Particle(m1);
419 Int_t mpdg = TMath::Abs(pm1->GetPdgCode());
420 // Check if mother is a p0, Simga0, or ???
421 if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
422 }
423
424 // Calculate the rapidity of the particle
425 Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
426 Double_t yy = rapidity(p, mm);
427
428 // Check if the rapidity of this particle is further out than any
429 // of the preceding particles
430 if (yy < y1) {
431 y1 = yy;
432 p1 = p;
433 }
434 if (yy > y2) {
435 y2 = yy;
436 p2 = p;
437 }
438 }
439 if (!p1 || !p2) return false;
440
441 // Calculate rapidities assuming protons
442 y1 = TMath::Abs(rapidity(p1, 0.938));
443 y2 = TMath::Abs(rapidity(p2, 0.938));
444
445 // Check if both or just one is a proton
446 Int_t pdg1 = p1->GetPdgCode();
447 Int_t pdg2 = p2->GetPdgCode();
448 Int_t arm = -99999;
449 if (pdg1 == 2212 && pdg2 == 2212)
450 arm = (y1 > y2 ? 0 : 1);
451 else if (pdg1 == 2212)
452 arm = 0;
453 else if (pdg2 == 2212)
454 arm = 1;
455 else
456 return false;
457
5bb5d1f6 458 // Rapidity shift
1f480471 459 Double_t m02s = 1 - 2 * p1->Energy() / fEnergy;
460 Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
461
462 if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
463 if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
464
465 return false;
466}
e1f47419 467//____________________________________________________________________
468Bool_t
e308a636 469AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd,
470 Double_t& cent,
471 UShort_t& qual) const
e1f47419 472{
473 //
474 // Read centrality from event
475 //
476 // Parameters:
477 // esd Event
478 // cent On return, the centrality or negative if not found
479 //
480 // Return:
481 // False on error, true otherwise
482 //
e308a636 483 cent = -1;
484 qual = 0;
e1f47419 485 AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
e308a636 486 if (!centObj) return true;
487
488 qual = centObj->GetQuality();
489 if (qual == 0x8) // Ignore ZDC outliers
490 cent = centObj->GetCentralityPercentileUnchecked("V0M");
491 else
492 cent = centObj->GetCentralityPercentile("V0M");
e1f47419 493
494 return true;
495}
496
e308a636 497//____________________________________________________________________
498Bool_t
499AliFMDMCEventInspector::CompareResults(Double_t vz, Double_t trueVz,
500 Double_t cent, Double_t b,
501 Int_t npart, Int_t nbin)
502{
503 fHVzComp->Fill(trueVz, vz);
504 fHBvsCent->Fill(b, cent);
505 fHCentVsPart->Fill(npart, cent);
506 fHCentVsBin->Fill(nbin, cent);
507
508 return true;
509}
e1f47419 510//
511// EOF
512//
513