+//____________________________________________________________________
+TF1*
+AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
+{
+ //
+ // Fit a composite particle signal to the passed energy loss
+ // distribution
+ //
+ // Parameters:
+ // dist Data to fit the function to
+ // sigman If larger than zero, the initial guess of the
+ // detector induced noise. If zero or less, then this
+ // parameter is ignored in the fit (fixed at 0)
+ //
+ // Return:
+ // The function fitted to the data
+ //
+
+ // Find the fit range
+ Int_t cutBin = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
+ Int_t maxBin = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
+ dist->GetNbinsX());
+ dist->GetXaxis()->SetRange(cutBin, maxBin);
+
+ // Get the bin with maximum
+ Int_t peakBin = dist->GetMaximumBin();
+ Double_t peakE = dist->GetBinLowEdge(peakBin);
+
+ // Get the low edge
+ // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+ Int_t minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
+ Double_t minE = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+ Double_t maxE = dist->GetBinCenter(peakBin+2*fMinusBins);
+
+ // Get the range in bins and the integral of that range
+ Int_t minEb = dist->GetXaxis()->FindBin(minE);
+ Int_t maxEb = dist->GetXaxis()->FindBin(maxE);
+ Double_t intg = dist->Integral(minEb, maxEb);
+ if (intg <= 0) {
+ ::Warning("Fit1Particle",
+ "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
+ dist->GetName(), minE, maxE, minEb, maxEb, intg);
+ return 0;
+ }
+
+ // Restore the range
+ dist->GetXaxis()->SetRange(1, maxBin);
+
+ // Define the function to fit
+ TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
+
+ // Set initial guesses, parameter names, and limits
+ seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
+ seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
+ seed->SetNpx(500);
+ seed->SetParLimits(kDelta, minE, fMaxRange);
+ seed->SetParLimits(kXi, 0.00, 0.1); // Was fMaxRange - too wide
+ seed->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
+ if (sigman <= 0) seed->FixParameter(kSigmaN, 0);
+ else seed->SetParLimits(kSigmaN, 0, fMaxRange);
+
+ // Do the fit, getting the result object
+ if (fDebug)
+ ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
+ /* TFitResultPtr r = */ dist->Fit(seed, FIT_OPTIONS, "", minE, maxE);
+
+ maxE = dist->GetXaxis()->GetXmax();
+#if 1
+ TF1* comp = new TF1("composite", landauGausComposite,
+ minE, maxE, kSigma+1+2);
+ comp->SetParNames("C", "#Delta_{p}", "#xi", "#sigma",
+ "C#prime", "#xi#prime");
+ comp->SetParameters(0.8 * seed->GetParameter(kC), // 0 Primary weight
+ seed->GetParameter(kDelta), // 1 Primary Delta
+ seed->GetParameter(kDelta)/10, // 2 primary Xi
+ seed->GetParameter(kDelta)/5, // 3 primary sigma
+ 1.20 * seed->GetParameter(kC), // 5 Secondary weight
+ seed->GetParameter(kXi)); // 7 secondary Xi
+ // comp->SetParLimits(kC, minE, fMaxRange); // C
+ comp->SetParLimits(kDelta, minE, fMaxRange); // Delta
+ comp->SetParLimits(kXi, 0.00, fMaxRange); // Xi
+ comp->SetParLimits(kSigma, 1e-5, fMaxRange); // Sigma
+ // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+ comp->SetParLimits(kSigma+2, 0.00, fMaxRange); // Xi'
+#else
+ TF1* comp = new TF1("composite", landauGausComposite,
+ minE, maxE, kSigma+1+4);
+ comp->SetParNames("C", "#Delta_{p}", "#xi", "#sigma",
+ "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
+ comp->SetParameters(0.8 * seed->GetParameter(kC), // 0 Primary weight
+ seed->GetParameter(kDelta), // 1 Primary Delta
+ seed->GetParameter(kDelta)/10, // 2 primary Xi
+ seed->GetParameter(kDelta)/5, // 3 primary sigma
+ 1.20 * seed->GetParameter(kC), // 5 Secondary weight
+ seed->GetParameter(kDelta), // 6 secondary Delta
+ seed->GetParameter(kXi), // 7 secondary Xi
+ seed->GetParameter(kSigma)); // 8 secondary sigma
+ // comp->SetParLimits(kC, minE, fMaxRange); // C
+ comp->SetParLimits(kDelta, minE, fMaxRange); // Delta
+ comp->SetParLimits(kXi, 0.00, fMaxRange); // Xi
+ comp->SetParLimits(kSigma, 1e-5, fMaxRange); // Sigma
+ // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+ comp->SetParLimits(kSigma+2, minE/10, fMaxRange); // Delta
+ comp->SetParLimits(kSigma+3, 0.00, fMaxRange); // Xi
+ comp->SetParLimits(kSigma+4, 1e-6, fMaxRange); // Sigma
+#endif
+ comp->SetLineColor(kRed+1);
+ comp->SetLineWidth(3);
+
+ // Do the fit, getting the result object
+ if (fDebug)
+ ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
+ /* TFitResultPtr r = */ dist->Fit(comp, FIT_OPTIONS, "", minE, maxE);
+
+#if 0
+ TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
+ part1->SetLineColor(kGreen+1);
+ part1->SetLineWidth(4);
+ part1->SetRange(minE, maxE);
+ part1->SetParameters(comp->GetParameter(0), // C
+ comp->GetParameter(1), // Delta
+ comp->GetParameter(2), // Xi
+ comp->GetParameter(3), // sigma
+ 0);
+ part1->Save(minE,maxE,0,0,0,0);
+ dist->GetListOfFunctions()->Add(part1);
+
+ TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
+ part2->SetLineColor(kBlue+1);
+ part2->SetLineWidth(4);
+ part2->SetRange(minE, maxE);
+ part2->SetParameters(comp->GetParameter(4), // C
+ comp->GetParameter(5), // Delta
+ comp->GetParameter(6), // Xi
+ comp->GetParameter(7), // sigma
+ 0);
+ part2->Save(minE,maxE,0,0,0,0);
+ dist->GetListOfFunctions()->Add(part2);
+#endif
+ return comp;
+}