Added missing correction object classes.
Removed dependence on old forward analysis code altogether - still need correction creators
}
// --- Make the task and add it to the manager ---------------------
- AliForwardMultiplicity* task = new AliForwardMultiplicity("FMD");
+ AliForwardMultiplicityTask* task = new AliForwardMultiplicityTask("FMD");
mgr->AddTask(task);
// --- Set parameters on the algorithms ----------------------------
+ // Whether to enable low flux specific code
+ task->SetEnableLowFlux(kFALSE);
// Set the number of SPD tracklets for which we consider the event a
// low flux event
task->GetEventInspector().SetLowFluxCut(1000);
// Set whether to use increasing bin sizes
task->GetEnergyFitter().SetUseIncreasingBins(true);
// Set whether to do fit the energy distributions
- task->GetEnergyFitter().SetDoFits(kTRUE);
+ task->GetEnergyFitter().SetDoFits(kFALSE);
+ // Set whether to make the correction object
+ task->GetEnergyFitter().SetDoMakeObject(kFALSE);
// Set the low cut used for energy
task->GetEnergyFitter().SetLowCut(0.4);
// Set the number of bins to subtract from maximum of distributions
// Set the minimum number of entries in the distribution before
// trying to fit to the data
task->GetEnergyFitter().SetMinEntries(1000);
- // Set the low cut used for sharing
- task->GetSharingFilter().SetLowCut(0.3);
+ // Set the low cut used for sharing - overrides settings in eloss fits
+ // task->GetSharingFilter().SetLowCut(0.4);
+ // Set the number of xi's (width of landau peak) to stop at
+ task->GetSharingFilter().SetNXi(1);
+ // Set the maximum number of particle to try to reconstruct
+ task->GetDensityCalculator().SetMaxParticles(2);
+ // Set the lower multiplicity cut. Overrides setting in energy loss fits.
+ // task->GetDensityCalculator().SetMultCut(0.4);
// Set the number of extra bins (beyond the secondary map border)
task->GetHistCollector().SetNCutBins(1);
// Set the correction cut, that is, when bins in the secondary map
// Set the overall debug level (1: some output, 3: a lot of output)
task->SetDebug(0);
// Set the debug level of a single algorithm
- task->GetEnergyFitter().SetDebug(3);
+ // task->GetEventInspector().SetDebug(4);
+ // --- Set limits on fits the energy -------------------------------
+ // Maximum relative error on parameters
+ AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
+ // Least weight to use
+ AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
+ // Maximum value of reduced chi^2
+ AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 5;
+
- // --- Set up the parameer manager ---------------------------------
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // --- Set up the parameter manager ---------------------------------
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
AliMCEventHandler* mcHandler =
dynamic_cast<AliMCEventHandler*>(mgr->GetMCtruthEventHandler());
Info("AddTaskFMD", "MC handler %p", mcHandler);
- if(mcHandler) {
- pars->SetRealData(kFALSE);
- pars->SetProcessPrimary(kTRUE);
- pars->SetProcessHits(kFALSE);
- }
- else {
- pars->SetRealData(kTRUE);
- pars->SetProcessPrimary(kFALSE);
- pars->SetProcessHits(kFALSE);
- }
- pars->Init();
+ // if(mcHandler) {
+ // pars->SetRealData(kFALSE);
+ // pars->SetProcessPrimary(kTRUE);
+ // pars->SetProcessHits(kFALSE);
+ // }
+ // else {
+ // pars->SetRealData(kTRUE);
+ // pars->SetProcessPrimary(kFALSE);
+ // pars->SetProcessHits(kFALSE);
+ // }
+ // pars->Init();
- // --- Makek the output container and connect it -------------------
+ // --- Make the output container and connect it --------------------
TString outputfile = AliAnalysisManager::GetCommonFileName();
- outputfile += Form(":%s",pars->GetDndetaAnalysisName());
+ // outputfile += ":PWG2forwardDnDeta";
+ // Form(":%s",pars->GetDndetaAnalysisName());
AliAnalysisDataContainer* histOut =
mgr->CreateContainer("Forward", TList::Class(),
AliAnalysisManager::kOutputContainer,outputfile);
--- /dev/null
+#include "AliFMDCorrDoubleHit.h"
+#include <TBrowser.h>
+#include <TH1D.h>
+#include <AliLog.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliFMDCorrDoubleHit::AliFMDCorrDoubleHit()
+ : fCorrections()
+{
+ fCorrections.SetOwner(kTRUE);
+ fCorrections.SetName("doubleHit");
+}
+//____________________________________________________________________
+AliFMDCorrDoubleHit::AliFMDCorrDoubleHit(const AliFMDCorrDoubleHit& o)
+ : TObject(o),
+ fCorrections(o.fCorrections)
+{
+}
+//____________________________________________________________________
+AliFMDCorrDoubleHit::~AliFMDCorrDoubleHit()
+{
+ fCorrections.Clear();
+}
+//____________________________________________________________________
+AliFMDCorrDoubleHit&
+AliFMDCorrDoubleHit::operator=(const AliFMDCorrDoubleHit& o)
+{
+ fCorrections = o.fCorrections;
+
+ return *this;
+}
+//____________________________________________________________________
+TH1D*
+AliFMDCorrDoubleHit::GetCorrection(UShort_t d, Char_t r) const
+{
+ Int_t idx = GetRingIndex(d, r);
+ if (idx < 0) return 0;
+
+ TObject* o = fCorrections.At(idx);
+ if (!o) {
+ AliWarning(Form("No double hit correction found for FMD%d%c", d, r));
+ return 0;
+ }
+ return static_cast<TH1D*>(o);
+}
+//____________________________________________________________________
+Int_t
+AliFMDCorrDoubleHit::GetRingIndex(UShort_t d, Char_t r) const
+{
+ switch (d) {
+ case 1: return 0;
+ case 2: return (r == 'I' || r == 'i' ? 1 : 2); break;
+ case 3: return (r == 'I' || r == 'i' ? 3 : 4); break;
+ }
+ AliWarning(Form("Index for FMD%d%c not found", d, r));
+ return -1;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrDoubleHit::SetCorrection(UShort_t d, Char_t r, TH1D* h)
+{
+ Int_t idx = GetRingIndex(d, r);
+ if (idx < 0) return kFALSE;
+
+ h->SetName(Form("FMD%d%c", d, r));
+ h->SetTitle(Form("Double hit correction for FMD%d%c", d,r));
+ h->SetXTitle("#eta");
+ h->SetYTitle("#sum_{i} N_{i,strips hit}(#eta)/"
+ "#sum_{i} N_{i,total hits}(#eta)");
+ // h->SetFillColor(Color(d,r));
+ h->SetFillStyle(3001);
+ h->SetDirectory(0);
+ h->SetStats(0);
+
+ fCorrections.AddAtAndExpand(h, idx);
+ return kTRUE;
+}
+//____________________________________________________________________
+void
+AliFMDCorrDoubleHit::Browse(TBrowser* b)
+{
+ b->Add(&fCorrections);
+}
+//____________________________________________________________________
+void
+AliFMDCorrDoubleHit::Print(Option_t* option) const
+{
+ std::cout << "Double hit correction" << std::endl;
+ fCorrections.Print(option);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+//
+// This class contains the secondary correction and the double hit
+// correction used in low-flux events.
+//
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRDOUBLEHIT_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRDOUBLEHIT_H
+#include <TObject.h>
+#include <TObjArray.h>
+class TH1D;
+
+/**
+ * This class contains double hit correction used in low-flux events.
+ *
+ *
+ * The double hit correction is given by
+ * @f[
+ * h_{r}(\eta) = \frac{\sum_i N_{i,strips hit}(\eta)}{
+ * \sum_i N_{i,total hits}(\eta)}
+ * @f]
+ *
+ * where @f$ N_{i,strips hit}(\eta)@f$ is the number of strips in the
+ * @f$\eta@f$ bin that had one or more hits in event @f$i@f$, and
+ * @f$N_{i,total hits}(\eta)@f$ is the total number hits in the
+ * @f$\eta@f$ bin.
+ *
+ * These are generated from Monte-Carlo truth information.
+ *
+ * @ingroup pwg2_forward_corr
+ *
+ */
+class AliFMDCorrDoubleHit : public TObject
+{
+public:
+ /**
+ * Default constructor
+ */
+ AliFMDCorrDoubleHit();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrDoubleHit(const AliFMDCorrDoubleHit& o);
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliFMDCorrDoubleHit();
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrDoubleHit& operator=(const AliFMDCorrDoubleHit& o);
+ /**
+ * @{
+ * @name Get corrections and parameters
+ */
+ /**
+ * Get the double hit correction @f$ h_{r}(\eta)@f$
+ *
+ * @param d Detector number
+ * @param r Ring identifier
+ *
+ * @return @f$ h_{r}(\eta)@f$
+ */
+ TH1D* GetCorrection(UShort_t d, Char_t r) const;
+ /* @} */
+
+ /**
+ * @{
+ * @name Set corrections and parameters
+ */
+ /**
+ * Set the double hit correction @f$ h_{r}(\eta)@f$. Note, that the
+ * object takes ownership of the passed pointer.
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param h @f$ h_{r}(\eta)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(UShort_t d, Char_t r, TH1D* h);
+ /* @} */
+
+ /**
+ * @{
+ * @name Auxiliary member functions
+ */
+ /**
+ * Declare this as a folder
+ *
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return true; }
+ /**
+ * Browse this object in the browser
+ *
+ * @param b
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Print this object
+ *
+ * @param option
+ */
+ void Print(Option_t* option="R") const; //*MENU*
+ /* @} */
+protected:
+ /**
+ * Get the index corresponding to the given ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Index (0 based) or negative in case of errors
+ */
+ Int_t GetRingIndex(UShort_t d, Char_t r) const;
+
+ TObjArray fCorrections; // Array of per-ring double hit corr.
+ ClassDef(AliFMDCorrDoubleHit,1); //
+};
+#endif
+// Local Variables:
+// mode: C++
+// End:
--- /dev/null
+#include "AliFMDCorrELossFit.h"
+#include "AliForwardUtil.h"
+#include <TF1.h>
+#include <TBrowser.h>
+#include <TVirtualPad.h>
+#include <THStack.h>
+#include <TH1D.h>
+#include <AliLog.h>
+#include <TList.h>
+#include <iostream>
+#include <iomanip>
+
+Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
+Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
+Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 20;
+
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit::ELossFit()
+ : fN(0),
+ fNu(0),
+ fChi2(0),
+ fC(0),
+ fDelta(0),
+ fXi(0),
+ fSigma(0),
+ fSigmaN(0),
+ fA(0),
+ fEC(0),
+ fEDelta(0),
+ fEXi(0),
+ fESigma(0),
+ fESigmaN(0),
+ fEA(0),
+ fQuality(0),
+ fDet(0),
+ fRing('\0'),
+ fBin(0)
+{}
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
+ : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
+ f.GetParameter(AliForwardUtil::ELossFitter::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)),
+ 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)),
+ fEA(0),
+ fQuality(quality),
+ fDet(0),
+ fRing('\0'),
+ fBin(0)
+{
+ if (fN <= 0) return;
+ 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[fN-1] = -9999;
+ fEA[fN-1] = -9999;
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality,UShort_t n,
+ Double_t chi2, UShort_t nu,
+ Double_t c, Double_t ec,
+ Double_t delta, Double_t edelta,
+ Double_t xi, Double_t exi,
+ Double_t sigma, Double_t esigma,
+ Double_t sigman, Double_t esigman,
+ Double_t* a, Double_t* ea)
+ : fN(n),
+ fNu(nu),
+ fChi2(chi2),
+ fC(c),
+ fDelta(delta),
+ fXi(xi),
+ fSigma(sigma),
+ fSigmaN(sigman),
+ fA(0),
+ fEC(ec),
+ fEDelta(edelta),
+ fEXi(exi),
+ fESigma(esigma),
+ fESigmaN(esigman),
+ fEA(0),
+ fQuality(quality),
+ fDet(0),
+ fRing('\0'),
+ fBin(0)
+{
+ if (fN <= 0) return;
+ fA = new Double_t[fN];
+ fEA = new Double_t[fN];
+ for (Int_t i = 0; i < fN-1; i++) {
+ fA[i] = a[i];
+ fEA[i] = ea[i];
+ }
+ fA[fN-1] = -9999;
+ fEA[fN-1] = -9999;
+}
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
+ : TObject(o),
+ fN(o.fN),
+ fNu(o.fNu),
+ fChi2(o.fChi2),
+ fC(o.fC),
+ fDelta(o.fDelta),
+ fXi(o.fXi),
+ fSigma(o.fSigma),
+ fSigmaN(o.fSigmaN),
+ fA(0),
+ fEC(o.fEC),
+ fEDelta(o.fEDelta),
+ fEXi(o.fEXi),
+ fESigma(o.fESigma),
+ fESigmaN(o.fESigmaN),
+ fEA(0),
+ fQuality(o.fQuality),
+ fDet(o.fDet),
+ fRing(o.fRing),
+ fBin(o.fBin)
+{
+ if (fN <= 0) return;
+ fA = new Double_t[fN];
+ fEA = new Double_t[fN];
+ for (Int_t i = 0; i < fN-1; i++) {
+ fA[i] = o.fA[i];
+ fEA[i] = o.fEA[i];
+ }
+ fA[fN-1] = -9999;
+ fEA[fN-1] = -9999;
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit&
+AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
+{
+ fN = o.fN;
+ fNu = o.fNu;
+ fChi2 = o.fChi2;
+ fC = o.fC;
+ fDelta = o.fDelta;
+ fXi = o.fXi;
+ fSigma = o.fSigma;
+ fSigmaN = o.fSigmaN;
+ fEC = o.fEC;
+ fEDelta = o.fEDelta;
+ fEXi = o.fEXi;
+ fESigma = o.fESigma;
+ fESigmaN = o.fESigmaN;
+ fQuality = o.fQuality;
+ fDet = o.fDet;
+ fRing = o.fRing;
+ fBin = o.fBin;
+ if (fA) delete [] fA;
+ if (fEA) delete [] fEA;
+ fA = 0;
+ fEA = 0;
+
+ if (fN <= 0) return *this;
+ fA = new Double_t[fN];
+ fEA = new Double_t[fN];
+ for (Int_t i = 0; i < fN; i++) {
+ fA[i] = o.fA[i];
+ fEA[i] = o.fEA[i];
+ }
+
+ return *this;
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit::~ELossFit()
+{
+ if (fA) delete[] fA;
+ if (fEA) delete[] fEA;
+}
+
+
+//____________________________________________________________________
+Int_t
+AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError,
+ Double_t leastWeight,
+ UShort_t maxN) const
+{
+ Int_t n = TMath::Min(maxN, UShort_t(fN-1));
+ Int_t m = 1;
+ // fN is one larger than we have data
+ for (Int_t i = 0; i < n-1; i++, m++) {
+ if (fA[i] < leastWeight) break;
+ if (fEA[i] / fA[i] > maxRelError) break;
+ }
+ return m;
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x,
+ UShort_t maxN) const
+{
+ return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
+ TMath::Min(maxN, UShort_t(fN)), fA);
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x,
+ UShort_t maxN) const
+{
+ UShort_t n = TMath::Min(maxN, UShort_t(fN-1));
+ Double_t num = 0;
+ Double_t den = 0;
+ 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);
+ num += i * a * f;
+ den += a * f;
+ }
+ if (den <= 0) return 1;
+ return num / den;
+}
+
+
+#define OUTPAR(N,V,E) \
+ std::setprecision(9) << \
+ std::setw(12) << N << ": " << \
+ std::setw(14) << V << " +/- " << \
+ std::setw(14) << E << " (" << \
+ std::setprecision(-1) << \
+ std::setw(5) << 100*(V>0?E/V:1) << "%)\n"
+
+
+//____________________________________________________________________
+Int_t
+AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
+{
+ const ELossFit* other = static_cast<const ELossFit*>(o);
+ if (this->fQuality < other->fQuality) return -1;
+ if (this->fQuality > other->fQuality) return +1;
+ Double_t chi2nu = (fNu == 0 ? 1000 : fChi2 / fNu);
+ Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
+ if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
+ if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
+ return 0;
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
+{
+ std::cout << GetName() << ":\n"
+ << " chi^2/nu = " << fChi2 << "/" << fNu << " = "
+ << (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("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]);
+ std::cout << std::flush;
+}
+
+//____________________________________________________________________
+const Char_t*
+AliFMDCorrELossFit::ELossFit::GetName() const
+{
+ return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
+{
+ // Draw this one
+ Draw(b ? b->GetDrawOption() : "comp");
+ gPad->SetLogy();
+ gPad->Update();
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
+{
+ TString opt(option);
+ opt.ToUpper();
+ bool comp = false;
+ if (opt.Contains("COMP")) {
+ opt.ReplaceAll("COMP","");
+ comp = true;
+ }
+ if (!opt.Contains("SAME")) {
+ gPad->Clear();
+ }
+
+ TObjArray cleanup;
+ TF1* tot = AliForwardUtil::MakeNLandauGaus(1,
+ fDelta, fXi,
+ fSigma, fSigmaN,
+ fN, fA,
+ 0.01, 10);
+ tot->SetLineColor(kBlack);
+ tot->SetLineWidth(2);
+ tot->SetLineStyle(1);
+ tot->SetTitle(GetName());
+ Double_t max = tot->GetMaximum();
+
+ if (!opt.Contains("SAME")) {
+ TH1* frame = new TH1F(GetName(),
+ Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
+ 100, 0, 10);
+ frame->SetMinimum(max/10000);
+ frame->SetMaximum(max*1.4);
+ frame->SetXTitle("#Delta / #Delta_{mip}");
+ frame->Draw();
+ opt.Append(" SAME");
+ }
+ tot->DrawCopy(opt.Data());
+ cleanup.Add(tot);
+
+ if (!comp) {
+ gPad->cd();
+ return;
+ }
+
+ Double_t min = max;
+ opt.Append(" same");
+ Int_t maxW = FindMaxWeight();
+ for (Int_t i=1; i <= fN; i++) {
+ TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]),
+ fDelta, fXi,
+ fSigma, fSigmaN,
+ i, 0.01, 10);
+ f->SetLineWidth(2);
+ f->SetLineStyle(i > maxW ? 2 : 1);
+ min = TMath::Min(f->GetMaximum(), min);
+ f->DrawCopy(opt.Data());
+ cleanup.Add(f);
+ }
+ min /= 100;
+ tot->GetHistogram()->SetMaximum(max);
+ tot->GetHistogram()->SetMinimum(min);
+ tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
+
+ gPad->cd();
+}
+
+
+//____________________________________________________________________
+#define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
+ Double_t maxRelError,
+ Double_t leastWeight)
+{
+ Int_t qual = 0;
+ if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
+ if (CHECKPAR(fDelta, fEDelta, maxRelError)) qual++;
+ if (CHECKPAR(fXi, fEXi, maxRelError)) qual++;
+ if (CHECKPAR(fSigma, fESigma, maxRelError)) qual++;
+ if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
+ qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
+ fQuality = qual;
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::AliFMDCorrELossFit()
+ : TObject(),
+ fRings(),
+ fEtaAxis(0,0,0),
+ fLowCut(0)
+{
+ fRings.SetOwner(kTRUE);
+ fEtaAxis.SetTitle("#eta");
+ fEtaAxis.SetName("etaAxis");
+ fRings.SetName("rings");
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
+ : TObject(o),
+ fRings(o.fRings),
+ fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
+ fLowCut(0)
+{
+ fEtaAxis.SetTitle("#eta");
+ fEtaAxis.SetName("etaAxis");
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit::~AliFMDCorrELossFit()
+{
+ fRings.Clear();
+}
+
+//____________________________________________________________________
+AliFMDCorrELossFit&
+AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
+{
+ fRings = o.fRings;
+ fLowCut = o.fLowCut;
+ SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
+
+ return *this;
+}
+//____________________________________________________________________
+Int_t
+AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
+{
+ if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
+ AliWarning("No eta axis defined");
+ return -1;
+ }
+ Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
+ if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
+ return bin;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
+{
+ TObjArray* ringArray = GetOrMakeRingArray(d, r);
+ if (!ringArray) {
+ AliError(Form("Failed to make ring array for FMD%d%c", d, r));
+ return kFALSE;
+ }
+ if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
+ AliError(Form("bin=%d is out of range [%d,%d]",
+ etaBin, 1, fEtaAxis.GetNbins()));
+ return kFALSE;
+ }
+ // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
+ ringArray->AddAtAndExpand(fit, etaBin);
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
+{
+ Int_t bin = FindEtaBin(eta);
+ if (bin <= 0) {
+ AliError(Form("eta=%f is out of range [%f,%f]",
+ eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
+ return kFALSE;
+ }
+
+ return SetFit(d, r, bin, fit);
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r,
+ Double_t eta,
+ Int_t quality,UShort_t n,
+ Double_t chi2, UShort_t nu,
+ Double_t c, Double_t ec,
+ Double_t delta, Double_t edelta,
+ Double_t xi, Double_t exi,
+ Double_t sigma, Double_t esigma,
+ Double_t sigman, Double_t esigman,
+ Double_t* a, Double_t* ea)
+{
+ ELossFit* e = new ELossFit(quality, n,
+ chi2, nu,
+ c, ec,
+ delta, edelta,
+ xi, exi,
+ sigma, esigma,
+ sigman, esigman,
+ a, ea);
+ if (!SetFit(d, r, eta, e)) {
+ delete e;
+ return kFALSE;
+ }
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta,
+ Int_t quality, const TF1& f)
+{
+ ELossFit* e = new ELossFit(quality, f);
+ if (!SetFit(d, r, eta, e)) {
+ delete e;
+ return kFALSE;
+ }
+ return kTRUE;
+}
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit*
+AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Int_t etabin) const
+{
+ TObjArray* ringArray = GetRingArray(d, r);
+ if (!ringArray) {
+ AliError(Form("Failed to make ring array for FMD%d%c", d, r));
+ return 0;
+ }
+ if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
+ AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
+ etabin, 1, fEtaAxis.GetNbins(), d, r));
+ return 0;
+ }
+ if (etabin > ringArray->GetEntriesFast()) {
+ AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
+ etabin, 1, ringArray->GetEntriesFast(), d, r));
+ return 0;
+ }
+ else if (etabin >= ringArray->GetEntriesFast()) {
+ AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, "
+ "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
+ etabin-1));
+ etabin--;
+ }
+ else if (!ringArray->At(etabin)) {
+ AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d",
+ etabin, d, r, etabin+1));
+ etabin++;
+ }
+ return static_cast<ELossFit*>(ringArray->At(etabin));
+}
+//____________________________________________________________________
+AliFMDCorrELossFit::ELossFit*
+AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Double_t eta) const
+{
+ Int_t etabin = FindEtaBin(eta);
+ return FindFit(d, r, etabin);
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
+ case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
+ }
+ if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
+ return static_cast<TObjArray*>(fRings.At(idx));
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
+ case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
+ }
+ if (idx < 0) return 0;
+ if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
+ TObjArray* a = new TObjArray(0);
+ // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
+ a->SetName(Form("FMD%d%c", d, r));
+ a->SetOwner();
+ fRings.AddAtAndExpand(a, idx);
+ }
+ return static_cast<TObjArray*>(fRings.At(idx));
+}
+
+namespace {
+ TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
+ Int_t color)
+ {
+ TH1D* h = new TH1D(Form("%s_%s", name, title),
+ Form("%s %s", name, title),
+ axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
+ h->SetDirectory(0);
+ h->SetMarkerStyle(20);
+ h->SetMarkerColor(color);
+ h->SetMarkerSize(0.5);
+ h->SetFillColor(color);
+ h->SetFillStyle(3001);
+ h->SetLineColor(color);
+ return h;
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::Draw(Option_t* option)
+{
+ TString opt(Form("nostack %s", option));
+ opt.ToLower();
+ Bool_t rel = (opt.Contains("rel"));
+ Bool_t err = (opt.Contains("err"));
+ if (rel) opt.ReplaceAll("rel","");
+ if (err) opt.ReplaceAll("err","");
+ Int_t nRings = fRings.GetEntriesFast();
+ UShort_t maxN = 0;
+ for (Int_t i = 0; i < nRings; i++) {
+ if (!fRings.At(i)) continue;
+ TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
+ Int_t nFits = a->GetEntriesFast();
+
+ for (Int_t j = 0; j < nFits; j++) {
+ ELossFit* fit = static_cast<ELossFit*>(a->At(j));
+ if (!fit) continue;
+ maxN = TMath::Max(maxN, UShort_t(fit->fN));
+ }
+ }
+ // AliInfo(Form("Maximum N is %d", maxN));
+ Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
+ TVirtualPad* pad = gPad;
+ pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
+
+ THStack* chi2nu;
+ THStack* c;
+ THStack* delta;
+ THStack* xi;
+ THStack* sigma;
+ THStack* sigman;
+ THStack* n;
+ TList stacks;
+ stacks.AddAt(chi2nu= new THStack("chi2", "#chi^{2}/#nu"), 0);
+ stacks.AddAt(c = new THStack("c", "C"), 1);
+ stacks.AddAt(delta = new THStack("delta", "#Delta_{mp}"), 2);
+ stacks.AddAt(xi = new THStack("xi", "#xi"), 3);
+ stacks.AddAt(sigma = new THStack("sigma", "#sigma"), 4);
+ stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"), 5);
+ stacks.AddAt(n = new THStack("n", "N"), 6);
+ for (Int_t i = 1; i <= maxN; i++) {
+ stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
+ }
+
+ for (Int_t i = 0; i < nRings; i++) {
+ if (!fRings.At(i)) continue;
+ TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
+ Int_t nFits = a->GetEntriesFast();
+ Int_t color = 0;
+ switch (i) {
+ case 0: color = kRed+2; break;
+ case 1: color = kGreen+2; break;
+ case 2: color = kGreen-2; break;
+ case 3: color = kBlue+2; break;
+ case 4: color = kBlue-2; break;
+ }
+
+ 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);
+ TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
+ TH1D* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
+ TH1D* hA[maxN];
+ const char* ho = (rel || !err ? "hist" : "e");
+ chi2nu->Add(hChi, "hist"); // 0
+ c ->Add(hC, ho); // 1
+ delta ->Add(hDelta, ho); // 2
+ xi ->Add(hXi, ho); // 3
+ sigma ->Add(hSigma, ho); // 4
+ sigman->Add(hSigmaN,ho); // 5
+ n ->Add(hN, "hist"); // 6
+ hChi->SetFillColor(color);
+ hChi->SetFillStyle(3001);
+ hN->SetFillColor(color);
+ hN->SetFillStyle(3001);
+
+ for (Int_t k = 1; k <= maxN; k++) {
+ hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
+ static_cast<THStack*>(stacks.At(6+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;
+ Int_t nW = f->FindMaxWeight();
+ hChi ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
+ hN ->SetBinContent(b, nW);
+
+ if (rel) {
+ hC ->SetBinContent(b, f->fC >0 ?f->fEC /f->fC :0);
+ hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta :0);
+ hXi ->SetBinContent(b, f->fXi >0 ?f->fEXi /f->fXi :0);
+ hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma :0);
+ hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
+ }
+ else {
+ hC ->SetBinContent(b, f->fC);
+ hDelta ->SetBinContent(b, f->fDelta);
+ hXi ->SetBinContent(b, f->fXi);
+ hSigma ->SetBinContent(b, f->fSigma);
+ hSigmaN->SetBinContent(b, f->fSigmaN);
+ 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++) {
+ if (rel)
+ hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
+ else {
+ hA[k]->SetBinContent(b, f->fA[k]);
+ hA[k]->SetBinError(b, f->fEA[k]);
+ }
+ }
+ }
+ }
+ Int_t nPad2 = (nPad+1) / 2;
+ for (Int_t i = 0; i < nPad; i++) {
+ Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
+ TVirtualPad* p = pad->cd(iPad);
+ p->SetLeftMargin(.15);
+ p->SetFillColor(0);
+ p->SetFillStyle(0);
+ p->SetGridx();
+ p->SetGridy();
+ if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
+
+ THStack* stack = static_cast<THStack*>(stacks.At(i));
+ // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
+ stack->Draw(opt.Data());
+
+ TString tit(stack->GetTitle());
+ if (rel && i != 0 && i != 6)
+ tit = Form("#delta %s/%s", tit.Data(), tit.Data());
+ TH1* hist = stack->GetHistogram();
+ TAxis* yaxis = hist->GetYaxis();
+ yaxis->SetTitle(tit.Data());
+ yaxis->SetTitleSize(0.15);
+ yaxis->SetLabelSize(0.08);
+ yaxis->SetTitleOffset(0.35);
+ yaxis->SetNdivisions(5);
+
+ TAxis* xaxis = stack->GetHistogram()->GetXaxis();
+ xaxis->SetTitle("#eta");
+ xaxis->SetTitleSize(0.15);
+ xaxis->SetLabelSize(0.08);
+ xaxis->SetTitleOffset(0.35);
+ xaxis->SetNdivisions(10);
+
+
+ if (!rel) {
+ switch (i) {
+ case 0: break; // chi^2/nu
+ case 1: break; // C
+ case 2: stack->SetMinimum(0.4); break; // Delta
+ case 3: stack->SetMinimum(0.02);break; // xi
+ case 4: stack->SetMinimum(0.05);break; // sigma
+ case 5: break; // sigmaN
+ case 6:
+ stack->SetMinimum(-.1);
+ stack->SetMaximum(maxN+.5);
+ break; // N
+ default: break; // Weights
+ }
+ }
+ stack->DrawClone(opt.Data());
+ }
+ pad->cd();
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::Print(Option_t* option) const
+{
+ TString opt(option);
+ opt.ToUpper();
+ Int_t nRings = fRings.GetEntriesFast();
+ bool recurse = opt.Contains("R");
+
+ std::cout << "Low cut in fit range: " << fLowCut << "\n"
+ << "Eta axis: " << fEtaAxis.GetNbins()
+ << " bins, range [" << fEtaAxis.GetXmin() << ","
+ << fEtaAxis.GetXmax() << "]" << std::endl;
+
+ for (Int_t i = 0; i < nRings; i++) {
+ if (!fRings.At(i)) continue;
+ TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
+ Int_t nFits = a->GetEntriesFast();
+
+ std::cout << a->GetName() << " [" << nFits << " entries]"
+ << (recurse ? ":\n" : "\t");
+ Int_t min = fEtaAxis.GetNbins()+1;
+ Int_t max = 0;
+ for (Int_t j = 0; j < nFits; j++) {
+ if (!a->At(j)) continue;
+
+ min = TMath::Min(j, min);
+ max = TMath::Max(j, max);
+
+ if (recurse) {
+ std::cout << "Bin # " << j << "\t";
+ ELossFit* fit = static_cast<ELossFit*>(a->At(j));
+ fit->Print(option);
+ }
+ }
+ if (!recurse)
+ std::cout << " bin range: " << std::setw(3) << min
+ << "-" << std::setw(3) << max << " " << std::setw(3)
+ << (max-min+1) << " bins" << std::endl;
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrELossFit::Browse(TBrowser* b)
+{
+ b->Add(&fRings);
+ b->Add(&fEtaAxis);
+}
+
+
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFMDCORRELOSSFIT_H
+#define ALIROOT_PWG2_FORWARD_ALIFMDCORRELOSSFIT_H
+#include <TObject.h>
+#include <TAxis.h>
+#include <TObjArray.h>
+class TF1;
+class TBrowser;
+
+/**
+ * Object holding the Energy loss fit 'correction'
+ *
+ * These are generated from Monte-Carlo or real ESDs.
+ *
+ * @ingroup pwg2_forward_corr
+ */
+class AliFMDCorrELossFit : public TObject
+{
+public:
+ /**
+ * POD structure to hold data from fits
+ *
+ */
+ struct ELossFit : public TObject
+ {
+ Int_t fN; // Number of peaks fitted
+ UShort_t fNu; // Number degrees of freedom
+ Double_t fChi2; // Chi square from fit
+ Double_t fC; // Scaling constant
+ Double_t fDelta; // Most probable value
+ Double_t fXi; // Width parameter of Landau
+ Double_t fSigma; // Sigma on folded gaussian
+ Double_t fSigmaN; // Sigma of detector noise
+ Double_t* fA; // [fN] Weights
+ Double_t fEC; // Error on C
+ Double_t fEDelta; // Error on Delta
+ Double_t fEXi; // Error on Xi
+ Double_t fESigma; // Error on sigma
+ Double_t fESigmaN;// Error on sigma (noise)
+ Double_t* fEA; // [fN] Error on weights
+ Int_t fQuality;// Assigned quality
+ UShort_t fDet; // Detector
+ Char_t fRing; // Ring
+ UShort_t fBin; // Eta bin
+
+ static Double_t fgMaxRelError; // Global default max relative error
+ static Double_t fgLeastWeight; // Global default least weight
+ static Double_t fgMaxChi2nu; // Global default maximum reduced chi^2
+ /**
+ * Default constructor
+ *
+ */
+ ELossFit();
+ /**
+ * Construct from a function
+ *
+ * @param quality Quality flag
+ * @param f Function
+ */
+ ELossFit(Int_t quality,const TF1& f);
+ /**
+ * Constructor with full parameter set
+ *
+ * @param d Detector number
+ * @param r Ring identifier
+ * @param quality Quality flag
+ * @param n @f$ N@f$ - Number of fitted peaks
+ * @param chi2 @f$ \chi^2 @f$
+ * @param nu @f$ \nu @f$ - number degrees of freedom
+ * @param c @f$ C2f$ - scale constant
+ * @param ec @f$ \delta C@f$ - error on @f$ C@f$
+ * @param delta @f$ \Delta2f$ - scale constant
+ * @param edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
+ * @param xi @f$ \xi2f$ - scale constant
+ * @param exi @f$ \delta\xi@f$ - error on @f$\xi@f$
+ * @param sigma @f$ \sigma@f$ - scale constant
+ * @param esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
+ * @param sigman @f$ \sigma_n@f$ - scale constant
+ * @param esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
+ * @param a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
+ * @f$ i=2,\ldots@f$
+ */
+ ELossFit(Int_t quality,UShort_t n,
+ Double_t chi2, UShort_t nu,
+ Double_t c, Double_t ec,
+ Double_t delta, Double_t edelta,
+ Double_t xi, Double_t exi,
+ Double_t sigma, Double_t esigma,
+ Double_t sigman, Double_t esigman,
+ Double_t* a, Double_t* ea);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ ELossFit(const ELossFit& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ ELossFit& operator=(const ELossFit& o);
+ /**
+ * Destructor
+ */
+ ~ELossFit();
+ /**
+ * @{
+ * @name Access to parameters
+ */
+ /**
+ * @return Number of peaks fitted
+ */
+ Int_t GetN() const { return fN; }
+ /**
+ * @return Number degrees of freedom
+ */
+ UShort_t GetNu() const { return fNu; }
+ /**
+ * @return Chi square from fit
+ */
+ Double_t GetChi2() const { return fChi2; }
+ /**
+ * @return Scaling constant
+ */
+ Double_t GetC() const { return fC; }
+ /**
+ * @return Most probable value
+ */
+ Double_t GetDelta() const { return fDelta; }
+ /**
+ * @return Width parameter of Landau
+ */
+ Double_t GetXi() const { return fXi; }
+ /**
+ * @return Sigma on folded gaussian
+ */
+ Double_t GetSigma() const { return fSigma; }
+ /**
+ * @return Sigma of detector noise
+ */
+ Double_t GetSigmaN() const { return fSigmaN; }
+ /**
+ * @return Weights
+ */
+ Double_t* GetAs() const { return fA; }
+ /**
+ * @return Weights
+ */
+ Double_t GetA(UShort_t i) const;
+ /**
+ * @return Error on C
+ */
+ Double_t GetEC() const { return fEC; }
+ /**
+ * @return Error on Delta
+ */
+ Double_t GetEDelta() const { return fEDelta; }
+ /**
+ * @return Error on Xi
+ */
+ Double_t GetEXi() const { return fEXi; }
+ /**
+ * @return Error on sigma
+ */
+ Double_t GetESigma() const { return fESigma; }
+ /**
+ * @return Error on sigma (noise)
+ */
+ Double_t GetESigmaN() const { return fESigmaN; }
+ /**
+ * @return Error on weights
+ */
+ Double_t* GetEAs() const { return fEA; }
+ /**
+ * @return Error on weights
+ */
+ Double_t GetEA(UShort_t i) const;
+ /**
+ * @return Assigned quality
+ */
+ Int_t GetQuality() const { return fQuality; }
+ /**
+ * @return Detector
+ */
+ UShort_t GetDet() const { return fDet; }
+ /**
+ * @return Ring
+ */
+ Char_t GetRing() const { return fRing; }
+ /**
+ * @return Eta bin
+ */
+ UShort_t GetBin() const { return fBin; }
+ /* @} */
+
+ /**
+ * @{
+ * @name Evaluation
+ */
+ /**
+ * Evaluate
+ * @f[
+ * f_N(x;\Delta,\xi,\sigma') =
+ * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
+ * @f]
+ *
+ * (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
+ * that fulfills the requirements
+ *
+ * @param x Where to evaluate
+ * @param maxN @f$ \max{N}@f$
+ *
+ * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
+ */
+ Double_t Evaluate(Double_t x,
+ UShort_t maxN=999) const;
+ /**
+ * Evaluate
+ * @f[
+ * f_W(x;\Delta,\xi,\sigma') =
+ * \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
+ * f_N(x;\Delta,\xi,\sigma')} =
+ * \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
+ * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
+ * @f]
+ * where @f$ n@f$ fulfills the requirements (see FindMaxWeight).
+ *
+ * If the denominator is zero, then 1 is returned.
+ *
+ * See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
+ * for more information on the evaluated functions.
+ *
+ * @param x Where to evaluate
+ * @param maxN @f$ \max{N}@f$
+ *
+ * @return @f$ f_W(x;\Delta,\xi,\sigma')@f$.
+ */
+ Double_t EvaluateWeighted(Double_t x,
+ UShort_t maxN=9999) const;
+ /**
+ * Find the maximum weight to use. The maximum weight is the
+ * largest i for which
+ *
+ * - @f$ i \leq \max{N}@f$
+ * - @f$ a_i > \min{a}@f$
+ * - @f$ \delta a_i/a_i > \delta_{max}@f$
+ *
+ * @param maxRelError @f$ \min{a}@f$
+ * @param leastWeight @f$ \delta_{max}@f$
+ * @param maxN @f$ \max{N}@f$
+ *
+ * @return The largest index @f$ i@f$ for which the above
+ * conditions hold. Will never return less than 1.
+ */
+ Int_t FindMaxWeight(Double_t maxRelError=2*fgMaxRelError,
+ Double_t leastWeight=fgLeastWeight,
+ UShort_t maxN=999) const;
+ /* @} */
+ /**
+ * @{
+ * @name TObject Sortable interface
+ */
+ /**
+ * Declare this object as sortable
+ *
+ * @return Always true
+ */
+ Bool_t IsSortable() const { return kTRUE; }
+ /**
+ * Compare to another ELossFit object.
+ *
+ * - +1, if this quality is better (larger) than other objects quality
+ * - -1, if this quality is worse (smaller) than other objects quality
+ * - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
+ * - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
+ * - 0 otherwise
+ *
+ * @param o Other object to compare to
+ */
+ Int_t Compare(const TObject* o) const;
+ /* @} */
+ /**
+ * @{
+ * name Auxiliary member functions
+ */
+ /**
+ * Information to standard output
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option) const; // *MENU*
+ /**
+ * Draw this fit
+ *
+ * @param option Options
+ * - COMP Draw components too
+ */
+ void Draw(Option_t* option="comp"); // *MENU*
+ /**
+ * Browse this object
+ *
+ * @param b Browser
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Get the name of this object
+ *
+ *
+ * @return
+ */
+ const Char_t* GetName() const;
+ /**
+ * Calculate the quality
+ */
+ void CalculateQuality(Double_t maxChi2nu=fgMaxChi2nu,
+ Double_t maxRelError=fgMaxRelError,
+ Double_t leastWeight=fgLeastWeight);
+ /* @} */
+ ClassDef(ELossFit,1); // Result of fit
+ };
+
+ /**
+ * Default constructor
+ */
+ AliFMDCorrELossFit();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrELossFit(const AliFMDCorrELossFit& o);
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDCorrELossFit();
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrELossFit& operator=(const AliFMDCorrELossFit& o);
+
+ /**
+ * @{
+ * @name Set fits
+ */
+ /**
+ * Set the fit parameters from a function
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta
+ * @param quality Quality flag
+ * @param f Function from fit
+ */
+ Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, Int_t quality,
+ const TF1& f);
+ /**
+ * Set the fit parameters from a function
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta
+ * @param f ELoss fit result - note, the object will take ownership
+ */
+ Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* f);
+ /**
+ * Set the fit parameters from a function
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta (bin number, 1->nBins)
+ * @param f ELoss fit result - note, the object will take ownership
+ */
+ Bool_t SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* f);
+ /**
+ * Set the fit parameters from a function
+ *
+ * @param d Detector number
+ * @param r Ring identifier
+ * @param quality Quality flag
+ * @param n @f$ N@f$ - Number of fitted peaks
+ * @param chi2 @f$ \chi^2 @f$
+ * @param nu @f$ \nu @f$ - number degrees of freedom
+ * @param c @f$ C2f$ - scale constant
+ * @param ec @f$ \delta C@f$ - error on @f$ C@f$
+ * @param delta @f$ \Delta2f$ - scale constant
+ * @param edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
+ * @param xi @f$ \xi2f$ - scale constant
+ * @param exi @f$ \delta\xi@f$ - error on @f$\xi@f$
+ * @param sigma @f$ \sigma@f$ - scale constant
+ * @param esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
+ * @param sigman @f$ \sigma_n@f$ - scale constant
+ * @param esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
+ * @param a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
+ * @f$ i=2,\ldots@f$
+ */
+ Bool_t SetFit(UShort_t d, Char_t r, Double_t eta,
+ Int_t quality,UShort_t n,
+ Double_t chi2, UShort_t nu,
+ Double_t c, Double_t ec,
+ Double_t delta, Double_t edelta,
+ Double_t xi, Double_t exi,
+ Double_t sigma, Double_t esigma,
+ Double_t sigman, Double_t esigman,
+ Double_t* a, Double_t* ea);
+ /* @} */
+
+ /**
+ * @{
+ * @name Set and get eta axis
+ */
+ /**
+ * Set the eta axis to use
+ *
+ * @param axis Eta axis
+ */
+ void SetEtaAxis(const TAxis& axis);
+ /**
+ * Set the eta axis to use
+ *
+ * @param axis Eta axis
+ */
+ void SetEtaAxis(Int_t nBins, Double_t min, Double_t max);
+ /**
+ * Get the eta axis used
+ *
+ * @return
+ */
+ const TAxis& GetEtaAxis() const { return fEtaAxis; }
+ /**
+ * Set the low cut used when fitting
+ *
+ * @param cut Cut value
+ */
+ void SetLowCut(Double_t cut) { fLowCut = cut; }
+ /**
+ * Get the low cut used when fitting
+ *
+ * @return Low cut used for fitting
+ */
+ Double_t GetLowCut() const { return fLowCut; }
+ /**
+ * Find the eta bin corresponding to the given eta
+ *
+ * @param eta Eta value
+ *
+ * @return Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
+ * eta, or 0 if out of range.
+ */
+ Int_t FindEtaBin(Double_t eta) const;
+ /* @} */
+
+ /**
+ * @{
+ * @name Find fits
+ */
+ /**
+ * Find the fit corresponding to the specified parameters
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta value
+ *
+ * @return Fit parameters or null in case of problems
+ */
+ ELossFit* FindFit(UShort_t d, Char_t r, Double_t eta) const;
+ /**
+ * Find the fit corresponding to the specified parameters
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param etabin Eta bin (1 based)
+ *
+ * @return Fit parameters or null in case of problems
+ */
+ ELossFit* FindFit(UShort_t d, Char_t r, Int_t etabin) const;
+ /* @} */
+
+ /**
+ * @{
+ * @name Miscellaneous
+ */
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or null in case of problems
+ */
+ TObjArray* GetRingArray(UShort_t d, Char_t r) const;
+ /**
+ * Signal that this is a folder
+ *
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return true; }
+ /**
+ * Browse this object
+ *
+ * @param b
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Draw this object
+ *
+ * @param option Options. Possible values are
+ * - err Plot error bars
+ */
+ void Draw(Option_t* option=""); //*MENU*
+ /**
+ * Print this object.
+ *
+ * @param option Options
+ * - R Print recursive
+ *
+ */
+ void Print(Option_t* option="R") const; //*MENU*
+ /* @} */
+protected:
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or newly created container
+ */
+ TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
+
+ TObjArray fRings; // Array of rings
+ TAxis fEtaAxis; // Eta axis used
+ Double_t fLowCut; // Low cut used when fitting
+
+ ClassDef(AliFMDCorrELossFit,2);
+};
+
+//____________________________________________________________________
+inline void
+AliFMDCorrELossFit::SetEtaAxis(Int_t nBins, Double_t min, Double_t max)
+{
+ fEtaAxis.Set(nBins, min, max);
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrELossFit::SetEtaAxis(const TAxis& e)
+{
+ fEtaAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
+}
+//____________________________________________________________________
+inline Double_t
+AliFMDCorrELossFit::ELossFit::GetA(UShort_t i) const
+{
+ if (i < 1) return 0;
+ if (i > fN) return 0;
+ if (i == 1) return 1;
+ return fA[i-2];
+}
+//____________________________________________________________________
+inline Double_t
+AliFMDCorrELossFit::ELossFit::GetEA(UShort_t i) const
+{
+ if (i < 1) return 0;
+ if (i > fN) return 0;
+ if (i == 1) return 1;
+ return fEA[i-2];
+}
+
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+#include "AliFMDCorrMergingEfficiency.h"
+#include <TBrowser.h>
+#include <TH1D.h>
+#include <AliLog.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliFMDCorrMergingEfficiency::AliFMDCorrMergingEfficiency()
+ : fRingArray(),
+ fVertexAxis(0,0,0)
+{
+ fRingArray.SetOwner(kTRUE);
+ fRingArray.SetName("rings");
+ fVertexAxis.SetName("vtxAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+
+}
+//____________________________________________________________________
+AliFMDCorrMergingEfficiency::AliFMDCorrMergingEfficiency(const
+ AliFMDCorrMergingEfficiency& o)
+ : TObject(o),
+ fRingArray(o.fRingArray),
+ fVertexAxis(o.fVertexAxis.GetNbins(), o.fVertexAxis.GetXmin(),
+ o.fVertexAxis.GetXmax())
+{
+ fVertexAxis.SetName("vtxAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+}
+//____________________________________________________________________
+AliFMDCorrMergingEfficiency::~AliFMDCorrMergingEfficiency()
+{
+ fRingArray.Clear();
+}
+//____________________________________________________________________
+AliFMDCorrMergingEfficiency&
+AliFMDCorrMergingEfficiency::operator=(const AliFMDCorrMergingEfficiency& o)
+{
+ fRingArray = o.fRingArray;
+ SetVertexAxis(o.fVertexAxis);
+
+ return *this;
+}
+//____________________________________________________________________
+TH1D*
+AliFMDCorrMergingEfficiency::GetCorrection(UShort_t d, Char_t r, Double_t v) const
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0) return 0;
+ return GetCorrection(d, r, UShort_t(b));
+}
+//____________________________________________________________________
+TH1D*
+AliFMDCorrMergingEfficiency::GetCorrection(UShort_t d, Char_t r, UShort_t b) const
+{
+ TObjArray* ringArray = GetRingArray(d, r);
+ if (!ringArray) return 0;
+
+ if (b <= 0 || b > ringArray->GetEntriesFast()) {
+ AliWarning(Form("vertex bin %d out of range [1,%d]",
+ b, ringArray->GetEntriesFast()));
+ return 0;
+ }
+
+ TObject* o = ringArray->At(b-1);
+ if (!o) {
+ AliWarning(Form("No secondary map found for FMD%d%c in vertex bin %d",
+ d,r,b));
+ return 0;
+ }
+ return static_cast<TH1D*>(o);
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDCorrMergingEfficiency::FindVertexBin(Double_t v) const
+{
+ if (fVertexAxis.GetNbins() <= 0) {
+ AliWarning("No vertex array defined");
+ return 0;
+ }
+ Int_t bin = const_cast<TAxis&>(fVertexAxis).FindBin(v);
+ if (bin <= 0 || bin > fVertexAxis.GetNbins()) {
+ AliWarning(Form("vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return 0;
+ }
+ return bin;
+}
+//____________________________________________________________________
+Int_t
+AliFMDCorrMergingEfficiency::GetRingIndex(UShort_t d, Char_t r) const
+{
+ switch (d) {
+ case 1: return 0;
+ case 2: return (r == 'I' || r == 'i' ? 1 : 2); break;
+ case 3: return (r == 'I' || r == 'i' ? 3 : 4); break;
+ }
+ AliWarning(Form("Index for FMD%d%c not found", d, r));
+ return -1;
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrMergingEfficiency::GetRingArray(UShort_t d, Char_t r) const
+{
+ Int_t idx = GetRingIndex(d,r);
+ if (idx < 0) return 0;
+
+ TObject* o = fRingArray.At(idx);
+ if (!o) {
+ AliWarning(Form("No array found for FMD%d%c", d, r));
+ return 0;
+ }
+
+ return static_cast<TObjArray*>(o);
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrMergingEfficiency::GetOrMakeRingArray(UShort_t d, Char_t r)
+{
+ Int_t idx = GetRingIndex(d,r);
+ if (idx < 0) return 0;
+
+ TObject* o = fRingArray.At(idx);
+ if (!o) {
+ TObjArray* a = new TObjArray(fVertexAxis.GetNbins());
+ a->SetName(Form("FMD%d%c", d, r));
+ a->SetOwner(kTRUE);
+ fRingArray.AddAtAndExpand(a, idx);
+ return a;
+ }
+
+ return static_cast<TObjArray*>(fRingArray.At(idx));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrMergingEfficiency::SetCorrection(UShort_t d, Char_t r,
+ UShort_t b, TH1D* h)
+{
+ TObjArray* ringArray = GetOrMakeRingArray(d, r);
+ if (!ringArray) return false;
+
+ if (b <= 0 || b > fVertexAxis.GetNbins()) {
+ AliWarning(Form("Vertex bin %3d out of range [1,%3d]",
+ b, fVertexAxis.GetNbins()));
+ return false;
+ }
+ h->SetName(Form("FMD%d%c_vtxbin%03d", d, r, b));
+ h->SetTitle(Form("Secondary map correction for FMD%d%c "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ d, r, b, fVertexAxis.GetBinLowEdge(b),
+ fVertexAxis.GetBinUpEdge(b)));
+ h->SetXTitle("#eta");
+ h->SetYTitle("dN_{ch}/d#eta / sum_i N_{ch,i}");
+ h->SetFillStyle(3001);
+ h->SetDirectory(0);
+ h->SetStats(0);
+ ringArray->AddAtAndExpand(h, b-1);
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrMergingEfficiency::SetCorrection(UShort_t d, Char_t r,
+ Double_t v, TH1D* h)
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0 || b > fVertexAxis.GetNbins()) {
+ AliWarning(Form("Vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return false;
+ }
+ return SetCorrection(d, r, UShort_t(b), h);
+}
+//____________________________________________________________________
+void
+AliFMDCorrMergingEfficiency::Browse(TBrowser* b)
+{
+ b->Add(&fRingArray);
+ b->Add(&fVertexAxis);
+}
+//____________________________________________________________________
+void
+AliFMDCorrMergingEfficiency::Print(Option_t* option) const
+{
+ std::cout << "Merging efficiency correction" << std::endl;
+ fRingArray.Print(option);
+ fVertexAxis.Print(option);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+//
+// This class contains the secondary correction and the double hit
+// correction used in low-flux events.
+//
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRMERGINGEFFICIENCY_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRMERGINGEFFICIENCY_H
+#include <TObject.h>
+#include <TObjArray.h>
+#include <TAxis.h>
+class TH1D;
+
+/**
+ * This class contains the merging efficiency correction.
+ *
+ * The secondary correction is given by
+ * @f[
+ * m_{r,v}(\eta) =
+ * @f]
+ *
+ * These are generated from Monte-Carlo truth and ESD information.
+ *
+ * @ingroup pwg2_forward_corr
+ */
+class AliFMDCorrMergingEfficiency : public TObject
+{
+public:
+ /**
+ * Default constructor
+ */
+ AliFMDCorrMergingEfficiency();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrMergingEfficiency(const AliFMDCorrMergingEfficiency& o);
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliFMDCorrMergingEfficiency();
+ /**
+ * @{
+ * @name Get corrections and parameters
+ */
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrMergingEfficiency& operator=(const AliFMDCorrMergingEfficiency& o);
+ /**
+ * Get the secondary correction @f$ c_{r,v}@f$
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ *
+ * @return The correction @f$ c_{r,v}@f$
+ */
+ TH1D* GetCorrection(UShort_t d, Char_t r, Double_t v) const;
+ /**
+ * Get the secondary correction @f$ c_{r,v}@f$
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ *
+ * @return The correction @f$ c_{r,v}@f$
+ */
+ TH1D* GetCorrection(UShort_t d, Char_t r, UShort_t b) const;
+ /**
+ * Get the vertex axis used
+ *
+ * @return vertex axis
+ */
+ const TAxis& GetVertexAxis() const { return fVertexAxis; }
+ /* @} */
+
+ /**
+ * @{
+ * @name Set corrections and parameters
+ */
+ /**
+ * Set the secondary map correction @f$ m_{r,v}(\eta)@f$.
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ * @param h @f$ m_{r,v}(\eta)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(UShort_t d, Char_t r, Double_t v, TH1D* h);
+ /**
+ * Set the secondary map correction @f$ m_{r,v}(\eta)@f$
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ * @param h @f$ m_{r,v}(\eta)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(UShort_t d, Char_t r, UShort_t b, TH1D* h);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(const TAxis& axis);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(Int_t nBins, Double_t min, Double_t max);
+ /* @} */
+
+ /**
+ * @{
+ * @name Auxiliary member functions
+ */
+ /**
+ * Declare this as a folder
+ *
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return true; }
+ /**
+ * Browse this object in the browser
+ *
+ * @param b
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Print this object
+ *
+ * @param option
+ */
+ void Print(Option_t* option="R") const; //*MENU*
+ /* @} */
+protected:
+ /**
+ * Find the vertex bin that corresponds to the passed vertex
+ *
+ * @param vertex The interaction points @f$z@f$-coordinate
+ *
+ * @return Vertex bin in @f$[1,N_{\box{vertex}}]@f$ or negative if
+ * out of range
+ */
+ Int_t FindVertexBin(Double_t vertex) const;
+ /**
+ * Get the index corresponding to the given ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Index (0 based) or negative in case of errors
+ */
+ Int_t GetRingIndex(UShort_t d, Char_t r) const;
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or null in case of problems
+ */
+ TObjArray* GetRingArray(UShort_t d, Char_t r) const;
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or newly created container
+ */
+ TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
+
+ TObjArray fRingArray; // Array of per-ring, per-vertex 2nd map
+ TAxis fVertexAxis; // The vertex axis
+ ClassDef(AliFMDCorrMergingEfficiency,1); //
+};
+
+//____________________________________________________________________
+inline void
+AliFMDCorrMergingEfficiency::SetVertexAxis(Int_t nBins, Double_t min,
+ Double_t max)
+{
+ fVertexAxis.Set(nBins, min, max);
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrMergingEfficiency::SetVertexAxis(const TAxis& e)
+{
+ fVertexAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
+}
+#endif
+// Local Variables:
+// mode: C++
+// End:
--- /dev/null
+#include "AliFMDCorrSecondaryMap.h"
+#include <TBrowser.h>
+#include <TH2D.h>
+#include <AliLog.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliFMDCorrSecondaryMap::AliFMDCorrSecondaryMap()
+ : fRingArray(),
+ fVertexAxis(0,0,0),
+ fEtaAxis(0,0,0)
+{
+ fRingArray.SetOwner(kTRUE);
+ fRingArray.SetName("rings");
+ fVertexAxis.SetName("vertexAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+ fEtaAxis.SetName("etaAxis");
+ fEtaAxis.SetTitle("#eta");
+
+}
+//____________________________________________________________________
+AliFMDCorrSecondaryMap::AliFMDCorrSecondaryMap(const
+ AliFMDCorrSecondaryMap& o)
+ : TObject(o),
+ fRingArray(o.fRingArray),
+ fVertexAxis(o.fVertexAxis.GetNbins(), o.fVertexAxis.GetXmin(),
+ o.fVertexAxis.GetXmax()),
+ fEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(),
+ o.fEtaAxis.GetXmax())
+{
+ fVertexAxis.SetName("vertexAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+ fEtaAxis.SetName("etaAxis");
+ fEtaAxis.SetTitle("v_{z} [cm]");
+}
+//____________________________________________________________________
+AliFMDCorrSecondaryMap::~AliFMDCorrSecondaryMap()
+{
+ fRingArray.Clear();
+}
+//____________________________________________________________________
+AliFMDCorrSecondaryMap&
+AliFMDCorrSecondaryMap::operator=(const AliFMDCorrSecondaryMap& o)
+{
+ fRingArray = o.fRingArray;
+ SetVertexAxis(o.fVertexAxis);
+ SetEtaAxis(o.fEtaAxis);
+
+ return *this;
+}
+//____________________________________________________________________
+TH2D*
+AliFMDCorrSecondaryMap::GetCorrection(UShort_t d, Char_t r, Double_t v) const
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0) return 0;
+ return GetCorrection(d, r, UShort_t(b));
+}
+//____________________________________________________________________
+TH2D*
+AliFMDCorrSecondaryMap::GetCorrection(UShort_t d, Char_t r, UShort_t b) const
+{
+ TObjArray* ringArray = GetRingArray(d, r);
+ if (!ringArray) return 0;
+
+ if (b <= 0 || b > ringArray->GetEntriesFast()) {
+ AliWarning(Form("vertex bin %d out of range [1,%d]",
+ b, ringArray->GetEntriesFast()));
+ return 0;
+ }
+
+ TObject* o = ringArray->At(b-1);
+ if (!o) {
+ AliWarning(Form("No secondary map found for FMD%d%c in vertex bin %d",
+ d,r,b));
+ return 0;
+ }
+ return static_cast<TH2D*>(o);
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDCorrSecondaryMap::FindVertexBin(Double_t v) const
+{
+ if (fVertexAxis.GetNbins() <= 0) {
+ AliWarning("No vertex array defined");
+ return 0;
+ }
+ Int_t bin = const_cast<TAxis&>(fVertexAxis).FindBin(v);
+ if (bin <= 0 || bin > fVertexAxis.GetNbins()) {
+ AliWarning(Form("vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return 0;
+ }
+ return bin;
+}
+//____________________________________________________________________
+Int_t
+AliFMDCorrSecondaryMap::GetRingIndex(UShort_t d, Char_t r) const
+{
+ switch (d) {
+ case 1: return 0;
+ case 2: return (r == 'I' || r == 'i' ? 1 : 2); break;
+ case 3: return (r == 'I' || r == 'i' ? 3 : 4); break;
+ }
+ AliWarning(Form("Index for FMD%d%c not found", d, r));
+ return -1;
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrSecondaryMap::GetRingArray(UShort_t d, Char_t r) const
+{
+ Int_t idx = GetRingIndex(d,r);
+ if (idx < 0) return 0;
+
+ TObject* o = fRingArray.At(idx);
+ if (!o) {
+ AliWarning(Form("No array found for FMD%d%c", d, r));
+ return 0;
+ }
+
+ return static_cast<TObjArray*>(o);
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrSecondaryMap::GetOrMakeRingArray(UShort_t d, Char_t r)
+{
+ Int_t idx = GetRingIndex(d,r);
+ if (idx < 0) return 0;
+
+ TObject* o = fRingArray.At(idx);
+ if (!o) {
+ TObjArray* a = new TObjArray(fVertexAxis.GetNbins());
+ a->SetName(Form("FMD%d%c", d, r));
+ a->SetOwner(kTRUE);
+ fRingArray.AddAtAndExpand(a, idx);
+ return a;
+ }
+
+ return static_cast<TObjArray*>(fRingArray.At(idx));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrSecondaryMap::SetCorrection(UShort_t d, Char_t r,
+ UShort_t b, TH2D* h)
+{
+ TObjArray* ringArray = GetOrMakeRingArray(d, r);
+ if (!ringArray) return false;
+
+ if (b <= 0 || b > fVertexAxis.GetNbins()) {
+ AliWarning(Form("Vertex bin %3d out of range [1,%3d]",
+ b, fVertexAxis.GetNbins()));
+ return false;
+ }
+ h->SetName(Form("FMD%d%c_vtxbin%03d", d, r, b));
+ h->SetTitle(Form("Secondary map correction for FMD%d%c "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ d, r, b, fVertexAxis.GetBinLowEdge(b),
+ fVertexAxis.GetBinUpEdge(b)));
+ h->SetXTitle("#eta");
+ h->SetYTitle("#phi [radians]");
+ h->SetZTitle("#sum_{i} N_{ch,i,primary} / #sum_{i} N_{ch,i,FMD}");
+ h->SetDirectory(0);
+ h->SetStats(0);
+ ringArray->AddAtAndExpand(h, b-1);
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrSecondaryMap::SetCorrection(UShort_t d, Char_t r,
+ Double_t v, TH2D* h)
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0 || b > fVertexAxis.GetNbins()) {
+ AliWarning(Form("Vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return false;
+ }
+ return SetCorrection(d, r, UShort_t(b), h);
+}
+//____________________________________________________________________
+void
+AliFMDCorrSecondaryMap::Browse(TBrowser* b)
+{
+ b->Add(&fRingArray);
+ b->Add(&fVertexAxis);
+ b->Add(&fEtaAxis);
+}
+//____________________________________________________________________
+void
+AliFMDCorrSecondaryMap::Print(Option_t* option) const
+{
+ std::cout << "Secondary correction map" << std::endl;
+ fRingArray.Print(option);
+ fVertexAxis.Print(option);
+ fEtaAxis.Print(option);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+//
+// This class contains the secondary correction and the double hit
+// correction used in low-flux events.
+//
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRSECONDARYMAP_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRSECONDARYMAP_H
+#include <TObject.h>
+#include <TObjArray.h>
+#include <TAxis.h>
+class TH2D;
+
+/**
+ * This class contains the secondary correction.
+ *
+ * The secondary correction is given by
+ * @f[
+ * c_{r,v}(\eta,\varphi) =
+ * \frac{\sum_i N_{ch,i,v,\mbox{primary}}(\eta,\varphi)}{
+ * \sum_i N_{ch,i,r,v,\mbox{FMD}}(\eta,\varphi)}
+ * @f]
+ * where @f$N_{ch,i,v,\mbox{primary}}(\eta,\varphi)@f$ is the
+ * is the number of primary charged particles that fall within
+ * the @f$(\eta,\varphi)@f$ bin in event @f$i@f$ with vertex @f$v@f$,
+ * and is the total (primary <i>and</i> secondary) charged particles
+ * that hit ring @f$r@f$ within @f$(\eta,\varphi)@f$ bin in event
+ * @f$i@f$ with vertex @f$v@f$.
+ *
+ * These are generated from Monte-Carlo truth information.
+ *
+ * @ingroup pwg2_forward_corr
+ */
+class AliFMDCorrSecondaryMap : public TObject
+{
+public:
+ /**
+ * Default constructor
+ */
+ AliFMDCorrSecondaryMap();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrSecondaryMap(const AliFMDCorrSecondaryMap& o);
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliFMDCorrSecondaryMap();
+ /**
+ * @{
+ * @name Get corrections and parameters
+ */
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrSecondaryMap& operator=(const AliFMDCorrSecondaryMap& o);
+ /**
+ * Get the secondary correction @f$ c_{r,v}@f$
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ *
+ * @return The correction @f$ c_{r,v}@f$
+ */
+ TH2D* GetCorrection(UShort_t d, Char_t r, Double_t v) const;
+ /**
+ * Get the secondary correction @f$ c_{r,v}@f$
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ *
+ * @return The correction @f$ c_{r,v}@f$
+ */
+ TH2D* GetCorrection(UShort_t d, Char_t r, UShort_t b) const;
+ /**
+ * Get the vertex axis used
+ *
+ * @return vertex axis
+ */
+ const TAxis& GetVertexAxis() const { return fVertexAxis; }
+ /**
+ * Get the eta axis used
+ *
+ * @return eta axis
+ */
+ const TAxis& GetEtaAxis() const { return fEtaAxis; }
+ /* @} */
+
+ /**
+ * @{
+ * @name Set corrections and parameters
+ */
+ /**
+ * Set the secondary map correction @f$ c_{r,v}(\eta,\varphi)@f$.
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ * @param h @f$ c_{r,v}(\eta,\varphi)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(UShort_t d, Char_t r, Double_t v, TH2D* h);
+ /**
+ * Set the secondary map correction @f$ c_{r,v}(\eta,\varphi)@f$
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param d Detector number (1-3)
+ * @param r Ring identifier (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ * @param h @f$ c_{r,v}(\eta,\varphi)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(UShort_t d, Char_t r, UShort_t b, TH2D* h);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(const TAxis& axis);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(Int_t nBins, Double_t min, Double_t max);
+ /**
+ * Set the eta axis to use
+ *
+ * @param axis Eta axis
+ */
+ void SetEtaAxis(const TAxis& axis);
+ /**
+ * Set the eta axis to use
+ *
+ * @param axis Eta axis
+ */
+ void SetEtaAxis(Int_t nBins, Double_t min, Double_t max);
+ /* @} */
+
+ /**
+ * @{
+ * @name Auxiliary member functions
+ */
+ /**
+ * Declare this as a folder
+ *
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return true; }
+ /**
+ * Browse this object in the browser
+ *
+ * @param b
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Print this object
+ *
+ * @param option
+ */
+ void Print(Option_t* option="R") const; //*MENU*
+ /* @} */
+protected:
+ /**
+ * Find the vertex bin that corresponds to the passed vertex
+ *
+ * @param vertex The interaction points @f$z@f$-coordinate
+ *
+ * @return Vertex bin in @f$[1,N_{\box{vertex}}]@f$ or negative if
+ * out of range
+ */
+ Int_t FindVertexBin(Double_t vertex) const;
+ /**
+ * Get the index corresponding to the given ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Index (0 based) or negative in case of errors
+ */
+ Int_t GetRingIndex(UShort_t d, Char_t r) const;
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or null in case of problems
+ */
+ TObjArray* GetRingArray(UShort_t d, Char_t r) const;
+ /**
+ * Get the ring array corresponding to the specified ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Pointer to ring array, or newly created container
+ */
+ TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
+
+ TObjArray fRingArray; // Array of per-ring, per-vertex 2nd map
+ TAxis fVertexAxis; // The vertex axis
+ TAxis fEtaAxis; // The eta axis
+ ClassDef(AliFMDCorrSecondaryMap,1); //
+};
+
+//____________________________________________________________________
+inline void
+AliFMDCorrSecondaryMap::SetVertexAxis(Int_t nBins, Double_t min, Double_t max)
+{
+ fVertexAxis.Set(nBins, min, max);
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrSecondaryMap::SetVertexAxis(const TAxis& e)
+{
+ fVertexAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrSecondaryMap::SetEtaAxis(Int_t nBins, Double_t min, Double_t max)
+{
+ fEtaAxis.Set(nBins, min, max);
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrSecondaryMap::SetEtaAxis(const TAxis& e)
+{
+ fEtaAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
+}
+#endif
+// Local Variables:
+// mode: C++
+// End:
--- /dev/null
+#include "AliFMDCorrVertexBias.h"
+#include <TBrowser.h>
+#include <TH2D.h>
+#include <AliLog.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliFMDCorrVertexBias::AliFMDCorrVertexBias()
+ : fVertexArray(),
+ fVertexAxis(0,0,0)
+{
+ fVertexArray.SetOwner(kTRUE);
+ fVertexArray.SetName("rings");
+ fVertexAxis.SetName("vtxAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+
+}
+//____________________________________________________________________
+AliFMDCorrVertexBias::AliFMDCorrVertexBias(const AliFMDCorrVertexBias& o)
+ : TObject(o),
+ fVertexArray(o.fVertexArray),
+ fVertexAxis(o.fVertexAxis.GetNbins(), o.fVertexAxis.GetXmin(),
+ o.fVertexAxis.GetXmax())
+{
+ fVertexAxis.SetName("vtxAxis");
+ fVertexAxis.SetTitle("v_{z} [cm]");
+}
+//____________________________________________________________________
+AliFMDCorrVertexBias::~AliFMDCorrVertexBias()
+{
+ fVertexArray.Clear();
+}
+//____________________________________________________________________
+AliFMDCorrVertexBias&
+AliFMDCorrVertexBias::operator=(const AliFMDCorrVertexBias& o)
+{
+ fVertexArray = o.fVertexArray;
+ SetVertexAxis(o.fVertexAxis);
+
+ return *this;
+}
+//____________________________________________________________________
+TH2D*
+AliFMDCorrVertexBias::GetCorrection(Char_t r, Double_t v) const
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0) return 0;
+ return GetCorrection(r, UShort_t(b));
+}
+//____________________________________________________________________
+TH2D*
+AliFMDCorrVertexBias::GetCorrection(Char_t r, UShort_t b) const
+{
+ TObjArray* vertexarray = GetVertexArray(b);
+ if (!vertexarray) return 0;
+
+ Int_t idx = -1;
+ switch (r) {
+ case 'i': case 'I': idx = 0; break;
+ case 'o': case 'O': idx = 1; break;
+ }
+ if (idx < 0) {
+ AliWarning(Form("Unknown ting type %c, not one of [iIoO]", r));
+ return 0;
+ }
+
+ TObject* o = vertexarray->At(idx);
+ if (!o) {
+ AliWarning(Form("No vertex bias found for ring type %c in vertex bin %d",
+ r,b));
+ return 0;
+ }
+ return static_cast<TH2D*>(o);
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDCorrVertexBias::FindVertexBin(Double_t v) const
+{
+ if (fVertexAxis.GetNbins() <= 0) {
+ AliWarning("No vertex array defined");
+ return 0;
+ }
+ Int_t bin = const_cast<TAxis&>(fVertexAxis).FindBin(v);
+ if (bin <= 0 || bin > fVertexAxis.GetNbins()) {
+ AliWarning(Form("vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return 0;
+ }
+ return bin;
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrVertexBias::GetVertexArray(UShort_t v) const
+{
+ if (v <= 0 || v > fVertexAxis.GetNbins()) {
+ AliWarning(Form("vertex bin %d out of range [1,%d]",
+ v, fVertexAxis.GetNbins()));
+ return 0;
+ }
+
+ TObject* o = fVertexArray.At(v-1);
+ if (!o) {
+ AliWarning(Form("No array found for vertex bin %d", v));
+ return 0;
+ }
+
+ return static_cast<TObjArray*>(o);
+}
+//____________________________________________________________________
+TObjArray*
+AliFMDCorrVertexBias::GetOrMakeVertexArray(UShort_t v)
+{
+ if (v <= 0 || v > fVertexAxis.GetNbins()) {
+ AliWarning(Form("vertex bin %d out of range [1,%d]",
+ v, fVertexAxis.GetNbins()));
+ return 0;
+ }
+
+ TObject* o = fVertexArray.At(v-1);
+ if (!o) {
+ TObjArray* a = new TObjArray(fVertexAxis.GetNbins());
+ a->SetName(Form("vertexbin%02d", v));
+ a->SetOwner(kTRUE);
+ fVertexArray.AddAtAndExpand(a, v-1);
+ return a;
+ }
+
+ return static_cast<TObjArray*>(fVertexArray.At(v-1));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrVertexBias::SetCorrection(Char_t r, UShort_t b, TH2D* h)
+{
+ TObjArray* vertexarray = GetOrMakeVertexArray(b);
+ if (!vertexarray) return false;
+
+
+ Int_t idx = -1;
+ switch (r) {
+ case 'i': case 'I': idx = 0; break;
+ case 'o': case 'O': idx = 1; break;
+ }
+ if (idx < 0) {
+ AliWarning(Form("Unknown ting type %c, not one of [iIoO]", r));
+ return false;
+ }
+ h->SetName(Form("FMDX%c", r));
+ h->SetTitle(Form("Vertex bias correction for %c rings "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ r, b, fVertexAxis.GetBinLowEdge(b),
+ fVertexAxis.GetBinUpEdge(b)));
+ h->SetXTitle("#eta");
+ h->SetYTitle("#phi [radians]");
+ h->SetZTitle("1/N_{t}#sum_{i}^{N_{tv}} N_{ch,i,primary} / "
+ "1/N_{v}#sum_{i}^{N_{v}} N_{ch,i,primary}");
+ h->SetDirectory(0);
+ h->SetStats(0);
+
+ vertexarray->AddAtAndExpand(h, idx);
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDCorrVertexBias::SetCorrection(Char_t r, Double_t v, TH2D* h)
+{
+ Int_t b = FindVertexBin(v);
+ if (b <= 0 || b > fVertexAxis.GetNbins()) {
+ AliWarning(Form("Vertex %+8.4f out of range [%+8.4f,%+8.4f]",
+ v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
+ return false;
+ }
+ return SetCorrection(r, UShort_t(b), h);
+}
+//____________________________________________________________________
+void
+AliFMDCorrVertexBias::Browse(TBrowser* b)
+{
+ b->Add(&fVertexArray);
+ b->Add(&fVertexAxis);
+}
+//____________________________________________________________________
+void
+AliFMDCorrVertexBias::Print(Option_t* option) const
+{
+ std::cout << "Vertex bias correction" << std::endl;
+ fVertexArray.Print(option);
+ fVertexAxis.Print(option);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+//
+// This class contains the secondary correction and the double hit
+// correction used in low-flux events.
+//
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRVERTEXBIAS_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS_ALIFMDCORRVERTEXBIAS_H
+#include <TObject.h>
+#include <TObjArray.h>
+#include <TAxis.h>
+class TH2D;
+
+/**
+ * This class contains the correction for the bias introduced by
+ * different vertex bins
+ *
+ * The correction is given by
+ * @f[
+ * b_{v}(\eta,\varphi) = \frac{1/N_{t}\sum_i^{N_{tv}} N_{ch,i,primary}}{
+ * 1/N_{v}\sum_i^{N_{v}} N_{ch,i,primary}}
+ * @f]
+ *
+ * where @f$N_{ch,i,primary}@f$ is the number of primary particles in
+ * the given @f$(\eta,\varphi)@f$, and where the denominator sum runs
+ * over all events with a vertex within the given vertex bin, and the
+ * sum in the numerator runs over only events that have a valid
+ * trigger and reconstructed vertex. @f$ N_{t}@f$ is the number of
+ * events with a valid trigger (but not necessarily a valid vertex).
+ * The vertex information used @f$v@f$ is in all cases the MC truth
+ * vertex
+ *
+ * These are generated from Monte-Carlo truth and ESD information.
+ *
+ * @ingroup pwg2_forward_corr
+ */
+class AliFMDCorrVertexBias : public TObject
+{
+public:
+ /**
+ * Default constructor
+ */
+ AliFMDCorrVertexBias();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrVertexBias(const AliFMDCorrVertexBias& o);
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliFMDCorrVertexBias();
+ /**
+ * @{
+ * @name Get corrections and parameters
+ */
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrVertexBias& operator=(const AliFMDCorrVertexBias& o);
+ /**
+ * Get the vertex bias correction @f$ b_{v}@f$
+ *
+ * @param r Ring type (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ *
+ * @return The correction @f$ b_{v}@f$
+ */
+ TH2D* GetCorrection(Char_t r, Double_t v) const;
+ /**
+ * Get the vertex bias correction @f$ b_{v}@f$
+ *
+ * @param r Ring type (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ *
+ * @return The correction @f$ b_{v}@f$
+ */
+ TH2D* GetCorrection(Char_t r, UShort_t b) const;
+ /**
+ * Get the vertex axis used
+ *
+ * @return vertex axis
+ */
+ const TAxis& GetVertexAxis() const { return fVertexAxis; }
+ /* @} */
+
+ /**
+ * @{
+ * @name Set corrections and parameters
+ */
+ /**
+ * Set the vertex bias correction @f$ b_{v}(\eta,\varphi)@f$.
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param r Ring type (I or O)
+ * @param v Primary interaction point @f$z@f$ coordinate
+ * @param h @f$ b_{v}(\eta,\varphi)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(Char_t r, Double_t v, TH2D* h);
+ /**
+ * Set the vertex bias correction @f$ b_{v}(\eta,\varphi)@f$
+ * Note, that the object takes ownership of the passed pointer.
+ *
+ * @param r Ring type (I or O)
+ * @param b Bin corresponding to the primary interaction point
+ * @f$z@f$ coordinate (1 based)
+ * @param h @f$ b_{v}(\eta,\varphi)@f$
+ *
+ * @return true if operation succeeded
+ */
+ Bool_t SetCorrection(Char_t r, UShort_t b, TH2D* h);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(const TAxis& axis);
+ /**
+ * Set the vertex axis to use
+ *
+ * @param axis Vertex axis
+ */
+ void SetVertexAxis(Int_t nBins, Double_t min, Double_t max);
+ /* @} */
+
+ /**
+ * @{
+ * @name Auxiliary member functions
+ */
+ /**
+ * Declare this as a folder
+ *
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return true; }
+ /**
+ * Browse this object in the browser
+ *
+ * @param b
+ */
+ void Browse(TBrowser* b);
+ /**
+ * Print this object
+ *
+ * @param option
+ */
+ void Print(Option_t* option="R") const; //*MENU*
+ /* @} */
+protected:
+ /**
+ * Find the vertex bin that corresponds to the passed vertex
+ *
+ * @param vertex The interaction points @f$z@f$-coordinate
+ *
+ * @return Vertex bin in @f$[1,N_{\box{vertex}}]@f$ or negative if
+ * out of range
+ */
+ Int_t FindVertexBin(Double_t vertex) const;
+ /**
+ * Get the vertex array corresponding to the specified ring
+ *
+ * @param v vertex bin (1 based)
+ *
+ * @return Pointer to vertex array, or null in case of problems
+ */
+ TObjArray* GetVertexArray(UShort_t v) const;
+ /**
+ * Get the vertex array corresponding to the specified ring
+ *
+ * @param v vertex bin (1 based)
+ *
+ * @return Pointer to vertex array, or newly created container
+ */
+ TObjArray* GetOrMakeVertexArray(UShort_t v);
+
+ TObjArray fVertexArray; // Array of per-ring, per-vertex 2nd map
+ TAxis fVertexAxis; // The vertex axis
+ ClassDef(AliFMDCorrVertexBias,1); //
+};
+
+//____________________________________________________________________
+inline void
+AliFMDCorrVertexBias::SetVertexAxis(Int_t nBins, Double_t min, Double_t max)
+{
+ fVertexAxis.Set(nBins, min, max);
+}
+//____________________________________________________________________
+inline void
+AliFMDCorrVertexBias::SetVertexAxis(const TAxis& e)
+{
+ fVertexAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
+}
+#endif
+// Local Variables:
+// mode: C++
+// End:
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
-#include "AliFMDAnaParameters.h"
+#include "AliForwardCorrectionManager.h"
+// #include "AliFMDAnaParameters.h"
#include "AliLog.h"
#include <TH2D.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
ClassImp(AliFMDCorrections)
#if 0
AliFMDCorrections::AliFMDCorrections()
: TNamed(),
fRingHistos(),
- fMultCut(0.3),
fDebug(0)
{}
AliFMDCorrections::AliFMDCorrections(const char* title)
: TNamed("fmdCorrections", title),
fRingHistos(),
- fMultCut(0.3),
fDebug(0)
{
fRingHistos.SetName(GetName());
AliFMDCorrections::AliFMDCorrections(const AliFMDCorrections& o)
: TNamed(o),
fRingHistos(),
- fMultCut(o.fMultCut),
fDebug(o.fDebug)
{
TIter next(&o.fRingHistos);
{
TNamed::operator=(o);
- fMultCut = o.fMultCut;
fDebug = o.fDebug;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
//____________________________________________________________________
Bool_t
AliFMDCorrections::Correct(AliForwardUtil::Histos& hists,
- Int_t vtxbin)
+ UShort_t vtxbin)
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ UShort_t uvb = vtxbin;
for (UShort_t d=1; d<=3; d++) {
UShort_t nr = (d == 1 ? 1 : 2);
for (UShort_t q=0; q<nr; q++) {
Char_t r = (q == 0 ? 'I' : 'O');
TH2D* h = hists.Get(d,r);
RingHistos* rh= GetRingHistos(d,r);
- TH2F* bg= pars->GetBackgroundCorrection(d, r, vtxbin);
- TH2F* ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
+ //TH2F* bg= pars->GetBackgroundCorrection(d, r, vtxbin);
+ //TH2F* ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
+ TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
+ TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
if (!bg) {
AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
- d, r, vtxbin));
+ d, r, uvb));
continue;
}
if (!ef) {
- AliWarning(Form("No event %s selection efficiency in vertex bin %d",
- "INEL", vtxbin));
+ AliWarning(Form("No event vertex bias correction in vertex bin %d",
+ uvb));
continue;
}
// Divide by the event selection efficiency
h->Divide(ef);
- if(!pars->SharingEffPresent()) {
- AliWarning("No sharing efficiencies");
+
+ // if(!pars->SharingEffPresent()) {
+ // AliWarning("No sharing efficiencies");
+ // continue;
+ // }
+ // TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin);
+ if (!fcm.GetMergingEfficiency()) {
+ AliWarning("No merging efficiencies");
+ continue;
+ }
+ TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
+ if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
+ AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
+ d, r, uvb));
continue;
}
- TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin);
for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
Float_t c = sf->GetBinContent(ieta);
}
}
+//____________________________________________________________________
+void
+AliFMDCorrections::Print(Option_t* /* option */) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ std::cout << ind << "AliFMDCorrections: " << GetName() << std::endl;
+}
+
//====================================================================
AliFMDCorrections::RingHistos::RingHistos()
: AliForwardUtil::RingHistos(),
*
* @return true on successs
*/
- virtual Bool_t Correct(AliForwardUtil::Histos& hists, Int_t vtxBin);
+ virtual Bool_t Correct(AliForwardUtil::Histos& hists, UShort_t vtxBin);
/**
* Scale the histograms to the total number of events
*
* @param dbg Debug level
*/
void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Internal data structure to keep track of the histograms
RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
TList fRingHistos; // List of histogram containers
- Double_t fMultCut; // Low cut on scaled energy loss
Int_t fDebug; // Debug level
ClassDef(AliFMDCorrections,1); // Calculate Nch density
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
-#include "AliFMDAnaParameters.h"
+#include "AliForwardCorrectionManager.h"
+// #include "AliFMDAnaParameters.h"
#include "AliLog.h"
#include <TH2D.h>
+#include <TProfile.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
ClassImp(AliFMDDensityCalculator)
#if 0
AliFMDDensityCalculator::AliFMDDensityCalculator()
: TNamed(),
fRingHistos(),
- fMultCut(0.3),
+ fMultCut(0),
fSumOfWeights(0),
fWeightedSum(0),
fCorrections(0),
+ fMaxParticles(5),
+ fAccI(0),
+ fAccO(0),
fDebug(0)
{}
AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
: TNamed("fmdDensityCalculator", title),
fRingHistos(),
- fMultCut(0.3),
+ fMultCut(0),
fSumOfWeights(0),
fWeightedSum(0),
fCorrections(0),
+ fMaxParticles(5),
+ fAccI(0),
+ fAccO(0),
fDebug(0)
{
fRingHistos.SetName(GetName());
fRingHistos.Add(new RingHistos(3, 'O'));
fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
200, 0, 20);
+ fSumOfWeights->SetFillColor(kRed+1);
+ fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
200, 0, 20);
+ fWeightedSum->SetFillColor(kBlue+1);
+ fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
fCorrections = new TH1D("corrections", "Distribution of corrections",
100, 0, 10);
+ fCorrections->SetFillColor(kBlue+1);
+ fCorrections->SetXTitle("correction");
+ fAccI = GenerateAcceptanceCorrection('I');
+ fAccO = GenerateAcceptanceCorrection('O');
}
//____________________________________________________________________
fSumOfWeights(o.fSumOfWeights),
fWeightedSum(o.fWeightedSum),
fCorrections(o.fCorrections),
+ fMaxParticles(o.fMaxParticles),
+ fAccI(o.fAccI),
+ fAccO(o.fAccO),
fDebug(o.fDebug)
{
TIter next(&o.fRingHistos);
{
TNamed::operator=(o);
- fMultCut = o.fMultCut;
- fDebug = o.fDebug;
+ fMultCut = o.fMultCut;
+ fDebug = o.fDebug;
+ fMaxParticles = o.fMaxParticles;
+ fAccI = o.fAccI;
+ fAccO = o.fAccO;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
+//____________________________________________________________________
+Double_t
+AliFMDDensityCalculator::GetMultCut() const
+{
+ if (fMultCut > 0) return fMultCut;
+
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ AliFMDCorrELossFit* fits = fcm.GetELossFit();
+ return fits->GetLowCut();
+}
+
//____________________________________________________________________
Bool_t
AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- Int_t vtxbin,
+ UShort_t vtxbin,
Bool_t lowFlux)
{
for (UShort_t d=1; d<=3; d++) {
Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
rh->fEvsN->Fill(mult,n);
+ rh->fEtaVsN->Fill(eta, n);
Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
fCorrections->Fill(c);
if (c > 0) n /= c;
rh->fEvsM->Fill(mult,n);
+ rh->fEtaVsM->Fill(eta, n);
+ rh->fCorr->Fill(eta, c);
h->Fill(eta,phi,n);
rh->fDensity->Fill(eta,phi,n);
Char_t r,
UShort_t /*s*/,
UShort_t /*t*/,
- Int_t /*v*/,
+ UShort_t /*v*/,
Float_t eta,
Bool_t lowFlux) const
{
- if (mult <= fMultCut) return 0;
+ if (mult <= GetMultCut()) return 0;
if (lowFlux) return 1;
+
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
+ if (!fit) {
+ AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
+ return 0;
+ }
+
+ Int_t m = fit->FindMaxWeight();
+ if (m < 1) {
+ AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
+ return 0;
+ }
+ UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
+ Double_t ret = fit->EvaluateWeighted(mult, n);
+
+ if (fDebug > 10) {
+ AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
+ }
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ fWeightedSum->Fill(ret);
+ fSumOfWeights->Fill(ret);
+ return ret;
+#if 0
Float_t mpv = pars->GetMPV(d,r,eta);
Float_t w = pars->GetSigma(d,r,eta);
Float_t w2 = pars->Get2MIPWeight(d,r,eta);
fWeightedSum->Fill(wsum);
fSumOfWeights->Fill(sum);
return (sum > 0) ? wsum / sum : 1;
+#endif
}
//_____________________________________________________________________
Float_t
-AliFMDDensityCalculator::Correction(UShort_t d, Char_t r, UShort_t /*s*/,
- UShort_t t, Int_t v, Float_t eta,
- Bool_t lowFlux) const
+AliFMDDensityCalculator::Correction(UShort_t d,
+ Char_t r,
+ UShort_t /*s*/,
+ UShort_t t,
+ UShort_t /*v*/,
+ Float_t eta,
+ Bool_t lowFlux) const
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
Float_t correction = AcceptanceCorrection(r,t);
if (lowFlux) {
- TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
+ // TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
+ TH1D* dblHitCor = 0;
+ if (fcm.GetDoubleHit())
+ dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
+
if (dblHitCor) {
- Float_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
- if (dblC > 0)
- correction *= dblC;
+ Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
+ if (dblC > 0) correction *= dblC;
}
- else
+ else {
AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
+ }
}
+#if 0
TH1F* deadCor = pars->GetFMDDeadCorrection(v);
if (deadCor) {
Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
if (deadC > 0)
correction *= deadC;
}
-
+#endif
+
return correction;
}
+//_____________________________________________________________________
+TH1D*
+AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
+{
+ const Double_t ic1[] = { 4.9895, 15.3560 };
+ const Double_t ic2[] = { 1.8007, 17.2000 };
+ const Double_t oc1[] = { 4.2231, 26.6638 };
+ const Double_t oc2[] = { 1.8357, 27.9500 };
+ const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
+ const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
+ Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
+ Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
+ Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
+ Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
+ Float_t basearc = 2 * TMath::Pi() / nSec;
+ Double_t rad = maxR - minR;
+ Float_t segment = rad / nStrips;
+ Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
+
+ // Numbers used to find end-point of strip.
+ // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
+ Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
+ Float_t dx = c2[0] - c1[0];
+ Float_t dy = c2[1] - c1[1];
+ Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
+
+ TH1D* ret = new TH1D(Form("acc%c", r),
+ Form("Acceptance correction for FMDx%c", r),
+ nStrips, -.5, nStrips-.5);
+ ret->SetXTitle("Strip");
+ ret->SetYTitle("#varphi acceptance");
+ ret->SetDirectory(0);
+ ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
+ ret->SetFillStyle(3001);
+
+ for (Int_t t = 0; t < nStrips; t++) {
+ Float_t radius = minR + t * segment;
+
+ // If the radius of the strip is smaller than the radius corresponding
+ // to the first corner we have a full strip length
+ if (radius <= cr) {
+ ret->SetBinContent(t+1, 1);
+ continue;
+ }
+
+ // Next, we should find the end-point of the strip - that is,
+ // the coordinates where the line from c1 to c2 intersects a circle
+ // with radius given by the strip.
+ // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
+ // Calculate the determinant
+ Float_t det = radius * radius * dr * dr - D*D;
+
+ if (det <= 0) {
+ // <0 means No intersection
+ // =0 means Exactly tangent
+ ret->SetBinContent(t+1, 1);
+ continue;
+ }
+
+ // Calculate end-point and the corresponding opening angle
+ Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
+ Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
+ Float_t th = TMath::ATan2(x, y);
+
+ ret->SetBinContent(t+1, th / basearc);
+ }
+ return ret;
+}
//_____________________________________________________________________
Float_t
AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
-
- //AliFMDRing fmdring(ring);
- //fmdring.Init();
- Float_t rad = pars->GetMaxR(r)-pars->GetMinR(r);
- Float_t slen = pars->GetStripLength(r,t);
- Float_t sblen = pars->GetBaseStripLength(r,t);
- Float_t nstrips = (r == 'I' ? 512 : 256);
- Float_t segment = rad / nstrips;
- Float_t radius = pars->GetMinR(r) + segment*t;
+ TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
+ return acc->GetBinContent(t+1);
+
+#if 0
+ const Double_t ic1[] = { 4.9895, 15.3560 };
+ const Double_t ic2[] = { 1.8007, 17.2000 };
+ const Double_t oc1[] = { 4.2231, 26.6638 };
+ const Double_t oc2[] = { 1.8357, 27.9500 };
+ const Double_t* c1 = (r == 'I' ? ic1 : oc1);
+ const Double_t* c2 = (r == 'I' ? ic2 : oc2);
+ Double_t minR = (r == 'I' ? 4.5213 : 15.4);
+ Double_t maxR = (r == 'I' ? 17.2 : 28.0);
+ Int_t nStrips = (r == 'I' ? 512 : 256);
+ Int_t nSec = (r == 'I' ? 20 : 40);
+ Float_t basearc = 2 * TMath::Pi() / nSec;
+ Double_t rad = maxR - minR;
+ Float_t segment = rad / nStrips;
+ Float_t radius = minR + t * segment;
+
+ // Old method - calculate full strip area and take ratio to extended
+ // strip area
+ Float_t baselen = basearc * radius;
+ Float_t slope = (c1[1] - c2[1]) / (c1[0] - c2[0]);
+ Float_t constant = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
+ Float_t d = (TMath::Power(TMath::Abs(radius*slope),2) +
+ TMath::Power(radius,2) - TMath::Power(constant,2));
+
+ // If below corners return 1
+ if (d >= 0) return 1;
+
+ Float_t x = ((-TMath::Sqrt(d) - slope * constant) /
+ (1+TMath::Power(slope, 2)));
+ Float_t y = slope*x + constant;
+
+ // If x is larger than corner x or y less than corner y, we have full
+ // length strip
+ if(x >= c1[0] || y <= c1[1]) return 1;
+
+ //One sector since theta is by definition half-hybrid
+ Float_t theta = TMath::ATan2(x,y);
+ Float_t arclen = radius * theta;
- Float_t basearea1 = 0.5*sblen*TMath::Power(radius,2);
- Float_t basearea2 = 0.5*sblen*TMath::Power((radius-segment),2);
+ // Calculate the area of a strip with no cut
+ Float_t basearea1 = 0.5 * baselen * TMath::Power(radius,2);
+ Float_t basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
Float_t basearea = basearea1 - basearea2;
-
- Float_t area1 = 0.5*slen*TMath::Power(radius,2);
- Float_t area2 = 0.5*slen*TMath::Power((radius-segment),2);
+
+ // Calculate the area of a strip with cut
+ Float_t area1 = 0.5 * arclen * TMath::Power(radius,2);
+ Float_t area2 = 0.5 * arclen * TMath::Power((radius-segment),2);
Float_t area = area1 - area2;
+ // Acceptance is ratio
return area/basearea;
+#endif
}
//____________________________________________________________________
d->Add(fWeightedSum);
d->Add(fSumOfWeights);
d->Add(fCorrections);
+ d->Add(fAccI);
+ d->Add(fAccO);
+
TIter next(&fRingHistos);
RingHistos* o = 0;
while ((o = static_cast<RingHistos*>(next()))) {
o->Output(d);
}
}
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::Print(Option_t*) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ std::cout << ind << "AliFMDDensityCalculator: " << GetName() << '\n'
+ << ind << " Multiplicity cut: " << fMultCut << '\n'
+ << ind << " Max(particles): " << fMaxParticles
+ << std::endl;
+}
//====================================================================
AliFMDDensityCalculator::RingHistos::RingHistos()
: AliForwardUtil::RingHistos(),
fEvsN(0),
fEvsM(0),
+ fEtaVsN(0),
+ fEtaVsM(0),
+ fCorr(0),
fDensity(0)
{}
: AliForwardUtil::RingHistos(d,r),
fEvsN(0),
fEvsM(0),
+ fEtaVsN(0),
+ fEtaVsM(0),
+ fCorr(0),
fDensity(0)
{
fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()),
- Form("Energy loss vs uncorrected inclusive "
+ Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
"N_{ch} in %s", fName.Data()),
- 100, -.5, 24.5, 100, -.5, 24.5);
+ 2500, -.5, 24.5, 2500, -.5, 24.5);
fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()),
- Form("Energy loss vs corrected inclusive N_{ch} in %s",
- fName.Data()), 100, -.5, 24.5, 100, -.5, 24.5);
+ Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
+ "N_{ch} in %s", fName.Data()),
+ 2500, -.5, 24.5, 2500, -.5, 24.5);
fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
fEvsN->Sumw2();
fEvsN->SetDirectory(0);
- fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
- fEvsN->SetYTitle("Inclusive N_{ch} (corrected)");
+ fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
+ fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
fEvsM->Sumw2();
fEvsM->SetDirectory(0);
+ fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
+ Form("Average inclusive N_{ch} vs #eta (uncorrected) "
+ "in %s", fName.Data()), 200, -4, 6);
+ fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
+ Form("Average inclusive N_{ch} vs #eta (corrected) "
+ "in %s", fName.Data()), 200, -4, 6);
+ fEtaVsN->SetXTitle("#eta");
+ fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
+ fEtaVsN->SetDirectory(0);
+ fEtaVsN->SetLineColor(Color());
+ fEtaVsN->SetFillColor(Color());
+ fEtaVsM->SetXTitle("#eta");
+ fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
+ fEtaVsM->SetDirectory(0);
+ fEtaVsM->SetLineColor(Color());
+ fEtaVsM->SetFillColor(Color());
+
+
+ fCorr = new TProfile(Form("%s_corr", fName.Data()),
+ Form("Average correction in %s", fName.Data()),
+ 200, -4, 6);
+ fCorr->SetXTitle("#eta");
+ fCorr->SetYTitle("#LT correction#GT");
+ fCorr->SetDirectory(0);
+ fCorr->SetLineColor(Color());
+ fCorr->SetFillColor(Color());
+
fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()),
- Form("Inclusive N_{ch} densitu in %s", fName.Data()),
+ Form("Inclusive N_{ch} density in %s", fName.Data()),
200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
0, 2*TMath::Pi());
fDensity->SetDirectory(0);
: AliForwardUtil::RingHistos(o),
fEvsN(o.fEvsN),
fEvsM(o.fEvsM),
+ fEtaVsN(o.fEtaVsN),
+ fEtaVsM(o.fEtaVsM),
+ fCorr(o.fCorr),
fDensity(o.fDensity)
{}
if (fEvsN) delete fEvsN;
if (fEvsM) delete fEvsM;
+ if (fEtaVsN) delete fEtaVsN;
+ if (fEtaVsM) delete fEtaVsM;
+ if (fCorr) delete fCorr;
if (fDensity) delete fDensity;
fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
+ fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
+ fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
+ fCorr = static_cast<TProfile*>(o.fCorr->Clone());
fDensity = static_cast<TH2D*>(o.fDensity->Clone());
return *this;
{
if (fEvsN) delete fEvsN;
if (fEvsM) delete fEvsM;
+ if (fEtaVsN) delete fEtaVsN;
+ if (fEtaVsM) delete fEtaVsM;
+ if (fCorr) delete fCorr;
if (fDensity) delete fDensity;
}
TList* d = DefineOutputList(dir);
d->Add(fEvsN);
d->Add(fEvsM);
+ d->Add(fEtaVsN);
+ d->Add(fEtaVsM);
+ d->Add(fCorr);
d->Add(fDensity);
}
class AliESDFMD;
class TH2D;
class TH1D;
+class TProfile;
/**
* This class calculates the inclusive charged particle density
*/
virtual Bool_t Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- Int_t vtxBin, Bool_t lowFlux);
+ UShort_t vtxBin, Bool_t lowFlux);
/**
* Scale the histograms to the total number of events
*
* @param dbg Debug level
*/
void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /**
+ * Maximum particle weight to use
+ *
+ * @param m
+ */
+ void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
+ /**
+ * Set the lower multiplicity cut. This overrides the setting in
+ * the energy loss fits.
+ *
+ * @param cut Cut to use
+ */
+ void SetMultCut(Double_t cut) { fMultCut = cut; }
+ /**
+ * Get the multiplicity cut. If the user has set fMultCut (via
+ * SetMultCut) then that value is used. If not, then the lower
+ * value of the fit range for the enery loss fits is returned.
+ *
+ * @return Lower cut on multiplicity
+ */
+ Double_t GetMultCut() const;
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Get the number of particles corresponding to the signal mult
*/
virtual Float_t NParticles(Float_t mult,
UShort_t d, Char_t r, UShort_t s, UShort_t t,
- Int_t v, Float_t eta, Bool_t lowFlux) const;
+ UShort_t v, Float_t eta, Bool_t lowFlux) const;
/**
* Get the inverse correction factor. This consist of
*
* @return
*/
virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
- Int_t v, Float_t eta, Bool_t lowFlux) const;
+ UShort_t v, Float_t eta, Bool_t lowFlux) const;
/**
* Get the acceptance correction for strip @a t in an ring of type @a r
*
* @return Inverse acceptance correction
*/
virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
-
+ /**
+ * Generate the acceptance corrections
+ *
+ * @param r Ring to generate for
+ *
+ * @return Newly allocated histogram of acceptance corrections
+ */
+ virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const;
/**
* Internal data structure to keep track of the histograms
*/
void ScaleHistograms(TList* dir, Int_t nEvents);
TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
+ TProfile* fEtaVsN; // Average uncorrected Nch vs eta
+ TProfile* fEtaVsM; // Average corrected Nch vs eta
+ TProfile* fCorr; // Average correction vs eta
TH2D* fDensity; // Distribution inclusive Nch
ClassDef(RingHistos,1);
};
TH1D* fSumOfWeights; //! Histogram
TH1D* fWeightedSum; //! Histogram
TH1D* fCorrections; //! Histogram
+ UShort_t fMaxParticles; // Maximum particle weight to use
+ TH1D* fAccI; // Acceptance correction for inner rings
+ TH1D* fAccO; // Acceptance correction for outer rings
Int_t fDebug; // Debug level
ClassDef(AliFMDDensityCalculator,1); // Calculate Nch density
// Class to do the energy correction of FMD ESD data
//
#include "AliFMDEnergyFitter.h"
+#include "AliForwardCorrectionManager.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TH1.h>
#include <TF1.h>
#include <TMath.h>
-#include "AliFMDAnaParameters.h"
#include <AliLog.h>
#include <TClonesArray.h>
#include <TFitResult.h>
#include <THStack.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
ClassImp(AliFMDEnergyFitter)
#if 0
; // This is for Emacs
#endif
+namespace {
+ const char* fgkEDistFormat = "%s_etabin%03d";
+}
//____________________________________________________________________
fMinEntries(100),
fFitRangeBinWidth(4),
fDoFits(false),
+ fDoMakeObject(false),
fEtaAxis(),
fMaxE(10),
fNEbins(300),
fUseIncreasingBins(true),
fMaxRelParError(.25),
fMaxChi2PerNDF(20),
+ fMinWeight(1e-7),
fDebug(0)
{}
fMinEntries(100),
fFitRangeBinWidth(4),
fDoFits(false),
+ fDoMakeObject(false),
fEtaAxis(0,0,0),
fMaxE(10),
fNEbins(300),
fUseIncreasingBins(true),
fMaxRelParError(.25),
- fMaxChi2PerNDF(20),
+ fMaxChi2PerNDF(20),
+ fMinWeight(1e-7),
fDebug(3)
{
fEtaAxis.SetName("etaAxis");
fMinEntries(o.fMinEntries),
fFitRangeBinWidth(o.fFitRangeBinWidth),
fDoFits(o.fDoFits),
+ fDoMakeObject(o.fDoMakeObject),
fEtaAxis(o.fEtaAxis),
fMaxE(o.fMaxE),
fNEbins(o.fNEbins),
fUseIncreasingBins(o.fUseIncreasingBins),
fMaxRelParError(o.fMaxRelParError),
fMaxChi2PerNDF(o.fMaxChi2PerNDF),
+ fMinWeight(o.fMinWeight),
fDebug(o.fDebug)
{
TIter next(&o.fRingHistos);
fMinEntries = o.fMinEntries;
fFitRangeBinWidth= o.fFitRangeBinWidth;
fDoFits = o.fDoFits;
+ fDoMakeObject = o.fDoMakeObject;
fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
fDebug = o.fDebug;
+ fMaxRelParError= o.fMaxRelParError;
+ fMaxChi2PerNDF = o.fMaxChi2PerNDF;
+ fMinWeight = o.fMinWeight;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
}
}
+
+ if (!fDoMakeObject) return;
+
+ MakeCorrectionsObject(d);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
+{
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+
+ AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
+ obj->SetEtaAxis(fEtaAxis);
+ obj->SetLowCut(fLowCut);
+
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError,
+ fMaxChi2PerNDF, fMinWeight);
+ }
+
+ TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
+ TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
+ AliInfo(Form("Object %s created in output - should be extracted and copied "
+ "to %s", oName.Data(), fileName.Data()));
+ d->Add(obj, oName.Data());
}
//____________________________________________________________________
while ((o = static_cast<RingHistos*>(next())))
o->fDebug = dbg;
}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::Print(Option_t*) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
+ << ind << " Low cut: " << fLowCut << " E/E_mip\n"
+ << ind << " Max(particles): " << fNParticles << '\n'
+ << ind << " Min(entries): " << fMinEntries << '\n'
+ << ind << " Fit range: "
+ << fFitRangeBinWidth << " bins\n"
+ << ind << " Make fits: "
+ << (fDoFits ? "yes\n" : "no\n")
+ << ind << " Make object: "
+ << (fDoMakeObject ? "yes\n" : "no\n")
+ << ind << " Max E/E_mip: " << fMaxE << '\n'
+ << ind << " N bins: " << fNEbins << '\n'
+ << ind << " Increasing bins: "
+ << (fUseIncreasingBins ?"yes\n" : "no\n")
+ << ind << " max(delta p/p): " << fMaxRelParError << '\n'
+ << ind << " max(chi^2/nu): " << fMaxChi2PerNDF << '\n'
+ << ind << " min(a_i): " << fMinWeight << std::endl;
+}
//====================================================================
AliFMDEnergyFitter::RingHistos::RingHistos()
fEmpty(0),
fEtaEDists(),
fList(0),
+ fFits("AliFMDCorrELossFit::ELossFit"),
fDebug(0)
{}
fEmpty(0),
fEtaEDists(),
fList(0),
+ fFits("AliFMDCorrELossFit::ELossFit"),
fDebug(0)
{
fEtaEDists.SetName("EDists");
fEmpty(o.fEmpty),
fEtaEDists(),
fList(0),
+ fFits("AliFMDCorrELossFit::ELossFit"),
fDebug(0)
{
+ fFits.Clear();
TIter next(&o.fEtaEDists);
TObject* obj = 0;
while ((obj = next())) fEtaEDists.Add(obj->Clone());
if (fEmpty) delete fEmpty;
fEtaEDists.Delete();
if (fList) fList->Delete();
+ fFits.Clear();
fEDist = static_cast<TH1D*>(o.fEDist->Clone());
fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
Bool_t incr)
{
TH1D* h = 0;
- TString name = Form("%s_etabin%03d", fName.Data(), ieta);
+ TString name = Form(fgkEDistFormat, fName.Data(), ieta);
TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
"(#eta bin %d)", fName.Data(), emin, emax, ieta);
if (!incr)
}
//____________________________________________________________________
TObjArray*
-AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
- Double_t lowCut,
- UShort_t nParticles,
- UShort_t minEntries,
- UShort_t minusBins,
- Double_t relErrorCut,
- Double_t chi2nuCut) const
+AliFMDEnergyFitter::RingHistos::Fit(TList* dir,
+ const TAxis& eta,
+ Double_t lowCut,
+ UShort_t nParticles,
+ UShort_t minEntries,
+ UShort_t minusBins,
+ Double_t relErrorCut,
+ Double_t chi2nuCut) const
{
// Get our ring list from the output container
TList* l = GetOutputList(dir);
// Normalize peak to 1
Double_t max = dist->GetMaximum();
+ if (max <= 0) continue;
dist->Scale(1/max);
// Check that we have enough entries
- if (dist->GetEntries() <= minEntries) continue;
+ if (dist->GetEntries() <= minEntries) {
+ AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
+ i, dist->GetName(), Int_t(dist->GetEntries()),
+ minEntries));
+ continue;
+ }
// Now fit
TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
relErrorCut,chi2nuCut);
if (!res) continue;
// dist->GetListOfFunctions()->Add(res);
-
+
// Store eta limits
low = TMath::Min(low,i+1);
high = TMath::Max(high,i+1);
// Fit the full-ring histogram
TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
if (total && total->GetEntries() >= minEntries) {
+
+ // Scale to the bin-width
+ total->Scale(1., "width");
+
+ // Normalize peak to 1
+ Double_t max = total->GetMaximum();
+ if (max > 0) total->Scale(1/max);
+
TF1* res = FitHist(total,lowCut,nParticles,minusBins,
relErrorCut,chi2nuCut);
if (res) {
TObjLink* lnk = dists->FirstLink();
while (lnk) {
TH1* h = static_cast<TH1*>(lnk->GetObject());
- if (h->GetEntries() <= 0 ||
- h->GetListOfFunctions()->GetEntries() <= 0) {
+ bool remove = false;
+ if (h->GetEntries() <= 0) {
+ // AliWarning(Form("No entries in %s - removing", h->GetName()));
+ remove = true;
+ }
+ else if (h->GetListOfFunctions()->GetEntries() <= 0) {
+ // AliWarning(Form("No fuctions associated with %s - removing",
+ // h->GetName()));
+ remove = true;
+ }
+ if (remove) {
TObjLink* keep = lnk->Next();
dists->Remove(lnk);
lnk = keep;
break;
}
// Give a warning if we're using fall-back
- if (ret == f.fFunctions.At(0))
+ if (ret == f.fFunctions.At(0)) {
AliWarning("Choosing fall-back 1 particle fit");
-
+ }
// Copy our result and return (the functions are owned by the fitter)
TF1* fret = new TF1(*ret);
return fret;
// Check that the reduced chi square isn't larger than 5
if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
- if (fDebug > 2)
+ if (fDebug > 2) {
AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
- chi2nuCut));
+ chi2nuCut)); }
ret = kFALSE;
}
Double_t re = e / v;
Double_t cut = relErrorCut * (i < kN ? 1 : 2);
if (re > cut) {
- if (fDebug > 2)
+ if (fDebug > 2) {
AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
n.Data(), r->ParName(i).c_str(),
r->ParName(i).c_str(), e, v,
100*(v == 0 ? 0 : e/v),
- 100*(cut)));
+ 100*(cut))); }
ret = kFALSE;
}
}
// Check that the last particle has a significant contribution
Double_t lastScale = r->Parameter(nPar-1);
if (lastScale <= 1e-7) {
- if (fDebug)
+ if (fDebug) {
AliWarning(Form("%s: %s=%9.6f<1e-7",
- n.Data(), r->ParName(nPar-1).c_str(), lastScale));
+ n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
ret = kFALSE;
}
}
}
+//__________________________________________________________________
+void
+AliFMDEnergyFitter::RingHistos::FindBestFits(TList* d,
+ AliFMDCorrELossFit& obj,
+ const TAxis& eta,
+ Double_t relErrorCut,
+ Double_t chi2nuCut,
+ Double_t minWeightCut)
+{
+ // Get our ring list from the output container
+ TList* l = GetOutputList(d);
+ if (!l) return;
+
+ // Get the energy distributions from the output container
+ TList* dists = static_cast<TList*>(l->FindObject("EDists"));
+ if (!dists) {
+ AliWarning(Form("Didn't find %s_EtaEDists in %s",
+ fName.Data(), l->GetName()));
+ l->ls();
+ return;
+ }
+ Int_t nBin = eta.GetNbins();
+
+ for (Int_t b = 1; b <= nBin; b++) {
+ TString n(Form(fgkEDistFormat, fName.Data(), b));
+ TH1D* dist = static_cast<TH1D*>(dists->FindObject(n));
+ // Ignore empty histograms altoghether
+ if (!dist || dist->GetEntries() <= 0) continue;
+
+ AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
+ relErrorCut,
+ chi2nuCut,
+ minWeightCut);
+ best->fDet = fDet;
+ best->fRing = fRing;
+ best->fBin = b; //
+ // Double_t eta = fAxis->GetBinCenter(b);
+ obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
+ }
+}
+
+//__________________________________________________________________
+AliFMDCorrELossFit::ELossFit*
+AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
+ Double_t relErrorCut,
+ Double_t chi2nuCut,
+ Double_t minWeightCut)
+{
+ TList* funcs = dist->GetListOfFunctions();
+ TIter next(funcs);
+ TF1* func = 0;
+ fFits.Clear();
+ Int_t i = 0;
+ // Info("FindBestFit", "%s", dist->GetName());
+ while ((func = static_cast<TF1*>(next()))) {
+ AliFMDCorrELossFit::ELossFit* fit =
+ new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
+ fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
+ }
+ fFits.Sort(false);
+ // fFits.Print();
+ return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
+}
+
+
//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::Output(TList* dir)
#include <TAxis.h>
#include <TList.h>
#include <TObjArray.h>
+#include <TClonesArray.h>
+#include "AliFMDCorrELossFit.h"
#include "AliForwardUtil.h"
class AliESDFMD;
class TFitResult;
* @param doFit Whether to do the fits or not
*/
void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
+ /**
+ * Set whether to make the corrections object on the output. Note,
+ * fits should be enable for this to have any effect.
+ *
+ * @param doMake If true (false is default), do make the corrections object.
+ */
+ void SetDoMakeObject(Bool_t doMake=kTRUE) { fDoMakeObject = doMake; }
/**
* Set how many particles we will try to fit at most to the data
*
* @param x Number of energy loss bins
*/
void SetNEbins(Int_t x) { fNEbins = x; }
- void SetMaxRelativeParameterError(Double_t e) { fMaxRelParError = e; }
- void SetMaxChi2PerNDF(Double_t c) { fMaxChi2PerNDF = c; }
+ void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
+ void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
+ void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
/**
* Set wheter to use increasing bin sizes
*
* @param dir Where the histograms are
*/
void Fit(TList* dir);
+ void MakeCorrectionsObject(TList* dir);
/**
* Define the output histograms. These are put in a sub list of the
* @param dbg Debug level
*/
void SetDebug(Int_t dbg=1);
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Internal data structure to keep track of the histograms
UShort_t minusBins,
Double_t relErrorCut,
Double_t chi2nuCut) const;
+ /**
+ * Find the best fits
+ *
+ * @param d Parent list
+ * @param obj Object to add fits to
+ * @param eta Eta axis
+ * @param relErrorCut Cut applied to relative error of parameter.
+ * Note, for multi-particle weights, the cut
+ * is loosend by a factor of 2
+ * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
+ * the reduced @f$\chi^2@f$
+ * @param minWeightCut Least valid @f$ a_i@f$
+ */
+ void FindBestFits(TList* d,
+ AliFMDCorrELossFit& obj,
+ const TAxis& eta,
+ Double_t relErrorCut,
+ Double_t chi2nuCut,
+ Double_t minWeightCut);
+ /**
+ * Find the best fit
+ *
+ * @param dist Histogram
+ * @param relErrorCut Cut applied to relative error of parameter.
+ * Note, for multi-particle weights, the cut
+ * is loosend by a factor of 2
+ * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
+ * the reduced @f$\chi^2@f$
+ * @param minWeightCut Least valid @f$ a_i@f$
+ *
+ * @return Best fit
+ */
+ AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist,
+ Double_t relErrorCut,
+ Double_t chi2nuCut,
+ Double_t minWeightCut);
/**
* Check the result of the fit. Returns true if
* - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
Int_t high,
Double_t val,
Double_t err) const;
- TH1D* fEDist; // Ring energy distribution
- TH1D* fEmpty; // Ring energy distribution for empty events
- TList fEtaEDists; // Energy distributions per eta bin.
- TList* fList;
- Int_t fDebug;
+ TH1D* fEDist; // Ring energy distribution
+ TH1D* fEmpty; // Ring energy distribution for empty events
+ TList fEtaEDists; // Energy distributions per eta bin.
+ TList* fList;
+ TClonesArray fFits;
+ Int_t fDebug;
ClassDef(RingHistos,1);
};
/**
UShort_t fNParticles; // Number of landaus to try to fit
UShort_t fMinEntries; // Minimum number of entries
UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
- Bool_t fDoFits; // Wheter to actually do the fits
+ Bool_t fDoFits; // Whether to actually do the fits
+ Bool_t fDoMakeObject; // Whether to make corrections object
TAxis fEtaAxis; // Eta axis
Double_t fMaxE; // Maximum energy loss to consider
Int_t fNEbins; // Number of energy loss bins
Bool_t fUseIncreasingBins; // Wheter to use increasing bin sizes
Double_t fMaxRelParError;// Relative error cut
Double_t fMaxChi2PerNDF; // chi^2/nu cit
+ Double_t fMinWeight; // Minimum weight value
Int_t fDebug; // Debug level
+
ClassDef(AliFMDEnergyFitter,1); //
};
--- /dev/null
+#include "AliFMDEnergyFitterTask.h"
+#include "AliLog.h"
+// #include "AliFMDAnaParameters.h"
+#include "AliESDEvent.h"
+#include "AliAODForwardMult.h"
+#include "AliForwardCorrectionManager.h"
+#include "AliAnalysisManager.h"
+#include <TH1.h>
+#include <TDirectory.h>
+#include <TTree.h>
+
+//====================================================================
+AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
+ : AliAnalysisTaskSE(),
+ fFirstEvent(true),
+ fEventInspector(),
+ fEnergyFitter(),
+ fList(0)
+{
+}
+
+//____________________________________________________________________
+AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name)
+ : AliAnalysisTaskSE(name),
+ fFirstEvent(true),
+ fEventInspector("event"),
+ fEnergyFitter("energy"),
+ fList(0)
+{
+ DefineOutput(1, TList::Class());
+}
+
+//____________________________________________________________________
+AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o)
+ : AliAnalysisTaskSE(o),
+ fFirstEvent(o.fFirstEvent),
+ fEventInspector(o.fEventInspector),
+ fEnergyFitter(o.fEnergyFitter),
+ fList(o.fList)
+{
+ DefineOutput(1, TList::Class());
+}
+
+//____________________________________________________________________
+AliFMDEnergyFitterTask&
+AliFMDEnergyFitterTask::operator=(const AliFMDEnergyFitterTask& o)
+{
+ AliAnalysisTaskSE::operator=(o);
+
+ fFirstEvent = o.fFirstEvent;
+ fEventInspector = o.fEventInspector;
+ fEnergyFitter = o.fEnergyFitter;
+ fList = o.fList;
+
+ return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
+{
+ fEventInspector.SetDebug(dbg);
+ fEnergyFitter.SetDebug(dbg);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::Init()
+{
+ fFirstEvent = true;
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::InitializeSubs()
+{
+
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // pars->Init(kTRUE);
+
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ fcm.Init(fEventInspector.GetCollisionSystem(),
+ fEventInspector.GetEnergy(),
+ fEventInspector.GetField(), 0);
+ TAxis eAxis(0,0,0);
+ TAxis vAxis(10,-10,10);
+ fEnergyFitter.Init(eAxis);
+ fEventInspector.Init(vAxis);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::UserCreateOutputObjects()
+{
+ fList = new TList;
+
+ fEventInspector.DefineOutput(fList);
+ fEnergyFitter.DefineOutput(fList);
+
+ PostData(1, fList);
+}
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::UserExec(Option_t*)
+{
+ // static Int_t cnt = 0;
+ // cnt++;
+ // Get the input data
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
+ if (!esd) {
+ AliWarning("No ESD event found for input event");
+ return;
+ }
+
+ // On the first event, initialize the parameters
+ if (fFirstEvent && esd->GetESDRun()) {
+ fEventInspector.ReadRunDetails(esd);
+
+ AliInfo(Form("Initializing with parameters from the ESD:\n"
+ " AliESDEvent::GetBeamEnergy() ->%f\n"
+ " AliESDEvent::GetBeamType() ->%s\n"
+ " AliESDEvent::GetCurrentL3() ->%f\n"
+ " AliESDEvent::GetMagneticField()->%f\n"
+ " AliESDEvent::GetRunNumber() ->%d\n",
+ esd->GetBeamEnergy(),
+ esd->GetBeamType(),
+ esd->GetCurrentL3(),
+ esd->GetMagneticField(),
+ esd->GetRunNumber()));
+
+
+
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // pars->SetParametersFromESD(esd);
+ // pars->PrintStatus();
+ fFirstEvent = false;
+
+ InitializeSubs();
+ }
+ Bool_t lowFlux = kFALSE;
+ UInt_t triggers = 0;
+ UShort_t ivz = 0;
+ Double_t vz = 0;
+ UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
+ if (found & AliFMDEventInspector::kNoEvent) return;
+ if (found & AliFMDEventInspector::kNoTriggers) return;
+ if (found & AliFMDEventInspector::kNoSPD) return;
+ if (found & AliFMDEventInspector::kNoFMD) return;
+ if (found & AliFMDEventInspector::kNoVertex) return;
+ if (found & AliFMDEventInspector::kBadVertex) return;
+
+ // Get FMD data
+ AliESDFMD* esdFMD = esd->GetFMDData();
+ // Do the energy stuff
+ if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
+ AliWarning("Energy fitter failed");
+ return;
+ }
+ PostData(1, fList);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::Terminate(Option_t*)
+{
+ TList* list = dynamic_cast<TList*>(GetOutputData(1));
+ if (!list) {
+ AliError(Form("No output list defined (%p)", GetOutputData(1)));
+ if (GetOutputData(1)) GetOutputData(1)->Print();
+ return;
+ }
+
+ // Get our histograms from the container
+ TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
+ TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
+ TH1I* hTriggers = 0;
+ if (!fEventInspector.FetchHistograms(list, hEventsTr,
+ hEventsTrVtx, hTriggers)) {
+ AliError(Form("Didn't get histograms from event selector "
+ "(hEventsTr=%p,hEventsTrVtx=%p)",
+ hEventsTr, hEventsTrVtx));
+ list->ls();
+ return;
+ }
+ fEnergyFitter.Fit(list);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitterTask::Print(Option_t*) const
+{
+}
+
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTERTASK_H
+#define ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTERTASK_H
+#include <AliAnalysisTaskSE.h>
+#include "AliForwardUtil.h"
+#include "AliFMDEventInspector.h"
+#include "AliFMDEnergyFitter.h"
+#include <AliESDFMD.h>
+#include <TH1I.h>
+class AliESDEvent;
+class TH2D;
+class TList;
+class TTree;
+
+
+/**
+ * @mainpage ALICE PWG2 Forward Multiplcity Analysis
+ */
+/**
+ * @defgroup pwg2_forward_analysis PWG2 Forward analysis
+ *
+ * Code to do the multiplicity analysis in the forward psuedo-rapidity
+ * regions
+ *
+ */
+/**
+ * Histogram and fit the energy loss distributions for the FMD
+ *
+ * @par Inputs:
+ * - AliESDEvent
+ *
+ * @par Outputs:
+ * - None
+ *
+ * @par Histograms
+ *
+ * @par Corrections used
+ * - None
+ *
+ * @ingroup pwg2_forward_analysis
+ *
+ */
+class AliFMDEnergyFitterTask : public AliAnalysisTaskSE
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param name Name of task
+ */
+ AliFMDEnergyFitterTask(const char* name);
+ /**
+ * Constructor
+ */
+ AliFMDEnergyFitterTask();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDEnergyFitterTask& operator=(const AliFMDEnergyFitterTask& o);
+ /**
+ * @{
+ * @name Interface methods
+ */
+ /**
+ * Initialize the task
+ *
+ */
+ virtual void Init();
+ /**
+ * Create output objects
+ *
+ */
+ virtual void UserCreateOutputObjects();
+ /**
+ * Process each event
+ *
+ * @param option Not used
+ */
+ virtual void UserExec(Option_t* option);
+ /**
+ * End of job
+ *
+ * @param option Not used
+ */
+ virtual void Terminate(Option_t* option);
+ /**
+ * @}
+ */
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
+ /**
+ * @{
+ * @name Access to sub-algorithms
+ */
+ /**
+ * Get reference to the EventInspector algorithm
+ *
+ * @return Reference to AliFMDEventInspector object
+ */
+ AliFMDEventInspector& GetEventInspector() { return fEventInspector; }
+ /**
+ * Get reference to the EnergyFitter algorithm
+ *
+ * @return Reference to AliFMDEnergyFitter object
+ */
+ AliFMDEnergyFitter& GetEnergyFitter() { return fEnergyFitter; }
+ /**
+ * @}
+ */
+ void SetDebug(Int_t dbg);
+protected:
+ /**
+ * Initialise the sub objects and stuff. Called on first event
+ *
+ */
+ virtual void InitializeSubs();
+
+ Bool_t fFirstEvent; // Whether the event is the first seen
+ AliFMDEventInspector fEventInspector; // Algorithm
+ AliFMDEnergyFitter fEnergyFitter; // Algorithm
+ TList* fList; // Output list
+
+ ClassDef(AliFMDEnergyFitterTask,1) // Forward multiplicity class
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliAODForwardMult.h"
+#include "AliForwardUtil.h"
#include <TH1.h>
#include <TList.h>
#include <TDirectory.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
//====================================================================
AliFMDEventInspector::AliFMDEventInspector()
fHEventsTr(0),
fHEventsTrVtx(0),
fHTriggers(0),
+ fHType(0),
fLowFluxCut(1000),
fMaxVzErr(0.1),
fList(0),
+ fEnergy(0),
+ fField(999),
+ fCollisionSystem(kUnknown),
fDebug(0)
{
}
fHEventsTr(0),
fHEventsTrVtx(0),
fHTriggers(0),
+ fHType(0),
fLowFluxCut(1000),
fMaxVzErr(0.1),
fList(0),
+ fEnergy(0),
+ fField(999),
+ fCollisionSystem(kUnknown),
fDebug(0)
{
}
fHEventsTr(o.fHEventsTr),
fHEventsTrVtx(o.fHEventsTrVtx),
fHTriggers(o.fHTriggers),
+ fHType(o.fHType),
fLowFluxCut(o.fMaxVzErr),
fMaxVzErr(o.fMaxVzErr),
fList(o.fList),
+ fEnergy(o.fEnergy),
+ fField(o.fField),
+ fCollisionSystem(o.fCollisionSystem),
fDebug(0)
{
}
if (fHEventsTr) delete fHEventsTr;
if (fHEventsTrVtx) delete fHEventsTrVtx;
if (fHTriggers) delete fHTriggers;
+ if (fHType) delete fHType;
if (fList) delete fList;
}
//____________________________________________________________________
fHEventsTr = o.fHEventsTr;
fHEventsTrVtx = o.fHEventsTrVtx;
fHTriggers = o.fHTriggers;
+ fHType = o.fHType;
fLowFluxCut = o.fLowFluxCut;
fMaxVzErr = o.fMaxVzErr;
fDebug = o.fDebug;
fList = (o.fList ? new TList : 0);
+ fEnergy = o.fEnergy;
+ fField = o.fField;
+ fCollisionSystem = o.fCollisionSystem;
if (fList) {
fList->SetName(GetName());
if (fHEventsTr) fList->Add(fHEventsTr);
if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
if (fHTriggers) fList->Add(fHTriggers);
+ if (fHType) fList->Add(fHType);
}
return *this;
}
fHTriggers->GetXaxis()->SetBinLabel(9, "spare1");
fHTriggers->GetXaxis()->SetBinLabel(10, "spare2");
fList->Add(fHTriggers);
+
+ fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)",
+ fLowFluxCut), 2, -.5, 1.5);
+ fHType->SetFillColor(kRed+1);
+ fHType->SetFillStyle(3001);
+ fHType->SetStats(0);
+ fHType->SetDirectory(0);
+ fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
+ fHType->GetXaxis()->SetBinLabel(2,"High-flux");
+ fList->Add(fHType);
+
}
//____________________________________________________________________
AliFMDEventInspector::Process(const AliESDEvent* event,
UInt_t& triggers,
Bool_t& lowFlux,
- Int_t& ivz,
+ UShort_t& ivz,
Double_t& vz)
{
// Check that we have an event
// Read trigger information from the ESD and store in AOD object
if (!ReadTriggers(event, triggers)) {
- if (fDebug > 2)
- AliWarning("Failed to read triggers from ESD");
+ if (fDebug > 2) {
+ AliWarning("Failed to read triggers from ESD"); }
return kNoTriggers;
}
// Check if this is a high-flux event
const AliMultiplicity* testmult = event->GetMultiplicity();
if (!testmult) {
- if (fDebug > 3)
- AliWarning("No central multiplicity object found");
+ if (fDebug > 3) {
+ AliWarning("No central multiplicity object found"); }
return kNoSPD;
}
lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
+
+ fHType->Fill(lowFlux ? 0 : 1);
// Check the FMD ESD data
if (!event->GetFMDData()) {
- if (fDebug > 3)
- AliWarning("No FMD data found in ESD");
+ if (fDebug > 3) {
+ AliWarning("No FMD data found in ESD"); }
return kNoFMD;
}
fHEventsTr->Fill(vz);
if (!vzOk) {
- if (fDebug > 3)
- AliWarning("Failed to read vertex from ESD");
+ if (fDebug > 3) {
+ AliWarning("Failed to read vertex from ESD"); }
return kNoVertex;
}
fHEventsTrVtx->Fill(vz);
// Get the vertex bin
- ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
- if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
- if (fDebug > 3)
+ ivz = fHEventsTr->GetXaxis()->FindBin(vz);
+ if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
+ if (fDebug > 3) {
AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
vz, fHEventsTr->GetXaxis()->GetXmin(),
- fHEventsTr->GetXaxis()->GetXmax()));
- ivz = -1;
+ fHEventsTr->GetXaxis()->GetXmax())); }
+ ivz = 0;
return kBadVertex;
}
return kOk;
// Get the vertex
const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
if (!vertex) {
- if (fDebug > 2)
- AliWarning("No SPD vertex found in ESD");
+ if (fDebug > 2) {
+ AliWarning("No SPD vertex found in ESD"); }
return kFALSE;
}
// Check that enough tracklets contributed
if(vertex->GetNContributors() <= 0) {
- if (fDebug > 2)
+ if (fDebug > 2) {
AliWarning(Form("Number of contributors to vertex is %d<=0",
- vertex->GetNContributors()));
+ vertex->GetNContributors())); }
vz = 0;
return kFALSE;
}
// Check that the uncertainty isn't too large
if (vertex->GetZRes() > fMaxVzErr) {
- if (fDebug > 2)
+ if (fDebug > 2) {
AliWarning(Form("Uncertaintity in Z of vertex is too large %f > %f",
- vertex->GetZRes(), fMaxVzErr));
+ vertex->GetZRes(), fMaxVzErr)); }
return kFALSE;
}
vz = vertex->GetZ();
return kTRUE;
}
+
+//____________________________________________________________________
+Bool_t
+AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
+{
+ fCollisionSystem =
+ AliForwardUtil::ParseCollisionSystem(esd->GetBeamType());
+ fEnergy =
+ AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,
+ 2 * esd->GetBeamEnergy());
+ fField =
+ AliForwardUtil::ParseMagneticField(esd->GetMagneticField());
+
+ if (fCollisionSystem == AliForwardUtil::kUnknown ||
+ fEnergy <= 0 ||
+ TMath::Abs(fField) > 10)
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliFMDEventInspector::Print(Option_t*) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
+ sNN.Strip(TString::kBoth, '0');
+ sNN.ReplaceAll("GeV", " GeV");
+ TString field(AliForwardUtil::MagneticFieldString(fField));
+ field.ReplaceAll("p", "+");
+ field.ReplaceAll("m", "-");
+ field.ReplaceAll("kG", " kG");
+
+ std::cout << ind << "AliFMDEventInspector: " << GetName() << '\n'
+ << ind << " Low flux cut: " << fLowFluxCut << '\n'
+ << ind << " Max(delta v_z): " << fMaxVzErr << " cm\n"
+ << ind << " System: "
+ << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
+ << ind << " CMS energy per nucleon: " << sNN << '\n'
+ << ind << " Field: " << field << std::endl;
+}
+
+
//
// EOF
//
kC,
kE
};
+ /**
+ * Collision systems
+ */
+ enum ECollisionSystem {
+ kUnknown,
+ kPP,
+ kPbPb
+ };
/**
* Constructor
*/
* Process the event
*
* @param event Input event
+ * @param first True on first event
* @param triggers On return, the triggers fired
* @param lowFlux On return, true if the event is considered a low-flux
* event (according to the setting of fLowFluxCut)
- * @param ivz On return, the found vertex bin (zero-based)
+ * @param ivz On return, the found vertex bin (1-based). A zero
+ * means outside of the defined vertex range
* @param vz On return, the z position of the interaction
*
* @return 0 (or kOk) on success, otherwise a bit mask of error codes
UInt_t Process(const AliESDEvent* event,
UInt_t& triggers,
Bool_t& lowFlux,
- Int_t& ivz,
+ UShort_t& ivz,
Double_t& vz);
/**
* Define the output histograms. These are put in a sub list of the
TH1I*& hEventsTr,
TH1I*& hEventsTrVtx,
TH1I*& hTriggers) const;
+ /**
+ * Read the collision system, collision energy, and L3 field setting
+ * from the ESD
+ *
+ * @param esd ESD to get information from
+ *
+ * @return true on success, false
+ */
+ Bool_t ReadRunDetails(const AliESDEvent* esd);
+ /**
+ * Get the collision system (one of the value in ECollisionSystem)
+ *
+ * @return Collision system
+ */
+ UShort_t GetCollisionSystem() const { return fCollisionSystem; }
+ /**
+ * Get the center of mass energy (per nucleon pair) in GeV
+ *
+ * @return center of mass energy (per nucleon pair) in GeV
+ */
+ UShort_t GetEnergy() const { return fEnergy; }
+ /**
+ * Get the magnetic field setting of the L3 magnet in kilo Gauss.
+ *
+ * @return The magnetic field setting
+ */
+ Short_t GetField() const { return fField; }
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Read the trigger information from the ESD event
TH1I* fHEventsTr; //! Histogram of events w/trigger
TH1I* fHEventsTrVtx; //! Events w/trigger and vertex
TH1I* fHTriggers; //! Triggers
+ TH1I* fHType; //! Type (low/high flux) of event
Int_t fLowFluxCut; // Low flux cut
Double_t fMaxVzErr; // Maximum error on v_z
TList* fList; //! Histogram container
+ UShort_t fEnergy; // CMS energy (per nucleon pair) [GeV]
+ Short_t fField; // L3 magnetic field [kG]
+ UShort_t fCollisionSystem; //
Int_t fDebug; // Debug level
ClassDef(AliFMDEventInspector,1); // Inspect the event
};
-
+// Calculate the total energy loss
#include "AliFMDHistCollector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
-#include "AliFMDAnaParameters.h"
+// #include "AliFMDAnaParameters.h"
+#include "AliForwardCorrectionManager.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TH1I.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TArrayI.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
ClassImp(AliFMDHistCollector)
#if 0
void
AliFMDHistCollector::Init(const TAxis& vtxAxis)
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
- Int_t nVz = vtxAxis.GetNbins();
+ UShort_t nVz = vtxAxis.GetNbins();
fFirstBins.Set(5*nVz);
fLastBins.Set(5*nVz);
// Find the eta bin ranges
- for (Int_t iVz = 0; iVz < nVz; iVz++) {
+ for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
// Find the first and last eta bin to use for each ring for
// each vertex bin. This is instead of using the methods
// provided by AliFMDAnaParameters
GetDetRing(iIdx, d, r);
// Get the background object
- TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
+ // TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
+ TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
Int_t nEta = bg->GetNbinsX();
Int_t first = nEta+1;
Int_t last = 0;
}
// Store the result for later use
- fFirstBins[iVz*5+iIdx] = first;
- fLastBins[iVz*5+iIdx] = last;
+ fFirstBins[(iVz-1)*5+iIdx] = first;
+ fLastBins[(iVz-1)*5+iIdx] = last;
} // for j
}
}
//____________________________________________________________________
void
-AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin,
+AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
Int_t& first, Int_t& last) const
{
first = 0;
last = 0;
- if (idx < 0) return;
- idx += vtxbin * 5;
+ if (idx < 0) return;
+ if (vtxbin <= 0) return;
+ idx += (vtxbin-1) * 5;
if (idx < 0 || idx >= fFirstBins.GetSize()) return;
//____________________________________________________________________
Int_t
-AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const
+AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
{
Int_t f, l;
GetFirstAndLast(idx,v,f,l);
//____________________________________________________________________
Int_t
-AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const
+AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
{
Int_t f, l;
GetFirstAndLast(idx,v,f,l);
//____________________________________________________________________
Int_t
AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
- Int_t bin, Int_t vtxbin) const
+ Int_t bin, UShort_t vtxbin) const
{
// Get the possibly overlapping histogram
Int_t other = -1;
}
//____________________________________________________________________
Int_t
-AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
+AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
{
UShort_t d = 0;
Char_t r = '\0';
//____________________________________________________________________
Bool_t
AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
- Int_t vtxbin,
+ UShort_t vtxbin,
TH2D& out)
{
for (UShort_t d=1; d<=3; d++) {
return kTRUE;
}
+//____________________________________________________________________
+void
+AliFMDHistCollector::Print(Option_t* /* option */) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
+ << ind << " # of cut bins: " << fNCutBins << '\n'
+ << ind << " Correction cut: " << fCorrectionCut << '\n'
+ << ind << " Bin ranges:\n" << ind << " v_z bin";
+ Int_t nVz = fFirstBins.fN / 5;
+ for (UShort_t iVz = 1; iVz <= nVz; iVz++)
+ std::cout << " | " << std::setw(7) << iVz;
+ std::cout << '\n' << ind << " / ring ";
+ for (UShort_t iVz = 1; iVz <= nVz; iVz++)
+ std::cout << "-+--------";
+ std::cout << std::endl;
+
+ for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
+ UShort_t d = 0;
+ Char_t r = 0;
+ GetDetRing(iIdx, d, r);
+
+ std::cout << ind << " FMD" << d << r << " ";
+ for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
+ Int_t first, last;
+ GetFirstAndLast(iIdx, iVz, first, last);
+ std::cout << " | " << std::setw(3) << first << "-"
+ << std::setw(3) << last;
+ }
+ std::cout << std::endl;
+ }
+}
//____________________________________________________________________
//
* Do the calculations
*
* @param hists Cache of histograms
- * @param vtxBin Vertex bin
+ * @param vtxBin Vertex bin (1 based)
* @param out Output histogram
*
* @return true on successs
*/
- virtual Bool_t Collect(AliForwardUtil::Histos& hists, Int_t vtxBin,
+ virtual Bool_t Collect(AliForwardUtil::Histos& hists, UShort_t vtxBin,
TH2D& out);
/**
* Set the number of extra bins (beyond the secondary map border)
* @param dbg Debug level
*/
void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Get the first and last eta bin to use for a given ring and vertex
*
* @param d Detector
* @param r Ring
- * @param vtxBin Vertex bin
+ * @param vtxBin Vertex bin (1 based)
* @param first On return, the first eta bin to use
* @param last On return, the last eta bin to use
*/
- virtual void GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxBin,
+ virtual void GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxBin,
Int_t& first, Int_t& last) const;
/**
* Get the first and last eta bin to use for a given ring and vertex
*
* @param idx Ring index as given by GetIdx
- * @param vtxBin Vertex bin
+ * @param vtxBin Vertex bin (1 based)
* @param first On return, the first eta bin to use
* @param last On return, the last eta bin to use
*/
- virtual void GetFirstAndLast(Int_t idx, Int_t vtxBin,
+ virtual void GetFirstAndLast(Int_t idx, UShort_t vtxBin,
Int_t& first, Int_t& last) const;
/**
* Get the first eta bin to use for a given ring and vertex
*
* @param d Detector
* @param r Ring
- * @param v vertex bin
+ * @param v vertex bin (1 based)
*
* @return First eta bin to use, or -1 in case of problems
*/
- Int_t GetFirst(UShort_t d, Char_t r, Int_t v) const;
+ Int_t GetFirst(UShort_t d, Char_t r, UShort_t v) const;
/**
* Get the first eta bin to use for a given ring and vertex
*
* @param idx Ring index as given by GetIdx
- * @param v vertex bin
+ * @param v vertex bin (1 based)
*
* @return First eta bin to use, or -1 in case of problems
*/
- Int_t GetFirst(Int_t idx, Int_t v) const;
+ Int_t GetFirst(Int_t idx, UShort_t v) const;
/**
* Get the last eta bin to use for a given ring and vertex
*
* @param d Detector
* @param r Ring
- * @param v vertex bin
+ * @param v vertex bin (1 based)
*
* @return Last eta bin to use, or -1 in case of problems
*/
- Int_t GetLast(UShort_t d, Char_t r, Int_t v) const;
+ Int_t GetLast(UShort_t d, Char_t r, UShort_t v) const;
/**
* Get the last eta bin to use for a given ring and vertex
*
* @param idx Ring index as given by GetIdx
- * @param v vertex bin
+ * @param v vertex bin (1 based)
*
* @return Last eta bin to use, or -1 in case of problems
*/
- Int_t GetLast(Int_t idx, Int_t v) const;
+ Int_t GetLast(Int_t idx, UShort_t v) const;
/**
* Get the detector and ring from the ring index
*
* @param d Detector
* @param r Ring
* @param e Eta bin
- * @param v Vertex bin
+ * @param v Vertex bin (1 based)
*
* @return Overlapping histogram index or -1
*/
- Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
/**
* Get the possibly overlapping histogram of eta bin @a e in
* detector and ring
*
* @param i Ring index
* @param e Eta bin
- * @param v Vertex bin
+ * @param v Vertex bin (1 based)
*
* @return Overlapping histogram index or -1
*/
- Int_t GetOverlap(Int_t i, Int_t e, Int_t v) const;
+ Int_t GetOverlap(Int_t i, Int_t e, UShort_t v) const;
/**
* Check if there's an overlapping histogram with this eta bin of
* the detector and ring
* @param d Detector
* @param r Ring
* @param e eta bin
- * @param v Vertex bin
+ * @param v Vertex bin (1 based)
*
* @return True if there's an overlapping histogram
*/
- Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
/**
* Check if there's an overlapping histogram with this eta bin of
* ring
*
* @return True if there's an overlapping histogram
*/
- Bool_t HasOverlap(Int_t i, Int_t e, Int_t v) const;
+ Bool_t HasOverlap(Int_t i, Int_t e, UShort_t v) const;
Int_t fNCutBins; // Number of additional bins to cut away
//____________________________________________________________________
inline void
-AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxbin,
+AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxbin,
Int_t& first, Int_t& last) const
{
GetFirstAndLast(GetIdx(d,r), vtxbin, first, last);
}
//____________________________________________________________________
inline Int_t
-AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, Int_t v) const
+AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, UShort_t v) const
{
return GetFirst(GetIdx(d,r), v);
}
//____________________________________________________________________
inline Int_t
-AliFMDHistCollector::GetLast(UShort_t d, Char_t r, Int_t v) const
+AliFMDHistCollector::GetLast(UShort_t d, Char_t r, UShort_t v) const
{
return GetLast(GetIdx(d, r), v);
}
//____________________________________________________________________
inline Bool_t
-AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const
+AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const
{
return GetOverlap(d,r,e,v) >= 0;
}
//____________________________________________________________________
inline Bool_t
-AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, Int_t v) const
+AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, UShort_t v) const
{
return GetOverlap(i,e,v) >= 0;
}
--- /dev/null
+#include "AliFMDMCDensityCalculator.h"
+#include <TMath.h>
+#include "AliForwardCorrectionManager.h"
+#include "AliFMDStripIndex.h"
+#include "AliMCEvent.h"
+// #include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+
+ClassImp(AliFMDMCDensityCalculator)
+#if 0
+; // For Emacs
+#endif
+
+
+//____________________________________________________________________
+AliFMDMCDensityCalculator&
+AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
+{
+ AliFMDDensityCalculator::operator=(o);
+ return *this;
+}
+
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCDensityCalculator::CalculateMC(const AliMCEvent& event,
+ AliForwardUtil::Histos& hists,
+ Double_t vz,
+ UShort_t vtxbin)
+{
+ Int_t nTracks = event.GetNumberOfTracks();
+ for (Int_t iTr = 0; iTr < nTracks; iTr++) {
+ AliMCParticle* particle =
+ static_cast<AliMCParticle*>(event.GetTrack(iTr));
+
+ // Check the returned particle
+ if (!particle) continue;
+
+ // Check if this charged and a primary
+ Bool_t isCharged = particle->Charge() != 0;
+ if (!isCharged) continue;
+
+ Int_t nTrRef = particle->GetNumberOfTrackReferences();
+ for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
+ AliTrackReference* ref = particle->GetTrackReference(iTrRef);
+
+ // Check existence
+ if (!ref) continue;
+
+ // Check that we hit an FMD element
+ if (ref->DetectorId() != AliTrackReference::kFMD)
+ continue;
+
+ // Get the detector coordinates
+ UShort_t d, s, t;
+ Char_t r;
+ AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
+
+ Double_t x = ref->X();
+ Double_t y = ref->Y();
+ Double_t z = ref->Z()-vz;
+ Double_t rr = TMath::Sqrt(x*x+y*y);
+ Double_t phi = TMath::ATan2(y,x);
+ Double_t theta= TMath::ATan2(rr,z);
+ Double_t eta = -TMath::Log(TMath::Tan(theta/2));
+
+ Float_t c = Correction(d,r,s,t,vtxbin,eta,false);
+ fCorrections->Fill(c);
+
+ TH2D* h = hists.Get(d,r);
+ h->Fill(eta,phi, 1 * c);
+ }
+ }
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCDensityCalculator::Calculate(const AliESDFMD&,
+ AliForwardUtil::Histos&,
+ UShort_t,
+ Bool_t)
+{
+ AliWarning("Method Calculate disabled for this class. If you need this, "
+ "make an AliFMDDensityCalculator object instead");
+ return kFALSE;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDMCDENSITYCALCULATOR_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDMCDENSITYCALCULATOR_H
+#include "AliFMDDensityCalculator.h"
+#include <TList.h>
+#include "AliForwardUtil.h"
+class AliMCEvent;
+class TH2D;
+class TH1D;
+
+/**
+ * This class calculates the inclusive charged particle density
+ * in each for the 5 FMD rings based on the MC truth.
+ *
+ * @par Input:
+ * - AliMCEvent MC truth event infromation
+ *
+ * @par Output:
+ * - None
+ *
+ * @par Corrections used:
+ * - None
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDMCDensityCalculator : public AliFMDDensityCalculator
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDMCDensityCalculator() : AliFMDDensityCalculator() {}
+ /**
+ * Constructor
+ *
+ * @param name Name of object
+ */
+ AliFMDMCDensityCalculator(const char* name)
+ : AliFMDDensityCalculator(name)
+ {}
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDMCDensityCalculator(const AliFMDMCDensityCalculator& o)
+ : AliFMDDensityCalculator(o)
+ {}
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDMCDensityCalculator() {}
+ /**
+ * Assignement operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDMCDensityCalculator& operator=(const AliFMDMCDensityCalculator&);
+ /**
+ * Do the calculations
+ *
+ * @param fmd AliESDFMD object (possibly) corrected for sharing
+ * @param hists Histogram cache
+ * @param vtxBin Vertex bin
+ * @param lowFlux Low flux flag.
+ *
+ * @return true on successs
+ */
+ virtual Bool_t Calculate(const AliESDFMD& fmd,
+ AliForwardUtil::Histos& hists,
+ UShort_t vtxBin, Bool_t lowFlux);
+ /**
+ * Calculate the charged particle density from the MC track references.
+ *
+ * @param event MC event
+ * @param hists Histograms to fill
+ * @param vz Interaction z coordinate @f$ v_z@f$
+ * @param vtxBin bin corresponding to @f$ v_z@f$
+ *
+ * @return true on success
+ */
+ virtual Bool_t CalculateMC(const AliMCEvent& event,
+ AliForwardUtil::Histos& hists,
+ Double_t vz,
+ UShort_t vtxBin);
+protected:
+
+ ClassDef(AliFMDMCDensityCalculator,1); // Calculate Nch density
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
#include <TList.h>
#include <TH1.h>
#include <TMath.h>
-#include "AliFMDAnaParameters.h"
+#include "AliForwardCorrectionManager.h"
+// #include "AliFMDAnaParameters.h"
#include <AliLog.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
ClassImp(AliFMDSharingFilter)
#if 0
AliFMDSharingFilter::AliFMDSharingFilter()
: TNamed(),
fRingHistos(),
- fLowCut(0.3),
+ fLowCut(0.),
fCorrectAngles(kFALSE),
+ fNXi(1),
fDebug(0)
{}
AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
: TNamed("fmdSharingFilter", title),
fRingHistos(),
- fLowCut(0.3),
+ fLowCut(0.),
fCorrectAngles(kFALSE),
+ fNXi(1),
fDebug(0)
{
fRingHistos.Add(new RingHistos(1, 'I'));
fRingHistos(),
fLowCut(o.fLowCut),
fCorrectAngles(o.fCorrectAngles),
+ fNXi(o.fNXi),
fDebug(o.fDebug)
{
TIter next(&o.fRingHistos);
fLowCut = o.fLowCut;
fCorrectAngles = o.fCorrectAngles;
+ fNXi = o.fNXi;
fDebug = o.fDebug;
fRingHistos.Delete();
{
output.Clear();
TIter next(&fRingHistos);
- RingHistos* o = 0;
+ RingHistos* o = 0;
+ Double_t lowCut = GetLowCut();
while ((o = static_cast<RingHistos*>(next())))
o->Clear();
Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
lowFlux,d,r,s,t,
usedPrev,usedThis);
- if (mergedEnergy > fLowCut) histos->Incr();
+ if (mergedEnergy > lowCut) histos->Incr();
histos->fAfter->Fill(mergedEnergy);
output.SetMultiplicity(d,r,s,t,mergedEnergy);
mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
return mult;
}
+//_____________________________________________________________________
+Double_t
+AliFMDSharingFilter::GetLowCut() const
+{
+ if (fLowCut > 0) return fLowCut;
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ AliFMDCorrELossFit* fits = fcm.GetELossFit();
+ return fits->GetLowCut();
+}
//_____________________________________________________________________
Double_t
AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+
// Get the high cut. The high cut is defined as the
// most-probably-value peak found from the energy distributions, minus
// 2 times the width of the corresponding Landau.
- Double_t mpv = pars->GetMPV(d,r,eta);
- Double_t w = pars->GetSigma(d,r,eta);
- if (mpv > 100)
- AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
- Double_t highCut = (mpv - 2 * w);
+ AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
+ Double_t delta = 1024;
+ Double_t xi = 1024;
+ if (!fit) {
+ AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
+ }
+ else {
+ delta = fit->fDelta;
+ xi = fit->fXi;
+ }
+
+ if (delta > 100) {
+ AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
+ }
+ Double_t highCut = (delta - fNXi * xi);
return highCut;
}
{
// IF the multiplicity is very large, or below the cut, reject it, and
// flag it as candidate for sharing
- if(mult > 12 || mult < fLowCut) {
+ Double_t lowCut = GetLowCut();
+ if(mult > 12 || mult < lowCut) {
usedThis = kFALSE;
usedPrev = kFALSE;
return 0;
//analyse and perform sharing on one strip
- // Calculate the total energy loss
+ // Get the high cut
Double_t highCut = GetHighCut(d, r, eta);
// If this signal is smaller than the next, and the next signal is smaller
// If the previous signal was larger than the low cut, and
// the previous was smaller than high cut, and the previous signal
// isn't marked as used, then add it's energy to this signal
- if (prevE > fLowCut &&
+ if (prevE > lowCut &&
prevE < highCut &&
!usedPrev)
totalE += prevE;
// If the next signal is larger than the low cut, and
// the next signal is smaller than the low cut, then add the next signal
// to this, and marked the next signal as used
- if(nextE > fLowCut && nextE < highCut ) {
+ if(nextE > lowCut && nextE < highCut ) {
totalE += nextE;
usedThis = kTRUE;
}
o->Output(d);
}
}
+//____________________________________________________________________
+void
+AliFMDSharingFilter::Print(Option_t* /*option*/) const
+{
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ std::cout << ind << "AliFMDSharingFilter: " << GetName() << '\n'
+ << ind << " Low cut: " << fLowCut << '\n'
+ << ind << " N xi factor: " << fNXi << '\n'
+ << ind << " Use corrected angles: "
+ << (fCorrectAngles ? "yes" : "no") << std::endl;
+}
//====================================================================
AliFMDSharingFilter::RingHistos::RingHistos()
*
* @param lowCut Low cut
*/
- void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
+ void SetLowCut(Double_t lowCut=0) { fLowCut = lowCut; }
+ /**
+ * Reset the low cut for sharing to use the fit range lower cut
+ *
+ */
+ void UnsetLowCut() { fLowCut = 0; }
/**
* Set the debug level. The higher the value the more output
*
* signals are always angle corrected.
*/
void UseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
+ /**
+ * Set the number of landau width to subtract from the most probably
+ * value to get the high cut for the merging algorithm.
+ *
+ * @param n Number of @f$ \xi@f$
+ */
+ void SetNXi(Short_t n) { fNXi = n; }
/**
* Filter the input AliESDFMD object
*
* @param dir Directory to add to
*/
void DefineOutput(TList* dir);
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* Internal data structure to keep track of the histograms
* 2 times the width of the corresponding Landau.
*/
virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta) const;
+ /**
+ * Get the low cut. Normally, the low cut is taken to be the lower
+ * value of the fit range used when generating the energy loss fits.
+ * However, if fLowCut is set (using SetLowCit) to a value greater
+ * than 0, then that value is used.
+ */
+ virtual Double_t GetLowCut() const;
TList fRingHistos; // List of histogram containers
Double_t fLowCut; // Low cut on sharing
Bool_t fCorrectAngles; // Whether to work on angle corrected signals
+ Short_t fNXi; // Number of xi's from Delta to stop merging
Int_t fDebug; // Debug level
ClassDef(AliFMDSharingFilter,1); //
--- /dev/null
+#include "AliForwardCorrectionManager.h"
+#include "AliForwardUtil.h"
+#include <TString.h>
+#include <AliLog.h>
+#include <TFile.h>
+#include <TSystem.h>
+
+
+//____________________________________________________________________
+AliForwardCorrectionManager* AliForwardCorrectionManager::fgInstance = 0;
+const char* AliForwardCorrectionManager::fgkSecondaryMapSkel = "secondary";
+const char* AliForwardCorrectionManager::fgkDoubleHitSkel = "doublehit";
+const char* AliForwardCorrectionManager::fgkELossFitsSkel = "elossfits";
+const char* AliForwardCorrectionManager::fgkVertexBiasSkel = "vertexbias";
+const char* AliForwardCorrectionManager::fgkMergingEffSkel = "merging";
+
+#define PREFIX "$(ALICE_ROOT)/PWG2/FORWARD/corrections/"
+
+//____________________________________________________________________
+AliForwardCorrectionManager& AliForwardCorrectionManager::Instance()
+{
+ if (!fgInstance) fgInstance= new AliForwardCorrectionManager;
+ return *fgInstance;
+}
+
+//____________________________________________________________________
+AliForwardCorrectionManager::AliForwardCorrectionManager()
+ : TObject(),
+ fInit(kFALSE),
+ fSys(0),
+ fSNN(0),
+ fField(999),
+ fELossFitsPath(PREFIX "ELossFits"),
+ fMergingEffPath(PREFIX "MergingEfficiency"),
+ fSecondaryMapPath(PREFIX "SecondaryMap"),
+ fDoubleHitPath(PREFIX "DoubleHit"),
+ fVertexBiasPath(PREFIX "VertexBias"),
+ fELossFit(0),
+ fSecondaryMap(0),
+ fDoubleHit(0),
+ fVertexBias(0),
+ fMergingEfficiency(0)
+{
+}
+//____________________________________________________________________
+AliForwardCorrectionManager::AliForwardCorrectionManager(const AliForwardCorrectionManager& o)
+ : TObject(o),
+ fInit(o.fInit),
+ fSys(o.fSys),
+ fSNN(o.fSNN),
+ fField(o.fField),
+ fELossFitsPath(o.fELossFitsPath),
+ fMergingEffPath(o.fMergingEffPath),
+ fSecondaryMapPath(o.fSecondaryMapPath),
+ fDoubleHitPath(o.fDoubleHitPath),
+ fVertexBiasPath(o.fVertexBiasPath),
+ fELossFit(o.fELossFit),
+ fSecondaryMap(o.fSecondaryMap),
+ fDoubleHit(o.fDoubleHit),
+ fVertexBias(o.fVertexBias),
+ fMergingEfficiency(o.fMergingEfficiency)
+
+{
+}
+//____________________________________________________________________
+AliForwardCorrectionManager&
+AliForwardCorrectionManager::operator=(const AliForwardCorrectionManager& o)
+{
+ fInit = o.fInit;
+ fSys = o.fSys;
+ fSNN = o.fSNN;
+ fField = o.fField;
+ fELossFitsPath = o.fELossFitsPath;
+ fMergingEffPath = o.fMergingEffPath;
+ fSecondaryMapPath = o.fSecondaryMapPath;
+ fDoubleHitPath = o.fDoubleHitPath;
+ fVertexBiasPath = o.fVertexBiasPath;
+ fELossFit = o.fELossFit;
+ fSecondaryMap = o.fSecondaryMap;
+ fDoubleHit = o.fDoubleHit;
+ fVertexBias = o.fVertexBias;
+ fMergingEfficiency= o.fMergingEfficiency;
+ return *this;
+}
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::Init(const char* cms,
+ Float_t sNN,
+ Float_t field,
+ Bool_t mc,
+ UInt_t what,
+ Bool_t force)
+{
+ UShort_t col = AliForwardUtil::ParseCollisionSystem(cms);
+ return Init(col,
+ AliForwardUtil::ParseCenterOfMassEnergy(col, sNN),
+ AliForwardUtil::ParseMagneticField(field),
+ mc, what, force);
+}
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::Init(UShort_t cms,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc,
+ UInt_t what,
+ Bool_t force)
+{
+ if (force) fInit = kFALSE;
+ if (fInit) return kTRUE;
+
+ Bool_t ret = kTRUE;
+ if (fSys == cms && TMath::Abs(fSNN - sNN) < 10 && fField == field) {
+ // We're already initialised for these settings - do nothing and return
+ fInit = kTRUE;
+ return ret;
+ }
+ // Set cached parameters
+ fSys = cms;
+ fSNN = sNN;
+ fField = field;
+
+ // Read secondary map if requested
+ if (what & kSecondaryMap) {
+ if (!ReadSecondaryMap(cms, sNN, field)) {
+ AliWarning(Form("Failed to read in secondary map for "
+ "cms=%d, sNN=%dGeV, field=%dkG", cms, sNN, field));
+ ret = kFALSE;
+ }
+ }
+ // Read double hit if requested
+ if (what & kDoubleHit) {
+ if (!ReadDoubleHit(cms, sNN, field)) {
+ AliWarning(Form("Failed to read in double hit correction for "
+ "cms=%d, sNN=%dGeV, field=%dkG", cms, sNN, field));
+ ret = kFALSE;
+ }
+ }
+ // Read energy loss fits if requested
+ if (what & kELossFits) {
+ if (!ReadELossFits(cms, sNN, field, mc)) {
+ AliWarning(Form("Failed to read in energy loss fits for "
+ "cms=%d, sNN=%dGeV, field=%dkG, %s",
+ cms, sNN, field, mc ? "MC" : "real"));
+ ret = kFALSE;
+ }
+ }
+ // Read event selection efficiencies if requested
+ if (what & kVertexBias) {
+ if (!ReadVertexBias(cms, sNN, field)) {
+ AliWarning(Form("Failed to read in event selection efficiency for "
+ "cms=%d, sNN=%dGeV, field=%dkG", cms, sNN, field));
+ ret = kFALSE;
+ }
+ }
+ // Read merging efficiencies if requested
+ if (what & kMergingEfficiency) {
+ if (!ReadMergingEfficiency(cms, sNN, field)) {
+ AliWarning(Form("Failed to read in hit merging efficiency for "
+ "cms=%d, sNN=%dGeV, field=%dkG",
+ cms, sNN, field));
+ ret = kFALSE;
+ }
+ }
+ fInit = kTRUE;
+ return ret;
+}
+
+//____________________________________________________________________
+TString
+AliForwardCorrectionManager::GetFileName(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const
+{
+ TString fname = "";
+ fname = GetObjectName(what);
+ fname.Append(Form("_%s_%04dGeV_%c%1dkG_%s.root",
+ AliForwardUtil::CollisionSystemString(sys),
+ sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field),
+ (mc ? "MC" : "real")));
+ return fname;
+}
+//____________________________________________________________________
+TString
+AliForwardCorrectionManager::GetFileName(ECorrection what) const
+{
+ if (!fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return "";
+ }
+ return GetFileName(what, fSys, fSNN, fField, false);
+}
+
+//____________________________________________________________________
+const Char_t*
+AliForwardCorrectionManager::GetFileDir(ECorrection what) const
+{
+ if (what & kSecondaryMap) return fSecondaryMapPath;
+ else if (what & kDoubleHit) return fDoubleHitPath;
+ else if (what & kELossFits) return fELossFitsPath;
+ else if (what & kVertexBias) return fVertexBiasPath;
+ else if (what & kMergingEfficiency) return fMergingEffPath;
+
+ AliWarning(Form("Unknown correction: %d", what));
+ return 0;
+}
+
+//____________________________________________________________________
+TString
+AliForwardCorrectionManager::GetFilePath(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const
+{
+ TString path = "";
+ const Char_t* dir = GetFileDir(what);
+ if (!dir) return path;
+
+ TString fname(GetFileName(what, sys, sNN, field, mc));
+ if (fname.IsNull()) return path;
+
+ path = gSystem->ConcatFileName(gSystem->ExpandPathName(dir), fname);
+
+ return path;
+}
+//____________________________________________________________________
+TString
+AliForwardCorrectionManager::GetFilePath(ECorrection what) const
+{
+ if (!fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return "";
+ }
+ return GetFilePath(what, fSys, fSNN, fField, false);
+}
+
+//____________________________________________________________________
+TFile*
+AliForwardCorrectionManager::GetFile(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc,
+ Bool_t rw,
+ Bool_t newfile) const
+{
+ TString path = GetFilePath(what, sys, sNN, field, mc);
+ if (path.IsNull()) return 0;
+
+ TString opt;
+ if (newfile) opt="RECREATE";
+ else {
+ if (gSystem->AccessPathName(path.Data(),
+ (rw ? kWritePermission : kReadPermission))) {
+ AliWarning(Form("file %s cannot be found or insufficient permissions",
+ path.Data()));
+ return 0;
+ }
+ opt=(rw ? "UPDATE" : "READ");
+ }
+ TFile* file = TFile::Open(path.Data(), opt.Data());
+ if (!file) {
+ AliWarning(Form("file %s cannot be opened in mode %s",
+ path.Data(), opt.Data()));
+ return 0;
+ }
+ return file;
+}
+//____________________________________________________________________
+TFile*
+AliForwardCorrectionManager::GetFile(ECorrection what) const
+{
+ if (!fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return 0;
+ }
+ return GetFile(what, fSys, fSNN, fField, false);
+}
+
+//____________________________________________________________________
+const Char_t*
+AliForwardCorrectionManager::GetObjectName(ECorrection what) const
+{
+ if (what & kSecondaryMap) return fgkSecondaryMapSkel;
+ else if (what & kDoubleHit) return fgkDoubleHitSkel;
+ else if (what & kELossFits) return fgkELossFitsSkel;
+ else if (what & kVertexBias) return fgkVertexBiasSkel;
+ else if (what & kMergingEfficiency) return fgkMergingEffSkel;
+ return 0;
+}
+
+//____________________________________________________________________
+TObject*
+AliForwardCorrectionManager::CheckObject(TFile* file, ECorrection what) const
+{
+ TObject* o = file->Get(GetObjectName(what));
+ if (!o) {
+ AliWarning(Form("Object %s not found in %s",
+ GetObjectName(what), file->GetName()));
+ file->Close();
+ return 0;
+ }
+ return o;
+}
+
+//____________________________________________________________________
+TObject*
+AliForwardCorrectionManager::GetObject(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const
+{
+ TFile* file = GetFile(what, sys, sNN, field, mc, false, false);
+ if (!file) return 0;
+
+ return CheckObject(file, what);
+}
+//____________________________________________________________________
+TObject*
+AliForwardCorrectionManager::GetObject(ECorrection what) const
+{
+ if (!fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return 0;
+ }
+ return GetObject(what, fSys, fSNN, fField, false);
+}
+
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::ReadSecondaryMap(UShort_t sys, UShort_t sNN,
+ Short_t field)
+{
+ if (fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return kFALSE;
+ }
+
+ TObject* o = GetObject(kSecondaryMap, sys, sNN, field, false);
+ if (!o) return kFALSE;
+
+ fSecondaryMap = dynamic_cast<AliFMDCorrSecondaryMap*>(o);
+ if (!fSecondaryMap) {
+ AliWarning(Form("Object %s (%p) is not an AliFMDCorrSecondaryMap object, "
+ "but %s", fgkSecondaryMapSkel, o, o->ClassName()));
+ return kFALSE;
+ }
+
+ // file->Close();
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::ReadDoubleHit(UShort_t sys, UShort_t sNN,
+ Short_t field)
+{
+ if (fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return kFALSE;
+ }
+
+ TObject* o = GetObject(kDoubleHit, sys, sNN, field, false);
+ if (!o) return kFALSE;
+
+ fDoubleHit = dynamic_cast<AliFMDCorrDoubleHit*>(o);
+ if (!fDoubleHit) {
+ AliWarning(Form("Object %s (%p) is not an AliFMDCorrDoubleHit object, "
+ "but %s", fgkDoubleHitSkel, o, o->ClassName()));
+ return kFALSE;
+ }
+
+ // file->Close();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::ReadELossFits(UShort_t sys, UShort_t sNN,
+ Short_t field, Bool_t mc)
+{
+ if (fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return kFALSE;
+ }
+
+ TObject* o = GetObject(kELossFits, sys, sNN, field, mc);
+ if (!o) return kFALSE;
+
+ fELossFit = dynamic_cast<AliFMDCorrELossFit*>(o);
+ if (!fELossFit) {
+ AliWarning(Form("Object %s (%p) is not an AliFMDCorrELossFit object, "
+ "but %s", fgkELossFitsSkel, o, o->ClassName()));
+ return kFALSE;
+ }
+
+ // file->Close();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::ReadVertexBias(UShort_t sys,
+ UShort_t sNN,
+ Short_t field)
+{
+ if (fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return kFALSE;
+ }
+
+ TObject* o = GetObject(kVertexBias, sys, sNN, field, false);
+ if (!o) return kFALSE;
+
+ fVertexBias = dynamic_cast<AliFMDCorrVertexBias*>(o);
+ if (!fVertexBias) {
+ AliWarning(Form("Object %s (%p) is not an AliFMDCorrVertexBias object, "
+ "but %s", fgkVertexBiasSkel, o, o->ClassName()));
+ return kFALSE;
+ }
+
+ // file->Close();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliForwardCorrectionManager::ReadMergingEfficiency(UShort_t sys,
+ UShort_t sNN,
+ Short_t field)
+{
+ if (fInit) {
+ AliWarning("Corrections manager initialised, do a forced Init(...)");
+ return kFALSE;
+ }
+
+ TObject* o = GetObject(kMergingEfficiency, sys, sNN, field, false);
+ if (!o) return kFALSE;
+
+ fMergingEfficiency = dynamic_cast<AliFMDCorrMergingEfficiency*>(o);
+ if (!fMergingEfficiency) {
+ AliWarning(Form("Object %s (%p) is not an AliFMDCorrMergingEfficiency "
+ "object, but %s", fgkMergingEffSkel, o, o->ClassName()));
+ return kFALSE;
+ }
+
+ // file->Close();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFORWARDCORRECTIONMANAGER_H
+#define ALIROOT_PWG2_FORWARD_ALIFORWARDCORRECTIONMANAGER_H
+#include <TObject.h>
+#include "AliFMDCorrELossFit.h"
+#include "AliFMDCorrSecondaryMap.h"
+#include "AliFMDCorrDoubleHit.h"
+#include "AliFMDCorrVertexBias.h"
+#include "AliFMDCorrMergingEfficiency.h"
+#include <TString.h>
+class TFile;
+
+class AliForwardCorrectionManager : public TObject
+{
+public:
+ /**
+ * Enumeration of things that can be read in
+ */
+ enum ECorrection {
+ kSecondaryMap = 0x01,
+ kELossFits = 0x02,
+ kVertexBias = 0x04,
+ kMergingEfficiency = 0x08,
+ kDoubleHit = 0x10,
+ kAll = (kSecondaryMap|
+ kELossFits|
+ kVertexBias|
+ kMergingEfficiency|
+ kDoubleHit)
+ };
+ /**
+ * Access to the singleton object
+ *
+ * @return Reference to the singleton object
+ */
+ static AliForwardCorrectionManager& Instance();
+ /**
+ * Read in corrections based on the parameters given
+ *
+ * @param collisionSystem Collision system
+ * @param cmsNN Center of mass energy per nuclean pair [GeV]
+ * @param field Magnetic field setting [kG]
+ * @param what What to read in.
+ * @param force Force (re-)reading of specified things
+ *
+ * @return
+ */
+ Bool_t Init(UShort_t collisionSystem,
+ UShort_t cmsNN,
+ Short_t field,
+ Bool_t mc=false,
+ UInt_t what=kAll,
+ Bool_t force=false);
+ Bool_t Init(const char* collisionSystem,
+ Float_t cmsNN,
+ Float_t field,
+ Bool_t mc=false,
+ UInt_t what=kAll,
+ Bool_t force=false);
+ /**
+ * Get the eta axis
+ *
+ * @return Eta axis or null
+ */
+ const TAxis* GetEtaAxis() const;
+ /**
+ * Get the vertex axis
+ *
+ * @return The vertex axis or null
+ */
+ const TAxis* GetVertexAxis() const;
+ /**
+ * Get the energy loss fit correction object.
+ *
+ * @return Get the energy loss fits corrections object or null pointer
+ */
+ AliFMDCorrELossFit* GetELossFit() const { return fELossFit; }
+ /**
+ * Get the secondary correction map
+ *
+ * @return Get the secondary correction map object or null
+ */
+ AliFMDCorrSecondaryMap* GetSecondaryMap() const { return fSecondaryMap; }
+ /**
+ * Get the double hit correction object
+ *
+ * @return Get the double hit correction object or null
+ */
+ AliFMDCorrDoubleHit* GetDoubleHit() const { return fDoubleHit; }
+ /**
+ * Get the vertex bias correction object
+ *
+ * @return Get the vertex bias correction object or null
+ */
+ AliFMDCorrVertexBias* GetVertexBias() const { return fVertexBias; }
+ AliFMDCorrMergingEfficiency* GetMergingEfficiency() const
+ {
+ return fMergingEfficiency;
+ }
+ /**
+ * @{
+ * @name Path, file, and object access utilities
+ */
+ /**
+ * Get the path to the specified object
+ *
+ * @param what Which stuff to get the path for
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ * @param mc Whether the correction objects should be valid for MC
+ *
+ * @return The full path or null
+ */
+ TString GetFileName(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const;
+ TString GetFileName(ECorrection what) const;
+ /**
+ * Get the path to the specified object
+ *
+ * @param what Which stuff to get the path for
+ *
+ * @return The full path or null
+ */
+ const Char_t* GetFileDir(ECorrection what) const;
+ /**
+ * Get the path to the specified object
+ *
+ * @param what Which stuff to get the path for
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ * @param mc Whether the correction objects should be valid for MC
+ *
+ * @return The full path or null
+ */
+ TString GetFilePath(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const;
+ TString GetFilePath(ECorrection what) const;
+ /**
+ * Open the file that contains the correction object specified
+ *
+ * @param what Which stuff to get the path for
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ * @param mc Whether the correction objects should be valid for MC
+ * @param rw Whether to open the file in read/write
+ * @param newfile Wheter to make the file if it doesn't exist
+ *
+ * @return The file that contains the correction object or null
+ */
+ TFile* GetFile(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc=false,
+ Bool_t rw=false,
+ Bool_t newfile=false) const;
+ TFile* GetFile(ECorrection what) const;
+ /**
+ * Get the object name corresponding to correction type
+ *
+ * @param what Correction
+ *
+ * @return Object name or null
+ */
+ const Char_t* GetObjectName(ECorrection what) const;
+ /**
+ * Check if the specified objet exists in the file, and return it
+ *
+ * @param file File to query
+ * @param what Correction type
+ *
+ * @return Object found, or null
+ */
+ TObject* CheckObject(TFile* file, ECorrection what) const;
+ /**
+ * Get the path to the specified object
+ *
+ * @param what Which stuff to get the path for
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ * @param mc Whether the correction objects should be valid for MC
+ *
+ * @return The full path or null
+ */
+ TObject* GetObject(ECorrection what,
+ UShort_t sys,
+ UShort_t sNN,
+ Short_t field,
+ Bool_t mc) const;
+ TObject* GetObject(ECorrection what) const;
+ /*
+ * @}
+ */
+private:
+ /**
+ * Default constructor
+ */
+ AliForwardCorrectionManager();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliForwardCorrectionManager(const AliForwardCorrectionManager& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliForwardCorrectionManager& operator=(const AliForwardCorrectionManager& o);
+ /**
+ * @{
+ * @name Read in corrections
+ */
+ /**
+ * Read in the secondary map
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t ReadSecondaryMap(UShort_t sys, UShort_t sNN, Short_t field);
+ /**
+ * Read in the double hit correction
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t ReadDoubleHit(UShort_t sys, UShort_t sNN, Short_t field);
+ /**
+ * Read in the energy loss fits
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ * @param mc Whether the correction objects should be valid for MC
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t ReadELossFits(UShort_t sys, UShort_t sNN, Short_t field, Bool_t mc);
+ /**
+ * Read in the event selection efficiency
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t ReadVertexBias(UShort_t sys, UShort_t sNN, Short_t field);
+ /**
+ * Read in the merging efficiency
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy [GeV]
+ * @param field Magnetic field in the L3 magnet [kG]
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t ReadMergingEfficiency(UShort_t sys, UShort_t sNN, Short_t field);
+ /*
+ * @}
+ */
+
+ /** Static singleton instance */
+ static AliForwardCorrectionManager* fgInstance;
+ Bool_t fInit; // whether we have been initialised
+ UShort_t fSys; // Collision System
+ UShort_t fSNN; // Collision energy per nucleon (GeV)
+ Short_t fField; // L3 magnetic field (kG)
+
+ /**
+ * @{
+ * @name Paths
+ */
+ TString fELossFitsPath; // Path to energy loss fit correction
+ TString fMergingEffPath; // Path to sharing efficiency correction
+ TString fSecondaryMapPath; // Path to secondary efficiency correction
+ TString fDoubleHitPath; // Path to double hit correction
+ TString fVertexBiasPath; // Path to event selection efficiency correction
+ /*
+ * @}
+ */
+ /**
+ * @{
+ * @name Object name
+ */
+ static const Char_t* fgkSecondaryMapSkel; // Name of correction object
+ static const Char_t* fgkDoubleHitSkel; // Name of correction object
+ static const Char_t* fgkELossFitsSkel; // Name of correction object
+ static const Char_t* fgkVertexBiasSkel; // Name of correction object
+ static const Char_t* fgkMergingEffSkel; // Name of correction object
+ /*
+ * @}
+ */
+ /**
+ * @{
+ * @name Correction objects
+ */
+ AliFMDCorrELossFit* fELossFit; // Energy loss fits
+ AliFMDCorrSecondaryMap* fSecondaryMap; // Secondary particle correction
+ AliFMDCorrDoubleHit* fDoubleHit; // Double hit correction (low flux)
+ AliFMDCorrVertexBias* fVertexBias; // Vertex bias correction
+ AliFMDCorrMergingEfficiency* fMergingEfficiency;
+ /*
+ * @}
+ */
+
+ ClassDef(AliForwardCorrectionManager,1) // Manager of corrections
+};
+//____________________________________________________________________
+inline const TAxis*
+AliForwardCorrectionManager::GetEtaAxis() const
+{
+ if (!fSecondaryMap) return 0;
+ return &(fSecondaryMap->GetEtaAxis());
+}
+//____________________________________________________________________
+inline const TAxis*
+AliForwardCorrectionManager::GetVertexAxis() const
+{
+ if (!fSecondaryMap) return 0;
+ return &(fSecondaryMap->GetVertexAxis());
+}
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
-#include "AliForwardMCCorrections.h"
+#include "AliForwardMCCorrectionsTask.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
-#include "AliFMDAnaParameters.h"
+// #include "AliFMDAnaParameters.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliFMDStripIndex.h"
#include <TH1.h>
+#include <TH3D.h>
#include <TDirectory.h>
#include <TTree.h>
namespace {
const char* GetEventName(Bool_t tr, Bool_t vtx)
{
- return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx", ""));
+ return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
}
const char* GetHitsName(UShort_t d, Char_t r)
{
}
//====================================================================
-AliForwardMCCorrections::AliForwardMCCorrections()
+AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
: AliAnalysisTaskSE(),
+ fHEvents(0),
fHEventsTr(0),
fHEventsTrVtx(0),
+ fHEventsVtx(0),
fHTriggers(0),
+ fPrimaryInnerAll(0),
+ fPrimaryOuterAll(0),
+ fPrimaryInnerTrVtx(0),
+ fPrimaryOuterTrVtx(0),
+ fHitsFMD1i(0),
+ fHitsFMD2i(0),
+ fHitsFMD2o(0),
+ fHitsFMD3i(0),
+ fHitsFMD3o(0),
+ fStripsFMD1i(0),
+ fStripsFMD2i(0),
+ fStripsFMD2o(0),
+ fStripsFMD3i(0),
+ fStripsFMD3o(0),
fVtxAxis(),
- fEtaAxis()
+ fEtaAxis(),
+ fList()
{
}
//____________________________________________________________________
-AliForwardMCCorrections::AliForwardMCCorrections(const char* name)
+AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
: AliAnalysisTaskSE(name),
+ fHEvents(0),
fHEventsTr(0),
- fHEventsTrVtx(0),
+ fHEventsTrVtx(0),
+ fHEventsVtx(0),
fHTriggers(0),
+ fPrimaryInnerAll(0),
+ fPrimaryOuterAll(0),
+ fPrimaryInnerTrVtx(0),
+ fPrimaryOuterTrVtx(0),
+ fHitsFMD1i(0),
+ fHitsFMD2i(0),
+ fHitsFMD2o(0),
+ fHitsFMD3i(0),
+ fHitsFMD3o(0),
+ fStripsFMD1i(0),
+ fStripsFMD2i(0),
+ fStripsFMD2o(0),
+ fStripsFMD3i(0),
+ fStripsFMD3o(0),
fVtxAxis(10,-10,10),
- fEtaAxis(200,-4,6)
+ fEtaAxis(200,-4,6),
+ fList()
{
DefineOutput(1, TList::Class());
}
//____________________________________________________________________
-AliForwardMCCorrections::AliForwardMCCorrections(const AliForwardMCCorrections& o)
+AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
: AliAnalysisTaskSE(o),
+ fHEvents(o.fHEvents),
fHEventsTr(o.fHEventsTr),
- fHEventsTrVtx(o.fHEventsTrVtx),
+ fHEventsTrVtx(o.fHEventsTrVtx),
+ fHEventsVtx(o.fHEventsVtx),
fHTriggers(o.fHTriggers),
- fVtxAxis(o.fVtxAxis),
- fEtaAxis(o.fEtaAxis)
+ fPrimaryInnerAll(o.fPrimaryInnerAll),
+ fPrimaryOuterAll(o.fPrimaryOuterAll),
+ fPrimaryInnerTrVtx(o.fPrimaryInnerTrVtx),
+ fPrimaryOuterTrVtx(o.fPrimaryOuterTrVtx),
+ fHitsFMD1i(o.fHitsFMD1i),
+ fHitsFMD2i(o.fHitsFMD2i),
+ fHitsFMD2o(o.fHitsFMD2o),
+ fHitsFMD3i(o.fHitsFMD3i),
+ fHitsFMD3o(o.fHitsFMD3o),
+ fStripsFMD1i(o.fStripsFMD1i),
+ fStripsFMD2i(o.fStripsFMD2i),
+ fStripsFMD2o(o.fStripsFMD2o),
+ fStripsFMD3i(o.fStripsFMD3i),
+ fStripsFMD3o(o.fStripsFMD3o),
+ fVtxAxis(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), o.fVtxAxis.GetXmax()),
+ fEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax()),
+ fList(o.fList)
{
}
//____________________________________________________________________
-AliForwardMCCorrections&
-AliForwardMCCorrections::operator=(const AliForwardMCCorrections& o)
+AliForwardMCCorrectionsTask&
+AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
{
fHEventsTr = o.fHEventsTr;
fHEventsTrVtx = o.fHEventsTrVtx;
fHTriggers = o.fHTriggers;
- fHData = o.fHData;
- fVtxAxis = o.fVtxAxis;
- fEtaAxis = o.fEtaAxis;
+ SetVertexAxis(o.fVtxAxis);
+ SetEtaAxis(o.fEtaAxis);
return *this;
}
//____________________________________________________________________
void
-AliForwardMCCorrections::Init()
+AliForwardMCCorrectionsTask::Init()
{
}
//____________________________________________________________________
void
-AliForwardMCCorrections::SetVertexAxis(Int_t nBin, Double_t min, Double_t max)
+AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
+ Double_t max)
{
if (max < min) max = -min;
if (min < max) {
}
//____________________________________________________________________
void
-AliForwardMCCorrections::SetVertexAxis(const TAxis& axis)
+AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
{
SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
}
//____________________________________________________________________
void
-AliForwardMCCorrections::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
+AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
{
if (max < min) max = -min;
if (min < max) {
}
//____________________________________________________________________
void
-AliForwardMCCorrections::SetVertexAxis(const TAxis& axis)
+AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
{
SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
}
-//____________________________________________________________________
-void
-AliForwardMCCorrections::InitializeSubs()
-{
-}
-
//____________________________________________________________________
TH3D*
-AliForwardMCCorrections::Make3D(const char* name, const char* title,
+AliForwardMCCorrectionsTask::Make3D(const char* name, const char* title,
Int_t nPhi) const
{
TH3D* ret = new TH3D(name, title,
}
//____________________________________________________________________
TH1D*
-AliForwardMCCorrections::Make1D(const char* name, const char* title) const
+AliForwardMCCorrectionsTask::Make1D(const char* name, const char* title) const
{
TH1D* ret = new TH1D(name, title,
fEtaAxis.GetNbins(),
//____________________________________________________________________
void
-AliForwardMCCorrections::UserCreateOutputObjects()
+AliForwardMCCorrectionsTask::UserCreateOutputObjects()
{
fList = new TList;
fList->SetName(GetName());
- fHEvents = new TH1I(GetEventsName(false,false),
+ fHEvents = new TH1I(GetEventName(false,false),
"Number of all events",
- fAxis.GetNbins(),
- fAxis.GetXmin(),
- fAxis.GetXmax());
+ fVtxAxis.GetNbins(),
+ fVtxAxis.GetXmin(),
+ fVtxAxis.GetXmax());
fHEvents->SetXTitle("v_{z} [cm]");
fHEvents->SetYTitle("# of events");
fHEvents->SetFillColor(kBlue+1);
fHEvents->SetDirectory(0);
fList->Add(fHEvents);
- fHEventsTr = new TH1I(GetEventsName(true, false),
+ fHEventsTr = new TH1I(GetEventName(true, false),
+ "Number of triggered events",
fVtxAxis.GetNbins(),
fVtxAxis.GetXmin(),
fVtxAxis.GetXmax());
fHEventsTr->SetDirectory(0);
fList->Add(fHEventsTr);
- fHEventsTrVtx = new TH1I(GetEventsName(true,true),
+ fHEventsTrVtx = new TH1I(GetEventName(true,true),
"Number of events w/trigger and vertex",
fVtxAxis.GetNbins(),
fVtxAxis.GetXmin(),
fHEventsTrVtx->SetFillStyle(3001);
fHEventsTrVtx->SetDirectory(0);
fList->Add(fHEventsTrVtx);
-
- fHEventsVtx = new TH1I(GetEventsName(false,true),
+
+ fHEventsVtx = new TH1I(GetEventName(false,true),
"Number of events w/vertex",
fVtxAxis.GetNbins(),
fVtxAxis.GetXmin(),
}
//____________________________________________________________________
void
-AliForwardMCCorrections::UserExec(Option_t*)
+AliForwardMCCorrectionsTask::UserExec(Option_t*)
{
// Get the input data - MC event
AliMCEvent* mcEvent = MCEvent();
- if (mcevent) {
+ if (mcEvent) {
AliWarning("No MC event found");
return;
}
// Get the event generate header
AliHeader* mcHeader = mcEvent->Header();
- AliGenEventHeader* genHeader = mcHeader->GetCenEventHeader();
+ AliGenEventHeader* genHeader = mcHeader->GenEventHeader();
// Get the generator vertex
TArrayF mcVertex;
- genHeader->PrimaryVertex();
+ genHeader->PrimaryVertex(mcVertex);
Double_t mcVtxZ = mcVertex.At(2);
// Check the MC vertex is in range
- Int_t mcVtxBin = fVtxAxis.FindBin(vZ);
- if (vZ < 1 || vZ > fVtxAxis.GetNbins()) {
+ Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ);
+ if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) {
#ifdef VERBOSE
AliWarning(Form("MC v_z=%f is out of rante [%,%f]",
- vZ, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
+ mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
#endif
return;
}
+ UInt_t triggers;
+ Bool_t gotTrigggers;
+ Bool_t gotInel;
+ Double_t vZ;
+ Bool_t gotVertex;
+#if 0
+ // Use event inspector instead
// Get the triggers
UInt_t triggers = 0;
Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers);
// Get the ESD vertex
Double_t vZ = -1000000;
Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ);
-
+#endif
+
// Fill the event count histograms
if (gotInel) fHEventsTr->Fill(mcVtxZ);
// Fill (eta,phi) of the particle into histograsm for b
Double_t eta = particle->Eta();
- Double_t phi = partcile->Phi();
+ Double_t phi = particle->Phi();
if (isCharged && isPrimary)
FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi);
// For the rest, ignore non-collisions, and events out of vtx range
- if (!gotInel || gotVtx) continue;
+ if (!gotInel || gotVertex) continue;
Int_t nTrRef = particle->GetNumberOfTrackReferences();
for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
AliTrackReference* trackRef = particle->GetTrackReference(iTrRef);
// Check existence
- if (!ret) continue;
+ if (!trackRef) continue;
// Check that we hit an FMD element
- if (ref->DetectorId() != AliTrackReference::kFMD)
+ if (trackRef->DetectorId() != AliTrackReference::kFMD)
continue;
// Get the detector coordinates
UShort_t d, s, t;
Char_t r;
- AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
+ AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
// Check if mother (?) is charged and that we haven't dealt with it
// already
// Increase counter of hits per strip
hitsByStrip(d,r,s,t) += 1;
- Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
- Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
+ // Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
+ // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
// Fill strip information into histograms
- FillStrip(mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
+ FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
// Set the last processed track number - marking it as done for
// this strip
//____________________________________________________________________
void
-AliForwardMCCorrections::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
+AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
Double_t vz, Double_t eta, Double_t phi)
{
if (gotInel && gotVtx) {
//____________________________________________________________________
void
-AliForwardMCCorrections::FillStrip(UShort_t d, Char_t r,
+AliForwardMCCorrectionsTask::FillStrip(UShort_t d, Char_t r,
Double_t vz, Double_t eta, Double_t phi,
Bool_t first)
{
// Number of hit strips per eta bin
- TH1D* hits = 0;
+ TH1D* strips = 0;
// All hits per ring, vertex in eta,phi bins. This takes the place
// of hHits_FMD<d><r>_vtx<v> and DoubleHits_FMD<d><r> (the later
// being just the 2D projection over the X bins)
- TH3D* strips = 0;
+ TH3D* hits = 0;
switch (d) {
case 1:
hits = fHitsFMD1i;
//____________________________________________________________________
TH2D*
-AliForwardMCCorrections::GetVertexProj(Int_t v, TH3D* src) const
+AliForwardMCCorrectionsTask::GetVertexProj(Int_t v, TH3D* src) const
{
Int_t nX = src->GetNbinsX();
if (v > nX || v < 1) return 0;
//____________________________________________________________________
void
-AliForwardMCCorrections::Terminate(Option_t*)
+AliForwardMCCorrectionsTask::Terminate(Option_t*)
{
TList* list = dynamic_cast<TList*>(GetOutputData(1));
if (!list) {
}
TH1I* eventsAll =
- static_cast<TH1I*>(list->FindObject(GetEventsName(false,false)));
+ static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
TH1I* eventsTr =
- static_cast<TH1I*>(list->FindObject(GetEventsName(true,false)));
+ static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
TH1I* eventsVtx =
- static_cast<TH1I*>(list->FindObject(GetEventsName(false,true)));
+ static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
TH1I* eventsTrVtx =
- static_cast<TH1I*>(list->FindObject(GetEventsName(true,true)));
+ static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) {
AliError(Form("Could not find all event histograms in %s "
"(a=%p,t=%p,v=%p,tv=%p)", list->GetName(),
// Calculate the over-all event selection efficiency
TH1D* selEff = new TH1D("selEff", "Event selection efficiency",
- fVtxArray.GetNbins(),
- fVtxArray.GetXmin(),
- fVtxArray.GetXmax());
+ fVtxAxis.GetNbins(),
+ fVtxAxis.GetXmin(),
+ fVtxAxis.GetXmax());
selEff->Sumw2();
selEff->SetDirectory(0);
selEff->SetFillColor(kMagenta+1);
for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) {
// Make a sub-list in the output
TList* vl = new TList;
- vl->SetName("vtx%02d", v);
+ vl->SetName(Form("vtx%02d", v));
list->Add(vl);
// Get event counts
- Int_t nEventsAll = eventsAll->GetEntries(v);
- Int_t nEventsTr = eventsTr->GetEntries(v);
- Int_t nEventsVtx = eventsVtx->GetEntries(v);
- Int_t nEventsTrVtx = eventsTrVtx->GetEntries(v);
+ Int_t nEventsAll = eventsAll->GetBinContent(v);
+ Int_t nEventsTr = eventsTr->GetBinContent(v);
+ Int_t nEventsVtx = eventsVtx->GetBinContent(v);
+ Int_t nEventsTrVtx = eventsTrVtx->GetBinContent(v);
// Project event histograms, set names, and store
TH2D* primIAllV = GetVertexProj(v, primIAll);
vl->Add(primITrVtxV);
vl->Add(primOTrVtxV);
- primIAllV->Scale(1. / nEventAll);
- primOAllV->Scale(1. / nEventAll);
- primITrVtxV->Scale(1. / nEventTr);
- primOTrVtxV->Scale(1. / nEventTr);
+ primIAllV->Scale(1. / nEventsAll);
+ primOAllV->Scale(1. / nEventsAll);
+ primITrVtxV->Scale(1. / nEventsTr);
+ primOTrVtxV->Scale(1. / nEventsTr);
// Calculate the vertex bias on the d^2N/detadphi
- TH2D* selBiasI = static_cast<TH2D*>(primITrVtxV->Clone("selBiasI%02d",v));
- TH2D* selBiasO = static_cast<TH2D*>(primOTrVtxV->Clone("selBiasO%02d",v));
+ TH2D* selBiasI =
+ static_cast<TH2D*>(primITrVtxV->Clone(Form("selBiasI%02d",v)));
+ TH2D* selBiasO =
+ static_cast<TH2D*>(primOTrVtxV->Clone(Form("selBiasO%02d",v)));
selBiasI->SetTitle(Form("Event selection bias in vertex bin %d", v));
selBiasO->SetTitle(Form("Event selection bias in vertex bin %d", v));
selBiasI->Divide(primIAllV);
// Do the double hit correction (only done once per ring in
// the vertex loop)
- TH1D* strips = 0;
+ TH1D* hStrips = 0;
switch (d) {
- case 1: strips = strips1i; break;
- case 2: strips = (q == 0 ? strips2i : strips2o); break;
- case 3: strips = (q == 0 ? strips3i : strips3o); break;
+ case 1: hStrips = strips1i; break;
+ case 2: hStrips = (q == 0 ? strips2i : strips2o); break;
+ case 3: hStrips = (q == 0 ? strips3i : strips3o); break;
}
TH2D* hits2D = GetVertexProj(v, hits3D);
doubleHit->SetFillColor(kGreen+1);
doubleHit->SetFillStyle(3001);
doubleHit->Sumw2();
- doubleHit->Divide(strips);
+ doubleHit->Divide(hStrips);
list->Add(doubleHit);
}
}
//____________________________________________________________________
void
-AliForwardMCCorrections::Print(Option_t*) const
+AliForwardMCCorrectionsTask::Print(Option_t*) const
{
}
class AliFMDAnaParameters;
class AliESDEvent;
class TH2D;
+class TH3D;
class TList;
class TTree;
* @ingroup pwg2_forward_analysis
*
*/
-class AliForwardMCCorrections : public AliAnalysisTaskSE
+class AliForwardMCCorrectionsTask : public AliAnalysisTaskSE
{
public:
/**
*
* @param name Name of task
*/
- AliForwardMCCorrections(const char* name);
+ AliForwardMCCorrectionsTask(const char* name);
/**
* Constructor
*/
- AliForwardMCCorrections();
+ AliForwardMCCorrectionsTask();
/**
* Copy constructor
*
* @param o Object to copy from
*/
- AliForwardMCCorrections(const AliForwardMCCorrections& o);
+ AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o);
/**
* Assignment operator
*
*
* @return Reference to this object
*/
- AliForwardMCCorrections& operator=(const AliForwardMCCorrections& o);
+ AliForwardMCCorrectionsTask& operator=(const AliForwardMCCorrectionsTask& o);
/**
* @{
* @name Interface methods
void SetEtaAxis(const TAxis& axis);
protected:
TH2D* GetVertexProj(Int_t v, TH3D* src) const;
-
+ TH3D* Make3D(const char* name, const char* title, Int_t nPhi) const;
+ TH1D* Make1D(const char* name, const char* title) const;
+ void FillPrimary(Bool_t gotInel, Bool_t gotVtx,
+ Double_t vz, Double_t eta, Double_t phi);
+ void FillStrip(UShort_t d, Char_t r,
+ Double_t vz, Double_t eta, Double_t phi,
+ Bool_t first);
TH1I* fHEvents; // All Events
TH1I* fHEventsTr; // Histogram of events w/trigger
TH1I* fHEventsTrVtx; // Events w/trigger and vertex
TList* fList; // Output list
- ClassDef(AliForwardMCCorrections,1) // Forward corrections class
+ ClassDef(AliForwardMCCorrectionsTask,1) // Forward corrections class
};
#endif
-#include "AliForwardMultiplicity.h"
+#include "AliForwardMultiplicityTask.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
-#include "AliFMDAnaParameters.h"
+// #include "AliFMDAnaParameters.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
+#include "AliForwardCorrectionManager.h"
+#include "AliAnalysisManager.h"
#include <TH1.h>
#include <TDirectory.h>
#include <TTree.h>
+#include <TROOT.h>
+#include <iostream>
+#include <iomanip>
//====================================================================
-AliForwardMultiplicity::AliForwardMultiplicity()
+AliForwardMultiplicityTask::AliForwardMultiplicityTask()
: AliAnalysisTaskSE(),
+ fEnableLowFlux(true),
fHData(0),
fFirstEvent(true),
fESDFMD(),
}
//____________________________________________________________________
-AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
+AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
: AliAnalysisTaskSE(name),
+ fEnableLowFlux(true),
fHData(0),
fFirstEvent(true),
fESDFMD(),
}
//____________________________________________________________________
-AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
+AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
: AliAnalysisTaskSE(o),
+ fEnableLowFlux(o.fEnableLowFlux),
fHData(o.fHData),
fFirstEvent(true),
fESDFMD(o.fESDFMD),
fHistCollector(o.fHistCollector),
fList(o.fList)
{
+ DefineOutput(1, TList::Class());
}
//____________________________________________________________________
-AliForwardMultiplicity&
-AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
+AliForwardMultiplicityTask&
+AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
{
AliAnalysisTaskSE::operator=(o);
+ fEnableLowFlux = o.fEnableLowFlux;
fHData = o.fHData;
fFirstEvent = o.fFirstEvent;
fEventInspector = o.fEventInspector;
//____________________________________________________________________
void
-AliForwardMultiplicity::SetDebug(Int_t dbg)
+AliForwardMultiplicityTask::SetDebug(Int_t dbg)
{
fEventInspector.SetDebug(dbg);
fEnergyFitter.SetDebug(dbg);
//____________________________________________________________________
void
-AliForwardMultiplicity::Init()
+AliForwardMultiplicityTask::Init()
{
fFirstEvent = true;
}
//____________________________________________________________________
void
-AliForwardMultiplicity::InitializeSubs()
+AliForwardMultiplicityTask::InitializeSubs()
{
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
- pars->Init(kTRUE);
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // pars->Init(kTRUE);
+
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ fcm.Init(fEventInspector.GetCollisionSystem(),
+ fEventInspector.GetEnergy(),
+ fEventInspector.GetField());
+ // Check that we have the energy loss fits, needed by
+ // AliFMDSharingFilter
+ // AliFMDDensityCalculator
+ if (!fcm.GetELossFit()) {
+ AliFatal(Form("No energy loss fits"));
+ return;
+ }
+ // Check that we have the double hit correction - (optionally) used by
+ // AliFMDDensityCalculator
+ if (!fcm.GetDoubleHit()) {
+ AliWarning("No double hit corrections");
+ }
+ // Check that we have the secondary maps, needed by
+ // AliFMDCorrections
+ // AliFMDHistCollector
+ if (!fcm.GetSecondaryMap()) {
+ AliFatal("No secondary corrections");
+ return;
+ }
+ // Check that we have the vertex bias correction, needed by
+ // AliFMDCorrections
+ if (!fcm.GetVertexBias()) {
+ AliFatal("No event vertex bias corrections");
+ return;
+ }
+ // Check that we have the merging efficiencies, optionally used by
+ // AliFMDCorrections
+ if (!fcm.GetMergingEfficiency()) {
+ AliWarning("No merging efficiencies");
+ }
+
+
+ const TAxis* pe = fcm.GetEtaAxis();
+ const TAxis* pv = fcm.GetVertexAxis();
+ if (!pe) AliFatal("No eta axis defined");
+ if (!pv) AliFatal("No vertex axis defined");
- TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
- TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
-
- fHistos.Init(e);
- fAODFMD.Init(e);
+ // TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
+ // TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
+ // pv->Dump();
+
+ fHistos.Init(*pe);
+ fAODFMD.Init(*pe);
fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
fHData->SetStats(0);
fHData->SetDirectory(0);
fList->Add(fHData);
- fEnergyFitter.Init(e);
- fEventInspector.Init(v);
- fHistCollector.Init(v);
+ fEnergyFitter.Init(*pe);
+ fEventInspector.Init(*pv);
+ fHistCollector.Init(*pv);
+
+ this->Print();
}
//____________________________________________________________________
void
-AliForwardMultiplicity::UserCreateOutputObjects()
+AliForwardMultiplicityTask::UserCreateOutputObjects()
{
fList = new TList;
AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
- AliAODHandler* ah = dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
+ AliAODHandler* ah =
+ dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
if (!ah) AliFatal("No AOD output handler set in analysis manager");
fSharingFilter.DefineOutput(fList);
fDensityCalculator.DefineOutput(fList);
fCorrections.DefineOutput(fList);
+
+ PostData(1, fList);
}
//____________________________________________________________________
void
-AliForwardMultiplicity::UserExec(Option_t*)
+AliForwardMultiplicityTask::UserExec(Option_t*)
{
+ // static Int_t cnt = 0;
+ // cnt++;
// Get the input data
AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
if (!esd) {
AliWarning("No ESD event found for input event");
return;
// On the first event, initialize the parameters
if (fFirstEvent && esd->GetESDRun()) {
- AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ fEventInspector.ReadRunDetails(esd);
+
AliInfo(Form("Initializing with parameters from the ESD:\n"
" AliESDEvent::GetBeamEnergy() ->%f\n"
" AliESDEvent::GetBeamType() ->%s\n"
esd->GetCurrentL3(),
esd->GetMagneticField(),
esd->GetRunNumber()));
- pars->SetParametersFromESD(esd);
- pars->PrintStatus();
+
+
+
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ // pars->SetParametersFromESD(esd);
+ // pars->PrintStatus();
fFirstEvent = false;
InitializeSubs();
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
- Int_t ivz = -1;
+ UShort_t ivz = 0;
Double_t vz = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
if (found & AliFMDEventInspector::kNoEvent) return;
if (found & AliFMDEventInspector::kBadVertex) return;
+ // We we do not want to use low flux specific code, we disable it here.
+ if (!fEnableLowFlux) lowFlux = false;
+
// Get FMD data
AliESDFMD* esdFMD = esd->GetFMDData();
// Apply the sharing filter (or hit merging or clustering if you like)
//____________________________________________________________________
void
-AliForwardMultiplicity::Terminate(Option_t*)
+AliForwardMultiplicityTask::Terminate(Option_t*)
{
TList* list = dynamic_cast<TList*>(GetOutputData(1));
if (!list) {
// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
- TH1D* norm = hData->ProjectionX("dNdeta", 0, 1, "");
+ TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
dNdeta->Divide(norm);
dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
"width");
list->Add(dNdeta);
+ list->Add(norm);
fEnergyFitter.Fit(list);
fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
//____________________________________________________________________
void
-AliForwardMultiplicity::MarkEventForStore() const
+AliForwardMultiplicityTask::MarkEventForStore() const
{
// Make sure the AOD tree is filled
AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
//____________________________________________________________________
void
-AliForwardMultiplicity::Print(Option_t*) const
+AliForwardMultiplicityTask::Print(Option_t* option) const
{
+ std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n"
+ << " Enable low flux code: " << (fEnableLowFlux ? "yes" : "no")
+ << std::endl;
+ gROOT->IncreaseDirLevel();
+ fEventInspector .Print(option);
+ fEnergyFitter .Print(option);
+ fSharingFilter .Print(option);
+ fDensityCalculator.Print(option);
+ fCorrections .Print(option);
+ fHistCollector .Print(option);
+ gROOT->DecreaseDirLevel();
}
//
* @ingroup pwg2_forward_analysis
*
*/
-class AliForwardMultiplicity : public AliAnalysisTaskSE
+class AliForwardMultiplicityTask : public AliAnalysisTaskSE
{
public:
/**
*
* @param name Name of task
*/
- AliForwardMultiplicity(const char* name);
+ AliForwardMultiplicityTask(const char* name);
/**
* Constructor
*/
- AliForwardMultiplicity();
+ AliForwardMultiplicityTask();
/**
* Copy constructor
*
* @param o Object to copy from
*/
- AliForwardMultiplicity(const AliForwardMultiplicity& o);
+ AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o);
/**
* Assignment operator
*
*
* @return Reference to this object
*/
- AliForwardMultiplicity& operator=(const AliForwardMultiplicity& o);
+ AliForwardMultiplicityTask& operator=(const AliForwardMultiplicityTask& o);
/**
* @{
* @name Interface methods
/**
* @}
*/
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
void Print(Option_t* option="") const;
-
+ /**
+ * Whether to enable low-flux code
+ *
+ * @param use IF true, enable low-flux code
+ */
+ void SetEnableLowFlux(Bool_t use=true) { fEnableLowFlux = use; }
/**
* @{
* @name Access to sub-algorithms
*/
virtual void MarkEventForStore() const;
+ Bool_t fEnableLowFlux;// Whether to use low-flux specific code
TH2D* fHData; // Summed 1/Nd^2N_{ch}/dphideta
Bool_t fFirstEvent; // Whether the event is the first seen
AliESDFMD fESDFMD; // Sharing corrected ESD object
TList* fList; // Output list
- ClassDef(AliForwardMultiplicity,1) // Forward multiplicity class
+ ClassDef(AliForwardMultiplicityTask,1) // Forward multiplicity class
};
#endif
#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;
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);
+ }
}
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-2];
- 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,
dist->GetXaxis()->SetRangeUser(0, fMaxRange);
// Define the function to fit
- TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
+ TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
// Set initial guesses, parameter names, and limits
landau1->SetParameters(1,0.5,0.07,0.1,sigman);
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,kN+n);
- landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
- landaun->SetLineColor(((n-2) % 10)+2); // start at red
- landaun->SetLineWidth(1);
- landaun->SetNpx(500);
- landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
-
- // Set the initial parameters from the seed fit
- landaun->SetParameter(kC, r->Parameter(kC)); // Constant
- landaun->SetParameter(kDelta, r->Parameter(kDelta)); // Delta
- landaun->SetParameter(kXi, r->Parameter(kXi)); // xi
- landaun->SetParameter(kSigma, r->Parameter(kSigma)); // sigma
- landaun->SetParameter(kSigmaN, r->Parameter(kSigmaN)); // sigma_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);
- // Fix the number parameter
- landaun->FixParameter(kN, n); // 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(kA+i-2, n == 2 ? 0.05 : 0.000001);
landaun->SetParLimits(kA+i-2, 0,1);
- landaun->SetParName(kA+i-2, Form("a_{%d}", i));
}
// Do the fit
class AliForwardUtil : public TObject
{
public:
+ //==================================================================
+ /**
+ * @{
+ * @nane Collision/run parameters
+ */
+ /**
+ * Defined collision types
+ */
+ enum ECollisionSystem {
+ kUnknown,
+ kPP,
+ kPbPb
+ };
+ //__________________________________________________________________
+ /**
+ * Parse a collision system spec given in a string. Known values are
+ *
+ * - "pp", "p-p" which returns kPP
+ * - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
+ * - Everything else gives kUnknown
+ *
+ * @param sys Collision system spec
+ *
+ * @return Collision system id
+ */
+ static UShort_t ParseCollisionSystem(const char* sys);
+ /**
+ * Get a string representation of the collision system
+ *
+ * @param sys Collision system
+ * - kPP -> "pp"
+ * - kPbPb -> "PbPb"
+ * - anything else gives "unknown"
+ *
+ * @return String representation of the collision system
+ */
+ static const char* CollisionSystemString(UShort_t sys);
+ //__________________________________________________________________
+ /**
+ * Parse the center of mass energy given as a float and return known
+ * values as a unsigned integer
+ *
+ * @param sys Collision system (needed for AA)
+ * @param cms Center of mass energy * total charge
+ *
+ * @return Center of mass energy per nucleon
+ */
+ static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
+ /**
+ * Get a string representation of the center of mass energy per nuclean
+ *
+ * @param sys Collision system
+ * @param sNN Center of mass energy per nucleon
+ *
+ * @return String representation of the center of mass energy per nuclean
+ */
+ static const char* CenterOfMassEnergyString(UShort_t cms);
+ //__________________________________________________________________
+ /**
+ * Parse the magnetic field (in kG) as given by a floating point number
+ *
+ * @param field Magnetic field in kG
+ *
+ * @return Short integer value of magnetic field in kG
+ */
+ static Short_t ParseMagneticField(Float_t field);
+ /**
+ * Get a string representation of the magnetic field
+ *
+ * @param field Magnetic field in kG
+ *
+ * @return String representation of the magnetic field
+ */
+ static const char* MagneticFieldString(Short_t field);
+ /* @} */
+
+ /**
+ * @{
+ * @name Energy stragling functions
+ */
//__________________________________________________________________
/**
* Number of steps to do in the Landau, Gaussiam convolution
*/
static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi,
Double_t sigma, Double_t sigma_n);
-
+
+ //------------------------------------------------------------------
+ /**
+ * Evaluate
+ * @f[
+ * f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
+ * @f]
+ * corresponding to @f$ i@f$ particles i.e., with the substitutions
+ * @f[
+ * \Delta \rightarrow \Delta_i = i(\Delta + \xi\log(i))\\
+ * \xi \rightarrow \xi_i = i \xi\\
+ * \sigma \rightarrow \sigma_i = \sqrt{i}\sigma\\
+ * \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
+ * @f]
+ *
+ * @param x Where to evaluate
+ * @param delta @f$ \Delta@f$
+ * @param xi @f$ \xi@f$
+ * @param sigma @f$ \sigma@f$
+ * @param sigma_n @f$ \sigma_n@f$
+ * @param i @f$ i@f$
+ *
+ * @return @f$ f_i@f$ evaluated
+ */
+ static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n, Int_t i);
+
+ //------------------------------------------------------------------
+ /**
+ * Numerically evaluate
+ * @f[
+ * \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
+ * @f]
+ * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter. The mapping
+ * of the parameters is given by
+ *
+ * - 0: @f$\Delta@f$
+ * - 1: @f$\xi@f$
+ * - 2: @f$\sigma@f$
+ * - 3: @f$\sigma_n@f$
+ *
+ * This is the partial derivative with respect to the parameter of
+ * the response function corresponding to @f$ i@f$ particles i.e.,
+ * with the substitutions
+ * @f[
+ * \Delta \rightarrow \Delta_i = i(\Delta + \xi\log(i))\\
+ * \xi \rightarrow \xi_i = i \xi\\
+ * \sigma \rightarrow \sigma_i = \sqrt{i}\sigma\\
+ * \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
+ * @f]
+ *
+ * @param x Where to evaluate
+ * @param ipar Parameter number
+ * @param dp @f$ \esilon\delta p_i@f$ for some value of @f$\epsilon@f$
+ * @param delta @f$ \Delta@f$
+ * @param xi @f$ \xi@f$
+ * @param sigma @f$ \sigma@f$
+ * @param sigma_n @f$ \sigma_n@f$
+ * @param i @f$ i@f$
+ *
+ * @return @f$ f_i@f$ evaluated
+ */
+ static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
+ Double_t delta, Double_t xi,
+ Double_t sigma, Double_t sigma_n, Int_t i);
+
//------------------------------------------------------------------
/**
* Evaluate
* @f[
- f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f(x;\Delta_i,\xi_i,\sigma'_i)
- @f]
+ * f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
+ * @f]
*
* where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
* Landau with a Gaussian (see LandauGaus). Note that
static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
Double_t sigma, Double_t sigma_n, Int_t n,
Double_t* a);
+ /**
+ * Generate a TF1 object of @f$ f_I@f$
+ *
+ * @param c Constant
+ * @param delta @f$ \Delta@f$
+ * @param xi @f$ \xi_1@f$
+ * @param sigma @f$ \sigma_1@f$
+ * @param sigma_n @f$ \sigma_n@f$
+ * @param i @f$ i@f$ - the number of particles
+ * @param xmin Least value of range
+ * @param xmax Largest value of range
+ *
+ * @return Newly allocated TF1 object
+ */
+ static TF1* 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);
+ /**
+ * Generate a TF1 object of @f$ f_N@f$
+ *
+ * @param c Constant
+ * @param delta @f$ \Delta@f$
+ * @param xi @f$ \xi_1@f$
+ * @param sigma @f$ \sigma_1@f$
+ * @param sigma_n @f$ \sigma_n@f$
+ * @param n @f$ N@f$ - how many particles to sum to
+ * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
+ * @f$ i > 1@f$
+ * @param xmin Least value of range
+ * @param xmax Largest value of range
+ *
+ * @return Newly allocated TF1 object
+ */
+ static TF1* 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);
+
//__________________________________________________________________
/**
* Structure to do fits to the energy loss spectrum
TObjArray fFitResults; // Array of fit results
TObjArray fFunctions; // Array of functions
};
+ /* @} */
- //__________________________________________________________________
+ //==================================================================
+ /**
+ * @{
+ * @name Convenience containers
+ */
/**
* Structure to hold histograms
*
ClassDef(RingHistos,1)
};
+ /* @} */
};
Bool_t hasOther = (other && other->GetListOfGraphs() &&
other->GetListOfGraphs()->GetEntries() > 0);
Bool_t hasHhd = (hhd && hhdsym);
- if (!hasOther || !hasHhd) return 0;
+ if (!hasOther && !hasHhd) return 0;
THStack* ratios = new THStack("ratios", "Ratios");
if (hasOther) {
Double_t vzMax=10,
Int_t rebin=5,
const char* title="",
- bool hhd=false,
+ bool hhd=true,
bool comp=true)
{
- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C");
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/DrawRes.C","g");
Int_t trgMask;
trgMask = 1;
trgs.Append("INEL");
}
-
+ TString tit(title);
+ tit.ReplaceAll("@", " ");
+
printf("--------------------------------------\n"
"Settings for this:\n"
" Input AOD: %s\n"
" HHD comp.: %s\n"
" Other comp.: %s\n"
"--------------------------------------\n",
- file, vzMin, vzMax, rebin, trgMask, trgs.Data(), energy, title,
+ file, vzMin, vzMax, rebin, trgMask, trgs.Data(), energy, tit.Data(),
hhd ? "yes" : "no", comp ? "yes" : "no");
DrawRes dr;
TStopwatch t;
t.Start();
- dr.Run(file, vzMin, vzMax, rebin, trgMask, energy, title, hhd, comp);
+ dr.Run(file, vzMin, vzMax, rebin, trgMask, energy, tit.Data(), hhd, comp);
t.Stop();
t.Print();
}
#!/bin/bash
+ana=$ALICE_ROOT/PWG2/FORWARD/analysis2
nev=10000
noanal=0
nodraw=0
if test $batch -gt 0 ; then
opts="-l -b -q -x"
redir="2>&1 | tee ${base}.log"
+ echo "redir=$redir"
fi
if test $noanal -lt 1 ; then
rm -f AnalysisResult.root AliAODs.root
if test $gdb -gt 0 ; then
export PROOF_WRAPPERCMD="gdb -batch -x $ALICE_ROOT/PWG2/FORWARD/analysis2/gdb_cmds --args"
fi
- aliroot $opts $ALICE_ROOT/PWG2/FORWARD/analysis2/Pass1.C\(\".\",$nev,$af\) \
- $redir
- rm -f event_stat.root EventStat_temp.root outputs_valid
+ echo "Running aliroot ${opts} ${ana}/Pass1.C\(\".\",$nev,$af\) $redir"
+ if test $batch -gt 0 ; then
+ aliroot $opts ${ana}/Pass1.C\(\".\",$nev,$af\) 2>&1 | tee ${base}.log
+ else
+ aliroot $opts ${ana}/Pass1.C\(\".\",$nev,$af\)
+ fi
+ rm -f event_stat.root \
+ EventStat_temp.root \
+ outputs_valid \
+ `printf %09d.stat $nev`
if test ! -f AnalysisResults.root || test ! -f AliAODs.root ; then
echo "Analysis failed"
exit 1
if test $nodraw -lt 1 ; then
rm -f result.root
if test "x$tit" = "x" ; then
- tit="$nev\ events,\ v_{z}#in[$vzmin,$vzmax],\ $type"
- else
- tit=`echo $tit | tr ' ' '\ '`
+ tit="$nev events, v_{z}#in[$vzmin,$vzmax], $type"
fi
- aliroot ${opts} $ALICE_ROOT/PWG2/FORWARD/analysis2/Pass2.C\(\"AliAODs.root\",\"$type\",$cms,$vzmin,$vzmax,$rebin,\"$tit\",$hhd,$comp\)
+ tit=`echo $tit | tr ' ' '@'`
+ aliroot ${opts} ${ana}/Pass2.C\(\"AliAODs.root\",\"$type\",$cms,$vzmin,$vzmax,$rebin,\"$tit\",$hhd,$comp\)
fi
--- /dev/null
+/**
+ * Flags for the analysis
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+/**
+ * Run a pass on ESD data to produce the energ loss fits
+ *
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void RunELossFitter(const char* esddir,
+ Int_t nEvents=1000,
+ Bool_t mc=false,
+ Bool_t proof=false)
+{
+ // --- Libraries to load -------------------------------------------
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+ // --- Check for proof mode, and possibly upload pars --------------
+ if (proof)
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
+
+ // --- Our data chain ----------------------------------------------
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C");
+ TChain* chain = MakeESDChain(esddir);
+ // If 0 or less events is select, choose all
+ if (nEvents <= 0) nEvents = chain->GetEntries();
+ Info("RunELossFitter", "Will analyse %d events", nEvents);
+
+ // --- Creating the manager and handlers ---------------------------
+ AliAnalysisManager *mgr = new AliAnalysisManager("Analysis Train",
+ "FMD analysis train");
+
+ AliESDInputHandler *esdHandler = new AliESDInputHandler();
+ esdHandler->SetInactiveBranches("AliESDACORDE "
+ "AliRawDataErrorLogs "
+ "CaloClusters "
+ "Cascades "
+ "EMCALCells "
+ "EMCALTrigger "
+ "Kinks "
+ "Cascades "
+ "MuonTracks "
+ "TrdTracks "
+ "CaloClusters "
+ "HLTGlobalTrigger");
+ mgr->SetInputEventHandler(esdHandler);
+
+ // AOD output handler
+ AliAODHandler* aodHandler = new AliAODHandler();
+ mgr->SetOutputEventHandler(aodHandler);
+ aodHandler->SetOutputFileName("AliAODs.root");
+
+ // --- Add tasks ---------------------------------------------------
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskFMD.C");
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ AddTaskPhysicsSelection(mc, kTRUE, kTRUE);
+ AliFMDEnergyFitterTask* task = new AliFMDEnergyFitterTask("fmdEnergyFitter");
+ mgr->AddTask(task);
+
+ // --- Make the output container and connect it --------------------
+ TString outputfile = "energyFits.root";
+ AliAnalysisDataContainer* histOut =
+ mgr->CreateContainer("Forward", TList::Class(),
+ AliAnalysisManager::kOutputContainer,outputfile);
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, histOut);
+
+ // --- Set parameters on the algorithms ----------------------------
+ // Set the number of SPD tracklets for which we consider the event a
+ // low flux event
+ task->GetEventInspector().SetLowFluxCut(1000);
+ // Set the maximum error on v_z [cm]
+ task->GetEventInspector().SetMaxVzErr(0.2);
+ // Set the eta axis to use - note, this overrides whatever is used
+ // by the rest of the algorithms - but only for the energy fitter
+ // algorithm.
+ task->GetEnergyFitter().SetEtaAxis(200, -4, 6);
+ // Set maximum energy loss to consider
+ task->GetEnergyFitter().SetMaxE(10);
+ // Set number of energy loss bins
+ task->GetEnergyFitter().SetNEbins(300);
+ // Set whether to use increasing bin sizes
+ task->GetEnergyFitter().SetUseIncreasingBins(true);
+ // Set whether to do fit the energy distributions
+ task->GetEnergyFitter().SetDoFits(kTRUE);
+ // Set whether to make the correction object
+ task->GetEnergyFitter().SetDoMakeObject(kTRUE);
+ // Set the low cut used for energy
+ task->GetEnergyFitter().SetLowCut(0.4);
+ // Set the number of bins to subtract from maximum of distributions
+ // to get the lower bound of the fit range
+ task->GetEnergyFitter().SetFitRangeBinWidth(4);
+ // Set the maximum number of landaus to try to fit (max 5)
+ task->GetEnergyFitter().SetNParticles(5);
+ // Set the minimum number of entries in the distribution before
+ // trying to fit to the data
+ task->GetEnergyFitter().SetMinEntries(1000);
+ // --- Set limits on fits the energy -------------------------------
+ // Maximum relative error on parameters
+ AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
+ // Least weight to use
+ AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
+ // Maximum value of reduced chi^2
+ AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 5;
+
+ // --- Run the analysis --------------------------------------------
+ TStopwatch t;
+ if (!mgr->InitAnalysis()) {
+ Error("RunManager", "Failed to initialize analysis train!");
+ return;
+ }
+ // Some informative output
+ mgr->PrintStatus();
+ // mgr->SetDebugLevel(3);
+ if (mgr->GetDebugLevel() < 1 && !proof) mgr->SetUseProgressBar(kTRUE);
+
+ // Run the train
+ t.Start();
+ Printf("=== RUNNING ANALYSIS ==================================");
+ mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+ t.Stop();
+ t.Print();
+}
+//
+// EOF
+//
UShort_t flags=kFull)
{
// --- Libraries to load -------------------------------------------
- gSystem->Load("libVMC");
- gSystem->Load("libTree");
-
- gSystem->Load("libSTEERBase");
-
- gSystem->Load("libESD") ;
- gSystem->Load("libAOD") ;
- gSystem->Load("libANALYSIS");
- gSystem->Load("libANALYSISalice");
-
- gSystem->Load("libPhysics");
- gSystem->Load("libPWG0base");
- gSystem->Load("libPWG0dep");
- gSystem->Load("libPWG2forward");
- gSystem->Load("libPWG2forward2");
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
// --- Check for proof mode, and possibly upload pars --------------
- if (flags & kProof) {
- TProof::Open("workers=2");
- const char* pkgs[] = { "STEERBase", "ESD", "AOD", "ANALYSIS",
- "ANALYSISalice", "PWG2forward", "PWG2forward2", 0};
- const char** pkg = pkgs;
- while (*pkg) {
- gProof->UploadPackage(Form("${ALICE_ROOT}/%s.par",*pkg));
- gProof->EnablePackage(*pkg);
- pkg++;
- }
- }
+ if (flags & kProof)
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
// --- Our data chain ----------------------------------------------
- TChain* chain = new TChain("esdTree");
-
- // --- Get list of ESDs --------------------------------------------
- // Open source directory, and make sure we go back to were we were
- TString oldDir(gSystem->WorkingDirectory());
- TSystemDirectory d(esddir, esddir);
- TList* files = d.GetListOfFiles();
- gSystem->ChangeDirectory(oldDir);
-
- // Sort list of files and check if we should add it
- files->Sort();
- TIter next(files);
- TSystemFile* file = 0;
- while ((file = static_cast<TSystemFile*>(next()))) {
- if (file->IsDirectory()) continue;
- TString name(file->GetName());
- if (!name.EndsWith(".root")) continue;
- if (!name.Contains("AliESDs")) continue;
- TString esd(Form("%s/%s", file->GetTitle(), name.Data()));
- Info("RunManager", "Adding %s to chain", esd.Data());
- chain->Add(esd);
- }
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C");
+ TChain* chain = MakeESDChain(esddir);
// If 0 or less events is select, choose all
if (nEvents <= 0) nEvents = chain->GetEntries();
-
// --- Creating the manager and handlers ---------------------------
AliAnalysisManager *mgr = new AliAnalysisManager("Analysis Train",
gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskFMD.C");
gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
AliAnalysisTask* task = AddTaskFMD();
- mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
+ // mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
task = AddTaskPhysicsSelection((flags & kMC), kTRUE, kTRUE);
- mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
+ // mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
// --- Run the analysis --------------------------------------------
TStopwatch t;
// Run the train
t.Start();
- if (!(flags & kTerminate))
+ if (!(flags & kTerminate)) {
+ Printf("=== RUNNING ANALYSIS ==================================");
mgr->StartAnalysis((flags & kProof) ? "proof" : "local", chain, nEvents);
+ }
else {
- mgr->ImportWrappers();
+ Printf("=== RUNNING TERMINATE =================================");
+ // mgr->ImportWrappers(0);
mgr->Terminate();
}
t.Stop();
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff