C++ TH1D::GetBinContent方法代码示例

TF1* GetFitFunc_ZYAM(TH1D* h)
{
                    TH1D* hcorrphi = (TH1D*)h->Clone(h->GetName());
                    double histminY = hcorrphi->GetBinContent(10);
                    double histminX = 1.0;

                    //hcorrphi->SetAxisRange(-0.01,1.2,"X"); 
                    TF1* fitfunc = new TF1(Form("fitfunc_%s",h->GetName()),"[0]+[1]*(x-[2])*(x-[2])",1.0,2.4); //std 0.6 1.55 vs pT ; 0.6 1.8 vs eta
                    fitfunc->SetParameters(histminY,0.0002,histminX);
                    fitfunc->SetParLimits(1,0,1000);
                    fitfunc->SetParLimits(2,1.0,2.4);
                   
/*
//                    hcorrphi->SetAxisRange(-0.01,1.2,"X"); 
                    TF1* fitfunc = new TF1(Form("fitfunc_%s",h->GetName()),"[0]+[1]*(x-[2])*(x-[2])",0.65,1.3); 
                    fitfunc->SetParameters(histminY,0.0002,histminX);
                    fitfunc->SetParLimits(1,0,1000);
                    fitfunc->SetParLimits(2,0.05,1000);
*/
                    for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%s",h->GetName()),"RNO");
                    return fitfunc;
}

Spectrum roothist_to_spectrum(const TH1D & hist, const TH2D * cov, bool use_hist_uncertainties){
    const int n = hist.GetNbinsX();
    if(cov){
        if(cov->GetDimension() != 2 || cov->GetNbinsX() != n || cov->GetNbinsY() != n){
            throw runtime_error(string("covariance histogram '") + cov->GetName() + "' has wrong dimension");
        }
    }
    Spectrum s(n);
    for(int i=0; i<n; ++i){
        double c = hist.GetBinContent(i+1);
        if(!isfinite(c)){
            throw runtime_error(string("non-finite entry in histogram '") + hist.GetName() + "'");
        }
        s[i] = c;
        if(use_hist_uncertainties){
            double e = hist.GetBinError(i+1);
            if(!isfinite(e)){
                throw runtime_error(string("non-finite error in histogram '") + hist.GetName() + "'");
            }
            s.cov()(i,i) += e*e;
        }
    }
    if(cov){
        for(int i=0; i<n; ++i){
            for(int j=0; j<n; ++j){
                double c_ij = cov->GetBinContent(i+1, j+1);
                double c_ji = cov->GetBinContent(j+1, i+1);
                if(!isfinite(c_ij)){
                    throw runtime_error(string("covariance histogram '") + cov->GetName() + "' does have non-finite entry");
                }
                if(fabs(c_ij - c_ji) > 1e-8 * max(fabs(c_ij), 1.0)){
                    cerr << "covariance histogram '" << cov->GetName() <<  "' is not symmetric: C(" << i << ","<< j << ") = " << c_ij << "; transposed: " << c_ji << "; diff = " << (c_ij - c_ji) << endl;
                }
                s.cov()(i,j) += c_ij;
            }
        }
    }
    return s;
}

double CalcChiSqr(TH1D* hist, Double_t coeff, Double_t expo, Int_t xMin, Int_t xMax)
{
  {
    TH1D* myHist = hist->Clone();
    int bins = myHist->GetSize()-2;
    double chiSqr = 0;
    for(int i = xMin+1; i < xMax+1; i++)
      {
	double x = myHist->GetBinCenter(i);
	double o = myHist->GetBinContent(i);
	double e = coeff*(TMath::Power(x,expo));
	//cout << x << "\t" << o << "\t" << e;
	if(o != 0)
	  chiSqr += ((o-e)*(o-e))/o;
	//cout << "\tchiSqr: " << chiSqr << endl;
      }
    return chiSqr/2;
    
  }
  
  
}

void fitpapvar(const char* infilename, const char* histname) {
  TFile* infile = new TFile(Form("%s",infilename),"read");
  TH1D* hcf = (TH1D*)infile->Get(histname);
  hcf->GetXaxis()->SetRangeUser(0,2);

   TF1 *fc2 = new TF1("fc2","[3]+[2]*TMath::Gaus(x,[0],[1])",0.4,1);
   fc2->SetParameters(1.5,0.3,-0.2,1);

  // TF1 *fc2 = new TF1("fc2","[1]+[0]*x*x",0,1);
  // fc2->SetParameters(-0.01,1.0);

  // TF1 *fc2 = new TF1("fc2","[2]+[1]*x*x*x*x*x",0,1);
  // fc2->SetParameters(-0.01,0.01,1.0);

  // TF1 *fc2 = new TF1("fc2","[3]+[2]*x*x+[1]*x*x*x+[0]*x*x*x*x",0.,1);
  // fc2->SetParameters(-0.01,0.01,1.0,1.0);

  hcf->Fit("fc2","r","",0.3,1);
  TH1D* hnew = new TH1D("hnew","hnew",hcf->GetNbinsX(),0,1);
  for(int i=1;i<=hcf->GetNbinsX();++i){
    // cout << i << endl;
    // cout << hcf->GetBinContent(i)/fc2->Eval(2.*i/hcf->GetNbinsX()) << endl;
    // cout << hcf->GetBinContent(i) << endl;
    // cout << fc2->Eval(2.*i/hcf->GetNbinsX()) << endl << endl;
    hnew->SetBinContent(i, hcf->GetBinContent(i)/fc2->Eval(1.*i/hcf->GetNbinsX()));
    hnew->SetBinError(i, hcf->GetBinError(i));

  }
  hnew->Draw("same");
  fc2->Draw("same");
  hnew->SetName(Form("divp4%s",histname));
  TFile* ofile = new TFile(Form("divp4%s",infilename),"update");
  hnew->Write();
  hcf->Write();
  ofile->Close();
}

