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

void oneDlinearcombination(TH1F* first, int firsttype, TH1F* second, int secondtype, TH1F* input, int inputtype, TH1F* finaloutput, int outputtype, TH1F* finaloutput2, int output2type){
	TH1F* output = (TH1F*) input->Clone();
	TH1F* output2 = (TH1F*) output->Clone();
	if(outputtype==kIntHist || outputtype==kSigHist){
		if(inputtype==kBSI25Hist && firsttype==kSigHist && secondtype==kBkgHist){
			output->Add(first, -25.0);
			output->Add(second, -1.0);
			output->Scale(0.2);
			if (outputtype == kSigHist){ delete output; output = (TH1F*) first->Clone(); }

			output2->Add(first, -25.0);
			output2->Add(second, -1.0);
			output2->Scale(0.2);
			if (output2type == kSigHist){ delete output2; output2 = (TH1F*) first->Clone(); }
		}
		if(inputtype==kBSI25Hist && firsttype==kBSIHist && secondtype==kBkgHist){
			for (int binx = 1; binx <= output->GetNbinsX(); binx++){
				double bsi = first->GetBinContent(binx);
				double bkg = second->GetBinContent(binx);
				double bsi25 = output->GetBinContent(binx);

				double weight=doLinearCombination(bsi25,25,bsi,1,bkg,outputtype);
				output->SetBinContent(binx,weight);
				if (finaloutput2 != 0){
					double weight2 = doLinearCombination(bsi25, 25, bsi, 1, bkg, output2type);
					output2->SetBinContent(binx, weight2);
				}
			}
		}		
		if(inputtype==kBSI25Hist && firsttype==kBSI10Hist && secondtype==kBkgHist){
			for (int binx = 1; binx <= output->GetNbinsX(); binx++){
				double bsi10 = first->GetBinContent(binx);
				double bkg = second->GetBinContent(binx);
				double bsi25 = output->GetBinContent(binx);

				double weight=doLinearCombination(bsi25,25,bsi10,10,bkg,outputtype);
				output->SetBinContent(binx,weight);
				if (finaloutput2 != 0){
					double weight2 = doLinearCombination(bsi25, 25, bsi10, 10, bkg, output2type);
					output2->SetBinContent(binx, weight2);
				}
			}
		}
		if(inputtype==kBSI10Hist && firsttype==kBSIHist && secondtype==kBkgHist){
			for (int binx = 1; binx <= output->GetNbinsX(); binx++){
				double bsi = first->GetBinContent(binx);
				double bkg = second->GetBinContent(binx);
				double bsi10 = output->GetBinContent(binx);

				double weight=doLinearCombination(bsi10,10,bsi,1,bkg,outputtype);
				output->SetBinContent(binx,weight);
				if (finaloutput2 != 0){
					double weight2 = doLinearCombination(bsi10, 10, bsi, 1, bkg, output2type);
					output2->SetBinContent(binx, weight2);
				}
			}
		}
		for (int binx = 1; binx <= output->GetNbinsX(); binx++){
			finaloutput->SetBinContent(binx,output->GetBinContent(binx));
			if (finaloutput2 != 0) finaloutput2->SetBinContent(binx,output2->GetBinContent(binx));
		}
	}
	else{cout<<"Option not yet supported. Exiting..."<<endl; assert(0);};
	delete output;
	delete output2;
}

