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

//--------------------------------------------------------------------------------------------------
// Get Total Number of Events in the sample
//--------------------------------------------------------------------------------------------------
Double_t getNormalizationWeight(TString filename, TString datasetName) {
  // Get Normalization Weight Factor

  //Get Number of Events in the Sample
  TFile *file = new TFile(filename.Data(),"READ");
  if (!file) {
    cout << "Could not open file " << filename << endl;
    return 0;
  }

  //TDirectory *dir = (TDirectory*)file->FindObjectAny("AnaFwkMod");
  //if (!dir) {
  //  cout << "Could not find directory AnaFwkMod"
  //       << " in file " << filename << endl;
  //  delete file;
  //  return 0;
  //}

  TH1D *hist = (TH1D*) gROOT->FindObject("hDAllEvents");
  if (!hist) {
    cout << "Could not find histogram hDAllEvents in directory AnaFwkMod"
         << " in file " << filename << endl;
    //delete dir;
    delete file;
    return 0;
  }
  Double_t NEvents = hist->Integral();
  cout << "Original events in the sample: " << NEvents << endl;

  //Get CrossSection
  mithep::SimpleTable xstab("$CMSSW_BASE/src/MitPhysics/data/xs.dat");
  Double_t CrossSection = xstab.Get(datasetName.Data());
  Double_t Weight = CrossSection / NEvents;
  // weight for data is always 1 (-1 to make a trick for fakes)
  if(CrossSection < 0) Weight = -1.0;

  // Only gg->H and qqH samples
  if(datasetName.Contains("-gf-") ||
     datasetName.Contains("-vbf-")){
    // Only high mass samples
    if(datasetName.Contains("h250") || datasetName.Contains("h300") || datasetName.Contains("h350") || datasetName.Contains("h400") || 
       datasetName.Contains("h450") || datasetName.Contains("h500") || datasetName.Contains("h550") || datasetName.Contains("h600") ||
       datasetName.Contains("h700") || datasetName.Contains("h800") || datasetName.Contains("h900") || datasetName.Contains("h1000")){
      TH1D *histBeforeWeight = (TH1D*) gROOT->FindObject("hDHKFactor_1");
      TH1D *histAfterWeight  = (TH1D*) gROOT->FindObject("hDHKFactor_2");
      cout << "Nevents before/after reweighting: " << histBeforeWeight->GetSumOfWeights() << " " << 
                                                      histAfterWeight->GetSumOfWeights() << endl;
      Weight = Weight * histBeforeWeight->GetSumOfWeights() / histAfterWeight->GetSumOfWeights();
    }
  }

  //delete dir;
  delete file;
  return Weight;

}

