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

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

//========================
void PlotRatio2D(Int_t iEff, Int_t iRap){

  gStyle->SetOptStat(0);
  gStyle->SetPaintTextFormat("5.2f");

  Char_t name[100];
  sprintf(name, "ratio2D_%sEff_%d", effName[iEff], iRap);
  TCanvas *c1 = new TCanvas(name, name);

  sprintf(name, "hRatio2D_%sEff_rap%d", effName[iEff], iRap);
  TH2D *hRatio = (TH2D *) hEffTP2D[iEff][iRap]->Clone(name);
  hRatio->Divide(hEffMCTruth2D[iEff][iRap]);
  hRatio->SetMinimum(0.);
  hRatio->SetMaximum(1.5);
  hRatio->Draw("colz text");

  sprintf(name, "Figures/ratio2D_%sEff_rap%d.pdf", effName[iEff], iRap);
  c1->Print(name);
}

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();
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
			}
			if (ynreco) {
				hRecoPtY->Fill(yreco,ptreco,weightF);
				hRecoPt->Fill(ptreco, weightF);
				hRecoRap->Fill(yreco, weightF); 
				hRecoLxy->Fill(lxy, weightF); 
			}
		}else { 
			if (doWeight) {weightF = getWeight(isPrompt, isPair, pt,eta);}
//				cout << "weightF = " <<weightF << endl;
			if (yngen)	{
				hGenPtY->Fill(eta,pt,weightF);
				hGenPt->Fill(pt,weightF);
				hGenRap->Fill(eta,weightF); 
				hGenLxy->Fill(lxy,weightF); 
			}
			if (ynreco) {
				hRecoPtY->Fill(etareco,ptreco,weightF);
				hRecoPt->Fill(ptreco,weightF);
				hRecoRap->Fill(etareco,weightF);
				hRecoLxy->Fill(lxy, weightF); 
			}
		}
	}// end of event loop

	cout << "hGenPt integral = " << hGenPt->Integral() << endl;
	cout << "hRecoPt integral = " << hRecoPt->Integral() << endl;
	cout << "hGenRap integral = " << hGenRap->Integral() << endl;
	cout << "hRecoRap integral = " << hRecoRap->Integral() << endl;
	cout << "hGenLxy integral = " << hGenLxy->Integral() << endl;
	cout << "hRecoLxy integral = " << hRecoLxy->Integral() << endl;

	//calculate efficiency
	hEffPtY->Divide(hRecoPtY,hGenPtY,1,1,"b");
	hEffPt->Divide(hRecoPt,hGenPt,1,1,"b");
	hEffRap->Divide(hRecoRap,hGenRap,1,1,"b");
	hEffLxy->Divide(hRecoLxy,hGenLxy,1,1,"b");

	//Legend
	TLegend *legUR = new TLegend(0.56,0.68,0.90,0.90,NULL,"brNDC");
	TLegend *legUL = new TLegend(0.17,0.68,0.51,0.90,NULL,"brNDC");
	TLegend *legBM = new TLegend(0.40,0.20,0.75,0.35,NULL,"brNDC");
	TLegend *legBR = new TLegend(0.56,0.20,0.90,0.42,NULL,"brNDC");
	SetLegendStyle(legUR);
	SetLegendStyle(legUL);
	SetLegendStyle(legBM);
	SetLegendStyle(legBR);

	TCanvas *c1 = new TCanvas("c1","c1",600,600); 
	
	//draw 2D Pt-y
	TPaletteAxis* pal; 
	hEffPtY->SetMinimum(0.0);
	hEffPtY->SetMaximum(1.0);
	if (isPair) hEffPtY->GetXaxis()->SetTitle("y");
	else hEffPtY->GetXaxis()->SetTitle("#eta");
	hEffPtY->GetXaxis()->CenterTitle();
	hEffPtY->GetYaxis()->SetTitle("p_{T} (GeV)");
	hEffPtY->Draw("colz");
	c1->Update();
	pal = (TPaletteAxis*)hEffPtY->GetListOfFunctions()->FindObject("palette"); 
	pal->SetX2NDC(0.92);
	c1->Modified();
	c1->Update();
	c1->SaveAs(Form("%s/EffPtY_isPtCut%d_doWeight%d_isEmbedded%d_%s_%s_%s.pdf",strDir.c_str(),(int)isPtCut,(int)doWeight,(int)isEmbedded,strPrompt.c_str(),strEmbd.c_str(),strPair.c_str()));
	c1->Clear();

