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