MC upgrades to displaced vertices
[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                                   Int_t&            npart, 
290                                   Int_t&            nbin,
291                                   Double_t&         phiR)
292 {
293   // 
294   // Process the event 
295   // 
296   // Parameters:
297   //   event     Input event 
298   //   triggers  On return, the triggers fired 
299   //   ivz       On return, the found vertex bin (1-based).  A zero
300   //                  means outside of the defined vertex range
301   //   vz        On return, the z position of the interaction
302   //   b         On return, the impact parameter - < 0 if not available
303   //   phiR      On return, the event plane angle - < 0 if not available 
304   // 
305   // Return:
306   //    0 (or kOk) on success, otherwise a bit mask of error codes 
307   //
308
309   // Check that we have an event 
310   if (!event) { 
311     AliWarning("No MC event found for input event");
312     return kNoEvent;
313   }
314
315   //Assign MC only triggers : True NSD etc.
316   AliHeader*               header          = event->Header();
317   AliGenEventHeader*       genHeader       = header->GenEventHeader();
318   AliCollisionGeometry*    colGeometry     = 
319     dynamic_cast<AliCollisionGeometry*>(genHeader);
320   AliGenPythiaEventHeader* pythiaHeader    = 
321     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
322   AliGenDPMjetEventHeader* dpmHeader       = 
323     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
324   AliGenGeVSimEventHeader* gevHeader       = 
325     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
326   AliGenHerwigEventHeader* herwigHeader    = 
327     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
328   // AliGenHijingEventHeader* hijingHeader    = 
329   //   dynamic_cast<AliGenHijingEventHeader*>(genHeader);
330   // AliGenHydjetEventHeader* hydjetHeader    = 
331   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
332   // AliGenEposEventHeader*   eposHeader      = 
333   //   dynamic_cast<AliGenEposEventHeader*>(genHeader);
334   
335   // Check if this is a single diffractive event 
336   Bool_t   sd    = kFALSE;
337   Double_t phi   = -1111;
338   npart          = 0;
339   nbin           = 0;
340   b              = -1;
341   if (colGeometry) { 
342     b     = colGeometry->ImpactParameter();
343     phi   = colGeometry->ReactionPlaneAngle();
344     npart = (colGeometry->ProjectileParticipants() + 
345              colGeometry->TargetParticipants());
346     nbin  = colGeometry->NN();
347   }
348   if(pythiaHeader) {
349     Int_t pythiaType = pythiaHeader->ProcessType();
350     // 92 and 93 are SD 
351     // 94 is DD 
352     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
353     b     = pythiaHeader->GetImpactParameter();
354     npart = 2; // Always 2 protons
355     nbin  = 1; // Always 1 binary collision 
356   }
357   if(dpmHeader) { // Also an AliCollisionGeometry 
358     Int_t processType = dpmHeader->ProcessType();
359     // 1 & 4 are ND 
360     // 5 & 6 are SD 
361     // 7 is DD 
362     if (processType == 5 || processType == 6)  sd = kTRUE;
363   }
364   if (gevHeader) { 
365     phi  = gevHeader->GetEventPlane();
366   }
367   if (herwigHeader) {
368     Int_t processType = herwigHeader->ProcessType();
369     // This is a guess 
370     if (processType == 5 || processType == 6)  sd = kTRUE;
371     npart = 2; // Always 2 protons
372     nbin  = 1; // Always 1 binary collision 
373   }
374   // if (hijingHeader) { 
375   // b    = hijingHeader->ImpactParameter();
376   // phi  = hijingHeader->ReactionPlaneAngle();
377   // }
378   // if (hydjetHeader) {
379   //   b    = hydjetHeader->ImpactParameter();
380   //   phi  = hydjetHeader->ReactionPlaneAngle();
381   // }
382   // if (eposHeader) {
383   //   b    = eposHeader->ImpactParameter();
384   // phi  = eposHeader->ReactionPlaneAngle();
385   // }
386
387   // Normalize event plane angle to [0,2pi]
388   if (phi <= -1111) phiR = -1;
389   else { 
390     while (true) {
391       if      (phi < 0)             phi += 2*TMath::Pi();
392       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
393       else                          break;
394     }
395   }
396
397   // Do a check on particles
398   sd = IsSingleDiffractive(event->Stack());
399
400   // Set NSD flag
401   if(!sd) {
402     triggers |= AliAODForwardMult::kMCNSD;
403     fHTriggers->Fill(kMCNSD+0.5);
404   }
405   
406   // Get the primary vertex from EG 
407   TArrayF vtx;
408   genHeader->PrimaryVertex(vtx);
409   vz = vtx[2];
410
411   fHVertex->Fill(vz);
412   fHPhiR->Fill(phiR);
413   fHB->Fill(b);
414   fHBvsPart->Fill(b, npart);
415   fHBvsBin->Fill(b, nbin);
416
417   // Check for the vertex bin 
418   ivz = fVtxAxis.FindBin(vz);
419   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
420     if (fDebug > 3) {
421       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
422                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
423     ivz = 0;
424     return kBadVertex;
425   }
426
427   
428   return kOk;
429 }
430
431 //____________________________________________________________________
432 namespace {
433   Double_t rapidity(TParticle* p, Double_t mass)
434   {
435     Double_t pt = p->Pt();
436     Double_t pz = p->Pz();
437     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
438     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
439     return .5 * TMath::Log((ee + pz) / (ee - pz)); 
440   }
441 }
442
443 //____________________________________________________________________
444 Bool_t
445 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
446                                             Double_t xiMin, 
447                                             Double_t xiMax) const
448 {
449   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
450   // 
451   // This is re-implemented here to be indendent of the PWG0 library. 
452   TParticle* p1 = 0;     // Particle with least y 
453   TParticle* p2 = 0;     // Particle with largest y 
454   Double_t   y1 =  1e10; // y of p1 
455   Double_t   y2 = -1e10; // y of p2 
456   
457   // Loop over primaries 
458   for (Int_t i = 0; i < stack->GetNprimary(); i++) { 
459     TParticle* p = stack->Particle(i);
460     if (!p) continue;
461     
462     Int_t pdg = TMath::Abs(p->GetPdgCode());
463     Int_t c1  = p->GetFirstDaughter();
464     Int_t s   = p->GetStatusCode();
465     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
466
467     // Select final state charm and beauty 
468     if (c1 > -1 || s != 1) mfl = 0; 
469
470     // Check if this is a primary, pi0, Sigma0, ???, or most
471     // significant digit is larger than or equal to 4
472     if (!(stack->IsPhysicalPrimary(i) || 
473           pdg == 111  || 
474           pdg == 3212 || 
475           pdg == 3124 || 
476           mfl >= 4)) continue;
477
478     Int_t m1 = p->GetFirstMother();
479     if (m1 > 0) { 
480       TParticle* pm1  = stack->Particle(m1);
481       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
482       // Check if mother is a p0, Simga0, or ???
483       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
484     }
485     
486     // Calculate the rapidity of the particle 
487     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
488     Double_t yy = rapidity(p, mm);
489     
490     // Check if the rapidity of this particle is further out than any
491     // of the preceding particles
492     if (yy < y1) { 
493       y1 = yy;
494       p1 = p;
495     }
496     if (yy > y2) { 
497       y2 = yy;
498       p2 = p;
499     }
500   }
501   if (!p1 || !p2) return false;
502   
503   // Calculate rapidities assuming protons 
504   y1            = TMath::Abs(rapidity(p1, 0.938));
505   y2            = TMath::Abs(rapidity(p2, 0.938));
506
507   // Check if both or just one is a proton
508   Int_t    pdg1 = p1->GetPdgCode();
509   Int_t    pdg2 = p2->GetPdgCode();
510   Int_t    arm  = -99999;
511   if (pdg1 == 2212 && pdg2 == 2212) 
512     arm = (y1 > y2 ? 0 : 1);
513   else if (pdg1 == 2212) 
514     arm = 0;
515   else if (pdg2 == 2212) 
516     arm = 1;
517   else 
518     return false;
519   
520   // Rapidity shift
521   Double_t m02s = 1 - 2 * p1->Energy() / fEnergy; 
522   Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
523   
524   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
525   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
526     
527   return false;
528 }
529 //____________________________________________________________________
530 Bool_t
531 AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, 
532                                        Double_t& cent, 
533                                        UShort_t& qual) const
534 {
535   // 
536   // Read centrality from event 
537   // 
538   // Parameters:
539   //    esd  Event 
540   //    cent On return, the centrality or negative if not found
541   // 
542   // Return:
543   //    False on error, true otherwise 
544   //
545   cent = -1;
546   qual = 0;
547   AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
548   if (!centObj) return true;
549 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
550   Bool_t isMC = am->GetMCtruthEventHandler() != 0;
551   
552   //std::cout<<fUseDisplacedVertices<<"  "<<isMC<<std::endl;
553   
554   if(fUseDisplacedVertices && isMC) {
555     
556     
557     AliMCEventHandler* mchandler = static_cast<AliMCEventHandler*>(am->GetMCtruthEventHandler());
558     AliMCEvent* mcevent          = mchandler->MCEvent();
559     
560     AliHeader*               header          = mcevent->Header();
561     AliGenEventHeader*       genHeader       = header->GenEventHeader();
562     AliCollisionGeometry*    colGeometry     = 
563       dynamic_cast<AliCollisionGeometry*>(genHeader);
564     Double_t b              = -1;
565     if (colGeometry)  
566       b     = colGeometry->ImpactParameter();
567     //std::cout<<"Hallo!!  "<<b<<std::endl;
568     cent = -1;
569     if(b<3.5 && b >0) cent = 2.5; //0-5%
570     if(b>3.5 && b<4.95) cent = 7.5; //5-10%
571     if(b>4.95 && b<6.98) cent = 15; //10-20%
572     if(b>6.98 && b<8.55) cent = 25; //20-30%
573     if(b>8.55 && b<9.88) cent = 35; //30-40%
574     if(b>9.88 && b<11.04) cent = 45; //40-50%
575     if(b>11.04) cent = 55; //50-60%
576     //cent = 10;
577     qual = 0;
578   }
579   else {
580     
581     qual = centObj->GetQuality();
582     if (qual == 0x8) // Ignore ZDC outliers 
583       cent = centObj->GetCentralityPercentileUnchecked("V0M");  
584     else
585       cent = centObj->GetCentralityPercentile("V0M");        
586   }
587   return true;
588 }
589
590 //____________________________________________________________________
591 Bool_t
592 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz, 
593                                        Double_t cent,  Double_t b,
594                                        Int_t    npart, Int_t    nbin)
595 {
596   fHVzComp->Fill(trueVz, vz);
597   fHBvsCent->Fill(b, cent);
598   fHCentVsPart->Fill(npart, cent);
599   fHCentVsBin->Fill(nbin, cent);
600
601   return true;
602 }  
603 //
604 // EOF
605 //
606