#include "TH3D.h"
#include "TRandom3.h"
+
ClassImp(AliCFUnfolding)
//______________________________________________________________
AliCFUnfolding::AliCFUnfolding() :
TNamed(),
- fResponse(0x0),
- fPrior(0x0),
- fEfficiency(0x0),
- fMeasured(0x0),
+ fResponseOrig(0x0),
+ fPriorOrig(0x0),
+ fEfficiencyOrig(0x0),
fMeasuredOrig(0x0),
fMaxNumIterations(0),
fNVariables(0),
fSmoothOption("iremn"),
fMaxConvergence(0),
fNRandomIterations(0),
- fOriginalPrior(0x0),
+ fResponse(0x0),
+ fPrior(0x0),
+ fEfficiency(0x0),
+ fMeasured(0x0),
fInverseResponse(0x0),
fMeasuredEstimate(0x0),
fConditional(0x0),
fCoordinates2N(0x0),
fCoordinatesN_M(0x0),
fCoordinatesN_T(0x0),
- fRandomizedDist(0x0),
+ fRandomResponse(0x0),
+ fRandomEfficiency(0x0),
+ fRandomMeasured(0x0),
fRandom3(0x0),
fDeltaUnfoldedP(0x0),
fDeltaUnfoldedN(0x0),
Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
) :
TNamed(name,title),
- fResponse((THnSparse*)response->Clone()),
- fPrior(0x0),
- fEfficiency((THnSparse*)efficiency->Clone()),
- fMeasured((THnSparse*)measured->Clone()),
+ fResponseOrig((THnSparse*)response->Clone()),
+ fPriorOrig(0x0),
+ fEfficiencyOrig((THnSparse*)efficiency->Clone()),
fMeasuredOrig((THnSparse*)measured->Clone()),
fMaxNumIterations(maxNumIterations),
fNVariables(nVar),
fSmoothOption("iremn"),
fMaxConvergence(0),
fNRandomIterations(maxNumIterations),
- fOriginalPrior(0x0),
+ fResponse((THnSparse*)response->Clone()),
+ fPrior(0x0),
+ fEfficiency((THnSparse*)efficiency->Clone()),
+ fMeasured((THnSparse*)measured->Clone()),
fInverseResponse(0x0),
fMeasuredEstimate(0x0),
fConditional(0x0),
fCoordinates2N(0x0),
fCoordinatesN_M(0x0),
fCoordinatesN_T(0x0),
- fRandomizedDist(0x0),
+ fRandomResponse((THnSparse*)response->Clone()),
+ fRandomEfficiency((THnSparse*)efficiency->Clone()),
+ fRandomMeasured((THnSparse*)measured->Clone()),
fRandom3(0x0),
fDeltaUnfoldedP(0x0),
fDeltaUnfoldedN(0x0),
if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
else {
fPrior = (THnSparse*) prior->Clone();
- fOriginalPrior = (THnSparse*)fPrior->Clone();
+ fPriorOrig = (THnSparse*)fPrior->Clone();
if (fPrior->GetNdimensions() != fNVariables)
AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
}
AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
}
- fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
- fRandomizedDist->SetTitle("Randomized");
+ fRandomResponse ->SetTitle("Randomized response matrix");
+ fRandomEfficiency->SetTitle("Randomized efficiency");
+ fRandomMeasured ->SetTitle("Randomized measured");
SetMaxConvergencePerDOF(maxConvergencePerDOF) ;
Init();
}
-
-//______________________________________________________________
-
-AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
- TNamed(c),
- fResponse((THnSparse*)c.fResponse->Clone()),
- fPrior((THnSparse*)c.fPrior->Clone()),
- fEfficiency((THnSparse*)c.fEfficiency->Clone()),
- fMeasured((THnSparse*)c.fMeasured->Clone()),
- fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
- fMaxNumIterations(c.fMaxNumIterations),
- fNVariables(c.fNVariables),
- fUseSmoothing(c.fUseSmoothing),
- fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
- fSmoothOption(c.fSmoothOption),
- fMaxConvergence(c.fMaxConvergence),
- fNRandomIterations(c.fNRandomIterations),
- fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
- fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
- fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
- fConditional((THnSparse*)c.fConditional->Clone()),
- fUnfolded((THnSparse*)c.fUnfolded->Clone()),
- fUnfoldedFinal((THnSparse*)c.fUnfoldedFinal->Clone()),
- fCoordinates2N(new Int_t(*c.fCoordinates2N)),
- fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
- fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)),
- fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()),
- fRandom3((TRandom3*)c.fRandom3->Clone()),
- fDeltaUnfoldedP((THnSparse*)c.fDeltaUnfoldedP),
- fDeltaUnfoldedN((THnSparse*)c.fDeltaUnfoldedN),
- fNCalcCorrErrors(c.fNCalcCorrErrors),
- fRandomSeed(c.fRandomSeed)
-{
- //
- // copy constructor
- //
-}
-
-//______________________________________________________________
-
-AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
- //
- // assignment operator
- //
-
- if (this!=&c) {
- TNamed::operator=(c);
- fResponse = (THnSparse*)c.fResponse->Clone() ;
- fPrior = (THnSparse*)c.fPrior->Clone() ;
- fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
- fMeasured = (THnSparse*)c.fMeasured->Clone() ;
- fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()),
- fMaxNumIterations = c.fMaxNumIterations ;
- fNVariables = c.fNVariables ;
- fMaxConvergence = c.fMaxConvergence ;
- fUseSmoothing = c.fUseSmoothing ;
- fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
- fSmoothOption = c.fSmoothOption ;
- fNRandomIterations = c.fNRandomIterations;
- fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
- fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
- fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
- fConditional = (THnSparse*)c.fConditional->Clone() ;
- fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
- fUnfoldedFinal = (THnSparse*)c.fUnfoldedFinal->Clone() ;
- fCoordinates2N = new Int_t(*c.fCoordinates2N) ;
- fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
- fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
- fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone();
- fRandom3 = (TRandom3*)c.fRandom3->Clone();
- fDeltaUnfoldedP = (THnSparse*)c.fDeltaUnfoldedP;
- fDeltaUnfoldedN = (THnSparse*)c.fDeltaUnfoldedN;
- fNCalcCorrErrors = c.fNCalcCorrErrors ;
- fRandomSeed = c.fRandomSeed ;
- }
- return *this;
-}
-
//______________________________________________________________
AliCFUnfolding::~AliCFUnfolding() {
if (fResponse) delete fResponse;
if (fPrior) delete fPrior;
if (fEfficiency) delete fEfficiency;
+ if (fEfficiencyOrig) delete fEfficiencyOrig;
if (fMeasured) delete fMeasured;
if (fMeasuredOrig) delete fMeasuredOrig;
if (fSmoothFunction) delete fSmoothFunction;
- if (fOriginalPrior) delete fOriginalPrior;
+ if (fPriorOrig) delete fPriorOrig;
if (fInverseResponse) delete fInverseResponse;
if (fMeasuredEstimate) delete fMeasuredEstimate;
if (fConditional) delete fConditional;
if (fCoordinates2N) delete [] fCoordinates2N;
if (fCoordinatesN_M) delete [] fCoordinatesN_M;
if (fCoordinatesN_T) delete [] fCoordinatesN_T;
- if (fRandomizedDist) delete fRandomizedDist;
+ if (fRandomResponse) delete fRandomResponse;
+ if (fRandomEfficiency) delete fRandomEfficiency;
+ if (fRandomMeasured) delete fRandomMeasured;
if (fRandom3) delete fRandom3;
if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
if (fDeltaUnfoldedN) delete fDeltaUnfoldedN;
+
}
//______________________________________________________________
fUnfolded->SetTitle("Unfolded");
// create the frame of the measurement estimate spectrum
fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
- // create the frame of the original measurement spectrum
- fMeasuredOrig = (THnSparse*) fMeasured->Clone();
+ // create the frame of the delta profiles
fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
fDeltaUnfoldedP->SetTitle("#Delta unfolded");
fDeltaUnfoldedP->Reset();
fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
- fDeltaUnfoldedP->SetTitle("");
+ fDeltaUnfoldedN->SetTitle("");
fDeltaUnfoldedN->Reset();
+
+
}
void AliCFUnfolding::CalculateCorrelatedErrors() {
- // Step 1: Create randomized distribution (fRandomizedDist) of each bin of
+ // Step 1: Create randomized distribution (fRandomXXXX) of each bin of
// the measured spectrum to calculate correlated errors.
// Poisson statistics: mean = measured value of bin
// Step 2: Unfold randomized distribution
// reset prior to original one
if (fPrior) delete fPrior ;
- fPrior = (THnSparse*) fOriginalPrior->Clone();
+ fPrior = (THnSparse*) fPriorOrig->Clone();
// create randomized distribution and stick measured spectrum to it
CreateRandomizedDist();
- if (fMeasured) delete fMeasured ;
- fMeasured = (THnSparse*) fRandomizedDist->Clone();
+
+ if (fResponse) delete fResponse ;
+ fResponse = (THnSparse*) fRandomResponse->Clone();
+ fResponse->SetTitle("Response");
+
+ if (fEfficiency) delete fEfficiency ;
+ fEfficiency = (THnSparse*) fRandomEfficiency->Clone();
+ fEfficiency->SetTitle("Efficiency");
+
+ if (fMeasured) delete fMeasured ;
+ fMeasured = (THnSparse*) fRandomMeasured->Clone();
fMeasured->SetTitle("Measured");
- //unfold fRandomizedDist
+ //unfold with randomized distributions
Unfold();
FillDeltaUnfoldedProfile();
}
// Get statistical errors for final unfolded spectrum
// ie. spread of each pt bin in fDeltaUnfoldedP
- Double_t sigma = 0.;
- Double_t dummy = 0.;
+ Double_t meanx2 = 0.;
+ Double_t mean = 0.;
+ Double_t checksigma = 0.;
+ Double_t entriesInBin = 0.;
for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
- dummy = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
- sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
+ fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
+ mean = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M);
+ meanx2 = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
+ entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
+ if(entriesInBin > 1.) checksigma = TMath::Sqrt((entriesInBin/(entriesInBin-1.))*TMath::Abs(meanx2-mean*mean));
+ //printf("mean %f, meanx2 %f, sigmacheck %f, nentries %f\n",mean, meanx2, checksigma,entriesInBin);
//AliDebug(2,Form("filling error %e\n",sigma));
- fUnfoldedFinal->SetBinError(fCoordinatesN_M,sigma);
+ fUnfoldedFinal->SetBinError(fCoordinatesN_M,checksigma);
}
// now errors are calculated
// This distribution is created several times, each time with a different random number
//
- Double_t random = 0.;
- Double_t measuredValue = 0.;
- Double_t measuredError = 0.;
-
- for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) {
- measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
- measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
+ for (Long_t iBin=0; iBin<fResponseOrig->GetNbins(); iBin++) {
+ Double_t val = fResponseOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
+ Double_t err = fResponseOrig->GetBinError(fCoordinatesN_M); //used as sigma
+ Double_t ran = fRandom3->Gaus(val,err);
// random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
- random = fRandom3->Gaus(measuredValue,measuredError);
- fRandomizedDist->SetBinContent(iBin,random);
+ fRandomResponse->SetBinContent(iBin,ran);
+ }
+ for (Long_t iBin=0; iBin<fEfficiencyOrig->GetNbins(); iBin++) {
+ Double_t val = fEfficiencyOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
+ Double_t err = fEfficiencyOrig->GetBinError(fCoordinatesN_M); //used as sigma
+ Double_t ran = fRandom3->Gaus(val,err);
+ // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
+ fRandomEfficiency->SetBinContent(iBin,ran);
+ }
+ for (Long_t iBin=0; iBin<fMeasuredOrig->GetNbins(); iBin++) {
+ Double_t val = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
+ Double_t err = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
+ Double_t ran = fRandom3->Gaus(val,err);
+ // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
+ fRandomMeasured->SetBinContent(iBin,ran);
}
}
// mean_{n+1} = (n*mean_n + value_{n+1}) / (n+1)
// sigma_{n+1} = sqrt { 1/(n+1) * [ n*sigma_n^2 + (n^2+n)*(mean_{n+1}-mean_n)^2 ] } (can this be optimized?)
- for (Long_t iBin=0; iBin<fUnfolded->GetNbins(); iBin++) {
- Double_t deltaInBin = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(iBin);
+ for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
+ Double_t deltaInBin = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(fCoordinatesN_M);
Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
//AliDebug(2,Form("%e %e ==> delta = %e\n",fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M),fUnfolded->GetBinContent(iBin),deltaInBin));
+ //printf("deltaInBin %f\n",deltaInBin);
+ //printf("pt %f\n",ptaxis->GetBinCenter(iBin+1));
+
Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
Double_t mean_nplus1 = mean_n ;
mean_nplus1 *= entriesInBin ;
mean_nplus1 += deltaInBin ;
mean_nplus1 /= (entriesInBin+1) ;
- Double_t sigma = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
- sigma *= sigma ;
- sigma *= entriesInBin ;
- sigma += ( (entriesInBin*entriesInBin+entriesInBin) * TMath::Power(mean_nplus1 - mean_n,2) ) ;
- sigma /= (entriesInBin+1) ;
- sigma = TMath::Sqrt(sigma) ;
+ Double_t meanx2_n = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M) ;
+ Double_t meanx2_nplus1 = meanx2_n ;
+ meanx2_nplus1 *= entriesInBin ;
+ meanx2_nplus1 += (deltaInBin*deltaInBin) ;
+ meanx2_nplus1 /= (entriesInBin+1) ;
//AliDebug(2,Form("sigma = %e\n",sigma));
+ fDeltaUnfoldedP->SetBinError(fCoordinatesN_M,meanx2_nplus1) ;
fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
- fDeltaUnfoldedP->SetBinError (fCoordinatesN_M,sigma) ;
fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
}
}
// --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
//
- fConditional = (THnSparse*) fResponse->Clone(); // output of this function
+ fConditional = (THnSparse*) fResponse->Clone(); // output of this function
Int_t* dim = new Int_t [fNVariables];
for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
fPrior->SetBinError (bin,0.); // put 0 everywhere
}
- fOriginalPrior = (THnSparse*)fPrior->Clone();
+ fPriorOrig = (THnSparse*)fPrior->Clone();
delete [] bin;
delete [] bins;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff