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

void ComputeWeights(){
  
  //*** mc pileup file
  //  TFile *fmc = TFile::Open("./PU_DYJetsToLL_Summer12_DR53X-PU_S10.root"); // 2012
  TFile *fmc = TFile::Open("./PU_DYJetsToLL_Summer12_DR53X-PU_RD1.root"); // 2012 RD ReReco
  //TFile *fmc = TFile::Open("./PU_DYJetsToLL_Summer11dr53X-PU_S13.root"); // 2011 53X ReReco
  TH1F *hmc  = (TH1F*)fmc->Get("hmc");
  TH1F *hmctrue  = (TH1F*)fmc->Get("hmctrue");
  cout << "Number of bins in MC histogram: " << hmc->GetNbinsX()<< endl;
  int nbinsmc = hmc->GetNbinsX();

  //*** data file -- observed PU
  //  TFile *fda   = TFile::Open("../pileup/pileup_2012D_minBiasXsec69400_corr_observed.root");
  //TFile *fda   = TFile::Open("../pileup/pileup_2012ABCD_22Jan2013ReReco_corr_observed.root");
  //TFile *fda   = TFile::Open("../pileup/pileup_2012ABCD_22Jan2013ReReco_pixelcorr_observed.root");
  //TFile *fda   = TFile::Open("../pileup/pileup_2011_minBiasXsec68000_pixelcorr_observed.root");

  // globe pileup files
  //  TFile *fda   = TFile::Open("/afs/cern.ch/work/m/malberti/HGG/LegacyPaper/CMSSW_6_1_1/src/h2gglobe/AnalysisScripts/aux/AB.json.69400.observed.pileup.root");
  //  TFile *fda   = TFile::Open("/afs/cern.ch/work/m/malberti/HGG/LegacyPaper/CMSSW_6_1_1/src/h2gglobe/AnalysisScripts/aux/C.json.69400.observed.pileup.root");
  TFile *fda   = TFile::Open("/afs/cern.ch/work/m/malberti/HGG/LegacyPaper/CMSSW_6_1_1/src/h2gglobe/AnalysisScripts/aux/D.json.69400.observed.pileup.root");

  TH1F  *hdata = (TH1F*)fda->Get("pileup");
  cout << "Number of bins in DATA histogram: " << hdata->GetNbinsX()<< endl;
  int nbinsdata = hdata->GetNbinsX();

  if (nbinsdata!=nbinsmc){
    TH1D  *htemp = (TH1D*)fda->Get("pileup");
    hdata = new TH1F("hdata","hdata",60,0,60);
    for (int ibin = 1; ibin < 61; ibin++){
      cout << hdata->GetBinContent(ibin) <<endl;
      hdata->SetBinContent(ibin,htemp->GetBinContent(ibin));
      cout << htemp->GetBinContent(ibin) <<endl;
    }
  }

  //*** compute weights
  TH1F *hweights = (TH1F*)hdata->Clone("hweights");
  hweights->Divide(hdata,hmc,1./hdata->GetSumOfWeights(),1./hmc->GetSumOfWeights());

  //  TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer12_DR53X-PU_RD1_minBiasXsec69400_corr_observed_2012ABCD.root","create");
  //  TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer12_DR53X-PU_RD1_minBiasXsec69400_pixelcorr_observed_2012ABCD.root","create");
  //  TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer11dr53X-PU_S13_minBiasXsec68000_pixelcorr_observed_2011.root","create");
  //  TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer11dr53X-PU_RD1_minBiasXsec69400_AB.root","create");
  //TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer11dr53X-PU_RD1_minBiasXsec69400_C.root","create");
  TFile *fout = new TFile("./PUweights_DYJetsToLL_Summer11dr53X-PU_RD1_minBiasXsec69400_D.root","create");
  hweights->Write("hweights");
  hdata->Write("hdata");
  hmc->Write("hmc");
  fout->Close();


//   //*** data file -- true PU
//   TFile *fdatrue = TFile::Open("../pileup/pileup_190456-208686_minBiasXsec69400_corr_observed.root");
//   TH1F *hdatatrue = (TH1F*)fdatrue->Get("pileup");
//   cout << hdatatrue->GetNbinsX()<< endl;

//   //*** compute weights
//   TH1F *hweightstrue = (TH1F*)hdatatrue->Clone("hweightstrue");
//   hweightstrue->Divide(hdatatrue,hmctrue,1./hdatatrue->GetSumOfWeights(),1./hmctrue->GetSumOfWeights());
 
//   TFile *fout2 = new TFile("./PUweights_DYJetsToLL_Summer12_DR53X-PU_S10_minBiasXsec69400_corr_true_Run2012ABCD.root","create");
//   hweightstrue->Write("hweights");
//   hdatatrue->Write("hdata");
//   hmctrue->Write("hmc");
//   fout2->Close();
}

