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

void QinvKstar (const char* filename, const char* histname) {

	TFile*  fileIn = new TFile(filename,"update");
	TH1D* hist = (TH1D*)fileIn->Get(histname);
	TH1D* histnew = new TH1D(Form("kstar_%s",histname),"",hist->GetNbinsX(),0,0.25);//(TH1D*)fileIn->Get(histname);

//   for (int ibins = 1; ibins <= hist->GetNbinsX()*2; ibins++ ) {
//     hist->SetBinContent(ibins, hist->GetBinContent(ibins));
//     //    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
//   }

	int bin=2;//1

	for (double kstar = 0.00375; kstar <= 0.25; kstar += 0.0025 ) {

//	for (double kstar = 0.00375; kstar <= 0.25; kstar += 0.0025 ) {
		//  for (double kstar = 0.005; kstar <= 0.25; kstar += 0.01 ) {
		//        hist->Fill(kstar,hist->GetBinContent(hist->FindBin(kstar*2)));

		histnew->SetBinContent(bin,hist->GetBinContent(hist->FindFixBin(kstar*2)));
		histnew->SetBinError(bin++,hist->GetBinError(hist->FindFixBin(kstar*2)));

		//    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
		//    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
	}

	TFile* fileOut = new TFile(Form("kstar_%s",filename), "update");
	hist->Write();
	histnew->Write();
}

// this is the event loop
void BeamGas::loop() {

  // DO ANY SETUP HERE
  nTracks_ = new TH1D("hNTracks", "N_{trk}", 100, 0., 100.);
  nTracksHQ_ = new TH1D("hNTracksHQ", "N_{trkHQ}", 100, 0., 100.);

  nTracksNoVtx_ = new TH1D("hNTracksNoVtx", "N_{trk}", 100, 0., 100.);
  nTracksHQNoVtx_ = new TH1D("hNTracksHQNoVtx", "N_{trkHQ}", 100, 0., 100.);

  trackPFraction_ = new TH1D("hTrackPFraction", "", 100, 0., 2.);

  reset();
  nextEvent();

  for (unsigned long i=0; i<maxEvents_; ++i, nextEvent()) {

    // occasional print out
    if (i%100000==0) {
      std::cout << "Processing " << i << "th event of " <<maxEvents_<< std::endl;
    }

    // YOUR CODE HERE
    if (event_->jet_N>0 && event_->jetE[0]>50. && event_->jetEta[0]<1.3) {

      // count tracks of different types
      unsigned nTracks(0), nTracksHQ(0);
      for (unsigned j=0; j<event_->track_N; ++j) {
	nTracks++;
      	if (event_->trackQuality[0] == 1) nTracksHQ++;
      }
      
      nTracks_->Fill(nTracks);
      nTracksHQ_->Fill(nTracksHQ);

      if (event_->vtx_N == 0) {
	nTracksNoVtx_->Fill(nTracks);
	nTracksHQNoVtx_->Fill(nTracksHQ);
      }

      // leading jet/track ratio
      if (event_->track_N>0) {
	trackPFraction_->Fill(event_->trackP[0]/event_->jetE[0]);
      }

    }

  }


  // SAVE HISTOGRAMS HERE
  nTracks_->Write("",TObject::kOverwrite);
  nTracksHQ_->Write("",TObject::kOverwrite);
  nTracksNoVtx_->Write("",TObject::kOverwrite);
  nTracksHQNoVtx_->Write("",TObject::kOverwrite);
  trackPFraction_->Write("",TObject::kOverwrite);
  
}

