fFMD3oMax(0),
fMaxWeights(0),
fLowCuts(0),
- fEtaLumping(5),
- fPhiLumping(5),
+ fEtaLumping(32),
+ fPhiLumping(4),
fDebug(0),
- fCuts()
+ fCuts(),
+ fUseRunningAverage(false)
{
//
// Constructor
fFMD3oMax(0),
fMaxWeights(0),
fLowCuts(0),
- fEtaLumping(5),
- fPhiLumping(5),
+ fEtaLumping(32),
+ fPhiLumping(4),
fDebug(0),
- fCuts()
+ fCuts(),
+ fUseRunningAverage(false)
{
//
// Constructor
fEtaLumping(o.fEtaLumping),
fPhiLumping(o.fPhiLumping),
fDebug(o.fDebug),
- fCuts(o.fCuts)
+ fCuts(o.fCuts),
+ fUseRunningAverage(o.fUseRunningAverage)
{
//
// Copy constructor
fEtaLumping = o.fEtaLumping;
fPhiLumping = o.fPhiLumping;
fCuts = o.fCuts;
+ fUseRunningAverage = o.fUseRunningAverage;
fRingHistos.Delete();
TIter next(&o.fRingHistos);
AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
UShort_t vtxbin,
- Bool_t lowFlux)
+ Bool_t lowFlux,
+ Double_t cent)
{
//
// Do the calculations
//
// Return:
// true on successs
- //
+
+
+
for (UShort_t d=1; d<=3; d++) {
UShort_t nr = (d == 1 ? 1 : 2);
for (UShort_t q=0; q<nr; q++) {
fRingHistos.ls();
return false;
}
+ rh->fPoisson.SetObject(d,r,vtxbin,cent);
rh->fPoisson.Reset(h);
// rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
if (mult == AliESDFMD::kInvalidMult || mult > 20) {
// rh->fEmptyStrips->Fill(eta,phi);
- rh->fPoisson.Fill(eta, phi, false);
+ rh->fPoisson.Fill(t , s, false);
rh->fEvsM->Fill(mult,0);
continue;
}
Bool_t hit = (n > 0.9 && c > 0);
if (hit) rh->fELossUsed->Fill(mult);
- rh->fPoisson.Fill(eta,phi,hit,1./c);
- // if (n > 0.9 && c > 0) rh->fELossUsed->Fill(mult);
- // if (n > 0.9 && c > 0) rh->fBasicHits->Fill(eta,phi, 1./c);
- // else rh->fEmptyStrips->Fill(eta,phi);
-
+ rh->fPoisson.Fill(t,s,hit,1./c);
h->Fill(eta,phi,n);
if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
} // for t
} // for s
-
-
- // --- Loop over poisson histograms ----------------------------
+
+ TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
+ if (!fUsePoisson) hclone->Reset();
+ if ( fUsePoisson) h->Reset();
+
TH2D* poisson = rh->fPoisson.Result();
- for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
- for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
- Double_t poissonV = poisson->GetBinContent(ieta, iphi);
- Double_t poissonE = poisson->GetBinError(ieta, iphi);
- Double_t eLossV = h->GetBinContent(ieta, iphi);
- Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
- Double_t phi = h->GetYaxis()->GetBinCenter(iphi);
-
- rh->fELossVsPoisson->Fill(eLossV, poissonV);
- if (!fUsePoisson) continue;
-
- h->SetBinContent(ieta,iphi,poissonV);
- h->SetBinError(ieta,iphi,poissonE);
+ for (Int_t t=0; t <= poisson->GetNbinsX(); t++) {
+ for (Int_t s=0; s<= poisson->GetNbinsY(); s++) {
+
+ Double_t poissonV = poisson->GetBinContent(t+1,s+1);
+ Double_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+ Double_t eta = fmd.Eta(d,r,s,t);
+ if (fUsePoisson)
+ h->Fill(eta,phi,poissonV);
+ else
+ hclone->Fill(eta,phi,poissonV);
rh->fDensity->Fill(eta, phi, poissonV);
-#if 0
- Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
- Double_t phi = h->GetYaxis()->GetBinCenter(iphi);
- Int_t jeta = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
- Int_t jphi = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
- Double_t empty = rh->fEmptyStrips->GetBinContent(jeta, jphi);
- Double_t total = rh->fTotalStrips->GetBinContent(jeta, jphi);
- Double_t hits = rh->fBasicHits->GetBinContent(ieta,iphi);
-
- // Mean in region of interest
- Double_t poissonM = (total <= 0 || empty <= 0 ? 0 :
- -TMath::Log(empty / total));
- //Full occupancy should give high correction not zero
- if(empty < 0.001 && total > 0 ) poissonM = -TMath::Log( 1. / total);
+ }
+ }
+
+ for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
+ for (Int_t iphi=1; iphi<= h->GetNbinsY(); iphi++) {
- // Note, that given filled=total-empty, and
- //
- // m = -log(empty/total)
- // = -log(1 - filled/total)
- //
- // v = m / (1 - exp(-m))
- // = -total/filled * (log(total-filled)-log(total))
- // = -total / (total-empty) * log(empty/total)
- // = total (log(total)-log(empty)) / (total-empty)
- //
- Double_t poissonV = hits;
- if(poissonM > 0)
- // Correct for counting statistics and weight by counts
- poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
- Double_t poissonE = TMath::Sqrt(hits);
- if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
+ Double_t poissonV = 0; //h->GetBinContent(,s+1);
+ Double_t eLossV = 0;
+ if(fUsePoisson) {
+ poissonV = h->GetBinContent(ieta,iphi);
+ eLossV = hclone->GetBinContent(ieta,iphi);
+ }
+ else {
+ poissonV = hclone->GetBinContent(ieta,iphi);
+ eLossV = h->GetBinContent(ieta,iphi);
+ }
rh->fELossVsPoisson->Fill(eLossV, poissonV);
- rh->fEmptyVsTotal->Fill(total, empty);
- if (fUsePoisson) {
- h->SetBinContent(ieta,iphi,poissonV);
- h->SetBinError(ieta,iphi,poissonE);
- rh->fDensity->Fill(eta, phi, poissonV);
- }
-#endif
}
}
-
+ delete hclone;
+
} // for q
} // for d
TIter next(&fRingHistos);
RingHistos* o = 0;
while ((o = static_cast<RingHistos*>(next()))) {
- o->fPoisson.SetEtaLumping(fEtaLumping);
- o->fPoisson.SetPhiLumping(fPhiLumping);
- o->fPoisson.Init();
+ // o->fPoisson.SetEtaLumping(fEtaLumping);
+ o->fPoisson.SetUseAverageOverEvents(fUseRunningAverage);
+ o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
o->fPoisson.GetOccupancy()->SetFillColor(o->Color());
o->fPoisson.GetMean()->SetFillColor(o->Color());
// o->fPoisson.GetOccupancy()->SetFillColor(o->Color());
void
AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
{
- fPoisson.Init();
+ fPoisson.Init(fDet,fRing,-1,-1);
}
//____________________________________________________________________
*/
virtual Bool_t Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- UShort_t vtxBin, Bool_t lowFlux);
+ UShort_t vtxBin, Bool_t lowFlux, Double_t cent=-1);
/**
* Scale the histograms to the total number of events
*
* @param dbg Debug level
*/
void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /**
+ * Set to use the running average in Poisson
+ *
+ * @param use use or not
+ */
+ void SetUseRunningAverage(Bool_t use) { fUseRunningAverage = use; }
/**
* Maximum particle weight to use
*
Int_t fPhiLumping; // How to lump phi bins for Poisson
Int_t fDebug; // Debug level
AliFMDMultCuts fCuts; // Cuts
+ Bool_t fUseRunningAverage; //Use running average for Poisson
ClassDef(AliFMDDensityCalculator,6); // Calculate Nch density
};
UInt_t triggers = 0;
UShort_t ivz = 0;
Double_t vz = 0;
- Double_t cent = 0;
+ Double_t cent = -1;
UShort_t nClusters = 0;
UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
ivz, vz, cent, nClusters);
fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
// Calculate the inclusive charged particle density
- if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
+ if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) {
AliWarning("Density calculator failed!");
return;
}
}
// Calculate the inclusive charged particle density
- if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
+ if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) {
// if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
AliWarning("Density calculator failed!");
return;
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator()
: TNamed(),
- fEtaLumping(5),
- fPhiLumping(5),
+ fEtaLumping(32),
+ fPhiLumping(4),
fTotal(0),
fEmpty(0),
fBasic(0),
fEmptyVsTotal(0),
fMean(0),
fOcc(0),
- fCorr(0)
+ fCorr(0),
+ fTotalList(),
+ fEmptyList(),
+ fRunningAverage(false)
{
//
// CTOR
}
//____________________________________________________________________
-AliPoissonCalculator::AliPoissonCalculator(const char*)
+AliPoissonCalculator::AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/)
: TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
- fEtaLumping(5),
- fPhiLumping(5),
+ fEtaLumping(32),
+ fPhiLumping(4),
fTotal(0),
fEmpty(0),
fBasic(0),
fEmptyVsTotal(0),
fMean(0),
fOcc(0),
- fCorr(0)
+ fCorr(0),
+ fTotalList(),
+ fEmptyList(),
+ fRunningAverage(false)
{
//
// CTOR
- //
+ //
+ fEmptyList.SetOwner();
+ fTotalList.SetOwner();
+
+
}
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
fEmptyVsTotal(0),
fMean(0),
fOcc(0),
- fCorr(0)
+ fCorr(0),
+ fTotalList(),
+ fEmptyList(),
+ fRunningAverage(o.fRunningAverage)
{
Init();
Reset(o.fBasic);
TNamed::operator=(o);
fEtaLumping = o.fEtaLumping;
fPhiLumping = o.fPhiLumping;
+ fRunningAverage = o.fRunningAverage;
CleanUp();
- Init(-1,-1);
+ Init();
Reset(o.fBasic);
return *this;
}
//____________________________________________________________________
void
-AliPoissonCalculator::Init(Int_t etaLumping, Int_t phiLumping)
+AliPoissonCalculator::Init(UShort_t d, Char_t r, Int_t etaLumping, Int_t phiLumping)
{
//
// Initialize
//
if (etaLumping > 0) SetEtaLumping(etaLumping);
if (phiLumping > 0) SetPhiLumping(phiLumping);
+ if(d > 0) {
+
+ Int_t nEtaF = (r == 'I' ? 512 : 256);
+ Int_t nEta = nEtaF / fEtaLumping;
+ Double_t etaMin = -0.5;
+ Double_t etaMax = nEtaF-0.5 ;
+ Int_t nPhiF = (r == 'I' ? 20 : 40);
+ Int_t nPhi = nPhiF / fPhiLumping;
+ Double_t phiMin = -0.5;
+ Double_t phiMax = nPhiF - 0.5;
+
+ fBasic = new TH2D("basic", "Basic number of hits",
+ nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
+ fBasic->SetDirectory(0);
+ fBasic->SetXTitle("#eta");
+ fBasic->SetYTitle("#varphi [radians]");
+ fBasic->Sumw2();
+
+ for(Int_t v = 1 ; v < 11; v++) { //CHC bins
+ for(Int_t centbin = 0 ; centbin < 13; centbin++) {
+
+ TH2D* hTotal = new TH2D(Form("totalFMD%d%c_vertex%d_cent%d",d,r,v,centbin),"Total number of bins/region",
+ nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
+ TH2D* hEmpty = new TH2D(Form("emptyFMD%d%c_vertex%d_cent%d",d,r,v,centbin), "Empty number of bins/region",
+ nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
+ hEmpty->Sumw2();
+ hTotal->Sumw2();
+ fEmptyList.Add(hEmpty);
+ fTotalList.Add(hTotal);
+
+ }
+ }
+ }
+ //Create diagnostics if void
if (fEmptyVsTotal) return;
Int_t n = fEtaLumping * fPhiLumping + 1;
fCorr->SetZTitle("Events");
fCorr->SetOption("colz");
fCorr->SetDirectory(0);
+
+
+}
+//____________________________________________________________________
+void AliPoissonCalculator::SetObject(UShort_t d, Char_t r, UShort_t v, Double_t cent) {
+
+ Int_t centbin = 0;
+ if(cent > 0) {
+ if(cent > 0 && cent <5) centbin = 1;
+ if(cent > 5 && cent <10) centbin = 2;
+ else if (cent>10) centbin = (Int_t)(cent/10.) + 2;
+ }
+
+ fTotal = static_cast<TH2D*>(fTotalList.FindObject(Form("totalFMD%d%c_vertex%d_cent%d",d,r,v,centbin)));
+ fEmpty = static_cast<TH2D*>(fEmptyList.FindObject(Form("emptyFMD%d%c_vertex%d_cent%d",d,r,v,centbin)));
+
+ return;
+
}
-
//____________________________________________________________________
void
AliPoissonCalculator::Output(TList* d)
// Reset histogram
//
if (!base) return;
- if (fBasic && fTotal && fEmpty) {
+ if (fBasic /* && fTotal && fEmpty*/) {
fBasic->Reset();
- fTotal->Reset();
- fEmpty->Reset();
+ if(!fRunningAverage) {
+ fTotal->Reset();
+ fEmpty->Reset();
+ }
return;
}
-
+ /*
Int_t nEtaF = base->GetNbinsX();
Int_t nEta = nEtaF / fEtaLumping;
Double_t etaMin = base->GetXaxis()->GetXmin();
Int_t nPhi = nPhiF / fPhiLumping;
Double_t phiMin = base->GetYaxis()->GetXmin();
Double_t phiMax = base->GetYaxis()->GetXmax();
+
- fTotal = new TH2D("total", "Total number of bins/region",
- nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
- fEmpty = new TH2D("empty", "Empty number of bins/region",
- nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
+
+
+
+ //fTotal = new TH2D("total", "Total number of bins/region",
+ // nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
+ //fEmpty = new TH2D("empty", "Empty number of bins/region",
+ //nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
fBasic = new TH2D("basic", "Basic number of hits",
nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
-
- fTotal->SetDirectory(0);
- fEmpty->SetDirectory(0);
+
+ //fTotal->SetDirectory(0);
+ //fEmpty->SetDirectory(0);
fBasic->SetDirectory(0);
- fTotal->SetXTitle("#eta");
- fEmpty->SetXTitle("#eta");
+ //fTotal->SetXTitle("#eta");
+ //fEmpty->SetXTitle("#eta");
fBasic->SetXTitle("#eta");
- fTotal->SetYTitle("#varphi [radians]");
- fEmpty->SetYTitle("#varphi [radians]");
+ //fTotal->SetYTitle("#varphi [radians]");
+ //fEmpty->SetYTitle("#varphi [radians]");
fBasic->SetYTitle("#varphi [radians]");
- fTotal->Sumw2();
- fEmpty->Sumw2();
+ //fTotal->Sumw2();
+ //fEmpty->Sumw2();
fBasic->Sumw2();
+ */
+
}
//____________________________________________________________________
void
-AliPoissonCalculator::Fill(Double_t eta, Double_t phi, Bool_t hit,
+AliPoissonCalculator::Fill(UShort_t strip, UShort_t sec, Bool_t hit,
Double_t weight)
{
//
// hit True if hit
// weight Weight if this
//
- fTotal->Fill(eta, phi);
- if (hit) fBasic->Fill(eta, phi, weight);
- else fEmpty->Fill(eta, phi);
+
+ fTotal->Fill(strip, sec);
+ if (hit) fBasic->Fill(strip, sec, weight);
+ else fEmpty->Fill(strip, sec);
}
//____________________________________________________________________
// Return:
// The result histogram (fBase overwritten)
//
+
+ // Double_t total = fEtaLumping * fPhiLumping;
+
for (Int_t ieta = 1; ieta <= fBasic->GetNbinsX(); ieta++) {
for (Int_t iphi = 1; iphi <= fBasic->GetNbinsY(); iphi++) {
Double_t eta = fBasic->GetXaxis()->GetBinCenter(ieta);
Double_t empty = fEmpty->GetBinContent(jeta, jphi);
Double_t total = fTotal->GetBinContent(jeta, jphi);
Double_t hits = fBasic->GetBinContent(ieta,iphi);
-
// Mean in region of interest
Double_t poissonM = CalculateMean(empty, total);
Double_t poissonC = CalculateCorrection(empty, total);
+
Double_t poissonV = hits * poissonM * poissonC;
Double_t poissonE = TMath::Sqrt(poissonV);
if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
Double_t corr = CalculateCorrection(empty, total);
fEmptyVsTotal->Fill(total, empty);
fMean->Fill(mean);
- fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
+ fOcc->Fill(100 * (1 - empty/total));
+ //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
fCorr->Fill(mean, corr);
}
}
#ifndef ALIPOISSONCALCULATOR_H
#define ALIPOISSONCALCULATOR_H
#include <TNamed.h>
+#include <TList.h>
class TH2D;
class TH1D;
class TBrowser;
-class TList;
+//class TList;
/**
* A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
* Constructor
*
*/
- AliPoissonCalculator(const char*);
+ AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/);
/**
* Copy constructor
*
* @param n Number of eta bins per region
*/
void SetEtaLumping(UShort_t n) { fEtaLumping = n; } //*MENU*
+ /**
+ * Set the actual object
+ *
+ * @param v vtxbin
+ */
+ void SetObject(UShort_t d, Char_t r, UShort_t v, Double_t cent); //*MENU*
/**
* Set the number of phi bins to group into a region
*
* @param etaLumping If larger than 0, set the eta lumping to this
* @param phiLumping If larger than 0, set the phi lumping to this
*/
- void Init(Int_t etaLumping=-1, Int_t phiLumping=-1);
+ void Init(UShort_t d=-1, Char_t r='I',Int_t etaLumping=-1, Int_t phiLumping=-1);
/**
* Output stuff to the passed list
*
* @param hit True if hit
* @param weight Weight if this
*/
- void Fill(Double_t eta, Double_t phi, Bool_t hit, Double_t weight=1);
+ void Fill(UShort_t strip, UShort_t sec, Bool_t hit, Double_t weight=1);
+ /**
+ * Take the average occupancy over many events or not
+ *
+ * @param use use or not
+ */
+ void SetUseAverageOverEvents(Bool_t use) {fRunningAverage = use;}
/**
* Calculate result and store in @a output
*
TH1D* fMean; // Mean calculated by poisson method
TH1D* fOcc; // Histogram of occupancies
TH2D* fCorr; // Correction as a function of mean
+ TList fTotalList; // List of total hists
+ TList fEmptyList; // List of empty hists
+ Bool_t fRunningAverage; // Whether to take average per event or not
ClassDef(AliPoissonCalculator,1) // Calculate N_ch using Poisson
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff