C++ TH1F::Clone方法代码示例

void beffAnalysis(const char* input, const char* output) {
  // input
  TFile* inputFile = TFile::Open(input);
  TTree* btagEff = (TTree*)inputFile->Get("btagEff");

  // output
  TFile* outputFile = TFile::Open(output,"RECREATE");

  // histogram with proper binning... cloned later on 
  Double_t binning[23] = {20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,120,140,160,180,200,1000};
  TH1F* ptSpectrum = new TH1F("PtSpectrum","PtSpectrum",22,binning);
  ptSpectrum->Sumw2();

  // produce the ratio plot for the 12 combinations of (CSVL,CSVM,CSVT),(Barrel,Endcap),(b,c,l)
  TClonesArray algorithms("TCut",3);
  new(algorithms[0]) TCut("CSVL","csv>0.244");
  new(algorithms[1]) TCut("CSVM","csv>0.679");
  new(algorithms[2]) TCut("CSVT","csv>0.898");

  TClonesArray etaRegions("TCut",2);
  new(etaRegions[0]) TCut("Barrel","abs(eta)<=1.2");
  new(etaRegions[1]) TCut("Endcaps","abs(eta)>1.2");

  TClonesArray flavor("TCut",3);
  new(flavor[0]) TCut("l","abs(flavor)!=4 && abs(flavor)!=5");
  new(flavor[1]) TCut("c","abs(flavor)==4");
  new(flavor[2]) TCut("b","abs(flavor)==5");

  for(int i=0; i< algorithms.GetEntries() ; ++i) {
    outputFile->mkdir(((TCut*)algorithms.At(i))->GetName());
    outputFile->cd(((TCut*)algorithms.At(i))->GetName());
    for(int j=0; j< etaRegions.GetEntries() ; ++j) {
      for(int k=0; k< flavor.GetEntries() ; ++k) {
        // histogram before tagging
        TH1F* pretag = ptSpectrum->Clone("pretag");
        btagEff->Draw("pt>>pretag",((*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // histogram after tagging
        TH1F* posttag = ptSpectrum->Clone("posttag");
        btagEff->Draw("pt>>posttag",((*(TCut*)algorithms.At(i))&&(*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // ratio
        TH1F* ratio = posttag->Clone(Form("h_eff_bTagOverGoodJet_pt%s_%s",((TCut*)flavor.At(k))->GetName(),
                                                                          ((TCut*)etaRegions.At(j))->GetName()));
        ratio->Divide(pretag);
        // cleanup
        delete pretag;
        delete posttag;
      }
    }
  }

  // cleanup
  algorithms.Delete();
  etaRegions.Delete();
  flavor.Delete();
  ptSpectrum->SetDirectory(0);
  outputFile->Write();
  outputFile->Close();
  inputFile->Close();
}

void applyEfficiencyToTauMuons() {

    TFile* _fPtSm = new TFile("smearedMuonPtFromTau.07.09.2013.root","read");
    TFile* _fEff = new TFile("mcWEffForTauMuonStudy_07_08_2013.root","read");
    TFile* _fPtSmEff = new TFile("smearedEffAppliedMuonPtFromTau.root","recreate");

    _fPtSm->cd();
    TH1F* hPtSm = (TH1F*)_fPtSm->Get("hSmearedMuonPt");
    _fEff->cd();
    TGraphAsymmErrors* grEff = (TGraphAsymmErrors*)_fEff->Get("pEffWAccCuts");

    double* yEff = grEff->GetY();
    TH1F*  hnew = (TH1F*)hPtSm->Clone("hnew");
    ///Apply efficiency bin by bin
    for(int igr=0; igr<grEff->GetN(); ++igr) {

        double ptNominal = hPtSm->GetBinContent(igr+1);
        std::cout << "pt nominal " << ptNominal << " eff " << yEff[igr] << " = " << yEff[igr]*ptNominal << std::endl;
        hnew->SetBinContent(igr+1,yEff[igr]*ptNominal);
    }

    _fPtSmEff->cd();
    hnew->Write();
    _fPtSmEff->Write();
    _fPtSmEff->Close();
}

void draw_lumi(){
	TString filenameData = "lumi_effData.root";
	TString filenameMC = "lumi_effMC.root";
	TString fileMC = "../newDoubleMuonTree/tnpZ_theTreeCleanedMC.root";
	//TString fileData = "doubleMuonTrees/tnpZ_theTreeCleaned.root";
        TString fileData = "/afs/cern.ch/user/q/quwang/work/Trigger/CMSSW_7_6_3_patch2/src/MuonAnalysis/TagAndProbe/test/zmumu/tnpZ_Data.root";//"doubleMuonTrees/tnpZ_theTreeCleaned.root";

	//TString fileMC = "newDoubleMuonTree/TnPtreeMC.root";
	//TString fileData = "doubleMuonTrees/tnpZ_Data_25ns_run2015D_doubleMuon.root";
	Eff( fileData, filenameData, "data_lumi");	
	//Eff( fileMC, filenameMC, "MC_lumi");	
	
	TFile * hData = new TFile("lumi_effData.root");
	TH1F * histD = (TH1F *)hData->Get("hlumi_eff");
	TFile * hMC = new TFile("lumi_effMC.root");
	TH1F * histM = (TH1F *)hMC->Get("hlumi_eff");
	TH1F * hlumi_DM = (TH1F *)histD->Clone("hlumi_DM");
	hlumi_DM->Sumw2();
	hlumi_DM->SetMinimum(0);
	hlumi_DM->Divide(histD, histM, 1, 1, "b");
	hlumi_DM->SetTitle("dZ eff Data/MC");
	hlumi_DM->GetXaxis()->SetTitle("lumi");
	hlumi_DM->GetYaxis()->SetTitle("eff Data/MC");

	TCanvas * C2 = new TCanvas();
	C2->Draw();
	hlumi_DM->Draw("colz");
	//C2->SaveAs("dZ_lumi_DM.pdf");

}

//--------------------------------------------------------------------------------------------------
// make fake rate histogram from numerator and denominators
//--------------------------------------------------------------------------------------------------
TH1F* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames,
                           int faketype, string histName) {

  string numeratorHistName = "";
  string denominatorHistName = "";
  string dirName = "";

  if (faketype == 11) {
    dirName = "ComputeElectronFakeRateMod";
    numeratorHistName = "hElectronNumeratorEt";
    denominatorHistName = "hElectronDenominatorEt";
  } else if (faketype == 13) {
    dirName = "ComputeMuonFakeRateMod";
    numeratorHistName = "hMuonNumeratorPt";
    denominatorHistName = "hMuonDenominatorPt";
  } else {
    cerr << "Error: faketype = " << faketype << " is not recognized." << endl;
    return 0;
  }

  TH1F *denominator = addAllSamples(datasetFiles, datasetNames, dirName, denominatorHistName);
  TH1F *numerator = addAllSamples(datasetFiles, datasetNames, dirName, numeratorHistName);

  //create fake rate hist from denominator binning
  TH1F *fakeRateHist = (TH1F*)numerator->Clone(histName.c_str());
  fakeRateHist->GetYaxis()->SetTitle("Fake Rate");
  //divide by the denominator to get fake rate;
  fakeRateHist->Divide(denominator);

  return fakeRateHist;
}

void HighlightZoom(TVirtualPad *pad, TObject *obj, Int_t xhb, Int_t yhb)
{
   TH1F *h = (TH1F *)obj;
   if(!h) return;

   TCanvas *c2 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("c2");
   static TH1 *hz = 0;
   if (!h->IsHighlight()) { // after highlight disabled
      if (c2) delete c2;
      if (hz) { delete hz; hz = 0; }
      return;
   }

   if (!c2) {
      c2 = new TCanvas("c2", "c2", 605, 0, 400, 400);
      c2->SetGrid();
      if (hz) hz->Draw(); // after reopen this canvas
   }
   if (!hz) {
      hz = (TH1 *)h->Clone("hz");
      hz->SetTitle(TString::Format("%s (zoomed)", hz->GetTitle()));
      hz->SetStats(kFALSE);
      hz->Draw();
      c2->Update();
      hz->SetHighlight(kFALSE);
   }

   Int_t zf = hz->GetNbinsX()*0.05; // zoom factor
   hz->GetXaxis()->SetRange(xhb-zf, xhb+zf);

   c2->Modified();
   c2->Update();
}

void TAGH_doubles_E()
{
    TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("TAGH_doubles");
    if(dir) dir->cd();

    TH1F *doubles =(TH1F*)gDirectory->Get("BeforeMergingDoubles/BM2_Energy");
    TH1F *total =(TH1F*)gDirectory->Get("BeforeMergingDoubles/BM1_Energy");

    TH1F *f_doubles = (TH1F*)doubles->Clone();
    f_doubles->Sumw2();
    f_doubles->Divide(doubles, total);

    if(gPad == NULL) {
        TCanvas *c1 = new TCanvas("c1","TAGH double-hit fraction",150,10,990,660);
        c1->cd(0);
        c1->Draw();
        c1->Update();
    }

    if(!gPad) return;
    TCanvas* c1 = gPad->GetCanvas();

    gStyle->SetOptStat("");
    f_doubles->SetTitle("TAGH double-hit fraction vs. energy");
    f_doubles->SetTitleSize(0.045, "XY");
    f_doubles->GetXaxis()->SetTitle("TAGH energy [GeV]");
    f_doubles->GetYaxis()->SetTitle("double-hit fraction");
    f_doubles->SetAxisRange(0.0,1.0,"Y");
    f_doubles->Draw();
}

/*
plotContainer* makeRatio (plotContainer& pcOS, plotContainer& pcSS,
                          vector<pair <TString, TCut> >& selections,
                          vector<string>& variablesList,
                          vector <string>& samples,
                          vector <float>& isoThr1,
                          vector <float>& isoThr2
                          )
*/
TH1F* makeRatio (plotContainer& pcOS, plotContainer& pcSS,
                 pair <TString, TCut>& selection,
                 string& varName,
                 string& sampleName,
                 float thr1,
                 float thr2
                )

{
  TString OS_selName = Form("OSaiso_%f_%f_", thr1, thr2);
  TString SS_selName = Form("SSaiso_%f_%f_", thr1, thr2);
  OS_selName = OS_selName + selection.first;
  SS_selName = SS_selName + selection.first;

  TH1F* hOS = pcOS.getHisto (varName, string(OS_selName.Data()), sampleName);
  TH1F* hSS = pcSS.getHisto (varName, string(SS_selName.Data()), sampleName);
  
  TString newName = "ratio_" + varName + "_" + selection.first + Form("_%f_%f", thr1, thr2) ;

  TH1F* hratio = (TH1F*)hOS->Clone(newName);
  hratio -> Divide (hSS);
  hratio -> SetMinimum(0.80);
  hratio -> SetMaximum(1.30);
  hratio -> GetYaxis()->SetTitle("OS/SS");

  return hratio;
}

TH1F* dist(TTree* treeElEl, TTree* treeMuMu, TTree* treeMuEl, const TString &variable, int nBin, float * bins, TCut precut, TCut precut_em, TCut fcut, const TString & name, const TString & tagName, int color, const TString & xtitle){

  TCut dencut_em =  precut_em + fcut;
  TCut dencut    =  precut    + fcut;

  TH1F* h_eff = new TH1F(Form("h_eff_%s",name.Data()),"h_eff",nBin,bins);

  TH1F* hden = new TH1F(Form("hden_%s",name.Data()),"hden",nBin,bins);
  TH1F* hden_ee = new TH1F(Form("hden_%s_ee",name.Data()),"hden_ee",nBin,bins);
  TH1F* hden_mm = new TH1F(Form("hden_%s_mm",name.Data()),"hden_mm",nBin,bins);
  TH1F* hden_em = new TH1F(Form("hden_%s_em",name.Data()),"hden_em",nBin,bins);

  treeElEl->Project(Form("hden_%s_ee",name.Data()),Form("%s",variable.Data()),dencut,"");
  treeMuMu->Project(Form("hden_%s_mm",name.Data()),Form("%s",variable.Data()),dencut,"");
  treeMuEl->Project(Form("hden_%s_em",name.Data()),Form("%s",variable.Data()),dencut_em,"");

  hden->Add(hden_ee,1);
  hden->Add(hden_mm,1);
  hden->Add(hden_em,1);

  //TCanvas*  tag_eff = new TCanvas(Form("eff_%s",name.Data()),Form("eff_%s",name.Data()),500,500);
  TH1F* h_effint = hden->Clone();
  h_effint->Scale(1/h_effint->Integral()); 
  h_effint->SetTitle(xtitle);

  h_effint->SetLineColor(color);
  h_effint->GetXaxis()->SetTitle(xtitle);
  h_effint->GetYaxis()->SetTitle("Normalized Entries");

  return h_effint;

}

void compnclusts(Int_t run)
{
  TFile* f = new TFile(Form("hodtest_%d.root",run));
  cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
  TH1F* h = nclust;

  TFile* f1 = new TFile(Form("%d_hbk.root",run));
  cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
  TH1F* h1;
  switch (run) {
  case 50017 :
  //    h1 = h212;   //A+
    break;
  default :
    h1 = h412;   //hnclusters
  }

  TCanvas *c1 = new TCanvas("c1", "Shower Cluster Map", 1000, 667);

  gPad->SetLogy();

  h1->SetFillColor(kGreen);
  h1->SetLineColor(kGreen);
  h1->SetFillStyle(1111);
  h1->Draw();

  h->SetFillColor(kBlue);
  h->SetLineWidth(2);
  h->SetFillStyle(0);
  h->Draw("same");

  TLatex l;
  l.SetTextSize(0.04);
  Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
  Float_t xmin = h1->GetXaxis()->GetXmin();
  Float_t xmax = h1->GetXaxis()->GetXmax();
  Float_t xt = xmin + 0.67*(xmax-xmin);

  l.SetTextColor(kGreen);
  l.DrawLatex(xt,0.095*maxy,"Engine");
  l.SetTextColor(kBlue);
  l.DrawLatex(xt,0.045*maxy,"hcana");

  // Difference between the histograms.

  TCanvas *c2 = new TCanvas("c2", "Cluster differences", 1000, 667);

  TH1F* dif = h->Clone();

  dif->Add(h,h1,1.,-1.);

  dif->SetTitle("Difference");
  dif->SetFillColor(kRed);
  dif->SetLineColor(kRed);
  dif->SetLineWidth(1);
  dif->SetFillStyle(1111);
  dif->Draw();

}

//_____________________________________________________________________________
Int_t ProofSimple::GetHistosFromFC(TCanvas *cv)
{
   // Check for the histograms in the files of a possible TFileCollection

   TIter nxo(fOutput);
   TFileCollection *fc = 0;
   Bool_t fc_found = kFALSE, hs_found = kFALSE;
   while ((fc = (TFileCollection *) nxo())) {
      if (strcmp(fc->ClassName(), "TFileCollection")) continue;
      fc_found = kTRUE;
      if (!fHist) {
         fHist = new TH1F*[fNhist];
         for (Int_t i = 0; i < fNhist; i++) { fHist[i] = 0; }
      } else {
         for (Int_t i = 0; i < fNhist; i++) { SafeDelete(fHist[i]); }
      }
      // Go through the list of files
      TIter nxf(fc->GetList());
      TFileInfo *fi = 0;
      while ((fi = (TFileInfo *) nxf())) {
         TFile *f = TFile::Open(fi->GetCurrentUrl()->GetUrl());
         if (f) {
            for (Int_t i = 0; i < fNhist; i++) {
               TString hn = TString::Format("h%d", i);
               TH1F *h = (TH1F *) f->Get(hn);
               if (h) {
                  hs_found = kTRUE;
                  if (!fHist[i]) {
                     fHist[i] = (TH1F *) h->Clone();
                     fHist[i]->SetDirectory(0);
                  } else {
                     fHist[i]->Add(h);
                  }
               } else {
                  Error("GetHistosFromFC", "histo '%s' not found in file '%s'",
                        hn.Data(), fi->GetCurrentUrl()->GetUrl());
               }
            }
            f->Close();
         } else {
            Error("GetHistosFromFC", "file '%s' could not be open", fi->GetCurrentUrl()->GetUrl());
         }
      }
      if (hs_found) break;
   }
   if (!fc_found) return -1;
   if (!hs_found) return -1;
   
   for (Int_t i = 0; i < fNhist; i++) {
      cv->cd(i+1);
      if (fHist[i]) {
         fHist[i]->DrawCopy();
      }
   }
   Info("GetHistosFromFC", "histograms read from %d files in TFileCollection '%s'",
                           fc->GetList()->GetSize(), fc->GetName());
   // Done
   return 0;
}

void peaks(Int_t np=10) {
   npeaks = TMath::Abs(np);
   TH1F *h = new TH1F("h","test",500,0,1000);
   //generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
   Int_t p;
   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 1;
      par[3*p+3] = 10+gRandom->Rndm()*980;
      par[3*p+4] = 3+2*gRandom->Rndm();
   }
   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
   f->SetNpx(1000);
   f->SetParameters(par);
   TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
   c1->Divide(1,2);
   c1->cd(1);
   h->FillRandom("f",200000);
   h->Draw();
   TH1F *h2 = (TH1F*)h->Clone("h2");
   //Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,2,"",0.10);
   printf("Found %d candidate peaks to fit\n",nfound);
   //Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,20,"same");
   if (hb) c1->Update();
   if (np <0) return;

   //estimate linear background using a fitting method
   c1->cd(2);
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h->Fit("fline","qn");
   //Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Double_t *xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
      par[3*npeaks+4] = 3;
      npeaks++;
   }
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   //we may have more than the default 25 parameters
   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h2->Fit("fit");
}

//--------------------------------------------------------------------------------------------------
// make fake rate histogram from numerator and denominators
//--------------------------------------------------------------------------------------------------
// TGraphAsymmErrors* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames, 
//                            int denominatortype, int faketype, int chargetype,  
//                            int eventType, string histName) {
TH1F* createFakeRatePtHist(vector<string> datasetFiles, vector<string> datasetNames, 
                           int denominatortype, int faketype, int chargetype,  
                           int eventType, string histName) {

  string numeratorHistName = "";
  string denominatorHistName = "";
  string dirName = "";

  string chargeTypeName = "";
  if (chargetype == -1)
    chargeTypeName = "Minus";
  else if (chargetype == 1)
    chargeTypeName = "Plus";

  string eventTypeName = "";
  if (eventType == -1)
    eventTypeName = "MinusW";
  else if (eventType == 1)
    eventTypeName = "PlusW";

  string denominatorTypeName = "";
  if (denominatortype == 1)
    denominatorTypeName = "Track";

  if (faketype == 11) {
    dirName = "ComputeElectronFakeRateMod";
    numeratorHistName = "h" + chargeTypeName + "Electron" + eventTypeName + "NumeratorPt";
    denominatorHistName = "h" + chargeTypeName + "Electron" + eventTypeName + denominatorTypeName + "DenominatorPt";
  } else if (faketype == 13) {
    dirName = "ComputeMuonFakeRateMod";
    numeratorHistName = "h" + chargeTypeName + "Muon" + eventTypeName + "NumeratorPt";
    denominatorHistName = "h" + chargeTypeName + "Muon" + eventTypeName + denominatorTypeName + "DenominatorPt";
  } else {
    cerr << "Error in createFakeRatePtHist: faketype = " << faketype << " is not recognized." << endl;
    assert(false);
  }

  TH1F *denominator = addAllSamples(datasetFiles, datasetNames, dirName, denominatorHistName);
  TH1F *numerator = addAllSamples(datasetFiles, datasetNames, dirName, numeratorHistName);

  //create fake rate hist from denominator binning
  TH1F *fakeRateHist = (TH1F*)numerator->Clone(histName.c_str());
  fakeRateHist->GetYaxis()->SetTitle("Fake Rate");
  //divide by the denominator to get fake rate;
  fakeRateHist->Divide(denominator);

  //TGraphAsymmErrors* fakeRate = (numerator, denominator);
  

  //for (int i=1 ; i <= fakeRateHist->GetXaxis()->GetNbins(); i++) {
    //cout << i << " : " << fakeRateHist->GetXaxis()->GetBinCenter(i) << " " << fakeRateHist->GetBinContent(i) << "  +- " << fakeRateHist->GetBinError(i) << endl;
    //if (fakeRateHist->GetBinContent(i) > 0 && fakeRateHist->GetBinError(i) < 0.01) {
      //fakeRateHist->SetBinError(i,0.01);
    //}
  //}

  return fakeRateHist;
}

void JetTagBin::inputData( TFile * file ) 
{
  TH1F * tempDataHist = (TH1F*)file->Get(("Data_"+jtBinName_).c_str());
  if ( verbose ) cout << "Loaded Data_" << jtBinName_ << " from the file." << endl;
  jtDataHist_ = (TH1F*) tempDataHist->Clone();
  //include check on binning with respect to other histograms
  if ( jtData_ ) delete jtData_;
  jtData_ = new RooDataHist(("Data_"+jtBinName_).c_str(),("Data_"+jtBinName_).c_str(),*jtBinVar_,Import(*jtDataHist_));
}

void displaySingleChannelWaveforms(TString fdat, int board, int channel) {

  gStyle->SetOptStat(0);

  TFile *f = new TFile(fdat);
  if (f->IsZombie()){
    cout << "Cannot open file: " << fdat << endl;
    return;
  }

  TBEvent *event = new TBEvent();
  TTree *t1041 = (TTree*)f->Get("t1041"); 
  TBranch *bevent = t1041->GetBranch("tbevent");
  bevent->SetAddress(&event);

  TCanvas * canv = new TCanvas("canv", "canv", 2000, 2000);
  canv->cd();

  TH1F * dummy = new TH1F("dummy", "dummy", 120, 0, 120);
  dummy->GetYaxis()->SetRangeUser(0, 2500);
  dummy->Draw();
        
  vector<TH1F*> waves;

  TH1F * wave = new TH1F("wave", "wave", 120, 0, 120);

  int nplots = 0;

  for (Int_t i = 0; i < t1041->GetEntries(); i++) {
    t1041->GetEntry(i);
    
    for (int j = 0; j < event->NPadeChan(); j++){
      PadeChannel pch = event->GetPadeChan(j);
      
      if((int)pch.GetBoardID() != board || (int)pch.GetChannelID() != channel) continue;
      
      pch.GetHist(wave);
      
      nplots++;
      TH1F * wavecopy = (TH1F*)wave->Clone("wave_"+TString(Form("%d", nplots)));
      waves.push_back(wavecopy);
      
    }
    
  }
  
  int nBigPeaks = 0;
  for(unsigned int ui = 0; ui < waves.size(); ui++) {
    if(waves[ui]->GetMaximum() > 400) {
      waves[ui]->SetLineColor(nBigPeaks+2);
      nBigPeaks++;
    }
    waves[ui]->Draw("same");
  }
  
}

//----------  Retrieve data histo  -----------------
TH1F* GetData(TFile* fin, string region, string varname)
{
    string cname = CHANNEL_NAME+string("/")+region+"/"+varname;
    TCanvas* c = (TCanvas*) fin->Get(cname.c_str());
    TList* l = c->GetListOfPrimitives();
    TPad* pad = (TPad*) l->At(0);
    string hname = "v:"+varname+"|r:"+region+string("|c:")+CHANNEL_NAME+string("|t:1DSumData");
    TH1F* h = (TH1F*) pad->GetPrimitive(hname.c_str());
    return (TH1F*) h->Clone();
}