void takeDirectlyFromMC(TFile* fin, TFile* fout, TString gentype) {

  TList* listOfDirs = fin->GetListOfKeys();
  for (auto k : *listOfDirs) {
    TString srname = k->GetName();
    if (!srname.Contains("sr")) continue;
    if (srname.Contains("base") || srname.Contains("incl") || srname.Contains("sb")) continue;
    if (gentype == "_1lepW" && !srname.EndsWith("2") && !srname.EndsWith("3")) continue;

    auto indir = (TDirectoryFile*) fin->Get(srname);
    auto outdir = (TDirectory*) fout->mkdir(srname);
    auto hlist = indir->GetListOfKeys();

    for (auto h : *hlist) {
      TString hname = h->GetName();
      if (!hname.BeginsWith("h_metbins" + gentype)) continue;
      TH1D* hin = (TH1D*) indir->Get(hname);
      outdir->cd();
      TH1D* hout = (TH1D*) hin->Clone(TString(hname).ReplaceAll(gentype, ""));
      for (int i = 1; i <= hout->GetNbinsX(); ++i) {
        // zero out negative yields
        if (hout->GetBinContent(i) < 0) {
          hout->SetBinContent(i, 0);
          hout->SetBinError(i, 0);
        }
      }
      hout->Write();

      if (yearSeparateSyst && (hname.EndsWith("Up") || hname.EndsWith("Dn"))) {
        for (int i = 1; i < 4; ++i) {
          TH1D* hcen_yi = (TH1D*) fbkgs[i]->Get(srname+"/h_metbins"+gentype);
          TH1D* hsys_yi = (TH1D*) fbkgs[i]->Get(srname+"/"+hname);
          if (hsys_yi && !hcen_yi) {
            cout << "Find " << srname+"/"+hname << " from " << fbkgs[i]->GetName() << " but not hcen " << srname+"/h_metbins"+gentype << " Should not happen?" << endl;
          }
          TH1D* hout_yi = (TH1D*) fin->Get(srname+"/h_metbins"+gentype)->Clone(TString(hname).Insert(hname.Length()-2, Form("%d", 15+i)).ReplaceAll(gentype, ""));
          if (hcen_yi) hout_yi->Add(hcen_yi, -1);
          if (hsys_yi) hout_yi->Add(hsys_yi);
          hout_yi->Write();
        }
      }
    }
    if (!outdir->Get("h_metbins")) {
      cout << "Didn't find yield hist for " << gentype << " in " << fin->GetName() << ":" << srname << "/. Faking a 0 one!" << endl;
      outdir->cd();
      // Get the MET binning from h_metbins, which shall always exist, and then set all bins to 0
      TH1D* hout = (TH1D*) fin->Get(srname + "/h_metbins")->Clone("h_metbins");
      for (int i = 1; i <= hout->GetNbinsX(); ++i) {
        hout->SetBinContent(i, 0);
        hout->SetBinError(i, 0);
      }
      hout->Write();
    }
  }
}

int produceEnergyHistos(CSPrereqs& csp)
{
    cout << endl << "Mapping " << csp.name << " TOF histogram to energy domain..." << endl;

    // open input file
    TFile* inputFile = new TFile(inputFileName.c_str(),"READ");
    if(!inputFile->IsOpen())
    {
        cerr << "Error: failed to open " << inputFileName << "  to fill histos." << endl;
        return 1;
    }

    // create output file
    TFile* outputFile = new TFile(outputFileName.c_str(),"RECREATE");

    for(auto& channelName : config.cs.DETECTOR_NAMES)
    {
        TDirectory* detectorDirectory = inputFile->GetDirectory(channelName.c_str());
        if(!detectorDirectory)
        {
            cerr << "Error: failed to find " << channelName << " directory in " << inputFileName << "." << endl;
            return 1;
        }

        TDirectory* outputDirectory = outputFile->mkdir(channelName.c_str(),channelName.c_str());

        // find TOF histograms in input file
        for(int i=0; i<config.target.TARGET_ORDER.size(); i++)
        {
            string targetName = config.target.TARGET_ORDER[i];
            string TOFHistoName = targetName + "TOFCorrected";

            detectorDirectory->cd();

            TH1D* TOF = (TH1D*)detectorDirectory->Get(TOFHistoName.c_str());
            if(!TOF)
            {
                cerr << "Error: failed to open " << TOFHistoName << " in " << inputFileName << endl;
                return 1;
            }

            outputDirectory->cd();

            TH1D* correctedEnergy = convertTOFtoEnergy(TOF, targetName + "Energy");
            TOF->Write();
            correctedEnergy->Write();
        }
    }

    outputFile->Close();
    inputFile->Close();

    return 0;
}

// Takes as input a given process with corresponding luminosity, as well as sample CME, cuts, directory location, and output file name - saves corresponding histogram in root file.
void processtree(TString sample, Double_t weight, TLorentzVector CME, TString filename, TString directory, TString param, TString cut, std::clock_t start, Int_t all){

    //Open File
    TFile *f = new TFile(directory);
    TTree *ttree = (TTree*)f->Get(sample);

    TString recohistname, paramhistname;
    if (all==0){
        recohistname = "_recoilmassHIST";
        paramhistname = "_parameterHIST";
    }
    if (all==1){
        recohistname = "_recoilmassHIST_all";
        paramhistname = "_parameterHIST_all";
    }

    TString histname = sample + recohistname;
    TH1D *recoilmassHIST   = new TH1D(histname,histname,CME.M()/2,0,CME.M());
    //recoilmassHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, "recoilmass", cut);

    histname = sample + paramhistname;
    TH1D *parameterHIST   = new TH1D(histname,histname,CME.M()/2,1,-1);
    //parameterHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, param, cut);

    //Normalize Hist
    recoilmassHIST->Scale(1/recoilmassHIST->GetMaximum());
    parameterHIST->Scale(1/parameterHIST->GetMaximum()) ;

    //Normalize Hist
//    recoilmassHIST->Scale(weight);
//    parameterHIST->Scale(weight) ;

    // Writes yield of each histogram for each cut to text file
    ofstream myfile;
    Double_t yield = recoilmassHIST->Integral(); 
    myfile.open ("yield.txt", ios::app);
        cout << "Yield for " << sample << " for Cut " << param << " is " << yield << endl;
        myfile << param << " " << all << " " << sample << " " << yield << "\n"; //write to file
    myfile.close();

    if (recoilmassHIST->Integral() == 0){
        recoilmassHIST->Fill(0);
    }
    if (parameterHIST->Integral() == 0){
        parameterHIST->Fill(0);
    }

    TFile g(filename, "update");
    recoilmassHIST->Write();
    parameterHIST->Write();
    g.Close();
}

void photonSpectra()
{
  TH1::SetDefaultSumw2();

  const Double_t PTBINS[] = {40, 50, 60, 80, 120, 300};
  const Int_t nPTBINS = 5;

  const TString DATA_FILE = "gammaJets_pA_Data.root";
  const TString MC_FILE = "gammaJets_pA_MC_allQCDPhoton.root";
  TFile *dataFile = TFile::Open(DATA_FILE);
  TTree *dataTree = (TTree*)dataFile->Get("photonTree");
  TFile *mcFile = TFile::Open(MC_FILE);
  TTree *mcTree = (TTree*)mcFile->Get("photonTree");

  TH1D *photonSpectra = new TH1D("photonSpectra","photonSpectra", nPTBINS, PTBINS);
  TH1D *correctedPhotonSpectra = (TH1D*)photonSpectra->Clone("correctedPhotonSpectra");

  const TCut photonEtaCut = "abs(photonTree.eta) < 1.44";
  const TCut isoCut = "ecalRecHitSumEtConeDR04 < 4.2  &&  hcalTowerSumEtConeDR04 < 2.2  &&  trkSumPtHollowConeDR04 < 2 && hadronicOverEm<0.1";
  const TCut genCut = "genCalIsoDR04<5 && abs(genMomId)<=22";
  const TCut sbCut = "(cc4+cr4+ct4PtCut20>10) && (cc4+cr4+ct4PtCut20<20) && hadronicOverEm<0.1";
  const TCut candidateCut = "sigmaIetaIeta<0.01";
  const TCut decayCut = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";

  TCut photonPtCut = "photonTree.corrPt>40 && photonTree.corrPt<300";

  dataTree->Project(photonSpectra->GetName(),"corrPt",isoCut && photonEtaCut && photonPtCut && candidateCut);

  for(int i = 0; i<nPTBINS; ++i)
  {
    photonPtCut = Form("photonTree.corrPt>%f && photonTree.corrPt<%f", PTBINS[i], PTBINS[i+1]  );

    TCut dataCandidateCut = isoCut && photonEtaCut && photonPtCut;
    TCut sidebandCut =  sbCut && photonEtaCut && photonPtCut;
    TCut mcSignalCut = dataCandidateCut && genCut;

    fitResult fitr = getPurity(dataTree, mcTree,
  			       dataCandidateCut, sidebandCut,
  			       mcSignalCut);
    Double_t purity = fitr.purity;
    correctedPhotonSpectra->SetBinContent(i+1,photonSpectra->GetBinContent(i+1)*purity);
  }

  // TCanvas *c1 = new TCanvas();
  // photonSpectra->Draw();
  // correctedPhotonSpectra->Draw("same");

  TFile *outFile = new TFile("photonSpectra.root", "RECREATE");
  outFile->cd();
  photonSpectra->Write();
  correctedPhotonSpectra->Write();
  outFile->Close();
}

void savePlots2(const std::string& canvasName, TH1D& dataPlot, TH1D& fitPlot,
		TH1D& signalPlot, TH1D& bkgPlot) {
  TCanvas c(canvasName.c_str(), canvasName.c_str()) ;
  TPad mainPad("mainPad", "mainPad", 0., 0.2, 1., 1.) ;
  mainPad.Draw() ;
  TPad pullPad("pullPad", "pullPad", 0., 0., 1., 0.2) ;
  pullPad.Draw() ;

  mainPad.cd() ;
  dataPlot.SetStats(true) ;
  //dataPlot.SetOptStat(111111111);
  dataPlot.SetLineWidth(2) ;
  dataPlot.Draw() ;
  fitPlot.SetLineWidth(2) ;
  fitPlot.SetLineColor(kBlue) ;
  fitPlot.Draw("same") ;
  signalPlot.SetLineWidth(2) ;
  signalPlot.SetLineColor(kGreen) ;
  signalPlot.Draw("same") ;
  bkgPlot.SetLineColor(kRed) ;
  bkgPlot.SetLineWidth(2) ;
  bkgPlot.Draw("same") ;

  TLegend leg(0.6, 0.7, 0.9, 0.9) ;
  leg.AddEntry(&dataPlot, "Data") ;
  leg.AddEntry(&fitPlot, "Fit") ;
  leg.AddEntry(&signalPlot, "Signal") ;
  leg.AddEntry(&bkgPlot, "Background") ;
  leg.SetFillStyle(0) ;
  leg.SetBorderSize(0) ;
  leg.Draw() ;

  pullPad.cd() ;
  pullPad.SetGridy() ;
  TH1D* h_residuals(NULL), *h_pulls(NULL), *h_pullDistribution(NULL) ;
  FitterTools::makePullPlot(dataPlot, fitPlot, h_residuals, h_pulls, h_pullDistribution) ;
  h_pulls->Draw() ;

  c.Write() ;

  h_pulls->Write() ;
  delete h_pulls ;
  h_residuals->Write() ;
  delete h_residuals ;
  h_pullDistribution->Write() ;
  delete h_pullDistribution ;

  dataPlot.Write() ;
  fitPlot.Write() ;
}

void MCefficiency(TString inputmc, TString outputfile, TString tfend, TString selmcgen, TString cut, TString weight, Float_t centmin, Float_t centmax)
{
  Float_t hiBinMin,hiBinMax;
  hiBinMin = centmin*2;
  hiBinMax = centmax*2;

  selmcgen = Form("%s&&hiBin>=%f&&hiBin<%f",selmcgen.Data(),hiBinMin,hiBinMax);
  cut = Form("%s&&hiBin>=%f&&hiBin<%f",cut.Data(),hiBinMin,hiBinMax);

  TFile* infMC = new TFile(inputmc.Data());
  TTree* ntMC = (TTree*)infMC->Get("ntDkpi");
  ntMC->AddFriend("ntHi");
  ntMC->AddFriend("ntSkim");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  ntGen->AddFriend("ntHi");
  ntGen->AddFriend("ntSkim");

  TH1D* hPtMC = new TH1D("hPtMC","",nPtBins,ptBins);
  TH1D* hPtGen = new TH1D("hPtGen","",nPtBins,ptBins);

  ntMC->Project("hPtMC","Dpt",TCut(weight)*(TCut(cut.Data())&&"(Dgen==23333)"));
  divideBinWidth(hPtMC);
  ntGen->Project("hPtGen","Gpt",TCut(weight)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);
  TH1D* hEff = (TH1D*)hPtMC->Clone("hEff");
  hEff->Divide(hPtGen);

  TFile* fout=new TFile(Form("%s_cent_%.0f_%.0f_%s.root",outputfile.Data(),centmin,centmax,tfend.Data()),"recreate");
  fout->cd();
  hEff->Write();
  fout->Close();
}

void
getLimits (TH1D *cutFlow, TDirectoryFile *dir)
{
  string histName = cutFlow->GetName ();
  TH1D *lowerLimit = (TH1D *) cutFlow->Clone ((histName + "LowerLimit").c_str ()),
       *upperLimit = (TH1D *) cutFlow->Clone ((histName + "UpperLimit").c_str ());

  TAxis *x = cutFlow->GetXaxis ();
  for (int i = 1; i <= x->GetNbins (); i++)
    {
      double events = cutFlow->GetBinContent (i), eventsLowerLimit, eventsUpperLimit;

      // The calculation of upper and lower limits is taken from the PDG Statistics chapter.  
      // "Poisson or binomal data", section 37.4.2.3 in 2013 version; Eqns. 37.71a, 37.71b.  
      // Here ALPHA is the confidence level, so ALPHA = 0.68 corresponds to a 68% C.L.
      // The PDG uses a different convention:  the C.L is 1 - \alpha.  
      // For example, the 68% CL upper limit on an observation of 0 events is:
      // 0.5 * TMath::ChisquareQuantile (0.68, 2) = 1.139  
      eventsLowerLimit = 0.5 * TMath::ChisquareQuantile (1 - ALPHA, 2 * events);
      eventsUpperLimit = 0.5 * TMath::ChisquareQuantile (ALPHA, 2 * (events + 1));
      lowerLimit->SetBinContent (i, eventsLowerLimit);
      upperLimit->SetBinContent (i, eventsUpperLimit);
    }
  dir->cd ();
  lowerLimit->Write ((histName + "LowerLimit").c_str ());
  upperLimit->Write ((histName + "UpperLimit").c_str ());
}

void Pileup_reweight_gen() {

   TFile *f1 = new TFile("MyDataPileupHistogram_true.root");
   TH1D *hPU_data = (TH1D*)f1->Get("pileup");  

   TFile *f2 = new TFile("PU_MC_histo.root");
   TH1D *hPU_MC = (TH1D*)f2->Get("PU_MC");

   hPU_data->Scale(1.0/hPU_data->Integral());
   hPU_MC->Scale(1.0/hPU_MC->Integral());

   TH1D *pileupWeights = new TH1D("pileupWeights", "Computed Weights", 52, 0., 52);


  std::vector<double> weights(hPU_data->GetNbinsX()+1, 1.0);

  for (int i = 1, n = weights.size(); i < n; ++i) {
        double nMC = hPU_MC->GetBinContent(i), nData = hPU_data->GetBinContent(i);
        weights[i-1] = (nMC > 0 ? nData/nMC : 1.0);
        pileupWeights->SetBinContent(i,weights[i-1]);
    }

  
     pileupWeights->Draw(); 

   TFile *outFile = new TFile("pileupWeights.root", "RECREATE");
   outFile->cd();
   pileupWeights->Write();
   outFile->Close();



}

void rebin (const char* filename, const char* histname, double rebin=2.0) {
        TFile* f = new TFile(filename,"update");
        TH1D* h = (TH1D*)f->Get(histname);
        h->Rebin(rebin);
        h->Scale(1.0/rebin);
        TFile* fileOut = new TFile(Form("rebin_%s",filename), "update");
        h->Write();
}