TCanvas *drawOneVariable(TTree *signalTree, TTree *backgroundTree1, TTree *backgroundTree2,TTree *backgroundTree3,
			 TCut signalCuts, TCut backgroundCuts,
			 TString var,
			 int nbins, double xlow, double xhigh,
			 TString sigLegend, TString bg1Legend, TString bg2Legend,TString bg3Legend,
			 TString comment)
{

  TString cname = "c_";
  cname += var;
  TCanvas *c1 = new TCanvas(cname,cname,10,10,600,600);
  c1->cd();

  TH1F *hsig = new TH1F(TString("hsig_")+var,"",nbins, xlow, xhigh);
  TH1F *hbg1 = new TH1F(TString("hbg1_")+var,"",nbins, xlow, xhigh);
  TH1F *hbg2 = new TH1F(TString("hbg2_")+var,"",nbins, xlow, xhigh);
  TH1F *hbg3 = new TH1F(TString("hbg3_")+var,"",nbins, xlow, xhigh);

  TString projectCommandSig = var+TString(">>hsig_")+var;
  TString projectCommandBg1 = var+TString(">>hbg1_")+var;
  TString projectCommandBg2 = var+TString(">>hbg2_")+var;
  TString projectCommandBg3 = var+TString(">>hbg3_")+var;
  
  if( !useSmallEventCount ){
    signalTree->Draw(projectCommandSig, "genWeight"*signalCuts);
  }else{
    printf("DEBUG MODE: using small event count\n");
    signalTree->Draw(projectCommandSig, "genWeight"*signalCuts, "", 100000);
  }
  TGaxis::SetMaxDigits(3);
  hsig->GetXaxis()->SetTitle(var);
  hsig->SetDirectory(0);

  if(backgroundTree1 != 0){
    if( !useSmallEventCount ){
      backgroundTree1->Draw(projectCommandBg1, "genWeight"*backgroundCuts);
    }else{
      printf("DEBUG MODE: using small event count\n");
      backgroundTree1->Draw(projectCommandBg1, "genWeight"*backgroundCuts, "", 100000);
    }
    hbg1->Scale(hsig->GetSumOfWeights() / hbg1->GetSumOfWeights());
    hbg1->SetDirectory(0);
  } else {
    delete hbg1;
    hbg1 = 0;
  }

  if(backgroundTree2 != 0){
    if( !useSmallEventCount ){
      backgroundTree2->Draw(projectCommandBg2, "genWeight"*backgroundCuts);
    }else{
      printf("DEBUG MODE: using small event count\n");
      backgroundTree2->Draw(projectCommandBg2, "genWeight"*backgroundCuts, "", 100000);
    }
    hbg2->Scale(hsig->GetSumOfWeights() / hbg2->GetSumOfWeights());
    hbg2->SetDirectory(0);
  } else {
    delete hbg2;
    hbg2 = 0;
  }


  if(backgroundTree3 != 0){
    if( !useSmallEventCount ){
      backgroundTree3->Draw(projectCommandBg3, "genWeight"*backgroundCuts);
    }else{
      printf("DEBUG MODE: using small event count\n");
      backgroundTree3->Draw(projectCommandBg3, "genWeight"*backgroundCuts, "", 100000);
    }
    hbg3->Scale(hsig->GetSumOfWeights() / hbg3->GetSumOfWeights());
    hbg3->SetDirectory(0);
  } else {
    delete hbg3;
    hbg3 = 0;
  }


  setHistogramAttributes(hsig, hbg1, hbg2, hbg3);

  c1->Clear();

  hsig->Draw("hist");
  if( hbg1 ){

    hbg1->Draw("same");
  }
  if( hbg2 ){

    hbg2->Draw("same");
  }
  if( hbg3 ){

    //hbg3->Draw("same");
  }
  TLegend *leg = new TLegend(0.55, 0.65, 0.95, 0.80); // 0.6 0.9
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->AddEntry(hsig, sigLegend, "lf");
  leg->AddEntry(hbg1, bg1Legend, "lf");
  leg->AddEntry(hbg2, bg2Legend, "lf");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TString TStrName(rootFileName);
  std::cout <<TStrName <<std::endl;  

  TFile * file = new TFile(era+"/"+TStrName+TString(".root"),"update");
  file->cd("");

  TH1D * inputEventsH = new TH1D("inputEventsH","",1,-0.5,0.5);
  TH1D * histWeightsH = new TH1D("histWeightsH","",1,-0.5,0.5);


  TH1D * metAll  = new TH1D("metAll","",400,0,4000);

  TH1D * hDiJetmet = new TH1D("hDiJetmet","",400,0,4000);
  TH1D * hDiJetmass = new TH1D("hDiJetmass","",400,0,4000);
  TH1D * hDiJet1mass = new TH1D("hDiJet1mass","",400,0,4000);
  TH1D * hDiJet2mass = new TH1D("hDiJet2mass","",400,0,4000);
  TH1D * hHT_ = new TH1D("hHT_","",400,0,4000);
  TH1D * PUweightsOfficialH = new TH1D("PUweightsOfficialH","PU weights w/ official reweighting",1000, 0, -1);

 // PILE UP REWEIGHTING - OPTIONS

  if (applyPUreweighting_vertices and applyPUreweighting_official) 
	{std::cout<<"ERROR: Choose only ONE PU reweighting method (vertices or official, not both!) " <<std::endl; exit(-1);}

  // reweighting with vertices

  // reading vertex weights
  TFile * fileDataNVert = new TFile(TString(cmsswBase)+"/src/"+dataBaseDir+"/"+vertDataFileName);
  TFile * fileMcNVert   = new TFile(TString(cmsswBase)+"/src/"+dataBaseDir+"/"+vertMcFileName);

  TH1F * vertexDataH = (TH1F*)fileDataNVert->Get(TString(vertHistName));
  TH1F * vertexMcH   = (TH1F*)fileMcNVert->Get(TString(vertHistName));

  float normVertexData = vertexDataH->GetSumOfWeights();
  float normVertexMc   = vertexMcH->GetSumOfWeights();

  vertexDataH->Scale(1/normVertexData);
  vertexMcH->Scale(1/normVertexMc);


 // reweighting official recipe 
  	// initialize pile up object
  PileUp * PUofficial = new PileUp();
  

  
  if (applyPUreweighting_official) {
    TFile * filePUdistribution_data = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/PileUpDistrib/Data_Pileup_2015D_Nov17.root","read"); 
    TFile * filePUdistribution_MC = new TFile (TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/PileUpDistrib/MC_Spring15_PU25_Startup.root", "read"); 
    TH1D * PU_data = (TH1D *)filePUdistribution_data->Get("pileup");
    TH1D * PU_mc = (TH1D *)filePUdistribution_MC->Get("pileup");
    PUofficial->set_h_data(PU_data); 
    PUofficial->set_h_MC(PU_mc);
  }

  

  int nFiles = 0;
  int nEvents = 0;
  int selEventsAllMuons = 0;
  int selEventsIdMuons = 0;
  int selEventsIsoMuons = 0;



  int nTotalFiles = 0;

  TLegend *leg = new TLegend(0.5,0.7,0.85, 0.90);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);


  //leg->AddEntry(&h1, "E near 1st peak", "l");
  //leg->AddEntry(&h2, "E near 2nd peak", "l");

  noise->Draw("ereco>>h1","wgt*((eshower==0)&&(ncore==0))");
  noise->Draw("ereco>>h2","wgt*((eshower==0)&&(ncore==1))");
  noise->Draw("ereco>>h3","wgt*((eshower==0)&&(ncore==2))");
  noise->Draw("ereco>>h4","wgt*((eshower==0)&&(ncore==3))");
  noise->Draw("ereco>>h5","wgt*((eshower==0)&&(ncore==4))");
  noise->Draw("ereco>>h6","wgt*((eshower==0)&&(ncore==5))");
  
  double s1 = h1->GetSumOfWeights();
  double s2 = h2->GetSumOfWeights();
  double s3 = h3->GetSumOfWeights();
  double s4 = h4->GetSumOfWeights();
  double s5 = h5->GetSumOfWeights();
  double s6 = h6->GetSumOfWeights();

  
  h1->Scale(100.0/s1);
  h2->Scale(100.0/s2);
  h3->Scale(100.0/s3);
  h4->Scale(100.0/s4);
  h5->Scale(100.0/s5);
  h6->Scale(100.0/s6);

  h1->SetMaximum(60.0);

//.........这里部分代码省略.........
	}
	
	// In case the full selection is applied:
	//       if( !(isTag1 && ele2_passID) && !(isTag2 && ele1_passID) ) continue;
	if( !(isTag1 && isIDProbePass2) && !(isTag2 && isIDProbePass1) ) continue;
	//       if( !(isTag1) && !(isTag2) ) continue;
	// Fill histogram
	hMass->Fill(dielectron->mass);
	
      } // end loop over dielectron candidates
    } // end loop over events
  
    delete infile;
    infile=0;
    eventTree=0;
    
    delete gen;
    delete info;
    delete dielectronArr;
  } // end loop over files

  //
  // Efficiency analysis
  //
  
