AliPhysics  ad6828d (ad6828d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliHFMultiTrials.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 2008-2019, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TMath.h>
17 #include <TPad.h>
18 #include <TCanvas.h>
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <TNtuple.h>
22 #include <TF1.h>
23 #include <TLatex.h>
24 #include <TFile.h>
25 #include "AliHFMassFitter.h"
26 #include "AliHFMassFitterVAR.h"
27 #include "AliHFMultiTrials.h"
28 
32 
33 
34 //_________________________________________________________________________
36  TNamed(),
37  fNumOfRebinSteps(4),
38  fRebinSteps(0x0),
39  fNumOfFirstBinSteps(1),
40  fNumOfLowLimFitSteps(6),
41  fLowLimFitSteps(0x0),
42  fNumOfUpLimFitSteps(6),
43  fUpLimFitSteps(0x0),
44  fNumOfnSigmaBinCSteps(11),
45  fnSigmaBinCSteps(0x0),
46  fnSigmaForBkgEval(3),
47  fSigmaGausMC(0.010),
48  fSigmaMCVariation(0.15),
49  fMassD(1.86484),
50  fSuffix(""),
51  fFitOption(0),
52  fUseExpoBkg(kTRUE),
53  fUseLinBkg(kTRUE),
54  fUsePol2Bkg(kTRUE),
55  fUsePol3Bkg(kTRUE),
56  fUsePol4Bkg(kTRUE),
57  fUsePol5Bkg(kFALSE),
58  fUsePowLawBkg(kFALSE),
59  fUsePowLawTimesExpoBkg(kFALSE),
60  fUseFixSigUpFreeMean(kTRUE),
61  fUseFixSigDownFreeMean(kTRUE),
62  fUseFreeS(kTRUE),
63  fUseFixedMeanFreeS(kTRUE),
64  fUseFixSigFreeMean(kTRUE),
65  fUseFixSigFixMean(kTRUE),
66  fSaveBkgVal(kFALSE),
67  fDrawIndividualFits(kFALSE),
68  fHistoRawYieldDistAll(0x0),
69  fHistoRawYieldTrialAll(0x0),
70  fHistoSigmaTrialAll(0x0),
71  fHistoMeanTrialAll(0x0),
72  fHistoChi2TrialAll(0x0),
73  fHistoSignifTrialAll(0x0),
74  fHistoBkgTrialAll(0x0),
75  fHistoBkgInBinEdgesTrialAll(0x0),
76  fHistoRawYieldDistBinCAll(0x0),
77  fHistoRawYieldTrialBinCAll(0x0),
78  fHistoRawYieldDist(0x0),
79  fHistoRawYieldTrial(0x0),
80  fHistoSigmaTrial(0x0),
81  fHistoMeanTrial(0x0),
82  fHistoChi2Trial(0x0),
83  fHistoSignifTrial(0x0),
84  fHistoBkgTrial(0x0),
85  fHistoBkgInBinEdgesTrial(0x0),
86  fHistoRawYieldDistBinC(0x0),
87  fHistoRawYieldTrialBinC(0x0),
88  fhTemplRefl(0x0),
89  fFixRefloS(0),
90  fNtupleMultiTrials(0x0),
91  fMinYieldGlob(0),
92  fMaxYieldGlob(0)
93 {
94  // constructor
95  Int_t rebinStep[4]={3,4,5,6};
96  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
97  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
98  Double_t nSigmasBC[11]={2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
99  ConfigureRebinSteps(4,rebinStep);
100  ConfigureLowLimFitSteps(6,minMassStep);
101  ConfigureUpLimFitSteps(6,maxMassStep);
102  ConfigurenSigmaBinCSteps(11,nSigmasBC);
103 }
104 
105 //________________________________________________________________________
107  // destructor
108  delete [] fRebinSteps;
109  delete [] fLowLimFitSteps;
110  delete [] fUpLimFitSteps;
111  if(fhTemplRefl) delete fhTemplRefl;
112 }
113 
114 //________________________________________________________________________
116  // creates output histograms
117 
118  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
119 
120  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
121  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
122 
124 
125  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
126  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
127  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
128  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
129  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
130  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
131  if(fSaveBkgVal) {
132  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
133  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
134  }
135 
136 
137  fHistoRawYieldDistBinCAll = new TH1F(Form("hRawYieldDistBinCAll%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
138  fHistoRawYieldTrialBinCAll = new TH2F(Form("hRawYieldTrialBinCAll%s",fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
139 
140  fHistoRawYieldDist = new TH1F*[nCases];
141  fHistoRawYieldTrial = new TH1F*[nCases];
142  fHistoSigmaTrial = new TH1F*[nCases];
143  fHistoMeanTrial = new TH1F*[nCases];
144  fHistoChi2Trial = new TH1F*[nCases];
145  fHistoSignifTrial = new TH1F*[nCases];
146  if(fSaveBkgVal) {
147  fHistoBkgTrial = new TH1F*[nCases];
148  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
149  }
150 
151  fHistoRawYieldDistBinC = new TH1F*[nCases];
152  fHistoRawYieldTrialBinC = new TH2F*[nCases];
153 
154  for(Int_t ib=0; ib<kNBkgFuncCases; ib++){
155  for(Int_t igs=0; igs<kNFitConfCases; igs++){
156  Int_t theCase=igs*kNBkgFuncCases+ib;
157  fHistoRawYieldDist[theCase]=new TH1F(Form("hRawYieldDist%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
158  fHistoRawYieldDistBinC[theCase]=new TH1F(Form("hRawYieldDistBinC%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
159  fHistoRawYieldTrial[theCase]=new TH1F(Form("hRawYieldTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Raw Yield",totTrials,-0.5,totTrials-0.5);
160  fHistoRawYieldTrialBinC[theCase]=new TH2F(Form("hRawYieldTrialBinC%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
161  fHistoSigmaTrial[theCase]=new TH1F(Form("hSigmaTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
162  fHistoMeanTrial[theCase]=new TH1F(Form("hMeanTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
163  fHistoChi2Trial[theCase]=new TH1F(Form("hChi2Trial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; #chi^{2}",totTrials,-0.5,totTrials-0.5);
164  fHistoSignifTrial[theCase]=new TH1F(Form("hSignifTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Significance",totTrials,-0.5,totTrials-0.5);
165  if(fSaveBkgVal) {
166  fHistoBkgTrial[theCase] = new TH1F(Form("hBkgTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background",totTrials,-0.5,totTrials-0.5);
167  fHistoBkgInBinEdgesTrial[theCase] = new TH1F(Form("hBkgInBinEdgesTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background in bin edges",totTrials,-0.5,totTrials-0.5);
168  }
169 
170  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
171  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
172  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
173  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
174  if(fSaveBkgVal) {
175  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
176  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
177  }
178 
179  }
180  }
181  fNtupleMultiTrials = new TNtuple(Form("ntuMultiTrial%s",fSuffix.Data()),Form("ntuMultiTrial%s",fSuffix.Data()),"rebin:firstb:minfit:maxfit:bkgfunc:confsig:confmean:chi2:signif:mean:emean:sigma:esigma:rawy:erawy",128000);
182  return kTRUE;
183 
184 }
185 
186 //________________________________________________________________________
187 Bool_t AliHFMultiTrials::DoMultiTrials(TH1D* hInvMassHisto, TPad* thePad){
188  // perform the multiple fits
189 
190  Bool_t hOK=CreateHistos();
191  if(!hOK) return kFALSE;
192 
193  Int_t itrial=0;
194  Int_t types=0;
195  Int_t itrialBC=0;
197 
198  fMinYieldGlob=999999.;
199  fMaxYieldGlob=0.;
200  Float_t xnt[15];
201 
202  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
203  Int_t rebin=fRebinSteps[ir];
204  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
205  TH1F* hRebinned=0x0;
206  if(fNumOfFirstBinSteps==1) hRebinned=RebinHisto(hInvMassHisto,rebin,-1);
207  else hRebinned=RebinHisto(hInvMassHisto,rebin,iFirstBin);
208  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
210  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
211  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
213  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
214  ++itrial;
215  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
216  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
217  if(typeb==kLinBkg && !fUseLinBkg) continue;
218  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
219  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
220  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
221  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
222  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
224  for(Int_t igs=0; igs<kNFitConfCases; igs++){
225  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
226  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
227  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
228  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
229  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
230  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
231  Int_t theCase=igs*kNBkgFuncCases+typeb;
232  Int_t globBin=itrial+theCase*totTrials;
233  for(Int_t j=0; j<15; j++) xnt[j]=0.;
234 
236  if(typeb<=kPol2Bkg){
237  fitter=new AliHFMassFitterVAR(hRebinned,hmin, hmax,1,typeb,types);
238  }else if(typeb==kPowBkg){
239  fitter=new AliHFMassFitterVAR(hRebinned,hmin, hmax,1,4,types);
240  }else if(typeb==kPowTimesExpoBkg){
241  fitter=new AliHFMassFitterVAR(hRebinned,hmin, hmax,1,5,types);
242  }else{
243  fitter=new AliHFMassFitterVAR(hRebinned,hmin, hmax,1,6,types);
244  if(typeb==kPol3Bkg) fitter->SetBackHighPolDegree(3);
245  if(typeb==kPol4Bkg) fitter->SetBackHighPolDegree(4);
246  if(typeb==kPol5Bkg) fitter->SetBackHighPolDegree(5);
247  }
248  fitter->SetReflectionSigmaFactor(0);
249  //if D0 Reflection
250  if(fhTemplRefl){
251  fitter=new AliHFMassFitterVAR(hRebinned,hmin,hmax,1,typeb,2);
253  fitter->SetFixReflOverS(fFixRefloS,kTRUE);
254  }
255  if(fFitOption==1) fitter->SetUseChi2Fit();
258  xnt[0]=rebin;
259  xnt[1]=iFirstBin;
260  xnt[2]=minMassForFit;
261  xnt[3]=maxMassForFit;
262  xnt[4]=typeb;
263  xnt[6]=0;
264  if(igs==kFixSigFreeMean){
265  fitter->SetFixGaussianSigma(fSigmaGausMC,kTRUE);
266  xnt[5]=1;
267  }else if(igs==kFixSigUpFreeMean){
269  xnt[5]=2;
270  }else if(igs==kFixSigDownFreeMean){
272  xnt[5]=3;
273  }else if(igs==kFreeSigFreeMean){
274  xnt[5]=0;
275  }else if(igs==kFixSigFixMean){
276  fitter->SetFixGaussianSigma(fSigmaGausMC,kTRUE);
277  fitter->SetFixGaussianMean(fMassD,kTRUE);
278  xnt[5]=1;
279  xnt[6]=1;
280  }else if(igs==kFreeSigFixMean){
281  fitter->SetFixGaussianMean(fMassD,kTRUE);
282  xnt[5]=0;
283  xnt[6]=1;
284  }
285  Bool_t out=kFALSE;
286  Double_t chisq=-1.;
287  Double_t sigma=0.;
288  Double_t esigma=0.;
289  Double_t pos=.0;
290  Double_t epos=.0;
291  Double_t ry=.0;
292  Double_t ery=.0;
293  Double_t significance=0.;
294  Double_t erSignif=0.;
295  Double_t bkg=0.;
296  Double_t erbkg=0.;
297  Double_t bkgBEdge=0;
298  Double_t erbkgBEdge=0;
299  TF1* fB1=0x0;
300  if(typeb<kNBkgFuncCases){
301  printf("****** START FIT OF HISTO %s WITH REBIN %d FIRST BIN %d MASS RANGE %f-%f BACKGROUND FIT FUNCTION=%d CONFIG SIGMA/MEAN=%d\n",hInvMassHisto->GetName(),rebin,iFirstBin,minMassForFit,maxMassForFit,typeb,igs);
302  out=fitter->MassFitter(0);
303  chisq=fitter->GetReducedChiSquare();
304  fitter->Significance(3,significance,erSignif);
305  sigma=fitter->GetSigma();
306  pos=fitter->GetMean();
307  esigma=fitter->GetSigmaUncertainty();
308  if(esigma<0.00001) esigma=0.0001;
309  epos=fitter->GetMeanUncertainty();
310  if(epos<0.00001) epos=0.0001;
311  ry=fitter->GetRawYield();
312  ery=fitter->GetRawYieldError();
313  fB1=fitter->GetBackgroundFullRangeFunc();
314  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
315  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
316  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
317  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
318  if(out && fDrawIndividualFits && thePad){
319  thePad->Clear();
320  fitter->DrawHere(thePad);
321  for (auto format : fInvMassFitSaveAsFormats) {
322  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
323  }
324  }
325  }
326  // else{
327  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
328  // if(out && thePad){
329  // thePad->Clear();
330  // hRebinned->Draw();
331  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
332  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
333  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
334  // fB1->SetLineColor(2);
335  // fB1->Draw("same");
336  // fSB->SetLineColor(4);
337  // fSB->Draw("same");
338  // thePad->Update();
339  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
340  // sigma=fSB->GetParameter(2);
341  // esigma=fSB->GetParError(2);
342  // if(esigma<0.00001) esigma=0.0001;
343  // pos=fSB->GetParameter(1);
344  // epos=fSB->GetParError(1);
345  // if(epos<0.00001) epos=0.0001;
346  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
347  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
348  // }
349  // }
350  xnt[7]=chisq;
351  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
352  xnt[8]=significance;
353  xnt[9]=pos;
354  xnt[10]=epos;
355  xnt[11]=sigma;
356  xnt[12]=esigma;
357  xnt[13]=ry;
358  xnt[14]=ery;
359  fHistoRawYieldDistAll->Fill(ry);
360  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
361  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
362  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
363  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
364  fHistoMeanTrialAll->SetBinContent(globBin,pos);
365  fHistoMeanTrialAll->SetBinError(globBin,epos);
366  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
367  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
368  fHistoSignifTrialAll->SetBinContent(globBin,significance);
369  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
370  if(fSaveBkgVal) {
371  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
372  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
373  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
374  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
375  }
376 
377  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
378  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
379  fHistoRawYieldDist[theCase]->Fill(ry);
380  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
381  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
382  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
383  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
384  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
385  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
386  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
387  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
388  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
389  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
390  if(fSaveBkgVal) {
391  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
392  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
393  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
394  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
395  }
396 
397  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
398  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
399  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
400  if(minMassBC>minMassForFit &&
401  maxMassBC<maxMassForFit &&
402  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
403  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
404  Double_t cnts,ecnts;
405  BinCount(hRebinned,fB1,1,minMassBC,maxMassBC,cnts,ecnts);
406  ++itrialBC;
407  fHistoRawYieldDistBinCAll->Fill(cnts);
408  fHistoRawYieldTrialBinCAll->SetBinContent(globBin,iStepBC+1,cnts);
409  fHistoRawYieldTrialBinCAll->SetBinError(globBin,iStepBC+1,ecnts);
410  fHistoRawYieldTrialBinC[theCase]->SetBinContent(itrial,iStepBC+1,cnts);
411  fHistoRawYieldTrialBinC[theCase]->SetBinError(itrial,iStepBC+1,ecnts);
412  fHistoRawYieldDistBinC[theCase]->Fill(cnts);
413  }
414  }
415  }
416  delete fitter;
417  // if(typeb>4) delete fB1;
418  fNtupleMultiTrials->Fill(xnt);
419  }
420  }
421  }
422  }
423  delete hRebinned;
424  }
425  }
426  return kTRUE;
427 }
428 
429 //________________________________________________________________________
431  // save histos in a root file for further analysis
432  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
433  TFile* outHistos=new TFile(fileName.Data(),option.Data());
434  fHistoRawYieldTrialAll->Write();
435  fHistoSigmaTrialAll->Write();
436  fHistoMeanTrialAll->Write();
437  fHistoChi2TrialAll->Write();
438  fHistoSignifTrialAll->Write();
439  if(fSaveBkgVal) {
440  fHistoBkgTrialAll->Write();
441  fHistoBkgInBinEdgesTrialAll->Write();
442  }
443  fHistoRawYieldDistBinCAll->Write();
444  fHistoRawYieldTrialBinCAll->Write();
445  for(Int_t ic=0; ic<nCases; ic++){
446  fHistoRawYieldTrial[ic]->Write();
447  fHistoSigmaTrial[ic]->Write();
448  fHistoMeanTrial[ic]->Write();
449  fHistoChi2Trial[ic]->Write();
450  fHistoSignifTrial[ic]->Write();
451  if(fSaveBkgVal) {
452  fHistoBkgTrial[ic]->Write();
453  fHistoBkgInBinEdgesTrial[ic]->Write();
454  }
455  fHistoRawYieldTrialBinC[ic]->Write();
456  fHistoRawYieldDistBinC[ic]->Write();
457  }
458  fNtupleMultiTrials->Write();
459  outHistos->Close();
460  delete outHistos;
461 }
462 
463 //________________________________________________________________________
464 void AliHFMultiTrials::DrawHistos(TCanvas* cry) const{
465  // draw histos
466  cry->Divide(2,2);
467  cry->cd(1);
468  fHistoSigmaTrialAll->Draw();
469  cry->cd(3);
470  fHistoChi2TrialAll->Draw();
471  cry->cd(2);
472  fHistoRawYieldTrialAll->Draw();
473  cry->cd(4);
474  fHistoRawYieldDistAll->Draw();
475  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
476  tmean->SetNDC();
477  tmean->Draw();
478  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
479  thrms->SetNDC();
480  thrms->Draw();
481  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
482  tmax->SetNDC();
483  tmax->Draw();
484  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
485  tmin->SetNDC();
486  tmin->Draw();
487  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
488  trms->SetNDC();
489  trms->Draw();
490 
491 }
492 //________________________________________________________________________
493 TH1F* AliHFMultiTrials::RebinHisto(TH1D* hOrig, Int_t reb, Int_t firstUse) const{
494  // Rebin histogram, from bin firstUse to lastUse
495  // Use all bins if firstUse=-1
496 
497  Int_t nBinOrig=hOrig->GetNbinsX();
498  Int_t firstBinOrig=1;
499  Int_t lastBinOrig=nBinOrig;
500  Int_t nBinOrigUsed=nBinOrig;
501  Int_t nBinFinal=nBinOrig/reb;
502  if(firstUse>=1){
503  firstBinOrig=firstUse;
504  nBinFinal=(nBinOrig-firstUse+1)/reb;
505  nBinOrigUsed=nBinFinal*reb;
506  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
507  }else{
508  Int_t exc=nBinOrigUsed%reb;
509  if(exc!=0){
510  nBinOrigUsed-=exc;
511  firstBinOrig+=exc/2;
512  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
513  }
514  }
515 
516  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
517  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
518  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
519  TH1F* hRebin=new TH1F(Form("%s-rebin%d_%d",hOrig->GetName(),reb,firstUse),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
520  Int_t lastSummed=firstBinOrig-1;
521  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
522  Float_t sum=0.;
523  Float_t sum2=0.;
524  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
525  sum+=hOrig->GetBinContent(lastSummed+1);
526  sum2+=hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1);
527  lastSummed++;
528  }
529  hRebin->SetBinContent(iBin,sum);
530  hRebin->SetBinError(iBin,TMath::Sqrt(sum2));
531  }
532  return hRebin;
533 }
534 
535 //________________________________________________________________________
536 void AliHFMultiTrials::BinCount(TH1F* h, TF1* fB, Int_t rebin, Double_t minMass, Double_t maxMass, Double_t& count, Double_t& ecount) const{
537  // compute yield with bin couting
538  Int_t minBinSum=h->FindBin(minMass);
539  Int_t maxBinSum=h->FindBin(maxMass);
540  Double_t cntSig=0.;
541  Double_t cntErr=0.;
542  for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
543  Double_t bkg=fB ? fB->Eval(h->GetBinCenter(iMB))/(Double_t)rebin : 0;
544  cntSig+=(h->GetBinContent(iMB)-bkg);
545  cntErr+=(h->GetBinError(iMB)*h->GetBinError(iMB));
546  }
547  count=cntSig;
548  ecount=TMath::Sqrt(cntErr);
549 }
550 //________________________________________________________________________
552  Int_t iCase){
553  //
554 
555  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
556  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
557  Double_t xc=hCutTmp->GetBinCenter(ib);
558  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
559  hCutTmp->SetBinContent(ib,0.);
560  hCutTmp->SetBinError(ib,0.);
561  }
562  }
563 
564  hCutTmp->Fit("pol2","E0","",hmin,hmax);
565  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
566  TF1* f3=new TF1("myPol3","pol3");
567  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
568  hCutTmp->Fit(f3,"E0","",hmin,hmax);
569  Double_t quickCount=0.;
570  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
571  Double_t xc=hCutTmp->GetBinCenter(ib);
572  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
573  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
574  }
575  }
576  TF1* fSB=new TF1("fSB","[0]*1./(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+[3]+[4]*x+[5]*x*x+[6]*x*x*x");
577  fSB->SetParameter(0,quickCount);
578  fSB->SetParameter(1,fMassD);
579  fSB->SetParameter(2,fSigmaGausMC);
580  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
581  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
582  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
583  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
584  else if(iCase==4){
585  fSB->FixParameter(1,fMassD);
586  fSB->FixParameter(2,fSigmaGausMC);
587  } else if(iCase==5){
588  fSB->FixParameter(1,fMassD);
589  }
590  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
591  // quality cuts
592  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
593  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
594 
595  delete hCutTmp;
596  return kTRUE;
597 }
598 //__________________________________________________________________________________
599 TH1F* AliHFMultiTrials::SetTemplateRefl(const TH1F *h) {
600  fhTemplRefl=(TH1F*)h->Clone("hTemplRefl");
601  return fhTemplRefl;
602 }
603 
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
void SetFixReflOverS(Double_t val, Bool_t fixpar=kTRUE)
TH1F ** fHistoBkgTrial
histo with chi2 from subsamples of trials
double Double_t
Definition: External.C:58
Double_t fMaxYieldGlob
minimum yield
Bool_t fUsePowLawBkg
switch for pol5 background
Definition: External.C:236
TH1F * fhTemplRefl
histo with bin counts from subsamples of trials
TH1F * fHistoSigmaTrialAll
histo with yield from all trials
TNtuple * fNtupleMultiTrials
Bool_t fSaveBkgVal
switch for FixSigFixMean
Double_t fMinYieldGlob
tree
TH1F * SetTemplateRefl(const TH1F *hTemplRefl)
TH1F * fHistoRawYieldDistAll
flag for drawing fits
Bool_t fUseLinBkg
switch for exponential background
TString fileName
void SetInitialGaussianMean(Double_t mean)
Double_t * fLowLimFitSteps
number of steps on the min. mass for fit
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1)
Bool_t fUseFixSigFixMean
switch for FixSigFreeMean
void SetBackHighPolDegree(Int_t deg)
void ConfigurenSigmaBinCSteps(Int_t nSteps, Double_t *values)
TString fSuffix
mass of D meson
TH1F * fHistoMeanTrialAll
histo with gauss sigma from all trials
Int_t types
Bool_t MassFitter(Bool_t draw=kTRUE)
TH1F ** fHistoRawYieldTrial
histo with yield from subsamples of trials
Double_t fMassD
relative variation of the sigma
TF1 * GetBackgroundFullRangeFunc()
Int_t * fRebinSteps
number of rebin steps
TH2F ** fHistoRawYieldTrialBinC
histo with bin counts from subsamples of trials
TH1F ** fHistoSigmaTrial
histo with yield from subsamples of trials
Int_t fNumOfLowLimFitSteps
number of steps in the first bin for rebin
void DrawHistos(TCanvas *cry) const
Double_t GetRawYieldError() const
Int_t typeb
Double_t * sigma
AliHFMassFitterVAR for the fit of invariant mass distribution of charmed mesons.
TH1F ** fHistoSignifTrial
histo with chi2 from subsamples of trials
TH1F ** fHistoRawYieldDist
histo with bin counts from all trials
Double_t * fnSigmaBinCSteps
number of steps on the bin counting
Int_t fFitOption
name to characterize analysis case
int Int_t
Definition: External.C:63
TH1F * fHistoSignifTrialAll
histo with chi2 from all trials
void ConfigureLowLimFitSteps(Int_t nSteps, Double_t *values)
Bool_t fDrawIndividualFits
switch for saving bkg values in nsigma
Double_t GetMeanUncertainty() const
Bool_t fUsePol2Bkg
switch for linear background
Bool_t fUsePol5Bkg
switch for pol4 background
float Float_t
Definition: External.C:68
Double_t GetMean() const
Bool_t fUseFreeS
switch for FixSigDownFreeMean
TH1F ** fHistoRawYieldDistBinC
histo with bkg in mass bin edges from subsamples of trials
Double_t fnSigmaForBkgEval
Double_t GetReducedChiSquare() const
Bool_t fUseFixedMeanFreeS
switch for FreeSigma
void SetReflectionSigmaFactor(Int_t constant)
Definition: External.C:212
Bool_t fUseFixSigUpFreeMean
switch for power law background
Int_t totTrials
TH1F * fHistoRawYieldDistBinCAll
histo with bkg in mass bin edges from all trials
Bool_t DoFitWithPol3Bkg(TH1F *histoToFit, Double_t hmin, Double_t hmax, Int_t theCase)
Double_t GetRawYield() const
TList * fitter
Definition: DrawAnaELoss.C:26
void SetInitialGaussianSigma(Double_t sigma)
change the default value of the mean
TH1F ** fHistoMeanTrial
histo with gauss sigma from subsamples of trials
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH1F * fHistoBkgInBinEdgesTrialAll
histo with bkg from all trials
std::set< std::string > fInvMassFitSaveAsFormats
TH1F * RebinHisto(TH1D *hOrig, Int_t reb, Int_t firstUse) const
Int_t fNumOfRebinSteps
saves the invariant mass fit canvases in the file formats listed in this vector (if empty...
void BinCount(TH1F *h, TF1 *fB, Int_t rebin, Double_t minMass, Double_t maxMass, Double_t &count, Double_t &ecount) const
TH2F * fHistoRawYieldTrialBinCAll
histo with bin counts from all trials
Double_t fSigmaMCVariation
sigma of D meson peak from MC
TH1F * SetTemplateReflections(const TH1 *h, TString option="templ", Double_t minRange=1.72, Double_t maxRange=2.05)
Bool_t fUsePol3Bkg
switch for pol2 background
void SetFixGaussianSigma(Double_t sigma=0.012, Bool_t fixpar=kTRUE)
Bool_t fUseFixSigDownFreeMean
switch for FixSigUpFreeMean
Bool_t fUsePol4Bkg
switch for pol3 background
TH1F * fHistoBkgTrialAll
histo with chi2 from all trials
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t minMass
Int_t rebin
TH1F ** fHistoChi2Trial
histo with gauss mean from subsamples of trials
TH1F * fHistoRawYieldTrialAll
histo with yield from all trials
void ConfigureRebinSteps(Int_t nSteps, Int_t *values)
void ConfigureUpLimFitSteps(Int_t nSteps, Double_t *values)
void SaveToRoot(TString fileName, TString option="recreate") const
Bool_t fUseFixSigFreeMean
switch for FixedMeanFreeS
TH1F ** fHistoBkgInBinEdgesTrial
histo with bkg from subsamples of trials
Double_t maxMass
Double_t GetSigma() const
bool Bool_t
Definition: External.C:53
Bool_t fUseExpoBkg
LL or chi2 fit.
Bool_t DoMultiTrials(TH1D *hInvMassHisto, TPad *thePad=0x0)
void SetFixGaussianMean(Double_t mean=1.865, Bool_t fixpar=kTRUE)
Double_t * fUpLimFitSteps
number of steps on the max. mass for fit
Definition: External.C:196
Float_t fFixRefloS
template of reflection contribution
Double_t minMassForFit
TH1F * fHistoChi2TrialAll
histo with gauss mean from all trials
Double_t maxMassForFit
Double_t GetSigmaUncertainty() const
Bool_t fUsePowLawTimesExpoBkg
switch for power law background