X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=CORRFW%2FAliCFUnfolding.cxx;h=2665d9d1ee8082fdf0bc9a055b61b59eb09ddd08;hb=0341a6c63cefaadc4682ba12674f14a890ddd6f4;hp=88cec41b363af2ab0022bf0980b5aefc6bdfd7d2;hpb=c0b10ad4da792e847326a70eafdfe34d92ee2956;p=u%2Fmrichter%2FAliRoot.git diff --git a/CORRFW/AliCFUnfolding.cxx b/CORRFW/AliCFUnfolding.cxx index 88cec41b363..2665d9d1ee8 100644 --- a/CORRFW/AliCFUnfolding.cxx +++ b/CORRFW/AliCFUnfolding.cxx @@ -13,21 +13,58 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//--------------------------------------------------------------------// -// // -// AliCFUnfolding Class // -// Class to handle general unfolding procedure // -// For the moment only bayesian unfolding is supported // -// The next steps are to add chi2 minimisation and weighting methods // -// // -// Author : renaud.vernet@cern.ch // -//--------------------------------------------------------------------// +//---------------------------------------------------------------------// +// // +// AliCFUnfolding Class // +// Class to handle general unfolding procedure // +// For the moment only bayesian unfolding is supported // +// The next steps are to add chi2 minimisation and weighting methods // +// // +// // +// // +// Use : // +// ------- // +// The Bayesian unfolding consists of several iterations. // +// At each iteration, an inverse response matrix is calculated, given // +// the measured spectrum, the a priori (guessed) spectrum, // +// the efficiency spectrum and the response matrix. // +// For each iteration, the unfolded spectrum is calculated using // +// the inverse response : the goal is to get an unfolded spectrum // +// similar (according to some criterion) to the a priori one. // +// If the difference is too big, another iteration is performed : // +// the a priori spectrum is updated to the unfolded one from the // +// previous iteration, and so on so forth, until the maximum number // +// of iterations or the similarity criterion is reached. // +// // +// Currently the similarity criterion is the Chi2 between the a priori // +// and the unfolded spectrum. // +// // +// Currently the user has to define the max. number of iterations // +// (::SetMaxNumberOfIterations) // +// and the chi2 below which the procedure will stop // +// (::SetMaxChi2 or ::SetMaxChi2PerDOF) // +// // +// An optional possibility is to smooth the unfolded spectrum at the // +// end of each iteration, either using a fit function // +// (only if #dimensions <=3) // +// or a simple averaging using the neighbouring bins values. // +// This is possible calling the function ::UseSmoothing // +// If no argument is passed to this function, then the second option // +// is used. // +// // +//---------------------------------------------------------------------// +// Author : renaud.vernet@cern.ch // +//---------------------------------------------------------------------// #include "AliCFUnfolding.h" #include "TMath.h" #include "TAxis.h" #include "AliLog.h" +#include "TF1.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TH3D.h" ClassImp(AliCFUnfolding) @@ -37,13 +74,15 @@ AliCFUnfolding::AliCFUnfolding() : TNamed(), fResponse(0x0), fPrior(0x0), - fOriginalPrior(0x0), fEfficiency(0x0), fMeasured(0x0), fMaxNumIterations(0), fNVariables(0), fMaxChi2(0), fUseSmoothing(kFALSE), + fSmoothFunction(0x0), + fSmoothOption(""), + fOriginalPrior(0x0), fInverseResponse(0x0), fMeasuredEstimate(0x0), fConditional(0x0), @@ -65,13 +104,15 @@ AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const In TNamed(name,title), fResponse((THnSparse*)response->Clone()), fPrior(0x0), - fOriginalPrior(0x0), fEfficiency((THnSparse*)efficiency->Clone()), fMeasured((THnSparse*)measured->Clone()), fMaxNumIterations(0), fNVariables(nVar), fMaxChi2(0), fUseSmoothing(kFALSE), + fSmoothFunction(0x0), + fSmoothOption(""), + fOriginalPrior(0x0), fInverseResponse(0x0), fMeasuredEstimate(0x0), fConditional(0x0), @@ -91,8 +132,18 @@ AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const In else { fPrior = (THnSparse*) prior->Clone(); fOriginalPrior = (THnSparse*)fPrior->Clone(); + if (fPrior->GetNdimensions() != fNVariables) + AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); } + + if (fEfficiency->GetNdimensions() != fNVariables) + AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions())); + if (fMeasured->GetNdimensions() != fNVariables) + AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions())); + if (fResponse->GetNdimensions() != 2*fNVariables) + AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions())); + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(),iVar)); AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar)); @@ -108,13 +159,15 @@ AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) : TNamed(c), fResponse((THnSparse*)c.fResponse->Clone()), fPrior((THnSparse*)c.fPrior->Clone()), - fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()), fEfficiency((THnSparse*)c.fEfficiency->Clone()), fMeasured((THnSparse*)c.fMeasured->Clone()), fMaxNumIterations(c.fMaxNumIterations), fNVariables(c.fNVariables), fMaxChi2(c.fMaxChi2), fUseSmoothing(c.fUseSmoothing), + fSmoothFunction((TF1*)c.fSmoothFunction->Clone()), + fSmoothOption(fSmoothOption), + fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()), fInverseResponse((THnSparse*)c.fInverseResponse->Clone()), fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()), fConditional((THnSparse*)c.fConditional->Clone()), @@ -140,13 +193,15 @@ AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) { TNamed::operator=(c); fResponse = (THnSparse*)c.fResponse->Clone() ; fPrior = (THnSparse*)c.fPrior->Clone() ; - fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ; fEfficiency = (THnSparse*)c.fEfficiency->Clone() ; fMeasured = (THnSparse*)c.fMeasured->Clone() ; fMaxNumIterations = c.fMaxNumIterations ; fNVariables = c.fNVariables ; fMaxChi2 = c.fMaxChi2 ; fUseSmoothing = c.fUseSmoothing ; + fSmoothFunction = (TF1*)c.fSmoothFunction->Clone(); + fSmoothOption = c.fSmoothOption ; + fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ; fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ; fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ; fConditional = (THnSparse*)c.fConditional->Clone() ; @@ -167,9 +222,10 @@ AliCFUnfolding::~AliCFUnfolding() { // if (fResponse) delete fResponse; if (fPrior) delete fPrior; - if (fOriginalPrior) delete fOriginalPrior; if (fEfficiency) delete fEfficiency; if (fMeasured) delete fMeasured; + if (fSmoothFunction) delete fSmoothFunction; + if (fOriginalPrior) delete fOriginalPrior; if (fInverseResponse) delete fInverseResponse; if (fMeasuredEstimate) delete fMeasuredEstimate; if (fConditional) delete fConditional; @@ -206,29 +262,45 @@ void AliCFUnfolding::Init() { void AliCFUnfolding::CreateEstMeasured() { // // This function creates a estimate (M) of the reconstructed spectrum - // given the a priori distribution (T) and the conditional matrix (COND) + // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND) // // --> P(M) = SUM { P(M|T) * P(T) } - // --> M(i) = SUM_k { COND(i,k) * T(k) } + // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)} // // This is needed to calculate the inverse response matrix // // clean the measured estimate spectrum - for (Long64_t i=0; iGetNbins(); i++) { + for (Long_t i=0; iGetNbins(); i++) { fMeasuredEstimate->GetBinContent(i,fCoordinatesN_M); fMeasuredEstimate->SetBinContent(fCoordinatesN_M,0.); + fMeasuredEstimate->SetBinError (fCoordinatesN_M,0.); } + THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); + priorTimesEff->Multiply(fEfficiency); + // fill it - for (Int_t iBin=0; iBinGetNbins(); iBin++) { + for (Long_t iBin=0; iBinGetNbins(); iBin++) { Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); + Double_t conditionalError = fConditional->GetBinError (iBin); GetCoordinates(); - Double_t priorValue = fPrior->GetBinContent(fCoordinatesN_T); - Double_t fill = fMeasuredEstimate->GetBinContent(fCoordinatesN_M) + conditionalValue * priorValue * fEfficiency->GetBinContent(fCoordinatesN_T); - if (fill>0.) fMeasuredEstimate->SetBinContent(fCoordinatesN_M,fill); + Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T); + Double_t priorTimesEffError = priorTimesEff->GetBinError (fCoordinatesN_T); + Double_t fill = conditionalValue * priorTimesEffValue ; + + if (fill>0.) { + fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill); + + // error calculation : gaussian error propagation (may be overestimated...) + Double_t err2 = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ; + err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ; + Double_t err = TMath::Sqrt(err2); + fMeasuredEstimate->SetBinError(fCoordinatesN_M,err); + } } + delete priorTimesEff ; } //______________________________________________________________ @@ -242,14 +314,32 @@ void AliCFUnfolding::CreateInvResponse() { // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) } // - for (Int_t iBin=0; iBinGetNbins(); iBin++) { + THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); + priorTimesEff->Multiply(fEfficiency); + + for (Long_t iBin=0; iBinGetNbins(); iBin++) { Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); + Double_t conditionalError = fConditional->GetBinError (iBin); GetCoordinates(); - Double_t priorValue = fPrior->GetBinContent(fCoordinatesN_T); - Double_t estimatedMeasured = fMeasuredEstimate->GetBinContent(fCoordinatesN_M); - Double_t fill = (estimatedMeasured>0. ? conditionalValue * priorValue * fEfficiency->GetBinContent(fCoordinatesN_T) / estimatedMeasured : 0. ) ; - if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) fInverseResponse->SetBinContent(fCoordinates2N,fill); - } + Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M); + Double_t estMeasuredError = fMeasuredEstimate->GetBinError (fCoordinatesN_M); + Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T); + Double_t priorTimesEffError = priorTimesEff ->GetBinError (fCoordinatesN_T); + Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ; + // error calculation : gaussian error propagation (may be overestimated...) + Double_t err = 0. ; + if (estMeasuredValue>0.) { + err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) + + TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) + + TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) ) + / TMath::Power(estMeasuredValue,2) ; + } + if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) { + fInverseResponse->SetBinContent(fCoordinates2N,fill); + fInverseResponse->SetBinError (fCoordinates2N,err ); + } + } + delete priorTimesEff ; } //______________________________________________________________ @@ -269,15 +359,21 @@ void AliCFUnfolding::Unfold() { CreateInvResponse(); CreateUnfolded(); chi2 = GetChi2(); - //printf("chi2 = %e\n",chi2); + AliDebug(1,Form("Chi2 at iteration %d is %e",iIterBayes,chi2)); if (fMaxChi2>0. && chi2Clone() ; // this should be changed (memory) } - AliInfo(Form("Finished at iteration %d : Chi2 is %e and you required it to be < %e",iIterBayes,chi2,fMaxChi2)); + AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2)); } //______________________________________________________________ @@ -291,20 +387,38 @@ void AliCFUnfolding::CreateUnfolded() { // clear the unfolded spectrum - for (Long64_t i=0; iGetNbins(); i++) { + for (Long_t i=0; iGetNbins(); i++) { fUnfolded->GetBinContent(i,fCoordinatesN_T); fUnfolded->SetBinContent(fCoordinatesN_T,0.); + fUnfolded->SetBinError (fCoordinatesN_T,0.); } - for (Int_t iBin=0; iBinGetNbins(); iBin++) { + for (Long_t iBin=0; iBinGetNbins(); iBin++) { Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N); + Double_t invResponseError = fInverseResponse->GetBinError (iBin); GetCoordinates(); - Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T); - Double_t fill = fUnfolded->GetBinContent(fCoordinatesN_T) + (effValue>0. ? invResponseValue*fMeasured->GetBinContent(fCoordinatesN_M)/effValue : 0.) ; - if (fill>0.) fUnfolded->SetBinContent(fCoordinatesN_T,fill); + Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T); + Double_t effError = fEfficiency->GetBinError (fCoordinatesN_T); + Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M); + Double_t measuredError = fMeasured ->GetBinError (fCoordinatesN_M); + Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ; + + if (fill>0.) { + fUnfolded->AddBinContent(fCoordinatesN_T,fill); + + // error calculation : gaussian error propagation (may be overestimated...) + Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ; + err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ; + err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ; + err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ; + Double_t err = TMath::Sqrt(err2); + fUnfolded->SetBinError(fCoordinatesN_T,err); + } } } +//______________________________________________________________ + void AliCFUnfolding::GetCoordinates() { // // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N) @@ -327,27 +441,28 @@ void AliCFUnfolding::CreateConditional() { fConditional = (THnSparse*) fResponse->Clone(); // output of this function fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator : // projection of the response matrix on the TRUE axis - - // set in fProjResponseInT zero everywhere - for (Int_t iBin=0; iBinGetNbins(); iBin++) { - fProjResponseInT->GetBinContent(iBin,fCoordinatesN_T); - fProjResponseInT->SetBinContent(fCoordinatesN_T,0.); - } - - // calculate the response projection on T axis - for (Int_t iBin=0; iBinGetNbins(); iBin++) { - Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N); - GetCoordinates(); - Double_t fill = fProjResponseInT->GetBinContent(fCoordinatesN_T) + responseValue ; - if (fill>0.) fProjResponseInT->SetBinContent(fCoordinatesN_T,fill); - } + Int_t* dim = new Int_t [fNVariables]; + for (Int_t iDim=0; iDimProjection(fNVariables,dim,"E"); //project + delete [] dim; // fill the conditional probability matrix - for (Int_t iBin=0; iBinGetNbins(); iBin++) { + for (Long_t iBin=0; iBinGetNbins(); iBin++) { Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N); + Double_t responseError = fResponse->GetBinError (iBin); GetCoordinates(); - Double_t fill = responseValue / fProjResponseInT->GetBinContent(fCoordinatesN_T) ; - if (fill>0. || fConditional->GetBinContent(fCoordinates2N)) fConditional->SetBinContent(fCoordinates2N,fill); + Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T); + Double_t projError = fProjResponseInT->GetBinError (fCoordinatesN_T); + + Double_t fill = responseValue / projValue ; + if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) { + fConditional->SetBinContent(fCoordinates2N,fill); + // gaussian error for the moment + Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ; + Double_t err = TMath::Sqrt(err2); + err /= TMath::Power(projValue,2) ; + fConditional->SetBinError (fCoordinates2N,err); + } } } @@ -359,9 +474,10 @@ Double_t AliCFUnfolding::GetChi2() { // Double_t chi2 = 0. ; - for (Int_t iBin=0; iBinGetNbins(); iBin++) { + for (Long_t iBin=0; iBinGetNbins(); iBin++) { Double_t priorValue = fPrior->GetBinContent(iBin); - chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ; +// chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ; + chi2 += (fUnfolded->GetBinError(iBin)>0. ? TMath::Power((fUnfolded->GetBinContent(iBin) - priorValue)/fUnfolded->GetBinError(iBin),2) / priorValue : 0.) ; } return chi2; } @@ -383,20 +499,37 @@ void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) { //______________________________________________________________ -void AliCFUnfolding::Smooth() { +Short_t AliCFUnfolding::Smooth() { // // Smoothes the unfolded spectrum - // Each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins) + // + // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins) + // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn // + if (fSmoothFunction) { + AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction)); + return SmoothUsingFunction(); + } + else return SmoothUsingNeighbours(); +} + +//______________________________________________________________ + +Short_t AliCFUnfolding::SmoothUsingNeighbours() { + // + // Smoothes the unfolded spectrum using neighouring bins + // + Int_t* numBins = new Int_t[fNVariables]; for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); //need a copy because fUnfolded will be updated during the loop, and this creates problems THnSparse* copy = (THnSparse*)fUnfolded->Clone(); - for (Int_t iBin=0; iBinGetNbins(); iBin++) { //loop on non-empty bins + for (Long_t iBin=0; iBinGetNbins(); iBin++) { //loop on non-empty bins Double_t content = copy->GetBinContent(iBin,fCoordinatesN_T); + Double_t error2 = TMath::Power(copy->GetBinError(iBin),2); // skip the under/overflow bins... Bool_t isOutside = kFALSE ; @@ -413,26 +546,80 @@ void AliCFUnfolding::Smooth() { for (Int_t iVar=0; iVar 1) { // must not be on low edge border fCoordinatesN_T[iVar]-- ; //get lower neighbouring bin - Double_t contentNeighbour = copy->GetBinContent(fCoordinatesN_T); - content += contentNeighbour; + content += copy->GetBinContent(fCoordinatesN_T); + error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2); neighbours++; fCoordinatesN_T[iVar]++ ; //back to initial coordinate } if (fCoordinatesN_T[iVar] < numBins[iVar]) { // must not be on up edge border fCoordinatesN_T[iVar]++ ; //get upper neighbouring bin - Double_t contentNeighbour = copy->GetBinContent(fCoordinatesN_T); - content += contentNeighbour ; + content += copy->GetBinContent(fCoordinatesN_T); + error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2); neighbours++; fCoordinatesN_T[iVar]-- ; //back to initial coordinate } } - content /= (1+neighbours) ; // make an average - fUnfolded->SetBinContent(fCoordinatesN_T,content); + // make an average + fUnfolded->SetBinContent(fCoordinatesN_T,content/(1.+neighbours)); + fUnfolded->SetBinError (fCoordinatesN_T,TMath::Sqrt(error2)/(1.+neighbours)); } delete [] numBins; delete copy; + return 0; } +//______________________________________________________________ + +Short_t AliCFUnfolding::SmoothUsingFunction() { + // + // Fits the unfolded spectrum using the function fSmoothFunction + // + + AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar())); + + for (Int_t iPar=0; iParGetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar))); + + Int_t fitResult = 0; + + switch (fNVariables) { + case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break; + case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue + case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break; + default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1; + } + + if (fitResult != 0) { + AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult)); + return 1; + } + + Int_t nDim = fNVariables; + Int_t* bins = new Int_t[nDim]; // number of bins for each variable + Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse + + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); + nBins *= bins[iVar]; + } + + Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates) + Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3)) + + // loop on the bins and update of fUnfolded + // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin + for (Long_t iBin=0; iBinGetAxis(iVar)->GetBinCenter(bin[iVar]); + } + Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ; + fUnfolded->SetBinContent(bin,functionValue); + fUnfolded->SetBinError (bin,functionValue*fUnfolded->GetBinError(bin)); + } + return 0; +} //______________________________________________________________ @@ -446,6 +633,9 @@ void AliCFUnfolding::CreateFlatPrior() { // create the frame of the THnSparse given (for example) the one from the efficiency map fPrior = (THnSparse*) fEfficiency->Clone(); + if (fNVariables != fPrior->GetNdimensions()) + AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); + Int_t nDim = fNVariables; Int_t* bins = new Int_t[nDim]; // number of bins for each variable Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse @@ -465,6 +655,7 @@ void AliCFUnfolding::CreateFlatPrior() { bin_tmp /= bins[iVar] ; } fPrior->SetBinContent(bin,1.); // put 1 everywhere + fPrior->SetBinError (bin,0.); // put 0 everywhere } fOriginalPrior = (THnSparse*)fPrior->Clone();