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