AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
charmCutsOptimization.C
Go to the documentation of this file.
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"
25 
26 #include <fstream>
27 
28 //global variables
33 Double_t sigma=0.0005;
36 
37 Double_t sigmaCut=0.035;//0.03;
38 Double_t errSgnCut=0.4;//0.35;
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;
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);
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;
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 
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 
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 
821  AliMultiDimVector* mvpur=cal->CalculatePurity();
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 
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 
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 
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;
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 
Int_t pdg
Float_t GetMaxSignificance(Int_t *cutIndices, Int_t ptbin) const
const char * filename
Definition: TestFCM.C:1
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1) const
write the canvas in a root file
double Double_t
Definition: External.C:58
Double_t PowerExpoBkgWoPeak(Double_t *x, Double_t *par)
Int_t GetNFinalPars() const
Int_t * GetRangeHistFromValues(AliMultiDimVector *vct, Int_t ptbin, Bool_t *valsgiven, Float_t *values, Int_t &nhistinrange)
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
TCanvas * canvas
Definition: DrawAnaELoss.C:28
TH1D * hSigma
Double_t PowerBkgWoPeak(Double_t *x, Double_t *par)
void SetInitialGaussianMean(Double_t mean)
Double_t sigma
Double_t mass
AliMultiDimVector * CalculatePurity() const
TSystem * gSystem
Int_t GetNHistFromValues(AliMultiDimVector *vct, Int_t ptbin, Float_t *values)
AliMultiDimVector * GetSignificanceError() const
Bool_t charmCutsOptimization(Bool_t isData=kTRUE, TString part="both", TString centr="no", Bool_t writefit=kTRUE, Int_t minentries=50, Double_t *rangefit=0x0, Bool_t useBinCounting=kTRUE)
TCanvas * c
Definition: TestFitELoss.C:172
Int_t GetNVariables() const
AliMultiDimVector * CalculatePurityError() const
Int_t GetNHistFromIndices(AliMultiDimVector *vct, Int_t ptbin, Int_t *indices)
ofstream outdetail
void Merge2Bins(Int_t b1, Int_t b2, TString pathin="./", Int_t decCh=2, TString part="both")
Int_t fitbtype
void SetElement(ULong64_t globadd, Float_t val)
virtual Bool_t RefitWithBkgOnly(Bool_t draw=kTRUE)
ofstream outcheck
ULong64_t GetGlobalAddressFromValues(const Float_t *values, Int_t ptbin) const
Float_t GetPtLimit(Int_t i) const
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
return total integral of the histogram
Bool_t GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Double_t GetMean() const
Int_t decCh
AliMultiDimVector * CalculateSOverBError() const
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)
Float_t GetCutValue(Int_t iVar, Int_t iCell) const
virtual void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
Double_t ExpoBkgWoPeak(Double_t *x, Double_t *par)
ULong64_t GetGlobalAddressFromIndices(const Int_t *ind, Int_t ptbin) const
Double_t nsigma
AliMultiDimVector * CalculateSOverB() const
void DrawPossibilities(Int_t ptbin, Bool_t *valsgiven, Float_t *values, TString path="./", Int_t decCh=2)
TList * fitter
Definition: DrawAnaELoss.C:26
void SetInitialGaussianSigma(Double_t sigma)
change the default value of the mean
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)
void DrawSigmas(TH2F *h2cuts)
Bool_t isMC
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)
Double_t errSgnCut
virtual Bool_t MassFitter(Bool_t draw=kTRUE)
ULong64_t GetNTotCells() const
Int_t GetNCutSteps(Int_t iVar) const
Double_t nSigmaMeanCut
Int_t rebin
TString GetAxisTitle(Int_t iVar) const
Float_t * GetCutValuesFromNHist(AliMultiDimVector *vct, Int_t ptbin, Int_t nhist)
const Int_t nbins
Double_t GetSigma() const
bool Bool_t
Definition: External.C:53
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
Int_t nptbins
Float_t GetElement(ULong64_t globadd) const
void SubtractBkg(Int_t nhisto)
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)
Bool_t reject
Double_t sigmaCut
virtual void WriteCanvas(TString userIDstring="", TString path="./", Double_t nsigma=3, Int_t writeFitInfo=1, Bool_t draw=kFALSE) const
write the TNtuple