Added D* (Raoul)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / charmCutsOptimization.C
1 #include <fstream>
2 #include <Riostream.h>
3 #include <TSystem.h>
4 #include <TMath.h>
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>
15 #include <TKey.h>
16 #include <TObjectTable.h>
17 #include <TDatabasePDG.h>
18 #include <TMath.h>
19 #include <TPaveText.h>
20 #include <TText.h>
21
22 #include <AliMultiDimVector.h>
23 #include "AliHFMassFitter.h"
24 #include <AliSignificanceCalculator.h>
25
26 #include <fstream>
27
28 //global variables
29 Double_t nsigma=3;
30 Int_t decCh=2;
31 Int_t fitbtype=5;
32 Int_t rebin=2;
33 Double_t sigma=0.0005;
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
41
42 ofstream outcheck;
43 ofstream outdetail;
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);
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);
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
49 //decCh:
50 //- 0 = kDplustoKpipi
51 //- 1 = kD0toKpi
52 //- 2 = kDStartoKpipi
53 //- 3 = kDstoKKpi
54 //- 4 = kD0toKpipipi
55 //- 5 = kLambdactopKpi
56
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
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
68   outcheck.open("output.dat",ios::out);
69   outdetail.open("discarddetails.dat",ios::out);
70
71   gStyle->SetFrameBorderMode(0);
72   gStyle->SetCanvasColor(0);
73   gStyle->SetFrameFillColor(0);
74   gStyle->SetTitleFillColor(0);
75   gStyle->SetStatColor(0);
76
77   //~/Lavoro/PbPb/tagli/SIGAOD33/mar02/cent3070/
78   TString filename="AnalysisResults.root",dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv";
79
80   TString hnamemass="hMass_",hnamesgn="hSig_",hnamebkg="hBkg_";
81
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:
95     listname+="Dstar0100";
96     mdvlistname+="Dstar0100";
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();
119   if (decCh==2) mass=(TDatabasePDG::Instance()->GetParticle(pdg)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); 
120
121   if(part!="both"){
122     listname.Append(part);
123     mdvlistname.Append(part);
124   }
125   if(centr!="no"){
126     listname.Append(centr);
127     mdvlistname.Append(centr);
128   }
129   cout<<"Mass = "<<mass<<endl;
130   
131   Int_t countFitFail=0,countSgnfFail=0,countNoHist=0,countBkgOnly=0;
132
133   outcheck<<"ptbin\tmdvGlobAddr\thistIndex\tSignif\tS\tB"<<endl;
134   outdetail<<"ptbin\tmdvGlobAddr\thistIndex\trelErrS\t\tmean_F-mass (mass = "<<mass<<")"<<endl;
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
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");
165     cst->SaveAs("hstat.png");
166   }else{
167     cout<<"Warning! fHistNEvents not found in "<<listname.Data()<<endl;
168   }
169
170   Bool_t isMC=kFALSE;
171   TH1F* htestIsMC=(TH1F*)histlist->FindObject("hSig_0");
172   if(htestIsMC) isMC=kTRUE;
173
174   Int_t nptbins=listamdv->GetEntries();
175   Int_t nhist=(histlist->GetEntries()-1);//-1 because of fHistNevents
176   if(isMC) nhist/=4; ///4 because hMass_, hSgn_,hBkg_,hRfl_
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
186   TH1F** hSigma=new TH1F*[nptbins];
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));
194
195   //Check wheter histograms are filled
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());
200
201       if(!h){
202         cout<<name<<" not found"<<endl;
203         continue;
204       }
205
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         }
213       }
214     }
215   
216
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     }
222   }
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());
234     outcheck<<"\n"<<mdvS->GetPtLimit(0)<<" < Pt <"<<mdvS->GetPtLimit(1)<<endl;
235
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
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
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]);
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);
275
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;
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);
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;
299         }
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;
305         }
306
307         MC(h,g,signal,errSignal,background,errBackground,signif,errSignif,sigmafit,status);
308
309       }
310
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);
321       }else{
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);
328         if(status==3){
329           mdvB->SetElement(ih,background);
330           mdvBerr->SetElement(ih,errBackground);
331         }
332
333       }
334
335     }
336   
337     fout->cd();
338     mdvS->Write();
339     mdvB->Write();
340     mdvSgnf->Write();
341     mdvSerr->Write();
342     mdvBerr->Write();
343     mdvSgnferr->Write();
344     hSigma[i]->Write();
345
346   }
347   
348
349   TCanvas *cinfo=new TCanvas("cinfo","Status");
350   cinfo->cd();
351   cinfo->SetGrid();
352   hstatus->Draw("htext0");
353   fout->cd();
354   hstatus->Write();
355   fout->Close();
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();
358
359 cout << "test1\n";
360   return kTRUE;
361 }
362
363
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
376   if(decCh != 2) min = h->GetBinLowEdge(1);
377   else min = TDatabasePDG::Instance()->GetParticle(211)->Mass();
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 }
488
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
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
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
537
538
539   Int_t nbin=h->GetNbinsX();
540   if(decCh != 2) Double_t min = h->GetBinLowEdge(1);
541   else Double_t min = TDatabasePDG::Instance()->GetParticle(211)->Mass();
542   Double_t max=h->GetBinLowEdge(nbin+1);
543
544   if(rangefit) {
545     min=rangefit[0];
546     max=rangefit[1];
547   }
548
549   reject = true;
550 //Bkg fit : exponential = A*exp(B*x) 
551   if(decCh != 2) TF1 *fBkgFit = new TF1("fBkgFit",ExpoBkgWoPeak,min,max,4);          
552   else { 
553    if (fitbtype == 5) TF1 *fBkgFit = new TF1("fBkgFit",PowerExpoBkgWoPeak,min,max,4); //Bkg fit : PowerExpo = A*(x-mpi)^(1/2)*exp(-B*(x-mpi))
554    else TF1 *fBkgFit = new TF1("fBkgFit",PowerBkgWoPeak,min,max,4); //Bkg fit : Power = A*(x-mpi)^B
555   }
556
557   fBkgFit->FixParameter(2,mass-nsigma*sigma);
558   fBkgFit->FixParameter(3,mass+nsigma*sigma);
559   TFitResultPtr r = h->Fit(fBkgFit,"RS+");
560   Int_t ok = r;
561
562   if(ok==-1){
563     cout<<"FIT NOT OK!"<<endl;
564     cout<<"\t 0\t xxx"<<"\t failed"<<endl;
565     status=0;
566     return kFALSE;
567   } 
568
569   reject = false;
570   if(decCh !=2) TF1 *fBkgFct = new TF1("fBkgFct",ExpoBkgWoPeak,min,max,4);  
571   else {
572    if (fitbtype == 5) TF1 *fBkgFct = new TF1("fBkgFct",PowerExpoBkgWoPeak,min,max,4);  
573    else TF1 *fBkgFct = new TF1("fBkgFct",PowerBkgWoPeak,min,max,4);
574   }
575   fBkgFct->SetLineStyle(2);
576   for(Int_t i=0; i<4; i++) fBkgFct->SetParameter(i,fBkgFit->GetParameter(i));
577   h->GetListOfFunctions()->Add(fBkgFct);
578   TH1F * hBkgFct = (TH1F*)fBkgFct->GetHistogram();
579
580
581   status = 1;    
582   Double_t binStartCount = h->FindBin(mass-nsigma*sigma);
583   Double_t binEndCount = h->FindBin(mass+nsigma*sigma);
584   Double_t counts=0., bkgcounts=0., errcounts=0., errbkgcounts=0.;
585   for (Int_t ibin = binStartCount; ibin<=binEndCount; ibin++) {
586     counts += h->GetBinContent( ibin );
587     errcounts += counts ;
588     Double_t center =  h->GetBinCenter(ibin);
589     bkgcounts += hBkgFct->GetBinContent( hBkgFct->FindBin(center) ); 
590     
591     errbkgcounts += bkgcounts ;
592   }
593   
594   bkg = bkgcounts;
595   errbkg = TMath::Sqrt( errbkgcounts );
596   sgn = counts - bkg ;
597   if(sgn<0) status = 2; // significance < 0
598   errsgn = TMath::Sqrt( counts + errbkg*errbkg );
599   sgnf = sgn / TMath::Sqrt( sgn + bkg );
600   errsgnf = TMath::Sqrt( sgnf*sgnf/(sgn+bkg)/(sgn+bkg)*(1/4.*errsgn*errsgn+errbkg*errbkg) + sgnf*sgnf/sgn/sgn*errsgn*errsgn );
601       cout << " Signal "<<sgn<<" +- "<<errsgn<<", bkg "<<bkg<<" +- "<<errbkg<<", significance "<<sgnf<<" +- "<<errsgnf<<endl;
602   
603   if(writefit) {
604
605     TString filename = Form("%sMassFit.root",h->GetName());
606
607     TFile* outputcv = new TFile(filename.Data(),"recreate");      
608
609     TCanvas* c = new TCanvas();
610
611     c->SetName(Form("%s",h->GetName()));
612
613     h->Draw();
614
615     TPaveText *pavetext=new TPaveText(0.4,0.7,0.85,0.9,"NDC");     
616     pavetext->SetBorderSize(0);
617     pavetext->SetFillStyle(0);
618     pavetext->AddText(Form("Signal = %4.2e #pm %4.2e",sgn,errsgn));
619     pavetext->AddText(Form("Bkg = %4.2e #pm %4.2e",bkg,errbkg));
620     pavetext->AddText(Form("Signif = %3.2f #pm %3.2f",sgnf,errsgnf));
621     c->cd();
622
623     pavetext->DrawClone();
624
625     outputcv->cd();
626
627     c->Write();
628
629
630     outputcv->Close();
631
632     delete outputcv;
633
634     delete c;
635   }
636
637
638   delete fBkgFit;
639
640   delete fBkgFct;
641
642
643   return kTRUE; 
644 }
645
646 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
647
648 // which=0 plot significance
649 //      =1 plot signal
650 //      =2 plot background
651 //      =3 plot S/B
652 // 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
653 // 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)
654
655
656 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){
657
658   gStyle->SetCanvasColor(0);
659   gStyle->SetFrameFillColor(0);
660   gStyle->SetTitleFillColor(0);
661   gStyle->SetOptStat(0);
662   //gStyle->SetOptTitle(0);
663   gStyle->SetFrameBorderMode(0);
664
665   TFile* fin=new TFile("outputSignifMaxim.root");
666   if(!fin->IsOpen()){
667     cout<<"outputSignifMaxim.root not found"<<endl;
668     return;
669   }
670
671   if(n>2){
672     cout<<"Error! cannot show "<<n+1<<" dimentions"<<endl;
673     return;
674   }
675
676   TString name,title,namebis,shorttitle;
677   switch (which){
678   case 0:
679     name="SgfmultiDimVectorPtBin";
680     title="Significance";
681     shorttitle="Sgnf";
682     break;
683   case 1:
684     name="SmultiDimVectorPtBin";
685     title="Signal";
686     shorttitle="S";
687     break;
688   case 2:
689     name="BmultiDimVectorPtBin";
690     title="Background";
691     shorttitle="B";
692     break;
693   case 3:
694     name="SmultiDimVectorPtBin";
695     namebis="BmultiDimVectorPtBin";
696     title="Signal over Background ";
697     shorttitle="SoB";
698     break;
699   // case 4:
700   //   name="errBmultiDimVectorPtBin";
701   //   title="Background (error)";
702   //   break;
703   }
704  
705   if(plotErrors && which!=3 && n==2){
706     name.Prepend("err");
707     title.Append(" Error") ;
708     shorttitle.Append("Err");
709   }
710
711   Int_t nptbins=50;
712   
713   for(Int_t ip=0;ip<=nptbins;ip++){
714     TString mdvname=Form("%s%d",name.Data(),ip);
715     AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname);
716     if(!mdv){
717       nptbins=ip;
718       cout<<"Number of pt bins "<<ip<<endl;
719       break;
720     }
721   }
722
723   cout<<"Projecting "<<title.Data()<<" with respect to the maximization variable(s) [chose]"<<endl;
724  
725   Int_t variable[2]; //no more than 2D
726   TString mdvname=Form("%s0",name.Data()), mdverrname="";//, mdvnamebis="", mdverrnamebis="";
727   AliMultiDimVector* mdv=(AliMultiDimVector*)fin->Get(mdvname);
728   AliMultiDimVector* mdverr=0x0;
729   if(!mdv){
730     cout<<mdvname.Data()<<" not found"<<endl;
731     return;
732   }
733   
734   Int_t nvarsopt=mdv->GetNVariables();
735   //Int_t nfixed=nvarsopt-n;
736  
737   Int_t fixedvars[nvarsopt];
738   Int_t allfixedvars[nvarsopt*nptbins];
739  
740
741
742   fstream writefixedvars;
743   if(readfromfile) {
744     //open file in read mode
745     writefixedvars.open("fixedvars.dat",ios::in);
746     Int_t longi=0;
747     while(writefixedvars){
748       writefixedvars>>allfixedvars[longi];
749       longi++;
750     }
751   }
752   else {
753     //open file in write mode
754     writefixedvars.open("fixedvars.dat",ios::out);
755   }
756
757   Bool_t freevars[nvarsopt];
758
759   //ask variables for projection
760   for(Int_t k=0;k<nvarsopt;k++){
761     cout<<k<<" "<<mdv->GetAxisTitle(k)<<endl;
762     freevars[k]=kTRUE;
763   }
764   cout<<"Choose "<<n<<" variable(s)"<<endl;
765   for(Int_t j=0;j<n;j++){
766     cout<<"var"<<j<<": ";
767     cin>>variable[j];
768     freevars[variable[j]]=kFALSE;
769   }
770   if(n==1) variable[1]=999;
771
772   TCanvas* cvpj=new TCanvas(Form("proj%d",variable[0]),Form("%s wrt %s",title.Data(),(mdv->GetAxisTitle(variable[0])).Data()));
773
774   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()));
775   TLegend *leg=new TLegend(0.7,0.2,0.9,0.6,"Pt Bin");
776   leg->SetBorderSize(0);
777   leg->SetFillStyle(0);
778
779   //set scale
780   fstream writerange;
781   Float_t axisrange[2*nptbins];
782   if(fixedrange) {
783     //open file in read mode
784     writerange.open("rangehistos.dat",ios::in);
785     Int_t longi=0;
786     while(writerange){
787       writerange>>axisrange[longi];
788       longi++;
789     }
790   }
791   else {
792     //open file in write mode
793     writerange.open("rangehistos.dat",ios::out);
794   }
795
796   for (Int_t i=0;i<nptbins;i++){   //loop on ptbins
797     cout<<"\nPtBin = "<<i<<endl;
798
799     //using AliSignificanceCalculator
800
801     TString nameS,nameB,nameerrS,nameerrB;
802     nameS.Form("SmultiDimVectorPtBin%d",i);
803     nameerrS.Form("errSmultiDimVectorPtBin%d",i);
804     nameB.Form("BmultiDimVectorPtBin%d",i);
805     nameerrB.Form("errBmultiDimVectorPtBin%d",i);
806  
807     AliMultiDimVector* mdvS=(AliMultiDimVector*)fin->Get(nameS.Data());
808     AliMultiDimVector* mdvB=(AliMultiDimVector*)fin->Get(nameB.Data());
809     AliMultiDimVector* mdvBerr=(AliMultiDimVector*)fin->Get(nameerrS.Data());
810     AliMultiDimVector* mdvSerr=(AliMultiDimVector*)fin->Get(nameerrB.Data());
811     if(!(mdvS && mdvB && mdvSerr && mdvBerr)){
812       cout<<"one of the multidimvector is not present"<<endl;
813       return;
814     }
815
816     AliSignificanceCalculator *cal=new AliSignificanceCalculator(mdvS,mdvB,mdvSerr,mdvBerr,1.,1.);
817
818     AliMultiDimVector* mvess=cal->GetSignificanceError();
819     AliMultiDimVector* mvpur=cal->CalculatePurity();
820     AliMultiDimVector* mvepur=cal->CalculatePurityError();
821
822     Int_t ncuts=mdvS->GetNVariables();
823     Int_t *maxInd=new Int_t[ncuts];
824     Float_t *cutvalues=new Float_t[ncuts];
825     //init
826     // for(Int_t ind=0;ind<ncuts;ind++)maxInd[ind]=0;
827
828     Float_t sigMax0=cal->GetMaxSignificance(maxInd,0); //look better into this!!
829
830     for(Int_t ic=0;ic<ncuts;ic++){
831       cutvalues[ic]=((AliMultiDimVector*)fin->Get(nameS.Data()))->GetCutValue(ic,maxInd[ic]);
832
833       //setting step of fixed variables
834       if(readfromfile){ //from file
835         fixedvars[ic]=allfixedvars[i+ic];
836       }
837
838       if(!readfromfile) { //using the values which maximize the significance
839         fixedvars[ic]=maxInd[ic];
840         //write to output fixedvars.dat
841         writefixedvars<<fixedvars[ic]<<"\t";
842       }
843     }
844     //output file: return after each pt bin
845     if(!readfromfile) writefixedvars<<endl;
846
847     printf("Maximum of significance for Ptbin %d found in bin:\n",i);
848     for(Int_t ic=0;ic<ncuts;ic++)printf(" %d  ",maxInd[ic]);
849     printf("\ncorresponding to cut:\n");
850     for(Int_t ic=0;ic<ncuts;ic++)printf(" %f  ",cutvalues[ic]);
851     
852     printf("\nSignificance = %f +- %f\n",sigMax0,mvess->GetElement(maxInd,0));
853     printf("Purity       = %f +- %f\n",mvpur->GetElement(maxInd,0),mvepur->GetElement(maxInd,i));
854     
855     if(which==3){
856       //mdv=0x0;
857       mdv=cal->CalculateSOverB();
858       if(!mdv)cout<<mdv->GetName()<<" null"<<endl;
859       //mdverr=0x0;
860       mdverr=cal->CalculateSOverBError();
861       if(!mdverr)cout<<mdverr->GetName()<<" null"<<endl;
862     }else{
863
864       //multidimvector
865       mdvname=Form("%s%d",name.Data(),i);   
866       mdv=(AliMultiDimVector*)fin->Get(mdvname);
867       if(!mdv)cout<<mdvname.Data()<<" not found"<<endl;
868       
869       //multidimvector of errors
870       mdverrname=Form("err%s%d",name.Data(),i);   
871       mdverr=(AliMultiDimVector*)fin->Get(mdverrname);
872       if(!mdverr)cout<<mdverrname.Data()<<" not found"<<endl;
873     }
874     printf("Global Address %d (%d)\n",(Int_t)mdv->GetGlobalAddressFromIndices(maxInd,0),(Int_t)mdv->GetNTotCells()*i+(Int_t)mdv->GetGlobalAddressFromIndices(maxInd,0));
875     TString ptbinrange=Form("%.1f < p_{t} < %.1f GeV/c",mdv->GetPtLimit(0),mdv->GetPtLimit(1));
876
877     Float_t maxval=0;
878
879     if(n==2) {
880       gStyle->SetPalette(1);
881       Int_t steps[2];
882       Int_t nstep[2]={mdv->GetNCutSteps(variable[0]),mdv->GetNCutSteps(variable[1])}; 
883   
884       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]));
885       if(fixedplane){
886         hproj=mdv->Project(variable[0],variable[1],fixedvars,0);
887         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()));
888       }else{
889         for(Int_t ist1=0;ist1<nstep[0];ist1++){
890           steps[0]=ist1;
891           Int_t fillbin1=ist1+1;
892           if(mdv->GetCutValue(variable[0],0)>mdv->GetCutValue(variable[0],mdv->GetNCutSteps(variable[0])-1))fillbin1=nstep[0]-ist1;
893           for(Int_t ist2=0;ist2<nstep[1];ist2++){
894             steps[1]=ist2;
895             Int_t fillbin2=ist2+1;
896             if(mdv->GetCutValue(variable[1],0)>mdv->GetCutValue(variable[1],mdv->GetNCutSteps(variable[1])-1))fillbin2=nstep[1]-ist2;
897             Int_t* varmaxim=mdv->FindLocalMaximum(maxval,variable,steps,n,0);
898             hproj->SetBinContent(fillbin1,fillbin2,maxval);
899             delete varmaxim;
900           }
901         }
902       }
903       if(fixedrange) {
904         hproj->SetMinimum(axisrange[2*i]);
905         hproj->SetMaximum(axisrange[2*i+1]);
906       } else{
907         writerange<<hproj->GetMinimum()<<"\t"<<hproj->GetMinimum()<<endl;
908       }
909       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));
910       cvpj->cd();
911       hproj->DrawClone("COLZtext");
912       cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName()));
913       delete hproj;
914     }
915
916     if(n==1){
917
918       Int_t nbins=mdv->GetNCutSteps(variable[0]);
919  
920       Double_t *x=new Double_t[nbins];
921       Double_t *y=new Double_t[nbins];
922       Double_t *errx=new Double_t[nbins];
923       Double_t *erry=new Double_t[nbins];
924
925       for(Int_t k=0;k<nbins;k++){ //loop on the steps (that is the bins of the graph)
926         //init
927         x[k]=0;y[k]=0;
928         errx[k]=0;erry[k]=0;
929         
930         x[k]=mdv->GetCutValue(variable[0],k);
931         errx[k]=mdv->GetCutStep(variable[0])/2.;
932         Int_t onevariable[1]={variable[0]};
933         Int_t onestep[1]={k};
934         ULong64_t gladd;
935
936         Float_t maxval;
937         Int_t* varmaxim=mdv->FindLocalMaximum(maxval,onevariable,onestep,n,0);
938         y[k]=maxval;
939         gladd=mdv->GetGlobalAddressFromIndices(varmaxim,0);
940         cout<<gladd<<endl;
941         delete varmaxim;
942       
943         erry[k]=mdverr->GetElement(gladd);
944         
945         cout<<mdv->GetAxisTitle(variable[0])<<" step "<<k<<" = "<<x[k]<<":"<<" y = "<<y[k]<<endl;
946       }
947       
948       cout<<"----------------------------------------------------------"<<endl;
949       TGraphErrors* gr=new TGraphErrors(nbins,x,y,errx,erry);
950       gr->SetMarkerStyle(20+i);
951       gr->SetMarkerColor(i+1);
952       gr->SetLineColor(i+1);
953       if(i>=9){
954         gr->SetMarkerColor(i+2);
955         gr->SetLineColor(i+2);
956       }
957       gr->SetMinimum(0);
958
959       gr->SetName(Form("g1%d",i));
960       mg->Add(gr,"P");
961       leg->AddEntry(gr,ptbinrange.Data(),"p");
962     }
963   }
964    
965   if(n==1){
966     cvpj->cd();
967     mg->Draw("A");
968     leg->Draw();
969     cvpj->SaveAs(Form("%s%s.png",shorttitle.Data(),cvpj->GetName()));
970   } else delete cvpj;
971 }
972
973 //draw sigma as a function of cuts
974
975 void DrawSigmas(TH2F* h2cuts){
976   TFile *fin=0x0;
977   TString fittype="ExpFit";
978   Int_t ntot=5;
979   if(fittype=="Pol2Fit") ntot=6;
980   Int_t ihfirst=0,ihlast=1; //change this (must think on it and remember what I wanted to do!)
981   for(Int_t ih=ihfirst;ih<ihlast;ih++){
982     fin=new TFile(Form("h%d%s.root",ih,fittype.Data()));
983     if(!fin) continue;
984     TCanvas *cv=(TCanvas*)fin->Get(Form("cv1%s%d",fittype.Data(),ih));
985     TF1* func=(TF1*)cv->FindObject("funcmass");
986     Int_t sigma=func->GetParameter(ntot-1);
987     //h2cuts->SetBinContent();
988     //to be finished
989   }
990 }
991
992 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
993 // Some methods to get the index of the histogram corresponding to a set of cuts
994
995 // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root
996 // ptbin = pt bin
997 // indices = array of the index of the cut variables (the dimension of the array must be equal to the number of variables maximized)
998
999 Int_t GetNHistFromIndices(AliMultiDimVector* vct,Int_t ptbin,Int_t* indices){
1000   cout<<"Calculating the index of the histogram corresponding to the following cut steps:"<<endl;
1001   for(Int_t i=0;i<vct->GetNVariables();i++){
1002     cout<<vct->GetAxisTitle(i)<<" --> "<<indices[i]<<endl;
1003   }
1004   cout<<"Pt bin "<<ptbin<<" from "<<vct->GetPtLimit(0)<<" to "<<vct->GetPtLimit(1)<<endl;
1005   cout<<"Info: total number of cells per multidim: "<<vct->GetNTotCells()<<endl;
1006   ULong64_t glindex=vct->GetGlobalAddressFromIndices(indices,0);
1007   cout<<"The histogram you want is:\t";
1008   return glindex+vct->GetNTotCells()*ptbin;
1009 }
1010
1011 // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root
1012 // ptbin = pt bin
1013 // values = array of the cut values (the dimension of the array must be equal to the number of variables maximized)
1014
1015 Int_t GetNHistFromValues(AliMultiDimVector* vct,Int_t ptbin,Float_t* values){
1016   cout<<"Calculating the index of the histogram corresponding to the following cut values:"<<endl;
1017   for(Int_t i=0;i<vct->GetNVariables();i++){
1018     cout<<vct->GetAxisTitle(i)<<" --> "<<values[i]<<endl;
1019   }
1020   cout<<"Pt bin "<<ptbin<<" from "<<vct->GetPtLimit(0)<<" to "<<vct->GetPtLimit(1)<<endl;
1021   cout<<"Info: total number of cells per multidim: "<<vct->GetNTotCells()<<endl;
1022   ULong64_t glindex=vct->GetGlobalAddressFromValues(values,0);
1023  
1024   cout<<"The histogram you want is:\t"<<glindex+vct->GetNTotCells()*ptbin<<endl;
1025   return glindex+vct->GetNTotCells()*ptbin;
1026 }
1027
1028 // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root
1029 // ptbin = pt bin
1030 // values = array of the cut values: the dimention can be <= number of variables maximized
1031 // valsgiven = array of dimention = to the number of variables optimized. For each variable put kTRUE if the value is given (in values), kFALSE otherwise
1032 // nhistinrange =  pass an integer which will contains the number of histogram returned (that is the dimention of the Int_t* returned)
1033
1034 //NB: Remember that the cut applied is the lower edge of the step where lower=looser
1035
1036 Int_t* GetRangeHistFromValues(AliMultiDimVector* vct,Int_t ptbin,Bool_t* valsgiven,Float_t* values,Int_t& nhistinrange){
1037   Int_t nvargiven=0;
1038   nhistinrange=1;
1039
1040   Int_t nvar4opt=vct->GetNVariables();
1041   Float_t allvals[nvar4opt];
1042  
1043   for (Int_t i=0;i<nvar4opt;i++) {
1044     if(valsgiven[i]==kTRUE) {
1045       allvals[i]=values[nvargiven];
1046       nvargiven++;
1047     }
1048     else {
1049       nhistinrange+=vct->GetNCutSteps(i);
1050       allvals[i]=vct->GetCutValue(i,0);
1051       //allvals[i]=vct->GetCutValue(0,i); //ivar,icell
1052     }
1053   }
1054   cout<<nhistinrange<<" index will be returned"<<endl;
1055   Int_t *rangeofhistos=new Int_t[nhistinrange];
1056
1057   if(nhistinrange==1){
1058     rangeofhistos[0]=GetNHistFromValues(vct,ptbin,allvals);
1059     cout<<"output"<<rangeofhistos[0]<<endl;
1060   }else{
1061     Int_t index[nvar4opt-nvargiven];
1062     Int_t k=0;
1063     for (Int_t i=0;i<nvar4opt;i++){
1064       if(valsgiven[i]==kFALSE) {
1065         //cout<<"kTRUE==>"<<i<<endl;
1066         index[k]=i;
1067         k++;
1068       }
1069     }
1070
1071     for(Int_t i=0;i<nvar4opt-nvargiven;i++){ //loop on number of free variables
1072       cout<<"Info: incrementing "<<vct->GetAxisTitle(index[i])<<endl;
1073       for(Int_t j=0;j<vct->GetNCutSteps(i);j++){ //loop on steps of each free variable
1074         allvals[index[i]]=vct->GetCutValue(index[i],j);
1075         rangeofhistos[i*vct->GetNCutSteps(i)+j]=GetNHistFromValues(vct,ptbin,allvals);
1076       }
1077     }
1078   }
1079   return rangeofhistos;
1080 }
1081
1082 // vct = the AliMultiDimVector correponding to ptbin: it can be from AnalysisResults.root or outputSignifMaxim.root
1083 // ptbin = pt bin
1084 // nhist = number of the histogram from which you want to have the cut values (returned)
1085
1086 Float_t* GetCutValuesFromNHist(AliMultiDimVector* vct,Int_t ptbin,Int_t nhist){
1087   ULong64_t totCells=vct->GetNTotCells();
1088   ULong64_t globadd=nhist-ptbin*totCells;
1089   const Int_t nvars=vct->GetNVariables();
1090   Float_t* cuts=new Float_t[nvars];
1091   Int_t ptinside;
1092   vct->GetCutValuesFromGlobalAddress(globadd,cuts,ptinside);
1093   return cuts;
1094 }
1095
1096 // ptbin = pt bin
1097 // values = array of the cut values: the dimention can be <= number of variables maximized
1098 // valsgiven = array of dimention = to the number of variables optimized. For each variable put kTRUE if the value is given (in values), kFALSE otherwise
1099 // 
1100
1101 void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString path="./",Int_t decCh=2){
1102   gStyle->SetFrameBorderMode(0);
1103   gStyle->SetCanvasColor(0);
1104   gStyle->SetFrameFillColor(0);
1105   gStyle->SetOptStat(0);
1106
1107   Int_t nhists;
1108   TString filename="AnalysisResults.root";
1109   TString dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv";
1110   TString centr="020";
1111
1112   TFile *fin=new TFile(Form("%s%s",path.Data(),filename.Data()));
1113   if(!fin){
1114     cout<<path.Data()<<filename.Data()<<" not found"<<endl;
1115     return;
1116   }
1117   TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname);
1118   if(!dir){
1119     cout<<"Directory "<<dirname<<" not found"<<endl;
1120     return;
1121   }
1122   switch (decCh) {
1123   case 0:
1124     listname+="Dplus";
1125     mdvlistname+="Dplus";
1126     
1127     break;
1128   case 1:
1129     listname+="D0";
1130     mdvlistname+="D0";
1131     
1132     break;
1133   case 2:
1134     listname+="Dstar0100";
1135     mdvlistname+="Dstar0100";
1136     
1137     break;
1138   case 3:
1139     listname+="Ds";
1140     mdvlistname+="Ds";
1141     
1142     break;
1143   case 4:
1144     listname+="D04";
1145     mdvlistname+="D04";
1146     
1147     break;
1148   case 5:
1149     listname+="Lc";
1150     mdvlistname+="Lc";
1151     
1152     break;
1153   default:
1154     cout<<decCh<<" is not allowed as decay channel "<<endl;
1155     return;
1156   }
1157   mdvlistname+=centr;
1158   listname+=centr;
1159
1160   TList* listamdv= (TList*)dir->Get(mdvlistname);
1161   if(!listamdv) {
1162     cout<<mdvlistname<<" doesn't exist"<<endl;
1163     return;
1164   }
1165
1166   AliMultiDimVector* vct=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",ptbin));
1167
1168   TFile* fin2;
1169   TString filehistname="";
1170   Int_t* indexes=GetRangeHistFromValues(vct,ptbin,valsgiven,values,nhists);
1171   for(Int_t i=0;i<nhists;i++){
1172
1173     fin2=new TFile(Form("%sh%dExpMassFit.root",path.Data(),indexes[i]));
1174     if(!fin2){
1175       cout<<"File "<<indexes[i]<<" not found!"<<endl;
1176       continue;
1177     }else{
1178       TCanvas *cv=(TCanvas*)fin2->Get(Form("c1h%dExp",indexes[i]));
1179       if(!cv){
1180         cout<<"Canvas c1h"<<indexes[i]<<"Exp not found among";
1181         fin2->ls();
1182         continue;
1183       }else{
1184         cv->Draw();
1185         cv->SaveAs(Form("h%dExpMassFit.png",indexes[i]));
1186         fin2=0x0;
1187       }
1188     }
1189   }
1190 }
1191
1192 void Merge2Bins(Int_t b1, Int_t b2,TString pathin="./",Int_t decCh=2,TString part="both"/*"A" anti-particle, "P" particle*/){
1193
1194   if(b2!=b1+1){
1195     printf("The bins to be merget must be consecutive. Check! [b1 = %d, b2= %d]\n",b1,b2);
1196     return;
1197   }
1198
1199   TFile *fin=new TFile(Form("%sAnalysisResults.root",pathin.Data()));
1200   if (!fin){
1201     cout<<"Input file not found"<<endl;
1202     return;
1203   }
1204
1205   TString dirname="PWG3_D2H_Significance",listname="coutputSig",mdvlistname="coutputmv";
1206
1207   switch (decCh) {
1208   case 0:
1209     listname+="Dplus";
1210     mdvlistname+="Dplus";
1211     break;
1212   case 1:
1213     listname+="D0";
1214     mdvlistname+="D0";
1215     break;
1216   case 2:
1217     listname+="Dstar0100";
1218     mdvlistname+="Dstar0100";
1219     break;
1220   case 3:
1221     listname+="Ds";
1222     mdvlistname+="Ds";
1223     break;
1224   case 4:
1225     listname+="D04";
1226     mdvlistname+="D04";
1227     break;
1228   case 5:
1229     listname+="Lc";
1230     mdvlistname+="Lc";
1231     break;
1232   default:
1233     cout<<decCh<<" is not allowed as decay channel "<<endl;
1234     return;
1235   }
1236
1237   if(part!="both"){
1238     listname.Append(part);
1239     mdvlistname.Append(part);
1240   }
1241
1242   TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname);
1243   if(!dir){
1244     cout<<"Directory "<<dirname<<" not found"<<endl;
1245     return;
1246   }
1247
1248   TList* histlist= (TList*)dir->Get(listname);
1249   if(!histlist) {
1250     cout<<listname<<" doesn't exist"<<endl;
1251     return;
1252   }
1253
1254   TList* listamdv= (TList*)dir->Get(mdvlistname);
1255   if(!listamdv) {
1256     cout<<mdvlistname<<" doesn't exist"<<endl;
1257     return;
1258   }
1259   if (!gSystem->AccessPathName(Form("merged%d%d",b1,b2))) gSystem->Exec(Form("mkdir merged%d%d",b1,b2));
1260   gSystem->Exec(Form("cd merged%d%d",b1,b2));
1261
1262   TFile* fout=new TFile("mergeAnalysisResults.root","recreate");
1263
1264   fout->mkdir(dirname);
1265   TList* listmdvout=new TList();
1266   listmdvout->SetName(listamdv->GetName());
1267   listmdvout->SetOwner();
1268   //listmdvout->SetTitle(listamdv->GetTitle());
1269   TList* histlistout=new TList();
1270   histlistout->SetName(histlist->GetName());
1271   histlistout->SetOwner();
1272   //histlistout->SetTitle(histlist->GetTitle());
1273
1274   AliMultiDimVector* mdvin1=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",b1));
1275   AliMultiDimVector* mdvin2=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",b2));
1276
1277   Int_t ntotHperbin = mdvin1->GetNTotCells();
1278   if(mdvin2->GetNTotCells() != ntotHperbin) {
1279     cout<<"Error! Number of histos in pt bin "<<b1<<" = "<<ntotHperbin<<" != Number of histos in pt bin "<<b2<<" = "<<mdvin2->GetNTotCells()<<endl;
1280     return;
1281   }
1282   Int_t nvar1=mdvin1->GetNVariables();
1283   if(nvar1 != mdvin2->GetNVariables()){
1284     cout<<"Error! Mismatch in number of variables"<<endl;
1285     return;
1286   }
1287   
1288   Float_t newptbins[2]={mdvin1->GetPtLimit(0),mdvin2->GetPtLimit(1)};
1289   Float_t loosercuts[nvar1], tightercuts[nvar1];
1290   TString axistitles[nvar1];
1291   Int_t ncells[nvar1];
1292
1293   for (Int_t ivar=0;ivar<nvar1;ivar++){
1294     loosercuts[ivar]=mdvin1->GetCutValue(ivar,0);
1295     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));
1296     tightercuts[ivar]=mdvin1->GetCutValue(ivar,mdvin1->GetNCutSteps(ivar)-1);
1297     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)));
1298     axistitles[ivar]=mdvin1->GetAxisTitle(ivar);
1299     cout<<axistitles[ivar]<<"\t";
1300     ncells[ivar]=mdvin1->GetNCutSteps(ivar);
1301   }
1302
1303   AliMultiDimVector* mdvout= new AliMultiDimVector(Form("multiDimVectorPtBin%d",b1),"MultiDimVector",1,newptbins,mdvin1->GetNVariables(),ncells,loosercuts,tightercuts,axistitles);
1304   cout<<"Info: writing mdv"<<endl;
1305   listmdvout->Add(mdvout);
1306
1307   Bool_t isMC=kFALSE;
1308   TH1F* htestIsMC=(TH1F*)histlist->FindObject("hSgn_0");
1309   if(htestIsMC) isMC=kTRUE;
1310   Int_t nptbins=listamdv->GetEntries();
1311   Int_t nhist=(histlist->GetEntries()-1);//-1 because of fHistNevents
1312   if(isMC) nhist/=4; ///4 because hMass_, hSgn_,hBkg_,hRfl_
1313
1314  cout<<"Merging bin from "<<mdvin1->GetPtLimit(0)<<" to "<<mdvin1->GetPtLimit(1)<<" and from "<<mdvin2->GetPtLimit(0)<<" to "<<mdvin2->GetPtLimit(1)<<endl;
1315  Int_t firsth1=b1*ntotHperbin,firsth2=b2*ntotHperbin; //firsth2 = (b1+1)*ntotHperbin
1316  Int_t lasth1=firsth1+ntotHperbin-1,lasth2=firsth2+ntotHperbin-1;
1317  cout<<"Histograms from "<<firsth1<<" to "<<lasth1<<" and "<<firsth2<<" to "<<lasth2<<endl;
1318
1319  //add the others mdv to the list
1320  Int_t cnt=0;
1321  for(Int_t i=0;i<nptbins;i++){
1322    if(i!=b1 && i!=b2){
1323      AliMultiDimVector* vcttmp=(AliMultiDimVector*)listamdv->FindObject(Form("multiDimVectorPtBin%d",i));
1324      if(i>b2) {
1325        vcttmp->SetName(Form("multiDimVectorPtBin%d",b2+cnt));
1326        cnt++;
1327      }
1328      listmdvout->Add(vcttmp);
1329    }
1330  }
1331
1332  histlistout->Add((TH1F*)histlist->FindObject("fHistNEvents"));
1333
1334  Int_t ih2=firsth2;
1335
1336  for(Int_t ih1=firsth1;ih1<lasth1;ih1++){
1337    TH1F* h1=(TH1F*)histlist->FindObject(Form("hMass_%d",ih1));
1338    if(!h1){
1339      cout<<"hMass_"<<ih1<<" not found!"<<endl;
1340      continue;
1341    }
1342    TH1F* h2=(TH1F*)histlist->FindObject(Form("hMass_%d",ih2));
1343    if(!h2){
1344      cout<<"hMass_"<<ih2<<" not found!"<<endl;
1345      continue;
1346    }
1347    //h1->SetName(Form("hMass_%d",cnt));
1348    h1->Add(h2);
1349    histlistout->Add(h1);
1350    ih2++;
1351    //cnt++;
1352    h1=0x0;
1353  }
1354
1355  cnt=0;
1356  for(Int_t j=0;j<ntotHperbin*nptbins;j++){
1357    if(!(j>=firsth1 && j<lasth2)){
1358      TH1F* htmp=(TH1F*)histlist->FindObject(Form("hMass_%d",j));
1359      if(j>=lasth2){
1360        //cout<<lasth1<<" + "<<cnt<<endl;
1361        htmp->SetName(Form("hMass_%d",lasth1+cnt));
1362        cnt++;
1363      }
1364      histlistout->Add(htmp);
1365    }
1366  }
1367
1368  fout->cd();
1369  ((TDirectoryFile*)fout->Get(dirname))->cd();
1370  listmdvout->Write(mdvlistname.Data(),TObject::kSingleKey);
1371  histlistout->Write(listname.Data(),TObject::kSingleKey);
1372  fout->Close();
1373 }
1374
1375 void SubtractBkg(Int_t nhisto){
1376
1377   gStyle->SetFrameBorderMode(0);
1378   gStyle->SetCanvasColor(0);
1379   gStyle->SetFrameFillColor(0);
1380   gStyle->SetOptStat(0);
1381
1382   TString fitType="Exp";
1383   TString filename=Form("h%d%sMassFit.root",nhisto,fitType.Data());
1384
1385   TFile* fin=new TFile(filename.Data());
1386   if(!fin->IsOpen()){
1387     cout<<filename.Data()<<" not found, exit"<<endl;
1388     return;
1389   }
1390
1391   TKey* key=((TKey*)((TList*)fin->GetListOfKeys())->At(fin->GetNkeys()-1));
1392   TCanvas* canvas=((TCanvas*)fin->Get(key->GetName()));
1393   if(!canvas){
1394     cout<<"Canvas not found"<<endl;
1395     return;
1396   }
1397   canvas->Draw();
1398
1399   TH1F* hfit=(TH1F*)canvas->FindObject("fhistoInvMass");
1400   if(!hfit){
1401     canvas->ls();
1402     cout<<"Histogram not found"<<endl;
1403     return;
1404   }
1405
1406   TF1* funcBkgRecalc=(TF1*)hfit->FindObject("funcbkgRecalc");
1407   if(!funcBkgRecalc){
1408     cout<<"Background fit function (final) not found"<<endl;
1409     return;
1410   }
1411
1412   TF1* funcBkgFullRange=(TF1*)hfit->FindObject("funcbkgFullRange");
1413   if(!funcBkgFullRange){
1414     cout<<"Background fit function (side bands) not found"<<endl;
1415     return;
1416   }
1417
1418   Int_t nbins=hfit->GetNbinsX();
1419   Double_t min=hfit->GetBinLowEdge(1), width=hfit->GetBinWidth(1);
1420   TH1F* hsubRecalc=(TH1F*)hfit->Clone("hsub");
1421   hsubRecalc->SetMarkerColor(kRed);
1422   hsubRecalc->SetLineColor(kRed);
1423   hsubRecalc->GetListOfFunctions()->Delete();
1424   TH1F* hsubFullRange=(TH1F*)hfit->Clone("hsub");
1425   hsubFullRange->SetMarkerColor(kGray+2);
1426   hsubFullRange->SetLineColor(kGray+2);
1427   hsubFullRange->GetListOfFunctions()->Delete();
1428   for(Int_t i=0;i<nbins;i++){
1429     //Double_t x=min+i*0.5*width;
1430     Double_t x1=min+i*width, x2=min+(i+1)*width;
1431     Double_t ycont=hfit->GetBinContent(i+1);
1432     Double_t y1=funcBkgRecalc->Integral(x1,x2) / width;//funcBkgRecalc->Eval(x);
1433     hsubRecalc->SetBinContent(i+1,ycont-y1);
1434     Double_t y2=funcBkgFullRange->Integral(x1,x2) / width;//funcBkgFullRange->Eval(x);
1435     hsubFullRange->SetBinContent(i+1,ycont-y2);
1436   }
1437
1438   TCanvas* c=new TCanvas("c","subtraction");
1439   c->cd();
1440   hsubRecalc->DrawClone();
1441   hsubFullRange->DrawClone("sames");
1442
1443   for(Int_t i=0;i<nbins;i++){
1444     if(hsubRecalc->GetBinContent(i+1)<0) hsubRecalc->SetBinContent(i+1,0);
1445     if(hsubFullRange->GetBinContent(i+1)<0) hsubFullRange->SetBinContent(i+1,0);
1446   }
1447
1448   TCanvas *cvnewfits=new TCanvas("cvnewfits", "new Fits",1200,600);
1449   cvnewfits->Divide(2,1);
1450
1451   AliHFMassFitter fitter1(hsubRecalc,min,min+nbins*width,1,1);
1452   fitter1.MassFitter(kFALSE);
1453   fitter1.DrawHere(cvnewfits->cd(1));
1454
1455   AliHFMassFitter fitter2(hsubFullRange,min,min+nbins*width,1,1);
1456   fitter2.MassFitter(kFALSE);
1457   fitter2.DrawHere(cvnewfits->cd(2));
1458
1459   canvas->SaveAs(Form("h%d%sMassFit.png",nhisto,fitType.Data()));
1460   c->SaveAs(Form("h%d%sSubtr.png",nhisto,fitType.Data()));
1461   cvnewfits->SaveAs(Form("h%d%sFitNew.png",nhisto,fitType.Data()));
1462 }
1463