In general, more diagnostics histograms.
Event inspector can now flag #_SPD_Cluster > 0 events
New implementation of sharing algorithm that zero's candidates
if not used.
dN/deta tasks now outputs summary stacks.
fHist(),
fTriggers(0),
fIpZ(fgkInvalidIpZ),
- fCentrality(-1)
+ fCentrality(-1),
+ fNClusters(0)
{
//
// Constructor
200, -4, 6, 20, 0, 2*TMath::Pi()),
fTriggers(0),
fIpZ(fgkInvalidIpZ),
- fCentrality(-1)
+ fCentrality(-1),
+ fNClusters(0)
{
//
// Constructor
// option Passed to TH1::Reset
//
fHist.Reset(option);
- fTriggers = 0;
- fIpZ = fgkInvalidIpZ;
+ fTriggers = 0;
+ fIpZ = fgkInvalidIpZ;
+ fNClusters = 0;
}
//____________________________________________________________________
void
static TObjString ipz;
static TObjString trg;
static TObjString cnt;
+ static TObjString ncl;
ipz = Form("ip_z=%fcm", fIpZ);
trg = GetTriggerString(fTriggers);
cnt = Form("%+6.1f%%", fCentrality);
+ ncl = Form("%d clusters", fNClusters);
b->Add(&fHist);
b->Add(&ipz);
b->Add(&trg);
b->Add(&cnt);
+ b->Add(&ncl);
}
//____________________________________________________________________
if ((mask & kC) != 0x0) trg.Append("C ");
if ((mask & kE) != 0x0) trg.Append("E ");
if ((mask & kMCNSD) != 0x0) trg.Append("MCNSD ");
+ if ((mask & kNClusterGt0) != 0x0) trg.Append("NCluster>0 ");
return trg.Data();
}
ret->GetXaxis()->SetBinLabel(kBinMCNSD, "NSD (MC truth)");
ret->GetXaxis()->SetBinLabel(kBinPileUp, "w/Pileup");
ret->GetXaxis()->SetBinLabel(kBinOffline, "w/Offline");
+ ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1");
ret->GetXaxis()->SetBinLabel(kWithVertex, "w/Vertex");
ret->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
ret->GetXaxis()->SetBinLabel(kAccepted, "Accepted by cut");
else if (s.CompareTo("A") == 0) trgMask |= AliAODForwardMult::kA;
else if (s.CompareTo("C") == 0) trgMask |= AliAODForwardMult::kC;
else if (s.CompareTo("E") == 0) trgMask |= AliAODForwardMult::kE;
+ else if (s.CompareTo("NCluster>0") == 0)
+ trgMask |= AliAODForwardMult::kNClusterGt0;
else
AliWarningGeneral("MakeTriggerMask",
Form("Unknown trigger %s", s.Data()));
if (IsTriggerBits(kPileUp)) hist->AddBinContent(kBinPileUp);
if (IsTriggerBits(kMCNSD)) hist->AddBinContent(kBinMCNSD);
if (IsTriggerBits(kOffline)) hist->AddBinContent(kBinOffline);
+ if (IsTriggerBits(kNClusterGt0)) hist->AddBinContent(kBinNClusterGt0);
}
// Check if we have an event of interest.
if (!IsTriggerBits(triggerMask|kB)) return false;
/** true NSD from MC */
kMCNSD = 0x400,
/** Offline MB triggered */
- kOffline = 0x800
+ kOffline = 0x800,
+ /** At least one SPD cluster */
+ kNClusterGt0 = 0x1000
};
/**
* Bin numbers in trigger histograms
kBinPileUp,
kBinMCNSD,
kBinOffline,
+ kBinNClusterGt0,
kWithTrigger,
kWithVertex,
kAccepted
*/
void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
/**
- * Check if bit(s) are set in the trigger mask
+ * Check if all bit(s) are set in the trigger mask. Note, this is
+ * an @e and between the bits. If you need an @e or you should use
+ * the member function IsTriggerOrBits
*
* @param bits Bits to test for
*
- * @return
+ * @return true if all enabled bits in the argument is also set in
+ * the trigger word
*/
Bool_t IsTriggerBits(UInt_t bits) const;
+ /**
+ * Check if any of bit(s) are enabled in the trigger word. This is
+ * an @e or between the selected bits. If you need and @a and you
+ * should use the member function IsTriggerBits;
+ *
+ * @param bits Bits to check for
+ *
+ * @return true if any of the enabled bits in the arguments are also
+ * enabled in the trigger mask
+ */
+ Bool_t IsTriggerOrBits(UInt_t bits) const;
/**
* Whether we have any trigger bits
*/
* @return
*/
Bool_t HasCentrality() const { return !(fCentrality < 0); }
-
+ /**
+ * Get the number of SPD clusters seen in @f$ |\eta|<1@f$
+ *
+ * @return Number of SPD clusters seen
+ */
+ UShort_t GetNClusters() const { return fNClusters; }
+ /**
+ * Set the number of SPD clusters seen in @f$ |\eta|<1@f$
+ *
+ * @param n Number of SPD clusters
+ */
+ void SetNClusters(UShort_t n) { fNClusters = n; }
/**
* Get the name of the object
*
*/
static UInt_t MakeTriggerMask(const char* what);
protected:
- Bool_t fIsMC; // Whether this is from MC
- TH2D fHist; // Histogram of d^2N_{ch}/(deta dphi) for this event
- UInt_t fTriggers; // Trigger bit mask
- Float_t fIpZ; // Z coordinate of the interaction point
- Float_t fCentrality; // Event centrality
+ Bool_t fIsMC; // Whether this is from MC
+ TH2D fHist; // Histogram of d^2N_{ch}/(deta dphi) for this event
+ UInt_t fTriggers; // Trigger bit mask
+ Float_t fIpZ; // Z coordinate of the interaction point
+ Float_t fCentrality; // Event centrality
+ UShort_t fNClusters; // Number of SPD clusters in |eta|<1
static const Float_t fgkInvalidIpZ; // Invalid IpZ value
- ClassDef(AliAODForwardMult,2); // AOD forward multiplicity
+ ClassDef(AliAODForwardMult,3); // AOD forward multiplicity
};
//____________________________________________________________________
{
return HasTrigger() && ((fTriggers & bits) == bits);
}
+//____________________________________________________________________
+inline Bool_t
+AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const
+{
+ return HasTrigger() && ((fTriggers & bits) != 0);
+}
#endif
// Local Variables:
#include "AliForwardUtil.h"
#include "AliAODForwardMult.h"
#include <TFile.h>
+#include <TStyle.h>
//____________________________________________________________________
AliBasedNdetaTask::AliBasedNdetaTask()
// high High cut
//
CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
- AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
fListOfCentralities->AddAtAndExpand(bin, at);
}
// Loop over centrality bins
TIter next(fListOfCentralities);
CentralityBin* bin = 0;
- while ((bin = static_cast<CentralityBin*>(next())))
+ gStyle->SetPalette(1);
+ THStack* dndetaStack = new THStack("dndeta", "dN/d#eta");
+ THStack* dndetaStackRebin = new THStack(Form("dndeta_rebin%02d", fRebin),
+ "dN_{ch}/d#eta");
+ THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta");
+ THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
+ "dN_{ch}/d#eta");
+
+ while ((bin = static_cast<CentralityBin*>(next()))) {
bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
- fTriggerMask, GetColor(), GetMarker());
+ fTriggerMask, GetMarker());
+ if (fCentAxis && bin->IsAllBin()) continue;
+ TH1* dndeta = bin->GetResult(0, false, "");
+ TH1* dndetaSym = bin->GetResult(0, true, "");
+ TH1* dndetaMC = bin->GetResult(0, false, "MC");
+ TH1* dndetaMCSym = bin->GetResult(0, true, "MC");
+ if (dndeta) dndetaStack->Add(dndeta);
+ if (dndetaSym) dndetaStack->Add(dndetaSym);
+ if (dndetaMC) dndetaMCStack->Add(dndetaMC);
+ if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
+ if (fRebin > 1) {
+ dndeta = bin->GetResult(fRebin, false, "");
+ dndetaSym = bin->GetResult(fRebin, true, "");
+ dndetaMC = bin->GetResult(fRebin, false, "MC");
+ dndetaMCSym = bin->GetResult(fRebin, true, "MC");
+ if (dndeta) dndetaStackRebin->Add(dndeta);
+ if (dndetaSym) dndetaStackRebin->Add(dndetaSym);
+ if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC);
+ if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
+ }
+ }
+ // Output the stack
+ fOutput->Add(dndetaStack);
+
+ // If available output rebinned stack
+ if (!dndetaStackRebin->GetHists() ||
+ dndetaStackRebin->GetHists()->GetEntries() <= 0) {
+ AliWarning("No rebinned histograms found");
+ delete dndetaStackRebin;
+ dndetaStackRebin = 0;
+ }
+ if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
+
+ // If available, output track-ref stack
+ if (!dndetaMCStack->GetHists() ||
+ dndetaMCStack->GetHists()->GetEntries() <= 0) {
+ AliWarning("No MC histograms found");
+ delete dndetaMCStack;
+ dndetaMCStack = 0;
+ }
+ if (dndetaMCStack) fOutput->Add(dndetaMCStack);
+
+ // If available, output rebinned track-ref stack
+ if (!dndetaMCStackRebin->GetHists() ||
+ dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
+ AliWarning("No rebinned MC histograms found");
+ delete dndetaMCStackRebin;
+ dndetaMCStackRebin = 0;
+ }
+ if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
// Output collision energy string
if (fSNNString) fOutput->Add(fSNNString->Clone());
return *this;
}
//____________________________________________________________________
+Int_t
+AliBasedNdetaTask::CentralityBin::GetColor() const
+{
+ if (IsAllBin()) return kRed+2;
+ Float_t fc = (fLow+double(fHigh-fLow)/2) / 100;
+ Int_t nCol = gStyle->GetNumberOfColors();
+ Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
+ Int_t col = gStyle->GetColorPalette(icol);
+ return col;
+}
+//____________________________________________________________________
const char*
AliBasedNdetaTask::CentralityBin::GetListName() const
{
return scaler;
}
+//________________________________________________________________________
+const char*
+AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
+ Bool_t sym,
+ const char* postfix) const
+{
+ static TString n;
+ n = Form("dndeta%s%s",GetName(), postfix);
+ if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
+ if (sym) n.Append("_mirror");
+ return n.Data();
+}
+//________________________________________________________________________
+TH1*
+AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
+ Bool_t sym,
+ const char* postfix) const
+{
+ if (!fOutput) {
+ AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
+ fLow, fHigh));
+ return 0;
+ }
+ TString n = GetResultName(rebin, sym, postfix);
+ TObject* o = fOutput->FindObject(n.Data());
+ if (!o) {
+ // AliWarning(Form("Object %s not found in output list", n.Data()));
+ return 0;
+ }
+ return static_cast<TH1*>(o);
+}
+
//________________________________________________________________________
void
AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
bool symmetrice,
Int_t rebin,
bool cutEdges,
- Int_t color,
Int_t marker)
{
//
copy->Scale(1., "width");
// --- Set some histogram attributes -------------------------------
- SetHistogramAttributes(dndeta, color, marker, Form("ALICE %s", GetName()));
- SetHistogramAttributes(accNorm, color, marker, Form("ALICE %s normalisation",
- GetName()));
+ TString post;
+ if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
+ SetHistogramAttributes(dndeta, GetColor(), marker,
+ Form("ALICE %s%s", GetName(), post.Data()));
+ SetHistogramAttributes(accNorm, GetColor(), marker,
+ Form("ALICE %s normalisation%s",
+ GetName(), post.Data()));
// --- Make symmetric extensions and rebinnings --------------------
if (symmetrice) fOutput->Add(Symmetrice(dndeta));
Bool_t corrEmpty,
Bool_t cutEdges,
Int_t triggerMask,
- Int_t color,
Int_t marker)
{
//
// --- Make result and store ---------------------------------------
MakeResult(fSum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
- scaler, symmetrice, rebin, cutEdges, color, marker);
+ scaler, symmetrice, rebin, cutEdges, marker);
// --- Process result from TrackRefs -------------------------------
- if (fSumMC) {
+ if (fSumMC)
MakeResult(fSumMC, "MC", rootProj, corrEmpty,
(scheme & kShape) ? shapeCorr : 0,
- scaler, symmetrice, rebin, cutEdges, color+2, marker);
- }
+ scaler, symmetrice, rebin, cutEdges, marker+2);
// Temporary stuff
- if (!IsAllBin()) return;
-
- TFile* forward = TFile::Open("forward.root", "READ");
- if (!forward) {
- AliWarning(Form("No forward.root file found"));
- return;
- }
-
- TH1D* shapeCorrProj = 0;
- if (shapeCorr) {
- shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
- shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
- shapeCorrProj->SetDirectory(0);
- fOutput->Add(shapeCorrProj);
- }
-
- TList* official = static_cast<TList*>(forward->Get("official"));
- if (official) {
- TH1F* histEta = static_cast<TH1F*>(official->FindObject("fHistEta"));
- if (histEta) {
- TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
- histEta->GetXaxis()->GetXmin(),
- histEta->GetXaxis()->GetXmax());
- for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
- oEta->SetBinContent(i, histEta->GetBinContent(i));
- oEta->SetBinError(i, histEta->GetBinError(i));
- }
- if (shapeCorrProj) oEta->Divide(shapeCorrProj);
- oEta->Scale(1./ntotal, "width");
- oEta->SetDirectory(0);
- oEta->SetMarkerStyle(marker+4);
- oEta->SetMarkerColor(color+5);
- fOutput->Add(oEta);
- fOutput->Add(Rebin(oEta, rebin, false));
- }
- else
- AliWarning(Form("Couldn't find histogram fHistEta in list %s",
- official->GetName()));
- }
- else
- AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
-
- TList* tracks = static_cast<TList*>(forward->Get("tracks"));
- if (tracks) {
- TH1F* histEta = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
- if (histEta) {
- TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
- histEta->GetXaxis()->GetXmin(),
- histEta->GetXaxis()->GetXmax());
- for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
- oEta->SetBinContent(i, histEta->GetBinContent(i));
- oEta->SetBinError(i, histEta->GetBinError(i));
- }
- if (shapeCorrProj) oEta->Divide(shapeCorrProj);
- oEta->Scale(1./ntotal, "width");
- oEta->SetDirectory(0);
- oEta->SetMarkerStyle(marker);
- oEta->SetMarkerColor(color+5);
- fOutput->Add(oEta);
- fOutput->Add(Rebin(oEta, rebin, false));
- }
- else
- AliWarning(Form("Couldn't find histogram fHistEta in list %s",
- tracks->GetName()));
- }
- else
- AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
+ // if (!IsAllBin()) return;
- forward->Close();
}
//
* @param symmetrice Whether to make symmetric extensions
* @param rebin Whether to rebin
* @param cutEdges Whether to cut edges when rebinning
- * @param color Color of markers
* @param marker Marker style
*/
virtual void MakeResult(const TH2D* sum,
bool symmetrice,
Int_t rebin,
bool cutEdges,
- Int_t color,
Int_t marker);
/**
* End of processing
* @param corrEmpty Whether to correct for empty bins
* @param cutEdges Whether to cut edges when rebinning
* @param triggerMask Trigger mask
- * @param color Base colour for markers
* @param marker Marker style
*/
virtual void End(TList* sums,
Bool_t corrEmpty,
Bool_t cutEdges,
Int_t triggerMask,
- Int_t color,
Int_t marker);
/**
* @{
*/
TH1I* GetTrigggers() { return fTriggers; }
/** @} */
+
+ Int_t GetColor() const;
+ TList* GetResults() const { return fOutput; }
+ const char* GetResultName(Int_t rebin, Bool_t sym,
+ const char* postfix="") const;
+ TH1* GetResult(Int_t rebin, Bool_t sym,
+ const char* postfix="") const;
+
protected:
/**
* Create sum histogram
//____________________________________________________________________
#define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f,
+ Bool_t includeSigma) const
+{
+ //
+ // Return
+ // Delta - f * (xi + sigma)
+ return fDelta - f * (fXi + (includeSigma ? fSigma : 0));
+}
+
//____________________________________________________________________
void
AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
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));
+ // 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));
+ // 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));
+ // 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));
+ // 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));
return static_cast<TObjArray*>(fRings.At(idx));
}
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Int_t etabin,
+ Double_t f, Bool_t showErrors,
+ Bool_t includeSigma) const
+{
+ ELossFit* fit = GetFit(d, r, etabin);
+ if (!fit) {
+ if (showErrors) {
+ AliWarning(Form("No fit for FMD%d%c @ etabin=%d", d, r, etabin));
+ }
+ return -1024;
+ }
+ return fit->GetLowerBound(f, includeSigma);
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Double_t eta,
+ Double_t f, Bool_t showErrors,
+ Bool_t includeSigma) const
+{
+ Int_t bin = FindEtaBin(eta);
+ if (bin <= 0) {
+ if (showErrors)
+ AliError(Form("eta=%f out of bounds for FMD%d%c", eta, d, r));
+ return -1024;
+ }
+ return GetLowerBound(d, r, bin, f, showErrors, includeSigma);
+}
+
+//____________________________________________________________________
namespace {
TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
Int_t color)
}
}
+
+
#define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O')
#define IDX2DET(I) (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3))
//____________________________________________________________________
* @return
*/
const Char_t* GetName() const;
+ /**
+ * Calculate the lower bound
+ *
+ * @param f Width factor
+ * @param includeSigma Whether to include sigma
+ *
+ * @return @f$ \Delta - f (\xi + \sigma)@f$
+ */
+ Double_t GetLowerBound(Double_t f, Bool_t includeSigma) const;
/**
* Calculate the quality
*/
ELossFit* GetFit(UShort_t d, Char_t r, Int_t etabin) const;
/* @} */
+ /**
+ * @{
+ * @name Lower bounds on fits
+ */
+ /**
+ * Get the lower validity bound of the fit
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param etaBin Eta bin (1-based)
+ * @param f Factor on xi (and sigma)
+ * @param showErrors Show errors
+ * @param includeSigma Whether to include sigma
+ *
+ * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
+ */
+ Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin,
+ Double_t f, Bool_t showErrors=true,
+ Bool_t includeSigma=true) const;
+ /**
+ * Get the lower validity bound of the fit
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta value
+ * @param f Factor on xi (and sigma)
+ * @param showErrors Show errors
+ * @param includeSigma Whether to include sigma
+ *
+ * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
+ */
+ Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta,
+ Double_t f, Bool_t showErrors=true,
+ Bool_t includeSigma=true) const;
+
+
/**
* @{
* @name Miscellaneous
#include "AliLog.h"
#include <TH2D.h>
#include <TROOT.h>
+#include <THStack.h>
#include <iostream>
#include <iomanip>
if (fUseSecondaryMap) {
TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
if (!bg) {
- AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
- d, r, uvb));
+ AliWarning(Form("No secondary correction for FMDM%d%c "
+ "in vertex bin %d", d, r, uvb));
continue;
}
// Divide by primary/total ratio
if (fUseAcceptance) {
TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
if (!ac) {
- AliWarning(Form("No acceptance correction for FMD%d%c in vertex bin %d",
- d, r, uvb));
+ AliWarning(Form("No acceptance correction for FMD%d%c in "
+ "vertex bin %d", d, r, uvb));
continue;
}
// Divide by the acceptance correction
h->Divide(ac);
}
-
-
-
-
if (fUseMergingEfficiency) {
if (!fcm.GetMergingEfficiency()) {
AliWarning("No merging efficiencies");
TIter next(&fRingHistos);
RingHistos* o = 0;
- while ((o = static_cast<RingHistos*>(next())))
+ THStack* sums = new THStack("sums", "Sums of ring results");
+ while ((o = static_cast<RingHistos*>(next()))) {
o->ScaleHistograms(d, nEvents);
+ TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
+ o->fDensity->GetNbinsY(),"e");
+ sum->Scale(1., "width");
+ sum->SetTitle(o->GetName());
+ sum->SetDirectory(0);
+ sum->SetYTitle("#sum N_{ch,primary}");
+ sums->Add(sum);
+ }
+ d->Add(sums);
+
}
//____________________________________________________________________
void
// d detector
// r ring
//
- fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r),
- Form("in FMD%d%c", d, r),
+ fDensity = new TH2D("primaryDensity",
+ "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)",
200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
0, 2*TMath::Pi());
fDensity->SetDirectory(0);
+ fDensity->SetMarkerColor(Color());
+ fDensity->Sumw2();
fDensity->SetXTitle("#eta");
fDensity->SetYTitle("#phi [radians]");
fDensity->SetZTitle("Primary N_{ch} density");
//____________________________________________________________________
void
-AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t /*nEvents*/)
+AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
{
//
// Scale the histograms to the total number of events
TList* l = GetOutputList(dir);
if (!l) return;
- //TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
- //if (density) density->Scale(1./nEvents);
+ TH1* density = GetOutputHist(l,"primaryDensity");
+ if (density) density->Scale(1./nEvents);
}
//____________________________________________________________________
*
* @param option Not used
*/
- void Print(Option_t* option="") const;
+ virtual void Print(Option_t* option="") const;
protected:
/**
* Internal data structure to keep track of the histograms
#include "AliLog.h"
#include <TH2D.h>
#include <TProfile.h>
+#include <THStack.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>
: TNamed(),
fRingHistos(),
fMultCut(0),
+ fNXi(1),
+ fIncludeSigma(true),
fSumOfWeights(0),
fWeightedSum(0),
fCorrections(0),
fMaxParticles(5),
fUsePoisson(false),
+ fUsePhiAcceptance(false),
fAccI(0),
fAccO(0),
fFMD1iMax(0),
fFMD2oMax(0),
fFMD3iMax(0),
fFMD3oMax(0),
+ fMaxWeights(0),
+ fLowCuts(0),
fEtaLumping(1),
fPhiLumping(1),
fDebug(0)
: TNamed("fmdDensityCalculator", title),
fRingHistos(),
fMultCut(0),
+ fNXi(1),
+ fIncludeSigma(true),
fSumOfWeights(0),
fWeightedSum(0),
fCorrections(0),
fMaxParticles(5),
fUsePoisson(false),
+ fUsePhiAcceptance(false),
fAccI(0),
fAccO(0),
fFMD1iMax(0),
fFMD2oMax(0),
fFMD3iMax(0),
fFMD3oMax(0),
+ fMaxWeights(0),
+ fLowCuts(0),
fEtaLumping(5),
fPhiLumping(1),
fDebug(0)
fAccI = GenerateAcceptanceCorrection('I');
fAccO = GenerateAcceptanceCorrection('O');
+
+ fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
+ 1, 0, 1, 1, 0, 1);
+ fMaxWeights->SetXTitle("#eta");
+ fMaxWeights->SetDirectory(0);
+
+ fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
+ fLowCuts->SetXTitle("#eta");
+ fLowCuts->SetDirectory(0);
+
+ for (Int_t i = 0; i < 5; i++) fMultCuts[i] = 0;
}
//____________________________________________________________________
: TNamed(o),
fRingHistos(),
fMultCut(o.fMultCut),
+ fNXi(o.fNXi),
+ fIncludeSigma(o.fIncludeSigma),
fSumOfWeights(o.fSumOfWeights),
fWeightedSum(o.fWeightedSum),
fCorrections(o.fCorrections),
fMaxParticles(o.fMaxParticles),
fUsePoisson(o.fUsePoisson),
+ fUsePhiAcceptance(o.fUsePhiAcceptance),
fAccI(o.fAccI),
fAccO(o.fAccO),
fFMD1iMax(o.fFMD1iMax),
fFMD2oMax(o.fFMD2oMax),
fFMD3iMax(o.fFMD3iMax),
fFMD3oMax(o.fFMD3oMax),
+ fMaxWeights(o.fMaxWeights),
+ fLowCuts(o.fLowCuts),
fEtaLumping(o.fEtaLumping),
fPhiLumping(o.fPhiLumping),
fDebug(o.fDebug)
TNamed::operator=(o);
fMultCut = o.fMultCut;
+ fNXi = o.fNXi;
+ fIncludeSigma = o.fIncludeSigma;
fDebug = o.fDebug;
fMaxParticles = o.fMaxParticles;
fUsePoisson = o.fUsePoisson;
+ fUsePhiAcceptance = o.fUsePhiAcceptance;
fAccI = o.fAccI;
fAccO = o.fAccO;
fFMD1iMax = o.fFMD1iMax;
fFMD2oMax = o.fFMD2oMax;
fFMD3iMax = o.fFMD3iMax;
fFMD3oMax = o.fFMD3oMax;
+ fMaxWeights = o.fMaxWeights;
+ fLowCuts = o.fLowCuts;
fEtaLumping = o.fEtaLumping;
fPhiLumping = o.fPhiLumping;
fRingHistos.Delete();
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
-
+
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::SetMultCuts(Double_t fmd1i,
+ Double_t fmd2i,
+ Double_t fmd2o,
+ Double_t fmd3i,
+ Double_t fmd3o)
+{
+ fMultCuts[0] = fmd1i;
+ fMultCuts[1] = fmd2i;
+ fMultCuts[2] = fmd2o;
+ fMultCuts[3] = fmd3i;
+ fMultCuts[4] = fmd3o;
+ fMultCut = (fmd1i+fmd2i+fmd2o+fmd3i+fmd3o) / 5;
+}
+
//____________________________________________________________________
Double_t
-AliFMDDensityCalculator::GetMultCut() const
+AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
+ Bool_t errors) const
{
//
// Get the multiplicity cut. If the user has set fMultCut (via
// Return:
// Lower cut on multiplicity
//
+ Int_t idx = (d == 1 ? 0 : 2*(d - 2) + 1 + ((r=='I' || r=='i') ? 0 : 1));
+ if (fMultCuts[idx] > 0) return fMultCuts[idx];
if (fMultCut > 0) return fMultCut;
AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
AliFMDCorrELossFit* fits = fcm.GetELossFit();
- return fits->GetLowCut();
+ if (fNXi < 0) return fits->GetLowCut();
+
+ return fits->GetLowerBound(d, r, ieta, fNXi, errors, fIncludeSigma);
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
+ Bool_t /*errors*/) const
+{
+ //
+ // 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
+ //
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+ AliFMDCorrELossFit* fits = fcm.GetELossFit();
+ Int_t iEta = fits->FindEtaBin(eta);
+
+ return GetMultCut(d, r, iEta);
}
//____________________________________________________________________
for (UShort_t s=0; s<ns; s++) {
for (UShort_t t=0; t<nt; t++) {
- Float_t mult = fmd.Multiplicity(d,r,s,t);
- Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
- Float_t eta = fmd.Eta(d,r,s,t);
+ Float_t mult = fmd.Multiplicity(d,r,s,t);
+ Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+ Float_t eta = fmd.Eta(d,r,s,t);
+ Double_t cut = 1024;
+ if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
rh->fTotalStrips->Fill(eta, phi);
- if (mult == AliESDFMD::kInvalidMult || mult > 20) {
+ if (mult == AliESDFMD::kInvalidMult || mult > 20) {
rh->fEmptyStrips->Fill(eta,phi);
+ rh->fEvsM->Fill(mult,0);
continue;
}
-
- Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
+
+ Double_t n = 0;
+ if (cut > 0 && mult > cut)
+ 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);
+ Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
fCorrections->Fill(c);
if (c > 0) n /= c;
rh->fEvsM->Fill(mult,n);
else rh->fEmptyStrips->Fill(eta,phi);
h->Fill(eta,phi,n);
- rh->fDensity->Fill(eta,phi,n);
+ if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
} // for t
} // for s
for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
Double_t eLossV = h->GetBinContent(ieta, iphi);
- Float_t eta = h->GetXaxis()->GetBinCenter(ieta);
- Float_t phi = h->GetYaxis()->GetBinCenter(iphi);
+ Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
+ Double_t phi = h->GetYaxis()->GetBinCenter(iphi);
Int_t jeta = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
Int_t jphi = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
Double_t empty = rh->fEmptyStrips->GetBinContent(jeta, jphi);
Double_t total = rh->fTotalStrips->GetBinContent(jeta, jphi);
Double_t hits = rh->fBasicHits->GetBinContent(ieta,iphi);
+
// Mean in region of interest
Double_t poissonM = (total <= 0 || empty <= 0 ? 0 :
-TMath::Log(empty / total));
+
+ // Note, that given filled=total-empty, and
+ //
+ // m = -log(empty/total)
+ // = -log(1 - filled/total)
+ //
+ // v = m / (1 - exp(-m))
+ // = -total/filled * (log(total-filled)-log(total))
+ // = -total / (total-empty) * log(empty/total)
+ // = total (log(total)-log(empty)) / (total-empty)
+ //
Double_t poissonV = hits;
if(poissonM > 0)
// Correct for counting statistics and weight by counts
- poissonV = (hits * poissonM) / (1 - TMath::Exp(-1*poissonM));
+ poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
Double_t poissonE = TMath::Sqrt(hits);
if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
if (fUsePoisson) {
h->SetBinContent(ieta,iphi,poissonV);
h->SetBinError(ieta,iphi,poissonE);
+ rh->fDensity->Fill(eta, phi, poissonV);
}
}
}
fFMD3iMax.Set(nEta);
fFMD3oMax.Set(nEta);
+ fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
+ fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
+ fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
+ fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
+ fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
+ fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
+
+ AliInfo(Form("Get eta axis with %d bins from %f to %f",
+ nEta, eta.GetXmin(), eta.GetXmax()));
+ fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
+ fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
+ fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
+ fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
+ fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
+ fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
+
for (Int_t i = 0; i < nEta; i++) {
- fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
- fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
- fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
- fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
- fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
+ Double_t w[5];
+ w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
+ w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
+ w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
+ w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
+ w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
+ Double_t l[5];
+ l[0] = GetMultCut(1, 'I', i+1, false);
+ l[1] = GetMultCut(2, 'I', i+1, false);
+ l[2] = GetMultCut(2, 'O', i+1, false);
+ l[3] = GetMultCut(3, 'I', i+1, false);
+ l[4] = GetMultCut(3, 'O', i+1, false);
+ for (Int_t j = 0; j < 5; j++) {
+ if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
+ if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
+ }
}
}
// Return:
// The number of particles
//
- if (mult <= GetMultCut()) return 0;
+ // if (mult <= GetMultCut()) return 0;
if (lowFlux) return 1;
AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
//
AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
- Float_t correction = AcceptanceCorrection(r,t);
+ Float_t correction = 1;
+ if (fUsePhiAcceptance) correction = AcceptanceCorrection(r,t);
if (lowFlux) {
TH1D* dblHitCor = 0;
if (fcm.GetDoubleHit())
TIter next(&fRingHistos);
RingHistos* o = 0;
- while ((o = static_cast<RingHistos*>(next())))
+ THStack* sums = new THStack("sums", "sums of ring signals");
+ while ((o = static_cast<RingHistos*>(next()))) {
o->ScaleHistograms(d, nEvents);
+ TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
+ o->fDensity->GetNbinsY(),"e");
+ sum->Scale(1., "width");
+ sum->SetTitle(o->GetName());
+ sum->SetDirectory(0);
+ sum->SetYTitle("#sum N_{ch,incl}");
+ sums->Add(sum);
+ }
+ d->Add(sums);
}
//____________________________________________________________________
// dir List to write in
//
TList* d = new TList;
+ d->SetOwner();
d->SetName(GetName());
dir->Add(d);
d->Add(fWeightedSum);
d->Add(fCorrections);
d->Add(fAccI);
d->Add(fAccO);
+ d->Add(fMaxWeights);
+ d->Add(fLowCuts);
TIter next(&fRingHistos);
RingHistos* o = 0;
for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
ind[gROOT->GetDirLevel()] = '\0';
std::cout << ind << ClassName() << ": " << GetName() << '\n'
+ << std::boolalpha
<< ind << " Multiplicity cut: " << fMultCut << '\n'
+ << ind << " # of (xi+sigma) factor: " << fNXi << '\n'
+ << ind << " Include sigma in cut: " << fIncludeSigma << '\n'
+ << ind << " Low cut method: "
+ << (fMultCut > 0 ? "fixed" :
+ (fNXi >= 0 ? "xi+sigma" : "fit range")) << '\n'
<< ind << " Max(particles): " << fMaxParticles << '\n'
+ << ind << " Poisson method: " << fUsePoisson << '\n'
+ << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n'
<< ind << " Eta lumping: " << fEtaLumping << '\n'
<< ind << " Phi lumping: " << fPhiLumping << '\n'
+ << std::noboolalpha
<< std::flush;
TString opt(option);
opt.ToLower();
// d detector
// r ring
//
- fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()),
- Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
- "N_{ch} in %s", fName.Data()),
- 2500, -.5, 24.5, 2500, -.5, 24.5);
- fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()),
- Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
- "N_{ch} in %s", fName.Data()),
- 2500, -.5, 24.5, 2500, -.5, 24.5);
+ fEvsN = new TH2D("elossVsNnocorr",
+ "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
+ 250, -.5, 24.5, 250, -.5, 24.5);
fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
fEvsN->Sumw2();
fEvsN->SetDirectory(0);
- fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
- fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
- fEvsM->Sumw2();
+
+ fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
+ fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
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 = new TProfile("etaVsNnocorr",
+ "Average inclusive N_{ch} vs #eta (uncorrected)",
+ 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 = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
+ fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
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 = new TProfile("corr", "Average correction", 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} density in %s", fName.Data()),
+ fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
0, 2*TMath::Pi());
fDensity->SetDirectory(0);
+ fDensity->Sumw2();
+ fDensity->SetMarkerColor(Color());
fDensity->SetXTitle("#eta");
fDensity->SetYTitle("#phi [radians]");
fDensity->SetZTitle("Inclusive N_{ch} density");
- fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()),
- Form("N_{ch} from energy loss vs from Poission %s",
- fName.Data()), 100, 0, 20, 100, 0, 20);
+ fELossVsPoisson = new TH2D("elossVsPoisson",
+ "N_{ch} from energy loss vs from Poission",
+ 100, 0, 20, 100, 0, 20);
fELossVsPoisson->SetDirectory(0);
fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
fELossVsPoisson->SetZTitle("Correlation");
- fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data()),
- Form("# of empty strips vs. total @ # strips in %s",
- fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5);
+ fEmptyVsTotal = new TH2D("emptyVsTotal",
+ "# of empty strips vs. total # strips",
+ 21, -.5, 20.5, 21, -0.5, 20.5);
fEmptyVsTotal->SetDirectory(0);
fEmptyVsTotal->SetXTitle("Total # strips");
fEmptyVsTotal->SetYTitle("Empty # strips");
TList* l = GetOutputList(dir);
if (!l) return;
- TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
+ TH1* density = GetOutputHist(l,"inclDensity");
if (density) density->Scale(1./nEvents);
}
* number of particles that has hit within a region.
*/
void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
+ /**
+ * Set whether to use the phi acceptance correction.
+ *
+ * @param u If true, use the phi acceptance (default is false)
+ */
+ void SetUsePhiAcceptance(Bool_t u) { fUsePhiAcceptance = u; }
/**
* 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; }
+ void SetMultCut(Double_t cut) { SetMultCuts(cut,cut,cut,cut,cut); }
+ /**
+ * Set the lower multiplicity cuts
+ *
+ * @param fmd1i Lower mulitplicyt cut for FMD1i
+ * @param fmd2i Lower mulitplicyt cut for FMD2i
+ * @param fmd2o Lower mulitplicyt cut for FMD2o
+ * @param fmd3i Lower mulitplicyt cut for FMD3i
+ * @param fmd3o Lower mulitplicyt cut for FMD3o
+ */
+ void SetMultCuts(Double_t fmd1i,
+ Double_t fmd2i,
+ Double_t fmd2o,
+ Double_t fmd3i,
+ Double_t fmd3o);
/**
* Set the luming factors used in the Poisson method
*
fEtaLumping = (eta < 1 ? 1 : eta);
fPhiLumping = (phi < 1 ? 1 : phi);
}
+ /**
+ * Set the number of landau width to subtract from the most probably
+ * value to get the low cut.
+ *
+ * @param n Number of @f$ \xi@f$
+ */
+ void SetNXi(Double_t nXi) { fNXi = nXi; }
+ /**
+ * Whether to include sigma in the number subtracted from the MPV to
+ * get the low cut
+ *
+ * @param u If true, then low cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$
+ */
+ void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; }
+ /**
+ * 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(UShort_t d, Char_t r, Double_t eta,
+ Bool_t errors=true) const;
/**
* Get the multiplicity cut. If the user has set fMultCut (via
* SetMultCut) then that value is used. If not, then the lower
*
* @return Lower cut on multiplicity
*/
- Double_t GetMultCut() const;
+ Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
+ Bool_t errors=true) const;
/**
* Print information
*
TList fRingHistos; // List of histogram containers
Double_t fMultCut; // Low cut on scaled energy loss
+ Double_t fNXi; // Delta-n(xi+sigma) factor
+ Bool_t fIncludeSigma; // Whether or not to include sigma in cut
TH1D* fSumOfWeights; // Histogram
TH1D* fWeightedSum; // Histogram
TH1D* fCorrections; // Histogram
UShort_t fMaxParticles; // Maximum particle weight to use
Bool_t fUsePoisson; // If true, then use poisson statistics
+ Bool_t fUsePhiAcceptance; // Whether to correct for corners
TH1D* fAccI; // Acceptance correction for inner rings
TH1D* fAccO; // Acceptance correction for outer rings
TArrayI fFMD1iMax; // Array of max weights
TArrayI fFMD2oMax; // Array of max weights
TArrayI fFMD3iMax; // Array of max weights
TArrayI fFMD3oMax; // Array of max weights
+ TH2D* fMaxWeights; // Histogram of max weights
+ TH2D* fLowCuts; // Histogram of low cuts
Int_t fEtaLumping; // How to lump eta bins for Poisson
Int_t fPhiLumping; // How to lump phi bins for Poisson
Int_t fDebug; // Debug level
+ Double_t fMultCuts[5]; // Per-ring multiplicity cuts
-
- ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density
+ ClassDef(AliFMDDensityCalculator,3); // Calculate Nch density
};
#endif
AliMCEvent* mcevent = MCEvent();
if(mcevent) {
AliHeader* header = mcevent->Header();
- AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
+ AliGenHijingEventHeader* hijingHeader =
+ dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
if(hijingHeader) {
Float_t b = hijingHeader->ImpactParameter();
if(b<fbLow || b>fbHigh) return;
InitializeSubs();
}
- Bool_t lowFlux = kFALSE;
- UInt_t triggers = 0;
- UShort_t ivz = 0;
- Double_t vz = 0;
- Double_t cent = 0;
- UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent);
+ Bool_t lowFlux = kFALSE;
+ UInt_t triggers = 0;
+ UShort_t ivz = 0;
+ Double_t vz = 0;
+ Double_t cent = 0;
+ UShort_t nClusters = 0;
+ UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
+ ivz, vz, cent, nClusters);
if (found & AliFMDEventInspector::kNoEvent) return;
if (found & AliFMDEventInspector::kNoTriggers) return;
if (found & AliFMDEventInspector::kNoSPD) return;
return;
}
-#if 0
- // 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;
- }
-#endif
AliInfo("Fitting energy loss spectra");
fEnergyFitter.Fit(list);
TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
list->GetName())));
PostData(2, list2);
-#if 0
- // Temporary code to save output to a file - should be disabled once
- // the plugin stuff works
- list->ls();
- TDirectory* savdir = gDirectory;
- TFile* tmp = TFile::Open("elossfits.root", "RECREATE");
- list->Write(list->GetName(), TObject::kSingleKey);
- tmp->Write();
- tmp->Close();
- savdir->cd();
-#endif
}
//____________________________________________________________________
: TNamed(),
fHEventsTr(0),
fHEventsTrVtx(0),
+ fHEventsAccepted(0),
fHTriggers(0),
fHType(0),
fHWords(0),
fField(999),
fCollisionSystem(kUnknown),
fDebug(0),
- fCentAxis(0)
+ fCentAxis(0),
+ fVtxAxis(10,-10,10)
{
//
// Constructor
: TNamed("fmdEventInspector", name),
fHEventsTr(0),
fHEventsTrVtx(0),
+ fHEventsAccepted(0),
fHTriggers(0),
fHType(0),
fHWords(0),
fField(999),
fCollisionSystem(kUnknown),
fDebug(0),
- fCentAxis(0)
+ fCentAxis(0),
+ fVtxAxis(10,-10,10)
{
//
// Constructor
: TNamed(o),
fHEventsTr(o.fHEventsTr),
fHEventsTrVtx(o.fHEventsTrVtx),
+ fHEventsAccepted(o.fHEventsAccepted),
fHTriggers(o.fHTriggers),
fHType(o.fHType),
fHWords(o.fHWords),
fField(o.fField),
fCollisionSystem(o.fCollisionSystem),
fDebug(0),
- fCentAxis(0)
+ fCentAxis(0),
+ fVtxAxis(o.fVtxAxis)
{
//
// Copy constructor
TNamed::operator=(o);
fHEventsTr = o.fHEventsTr;
fHEventsTrVtx = o.fHEventsTrVtx;
+ fHEventsAccepted = o.fHEventsAccepted;
fHTriggers = o.fHTriggers;
fHType = o.fHType;
fHWords = o.fHWords;
fEnergy = o.fEnergy;
fField = o.fField;
fCollisionSystem = o.fCollisionSystem;
+ fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(),
+ o.fVtxAxis.GetXmax());
if (fList) {
fList->SetName(GetName());
if (fHEventsTr) fList->Add(fHEventsTr);
// ----- 92 number --------- ---- 1 ---
TArrayD limits(93);
for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
+
+ fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
fCentAxis = new TAxis(limits.GetSize()-1, limits.GetArray());
fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
- vtxAxis.GetNbins(),
- vtxAxis.GetXmin(),
- vtxAxis.GetXmax());
+ 4*vtxAxis.GetNbins(),
+ 2*vtxAxis.GetXmin(),
+ 2*vtxAxis.GetXmax());
fHEventsTr->SetXTitle("v_{z} [cm]");
fHEventsTr->SetYTitle("# of events");
fHEventsTr->SetFillColor(kRed+1);
// fHEventsTr->Sumw2();
fList->Add(fHEventsTr);
- fHEventsTrVtx = new TH1I("nEventsTrVtx",
- "Number of events w/trigger and vertex",
- vtxAxis.GetNbins(),
- vtxAxis.GetXmin(),
- vtxAxis.GetXmax());
- fHEventsTrVtx->SetXTitle("v_{z} [cm]");
- fHEventsTrVtx->SetYTitle("# of events");
+ fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
+ fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex");
fHEventsTrVtx->SetFillColor(kBlue+1);
- fHEventsTrVtx->SetFillStyle(3001);
fHEventsTrVtx->SetDirectory(0);
// fHEventsTrVtx->Sumw2();
fList->Add(fHEventsTrVtx);
+ fHEventsAccepted = new TH1I("nEventsAccepted",
+ "Number of events w/trigger and vertex in range",
+ 2*vtxAxis.GetNbins(),
+ 2*vtxAxis.GetXmin(),
+ 2*vtxAxis.GetXmax());
+ fHEventsAccepted->SetXTitle("v_{z} [cm]");
+ fHEventsAccepted->SetYTitle("# of events");
+ fHEventsAccepted->SetFillColor(kGreen+1);
+ fHEventsAccepted->SetFillStyle(3001);
+ fHEventsAccepted->SetDirectory(0);
+ // fHEventsAccepted->Sumw2();
+ fList->Add(fHEventsAccepted);
+
fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
fHTriggers->SetFillColor(kRed+1);
Bool_t& lowFlux,
UShort_t& ivz,
Double_t& vz,
- Double_t& cent)
+ Double_t& cent,
+ UShort_t& nClusters)
{
//
// Process the event
}
// --- Read trigger information from the ESD and store in AOD object
- if (!ReadTriggers(event, triggers)) {
+ if (!ReadTriggers(event, triggers, nClusters)) {
if (fDebug > 2) {
AliWarning("Failed to read triggers from ESD"); }
return kNoTriggers;
fHEventsTrVtx->Fill(vz);
// --- Get the vertex bin ------------------------------------------
- ivz = fHEventsTr->GetXaxis()->FindBin(vz);
- if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
+ ivz = fVtxAxis.FindBin(vz);
+ if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
if (fDebug > 3) {
AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
- vz, fHEventsTr->GetXaxis()->GetXmin(),
- fHEventsTr->GetXaxis()->GetXmax())); }
+ vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
+ }
ivz = 0;
return kBadVertex;
}
+ fHEventsAccepted->Fill(vz);
// --- Check the FMD ESD data --------------------------------------
if (!event->GetFMDData()) {
//____________________________________________________________________
Bool_t
-AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
+AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers,
+ UShort_t& nClusters)
{
//
// Read the trigger information from the ESD event
// then the AliPhysicsSelection object would overcount by a
// factor of 2! :-(
Bool_t offline = ih->IsEventSelected();
+ nClusters = 0;
if (offline) {
triggers |= AliAODForwardMult::kOffline;
triggers |= AliAODForwardMult::kInel;
// Check if we have one or more tracklets
// in the range -1 < eta < 1 to set the INEL>0
// trigger flag.
+ //
+ // Also count tracklets as a single cluster
Int_t n = spdmult->GetNumberOfTracklets();
for (Int_t j = 0; j < n; j++) {
if(TMath::Abs(spdmult->GetEta(j)) < 1) {
triggers |= AliAODForwardMult::kInelGt0;
- break;
+ nClusters++;
}
}
+ n = spdmult->GetNumberOfSingleClusters();
+ for (Int_t j = 0; j < n; j++) {
+ Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
+ if (TMath::Abs(eta) < 1) nClusters++;
+ }
}
+ if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
}
// Analyse some trigger stuff
triggers |= AliAODForwardMult::kNSD;
- //Check pileup
+ // Check for multiple vertices (pile-up) with at least 3
+ // contributors and at least 0.8cm from the primary vertex
Bool_t pileup = esd->IsPileupFromSPD(3,0.8);
if (pileup) {
triggers |= AliAODForwardMult::kPileUp;
// @c true on success, @c false otherwise
//
vz = 0;
+#if 1
+ // This is the code used by the 1st physics people
+ const AliESDVertex* vertex = esd->GetPrimaryVertex();
+ if (!vertex || !vertex->GetStatus()) {
+ if (fDebug > 2) {
+ AliWarning(Form("No primary vertex (%p) or bad status %d",
+ vertex, (vertex ? vertex->GetStatus() : -1)));
+ }
+ return false;
+ }
+ const AliESDVertex* vertexSPD = esd->GetPrimaryVertexSPD();
+ if (!vertexSPD || !vertexSPD->GetStatus()) {
+ if (fDebug > 2) {
+ AliWarning(Form("No primary SPD vertex (%p) or bad status %d",
+ vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1)));
+ }
+ return false;
+ }
+
+ // if vertex is from SPD vertexZ, require more stringent cuts
+ if (vertex->IsFromVertexerZ()) {
+ if (vertex->GetDispersion() > fMaxVzErr ||
+ vertex->GetZRes() > 1.25 * fMaxVzErr) {
+ if (fDebug > 2) {
+ AliWarning(Form("Dispersion %f > %f or resolution %f > %f",
+ vertex->GetDispersion(), fMaxVzErr,
+ vertex->GetZRes(), 1.25 * fMaxVzErr));
+ }
+ return false;
+ }
+ }
+ vz = vertex->GetZ();
+ return true;
+#else
// Get the vertex
const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
if (!vertex) {
// Get the z coordiante
vz = vertex->GetZ();
return kTRUE;
+#endif
}
//____________________________________________________________________
* @ingroup pwg2_forward_aod
*/
#include <TNamed.h>
+#include <TAxis.h>
class AliESDEvent;
class TH2D;
class TH1D;
* @param vz On return, the z position of the interaction
* @param cent On return, the centrality (in percent) or < 0
* if not found
+ * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$
*
* @return 0 (or kOk) on success, otherwise a bit mask of error codes
*/
Bool_t& lowFlux,
UShort_t& ivz,
Double_t& vz,
- Double_t& cent);
+ Double_t& cent,
+ UShort_t& nClusters);
/**
* Define the output histograms. These are put in a sub list of the
* passed list. The histograms are merged before the parent task calls
*
* @param esd ESD event
* @param triggers On return, contains the trigger bits
+ * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$
*
* @return @c true on success, @c false otherwise
*/
- Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers);
+ Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers,
+ UShort_t& nClusters);
/**
* Read the vertex information from the ESD event
*
TH1I* fHEventsTr; //! Histogram of events w/trigger
TH1I* fHEventsTrVtx; //! Events w/trigger and vertex
+ TH1I* fHEventsAccepted; //! Events w/trigger and vertex in range
TH1I* fHTriggers; //! Triggers
TH1I* fHType; //! Type (low/high flux) of event
TH1I* fHWords; //! Trigger words
UShort_t fCollisionSystem; // Collision system
Int_t fDebug; // Debug level
TAxis* fCentAxis; // Centrality axis used in histograms
+ TAxis fVtxAxis;
ClassDef(AliFMDEventInspector,3); // Inspect the event
};
; // For Emacs
#endif
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector()
+ : fNCutBins(0),
+ fCorrectionCut(0),
+ fFirstBins(),
+ fLastBins(),
+ fDebug(0),
+ fList(0),
+ fSumRings(0),
+ fCoverage(0),
+ fMergeMethod(kStraightMean),
+ fFiducialMethod(kByCut)
+{}
+
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const char* title)
+ : TNamed("fmdHistCollector", title),
+ fNCutBins(2),
+ fCorrectionCut(0.5),
+ fFirstBins(1),
+ fLastBins(1),
+ fDebug(0),
+ fList(0),
+ fSumRings(0),
+ fCoverage(0),
+ fMergeMethod(kStraightMean),
+ fFiducialMethod(kByCut)
+{
+}
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
+ : TNamed(o),
+ fNCutBins(o.fNCutBins),
+ fCorrectionCut(o.fCorrectionCut),
+ fFirstBins(o.fFirstBins),
+ fLastBins(o.fLastBins),
+ fDebug(o.fDebug),
+ fList(o.fList),
+ fSumRings(o.fSumRings),
+ fCoverage(o.fCoverage),
+ fMergeMethod(o.fMergeMethod),
+ fFiducialMethod(o.fFiducialMethod)
+{}
+
//____________________________________________________________________
AliFMDHistCollector::~AliFMDHistCollector()
{
fCoverage = o.fCoverage;
fMergeMethod = o.fMergeMethod;
fFiducialMethod = o.fFiducialMethod;
+
return *this;
}
// Store the result for later use
fFirstBins[(iVz-1)*5+iIdx] = first;
fLastBins[(iVz-1)*5+iIdx] = last;
- TH2D* obg = static_cast<TH2D*>(bg->Clone());
+ TH2D* obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
obg->SetDirectory(0);
obg->Reset();
vtxList->Add(obg);
fList->SetOwner();
fList->SetName(GetName());
dir->Add(fList);
+
}
+
//____________________________________________________________________
Int_t
AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
//____________________________________________________________________
Bool_t
-AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
+AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
+ AliForwardUtil::Histos& sums,
UShort_t vtxbin,
TH2D& out)
{
TH2D* h = hists.Get(d,r);
TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
Int_t i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
+ TH2D* o = sums.Get(d, r);
- // Outer rings have better phi segmentation - rebin to same as inner.
- if (q == 1) t->RebinY(2);
-
+ // Get valid range
Int_t first = 0;
Int_t last = 0;
GetFirstAndLast(d, r, vtxbin, first, last);
+
+ // Zero outside valid range
+ for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) {
+ // Lower range
+ for (Int_t iEta = 1; iEta < first; iEta++) {
+ t->SetBinContent(iEta,iPhi,0);
+ t->SetBinError(iEta,iPhi,0);
+ }
+ for (Int_t iEta = last+1; iEta <= t->GetNbinsX(); iEta++) {
+ t->SetBinContent(iEta,iPhi,0);
+ t->SetBinError(iEta,iPhi,0);
+ }
+ }
+ for (Int_t iEta = first; iEta <= last; iEta++)
+ t->SetBinContent(iEta,0,1);
+ // Add to our per-ring sum
+ o->Add(t);
+
+ // Outer rings have better phi segmentation - rebin to same as inner.
+ if (q == 1) t->RebinY(2);
// Now update profile output
for (Int_t iEta = first; iEta <= last; iEta++) {
Float_t noc = overlap >= 0? 0.5 : 1;
out.SetBinContent(iEta, 0, ooc + noc);
- for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
+ // Should we loop over h or t Y bins - I think it's t
+ for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) {
Double_t c = t->GetBinContent(iEta,iPhi);
Double_t e = t->GetBinError(iEta,iPhi);
Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
case kLeastError: std::cout << "least error\n"; break;
}
- std::cout << ind << " Bin ranges:\n" << ind << " v_z bin";
+ std::cout << ind << " Bin ranges:\n" << ind << " rings ";
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 << " ";
+ }
+ std::cout << '\n' << ind << " /vz_bin ";
+ for (Int_t iIdx = 0; iIdx < 5; iIdx++)
+ std::cout << "-+--------";
+ std::cout << std::endl;
+
+ for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
+ std::cout << " " << std::setw(7) << iVz << " ";
+ 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 << "-"
/**
* Constructor
*/
- AliFMDHistCollector()
- : fNCutBins(0),
- fCorrectionCut(0),
- fFirstBins(),
- fLastBins(),
- fDebug(0),
- fList(0),
- fSumRings(0),
- fCoverage(0),
- fMergeMethod(kStraightMean),
- fFiducialMethod(kByCut)
- {}
+ AliFMDHistCollector();
/**
* Constructor
*
* @param title Name of object
*/
- AliFMDHistCollector(const char* title)
- : TNamed("fmdHistCollector", title),
- fNCutBins(2),
- fCorrectionCut(0.5),
- fFirstBins(1),
- fLastBins(1),
- fDebug(0),
- fList(0),
- fSumRings(0),
- fCoverage(0),
- fMergeMethod(kStraightMean),
- fFiducialMethod(kByCut)
- {}
+ AliFMDHistCollector(const char* title);
/**
* Copy constructor
*
* @param o Object to copy from
*/
- AliFMDHistCollector(const AliFMDHistCollector& o)
- : TNamed(o),
- fNCutBins(o.fNCutBins),
- fCorrectionCut(o.fCorrectionCut),
- fFirstBins(o.fFirstBins),
- fLastBins(o.fLastBins),
- fDebug(o.fDebug),
- fList(o.fList),
- fSumRings(o.fSumRings),
- fCoverage(o.fCoverage),
- fMergeMethod(o.fMergeMethod),
- fFiducialMethod(o.fFiducialMethod)
- {}
+ AliFMDHistCollector(const AliFMDHistCollector& o);
/**
* Destructor
* Do the calculations
*
* @param hists Cache of histograms
+ * @param sums Cache to sum ring histograms in
* @param vtxBin Vertex bin (1 based)
* @param out Output histogram
*
* @return true on successs
*/
- virtual Bool_t Collect(AliForwardUtil::Histos& hists, UShort_t vtxBin,
- TH2D& out);
+ virtual Bool_t Collect(const AliForwardUtil::Histos& hists,
+ AliForwardUtil::Histos& sums,
+ UShort_t vtxBin,
+ TH2D& out);
/**
* Output diagnostic histograms to directory
*
void MergeBins(Double_t c, Double_t e,
Double_t oc, Double_t oe,
Double_t& rc, Double_t& re) const;
+
Int_t fNCutBins; // Number of additional bins to cut away
Float_t fCorrectionCut; // Cut-off on secondary corrections
TH2D* fCoverage; // Sum per ring (on y-axis)
MergeMethod fMergeMethod; // Merge methiod for overlapping bins
FiducialMethod fFiducialMethod; // Fidicual method
+
ClassDef(AliFMDHistCollector,1); // Calculate Nch density
};
// Reference to this object
//
AliFMDCorrector::operator=(o);
-
+ fSecondaryForMC = o.fSecondaryForMC;
return *this;
}
// Return:
// true on successs
//
+ if ((!fUseSecondaryMap || !fSecondaryForMC) && fUseVertexBias)
+ return kTRUE;
+
AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
UShort_t uvb = vtxbin;
TH2D* h = hists.Get(d,r);
- if (fUseSecondaryMap) {
+ if (fUseSecondaryMap && fSecondaryForMC) {
TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
if (!bg) {
AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
fComps->SetName("esd_mc_comparison");
d->Add(fComps);
}
+//____________________________________________________________________
+void
+AliFMDMCCorrector::Print(Option_t* option) const
+{
+ //
+ // Print information
+ // Parameters:
+ // option Not used
+ //
+ char ind[gROOT->GetDirLevel()+1];
+ for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+ ind[gROOT->GetDirLevel()] = '\0';
+ AliFMDCorrector::Print(option);
+ std::cout << std::boolalpha
+ << ind << " Use sec. map on MC: " << fSecondaryForMC
+ << std::noboolalpha << std::endl;
+}
//____________________________________________________________________
//
fFMD2o(0),
fFMD3i(0),
fFMD3o(0),
- fComps(0)
+ fComps(0),
+ fSecondaryForMC(true)
{}
/**
* Constructor
fFMD2o(0),
fFMD3i(0),
fFMD3o(0),
- fComps(0)
+ fComps(0),
+ fSecondaryForMC(true)
{}
/**
* Copy constructor
fFMD2o(o.fFMD2o),
fFMD3i(o.fFMD3i),
fFMD3o(o.fFMD3o),
- fComps(0)
+ fComps(0),
+ fSecondaryForMC(o.fSecondaryForMC)
{}
/**
* Destructor
* @return Reference to this object
*/
AliFMDMCCorrector& operator=(const AliFMDMCCorrector&);
+ /**
+ * If set, then do not do the secondary correction for MC data
+ *
+ * @param use
+ */
+ void SetSecondaryForMC(Bool_t use) { fSecondaryForMC = use; }
/**
* Initialize this object
*
* @param dir List to write in
*/
void DefineOutput(TList* dir);
+
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
+ void Print(Option_t* option="") const;
protected:
/**
* MAke comparison profiles
TProfile2D* fFMD3i; // Comparison
TProfile2D* fFMD3o; // Comparison
TList* fComps; // List of comparisons
-
+ Bool_t fSecondaryForMC; // Whether to correct MC data
+
ClassDef(AliFMDMCCorrector,1); // Calculate Nch density
};
else
return false;
- // Invariant masses
+ // Rapidity shift
Double_t m02s = 1 - 2 * p1->Energy() / fEnergy;
Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
// - AliESDFMD object - copy of input, but with signals merged
//
// Corrections used:
-// - None
+// - ELoss fits
//
// Histograms:
// - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
#include <TH1.h>
#include <TMath.h>
#include "AliFMDStripIndex.h"
+#include "AliFMDFloatMap.h"
#include <AliLog.h>
#include <TROOT.h>
#include <iostream>
fFMD2o(0),
fFMD3i(0),
fFMD3o(0),
- fSumEta(0)
+ fSumEta(0),
+ fOperComp(0),
+ fThetaVsNr(0),
+ fOnlyPrimary(false)
{
//
// Constructor
fSumEta->SetMarkerStyle(22);
fSumEta->SetFillColor(0);
fSumEta->SetFillStyle(0);
+
+ fOper = new AliFMDFloatMap(0,0,0,0);
+ fOperComp = new TH2I("operComp", "Operation vs # track refs",
+ kMergedInto, kNone-.5, kMergedInto+.5,
+ 20, -.5, 19.5);
+ fOperComp->SetXTitle("Operation");
+ fOperComp->SetYTitle("# of track refs in sector");
+ fOperComp->SetZTitle("Observations");
+ fOperComp->GetXaxis()->SetBinLabel(kNone, "None");
+ fOperComp->GetXaxis()->SetBinLabel(kCandidate, "Candidate");
+ fOperComp->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
+ fOperComp->GetXaxis()->SetBinLabel(kMergedInto, "Merged into");
+ fOperComp->SetDirectory(0);
+ fThetaVsNr = new TH2D("thetaVsNr", "#theta of track vs # track references",
+ 360, 0, 360, 20, -.5, 19.5);
+ fThetaVsNr->SetXTitle("#theta [degrees]");
+ fThetaVsNr->SetYTitle("# of track references");
+ fThetaVsNr->SetDirectory(0);
}
//____________________________________________________________________
fFMD2o(o.fFMD2o),
fFMD3i(o.fFMD3i),
fFMD3o(o.fFMD3o),
- fSumEta(o.fSumEta)
+ fSumEta(o.fSumEta),
+ fOperComp(o.fOperComp),
+ fThetaVsNr(o.fThetaVsNr),
+ fOnlyPrimary(o.fOnlyPrimary)
{
//
// Copy constructor
// Reference to this
//
AliFMDSharingFilter::operator=(o);
+ fOnlyPrimary = o.fOnlyPrimary;
return *this;
}
Char_t r,
UShort_t s,
UShort_t t,
+ UShort_t nR,
+ Double_t theta,
AliESDFMD& output) const
{
//
// r Ring
// s Sector
// t Strip
+ // nR Number of references to this particle in this sector
// output Output ESD object
//
Double_t old = output.Multiplicity(d,r,s,t);
if (old == AliESDFMD::kInvalidMult) old = 0;
+ if (fOper) fOperComp->Fill(fOper->operator()(d,r,s,t), nR);
+ if (theta < 0) theta += 2*TMath::Pi();
+ theta *= 180. / TMath::Pi();
+ fThetaVsNr->Fill(theta, nR);
output.SetMultiplicity(d,r,s,t,old+1);
}
//
output.Clear();
+ // Increase event count - stored in
+ // underflow bin
+ fSumEta->AddBinContent(0);
+
// Copy eta values to output
for (UShort_t ed = 1; ed <= 3; ed++) {
UShort_t nq = (ed == 1 ? 1 : 2);
if (!isCharged) continue;
Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
+ // Pseudo rapidity and azimuthal angle
+ Double_t eta = particle->Eta();
+ Double_t phi = particle->Phi();
+
// Fill 'dn/deta' histogram
if (isPrimary && iTr < nPrim) {
- fSumEta->Fill(particle->Eta());
- primary->Fill(particle->Eta(), particle->Phi());
+ // Avoid under count - used to store event count
+ if (eta >= fSumEta->GetXaxis()->GetXmin()) fSumEta->Fill(eta);
+ primary->Fill(eta, phi);
}
+ // Bail out if we're only processing primaries - perhaps we should
+ // track back to the original primary?
+ if (fOnlyPrimary && !isPrimary) continue;
+
Int_t nTrRef = particle->GetNumberOfTrackReferences();
Int_t longest = -1;
Double_t angle = 0;
UShort_t oD = 0, oS = 1024, oT = 1024;
Char_t oR = '\0';
+ UShort_t nC = 0;
+ Double_t oTheta = 0;
for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
AliTrackReference* ref = particle->GetTrackReference(iTrRef);
if (ref->DetectorId() != AliTrackReference::kFMD)
continue;
+ // Count number of track refs in this sector
+ nC++;
+
// Get the detector coordinates
UShort_t d, s, t;
Char_t r;
AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
- if (oD > 0 && oR != '\0' && (d != oD || r != oR)) {
+ // If this is a new detector/ring, then reset the other one
+ if (oD > 0 && oR != '\0' && oS != 1024 &&
+ (d != oD || r != oR || s != oS)) {
longest = -1;
- StoreParticle(oD, oR, oS, oT, output);
+ angle = 0;
+ StoreParticle(oD, oR, oS, oT, nC, oTheta, output);
+ nC = 0;
+ oD = 0;
+ oR = '\0';
+ oS = 1024;
+ oT = 1024;
}
oD = d;
longest = iTrRef;
angle = ang;
}
+ oTheta = theta;
} // Loop over track references
if (longest < 0) continue;
Char_t r;
AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
- StoreParticle(d,r,s,t,output);
+ StoreParticle(d,r,s,t,nC,particle->Theta(),output);
} // Loop over tracks
return kTRUE;
}
cd->Add(fFMD3i);
cd->Add(fFMD3o);
dir->Add(fSumEta);
+ cd->Add(fOperComp);
+ cd->Add(fThetaVsNr);
}
//____________________________________________________________________
AliWarning(Form("No mcSumEta histogram found in output list"));
return;
}
- sumEta->Scale(1. / nEvents, "width");
+ Double_t n = nEvents; // sumEta->GetBinContent(0);
+ sumEta->Scale(1. / n, "width");
}
//____________________________________________________________________
for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
ind[gROOT->GetDirLevel()] = '\0';
AliFMDSharingFilter::Print(option);
+ std::cout << std::boolalpha
+ << ind << " Only primary tracks: " << fOnlyPrimary
+ << std::noboolalpha << std::endl;
}
//____________________________________________________________________
fFMD2o(0),
fFMD3i(0),
fFMD3o(0),
- fSumEta(0)
+ fSumEta(0),
+ fOperComp(0),
+ fThetaVsNr(0),
+ fOnlyPrimary(false)
{}
/**
* Constructor
* @return Reference to this
*/
AliFMDMCSharingFilter& operator=(const AliFMDMCSharingFilter& o);
+
+ /**
+ * If set, then only process primary tracks
+ *
+ * @param use
+ */
+ void SetOnlyPrimary(Bool_t use) { fOnlyPrimary = use; }
/**
* Filter the input kinematics and track references, using
* some of the ESD information
* @param output Output ESD object
*/
void StoreParticle(UShort_t d, Char_t r, UShort_t s, UShort_t t,
- AliESDFMD& output) const;
+ UShort_t nr, Double_t theta, AliESDFMD& output) const;
TH2D* fFMD1i; // ESD-MC correlation
TH2D* fFMD2i; // ESD-MC correlation
TH2D* fFMD2o; // ESD-MC correlation
TH2D* fFMD3i; // ESD-MC correlation
TH2D* fFMD3o; // ESD-MC correlation
TH1D* fSumEta; // MC dN/deta
-
+ TH2I* fOperComp; // Operation vs # trackrefs
+ TH2D* fThetaVsNr; // Theta vs # trackrefs
+ Bool_t fOnlyPrimary; // Only process primary tracks
ClassDef(AliFMDMCSharingFilter,1); //
};
#include "AliFMDCorrELossFit.h"
#include <AliLog.h>
#include <TROOT.h>
+#include <THStack.h>
#include <iostream>
#include <iomanip>
; // This is for Emacs
#endif
+#define DBG(L,M) \
+ do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false)
+#define DBGL(L,M) \
+ do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false)
+
+
//____________________________________________________________________
AliFMDSharingFilter::AliFMDSharingFilter()
fLowCut(0.),
fCorrectAngles(kFALSE),
fNXi(1),
+ fIncludeSigma(true),
+ fSummed(0),
+ fHighCuts(0),
+ fOper(0),
fDebug(0)
{
//
fLowCut(0.),
fCorrectAngles(kFALSE),
fNXi(1),
+ fIncludeSigma(true),
+ fSummed(0),
+ fHighCuts(0),
+ fOper(0),
fDebug(0)
{
//
fLowCut(o.fLowCut),
fCorrectAngles(o.fCorrectAngles),
fNXi(o.fNXi),
+ fIncludeSigma(o.fIncludeSigma),
+ fSummed(o.fSummed),
+ fHighCuts(o.fHighCuts),
+ fOper(o.fOper),
fDebug(o.fDebug)
{
//
fCorrectAngles = o.fCorrectAngles;
fNXi = o.fNXi;
fDebug = o.fDebug;
+ fOper = o.fOper;
+ fSummed = o.fSummed;
+ fHighCuts = o.fHighCuts;
+ fIncludeSigma = o.fIncludeSigma;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
+//____________________________________________________________________
+void
+AliFMDSharingFilter::Init()
+{
+ // Initialise
+ 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.
+ AliFMDCorrELossFit* fits = fcm.GetELossFit();
+ const TAxis& eAxis = fits->GetEtaAxis();
+
+ UShort_t nEta = eAxis.GetNbins();
+ fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
+ fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
+ fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
+ fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
+ fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
+ fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
+
+ UShort_t ybin = 0;
+ 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');
+ ybin++;
+ for (UShort_t e = 1; e <= nEta; e++) {
+ Double_t eta = eAxis.GetBinCenter(e);
+ Double_t cut = GetHighCut(d, r, eta, false);
+ if (cut <= 0) continue;
+ fHighCuts->SetBinContent(e, ybin, cut);
+ }
+ }
+ }
+}
+
//____________________________________________________________________
Bool_t
AliFMDSharingFilter::Filter(const AliESDFMD& input,
output.Clear();
TIter next(&fRingHistos);
RingHistos* o = 0;
- Double_t lowCut = GetLowCut();
while ((o = static_cast<RingHistos*>(next())))
o->Clear();
- // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ if (fOper) fOper->Reset(0);
+ Int_t nNone = 0;
+ Int_t nCandidate = 0;
+ Int_t nMerged = 0;
+ Int_t nSummed = 0;
+ Status status[512];
+
for(UShort_t d = 1; d <= 3; d++) {
Int_t nRings = (d == 1 ? 1 : 2);
for (UShort_t q = 0; q < nRings; q++) {
- Char_t r = (q == 0 ? 'I' : 'O');
- UShort_t nsec = (q == 0 ? 20 : 40);
- UShort_t nstr = (q == 0 ? 512 : 256);
-
+ Char_t r = (q == 0 ? 'I' : 'O');
+ UShort_t nsec = (q == 0 ? 20 : 40);
+ UShort_t nstr = (q == 0 ? 512 : 256);
RingHistos* histos = GetRingHistos(d, r);
for(UShort_t s = 0; s < nsec; s++) {
- Bool_t usedThis = kFALSE;
- Bool_t usedPrev = kFALSE;
-
+ for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate;
+#ifdef USE_OLDER_MERGING
+ Bool_t usedThis = kFALSE;
+ Bool_t usedPrev = kFALSE;
+#endif
for(UShort_t t = 0; t < nstr; t++) {
output.SetMultiplicity(d,r,s,t,0.);
Float_t mult = SignalInStrip(input,d,r,s,t);
+ // Get the pseudo-rapidity
+ Double_t eta = input.Eta(d,r,s,t);
+ Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
+ if (s == 0) output.SetEta(d,r,s,t,eta);
// Keep dead-channel information.
if(mult == AliESDFMD::kInvalidMult)
output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
// If no signal or dead strip, go on.
- if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
+ if (mult == AliESDFMD::kInvalidMult || mult == 0) {
+ if (mult == 0) histos->fSum->Fill(eta,phi,mult);
+ status[t] = kNone;
+ continue;
+ }
- // Get the pseudo-rapidity
- Double_t eta = input.Eta(d,r,s,t);
-
// Fill the diagnostics histogram
histos->fBefore->Fill(mult);
// Get next and previous signal - if any
Double_t prevE = 0;
Double_t nextE = 0;
+ Status prevStatus = (t == 0 ? kNone : status[t-1]);
+ Status thisStatus = status[t];
+ Status nextStatus = (t == nstr-1 ? kNone : status[t+1]);
if (t != 0) {
prevE = SignalInStrip(input,d,r,s,t-1);
if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
nextE = SignalInStrip(input,d,r,s,t+1);
if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
}
+ if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult);
+#ifdef USE_OLDER_MERGING
Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
lowFlux,d,r,s,t,
usedPrev,usedThis);
- if (mergedEnergy > lowCut) histos->Incr();
+ status[t] = (usedPrev ? kMergedWithOther : kNone);
+ if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone);
+#else
+ Double_t mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE,
+ eta, lowFlux,
+ d, r, s, t,
+ prevStatus,
+ thisStatus,
+ nextStatus);
+ if (t != 0) status[t-1] = prevStatus;
+ if (t != nstr-1) status[t+1] = nextStatus;
+ status[t] = thisStatus;
+
+#endif
+ // If we're processing on non-angle corrected data, we
+ // should do the angle correction here
+ if (!fCorrectAngles)
+ mergedEnergy = AngleCorrect(mergedEnergy, eta);
+ if (mergedEnergy > 0) histos->Incr();
+
+ if (t != 0)
+ histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1),
+ mergedEnergy);
+ histos->fBeforeAfter->Fill(mult, mergedEnergy);
histos->fAfter->Fill(mergedEnergy);
+ histos->fSum->Fill(eta,phi,mergedEnergy);
output.SetMultiplicity(d,r,s,t,mergedEnergy);
- output.SetEta(d,r,s,t,eta);
} // for strip
+ for (UShort_t t = 0; t < nstr; t++) {
+ if (fOper) fOper->operator()(d, r, s, t) = status[t];
+ switch (status[t]) {
+ case kNone: nNone++; break;
+ case kCandidate: nCandidate++; break;
+ case kMergedWithOther: nMerged++; break;
+ case kMergedInto: nSummed++; break;
+ }
+ }
} // for sector
} // for ring
} // for detector
+ fSummed->Fill(kNone, nNone);
+ fSummed->Fill(kCandidate, nCandidate);
+ fSummed->Fill(kMergedWithOther, nMerged);
+ fSummed->Fill(kMergedInto, nSummed);
+ DBGL(1, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d",
+ nNone, nCandidate, nMerged, nSummed));
next.Reset();
while ((o = static_cast<RingHistos*>(next())))
o->Finish();
// The energy signal
//
Double_t mult = input.Multiplicity(d,r,s,t);
- if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
-
- if (fCorrectAngles && !input.IsAngleCorrected())
- mult = AngleCorrect(mult, input.Eta(d,r,s,t));
- else if (!fCorrectAngles && input.IsAngleCorrected())
- mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
+ // In case of
+ // - bad value (invalid or 0)
+ // - we want angle corrected and data is
+ // - we don't want angle corrected and data isn't
+ // just return read value
+ if (mult == AliESDFMD::kInvalidMult ||
+ mult == 0 ||
+ (fCorrectAngles && input.IsAngleCorrected()) ||
+ (!fCorrectAngles && !input.IsAngleCorrected()))
+ return mult;
+
+ // If we want angle corrected data, correct it,
+ // otherwise de-correct it
+ if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
+ else mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
return mult;
}
//_____________________________________________________________________
Double_t
-AliFMDSharingFilter::GetLowCut() const
+AliFMDSharingFilter::GetLowCut(UShort_t, Char_t, Double_t) const
{
//
// Get the low cut. Normally, the low cut is taken to be the lower
//_____________________________________________________________________
Double_t
-AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
+AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r,
+ Double_t eta, Bool_t errors) const
{
//
// Get the high cut. The high cut is defined as the
// 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.
- 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));
+ AliFMDCorrELossFit* fits = fcm.GetELossFit();
+
+ return fits->GetLowerBound(d, r, eta, fNXi, errors, fIncludeSigma);
+}
+
+//_____________________________________________________________________
+namespace {
+ const char* status2String(AliFMDSharingFilter::Status s)
+ {
+ switch(s) {
+ case AliFMDSharingFilter::kNone: return "none ";
+ case AliFMDSharingFilter::kCandidate: return "candidate";
+ case AliFMDSharingFilter::kMergedWithOther: return "merged ";
+ case AliFMDSharingFilter::kMergedInto: return "summed ";
+ }
+ return "unknown ";
}
- else {
- delta = fit->fDelta;
- xi = fit->fXi;
+}
+
+//_____________________________________________________________________
+Double_t
+AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE,
+ Double_t prevE,
+ Double_t nextE,
+ Double_t eta,
+ Bool_t lowFlux,
+ UShort_t d,
+ Char_t r,
+ UShort_t s,
+ UShort_t t,
+ Status& prevStatus,
+ Status& thisStatus,
+ Status& nextStatus) const
+{
+ Double_t lowCut = GetLowCut(d, r, eta);
+
+ DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)",
+ d, r, s, t,
+ thisE, status2String(thisStatus),
+ prevE, status2String(prevStatus),
+ nextE, status2String(nextStatus)));
+
+ // If below cut, then modify zero signal and make sure the next
+ // strip is considered a candidate.
+ if (thisE < lowCut || thisE > 20) {
+ thisStatus = kNone;
+ DBGL(3,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE));
+ if (prevStatus == kCandidate) prevStatus = kNone;
+ return 0;
+ }
+ // It this strip was merged with the previous strip, then
+ // make the next strip a candidate and zero the value in this strip.
+ if (thisStatus == kMergedWithOther) {
+ DBGL(3,Form(" Merged with other, 0'ed"));
+ return 0;
}
- if (delta > 100) {
- AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
+ // Get the upper cut on the sharing
+ Double_t highCut = GetHighCut(d, r, eta ,false);
+ if (highCut <= 0) {
+ thisStatus = kNone;
+ return 0;
+ }
+
+ // Only for low-flux events
+ if (lowFlux) {
+ // If this is smaller than the next, and the next is larger
+ // then the cut, mark both as candidates for sharing.
+ // They should be be merged on the next iteration
+ if (thisE < nextE && nextE > highCut) {
+ thisStatus = kCandidate;
+ nextStatus = kCandidate;
+ return 0;
+ }
+ }
+
+ // Variable to sum signal in
+ Double_t totalE = thisE;
+
+ // If the previous signal was between the two cuts, and is still
+ // marked as a candidate , then merge it in here.
+ if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) {
+ totalE += prevE;
+ prevStatus = kMergedWithOther;
+ DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f",
+ lowCut, prevE, highCut, totalE));
+ }
+
+ // If the next signal is between the two cuts, then merge it here
+ if (nextE > lowCut && nextE < highCut) {
+ totalE += nextE;
+ nextStatus = kMergedWithOther;
+ DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE));
+ }
+
+ // If the signal is bigger than the high-cut, then return the value
+ if (totalE > highCut) {
+ if (totalE > thisE)
+ thisStatus = kMergedInto;
+ else
+ thisStatus = kNone;
+ DBGL(3, Form(" %9f>%f9 (was %9f) -> %9f %s",
+ totalE, highCut, thisE, totalE,status2String(thisStatus)));
+ return totalE;
+ }
+ else {
+ // This is below cut, so flag next as a candidate
+ DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut));
+ nextStatus = kCandidate;
+ }
+
+ // If the total signal is smaller than low cut then zero this and kill this
+ if (totalE < lowCut) {
+ DBGL(3, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE));
+ thisStatus = kNone;
+ return 0;
}
- Double_t highCut = (delta - fNXi * xi);
- return highCut;
+
+ // If total signal not above high cut or lower than low cut,
+ // mark this as a candidate for merging into the next, and zero signal
+ DBGL(3, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate",
+ lowCut, totalE, highCut, thisE));
+ thisStatus = kCandidate;
+ return 0;
}
//_____________________________________________________________________
// IF the multiplicity is very large, or below the cut, reject it, and
// flag it as candidate for sharing
- Double_t lowCut = GetLowCut();
+ Double_t lowCut = GetLowCut(d,r,eta);
if(mult > 12 || mult < lowCut) {
usedThis = kFALSE;
usedPrev = kFALSE;
//analyse and perform sharing on one strip
// Get the high cut
- Double_t highCut = GetHighCut(d, r, eta);
+ Double_t highCut = GetHighCut(d, r, eta, false);
+ if (highCut <= 0) {
+ usedThis = kFALSE;
+ usedPrev = kTRUE;
+ return 0;
+ }
// If this signal is smaller than the next, and the next signal is smaller
// than then the high cut, and it's a low-flux event, then mark this and
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
+ // the next signal is smaller than the high cut, then add the next signal
// to this, and marked the next signal as used
if(nextE > lowCut && nextE < highCut ) {
totalE += nextE;
// If we're processing on non-angle corrected data, we should do the
// angle correction here
- if (!fCorrectAngles)
- totalE = AngleCorrect(totalE, eta);
+ // if (!fCorrectAngles)
+ // totalE = AngleCorrect(totalE, eta);
// Fill post processing histogram
// if(totalE > fLowCut)
TIter next(&fRingHistos);
RingHistos* o = 0;
- while ((o = static_cast<RingHistos*>(next())))
+ THStack* sums = new THStack("sums", "Sum of ring signals");
+ while ((o = static_cast<RingHistos*>(next()))) {
o->ScaleHistograms(d, nEvents);
+ TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
+ sum->Scale(1., "width");
+ sum->SetTitle(o->GetName());
+ sum->SetDirectory(0);
+ sum->SetYTitle("#sum #Delta/#Delta_{mip}");
+ sums->Add(sum);
+ }
+ d->Add(sums);
}
//____________________________________________________________________
TList* d = new TList;
d->SetName(GetName());
dir->Add(d);
+
+ fSummed = new TH2I("operations", "Strip operations",
+ kMergedInto, kNone-.5, kMergedInto+.5,
+ 51201, -.5, 51200.5);
+ fSummed->SetXTitle("Operation");
+ fSummed->SetYTitle("# of strips");
+ fSummed->SetZTitle("Events");
+ fSummed->GetXaxis()->SetBinLabel(kNone, "None");
+ fSummed->GetXaxis()->SetBinLabel(kCandidate, "Candidate");
+ fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
+ fSummed->GetXaxis()->SetBinLabel(kMergedInto, "Merged into");
+ fSummed->SetDirectory(0);
+ d->Add(fSummed);
+
+ fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
+ fHighCuts->SetXTitle("#eta");
+ fHighCuts->SetDirectory(0);
+ d->Add(fHighCuts);
+
+
TIter next(&fRingHistos);
RingHistos* o = 0;
while ((o = static_cast<RingHistos*>(next()))) {
for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
ind[gROOT->GetDirLevel()] = '\0';
std::cout << ind << ClassName() << ": " << GetName() << '\n'
+ << std::boolalpha
+ << ind << " Debug: " << fDebug << "\n"
<< ind << " Low cut: " << fLowCut << '\n'
<< ind << " N xi factor: " << fNXi << '\n'
- << ind << " Use corrected angles: "
- << (fCorrectAngles ? "yes" : "no") << std::endl;
+ << ind << " Include sigma in cut: " << fIncludeSigma << '\n'
+ << ind << " Use corrected angles: " << fCorrectAngles
+ << std::noboolalpha << std::endl;
}
//====================================================================
: AliForwardUtil::RingHistos(),
fBefore(0),
fAfter(0),
+ fBeforeAfter(0),
+ fNeighborsBefore(0),
+ fNeighborsAfter(0),
+ fSum(0),
fHits(0),
fNHits(0)
{
: AliForwardUtil::RingHistos(d,r),
fBefore(0),
fAfter(0),
+ fBeforeAfter(0),
+ fNeighborsBefore(0),
+ fNeighborsAfter(0),
+ fSum(0),
fHits(0),
fNHits(0)
{
// d detector
// r ring
//
- fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
- Form("Energy loss in %s (reconstruction)", fName.Data()),
+ fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)",
1000, 0, 25);
- fAfter = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
- Form("Energy loss in %s (sharing corrected)",
- fName.Data()), 1000, 0, 25);
fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
- fBefore->SetFillColor(kRed+1);
+ fBefore->SetFillColor(Color());
fBefore->SetFillStyle(3001);
- // fBefore->Sumw2();
fBefore->SetDirectory(0);
- fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
- fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
- fAfter->SetFillColor(kBlue+1);
- fAfter->SetFillStyle(3001);
- // fAfter->Sumw2();
+
+ fAfter = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
+ fAfter->SetTitle("Energy loss in %s (sharing corrected)");
+ fAfter->SetFillColor(Color()+2);
fAfter->SetDirectory(0);
- fHits = new TH1D(Form("%s_Hits", fName.Data()),
- Form("Number of hits in %s", fName.Data()),
- 200, 0, 200000);
+ Double_t max = 15;
+ Double_t min = -1;
+ Int_t n = int((max-min) / (max / 300));
+ fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation",
+ n, min, max, n, min, max);
+ fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
+ fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
+ fBeforeAfter->SetZTitle("Correlation");
+ fBeforeAfter->SetDirectory(0);
+
+ fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
+ fNeighborsBefore->SetTitle("Correlation of neighbors befire");
+ fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
+ fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
+ fNeighborsBefore->SetDirectory(0);
+
+ fNeighborsAfter =
+ static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
+ fNeighborsAfter->SetTitle("Correlation of neighbors after");
+ fNeighborsAfter->SetDirectory(0);
+
+ fSum = new TH2D("summed", "Summed signal", 200, -4, 6,
+ (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi());
+ fSum->SetDirectory(0);
+ fSum->Sumw2();
+ fSum->SetMarkerColor(Color());
+ // fSum->SetFillColor(Color());
+ fSum->SetXTitle("#eta");
+ fSum->SetYTitle("#varphi [radians]");
+ fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) ");
+
+ fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
fHits->SetDirectory(0);
fHits->SetXTitle("# of hits");
fHits->SetFillColor(kGreen+1);
: AliForwardUtil::RingHistos(o),
fBefore(o.fBefore),
fAfter(o.fAfter),
+ fBeforeAfter(o.fBeforeAfter),
+ fNeighborsBefore(o.fNeighborsBefore),
+ fNeighborsAfter(o.fNeighborsAfter),
+ fSum(o.fSum),
fHits(o.fHits),
fNHits(o.fNHits)
{
if (fAfter) delete fAfter;
if (fHits) delete fHits;
- fBefore = static_cast<TH1D*>(o.fBefore->Clone());
- fAfter = static_cast<TH1D*>(o.fAfter->Clone());
- fHits = static_cast<TH1D*>(o.fHits->Clone());
-
+ fBefore = static_cast<TH1D*>(o.fBefore->Clone());
+ fAfter = static_cast<TH1D*>(o.fAfter->Clone());
+ fBeforeAfter = static_cast<TH2D*>(o.fBeforeAfter->Clone());
+ fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
+ fNeighborsAfter = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
+ fHits = static_cast<TH1D*>(o.fHits->Clone());
+ fSum = static_cast<TH2D*>(o.fSum->Clone());
+
return *this;
}
//____________________________________________________________________
//
// Destructor
//
- if (fBefore) delete fBefore;
- if (fAfter) delete fAfter;
- if (fHits) delete fHits;
}
//____________________________________________________________________
void
TList* l = GetOutputList(dir);
if (!l) return;
- TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
- fName.Data())));
- TH1D* after = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
- fName.Data())));
+ TH1D* before = static_cast<TH1D*>(l->FindObject("esdELoss"));
+ TH1D* after = static_cast<TH1D*>(l->FindObject("anaELoss"));
+ TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
if (before) before->Scale(1./nEvents);
if (after) after->Scale(1./nEvents);
+ if (summed) summed->Scale(1./nEvents);
+
}
//____________________________________________________________________
// dir where to store
//
TList* d = DefineOutputList(dir);
+
d->Add(fBefore);
d->Add(fAfter);
+ d->Add(fBeforeAfter);
+ d->Add(fNeighborsBefore);
+ d->Add(fNeighborsAfter);
d->Add(fHits);
+ d->Add(fSum);
+
dir->Add(d);
}
class TAxis;
class TList;
class TH2;
+class AliFMDFloatMap;
/**
* Class to do the sharing correction. That is, a filter that merges
class AliFMDSharingFilter : public TNamed
{
public:
+ enum Status {
+ kNone = 1,
+ kCandidate = 2,
+ kMergedWithOther = 3,
+ kMergedInto = 4
+ };
/**
* Destructor
*/
*/
AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o);
+ /**
+ * Initialize
+ *
+ */
+ void Init();
/**
* Set the low cut used for sharing
*
* otherwise use de-corrected signals. In the final output, the
* signals are always angle corrected.
*/
- void UseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
+ void SetUseAngleCorrectedSignals(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; }
+ void SetNXi(Double_t n) { fNXi = n; }
+ /**
+ * Whether to include sigma in the number subtracted from the MPV to
+ * get the high cut
+ *
+ * @param u If true, then high cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$
+ */
+ void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; }
/**
* Filter the input AliESDFMD object
*
void ScaleHistograms(const TList* dir, Int_t nEvents);
TH1D* fBefore; // Distribution of signals before filter
TH1D* fAfter; // Distribution of signals after filter
+ TH2D* fBeforeAfter; // Correlation of before and after
+ TH2D* fNeighborsBefore; // Correlation of neighbors
+ TH2D* fNeighborsAfter; // Correlation of neighbors
+ TH2D* fSum; // Summed signal
TH1D* fHits; // Distribution of hit strips.
Int_t fNHits; // Number of hit strips per event
ClassDef(RingHistos,1);
UShort_t t,
Bool_t& usedPrev,
Bool_t& usedThis) const;
+ Double_t MultiplicityOfStrip(Double_t thisE,
+ Double_t prevE,
+ Double_t nextE,
+ Double_t eta,
+ Bool_t lowFlux,
+ UShort_t d,
+ Char_t r,
+ UShort_t s,
+ UShort_t t,
+ Status& prevStatus,
+ Status& thisStatus,
+ Status& nextStatus) const;
/**
* Angle correct the signal
*
* @return Angle un-corrected signal
*/
Double_t DeAngleCorrect(Double_t mult, Double_t eta) const;
- /**
+ /**
* 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.
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta value
+ * @param errors If false, do not show errors
+ *
+ * @return 0 or less on failure, otherwise @f$\Delta_{mp}-n\xi@f$
*/
- virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta) const;
+ virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta,
+ Bool_t errors=true) 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.
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta value
+ *
+ * @return
*/
- virtual Double_t GetLowCut() const;
+ virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) 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
+ Double_t fNXi; // Number of xi's from Delta to stop merging
+ Bool_t fIncludeSigma; // Whether to include sigma in cut
+ TH2* fSummed; // Operations histogram
+ TH2* fHighCuts; // High cuts used
+ AliFMDFloatMap* fOper; // Operation done per strip
Int_t fDebug; // Debug level
- ClassDef(AliFMDSharingFilter,1); //
+ ClassDef(AliFMDSharingFilter,2); //
};
#endif
#include "AliInputEventHandler.h"
#include "AliStack.h"
#include "AliMCEvent.h"
-//#include "AliFMDStripIndex.h"
+#include "AliAODForwardMult.h"
+#include "AliFMDStripIndex.h"
#include <TH1.h>
#include <TH3D.h>
#include <TH2D.h>
//====================================================================
AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
: AliAnalysisTaskSE(),
+ fInspector(),
+ fFirstEvent(true),
fHEvents(0),
fHEventsTr(0),
fHEventsTrVtx(0),
- fHEventsVtx(0),
fHTriggers(0),
fPrimaryInnerAll(0),
fPrimaryOuterAll(0),
//____________________________________________________________________
AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
- : AliAnalysisTaskSE(name),
+ : AliAnalysisTaskSE(name),
+ fInspector("eventInspector"),
+ fFirstEvent(true),
fHEvents(0),
fHEventsTr(0),
fHEventsTrVtx(0),
- fHEventsVtx(0),
fHTriggers(0),
fPrimaryInnerAll(0),
fPrimaryOuterAll(0),
// name Name of task
//
DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
}
//____________________________________________________________________
AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
: AliAnalysisTaskSE(o),
+ fInspector(o.fInspector),
+ fFirstEvent(o.fFirstEvent),
fHEvents(o.fHEvents),
fHEventsTr(o.fHEventsTr),
fHEventsTrVtx(o.fHEventsTrVtx),
- fHEventsVtx(o.fHEventsVtx),
fHTriggers(o.fHTriggers),
fPrimaryInnerAll(o.fPrimaryInnerAll),
fPrimaryOuterAll(o.fPrimaryOuterAll),
// Return:
// Reference to this object
//
+ fInspector = o.fInspector;
+ fFirstEvent = o.fFirstEvent;
fHEventsTr = o.fHEventsTr;
fHEventsTrVtx = o.fHEventsTrVtx;
fHTriggers = o.fHTriggers;
//
//
fList = new TList;
- fList->SetName(GetName());
+ fList->SetOwner();
+ fList->SetName(Form("%sSums", GetName()));
fHEvents = new TH1I(GetEventName(false,false),
"Number of all events",
fHEventsTrVtx->SetDirectory(0);
fList->Add(fHEventsTrVtx);
- fHEventsVtx = new TH1I(GetEventName(false,true),
- "Number of events w/vertex",
- fVtxAxis.GetNbins(),
- fVtxAxis.GetXmin(),
- fVtxAxis.GetXmax());
- fHEventsVtx->SetXTitle("v_{z} [cm]");
- fHEventsVtx->SetYTitle("# of events");
- fHEventsVtx->SetFillColor(kBlue+1);
- fHEventsVtx->SetFillStyle(3001);
- fHEventsVtx->SetDirectory(0);
- fList->Add(fHEventsVtx);
-
-
fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
fHTriggers->SetFillColor(kRed+1);
fHTriggers->SetFillStyle(3001);
strips->Add(fStripsFMD3o);
fList->Add(strips);
+ fInspector.Init(fVtxAxis);
+ fInspector.DefineOutput(fList);
+
PostData(1, fList);
}
//____________________________________________________________________
AliWarning("No ESD event found for input event");
return;
}
-
- // Get the particle stack
- AliStack* stack = mcEvent->Stack();
- // Get the event generate header
- AliHeader* mcHeader = mcEvent->Header();
- AliGenEventHeader* genHeader = mcHeader->GenEventHeader();
-
- // Get the generator vertex
- TArrayF mcVertex;
- genHeader->PrimaryVertex(mcVertex);
- Double_t mcVtxZ = mcVertex.At(2);
-
- // Check the MC vertex is in range
- 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]",
- mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
-#endif
- return;
+ //--- Read run information -----------------------------------------
+ if (fFirstEvent && esd->GetESDRun()) {
+ fInspector.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()));
+
+ fFirstEvent = false;
}
- // UInt_t triggers;
- // Bool_t gotTrigggers = false;
- Bool_t gotInel = false;
- // Double_t vZ;
- Bool_t gotVertex = false;
-#if 0
- // Use event inspector instead
- // Get the triggers
- UInt_t triggers = 0;
- Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers);
- Bool_t gotInel = triggers & AliAODForwardMult::kInel;
-
- // Get the ESD vertex
- Double_t vZ = -1000000;
- Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ);
-#endif
-
+ // Get the particle stack
+ AliStack* stack = mcEvent->Stack();
+ // Some variables
+ UInt_t triggers; // Trigger bits
+ Bool_t lowFlux; // Low flux flag
+ UShort_t iVz; // Vertex bin from ESD
+ Double_t vZ; // Z coordinate from ESD
+ Double_t cent; // Centrality
+ UShort_t iVzMc; // Vertex bin from MC
+ Double_t vZMc; // Z coordinate of IP vertex from MC
+ Double_t b; // Impact parameter
+ Int_t nPart; // Number of participants
+ Int_t nBin; // Number of binary collisions
+ Double_t phiR; // Reaction plane from MC
+ UShort_t nClusters;// Number of SPD clusters
+ // Process the data
+ UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
+ cent, nClusters);
+ /*UInt_t retMC =*/ fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
+ b, nPart, nBin, phiR);
+
+ // Bool_t isInelMC = true; // (triggers & AliAODForwardMult::kB);
+ // Bool_t isNSDMC = (triggers & AliAODForwardMult::kMCNSD);
+ Bool_t isInel = triggers & AliAODForwardMult::kInel;
+ // Bool_t isNSD = triggers & AliAODForwardMult::kNSD;
+ Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk;
// Fill the event count histograms
- if (gotInel) fHEventsTr->Fill(mcVtxZ);
- if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ);
- if (gotVertex) fHEventsVtx->Fill(mcVtxZ);
- fHEvents->Fill(mcVtxZ);
+ if (isInel) fHEventsTr->Fill(vZMc);
+ if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
+ fHEvents->Fill(vZMc);
// Cache of the hits
AliFMDFloatMap hitsByStrip;
Double_t phi = particle->Phi();
if (isCharged && isPrimary)
- FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi);
+ FillPrimary(isInel, hasVtx, vZMc, eta, phi);
// For the rest, ignore non-collisions, and events out of vtx range
- if (!gotInel || gotVertex) continue;
+ if (!isInel || !hasVtx) continue;
Int_t nTrRef = particle->GetNumberOfTrackReferences();
for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
// Get the detector coordinates
UShort_t d = 0, s = 0, t = 0;
Char_t r = '\0';
- // AliFMDStripIndex::Unpack(trackRef->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
// Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
// Fill strip information into histograms
- FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
+ FillStrip(d, r, vZMc, eta, phi, hitsByStrip(d,r,s,t) == 1);
// Set the last processed track number - marking it as done for
// this strip
//____________________________________________________________________
void
-AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
- Double_t vz, Double_t eta, Double_t phi)
+AliForwardMCCorrectionsTask::FillPrimary(Bool_t isInel, Bool_t hasVtx,
+ Double_t vz,
+ Double_t eta, Double_t phi)
{
//
// Fill in primary information
//
// Parameters:
- // gotInel Got INEL trigger from ESD
+ // isInel Got INEL trigger from ESD
// gotVtx Got vertex Z from ESD
// vz @f$z@f$ coordinate of interation point
// eta Pseudo rapidity
// phi Azimuthal angle
//
- if (gotInel && gotVtx) {
+ if (isInel && hasVtx) {
// This takes the place of hPrimary_FMD_<r>_vtx<v> and
// Analysed_FMD<r>_vtx<v>
fPrimaryInnerTrVtx->Fill(vz,eta,phi);
// option Not used
//
- TList* list = dynamic_cast<TList*>(GetOutputData(1));
- if (!list) {
+ fList = dynamic_cast<TList*>(GetOutputData(1));
+ if (!fList) {
AliError("No output list defined");
return;
}
- TList* primaries = static_cast<TList*>(list->FindObject("primaries"));
- TList* hits = static_cast<TList*>(list->FindObject("hits"));
- TList* strips = static_cast<TList*>(list->FindObject("strips"));
+ TList* primaries = static_cast<TList*>(fList->FindObject("primaries"));
+ TList* hits = static_cast<TList*>(fList->FindObject("hits"));
+ TList* strips = static_cast<TList*>(fList->FindObject("strips"));
if (!primaries || !hits || !strips) {
AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)",
- list->GetName(), primaries, hits, strips));
+ fList->GetName(), primaries, hits, strips));
return;
}
TH1I* eventsAll =
- static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
+ static_cast<TH1I*>(fList->FindObject(GetEventName(false,false)));
TH1I* eventsTr =
- static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
+ static_cast<TH1I*>(fList->FindObject(GetEventName(true,false)));
TH1I* eventsVtx =
- static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
+ static_cast<TH1I*>(fList->FindObject(GetEventName(false,true)));
TH1I* eventsTrVtx =
- static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
+ static_cast<TH1I*>(fList->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(),
+ "(a=%p,t=%p,v=%p,tv=%p)", fList->GetName(),
eventsAll, eventsTr, eventsVtx, eventsTrVtx));
return;
}
static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',true)));
if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) {
AliError(Form("Could not find all primary particle histograms in %s "
- "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(),
+ "(ai=%p,ao=%p,tvi=%p,tvo=%p)", fList->GetName(),
primIAll, primOAll, primITrVtx, primOTrVtx));
return;
}
return;
}
+ // Output list
+ TList* output = new TList;
+ output->SetOwner();
+ output->SetName(Form("%sResults", GetName()));
+
// Calculate the over-all event selection efficiency
TH1D* selEff = new TH1D("selEff", "Event selection efficiency",
fVtxAxis.GetNbins(),
selEff->SetFillStyle(3001);
selEff->Add(eventsAll);
selEff->Divide(eventsTrVtx);
- list->Add(selEff);
+ output->Add(selEff);
// Loop over vertex bins and do vertex dependendt stuff
for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) {
- // Make a sub-list in the output
+ // Make a sub-fList in the output
TList* vl = new TList;
+ vl->SetOwner();
vl->SetName(Form("vtx%02d", v));
- list->Add(vl);
+ output->Add(vl);
// Get event counts
Int_t nEventsAll = Int_t(eventsAll->GetBinContent(v));
doubleHit->SetFillStyle(3001);
doubleHit->Sumw2();
doubleHit->Divide(hStrips);
- list->Add(doubleHit);
- }
- }
- }
+ output->Add(doubleHit);
+ } // for q
+ } // for d
+ } // for v
+ PostData(2, output);
}
//____________________________________________________________________
*/
#include <AliAnalysisTaskSE.h>
#include <AliESDFMD.h>
+#include "AliFMDMCEventInspector.h"
#include <TH1I.h>
-class AliFMDAnaParameters;
class AliESDEvent;
class TH2D;
class TH3D;
void FillStrip(UShort_t d, Char_t r,
Double_t vz, Double_t eta, Double_t phi,
Bool_t first);
+ AliFMDMCEventInspector fInspector; // Event inspector
+ Bool_t fFirstEvent; // First event flag
TH1I* fHEvents; // All Events
TH1I* fHEventsTr; // Histogram of events w/trigger
TH1I* fHEventsTrVtx; // Events w/trigger and vertex
- TH1I* fHEventsVtx; // Events w/vertex
TH1I* fHTriggers; // Triggers
TH3D* fPrimaryInnerAll; // Distribution of primaries - all events
TH3D* fPrimaryOuterAll; // Distribution of primaries - all events
fMCESDFMD(),
fMCHistos(),
fMCAODFMD(),
+ fRingSums(),
+ fMCRingSums(),
fPrimary(0),
fEventInspector(),
fSharingFilter(),
fMCESDFMD(),
fMCHistos(),
fMCAODFMD(kTRUE),
+ fRingSums(),
+ fMCRingSums(),
fPrimary(0),
fEventInspector("event"),
fSharingFilter("sharing"),
fMCESDFMD(o.fMCESDFMD),
fMCHistos(o.fMCHistos),
fMCAODFMD(o.fMCAODFMD),
+ fRingSums(o.fRingSums),
+ fMCRingSums(o.fMCRingSums),
fPrimary(o.fPrimary),
fEventInspector(o.fEventInspector),
fSharingFilter(o.fSharingFilter),
fAODFMD = o.fAODFMD;
fMCHistos = o.fMCHistos;
fMCAODFMD = o.fMCAODFMD;
+ fRingSums = o.fRingSums;
+ fMCRingSums = o.fMCRingSums;
fPrimary = o.fPrimary;
fList = o.fList;
fCorrections.SetDebug(dbg);
fHistCollector.SetDebug(dbg);
}
+//____________________________________________________________________
+void
+AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
+{
+ fSharingFilter.SetOnlyPrimary(use);
+ fCorrections.SetSecondaryForMC(!use);
+}
+
//____________________________________________________________________
void
AliForwardMCMultiplicityTask::InitializeSubs()
fAODFMD.Init(*pe);
fMCHistos.Init(*pe);
fMCAODFMD.Init(*pe);
+ fRingSums.Init(*pe);
+ fMCRingSums.Init(*pe);
fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
fHData->SetStats(0);
fList->Add(fHData);
+ TList* rings = new TList;
+ rings->SetName("ringSums");
+ rings->SetOwner();
+ fList->Add(rings);
+
+ rings->Add(fRingSums.Get(1, 'I'));
+ rings->Add(fRingSums.Get(2, 'I'));
+ rings->Add(fRingSums.Get(2, 'O'));
+ rings->Add(fRingSums.Get(3, 'I'));
+ rings->Add(fRingSums.Get(3, 'O'));
+ fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+ fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+ fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+ fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+ fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
+ TList* mcRings = new TList;
+ mcRings->SetName("mcRingSums");
+ mcRings->SetOwner();
+ fList->Add(mcRings);
+
+ mcRings->Add(fMCRingSums.Get(1, 'I'));
+ mcRings->Add(fMCRingSums.Get(2, 'I'));
+ mcRings->Add(fMCRingSums.Get(2, 'O'));
+ mcRings->Add(fMCRingSums.Get(3, 'I'));
+ mcRings->Add(fMCRingSums.Get(3, 'O'));
+ fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+ fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+ fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+ fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+ fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
+
+
fEventInspector.Init(*pv);
+ fSharingFilter.Init();
fDensityCalculator.Init(*pe);
fCorrections.Init(*pe);
fHistCollector.Init(*pv,*pe);
fMCAODFMD.Clear();
fPrimary->Reset();
- Bool_t lowFlux = kFALSE;
- UInt_t triggers = 0;
- UShort_t ivz = 0;
- Double_t vz = 0;
- Double_t cent = 0;
- UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent);
+ Bool_t lowFlux = kFALSE;
+ UInt_t triggers = 0;
+ UShort_t ivz = 0;
+ Double_t vz = 0;
+ Double_t cent = 0;
+ UShort_t nClusters = 0;
+ UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
+ ivz, vz, cent, nClusters);
UShort_t ivzMC = 0;
Double_t vzMC = 0;
Double_t phiR = 0;
fAODFMD.SetSNN(fEventInspector.GetEnergy());
fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
fAODFMD.SetCentrality(cent);
+ fAODFMD.SetNClusters(nClusters);
fMCAODFMD.SetTriggerBits(triggers);
fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
fMCAODFMD.SetCentrality(cent);
+ fMCAODFMD.SetNClusters(nClusters);
//All events should be stored - HHD
//if (isAccepted) MarkEventForStore();
}
fCorrections.CompareResults(fHistos, fMCHistos);
- if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
+ if (!fHistCollector.Collect(fHistos, fRingSums,
+ ivz, fAODFMD.GetHistogram())) {
AliWarning("Histogram collector failed");
return;
}
- if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
+ if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
+ ivz, fMCAODFMD.GetHistogram())) {
AliWarning("MC Histogram collector failed");
return;
}
list->Add(dNdeta);
list->Add(norm);
+ MakeRingdNdeta(list, "ringSums", list, "ringResults");
+ MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
+
fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
/**
* @}
*/
+ /**
+ * Process only primary MC tracks
+ *
+ * @param use
+ */
+ void SetOnlyPrimary(Bool_t use);
/**
* @{
* @name Access to sub-algorithms
AliESDFMD fMCESDFMD; // MC 'Sharing corrected' ESD object
AliForwardUtil::Histos fMCHistos; // MC Cache histograms
AliAODForwardMult fMCAODFMD; // MC Output object
+ AliForwardUtil::Histos fRingSums; // Cache histograms
+ AliForwardUtil::Histos fMCRingSums; // Cache histograms
TH2D* fPrimary; // Per event primary particles
AliFMDMCEventInspector fEventInspector; // Algorithm
#include "AliESDEvent.h"
#include <TROOT.h>
#include <TAxis.h>
+#include <THStack.h>
#include <iostream>
#include <iomanip>
ah->SetFillAOD(kTRUE);
}
+//____________________________________________________________________
+void
+AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
+ const char* inName,
+ TList* output,
+ const char* outName,
+ Int_t style) const
+{
+ // Make dN/deta for each ring found in the input list.
+ //
+ // A stack of all the dN/deta is also made for easy drawing.
+ //
+ // Note, that the distributions are normalised to the number of
+ // observed events only - they should be corrected for
+ if (!input) return;
+ TList* list = static_cast<TList*>(input->FindObject(inName));
+ if (!list) {
+ AliWarning(Form("No list %s found in %s", inName, input->GetName()));
+ return;
+ }
+
+ TList* out = new TList;
+ out->SetName(outName);
+ out->SetOwner();
+ output->Add(out);
+
+ THStack* dndetaRings = new THStack("all", "dN/d#eta per ring");
+ const char* names[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
+ const char** ptr = names;
+
+ while (*ptr) {
+ TList* thisList = new TList;
+ thisList->SetOwner();
+ thisList->SetName(*ptr);
+ out->Add(thisList);
+
+ TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
+ if (!h) {
+ AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
+ ptr++;
+ continue;
+ }
+ TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
+ copy->SetDirectory(0);
+ thisList->Add(copy);
+
+ TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
+ for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
+ for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
+ Double_t c = copy->GetBinContent(i, j);
+ Double_t e = copy->GetBinError(i, j);
+ Double_t a = norm->GetBinContent(i);
+ copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
+ copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
+ }
+ }
+
+ TH1D* res =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
+ h->GetNbinsY(),"e"));
+ TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
+ res->SetTitle(*ptr);
+ res->Scale(1., "width");
+ copy->Scale(1., "width");
+ proj->Scale(1. / norm->GetMaximum(), "width");
+ norm->Scale(1. / norm->GetMaximum());
+
+ res->SetMarkerStyle(style);
+ norm->SetDirectory(0);
+ res->SetDirectory(0);
+ proj->SetDirectory(0);
+ thisList->Add(norm);
+ thisList->Add(res);
+ thisList->Add(proj);
+ dndetaRings->Add(res);
+ ptr++;
+ }
+ out->Add(dndetaRings);
+}
+
//____________________________________________________________________
void
AliForwardMultiplicityBase::Print(Option_t* option) const
*
*/
virtual void MarkEventForStore() const;
-
+ /**
+ * Make Ring @f$ dN/d\eta @f$ histogram and a stack
+ *
+ * @param input List with summed signals
+ * @param output Output list
+ * @param stackName Stack name
+ */
+ virtual void MakeRingdNdeta(const TList* input,
+ const char* inName,
+ TList* output,
+ const char* outName,
+ Int_t style=20) const;
Bool_t fEnableLowFlux;// Whether to use low-flux specific code
Bool_t fFirstEvent; // Whether the event is the first seen
private:
fESDFMD(),
fHistos(),
fAODFMD(),
+ fRingSums(),
fEventInspector(),
fSharingFilter(),
fDensityCalculator(),
fESDFMD(),
fHistos(),
fAODFMD(false),
+ fRingSums(),
fEventInspector("event"),
fSharingFilter("sharing"),
fDensityCalculator("density"),
fESDFMD(o.fESDFMD),
fHistos(o.fHistos),
fAODFMD(o.fAODFMD),
+ fRingSums(o.fRingSums),
fEventInspector(o.fEventInspector),
fSharingFilter(o.fSharingFilter),
fDensityCalculator(o.fDensityCalculator),
fHistCollector = o.fHistCollector;
fHistos = o.fHistos;
fAODFMD = o.fAODFMD;
+ fRingSums = o.fRingSums;
fList = o.fList;
return *this;
fHistos.Init(*pe);
fAODFMD.Init(*pe);
+ fRingSums.Init(*pe);
fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
fHData->SetStats(0);
fHData->SetDirectory(0);
fList->Add(fHData);
+ TList* rings = new TList;
+ rings->SetName("ringSums");
+ rings->SetOwner();
+ fList->Add(rings);
+
+ rings->Add(fRingSums.Get(1, 'I'));
+ rings->Add(fRingSums.Get(2, 'I'));
+ rings->Add(fRingSums.Get(2, 'O'));
+ rings->Add(fRingSums.Get(3, 'I'));
+ rings->Add(fRingSums.Get(3, 'O'));
+ fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+ fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+ fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+ fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+ fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
fEventInspector.Init(*pv);
+ fSharingFilter.Init();
fDensityCalculator.Init(*pe);
fCorrections.Init(*pe);
fHistCollector.Init(*pv,*pe);
fESDFMD.Clear();
fAODFMD.Clear();
- Bool_t lowFlux = kFALSE;
- UInt_t triggers = 0;
- UShort_t ivz = 0;
- Double_t vz = 0;
- Double_t cent = -1;
- UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent);
+ Bool_t lowFlux = kFALSE;
+ UInt_t triggers = 0;
+ UShort_t ivz = 0;
+ Double_t vz = 0;
+ Double_t cent = -1;
+ UShort_t nClusters = 0;
+ UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
+ ivz, vz, cent, nClusters);
if (found & AliFMDEventInspector::kNoEvent) return;
if (found & AliFMDEventInspector::kNoTriggers) return;
fAODFMD.SetSNN(fEventInspector.GetEnergy());
fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
fAODFMD.SetCentrality(cent);
+ fAODFMD.SetNClusters(nClusters);
MarkEventForStore();
if (found & AliFMDEventInspector::kNoSPD) return;
return;
}
- if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
+ if (!fHistCollector.Collect(fHistos, fRingSums,
+ ivz, fAODFMD.GetHistogram())) {
AliWarning("Histogram collector failed");
return;
}
list->Add(dNdeta);
list->Add(norm);
+ MakeRingdNdeta(list, "ringSums", list, "ringResults");
+
fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
AliESDFMD fESDFMD; // Sharing corrected ESD object
AliForwardUtil::Histos fHistos; // Cache histograms
AliAODForwardMult fAODFMD; // Output object
+ AliForwardUtil::Histos fRingSums; // Cache histograms
AliFMDEventInspector fEventInspector; // Algorithm
AliFMDSharingFilter fSharingFilter; // Algorithm
//
if (!d) return 0;
TList* list = new TList;
+ list->SetOwner();
list->SetName(fName.Data());
d->Add(list);
return list;
static Color_t RingColor(UShort_t d, Char_t r)
{
return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
- + ((r == 'I' || r == 'i') ? 2 : -2));
+ + ((r == 'I' || r == 'i') ? 2 : -3));
}
//==================================================================
/**
{
return AliForwardUtil::RingColor(fDet, fRing);
}
+ const char* GetName() const { return fName.Data(); }
UShort_t fDet; // Detector
Char_t fRing; // Ring
TString fName; // Name
#include <TH1D.h>
#include <THStack.h>
#include <TList.h>
+#include <TFile.h>
#include <AliAnalysisManager.h>
#include <AliAODEvent.h>
#include <AliAODHandler.h>
bin->ProcessPrimary(forward, fTriggerMask, fVtxMin, fVtxMax, primary);
}
+//________________________________________________________________________
+void
+AliForwarddNdetaTask::Terminate(Option_t *option)
+{
+ //
+ // Called at end of event processing..
+ //
+ // This is called once in the master
+ //
+ // Parameters:
+ // option Not used
+ AliBasedNdetaTask::Terminate(option);
+
+ THStack* truth = new THStack("dndetaTruth", "dN/d#eta MC Truth");
+ THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin),
+ "dN/d#eta MC Truth");
+
+ TIter next(fListOfCentralities);
+ CentralityBin* bin = 0;
+ while ((bin = static_cast<CentralityBin*>(next()))) {
+ if (fCentAxis && bin->IsAllBin()) continue;
+
+ TList* results = bin->GetResults();
+ if (!results) continue;
+
+ TH1* dndeta = static_cast<TH1*>(results->FindObject("dndetaTruth"));
+ TH1* dndetaRebin =
+ static_cast<TH1*>(results->FindObject(Form("dndetaTruth_rebin%02d",
+ fRebin)));
+ if (dndeta) truth->Add(dndeta);
+ if (dndetaRebin) truthRebin->Add(dndetaRebin);
+ }
+ // If available output rebinned stack
+ if (!truth->GetHists() ||
+ truth->GetHists()->GetEntries() <= 0) {
+ AliWarning("No MC truth histograms found");
+ delete truth;
+ truth = 0;
+ }
+ if (truth) fOutput->Add(truth);
+
+ // If available output rebinned stack
+ if (!truthRebin->GetHists() ||
+ truthRebin->GetHists()->GetEntries() <= 0) {
+ AliWarning("No rebinned MC truth histograms found");
+ delete truthRebin;
+ truthRebin = 0;
+ }
+ if (truthRebin) fOutput->Add(truthRebin);
+
+}
//____________________________________________________________________
TH2D*
AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult*
forward,
Int_t triggerMask,
- Double_t /*vzMin*/,
- Double_t /*vzMax*/,
+ Double_t vzMin,
+ Double_t vzMax,
const TH2D* primary)
{
// Check the centrality class unless this is the 'all' bin
// translate real trigger mask to MC trigger mask
Int_t mask = AliAODForwardMult::kB;
- if (triggerMask == AliAODForwardMult::kNSD)
- mask = AliAODForwardMult::kMCNSD;
+ if (triggerMask == AliAODForwardMult::kNSD) {
+ mask ^= AliAODForwardMult::kNSD;
+ mask = AliAODForwardMult::kMCNSD;
+ }
// Now use our normal check, but with the new mask, except
- if (!forward->CheckEvent(mask, -1000, -1000, 0, 0, 0)) return;
+ vzMin = vzMax = -10000; // ignore vertex
+ if (!forward->CheckEvent(mask, vzMin, vzMax, 0, 0, 0)) return;
fSumPrimary->Add(primary);
+ Int_t n = fSumPrimary->GetBinContent(0,0);
+ fSumPrimary->SetBinContent(0,0, ++n);
}
//________________________________________________________________________
Bool_t corrEmpty,
Bool_t cutEdges,
Int_t triggerMask,
- Int_t color,
Int_t marker)
{
AliInfo(Form("In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d",
shapeCorr, trigEff,
symmetrice, rebin,
rootProj, corrEmpty, cutEdges,
- triggerMask, color, marker);
+ triggerMask, marker);
fSumPrimary = static_cast<TH2D*>(fSums->FindObject("truth"));
if (!fSumPrimary) {
AliWarning("No MC truth histogram found");
- fSums->ls();
- return;
}
- Int_t n = (triggerMask == AliAODForwardMult::kNSD ?
- Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) :
- Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
- AliInfo(Form("Normalising MC truth to %d", n));
-
- TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
+ else {
+#if 0
+ Int_t n = fSumPrimary->GetBinContent(0,0);
+#else
+ Int_t n = (triggerMask == AliAODForwardMult::kNSD ?
+ Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) :
+ Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
+#endif
+ AliInfo(Form("Normalising MC truth to %d", n));
+
+ TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
fSumPrimary->GetNbinsY(),"e");
- dndetaTruth->Scale(1./n, "width");
+ dndetaTruth->SetDirectory(0);
+ dndetaTruth->Scale(1./n, "width");
- SetHistogramAttributes(dndetaTruth, color-2, marker+4, "Monte-Carlo truth");
- fOutput->Add(dndetaTruth);
- fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
+ SetHistogramAttributes(dndetaTruth, GetColor(), 30, "Monte-Carlo truth");
+
+ fOutput->Add(dndetaTruth);
+ fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
+ }
+
+ if (!IsAllBin()) return;
+ TFile* file = TFile::Open("forward.root", "READ");
+ if (!file) return;
+
+ TList* forward = static_cast<TList*>(file->Get("Forward"));
+ if (!forward) {
+ AliError("List Forward not found in forward.root");
+ return;
+ }
+ TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
+ if (!rings) {
+ AliError("List ringResults not found in forward.root");
+ return;
+ }
+ THStack* res = static_cast<THStack*>(rings->FindObject("all"));
+ if (!res) {
+ AliError(Form("Stack all not found in %s", rings->GetName()));
+ return;
+ }
+ if (!fTriggers) {
+ AliError("Triggers histogram not set");
+ return;
+ }
+ Double_t ntotal = 0;
+ Double_t epsilonT = trigEff;
+ // TEMPORARY FIX
+ if (triggerMask == AliAODForwardMult::kNSD) {
+ // This is a local change
+ epsilonT = 0.92;
+ AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
+ }
+ AliInfo("Adding per-ring histograms to output");
+ Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
+ TIter next(res->GetHists());
+ TH1* hist = 0;
+ while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
+ res->SetName("dndetaRings");
+ fOutput->Add(res);
}
//________________________________________________________________________
* @param option Not used
*/
virtual void UserExec(Option_t* option);
+ /**
+ * Called at end of event processing.
+ *
+ * This is called once in the master
+ *
+ * @param option Not used
+ */
+ virtual void Terminate(Option_t* option);
protected:
/**
* Copy constructor
Bool_t corrEmpty,
Bool_t cutEdges,
Int_t triggerMask,
- Int_t color,
Int_t marker);
protected:
TH2D* fSumPrimary; // Sum of primary histograms
* @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
* @date Wed Mar 23 14:07:10 2011
*
- * @brief Script to visualise the dN/deta
+ * @brief Script to visualise the dN/deta for pp and PbPb
*
* This script is independent of any AliROOT code - and should stay
* that way.
*
+ * The script is <emph>very</emph> long - sigh - the joy of drawing
+ * things nicely in ROOT
*
* @ingroup pwg2_forward_dndeta
*/
sys,snn,trg,
centLow,centHigh,onlya));
if (!ret) {
+#if 0
Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
"trigger=%d %d%%-%d%% central %s",
sys, snn, trg, centLow, centHigh,
onlya ? " (ALICE results)" : "all");
+#endif
return 0;
}
thisOther = reinterpret_cast<TMultiGraph*>(ret);
Error("FetchResults", "Couldn't find list 'all' in %s",
list->GetName());
else
- FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
+ FetchResults(all, name, FetchOthers(0,0), -1000, 0, max, rmax, amax);
return;
}
- Int_t nCol = gStyle->GetNumberOfColors();
for (UShort_t i = 0; i < n; i++) {
- UShort_t centLow = fCentAxis->GetXbins()->At(i);
- UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
+ UShort_t centLow = fCentAxis->GetBinLowEdge(i+1);
+ UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
TString lname = Form("cent%03d_%03d", centLow, centHigh);
TList* thisCent = static_cast<TList*>(list->FindObject(lname));
-
- Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
- Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
- Int_t col = gStyle->GetColorPalette(icol);
- Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d",
- centLow, centHigh, icol, nCol, fc, col);
+ Int_t col = GetCentralityColor(i+1);
TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
if (!thisCent)
}
}
//__________________________________________________________________
+ Int_t GetCentralityColor(Int_t bin) const
+ {
+ UShort_t centLow = fCentAxis->GetBinLowEdge(bin);
+ UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
+ Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
+ Int_t nCol = gStyle->GetNumberOfColors();
+ Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
+ Int_t col = gStyle->GetColorPalette(icol);
+ //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
+ return col;
+ }
+ //__________________________________________________________________
void SetAttributes(TH1* h, Int_t color)
{
if (!h) return;
//__________________________________________________________________
void ModifyTitle(TNamed* h, const char* centTxt)
{
- if (!centTxt || !h) return;
- h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
+ return;
+ // if (!centTxt || !h) return;
+ // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
}
//__________________________________________________________________
TH1* tracks = FetchResult(list, "tracks");
if (tracks) tracks->SetTitle("ALICE Tracks");
SetAttributes(dndeta, color);
- SetAttributes(dndetaMC, color+2);
+ SetAttributes(dndetaMC, fCentAxis ? color : color+2);
SetAttributes(dndetaTruth,color);
SetAttributes(dndetaSym, color);
- SetAttributes(dndetaMCSym,color+2);
- SetAttributes(tracks, color+3);
+ SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
+ SetAttributes(tracks, fCentAxis ? color : color+3);
+ if (dndetaMC && fCentAxis)
+ dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
+ if (dndetaMCSym && fCentAxis)
+ dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
+ if (dndetaTruth && fCentAxis)
+ dndetaTruth->SetMarkerStyle(34);
ModifyTitle(dndeta, centTxt);
ModifyTitle(dndetaMC, centTxt);
ModifyTitle(dndetaTruth,centTxt);
max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
max = TMath::Max(max, AddHistogram(fResults, tracks));
-
+
+ THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
+ if (rings) {
+ TIter next(rings->GetHists());
+ TH1* hist = 0;
+ while ((hist = static_cast<TH1*>(next())))
+ max = TMath::Max(max, AddHistogram(fResults, hist));
+ }
// Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
// dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
fRatios->Add(Ratio(dndeta, tracks, rmax));
}
+ if (dndetaMC) {
+ fRatios->Add(Ratio(dndeta, dndetaMC, rmax));
+ fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
+ }
if (dndetaTruth) {
fRatios->Add(Ratio(dndeta, dndetaTruth, rmax));
fRatios->Add(Ratio(dndetaSym, dndetaTruth, rmax));
- fRatios->Add(Ratio(dndetaMC, dndetaTruth, rmax));
- fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
}
}
//__________________________________________________________________
c->SetBorderSize(0);
c->SetBorderMode(0);
-#if 1
- Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s",
- y1, y2, y2, (fShowRatios ? "ratios" : ""),
- (fShowLeftRight ? "right/left" : ""));
- Info("Plot", "Maximum is %f", max);
-#endif
PlotResults(max, y1);
c->cd();
Double_t x1, Double_t y1, Double_t x2, Double_t y2)
{
TLegend* l = new TLegend(x1,y1,x2,y2);
- l->SetNColumns(2);
+ l->SetNColumns(fCentAxis ? 1 : 2);
l->SetFillColor(0);
l->SetFillStyle(0);
l->SetBorderSize(0);
l->SetTextFont(132);
TIter next(stack->GetHists());
- TObject* hist = 0;
- while ((hist = next())) {
+ TH1* hist = 0;
+ TObjArray unique;
+ unique.SetOwner();
+ while ((hist = static_cast<TH1*>(next()))) {
TString n(hist->GetTitle());
if (n.Contains("mirrored")) continue;
- l->AddEntry(hist, hist->GetTitle(), "pl");
+ if (unique.FindObject(n.Data())) continue;
+ TObjString* s = new TObjString(n);
+ s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
+ ((hist->GetMarkerColor() & 0xFFFF) << 0));
+ unique.Add(s);
+ // l->AddEntry(hist, hist->GetTitle(), "pl");
}
if (mg) {
TIter nexto(mg->GetListOfGraphs());
- while ((hist = nexto())) {
- TString n(hist->GetTitle());
+ TGraph* g = 0;
+ while ((g = static_cast<TGraph*>(nexto()))) {
+ TString n(g->GetTitle());
if (n.Contains("mirrored")) continue;
- l->AddEntry(hist, hist->GetTitle(), "pl");
+ if (unique.FindObject(n.Data())) continue;
+ TObjString* s = new TObjString(n);
+ s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
+ ((g->GetMarkerColor() & 0xFFFF) << 0));
+ unique.Add(s);
+ // l->AddEntry(hist, hist->GetTitle(), "pl");
}
}
+ unique.ls();
+ TIter nextu(&unique);
+ TObject* s = 0;
+ Int_t i = 0;
+ while ((s = nextu())) {
+ TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
+ s->GetName(), "lp");
+ Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
+ Int_t color = (s->GetUniqueID() >> 0) & 0xFFFF;
+ dd->SetLineColor(kBlack);
+ if (fCentAxis) dd->SetMarkerColor(kBlack);
+ else dd->SetMarkerColor(color);
+ dd->SetMarkerStyle(style);
+ }
TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
d1->SetLineColor(kBlack);
d1->SetMarkerColor(kBlack);
l->Draw();
}
//__________________________________________________________________
+ /**
+ * Build centrality legend
+ *
+ * @param axis Stack to include
+ * @param x1 Lower X coordinate in the range [0,1]
+ * @param y1 Lower Y coordinate in the range [0,1]
+ * @param x2 Upper X coordinate in the range [0,1]
+ * @param y2 Upper Y coordinate in the range [0,1]
+ */
+ void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
+ {
+ if (!fCentAxis) return;
+
+ TLegend* l = new TLegend(x1,y1,x2,y2);
+ l->SetNColumns(1);
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->SetBorderSize(0);
+ l->SetTextFont(132);
+
+ Int_t n = fCentAxis->GetNbins();
+ for (Int_t i = 1; i <= n; i++) {
+ Double_t low = fCentAxis->GetBinLowEdge(i);
+ Double_t upp = fCentAxis->GetBinUpEdge(i);
+ TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
+ Form("%3d%% - %3d%%",
+ int(low), int(upp)), "pl");
+ e->SetMarkerColor(GetCentralityColor(i));
+ }
+ l->Draw();
+ }
+ //__________________________________________________________________
/**
* Plot the results
*
}
// Make a legend in the main result pad
- BuildLegend(fResults, fOthers, .15,p1->GetBottomMargin()+.01,.90,.35);
-#if 0
- TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
- l->SetNColumns(2);
- l->SetFillColor(0);
- l->SetFillStyle(0);
- l->SetBorderSize(0);
- l->SetTextFont(132);
-#endif
+ BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,
+ .35, 1-p1->GetTopMargin()-.01-.1);
+ Double_t x1 = (fCentAxis ? .7 : .15);
+ Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
+ Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01);
+ Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
+
+ BuildLegend(fResults, fOthers, x1, y1, x2, y2);
// Put a title on top
fTitle.ReplaceAll("@", " ");
TString eS;
UShort_t snn = fSNNString->GetUniqueID();
const char* sys = fSysString->GetTitle();
+ if (snn == 2750) snn = 2760;
if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
else eS = Form("%3dGeV", snn);
- TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
+ TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s",
sys,
+ (fCentAxis ? "_{NN}" : ""),
eS.Data(),
+ fCentAxis ? "by centrality" :
fTrigString->GetTitle()));
tt->SetNDC();
tt->SetTextFont(132);
// Make a legend
- TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
- l2->SetNColumns(2);
- l2->SetFillColor(0);
- l2->SetFillStyle(0);
- l2->SetBorderSize(0);
- l2->SetTextFont(132);
+ Double_t xx1 = (fCentAxis ? .7 : .15);
+ Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01 : .90);
+ Double_t yy1 = p3->GetBottomMargin()+.01;
+ Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
+ BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
+ // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
+ // l2->SetNColumns(2);
+ // l2->SetFillColor(0);
+ // l2->SetFillStyle(0);
+ // l2->SetBorderSize(0);
+ // l2->SetTextFont(132);
// Make a nice band from 0.9 to 1.1
TGraphErrors* band = new TGraphErrors(2);
if (!list) return 0;
TH1* ret = static_cast<TH1*>(list->FindObject(name));
+#if 0
if (!ret) {
// all->ls();
Warning("GetResult", "Histogram %s not found", name);
}
+#endif
return ret;
}
// Rebin if needed
Rebin(hist);
- // Info("AddHistogram", "Adding %s to %s",
- // hist->GetName(), stack->GetName());
stack->Add(hist, option);
- // stack->ls();
return hist->GetMaximum();
}
//__________________________________________________________________
sym = Symmetrice(hist);
stack->Add(sym, option);
- // Info("AddHistogram", "Adding %s and %s to %s",
- // hist->GetName(), sym->GetName(), stack->GetName());
return hist->GetMaximum();
}
{
if (!o1 || !o2) return 0;
- TH1* r = 0;
+ TH1* r = 0;
+ const TAttMarker* m1 = 0;
+ const TAttMarker* m2 = 0;
const TH1* h1 = dynamic_cast<const TH1*>(o1);
if (h1) {
- // Info("Ratio", "First is a TH1");
+ m1 = h1;
const TH1* h2 = dynamic_cast<const TH1*>(o2);
if (h2) {
- // Info("Ratio", "Second is a TH1");
- r = RatioHH(h1,h2,max);
+ m2 = h2;
+ r = RatioHH(h1,h2);
}
else {
const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
if (g2) {
- // Info("Ratio", "Second os a TGraph");
- r = RatioHG(h1,g2,max);
+ m2 = g2;
+ r = RatioHG(h1,g2);
}
}
}
else {
const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
if (g1) {
- // Info("Ratio", "First is a TGraphAsymmErrors");
+ m1 = g1;
const TGraphAsymmErrors* g2 =
dynamic_cast<const TGraphAsymmErrors*>(o2);
if (g2) {
- // Info("Ratio", "Second is a TGraphAsymmErrors");
- r = RatioGG(g1, g2, max);
+ m2 = g2;
+ r = RatioGG(g1, g2);
}
}
}
delete r;
r = 0;
}
- // for (Int_t bin = 1; bin <= r->GetNbinsX(); bin++)
- // if (r->GetBinContent(bin) != 0) return r;
-
+ if (r) {
+ r->SetMarkerStyle(m1->GetMarkerStyle());
+ r->SetMarkerColor(m2->GetMarkerColor());
+ r->SetMarkerSize(0.9*m1->GetMarkerSize());
+ r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
+ r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
+ r->SetDirectory(0);
+ max = TMath::Max(RatioMax(r), max);
+ }
+
return r;
}
//__________________________________________________________________
*
* @param h Numerator
* @param g Divisor
- * @param max Maximum diviation from 1
*
* @return h/g
*/
- TH1* RatioHG(const TH1* h, const TGraph* g, Double_t& max) const
+ TH1* RatioHG(const TH1* h, const TGraph* g) const
{
if (!h || !g) return 0;
TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
- ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
- ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
ret->Reset();
- ret->SetMarkerStyle(h->GetMarkerStyle());
- ret->SetMarkerColor(g->GetMarkerColor());
- ret->SetMarkerSize(0.9*h->GetMarkerSize());
- // ret->SetMarkerStyle(g->GetMarkerStyle());
- // ret->SetMarkerColor(h->GetMarkerColor());
- // ret->SetMarkerSize(0.9*g->GetMarkerSize());
- ret->SetDirectory(0);
Double_t xlow = g->GetX()[0];
Double_t xhigh = g->GetX()[g->GetN()-1];
if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
ret->SetBinContent(i, c / f);
ret->SetBinError(i, h->GetBinError(i) / f);
}
- if (ret->GetEntries() > 0)
- max = TMath::Max(RatioMax(ret), max);
return ret;
}
//__________________________________________________________________
*
* @param h1 First histogram (numerator)
* @param h2 Second histogram (denominator)
- * @param max Maximum diviation from 1
*
* @return h1 / h2
*/
- TH1* RatioHH(const TH1* h1, const TH1* h2, Double_t& max) const
+ TH1* RatioHH(const TH1* h1, const TH1* h2) const
{
if (!h1 || !h2) return 0;
- TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
- h1->GetName(),
- h2->GetName())));
- t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
+ TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
t1->Divide(h2);
- // t1->SetMarkerColor(h1->GetMarkerColor());
- // t1->SetMarkerStyle(h2->GetMarkerStyle());
- t1->SetMarkerColor(h2->GetMarkerColor());
- t1->SetMarkerStyle(h1->GetMarkerStyle());
- t1->SetMarkerSize(0.9*h1->GetMarkerSize());
- t1->SetDirectory(0);
- max = TMath::Max(RatioMax(t1), max);
return t1;
}
//__________________________________________________________________
*
* @param g1 Numerator
* @param g2 Denominator
- * @param max Maximum diviation from 1
*
* @return g1 / g2 in a histogram
*/
TH1* RatioGG(const TGraphAsymmErrors* g1,
- const TGraphAsymmErrors* g2, Double_t& max) const
+ const TGraphAsymmErrors* g2) const
{
Int_t nBins = g1->GetN();
TArrayF bins(nBins+1);
if (i == 0) dx = dxi;
else if (dxi != dx) dx = 0;
}
- TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
- TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
TH1* h = 0;
if (dx != 0) {
- h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
+ h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
}
else {
- h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
+ h = new TH1F("tmp", "tmp", nBins, bins.fArray);
}
- h->SetMarkerStyle(g1->GetMarkerStyle());
- h->SetMarkerColor(g2->GetMarkerColor());
- h->SetMarkerSize(0.9*g1->GetMarkerSize());
- // h->SetMarkerStyle(g2->GetMarkerStyle());
- // h->SetMarkerColor(g1->GetMarkerColor());
- // h->SetMarkerSize(0.9*g2->GetMarkerSize());
- h->SetDirectory(0);
Double_t low = g2->GetX()[0];
Double_t high = g2->GetX()[g2->GetN()-1];
h->SetBinContent(i+1, c1 / c2);
h->SetBinError(i+1, e1 / c2);
}
- max = TMath::Max(RatioMax(h), max);
return h;
}
/* @} */
Int_t last = p->fMasterAxis->GetLast();
Double_t x1 = p->fMasterAxis->GetBinCenter(first);
Double_t x2 = p->fMasterAxis->GetBinCenter(last);
- //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
if (p->fSlave1Axis) {
Int_t i1 = p->fSlave1Axis->FindBin(x1);
// Whether to enable low flux specific code
task->SetEnableLowFlux(kFALSE);
+ // Would like to use dynamic cast but CINT interprets that as a
+ // static cast - sigh!
+ Bool_t mc = false;
+ if (task->IsA() == AliForwardMCMultiplicityTask::Class())
+ mc = true;
+
+#if 0
+ if (mc) {
+ AliForwardMCMultiplicityTask* mcTask =
+ static_cast<AliForwardMCMultiplicityTask*>(task);
+ mcTask->SetOnlyPrimary(true);
+ }
+#endif
+ Double_t nXi = mc ? 1 : .5;
+ Bool_t includeSigma = true;
+
// --- Event inspector ---------------------------------------------
// Set the number of SPD tracklets for which we consider the event a
// low flux event
// --- Sharing filter ----------------------------------------------
// Set the low cut used for sharing - overrides settings in eloss fits
- task->GetSharingFilter().SetLowCut(0.3);
+ task->GetSharingFilter().SetLowCut(0.15);
// Set the number of xi's (width of landau peak) to stop at
- task->GetSharingFilter().SetNXi(1);
+ task->GetSharingFilter().SetNXi(nXi);
+ // Set whether or not to include sigma in cut
+ task->GetSharingFilter().SetIncludeSigma(includeSigma);
+ // Enable use of angle corrected signals in the algorithm
+ task->GetSharingFilter().SetUseAngleCorrectedSignals(true);
// --- Density calculator
// Set the maximum number of particle to try to reconstruct
- task->GetDensityCalculator().SetMaxParticles(3);
+ task->GetDensityCalculator().SetMaxParticles(10);
// Wet whether to use poisson statistics to estimate N_ch
task->GetDensityCalculator().SetUsePoisson(false);
// Set the lower multiplicity cut. Overrides setting in energy loss fits.
- task->GetDensityCalculator().SetMultCut(0.3); //was 0.3
+ task->GetDensityCalculator().SetMultCut(-1); //was 0.3
+ // Set the lower per-ring multiplicity cuts
+ task->GetDensityCalculator().SetMultCuts(-1,-1,-1,-1,-1);
+ // USe this many times xi+sigma below MPV
+ task->GetDensityCalculator().SetNXi(nXi);
+ // Set whether or not to include sigma in cut
+ task->GetDensityCalculator().SetIncludeSigma(includeSigma);
+ // Set whether or not to use the phi acceptance
+ task->GetDensityCalculator().SetUsePhiAcceptance(true);
// --- Corrector ---------------------------------------------------
// Whether to use the secondary map correction
// 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->GetEventInspector().SetDebug(4);
+ task->GetSharingFilter().SetDebug(0);
// --- Set limits on fits the energy -------------------------------
// Maximum relative error on parameters
--- /dev/null
+#ifdef BUILD
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliMultiplicity.h"
+#include "AliFMDMCEventInspector.h"
+#include "AliAODForwardMult.h"
+#include "AliLog.h"
+#include <TH1I.h>
+#include <TH2D.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TParameter.h>
+#include <TStopwatch.h>
+#include <TROOT.h>
+#include <THStack.h>
+#include <TStyle.h>
+
+//====================================================================
+/**
+ * Task to evaluate trigger bias in pp
+ *
+ */
+class EvaluateTrigger : public AliAnalysisTaskSE
+{
+public:
+ enum {
+ kNone = 0x0,
+ kESD = 0x1,
+ kMC = 0x2
+ };
+ /**
+ * Constructor
+ */
+ EvaluateTrigger()
+ : AliAnalysisTaskSE(),
+ fInel(),
+ fInelGt0(),
+ fNSD(),
+ fInspector(),
+ fFirstEvent(true),
+ fData(0),
+ fTriggers(0),
+ fTrackletRequirement(kESD),
+ fVertexRequirement(kESD),
+ fVertexAxis(0, 0, 0),
+ fVertexESD(0),
+ fVertexMC(0),
+ fM(0)
+ {}
+ /**
+ * Constructor
+ */
+ EvaluateTrigger(const char* /*name*/)
+ : AliAnalysisTaskSE("evaluateTrigger"),
+ fInel(AliAODForwardMult::kInel),
+ fInelGt0(AliAODForwardMult::kInelGt0),
+ fNSD(AliAODForwardMult::kNSD),
+ fInspector("eventInspector"),
+ fFirstEvent(true),
+ fData(0),
+ fTriggers(0),
+ fTrackletRequirement(kESD),
+ fVertexRequirement(kESD),
+ fVertexAxis(10, -10, 10),
+ fVertexESD(0),
+ fVertexMC(0),
+ fM(0)
+ {
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
+ }
+ void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
+ void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
+ void SetVertexAxis(Int_t nBins, Double_t low, Double_t high)
+ {
+ fVertexAxis.Set(nBins, low, high);
+ }
+ /**
+ * Intialize
+ */
+ void Init() {}
+ /**
+ * Create output objects
+ */
+ void UserCreateOutputObjects()
+ {
+ fList = new TList;
+ fList->SetOwner();
+ fList->SetName(GetName());
+
+ Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 1000 };
+ Int_t nM = 10;
+ TAxis mAxis(nM, mb);
+ TAxis eAxis(200, -4, 6);
+ TAxis pAxis(40, 0, 2*TMath::Pi());
+
+ fData = new TH2D("data", "Cache",
+ eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
+ pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
+ fData->SetDirectory(0);
+ fData->SetXTitle("#eta");
+ fData->SetYTitle("#varphi [radians]");
+ fData->SetZTitle("N_{ch}(#eta,#varphi)");
+ fData->Sumw2();
+
+ fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}",nM+1, -0.5, nM+.5);
+ fM->SetXTitle("N_{ch}|_{|#eta|<1}");
+ fM->SetYTitle("Events");
+ fM->SetFillColor(kRed+1);
+ fM->SetFillStyle(3001);
+ fM->SetDirectory(0);
+ fList->Add(fM);
+
+ for (Int_t i = 0; i <= nM; i++) {
+ TString lbl;
+ if (i == 0) lbl = "all";
+ else if (i == nM) lbl = Form("%d+",int(mAxis.GetBinLowEdge(i)+.5));
+ else lbl = Form("<%d",int(mAxis.GetBinUpEdge(i)+.5));
+ fM->GetXaxis()->SetBinLabel(i+1, lbl);
+ }
+
+ fTriggers = new TH1I("triggers", "Triggers", 6, -.5, 5.5);
+ fTriggers->SetDirectory(0);
+ fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
+ fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
+ fTriggers->GetXaxis()->SetBinLabel(3, "INEL>0 (MC)");
+ fTriggers->GetXaxis()->SetBinLabel(4, "INEL>0 (ESD)");
+ fTriggers->GetXaxis()->SetBinLabel(5, "NSD (MC)");
+ fTriggers->GetXaxis()->SetBinLabel(6, "NSD (ESD)");
+ fTriggers->SetFillColor(kYellow+1);
+ fTriggers->SetFillStyle(3001);
+ fList->Add(fTriggers);
+
+ fVertexESD = new TH1D("vertexESD", "ESD vertex distribution",
+ fVertexAxis.GetNbins(),
+ fVertexAxis.GetXmin(),
+ fVertexAxis.GetXmax());
+ fVertexESD->SetDirectory(0);
+ fVertexESD->SetFillColor(kRed+1);
+ fVertexESD->SetFillStyle(3001);
+ fVertexESD->SetXTitle("v_{z} [cm]");
+ fVertexESD->SetYTitle("P(v_{z})");
+ fList->Add(fVertexESD);
+
+ fVertexMC = new TH1D("vertexMC", "MC vertex distribution",
+ fVertexAxis.GetNbins(),
+ fVertexAxis.GetXmin(),
+ fVertexAxis.GetXmax());
+ fVertexMC->SetDirectory(0);
+ fVertexMC->SetFillColor(kBlue+1);
+ fVertexMC->SetFillStyle(3001);
+ fVertexMC->SetXTitle("v_{z} [cm]");
+ fVertexMC->SetYTitle("P(v_{z})");
+ fList->Add(fVertexMC);
+
+ fInel.CreateObjects(fList, fM, fData);
+ fInelGt0.CreateObjects(fList, fM, fData);
+ fNSD.CreateObjects(fList, fM, fData);
+
+
+ fInspector.DefineOutput(fList);
+ fInspector.Init(fVertexAxis);
+
+ PostData(1, fList);
+ }
+ /**
+ * Event processing
+ */
+ void UserExec(Option_t*)
+ {
+ // Get the input data - MC event
+ AliMCEvent* mcEvent = MCEvent();
+ if (!mcEvent) {
+ AliWarning("No MC event found");
+ return;
+ }
+
+ // Get the input data - ESD event
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ if (!esd) {
+ AliWarning("No ESD event found for input event");
+ return;
+ }
+
+ if (fFirstEvent && esd->GetESDRun()) {
+ fInspector.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()));
+
+ fFirstEvent = false;
+ }
+
+ // Get the particle stack
+ AliStack* stack = mcEvent->Stack();
+
+ // Some variables
+ UInt_t triggers; // Trigger bits
+ Bool_t lowFlux; // Low flux flag
+ UShort_t iVz; // Vertex bin from ESD
+ Double_t vZ; // Z coordinate from ESD
+ Double_t cent; // Centrality
+ UShort_t iVzMc; // Vertex bin from MC
+ Double_t vZMc; // Z coordinate of IP vertex from MC
+ Double_t b; // Impact parameter
+ Int_t nPart; // Number of participants
+ Int_t nBin; // Number of binary collisions
+ Double_t phiR; // Reaction plane from MC
+
+ // Process the data
+ Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent);
+ Int_t retMC = fInspector.ProcessMC(mcEvent, triggers, iVzMc,
+ vZMc, b, nPart, nBin, phiR);
+
+ Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
+ Bool_t hasMCVtx = retMC == AliFMDEventInspector::kOk;
+ if (hasESDVtx) fVertexESD->Fill(vZ);
+ if (hasMCVtx) fVertexMC->Fill(vZMc);
+
+ Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
+ Bool_t isMcNSD = (triggers & AliAODForwardMult::kMCNSD);
+
+ Int_t mESD = 0;
+ const AliMultiplicity* spdmult = esd->GetMultiplicity();
+ if (!spdmult) {
+ AliWarning("No SPD multiplicity");
+ }
+ else {
+ // Check if we have one or more tracklets
+ // in the range -1 < eta < 1 to set the INEL>0
+ // trigger flag.
+ Int_t n = spdmult->GetNumberOfTracklets();
+ for (Int_t j = 0; j < n; j++)
+ if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
+ }
+
+ // Reset cache
+ fData->Reset();
+ Int_t mMC = 0; // Number of particles in |eta|<1
+
+ // Loop over all tracks
+ Int_t nTracks = mcEvent->GetNumberOfTracks();
+ for (Int_t iTr = 0; iTr < nTracks; iTr++) {
+ AliMCParticle* particle =
+ static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
+
+ // Check the returned particle
+ if (!particle) continue;
+
+ // Check if this charged and a primary
+ Bool_t isCharged = particle->Charge() != 0;
+ Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
+
+ if (!isCharged || !isPrimary) continue;
+
+
+ // Fill (eta,phi) of the particle into histograsm for b
+ Double_t eta = particle->Eta();
+ Double_t phi = particle->Phi();
+
+ fData->Fill(eta, phi);
+ if (TMath::Abs(eta) <= 1) mMC++;
+ }
+ Int_t m = mESD;
+ if (fTrackletRequirement == kMC) m = mMC;
+ fM->Fill(m);
+
+ bool isMcInelGt0 = isMcInel && (mMC > 0);
+
+ bool hasVertex = true;
+ if (fVertexRequirement & kMC) hasVertex = hasVertex && hasMCVtx;
+ if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
+
+ if (isMcInel) {
+ fTriggers->Fill(0);
+ bool triggered = (triggers & AliAODForwardMult::kInel);
+ if (triggered) fTriggers->Fill(1);
+ fInel.AddEvent(triggered, hasVertex, m, fData);
+ }
+ if (isMcInelGt0) {
+ fTriggers->Fill(2);
+ bool triggered = (triggers & AliAODForwardMult::kInelGt0);
+ if (triggered) fTriggers->Fill(3);
+ fInelGt0.AddEvent(triggered, hasVertex, m, fData);
+ }
+ if (isMcNSD) {
+ fTriggers->Fill(4);
+ bool triggered = (triggers & AliAODForwardMult::kNSD);
+ if (triggered) fTriggers->Fill(5);
+ fNSD.AddEvent(triggered, hasVertex, m, fData);
+ }
+ PostData(1, fList);
+ }
+ /**
+ * End of job processing
+ */
+ void Terminate(Option_t*)
+ {
+ fList = dynamic_cast<TList*>(GetOutputData(1));
+ if (!fList) {
+ AliError(Form("No output list defined (%p)", GetOutputData(1)));
+ if (GetOutputData(1)) GetOutputData(1)->Print();
+ return;
+ }
+
+
+ TList* output = new TList;
+ output->SetName(GetName());
+ output->SetOwner();
+
+ fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
+ fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
+ fM = static_cast<TH1D*>(fList->FindObject("m"));
+ if (fVertexMC) {
+ TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
+ vtxMC->SetDirectory(0);
+ vtxMC->Scale(1. / vtxMC->GetEntries());
+ output->Add(vtxMC);
+ }
+ if (fVertexESD) {
+ TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
+ vtxESD->SetDirectory(0);
+ vtxESD->Scale(1. / vtxESD->GetEntries());
+ output->Add(vtxESD);
+ }
+ if (fM) {
+ TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
+ m->SetDirectory(0);
+ m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
+ m->Scale(1. / m->GetBinContent(1));
+ output->Add(m);
+ }
+
+ fInel.Finish(fList, output);
+ fInelGt0.Finish(fList, output);
+ fNSD.Finish(fList, output);
+
+ PostData(2, output);
+ }
+
+protected:
+ //__________________________________________________________________
+ /**
+ * Structure to hold per trigger type information
+ */
+ struct TriggerType : public TNamed
+ {
+ //________________________________________________________________
+ /**
+ * Structure to hold per multiplicity bin information
+ */
+ struct MBin : public TNamed
+ {
+ TH2D* fTruth;
+ TH2D* fTriggered;
+ TH2D* fAccepted;
+ TH1I* fCounts;
+ UShort_t fLow;
+ UShort_t fHigh;
+ /**
+ * Constructor
+ */
+ MBin()
+ : fTruth(0), fTriggered(0), fAccepted(0),
+ fCounts(0), fLow(0), fHigh(1000) {}
+ /**
+ * Constructor
+ *
+ * @param p Parent list
+ * @param low Low cut
+ * @param high High cut
+ * @param eAxis Eta axis
+ * @param pAxis Phi axis
+ */
+ MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist)
+ : fTruth(0),
+ fTriggered(0),
+ fAccepted(0),
+ fCounts(0),
+ fLow(low),
+ fHigh(high)
+ {
+ if (low >= high) SetName("all");
+ else SetName(Form("m%03d_%03d", fLow, fHigh));
+ TList* l = new TList;
+ l->SetOwner();
+ l->SetName(GetName());
+ p->Add(l);
+
+ fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
+ fTruth->SetTitle("MC truth");
+ fTruth->SetDirectory(0);
+ fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
+ fTruth->Sumw2();
+ fTruth->Reset();
+ l->Add(fTruth);
+
+ fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
+ fTriggered->SetTitle("Triggered");
+ fTriggered->SetDirectory(0);
+ fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
+ fTriggered->Sumw2();
+ fTriggered->Reset();
+ l->Add(fTriggered);
+
+ fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
+ fAccepted->SetTitle("Accepted");
+ fAccepted->SetDirectory(0);
+ fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
+ fAccepted->Sumw2();
+ fAccepted->Reset();
+ l->Add(fAccepted);
+
+ fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
+ fCounts->SetDirectory(0);
+ fCounts->GetXaxis()->SetBinLabel(1, "Truth");
+ fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
+ fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
+ fCounts->SetYTitle("# events");
+ l->Add(fCounts);
+ }
+ /**
+ * Add event observation
+ *
+ * @param triggered Whether the event was triggered
+ * @param event Data for this event
+ */
+ void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event)
+ {
+ fCounts->Fill(0);
+ fTruth->Add(event);
+ if (triggered) {
+ fCounts->Fill(1);
+ fTriggered->Add(event);
+ if (hasVtx) {
+ fCounts->Fill(2);
+ fAccepted->Add(event);
+ }
+ }
+ }
+ /**
+ * End of job processing
+ *
+ * @param p Parent list
+ * @param o Output parent list
+ * @param stack Stack of histograms
+ *
+ * @return Trigger efficiency
+ */
+ Double_t Finish(const TList* p, TList* o, THStack* stack)
+ {
+ TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
+ if (!l) {
+ Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
+ return 0;
+ }
+ fTruth = static_cast<TH2D*>(l->FindObject("truth"));
+ fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
+ fAccepted = static_cast<TH2D*>(l->FindObject("accepted"));
+ fCounts = static_cast<TH1I*>(l->FindObject("counts"));
+
+ Int_t nTruth = fCounts->GetBinContent(1);
+ Int_t nTriggered = fCounts->GetBinContent(2);
+ Int_t nAccepted = fCounts->GetBinContent(3);
+ Double_t eff = 0;
+ if (nTruth > 0) eff = double(nTriggered) / nTruth;
+ else if (nTriggered == nTruth) eff = 1;
+
+ if (nTruth > 0) fTruth->Scale(1. / nTruth);
+ if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
+ if (nAccepted > 0) fAccepted->Scale(1. / nAccepted);
+
+ if (fLow >= fHigh)
+ Info("Finish", "%-6s [all] E_X=N_T/N_X=%9d/%-9d=%f "
+ "E_V=N_A/N_T=%9d/%-9d=%f",
+ p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered,
+ (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
+ else
+ Info("Finish", "%-6s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
+ "E_V=N_A/N_T=%9d/%-9d=%f",
+ p->GetName(), fLow, fHigh, nTriggered, nTruth, eff,
+ nAccepted, nTriggered,
+ (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
+
+ TList* out = new TList;
+ out->SetName(GetName());
+ out->SetOwner();
+ o->Add(out);
+
+ out->Add(fTruth);
+ out->Add(fTriggered);
+ out->Add(new TParameter<double>("eff", eff));
+
+ TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
+ bias->Divide(fTruth);
+ bias->SetDirectory(0);
+ bias->SetZTitle("Trigger bias (accepted/truth)");
+ out->Add(bias);
+
+ TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
+ truth_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
+ fLow, fHigh));
+ truth_px->Scale(1. / fTruth->GetNbinsY());
+ truth_px->SetDirectory(0);
+ truth_px->SetLineColor(kRed+1);
+ truth_px->SetMarkerColor(kRed+1);
+ truth_px->SetFillColor(kRed+1);
+ truth_px->SetFillStyle(0);
+ out->Add(truth_px);
+
+ TH1D* triggered_px =
+ static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
+ triggered_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
+ fLow, fHigh));
+ triggered_px->Scale(1. / fTriggered->GetNbinsY());
+ triggered_px->SetDirectory(0);
+ triggered_px->SetLineColor(kGreen+1);
+ triggered_px->SetMarkerColor(kGreen+1);
+ triggered_px->SetFillColor(kGreen+1);
+ triggered_px->SetFillStyle(0);
+ out->Add(triggered_px);
+
+ TH1D* accepted_px =
+ static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
+ accepted_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
+ fLow, fHigh));
+ accepted_px->Scale(1. / fAccepted->GetNbinsY());
+ accepted_px->SetLineColor(kBlue+1);
+ accepted_px->SetMarkerColor(kBlue+1);
+ accepted_px->SetFillColor(kBlue+1);
+ accepted_px->SetDirectory(0);
+ out->Add(accepted_px);
+
+ THStack* data = new THStack("data", "Data distributions");
+ data->Add(truth_px);
+ data->Add(triggered_px);
+ data->Add(accepted_px);
+ out->Add(data);
+
+#if 0
+ TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
+ bias_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
+ fLow, fHigh));
+ bias_px->Scale(1. / bias->GetNbinsY());
+#else
+ TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
+ bias_px->Divide(truth_px);
+ bias_px->SetYTitle("Trigger bias (triggered/truth)");
+#endif
+ bias_px->SetDirectory(0);
+ bias_px->SetMarkerStyle(20);
+ bias_px->SetFillStyle(0);
+ bias_px->SetMinimum(0);
+ out->Add(bias_px);
+
+ stack->Add(bias_px);
+
+ return eff;
+ }
+ ClassDef(MBin,1);
+ };
+
+ /**
+ * Constructor
+ */
+ TriggerType()
+ : TNamed(),
+ fMask(0),
+ fM(0),
+ fBins(0)
+ {}
+ //--- Constructor ------------------------------------------------
+ /**
+ * Constructor
+ *
+ * @param mask Trigger mask
+ */
+ TriggerType(UShort_t mask)
+ : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
+ fMask(mask),
+ fM(0),
+ fBins(0)
+ {
+ }
+ /**
+ * Create objects
+ *
+ * @param list PArent list
+ * @param mAxis Multiplicity axis
+ * @param eAxis Eta axis
+ * @param pAxis Phi axis
+ */
+ void CreateObjects(TList* list,
+ const TH1D* mHist,
+ const TH2D* dHist)
+ {
+ TList* ours = new TList;
+ ours->SetName(GetName());
+ ours->SetOwner();
+ list->Add(ours);
+
+ fM = static_cast<TH1D*>(mHist->Clone("m"));
+ fM->SetDirectory(0);
+ ours->Add(fM);
+
+ fBins = new TObjArray;
+ fBins->AddAt(new MBin(ours, 0, 0, dHist), 0);
+ for (Int_t i = 1; i <= fM->GetNbinsX(); i++) {
+ Double_t low = fM->GetXaxis()->GetBinLowEdge(i);
+ Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
+
+ fBins->AddAt(new MBin(ours, low, high, dHist), i);
+ }
+ }
+ /**
+ * Find bin corresponding to m
+ *
+ * @param m Multiplicity
+ *
+ * @return Bin.
+ */
+ MBin* FindBin(UShort_t m)
+ {
+ Int_t b = fM->GetXaxis()->FindBin(m);
+ if (b <= 0) return 0;
+ if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
+ return static_cast<MBin*>(fBins->At(b));
+ }
+ /**
+ * Add event observation
+ *
+ * @param triggered IF this is triggered
+ * @param m Multiplicity
+ * @param data Observation
+ */
+ void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
+ {
+ fM->AddBinContent(1);
+ fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
+
+ MBin* all = static_cast<MBin*>(fBins->At(0));
+ all->AddEvent(triggered, hasVtx, data);
+
+ MBin* bin = FindBin(m);
+ bin->AddEvent(triggered, hasVtx, data);
+ }
+ /**
+ * End of job processing
+ *
+ * @param p Parent list
+ * @param o Parent output list
+ */
+ void Finish(const TList* p, TList* o)
+ {
+ TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
+ if (!l) {
+ Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
+ return;
+ }
+
+ TList* ours = new TList;
+ ours->SetName(GetName());
+ ours->SetOwner();
+ o->Add(ours);
+
+ fM = static_cast<TH1D*>(l->FindObject("m"));
+ if (!fM) {
+ Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
+ return;
+ }
+ TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
+ m->SetDirectory(0);
+ m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
+ if (m->GetBinContent(1) > 0)
+ m->Scale(1. / m->GetBinContent(1));
+ ours->Add(m);
+
+ Int_t nBin = fM->GetNbinsX();
+ TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
+ effs->SetYTitle("#epsilon_{X}");
+ effs->SetFillColor(kRed+1);
+ effs->SetDirectory(0);
+ effs->SetMinimum(0);
+
+ gStyle->SetPalette(1);
+ Int_t nCol = gStyle->GetNumberOfColors();
+ THStack* stack = new THStack("biases", "Trigger biases");
+ for (Int_t i = 0; i <= nBin; i++) {
+ MBin* bin = static_cast<MBin*>(fBins->At(i));
+ effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
+ TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
+ Int_t col = kBlack;
+ if (i != 0) {
+ Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
+ col = gStyle->GetColorPalette(icol);
+ }
+ h->SetMarkerColor(col);
+ h->SetFillColor(col);
+ h->SetLineColor(col);
+ }
+
+ ours->Add(stack);
+ ours->Add(effs);
+ }
+ UShort_t fMask;
+ TH1D* fM;
+ TObjArray* fBins;
+ ClassDef(TriggerType,1);
+ };
+ TriggerType fInel;
+ TriggerType fInelGt0;
+ TriggerType fNSD;
+ AliFMDMCEventInspector fInspector;
+ TList* fList;
+ Bool_t fFirstEvent;
+ TH2D* fData;
+ TH1I* fTriggers;
+ UInt_t fTrackletRequirement;
+ UInt_t fVertexRequirement;
+ TAxis fVertexAxis;
+ TH1D* fVertexESD;
+ TH1D* fVertexMC;
+ TH1D* fM;
+ ClassDef(EvaluateTrigger,1);
+};
+#else
+//====================================================================
+void MakeEvaluateTriggers(const char* esddir,
+ Int_t nEvents = -1,
+ UInt_t vtx = 0x1,
+ UInt_t trk = 0x1,
+ UInt_t vz = 10,
+ Int_t proof = 0)
+{
+ // --- Libraries to load -------------------------------------------
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+ // --- Check for proof mode, and possibly upload pars --------------
+ if (proof> 0) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
+ if (!LoadPars(proof)) {
+ Error("MakeAOD", "Failed to load PARs");
+ return;
+ }
+ }
+
+ // --- Our data chain ----------------------------------------------
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
+ TChain* chain = MakeChain("ESD", esddir,true);
+ // If 0 or less events is select, choose all
+ if (nEvents <= 0) nEvents = chain->GetEntries();
+
+ // --- Set the macro path ------------------------------------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+ "$ALICE_ROOT/ANALYSIS/macros",
+ gROOT->GetMacroPath()));
+
+ // --- Creating the manager and handlers ---------------------------
+ AliAnalysisManager *mgr = new AliAnalysisManager("Triggers",
+ "Forward multiplicity");
+ AliAnalysisManager::SetCommonFileName("triggers.root");
+
+ // --- ESD input handler -------------------------------------------
+ AliESDInputHandler *esdHandler = new AliESDInputHandler();
+ mgr->SetInputEventHandler(esdHandler);
+
+ // --- Monte Carlo handler -----------------------------------------
+ AliMCEventHandler* mcHandler = new AliMCEventHandler();
+ mgr->SetMCtruthEventHandler(mcHandler);
+ mcHandler->SetReadTR(true);
+
+ // --- Add tasks ---------------------------------------------------
+ // Physics selection
+ gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));
+
+#if 0
+ // --- Fix up physics selection to give proper A,C, and E triggers -
+ AliInputEventHandler* ih =
+ static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+ AliPhysicsSelection* ps =
+ static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+ // Ignore trigger class when selecting events. This mean that we
+ // get offline+(A,C,E) events too
+ ps->SetSkipTriggerClassSelection(true);
+#endif
+
+ // --- compile our code --------------------------------------------
+ gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
+ "-I${ALICE_ROOT}/ANALYSIS "
+ "-I${ALICE_ROOT}/include -DBUILD=1");
+ gROOT->LoadMacro("${ALICE_ROOT}/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
+
+ // --- Make our object ---------------------------------------------
+ EvaluateTrigger* task = new EvaluateTrigger("triggers");
+ mgr->AddTask(task);
+ task->SetVertexRequirement(vtx);
+ task->SetTrackletRequirement(trk);
+ task->SetVertexAxis(10, -vz, vz);
+
+ // --- create containers for input/output --------------------------
+ AliAnalysisDataContainer *sums =
+ mgr->CreateContainer("triggerSums", TList::Class(),
+ AliAnalysisManager::kOutputContainer,
+ AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *output =
+ mgr->CreateContainer("triggerResults", TList::Class(),
+ AliAnalysisManager::kParamContainer,
+ AliAnalysisManager::GetCommonFileName());
+
+ // --- connect input/output ----------------------------------------
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, sums);
+ mgr->ConnectOutput(task, 2, output);
+
+ // --- Run the analysis --------------------------------------------
+ TStopwatch t;
+ if (!mgr->InitAnalysis()) {
+ Error("MakeAOD", "Failed to initialize analysis train!");
+ return;
+ }
+ // Skip terminate if we're so requested and not in Proof or full mode
+ mgr->SetSkipTerminate(false);
+ // Some informative output
+ mgr->PrintStatus();
+ if (proof) mgr->SetDebugLevel(3);
+ if (mgr->GetDebugLevel() < 1 && !proof)
+ mgr->SetUseProgressBar(kTRUE,100);
+
+ // Run the train
+ t.Start();
+ Printf("=== RUNNING ANALYSIS on %9d events ==========================",
+ nEvents);
+ mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+ t.Stop();
+ t.Print();
+}
+#endif
mp->Add(AliceCentralInelGt7000());
}
}
+#if 0
else
Warning("GetData", "No other results for sys=%d, energy=%d",
sys, energy);
+#endif
}
else if (sys == 2) {
// Nothing for PbPb so far
en = Form(", #sqrt{s_{NN}}=%dGeV", energy);
else
en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000);
- Warning("GetData", "No other data for PbP b yet");
+ // Warning("GetData", "No other data for PbPb yet");
}
else
Warning("GetData", "Unknown system %d", sys);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff