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