#define gammaJetAnalysis_cxx
#include "gammaJetAnalysis.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "commonFunctions.h"
#include <TRandom.h>

void gammaJetAnalysis::Loop(Long64_t eventMax)  {
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  for (Long64_t e=0; e<nentries && e<eventMax; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    cout << "Event= " << e << " Weight= " << Weight << endl;
  }
  return;
}


void gammaJetAnalysis::ToggleBranchesOnOff(string command) {
  /* To speed up, we may wish to deactivate most branches, 
     and keep reading only those that we need. */
  assert(command=="on" || command=="off" || command=="onlyNeeded");
  if (fChain == 0) return;
  
  if (command=="on") 
    fChain->SetBranchStatus("*",1);
  else if (command=="off")
    fChain->SetBranchStatus("*",0);
  else if (command=="onlyNeeded") {
    fChain->SetBranchStatus("*",0);
    fChain->SetBranchStatus("Weight",1);
    //count variables
    fChain->SetBranchStatus("NPhot_true",1);
    fChain->SetBranchStatus("NPhot",1);
    fChain->SetBranchStatus("NJets_true",1);
    fChain->SetBranchStatus("NJets",1);
    //eta phi pT of photon
    fChain->SetBranchStatus("PhotEta_true",1);
    fChain->SetBranchStatus("PhotPhi_true",1);
    fChain->SetBranchStatus("PhotEta",1);
    fChain->SetBranchStatus("PhotPhi",1);
    fChain->SetBranchStatus("PhotPt",1);
    fChain->SetBranchStatus("PhotPt_true",1);
    //eta phi pT of jet
    fChain->SetBranchStatus("JetEta_true",1);
    fChain->SetBranchStatus("JetPhi_true",1);
    fChain->SetBranchStatus("JetEta",1);
    fChain->SetBranchStatus("JetPhi",1);
    fChain->SetBranchStatus("JetPt",1);
    fChain->SetBranchStatus("JetEMPt",1);
    fChain->SetBranchStatus("JetPt_true",1);
    //quality variables
    fChain->SetBranchStatus("PhotGood",1);
    fChain->SetBranchStatus("PhotEtcone",1);
    fChain->SetBranchStatus("PhotNTrack",1);
  }
  return;
}


void gammaJetAnalysis::SetWeight(double weight, string outputFilename, bool alsoFilterEvents, vector<Long64_t> events)  {
  //produce a new TTree, based on the fChain here, and fill the Weight brunch with something (it's initially filled with 0).  The new tree is then saved in the filename specified.
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  if(events.size() == 0) { //reweight all events.
    TFile fout(outputFilename.c_str(),"recreate");
    TTree* newTree;
    newTree = fChain->CloneTree(0);
    for (Long64_t e=0; e<nentries; e += prescale) { 
      Long64_t ientry = LoadTree(e); 
      if (ientry < 0) break;
      fChain->GetEntry(e);
      Weight = weight;
      if (alsoFilterEvents) 
	if (!PassesEventFilter()) continue;
      newTree->Fill();
    }
    fout.cd();
    newTree->Write();
  }
  else { //reweight just specific events
    TFile fout(outputFilename.c_str(),"recreate");
    TTree* newTree;
    newTree = fChain->CloneTree(0);
    for (unsigned int i=0; i<events.size(); i++) {
      Long64_t e=events[i];
      Long64_t ientry = LoadTree(e); 
      if (ientry < 0) break;
      fChain->GetEntry(e);
      Weight = weight;
      if (alsoFilterEvents) 
	if (!PassesEventFilter()) continue;
      newTree->Fill();
    }
    fout.cd();
    newTree->Write();
  }
        
  return;
}


bool gammaJetAnalysis::PassesEventFilter() {
  /* those are not cuts for the analysis.  it's a prefilter, used to skim ntuples smaller, and reduce the time it will later take to loop through them.  Basically, it's a very crude cut, removing events that for any consideration will not be selected later by the PassesEventCuts() */
  if (NPhot==0) return false;

  //require at least one tight photon
  unsigned int goodPhotons = 0;
  for (unsigned int iPhoton=0; iPhoton<NPhot; iPhoton++)
    if ((*PhotGood)[iPhoton]==1) goodPhotons++;
  if (goodPhotons==0) return false;

  //the following is to keep subsets of the GJ signal (runs 108001-108004), in such a way that then we will stitch them together continuously, without double-counting some pT regions.
  double hardestPt = (*JetPt_true)[0]/GeV;//can be jet or photon
  if (RunNumber==108001) {
    if (hardestPt<20 || hardestPt>70) return false;
  }
  else if (RunNumber==108002) {
    if (hardestPt<70 || hardestPt>150) return false;
  }
  else if (RunNumber==108003) {
    if (hardestPt<150 || hardestPt>220) return false;
  }
  else if (RunNumber==108004) {
    if (hardestPt<220) return false;
  }
  
  return true;
}


void gammaJetAnalysis::FilterOutEvents(string outputFilename)  {
  //produce a new TTree, based on the fChain here, and fill it with a subset of events.
  bool debug=false;
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  if(debug) cout << "DEBUG " << nentries << endl << flush;

  TFile fout(outputFilename.c_str(),"recreate");  

  TTree* newTree;
  if(debug) cout << "DEBUG we made newTree" << endl << flush;

  newTree = fChain->CloneTree(0);
  if(debug) cout << "DEBUG we cloned fChain" << endl << flush;
  for (Long64_t e=0; e<nentries; e++) { 
    Long64_t ientry = LoadTree(e); 
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if(debug) cout << "DEBUG will check event " << e << endl << flush;
    if(!PassesEventFilter()) continue;
    newTree->Fill();
  }
  fout.cd();
  newTree->Write();
  fout.Close();
  return;
}




double gammaJetAnalysis::GetWeightOfEvent() {
  //this function assumes an event has been read in from the ntuple. 
  double answer = Weight;
  answer *= luminosity; //integrated luminosity in pb^{-1}, since we have assigned weighs assuming luminosity 1 pb^{-1}.
  answer *= ((double)prescale); //if we read every 3rd event, then give each event 3 times its normal weight.
  return answer;
}


double gammaJetAnalysis::GetJetRecoPt(int iJet) {
  assert (iJet>=0);
  assert (iJet<NJets);
  double pt=0;
  if (!useEMscaleAsRecoJetPt) pt = (*JetPt)[iJet];
  if (useEMscaleAsRecoJetPt)  pt = (*JetEMPt)[iJet];
  pt = pt/GeV;
  return pt;
}


