AliPhysics  e34b7ac (e34b7ac)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliForwardUtil.cxx
Go to the documentation of this file.
1 //
2 // Utilities used in the forward multiplcity analysis
3 //
4 //
5 #include "AliForwardUtil.h"
6 //#include <ARVersion.h>
7 #include <AliAnalysisManager.h>
8 #include "AliAODForwardMult.h"
9 #include <AliLog.h>
10 #include <AliInputEventHandler.h>
11 #include <AliAODInputHandler.h>
12 #include <AliAODHandler.h>
13 #include <AliAODEvent.h>
14 #include <AliESDEvent.h>
15 #include <AliAnalysisTaskSE.h>
16 #include <AliPhysicsSelection.h>
17 #include <AliTriggerAnalysis.h>
18 #include <AliMultiplicity.h>
19 #include "AliMultEstimator.h" // <-- Sigh, needed by AliMultSelection
20 #include "AliMultVariable.h" // <-- Sigh, needed by AliMultSelection
21 #include "AliMultInput.h" // <-- Sigh, needed by AliMultSelection
22 #include "AliMultSelection.h"
23 #include "AliMultSelectionCuts.h"
24 #include <TParameter.h>
25 #include <TH2D.h>
26 #include <TH1I.h>
27 #include <TF1.h>
28 #include <TFitResult.h>
29 #include <TMath.h>
30 #include <TError.h>
31 #include <TROOT.h>
32 #include <TVector3.h>
33 #define FIT_OPTIONS "RNS"
34 
35 //====================================================================
37 {
38 #ifdef ALIROOT_SVN_REVISION
39  return ALIROOT_SVN_REVISION;
40 #elif defined(ALIROOT_REVISION)
41  static ULong_t ret = 0;
42  if (ret != 0) return ret;
43 
44  // Select first 32bits of the 40byte long check-sum
45  TString rev(ALIROOT_REVISION, 8);
46  for (ULong_t i = 0; i < 8; i++) {
47  ULong_t p = 0;
48  switch (rev[i]) {
49  case '0': p = 0; break;
50  case '1': p = 1; break;
51  case '2': p = 2; break;
52  case '3': p = 3; break;
53  case '4': p = 4; break;
54  case '5': p = 5; break;
55  case '6': p = 6; break;
56  case '7': p = 7; break;
57  case '8': p = 8; break;
58  case '9': p = 9; break;
59  case 'a': case 'A': p = 10; break;
60  case 'b': case 'B': p = 11; break;
61  case 'c': case 'C': p = 12; break;
62  case 'd': case 'D': p = 13; break;
63  case 'e': case 'E': p = 14; break;
64  case 'f': case 'F': p = 15; break;
65  }
66  ret |= (p << (32-4*(i+1)));
67  }
68  return ret;
69 #else
70  return 0;
71 #endif
72 }
73 //____________________________________________________________________
75 {
76  // Do something here when we switch to git - sigh!
77 #if !defined(ALIROOT_SVN_BRANCH) && !defined(ALIROOT_BRANCH)
78  return 0;
79 #endif
80  static ULong_t ret = 0;
81  if (ret != 0) return ret;
82  TString str;
83  TString top;
84 #ifdef ALIROOT_SVN_BRANCH
85  str = ALIROOT_SVN_BRANCH;
86  top = "trunk";
87 #elif defined(ALIROOT_BRANCH)
88  str = ALIROOT_BRANCH;
89  top = "master";
90 #endif
91  if (str.IsNull()) return 0xFFFFFFFF;
92  if (str[0] == 'v') str.Remove(0,1);
93  if (str.EqualTo(top)) return ret = 0xFFFFFFFF;
94 
95  TObjArray* tokens = str.Tokenize("-");
96  TObjString* pMajor = tokens->GetEntries()>0 ?
97  (static_cast<TObjString*>(tokens->At(0))) : 0;
98  TObjString* pMinor = tokens->GetEntries()>1 ?
99  (static_cast<TObjString*>(tokens->At(1))) : 0;
100  TObjString* pRelea = tokens->GetEntries() > 2 ?
101  static_cast<TObjString*>(tokens->At(2)) : 0;
102  TObjString* pAn = tokens->GetEntries() > 3 ?
103  static_cast<TObjString*>(tokens->At(3)) : 0;
104  TString sMajor,sMinor,sRelea;
105  if (pMajor) sMajor = pMajor->String().Strip(TString::kLeading, '0');
106  if (pMinor) sMinor = pMinor->String().Strip(TString::kLeading, '0');
107  if (pRelea) sRelea = pRelea->String().Strip(TString::kLeading, '0');
108  //
109  ret = (((sMajor.Atoi() & 0xFF) << 12) |
110  ((sMinor.Atoi() & 0xFF) << 8) |
111  ((sRelea.Atoi() & 0xFF) << 4) |
112  (pAn ? 0xAA : 0));
113 
114  return ret;
115 }
116 
117 //====================================================================
118 UShort_t
120 {
121  //
122  // Parse a collision system spec given in a string. Known values are
123  //
124  // - "ppb", "p-pb", "pa", "p-a" which returns kPPb
125  // - "pp", "p-p" which returns kPP
126  // - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
127  // - Everything else gives kUnknown
128  //
129  // Parameters:
130  // sys Collision system spec
131  //
132  // Return:
133  // Collision system id
134  //
135  TString s(sys);
136  s.ToLower();
137  // we do pA first to avoid pp catch on ppb string (AH)
138  if (s.Contains("p-pb") || s.Contains("ppb")) return AliForwardUtil::kPPb;
139  if (s.Contains("p-a") || s.Contains("pa")) return AliForwardUtil::kPPb;
140  if (s.Contains("pb-p") || s.Contains("pbp")) return AliForwardUtil::kPbp;
141  if (s.Contains("a-p") || s.Contains("ap")) return AliForwardUtil::kPbp;
142  if (s.Contains("p-p") || s.Contains("pp")) return AliForwardUtil::kPP;
143  if (s.Contains("pb-pb") || s.Contains("pbpb")) return AliForwardUtil::kPbPb;
144  if (s.Contains("a-a") || s.Contains("aa")) return AliForwardUtil::kPbPb;
146 }
147 //____________________________________________________________________
148 const char*
150 {
151  //
152  // Get a string representation of the collision system
153  //
154  // Parameters:
155  // sys Collision system
156  // - kPP -> "pp"
157  // - kPbPb -> "PbPb"
158  // - kPPb -> "pPb"
159  // - kPbp -> "Pbp"
160  // - anything else gives "unknown"
161  //
162  // Return:
163  // String representation of the collision system
164  //
165  switch (sys) {
166  case AliForwardUtil::kPP: return "pp";
167  case AliForwardUtil::kPbPb: return "PbPb";
168  case AliForwardUtil::kPPb: return "pPb";
169  case AliForwardUtil::kPbp: return "Pbp";
170  }
171  return "unknown";
172 }
173 //____________________________________________________________________
174 Float_t
176 {
177  const Double_t pMass = 9.38271999999999995e-01;
178  const Double_t nMass = 9.39564999999999984e-01;
179  Double_t beamE = z * beam / 2;
180  Double_t beamM = z * pMass + (a - z) * nMass;
181  Double_t beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
182  Double_t beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
183  return beamY;
184 }
185 //____________________________________________________________________
186 Float_t
188  UShort_t z1,
189  UShort_t a1,
190  Short_t z2,
191  Short_t a2)
192 {
193  // Calculate the center of mass energy given target/projectile
194  // mass and charge numbers
195  if (z2 < 0) z2 = z1;
196  if (a2 < 0) a2 = a1;
197  return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
198 }
199 //____________________________________________________________________
200 Float_t
202  UShort_t a1,
203  Short_t z2,
204  Short_t a2)
205 {
206  // Calculate the center of mass rapidity (shift) given target/projectile
207  // mass and charge numbers
208  if (z2 < 0) z2 = z1;
209  if (a2 < 0) a2 = a1;
210  if (z2 == z1 && a2 == a1) return 0;
211  return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
212 }
213 
214 namespace {
215  UShort_t CheckSNN(Float_t energy)
216  {
217  if (TMath::Abs(energy - 900.) < 10) return 900;
218  if (TMath::Abs(energy - 2400.) < 10) return 2400;
219  if (TMath::Abs(energy - 2760.) < 20) return 2760;
220  if (TMath::Abs(energy - 4400.) < 10) return 4400;
221  if (TMath::Abs(energy - 5000.) < 10) return 5000;
222  if (TMath::Abs(energy - 5022.) < 10) return 5023;
223  if (TMath::Abs(energy - 5125.) < 30) return 5100;
224  if (TMath::Abs(energy - 5500.) < 40) return 5500;
225  if (TMath::Abs(energy - 7000.) < 10) return 7000;
226  if (TMath::Abs(energy - 8000.) < 10) return 8000;
227  if (TMath::Abs(energy - 10000.) < 10) return 10000;
228  if (TMath::Abs(energy - 13000.) < 10) return 13000;
229  if (TMath::Abs(energy - 14000.) < 10) return 14000;
230  return 0;
231  }
232 }
233 //____________________________________________________________________
234 UShort_t
236 {
237  //
238  // Parse the center of mass energy given as a float and return known
239  // values as a unsigned integer
240  //
241  // Parameters:
242  // sys Collision system (needed for AA)
243  // beam Center of mass energy * total charge
244  //
245  // Return:
246  // Center of mass energy per nucleon
247  //
248  Float_t energy = beam;
249  // Below no longer needed apparently
250  // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
251  if (sys == AliForwardUtil::kPPb)
252  energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
253  if (sys == AliForwardUtil::kPbp)
254  energy = CenterOfMassEnergy(beam, 1, 1, 82, 208);
255  else if (sys == AliForwardUtil::kPbPb)
256  energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
257  UShort_t ret = CheckSNN(energy);
258  if (ret > 1) return ret;
259  if (sys == AliForwardUtil::kPbPb ||
260  sys == AliForwardUtil::kPPb ||
261  sys == AliForwardUtil::kPbp) {
262  ret = CheckSNN(beam);
263  }
264  return ret;
265 }
266 //____________________________________________________________________
267 const char*
269 {
270  //
271  // Get a string representation of the center of mass energy per nuclean
272  //
273  // Parameters:
274  // cms Center of mass energy per nucleon
275  //
276  // Return:
277  // String representation of the center of mass energy per nuclean
278  //
279  return Form("%04dGeV", cms);
280 }
281 //____________________________________________________________________
282 Short_t
284 {
285  //
286  // Parse the magnetic field (in kG) as given by a floating point number
287  //
288  // Parameters:
289  // field Magnetic field in kG
290  //
291  // Return:
292  // Short integer value of magnetic field in kG
293  //
294  if (TMath::Abs(v - 5.) < 1 ) return +5;
295  if (TMath::Abs(v + 5.) < 1 ) return -5;
296  if (TMath::Abs(v) < 1) return 0;
297  return 999;
298 }
299 //____________________________________________________________________
300 const char*
302 {
303  //
304  // Get a string representation of the magnetic field
305  //
306  // Parameters:
307  // field Magnetic field in kG
308  //
309  // Return:
310  // String representation of the magnetic field
311  //
312  return Form("%01dkG", f);
313 }
314 //_____________________________________________________________________
316 {
317  // Check if AOD is the output event
318  if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
319 
320  AliAODEvent* ret = task->AODEvent();
321  if (ret) return ret;
322 
323  // Check if AOD is the input event
324  ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
325  if (!ret) ::Warning("GetAODEvent", "No AOD event found");
326 
327  return ret;
328 }
329 //_____________________________________________________________________
331 {
332  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
333  if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
334  // ::Info("CheckForAOD", "Found AOD Input handler");
335  return 1;
336  }
337  if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
338  // ::Info("CheckForAOD", "Found AOD Output handler");
339  return 2;
340  }
341 
342  ::Warning("CheckForAOD",
343  "Neither and input nor output AOD handler is specified");
344  return 0;
345 }
346 //_____________________________________________________________________
347 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
348 {
349  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
350  if (!cls) {
351  AliAnalysisTask* t = am->GetTask(clsOrName);
352  if (!t) {
353  ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
354  return false;
355  }
356  ::Info("CheckForTask", "Found task %s", clsOrName);
357  return true;
358  }
359  TClass* dep = gROOT->GetClass(clsOrName);
360  if (!dep) {
361  ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
362  return false;
363  }
364  TIter next(am->GetTasks());
365  TObject* o = 0;
366  while ((o = next())) {
367  if (o->IsA()->InheritsFrom(dep)) {
368  ::Info("CheckForTask", "Found task of class %s: %s",
369  clsOrName, o->GetName());
370  return true;
371  }
372  }
373  ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
374  return false;
375 }
376 
377 //_____________________________________________________________________
379 {
380  TParameter<int>* ret = new TParameter<int>(name, value);
381  ret->SetMergeMode('f');
382  ret->SetUniqueID(value);
383  return ret;
384 }
385 //_____________________________________________________________________
386 TObject* AliForwardUtil::MakeParameter(const char* name, Int_t value)
387 {
388  TParameter<int>* ret = new TParameter<int>(name, value);
389  ret->SetMergeMode('f');
390  ret->SetUniqueID(value);
391  return ret;
392 }
393 //_____________________________________________________________________
395 {
396  TParameter<Long_t>* ret = new TParameter<Long_t>(name, value);
397  ret->SetMergeMode('f');
398  ret->SetUniqueID(value);
399  return ret;
400 }
401 //_____________________________________________________________________
403 {
404  TParameter<double>* ret = new TParameter<double>(name, value);
405  // Float_t v = value;
406  // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
407  ret->SetMergeMode('f');
408  // ret->SetUniqueID(*tmp);
409  return ret;
410 }
411 //_____________________________________________________________________
413 {
414  TParameter<bool>* ret = new TParameter<bool>(name, value);
415  ret->SetMergeMode('f');
416  ret->SetUniqueID(value);
417  return ret;
418 }
419 
420 //_____________________________________________________________________
422 {
423  if (!o) return;
424  TParameter<int>* p = static_cast<TParameter<int>*>(o);
425  if (p->TestBit(BIT(19)))
426  value = p->GetVal();
427  else
428  value = o->GetUniqueID();
429 }
430 //_____________________________________________________________________
432 {
433  if (!o) return;
434  TParameter<int>* p = static_cast<TParameter<int>*>(o);
435  if (p->TestBit(BIT(19)))
436  value = p->GetVal();
437  else
438  value = o->GetUniqueID();
439 }
440 //_____________________________________________________________________
442 {
443  if (!o) return;
444  TParameter<Long_t>* p = static_cast<TParameter<Long_t>*>(o);
445  if (p->TestBit(BIT(19)))
446  value = p->GetVal();
447  else
448  value = o->GetUniqueID();
449 }
450 //_____________________________________________________________________
452 {
453  if (!o) return;
454  TParameter<double>* p = static_cast<TParameter<double>*>(o);
455  if (p->TestBit(BIT(19)))
456  value = p->GetVal(); // o->GetUniqueID();
457  else {
458  UInt_t i = o->GetUniqueID();
459  Float_t v = *reinterpret_cast<Float_t*>(&i);
460  value = v;
461  }
462 }
463 //_____________________________________________________________________
465 {
466  if (!o) return;
467  TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
468  if (p->TestBit(BIT(19)))
469  value = p->GetVal(); // o->GetUniqueID();
470  else
471  value = o->GetUniqueID();
472 }
473 
474 //_____________________________________________________________________
476 {
477  // Get max R of ring
478  //
479  // New implementation has only one branch
480  const Double_t minR[] = { 4.5213, 15.4 };
481  const Double_t maxR[] = { 17.2, 28.0 };
482  const Int_t nStr[] = { 512, 256 };
483 
484  Int_t q = (ring == 'I' || ring == 'i') ? 0 : 1;
485  Double_t rad = maxR[q] - minR[q];
486  Double_t segment = rad / nStr[q];
487  Double_t r = minR[q] + segment*strip;
488 
489  return r;
490 }
491 
492 //_____________________________________________________________________
494 {
495  Int_t hybrid = sec / 2;
496  Int_t q = (ring == 'I' || ring == 'i') ? 0 : 1;
497  Int_t r = q == 0 ? 1 : 0;
498  const Double_t zs[][2] = { { 320.266+1.5, kInvalidValue },
499  { 83.666, 74.966+.5 },
500  { -63.066+.5, -74.966 } };
501  if (det > 3 || zs[det-1][q] == kInvalidValue) {
502  ::Warning("GetSectorZ", "Unknown sub-detector FMD%d%c", det, ring);
503  return kInvalidValue;
504  }
505 
506  Double_t z = zs[det-1][q];
507  switch (det) {
508  case 1: if ((hybrid % 2) == 1) z -= .5; break;
509  case 2: if ((hybrid % 2) == r) z -= .5; break;
510  case 3: if ((hybrid % 2) == q) z -= .5; break;
511  }
512 
513  return z;
514 }
515 
516 //_____________________________________________________________________
518 {
519  UShort_t nSec = (ring == 'I' || ring == 'i') ? 20 : 40;
520  Double_t base = float(sec+.5) / nSec * TMath::TwoPi();
521  switch (d) {
522  case 1: base += TMath::Pi()/2; break;
523  case 2: break;
524  case 3: base = TMath::Pi() - base; break;
525  default:
526  ::Warning("GetSectorPhi", "Unknown detector %d", d);
527  return kInvalidValue;
528  }
529  if (base < 0) base += TMath::TwoPi();
530  if (base > TMath::TwoPi()) base -= TMath::TwoPi();
531  return base;
532 }
533 
534 //_____________________________________________________________________
536  UShort_t sec, UShort_t strip,
537  Double_t zvtx)
538 {
539  // Calculate eta from strip with vertex (redundant with
540  // AliESDFMD::Eta but support displaced vertices)
541  //
542  // Slightly more optimized version that uses less branching
543 
544  // Get R of the strip
545  Double_t r = GetStripR(ring, strip);
546  Double_t z = GetSectorZ(det, ring, sec);
547  Double_t theta = TMath::ATan2(r,z-zvtx);
548  Double_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
549 
550  return eta;
551 }
552 
553 //_____________________________________________________________________
555  UShort_t sec, UShort_t strip,
556  const TVector3& ip,
557  TVector3& pos)
558 {
559  Double_t rD = GetStripR(ring, strip);
560  Double_t phiD = GetSectorPhi(det,ring, sec);
561  Double_t zD = GetSectorZ(det, ring, sec);
562  if (phiD == kInvalidValue || zD == kInvalidValue) {
564  return false;
565  }
566  Double_t xD = rD*TMath::Cos(phiD);
567  Double_t yD = rD*TMath::Sin(phiD);
568  Double_t iX = ip.X(); if (iX > 100) iX = 0; // No X
569  Double_t iY = ip.Y(); if (iY > 100) iY = 0; // No Y
570  Double_t dX = xD-iX;
571  Double_t dY = yD-iY;
572  Double_t dZ = zD-ip.Z();
573  pos.SetXYZ(dX, dY, dZ);
574  return true;
575 }
576 //_____________________________________________________________________
578  UShort_t sec, UShort_t strip,
579  const TVector3& ip,
580  Double_t& eta, Double_t& phi)
581 {
582  TVector3 pos;
583  if (!GetXYZ(det, ring, sec, strip, ip, pos)) {
584  ::Warning("GetEtaPhi", "Invalid position for FMD%d%c[%2d,%3d]=(%f,%f,%f)",
585  det, ring, sec, strip, pos.X(), pos.Y(), pos.Z());
586  eta = kInvalidValue;
587  phi = kInvalidValue;
588  return false;
589  }
590  Double_t r = TMath::Sqrt(TMath::Power(pos.X(),2)+
591  TMath::Power(pos.Y(),2));
592  Double_t theta = TMath::ATan2(r, pos.Z());
593  Double_t tant = TMath::Tan(theta/2);
594  if (TMath::Abs(theta) < 1e-9) {
595  ::Warning("GetEtaPhi","tan(theta/2)=%f very small", tant);
596  eta = kInvalidValue;
597  phi = kInvalidValue;
598  return false;
599  }
600  phi = TMath::ATan2(pos.Y(), pos.X());
601  eta = -TMath::Log(tant);
602  if (phi < 0) phi += TMath::TwoPi();
603  if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
604 
605  return true;
606 }
607 
608 //_____________________________________________________________________
610  UShort_t strip,
611  Double_t& eta, Double_t& phi ,
612  Double_t ipX, Double_t ipY)
613 {
614  Double_t rs = GetStripR(r, strip);
615  Double_t sx = rs*TMath::Cos(phi);
616  Double_t sy = rs*TMath::Sin(phi);
617  Double_t dx = sx-ipX;
618  Double_t dy = sy-ipY;
619  Double_t rv = TMath::Sqrt(TMath::Power(dx,2) + TMath::Power(dy,2));
620  Double_t the = 2*TMath::ATan(TMath::Exp(-eta));
621  Double_t tth = TMath::Tan(the);
622  if (TMath::Abs(tth) < 1e-9) {
623  ::Warning("GetEtaPhiFromStrip",
624  "eta=%f -> theta=%f tan(theta)=%f invalid (no change)",
625  eta, the, tth);
626  return false;
627  }
628  Double_t z = rs / tth;
629  // Printf("IP(x,y)=%f,%f S(x,y)=%f,%f D(x,y)=%f,%f R=%f theta=%f "
630  // "tan(theta)=%f z=%f",
631  // ipX, ipY, sx, sy, dx, dy, rv, the, TMath::Tan(the), z);
632  eta = -TMath::Log(TMath::Tan(TMath::ATan2(rv,z)/2));
633  phi = TMath::ATan2(dy,dx);
634  if (phi < 0) phi += TMath::TwoPi();
635 
636  return true;
637 }
638 
639 
640 //_____________________________________________________________________
642  Double_t phi,
643  Double_t xvtx, Double_t yvtx)
644 {
645  // Calculate eta from strip with vertex (redundant with
646  // AliESDFMD::Eta but support displaced vertices)
647 
648  // Unknown x,y -> no change
649  if (yvtx > 999 || xvtx > 999) return phi;
650 
651  //Get max R of ring
652  Double_t r = GetStripR(ring, strip);
653  Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
654  Double_t pha = (TMath::Abs(yvtx) < 1e-12 ? 0 : TMath::ATan2(xvtx, yvtx));
655  Double_t cha = amp * TMath::Cos(phi+pha);
656  phi += cha;
657  if (phi < 0) phi += TMath::TwoPi();
658  if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
659  return phi;
660 }
661 //====================================================================
662 TAxis*
664 {
665  TArrayD bins;
666  MakeFullIpZAxis(nCenter, bins);
667  TAxis* a = new TAxis(bins.GetSize()-1,bins.GetArray());
668  return a;
669 }
670 //____________________________________________________________________
671 void
673 {
674  // Custom vertex axis that will include satellite vertices
675  // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10
676  // Nominal vertices are usually in -10 to 10 and we should have
677  // 10 bins in that range. That gives us a total of
678  //
679  // 10+10+10=30 bins
680  //
681  // or 31 bin boundaries
682  if (nCenter % 2 == 1)
683  // Number of central bins is odd - make it even
684  nCenter--;
685  const Double_t mCenter = 20;
686  const Int_t nSat = 10;
687  const Int_t nBins = 2*nSat + nCenter;
688  const Int_t mBin = nBins / 2;
689  Double_t dCenter = 2*mCenter / nCenter;
690  bins.Set(nBins+1);
691  bins[mBin] = 0;
692  for (Int_t i = 1; i <= nCenter/2; i++) {
693  // Assign from the middle out
694  Double_t v = i * dCenter;
695  // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-i,mBin+i);
696  bins[mBin-i] = -v;
697  bins[mBin+i] = +v;
698  }
699  for (Int_t i = 1; i <= nSat; i++) {
700  Double_t v = (i+.5) * 37.5;
701  Int_t o = nCenter/2+i;
702  // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-o,mBin+o);
703  bins[mBin-o] = -v;
704  bins[mBin+o] = +v;
705  }
706 }
707 //____________________________________________________________________
708 void
710  Int_t minOrder,
711  Int_t maxOrder,
712  TArrayD& bins)
713 {
714  Double_t dO = Double_t(maxOrder-minOrder) / nBins;
715  bins.Set(nBins+1);
716  for (Int_t i = 0; i <= nBins; i++) bins[i] = TMath::Power(10, i * dO+minOrder);
717 }
718 
719 //____________________________________________________________________
720 void
722 {
723  Int_t ind = gROOT->GetDirLevel();
724  if (ind > 0)
725  // Print indention
726  std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
727 
728  TString t = TString::Format("%s %s", o.GetName(), o.ClassName());
729  const Int_t maxN = 75;
730  std::cout << "--- " << t << " " << std::setfill('-')
731  << std::setw(maxN-ind-5-t.Length()) << "-" << std::endl;
732 }
733 //____________________________________________________________________
734 void
735 AliForwardUtil::PrintName(const char* name)
736 {
737  Int_t ind = gROOT->GetDirLevel();
738  if (ind > 0)
739  // Print indention
740  std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
741 
742  // Now print field name
743  const Int_t maxN = 29;
744  Int_t width = maxN - ind;
745  TString n(name);
746  if (n.Length() > width-1) {
747  // Truncate the string, and put in "..."
748  n.Remove(width-4);
749  n.Append("...");
750  }
751  n.Append(":");
752  std::cout << std::setfill(' ') << std::left << std::setw(width)
753  << n << std::right << std::flush;
754 }
755 //____________________________________________________________________
756 void
757 AliForwardUtil::PrintField(const char* name, const char* value, ...)
758 {
759  PrintName(name);
760 
761  // Now format the field value
762  va_list ap;
763  va_start(ap, value);
764  static char buf[512];
765  vsnprintf(buf, 511, value, ap);
766  buf[511] = '\0';
767  va_end(ap);
768 
769  std::cout << buf << std::endl;
770 }
771 
772 
773 //====================================================================
774 Float_t AliForwardUtil::GetCentrality(const AliVEvent& event,
775  const TString& method,
776  Int_t& qual,
777  Bool_t verbose)
778 {
779  qual = 0xFFFF;
780  Float_t cent = GetCentralityMult(event,method,qual,verbose);
781  if (qual < 0xFFF) return cent;
782 
783  // No selection object found, try compat
784  return GetCentralityCompat(event,method,qual,verbose);
785 }
786 //____________________________________________________________________
788  const TString& method,
789  Int_t& qual,
790  Bool_t verbose)
791 {
792  TObject* o = event.FindListObject("MultSelection");
793  if (!o) {
794  if (verbose)
795  ::Warning("AliForwardUtil::GetCentralityMult",
796  "No MultSelection object found in event");
797  return -1;
798  }
799  AliMultSelection* sel = static_cast<AliMultSelection*>(o);
800  if (!sel->GetEstimatorList() ||
801  sel->GetEstimatorList()->GetEntries() <= 0){
802  if (verbose) {
803  ::Warning("AliForwardUtil::GetCentralityMult",
804  "No list of estimators, falling back to compat");
805  sel->PrintInfo();
806  }
807  return -1;
808  }
809  AliMultEstimator* est = sel->GetEstimator(method);
810  // if (verbose) sel->GetEstimatorList()->ls();
811  if (!est) {
812  if (verbose) {
813  ::Warning("AliForwardUtil::GetCentralityMult",
814  "Unknown estimator: %s", method.Data());
815  sel->GetEstimatorList()->Print();
816  }
817  return -1;
818  }
819  // 198 -> 1: beyond anchor
820  // 199 -> 2: no calib
821  // 200 -> 3: not desired trigger
822  // 201 -> 4: not INEL>0 with tracklets
823  // 202 -> 5: vertex Z not within 10cm
824  // 203 -> 6: tagged as pileup (SPD)
825  // 204 -> 7: inconsistent SPD/tracking vertex
826  // 205 -> 8: rejected by tracklets-vs-clusters
827  Float_t cent = est->GetPercentile();
828  qual = sel->GetEvSelCode();
829  if (qual == AliMultSelectionCuts::kNoCalib) cent = -1;
830  if (verbose)
831  ::Info("AliForwardUtil::GetCentralityMult",
832  "Got centrality %5.1f%% (%d)", /*" - old %5.1f%% (%d)",*/
833  cent, qual/*, old, oldQual*/);
834  return cent;
835 
836 }
837 
838 //____________________________________________________________________
840  const TString& method,
841  Int_t& qual,
842  Bool_t verbose)
843 {
844  if (event.IsA()->InheritsFrom(AliESDEvent::Class()))
845  return GetCentralityCompat(static_cast<const AliESDEvent&>(event),
846  method,qual,verbose);
847  if (event.IsA()->InheritsFrom(AliAODEvent::Class()))
848  return GetCentralityCompat(static_cast<const AliAODEvent&>(event),
849  method,qual,verbose);
850  return -1;
851 }
852 //____________________________________________________________________
854  const TString& method,
855  Int_t& qual,
856  Bool_t verbose)
857 {
858  AliCentrality* centObj = const_cast<AliESDEvent&>(event).GetCentrality();
859  if (!centObj) {
860  if (verbose)
861  ::Warning("AliForwardUtil::GetCentralityCompat",
862  "No centrality object found in ESD");
863  return -1;
864  }
865  Float_t cent = centObj->GetCentralityPercentile(method);
866  qual = centObj->GetQuality();
867  if (verbose)
868  ::Info("AliForwardUtil::GetCentralityCompat<ESD>",
869  "Got centrality %5.1f%% (%d)", cent, qual);
870  if (qual & 0x1) qual = AliMultSelectionCuts::kRejVzCut;
871  if (qual & 0x2) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
872  if (qual & 0x4) qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
873  if (qual & 0x8) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
874  return cent;
875 }
876 //____________________________________________________________________
878  const TString& method,
879  Int_t& qual,
880  Bool_t verbose)
881 {
882  AliAODHeader* hdr = dynamic_cast<AliAODHeader*>(event.GetHeader());
883  if (!hdr) {
884  if (verbose)
885  ::Warning("AliForwardUtil::GetCentralityCompat","Not a standard AOD");
886  return -1;
887  }
888  AliCentrality* cP = hdr->GetCentralityP();
889  if (!cP) {
890  if (verbose)
891  ::Warning("AliForwardUtil::GetCentralityCompat",
892  "No centrality found in AOD");
893  return -1;
894  }
895  Float_t cent = cP->GetCentralityPercentile(method);
896  qual = cP->GetQuality();
897  if (qual & 0x1) qual = AliMultSelectionCuts::kRejVzCut;
898  if (qual & 0x2) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
899  if (qual & 0x4) qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
900  if (qual & 0x8) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
901 
902  if (verbose)
903  ::Info("AliForwardUtil::GetCentralityCompat<AOD>",
904  "Got centrality %5.1f%% (%d)", cent, qual);
905  return cent;
906 }
907 //====================================================================
909 {
910  //
911  // Destructor
912  //
913 }
914 
915 //____________________________________________________________________
916 void
918 {
919  if (fFMD1i) delete fFMD1i;
920  if (fFMD2i) delete fFMD2i;
921  if (fFMD2o) delete fFMD2o;
922  if (fFMD3i) delete fFMD3i;
923  if (fFMD3o) delete fFMD3o;
924  fFMD1i = 0;
925  fFMD2i = 0;
926  fFMD2o = 0;
927  fFMD3i = 0;
928  fFMD3o = 0;
929  TObject::Delete(opt);
930 }
931 
932 //____________________________________________________________________
933 TH2D*
935 {
936  //
937  // Make a histogram
938  //
939  // Parameters:
940  // d Detector
941  // r Ring
942  // etaAxis Eta axis to use
943  //
944  // Return:
945  // Newly allocated histogram
946  //
947  Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
948  TH2D* hist = 0;
949  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
950  hist = new TH2D(Form("FMD%d%c_cache", d, r),
951  Form("FMD%d%c cache", d, r),
952  etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(),
953  ns, 0, TMath::TwoPi());
954  else
955  hist = new TH2D(Form("FMD%d%c_cache", d, r),
956  Form("FMD%d%c cache", d, r),
957  etaAxis.GetNbins(), etaAxis.GetXmin(),
958  etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
959  hist->SetXTitle("#eta");
960  hist->SetYTitle("#phi [radians]");
961  hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
962  hist->Sumw2();
963  hist->SetDirectory(0);
964 
965  return hist;
966 }
967 //____________________________________________________________________
968 void
970 {
971  TAxis* xAxis = hist->GetXaxis();
972  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
973  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
974  else
975  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
976  hist->Rebuild();
977 }
978 
979 
980 //____________________________________________________________________
981 void
983 {
984  //
985  // Initialize the object
986  //
987  // Parameters:
988  // etaAxis Eta axis to use
989  //
990  fFMD1i = Make(1, 'I', etaAxis);
991  fFMD2i = Make(2, 'I', etaAxis);
992  fFMD2o = Make(2, 'O', etaAxis);
993  fFMD3i = Make(3, 'I', etaAxis);
994  fFMD3o = Make(3, 'O', etaAxis);
995 }
996 //____________________________________________________________________
997 void
999 {
1000  //
1001  // Initialize the object
1002  //
1003  // Parameters:
1004  // etaAxis Eta axis to use
1005  //
1006  if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
1007  if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
1008  if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
1009  if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
1010  if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
1011 }
1012 
1013 //____________________________________________________________________
1014 void
1016 {
1017  //
1018  // Clear data
1019  //
1020  // Parameters:
1021  // option Not used
1022  //
1023  if (fFMD1i) { fFMD1i->Reset(option); fFMD1i->ResetBit(kSkipRing); }
1024  if (fFMD2i) { fFMD2i->Reset(option); fFMD2i->ResetBit(kSkipRing); }
1025  if (fFMD2o) { fFMD2o->Reset(option); fFMD2o->ResetBit(kSkipRing); }
1026  if (fFMD3i) { fFMD3i->Reset(option); fFMD3i->ResetBit(kSkipRing); }
1027  if (fFMD3o) { fFMD3o->Reset(option); fFMD3o->ResetBit(kSkipRing); }
1028 }
1029 
1030 //____________________________________________________________________
1031 TH2D*
1033 {
1034  //
1035  // Get the histogram for a particular detector,ring
1036  //
1037  // Parameters:
1038  // d Detector
1039  // r Ring
1040  //
1041  // Return:
1042  // Histogram for detector,ring or nul
1043  //
1044  switch (d) {
1045  case 1: return fFMD1i;
1046  case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1047  case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1048  }
1049  return 0;
1050 }
1051 //====================================================================
1052 TList*
1054 {
1055  //
1056  // Define the outout list in @a d
1057  //
1058  // Parameters:
1059  // d Where to put the output list
1060  //
1061  // Return:
1062  // Newly allocated TList object or null
1063  //
1064  if (!d) return 0;
1065  TList* list = new TList;
1066  list->SetOwner();
1067  list->SetName(fName.Data());
1068  d->Add(list);
1069  return list;
1070 }
1071 //____________________________________________________________________
1072 TList*
1074 {
1075  //
1076  // Get our output list from the container @a d
1077  //
1078  // Parameters:
1079  // d where to get the output list from
1080  //
1081  // Return:
1082  // The found TList or null
1083  //
1084  if (!d) return 0;
1085  TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1086  return list;
1087 }
1088 
1089 //____________________________________________________________________
1090 TH1*
1091 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1092 {
1093  //
1094  // Find a specific histogram in the source list @a d
1095  //
1096  // Parameters:
1097  // d (top)-container
1098  // name Name of histogram
1099  //
1100  // Return:
1101  // Found histogram or null
1102  //
1103  return static_cast<TH1*>(d->FindObject(name));
1104 }
1105 
1106 //====================================================================
1108  const char* format, ...)
1109  : fMsg("")
1110 {
1111  if (lvl < msgLvl) return;
1112  va_list ap;
1113  va_start(ap, format);
1114  Format(fMsg, format, ap);
1115  va_end(ap);
1116  Output(+1, fMsg);
1117 }
1118 //____________________________________________________________________
1120 {
1121  if (fMsg.IsNull()) return;
1122  Output(-1, fMsg);
1123 }
1124 //____________________________________________________________________
1125 void
1127  const char* format, ...)
1128 {
1129  if (lvl < msgLvl) return;
1130  TString msg;
1131  va_list ap;
1132  va_start(ap, format);
1133  Format(msg, format, ap);
1134  va_end(ap);
1135  Output(0, msg);
1136 }
1137 
1138 //____________________________________________________________________
1139 void
1140 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1141 {
1142  static char buf[512];
1143  Int_t n = gROOT->GetDirLevel() + 2;
1144  for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1145  vsnprintf(&(buf[n]), 511-n, format, ap);
1146  buf[511] = '\0';
1147  out = buf;
1148 }
1149 //____________________________________________________________________
1150 void
1152 {
1153  msg[0] = (in > 0 ? '>' : in < 0 ? '<' : '=');
1154  AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1155  if (in > 0) gROOT->IncreaseDirLevel();
1156  else if (in < 0) gROOT->DecreaseDirLevel();
1157 }
1158 
1159 
1160 
1161 //
1162 // EOF
1163 //
return jsonbuilder str().c_str()
static Short_t ParseMagneticField(Float_t field)
void Clear(Option_t *option="")
double Double_t
Definition: External.C:58
static TAxis * MakeFullIpZAxis(Int_t nCenter=20)
static const char * CenterOfMassEnergyString(UShort_t cms)
static void Message(Int_t lvl, Int_t msgLvl, const char *format,...)
static Double_t GetPhiFromStrip(Char_t ring, UShort_t strip, Double_t phi, Double_t xvtx, Double_t yvtx)
static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms)
void ReInit(const TAxis &etaAxis)
static Bool_t GetEtaPhi(UShort_t det, Char_t ring, UShort_t sec, UShort_t str, const TVector3 &ip, Double_t &eta, Double_t &phi)
static Double_t GetSectorZ(UShort_t det, Char_t ring, UShort_t sec)
char Char_t
Definition: External.C:18
TList * list
static void Format(TString &out, const char *format, va_list ap)
static Bool_t GetXYZ(UShort_t det, Char_t ring, UShort_t sec, UShort_t str, const TVector3 &ip, TVector3 &pos)
static Bool_t GetEtaPhiFromStrip(Char_t r, UShort_t strip, Double_t &eta, Double_t &phi, Double_t ipX, Double_t ipY)
static AliAODEvent * GetAODEvent(AliAnalysisTaskSE *task)
#define ALIROOT_SVN_REVISION
static void PrintName(const char *name)
Per-event per bin.
static const char * MagneticFieldString(Short_t field)
static TH2D * Make(UShort_t d, Char_t r, const TAxis &etaAxis)
int Int_t
Definition: External.C:63
static Double_t GetSectorPhi(UShort_t det, Char_t ring, UShort_t sec)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
static Float_t GetCentralityCompat(const AliVEvent &event, const TString &method, Int_t &qual, Bool_t verbose)
Various utilities used in PWGLF/FORWARD.
static Double_t GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx)
static void Output(int in, TString &msg)
Definition: External.C:228
static void GetParameter(TObject *o, UShort_t &value)
unsigned long ULong_t
Definition: External.C:38
static void PrintTask(const TObject &o)
TH2D * Get(UShort_t d, Char_t r) const
short Short_t
Definition: External.C:23
energy
static TObject * MakeParameter(const char *name, UShort_t value)
TH1 * Make(UShort_t d, Char_t r)
Definition: TestEtaPhi.C:3
static Float_t BeamRapidity(Float_t beam, UShort_t z, UShort_t a)
static UShort_t ParseCollisionSystem(const char *sys)
TH1 * GetOutputHist(const TList *d, const char *name) const
static Double_t GetStripR(Char_t ring, UShort_t strip)
static UShort_t CheckForAOD()
static Float_t GetCentrality(const AliVEvent &event, const TString &method, Int_t &qual, Bool_t verbose=false)
void Init(const TAxis &etaAxis)
static ULong_t AliROOTBranch()
static void MakeLogScale(Int_t nBins, Int_t minOrder, Int_t maxOrder, TArrayD &bins)
static void PrintField(const char *name, const char *value,...)
unsigned short UShort_t
Definition: External.C:28
TList * DefineOutputList(TList *d) const
const char Option_t
Definition: External.C:48
void Delete(Option_t *opt="")
static Float_t CenterOfMassEnergy(Float_t beam, UShort_t z1, UShort_t a1, Short_t z2=-1, Short_t a2=-1)
static Float_t CenterOfMassRapidity(UShort_t z1, UShort_t a1, Short_t z2=-1, Short_t a2=-1)
static const char * CollisionSystemString(UShort_t sys)
bool Bool_t
Definition: External.C:53
static Bool_t CheckForTask(const char *clsOrName, Bool_t cls=true)
DebugGuard(Int_t lvl, Int_t msgLvl, const char *format,...)
static ULong_t AliROOTRevision()
static Float_t GetCentralityMult(const AliVEvent &event, const TString &method, Int_t &qual, Bool_t verbose=false)
static void RebinEta(TH2D *hist, const TAxis &etaAxis)
Definition: External.C:196
TList * GetOutputList(const TList *d) const