void merge(){
	TVectorD Nevent;	Nevent.ResizeTo(nbin);  Nevent.Zero();
     //   TVectorD totmultall;	totmultall.ResizeTo(nbin);      totmultall.Zero();
     //   TVectorD avgmultall;	avgmultall.ResizeTo(nbin);      avgmultall.Zero();
        TVectorD tottrk;	tottrk.ResizeTo(nbin);      tottrk.Zero();
//	TVectorD totptall;      totptall.ResizeTo(nbin);    totptall.Zero();
//	TVectorD totetaall;      totetaall.ResizeTo(nbin);    totetaall.Zero();
        TVectorD avgtrk;	avgtrk.ResizeTo(nbin);      avgtrk.Zero();
//	TVectorD avgmult;       avgmult.ResizeTo(nbin);
            
        TH2F* s[nbin];
        TH2F* b[nbin];
        for(int ibin=0;ibin<nbin;ibin++){
            s[ibin] = new TH2F(Form("s_%d",ibin),Form("signal",ibin),detastep,detamin,detamax,dphistep,dphimin,dphimax);
            s[ibin]->Sumw2();
            b[ibin]  = new TH2F(Form("b_%d",ibin), "background",detastep,detamin,detamax,dphistep,dphimin,dphimax);
            b[ibin]->Sumw2();
        }
        TFile *fout = new TFile(Form("Anav3_merged.root"),"Recreate");
        TFile *f[nFileAll];
        for(int ifile=; ifile<63; ifile++){
                f[ifile] = TFile::Open(Form("%s/Anav3_%d.root",outdir.Data(),ifile));
                if(!f[ifile]) continue;
		TVectorD* Nevent_t =  (TVectorD*)f[ifile]->Get(Form("Nevent"));
	//	TVectorD* totmultall_t =  (TVectorD*)f[ifile]->Get(Form("totmultall"));
		TVectorD* tottrk_t =  (TVectorD*)f[ifile]->Get(Form("tottrk"));
	//	TVectorD* totptall_t =  (TVectorD*)f[ifile]->Get(Form("totptall"));
//		TVectorD* totetaall_t =  (TVectorD*)f[ifile]->Get(Form("totetaall"));
		for(int ibin=0;ibin<nbin;ibin++){
			//	totptall[ibin] += (*totptall_t)[ibin];
			//	totetaall[ibin] += (*totetaall_t)[ibin];
			        Nevent[ibin] += (*Nevent_t)[ibin];
			  //      totmultall[ibin] += (*totmultall_t)[ibin];	
			        tottrk[ibin] += (*tottrk_t)[ibin];	
                                TH1D* s_t = (TH1D*)f[ifile]->Get(Form("signal_%d",ibin));
                                TH1D* b_t = (TH1D*)f[ifile]->Get(Form("background_%d",ibin));
                                s[ibin]->Add(s_t);
                                b[ibin]->Add(b_t);
		}
		f[ifile]->Close();
        }
    for(int ibin=0;ibin<nbin;ibin++){
    avgtrk[ibin] = tottrk[ibin]/Nevent[ibin];
    TH1D* hr = (TH1D*)s[ibin]->Clone("hr");
    hr->Divide(b[ibin]);
    hr->Scale(b[ibin]->GetBinContent(b[ibin]->FindBin(0,0)));

    fout->cd();
    Nevent.Write("Nevent");
    avgtrk.Write("avgtrk");
    TDirectory *dir0 = (TDirectory*)fout->mkdir(Form("D_%d",ibin));
    dir0->cd();
    s[ibin]->Write("s");
    b[ibin]->Write("b");
    hr->Write();
    }
}