void gammaJetAnalysis::FillHistogram(TH1F& hist, const double& value) {
  bool quantize=false;
  if (!quantize) {
    hist.Fill(value,GetWeightOfEvent());
    return;
  }
  if (quantize) {
    double expected = GetWeightOfEvent();
    int observed = randomEngine.Poisson(expected);
    for (int i = 0; i<observed; ++i) {
      hist.Fill(value,1);
    }
  }
}


void gammaJetAnalysis::FillHistogram(TH2F& hist, const double& valuex, const double& valuey) {
  bool quantize=false;
  if (!quantize) {
    hist.Fill(valuex,valuey,GetWeightOfEvent());
    return;
  }
  if (quantize) {
    double expected = GetWeightOfEvent();
    int observed = randomEngine.Poisson(expected);
    for (int i = 0; i<observed; ++i) {
      hist.Fill(valuex,valuey,1);
    }
  }
}


pair<int,double> gammaJetAnalysis::GetNearestGJToThisJG(string object, unsigned int iObject, string trueOrReco) {
  //to find the nearest true jet to a true photon N give ("photon", N , "true").
  //to find the nearest reco photon to a reco jet N give ("jet", N , "reco").
  //returns the index of the nearest guy and its energy, in the pair.
  bool debug=false;
 
  assert(object=="jet"||object=="photon");
  assert(trueOrReco=="true"||trueOrReco=="reco");

  if(debug) {
    if (object=="jet") cout << "- Looking for nearest " << trueOrReco << " photon to the " << object << " " << iObject << endl << flush;
    if (object=="photon") cout << "- Looking for nearest " << trueOrReco << " jet to the " << object << " " << iObject << endl << flush;
  }

  if (object=="photon") {
    if (trueOrReco=="true") assert (iObject<NPhot_true);
    if (trueOrReco=="reco") assert (iObject<NPhot);
  }
  if (object=="jet") {
    if (trueOrReco=="true") assert (iObject<NJets_true);
    if (trueOrReco=="reco") assert (iObject<NJets);
  }

  unsigned int numberOfTheOtherGuys=0;
  if (object=="photon") {
    if (trueOrReco=="true") numberOfTheOtherGuys=NJets_true;
    if (trueOrReco=="reco") numberOfTheOtherGuys=NJets;
  }  
  if (object=="jet") {
    if (trueOrReco=="true") numberOfTheOtherGuys=NPhot_true;
    if (trueOrReco=="reco") numberOfTheOtherGuys=NPhot;
  }  

  //check if there is any of the other guys at all of the kind we are looking for.
  if (numberOfTheOtherGuys==0) return (pair<int,double>(-1,-1));

  double ourGuyEta, ourGuyPhi;
  if (object=="photon") {
    if (trueOrReco=="true") {
      ourGuyEta=(*PhotEta_true)[iObject];
      ourGuyPhi=(*PhotPhi_true)[iObject];
    }
    if (trueOrReco=="reco") {
      ourGuyEta=(*PhotEta)[iObject];
      ourGuyPhi=(*PhotPhi)[iObject];
    }
  }
  if (object=="jet") {
    if (trueOrReco=="true") {
      ourGuyEta=(*JetEta_true)[iObject];
      ourGuyPhi=(*JetPhi_true)[iObject];
    }
    if (trueOrReco=="reco") {
      ourGuyEta=(*JetEta)[iObject];
      ourGuyPhi=(*JetPhi)[iObject];
    }
  }

  double DRmin = 1e6;
  double otherGuyEta=0;
  double otherGuyPhi=0;
  double dr=1e6;
  int nearestOtherGuy = -1;
  for (unsigned int iOtherGuy=0; iOtherGuy<numberOfTheOtherGuys; ++iOtherGuy) {
    if (object=="photon") {
      if (trueOrReco=="true") {
	otherGuyEta=(*JetEta_true)[iOtherGuy];
	otherGuyPhi=(*JetPhi_true)[iOtherGuy];
      }
      if (trueOrReco=="reco") {
	otherGuyEta=(*JetEta)[iOtherGuy];
	otherGuyPhi=(*JetPhi)[iOtherGuy];
      }
    }
    if (object=="jet") {
      if (trueOrReco=="true") {
	otherGuyEta=(*PhotEta_true)[iOtherGuy];
	otherGuyPhi=(*PhotPhi_true)[iOtherGuy];
      }
      if (trueOrReco=="reco") {
	otherGuyEta=(*PhotEta)[iOtherGuy];
	otherGuyPhi=(*PhotPhi)[iOtherGuy];
      }
    }

    dr = deltaR(ourGuyEta,otherGuyEta,ourGuyPhi,otherGuyPhi);
    if (dr < DRmin)  {
      DRmin = dr;
      nearestOtherGuy = iOtherGuy;
      if(debug) cout << "Nearest is: " << nearestOtherGuy << " at " << DRmin << endl << flush;
    }
  }
  return pair<int,double>(nearestOtherGuy,DRmin);
}




int gammaJetAnalysis::GetLeadingTightPhoton() {
  //tight here means ultra tight, namely being "good=tight" and isolated as we demand.
  //We want to loop through events that are tight and pass our cuts, and then return the leading among those.
  int leadingTightPhoton = -1;
  double ptmax = 0;
  for (unsigned int iPhot = 0; iPhot<NPhot; iPhot++) {
    double pt = (*PhotPt)[iPhot]/GeV;
    //make sure it's the kind of object we want, namely a tight photon.
    if ((*PhotGood)[iPhot]!=1) continue;
    double pEta = (*PhotEta)[iPhot];  
    bool photonInCrack = (fabs(pEta)>1.37 && fabs(pEta)<1.52);
    if (photonInCrack) continue ;
    if ((*PhotEtcone)[iPhot] > etConeCut*(*PhotPt)[iPhot]) continue;
    if ((*PhotNTrack)[iPhot] > nTracksCut) continue;
  
    if (pt > ptmax) {
      ptmax = pt;
      leadingTightPhoton = iPhot;
    }
  }
  return leadingTightPhoton;
}



bool gammaJetAnalysis::TrustJetIsReallyJet(unsigned int iJet,string trueOrReco) {
  //return 'true' if the jet does not coincide in DeltaR with any photon of significant energy.  If it coincides, then it's most likely a photon, which is double-written in both photon and jet containers, because it satisfies the definition of jet too, since the jet ID requirements are looser.
  assert (trueOrReco=="true"||trueOrReco=="reco");
  if (trueOrReco=="true") assert(iJet < NJets_true);
  if (trueOrReco=="reco") assert(iJet < NJets);

  double answer=true;

  pair<int,double> pairNearestPhotonAndDRmin = GetNearestGJToThisJG("jet",iJet,trueOrReco);
  double DRmin = pairNearestPhotonAndDRmin.second;
  if (DRmin>=0 && DRmin<0.4) {
    //it coincides with a photon in DeltaR
    //But does the photon have comparable momentum?
    int nearestPhoton = pairNearestPhotonAndDRmin.first;
    double photonPt=1;
    double jetPt=1;
    if (trueOrReco=="true") {
      jetPt=(*JetPt_true)[iJet] /GeV;
      photonPt=(*PhotPt_true)[nearestPhoton] /GeV;
    }
    if (trueOrReco=="reco") {
      //use always EM scale energy for this criterion
      jetPt=(*JetEMPt)[iJet] /GeV;
      photonPt=(*PhotPt)[nearestPhoton] /GeV;
    }
    assert(jetPt > 0);
    double photonOverJetPt = photonPt/jetPt;
    double threshold=0;
    if (trueOrReco=="reco") threshold = 0.55;
    if (trueOrReco=="true") threshold = 0.50;

    if (photonOverJetPt > threshold)
      answer = false;
  }
  
  return answer;
}



int gammaJetAnalysis::GetLeading(string object, string trueOrReco) {
  bool debug=false;
  if(debug) cout << "Looking for leading " << object << " among " << trueOrReco << endl << flush;
  assert(trueOrReco=="reco"||trueOrReco=="true");
  assert(object=="photon"||object=="jet");
  
  unsigned int Nobjects=0;
  if (object=="photon" && trueOrReco=="reco") Nobjects = NPhot;
  if (object=="jet" && trueOrReco=="reco")    Nobjects = NJets;
  if (object=="photon" && trueOrReco=="true") Nobjects = NPhot_true;
  if (object=="jet" && trueOrReco=="true")    Nobjects = NJets_true;

  if(debug) cout << "Nobjects= " << Nobjects << endl << flush;
	      
  int leadingGuy = -1;
  double ptmax = 0;
  for (unsigned int i = 0; i<Nobjects; i++) {
    double pt = 0;
    if (object=="photon" && trueOrReco=="reco") pt=(*PhotPt)[i]/GeV;
    if (object=="jet" && trueOrReco=="reco")    pt=GetJetRecoPt(i);
    if (object=="photon" && trueOrReco=="true") pt=(*PhotPt_true)[i]/GeV;
    if (object=="jet" && trueOrReco=="true")    pt=(*JetPt_true)[i]/GeV;

    if(debug) cout << "i= " << i << " pT= " << pt << endl << flush;

    //make sure it's the kind of object we want.  Photons are alse recorded as jets.
    if (object=="jet") {
      if (!TrustJetIsReallyJet(i,trueOrReco)) {
	if(debug) cout << "This 'jet' coincides with a photon " << endl << flush;
	continue; //it's actually a photon, not a jet.  Go to next jet
      }
    }

    if (pt > ptmax) {
      ptmax = pt;
      leadingGuy = i;
    }
    if(debug) cout << "So far leading " << object << " is " << leadingGuy << " with ptmax= " << ptmax << endl << flush;
  }
  return leadingGuy;
}


int gammaJetAnalysis::GetSubleading(string object, string trueOrReco) {
  assert(trueOrReco=="true"||trueOrReco=="reco");
  assert(object=="photon"||object=="jet");
  unsigned int Nobjects=0;
  if (object=="photon" && trueOrReco=="reco") Nobjects = NPhot;
  if (object=="jet" && trueOrReco=="reco")    Nobjects = NJets;
  if (object=="photon" && trueOrReco=="true") Nobjects = NPhot_true;
  if (object=="jet" && trueOrReco=="true")    Nobjects = NJets_true;

  int leadingGuy = GetLeading(object, trueOrReco);
  int subleadingGuy = -1;
  double ptmax = 0;
  double subleadingPt = 0;
  for (unsigned int i = 0; i<Nobjects; i++) {
    if (i == leadingGuy) continue;
    double pt = 0;
    if (object=="photon" && trueOrReco=="reco") pt=(*PhotPt)[i]/GeV;
    if (object=="jet" && trueOrReco=="reco")    pt=GetJetRecoPt(i);
    if (object=="photon" && trueOrReco=="true") pt=(*PhotPt_true)[i]/GeV;
    if (object=="jet" && trueOrReco=="true")    pt=(*JetPt_true)[i]/GeV;

    //make sure it's the kind of object we want.  Photons are alse recorded as jets.
    if (object=="jet")
      if (!TrustJetIsReallyJet(i,trueOrReco)) 
	continue; //it's actually a photon, not a jet.  Go to next jet.
    

    if (pt > ptmax) {
      ptmax = pt;
      subleadingGuy = i;
    }
  }
  return subleadingGuy;
}



int gammaJetAnalysis::GetMatchingTruePhotonToPhoton(unsigned int recoPhoton) {
  assert(recoPhoton<NPhot);
  int nearestTruePhoton = -1;
  double recoEta = (*PhotEta)[recoPhoton];
  double recoPhi = (*PhotPhi)[recoPhoton];

  double drMin = 1e6; //something big.
  for ( int i = 0; i < NPhot_true; ++i) {
    double trueEta = (*PhotEta_true)[i];
    double truePhi = (*PhotPhi_true)[i];
    double dr = deltaR(trueEta,recoEta,truePhi,recoPhi);
    if ( dr < drMin ) {
      drMin = dr;
      nearestTruePhoton = i;
    }
  }
  if (drMin < 0.4) 
    return nearestTruePhoton;
  else
    return (-1);
}


int gammaJetAnalysis::GetMatchingTrueJetToJet(unsigned int recoJet) {
  assert(recoJet<NJets);
  int nearestTrueJet = -1;
  double recoEta = (*JetEta)[recoJet];
  double recoPhi = (*JetPhi)[recoJet];

  double drMin = 1e6; //something big.
  for ( int i = 0; i < NJets_true; ++i) {
    double trueEta = (*JetEta_true)[i];
    double truePhi = (*JetPhi_true)[i];
    if(!TrustJetIsReallyJet(i,"true"))
      continue; //true "jet" i is actually a true photon, not a true jet.  Go to next true jet.
    double dr = deltaR(trueEta,recoEta,truePhi,recoPhi);
    if ( dr < drMin ) {
      drMin = dr;
      nearestTrueJet = i;
    }
  }
  if (drMin < 0.4) 
    return nearestTrueJet;
  else
    return (-1);
}


void gammaJetAnalysis::SetCuts(const double& etConeCut_, const int& nTracksCut_, const double& ptJ2cut_, const double& phiCut_) {
  assert (etConeCut_ >= 0);
  assert (nTracksCut_ >= 0);
  assert (ptJ2cut_ >= 0);
  assert (phiCut_ >=0 && phiCut_ <= pi);

  etConeCut=etConeCut_;
  nTracksCut=nTracksCut_;
  ptJ2cut=ptJ2cut_;
  phiCut=phiCut_;
  return;
}


void gammaJetAnalysis::CopyParametersFrom(gammaJetAnalysis* otherAnalysis) {
  SetCuts(otherAnalysis->etConeCut,otherAnalysis->nTracksCut,otherAnalysis->ptJ2cut,otherAnalysis->phiCut);
  luminosity = otherAnalysis->luminosity;
  prescale = otherAnalysis->prescale;
  useEMscaleAsRecoJetPt = otherAnalysis->useEMscaleAsRecoJetPt;
  return;
}


TH1F gammaJetAnalysis::CountTightPhotons() {
  TH1F hist("hist","",10,0,10);
  hist.Sumw2();
  hist.SetTitle(";Tight Photons in event;Events");

  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);

    int tightPhotons=0;
    for (int iPhot = 0 ; iPhot < NPhot; iPhot++) {
      if ((*PhotGood)[iPhot]!=1) continue;
      double pEta = (*PhotEta)[iPhot];  
      bool photonInCrack = (fabs(pEta)>1.37 && fabs(pEta)<1.52);
      if (photonInCrack) continue ;
      if ((*PhotEtcone)[iPhot] > etConeCut*(*PhotPt)[iPhot]) continue;
      if ((*PhotNTrack)[iPhot] > nTracksCut) continue;
      tightPhotons++;
    }
    hist.Fill(tightPhotons,GetWeightOfEvent());
  }
  return hist;
}


TH1F gammaJetAnalysis::GetHistoIsTheLeadingPhotonTheSameAsTheLeadingTightPhoton() {
  //tight here means ultra tight, namely being "good=tight" and isolated as we demand.
  TH1F hist("hist","",2,-1,1);
  hist.Sumw2();
  hist.SetTitle("How Often is #gamma_1 = leading tight #gamma?;Left bin are the events where #gamma_{1}#neq leading tight #gamma;Events");

  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);

    int leadingPhoton = GetLeading("photon","reco");
    int leadingTightPhoton = GetLeadingTightPhoton();
    if (leadingPhoton==-1 || leadingTightPhoton==-1) continue;

    if (leadingPhoton==leadingTightPhoton)  hist.Fill(0.5,GetWeightOfEvent());
    if (leadingPhoton!=leadingTightPhoton)  hist.Fill(-0.5,GetWeightOfEvent());
  }
  return hist;
}



bool gammaJetAnalysis::PassesEventCuts() {
  /* Apply the basic cuts, when we loop through events to fill histograms etc. */
  bool passes=false;

  //jet trigger
  /* I'm disabling this criterion, because the background 
    samples with JetFilter don't contain trigger information.
    For them, the Bitmaps are 0, so all events would fail. */
  //   bool passesJetTrigger = (EFJetBitmap != 0);
  //   bool passesPhotonTrigger = (EFPhBitmap != 0);
  //   bool passesTrigger = (passesJetTrigger || passesPhotonTrigger);
  //   if (!passesTrigger) return false;


  //require at least one one good (=tight) photon
  unsigned int goodPhotons = 0;
  for (unsigned int iPhoton=0; iPhoton<NPhot; iPhoton++)
    if ((*PhotGood)[iPhoton]==1) goodPhotons++;
  if (goodPhotons==0) return false;

  //require that the leading photon is a good (=tight) photon.
  int leadingPhoton = GetLeading("photon","reco");
  if (leadingPhoton == -1) return false; //there is no photon at all.
  if ((*PhotGood)[leadingPhoton]!=1) return false;

  //require that it is not in the crack
  double pEta = (*PhotEta)[leadingPhoton];  
  bool photonInCrack = (fabs(pEta)>1.37 && fabs(pEta)<1.52);
  if (photonInCrack) return false;

  //require that it's also isolated of energy around it.
  if ((*PhotEtcone)[leadingPhoton] > etConeCut*(*PhotPt)[leadingPhoton]) return false;
  
  //require that it's also isolated of tracks around it.
  if ((*PhotNTrack)[leadingPhoton]>nTracksCut) return false;

  //require at least one jet (that is not a photon)
  int leadingJet = GetLeading("jet","reco");
  bool hasLeadingJet = (leadingJet!=-1);
  if (!hasLeadingJet) return false;

  //require that the subleading jet (if any) isn't too strong
  //at the same time, apply the same pT cut on any sybleading photon (even loose; anything that is in the photon container)
  int subleadingJet = GetSubleading("jet","reco");
  int subleadingPhoton = GetSubleading("photon","reco");
  if (subleadingJet != -1)
    if (GetJetRecoPt(subleadingJet) > ptJ2cut) return false;
  if (subleadingPhoton != -1)
    if ((*PhotPt)[subleadingPhoton]/GeV > ptJ2cut) return false;

  //check the leading jet
  //eta cut
  double jEta = (*JetEta)[leadingJet];
  bool jetInCrack = ( fabs(jEta)>1.3 && fabs(jEta)<1.8 );
  if (jetInCrack) return false;

  bool jetWithinEtaCut = (fabs(jEta)<2.5);
  if (!jetWithinEtaCut) return false;

  double jPhi = (*JetPhi)[leadingJet];
  double pPhi = (*PhotPhi)[leadingPhoton];
  double DR   = deltaR(jEta,pEta,jPhi,pPhi);
  double Dphi = deltaPhi(jPhi,pPhi);
  bool backToBack = (Dphi > phiCut);
  if(!backToBack) return false;
  
  return true;
}



