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