PWG2/FORWARD/analysis2/AliForwardUtil.cxx

   1 #include "AliForwardUtil.h"
   2 #include <AliAnalysisManager.h>
   3 #include "AliAODForwardMult.h"
   4 #include <AliLog.h>
   5 #include <AliInputEventHandler.h>
   6 #include <AliESDEvent.h>
   7 #include <AliPhysicsSelection.h>
   8 #include <AliTriggerAnalysis.h>
   9 #include <AliMultiplicity.h>
  10 #include <TH2D.h>
  11 #include <TH1I.h>
  12 #include <TF1.h>
  13 #include <TFitResult.h>
  14 #include <TMath.h>
  15 #include <TError.h>
  16
  17 //====================================================================
  18 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
  19 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
  20 namespace {
  21   /**
  22    * The shift of the most probable value for the ROOT function TMath::Landau
  23    */
  24   const Double_t  mpshift  = -0.22278298;
  25   /**
  26    * Integration normalisation
  27    */
  28   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
  29
  30   /**
  31    * Utility function to use in TF1 defintition
  32    */
  33   Double_t landauGaus1(Double_t* xp, Double_t* pp)
  34   {
  35     Double_t x        = xp[0];
  36     Double_t constant = pp[0];
  37     Double_t delta    = pp[1];
  38     Double_t xi       = pp[2];
  39     Double_t sigma    = pp[3];
  40     Double_t sigma_n  = pp[4];
  41
  42     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
  43   }
  44
  45   /**
  46    * Utility function to use in TF1 defintition
  47    */
  48   Double_t landauGausN(Double_t* xp, Double_t* pp)
  49   {
  50     Double_t  x        = xp[0];
  51     Double_t  constant = pp[0];
  52     Double_t  delta    = pp[1];
  53     Double_t  xi       = pp[2];
  54     Double_t  sigma    = pp[3];
  55     Double_t  sigma_n  = pp[4];
  56     Int_t     n        = Int_t(pp[5]);
  57     Double_t* a        = &(pp[6]);
  58
  59     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
  60                                                   n, a);
  61   }
  62
  63
  64 }
  65 //____________________________________________________________________
  66 Double_t
  67 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
  68 {
  69   return TMath::Landau(x, delta - xi * mpshift, xi);
  70 }
  71 //____________________________________________________________________
  72 Double_t
  73 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
  74                            Double_t sigma, Double_t sigma_n)
  75 {
  76   Double_t deltap = delta - xi * mpshift;
  77   Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
  78   Double_t sigma1 = TMath::Sqrt(sigma2);
  79   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
  80   Double_t xhigh  = x - fgConvolutionNSigma * sigma1;
  81   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
  82   Double_t sum    = 0;
  83
  84   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
  85     Double_t x1 = xlow + (i - .5) * step;
  86     Double_t x2 = xlow - (i - .5) * step;
  87
  88     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
  89     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
  90   }
  91   return step * sum * invSq2pi / sigma1;
  92 }
  93
  94 //____________________________________________________________________
  95 Double_t
  96 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
  97                             Double_t sigma, Double_t sigma_n, Int_t n,
  98                             Double_t* a)
  99 {
 100   Double_t result = LandauGaus(x, delta, xi, sigma, sigma_n);
 101   for (Int_t i = 2; i <= n; i++) {
 102     Double_t delta_i =  i * (delta + xi * TMath::Log(i));
 103     Double_t xi_i    =  i * xi;
 104     Double_t sigma_i =  TMath::Sqrt(Double_t(n)*sigma);
 105     Double_t a_i     =  a[i];
 106     result           += a_i * AliForwardUtil::LandauGaus(x, delta_i, xi_i,
 107                                                          sigma_i, sigma_n);
 108   }
 109   return result;
 110 }
 111
 112 //====================================================================
 113 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
 114                                          Double_t maxRange,
 115                                          UShort_t minusBins)
 116   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
 117     fFitResults(0), fFunctions(0)
 118 {
 119   fFitResults.SetOwner();
 120   fFunctions.SetOwner();
 121 }
 122 //____________________________________________________________________
 123 AliForwardUtil::ELossFitter::~ELossFitter()
 124 {
 125   fFitResults.Delete();
 126   fFunctions.Delete();
 127 }
 128 //____________________________________________________________________
 129 void
 130 AliForwardUtil::ELossFitter::Clear()
 131 {
 132   fFitResults.Clear();
 133   fFunctions.Clear();
 134 }
 135 //____________________________________________________________________
 136 TF1*
 137 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
 138 {
 139   // Clear the cache
 140   Clear();
 141
 142   // Find the fit range
 143   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
 144
 145   // Normalize peak to 1
 146   Double_t max = dist->GetMaximum();
 147   dist->Scale(1/max);
 148
 149   // Get the bin with maximum
 150   Int_t    maxBin = dist->GetMaximumBin();
 151   Double_t maxE   = dist->GetBinLowEdge(maxBin);
 152
 153   // Get the low edge
 154   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
 155   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
 156   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
 157   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
 158
 159   // Restore the range
 160   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
 161
 162   // Define the function to fit
 163   TF1*          landau1 = new TF1("landau1", landauGaus1, minE, maxEE, 5);
 164
 165   // Set initial guesses, parameter names, and limits
 166   landau1->SetParameters(5,.5,.07,1,sigman);
 167   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
 168   landau1->SetParLimits(1, minE, fMaxRange);
 169   landau1->SetParLimits(2, 0.00, fMaxRange);
 170   landau1->SetParLimits(3, 0.01, fMaxRange);
 171   if (sigman <= 0)  landau1->FixParameter(4, 0);
 172   else              landau1->SetParLimits(4, 0, fMaxRange);
 173
 174   // Do the fit, getting the result object
 175   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
 176
 177   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
 178   fFunctions.AddAtAndExpand(landau1, 0);
 179
 180   return landau1;
 181 }
 182 //____________________________________________________________________
 183 TF1*
 184 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
 185                                           Double_t sigman)
 186 {
 187   // Get the seed fit result
 188   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
 189   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
 190   if (!r || !f) {
 191     f = Fit1Particle(dist, sigman);
 192     r = static_cast<TFitResult*>(fFitResults.At(0));
 193     if (!r || !f) {
 194       ::Warning("FitNLandau", "No first shot at landau fit");
 195       return 0;
 196     }
 197   }
 198
 199   // Get some parameters from seed fit
 200   Double_t delta1  = r->Parameter(1);
 201   Double_t xi1     = r->Parameter(2);
 202   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
 203   Double_t minE    = f->GetXmin();
 204
 205   // Make the fit function
 206   TF1* landaun     = new TF1(Form("landau%d", n), &landauGausN,minE,maxEi,5+n);
 207   landaun->SetLineStyle((n % 10)+1);
 208   landaun->SetLineWidth(2);
 209   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 210
 211   // Set the initial parameters from the seed fit
 212   landaun->SetParameter(0, r->Parameter(0)); // Constant
 213   landaun->SetParameter(1, r->Parameter(1)); // Delta
 214   landaun->SetParameter(2, r->Parameter(2)); // xi
 215   landaun->SetParameter(3, r->Parameter(3)); // sigma
 216   landaun->SetParameter(4, r->Parameter(4)); // sigma_n
 217   landaun->SetParLimits(1, minE, fMaxRange); // Delta
 218   landaun->SetParLimits(2, 0.00, fMaxRange); // xi
 219   landaun->SetParLimits(3, 0.01, fMaxRange); // sigma
 220   landaun->SetParLimits(4, 0.00, fMaxRange); // sigma_n
 221   // Fix the number parameter
 222   landaun->FixParameter(5, n);               // N
 223
 224   // Set the range and name of the scale parameters
 225   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 226     landaun->SetParameter(5+i-2, n == 2 ? 0.05 : 0);
 227     landaun->SetParLimits(5+i-2, 0,1);
 228     landaun->SetParName(5+i-2, Form("a_{%d}", i));
 229   }
 230
 231   // Do the fit
 232   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
 233
 234   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
 235   fFunctions.AddAtAndExpand(landaun, n-1);
 236
 237   return landaun;
 238 }
 239
 240 //====================================================================
 241 AliForwardUtil::Histos::~Histos()
 242 {
 243   if (fFMD1i) delete fFMD1i;
 244   if (fFMD2i) delete fFMD2i;
 245   if (fFMD2o) delete fFMD2o;
 246   if (fFMD3i) delete fFMD3i;
 247   if (fFMD3o) delete fFMD3o;
 248 }
 249
 250 //____________________________________________________________________
 251 TH2D*
 252 AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
 253                              const TAxis& etaAxis) const
 254 {
 255   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
 256   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r),
 257                         Form("FMD%d%c cache", d, r),
 258                         etaAxis.GetNbins(), etaAxis.GetXmin(),
 259                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
 260   hist->SetXTitle("#eta");
 261   hist->SetYTitle("#phi [radians]");
 262   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 263   hist->Sumw2();
 264   hist->SetDirectory(0);
 265
 266   return hist;
 267 }
 268 //____________________________________________________________________
 269 void
 270 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
 271 {
 272   fFMD1i = Make(1, 'I', etaAxis);
 273   fFMD2i = Make(2, 'I', etaAxis);
 274   fFMD2o = Make(2, 'O', etaAxis);
 275   fFMD3i = Make(3, 'I', etaAxis);
 276   fFMD3o = Make(3, 'O', etaAxis);
 277 }
 278 //____________________________________________________________________
 279 void
 280 AliForwardUtil::Histos::Clear(Option_t* option)
 281 {
 282   fFMD1i->Reset(option);
 283   fFMD2i->Reset(option);
 284   fFMD2o->Reset(option);
 285   fFMD3i->Reset(option);
 286   fFMD3o->Reset(option);
 287 }
 288
 289 //____________________________________________________________________
 290 TH2D*
 291 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
 292 {
 293   switch (d) {
 294   case 1: return fFMD1i;
 295   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
 296   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
 297   }
 298   return 0;
 299 }
 300 //====================================================================
 301 TList*
 302 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
 303 {
 304   if (!d) return 0;
 305   TList* list = new TList;
 306   list->SetName(fName.Data());
 307   d->Add(list);
 308   return list;
 309 }
 310 //____________________________________________________________________
 311 TList*
 312 AliForwardUtil::RingHistos::GetOutputList(TList* d) const
 313 {
 314   if (!d) return 0;
 315   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
 316   return list;
 317 }
 318
 319 //____________________________________________________________________
 320 TH1*
 321 AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
 322 {
 323   return static_cast<TH1*>(d->FindObject(name));
 324 }
 325
 326 //
 327 // EOF
 328 //