TH2F gammaJetAnalysis::GetJetPtVsTruePt() {
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  TH2F hist("hist","",12,limits,200,0,2);
  hist.Sumw2();
  hist.SetTitle(";p^{Jet,True}_{T} [GeV];p^{Jet,Reco}_{T}/p^{Jet,True}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingJet = GetLeading("jet","reco");
    double jetPtReco = GetJetRecoPt(leadingJet);
    
    int correspondingTrueJet = GetMatchingTrueJetToJet(leadingJet);
    if (correspondingTrueJet==-1) continue;
    double jetPtTrue = (*JetPt_true)[correspondingTrueJet] /GeV;
    hist.Fill(jetPtTrue,jetPtReco/jetPtTrue,GetWeightOfEvent());
  }
  return hist;
}



TH2F gammaJetAnalysis::GetPhotonPtVsPhotonTruePt() {
  Double_t limits[] = {0,10,20,30,40,50,60,70,80,90,100,120,140,160,180,200,250,300,350,400,450,500};
  TH2F hist("hist","",21,limits,200,0,2);
  hist.Sumw2();
  hist.SetTitle(";p^{#gamma,True}_{T} [GeV];p^{#gamma,Reco}_{T}/p^{#gamma,True}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    double photonPtReco = (*PhotPt)[leadingPhoton] /GeV;

    int matchingTruePhoton = GetMatchingTruePhotonToPhoton(leadingPhoton);
    double photonPtTrue = (*PhotPt_true)[matchingTruePhoton] /GeV;
    hist.Fill(photonPtTrue,photonPtReco/photonPtTrue,GetWeightOfEvent());
  }
  return hist;
}



TH1F  gammaJetAnalysis::GetHistoPtPhotonOverPtJetWhenTheTwoOverlap(string trueOrReco) {
  //this is useful to know where to cut, namely where the threshold should be in TrustJetIsReallyJet();
  assert(trueOrReco=="true"||trueOrReco=="reco");
  TH1F hist("hist","",100,0,2);
  hist.Sumw2();
  if (trueOrReco=="true") hist.SetTitle(";p^{#gamma,True}_{T}/p^{jet,True}_{T} when matching;jets");
  if (trueOrReco=="reco") hist.SetTitle(";p^{#gamma}_{T}/p^{jet}_{T} when matching;jets");

  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);

    unsigned int njets=0;
    if (trueOrReco=="true") njets=NJets_true;
    if (trueOrReco=="reco") njets=NJets;
    for (unsigned int iJet=0; iJet<njets ; ++iJet) {
      pair<int,double> pairPhotonDRmin = GetNearestGJToThisJG("jet",iJet,trueOrReco);
      double DRmin=pairPhotonDRmin.second;
      if (DRmin>=0 && DRmin<0.4) {
	double jetPt = 1;
	if (trueOrReco=="true") jetPt=(*JetPt_true)[iJet]/GeV;
	if (trueOrReco=="reco") jetPt=GetJetRecoPt(iJet);
	double photonPt = 1;
	if (trueOrReco=="true")	photonPt=(*PhotPt_true)[pairPhotonDRmin.first]/GeV;
	if (trueOrReco=="reco")	photonPt=(*PhotPt)[pairPhotonDRmin.first]/GeV;
	if (jetPt>0)
	  hist.Fill(photonPt/jetPt,GetWeightOfEvent());
      }
    }//next jet
  }//next event
  return hist;
}


