//
#include "AliFMDCorrELossFit.h"
#include "AliForwardUtil.h"
+#include "AliLandauGaus.h"
#include <TF1.h>
#include <TGraph.h>
#include <TBrowser.h>
#include <THStack.h>
#include <TLatex.h>
#include <TLegend.h>
+#include <TLine.h>
#include <TH1D.h>
#include <AliLog.h>
#include <TMath.h>
}
//____________________________________________________________________
AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
- : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
- Int_t(f.GetParameter(AliForwardUtil::ELossFitter::kN)) :
+ : fN(f.GetNpar() > AliLandauGaus::kN ?
+ Int_t(f.GetParameter(AliLandauGaus::kN)) :
1),
fNu(f.GetNDF()),
fChi2(f.GetChisquare()),
- fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
- fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
- fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
- fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
- fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
+ fC(f.GetParameter(AliLandauGaus::kC)),
+ fDelta(f.GetParameter(AliLandauGaus::kDelta)),
+ fXi(f.GetParameter(AliLandauGaus::kXi)),
+ fSigma(f.GetParameter(AliLandauGaus::kSigma)),
+ fSigmaN(f.GetParameter(AliLandauGaus::kSigmaN)),
fA(0),
- fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
- fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
- fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
- fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
- fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
+ fEC(f.GetParError(AliLandauGaus::kC)),
+ fEDelta(f.GetParError(AliLandauGaus::kDelta)),
+ fEXi(f.GetParError(AliLandauGaus::kXi)),
+ fESigma(f.GetParError(AliLandauGaus::kSigma)),
+ fESigmaN(f.GetParError(AliLandauGaus::kSigmaN)),
fEA(0),
fQuality(quality),
fDet(0),
fA = new Double_t[fN];
fEA = new Double_t[fN];
for (Int_t i = 0; i < fN-1; i++) {
- fA[i] = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
- fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
+ fA[i] = f.GetParameter(AliLandauGaus::kA+i);
+ fEA[i] = f.GetParError(AliLandauGaus::kA+i);
}
fA[fN-1] = -9999;
fEA[fN-1] = -9999;
// \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
// @f]
//
- // (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
+ // (see AliLandauGaus::NLandauGaus) for the maximum @f$ N @f$
// that fulfills the requirements
//
// Parameters:
// Return:
// @f$ f_N(x;\Delta,\xi,\sigma')@f$
//
- return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
- TMath::Min(maxN, UShort_t(fN)), fA);
+ return AliLandauGaus::Fn(x, fDelta, fXi, fSigma, fSigmaN,
+ TMath::Min(maxN, UShort_t(fN)), fA);
}
//____________________________________________________________________
//
// If the denominator is zero, then 1 is returned.
//
- // See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
+ // See also AliLandauGaus::Fi and AliLandauGaus::NLandauGaus
// for more information on the evaluated functions.
//
// Parameters:
for (Int_t i = 1; i <= n; i++) {
Double_t a = (i == 1 ? 1 : fA[i-1]);
if (fA[i-1] < 0) break;
- Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
+ Double_t f = AliLandauGaus::Fi(x,fDelta,fXi,fSigma,fSigmaN,i);
num += i * a * f;
den += a * f;
}
<< (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
<< " Quality: " << fQuality << "\n"
<< " NParticles: " << fN << " (" << FindMaxWeight() << ")\n"
- << OUTPAR("Delta", fDelta, fEDelta)
- << OUTPAR("xi", fXi, fEXi)
- << OUTPAR("sigma", fSigma, fESigma)
+ << OUTPAR("Delta", fDelta, fEDelta)
+ << OUTPAR("xi", fXi, fEXi)
+ << OUTPAR("sigma", fSigma, fESigma)
<< OUTPAR("sigma_n", fSigmaN, fESigmaN);
for (Int_t i = 0; i < fN-1; i++)
std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
TF1* ret = 0;
if (i <= 0)
- ret = AliForwardUtil::MakeNLandauGaus(fC * 1, fDelta, fXi,
- fSigma, fSigmaN, maxW/*fN*/,
- fA, lowX, upX);
+ ret = AliLandauGaus::MakeFn(fC * 1, fDelta, fXi,
+ fSigma, fSigmaN, maxW/*fN*/, fA, lowX, upX);
+ else if (i == 1)
+ ret = AliLandauGaus::MakeF1(fC, fDelta, fXi, fSigma, fSigmaN, lowX, upX);
else if (i <= maxW)
- ret = AliForwardUtil::MakeILandauGaus(fC*(i == 1 ? 1 : fA[i-2]),
- fDelta, fXi,
- fSigma, fSigmaN, i, lowX, upX);
-
+ ret = AliLandauGaus::MakeFi(fC*(i == 1 ? 1 : fA[i-2]),
+ fDelta, fXi, fSigma, fSigmaN, i, lowX, upX);
+
return ret;
}
//____________________________________________________________________
TF1* f = 0;
TGraph* g = 0;
try {
- if (!(f = GetF1()))
+ if (!(f = GetF1(1,5))) // First landau up to Delta/Delta_{mip}=5
throw TString("Didn't TF1 object");
if (!(g = new TGraph(f, "i")))
throw TString("Failed to integrate function");
bool good = false;
bool vals = false;
bool legd = false;
+ bool peak = false;
if (opt.Contains("COMP")) {
opt.ReplaceAll("COMP","");
comp = true;
if (!opt.Contains("SAME")) {
gPad->Clear();
}
-
+ if (opt.Contains("PEAK")) {
+ peak = true;
+ }
TLegend* l = 0;
if (legd) {
l = new TLegend(.3, .5, .59, .94);
}
TObjArray cleanup;
Int_t maxW = FindMaxWeight();
- TF1* tot = AliForwardUtil::MakeNLandauGaus(fC * 1,
- fDelta, fXi,
- fSigma, fSigmaN,
- maxW/*fN*/, fA,
- 0.01, 10);
+ TF1* tot = AliLandauGaus::MakeFn(fC * 1, fDelta, fXi, fSigma, fSigmaN,
+ maxW/*fN*/, fA, 0.01, 10);
tot->SetLineColor(kBlack);
tot->SetLineWidth(2);
tot->SetLineStyle(1);
}
+ if (peak) {
+ TLine* pl = new TLine(fDelta, 0.01*max, fDelta, 1.5*max);
+ pl->SetLineStyle(2);
+ pl->SetLineColor(kBlack);
+ pl->Draw();
+ }
if (!comp) {
gPad->cd();
return;
opt.Append(" same");
for (Int_t i=1; i <= fN; i++) {
if (good && i > maxW) break;
- TF1* f = AliForwardUtil::MakeILandauGaus(fC*(i == 1 ? 1 : fA[i-2]),
+ TF1* f = AliLandauGaus::MakeFi(fC*(i == 1 ? 1 : fA[i-2]),
fDelta, fXi,
fSigma, fSigmaN,
i, 0.01, 10);
return *this;
}
#define CACHE(BIN,IDX,OFFSET) fCache[IDX*OFFSET+BIN-1]
+#define CACHEDR(BIN,D,R,OFFSET) \
+ CACHE(BIN,(D == 1 ? 0 : (D - 2) * 2 + 1 + (R=='I' || R=='i' ? 0 : 1)),OFFSET)
//____________________________________________________________________
void
AliFMDCorrELossFit::CacheBins(UShort_t minQuality) const
{
+ AliLandauGaus::EnableSigmaShift(TestBit(kHasShift));
if (fCache.GetSize() > 0) return;
Int_t nRings = fRings.GetEntriesFast();
}
DMSG(fDebug, 10, "Got ringArray %s for FMD%d%c", ringArray->GetName(), d, r);
if (fCache.GetSize() <= 0) CacheBins(minQ);
+#if 0
Int_t idx = (d == 1 ? 0 :
(d - 2) * 2 + 1 + (r=='I' || r=='i' ? 0 : 1));
- Int_t bin = CACHE(etabin, idx, fEtaAxis.GetNbins());
+#endif
+ Int_t bin = CACHEDR(etabin, d, r, fEtaAxis.GetNbins());
if (bin < 0 || bin > ringArray->GetEntriesFast()) return 0;
TList*
AliFMDCorrELossFit::GetStacks(Bool_t err,
Bool_t rel,
- Bool_t /*good*/,
+ Bool_t good,
UShort_t maxN) const
{
// Get a list of THStacks
if (!fRings.At(i)) continue;
TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
Int_t nFits = a->GetEntriesFast();
- Int_t color = AliForwardUtil::RingColor(IDX2DET(i), IDX2RING(i));
-
- TH1D* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
- TH1D* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
- TH1D* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
- TH1D* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
- TH1D* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
+ UShort_t d = IDX2DET(i);
+ Char_t r = IDX2RING(i);
+ Int_t color = AliForwardUtil::RingColor(d, r);
+
+ TH1* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
+ TH1* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
+ TH1* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
+ TH1* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
+ TH1* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
// TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
- TH1D* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
- TH1D* hA[maxN];
+ TH1* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
+ TH1* hA[maxN];
const char* ho = (rel || !err ? "hist" : "e");
sChi2nu->Add(hChi, "hist"); // 0
sC ->Add(hC, ho); // 1
hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
static_cast<THStack*>(stacks->At(kN+k))->Add(hA[k-1], ho);
}
-
- for (Int_t j = 0; j < nFits; j++) {
- ELossFit* f = static_cast<ELossFit*>(a->At(j));
- if (!f) continue;
-
- Int_t b = f->fBin;
- Double_t vChi2nu = (rel ? f->fQuality
- : (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
- Double_t vC = (rel ? (f->fC >0 ?f->fEC /f->fC :0)
- : f->fC);
- Double_t vDelta = (rel ? (f->fDelta >0 ?f->fEDelta /f->fDelta :0)
- : f->fDelta);
- Double_t vXi = (rel ? (f->fXi >0 ?f->fEXi /f->fXi :0)
- : f->fXi);
- Double_t vSigma = (rel ? (f->fSigma >0 ?f->fESigma /f->fSigma :0)
- : f->fSigma);
- Int_t nW = (rel ? CACHE(j,i,fEtaAxis.GetNbins()) :
- f->FindMaxWeight());
- // Double_t sigman = (rel ? (f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0)
- // : f->SigmaN);
- hChi ->SetBinContent(b, vChi2nu);
- hN ->SetBinContent(b, nW);
- hC ->SetBinContent(b, vC);
- hDelta ->SetBinContent(b, vDelta);
- hXi ->SetBinContent(b, vXi);
- hSigma ->SetBinContent(b, vSigma);
- // if (vChi2nu > 1e-12) {
- // min[kChi2nu] = TMath::Min(min[kChi2nu], vChi2nu);
- // max[kChi2nu] = TMath::Max(max[kChi2nu], vChi2nu);
- // }
- // if (vC > 1e-12) {
- // min[kC] = TMath::Min(min[kC], vC);
- // max[kC] = TMath::Max(max[kC], vC);
- // }
- // if (vDelta > 1e-12) {
- // min[kDelta] = TMath::Min(min[kDelta], vDelta);
- // max[kDelta] = TMath::Max(max[kDelta], vDelta);
- // }
- // if (vXi > 1e-12) {
- // min[kXi] = TMath::Min(min[kXi], vXi);
- // max[kXi] = TMath::Max(max[kXi], vXi);
- // }
- // if (vSigma > 1e-12) {
- // min[kSigma] = TMath::Min(min[kSigma], vSigma);
- // max[kSigma] = TMath::Max(max[kSigma], vSigma);
- // }
- // if (nW > 1e-12) {
- // min[kN] = TMath::Min(min[kN], Double_t(nW));
- // max[kN] = TMath::Max(max[kN], Double_t(nW));
- // }
- // hSigmaN->SetBinContent(b, sigman);
- if (!rel) {
- hC ->SetBinError(b, f->fEC);
- hDelta ->SetBinError(b, f->fEDelta);
- hXi ->SetBinError(b, f->fEXi);
- hSigma ->SetBinError(b, f->fESigma);
- // hSigmaN->SetBinError(b, f->fESigmaN);
+
+ if (good) {
+ for (Int_t j = 1; j <= fEtaAxis.GetNbins(); j++) {
+ ELossFit* f = FindFit(d, r, j, 20);
+ if (!f) continue;
+
+ UpdateStackHist(f, rel, j, hChi, hN, hC, hDelta, hXi, hSigma, maxN, hA);
}
- for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
- Double_t vA = (rel ? (f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0)
- : f->fA[k]);
- hA[k]->SetBinContent(b, vA);
- if (!rel) hA[k]->SetBinError(b, f->fEA[k]);
- // if (vA > 1e-12) {
- // min[kN+1+k] = TMath::Min(min[kN+1+k], vA);
- // max[kN+1+k] = TMath::Max(max[kN+1+k], vA);
- // }
+ }
+ else {
+ for (Int_t j = 0; j < nFits; j++) {
+ ELossFit* f = static_cast<ELossFit*>(a->At(j));
+ if (!f) continue;
+
+ UpdateStackHist(f, rel, CACHE(j,i,fEtaAxis.GetNbins()),
+ hChi, hN, hC, hDelta, hXi, hSigma, maxN, hA);
}
}
}
return stacks;
}
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::UpdateStackHist(ELossFit* f, Bool_t rel,
+ Int_t used,
+ TH1* hChi, TH1* hN,
+ TH1* hC, TH1* hDelta,
+ TH1* hXi, TH1* hSigma,
+ Int_t maxN, TH1** hA) const
+{
+ Int_t b =f->fBin;
+ Double_t chi2nu =(rel ? f->fQuality : (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
+ Double_t c =(rel ? (f->fC >0 ?f->fEC /f->fC :0) : f->fC);
+ Double_t delta =(rel ? (f->fDelta >0 ?f->fEDelta /f->fDelta :0) : f->fDelta);
+ Double_t xi =(rel ? (f->fXi >0 ?f->fEXi /f->fXi :0) : f->fXi);
+ Double_t sigma =(rel ? (f->fSigma >0 ?f->fESigma /f->fSigma :0) : f->fSigma);
+ Int_t w =(rel ? used : f->FindMaxWeight());
+ // Double_t sigman = (rel ? (f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0)
+ // : f->SigmaN);
+ hChi ->SetBinContent(b, chi2nu);
+ hN ->SetBinContent(b, w);
+ hC ->SetBinContent(b, c);
+ hDelta ->SetBinContent(b, delta);
+ hXi ->SetBinContent(b, xi);
+ hSigma ->SetBinContent(b, sigma);
+
+ if (!rel) {
+ hC ->SetBinError(b, f->fEC);
+ hDelta ->SetBinError(b, f->fEDelta);
+ hXi ->SetBinError(b, f->fEXi);
+ hSigma ->SetBinError(b, f->fESigma);
+ // hSigmaN->SetBinError(b, f->fESigmaN);
+ }
+ for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
+ Double_t a = (rel ? (f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0) : f->fA[k]);
+ hA[k]->SetBinContent(b, a);
+ if (!rel) hA[k]->SetBinError(b, f->fEA[k]);
+ }
+}
//____________________________________________________________________
void
THStack* stack = static_cast<THStack*>(stacks->At(i));
-
+ if (!stack->GetHists() || stack->GetHists()->GetEntries() <= 0) {
+ AliWarningF("No histograms in %s", stack->GetName());
+ continue;
+ }
// Double_t powMax = TMath::Log10(max[i]);
// Double_t powMin = min[i] <= 0 ? powMax : TMath::Log10(min[i]);
// if (powMax-powMin > 2. && min[i] != 0) p->SetLogy();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff