Multiple fixes:
[u/mrichter/AliRoot.git] / PWGLF / 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 "AliStack.h"
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 "AliAnalysisManager.h"
33 #include "AliMCEventHandler.h"
34 #include "AliHeader.h"
35 #include <TList.h>
36 #include <TH2F.h>
37 #include <TParticle.h>
38 #include <TMath.h>
39
40 //====================================================================
41 AliFMDMCEventInspector::AliFMDMCEventInspector()
42   : AliFMDEventInspector(), 
43     fHVertex(0),
44     fHPhiR(0), 
45     fHB(0),
46     fHBvsPart(0),
47     fHBvsBin(0),
48     fHBvsCent(0),
49     fHVzComp(0),
50     fHCentVsPart(0),
51     fHCentVsBin(0),
52     fProduction("")
53 {
54   // 
55   // Constructor 
56   //
57 }
58
59 //____________________________________________________________________
60 AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
61   : AliFMDEventInspector("fmdEventInspector"), 
62     fHVertex(0),
63     fHPhiR(0), 
64     fHB(0),
65     fHBvsPart(0),
66     fHBvsBin(0),
67     fHBvsCent(0),
68     fHVzComp(0),
69     fHCentVsPart(0),
70     fHCentVsBin(0),
71     fProduction("")
72 {
73   // 
74   // Constructor 
75   // 
76   // Parameters:
77   //   name Name of object
78   //
79 }
80
81 //____________________________________________________________________
82 AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
83   : AliFMDEventInspector(o), 
84     fHVertex(0),
85     fHPhiR(0), 
86     fHB(0),
87     fHBvsPart(0),
88     fHBvsBin(0),
89     fHBvsCent(0),
90     fHVzComp(0),
91     fHCentVsPart(0),
92     fHCentVsBin(0),
93     fProduction("")
94 {
95   // 
96   // Copy constructor 
97   // 
98   // Parameters:
99   //   o Object to copy from 
100   //
101 }
102
103 //____________________________________________________________________
104 AliFMDMCEventInspector::~AliFMDMCEventInspector()
105 {
106   // 
107   // Destructor 
108   //
109 }
110 //____________________________________________________________________
111 AliFMDMCEventInspector&
112 AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
113 {
114   // 
115   // Assignement operator
116   // 
117   // Parameters:
118   //   o Object to assign from 
119   // 
120   // Return:
121   //    Reference to this object
122   //
123   AliFMDEventInspector::operator=(o);
124   return *this;
125 }
126
127 //____________________________________________________________________
128 void
129 AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
130 {
131   // 
132   // Initialize the object 
133   // 
134   // Parameters:
135   //   vtxAxis Vertex axis in use 
136   //
137   AliFMDEventInspector::Init(vtxAxis);
138
139   Int_t    maxPart = 450;
140   Int_t    maxBin  = 225;
141   Int_t    maxB    = 25;
142   Int_t    nVtx = vtxAxis.GetNbins();
143   Double_t lVtx = vtxAxis.GetXmin();
144   Double_t hVtx = vtxAxis.GetXmax();
145   fHVertex = new TH1F("vertex", "True vertex distribution", nVtx, lVtx, hVtx);
146   fHVertex->SetXTitle("v_{z} [cm]");
147   fHVertex->SetYTitle("# of events");
148   fHVertex->SetFillColor(kGreen+1);
149   fHVertex->SetFillStyle(3001);
150   fHVertex->SetDirectory(0);
151   // fHVertex->Sumw2();
152   fList->Add(fHVertex);
153
154   fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
155   fHPhiR->SetXTitle("#Phi_{R} [radians]");
156   fHPhiR->SetYTitle("# of events");
157   fHPhiR->SetFillColor(kGreen+1);
158   fHPhiR->SetFillStyle(3001);
159   fHPhiR->SetDirectory(0);
160   fList->Add(fHPhiR);
161
162   fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
163   fHB->SetXTitle("b [fm]");
164   fHB->SetYTitle("# of events");
165   fHB->SetFillColor(kGreen+1);
166   fHB->SetFillStyle(3001);
167   fHB->SetDirectory(0);
168   fList->Add(fHB);
169
170   fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
171                        5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
172   fHBvsPart->SetXTitle("b [fm]");
173   fHBvsPart->SetYTitle("# of participants");
174   fHBvsPart->SetZTitle("Events");
175   fHBvsPart->SetDirectory(0);
176   fList->Add(fHBvsPart);
177
178   fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
179                        5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
180   fHBvsBin->SetXTitle("b [fm]");
181   fHBvsBin->SetYTitle("# of binary collisions");
182   fHBvsBin->SetZTitle("Events");
183   fHBvsBin->SetDirectory(0);
184   fList->Add(fHBvsBin);
185   
186   fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
187                        5*maxB, 0, maxB, fCentAxis->GetNbins(), 
188                        fCentAxis->GetXbins()->GetArray());
189   fHBvsCent->SetXTitle("b [fm]");
190   fHBvsCent->SetYTitle("Centrality [%]");
191   fHBvsCent->SetZTitle("Event");
192   fHBvsCent->SetDirectory(0);
193   fList->Add(fHBvsCent);
194   
195   
196   fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
197                       10*nVtx, lVtx, hVtx, 10*nVtx, lVtx, hVtx);
198   fHVzComp->SetXTitle("True v_{z} [cm]");
199   fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
200   fHVzComp->SetZTitle("Events");
201   fHVzComp->SetDirectory(0);
202   fList->Add(fHVzComp);
203
204   fHCentVsPart = new TH2F("centralityVsParticipans", 
205                           "# of participants vs Centrality",
206                           maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(), 
207                           fCentAxis->GetXbins()->GetArray());
208   fHCentVsPart->SetXTitle("Participants");
209   fHCentVsPart->SetYTitle("Centrality [%]");
210   fHCentVsPart->SetZTitle("Event");
211   fHCentVsPart->SetDirectory(0);
212   fList->Add(fHCentVsPart);
213
214   fHCentVsBin = new TH2F("centralityVsBinary", 
215                          "# of binary collisions vs Centrality",
216                          maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(), 
217                          fCentAxis->GetXbins()->GetArray());
218   fHCentVsBin->SetXTitle("Binary collisions");
219   fHCentVsBin->SetYTitle("Centrality [%]");
220   fHCentVsBin->SetZTitle("Event");
221   fHCentVsBin->SetDirectory(0);
222   fList->Add(fHCentVsBin);
223 }
224
225 //____________________________________________________________________
226 void
227 AliFMDMCEventInspector::StoreInformation(Int_t runNo)
228 {
229   // Store information about running conditions in the output list 
230   if (!fList) return;
231   AliFMDEventInspector::StoreInformation(runNo);
232   TNamed* mc = new TNamed("mc", fProduction.Data());
233   mc->SetUniqueID(1);
234   fList->Add(mc);
235 }
236
237 namespace
238 {
239   TString& AppendField(TString& s, bool test, const char* f)
240   {
241     if (!test) return s;
242     if (!s.IsNull()) s.Append(",");
243     s.Append(f);
244     return s;
245   }
246 }
247
248 //____________________________________________________________________
249 void
250 AliFMDMCEventInspector::ReadProductionDetails(AliMCEvent* event)
251 {
252   //Assign MC only triggers : True NSD etc.
253   AliHeader*               header          = event->Header();
254   AliGenEventHeader*       genHeader       = header->GenEventHeader();
255   AliCollisionGeometry*    colGeometry     = 
256     dynamic_cast<AliCollisionGeometry*>(genHeader);
257   AliGenPythiaEventHeader* pythiaHeader    = 
258     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
259   AliGenDPMjetEventHeader* dpmHeader       = 
260     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
261   AliGenGeVSimEventHeader* gevHeader       = 
262     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
263   AliGenHerwigEventHeader* herwigHeader    = 
264     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
265   AliGenHijingEventHeader* hijingHeader    = 
266     dynamic_cast<AliGenHijingEventHeader*>(genHeader);
267   // AliGenHydjetEventHeader* hydjetHeader    = 
268   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
269   AliGenEposEventHeader*   eposHeader      = 
270     dynamic_cast<AliGenEposEventHeader*>(genHeader);
271
272   AppendField(fProduction, colGeometry,  "Geometry");
273   AppendField(fProduction, pythiaHeader, "Pythia");
274   AppendField(fProduction, dpmHeader,    "DPM");
275   AppendField(fProduction, gevHeader,    "GevSim");
276   AppendField(fProduction, herwigHeader, "Herwig");
277   AppendField(fProduction, hijingHeader, "Hijing");
278   // AppendField(fProduction, hydjetHeader, "Hydjet");
279   AppendField(fProduction, eposHeader,   "EPOS");
280 }
281
282 //____________________________________________________________________
283 UInt_t
284 AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event, 
285                                   UInt_t&           triggers,
286                                   UShort_t&         ivz, 
287                                   Double_t&         vz,
288                                   Double_t&         b,
289                                   Double_t&         c,
290                                   Int_t&            npart, 
291                                   Int_t&            nbin,
292                                   Double_t&         phiR)
293 {
294   // 
295   // Process the event 
296   // 
297   // Parameters:
298   //   event     Input event 
299   //   triggers  On return, the triggers fired 
300   //   ivz       On return, the found vertex bin (1-based).  A zero
301   //                  means outside of the defined vertex range
302   //   vz        On return, the z position of the interaction
303   //   b         On return, the impact parameter - < 0 if not available
304   //   c         On return, centrality estimate - < 0 if not available 
305   //   phiR      On return, the event plane angle - < 0 if not available 
306   // 
307   // Return:
308   //    0 (or kOk) on success, otherwise a bit mask of error codes 
309   //
310
311   // Check that we have an event 
312   if (!event) { 
313     AliWarning("No MC event found for input event");
314     return kNoEvent;
315   }
316
317   //Assign MC only triggers : True NSD etc.
318   AliHeader*               header          = event->Header();
319   AliGenEventHeader*       genHeader       = header->GenEventHeader();
320   AliCollisionGeometry*    colGeometry     = 
321     dynamic_cast<AliCollisionGeometry*>(genHeader);
322   AliGenPythiaEventHeader* pythiaHeader    = 
323     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
324   AliGenDPMjetEventHeader* dpmHeader       = 
325     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
326   AliGenGeVSimEventHeader* gevHeader       = 
327     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
328   AliGenHerwigEventHeader* herwigHeader    = 
329     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
330   // AliGenHijingEventHeader* hijingHeader    = 
331   //   dynamic_cast<AliGenHijingEventHeader*>(genHeader);
332   // AliGenHydjetEventHeader* hydjetHeader    = 
333   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
334   // AliGenEposEventHeader*   eposHeader      = 
335   //   dynamic_cast<AliGenEposEventHeader*>(genHeader);
336   
337   // Check if this is a single diffractive event 
338   Bool_t   sd    = kFALSE;
339   Double_t phi   = -1111;
340   npart          = 0;
341   nbin           = 0;
342   b              = -1;
343   c              = -1;
344   if (colGeometry) { 
345     b     = colGeometry->ImpactParameter();
346     phi   = colGeometry->ReactionPlaneAngle();
347     npart = (colGeometry->ProjectileParticipants() + 
348              colGeometry->TargetParticipants());
349     nbin  = colGeometry->NN();
350   }
351   if(pythiaHeader) {
352     Int_t pythiaType = pythiaHeader->ProcessType();
353     // 92 and 93 are SD 
354     // 94 is DD 
355     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
356     b     = pythiaHeader->GetImpactParameter();
357     npart = 2; // Always 2 protons
358     nbin  = 1; // Always 1 binary collision 
359   }
360   if (b >= 0) { 
361     if      (b <  3.5)  c = 2.5; //0-5%
362     else if (b <  4.95) c = 7.5; //5-10%
363     else if (b <  6.98) c = 15; //10-20%
364     else if (b <  8.55) c = 25; //20-30%
365     else if (b <  9.88) c = 35; //30-40%
366     else if (b < 11.04) c = 45; //40-50%
367     else                c = 55; //50-60%
368   }
369   if(dpmHeader) { // Also an AliCollisionGeometry 
370     Int_t processType = dpmHeader->ProcessType();
371     // 1 & 4 are ND 
372     // 5 & 6 are SD 
373     // 7 is DD 
374     if (processType == 5 || processType == 6)  sd = kTRUE;
375   }
376   if (gevHeader) { 
377     phi  = gevHeader->GetEventPlane();
378   }
379   if (herwigHeader) {
380     Int_t processType = herwigHeader->ProcessType();
381     // This is a guess 
382     if (processType == 5 || processType == 6)  sd = kTRUE;
383     npart = 2; // Always 2 protons
384     nbin  = 1; // Always 1 binary collision 
385   }
386   // if (hijingHeader) { 
387   // b    = hijingHeader->ImpactParameter();
388   // phi  = hijingHeader->ReactionPlaneAngle();
389   // }
390   // if (hydjetHeader) {
391   //   b    = hydjetHeader->ImpactParameter();
392   //   phi  = hydjetHeader->ReactionPlaneAngle();
393   // }
394   // if (eposHeader) {
395   //   b    = eposHeader->ImpactParameter();
396   // phi  = eposHeader->ReactionPlaneAngle();
397   // }
398
399   // Normalize event plane angle to [0,2pi]
400   if (phi <= -1111) phiR = -1;
401   else { 
402     while (true) {
403       if      (phi < 0)             phi += 2*TMath::Pi();
404       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
405       else                          break;
406     }
407   }
408
409   // Do a check on particles
410   sd = IsSingleDiffractive(event->Stack());
411
412   // Set NSD flag
413   if(!sd) {
414     triggers |= AliAODForwardMult::kMCNSD;
415     fHTriggers->Fill(kMCNSD+0.5);
416   }
417   
418   // Get the primary vertex from EG 
419   TArrayF vtx;
420   genHeader->PrimaryVertex(vtx);
421   vz = vtx[2];
422
423   fHVertex->Fill(vz);
424   fHPhiR->Fill(phiR);
425   fHB->Fill(b);
426   fHBvsPart->Fill(b, npart);
427   fHBvsBin->Fill(b, nbin);
428
429   if(fUseDisplacedVertices) {
430     // Put the vertex at fixed locations 
431     Double_t zvtx  = vz;
432     Double_t ratio = zvtx/37.5;
433     if(ratio > 0) ratio = ratio + 0.5;
434     if(ratio < 0) ratio = ratio - 0.5;
435     Int_t ratioInt = Int_t(ratio);
436     zvtx = 37.5*((Double_t)ratioInt);
437     if(TMath::Abs(zvtx) > 999) 
438       return kBadVertex;
439   }
440
441   // Check for the vertex bin 
442   ivz = fVtxAxis.FindBin(vz);
443   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
444     if (fDebug > 3) {
445       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
446                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
447     ivz = 0;
448     return kBadVertex;
449   }
450
451   
452   return kOk;
453 }
454
455 //____________________________________________________________________
456 namespace {
457   Double_t rapidity(TParticle* p, Double_t mass)
458   {
459     Double_t pt = p->Pt();
460     Double_t pz = p->Pz();
461     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
462     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
463     return .5 * TMath::Log((ee + pz) / (ee - pz)); 
464   }
465 }
466
467 //____________________________________________________________________
468 Bool_t
469 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
470                                             Double_t xiMin, 
471                                             Double_t xiMax) const
472 {
473   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
474   // 
475   // This is re-implemented here to be indendent of the PWG0 library. 
476   TParticle* p1 = 0;     // Particle with least y 
477   TParticle* p2 = 0;     // Particle with largest y 
478   Double_t   y1 =  1e10; // y of p1 
479   Double_t   y2 = -1e10; // y of p2 
480   
481   // Loop over primaries 
482   for (Int_t i = 0; i < stack->GetNprimary(); i++) { 
483     TParticle* p = stack->Particle(i);
484     if (!p) continue;
485     
486     Int_t pdg = TMath::Abs(p->GetPdgCode());
487     Int_t c1  = p->GetFirstDaughter();
488     Int_t s   = p->GetStatusCode();
489     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
490
491     // Select final state charm and beauty 
492     if (c1 > -1 || s != 1) mfl = 0; 
493
494     // Check if this is a primary, pi0, Sigma0, ???, or most
495     // significant digit is larger than or equal to 4
496     if (!(stack->IsPhysicalPrimary(i) || 
497           pdg == 111  || 
498           pdg == 3212 || 
499           pdg == 3124 || 
500           mfl >= 4)) continue;
501
502     Int_t m1 = p->GetFirstMother();
503     if (m1 > 0) { 
504       TParticle* pm1  = stack->Particle(m1);
505       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
506       // Check if mother is a p0, Simga0, or ???
507       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
508     }
509     
510     // Calculate the rapidity of the particle 
511     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
512     Double_t yy = rapidity(p, mm);
513     
514     // Check if the rapidity of this particle is further out than any
515     // of the preceding particles
516     if (yy < y1) { 
517       y1 = yy;
518       p1 = p;
519     }
520     if (yy > y2) { 
521       y2 = yy;
522       p2 = p;
523     }
524   }
525   if (!p1 || !p2) return false;
526   
527   // Calculate rapidities assuming protons 
528   y1            = TMath::Abs(rapidity(p1, 0.938));
529   y2            = TMath::Abs(rapidity(p2, 0.938));
530
531   // Check if both or just one is a proton
532   Int_t    pdg1 = p1->GetPdgCode();
533   Int_t    pdg2 = p2->GetPdgCode();
534   Int_t    arm  = -99999;
535   if (pdg1 == 2212 && pdg2 == 2212) 
536     arm = (y1 > y2 ? 0 : 1);
537   else if (pdg1 == 2212) 
538     arm = 0;
539   else if (pdg2 == 2212) 
540     arm = 1;
541   else 
542     return false;
543   
544   // Rapidity shift
545   Double_t m02s = 1 - 2 * p1->Energy() / fEnergy; 
546   Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
547   
548   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
549   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
550     
551   return false;
552 }
553
554 //____________________________________________________________________
555 Bool_t
556 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz, 
557                                        Double_t cent,  Double_t b,
558                                        Int_t    npart, Int_t    nbin)
559 {
560   fHVzComp->Fill(trueVz, vz);
561   fHBvsCent->Fill(b, cent);
562   fHCentVsPart->Fill(npart, cent);
563   fHCentVsBin->Fill(nbin, cent);
564
565   return true;
566 }  
567 //
568 // EOF
569 //
570