AliPhysics  97344c9 (97344c9)
 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  }
322  }
323  // else{
324  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
325  // if(out && thePad){
326  // thePad->Clear();
327  // hRebinned->Draw();
328  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
329  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
330  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
331  // fB1->SetLineColor(2);
332  // fB1->Draw("same");
333  // fSB->SetLineColor(4);
334  // fSB->Draw("same");
335  // thePad->Update();
336  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
337  // sigma=fSB->GetParameter(2);
338  // esigma=fSB->GetParError(2);
339  // if(esigma<0.00001) esigma=0.0001;
340  // pos=fSB->GetParameter(1);
341  // epos=fSB->GetParError(1);
342  // if(epos<0.00001) epos=0.0001;
343  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
344  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
345  // }
346  // }
347  xnt[7]=chisq;
348  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
349  xnt[8]=significance;
350  xnt[9]=pos;
351  xnt[10]=epos;
352  xnt[11]=sigma;
353  xnt[12]=esigma;
354  xnt[13]=ry;
355  xnt[14]=ery;
356  fHistoRawYieldDistAll->Fill(ry);
357  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
358  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
359  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
360  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
361  fHistoMeanTrialAll->SetBinContent(globBin,pos);
362  fHistoMeanTrialAll->SetBinError(globBin,epos);
363  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
364  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
365  fHistoSignifTrialAll->SetBinContent(globBin,significance);
366  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
367  if(fSaveBkgVal) {
368  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
369  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
370  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
371  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
372  }
373 
374  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
375  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
376  fHistoRawYieldDist[theCase]->Fill(ry);
377  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
378  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
379  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
380  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
381  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
382  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
383  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
384  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
385  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
386  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
387  if(fSaveBkgVal) {
388  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
389  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
390  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
391  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
392  }
393 
394  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
395  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
396  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
397  if(minMassBC>minMassForFit &&
398  maxMassBC<maxMassForFit &&
399  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
400  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
401  Double_t cnts,ecnts;
402  BinCount(hRebinned,fB1,1,minMassBC,maxMassBC,cnts,ecnts);
403  ++itrialBC;
404  fHistoRawYieldDistBinCAll->Fill(cnts);
405  fHistoRawYieldTrialBinCAll->SetBinContent(globBin,iStepBC+1,cnts);
406  fHistoRawYieldTrialBinCAll->SetBinError(globBin,iStepBC+1,ecnts);
407  fHistoRawYieldTrialBinC[theCase]->SetBinContent(itrial,iStepBC+1,cnts);
408  fHistoRawYieldTrialBinC[theCase]->SetBinError(itrial,iStepBC+1,ecnts);
409  fHistoRawYieldDistBinC[theCase]->Fill(cnts);
410  }
411  }
412  }
413  delete fitter;
414  // if(typeb>4) delete fB1;
415  fNtupleMultiTrials->Fill(xnt);
416  }
417  }
418  }
419  }
420  delete hRebinned;
421  }
422  }
423  return kTRUE;
424 }
425 
426 //________________________________________________________________________
428  // save histos in a root file for further analysis
429  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
430  TFile* outHistos=new TFile(fileName.Data(),option.Data());
431  fHistoRawYieldTrialAll->Write();
432  fHistoSigmaTrialAll->Write();
433  fHistoMeanTrialAll->Write();
434  fHistoChi2TrialAll->Write();
435  fHistoSignifTrialAll->Write();
436  if(fSaveBkgVal) {
437  fHistoBkgTrialAll->Write();
438  fHistoBkgInBinEdgesTrialAll->Write();
439  }
440  fHistoRawYieldDistBinCAll->Write();
441  fHistoRawYieldTrialBinCAll->Write();
442  for(Int_t ic=0; ic<nCases; ic++){
443  fHistoRawYieldTrial[ic]->Write();
444  fHistoSigmaTrial[ic]->Write();
445  fHistoMeanTrial[ic]->Write();
446  fHistoChi2Trial[ic]->Write();
447  fHistoSignifTrial[ic]->Write();
448  if(fSaveBkgVal) {
449  fHistoBkgTrial[ic]->Write();
450  fHistoBkgInBinEdgesTrial[ic]->Write();
451  }
452  fHistoRawYieldTrialBinC[ic]->Write();
453  fHistoRawYieldDistBinC[ic]->Write();
454  }
455  fNtupleMultiTrials->Write();
456  outHistos->Close();
457  delete outHistos;
458 }
459 
460 //________________________________________________________________________
461 void AliHFMultiTrials::DrawHistos(TCanvas* cry) const{
462  // draw histos
463  cry->Divide(2,2);
464  cry->cd(1);
465  fHistoSigmaTrialAll->Draw();
466  cry->cd(3);
467  fHistoChi2TrialAll->Draw();
468  cry->cd(2);
469  fHistoRawYieldTrialAll->Draw();
470  cry->cd(4);
471  fHistoRawYieldDistAll->Draw();
472  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
473  tmean->SetNDC();
474  tmean->Draw();
475  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
476  thrms->SetNDC();
477  thrms->Draw();
478  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
479  tmax->SetNDC();
480  tmax->Draw();
481  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
482  tmin->SetNDC();
483  tmin->Draw();
484  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
485  trms->SetNDC();
486  trms->Draw();
487 
488 }
489 //________________________________________________________________________
490 TH1F* AliHFMultiTrials::RebinHisto(TH1D* hOrig, Int_t reb, Int_t firstUse) const{
491  // Rebin histogram, from bin firstUse to lastUse
492  // Use all bins if firstUse=-1
493 
494  Int_t nBinOrig=hOrig->GetNbinsX();
495  Int_t firstBinOrig=1;
496  Int_t lastBinOrig=nBinOrig;
497  Int_t nBinOrigUsed=nBinOrig;
498  Int_t nBinFinal=nBinOrig/reb;
499  if(firstUse>=1){
500  firstBinOrig=firstUse;
501  nBinFinal=(nBinOrig-firstUse+1)/reb;
502  nBinOrigUsed=nBinFinal*reb;
503  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
504  }else{
505  Int_t exc=nBinOrigUsed%reb;
506  if(exc!=0){
507  nBinOrigUsed-=exc;
508  firstBinOrig+=exc/2;
509  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
510  }
511  }
512 
513  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
514  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
515  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
516  TH1F* hRebin=new TH1F(Form("%s-rebin%d_%d",hOrig->GetName(),reb,firstUse),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
517  Int_t lastSummed=firstBinOrig-1;
518  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
519  Float_t sum=0.;
520  Float_t sum2=0.;
521  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
522  sum+=hOrig->GetBinContent(lastSummed+1);
523  sum2+=hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1);
524  lastSummed++;
525  }
526  hRebin->SetBinContent(iBin,sum);
527  hRebin->SetBinError(iBin,TMath::Sqrt(sum2));
528  }
529  return hRebin;
530 }
531 
532 //________________________________________________________________________
533 void AliHFMultiTrials::BinCount(TH1F* h, TF1* fB, Int_t rebin, Double_t minMass, Double_t maxMass, Double_t& count, Double_t& ecount) const{
534  // compute yield with bin couting
535  Int_t minBinSum=h->FindBin(minMass);
536  Int_t maxBinSum=h->FindBin(maxMass);
537  Double_t cntSig=0.;
538  Double_t cntErr=0.;
539  for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
540  Double_t bkg=fB ? fB->Eval(h->GetBinCenter(iMB))/(Double_t)rebin : 0;
541  cntSig+=(h->GetBinContent(iMB)-bkg);
542  cntErr+=(h->GetBinError(iMB)*h->GetBinError(iMB));
543  }
544  count=cntSig;
545  ecount=TMath::Sqrt(cntErr);
546 }
547 //________________________________________________________________________
549  Int_t iCase){
550  //
551 
552  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
553  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
554  Double_t xc=hCutTmp->GetBinCenter(ib);
555  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
556  hCutTmp->SetBinContent(ib,0.);
557  hCutTmp->SetBinError(ib,0.);
558  }
559  }
560 
561  hCutTmp->Fit("pol2","E0","",hmin,hmax);
562  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
563  TF1* f3=new TF1("myPol3","pol3");
564  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
565  hCutTmp->Fit(f3,"E0","",hmin,hmax);
566  Double_t quickCount=0.;
567  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
568  Double_t xc=hCutTmp->GetBinCenter(ib);
569  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
570  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
571  }
572  }
573  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");
574  fSB->SetParameter(0,quickCount);
575  fSB->SetParameter(1,fMassD);
576  fSB->SetParameter(2,fSigmaGausMC);
577  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
578  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
579  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
580  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
581  else if(iCase==4){
582  fSB->FixParameter(1,fMassD);
583  fSB->FixParameter(2,fSigmaGausMC);
584  } else if(iCase==5){
585  fSB->FixParameter(1,fMassD);
586  }
587  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
588  // quality cuts
589  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
590  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
591 
592  delete hCutTmp;
593  return kTRUE;
594 }
595 //__________________________________________________________________________________
596 TH1F* AliHFMultiTrials::SetTemplateRefl(const TH1F *h) {
597  fhTemplRefl=(TH1F*)h->Clone("hTemplRefl");
598  return fhTemplRefl;
599 }
600 
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
TH1F * RebinHisto(TH1D *hOrig, Int_t reb, Int_t firstUse) const
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