Fixes for pA indenfication of events
[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 (fDebug && !colGeometry) { 
352     AliWarningF("Collision header of class %s is not a CollisionHeader", 
353                 genHeader->ClassName());
354   }
355     
356   if(pythiaHeader) {
357     Int_t pythiaType = pythiaHeader->ProcessType();
358     // 92 and 93 are SD 
359     // 94 is DD 
360     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
361     b     = pythiaHeader->GetImpactParameter();
362     npart = 2; // Always 2 protons
363     nbin  = 1; // Always 1 binary collision 
364   }
365   if (b >= 0) { 
366     if      (b <  3.5)  c = 2.5; //0-5%
367     else if (b <  4.95) c = 7.5; //5-10%
368     else if (b <  6.98) c = 15; //10-20%
369     else if (b <  8.55) c = 25; //20-30%
370     else if (b <  9.88) c = 35; //30-40%
371     else if (b < 11.04) c = 45; //40-50%
372     else                c = 55; //50-60%
373   }
374   if(dpmHeader) { // Also an AliCollisionGeometry 
375     Int_t processType = dpmHeader->ProcessType();
376     // 1 & 4 are ND 
377     // 5 & 6 are SD 
378     // 7 is DD 
379     if (processType == 5 || processType == 6)  sd = kTRUE;
380   }
381   if (gevHeader) { 
382     phi  = gevHeader->GetEventPlane();
383   }
384   if (herwigHeader) {
385     Int_t processType = herwigHeader->ProcessType();
386     // This is a guess 
387     if (processType == 5 || processType == 6)  sd = kTRUE;
388     npart = 2; // Always 2 protons
389     nbin  = 1; // Always 1 binary collision 
390   }
391   // if (hijingHeader) { 
392   // b    = hijingHeader->ImpactParameter();
393   // phi  = hijingHeader->ReactionPlaneAngle();
394   // }
395   // if (hydjetHeader) {
396   //   b    = hydjetHeader->ImpactParameter();
397   //   phi  = hydjetHeader->ReactionPlaneAngle();
398   // }
399   // if (eposHeader) {
400   //   b    = eposHeader->ImpactParameter();
401   // phi  = eposHeader->ReactionPlaneAngle();
402   // }
403
404   // Normalize event plane angle to [0,2pi]
405   if (phi <= -1111) phiR = -1;
406   else { 
407     while (true) {
408       if      (phi < 0)             phi += 2*TMath::Pi();
409       else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
410       else                          break;
411     }
412   }
413
414   // Do a check on particles
415   sd = IsSingleDiffractive(event->Stack());
416
417   // Set NSD flag
418   if(!sd) {
419     triggers |= AliAODForwardMult::kMCNSD;
420     fHTriggers->Fill(kMCNSD+0.5);
421   }
422   
423   // Get the primary vertex from EG 
424   TArrayF vtx;
425   genHeader->PrimaryVertex(vtx);
426   vz = vtx[2];
427
428   if (fDebug) { 
429     AliInfoF("vz=%f, phiR=%f, b=%f, npart=%d, nbin=%d", 
430              vz, phiR, b, npart, nbin);
431   }
432   fHVertex->Fill(vz);
433   fHPhiR->Fill(phiR);
434   fHB->Fill(b);
435   fHBvsPart->Fill(b, npart);
436   fHBvsBin->Fill(b, nbin);
437
438   if(fUseDisplacedVertices) {
439     // Put the vertex at fixed locations 
440     Double_t zvtx  = vz;
441     Double_t ratio = zvtx/37.5;
442     if(ratio > 0) ratio = ratio + 0.5;
443     if(ratio < 0) ratio = ratio - 0.5;
444     Int_t ratioInt = Int_t(ratio);
445     zvtx = 37.5*((Double_t)ratioInt);
446     if(TMath::Abs(zvtx) > 999) 
447       return kBadVertex;
448   }
449
450   // Check for the vertex bin 
451   ivz = fVtxAxis.FindBin(vz);
452   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
453     if (fDebug > 3) {
454       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
455                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
456     ivz = 0;
457     return kBadVertex;
458   }
459
460   
461   return kOk;
462 }
463
464 //____________________________________________________________________
465 namespace {
466   Double_t rapidity(TParticle* p, Double_t mass)
467   {
468     Double_t pt = p->Pt();
469     Double_t pz = p->Pz();
470     Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
471     if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
472     return .5 * TMath::Log((ee + pz) / (ee - pz)); 
473   }
474 }
475
476 //____________________________________________________________________
477 Bool_t
478 AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
479                                             Double_t xiMin, 
480                                             Double_t xiMax) const
481 {
482   // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
483   // 
484   // This is re-implemented here to be indendent of the PWG0 library. 
485   TParticle* p1 = 0;     // Particle with least y 
486   TParticle* p2 = 0;     // Particle with largest y 
487   Double_t   y1 =  1e10; // y of p1 
488   Double_t   y2 = -1e10; // y of p2 
489   
490   // Loop over primaries 
491   for (Int_t i = 0; i < stack->GetNprimary(); i++) { 
492     TParticle* p = stack->Particle(i);
493     if (!p) continue;
494     
495     Int_t pdg = TMath::Abs(p->GetPdgCode());
496     Int_t c1  = p->GetFirstDaughter();
497     Int_t s   = p->GetStatusCode();
498     Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
499
500     // Select final state charm and beauty 
501     if (c1 > -1 || s != 1) mfl = 0; 
502
503     // Check if this is a primary, pi0, Sigma0, ???, or most
504     // significant digit is larger than or equal to 4
505     if (!(stack->IsPhysicalPrimary(i) || 
506           pdg == 111  || 
507           pdg == 3212 || 
508           pdg == 3124 || 
509           mfl >= 4)) continue;
510
511     Int_t m1 = p->GetFirstMother();
512     if (m1 > 0) { 
513       TParticle* pm1  = stack->Particle(m1);
514       Int_t      mpdg = TMath::Abs(pm1->GetPdgCode());
515       // Check if mother is a p0, Simga0, or ???
516       if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
517     }
518     
519     // Calculate the rapidity of the particle 
520     Double_t mm = (pdg != 3124 ? p->GetMass() : 1.5195);
521     Double_t yy = rapidity(p, mm);
522     
523     // Check if the rapidity of this particle is further out than any
524     // of the preceding particles
525     if (yy < y1) { 
526       y1 = yy;
527       p1 = p;
528     }
529     if (yy > y2) { 
530       y2 = yy;
531       p2 = p;
532     }
533   }
534   if (!p1 || !p2) return false;
535   
536   // Calculate rapidities assuming protons 
537   y1            = TMath::Abs(rapidity(p1, 0.938));
538   y2            = TMath::Abs(rapidity(p2, 0.938));
539
540   // Check if both or just one is a proton
541   Int_t    pdg1 = p1->GetPdgCode();
542   Int_t    pdg2 = p2->GetPdgCode();
543   Int_t    arm  = -99999;
544   if (pdg1 == 2212 && pdg2 == 2212) 
545     arm = (y1 > y2 ? 0 : 1);
546   else if (pdg1 == 2212) 
547     arm = 0;
548   else if (pdg2 == 2212) 
549     arm = 1;
550   else 
551     return false;
552   
553   // Rapidity shift
554   Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0); 
555   Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
556   
557   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
558   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
559     
560   return false;
561 }
562
563 //____________________________________________________________________
564 Bool_t
565 AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz, 
566                                        Double_t cent,  Double_t b,
567                                        Int_t    npart, Int_t    nbin)
568 {
569   fHVzComp->Fill(trueVz, vz);
570   fHBvsCent->Fill(b, cent);
571   fHCentVsPart->Fill(npart, cent);
572   fHCentVsBin->Fill(nbin, cent);
573
574   return true;
575 }  
576 //
577 // EOF
578 //
579