//.........这里部分代码省略.........
	// -------------------------------
	// fit range without data entries
	// -------------------------------
	// NB: needed if loaded data is interpreted as points with nothing 
	//     between (like kidonakis)
	if(points){
	  // perform a linear fit wrt previous point
	  // get the two points (this bin and the previous one)
	  int binPrev= (bin==1&&variable=="topPt") ? 0 : hist->FindBin(Xvalues[bin-2]+0.5*binwidth);
	  double x2=-1;
	  double x1= 0;
	  double y2=-1;
	  double y1= 0.;
	  rawHist->GetPoint(bin-1, x2, y2);
	  x2+=0.5*binwidth;
	  if(bin==1&&variable=="topPt"){
	    y1=0;
	    x1=0;
	  }
	  else{
	    rawHist->GetPoint(bin-2, x1, y1);
	    x1+=0.5*binwidth;
	  }
	  // calculate linear funtion
	  double a=(y2-y1)/(x2-x1);
	  double b=y1-a*x1;
	  TF1* linInterpol=new TF1("linInterpol"+getTStringFromInt(bin), "[0]*x+[1]", x1, x2);
	  linInterpol->SetParameter(0,a);
	  linInterpol->SetParameter(1,b);
	  double xlow  = (bin==1&&variable=="topPt") ? 0. : hist->GetBinLowEdge(binPrev+1);
	  double xhigh = hist->GetBinLowEdge(bin2+1);
	  std::cout << " lin. interpolation [ (" << x1 << "," << y1 << ") .. (" << x2 << "," << y2 << ") ]: " << a << "*x+" << b << std::endl;
	  linInterpol->SetRange(xlow, xhigh);
	  hist->Add(linInterpol);
	}
      }
    }
    // check theory curve
    std::cout << std::endl << "analyze theory curve:" << std::endl;
    double integralRawTheory=hist->Integral(0.,hist->GetNbinsX()+1);
    std::cout << "Integral: " << integralRawTheory << std::endl;
    if(integralRawTheory==0.){
      std::cout << "ERROR: Integral can not be 0!" << std::endl;
      exit(0);
    }
    std::cout << "binwidth: " << binwidth << std::endl;
    std::cout << "Integral*binwidth: " << integralRawTheory*binwidth << std::endl;
    std::cout << "binRange: " << xmin << ".." << xmax << std::endl;
    std::cout << "   -> reco range: " << low << ".." << high << std::endl;
    std::vector<double> recoBins_=bins_[variable];
    for(unsigned int i=0; i<recoBins_.size(); ++i){
      i==0 ? std::cout << "(" : std::cout << ", ";
      std::cout << recoBins_[i];
      if(i==recoBins_.size()-1) std::cout << ")" << std::endl;
    }
    // check if you need to divide by binwidth
    if(std::abs(1.-integralRawTheory)<0.03){
      std::cout << "Integral is approx. 1 -> need to divide by binwidth!" << std::endl;
      hist->Scale(1./binwidth);
    }
    else if(std::abs(1-integralRawTheory*binwidth)<0.03){
      std::cout << "Integral*binwidth is approx. 1 -> no scaling needed!" << std::endl;
    }
    else{
      std::cout << "need to normalize and divide by binwidth";
      hist->Scale(1./(binwidth*integralRawTheory));

void ReactorNuAnalysis()
{
	// -------- VARIABLE DEFINITIONS and SETUP --------
	// Style settings
	gStyle->SetOptStat("");
	
	// Constants
	Double_t NuSpectrumMinE = 0.0;  // [MeV]
	Double_t NuSpectrumMaxE = 10.0; // [MeV]
	Int_t seed = 43534;
	const Double_t Gfermi = 1.16637e-11;																			 // [MeV^-2]
	const Double_t Sin2ThetaW = 0.2387;
	const Double_t InverseMeVtoCm = 1.97e-11;
	const Double_t elementaryCharge = 1.602176565e-19;
	
	// Spectrum files (these are just samples for now)
	const char ReactorNeutronBackgroundFile[] = "SampleData.txt";
	
	// Define constants relevant for this run
	Double_t OnOffTimeRatio = 1.0/1.0;
	Double_t time = 100 * 24.0*3600.0;																				 // [sec]
	Double_t detMass = 5000.0;																								 // [g]
	Double_t distance = 200.0;																								 // [cm]
	Double_t activity = 6.0/200.0/elementaryCharge;														 // [sec^-1]
	const Int_t nNeutrons = 14;
	const Int_t nProtons = 14;
	const Double_t Qweak = nNeutrons - (1 - 4*Sin2ThetaW) * nProtons;
	const Double_t NucleonMass = (nNeutrons * 939.565) + (nProtons * 938.272);			 // [MeV]
	
	// Define the histograms
	TH1F* RecoilEvtHistogram = new TH1F("RecoilEvtHistogram","^{28}Si recoil energy spectrum;Energy [eV];Events / 10 eV",50,0.0,500.0);
	TH1F* NeutrinoEvtHistogram = new TH1F("NeutrinoEvtHistogram","#bar{#nu}_{e} energy spectrum;Energy [MeV];Events / 0.5 MeV",50,0.0,20);
	TH1F* TheoryRecoilEvtHistogram = new TH1F("TheoryRecoilEvtHistogram","High-statistics (\"theoretical\") Coherent Recoil Spectrum;Energy [eV];Events / 10 eV",50,0.0,500.0);
	TH1F* TheoryNeutrinoEvtHistogram = new TH1F("TheoryNeutrinoEvtHistogram","High-statistics (\"theoretical\") Neutrino Energy Spectrum;Energy [MeV] / 0.5 MeV;Events",50,0.0,20);
	TH1F* TheoryRecoilEvtHistogramFit = new TH1F("TheoryRecoilEvtHistogram","MC Fit;Events",50,0.0,0.002);
	TH1F* EMNoLukeBackground = new TH1F("EMNoLukeBackground","Background from EM recoils (before Luke Effect amplification;Energy [MeV];Events)",50,0.0,0.002);
	TH1F* ReactorNeutronBackground = new TH1F("ReactorNeutronBackground","Background from reactor neutrons;Energy [MeV];Events",50,0.0,0.002);
	TH1F* ReactorOnCosmoNeutronBackground = new TH1F("ReactorOnCosmoNeutronBackground","Reactor-on background from muon-induced neutrons;Energy [MeV];Events)",50,0.0,0.002);
	TH1F* ReactorOffCosmoNeutronBackground = new TH1F("ReactorOffCosmoNeutronBackground","Reactor-off background from muon-induced neutrons;Energy [MeV];Events)",50,0.0,0.002);	
	TH1F* ReactorOnHisto = new TH1F("ReactorOnHisto","Recoil Spectrum for Reactor-On Data;Energy [MeV];Events",50,0.0,0.002);
	TH1F* ReactorOffHisto = new TH1F("ReactorOffHisto","Recoil Spectrum for Reactor-Off Data;Energy [MeV];Events",50,0.0,0.002);
	TH1F* BackgroundSubtractedSignal = new TH1F("BackgroundSubtractedSignal","Recoil Spectrum for Reactor-Off Data;Energy [MeV];Events",50,0.0,0.002);
	
	// Energy spectra
	// (Spectral parameterizations from arXiv:1101.2663v3)
	TF1* NeutrinoEnergySpectrum = new TF1("NeutrinoSpectrum", "TMath::Exp([0] + [1]*x + [2]*TMath::Power(x,2) + [3]*TMath::Power(x,3) + [4]*TMath::Power(x,4) + [5]*TMath::Power(x,5))", NuSpectrumMinE, NuSpectrumMaxE);
	NeutrinoEnergySpectrum->SetParameters(3.217, -3.111, 1.395, -0.369, 0.04445, -0.002053);
	TF1* IntegratedRecoilSpectrum = new TF1("IntegratedRecoilSpectrum", "TMath::Power([0]*[1]*[3],2)/(4*TMath::Pi()) * [2] * (x/(1 + [2]/(2*x))) * (1 - ([2]/(4*x*x)) * (x/(1 + [2]/(2*x))))", NuSpectrumMinE, NuSpectrumMaxE);
	IntegratedRecoilSpectrum->SetParameters(Gfermi, Qweak, NucleonMass, InverseMeVtoCm);
	TF1* RecoilSpectrum = new TF1("RecoilSpectrum","IntegratedRecoilSpectrum * NeutrinoSpectrum", 0.0, 10.0);
	TF1* DiffRecoilSpectrumAtConstE = new TF1("DiffRecoilSpectrumAtConstE", "(x<[4])*TMath::Power([0]*[1],2)/(4*TMath::Pi()) * [2] * (1 - ([2] * x)/(2 * TMath::Power([3],2))) + (x>[4])*0", NuSpectrumMinE, 0.01);
	DiffRecoilSpectrumAtConstE->SetParameters(Gfermi, Qweak, NucleonMass);
	
	
	
	
	// -------- HISTOGRAM FILLING --------
	// Fill "experimental" histograms
	Int_t nEvt = GenerateNumOfNuRecoils(time, detMass, distance, activity, nNeutrons, nProtons, RecoilSpectrum, NuSpectrumMinE, NuSpectrumMaxE, seed);
	FillNuRecoilSpectrum(RecoilEvtHistogram, NeutrinoEvtHistogram, nEvt, nNeutrons, nProtons, RecoilSpectrum, DiffRecoilSpectrumAtConstE, NuSpectrumMinE, NuSpectrumMaxE, seed+4);
	FillRecoilSpectrumFromFile(ReactorNeutronBackground, 100.0, 10.0, ReactorNeutronBackgroundFile, seed+1);      // Testing purposes only.  CHANGE ME!!
	FillRecoilSpectrumFromFile(ReactorOnCosmoNeutronBackground, 50.0, 10.0, ReactorNeutronBackgroundFile, seed+2);		// Testing purposes only.  CHANGE ME!!
	FillRecoilSpectrumFromFile(ReactorOffCosmoNeutronBackground, 50.0, 10.0, ReactorNeutronBackgroundFile, seed+3);		// Testing purposes only.  CHANGE ME!!
	
	cout << "nEvt: " << nEvt << endl;
	
	// Fill high-statistics "theoretical" histograms
	FillNuRecoilSpectrum(TheoryRecoilEvtHistogram, TheoryNeutrinoEvtHistogram, 10000, nNeutrons, nProtons, RecoilSpectrum, DiffRecoilSpectrumAtConstE, NuSpectrumMinE, NuSpectrumMaxE, seed+5);
	TheoryNeutrinoEvtHistogram->Scale(nEvt/TheoryNeutrinoEvtHistogram->GetEntries());
	TheoryRecoilEvtHistogram->Scale(nEvt/TheoryNeutrinoEvtHistogram->GetEntries());
	
	// Combine the histograms into total
	ReactorOnHisto->Add(RecoilEvtHistogram);
	ReactorOnHisto->Add(ReactorNeutronBackground);
	ReactorOnHisto->Add(ReactorOnCosmoNeutronBackground);
	
	ReactorOffHisto->Add(ReactorOffCosmoNeutronBackground);
	
	
	
	
	// -------- HYPOTHESIS TESTING --------
	cout << "p-value between Reactor-On and Reactor-Off Data: " << ReactorOnHisto->Chi2Test(ReactorOffHisto) << endl;
	cout << "p-value between the simulated data and the Monte Carlo histogram: " << RecoilEvtHistogram->Chi2Test(TheoryRecoilEvtHistogram) << endl;
	
	
	
	
	// -------- BACKGROUND SUBTRACTION and FITTING --------
	// Normalize reactor-off data to reactor-on data by exposure
	BackgroundSubtractedSignal->Add(ReactorOnHisto, ReactorOffHisto, 1.0, -1.0*OnOffTimeRatio);
	
	// Use TFractionFitter to do fitting
	TObjArray *FractionFitData = new TObjArray(2);
	FractionFitData->Add(TheoryRecoilEvtHistogram);
	FractionFitData->Add(ReactorOffHisto);
	TFractionFitter* ffit = new TFractionFitter(ReactorOnHisto, FractionFitData);
	ffit->Constrain(0,0.1,10.0);
	ffit->Constrain(1,1.0,1.0);
	Int_t status = ffit->Fit();
//.........这里部分代码省略.........

void DeltaEta_norm(){


  TLegend* leg = new TLegend(0.70,0.55,0.89,0.89);
  TFile *data =  TFile::Open("./50ns_betaV4/DoubleEG.root","READ");
  TFile *sig_M1 = TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M1.root","READ");
  TFile *sig_M10 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M10.root","READ");  
  TFile *sig_M100 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M100.root","READ");
  TFile *sig_M1000 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M1000.root","READ");
  TFile *bkg1 =  TFile::Open("./50ns_betaV4/NewWeightDiPhoton.root","READ"); 
  TFile *bkg2 =  TFile::Open("./50ns_betaV4/NewWeightDYJetsToLL.root","READ");
  TFile *bkg3 =  TFile::Open("./50ns_betaV4/NewWeightGJets.root","READ");  
  TFile *bkg4 =  TFile::Open("./50ns_betaV4/NewWeightGluGluHToGG.root","READ");
  TFile *bkg5 =  TFile::Open("./50ns_betaV4/NewWeightQCD.root","READ");
  TFile *bkg6 =  TFile::Open("./50ns_betaV4/NewWeightVH.root","READ");


 
  TTree *tree_data = (TTree*) data->Get("DiPhotonTree");
  TTree *tree_sig_M1 = (TTree*) sig_M1->Get("WeightDiPhotonTree"); 
  TTree *tree_sig_M10 = (TTree*) sig_M10->Get("WeightDiPhotonTree");
  TTree *tree_sig_M100 = (TTree*) sig_M100->Get("WeightDiPhotonTree");
  TTree *tree_sig_M1000 = (TTree*) sig_M1000->Get("WeightDiPhotonTree");
  TTree *tree_bkg1 = (TTree*) bkg1->Get("WeightDiPhotonTree");
  TTree *tree_bkg2 = (TTree*) bkg2->Get("WeightDiPhotonTree");
  TTree *tree_bkg3 = (TTree*) bkg3->Get("WeightDiPhotonTree");
  TTree *tree_bkg4 = (TTree*) bkg4->Get("WeightDiPhotonTree");
  TTree *tree_bkg5 = (TTree*) bkg5->Get("WeightDiPhotonTree");
  TTree *tree_bkg6 = (TTree*) bkg6->Get("WeightDiPhotonTree");
  
  TCanvas *c1n = new TCanvas("c1n","",1200,400);
  
  c1n->Divide(2,1);  
  c1n->cd(1);
  gPad->SetLogy();  
  
   
 
  tree_data->Draw("(fabs(eta1-eta2))>>hdata(30,0,3)");
  TH1F *hdata =(TH1F*)gPad->GetPrimitive("hdata");
  tree_sig_M1->Draw("(fabs(eta1-eta2))>>h1(30,0,3)","weightNew");
  TH1F *h1 =(TH1F*)gPad->GetPrimitive("h1"); 
  tree_sig_M10->Draw("(fabs(eta1-eta2))>>h2(30,0,3)","weightNew");  
  TH1F *h2 =(TH1F*)gPad->GetPrimitive("h2"); 
  tree_sig_M100->Draw("(fabs(eta1-eta2))>>h3(30,0,3)","weightNew");  
  TH1F *h3 =(TH1F*)gPad->GetPrimitive("h3");
  tree_sig_M1000->Draw("(fabs(eta1-eta2))>>h4(30,0,3)","weightNew");  
  TH1F *h4 =(TH1F*)gPad->GetPrimitive("h4");
  tree_bkg1->Draw("(fabs(eta1-eta2))>>hbkg1(30,0,3)","weight*weightNew");
  TH1F *hbkg1 =(TH1F*)gPad->GetPrimitive("hbkg1");
  tree_bkg2->Draw("(fabs(eta1-eta2))>>hbkg2(30,0,3)","weight*weightNew");  
  TH1F *hbkg2 =(TH1F*)gPad->GetPrimitive("hbkg2");
  tree_bkg3->Draw("(fabs(eta1-eta2))>>hbkg3(30,0,3)","weight*weightNew");  
  TH1F *hbkg3 =(TH1F*)gPad->GetPrimitive("hbkg3");
  tree_bkg4->Draw("(fabs(eta1-eta2))>>hbkg4(30,0,3)","weight*weightNew");  
  TH1F *hbkg4 =(TH1F*)gPad->GetPrimitive("hbkg4");
  tree_bkg5->Draw("(fabs(eta1-eta2))>>hbkg5(30,0,3)","weight*weightNew");
  TH1F *hbkg5 =(TH1F*)gPad->GetPrimitive("hbkg5");
  tree_bkg6->Draw("(fabs(eta1-eta2))>>hbkg6(30,0,3)","weight*weightNew");  
  TH1F *hbkg6 =(TH1F*)gPad->GetPrimitive("hbkg6");
 
  

  
  h1->SetLineColor(kRed-4);
  h2->SetLineColor(kRed+2);  
  h3->SetLineColor(kPink+2);
  h4->SetLineColor(kPink+4);
  h1->SetLineWidth(2);
  h2->SetLineWidth(2);  
  h3->SetLineWidth(2);
  h4->SetLineWidth(2);

  THStack *hs= new THStack("hs","");
  hbkg6->SetFillColor(kGreen+3);  
  hbkg4->SetFillColor(kGreen+4);
  hbkg1->SetFillColor(kBlue);
  hbkg2->SetFillColor(kBlue+2);
  hbkg3->SetFillColor(kMagenta+4);
  hbkg5->SetFillColor(kMagenta+2);

  hbkg1->SetLineColor(kBlack);
  hbkg2->SetLineColor(kBlack);
  hbkg3->SetLineColor(kBlack);
  hbkg4->SetLineColor(kBlack);
  hbkg5->SetLineColor(kBlack);
  hbkg6->SetLineColor(kBlack);

  hs->Add(hbkg6);
  hs->Add(hbkg4);
  hs->Add(hbkg1);
  hs->Add(hbkg2);
  hs->Add(hbkg3); 
  hs->Add(hbkg5);
  
  TH1F *hsum = (TH1F*)hbkg1->Clone("hsum");
  hsum->Add(hbkg2);
  hsum->Add(hbkg3);
  hsum->Add(hbkg4);
  hsum->Add(hbkg5);
//.........这里部分代码省略.........

/*************************************************************************************
 * createMFOutfile: moves the MLB distributions into an output file for use in 
 *                  R. Nally's MassFit.C code.
 *  input:  the main() arguments array
 *  output: writes to an output file the histograms, in a MassFit.C-readable format 
 *
 *  Structure-wise: can be implemented into class easily.
 ***********************************/
void createMFOutfile(const char* argv[]) {
  TFile *f = new TFile(argv[2]);

  //create the output file we'd like to write histograms to
  TFile *output = new TFile(outfileName, "RECREATE");
  
  //get the filesystem information from the file
  f->cd();
  TList *alokDirs = (TList*) f->GetListOfKeys();

  //create a list of "All" histograms and initialize (TODO)
  TList *allHists = new TList;

  //loop through the directories in the input file
  for(int idir=0; alokDirs->At(idir-1) != alokDirs->Last(); idir++) {
    //if it's not mlb, we don't care
    if(!TString(alokDirs->At(idir)->GetName()).Contains("_Mlb")) continue;
    
    //get the file directory information, its histograms
    TDirectory *cDir  = (TDirectory*) f->Get(alokDirs->At(idir)->GetName());
    TList *alokHistos = (TList*)      cDir->GetListOfKeys();

    //loop through the histograms in the current directory 
    for(int ihisto=0; alokHistos->At(ihisto-1) != alokHistos->Last(); ihisto++) {
      // don't consider the graph objects
      if(TString(alokHistos->At(ihisto)->GetName()).Contains("Graph_")) continue; 

      // clone the histogram, give it a new name
      TString cloneName = formatName(alokHistos->At(ihisto)->GetName(),nominalWidth);
      TH1F *thisto = (TH1F*) cDir->Get(alokHistos->At(ihisto)->GetName());
      TH1F *tclone = (TH1F*) thisto->Clone(cloneName);

      // open the outfile and write the histo in
      output->cd();
      TList *outkeys = (TList*) output->GetListOfKeys();

      // if the histogram already exists, add thisto to the existing one
      // messy, but removes the need for magic
      for(int iout=0; outkeys->At(iout-1) != outkeys->Last(); iout++) {
        if(outkeys->At(iout)->GetName() == cloneName) {
          cout<<" - found another histogram in output with name "<<cloneName<<endl;
          TH1F *theHisto = (TH1F*) output->Get(cloneName);
          cout<<" - got the same histogram from the output file"<<endl;
          TH1F *tHclone  = (TH1F*) theHisto->Clone(cloneName);
          cout<<" - cloned the histogram"<<endl;

          cout<<" - adding in clone"<<endl;
          tclone->Add(tHclone);

          cout<<" - deleting the original histogram from the output file"<<endl;
          output->Delete(cloneName + TString(";1"));
          cout<<" - deleted thing from output file"<<endl;
          cout<<" - tclone looks like "<<tclone<<endl;
        }
      }
      cout<<" - writing the tclone to file"<<endl;
      tclone->Write();

      // reopen the input root file and clean up
      f->cd();
      delete thisto;
      delete tclone;
    }
  }

  f->Close();
  output->cd();
  output->Close();

  // if we want to interpolate, start making more histograms
  if(interpolate) {
    // Have to reopen the outfile because ROOT is odd
    TFile *output2 = new TFile(outfileName, "UPDATE");
    output2->cd();

    std::pair<double,TString> maxPair = *moreFiles.begin();
    double maxWidth = maxPair.first;

    // check that it makes sense to interpolate with our settings
    if(maxWidth <= nominalWidth) {
      cout<<"\n\nERROR: Max width is less than or equal to the nominal width! Exiting."<<endl;
      exit(EXIT_FAILURE);
    }

    // open the input file and loop through the relevant directories
    TFile *maxFile = new TFile(maxPair.second);
    for(int i=0; i<lepsSize; i++) {
      // change directory and get the name of the folder we want to access
      maxFile->cd();
      char a[128];
      sprintf(a, "mlbwa_%s_Mlb", leps[i]);
      TDirectory *tDir = (TDirectory*) maxFile->Get(a);
//.........这里部分代码省略.........

void quickLifetimeFit(string filename, bool isB0)
{
    setTDRStyle();

    TFile *_file0 = TFile::Open(filename.c_str());

    TTree *tree = (TTree*)gDirectory->Get("fittree");

    TCanvas *c = new TCanvas("c","c",600,800);
    c->Divide(2,3);
    int canvasctr(0);

    string strSigWindow, strBgrWindow;
    double SoverSB;

    if (isB0)
    {
	strSigWindow = "mass>5.25&&mass<5.31";
	strBgrWindow = "mass>5.35&&mass<5.55";
	SoverSB = 0.90; // from another fit
    }
    else
    {
	strSigWindow = "mass>5.60&&mass<5.64";
	strBgrWindow = "mass>5.75&&mass<6.30";
	SoverSB = 0.60; // from another fit
    }

    c->cd(++canvasctr);
    tree->Draw("mass>>hmsigwindow", strSigWindow.c_str());
    TH1F *hmsigwindow = (TH1F*)gDirectory->GetList()->FindObject("hmsigwindow");
    const int nSigWindow = hmsigwindow->GetEntries();

    c->cd(++canvasctr);
    tree->Draw("mass>>hbgrwindow", strBgrWindow.c_str());
    TH1F *hbgrwindow = (TH1F*)gDirectory->GetList()->FindObject("hbgrwindow");
    const int nBgrWindow = hbgrwindow->GetEntries();

    const double nSig = nSigWindow * SoverSB;
    const double nBgr = nSigWindow * (1.-SoverSB);
    const double scale = (double)nBgr / (double)nBgrWindow;

    c->cd(++canvasctr);
    tree->Draw("t>>hsigt(34,-2e-12,15e-12)", strSigWindow.c_str());
    TH1F *hsigt = (TH1F*)gDirectory->GetList()->FindObject("hsigt");
    gPad->SetLogy();

    c->cd(++canvasctr);
    tree->Draw("t>>hbgrt(34,-2e-12,15e-12)", strBgrWindow.c_str());
    TH1F *hbgrt = (TH1F*)gDirectory->GetList()->FindObject("hbgrt");
    hbgrt->Scale(scale);
    gPad->SetLogy();

    c->cd(++canvasctr);
    hsigt->Draw();
    hsigt->SetLineColor(4);
    hbgrt->Draw("same");
    hbgrt->SetLineColor(2);
    gPad->SetLogy();

    c->cd(++canvasctr);
    TH1F *hdiff = new TH1F("hdiff","hdiff",34,-2e-12,15e-12);
    hdiff->Add(hsigt,hbgrt,1,-1);
    hdiff->Draw();
    gPad->SetLogy();

    c->SaveAs("quickFit.pdf");

    TCanvas *c2 = new TCanvas("c2","c2",1000,800);
    hdiff->Draw();
    gPad->SetLogy();

}

void 
//HHH_TT_X_1or2tag(bool scaled=true, bool log=true, float min=0.1, float max=-1., string inputfile="root/$HISTFILE", const char* directory="tauTau_$CATEGORY")
HHH_TT_X_1or2tag(bool scaled=true, bool log=true, float min=0.1, float max=-1., string inputfile="root/$HISTFILE", const char* directory="tauTau_$CATEGORY")
{
  // defining the common canvas, axes pad styles
  SetStyle(); gStyle->SetLineStyleString(11,"20 10");

  // determine category tag
  const char* category = ""; const char* category_extra = ""; const char* category_extra2 = "";
  if(std::string(directory) == std::string("tauTau_2jet0tag")){ category = "#tau_{h}#tau_{h}";           }
  if(std::string(directory) == std::string("tauTau_2jet0tag")){ category_extra= "2-jet 0 b-tag";           }
  if(std::string(directory) == std::string("tauTau_2jet1tag"  )){ category = "#tau_{h}#tau_{h}";           }
  if(std::string(directory) == std::string("tauTau_2jet1tag"  )){ category_extra= "2-jet 1 b-tag";     }
  if(std::string(directory) == std::string("tauTau_2jet2tag"  )){ category = "#tau_{h}#tau_{h}";           }
  if(std::string(directory) == std::string("tauTau_2jet2tag"  )){ category_extra = "2-jet 2 b-tag";              }

  const char* dataset;
#ifdef MSSM
  if(std::string(inputfile).find("7TeV")!=std::string::npos){dataset = "#scale[1.5]{CMS}  h,H,A#rightarrow#tau#tau                                 4.9 fb^{-1} (7 TeV)";}
  if(std::string(inputfile).find("8TeV")!=std::string::npos){
    if(std::string(directory).find("btag")!=std::string::npos){
      dataset = "#scale[1.5]{CMS}  h,H,A#rightarrow#tau#tau                                18.3 fb^{-1} (8 TeV)";
    }
    else{
      dataset = "#scale[1.5]{CMS}  h,H,A#rightarrow#tau#tau                                19.7 fb^{-1} (8 TeV)";
    }
  }
#else
  if(std::string(inputfile).find("8TeV")!=std::string::npos){dataset = "CMS, 19.7 fb^{-1} at 8 TeV";}
#endif
  
  // open example histogram file
  TFile* input = new TFile(inputfile.c_str());
#ifdef MSSM
  TFile* input2 = new TFile((inputfile+"_$MH_$TANB").c_str());
#endif
  TH1F* Fakes  = refill((TH1F*)input->Get(TString::Format("%s/QCD"     , directory)), "QCD"); InitHist(Fakes, "", "", TColor::GetColor(250,202,255), 1001);
//  TH1F* EWK1   = refill((TH1F*)input->Get(TString::Format("%s/W"       , directory)), "W"  ); InitHist(EWK1 , "", "", TColor::GetColor(222,90,106), 1001);
  TH1F* EWK2   = refill((TH1F*)input->Get(TString::Format("%s/ZLL"      , directory)), "ZLL" ); InitHist(EWK2 , "", "", TColor::GetColor(222,90,106), 1001);
//TH1F* EWK3   = refill((TH1F*)input->Get(TString::Format("%s/ZL"      , directory)), "ZL" ); InitHist(EWK3 , "", "", TColor::GetColor(222,90,106), 1001);
  TH1F* EWK    = refill((TH1F*)input->Get(TString::Format("%s/VV"      , directory)), "VV" ); InitHist(EWK  , "", "", TColor::GetColor(222,90,106), 1001);
  TH1F* ttbar  = refill((TH1F*)input->Get(TString::Format("%s/TT"      , directory)), "TT" ); InitHist(ttbar, "", "", TColor::GetColor(155,152,204), 1001);
  TH1F* Ztt    = refill((TH1F*)input->Get(TString::Format("%s/ZTT"     , directory)), "ZTT"); InitHist(Ztt  , "", "", TColor::GetColor(248,206,104), 1001);
#ifdef MSSM
  TH1F* ggHTohhTo2Tau2B    = refill((TH1F*)input2->Get(TString::Format("%s/ggHTohhTo2Tau2B$MH" , directory)), "ggHTohhTo2Tau2B"); InitSignal(ggHTohhTo2Tau2B); ggHTohhTo2Tau2B->Scale($TANB*SIGNAL_SCALE);
 /* if(std::string(directory)=="tauTau_2jet1tag") ggHTohhTo2Tau2B->Scale(0.5);
  if(std::string(directory)=="tauTau_2jet2tag") ggHTohhTo2Tau2B->Scale(0.3);
  TH1F* ggH_SM125 = refill((TH1F*)input->Get(TString::Format("%s/ggH_SM125",directory)),"ggH_SM125");InitSignal(ggH_SM125);ggH_SM125->Scale(SIGNAL_SCALE);
  TH1F* qqH_SM125 = refill((TH1F*)input->Get(TString::Format("%s/qqH_SM125",directory)),"qqH_SM125");InitSignal(qqH_SM125);qqH_SM125->Scale(SIGNAL_SCALE);
  TH1F* VH_SM125 = refill((TH1F*)input->Get(TString::Format("%s/VH_SM125",directory)),"VH_SM125");InitSignal(VH_SM125);VH_SM125->Scale(SIGNAL_SCALE);
//  TH1F* WHToBB_SM125 = refill((TH1F*)input->Get(TString::Format("%s/WHToBB_SM125",directory)),"WHToBB_SM125");InitSignal(WHToBB_SM125);WHToBB_SM125->Scale(SIGNAL_SCALE);
  TH1F* ZHToBB_SM125 = refill((TH1F*)input->Get(TString::Format("%s/ZHToBB_SM125",directory)),"ZHToBB_SM125");InitSignal(ZHToBB_SM125);ZHToBB_SM125->Scale(SIGNAL_SCALE);
*/
/*
  TH1F* ggAToZhToLLTauTau = refill((TH1F*)input2->Get(TString::Format("%s/ggAToZhToLLTauTau$MH",directory)),"ggAToZhToLLTauTau"); InitSignal(ggAToZhToLLTauTau);
  TH1F* ggAToZhToLLBB = refill((TH1F*)input2->Get(TString::Format("%s/ggAToZhToLLBB$MH",directory)),"ggAToZhToLLBB"); InitSignal(ggAToZhToLLBB);
  TH1F* bbH    = refill((TH1F*)input2->Get(TString::Format("%s/bbH$MH" , directory)), "bbH"); InitSignal(bbH);
*/
#endif
#ifdef ASIMOV
  TH1F* data   = refill((TH1F*)input->Get(TString::Format("%s/data_obs_asimov", directory)), "data", true);
#else
  TH1F* data   = refill((TH1F*)input->Get(TString::Format("%s/data_obs", directory)), "data",true);
#endif
  InitHist(data, "#bf{m_{H} [GeV]}", "#bf{dN/dm_{H} [1/GeV]}"); InitData(data);

  TH1F* ref=(TH1F*)Fakes->Clone("ref");
//  ref->Add(EWK1 );
  ref->Add(EWK2 );
//ref->Add(EWK3 );
  ref->Add(EWK );
  ref->Add(ttbar);
  ref->Add(Ztt  );


  double unscaled[8];
  unscaled[0] = Fakes->Integral();
  unscaled[1] = EWK  ->Integral();
//  unscaled[1]+= EWK1 ->Integral();
  unscaled[1]+= EWK2 ->Integral();
//unscaled[1]+= EWK3 ->Integral();
  unscaled[2] = ttbar->Integral();
  unscaled[3] = Ztt  ->Integral();
#ifdef MSSM
  unscaled[4] = ggHTohhTo2Tau2B  ->Integral();
  /*unscaled[5] = ggH_SM125->Integral();
  unscaled[5]+= qqH_SM125->Integral();
  unscaled[5]+= VH_SM125->Integral();
*/
/*
  unscaled[5] = ggAToZhToLLTauTau->Integral();
  unscaled[6] = ggAToZhToLLBB->Integral();
  unscaled[7] = bbH  ->Integral();
*/
#endif

  if(scaled){

/*    Fakes = refill(shape_histos(Fakes, datacard, "QCD"), "QCD");
    EWK1  = refill(shape_histos(EWK1, datacard, "W"), "W");
//.........这里部分代码省略.........

void shapehist() {
    bool makeplot  = true;

    double bins[] = {250., 266, 285., 305., 330., 360., 400., 450., 500., 550., 1000.};
    int nbins = 10;

    TH1F* hisdt = new TH1F("data_obs"                 , "", nbins, bins);
    TH1F* hisdm = new TH1F("DM"                       , "", nbins, bins);
    TH1F* hiszj = new TH1F("ZJets"                    , "", nbins, bins);
    TH1F* hisdi = new TH1F("Dibosons"                 , "", nbins, bins);
    TH1F* histt = new TH1F("ttbar"                    , "", nbins, bins);
    TH1F* hisst = new TH1F("singletop"                , "", nbins, bins);
    TH1F* hisqc = new TH1F("QCD"                      , "", nbins, bins);
    TH1F* hiswj = new TH1F("WJets"                    , "", nbins, bins);
    TH1F* hiszn = new TH1F("Znunu"                    , "", nbins, bins);
    TH1F* hist1 = new TH1F("hist1"                    , "", nbins, bins);
    TH1F* hist2 = new TH1F("hist2"                    , "", nbins, bins);
    TH1F* hiswu = new TH1F("WJets_Stat_WJETSUp"       , "", nbins, bins);
    TH1F* hiswd = new TH1F("WJets_Stat_WJETSDown"     , "", nbins, bins);
    TH1F* hiszu = new TH1F("Znunu_Stat_ZNUNUUp"       , "", nbins, bins);
    TH1F* hiszd = new TH1F("Znunu_Stat_ZNUNUDown"     , "", nbins, bins);

    hisdt->Sumw2();
    hisdm->Sumw2();
    hiszn->Sumw2();
    hiswj->Sumw2();
    hiszj->Sumw2();
    hisdi->Sumw2();
    histt->Sumw2();
    hisst->Sumw2();
    hisqc->Sumw2();
    hiszu->Sumw2();
    hiszd->Sumw2();
    hiswu->Sumw2();
    hiswd->Sumw2();
    hist1->Sumw2();
    hist2->Sumw2();

    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/met/ztree.root",    "tree/tree", 1, 1.0 , hiszn);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/bkgz/ztree.root",   "tree/tree", 1, 19.7, hist1);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/bkgnowz/ztree.root","tree/tree", 1, 19.7, hist1);

    hiszn->Add(hist1, -1.0);

    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/met/wtree.root",    "tree/tree", 2, 1.0 , hiswj);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/bkgw/wtree.root",   "tree/tree", 2, 19.7, hist2);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/bkgnowz/wtree.root","tree/tree", 2, 19.7, hist2);

    hiswj->Add(hist2, -1.0);

    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/dmVM100/reducedtree.root"  ,  "tree/tree", 0, 19.7     , hisdm);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/met/reducedtree.root"      ,  "tree/tree", 0, 1.0      , hisdt);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/zjets/reducedtree.root"    ,  "tree/tree", 0, 19.7     , hiszj);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/dibosons/reducedtree.root" ,  "tree/tree", 0, 19.7     , hisdi);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/ttbar/reducedtree.root"    ,  "tree/tree", 0, 19.7     , histt);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/singletop/reducedtree.root",  "tree/tree", 0, 19.7     , hisst);
    fillhist("/home/avartak/CMS/MonoX/CMSSW_5_3_12/src/MonoXAnalysis/AnalysisStep/trees/qcd/reducedtree.root"      ,  "tree/tree", 0, 19.7*1.3 , hisqc);

    for (int i = 1; i <= nbins; i++) {
        hiszu->SetBinContent(i, hiszn->GetBinContent(i) + hiszn->GetBinError(i));
        hiszd->SetBinContent(i, hiszn->GetBinContent(i) - hiszn->GetBinError(i));
        hiswu->SetBinContent(i, hiswj->GetBinContent(i) + hiswj->GetBinError(i));
        hiswd->SetBinContent(i, hiswj->GetBinContent(i) - hiswj->GetBinError(i));
    }

    for (int i = 1; i <= nbins; i++) {
        if (hiszj->GetBinContent(i) == 0) hiszj->SetBinContent(i, 0.001);
        if (hisdi->GetBinContent(i) == 0) hisdi->SetBinContent(i, 0.001);
        if (histt->GetBinContent(i) == 0) histt->SetBinContent(i, 0.001);
        if (hisst->GetBinContent(i) == 0) hisst->SetBinContent(i, 0.001);
        if (hisqc->GetBinContent(i) == 0) hisqc->SetBinContent(i, 0.001);

        hiszj->SetBinError(i, 0.5*hiszj->GetBinContent(i));
        hisdi->SetBinError(i, 0.5*hisdi->GetBinContent(i));
        histt->SetBinError(i, 0.5*histt->GetBinContent(i));
        hisst->SetBinError(i, 0.5*hisst->GetBinContent(i));
        hisqc->SetBinError(i, 0.5*hisqc->GetBinContent(i));
    }


    std::cout << "data_obs  : " << hisdt->Integral() << std::endl;
    std::cout << "DM        : " << hisdm->Integral() << std::endl;
    std::cout << "Znunu     : " << hiszn->Integral() << std::endl;
    std::cout << "WJets     : " << hiswj->Integral() << std::endl;
    std::cout << "ZJets     : " << hiszj->Integral() << std::endl;
    std::cout << "Dibosons  : " << hisdi->Integral() << std::endl;
    std::cout << "ttbar     : " << histt->Integral() << std::endl;
    std::cout << "singletop : " << hisst->Integral() << std::endl;
    std::cout << "QCD       : " << hisqc->Integral() << std::endl;

    if (makeplot) {
    hiszj->Add(hisdi);
    hiszj->Add(histt);
    hiszj->Add(hisst);
    hiszj->Add(hisqc);

    hiswj->Add(hiszj);
    hiszn->Add(hiswj);
    hiszu->Add(hiswj);
    hiszd->Add(hiswj);
//.........这里部分代码省略.........

void macro_MakeQcdClosureTest()
{
  // parameters //////////////////////////////////////////////////////////////
  //TFile input("./emuSpec_19619pb-1.root", "open");
  TFile input("test_19619pb-1.root", "open");
  input.cd();

  TParameter<float> *lumi = (TParameter<float> *)input.Get("lumi");

  const int nBins = 75;
  const bool usePu = 1;
  const bool useWeight = 1;
  const int qcdEst = 1; // estimation method of QCD contribution. none(0), from SS spectrum(1), from fake rate(2)

  int eRegion = 2; // electron region EB(0), EE(1), EB+EE(2)

  bool plotSign[3];
  plotSign[0] = 1;  // all
  plotSign[1] = 1;  // SS same sign
  plotSign[2] = 1;  // OS opposite sign

  bool plotType[2];
  plotType[0] = 1;  // emu spectrum
  plotType[1] = 1;  // cumulative emu spectrum

  const bool plotPull = 0; // plot (data-bkg)/bkg
  const bool plotPullBelowSpec = 0; // plot (data-bkg)/bkg below spectrum
  const bool logPlotX = 0;
  const bool logPlotY = 1;
  const bool prelim = 1;
  const bool groupedPlot = 0;
  const bool overflowBin = 1;

  float xRangeMin = 60.;
  float xRangeMax = 1200.;
  //float xRangeMin = 0.;
  //float xRangeMax = 1500.;
  float yRangeMin[6] = {0.002, 0.002, 0.002, 0.4, 0.4, 0.4};
  float yRangeMax[6] = {30, 10, 30, 3000, 1000, 3000};
  float yRangeMinRatio[3] = {-0.7, -0.7, -0.7};
  float yRangeMaxRatio[3] = {0.7, 0.7, 0.7};
  float fitMin = xRangeMin;
  float fitMax = 1100.; // set to highest bin with a data point
  float xRangeMinRatio = fitMin;
  float xRangeMaxRatio = fitMax;

  // output file formats
  const bool savePull = 0;
  const bool saveSpec = 0;
  const bool saveCumSpec = 0;
  const bool saveAsPdf = 0;
  const bool saveAsPng = 1;
  const bool saveAsRoot = 0;
  const char *fileNameExtra = "";
  //const char *fileNameExtra = "madgraphTTbar_";
  const char *plotDir = "./plottemp/";

  // plot style
  int wjetColour=  TColor::GetColor("#ffd324");
  int jetBkgColour = TColor::GetColor("#ffff66"); 

  int font = 42; //62
  ////////////////////////////////////////////////////////////////////////////

  // systematic errors
  float systErrLumi = ((TParameter<float> *)input.Get("systErrLumi"))->GetVal();
  systErrLumi = 0.; // since we normalize to the Z peak
  float systErrEff = ((TParameter<float> *)input.Get("systErrEff"))->GetVal(); // muon err & ele err
  THashList *systErrMCs = (THashList *)input.Get("systErrMCs");
  vector<float> systErrMC;
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcTtbar"))->GetVal());  // NNLO ttbar
  //systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcTtbar700to1000"))->GetVal());  // NLO ttbar700to1000
  //systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcTtbar1000up"))->GetVal());  // NLO ttbar1000up
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcDyTauTau"))->GetVal()); //z->tt
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcWW"))->GetVal()); //WW
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcWZ"))->GetVal()); //WZ
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcZZ"))->GetVal()); //ZZ
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcTW"))->GetVal()); //tW
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcDyMuMu"))->GetVal()); //Z->mm
  systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcDyEE"))->GetVal()); //Z->ee
  if (qcdEst == 2) systErrMC.push_back(0.4); // qcd error
  else systErrMC.push_back(((TParameter<float> *)systErrMCs->FindObject("systErrMcWJets"))->GetVal());  //WJets

  // to keep the histogram when the file is closed
  TH1::AddDirectory(kFALSE);
  TH1::SetDefaultSumw2(kTRUE);

  TString histoSign[3] = {"", "SS_", "OS_"};
  TString xAxisTitle[3] = {"m(e#mu)", "m(e^{#pm}#mu^{#pm})", "m(e^{#pm}#mu^{#mp})"};
  TString nameSuffix[2] = {"", "cumul"};
  TString titleSuffix[2] = {"", " - Cumulative"};

  vector<TH1F *> emuMass_wjets;
  vector<TH1F *> emuMass_qcd;
  vector<TH1F *> emuMass_qcdFromFake;

  // define the binning
  vector<float> binning;
  if (logPlotX) {
    //for (float bin = 0.; bin < 100.; bin += 5.)
//.........这里部分代码省略.........

void DrawQCDClosure(TString VAR,TString XTITLE)
{
  gROOT->ForceStyle();
  TString FileName[7] = {"QCD_HT200to300","QCD_HT300to500","QCD_HT500to700","QCD_HT700to1000","QCD_HT1000to1500","QCD_HT1500to2000","QCD_HT2000toInf"};
  float XSEC[7]       = {1.74e+6,3.67e+5,2.94e+4,6.524e+03,1.064e+03,121.5,2.542e+01};
  TFile *inf[7];
  TH1F  *h[7],*h1[7];

  TCanvas *can = new TCanvas("can_QCDClosure_"+VAR,"can_QCDClosure_"+VAR,900,600);
  can->cd(1);
  can->SetBottomMargin(0.3);
  can->SetRightMargin(0.15);

  for(int i=0;i<7;i++) {
    inf[i] = TFile::Open("Histo_"+FileName[i]+".root");
    TH1F *hTriggerPass = (TH1F*)inf[i]->Get("hadtopL/TriggerPass");
    h[i]   = (TH1F*)inf[i]->Get("hadtopL/h_"+VAR);
    h1[i]   = (TH1F*)inf[i]->Get("hadtop/h_"+VAR);
    h[i]->Sumw2();
    h1[i]->Sumw2();
    h[i]->Rebin(5);
    h1[i]->Rebin(5);
    h[i]->Scale(XSEC[i]/hTriggerPass->GetBinContent(1));
    h1[i]->Scale(XSEC[i]/hTriggerPass->GetBinContent(1));
    cout<<hTriggerPass->GetBinContent(1)<<endl;
  }

  TH1F *hQCD  = (TH1F*)h[0]->Clone("hQCD");
  TH1F *hQCD1 = (TH1F*)h1[0]->Clone("hQCD1");
  for(int i=0;i<7;i++) {
    hQCD->Add(h[i]);
    hQCD1->Add(h1[i]);
  } 
  hQCD->SetFillColor(kGray);

  hQCD->Scale(1./hQCD->Integral());
  hQCD1->Scale(1./hQCD1->Integral());
  hQCD->GetXaxis()->SetLabelSize(0.0);
  double max = 1.1*TMath::Max(hQCD->GetBinContent(hQCD->GetMaximumBin()),hQCD1->GetBinContent(hQCD1->GetMaximumBin()));
  hQCD->SetMinimum(1e-5);
  hQCD->SetMaximum(max);
  hQCD->Draw("hist");
  hQCD1->Draw("sameE");
  gPad->RedrawAxis();

  TLegend *leg = new TLegend(0.86,0.65,0.99,0.9);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);
  leg->AddEntry(hQCD,"Control","F");
  leg->AddEntry(hQCD1,"Signal","LP");
  leg->Draw();

  TH1F *hRatio = (TH1F*)hQCD1->Clone("Ratio");
  hRatio->Divide(hQCD);

  TPad* pad = new TPad("pad", "pad", 0., 0., 1., 1.);
  pad->SetTopMargin(0.7);
  pad->SetRightMargin(0.15);
  pad->SetFillColor(0);
  pad->SetFillStyle(0);
  pad->Draw();
  pad->cd(0);
  gPad->SetGridy();
  hRatio->GetXaxis()->SetTitle(XTITLE);
  hRatio->GetYaxis()->SetNdivisions(505);
  hRatio->GetYaxis()->SetRangeUser(0,2);
  hRatio->GetYaxis()->SetLabelSize(0.04);
  hRatio->Draw();
}

//.........这里部分代码省略.........
  wjets->Write();
  ewk->SetName("VV");
  ewk->SetTitle("VV");
  ewk->Write();
  vbfh->SetName("qqH");
  vbfh->SetTitle("qqH");
  vbfh->Write();
  ggh->SetName("ggH");
  ggh->SetTitle("ggH");
  ggh->Write();
  assoh->SetName("assoH");
  assoh->SetTitle("assoH");
  assoh->Write();
  smhh->SetName("ggHH");
  smhh->SetTitle("ggHH");
  smhh->Write();
  hh_0->SetName("lam0");
  hh_0->SetTitle("lam0");
  hh_0->Write();
  hh_m5->SetName("lamm5");
  hh_m5->SetTitle("lamm5");
  hh_m5->Write();
  hh_m1->SetName("lamm1");
  hh_m1->SetTitle("lamm1");
  hh_m1->Write();
  hh_p5->SetName("lamp5");
  hh_p5->SetTitle("lamp5");
  hh_p5->Write();
  hh_vbf->SetName("qqHH");
  hh_vbf->SetTitle("qqHH");
  hh_vbf->Write();
  outDC->Close();
  //stack some  histtograms together
  vbfh->Add(ggh);
  vbfh->Add(assoh);
  cout << "Single H    "  << vbfh->IntegralAndError(0,vbfh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  wjets->Add(ewk); 
  cout << "Total  Electroweak    "  << wjets->IntegralAndError(0,wjets->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  //-----------------------------------------------------------------------
  smhh->Scale(sigscale);
  //Draw the histograms
  TCanvas *canv = MakeCanvas("canv", "histograms", 800, 600);
  canv->cd();
  wjets->Add(ttbar);
  Ztt->Add(wjets);
  vbfh->Add(Ztt);
  //Error band stat
  TH1F* errorBand = (TH1F*)vbfh->Clone("errorBand");
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(13);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  //  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
  //     if(errorBand->GetBinContent(idx)>0){
  //       std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
  //       break;
  //     }
  //}
  vbfh->SetMaximum(1.0*std::max(maximum(vbfh, 0), maximum(smhh, 0)));
  //blind(data,75,150);
  //data->Draw("e");
  vbfh->Draw("hist");
  Ztt->Draw("histsame");
  wjets->Draw("histsame");

void drawVBFSpec(TTree **iTree,TH1F **iH,TH1F **iHSS,TH1F **iHMT,TH1F **iHLIS,TH1F **iHTIS,TH1F **iHNMT,TH1F **iHMTSS,TH1F **iHTemp,int iN,std::string iVar,std::string iCut,bool is2012,TFile * iFile = 0,std::string iDirName="",bool *iUseScale,TH1F** iHHigh,TH1F** iHLow) { 
  std::string lVBF = "*(njetingap  == 0 && mjj >  500 && jdeta > 3.5)";// && nbtag == 0)"; 
  for(int i0 = 0; i0 < iN; i0++) {
      //iHNMT [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 > -10)*(iso_1 < 0.1 && iso_2  >  0.795)               "+lVBF," nm_{T}");
      //iHMTSS[i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 > -10)*(iso_1 < 0.3 && iso_2  > -0.50 )               "+lVBF," nm_{T} SS");
    iH    [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)*(iso_1 < 0.1 &&  "+fIso+" )                 "+lVBF," Main");
    iHSS  [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 < 20)*(iso_1 < 0.3 &&  "+fIso+"  && iso_1 > 0.1 )"+lVBF," Same Sign");
    iHMT  [i0]   = draw(iVar,iTree[i0],i0,iCut+"*              (mtMVA_1 > 70 && mtMVA_1 < 120)*(iso_1 < 0.1 &&  "+fIso+" && m_sv < 145)"+lVBF," m_{T} > 70 GeV");
    iHNMT [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 >  0)                 *(iso_1 < 0.3 &&  "+fIso+" )"+lVBF," No m");
    iHMTSS[i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 > 70 && mtMVA_1 < 120)*(iso_1 < 0.1 &&  "+fIso+" )"+lVBF," XNomSS");
    iHTemp[i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*"+fFreeWeights[i0]                                                   ," XFree");
    iHLIS [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 < 20)*(iso_1 < 0.3 &&  "+fIso+" & iso_1 > 0.1)"  ," X2 Jet Loose Iso SS Cut");
    iHTIS [i0]   = draw(iVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 < 20)*(iso_1 < 0.1 &&  "+fIso+" )"               ," Y2 Jet Tight Iso SS Cut");
    cout << " +===> " << iHLIS[i0]->Integral() << " -- " << iHTIS[i0]->Integral() << endl;
    if(iFile == 0    ) continue;
    iHHigh[i0]   = 0;
    iHLow[i0]    = 0;
    if(!iUseScale[i0]) continue;
    iHHigh[i0]   = draw(iVar+"High",iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)"," Main High");
    iHLow [i0]   = draw(iVar+"Low" ,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)"," Main Low");
  }
  //Comput MT Scale Factor
  TH1F *lMTMC  = (TH1F*) iHMT  [0]->Clone("mTTmp");  clear(lMTMC);
  TH1F *lMTMCS = (TH1F*) iHMTSS[0]->Clone("mTTmpS"); clear(lMTMCS);
  for(int i0 = 0; i0 < iN-1; i0++) {
    if(i0 == fWId || i0 == fQCDId) continue;
    lMTMC ->Add(iHMT[i0]);
    lMTMCS->Add(iHMTSS[i0]);
  }
  double lDataInt   = iHMT[iN-1]  ->Integral(-1,1000) - lMTMC ->Integral(-1,1000); 
  double lDataSSInt = iHMTSS[iN-1]->Integral(-1,1000) - lMTMCS->Integral(-1,1000); 
  double lWSF       = float(iHMT  [fWId]->Integral(-1,1000))/lDataInt;
  double lWSFSS     = float(iHMTSS[fWId]->Integral(-1,1000))/lDataSSInt;
  cout << "===> W Boson Scale Factor : " << lWSF << " -W- " << lWSFSS  << " -Data- "<< lDataInt << " - " << iHMT[iN-1]->Integral(-1,1000) << "  MC- " << lMTMC->Integral(-1,1000) << endl;
  if(lWSF   == 0) lWSF   = 1.;
  if(lWSFSS == 0) lWSFSS = 1.;
  iH    [fWId]->Scale(1./lWSF);
  iHSS  [fWId]->Scale(1./lWSF);
  iHMT  [fWId]->Scale(1./lWSF);
  iHMTSS[fWId]->Scale(1./lWSF);
  //!!!!! Do we need to do this
  //iHLIS [fWId]->Scale(1./lWSF);
  //iHTIS [fWId]->Scale(1./lWSF);
  
  //Compute QCD Shape
  TH1F *lSS   = (TH1F*) iHSS  [0]->Clone("SSTmp0"); clear(lSS);
  TH1F *lSSL  = (TH1F*) iHLIS [0]->Clone("SSTmp1"); clear(lSSL);
  TH1F *lSST  = (TH1F*) iHTIS [0]->Clone("SSTmp2"); clear(lSST);
  TH1F *lSSMT = (TH1F*) iHMTSS[0]->Clone("SSTmp3"); clear(lSSMT);
  for(int i0 = 0; i0 < iN-1; i0++) {
    if(i0 == fQCDId) continue;
    lSS  ->Add(iHSS  [i0]);
    lSSL ->Add(iHLIS [i0]);
    lSST ->Add(iHTIS [i0]);
    lSSMT->Add(iHMTSS[i0]);
  }
  iHSS  [fQCDId]->Add(lSS   ,-1); 
  iHMTSS[fQCDId]->Add(lSSMT ,-1); 
  iHLIS [fQCDId]->Add(lSSL  ,-1); 
  iHTIS [fQCDId]->Add(lSST  ,-1); 
  clear(iHMT[fQCDId]);
  for(int i0 = 0; i0 < iHSS  [fQCDId]->GetNbinsX()+1; i0++) if(iHSS  [fQCDId]->GetBinContent(i0) < 0) iHSS  [fQCDId]->SetBinContent(i0,0);
  for(int i0 = 0; i0 < iHMTSS[fQCDId]->GetNbinsX()+1; i0++) if(iHMTSS[fQCDId]->GetBinContent(i0) < 0) iHMTSS[fQCDId]->SetBinContent(i0,0);
  for(int i0 = 0; i0 < iHLIS [fQCDId]->GetNbinsX()+1; i0++) if(iHLIS [fQCDId]->GetBinContent(i0) < 0) iHLIS [fQCDId]->SetBinContent(i0,0);
  for(int i0 = 0; i0 < iHTIS [fQCDId]->GetNbinsX()+1; i0++) if(iHTIS [fQCDId]->GetBinContent(i0) < 0) iHTIS [fQCDId]->SetBinContent(i0,0);
  iHTIS[fQCDId]->Scale(1.06);
  double lTightLooseRatio = iHTIS[fQCDId]->Integral()/iHLIS[fQCDId]->Integral(); 
  cout << "===> Tight/Loose QCD Ratio " << lTightLooseRatio << endl;
  iHSS   [fQCDId]->Scale(lTightLooseRatio);
  iHTIS  [fQCDId]->Scale(iHSS[fQCDId]->Integral()/iHTIS[fQCDId]->Integral());
  iH     [fQCDId] = iHSS  [fQCDId];
  //Generic Template Definition for the loose cuts
  for(int i0 = 0; i0 < iN-1; i0++) {
    if(i0 == fQCDId || i0 == 0) continue; //Keep full selection for Z and QCD
    double pInt = iH[i0]->Integral();
    iH[i0] = iHTemp[i0];
    iH[i0]->Scale(pInt/iH[i0]->Integral());
  }
  //Make Data card
  if(iFile !=0) makeDataCard(iFile,iH,iHHigh,iHLow,iN,iDirName,fString);
  //Additional Crap
  //iHNMT [fWId]->Scale(1./lWSF);
  //iHMTSS [fQCDId]->Scale(lTightLooseRatio);
  //iHNMT  [fQCDId] = iHMTSS[fQCDId];

  //Blind
  //for(int i0 = 0; i0 < iH[iN-1]->GetNbinsX()+1; i0++) if(iH[iN-1]->GetXaxis()->GetBinCenter(i0) > 60 && iH[iN-1]->GetXaxis()->GetBinCenter(i0) < 130) iH[iN-1]->SetBinContent(i0,0);
  //Draw the plot
  draw(iH    ,iN,iVar+"VBFA",iVar);
  draw(iHSS  ,iN,iVar+"VBFB",iVar);
  draw(iHMT  ,iN,iVar+"VBFC",iVar);
  //draw(iHLIS ,iN,iVar+"VBFD",iVar);
  //draw(iHTIS ,iN,iVar+"VBFE",iVar);
  //draw(iHMTSS,iN,iVar+"VBFF",iVar);
  //draw(iHNMT ,iN,iVar+"VBFG",iVar);
}

void drawSpec(TTree **iTree,TH1F **iH,TH1F **iHSS,TH1F **iHMT,TH1F **iHNMT,TH1F **iHSSMT,TH1F **iHMTS,int iN,std::string iVar,std::string iCut,std::string iName,TFile * iFile = 0,std::string iDirName="",bool *iUseScale,TH1F** iHHigh,TH1F** iHLow) { 
  std::string lVar = iVar;
  TH1F *lQCDShape = 0;
  for(int i0 = 0; i0 < iN; i0++) {
    lVar = iVar;
    iH   [i0]   = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)"," Main");
    iHSS [i0]   = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 < 20)"," Same Sign");
    iHMT [i0]   = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 > 70)"," OS m_{T} > 70 GeV");
    iHMTS[i0]   = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 > 70)"," SS m_{T} > 70 GeV");
    cout << "====> " << fString[i0] << " -- " << iH[i0]->Integral() << " -- " << fWId << " -- " << fQCDId << endl;
    if(i0 == fQCDId && iName.find("A") == std::string::npos) lQCDShape =  draw(lVar,iTree[i0],iN,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 < 20)"," Same Sign Shape");
    
    //Additional Useless plots
    //iHNMT[i0]   = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 >  0)"," No m_{T} Cut");
    //iHSSMT[i0]  = draw(lVar,iTree[i0],i0,iCut+"*(q_1*q_2 > 0)*(mtMVA_1 >  0)"," Same Sign");
    
    //Dealing with datacards
    if(iFile == 0    ) continue;
    iHHigh[i0]   = 0;
    iHLow[i0]    = 0;
    if(!iUseScale[i0]) continue;
    iHHigh[i0]   = draw(iVar+"High",iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)"," Main High");
    iHLow [i0]   = draw(iVar+"Low" ,iTree[i0],i0,iCut+"*(q_1*q_2 < 0)*(mtMVA_1 < 20)"," Main Low");
  }

  //Comput MT Scale Factor => separately for SS and OS
  TH1F *lMTMC  = (TH1F*) iHMT [0]->Clone("mTTmp");  clear(lMTMC);
  TH1F *lMTMCS = (TH1F*) iHMTS[0]->Clone("mTTmpS"); clear(lMTMCS);
  for(int i0 = 0; i0 < iN-1; i0++) {
    if(i0 == fWId || i0 == fQCDId) continue;
    lMTMC ->Add(iHMT [i0]);
    lMTMCS->Add(iHMTS[i0]);
  }
  double lDataInt   = iHMT [iN-1]->Integral(-1,1000) - lMTMC ->Integral(-1,1000); 
  double lWSF       = float(iHMT[fWId]->Integral(-1,1000))/lDataInt;
  double lDataSSInt = iHMTS[iN-1]->Integral(-1,1000) - lMTMCS->Integral(-1,1000); 
  double lWSFSS     = float(iHMTS[fWId]->Integral(-1,1000))/lDataSSInt;
  cout << "===> W Boson    Scale Factor : " << lWSF   << " -W- " << iHMT[fWId]->Integral(-1,1000) << " -- " << iH[fWId]->Integral(-1,1000)   << endl;
  cout << "===> W Boson SS Scale Factor : " << lWSFSS << " -W- " << iHMT[fWId]->Integral(-1,1000) << endl; 
  if(lWSF   == 0) lWSF   = 1.;
  if(lWSFSS == 0) lWSFSS = 1.;
  iH    [fWId]->Scale(1./lWSF);
  iHSS  [fWId]->Scale(1./lWSFSS);
  iHMT  [fWId]->Scale(1./lWSF);//No ss
  iHMTS [fWId]->Scale(1./lWSFSS);

  //Compute QCD Yield
  TH1F *lSS   = (TH1F*) iHSS[0]->Clone("SSTmp");   clear(lSS);
  for(int i0 = 0; i0 < iN-1; i0++) {if(i0 == fQCDId) continue;    lSS  ->Add(iHSS  [i0]);  }
  iHSS  [fQCDId]->Add(lSS  ,-1); 
  //Scale the Shape
  if(lQCDShape != 0) lQCDShape->Scale(iHSS[fQCDId]->Integral()/lQCDShape->Integral());
  if(lQCDShape != 0) iHSS  [fQCDId] = lQCDShape;
  
  for(int i0 = 0; i0 < iHSS  [fQCDId]->GetNbinsX()+1; i0++) if(iHSS  [fQCDId]->GetBinContent(i0) < 0) iHSS  [fQCDId]->SetBinContent(i0,0);
  TH1F *lXSS = (TH1F*) iHSS[fQCDId]->Clone("SSF");  lXSS->Scale(1.06);
  iH   [fQCDId] = lXSS;//iHSS[fQCDId];
  
  //Dealing with Addtional Plots
  //iHNMT [fWId]->Scale(1./lWSF); //No ss
  //iHSSMT[fWId]->Scale(1./lWSFSS);
  //TH1F *lSSMT = (TH1F*) iHSS[0]->Clone("SSTmpmT"); clear(lSSMT);
  //for(int i0 = 0; i0 < iN-1; i0++) {if(i0 == fQCDId) continue;    lSSMT->Add(iHSSMT[i0]);  }
  //iHSSMT[fQCDId]->Add(lSSMT,-1); 
  clear(iHMT [fQCDId]);
  //clear(iHNMT[fQCDId]);
  //for(int i0 = 0; i0 < iHSSMT[fQCDId]->GetNbinsX()+1; i0++) if(iHSSMT[fQCDId]->GetBinContent(i0) < 0) iHSSMT[fQCDId]->SetBinContent(i0,0);
  //iHSSMT[fQCDId]->Scale(1.06);
  //iHNMT[fQCDId] = iHSSMT[fQCDId];
  
  //Make DataCards ==> If Asked
  if(iFile !=0) makeDataCard(iFile,iH,iHHigh,iHLow,iN,iDirName,fString);
  
  //Blind
  //for(int i0 = 0; i0 < iH[iN-1]->GetNbinsX()+1; i0++) if(iH[iN-1]->GetXaxis()->GetBinCenter(i0) > 100 && iH[iN-1]->GetXaxis()->GetBinCenter(i0) < 150) iH[iN-1]->SetBinContent(i0,0);
 
  //Draw the plot
  draw(iH   ,iN,iVar+"A"+iName,iVar);
  draw(iHSS ,iN,iVar+"B"+iName,iVar);
  draw(iHMT ,iN,iVar+"C"+iName,iVar);
  //draw(iHNMT,iN,iVar+"D"+iName,iVar);
  //draw(iHMTS,iN,iVar+"E"+iName,iVar);
}

void 
HTT_EM_X(bool scaled=true, bool log=true, float min=0.1, float max=-1., TString datacard="htt_em_0_7TeV", string inputfile="root/$HISTFILE", const char* directory="emu_$CATEGORY")
{
  // define common canvas, axes pad styles
  SetStyle(); gStyle->SetLineStyleString(11,"20 10");

  // determine category tag
  const char* category = ""; const char* category_extra = ""; const char* category_extra2 = "";
  if(std::string(directory) == std::string("emu_0jet_low"             )){ category = "e#mu, 0 jet";          }    
  if(std::string(directory) == std::string("emu_0jet_low"             )){ category_extra = "p_{T}(lep1) low";          }    
  if(std::string(directory) == std::string("emu_0jet_high"            )){ category = "e#mu, 0 jet";          }    
  if(std::string(directory) == std::string("emu_0jet_high"            )){ category_extra = "p_{T}(lep1) high";         }    
  if(std::string(directory) == std::string("emu_1jet_low"          )){ category = "e#mu, 1 jet";          }    
  if(std::string(directory) == std::string("emu_1jet_low"          )){ category_extra = "p_{T}(lep1) low";       }    
  if(std::string(directory) == std::string("emu_1jet_high"          )){ category = "e#mu, 1 jet";          }    
  if(std::string(directory) == std::string("emu_1jet_high"          )){ category_extra = "p_{T}(lep1) high";       }    
  if(std::string(directory) == std::string("emu_vbf_loose"            )){ category = "e#mu, 2 jet";          }    
  if(std::string(directory) == std::string("emu_vbf_loose"            )){ category_extra = "VBF, loose";              }    
  if(std::string(directory) == std::string("emu_vbf_tight"            )){ category = "e#mu, 2 jet";          }    
  if(std::string(directory) == std::string("emu_vbf_tight"            )){ category_extra = "VBF, tight";              }    
  if(std::string(directory) == std::string("emu_nobtag"               )){ category = "e#mu";          }    
  if(std::string(directory) == std::string("emu_nobtag"               )){ category_extra = "No B-Tag";                        }    
  if(std::string(directory) == std::string("emu_btag"                 )){ category = "e#mu";          }    
  if(std::string(directory) == std::string("emu_btag"                 )){ category_extra = "B-Tag";                           }

  const char* dataset;
  if(std::string(inputfile).find("7TeV")!=std::string::npos){dataset = "CMS Preliminary,  H#rightarrow#tau#tau, 4.9 fb^{-1} at 7 TeV";}
  if(std::string(inputfile).find("8TeV")!=std::string::npos){dataset = "CMS Preliminary,  H#rightarrow#tau#tau, 19.7 fb^{-1} at 8 TeV";}
  
  TFile* input = new TFile(inputfile.c_str());
#ifdef MSSM
  TFile* input2 = new TFile((inputfile+"_$MA_$TANB").c_str());
#endif
  TH1F* Fakes  = refill((TH1F*)input->Get(TString::Format("%s/Fakes"     , directory)), "Fakes"  ); InitHist(Fakes, "", "", kMagenta-10, 1001);
  TH1F* EWK    = refill((TH1F*)input->Get(TString::Format("%s/EWK"       , directory)), "EWK"    ); InitHist(EWK  , "", "", kRed    + 2, 1001);
  TH1F* ttbar  = refill((TH1F*)input->Get(TString::Format("%s/ttbar"     , directory)), "ttbar"  ); InitHist(ttbar, "", "", kBlue   - 8, 1001);
  TH1F* Ztt    = refill((TH1F*)input->Get(TString::Format("%s/Ztt"       , directory)), "Ztt"    ); InitHist(Ztt  , "", "", kOrange - 4, 1001);
#ifdef MSSM
  TH1F* ggH  = refill((TH1F*)input2->Get(TString::Format("%s/ggH$MA"     , directory)), "ggH"    ); InitSignal(ggH    ); ggH    ->Scale($TANB);
  TH1F* bbH  = refill((TH1F*)input2->Get(TString::Format("%s/bbH$MA"     , directory)), "bbH"    ); InitSignal(bbH    ); bbH    ->Scale($TANB);
#else
#ifdef HWW_BG
  TH1F* ggH_hww= refill((TH1F*)input->Get(TString::Format("%s/ggH_hww125", directory)), "ggH_hww"); InitHist(ggH_hww  , "", "", kGreen + 2, 1001); ggH_hww->Scale(SIGNAL_SCALE);
  TH1F* qqH_hww= refill((TH1F*)input->Get(TString::Format("%s/qqH_hww125", directory)), "qqH_hww"); InitHist(qqH_hww  , "", "", kGreen + 2, 1001); qqH_hww->Scale(SIGNAL_SCALE);
#else
 TH1F* ggH_hww= refill((TH1F*)input->Get(TString::Format("%s/ggH_hww125", directory)), "ggH_hww"); InitSignal(ggH_hww); ggH_hww->Scale(SIGNAL_SCALE);
  TH1F* qqH_hww= refill((TH1F*)input->Get(TString::Format("%s/qqH_hww125", directory)), "qqH_hww"); InitSignal(qqH_hww); qqH_hww->Scale(SIGNAL_SCALE);
#endif
#ifndef DROP_SIGNAL
  TH1F* ggH    = refill((TH1F*)input->Get(TString::Format("%s/ggH125"    , directory)), "ggH"    ); InitSignal(ggH    ); ggH    ->Scale(SIGNAL_SCALE);
  TH1F* qqH    = refill((TH1F*)input->Get(TString::Format("%s/qqH125"    , directory)), "qqH"    ); InitSignal(qqH    ); qqH    ->Scale(SIGNAL_SCALE);
  TH1F* VH     = refill((TH1F*)input->Get(TString::Format("%s/VH125"     , directory)), "VH"     ); InitSignal(VH     ); VH     ->Scale(SIGNAL_SCALE);
#ifndef HWW_BG
  ggH->Add(ggH_hww);
  qqH->Add(qqH_hww);
#endif
#endif
#endif
#ifdef ASIMOV
  TH1F* data   = refill((TH1F*)input->Get(TString::Format("%s/data_obs_asimov", directory)), "data", true);
#else
  TH1F* data   = refill((TH1F*)input->Get(TString::Format("%s/data_obs", directory)), "data", true);
#endif
  InitHist(data, "#bf{m_{#tau#tau} [GeV]}", "#bf{dN/dm_{#tau#tau} [1/GeV]}"); InitData(data);

  TH1F* ref=(TH1F*)Fakes->Clone("ref");
  ref->Add(EWK  );
  ref->Add(ttbar);
  ref->Add(Ztt  );
#ifdef HWW_BG
  ref->Add(qqH_hww);
  ref->Add(ggH_hww);
#endif

#ifdef HWW_BG
  const int numhistos = 9;
  const int offset = 2;
#else
  const int numhistos = 7;
  const int offset = 0;
#endif
  double unscaled[numhistos];
  unscaled[0] = Fakes->Integral();
  unscaled[1] = EWK  ->Integral();
  unscaled[2] = ttbar->Integral();
  unscaled[3] = Ztt  ->Integral();
#ifdef MSSM
  unscaled[4] = ggH  ->Integral();
  unscaled[5] = bbH  ->Integral();
  unscaled[6] = 0;
#else
#ifdef HWW_BG
  unscaled[4] = ggH_hww  ->Integral();
  unscaled[5] = qqH_hww  ->Integral();
#endif
#ifndef DROP_SIGNAL
  unscaled[4+offset] = ggH  ->Integral();
  unscaled[5+offset] = qqH  ->Integral();
  unscaled[6+offset] = VH   ->Integral();
#endif
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  cout<<"GluGlu = "<<GluGlu<<" +/- "<<DeltaF4<<endl;
  double QCD = hbkg5->Integral();
  cout<<"QCD = "<<QCD<<" +/- "<<DeltaF5<<endl;
 double VHiggs = hbkg6->Integral();
  cout<<"VH = "<<VHiggs<<" +/- "<<DeltaF6<<endl;


 
  h1->SetLineColor(kRed-4);
  h2->SetLineColor(kRed+2);  
  h3->SetLineColor(kPink+2);
  h4->SetLineColor(kPink+4);
  h5->SetLineColor(kMagenta+2); //only for 25ns samples
  h1->SetLineWidth(2);
  h2->SetLineWidth(2);  
  h3->SetLineWidth(2);
  h4->SetLineWidth(2);
  h5->SetLineWidth(2); //only for 25ns samples
  
  
  THStack *hs=new THStack("hs","");
  hbkg6->SetFillColor(kGreen+3);  
  hbkg4->SetFillColor(kGreen+4);
  hbkg1->SetFillColor(kBlue);
  hbkg2->SetFillColor(kBlue+2);
  hbkg3->SetFillColor(kMagenta+4);
  hbkg5->SetFillColor(kMagenta+2);
  
  hbkg1->SetLineColor(kBlack);
  hbkg2->SetLineColor(kBlack);
  hbkg3->SetLineColor(kBlack);
  hbkg4->SetLineColor(kBlack);
  hbkg5->SetLineColor(kBlack);
  hbkg6->SetLineColor(kBlack);
  
  hs->Add(hbkg6);
  hs->Add(hbkg4);
  hs->Add(hbkg1);
  hs->Add(hbkg2);
  hs->Add(hbkg3);
  hs->Add(hbkg5);
  
  TH1F *hsum = (TH1F*)hbkg1->Clone("hsum"); 
  hsum->Add(hbkg2);
  hsum->Add(hbkg3);
  hsum->Add(hbkg4);
  hsum->Add(hbkg5);
  hsum->Add(hbkg6);
  hs->SetMinimum(0.001);
  hs->SetMaximum(5000);
  hs->SetTitle("");
  hs->Draw("HIST"); 
  hsum->SetMarkerStyle(1);
  hsum->SetFillColor(kGray+3);
  hsum->SetFillStyle(3001);
  hsum->Draw("same e2");
  h2->Draw("same hist");
  h3->Draw("same hist"); 
  h4->Draw("same hist");
  h1->Draw("same hist");
  h5->Draw("same hist"); //only for 25ns samples
  if(DataFlag==1){
    hdata->Draw("same E1");
  }

  hs->GetXaxis()->SetTitle("MET [GeV]");
  hs->GetYaxis()->SetTitleOffset(1.2);
  hs->GetYaxis()->SetTitle("Events/15 GeV");  
  gPad->Modified();



  /*leg->AddEntry(h1,"m_{#chi} = 1 GeV","l");
  leg->AddEntry(h2,"m_{#chi} = 10 GeV","l");
  leg->AddEntry(h3,"m_{#chi} = 100 GeV","l");      
  leg->AddEntry(h4,"m_{#chi} = 1000 GeV","l");*/
  leg->AddEntry(h1,"m_{Z'} = 600 GeV","l");
  leg->AddEntry(hbkg1,"#gamma #gamma","f");         
  leg->AddEntry(h2,"m_{Z'} = 800 GeV","l");
  leg->AddEntry(hbkg2,"Drell Yann","f");       
  leg->AddEntry(h3,"m_{Z'} = 1000 GeV","l");      
  leg->AddEntry(hbkg3,"#gamma + Jets","f");      
  leg->AddEntry(h4,"m_{Z'} = 1200 GeV","l");     
  leg->AddEntry(hbkg4,"ggH","f");      
  leg->AddEntry(h5,"m_{Z'} = 1400 GeV","l"); //only for 25ns samples    
  leg->AddEntry(hbkg5,"QCD","f");      
  leg->AddEntry(hbkg6,"VH","f");  
  leg->Draw("same");
  
  gStyle->SetOptStat(0);  
 

  /*    if(jets<3){
    canvas->SaveAs(Form("met_cuts_M%d_%djets.pdf",sigMass,(jets-1))); 
    canvas->SaveAs(Form("met_cuts_M%d_%djets.png",sigMass,(jets-1))); 
  }else{
    canvas->SaveAs(Form("met_cuts_M%d_noJetsCut.pdf",sigMass)); 
    canvas->SaveAs(Form("met_cuts_M%d_noJetsCut.png",sigMass)); 
    }*/
}