Fix some documentation issues
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardUtil.cxx
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 <TParameter.h>
20 #include <TH2D.h>
21 #include <TH1I.h>
22 #include <TF1.h>
23 #include <TFitResult.h>
24 #include <TMath.h>
25 #include <TError.h>
26 #include <TROOT.h>
27 #define FIT_OPTIONS "RNQS"
28
29 //====================================================================
30 ULong_t AliForwardUtil::AliROOTRevision()
31 {
32 #ifdef ALIROOT_SVN_REVISION
33   return ALIROOT_SVN_REVISION;
34 #else 
35   return 0;
36 #endif
37 }
38 //____________________________________________________________________
39 ULong_t AliForwardUtil::AliROOTBranch()
40 {
41   // Do something here when we switch to git - sigh!
42 #ifdef ALIROOT_SVN_BRANCH
43   static ULong_t ret = 0;
44   if (ret != 0) return ret;
45   
46   TString str(ALIROOT_SVN_BRANCH);
47   if (str[0] == 'v') str.Remove(0,1);
48   if (str.EqualTo("trunk")) return ret = 0xFFFFFFFF;
49
50   TObjArray*   tokens = str.Tokenize("-");
51   TObjString*  pMajor = static_cast<TObjString*>(tokens->At(0));
52   TObjString*  pMinor = static_cast<TObjString*>(tokens->At(1));
53   TObjString*  pRelea = static_cast<TObjString*>(tokens->At(2));
54   TObjString* pAn     = (tokens->GetEntries() > 3 ? 
55     static_cast<TObjString*>(tokens->At(3)) : 0);
56   TString sMajor = pMajor->String().Strip(TString::kLeading, '0');
57   TString sMinor = pMinor->String().Strip(TString::kLeading, '0');
58   TString sRelea = pRelea->String().Strip(TString::kLeading, '0');
59   
60   ret = (((sMajor.Atoi() & 0xFF) << 12) |
61     ((sMinor.Atoi() & 0xFF) <<  8) |
62     ((sRelea.Atoi() & 0xFF) <<  4) |
63     (pAn ? 0xAA : 0));
64   
65   return ret;
66 #else 
67   return 0;
68 #endif
69 }
70
71 //====================================================================
72 UShort_t
73 AliForwardUtil::ParseCollisionSystem(const char* sys)
74 {
75   // 
76   // Parse a collision system spec given in a string.   Known values are 
77   // 
78   //  - "ppb", "p-pb", "pa", "p-a"  which returns kPPb
79   //  - "pp", "p-p"                 which returns kPP 
80   //  - "PbPb", "Pb-Pb", "A-A",     which returns kPbPb 
81   //  - Everything else gives kUnknown 
82   // 
83   // Parameters:
84   //    sys Collision system spec 
85   // 
86   // Return:
87   //    Collision system id 
88   //
89   TString s(sys);
90   s.ToLower();
91   // we do pA first to avoid pp catch on ppb string (AH)
92   if (s.Contains("p-pb")  || s.Contains("ppb"))   return AliForwardUtil::kPPb;
93   if (s.Contains("p-a")   || s.Contains("pa"))    return AliForwardUtil::kPPb;
94   if (s.Contains("a-p")   || s.Contains("ap"))    return AliForwardUtil::kPPb;
95   if (s.Contains("p-p")   || s.Contains("pp"))    return AliForwardUtil::kPP; 
96   if (s.Contains("pb-pb") || s.Contains("pbpb"))  return AliForwardUtil::kPbPb;
97   if (s.Contains("a-a")   || s.Contains("aa"))    return AliForwardUtil::kPbPb;
98   return AliForwardUtil::kUnknown;
99 }
100 //____________________________________________________________________
101 const char*
102 AliForwardUtil::CollisionSystemString(UShort_t sys)
103 {
104   // 
105   // Get a string representation of the collision system 
106   // 
107   // Parameters:
108   //    sys  Collision system 
109   // - kPP -> "pp"
110   // - kPbPb -> "PbPb" 
111   // - anything else gives "unknown"
112   // 
113   // Return:
114   //    String representation of the collision system 
115   //
116   switch (sys) { 
117   case AliForwardUtil::kPP:   return "pp";
118   case AliForwardUtil::kPbPb: return "PbPb";
119   case AliForwardUtil::kPPb:  return "pPb";
120   }
121   return "unknown";
122 }
123 //____________________________________________________________________
124 Float_t
125 AliForwardUtil::BeamRapidity(Float_t beam, UShort_t z, UShort_t a)
126 {
127   const Double_t pMass = 9.38271999999999995e-01;
128   const Double_t nMass = 9.39564999999999984e-01;
129   Double_t       beamE = z * beam / 2;
130   Double_t       beamM = z * pMass + (a - z) * nMass;
131   Double_t       beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
132   Double_t       beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
133   return beamY;
134 }
135 //____________________________________________________________________
136 Float_t
137 AliForwardUtil::CenterOfMassEnergy(Float_t beam, 
138                                    UShort_t z1, 
139                                    UShort_t a1, 
140                                    Short_t z2, 
141                                    Short_t a2) 
142 {
143   // Calculate the center of mass energy given target/projectile 
144   // mass and charge numbers
145   if (z2 < 0) z2 = z1;
146   if (a2 < 0) a2 = a1;
147   return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
148 }
149 //____________________________________________________________________
150 Float_t
151 AliForwardUtil::CenterOfMassRapidity(UShort_t z1, 
152                                      UShort_t a1, 
153                                      Short_t z2, 
154                                      Short_t a2) 
155 {
156   // Calculate the center of mass rapidity (shift) given target/projectile 
157   // mass and charge numbers
158   if (z2 < 0) z2 = z1;
159   if (a2 < 0) a2 = a1;
160   if (z2 == z1 && a2 == a1) return 0;
161   return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
162 }
163
164 namespace {
165   UShort_t CheckSNN(Float_t energy)
166   {
167     if (TMath::Abs(energy - 900.)   < 10)  return 900;
168     if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
169     if (TMath::Abs(energy - 2760.)  < 20)  return 2760;
170     if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
171     if (TMath::Abs(energy - 5022.)  < 10)  return 5023;
172     if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
173     if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
174     if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
175     if (TMath::Abs(energy - 10000.) < 10)  return 10000;
176     if (TMath::Abs(energy - 14000.) < 10)  return 14000;
177     return 0;
178   }
179 }
180 //____________________________________________________________________
181 UShort_t
182 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t beam)
183 {
184   // 
185   // Parse the center of mass energy given as a float and return known 
186   // values as a unsigned integer
187   // 
188   // Parameters:
189   //    sys   Collision system (needed for AA)
190   //    beam  Center of mass energy * total charge 
191   // 
192   // Return:
193   //    Center of mass energy per nucleon
194   //
195   Float_t energy = beam; 
196   // Below no longer needed apparently
197   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
198   if (sys == AliForwardUtil::kPPb) 
199     energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
200   else if (sys == AliForwardUtil::kPbPb) 
201     energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
202   UShort_t ret = CheckSNN(energy);
203   if (ret > 1) return ret;
204   if (sys == AliForwardUtil::kPbPb || sys == AliForwardUtil::kPPb) {
205     ret = CheckSNN(beam);
206   }
207   return ret;
208 }
209 //____________________________________________________________________
210 const char* 
211 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
212 {
213   // 
214   // Get a string representation of the center of mass energy per nuclean
215   // 
216   // Parameters:
217   //    cms  Center of mass energy per nucleon
218   // 
219   // Return:
220   //    String representation of the center of mass energy per nuclean
221   //
222   return Form("%04dGeV", cms);
223 }
224 //____________________________________________________________________
225 Short_t
226 AliForwardUtil::ParseMagneticField(Float_t v)
227 {
228   // 
229   // Parse the magnetic field (in kG) as given by a floating point number
230   // 
231   // Parameters:
232   //    field  Magnetic field in kG 
233   // 
234   // Return:
235   //    Short integer value of magnetic field in kG 
236   //
237   if (TMath::Abs(v - 5.) < 1 ) return +5;
238   if (TMath::Abs(v + 5.) < 1 ) return -5;
239   if (TMath::Abs(v) < 1)       return 0;
240   return 999;
241 }
242 //____________________________________________________________________
243 const char* 
244 AliForwardUtil::MagneticFieldString(Short_t f)
245 {
246   // 
247   // Get a string representation of the magnetic field
248   // 
249   // Parameters:
250   //    field Magnetic field in kG
251   // 
252   // Return:
253   //    String representation of the magnetic field
254   //
255   return Form("%01dkG", f);
256 }
257 //_____________________________________________________________________
258 AliAODEvent* AliForwardUtil::GetAODEvent(AliAnalysisTaskSE* task)
259 {
260   // Check if AOD is the output event
261   if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
262
263   AliAODEvent* ret = task->AODEvent();
264   if (ret) return ret; 
265   
266   // Check if AOD is the input event 
267   ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
268   if (!ret) ::Warning("GetAODEvent", "No AOD event found");
269   
270   return ret; 
271 }
272 //_____________________________________________________________________
273 UShort_t AliForwardUtil::CheckForAOD()
274 {
275   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
276   if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
277     // ::Info("CheckForAOD", "Found AOD Input handler");
278     return 1;
279   }
280   if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
281     // ::Info("CheckForAOD", "Found AOD Output handler");
282     return 2;
283   }
284
285   ::Warning("CheckForAOD", 
286             "Neither and input nor output AOD handler is specified");
287   return 0;
288 }
289 //_____________________________________________________________________
290 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
291 {
292   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
293   if (!cls) { 
294     AliAnalysisTask* t = am->GetTask(clsOrName);
295     if (!t) { 
296       ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
297       return false;
298     }
299     ::Info("CheckForTask", "Found task %s", clsOrName);
300     return true;
301   }
302   TClass* dep = gROOT->GetClass(clsOrName);
303   if (!dep) { 
304     ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
305     return false;
306   }
307   TIter next(am->GetTasks());
308   TObject* o = 0;
309   while ((o = next())) { 
310     if (o->IsA()->InheritsFrom(dep)) {
311       ::Info("CheckForTask", "Found task of class %s: %s", 
312              clsOrName, o->GetName());
313       return true;
314     }
315   }
316   ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
317   return false;
318 }
319
320 //_____________________________________________________________________
321 TObject* AliForwardUtil::MakeParameter(const Char_t* name, UShort_t value)
322 {
323   TParameter<int>* ret = new TParameter<int>(name, value);
324   ret->SetMergeMode('f');
325   ret->SetUniqueID(value);
326   return ret;
327 }
328 //_____________________________________________________________________
329 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Int_t value)
330 {
331   TParameter<int>* ret = new TParameter<int>(name, value);
332   ret->SetMergeMode('f');
333   ret->SetUniqueID(value);
334   return ret;
335 }
336 //_____________________________________________________________________
337 TObject* AliForwardUtil::MakeParameter(const Char_t* name, ULong_t value)
338 {
339   TParameter<Long_t>* ret = new TParameter<Long_t>(name, value);
340   ret->SetMergeMode('f');
341   ret->SetUniqueID(value);
342   return ret;
343 }
344 //_____________________________________________________________________
345 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Double_t value)
346 {
347   TParameter<double>* ret = new TParameter<double>(name, value);
348   // Float_t v = value;
349   // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
350   ret->SetMergeMode('f');
351   // ret->SetUniqueID(*tmp);
352   return ret;
353 }
354 //_____________________________________________________________________
355 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
356 {
357   TParameter<bool>* ret = new TParameter<bool>(name, value);
358   ret->SetMergeMode('f');
359   ret->SetUniqueID(value);
360   return ret;
361 }
362
363 //_____________________________________________________________________
364 void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
365 {
366   if (!o) return;
367   TParameter<int>* p = static_cast<TParameter<int>*>(o);
368   if (p->TestBit(BIT(19)))
369     value = p->GetVal(); 
370   else
371     value = o->GetUniqueID();
372 }
373 //_____________________________________________________________________
374 void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
375 {
376   if (!o) return;
377   TParameter<int>* p = static_cast<TParameter<int>*>(o);
378   if (p->TestBit(BIT(19)))
379     value = p->GetVal(); 
380   else
381     value = o->GetUniqueID();
382 }
383 //_____________________________________________________________________
384 void AliForwardUtil::GetParameter(TObject* o, ULong_t& value)
385 {
386   if (!o) return;
387   TParameter<Long_t>* p = static_cast<TParameter<Long_t>*>(o);
388   if (p->TestBit(BIT(19)))
389     value = p->GetVal(); 
390   else
391     value = o->GetUniqueID();
392 }
393 //_____________________________________________________________________
394 void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
395 {
396   if (!o) return;
397   TParameter<double>* p = static_cast<TParameter<double>*>(o);
398   if (p->TestBit(BIT(19)))
399     value = p->GetVal(); // o->GetUniqueID();
400   else {
401     UInt_t  i = o->GetUniqueID();
402     Float_t v = *reinterpret_cast<Float_t*>(&i);
403     value = v;
404   }
405 }
406 //_____________________________________________________________________
407 void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
408 {
409   if (!o) return;
410   TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
411   if (p->TestBit(BIT(19)))
412     value = p->GetVal(); // o->GetUniqueID();
413   else
414     value = o->GetUniqueID();
415 }
416   
417 #if 0
418 //_____________________________________________________________________
419 Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
420 {
421   // Get max R of ring
422   // 
423   // Optimized version that has a cache 
424   static TArrayD inner;
425   static TArrayD outer; 
426   if (inner.GetSize() <= 0 || outer.GetSize() <= 0) {
427     const Double_t minR[] = {  4.5213, 15.4 };
428     const Double_t maxR[] = { 17.2,    28.0 };
429     const Int_t    nStr[] = { 512,     256  };
430     for (Int_t q = 0; q < 2; q++) { 
431       TArrayD& a = (q == 0 ? inner : outer);
432       a.Set(nStr[q]);
433
434       for (Int_t it = 0; it < nStr[q]; it++) {
435         Double_t   rad     = maxR[q] - minR[q];
436         Double_t   segment = rad / nStr[q];
437         Double_t   r       = minR[q] + segment*strip;
438         a[it]              = r;
439       }
440     }
441   }
442   if (ring == 'I' || ring == 'i') return inner.At(strip);
443   return outer.At(strip);
444 }
445 #else
446 //_____________________________________________________________________
447 Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
448 {
449   // Get max R of ring
450   // 
451   // New implementation has only one branch
452   const Double_t minR[] = {  4.5213, 15.4 };
453   const Double_t maxR[] = { 17.2,    28.0 };
454   const Int_t    nStr[] = { 512,     256  };
455
456   Int_t      q       = (ring == 'I' || ring == 'i') ? 0 : 1;  
457   Double_t   rad     = maxR[q] - minR[q];
458   Double_t   segment = rad / nStr[q];
459   Double_t   r       = minR[q] + segment*strip;
460
461   return r;
462 }
463 #endif
464
465 #if 1
466 //_____________________________________________________________________
467 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, 
468                                          UShort_t sec, UShort_t strip, 
469                                          Double_t zvtx)
470 {
471   // Calculate eta from strip with vertex (redundant with
472   // AliESDFMD::Eta but support displaced vertices)
473   //
474   // Slightly more optimized version that uses less branching 
475   
476   // Get R of the strip
477   Double_t   r         = GetStripR(ring, strip);
478   Int_t      hybrid    = sec / 2;
479   Int_t      q        = (ring == 'I' || ring == 'i') ? 0 : 1;
480
481   const Double_t zs[][2] = { { 320.266, -999999 }, 
482                             {  83.666,  74.966 },
483                             { -63.066, -74.966 } };
484   if (det > 3 || zs[det-1][q] == -999999) return -999999;
485
486   Double_t z = zs[det-1][q];
487   if ((hybrid % 2) == 0) z -= .5;
488   
489   Double_t   theta = TMath::ATan2(r,z-zvtx);
490   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
491   
492   return eta;
493 }
494 #else
495 //_____________________________________________________________________
496 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, 
497                                          UShort_t sec, UShort_t strip, 
498                                          Double_t zvtx)
499 {
500   // Calculate eta from strip with vertex (redundant with
501   // AliESDFMD::Eta but support displaced vertices)
502   
503   //Get max R of ring
504   Double_t   r         = GetStripR(ring, strip);
505   Int_t      hybrid    = sec / 2;
506   Bool_t     inner     = (ring == 'I' || ring == 'i');
507   Double_t   z         = 0;
508
509
510   switch (det) { 
511   case 1: z = 320.266;                     break;
512   case 2: z = (inner ?  83.666 :  74.966); break;
513   case 3: z = (inner ? -63.066 : -74.966); break; 
514   default: return -999999;
515   }
516   if ((hybrid % 2) == 0) z -= .5;
517   
518   Double_t   theta = TMath::ATan2(r,z-zvtx);
519   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
520   
521   return eta;
522 }
523 #endif
524
525 //_____________________________________________________________________
526 Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip, 
527                                          Double_t phi,
528                                          Double_t xvtx, Double_t yvtx)
529 {
530   // Calculate eta from strip with vertex (redundant with
531   // AliESDFMD::Eta but support displaced vertices)
532
533   // Unknown x,y -> no change
534   if (yvtx > 999 || xvtx > 999) return phi;
535   
536   //Get max R of ring
537   Double_t r   = GetStripR(ring, strip);
538   Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
539   Double_t pha = (TMath::Abs(yvtx) < 1e-12  ? 0 : TMath::ATan2(xvtx, yvtx));
540   Double_t cha = amp * TMath::Cos(phi+pha);
541   phi += cha;
542   if (phi < 0)              phi += TMath::TwoPi();
543   if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
544   return phi;
545 }
546 //====================================================================
547 TAxis*
548 AliForwardUtil::MakeFullIpZAxis(Int_t nCenter)
549 {
550   // Custom vertex axis that will include satellite vertices 
551   // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10 
552   // Nominal vertices are usually in -10 to 10 and we should have 
553   // 10 bins in that range.  That gives us a total of 
554   //
555   //   10+10+10=30 bins 
556   // 
557   // or 31 bin boundaries 
558   if (nCenter % 2 == 1) 
559     // Number of central bins is odd - make it even
560     nCenter--;
561   const Double_t mCenter = 20;
562   const Int_t    nSat    = 10;
563   const Int_t    nBins   = 2*nSat + nCenter;
564   const Int_t    mBin    = nBins / 2;
565   Double_t       dCenter = 2*mCenter / nCenter;
566   TArrayD        bins(nBins+1); 
567   bins[mBin] = 0;
568   for (Int_t i = 1; i <= nCenter/2; i++) { 
569     // Assign from the middle out 
570     Double_t  v  = i * dCenter;
571     // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-i,mBin+i);
572     bins[mBin-i] = -v;
573     bins[mBin+i] = +v;
574   }
575   for (Int_t i = 1; i <= nSat; i++) { 
576     Double_t v = (i+.5) * 37.5;
577     Int_t    o = nCenter/2+i;
578     // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-o,mBin+o);
579     bins[mBin-o] = -v;
580     bins[mBin+o] = +v;
581   }
582   TAxis* a = new TAxis(nBins,bins.GetArray());
583   return a;
584 }
585 void 
586 AliForwardUtil::PrintTask(const TObject& o)
587 {
588   Int_t ind = gROOT->GetDirLevel();
589   if (ind > 0) 
590     // Print indention 
591     std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
592
593   TString t = TString::Format("%s %s", o.GetName(), o.ClassName());
594   const Int_t maxN = 75;
595   std::cout << "--- " << t << " " << std::setfill('-') 
596             << std::setw(maxN-ind-5-t.Length()) << "-" << std::endl;
597 }
598 void
599 AliForwardUtil::PrintName(const char* name)
600 {
601   Int_t ind = gROOT->GetDirLevel();
602   if (ind > 0) 
603     // Print indention 
604     std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
605     
606   // Now print field name 
607   const Int_t maxN  = 29;
608   Int_t       width = maxN - ind;
609   TString     n(name);
610   if (n.Length() > width-1) {
611     // Truncate the string, and put in "..."
612     n.Remove(width-4);
613     n.Append("...");
614   }
615   n.Append(":");
616   std::cout << std::setfill(' ') << std::left << std::setw(width) 
617             << n << std::right << std::flush;
618 }
619 void
620 AliForwardUtil::PrintField(const char* name, const char* value, ...)
621 {
622   PrintName(name);
623
624   // Now format the field value 
625   va_list ap;
626   va_start(ap, value);
627   static char buf[512];
628   vsnprintf(buf, 511, value, ap);
629   buf[511] = '\0';
630   va_end(ap);
631
632   std::cout << buf << std::endl;
633 }
634
635 //====================================================================
636 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
637 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
638 namespace {
639   // 
640   // The shift of the most probable value for the ROOT function TMath::Landau 
641   //
642   const Double_t  mpshift  = -0.22278298;
643   // 
644   // Integration normalisation 
645   //
646   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
647
648   // 
649   // Utility function to use in TF1 defintition 
650   //
651   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
652   {
653     Double_t x        = xp[0];
654     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
655     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
656     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
657     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
658     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
659
660     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
661   }
662
663   Double_t landauGausComposite(Double_t* xp, Double_t* pp)
664   {
665     Double_t x           = xp[0];
666     Double_t cP          = pp[AliForwardUtil::ELossFitter::kC];
667     Double_t deltaP      = pp[AliForwardUtil::ELossFitter::kDelta];
668     Double_t xiP         = pp[AliForwardUtil::ELossFitter::kXi];
669     Double_t sigmaP      = pp[AliForwardUtil::ELossFitter::kSigma];
670     Double_t cS          = pp[AliForwardUtil::ELossFitter::kSigma+1];
671     Double_t deltaS      = deltaP; // pp[AliForwardUtil::ELossFitter::kSigma+2];
672     Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+2/*3*/];
673     Double_t sigmaS      = sigmaP; // pp[AliForwardUtil::ELossFitter::kSigma+4];
674
675     return (cP * AliForwardUtil::LandauGaus(x,deltaP,xiP,sigmaP,0) + 
676             cS * AliForwardUtil::LandauGaus(x,deltaS,xiS,sigmaS,0));
677   }
678     
679   // 
680   // Utility function to use in TF1 defintition 
681   //
682   Double_t landauGausN(Double_t* xp, Double_t* pp) 
683   {
684     Double_t  x        = xp[0];
685     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
686     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
687     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
688     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
689     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
690     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
691     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
692
693     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
694                                                   n, a);
695   }
696   // 
697   // Utility function to use in TF1 defintition 
698   //
699   Double_t landauGausI(Double_t* xp, Double_t* pp) 
700   {
701     Double_t x         = xp[0];
702     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
703     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
704     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
705     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
706     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
707     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
708
709     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
710   }
711
712
713 }
714 //____________________________________________________________________
715 Double_t 
716 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
717 {
718   // 
719   // Calculate the shifted Landau
720   // @f[
721   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
722   // @f]
723   //
724   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
725   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
726   // @f$\Delta=0,\xi=1@f$. 
727   // 
728   // Parameters:
729   //    x      Where to evaluate @f$ f'_{L}@f$ 
730   //    delta  Most probable value 
731   //    xi     The 'width' of the distribution 
732   //
733   // Return:
734   //    @f$ f'_{L}(x;\Delta,\xi) @f$
735   //
736   return TMath::Landau(x, delta - xi * mpshift, xi);
737 }
738 //____________________________________________________________________
739 Double_t 
740 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
741                            Double_t sigma, Double_t sigmaN)
742 {
743   // 
744   // Calculate the value of a Landau convolved with a Gaussian 
745   // 
746   // @f[ 
747   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
748   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
749   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
750   // @f]
751   // 
752   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
753   // energy loss, @f$ \xi@f$ the width of the Landau, and 
754   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
755   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
756   // noise in the detector.  
757   //
758   // Note that this function uses the constants fgConvolutionSteps and
759   // fgConvolutionNSigma
760   // 
761   // References: 
762   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
763   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
764   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
765   // 
766   // Parameters:
767   //    x         where to evaluate @f$ f@f$
768   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
769   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
770   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
771   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
772   // 
773   // Return:
774   //    @f$ f@f$ evaluated at @f$ x@f$.  
775   //
776   Double_t deltap = delta - xi * mpshift;
777   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
778   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
779   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
780   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
781   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
782   Double_t sum    = 0;
783   
784   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
785     Double_t x1 = xlow  + (i - .5) * step;
786     Double_t x2 = xhigh - (i - .5) * step;
787     
788     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
789     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
790   }
791   return step * sum * invSq2pi / sigma1;
792 }
793
794 //____________________________________________________________________
795 Double_t 
796 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
797                             Double_t sigma, Double_t sigmaN, Int_t i)
798 {
799   // 
800   // Evaluate 
801   // @f[ 
802   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
803   // @f] 
804   // corresponding to @f$ i@f$ particles i.e., with the substitutions 
805   // @f{eqnarray*}{ 
806   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
807   //    \xi       \rightarrow \xi_i       &=& i \xi
808   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
809   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
810   // @f} 
811   // 
812   // Parameters:
813   //    x        Where to evaluate 
814   //    delta    @f$ \Delta@f$ 
815   //    xi       @f$ \xi@f$ 
816   //    sigma    @f$ \sigma@f$ 
817   //    sigma_n  @f$ \sigma_n@f$
818   //    i        @f$ i @f$
819   // 
820   // Return:
821   //    @f$ f_i @f$ evaluated
822   //  
823   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
824   Double_t xiI    =  i * xi;
825   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
826   if (sigmaI < 1e-10) { 
827     // Fall back to landau 
828     return AliForwardUtil::Landau(x, deltaI, xiI);
829   }
830   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
831 }
832
833 //____________________________________________________________________
834 Double_t 
835 AliForwardUtil::IdLandauGausdPar(Double_t x, 
836                                  UShort_t par,   Double_t dPar, 
837                                  Double_t delta, Double_t xi, 
838                                  Double_t sigma, Double_t sigmaN, 
839                                  Int_t    i)
840 {
841   // 
842   // Numerically evaluate 
843   // @f[ 
844   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
845   // @f] 
846   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
847   // of the parameters is given by 
848   //
849   // - 0: @f$\Delta@f$ 
850   // - 1: @f$\xi@f$ 
851   // - 2: @f$\sigma@f$ 
852   // - 3: @f$\sigma_n@f$ 
853   //
854   // This is the partial derivative with respect to the parameter of
855   // the response function corresponding to @f$ i@f$ particles i.e.,
856   // with the substitutions
857   // @f[ 
858   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
859   //    \xi       \rightarrow \xi_i       = i \xi
860   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
861   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
862   // @f] 
863   // 
864   // Parameters:
865   //    x        Where to evaluate 
866   //    ipar     Parameter number 
867   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
868   //    delta    @f$ \Delta@f$ 
869   //    xi       @f$ \xi@f$ 
870   //    sigma    @f$ \sigma@f$ 
871   //    sigma_n  @f$ \sigma_n@f$
872   //    i        @f$ i@f$
873   // 
874   // Return:
875   //    @f$ f_i@f$ evaluated
876   //  
877   if (dPar == 0) return 0;
878   Double_t dp      = dPar;
879   Double_t d2      = dPar / 2;
880   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
881   Double_t xiI     =  i * xi;
882   Double_t si      =  TMath::Sqrt(Double_t(i));
883   Double_t sigmaI  =  si*sigma;
884   Double_t y1      = 0;
885   Double_t y2      = 0;
886   Double_t y3      = 0;
887   Double_t y4      = 0;
888   switch (par) {
889   case 0: 
890     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
891     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
892     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
893     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
894     break;
895   case 1: 
896     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
897     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
898     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
899     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
900     break;
901   case 2: 
902     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
903     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
904     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
905     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
906     break;
907   case 3: 
908     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
909     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
910     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
911     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
912     break;
913   default:
914     return 0;
915   } 
916   
917   Double_t d0  = y1 - y4;
918   Double_t d1  = 2 * (y2 - y3);
919   
920   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
921    
922   return g;
923 }
924
925 //____________________________________________________________________
926 Double_t 
927 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
928                             Double_t sigma, Double_t sigmaN, Int_t n, 
929                             const Double_t* a)
930 {
931   // 
932   // Evaluate 
933   // @f[ 
934   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
935   // @f] 
936   // 
937   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
938   // Landau with a Gaussian (see LandauGaus).  Note that 
939   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
940   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
941   //  
942   // References: 
943   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
944   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
945   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
946   // 
947   // Parameters:
948   //    x        Where to evaluate @f$ f_N@f$
949   //    delta    @f$ \Delta_1@f$ 
950   //    xi       @f$ \xi_1@f$
951   //    sigma    @f$ \sigma_1@f$ 
952   //    sigma_n  @f$ \sigma_n@f$ 
953   //    n        @f$ N@f$ in the sum above.
954   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
955   //                 @f$ i > 1@f$ 
956   // 
957   // Return:
958   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
959   //
960   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
961   for (Int_t i = 2; i <= n; i++) 
962     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
963   return result;
964 }
965 namespace { 
966   const Int_t kColors[] = { kRed+1, 
967                             kPink+3, 
968                             kMagenta+2, 
969                             kViolet+2, 
970                             kBlue+1, 
971                             kAzure+3, 
972                             kCyan+1, 
973                             kTeal+2, 
974                             kGreen+2, 
975                             kSpring+3, 
976                             kYellow+2, 
977                             kOrange+2 };
978 }
979
980 //____________________________________________________________________
981 TF1*
982 AliForwardUtil::MakeNLandauGaus(Double_t  c, 
983                                 Double_t  delta, Double_t xi, 
984                                 Double_t  sigma, Double_t sigmaN, Int_t n, 
985                                 const Double_t* a, 
986                                 Double_t  xmin, Double_t xmax)
987 {
988   // 
989   // Generate a TF1 object of @f$ f_N@f$ 
990   // 
991   // Parameters:
992   //    c         Constant                             
993   //    delta     @f$ \Delta@f$                        
994   //    xi            @f$ \xi_1@f$                     
995   //    sigma     @f$ \sigma_1@f$                      
996   //    sigma_n   @f$ \sigma_n@f$                      
997   //    n             @f$ N@f$ - how many particles to sum to
998   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
999   //                  @f$ i > 1@f$ 
1000   //    xmin      Least value of range  
1001   //    xmax      Largest value of range
1002   // 
1003   // Return:
1004   //    Newly allocated TF1 object
1005   //
1006   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
1007   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
1008   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
1009   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
1010   landaun->SetLineWidth(2);
1011   landaun->SetNpx(500);
1012   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
1013
1014   // Set the initial parameters from the seed fit 
1015   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
1016   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
1017   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
1018   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
1019   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
1020   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);       
1021
1022   // Set the range and name of the scale parameters 
1023   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
1024     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
1025     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
1026   }
1027   return landaun;
1028 }
1029 //____________________________________________________________________
1030 TF1*
1031 AliForwardUtil::MakeILandauGaus(Double_t  c, 
1032                                 Double_t  delta, Double_t xi, 
1033                                 Double_t  sigma, Double_t sigmaN, Int_t i, 
1034                                 Double_t  xmin, Double_t xmax)
1035 {
1036   // 
1037   // Generate a TF1 object of @f$ f_I@f$ 
1038   // 
1039   // Parameters:
1040   //    c        Constant
1041   //    delta    @f$ \Delta@f$ 
1042   //    xi       @f$ \xi_1@f$          
1043   //    sigma    @f$ \sigma_1@f$               
1044   //    sigma_n  @f$ \sigma_n@f$               
1045   //    i            @f$ i@f$ - the number of particles
1046   //    xmin     Least value of range
1047   //    xmax     Largest value of range
1048   // 
1049   // Return:
1050   //    Newly allocated TF1 object
1051   //
1052   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
1053   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
1054   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
1055   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
1056   landaui->SetLineWidth(1);
1057   landaui->SetNpx(500);
1058   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
1059
1060   // Set the initial parameters from the seed fit 
1061   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
1062   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
1063   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
1064   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
1065   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
1066   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);       
1067
1068   return landaui;
1069 }
1070
1071 //====================================================================
1072 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
1073                                          Double_t maxRange, 
1074                                          UShort_t minusBins) 
1075   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
1076     fFitResults(0), fFunctions(0), fDebug(false)
1077 {
1078   // 
1079   // Constructor 
1080   // 
1081   // Parameters:
1082   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
1083   //    maxRange   Maximum range to fit to 
1084   //    minusBins  The number of bins below maximum to use 
1085   //
1086   fFitResults.SetOwner();
1087   fFunctions.SetOwner();
1088 }
1089 //____________________________________________________________________
1090 AliForwardUtil::ELossFitter::~ELossFitter()
1091 {
1092   // 
1093   // Destructor
1094   // 
1095   //
1096   fFitResults.Delete();
1097   fFunctions.Delete();
1098 }
1099 //____________________________________________________________________
1100 void
1101 AliForwardUtil::ELossFitter::Clear()
1102 {
1103   // 
1104   // Clear internal arrays 
1105   // 
1106   //
1107   fFitResults.Clear();
1108   fFunctions.Clear();
1109 }
1110 //____________________________________________________________________
1111 TF1*
1112 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
1113 {
1114   // 
1115   // Fit a 1-particle signal to the passed energy loss distribution 
1116   // 
1117   // Note that this function clears the internal arrays first 
1118   // 
1119   // Parameters:
1120   //    dist    Data to fit the function to 
1121   //    sigman If larger than zero, the initial guess of the
1122   //               detector induced noise. If zero or less, then this 
1123   //               parameter is ignored in the fit (fixed at 0)
1124   // 
1125   // Return:
1126   //    The function fitted to the data 
1127   //
1128
1129   // Clear the cache 
1130   Clear();
1131   
1132   // Find the fit range 
1133   // Find the fit range 
1134   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1135   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1136                                 dist->GetNbinsX());
1137   dist->GetXaxis()->SetRange(cutBin, maxBin);
1138   // dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
1139   
1140   // Get the bin with maximum 
1141   Int_t    peakBin = dist->GetMaximumBin();
1142   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1143   
1144   // Get the low edge 
1145   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1146   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1147   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1148   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1149
1150   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1151   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1152   Double_t intg  = dist->Integral(minEb, maxEb);
1153   if (intg <= 0) {
1154     ::Warning("Fit1Particle", 
1155               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
1156               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1157     return 0;
1158   }
1159     
1160   // Restore the range 
1161   dist->GetXaxis()->SetRange(1, maxBin);
1162   
1163   // Define the function to fit 
1164   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
1165
1166   // Set initial guesses, parameter names, and limits  
1167   landau1->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1168   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1169   landau1->SetNpx(500);
1170   if (peakE >= minE && peakE <= fMaxRange) {
1171     // printf("Fit1: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1172     landau1->SetParLimits(kDelta, minE, fMaxRange);
1173   }
1174   if (peakE/10 >= 0 && peakE <= 0.1) {
1175     // printf("Fit1: Set par limits on xi: %f, %f\n", 0., 0.1);
1176     landau1->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1177   }
1178   if (peakE/5 >= 0 && peakE/5 <= 0.1) {
1179     // printf("Fit1: Set par limits on sigma: %f, %f\n", 0., 0.1);
1180     landau1->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1181   }
1182   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
1183   else {
1184     // printf("Fit1: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1185     landau1->SetParLimits(kSigmaN, 0, fMaxRange);
1186   }
1187
1188   // Do the fit, getting the result object 
1189   if (fDebug) 
1190     ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
1191   TFitResultPtr r = dist->Fit(landau1, FIT_OPTIONS, "", minE, maxE);
1192   if (!r.Get()) { 
1193     ::Warning("Fit1Particle", 
1194               "No fit returned when processing %s in the range [%f,%f] "
1195               "options %s", dist->GetName(), minE, maxE, FIT_OPTIONS);
1196     return 0;
1197   }
1198   // landau1->SetRange(minE, fMaxRange);
1199   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
1200   fFunctions.AddAtAndExpand(landau1, 0);
1201
1202   return landau1;
1203 }
1204 //____________________________________________________________________
1205 TF1*
1206 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
1207                                           Double_t sigman)
1208 {
1209   // 
1210   // Fit a N-particle signal to the passed energy loss distribution 
1211   //
1212   // If there's no 1-particle fit present, it does that first 
1213   // 
1214   // Parameters:
1215   //    dist   Data to fit the function to 
1216   //    n      Number of particle signals to fit 
1217   //    sigman If larger than zero, the initial guess of the
1218   //               detector induced noise. If zero or less, then this 
1219   //               parameter is ignored in the fit (fixed at 0)
1220   // 
1221   // Return:
1222   //    The function fitted to the data 
1223   //
1224
1225   // Get the seed fit result 
1226   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
1227   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
1228   if (!r || !f) { 
1229     f = Fit1Particle(dist, sigman);
1230     r = static_cast<TFitResult*>(fFitResults.At(0));
1231     if (!r || !f) { 
1232       ::Warning("FitNLandau", "No first shot at landau fit");
1233       return 0;
1234     }
1235   }
1236
1237   // Get some parameters from seed fit 
1238   Double_t delta1  = r->Parameter(kDelta);
1239   Double_t xi1     = r->Parameter(kXi);
1240   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
1241   Double_t minE    = f->GetXmin();
1242
1243   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1244   Int_t    maxEb = dist->GetXaxis()->FindBin(maxEi);
1245   Double_t intg  = dist->Integral(minEb, maxEb);
1246   if (intg <= 0) {
1247     ::Warning("FitNParticle",
1248               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
1249               dist->GetName(), minE, maxEi, minEb, maxEb, intg);
1250     return 0;
1251   }
1252
1253   // Array of weights 
1254   TArrayD a(n-1);
1255   for (UShort_t i = 2; i <= n; i++) 
1256     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
1257   // Make the fit function 
1258   TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
1259                                  r->Parameter(kDelta),
1260                                  r->Parameter(kXi),
1261                                  r->Parameter(kSigma),
1262                                  r->Parameter(kSigmaN),
1263                                  n, a.fArray, minE, maxEi);
1264   if (minE      <= r->Parameter(kDelta) &&
1265       fMaxRange >= r->Parameter(kDelta)) {
1266     // Protect against warning from ParameterSettings
1267     // printf("FitN: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1268     landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
1269   }
1270   if (r->Parameter(kXi) >= 0 && r->Parameter(kXi) <= 0.1) {
1271     // printf("FitN: Set par limits on xi: %f, %f\n", 0., 0.1);
1272     landaun->SetParLimits(kXi,     0.00, 0.1);  // was fMaxRange - too wide
1273   }
1274   if (r->Parameter(kSigma) >= 1e-5 && r->Parameter(kSigma) <= 0.1) {
1275     // printf("FitN: Set par limits on sigma: %f, %f\n", 1e-5, 0.1);
1276     landaun->SetParLimits(kSigma,  1e-5, 0.1);  // was fMaxRange - too wide
1277   }
1278   // Check if we're using the noise sigma 
1279   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
1280   else {
1281     // printf("FitN: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1282     landaun->SetParLimits(kSigmaN, 0, fMaxRange);
1283   }
1284
1285   // Set the range and name of the scale parameters 
1286   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
1287     if (a[i-2] >= 0 && a[i-2] <= 1) {
1288       // printf("FitN: Set par limits on a_%d: %f, %f\n", i, 0., 1.);
1289       landaun->SetParLimits(kA+i-2, 0,1);
1290     }
1291   }
1292
1293   // Do the fit 
1294   if (fDebug) 
1295     ::Info("FitNParticle", "Fitting in the range %f,%f (%d)", minE, maxEi, n);
1296   TFitResultPtr tr = dist->Fit(landaun, FIT_OPTIONS, "", minE, maxEi);
1297   
1298   // landaun->SetRange(minE, fMaxRange);
1299   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
1300   fFunctions.AddAtAndExpand(landaun, n-1);
1301   
1302   return landaun;
1303 }  
1304 //____________________________________________________________________
1305 TF1*
1306 AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
1307 {
1308   // 
1309   // Fit a composite particle signal to the passed energy loss
1310   // distribution
1311   // 
1312   // Parameters:
1313   //    dist    Data to fit the function to 
1314   //    sigman If larger than zero, the initial guess of the
1315   //               detector induced noise. If zero or less, then this 
1316   //               parameter is ignored in the fit (fixed at 0)
1317   // 
1318   // Return:
1319   //    The function fitted to the data 
1320   //
1321
1322   // Find the fit range 
1323   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1324   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1325                                 dist->GetNbinsX());
1326   dist->GetXaxis()->SetRange(cutBin, maxBin);
1327   
1328   // Get the bin with maximum 
1329   Int_t    peakBin = dist->GetMaximumBin();
1330   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1331   
1332   // Get the low edge 
1333   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1334   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1335   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1336   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1337
1338   // Get the range in bins and the integral of that range 
1339   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1340   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1341   Double_t intg  = dist->Integral(minEb, maxEb);
1342   if (intg <= 0) {
1343     ::Warning("Fit1Particle", 
1344               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
1345               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1346     return 0;
1347   }
1348     
1349   // Restore the range 
1350   dist->GetXaxis()->SetRange(1, maxBin);
1351   
1352   // Define the function to fit 
1353   TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
1354
1355   // Set initial guesses, parameter names, and limits  
1356   seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1357   seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1358   seed->SetNpx(500);
1359   seed->SetParLimits(kDelta, minE, fMaxRange);
1360   seed->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1361   seed->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1362   if (sigman <= 0)  seed->FixParameter(kSigmaN, 0);
1363   else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
1364
1365   // Do the fit, getting the result object 
1366   if (fDebug) 
1367     ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
1368   /* TFitResultPtr r = */ dist->Fit(seed, FIT_OPTIONS, "", minE, maxE);
1369
1370   maxE = dist->GetXaxis()->GetXmax();
1371 #if 1
1372   TF1* comp = new TF1("composite", landauGausComposite, 
1373                       minE, maxE, kSigma+1+2);
1374   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
1375                     "C#prime", "#xi#prime");
1376   comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight 
1377                       seed->GetParameter(kDelta),    // 1 Primary Delta
1378                       seed->GetParameter(kDelta)/10, // 2 primary Xi
1379                       seed->GetParameter(kDelta)/5,  // 3 primary sigma
1380                       1.20 * seed->GetParameter(kC), // 5 Secondary weight
1381                       seed->GetParameter(kXi));      // 7 secondary Xi
1382   // comp->SetParLimits(kC,       minE, fMaxRange); // C
1383   comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
1384   comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi 
1385   comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
1386   // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
1387   comp->SetParLimits(kSigma+2,    0.00, fMaxRange); // Xi'
1388 #else
1389   TF1* comp = new TF1("composite", landauGausComposite, 
1390                       minE, maxE, kSigma+1+4);
1391   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
1392                     "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
1393   comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight 
1394                       seed->GetParameter(kDelta),    // 1 Primary Delta
1395                       seed->GetParameter(kDelta)/10, // 2 primary Xi
1396                       seed->GetParameter(kDelta)/5,  // 3 primary sigma
1397                       1.20 * seed->GetParameter(kC), // 5 Secondary weight
1398                       seed->GetParameter(kDelta),    // 6 secondary Delta
1399                       seed->GetParameter(kXi),       // 7 secondary Xi
1400                       seed->GetParameter(kSigma));   // 8 secondary sigma
1401   // comp->SetParLimits(kC,       minE, fMaxRange); // C
1402   comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
1403   comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi 
1404   comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
1405   // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
1406   comp->SetParLimits(kSigma+2,    minE/10, fMaxRange); // Delta
1407   comp->SetParLimits(kSigma+3,    0.00,    fMaxRange); // Xi 
1408   comp->SetParLimits(kSigma+4,    1e-6,    fMaxRange); // Sigma
1409 #endif                
1410   comp->SetLineColor(kRed+1);
1411   comp->SetLineWidth(3);
1412   
1413   // Do the fit, getting the result object 
1414   if (fDebug) 
1415     ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
1416   /* TFitResultPtr r = */ dist->Fit(comp, FIT_OPTIONS, "", minE, maxE);
1417
1418 #if 0
1419   TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
1420   part1->SetLineColor(kGreen+1);
1421   part1->SetLineWidth(4);
1422   part1->SetRange(minE, maxE);
1423   part1->SetParameters(comp->GetParameter(0), // C 
1424                        comp->GetParameter(1), // Delta
1425                        comp->GetParameter(2), // Xi
1426                        comp->GetParameter(3), // sigma
1427                        0);
1428   part1->Save(minE,maxE,0,0,0,0);
1429   dist->GetListOfFunctions()->Add(part1);
1430
1431   TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
1432   part2->SetLineColor(kBlue+1);
1433   part2->SetLineWidth(4);
1434   part2->SetRange(minE, maxE);
1435   part2->SetParameters(comp->GetParameter(4), // C 
1436                        comp->GetParameter(5), // Delta
1437                        comp->GetParameter(6), // Xi
1438                        comp->GetParameter(7), // sigma
1439                        0);
1440   part2->Save(minE,maxE,0,0,0,0);
1441   dist->GetListOfFunctions()->Add(part2);
1442 #endif
1443   return comp;
1444 }
1445
1446 //====================================================================
1447 AliForwardUtil::Histos::~Histos()
1448 {
1449   // 
1450   // Destructor
1451   //
1452 }
1453
1454 //____________________________________________________________________
1455 void
1456 AliForwardUtil::Histos::Delete(Option_t* opt)
1457 {
1458   if (fFMD1i) delete fFMD1i;
1459   if (fFMD2i) delete fFMD2i;
1460   if (fFMD2o) delete fFMD2o;
1461   if (fFMD3i) delete fFMD3i;
1462   if (fFMD3o) delete fFMD3o;
1463   fFMD1i = 0;
1464   fFMD2i = 0;
1465   fFMD2o = 0;
1466   fFMD3i = 0;
1467   fFMD3o = 0;
1468   TObject::Delete(opt);
1469 }
1470
1471 //____________________________________________________________________
1472 TH2D*
1473 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, const TAxis& etaAxis)
1474 {
1475   // 
1476   // Make a histogram 
1477   // 
1478   // Parameters:
1479   //    d        Detector
1480   //    r        Ring 
1481   //    etaAxis  Eta axis to use
1482   // 
1483   // Return:
1484   //    Newly allocated histogram 
1485   //
1486   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
1487   TH2D* hist = 0;
1488   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1489     hist = new TH2D(Form("FMD%d%c_cache", d, r), 
1490                     Form("FMD%d%c cache", d, r),
1491                     etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(), 
1492                     ns, 0, TMath::TwoPi());
1493   else
1494     hist = new TH2D(Form("FMD%d%c_cache", d, r), 
1495                     Form("FMD%d%c cache", d, r),
1496                     etaAxis.GetNbins(), etaAxis.GetXmin(), 
1497                     etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
1498   hist->SetXTitle("#eta");
1499   hist->SetYTitle("#phi [radians]");
1500   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
1501   hist->Sumw2();
1502   hist->SetDirectory(0);
1503
1504   return hist;
1505 }
1506 //____________________________________________________________________
1507 void
1508 AliForwardUtil::Histos::RebinEta(TH2D* hist, const TAxis& etaAxis)
1509 {
1510   TAxis* xAxis = hist->GetXaxis();
1511   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1512     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
1513   else
1514     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
1515   hist->Rebuild();
1516 }
1517   
1518
1519 //____________________________________________________________________
1520 void
1521 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
1522 {
1523   // 
1524   // Initialize the object 
1525   // 
1526   // Parameters:
1527   //    etaAxis Eta axis to use 
1528   //
1529   fFMD1i = Make(1, 'I', etaAxis);
1530   fFMD2i = Make(2, 'I', etaAxis);
1531   fFMD2o = Make(2, 'O', etaAxis);
1532   fFMD3i = Make(3, 'I', etaAxis);
1533   fFMD3o = Make(3, 'O', etaAxis);
1534 }
1535 //____________________________________________________________________
1536 void
1537 AliForwardUtil::Histos::ReInit(const TAxis& etaAxis)
1538 {
1539   // 
1540   // Initialize the object 
1541   // 
1542   // Parameters:
1543   //    etaAxis Eta axis to use 
1544   //
1545   if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
1546   if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
1547   if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
1548   if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
1549   if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
1550 }
1551
1552 //____________________________________________________________________
1553 void
1554 AliForwardUtil::Histos::Clear(Option_t* option)
1555 {
1556   // 
1557   // Clear data 
1558   // 
1559   // Parameters:
1560   //    option Not used 
1561   //
1562   if (fFMD1i) fFMD1i->Reset(option);
1563   if (fFMD2i) fFMD2i->Reset(option);
1564   if (fFMD2o) fFMD2o->Reset(option);
1565   if (fFMD3i) fFMD3i->Reset(option);
1566   if (fFMD3o) fFMD3o->Reset(option);
1567 }
1568
1569 //____________________________________________________________________
1570 TH2D*
1571 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
1572 {
1573   // 
1574   // Get the histogram for a particular detector,ring
1575   // 
1576   // Parameters:
1577   //    d Detector 
1578   //    r Ring 
1579   // 
1580   // Return:
1581   //    Histogram for detector,ring or nul 
1582   //
1583   switch (d) { 
1584   case 1: return fFMD1i;
1585   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1586   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1587   }
1588   return 0;
1589 }
1590 //====================================================================
1591 TList*
1592 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1593 {
1594   // 
1595   // Define the outout list in @a d
1596   // 
1597   // Parameters:
1598   //    d Where to put the output list
1599   // 
1600   // Return:
1601   //    Newly allocated TList object or null
1602   //
1603   if (!d) return 0;
1604   TList* list = new TList;
1605   list->SetOwner();
1606   list->SetName(fName.Data());
1607   d->Add(list);
1608   return list;
1609 }
1610 //____________________________________________________________________
1611 TList*
1612 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1613 {
1614   // 
1615   // Get our output list from the container @a d
1616   // 
1617   // Parameters:
1618   //    d where to get the output list from 
1619   // 
1620   // Return:
1621   //    The found TList or null
1622   //
1623   if (!d) return 0;
1624   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1625   return list;
1626 }
1627
1628 //____________________________________________________________________
1629 TH1*
1630 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1631 {
1632   // 
1633   // Find a specific histogram in the source list @a d
1634   // 
1635   // Parameters:
1636   //    d     (top)-container 
1637   //    name  Name of histogram
1638   // 
1639   // Return:
1640   //    Found histogram or null
1641   //
1642   return static_cast<TH1*>(d->FindObject(name));
1643 }
1644
1645 //====================================================================
1646 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl, 
1647                                        const char* format, ...)
1648   : fMsg("")
1649 {
1650   if (lvl < msgLvl) return; 
1651   va_list ap;
1652   va_start(ap, format);
1653   Format(fMsg, format, ap);
1654   va_end(ap);
1655   Output(+1, fMsg);
1656 }
1657 //____________________________________________________________________
1658 AliForwardUtil::DebugGuard::~DebugGuard()
1659 {
1660   if (fMsg.IsNull()) return;
1661   Output(-1, fMsg);
1662 }
1663 //____________________________________________________________________
1664 void
1665 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl, 
1666                                     const char* format, ...)
1667 {
1668   if (lvl < msgLvl) return; 
1669   TString msg;
1670   va_list ap;
1671   va_start(ap, format);
1672   Format(msg, format, ap);
1673   va_end(ap);
1674   Output(0, msg);
1675 }
1676
1677 //____________________________________________________________________
1678 void
1679 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1680 {
1681   static char buf[512];
1682   Int_t n = gROOT->GetDirLevel() + 2;
1683   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1684   vsnprintf(&(buf[n]), 511-n, format, ap);
1685   buf[511] = '\0';
1686   out = buf;  
1687 }
1688 //____________________________________________________________________
1689 void
1690 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1691 {
1692   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1693   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1694   if      (in > 0) gROOT->IncreaseDirLevel();
1695   else if (in < 0) gROOT->DecreaseDirLevel();
1696 }
1697
1698
1699
1700 //
1701 // EOF
1702 //