]>
Commit | Line | Data |
---|---|---|
1282b9e5 | 1 | TCanvas *cCompareSimple,*cCompareCum; |
2 | TCanvas *cSubtraction,*cCompareSubtractSimple,*cCompareSubtractCum; | |
3 | const Int_t timeWait=500; | |
4 | ||
5 | void AnalyzeCharmFractionHists(Int_t ptbin,Int_t rebin){ | |
6 | ||
7 | // | |
8 | // Macro for analyzing the impact parameter distribution of the D0 candidates. | |
9 | // | |
10 | // andrea.rossi@ts.infn.it | |
11 | // | |
12 | //========================================================================== | |
13 | ||
14 | ||
15 | DrawHistos("d0D0_PureBack.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin); | |
16 | SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin); | |
17 | return; | |
18 | ||
19 | } | |
20 | ||
21 | TH1F* CumulativeHist(const TH1F *h1,TString nameHist=0x0){ | |
22 | //ASSUME SAME BINNING | |
23 | TString histname="hCumulative"; | |
24 | if(!nameHist.IsNull())histname.Append(nameHist.Data()); | |
25 | Int_t nbins; | |
26 | Double_t minX,maxX,binwidth,cumul=0.,errcumul=0.; | |
27 | nbins=h1->GetNbinsX(); | |
28 | minX=h1->GetBinCenter(1); | |
29 | maxX=h1->GetBinCenter(nbins); | |
30 | binwidth=h1->GetBinWidth(1); | |
31 | ||
32 | TH1F *h1copy=new TH1F(); | |
33 | *h1copy=*h1; | |
34 | //h1copy->Sumw2(); | |
35 | ||
36 | TH1F *hCumul=new TH1F(histname.Data(),histname.Data(),nbins,minX-binwidth/2.,maxX+binwidth/2.); | |
37 | ||
38 | ||
39 | for(Int_t j=1;j<=nbins;j++){ | |
40 | cumul+=h1copy->GetBinContent(j); | |
41 | hCumul->SetBinContent(j,cumul); | |
42 | } | |
43 | hCumul->Sumw2(); | |
44 | hCumul->Scale(1./h1copy->Integral(1,nbins)); | |
45 | ||
46 | delete h1copy; | |
47 | return hCumul; | |
48 | } | |
49 | ||
50 | ||
51 | ///////////////////// HISTOGRAMS COMPARISON ////////////////////////////////// | |
52 | // | |
53 | // Compare the desired histograms. The default case is the comparison | |
54 | // between the d0 distribution for background (from MC) D0 candidates with the | |
55 | // d0 distribution for "side bands" d0 candidates | |
56 | // | |
57 | /////////////////////////////////////////////////////////////////////////////////// | |
58 | ||
59 | void DrawHistos(Int_t pt,Int_t rebin=1){ | |
60 | TString hist="hd0D0"; | |
61 | DrawHistos("","",hist,pt,rebin); | |
62 | return; | |
63 | } | |
64 | void DrawHistos(TString file1,TString file2,TString hist,Int_t pt,Int_t rebin=1){ | |
65 | ||
66 | if(pt>=0){ | |
67 | hist.Append("pt"); | |
68 | hist+=pt; | |
69 | } | |
70 | TString nameH=hist; | |
71 | Double_t integr1,integr2; | |
72 | if(file1.IsNull())file1="d0D0_PureBack.root"; | |
73 | if(file2.IsNull())file2="d0D0NoMCSel_SideBand.root"; | |
74 | TFile *fSide=TFile::Open(file2.Data()); | |
75 | TFile *fmer=TFile::Open(file1.Data()); | |
76 | TH1F *h1=(TH1F*)fmer->Get(hist.Data()); | |
77 | nameH.Append("_Gaus"); | |
78 | h1->SetName(nameH.Data()); | |
79 | nameH=hist; | |
80 | TH1F *h2=(TH1F*)fSide->Get(hist.Data()); | |
81 | nameH.Append("_SideBands"); | |
82 | h2->SetName(nameH.Data()); | |
83 | ||
84 | Bool_t satisfied=kFALSE; | |
85 | ||
86 | while(!satisfied){ | |
87 | h1->Rebin(rebin); | |
88 | h2->Rebin(rebin); | |
89 | TH1F *hCumul1=CumulativeHist(h1,"_Gaus"); | |
90 | hCumul1->SetDrawOption("E4"); | |
91 | TH1F *hCumul2=CumulativeHist(h2,"_SideBands"); | |
92 | hCumul2->SetDrawOption("E4"); | |
93 | ||
94 | ||
95 | TH1F *hDivCumul=new TH1F(); | |
96 | *hDivCumul=*hCumul2; | |
97 | // hDivCumul->Sumw2(); | |
98 | hDivCumul->Divide(hCumul1,hCumul2); | |
99 | hDivCumul->SetName("RatioCumulGausSdBands"); | |
100 | hDivCumul->SetTitle("Ratio Cumulative Gaus Over SideBands"); | |
101 | ||
102 | TH1F *hGausComp=new TH1F(); | |
103 | *hGausComp=*h1; | |
104 | nameH=h1->GetName(); | |
105 | nameH.Append("Compare"); | |
106 | hGausComp->SetName(nameH.Data()); | |
107 | hGausComp->SetLineColor(kRed); | |
108 | hGausComp->Sumw2(); | |
109 | integr1=hGausComp->Integral(); | |
110 | ||
111 | TH1F *hSideComp=new TH1F(); | |
112 | *hSideComp=*h2; | |
113 | nameH=h2->GetName(); | |
114 | nameH.Append("Compare"); | |
115 | hSideComp->SetName(nameH.Data()); | |
116 | hSideComp->SetLineColor(kBlue); | |
117 | hSideComp->Sumw2(); | |
118 | integr2=hSideComp->Integral(); | |
119 | if(integr2>integr1)hSideComp->Scale(integr1/integr2); | |
120 | else hGausComp->Scale(integr2/integr1); | |
121 | ||
122 | TH1F *hDivision=new TH1F(); | |
123 | *hDivision=*h1; | |
124 | hDivision->SetName("RatioGausSdBands"); | |
125 | hDivision->SetTitle("Ratio Gaus Over Sd Bands"); | |
126 | hDivision->Sumw2(); | |
127 | hDivision->Divide(hGausComp,hSideComp); | |
128 | ||
129 | ||
130 | cCompareSimple=new TCanvas("cCompareSimple","Comparison of Under Gaus and Under Side Band background",700,700); | |
131 | cCompareSimple->Divide(1,2); | |
132 | cCompareSimple->cd(1); | |
133 | hSideComp->Draw(); | |
134 | hGausComp->Draw("Sames"); | |
135 | cCompareSimple->Update(); | |
136 | TPaveStats *p=hGausComp->FindObject("stats"); | |
137 | p->SetTextColor(kRed); | |
138 | p=(TPaveStats*)hSideComp->FindObject("stats"); | |
139 | p->SetTextColor(kBlue); | |
140 | cCompareSimple->cd(2); | |
141 | hDivision->Draw(); | |
142 | cCompareSimple->Update(); | |
143 | ||
144 | cCompareCum=new TCanvas("cCompareCum","Comparison of cumulative function",700,700); | |
145 | ||
146 | cCompareCum->Divide(1,2); | |
147 | ||
148 | cCompareCum->cd(1); | |
149 | ||
150 | hCumul1->SetLineColor(kRed); | |
151 | hCumul2->SetLineColor(kBlue); | |
152 | hCumul1->SetFillColor(kRed); | |
153 | hCumul2->SetFillColor(kBlue); | |
154 | hCumul1->Draw(); | |
155 | hCumul2->Draw("Sames"); | |
156 | cCompareCum->cd(2); | |
157 | hDivCumul->Draw(); | |
158 | cCompareCum->Update(); | |
159 | p=(TPaveStats*)hCumul1->FindObject("stats"); | |
160 | p->SetTextColor(kRed); | |
161 | p=(TPaveStats*)hCumul2->FindObject("stats"); | |
162 | p->SetTextColor(kBlue); | |
163 | cCompareCum->Update(); | |
164 | // cCompareCum->cd(2); | |
165 | //hCumul2->Draw(); | |
166 | // cCompareCum->Update(); | |
167 | ||
168 | for(Int_t timewait=0;timewait<timeWait;timewait++){ | |
169 | cCompareCum->Update(); | |
170 | cCompareSimple->Update(); | |
171 | } | |
172 | cout<<"Are you satisfied?"<<endl; | |
173 | cin>>satisfied; | |
174 | ||
175 | if(!satisfied){ | |
176 | ||
177 | ||
178 | delete hCumul1; | |
179 | delete hCumul2; | |
180 | delete hDivCumul; | |
181 | delete cCompareCum; | |
182 | delete hSideComp; | |
183 | delete hGausComp; | |
184 | delete cCompareSimple; | |
185 | delete hDivision; | |
186 | ||
187 | cout<<"Rebin"<<endl; | |
188 | cin>>rebin; | |
189 | if(rebin<0)return; | |
190 | } | |
191 | } | |
192 | ||
193 | return; | |
194 | ||
195 | } | |
196 | ||
197 | /////////////////////////////////////// BACKGROUND SUBTRACTION EXAMPLE /////////////////// | |
198 | // | |
199 | // | |
200 | // Performs: Signal Imp.Par Distr - B* Normalized Background Imp.Par Distr.(from Side Bands) | |
201 | // and compare it with the MC Signal Imp. Par distribution. | |
202 | // The amount of Background B is taken as the true one from the difference between | |
203 | // the integrals of the "Observed Signal" and the MC Signal d0 distributions. | |
204 | // | |
205 | ////////////////////////////////////////////////////////////////////////////////////////////// | |
206 | ||
207 | void SubtractHist(Int_t pt,Int_t rebin=1){ | |
208 | SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",pt,rebin); | |
209 | ||
210 | return; | |
211 | } | |
212 | ||
213 | void SubtractHist(TString fileSignal,TString fileNoMC,TString fileNoMCSB,TString hist,Int_t pt,Int_t rebin){ | |
214 | ||
215 | if(pt>=0){ | |
216 | hist.Append("pt"); | |
217 | hist+=pt; | |
218 | } | |
219 | ||
220 | TString nameH=hist; | |
221 | Double_t integr1,integr2; | |
222 | /* if(fileNoMC.IsNull())fileNoMC="d0D0merged.root"; | |
223 | if(fileNoMCSB.IsNull())fileNoMCSB="d0D0SideBmerged.root";*/ | |
224 | if(gSystem->AccessPathName(fileSignal.Data())){ | |
225 | Printf("Wrong signal file! \n"); | |
226 | return; | |
227 | } | |
228 | if(gSystem->AccessPathName(fileNoMC.Data())){ | |
229 | Printf("Wrong d0distr under inv mass peak file! \n"); | |
230 | return; | |
231 | } | |
232 | if(gSystem->AccessPathName(fileNoMCSB.Data())){ | |
233 | Printf("Wrong d0distr Side band file! \n"); | |
234 | return; | |
235 | } | |
236 | ||
237 | TFile *fSide=TFile::Open(fileNoMCSB.Data()); | |
238 | TFile *fNoMCSignal=TFile::Open(fileNoMC.Data()); | |
239 | ||
240 | TH1F *h1=(TH1F*)fNoMCSignal->Get(hist.Data()); | |
241 | nameH.Append("_SelSignal"); | |
242 | h1->SetName(nameH.Data()); | |
243 | nameH=hist; | |
244 | h1->Sumw2(); | |
245 | ||
246 | ||
247 | TH1F *h2=(TH1F*)fSide->Get(hist.Data()); | |
248 | nameH.Append("_SideBands"); | |
249 | h2->SetName(nameH.Data()); | |
250 | h2->Sumw2(); | |
251 | ||
252 | ||
253 | Double_t integrGaus,integrSB,integrSign; | |
254 | integrGaus=h1->Integral(); | |
255 | integrSB=h2->Integral(); | |
256 | ||
257 | ||
258 | TFile *fSignal=TFile::Open(fileSignal.Data()); | |
259 | TH1F *hSign=(TH1F*)fSignal->Get(hist.Data()); | |
260 | nameH=hist; | |
261 | nameH.Append("_MCSignal"); | |
262 | hSign->SetName(nameH.Data()); | |
263 | hSign->Sumw2(); | |
264 | ||
265 | Double_t integrGaus,integrSB,integrSign; | |
266 | integrGaus=h1->Integral(); | |
267 | integrSB=h2->Integral(); | |
268 | integrSign=hSign->Integral(); | |
269 | ||
270 | ||
271 | Bool_t satisfied=kFALSE; | |
272 | ||
273 | while(!satisfied){ | |
274 | ||
275 | h1->Rebin(rebin); | |
276 | h2->Rebin(rebin); | |
277 | hSign->Rebin(rebin); | |
278 | ||
279 | TH1F *hGausSubtract=new TH1F(); | |
280 | *hGausSubtract=*h1; | |
281 | nameH=hist; | |
282 | hist.Append("_Subtracted"); | |
283 | hGausSubtract->SetName(hist.Data()); | |
284 | //hGausSubtract->Sumw2(); | |
285 | hGausSubtract->Add(h1,h2,1.,-1./integrSB*(integrGaus-integrSign)); | |
286 | /* hGausSubtract->Draw(); | |
287 | return;*/ | |
288 | cSubtraction=new TCanvas("cSubtraction","cSubtraction",600,700); | |
289 | cSubtraction->cd(); | |
290 | hGausSubtract->SetLineColor(kRed); | |
291 | hGausSubtract->Draw("E0"); | |
292 | hSign->SetLineColor(kBlue); | |
293 | hSign->Draw("sames"); | |
294 | cSubtraction->Update(); | |
295 | TPaveStats *p=hGausSubtract->FindObject("stats"); | |
296 | p->SetTextColor(kRed); | |
297 | cSubtraction->Update(); | |
298 | p=(TPaveStats*)hSign->FindObject("stats"); | |
299 | p->SetTextColor(kBlue); | |
300 | cSubtraction->Update(); | |
301 | /*TString strText="Side Band integral: "; | |
302 | strText+=integrSB; | |
303 | TLatex *lat=new TLatex(500.,500.,strText.Data());// *text=p->AddText(strText.Data()); | |
304 | lat->Draw(); | |
305 | TPaveText *pvInfo = new TPaveText(72.74276,79.18782,92.84497,95.43147,"brNDC"); | |
306 | pvInfo->SetName("pvInfo"); | |
307 | pvInfo->SetBorderSize(2); | |
308 | pvInfo->SetFillColor(19); | |
309 | pvInfo->SetTextAlign(12); | |
310 | TText *text=p->AddText(strText.Data()); | |
311 | pvInfo->Draw(); | |
312 | hGausSubtract->GetListOfFunctions()->Add(pvInfo); | |
313 | // pvInfo->SetParent(hGausSubtract->GetListOfFunctions()); | |
314 | ||
315 | // strText="MCSignal integral: "; | |
316 | //strText+=integrSign; | |
317 | text=p->AddText(strText.Data()); | |
318 | strText="Sel Signal integral: "; | |
319 | strText+=integrGaus; | |
320 | text=p->AddText(strText.Data()); | |
321 | cSubtraction->Modified();*/ | |
322 | ||
323 | ||
324 | TH1F *hCumulGausSubtract=CumulativeHist(hGausSubtract,"_BackSubtr"); | |
325 | hCumulGausSubtract->SetDrawOption("E4"); | |
326 | TH1F *hCumulSign=CumulativeHist(hSign,"_MCSignal"); | |
327 | hCumulSign->SetDrawOption("E4"); | |
328 | ||
329 | ||
330 | TH1F *hDivCumul=new TH1F(); | |
331 | *hDivCumul=*hCumulSign; | |
332 | hDivCumul->Divide(hCumulGausSubtract,hCumulSign); | |
333 | hDivCumul->SetName("RatioCumulBackSubtr_MCSignal"); | |
334 | hDivCumul->SetTitle("Ratio Cumulative BackSubtracted Over MCSignal"); | |
335 | ||
336 | TH1F *hGausSubComp=new TH1F(); | |
337 | *hGausSubComp=*hGausSubtract; | |
338 | nameH=hGausSubtract->GetName(); | |
339 | nameH.Append("Compare"); | |
340 | hGausSubComp->SetName(nameH.Data()); | |
341 | hGausSubComp->SetLineColor(kRed); | |
342 | integr1=hGausSubComp->Integral(); | |
343 | ||
344 | TH1F *hSignComp=new TH1F(); | |
345 | *hSignComp=*hSign; | |
346 | nameH=hSign->GetName(); | |
347 | nameH.Append("Compare"); | |
348 | hSignComp->SetName(nameH.Data()); | |
349 | hSignComp->SetLineColor(kBlue); | |
350 | integr2=hSignComp->Integral(); | |
351 | if(integr2>integr1)hSignComp->Scale(integr1/integr2); | |
352 | else hGausSubComp->Scale(integr2/integr1); | |
353 | ||
354 | TH1F *hDivision=new TH1F(); | |
355 | *hDivision=*hGausSubtract; | |
356 | hDivision->SetName("RatioBackSubtr_Signal"); | |
357 | hDivision->SetTitle("Ratio BackSubtracted Over MCSignal"); | |
358 | hDivision->Divide(hGausSubComp,hSignComp); | |
359 | ||
360 | ||
361 | cCompareSubtractSimple=new TCanvas("cCompareSubtractSimple","Comparison of BackSubtracted and MCSignal",700,700); | |
362 | cCompareSubtractSimple->Divide(1,2); | |
363 | cCompareSubtractSimple->cd(1); | |
364 | hSignComp->Draw(); | |
365 | hGausSubComp->Draw("Sames"); | |
366 | cCompareSubtractSimple->cd(2); | |
367 | hDivision->SetLineColor(1); | |
368 | hDivision->Draw(); | |
369 | cCompareSubtractSimple->Update(); | |
370 | p=(TPaveStats*)hGausSubComp->FindObject("stats"); | |
371 | p->SetTextColor(kRed); | |
372 | p=(TPaveStats*)hSignComp->FindObject("stats"); | |
373 | p->SetTextColor(kBlue); | |
374 | cCompareSubtractSimple->Update(); | |
375 | ||
376 | ||
377 | cCompareSubtractCum=new TCanvas("cCompareSubtractCumulative","Comparison of cumulative functions",700,700); | |
378 | ||
379 | cCompareSubtractCum->Divide(1,2); | |
380 | ||
381 | cCompareSubtractCum->cd(1); | |
382 | ||
383 | hCumulGausSubtract->SetLineColor(kRed); | |
384 | hCumulSign->SetLineColor(kBlue); | |
385 | hCumulGausSubtract->SetFillColor(kRed); | |
386 | hCumulSign->SetFillColor(kBlue); | |
387 | hCumulGausSubtract->Draw(); | |
388 | hCumulSign->Draw("Sames"); | |
389 | cCompareSubtractCum->cd(2); | |
390 | hDivCumul->SetLineColor(1); | |
391 | hDivCumul->Draw(); | |
392 | ||
393 | cCompareSubtractCum->Update(); | |
394 | p=(TPaveStats*)hCumulGausSubtract->FindObject("stats"); | |
395 | p->SetTextColor(kRed); | |
396 | p=(TPaveStats*)hCumulSign->FindObject("stats"); | |
397 | p->SetTextColor(kBlue); | |
398 | cCompareSubtractCum->Update(); | |
399 | ||
400 | ||
401 | for(Int_t timewait=0;timewait<timeWait;timewait++){ | |
402 | cCompareSubtractCum->Update(); | |
403 | cCompareSubtractSimple->Update(); | |
404 | cSubtraction->Update(); | |
405 | } | |
406 | ||
407 | cout<<"Are you satisfied?"<<endl; | |
408 | cin>>satisfied; | |
409 | ||
410 | if(!satisfied){ | |
411 | cout<<"Rebin"<<endl; | |
412 | cin>>rebin; | |
413 | delete hCumulGausSubtract; | |
414 | delete hCumulSign; | |
415 | delete hDivCumul; | |
416 | delete cCompareSubtractCum; | |
417 | delete hSignComp; | |
418 | delete hGausSubComp; | |
419 | delete cCompareSubtractSimple; | |
420 | delete hDivision; | |
421 | } | |
422 | } | |
423 | ||
424 | printf("Side Band integral: %d \n",integrSB); | |
425 | printf("Selected Signal integral: %d \n",integrGaus); | |
426 | printf("MC Signal integral: %d \n",integrSign); | |
427 | printf(" -> Background = %d \n",integrGaus-integrSign); | |
428 | ||
429 | return; | |
430 | ||
431 | } | |
432 | ||
433 | //////////////////////////////// FIT DISTRIBUTION WITH GAUSSIAN + EXP TAILS /////////////////// | |
434 | // | |
435 | // Fit a istogram with a gaus(mean,sigma) + a*exp(mean2,lambda) function | |
436 | // | |
437 | ///////////////////////////////////////////////////////////////////////////////////////////////// | |
438 | void DoFit(TString filename,TString histo,TString gausSide,Int_t ptbin,Int_t rebin=-1){ | |
439 | ||
440 | TString histname=histo; | |
441 | TString fileout=histo; | |
442 | ||
443 | if(ptbin>=0){ | |
444 | histname.Append("pt"); | |
445 | histname+=ptbin; | |
446 | fileout.Append("pt"); | |
447 | fileout+=ptbin; | |
448 | } | |
449 | else fileout.Append("ptAll"); | |
450 | ||
451 | TFile *fIN=TFile::Open(filename.Data()); | |
452 | ||
453 | TH1F *h=(TH1F*)fIN->Get(histname.Data()); | |
454 | ||
455 | Int_t ok=0; | |
456 | Int_t binL=1,binR=h->GetNbinsX(); | |
457 | ||
458 | while(ok!=1){ | |
459 | if(rebin>0)h->Rebin(rebin); | |
460 | ||
461 | Double_t integr=h->Integral(binL,binR); | |
462 | Double_t binwidth=h->GetBinWidth(10); | |
463 | Double_t minX,maxX; | |
464 | Int_t nbin=h->GetNbinsX(); | |
465 | ||
466 | ||
467 | TCanvas *c1=new TCanvas("c1","c1",700,600); | |
468 | c1->cd(); | |
469 | // h->SetFillColor(kRed); | |
470 | h->Draw(""); | |
471 | // c1->SetLogy(); | |
472 | c1->Update(); | |
473 | ||
474 | TF1 *f=CreateFunc(integr,binwidth); | |
475 | h->Fit(f->GetName(),"RI","",h->GetBinCenter(binL)-binwidth/2.,h->GetBinCenter(binR)+binwidth/2.); | |
476 | c1->Update(); | |
477 | ||
478 | cout<<"Is it ok?"<<endl; | |
479 | cin>>ok; | |
480 | if(ok==1){ | |
481 | ||
482 | fileout.Append("_"); | |
483 | fileout.Append(gausSide.Data()); | |
484 | fileout.Append("Fit.root"); | |
485 | ||
486 | TFile *fOUT=new TFile(fileout.Data(),"RECREATE"); | |
487 | fOUT->cd(); | |
488 | c1->Write(); | |
489 | h->Write(); | |
490 | fOUT->Close(); | |
491 | delete c1; | |
492 | delete f; | |
493 | } | |
494 | else if(ok==0){ | |
495 | cout<<"rebin value= "; | |
496 | cin>>rebin; | |
497 | cout<<endl; | |
498 | cout<<"Change Interval?"<<endl; | |
499 | cin>>ok; | |
500 | if(ok){ | |
501 | cout<<"Lower value= "; | |
502 | cin>>minX; | |
503 | cout<<endl; | |
504 | cout<<"Upper value= "; | |
505 | cin>>maxX; | |
506 | cout<<endl; | |
507 | ||
508 | binL=h->FindBin(minX); | |
509 | binR=h->FindBin(maxX); | |
510 | } | |
511 | ok=0; | |
512 | delete f; | |
513 | delete c1; | |
514 | } | |
515 | else { | |
516 | ok=1; | |
517 | delete f; | |
518 | delete c1; | |
519 | } | |
520 | } | |
521 | ||
522 | return; | |
523 | } | |
524 | ||
525 | ||
526 | TF1* CreateFunc(Double_t integral,Double_t binw){ | |
527 | TF1 *func=new TF1("func","[5]*[6]*((1.-[2])*1./2./[1]*TMath::Exp(-TMath::Abs((x-[0])/[1]))+[2]/TMath::Sqrt(2*TMath::Pi())/[4]*TMath::Exp(-1*(x-[3])*(x-[3])/2./[4]/[4]))"); | |
528 | ||
529 | func->FixParameter(5,integral); | |
530 | func->SetParName(5,"histInt"); | |
531 | func->FixParameter(6,binw); | |
532 | func->SetParName(6,"binW"); | |
533 | func->SetParLimits(0,-100,100.); | |
534 | func->SetParName(0,"expMean"); | |
535 | func->SetParLimits(1,0.001,1000); | |
536 | func->SetParName(1,"expDecay"); | |
537 | func->SetParLimits(2,0.00001,1.); | |
538 | func->SetParName(2,"fracGaus"); | |
539 | func->SetParLimits(3,-100,100.); | |
540 | func->SetParName(3,"gausMean"); | |
541 | func->SetParLimits(4,5.,200.); | |
542 | func->SetParName(4,"gausSigma"); | |
543 | ||
544 | func->SetParameter(1,50.); | |
545 | func->SetParameter(4,50.); | |
546 | ||
547 | return func; | |
548 | ||
549 | } | |
550 | ||
551 | ////////////////////////////////// CHI2 TEST: //////////////////////// | |
552 | // IN PROGRESS | |
553 | // | |
554 | /////////////////////////////////////////////////////////////////////////////// | |
555 | Double_t GetChi2(const TH1F *h1,const TH1F *h2,TH1F *hchi2,Int_t &ndof){//ASSUME SAME BINNING | |
556 | ||
557 | Int_t int1=h1->Integral(); | |
558 | Int_t int2=h2->Integral(); | |
559 | Int_t nm,nM; | |
560 | TH1F **h=new TH1F*[2]; | |
561 | TH1F *hchisq=new TH1F("hChi2","Chi2 histo",1000,0.,10.); | |
562 | ||
563 | if(int2>int1){ | |
564 | nM=int2; | |
565 | nm=int1; | |
566 | h[0]=h2; | |
567 | h[1]=h1; | |
568 | } | |
569 | else { | |
570 | nM=int1; | |
571 | nm=int2; | |
572 | h[0]=h1; | |
573 | h[1]=h2; | |
574 | } | |
575 | ndof=0; | |
576 | Int_t nbin=h[0]->GetNbinsX();//ASSUME SAME BINNING | |
577 | ||
578 | Double_t chi2=0.,chi2Sum=0.; | |
579 | Double_t f1,fTrue; | |
580 | for(Int_t j=1;j<=nbin;j++){//THE BINNING IS NOT TAKEN INTO ACCOUNT | |
581 | ||
582 | if(h[1]->GetBinContent(j)==0||h[0]->GetBinContent(j)==0)continue; | |
583 | fTrue=h[0]->GetBinContent(j)/nM; | |
584 | chi2=(h[1]->GetBinContent(j)-fTrue*nm)*(h[1]->GetBinContent(j)-fTrue*nm)/h[1]->GetBinContent(j); | |
585 | hchisq->Fill(chi2); | |
586 | chi2Sum+=chi2; | |
587 | ndof++; | |
588 | } | |
589 | ||
590 | *hchi2=*hchisq; | |
591 | delete hchisq; | |
592 | ||
593 | return chi2Sum; | |
594 | ||
595 | ||
596 | } | |
597 | ||
598 | void TestChi2(){ | |
599 | ||
600 | ||
601 | TH1F *h1=new TH1F("h1","h1",50,-10.,10.); | |
602 | TH1F *h2=new TH1F("h2","h2",50,-10.,10.); | |
603 | ||
604 | h1->FillRandom("gaus",50000); | |
605 | h2->FillRandom("gaus",80000); | |
606 | TH1F *hchi2=new TH1F(); | |
607 | Int_t ndof; | |
608 | Double_t chi2=GetChi2(h1,h2,hchi2,ndof); | |
609 | hchi2->Draw(); | |
610 | cout<<"The chi2 is "<<chi2; | |
611 | } |