C++ TH2D::Clone方法代码示例

TH1* getHisto(char * filename, char* histoName, char * dirName, int nBin, double lumi)
{
  TH1 * hpt_=0;
  TFile *file0 = TFile::Open(filename);
  if(!file0) return hpt_;
  TDirectory *dir;
  TH2D * hMuPt;

  if(dirName == "0") {
    hMuPt = (TH2D*) file0->Get(histoName);
  } else {
    dir = (TDirectory*) file0->Get(dirName);
    if(!dir) return hpt_;
    hMuPt = (TH2D*) dir->Get(histoName);
  }

  if(hMuPt) {
  hpt_ = (TH1*) hMuPt->Clone();
  hpt_->Sumw2();
  hpt_->Scale(1./lumi); // this take into into account the luminosity
  hpt_->SetLineWidth(2);
  hpt_->Rebin(nBin);
  double nBinX=hpt_->GetNbinsX();

  // overFlow
  hpt_->SetBinContent((int)nBinX,hpt_->GetBinContent((int)nBinX)+hpt_->GetBinContent((int)nBinX+1));
  hpt_->SetDirectory(0);
  file0->Close();
  hpt_->SetLineWidth(3);
  }
  return hpt_;
}

void runcorr(int filenum)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/oldstuff/sortedforests.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }


  corrana(listoffiles[filenum].data());
  
  double pttriglow[] = {1,4};
  double pttrighigh[] = {2,8};
  double ptasslow[] = {1,2};
  double ptasshigh[] = {2,3};
  int centmin[] = {0,4,8,12,16,20,24,28,32};
  int centmax[] = {40,8,12,16,20,24,28,32,36};
  TFile * outf = new TFile(Form("corrhists_%d.root",filenum),"recreate");
  
  // for(int i = 0 ; i < 9 ; ++i)
  // for(int i = 0 ; i < 2 ; ++i)
  for(int i = 0 ; i < 1 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    // for(int cent = 0 ; cent < 3 ; ++cent)
    // for(int cent = 0 ; cent < 9 ; ++cent)
    for(int cent = 0 ; cent < 1 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      TH2D * ttsig = TrackTrackSignal(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      TH2D * ttbak = TrackTrackBackground(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      // TCanvas * c1 = new TCanvas();
      // ttsig->Draw("surf1");
      // TCanvas * c2 = new TCanvas();
      // ttbak->Draw("surf1");
      // TCanvas * c3 = new TCanvas();
      TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[i],(int)pttrighigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      ttcorr->Divide(ttbak);
      ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(0,0)));
      ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(-4.0),ttcorr->GetXaxis()->FindBin(4.0));
      ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      // ttcorr->Draw("surf1");
    }
  }
  
  outf->Write();
  outf->Close();
}

TH2D* GetJetCorrFunc2D_ZYAM(int itrg, int jass)
{
                    TH2D* hcorr = (TH2D*)GetRawCorrFunc2D_ratio(itrg,jass);
                    TH2D* hcorr_clone = (TH2D*)hcorr->Clone(Form("corr_clone_itrg%d_jass%d",itrg,jass));
                    hcorr_clone->Reset();
		    for(int ietabin=1;ietabin<=hcorr->GetNbinsX();ietabin++)
		    {
                      TH1D* hcorrphi = (TH1D*)hcorr->ProjectionY(Form("corrphi_%d",ietabin),ietabin,ietabin,"e");  
                     // hcorrphi->SetAxisRange(0.35,1.57,"X");
                      double histminX = hcorrphi->GetBinCenter(hcorrphi->GetMinimumBin());
                      double histminY = hcorrphi->GetBinContent(hcorrphi->GetMinimumBin());
	   	      fitfunc = new TF1(Form("fitfunc_%d",ietabin),"[0]+[1]*(x-[2])*(x-[2])",0.6,1.2);		    
                      fitfunc->SetParameters(histminY,1,histminX);
                      fitfunc->SetParLimits(1,0,0.10000);
                  //    fitfunc->SetParLimits(2,0.35,1.57);
		      for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%d",ietabin),"RNO");                     		  
		      float level = fitfunc->GetParameter(0);
		      for(int iphibin=1;iphibin<=hcorr->GetNbinsY();iphibin++) 
                        hcorr_clone->SetBinContent(ietabin,iphibin,hcorr->GetBinContent(ietabin,iphibin)-level);
		    }
                    
                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");		    
                    return hcorr_clone;
}

//loads phi-pT efficiency from root file
TH2D *GetEfficiencyPhi(TString effLoc){
  TFile *fIn = 0;
  TH2D *hPhiPt = 0; 

  if(!fIn)fIn = TFile::Open(effLoc.Data());
  if(!fIn)Printf("%s%d no input data",(char*)__FILE__,__LINE__);
  if(!hPhiPt)hPhiPt = (TH2D*)fIn->Get("h2TrackPtPhiNorm");
  if(!hPhiPt) cout<<"Could not load h2TrackPtPhiNorm"<<endl; 
  if(!hPhiPt)Printf("%s%d no phi-pt efficiency spectrum available",(char*)__FILE__,__LINE__);

  gROOT->cd();
  TH2D *hPhiPtClone = (TH2D*)hPhiPt->Clone(hPhiPt->GetName()); 
  fIn->Close();
  return hPhiPtClone;

}