void NTupleSS() {
    double Lumi = 19700.;
    double signal_weight = 2.39163;//Lumi/(3263180./1177.3);
    double ttbar_weight  = Lumi / (4246440. / 23.64);
    double tautau_weight = Lumi / (3297032. / 1966.7);
    double tbarw_weight  = Lumi / (493460. / 11.1);
    double tw_weight     = Lumi / (497658. / 11.1);
    double ww_weight     = Lumi / (10000430. / 54.94);
    double wz_weight     = Lumi / (10000280. / 33.21);
    double zz_weight     = Lumi / (9799908. / 17.7);

    TH1D* Data = GetDataPhiStar();
    TH1D* h_ee    = GetBGPhiStar(File_Signal_reco,    signal_weight);
    TH1D* h_zz    = GetBGPhiStar(File_zz,    zz_weight);
    TH1D* h_tt    = GetBGPhiStar(File_tt,    ttbar_weight);
    TH1D* h_tautau = GetBGPhiStar(File_tautau, tautau_weight);
    TH1D* h_tbarw = GetBGPhiStar(File_tbarw, tbarw_weight);
    TH1D* h_tw    = GetBGPhiStar(File_tw,    tw_weight);
    TH1D* h_ww    = GetBGPhiStar(File_ww,    ww_weight);
    TH1D* h_wz    = GetBGPhiStar(File_wz,    wz_weight);
    //  TH1D* h_zz    =GetBGPhiStar(File_zz,    zz_weight);

    double data_sel = Data->GetSumOfWeights();
    double ee_sel = h_ee->GetSumOfWeights();
    double tt_sel = h_tt->GetSumOfWeights();
    double tautau_sel = h_tautau->GetSumOfWeights();
    double tbarw_sel = h_tbarw->GetSumOfWeights();
    double tw_sel = h_tw->GetSumOfWeights();
    double ww_sel = h_ww->GetSumOfWeights();
    double wz_sel = h_wz->GetSumOfWeights();
    double zz_sel = h_zz->GetSumOfWeights();

    // double t_bg=ee_sel+tt_sel+tautau_sel+tbarw_sel+tw_sel+ww_sel+wz_sel+zz_sel;
    double t_bg = ee_sel + tt_sel + tautau_sel + tbarw_sel + tw_sel + ww_sel +
                  wz_sel + zz_sel;
    cout << "total left:" << data_sel - t_bg;
    cout << "Data: " << data_sel << "  ee: " << ee_sel << "  tt: " << tt_sel <<
         "  tautau: " << tautau_sel << "  tbarw: " << tbarw_sel << "  tw: " << tw_sel <<
         " singletop: " << tbarw_sel + tw_sel << "  ww: " << ww_sel << "  wz: " << wz_sel
         << "  zz: " << zz_sel << endl;
    cout << "ratio:"          << "  ee: " << ee_sel / data_sel << "  tt: " << tt_sel
         / data_sel << "  tautau: " << tautau_sel / data_sel << "  tbarw: " << tbarw_sel
         / data_sel << "  tw: " << tw_sel / data_sel << " singletop: " <<
         (tbarw_sel + tw_sel) / data_sel << "  ww: " << ww_sel / data_sel << "  wz: " <<
         wz_sel / data_sel << "  zz: " << zz_sel / data_sel << "data: " <<
         (data_sel - t_bg) / data_sel << endl;

    TH1D* h_qcd  = (TH1D*)Data->Clone();
    h_qcd->Add(h_ee, -1.0);
    h_qcd->Add(h_tt, -1.0);
    h_qcd->Add(h_tautau, -1.0);
    h_qcd->Add(h_tbarw, -1.0);
    h_qcd->Add(h_tw, -1.0);
    h_qcd->Add(h_ww, -1.0);
    h_qcd->Add(h_wz, -1.0);
    h_qcd->Add(h_zz, -1.0);

    h_qcd->Draw();
    cout << "norm=" << h_qcd->GetSumOfWeights() << endl;;
}

double getweight(TFile *file, double xsec) {
 
  TDirectory *dir = (TDirectory*)file->FindObjectAny("AnaFwkMod");
  TH1D *hallevts = (TH1D*)dir->Get("hDTotalMCWeight");
  
  return xsec/hallevts->GetSumOfWeights();
  
}

