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