TF1* fitDstar(TTree* nt, TTree* ntMC, Float_t ptmin, Bool_t plotgenmatch)
{
  TCanvas* c = new TCanvas(Form("c_5p_%.0f",ptmin),"",600,600);
  TH1D* h = new TH1D(Form("h_5p_%.0f",ptmin),"",BINNUM,BINMIN,BINMAX);
  TH1D* hMCSignal = new TH1D(Form("hMCSignal_5p_%.0f",ptmin),"",BINNUM,BINMIN,BINMAX);
  TH1D* hMCSignalplot = new TH1D(Form("hMCSignalplot_5p_%.0f",ptmin),"",BINNUM,BINMIN,BINMAX);
  TH1D* hMCSwapped = new TH1D(Form("hMCSwapped_5p_%.0f",ptmin),"",BINNUM,BINMIN,BINMAX);
  TH1D* hMCSwappedplot = new TH1D(Form("hMCSwappedplot_5p_%.0f",ptmin),"",BINNUM,BINMIN,BINMAX);

  TF1* f = new TF1(Form("f_5p_%.0f",ptmin),"[0]*([4]*([6]*([12]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[12])*Gaus(x,[1],[11])/(sqrt(2*3.14159)*[11]))+(1-[6])*Gaus(x,[1],[5])/(sqrt(2*3.14159)*[5]))+(1-[4])*Gaus(x,[1],[3])/(sqrt(2*3.14159)*[3]))+[10]*((1-exp((0.13957-x)/[7]))*pow(x/0.13957,[8])+[9]*(x/0.13957-1))",BINMIN,BINMAX);
  nt->Project(Form("h_5p_%.0f",ptmin),"Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightdata[isData].Data(),seldata5p[isData].Data(),triggerselectiondata[isData].Data(),ptmin));
  ntMC->Project(Form("hMCSignal_5p_%.0f",ptmin),"Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightmc[isData].Data(),selmc5p[isData].Data(),triggerselectionmc[isData].Data(),ptmin));
  ntMC->Project(Form("hMCSwapped_5p_%.0f",ptmin),"Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightmc[isData].Data(),selswp5p[isData].Data(),triggerselectionmc[isData].Data(),ptmin));
  for(int ibin=0;ibin<BINNUM;ibin++) hMCSignalplot->SetBinContent(ibin+1,hMCSignal->GetBinContent(ibin+1));
  for(int ibin=0;ibin<BINNUM;ibin++) hMCSwappedplot->SetBinContent(ibin+1,hMCSwapped->GetBinContent(ibin+1));
  f->FixParameter(4,1.);
  f->FixParameter(1,0.145491);
  f->FixParameter(10,0);

  f->SetParLimits(0,0,1.e+5);
  f->SetParLimits(6,0,1.);
  f->SetParLimits(12,0,1.);
  f->SetParLimits(2,3.e-4,1.e-3);
  f->SetParameter(2,5.e-4);
  f->SetParLimits(11,1.6e-4,3.e-4);//1.5e-4 keyong
  f->SetParameter(11,2.e-4);
  f->SetParLimits(5,1.e-3,1.6e-3);
  f->SetParameter(5,1.e-3);
  hMCSignal->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);
  hMCSignal->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);
  f->ReleaseParameter(1);
  f->SetParLimits(1,minmass,maxmass);
  hMCSignal->Fit(Form("f_5p_%.0f",ptmin),"LL","",minmass,maxmass);
  hMCSignal->Fit(Form("f_5p_%.0f",ptmin),"LL","",minmass,maxmass);

  f->FixParameter(1,f->GetParameter(1));
  f->FixParameter(2,f->GetParameter(2));
  f->FixParameter(5,f->GetParameter(5));
  f->FixParameter(11,f->GetParameter(11));
  f->FixParameter(6,f->GetParameter(6));
  f->FixParameter(12,f->GetParameter(12));
  f->FixParameter(4,0);
  f->SetParLimits(3,2.e-4,2.e-3);
  f->SetParameter(3,1.e-3);

  hMCSwapped->Fit(Form("f_5p_%.0f",ptmin),"L q","",BINMIN,BINMAX);
  hMCSwapped->Fit(Form("f_5p_%.0f",ptmin),"L q","",BINMIN,BINMAX);
  hMCSwapped->Fit(Form("f_5p_%.0f",ptmin),"L q","",minmass,maxmass);
  hMCSwapped->Fit(Form("f_5p_%.0f",ptmin),"L q","",minmass,maxmass);

  f->FixParameter(4,hMCSignal->Integral(0,1000)/(hMCSwapped->Integral(0,1000)+hMCSignal->Integral(0,1000)));
  f->FixParameter(3,f->GetParameter(3));

  f->SetParLimits(7,5.e-4,1.e-2);
  f->SetParameter(7,1.6e-3);
  f->SetParLimits(8,0.,15.);
  f->SetParameter(8,0.35);
  f->SetParLimits(9,-2.e+1,2.e+1);
  f->SetParameter(9,13.);

  f->ReleaseParameter(10);
  f->SetParLimits(10,0,1.e+6);

  h->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);
  h->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);

  f->ReleaseParameter(1);
  f->SetParLimits(1,minmass,maxmass);
  f->SetParameter(1,f->GetParameter(1));
  h->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);
  h->Fit(Form("f_5p_%.0f",ptmin),"LL","",BINMIN,BINMAX);

  TF1* background = new TF1(Form("background_5p_%.0f",ptmin),"[3]*((1-exp((0.13957-x)/[0]))*pow(x/0.13957,[1])+[2]*(x/0.13957-1))");
  background->SetParameters(f->GetParameter(7),f->GetParameter(8),f->GetParameter(9),f->GetParameter(10));
  background->SetRange(BINMIN,BINMAX);
  background->SetLineColor(4);
  background->SetLineWidth(3);
  background->SetLineStyle(2);

  TF1* mass = new TF1(Form("fmass_5p_%.0f",ptmin),"[0]*[3]*([5]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[6])/(sqrt(2*3.14159)*[6]))+(1-[5])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
  mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(4),f->GetParameter(5),f->GetParameter(6),f->GetParameter(11),f->GetParameter(12));
  mass->SetParError(0,f->GetParError(0));
  mass->SetParError(1,f->GetParError(1));
  mass->SetParError(2,f->GetParError(2));
  mass->SetParError(3,f->GetParError(4));
  mass->SetParError(4,f->GetParError(5));
  mass->SetParError(5,f->GetParError(6));
  mass->SetRange(BINMIN,BINMAX);
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);

  TF1* massSwap = new TF1(Form("fmassSwap_5p_%.0f",ptmin),"[0]*(1-[2])*Gaus(x,[1],[3])/(sqrt(2*3.14159)*[3])");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(4),f->GetParameter(3));
  massSwap->SetRange(BINMIN,BINMAX);
  massSwap->SetFillColor(kGreen+4);
  massSwap->SetFillStyle(3005);
  massSwap->SetLineColor(kGreen+4);
//.........这里部分代码省略.........

