97 TH1 *hfit = (
TH1*)gFitter->GetObjectFit();
98 TF1 *f1 = (TF1*)gFitter->GetUserFunc();
101 npar = f1->GetNpar();
107 for (
Int_t i=0;i<nPoints;i++) {
114 for (
Int_t i=0;i<nPoints;i++) {
115 for (
Int_t j=0;j<nPoints;j++) {
116 f += df[i]*df[j]*gHistInvEMatrix->GetBinContent(i+1,j+1);
120 f1->SetNumberFitPoints(npfit);
125 return par[0]/(dm*dm+par[2]*par[2]);
138 if(rnd->Rndm()>bgr) {
178 if(rnd->Rndm()>frac) {
179 return rnd->Gaus(mTrue+smallBias,smallSigma);
181 return rnd->Gaus(mTrue+wideBias,wideSigma);
197 Double_t const luminosityData=100000;
198 Double_t const luminosityMC=1000000;
201 Int_t const nDet=250;
202 Int_t const nGen=100;
211 TH1D *histMgenMC=
new TH1D(
"MgenMC",
";mass(gen)",nGen,xminGen,xmaxGen);
212 TH1D *histMdetMC=
new TH1D(
"MdetMC",
";mass(det)",nDet,xminDet,xmaxDet);
213 TH2D *histMdetGenMC=
new TH2D(
"MdetgenMC",
";mass(det);mass(gen)",
214 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
215 Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);
216 for(
Int_t i=0;i<neventMC;i++) {
229 histMgenMC->
Fill(mGen,luminosityData/luminosityMC);
231 histMdetMC->
Fill(mDet,luminosityData/luminosityMC);
246 histMdetGenMC->
Fill(mDet,mGen,luminosityData/luminosityMC);
253 TH2D *histMdetGenSysMC=
new TH2D(
"MdetgenSysMC",
";mass(det);mass(gen)",
254 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
255 neventMC=rnd->Poisson(luminosityMC*crossSection);
256 for(
Int_t i=0;i<neventMC;i++) {
262 histMdetGenSysMC->
Fill(mDet,mGen,luminosityData/luminosityMC);
268 TH1D *histMgenData=
new TH1D(
"MgenData",
";mass(gen)",nGen,xminGen,xmaxGen);
269 TH1D *histMdetData=
new TH1D(
"MdetData",
";mass(det)",nDet,xminDet,xmaxDet);
270 Int_t neventData=rnd->Poisson(luminosityData*crossSection);
271 for(
Int_t i=0;i<neventData;i++) {
278 histMgenData->
Fill(mGen);
281 histMdetData->
Fill(mDet);
286 TH1D *histDensityGenData=
new TH1D(
"DensityGenData",
";mass(gen)",
287 nGen,xminGen,xmaxGen);
288 TH1D *histDensityGenMC=
new TH1D(
"DensityGenMC",
";mass(gen)",
289 nGen,xminGen,xmaxGen);
290 for(
Int_t i=1;i<=nGen;i++) {
308 if(unfold.SetInput(histMdetData)>=10000) {
309 std::cout<<
"Unfolding result may be wrong\n";
325 #ifdef VERBOSE_LCURVE_SCAN
331 iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
334 #ifdef VERBOSE_LCURVE_SCAN
344 TH2D *histMdetGenSys1=
new TH2D(
"Mdetgensys1",
";mass(det);mass(gen)",
345 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
346 for(
Int_t i=0;i<=nDet+1;i++) {
348 for(
Int_t j=0;j<=nGen+1;j++) {
361 std::cout<<
"tau="<<unfold.GetTau()<<
"\n";
362 std::cout<<
"chi**2="<<unfold.GetChi2A()<<
"+"<<unfold.GetChi2L()
363 <<
" / "<<unfold.GetNdf()<<
"\n";
364 std::cout<<
"chi**2(sys)="<<unfold.GetChi2Sys()<<
"\n";
371 logTauX->
GetKnot(iBest,t[0],x[0]);
372 logTauY->
GetKnot(iBest,t[0],y[0]);
373 TGraph *bestLcurve=
new TGraph(1,x,y);
374 TGraph *bestLogTauLogChi2=
new TGraph(1,t,x);
380 TH1 *histMunfold=unfold.GetOutput(
"Unfolded");
383 TH1 *histMdetFold=unfold.GetFoldedOutput(
"FoldedBack");
391 TH2 *histEmatTotal=unfold.GetEmatrixTotal(
"EmatTotal");
394 TH1D *histTotalError=
395 new TH1D(
"TotalError",
";mass(gen)",nGen,xminGen,xmaxGen);
396 for(
Int_t bin=1;bin<=nGen;bin++) {
407 TH1 *histRhoi=unfold.GetRhoItotal(
"rho_I",
421 gFitter->SetFCN(chisquare_corr);
423 TF1 *bw=
new TF1(
"bw",bw_func,xminGen,xmaxGen,3);
424 bw->SetParameter(0,1000.);
425 bw->SetParameter(1,3.8);
426 bw->SetParameter(2,0.2);
430 histMunfold->
Fit(bw,
"UE");
444 histMdetGenMC->
Draw(
"BOX");
452 histTotalError->
Draw(
"E");
454 histMunfold->
Draw(
"SAME E1");
456 histDensityGenData->
Draw(
"SAME");
457 histDensityGenMC->
Draw(
"SAME HIST");
465 histMdetFold->
Draw();
466 histMdetMC->
Draw(
"SAME HIST");
468 TH1 *histInput=unfold.GetInput(
"Minput",
";mass(det)");
471 histInput->
Draw(
"SAME");
480 bestLogTauLogChi2->SetMarkerColor(
kRed);
481 bestLogTauLogChi2->Draw(
"*");
486 bestLcurve->SetMarkerColor(
kRed);
487 bestLcurve->Draw(
"*");
489 output.SaveAs(
"testUnfold1.ps");
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Random number generator class based on M.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
R__EXTERN Int_t gErrorIgnoreLevel
R__EXTERN TStyle * gStyle
Base class for spline implementation containing the Draw/Paint methods //.
virtual Int_t GetNbinsX() const
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
static function.
virtual Double_t GetBinWidth(Int_t bin) const
return bin width for 1D historam Better to use h1.GetXaxis().GetBinWidth(bin)
This is the base class for the ROOT Random number generators.
virtual void SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetLineColor(Color_t lcolor)
virtual Double_t GetBinCenter(Int_t bin) const
return bin center for 1D historam Better to use h1.GetXaxis().GetBinCenter(bin)
Service class for 2-Dim histogram classes.
virtual void Draw(Option_t *option="")
Draw this histogram with options.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
static void RejectPoint(Bool_t reject=kTRUE)
Static function to set the global flag to reject points the fgRejectPoint global flag is tested by al...
1-D histogram with a double per channel (see TH1 documentation)}
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Abstract Base Class for Fitting.
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Double_t Sqrt(Double_t x)
static void output(int code)
Int_t Fill(Double_t)
Invalid Fill method.
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
2-D histogram with a double per channel (see TH1 documentation)}