AliPhysics  608b256 (608b256)
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("xe-xe") || s.Contains("xexe")) return AliForwardUtil::kXeXe;
145  if (s.Contains("xe54-xe54")) return AliForwardUtil::kXeXe;
146  if (s.Contains("a-a") || s.Contains("aa")) return AliForwardUtil::kPbPb;
148 }
149 //____________________________________________________________________
150 UShort_t
152  Int_t b2a, Int_t b2z,
153  const char* sys)
154 {
155  if (b1a == 129 && b1z == 54 && b2a == 129 && b2z == 54)
156  return AliForwardUtil::kXeXe;
157  if (b1a == 208 && b1z == 82 && b2a == 208 && b2z == 82)
158  return AliForwardUtil::kPbPb;
159  if (b1a == 208 && b1z == 82 && b2a == 1 && b2z == 1)
160  return AliForwardUtil::kPbp;
161  if (b1a == 1 && b1z == 1 && b2a == 208 && b2z == 82)
162  return AliForwardUtil::kPPb;
163  if (b1a == 1 && b1z == 1 && b2a == 1 && b2z == 1)
164  return AliForwardUtil::kPP;
165  return ParseCollisionSystem(sys);
166 }
167 
168 //____________________________________________________________________
169 const char*
171 {
172  //
173  // Get a string representation of the collision system
174  //
175  // Parameters:
176  // sys Collision system
177  // - kPP -> "pp"
178  // - kPbPb -> "PbPb"
179  // - kPPb -> "pPb"
180  // - kPbp -> "Pbp"
181  // - anything else gives "unknown"
182  //
183  // Return:
184  // String representation of the collision system
185  //
186  switch (sys) {
187  case AliForwardUtil::kPP: return "pp";
188  case AliForwardUtil::kPbPb: return "PbPb";
189  case AliForwardUtil::kPPb: return "pPb";
190  case AliForwardUtil::kPbp: return "Pbp";
191  case AliForwardUtil::kXeXe: return "XeXe";
192  }
193  return "unknown";
194 }
195 //____________________________________________________________________
196 Float_t
198 {
199  const Double_t pMass = 9.38271999999999995e-01;
200  const Double_t nMass = 9.39564999999999984e-01;
201  Double_t beamE = z * beam / 2;
202  Double_t beamM = z * pMass + (a - z) * nMass;
203  Double_t beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
204  Double_t beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
205  return beamY;
206 }
207 //____________________________________________________________________
208 Float_t
210  UShort_t z1,
211  UShort_t a1,
212  Short_t z2,
213  Short_t a2)
214 {
215  // Calculate the center of mass energy given target/projectile
216  // mass and charge numbers
217  if (z2 < 0) z2 = z1;
218  if (a2 < 0) a2 = a1;
219  return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
220 }
221 //____________________________________________________________________
222 Float_t
224  UShort_t a1,
225  Short_t z2,
226  Short_t a2)
227 {
228  // Calculate the center of mass rapidity (shift) given target/projectile
229  // mass and charge numbers
230  if (z2 < 0) z2 = z1;
231  if (a2 < 0) a2 = a1;
232  if (z2 == z1 && a2 == a1) return 0;
233  return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
234 }
235 
236 namespace {
237  UShort_t CheckSNN(Float_t energy)
238  {
239  if (TMath::Abs(energy - 900.) < 10) return 900;
240  if (TMath::Abs(energy - 2400.) < 10) return 2400;
241  if (TMath::Abs(energy - 2760.) < 20) return 2760;
242  if (TMath::Abs(energy - 4400.) < 10) return 4400;
243  if (TMath::Abs(energy - 5000.) < 10) return 5000;
244  if (TMath::Abs(energy - 5022.) < 10) return 5023;
245  if (TMath::Abs(energy - 5125.) < 30) return 5100;
246  if (TMath::Abs(energy - 5200.) < 50) return 5200;
247  if (TMath::Abs(energy - 5440.) < 20) return 5440;
248  if (TMath::Abs(energy - 5500.) < 40) return 5500;
249  if (TMath::Abs(energy - 7000.) < 10) return 7000;
250  if (TMath::Abs(energy - 8000.) < 10) return 8000;
251  if (TMath::Abs(energy - 10000.) < 10) return 10000;
252  if (TMath::Abs(energy - 13000.) < 10) return 13000;
253  if (TMath::Abs(energy - 14000.) < 10) return 14000;
254  return 0;
255  }
256 }
257 //____________________________________________________________________
258 UShort_t
260 {
261  //
262  // Parse the center of mass energy given as a float and return known
263  // values as a unsigned integer
264  //
265  // Parameters:
266  // sys Collision system (needed for AA)
267  // beam Center of mass energy * total charge
268  //
269  // Return:
270  // Center of mass energy per nucleon
271  //
272  Float_t energy = beam;
273  // Below no longer needed apparently
274  // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
275  if (sys == AliForwardUtil::kPPb)
276  energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
277  if (sys == AliForwardUtil::kPPb)
278  energy = CenterOfMassEnergy(beam, 54, 129, 1, 1);
279  if (sys == AliForwardUtil::kPbp)
280  energy = CenterOfMassEnergy(beam, 1, 1, 82, 208);
281  else if (sys == AliForwardUtil::kPbPb)
282  energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
283  else if (sys == AliForwardUtil::kXeXe)
284  energy = CenterOfMassEnergy(beam, 54, 129, 54, 129);
285  UShort_t ret = CheckSNN(energy);
286  if (ret > 1) return ret;
287  if (sys == AliForwardUtil::kPbPb ||
288  sys == AliForwardUtil::kPPb ||
289  sys == AliForwardUtil::kPbp) {
290  ret = CheckSNN(beam);
291  }
292  return ret;
293 }
294 //____________________________________________________________________
295 const char*
297 {
298  //
299  // Get a string representation of the center of mass energy per nuclean
300  //
301  // Parameters:
302  // cms Center of mass energy per nucleon
303  //
304  // Return:
305  // String representation of the center of mass energy per nuclean
306  //
307  return Form("%04dGeV", cms);
308 }
309 //____________________________________________________________________
310 Short_t
312 {
313  //
314  // Parse the magnetic field (in kG) as given by a floating point number
315  //
316  // Parameters:
317  // field Magnetic field in kG
318  //
319  // Return:
320  // Short integer value of magnetic field in kG
321  //
322  if (TMath::Abs(v - 5.) < 1 ) return +5;
323  if (TMath::Abs(v - 2.5) < 1 ) return +2;
324  if (TMath::Abs(v - 2.5) < 1 ) return -2;
325  if (TMath::Abs(v + 5.) < 1 ) return -5;
326  if (TMath::Abs(v) < 1) return 0;
327  ::Warning("ParseMagneticfield", "Magnetic field value: %f not known", v);
328  return 999;
329 }
330 //____________________________________________________________________
331 const char*
333 {
334  //
335  // Get a string representation of the magnetic field
336  //
337  // Parameters:
338  // field Magnetic field in kG
339  //
340  // Return:
341  // String representation of the magnetic field
342  //
343  return Form("%01dkG", f);
344 }
345 //_____________________________________________________________________
347 {
348  // Check if AOD is the output event
349  if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
350 
351  AliAODEvent* ret = task->AODEvent();
352  if (ret) return ret;
353 
354  // Check if AOD is the input event
355  ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
356  if (!ret) ::Warning("GetAODEvent", "No AOD event found");
357 
358  return ret;
359 }
360 //_____________________________________________________________________
362 {
363  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
364  if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
365  // ::Info("CheckForAOD", "Found AOD Input handler");
366  return 1;
367  }
368  if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
369  // ::Info("CheckForAOD", "Found AOD Output handler");
370  return 2;
371  }
372 
373  ::Warning("CheckForAOD",
374  "Neither and input nor output AOD handler is specified");
375  return 0;
376 }
377 //_____________________________________________________________________
378 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
379 {
380  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
381  if (!cls) {
382  AliAnalysisTask* t = am->GetTask(clsOrName);
383  if (!t) {
384  ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
385  return false;
386  }
387  ::Info("CheckForTask", "Found task %s", clsOrName);
388  return true;
389  }
390  TClass* dep = gROOT->GetClass(clsOrName);
391  if (!dep) {
392  ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
393  return false;
394  }
395  TIter next(am->GetTasks());
396  TObject* o = 0;
397  while ((o = next())) {
398  if (o->IsA()->InheritsFrom(dep)) {
399  ::Info("CheckForTask", "Found task of class %s: %s",
400  clsOrName, o->GetName());
401  return true;
402  }
403  }
404  ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
405  return false;
406 }
407 
408 //_____________________________________________________________________
410 {
411  TParameter<int>* ret = new TParameter<int>(name, value);
412  ret->SetMergeMode('f');
413  ret->SetUniqueID(value);
414  return ret;
415 }
416 //_____________________________________________________________________
417 TObject* AliForwardUtil::MakeParameter(const char* name, Int_t value)
418 {
419  TParameter<int>* ret = new TParameter<int>(name, value);
420  ret->SetMergeMode('f');
421  ret->SetUniqueID(value);
422  return ret;
423 }
424 //_____________________________________________________________________
426 {
427  TParameter<Long_t>* ret = new TParameter<Long_t>(name, value);
428  ret->SetMergeMode('f');
429  ret->SetUniqueID(value);
430  return ret;
431 }
432 //_____________________________________________________________________
434 {
435  TParameter<double>* ret = new TParameter<double>(name, value);
436  // Float_t v = value;
437  // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
438  ret->SetMergeMode('f');
439  // ret->SetUniqueID(*tmp);
440  return ret;
441 }
442 //_____________________________________________________________________
444 {
445  TParameter<bool>* ret = new TParameter<bool>(name, value);
446  ret->SetMergeMode('f');
447  ret->SetUniqueID(value);
448  return ret;
449 }
450 
451 //_____________________________________________________________________
453 {
454  if (!o) return;
455  TParameter<int>* p = static_cast<TParameter<int>*>(o);
456  if (p->TestBit(BIT(19)))
457  value = p->GetVal();
458  else
459  value = o->GetUniqueID();
460 }
461 //_____________________________________________________________________
463 {
464  if (!o) return;
465  TParameter<int>* p = static_cast<TParameter<int>*>(o);
466  if (p->TestBit(BIT(19)))
467  value = p->GetVal();
468  else
469  value = o->GetUniqueID();
470 }
471 //_____________________________________________________________________
473 {
474  if (!o) return;
475  TParameter<Long_t>* p = static_cast<TParameter<Long_t>*>(o);
476  if (p->TestBit(BIT(19)))
477  value = p->GetVal();
478  else
479  value = o->GetUniqueID();
480 }
481 //_____________________________________________________________________
483 {
484  if (!o) return;
485  TParameter<double>* p = static_cast<TParameter<double>*>(o);
486  if (p->TestBit(BIT(19)))
487  value = p->GetVal(); // o->GetUniqueID();
488  else {
489  UInt_t i = o->GetUniqueID();
490  Float_t v = *reinterpret_cast<Float_t*>(&i);
491  value = v;
492  }
493 }
494 //_____________________________________________________________________
496 {
497  if (!o) return;
498  TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
499  if (p->TestBit(BIT(19)))
500  value = p->GetVal(); // o->GetUniqueID();
501  else
502  value = o->GetUniqueID();
503 }
504 
505 //_____________________________________________________________________
507 {
508  // Get max R of ring
509  //
510  // New implementation has only one branch
511  const Double_t minR[] = { 4.5213, 15.4 };
512  const Double_t maxR[] = { 17.2, 28.0 };
513  const Int_t nStr[] = { 512, 256 };
514 
515  Int_t q = (ring == 'I' || ring == 'i') ? 0 : 1;
516  Double_t rad = maxR[q] - minR[q];
517  Double_t segment = rad / nStr[q];
518  Double_t r = minR[q] + segment*strip;
519 
520  return r;
521 }
522 
523 //_____________________________________________________________________
525 {
526  Int_t hybrid = sec / 2;
527  Int_t q = (ring == 'I' || ring == 'i') ? 0 : 1;
528  Int_t r = q == 0 ? 1 : 0;
529  const Double_t zs[][2] = { { 320.266+1.5, kInvalidValue },
530  { 83.666, 74.966+.5 },
531  { -63.066+.5, -74.966 } };
532  if (det > 3 || zs[det-1][q] == kInvalidValue) {
533  ::Warning("GetSectorZ", "Unknown sub-detector FMD%d%c", det, ring);
534  return kInvalidValue;
535  }
536 
537  Double_t z = zs[det-1][q];
538  switch (det) {
539  case 1: if ((hybrid % 2) == 1) z -= .5; break;
540  case 2: if ((hybrid % 2) == r) z -= .5; break;
541  case 3: if ((hybrid % 2) == q) z -= .5; break;
542  }
543 
544  return z;
545 }
546 
547 //_____________________________________________________________________
549 {
550  UShort_t nSec = (ring == 'I' || ring == 'i') ? 20 : 40;
551  Double_t base = float(sec+.5) / nSec * TMath::TwoPi();
552  switch (d) {
553  case 1: base += TMath::Pi()/2; break;
554  case 2: break;
555  case 3: base = TMath::Pi() - base; break;
556  default:
557  ::Warning("GetSectorPhi", "Unknown detector %d", d);
558  return kInvalidValue;
559  }
560  if (base < 0) base += TMath::TwoPi();
561  if (base > TMath::TwoPi()) base -= TMath::TwoPi();
562  return base;
563 }
564 
565 //_____________________________________________________________________
567  UShort_t sec, UShort_t strip,
568  Double_t zvtx)
569 {
570  // Calculate eta from strip with vertex (redundant with
571  // AliESDFMD::Eta but support displaced vertices)
572  //
573  // Slightly more optimized version that uses less branching
574 
575  // Get R of the strip
576  Double_t r = GetStripR(ring, strip);
577  Double_t z = GetSectorZ(det, ring, sec);
578  Double_t theta = TMath::ATan2(r,z-zvtx);
579  Double_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
580 
581  return eta;
582 }
583 
584 //_____________________________________________________________________
586  UShort_t sec, UShort_t strip,
587  const TVector3& ip,
588  TVector3& pos)
589 {
590  Double_t rD = GetStripR(ring, strip);
591  Double_t phiD = GetSectorPhi(det,ring, sec);
592  Double_t zD = GetSectorZ(det, ring, sec);
593  if (phiD == kInvalidValue || zD == kInvalidValue) {
595  return false;
596  }
597  Double_t xD = rD*TMath::Cos(phiD);
598  Double_t yD = rD*TMath::Sin(phiD);
599  Double_t iX = ip.X(); if (iX > 100) iX = 0; // No X
600  Double_t iY = ip.Y(); if (iY > 100) iY = 0; // No Y
601  Double_t dX = xD-iX;
602  Double_t dY = yD-iY;
603  Double_t dZ = zD-ip.Z();
604  pos.SetXYZ(dX, dY, dZ);
605  return true;
606 }
607 //_____________________________________________________________________
609  UShort_t sec, UShort_t strip,
610  const TVector3& ip,
611  Double_t& eta, Double_t& phi)
612 {
613  TVector3 pos;
614  if (!GetXYZ(det, ring, sec, strip, ip, pos)) {
615  ::Warning("GetEtaPhi", "Invalid position for FMD%d%c[%2d,%3d]=(%f,%f,%f)",
616  det, ring, sec, strip, pos.X(), pos.Y(), pos.Z());
617  eta = kInvalidValue;
618  phi = kInvalidValue;
619  return false;
620  }
621  Double_t r = TMath::Sqrt(TMath::Power(pos.X(),2)+
622  TMath::Power(pos.Y(),2));
623  Double_t theta = TMath::ATan2(r, pos.Z());
624  Double_t tant = TMath::Tan(theta/2);
625  if (TMath::Abs(theta) < 1e-9) {
626  ::Warning("GetEtaPhi","tan(theta/2)=%f very small", tant);
627  eta = kInvalidValue;
628  phi = kInvalidValue;
629  return false;
630  }
631  phi = TMath::ATan2(pos.Y(), pos.X());
632  eta = -TMath::Log(tant);
633  if (phi < 0) phi += TMath::TwoPi();
634  if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
635 
636  return true;
637 }
638 
639 //_____________________________________________________________________
641  UShort_t strip,
642  Double_t& eta, Double_t& phi ,
643  Double_t ipX, Double_t ipY)
644 {
645  Double_t rs = GetStripR(r, strip);
646  Double_t sx = rs*TMath::Cos(phi);
647  Double_t sy = rs*TMath::Sin(phi);
648  Double_t dx = sx-ipX;
649  Double_t dy = sy-ipY;
650  Double_t rv = TMath::Sqrt(TMath::Power(dx,2) + TMath::Power(dy,2));
651  Double_t the = 2*TMath::ATan(TMath::Exp(-eta));
652  Double_t tth = TMath::Tan(the);
653  if (TMath::Abs(tth) < 1e-9) {
654  ::Warning("GetEtaPhiFromStrip",
655  "eta=%f -> theta=%f tan(theta)=%f invalid (no change)",
656  eta, the, tth);
657  return false;
658  }
659  Double_t z = rs / tth;
660  // Printf("IP(x,y)=%f,%f S(x,y)=%f,%f D(x,y)=%f,%f R=%f theta=%f "
661  // "tan(theta)=%f z=%f",
662  // ipX, ipY, sx, sy, dx, dy, rv, the, TMath::Tan(the), z);
663  eta = -TMath::Log(TMath::Tan(TMath::ATan2(rv,z)/2));
664  phi = TMath::ATan2(dy,dx);
665  if (phi < 0) phi += TMath::TwoPi();
666 
667  return true;
668 }
669 
670 
671 //_____________________________________________________________________
673  Double_t phi,
674  Double_t xvtx, Double_t yvtx)
675 {
676  // Calculate eta from strip with vertex (redundant with
677  // AliESDFMD::Eta but support displaced vertices)
678 
679  // Unknown x,y -> no change
680  if (yvtx > 999 || xvtx > 999) return phi;
681 
682  //Get max R of ring
683  Double_t r = GetStripR(ring, strip);
684  Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
685  Double_t pha = (TMath::Abs(yvtx) < 1e-12 ? 0 : TMath::ATan2(xvtx, yvtx));
686  Double_t cha = amp * TMath::Cos(phi+pha);
687  phi += cha;
688  if (phi < 0) phi += TMath::TwoPi();
689  if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
690  return phi;
691 }
692 //====================================================================
693 TAxis*
695 {
696  TArrayD bins;
697  MakeFullIpZAxis(nCenter, bins);
698  TAxis* a = new TAxis(bins.GetSize()-1,bins.GetArray());
699  return a;
700 }
701 //____________________________________________________________________
702 void
704 {
705  // Custom vertex axis that will include satellite vertices
706  // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10
707  // Nominal vertices are usually in -10 to 10 and we should have
708  // 10 bins in that range. That gives us a total of
709  //
710  // 10+10+10=30 bins
711  //
712  // or 31 bin boundaries
713  if (nCenter % 2 == 1)
714  // Number of central bins is odd - make it even
715  nCenter--;
716  const Double_t mCenter = 20;
717  const Int_t nSat = 10;
718  const Int_t nBins = 2*nSat + nCenter;
719  const Int_t mBin = nBins / 2;
720  Double_t dCenter = 2*mCenter / nCenter;
721  bins.Set(nBins+1);
722  bins[mBin] = 0;
723  for (Int_t i = 1; i <= nCenter/2; i++) {
724  // Assign from the middle out
725  Double_t v = i * dCenter;
726  // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-i,mBin+i);
727  bins[mBin-i] = -v;
728  bins[mBin+i] = +v;
729  }
730  for (Int_t i = 1; i <= nSat; i++) {
731  Double_t v = (i+.5) * 37.5;
732  Int_t o = nCenter/2+i;
733  // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-o,mBin+o);
734  bins[mBin-o] = -v;
735  bins[mBin+o] = +v;
736  }
737 }
738 //____________________________________________________________________
739 void
741  Int_t minOrder,
742  Int_t maxOrder,
743  TArrayD& bins)
744 {
745  Double_t dO = Double_t(maxOrder-minOrder) / nBins;
746  bins.Set(nBins+1);
747  for (Int_t i = 0; i <= nBins; i++) bins[i] = TMath::Power(10, i * dO+minOrder);
748 }
749 
750 //____________________________________________________________________
751 void
753 {
754  Int_t ind = gROOT->GetDirLevel();
755  if (ind > 0)
756  // Print indention
757  std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
758 
759  TString t = TString::Format("%s %s", o.GetName(), o.ClassName());
760  const Int_t maxN = 75;
761  std::cout << "--- " << t << " " << std::setfill('-')
762  << std::setw(maxN-ind-5-t.Length()) << "-" << std::endl;
763 }
764 //____________________________________________________________________
765 void
766 AliForwardUtil::PrintName(const char* name)
767 {
768  Int_t ind = gROOT->GetDirLevel();
769  if (ind > 0)
770  // Print indention
771  std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
772 
773  // Now print field name
774  const Int_t maxN = 29;
775  Int_t width = maxN - ind;
776  TString n(name);
777  if (n.Length() > width-1) {
778  // Truncate the string, and put in "..."
779  n.Remove(width-4);
780  n.Append("...");
781  }
782  n.Append(":");
783  std::cout << std::setfill(' ') << std::left << std::setw(width)
784  << n << std::right << std::flush;
785 }
786 //____________________________________________________________________
787 void
788 AliForwardUtil::PrintField(const char* name, const char* value, ...)
789 {
790  PrintName(name);
791 
792  // Now format the field value
793  va_list ap;
794  va_start(ap, value);
795  static char buf[512];
796  vsnprintf(buf, 511, value, ap);
797  buf[511] = '\0';
798  va_end(ap);
799 
800  std::cout << buf << std::endl;
801 }
802 
803 
804 //====================================================================
805 Float_t AliForwardUtil::GetCentrality(const AliVEvent& event,
806  const TString& method,
807  Int_t& qual,
808  Bool_t verbose)
809 {
810  qual = 0xFFFF;
811  Float_t cent = GetCentralityMult(event,method,qual,verbose);
812  if (qual < 0xFFF) return cent;
813 
814  // No selection object found, try compat
815  return GetCentralityCompat(event,method,qual,verbose);
816 }
817 //____________________________________________________________________
819  const TString& method,
820  Int_t& qual,
821  Bool_t verbose)
822 {
823  TObject* o = event.FindListObject("MultSelection");
824  if (!o) {
825  if (verbose)
826  ::Warning("AliForwardUtil::GetCentralityMult",
827  "No MultSelection object found in event");
828  return -1;
829  }
830  AliMultSelection* sel = static_cast<AliMultSelection*>(o);
831  if (!sel->GetEstimatorList() ||
832  sel->GetEstimatorList()->GetEntries() <= 0){
833  if (verbose) {
834  ::Warning("AliForwardUtil::GetCentralityMult",
835  "No list of estimators, falling back to compat");
836  sel->PrintInfo();
837  }
838  return -1;
839  }
840  AliMultEstimator* est = sel->GetEstimator(method);
841  // if (verbose) sel->GetEstimatorList()->ls();
842  if (!est) {
843  if (verbose) {
844  ::Warning("AliForwardUtil::GetCentralityMult",
845  "Unknown estimator: %s", method.Data());
846  sel->GetEstimatorList()->Print();
847  }
848  return -1;
849  }
850  // 198 -> 1: beyond anchor
851  // 199 -> 2: no calib
852  // 200 -> 3: not desired trigger
853  // 201 -> 4: not INEL>0 with tracklets
854  // 202 -> 5: vertex Z not within 10cm
855  // 203 -> 6: tagged as pileup (SPD)
856  // 204 -> 7: inconsistent SPD/tracking vertex
857  // 205 -> 8: rejected by tracklets-vs-clusters
858  Float_t cent = est->GetPercentile();
859  qual = sel->GetEvSelCode();
860  if (qual == AliMultSelectionCuts::kNoCalib) cent = -1;
861  if (verbose)
862  ::Info("AliForwardUtil::GetCentralityMult",
863  "Got centrality %5.1f%% (%d)", /*" - old %5.1f%% (%d)",*/
864  cent, qual/*, old, oldQual*/);
865  return cent;
866 
867 }
868 
869 //____________________________________________________________________
871  const TString& method,
872  Int_t& qual,
873  Bool_t verbose)
874 {
875  if (event.IsA()->InheritsFrom(AliESDEvent::Class()))
876  return GetCentralityCompat(static_cast<const AliESDEvent&>(event),
877  method,qual,verbose);
878  if (event.IsA()->InheritsFrom(AliAODEvent::Class()))
879  return GetCentralityCompat(static_cast<const AliAODEvent&>(event),
880  method,qual,verbose);
881  return -1;
882 }
883 //____________________________________________________________________
885  const TString& method,
886  Int_t& qual,
887  Bool_t verbose)
888 {
889  AliCentrality* centObj = const_cast<AliESDEvent&>(event).GetCentrality();
890  if (!centObj) {
891  if (verbose)
892  ::Warning("AliForwardUtil::GetCentralityCompat",
893  "No centrality object found in ESD");
894  return -1;
895  }
896  Float_t cent = centObj->GetCentralityPercentileUnchecked(method);
897  qual = centObj->GetQuality();
898  if (verbose)
899  ::Info("AliForwardUtil::GetCentralityCompat<ESD>",
900  "Got centrality %5.1f%% (%d)", cent, qual);
901  if (qual & 0x1)
902  qual = AliMultSelectionCuts::kRejVzCut;
903  if (qual & 0x2)
904  qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
905  if (qual & 0x4)
906  qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
907  if (qual & 0x8)
908  qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
909  return cent;
910 }
911 //____________________________________________________________________
913  const TString& method,
914  Int_t& qual,
915  Bool_t verbose)
916 {
917  AliAODHeader* hdr = dynamic_cast<AliAODHeader*>(event.GetHeader());
918  if (!hdr) {
919  if (verbose)
920  ::Warning("AliForwardUtil::GetCentralityCompat","Not a standard AOD");
921  return -1;
922  }
923  AliCentrality* cP = hdr->GetCentralityP();
924  if (!cP) {
925  if (verbose)
926  ::Warning("AliForwardUtil::GetCentralityCompat",
927  "No centrality found in AOD");
928  return -1;
929  }
930  Float_t cent = cP->GetCentralityPercentile(method);
931  qual = cP->GetQuality();
932  if (qual & 0x1) qual = AliMultSelectionCuts::kRejVzCut;
933  if (qual & 0x2) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
934  if (qual & 0x4) qual = AliMultSelectionCuts::kRejConsistencySPDandTrackVertices;
935  if (qual & 0x8) qual = AliMultSelectionCuts::kRejTrackletsVsClusters;
936 
937  if (verbose)
938  ::Info("AliForwardUtil::GetCentralityCompat<AOD>",
939  "Got centrality %5.1f%% (%d)", cent, qual);
940  return cent;
941 }
942 //====================================================================
944 {
945  //
946  // Destructor
947  //
948 }
949 
950 //____________________________________________________________________
951 void
953 {
954  if (fFMD1i) delete fFMD1i;
955  if (fFMD2i) delete fFMD2i;
956  if (fFMD2o) delete fFMD2o;
957  if (fFMD3i) delete fFMD3i;
958  if (fFMD3o) delete fFMD3o;
959  fFMD1i = 0;
960  fFMD2i = 0;
961  fFMD2o = 0;
962  fFMD3i = 0;
963  fFMD3o = 0;
964  TObject::Delete(opt);
965 }
966 
967 //____________________________________________________________________
968 TH2D*
970 {
971  //
972  // Make a histogram
973  //
974  // Parameters:
975  // d Detector
976  // r Ring
977  // etaAxis Eta axis to use
978  //
979  // Return:
980  // Newly allocated histogram
981  //
982  Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
983  TH2D* hist = 0;
984  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
985  hist = new TH2D(Form("FMD%d%c_cache", d, r),
986  Form("FMD%d%c cache", d, r),
987  etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(),
988  ns, 0, TMath::TwoPi());
989  else
990  hist = new TH2D(Form("FMD%d%c_cache", d, r),
991  Form("FMD%d%c cache", d, r),
992  etaAxis.GetNbins(), etaAxis.GetXmin(),
993  etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
994  hist->SetXTitle("#eta");
995  hist->SetYTitle("#phi [radians]");
996  hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
997  hist->Sumw2();
998  hist->SetDirectory(0);
999 
1000  return hist;
1001 }
1002 //____________________________________________________________________
1003 void
1005 {
1006  TAxis* xAxis = hist->GetXaxis();
1007  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1008  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
1009  else
1010  xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
1011  hist->Rebuild();
1012 }
1013 
1014 
1015 //____________________________________________________________________
1016 void
1018 {
1019  //
1020  // Initialize the object
1021  //
1022  // Parameters:
1023  // etaAxis Eta axis to use
1024  //
1025  fFMD1i = Make(1, 'I', etaAxis);
1026  fFMD2i = Make(2, 'I', etaAxis);
1027  fFMD2o = Make(2, 'O', etaAxis);
1028  fFMD3i = Make(3, 'I', etaAxis);
1029  fFMD3o = Make(3, 'O', etaAxis);
1030 }
1031 //____________________________________________________________________
1032 void
1034 {
1035  //
1036  // Initialize the object
1037  //
1038  // Parameters:
1039  // etaAxis Eta axis to use
1040  //
1041  if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
1042  if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
1043  if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
1044  if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
1045  if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
1046 }
1047 
1048 //____________________________________________________________________
1049 void
1051 {
1052  //
1053  // Clear data
1054  //
1055  // Parameters:
1056  // option Not used
1057  //
1058  if (fFMD1i) { fFMD1i->Reset(option); fFMD1i->ResetBit(kSkipRing); }
1059  if (fFMD2i) { fFMD2i->Reset(option); fFMD2i->ResetBit(kSkipRing); }
1060  if (fFMD2o) { fFMD2o->Reset(option); fFMD2o->ResetBit(kSkipRing); }
1061  if (fFMD3i) { fFMD3i->Reset(option); fFMD3i->ResetBit(kSkipRing); }
1062  if (fFMD3o) { fFMD3o->Reset(option); fFMD3o->ResetBit(kSkipRing); }
1063 }
1064 
1065 //____________________________________________________________________
1066 TH2D*
1068 {
1069  //
1070  // Get the histogram for a particular detector,ring
1071  //
1072  // Parameters:
1073  // d Detector
1074  // r Ring
1075  //
1076  // Return:
1077  // Histogram for detector,ring or nul
1078  //
1079  switch (d) {
1080  case 1: return fFMD1i;
1081  case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1082  case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1083  }
1084  return 0;
1085 }
1086 //====================================================================
1087 TList*
1089 {
1090  //
1091  // Define the outout list in @a d
1092  //
1093  // Parameters:
1094  // d Where to put the output list
1095  //
1096  // Return:
1097  // Newly allocated TList object or null
1098  //
1099  if (!d) return 0;
1100  TList* list = new TList;
1101  list->SetOwner();
1102  list->SetName(fName.Data());
1103  d->Add(list);
1104  return list;
1105 }
1106 //____________________________________________________________________
1107 TList*
1109 {
1110  //
1111  // Get our output list from the container @a d
1112  //
1113  // Parameters:
1114  // d where to get the output list from
1115  //
1116  // Return:
1117  // The found TList or null
1118  //
1119  if (!d) return 0;
1120  TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1121  return list;
1122 }
1123 
1124 //____________________________________________________________________
1125 TH1*
1126 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1127 {
1128  //
1129  // Find a specific histogram in the source list @a d
1130  //
1131  // Parameters:
1132  // d (top)-container
1133  // name Name of histogram
1134  //
1135  // Return:
1136  // Found histogram or null
1137  //
1138  return static_cast<TH1*>(d->FindObject(name));
1139 }
1140 
1141 
1142 //====================================================================
1144  const char* format, ...)
1145  : fMsg("")
1146 {
1147  if (lvl < msgLvl) return;
1148  va_list ap;
1149  va_start(ap, format);
1150  Format(fMsg, format, ap);
1151  va_end(ap);
1152  Output(+1, fMsg);
1153 }
1154 //____________________________________________________________________
1156 {
1157  if (fMsg.IsNull()) return;
1158  Output(-1, fMsg);
1159 }
1160 //____________________________________________________________________
1161 void
1163  const char* format, ...)
1164 {
1165  if (lvl < msgLvl) return;
1166  TString msg;
1167  va_list ap;
1168  va_start(ap, format);
1169  Format(msg, format, ap);
1170  va_end(ap);
1171  Output(0, msg);
1172 }
1173 
1174 //____________________________________________________________________
1175 void
1176 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1177 {
1178  static char buf[512];
1179  Int_t n = gROOT->GetDirLevel() + 2;
1180  for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1181  vsnprintf(&(buf[n]), 511-n, format, ap);
1182  buf[511] = '\0';
1183  out = buf;
1184 }
1185 //____________________________________________________________________
1186 void
1188 {
1189  msg[0] = (in > 0 ? '>' : in < 0 ? '<' : '=');
1190  AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1191  if (in > 0) gROOT->IncreaseDirLevel();
1192  else if (in < 0) gROOT->DecreaseDirLevel();
1193 }
1194 //====================================================================
1196  : save(gErrorIgnoreLevel)
1197 {
1198  gErrorIgnoreLevel = lvl;
1199  AliLog::SetModuleDebugLevel("ROOT", lvl);
1200  AliLog::SetGlobalDebugLevel(lvl);
1201 }
1202 //____________________________________________________________________
1204 {
1205  gErrorIgnoreLevel = save;
1206  AliLog::SetModuleDebugLevel("ROOT", save);
1207  AliLog::SetGlobalDebugLevel(save);
1208 }
1209 
1210 //
1211 // EOF
1212 //
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)
TString format
file names tag, basically the trigger and calorimeter combination
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:45
static Double_t GetSectorZ(UShort_t det, Char_t ring, UShort_t sec)
char Char_t
Definition: External.C:18
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