void fitDstar5pMinpt(Bool_t genmatchpoint=true)
{
  gStyle->SetTextSize(0.05);
  gStyle->SetTextFont(42);
  gStyle->SetPadRightMargin(0.043);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadBottomMargin(0.145);
  gStyle->SetTitleX(.0f);

  void clean0(TH1D* h);  
  TF1* fitDstar(TTree* nt, TTree* ntMC, Float_t ptmin, Bool_t plotgenmatch);

  TFile* infData = new TFile(infnameData5p[isData].Data());
  TFile* infMC = new TFile(infnameMC5p[isData].Data());

  TTree* ntData = (TTree*)infData->Get("ntDD0kpipipipi");
  TTree* ntMC = (TTree*)infMC->Get("ntDD0kpipipipi");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  
  ntData->AddFriend("ntHlt");
  if(isData!=Data_MB&&isData!=Data) ntData->AddFriend("ntHi");
  ntMC->AddFriend("ntHlt");
  ntMC->AddFriend("ntHi");
  ntGen->AddFriend("ntHlt");
  ntGen->AddFriend("ntHi");

  Float_t aZero[nBins];
  for(int i=0;i<nBins;i++) aZero[i]=0;
  Float_t aPt[nBins],aPtErr[nBins],aGen[nBins],aGenErr[nBins];
  TH1F* hPt = new TH1F("hPt","",nBins,ptBinsPlus);
  TH1F* hGen = new TH1F("hGen","",nBins,ptBinsPlus);

  for(int i=0;i<nBins;i++)
    {
      TF1* fData = fitDstar(ntData,ntMC,ptBins[i],genmatchpoint);
      Float_t yieldData = fData->Integral(BINMIN,BINMAX)/BINWID;
      Float_t yieldDataErr = fData->Integral(BINMIN,BINMAX)/BINWID*fData->GetParError(0)/fData->GetParameter(0);      
      aPt[i] = yieldData;
      aPtErr[i] = yieldDataErr;
      hPt->SetBinContent(i+1,aPt[i]);
      hPt->SetBinError(i+1,aPtErr[i]);
      TH1D* hGenFill = new TH1D(Form("hGenFill_%.0f",ptBins[i]),"",1,ptBins[i],1.e+3);
      hGenFill->Sumw2();
      ntGen->Project(Form("hGenFill_%.0f",ptBins[i]),"Gpt",TCut("%s",weightmc[isData].Data())*Form("%s&&%s",selgen5p.Data(),triggerselectionmc[isData].Data()));
      aGen[i] = hGenFill->GetBinContent(1);
      aGenErr[i] = hGenFill->GetBinError(1);
      hGen->SetBinContent(i+1,aGen[i]);
      hGen->SetBinError(i+1,aGenErr[i]);
    }

  TGraphErrors* gPt = new TGraphErrors(nBins,ptBins,aPt,aZero,aPtErr);
  gPt->SetName("gPt");
  TGraphErrors* gGen = new TGraphErrors(nBins,ptBins,aGen,aZero,aGenErr);
  gGen->SetName("gGen");

  TFile* outputfile = new TFile(Form("outputfiles/output_5p_%s_Minpt.root",texData[isData].Data()),"recreate");
  outputfile->cd();
  gPt->Write();
  gGen->Write();
  hPt->Write();
  hGen->Write();
  outputfile->Close();
}

//.........这里部分代码省略.........
      ppos[i] = new TH1D(Form("ppos%d",i),"",nBinAj,ajBins);
      ppos[i]->SetLineColor(1);     
      ppos[i]->SetMarkerColor(colors[i]);
      ppos[i]->SetFillColor(colors[i]);
      ppos[i]->SetFillStyle(1001);
      pneg[i] = new TH1D(Form("pneg%d",i),"",nBinAj,ajBins);
      pneg[i]->SetLineColor(1);     
      pneg[i]->SetMarkerColor(colors[i]);
      pneg[i]->SetFillColor(colors[i]);
      pneg[i]->SetFillStyle(1001);

      // =================================
      // Caculated Stat Error of the Mean
      // =================================
      cout << "Stat Error for pt bin " << i << ": ";
      for ( int iaj = 0 ; iaj< nBinAj ; iaj++) {
         he[iaj] = new TH1D(Form("he%d_aj%d",i,iaj),"",100,-200,200);
         TCut ajCut = Form("Aj>%f && Aj<%f", ajBins[iaj],ajBins[iaj+1]);
         t->Draw(Form("((metxMerged%d))>>he%d_aj%d",i,i,iaj), "weight" * evtCut&&myCut&&ajCut);
         float theError = he[iaj]->GetRMS()/ (sqrt(he[iaj]->GetEntries()));
	 cout << theError << " ";
	 pe[i]->SetBinError(iaj+1, theError);
      }
      cout << endl;
   }

   // Stack
   for (int i=nBin-1;i>=0;i--)
   {
      for(int iaj = 0 ; iaj< nBinAj ; iaj++) {
         double posVal=0, negVal=0;
         double posValErr=0, negValErr=0;
         if (i!=nBin-1) {
            posVal = ppos[i+1]->GetBinContent(iaj+1);
            posValErr = ppos[i+1]->GetBinError(iaj+1);
            negVal = pneg[i+1]->GetBinContent(iaj+1);
            negValErr = pneg[i+1]->GetBinError(iaj+1);
         }
         if (pe[i]->GetBinContent(iaj+1)<0) {
            negVal+=pe[i]->GetBinContent(iaj+1);
            negValErr=pe[i]->GetBinError(iaj+1);
            posValErr=0;
         } else if (pe[i]->GetBinContent(iaj+1)>0) {
            posVal+=pe[i]->GetBinContent(iaj+1);
            posValErr=pe[i]->GetBinError(iaj+1);
            negValErr=0;
         } 
         ppos[i]->SetBinContent(iaj+1,posVal);
         ppos[i]->SetBinError(iaj+1,posValErr);
         pneg[i]->SetBinContent(iaj+1,negVal);
         pneg[i]->SetBinError(iaj+1,negValErr);
      }
   }
   
   TH1D *pall;
   TH1D *pallE;
   TH1D *h1 = new TH1D(Form("hAll1"),"",nBinAj,ajBins);
   TH1D *h2 = new TH1D(Form("hAll2"),"",nBinAj,ajBins);

   h1->Sumw2();
   h2->Sumw2();
   t->Draw(Form("Aj>>hAll1"), "weight"*(evtCut&&myCut));
   t->Draw(Form("Aj>>hAll2"), Form("((-weight*metxMergedAll))")*(evtCut&&myCut));
   pall=(TH1D*)h2->Clone();
   pall->SetName("pall");
   pall->Divide(h1);

