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

void Wpt_ZmassCompEtaBins_Gaus(const TString Mode,//Channel - Muon or Electron
    const TString corrName,
    const TString outputDir 
    )
{
  TString plotTitle;
  TString mu_etaRange[6];
  mu_etaRange[0] = "-2.1 #geq eta < -1.4";
  mu_etaRange[1] = "-1.4 #geq eta < -0.7";
  mu_etaRange[2] = "-0.7 #geq eta < 0";
  mu_etaRange[3] = "0 #geq eta < 0.7";
  mu_etaRange[4] = "0.7 #geq eta < 1.4";
  mu_etaRange[5] = "1.4 #geq eta < 2.1";

  TH1D* makeDiffHist(TH1D* h1, TH1D* h2, const TString name);
  const TString format("png"); 
  Int_t ratioColor = kGray+2;

  TFile *fname_MC;
  TFile *fname_RD;

  gSystem->mkdir(outputDir,kTRUE);
  CPlot::sOutDir = outputDir;

  fname_MC = new TFile("Muon2012LoPU/Muon_DYToMuMu_S8.root");
  fname_RD = new TFile("Muon2012LoPU/Muon_RD_LowPU.root");

  if(Mode=="Electron")
  {
    fname_MC = new TFile("Electron2012LoPU/Ele_DYToEE_S8.root");
    fname_RD = new TFile("Electron2012LoPU/Ele_RD_LowPU.root");
  }

  CPlot *plotMllEtaBinP;
  CPlot *plotMllEtaBinM;
  CPlot *plotMllEtameanP;
  CPlot *plotMllEtameanM;
  CPlot *plotMllEtawidthP;
  CPlot *plotMllEtawidthM;

  TH1D *hMCetaBinP[ScaleBins];
  TH1D *hRDetaBinP[ScaleBins];
  TH1D *hMCetaBinM[ScaleBins];
  TH1D *hRDetaBinM[ScaleBins];

  RooDataHist *h1_MCetaBinP;
  RooDataHist *h1_RDetaBinP;
  RooDataHist *h1_MCetaBinM;
  RooDataHist *h1_RDetaBinM;

  TH1D *hMCmeanp = new TH1D("hMCmeanp","hMCmeanp",ScaleBins,-2.1,2.1);hMCmeanp->Sumw2();
  TH1D *hRDmeanp = new TH1D("hRDmeanp","hRDmeanp",ScaleBins,-2.1,2.1);hRDmeanp->Sumw2();
  TH1D *hMCmeanm = new TH1D("hMCmeanm","hMCmeanm",ScaleBins,-2.1,2.1);hMCmeanm->Sumw2();
  TH1D *hRDmeanm = new TH1D("hRDmeanm","hRDmeanm",ScaleBins,-2.1,2.1);hRDmeanm->Sumw2();

  TH1D *hMCwidthp = new TH1D("hMCwidthp","hMCwidthp",ScaleBins,-2.1,2.1);hMCwidthp->Sumw2();
  TH1D *hRDwidthp = new TH1D("hRDwidthp","hRDwidthp",ScaleBins,-2.1,2.1);hRDwidthp->Sumw2();
  TH1D *hMCwidthm = new TH1D("hMCwidthm","hMCwidthm",ScaleBins,-2.1,2.1);hMCwidthm->Sumw2();
  TH1D *hRDwidthm = new TH1D("hRDwidthm","hRDwidthm",ScaleBins,-2.1,2.1);hRDwidthm->Sumw2();

  char histName[50];
  char tmpName[50];

  TCanvas *myCan;
  myCan = MakeCanvas("myCan","myCan",800,600);

  myCan->SetPad(0,0,1.0,1.0);
  myCan->SetTopMargin(0.11);
  myCan->SetBottomMargin(0.15);
  myCan->SetLeftMargin(0.15);  
  myCan->SetRightMargin(0.05);  
  myCan->SetTickx(1);
  myCan->SetTicky(1);  

  //=============================
  //Read Zmass histograms for each pt bin
  //=============================
  for(int i(0);i<ScaleBins;i++){
    sprintf(tmpName,"h1_Zmass_muEtaP_%d",i);
    if(corrName=="CorrTotalRegion")
      sprintf(tmpName,"h1_ZmassCorr_muEtaP_%d",i);

    if(outputDir=="Wpt_ZmassPlotsEtaBins_noOverLap_Gaus")
    {
      sprintf(tmpName,"h1_Zmass_noOverLap_muEtaP_%d",i);
      if(corrName=="CorrTotalRegion")
	sprintf(tmpName,"h1_ZmassCorr_noOverLap_muEtaP_%d",i);
    }

    if(outputDir=="Wpt_ZmassPlotsEtaBins_LeadingLept_noOverLap_Gaus")
    {
      sprintf(tmpName,"h1_Zmass_LeadingLept_noOverLap_muEtaP_%d",i);
      if(corrName=="CorrTotalRegion")
	sprintf(tmpName,"h1_ZmassCorr_LeadingLept_noOverLap_muEtaP_%d",i);
    }

    if(outputDir=="Wpt_ZmassPlotsEtaBins_LeadingLept_Gaus")
    {
      sprintf(tmpName,"h1_Zmass_LeadingLept_muEtaP_%d",i);
      if(corrName=="CorrTotalRegion")
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ntReference->SetAlias("var","matchPt2/pt2");
    ntMix->SetAlias("var","matchPt2/pt2");
  }


  nt->Draw("var>>hFull",dijet&&noise&&jetID&&centHF&&vtx); 
  nt->Draw("var>>hB",side&&noise&&jetID&&centHF&&vtx);
  nt->Draw("pt1>>hNorm",lead&&noise&&jetID&&centHF&&vtx);
   
  ntMix->Draw("var>>hDataMixFull",weight*(dijet&&jetID&&centHF&&vtx));
  ntMix->Draw("var>>hDataMixB",weight*(side&&jetID&&centHF&&vtx));
  ntMix->Draw("pt1>>hNormDataMix",weight*(lead&&jetID&&centHF&&vtx));
  
  ntReference->Draw("var>>hReferenceFull",dijet&&noise&&jetID&&centHF&&vtx);
  ntReference->Draw("var>>hReferenceB",side&&noise&&jetID&&centHF&&vtx);
  ntReference->Draw("pt1>>hNormReference",lead&&noise&&jetID&&centHF&&vtx);

  hDataMix->Add(hDataMixFull);
  h->Add(hFull);
  hReference->Add(hReferenceFull);

  if(subtract){
    hDataMix->Add(hDataMixB,-1);
    h->Add(hB,-1);
    hReference->Add(hReferenceB,-1);
  }

  hB->SetFillStyle(3005);
  hB->SetFillColor(15);

  // calculate the statistical error and normalize
  h->SetLineColor(dataColor);
  h->SetMarkerColor(dataColor);
  h->Sumw2();
  if(normLead){
    h->Scale(1./hNorm->Integral());
    hB->Scale(1./hNorm->Integral());
  }else{
    hB->Scale(1./h->Integral());
    h->Scale(1./h->Integral());
  }
  h->SetMarkerStyle(20);

  if(hReference->Integral() > 0){
    hReference->Scale(1./hNormReference->Integral());
  }

  hReference->SetLineColor(kBlue);
  hReference->SetFillColor(kAzure-8);
  hReference->SetFillStyle(3005);

  if(normLead){
    hDataMixB->Scale(1./hNormDataMix->Integral());
    hDataMix->Scale(1./hNormDataMix->Integral());
  }else{
    hDataMixB->Scale(1./hDataMix->Integral());
    hDataMix->Scale(1./hDataMix->Integral());
  }
  hDataMix->SetLineColor(mixColor);
  hDataMix->SetFillColor(mixColor);
  hDataMix->SetFillStyle(3004);
  
  hDataMix->SetMarkerSize(0);
  hDataMix->SetStats(0);

  hDataMix->GetXaxis()->SetLabelSize(22);

void fitWe(const TString  outputDir="test",   // output directory
           const Double_t lumi=18.7,        // integrated luminosity (/fb)
	   const Int_t    Ecm=8,         // center-of-mass energy
	   const Int_t    doPU=1         // option for PU-reweighting
) {
  gBenchmark->Start("fitWe");

  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //==============================================================================================================   
  
  // MET histogram binning and range
  const Int_t    NBINS  = 50;
  const Double_t METMAX = 100;

  const Double_t PT_CUT  = 25;
  const Double_t ETA_CUT = 2.5;

  // file format for output plots
  const TString format("png"); 
    
  // file name with recoil correction
  TString recoilfname("../RecoilSyst/ZeeData/fits.root");
  
  // file name(s) with PU weights
  TString pufname("");
  if(doPU>0) {
    if(doPU==1) { pufname = "Utils/PileupReweighting.Summer11DYmm_To_Run2011A.root"; }
    else        { cout << "Invalid option for PU re-weighting! Aborting..." << endl; assert(0); }
  }
  
  //
  // input ntuple file names
  //
  enum { eData, eWenu, eEWK };  // data type enum
  vector<TString> fnamev;
  vector<Int_t>   typev;
  
  fnamev.push_back("../Selection/Wenu/ntuples/data_m23_select.root");   typev.push_back(eData);
  fnamev.push_back("../Selection/Wenu/ntuples/we_select.root");         typev.push_back(eWenu);
  fnamev.push_back("../Selection/Wenu/ntuples/ewk_select.root");        typev.push_back(eEWK);
  fnamev.push_back("../Selection/Wenu/ntuples/top_select.root");        typev.push_back(eEWK);


  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================  
  
  // Create output directory
  gSystem->mkdir(outputDir,kTRUE);
  CPlot::sOutDir = outputDir;  
  
  // Get pile-up weights
  TFile *pufile    = 0;
  TH1D  *puWeights = 0;
  if(doPU>0) {
    pufile = new TFile(pufname);
    assert(pufile);
    puWeights = (TH1D*)pufile->Get("puWeights");
    assert(puWeights);
  }
  
  // Access recoil corrections
  //RecoilCorrector recoilCorr(recoilfname);
  
  //
  // Declare MET histograms
  //
  TH1D *hDataMet  = new TH1D("hDataMet", "",NBINS,0,METMAX); hDataMet->Sumw2();
  TH1D *hDataMetm = new TH1D("hDataMetm","",NBINS,0,METMAX); hDataMetm->Sumw2();  
  TH1D *hDataMetp = new TH1D("hDataMetp","",NBINS,0,METMAX); hDataMetp->Sumw2();
  TH1D *hWenuMet  = new TH1D("hWenuMet", "",NBINS,0,METMAX); hWenuMet->Sumw2();
  TH1D *hWenuMetp = new TH1D("hWenuMetp","",NBINS,0,METMAX); hWenuMetp->Sumw2();
  TH1D *hWenuMetm = new TH1D("hWenuMetm","",NBINS,0,METMAX); hWenuMetm->Sumw2();
  TH1D *hEWKMet   = new TH1D("hEWKMet",  "",NBINS,0,METMAX); hEWKMet->Sumw2();
  TH1D *hEWKMetp  = new TH1D("hEWKMetp", "",NBINS,0,METMAX); hEWKMetp->Sumw2();
  TH1D *hEWKMetm  = new TH1D("hEWKMetm", "",NBINS,0,METMAX); hEWKMetm->Sumw2();

  //
  // Declare variables to read in ntuple
  //
  UInt_t  runNum, lumiSec, evtNum;
  UInt_t  npv, npu;
  Float_t genWPt, genWPhi;
  Float_t scale1fb;
  Float_t met, metPhi, sumEt, mt, u1, u2;
  Int_t   q;
  LorentzVector *lep=0;
  LorentzVector *sc=0;
    
  TFile *infile=0;
  TTree *intree=0;

  //
  // Loop over files
  //
  for(UInt_t ifile=0; ifile<fnamev.size(); ifile++) {
    
    // Read input file and get the TTrees
    cout << "Processing " << fnamev[ifile] << "..." << endl;
//.........这里部分代码省略.........

void limit() {
//This program demonstrates the computation of 95 % C.L. limits.
//It uses a set of randomly created histograms.
//
//Author:  [email protected] on 21.08.02

// Create a new canvas.
  TCanvas *c1 = new TCanvas("c1","Dynamic Filling Example",200,10,700,500);
  c1->SetFillColor(42);
  
// Create some histograms
  TH1D* background = new TH1D("background","The expected background",30,-4,4);
  TH1D* signal     = new TH1D("signal","the expected signal",30,-4,4);
  TH1D* data       = new TH1D("data","some fake data points",30,-4,4);
  background->SetFillColor(48);
  signal->SetFillColor(41);
  data->SetMarkerStyle(21);
  data->SetMarkerColor(kBlue);
  background->Sumw2(); // needed for stat uncertainty
  signal->Sumw2(); // needed for stat uncertainty
  
// Fill histograms randomly
  TRandom2 r;
  Float_t bg,sig,dt;
  for (Int_t i = 0; i < 25000; i++) {
     bg  = r.Gaus(0,1);
     sig = r.Gaus(1,.2);
     background->Fill(bg,0.02);
     signal->Fill(sig,0.001);
  }
  for (Int_t i = 0; i < 500; i++) {
     dt = r.Gaus(0,1);
     data->Fill(dt);
  }
  THStack *hs = new THStack("hs","Signal and background compared to data...");
  hs->Add(background);
  hs->Add(signal);
  hs->Draw("hist");
  data->Draw("PE1,Same");
  c1->Modified();
  c1->Update();
  c1->GetFrame()->SetFillColor(21);
  c1->GetFrame()->SetBorderSize(6);
  c1->GetFrame()->SetBorderMode(-1);
  c1->Modified();
  c1->Update();
  gSystem->ProcessEvents();

// Compute the limits
  cout << "Computing limits... " << endl;
  TLimitDataSource* mydatasource = new TLimitDataSource(signal,background,data);
  TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
  cout << "CLs    : "   << myconfidence->CLs()  << endl;
  cout << "CLsb   : "   << myconfidence->CLsb() << endl;
  cout << "CLb    : "   << myconfidence->CLb()  << endl;
  cout << "< CLs >  : " << myconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << myconfidence->GetExpectedCLb_b()  << endl;

// Add stat uncertainty
  cout << endl << "Computing limits with stat systematics... " << endl;
  TConfidenceLevel *mystatconfidence = TLimit::ComputeLimit(mydatasource,50000,true);
  cout << "CLs    : "   << mystatconfidence->CLs()  << endl;
  cout << "CLsb   : "   << mystatconfidence->CLsb() << endl;
  cout << "CLb    : "   << mystatconfidence->CLb()  << endl;
  cout << "< CLs >  : " << mystatconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << mystatconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << mystatconfidence->GetExpectedCLb_b()  << endl;

// Add some systematics
  cout << endl << "Computing limits with systematics... " << endl;
  TVectorD errorb(2);
  TVectorD errors(2);
  TObjArray* names = new TObjArray();
  TObjString name1("bg uncertainty");
  TObjString name2("sig uncertainty");
  names->AddLast(&name1);
  names->AddLast(&name2);
  errorb[0]=0.05; // error source 1: 5%
  errorb[1]=0;    // error source 2: 0%
  errors[0]=0;    // error source 1: 0%
  errors[1]=0.01; // error source 2: 1%
  TLimitDataSource* mynewdatasource  = new TLimitDataSource();
  mynewdatasource->AddChannel(signal,background,data,&errors,&errorb,names);
  TConfidenceLevel *mynewconfidence = TLimit::ComputeLimit(mynewdatasource,50000,true);
  cout << "CLs    : " << mynewconfidence->CLs()  << endl;
  cout << "CLsb   : " << mynewconfidence->CLsb() << endl;
  cout << "CLb    : " << mynewconfidence->CLb()  << endl;
  cout << "< CLs >  : " << mynewconfidence->GetExpectedCLs_b()  << endl;
  cout << "< CLsb > : " << mynewconfidence->GetExpectedCLsb_b() << endl;
  cout << "< CLb >  : " << mynewconfidence->GetExpectedCLb_b()  << endl;

// show canonical -2lnQ plots in a new canvas
// - The histogram of -2lnQ for background hypothesis (full)
// - The histogram of -2lnQ for signal and background hypothesis (dashed)
  TCanvas *c2 = new TCanvas("c2");
  myconfidence->Draw();
  
// clean up (except histograms and canvas)
  delete myconfidence;
//.........这里部分代码省略.........

void MatchComparison()
{
	//
	// Initialize AliRun manager
	//
	
	//
	// Initialize run loader and load Kinematics
	//
	AliRunLoader *runLoader = AliRunLoader::Open("galice.root");
	if (!runLoader) return;
	runLoader->LoadgAlice();
	gAlice = runLoader->GetAliRun();
	runLoader->LoadKinematics();
	
	//
	// Initialize histograms with their error computation
	//
	TH1D *hgood  = new TH1D("hgood", "Well matched tracks", 40, 0.0,  40.0);
	TH1D *hfake  = new TH1D("hfake", "Fake matched tracks", 40, 0.0,  40.0);
	TH1D *htrue  = new TH1D("htrue", "True matches"       , 40, 0.0,  40.0);
	TH1D *hfound = new TH1D("hfound","Found matches"      , 40, 0.0,  40.0);
	hgood->Sumw2();
	hfake->Sumw2();
	htrue->Sumw2();
	hfound->Sumw2();
	
	//
	// Open file containing true matches,
	// retrieve the Tree and link to a cursor.
	//
	TFile *fileTrue = TFile::Open("true-matches.root");
	match_t trueMatch;
	
	//
	// Open file of found matches,
	// link the modified ESD container.
	//
	TFile *fileFound = TFile::Open("matchESD.root");
	TTree *treeFound = (TTree*)fileFound->Get("esdTree");
	AliESDEvent* esd = new AliESDEvent();
	esd->ReadFromTree(treeFound);
	Long64_t nEvents = treeFound->GetEntries();
	
	//
	// Loop on all events
	//
	Int_t im, it, ic, nTrueMatches, nTracks;
	Int_t label, trkLabel, cluLabel;
	for (Long64_t iev = 0; iev < nEvents; iev++) {
		
		// get true matches tree of given event
		TTree *treeTrue = (TTree*)fileTrue->Get(Form("tm_%d", iev));
		treeTrue->SetBranchAddress("matches", &trueMatch);
		nTrueMatches = treeTrue->GetEntries();
		
		// set TTree pointers to selected event
		runLoader->GetEvent(iev);
		treeFound->GetEntry(iev);
		AliStack *stack = runLoader->Stack();
		nTracks = esd->GetNumberOfTracks();
		
		// read all true pairs
		for (im = 0; im < nTrueMatches; im++) {
			treeTrue->GetEntry(im);
			AliESDtrack *track = esd->GetTrack(trueMatch.indexT);
			if (!track) continue;
			
			label = TMath::Abs(track->GetLabel());
			TParticle *p = stack->Particle(label);
			htrue->Fill(p->Pt());
			cout <<"filling true"<< endl;
		}
		
		// compare found matches
		for (Int_t it = 0; it < nTracks; it++) {
			AliESDtrack *track = esd->GetTrack(it);
			ic = track->GetEMCALcluster();
			if (ic == AliEMCALTracker::kUnmatched) continue;
			ic = TMath::Abs(ic);
			AliESDCaloCluster *cl = esd->GetCaloCluster(ic);
			if (!cl) continue;
			if (!cl->IsEMCAL()) continue ;
			trkLabel = TMath::Abs(track->GetLabel());
			cluLabel = cl->GetLabel();
			if (trkLabel == cluLabel && trkLabel >= 0) {
				TParticle *p = stack->Particle(TMath::Abs(trkLabel));
				hgood->Fill(p->Pt());
				hfound->Fill(p->Pt());
				cout <<"filling GOOD, pt:" << p->Pt()<< endl;
			}
			else  {
				TParticle *p = stack->Particle(TMath::Abs(trkLabel));
				hfake->Fill(p->Pt());
				hfound->Fill(p->Pt());
				cout <<"filling FAKE" << endl;
			}
		}
	}
	
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	h1->GetYaxis()->SetLabelSize( 0.03);
	h1->GetXaxis()->CenterTitle();
	h1->GetYaxis()->CenterTitle();
	h1->Draw();

	h1->SetMinimum( 0);
	yy = h1->GetMaximum();
	ym = yy*1.05;
	h1->SetMaximum( ym);
	y1 = yy*0.9;
	y2 = yy;

	l1 = new TLine( prompt[0], 0,prompt[0], ym);
	l1->SetLineWidth(3);
	l1->Draw();
	l2 = new TLine( prompt[1], 0,prompt[1], ym);
	l2->SetLineWidth(3);
	l2->Draw();
	l3 = new TLine( acc[0], 0,acc[0], ym);
	l3->SetLineStyle(2);
	l3->SetLineWidth(3);
	l3->Draw();
	l4 = new TLine( acc[1], 0,acc[1], ym);
	l4->SetLineStyle(2);
	l4->SetLineWidth(3);
	l4->Draw();
	l5 = new TLine( acc[2], 0,acc[2], ym);
	l5->SetLineStyle(2);
	l5->SetLineWidth(3);
	l5->Draw();
	l6 = new TLine( acc[3], 0,acc[3], ym);
	l6->SetLineStyle(2);
	l6->SetLineWidth(3);
	l6->Draw();

	// Photon Missing Mass
	c1->cd( 2);

	name = "COMP_PhotonMmissP_v_TChanPhotP";
	TH2D *h2P = (TH2D*)full.Get( name);
	name = "COMP_PhotonMmissR_v_TChanPhotR";
	TH2D *h2R = (TH2D*)full.Get( name);

	bin_lo = h2P->GetXaxis()->FindBin( chan_lo);
	bin_hi = h2P->GetXaxis()->FindBin( chan_hi);
	TH1D *hP = h2P->ProjectionY( "PhotonMmissP", bin_lo, bin_hi);
	TH1D *hR = h2R->ProjectionY( "PhotonMmissR", bin_lo, bin_hi);
	TH1D *hS = (TH1D*)hP->Clone( "subt");

	// TDC Spectrum
	hP->SetTitle();
	hP->SetLineWidth( 2);
	hP->GetXaxis()->SetTitleOffset( 1.1);
	hP->GetXaxis()->SetTitle( "Proton Missing Mass (MeV)");
	hP->GetXaxis()->SetLabelSize( 0.03);
	hP->GetXaxis()->CenterTitle();
	hP->GetXaxis()->SetRangeUser( 800, 1200);
	hP->Draw();

	hR->SetLineWidth( 2);
	hR->SetLineStyle( 2);
	hR->Scale( 0.0833);
	hR->Draw( "same");

	tl = new TLegend( 0.2, 0.7, 0.4, 0.8);
	tl->SetFillStyle( 0);
	tl->SetBorderSize( 0);
	tl->SetTextSize( 0.05);
	tl->SetTextFont( 132);
	tl->AddEntry( hP, "prompt");
	tl->AddEntry( hR, "random");
	tl->Draw();

	c1->cd( 3);
	hS->Sumw2();
	hS->Add( hR, -1.0);

	hS->SetTitle();
	hS->GetXaxis()->SetTitleOffset( 1.1);
	hS->GetXaxis()->SetTitle( "Proton Missing Mass (MeV)");
	hS->GetXaxis()->SetLabelSize( 0.03);
	hS->GetXaxis()->CenterTitle();
	hS->SetLineWidth( 2);
	hS->SetMarkerStyle( 20);
	hS->SetMarkerSize( 0.5);
	hS->GetXaxis()->SetRangeUser( 800, 1200);
	hS->Draw();

	tl = new TLegend( 0.2, 0.7, 0.4, 0.8);
	tl->SetFillStyle( 0);
	tl->SetBorderSize( 0);
	tl->SetTextSize( 0.05);
	tl->SetTextFont( 132);
	tl->AddEntry( hS, "subtracted", "p");
	tl->Draw();

	name = Form( "plots/Compton/random_subt_%d", ebin);
	name.Append( ".pdf");
	if ( print == kTRUE) c1->Print( name);
}

void MVA_allsrcComp()
{
  //InpFile = new TFile("../srcbatch/scTTMuon/TTsemiMuon.root");
  //InpFile = new TFile("../srcbatch/scTTElectron/TTsemiElectron.root");
  //InpFile = new TFile("../srcbatch/MonoJet_MET/MonoJet.root");
  InpFile = new TFile("../srcbatch/mSUGRA_MET/mSUGRA.root");
  //InpFile = new TFile("../srcbatch/scTT2Hadron/TT2Hadron.root");
  //InpFile = new TFile("../srcbatch/scWToENu/WToENu.root");
  //InpFile = new TFile("../srcbatch/scWToMuNu/WToMuNu.root");
  //InpFile = new TFile("../srcbatch/scDYElEl/DYElEl.root");
  //InpFile = new TFile("../srcbatch/scDYMuMu/DYMuMu.root");
  //InpFile = new TFile("../srcbatch/scZZ4Mu/ZZ4Mu.root");
  //InpFile = new TFile("../srcbatch/scZZ4El/ZZ4El.root");
  //InpFile = new TFile("../srcbatch/scZZMuEl/ZZMuEl.root");
  //InpFile = new TFile("../srcbatch/ZZ2Mu2Nu_MET/ZZ2Mu2Nu.root");
  //InpFile = new TFile("../srcbatch/ZZ2El2Nu_MET/ZZ2El2Nu.root");
  //InpFile = new TFile("../srcbatch/scWW/WW.root");
  //InpFile = new TFile("../srcbatch/scWZ/WZ.root");
  //InpFile = new TFile("../srcbatch/scMonoPhoton/MonoPhoton.root");
  
  //InpFileall = new TFile("../allbatch/TTsemiMuon_MET/TTsemiMuon.root");
  //InpFileall = new TFile("../allbatch/TTsemiElectron_MET/TTsemiElectron.root");
  //InpFileall = new TFile("../allbatch/MonoJet_MET/MonoJet.root");
  InpFileall = new TFile("../allbatch/mSUGRA_MET/mSUGRA.root");
  //InpFileall = new TFile("../allbatch/TT2Hadron_MET/TT2Hadron.root");
  //InpFileall = new TFile("../allbatch/WToENu_MET/WToENu.root");
  //InpFileall = new TFile("../allbatch/WToMuNu_MET/WToMuNu.root");
  //InpFileall = new TFile("../allbatch/DYElEl_MET/DYElEl.root");
  //InpFileall = new TFile("../allbatch/DYMuMu_MET/DYMuMu.root");
  //InpFileall = new TFile("../allbatch/ZZ4Mu_MET/ZZ4Mu.root");
  //InpFileall = new TFile("../allbatch/ZZ4El_MET/ZZ4El.root");
  //InpFileall = new TFile("../allbatch/ZZMuEl_MET/ZZMuEl.root");
  //InpFileall = new TFile("../allbatch/ZZ2Mu2Nu_MET/ZZ2Mu2Nu.root");
  //InpFileall = new TFile("../allbatch/ZZ2El2Nu_MET/ZZ2El2Nu.root");
  //InpFileall = new TFile("../allbatch/WW_MET/WW.root");
  //InpFileall = new TFile("../allbatch/WZ_MET/WZ.root");
  //InpFileall = new TFile("../allbatch/MonoPhoton_MET/MonoPhoton.root");

  char ylabel[100];
  char histName[100];
  char histName_org[100];
  
  TH1D *pfMet;
  TH1D *NoPUmet;
  TH1D *MVAmet;
  TH1D *NoPUmetall;
  TH1D *MVAmetall;
 
    sprintf(histName,"pfMet");
    sprintf(histName_org,"pfMET_0");
    pfMet =(TH1D*)InpFile->Get(histName_org)->Clone(histName);pfMet->Sumw2();
    
    sprintf(histName,"NoPUmet");
    sprintf(histName_org,"NoPuMET_0");
    NoPUmet =(TH1D*)InpFile->Get(histName_org)->Clone(histName);NoPUmet->Sumw2();
    
    sprintf(histName,"NoPUmetall");
    sprintf(histName_org,"NoPuMET_0");
    NoPUmetall =(TH1D*)InpFileall->Get(histName_org)->Clone(histName);NoPUmetall->Sumw2();
    
    sprintf(histName,"MVAmet");
    sprintf(histName_org,"MVaMET_0");
    MVAmet =(TH1D*)InpFile->Get(histName_org)->Clone(histName);MVAmet->Sumw2();

    sprintf(histName,"MVAmetall");
    sprintf(histName_org,"MVaMET_0");
    MVAmetall =(TH1D*)InpFileall->Get(histName_org)->Clone(histName);MVAmetall->Sumw2();

    pfMet->SetLineWidth(2);
    pfMet->SetLineColor(kBlack);
    pfMet->SetMarkerStyle(26);
    pfMet->SetMarkerColor(kBlack);
    NoPUmet->SetLineWidth(2);
    NoPUmet->SetLineColor(kRed);
    NoPUmet->SetMarkerStyle(kFullCircle);
    NoPUmet->SetMarkerColor(kRed);
    NoPUmetall->SetLineWidth(2);
    NoPUmetall->SetLineColor(kRed);
    NoPUmetall->SetMarkerStyle(kCircle);
    NoPUmetall->SetMarkerColor(kRed);
    MVAmet->SetLineWidth(2);
    MVAmet->SetLineColor(kBlue);
    MVAmet->SetMarkerStyle(kFullSquare);
    MVAmet->SetMarkerColor(kBlue);
    MVAmetall->SetLineWidth(2);
    MVAmetall->SetLineColor(kBlue);
    MVAmetall->SetMarkerStyle(25);
    MVAmetall->SetMarkerColor(kBlue);

    sprintf(ylabel,"MET-genMEtTrue",MVAmet->GetBinWidth(1));
    MVAmet->GetYaxis()->SetTitle(ylabel);
    MVAmet->GetYaxis()->SetTitleOffset(1.2);
    MVAmet->GetYaxis()->SetTitleSize(0.04);
    MVAmet->GetYaxis()->SetLabelSize(0.04);
    MVAmet->GetXaxis()->SetTitle("genMEt");
    MVAmet->GetXaxis()->SetTitleOffset(1.0);
    //MVAmet->GetXaxis()->SetTitleSize(0.13);
   
    //MVAmet->GetYaxis()->SetRangeUser(-15,30);//WJet TTsemiPU012
    //MVAmet->GetYaxis()->SetRangeUser(-25,30);// TTsemiPU3
//.........这里部分代码省略.........

//.........这里部分代码省略.........
 
cout << "got samples" << endl;

THStack *hs = new THStack("hs","test");

  hs->Add(QCD_all);
  hs->Add(wjets);
  hs->Add(WZ); 
  hs->Add(WW);
  hs->Add(ZZ);
  hs->Add(DY1);
  hs->Add(DY2);
  hs->Add(T_tW); 
  hs->Add(Tbar_tW);
  hs->Add(tt);
  hs->Add(ttgamma);
  

TH1D* allMC = (TH1D*)ttgamma->Clone("ratio");
  allMC->Add(tt);
  allMC->Add(wjets);
  allMC->Add(DY1);
  allMC->Add(DY2);
  allMC->Add(T_tW);
  allMC->Add(Tbar_tW);
  allMC->Add(ZZ);
  allMC->Add(WW);
  allMC->Add(WZ);
  allMC->Add(QCD_all);
  
TH1D* dataEff = new TH1D("data eff","data eff",9,0,9);
TH1D* mcEff = new TH1D("mc eff","mc eff",9,0,9);

mcEff->Sumw2();
dataEff->Sumw2();

for(int q =1; q<ttgamma->GetNbinsX(); q++){
mcEff->GetXaxis()->SetBinLabel(q, step[q]);
dataEff->SetBinContent(q, data->GetBinContent(q+1)/data->GetBinContent(q));
dataEff->SetBinError(q, dataEff->GetBinContent(q)*sqrt(pow(data->GetBinError(q+1)/data->GetBinContent(q+1),2)+pow(data->GetBinError(q)/data->GetBinContent(q),2)));
mcEff->SetBinContent(q, allMC->GetBinContent(q+1)/allMC->GetBinContent(q));

}  

	
	TCanvas *c2 = new TCanvas("cutflow eff","cutflow eff",600, 500);
  	
	mcEff->SetLineColor(kRed);
	mcEff->Draw();
	dataEff->Draw("Esame");
	
	TLegend *tleg3;
	tleg3 = new TLegend(0.7,0.7,0.8,0.9);
	tleg3->SetTextSize(0.04);
	tleg3->SetBorderSize(0);
	tleg3->SetFillColor(10);
	tleg3->AddEntry(dataEff , "data", "l");
	tleg3->AddEntry(mcEff , "mc", "l");
	
	tleg3->Draw("same");
	
	c2->SaveAs("plots/cutFlow/cutEff_"+Variable+".pdf");
	c2->SaveAs("plots/cutFlow/cutEff_"+Variable+".png");
	delete c2;

  //draw histos to files

//.........这里部分代码省略.........
    //RooFitResult* filters = model.fitTo(dh, "mhe");
    //RooFitResult* filters = cb.fitTo(dh,Range(0,200),"qr");
    //RooFitResult* filters = bw.fitTo(dh,Range(0,100),"qr");
    //RooFitResult* filters = bag.fitTo(dh,"qr");


    TCanvas* canvas = new TCanvas("ZmassHisto","ZmassHisto",0, 0, 1000,700) ;
    canvas->cd() ; //gPad->SetLeftMargin(0.15);
    //gPad->SetLogy();
    RooPlot* frame = x.frame(Title("e #gamma invariant mass fit")) ;
    dh.plotOn(frame,MarkerSize(0.5),Name("data_hist"));  //this will show histogram data points on canvas

    dh.statOn(frame,Layout(0.15,0.37,0.85),What("N")) ;


    //dh.statOn(frame);  //this will display hist stat on canvas

    //sum.plotOn(frame,LineColor(4));//this will show fit overlay on canvas
    //sum.paramOn(frame); //this will display the fit parameters on canvas

    //bag.plotOn(frame,LineColor(2));//this will show fit overlay on canvas 
    //bag.paramOn(frame); //this will display the fit parameters on canvas

    model.plotOn(frame,Components(bag),LineColor(2));//this will show fit overlay on canvas 
    model.plotOn(frame,Components(signal),LineColor(3));//this will show fit overlay on canvas 
    model.plotOn(frame,LineColor(kBlue),Name("main_curve"));//this will show fit overlay on canvas 

    //	model.paramOn(frame); //this will display the fit parameters on canvas
    model.paramOn(frame, Layout(0.6, 0.95, 0.92)); //this will display the fit parameters on canvas
    //model.paramOn(frame, Layout(0.6, 0.99, 0.75));
    //	model.plotOn(frame,LineColor(4));//this will show fit overlay on canvas 
    //	model.paramOn(frame); //this will display the fit parameters on canvas

    RooHist* histogram = frame->getHist("data_hist");
    RooCurve* curve = frame->getCurve("main_curve");
    TH1D* hresidual  = residualHist(histogram,curve);
    hresidual->Sumw2();
    canvas->Divide( 1, 2, .1, .1 );
    TPad* padHisto = (TPad*) canvas->cd(1);
    TPad* padResid = (TPad*) canvas->cd(2);
    double small = 0.1;
    padHisto->SetPad( 0., r , 1., 1. );
    padHisto->SetBottomMargin( small );
    padResid->SetPad( 0., 0., 1., r  );
    padResid->SetBottomMargin( 0.3  );
    padResid->SetTopMargin   ( small );
    padHisto->cd();


    //float fitvalue = frame.GetX()[60];
    //Double_t nX = x.getVal();
    //cout << "nX  : " << nX << endl;
    //results->SetBinContent(results->GetXaxis()->FindBin(60), nX);
    //float fitvalue = results->GetBinContent(59);
    //cout << "59.5daki fit value  : " << fitvalue << endl;

    //cb.plotOn(frame,LineColor(2));//this will show fit overlay on canvas 
    //cb.paramOn(frame); //this will display the fit parameters on canvas

    //bw.plotOn(frame,LineColor(4));//this will show fit overlay on canvas 
    //bw.paramOn(frame); //this will display the fit parameters on canvas


    //gPad->SetLogy();
    cout << "chisquare : " << frame->chiSquare() << endl ;

    //cout << "Total Number of events: " << Z_mass->Integral() << endl;
    //cout << "Number of signal events: " << fsig.getVal() *  Z_mass->Integral() << endl;
    //cout << "Number of background events: " << (1 - fsig.getVal()) * Z_mass->Integral() << endl;

    //Draw all frames on a canvas

    //TPaveLabel *label1 = new TPaveLabel(1,500,15,700,"Chisquare:");
    frame->GetXaxis()->SetTitle("Invariant mass w/ NoElectronVeto photon (in GeV/c^{2})");  
    //frame->GetXaxis()->SetTitle("Invariant mass w/ TIGHT photon (in GeV/c^{2})");  
    //frame->GetXaxis()->SetTitleOffset(1.2);
    frame->Draw();
    //float binsize = Z_mass->GetBinWidth(1);
    //char Bsize[50]; 
    //sprintf(Bsize,"Events per %2.2f",binsize);
    //frame->GetYaxis()->SetTitle(Bsize);  
    //frame->GetYaxis()->SetTitleOffset(1.2);
    //results->Sumw2();
    padResid->cd();
    hresidual->Draw();
    Lines( hresidual );
    hresidual->Draw( "SAME" );

    //frame->Draw();
    //results->Draw();

    //canvas->Update();
    //title->Draw("same");
    //hresidual->Draw();

    canvas->SaveAs("pdf_TagAndProbe/ResidualNumerator_binned_negligible_errorfit.pdf");
    //canvas->SaveAs("pdf_TagAndProbe/ResidualDenominator_binned_negligible_errorfit.pdf");


}

int resbosComp(const TString BaseName)
{
  TString tmpTStr;
  CPlot* pltComp;
  CPlot* pltCompLog;

  TFile *F_29;
  TFile *F_30;
  TFile *F_31;
  TFile *F_32;
  TFile *F_33;
  TFile *F_34;
  TFile *F_35;
  TFile *F_Data;

  F_29 = new TFile("../"+BaseName+"/29lResbos.root");
  F_30 = new TFile("../"+BaseName+"/30lResbos.root");
  F_31 = new TFile("../"+BaseName+"/31lResbos.root");
  F_32 = new TFile("../"+BaseName+"/32lResbos.root");
  F_33 = new TFile("../"+BaseName+"/33lResbos.root");
  F_34 = new TFile("../"+BaseName+"/34lResbos.root");
  F_35 = new TFile("../"+BaseName+"/35lResbos.root");
  F_Data = new TFile("../RstUnfold/Result_"+BaseName+".root");

  TString resultDir = "Result";
  gSystem->mkdir(resultDir,kTRUE);
  TFile f_out(resultDir+"/Resbos_"+BaseName+".root","recreate");

  const int nBins = 14;
  double WptBins[nBins]    = {0.0,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  double WptLogBins[nBins] = {1.0,7.5,12.5,17.5,24,30,40,50,70,110,150,190,250,600};
  
  TH1D *hResbos29 = new TH1D("hResbos29","hResbos29",nBins-1,WptBins);hResbos29->Sumw2();
  TH1D *hResbos30 = new TH1D("hResbos30","hResbos30",nBins-1,WptBins);hResbos30->Sumw2();
  TH1D *hResbos31 = new TH1D("hResbos31","hResbos31",nBins-1,WptBins);hResbos31->Sumw2();
  TH1D *hResbos32 = new TH1D("hResbos32","hResbos32",nBins-1,WptBins);hResbos32->Sumw2();
  TH1D *hResbos33 = new TH1D("hResbos33","hResbos33",nBins-1,WptBins);hResbos33->Sumw2();
  TH1D *hResbos34 = new TH1D("hResbos34","hResbos34",nBins-1,WptBins);hResbos34->Sumw2();
  TH1D *hResbos35 = new TH1D("hResbos35","hResbos35",nBins-1,WptBins);hResbos35->Sumw2();
  TH1D *hData = new TH1D("hData","hData",nBins-1,WptBins);hData->Sumw2();
  
  TH1D *hResbosLog29 = new TH1D("hResbosLog29","hResbosLog29",13,WptLogBins);hResbosLog29->Sumw2();
  TH1D *hResbosLog30 = new TH1D("hResbosLog30","hResbosLog30",13,WptLogBins);hResbosLog30->Sumw2();
  TH1D *hResbosLog31 = new TH1D("hResbosLog31","hResbosLog31",13,WptLogBins);hResbosLog31->Sumw2();
  TH1D *hResbosLog32 = new TH1D("hResbosLog32","hResbosLog32",13,WptLogBins);hResbosLog32->Sumw2();
  TH1D *hResbosLog33 = new TH1D("hResbosLog33","hResbosLog33",13,WptLogBins);hResbosLog33->Sumw2();
  TH1D *hResbosLog34 = new TH1D("hResbosLog34","hResbosLog34",13,WptLogBins);hResbosLog34->Sumw2();
  TH1D *hResbosLog35 = new TH1D("hResbosLog35","hResbosLog35",13,WptLogBins);hResbosLog35->Sumw2();
  TH1D *hDataLog = new TH1D("hDataLog","hDataLog",13,WptLogBins);hDataLog->Sumw2();
  
  TH1D* lResbos29;
  TH1D* lResbos30;
  TH1D* lResbos31;
  TH1D* lResbos32;
  TH1D* lResbos33;
  TH1D* lResbos34;
  TH1D* lResbos35;
  TH1D* lData;

  lResbos29 =(TH1D*)F_29->Get("lResbos")->Clone();
  lResbos30 =(TH1D*)F_30->Get("lResbos")->Clone();
  lResbos31 =(TH1D*)F_31->Get("lResbos")->Clone();
  lResbos32 =(TH1D*)F_32->Get("lResbos")->Clone();
  lResbos33 =(TH1D*)F_33->Get("lResbos")->Clone();
  lResbos34 =(TH1D*)F_34->Get("lResbos")->Clone();
  lResbos35 =(TH1D*)F_35->Get("lResbos")->Clone();
  lData =(TH1D*)F_Data->Get("BornEffCorr")->Clone();

  lData->Scale(1./18.429);
  
  for( int ipt(1);ipt<nBins;ipt++)
  {
    hResbos29->SetBinContent(ipt,lResbos29->GetBinContent(ipt));
    hResbos30->SetBinContent(ipt,lResbos30->GetBinContent(ipt));
    hResbos31->SetBinContent(ipt,lResbos31->GetBinContent(ipt));
    hResbos32->SetBinContent(ipt,lResbos32->GetBinContent(ipt));
    hResbos33->SetBinContent(ipt,lResbos33->GetBinContent(ipt));
    hResbos34->SetBinContent(ipt,lResbos34->GetBinContent(ipt));
    hResbos35->SetBinContent(ipt,lResbos35->GetBinContent(ipt));
    hData->SetBinContent(ipt,lData->GetBinContent(ipt));
    hData->SetBinError(ipt,lData->GetBinError(ipt));
    
    hResbos29->SetBinError(ipt,lResbos29->GetBinError(ipt));
    hResbos30->SetBinError(ipt,lResbos30->GetBinError(ipt));
    hResbos31->SetBinError(ipt,lResbos31->GetBinError(ipt));
    hResbos32->SetBinError(ipt,lResbos32->GetBinError(ipt));
    hResbos33->SetBinError(ipt,lResbos33->GetBinError(ipt));
    hResbos34->SetBinError(ipt,lResbos34->GetBinError(ipt));
    hResbos35->SetBinError(ipt,lResbos35->GetBinError(ipt));
    
    hResbosLog29->SetBinContent(ipt,lResbos29->GetBinContent(ipt));
    hResbosLog30->SetBinContent(ipt,lResbos30->GetBinContent(ipt));
    hResbosLog31->SetBinContent(ipt,lResbos31->GetBinContent(ipt));
    hResbosLog32->SetBinContent(ipt,lResbos32->GetBinContent(ipt));
    hResbosLog33->SetBinContent(ipt,lResbos33->GetBinContent(ipt));
    hResbosLog34->SetBinContent(ipt,lResbos34->GetBinContent(ipt));
    hResbosLog35->SetBinContent(ipt,lResbos35->GetBinContent(ipt));
    hDataLog->SetBinContent(ipt,lData->GetBinContent(ipt));
    hDataLog->SetBinError(ipt,lData->GetBinError(ipt));
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  legendTot->SetTextFont(72);
  legendTot->SetTextSize(0.035);
  legendTot->SetFillColor(0);
  char nev[50];

  TCanvas * canvas = new TCanvas ( "diffxSec", "differential xSec", 1200, 400 );
  gStyle->SetOptStat(0);
  //  canvas->UseCurrentStyle();
  canvas->Divide(3,1);
  canvas->cd(1);
  histo1->SetNormFactor(histo1Entries*invLuminosityEff1*luminosityFactor);
  histo1->Draw();
  sprintf(nev,"number of events in 100 fb^{-1}:");
  sprintf(nev,"%.2f",histo1Entries*invLuminosityEff1*luminosityFactor);
  legend1->AddEntry(histo1,nev,"");
  legend1->Draw();
  canvas->cd(2);
  histo2->SetNormFactor(histo2Entries*invLuminosityEff2*luminosityFactor);
  histo2->Draw();
  legend2->Clear();
  sprintf(nev,"number of events in 100 fb^{-1}:");
  sprintf(nev,"%.2f",histo2Entries*invLuminosityEff2*luminosityFactor);
  legend2->AddEntry(histo2,nev,"");
  legend2->Draw();
  canvas->cd(3);
  histo3->SetNormFactor(histo3Entries*invLuminosityEff3*luminosityFactor);
  histo3->Draw();
  legend3->Clear();
  sprintf(nev,"number of events in 100 fb^{-1}:");
  sprintf(nev,"%.2f",histo3Entries*invLuminosityEff3*luminosityFactor);
  legend3->AddEntry(histo3,nev,"");
  legend3->Draw();

  histoTot->Sumw2();
  histoTot->Add(histo1);
  histoTot->Add(histo2);
  histoTot->Add(histo3);

  TCanvas * canvasTot = new TCanvas ( "inclusivediffxSec", "differential xSec", 1200, 400 );
  gStyle->SetOptStat(0);
  canvasTot->UseCurrentStyle();
  canvasTot->cd();
  histo1->SetTitle("4-body mass (100 fb^{-1})");
  histo1->GetXaxis()->SetTitle("m_{llbb} (GeV)");
  histo1->GetYaxis()->SetTitle("d#sigma/dm_{llbb} events/10GeV");
  histo1->SetMinimum(0.001);
  histo1->SetLineColor(51);
  histo1->Draw();
  histo2->SetLineColor(56);
  histo2->Draw("same");
  histo3->SetLineColor(60);
  histo3->Draw("same");
  histoTot->Draw("same");

  legendTot->Clear();
  sprintf(nev,"Zb\\bar{b}+0jets");
  legendTot->AddEntry(histo1,nev,"l");
  sprintf(nev,"Zb\\bar{b}+1jets");
  legendTot->AddEntry(histo2,nev,"l");
  sprintf(nev,"Zb\\bar{b}+2jets");
  legendTot->AddEntry(histo3,nev,"l");
  sprintf(nev,"inclusive");
  legendTot->AddEntry(histoTot,nev,"l");
  legendTot->Draw();

void drawQGFraction(){

	const double PI = 3.14159;

	TChain *mix = new TChain("mixing_tree");
	mix->Add("/data/kurtjung/JetTrackCorr_skims/2p76TeV_MC_Pythia6/MergedPythia_withPartonFlavor.root");
	//mix->Add("/data/kurtjung/JetTrackCorr_skims/5TeV_MC_Pythia6/*");

	//double xsecs[11] = {5.335E-01, 3.378E-02, 3.778E-03, 4.412E-04, 6.147E-05, 1.018E-05, 2.477E-06, 6.160E-07, 1.088E-07, 3.216E-08, 0}; //pythia6 5.02 tev weights
    double xsecs[11] = {2.043e-01, 1.075E-02, 1.025E-03, 9.865E-05, 1.129E-05, 1.465E-06, 2.837E-07, 5.323E-08, 5.934e-09, 8.125e-10, 0}; //2.76 tev weights
	int recalculatedEntries[10] = {0,0,0,0,0,0,0,0,0,0};
	double pthatbins[11] = {15,30,50,80,120,170,220,280,370,460,9999};

	TFile *fout = new TFile("QGFrac_pythia6_2p76TeV.root","recreate");

	TH1D *quarkFracIncl = new TH1D("quarkFracIncl","",20,120,500); quarkFracIncl->Sumw2();
	TH1D *glueFracIncl = new TH1D("glueFracIncl","",20,120,500); glueFracIncl->Sumw2();
	TH1D *inclJets = new TH1D("inclJets","",20,120,500); inclJets->Sumw2();

	TH1D *quarkFracLead = new TH1D("quarkFracLead","",20,120,500); quarkFracLead->Sumw2();
	TH1D *glueFracLead = new TH1D("glueFracLead","",20,120,500); glueFracLead->Sumw2();
	TH1D *leadJets = new TH1D("leadJets","",20,120,500); leadJets->Sumw2();


	Int_t HBHENoiseFilterResultRun2Loose, pPAprimaryVertexFilter;
	vector<float> *calo_corrpt=0, *calo_jtphi=0;
	vector<int> *calo_refparton_flavor=0;
	float pthat;

	mix->SetBranchAddress("calo_corrpt",&calo_corrpt);
	mix->SetBranchAddress("calo_jtphi",&calo_jtphi);
	mix->SetBranchAddress("calo_refparton_flavor",&calo_refparton_flavor);
	mix->SetBranchAddress("pthat",&pthat);

	mix->SetBranchAddress("HBHENoiseFilterResultRun2Loose",&HBHENoiseFilterResultRun2Loose);
	mix->SetBranchAddress("pPAprimaryVertexFilter",&pPAprimaryVertexFilter);


	TH1D *pthatHisto = new TH1D("pthatHisto","",10,pthatbins);
	mix->Project("pthatHisto","pthat");
	for(int i=0; i<10; i++){
		recalculatedEntries[i] = pthatHisto->GetBinContent(i+1);
		cout << "entries between pthat " << pthatbins[i] << " and " << pthatbins[i+1] << ": " << recalculatedEntries[i] << endl;
		cout << "weight: " << (xsecs[i]-xsecs[i+1])/recalculatedEntries[i] <<endl;
	}

	int totEntries = mix->GetEntries();
	cout << "entries: "<< totEntries << endl;
	totEntries=100000;
	for(int ievt=0; ievt<totEntries; ievt++){
		mix->GetEntry(ievt);
		if(ievt && ievt%10000==0) cout << "entry: " << ievt << endl;

		//if(!HBHENoiseFilterResultRun2Loose || !pPAprimaryVertexFilter) continue;

		int ibin=0;
		double weight=0.;
		while(pthat>pthatbins[ibin]) ibin++;
		ibin--;
		weight = (xsecs[ibin]-xsecs[ibin+1])/recalculatedEntries[ibin];
		if(weight>1){ 
			cout << "xsec: "<< xsecs[ibin]-xsecs[ibin+1] << " entries: " << recalculatedEntries[ibin] << endl;
			cout << "pthat: "<< pthat << " bin " << ibin << endl;
		}
		
		for(unsigned int ijet=0; ijet<calo_corrpt->size(); ijet++){
			
			if(abs(calo_refparton_flavor->at(ijet))>0 && abs(calo_refparton_flavor->at(ijet))<6) quarkFracIncl->Fill(calo_corrpt->at(ijet), weight);
			if(abs(calo_refparton_flavor->at(ijet))==21){ glueFracIncl->Fill(calo_corrpt->at(ijet), weight); }
			inclJets->Fill(calo_corrpt->at(ijet), weight);

			if(ijet==0 && calo_corrpt->size()>1){
				double dphi = abs(calo_jtphi->at(0) - calo_jtphi->at(1));
				if(dphi>(7*PI/8.) && dphi<(9*PI/8.) && calo_corrpt->at(1)>50) {

					if(abs(calo_refparton_flavor->at(ijet))>0 && abs(calo_refparton_flavor->at(ijet))<6) quarkFracLead->Fill(calo_corrpt->at(ijet), weight);
					if(abs(calo_refparton_flavor->at(ijet))==21) glueFracLead->Fill(calo_corrpt->at(ijet), weight);
					leadJets->Fill(calo_corrpt->at(ijet), weight);

				}
			}
		}
	}

	//quarkFracLead->Divide(leadJets);
	//glueFracLead->Divide(leadJets);

	//quarkFracIncl->Divide(inclJets);
	//glueFracIncl->Divide(inclJets);

	fout->cd();

	formatHisto(quarkFracLead,1);
	formatHisto(glueFracLead,2);
	formatHisto(quarkFracIncl,4);
	formatHisto(glueFracIncl,8);

	quarkFracLead->Write();
	glueFracLead->Write();
	quarkFracIncl->Write();
//.........这里部分代码省略.........

TH1D* getSample(TString sample, int rebinFact, TString Obj, TString dir, bool ge4){
		
	TFile* file = new TFile();
	
	//cout << "file: " << "rootFilesV3/"+dir+"/"+ sample + "_5814pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" << endl;
	
	if(dir == "central" || (dir == "Scale_up_tt" && sample != "TTJet") || (dir == "Scale_down_tt" && sample != "TTJet") || ((dir == "Scale_up" && sample != "WJetsToLNu") && (dir == "Scale_up"  && sample != "DYJetsToLL"))  || ((dir == "Scale_down" && sample != "WJetsToLNu") && (dir == "Scale_down"  && sample != "DYJetsToLL")) || (dir == "Match_up_tt" && sample != "TTJet") || (dir == "Match_down_tt" && sample != "TTJet") || ((dir == "Match_up" && sample != "WJetsToLNu") && (dir == "Match_up"  && sample != "DYJetsToLL"))  || ((dir == "Match_down" && sample != "WJetsToLNu") && (dir == "Match_down"  && sample != "DYJetsToLL")) || dir == "UnclusteredEnUp" || dir == "UnclusteredEnDown" || dir == "JetEnUp" || dir == "JetEnDown" || dir == "JetResUp" || dir == "JetResDown" || dir == "TauEnUp" || dir == "TauEnDown" || dir == "MuonEnUp" || dir == "MuonEnDown" || dir == "ElectronEnUp" || dir == "ElectronEnDown"){
	file = new TFile("rootFilesV4/central/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	}else if(dir == "JES_up")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_plusJES.root");
	else if(dir == "JES_down")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_minusJES.root");
	else if(dir == "BJet_up")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_plusBjet.root");
	else if(dir == "BJet_down")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_minusBJet.root");
	else if(dir == "PU_up")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_PU_72765mb.root");
	else if(dir == "PU_down")
	file = new TFile("rootFilesV4/"+ dir +"/"+ sample + "_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_PU_65835mb.root");
	else if(dir == "Scale_up_tt" && sample == "TTJet")
	file = new TFile("rootFilesV4/central/TTJets-scaleup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Scale_down_tt" && sample == "TTJet")
	file = new TFile("rootFilesV4/central/TTJets-scaledown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Scale_up" && sample == "WJetsToLNu")
	file = new TFile("rootFilesV4/central/WJets-scaleup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Scale_up" && sample == "DYJetsToLL")
	file = new TFile("rootFilesV4/central/ZJets-scaleup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Scale_down" && sample == "WJetsToLNu")
	file = new TFile("rootFilesV4/central/WJets-scaledown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Scale_down" && sample == "DYJetsToLL")
	file = new TFile("rootFilesV4/central/ZJets-scaledown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_up_tt" && sample == "TTJet")
	file = new TFile("rootFilesV4/central/TTJets-matchingup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_down_tt" && sample == "TTJet")
	file = new TFile("rootFilesV4/central/TTJets-matchingdown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_up" && sample == "WJetsToLNu")
	file = new TFile("rootFilesV4/central/WJets-matchingup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_up" && sample == "DYJetsToLL")
	file = new TFile("rootFilesV4/central/ZJets-matchingup_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_down" && sample == "WJetsToLNu")
	file = new TFile("rootFilesV4/central/WJets-matchingdown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(dir == "Match_down" && sample == "DYJetsToLL")
	file = new TFile("rootFilesV4/central/ZJets-matchingdown_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root");
	else if(sample == "TTJet_MCNLO")
	file = new TFile("rootFilesV4/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_MCatNLO.root");
	else if(sample == "TTJet_POWHEG")
	file = new TFile("rootFilesV4/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET_POWHEG.root");
	
	TH1D* plot; 
	TH1D* plot2;
	TH1D* plot3;
	TH1D* plot4;
	TH1D* plot5;
	
	
	cout<< "folder: " << "TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_2btags" << endl;
	
	plot = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_2btags");
	plot2 = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_3btags");
	plot3 = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_4orMoreBtags");
	plot4 = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_0btag");
	plot5 = (TH1D*) file->Get("TTbar_plus_X_analysis/MuPlusJets/Ref selection/"+Obj+"/muon_absolute_eta_1btag");		
	plot->Sumw2();
	plot2->Sumw2();
	plot3->Sumw2();
	plot4->Sumw2();
	plot5->Sumw2();
	
	plot->Add(plot2);
	plot->Add(plot3);
	
	if(ge4 ==true){
	plot->Add(plot4);
	plot->Add(plot5);
	}
	
        if(dir == "central"){
	//plot->SetFillColor(kRed+1);
        plot->SetLineColor(kRed+1);
	}else if(dir == "JES_up"){
	plot->SetLineColor(kGreen-3);	  
  	//plot->SetFillColor(kGreen-3);
	}else if(dir == "JES_down"){
	//plot->SetFillColor(kAzure-2);
	plot->SetLineColor(kAzure-2);
	}
	
	plot->Scale(1/plot->Integral());
	plot->Rebin(rebinFact);
	
	plot->SetDirectory(gROOT);
	
	file->Close();
	return plot;

}

TH1D * fillPlot2012_gammaJets::Plot(string var, string name, int nbin, double min, double max, int signal) {
  
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
  
  TH1D * tempplot = new TH1D(name.c_str(),name.c_str(),nbin,min,max);
  tempplot->Sumw2();

  ofstream outfile;

  if (writetxt != "") {
    string filename(writetxt);
    outfile.open(filename.c_str()); 
  }

  // photonID MVA
  tmvaReaderID_Single_Endcap=0;
  tmvaReaderID_Single_Barrel=0;
  if (!tmvaReaderID_Single_Barrel || !tmvaReaderID_Single_Endcap) SetAllMVA();  

  // Loop over entries
  int enMax = nentries; 
  for (Long64_t jentry=0; jentry<enMax;jentry++) {
    Long64_t ientry = LoadTree(jentry);

    if (jentry%100000==0) cout << jentry << endl;

    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;

    // patological events
    if (npu>=60) continue;    

    // test to remove possible duplicated events between samples
    // if ( signal<6 && (ISRGamma>0 || FSRGamma>0) ) continue; // chiara

    // first vertex must be good
    if (vtxId<0) continue;

    // HLT selection - for data only
    if ( signal==100 && hltcut==30 && !isHLT_30() )  continue;
    if ( signal==100 && hltcut==50 && !isHLT_50() )  continue;
    if ( signal==100 && hltcut==75 && !isHLT_75() )  continue;
    if ( signal==100 && hltcut==90 && !isHLT_90() )  continue;

    // vector with index of fully selected gammas
    vector<int> fullSelected;

    // apply the full cut based photonID selection (medium WP)
    // https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
    float EA_chargedH[7] = { 0.012, 0.010, 0.014, 0.012, 0.016, 0.020, 0.012};
    float EA_neutralH[7] = { 0.030, 0.057, 0.039, 0.015, 0.024, 0.039, 0.072};
    float EA_photons[7]  = { 0.148, 0.130, 0.112, 0.216, 0.262, 0.260, 0.266};

    int theNumber = nPhot_presel;
    if (nPhot_presel>100) theNumber = 100;
    for (int theGamma=0; theGamma<theNumber; theGamma++) {

      // for effective area calculation
      int theEAregion = effectiveAreaRegion(etaPhot_presel[theGamma]); 
      if (theEAregion>6) continue;

      // ---------------------------------------------------------------------------
      // preselection to really emulate the single gamma HLT requirements (current preselection is H->gg one)
      bool isReallyPresel = true;      
      float preselECAL    = pid_jurECAL03_presel[theGamma]  - 0.012*ptPhot_presel[theGamma];  
      float preselHCAL    = pid_twrHCAL03_presel[theGamma]  - 0.005*ptPhot_presel[theGamma]; 
      float preselTracker = pid_hlwTrack03_presel[theGamma] - 0.002*ptPhot_presel[theGamma]; 
      if ( preselECAL > 5.5)    isReallyPresel = false;
      if ( preselHCAL > 3.5)    isReallyPresel = false;
      if ( preselTracker > 3.5) isReallyPresel = false;
      if ( theEAregion<2) {  // EB
	if ( pid_HoverE_presel[theGamma]>0.15 )   isReallyPresel = false;       
	if ( sEtaEtaPhot_presel[theGamma]>0.024 ) isReallyPresel = false;
      } else {     // EE
	if(pid_HoverE_presel[theGamma]>0.10)      isReallyPresel = false;       
	if (sEtaEtaPhot_presel[theGamma]>0.040)   isReallyPresel = false;
      }
      if( !isReallyPresel ) continue;

      // ---------------------------------------------------------------------------
      /*
      // full 2012 egamma pog cut based selection
      bool isFullySel = true;

      // isolations corrected with effective areas
      float rhoCorrCharged = pid_pfIsoCharged03ForCiC_presel[theGamma] - rhoAllJets*EA_chargedH[theEAregion];   
      float rhoCorrNeutral = pid_pfIsoNeutrals03ForCiC_presel[theGamma] - rhoAllJets*EA_neutralH[theEAregion];
      float rhoCorrPhoton  = pid_pfIsoPhotons03ForCiC_presel[theGamma] - rhoAllJets*EA_photons[theEAregion];
      if (rhoCorrCharged<0) rhoCorrCharged = 0.;
      if (rhoCorrNeutral<0) rhoCorrNeutral = 0.;
      if (rhoCorrPhoton<0)  rhoCorrPhoton  = 0.;

      // if(pid_hasMatchedPromptElephot[theGamma]) isFullySel = false;   // already applied at preselection level
      if (pid_HoverE_presel[theGamma]>0.05)     isFullySel = false;       
      if (theEAregion<2) {  // EB
	if (sEtaEtaPhot_presel[theGamma]>0.011) isFullySel = false;
	if (rhoCorrCharged > 1.5)               isFullySel = false;
	if (rhoCorrNeutral > 1.0 + 0.04*ptPhot_presel[theGamma])  isFullySel = false;
	if (rhoCorrPhoton  > 0.7 + 0.005*ptPhot_presel[theGamma]) isFullySel = false;
//.........这里部分代码省略.........

void QCDxcheck_shape() {

    setTDRStyle();

//MET will need choice of variable at the top
//  TString Variable = "MET";
//  int Nbins = 6;
//  TString bins[6] = {"_bin_0-25", "_bin_25-45", "_bin_45-70", "_bin_70-100", "_bin_100-150", "_bin_150-inf"};
//  double xbins[7] = {0,25,45,70,100,150, 250};
//  TString Xtitle = "E_{T}^{miss} ";

    TString Variable ="WPT";
    int Nbins = 6;
    TString bins[6] = {"_bin_0-40", "_bin_40-70", "_bin_70-100", "_bin_100-130", "_bin_130-170", "_bin_170-inf"};
    double xbins[7] = {1,40,70,100,130,170,220};
    TString Xtitle = "p_{T}(W)";

//  TString Variable ="MT";
//  int Nbins = 5;
//  TString bins[5] = {"_bin_0-30", "_bin_30-50", "_bin_50-80", "_bin_80-100", "_bin_100-inf"};
//  double xbins[6] = {1,30,50,80,100,150};
//  TString Xtitle = "M(W)_{T}";

//  TString Variable ="ST";
//  int Nbins = 7;
//  TString bins[7] = {"_bin_0-350","_bin_350-400", "_bin_400-450", "_bin_450-500", "_bin_500-580", "_bin_580-700", "_bin_700-inf"};
//  double xbins[8] = {1,350,400,450,500,580,700,1000};
//  TString Xtitle = "S_{T}";

// TString Variable ="HT";
// int Nbins = 7;
// TString bins[7] = {"_bin_0-240", "_bin_240-280", "_bin_280-330", "_bin_330-380", "_bin_380-450", "_bin_450-600", "_bin_600-inf"};
// double xbins[8] = {1,240,280,330,380,450,600,800};
// TString Xtitle = "H_{T}";

    //differential histo
    TH1D *qcd_data = new TH1D("qcd data", "", Nbins, xbins);


    for(int i = 0; i < Nbins; i++) {

        TString Obj;

        if(Variable == "MET") {
            Obj = "Binned_"+Variable+"_Analysis/patType1CorrectedPFMet"+bins[i]+"/muon_absolute_eta_";
        } else if(Variable == "HT") {
            Obj = "Binned_"+Variable+"_Analysis/HT"+bins[i]+"/muon_absolute_eta_";
        } else {
            Obj = "Binned_"+Variable+"_Analysis/"+Variable+"_with_patType1CorrectedPFMet"+bins[i]+"/muon_absolute_eta_";
        }

        double error;

        double content = getQCD(Variable, Obj, 1, &error);
        qcd_data->SetBinContent(i+1, content);
        qcd_data->SetBinError(i+1, error);

    }

    qcd_data->Sumw2();
    qcd_data->Scale(1./qcd_data->Integral());

    for(int i = 0; i < Nbins; i++) {
        cout << qcd_data->GetBinContent(i+1) << ", ";

    }

    cout << " " << endl;

//get fit results
    TFile* fit_file = new TFile("outFiles/diffResults_"+Variable+"_bin_0-40.root");
    TH1D*  fit = (TH1D*) fit_file->Get("central_dir/central_qcd_fit");
    fit->SetLineColor(kRed);
    fit->SetMarkerColor(kRed);
    fit->SetMarkerSize(0.1);
    fit->Sumw2();
    fit->Scale(1./fit->Integral());

    TFile* fit_200f = new TFile("outFiles/diffResults_"+Variable+".root");
    TH1D*  fit_200 = (TH1D*) fit_200f->Get("central_dir/central_qcd_fit");
    fit_200->SetLineColor(kGreen);
    fit_200->SetMarkerColor(kGreen);
    fit_200->SetMarkerSize(0.1);
    fit_200->Sumw2();
    fit_200->Scale(1./fit_200->Integral());

    TFile* fit_noconf = new TFile("outFiles/diffResults_"+Variable+"_bin_100-130.root");
    TH1D*  fit_nocon = (TH1D*) fit_noconf->Get("central_dir/central_qcd_fit");
    fit_nocon->SetLineColor(kMagenta);
    fit_nocon->SetMarkerColor(kMagenta);
    fit_nocon->SetMarkerSize(0.1);
    fit_nocon->Sumw2();
    fit_nocon->Scale(1./fit_nocon->Integral());

    TFile* fit_dataf = new TFile("outFiles/diffResults_"+Variable+"_bins.root");
    TH1D*  fit_data = (TH1D*) fit_dataf->Get("central_dir/central_qcd_fit");
    fit_data->SetLineColor(kYellow);
    fit_data->SetMarkerColor(kYellow);
    fit_data->Sumw2();
    fit_data->Scale(1./fit_data->Integral());
//.........这里部分代码省略.........