]>
Commit | Line | Data |
---|---|---|
3ed21f0a | 1 | #include <fstream> |
2 | #include <Riostream.h> | |
3 | #include <TH1F.h> | |
4 | #include <TF1.h> | |
5 | #include <TFile.h> | |
6 | #include <TCanvas.h> | |
7 | #include <TClonesArray.h> | |
8 | #include <TStyle.h> | |
9 | #include <TLegend.h> | |
10 | #include <TGraphErrors.h> | |
11 | #include <TGraph.h> | |
12 | #include <TMultiGraph.h> | |
13 | #include <TObjectTable.h> | |
14 | #include <TDatabasePDG.h> | |
15 | ||
16 | #include <AliMultiDimVector.h> | |
c24d2d9f | 17 | #include "AliHFMassFitter.h" |
3ed21f0a | 18 | #include <AliSignificanceCalculator.h> |
19 | ||
20 | #include <fstream> | |
21 | ||
22 | //decCh: | |
23 | //- 0 = kDplustoKpipi | |
24 | //- 1 = kD0toKpi | |
25 | //- 2 = kDstartoKpipi | |
26 | //- 3 = kDstoKKpi | |
27 | //- 4 = kD0toKpipipi | |
28 | //- 5 = kLambdactopKpi | |
29 | ||
4bcb0ebb | 30 | //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) |
31 | ||
32 | Bool_t charmCutsOptimization(Double_t nsigma=2,Int_t decCh=1,Int_t fitbtype=0,Int_t rebin=2,Bool_t writefit=kFALSE,Double_t sigma=0.012,Int_t minentries=50,Double_t *rangefit=0x0,TString hname="hMass_"){ | |
3ed21f0a | 33 | |
34 | TString filename="AnalysisResults.root",dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv"; | |
35 | ||
36 | ||
37 | Int_t pdg; | |
38 | Double_t mass; | |
39 | ||
40 | switch (decCh) { | |
41 | case 0: | |
42 | listname+="Dplus"; | |
43 | mdvlistname+="Dplus"; | |
44 | pdg=411; | |
45 | break; | |
46 | case 1: | |
47 | listname+="D0"; | |
48 | mdvlistname+="D0"; | |
49 | pdg=421; | |
50 | break; | |
51 | case 2: | |
52 | listname+="Dstar"; | |
53 | mdvlistname+="Dstar"; | |
54 | pdg=413; | |
55 | break; | |
56 | case 3: | |
57 | listname+="Ds"; | |
58 | mdvlistname+="Ds"; | |
59 | pdg=431; | |
60 | break; | |
61 | case 4: | |
62 | listname+="D04"; | |
63 | mdvlistname+="D04"; | |
64 | pdg=421; | |
65 | break; | |
66 | case 5: | |
67 | listname+="Lc"; | |
68 | mdvlistname+="Lc"; | |
69 | pdg=4122; | |
70 | break; | |
71 | default: | |
72 | cout<<decCh<<" is not allowed as decay channel "<<endl; | |
73 | return kFALSE; | |
74 | } | |
75 | mass=TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
76 | ||
77 | cout<<"Mass = "<<mass<<endl; | |
78 | ||
79 | gStyle->SetCanvasColor(0); | |
80 | gStyle->SetFrameFillColor(0); | |
81 | gStyle->SetTitleFillColor(0); | |
82 | gStyle->SetStatColor(0); | |
83 | ||
84 | Int_t countFitFail=0,countSgnfFail=0,countNoHist=0,countBkgOnly=0; | |
85 | ofstream outcheck("output.dat"); | |
c24d2d9f | 86 | ofstream outdetail("discarddetails.dat"); |
3ed21f0a | 87 | |
88 | outcheck<<"ptbin\tmdvGlobAddr\thistIndex\tSignif\tS\tB"<<endl; | |
c24d2d9f | 89 | outdetail<<"ptbin\tmdvGlobAddr\thistIndex\trelErrS\t\tmean_F-mass"<<endl; |
3ed21f0a | 90 | TFile *fin=new TFile(filename.Data()); |
91 | if(!fin->IsOpen()){ | |
92 | cout<<"File "<<filename.Data()<<" not found"<<endl; | |
93 | return kFALSE; | |
94 | } | |
95 | ||
96 | TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname); | |
97 | if(!dir){ | |
98 | cout<<"Directory "<<dirname<<" not found"<<endl; | |
99 | return kFALSE; | |
100 | } | |
101 | ||
102 | TList* histlist= (TList*)dir->Get(listname); | |
103 | if(!histlist) { | |
104 | cout<<listname<<" doesn't exist"<<endl; | |
105 | return kFALSE; | |
106 | } | |
107 | ||
108 | TList* listamdv= (TList*)dir->Get(mdvlistname); | |
109 | if(!listamdv) { | |
110 | cout<<mdvlistname<<" doesn't exist"<<endl; | |
111 | return kFALSE; | |
112 | } | |
113 | ||
6b64765b | 114 | TH1F* hstat=(TH1F*)histlist->FindObject("fHistNEvents"); |
115 | TCanvas *cst=new TCanvas("hstat","Summary of statistics"); | |
116 | if(hstat) { | |
117 | cst->cd(); | |
118 | cst->SetGrid(); | |
119 | hstat->Draw("htext0"); | |
120 | hstat->SaveAs("hstat.png"); | |
121 | }else{ | |
122 | cout<<"Warning! fHistNEvents not found in "<<listname.Data()<<endl; | |
123 | } | |
124 | ||
125 | Bool_t isMC=kFALSE; | |
126 | TH1F* htestIsMC=(TH1F*)histlist->FindObject("hSgn_0"); | |
127 | if(htestIsMC) isMC=kTRUE; | |
128 | ||
3ed21f0a | 129 | Int_t nptbins=listamdv->GetEntries(); |
6b64765b | 130 | Int_t nhist=(histlist->GetEntries()-1);//-1 because of fHistNevents |
131 | if(isMC) nhist/=4; ///4 because hMass_, hSgn_,hBkg_,hRfl_ | |
3ed21f0a | 132 | Int_t count=0; |
133 | Int_t *indexes= new Int_t[nhist]; | |
134 | //initialize indexes[i] to -1 | |
135 | for(Int_t i=0;i<nhist;i++){ | |
136 | indexes[i]=-1; | |
137 | } | |
138 | ||
139 | TFile* fout=new TFile(Form("outputSignifMaxim.root"),"recreate"); | |
140 | ||
141 | //Check wheter histograms are filled | |
142 | for(Int_t i=0;i<nhist;i++){ | |
143 | TString name=Form("%s%d",hname.Data(),i); | |
144 | TH1F* h=(TH1F*)histlist->FindObject(name.Data()); | |
145 | ||
146 | if(!h){ | |
147 | cout<<name<<" not found"<<endl; | |
148 | continue; | |
149 | } | |
150 | ||
151 | if(h->GetEntries()>minentries){ | |
152 | //cout<<"Entries = "<<h->GetEntries()<<endl; | |
153 | if (h->Integral() > minentries){ | |
154 | cout<<i<<") Integral = "<<h->Integral()<<endl; | |
b7194c48 | 155 | indexes[i]=i; |
3ed21f0a | 156 | count++; |
157 | } | |
158 | } | |
159 | } | |
160 | ||
161 | cout<<"There are "<<count<<" histogram with more than "<<minentries<<" entries"<<endl; | |
162 | if(count==0) { | |
163 | cout<<"No histogram to draw..."<<endl; | |
164 | return kFALSE; | |
165 | } | |
166 | ||
167 | //create multidimvectors | |
168 | ||
169 | //for(Int_t i=0;i<1;i++){ | |
170 | for(Int_t i=0;i<nptbins;i++){ | |
171 | ||
172 | //multidimvectors for signal | |
173 | AliMultiDimVector *mdvS=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",i)); | |
174 | TString name=mdvS->GetName(),nameErr="err",setname=""; | |
175 | ||
176 | setname=Form("S%s",name.Data()); | |
177 | mdvS->SetName(setname.Data()); | |
4bcb0ebb | 178 | outcheck<<"\n"<<mdvS->GetPtLimit(0)<<" < Pt <"<<mdvS->GetPtLimit(1)<<endl; |
3ed21f0a | 179 | AliMultiDimVector *mdvSerr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); |
180 | setname=Form("%sS%s",nameErr.Data(),name.Data()); | |
181 | mdvSerr->SetName(setname.Data()); | |
182 | ||
183 | //multidimvectors for background | |
184 | setname=Form("B%s",name.Data()); | |
185 | AliMultiDimVector *mdvB=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
186 | ||
187 | AliMultiDimVector *mdvBerr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
188 | setname=Form("%sB%s",nameErr.Data(),name.Data()); | |
189 | mdvBerr->SetName(setname.Data()); | |
190 | ||
191 | //multidimvectors for significance | |
192 | setname=Form("Sgf%s",name.Data()); | |
193 | AliMultiDimVector *mdvSgnf=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
194 | ||
195 | AliMultiDimVector *mdvSgnferr=(AliMultiDimVector*)mdvS->Clone(setname.Data()); | |
196 | setname=Form("%sSgf%s",nameErr.Data(),name.Data()); | |
197 | mdvSgnferr->SetName(setname.Data()); | |
198 | ||
199 | Int_t nhistforptbin=mdvS->GetNTotCells(); | |
200 | //Int_t nvarsopt=mdvS->GetNVariables(); | |
201 | ||
202 | cout<<"nhistforptbin = "<<nhistforptbin<<endl; | |
203 | ||
204 | //loop on all histograms and do AliHFMassFitter | |
205 | //for(Int_t ih=0;ih<1;ih++){ | |
206 | for(Int_t ih=0;ih<nhistforptbin;ih++){ | |
207 | printf("Analyzing indexes[%d] = %d \n",ih+i*nhistforptbin,indexes[ih+i*nhistforptbin]); | |
208 | ||
209 | if(indexes[ih+i*nhistforptbin] != -1){ | |
210 | TString name=Form("%s%d",hname.Data(),indexes[ih+i*nhistforptbin]); | |
211 | TH1F* h=(TH1F*)histlist->FindObject(name.Data()); | |
212 | ||
213 | Int_t nbin=((TH1F*)histlist->FindObject(name))->GetNbinsX(); | |
214 | Double_t min=((TH1F*)histlist->FindObject(name))->GetBinLowEdge(7); | |
215 | Double_t max=((TH1F*)histlist->FindObject(name))->GetBinLowEdge(nbin-5)+((TH1F*)histlist->FindObject(name))->GetBinWidth(nbin-5); | |
216 | if(rangefit) { | |
217 | min=rangefit[0]; | |
218 | max=rangefit[1]; | |
219 | } | |
220 | ||
c24d2d9f | 221 | AliHFMassFitter fitter(h,min, max,rebin,fitbtype); |
3ed21f0a | 222 | fitter.SetInitialGaussianMean(mass); |
c24d2d9f | 223 | fitter.SetInitialGaussianSigma(sigma); |
3ed21f0a | 224 | |
225 | if(ih==0) fitter.InitNtuParam(Form("ntuPtbin%d",i)); | |
226 | // fitter.SetHisto(h); | |
227 | // fitter.SetRangeFit(min,max); | |
228 | //fitter.SetRangeFit(1.68,2.05); | |
229 | ||
230 | //fitter.SetType(fitbtype,0); | |
231 | ||
232 | Bool_t ok=fitter.MassFitter(kFALSE); | |
c24d2d9f | 233 | if(!ok){ |
234 | cout<<"FIT NOT OK!"<<endl; | |
235 | // ok=fitter.RefitWithBkgOnly(kFALSE); | |
236 | // if (ok){ //onlybkg | |
237 | countBkgOnly++; | |
238 | // Double_t bkg=0,errbkg=0.; | |
239 | // fitter.Background(nsigma,bkg,errbkg); | |
240 | outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t 0\t xxx"<<"\t bkgonly"<<endl; | |
241 | // mdvSgnf->SetElement(ih,0); | |
242 | // mdvSgnferr->SetElement(ih,0); | |
243 | // mdvS->SetElement(ih,0); | |
244 | // mdvSerr->SetElement(ih,0); | |
245 | // mdvB->SetElement(ih,bkg); | |
246 | // mdvBerr->SetElement(ih,errbkg); | |
247 | // }else{ //bkg fit failed | |
248 | // cout<<"Setting to 0"<<endl; | |
249 | // outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t fit failed also with only bkg"<<endl; | |
250 | // countFitFail++; | |
251 | // mdvSgnf->SetElement(ih,0); | |
252 | // mdvSgnferr->SetElement(ih,0); | |
253 | // mdvS->SetElement(ih,0); | |
254 | // mdvSerr->SetElement(ih,0); | |
255 | // mdvB->SetElement(ih,0); | |
256 | // mdvBerr->SetElement(ih,0); | |
257 | // } | |
3ed21f0a | 258 | }else{ //fit ok! |
259 | ||
4bcb0ebb | 260 | if(writefit) fitter.WriteCanvas(Form("h%d",ih+i*nhistforptbin),"./",nsigma); |
3ed21f0a | 261 | Double_t signif=0, signal=0, background=0, errSignif=0, errSignal=0, errBackground=0; |
262 | fitter.Signal(nsigma,signal,errSignal); | |
263 | fitter.Background(nsigma,background,errBackground); | |
264 | Double_t meanfit=fitter.GetMean(); | |
265 | Double_t sigmafit=fitter.GetSigma(); | |
266 | ||
267 | if(sigmafit > 0.03){ | |
268 | //refit | |
269 | fitter.Reset(); | |
270 | fitter.SetHisto(h); | |
c24d2d9f | 271 | fitter.SetRangeFit(1.8,1.93); //WARNING!! this is candidate dependant!! (change) |
272 | fitter.SetInitialGaussianMean(mass); | |
273 | fitter.SetInitialGaussianSigma(sigma); | |
3ed21f0a | 274 | ok=fitter.MassFitter(kFALSE); |
275 | if(ok){ | |
276 | meanfit=fitter.GetMean(); | |
277 | sigmafit=fitter.GetSigma(); | |
278 | fitter.Signal(nsigma,signal,errSignal); | |
279 | fitter.Background(nsigma,background,errBackground); | |
280 | } | |
281 | } //sigma check done and fit recalc | |
282 | ||
283 | if(ok==kTRUE && sigmafit < 0.03 && signal > 0 && background > 0){ | |
284 | fitter.Significance(nsigma,signif,errSignif); | |
285 | if(signif >0){ | |
c24d2d9f | 286 | if(errSignal/signal < 0.35 && TMath::Abs(meanfit-mass)<0.01){ |
3ed21f0a | 287 | outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t"<<signif<<" +- "<<errSignif<<"\t"<<signal<<" +- "<<errSignal<<"\t"<<background<<" +- "<<errBackground<<endl; |
288 | mdvSgnf->SetElement(ih,signif); | |
289 | mdvSgnferr->SetElement(ih,errSignif); | |
290 | mdvS->SetElement(ih,signal); | |
291 | mdvSerr->SetElement(ih,errSignal); | |
292 | mdvB->SetElement(ih,background); | |
293 | mdvBerr->SetElement(ih,errBackground); | |
294 | ||
295 | }else{ | |
296 | ok=fitter.RefitWithBkgOnly(kFALSE); | |
297 | if (ok){ | |
298 | countBkgOnly++; | |
299 | Double_t bkg=0,errbkg=0.; | |
c24d2d9f | 300 | Double_t nsigmarange[2]={mass-nsigma*sigma,mass+nsigma*sigma}; |
301 | fitter.Background(nsigmarange[0],nsigmarange[1],bkg,errbkg); | |
3ed21f0a | 302 | outcheck<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t 0\t "<<bkg <<"\t bkgonly"<<endl; |
c24d2d9f | 303 | outdetail<<i<<"\t\t "<<ih<<"\t\t"<<indexes[ih+i*nhistforptbin]<<"\t"<<errSignal/signal<<"\t\t "<<meanfit-mass<<endl; |
3ed21f0a | 304 | mdvSgnf->SetElement(ih,0); |
305 | mdvSgnferr->SetElement(ih,0); | |
306 | mdvS->SetElement(ih,0); | |
307 | mdvSerr->SetElement(ih,0); | |
308 | mdvB->SetElement(ih,bkg); | |
309 | mdvBerr->SetElement(ih,errbkg); | |
310 | } | |
311 | }//only bkg | |
312 | }//check signif>0 | |
313 | else{ | |
314 | countSgnfFail++; | |
315 | cout<<"Setting to 0 (fitter results meaningless)"<<endl; | |
316 | outcheck<<"\t S || B || sgnf negative"; | |
317 | outcheck<<endl; | |
318 | mdvSgnf->SetElement(ih,0); | |
319 | mdvSgnferr->SetElement(ih,0); | |
320 | mdvS->SetElement(ih,0); | |
321 | mdvSerr->SetElement(ih,0); | |
322 | mdvB->SetElement(ih,0); | |
323 | mdvBerr->SetElement(ih,0); | |
324 | } | |
325 | } //end fit ok! | |
326 | } | |
327 | }else{ //check on histo | |
328 | ||
329 | countNoHist++; | |
330 | cout<<"Setting to 0 (indexes = -1)"<<endl; | |
331 | outcheck<<"\t histo not accepted for fit"; | |
332 | outcheck<<endl; | |
333 | mdvSgnf->SetElement(ih,0); | |
334 | mdvSgnferr->SetElement(ih,0); | |
335 | mdvS->SetElement(ih,0); | |
336 | mdvSerr->SetElement(ih,0); | |
337 | mdvB->SetElement(ih,0); | |
338 | mdvBerr->SetElement(ih,0); | |
339 | ||
340 | } | |
341 | ||
342 | ||
343 | cout<<mdvS->GetElement(ih)<<"\t"<<mdvB->GetElement(ih)<<endl; | |
344 | ||
345 | } | |
346 | ||
347 | fout->cd(); | |
348 | mdvS->Write(); | |
349 | mdvB->Write(); | |
350 | mdvSgnf->Write(); | |
351 | ||
352 | mdvSerr->Write(); | |
353 | mdvBerr->Write(); | |
354 | mdvSgnferr->Write(); | |
355 | ||
356 | } | |
357 | ||
358 | ||
359 | fout->Close(); | |
360 | ||
6b64765b | 361 | 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; |
3ed21f0a | 362 | outcheck.close(); |
363 | return kTRUE; | |
364 | } | |
365 | ||
366 | ||
367 | //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: | |
368 | ||
369 | // which=0 plot significance | |
370 | // =1 plot signal | |
371 | // =2 plot background | |
c24d2d9f | 372 | // =3 plot S/B |
6b64765b | 373 | // 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 |
374 | // 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 | 375 | |
6b64765b | 376 | |
377 | void showMultiDimVector(Int_t n=2,Int_t which=0, Bool_t maximize=kTRUE,Bool_t readfromfile=kFALSE){ | |
3ed21f0a | 378 | |
379 | gStyle->SetCanvasColor(0); | |
380 | gStyle->SetFrameFillColor(0); | |
381 | gStyle->SetTitleFillColor(0); | |
382 | gStyle->SetOptStat(0); | |
383 | //gStyle->SetOptTitle(0); | |
384 | ||
6b64765b | 385 | if((maximize && readfromfile) || (!maximize && !readfromfile)){ |
386 | cout<<"Error! maximize & readfromfile cannot be both kTRUE or kFALSE"<<endl; | |
387 | return; | |
388 | } | |
389 | ||
3ed21f0a | 390 | TFile* fin=new TFile("outputSignifMaxim.root"); |
391 | if(!fin->IsOpen()){ | |
392 | cout<<"outputSignifMaxim.root not found"<<endl; | |
393 | return; | |
394 | } | |
395 | ||
396 | if(n>2){ | |
397 | cout<<"Error! cannot show "<<n+1<<" dimentions"<<endl; | |
398 | return; | |
399 | } | |
400 | ||
c24d2d9f | 401 | TString name,title,namebis,shorttitle; |
3ed21f0a | 402 | switch (which){ |
403 | case 0: | |
404 | name="SgfmultiDimVectorPtBin"; | |
405 | title="Significance"; | |
c24d2d9f | 406 | shorttitle="Sgnf"; |
3ed21f0a | 407 | break; |
408 | case 1: | |
409 | name="SmultiDimVectorPtBin"; | |
410 | title="Signal"; | |
c24d2d9f | 411 | shorttitle="S"; |
3ed21f0a | 412 | break; |
413 | case 2: | |
414 | name="BmultiDimVectorPtBin"; | |
415 | title="Background"; | |
c24d2d9f | 416 | shorttitle="B"; |
3ed21f0a | 417 | break; |
418 | case 3: | |
c24d2d9f | 419 | name="SmultiDimVectorPtBin"; |
420 | namebis="BmultiDimVectorPtBin"; | |
421 | title="Signal over Background "; | |
422 | shorttitle="SoB"; | |
3ed21f0a | 423 | break; |
c24d2d9f | 424 | // case 4: |
425 | // name="errBmultiDimVectorPtBin"; | |
426 | // title="Background (error)"; | |
427 | // break; | |
3ed21f0a | 428 | } |
429 | ||
4bcb0ebb | 430 | Int_t nptbins=25; |
3ed21f0a | 431 | |
4bcb0ebb | 432 | for(Int_t ip=0;ip<=nptbins;ip++){ |
3ed21f0a | 433 | TString mdvname=Form("%s%d",name.Data(),ip); |
434 | AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname); | |
435 | if(!mdv){ | |
436 | nptbins=ip; | |
437 | cout<<"Number of pt bins "<<ip<<endl; | |
438 | break; | |
439 | } | |
440 | } | |
441 | ||
442 | cout<<"Projecting "<<title.Data()<<" with respect to the maximization variable(s) [chose]"<<endl; | |
443 | ||
444 | Int_t variable[2]; //no more than 2D | |
c24d2d9f | 445 | TString mdvname=Form("%s0",name.Data()), mdverrname="";//, mdvnamebis="", mdverrnamebis=""; |
3ed21f0a | 446 | AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname); |
447 | AliMultiDimVector* mdverr=0x0; | |
448 | if(!mdv){ | |
449 | cout<<mdvname.Data()<<" not found"<<endl; | |
450 | return; | |
451 | } | |
c24d2d9f | 452 | |
3ed21f0a | 453 | |
454 | Int_t nvarsopt=mdv->GetNVariables(); | |
6b64765b | 455 | //Int_t nfixed=nvarsopt-n; |
456 | Int_t fixedvars[nvarsopt]; | |
457 | Int_t allfixedvars[nvarsopt*nptbins]; | |
3ed21f0a | 458 | |
6b64765b | 459 | fstream writefixedvars; |
460 | if(readfromfile) { | |
461 | //open file in read mode | |
462 | writefixedvars.open("fixedvars.dat",ios::in); | |
463 | Int_t longi=0; | |
464 | while(writefixedvars){ | |
465 | writefixedvars>>allfixedvars[longi]; | |
466 | longi++; | |
467 | } | |
468 | } | |
469 | else { | |
470 | //open file in write mode | |
471 | writefixedvars.open("fixedvars.dat",ios::out); | |
3ed21f0a | 472 | } |
473 | ||
474 | //ask variables for projection | |
475 | for(Int_t k=0;k<nvarsopt;k++){ | |
3ed21f0a | 476 | cout<<k<<" "<<mdv->GetAxisTitle(k)<<endl; |
477 | } | |
478 | cout<<"Choose "<<n<<" variable(s)"<<endl; | |
479 | for(Int_t j=0;j<n;j++){ | |
480 | cout<<"var"<<j<<": "; | |
481 | cin>>variable[j]; | |
482 | } | |
483 | if(n==1) variable[1]=999; | |
484 | ||
3ed21f0a | 485 | TCanvas* cvpj=new TCanvas(Form("proj%d",variable[0]),Form("%s wrt %s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data())); |
6b64765b | 486 | |
3ed21f0a | 487 | 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())); |
488 | TLegend *leg=new TLegend(0.7,0.2,0.9,0.6,"Pt Bin"); | |
489 | leg->SetBorderSize(0); | |
490 | leg->SetFillStyle(0); | |
491 | ||
492 | for (Int_t i=0;i<nptbins;i++){ //loop on ptbins | |
493 | cout<<"\nPtBin = "<<i<<endl; | |
494 | ||
495 | //using AliSignificanceCalculator | |
496 | ||
497 | TString nameS,nameB,nameerrS,nameerrB; | |
498 | nameS.Form("SmultiDimVectorPtBin%d",i); | |
499 | nameerrS.Form("errSmultiDimVectorPtBin%d",i); | |
500 | nameB.Form("BmultiDimVectorPtBin%d",i); | |
501 | nameerrB.Form("errBmultiDimVectorPtBin%d",i); | |
502 | ||
503 | AliMultiDimVector* mdvS=(AliMultiDimVector*)fin->Get(nameS.Data()); | |
504 | AliMultiDimVector* mdvB=(AliMultiDimVector*)fin->Get(nameB.Data()); | |
505 | AliMultiDimVector* mdvBerr=(AliMultiDimVector*)fin->Get(nameerrS.Data()); | |
506 | AliMultiDimVector* mdvSerr=(AliMultiDimVector*)fin->Get(nameerrB.Data()); | |
507 | if(!(mdvS && mdvB && mdvSerr && mdvBerr)){ | |
508 | cout<<"one of the multidimvector is not present"<<endl; | |
509 | return; | |
510 | } | |
511 | ||
512 | AliSignificanceCalculator *cal=new AliSignificanceCalculator(mdvS,mdvB,mdvSerr,mdvBerr,1.,1.); | |
513 | ||
3ed21f0a | 514 | AliMultiDimVector* mvess=cal->GetSignificanceError(); |
3ed21f0a | 515 | AliMultiDimVector* mvpur=cal->CalculatePurity(); |
516 | AliMultiDimVector* mvepur=cal->CalculatePurityError(); | |
6b64765b | 517 | |
3ed21f0a | 518 | Int_t ncuts=mdvS->GetNVariables(); |
519 | Int_t *maxInd=new Int_t[ncuts]; | |
520 | Float_t *cutvalues=new Float_t[ncuts]; | |
6b64765b | 521 | //init |
3ed21f0a | 522 | for(Int_t ind=0;ind<ncuts;ind++)maxInd[ind]=0; |
6b64765b | 523 | |
3ed21f0a | 524 | Float_t sigMax0=cal->GetMaxSignificance(maxInd,0); |
525 | for(Int_t ic=0;ic<ncuts;ic++){ | |
526 | cutvalues[ic]=((AliMultiDimVector*)fin->Get(nameS.Data()))->GetCutValue(ic,maxInd[ic]); | |
6b64765b | 527 | |
528 | //setting step of fixed variables | |
529 | if(readfromfile){ //from file | |
530 | fixedvars[ic]=allfixedvars[i+ic]; | |
531 | } | |
532 | ||
533 | if(maximize) { //using the values which maximize the significance | |
534 | fixedvars[ic]=maxInd[ic]; | |
535 | //write to output fixedvars.dat | |
536 | writefixedvars<<fixedvars[ic]<<"\t"; | |
537 | } | |
3ed21f0a | 538 | } |
6b64765b | 539 | //output file: return after each pt bin |
540 | if(maximize) writefixedvars<<endl; | |
541 | ||
3ed21f0a | 542 | printf("Maximum of significance for Ptbin %d found in bin:\n",i); |
543 | for(Int_t ic=0;ic<ncuts;ic++){ | |
544 | printf(" %d\n",maxInd[ic]); | |
545 | printf("corresponding to cut:\n"); | |
546 | printf(" %f\n",cutvalues[ic]); | |
547 | } | |
3ed21f0a | 548 | |
6b64765b | 549 | printf("Significance = %f +- %f\n",sigMax0,mvess->GetElement(maxInd,0)); |
550 | printf("Purity = %f +- %f\n",mvpur->GetElement(maxInd,0),mvepur->GetElement(maxInd,i)); | |
3ed21f0a | 551 | |
c24d2d9f | 552 | if(which==3){ |
553 | //mdv=0x0; | |
554 | mdv=cal->CalculateSOverB(); | |
555 | if(!mdv)cout<<mdv->GetName()<<" null"<<endl; | |
556 | //mdverr=0x0; | |
557 | mdverr=cal->CalculateSOverBError(); | |
558 | if(!mdverr)cout<<mdverr->GetName()<<" null"<<endl; | |
559 | }else{ | |
560 | ||
561 | //multidimvector | |
562 | mdvname=Form("%s%d",name.Data(),i); | |
563 | mdv=(AliMultiDimVector*)fin->Get(mdvname); | |
564 | if(!mdv)cout<<mdvname.Data()<<" not found"<<endl; | |
565 | ||
566 | //multidimvector of errors | |
567 | mdverrname=Form("err%s%d",name.Data(),i); | |
568 | mdverr=(AliMultiDimVector*)fin->Get(mdverrname); | |
569 | if(!mdverr)cout<<mdverrname.Data()<<" not found"<<endl; | |
570 | } | |
3ed21f0a | 571 | TString ptbinrange=Form("%.0f < p_{t} < %.0f GeV/c",mdv->GetPtLimit(0),mdv->GetPtLimit(1)); |
572 | ||
573 | if(n==2) { | |
574 | gStyle->SetPalette(1); | |
575 | TH2F* hproj=mdv->Project(variable[0],variable[1],fixedvars,0); | |
576 | hproj->SetTitle(Form("%s wrt %s vs %s (Ptbin%d);%s;%s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data(),i,(mdv->GetAxisTitle(variable[0])).Data(),mdv->GetAxisTitle(variable[1]).Data())); | |
c24d2d9f | 577 | |
3ed21f0a | 578 | 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)); |
579 | cvpj->cd(); | |
6b64765b | 580 | hproj->DrawClone("COLZtext"); |
c24d2d9f | 581 | cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName())); |
3ed21f0a | 582 | delete hproj; |
583 | } | |
584 | ||
585 | if(n==1){ | |
6b64765b | 586 | |
3ed21f0a | 587 | Int_t nbins=mdv->GetNCutSteps(variable[0]); |
588 | ||
589 | Double_t *x=new Double_t[nbins]; | |
590 | Double_t *y=new Double_t[nbins]; | |
591 | Double_t *errx=new Double_t[nbins]; | |
592 | Double_t *erry=new Double_t[nbins]; | |
593 | ||
594 | for(Int_t k=0;k<nbins;k++){ //loop on the steps (that is the bins of the graph) | |
6b64765b | 595 | //init |
3ed21f0a | 596 | x[k]=0;y[k]=0; |
597 | errx[k]=0;erry[k]=0; | |
598 | ||
6b64765b | 599 | fixedvars[variable[0]]=k; //variable[0] is the index of the variable on which we project. This variable must increase, the others have been set before |
3ed21f0a | 600 | |
601 | x[k]=mdv->GetCutValue(variable[0],k); | |
602 | errx[k]=mdv->GetCutStep(variable[0])/2.; | |
603 | ||
604 | y[k]=mdv->GetElement(fixedvars,0); | |
605 | erry[k]=mdverr->GetElement(fixedvars,0); | |
c24d2d9f | 606 | if(which==3){ |
3ed21f0a | 607 | |
c24d2d9f | 608 | } |
6b64765b | 609 | cout<<mdv->GetAxisTitle(variable[0])<<" step "<<k<<" = "<<x[k]<<":"<<" y = "<<y[k]<<endl; |
3ed21f0a | 610 | } |
6b64765b | 611 | |
612 | cout<<"----------------------------------------------------------"<<endl; | |
3ed21f0a | 613 | TGraphErrors* gr=new TGraphErrors(nbins,x,y,errx,erry); |
614 | gr->SetMarkerStyle(20+i); | |
615 | gr->SetMarkerColor(i+1); | |
616 | gr->SetMinimum(0); | |
617 | ||
618 | gr->SetName(Form("g1%d",i)); | |
619 | mg->Add(gr,"P"); | |
620 | leg->AddEntry(gr,ptbinrange.Data(),"p"); | |
621 | } | |
622 | } | |
623 | ||
624 | if(n==1){ | |
625 | cvpj->cd(); | |
626 | mg->Draw("A"); | |
627 | leg->Draw(); | |
c24d2d9f | 628 | cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName())); |
6b64765b | 629 | } else delete cvpj; |
4bcb0ebb | 630 | |
3ed21f0a | 631 | } |
632 |