//   printf("Number of regular candidates:      %15.0f\n", hMass->GetSumOfWeights());
  printf("Total events in ntuple                                       %15d\n",eventsInNtuple);
  printf("    events after event level trigger cut                     %15d\n",eventsAfterTrigger);
  printf("\nTotal candidates (no cuts)                                   %15d\n",totalCand);
  printf("        candidates in 60-120 mass window                     %15d\n",totalCandInMassWindow);
  printf("        candidates witheta 0-1.4442, 1.566-2.5               %15d\n",totalCandInEtaAcceptance);
  printf("        candidates, both electrons above 10 GeV              %15d\n",totalCandEtAbove10GeV);
  printf("        candidates matched to GEN level (if MC)              %15d\n",totalCandMatchedToGen);
  printf("\nNumber of probes, total                                      %15.0f\n", hMassTotal->GetSumOfWeights());
  printf("Number of probes, passed                                     %15.0f\n", hMassPass->GetSumOfWeights());
  printf("Number of probes, failed                                     %15.0f\n", hMassFail->GetSumOfWeights());


  // Human-readbale text file to store measured efficiencies
  TString reslog = tagAndProbeDir+TString("/efficiency_TnP_")+label+TString(".txt");
  ofstream effOutput;
  effOutput.open(reslog);
  // Print into the results file the header.
  effOutput << "Efficiency calculation method: " << calcMethodString.Data() << endl;
  effOutput << "Efficiency type to measure: " << effTypeString.Data() << endl;
  effOutput << "SC ET binning: " << etBinningString.Data() << endl;
  effOutput << "SC eta binning: " << etaBinningString.Data() << endl;
  effOutput << "Sample: " << sampleTypeString.Data() << endl;
  effOutput << "Files processed: " << endl;
  for(UInt_t i=0; i<ntupleFileNames.size(); i++)
    effOutput << "   " << ntupleFileNames[i].Data() << endl;
  
  // ROOT file to store measured efficiencies in ROOT format
  TString resroot = tagAndProbeDir+TString("/efficiency_TnP_")+label+TString(".root");
  TFile *resultsRootFile = new TFile(resroot,"recreate");

  // Fit log 
  TString fitlogname = tagAndProbeDir+TString("/efficiency_TnP_")+label+TString("_fitlog.dat");
  ofstream fitLog;
  fitLog.open(fitlogname);

  //
  //  Find efficiency
  //
  bool useTemplates = false;
  if(sample == DATA && effType == ID && (calcMethod == COUNTnFIT || FITnFIT) )