int genbod_diff(const char* filename1, const char* filename2, bool ifSave = false)
{
  TFile* f1 = new TFile(filename1,"READ");
  TFile* f2 = new TFile(filename2,"READ");

  TH2D* numF1 = ((TH2D*) f1->Get("cnumepNonIdEP"));
  TH2D* numF2 = ((TH2D*) f2->Get("cnumepNonIdEP"));
  TH2D* denF1 = ((TH2D*) f1->Get("cdenepNonIdEP"));
  TH2D* denF2 = ((TH2D*) f2->Get("cdenepNonIdEP"));
  TH2D* numB1 = ((TH2D*) f1->Get("cnumepNonIdEPTrue"));
  TH2D* numB2 = ((TH2D*) f2->Get("cnumepNonIdEPTrue"));
  TH2D* denB1 = ((TH2D*) numF1->Clone());
  TH2D* denB2 = ((TH2D*) numF2->Clone());
  TH2D* numC1 = ((TH2D*) numB1->Clone());
  TH2D* numC2 = ((TH2D*) numB2->Clone());
  TH2D* denC1 = ((TH2D*) f1->Get("cdenepNonIdEPTrue"));
  TH2D* denC2 = ((TH2D*) f2->Get("cdenepNonIdEPTrue"));

  double eta = 1.1;
  gStyle->SetOptStat(111);
  //TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,800,600);
   TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,1600,1200);
  canv->Divide(3,3);
  //TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,800,600);
  TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,1600,1200);
  canvND->Divide(2,2);

  //////   Pure Corr fun
  TH2D* numdiff = (TH2D*)numF1->Clone();
  numdiff->Divide(numF2);
  canvND->cd(1); numdiff->DrawCopy("colz");
  TH2D* dendiff = (TH2D*)denF1->Clone();
  dendiff->Divide(denF2);
  canvND->cd(2); dendiff->DrawCopy("colz");
  double sF1 = numF1->GetEntries()/denF1->GetEntries();
  numF1->Divide(denF1);
  numF1->Scale(1./sF1);
  numF1->GetYaxis()->SetRangeUser(-eta,eta);
  double sF2 = numF2->GetEntries()/denF2->GetEntries();
  numF2->Divide(denF2);
  numF2->Scale(1./sF2);
  numF2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(1);
  TH2D* corrF1 = (TH2D*) numF1->Clone();
  corrF1->Draw("surf1");
  DrawNicely(corrF1, 0,0, filename1);
  canv->cd(2);
  TH2D* corrF2 = (TH2D*) numF2->Clone();
  corrF2->Draw("surf1");
  DrawNicely(corrF2, 0,0, filename2);
  canv->cd(3);
  numF1->Divide(numF2);
  DrawNicely(numF1, 0,0, "Femtoscopic component");
  numF1->Draw("colz");
  
  ////// Pure Bkg
  TH2D* numdiffTrue = (TH2D*)numB1->Clone();
  numdiffTrue->Divide(numB2);
  canvND->cd(3); numdiffTrue->DrawCopy("colz");
  double sB1 = numB1->GetEntries()/denB1->GetEntries();
  numB1->Divide(denB1);
  numB1->Scale(1./sB1);
  numB1->GetYaxis()->SetRangeUser(-eta,eta);
  double sB2 = numB2->GetEntries()/denB2->GetEntries();
  numB2->Divide(denB2);
  numB2->Scale(1./sB2);
  numB2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(4);
  TH2D* corrB1 = (TH2D*) numB1->Clone();
  corrB1->Draw("surf1");
  DrawNicely(corrB1, 0,0, filename1);
  canv->cd(5);
  TH2D* corrB2 = (TH2D*) numB2->Clone();
  corrB2->Draw("surf1");
  DrawNicely(corrB2, 0,0, filename2);
  canv->cd(6);
  numB1->Divide(numB2);
  DrawNicely(numB1, 0,0, "Background component");
  numB1->DrawCopy("colz");
  
  //////   Norm Corr fun
  TH2D* dendiffTrue = (TH2D*)denC1->Clone();
  dendiffTrue->Divide(denC2);
  canvND->cd(4); dendiffTrue->DrawCopy("colz");
  double sC1 = numC1->GetEntries()/denC1->GetEntries();
  numC1->Divide(denC1);
  numC1->Scale(1./sC1);
  numC1->GetYaxis()->SetRangeUser(-eta,eta);
  double sC2 = numC2->GetEntries()/denC2->GetEntries();
  numC2->Divide(denC2);
  numC2->Scale(1./sC2);
  numC2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(7);
  TH2D* corrC1 = (TH2D*) numC1->Clone();
  corrC1->Draw("surf1");
  DrawNicely(corrC1, 0,0, filename1);
  canv->cd(8);
  TH2D* corrC2 = (TH2D*) numC2->Clone();
  corrC2->Draw("surf1");
  DrawNicely(corrC2, 0,0, filename2);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if(nTauGood!=0)  continue;
      CutFlow->Fill(1.);
      // if(fabs(LepGood_pdgId[0]* LepGood_pdgId[1])<122)   continue;
      CutFlow->Fill(2.);
      if(lep.Pt()<5||fabs(lep.Eta())>1.5||lep.Pt()>60) continue;
      CutFlow->Fill(3.);
      if(fabs(LepGood_pdgId[0])==11&&lep.Pt()<7) continue;
      CutFlow->Fill(4.);
      // check that a cut was ahead in nLepton !!!!!!!!
      TLorentzVector secondLep;
      secondLep.SetPtEtaPhiM(LepGood_pt[1], LepGood_eta[1], LepGood_phi[1],0.1);
      if(secondLep.Pt()<3||fabs(secondLep.Eta())>1.5||secondLep.Pt()>60) continue;
      CutFlow->Fill(5.);
      // lepton cuts
      if(LepGood_relIso[0]*lep.Pt()>10.) continue;
      if(LepGood_relIso[1]*secondLep.Pt()>5.) continue;
      if(LepGood_relIso[1]>.5) continue;
      CutFlow->Fill(6.);
      // JET CUTS
      if(nJet>0){
	if (Jet_pt[0] < 150) continue;
	if(fabs(Jet_eta[0]) > 2.4) continue;
      }
      if(nJet>2){
	if (Jet_pt[2] > 60) continue;
      }
      CutFlow->Fill(7.);
      if(met.Pt()/HT30<2./3.) continue;
      if(nb40!=0) continue;    
      CutFlow->Fill(9.);
    // preselection syncronized to  ++++++++++++++++++++

      if(MT>60&&MT<100) continue;
 
      if (weight==1) puWeight=1; // to not effect data
      float pairmass = DiTau_InvMass(met,lep,secondLep,0);     
      if(lep.Pt()>25||secondLep.Pt()>15) continue;

      if(fabs(LepGood_dz[1]) > 0.02||fabs(LepGood_dz[0]) > 0.02 ||fabs(LepGood_dxy[1])> 0.02  || fabs(LepGood_dxy[0])> 0.02 ) continue;


      // define some over and underflow 
      if (pairmass>2000) pairmass=1999.;
      if (pairmass<0) pairmass=0.;

      LepTwoZmass->Fill(pairmass,weight*puWeight);
      //      cout << LepGood_pdgId[0]* LepGood_pdgId[1] << " Ids "<<endl;

      if(LepGood_pdgId[0]*LepGood_pdgId[1]>121) //same sign
	  {
	 
	      if((pairmass>160||pairmass<20)){
	        
	    
	      CutFlow->Fill(10., weight*puWeight);
	      scanB->Fill(GenSusyMStop,GenSusyMNeutralino);
	      // if (weight<0) weight=1;
	      Sig->Fill(LepGood_sip3d[1],weight*puWeight);
	      Sig0->Fill(LepGood_sip3d[0],weight*puWeight);
	      IP3D->Fill(LepGood_ip3d[0],weight*puWeight);
	      Pt->Fill(secondLep.Pt(),weight*puWeight);
	      Iso->Fill(LepGood_relIso[1]*secondLep.Pt(),weight*puWeight);
	      Dxy->Fill(fabs(LepGood_dxy[1]),weight*puWeight);
	      Dz->Fill(fabs(LepGood_dz[1]),weight*puWeight);
	      DxyDz->Fill(fabs(LepGood_dxy[1]),fabs(LepGood_dz[1]),weight*puWeight);

	      if(nb30L==0)
		{	
		  scanC->Fill(GenSusyMStop,GenSusyMNeutralino);
		}
	      for (int i=0;i<nJet;i++)
		{
		  if(Jet_btagCSV[i]>0.244 ) {
		   LooseBPt->Fill(Jet_pt[i],weight*puWeight );		    
		  }
		  //	   QG->Fill(Jet_quarkGluonID[i]);
		}

	  }


	else // same sign leptons
	  {	  
	   
	  }
      }
          
   }


   scanA->Divide(scan);
   scanB->Divide(scan);
   scanC->Divide(scan);

   //  CutFlow->DrawCopy();

   myFile->Write();
   return (myFile);
      // if (Cut(ientry) < 0) continue;
}