int main(){
	TFile * data = TFile::Open("TreesMu_Data_plots.root");
	TH1D * dataMtop = (TH1D*)data->Get("antiEtaFwDTrue_allW/antiEtaFwDTrue_allWcosTheta");
	
	TFile * tt = TFile::Open("TreesMu_TTBar_RW.root");
	
	TH1D * ttMtop = (TH1D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta");
	
	TH2D * ttMtop2D = (TH2D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta2D");
	
  std::pair<TF1, WeightFunctionCreator*> WeightFuncUD(WeightFunctionCreator::getWeightFunction("WeightFuncUDF", F0, FL));
  std::pair<TF1, WeightFunctionCreator*> WeightFuncDU(WeightFunctionCreator::getWeightFunction("WeightFuncDUF", F0, FL));
  cout<<F0 + F0Sys<<"\t"<<FL - FLSys<<endl;
	WeightFuncUD.first.SetParameters(F0 + F0Sys, FL - FLSys);
	WeightFuncDU.first.SetParameters(F0 - F0Sys, FL + FLSys);
	
	TH1D * ttUD = myReweightor(ttMtop2D,WeightFuncUD,"WeightFuncUD_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttUD->Integral()<<endl;
	ttUD->Scale(ttMtop2D->Integral()/ttUD->Integral());
	ttUD->Add(ttMtop);
	TH1D * wUD = (TH1D*)ttUD->Clone("wUD");
	wUD->Scale(-1.);
	wUD->Add(dataMtop);
	
	TH1D * ttDU = myReweightor(ttMtop2D,WeightFuncDU,"WeightFuncDU_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttDU->Integral()<<endl;
	ttDU->Scale(ttMtop2D->Integral()/ttDU->Integral());
	ttDU->Add(ttMtop);
	TH1D * wDU = (TH1D*)ttDU->Clone("wDU");
	wDU->Scale(-1.);
	wDU->Add(dataMtop);
	
	TFile * out = new TFile("file.root","recreate");
	out->cd();
	WeightFuncUD.first.Write();
	WeightFuncDU.first.Write();
	ttUD->Write();
	ttDU->Write();	
	wUD->Write();
	wDU->Write();	
	out->Close();
	return 1;
}

void Run_ZgSelectorMC(Int_t eleID, Int_t phoID,Float_t DelRCut){

   gROOT->LoadMacro("Zg_SelectorMC.C+");

   FILE *Table = fopen("Input_table_Sig", "r");
   char filename[200];
   char tmp[100];
   int  flag = 1;
   int  nfile(0);
   double EvtWeight[40];
   double nevt[40];
   double luminosity(1.0);
   TString ProcessTag[40];

   //Load Input table 
   cout<<"<<< File Table >>>"<<endl;
   while (flag != -1){
      flag=fscanf(Table, "%s", tmp);
      if (flag == -1) continue;
      ProcessTag[nfile] = tmp;
      flag=fscanf(Table, "%s", filename);
      flag=fscanf(Table, "%s", tmp);
      double cross=atof(tmp);
      flag=fscanf(Table, "%s", tmp);
      double evt=atof(tmp);

      EvtWeight[nfile] = 0;
      EvtWeight[nfile] = cross/evt*luminosity;
      if (nfile == 0) nevt[nfile] = evt;
      else nevt[nfile] = evt + nevt[nfile-1];
      cout<<">>> "<<nfile<<" : "<<filename<<endl;
      cout<<"    "<<ProcessTag[nfile]<<"  "<<cross<<"  "<<evt<<"  "<<nevt[nfile]<<"  "<<EvtWeight[nfile]<<endl;
      nfile += 1;
   }
   fclose(Table);

   cout<<"  "<<endl;
   TString TSaveFileName = Form("Test_eID%d_pID%d_dR07",eleID,phoID);
   TH1D *hInfo  = new TH1D("hInfo","Event Information",15,0,15);
   hInfo->SetBinContent(1,nfile);
   hInfo->SetBinContent(2,eleID);
   hInfo->SetBinContent(3,phoID);

   TFile *f1 = new TFile(TSaveFileName+".root","RECREATE");
   f1->cd();
   hInfo->Write();
   f1->Close();

   for (Int_t iSample=0;iSample<nfile;iSample++) {
     cout<<"Now Running....sample:"<<iSample<<", Weight:"<<EvtWeight[iSample]<<endl;
     Zg_SelectorMC t(iSample);
     t.Loop(iSample,phoID,eleID,EvtWeight[iSample],ProcessTag[iSample],TSaveFileName,DelRCut);
   }
   cout<<"Done!"<<endl;
}

void fixLeptonBJet(char* filename)
{
    TH1D *lb;
    TFile f(filename, "UPDATE");
    f.GetObject("HypLeptonBAjetMass", lb);
    if (lb) {
        cout << "HypLeptonBAjetMass found, saving as HypLeptonBjetMass\n";
        lb->SetName("HypLeptonBjetMass");
        lb->Write();
    }
    f.Close();
}