#include <TMath.h>
#include <TError.h>
+//====================================================================
+UShort_t
+AliForwardUtil::ParseCollisionSystem(const char* sys)
+{
+ TString s(sys);
+ s.ToLower();
+ if (s.Contains("p-p") || s.Contains("pp")) return AliForwardUtil::kPP;
+ if (s.Contains("pb-pb") || s.Contains("pbpb")) return AliForwardUtil::kPbPb;
+ if (s.Contains("a-a") || s.Contains("aa")) return AliForwardUtil::kPbPb;
+ return AliForwardUtil::kUnknown;
+}
+//____________________________________________________________________
+const char*
+AliForwardUtil::CollisionSystemString(UShort_t sys)
+{
+ switch (sys) {
+ case AliForwardUtil::kPP: return "pp";
+ case AliForwardUtil::kPbPb: return "PbPb";
+ }
+ return "unknown";
+}
+//____________________________________________________________________
+UShort_t
+AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t v)
+{
+ Float_t energy = v;
+ if (sys != AliForwardUtil::kPP) energy = energy / 208 * 82;
+ if (TMath::Abs(energy - 900.) < 10) return 900;
+ if (TMath::Abs(energy - 2400.) < 10) return 2400;
+ if (TMath::Abs(energy - 2750.) < 10) return 2750;
+ if (TMath::Abs(energy - 5500.) < 40) return 5500;
+ if (TMath::Abs(energy - 7000.) < 10) return 7000;
+ if (TMath::Abs(energy - 10000.) < 10) return 10000;
+ if (TMath::Abs(energy - 14000.) < 10) return 14000;
+ return 0;
+}
+//____________________________________________________________________
+const char*
+AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
+{
+ return Form("%04dGeV", cms);
+}
+//____________________________________________________________________
+Short_t
+AliForwardUtil::ParseMagneticField(Float_t v)
+{
+ if (TMath::Abs(v - 5.) < 1 ) return +5;
+ if (TMath::Abs(v + 5.) < 1 ) return -5;
+ if (TMath::Abs(v) < 1) return 0;
+ return 999;
+}
+//____________________________________________________________________
+const char*
+AliForwardUtil::MagneticFieldString(Short_t f)
+{
+ return Form("%01dkG", f);
+}
+
+
//====================================================================
Int_t AliForwardUtil::fgConvolutionSteps = 100;
Double_t AliForwardUtil::fgConvolutionNSigma = 5;
Double_t landauGaus1(Double_t* xp, Double_t* pp)
{
Double_t x = xp[0];
- Double_t constant = pp[0];
- Double_t delta = pp[1];
- Double_t xi = pp[2];
- Double_t sigma = pp[3];
- Double_t sigma_n = pp[4];
+ Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
+ Double_t delta = pp[AliForwardUtil::ELossFitter::kDelta];
+ Double_t xi = pp[AliForwardUtil::ELossFitter::kXi];
+ Double_t sigma = pp[AliForwardUtil::ELossFitter::kSigma];
+ Double_t sigma_n = pp[AliForwardUtil::ELossFitter::kSigmaN];
return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
}
Double_t landauGausN(Double_t* xp, Double_t* pp)
{
Double_t x = xp[0];
- Double_t constant = pp[0];
- Double_t delta = pp[1];
- Double_t xi = pp[2];
- Double_t sigma = pp[3];
- Double_t sigma_n = pp[4];
- Int_t n = Int_t(pp[5]);
- Double_t* a = &(pp[6]);
+ Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
+ Double_t delta = pp[AliForwardUtil::ELossFitter::kDelta];
+ Double_t xi = pp[AliForwardUtil::ELossFitter::kXi];
+ Double_t sigma = pp[AliForwardUtil::ELossFitter::kSigma];
+ Double_t sigma_n = pp[AliForwardUtil::ELossFitter::kSigmaN];
+ Int_t n = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
+ Double_t* a = &(pp[AliForwardUtil::ELossFitter::kA]);
return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
n, a);
}
+ /**
+ * Utility function to use in TF1 defintition
+ */
+ Double_t landauGausI(Double_t* xp, Double_t* pp)
+ {
+ Double_t x = xp[0];
+ Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
+ Double_t delta = pp[AliForwardUtil::ELossFitter::kDelta];
+ Double_t xi = pp[AliForwardUtil::ELossFitter::kXi];
+ Double_t sigma = pp[AliForwardUtil::ELossFitter::kSigma];
+ Double_t sigma_n = pp[AliForwardUtil::ELossFitter::kSigmaN];
+ Int_t i = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
+
+ return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
+ }
}
{
Double_t deltap = delta - xi * mpshift;
Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
- Double_t sigma1 = TMath::Sqrt(sigma2);
+ Double_t sigma1 = sigma_n == 0 ? sigma : TMath::Sqrt(sigma2);
Double_t xlow = x - fgConvolutionNSigma * sigma1;
- Double_t xhigh = x - fgConvolutionNSigma * sigma1;
+ Double_t xhigh = x + fgConvolutionNSigma * sigma1;
Double_t step = (xhigh - xlow) / fgConvolutionSteps;
Double_t sum = 0;
for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
- Double_t x1 = xlow + (i - .5) * step;
- Double_t x2 = xlow - (i - .5) * step;
+ Double_t x1 = xlow + (i - .5) * step;
+ Double_t x2 = xhigh - (i - .5) * step;
sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
return step * sum * invSq2pi / sigma1;
}
+//____________________________________________________________________
+Double_t
+AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n, Int_t i)
+{
+ Double_t delta_i = (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
+ Double_t xi_i = i * xi;
+ Double_t sigma_i = (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
+ if (sigma_i < 1e-10) {
+ // Fall back to landau
+ return AliForwardUtil::Landau(x, delta_i, xi_i);
+ }
+ return AliForwardUtil::LandauGaus(x, delta_i, xi_i, sigma_i, sigma_n);
+}
+
+//____________________________________________________________________
+Double_t
+AliForwardUtil::IdLandauGausdPar(Double_t x,
+ UShort_t par, Double_t dPar,
+ Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n,
+ Int_t i)
+{
+ if (dPar == 0) return 0;
+ Double_t dp = dPar;
+ Double_t d2 = dPar / 2;
+ Double_t delta_i = i * (delta + xi * TMath::Log(i));
+ Double_t xi_i = i * xi;
+ Double_t si = TMath::Sqrt(Double_t(i));
+ Double_t sigma_i = si*sigma;
+ Double_t y1 = 0;
+ Double_t y2 = 0;
+ Double_t y3 = 0;
+ Double_t y4 = 0;
+ switch (par) {
+ case 0:
+ y1 = ILandauGaus(x, delta_i+i*dp, xi_i, sigma_i, sigma_n, i);
+ y2 = ILandauGaus(x, delta_i+i*d2, xi_i, sigma_i, sigma_n, i);
+ y3 = ILandauGaus(x, delta_i-i*d2, xi_i, sigma_i, sigma_n, i);
+ y4 = ILandauGaus(x, delta_i-i*dp, xi_i, sigma_i, sigma_n, i);
+ break;
+ case 1:
+ y1 = ILandauGaus(x, delta_i, xi_i+i*dp, sigma_i, sigma_n, i);
+ y2 = ILandauGaus(x, delta_i, xi_i+i*d2, sigma_i, sigma_n, i);
+ y3 = ILandauGaus(x, delta_i, xi_i-i*d2, sigma_i, sigma_n, i);
+ y4 = ILandauGaus(x, delta_i, xi_i-i*dp, sigma_i, sigma_n, i);
+ break;
+ case 2:
+ y1 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*dp, sigma_n, i);
+ y2 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*d2, sigma_n, i);
+ y3 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*d2, sigma_n, i);
+ y4 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*dp, sigma_n, i);
+ break;
+ case 3:
+ y1 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+dp, i);
+ y2 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+d2, i);
+ y3 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-d2, i);
+ y4 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-dp, i);
+ break;
+ default:
+ return 0;
+ }
+
+ Double_t d0 = y1 - y4;
+ Double_t d1 = 2 * (y2 - y3);
+
+ Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
+
+ return g;
+}
+
//____________________________________________________________________
Double_t
AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
Double_t sigma, Double_t sigma_n, Int_t n,
Double_t* a)
{
- Double_t result = LandauGaus(x, delta, xi, sigma, sigma_n);
- for (Int_t i = 2; i <= n; i++) {
- Double_t delta_i = i * (delta + xi * TMath::Log(i));
- Double_t xi_i = i * xi;
- Double_t sigma_i = TMath::Sqrt(Double_t(n)*sigma);
- Double_t a_i = a[i];
- result += a_i * AliForwardUtil::LandauGaus(x, delta_i, xi_i,
- sigma_i, sigma_n);
- }
+ Double_t result = ILandauGaus(x, delta, xi, sigma, sigma_n, 1);
+ for (Int_t i = 2; i <= n; i++)
+ result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
return result;
}
+namespace {
+ const Int_t kColors[] = { kRed+1,
+ kPink+3,
+ kMagenta+2,
+ kViolet+2,
+ kBlue+1,
+ kAzure+3,
+ kCyan+1,
+ kTeal+2,
+ kGreen+2,
+ kSpring+3,
+ kYellow+2,
+ kOrange+2 };
+}
+
+//____________________________________________________________________
+TF1*
+AliForwardUtil::MakeNLandauGaus(Double_t c,
+ Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n, Int_t n,
+ Double_t* a,
+ Double_t xmin, Double_t xmax)
+{
+ Int_t npar = AliForwardUtil::ELossFitter::kN+n;
+ TF1* landaun = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
+ // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
+ landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
+ landaun->SetLineWidth(2);
+ landaun->SetNpx(500);
+ landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
+
+ // Set the initial parameters from the seed fit
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kC, c);
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta, delta);
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kXi, xi);
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma, sigma);
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n);
+ landaun->FixParameter(AliForwardUtil::ELossFitter::kN, n);
+
+ // Set the range and name of the scale parameters
+ for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
+ landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
+ landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
+ }
+ return landaun;
+}
+//____________________________________________________________________
+TF1*
+AliForwardUtil::MakeILandauGaus(Double_t c,
+ Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n, Int_t i,
+ Double_t xmin, Double_t xmax)
+{
+ Int_t npar = AliForwardUtil::ELossFitter::kN+1;
+ TF1* landaui = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
+ // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
+ landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
+ landaui->SetLineWidth(1);
+ landaui->SetNpx(500);
+ landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
+
+ // Set the initial parameters from the seed fit
+ landaui->SetParameter(AliForwardUtil::ELossFitter::kC, c);
+ landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta, delta);
+ landaui->SetParameter(AliForwardUtil::ELossFitter::kXi, xi);
+ landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma, sigma);
+ landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n);
+ landaui->FixParameter(AliForwardUtil::ELossFitter::kN, i);
+
+ return landaui;
+}
//====================================================================
AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
// Find the fit range
dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
- // Normalize peak to 1
- Double_t max = dist->GetMaximum();
- dist->Scale(1/max);
-
// Get the bin with maximum
Int_t maxBin = dist->GetMaximumBin();
Double_t maxE = dist->GetBinLowEdge(maxBin);
dist->GetXaxis()->SetRangeUser(0, fMaxRange);
// Define the function to fit
- TF1* landau1 = new TF1("landau1", landauGaus1, minE, maxEE, 5);
+ TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
// Set initial guesses, parameter names, and limits
- landau1->SetParameters(5,.5,.07,1,sigman);
+ landau1->SetParameters(1,0.5,0.07,0.1,sigman);
landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
- landau1->SetParLimits(1, minE, fMaxRange);
- landau1->SetParLimits(2, 0.00, fMaxRange);
- landau1->SetParLimits(3, 0.01, fMaxRange);
- if (sigman <= 0) landau1->FixParameter(4, 0);
- else landau1->SetParLimits(4, 0, fMaxRange);
+ landau1->SetNpx(500);
+ landau1->SetParLimits(kDelta, minE, fMaxRange);
+ landau1->SetParLimits(kXi, 0.00, fMaxRange);
+ landau1->SetParLimits(kSigma, 0.01, 0.1);
+ if (sigman <= 0) landau1->FixParameter(kSigmaN, 0);
+ else landau1->SetParLimits(kSigmaN, 0, fMaxRange);
// Do the fit, getting the result object
TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
-
+ landau1->SetRange(minE, fMaxRange);
fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
fFunctions.AddAtAndExpand(landau1, 0);
}
// Get some parameters from seed fit
- Double_t delta1 = r->Parameter(1);
- Double_t xi1 = r->Parameter(2);
+ Double_t delta1 = r->Parameter(kDelta);
+ Double_t xi1 = r->Parameter(kXi);
Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
Double_t minE = f->GetXmin();
+ // Array of weights
+ TArrayD a(n-1);
+ for (UShort_t i = 2; i <= n; i++)
+ a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
// Make the fit function
- TF1* landaun = new TF1(Form("landau%d", n), &landauGausN,minE,maxEi,5+n);
- landaun->SetLineStyle((n % 10)+1);
- landaun->SetLineWidth(2);
- landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
-
- // Set the initial parameters from the seed fit
- landaun->SetParameter(0, r->Parameter(0)); // Constant
- landaun->SetParameter(1, r->Parameter(1)); // Delta
- landaun->SetParameter(2, r->Parameter(2)); // xi
- landaun->SetParameter(3, r->Parameter(3)); // sigma
- landaun->SetParameter(4, r->Parameter(4)); // sigma_n
- landaun->SetParLimits(1, minE, fMaxRange); // Delta
- landaun->SetParLimits(2, 0.00, fMaxRange); // xi
- landaun->SetParLimits(3, 0.01, fMaxRange); // sigma
- landaun->SetParLimits(4, 0.00, fMaxRange); // sigma_n
- // Fix the number parameter
- landaun->FixParameter(5, n); // N
+ TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
+ r->Parameter(kDelta),
+ r->Parameter(kXi),
+ r->Parameter(kSigma),
+ r->Parameter(kSigmaN),
+ n,a.fArray,minE,maxEi);
+ landaun->SetParLimits(kDelta, minE, fMaxRange); // Delta
+ landaun->SetParLimits(kXi, 0.00, fMaxRange); // xi
+ landaun->SetParLimits(kSigma, 0.01, 1); // sigma
+ // Check if we're using the noise sigma
+ if (sigman <= 0) landaun->FixParameter(kSigmaN, 0);
+ else landaun->SetParLimits(kSigmaN, 0, fMaxRange);
// Set the range and name of the scale parameters
for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
- landaun->SetParameter(5+i-2, n == 2 ? 0.05 : 0);
- landaun->SetParLimits(5+i-2, 0,1);
- landaun->SetParName(5+i-2, Form("a_{%d}", i));
+ landaun->SetParLimits(kA+i-2, 0,1);
}
// Do the fit
TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
+ landaun->SetRange(minE, fMaxRange);
fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
fFunctions.AddAtAndExpand(landaun, n-1);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff