PWGHF/vertexingHF/macros/FitMassSpectra.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TInterpreter.h>
   3 #include <TString.h>
   4 #include <TObjString.h>
   5 #include <TObjArray.h>
   6 #include <TMath.h>
   7 #include <TFile.h>
   8 #include <TCanvas.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH1D.h>
  12 #include <TF1.h>
  13 #include <TStyle.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TDatabasePDG.h>
  17 #include <TGraph.h>
  18
  19 #include "AliAODRecoDecayHF.h"
  20 #include "AliRDHFCuts.h"
  21 #include "AliRDHFCutsDplustoKpipi.h"
  22 #include "AliRDHFCutsDStartoKpipi.h"
  23 #include "AliRDHFCutsDstoKKpi.h"
  24 #include "AliRDHFCutsD0toKpi.h"
  25 #include "AliHFMassFitter.h"
  26 #include "AliNormalizationCounter.h"
  27 #endif
  28
  29
  30 // MACRO to perform fits to D meson invariant mass spectra
  31 // and store raw yields and cut object into a root output file
  32 // Origin: F. Prino (prino@to.infn.it)
  33 // D0: C. Bianchin (cbianchi@pd.infn.it)
  34
  35
  36
  37 //
  38 enum {kD0toKpi, kDplusKpipi, kDStarD0pi, kDsKKpi};
  39 enum {kBoth, kParticleOnly, kAntiParticleOnly};
  40 enum {kExpo=0, kThrExpo=5, kLinear, kPol2};
  41 enum {kGaus=0, kDoubleGaus};
  42
  43
  44 // Common variables: to be configured by the user
  45 const Int_t nPtBins=6;
  46 Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.};
  47 Int_t rebin[nPtBins]={2,4,4,4,4,4};
  48 Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo
  49
  50 TString suffix="Loose_SecondSet1236_ForCF08";
  51
  52
  53 //const Int_t nPtBins=7;//6;
  54 //Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.};
  55 //Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts
  56 //Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts
  57 Int_t typeb=kExpo;
  58 Int_t types=kGaus;
  59 Int_t optPartAntiPart=kBoth;
  60 Int_t factor4refl=0;
  61 Double_t minMassForFit=1.7;
  62 Double_t maxMassForFit=2.1;
  63 Float_t massRangeForCounting=0.05; // GeV
  64 TH2F* hPtMass=0x0;
  65 Double_t nEventsForNorm=0.;
  66
  67 //for D0only
  68 Bool_t cutsappliedondistr=kTRUE;
  69 const Int_t nsamples=2;//3;
  70 //Int_t nevents[nsamples]={8.5635859e+07 /*LHC10b+cpart*/,8.9700624e+07/*LHC10dpart*/};
  71 //Int_t nevents[nsamples]={9.0374946e+07 /*LHC10b+c*/,9.7593785e+07/*LHC10d*/};
  72 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/};
  73 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/,9.0374946e+07 /*LHC10b+c*/};
  74 Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};
  75 //Int_t nevents[nsamples]={1. /*LHC10dnewTPCpid*/,1 /*LHC10dnewTPCpidrescale*/};
  76 // Functions
  77
  78 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass);
  79 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass);
  80 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass);
  81 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass);
  82
  83 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);
  84
  85
  86
  87
  88 void FitMassSpectra(Int_t analysisType=kDplusKpipi,
  89                     TString fileNameb="",
  90                     TString fileNamec="",
  91                     TString fileNamed="",
  92                     TString fileNamee=""
  93                ){
  94   //
  95
  96   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
  97   gStyle->SetOptTitle(0);
  98
  99   TObjArray* listFiles=new TObjArray();
 100   if(fileNameb!="") listFiles->AddLast(new TObjString(fileNameb.Data()));
 101   if(fileNamec!="") listFiles->AddLast(new TObjString(fileNamec.Data()));
 102   if(fileNamed!="") listFiles->AddLast(new TObjString(fileNamed.Data()));
 103   if(fileNamee!="") listFiles->AddLast(new TObjString(fileNamee.Data()));
 104   if(listFiles->GetEntries()==0){
 105     printf("Missing file names in input\n");
 106     return;
 107   }
 108
 109   TH1F** hmass=new TH1F*[nPtBins];
 110   for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0;
 111
 112   Float_t massD, massD0_fDstar;
 113   Bool_t retCode;
 114   if(analysisType==kD0toKpi){
 115     retCode=LoadD0toKpiHistos(listFiles,hmass);
 116     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 117   }
 118   else if(analysisType==kDplusKpipi){
 119     retCode=LoadDplusHistos(listFiles,hmass);
 120     massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 121   }
 122  else if(analysisType==kDStarD0pi){
 123     retCode=LoadDstarD0piHistos(listFiles,hmass);
 124     massD=TDatabasePDG::Instance()->GetParticle(413)->Mass();
 125     massD0_fDstar=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 126     massD =massD-massD0_fDstar;
 127   }
 128   else if(analysisType==kDsKKpi){
 129     retCode=LoadDsHistos(listFiles,hmass);
 130     massD=TDatabasePDG::Instance()->GetParticle(431)->Mass();
 131   }
 132   else{
 133     printf("Wrong analysis type parameter\n");
 134     return;
 135   }
 136   if(!retCode){
 137     printf("ERROR in reading input files\n");
 138     return;
 139   }
 140
 141
 142
 143   TH1D* hCntSig1=new TH1D("hCntSig1","hCntSig1",nPtBins,ptlims);
 144   TH1D* hCntSig2=new TH1D("hCntSig2","hCntSig2",nPtBins,ptlims);
 145   TH1D* hNDiffCntSig1=new TH1D("hNDiffCntSig1","hNDiffCntSig1",nPtBins,ptlims);
 146   TH1D* hNDiffCntSig2=new TH1D("hNDiffCntSig2","hNDiffCntSig2",nPtBins,ptlims);
 147   TH1D* hSignal=new TH1D("hSignal","hSignal",nPtBins,ptlims);
 148   TH1D* hRelErrSig=new TH1D("hRelErrSig","hRelErrSig",nPtBins,ptlims);
 149   TH1D* hInvSignif=new TH1D("hInvSignif","hInvSignif",nPtBins,ptlims);
 150   TH1D* hBackground=new TH1D("hBackground","hBackground",nPtBins,ptlims);
 151   TH1D* hBackgroundNormSigma=new TH1D("hBackgroundNormSigma","hBackgroundNormSigma",nPtBins,ptlims);
 152   TH1D* hSignificance=new TH1D("hSignificance","hSignificance",nPtBins,ptlims);
 153   TH1D* hMass=new TH1D("hMass","hMass",nPtBins,ptlims);
 154   TH1D* hSigma=new TH1D("hSigma","hSigma",nPtBins,ptlims);
 155
 156
 157   Int_t nMassBins=hmass[0]->GetNbinsX();
 158   Double_t hmin=TMath::Max(minMassForFit,hmass[0]->GetBinLowEdge(2));
 159   Double_t hmax=TMath::Min(maxMassForFit,hmass[0]->GetBinLowEdge(nMassBins-2)+hmass[0]->GetBinWidth(nMassBins-2));
 160   Float_t minBinSum=hmass[0]->FindBin(massD-massRangeForCounting);
 161   Float_t maxBinSum=hmass[0]->FindBin(massD+massRangeForCounting);
 162   Int_t iPad=1;
 163
 164   TF1* funBckStore1=0x0;
 165   TF1* funBckStore2=0x0;
 166   TF1* funBckStore3=0x0;
 167
 168   AliHFMassFitter** fitter=new AliHFMassFitter*[nPtBins];
 169   Double_t arrchisquare[nPtBins];
 170   TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800);
 171   Int_t nx=3, ny=2;
 172   if(nPtBins>6){
 173     nx=4;
 174     ny=3;
 175   }
 176   if(nPtBins>12){
 177     nx=5;
 178     ny=4;
 179   }
 180   c1->Divide(nx,ny);
 181
 182   Double_t sig,errsig,s,errs,b,errb;
 183   for(Int_t iBin=0; iBin<nPtBins; iBin++){
 184     c1->cd(iPad++);
 185     Int_t origNbins=hmass[iBin]->GetNbinsX();
 186     TH1F* hRebinned=RebinHisto(hmass[iBin],rebin[iBin],firstUsedBin[iBin]);
 187     hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
 188     hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
 189     fitter[iBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
 190     rebin[iBin]=origNbins/fitter[iBin]->GetBinN();
 191     fitter[iBin]->SetReflectionSigmaFactor(factor4refl);
 192     fitter[iBin]->SetInitialGaussianMean(massD);
 193     if(analysisType==kDStarD0pi) fitter[iBin]->SetInitialGaussianSigma(0.00065);
 194     Bool_t out=fitter[iBin]->MassFitter(0);
 195     if(!out) {
 196       fitter[iBin]->GetHistoClone()->Draw();
 197       continue;
 198     }
 199     Double_t mass=fitter[iBin]->GetMean();
 200     Double_t sigma=fitter[iBin]->GetSigma();
 201     arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare();
 202     TF1* fB1=fitter[iBin]->GetBackgroundFullRangeFunc();
 203     TF1* fB2=fitter[iBin]->GetBackgroundRecalcFunc();
 204     TF1* fM=fitter[iBin]->GetMassFunc();
 205     if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
 206     if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
 207     if(iBin==0 && fM) funBckStore3=new TF1(*fM);
 208
 209     fitter[iBin]->DrawHere(gPad);
 210     fitter[iBin]->Signal(3,s,errs);
 211     fitter[iBin]->Background(3,b,errb);
 212     fitter[iBin]->Significance(3,sig,errsig);
 213     Double_t ry=fitter[iBin]->GetRawYield();
 214     Double_t ery=fitter[iBin]->GetRawYieldError();
 215     Float_t cntSig1=0.;
 216     Float_t cntSig2=0.;
 217     Float_t cntErr=0.;
 218     for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
 219       Float_t bkg1=fB1 ? fB1->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
 220       Float_t bkg2=fB2 ? fB2->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
 221       cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1);
 222       cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2);
 223       cntErr+=(hmass[iBin]->GetBinContent(iMB));
 224     }
 225     hCntSig1->SetBinContent(iBin+1,cntSig1);
 226     hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr));
 227     hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s);
 228     hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
 229     hCntSig2->SetBinContent(iBin+1,cntSig2);
 230     hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s);
 231     hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
 232     hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr));
 233     hSignal->SetBinContent(iBin+1,ry);
 234     hSignal->SetBinError(iBin+1,ery);
 235     hRelErrSig->SetBinContent(iBin+1,errs/s);
 236     hInvSignif->SetBinContent(iBin+1,1/sig);
 237     hInvSignif->SetBinError(iBin+1,errsig/(sig*sig));
 238     hBackground->SetBinContent(iBin+1,b); //consider sigma
 239     hBackground->SetBinError(iBin+1,errb);
 240     hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma
 241     hBackgroundNormSigma->SetBinError(iBin+1,errb);
 242     hSignificance->SetBinContent(iBin+1,sig);
 243     hSignificance->SetBinError(iBin+1,errsig);
 244     hMass->SetBinContent(iBin+1,mass);
 245     hMass->SetBinError(iBin+1,0.0001);
 246     hSigma->SetBinContent(iBin+1,sigma);
 247     hSigma->SetBinError(iBin+1,0.0001);
 248
 249   }
 250
 251   /*
 252   c1->cd(1); // is some cases the fitting function of 1st bin get lost
 253   funBckStore1->Draw("same");
 254   funBckStore2->Draw("same");
 255   funBckStore3->Draw("same");
 256   */
 257
 258   TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
 259   cpar->Divide(2,1);
 260   cpar->cd(1);
 261   hMass->SetMarkerStyle(20);
 262   hMass->Draw("PE");
 263   hMass->GetXaxis()->SetTitle("Pt (GeV/c)");
 264   hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
 265   cpar->cd(2);
 266   hSigma->SetMarkerStyle(20);
 267   hSigma->Draw("PE");
 268   hSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
 269   hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");
 270
 271   TCanvas* csig=new TCanvas("csig","Results",1200,600);
 272   csig->Divide(3,1);
 273   csig->cd(1);
 274   hSignal->SetMarkerStyle(20);
 275   hSignal->SetMarkerColor(4);
 276   hSignal->SetLineColor(4);
 277   hSignal->GetXaxis()->SetTitle("Pt (GeV/c)");
 278   hSignal->GetYaxis()->SetTitle("Signal");
 279   hSignal->Draw("P");
 280   hCntSig1->SetMarkerStyle(26);
 281   hCntSig1->SetMarkerColor(2);
 282   hCntSig1->SetLineColor(2);
 283   hCntSig1->Draw("PSAME");
 284   hCntSig2->SetMarkerStyle(29);
 285   hCntSig2->SetMarkerColor(kGray+1);
 286   hCntSig2->SetLineColor(kGray+1);
 287   hCntSig2->Draw("PSAME");
 288   TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
 289   leg->SetFillColor(0);
 290   TLegendEntry* ent=leg->AddEntry(hSignal,"From Fit","PL");
 291   ent->SetTextColor(hSignal->GetMarkerColor());
 292   ent=leg->AddEntry(hCntSig1,"From Counting1","PL");
 293   ent->SetTextColor(hCntSig1->GetMarkerColor());
 294   ent=leg->AddEntry(hCntSig2,"From Counting2","PL");
 295   ent->SetTextColor(hCntSig2->GetMarkerColor());
 296   leg->Draw();
 297   csig->cd(2);
 298   hBackground->SetMarkerStyle(20);
 299   hBackground->Draw("P");
 300   hBackground->GetXaxis()->SetTitle("Pt (GeV/c)");
 301   hBackground->GetYaxis()->SetTitle("Background");
 302   csig->cd(3);
 303   hSignificance->SetMarkerStyle(20);
 304   hSignificance->Draw("P");
 305   hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)");
 306   hSignificance->GetYaxis()->SetTitle("Significance");
 307
 308   TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
 309   cDiffS->Divide(2,1);
 310   cDiffS->cd(1);
 311   hRelErrSig->SetMarkerStyle(20); //fullcircle
 312   hRelErrSig->SetTitleOffset(1.2);
 313   hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
 314   hRelErrSig->Draw("P");
 315   hInvSignif->SetMarkerStyle(21); //fullsquare
 316   hInvSignif->SetMarkerColor(kMagenta+1);
 317   hInvSignif->SetLineColor(kMagenta+1);
 318   hInvSignif->Draw("PSAMES");
 319   TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
 320   leg2->SetFillColor(0);
 321   TLegendEntry* ent2=leg2->AddEntry(hRelErrSig,"From Fit","P");
 322   ent2->SetTextColor(hRelErrSig->GetMarkerColor());
 323   ent2=leg2->AddEntry(hInvSignif,"1/Significance","PL");
 324   ent2->SetTextColor(hInvSignif->GetMarkerColor());
 325   leg2->Draw();
 326
 327   cDiffS->cd(2);
 328   hNDiffCntSig1->SetMarkerStyle(26);
 329   hNDiffCntSig1->SetMarkerColor(2);
 330   hNDiffCntSig1->SetLineColor(2);
 331   hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
 332   hNDiffCntSig1->Draw("P");
 333   hNDiffCntSig2->SetMarkerStyle(29);
 334   hNDiffCntSig2->SetMarkerColor(kGray+1);
 335   hNDiffCntSig2->SetLineColor(kGray+1);
 336   hNDiffCntSig2->Draw("PSAME");
 337   TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
 338   leg1->SetFillColor(0);
 339   TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1,"From Counting1","PL");
 340   ent1->SetTextColor(hNDiffCntSig1->GetMarkerColor());
 341   ent1=leg1->AddEntry(hNDiffCntSig2,"From Counting2","PL");
 342   ent1->SetTextColor(hNDiffCntSig2->GetMarkerColor());
 343   leg1->Draw();
 344
 345   TGraph* grReducedChiSquare=new TGraph(nPtBins,ptlims,arrchisquare);
 346   grReducedChiSquare->SetName("grReducedChiSquare");
 347   grReducedChiSquare->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
 348   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
 349   cChi2->cd();
 350   grReducedChiSquare->SetMarkerStyle(21);
 351   grReducedChiSquare->Draw("AP");
 352
 353   TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
 354   cbkgNormSigma->cd();
 355   hBackgroundNormSigma->SetMarkerStyle(20);
 356   hBackgroundNormSigma->Draw("P");
 357   hBackgroundNormSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
 358   hBackgroundNormSigma->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
 359   hBackgroundNormSigma->Draw();
 360
 361
 362   TString partname="Both";
 363   if(optPartAntiPart==kParticleOnly) {
 364     if(analysisType==kD0toKpi) partname="D0";
 365     if(analysisType==kDplusKpipi) partname="Dplus";
 366     if(analysisType==kDsKKpi) partname="Dsplus";
 367   }
 368   if(optPartAntiPart==kAntiParticleOnly) {
 369     if(analysisType==kD0toKpi) partname="D0bar";
 370     if(analysisType==kDplusKpipi) partname="Dminus";
 371     if(analysisType==kDsKKpi) partname="Dsminus";
 372   }
 373
 374   printf("Events for norm = %f\n",nEventsForNorm);
 375   TH1F* hNEvents=new TH1F("hNEvents","",1,0.,1.);
 376   hNEvents->SetBinContent(1,nEventsForNorm);
 377
 378   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update");
 379   outf->cd();
 380   hNEvents->Write();
 381   hMass->Write();
 382   hSigma->Write();
 383   hCntSig1->Write();
 384   hCntSig2->Write();
 385   hNDiffCntSig1->Write();
 386   hNDiffCntSig2->Write();
 387   hSignal->Write();
 388   hRelErrSig->Write();
 389   hInvSignif->Write();
 390   hBackground->Write();
 391   hBackgroundNormSigma->Write();
 392   hSignificance->Write();
 393   grReducedChiSquare->Write();
 394   hPtMass->Write();
 395   outf->Close();
 396 }
 397
 398
 399 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){
 400
 401   Int_t nFiles=listFiles->GetEntries();
 402   TList **hlist=new TList*[nFiles];
 403   AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
 404
 405   Int_t nReadFiles=0;
 406   for(Int_t iFile=0; iFile<nFiles; iFile++){
 407     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 408     TFile *f=TFile::Open(fName.Data());
 409     if(!f){
 410       printf("ERROR: file %s does not exist\n",fName.Data());
 411       continue;
 412     }
 413     printf("Open File %s\n",f->GetName());
 414     TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus");
 415     if(!dir){
 416       printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data());
 417       continue;
 418     }
 419     hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
 420     TList *listcut = (TList*)dir->Get(Form("coutputDplusCuts%s",suffix.Data()));
 421     dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
 422     if(nReadFiles>0){
 423       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 424       if(!sameCuts){
 425         printf("ERROR: Cut objects do not match\n");
 426         return kFALSE;
 427       }
 428     }
 429     AliNormalizationCounter* c=(AliNormalizationCounter*)dir->Get(Form("coutputDplusNorm%s",suffix.Data()));
 430     printf("Events for normalization = %f\n",c->GetNEventsForNorm());
 431     nEventsForNorm+=c->GetNEventsForNorm();
 432     nReadFiles++;
 433   }
 434   if(nReadFiles<nFiles){
 435     printf("WARNING: not all requested files have been found\n");
 436   }
 437
 438   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 439   printf("Number of pt bins for cut object = %d\n",nPtBins);
 440   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 441   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 442
 443   Int_t iFinBin=0;
 444   for(Int_t i=0;i<nPtBinsCuts;i++){
 445     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 446     if(iFinBin>nPtBins) break;
 447     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 448        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 449       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 450         TString histoName;
 451         if(optPartAntiPart==kBoth) histoName.Form("hMassPt%dTC",i);
 452         else if(optPartAntiPart==kParticleOnly) histoName.Form("hMassPt%dTCPlus",i);
 453         else if(optPartAntiPart==kAntiParticleOnly) histoName.Form("hMassPt%dTCMinus",i);
 454         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
 455         if(!htemp){
 456           printf("ERROR: Histogram %s not found\n",histoName.Data());
 457           return kFALSE;
 458         }
 459         if(!hMass[iFinBin]){
 460           hMass[iFinBin]=new TH1F(*htemp);
 461         }else{
 462           hMass[iFinBin]->Add(htemp);
 463         }
 464       }
 465     }
 466   }
 467   TString partname="Both";
 468   if(optPartAntiPart==kParticleOnly) partname="Dplus";
 469   if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
 470
 471   for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 472     TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC");
 473     if(iFile==0){
 474       hPtMass=new TH2F(*htemp2);
 475     }else{
 476       hPtMass->Add(htemp2);
 477     }
 478   }
 479
 480   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 481   outf->cd();
 482   dcuts[0]->Write();
 483   outf->Close();
 484
 485   return kTRUE;
 486
 487 }
 488
 489
 490 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass){
 491
 492   Int_t nFiles=listFiles->GetEntries();
 493   TList **hlist=new TList*[nFiles];
 494   AliRDHFCutsDstoKKpi** dcuts=new AliRDHFCutsDstoKKpi*[nFiles];
 495
 496   Int_t nReadFiles=0;
 497   for(Int_t iFile=0; iFile<nFiles; iFile++){
 498     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 499     TFile *f=TFile::Open(fName.Data());
 500     if(!f){
 501       printf("ERROR: file %s does not exist\n",fName.Data());
 502       continue;
 503     }
 504     printf("Open File %s\n",f->GetName());
 505     TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDs");
 506     if(!dir){
 507       printf("ERROR: directory PWG3_D2H_InvMassDs not found in %s\n",fName.Data());
 508       continue;
 509     }
 510     hlist[nReadFiles]=(TList*)dir->Get("coutputDs");
 511     TList *listcut = (TList*)dir->Get("coutputDsCuts");
 512     dcuts[nReadFiles]=(AliRDHFCutsDstoKKpi*)listcut->At(0);
 513     cout<< dcuts[nReadFiles]<<endl;
 514     if(nReadFiles>0){
 515       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 516       if(!sameCuts){
 517         printf("ERROR: Cut objects do not match\n");
 518         return kFALSE;
 519       }
 520     }
 521     nReadFiles++;
 522   }
 523   if(nReadFiles<nFiles){
 524     printf("WARNING: not all requested files have been found\n");
 525   }
 526
 527   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 528   printf("Number of pt bins for cut object = %d\n",nPtBins);
 529   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 530   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 531
 532   Int_t iFinBin=0;
 533   for(Int_t i=0;i<nPtBinsCuts;i++){
 534     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 535     if(iFinBin>nPtBins) break;
 536     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 537        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 538       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 539         TString histoName;
 540         if(optPartAntiPart==kBoth) histoName.Form("hMassAllPt%dphi",i);
 541         else if(optPartAntiPart==kParticleOnly){
 542           printf("Particle/Antiparticle not yet enabled for Ds");
 543           histoName.Form("hMassAllPt%dphi",i);
 544         }
 545         else if(optPartAntiPart==kAntiParticleOnly){
 546           printf("Particle/Antiparticle not yet enabled for Ds");
 547           histoName.Form("hMassAllPt%dphi",i);
 548         }
 549         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
 550         if(!htemp){
 551           printf("ERROR: Histogram %s not found\n",histoName.Data());
 552           return kFALSE;
 553         }
 554         if(!hMass[iFinBin]){
 555           hMass[iFinBin]=new TH1F(*htemp);
 556         }else{
 557           hMass[iFinBin]->Add(htemp);
 558         }
 559       }
 560     }
 561   }
 562   TString partname="Both";
 563   if(optPartAntiPart==kParticleOnly) partname="Both";
 564   if(optPartAntiPart==kAntiParticleOnly) partname="Both";
 565
 566   for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 567     TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassPhi");
 568     if(iFile==0){
 569       hPtMass=new TH2F(*htemp2);
 570     }else{
 571       hPtMass->Add(htemp2);
 572     }
 573   }
 574
 575   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 576   outf->cd();
 577   dcuts[0]->Write();
 578   outf->Close();
 579
 580   return kTRUE;
 581
 582 }
 583
 584
 585 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){
 586   //
 587   Int_t nFiles=listFiles->GetEntries();
 588   TList **hlist=new TList*[nFiles];
 589   AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
 590
 591   Int_t nReadFiles=0;
 592   for(Int_t iFile=0; iFile<nFiles; iFile++){
 593     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 594     TFile *f=TFile::Open(fName.Data());
 595     if(!f){
 596       printf("ERROR: file %s does not exist\n",fName.Data());
 597       continue;
 598     }
 599     printf("Open File %s\n",f->GetName());
 600
 601     TString dirname="PWG3_D2H_D0InvMass";
 602     if(optPartAntiPart==kParticleOnly) dirname+="D0";
 603     else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
 604     if(cutsappliedondistr) dirname+="C";
 605     TDirectory *dir = (TDirectory*)f->Get(dirname);
 606     if(!dir){
 607       printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
 608       continue;
 609     }
 610     TString listmassname="coutputmassD0Mass";
 611     if(optPartAntiPart==kParticleOnly) listmassname+="D0";
 612     else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
 613     if(cutsappliedondistr) listmassname+="C";
 614
 615     hlist[nReadFiles]=(TList*)dir->Get(listmassname);
 616
 617     TString cutsobjname="cutsD0";
 618     if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
 619     else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
 620     if(cutsappliedondistr) cutsobjname+="C";
 621
 622     dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
 623     if(nReadFiles>0){
 624       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 625       if(!sameCuts){
 626         printf("ERROR: Cut objects do not match\n");
 627         return kFALSE;
 628       }
 629     }
 630     nReadFiles++;
 631   }
 632   if(nReadFiles<nFiles){
 633     printf("WARNING: not all requested files have been found\n");
 634     if (nReadFiles==0) {
 635       printf("ERROR: Any file/dir found\n");
 636       return kFALSE;
 637     }
 638   }
 639
 640   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 641   printf("Number of pt bins for cut object = %d\n",nPtBins);
 642   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 643   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 644
 645   Int_t iFinBin=0;
 646   for(Int_t i=0;i<nPtBinsCuts;i++){
 647     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 648     if(iFinBin>nPtBins) break;
 649     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 650        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 651       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 652         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histMass_%d",i));
 653         if(!hMass[iFinBin]){
 654           hMass[iFinBin]=new TH1F(*htemp);
 655         }else{
 656           hMass[iFinBin]->Add(htemp);
 657         }
 658       }
 659     }
 660   }
 661   TString partname="Both";
 662   if(optPartAntiPart==kParticleOnly) partname="D0";
 663   if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
 664
 665   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 666   outf->cd();
 667   dcuts[0]->Write();
 668   outf->Close();
 669   return kTRUE;
 670 }
 671
 672 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass){
 673   //
 674   Int_t nFiles=listFiles->GetEntries();
 675   TList **hlist=new TList*[nFiles];
 676   AliRDHFCutsDStartoKpipi** dcuts=new AliRDHFCutsDStartoKpipi*[nFiles];
 677   Int_t nReadFiles=0;
 678   for(Int_t iFile=0; iFile<nFiles; iFile++){
 679     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 680     TFile *f=TFile::Open(fName.Data());
 681     if(!f){
 682       printf("ERROR: file %s does not exist\n",fName.Data());
 683       continue;
 684     }
 685     printf("Open File %s\n",f->GetName());
 686     TString dirname="PWG3_D2H_DStarSpectra";
 687     TDirectory *dir = (TDirectory*)f->Get(dirname);
 688     if(!dir){
 689       printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
 690       continue;
 691     }
 692     TString listmassname="DStarPID00";
 693
 694     hlist[nReadFiles]=(TList*)dir->Get(listmassname);
 695     TString cutsobjname="DStartoKpipiCuts";
 696     dcuts[nReadFiles]=(AliRDHFCutsDStartoKpipi*)dir->Get(cutsobjname);
 697     if(nReadFiles>0){
 698       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 699       if(!sameCuts){
 700         printf("ERROR: Cut objects do not match\n");
 701         return kFALSE;
 702       }
 703     }
 704     nReadFiles++;
 705   }
 706   if(nReadFiles<nFiles){
 707     printf("WARNING: not all requested files have been found\n");
 708     if (nReadFiles==0) {
 709       printf("ERROR: Any file/dir found\n");
 710       return kFALSE;
 711     }
 712   }
 713   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 714   printf("Number of pt bins for cut object = %d\n",nPtBins);
 715   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 716   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 717
 718
 719
 720
 721   Int_t iFinBin=0;
 722   for(Int_t i=0;i<nPtBinsCuts;i++){
 723     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 724     if(iFinBin>nPtBins) break;
 725     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 726        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 727       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 728         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histDeltaMass_%d",i));
 729         if(!hMass[iFinBin]){
 730           hMass[iFinBin]=new TH1F(*htemp);
 731         }else{
 732           hMass[iFinBin]->Add(htemp);
 733         }
 734       }
 735     }
 736   }
 737   TString partname="Both";
 738   if(optPartAntiPart==kParticleOnly) partname="DStar";
 739   if(optPartAntiPart==kAntiParticleOnly) partname="DStarbar";
 740
 741   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 742   outf->cd();
 743   dcuts[0]->Write();
 744   outf->Close();
 745   return kTRUE;
 746 }
 747
 748 void CompareFitTypes(TString* paths, TString* legtext,Int_t ncmp=3,TString* filenameYield=0x0){
 749   //read ncmp RawYield.roots and draw them together
 750   //set the global variable nevents before running
 751   //arguments:
 752   // - paths= vector of ncmp dimension with the paths of the file RawYield.root
 753   // - legtext= vector of ncmp dimension with the label for the legend
 754   // - ncmp= number of files to compare (default is 3)
 755   // - filenameYield= optional ncmp-dimensional array with the difference between the names of the files to be compared (useful if the 2 files are in the same directory but have different names)
 756
 757   gStyle->SetOptStat(0);
 758   gStyle->SetOptTitle(0);
 759   gStyle->SetCanvasColor(0);
 760   gStyle->SetFrameFillColor(0);
 761   gStyle->SetTitleFillColor(0);
 762   gStyle->SetFrameBorderMode(0);
 763
 764   if(!filenameYield) filenameYield=new TString[ncmp];
 765
 766   for(Int_t k=0;k<ncmp;k++){
 767     if(!filenameYield) filenameYield[k]="RawYield.root";
 768     filenameYield[k].Prepend(paths[k]);
 769   }
 770
 771   TCanvas* cSig=new TCanvas("cSig","Results",1200,600);
 772   cSig->Divide(3,1);
 773   TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600);
 774   TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
 775   cDiffS->Divide(2,1);
 776   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
 777
 778   TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
 779   leg1->SetFillColor(0);
 780   TLegend* leg2=(TLegend*)leg1->Clone();
 781   TLegend* leg3=(TLegend*)leg1->Clone();
 782   TLegend* leg4=new TLegend(0.4,0.6,0.8,0.89);
 783   leg4->SetFillColor(0);
 784
 785   for(Int_t iTypes=0;iTypes<ncmp;iTypes++){
 786     TFile* fin=new TFile(filenameYield[iTypes]);
 787     if(!fin){
 788       printf("WARNING: %s not found",filenameYield[iTypes].Data());
 789       continue;
 790     }
 791
 792     TH1F* hSignal=(TH1F*)fin->Get("hSignal");
 793     TH1F* hBackground=(TH1F*)fin->Get("hBackground");
 794     TH1F* hBackgroundNormSigma=(TH1F*)fin->Get("hBackgroundNormSigma");
 795     TH1F* hSignificance=(TH1F*)fin->Get("hSignificance");
 796     hSignal->SetName(Form("%s%d",hSignal->GetName(),iTypes));
 797     hBackground->SetName(Form("%s%d",hBackground->GetName(),iTypes));
 798     hBackgroundNormSigma->SetName(Form("%s%d",hBackgroundNormSigma->GetName(),iTypes));
 799     hSignificance->SetName(Form("%s%d",hSignificance->GetName(),iTypes));
 800
 801     hSignal->SetMarkerColor(iTypes+2);
 802     hSignal->SetLineColor(iTypes+2);
 803     hBackground->SetMarkerColor(iTypes+2);
 804     hBackground->SetLineColor(iTypes+2);
 805     hBackgroundNormSigma->SetMarkerColor(iTypes+2);
 806     hBackgroundNormSigma->SetLineColor(iTypes+2);
 807     hSignificance->SetMarkerColor(iTypes+2);
 808     hSignificance->SetLineColor(iTypes+2);
 809
 810     TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL");
 811     ent4->SetTextColor(hSignal->GetMarkerColor());
 812
 813     if(ncmp==nsamples){
 814       printf("Info: Normalizing signal, background and significance to the number of events\n");
 815       hSignal->Scale(1./nevents[iTypes]);
 816       hBackground->Scale(1./nevents[iTypes]);
 817       hBackgroundNormSigma->Scale(1./nevents[iTypes]);
 818       hSignificance->Scale(1./TMath::Sqrt(nevents[iTypes]));
 819     }
 820
 821     if(iTypes==0){
 822       cSig->cd(1);
 823       hSignal->DrawClone("P");
 824       cSig->cd(2);
 825       hBackground->DrawClone("P");
 826       cSig->cd(3);
 827       hSignificance->DrawClone("P");
 828       cBkgN->cd();
 829       hBackgroundNormSigma->DrawClone("P");
 830     } else{
 831       cSig->cd(1);
 832       hSignal->DrawClone("Psames");
 833       cSig->cd(2);
 834       hBackground->DrawClone("Psames");
 835       cSig->cd(3);
 836       hSignificance->DrawClone("Psames");
 837       cBkgN->cd();
 838       hBackgroundNormSigma->DrawClone("Psames");
 839     }
 840
 841     TH1F* hRelErrSig=(TH1F*)fin->Get("hRelErrSig");
 842     TH1F* hInvSignif=(TH1F*)fin->Get("hInvSignif");
 843     hRelErrSig->SetName(Form("%s%d",hRelErrSig->GetName(),iTypes));
 844     hInvSignif->SetName(Form("%s%d",hInvSignif->GetName(),iTypes));
 845
 846     hRelErrSig->SetMarkerColor(iTypes+2);
 847     hRelErrSig->SetLineColor(iTypes+2);
 848     hInvSignif->SetMarkerColor(iTypes+2);
 849     hInvSignif->SetLineColor(iTypes+2);
 850
 851     TLegendEntry* ent1=leg1->AddEntry(hRelErrSig,Form("From Fit (%s)",legtext[iTypes].Data()),"P");
 852     ent1->SetTextColor(hRelErrSig->GetMarkerColor());
 853     ent1=leg1->AddEntry(hInvSignif,Form("1/Significance (%s)",legtext[iTypes].Data()),"PL");
 854     ent1->SetTextColor(hInvSignif->GetMarkerColor());
 855
 856     cDiffS->cd(1);
 857     if(iTypes==0){
 858       hRelErrSig->DrawClone("P");
 859       hInvSignif->DrawClone();
 860     } else{
 861       hRelErrSig->DrawClone("Psames");
 862       hInvSignif->DrawClone("sames");
 863     }
 864
 865     TH1F* hNDiffCntSig1=(TH1F*)fin->Get("hNDiffCntSig1");
 866     TH1F* hNDiffCntSig2=(TH1F*)fin->Get("hNDiffCntSig2");
 867     hNDiffCntSig1->SetName(Form("%s%d",hNDiffCntSig1->GetName(),iTypes));
 868     hNDiffCntSig2->SetName(Form("%s%d",hNDiffCntSig2->GetName(),iTypes));
 869
 870     hNDiffCntSig1->SetMarkerColor(iTypes+2);
 871     hNDiffCntSig1->SetLineColor(iTypes+2);
 872     hNDiffCntSig2->SetMarkerColor(iTypes+2);
 873     hNDiffCntSig2->SetLineColor(iTypes+2);
 874     TLegendEntry* ent2=leg2->AddEntry(hNDiffCntSig1,Form("From Counting1 (%s)",legtext[iTypes].Data()),"PL");
 875     ent2->SetTextColor(hNDiffCntSig1->GetMarkerColor());
 876     ent2=leg2->AddEntry(hNDiffCntSig2,Form("From Counting2 (%s)",legtext[iTypes].Data()),"PL");
 877     ent2->SetTextColor(hNDiffCntSig2->GetMarkerColor());
 878
 879     cDiffS->cd(2);
 880     if(iTypes==0){
 881       hNDiffCntSig1->DrawClone();
 882       hNDiffCntSig2->DrawClone();
 883     }else{
 884      hNDiffCntSig1->DrawClone("sames");
 885      hNDiffCntSig2->DrawClone("sames");
 886     }
 887
 888     TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare");
 889     grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes));
 890
 891     grReducedChiSquare->SetMarkerColor(iTypes+2);
 892     grReducedChiSquare->SetLineColor(iTypes+2);
 893     TLegendEntry* ent3=leg3->AddEntry(grReducedChiSquare,Form("%s",legtext[iTypes].Data()),"PL");
 894     ent3->SetTextColor(grReducedChiSquare->GetMarkerColor());
 895
 896     cChi2->cd();
 897     if(iTypes==0) grReducedChiSquare->DrawClone("AP");
 898     else grReducedChiSquare->DrawClone("P");
 899   }
 900
 901   cSig->cd(1);
 902   leg4->Draw();
 903
 904   cDiffS->cd(1);
 905   leg1->Draw();
 906
 907   cDiffS->cd(2);
 908   leg2->Draw();
 909
 910   cChi2->cd();
 911   leg3->Draw();
 912
 913   TFile* fout=new TFile("ComparisonRawYield.root","RECREATE");
 914   fout->cd();
 915   cDiffS->Write();
 916   cChi2->Write();
 917   fout->Close();
 918 }
 919
 920
 921 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){
 922   // Rebin histogram, from bin firstUse to lastUse
 923   // Use all bins if firstUse=-1
 924
 925   Int_t nBinOrig=hOrig->GetNbinsX();
 926   Int_t firstBinOrig=1;
 927   Int_t lastBinOrig=nBinOrig;
 928   Int_t nBinOrigUsed=nBinOrig;
 929   Int_t nBinFinal=nBinOrig/reb;
 930   if(firstUse>=1){
 931     firstBinOrig=firstUse;
 932     nBinFinal=(nBinOrig-firstUse+1)/reb;
 933     nBinOrigUsed=nBinFinal*reb;
 934     lastBinOrig=firstBinOrig+nBinOrigUsed-1;
 935   }else{
 936     Int_t exc=nBinOrigUsed%reb;
 937     if(exc!=0){
 938       nBinOrigUsed-=exc;
 939       firstBinOrig+=exc/2;
 940       lastBinOrig=firstBinOrig+nBinOrigUsed-1;
 941     }
 942   }
 943
 944   printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
 945   Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
 946   Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
 947   TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
 948   Int_t lastSummed=firstBinOrig-1;
 949   for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
 950     Float_t sum=0.;
 951     for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
 952       sum+=hOrig->GetBinContent(lastSummed+1);
 953       lastSummed++;
 954     }
 955     hRebin->SetBinContent(iBin,sum);
 956   }
 957   return hRebin;
 958 }