TH1F  gammaJetAnalysis::GetHistoSpectrum(string quantity, string trueOrReco) {
  bool debug=false;
  assert (quantity=="pTaverage"||quantity=="j1Pt"||quantity=="gamma1Pt"||quantity=="hardestGammaOrJetPt"||quantity=="hardestInJetContainer");
  assert (trueOrReco=="true"||trueOrReco=="reco");
  TH1F hist("hist","",100,0,500);
  hist.Sumw2();
  if (quantity=="pTaverage") {
    if (trueOrReco=="true") hist.SetTitle(";p^{average,True}_{T} [GeV];Events");
    if (trueOrReco=="reco") hist.SetTitle(";p^{average}_{T} [GeV];Events");
  }
  if (quantity=="j1Pt") {
    if (trueOrReco=="true") hist.SetTitle(";p^{j1,True}_{T} [GeV];Events");
    if (trueOrReco=="reco") hist.SetTitle(";p^{j1}_{T} [GeV];Events");
  }
  if (quantity=="gamma1Pt") {
    if (trueOrReco=="true") hist.SetTitle(";p^{#gamma1,True}_{T} [GeV];Events");
    if (trueOrReco=="reco") hist.SetTitle(";p^{#gamma1}_{T} [GeV];Events");
  }
  if (quantity=="hardestGammaOrJetPt") {
    if (trueOrReco=="true") hist.SetTitle(";max(p^{#gamma1}_{T,True},p^{jet1}_{T,True}) [GeV];Events");
    if (trueOrReco=="reco") hist.SetTitle(";max(p^{#gamma1}_{T},p^{jet1}_{T}) [GeV];Events");
  }
  if (quantity=="hardestInJetContainer") {
    if (trueOrReco=="true") hist.SetTitle(";(*JetPt_true)[0] [GeV];Events");
    if (trueOrReco=="reco") hist.SetTitle(";(*JetPt)[0] [GeV];Events");
  }

  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    //    if (!PassesEventCuts()) continue;
    
    if(debug) cout << "-------- event " << e << " --------" << endl << flush;

    if (quantity=="pTaverage") {
      int leadingPhoton = -1;
      double photonPt = 0;
      int leadingJet = -1;
      double jetPt = 0;
      if (trueOrReco=="reco") {
	leadingPhoton = GetLeading("photon","reco");
	if (leadingPhoton==-1) continue;      
	photonPt = (*PhotPt)[leadingPhoton] /GeV;
	leadingJet = GetLeading("jet","reco");
	if (leadingJet==-1) continue;
	jetPt = GetJetRecoPt(leadingJet);
      }
      if (trueOrReco=="true") {
	leadingPhoton = GetLeading("photon","reco");
	if (leadingPhoton==-1) continue;      
	leadingPhoton=GetMatchingTruePhotonToPhoton(leadingPhoton);
	if (leadingPhoton==-1) continue;
	photonPt = (*PhotPt_true)[leadingPhoton] /GeV;
	leadingJet = GetLeading("jet","reco");
	if (leadingJet==-1) continue;
	leadingJet=GetMatchingTrueJetToJet(leadingJet);
	if (leadingJet==-1) continue;
	jetPt = (*JetPt_true)[leadingJet] /GeV;
      }
      double pTaverage =  (photonPt+jetPt)/2;
      hist.Fill(pTaverage,GetWeightOfEvent());
    }
    if (quantity=="gamma1Pt") {
      int leadingPhoton = -1;
      double photonPt = 0;
      if (trueOrReco=="reco") {
	leadingPhoton=GetLeading("photon","reco");
	if (leadingPhoton==-1) continue;
	photonPt = (*PhotPt)[leadingPhoton] /GeV;
      }
      if (trueOrReco=="true") {
	leadingPhoton=GetLeading("photon","true");
	if (leadingPhoton==-1) continue;
	photonPt = (*PhotPt_true)[leadingPhoton] /GeV;
      }
      hist.Fill(photonPt,GetWeightOfEvent());
    }
    if (quantity=="j1Pt") {
      int leadingJet = -1;
      double jetPt = 0;
      if (trueOrReco=="reco") {
	leadingJet=GetLeading("jet","reco");
	if (leadingJet==-1) continue;
	jetPt = GetJetRecoPt(leadingJet);
      }
      if (trueOrReco=="true") {
	leadingJet=GetLeading("jet","true");
	if (leadingJet==-1) continue;
	jetPt = (*JetPt_true)[leadingJet] /GeV;
      }
      hist.Fill(jetPt,GetWeightOfEvent());
    }
    if (quantity=="hardestGammaOrJetPt") {
      int leadingPhoton = -1;
      double photonPt = 0;
      int leadingJet = -1;
      double jetPt = 0;
      if (trueOrReco=="reco") {
	leadingPhoton = GetLeading("photon","reco");
	leadingJet = GetLeading("jet","reco");
	if (leadingJet == -1 && leadingPhoton == -1) continue;
	if (leadingJet == -1)    
	  jetPt = -1;
	else
	  jetPt = GetJetRecoPt(leadingJet);
	if (leadingPhoton == -1) 
	  photonPt = -1;
	else
	  photonPt = (*PhotPt)[leadingPhoton] /GeV;
      }
      hist.Fill(max(photonPt,jetPt),GetWeightOfEvent());
    }
    if (quantity=="hardestInJetContainer") {
      if (trueOrReco=="true") 
	hist.Fill((*JetPt_true)[0]/GeV,GetWeightOfEvent()); 
      if (trueOrReco=="reco") 
	hist.Fill((*JetPt)[0]/GeV,GetWeightOfEvent()); 
    }
  }//next event
  return hist;
}


TH2F  gammaJetAnalysis::GetHistoJetPtVsGammaPt(string trueOrReco) {
  assert(trueOrReco=="true"||trueOrReco=="reco");
  TH2F hist("hist","",100,0,200,100,0,200);
  hist.Sumw2();
  if (trueOrReco=="true") hist.SetTitle(";p^{#gamma,True}_{T} [GeV];p^{Jet,True}_{T} [GeV];Events");
  if (trueOrReco=="reco") hist.SetTitle(";p^{#gamma}_{T} [GeV];p^{Jet}_{T} [GeV];Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = -1;
    leadingPhoton = GetLeading("photon","reco");
    if (leadingPhoton==-1) continue;
    if (trueOrReco=="true") leadingPhoton=GetMatchingTruePhotonToPhoton(leadingPhoton);


    int leadingJet = -1;
    leadingJet = GetLeading("jet", "reco");
    if (leadingJet==-1) continue;
    if (trueOrReco=="true") leadingJet=GetMatchingTrueJetToJet(leadingJet);
	
    double photonPt = 0;
    if (trueOrReco=="reco") photonPt=(*PhotPt)[leadingPhoton] /GeV;
    if (trueOrReco=="true") photonPt=(*PhotPt_true)[leadingPhoton] /GeV;
    
    double jetPt = 0;
    if (trueOrReco=="reco") jetPt=GetJetRecoPt(leadingJet);
    if (trueOrReco=="true") jetPt=(*JetPt_true)[leadingJet] /GeV;

    hist.Fill(photonPt,jetPt,GetWeightOfEvent());
  }
  return hist;
}



TH1F  gammaJetAnalysis::GetHistoBalance() {
  TH1F hist("hist","",200,-1,1);
  hist.Sumw2();
  hist.SetTitle(";(p^{Jet}_{T}-p^{#gamma}_{T})/(p^{Jet}_{T}+p^{#gamma}_{T});events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    int leadingJet = GetLeading("jet","reco");
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double balance = (jetPt - photonPt)/(jetPt + photonPt);
    hist.Fill(balance,GetWeightOfEvent());
  }
  return hist;
}


TH2F  gammaJetAnalysis::GetHistoJet2PtVsBalance() {
  TH2F hist("hist","",200,-1,1,100,0,200);
  hist.Sumw2();
  hist.SetTitle(";(p^{Jet}_{T}-p^{#gamma}_{T})/(p^{Jet}_{T}+p^{#gamma}_{T});p^{jet2}_{T} [GeV];Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    int leadingJet = GetLeading("jet","reco");
    int jet2 = GetSubleading("jet","reco");
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double balance = (jetPt - photonPt)/(jetPt + photonPt);
    double jet2Pt = 0;
    if (jet2 != -1) jet2Pt = GetJetRecoPt(jet2);
    hist.Fill(balance,jet2Pt,GetWeightOfEvent());
  }
  return hist;
}


TH2F    gammaJetAnalysis::GetHistoJESvsJetPt(string option) {
  assert (option=="mcWeights"||option=="pseudodata");
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  TH2F hist("hist","",12,limits,300,0,3);
  hist.Sumw2();
  hist.SetTitle(";p^{Jet}_{T} [GeV];p^{#gamma}_{T}/p^{Jet}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    int leadingJet = GetLeading("jet","reco");
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double correctionFactor = photonPt/jetPt;
    hist.Fill(jetPt,correctionFactor,GetWeightOfEvent());
  }
  if (option=="mcWeights") {
    return hist;
  }
  else { //if option=="pseudodata"
    //simulate data (quantize)
    TH2F hist2 = hist;
    hist2.Reset();
    int observedEntries = randomEngine.Poisson(hist.Integral());
    hist2.FillRandom(&hist,observedEntries);
    return hist2;
  }
}


