Coverity 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"
24#include "AliGenPythiaEventHeader.h"
25#include "AliGenDPMjetEventHeader.h"
26#include "AliGenHijingEventHeader.h"
27// #include "AliGenHydjetEventHeader.h"
28#include "AliGenEposEventHeader.h"
29#include "AliGenHerwigEventHeader.h"
30#include "AliGenGeVSimEventHeader.h"
31#include "AliHeader.h"
32#include <TList.h>
e308a636 33#include <TH2F.h>
34
e1f47419 35//====================================================================
36AliFMDMCEventInspector::AliFMDMCEventInspector()
37 : AliFMDEventInspector(),
38 fHVertex(0),
39 fHPhiR(0),
e308a636 40 fHB(0),
41 fHBvsPart(0),
42 fHBvsBin(0),
43 fHBvsCent(0),
44 fHVzComp(0),
45 fHCentVsPart(0),
46 fHCentVsBin(0)
e1f47419 47{
48 //
49 // Constructor
50 //
51}
52
53//____________________________________________________________________
54AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
55 : AliFMDEventInspector("fmdEventInspector"),
56 fHVertex(0),
57 fHPhiR(0),
e308a636 58 fHB(0),
59 fHBvsPart(0),
60 fHBvsBin(0),
61 fHBvsCent(0),
62 fHVzComp(0),
63 fHCentVsPart(0),
64 fHCentVsBin(0)
e1f47419 65{
66 //
67 // Constructor
68 //
69 // Parameters:
70 // name Name of object
71 //
72}
73
74//____________________________________________________________________
75AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
76 : AliFMDEventInspector(o),
77 fHVertex(0),
78 fHPhiR(0),
e308a636 79 fHB(0),
80 fHBvsPart(0),
81 fHBvsBin(0),
82 fHBvsCent(0),
83 fHVzComp(0),
84 fHCentVsPart(0),
85 fHCentVsBin(0)
e1f47419 86{
87 //
88 // Copy constructor
89 //
90 // Parameters:
91 // o Object to copy from
92 //
93}
94
95//____________________________________________________________________
96AliFMDMCEventInspector::~AliFMDMCEventInspector()
97{
98 //
99 // Destructor
100 //
101}
102//____________________________________________________________________
103AliFMDMCEventInspector&
104AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
105{
106 //
107 // Assignement operator
108 //
109 // Parameters:
110 // o Object to assign from
111 //
112 // Return:
113 // Reference to this object
114 //
115 AliFMDEventInspector::operator=(o);
116 return *this;
117}
118
119//____________________________________________________________________
120void
121AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
122{
123 //
124 // Initialize the object
125 //
126 // Parameters:
127 // vtxAxis Vertex axis in use
128 //
129 AliFMDEventInspector::Init(vtxAxis);
130
e308a636 131 Int_t maxPart = 450;
132 Int_t maxBin = 225;
133 Int_t maxB = 25;
134 Int_t nVtx = vtxAxis.GetNbins();
135 Double_t lVtx = vtxAxis.GetXmin();
136 Double_t hVtx = vtxAxis.GetXmax();
137 fHVertex = new TH1F("vertex", "True vertex distribution", nVtx, lVtx, hVtx);
e1f47419 138 fHVertex->SetXTitle("v_{z} [cm]");
139 fHVertex->SetYTitle("# of events");
140 fHVertex->SetFillColor(kGreen+1);
141 fHVertex->SetFillStyle(3001);
142 fHVertex->SetDirectory(0);
143 // fHVertex->Sumw2();
144 fList->Add(fHVertex);
145
146 fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
147 fHPhiR->SetXTitle("#Phi_{R} [radians]");
148 fHPhiR->SetYTitle("# of events");
149 fHPhiR->SetFillColor(kGreen+1);
150 fHPhiR->SetFillStyle(3001);
151 fHPhiR->SetDirectory(0);
152 fList->Add(fHPhiR);
153
e308a636 154 fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
e1f47419 155 fHB->SetXTitle("b [fm]");
156 fHB->SetYTitle("# of events");
157 fHB->SetFillColor(kGreen+1);
158 fHB->SetFillStyle(3001);
159 fHB->SetDirectory(0);
160 fList->Add(fHB);
e308a636 161
162 fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
163 5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
164 fHBvsPart->SetXTitle("b [fm]");
165 fHBvsPart->SetYTitle("# of participants");
166 fHBvsPart->SetZTitle("Events");
167 fHBvsPart->SetDirectory(0);
168 fList->Add(fHBvsPart);
169
170 fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
171 5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
172 fHBvsBin->SetXTitle("b [fm]");
173 fHBvsBin->SetYTitle("# of binary collisions");
174 fHBvsBin->SetZTitle("Events");
175 fHBvsBin->SetDirectory(0);
176 fList->Add(fHBvsBin);
e1f47419 177
e308a636 178 fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
179 5*maxB, 0, maxB, fCentAxis->GetNbins(),
180 fCentAxis->GetXbins()->GetArray());
181 fHBvsCent->SetXTitle("b [fm]");
182 fHBvsCent->SetYTitle("Centrality [%]");
183 fHBvsCent->SetZTitle("Event");
184 fHBvsCent->SetDirectory(0);
185 fList->Add(fHBvsCent);
186
187
188 fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
189 10*nVtx, lVtx, hVtx, 10*nVtx, lVtx, hVtx);
190 fHVzComp->SetXTitle("True v_{z} [cm]");
191 fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
192 fHVzComp->SetZTitle("Events");
193 fHVzComp->SetDirectory(0);
194 fList->Add(fHVzComp);
195
196 fHCentVsPart = new TH2F("centralityVsParticipans",
197 "# of participants vs Centrality",
198 maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(),
199 fCentAxis->GetXbins()->GetArray());
200 fHCentVsPart->SetXTitle("Participants");
201 fHCentVsPart->SetYTitle("Centrality [%]");
202 fHCentVsPart->SetZTitle("Event");
203 fHCentVsPart->SetDirectory(0);
204 fList->Add(fHCentVsPart);
205
206 fHCentVsBin = new TH2F("centralityVsBinary",
207 "# of binary collisions vs Centrality",
208 maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(),
209 fCentAxis->GetXbins()->GetArray());
210 fHCentVsBin->SetXTitle("Binary collisions");
211 fHCentVsBin->SetYTitle("Centrality [%]");
212 fHCentVsBin->SetZTitle("Event");
213 fHCentVsBin->SetDirectory(0);
214 fList->Add(fHCentVsBin);
e1f47419 215}
216
217//____________________________________________________________________
218UInt_t
219AliFMDMCEventInspector::ProcessMC(AliMCEvent* event,
220 UInt_t& triggers,
221 UShort_t& ivz,
222 Double_t& vz,
223 Double_t& b,
e308a636 224 Int_t& npart,
225 Int_t& nbin,
e1f47419 226 Double_t& phiR)
227{
228 //
229 // Process the event
230 //
231 // Parameters:
232 // event Input event
233 // triggers On return, the triggers fired
234 // ivz On return, the found vertex bin (1-based). A zero
235 // means outside of the defined vertex range
236 // vz On return, the z position of the interaction
237 // b On return, the impact parameter - < 0 if not available
238 // phiR On return, the event plane angle - < 0 if not available
239 //
240 // Return:
241 // 0 (or kOk) on success, otherwise a bit mask of error codes
242 //
243
244 // Check that we have an event
245 if (!event) {
246 AliWarning("No MC event found for input event");
247 return kNoEvent;
248 }
249
250 //Assign MC only triggers : True NSD etc.
251 AliHeader* header = event->Header();
252 AliGenEventHeader* genHeader = header->GenEventHeader();
e308a636 253 AliCollisionGeometry* colGeometry =
254 dynamic_cast<AliCollisionGeometry*>(genHeader);
e1f47419 255 AliGenPythiaEventHeader* pythiaHeader =
256 dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
257 AliGenDPMjetEventHeader* dpmHeader =
258 dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
259 AliGenGeVSimEventHeader* gevHeader =
260 dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
e1f47419 261 AliGenHerwigEventHeader* herwigHeader =
262 dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
e308a636 263 // AliGenHijingEventHeader* hijingHeader =
264 // dynamic_cast<AliGenHijingEventHeader*>(genHeader);
e1f47419 265 // AliGenHydjetEventHeader* hydjetHeader =
266 // dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
e308a636 267 // AliGenEposEventHeader* eposHeader =
268 // dynamic_cast<AliGenEposEventHeader*>(genHeader);
e1f47419 269
270 // Check if this is a single diffractive event
e308a636 271 Bool_t sd = kFALSE;
272 Double_t phi = -1111;
273 npart = 0;
274 nbin = 0;
275 b = -1;
276 if (colGeometry) {
277 b = colGeometry->ImpactParameter();
278 phi = colGeometry->ReactionPlaneAngle();
279 npart = (colGeometry->ProjectileParticipants() +
280 colGeometry->TargetParticipants());
281 nbin = colGeometry->NN();
282 }
e1f47419 283 if(pythiaHeader) {
284 Int_t pythiaType = pythiaHeader->ProcessType();
285 if (pythiaType==92 || pythiaType==93) sd = kTRUE;
e308a636 286 b = pythiaHeader->GetImpactParameter();
287 npart = 2; // Always 2 protons
288 nbin = 1; // Always 1 binary collision
e1f47419 289 }
e308a636 290 if(dpmHeader) { // Also an AliCollisionGeometry
e1f47419 291 Int_t processType = dpmHeader->ProcessType();
292 if (processType == 5 || processType == 6) sd = kTRUE;
e1f47419 293 }
294 if (gevHeader) {
295 phi = gevHeader->GetEventPlane();
296 }
e1f47419 297 if (herwigHeader) {
298 Int_t processType = herwigHeader->ProcessType();
299 // This is a guess
300 if (processType == 5 || processType == 6) sd = kTRUE;
e308a636 301 npart = 2; // Always 2 protons
302 nbin = 1; // Always 1 binary collision
e1f47419 303 }
e308a636 304 // if (hijingHeader) {
305 // b = hijingHeader->ImpactParameter();
306 // phi = hijingHeader->ReactionPlaneAngle();
307 // }
e1f47419 308 // if (hydjetHeader) {
309 // b = hydjetHeader->ImpactParameter();
310 // phi = hydjetHeader->ReactionPlaneAngle();
311 // }
e308a636 312 // if (eposHeader) {
313 // b = eposHeader->ImpactParameter();
314 // phi = eposHeader->ReactionPlaneAngle();
315 // }
e1f47419 316
317 // Normalize event plane angle to [0,2pi]
318 if (phi <= -1111) phiR = -1;
319 else {
320 while (true) {
321 if (phi < 0) phi += 2*TMath::Pi();
322 else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
323 else break;
324 }
325 }
326
327 // Set NSD flag
328 if(!sd) {
329 triggers |= AliAODForwardMult::kMCNSD;
330 fHTriggers->Fill(kMCNSD+0.5);
331 }
332
333 // Get the primary vertex from EG
334 TArrayF vtx;
335 genHeader->PrimaryVertex(vtx);
336 vz = vtx[2];
337
338 fHVertex->Fill(vz);
339 fHPhiR->Fill(phiR);
340 fHB->Fill(b);
e308a636 341 fHBvsPart->Fill(b, npart);
342 fHBvsBin->Fill(b, nbin);
e1f47419 343
344 // Check for the vertex bin
345 ivz = fHEventsTr->GetXaxis()->FindBin(vz);
346 if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
347 if (fDebug > 3) {
348 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
349 vz, fHEventsTr->GetXaxis()->GetXmin(),
350 fHEventsTr->GetXaxis()->GetXmax())); }
351 ivz = 0;
352 return kBadVertex;
353 }
354
355
356 return kOk;
357}
358//____________________________________________________________________
359Bool_t
e308a636 360AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd,
361 Double_t& cent,
362 UShort_t& qual) const
e1f47419 363{
364 //
365 // Read centrality from event
366 //
367 // Parameters:
368 // esd Event
369 // cent On return, the centrality or negative if not found
370 //
371 // Return:
372 // False on error, true otherwise
373 //
e308a636 374 cent = -1;
375 qual = 0;
e1f47419 376 AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
e308a636 377 if (!centObj) return true;
378
379 qual = centObj->GetQuality();
380 if (qual == 0x8) // Ignore ZDC outliers
381 cent = centObj->GetCentralityPercentileUnchecked("V0M");
382 else
383 cent = centObj->GetCentralityPercentile("V0M");
e1f47419 384
385 return true;
386}
387
e308a636 388//____________________________________________________________________
389Bool_t
390AliFMDMCEventInspector::CompareResults(Double_t vz, Double_t trueVz,
391 Double_t cent, Double_t b,
392 Int_t npart, Int_t nbin)
393{
394 fHVzComp->Fill(trueVz, vz);
395 fHBvsCent->Fill(b, cent);
396 fHCentVsPart->Fill(npart, cent);
397 fHCentVsBin->Fill(nbin, cent);
398
399 return true;
400}
e1f47419 401//
402// EOF
403//
404