AliPhysics  7f2a7c4 (7f2a7c4)
 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->GetCentralityPercentileUnchecked(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)
871  qual = AliMultSelectionCuts::kRejVzCut;
872  if (qual & 0x2)
873  qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
874  if (qual & 0x4)
875  qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
876  if (qual & 0x8)
877  qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
878  return cent;
879 }
880 //____________________________________________________________________
882  const TString& method,
883  Int_t& qual,
884  Bool_t verbose)
885 {
886  AliAODHeader* hdr = dynamic_cast<AliAODHeader*>(event.GetHeader());
887  if (!hdr) {
888  if (verbose)
889  ::Warning("AliForwardUtil::GetCentralityCompat","Not a standard AOD");
890  return -1;
891  }
892  AliCentrality* cP = hdr->GetCentralityP();
893  if (!cP) {
894  if (verbose)
895  ::Warning("AliForwardUtil::GetCentralityCompat",
896  "No centrality found in AOD");
897  return -1;
898  }
899  Float_t cent = cP->GetCentralityPercentile(method);
900  qual = cP->GetQuality();
901  if (qual & 0x1) qual = AliMultSelectionCuts::kRejVzCut;
902  if (qual & 0x2) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
903  if (qual & 0x4) qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
904  if (qual & 0x8) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
905 
906  if (verbose)
907  ::Info("AliForwardUtil::GetCentralityCompat<AOD>",
908  "Got centrality %5.1f%% (%d)", cent, qual);
909  return cent;
910 }
911 //====================================================================
913 {
914  //
915  // Destructor
916  //
917 }
918 
919 //____________________________________________________________________
920 void
922 {
923  if (fFMD1i) delete fFMD1i;
924  if (fFMD2i) delete fFMD2i;
925  if (fFMD2o) delete fFMD2o;
926  if (fFMD3i) delete fFMD3i;
927  if (fFMD3o) delete fFMD3o;
928  fFMD1i = 0;
929  fFMD2i = 0;
930  fFMD2o = 0;
931  fFMD3i = 0;
932  fFMD3o = 0;
933  TObject::Delete(opt);
934 }
935 
936 //____________________________________________________________________
937 TH2D*
939 {
940  //
941  // Make a histogram
942  //
943  // Parameters:
944  // d Detector
945  // r Ring
946  // etaAxis Eta axis to use
947  //
948  // Return:
949  // Newly allocated histogram
950  //
951  Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
952  TH2D* hist = 0;
953  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
954  hist = new TH2D(Form("FMD%d%c_cache", d, r),
955  Form("FMD%d%c cache", d, r),
956  etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(),
957  ns, 0, TMath::TwoPi());
958  else
959  hist = new TH2D(Form("FMD%d%c_cache", d, r),
960  Form("FMD%d%c cache", d, r),
961  etaAxis.GetNbins(), etaAxis.GetXmin(),
962  etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
963  hist->SetXTitle("#eta");
964  hist->SetYTitle("#phi [radians]");
965  hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
966  hist->Sumw2();
967  hist->SetDirectory(0);
968 
969  return hist;
970 }
971 //____________________________________________________________________
972 void
974 {
975  TAxis* xAxis = hist->GetXaxis();
976  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
977  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
978  else
979  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
980  hist->Rebuild();
981 }
982 
983 
984 //____________________________________________________________________
985 void
987 {
988  //
989  // Initialize the object
990  //
991  // Parameters:
992  // etaAxis Eta axis to use
993  //
994  fFMD1i = Make(1, 'I', etaAxis);
995  fFMD2i = Make(2, 'I', etaAxis);
996  fFMD2o = Make(2, 'O', etaAxis);
997  fFMD3i = Make(3, 'I', etaAxis);
998  fFMD3o = Make(3, 'O', etaAxis);
999 }
1000 //____________________________________________________________________
1001 void
1003 {
1004  //
1005  // Initialize the object
1006  //
1007  // Parameters:
1008  // etaAxis Eta axis to use
1009  //
1010  if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
1011  if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
1012  if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
1013  if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
1014  if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
1015 }
1016 
1017 //____________________________________________________________________
1018 void
1020 {
1021  //
1022  // Clear data
1023  //
1024  // Parameters:
1025  // option Not used
1026  //
1027  if (fFMD1i) { fFMD1i->Reset(option); fFMD1i->ResetBit(kSkipRing); }
1028  if (fFMD2i) { fFMD2i->Reset(option); fFMD2i->ResetBit(kSkipRing); }
1029  if (fFMD2o) { fFMD2o->Reset(option); fFMD2o->ResetBit(kSkipRing); }
1030  if (fFMD3i) { fFMD3i->Reset(option); fFMD3i->ResetBit(kSkipRing); }
1031  if (fFMD3o) { fFMD3o->Reset(option); fFMD3o->ResetBit(kSkipRing); }
1032 }
1033 
1034 //____________________________________________________________________
1035 TH2D*
1037 {
1038  //
1039  // Get the histogram for a particular detector,ring
1040  //
1041  // Parameters:
1042  // d Detector
1043  // r Ring
1044  //
1045  // Return:
1046  // Histogram for detector,ring or nul
1047  //
1048  switch (d) {
1049  case 1: return fFMD1i;
1050  case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1051  case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1052  }
1053  return 0;
1054 }
1055 //====================================================================
1056 TList*
1058 {
1059  //
1060  // Define the outout list in @a d
1061  //
1062  // Parameters:
1063  // d Where to put the output list
1064  //
1065  // Return:
1066  // Newly allocated TList object or null
1067  //
1068  if (!d) return 0;
1069  TList* list = new TList;
1070  list->SetOwner();
1071  list->SetName(fName.Data());
1072  d->Add(list);
1073  return list;
1074 }
1075 //____________________________________________________________________
1076 TList*
1078 {
1079  //
1080  // Get our output list from the container @a d
1081  //
1082  // Parameters:
1083  // d where to get the output list from
1084  //
1085  // Return:
1086  // The found TList or null
1087  //
1088  if (!d) return 0;
1089  TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1090  return list;
1091 }
1092 
1093 //____________________________________________________________________
1094 TH1*
1095 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1096 {
1097  //
1098  // Find a specific histogram in the source list @a d
1099  //
1100  // Parameters:
1101  // d (top)-container
1102  // name Name of histogram
1103  //
1104  // Return:
1105  // Found histogram or null
1106  //
1107  return static_cast<TH1*>(d->FindObject(name));
1108 }
1109 
1110 
1111 //====================================================================
1113  const char* format, ...)
1114  : fMsg("")
1115 {
1116  if (lvl < msgLvl) return;
1117  va_list ap;
1118  va_start(ap, format);
1119  Format(fMsg, format, ap);
1120  va_end(ap);
1121  Output(+1, fMsg);
1122 }
1123 //____________________________________________________________________
1125 {
1126  if (fMsg.IsNull()) return;
1127  Output(-1, fMsg);
1128 }
1129 //____________________________________________________________________
1130 void
1132  const char* format, ...)
1133 {
1134  if (lvl < msgLvl) return;
1135  TString msg;
1136  va_list ap;
1137  va_start(ap, format);
1138  Format(msg, format, ap);
1139  va_end(ap);
1140  Output(0, msg);
1141 }
1142 
1143 //____________________________________________________________________
1144 void
1145 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1146 {
1147  static char buf[512];
1148  Int_t n = gROOT->GetDirLevel() + 2;
1149  for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1150  vsnprintf(&(buf[n]), 511-n, format, ap);
1151  buf[511] = '\0';
1152  out = buf;
1153 }
1154 //____________________________________________________________________
1155 void
1157 {
1158  msg[0] = (in > 0 ? '>' : in < 0 ? '<' : '=');
1159  AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1160  if (in > 0) gROOT->IncreaseDirLevel();
1161  else if (in < 0) gROOT->DecreaseDirLevel();
1162 }
1163 //====================================================================
1165  : save(gErrorIgnoreLevel)
1166 {
1167  gErrorIgnoreLevel = lvl;
1168  AliLog::SetModuleDebugLevel("ROOT", lvl);
1169  AliLog::SetGlobalDebugLevel(lvl);
1170 }
1171 //____________________________________________________________________
1173 {
1174  gErrorIgnoreLevel = save;
1175  AliLog::SetModuleDebugLevel("ROOT", save);
1176  AliLog::SetGlobalDebugLevel(save);
1177 }
1178 
1179 //
1180 // EOF
1181 //
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)
energy
Definition: HFPtSpectrum.C:44
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.
Int_t method
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
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