//.........这里部分代码省略.........
    reco->SetBranchAddress("rmin",&rmin);
    reco->SetBranchAddress("jtpt",&maxJetPt);
    reco->SetBranchAddress("trkStatus",&trkStatus);
    reco->SetBranchAddress("nEv",&nEv); 
    for(int i = 0; i<reco->GetEntries(); i++)
    {
      reco->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      if(trkStatus<-100) continue;
      mrecoPre[0]->Fill(pt,weight);
      mrecoPre2[1]->Fill(eta,phi,weight); 
      mrecoPre[2]->Fill(centPU,weight);
      mrecoPre[3]->Fill(maxJetPt,weight);
      mrecoPre[4]->Fill(eta,weight); 
      mrecoPre[5]->Fill(rmin,weight);
      mrecoPre2[7]->Fill(eta,pt,weight);
    }
  
    //Secondary calculation (no iterations)
    std::cout << "Quickly calculating the Secondary Rate from the reco tree (No further iterations needed)" << std::endl;
    TH2D * Secondary_Matched = new TH2D("Secondary_Matched",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax,24,-2.4,2.4); 
    TH2D * Secondary_Secondaries = new TH2D("Secondary_Secondaries",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax,24,-2.4,2.4); 
    for(int i = 0; i<reco->GetEntries(); i++)
    {
      reco->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      if(trkStatus>-100)
      {
        Secondary_Matched->Fill(pt,eta,weight);
        if(trkStatus==-99) Secondary_Secondaries->Fill(pt,eta,weight);
      }
    }
    TH2D * Secondary = (TH2D*)Secondary_Secondaries->Clone("SecondaryRate");
    Secondary->Divide(Secondary_Matched);
    Secondary->SetDirectory(histFile);
    Secondary_Matched->SetDirectory(histFile);
    Secondary_Secondaries->SetDirectory(histFile);
    //end Secondary Reco calculation

    //multiReco calculation (no iterations)
    std::cout << "Quickly calculating the Multiple Reconstruction Rate from the gen tree (No further iterations needed)" << std::endl;
    TH1D * MultiGen = new TH1D("MultiGen",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax); 
    TH1D * MultiReco = new TH1D("MultiMatchedReco",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax); 
    for(int i = 0; i<gen->GetEntries(); i++)
    {
      gen->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      if(pNRec>-1)
      {
        MultiGen->Fill(pt,weight);
        if(pNRec>1) MultiReco->Fill(pt,(pNRec-1)*weight);
      }
    }
    TH1D * Multi = (TH1D*)MultiReco->Clone("MultipleRecoRate");
    Multi->Divide(MultiGen);
    Multi->SetDirectory(histFile);
    MultiReco->SetDirectory(histFile);
    MultiGen->SetDirectory(histFile);
    //end Multiple Reco calculation

    skim->Close();
    histFile->Write(); 
  }
	  
  //redundant for first step, but needed if the gen file was made and saved previously 
  std::cout << "Loading appropriate information for denominator (gen for efficiency, reco for fake)..." << std::endl;

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

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

