*
* @ingroup pwglf_forward_dndeta
*/
+
AliAnalysisTask*
AddTaskForwardMultDists(const char* trig = "V0AND",
Double_t vzMin = -4,
Double_t vzMax = +4,
- UShort_t maxN = 150,
- UShort_t nDiv = 1,
Bool_t usePhiAcc = true,
Bool_t useAsymm = false)
{
task->SetIpZRange(vzMin, vzMax);
// Set the trigger mask to use (INEL,INEL>0,NSD)
task->SetTriggerMask(trig);
- // Set the maximum number of charged particles
- task->SetMaxN(maxN);
- // Set number of divisions per particle number
- task->SetNDivisions(nDiv);
// Set whether to use stored phi acceptance
task->SetUsePhiAcc(usePhiAcc);
// add the task
mgr->AddTask(task);
+ // Variable size axis objects
+ // for |eta|<0.5 from CMD
+ AliForwardMultDists::BinSpec b05(-0.5, 0.5, -0.5);
+ b05.Push(21, 1);
+ b05.Push(1, 3);
+ b05.Push(1, 5);
+ b05.Push(3, 5); // <-- Extra
+
+ // for |eta|<0.5 from ALICE
+ AliForwardMultDists::BinSpec a05(-0.5, 0.5, -0.5);
+ a05.Push(21, 1);
+ a05.Push(3, 2);
+
+ // For |eta|<1 from CMS
+ AliForwardMultDists::BinSpec b10(-1.0, 1.0, -0.5);
+ b10.Push(35, 1);
+ b10.Push(1, 3);
+ b10.Push(4, 5);
+ b10.Push(2, 5); // <-- Extra
+
+ // For |eta|<1 from ALICE
+ AliForwardMultDists::BinSpec a10(-1.0, 1.0, -0.5);
+ a10.Push(41, 1);
+ a10.Push(1, 2);
+
+ // For |eta|<1.3 from ALICE
+ AliForwardMultDists::BinSpec a13(-1.3, 1.3, -0.5);
+ a13.Push(41, 1);
+ a13.Push(7, 2);
+
+ // For |eta|<1.5 from CMS
+ AliForwardMultDists::BinSpec b15(-1.5, 1.5, -0.5);
+ b15.Push(46, 1);
+ b15.Push(1, 2);
+ b15.Push(1, 3);
+ b15.Push(4, 5);
+ b15.Push(4, 5); // <-- Extra
+
+ // for |eta|<2.0 from CMS
+ AliForwardMultDists::BinSpec b20(-2.0, 2.0, -0.5);
+ b20.Push(62, 1);
+ b20.Push(2, 2);
+ b20.Push(1, 4);
+ b20.Push(3, 10);
+ b20.Push(2, 10); // <-- Extra
+
+ // for |eta|<2 from CMS
+ AliForwardMultDists::BinSpec b24(-2.4, 2.4, -0.5);
+ b24.Push(57, 1);
+ b24.Push(1, 2);
+ b24.Push(1, 3);
+ b24.Push(3, 10);
+ b24.Push(4, 10); // <-- Extra
+
// Add bins
- Double_t limits[] = { 1, 1.5, 2, 3, 0 };
- Double_t* plimit = limits;
- while (*plimit) {
- Double_t eta = *plimit;
- task->AddBin(-eta, +eta);
+ AliForwardMultDists::BinSpec* bs[] = { &b05, &b10, &b15, &b20, &b24, 0 };
+ AliForwardMultDists::BinSpec** pb = bs;
+ while (*pb) {
+ AliForwardMultDists::BinSpec* b = *pb;
+ task->AddBin(*b);
if (useAsymm) {
- task->AddBin(-eta, 0);
- task->AddBin(0, +eta);
+ task->AddBin(b->fEtaMin, 0, b->Axis());
+ task->AddBin(0, b->fEtaMax, b->Axis());
}
- plimit++;
+ pb++;
}
// --- create containers for input/output --------------------------
#include <TVector2.h>
#include <THStack.h>
#include <TParameter.h>
+#include <TError.h>
#include <iostream>
#include <iomanip>
+//____________________________________________________________________
+AliForwardMultDists::BinSpec::BinSpec(Double_t etaMin,
+ Double_t etaMax,
+ Double_t nchLow)
+ : fEtaMin(etaMin), fEtaMax(etaMax), fLow(nchLow), fN(0), fD(0), fAxis(1,0,1)
+{
+}
+//____________________________________________________________________
+void
+AliForwardMultDists::BinSpec::Push(UShort_t n, Double_t d)
+{
+ Int_t l = fN.GetSize();
+ fN.Set(l+1); // Not terribly efficient, but that's ok because we do
+ fD.Set(l+1); // this infrequently
+ fN[l] = n;
+ fD[l] = d;
+}
+//____________________________________________________________________
+const TAxis&
+AliForwardMultDists::BinSpec::Axis() const
+{
+ if (fAxis.GetNbins() > 1) return fAxis;
+ if (fN.GetSize() <= 0) return fAxis;
+ if (fN.GetSize() == 1) {
+ // Exactly one spec,
+ fAxis.Set(fN[0], fLow, fN[0]*fD[0]+fLow);
+ }
+ else {
+ Int_t n = fN.GetSum();
+ Int_t l = fN.GetSize();
+ Int_t i = 0;
+ TArrayD b(n+1);
+ b[0] = fLow;
+ for (Int_t j = 0; j < l; j++) {
+ for (Int_t k = 0; k < fN[j]; k++) {
+ b[i+1] = b[i] + fD[j];
+ i++;
+ }
+ }
+ if (i != n) {
+ ::Warning("Axis", "Assigned bins, and number of bins not equal");
+ n = TMath::Min(i, n);
+ }
+ fAxis.Set(n, b.GetArray());
+ }
+ return fAxis;
+}
+
+
+
//____________________________________________________________________
AliForwardMultDists::AliForwardMultDists()
: AliAnalysisTaskSE(),
fForwardCache(0),
fCentralCache(0),
fMCCache(0),
- fMaxN(200),
- fNDivisions(1),
fUsePhiAcc(true)
{}
fForwardCache(0),
fCentralCache(0),
fMCCache(0),
- fMaxN(200),
- fNDivisions(1),
fUsePhiAcc(true)
{
/**
fForwardCache(o.fForwardCache),
fCentralCache(o.fCentralCache),
fMCCache(o.fMCCache),
- fMaxN(o.fMaxN),
- fNDivisions(o.fNDivisions),
fUsePhiAcc(o.fUsePhiAcc)
{}
//____________________________________________________________________
fForwardCache = o.fForwardCache;
fCentralCache = o.fCentralCache;
fMCCache = o.fMCCache;
- fMaxN = o.fMaxN;
- fNDivisions = o.fNDivisions;
fUsePhiAcc = o.fUsePhiAcc;
return *this;
}
-//____________________________________________________________________
-void
-AliForwardMultDists::SetMaxN(UShort_t maxN,
- UShort_t nDivisions)
-{
- const UShort_t one = 1;
- fMaxN = maxN;
- fNDivisions = TMath::Max(nDivisions, one);
-}
-
//____________________________________________________________________
void
AliForwardMultDists::UserCreateOutputObjects()
EtaBin* bin = 0;
TIter next(&fBins);
while ((bin = static_cast<EtaBin*>(next()))) {
- bin->Terminate(in, out, fMaxN);
+ bin->Terminate(in, out);
sta = oth;
if (bin->IsSymmetric()) sta = sym;
fList->Add(AliForwardUtil::MakeParameter("trigger", ULong_t(fTriggerMask)));
fList->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
fList->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
- fList->Add(AliForwardUtil::MakeParameter("maxN", UShort_t(fMaxN)));
fList->Add(AliForwardUtil::MakeParameter("count", UShort_t(1)));
}
TIter next(&fBins);
EtaBin* bin = 0;
while ((bin = static_cast<EtaBin*>(next()))) {
- bin->SetupForData(fList, hist, fMaxN, fNDivisions, useMC);
+ bin->SetupForData(fList, hist, useMC);
}
}
}
//____________________________________________________________________
void
-AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax)
+AliForwardMultDists::AddBin(const BinSpec& spec)
+{
+ AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
+}
+
+//____________________________________________________________________
+void
+AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
+ const TAxis& mAxis)
{
/**
* Add an @f$\eta@f$ bin
etaLow = -TMath::Abs(etaLow);
etaMax = +TMath::Abs(etaLow);
}
- EtaBin* bin = new EtaBin(etaLow, etaMax);
+ EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
+ // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
+ fBins.Add(bin);
+}
+//____________________________________________________________________
+void
+AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
+ UShort_t nMax, UShort_t nDiv)
+{
+ /**
+ * Add an @f$\eta@f$ bin
+ *
+ * @param etaLow Low cut on @f$\eta@f$
+ * @param etaMax High cut on @f$\eta@f$
+ */
+ // Symmetric bin
+ if (etaMax >= kInvalidEta) {
+ etaLow = -TMath::Abs(etaLow);
+ etaMax = +TMath::Abs(etaLow);
+ }
+ TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
+ EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
// AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
fBins.Add(bin);
}
<< " Trigger mask: "
<< AliAODForwardMult::GetTriggerString(fTriggerMask) << "\n"
<< " IpZ range: " << fMinIpZ <<" - "<< fMaxIpZ << "\n"
- << " Max Nch: " << fMaxN << "\n"
<< " Use phi acceptance: " << fUsePhiAcc << "\n"
<< " Bins:" << std::endl;
TIter next(&fBins);
//====================================================================
AliForwardMultDists::EtaBin::EtaBin()
: fName(""),
+ fMAxis(1,0,0),
+ fTAxis(1,0,0),
fMinEta(0),
fMaxEta(0),
fMinBin(0),
}
//____________________________________________________________________
-AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta)
+AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta,
+ const TAxis& mAxis)
: fName(""),
+ fMAxis(1,0,0),
+ fTAxis(1,0,0),
fMinEta(minEta),
fMaxEta(maxEta),
fMinBin(0),
fName.ReplaceAll("-", "m");
fName.ReplaceAll("+", "p");
fName.ReplaceAll(".", "d");
+
+ // Copy to other our object
+ mAxis.Copy(fMAxis);
+
+ if (mAxis.GetXbins() && mAxis.GetXbins()->GetArray()) {
+ const TArrayD& mA = *(mAxis.GetXbins());
+ TArrayD tA(mA.GetSize());
+ Int_t j = 0;
+ Double_t min = mA[0];
+ tA[0] = min;
+ for (Int_t i = 1; i < mA.GetSize(); i++) {
+ Double_t d = mA[i] - min;
+ if (d < 1)
+ // Not full integer bin
+ continue;
+ tA[j+1] = tA[j] + Int_t(d);
+ min = tA[j+1];
+ j++;
+ }
+ fTAxis.Set(j, tA.GetArray());
+ }
+ else {
+ // Rounded down maximum and minimum
+ Int_t max = mAxis.GetXmax();
+ Int_t min = mAxis.GetXmin();
+ fTAxis.Set((max-min)+1, min-.5, max+.5);
+ }
}
//____________________________________________________________________
AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
: TObject(o),
fName(o.fName),
+ fMAxis(o.fMAxis),
+ fTAxis(o.fTAxis),
fMinEta(o.fMinEta),
fMaxEta(o.fMaxEta),
fMinBin(o.fMinBin),
if (&o == this) return *this;
fName = o.fName;
+ fMAxis = o.fMAxis;
+ fTAxis = o.fTAxis;
fMinEta = o.fMinEta;
fMaxEta = o.fMaxEta;
fMinBin = o.fMinBin;
return p;
}
+//____________________________________________________________________
+TH1*
+AliForwardMultDists::EtaBin::CreateH1(const char* name,
+ const char* title,
+ const TAxis& xAxis)
+{
+ TH1* ret = 0;
+
+ if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
+ ret = new TH1D(name,title,xAxis.GetNbins(), xAxis.GetXbins()->GetArray());
+ else
+ ret = new TH1D(name,title,xAxis.GetNbins(),xAxis.GetXmin(),xAxis.GetXmax());
+ return ret;
+}
+//____________________________________________________________________
+TH2*
+AliForwardMultDists::EtaBin::CreateH2(const char* name,
+ const char* title,
+ const TAxis& xAxis,
+ const TAxis& yAxis)
+{
+ TH2* ret = 0;
+
+ if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray()) {
+ if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
+ // Variable variable
+ ret = new TH2D(name,title,
+ xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
+ yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
+ else
+ // variable fixed
+ ret = new TH2D(name,title,
+ xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
+ yAxis.GetNbins(),yAxis.GetXmin(),yAxis.GetXmax());
+ }
+ else {
+ if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
+ // fixed variable
+ ret = new TH2D(name,title,
+ xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
+ yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
+ else
+ // fixed fixed
+ ret = new TH2D(name,title,
+ xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
+ yAxis.GetNbins(), yAxis.GetXmin(), yAxis.GetXmax());
+ }
+ return ret;
+}
+
//____________________________________________________________________
void
AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
- UShort_t max, UShort_t nDiv,
Bool_t useMC)
{
/**
fMinBin = hist.GetXaxis()->FindBin(fMinEta);
fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
- TString eta(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
- fSum = new TH1D("rawDist", Form("Raw P(#it{N}_{ch}) in %s", eta.Data()),
- (max+1)*nDiv, -.5, max+.5);
+ TString t(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
+ fSum = CreateH1("rawDist",Form("Raw P(#it{N}_{ch}) in %s",t.Data()),fMAxis);
fSum->SetDirectory(0);
fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
fSum->Sumw2();
l->Add(fSum);
- fCorr = new TH2D("corr", Form("C_{SPD,FMD} in %s", eta.Data()),
- max+2, -1.5, max+.5, max+2, -1.5, max+.5);
+ fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
fCorr->SetDirectory(0);
fCorr->GetXaxis()->SetTitle("Forward");
fCorr->GetYaxis()->SetTitle("Central");
l->Add(fCoverage);
if (!useMC) return;
- fResponse = new TH2D("response", Form("Reponse matrix in %s", eta.Data()),
- nDiv*(max+1), -0.5, max+0.5, max+1, -.5, max+.5);
+ fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
+ fMAxis, fTAxis);
fResponse->SetDirectory(0);
fResponse->SetYTitle("MC");
fResponse->SetXTitle("Signal");
fResponse->SetOption("colz");
l->Add(fResponse);
- fTruth = new TH1D("truth", Form("True P(#it{N}_{ch}) in %s", eta.Data()),
- (max+1), -0.5, max+0.5);
+ fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
fTruth->SetLineColor(kBlack);
fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s",
- eta.Data()));
+ t.Data()));
fTruthAccepted->SetLineColor(kGray+2);
fTruthAccepted->SetFillColor(kOrange+2);
fTruthAccepted->SetDirectory(0);
// Fill with weight
Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
fSum->Fill(sum, w);
- fCorr->Fill(fsum, csum);
+ fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
Int_t nTotal = fMaxBin - fMinBin + 1;
fCoverage->Fill(100*float(covered) / nTotal);
}
//____________________________________________________________________
void
-AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
+AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
{
/**
* Called at the end of the final processing of the job on the
while (*phist) {
TH1* h = *phist;
if (h) {
+ Int_t maxN = h->GetNbinsX();
Double_t intg = h->Integral(1, maxN);
h->Scale(1. / intg, "width");
}
kVertex = 4,
kTriggerVertex = 5
};
+ /**
+ * Structure to define @f$\eta@f$ bins with an @f$ N_{ch}@f$ axis
+ *
+ */
+ struct BinSpec
+ {
+ /**
+ * Constructor
+ *
+ * @param etaMin Low cut on @f$\eta@f$
+ * @param etaMax High cut on @f$\eta@f$
+ * @param nchLow Lowest @f$ N_{ch}@f$ (e.g., -0.5);
+ */
+ BinSpec(Double_t etaMin, Double_t etaMax, Double_t nchLow);
+ /**
+ * Push @a n bins of with @a d in @f$ N_{ch}@f$ onto the axis
+ *
+ * If only a single push is done, then we will get an axis of
+ * equally sized bins (@a d) from @f$ l@f$ edge to @f$ nd+low@f$ -
+ * e.g., if one only does
+ *
+ * @code
+ * BinSpec b(e1, e2, -.5);
+ * b.Push(10, 1);
+ * @endcode
+ *
+ * One ends up with 10 bins from -0.5 to 9.5.
+ *
+ * @param n Number of bins to push
+ * @param d Bin width of each of the bins
+ */
+ void Push(UShort_t n, Double_t d);
+ /**
+ * Get the axis computed from the setup using Push
+ *
+ * @return Reference to the axis
+ */
+ const TAxis& Axis() const;
+ Double_t fEtaMin;
+ Double_t fEtaMax;
+ Double_t fLow;
+ TArrayI fN;
+ TArrayD fD;
+ mutable TAxis fAxis;
+ };
+
/**
* Default constructor
*/
* @param cache
*/
static void ProjectX(const TH2* input, TH1* cache);
+ /**
+ * Add an @f$\eta@f$ bin
+ *
+ * @param spec Bin specification
+ */
+ void AddBin(const BinSpec& spec);
/**
* Add an @f$\eta@f$ bin
*
* @param etaLow Low cut on @f$\eta@f$
* @param etaMax High cut on @f$\eta@f$
+ * @param nAxis Axis to use for measured @f$ N_{ch}@f$
*/
- void AddBin(Double_t etaLow, Double_t etaMax=kInvalidEta);
+ void AddBin(Double_t etaLow, Double_t etaMax, const TAxis& nAxis);
/**
- * Set the maximum @f$N_{ch}@f$ to consider
+ * Add an @f$\eta@f$ bin
*
- * @param maxN Maximum
- * @param nBins (Optionally) the number of bins per particle number
+ * @param etaLow Low cut on @f$\eta@f$
+ * @param etaMax High cut on @f$\eta@f$
+ * @param nMax Maximum @f$ N_{ch}@f$
+ * @param nDiv Number of subdivisions per @f$ N_{ch}@f$
*/
- void SetMaxN(UShort_t maxN, UShort_t nBins=-1);
- void SetNDivisions(UInt_t nDiv) { fNDivisions = nDiv; }
- /**
+ void AddBin(Double_t etaLow, Double_t etaMax, UShort_t nMax, UShort_t nDiv);
+ /**
* Set the range of valid interaction points
*
* @param z1 Least Z coordinate
*
* @param minEta Least @f$\eta@f$ to consider
* @param maxEta Largest @f$\eta@f$ to consider
+ * @param mAxis The @f$ N_{ch}@f$ axis to use for measured data
+ * @param tAxis The @f$ N_{ch}@f$ axis to use for truth data
*/
- EtaBin(Double_t minEta, Double_t maxEta);
+ EtaBin(Double_t minEta, Double_t maxEta, const TAxis& mAxis);
/**
* Copy constructor
*
* @return Container, or null
*/
TList* FindParent(TList* l, Bool_t create=true) const;
+ /**
+ * Create a 1D histogram with specified axis
+ *
+ * @param name Name of histogram
+ * @param title Title of histogram
+ * @param xAxis X-axis to use
+ *
+ * @return Created histogram
+ */
+ static TH1* CreateH1(const char* name, const char* title,
+ const TAxis& xAxis);
+ /**
+ * Create a 2D histogram with specified axis
+ *
+ * @param name Name of histogram
+ * @param title Title of histogram
+ * @param xAxis X-axis to use
+ * @param yAxis Y-axis to use
+ *
+ * @return Created histogram
+ */
+ static TH2* CreateH2(const char* name, const char* title,
+ const TAxis& xAxis, const TAxis& yAxis);
/**
* Set-up internal structures on first event.
*
* @param list List to add information to
* @param hist Template histogram
- * @param max Maximum number of particles
- * @param nDiv Number of divisions per charged particle bin
* @param useMC Whether to set-up for MC input
*/
- void SetupForData(TList* list, const TH2& hist, UShort_t max,
- UShort_t nDiv, Bool_t useMC);
+ void SetupForData(TList* list, const TH2& hist, Bool_t useMC);
/**
* Process a single event
*
*
* @param in Input list
* @param out Output list
- * @param maxN Maximum number of @f$N_{ch}@f$ to consider
*/
- void Terminate(TList* in, TList* out, UShort_t maxN);
+ void Terminate(TList* in, TList* out);
TString fName; // Name of this bin
+ TAxis fMAxis; // Axis used for measured Nch
+ TAxis fTAxis; // Axis used for true Nch
Double_t fMinEta; // Least @f$\eta@f$ to consider
Double_t fMaxEta; // Largest @f$\eta@f$ to consider
Int_t fMinBin; // Least @f$\eta@f$ bin to consider
TH1* fTruthAccepted; // `true' distribution for accepted events
TH1* fCoverage; // How much was covered
- ClassDef(EtaBin,1);
+ ClassDef(EtaBin,2);
};
TList fBins; // List of bins
TList* fSymmetric; // Bins symmetric around 0
TH1* fForwardCache; // Projection cache
TH1* fCentralCache; // Projection cache
TH1* fMCCache; // Projection cache
- UShort_t fMaxN; // Maximum of @f$N_{ch}@f$
- UShort_t fNDivisions; // Number of particle number sub-divions
Bool_t fUsePhiAcc; // If true, scale by phi acceptance
ClassDef(AliForwardMultDists,1);
// Get some info from the input collection
UShort_t sys;
UShort_t sNN;
- UShort_t maxN;
ULong_t trig;
Double_t minZ;
Double_t maxZ;
GetParameter(mTop, "sys", sys);
GetParameter(mTop, "snn", sNN);
- GetParameter(mTop, "maxN", maxN);
GetParameter(mTop, "trigger", trig);
GetParameter(mTop, "minIpZ", minZ);
GetParameter(mTop, "maxIpZ", maxZ);
if (mId == 0xDeadBeef) return;
// Store information
- SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,maxN);
+ SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
+ corrFile.IsNull());
// Load other data
TString savPath(gROOT->GetMacroPath());
SaveSummarize();
}
+ /**
+ * Append an & to a string and the next term.
+ *
+ * @param trg Output string
+ * @param what Term
+ */
static void AppendAnd(TString& trg, const TString& what)
{
if (!trg.IsNull()) trg.Append(" & ");
* @param trigger Trigger mask
* @param minIpZ Least z coordinate of interaction point
* @param maxIpZ Largest z coordinate of interaction point
- * @param maxN Largest @f$N_{ch}@f$ to consider
*/
void SaveInformation(TDirectory* dir,
const TString& method,
UInt_t trigger,
Double_t minIpZ,
Double_t maxIpZ,
- UShort_t maxN) const
+ Bool_t self) const
{
dir->cd();
+ TParameter<bool>* pM = new TParameter<bool>("self", self);
+ pM->SetBit(BIT(19));
+ pM->Write();
+
TNamed* outMeth = new TNamed("method", method.Data());
outMeth->SetUniqueID(mId);
outMeth->Write();
TParameter<double>* pR = new TParameter<double>("regParam", regParam);
+ pR->SetBit(BIT(19));
pR->Write();
TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
pT->Write();
TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
+ pY->SetBit(BIT(19));
pY->Write();
TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
+ pZ->SetBit(BIT(19));
pZ->Write();
-
- TParameter<int>* pN = new TParameter<int>("maxN", maxN);
- pN->Write();
}
/**
* Save a script to do a summary of this step
"Unfolded P(#it{N}_{ch})");
THStack* allCorrected = new THStack("corrected",
"Corrected P(#it{N}_{ch})");
+ THStack* allRatio = (sys != 1 ? 0 :
+ new THStack("ratios", "Ratios to other"));
TMultiGraph* allALICE = (sys != 1 ? 0 :
new TMultiGraph("alice", "ALICE Published"));
TMultiGraph* allCMS = (sys != 1 ? 0 :
THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
if (!binS) continue;
+ TH1* result = 0;
Bin2Stack(binS, i, allMeasured, allTruth, allTruthA,
- allUnfolded, allCorrected);
- Other2Stack(o->GetName(), i, sNN, allALICE, allCMS);
+ allUnfolded, allCorrected, result);
+
+ TGraph* alice = 0;
+ TGraph* cms = 0;
+ Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
+ Ratio2Stack(i, result, alice, cms, allRatio);
i++;
}
dir->Add(allMeasured);
else
delete allCMS;
}
+ if (allRatio && allRatio->GetHists()) {
+ if (allRatio->GetHists()->GetEntries() > 0)
+ dir->Add(allRatio);
+ else
+ delete allRatio;
+ }
}
/**
* Process a single eta bin
THStack* truth,
THStack* accepted,
THStack* unfolded,
- THStack* corrected)
+ THStack* corrected,
+ TH1*& result)
{
Int_t open, closed;
Double_t factor;
cln->SetMarkerStyle(sty);
cln->SetMarkerSize(size);
cln->Scale(factor); // Scale by 10^i
+ if (col == kColorCorrected) result = cln;
// Make sure we do not get the old legend
TObject* tst = cln->FindObject("legend");
* @param allCMS Stack of CMS data
*/
void Other2Stack(const TString& name, Int_t i,
- UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS)
+ UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
+ TGraph*& alice, TGraph*& cms)
{
if (!allALICE && !allCMS) return;
tmp.ReplaceAll("p", "+");
tmp.ReplaceAll("m", "-");
tmp.ReplaceAll("_", " ");
+ tmp.ReplaceAll("d", ".");
TObjArray* tokens = tmp.Tokenize(" ");
if (!tokens || tokens->GetEntriesFast() < 2) {
Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
g->SetMarkerColor(kColorALICE);
g->SetMarkerSize(size);
allALICE->Add(g, "p same");
+ alice = g;
}
}
if (allCMS) {
g->SetMarkerColor(kColorCMS);
g->SetMarkerSize(size);
allCMS->Add(g, "p same");
+ cms = g;
+ }
+ }
+ }
+ /**
+ * Create ratios to other data
+ *
+ * @param i
+ * @param res
+ * @param alice
+ * @param cms
+ * @param all
+ */
+ void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
+ {
+ if (!all || !res || !(alice || cms)) return;
+
+ Int_t off = 5*ib;
+ TGraph* gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
+ TGraph** pg = gs;
+ while (*pg) {
+ TGraph* g = *pg;
+ const char* n = (g == alice ? "ALICE" : "CMS");
+
+ TH1* r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
+ TString tit(r->GetTitle());
+ tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
+ r->SetTitle(tit);
+ r->SetMarkerColor(g->GetMarkerColor());
+ r->SetLineColor(g->GetLineColor());
+
+ TObject* tst = r->FindObject("legend");
+ if (tst) r->GetListOfFunctions()->Remove(tst);
+
+ for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
+ Double_t c = r->GetBinContent(i);
+ Double_t e = r->GetBinError(i);
+ Double_t o = g->Eval(r->GetBinCenter(i));
+ if (o < 1e-12) {
+ r->SetBinContent(i, 0);
+ r->SetBinError(i, 0);
+ continue;
+ }
+ r->SetBinContent(i, (c - o) / o + off);
+ r->SetBinError(i, e / o);
}
+ all->Add(r);
+ pg++;
}
+ TLegend* leg = StackLegend(all);
+ if (!leg) return;
+
+ TString txt = res->GetTitle();
+ txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
+ if (ib == 0) txt.Append(" "); // (#times1)");
+ // else if (ib == 1) txt.Append(" (#times10)");
+ else txt.Append(Form(" (+%d)", off));
-
+ TObject* dummy = 0;
+ TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
+ e->SetMarkerStyle(res->GetMarkerStyle());
+ e->SetMarkerSize(res->GetMarkerSize());
+ e->SetMarkerColor(kBlack);
+ e->SetFillColor(0);
+ e->SetFillStyle(0);
+ e->SetLineColor(kBlack);
}
+
/**
* Get or create a stack legend. This is done by adding a TLegend
* object to the list of functions for the first histogram in the
TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
Int_t factor)
{
- TString oC = Form("OtherPNch::GetData(%hu,%f,%hu)",
+ TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);",
type, eta, sNN);
TGraphAsymmErrors* g =
reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
if (!g) {
- Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
- type, eta, sNN);
+ // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
+ // type, eta, sNN);
return 0;
}
-
+
for (Int_t j = 0; j < g->GetN(); j++) {
g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
#pragma link C++ class AliForwardMultiplicityDistribution::Bin+;
#pragma link C++ class AliForwardMultDists+;
#pragma link C++ class AliForwardMultDists::EtaBin+;
+#pragma link C++ class AliForwardMultDists::BinSpec+;
#else
# error Not for compilation