#include <TList.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
-// #include "AliFMDCorrDoubleHit.h"
+#include "AliFMDCorrSecondaryMap.h"
#include "AliFMDCorrVertexBias.h"
#include "AliFMDCorrMergingEfficiency.h"
#include "AliFMDCorrAcceptance.h"
fRingHistos(),
fUseSecondaryMap(true),
fUseVertexBias(false),
- fUseAcceptance(true),
+ fUseAcceptance(false),
fUseMergingEfficiency(false),
fDebug(0)
{
// Constructor
- DGUARD(fDebug, 0, "Default CTOR of AliFMDCorrector");
+ DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
}
//____________________________________________________________________
fRingHistos(),
fUseSecondaryMap(true),
fUseVertexBias(false),
- fUseAcceptance(true),
+ fUseAcceptance(false),
fUseMergingEfficiency(false),
fDebug(0)
{
//
// 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);
//____________________________________________________________________
void
-AliFMDCorrector::Init(const TAxis&)
+AliFMDCorrector::SetupForData(const TAxis&)
{
//
// Initialize this object
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);
}
//____________________________________________________________________
void
-AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
+AliFMDCorrector::Terminate(const TList* dir, TList* output, Int_t nEvents)
{
//
// Scale the histograms to the total number of events
TList* d = static_cast<TList*>(dir->FindObject(GetName()));
if (!d) return;
+ TList* out = new TList;
+ out->SetName(d->GetName());
+ out->SetOwner();
+
TIter next(&fRingHistos);
RingHistos* o = 0;
THStack* sums = new THStack("sums", "Sums of ring results");
while ((o = static_cast<RingHistos*>(next()))) {
- o->ScaleHistograms(d, nEvents);
+ o->Terminate(d, nEvents);
+ if (!o->fDensity) {
+ Warning("Terminate", "No density from %s", o->GetName());
+ continue;
+ }
TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
o->fDensity->GetNbinsY(),"e");
sum->Scale(1., "width");
sum->SetYTitle("#sum N_{ch,primary}");
sums->Add(sum);
}
- d->Add(sums);
-
+ out->Add(sums);
+ output->Add(out);
}
//____________________________________________________________________
void
-AliFMDCorrector::DefineOutput(TList* dir)
+AliFMDCorrector::CreateOutputObjects(TList* dir)
{
//
// Output diagnostic histograms to directory
TIter next(&fRingHistos);
RingHistos* o = 0;
while ((o = static_cast<RingHistos*>(next()))) {
- o->Output(d);
+ o->CreateOutputObjects(d);
}
}
+#define PF(N,V,...) \
+ AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
+#define PFB(N,FLAG) \
+ do { \
+ AliForwardUtil::PrintName(N); \
+ std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
+ } while(false)
+#define PFV(N,VALUE) \
+ do { \
+ AliForwardUtil::PrintName(N); \
+ std::cout << (VALUE) << std::endl; } while(false)
+
//____________________________________________________________________
void
AliFMDCorrector::Print(Option_t* /* option */) const
// Parameters:
// option Not used
//
- char ind[gROOT->GetDirLevel()+1];
- for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
- ind[gROOT->GetDirLevel()] = '\0';
- std::cout << ind << ClassName() << ": " << GetName() << "\n"
- << std::boolalpha
- << ind << " Use secondary maps: " << fUseSecondaryMap << "\n"
- << ind << " Use vertex bias: " << fUseVertexBias << "\n"
- << ind << " Use acceptance: " << fUseAcceptance << "\n"
- << ind << " Use merging efficiency: " << fUseMergingEfficiency
- << std::noboolalpha << std::endl;
+ AliForwardUtil::PrintTask(*this);
+ gROOT->IncreaseDirLevel();
+ PFB("Use secondary maps", fUseSecondaryMap );
+ PFB("Use vertex bias", fUseVertexBias );
+ PFB("Use acceptance", fUseAcceptance );
+ PFB("Use merging efficiency", fUseMergingEfficiency);
+ gROOT->DecreaseDirLevel();
}
//====================================================================
//
// Destructor
//
- if (fDensity) delete fDensity;
+ // if (fDensity) delete fDensity;
}
//____________________________________________________________________
void
-AliFMDCorrector::RingHistos::Output(TList* dir)
+AliFMDCorrector::RingHistos::CreateOutputObjects(TList* dir)
{
//
// Make output
//____________________________________________________________________
void
-AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
+AliFMDCorrector::RingHistos::Terminate(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,"primaryDensity");
+ TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
if (density) density->Scale(1./nEvents);
+ fDensity = density;
}
//____________________________________________________________________