TH2F    gammaJetAnalysis::GetHistoJESvsJetPt_true() {
  Double_t limits[] = {0,10,20,30,40,50,60,70,80,90,100,120,140,160,180,200,250,300,350,400,450,500};
  TH2F hist("hist","",21,limits,300,0,3);
  hist.Sumw2();
  hist.SetTitle(";p^{Jet,True}_{T} [GeV];p^{#gamma,True}_{T}/p^{Jet,True}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    leadingPhoton = GetMatchingTruePhotonToPhoton(leadingPhoton);
    int leadingJet = GetLeading("jet","reco");
    leadingJet = GetMatchingTrueJetToJet(leadingJet);

    double photonPt = (*PhotPt_true)[leadingPhoton] /GeV; 
    double jetPt = (*JetPt_true)[leadingJet] /GeV;
    double correctionFactor = photonPt/jetPt;
    hist.Fill(jetPt,correctionFactor,GetWeightOfEvent());
  }
  return hist;
}


TH2F    gammaJetAnalysis::GetHistoPhotonPtOverJetPtVsAveragePt(string option) {
  assert (option=="mcWeights"||option=="pseudodata");
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  TH2F hist("hist","",12,limits,75,0,3);
  hist.Sumw2();
  hist.SetTitle(";p^{average}_{T} [GeV];p^{#gamma}_{T}/p^{Jet}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    int leadingJet = GetLeading("jet","reco");
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double averagePt = (jetPt + photonPt)/2;
    double ratio = photonPt/jetPt;
    hist.Fill(averagePt,ratio,GetWeightOfEvent());
  }
  if (option=="mcWeights") {
    return hist;
  }
  else { //if option=="pseudodata"
    //simulate data (quantize)
    TH2F hist2 = hist;
    hist2.Reset();
    int observedEntries = randomEngine.Poisson(hist.Integral());
    hist2.FillRandom(&hist,observedEntries);
    return hist2;
  }
}


TH2F    gammaJetAnalysis::GetHistoGRecoOverJTrueVsAveragePt() {
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  TH2F hist("histGRecoOverJTrueVsAveragePt","",12,limits,300,0,3);
  hist.Sumw2();
  hist.SetTitle(";p^{average}_{T} [GeV];p^{#gamma}_{T}/p^{Jet,True}_{T};Events");
  //  hist.SetMarkerStyle(1);
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    if (leadingPhoton==-1) continue;
    int leadingJet = GetLeading("jet","reco");
    if (leadingJet==-1) continue;
    int correspondingTrueJet = GetMatchingTrueJetToJet(leadingJet);
    if (correspondingTrueJet==-1) continue;
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double jetPt_true = (*JetPt_true)[correspondingTrueJet] /GeV;
    double ratio = photonPt/jetPt_true;
    double averagePt = (jetPt + photonPt)/2;

    double jes = photonPt/jetPt;
    if ( jes>(1-2*2/sqrt(jetPt)) && jes<(1+2*2/sqrt(jetPt)) ) //mostly events with JES within 2 sigma of resolution have a role in determining the JES, due to the way the iterative fitting is done.  If it's out in the tails, the fitter doesn't use those events.  So, don't  fill the histogram with those events.  In this way, when we later find the systematic error, we will limit its definition to the fraction *of events that play a role in the JES* which are not in good balance at semi-true level.  To define this cut, we assumed that sigmaJES = JES*sigmaPtJetReco/PtJetReco, where sigmaPtJetReco/PtJetReco is roughly 2/sqrt(pTjetReco).
      hist.Fill(averagePt,ratio,GetWeightOfEvent());
  }
  return hist;
}


void gammaJetAnalysis::GetBothHistosGoverJandJESVsAverage(TH2F& histJESVsAverage, TH2F& histActualBalance, string option) {
  //this is a merger of GetHistoPhotonPtOverJetPtVsAveragePt() and GetHistoGRecoOverJTrueVsAveragePt(), to gain speed by looping through events once instead of twice.
  assert (option=="mcWeights"||option=="pseudodata");
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  histJESVsAverage = TH2F("histJES","",12,limits,75,0,3);
  histJESVsAverage.Sumw2();
  histJESVsAverage.SetTitle(";p^{average}_{T} [GeV];p^{#gamma}_{T}/p^{Jet}_{T};Events");
  //  histJESVsAverage.SetMarkerStyle(1);
  
  histActualBalance = TH2F("histGRecoOverJTrueVsAveragePt","",12,limits,300,0,3);
  histActualBalance.Sumw2();
  histActualBalance.SetTitle(";p^{average}_{T} [GeV];p^{#gamma}_{T}/p^{Jet,True}_{T};Events");
  //  histActualBalance.SetMarkerStyle(1);


  if (fChain == 0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    if (leadingPhoton==-1) continue;
    int leadingJet = GetLeading("jet","reco");
    if (leadingJet==-1) continue;
    
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double averagePt = (jetPt + photonPt)/2;
    double jes = photonPt/jetPt;
    histJESVsAverage.Fill(averagePt,jes,GetWeightOfEvent());

    if ( !(jes>(1-2*2/sqrt(jetPt)) && jes<(1+2*2/sqrt(jetPt))) ) continue; //mostly events within 2 sigma of resolution have a role in determining the JES, due to the way the iterative fitting is done.  If it's out in the tails, the fitter doesn't use those events.  So, don't  fill the histogram with those events.  In this way, when we later find the systematic error, we will limit its definition to the fraction *of events that play a role in the JES* which are not in good balance at semi-true level.  To define this cut, we assumed that sigmaJES = JES*sigmaPtJetReco/PtJetReco, where sigmaPtJetReco/PtJetReco is roughly 2/sqrt(pTjetReco).
      
    int correspondingTrueJet = GetMatchingTrueJetToJet(leadingJet);
    if (correspondingTrueJet==-1) continue;
    double jetPt_true = (*JetPt_true)[correspondingTrueJet] /GeV;
    double ratio = photonPt/jetPt_true;
    histActualBalance.Fill(averagePt,ratio,GetWeightOfEvent());
  }

  if (option=="pseudodata") {
    //  simulate data (quantize) for histJESVsAverage
    TH2F hist2 = histJESVsAverage;
    hist2.Reset();
    int observedEntries = randomEngine.Poisson(histJESVsAverage.Integral());
    hist2.FillRandom(&histJESVsAverage,observedEntries);
    histJESVsAverage.Reset();
    histJESVsAverage=hist2;
  }

  return;
}



