fDebug(false)
{
// Default constructor
- DGUARD(0,0,"MC track density default construction");
+ DGUARD(fDebug, 3,"Default CTOR of AliBasMCTrackDensity");
}
//____________________________________________________________________
fDebug(false)
{
// Normal constructor constructor
- DGUARD(0,0,"MC track density named construction: %s", name);
+ DGUARD(fDebug, 3,"Named CTOR of AliBasMCTrackDensity: %s", name);
}
//____________________________________________________________________
fDebug(o.fDebug)
{
// Normal constructor constructor
- DGUARD(0,0,"MC track density copy construction");
+ DGUARD(fDebug, 3,"Copy CTOR of AliBasMCTrackDensity");
}
//____________________________________________________________________
//
// Constructor
//
- DGUARD(0,0,"Default construction of AliBasedNdetaTask");
+ DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask");
}
//____________________________________________________________________
//
// Constructor
//
- DGUARD(0,0,"Named construction of AliBasedNdetaTask: %s", name);
+ DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name);
fListOfCentralities = new TObjArray(1);
// Set the normalisation scheme
fTriggerString(o.fTriggerString),
fFinalMCCorrFile(o.fFinalMCCorrFile)
{
- DGUARD(0,0,"Copy construction of AliBasedNdetaTask");
+ DGUARD(fDebug, 3,"Copy CTOR of AliBasedNdetaTask");
}
//____________________________________________________________________
}
}
}
+//________________________________________________________________________
+void
+AliBasedNdetaTask::ScaleToCoverage(TH1D* copy, const TH1D* norm)
+{
+ // Normalize to the acceptance -
+ // dndeta->Divide(accNorm);
+ for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
+ Double_t a = norm->GetBinContent(i);
+ if (a <= 0) {
+ copy->SetBinContent(i,0);
+ copy->SetBinError(i,0);
+ continue;
+ }
+ Double_t c = copy->GetBinContent(i);
+ Double_t e = copy->GetBinError(i);
+ copy->SetBinContent(i, c / a);
+ copy->SetBinError(i, e / a);
+ }
+}
//________________________________________________________________________
TH1D*
Int_t first = firstbin;
Int_t last = lastbin;
- if (first < 0) first = 0;
- else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
+ if (first < 0) first = 1;
+ else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1;
if (last < 0) last = yaxis->GetNbins();
- else if (last > yaxis->GetNbins()+1) last = yaxis->GetNbins();
+ else if (last >= yaxis->GetNbins()+2) last = yaxis->GetNbins()+1;
if (last-first < 0) {
AliWarningGeneral("AliBasedNdetaTask",
Form("Nothing to project [%d,%d]", first, last));
retCopy->SetMarkerStyle(marker);
retCopy->SetDirectory(0);
- TH1D* norm = ProjectX(fSum, "norm", 0, 0, rootProj, corrEmpty, false);
- TH1D* norm0 = ProjectX(fSum0, "norm0", 0, 0, rootProj, corrEmpty, false);
- TH1D* normAll = ProjectX(ret, "normAll", 0, 0, rootProj, corrEmpty, false);
+ Int_t nY = fSum->GetNbinsY();
+ Int_t o = 0; // nY+1;
+ TH1D* norm = ProjectX(fSum, "norm", o, o, rootProj, corrEmpty, false);
+ TH1D* norm0 = ProjectX(fSum0, "norm0", o, o, rootProj, corrEmpty, false);
+ TH1D* normAll = ProjectX(ret, "normAll", o, o, rootProj, corrEmpty, false);
norm->SetDirectory(0);
norm0->SetDirectory(0);
normAll->SetDirectory(0);
ScaleToCoverage(sum0Copy, norm0);
ScaleToCoverage(retCopy, normAll);
- Int_t nY = fSum->GetNbinsY();
TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1, nY,rootProj,corrEmpty);
TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1, nY,rootProj,corrEmpty);
TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1, nY,rootProj,corrEmpty);
sum0CopyPx->SetDirectory(0);
retCopyPx->SetDirectory(0);
+ TH1D* phi = ProjectX(fSum, "phi", nY+1, nY+1,rootProj,corrEmpty);
+ TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1, nY+1,rootProj,corrEmpty);
+ TH1D* phiAll = ProjectX(ret, "phiAll", nY+1, nY+1,rootProj,corrEmpty);
+ phi->SetDirectory(0);
+ phi0->SetDirectory(0);
+ phiAll->SetDirectory(0);
+
// Scale our 1D histograms
sumCopyPx->Scale(1., "width");
sum0CopyPx->Scale(1., "width");
out->Add(norm);
out->Add(norm0);
out->Add(normAll);
+ out->Add(phi);
+ out->Add(phi0);
+ out->Add(phiAll);
AliInfoF("Returning (1/%f * %s + 1/%f * %s), "
"1/%f * %d + 1/%f * %d = %d",
//
// Constructor
//
- DGUARD(0,0,"Centrality bin default construction");
+ DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin");
}
//____________________________________________________________________
AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
// low Lower centrality cut in percent
// high Upper centrality cut in percent
//
- DGUARD(0,0,"Named Centrality bin construction: %s [%3d,%3d]",name,low,high);
+ DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: %s [%3d,%3d]",
+ name,low,high);
if (low <= 0 && high <= 0) {
fLow = 0;
fHigh = 0;
// Parameters:
// other Object to copy from
//
- DGUARD(0,0,"Copy Centrality bin construction");
+ DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin");
}
//____________________________________________________________________
AliBasedNdetaTask::CentralityBin::~CentralityBin()
//
// Destructor
//
- DGUARD(fDebug,3,"Centrality bin desctruction");
- if (fSums) fSums->Delete();
- if (fOutput) fOutput->Delete();
+ DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin");
+
+ // if (fSums) fSums->Delete();
+ // if (fOutput) fOutput->Delete();
}
//____________________________________________________________________
DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
TH2D* copy = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
GetName(), postfix)));
- TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0,
+ Int_t nY = sum->GetNbinsY();
+ Int_t o = (corrEmpty ? 0 : nY+1);
+ TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), o, o,
rootProj, corrEmpty, false);
accNorm->SetDirectory(0);
else AliInfo("No shape correction specified, or disabled");
// --- Normalize to the coverage -----------------------------------
- ScaleToCoverage(copy, accNorm);
-
- // --- Event-level normalization -----------------------------------
- copy->Scale(scaler);
+ if (corrEmpty) {
+ ScaleToCoverage(copy, accNorm);
+ // --- Event-level normalization ---------------------------------
+ copy->Scale(scaler);
+ }
// --- Project on X axis -------------------------------------------
TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
- 1, sum->GetNbinsY(), rootProj, corrEmpty);
+ 1, nY, rootProj, corrEmpty);
dndeta->SetDirectory(0);
// Event-level normalization
+ if (!corrEmpty) {
+ ScaleToCoverage(dndeta, accNorm);
+ dndeta->Scale(scaler);
+ }
dndeta->Scale(1., "width");
copy->Scale(1., "width");
* @param norm Coverage histogram
*/
static void ScaleToCoverage(TH2D* copy, const TH1D* norm);
+ /**
+ * Scale the copy of the 1D histogram by coverage in supplied 1D histogram
+ *
+ * @param copy Data to scale
+ * @param norm Coverage histogram
+ */
+ static void ScaleToCoverage(TH1D* copy, const TH1D* norm);
/**
* Set histogram graphical options, etc.
*
// Parameters:
// name Name of task
//
- DGUARD(fDebug,0,"Default construction of AliCentralMCCorrectionsTask");
+ DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
}
//____________________________________________________________________
// Parameters:
// name Name of task
//
- DGUARD(fDebug,0,"Named construction of AliCentralMCCorrectionsTask: %s",name);
+ DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class());
fBranchNames =
// Parameters:
// o Object to copy from
//
- DGUARD(fDebug,0,"Copy construction of AliCentralMCCorrectionsTask");
+ DGUARD(fDebug, 3,"Copy CTOR of AliCentralMCCorrectionsTask");
}
//____________________________________________________________________
UInt_t triggers = 0; // Trigger bits
Bool_t lowFlux = true; // Low flux flag
UShort_t iVz = 0; // Vertex bin from ESD
- Double_t vZ = 1000; // Z coordinate from ESD
+ TVector3 ip; // Z coordinate from ESD
Double_t cent = -1; // Centrality
UShort_t iVzMc = 0; // Vertex bin from MC
Double_t vZMc = 1000; // Z coordinate of IP vertex from MC
Double_t phiR = 1000; // Reaction plane from MC
UShort_t nClusters = 0; // Number of SPD clusters
// Process the data
- UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
+ UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip,
cent, nClusters);
fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
b, cMC, nPart, nBin, phiR);
//
// Constructor
//
- DGUARD(fDebug,0,"Named CTOR of AliCentralMultiplicityTask: %s", name);
+ DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
DefineOutput(1, TList::Class());
fBranchNames =
"ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
//
// Constructor
//
- DGUARD(fDebug,0,"Default CTOR of AliCentralMultiplicityTask");
+ DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
}
//____________________________________________________________________
AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
//
// Copy constructor
//
- DGUARD(fDebug,0,"COPY CTOR of AliCentralMultiplicityTask");
+ DGUARD(fDebug, 3,"COPY CTOR of AliCentralMultiplicityTask");
}
//____________________________________________________________________
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip;
Double_t cent = -1;
UShort_t nClusters = 0;
UInt_t found = fInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ ivz, ip, cent, nClusters);
// No event or no trigger
if (found & AliFMDEventInspector::kNoEvent) return;
if (!hAcceptance) AliFatal("No acceptance!");
if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
-
- for(Int_t nx = 1; nx <= aodHist.GetNbinsX(); nx++) {
- Float_t accCor = hAcceptance->GetBinContent(nx);
- Float_t accErr = hAcceptance->GetBinError(nx);
+
+ Int_t nY = aodHist.GetNbinsY();
+ for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) {
+ Float_t accCor = hAcceptance->GetBinContent(ix);
+ Float_t accErr = hAcceptance->GetBinError(ix);
Bool_t fiducial = true;
- if (nx < fEtaMin[vtxbin-1] || nx > fEtaMax[vtxbin-1])
+ if (ix < fEtaMin[vtxbin-1] || ix > fEtaMax[vtxbin-1])
fiducial = false;
// Bool_t etabinSeen = kFALSE;
- for(Int_t ny = 1; ny <= aodHist.GetNbinsY(); ny++) {
+ for(Int_t iy = 1; iy <= nY; iy++) {
#if 1
if (!fiducial) {
- aodHist.SetBinContent(nx, ny, 0);
- aodHist.SetBinError(nx, ny, 0);
+ aodHist.SetBinContent(ix, iy, 0);
+ aodHist.SetBinError(ix, iy, 0);
continue;
}
#endif
// Get currrent value
- Float_t aodValue = aodHist.GetBinContent(nx,ny);
- Float_t aodErr = aodHist.GetBinError(nx,ny);
+ Float_t aodValue = aodHist.GetBinContent(ix,iy);
+ Float_t aodErr = aodHist.GetBinError(ix,iy);
#if 0 // This is done once in the FindEtaBins function
// Set underflow bin
Float_t secCor = 0;
- if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
+ if(hSecMap) secCor = hSecMap->GetBinContent(ix,iy);
if (secCor > 0.5) etabinSeen = kTRUE;
#endif
if (aodValue < 0.000001) {
- aodHist.SetBinContent(nx,ny, 0);
- aodHist.SetBinError(nx,ny, 0);
+ aodHist.SetBinContent(ix,iy, 0);
+ aodHist.SetBinError(ix,iy, 0);
continue;
}
if (!fUseAcceptance) continue;
Float_t aodNew = aodValue / accCor ;
Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
TMath::Power(accErr/accCor,2) );
- aodHist.SetBinContent(nx,ny, aodNew);
+ aodHist.SetBinContent(ix,iy, aodNew);
//test
- aodHist.SetBinError(nx,ny,error);
- aodHist.SetBinError(nx,ny,aodErr);
+ aodHist.SetBinError(ix,iy,error);
+ aodHist.SetBinError(ix,iy,aodErr);
}
//Filling underflow bin if we eta bin is in range
- if (fiducial) aodHist.SetBinContent(nx,0, 1.);
- // if (etabinSeen) aodHist.SetBinContent(nx,0, 1.);
+ if (fiducial) {
+ aodHist.SetBinContent(ix,0, 1.);
+ aodHist.SetBinContent(ix,nY+1, 1.);
+ }
+ // if (etabinSeen) aodHist.SetBinContent(ix,0, 1.);
}
}
for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
std::cout << " " << std::setw(2) << v << " "
<< std::setw(5) << vaxis.GetBinLowEdge(v) << "-"
- << std::setw(5) << vaxis.GetBinUpEdge(v) << " | "
- << std::setw(3) << fEtaMin[v-1] << "-"
- << std::setw(3) << fEtaMax[v-1] << std::endl;
+ << std::setw(5) << vaxis.GetBinUpEdge(v) << " | ";
+ if (fEtaMin.GetSize() <= 0)
+ std::cout << " ? - ?";
+ else
+ std::cout << std::setw(3) << fEtaMin[v-1] << "-"
+ << std::setw(3) << fEtaMax[v-1];
+ std::cout << std::endl;
}
}
Form("Failed to write %s to disk (%d)", ofName.Data(),ret));
return false;
}
- output->ls();
+ // output->ls();
output->Close();
+#if 0
TString cName(obj->IsA()->GetName());
AliInfoGeneral("Manager",
Form("Wrote %s object %s to %s\n",
if (!full) {
TString dName(GetFileDir(what));
AliInfoGeneral("Manager",
- Form("%s should be copied to %s\n"
+ Form("\n %s should be copied to %s\n"
"Do for example\n\t"
"aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/"
"MoveCorrections.C\\(%d\\)\nor\n\t"
}
+#endif
return true;
}
AliCentraldNdetaTask::AliCentraldNdetaTask()
: AliBasedNdetaTask()
{
- DGUARD(fDebug,0,"Default construction of AliCentraldNdetaTask");
+ DGUARD(fDebug,3,"Default CTOR of AliCentraldNdetaTask");
}
//____________________________________________________________________
AliCentraldNdetaTask::AliCentraldNdetaTask(const char*)
: AliBasedNdetaTask("Central")
{
- DGUARD(fDebug,0,"Named construction of AliCentraldNdetaTask");
+ DGUARD(fDebug,3,"Named CTOR of AliCentraldNdetaTask");
fSymmetrice = false;
fCorrEmpty = false;
SetTitle("Central");
#include <TObjArray.h>
#include <TAxis.h>
class TH2D;
+class TH1D;
/**
* This class contains the acceptance correction due to dead channels
* @return The correction @f$ a_{r,v}@f$
*/
TH2D* GetCorrection(UShort_t d, Char_t r, UShort_t b) const;
+ /**
+ * Get the phi acceptance
+ *
+ * @param d
+ * @param r
+ * @param v
+ *
+ * @return
+ */
+ TH1D* GetPhiAcceptance(UShort_t d, Char_t r, Double_t v) const;
+ /**
+ * Get the phi acceptance
+ *
+ * @param d
+ * @param r
+ * @param v
+ *
+ * @return
+ */
+ TH1D* GetPhiAcceptance(UShort_t d, Char_t r, UShort_t b) const;
/**
* Get the vertex axis used
*
* @return vertex axis
*/
const TAxis& GetVertexAxis() const { return fVertexAxis; }
+ /**
+ * @return true if the overflow bins along eta are set to the ratio
+ * of OK strips to all strips for a given eta bin.
+ */
+ Bool_t HasOverflow() const { return fHasOverflow; }
/* @} */
/**
* @param max Maximum
*/
void SetVertexAxis(Int_t nBins, Double_t min, Double_t max);
+ /**
+ * Set that we have (or not) the overflow bin present
+ *
+ * @param present If true, assume the overflow bins are filled
+ */
+ void SetHasOverflow(Bool_t present=true) { fHasOverflow = present; }
/* @} */
/**
*
* @return Pointer to ring array, or null in case of problems
*/
- TObjArray* GetRingArray(UShort_t d, Char_t r) const;
+ TObjArray* GetRingArray(const TObjArray& m, UShort_t d, Char_t r) const;
/**
* Get the ring array corresponding to the specified ring
*
*
* @return Pointer to ring array, or newly created container
*/
- TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
+ TObjArray* GetOrMakeRingArray(TObjArray& m, UShort_t d, Char_t r) const;
+ TObject* GetObject(const TObjArray& m, UShort_t d,
+ Char_t r, UShort_t b) const;
+
+ void FillCache() const;
- TObjArray fRingArray; // Array of per-ring, per-vertex 2nd map
+ TObjArray fRingArray; // Array of per-ring, per-vertex 2d map
+ mutable TObjArray* fCache; //! Array of per-ring, per-vertex 1d factors
TAxis fVertexAxis; // The vertex axis
- ClassDef(AliFMDCorrAcceptance,1); // Acceptance correction due to dead areas
+ Bool_t fHasOverflow; // Whether we have the overflow bin set
+ ClassDef(AliFMDCorrAcceptance,2); // Acceptance correction due to dead areas
};
//____________________________________________________________________
fDebug(0)
{
// Constructor
- DGUARD(fDebug, 0, "Default CTOR of AliFMDCorrector");
+ DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
}
//____________________________________________________________________
//
// Parameters:
// title Title
- DGUARD(fDebug, 0, "Named CTOR of AliFMDCorrector: %s", title);
+ DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
fRingHistos.SetName(GetName());
fRingHistos.Add(new RingHistos(1, 'I'));
fRingHistos.Add(new RingHistos(2, 'I'));
//
// Parameters:
// o Object to copy from
- DGUARD(fDebug, 0, "Copy CTOR of AliFMDCorrector");
+ DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
TIter next(&o.fRingHistos);
TObject* obj = 0;
while ((obj = next())) fRingHistos.Add(obj);
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
+//____________________________________________________________________
+void
+AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
+ Bool_t alsoUnderOver) const
+{
+ //
+ // Implement TH1::Divide but
+ // - Assume compatible histograms
+ // - Unless third argument is true, do not divide over/under flow bins
+ //
+ if (!num || !denom) return;
+
+ Int_t first = (alsoUnderOver ? 0 : 1);
+ Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
+ Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
+
+ for (Int_t ix = first; ix <= lastX; ix++) {
+ for (Int_t iy = first; iy <= lastY; iy++) {
+ Int_t bin = num->GetBin(ix,iy);
+ Double_t c0 = num->GetBinContent(bin);
+ Double_t c1 = denom->GetBinContent(bin);
+ if (!c1) {
+ num->SetBinContent(bin,0);
+ num->SetBinError(bin, 0);
+ continue;
+ }
+ Double_t w = c0 / c1;
+ Double_t e0 = num->GetBinError(bin);
+ Double_t e1 = denom->GetBinError(bin);
+ Double_t c12 = c1*c1;
+ Double_t e2 = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
+
+ num->SetBinContent(bin, w);
+ num->SetBinError(bin, TMath::Sqrt(e2));
+ }
+ }
+}
//____________________________________________________________________
Bool_t
AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
continue;
}
// Divide by primary/total ratio
- h->Divide(bg);
+ DivideMap(h, bg, false);
}
if (fUseVertexBias) {
TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
continue;
}
// Divide by the event selection efficiency
- h->Divide(ef);
+ DivideMap(h, ef, false);
}
if (fUseAcceptance) {
TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
"vertex bin %d", d, r, uvb));
continue;
}
- // Divide by the acceptance correction
- h->Divide(ac);
+ // Fill overflow bin with ones
+ for (Int_t i = 1; i <= h->GetNbinsX(); i++)
+ h->SetBinContent(i, h->GetNbinsY()+1, 1);
+
+ // Divide by the acceptance correction -
+ DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
}
if (fUseMergingEfficiency) {
for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
Float_t c = sf->GetBinContent(ieta);
Float_t ec = sf->GetBinError(ieta);
-
- if (c == 0) continue;
-
+
for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
+ if (c == 0) {
+ h->SetBinContent(ieta,iphi,0);
+ h->SetBinError(ieta,iphi,0);
+ continue;
+ }
+
Double_t m = h->GetBinContent(ieta, iphi) / c;
Double_t em = h->GetBinError(ieta, iphi);
}
}
}
- //HHD
- /*
- TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
- TH2D hRing("hring","hring",bg->GetNbinsX(),
- bg->GetXaxis()->GetXmin(),
- bg->GetXaxis()->GetXmax(),
- bg->GetNbinsY(),
- bg->GetYaxis()->GetXmin(),
- bg->GetYaxis()->GetXmax());
-
- Int_t edgebin[4] = {0,0,0,0};
- for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
- for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
- Float_t bgcor = bg->GetBinContent(ii,jj);
- if(bgcor<0.1) continue;
- if(edgebin[0] == 0) edgebin[0] = ii;
- if(edgebin[0] == ii) continue;
- if(edgebin[0] > 0 && edgebin[1] == 0) edgebin[1] = ii;
- if(edgebin[0]>0 && edgebin[1]>0) break;
- }
- }
- for(Int_t ii = bg->GetNbinsX(); ii >= 1; ii--) {
- for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
- Float_t bgcor = bg->GetBinContent(ii,jj);
- if(bgcor<0.1) continue;
- if(edgebin[2] == 0) edgebin[2] = ii;
- if(edgebin[2] == ii) continue;
- if(edgebin[2] > 0 && edgebin[3] == 0) edgebin[3] = ii;
- if(edgebin[2]>0 && edgebin[3]>0) break;
- }
- }
- for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
- for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
- Float_t data = h->GetBinContent(ii,jj);
- if(data <0.000001) continue;
- if(edgebin[0] == ii || edgebin[1] == ii || edgebin[2] == ii || edgebin[3] == ii) continue;
- hRing.SetBinContent(ii,jj,data);
- hRing.SetBinError(ii,jj,h->GetBinError(ii,jj));
- }
- }
-
- //std::cout<<edgebin[0]<<" "<<edgebin[1]<<" "<<edgebin[2]<<" "<<edgebin[3]<<std::endl;
- */
rh->fDensity->Add(h);
}
#include <TList.h>
#include "AliForwardUtil.h"
class TH2D;
-
+class TH2;
/**
* @defgroup pwglf_forward_algo Algorithms
*
* @return Ring histogram container
*/
RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
-
+ /**
+ * Divide a map with another map. This is a reimplementation of
+ * TH1::Divide, but we assume compatible histograms, and the under-
+ * and overflow bins are only divided if the third argument is true.
+ *
+ * @param num Numerator. On return contains the result
+ * @param denom Denominator
+ * @param alsoUnderOver If true, also divide under/overflow bins
+ */
+ void DivideMap(TH2* num, const TH2* denom, Bool_t alsoUnderOver=false) const;
TList fRingHistos; // List of histogram containers
Bool_t fUseSecondaryMap; // Whether to do correction for secondaries
Bool_t fUseVertexBias; // Whether to do correction for vertex bias
Bool_t fUseMergingEfficiency; // Whether to use the merging efficiency
Int_t fDebug; // Debug level
- ClassDef(AliFMDCorrector,2); // Correct the inclusive d2N/detadphi
+ ClassDef(AliFMDCorrector,3); // Correct the inclusive d2N/detadphi
};
#endif
#include <TProfile.h>
#include <THStack.h>
#include <TROOT.h>
+#include <TVector3.h>
#include <iostream>
#include <iomanip>
; // For Emacs
#endif
+//____________________________________________________________________
+const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator";
+
//____________________________________________________________________
AliFMDDensityCalculator::AliFMDDensityCalculator()
: TNamed(),
fPhiLumping(4),
fDebug(0),
fCuts(),
- fRecalculateEta(false)
+ fRecalculateEta(false),
+ fRecalculatePhi(false)
{
//
// Constructor
//
- DGUARD(fDebug, 0, "Default CTOR of FMD density calculator");
+ DGUARD(fDebug, 3, "Default CTOR of FMD density calculator");
}
//____________________________________________________________________
AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
- : TNamed("fmdDensityCalculator", title),
+ : TNamed(fgkFolderName, title),
fRingHistos(),
fSumOfWeights(0),
fWeightedSum(0),
fPhiLumping(4),
fDebug(0),
fCuts(),
- fRecalculateEta(false)
+ fRecalculateEta(false),
+ fRecalculatePhi(false)
{
//
// Constructor
// Parameters:
// name Name of object
//
- DGUARD(fDebug, 0, "Named CTOR of FMD density calculator: %s", title);
+ DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title);
fRingHistos.SetName(GetName());
fRingHistos.SetOwner();
fRingHistos.Add(new RingHistos(1, 'I'));
fPhiLumping(o.fPhiLumping),
fDebug(o.fDebug),
fCuts(o.fCuts),
- fRecalculateEta(o.fRecalculateEta)
+ fRecalculateEta(o.fRecalculateEta),
+ fRecalculatePhi(o.fRecalculatePhi)
{
//
// Copy constructor
// Parameters:
// o Object to copy from
//
- DGUARD(fDebug, 0, "Copy CTOR of FMD density calculator");
+ DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator");
TIter next(&o.fRingHistos);
TObject* obj = 0;
while ((obj = next())) fRingHistos.Add(obj);
fPhiLumping = o.fPhiLumping;
fCuts = o.fCuts;
fRecalculateEta = o.fRecalculateEta;
+ fRecalculatePhi = o.fRecalculatePhi;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
// Intialize this sub-algorithm
//
// Parameters:
- // etaAxis Not used
+ // etaAxis Eta axis
DGUARD(fDebug, 1, "Initialize FMD density calculator");
CacheMaxWeights(axis);
Bool_t
AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- UShort_t vtxbin,
Bool_t lowFlux,
Double_t /*cent*/,
- Double_t zvtx)
+ const TVector3& ip)
{
//
// Do the calculations
// Return:
// true on successs
DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
-
+
+ // --- Loop over detectors -----------------------------------------
for (UShort_t d=1; d<=3; d++) {
UShort_t nr = (d == 1 ? 1 : 2);
for (UShort_t q=0; q<nr; q++) {
}
// rh->fPoisson.SetObject(d,r,vtxbin,cent);
rh->fPoisson.Reset(0);
+ rh->fTotal->Reset();
+ rh->fGood->Reset();
// rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
+ // --- Loop over sectors and strips ----------------------------
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);
-
-
- if(fRecalculateEta)
- eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
-
+ Float_t mult = fmd.Multiplicity(d,r,s,t);
+ Float_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
+ Float_t eta = fmd.Eta(d,r,s,t);
+ Double_t oldPhi = phi;
+
+ // --- Re-calculate eta - needed for satelittes ------------
+ if( fRecalculateEta)
+ eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
+
+ // --- Check this strip ------------------------------------
+ rh->fTotal->Fill(eta);
if (mult == AliESDFMD::kInvalidMult || mult > 20) {
- rh->fPoisson.Fill(t , s, false);
- rh->fEvsM->Fill(mult,0);
+ // Do not count invalid stuff
+ // rh->fPoisson.Fill(t , s, false);
+ rh->fEvsN->Fill(mult,-1);
+ rh->fEvsM->Fill(mult,-1);
continue;
}
+ // --- Automatic calculation of acceptance -----------------
+ rh->fGood->Fill(eta);
+
+ // --- If we asked to re-calculate phi for (x,y) IP --------
+ if (fRecalculatePhi) {
+ oldPhi = phi;
+ phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
+ }
+
+ // --- Apply phi corner correction to eloss ----------------
if (fUsePhiAcceptance == kPhiCorrectELoss)
mult *= AcceptanceCorrection(r,t);
+ // --- Get the low multiplicity cut ------------------------
Double_t cut = 1024;
if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
+ // --- Now caluculate Nch for this strip using fits --------
Double_t n = 0;
- if (cut > 0 && mult > cut)
- n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
-
+ if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
rh->fELoss->Fill(mult);
rh->fEvsN->Fill(mult,n);
rh->fEtaVsN->Fill(eta, n);
- Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
+ // --- Calculate correction if needed ----------------------
+ Double_t c = Correction(d,r,t,eta,lowFlux);
fCorrections->Fill(c);
if (c > 0) n /= c;
rh->fEvsM->Fill(mult,n);
rh->fEtaVsM->Fill(eta, n);
rh->fCorr->Fill(eta, c);
+ // --- Accumulate Poisson statistics -----------------------
Bool_t hit = (n > 0.9 && c > 0);
- if (hit) rh->fELossUsed->Fill(mult);
+ if (hit) {
+ rh->fELossUsed->Fill(mult);
+ if (fRecalculatePhi) {
+ rh->fPhiBefore->Fill(oldPhi);
+ rh->fPhiAfter->Fill(phi);
+ }
+ }
rh->fPoisson.Fill(t,s,hit,1./c);
h->Fill(eta,phi,n);
+
+ // --- If we use ELoss fits, apply now ---------------------
if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
} // for t
} // for s
-
+
+ // --- Automatic acceptance - Calculate as an efficiency -------
+ rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
+
+ // --- Make a copy and reset as needed -------------------------
TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
if (!fUsePoisson) hclone->Reset();
if ( fUsePoisson) h->Reset();
+ // --- Store Poisson result ------------------------------------
TH2D* poisson = rh->fPoisson.Result();
for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
Double_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
Double_t eta = fmd.Eta(d,r,s,t);
if(fRecalculateEta)
- eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
+ eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
if (fUsePoisson)
h->Fill(eta,phi,poissonV);
else
}
}
+ // --- Make diagnostics - eloss vs poisson ---------------------
+ Int_t nY = h->GetNbinsY();
for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
- for (Int_t iphi=1; iphi<= h->GetNbinsY(); iphi++) {
+ // Set the overflow bin to contain the phi acceptance
+ Double_t phiAcc = rh->fGood->GetBinContent(ieta);
+ Double_t phiAccE = rh->fGood->GetBinError(ieta);
+ h->SetBinContent(ieta, nY+1, phiAcc);
+ h->SetBinError(ieta, nY+1, phiAccE);
+ Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
+ rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
+#if 0
+ if (phiAcc > 1e-12) {
+ Info("", "FMD%d%c, eta=%3d phi acceptance: %f (%f)",
+ d, r, ieta, phiAcc, h->GetBinContent(ieta, nY+1));
+ }
+#endif
+ for (Int_t iphi=1; iphi<= nY; iphi++) {
Double_t poissonV = 0; //h->GetBinContent(,s+1);
Double_t eLossV = 0;
}
rh->fELossVsPoisson->Fill(eLossV, poissonV);
+ rh->fDiffELossPoisson->Fill(poissonV < 1e-12 ? 0 :
+ (eLossV - poissonV) / poissonV);
+
}
}
delete hclone;
AliFMDDensityCalculator::NParticles(Float_t mult,
UShort_t d,
Char_t r,
- UShort_t /*s*/,
- UShort_t /*t*/,
- UShort_t /*v*/,
Float_t eta,
Bool_t lowFlux) const
{
Float_t
AliFMDDensityCalculator::Correction(UShort_t d,
Char_t r,
- UShort_t /*s*/,
UShort_t t,
- UShort_t /*v*/,
Float_t eta,
Bool_t lowFlux) const
{
// TNamed* sigma = new TNamed("sigma",
// (fIncludeSigma ? "included" : "excluded"));
-#if 0
- TNamed* maxP = new TNamed("maxParticle", Form("%d", fMaxParticles));
- TNamed* method = new TNamed("method",
- (fUsePoisson ? "Poisson" : "Energy loss"));
- TNamed* phiA = new TNamed("phiAcceptance",
- (fUsePhiAcceptance == 0 ? "disabled" :
- fUsePhiAcceptance == 1 ? "particles" :
- "energy loss"));
- TNamed* etaL = new TNamed("etaLumping", Form("%d", fEtaLumping));
- TNamed* phiL = new TNamed("phiLumping", Form("%d", fPhiLumping));
- // TParameter<double>* nxi = new TParameter<double>("nXi", fNXi);
-#else
TObject* maxP = AliForwardUtil::MakeParameter("maxParticle", fMaxParticles);
TObject* method = AliForwardUtil::MakeParameter("method", fUsePoisson);
TObject* phiA = AliForwardUtil::MakeParameter("phiAcceptance",
fUsePhiAcceptance);
TObject* etaL = AliForwardUtil::MakeParameter("etaLumping", fEtaLumping);
TObject* phiL = AliForwardUtil::MakeParameter("phiLumping", fPhiLumping);
-#endif
+ TObject* reEt = AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta);
+ TObject* rePh = AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi);
+
// d->Add(sigma);
d->Add(maxP);
d->Add(method);
d->Add(phiA);
d->Add(etaL);
d->Add(phiL);
+ d->Add(reEt);
+ d->Add(rePh);
// d->Add(nxi);
- fCuts.Output(d,0);
+ fCuts.Output(d,"lCuts");
TIter next(&fRingHistos);
RingHistos* o = 0;
<< std::boolalpha
<< 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'
+ << ind << " Recalculate eta: " << fRecalculateEta << '\n'
+ << ind << " Recalculate phi: " << fRecalculatePhi << '\n'
+ << ind << " Use phi acceptance: "
<< std::noboolalpha
<< std::flush;
+ switch (fUsePhiAcceptance) {
+ case kPhiNoCorrect: std::cout << "none\n"; break;
+ case kPhiCorrectNch: std::cout << "correct Nch\n"; break;
+ case kPhiCorrectELoss: std::cout << "correct energy loss\n"; break;
+ }
std::cout << ind << " Lower cut:" << std::endl;
fCuts.Print();
TString opt(option);
//====================================================================
AliFMDDensityCalculator::RingHistos::RingHistos()
: AliForwardUtil::RingHistos(),
+ fList(0),
fEvsN(0),
fEvsM(0),
fEtaVsN(0),
fCorr(0),
fDensity(0),
fELossVsPoisson(0),
+ fDiffELossPoisson(0),
fPoisson(),
fELoss(0),
fELossUsed(0),
- fMultCut(0)
+ fMultCut(0),
+ fTotal(0),
+ fGood(0),
+ fPhiAcc(0),
+ fPhiBefore(0),
+ fPhiAfter(0)
{
//
// Default CTOR
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
: AliForwardUtil::RingHistos(d,r),
+ fList(0),
fEvsN(0),
fEvsM(0),
fEtaVsN(0),
fCorr(0),
fDensity(0),
fELossVsPoisson(0),
+ fDiffELossPoisson(0),
fPoisson("ignored"),
fELoss(0),
fELossUsed(0),
- fMultCut(0)
+ fMultCut(0),
+ fTotal(0),
+ fGood(0),
+ fPhiAcc(0),
+ fPhiBefore(0),
+ fPhiAfter(0)
{
//
// Constructor
//
fEvsN = new TH2D("elossVsNnocorr",
"#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
- 250, -.5, 24.5, 250, -.5, 24.5);
+ 250, -.5, 24.5, 251, -1.5, 24.5);
fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
fEvsN->Sumw2();
200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
0, 2*TMath::Pi());
fDensity->SetDirectory(0);
- fDensity->Sumw2();
- fDensity->SetMarkerColor(Color());
+ fDensity->Sumw2(); fDensity->SetMarkerColor(Color());
fDensity->SetXTitle("#eta");
fDensity->SetYTitle("#phi [radians]");
fDensity->SetZTitle("Inclusive N_{ch} density");
fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
fELossVsPoisson->SetZTitle("Correlation");
+ fDiffELossPoisson = new TH1D("diffElossPoisson",
+ "(N_{ch,#DeltaE}-N_{ch,Poisson})/N_{ch,Poisson}",
+ 100, -1, 1);
+ fDiffELossPoisson->SetDirectory(0);
+ fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
+ fDiffELossPoisson->SetYTitle("Frequency");
+ fDiffELossPoisson->SetMarkerColor(Color());
+ fDiffELossPoisson->SetFillColor(Color());
+ fDiffELossPoisson->Sumw2();
+
fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
600, 0, 15);
fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
fELossUsed->SetFillStyle(3002);
fELossUsed->SetLineStyle(1);
fELossUsed->SetDirectory(0);
+
+ fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
+ (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
+ fPhiBefore->SetDirectory(0);
+ fPhiBefore->SetXTitle("#phi");
+ fPhiBefore->SetYTitle("Events");
+ fPhiBefore->SetMarkerColor(Color());
+ fPhiBefore->SetLineColor(Color());
+ fPhiBefore->SetFillColor(Color());
+ fPhiBefore->SetMarkerStyle(20);
+
+ fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
+ fPhiAfter->SetTitle("#phi distribution (after re-calc)");
+ fPhiAfter->SetDirectory(0);
}
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
- : AliForwardUtil::RingHistos(o),
+ : AliForwardUtil::RingHistos(o),
+ fList(o.fList),
fEvsN(o.fEvsN),
fEvsM(o.fEvsM),
fEtaVsN(o.fEtaVsN),
fCorr(o.fCorr),
fDensity(o.fDensity),
fELossVsPoisson(o.fELossVsPoisson),
+ fDiffELossPoisson(o.fDiffELossPoisson),
fPoisson(o.fPoisson),
fELoss(o.fELoss),
fELossUsed(o.fELossUsed),
- fMultCut(o.fMultCut)
+ fMultCut(o.fMultCut),
+ fTotal(o.fTotal),
+ fGood(o.fGood),
+ fPhiAcc(o.fPhiAcc),
+ fPhiBefore(o.fPhiBefore),
+ fPhiAfter(o.fPhiAfter)
{
//
// Copy constructor
if (&o == this) return *this;
AliForwardUtil::RingHistos::operator=(o);
- if (fEvsN) delete fEvsN;
- if (fEvsM) delete fEvsM;
- if (fEtaVsN) delete fEtaVsN;
- if (fEtaVsM) delete fEtaVsM;
- if (fCorr) delete fCorr;
- if (fDensity) delete fDensity;
- if (fELossVsPoisson) delete fELossVsPoisson;
-
- fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
- fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
- fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
- fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
- fCorr = static_cast<TProfile*>(o.fCorr->Clone());
- fDensity = static_cast<TH2D*>(o.fDensity->Clone());
- fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
- fPoisson = o.fPoisson;
- fELoss = static_cast<TH1D*>(o.fELoss->Clone());
- fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
-
+ if (fEvsN) delete fEvsN;
+ if (fEvsM) delete fEvsM;
+ if (fEtaVsN) delete fEtaVsN;
+ if (fEtaVsM) delete fEtaVsM;
+ if (fCorr) delete fCorr;
+ if (fDensity) delete fDensity;
+ if (fELossVsPoisson) delete fELossVsPoisson;
+ if (fDiffELossPoisson) delete fDiffELossPoisson;
+ if (fTotal) delete fTotal;
+ if (fGood) delete fGood;
+ if (fPhiAcc) delete fPhiAcc;
+ if (fPhiBefore) delete fPhiBefore;
+ if (fPhiAfter) delete fPhiAfter;
+
+ fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
+ fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
+ fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
+ fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
+ fCorr = static_cast<TProfile*>(o.fCorr->Clone());
+ fDensity = static_cast<TH2D*>(o.fDensity->Clone());
+ fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
+ fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
+ fPoisson = o.fPoisson;
+ fELoss = static_cast<TH1D*>(o.fELoss->Clone());
+ fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
+ fTotal = static_cast<TH1D*>(o.fTotal->Clone());
+ fGood = static_cast<TH1D*>(o.fGood->Clone());
+ fPhiAcc = static_cast<TH2D*>(o.fPhiAcc->Clone());
+ fPhiBefore = static_cast<TH1D*>(o.fPhiBefore->Clone());
+ fPhiAfter = static_cast<TH1D*>(o.fPhiAfter->Clone());
return *this;
}
//____________________________________________________________________
//____________________________________________________________________
void
-AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
+AliFMDDensityCalculator::RingHistos::Init(const TAxis& eAxis)
{
+ // Initialize
+ // This is called on first event
fPoisson.Init(-1,-1);
+ fTotal = new TH1D("total", "Total # of strips per #eta",
+ eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
+ fTotal->SetDirectory(0);
+ fTotal->SetXTitle("#eta");
+ fTotal->SetYTitle("# of strips");
+ fGood = static_cast<TH1D*>(fTotal->Clone("good"));
+ fGood->SetTitle("# of good strips per #eta");
+ fGood->SetDirectory(0);
+
+ fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}",
+ eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
+ 10, -10, 10);
+ fPhiAcc->SetXTitle("#eta");
+ fPhiAcc->SetYTitle("v_{z} [cm]");
+ fPhiAcc->SetZTitle("#phi acceptance");
+ fPhiAcc->SetDirectory(0);
+
+ if (fList) fList->Add(fPhiAcc);
}
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::Output(TList* dir)
{
//
- // Make output
+ // Make output. This is called as part of SlaveBegin
//
// Parameters:
// dir Where to put it
d->Add(fCorr);
d->Add(fDensity);
d->Add(fELossVsPoisson);
+ d->Add(fDiffELossPoisson);
fPoisson.Output(d);
fPoisson.GetOccupancy()->SetFillColor(Color());
fPoisson.GetMean()->SetFillColor(Color());
fPoisson.GetOccupancy()->SetFillColor(Color());
d->Add(fELoss);
d->Add(fELossUsed);
+ d->Add(fPhiBefore);
+ d->Add(fPhiAfter);
Bool_t inner = (fRing == 'I' || fRing == 'i');
Int_t nStr = inner ? 512 : 256;
fPoisson.Define(x, y);
d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
+ fList = d;
}
//____________________________________________________________________
#include <TNamed.h>
#include <TList.h>
#include <TArrayI.h>
+#include <TVector3.h>
#include "AliForwardUtil.h"
#include "AliFMDMultCuts.h"
#include "AliPoissonCalculator.h"
/** Correct the energy loss */
kPhiCorrectELoss
};
+ /**
+ * Folder name
+ */
+ static const char* fgkFolderName;
/**
* Constructor
*/
*
* @param fmd AliESDFMD object (possibly) corrected for sharing
* @param hists Histogram cache
- * @param vtxBin Vertex bin
* @param lowFlux Low flux flag.
* @param cent Centrality
* @param vz Vertex Z position
*/
virtual Bool_t Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- UShort_t vtxBin,
Bool_t lowFlux,
Double_t cent=-1,
- Double_t vz=0);
+ const TVector3& ip=TVector3(1024,1024,0));
/**
* Scale the histograms to the total number of events
*
*
*/
void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; }
+ /**
+ * In case of a displaced vertices recalculate eta and angle correction
+ *
+ * @param use recalculate or not
+ *
+ */
+ void SetRecalculatePhi(Bool_t use) { fRecalculatePhi = use; }
/**
* Set whether to use the phi acceptance correction.
*
* @param mult Signal
* @param d Detector
* @param r Ring
- * @param s Sector
- * @param t Strip (not used)
- * @param v Vertex bin
* @param eta Pseudo-rapidity
* @param lowFlux Low-flux flag
*
* @return The number of particles
*/
- virtual Float_t NParticles(Float_t mult,
- UShort_t d, Char_t r, UShort_t s, UShort_t t,
- UShort_t v, Float_t eta, Bool_t lowFlux) const;
+ virtual Float_t NParticles(Float_t mult,
+ UShort_t d,
+ Char_t r,
+ Float_t eta,
+ Bool_t lowFlux) const;
/**
* Get the inverse correction factor. This consist of
*
*
* @return
*/
- virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
- UShort_t v, Float_t eta, Bool_t lowFlux) const;
+ virtual Float_t Correction(UShort_t d, Char_t r, UShort_t t,
+ Float_t eta, Bool_t lowFlux) const;
/**
* Get the acceptance correction for strip @a t in an ring of type @a r
*
* @param nEvents Number of events
*/
void ScaleHistograms(TList* dir, Int_t nEvents);
+ TList* fList;
TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
TProfile* fEtaVsN; // Average uncorrected Nch vs eta
TProfile* fCorr; // Average correction vs eta
TH2D* fDensity; // Distribution inclusive Nch
TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
+ TH1D* fDiffELossPoisson;// Relative difference to Poisson
AliPoissonCalculator fPoisson; // Calculate density using Poisson method
TH1D* fELoss; // Energy loss as seen by this
TH1D* fELossUsed; // Energy loss in strips with signal
Double_t fMultCut; // If set, use this
-
- ClassDef(RingHistos,6);
+ TH1D* fTotal; // Total number of strips per eta
+ TH1D* fGood; // Number of good strips per eta
+ TH2D* fPhiAcc; // Phi acceptance vs IpZ
+ TH1D* fPhiBefore; // Phi before re-calce
+ TH1D* fPhiAfter; // Phi after re-calc
+ ClassDef(RingHistos,8);
};
/**
* Get the ring histogram container
Int_t fPhiLumping; // How to lump phi bins for Poisson
Int_t fDebug; // Debug level
AliFMDMultCuts fCuts; // Cuts
- Bool_t fRecalculateEta; // //Whether to recalculate eta and angle correction (disp vtx)
+ Bool_t fRecalculateEta; // Whether to recalc eta and angle correction (disp vtx)
+ Bool_t fRecalculatePhi; // Whether to correct for (X,Y) offset
ClassDef(AliFMDDensityCalculator,7); // Calculate Nch density
};
//
// Constructor
//
- DGUARD(fDebug,0,"Default CTOR of AliFMDEnergyFitterTask");
+ DGUARD(fDebug, 3,"Default CTOR of AliFMDEnergyFitterTask");
}
//____________________________________________________________________
// Parameters:
// name Name of task
//
- DGUARD(fDebug,0,"Named CTOR of AliFMDEnergyFitterTask: %s", name);
+ DGUARD(fDebug, 3,"Named CTOR of AliFMDEnergyFitterTask: %s", name);
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class());
fBranchNames =
// Parameters:
// o Object to copy from
//
- DGUARD(fDebug,0,"COPY CTOR of AliFMDEnergyFitterTask");
+ DGUARD(fDebug, 3,"COPY CTOR of AliFMDEnergyFitterTask");
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class());
}
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip;
Double_t cent = 0;
UShort_t nClusters = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ ivz, ip, cent, nClusters);
if (found & AliFMDEventInspector::kNoEvent) return;
if (found & AliFMDEventInspector::kNoTriggers) return;
if (found & AliFMDEventInspector::kNoSPD) return;
#include "AliVVZERO.h"
//====================================================================
+const char* AliFMDEventInspector::fgkFolderName = "fmdEventInspector";
+
+//____________________________________________________________________
AliFMDEventInspector::AliFMDEventInspector()
: TNamed(),
fHEventsTr(0),
fHTriggers(0),
fHTriggerCorr(0),
fHType(0),
- fHWords(0),
- fHCent(0),
- fHCentVsQual(0),
- fHStatus(0),
- fLowFluxCut(1000),
- fMaxVzErr(0.2),
- fList(0),
- fEnergy(0),
+ fHWords(0),
+ fHCent(0),
+ fHCentVsQual(0),
+ fHStatus(0),
+ fLowFluxCut(1000),
+ fMaxVzErr(0.2),
+ fList(0),
+ fEnergy(0),
fField(999),
- fCollisionSystem(kUnknown),
+ fCollisionSystem(kUnknown),
fDebug(0),
fCentAxis(0),
fVtxAxis(10,-10,10),
- fUseFirstPhysicsVertex(true),
- fUseV0AND(false),
- fMinPileupContrib(3),
- fMinPileupDistance(0.8),
+ fUseFirstPhysicsVertex(false),
+ fUseV0AND(false),
+ fMinPileupContrib(3),
+ fMinPileupDistance(0.8),
fUseDisplacedVertices(false),
fDisplacedVertex(),
fCollWords(),
//____________________________________________________________________
AliFMDEventInspector::AliFMDEventInspector(const char* name)
- : TNamed("fmdEventInspector", name),
+ : TNamed(fgkFolderName, name),
fHEventsTr(0),
fHEventsTrVtx(0),
fHEventsAccepted(0),
fDebug(0),
fCentAxis(0),
fVtxAxis(10,-10,10),
- fUseFirstPhysicsVertex(true),
+ fUseFirstPhysicsVertex(false),
fUseV0AND(false),
fMinPileupContrib(3),
fMinPileupDistance(0.8),
AliFMDEventInspector::FetchHistograms(const TList* d,
TH1I*& hEventsTr,
TH1I*& hEventsTrVtx,
+ TH1I*& hEventsAcc,
TH1I*& hTriggers) const
{
//
DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
hEventsTr = 0;
hEventsTrVtx = 0;
- hTriggers = 0;
+ hEventsAcc = 0;
+ hTriggers = 0;
TList* dd = dynamic_cast<TList*>(d->FindObject(GetName()));
if (!dd) return kFALSE;
hEventsTr = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
+ hEventsAcc = dynamic_cast<TH1I*>(dd->FindObject("nEventsAccepted"));
hTriggers = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
- if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
+ if (!hEventsTr ||
+ !hEventsTrVtx ||
+ !hEventsAcc ||
+ !hTriggers) return kFALSE;
return kTRUE;
}
//____________________________________________________________________
UInt_t& triggers,
Bool_t& lowFlux,
UShort_t& ivz,
- Double_t& vz,
+ TVector3& ip,
Double_t& cent,
UShort_t& nClusters)
{
if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
}
- // --- Get the vertex information ----------------------------------
- Double_t vx = 0;
- Double_t vy = 0;
- vz = 0;
-
- Bool_t vzOk = ReadVertex(*event, vz,vx,vy);
-
- fHEventsTr->Fill(vz);
+ // --- Get the interaction point -----------------------------------
+ Bool_t vzOk = ReadVertex(*event, ip);
+ fHEventsTr->Fill(ip.Z());
if (!vzOk) {
if (fDebug > 3) {
AliWarning("Failed to read vertex from ESD"); }
fHStatus->Fill(6);
return kNoVertex;
}
- fHEventsTrVtx->Fill(vz);
+ fHEventsTrVtx->Fill(ip.Z());
// --- Get the vertex bin ------------------------------------------
- ivz = fVtxAxis.FindBin(vz);
+ ivz = fVtxAxis.FindBin(ip.Z());
if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
if (fDebug > 3) {
AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
- vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
+ ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
}
ivz = 0;
fHStatus->Fill(7);
return kBadVertex;
}
- fHEventsAccepted->Fill(vz);
- fHEventsAcceptedXY->Fill(vx,vy);
+ fHEventsAccepted->Fill(ip.Z());
+ fHEventsAcceptedXY->Fill(ip.X(),ip.Y());
// --- Check the FMD ESD data --------------------------------------
if (!event->GetFMDData()) {
//____________________________________________________________________
Bool_t
-AliFMDEventInspector::ReadVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy)
+AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, TVector3& ip)
{
//
// Read the vertex information from the ESD event
// @c true on success, @c false otherwise
//
DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
- vz = 0;
- vx = 1024;
- vy = 1024;
+ ip.SetXYZ(1024, 1024, 0);
if(fUseDisplacedVertices) {
Double_t zvtx = fDisplacedVertex.GetVertexZ();
if(TMath::Abs(zvtx) < 999) {
- vz = zvtx;
+ ip.SetZ(zvtx);
return true;
}
return false;
}
- if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, vz, vx, vy);
+ if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, ip);
- return CheckVertex(esd, vz, vx, vy);
+ return CheckVertex(esd, ip);
}
//____________________________________________________________________
Bool_t
AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy) const
+ TVector3& ip) const
{
// This is the code used by the 1st physics people
const AliESDVertex* vertex = esd.GetPrimaryVertex();
return false;
}
}
- vz = vertex->GetZ();
+ ip.SetZ(vertex->GetZ());
if(!vertex->IsFromVertexerZ()) {
- vx = vertex->GetX();
- vy = vertex->GetY();
+ ip.SetX(vertex->GetX());
+ ip.SetY(vertex->GetY());
}
return true;
}
//____________________________________________________________________
Bool_t
AliFMDEventInspector::CheckVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy) const
+ TVector3& ip) const
{
// Use standard SPD vertex (perhaps preferable for Pb+Pb)
// Get the vertex
if(vertex->GetNContributors() <= 0) {
DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
vertex->GetNContributors());
- vz = 0;
+ ip.SetZ(0);
return false;
}
// Check that the uncertainty isn't too large
}
// Get the z coordiante
- vz = vertex->GetZ();
+ ip.SetZ(vertex->GetZ());
const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
+
if(!vertexXY->IsFromVertexerZ()) {
- vx = vertexXY->GetX();
- vy = vertexXY->GetY();
+ ip.SetX(vertexXY->GetX());
+ ip.SetY(vertexXY->GetY());
}
return true;
}
class TH2F;
class TH2I;
class TAxis;
+class TVector3;
// class TList;
/**
kPP,
kPbPb
};
+ /**
+ * Folder name
+ */
+ static const char* fgkFolderName;
/**
* Constructor
*/
* event (according to the setting of fLowFluxCut)
* @param ivz On return, the found vertex bin (1-based). A zero
* means outside of the defined vertex range
- * @param vz On return, the z position of the interaction
+ * @param vx On return, the x position of the interaction point
+ * @param vy On return, the y position of the interaction point
+ * @param vz On return, the z position of the interaction point
* @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$
UInt_t& triggers,
Bool_t& lowFlux,
UShort_t& ivz,
- Double_t& vz,
+ TVector3& ip,
Double_t& cent,
UShort_t& nClusters);
/**
* @param d Input
* @param hEventsTr On return, pointer to histogram, or null
* @param hEventsTrVtx On return, pointer to histogram, or null
+ * @param hEventsAcc On return, pointer to histogram, or null
* @param hTriggers On return, pointer to histogram, or null
*
* @return true on success, false otherwise
Bool_t FetchHistograms(const TList* d,
TH1I*& hEventsTr,
TH1I*& hEventsTrVtx,
+ TH1I*& hEventsAcc,
TH1I*& hTriggers) const;
/**
* Read the collision system, collision energy, and L3 field setting
* Read the vertex information from the ESD event
*
* @param esd ESD event
- * @param vz On return, the vertex Z position
- * @param vx On return, the vertex X position
- * @param vy On return, the vertex Y position
+ * @param ip On return, the coordinates of the IP
*
* @return @c true on success, @c false otherwise
*/
- Bool_t ReadVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy);
+ Bool_t ReadVertex(const AliESDEvent& esd, TVector3& ip);
/**
* Check the vertex using the method used in PWG-UD. That is
*
* - Check that the dispersion and resolution is OK
*
* @param esd Data
- * @param vz On return, the Z coordinate of the IP
- * @param vx On return, possibly the X coordinate of the IP
- * @param vy On return, possibly the Y coordinate of the IP
+ * @param ip On return, the coordinates of the IP
*
* @return true if the vertex was found and met the requirements
*/
virtual Bool_t CheckPWGUDVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy) const;
+ TVector3& ip) const;
/**
* Check the vertex. That is
*
* - Check that the reslution is OK
*
* @param esd Data
- * @param vz On return, the Z coordinate of the IP
- * @param vx On return, possibly the X coordinate of the IP
- * @param vy On return, possibly the Y coordinate of the IP
+ * @param ip On return, the coordinates of the IP
*
* @return true if the vertex was found and met the requirements
*/
- virtual Bool_t CheckVertex(const AliESDEvent& esd,
- Double_t& vz,
- Double_t& vx,
- Double_t& vy) const;
+ virtual Bool_t CheckVertex(const AliESDEvent& esd, TVector3& ip) const;
/**
* Read centrality from event
*
fSkipFMDRings(0),
fBgAndHitMaps(false)
{
- DGUARD(fDebug, 0, "Default CTOR of AliFMDHistCollector");
+ DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector");
}
//____________________________________________________________________
fSkipFMDRings(0),
fBgAndHitMaps(false)
{
- DGUARD(fDebug, 0, "Named CTOR of AliFMDHistCollector: %s", title);
+ DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title);
}
//____________________________________________________________________
AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
fSkipFMDRings(o.fSkipFMDRings),
fBgAndHitMaps(o.fBgAndHitMaps)
{
- DGUARD(fDebug, 0, "Copy CTOR of AliFMDHistCollector");
+ DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
}
//____________________________________________________________________
re = oe;
}
break;
+ case kSum:
+ rc = c + oc;
+ re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature
+ break;
default:
AliError("No method for defining content of overlapping bins defined");
return;
GetFirstAndLast(d, r, vtxbin, first, last);
// Zero outside valid range
- for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) {
+ Int_t nY = t->GetNbinsY();
+ for (Int_t iPhi = 0; iPhi <= nY+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 = first; iEta <= last; iEta++)
+ for (Int_t iEta = first; iEta <= last; iEta++) {
+ // Double_t phiAcc = t->GetBinContent(iEta, nY+1);
+ // if (phiAcc > 1e-12 && phiAcc < 1)
+ // Info("", "FMD%d%c %3d phi acceptance: %f",d,r,iEta,phiAcc);
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);
+ nY = t->GetNbinsY();
// Now update profile output
for (Int_t iEta = first; iEta <= last; iEta++) {
// Get the possibly overlapping histogram
Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
- // Fill eta acceptance for this event into the phi underlow bin
+ // Get factor
+ Float_t fac = (fMergeMethod != kSum && overlap >= 0 ? .5 : 1);
+
+ // Fill eta acceptance for this event into the phi underflow bin
Float_t ooc = out.GetBinContent(iEta,0);
- Float_t noc = overlap >= 0? 0.5 : 1;
- out.SetBinContent(iEta, 0, ooc + noc);
+ out.SetBinContent(iEta, 0, ooc + fac);
+
+ // Fill phi acceptance for this event into the phi overflow bin
+ Float_t oop = out.GetBinContent(iEta,nY+1);
+ Float_t nop = t->GetBinContent(iEta,nY+1);
+#if 0
+ Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f",
+ iEta, fac, oop, nop, fac*(oop+nop));
+#endif
+ out.SetBinContent(iEta, nY+1, fac * nop + oop);
// Should we loop over h or t Y bins - I think it's t
- for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) {
+ for (Int_t iPhi = 1; iPhi <= nY; iPhi++) {
Double_t c = t->GetBinContent(iEta,iPhi);
Double_t e = t->GetBinError(iEta,iPhi);
Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
ind[gROOT->GetDirLevel()] = '\0';
std::cout << ind << ClassName() << ": " << GetName() << '\n'
<< ind << " # of cut bins: " << fNCutBins << '\n'
- << ind << " Correction cut: " << fCorrectionCut << '\n'
+ << ind << " Fiducal method: "
+ << (fFiducialMethod == kByCut ? "cut" : "distance") << "\n"
+ << ind << " Fiducial cut: " << fCorrectionCut << "\n"
<< ind << " Merge method: ";
switch (fMergeMethod) {
case kStraightMean: std::cout << "straight mean\n"; break;
case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
case kWeightedMean: std::cout << "weighted mean\n"; break;
case kLeastError: std::cout << "least error\n"; break;
+ case kSum: std::cout << "straight sum\n"; break;
}
std::cout << ind << " Bin ranges:\n" << ind << " rings ";
* c_2 & \mbox{otherwise}\end{array}\right.
* @f]
*/
- kLeastError
+ kLeastError,
+ /**
+ * Just sum the signals
+ */
+ kSum
};
/**
* How to obtain the fiducial cuts
Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
Float_t eta = fmd.Eta(d,r,s,t);
- Float_t corr = Correction(d, r, s, t, 0, eta, false);
+ Float_t corr = Correction(d, r, t, eta, false);
Float_t sig = (corr <= 0 ? 0 : mult / corr);
h->Fill(eta,phi,sig);
}
* @return Reference to this object
*/
State& operator=(const State& o);
- } fState; // State
+ } fState; //! State
UShort_t fMaxConsequtiveStrips; // Max 'cluster' size
TH2D* fNcr; // Number of clusters per track
AliESDFMD* fOutput; //! Output ESD object
- ClassDef(AliFMDMCTrackDensity,3); // Calculate track-ref density
+ ClassDef(AliFMDMCTrackDensity,4); // Calculate track-ref density
};
#endif
//
//
#include "AliFMDSharingFilter.h"
+#include "AliFMDStripIndex.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
fHCuts(),
fUseSimpleMerging(false),
fThreeStripSharing(true),
- fRecalculateEta(false)
+ fRecalculateEta(false),
+ fExtraDead(0)
{
//
// Default Constructor - do not use
fHCuts(),
fUseSimpleMerging(false),
fThreeStripSharing(true),
- fRecalculateEta(false)
+ fRecalculateEta(false),
+ fExtraDead(51200)
{
//
// Constructor
fHCuts.SetNXi(1);
fHCuts.SetIncludeSigma(1);
fLCuts.SetMultCuts(.15);
+
+ fExtraDead.Reset(-1);
}
//____________________________________________________________________
fHCuts(o.fHCuts),
fUseSimpleMerging(o.fUseSimpleMerging),
fThreeStripSharing(o.fThreeStripSharing),
- fRecalculateEta(o.fRecalculateEta)
+ fRecalculateEta(o.fRecalculateEta),
+ fExtraDead(o.fExtraDead)
{
//
// Copy constructor
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
+//____________________________________________________________________
+void
+AliFMDSharingFilter::AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t)
+{
+ if (d < 1 || d > 3) {
+ Warning("AddDead", "Invalid detector FMD%d", d);
+ return;
+ }
+ Bool_t inner = (r == 'I' || r == 'i');
+ if (d == 1 && !inner) {
+ Warning("AddDead", "Invalid ring FMD%d%c", d, r);
+ return;
+ }
+ if ((inner && s >= 20) || (!inner && s >= 40)) {
+ Warning("AddDead", "Invalid sector FMD%d%c[%02d]", d, r, s);
+ return;
+ }
+ if ((inner && t >= 512) || (!inner && t >= 256)) {
+ Warning("AddDead", "Invalid strip FMD%d%c[%02d,%03d]", d, r, s, t);
+ return;
+ }
+
+ Int_t id = AliFMDStripIndex::Pack(d, r, s, t);
+ Int_t i = 0;
+ for (i = 0; i < fExtraDead.GetSize(); i++) {
+ Int_t j = fExtraDead.At(i);
+ if (j == id) return; // Already there
+ if (j < 0) break; // Free slot
+ }
+ if (i >= fExtraDead.GetSize()) {
+ Warning("AddDead", "No free slot to add FMD%d%c[%02d,%03d] at",
+ d, r, s, t);
+ return;
+ }
+ fExtraDead[i] = id;
+}
+//____________________________________________________________________
+void
+AliFMDSharingFilter::AddDeadRegion(UShort_t d, Char_t r,
+ UShort_t s1, UShort_t s2,
+ UShort_t t1, UShort_t t2)
+{
+ // Add a dead region spanning from FMD<d><r>[<s1>,<t1>] to
+ // FMD<d><r>[<s2>,<t2>] (both inclusive)
+ for (Int_t s = s1; s <= s2; s++)
+ for (Int_t t = t1; t <= t2; t++)
+ AddDead(d, r, s, t);
+}
+//____________________________________________________________________
+Bool_t
+AliFMDSharingFilter::IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const
+{
+ Int_t id = AliFMDStripIndex::Pack(d, r, s, t);
+ for (Int_t i = 0; i < fExtraDead.GetSize(); i++) {
+ Int_t j = fExtraDead.At(i);
+ if (j == id) {
+ //Info("IsDead", "FMD%d%c[%02d,%03d] marked as dead here", d, r, s, t);
+ return true;
+ }
+ if (j < 0) break; // High water mark
+ }
+ return false;
+}
//____________________________________________________________________
void
AliFMDSharingFilter::Init(const TAxis& axis)
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 simple merging
Bool_t used = kFALSE;
Double_t eTotal = -1;
}
// Keep dead-channel information.
- if(mult == AliESDFMD::kInvalidMult)
+ if(mult == AliESDFMD::kInvalidMult || IsDead(d,r,s,t))
output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
// If no signal or dead strip, go on.
}
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);
- 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);
+ 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
+ status[t] = thisStatus;
// If we're processing on non-angle corrected data, we
// should do the angle correction here
} // End of non-simple
d->Add(AliForwardUtil::MakeParameter("lowSignal",
fZeroSharedHitsBelowThreshold));
d->Add(AliForwardUtil::MakeParameter("simple", fUseSimpleMerging));
+
+ TObjArray* extraDead = new TObjArray;
+ extraDead->SetOwner();
+ extraDead->SetName("extraDead");
+ for (Int_t i = 0; i < fExtraDead.GetSize(); i++) {
+ if (fExtraDead.At(i) < 0) break;
+ UShort_t dd, s, t;
+ Char_t r;
+ Int_t id = fExtraDead.At(i);
+ AliFMDStripIndex::Unpack(id, dd, r, s, t);
+ extraDead->Add(AliForwardUtil::MakeParameter(Form("FMD%d%c[%02d,%03d]",
+ dd, r, s, t), id));
+ }
+ d->Add(extraDead);
fLCuts.Output(d,"lCuts");
fHCuts.Output(d,"hCuts");
* @param c Cuts object
*/
void SetHCuts(const AliFMDMultCuts& c) { fHCuts = c; }
+
+ void AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t);
+ void AddDeadRegion(UShort_t d, Char_t r, UShort_t s1, UShort_t s2,
+ UShort_t t1, UShort_t t2);
protected:
/**
* Internal data structure to keep track of the histograms
* @param nextStatus Next status
*
* @return The filtered signal in the strip
- */ Double_t MultiplicityOfStrip(Double_t thisE,
+ */
+ Double_t MultiplicityOfStrip(Double_t thisE,
Double_t prevE,
Double_t nextE,
Double_t eta,
*/
virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const;
+ virtual Bool_t IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const;
TList fRingHistos; // List of histogram containers
// Double_t fLowCut; // Low cut on sharing
Bool_t fCorrectAngles; // Whether to work on angle corrected signals
Bool_t fUseSimpleMerging; //enable simple sharing by HHD
Bool_t fThreeStripSharing; //In case of simple sharing allow 3 strips
Bool_t fRecalculateEta; //Whether to recalculate eta and angle correction (disp vtx)
- ClassDef(AliFMDSharingFilter,4); //
+ TArrayI fExtraDead; // List of extra dead channels
+ ClassDef(AliFMDSharingFilter,5); //
};
#endif
AliError(Form("Failed to write %s to disk (%d)", ofName.Data(), ret));
return false;
}
- output->ls();
+ // output->ls();
output->Close();
+#if 0
TString cName(obj->IsA()->GetName());
AliInfo(Form("Wrote %s object %s to %s",
cName.Data(), oName.Data(), ofName.Data()));
if (!full) {
TString dName(GetFileDir(what));
- AliInfo(Form("%s should be copied to %s"
- "Do for example\n\n\t"
- "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/"
- "MoveCorrections.C\\(%d\\)\nor\n\t"
- "cp %s %s/\n",
- ofName.Data(),dName.Data(),
- what, ofName.Data(),
- gSystem->ExpandPathName(dName.Data())));
+ AliInfoF("\n %s should be copied to %s"
+ "Do for example\n\n\t"
+ "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/"
+ "MoveCorrections.C\\(%d\\)\nor\n\t"
+ "cp %s %s/\n",
+ ofName.Data(),dName.Data(),
+ what, ofName.Data(),
+ gSystem->ExpandPathName(dName.Data()));
}
+#endif
return true;
}
UInt_t triggers = 0; // Trigger bits
Bool_t lowFlux = true; // Low flux flag
UShort_t iVz = 0; // Vertex bin from ESD
- Double_t vZ = 1000; // Z coordinate from ESD
+ TVector3 ip(1024,1024,1000);
Double_t cent = -1; // Centrality
UShort_t iVzMc = 0; // Vertex bin from MC
Double_t vZMc = 1000; // Z coordinate of IP vertex from MC
Double_t phiR = 100; // Reaction plane from MC
UShort_t nClusters = 0; // Number of SPD clusters
// Process the data
- UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
+ UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip,
cent, nClusters);
fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
b, cMC, nPart, nBin, phiR);
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip(1024, 1024, 0);
Double_t cent = -1;
UShort_t nClusters = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ ivz, ip, cent, nClusters);
UShort_t ivzMC = 0;
Double_t vzMC = 0;
Double_t phiR = 0;
// UInt_t foundMC =
fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
npart, nbin, phiR);
- fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
+ fEventInspector.CompareResults(ip.Z(), vzMC, cent, b, npart, nbin);
//Store all events
MarkEventForStore();
if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
if (isAccepted) {
- fAODFMD.SetIpZ(vz);
- fMCAODFMD.SetIpZ(vz);
+ fAODFMD.SetIpZ(ip.Z());
+ fMCAODFMD.SetIpZ(ip.Z());
}
if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
AliESDFMD* esdFMD = esd->GetFMDData();
// Apply the sharing filter (or hit merging or clustering if you like)
- if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
+ if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())) {
AliWarning("Sharing filter failed!");
return;
}
- if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
+ if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),fMCESDFMD,fPrimary)) {
AliWarning("MC Sharing filter failed!");
return;
}
fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
// Calculate the inclusive charged particle density
- if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
+ if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip.Z())) {
AliWarning("Density calculator failed!");
return;
}
if (GetOutputData(1)) GetOutputData(1)->Print();
return;
}
-
- // Get our histograms from the container
- TH1I* hEventsTr = 0;
- TH1I* hEventsTrVtx = 0;
- 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;
- }
-
- TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
- if (!hData) {
- AliError(Form("Couldn't get our summed histogram from output "
- "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
- list->ls();
- return;
- }
-
- // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
- TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
- TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
- dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
- dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
- dNdeta->Divide(norm);
- dNdeta->SetStats(0);
- dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
- "width");
- list->Add(dNdeta);
- list->Add(norm);
+ Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
+ MakeSimpledNdeta(list, list, nTr, nTrVtx, nAcc);
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()));
+ fSharingFilter.ScaleHistograms(list,Int_t(nTr));
+ fDensityCalculator.ScaleHistograms(list,Int_t(nTrVtx));
+ fCorrections.ScaleHistograms(list,Int_t(nTrVtx));
}
fFirstEvent(true),
fCorrManager(0)
{
- DGUARD(0,0,"Named Construction of AliForwardMultiplicityBase %s",name);
+ DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name);
// Set our persistent pointer
fCorrManager = &AliForwardCorrectionManager::Instance();
fBranchNames =
ah->SetFillAOD(kTRUE);
}
+//____________________________________________________________________
+Bool_t
+AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input,
+ TList* output,
+ Double_t& nTr,
+ Double_t& nTrVtx,
+ Double_t& nAcc)
+{
+ // Get our histograms from the container
+ TH1I* hEventsTr = 0;
+ TH1I* hEventsTrVtx = 0;
+ TH1I* hEventsAcc = 0;
+ TH1I* hTriggers = 0;
+ if (!GetEventInspector().FetchHistograms(input,
+ hEventsTr,
+ hEventsTrVtx,
+ hEventsAcc,
+ hTriggers)) {
+ AliError(Form("Didn't get histograms from event selector "
+ "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
+ hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
+ input->ls();
+ return false;
+ }
+ nTr = hEventsTr->Integral();
+ nTrVtx = hEventsTrVtx->Integral();
+ nAcc = hEventsAcc->Integral();
+ Double_t vtxEff = nTrVtx / nTr;
+ TH2D* hData = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
+ if (!hData) {
+ AliError(Form("Couldn't get our summed histogram from output "
+ "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
+ input->ls();
+ return false;
+ }
+
+ Int_t nY = hData->GetNbinsY();
+ TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, nY, "e");
+ TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1, nY, "e");
+ TH1D* norm = hData->ProjectionX("norm", 0, 0, "");
+ TH1D* phi = hData->ProjectionX("phi", nY+1, nY+1, "");
+ dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
+ dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ dNdeta->SetMarkerColor(kRed+1);
+ dNdeta->SetMarkerStyle(20);
+ dNdeta->SetDirectory(0);
+
+ dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
+ dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ dNdeta_->SetMarkerColor(kMagenta+1);
+ dNdeta_->SetMarkerStyle(21);
+ dNdeta_->SetDirectory(0);
+
+ norm->SetTitle("Normalization to #eta coverage");
+ norm->SetYTitle("#eta coverage");
+ norm->SetMarkerColor(kBlue+1);
+ norm->SetMarkerStyle(21);
+ norm->SetFillColor(kBlue+1);
+ norm->SetFillStyle(3005);
+ norm->SetDirectory(0);
+
+ phi->SetTitle("Normalization to #phi acceptance");
+ phi->SetYTitle("#phi acceptance");
+ phi->SetMarkerColor(kGreen+1);
+ phi->SetMarkerStyle(20);
+ phi->SetFillColor(kGreen+1);
+ phi->SetFillStyle(3004);
+ // phi->Scale(1. / nAcc);
+ phi->SetDirectory(0);
+
+ // dNdeta->Divide(norm);
+ dNdeta->Divide(phi);
+ dNdeta->SetStats(0);
+ dNdeta->Scale(vtxEff, "width");
+
+ dNdeta_->Divide(norm);
+ dNdeta_->SetStats(0);
+ dNdeta_->Scale(vtxEff, "width");
+
+ output->Add(dNdeta);
+ output->Add(dNdeta_);
+ output->Add(norm);
+ output->Add(phi);
+
+ return true;
+}
+
+
//____________________________________________________________________
void
AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
ptr++;
continue;
}
- TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
- copy->SetDirectory(0);
- thisList->Add(copy);
+ TH2D* sumPhi = static_cast<TH2D*>(h->Clone("sum_phi"));
+ sumPhi->SetDirectory(0);
+ thisList->Add(sumPhi);
+
+ TH2D* sumEta = static_cast<TH2D*>(h->Clone("sum_eta"));
+ sumEta->SetDirectory(0);
+ thisList->Add(sumEta);
- 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);
- }
- }
+ Int_t nY = sumEta->GetNbinsY();
+ TH1D* etaCov =static_cast<TH1D*>(h->ProjectionX("etaCov", 0, 0, ""));
+ TH1D* phiAcc =static_cast<TH1D*>(h->ProjectionX("phiAcc", nY+1, nY+1, ""));
- 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");
+ etaCov->SetTitle("Normalization to #eta coverage");
+ etaCov->SetYTitle("#eta coverage");
+ etaCov->SetMarkerColor(kBlue+1);
+ etaCov->SetFillColor(kBlue+1);
+ etaCov->SetFillStyle(3005);
+ etaCov->SetDirectory(0);
- if(norm->GetMaximum() > 0) {
- proj->Scale(1. / norm->GetMaximum(), "width");
- norm->Scale(1. / norm->GetMaximum());
+ phiAcc->SetTitle("Normalization to #phi acceptance");
+ phiAcc->SetYTitle("#phi acceptance");
+ phiAcc->SetMarkerColor(kGreen+1);
+ phiAcc->SetFillColor(kGreen+1);
+ phiAcc->SetFillStyle(3004);
+ // phiAcc->Scale(1. / nAcc);
+ phiAcc->SetDirectory(0);
+
+ // Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1);
+ for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) {
+ for (Int_t j = 1; j <= nY; j++) {
+ Double_t c = sumEta->GetBinContent(i, j);
+ Double_t e = sumEta->GetBinError(i, j);
+ Double_t a = etaCov->GetBinContent(i);
+ Double_t p = phiAcc->GetBinContent(i);
+ // Double_t t = p; // * a
+ sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a);
+ sumEta->SetBinError( i, j, a <= 0 ? 0 : e / a);
+ sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p);
+ sumPhi->SetBinError( i, j, p <= 0 ? 0 : e / p);
+ }
}
-
- res->SetMarkerStyle(style);
- norm->SetDirectory(0);
- res->SetDirectory(0);
- proj->SetDirectory(0);
- thisList->Add(norm);
- thisList->Add(res);
- thisList->Add(proj);
- dndetaRings->Add(res);
+ // etaCov->Scale(s);
+ // phiAcc->Scale(s);
+
+ TH1D* resPhi =static_cast<TH1D*>(sumPhi->ProjectionX("dndeta_phi",
+ 1,nY,"e"));
+ resPhi->SetMarkerStyle(style);
+ resPhi->SetDirectory(0);
+ resPhi->Scale(1, "width");
+
+ TH1D* resEta =static_cast<TH1D*>(sumEta->ProjectionX("dndeta_eta",
+ 1,nY,"e"));
+ resEta->SetMarkerStyle(style);
+ resEta->SetDirectory(0);
+ resEta->Scale(1, "width");
+
+ thisList->Add(resEta);
+ thisList->Add(etaCov);
+ thisList->Add(resPhi);
+ thisList->Add(phiAcc);
+ dndetaRings->Add(resPhi);
ptr++;
}
out->Add(dndetaRings);
}
-
//____________________________________________________________________
void
AliForwardMultiplicityBase::Print(Option_t* option) const
*
*/
virtual void MarkEventForStore() const;
+ /**
+ * Calculate a simple dN/deta from all accepted events
+ *
+ * @param input Input list
+ * @param output Output list
+ * @param nTr On return, number of triggers
+ * @param nTrVtx On return, number of trigger+vertex events
+ * @param nAcc On return, number of accepted events
+ *
+ * @return true on success
+ */
+ virtual Bool_t MakeSimpledNdeta(const TList* input,
+ TList* output,
+ Double_t& nTr,
+ Double_t& nTrVtx,
+ Double_t& nAcc);
/**
* Make Ring @f$ dN/d\eta @f$ histogram and a stack
*
//
// Constructor
//
- DGUARD(0,0,"Default construction of AliForwardMultiplicityTask");
+ DGUARD(fDebug, 3,"Default CTOR of AliForwardMultiplicityTask");
}
//____________________________________________________________________
// Parameters:
// name Name of task
//
- DGUARD(0,0,"named construction of AliForwardMultiplicityTask: %s", name);
+ DGUARD(fDebug, 3,"named CTOR of AliForwardMultiplicityTask: %s", name);
DefineOutput(1, TList::Class());
}
// Parameters:
// o Object to copy from
//
- DGUARD(0,0,"Copy construction of AliForwardMultiplicityTask");
+ DGUARD(fDebug, 3,"Copy CTOR of AliForwardMultiplicityTask");
DefineOutput(1, TList::Class());
}
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip;
Double_t cent = -1;
UShort_t nClusters = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ ivz, ip, cent, nClusters);
if (found & AliFMDEventInspector::kNoEvent) return;
if (found & AliFMDEventInspector::kNoTriggers) return;
if (triggers & AliAODForwardMult::kPileUp) return;
- fAODFMD.SetIpZ(vz);
+ fAODFMD.SetIpZ(ip.Z());
if (found & AliFMDEventInspector::kBadVertex) return;
// Get FMD data
AliESDFMD* esdFMD = esd->GetFMDData();
// // Apply the sharing filter (or hit merging or clustering if you like)
- if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
+ if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) {
AliWarning("Sharing filter failed!");
return;
}
// Calculate the inclusive charged particle density
- if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent,vz)) {
+ if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
// if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
AliWarning("Density calculator failed!");
return;
}
if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
- if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos))
+ if (!fEventPlaneFinder.FindEventplane(esd, fAODEP,
+ &(fAODFMD.GetHistogram()), &fHistos))
AliWarning("Eventplane finder failed!");
}
if (GetOutputData(1)) GetOutputData(1)->Print();
return;
}
-
- // Get our histograms from the container
- TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
- TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
- TH1I* hTriggers = 0;
- if (!fEventInspector.FetchHistograms(list, hEventsTr,
- hEventsTrVtx, hTriggers)) {
- AliError(Form("Didn't get histograms from event selector "
- "(hEventsTr=%p,hEventsTrVtx=%p)",
- hEventsTr, hEventsTrVtx));
- list->ls();
- return;
- }
-
- TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
- if (!hData) {
- AliError(Form("Couldn't get our summed histogram from output "
- "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
- list->ls();
- return;
- }
-
- // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
- TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
- TH1D* norm = hData->ProjectionX("norm", 0, 0, "");
- dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
- dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
- dNdeta->Divide(norm);
- dNdeta->SetStats(0);
- dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
- "width");
- list->Add(dNdeta);
- list->Add(norm);
+ Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
+ MakeSimpledNdeta(list, list, nTr, nTrVtx, nAcc);
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()));
+ fSharingFilter.ScaleHistograms(list,Int_t(nTr));
+ fDensityCalculator.ScaleHistograms(list,Int_t(nTrVtx));
+ fCorrections.ScaleHistograms(list,Int_t(nTrVtx));
}
//
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip;
Double_t cent = -1;
UShort_t nClusters = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ ivz, ip, cent, nClusters);
Bool_t ok = true;
if (found & AliFMDEventInspector::kNoEvent) ok = false;
AliFMDEventInspector::CodeString(found));
return;
}
- DMSG(fDebug,2,"Event triggers: %s", AliAODForwardMult::GetTriggerString(triggers));
+ DMSG(fDebug,2,"Event triggers: %s",
+ AliAODForwardMult::GetTriggerString(triggers));
// We we do not want to use low flux specific code, we disable it here.
if (!fEnableLowFlux) lowFlux = false;
}
// // Apply the sharing filter (or hit merging or clustering if you like)
- if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
+ if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) {
AliWarning("Sharing filter failed!");
return;
}
// Get our histograms from the container
TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
+ TH1I* hEventsAcc = 0;
TH1I* hTriggers = 0;
- if (!fEventInspector.FetchHistograms(list, hEventsTr,
- hEventsTrVtx, hTriggers)) {
+ if (!fEventInspector.FetchHistograms(list,
+ hEventsTr,
+ hEventsTrVtx,
+ hEventsAcc,
+ hTriggers)) {
AliError(Form("Didn't get histograms from event selector "
- "(hEventsTr=%p,hEventsTrVtx=%p)",
- hEventsTr, hEventsTrVtx));
+ "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p)",
+ hEventsTr, hEventsTrVtx,hEventsAcc));
return;
}
}
//_____________________________________________________________________
-Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx)
+Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
{
- //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta but support displaced vertices)
-
//Get max R of ring
- Double_t maxR = 0;
- Double_t minR = 0;
- Bool_t inner = false;
- switch (ring) {
- case 'i': case 'I': maxR = 17.2; minR = 4.5213; inner = true; break;
- case 'o': case 'O': maxR = 28.0; minR = 15.4; inner = false; break;
- default:
- return -99999;
- }
+ const Double_t iMinR = 4.5213;
+ const Double_t iMaxR = 17.2;
+ const Double_t oMinR = 15.4;
+ const Double_t oMaxR = 28.0;
+
+ Double_t minR = (ring == 'I' || ring == 'i') ? iMinR : oMinR;
+ Double_t maxR = (ring == 'I' || ring == 'i') ? iMaxR : oMaxR;
+ Double_t nStrips = (ring == 'I' || ring == 'i') ? 512 : 256;
+ Double_t rad = maxR - minR;
+ Double_t segment = rad / nStrips;
+ Double_t r = minR + segment*strip;
- Double_t rad = maxR- minR;
- Double_t nStrips = (ring == 'I' ? 512 : 256);
- Double_t segment = rad / nStrips;
- Double_t r = minR + segment*strip;
- Int_t hybrid = sec / 2;
+ return r;
+}
+
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
+ UShort_t sec, UShort_t strip,
+ Double_t zvtx)
+{
+ // Calculate eta from strip with vertex (redundant with
+ // AliESDFMD::Eta but support displaced vertices)
- Double_t z = 0;
+ //Get max R of ring
+ Double_t r = GetStripR(ring, strip);
+ Int_t hybrid = sec / 2;
+ Bool_t inner = (ring == 'I' || ring == 'i');
+ Double_t z = 0;
switch (det) {
- case 1: z = 320.266; break;
- case 2: z = (inner ? 83.666 : 74.966); break;
+ case 1: z = 320.266; break;
+ case 2: z = (inner ? 83.666 : 74.966); break;
case 3: z = (inner ? -63.066 : -74.966); break;
default: return -999999;
}
return eta;
}
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip,
+ Double_t phi,
+ Double_t xvtx, Double_t yvtx)
+{
+ // Calculate eta from strip with vertex (redundant with
+ // AliESDFMD::Eta but support displaced vertices)
+
+ // Unknown x,y -> no change
+ if (yvtx > 999 || xvtx > 999) return phi;
+
+ //Get max R of ring
+ Double_t r = GetStripR(ring, strip);
+ Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
+ Double_t pha = (TMath::Abs(yvtx) < 1e-12 ? 0 : TMath::ATan2(xvtx, yvtx));
+ Double_t cha = amp * TMath::Cos(phi+pha);
+ phi += cha;
+ if (phi < 0) phi += TMath::TwoPi();
+ if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
+ return phi;
+}
//====================================================================
Int_t AliForwardUtil::fgConvolutionSteps = 100;
* @return Short integer value of magnetic field in kG
*/
static Short_t ParseMagneticField(Float_t field);
+ /**
+ * Get the radius of a strip.
+ *
+ * @param ring Ring identifier 'I' or 'O'
+ * @param strip Strip number
+ *
+ * @return Radial distance from beam of the strip
+ */
+ static Double_t GetStripR(Char_t ring, UShort_t strip);
/**
* Get eta from strip
*
*
* @return eta
*/
- static Double_t GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx) ;
+ static Double_t GetEtaFromStrip(UShort_t det, Char_t ring,
+ UShort_t sec, UShort_t strip, Double_t zvtx);
+ /**
+ * Get the azimuthal angle of a strip
+ *
+ * @param ring Ring identifier 'I' or 'O'
+ * @param strip Strip number
+ * @param phi Straight forward strip phi
+ * @param xvtx Ip X coordinate
+ * @param yvtx Ip Y coordinate
+ *
+ * @return The phi angle correctef for (X,Y) off set.
+ */
+ static Double_t GetPhiFromStrip(Char_t ring, UShort_t strip,
+ Double_t phi, Double_t xvtx, Double_t yvtx);
/**
* Get a string representation of the magnetic field
*
fMean->SetDirectory(0);
fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
- 1000, 0, 100);
+ 101, -.5, 100.5);
fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
fOcc->SetYTitle("Events");
fOcc->SetFillColor(kBlue+1);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff