T0 QA Task
[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>
33//====================================================================
34AliFMDMCEventInspector::AliFMDMCEventInspector()
35 : AliFMDEventInspector(),
36 fHVertex(0),
37 fHPhiR(0),
38 fHB(0)
39{
40 //
41 // Constructor
42 //
43}
44
45//____________________________________________________________________
46AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
47 : AliFMDEventInspector("fmdEventInspector"),
48 fHVertex(0),
49 fHPhiR(0),
50 fHB(0)
51{
52 //
53 // Constructor
54 //
55 // Parameters:
56 // name Name of object
57 //
58}
59
60//____________________________________________________________________
61AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
62 : AliFMDEventInspector(o),
63 fHVertex(0),
64 fHPhiR(0),
65 fHB(0)
66{
67 //
68 // Copy constructor
69 //
70 // Parameters:
71 // o Object to copy from
72 //
73}
74
75//____________________________________________________________________
76AliFMDMCEventInspector::~AliFMDMCEventInspector()
77{
78 //
79 // Destructor
80 //
81}
82//____________________________________________________________________
83AliFMDMCEventInspector&
84AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
85{
86 //
87 // Assignement operator
88 //
89 // Parameters:
90 // o Object to assign from
91 //
92 // Return:
93 // Reference to this object
94 //
95 AliFMDEventInspector::operator=(o);
96 return *this;
97}
98
99//____________________________________________________________________
100void
101AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
102{
103 //
104 // Initialize the object
105 //
106 // Parameters:
107 // vtxAxis Vertex axis in use
108 //
109 AliFMDEventInspector::Init(vtxAxis);
110
111 fHVertex = new TH1F("vertex", "True vertex distribution",
112 vtxAxis.GetNbins(),
113 vtxAxis.GetXmin(),
114 vtxAxis.GetXmax());
115 fHVertex->SetXTitle("v_{z} [cm]");
116 fHVertex->SetYTitle("# of events");
117 fHVertex->SetFillColor(kGreen+1);
118 fHVertex->SetFillStyle(3001);
119 fHVertex->SetDirectory(0);
120 // fHVertex->Sumw2();
121 fList->Add(fHVertex);
122
123 fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
124 fHPhiR->SetXTitle("#Phi_{R} [radians]");
125 fHPhiR->SetYTitle("# of events");
126 fHPhiR->SetFillColor(kGreen+1);
127 fHPhiR->SetFillStyle(3001);
128 fHPhiR->SetDirectory(0);
129 fList->Add(fHPhiR);
130
131 fHB = new TH1F("b", "Impact parameter", 125, 0, 25);
132 fHB->SetXTitle("b [fm]");
133 fHB->SetYTitle("# of events");
134 fHB->SetFillColor(kGreen+1);
135 fHB->SetFillStyle(3001);
136 fHB->SetDirectory(0);
137 fList->Add(fHB);
138
139}
140
141//____________________________________________________________________
142UInt_t
143AliFMDMCEventInspector::ProcessMC(AliMCEvent* event,
144 UInt_t& triggers,
145 UShort_t& ivz,
146 Double_t& vz,
147 Double_t& b,
148 Double_t& phiR)
149{
150 //
151 // Process the event
152 //
153 // Parameters:
154 // event Input event
155 // triggers On return, the triggers fired
156 // ivz On return, the found vertex bin (1-based). A zero
157 // means outside of the defined vertex range
158 // vz On return, the z position of the interaction
159 // b On return, the impact parameter - < 0 if not available
160 // phiR On return, the event plane angle - < 0 if not available
161 //
162 // Return:
163 // 0 (or kOk) on success, otherwise a bit mask of error codes
164 //
165
166 // Check that we have an event
167 if (!event) {
168 AliWarning("No MC event found for input event");
169 return kNoEvent;
170 }
171
172 //Assign MC only triggers : True NSD etc.
173 AliHeader* header = event->Header();
174 AliGenEventHeader* genHeader = header->GenEventHeader();
175 AliGenPythiaEventHeader* pythiaHeader =
176 dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
177 AliGenDPMjetEventHeader* dpmHeader =
178 dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
179 AliGenGeVSimEventHeader* gevHeader =
180 dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
181 AliGenHijingEventHeader* hijingHeader =
182 dynamic_cast<AliGenHijingEventHeader*>(genHeader);
183 AliGenHerwigEventHeader* herwigHeader =
184 dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
185 // AliGenHydjetEventHeader* hydjetHeader =
186 // dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
187 AliGenEposEventHeader* eposHeader =
188 dynamic_cast<AliGenEposEventHeader*>(genHeader);
189
190 // Check if this is a single diffractive event
191 Bool_t sd = kFALSE;
192 Double_t phi = -1111;
193 b = -1;
194 if(pythiaHeader) {
195 Int_t pythiaType = pythiaHeader->ProcessType();
196 if (pythiaType==92 || pythiaType==93) sd = kTRUE;
197 b = pythiaHeader->GetImpactParameter();
198 }
199 if(dpmHeader) {
200 Int_t processType = dpmHeader->ProcessType();
201 if (processType == 5 || processType == 6) sd = kTRUE;
202 b = dpmHeader->ImpactParameter();
203 phi = dpmHeader->ReactionPlaneAngle();
204
205 }
206 if (gevHeader) {
207 phi = gevHeader->GetEventPlane();
208 }
209 if (hijingHeader) {
210 b = hijingHeader->ImpactParameter();
211 phi = hijingHeader->ReactionPlaneAngle();
212 }
213 if (herwigHeader) {
214 Int_t processType = herwigHeader->ProcessType();
215 // This is a guess
216 if (processType == 5 || processType == 6) sd = kTRUE;
217 }
218 // if (hydjetHeader) {
219 // b = hydjetHeader->ImpactParameter();
220 // phi = hydjetHeader->ReactionPlaneAngle();
221 // }
222 if (eposHeader) {
223 b = eposHeader->ImpactParameter();
224 phi = eposHeader->ReactionPlaneAngle();
225 }
226
227 // Normalize event plane angle to [0,2pi]
228 if (phi <= -1111) phiR = -1;
229 else {
230 while (true) {
231 if (phi < 0) phi += 2*TMath::Pi();
232 else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
233 else break;
234 }
235 }
236
237 // Set NSD flag
238 if(!sd) {
239 triggers |= AliAODForwardMult::kMCNSD;
240 fHTriggers->Fill(kMCNSD+0.5);
241 }
242
243 // Get the primary vertex from EG
244 TArrayF vtx;
245 genHeader->PrimaryVertex(vtx);
246 vz = vtx[2];
247
248 fHVertex->Fill(vz);
249 fHPhiR->Fill(phiR);
250 fHB->Fill(b);
251
252 // Check for the vertex bin
253 ivz = fHEventsTr->GetXaxis()->FindBin(vz);
254 if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
255 if (fDebug > 3) {
256 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
257 vz, fHEventsTr->GetXaxis()->GetXmin(),
258 fHEventsTr->GetXaxis()->GetXmax())); }
259 ivz = 0;
260 return kBadVertex;
261 }
262
263
264 return kOk;
265}
266//____________________________________________________________________
267Bool_t
268AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
269{
270 //
271 // Read centrality from event
272 //
273 // Parameters:
274 // esd Event
275 // cent On return, the centrality or negative if not found
276 //
277 // Return:
278 // False on error, true otherwise
279 //
280 AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
281 if (centObj) {
282 // AliInfo(Form("Got centrality object %p with quality %d",
283 // centObj, centObj->GetQuality()));
284 // centObj->Print();
285 if (centObj->GetQuality() == 0x8)
286 cent = centObj->GetCentralityPercentileUnchecked("V0M");
287 else
288 cent = centObj->GetCentralityPercentile("V0M");
289 }
290 // AliInfo(Form("Centrality is %f", cent));
291 fHCent->Fill(cent);
292
293 return true;
294}
295
296
297//
298// EOF
299//
300