void gammaJetAnalysis::GetBothHistosGoverJandJESVsJetPt(TH2F& histJESVsJetPt, TH2F& histActualBalance, string option) {
  //this is a merger of GetHistoPhotonPtOverJetPtVsAveragePt() and GetHistoGRecoOverJTrueVsAveragePt(), to gain speed by looping through events once instead of twice.
  assert (option=="mcWeights"||option=="pseudodata");
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  histJESVsJetPt = TH2F("histJES","",12,limits,75,0,3);
  histJESVsJetPt.Sumw2();
  histJESVsJetPt.SetTitle(";p^{Jet}_{T} [GeV];p^{#gamma}_{T}/p^{Jet}_{T};Events");
  //  histJESVsJetPt.SetMarkerStyle(1);
  
  histActualBalance = TH2F("histGRecoOverJTrueVsJetPt","",12,limits,300,0,3);
  histActualBalance.Sumw2();
  histActualBalance.SetTitle(";p^{Jet}_{T} [GeV];p^{#gamma}_{T}/p^{Jet,True}_{T};Events");
  //  histActualBalance.SetMarkerStyle(1);


  if (fChain == 0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    if (!PassesEventCuts()) continue;

    int leadingPhoton = GetLeading("photon","reco");
    if (leadingPhoton==-1) continue;
    int leadingJet = GetLeading("jet","reco");
    if (leadingJet==-1) continue;
    
    double photonPt = (*PhotPt)[leadingPhoton] /GeV; 
    double jetPt = GetJetRecoPt(leadingJet);
    double averagePt = (jetPt + photonPt)/2;
    double jes = photonPt/jetPt;
    histJESVsJetPt.Fill(jetPt,jes,GetWeightOfEvent());

    //    if ( !(jes>(1-2*2/sqrt(jetPt)) && jes<(1+2*2/sqrt(jetPt))) ) continue; //mostly events within 2 sigma of resolution have a role in determining the JES, due to the way the iterative fitting is done.  If it's out in the tails, the fitter doesn't use those events.  So, don't  fill the histogram with those events.  In this way, when we later find the systematic error, we will limit its definition to the fraction *of events that play a role in the JES* which are not in good balance at semi-true level.  To define this cut, we assumed that sigmaJES = JES*sigmaPtJetReco/PtJetReco, where sigmaPtJetReco/PtJetReco is roughly 2/sqrt(pTjetReco).
      
    int correspondingTrueJet = GetMatchingTrueJetToJet(leadingJet);
    if (correspondingTrueJet==-1) continue;
    double jetPt_true = (*JetPt_true)[correspondingTrueJet] /GeV;
    double ratio = photonPt/jetPt_true;
    histActualBalance.Fill(jetPt,ratio,GetWeightOfEvent());
  }

  if (option=="pseudodata") {
    //  simulate data (quantize) for histJESVsJetPt
    TH2F hist2 = histJESVsJetPt;
    hist2.Reset();
    int observedEntries = randomEngine.Poisson(histJESVsJetPt.Integral());
    hist2.FillRandom(&histJESVsJetPt,observedEntries);
    histJESVsJetPt.Reset();
    histJESVsJetPt=hist2;
  }

  return;
}




/////////////////////////////////////////////////////////////////////////////////////
//QCD test functions, used to examine the effects of JES correction to the jet response.
/////////////////////////////////////////////////////////////////////////////////////

TH2F   gammaJetAnalysis::GetHistoResponseToJets(TGraphErrors* graphJESVsJetPt) {
  Double_t limits[] = {0,10,20,40,60,80,100,150,200,250,300,400,500};
  TH2F hist("hist","",12,limits,200,0,2);
  hist.Sumw2();
  hist.SetTitle(";p^{Jet,True}_{T} [GeV];p^{Jet,Reco}_{T}/p^{Jet,True}_{T};Events");
  
  if (fChain == 0) return hist;
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t e=0; e<nentries; e += prescale) { 
    Long64_t ientry = LoadTree(e); //ientry is the entry number in the sub-tree of the chain where the e^th event is found.
    if (ientry < 0) break;
    fChain->GetEntry(e);
    //    if (!PassesEventCuts()) continue;
    
    for (unsigned int iJet=0; iJet<NJets; iJet++) {
      double jetEta  = (*JetEta)[iJet];
      if (fabs(jetEta)>=2.5) continue; //we don't derive JES for forward jets and we don't want to examine how they behave.
      //check if it's not a jet but a photon
      if (!TrustJetIsReallyJet(iJet,"reco"))
	continue; //it's actually a photon, not a jet.  Go to next jet
      double jetPtReco = GetJetRecoPt(iJet);
      if (graphJESVsJetPt != NULL) {
	//apply the JES correction
	//find which point corresponds to the given jetPtReco before correction.
	if (jetPtReco < 500) {
	  double lowEdge=0;
	  double pTjetAtCenterOfBin=0;
	  double highEdge=0;
	  double deltaPt=0;
	  double jesCorrection=1;
	  for (int point=0; point<graphJESVsJetPt->GetN(); ++point) {
	    graphJESVsJetPt->GetPoint(point,pTjetAtCenterOfBin,jesCorrection);
	    deltaPt = graphJESVsJetPt->GetErrorX(point);
	    lowEdge=pTjetAtCenterOfBin - deltaPt;
	    highEdge=pTjetAtCenterOfBin + deltaPt;
	    if (jetPtReco >= lowEdge && jetPtReco < highEdge)
	      break;
	  }
	  jetPtReco *= jesCorrection;
	}
      }
      int correspondingTrueJet = GetMatchingTrueJetToJet(iJet);
      if (correspondingTrueJet==-1) continue;
      double jetPtTrue = (*JetPt_true)[correspondingTrueJet] /GeV;
      hist.Fill(jetPtTrue,jetPtReco/jetPtTrue,GetWeightOfEvent());
    }
  }
  return hist;  
}