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