//.........这里部分代码省略.........
  gDirectory->pwd();

  gDirectory->ReadAll(); // load histos

  TList * lst = gDirectory->GetList();
  cout << lst->GetName() << endl;
  cout << lst->GetTitle() << endl;
  cout << "size    " << lst->GetSize() << endl;
  cout << "entries " << lst->GetEntries() << endl;
  cout << "last    " << lst->LastIndex() << endl;

  TIterator *iter = lst->MakeIterator();
  int ii = 0;
  TObject *obj;
  TH1D *h;
  TH1D *h0;
  TH2D *H;
  TH2D *H0;

  while( obj = iter->Next() ){
    ii++;
    cout << setw(4) << ii << ": ";
    cout << obj->ClassName() << " ";
    cout << obj->InheritsFrom("TH1D") << " ";
    cout << obj->GetName() << " \"";
    cout << obj->GetTitle() << "\"";
    cout << endl;
    //    if( obj->ClassName() == "TH1D" ){
    if( obj->InheritsFrom("TH1D") ){
      h = (TH1D*) obj;
      cout << "       1D";
      cout << h->GetNbinsX() << " bins, ";
      cout << h->GetEntries() << " entries, ";
      cout << h->GetSumOfWeights() << " inside, ";
      cout << h->GetBinContent(0) << " under, ";
      cout << h->GetBinContent(h->GetNbinsX()+1) << " over";
      cout << endl;

      f0.cd(); // output file

      //      TH1D* h0 = (TH1D*) h->Clone();
      h0 = h; // copy
      h0->Write(); // write to file f0
      
      f1.cd(); // back to file 1 for the loop
    }
    else{

      if( obj->InheritsFrom("TH2D") ){

	H = (TH2D*) obj;
	cout << "       2D";
	cout << H->GetNbinsX() << " bins, ";
	cout << H->GetEntries() << " entries, ";
	cout << H->GetSumOfWeights() << " inside, ";
	cout << H->GetBinContent(0) << " under, ";
	cout << H->GetBinContent(H->GetNbinsX()+1) << " over";
	cout << endl;
	
	f0.cd(); // output file
	
	H0 = H; // copy
	H0->Write(); // write to file f0
	
	f1.cd(); // back to file 1 for the loop
      }

void compIsol(int varnum)
{
  TString varname, vartitle;
  if(varnum==1) {
    varname="hHEIsoPt";
    vartitle="H/E Isolation";
  }
  if(varnum==2) {
    varname="hHcalIsoPt";
    vartitle="HCAL Isolation";
  }
  if(varnum==3) {
    varname="hEcalIsoPt";
    vartitle="ECAL Isolation";
  }
  if(varnum==4) {
    varname="hTrkIsoPt";
    vartitle="Track Isolation";
  }
  if(varnum==5) {
    varname="hJetEt";
    vartitle="Jet E_{T} [GeV]";
  }
  if(varnum==6) {
    varname="hNJets";
    vartitle="Number of Jets";
  }
  
  
  TH1D* hBorn;
  TH1D* hADD;
  if(varnum<=4) {
    TH2D* hBorn2d = dynamic_cast<TH2D*>(GetPlot(1, varname));
    TH2D* hADD2d = dynamic_cast<TH2D*>(GetPlot(3, varname));
    hBorn = hBorn2d->ProjectionX(varname+"born");
    hADD = hADD2d->ProjectionX(varname+"add");
    cout << "asdf" << endl;
  } else {
    hBorn = dynamic_cast<TH1D*>(GetPlot(1, varname));
    hADD = dynamic_cast<TH1D*>(GetPlot(3, varname));
  }
  
  hBorn->Scale(1/hBorn->GetSumOfWeights());
  hADD->Scale(1/hADD->GetSumOfWeights());

  hBorn->SetLineColor(kBlue+2);
  hADD->SetLineColor(kCyan+4);
  hBorn->SetFillColor(0);
  hADD->SetFillColor(0);
  hBorn->SetLineStyle(1);
  hADD->SetLineStyle(2);
  
  //  hBorn->SetStats(0);
  //  hADD->SetStats(0);

  gPad->SetLogy(1);
  
  hBorn->GetXaxis()->SetTitle(vartitle);
  hBorn->GetYaxis()->SetTitle("Normalized Units");
  hBorn->SetTitle("Leading Photons in Pythia/Sherpa Samples");

  hBorn->Draw();
  hADD->Draw("sames");

  hBorn->SetMaximum(1);

  gPad->Update();
  TPaveStats *st1=(TPaveStats*)hBorn->GetListOfFunctions()->FindObject("stats");
  TPaveStats *st2=(TPaveStats*)hADD->GetListOfFunctions()->FindObject("stats");
  st1->SetName("Born");
  st2->SetName("ADD");
  st1->SetOptStat(101100);
  st2->SetOptStat(101100);
  st1->SetX1NDC(.25);
  st1->SetX2NDC(.55);
  st1->SetY1NDC(.56);
  st1->SetY2NDC(.80);
  st2->SetX1NDC(.56);
  st2->SetX2NDC(.86);
  st2->SetY1NDC(.56);
  st2->SetY2NDC(.80);

  TLegend* leg=new TLegend(.25,.78,.55,.85);
  leg->AddEntry(hBorn, "Pythia (Born)", "l");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->SetTextSize(0.045);
  leg->Draw();

  TLegend* leg=new TLegend(.56,.78,.86,.85);
  leg->AddEntry(hADD, "Sherpa (ADD)", "l");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->SetTextSize(0.045);
  leg->Draw();

  gPad->Update();
  

  return;
//.........这里部分代码省略.........

//.........这里部分代码省略.........

	TFile * file = new TFile(TStrName+TString(".root"),"recreate");
	file->cd("");
	*/
	std::string initNtupleName("initroottree/AC1B");
	TH1D * inputEventsH = new TH1D("inputEventsH","",1,-0.5,0.5);
	TH1D * histWeightsH = new TH1D("histWeightsH","",1,-0.5,0.5);
	TH1D * histWeightsSkimmedH = new TH1D("histWeightsSkimmedH","",1,-0.5,0.5);

	// Histograms after selecting unique dimuon pair
	TH1D * massSelH = new TH1D("massSelH","",200,0,200);
	TH1D * metSelH  = new TH1D("metSelH","",200,0,400);

	TH1D * hDiJetmet = new TH1D("hDiJetmet","",200,0,400);
	TH1D * hDiJetmass = new TH1D("hDiJetmass","",500,0,1000);
	TH1D * hHT_ = new TH1D("hHT_","",800,0,1600);
	TH1D * metAll  = new TH1D("metAll","",200,0,400);

	TH1D * muon1PtH = new TH1D("muon1PtH","",200,0,400);
	TH1D * muon2PtH = new TH1D("muon2PtH","",200,0,400);

	TH1F * NumberOfVerticesH = new TH1F("NumberOfVerticesH","",51,-0.5,50.5);


	//*****  create eta histogram with eta ranges associated to their names (eg. endcap, barrel)   ***** //

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


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

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

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


	PileUp * PUofficial = new PileUp();

	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);


	TFile *f10= new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfDataBarrel);  // mu SF barrel data
	TFile *f11 = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfDataEndcap); // mu SF endcap data
	TFile *f12= new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfMcBarrel);  // mu SF barrel MC
	TFile *f13 = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfMcEndcap); // mu SF endcap MC 

	TGraphAsymmErrors *hEffBarrelData = (TGraphAsymmErrors*)f10->Get("ZMassBarrel");
	TGraphAsymmErrors *hEffEndcapData = (TGraphAsymmErrors*)f11->Get("ZMassEndcap");
	TGraphAsymmErrors *hEffBarrelMC = (TGraphAsymmErrors*)f12->Get("ZMassBarrel");
	TGraphAsymmErrors *hEffEndcapMC = (TGraphAsymmErrors*)f13->Get("ZMassEndcap");

	double * dataEffBarrel = new double[10];
	double * dataEffEndcap = new double[10];
	double * mcEffBarrel = new double[10];
	double * mcEffEndcap = new double[10];

	dataEffBarrel = hEffBarrelData->GetY();

void
//Impose (TDirectory * target, TList * sourcelist, string & np_title_, vector<string> titles,vector<float> xsecs)
Impose (TList * sourcelist, string & np_title_, vector<string> titles,vector<float> xsecs, TString &variable)
{
	cout << "	" << "========================================================" << endl;
	cout << "	" << "This is a macro to superimpose plots of different root files." << endl;
	cout << "	" << "Only TH2Dobjects are superimposed." << endl;
	float Lumi=1;


	Lumi = 15712.;
	bool norm_=false;
        int MaxEventsBin = 10;
	cout<<titles[0]<<"   "<<titles.size()<<endl;

	//not really useful if plots already weighted to lumi - usefull is plots are in a.u.
	vector <float > lumiweights;
	lumiweights.clear();
//	for (unsigned int kk=0; kk<signal_names.size();kk++){
//		cout<<" HERE is some signal  ===============================================  "<<signal_names[kk]<<"  "<<signalnames[kk]<<endl;
//	}

	TH2D* allbkg, *htt,*hstop,*hwj,*hdyj,*hrare,*hdib,*hqcd,*httx, *hrest;

	TFile *first_source = (TFile *) sourcelist->First ();
	first_source->cd ("mutau");

	TH1D* eventCount = (TH1D*)first_source->Get("mutau/histWeightsH");
	//TH1D* eventCount = (TH1D*)first_source->Get("mutau/inputEventsH");
	//TH1D* hxsec = (TH1D*)first_source->Get("mutau/xsec");
	float nGen = eventCount->GetSumOfWeights();
	float xsec = 1;//hxsec->GetMean();
	float norm = xsec*Lumi/nGen;

	norm =1;
	lumiweights.push_back(float(norm));


	//cout<< " for first source file, there where "<<nGen<<" events with xsec "<<xsec<<" weight "<<lumiweights[0]<<endl;//" weight "<<float(xsec*Lumi/nGen)<<"  norm "<<norm<<endl;

	TDirectory *current_sourcedir = gDirectory;
	//gain time, do not add the objects in the list in memory
	Bool_t status = TH1::AddDirectoryStatus ();
	TH1::AddDirectory (kFALSE);
	// loop over all keys in this directory
	TChain *globChain = 0;
	TIter nextkey (current_sourcedir->GetListOfKeys ());
	//TIter nextkey (((TDirectory *) current_sourcedir->Get ("ana"))->GetListOfKeys ());
	TKey *key, *oldkey = 0;
	while ((key = (TKey *) nextkey ())) {
	//variable="met_MTsum_16";

	int count=0;
		count++;
		//if (count>20) break;
		//keep only the highest cycle number for each key
		//        if (oldkey && !strcmp (oldkey->GetName (), key->GetName ()))
		//            continue;

		// read object from first source file and create a canvas
		// first_source->cd (path);
		first_source->cd ("mutau");
		TObject *obj = key->ReadObj ();

		//string nn = obj->GetName();
		// if (std::string::npos == nn.find("Cut")) continue;
		//cout<<obj->GetName()<<endl;



		string nn = obj->GetName();
		bool flcont = true;
		bool NormTT = false;
		if (string::npos != nn.find("_9") || string::npos != nn.find("_10") || string::npos != nn.find("_11") || string::npos != nn.find("_12") || string::npos != nn.find("_13") || string::npos != nn.find("_14") || string::npos != nn.find("_15") || string::npos != nn.find("_16") || string::npos != nn.find("_18") || string::npos != nn.find("_18") || string::npos != nn.find("_18")) NormTT = true;
		//if ( string::npos != nn.find("_11") || string::npos != nn.find("_12") || string::npos != nn.find("_13") || string::npos != nn.find("_14")) NormTT = true;
		NormTT=true;
		//if ( string::npos == nn.find("CutFlowUnW")) flcont=false;
		//if (string::npos == nn.find(""+variable) ) flcont=false;
		if (string::npos == nn.find(variable) ) continue;
		//		if (!flcont) continue;
		if (obj->IsA ()->InheritsFrom ("TTree") ) continue;
		//	if (obj->IsA ()->InheritsFrom ("TH1") ) continue;
		if (obj->IsA ()->InheritsFrom ("TH2") ) {

			cout<<"=================================================== OK for variable "<<variable<<endl;


			TH2D* hh[1500];
			TH2D* hsignal[1500];
			TH2D* h1 = (TH2D*) obj;

			ModifyHist (h1,1,Lumi,lumiweights[0],titles[0],norm_);
			TFile *nextsource = (TFile *) sourcelist->After (first_source);

			int cl, countsignal;
			h1->SetStats(000000);
			cl=1;
			countsignal=1;

			hh[cl]=h1;
//.........这里部分代码省略.........

mainClass(int luminosity){//constructor
//Importnat
//make sure this initialization of the 
//maps is the same as that in main.cpp
     cutname[0]="RA2nocut";
    cutname[1]="RA2Asys";
    cutname[2]="RA2Inc3Jetcut";
    cutname[3]="RA2HT500cut";
    cutname[4]="RA2MHT200cut";
    cutname[5]="RA2delphicut";
    cutname[6]="RA2noleptoncut";
    cutname[7]="noPhotoncut";
    cutname[8]="RA2Inc4Jetcut";
    cutname[9]="RA2Inc5Jetcut";
    cutname[10]="RA2Inc6Jetcut";
    cutname[11]="RA2allbutHT2500cut";
    cutname[12]="RA2allbutMHT1000cut";
    cutname[13]="RA2allcut";
    cutname[14]="RA2noleptoncutMHT1000";
    cutname[15]="RA2noleptoncutBtag2";
    cutname[16]="RA2noleptoncutBtag2MHT1000";
    cutname[17]="RA2Inc4JetcutMHT1000";
    cutname[18]="RA2Inc4JetcutBtag2";
    cutname[19]="RA2Inc4JetcutBtag2MHT1000";
    cutname[20]="RA2Inc5JetcutMHT1000";
    cutname[21]="RA2Inc5JetcutBtag2";
    cutname[22]="RA2Inc5JetcutBtag2MHT1000";
    cutname[23]="RA2Inc6JetcutMHT1000";
    cutname[24]="RA2Inc6JetcutBtag2";
    cutname[25]="RA2Inc6JetcutBtag2MHT1000";


    sigtype[0]="allEvents";
    sigtype[1]="glgl";

    BJtype[0]="allEvents";
    BJtype[1]="W";
    BJtype[2]="Wlv";
    BJtype[3]="Wjj";
    BJtype[4]="Z";
    BJtype[5]="Zll";
    BJtype[6]="Zvv";
    BJtype[7]="Zjj";
    BJtype[8]="photon";
    BJtype[9]="H";

   TTtype[0]="allEvents";
   TTtype[1]="TTbar";
   TTtype[2]="TTSingLep";
   TTtype[3]="TTdiLep";
   TTtype[4]="TThadronic";



  //KH
  histname[0]="weight";
  histname[1]="HT";
  histname[2]="MHT";
  histname[3]="NJet";
  histname[4]="NBtagLoose";
  histname[5]="NBtagTight";
  histname[6]="BtagLoose1Pt";
histname[7]="BtagLoose1Eta";
histname[8]="BtagLoose1Phi";
histname[9]="BtagLoose2Pt";
histname[10]="BtagLoose2Eta";
histname[11]="BtagLoose2Phi";
histname[12]="BtagTight1Pt";
histname[13]="BtagTight1Eta";
histname[14]="BtagTight1Phi";
histname[15]="BtagTight2Pt";
histname[16]="BtagTight2Eta";
histname[17]="BtagTight2Phi";

  ///end of initialization of the maps

yieldmap.clear();



//Signal Section//Signal Section//Signal Section//Signal Section//Signal Section//Signal Section//Signal Section//Signal Section

  //build a vector of scale factors
  //first load the cross sections into a vector
//Sig_xs_vec.push_back(0.757); /// v1
//Sig_xs_vec.push_back(1.12); // v2
//Sig_xs_vec.push_back(1.15); // v3
//Sig_xs_vec.push_back(1.14); // M(Stop,LSP)=(450,410) and also M(Stop,LSP)=(450,440)
//Sig_xs_vec.push_back(2.18); // M(Stop,LSP)=(400,390) and also M(Stop,LSP)=(400,360)
//Sig_xs_vec.push_back(4.41); // M(Stop,LSP)=(350,340) and also M(Stop,LSP)=(350,310)
//Sig_xs_vec.push_back(0.009635); //STOCv4
Sig_xs_vec.push_back(1.58); //StauC

  double Sig_numberofevents =0;//this will use GetSumOfWeights() 
  const int Sig_nHT = 1;   // Total number of HT bin samples
  const int nHist = 18; // Number of histograms in each TDirectory


  for(int i=1; i<=Sig_nHT ; i++){
//sprintf(tempname,"../Results/results_PhaseII4_Stop_CharmLSP_14TEV_140PileUp_00.root");
//.........这里部分代码省略.........

void ProduceDatacards_em(int nBins = 35,
			 float xmin = 0,
			 float xmax = 350,
			 bool pileUp = false) {

  double lumi = 2092;
  double normSS = 1.06;
    
  // sample names
  TString sampleNames[20] = {
    "MuEG_2015D", // data
    "DYJetsToLL_M-50_MG", // isZTT (ZTT)
    "DYJetsToLL_M-50_MG", // !isZTT (ZLL)
    "WJetsToLNu_MG", 
    "TTPowHeg",  
    "ST_t-channel_top_4f_leptonDecays",  
    "ST_t-channel_antitop_4f_leptonDecays",
    "ST_tW_antitop_5f_inclusiveDecays",
    "ST_tW_top_5f_inclusiveDecays",
    "VVTo2L2Nu",
    "WWToLNuQQ",
    "WZTo2L2Q",
    "WZTo1L1Nu2Q",
    "WZTo1L3Nu",
    "WZJets",
    "ZZTo4L",
    "ZZTo2L2Q",
    "",
    "",
    ""};

  double xsec[20] = {1, // data (0)
		     6025.2, // DY (1)
		     6025.2, // DY (2)
		     61526.7, // WJets (3)
		     831.76, // TT (4)
		     136.95, // ST_t-channel_top (5)
		     80.95,  // ST_t-channel_antitop (6)
		     35.6, // ST_tW_antitop (7)
		     35.6, // ST_tW_top_5f (8)
		     11.95,  // VV (9)
		     49.997, // WWToLNuQQ (10)
		     5.595, // WZTo2L2Q (11)
		     10.71, // WZTo1L1Nu2Q (12)
		     3.05, // WZTo1L3Nu (13)
		     5.26, // WZJets (3L1Nu) (14)
		     1.212, // ZZTo4L (15)
		     3.22, // ZZTo2L2Q (16)
		     0, //
		     0, //
		     0}; // 
  
  TString cuts[20];
  TString cutsSS[20];
  for (int i=0; i<20; ++i) {
    cuts[i] = "mcweight*(os>0.5)";
    cutsSS[i] = "mcweight*(os<0.5)";
  }

  cuts[0] = "os>0.5";
  cuts[1] = "mcweight*(os>0.5&&isZTT)";
  cuts[2] = "mcweight*(os>0.5&&!isZTT)";

  cutsSS[0] = "os<0.5";
  cutsSS[1] = "mcweight*(os<0.5&&isZTT)";
  cutsSS[2] = "mcweight*(os<0.5&&!isZTT)";
  
  if (pileUp) {
    for (int i=1; i<20; ++i) {
      cuts[i] = "puweight*" + cuts[i];
      cutsSS[i] = "puweight*" + cutsSS[i];
    }
  }

  TH1D * hist[20];
  TH1D * histSS[20];

  for (int i=0; i<17; ++i) {
    std::cout << sampleNames[i] << std::endl;
    TFile * file = new TFile(sampleNames[i]+".root");
    TH1D * histWeightsH = (TH1D*)file->Get("histWeightsH");
    TTree * tree = (TTree*)file->Get("TauCheck");
    double norm = xsec[i]*lumi/histWeightsH->GetSumOfWeights();
    TString histName = sampleNames[i] + "_mvis";
    TString histNameSS = sampleNames[i] + "_mvis_os";
    hist[i] = new TH1D(histName,"",nBins,xmin,xmax);
    histSS[i] = new TH1D(histNameSS,"",nBins,xmin,xmax);
    hist[i]->Sumw2();
    histSS[i]->Sumw2();
    tree->Draw("m_vis>>"+histName,cuts[i]);
    tree->Draw("m_vis>>"+histNameSS,cutsSS[i]);
    if (i>0) {
      for (int iB=1; iB<=nBins; ++iB) {
	double x = hist[i]->GetBinContent(iB);
	double e = hist[i]->GetBinError(iB);
    	hist[i]->SetBinContent(iB,norm*x);
    	hist[i]->SetBinError(iB,norm*e);
	double xSS = histSS[i]->GetBinContent(iB);
	double eSS = histSS[i]->GetBinError(iB);
    	histSS[i]->SetBinContent(iB,norm*xSS);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	else if (argv[2] == st2) {ChiMass=200;mIntermediate=500;}
	*/
      if (isData) XSec=1.;
	ChiMass=0.0;
    }

  if (XSec<0&& !isData) {cout<<" Something probably wrong with the xsecs...please check  - the input was "<<argv[2]<<endl;return 0;}

  std::vector<unsigned int> allRuns; allRuns.clear();

  cout<<" ChiMass is "<<ChiMass<<"  "<<mIntermediate<<endl;
  bool doThirdLeptVeto=true;
  bool doMuVeto=true;

  //CutList[CutNumb]=CutListt[CutNumb];
  char ff[100];

	
  sprintf(ff,"%s/%s",argv[3],argv[2]);

  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);

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

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


  vertexDataH->Scale(1/normVertexData);
  vertexMcH->Scale(1/normVertexMc);
 
  PileUp * PUofficial = new PileUp();

  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);


  TFile *f10= new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfDataBarrel);  // mu SF barrel data
  TFile *f11 = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfDataEndcap); // mu SF endcap data
  TFile *f12= new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfMcBarrel);  // mu SF barrel MC
  TFile *f13 = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/"+muonSfMcEndcap); // mu SF endcap MC 
  
  TGraphAsymmErrors *hEffBarrelData = (TGraphAsymmErrors*)f10->Get("ZMassBarrel");
  TGraphAsymmErrors *hEffEndcapData = (TGraphAsymmErrors*)f11->Get("ZMassEndcap");
  TGraphAsymmErrors *hEffBarrelMC = (TGraphAsymmErrors*)f12->Get("ZMassBarrel");
  TGraphAsymmErrors *hEffEndcapMC = (TGraphAsymmErrors*)f13->Get("ZMassEndcap");
  
  double * dataEffBarrel = new double[10];
  double * dataEffEndcap = new double[10];
  double * mcEffBarrel = new double[10];
  double * mcEffEndcap = new double[10];
  
  dataEffBarrel = hEffBarrelData->GetY();

//.........这里部分代码省略.........
  //  type[8]="photon";
//    type[9]="H";
   type[10]="TTbar";
   type[11]="TTSingLep";
   type[12]="TTdiLep";
   type[13]="TThadronic";



  //KH
  histname[0]="weight";
  histname[1]="METAsys";
  histname[2]="MET";
  histname[3]="NJet";
  histname[4]="j1Pt";
  histname[5]="Jet1Eta";
  histname[6]="Jet1Phi";
histname[7]="j2Pt";
histname[8]="Jet2Eta";
histname[9]="Jet2Phi";
histname[10]="j3Pt";
histname[11]="Jet3Eta";
histname[12]="Jet3Phi";
histname[13]="DelPhij1j2";
histname[14]="NLep";
histname[15]="NElec";
histname[16]="NMuon";
histname[17]="NTau";
  ///end of initialization of the maps

  //build a vector of scale factors
  //first load the cross sections into a vector
  xs_vec.push_back(530.89358);
  xs_vec.push_back(42.55351);
  xs_vec.push_back(4.48209);
  xs_vec.push_back(0.52795);
  xs_vec.push_back(0.05449);

  double numberofevents =0;
  const int ttnHT = 5;   // Total number of HT bin samples
  const int nHist = 18; // Number of histograms in each TDirectory


  for(int i=1; i<=ttnHT ; i++){
    sprintf(tempname,"../Results/results_PhaseI_TT_14TEV_HT%d_NoPileUp_00.root",i);
    file = new TFile(tempname, "R");
    sprintf(tempname,"allEvents/nocut/MET_nocut_allEvents");
    tempvalue = (luminosity*xs_vec[i-1])/((* (TH1D* ) file->Get(tempname)).GetEntries());
    scalevec.push_back(tempvalue);
  }//end of loop over HTbins 
  std::cout << "normalization scale factor determination done" << std::endl;
for(int i=1; i<=ttnHT; i++){
sprintf(tempname,"../Results/results_PhaseI_TT_14TEV_HT%d_NoPileUp_00.root",i);
inputfilevec.push_back(TFile::Open(tempname,"R"));
}

tempstack = new THStack("stack","Binned Sample Stack");
file = new TFile("stack.root","RECREATE");
 for(map<int , string >::iterator itt=type.begin(); itt!=type.end();itt++){        // loop over different event types
    cdtoitt = file->mkdir((itt->second).c_str());
    cdtoitt->cd();
    for(map<int , string >::iterator it=cutname.begin(); it!=cutname.end();it++){   // loop over different cutnames
      cdtoit =  cdtoitt->mkdir((it->second).c_str());
      cdtoit->cd();
      for(int j=0; j<histname.size(); j++){                                        // loop over different histograms
        for(int i=0; i<ttnHT ; i++){                                                  // loop over different HT bins
sprintf(tempname,"%s/%s/%s_%s_%s",(itt->second).c_str(),(it->second).c_str(),(histname[j]).c_str(),(it->second).c_str(),(itt->second).c_str());
temphist = (TH1D *) inputfilevec.at(i)->Get(tempname)->Clone();
temphist->Scale(scalevec[i]);
if(histname[j]=="MET"){numberofevents+=(double)temphist->GetSumOfWeights();} //all the histograms in one directory have the same number of events
//if(histname[j]=="MET"){cout << " temphist->GetSumOfWeights() " << temphist->GetSumOfWeights() << endl;}
/*if(i==0){
cout << "" << endl;
cout << "type: " << (itt->second).c_str() << ",  cutname: " << (it->second).c_str()<< ", histname: " << histname[j].c_str() << ", bin#: " << i << endl;
cout << "temphist->GetEntries():  " << temphist->GetEntries() << endl;
cout << "temphist->GetSumOfWeights(): " << temphist->GetSumOfWeights() << endl;
cout << " ===============================================================" << endl;
}
*/
temphist->SetFillColor(i+2);
          tempstack->Add(temphist);
               }//end of loop over HTbins 1..7
if(histname[j]=="MET"){
cout << " +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
cout << "type: " << (itt->second).c_str() << ",  cutname: " << (it->second).c_str() << endl;
cout << "Number of events:  " << numberofevents << endl;
cout << " +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
}
numberofevents=0;
        sprintf(tempname,"%s_%s_%s",histname[j].c_str(),(it->second).c_str(),(itt->second).c_str());
        tempstack->Write(tempname);
        delete tempstack;
        tempstack = new THStack("stack","Binned Sample Stack");
      }//end of loop over histograms
    }//end of loop over cutnames
  }//end of loop over event types
file->Close();


}//end of the constructor