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()),
+ fResponseOrig((THnSparse*)c.fResponseOrig->Clone()),
+ fPriorOrig((THnSparse*)c.fPriorOrig->Clone()),
+ fEfficiencyOrig((THnSparse*)c.fEfficiencyOrig->Clone()),
fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
fMaxNumIterations(c.fMaxNumIterations),
fNVariables(c.fNVariables),
fSmoothOption(c.fSmoothOption),
fMaxConvergence(c.fMaxConvergence),
fNRandomIterations(c.fNRandomIterations),
- fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
+ fResponse((THnSparse*)c.fResponse->Clone()),
+ fPrior((THnSparse*)c.fPrior->Clone()),
+ fEfficiency((THnSparse*)c.fEfficiency->Clone()),
+ fMeasured((THnSparse*)c.fMeasured->Clone()),
fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
fConditional((THnSparse*)c.fConditional->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()),
+ fRandomResponse((THnSparse*)c.fRandomResponse->Clone()),
+ fRandomEfficiency((THnSparse*)c.fRandomEfficiency->Clone()),
+ fRandomMeasured((THnSparse*)c.fRandomMeasured->Clone()),
fRandom3((TRandom3*)c.fRandom3->Clone()),
fDeltaUnfoldedP((THnSparse*)c.fDeltaUnfoldedP),
fDeltaUnfoldedN((THnSparse*)c.fDeltaUnfoldedN),
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()),
+ fResponseOrig = (THnSparse*)c.fResponseOrig->Clone() ;
+ fPriorOrig = (THnSparse*)c.fPriorOrig->Clone() ;
+ fEfficiencyOrig = (THnSparse*)c.fEfficiencyOrig->Clone() ;
+ fMeasuredOrig = (THnSparse*)c.fMeasuredOrig->Clone() ;
fMaxNumIterations = c.fMaxNumIterations ;
fNVariables = c.fNVariables ;
- fMaxConvergence = c.fMaxConvergence ;
fUseSmoothing = c.fUseSmoothing ;
- fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
+ fSmoothFunction = (TF1*)c.fSmoothFunction->Clone() ;
fSmoothOption = c.fSmoothOption ;
- fNRandomIterations = c.fNRandomIterations;
- fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
+ fMaxConvergence = c.fMaxConvergence ;
+ fNRandomIterations = c.fNRandomIterations ;
+ fResponse = (THnSparse*)c.fResponse->Clone() ;
+ fPrior = (THnSparse*)c.fPrior->Clone() ;
+ fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
+ fMeasured = (THnSparse*)c.fMeasured->Clone() ;
fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
- fConditional = (THnSparse*)c.fConditional->Clone() ;
+ fConditional = (THnSparse*)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;
+ fRandomResponse = (THnSparse*)c.fRandomResponse->Clone() ;
+ fRandomEfficiency = (THnSparse*)c.fRandomEfficiency->Clone() ;
+ fRandomMeasured = (THnSparse*)c.fRandomMeasured->Clone() ;
+ fRandom3 = (TRandom3*)c.fRandom3->Clone() ;
+ fDeltaUnfoldedP = (THnSparse*)c.fDeltaUnfoldedP->Clone() ;
+ fDeltaUnfoldedN = (THnSparse*)c.fDeltaUnfoldedN->Clone() ;
fNCalcCorrErrors = c.fNCalcCorrErrors ;
fRandomSeed = c.fRandomSeed ;
}
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();
}
// 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)
+ 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)
- random = fRandom3->Gaus(measuredValue,measuredError);
- fRandomizedDist->SetBinContent(iBin,random);
+ fRandomMeasured->SetBinContent(iBin,ran);
}
}
// --> 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;
class AliCFUnfolding : public TNamed {
public :
- enum EUnfoldingErrorStatus {
- kNotDone = BIT(0), // Doing the normal unfolding, before errors are calculated
- kBeingDone = BIT(1), // In the process of calculating errors
- kDone = BIT(2) // Errors have been calculated
- };
AliCFUnfolding();
AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
void Unfold();
- THnSparse* GetResponse() const {return fResponse;}
- THnSparse* GetInverseResponse() const {return fInverseResponse;}
- THnSparse* GetPrior() const {return fPrior;}
- THnSparse* GetOriginalPrior() const {return fOriginalPrior;}
- THnSparse* GetEfficiency() const {return fEfficiency;}
- THnSparse* GetUnfolded() const {return fUnfoldedFinal;}
- THnSparse* GetEstMeasured() const {return fMeasuredEstimate;}
- THnSparse* GetMeasured() const {return fMeasured;}
- THnSparse* GetConditional() const {return fConditional;}
- TF1* GetSmoothFunction() const {return fSmoothFunction;}
- THnSparse* GetDeltaUnfoldedProfile() const {return fDeltaUnfoldedP;}
- Int_t GetDOF(); // Returns number of degrees of freedom
+ const THnSparse* GetResponse() const {return fResponseOrig;}
+ const THnSparse* GetEfficiency() const {return fEfficiencyOrig;}
+ const THnSparse* GetMeasured() const {return fMeasuredOrig;}
+ const THnSparse* GetOriginalPrior() const {return fPriorOrig;}
+ THnSparse* GetInverseResponse() const {return fInverseResponse;}
+ THnSparse* GetPrior() const {return fPrior;}
+ THnSparse* GetUnfolded() const {return fUnfoldedFinal;}
+ THnSparse* GetEstMeasured() const {return fMeasuredEstimate;}
+ THnSparse* GetConditional() const {return fConditional;}
+ TF1* GetSmoothFunction() const {return fSmoothFunction;}
+ THnSparse* GetDeltaUnfoldedProfile() const {return fDeltaUnfoldedP;}
+ Int_t GetDOF(); // Returns number of degrees of freedom
static Short_t SmoothUsingNeighbours(THnSparse*); // smoothes the unfolded spectrum using the neighbouring cells
private :
+ //
// user-related settings
- THnSparse *fResponse; // Response matrix : dimensions must be 2N = 2 x (number of variables)
- // dimensions 0 -> N-1 must be filled with reconstructed values
- // dimensions N -> 2N-1 must be filled with generated values
- THnSparse *fPrior; // This is the assumed generated distribution : dimensions must be N = number of variables
- // it will be used at the first step
- // then will be updated automatically at each iteration
- THnSparse *fEfficiency; // Efficiency map : dimensions must be N = number of variables
- // this map must be filled only with "true" values of the variables (should not include resolution effects)
- THnSparse *fMeasured; // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)
- THnSparse *fMeasuredOrig; // Measured spectrum to be unfolded : dimensions must be N = number of variables (not modified)
- Int_t fMaxNumIterations; // Maximum number of iterations to be performed
- Int_t fNVariables; // Number of variables used in analysis spectra (pt, y, ...)
- Bool_t fUseSmoothing; // Smooth the unfolded sectrum at each iteration; default is kFALSE
- TF1 *fSmoothFunction; // Function used to smooth the unfolded spectrum
- Option_t *fSmoothOption; // Option to use during the fit (with fSmoothFunction) ; default is "iremn"
-
- /* correlated error calculation */
- Double_t fMaxConvergence; // Convergence criterion in case of correlated error calculation
- Int_t fNRandomIterations; // Number of random distributed measured spectra
-
+ //
+ const THnSparse *fResponseOrig; // Response matrix : dimensions must be 2N = 2 x (number of variables)
+ // dimensions 0 -> N-1 must be filled with reconstructed values
+ // dimensions N -> 2N-1 must be filled with generated values
+ const THnSparse *fPriorOrig; // This is the assumed generated distribution : dimensions must be N = number of variables
+ // it will be used at the first step
+ // then will be updated automatically at each iteration
+ const THnSparse *fEfficiencyOrig; // Efficiency map : dimensions must be N = number of variables (modified)
+ // this map must be filled only with "true" values of the variables (should not do "bin smearing")
+ const THnSparse *fMeasuredOrig; // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)
+
+ Int_t fMaxNumIterations; // Maximum number of iterations to be performed
+ Int_t fNVariables; // Number of variables used in analysis spectra (pt, y, ...)
+ Bool_t fUseSmoothing; // Smooth the unfolded sectrum at each iteration; default is kFALSE
+ TF1 *fSmoothFunction; // Function used to smooth the unfolded spectrum
+ Option_t *fSmoothOption; // Option to use during the fit (with fSmoothFunction) ; default is "iremn"
+
+ //
// internal settings
- THnSparse *fOriginalPrior; // This is the original prior distribution : will not be modified
- THnSparse *fInverseResponse; // Inverse response matrix
- THnSparse *fMeasuredEstimate; // Estimation of the measured (M) spectrum given the a priori (T) distribution
- THnSparse *fConditional; // Matrix holding the conditional probabilities P(M|T)
- THnSparse *fUnfolded; // Unfolded spectrum (modified before and during error calculation)
- THnSparse *fUnfoldedFinal; // Final unfolded spectrum
- Int_t *fCoordinates2N; // Coordinates in 2N (measured,true) space
- Int_t *fCoordinatesN_M; // Coordinates in measured space
- Int_t *fCoordinatesN_T; // Coordinates in true space
+ //
+ Double_t fMaxConvergence; // Convergence criterion in case of correlated error calculation
+ Int_t fNRandomIterations; // Number of random distributed measured spectra
+ THnSparse *fResponse; // Copy of the original response matrix (modified)
+ THnSparse *fPrior; // Copy of the original prior spectrum (modified)
+ THnSparse *fEfficiency; // Copy of original efficiency (modified)
+ THnSparse *fMeasured; // Copy of the original measureed spectrum (modified)
+ THnSparse *fInverseResponse; // Inverse response matrix
+ THnSparse *fMeasuredEstimate; // Estimation of the measured (M) spectrum given the a priori (T) distribution
+ THnSparse *fConditional; // Matrix holding the conditional probabilities P(M|T)
+ THnSparse *fUnfolded; // Unfolded spectrum (modified before and during error calculation)
+ THnSparse *fUnfoldedFinal; // Final unfolded spectrum
+ Int_t *fCoordinates2N; // Coordinates in 2N (measured,true) space
+ Int_t *fCoordinatesN_M; // Coordinates in measured space
+ Int_t *fCoordinatesN_T; // Coordinates in true space
/* correlated error calculation */
- THnSparse *fRandomizedDist; // Randomized distribution for each bin of the measured spectrum to calculate correlated errors
+ THnSparse *fRandomResponse; // Randomized distribution for each bin of the response matrix to calculate correlated errors
+ THnSparse *fRandomEfficiency; // Randomized distribution for each bin of the efficiency spectrum to calculate correlated errors
+ THnSparse *fRandomMeasured; // Randomized distribution for each bin of the measured spectrum to calculate correlated errors
TRandom3 *fRandom3; // Object to get random number following Poisson distribution
THnSparse *fDeltaUnfoldedP; // Profile of the delta-unfolded distribution
THnSparse *fDeltaUnfoldedN; // Entries of the delta-unfolded distribution (count for each bin)