]>
Commit | Line | Data |
---|---|---|
3ed21f0a | 1 | #include <fstream> |
2 | #include <Riostream.h> | |
82487ae7 | 3 | #include <TSystem.h> |
5b8639e7 | 4 | #include <TMath.h> |
3ed21f0a | 5 | #include <TH1F.h> |
6 | #include <TF1.h> | |
7 | #include <TFile.h> | |
8 | #include <TCanvas.h> | |
9 | #include <TClonesArray.h> | |
10 | #include <TStyle.h> | |
11 | #include <TLegend.h> | |
12 | #include <TGraphErrors.h> | |
13 | #include <TGraph.h> | |
14 | #include <TMultiGraph.h> | |
82487ae7 | 15 | #include <TKey.h> |
3ed21f0a | 16 | #include <TObjectTable.h> |
17 | #include <TDatabasePDG.h> | |
0adeb69c | 18 | #include <TMath.h> |
19 | #include <TPaveText.h> | |
20 | #include <TText.h> | |
3ed21f0a | 21 | |
22 | #include <AliMultiDimVector.h> | |
c24d2d9f | 23 | #include "AliHFMassFitter.h" |
3ed21f0a | 24 | #include <AliSignificanceCalculator.h> |
25 | ||
26 | #include <fstream> | |
27 | ||
82487ae7 | 28 | //global variables |
29 | Double_t nsigma=3; | |
5b8639e7 | 30 | Int_t decCh=2; |
31 | Int_t fitbtype=5; | |
82487ae7 | 32 | Int_t rebin=2; |
5b8639e7 | 33 | Double_t sigma=0.0005; |
82487ae7 | 34 | Int_t pdg; |
35 | Double_t mass; | |
36 | ||
37 | Double_t sigmaCut=0.035;//0.03; | |
38 | Double_t errSgnCut=0.4;//0.35; | |
39 | Double_t nSigmaMeanCut=4.;//3.; | |
40 | ||
0adeb69c | 41 | |
42 | ofstream outcheck; | |
43 | ofstream outdetail; | |
82487ae7 | 44 | |
45 | 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); | |
0adeb69c | 46 | Bool_t BinCounting(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, Int_t& status); |
82487ae7 | 47 | Bool_t MC(TH1F* hs,TH1F* hb, Double_t& sgn, Double_t& errsgn, Double_t& bkg, Double_t& errbkg, Double_t& sgnf, Double_t& errsgnf, Double_t& sigmaused, Int_t& status); |
48 | ||
3ed21f0a | 49 | //decCh: |
50 | //- 0 = kDplustoKpipi | |
51 | //- 1 = kD0toKpi | |
5b8639e7 | 52 | //- 2 = kDStartoKpipi |
3ed21f0a | 53 | //- 3 = kDstoKKpi |
54 | //- 4 = kD0toKpipipi | |
55 | //- 5 = kLambdactopKpi | |
56 | ||
4bcb0ebb | 57 | //Note: writefit=kTRUE writes the root files with the fit performed but it also draw all the canvas, so if your computer is not powerfull enough I suggest to run it in batch mode (root -b) |
58 | ||
5b8639e7 | 59 | // Method using this code: |
60 | // root -b | |
61 | // .L macros/LoadLibraries | |
62 | // Loadlibraries() | |
63 | // .L macros/CharmCutsOptimization | |
64 | // charmCutsOptimization(/*options*/) | |
65 | ||
66 | Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-particle, "P" particle*/,TString centr="no",Bool_t writefit=kTRUE,Int_t minentries=50,Double_t *rangefit=0x0, Bool_t useBinCounting=kTRUE){ | |
67 | ||
0adeb69c | 68 | outcheck.open("output.dat",ios::out); |
69 | outdetail.open("discarddetails.dat",ios::out); | |
cb3b714e | 70 | |
71 | gStyle->SetFrameBorderMode(0); | |
72 | gStyle->SetCanvasColor(0); | |
73 | gStyle->SetFrameFillColor(0); | |
74 | gStyle->SetTitleFillColor(0); | |
75 | gStyle->SetStatColor(0); | |
76 | ||
0adeb69c | 77 | //~/Lavoro/PbPb/tagli/SIGAOD33/mar02/cent3070/ |
3ed21f0a | 78 | TString filename="AnalysisResults.root",dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv"; |
79 | ||
82487ae7 | 80 | TString hnamemass="hMass_",hnamesgn="hSig_",hnamebkg="hBkg_"; |
3ed21f0a | 81 | |
3ed21f0a | 82 | |
83 | switch (decCh) { | |
84 | case 0: | |
85 | listname+="Dplus"; | |
86 | mdvlistname+="Dplus"; | |
87 | pdg=411; | |
88 | break; | |
89 | case 1: | |
90 | listname+="D0"; | |
91 | mdvlistname+="D0"; | |
92 | pdg=421; | |
93 | break; | |
94 | case 2: | |
5b8639e7 | 95 | listname+="Dstar0100"; |
96 | mdvlistname+="Dstar0100"; | |
3ed21f0a | 97 | pdg=413; |
98 | break; | |
99 | case 3: | |
100 | listname+="Ds"; | |
101 | mdvlistname+="Ds"; | |
102 | pdg=431; | |
103 | break; | |
104 | case 4: | |
105 | listname+="D04"; | |
106 | mdvlistname+="D04"; | |
107 | pdg=421; | |
108 | break; | |
109 | case 5: | |
110 | listname+="Lc"; | |
111 | mdvlistname+="Lc"; | |
112 | pdg=4122; | |
113 | break; | |
114 | default: | |
115 | cout<<decCh<<" is not allowed as decay channel "<<endl; | |
116 | return kFALSE; | |
117 | } | |
118 | mass=TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
5b8639e7 | 119 | if (decCh==2) mass=(TDatabasePDG::Instance()->GetParticle(pdg)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); |
3ed21f0a | 120 | |
cb3b714e | 121 | if(part!="both"){ |
122 | listname.Append(part); | |
123 | mdvlistname.Append(part); | |
124 | } | |
82487ae7 | 125 | if(centr!="no"){ |
126 | listname.Append(centr); | |
127 | mdvlistname.Append(centr); | |
128 | } | |
3ed21f0a | 129 | cout<<"Mass = "<<mass<<endl; |
130 | ||
3ed21f0a | 131 | Int_t countFitFail=0,countSgnfFail=0,countNoHist=0,countBkgOnly=0; |
3ed21f0a | 132 | |
133 | outcheck<<"ptbin\tmdvGlobAddr\thistIndex\tSignif\tS\tB"<<endl; | |
82487ae7 | 134 | outdetail<<"ptbin\tmdvGlobAddr\thistIndex\trelErrS\t\tmean_F-mass (mass = "<<mass<<")"<<endl; |
3ed21f0a | 135 | TFile *fin=new TFile(filename.Data()); |
136 | if(!fin->IsOpen()){ | |
137 | cout<<"File "<<filename.Data()<<" not found"<<endl; | |
138 | return kFALSE; | |
139 | } | |
140 | ||
141 | TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname); | |
142 | if(!dir){ | |
143 | cout<<"Directory "<<dirname<<" not found"<<endl; | |
144 | return kFALSE; | |
145 | } | |
146 | ||
147 | TList* histlist= (TList*)dir->Get(listname); | |
148 | if(!histlist) { | |
149 | cout<<listname<<" doesn't exist"<<endl; | |
150 | return kFALSE; | |
151 | } | |
152 | ||
153 | TList* listamdv= (TList*)dir->Get(mdvlistname); | |
154 | if(!listamdv) { | |
155 | cout<<mdvlistname<<" doesn't exist"<<endl; | |
156 | return kFALSE; | |
157 | } | |
158 | ||
6b64765b | 159 | TH1F* hstat=(TH1F*)histlist->FindObject("fHistNEvents"); |
160 | TCanvas *cst=new TCanvas("hstat","Summary of statistics"); | |
161 | if(hstat) { | |
162 | cst->cd(); | |
163 | cst->SetGrid(); | |
164 | hstat->Draw("htext0"); | |
82487ae7 | 165 | cst->SaveAs("hstat.png"); |
6b64765b | 166 | }else{ |
167 | cout<<"Warning! fHistNEvents not found in "<<listname.Data()<<endl; | |
168 | } | |
169 | ||
170 | Bool_t isMC=kFALSE; | |
82487ae7 | 171 | TH1F* htestIsMC=(TH1F*)histlist->FindObject("hSig_0"); |
6b64765b | 172 | if(htestIsMC) isMC=kTRUE; |
173 | ||
3ed21f0a | 174 | Int_t nptbins=listamdv->GetEntries(); |
6b64765b | 175 | Int_t nhist=(histlist->GetEntries()-1);//-1 because of fHistNevents |
176 | if(isMC) nhist/=4; ///4 because hMass_, hSgn_,hBkg_,hRfl_ | |
3ed21f0a | 177 | Int_t count=0; |
178 | Int_t *indexes= new Int_t[nhist]; | |
179 | //initialize indexes[i] to -1 | |
180 | for(Int_t i=0;i<nhist;i++){ | |
181 | indexes[i]=-1; | |
182 | } | |
183 | ||
184 | TFile* fout=new TFile(Form("outputSignifMaxim.root"),"recreate"); | |
185 | ||
cb3b714e | 186 | TH1F** hSigma=new TH1F*[nptbins]; |
82487ae7 | 187 | TH1F* hstatus=new TH1F("hstatus","Flag status",6,-0.5,5.5); |
188 | hstatus->GetXaxis()->SetBinLabel(1,"fit fail"); | |
189 | hstatus->GetXaxis()->SetBinLabel(2,"fit ok and good results"); | |
190 | hstatus->GetXaxis()->SetBinLabel(3,"quality requirements not satisfied"); | |
191 | hstatus->GetXaxis()->SetBinLabel(4,"only bkg fit ok"); | |
192 | hstatus->GetXaxis()->SetBinLabel(5,"negative signif"); | |
193 | hstatus->GetXaxis()->SetBinLabel(6,Form("< %d entries",minentries)); | |
cb3b714e | 194 | |
3ed21f0a | 195 | //Check wheter histograms are filled |
82487ae7 | 196 | if(isData){ |
197 | for(Int_t i=0;i<nhist;i++){ | |
198 | TString name=Form("%s%d",hnamemass.Data(),i); | |
199 | TH1F* h=(TH1F*)histlist->FindObject(name.Data()); | |
3ed21f0a | 200 | |
82487ae7 | 201 | if(!h){ |
202 | cout<<name<<" not found"<<endl; | |
203 | continue; | |
204 | } | |
3ed21f0a | 205 | |
82487ae7 | 206 | if(h->GetEntries()>minentries){ |
207 | //cout<<"Entries = "<<h->GetEntries()<<endl; | |
208 | if (h->Integral() > minentries){ | |
209 | cout<<i<<") Integral = "<<h->Integral()<<endl; | |
210 | indexes[i]=i; | |
211 | count++; | |
212 | } | |
3ed21f0a | 213 | } |
214 | } | |
82487ae7 | 215 | |
3ed21f0a | 216 | |
82487ae7 | 217 | cout<<"There are "<<count<<" histogram with more than "<<minentries<<" entries"<<endl; |
218 | if(count==0) { | |
219 | cout<<"No histogram to draw..."<<endl; | |
220 | return kFALSE; | |
221 | } | |
3ed21f0a | 222 | } |
3ed21f0a | 223 | //create multidimvectors |
224 | ||
225 | //for(Int_t i=0;i<1;i++){ | |
226 | for(Int_t i=0;i<nptbins;i++){ | |
227 | ||
228 | //multidimvectors for signal | |
229 | AliMultiDimVector *mdvS=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",i)); | |
230 | TString name=mdvS->GetName(),nameErr="err",setname=""; | |
231 | ||
232 | setname=Form("S%s",name.Data()); | |
233 | mdvS->SetName(setname.Data()); | |
4bcb0ebb | 234 | outcheck<<"\n"<<mdvS->GetPtLimit(0)<<" < Pt <"<<mdvS->GetPtLimit(1)<<endl; |
cb3b714e | 235 | |
3ed21f0a | 236 | AliMultiDimVector *mdvSerr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); |
237 | setname=Form("%sS%s",nameErr.Data(),name.Data()); | |
238 | mdvSerr->SetName(setname.Data()); | |
239 | ||
240 | //multidimvectors for background | |
241 | setname=Form("B%s",name.Data()); | |
242 | AliMultiDimVector *mdvB=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
243 | ||
244 | AliMultiDimVector *mdvBerr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
245 | setname=Form("%sB%s",nameErr.Data(),name.Data()); | |
246 | mdvBerr->SetName(setname.Data()); | |
247 | ||
248 | //multidimvectors for significance | |
249 | setname=Form("Sgf%s",name.Data()); | |
250 | AliMultiDimVector *mdvSgnf=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
251 | ||
252 | AliMultiDimVector *mdvSgnferr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
253 | setname=Form("%sSgf%s",nameErr.Data(),name.Data()); | |
254 | mdvSgnferr->SetName(setname.Data()); | |
255 | ||
cb3b714e | 256 | hSigma[i]=new TH1F(Form("hSigmapt%d",i),Form("Sigma distribution pt bin %d (%.1f < pt < %.1f)",i,mdvSgnf->GetPtLimit(0),mdvSgnf->GetPtLimit(1)), 200,0.,0.05); |
257 | ||
3ed21f0a | 258 | Int_t nhistforptbin=mdvS->GetNTotCells(); |
259 | //Int_t nvarsopt=mdvS->GetNVariables(); | |
260 | ||
261 | cout<<"nhistforptbin = "<<nhistforptbin<<endl; | |
262 | ||
263 | //loop on all histograms and do AliHFMassFitter | |
264 | //for(Int_t ih=0;ih<1;ih++){ | |
265 | for(Int_t ih=0;ih<nhistforptbin;ih++){ | |
266 | printf("Analyzing indexes[%d] = %d \n",ih+i*nhistforptbin,indexes[ih+i*nhistforptbin]); | |
82487ae7 | 267 | Int_t status=-1; |
268 | if(isData && indexes[ih+i*nhistforptbin] == -1) { | |
269 | status=5; | |
270 | mdvSgnferr->SetElement(ih,0); | |
271 | mdvS->SetElement(ih,0); | |
272 | mdvSerr->SetElement(ih,0); | |
273 | mdvB->SetElement(ih,0); | |
274 | mdvBerr->SetElement(ih,0); | |
3ed21f0a | 275 | |
82487ae7 | 276 | continue; |
277 | } | |
278 | outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]; | |
279 | TString name; | |
280 | TH1F* h=0x0; | |
281 | TH1F* g=0x0; | |
282 | Double_t signif=0, signal=0, background=0, errSignif=0, errSignal=0, errBackground=0,sigmafit=0; | |
283 | ||
284 | if(isData){ | |
285 | name=Form("%s%d",hnamemass.Data(),indexes[ih+i*nhistforptbin]); | |
286 | h=(TH1F*)histlist->FindObject(name.Data()); | |
287 | if(!h)continue; | |
0adeb69c | 288 | if(useBinCounting) { |
289 | if (h->GetEntries() >= minentries) | |
290 | BinCounting(h,rangefit,writefit,signal,errSignal,background,errBackground,signif,errSignif,status); | |
291 | } else | |
292 | Data(h,rangefit,writefit,signal,errSignal,background,errBackground,signif,errSignif,sigmafit,status); | |
82487ae7 | 293 | }else{ |
294 | name=Form("%s%d",hnamesgn.Data(),ih+i*nhistforptbin); | |
295 | h=(TH1F*)histlist->FindObject(name.Data()); | |
296 | if(!h){ | |
297 | cout<<name.Data()<<" not found"<<endl; | |
298 | continue; | |
3ed21f0a | 299 | } |
82487ae7 | 300 | name=Form("%s%d",hnamebkg.Data(),ih+i*nhistforptbin); |
301 | g=(TH1F*)histlist->FindObject(name.Data()); | |
302 | if(!g){ | |
303 | cout<<name.Data()<<" not found"<<endl; | |
304 | continue; | |
cb3b714e | 305 | } |
5b8639e7 | 306 | |
82487ae7 | 307 | MC(h,g,signal,errSignal,background,errBackground,signif,errSignif,sigmafit,status); |
5b8639e7 | 308 | |
82487ae7 | 309 | } |
5b8639e7 | 310 | |
82487ae7 | 311 | hstatus->Fill(status); |
312 | ||
313 | if(status==1){ | |
314 | mdvSgnf->SetElement(ih,signif); | |
315 | mdvSgnferr->SetElement(ih,errSignif); | |
316 | mdvS->SetElement(ih,signal); | |
317 | mdvSerr->SetElement(ih,errSignal); | |
318 | mdvB->SetElement(ih,background); | |
319 | mdvBerr->SetElement(ih,errBackground); | |
320 | hSigma[i]->Fill(sigmafit); | |
82487ae7 | 321 | }else{ |
cb3b714e | 322 | mdvSgnf->SetElement(ih,0); |
323 | mdvSgnferr->SetElement(ih,0); | |
324 | mdvS->SetElement(ih,0); | |
325 | mdvSerr->SetElement(ih,0); | |
326 | mdvB->SetElement(ih,0); | |
327 | mdvBerr->SetElement(ih,0); | |
82487ae7 | 328 | if(status==3){ |
329 | mdvB->SetElement(ih,background); | |
330 | mdvBerr->SetElement(ih,errBackground); | |
331 | } | |
3ed21f0a | 332 | |
82487ae7 | 333 | } |
3ed21f0a | 334 | |
cb3b714e | 335 | } |
82487ae7 | 336 | |
cb3b714e | 337 | fout->cd(); |
338 | mdvS->Write(); | |
339 | mdvB->Write(); | |
340 | mdvSgnf->Write(); | |
cb3b714e | 341 | mdvSerr->Write(); |
342 | mdvBerr->Write(); | |
343 | mdvSgnferr->Write(); | |
344 | hSigma[i]->Write(); | |
5b8639e7 | 345 | |
cb3b714e | 346 | } |
82487ae7 | 347 | |
348 | ||
349 | TCanvas *cinfo=new TCanvas("cinfo","Status"); | |
350 | cinfo->cd(); | |
351 | cinfo->SetGrid(); | |
352 | hstatus->Draw("htext0"); | |
82487ae7 | 353 | fout->cd(); |
354 | hstatus->Write(); | |
cb3b714e | 355 | fout->Close(); |
cb3b714e | 356 | outcheck<<"\nSummary:\n - Total number of histograms: "<<nhist<<"\n - "<<count<<" histograms with more than "<<minentries<<" entries; \n - Too few entries in histo "<<countNoHist<<" times;\n - Fit failed "<<countFitFail<<" times \n - no sense Signal/Background/Significance "<<countSgnfFail<<" times\n - only background "<<countBkgOnly<<" times"<<endl; |
357 | outcheck.close(); | |
5b8639e7 | 358 | |
359 | cout << "test1\n"; | |
cb3b714e | 360 | return kTRUE; |
3ed21f0a | 361 | } |
362 | ||
5b8639e7 | 363 | |
82487ae7 | 364 | //this function fit the hMass histograms |
365 | //status = 0 -> fit fail | |
366 | // 1 -> fit ok and good results | |
367 | // 2 -> quality requirements not satisfied, try to fit with bkg only | |
368 | // 3 -> only bkg fit ok | |
369 | // 4 -> negative signif | |
370 | // 5 -> not enough entries in the hisotgram | |
371 | 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){ | |
372 | Int_t nbin=h->GetNbinsX(); | |
373 | Double_t min=h->GetBinLowEdge(7); | |
374 | Double_t max=h->GetBinLowEdge(nbin-5)+h->GetBinWidth(nbin-5); | |
375 | ||
5b8639e7 | 376 | if(decCh != 2) min = h->GetBinLowEdge(1); |
377 | else min = TDatabasePDG::Instance()->GetParticle(211)->Mass(); | |
82487ae7 | 378 | max=h->GetBinLowEdge(nbin+1); |
379 | ||
380 | if(rangefit) { | |
381 | min=rangefit[0]; | |
382 | max=rangefit[1]; | |
383 | } | |
384 | ||
385 | AliHFMassFitter fitter(h,min, max,rebin,fitbtype); | |
386 | fitter.SetInitialGaussianMean(mass); | |
387 | fitter.SetInitialGaussianSigma(sigma); | |
388 | ||
389 | //if(ih==0) fitter.InitNtuParam(Form("ntuPtbin%d",i)); | |
390 | // fitter.SetHisto(h); | |
391 | // fitter.SetRangeFit(min,max); | |
392 | //fitter.SetRangeFit(1.68,2.05); | |
393 | ||
394 | //fitter.SetType(fitbtype,0); | |
395 | ||
396 | Bool_t ok=fitter.MassFitter(kFALSE); | |
397 | if(!ok){ | |
398 | cout<<"FIT NOT OK!"<<endl; | |
399 | //countBkgOnly++; | |
400 | //outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t 0\t xxx"<<"\t bkgonly"<<endl; | |
401 | outcheck<<"\t 0\t xxx"<<"\t failed"<<endl; | |
402 | status=0; | |
403 | return kFALSE; | |
404 | }else{ //fit ok! | |
405 | ||
406 | if(writefit) fitter.WriteCanvas(h->GetName(),"./",nsigma); | |
407 | fitter.Signal(nsigma,sgn,errsgn); | |
408 | fitter.Background(nsigma,bkg,errbkg); | |
409 | Double_t meanfit=fitter.GetMean(); | |
410 | sigmafit=fitter.GetSigma(); | |
411 | ||
412 | ||
413 | //if(ok==kTRUE && ( (sigmafit < 0.03) || (sigmafit < 0.04 && mdvS->GetPtLimit(0)>8.)) && sgn > 0 && bkg > 0){ | |
414 | if(ok==kTRUE && ( (sigmafit < sigmaCut) ) && sgn > 0 && bkg > 0){ | |
415 | Double_t errmeanfit=fitter.GetMassFunc()->GetParError(fitter.GetNFinalPars()-2); | |
416 | fitter.Significance(nsigma,sgnf,errsgnf); | |
417 | if(sgnf >0){ | |
418 | ||
419 | if(errsgn/sgn < errSgnCut && /*TMath::Abs(meanfit-mass)<0.015*/TMath::Abs(meanfit-mass)/errmeanfit < nSigmaMeanCut){ | |
420 | //outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t"<<signif<<" +- "<<errSignif<<"\t"<<sgn<<" +- "<<errsgn<<"\t"<<bkg<<" +- "<<errbkg<<endl; | |
421 | outcheck<<"\t\t "<<sgnf<<" +- "<<errsgnf<<"\t"<<sgn<<" +- "<<errsgn<<"\t"<<bkg<<" +- "<<errbkg<<endl; | |
422 | status=1; | |
423 | ||
424 | }else{ | |
425 | status=2; | |
426 | //outdetail<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t"<<errsgn/sgn<<"\t\t "<<(meanfit-mass)/errmeanfit<<endl; | |
427 | outdetail<<"\t\t "<<errsgn/sgn<<"\t\t "<<(meanfit-mass)/errmeanfit<<endl; | |
428 | ok=fitter.RefitWithBkgOnly(kFALSE); | |
429 | if (ok){ | |
430 | status=3; | |
431 | //countBkgOnly++; | |
432 | Double_t bkg=0,errbkg=0.; | |
433 | Double_t nsigmarange[2]={mass-nsigma*sigma,mass+nsigma*sigma}; | |
434 | fitter.Background(nsigmarange[0],nsigmarange[1],bkg,errbkg); | |
435 | //outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t 0\t "<<bkg <<"\t bkgonly"<<endl; | |
436 | outcheck<<"\t\t 0\t "<<bkg <<"\t bkgonly"<<endl; | |
437 | }else{ | |
438 | //outdetail<<i<<"\t\t "<<ih<<"\t\tnot able to refit with bkg obly"<<endl; | |
439 | outdetail<<"\t\t \t\tnot able to refit with bkg obly"<<endl; | |
440 | status=0; | |
441 | return kFALSE; | |
442 | } | |
443 | }//only bkg | |
444 | }//check signif>0 | |
445 | else{ | |
446 | status=4; | |
447 | //countSgnfFail++; | |
448 | cout<<"Setting to 0 (fitter results meaningless)"<<endl; | |
449 | outcheck<<"\t S || B || sgnf negative"; | |
450 | ||
451 | return kFALSE; | |
452 | } | |
453 | } //end fit ok! | |
454 | } | |
455 | outcheck<<endl; | |
456 | return kTRUE; | |
457 | } | |
458 | ||
459 | ||
460 | ||
461 | //this function counts the entries in hSgn and hBgk | |
462 | Bool_t MC(TH1F* hs,TH1F* hb, Double_t& sgn, Double_t& errsgn, Double_t& bkg, Double_t& errbkg, Double_t& sgnf, Double_t& errsgnf, Double_t& sigmaused, Int_t& status){ | |
463 | ||
464 | //do we want to use a fixed sigma or take the standard deviation of the signal histogram? | |
465 | sigmaused=hs->GetRMS(); | |
466 | //sigmaused=sigma; | |
467 | ||
468 | Double_t nsigmarange[2]={mass-nsigma*sigmaused,mass+nsigma*sigmaused}; | |
469 | cout<<"from "<<nsigmarange[0]<<" to "<<nsigmarange[1]<<endl; | |
470 | ||
471 | Int_t binnsigmarange[2]={hs->FindBin(nsigmarange[0]),hs->FindBin(nsigmarange[1])};//for bkg histo it's the same | |
472 | cout<<"bins "<<binnsigmarange[0]<<" e "<<binnsigmarange[1]<<endl; | |
473 | ||
474 | sgn=hs->Integral(binnsigmarange[0],binnsigmarange[1]); | |
475 | errsgn=TMath::Sqrt(sgn); | |
476 | bkg=hb->Integral(binnsigmarange[0],binnsigmarange[1]); | |
477 | errbkg=TMath::Sqrt(bkg); | |
478 | if(sgn+bkg>0.) sgnf=sgn/TMath::Sqrt(sgn+bkg); | |
479 | else { | |
480 | status=2; | |
481 | return kFALSE; | |
482 | } | |
483 | errsgnf=TMath::Sqrt(sgnf*sgnf/(sgn+bkg)/(sgn+bkg)*(1/4.*errsgn*errsgn+errbkg*errbkg)+sgnf*sgnf/sgn/sgn*errsgn*errsgn); | |
484 | status=1; | |
485 | return kTRUE; | |
486 | ||
487 | } | |
3ed21f0a | 488 | |
0adeb69c | 489 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: |
490 | // par[0], par[1] expo params, par[2], par[3] exclusion range | |
491 | Bool_t reject = true; | |
492 | Double_t ExpoBkgWoPeak(Double_t *x, Double_t *par){ | |
493 | ||
494 | if( reject && x[0]>par[2] && x[0]<par[3] ){ | |
495 | TF1::RejectPoint(); | |
496 | return 0; | |
497 | } | |
498 | return par[0] + TMath::Exp(par[1]*x[0]) ; | |
499 | ||
500 | } | |
501 | ||
5b8639e7 | 502 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: |
503 | // par[0], par[1] expo params, par[2], par[3] exclusion range | |
504 | ||
505 | Double_t PowerBkgWoPeak(Double_t *x, Double_t *par){ | |
506 | ||
507 | if( reject && x[0]>par[2] && x[0]<par[3] ){ | |
508 | TF1::RejectPoint(); | |
509 | return 0; | |
510 | } | |
511 | ||
512 | Double_t xminusmpi = x[0]-TDatabasePDG::Instance()->GetParticle(211)->Mass(); | |
513 | return par[0]*TMath::Power(TMath::Abs(xminusmpi),par[1]) ; | |
514 | } | |
515 | ||
516 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: | |
517 | // par[0], par[1] expo params, par[2], par[3] exclusion range | |
518 | ||
519 | Double_t PowerExpoBkgWoPeak(Double_t *x, Double_t *par){ | |
520 | ||
521 | if( reject && x[0]>par[2] && x[0]<par[3] ){ | |
522 | TF1::RejectPoint(); | |
523 | return 0; | |
524 | } | |
525 | Double_t xminusmpi = x[0]-TDatabasePDG::Instance()->GetParticle(211)->Mass(); | |
526 | return par[0]*TMath::Sqrt(xminusmpi)*TMath::Exp(-1.*par[1]*xminusmpi); | |
527 | } | |
528 | ||
529 | ||
0adeb69c | 530 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: |
531 | //this function fit the hMass histograms | |
532 | //status = 0 -> fit fail | |
533 | // 1 -> fit ok and good results | |
534 | // 2 -> negative signif | |
535 | Bool_t BinCounting(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, Int_t& status){ | |
536 | ||
5b8639e7 | 537 | |
a7665a76 | 538 | Double_t min, max; |
0adeb69c | 539 | Int_t nbin=h->GetNbinsX(); |
a7665a76 | 540 | if(decCh != 2) min = h->GetBinLowEdge(1); |
541 | else min = TDatabasePDG::Instance()->GetParticle(211)->Mass(); | |
542 | max=h->GetBinLowEdge(nbin+1); | |
0adeb69c | 543 | |
544 | if(rangefit) { | |
545 | min=rangefit[0]; | |
546 | max=rangefit[1]; | |
547 | } | |
548 | ||
0adeb69c | 549 | reject = true; |
5b8639e7 | 550 | //Bkg fit : exponential = A*exp(B*x) |
a7665a76 | 551 | TF1 *fBkgFit; |
552 | if(decCh != 2) fBkgFit = new TF1("fBkgFit",ExpoBkgWoPeak,min,max,4); | |
553 | else { | |
554 | if (fitbtype == 5) fBkgFit = new TF1("fBkgFit",PowerExpoBkgWoPeak,min,max,4); //Bkg fit : PowerExpo = A*(x-mpi)^(1/2)*exp(-B*(x-mpi)) | |
555 | else fBkgFit = new TF1("fBkgFit",PowerBkgWoPeak,min,max,4); //Bkg fit : Power = A*(x-mpi)^B | |
5b8639e7 | 556 | } |
557 | ||
0adeb69c | 558 | fBkgFit->FixParameter(2,mass-nsigma*sigma); |
559 | fBkgFit->FixParameter(3,mass+nsigma*sigma); | |
560 | TFitResultPtr r = h->Fit(fBkgFit,"RS+"); | |
561 | Int_t ok = r; | |
562 | ||
0adeb69c | 563 | if(ok==-1){ |
564 | cout<<"FIT NOT OK!"<<endl; | |
565 | cout<<"\t 0\t xxx"<<"\t failed"<<endl; | |
566 | status=0; | |
567 | return kFALSE; | |
568 | } | |
0adeb69c | 569 | |
f79ea837 | 570 | reject = false; |
a7665a76 | 571 | TF1 *fBkgFct; |
572 | if(decCh !=2) fBkgFct = new TF1("fBkgFct",ExpoBkgWoPeak,min,max,4); | |
5b8639e7 | 573 | else { |
a7665a76 | 574 | if (fitbtype == 5) fBkgFct = new TF1("fBkgFct",PowerExpoBkgWoPeak,min,max,4); |
575 | else fBkgFct = new TF1("fBkgFct",PowerBkgWoPeak,min,max,4); | |
5b8639e7 | 576 | } |
f79ea837 | 577 | fBkgFct->SetLineStyle(2); |
578 | for(Int_t i=0; i<4; i++) fBkgFct->SetParameter(i,fBkgFit->GetParameter(i)); | |
579 | h->GetListOfFunctions()->Add(fBkgFct); | |
580 | TH1F * hBkgFct = (TH1F*)fBkgFct->GetHistogram(); | |
581 | ||
5b8639e7 | 582 | |
f79ea837 | 583 | status = 1; |
584 | Double_t binStartCount = h->FindBin(mass-nsigma*sigma); | |
585 | Double_t binEndCount = h->FindBin(mass+nsigma*sigma); | |
586 | Double_t counts=0., bkgcounts=0., errcounts=0., errbkgcounts=0.; | |
587 | for (Int_t ibin = binStartCount; ibin<=binEndCount; ibin++) { | |
588 | counts += h->GetBinContent( ibin ); | |
589 | errcounts += counts ; | |
590 | Double_t center = h->GetBinCenter(ibin); | |
5b8639e7 | 591 | bkgcounts += hBkgFct->GetBinContent( hBkgFct->FindBin(center) ); |
592 | ||
f79ea837 | 593 | errbkgcounts += bkgcounts ; |
594 | } | |
5b8639e7 | 595 | |
f79ea837 | 596 | bkg = bkgcounts; |
597 | errbkg = TMath::Sqrt( errbkgcounts ); | |
598 | sgn = counts - bkg ; | |
599 | if(sgn<0) status = 2; // significance < 0 | |
600 | errsgn = TMath::Sqrt( counts + errbkg*errbkg ); | |
601 | sgnf = sgn / TMath::Sqrt( sgn + bkg ); | |
602 | errsgnf = TMath::Sqrt( sgnf*sgnf/(sgn+bkg)/(sgn+bkg)*(1/4.*errsgn*errsgn+errbkg*errbkg) + sgnf*sgnf/sgn/sgn*errsgn*errsgn ); | |
5b8639e7 | 603 | cout << " Signal "<<sgn<<" +- "<<errsgn<<", bkg "<<bkg<<" +- "<<errbkg<<", significance "<<sgnf<<" +- "<<errsgnf<<endl; |
f79ea837 | 604 | |
605 | if(writefit) { | |
5b8639e7 | 606 | |
f79ea837 | 607 | TString filename = Form("%sMassFit.root",h->GetName()); |
5b8639e7 | 608 | |
f79ea837 | 609 | TFile* outputcv = new TFile(filename.Data(),"recreate"); |
5b8639e7 | 610 | |
f79ea837 | 611 | TCanvas* c = new TCanvas(); |
5b8639e7 | 612 | |
f79ea837 | 613 | c->SetName(Form("%s",h->GetName())); |
5b8639e7 | 614 | |
f79ea837 | 615 | h->Draw(); |
5b8639e7 | 616 | |
f79ea837 | 617 | TPaveText *pavetext=new TPaveText(0.4,0.7,0.85,0.9,"NDC"); |
618 | pavetext->SetBorderSize(0); | |
619 | pavetext->SetFillStyle(0); | |
620 | pavetext->AddText(Form("Signal = %4.2e #pm %4.2e",sgn,errsgn)); | |
621 | pavetext->AddText(Form("Bkg = %4.2e #pm %4.2e",bkg,errbkg)); | |
622 | pavetext->AddText(Form("Signif = %3.2f #pm %3.2f",sgnf,errsgnf)); | |
623 | c->cd(); | |
5b8639e7 | 624 | |
f79ea837 | 625 | pavetext->DrawClone(); |
5b8639e7 | 626 | |
f79ea837 | 627 | outputcv->cd(); |
5b8639e7 | 628 | |
f79ea837 | 629 | c->Write(); |
5b8639e7 | 630 | |
631 | ||
f79ea837 | 632 | outputcv->Close(); |
5b8639e7 | 633 | |
f79ea837 | 634 | delete outputcv; |
5b8639e7 | 635 | |
f79ea837 | 636 | delete c; |
0adeb69c | 637 | } |
f79ea837 | 638 | |
5b8639e7 | 639 | |
0adeb69c | 640 | delete fBkgFit; |
5b8639e7 | 641 | |
0adeb69c | 642 | delete fBkgFct; |
643 | ||
5b8639e7 | 644 | |
0adeb69c | 645 | return kTRUE; |
646 | } | |
647 | ||
3ed21f0a | 648 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: |
649 | ||
650 | // which=0 plot significance | |
651 | // =1 plot signal | |
652 | // =2 plot background | |
c24d2d9f | 653 | // =3 plot S/B |
6b64765b | 654 | // maximize = kTRUE (default) if you want to fix the step of the variables not shown to the value that maximize the significance. Note that these values are saved in fixedvars.dat |
655 | // readfromfile = kTRUE (default is kFALSE) if you want to read the value fixed in a previous run of this function (e.g. significance or signal maximization) | |
3ed21f0a | 656 | |
6b64765b | 657 | |
0adeb69c | 658 | void showMultiDimVector(Int_t n=2,Int_t which=0, Bool_t plotErrors=kFALSE,Bool_t readfromfile=kFALSE, Bool_t fixedrange=kFALSE, Bool_t fixedplane=kFALSE){ |
3ed21f0a | 659 | |
660 | gStyle->SetCanvasColor(0); | |
661 | gStyle->SetFrameFillColor(0); | |
662 | gStyle->SetTitleFillColor(0); | |
663 | gStyle->SetOptStat(0); | |
664 | //gStyle->SetOptTitle(0); | |
cb3b714e | 665 | gStyle->SetFrameBorderMode(0); |
3ed21f0a | 666 | |
667 | TFile* fin=new TFile("outputSignifMaxim.root"); | |
668 | if(!fin->IsOpen()){ | |
669 | cout<<"outputSignifMaxim.root not found"<<endl; | |
670 | return; | |
671 | } | |
672 | ||
673 | if(n>2){ | |
674 | cout<<"Error! cannot show "<<n+1<<" dimentions"<<endl; | |
675 | return; | |
676 | } | |
677 | ||
c24d2d9f | 678 | TString name,title,namebis,shorttitle; |
3ed21f0a | 679 | switch (which){ |
680 | case 0: | |
681 | name="SgfmultiDimVectorPtBin"; | |
682 | title="Significance"; | |
c24d2d9f | 683 | shorttitle="Sgnf"; |
3ed21f0a | 684 | break; |
685 | case 1: | |
686 | name="SmultiDimVectorPtBin"; | |
687 | title="Signal"; | |
c24d2d9f | 688 | shorttitle="S"; |
3ed21f0a | 689 | break; |
690 | case 2: | |
691 | name="BmultiDimVectorPtBin"; | |
692 | title="Background"; | |
c24d2d9f | 693 | shorttitle="B"; |
3ed21f0a | 694 | break; |
695 | case 3: | |
c24d2d9f | 696 | name="SmultiDimVectorPtBin"; |
697 | namebis="BmultiDimVectorPtBin"; | |
698 | title="Signal over Background "; | |
699 | shorttitle="SoB"; | |
3ed21f0a | 700 | break; |
c24d2d9f | 701 | // case 4: |
702 | // name="errBmultiDimVectorPtBin"; | |
703 | // title="Background (error)"; | |
704 | // break; | |
3ed21f0a | 705 | } |
706 | ||
cb3b714e | 707 | if(plotErrors && which!=3 && n==2){ |
708 | name.Prepend("err"); | |
709 | title.Append(" Error") ; | |
710 | shorttitle.Append("Err"); | |
711 | } | |
712 | ||
713 | Int_t nptbins=50; | |
3ed21f0a | 714 | |
4bcb0ebb | 715 | for(Int_t ip=0;ip<=nptbins;ip++){ |
3ed21f0a | 716 | TString mdvname=Form("%s%d",name.Data(),ip); |
717 | AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname); | |
718 | if(!mdv){ | |
719 | nptbins=ip; | |
720 | cout<<"Number of pt bins "<<ip<<endl; | |
721 | break; | |
722 | } | |
723 | } | |
724 | ||
725 | cout<<"Projecting "<<title.Data()<<" with respect to the maximization variable(s) [chose]"<<endl; | |
726 | ||
727 | Int_t variable[2]; //no more than 2D | |
c24d2d9f | 728 | TString mdvname=Form("%s0",name.Data()), mdverrname="";//, mdvnamebis="", mdverrnamebis=""; |
3ed21f0a | 729 | AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname); |
730 | AliMultiDimVector* mdverr=0x0; | |
731 | if(!mdv){ | |
732 | cout<<mdvname.Data()<<" not found"<<endl; | |
733 | return; | |
734 | } | |
c24d2d9f | 735 | |
3ed21f0a | 736 | Int_t nvarsopt=mdv->GetNVariables(); |
6b64765b | 737 | //Int_t nfixed=nvarsopt-n; |
0adeb69c | 738 | |
6b64765b | 739 | Int_t fixedvars[nvarsopt]; |
740 | Int_t allfixedvars[nvarsopt*nptbins]; | |
0adeb69c | 741 | |
742 | ||
3ed21f0a | 743 | |
6b64765b | 744 | fstream writefixedvars; |
745 | if(readfromfile) { | |
746 | //open file in read mode | |
747 | writefixedvars.open("fixedvars.dat",ios::in); | |
748 | Int_t longi=0; | |
749 | while(writefixedvars){ | |
750 | writefixedvars>>allfixedvars[longi]; | |
751 | longi++; | |
752 | } | |
753 | } | |
754 | else { | |
755 | //open file in write mode | |
756 | writefixedvars.open("fixedvars.dat",ios::out); | |
3ed21f0a | 757 | } |
758 | ||
0adeb69c | 759 | Bool_t freevars[nvarsopt]; |
760 | ||
3ed21f0a | 761 | //ask variables for projection |
762 | for(Int_t k=0;k<nvarsopt;k++){ | |
3ed21f0a | 763 | cout<<k<<" "<<mdv->GetAxisTitle(k)<<endl; |
0adeb69c | 764 | freevars[k]=kTRUE; |
3ed21f0a | 765 | } |
766 | cout<<"Choose "<<n<<" variable(s)"<<endl; | |
767 | for(Int_t j=0;j<n;j++){ | |
768 | cout<<"var"<<j<<": "; | |
0adeb69c | 769 | cin>>variable[j]; |
770 | freevars[variable[j]]=kFALSE; | |
3ed21f0a | 771 | } |
772 | if(n==1) variable[1]=999; | |
773 | ||
3ed21f0a | 774 | TCanvas* cvpj=new TCanvas(Form("proj%d",variable[0]),Form("%s wrt %s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data())); |
6b64765b | 775 | |
3ed21f0a | 776 | TMultiGraph* mg=new TMultiGraph(Form("proj%d",variable[0]),Form("%s wrt %s;%s;%s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data(),(mdv->GetAxisTitle(variable[0])).Data(),title.Data())); |
777 | TLegend *leg=new TLegend(0.7,0.2,0.9,0.6,"Pt Bin"); | |
778 | leg->SetBorderSize(0); | |
779 | leg->SetFillStyle(0); | |
780 | ||
82487ae7 | 781 | //set scale |
782 | fstream writerange; | |
783 | Float_t axisrange[2*nptbins]; | |
784 | if(fixedrange) { | |
785 | //open file in read mode | |
786 | writerange.open("rangehistos.dat",ios::in); | |
787 | Int_t longi=0; | |
788 | while(writerange){ | |
789 | writerange>>axisrange[longi]; | |
790 | longi++; | |
791 | } | |
792 | } | |
793 | else { | |
794 | //open file in write mode | |
795 | writerange.open("rangehistos.dat",ios::out); | |
796 | } | |
797 | ||
3ed21f0a | 798 | for (Int_t i=0;i<nptbins;i++){ //loop on ptbins |
799 | cout<<"\nPtBin = "<<i<<endl; | |
800 | ||
801 | //using AliSignificanceCalculator | |
802 | ||
803 | TString nameS,nameB,nameerrS,nameerrB; | |
804 | nameS.Form("SmultiDimVectorPtBin%d",i); | |
805 | nameerrS.Form("errSmultiDimVectorPtBin%d",i); | |
806 | nameB.Form("BmultiDimVectorPtBin%d",i); | |
807 | nameerrB.Form("errBmultiDimVectorPtBin%d",i); | |
808 | ||
809 | AliMultiDimVector* mdvS=(AliMultiDimVector*)fin->Get(nameS.Data()); | |
810 | AliMultiDimVector* mdvB=(AliMultiDimVector*)fin->Get(nameB.Data()); | |
811 | AliMultiDimVector* mdvBerr=(AliMultiDimVector*)fin->Get(nameerrS.Data()); | |
812 | AliMultiDimVector* mdvSerr=(AliMultiDimVector*)fin->Get(nameerrB.Data()); | |
813 | if(!(mdvS && mdvB && mdvSerr && mdvBerr)){ | |
814 | cout<<"one of the multidimvector is not present"<<endl; | |
815 | return; | |
816 | } | |
817 | ||
818 | AliSignificanceCalculator *cal=new AliSignificanceCalculator(mdvS,mdvB,mdvSerr,mdvBerr,1.,1.); | |
819 | ||
3ed21f0a | 820 | AliMultiDimVector* mvess=cal->GetSignificanceError(); |
3ed21f0a | 821 | AliMultiDimVector* mvpur=cal->CalculatePurity(); |
822 | AliMultiDimVector* mvepur=cal->CalculatePurityError(); | |
6b64765b | 823 | |
3ed21f0a | 824 | Int_t ncuts=mdvS->GetNVariables(); |
825 | Int_t *maxInd=new Int_t[ncuts]; | |
826 | Float_t *cutvalues=new Float_t[ncuts]; | |
6b64765b | 827 | //init |
0adeb69c | 828 | // for(Int_t ind=0;ind<ncuts;ind++)maxInd[ind]=0; |
6b64765b | 829 | |
82487ae7 | 830 | Float_t sigMax0=cal->GetMaxSignificance(maxInd,0); //look better into this!! |
0adeb69c | 831 | |
3ed21f0a | 832 | for(Int_t ic=0;ic<ncuts;ic++){ |
833 | cutvalues[ic]=((AliMultiDimVector*)fin->Get(nameS.Data()))->GetCutValue(ic,maxInd[ic]); | |
6b64765b | 834 | |
835 | //setting step of fixed variables | |
836 | if(readfromfile){ //from file | |
837 | fixedvars[ic]=allfixedvars[i+ic]; | |
838 | } | |
839 | ||
0adeb69c | 840 | if(!readfromfile) { //using the values which maximize the significance |
841 | fixedvars[ic]=maxInd[ic]; | |
842 | //write to output fixedvars.dat | |
843 | writefixedvars<<fixedvars[ic]<<"\t"; | |
6b64765b | 844 | } |
3ed21f0a | 845 | } |
6b64765b | 846 | //output file: return after each pt bin |
0adeb69c | 847 | if(!readfromfile) writefixedvars<<endl; |
6b64765b | 848 | |
3ed21f0a | 849 | printf("Maximum of significance for Ptbin %d found in bin:\n",i); |
0adeb69c | 850 | for(Int_t ic=0;ic<ncuts;ic++)printf(" %d ",maxInd[ic]); |
851 | printf("\ncorresponding to cut:\n"); | |
852 | for(Int_t ic=0;ic<ncuts;ic++)printf(" %f ",cutvalues[ic]); | |
853 | ||
854 | printf("\nSignificance = %f +- %f\n",sigMax0,mvess->GetElement(maxInd,0)); | |
6b64765b | 855 | printf("Purity = %f +- %f\n",mvpur->GetElement(maxInd,0),mvepur->GetElement(maxInd,i)); |
0adeb69c | 856 | |
c24d2d9f | 857 | if(which==3){ |
858 | //mdv=0x0; | |
859 | mdv=cal->CalculateSOverB(); | |
860 | if(!mdv)cout<<mdv->GetName()<<" null"<<endl; | |
861 | //mdverr=0x0; | |
862 | mdverr=cal->CalculateSOverBError(); | |
863 | if(!mdverr)cout<<mdverr->GetName()<<" null"<<endl; | |
864 | }else{ | |
865 | ||
866 | //multidimvector | |
867 | mdvname=Form("%s%d",name.Data(),i); | |
868 | mdv=(AliMultiDimVector*)fin->Get(mdvname); | |
869 | if(!mdv)cout<<mdvname.Data()<<" not found"<<endl; | |
870 | ||
871 | //multidimvector of errors | |
872 | mdverrname=Form("err%s%d",name.Data(),i); | |
873 | mdverr=(AliMultiDimVector*)fin->Get(mdverrname); | |
874 | if(!mdverr)cout<<mdverrname.Data()<<" not found"<<endl; | |
875 | } | |
0adeb69c | 876 | printf("Global Address %d (%d)\n",(Int_t)mdv->GetGlobalAddressFromIndices(maxInd,0),(Int_t)mdv->GetNTotCells()*i+(Int_t)mdv->GetGlobalAddressFromIndices(maxInd,0)); |
cb3b714e | 877 | TString ptbinrange=Form("%.1f < p_{t} < %.1f GeV/c",mdv->GetPtLimit(0),mdv->GetPtLimit(1)); |
3ed21f0a | 878 | |
0adeb69c | 879 | Float_t maxval=0; |
880 | ||
3ed21f0a | 881 | if(n==2) { |
882 | gStyle->SetPalette(1); | |
0adeb69c | 883 | Int_t steps[2]; |
884 | Int_t nstep[2]={mdv->GetNCutSteps(variable[0]),mdv->GetNCutSteps(variable[1])}; | |
885 | ||
886 | TH2F* hproj=new TH2F(Form("hproj%d",i),Form("%s wrt %s vs %s (Ptbin%d %.1f<pt<%.1f);%s;%s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data(),i,mdv->GetPtLimit(0),mdv->GetPtLimit(1),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data()),nstep[0],mdv->GetMinLimit(variable[0]),mdv->GetMaxLimit(variable[0]),nstep[1],mdv->GetMinLimit(variable[1]),mdv->GetMaxLimit(variable[1])); | |
887 | if(fixedplane){ | |
888 | hproj=mdv->Project(variable[0],variable[1],fixedvars,0); | |
889 | hproj->SetTitle(Form("%s wrt %s vs %s (Ptbin%d %.1f<pt<%.1f);%s;%s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data(),i,mdv->GetPtLimit(0),mdv->GetPtLimit(1),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data())); | |
890 | }else{ | |
891 | for(Int_t ist1=0;ist1<nstep[0];ist1++){ | |
892 | steps[0]=ist1; | |
893 | Int_t fillbin1=ist1+1; | |
894 | if(mdv->GetCutValue(variable[0],0)>mdv->GetCutValue(variable[0],mdv->GetNCutSteps(variable[0])-1))fillbin1=nstep[0]-ist1; | |
895 | for(Int_t ist2=0;ist2<nstep[1];ist2++){ | |
896 | steps[1]=ist2; | |
897 | Int_t fillbin2=ist2+1; | |
898 | if(mdv->GetCutValue(variable[1],0)>mdv->GetCutValue(variable[1],mdv->GetNCutSteps(variable[1])-1))fillbin2=nstep[1]-ist2; | |
899 | Int_t* varmaxim=mdv->FindLocalMaximum(maxval,variable,steps,n,0); | |
900 | hproj->SetBinContent(fillbin1,fillbin2,maxval); | |
901 | delete varmaxim; | |
902 | } | |
903 | } | |
904 | } | |
82487ae7 | 905 | if(fixedrange) { |
906 | hproj->SetMinimum(axisrange[2*i]); | |
907 | hproj->SetMaximum(axisrange[2*i+1]); | |
908 | } else{ | |
909 | writerange<<hproj->GetMinimum()<<"\t"<<hproj->GetMinimum()<<endl; | |
910 | } | |
3ed21f0a | 911 | TCanvas* cvpj=new TCanvas(Form("proj%d%dpt%d",variable[0],variable[1],i),Form("%s wrt %s vs %s (Ptbin%d)",title.Data(),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data(),i)); |
912 | cvpj->cd(); | |
6b64765b | 913 | hproj->DrawClone("COLZtext"); |
c24d2d9f | 914 | cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName())); |
3ed21f0a | 915 | delete hproj; |
916 | } | |
917 | ||
918 | if(n==1){ | |
6b64765b | 919 | |
3ed21f0a | 920 | Int_t nbins=mdv->GetNCutSteps(variable[0]); |
921 | ||
922 | Double_t *x=new Double_t[nbins]; | |
923 | Double_t *y=new Double_t[nbins]; | |
924 | Double_t *errx=new Double_t[nbins]; | |
925 | Double_t *erry=new Double_t[nbins]; | |
926 | ||
927 | for(Int_t k=0;k<nbins;k++){ //loop on the steps (that is the bins of the graph) | |
6b64765b | 928 | //init |
3ed21f0a | 929 | x[k]=0;y[k]=0; |
930 | errx[k]=0;erry[k]=0; | |
0adeb69c | 931 | |
3ed21f0a | 932 | x[k]=mdv->GetCutValue(variable[0],k); |
933 | errx[k]=mdv->GetCutStep(variable[0])/2.; | |
0adeb69c | 934 | Int_t onevariable[1]={variable[0]}; |
935 | Int_t onestep[1]={k}; | |
936 | ULong64_t gladd; | |
937 | ||
938 | Float_t maxval; | |
939 | Int_t* varmaxim=mdv->FindLocalMaximum(maxval,onevariable,onestep,n,0); | |
940 | y[k]=maxval; | |
941 | gladd=mdv->GetGlobalAddressFromIndices(varmaxim,0); | |
942 | cout<<gladd<<endl; | |
943 | delete varmaxim; | |
944 | ||
945 | erry[k]=mdverr->GetElement(gladd); | |
3ed21f0a | 946 | |
6b64765b | 947 | cout<<mdv->GetAxisTitle(variable[0])<<" step "<<k<<" = "<<x[k]<<":"<<" y = "<<y[k]<<endl; |
3ed21f0a | 948 | } |
0adeb69c | 949 | |
6b64765b | 950 | cout<<"----------------------------------------------------------"<<endl; |
3ed21f0a | 951 | TGraphErrors* gr=new TGraphErrors(nbins,x,y,errx,erry); |
952 | gr->SetMarkerStyle(20+i); | |
953 | gr->SetMarkerColor(i+1); | |
cb3b714e | 954 | gr->SetLineColor(i+1); |
955 | if(i>=9){ | |
956 | gr->SetMarkerColor(i+2); | |
957 | gr->SetLineColor(i+2); | |
958 | } | |
3ed21f0a | 959 | gr->SetMinimum(0); |
960 | ||
961 | gr->SetName(Form("g1%d",i)); | |
962 | mg->Add(gr,"P"); | |
963 | leg->AddEntry(gr,ptbinrange.Data(),"p"); | |
964 | } | |
965 | } | |
966 | ||
967 | if(n==1){ | |
968 | cvpj->cd(); | |
969 | mg->Draw("A"); | |
970 | leg->Draw(); | |
c24d2d9f | 971 | cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName())); |
6b64765b | 972 | } else delete cvpj; |
cb3b714e | 973 | } |
974 | ||
82487ae7 | 975 | //draw sigma as a function of cuts |
976 | ||
977 | void DrawSigmas(TH2F* h2cuts){ | |
978 | TFile *fin=0x0; | |
979 | TString fittype="ExpFit"; | |
980 | Int_t ntot=5; | |
981 | if(fittype=="Pol2Fit") ntot=6; | |
982 | Int_t ihfirst=0,ihlast=1; //change this (must think on it and remember what I wanted to do!) | |
983 | for(Int_t ih=ihfirst;ih<ihlast;ih++){ | |
984 | fin=new TFile(Form("h%d%s.root",ih,fittype.Data())); | |
985 | if(!fin) continue; | |
986 | TCanvas *cv=(TCanvas*)fin->Get(Form("cv1%s%d",fittype.Data(),ih)); | |
987 | TF1* func=(TF1*)cv->FindObject("funcmass"); | |
988 | Int_t sigma=func->GetParameter(ntot-1); | |
989 | //h2cuts->SetBinContent(); | |
990 | //to be finished | |
991 | } | |
992 | } | |
993 | ||
cb3b714e | 994 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: |
995 | // Some methods to get the index of the histogram corresponding to a set of cuts | |
996 | ||
997 | // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root | |
998 | // ptbin = pt bin | |
999 | // indices = array of the index of the cut variables (the dimension of the array must be equal to the number of variables maximized) | |
4bcb0ebb | 1000 | |
cb3b714e | 1001 | Int_t GetNHistFromIndices(AliMultiDimVector* vct,Int_t ptbin,Int_t* indices){ |
1002 | cout<<"Calculating the index of the histogram corresponding to the following cut steps:"<<endl; | |
1003 | for(Int_t i=0;i<vct->GetNVariables();i++){ | |
1004 | cout<<vct->GetAxisTitle(i)<<" --> "<<indices[i]<<endl; | |
1005 | } | |
1006 | cout<<"Pt bin "<<ptbin<<" from "<<vct->GetPtLimit(0)<<" to "<<vct->GetPtLimit(1)<<endl; | |
1007 | cout<<"Info: total number of cells per multidim: "<<vct->GetNTotCells()<<endl; | |
1008 | ULong64_t glindex=vct->GetGlobalAddressFromIndices(indices,0); | |
1009 | cout<<"The histogram you want is:\t"; | |
1010 | return glindex+vct->GetNTotCells()*ptbin; | |
3ed21f0a | 1011 | } |
1012 | ||
cb3b714e | 1013 | // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root |
1014 | // ptbin = pt bin | |
1015 | // values = array of the cut values (the dimension of the array must be equal to the number of variables maximized) | |
1016 | ||
1017 | Int_t GetNHistFromValues(AliMultiDimVector* vct,Int_t ptbin,Float_t* values){ | |
1018 | cout<<"Calculating the index of the histogram corresponding to the following cut values:"<<endl; | |
1019 | for(Int_t i=0;i<vct->GetNVariables();i++){ | |
1020 | cout<<vct->GetAxisTitle(i)<<" --> "<<values[i]<<endl; | |
1021 | } | |
1022 | cout<<"Pt bin "<<ptbin<<" from "<<vct->GetPtLimit(0)<<" to "<<vct->GetPtLimit(1)<<endl; | |
1023 | cout<<"Info: total number of cells per multidim: "<<vct->GetNTotCells()<<endl; | |
1024 | ULong64_t glindex=vct->GetGlobalAddressFromValues(values,0); | |
1025 | ||
1026 | cout<<"The histogram you want is:\t"<<glindex+vct->GetNTotCells()*ptbin<<endl; | |
1027 | return glindex+vct->GetNTotCells()*ptbin; | |
1028 | } | |
1029 | ||
1030 | // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root | |
1031 | // ptbin = pt bin | |
1032 | // values = array of the cut values: the dimention can be <= number of variables maximized | |
1033 | // valsgiven = array of dimention = to the number of variables optimized. For each variable put kTRUE if the value is given (in values), kFALSE otherwise | |
1034 | // nhistinrange = pass an integer which will contains the number of histogram returned (that is the dimention of the Int_t* returned) | |
1035 | ||
1036 | //NB: Remember that the cut applied is the lower edge of the step where lower=looser | |
1037 | ||
1038 | Int_t* GetRangeHistFromValues(AliMultiDimVector* vct,Int_t ptbin,Bool_t* valsgiven,Float_t* values,Int_t& nhistinrange){ | |
1039 | Int_t nvargiven=0; | |
82487ae7 | 1040 | nhistinrange=1; |
cb3b714e | 1041 | |
1042 | Int_t nvar4opt=vct->GetNVariables(); | |
1043 | Float_t allvals[nvar4opt]; | |
1044 | ||
1045 | for (Int_t i=0;i<nvar4opt;i++) { | |
1046 | if(valsgiven[i]==kTRUE) { | |
1047 | allvals[i]=values[nvargiven]; | |
1048 | nvargiven++; | |
1049 | } | |
1050 | else { | |
1051 | nhistinrange+=vct->GetNCutSteps(i); | |
1052 | allvals[i]=vct->GetCutValue(i,0); | |
82487ae7 | 1053 | //allvals[i]=vct->GetCutValue(0,i); //ivar,icell |
cb3b714e | 1054 | } |
1055 | } | |
1056 | cout<<nhistinrange<<" index will be returned"<<endl; | |
cb3b714e | 1057 | Int_t *rangeofhistos=new Int_t[nhistinrange]; |
82487ae7 | 1058 | |
1059 | if(nhistinrange==1){ | |
1060 | rangeofhistos[0]=GetNHistFromValues(vct,ptbin,allvals); | |
1061 | cout<<"output"<<rangeofhistos[0]<<endl; | |
1062 | }else{ | |
1063 | Int_t index[nvar4opt-nvargiven]; | |
1064 | Int_t k=0; | |
1065 | for (Int_t i=0;i<nvar4opt;i++){ | |
1066 | if(valsgiven[i]==kFALSE) { | |
1067 | //cout<<"kTRUE==>"<<i<<endl; | |
1068 | index[k]=i; | |
1069 | k++; | |
1070 | } | |
1071 | } | |
1072 | ||
1073 | for(Int_t i=0;i<nvar4opt-nvargiven;i++){ //loop on number of free variables | |
1074 | cout<<"Info: incrementing "<<vct->GetAxisTitle(index[i])<<endl; | |
1075 | for(Int_t j=0;j<vct->GetNCutSteps(i);j++){ //loop on steps of each free variable | |
1076 | allvals[index[i]]=vct->GetCutValue(index[i],j); | |
1077 | rangeofhistos[i*vct->GetNCutSteps(i)+j]=GetNHistFromValues(vct,ptbin,allvals); | |
1078 | } | |
cb3b714e | 1079 | } |
1080 | } | |
1081 | return rangeofhistos; | |
1082 | } | |
1083 | ||
1084 | // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root | |
1085 | // ptbin = pt bin | |
1086 | // nhist = number of the histogram from which you want to have the cut values (returned) | |
1087 | ||
1088 | Float_t* GetCutValuesFromNHist(AliMultiDimVector* vct,Int_t ptbin,Int_t nhist){ | |
1089 | ULong64_t totCells=vct->GetNTotCells(); | |
1090 | ULong64_t globadd=nhist-ptbin*totCells; | |
1091 | const Int_t nvars=vct->GetNVariables(); | |
1092 | Float_t* cuts=new Float_t[nvars]; | |
1093 | Int_t ptinside; | |
1094 | vct->GetCutValuesFromGlobalAddress(globadd,cuts,ptinside); | |
1095 | return cuts; | |
1096 | } | |
1097 | ||
1098 | // ptbin = pt bin | |
1099 | // values = array of the cut values: the dimention can be <= number of variables maximized | |
1100 | // valsgiven = array of dimention = to the number of variables optimized. For each variable put kTRUE if the value is given (in values), kFALSE otherwise | |
1101 | // | |
1102 | ||
5b8639e7 | 1103 | void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString path="./",Int_t decCh=2){ |
cb3b714e | 1104 | gStyle->SetFrameBorderMode(0); |
1105 | gStyle->SetCanvasColor(0); | |
1106 | gStyle->SetFrameFillColor(0); | |
1107 | gStyle->SetOptStat(0); | |
1108 | ||
1109 | Int_t nhists; | |
1110 | TString filename="AnalysisResults.root"; | |
1111 | TString dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv"; | |
82487ae7 | 1112 | TString centr="020"; |
1113 | ||
cb3b714e | 1114 | TFile *fin=new TFile(Form("%s%s",path.Data(),filename.Data())); |
1115 | if(!fin){ | |
1116 | cout<<path.Data()<<filename.Data()<<" not found"<<endl; | |
1117 | return; | |
1118 | } | |
1119 | TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname); | |
1120 | if(!dir){ | |
1121 | cout<<"Directory "<<dirname<<" not found"<<endl; | |
1122 | return; | |
1123 | } | |
1124 | switch (decCh) { | |
1125 | case 0: | |
1126 | listname+="Dplus"; | |
1127 | mdvlistname+="Dplus"; | |
1128 | ||
1129 | break; | |
1130 | case 1: | |
1131 | listname+="D0"; | |
1132 | mdvlistname+="D0"; | |
1133 | ||
1134 | break; | |
1135 | case 2: | |
5b8639e7 | 1136 | listname+="Dstar0100"; |
1137 | mdvlistname+="Dstar0100"; | |
cb3b714e | 1138 | |
1139 | break; | |
1140 | case 3: | |
1141 | listname+="Ds"; | |
1142 | mdvlistname+="Ds"; | |
1143 | ||
1144 | break; | |
1145 | case 4: | |
1146 | listname+="D04"; | |
1147 | mdvlistname+="D04"; | |
1148 | ||
1149 | break; | |
1150 | case 5: | |
1151 | listname+="Lc"; | |
1152 | mdvlistname+="Lc"; | |
1153 | ||
1154 | break; | |
1155 | default: | |
1156 | cout<<decCh<<" is not allowed as decay channel "<<endl; | |
1157 | return; | |
1158 | } | |
82487ae7 | 1159 | mdvlistname+=centr; |
1160 | listname+=centr; | |
cb3b714e | 1161 | |
1162 | TList* listamdv= (TList*)dir->Get(mdvlistname); | |
1163 | if(!listamdv) { | |
1164 | cout<<mdvlistname<<" doesn't exist"<<endl; | |
1165 | return; | |
1166 | } | |
1167 | ||
1168 | AliMultiDimVector* vct=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",ptbin)); | |
1169 | ||
1170 | TFile* fin2; | |
1171 | TString filehistname=""; | |
1172 | Int_t* indexes=GetRangeHistFromValues(vct,ptbin,valsgiven,values,nhists); | |
1173 | for(Int_t i=0;i<nhists;i++){ | |
1174 | ||
1175 | fin2=new TFile(Form("%sh%dExpMassFit.root",path.Data(),indexes[i])); | |
1176 | if(!fin2){ | |
1177 | cout<<"File "<<indexes[i]<<" not found!"<<endl; | |
1178 | continue; | |
1179 | }else{ | |
1180 | TCanvas *cv=(TCanvas*)fin2->Get(Form("c1h%dExp",indexes[i])); | |
1181 | if(!cv){ | |
1182 | cout<<"Canvas c1h"<<indexes[i]<<"Exp not found among"; | |
1183 | fin2->ls(); | |
1184 | continue; | |
1185 | }else{ | |
1186 | cv->Draw(); | |
1187 | cv->SaveAs(Form("h%dExpMassFit.png",indexes[i])); | |
1188 | fin2=0x0; | |
1189 | } | |
1190 | } | |
1191 | } | |
1192 | } | |
82487ae7 | 1193 | |
5b8639e7 | 1194 | void Merge2Bins(Int_t b1, Int_t b2,TString pathin="./",Int_t decCh=2,TString part="both"/*"A" anti-particle, "P" particle*/){ |
82487ae7 | 1195 | |
1196 | if(b2!=b1+1){ | |
1197 | printf("The bins to be merget must be consecutive. Check! [b1 = %d, b2= %d]\n",b1,b2); | |
1198 | return; | |
1199 | } | |
1200 | ||
1201 | TFile *fin=new TFile(Form("%sAnalysisResults.root",pathin.Data())); | |
1202 | if (!fin){ | |
1203 | cout<<"Input file not found"<<endl; | |
1204 | return; | |
1205 | } | |
1206 | ||
1207 | TString dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv"; | |
1208 | ||
1209 | switch (decCh) { | |
1210 | case 0: | |
1211 | listname+="Dplus"; | |
1212 | mdvlistname+="Dplus"; | |
1213 | break; | |
1214 | case 1: | |
1215 | listname+="D0"; | |
1216 | mdvlistname+="D0"; | |
1217 | break; | |
1218 | case 2: | |
5b8639e7 | 1219 | listname+="Dstar0100"; |
1220 | mdvlistname+="Dstar0100"; | |
82487ae7 | 1221 | break; |
1222 | case 3: | |
1223 | listname+="Ds"; | |
1224 | mdvlistname+="Ds"; | |
1225 | break; | |
1226 | case 4: | |
1227 | listname+="D04"; | |
1228 | mdvlistname+="D04"; | |
1229 | break; | |
1230 | case 5: | |
1231 | listname+="Lc"; | |
1232 | mdvlistname+="Lc"; | |
1233 | break; | |
1234 | default: | |
1235 | cout<<decCh<<" is not allowed as decay channel "<<endl; | |
1236 | return; | |
1237 | } | |
1238 | ||
1239 | if(part!="both"){ | |
1240 | listname.Append(part); | |
1241 | mdvlistname.Append(part); | |
1242 | } | |
1243 | ||
1244 | TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname); | |
1245 | if(!dir){ | |
1246 | cout<<"Directory "<<dirname<<" not found"<<endl; | |
1247 | return; | |
1248 | } | |
1249 | ||
1250 | TList* histlist= (TList*)dir->Get(listname); | |
1251 | if(!histlist) { | |
1252 | cout<<listname<<" doesn't exist"<<endl; | |
1253 | return; | |
1254 | } | |
1255 | ||
1256 | TList* listamdv= (TList*)dir->Get(mdvlistname); | |
1257 | if(!listamdv) { | |
1258 | cout<<mdvlistname<<" doesn't exist"<<endl; | |
1259 | return; | |
1260 | } | |
1261 | if (!gSystem->AccessPathName(Form("merged%d%d",b1,b2))) gSystem->Exec(Form("mkdir merged%d%d",b1,b2)); | |
1262 | gSystem->Exec(Form("cd merged%d%d",b1,b2)); | |
1263 | ||
1264 | TFile* fout=new TFile("mergeAnalysisResults.root","recreate"); | |
1265 | ||
1266 | fout->mkdir(dirname); | |
1267 | TList* listmdvout=new TList(); | |
1268 | listmdvout->SetName(listamdv->GetName()); | |
1269 | listmdvout->SetOwner(); | |
1270 | //listmdvout->SetTitle(listamdv->GetTitle()); | |
1271 | TList* histlistout=new TList(); | |
1272 | histlistout->SetName(histlist->GetName()); | |
1273 | histlistout->SetOwner(); | |
1274 | //histlistout->SetTitle(histlist->GetTitle()); | |
1275 | ||
1276 | AliMultiDimVector* mdvin1=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",b1)); | |
1277 | AliMultiDimVector* mdvin2=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",b2)); | |
1278 | ||
1279 | Int_t ntotHperbin = mdvin1->GetNTotCells(); | |
1280 | if(mdvin2->GetNTotCells() != ntotHperbin) { | |
1281 | cout<<"Error! Number of histos in pt bin "<<b1<<" = "<<ntotHperbin<<" != Number of histos in pt bin "<<b2<<" = "<<mdvin2->GetNTotCells()<<endl; | |
1282 | return; | |
1283 | } | |
1284 | Int_t nvar1=mdvin1->GetNVariables(); | |
1285 | if(nvar1 != mdvin2->GetNVariables()){ | |
1286 | cout<<"Error! Mismatch in number of variables"<<endl; | |
1287 | return; | |
1288 | } | |
1289 | ||
1290 | Float_t newptbins[2]={mdvin1->GetPtLimit(0),mdvin2->GetPtLimit(1)}; | |
1291 | Float_t loosercuts[nvar1], tightercuts[nvar1]; | |
1292 | TString axistitles[nvar1]; | |
1293 | Int_t ncells[nvar1]; | |
1294 | ||
1295 | for (Int_t ivar=0;ivar<nvar1;ivar++){ | |
1296 | loosercuts[ivar]=mdvin1->GetCutValue(ivar,0); | |
1297 | if(loosercuts[ivar] - mdvin2->GetCutValue(ivar,0) < 1e-8) printf("Warning! The loose cut %s is different between the 2: %f and %f\n",mdvin1->GetAxisTitle(ivar).Data(),loosercuts[ivar],mdvin2->GetCutValue(ivar,0)); | |
1298 | tightercuts[ivar]=mdvin1->GetCutValue(ivar,mdvin1->GetNCutSteps(ivar)-1); | |
1299 | if(tightercuts[ivar] - mdvin2->GetCutValue(ivar,mdvin1->GetNCutSteps(ivar)-1) < 1e-8) printf("Warning! The tight cut %s is different between the 2: %f and %f\n",mdvin1->GetAxisTitle(ivar).Data(),tightercuts[ivar],mdvin2->GetCutValue(ivar,mdvin2->GetNCutSteps(ivar))); | |
1300 | axistitles[ivar]=mdvin1->GetAxisTitle(ivar); | |
1301 | cout<<axistitles[ivar]<<"\t"; | |
1302 | ncells[ivar]=mdvin1->GetNCutSteps(ivar); | |
1303 | } | |
1304 | ||
1305 | AliMultiDimVector* mdvout= new AliMultiDimVector(Form("multiDimVectorPtBin%d",b1),"MultiDimVector",1,newptbins,mdvin1->GetNVariables(),ncells,loosercuts,tightercuts,axistitles); | |
1306 | cout<<"Info: writing mdv"<<endl; | |
1307 | listmdvout->Add(mdvout); | |
1308 | ||
1309 | Bool_t isMC=kFALSE; | |
1310 | TH1F* htestIsMC=(TH1F*)histlist->FindObject("hSgn_0"); | |
1311 | if(htestIsMC) isMC=kTRUE; | |
1312 | Int_t nptbins=listamdv->GetEntries(); | |
1313 | Int_t nhist=(histlist->GetEntries()-1);//-1 because of fHistNevents | |
1314 | if(isMC) nhist/=4; ///4 because hMass_, hSgn_,hBkg_,hRfl_ | |
1315 | ||
1316 | cout<<"Merging bin from "<<mdvin1->GetPtLimit(0)<<" to "<<mdvin1->GetPtLimit(1)<<" and from "<<mdvin2->GetPtLimit(0)<<" to "<<mdvin2->GetPtLimit(1)<<endl; | |
1317 | Int_t firsth1=b1*ntotHperbin,firsth2=b2*ntotHperbin; //firsth2 = (b1+1)*ntotHperbin | |
1318 | Int_t lasth1=firsth1+ntotHperbin-1,lasth2=firsth2+ntotHperbin-1; | |
1319 | cout<<"Histograms from "<<firsth1<<" to "<<lasth1<<" and "<<firsth2<<" to "<<lasth2<<endl; | |
1320 | ||
1321 | //add the others mdv to the list | |
1322 | Int_t cnt=0; | |
1323 | for(Int_t i=0;i<nptbins;i++){ | |
1324 | if(i!=b1 && i!=b2){ | |
1325 | AliMultiDimVector* vcttmp=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",i)); | |
1326 | if(i>b2) { | |
1327 | vcttmp->SetName(Form("multiDimVectorPtBin%d",b2+cnt)); | |
1328 | cnt++; | |
1329 | } | |
1330 | listmdvout->Add(vcttmp); | |
1331 | } | |
1332 | } | |
1333 | ||
1334 | histlistout->Add((TH1F*)histlist->FindObject("fHistNEvents")); | |
1335 | ||
1336 | Int_t ih2=firsth2; | |
1337 | ||
1338 | for(Int_t ih1=firsth1;ih1<lasth1;ih1++){ | |
1339 | TH1F* h1=(TH1F*)histlist->FindObject(Form("hMass_%d",ih1)); | |
1340 | if(!h1){ | |
1341 | cout<<"hMass_"<<ih1<<" not found!"<<endl; | |
1342 | continue; | |
1343 | } | |
1344 | TH1F* h2=(TH1F*)histlist->FindObject(Form("hMass_%d",ih2)); | |
1345 | if(!h2){ | |
1346 | cout<<"hMass_"<<ih2<<" not found!"<<endl; | |
1347 | continue; | |
1348 | } | |
1349 | //h1->SetName(Form("hMass_%d",cnt)); | |
1350 | h1->Add(h2); | |
1351 | histlistout->Add(h1); | |
1352 | ih2++; | |
1353 | //cnt++; | |
1354 | h1=0x0; | |
1355 | } | |
1356 | ||
1357 | cnt=0; | |
1358 | for(Int_t j=0;j<ntotHperbin*nptbins;j++){ | |
1359 | if(!(j>=firsth1 && j<lasth2)){ | |
1360 | TH1F* htmp=(TH1F*)histlist->FindObject(Form("hMass_%d",j)); | |
1361 | if(j>=lasth2){ | |
1362 | //cout<<lasth1<<" + "<<cnt<<endl; | |
1363 | htmp->SetName(Form("hMass_%d",lasth1+cnt)); | |
1364 | cnt++; | |
1365 | } | |
1366 | histlistout->Add(htmp); | |
1367 | } | |
1368 | } | |
1369 | ||
1370 | fout->cd(); | |
1371 | ((TDirectoryFile*)fout->Get(dirname))->cd(); | |
1372 | listmdvout->Write(mdvlistname.Data(),TObject::kSingleKey); | |
1373 | histlistout->Write(listname.Data(),TObject::kSingleKey); | |
1374 | fout->Close(); | |
1375 | } | |
1376 | ||
1377 | void SubtractBkg(Int_t nhisto){ | |
1378 | ||
1379 | gStyle->SetFrameBorderMode(0); | |
1380 | gStyle->SetCanvasColor(0); | |
1381 | gStyle->SetFrameFillColor(0); | |
1382 | gStyle->SetOptStat(0); | |
1383 | ||
1384 | TString fitType="Exp"; | |
1385 | TString filename=Form("h%d%sMassFit.root",nhisto,fitType.Data()); | |
1386 | ||
1387 | TFile* fin=new TFile(filename.Data()); | |
1388 | if(!fin->IsOpen()){ | |
1389 | cout<<filename.Data()<<" not found, exit"<<endl; | |
1390 | return; | |
1391 | } | |
1392 | ||
1393 | TKey* key=((TKey*)((TList*)fin->GetListOfKeys())->At(fin->GetNkeys()-1)); | |
1394 | TCanvas* canvas=((TCanvas*)fin->Get(key->GetName())); | |
1395 | if(!canvas){ | |
1396 | cout<<"Canvas not found"<<endl; | |
1397 | return; | |
1398 | } | |
1399 | canvas->Draw(); | |
1400 | ||
1401 | TH1F* hfit=(TH1F*)canvas->FindObject("fhistoInvMass"); | |
1402 | if(!hfit){ | |
1403 | canvas->ls(); | |
1404 | cout<<"Histogram not found"<<endl; | |
1405 | return; | |
1406 | } | |
1407 | ||
1408 | TF1* funcBkgRecalc=(TF1*)hfit->FindObject("funcbkgRecalc"); | |
1409 | if(!funcBkgRecalc){ | |
1410 | cout<<"Background fit function (final) not found"<<endl; | |
1411 | return; | |
1412 | } | |
1413 | ||
1414 | TF1* funcBkgFullRange=(TF1*)hfit->FindObject("funcbkgFullRange"); | |
1415 | if(!funcBkgFullRange){ | |
1416 | cout<<"Background fit function (side bands) not found"<<endl; | |
1417 | return; | |
1418 | } | |
1419 | ||
1420 | Int_t nbins=hfit->GetNbinsX(); | |
1421 | Double_t min=hfit->GetBinLowEdge(1), width=hfit->GetBinWidth(1); | |
1422 | TH1F* hsubRecalc=(TH1F*)hfit->Clone("hsub"); | |
1423 | hsubRecalc->SetMarkerColor(kRed); | |
1424 | hsubRecalc->SetLineColor(kRed); | |
1425 | hsubRecalc->GetListOfFunctions()->Delete(); | |
1426 | TH1F* hsubFullRange=(TH1F*)hfit->Clone("hsub"); | |
1427 | hsubFullRange->SetMarkerColor(kGray+2); | |
1428 | hsubFullRange->SetLineColor(kGray+2); | |
1429 | hsubFullRange->GetListOfFunctions()->Delete(); | |
1430 | for(Int_t i=0;i<nbins;i++){ | |
1431 | //Double_t x=min+i*0.5*width; | |
1432 | Double_t x1=min+i*width, x2=min+(i+1)*width; | |
1433 | Double_t ycont=hfit->GetBinContent(i+1); | |
1434 | Double_t y1=funcBkgRecalc->Integral(x1,x2) / width;//funcBkgRecalc->Eval(x); | |
1435 | hsubRecalc->SetBinContent(i+1,ycont-y1); | |
1436 | Double_t y2=funcBkgFullRange->Integral(x1,x2) / width;//funcBkgFullRange->Eval(x); | |
1437 | hsubFullRange->SetBinContent(i+1,ycont-y2); | |
1438 | } | |
1439 | ||
1440 | TCanvas* c=new TCanvas("c","subtraction"); | |
1441 | c->cd(); | |
1442 | hsubRecalc->DrawClone(); | |
1443 | hsubFullRange->DrawClone("sames"); | |
1444 | ||
1445 | for(Int_t i=0;i<nbins;i++){ | |
1446 | if(hsubRecalc->GetBinContent(i+1)<0) hsubRecalc->SetBinContent(i+1,0); | |
1447 | if(hsubFullRange->GetBinContent(i+1)<0) hsubFullRange->SetBinContent(i+1,0); | |
1448 | } | |
1449 | ||
1450 | TCanvas *cvnewfits=new TCanvas("cvnewfits", "new Fits",1200,600); | |
1451 | cvnewfits->Divide(2,1); | |
1452 | ||
1453 | AliHFMassFitter fitter1(hsubRecalc,min,min+nbins*width,1,1); | |
1454 | fitter1.MassFitter(kFALSE); | |
1455 | fitter1.DrawHere(cvnewfits->cd(1)); | |
1456 | ||
1457 | AliHFMassFitter fitter2(hsubFullRange,min,min+nbins*width,1,1); | |
1458 | fitter2.MassFitter(kFALSE); | |
1459 | fitter2.DrawHere(cvnewfits->cd(2)); | |
1460 | ||
1461 | canvas->SaveAs(Form("h%d%sMassFit.png",nhisto,fitType.Data())); | |
1462 | c->SaveAs(Form("h%d%sSubtr.png",nhisto,fitType.Data())); | |
1463 | cvnewfits->SaveAs(Form("h%d%sFitNew.png",nhisto,fitType.Data())); | |
1464 | } | |
5b8639e7 | 1465 |