74a9b0119b843f1a7543c4fe67e4804fc6f7d44c
[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 //====================================================================
34 AliFMDMCEventInspector::AliFMDMCEventInspector()
35   : AliFMDEventInspector(), 
36     fHVertex(0),
37     fHPhiR(0), 
38     fHB(0)
39 {
40   // 
41   // Constructor 
42   //
43 }
44
45 //____________________________________________________________________
46 AliFMDMCEventInspector::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 //____________________________________________________________________
61 AliFMDMCEventInspector::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 //____________________________________________________________________
76 AliFMDMCEventInspector::~AliFMDMCEventInspector()
77 {
78   // 
79   // Destructor 
80   //
81 }
82 //____________________________________________________________________
83 AliFMDMCEventInspector&
84 AliFMDMCEventInspector::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 //____________________________________________________________________
100 void
101 AliFMDMCEventInspector::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 //____________________________________________________________________
142 UInt_t
143 AliFMDMCEventInspector::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 //____________________________________________________________________
267 Bool_t
268 AliFMDMCEventInspector::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