Various fixes to deal with centrality
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCEventInspector.cxx
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 #include <TH2F.h>
34
35 //====================================================================
36 AliFMDMCEventInspector::AliFMDMCEventInspector()
37   : AliFMDEventInspector(), 
38     fHVertex(0),
39     fHPhiR(0), 
40     fHB(0),
41     fHBvsPart(0),
42     fHBvsBin(0),
43     fHBvsCent(0),
44     fHVzComp(0),
45     fHCentVsPart(0),
46     fHCentVsBin(0)
47 {
48   // 
49   // Constructor 
50   //
51 }
52
53 //____________________________________________________________________
54 AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
55   : AliFMDEventInspector("fmdEventInspector"), 
56     fHVertex(0),
57     fHPhiR(0), 
58     fHB(0),
59     fHBvsPart(0),
60     fHBvsBin(0),
61     fHBvsCent(0),
62     fHVzComp(0),
63     fHCentVsPart(0),
64     fHCentVsBin(0)
65 {
66   // 
67   // Constructor 
68   // 
69   // Parameters:
70   //   name Name of object
71   //
72 }
73
74 //____________________________________________________________________
75 AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
76   : AliFMDEventInspector(o), 
77     fHVertex(0),
78     fHPhiR(0), 
79     fHB(0),
80     fHBvsPart(0),
81     fHBvsBin(0),
82     fHBvsCent(0),
83     fHVzComp(0),
84     fHCentVsPart(0),
85     fHCentVsBin(0)
86 {
87   // 
88   // Copy constructor 
89   // 
90   // Parameters:
91   //   o Object to copy from 
92   //
93 }
94
95 //____________________________________________________________________
96 AliFMDMCEventInspector::~AliFMDMCEventInspector()
97 {
98   // 
99   // Destructor 
100   //
101 }
102 //____________________________________________________________________
103 AliFMDMCEventInspector&
104 AliFMDMCEventInspector::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 //____________________________________________________________________
120 void
121 AliFMDMCEventInspector::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
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);
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
154   fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
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);
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);
177   
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);
215 }
216
217 //____________________________________________________________________
218 UInt_t
219 AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event, 
220                                   UInt_t&           triggers,
221                                   UShort_t&         ivz, 
222                                   Double_t&         vz,
223                                   Double_t&         b,
224                                   Int_t&            npart, 
225                                   Int_t&            nbin,
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();
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   // Check if this is a single diffractive event 
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   }
283   if(pythiaHeader) {
284     Int_t pythiaType = pythiaHeader->ProcessType();
285     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
286     b     = pythiaHeader->GetImpactParameter();
287     npart = 2; // Always 2 protons
288     nbin  = 1; // Always 1 binary collision 
289   }
290   if(dpmHeader) { // Also an AliCollisionGeometry 
291     Int_t processType = dpmHeader->ProcessType();
292     if (processType == 5 || processType == 6)  sd = kTRUE;
293   }
294   if (gevHeader) { 
295     phi  = gevHeader->GetEventPlane();
296   }
297   if (herwigHeader) {
298     Int_t processType = herwigHeader->ProcessType();
299     // This is a guess 
300     if (processType == 5 || processType == 6)  sd = kTRUE;
301     npart = 2; // Always 2 protons
302     nbin  = 1; // Always 1 binary collision 
303   }
304   // if (hijingHeader) { 
305   // b    = hijingHeader->ImpactParameter();
306   // phi  = hijingHeader->ReactionPlaneAngle();
307   // }
308   // if (hydjetHeader) {
309   //   b    = hydjetHeader->ImpactParameter();
310   //   phi  = hydjetHeader->ReactionPlaneAngle();
311   // }
312   // if (eposHeader) {
313   //   b    = eposHeader->ImpactParameter();
314   // phi  = eposHeader->ReactionPlaneAngle();
315   // }
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);
341   fHBvsPart->Fill(b, npart);
342   fHBvsBin->Fill(b, nbin);
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 //____________________________________________________________________
359 Bool_t
360 AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, 
361                                        Double_t& cent, 
362                                        UShort_t& qual) const
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   //
374   cent = -1;
375   qual = 0;
376   AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
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");        
384
385   return true;
386 }
387
388 //____________________________________________________________________
389 Bool_t
390 AliFMDMCEventInspector::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 }  
401 //
402 // EOF
403 //
404