TH2D* Plot2D (std::string var,std::vector<TFile*>& tfiles,std::vector<double>& weights) {

    double weightTot = 0.0;
    for ( int i = 0; i<weights.size(); i++) {
        weightTot+=weights[i];
    }

    std::string name = var;
    name+="NEW";
    TH2D* hVar = (TH2D*)tfiles[0]->Get(var.c_str());
    TH2D* HVar = (TH2D*)hVar->Clone(name.c_str());
    //HVar->Sumw2();

    for (int i=1; i<tfiles.size(); i++) {
        TH2D* htempVar = (TH2D*)tfiles[i]->Get(var.c_str());
        HVar->Add(htempVar,weights[i]/weightTot);
    }

    return HVar;

}

void mcruncorr(int filenum = 0)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/franksorted.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }

  // int filenum = 13;
  
  cout<<"opening: "<<listoffiles[filenum].data()<<endl;
  // alicecorrana("/mnt/hadoop/cms/store/user/velicanu/mergedv1_sortedforest/mergesortv1_2.root");
  mccorrana(listoffiles[filenum].data());
  
  // double leadingjetptlow[] =     {0  ,70 ,70 ,120,90 ,60 ,120,90 ,60 };
  // double leadingjetpthigh[] =    {300,90 ,300,300,300,300,300,300,300};
  // double subleadingjetptlow[] =  {70 ,0  ,50 ,50 ,50 ,50 ,50 ,50 ,50 };
  // double subleadingjetpthigh[] = {90 ,300,300,300,300,300,300,300,300};
  // double ptasslow[] = {1,1,1,2,2,2,3,3,3};
  // double ptasshigh[] = {2,2,2,3,3,3,4,4,4};
  // int centmin[] = {0,4,8,16,24};
  // int centmax[] = {4,8,16,24,28};
  
  
  
  double leadingjetptlow[] =     {0  ,70 ,0  ,70 ,0  ,70 ,90 ,90 ,90 };
  double leadingjetpthigh[] =    {300,90 ,300,90 ,300,90 ,300,300,300};
  double subleadingjetptlow[] =  {70 ,0  ,70 ,0  ,70 ,0  ,50 ,50 ,50 };
  double subleadingjetpthigh[] = {90 ,300,300,300,300,300,300,300,300};
  double ptasslow[] =            {1  ,1  ,2  ,2  ,3  ,3  ,1  ,2  ,3  };
  double ptasshigh[] =           {2  ,2  ,3  ,3  ,4  ,4  ,2  ,3  ,4  };
  int centmin[] = {0,4,8,12,16,20};
  int centmax[] = {4,8,12,16,20,24};
  TFile * outf = new TFile(Form("frank_pf3_%d.root",filenum),"recreate");
  
  for(int i = 6 ; i < 9 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    // for(int cent = 0 ; cent < 3 ; ++cent)
    // for(int cent = 0 ; cent < 9 ; ++cent)
    for(int cent = 0 ; cent < 6 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      TH2D * ljtsig = JetTrackSignal(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * ljtbak = JetTrackBackground(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * ljtcorr = (TH2D*)ljtsig->Clone(Form("corr_leadingjet%d_%d_ass%d_%d_cmin%d_cmax%d",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      ljtcorr->Divide(ljtbak);
      ljtcorr->Scale(ljtbak->GetBinContent(ljtbak->FindBin(0,0)));
      ljtcorr->GetXaxis()->SetRange(ljtcorr->GetXaxis()->FindBin(-1.6),ljtcorr->GetXaxis()->FindBin(1.6));
      ljtcorr->GetYaxis()->SetRange(ljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),ljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      
        TH2D * sljtsig = JetTrackSignal(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * sljtbak = JetTrackBackground(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * sljtcorr = (TH2D*)sljtsig->Clone(Form("corr_subleadingjet%d_%d_ass%d_%d_cmin%d_cmax%d",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      sljtcorr->Divide(sljtbak);
      sljtcorr->Scale(sljtbak->GetBinContent(sljtbak->FindBin(0,0)));
      sljtcorr->GetXaxis()->SetRange(sljtcorr->GetXaxis()->FindBin(-1.6),sljtcorr->GetXaxis()->FindBin(1.6));
      sljtcorr->GetYaxis()->SetRange(sljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),sljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      
      
    }
  }
  outf->Write();
  outf->Close();
}

//.........这里部分代码省略.........
  inTree->SetBranchAddress("sigmaIphiIphi",sigmaIphiIphi);
  inTree->SetBranchAddress("swissCrx",swissCrx);
  inTree->SetBranchAddress("seedTime",seedTime);

  inTree->SetBranchAddress("L1_SingleEG2_BptxAND", &L1_SingleEG2_BptxAND);
  inTree->SetBranchAddress("L1_SingleEG5_BptxAND", &L1_SingleEG5_BptxAND);
  inTree->SetBranchAddress("L1_SingleEG8_BptxAND", &L1_SingleEG8_BptxAND);
  inTree->SetBranchAddress("L1_SingleEG12", &L1_SingleEG12);


  Int_t nGen;
  Float_t gen_pt[MAXJETS], gen_eta[MAXJETS];//, gen_phi[MAXJETS];
  Float_t gen_iso[MAXJETS];
  Int_t gen_id[MAXJETS], gen_momId[MAXJETS];

  inTree->SetBranchAddress("nGen",&nGen);
  inTree->SetBranchAddress("gen_pt",gen_pt);
  inTree->SetBranchAddress("gen_eta",gen_eta);
  //inTree->SetBranchAddress("gen_phi",gen_phi);
  inTree->SetBranchAddress("gen_iso",gen_iso);
  inTree->SetBranchAddress("gen_id",gen_id);
  inTree->SetBranchAddress("gen_momId",gen_momId);


  TFile *outFile = new TFile(Form("hist_hydjet_photons_reco.root"),"RECREATE");
  outFile->cd();

  const int nBins = 100;
  const double maxPt = 100;

  TH1D *l1Pt = new TH1D("l1Pt",";L1 p_{T} (GeV)",nBins,0,maxPt);
  TH1D *fPt[3];
  fPt[0] = new TH1D("fPt",";offline p_{T} (GeV)",nBins,0,maxPt);
  fPt[1] = (TH1D*)fPt[0]->Clone("fPt_cen");
  fPt[2] = (TH1D*)fPt[0]->Clone("fPt_periph");
  TH1D *accepted[THRESHOLDS][3];
  TH1D *isoDistribution = new TH1D("isoDistribution",";isolation energy (GeV)",100, -5, 45);
  TH1D *badIsoDistribution = (TH1D*)isoDistribution->Clone("badIsoDistribution");

  TH1D *jetSpectra = new TH1D("jetSpectra","l1 jet (GeV)",64, 0, 256);
  TH1D *badJetSpectra = (TH1D*)jetSpectra->Clone("badJetSpectra");
  TH1D *goodJetSpectra = (TH1D*)jetSpectra->Clone("goodJetSpectra");


  for(int i = 0; i < THRESHOLDS; ++i)
    for(int j = 0; j < 3; ++j)
    {
      accepted[i][j] = new TH1D(Form("accepted_pt%d_%d",(int)L1_THRESHOLD[i],j),";offline p_{T}",nBins,0,maxPt);
    }

  TH2D *corr = new TH2D("corr",";offline p_{T};l1 p_{T}",nBins,0,maxPt,nBins,0,maxPt);
  TH2D *matching = new TH2D("matching",";#Delta #eta;#Delta #phi",100,-5,5,100,0, TMath::Pi() );
  TH2D *matched_bad = (TH2D*)matching->Clone("matched_bad");
  TH2D *absolute_position = new TH2D("absolute_position",";#eta;#phi",50,-5,5,50,-TMath::Pi(), TMath::Pi());

  Long64_t entries = inTree->GetEntries();
  for(Long64_t j = 0; j < entries; ++j)
  {
    if(j % 10000 == 0)
      printf("%lld / %lld\n",j,entries);

    inTree->GetEntry(j);

    double maxl1eta = -999;
    double maxl1phi = -999;
    double maxl1pt = -1;

void pidFitter::processSpecies( string species, int charge, reporter * rp ){
	cout << "[pidFitter." << __FUNCTION__ << "]" << endl;

	string hName = "nSig_" + sName( species, charge );
	// get the pt Binning
	vector<double>pBins = config->getDoubleVector( "binning.p" );

	string useNode = "";

	TH3* h3 = book->get3D( hName );
	int nFits = config->getInt( species + "_Fit:nFits", 1 );
	cout << "Number of Fit Categories: " << nFits << endl;
	for ( int iFit = 1; iFit <= nFits; iFit++ ){

		string optPath = species + "_Fit.opt" + config->getString( species + "_Fit.fit"+ts(iFit)+":options");
		int fBin = config->getInt( species + "_Fit.fit"+ts(iFit)+":min", 1 );
		int lBin = config->getInt( species + "_Fit.fit"+ts(iFit)+":max", pBins.size() );	

		cout << "Fitting P bins ( " << fBin << " --> " << lBin << " ) " << endl;
		for ( int i = fBin; i <= lBin; i ++ ){

			// get the Pt range for title etc.
			double pLow = h3->GetZaxis()->GetBinLowEdge( i );
			double pHi = h3->GetZaxis()->GetBinLowEdge( i + 1 );	

			// look at one Pt bin at a time
			h3->GetZaxis()->SetRange( i, i );
			
			// Get the 2D projection we want
			TH2D* proj;
			proj = (TH2D*)h3->Project3D( "xy" );

			string name = proj->GetName();
			TH2D* cut = (TH2D*) proj->Clone( (name + "cut").c_str() );
			rp->newPage( 1, 2 );
			proj->SetTitle( ( ts( pLow, 4 ) + " #leq " + " P #leq" + ts( pHi, 4 ) ).c_str()  );
			gPad->SetLogz( 1 );
			proj->Draw( "colz" );
			//nProj->Draw("colz");

			// process the square cuts
			makeSquareCuts( cut, optPath + ".squareCut" );

			// draw the distribution after square cuts
			rp->cd( 1, 2 );
			gPad->SetLogz( 1 );
			cut->SetTitle( ( "After 1D Cuts : " + ts( pLow, 4 ) + " #leq " + " P #leq" + ts( pHi, 4 ) ).c_str()  );
			cut->Draw( "colz" );
			rp->savePage( );

			// fit using the dkl algorithm
			if ( config->nodeExists( optPath + ".dkl" ) )
				runDkl( cut, rp, optPath );
			// runs the 2d gaussian fit
			else if ( config->nodeExists( optPath + ".mgf" ) )
				runMultiGauss( cut, rp, species, optPath, i );
			
			
		}

	} 
	
	
}

void BasicCase_tests::MyTest()
{
	float fnum = 2.00001f;
// 	CPPUNIT_FAIL("zxczc");
	CPPUNIT_ASSERT_DOUBLES_EQUAL( fnum, 2.0f, 0.0005 );

	std::string pattern_string("%%d pattern");
	std::string test_string("test pattern");
	std::string replace_string("test");

	TROOT troot("TreeAnalysis","compiled analysisDST macros");
	TApplication app();

// 	h->Draw();
// 	h->Print();
// 	hd->Print();

// 	h->PrintErrors();
// 	PrintErrors(hd);

	TH2DA * h2 = (TH2DA *)h->Clone("new h");
// 	h2->PrintErrors();

	h2->Divide(h);
// 	h2->PrintErrors();
	CPPUNIT_ASSERT_DOUBLES_EQUAL( sqrt(2), h2->GetBinErrorH(1, 1), 0.0005 );
	CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, h2->GetBinContent(1, 1), 0.0005 );

	TH2D * hd2 = (TH2D*)hd->Clone("new hd");
	hd2->Divide(hd);
// 	PrintErrors(hd2);
	CPPUNIT_ASSERT_DOUBLES_EQUAL( sqrt(2), hd2->GetBinError(1, 1), 0.0005 );

	TH2DA * h3 = (TH2DA *)h->Clone("h3: new h");
	TH2DA * h4 = (TH2DA *)h->Clone("h4: new h");

	h3->Reset();
	h4->Reset();

	fill_flat(h3);
	fill_rising(h4);

// 	h3->PrintErrors();
// 	h4->PrintErrors();

	TH2DA * h5 = TH2DA::Efficiency("h5_eff", h3, h4);
// 	h5->PrintErrors();

// 	PR(h5->GetTotalErrorU());
// 	PR(h5->Integral());
// 	PR(h5->GetTotalErrorL());

	CPPUNIT_ASSERT_DOUBLES_EQUAL( (h5->GetTotalErrorH() - h5->GetTotalErrorL()) == 0.0, false, 0.0005 );

// 	TFile * f = new TFile("/tmp/test.root", "RECREATE");
// 	if (f->IsOpen())
// 	{
// 		f->cd();
// 		h5->Write();
// 		h5->GetAsymErrorsGraphX(4)->Write("grx");
// 		h5->GetAsymErrorsGraphY(4)->Write("gry");
// 		f->Close();
// 	}
}

void SimplePlot2D(){

  //string filename = "data/scan_CH1-64_unmasked.root";
  //string filename = "data/scan_CH1-50_masked.root";
  string filename = "test.root";
  string varXname = "VMM #";
  string varYname = "CH #";
  
  // delay count stuff
  int CH = 21;
  double delays[5];
  double count_tot[5];
  double count_right[5];
  for(int i = 0; i < 5; i++){
    delays[i] = double(i)*5.;
    count_tot[i] = 0.;
    count_right[i] = 0.;
  }
  
      

  ///////////////////////////////////////////////////////
  
  TChain* tree = new TChain("MMFE8","MMFE8");

  tree->AddFile(filename.c_str());

  MMFE8Base* base = new MMFE8Base(tree);

  int N = tree->GetEntries();

  TH2D* hist = new TH2D("hist","hist", 8, 0.5, 8.5, 64, 0.5,64.5);
  TH2D* histN = (TH2D*) hist->Clone("norm");
  TH2D* histchch = new TH2D("histchch","histchch", 64, 0.5, 64.5, 64, 0.5,64.5);
  TH2D* histDelay = new TH2D("histN","histN", 31, 9.5, 40.5, 5,-0.5, 4.5);	
  TH2D* histDelayD = new TH2D("histD","histD", 31, 9.5, 40.5, 5,-0.5, 4.5);

  for(int i = 0; i < N; i++){
    base->GetEntry(i);
    
    if(base->CHpulse == CH){
      //count_tot[base->Delay] += 1.;
      count_tot[(base->TPDAC-80)/20] += 1.;
      if(base->CHpulse == base->CHword)
	//count_right[base->Delay] += base->TDO;
      	count_right[(base->TPDAC-80)/20] += base->PDO;
    }

    //histDelayD->Fill(base->CHpulse,base->Delay);
    histDelayD->Fill(base->CHpulse,(base->TPDAC-80)/20);
    if(base->CHpulse == base->CHword)
      //histDelay->Fill(base->CHpulse,base->Delay,base->TDO);
      histDelay->Fill(base->CHpulse,(base->TPDAC-80)/20,base->PDO);

    if((base->CHpulse != base->CHword || true) &&
       base->VMM == 6)
      histchch->Fill(base->CHpulse,base->CHword);

    if(base->CHpulse != base->CHword)
      continue;

    hist->Fill(base->VMM,base->CHpulse,base->PDO);
    histN->Fill(base->VMM,base->CHpulse);
  
  }

  for(int x = 0; x < 8; x++){
    for(int y = 0; y < 64; y++){
      double v = hist->GetBinContent(x+1,y+1);
      double N = histN->GetBinContent(x+1,y+1);
      hist->SetBinContent(x+1,y+1,v/max(int(N),1));
    }
  }
  
  TLatex l;
  //l.NDC();

  TCanvas* can = new TCanvas("can","can",600,500);
  can->SetTopMargin(0.05);
  can->SetLeftMargin(0.12);
  can->SetRightMargin(0.15);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  can->Draw();
  can->SetGridx();
  can->SetGridy();
  
  can->cd();

  hist->Draw("COLZ");
  
  hist->GetXaxis()->SetTitle(varXname.c_str());
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->SetTitle(varYname.c_str());
  hist->GetYaxis()->CenterTitle();
  hist->GetYaxis()->SetTitleOffset(1.4);
  hist->GetYaxis()->CenterTitle();
  //hist->GetYaxis()->SetRangeUser(0.,hist->GetMaximum()*1.1) ;

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

void stdruncorr(int filenum = 0)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/mcsorted.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }

  
  cout<<"opening: "<<listoffiles[filenum].data()<<endl;
  stdcorrana(listoffiles[filenum].data());
  double leadingjetptlow[] =     {100,100,100,100,100,100,100,100,100};
  double leadingjetpthigh[] =    {120,120,120,300,300,300,300,300,300};
  double subleadingjetptlow[] =  {50 ,50 ,50 ,50 ,50 ,50 ,50 ,50 ,50 };
  double subleadingjetpthigh[] = {120,120,120,300,300,300,300,300,300};
  double ptasslow[] =            {2  ,3  ,5  ,2  ,3  ,3  ,1  ,2  ,3  };
  double ptasshigh[] =           {3  ,5  ,8  ,3  ,4  ,4  ,2  ,3  ,4  };


  int centmin[] = {0,4,12,20};
  int centmax[] = {4,12,20,40};
  
  float ajmin[] = { 0.00, 0.13, 0.24, 0.35 };
  float ajmax[] = { 0.13, 0.24, 0.35, 1.00 };
 
  
  TFile * outf = new TFile(Form("stdmcv2_%d.root",filenum),"recreate");
  
  for(int i = 0 ; i < 3 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    for(int cent = 0 ; cent < 4 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      for(int aj = 0 ; aj < 4 ; ++aj)
      {
        cout<<"aj iteration "<<aj<<endl;
        
        TH2D * ljtsig = JetTrackSignal    (0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * ljtbak = JetTrackBackground(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * ljtcorr = (TH2D*)ljtsig->Clone(Form("corr_leadingjet%d_%d_ass%d_%d_cmin%d_cmax%d_ajmin%2.2f_ajmax%2.2f",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent],ajmin[aj],ajmax[aj]));
        ljtcorr->Divide(ljtbak);
        ljtcorr->Scale(ljtbak->GetBinContent(ljtbak->FindBin(0,0)));
        ljtcorr->GetXaxis()->SetRange(ljtcorr->GetXaxis()->FindBin(-1.6),ljtcorr->GetXaxis()->FindBin(1.6));
        ljtcorr->GetYaxis()->SetRange(ljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),ljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
        
        TH2D * sljtsig = JetTrackSignal    (1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * sljtbak = JetTrackBackground(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * sljtcorr = (TH2D*)sljtsig->Clone(Form("corr_subleadingjet%d_%d_ass%d_%d_cmin%d_cmax%d_ajmin%2.2f_ajmax%2.2f",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent],ajmin[aj],ajmax[aj]));
        sljtcorr->Divide(sljtbak);
        sljtcorr->Scale(sljtbak->GetBinContent(sljtbak->FindBin(0,0)));
        sljtcorr->GetXaxis()->SetRange(sljtcorr->GetXaxis()->FindBin(-4.0),sljtcorr->GetXaxis()->FindBin(4.0));
        sljtcorr->GetYaxis()->SetRange(sljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),sljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      }
    }
  }
  outf->Write();
  outf->Close();
}

void runcorr(int condor_iter, int trackqual)
{

  const int nptbins = 2;
  const int ncentbins = 1;
  const int najbins = 1;
  
  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  // ifstream infile("sortedforests.txt");
  // ifstream infile("doeproposalforests.txt");
  // ifstream infile("ppb_unmerged_minbias_forests.txt");
  ifstream infile("HIRun2013-PromptReco-v1-HLT_PAPixelTracks_Multiplicity190_v1-forest-v2.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }
  bool dostdhists = (condor_iter%(nptbins * ncentbins * najbins) == 0);
  int ptbin = condor_iter % nptbins;
  int centbin = (condor_iter / nptbins) % ncentbins;
  int ajbin = (condor_iter / (nptbins * ncentbins)) % najbins;
  int filenum = (condor_iter / (nptbins * ncentbins * najbins));
  cout << "ipt: " << ptbin << " icent: " << centbin << " iaj: " << ajbin << " filenum: " << filenum << " dostdhists: " << dostdhists << " condor_iter "<<condor_iter<< endl;
 

  int nmin = 190 , nmax = 1000;
  //! for first iteration of forest production
  /*
  if(filenum==0) { nmin = 110 ; nmax = 1000; }
  if(filenum==1) { nmin = 90  ; nmax = 110 ; }
  if(filenum>1 ) { nmin = 0   ; nmax = 35  ; }
  if(filenum>9 ) { nmin = 35  ; nmax = 90  ; }
  */
  //! for second iteration of forest production
  // /*
  // if(filenum<24 ) { nmin = 90  ; nmax = 110 ; }
  // if(filenum<22 ) { nmin = 35  ; nmax = 90  ; }
  // if(filenum<12 ) { nmin = 110 ; nmax = 1000; }
  // if(filenum<10 ) { nmin = 0   ; nmax = 35  ; }
  // */
  //! for second iteration of forest production
  /*
  if(filenum<26 ) { nmin = 90  ; nmax = 110 ; }
  if(filenum<23 ) { nmin = 35  ; nmax = 90  ; }
  if(filenum<13)  { nmin = 110 ; nmax = 1000; }
  if(filenum<10)  { nmin = 0   ; nmax = 35  ; }
  */
  corrana(listoffiles[filenum].data(),trackqual);
  
  
  double pttriglow[] =  {1    ,1    ,1 ,3, 1};
  double pttrighigh[] = {3    ,3    ,3 ,4, 3};
  double ptasslow[] =   {0.25 ,2 ,0.25 ,3, 1};
  double ptasshigh[] =  {0.5  ,3    ,3 ,4, 3};
  
  
  int centmin[] = {1,0,4,8,12,16,20,24,28,32};
  int centmax[] = {1,41,8,12,16,20,24,28,32,36};
  TFile * outf = new TFile(Form("corrhists_trkhfminus_trkqaul%d_nmin%d_nmax%d_ptmin%d_ptmax%d_hfmin%d_hfmax%d_%d.root",trackqual,nmin,nmax,(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],filenum),"recreate");
  
  int i = 0;
  int cent = 0;

  cout<<"cent iteration "<<cent<<endl;
  TH2D * ttsig = HFTrackSignal(pttriglow[ptbin],pttrighigh[ptbin],ptasslow[ptbin],ptasshigh[ptbin],centmin[cent],centmax[cent],nmin,nmax);
  TH2D * ttbak = HFTrackBackground(pttriglow[ptbin],pttrighigh[ptbin],ptasslow[ptbin],ptasshigh[ptbin],centmin[cent],centmax[cent],nmin,nmax);

  // TH1I * hntottrig = new TH1I(Form("nttottrig_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],centmin[cent],centmax[cent]),"",1,0.5,1.5);
  // int myntottrig = GetNTotTrig();
  // hntottrig->Fill(1,myntottrig);
  // cout<<"ntottrig: "<<myntottrig<<endl;

  TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],centmin[cent],centmax[cent]));
  ttcorr->Divide(ttbak);
  ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(4,0)));
  ttcorr->Scale(1/0.0594998609); //! bin width
  ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(1),ttcorr->GetXaxis()->FindBin(7));
  ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
  
  outf->Write();
  outf->Close();
}

//.........这里部分代码省略.........
  TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
  if (!fSave) {
    std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
    return -1;
  }
  
  const Double_t pBoundLow = 0.1;
  const Double_t pBoundUp = 5;

  const Int_t nBins1 = TMath::Ceil(180 / downScaleFactor);
  const Int_t nBins2 = TMath::Ceil(100 / downScaleFactor);
  const Int_t nBins3 = TMath::Ceil(60 / downScaleFactor);
  
  const Int_t nPbinsForMap = nBins1 + nBins2 + nBins3;
  Double_t binsPforMap[nPbinsForMap + 1];
  
  Double_t binWidth1 = (1.0 - pBoundLow) / nBins1;
  Double_t binWidth2 = (2.0 - 1.0 ) / nBins2;
  Double_t binWidth3 = (pBoundUp - 2.0) / nBins3;
  
  for (Int_t i = 0; i < nBins1; i++)  {
    binsPforMap[i] = pBoundLow + i * binWidth1;
  }
  for (Int_t i = nBins1, j = 0; i < nBins1 + nBins2; i++, j++)  {
    binsPforMap[i] = 1.0 + j * binWidth2;
  }
  for (Int_t i = nBins1 + nBins2, j = 0; i < nBins1 + nBins2 + nBins3; i++, j++)  {
    binsPforMap[i] = 2.0 + j * binWidth3;
  }
  binsPforMap[nPbinsForMap] = pBoundUp;

  TH2D* hPull = new TH2D("hPull", "Pull vs. p_{TPC} integrated over tan(#Theta);p_{TPC} (GeV/c);Pull", nPbinsForMap, binsPforMap, 
                         plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
  TH2D* hPullAdditionalCorr = (TH2D*)hPull->Clone("hPullAdditionalCorr");
  hPullAdditionalCorr->SetTitle("Pull vs. p_{TPC} integrated over tan(#Theta) with additional dEdx correction w.r.t. tan(#Theta)");
  /*
  const Int_t nThetaHistos = 3;
  TH2D* hPullTheta[nThetaHistos];
  TH2D* hPullAdditionalCorrTheta[nThetaHistos];
  Double_t tThetaLow[nThetaHistos] = { 0.0, 0.4, 0.9 };
  Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.5, 1.0 };
  */
  const Int_t nThetaHistos = 10;
  TH2D* hPullTheta[nThetaHistos];
  TH2D* hPullAdditionalCorrTheta[nThetaHistos];
  Double_t tThetaLow[nThetaHistos] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
  Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
  

  for (Int_t i = 0; i < nThetaHistos; i++)    {
    hPullTheta[i] = new TH2D(Form("hPullTheta_%d", i),
                             Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f;p_{TPC} (GeV/c);Pull", tThetaLow[i], tThetaHigh[i]),
                             nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);

    hPullAdditionalCorrTheta[i] =
      new TH2D(Form("hPullAdditionalCorrTheta_%d", i),
               Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f with additional dEdx correction w.r.t. tan(#Theta);p_{TPC} (GeV/c);Pull",
                    tThetaLow[i], tThetaHigh[i]),
               nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
  }

  
  
  
  
  

void sysError(
    TString inFileName="jfh_HCPR_J50U_Cent30to100_Aj0to100_SubEtaRefl.root",
    Int_t compMode = 0, // Compare mode: 0 reco-genSig, 1 reco-genAll, 3 genAll-genSig, 4 calo_genp-allGen
    Int_t sysMode = 0, // Plot mode: 0 for simple plot, 1 for difference
    TString outdir = ".",
    TString title = "test"
    ) {
  // ===============================================
  // Inputs
  // ===============================================
  TFile *f = new TFile(inFileName);
  TString inFileNameStrip(inFileName); inFileNameStrip.ReplaceAll(".root","");
  TString inFileNameGen(inFileName);
  if (compMode==0) { 
    inFileNameGen.ReplaceAll("djcalo","djcalo_genp");
    inFileNameGen.ReplaceAll("HydjetAll","HydjetSig");
  }
  if (compMode==1) inFileNameGen.ReplaceAll("djcalo","djcalo_genp");
  else if (compMode==3) inFileNameGen.ReplaceAll("HydjetAll","HydjetSig");
  else if (compMode==5) inFileNameGen.ReplaceAll("djcalo_genp","djgen");
  TFile *fgen = new TFile(inFileNameGen);
  TString inFileNameStripGen(inFileNameGen); inFileNameStripGen.ReplaceAll(".root","");
  cout << "==========================================================" << endl;
  cout << "Compare: " << inFileName << endl
       << "         vs " << endl
       << "         " << inFileNameGen << endl;
  cout << "==========================================================" << endl;

  // ===============================================
  // Setup
  // ===============================================
  TString tag=Form("sysError_%s_%s_%d_%d",inFileNameStrip.Data(),title.Data(),compMode,sysMode);

  // ===============================================
  // Analyze
  // ===============================================
  TH2D * hPtPNDR = (TH2D*) f->Get("hPtPNDR");
  TH2D * hPtPADR = (TH2D*) f->Get("hPtPADR");
  TH2D * hPtPNDRBg = (TH2D*) f->Get("hPtPNDRBg");
  TH2D * hPtPADRBg = (TH2D*) f->Get("hPtPADRBg");
  TH2D * hPtPNDRSub = (TH2D*)hPtPNDR->Clone(tag+"hPtPNDRSub");
  TH2D * hPtPADRSub = (TH2D*)hPtPADR->Clone(tag+"hPtPADRSub");
  hPtPNDRSub->Add(hPtPNDR,hPtPNDRBg,1,-1);
  hPtPADRSub->Add(hPtPADR,hPtPADRBg,1,-1);

  TH2D * hPtPNDRGen = (TH2D*) fgen->Get("hPtPNDR");
  TH2D * hPtPADRGen = (TH2D*) fgen->Get("hPtPADR");
  TH2D * hPtPNDRBgGen = (TH2D*) fgen->Get("hPtPNDRBg");
  TH2D * hPtPADRBgGen = (TH2D*) fgen->Get("hPtPADRBg");
  TH2D * hPtPNDRSubGen = (TH2D*)hPtPNDRGen->Clone(tag+"hPtPNDRSub");
  TH2D * hPtPADRSubGen = (TH2D*)hPtPADRGen->Clone(tag+"hPtPADRSub");
  hPtPNDRSubGen->Add(hPtPNDRGen,hPtPNDRBgGen,1,-1);
  hPtPADRSubGen->Add(hPtPADRGen,hPtPADRBgGen,1,-1);

  // ===============================================
  // Draw
  // ===============================================
  // Get Pt info
  Int_t numPtBins=hPtPNDR->GetNbinsX();
  TH1D * hPt = (TH1D*)hPtPNDR->ProjectionX("hPt");
  /*
  cout << "Pt bins: " << numPtBins << endl;
  for (Int_t i=0; i<numPtBins+2; ++i) {
    cout << "Pt Bin " << i << " Low Edge: " << hPt->GetBinLowEdge(i) << endl;
  }
  */

  // What pt bins to draw
  const Int_t numPtBinsDraw=3;

  TCanvas * c6 = new TCanvas("c"+tag,"c"+tag,1400,500);
  c6->Divide(3,1);
  for (Int_t i=0; i<numPtBinsDraw; ++i) {
    Int_t iBeg,iEnd;
    if (i==0) { iBeg=2; iEnd=3;}
    if (i==1) { iBeg=4; iEnd=4;}
    if (i==2) { iBeg=5; iEnd=numPtBins;}
    cout << "Bin: " << iBeg <<  " to " << iEnd << endl;
    TH1D * hNr = (TH1D*)hPtPNDRSub->ProjectionY(tag+Form("hPNDRSub_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hAw = (TH1D*)hPtPADRSub->ProjectionY(tag+Form("hPADRSub_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hNrGen = (TH1D*)hPtPNDRSubGen->ProjectionY(tag+Form("hPNDRSubGen_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hAwGen = (TH1D*)hPtPADRSubGen->ProjectionY(tag+Form("hPADRSubGen_%d_%d",iBeg,iEnd),iBeg,iEnd);
    if (sysMode==1) {
      hNr->Add(hNrGen,-1);
      hAw->Add(hAwGen,-1);
    }
    if (sysMode==2) {
      hNr->Divide(hNrGen);
      hAw->Divide(hAwGen);
    }
    // Print
    cout << Form("%.1f < P_{T} < %.1f GeV: ",hPt->GetBinLowEdge(iBeg),hPt->GetBinLowEdge(iEnd+1))
      << " SigSubBkg Integral - Nr: " << hNr->Integral() << " Aw: " << hAw->Integral() << endl
      << " Gen - Nr: " << hNrGen->Integral() << " Aw: " << hAwGen->Integral() << endl;
    // Styles
    hNr->SetMarkerColor(kRed);
    hNr->SetLineColor(kRed);
    hAw->SetMarkerColor(kBlue);
    hAw->SetLineColor(kBlue);
    hAwGen->SetLineStyle(2);
//.........这里部分代码省略.........

void signalEff2012_PDF_details(const TString what="CTEQMSTW",const TString sample="T1bbbb",const int minnjets=3,const bool useisr=false) {

  TString njetsstring=".";

  if (minnjets==3) {
    //do nothing
  }
  else if (minnjets==5) {
    njetsstring = ".minnjets5.";
  }
  else {
    cout<<" minnjets = "<<minnjets<<" is not ok"<<endl;
    return;
  }

  TString stub1="eventcounts2x2.mergebbins";
  if (sample.Contains("pMSSM") ||sample.Contains("T1ttcc") ||sample.Contains("14TeV")||sample.Contains("TChi")) stub1="eventcounts.mergebbins";

  TString stub2=stub1;
  stub2+=".withpdfs";

  if ( useisr) {
    stub1+=".Isr0";
    stub2+=".Isr0";
  }

  //  TString nominalstub="CSVM_PF2PATjets_JES0_JER0_PFMETTypeI_METunc0_PUunc0_BTagEff05_HLTEff0."; //old
  TString nominalstub="JES0_JERbias_PFMETTypeI_METunc0_PUunc0_hpt20."; //new, and with JERbias
  TString f0file = stub1+njetsstring+nominalstub+sample+".root";
  TString fpdffile = stub2+njetsstring+nominalstub+sample+".root";

  TFile f0(f0file);
  TFile fpdf(fpdffile);

  vector<TH2D*> effratios;

   for (int ih = 0; ih<fpdf.GetListOfKeys()->GetEntries(); ih++) {
     TString histname = fpdf.GetListOfKeys()->At(ih)->GetName();
     if (! (histname.BeginsWith("events_")&&histname.Contains(what))) continue;

     TString histnametotal = histname;
     histnametotal.ReplaceAll("events_","eventstotal_");

     TString histname0=histname;
     histname0.ReplaceAll("_"+what+"0","");

     TString histname0total = histnametotal;
     histname0total.ReplaceAll("_"+what+"0","");


     TH2D* h0c = (TH2D*) f0.Get(histname0);
     TH2D* h0t = (TH2D*) f0.Get(histname0total);
     
     TH2D* hpdfc = (TH2D*) fpdf.Get(histname);
     TH2D* hpdft = (TH2D*) fpdf.Get(histnametotal);
     
     //now calculate efficiency
     TH2D* h0r = (TH2D*) h0c->Clone("h0r_"+histname);
     h0r->Reset();
     h0r->Divide(h0c,h0t);

     TH2D* hpdfr = (TH2D*) hpdfc->Clone("hpdfr_"+histname);
     hpdfr->Reset();
     hpdfr->Divide(hpdfc,hpdft);

     //and then the ratio of efficiencies
     TString rname = histname;
     rname.ReplaceAll("events_","effratio_");
     TH2D*  heffratio = (TH2D*) h0r->Clone(rname);
     heffratio->Reset();
     heffratio->Divide(hpdfr,h0r);
     effratios.push_back(heffratio);

   }

   TFile fout("signalsyst_PDF_"+what+"."+sample+njetsstring+"root","recreate");
   for (unsigned int ii=0; ii<effratios.size();ii++) {
     effratios.at(ii)->Write();
   }
   fout.Close();


   /* quick and dirty version
  TH2D* h0c = (TH2D*) f0.Get("events_b1_HT1000to100000_MET350to100000_2x2");
  TH2D* h0t = (TH2D*) f0.Get("eventstotal_b1_HT1000to100000_MET350to100000_2x2");

  TH2D* h0r = (TH2D*) h0c->Clone("h0r");
  h0r->Reset();
  h0r->Divide(h0c,h0t);

  TH2D* hpdfc = (TH2D*) fpdf.Get("events_b1_HT1000to100000_MET350to100000_CTEQMSTW0_2x2");
  TH2D* hpdft = (TH2D*) fpdf.Get("eventstotal_b1_HT1000to100000_MET350to100000_CTEQMSTW0_2x2");

  TH2D* hpdfr = (TH2D*) hpdfc->Clone("hpdfr");
  hpdfr->Reset();
  hpdfr->Divide(hpdfc,hpdft);


  TH2D*  heffratio = (TH2D*) h0r->Clone("heffratio");
  heffratio->Reset();
//.........这里部分代码省略.........