void HwwNoteFigs(TString filename = "../histo_tmva_new-ntuples-1_160_0.root",
                 int mass = 170, int ReBin = 4, int ysel = 0)
{
	TFile* infile = new TFile(filename.Data(),"READ");
	setTDRStyle(0);
	
	//BDTD dists
	TCanvas* c1 = new TCanvas("c1","c1",0,-200,500,500);
	
        TH1F* histoBg =  histo4->Clone();
        histoBg->Add(histo3);
        histoBg->Add(histo2);
        histoBg->Add(histo1);

    	TH1F* histo_4 =  histo4->Clone();
    	histo3->Add(histo_4);
    	TH1F* histo_3 =  histo3->Clone();
    	histo2->Add(histo_3);
    	TH1F* histo_2 =  histo2->Clone();
    	histo1->Add(histo_2);

        histo1->Rebin(ReBin);
        histo1->SetFillColor(kBlue);
        histo1->SetFillStyle(1001);
        histo1->SetLineStyle(0);
        histo1->SetLineWidth(0);

        histo2->Rebin(ReBin);
        histo2->SetFillColor(kMagenta);
        histo2->SetFillStyle(1001);
        histo2->SetLineStyle(0);
        histo2->SetLineWidth(0);

        histo3->Rebin(ReBin);
        histo3->SetFillColor(kGreen);
        histo3->SetFillStyle(1001);
        histo3->SetLineStyle(0);
        histo3->SetLineWidth(0);
        
        histo4->Rebin(ReBin);
        histo4->SetFillColor(kCyan);
        histo4->SetFillStyle(1001);
        histo4->SetLineStyle(0);
        histo4->SetLineWidth(0);

        char YTitle[300];
        sprintf(YTitle,"events / bin");
        char XTitle[300];
        sprintf(XTitle,"BDT Output");
    	histo0->SetYTitle(YTitle);
    	histo1->SetYTitle(YTitle);
    	histo2->SetYTitle(YTitle);
    	histo3->SetYTitle(YTitle);
    	histo4->SetYTitle(YTitle);

    	histo0->SetXTitle(XTitle);
    	histo1->SetXTitle(XTitle);
    	histo2->SetXTitle(XTitle);
    	histo3->SetXTitle(XTitle);
    	histo4->SetXTitle(XTitle);
    	histo1->SetTitleSize(0.05, "X");
    	histo1->GetXaxis()->SetTitleFont(62);
    	histo1->GetXaxis()->SetLabelFont(61);
    	histo1->GetYaxis()->SetLabelFont(61); 
    	histo1->GetYaxis()->SetTitleOffset(1.3);
    	histo1->SetLabelSize(0.04, "Y");
    	histo1->SetLabelSize(0.04, "X");

	int min = histoBg->FindBin(-1.5);
	int max = histoBg->FindBin(1.5);

	histoBg->GetXaxis()->SetRange(min,max);
	
	histoBg->SetMarkerStyle(20);
	histoBg->SetMarkerSize(1.0);
	histoBg->GetYaxis()->SetTitleOffset(1.40);
	
	histo0->SetMarkerStyle(21);
	histo0->SetMarkerSize(1.0);
	
	histoBg->Rebin(ReBin);
	histo0->Rebin(ReBin);

	histoBg->SetLineColor(4);
	histo0->SetLineColor(1);

	scaleHist(histoBg);
	scaleHist(histo0);
	cout << "bg events: " << histoBg->GetSumOfWeights() << endl;
	cout << "si events: " << histo0->GetSumOfWeights() << endl;

	histoBg->SetYTitle("Events");
    	histo1->SetMinimum(0.01);
	if(ysel == 0)  {
	  histo1->Draw("hist");
	}
	else          {
	  histo0->Draw("E");
	  histo1->Draw("hist,same");
    	}
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  hDPhiPtllJetWWLevelLL    ->Fill(dphilljet,    totalW*fakeW);
	  //}

	  if ( passIDISO1 && passIDISO2 ) { //&&  isSemi == 1) { //&& pfmet > 20 )  {
	    
	    hWTriggerA ->Fill(1,      totalW);
	    hPtLepton1WWLevelA      ->Fill(pt1,       totalW);
	    hPtLepton2WWLevelA      ->Fill(pt2,       totalW);
	    hPtDiLeptonWWLevelA     ->Fill(ptll,      totalW);
	    hMinvWWLevelA           ->Fill(mll,       totalW);
	    hMtWWLevelA             ->Fill(mth,       totalW);
	    hNJets30WWLevelA        ->Fill(njet,      totalW);
	    hpfMetWWLevelA          ->Fill(pfType1Met,     totalW);
	    //hppfMetWWLevelA         ->Fill(ppfmet,    totalW);
	    //hchMetWWLevelA          ->Fill(chmet,     totalW);
	    //hpchMetWWLevelA         ->Fill(pchmet,    totalW);
	    hpminMetWWLevelA        ->Fill(mpmet,     totalW);
	    hDeltaRLeptonsWWLevelA  ->Fill(drll,      totalW);
	    hDeltaPhiLeptonsWWLevelA->Fill(dphill,    totalW);
	    hDPhiPtllJetWWLevelA    ->Fill(dphilljet, totalW);
	    
	  }
	
	  if ( (passIDISO1 && !passIDISO2) || (!passIDISO1 && passIDISO2) ) {
	    
	    hWTriggerB ->Fill(1,      totalW);
	    hpfMetWWLevelBNoMET          ->Fill(pfType1Met,     totalW); 
	    
	    hPtLepton1WWLevelB      ->Fill(pt1,       totalW);
	    hPtLepton2WWLevelB      ->Fill(pt2,       totalW);
	    hPtDiLeptonWWLevelB     ->Fill(ptll,      totalW);
	    hMinvWWLevelB           ->Fill(mll,       totalW);
	    hMtWWLevelB             ->Fill(mth,       totalW);
	    hNJets30WWLevelB        ->Fill(njet,      totalW);
	    hpfMetWWLevelB          ->Fill(pfType1Met,     totalW);
	    //hppfMetWWLevelB         ->Fill(ppfmet,    totalW);
	    //hchMetWWLevelB          ->Fill(chmet,     totalW);
	    //hpchMetWWLevelB         ->Fill(pchmet,    totalW);
	    hpminMetWWLevelB        ->Fill(mpmet,     totalW);
	    hDeltaRLeptonsWWLevelB  ->Fill(drll,      totalW);
	    hDeltaPhiLeptonsWWLevelB->Fill(dphill,    totalW);
	    hDPhiPtllJetWWLevelB    ->Fill(dphilljet, totalW);
	  }	
	
	  if ( !passIDISO1 && !passIDISO2)  { // && pfmet > 20 ) { 
	    
	    hWTriggerC ->Fill(1,      totalW);
	    hPtLepton1WWLevelC      ->Fill(pt1,       totalW);
	    hPtLepton2WWLevelC      ->Fill(pt2,       totalW);
	    hPtDiLeptonWWLevelC     ->Fill(ptll,      totalW);
	    hMinvWWLevelC           ->Fill(mll,       totalW);
	    hMtWWLevelC             ->Fill(mth,       totalW);
	    hNJets30WWLevelC        ->Fill(njet,      totalW);
	    hpfMetWWLevelC          ->Fill(pfType1Met,     totalW);
	    //hppfMetWWLevelC         ->Fill(ppfmet,    totalW);
	    //hchMetWWLevelC          ->Fill(chmet,     totalW);
	    //hpchMetWWLevelC         ->Fill(pchmet,    totalW);
	    hpminMetWWLevelC        ->Fill(mpmet,     totalW);
	    hDeltaRLeptonsWWLevelC  ->Fill(drll,      totalW);
	    hDeltaPhiLeptonsWWLevelC->Fill(dphill,    totalW);
	    hDPhiPtllJetWWLevelC    ->Fill(dphilljet, totalW);	    
	  
	  }
		
      }
      
    }
    
  }


	  
  verbose = true;

  // Print
  //----------------------------------------------------------------------------
  if (verbose) {

    float norm = 1933235;
  

  //// OUTPUT 

  cout << "LL:  "  << hWTriggerLL->GetSumOfWeights() << endl;
  cout << "TT:  "  << hWTriggerA->GetSumOfWeights() << endl;
  cout << "TnoT:  "  << hWTriggerB->GetSumOfWeights() << endl;
  cout << "noTnoT:  "  << hWTriggerC->GetSumOfWeights() << endl;
  }
 
   

  // Save the histograms
  //----------------------------------------------------------------------------
  output->cd();
  output->Write("", TObject::kOverwrite);
  output->Close();
  
  

}

