#include "AliFMDCorrDoubleHit.h"
#include "AliFMDCorrELossFit.h"
#include "AliLog.h"
+#include "AliForwardUtil.h"
#include <TH2D.h>
#include <TProfile.h>
#include <THStack.h>
fPhiLumping(4),
fDebug(0),
fCuts(),
- fRecalculateEta(false),
fRecalculatePhi(false),
fMinQuality(10),
fCache(),
fDoTiming(false),
- fHTiming(0)
+ fHTiming(0),
+ fMaxOutliers(0.05),
+ fOutlierCut(0.50)
{
//
// Constructor
fPhiLumping(4),
fDebug(0),
fCuts(),
- fRecalculateEta(false),
fRecalculatePhi(false),
fMinQuality(10),
fCache(),
fDoTiming(false),
- fHTiming(0)
+ fHTiming(0),
+ fMaxOutliers(0.05),
+ fOutlierCut(0.50)
{
//
// Constructor
fPhiLumping(o.fPhiLumping),
fDebug(o.fDebug),
fCuts(o.fCuts),
- fRecalculateEta(o.fRecalculateEta),
fRecalculatePhi(o.fRecalculatePhi),
fMinQuality(o.fMinQuality),
fCache(o.fCache),
fDoTiming(o.fDoTiming),
- fHTiming(o.fHTiming)
+ fHTiming(o.fHTiming),
+ fMaxOutliers(o.fMaxOutliers),
+ fOutlierCut(o.fOutlierCut)
{
//
// Copy constructor
fEtaLumping = o.fEtaLumping;
fPhiLumping = o.fPhiLumping;
fCuts = o.fCuts;
- fRecalculateEta = o.fRecalculateEta;
fRecalculatePhi = o.fRecalculatePhi;
fMinQuality = o.fMinQuality;
fCache = o.fCache;
fDoTiming = o.fDoTiming;
fHTiming = o.fHTiming;
+ fMaxOutliers = o.fMaxOutliers;
+ fOutlierCut = o.fOutlierCut;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
+namespace {
+ Double_t Rng2Cut(UShort_t d, Char_t r, Int_t xbin, TH2* h)
+ {
+ Double_t ret = 1024;
+ if (xbin < 1 && xbin > h->GetXaxis()->GetNbins()) return ret;
+ Int_t ybin = 0;
+ switch(d) {
+ case 1: ybin = 1; break;
+ case 2: ybin = (r=='i' || r=='I') ? 2 : 3; break;
+ case 3: ybin = (r=='i' || r=='I') ? 4 : 5; break;
+ default: return ret;
+ }
+ ret = h->GetBinContent(xbin,ybin);
+ return ret;
+ }
+}
+
//____________________________________________________________________
Double_t
AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
- Bool_t errors) const
+ Bool_t /*errors*/) const
{
//
// Get the multiplicity cut. If the user has set fMultCut (via
// Return:
// Lower cut on multiplicity
//
- return fCuts.GetMultCut(d,r,ieta,errors);
+ return Rng2Cut(d, r, ieta, fLowCuts);
+ // return fCuts.GetMultCut(d,r,ieta,errors);
}
//____________________________________________________________________
Double_t
AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
- Bool_t errors) const
+ Bool_t /*errors*/) const
{
//
// Get the multiplicity cut. If the user has set fMultCut (via
// Return:
// Lower cut on multiplicity
//
- return fCuts.GetMultCut(d,r,eta,errors);
+ Int_t ieta = fLowCuts->GetXaxis()->FindBin(eta);
+ return Rng2Cut(d, r, ieta, fLowCuts);
+ // return fCuts.GetMultCut(d,r,eta,errors);
}
-
+
+#ifndef NO_TIMING
+# define START_TIMER(T) if (fDoTiming) T.Start(true)
+# define GET_TIMER(T,V) if (fDoTiming) V = T.CpuTime()
+# define ADD_TIMER(T,V) if (fDoTiming) V += T.CpuTime()
+#else
+# define START_TIMER(T) do {} while (false)
+# define GET_TIMER(T,V) do {} while (false)
+# define ADD_TIMER(T,V) do {} while (false)
+#endif
+
//____________________________________________________________________
Bool_t
AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
// Copy to cache : fraction of sum 3.9% of total 2.2%
// Poisson calculation : fraction of sum 18.7% of total 10.6%
// Diagnostics : fraction of sum 3.7% of total 2.1%
- Double_t reEtaTime = 0;
Double_t nPartTime = 0;
Double_t corrTime = 0;
Double_t rePhiTime = 0;
Double_t copyTime = 0;
Double_t poissonTime= 0;
Double_t diagTime = 0;
- if (fDoTiming) totalT.Start(true);
+ START_TIMER(totalT);
Double_t etaCache[20*512]; // Same number of strips per ring
Double_t phiCache[20*512]; // whether it is inner our outer.
Double_t oldPhi = phi;
// --- Re-calculate eta - needed for satelittes ------------
- if (fDoTiming) timer.Start(true);
- if (eta == AliESDFMD::kInvalidEta || fRecalculateEta)
- eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
- if (fDoTiming) reEtaTime += timer.CpuTime();
+ START_TIMER(timer);
etaCache[s*nt+t] = eta;
// --- Check this strip ------------------------------------
rh->fGood->Fill(eta);
// --- If we asked to re-calculate phi for (x,y) IP --------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
if (fRecalculatePhi) {
oldPhi = phi;
phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
}
phiCache[s*nt+t] = phi;
- if (fDoTiming) rePhiTime += timer.CpuTime();
+ ADD_TIMER(timer,rePhiTime);
// --- Apply phi corner correction to eloss ----------------
if (fUsePhiAcceptance == kPhiCorrectELoss)
d, r, s, t, eta);
// --- Now caluculate Nch for this strip using fits --------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
Double_t n = 0;
if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
rh->fELoss->Fill(mult);
// rh->fEvsN->Fill(mult,n);
// rh->fEtaVsN->Fill(eta, n);
- if (fDoTiming) nPartTime += timer.CpuTime();
+ ADD_TIMER(timer,nPartTime);
// --- Calculate correction if needed ----------------------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
// Temporary stuff - remove Correction call
Double_t c = 1;
if (fUsePhiAcceptance == kPhiCorrectNch)
c = AcceptanceCorrection(r,t);
// Double_t c = Correction(d,r,t,eta,lowFlux);
- if (fDoTiming) corrTime += timer.CpuTime();
+ ADD_TIMER(timer,corrTime);
fCorrections->Fill(c);
if (c > 0) n /= c;
// rh->fEvsM->Fill(mult,n);
rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
// --- Make a copy and reset as needed -------------------------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
TH2D* hclone = fCache.Get(d,r);
// hclone->Reset();
// TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
// hclone->Add(h);
h->Reset();
}
- if (fDoTiming) copyTime += timer.CpuTime();
+ ADD_TIMER(timer,copyTime);
// --- Store Poisson result ------------------------------------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
TH2D* poisson = rh->fPoisson.Result();
for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
Double_t eta = etaCache[s*nt+t];
// Double_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
// Double_t eta = fmd.Eta(d,r,s,t);
- // if (fRecalculateEta)
- // eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
- // if (fRecalculatePhi)
- // phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
if (fUsePoisson) {
h->Fill(eta,phi,poissonV);
rh->fDensity->Fill(eta, phi, poissonV);
hclone->Fill(eta,phi,poissonV);
}
}
- if (fDoTiming) poissonTime += timer.CpuTime();
+ ADD_TIMER(timer,poissonTime);
// --- Make diagnostics - eloss vs poisson ---------------------
- if (fDoTiming) timer.Start(true);
+ START_TIMER(timer);
Int_t nY = h->GetNbinsY();
+ Int_t nIn = 0; // Count non-outliers
+ Int_t nOut = 0; // Count outliers
for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
// Set the overflow bin to contain the phi acceptance
Double_t phiAcc = rh->fGood->GetBinContent(ieta);
eLossV = h->GetBinContent(ieta,iphi);
}
- rh->fELossVsPoisson->Fill(eLossV, poissonV);
- rh->fDiffELossPoisson->Fill(poissonV < 1e-12 ? 0 :
- (eLossV - poissonV) / poissonV);
+ if (poissonV < 1e-12 && eLossV < 1e-12)
+ // we do not care about trivially empty bins
+ continue;
+ Bool_t outlier = CheckOutlier(eLossV, poissonV, fOutlierCut);
+ Double_t rel = eLossV < 1e-12 ? 0 : (poissonV - eLossV) / eLossV;
+ if (outlier) {
+ rh->fELossVsPoissonOut->Fill(eLossV, poissonV);
+ rh->fDiffELossPoissonOut->Fill(rel);
+ nOut++;
}
- }
- if (fDoTiming) diagTime += timer.CpuTime();
+ else {
+ rh->fELossVsPoisson->Fill(eLossV, poissonV);
+ rh->fDiffELossPoisson->Fill(rel);
+ nIn++;
+ } // if (outlier)
+ } // for (iphi)
+ } // for (ieta)
+ Int_t nTotal = (nIn+nOut);
+ Double_t outRatio = (nTotal > 0 ? Double_t(nOut) / nTotal : 0);
+ rh->fOutliers->Fill(outRatio);
+ if (outRatio < fMaxOutliers) rh->fPoisson.FillDiagnostics();
+ else h->SetBit(AliForwardUtil::kSkipRing);
+ ADD_TIMER(timer,diagTime);
// delete hclone;
} // for q
} // for d
if (fDoTiming) {
- fHTiming->Fill(1,reEtaTime);
+ // fHTiming->Fill(1,reEtaTime);
fHTiming->Fill(2,nPartTime);
fHTiming->Fill(3,corrTime);
fHTiming->Fill(4,rePhiTime);
return kTRUE;
}
+//_____________________________________________________________________
+Bool_t
+AliFMDDensityCalculator::CheckOutlier(Double_t eloss,
+ Double_t poisson,
+ Double_t cut) const
+{
+ if (eloss < 1e-6) return true;
+ Double_t diff = TMath::Abs(poisson - eloss) / eloss;
+ return diff > cut;
+}
//_____________________________________________________________________
Int_t
AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
AliFMDCorrELossFit::ELossFit::fgLeastWeight,
fMaxParticles);
}
+//_____________________________________________________________________
+Int_t
+AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
+ UShort_t d, Char_t r, Double_t eta) const
+{
+ //
+ // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
+ //
+ // Parameters:
+ // cor Correction
+ // d Detector
+ // r Ring
+ // eta Eta
+ //
+ DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
+ if(!cor) return -1;
+
+ AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,eta, -1);
+ if (!fit) {
+ // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
+ return -1;
+ }
+ return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError,
+ AliFMDCorrELossFit::ELossFit::fgLeastWeight,
+ fMaxParticles);
+}
//_____________________________________________________________________
void
AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
const AliFMDCorrELossFit* cor = fcm.GetELossFit();
cor->CacheBins(fMinQuality);
+ cor->Print(fDebug > 5 ? "RCS" : "C");
TAxis eta(axis.GetNbins(),
axis.GetXmin(),
fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
-
+
for (Int_t i = 0; i < nEta; i++) {
+ Double_t leta = fMaxWeights->GetXaxis()->GetBinCenter(i+1);
Double_t w[5];
- w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
- w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
- w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
- w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
- w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
- Double_t l[5];
- l[0] = GetMultCut(1, 'I', i+1, false);
- l[1] = GetMultCut(2, 'I', i+1, false);
- l[2] = GetMultCut(2, 'O', i+1, false);
- l[3] = GetMultCut(3, 'I', i+1, false);
- l[4] = GetMultCut(3, 'O', i+1, false);
- for (Int_t j = 0; j < 5; j++) {
+ w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', leta);
+ w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', leta);
+ w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', leta);
+ w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', leta);
+ w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', leta);
+ for (Int_t j = 0; j < 5; j++)
if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
- if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
- }
}
+
+ // Cache cuts in histogram
+ fCuts.FillHistogram(fLowCuts);
}
//_____________________________________________________________________
d->Add(fMaxWeights);
d->Add(fLowCuts);
- // TNamed* sigma = new TNamed("sigma",
- // (fIncludeSigma ? "included" : "excluded"));
- TObject* maxP = AliForwardUtil::MakeParameter("maxParticle", fMaxParticles);
- TObject* method = AliForwardUtil::MakeParameter("method", fUsePoisson);
- TObject* phiA = AliForwardUtil::MakeParameter("phiAcceptance",
- fUsePhiAcceptance);
- TObject* etaL = AliForwardUtil::MakeParameter("etaLumping", fEtaLumping);
- TObject* phiL = AliForwardUtil::MakeParameter("phiLumping", fPhiLumping);
- TObject* reEt = AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta);
- TObject* rePh = AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi);
-
TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
nFiles->SetMergeMode('+');
// d->Add(sigma);
- d->Add(maxP);
- d->Add(method);
- d->Add(phiA);
- d->Add(etaL);
- d->Add(phiL);
- d->Add(reEt);
- d->Add(rePh);
+ d->Add(AliForwardUtil::MakeParameter("maxParticle", fMaxParticles));
+ d->Add(AliForwardUtil::MakeParameter("minQuality", fMinQuality));
+ d->Add(AliForwardUtil::MakeParameter("method", fUsePoisson));
+ d->Add(AliForwardUtil::MakeParameter("phiAcceptance",fUsePhiAcceptance));
+ d->Add(AliForwardUtil::MakeParameter("etaLumping", fEtaLumping));
+ d->Add(AliForwardUtil::MakeParameter("phiLumping", fPhiLumping));
+ d->Add(AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi));
+ d->Add(AliForwardUtil::MakeParameter("maxOutliers", fMaxOutliers));
+ d->Add(AliForwardUtil::MakeParameter("outlierCut", fOutlierCut));
d->Add(nFiles);
// d->Add(nxi);
fCuts.Output(d,"lCuts");
}
PFV("Max(particles)", fMaxParticles );
- PFV("Poisson method", fUsePoisson );
+ PFV("Min(quality)", fMinQuality );
+ PFB("Poisson method", fUsePoisson );
PFV("Eta lumping", fEtaLumping );
PFV("Phi lumping", fPhiLumping );
- PFV("Recalculate eta", fRecalculateEta );
- PFV("Recalculate phi", fRecalculatePhi );
- PFV("Use phi acceptance", phiM);
+ PFB("Recalculate phi", fRecalculatePhi );
+ PFB("Use phi acceptance", phiM);
PFV("Lower cut", "");
fCuts.Print();
fDensity(0),
fELossVsPoisson(0),
fDiffELossPoisson(0),
+ fELossVsPoissonOut(0),
+ fDiffELossPoissonOut(0),
+ fOutliers(0),
fPoisson(),
fELoss(0),
fELossUsed(0),
fDensity(0),
fELossVsPoisson(0),
fDiffELossPoisson(0),
+ fELossVsPoissonOut(0),
+ fDiffELossPoissonOut(0),
+ fOutliers(0),
fPoisson("ignored"),
fELoss(0),
fELossUsed(0),
fDensity->SetYTitle("#phi [radians]");
fDensity->SetZTitle("Inclusive N_{ch} density");
+ // --- Create increasing sized bins --------------------------------
+ TArrayD bins;
+ // bins, lowest order, higest order, return array
+ const char* nchP = "N_{ch}^{Poisson}";
+ const char* nchE = "N_{ch}^{#Delta}";
+ AliForwardUtil::MakeLogScale(300, 0, 2, bins);
fELossVsPoisson = new TH2D("elossVsPoisson",
- "N_{ch} from energy loss vs from Poission",
- 500, 0, 100, 500, 0, 100);
+ "N_{ch} from energy loss vs from Poisson",
+ bins.GetSize()-1, bins.GetArray(),
+ bins.GetSize()-1, bins.GetArray());
fELossVsPoisson->SetDirectory(0);
- fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
- fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
+ fELossVsPoisson->SetXTitle(nchE);
+ fELossVsPoisson->SetYTitle(nchP);
fELossVsPoisson->SetZTitle("Correlation");
+ fELossVsPoissonOut =
+ static_cast<TH2D*>(fELossVsPoisson
+ ->Clone(Form("%sOutlier",
+ fELossVsPoisson->GetName())));
+ fELossVsPoissonOut->SetDirectory(0);
+ fELossVsPoissonOut->SetMarkerStyle(20);
+ fELossVsPoissonOut->SetMarkerSize(0.3);
+ fELossVsPoissonOut->SetMarkerColor(kBlack);
+ fELossVsPoissonOut->SetTitle(Form("%s for outliers",
+ fELossVsPoisson->GetTitle()));
fDiffELossPoisson = new TH1D("diffElossPoisson",
- "(N_{ch,#DeltaE}-N_{ch,Poisson})/N_{ch,Poisson}",
+ Form("(%s-%s)/%s", nchP, nchE, nchE),
100, -1, 1);
fDiffELossPoisson->SetDirectory(0);
fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
fDiffELossPoisson->SetFillStyle(3001);
fDiffELossPoisson->Sumw2();
+ fDiffELossPoissonOut =
+ static_cast<TH1D*>(fDiffELossPoisson
+ ->Clone(Form("%sOutlier",fDiffELossPoisson->GetName())));
+ fDiffELossPoissonOut->SetDirectory(0);
+ fDiffELossPoissonOut->SetTitle(Form("%s for outliers",
+ fDiffELossPoisson->GetTitle()));
+ fDiffELossPoissonOut->SetMarkerColor(Color()-2);
+ fDiffELossPoissonOut->SetFillColor(Color()-2);
+ fDiffELossPoissonOut->SetFillStyle(3002);
+
+ fOutliers = new TH1D("outliers", "Fraction of outliers", 100, 0, 1);
+ fOutliers->SetDirectory(0);
+ fOutliers->SetXTitle("N_{outlier}/(N_{outlier}+N_{inside})");
+ fOutliers->SetYTitle("#sum_{events}#sum_{bins}");
+ fOutliers->SetFillColor(Color());
+ fOutliers->SetFillStyle(3001);
+ fOutliers->SetLineColor(kBlack);
+
fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
640, -1, 15);
fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
fELoss->SetFillStyle(3003);
fELoss->SetLineColor(kBlack);
fELoss->SetLineStyle(2);
- fELoss->SetLineWidth(2);
+ fELoss->SetLineWidth(1);
fELoss->SetDirectory(0);
fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
fDensity(o.fDensity),
fELossVsPoisson(o.fELossVsPoisson),
fDiffELossPoisson(o.fDiffELossPoisson),
+ fELossVsPoissonOut(o.fELossVsPoissonOut),
+ fDiffELossPoissonOut(o.fDiffELossPoissonOut),
+ fOutliers(o.fOutliers),
fPoisson(o.fPoisson),
fELoss(o.fELoss),
fELossUsed(o.fELossUsed),
fDensity = static_cast<TH2D*>(o.fDensity->Clone());
fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
+ fELossVsPoissonOut = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
+ fDiffELossPoissonOut = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
+ fOutliers = static_cast<TH1D*>(o.fOutliers->Clone());
fPoisson = o.fPoisson;
fELoss = static_cast<TH1D*>(o.fELoss->Clone());
fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
d->Add(fCorr);
d->Add(fDensity);
d->Add(fELossVsPoisson);
+ d->Add(fELossVsPoissonOut);
d->Add(fDiffELossPoisson);
+ d->Add(fDiffELossPoissonOut);
+ d->Add(fOutliers);
fPoisson.Output(d);
fPoisson.GetOccupancy()->SetFillColor(Color());
fPoisson.GetMean()->SetFillColor(Color());
d->Add(fPhiBefore);
d->Add(fPhiAfter);
- Bool_t inner = (fRing == 'I' || fRing == 'i');
- Int_t nStr = inner ? 512 : 256;
- Int_t nSec = inner ? 20 : 40;
- TAxis x(nStr, -.5, nStr-.5);
- TAxis y(nSec, -.5, nSec-.5);
+ TAxis x(NStrip(), -.5, NStrip()-.5);
+ TAxis y(NSector(), -.5, NSector()-.5);
x.SetTitle("strip");
y.SetTitle("sector");
fPoisson.Define(x, y);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff