AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFMDMCEventInspector.cxx
Go to the documentation of this file.
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 "AliCollisionGeometry.h"
19 #include "AliLog.h"
20 #include "AliAODForwardMult.h"
21 #include "AliForwardUtil.h"
22 #include "AliCentrality.h"
23 #include "AliESDEvent.h"
24 #include "AliMCEvent.h"
25 #include "AliStack.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliGenDPMjetEventHeader.h"
28 #include "AliGenHijingEventHeader.h"
29 // #include "AliGenHydjetEventHeader.h"
30 #include "AliGenEposEventHeader.h"
31 #include "AliGenHerwigEventHeader.h"
32 #include "AliGenGeVSimEventHeader.h"
33 #include "AliAnalysisManager.h"
34 #include "AliMCEventHandler.h"
35 #include "AliHeader.h"
36 #include <TList.h>
37 #include <TH2F.h>
38 #include <TParticle.h>
39 #include <TMath.h>
40 #include <TParameter.h>
41 
42 //====================================================================
45  fHVertex(0),
46  fHVertexXY(0),
47  fHPhiR(0),
48  fHB(0),
49  fHMcC(0),
50  fHBvsPart(0),
51  fHBvsBin(0),
52  fHBvsCent(0),
53  fHVzComp(0),
54  fHCentVsPart(0),
55  fHCentVsBin(0),
56  fHCentVsMcC(0),
57  fProduction("")
58 {
59  //
60  // Constructor
61  //
62 }
63 
64 //____________________________________________________________________
66  : AliFMDEventInspector("fmdEventInspector"),
67  fHVertex(0),
68  fHVertexXY(0),
69  fHPhiR(0),
70  fHB(0),
71  fHMcC(0),
72  fHBvsPart(0),
73  fHBvsBin(0),
74  fHBvsCent(0),
75  fHVzComp(0),
76  fHCentVsPart(0),
77  fHCentVsBin(0),
78  fHCentVsMcC(0),
79  fProduction("")
80 {
81  //
82  // Constructor
83  //
84  // Parameters:
85  // name Name of object
86  //
87  fMC = true;
88 }
89 
90 //____________________________________________________________________
93  fHVertex(0),
94  fHVertexXY(0),
95  fHPhiR(0),
96  fHB(0),
97  fHMcC(0),
98  fHBvsPart(0),
99  fHBvsBin(0),
100  fHBvsCent(0),
101  fHVzComp(0),
102  fHCentVsPart(0),
103  fHCentVsBin(0),
104  fHCentVsMcC(0),
105  fProduction("")
106 {
107  //
108  // Copy constructor
109  //
110  // Parameters:
111  // o Object to copy from
112  //
113 }
114 
115 //____________________________________________________________________
117 {
118  //
119  // Destructor
120  //
121 }
122 //____________________________________________________________________
125 {
126  //
127  // Assignement operator
128  //
129  // Parameters:
130  // o Object to assign from
131  //
132  // Return:
133  // Reference to this object
134  //
136  return *this;
137 }
138 
139 //____________________________________________________________________
140 void
142 {
143  //
144  // Initialize the object
145  //
146  // Parameters:
147  // vtxAxis Vertex axis in use
148  //
149 
150  // We temporary disable the displaced vertices so we can initialize
151  // the routine ourselves.
152  Bool_t saveDisplaced = AllowDisplaced();
153  if (saveDisplaced) SetVertexMethod(kNormal);
155  if (saveDisplaced) SetVertexMethod(kDisplaced);
156 
157  Int_t maxPart = 450;
158  Int_t maxBin = 225;
159  Int_t maxB = 25;
160  fHVertex = 0;
161  if (vtxAxis.GetXbins() && vtxAxis.GetXbins()->GetArray())
162  fHVertex = new TH1F("vertex", "True vertex distribution",
163  vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray());
164  else
165  fHVertex = new TH1F("vertex", "True vertex distribution",
166  vtxAxis.GetNbins(), vtxAxis.GetXmin(),
167  vtxAxis.GetXmax());
168  fHVertex->SetXTitle("v_{z} [cm]");
169  fHVertex->SetYTitle("# of events");
170  fHVertex->SetFillColor(kGreen+1);
171  fHVertex->SetFillStyle(3001);
172  fHVertex->SetDirectory(0);
173  // fHVertex->Sumw2();
174  fList->Add(fHVertex);
175 
176  fHVertexXY = new TH2F("vertexXY", "True vertex distribution",
177  100, -1, 1, 100, -1, 1);
178  fHVertexXY->SetDirectory(0);
179  fHVertexXY->SetXTitle("x [cm]");
180  fHVertexXY->SetYTitle("y [cm]");
181  fList->Add(fHVertexXY);
182 
183  fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
184  fHPhiR->SetXTitle("#Phi_{R} [radians]");
185  fHPhiR->SetYTitle("# of events");
186  fHPhiR->SetFillColor(kGreen+1);
187  fHPhiR->SetFillStyle(3001);
188  fHPhiR->SetDirectory(0);
189  fList->Add(fHPhiR);
190 
191  fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
192  fHB->SetXTitle("b [fm]");
193  fHB->SetYTitle("# of events");
194  fHB->SetFillColor(kGreen+1);
195  fHB->SetFillStyle(3001);
196  fHB->SetDirectory(0);
197  fList->Add(fHB);
198 
199  fHMcC = static_cast<TH1F*>(fHCent->Clone("mcC"));
200  fHMcC->SetFillColor(kCyan+2);
201  fHMcC->SetDirectory(0);
202  fList->Add(fHMcC);
203 
204  fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
205  5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
206  fHBvsPart->SetXTitle("b [fm]");
207  fHBvsPart->SetYTitle("# of participants");
208  fHBvsPart->SetZTitle("Events");
209  fHBvsPart->SetDirectory(0);
210  fList->Add(fHBvsPart);
211 
212  fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
213  5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
214  fHBvsBin->SetXTitle("b [fm]");
215  fHBvsBin->SetYTitle("# of binary collisions");
216  fHBvsBin->SetZTitle("Events");
217  fHBvsBin->SetDirectory(0);
218  fList->Add(fHBvsBin);
219 
220  fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
221  5*maxB, 0, maxB, fCentAxis->GetNbins(),
222  fCentAxis->GetXbins()->GetArray());
223  fHBvsCent->SetXTitle("b [fm]");
224  fHBvsCent->SetYTitle("Centrality [%]");
225  fHBvsCent->SetZTitle("Event");
226  fHBvsCent->SetDirectory(0);
227  fList->Add(fHBvsCent);
228 
229 
230  fHVzComp = 0;
231  if (vtxAxis.GetXbins() && vtxAxis.GetXbins()->GetArray())
232  fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
233  vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray(),
234  vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray());
235  else
236  fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
237  10*vtxAxis.GetNbins(), vtxAxis.GetXmin(),
238  vtxAxis.GetXmax(), 10*vtxAxis.GetNbins(),
239  vtxAxis.GetXmin(), vtxAxis.GetXmax());
240  fHVzComp->SetXTitle("True v_{z} [cm]");
241  fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
242  fHVzComp->SetZTitle("Events");
243  fHVzComp->SetDirectory(0);
244  fList->Add(fHVzComp);
245 
246  fHCentVsPart = new TH2F("centralityVsParticipans",
247  "# of participants vs Centrality",
248  maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(),
249  fCentAxis->GetXbins()->GetArray());
250  fHCentVsPart->SetXTitle("Participants");
251  fHCentVsPart->SetYTitle("Centrality [%]");
252  fHCentVsPart->SetZTitle("Event");
253  fHCentVsPart->SetDirectory(0);
254  fList->Add(fHCentVsPart);
255 
256  fHCentVsBin = new TH2F("centralityVsBinary",
257  "# of binary collisions vs Centrality",
258  maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(),
259  fCentAxis->GetXbins()->GetArray());
260  fHCentVsBin->SetXTitle("Binary collisions");
261  fHCentVsBin->SetYTitle("Centrality [%]");
262  fHCentVsBin->SetZTitle("Event");
263  fHCentVsBin->SetDirectory(0);
264  fList->Add(fHCentVsBin);
265 
266  Int_t nC = fHCent->GetNbinsX();
267  Double_t cL = fHCent->GetXaxis()->GetXmin();
268  Double_t cH = fHCent->GetXaxis()->GetXmax();
269  fHCentVsMcC = new TH2F("centralityRecoVsMC",
270  "Centrality from reconstruction vs MC derived",
271  nC, cL, cH, nC, cL, cH);
272  fHCentVsMcC->SetDirectory(0);
273  fHCentVsMcC->SetStats(0);
274  fHCentVsMcC->SetXTitle("Centralty from Reco [%]");
275  fHCentVsMcC->SetYTitle("Centralty derived from Impact Par. [%]");
276  fHCentVsMcC->SetZTitle("Events");
277  fList->Add(fHCentVsMcC);
278 
280 }
281 
282 namespace
283 {
284  TString& AppendField(TString& s, bool test, const char* f)
285  {
286  if (!test) return s;
287  if (!s.IsNull()) s.Append(",");
288  s.Append(f);
289  return s;
290  }
291 }
292 
293 //____________________________________________________________________
294 void
296 {
297  //Assign MC only triggers : True NSD etc.
298  AliHeader* header = event->Header();
299  AliGenEventHeader* genHeader = header->GenEventHeader();
300  AliCollisionGeometry* colGeometry =
301  dynamic_cast<AliCollisionGeometry*>(genHeader);
302  AliGenPythiaEventHeader* pythiaHeader =
303  dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
304  AliGenDPMjetEventHeader* dpmHeader =
305  dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
306  AliGenGeVSimEventHeader* gevHeader =
307  dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
308  AliGenHerwigEventHeader* herwigHeader =
309  dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
310  AliGenHijingEventHeader* hijingHeader =
311  dynamic_cast<AliGenHijingEventHeader*>(genHeader);
312  // AliGenHydjetEventHeader* hydjetHeader =
313  // dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
314  AliGenEposEventHeader* eposHeader =
315  dynamic_cast<AliGenEposEventHeader*>(genHeader);
316 
317  AppendField(fProduction, colGeometry, "Geometry");
318  AppendField(fProduction, pythiaHeader, "Pythia");
319  AppendField(fProduction, dpmHeader, "DPM");
320  AppendField(fProduction, gevHeader, "GevSim");
321  AppendField(fProduction, herwigHeader, "Herwig");
322  AppendField(fProduction, hijingHeader, "Hijing");
323  // AppendField(fProduction, hydjetHeader, "Hydjet");
324  AppendField(fProduction, eposHeader, "EPOS");
325 }
326 
327 //____________________________________________________________________
328 UInt_t
330  UInt_t& triggers,
331  UShort_t& ivz,
332  TVector3& ip,
333  Double_t& b,
334  Double_t& c,
335  Int_t& npart,
336  Int_t& nbin,
337  Double_t& phiR)
338 {
339  //
340  // Process the event
341  //
342  // Parameters:
343  // event Input event
344  // triggers On return, the triggers fired
345  // ivz On return, the found vertex bin (1-based). A zero
346  // means outside of the defined vertex range
347  // vz On return, the z position of the interaction
348  // b On return, the impact parameter - < 0 if not available
349  // c On return, centrality estimate - < 0 if not available
350  // phiR On return, the event plane angle - < 0 if not available
351  //
352  // Return:
353  // 0 (or kOk) on success, otherwise a bit mask of error codes
354  //
355 
356  // Check that we have an event
357  if (!event) {
358  AliWarning("No MC event found for input event");
359  return kNoEvent;
360  }
361  // MC events are _always_ B collisions
362  triggers |= AliAODForwardMult::kB;
363  //Assign MC only triggers : True NSD etc.
364  AliHeader* header = event->Header();
365  AliGenEventHeader* genHeader = header->GenEventHeader();
366  AliCollisionGeometry* colGeometry =
367  dynamic_cast<AliCollisionGeometry*>(genHeader);
368  AliGenPythiaEventHeader* pythiaHeader =
369  dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
370  AliGenDPMjetEventHeader* dpmHeader =
371  dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
372  AliGenGeVSimEventHeader* gevHeader =
373  dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
374  AliGenHerwigEventHeader* herwigHeader =
375  dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
376  // AliGenHijingEventHeader* hijingHeader =
377  // dynamic_cast<AliGenHijingEventHeader*>(genHeader);
378  // AliGenHydjetEventHeader* hydjetHeader =
379  // dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
380  // AliGenEposEventHeader* eposHeader =
381  // dynamic_cast<AliGenEposEventHeader*>(genHeader);
382 
383  // Check if this is a single diffractive event
384  Bool_t sd = false;
385  Double_t phi = -1111;
386  Bool_t egSD = false;
387  npart = 0;
388  nbin = 0;
389  b = -1;
390  c = -1;
391  phiR = -1111;
392  if (colGeometry) {
393  b = colGeometry->ImpactParameter();
394  phi = colGeometry->ReactionPlaneAngle();
395  npart = (colGeometry->ProjectileParticipants() +
396  colGeometry->TargetParticipants());
397  nbin = colGeometry->NN();
398  }
399  if (fDebug && !colGeometry) {
400  AliWarningF("Collision header of class %s is not a CollisionHeader",
401  genHeader->ClassName());
402  }
403 
404  if(pythiaHeader) {
405  Int_t pythiaType = pythiaHeader->ProcessType();
406  egSD = true; // We have SD flag in EG
407  if (pythiaType <= 100) {
408  // Pythia6
409  // 92 and 93 are SD
410  // 94 is DD
411  if (pythiaType==92 || pythiaType==93) sd = true;
412  }
413  else {
414  // Pythia8
415  if (pythiaType==103 || pythiaType==104) sd = true;
416  }
417  b = pythiaHeader->GetImpactParameter();
418  npart = 2; // Always 2 protons
419  nbin = 1; // Always 1 binary collision
420  }
421  if (b >= 0) {
422 #if 0
423  if (b < 3.5) c = 2.5; //0-5%
424  else if (b < 4.95) c = 7.5; //5-10%
425  else if (b < 6.98) c = 15; //10-20%
426  else if (b < 8.55) c = 25; //20-30%
427  else if (b < 9.88) c = 35; //30-40%
428  else if (b < 11.04) c = 45; //40-50%
429  else c = 55; //50-60%
430 #else
431  // Updated 4th of November 2014 from
432  // cern.ch/twiki/bin/view/ALICE/CentStudies#Tables_with_centrality_bins_AN1
433  Float_t np=0;
434  Float_t nc=0;
435  if (0.00 >= b && b < 1.57) { c=0.5; np=403.8; nc=1861; }
436  else if (1.57 >= b && b < 2.22) { c=1.5; np=393.6; nc=1766; }
437  else if (2.22 >= b && b < 2.71) { c=2.5; np=382.9; nc=1678; }
438  else if (2.71 >= b && b < 3.13) { c=3.5; np=372; nc=1597; }
439  else if (3.13 >= b && b < 3.50) { c=4.5; np=361.1; nc=1520; }
440  else if (3.50 >= b && b < 4.94) { c=7.5; np=329.4; nc=1316; }
441  else if (4.94 >= b && b < 6.05) { c=12.5; np=281.2; nc=1032; }
442  else if (6.05 >= b && b < 6.98) { c=17.5; np=239; nc=809.8; }
443  else if (6.98 >= b && b < 7.81) { c=22.5; np=202.1; nc=629.6; }
444  else if (7.81 >= b && b < 8.55) { c=27.5; np=169.5; nc=483.7; }
445  else if (8.55 >= b && b < 9.23) { c=32.5; np=141; nc=366.7; }
446  else if (9.23 >= b && b < 9.88) { c=37.5; np=116; nc=273.4; }
447  else if (9.88 >= b && b < 10.47) { c=42.5; np=94.11; nc=199.4; }
448  else if (10.47 >= b && b < 11.04) { c=47.5; np=75.3; nc=143.1; }
449  else if (11.04 >= b && b < 11.58) { c=52.5; np=59.24; nc=100.1; }
450  else if (11.58 >= b && b < 12.09) { c=57.5; np=45.58; nc=68.46; }
451  else if (12.09 >= b && b < 12.58) { c=62.5; np=34.33; nc=45.79; }
452  else if (12.58 >= b && b < 13.05) { c=67.5; np=25.21; nc=29.92; }
453  else if (13.05 >= b && b < 13.52) { c=72.5; np=17.96; nc=19.08; }
454  else if (13.52 >= b && b < 13.97) { c=77.5; np=12.58; nc=12.07; }
455  else if (13.97 >= b && b < 14.43) { c=82.5; np=8.812; nc=7.682; }
456  else if (14.43 >= b && b < 14.96) { c=87.5; np=6.158; nc=4.904; }
457  else if (14.96 >= b && b < 15.67) { c=92.5; np=4.376; nc=3.181; }
458  else if (15.67 >= b && b < 20.00) { c=97.5; np=3.064; nc=1.994; }
459  // Be careful to round off
460  if (npart <= 0) npart = Int_t(np+.5);
461  if (nbin <= 0) nbin = Int_t(nc+.5)/2;
462 #endif
463  }
464  if(dpmHeader) { // Also an AliCollisionGeometry
465  Int_t processType = dpmHeader->ProcessType();
466  egSD = true; // We have SD flag in EG
467  // 1 & 4 are ND
468  // 5 & 6 are SD
469  // 7 is DD
470 #if 0
471  if (processType == 5 || processType == 6) sd = kTRUE;
472 #else
473  // The below - or rather a different implementation with some
474  // errors - was proposed by Cvetan - I don't think it's right
475  // though. See also
476  //
477  // https://cern.ch/twiki/pub/ALICE/PAPaperCentrality/normalization.pdf
478  // https://cern.ch/twiki/bin/view/ALICE/PAMCProductionStudies
479  //
480  Int_t nsd1=0, nsd2=0, ndd=0;
481  Int_t npP = dpmHeader->ProjectileParticipants();
482  Int_t npT = dpmHeader->TargetParticipants();
483  // Get the numbeer of single and double diffractive participants
484  dpmHeader->GetNDiffractive(nsd1,nsd2,ndd);
485  // Check if all partipants are single/double diffractive
486  if ((ndd == 0) && ((npP == nsd1) || (npT == nsd2))) sd = true;
487  // else if (ndd == (npP + npT)) dd = true;
488  // Printf("Projectile: %3d (%3d) Target: %3d (%3d) DD: %3d Process: %d",
489  // npP, nsd1, npT, nsd2, ndd, type);
490 #endif
491  }
492  if (gevHeader) {
493  phi = gevHeader->GetEventPlane();
494  }
495  if (herwigHeader) {
496  Int_t processType = herwigHeader->ProcessType();
497  egSD = true; // We have SD flag in EG
498  // This is a guess
499  if (processType == 5 || processType == 6) sd = kTRUE;
500  npart = 2; // Always 2 protons
501  nbin = 1; // Always 1 binary collision
502  }
503  // if (hijingHeader) {
504  // b = hijingHeader->ImpactParameter();
505  // phi = hijingHeader->ReactionPlaneAngle();
506  // }
507  // if (hydjetHeader) {
508  // b = hydjetHeader->ImpactParameter();
509  // phi = hydjetHeader->ReactionPlaneAngle();
510  // }
511  // if (eposHeader) {
512  // b = eposHeader->ImpactParameter();
513  // phi = eposHeader->ReactionPlaneAngle();
514  // }
515 
516  // Normalize event plane angle to [0,2pi]
517  if (phi <= -1111) phiR = -1;
518  else {
519  while (true) {
520  if (phi < 0) phi += 2*TMath::Pi();
521  else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
522  else break;
523  }
524  phiR = phi;
525  }
526 
527  // Do a check on particles - but only if the EG does not give it or
528  // the impact parameter is large enough
529  if (!egSD && (b < 0 || b > 15)) sd = IsSingleDiffractive(event->Stack());
530 
531  // Set NSD flag
532  if(!sd) {
533  triggers |= AliAODForwardMult::kMCNSD;
534  fHTriggers->Fill(kMCNSD+0.5);
535  }
536 
537  // Get the primary vertex from EG
538  TArrayF vtx;
539  genHeader->PrimaryVertex(vtx);
540  ip.SetXYZ(vtx[0], vtx[1], vtx[2]);
541 
542  DMSG(fDebug, 1, "ip=(%f,%f,%f), phiR=%f, b=%f, npart=%d, nbin=%d",
543  vtx[0], vtx[1], vtx[2], phiR, b, npart, nbin);
544 
545  fHVertex->Fill(vtx[2]);
546  fHVertexXY->Fill(vtx[0], vtx[1]);
547  fHPhiR->Fill(phiR);
548  fHB->Fill(b);
549  fHMcC->Fill(c);
550  fHBvsPart->Fill(b, npart);
551  fHBvsBin->Fill(b, nbin);
552 
553  if(AllowDisplaced()) {
554 #if 0
555  // Put the vertex at fixed locations
556  Double_t zvtx = vz;
557  Double_t ratio = zvtx/37.5;
558  if(ratio > 0) ratio = ratio + 0.5;
559  if(ratio < 0) ratio = ratio - 0.5;
560  Int_t ratioInt = Int_t(ratio);
561  zvtx = 37.5*((Double_t)ratioInt);
562  if(TMath::Abs(zvtx) > 999)
563  return kBadVertex;
564 #endif
565  if (!fDisplacedVertex.ProcessMC(event))
566  return kBadVertex;
568  ip.SetZ(fDisplacedVertex.GetVertexZ());
569  }
570 
571  // Check for the vertex bin
572  ivz = fVtxAxis.FindBin(ip.Z());
573  if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
574  if (fDebug > 3) {
575  AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
576  ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); }
577  ivz = 0;
578  return kBadVertex;
579  }
580 
581 
582  return kOk;
583 }
584 //____________________________________________________________________
585 Bool_t
587  Double_t& cent,
588  UShort_t& qual) const
589 {
590  //
591  // Read centrality from event
592  //
593  // Parameters:
594  // esd Event
595  // cent On return, the centrality or negative if not found
596  //
597  // Return:
598  // False on error, true otherwise
599  //
600  Bool_t ret = AliFMDEventInspector::ReadCentrality(esd, cent, qual);
601  if (qual != 0) {
602  AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
603  if (!centObj) return ret;
604 
605  // For MC, we allow `bad' centrality selections
606  cent = centObj->GetCentralityPercentileUnchecked(fCentMethod);
607  }
608  return ret;
609 }
610 
611 //____________________________________________________________________
612 namespace {
613  Double_t rapidity(TParticle* p, Double_t mass)
614  {
615  Double_t pt = p->Pt();
616  Double_t pz = p->Pz();
617  Double_t ee = TMath::Sqrt(pt*pt+pz*pz+mass*mass);
618  if (TMath::Abs(ee - TMath::Abs(pz)) < 1e-9) return TMath::Sign(1e30, pz);
619  return .5 * TMath::Log((ee + pz) / (ee - pz));
620  }
621 }
622 
623 //____________________________________________________________________
624 Bool_t
626  Double_t xiMin,
627  Double_t xiMax) const
628 {
629  // Re-implementation of AliPWG0Helper::IsHadronLevelSingleDiffrative
630  //
631  // This is re-implemented here to be indendent of the PWG0 library.
632  TParticle* p1 = 0; // Particle with least y
633  TParticle* p2 = 0; // Particle with largest y
634  Double_t y1 = 1e10; // y of p1
635  Double_t y2 = -1e10; // y of p2
636 
637  // Loop over primaries
638  for (Int_t i = 0; i < stack->GetNprimary(); i++) {
639  TParticle* p = stack->Particle(i);
640  if (!p) continue;
641 
642  Int_t pdg = TMath::Abs(p->GetPdgCode());
643  Int_t c1 = p->GetFirstDaughter();
644  Int_t s = p->GetStatusCode();
645  Int_t mfl = Int_t(pdg/TMath::Power(10,Int_t(TMath::Log10(pdg))));
646 
647  // Select final state charm and beauty
648  if (c1 > -1 || s != 1) mfl = 0;
649 
650  // Check if this is a primary, pi0, Sigma0, ???, or most
651  // significant digit is larger than or equal to 4
652  if (!(stack->IsPhysicalPrimary(i) ||
653  pdg == 111 ||
654  pdg == 3212 ||
655  pdg == 3124 ||
656  mfl >= 4)) continue;
657 
658  Int_t m1 = p->GetFirstMother();
659  if (m1 > 0) {
660  TParticle* pm1 = stack->Particle(m1);
661  Int_t mpdg = TMath::Abs(pm1->GetPdgCode());
662  // Check if mother is a p0, Simga0, or ???
663  if (mpdg == 111 || mpdg == 3124 || mpdg == 3212) continue;
664  }
665 
666  // Calculate the rapidity of the particle
667  Double_t mm = (pdg != 3124 && p->GetPDG() ? p->GetMass() : 1.5195);
668  Double_t yy = rapidity(p, mm);
669 
670  // Check if the rapidity of this particle is further out than any
671  // of the preceding particles
672  if (yy < y1) {
673  y1 = yy;
674  p1 = p;
675  }
676  if (yy > y2) {
677  y2 = yy;
678  p2 = p;
679  }
680  }
681  if (!p1 || !p2) return false;
682 
683  // Calculate rapidities assuming protons
684  y1 = TMath::Abs(rapidity(p1, 0.938));
685  y2 = TMath::Abs(rapidity(p2, 0.938));
686 
687  // Check if both or just one is a proton
688  Int_t pdg1 = p1->GetPdgCode();
689  Int_t pdg2 = p2->GetPdgCode();
690  Int_t arm = -99999;
691  if (pdg1 == 2212 && pdg2 == 2212)
692  arm = (y1 > y2 ? 0 : 1);
693  else if (pdg1 == 2212)
694  arm = 0;
695  else if (pdg2 == 2212)
696  arm = 1;
697  else
698  return false;
699 
700  // Rapidity shift
701  Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0);
702  Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
703 
704  if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
705  if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
706 
707  return false;
708 }
709 
710 //____________________________________________________________________
711 Bool_t
713  Double_t cent, Double_t mcC,
714  Double_t b,
715  Int_t npart, Int_t nbin)
716 {
717  fHVzComp->Fill(trueVz, vz);
718  fHBvsCent->Fill(b, cent);
719  fHCentVsPart->Fill(npart, cent);
720  fHCentVsBin->Fill(nbin, cent);
721  fHCentVsMcC->Fill(cent, mcC);
722 
723  return true;
724 }
725 
726 
727 //
728 // EOF
729 //
730 
Int_t pdg
double Double_t
Definition: External.C:58
UShort_t fEnergy
Histogram container.
Definition: External.C:236
virtual Bool_t CompareResults(Double_t vz, Double_t trueVz, Double_t cent, Double_t mcC, Double_t b, Int_t npart, Int_t nbin)
virtual Bool_t ReadCentrality(const AliESDEvent &esd, Double_t &cent, UShort_t &qual) const
Double_t mass
TH1I * fHTriggers
XY vtx with trigger and Z vertex in range.
#define DMSG(L, N, F,...)
TCanvas * c
Definition: TestFitELoss.C:172
AliFMDEventInspector & operator=(const AliFMDEventInspector &)
AliStack * stack
Int_t cH
Definition: Combine.C:26
Per-event per bin.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Various utilities used in PWGLF/FORWARD.
void SetVertexMethod(EVtxType t)
Bool_t ProcessMC(const AliMCEvent *mcevent)
rapidity
Definition: HFPtSpectrum.C:47
virtual Bool_t ReadCentrality(const AliESDEvent &esd, Double_t &cent, UShort_t &qual) const
virtual void ReadProductionDetails(AliMCEvent *event)
TString fCentMethod
Configured background words.
void SetupForData(const TAxis &vtxAxis)
AliDisplacedVertexSelection fDisplacedVertex
AliFMDMCEventInspector & operator=(const AliFMDMCEventInspector &)
UInt_t ProcessMC(AliMCEvent *event, UInt_t &triggers, UShort_t &ivz, TVector3 &ip, Double_t &b, Double_t &c, Int_t &npart, Int_t &nbin, Double_t &phiR)
Bool_t IsSingleDiffractive(AliStack *stack, Double_t xiMin=0, Double_t xiMax=1./81) const
unsigned short UShort_t
Definition: External.C:28
void test(int runnumber=195345)
TH1F * fHCent
Trigger words.
bool Bool_t
Definition: External.C:53
void SetupForData(TList *l, const char *name=0, Bool_t mc=false)
Bool_t AllowDisplaced() const
virtual void SetupForData(const TAxis &vtxAxis)