/*************************************************************************************
 * getYields: Goes through each Counts histogram for every process and subprocess
 *            (i.e., EE & Drell-Yan, E & WJets, etc.) and properly adds the yields
 *            for data and MC. Reports ratios of Data:MC as well.
 *  input:  the main() arguments array
 *  output: writes to stdout the LaTeX-formatted table of yields (pipe the output to 
 *          save)
 ***********************************/
void getYields(const char* argv[]) { 
  //open mlbwidth output file
  TFile *f = new TFile(argv[2]);
  f->cd();

  //very inefficient, but for now it works
  //loop over all procs, leps and figure out the event counts
  for(int i=0; i<lepsSize; i++) {
    char a[128];
    sprintf(a, "mlbwa_%s_Count", leps[i]);
    TDirectory *tDir = (TDirectory*) f->Get(a);
    TList *alok = tDir->GetListOfKeys();

    for(int j=0; j<procsSize; j++) {
      for(int k=0; alok->At(k)->GetName() != alok->Last()->GetName(); k++) {
        if(TString(alok->At(k)->GetName()).Contains(TRegexp(procs[j]))) {
          char b[128];
          sprintf(b, "mlbwa_%s_Count/%s", leps[i], alok->At(k)->GetName());

          TH1F *h = (TH1F*) f->Get(b);
          eCounts[i][j] += h->GetSumOfWeights();
          eErrors[i][j] =  sqrt(pow(eErrors[i][j],2) + pow(h->GetBinError(2),2));

          delete h;
        }
      }
    }

    char d[128];
    sprintf(d, "mlbwa_%s_Count/%s", leps[i], alok->Last()->GetName());
    if(d == "") { exit(EXIT_FAILURE); }

    TH1F *tth = (TH1F*) f->Get(d);
    eCounts[i][procsSize] = tth->GetEntries();
    double integral = tth->GetSumOfWeights();
    eErrors[i][procsSize] = tth->GetBinError(2)*eCounts[i][procsSize]/integral;
    delete tth;
  }

  //Get a string to tell us how many columns we want (size leps + 1)
  char cols[] = "c";
  for(int i=0; i<lepsSize; i++) {
    strcat(cols, "c");
  }

  //print out LaTeX:
  //formatting
  cout<<"\\documentclass[12pt,a4paper,titlepage]{article}"<<endl;
  cout<<"\\usepackage[utf8]{inputenc}"<<endl;
  cout<<"\\usepackage{amsmath}"<<endl;
  cout<<"\\usepackage{amsfonts}"<<endl;
  cout<<"\\usepackage{amssymb}"<<endl;
  cout<<"\\usepackage{hyperref}\n\n"<<endl;
  cout<<"\\usepackage[margin=0.0025in]{geometry}"<<endl;
  cout<<"\\begin{document}"<<endl;
  cout<<"\\begin{tabular}{l|"<<cols<<"} \\\\"<<endl;

  //Print out the table header
  cout<<"Sample & ";
  for(int i=0; i<lepsLaTeXSize; i++) {
    cout<<lepsLaTeX[i]<<" & ";
  }
  cout<<"Sum \\\\ \\hline\\hline"<<endl;

  //Printing out the yields for each process
  for(int i=0; i<procsSize; i++) {
      cout<<yieldLaTeX[i]<<" & ";
    for(int j=0; j<lepsSize; j++) {
      cout<<GetLatex(j,i)<<" & ";
    }
    cout<<GetLatexSum(true,i)<<"\\\\"<<(i==procsSize-1 ? "\\hline" : "")<<endl;
  }
 
  //Print out the total MC yield
  cout<<"Total MC & ";
  for(int i=0; i<lepsSize; i++) {
    cout<<GetLatexSum(false,i)<<" & ";
  }
  cout<<GetLatexSum(false,-1)<<"\\\\"<<endl;

  //Print out the data yield
  cout<<"Data & ";
  for(int i=0; i<lepsSize; i++) {
    cout<<GetLatex(i,lepsSize+1)<<" & ";
  }
  cout<<GetLatexSum(true,lepsSize+1)<<"\\\\\\hline\\hline"<<endl;
  
  //Print out the data:MC ratio
  cout<<"Data/MC & ";
  for(int i=1; i<=procsSize;i++) {
    cout<<GetLatexRatio(i)<<(i==procsSize ? "" : " & ");
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  gPad->SetRightMargin(0.15);
  gStyle->SetPalette(1);
  hResponse->GetXaxis()->SetMoreLogLabels();
  hResponse->GetXaxis()->SetNoExponent();
  hResponse->GetYaxis()->SetNoExponent();
  plotResponse.Draw(e);

  // Create the plot of the inverted response matrix
  TH2F *hInvResponse = new TH2F("hInvResponse","",DYTools::nMassBins, DYTools::massBinLimits,
			     DYTools::nMassBins, DYTools::massBinLimits);
  for(int i=1; i<=DYTools::nMassBins; i++)
    for(int j=1; j<=DYTools::nMassBins; j++)
      hInvResponse->SetBinContent(i,j,DetInvertedResponse(i-1,j-1));
  TCanvas *f = MakeCanvas("invResponse","invResponse",xsize,ysize);
  canvasV.push_back(f);
  CPlot plotInvResponse("inverted response","",
			"reconstructed  m(ee) [GeV/c^{2}]",
			"generator post-FSR m(ee) [GeV/c^{2}]");
  plotInvResponse.AddHist2D(hInvResponse,"COLZ");
  f->cd();
  plotInvResponse.SetLogx();
  plotInvResponse.SetLogy();
  gPad->SetRightMargin(0.15);
  gStyle->SetPalette(1);
  hInvResponse->GetXaxis()->SetMoreLogLabels();
  hInvResponse->GetXaxis()->SetNoExponent();
  hInvResponse->GetYaxis()->SetNoExponent();
  plotInvResponse.Draw(f);

  // Create a plot of detector resolution without mass binning
  TCanvas *g = MakeCanvas("massDiff","massDiff",600,600);
  canvasV.push_back(g);
  CPlot plotMassDiff("massDiff","","reco mass - gen post-FSR mass [GeV/c^{2}]","a.u.");
  hMassDiffBB->Scale(1.0/hMassDiffBB->GetSumOfWeights());
  hMassDiffEB->Scale(1.0/hMassDiffEB->GetSumOfWeights());
  hMassDiffEE->Scale(1.0/hMassDiffEE->GetSumOfWeights());
  plotMassDiff.AddHist1D(hMassDiffBB,"EB-EB","hist",kBlack);
  plotMassDiff.AddHist1D(hMassDiffEB,"EE-EB","hist",kBlue);
  plotMassDiff.AddHist1D(hMassDiffEE,"EE-EE","hist",kRed);
  plotMassDiff.TransLegend(0.1, 0.0);
  plotMassDiff.Draw(g);



  // Create a plot of reco - gen post-FSR mass difference for several mass bins
  TCanvas *h = MakeCanvas("massDiffV","massDiffV",600,600);
  canvasV.push_back(h);
  CPlot plotMassDiffV("massDiffV","","reco mass - gen post-FSR mass [GeV/c^{2}]","a.u.");
  hMassDiffV[7]->Scale(1.0/hMassDiffV[7]->GetSumOfWeights());
  hMassDiffV[6]->Scale(1.0/hMassDiffV[6]->GetSumOfWeights());
  hMassDiffV[5]->Scale(1.0/hMassDiffV[5]->GetSumOfWeights());
  hMassDiffV[4]->Scale(1.0/hMassDiffV[4]->GetSumOfWeights());
  plotMassDiffV.AddHist1D(hMassDiffV[4],"50-60 GeV/c^{2}","hist",kBlack);
  plotMassDiffV.AddHist1D(hMassDiffV[5],"60-76 GeV/c^{2}","hist",kBlue);
  plotMassDiffV.AddHist1D(hMassDiffV[6],"76-86 GeV/c^{2}","hist",kRed);
  plotMassDiffV.AddHist1D(hMassDiffV[7],"86-96 GeV/c^{2}","hist",kMagenta);
  plotMassDiffV.SetXRange(-20,50);
  plotMassDiffV.Draw(h);

  // Create graph of bin-to-bin corrections
  TGraphAsymmErrors *grCorrFactor 
    = new TGraphAsymmErrors(massBinCentral, DetCorrFactor,
			    massBinHalfWidth, massBinHalfWidth,
			    DetCorrFactorErrNeg, DetCorrFactorErrPos);
  TCanvas *c11 = MakeCanvas("corrFactor","corrFactor",800,600);
  canvasV.push_back(c11);

//.........这里部分代码省略.........
    if( (it->first).find("Ztautau")!=string::npos ) {
      h1->SetFillColor(kRed);
      leg->AddEntry(h1,"Z+jets, genuine-#tau","F");
    }
    if( (it->first).find("TTbar")!=string::npos ) {
      h1->SetFillColor(kBlue);
      leg->AddEntry(h1,"t#bar{t}+jets","F");
    }
    if( (it->first).find("SingleTop")!=string::npos ) {
      h1->SetFillColor(29);
      leg->AddEntry(h1,"single-t","F");
    }
    if( (it->first).find("Wjets")!=string::npos ) {
      h1->SetFillColor(kGreen);
      leg->AddEntry(h1,"W+jets","F");
    }
    if( (it->first).find("DiBoson")!=string::npos ) {
      h1->SetFillColor(38);
      leg->AddEntry(h1,"Di-Boson","F");
    }
    if( (it->first).find("QCD")!=string::npos ) {
      h1->SetFillColor(42);
      leg->AddEntry(h1,"QCD-multijets","F");
    }
    if((it->first).find("qqH115")!=string::npos){
      hSgn1 = (TH1F*)h1->Clone("hSgn1");
      hSgn1->Scale(magnifySgn_);
      hSgn1->SetLineWidth(2);
      h1->SetFillColor(kBlack);
      h1->SetFillStyle(3004);
      h1->SetLineColor(kBlack);
      leg->AddEntry(h1,Form("VBF H(115) X %.0f",magnifySgn_),"F");
    }
    if((it->first).find("ggH115")!=string::npos){
      hSgn2 = (TH1F*)h1->Clone("hSgn2");
      hSgn2->Scale(magnifySgn_);
      hSgn2->SetLineWidth(2);
      h1->SetFillColor(kBlack);
      h1->SetFillStyle(3005);
      h1->SetLineColor(kBlack);
      leg->AddEntry(h1,Form("GGF H(115) X %.0f",magnifySgn_),"F");
    }
    if((it->first).find("Data")!=string::npos){
      hData = (TH1F*)h1->Clone("hData");
      hData->Sumw2();
      hData->SetMarkerStyle(20);
      hData->SetMarkerSize(1.2);
      hData->SetMarkerColor(kBlack);
      hData->SetLineColor(kBlack);
      hData->SetXTitle(XTitle_+" ("+Unities_+")");
      hData->SetYTitle(Form(" Events/(%.0f %s)", hData->GetBinWidth(1), Unities_.Data() ) );
      hData->SetTitleSize(0.04,"X");
      hData->SetTitleSize(0.05,"Y");
      hData->SetTitleOffset(0.95,"Y");
      leg->AddEntry(hData,"DATA","P");
    }

    if(iter==0) hSiml=(TH1F*)h1->Clone("hSiml");
    else if((it->first).find("Data")==string::npos) hSiml->Add(h1);

    if((it->first).find("Data")==string::npos) aStack->Add(h1);

    if(VERBOSE) cout<<(it->first) << " ==> " 
		   << h1->Integral() << " +/- " 
		   << TMath::Sqrt(h1->GetEntries())*(h1->Integral()/h1->GetEntries())
		   << endl;
  }

  // all signal summed together:
  hSgn = (TH1F*)hSgn1->Clone("hSgn");
  hSgn->Add(hSgn1,hSgn2,1,1);

  if(VERBOSE) cout<< "S/sqrt(B) ==> " 
		  << hSgn->Integral()/ TMath::Sqrt(hSiml->Integral()) << " +/- " 
		  << (1./2)*TMath::Sqrt(hSiml->GetEntries())*(hSiml->GetSumOfWeights())/hSiml->Integral()*( hSgn->Integral()/ TMath::Sqrt(hSiml->Integral())  )
		  << endl;

  hData->SetXTitle(XTitle_+" ("+Unities_+")");
  hData->SetYTitle(Form(" Events/(%.1f %s)", hData->GetBinWidth(1), Unities_.Data() ) );
  hData->SetTitleSize(0.04,"X");
  hData->SetTitleSize(0.05,"Y");
  hData->SetTitleOffset(0.95,"Y");

  hData->Draw("P");
  aStack->Draw("HISTSAME");
  hData->Draw("PSAME");
  //hSgn1->Draw("HISTSAME");
  //hSgn2->Draw("HISTSAME");

  TH1F* hStack = (TH1F*)aStack->GetHistogram();
  hStack->SetXTitle(XTitle_+" ("+Unities_+")");
  hStack->SetYTitle(Form(" Events/(%.0f %s)", hStack->GetBinWidth(1), Unities_.Data() ) );
  hStack->SetTitleSize(0.04,"X");
  hStack->SetTitleSize(0.05,"Y");
  hStack->SetTitleOffset(0.95,"Y");
  //aStack->GetYaxis()->SetRangeUser(0.0,std::max(hData->GetMaximum(),hStack->GetMaximum())*1.1);
  leg->Draw();


}