//.........这里部分代码省略.........
  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) ;

  TCanvas* canN = new TCanvas("canN","canN",600,500);
  canN->SetTopMargin(0.05);
  canN->SetLeftMargin(0.12);
  canN->SetRightMargin(0.15);

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

  histN->Draw("COLZ");
  
  histN->GetXaxis()->SetTitle(varXname.c_str());
  histN->GetXaxis()->CenterTitle();
  histN->GetYaxis()->SetTitle(varYname.c_str());
  histN->GetYaxis()->CenterTitle();
  histN->GetYaxis()->SetTitleOffset(1.4);
  histN->GetYaxis()->CenterTitle();
			     
  TCanvas* canchch = new TCanvas("canchch","canchch",600,500);
  canchch->SetTopMargin(0.05);
  canchch->SetLeftMargin(0.12);
  canchch->SetRightMargin(0.15);

  canchch->Draw();
  canchch->SetGridx();
  canchch->SetGridy();
  
  canchch->cd();
  histchch->Draw("COLZ");
  
  histchch->GetXaxis()->SetTitle("CH pulsed");
  histchch->GetXaxis()->CenterTitle();
  histchch->GetYaxis()->SetTitle("CH data");
  histchch->GetYaxis()->CenterTitle();
  histchch->GetYaxis()->SetTitleOffset(1.4);
  histchch->GetYaxis()->CenterTitle();
  histchch->GetZaxis()->SetTitle("Number of data events");
  histchch->GetZaxis()->SetTitleOffset(1.4);
  histchch->GetZaxis()->CenterTitle();
 
  l.DrawLatex(.54,65.2,"VMM 2");

  TCanvas* can_delay = new TCanvas("can_delay","can_delay",600,500);
  can_delay->Draw();
  can_delay->cd();

  for(int i = 0; i < 5; i++)
    count_tot[i] = count_right[i]/count_tot[i];
  TGraph* gr = new TGraph(5,delays,count_tot);
  gr->SetMarkerSize(4);
  gr->SetMarkerStyle(5);
  gr->Draw("AP");

  histDelay->Divide(histDelayD);

  TCanvas* canDelay = new TCanvas("canDelay","canDelay",600,500);
  canDelay->SetTopMargin(0.05);
  canDelay->SetLeftMargin(0.12);
  canDelay->SetRightMargin(0.15);

  canDelay->Draw();
  canDelay->SetGridx();
  canDelay->SetGridy();
  
  canDelay->cd();
  histDelay->Draw("COLZ");
  
  histDelay->GetXaxis()->SetTitle("CH pulsed");
  histDelay->GetXaxis()->CenterTitle();
  histDelay->GetYaxis()->SetTitle("Delay Count");
  histDelay->GetYaxis()->CenterTitle();
  histDelay->GetYaxis()->SetTitleOffset(1.4);
  histDelay->GetYaxis()->CenterTitle();
  histDelay->GetZaxis()->SetTitle("Fraction zeroes");
  histDelay->GetZaxis()->SetTitleOffset(1.4);
  histDelay->GetZaxis()->CenterTitle();
}