void fit(const char *run="428211_429133_5s",
	 int key=1, int bmin=10, bool draw=true, bool pa=false) {
  int minentries=1000;
  gSystem->Exec( Form("mkdir -p %s/SEN%03d",run,key/128) );
  gSystem->Exec( Form("mkdir -p %s/SEN%03d",run,key/128) );
  int state = findstate(key);
  printf("state %d\n",state);
  // data
  TString inname = Form("%s/adc/HI_KEY%05d.root",run,key);
  TString outname = Form("HI_KEY%05d",key);
  TFile *file = new TFile( inname.Data() );
  cout << inname.Data() << endl;

  TH1D *out = (TH1D*) file->Get("out");
  double xfit_min=out->GetBinLowEdge(bmin);
  double xfit_max=122.5;
  int bmax = out->GetXaxis()->FindBin(xfit_max);
  int entries = out->Integral(bmin,bmax);
  if(entries<minentries) {
    cout << "not enough entries: ";
    cout << entries << endl;
    return;
  }

  // fit
  TCanvas *main = new TCanvas("main","main");
  int pkt=0;
  if(pa)
    pkt = (key%(8*4*12*64))/(4*12*64);
  TF1 *fitH = GetFit( Form("%s/SEN%03d/%s.dat",run,key/128,outname.Data()) ,pkt,xfit_min);
  out->Fit(fitH,"MELIR","",xfit_min,xfit_max);
  TF1 *MIPH1 = GetMIP(fitH,1,kCyan-3);
  TF1 *MIPH2 = GetMIP(fitH,2,kGreen-3);
  TF1 *MIPH3 = GetMIP(fitH,3,kOrange-3);
  TF1 *MIPH4 = GetMIP(fitH,4,kMagenta-3);
  TF1 *BGR = GetBGR(fitH,xfit_min);

  double amp = fitH->GetParameter(0);
  double eamp= fitH->GetParError(0);
  double lda = fitH->GetParameter(1);
  double elda= fitH->GetParError(1);
  double sg1 = fitH->GetParameter(2);
  double esg1= fitH->GetParError(2);
  double fr2 = fitH->GetParameter(3);
  double efr2= fitH->GetParError(3);
  double fr3 = fitH->GetParameter(4);
  double efr3= fitH->GetParError(4);
  double fr4 = fitH->GetParameter(5);
  double efr4= fitH->GetParError(5);
  double fr1 = 1 - fr2 - fr3 - fr4;
  double ncs = fitH->GetChisquare()/fitH->GetNDF();
  double ba = fitH->GetParameter(6);
  double eba= fitH->GetParError(6);
  double bsl = fitH->GetParameter(7);
  double ebsl= fitH->GetParError(7);

  // saving fit
  ofstream outfit;
  outfit.open( Form("%s/SEN%03d/%s.dat",run,key/128,outname.Data()) );
  outfit << amp << " " << eamp << endl;
  outfit << lda << " " << elda << endl;
  outfit << sg1 << " " << esg1 << endl;
  outfit << fr2 << " " << efr2 << endl;
  outfit << fr3 << " " << efr3 << endl;
  outfit << fr4 << " " << efr4 << endl;
  outfit << ba  << " " << eba  << endl;
  outfit << bsl << " " << ebsl << endl;
  outfit << ncs << endl;
  outfit.close();
  cout << "Parameters saved to ";
  cout << outname.Data() << ".dat" << endl;

  // draw
  if(!draw) return;

  gStyle->SetOptFit(0);
  gStyle->SetOptStat(0);
  out->Draw("HE");
  double ymax = out->GetBinContent( out->FindBin(xfit_min) )*1.5;
  out->GetYaxis()->SetRangeUser(0.5,ymax);
  out->GetXaxis()->SetRangeUser(-5,125);
  out->Sumw2();
  out->SetLineColor(kBlack);
  out->SetMarkerStyle(20);
  out->SetTitle("");
  out->GetXaxis()->SetTitle("ADC-PED (a.u.)");
  BGR->Draw("SAME");
  MIPH1->Draw("SAME");
  MIPH2->Draw("SAME");
  MIPH3->Draw("SAME");
  MIPH4->Draw("SAME");
  fitH->SetRange(xfit_min,xfit_max);
  fitH->Draw("SAME");
  TLatex *text = new TLatex();
  text->DrawLatex(0, (1.03*(ymax)), inname.Data() );
  text->DrawLatex(30, (0.83*(ymax)), Form("Entries  %d",entries) );
  text->DrawLatex(30, (0.73*(ymax)), Form("State  %d",state) );
  text->DrawLatex(30, (0.53*(ymax)), Form("#lambda  %.1f #pm %.1f",lda,elda) );
  text->DrawLatex(30, (0.43*(ymax)), Form("#sigma  %.1f #pm %.1f",sg1,esg1) );
  text->SetTextColor(kRed-3);
//.........这里部分代码省略.........

int purityRandNorm(TH1D* h_template, TString name , TFile * fData, TFile * fZinv, TFile * fDY, int & lastmt2val_hybrid) {

  //cout<<"purityRandNorm for template "<<name<<endl;
  //h_template->Print("all");
  int lastbin_hybrid = 0;
  lastmt2val_hybrid = 200;
  TString name_emu = name + "emu";
  TString name_zinv = name;
  name_zinv.ReplaceAll("crdy", "sr");
  TH1D* hEMU   = (TH1D*) fData->Get(name_emu);    
  TH1D* hDY    = (TH1D*) fDY->Get(name);    
  TH1D* hZinv  = (TH1D*) fZinv->Get(name_zinv); 
  if (h_template == 0) {
    cout<<"ZinvMaker::purityAndRatio  :  could not find input template"<<endl;
    return lastbin_hybrid;
  }
  if (hDY == 0 || hZinv == 0) {
    cout<<"ZinvMaker::purityAndRatio  :  could not find DY or Zinv MC histogram"<<endl;
    return lastbin_hybrid;
  }
  if (hEMU) h_template->Add(hEMU, -1*rSFOF);

  // find the last bin
  for ( int ibin=1; ibin <= hDY->GetNbinsX(); ++ibin ) {
    float integratedYield = 0;
    integratedYield = hDY->Integral(ibin,-1);
    if (integratedYield < hybrid_nevent_threshold) {
      if (ibin == 1) lastbin_hybrid = 1;
      else {
	lastbin_hybrid = ibin-1;
	lastmt2val_hybrid = hDY->GetBinLowEdge(ibin);
      }
      break;
    }
  }

  // multiply R(Znn/Zll)
  TH1D* ratio =   (TH1D*) hZinv->Clone("ratio");
  ratio->Divide(hDY);

  h_template->Multiply(ratio);

  // Get the integrals to normalize the Zinv tails
  // and the uncertainties on the CR yield (dominated by data stats in the last N bins)
  double integratedYieldErrZinv = 0;
  float integratedYieldZinv = hZinv->IntegralAndError(lastbin_hybrid, -1., integratedYieldErrZinv);
  float relativeErrorZinv = integratedYieldErrZinv/integratedYieldZinv;
  double integratedYieldErr = 0;
  float integratedYield = h_template->IntegralAndError(lastbin_hybrid,-1,integratedYieldErr); 
  float relativeError = integratedYieldErr/integratedYield;

  // Hybridize the template: last N bins have a common stat uncertainty, and they follow the Zinv MC shape
  for ( int ibin=1; ibin <= hZinv->GetNbinsX(); ++ibin ) {

    if (ibin < lastbin_hybrid) continue;

    float kMT2 = hZinv->GetBinContent(ibin) / integratedYieldZinv;
    float err = sqrt( relativeError*relativeError + relativeErrorZinv*relativeErrorZinv);
    h_template->SetBinContent(ibin, integratedYield * kMT2);
    h_template->SetBinError(ibin, integratedYield * kMT2 * err );

  }

  // Normalize it: we just need a shape after all
  h_template->Scale(1./h_template->Integral());

  //h_template->Print("all");


  return lastbin_hybrid;
}

void MiniTreeSignalProducerVV13TeV(int samplemin, int samplemax){

 string dir = "/usr/users/dschaefer/CMSSW_7_4_7/src/DijetCombineLimitCode/";
  double mgg, mjj,evWeight, mtot, normWeight;
 int categories;

 evWeight = 1.0;
 normWeight = 1;

 for (int iSample = samplemin; iSample < samplemax; iSample++){
   
   string inFile("WprimeToWZ");
   if (iSample == 1) inFile = string("RS1WW");
   if (iSample == 2) inFile = string("RS1ZZ"); // Fake ZZ signal
   if (iSample == 3) inFile = string("QstarQW");
   if (iSample == 4) inFile = string("QstarQZ");
   if (iSample == 5) inFile = string("BulkWW");
   if (iSample == 6) inFile = string("BulkZZ");
   if (iSample == 7) inFile = string("ZprimeWW");
   if (iSample == 8) inFile = string("WprimeWZ");

   string outFile("dijetVV_13TeV_WZ");
   if (iSample == 1) outFile = string("dijetVV_13TeV_RS1WW");
   if (iSample == 2) outFile = string("dijetVV_13TeV_RS1ZZ");
   if (iSample == 3) outFile = string("dijetVV_13TeV_QstarQW");
   if (iSample == 4) outFile = string("dijetVV_13TeV_QstarQZ");
   if (iSample == 5) outFile = string("dijetVV_13TeV_BulkWW");
   if (iSample == 6) outFile = string("dijetVV_13TeV_BulkZZ");
   if (iSample == 7) outFile = string("dijetVV_13TeV_ZprimeWW");
   if (iSample == 8) outFile = string("dijetVV_13TeV_WZ");
   
   int massrange=36;

   for (int iMass = 0; iMass<massrange; iMass++){

     string sInFile = dir+"input/" + inFile + "_13TeV_" + Form("10k_OUT%dGeV.root", 1000+iMass*100);
     cout << sInFile.c_str() << endl;
     TFile file0(sInFile.c_str(), "read");

     string sOutFile = dir+"MiniTrees/Signal_VV_13TeV/" + outFile + Form("OUT%d_miniTree.root", 1000+iMass*100);
     TFile f1(sOutFile.c_str(), "recreate");
     f1.cd();

     TTree *TCVARS = new TTree("TCVARS", "VV selection");
     TCVARS->Branch("mgg13TeV",&mgg,"mgg/D");

     TCVARS->Branch("evWeight",&evWeight,"evWeight/D");
     TCVARS->Branch("normWeight",&normWeight,"normWeight/D");
     
     TCVARS->Branch("categories",&categories,"categories/I");

  
     double dMass = 1000.+iMass*100.;

     
     for (int iCat = 2; iCat < 8; iCat++){
       TH1D* hMass = (TH1D*) file0.Get("DijetMassHighPuriVV;1");
       if (iCat == 1) hMass = (TH1D*) file0.Get("DijetMassLowPuriVV;1");
       
       if (iCat == 2) hMass = (TH1D*) file0.Get("DijetMassHighPuriWW;1");
       if (iCat == 3) hMass = (TH1D*) file0.Get("DijetMassLowPuriWW;1");
       
       if (iCat == 4) hMass = (TH1D*) file0.Get("DijetMassHighPuriWZ;1");
       if (iCat == 5) hMass = (TH1D*) file0.Get("DijetMassLowPuriWZ;1");
       
       if (iCat == 6) hMass = (TH1D*) file0.Get("DijetMassHighPuriZZ;1");
       if (iCat == 7) hMass = (TH1D*) file0.Get("DijetMassLowPuriZZ;1");
    
       if (iCat == 8) hMass = (TH1D*) file0.Get("DijetMassHighPuriqV;1");
       if (iCat == 9) hMass = (TH1D*) file0.Get("DijetMassLowPuriqV;1");

       if (iCat == 10) hMass = (TH1D*) file0.Get("DijetMassHighPuriqW;1");
       if (iCat == 11) hMass = (TH1D*) file0.Get("DijetMassLowPuriqW;1");

       if (iCat == 12) hMass = (TH1D*) file0.Get("DijetMassHighPuriqZ;1");
       if (iCat == 13) hMass = (TH1D*) file0.Get("DijetMassLowPuriqZ;1");

       if (iCat == 14) hMass = (TH1D*) file0.Get("DijetMassNoPuriVV;1");
       if (iCat == 15) hMass = (TH1D*) file0.Get("DijetMassNoPuriWW;1");
       if (iCat == 16) hMass = (TH1D*) file0.Get("DijetMassNoPuriWZ;1");
       if (iCat == 17) hMass = (TH1D*) file0.Get("DijetMassNoPuriZZ;1");
       if (iCat == 18) hMass = (TH1D*) file0.Get("DijetMassNoPuriqV;1");
       if (iCat == 19) hMass = (TH1D*) file0.Get("DijetMassNoPuriqW;1");
       if (iCat == 20) hMass = (TH1D*) file0.Get("DijetMassNoPuriqZ;1");
       // TH1D* hMass = (TH1D*) file0.Get("DijetMassHighPuriVV"); // WW high purity
//        if (iCat == 1) hMass = (TH1D*) file0.Get("DijetMassLowPuriVV"); // WW low purity
//        if (iCat == 2) hMass = (TH1D*) file0.Get("DijetMassHighPuriWW"); // WW high purity
//        if (iCat == 3) hMass = (TH1D*) file0.Get("DijetMassLowPuriWW"); // WW low purity
//        if (iCat == 4) hMass = (TH1D*) file0.Get("DijetMassHighPuriWZ"); // WZ high purity
//        if (iCat == 5) hMass = (TH1D*) file0.Get("DijetMassLowPuriWZ"); // WZ low purity
//        if (iCat == 6) hMass = (TH1D*) file0.Get("DijetMassHighPuriZZ"); // ZZ high purity
//        if (iCat == 7) hMass = (TH1D*) file0.Get("DijetMassLowPuriZZ"); // ZZ low purity
       
       if(!hMass) continue;
       
       TAxis* Axis =   hMass->GetXaxis();
       for (int i = 1 ; i < hMass->GetNbinsX()+1; i++){
	 //if (hMass->GetBinCenter(i) < dMass*0.75 || hMass->GetBinCenter(i) > dMass*1.25) continue;
	 int N = abs(hMass->GetBinContent(i));
	 if (i%1000 == 0) cout << "i = " << i << " N = " << N << endl;
//.........这里部分代码省略.........

TCanvas *drawRatioPlot(TH1D *prediction, TH1D *sr, TH1D *data, TString xTitle, TString filename, double ecaloCut){

  gStyle -> SetPadLeftMargin(0.20);

  data->SetMarkerStyle(20);
  data->SetMarkerColor(kGreen);
  data->SetLineColor(kGreen);


  
  TCanvas *c = new TCanvas("c"+filename,"c",0,0,500,500);

  float y = 0.3;

  TPad *pad1     = new TPad("pad1", "Control Plots 1", 0.01, y, 0.99, 0.99);
  TPad *padRatio = new TPad("rp1", "Ratio1", 0.01, 0.01, 0.99, y-0.01);

  pad1->SetNumber(100);
  pad1->SetTicks(0,1);
  cout<<"number pad1     = "<<pad1->GetNumber()<<endl;
  cout<<"number padRatio = "<<padRatio->GetNumber()<<endl;

  TH1D *ratio = 0;

  //ratio = (TH1D*) sr->Clone();
  //ratio->Divide(prediction);

  ratio = (TH1D*) prediction->Clone();
  ratio->Divide(data);



  for(int i=1; i<=ratio->GetNbinsX();i++){

    if(ratio->GetBinContent(i) != 0){
      cout<<"N in CR in "<<i<<". bin ="<<prediction->GetBinContent(i)<<endl;
      cout<<"N in SR in "<<i<<". bin ="<<sr->GetBinContent(i)<<endl;
      cout<<"Rel. difference in "<<i<<". bin ="<<(1./ratio->GetBinContent(i)-1.)*100<<"%"<<endl; 
    }
    else if(sr->GetBinContent(i) == 0 && prediction->GetBinContent(i) !=0)    cout<<"Scaling Factor in "<<i<<". bin <"<<prediction->GetBinContent(i)/1.15<<" +/- "<<ratio->GetBinError(i)<<endl;
    else if(sr->GetBinContent(i) != 0 && prediction->GetBinContent(i) ==0)    cout<<"Scaling Factor in "<<i<<". bin <"<<(sr->GetEntries()/prediction->GetEntries())/sr->GetBinContent(i)<<" +/- "<<ratio->GetBinError(i)<<endl;

  }

  ratio->GetYaxis()->SetTitle("#frac{CR (MC)}{CR (data)}");
  ratio->SetTitle("");
  ratio->SetLabelSize(0.1,"X");
  ratio->SetLabelSize(0.1,"Y");
  ratio->SetTitleOffset(0.5,"Y");
  ratio->SetTitleSize(0.15,"Y");


  padRatio->cd();
  //ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->Draw("e");

  // Draw line at one!
  float xmin = ratio->GetXaxis()->GetXmin();
  float xmax = ratio->GetXaxis()->GetXmax();
  TLine *line = new TLine(xmin,1,xmax,1);
  line->SetLineWidth(2);
  line->Draw("same");
  padRatio->Modified();

  TLegend *leg = new TLegend(0.5,0.7,0.9,0.9);
  leg->AddEntry(sr,"SR (MC)","l"); 
  leg->AddEntry(prediction,"lepton CR (MC)","pel"); 
 
  pad1->cd();
  pad1->SetLogy();
  //  pad1->SetLogx();

  sr->SetLineColor(kRed);
  sr->SetMarkerColor(kRed);
  sr->SetMarkerStyle(20);
  sr->SetTitle("");
  prediction->SetMarkerStyle(20);
  prediction->SetTitle("");
  prediction->GetXaxis()->SetTitle(xTitle);
  sr->GetXaxis()->SetTitle(xTitle);

  
  

  prediction->SetTitleSize(0.07,"X");
  prediction->GetXaxis()->SetTitleOffset(0.7);
  sr->SetTitleSize(0.07,"X");
  sr->GetXaxis()->SetTitleOffset(0.7);

  double maximum = 0;
  double minimum = 1000000;
  if(sr->GetMinimum()!=0 && sr->GetMinimum()<minimum){
    minimum=sr->GetMinimum()*0.5;
  }
  if(prediction->GetMinimum()!=0 && prediction->GetMinimum()<minimum){
    minimum=prediction->GetMinimum()*0.5;
  }
  if(data->GetMinimum()!=0 && data->GetMinimum()<minimum){
    minimum=data->GetMinimum()*0.5;
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    hZ_EffVBFC_D->Scale(weight);
    hZ_EffVBFC_N->Scale(weight);

    hZ_DY_NoWeight->Scale(weight);
    hZ_DY_Weight->Scale(weight);
    
    // add to output histograms
    if (dataset.isData) {
      hZ_Data_C_DPhi->Add(hZ_C_DPhi);
    }
    else if (isDY) {
      hZ_DY_C_DPhi->Add(hZ_C_DPhi);

      hZ_DY_NoVBFNoWeight->Add(hZ_DY_NoWeight);
      hZ_DY_NoVBFWeight->Add(hZ_DY_Weight);

      if(dataset.name == "DYJetsToLL_NoTrig" || dataset.name == "DYJetsToLL_PtZ-100_NoTrig") {
        hZ_QCD_EffVBFS_D->Add(hZ_EffVBFS_D);
        hZ_QCD_EffVBFS_N->Add(hZ_EffVBFS_N);
        hZ_QCD_EffVBFC_D->Add(hZ_EffVBFC_D);
        hZ_QCD_EffVBFC_N->Add(hZ_EffVBFC_N);
      }
      if(dataset.name == "DYJetsToLL_EWK_NoTrig") {
        hZ_EWK_EffVBFS_D->Add(hZ_EffVBFS_D);
        hZ_EWK_EffVBFS_N->Add(hZ_EffVBFS_N);
        hZ_EWK_EffVBFC_D->Add(hZ_EffVBFC_D);
        hZ_EWK_EffVBFC_N->Add(hZ_EffVBFC_N);
      }
    }
    else {
      hZ_BG_C_DPhi->Add(hZ_C_DPhi);
    }

    std::cout << "  N ctrl (dphi<1.0) : " << hZ_C_DPhi->GetBinContent(1) << " +/- " << hZ_C_DPhi->GetBinError(1) << std::endl;  
    

    delete hZ_C_DPhi;
    delete hZ_EffVBFS_D;
    delete hZ_EffVBFS_N;
    delete hZ_EffVBFC_D;
    delete hZ_EffVBFC_N;

    delete hZ_DY_NoWeight;
    delete hZ_DY_Weight;

    ofile->cd();

    // cut flow histograms
    std::string hname = std::string("hZ_CutFlow_")+dataset.name;
    TH1D* hZ_CutFlow = new TH1D(hname.c_str(), "", nCutsZMuMu, 0., nCutsZMuMu);

    for (unsigned c=0; c<nCutsZMuMu; ++c) {

      TCut cut;

      if(c == nCutsZMuMu-1) cut = otherWeight * trigCorr * (cutD + cuts.cutflowZMuMu(c));
      else cut = otherWeight * (cutD + cuts.cutflowZMuMu(c));

      if(!(dataset.isData)) cut *= muTightWeight;

      TH1D* h = new TH1D("h","", 1, 0., 1.);
      tree->Draw("0.5>>h", cut);

      hZ_CutFlow->SetBinContent(c+1, h->GetBinContent(1));
      hZ_CutFlow->SetBinError(c+1, h->GetBinError(1));

//.........这里部分代码省略.........
  TH1D *h2Cut = new TH1D("h2Cut","",nBin,m);

  TH1D *hStat = new TH1D("hStat","",nBin,m);
  TH1D *hNpartSum = new TH1D("hNpartSum","",nBin,m);
  TH1D *hStat2 = new TH1D("hStat2","",nBin,m);
  TH1D *hNpartSum2 = new TH1D("hNpartSum2","",nBin,m);


  Float_t bin=0;
  Float_t et1=0;
  nt->SetBranchAddress("bin",&bin);
  nt->SetBranchAddress("pt1",&et1);
  
  for (int i=0;i<nt->GetEntries();i++)
  {
     nt->GetEntry(i);
     if (et1<threshold1) continue;
     
     if (et1>threshold2) {
        hNpartSum2->Fill(bin,npartValue[(int)bin]);
        hStat2->Fill(bin); 
        if (et1>threshold1) {
        hNpartSum->Fill(bin,npartValue[(int)bin]);
        hStat->Fill(bin); 
	}
     }	 
  }

  hNpartSum->Divide(hStat);
  hNpartSum2->Divide(hStat2);
  
  for (int i=1;i<nBin;i++)
  {
     cout <<hNpartSum->GetBinContent(i+1)<<endl;
     npart[i]=hNpartSum->GetBinContent(i+1);
     cout <<hNpartSum2->GetBinContent(i+1)<<endl;
     npart2[i]=hNpartSum2->GetBinContent(i+1);
     
  }

  nt->Draw("bin>>h",Form("abs(dphi)>%f&&%s",dphiCut,cut1.Data()));
  nt->Draw("bin>>hCut",Form("%s",cut1.Data()));
  TGraphAsymmErrors *g = calcEff(hCut,h,npart);
  g->SetMarkerSize(1.25);

  cout <<cut2.Data()<<endl;
  nt->Draw("bin>>h2",Form("abs(dphi)>%f&&%s",dphiCut2,cut2.Data()));
  nt->Draw("bin>>h2Cut",Form("%s",cut2.Data()));
  TGraphAsymmErrors *g2 = calcEff(h2Cut,h2,npart2);
  g2->SetMarkerSize(1.25);


  ntPythia->Draw("bin>>h",Form("abs(dphi)>%f&&%s",dphiCut2,cut2.Data()));
  ntPythia->Draw("bin>>hCut",Form("%s",cut2.Data()));
  TGraphAsymmErrors *gPythia = calcEffpythia(hCut,h,npart);
  gPythia->SetMarkerSize(1.7);
  
  if(useWeight){
    ntMix->Draw("bin>>h",Form("weight*(abs(dphi)>%f&&%s)",dphiCut2,cut2.Data()));
    ntMix->Draw("bin>>hCut",Form("weight*(%s)",cut2.Data()));
  }else{
    ntMix->Draw("bin>>h",Form("(abs(dphi)>%f&&%s)",dphiCut2,cut2.Data()));
    ntMix->Draw("bin>>hCut",Form("(%s)",cut2.Data()));
  }
  TGraphAsymmErrors *gMix = calcEff(hCut,h,npart);
  gMix->SetMarkerSize(1.25);

void tStudent()
{

   //gSystem->Load("libMathMore");
   // this is the way to force load of MathMore in Cling
   ROOT::Math::MathMoreLibrary::Load();

   int n=100;
   double a=-5.;
   double b=5.;
   //double r  = 3;
   TF1* pdf = new TF1("pdf", "ROOT::Math::tdistribution_pdf(x,3.0)", a,b);
   TF1* cum = new TF1("cum", "ROOT::Math::tdistribution_cdf(x,3.0)", a,b);

   TH1D* quant = new TH1D("quant", "", 9, 0, 0.9);

   for(int i=1; i < 10; i++)
      quant->Fill((i-0.5)/10.0, ROOT::Math::tdistribution_quantile((1.0*i)/10, 3.0 ) );

   double xx[10];
   xx[0] = -1.5;
   for(int i=1; i<9; i++)
      xx[i]= quant->GetBinContent(i);
   xx[9] = 1.5;
   TH1D* pdfq[10];
   //int nbin = n/10.0;
   for(int i=0; i < 9; i++) {
      int nbin = n * (xx[i+1]-xx[i])/3.0 + 1.0;
      TString name = "pdf";
      name += i;
      pdfq[i]= new TH1D(name, "", nbin,xx[i],xx[i+1] );
      for(int j=1; j<nbin; j++) {
         double x= j*(xx[i+1]-xx[i])/nbin + xx[i];
         pdfq[i]->SetBinContent(j, ROOT::Math::tdistribution_pdf(x,3));
      }
   }

   TCanvas *Canvas = new TCanvas("DistCanvas", "Student Distribution graphs", 10, 10, 800, 700);
   pdf->SetTitle("Student t distribution function");
   cum->SetTitle("Cumulative for Student t");
   quant->SetTitle("10-quantiles  for Student t");
   Canvas->Divide(2, 2);
   Canvas->cd(1);
   pdf->SetLineWidth(2);
   pdf->DrawCopy();
   Canvas->cd(2);
   cum->SetLineWidth(2);
   cum->SetLineColor(kRed);
   cum->Draw();
   Canvas->cd(3);
   quant->Draw();
   quant->SetLineWidth(2);
   quant->SetLineColor(kBlue);
   quant->SetStats(0);
   Canvas->cd(4);
   pdfq[0]->SetTitle("Student t & its quantiles");
   pdf->SetTitle("");
   pdf->Draw();
   //pdfq[0]->SetAxisRange(-1.5, 0, 1.5,1.0);
   pdfq[0]->SetTitle("Student t & its quantiles");
   for(int i=0; i < 9; i++) {
      pdfq[i]->SetStats(0);
      pdfq[i]->SetFillColor(i+1);
      pdfq[i]->Draw("same");
   }
   Canvas->Modified();
   Canvas->cd();
}

//.........这里部分代码省略.........
    if(!hDY || !hZinv || !hData){
      cout<<"could not find histogram "<<fullhistname<<endl;
      continue;
    }
    
    if (hDY->GetNbinsX() != hZinv->GetNbinsX() ) {
      cout<<"different binning for histograms "<<fullhistname<<endl;
      continue;
    }
    
    // Make directory and plot(s) in the output file
    TDirectory* dir = 0;
    dir = (TDirectory*)outfile->Get(directory.Data());
    if (dir == 0) {
      dir = outfile->mkdir(directory.Data());
    } 
    dir->cd();

    cout<<"Looking at histo "<<fullhistname<<endl;

    int lastbin_hybrid = 1;
    int lastmt2val_hybrid = 200;
    int ht_LOW = 0, ht_HI = 0, njets_LOW = 0, njets_HI = 0, nbjets_LOW = 0, nbjets_HI = 0;
    TH1D* h_MT2Template = (TH1D*) hZinv->Clone("h_MT2Template");  
    TString inclusiveTemplateName = "";

    TH1D *h_ht_LOW, *h_ht_HI, *h_njets_LOW, *h_njets_HI, *h_nbjets_LOW, *h_nbjets_HI;
    if (doHybridInclusiveTemplate) {
      //  Inclusive template: use inclusive template corresponding to this region
      
      //Get variable boundaries for signal region.
      h_ht_LOW = (TH1D*) fData->Get(directory+"/h_ht_LOW");
      h_ht_HI  = (TH1D*) fData->Get(directory+"/h_ht_HI");
      if (h_ht_LOW) ht_LOW = h_ht_LOW->GetBinContent(1);
      if (h_ht_HI)  ht_HI = h_ht_HI->GetBinContent(1);
      h_njets_LOW = (TH1D*) fData->Get(directory+"/h_njets_LOW");
      h_njets_HI  = (TH1D*) fData->Get(directory+"/h_njets_HI");
      if (h_njets_LOW) njets_LOW = h_njets_LOW->GetBinContent(1);
      if (h_njets_HI)  njets_HI = h_njets_HI->GetBinContent(1);
      h_nbjets_LOW = (TH1D*) fData->Get(directory+"/h_nbjets_LOW");
      h_nbjets_HI  = (TH1D*) fData->Get(directory+"/h_nbjets_HI");
      if (h_nbjets_LOW) nbjets_LOW = h_nbjets_LOW->GetBinContent(1);
      if (h_nbjets_HI)  nbjets_HI = h_nbjets_HI->GetBinContent(1);
      
      //Determine which inclusive template to use. If none works, this reverts to HybridSimple, taking template from its own TopoRegion 
      // Start from the Aggregate Regions (hardcoded, since they can partially overlap with the standard regions)
      if (srname == "20") inclusiveTemplateName = "crdy20/h_mt2bins";  // self (2j, HT1200)
      else if (srname == "base") inclusiveTemplateName = "crdybase/h_mt2bins";
      else if (srname == "baseVL") inclusiveTemplateName = "crdybaseVL/h_mt2bins";
      else if (srname == "baseL") inclusiveTemplateName = "crdybaseL/h_mt2bins";
      else if (srname == "baseM") inclusiveTemplateName = "crdybaseM/h_mt2bins";
      else if (srname == "baseH") inclusiveTemplateName = "crdybaseH/h_mt2bins";
      else if (srname == "baseUH") inclusiveTemplateName = "crdybaseUH/h_mt2bins";
      else if (srname == "21") inclusiveTemplateName = "crdy21/h_mt2bins"; // self (2j, HT1500)
      else if (srname == "22") inclusiveTemplateName = "crdy22/h_mt2bins"; // self (4j, HT1200)
      else if (srname == "23") inclusiveTemplateName = "crdy21/h_mt2bins"; // from 21 
      else if (srname == "24") inclusiveTemplateName = "crdy24/h_mt2bins"; // self (7J, HT1200)
      else if (srname == "25") inclusiveTemplateName = "crdy21/h_mt2bins"; // from 21
      else if (srname == "26") inclusiveTemplateName = "crdy20/h_mt2bins"; // from 20
      else if (srname == "27") inclusiveTemplateName = "crdy21/h_mt2bins"; // from 21
      else if (srname == "28") inclusiveTemplateName = "crdy20/h_mt2bins3J"; // need a 3J histogram within SR20
      //      else if (srname == "28") inclusiveTemplateName = "crdy20/h_mt2bins"; // test
      else if (srname == "29") inclusiveTemplateName = "crdy21/h_mt2bins"; // from 21
      else if (srname == "30") inclusiveTemplateName = "crdy24/h_mt2bins"; // from 24
      else if (srname == "31") inclusiveTemplateName = "crdy21/h_mt2bins"; // from 21
      // Now the standard regions