PWG3/vertexingHF/macros/charmCutsOptimization.C

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