int genbod_drawNice2(const char* filename1, const char* filename2, bool ifSave = false)
{
  TFile* f1 = new TFile(filename1,"READ");
  TFile* f2 = new TFile(filename2,"READ");

  TH2D* numF1 = ((TH2D*) f1->Get("cnumepNonIdEP"));
  TH2D* numF2 = ((TH2D*) f2->Get("cnumepNonIdEP"));
  TH2D* denF1 = ((TH2D*) f1->Get("cdenepNonIdEP"));
  TH2D* denF2 = ((TH2D*) f2->Get("cdenepNonIdEP"));
  TH2D* numB1 = ((TH2D*) f1->Get("cnumepNonIdEPTrue"));
  TH2D* numB2 = ((TH2D*) f2->Get("cnumepNonIdEPTrue"));
  TH2D* denB1 = ((TH2D*) numF1->Clone());
  TH2D* denB2 = ((TH2D*) numF2->Clone());
  TH2D* numC1 = ((TH2D*) numB1->Clone());
  TH2D* numC2 = ((TH2D*) numB2->Clone());
  TH2D* denC1 = ((TH2D*) f1->Get("cdenepNonIdEPTrue"));
  TH2D* denC2 = ((TH2D*) f2->Get("cdenepNonIdEPTrue"));

  double eta = 1.1;
  gStyle->SetOptStat(000);//111);

  int drawDifference = 0;

  const char* first = "";//Cons. laws only";
  const char* second = "";//Minijets";

  int drawWhich = 0; // 0 - cons. laws, 1 - femto, 2 - full corr fun

  //TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,800,600);
  TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,(2+drawDifference)*550,400);
  canv->Divide(2+drawDifference,1);

  //TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,800,600);
  // TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,1600,1200);
  // canvND->Divide(2,2);

  if(drawWhich==1){
    //////   Pure Corr fun
    TH2D* numdiff = (TH2D*)numF1->Clone();
    numdiff->Divide(numF2);
    //canvND->cd(1); numdiff->DrawCopy("colz");
    TH2D* dendiff = (TH2D*)denF1->Clone();
    dendiff->Divide(denF2);
    //canvND->cd(2); dendiff->DrawCopy("colz");
    double sF1 = numF1->GetEntries()/denF1->GetEntries();
    numF1->Divide(denF1);
    numF1->Scale(1./sF1);
    numF1->GetYaxis()->SetRangeUser(-eta,eta);
    double sF2 = numF2->GetEntries()/denF2->GetEntries();
    numF2->Divide(denF2);
    numF2->Scale(1./sF2);
    numF2->GetYaxis()->SetRangeUser(-eta,eta);
    canv->cd(1);
    TH2D* corrF1 = (TH2D*) numF1->Clone();
    corrF1->Draw("surf1");
    DrawNicely(corrF1, 0,0, filename1);
    canv->cd(2);
    TH2D* corrF2 = (TH2D*) numF2->Clone();
    corrF2->Draw("surf1");
    DrawNicely(corrF2, 0,0, filename2);
    if(drawDifference){
      canv->cd(3);
      numF1->Divide(numF2);
      DrawNicely(numF1, 0,0, "Femtoscopic component");
      numF1->Draw("colz");
    }
  }
  if(drawWhich==0){
    ////// Pure Bkg
    TH2D* numdiffTrue = (TH2D*)numB1->Clone();
    numdiffTrue->Divide(numB2);
    //canvND->cd(3); numdiffTrue->DrawCopy("colz");
    double sB1 = numB1->GetEntries()/denB1->GetEntries();
    numB1->Divide(denB1);
    numB1->Scale(1./sB1);
    numB1->GetYaxis()->SetRangeUser(-eta,eta);
    double sB2 = numB2->GetEntries()/denB2->GetEntries();
    numB2->Divide(denB2);
    numB2->Scale(1./sB2);
    numB2->GetYaxis()->SetRangeUser(-eta,eta);
    canv->cd(1);
    TH2D* corrB1 = (TH2D*) numB1->Clone();
    corrB1->Draw("surf1");
    DrawNicely(corrB1, 0.7,1.4, first);//filename1);
    canv->cd(2);
    TH2D* corrB2 = (TH2D*) numB2->Clone();
    corrB2->Draw("surf1");
    DrawNicely(corrB2, 0.7,1.4, second);//filename2);
    if(drawDifference){
      canv->cd(3);
      numB1->Divide(numB2);
      DrawNicely(numB1, 0,0, "Background component");
      numB1->DrawCopy("colz");
    }
  }
  if(drawWhich==2){
    //////   Norm Corr fun
    TH2D* dendiffTrue = (TH2D*)denC1->Clone();
    dendiffTrue->Divide(denC2);
    //canvND->cd(4); dendiffTrue->DrawCopy("colz");
//.........这里部分代码省略.........