TPoll/src/RecoFirst.cc

// -*- C++ -*-
//
// Package:    RecoFirst
// Class:      RecoFirst
// 
/**\class RecoFirst RecoFirst.cc recoFirst/RecoFirst/src/RecoFirst.cc

 Description: Displaced vertex method for long lived reconstruction

 Implementation:
    From a MC sample find a RECO electron or photon
    Find a matching MC lepton
    See if MC lepton is a (grand+)daughter of signal (Z or Z', say)
    Construct displaced vertex and hence reconstructed mass
    Apply iterative algorithm to refine solution
 
   This analysis works with both RECO and AOD data
   The reconstructed particles are generically referred to as RECO
   To distinguish reconstructed particles vs MC
*/
//
// Original Author:  Tony Poll
//         Created:  Thu Jul  7 12:48:08 BST 2011
// $Id$
//
//


// system include files
#include <memory>
#include <typeinfo>

// user include files

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h" 

#include "DataFormats/EgammaCandidates/interface/Electron.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"

#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronCore.h"
#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
#include <DataFormats/CaloRecHit/interface/CaloCluster.h> 
#include <DataFormats/EgammaReco/interface/BasicCluster.h> 
//#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <DataFormats/Common/interface/Ref.h> 
#include <DataFormats/DetId/interface/DetId.h> 

#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
#include <DataFormats/EgammaReco/interface/BasicCluster.h> 
#include <DataFormats/CaloRecHit/interface/CaloRecHit.h> 

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
 
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"
#include "DataFormats/EcalRawData/interface/EcalDCCHeaderBlock.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/Candidate/interface/Candidate.h"

#include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
 
#include "CaloOnlineTools/EcalTools/interface/EcalFedMap.h"

#include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
#include <CommonTools/Utils/interface/Angle.h>

#include "RecoFirst.h"
#include "TVector3.h"
#include "TFile.h"
#include "TTree.h"
#include "TText.h"
#include "TGraphAsymmErrors.h"

#include "CLHEP/Geometry/Transform3D.h"


//
// class declaration
//
/*
class RecoFirst : public edm::EDAnalyzer {
   public:
      explicit RecoFirst(const edm::ParameterSet&);
      ~RecoFirst();

      static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);


   private:
      virtual void beginJob() ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      virtual void beginRun(edm::Run const&, edm::EventSetup const&);
      virtual void endRun(edm::Run const&, edm::EventSetup const&);
      virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
      virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);

      // ----------member data ---------------------------
};
*/

//
// constants, enums and typedefs
//
        using namespace edm;
        using namespace std;

        edm::Handle<EcalRecHitCollection> EBhits_;
        edm::Handle<EcalRecHitCollection> EEhits_;

        edm::Handle<vector<reco::GsfElectron> > ElectronCollection;
        edm::Handle<vector<reco::Photon> > PhotonCollection;

        float speedOfLight_ = 30.0; // speed  of light in units of cm/ns
        float decayLengthAccuracy = 0.01; // difference between one iteration & the next when calculating decay length (%)
    float timeAccuracy = 0.02;  // Iterate until differences in measured and calculated flight time are less than this (ns)

    const int pdgIdElectron = 11;
        const int pdgIdMuon = 13;
        const int pdgIdTau = 15;
        const int pdgIdPhoton = 22;

    // Number of each long lived daughter decay types
    int numberOfElectrons = 0;
    int numberOfPhotons = 0;
    int numberOfElectronPhotonPairs = 0;
    int numberOfDaughters = 0;
    int lowCrystalHitEnergy = 0;
    int negativeHitTime = 0;
    int noDisplacedVertexConvergence = 0;
    int ValidatedReconstructions = 0;
    int nonPhysicalDecayLength = 0;
    int closeElectronsWithNegativeHitTime = 0;
    int negativeFlightLength = 0;

    int numberOfPrimaryVertexEvents = 0;
    int ValidatedEvents = 0;
    int noParticles = 0;
    int numberOfZParticles = 0;
    int numberOfZeventsAboveMinimumMass =0;
    int numberOfZeventsWithTwoMcDaughters = 0;
    int numberMcDaughterElectronsNotFinalState = 0;
    int foundDaughters = 0;
    int zDecaysWithEnergeticDaughters = 0;
    int badRecoDeltaRMatch = 0;
    int goodRecoElectron = 0;
    int goodRecoPhoton = 0;
        int noRecoMatch = 0;
    int nullSuperClusters = 0;
    int zeroEnergyLepton = 0;
    int zWithTwoRecoDaughters = 0;
    int numberOfZsWithoutTwoPositiveTimes = 0;
    int numberOfDisplacedVertices = 0;
    int lookingForLeptons = 0;
    int particleHits = 0;
    int recoPairsGoingIntoVeto = 0;
    int recoPairsPassedVeto = 0;
    int failedIsolation = 0;
    int failedBarrelEndcapGap =0;
    int failedEtaCut = 0;
    int failedMinEt = 0;
    int massReconstructions = 0;
    int middleReconstructedMass = 0;
    int massCalcCounter = 0;
    int countValidations = 0;

    // Array for Error Plots
    int ErrorsCount = 0;
    const int n = 10000;
    float bestMass [n];
    float bestFlightLength [n];
    float lowMass [n];
    float lowFlightLength [n];
    float highMass [n];
    float highFlightLength [n];

    // Values for measured - error, measured, and measured + error values, 
    int low = 0;
    int medium = 1;
    int high = 2;
        
    int ElectronDecays = 0;
    int bremElectrons = 0;
    int MuonDecays = 0;
    int TauDecays = 0;
    int PhotonDecays = 0;
    int minimumAngleBetweenElectronsRejects = 0;

    float crystalMinEnergy = 4.0;      // Below 4 Gev hit time has too much error
        // float crystalMinEnergy = 0.0;      // Below 4 Gev hit time has too much error
    float photonMinPt = 60.0;          // Trigger, as of ~ 18/05/2011, rejects diphotons with Pt < 60GeV
    // float minimumZMass = 30.0;         // Don't consider off-shell MC Z events
    float minimumZMass = 0.0;         // Don't consider Pythia MC (gamma* + Z) 'particles'
    // float minimumZMass = 80.0;         // Don't consider Pythia MC (gamma* + Z) 'particles'
    // float minimumAngleBetweenElectrons = 0.05; // Very small angles between electrons gives reconstructed mass ~ zero
    float minimumAngleBetweenElectrons = 0.1; // Very small angles between electrons gives reconstructed mass ~ zero
    float bremFraction = 0.2;   // A high bremFraction means ECAL hit energy << parent energy, which results in poor mass reconstruction

    float deltaConstructedMass = 0.10;  // Iterate until changes are small
    // float deltaConstructedMass = 0.01;  // Iterate until changes are small

    float dummy = -10; 

  // Photon selection cuts
    float minimumMatchDeltaR = 0.1;     // Reject MC to reco match if difference is greater than this value
    // float maxEta = 50.0;             // Photon eta cut
    float maxEta = 2.5;                 // Photon eta cut
    float isolationMinimum = 0.0;       // isolation: deltaR between two photons
    // float isolationMinimum = 0.45;   // isolation: deltaR between two photons
    // float photonMinEt = 1.0;         // Electron/Photon cut
    float photonMinEt = 30.0;           // Electron/Photon cut
    // float minDeltaR = 1.0;              // MC to Reco electron deltaR
    float minDeltaR = 0.5;              // MC to Reco electron deltaR
    float minDeltaOmega = 0.01;         // Loop displaced vertex calc until change in angle betweeen electrons is small

    float barrelLimit = 1.44;
    float endcapLimit = 1.57;

    float barrelTimingError = 0.27;     //ns from CMS Note 2010/012
    float endcapTimingError = 0.18;

        math::XYZPoint PrimaryVertex_;
    float AngleBetweenHits;

        float pi = acos (-1);

        // For histograms
        map<string, TH1*> hists_;
        map<string, TH2*> hists2d_;
    map<string, TGraphAsymmErrors*> TGraphAsymmErrors_;

    //bool minimum = true;
    bool minimum = false;
    bool measured = true;
    // bool measured = false;
    // bool maximum = true;
    bool maximum = false;

    // bool Z-ee = true;
    bool Zee = false;
    // bool ZPrimeee = true;
    bool ZPrimeee = false;
    bool LLee = true;
    // bool LLee = false;

        bool Debug1 = true;
        // bool Debug1 = false;
        // bool Debug2 = true;
        bool Debug2 = false;

//
// static data member definitions
//

//
// constructors and destructor
//
RecoFirst::RecoFirst(const edm::ParameterSet& iConfig):
   //now do what ever initialization is needed
        PhotonSrc_                      (iConfig.getParameter<InputTag>("Photons")),
        electronSrc_        (iConfig.getParameter<InputTag>("gsfElectrons")),
        EBRecHitCollection_ (iConfig.getParameter<InputTag>("EcalRecHitCollectionEB")),
        EERecHitCollection_ (iConfig.getParameter<InputTag>("EcalRecHitCollectionEE")),
        generatorTag_ (iConfig.getParameter<InputTag>("generatorSrc")),
        signalPDGId_ (iConfig.getParameter<int>("signalPDGId"))
{
}


RecoFirst::~RecoFirst()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called for each event  ------------
void
RecoFirst::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;

    if(Debug2){cout << "Starting to read event" << endl;}

        // Get Ecal hits for this event
        iEvent.getByLabel(EBRecHitCollection_, EBhits_);
    if (EBhits_.failedToGet())
    {
        std::cout << "WARNING: cannot access EBhits_" << std::endl;
        return;
    }

        iEvent.getByLabel(EERecHitCollection_, EEhits_);
    if (EEhits_.failedToGet())
    {
        std::cout << "WARNING: cannot access EEhits_" << std::endl;
        return;
    }
        if(Debug2){cout << "EBhits: " << EBhits_->size() << ", EEhits: " << EEhits_->size() << endl;}

    // Get electron collection
        iEvent.getByLabel(electronSrc_, ElectronCollection);    // Indirectly pointed to in recofirst_cfi.py
    if (ElectronCollection.failedToGet())
    {
        std::cout << "WARNING: cannot access ElectronCollection" << std::endl;
        return;
    }
        if(Debug2) {cout << "Number of electrons: " << ElectronCollection->size() << endl;}
    // Get photon collection
        iEvent.getByLabel(PhotonSrc_, PhotonCollection);        // Indirectly pointed to in recofirst_cfi.py
      if (PhotonCollection.failedToGet())
    {
        std::cout << "WARNING: cannot access PhotonCollection" << std::endl;
        return;
    }
    if(Debug2) {cout << "Number of photons: " << PhotonCollection->size() << endl;}

    // Get primary vertex
    getPrimaryVertex(iEvent);

    // Confirm algorithms accurately reconstructed parent flight length and mass
    // *************************************************************************
    // By comparing RECO particle reconstructions (mass & flight length) to 'true' MC particle values
    // If the two RECO electrons are close to the MC electrons (i.e. deltaR is small)
    // Plot (reconstructed mass vs true MC mass) and (reconstructed flight length vs true MC flight length)
    
    std::vector<reco::GenParticle> goodMcElectrons;

    // Create a particle, with initial dummy values, to hold the Monte Calro signal 
    myParticle McSignal (dummy);
       
    // Get MC particles
        edm::Handle<View<reco::GenParticle> > Particles;
        iEvent.getByLabel(generatorTag_, Particles);    // 'generatorTag_' defined in recofirst_cfi.py
        if (Particles.failedToGet())
        { 
        noParticles++;
                if (Debug2) {cout << "No MC particles in this event" << endl; }
                // return;  // Do not exit as real physics data will not have MC particles
        }
    else    // Only do this if the sample contains MC particles
    {
        // Match RECO electrons/photons to MC electrons that are descendants of our signal (Z, Z' or LL)
        // Rather than attempt to match to all the MC particles every time, which would be slow,
        // create a collection of MC electrons that are descendants of our signal
        bool signalFound = false;
        for (unsigned imc = 0; imc < Particles->size(); ++imc)  // Get the MC particles
        { 
            // Looking for MC electron/positron in final state i.e. Pythia status == 1
            if((abs((*Particles)[imc].pdgId()) == pdgIdElectron) && ((*Particles)[imc].status()==1))
            {
                if(Debug2){cout << "Got a MC final state electron" << endl;}
                reco::GenParticle decayParticle;
                decayParticle = (*Particles)[imc];
                // Now navigate up the decay chain looking for our signal particle
                for (unsigned i = 0; i < decayParticle.numberOfMothers(); ++i)
                {
                    const reco::Candidate* mother = decayParticle.mother(i);
                    // 'signalPDGId_' defined in recofirst_cfi.py
                    while ((mother->pdgId() != signalPDGId_) && (mother->numberOfMothers() != 0 ))
                    {
                        mother = mother->mother(0);
                    }
                    // if ((mother->pdgId() == signalPDGId_) && (mother->status() == 3))
                    if ((mother->pdgId() == signalPDGId_))
                    {
                        signalFound = true;
                        goodMcElectrons.push_back(decayParticle);  // Store good MC electrons
                        McSignal.mass[medium] = mother->p4().M();
                        McSignal.energy = mother->energy();
                        if(Debug2){cout << "Signal mass: " << McSignal.mass[medium] << endl;}
                        McSignal.flightLength[medium] = getFlightLength(mother);
                       // Plot MC mass & flight length
                        hists_["t6"]->Fill(McSignal.mass[medium]);      // True (MC) mass
                        hists_["t3"]->Fill(McSignal.flightLength[medium]);      // True (MC) flight length
                    }
                    else if(Debug2) {cout << "Got a bad signal MC mother: " << i << " with PDGId: " << mother->pdgId() << endl;}
                }
            }
        }
    }

    // Create an instance for the two daughter particles and the reconstructed particle
    // Initialised with dummy values: usually -10, as this is easily spotted when not set with calculated values
    myParticle firstDaughterRECO (dummy);
    myParticle secondDaughterRECO (dummy);
    myParticle reconstructedRECOParent(dummy);
    
    // Get some statistics: number of electrons, photons, those without a supercluster, low ECAL crystal energy, negative ECAL hit time etc
    if(Debug2){cout << "Number of electrons: " << ElectronCollection->size() << endl;}
    if(Debug2){cout << "Number of photons: " << PhotonCollection->size() << endl;}
    for(std::vector<reco::GsfElectron>::const_iterator firstElectron = ElectronCollection->begin(); firstElectron != ElectronCollection->end(); ++firstElectron)
    {
        getDaughterStats(firstDaughterRECO, firstElectron);
    }
    for(std::vector<reco::Photon>::const_iterator firstPhoton = PhotonCollection->begin(); firstPhoton != PhotonCollection->end(); ++firstPhoton)
    {   
        getDaughterStats(firstDaughterRECO, firstPhoton);    
    }
    
    // Search for resonances
    // Find all electron/photon pairs and reconstruct an invariant mass
    // Find the firstDaughterRECO electron
    for(std::vector<reco::GsfElectron>::const_iterator firstElectron = ElectronCollection->begin(); firstElectron != ElectronCollection->end(); ++firstElectron)
    {
        numberOfElectrons++;
        getDaughterParameters(firstDaughterRECO, firstElectron);
        if(Debug2){cout << "Electron firstDaughterRECO.measuredHitTime: " << firstDaughterRECO.measuredHitTime[medium] << endl;}
        // For an electron/electron pair find a second electron
        for(std::vector<reco::GsfElectron>::const_iterator secondElectron = firstElectron + 1; secondElectron != ElectronCollection->end(); ++secondElectron)
        {   
            numberOfElectronPhotonPairs++;
            getDaughterParameters(secondDaughterRECO, secondElectron);
            if(Debug2){cout << "Electron secondDaughterRECO.measuredHitTime: " << secondDaughterRECO.measuredHitTime[medium] << endl;}
            reconstructMass(reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);
            validate(McSignal, goodMcElectrons, reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);
        }
        //  For an electron/photon pair find a photon
        for(std::vector<reco::Photon>::const_iterator secondPhoton = PhotonCollection->begin(); secondPhoton != PhotonCollection->end(); ++secondPhoton)
        {   
            numberOfElectronPhotonPairs++;
            getDaughterParameters(secondDaughterRECO, secondPhoton);
            if(Debug2){cout << "Photon secondDaughterRECO.measuredHitTime: " << secondDaughterRECO.measuredHitTime[medium] << endl;}
            reconstructMass(reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);    
            validate(McSignal, goodMcElectrons, reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);               
        }
    }
 
    // Find photon/photon pairs
    // Find the firstDaughterRECO photon
    for(std::vector<reco::Photon>::const_iterator firstPhoton = PhotonCollection->begin(); firstPhoton != PhotonCollection->end(); ++firstPhoton)
    {   
        numberOfPhotons++;
        getDaughterParameters(firstDaughterRECO, firstPhoton);
        if(Debug2){cout << "Photon firstDaughterRECO.measuredHitTime: " << firstDaughterRECO.measuredHitTime[medium] << endl;}
        //  For a photon/photon pair find the secondDaughterRECO photon
        for(std::vector<reco::Photon>::const_iterator secondPhoton = firstPhoton + 1; secondPhoton != PhotonCollection->end(); ++secondPhoton)
        {
            numberOfElectronPhotonPairs++;
            getDaughterParameters(secondDaughterRECO, secondPhoton);
            if(Debug2){cout << "Photon secondDaughterRECO.measuredHitTime: " << secondDaughterRECO.measuredHitTime[medium] << endl;}
            reconstructMass(reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);
            validate(McSignal, goodMcElectrons, reconstructedRECOParent, firstDaughterRECO, secondDaughterRECO);
        }
    }
}


// Compare calculated values against true MC flight length and mass
// Check both daughters come from the same parent (In LL samples there will be two LL decays)
// Check each daughter is close (i.e. small deltaR) to MC electron or photon
void 
RecoFirst::validate(myParticle McSignal, std::vector<reco::GenParticle> McElectrons, myParticle reconstructedRECOParent, myParticle firstDaughterRECO, myParticle secondDaughterRECO)
{
    countValidations++;

    if(reconstructedRECOParent.mass[medium] != dummy) // If mass is not reconstructed it retains its initialisation negative value
    {
        
        // If MC electrons are close (small deltaR) to RECO electrons, then 
        // plot MC vs reconstructedRECOParent mass
        // plot MC flight length vs reconstructedRECOParent flight length
        if(Debug2){cout << "McSignalMass: " << McSignal.mass[medium] << endl;}
        // Note: An electron that decays from the signal may decay into a subsequent electron, while emitting brem radiation.
        // Therefore the best match electron may or may not be either the initial or final electron in the decay chain
        for (std::vector<reco::GenParticle>::const_iterator McElectron1 = McElectrons.begin(); McElectron1 != McElectrons.end(); McElectron1++)
        {
            float parentPt1 = parentPt(McElectron1->mother());  // Navigate up the decay chain until the signal is found, then get its Pt
            if (parentPt1 == dummy)    continue;   // Only proceed if we have a real Pt
            for (std::vector<reco::GenParticle>::const_iterator McElectron2 = McElectron1 + 1; McElectron2 != McElectrons.end(); McElectron2++)
            {
                // Only proceed is McElectron2 is not a daughter of McElectron1, or McElectron1 is not a daughter of McElectron2
                // We want two electrons from separate electron decay chains (this is different from two electrons from the same signal parent).
                // Note: You must dereference the McElectron as it is a pointing into a vector (which itself is an array of pointers)
                if (sameDecayChain(*McElectron1, *McElectron2)) continue;
                // Only proceed if the two MC particles share the same parent, as confirmed by its pt
                // if (Debug1){cout << "typeid(McElectron1).name(): " << typeid(McElectron1).name() << endl;}
                if (parentPt1 != parentPt(McElectron2->mother())) continue;   // Same parent?
                    
                if (Debug2) {cout << "Got two MC electrons with the same signal parent" << endl;}
                // Is the first RECO daughter close to one MC electron
                // and the second RECO daughter close to the other MC electron?
                // float deltaRValue = deltaR(McElectron1, firstDaughterRECO);
                float deltaPhi1 = calcDeltaPhi(McElectron1->phi(), firstDaughterRECO.phi);
                float deltaEta1 = calcDeltaEta(McElectron1->eta(), firstDaughterRECO.eta);
                float deltaPhi2 = calcDeltaPhi(McElectron2->phi(), secondDaughterRECO.phi);
                float deltaEta2 = calcDeltaEta(McElectron2->eta(), secondDaughterRECO.eta);
                
                float deltaPhi3 = calcDeltaPhi(McElectron2->phi(), firstDaughterRECO.phi);
                float deltaEta3 = calcDeltaEta(McElectron2->eta(), firstDaughterRECO.eta);
                float deltaPhi4 = calcDeltaPhi(McElectron1->phi(), secondDaughterRECO.phi);
                float deltaEta4 = calcDeltaEta(McElectron1->eta(), secondDaughterRECO.eta);

                // Only proceed if both RECO electrons are close to different MC electrons
                if ( !(  (  // firstDaughterRECO close to McElectron1 ?
                            (sqrt ((deltaPhi1 * deltaPhi1) + (deltaEta1 * deltaEta1)) < minDeltaR) &&
                            // secondDaughterRECO close to McElectron2 ?
                            (sqrt ((deltaPhi2 * deltaPhi2) + (deltaEta2 * deltaEta2)) < minDeltaR) )
                    ||      // OR
                        (   // firstDaughterRECO close to McElectron2 ?
                            (sqrt ((deltaPhi3 * deltaPhi3) + (deltaEta3 * deltaEta3)) < minDeltaR) &&
                            // second DaughterRECO close to McElectron1 ?
                            (sqrt ((deltaPhi4 * deltaPhi4) + (deltaEta4 * deltaEta4)) < minDeltaR) )
                    )) continue;    // Skip loop if the two RECO electrons don't match separate MC electrons
                {
                    ValidatedReconstructions++;
                    if ((firstDaughterRECO.measuredHitTime[medium] != dummy) &&  // Only plot measured hit times
                        (secondDaughterRECO.measuredHitTime[medium] != dummy))
                    {
                        closeElectronsWithNegativeHitTime++;
                    }
                    if (reconstructedRECOParent.flightLength[medium] < 0.0)
                    {
                        negativeFlightLength++;
                    }

                    // Find (true) angle between the two MC electrons
                    reco::GenParticle::Vector McParticleMom1 = McElectron1->momentum();
                    reco::GenParticle::Vector McParticleMom2 = McElectron2->momentum();
                    double McIncludedAngle = acos(McParticleMom1.Dot(McParticleMom2)/McParticleMom1.R()/McParticleMom2.R());
                    if (Debug2) {cout << "McIncludedAngle: " << McIncludedAngle << endl;}

                    // Calculate the expected MC hit time. I.e. LL flight time + electron flight time
                    float McSignalGamma = McSignal.energy / McSignal.mass[medium];
                    float McSignalVelocity = speedOfLight_ * sqrt( 1 - 1/(McSignalGamma * McSignalGamma));
                    if (Debug1){cout << "McSignalVelocity: " << McSignalVelocity << endl;}
                    float McSignalFlightTime = McSignal.flightLength[medium] / McSignalVelocity;
                    if (Debug1){cout << "McSignalFlightTime: " << McSignalFlightTime << endl;}
                    float McElectronFlightTime1 = sqrt( (McParticleMom1.X() * McParticleMom1.X()) + (McParticleMom1.Y() * McParticleMom1.Y()) + (McParticleMom1.Z() * McParticleMom1.Z()) )/speedOfLight_;
                    float McElectronFlightTime2 = sqrt( (McParticleMom2.X() * McParticleMom2.X()) + (McParticleMom2.Y() * McParticleMom2.Y()) + (McParticleMom2.Z() * McParticleMom2.Z()) )/speedOfLight_;
                    McSignal.calculatedHitTime[1] = McSignalFlightTime + McElectronFlightTime1;   // MC LL flight time + MC electron1 flight time
                    McSignal.calculatedHitTime[2] = McSignalFlightTime + McElectronFlightTime2;   // MC LL flight time + MC electron2 flight time
                    
                    hists_["t35"]->Fill(McSignalFlightTime);
                    hists_["t36"]->Fill(McElectronFlightTime1);
                    hists_["t36"]->Fill(McElectronFlightTime2);

                    // We need the time relative to a photon travelling at velocity from the primary vertex
                    // Find the XYZ coords of the hit for each electron
                    float x1 = McElectron1->vx() + McElectron1->px();
                    float y1 = McElectron1->vy() + McElectron1->py();
                    float z1 = McElectron1->vz() + McElectron1->pz();
                    float x2 = McElectron2->vx() + McElectron2->px();
                    float y2 = McElectron2->vy() + McElectron2->py();
                    float z2 = McElectron2->vz() + McElectron2->pz();

                    // Time for a photon to travel from IP to hit position
                    float time1 = (sqrt(x1*x1 + y1*y1 + z1*z1))/speedOfLight_;
                    float time2 = (sqrt(x2*x2 + y2*y2 + z2*z2))/speedOfLight_;
                    hists_["t37"]->Fill(time1);
                    hists_["t37"]->Fill(time2);

                    // Plot the MC vs measured times (relative to time=0 for a relativistic particle)
                    // Want to plot hit times for correct MC & RECO electron pairs
                    // firstDaughterRECO close to McElectron1 ?
                    if (sqrt ((deltaPhi1 * deltaPhi1) + (deltaEta1 * deltaEta1)) < minDeltaR)
                    {
                        if (Debug1){cout << "First pair match" << endl;}
                        hists_["t33"]->Fill(McSignal.calculatedHitTime[1] - time1, firstDaughterRECO.measuredHitTime[medium]);
                        hists_["t33"]->Fill(McSignal.calculatedHitTime[2] - time2, secondDaughterRECO.measuredHitTime[medium]);
                    }
                    else
                    {
                        if (Debug1){cout << "Second pair match" << endl;}
                        hists_["t33"]->Fill(McSignal.calculatedHitTime[2] - time2, firstDaughterRECO.measuredHitTime[medium]);
                        hists_["t33"]->Fill(McSignal.calculatedHitTime[1] - time1, secondDaughterRECO.measuredHitTime[medium]);
                    }

                    // Plot actual hit time i.e. don't subtract time for photon to travel from IP to ECAL hit
                    float x = firstDaughterRECO.endPoint.x();
                        float y = firstDaughterRECO.endPoint.y();
                        float z = firstDaughterRECO.endPoint.z();
                    float timeAdjustment = sqrt(x*x + y*y + z*z)/speedOfLight_;
                    firstDaughterRECO.calculatedHitTime[medium] = firstDaughterRECO.measuredHitTime[medium] + timeAdjustment;

                    x = secondDaughterRECO.endPoint.x();
                        y = secondDaughterRECO.endPoint.y();
                        z = secondDaughterRECO.endPoint.z();
                    timeAdjustment = sqrt(x*x + y*y + z*z)/speedOfLight_;
                    secondDaughterRECO.calculatedHitTime[medium] = secondDaughterRECO.measuredHitTime[medium] + timeAdjustment;

                    // Plot the MC vs measured absolute times (No adjustment for time=0 for a relativistic particle)
                    // Want to plot hit times for correct MC & RECO electron pairs
                    // firstDaughterRECO close to McElectron1 ?
                    if (sqrt ((deltaPhi1 * deltaPhi1) + (deltaEta1 * deltaEta1)) < minDeltaR)
                    {
                        if (Debug1){cout << "First pair match" << endl;}
                        hists_["t34"]->Fill(McSignal.calculatedHitTime[1], firstDaughterRECO.calculatedHitTime[medium]);
                        hists_["t34"]->Fill(McSignal.calculatedHitTime[2], secondDaughterRECO.calculatedHitTime[medium]);
                    }
                    else
                    {
                        if (Debug1){cout << "Second pair match" << endl;}
                        hists_["t34"]->Fill(McSignal.calculatedHitTime[2], firstDaughterRECO.calculatedHitTime[medium]);
                        hists_["t34"]->Fill(McSignal.calculatedHitTime[1], secondDaughterRECO.calculatedHitTime[medium]);
                    }


                    if (Debug2) {cout << "McSignal flight Length: " << McSignal.flightLength[medium] << endl;}

                    //If both RECO electrons are close to MC electrons, then compare true vs reconstructed mass & flight length
                    if ((reconstructedRECOParent.flightLength[medium] >= 0.0) // flight length < 0 is non-physical
                        && (firstDaughterRECO.measuredHitTime[medium] != dummy) // Only plot measured hit times
                        && (secondDaughterRECO.measuredHitTime[medium] != dummy)
                        && (reconstructedRECOParent.omega[medium] > minimumAngleBetweenElectrons ) // Large errors if angle between electrons is small
                        && ((McSignal.energy - (firstDaughterRECO.energy + secondDaughterRECO.energy))/McSignal.energy < (bremFraction) ) // Reject events with high brem
                        )    
                    {
                        hists_["t4"]->Fill(McSignal.flightLength[medium], reconstructedRECOParent.flightLength[medium]);   // True (MC) mass vs constructed flight length
                        hists_["t12"]->Fill(McSignal.mass[medium], reconstructedRECOParent.mass[medium]);   // True (MC) mass vs constructed mass
                        hists_["t13"]->Fill(reconstructedRECOParent.flightLength[medium] - McSignal.flightLength[medium], reconstructedRECOParent.mass[medium] - McSignal.mass[medium]);   // Reconstructed flight length residual vs Reconstructed mass residual
                        hists_["t14"]->Fill(reconstructedRECOParent.flightLength[medium]);   // Reconstructed flight length for 'good' electron pairs

                        // if (firstDaughterRECO.measuredHitTime[medium] != dummy)  // Only plot measured hit times
                        {
                            hists_["t15"]->Fill(firstDaughterRECO.measuredHitTime[medium]);   // ECAL hit time for 'good' electron 
                        }
                        // if (secondDaughterRECO.measuredHitTime[medium] != dummy)  // Only plot measured hit times
                        {
                            hists_["t15"]->Fill(secondDaughterRECO.measuredHitTime[medium]);   // ECAL hit time for 'good' electron 
                        }
                        hists_["t16"]->Fill(McSignal.mass[medium], reconstructedRECOParent.mass[medium] - McSignal.mass[medium]);   // Reconstructed flight length residual vs Reconstructed mass residual
                        hists_["t17"]->Fill(reconstructedRECOParent.mass[medium]);   // Reconstructed mass for electron hits matched to two MC electrons with same signal parent
                        hists_["t18"]->Fill(reconstructedRECOParent.flightLength[medium]);   // Reconstructed flight length for electron hits matched to two MC electrons with same signal parent
                        hists_["t19"]->Fill(McSignal.flightLength[medium], reconstructedRECOParent.flightLength[medium]);   // MC vs Reconstructed flight length for electron hits matched to two MC electrons with same signal parent
                        hists_["t23"]->Fill(reconstructedRECOParent.flightLength[medium], reconstructedRECOParent.omega[medium]);     // How included angle depends upon flight length
                        hists_["t24"]->Fill(McSignal.flightLength[medium], reconstructedRECOParent.flightLength[medium] - McSignal.flightLength[medium]); // How constructed flight length error depends upon MC flight length
                        hists_["t25"]->Fill(reconstructedRECOParent.omega[medium], reconstructedRECOParent.mass[medium]); // How reconstructed mass depends upon included angle
                        hists_["t26"]->Fill(sqrt(firstDaughterRECO.energy * secondDaughterRECO.energy), reconstructedRECOParent.mass[medium]); // How reconstructed mass depends upon energy of daughters
                        hists_["t27"]->Fill((firstDaughterRECO.energy * secondDaughterRECO.energy), (1-cos(reconstructedRECOParent.omega[medium]))); // How constructed included angle depends upon energy of daughters
                        hists_["t28"]->Fill(McIncludedAngle, McSignal.mass[medium]); // MC (true) included angle vs MC (true) mass
                        hists_["t29"]->Fill(McIncludedAngle, reconstructedRECOParent.omega[medium]); //  MC (true) included angle vs reconstructed included angle
                        hists_["t30"]->Fill(McSignal.flightLength[medium], McIncludedAngle); //  MC (true) flight length vs MC (true) included angle
                        hists_["t31"]->Fill(McIncludedAngle); // MC (true) included angle
                        hists_["t32"]->Fill(reconstructedRECOParent.omega[medium]); //  Constructed included angle

                        if (Debug2)
                        {
                            // If true mass >> reconstructed mass - to highlight serious mis-constructions
                            if (McSignal.mass[medium]/reconstructedRECOParent.mass[medium] > 1.4)
                            {   
                                cout << "Reconstructed MC mass ratio: " << McSignal.mass[medium]/reconstructedRECOParent.mass[medium] << endl;
                                cout << "MC signal parent mass : " << McSignal.mass[medium] << endl;
                                cout << "Reconstructed MC mass: " << reconstructedRECOParent.mass[medium] << endl;
                            }
                        }
                    }
                    return;
                }
            }
        }
    }
}


// Return the Pt of a McParticle's parent
float
RecoFirst::parentPt(const reco::Candidate* McParticle)
{
    float McParticlePt = dummy;
    while ((McParticle->pdgId() != signalPDGId_) && (McParticle->numberOfMothers() != 0))
    {
        McParticle = McParticle->mother(0);
    }
    if (McParticle->pdgId() == signalPDGId_)
    {
        McParticlePt = McParticle->pt();
    }
    if (Debug2){cout << "McParticle->pt: " << McParticlePt << endl;}
    return McParticlePt;
}

// True if the electrons are from separate decay chains
// Compare electrons on their Pt.
bool 
RecoFirst::sameDecayChain(const reco::GenParticle& McElectron1, const reco::GenParticle& McElectron2)
{
    if (McElectron1.pt() == McElectron2.pt()) return true;

    const reco::GenParticle& startMcElectron1 = McElectron1;
    const reco::GenParticle& startMcElectron2 = McElectron2;

    while ((McElectron1.pdgId() != signalPDGId_) && (McElectron1.numberOfMothers() != 0))
    {
        McElectron1 = McElectron1.mother(0);
        if (McElectron1.pt() == startMcElectron2.pt()) return true;

    }

    while ((McElectron2.pdgId() != signalPDGId_) && (McElectron2.numberOfMothers() != 0))
        {
            McElectron2 = McElectron2.mother(0);
            if (McElectron2.pt() == startMcElectron1.pt()) return true;
        }
    }

    return false;
}


// Primary Vertex setup
void 
RecoFirst::getPrimaryVertex(const edm::Event& iEvent)
{
    // Get primary vertex
    edm::Handle<reco::VertexCollection> primaryVertex;
    iEvent.getByLabel("offlinePrimaryVertices", primaryVertex);
    if(primaryVertex.failedToGet())
    {
        std::cout << "WARNING: cannot access PhotonCollection" << std::endl;
        return;
    }
    
    numberOfPrimaryVertexEvents++;
    if(Debug2) {cout << "Number of primary vertice(s): " << primaryVertex->size() << endl;}
    // If more than one primary use the first
    // TODO check which primary vertex should be used. Currently using the first which may be (?) the primary vertex with highest Pt
    reco::Vertex pv = primaryVertex->front();
    PrimaryVertex_ = pv.position();
}

float 
RecoFirst::getFlightLength(const reco::Candidate* mother)
{
    float flightLength = dummy; // A negative value that can be used for debugging, if required.
    float dx=100000,dy=100000,dz=100000;
    // We have the mother,
        dx=mother->vertex().X();
        dy=mother->vertex().Y();
        dz=mother->vertex().Z();
    if(Debug2){cout << "mother dx: " << dx << ", dy: " << dy << ", dz: " << dz << endl;}

        // Calculate the x, y & z distances between the mother & daughter vertices
    // If ( (particle is Pythia status == 3, or (not an electron or photon) ), and has a daughter: then get the daughter
    while ( ( (mother->status() == 3 ) || !( (abs(mother->pdgId()) == pdgIdElectron) ||  (mother->pdgId() == pdgIdPhoton) ) ) 
        &&  (mother->numberOfDaughters() != 0) )
    {
        mother = mother->daughter(0);
    }
    dx-=mother->vertex().X();
        dy-=mother->vertex().Y();
        dz-=mother->vertex().Z();
    if(Debug2){cout << "mother - daughter dx: " << dx << ", dy: " << dy << ", dz: " << dz << endl;}    
    flightLength = sqrt(dx*dx + dy*dy + dz*dz);
    if(Debug2){cout << "MC flight length: " << flightLength << endl;}
    return flightLength;
}

// Reconstruct the long lived particle's invariant mass from the 2 daughter particles
// This is an iterative algorithm
// Find an initial displaced vertex, and hence an initial value for LL mass assuming the LL particle's velocity = speed of light
// Using the LL mass and energy, calculate its gamma value, and hence velocity
// With the new LL velocity and the displaced vertex calculate the time of flight of the LL.
// With the electron path length and velocity c calculate the time of flight of the lectron.
// Compare the LL plus electron flight times with the measured ECAL hit time
// Now iterate on LL flight length (displaced vertex) and direction until the calculated hit time is close to that measured
// Each time the displaced vertex changes the electron path also changes to retain the same hit position on the ECAL
void 
RecoFirst::reconstructMass(myParticle &Parent, myParticle &firstDaughterRECO, myParticle &secondDaughterRECO)   // pass particles 'by reference' so parameters are updated
{
    if ((firstDaughterRECO.measuredHitTime[medium] == dummy) || (secondDaughterRECO.measuredHitTime[medium] == dummy))
    {
        if(Debug2){cout << "One or both electrons do not have a good ECAL hit time" << endl;}
        return; // Exit if either hit time has not been set
    }
    Parent.flightLength[medium] = dummy;  // Re-set values each time, else it retains previous values
    Parent.mass[medium] = dummy;
    // Only reconstruct mass for particles above the minimum Et
    if ((firstDaughterRECO.et < photonMinEt) || (secondDaughterRECO.et < photonMinEt) )
    {
        if(Debug2){cout << "Low daughter Et" << endl;}
        return;  // Exit
    }

    massReconstructions++;
    if(Debug2){cout << "Reconstructing mass" << endl;}
    // Find the displaced vertex for the two particles

    // Cannot find a displaced vertex if the timed flight length for (LL + daughter) < daughter direct flight length for either daughter
    if((firstDaughterRECO.timedDistanceToHit[medium] < firstDaughterRECO.PrimaryVertexToHit) || (secondDaughterRECO.timedDistanceToHit[medium] < secondDaughterRECO.PrimaryVertexToHit))
    {
        if(Debug2){cout << "firstDaughterRECO.timedDistanceToHit[" << 1 << "]: " << firstDaughterRECO.timedDistanceToHit[medium] << ", firstDaughterRECO.PrimaryVertexToHit: " << firstDaughterRECO.PrimaryVertexToHit << endl
                << "secondDaughterRECO.timedDistanceToHit[" << 1 << "]: " << secondDaughterRECO.timedDistanceToHit[medium] << ", secondDaughterRECO.PrimaryVertexToHit: " << secondDaughterRECO.PrimaryVertexToHit << endl 
                << "Cannot construct displaced vertex" << endl;}
                    
        numberOfZsWithoutTwoPositiveTimes++;
        return; // Exit without making any plots
    }
    // else: calculate an initial value for the displaced vertex (i.e. the parent particle's flight length)
        
    if(Debug2){cout << "Looking for displaced vertex" << endl;}
    Parent.flightLength[medium] = findDisplacedVertexLength(firstDaughterRECO, 1, secondDaughterRECO, 1);
    if (Parent.flightLength[medium] == dummy)    // displaced vertex algorithm never converged.
    {
        return;     // Without setting any values for reconstructed mass
    }
    numberOfDisplacedVertices++;
    if(Debug2){cout << "Initial flight length: " <<  Parent.flightLength[medium] << endl;}
        
    // Find the parent's mass, energy & hence gamma
    // m = sqrt[2*E(1) * E(2)*(1-cos(theta)] where 'theta' is the included angle between the 2 electron hits
    // firstDaughterRECO.omega: Angle between firstDaughterRECO electron hit, primary vertex and displaced vertex
    // secondDaughterRECO.omega:  Angle between secondDaughterRECO electron hit, primary vertex and displaced vertex
    // if(Debug2){cout << "firstDaughterRECO omega[" << 1 << "]: " << firstDaughterRECO.omega[medium] << ", secondDaughterRECO omega[" << 1 << "]: " << secondDaughterRECO.omega[medium] << endl;}

    if(Debug2){cout << "firstDaughterRECO.flightLength[medium]: " << firstDaughterRECO.flightLength[medium] << ", Parent.flightLength[medium]: " << Parent.flightLength[medium] << ",firstDaughterRECO.PrimaryVertexToHit: " << firstDaughterRECO.PrimaryVertexToHit << endl;}
    if(Debug2){cout << "secondDaughterRECO.flightLength[medium]: " << secondDaughterRECO.flightLength[medium] << ", Parent.flightLength[medium]: " << Parent.flightLength[medium] << ",secondDaughterRECO.PrimaryVertexToHit: " << secondDaughterRECO.PrimaryVertexToHit << endl;}
 
    // angleBetweenHits = 2*pi - (angle between PV->DV->Hit1) - (angle between PV->DV->Hit2)

    float angleBetweenPrimaryVertexDisplacedVertexHitOne = getAngle(firstDaughterRECO.flightLength[medium],Parent.flightLength[medium],firstDaughterRECO.PrimaryVertexToHit);
    float angleBetweenPrimaryVertexDisplacedVertexHitTwo = getAngle(secondDaughterRECO.flightLength[medium],Parent.flightLength[medium],secondDaughterRECO.PrimaryVertexToHit);
    if(Debug2){cout << "initial angleBetweenPrimaryVertexDisplacedVertexHitOne: " << angleBetweenPrimaryVertexDisplacedVertexHitOne << endl;}
    if(Debug2){cout << "initial angleBetweenPrimaryVertexDisplacedVertexHitTwo: " << angleBetweenPrimaryVertexDisplacedVertexHitTwo << endl;}

    float AngleBetweenHits = (2*pi) - angleBetweenPrimaryVertexDisplacedVertexHitOne -angleBetweenPrimaryVertexDisplacedVertexHitTwo;  
    if(Debug2){cout << "Initial AngleBetweenHits: " << AngleBetweenHits << endl;}
    // Reconstructed mass is ~ zero for very small angles
    if (AngleBetweenHits < minimumAngleBetweenElectrons)
    {
        minimumAngleBetweenElectronsRejects++;
        return;
    }
    if (Debug2) {cout << "firstDaughterRECO.energy: " << firstDaughterRECO.energy << ", secondDaughterRECO.energy: " << secondDaughterRECO.energy << endl;}
    
    Parent.mass[medium] = sqrt(2.0 * firstDaughterRECO.energy * secondDaughterRECO.energy * (1 - cos(AngleBetweenHits)));
    Parent.energy = firstDaughterRECO.energy + secondDaughterRECO.energy;
    if (Debug2) {cout << "Parent.mass[medium]: " << Parent.mass[medium] << ", Parent.energy: " << Parent.energy << endl;}

    
    // Parent.gamma[medium] = Parent.energy / Parent.mass[medium];
  
    // Now iterate: modify the constructed parent's flight length to allow for the current calculated mass of the constructed parent
    // The angles between the hits will initially remain constant as we modify the displaced vertex length (for the mass of the LL) and electron path length
    
    // firstDaughterRECO.flightLength[medium] = firstDaughterRECO.timedDistanceToHit[medium] - Parent.flightLength[medium];
    // secondDaughterRECO.flightLength[medium] = secondDaughterRECO.timedDistanceToHit[medium] - Parent.flightLength[medium];
    
    int count = 0;
    float oldConstructedParentMass = dummy; // arbitrary initial value
    float deltaTime = 1.0;  // arbitrary initial value
    float firstDeltaTime;
    float secondDeltaTime;
    // float deltaOmega; // Incremental change to omega for both daughters

    // Calculate initial value for reconstructed mass
    // Initial approximation assumed a LL mass of zero and velocity of c
    // Parent.mass[medium] = sqrt(2.0 * firstDaughterRECO.energy * secondDaughterRECO.energy * (1 - cos(AngleBetweenHits)));
    
    if(Debug2){cout << endl << "Finding new reconstructed mass" << endl;}
    // while ( (fabs(oldConstructedParentMass - Parent.mass[medium] ) > deltaConstructedMass)  && (count < 100))
    float deltaOmega = (firstDaughterRECO.omega[medium] + secondDaughterRECO.omega[medium])/5;// Only a small change required

    while ( (fabs(deltaTime ) > 0.2)  && (count < 100))
    {
        oldConstructedParentMass = Parent.mass[medium];
        count++;
        // Now modify result for a LL mass > 0, and LL velocity < c
        // Use the LL energy & mass to calculate its gamma, and hence velocity
        if (Parent.mass[medium] != 0)
        {
            Parent.gamma[medium] = Parent.energy / Parent.mass[medium];
        }
        else
        {
            Parent.gamma[medium] = 1000; // arbitrary large value
        }
        
        if (Debug2) {cout << "count: " << count << ", LLGamma: " << Parent.gamma[medium] << endl;}
        // Find the LLvelocity
        if (Parent.gamma[medium] < 1000)
        {
            Parent.velocity = speedOfLight_*sqrt(1 - 1/(Parent.gamma[medium] * Parent.gamma[medium]));
            if(Debug2){cout << "Constructed parent velocity: " << Parent.velocity << endl;}
        }
        else
        {  
            Parent.velocity = speedOfLight_; // If gamma = infinity, velocity = speedOfLight_
        }
        // Using the LL velocity calculate the expected flight times and compare with the measured time
        // Calculate the hit time for each electron, being constructed parent flight time + electron flight time
        // For firstDaughterRECO electron
        float firstConstructedTimeOfFlight = (Parent.flightLength[medium]/Parent.velocity) + (firstDaughterRECO.flightLength[medium]/speedOfLight_);
        if(Debug2){cout << "firstConstructedTimeOfFlight: " << firstConstructedTimeOfFlight << ", firstDaughterRECO.measuredHitTime[medium]: " << (firstDaughterRECO.PrimaryVertexToHit/speedOfLight_) + firstDaughterRECO.measuredHitTime[medium] << endl;}
        float secondConstructedTimeOfFlight = (Parent.flightLength[medium]/Parent.velocity) + (secondDaughterRECO.flightLength[medium]/speedOfLight_);
        if(Debug2){cout << "secondConstructedTimeOfFlight: " << secondConstructedTimeOfFlight << ",secondDaughterRECO.measuredHitTime[medium]: " << (secondDaughterRECO.PrimaryVertexToHit/speedOfLight_) + secondDaughterRECO.measuredHitTime[medium] << endl;}
        
        // Find difference in calculated and measured flight times
        // Note - the measured time is the excess over the time a relativistic particle would take.
        firstDeltaTime = firstConstructedTimeOfFlight - (firstDaughterRECO.PrimaryVertexToHit/speedOfLight_ + firstDaughterRECO.measuredHitTime[medium]);
        secondDeltaTime = secondConstructedTimeOfFlight - (secondDaughterRECO.PrimaryVertexToHit/speedOfLight_ + secondDaughterRECO.measuredHitTime[medium]);
        // deltaTime = fabs(firstDeltaTime) + fabs(secondDeltaTime);
        deltaTime = fabs(firstDeltaTime) + fabs(secondDeltaTime);
        if(Debug2){cout << "firstDeltaTime: " << firstDeltaTime << ", secondDeltaTime: " << secondDeltaTime << ", deltaTime: " << deltaTime << endl;}
        
        // Modify displaced vertex length based upon the average deltaTime for both electrons
        if ((firstDeltaTime * secondDeltaTime) > 0)     // If delta times are of the same sign
        {
            // A contrieved algorithm to modify parent flight length. Maybe a better algorithm could be used
            if(Debug2){cout << "firstDaughterRECO.measuredHitTime[medium]: " << firstDaughterRECO.measuredHitTime[medium] << endl;}
            if(Debug2){cout << "secondDaughterRECO.measuredHitTime[medium]: " <<secondDaughterRECO.measuredHitTime[medium] << endl;}
            if(Debug2){cout << "Old Parent.flightLength[medium]: " << Parent.flightLength[medium] << endl;}

            // Adjust parent flight length by ratio of (difference of measured time) to (total flight time)
            // Factor of "2*" gives faster convergence
            // If the calculated time > measured, then the slow parent flight length must be shorter, and the fast daughters' length longer
            // Parent.flightLength[medium] = Parent.flightLength[medium] * (1 - (10* deltaTime)/(firstDaughterRECO.PrimaryVertexToHit/speedOfLight_ + firstDaughterRECO.measuredHitTime[medium] + secondDaughterRECO.PrimaryVertexToHit/speedOfLight_ + secondDaughterRECO.measuredHitTime[medium]));
            Parent.flightLength[medium] = Parent.flightLength[medium] * (1 - (firstDeltaTime + secondDeltaTime));
        }
           
        if (Debug2) {cout << "New displacedVertexLength: " << Parent.flightLength[medium] << endl;}
        
        // Now modify the direction of the parent flight path between the two particle hits
        // if one deltaTime is +ve & the other -ve then slightly rotate the parentMC path around the primary vertex
        
        if ((firstDeltaTime * secondDeltaTime) < 0.0)     // Only if times are of the opposite sign
        {
            if(Debug2){cout << "Opposite sign delta times" << endl;}
            deltaOmega = deltaOmega/2;// Only a small change required
            // float deltaOmega = 0.1;   // Only a small change required
            // Algorithm to modify angle. delta time will be small wrt actual times of flight
            // float deltaOmega = firstDaughterRECO.omega[medium] * (1 - (deltaTime/(firstDaughterRECO.PrimaryVertexToHit/speedOfLight_ + firstDaughterRECO.measuredHitTime[medium] + secondDaughterRECO.PrimaryVertexToHit/speedOfLight_ + secondDaughterRECO.measuredHitTime[medium])));
            if(Debug2){cout << "Old firstDaughterRECO omega[" << 1 << "]: " << firstDaughterRECO.omega[medium] << ", Old secondDaughterRECO omega[" << 1 << "]: " << secondDaughterRECO.omega[medium] << ", deltaOmega: " << deltaOmega << endl;}
         
            if((firstDeltaTime > 0.0) && (secondDeltaTime < 0.0))
            {    
                // If constructed flight time > true flight time, then increase electron path length i.e path of faster particle
                // angle between ECAL hit to primary vertex to displaced vertex for daughter particles
                if (Debug2){cout << "here first" << endl;}
                firstDaughterRECO.omega[medium] = firstDaughterRECO.omega[medium] - deltaOmega;
                secondDaughterRECO.omega[medium] = secondDaughterRECO.omega[medium] + deltaOmega;
            }
            else if((firstDeltaTime < 0.0) && (secondDeltaTime > 0.0))
            {
                if (Debug2){cout << "here second" << endl;}  
                // If constructed flight time < true flight time, then decrease electron path length i.e path of faster particle
                firstDaughterRECO.omega[medium] = firstDaughterRECO.omega[medium] + deltaOmega;
                secondDaughterRECO.omega[medium] = secondDaughterRECO.omega[medium] - deltaOmega;
            }
            if(Debug2){cout << "New firstDaughterRECO omega[" << 1 << "]: " << firstDaughterRECO.omega[medium] << ", New secondDaughterRECO omega[" << 1 << "]: " << secondDaughterRECO.omega[medium] << endl;}
        }
        
        
        // Adjust path lengths for both electrons for the change in displaced vertex and angle omega
        firstDaughterRECO.flightLength[medium] = getOppositeLength (Parent.flightLength[medium], firstDaughterRECO.PrimaryVertexToHit, firstDaughterRECO.omega[medium]);
        secondDaughterRECO.flightLength[medium] = getOppositeLength (Parent.flightLength[medium], secondDaughterRECO.PrimaryVertexToHit, secondDaughterRECO.omega[medium]);
        if(Debug2){cout << "firstDaughterRECO.flightLength[medium]: " <<  firstDaughterRECO.flightLength[medium] << ",secondDaughterRECO.flightLength[medium]: " << secondDaughterRECO.flightLength[medium] << endl;}
        
        // Now recalculate the angle between hits
        angleBetweenPrimaryVertexDisplacedVertexHitOne = getAngle(firstDaughterRECO.flightLength[medium],Parent.flightLength[medium],firstDaughterRECO.PrimaryVertexToHit);
        angleBetweenPrimaryVertexDisplacedVertexHitTwo = getAngle(secondDaughterRECO.flightLength[medium],Parent.flightLength[medium],secondDaughterRECO.PrimaryVertexToHit);
        if(Debug2){cout << "angleBetweenPrimaryVertexDisplacedVertexHitOne: " << angleBetweenPrimaryVertexDisplacedVertexHitOne << endl;}
        if(Debug2){cout << "angleBetweenPrimaryVertexDisplacedVertexHitTwo: " << angleBetweenPrimaryVertexDisplacedVertexHitTwo << endl;}

        AngleBetweenHits = (2*pi) - angleBetweenPrimaryVertexDisplacedVertexHitOne - angleBetweenPrimaryVertexDisplacedVertexHitTwo;  

        if (Debug2) {cout << "AngleBetweenHits : " << AngleBetweenHits << endl;}

        // Recalculate reconstructed mass
        // Initial approximation assumed a LL mass of zero and velocity of c
        Parent.mass[medium] = sqrt(2.0 * firstDaughterRECO.energy * secondDaughterRECO.energy * (1 - cos(AngleBetweenHits)));
        Parent.omega[medium] = AngleBetweenHits;
        if (Debug2) {cout << "Revised reconstructed mass[medium]: " << Parent.mass[medium] << endl;}

        // Store values for Graph with asymm errors
        // bestMass[ErrorsCount] = Parent.mass[medium];
        // bestFlightLength[ErrorsCount] = Parent.flightLength[medium];

        ErrorsCount++;
        if (Debug2 && (count >= 100)) {cout << "Count: " << count << endl;}
    }

    if (Debug2){cout << "Parent.flightLength[medium]: " << Parent.flightLength[medium] << endl;}
    if (Parent.flightLength[medium] < 0.0)
    {
            nonPhysicalDecayLength++;
    }
    if (Parent.flightLength[medium] >= 0.0)  // Flight length < 0 is non-physical
    {
        if (Debug2){cout << "Parent.mass[medium]: " << Parent.mass[medium] << endl;}
        if (Debug2){cout << "Parent.flightLength[medium]: " << Parent.flightLength[medium] << endl;}

        if(Debug2){cout << "massCalcCounter: " << massCalcCounter << endl;}
        massCalcCounter++;
        if (AngleBetweenHits > minimumAngleBetweenElectrons) // Very small angles between electrons has large errors
        {
            hists_["t1"]->Fill(Parent.mass[medium]);
            if(Debug2){cout << "t1 OK" << endl;}
            hists_["t2"]->Fill(Parent.flightLength[medium]);
            if(Debug2){cout << "t2 OK" << endl;}   
            hists_["t5"]->Fill(Parent.mass[medium], Parent.flightLength[medium]);    // Reconstructed mass vs constructed LL flight length
            if(Debug2){cout << "t5 OK" << endl;}
            hists_["t9"]->Fill(Parent.mass[medium], Parent.flightLength[medium]);    // Profile plot of reconstructed mass vs constructed LL flight length
            if(Debug2){cout << "t9 OK" << endl;}
            //hists_["t11"]->Fill(Parent.mass[medium], Parent.mass[medium] - Parent.mass[medium]);    // Profile plot of True (MC) mass vs True mass - constructed mass
            // if(Debug2){cout << "t11 OK" << endl;}

            // Plot daughter's hit time vs constructed flight length
            if(Debug2){cout << "firstDaughterRECO.measuredHitTime[" << 1 << "]: " << firstDaughterRECO.measuredHitTime[medium] << ", parentMC.flightLength[" << 1 << "]: " << Parent.flightLength[medium] << endl;}
            if (firstDaughterRECO.measuredHitTime[medium] != dummy)  // Only plot measured hit times
            {
                hists_["t7"]->Fill(firstDaughterRECO.measuredHitTime[medium], Parent.flightLength[medium]);    // firstDaughterRECO.measuredHitTime vs parentMC constructed flight length
            }
            if (secondDaughterRECO.measuredHitTime[medium] != dummy)  // Only plot measured hit times
            {
                hists_["t7"]->Fill(secondDaughterRECO.measuredHitTime[medium], Parent.flightLength[medium]);    // secondDaughterRECO.measuredHitTime vs parentMC constructed flight length
            }
        }
    }
    if(Debug2) {cout    << "Constructed displaced vertex[medium]: " << Parent.flightLength[medium] << " cm" << endl; }
    return;
}

// Take two hits. Return the displaced vertex shared by the two hits
float 
RecoFirst::findDisplacedVertexLength(myParticle &firstDaughterRECO, int i, myParticle &secondDaughterRECO, int j)
    // Using polar coords for the ellipses. The Primary Vertex is a shared focus. Each hit is the other focus for each ellipse. 
    // The flight time for each hit is the twice major axis for each ellipse.
    // Equations for ellipse using polar coords from one focus:
    // r(theta) = a*(1-eta*eta)/(1 +- eta * cos(theta))
    // We will always use r(theta) = a*(1-eta*eta)/(1 + eta * cos(theta)) This is correct, and has been validated with an Excel plot
    // For each ellipse:
    // eta = f/a, 
    // with f = 1/2 * distance between foci (i.e. half the distance between PV & hit)
    // a = major axis, i.e. half of flight time distance ... flight time distance = f + a + (a-f) = 2a
    // As an initial value assume displaced vertex is angled half way between the two ellipses, then theta1 = alhpa/2
{
        float decayLength1 = 1.0; // Arbitrary initial value for decay length for first ellipse
        float decayLength2 = 10.0; // Arbitrary initial value for decay length for second ellipse
    // We want DecayLength1 = DecayLength2 at a shared displaced vertex

    float InitialAngleBetweenHits = getInitialAngleBetweenLeptons(firstDaughterRECO, secondDaughterRECO);

    // Reject cases where the two electrons are very close together
    if (fabs(InitialAngleBetweenHits) < minimumAngleBetweenElectrons)
    {
        decayLength1 = dummy;   // return a negative decay length
        firstDaughterRECO.omega[i] = dummy;
        secondDaughterRECO.omega[j] = dummy;
        return dummy;
    }
        
    // We want DecayLength1 = DecayLength2 at a shared displaced vertex
        float omega1 = InitialAngleBetweenHits/2.0; // Give omega1 & omega2 aribitrary starting values (in radians)
        float omega2 = InitialAngleBetweenHits/2.0; 
        int count = 0;
        float deltaOmega1 = omega1/2.0;

        // For first ellipse
        float a1 = firstDaughterRECO.timedDistanceToHit[i]/2.0;
        float f1 = firstDaughterRECO.PrimaryVertexToHit/2.0;
        float eta1 = f1/a1;
        if (Debug2) {cout << "Major axis 1: " << a1 << ", foci1: " << f1 << ", eta1: " << eta1 << endl;}

    // For second ellipse
        float a2 = secondDaughterRECO.timedDistanceToHit[j]/2.0;
        float f2 = secondDaughterRECO.PrimaryVertexToHit/2.0;
        float eta2 = f2/a2;
        if (Debug2) {cout << "Major axis 2: " << a2 << ", foci2: " << f2 << ", eta2: " << eta2 << endl;}

    // Loop until ratio of difference in decay lengths < decayLengthAccuracy i.e. 1%
        while ( (fabs((decayLength1-decayLength2)/(decayLength1 + decayLength2)) > decayLengthAccuracy) && (count < 100) && (deltaOmega1 > minDeltaOmega))
        {
                if (Debug2) {cout << "count: " << count << ", omega1: " << omega1 << ", omega2: " << omega2 << endl;}
                
        decayLength1 = a1*(1-eta1*eta1)/(1 + eta1*cos(omega1));
        if (Debug2){cout << "a1: " << a1 << ", eta1: " << eta1 << ", decayLength1 = a1*(1-eta1*eta1)/(1 + eta1*cos(omega1)): " << decayLength1 << endl;}
       
            if (Debug2) {cout << "omega1: " << omega1 << ", displaced vertex 1: " << decayLength1 << endl;}
        
        decayLength2 = a2*(1-eta2*eta2)/(1 + eta2*cos(omega2));
        if (Debug2){cout << "a2: " << a2 << ", eta2: " << eta2 <<", decayLength2 = a2*(1-eta2*eta2)/(1 + eta2*cos(omega2)):" << decayLength2 << endl;}
        
                if (Debug2) {cout << "omega2: " << omega2 << ", displaced vertex 2: " << decayLength2 << endl;}

                // compare decayLength1 & decayLength2
                if (decayLength1 > decayLength2)
                {
                omega1 = omega1 - deltaOmega1;  
        //omega1 = omega1 + deltaOmega1;
                    if (Debug2) {cout << "Length1 > Length2, so omega1: " << omega1 << ", deltaOmega1: " << deltaOmega1 << endl;}
                }
                if (decayLength1 < decayLength2)
                {
            omega1 = omega1 + deltaOmega1;
                        // omega1 = omega1 - deltaOmega1;
                    if (Debug2) {cout << "Length1 < Length2, so omega1: " << omega1 << ", deltaOmega1: " << deltaOmega1 << endl;}
                }
                omega2 = InitialAngleBetweenHits - omega1;
                deltaOmega1 = deltaOmega1/2.0;
                if (Debug2) {cout << "deltaOmega1: " << deltaOmega1 << ", omega2: " << omega2 << endl;}

                // else if decayLength1 = decayLength2, do nothing
                ++count;
        }
        if (Debug2) { cout << "Final calculated decay length is: " << decayLength1 << ", omega1: " << omega1 << ", omega2: " << omega2<< endl << endl;}
    firstDaughterRECO.omega[i] = omega1;
    secondDaughterRECO.omega[j] = omega2;
    firstDaughterRECO.flightLength[i] = getOppositeLength(decayLength1, firstDaughterRECO.PrimaryVertexToHit, omega1);
    secondDaughterRECO.flightLength[j] = getOppositeLength(decayLength1, secondDaughterRECO.PrimaryVertexToHit, omega2);
    if(Debug2){cout << "firstDaughterRECO.omega: " << omega1 << ", firstDaughterRECO.flightLength[" << i << "]: " << firstDaughterRECO.flightLength[i] << endl;}
    if(Debug2){cout << "secondDaughterRECO.omega: " << omega2 << ", secondDaughterRECO.flightLength[" << j << "]: " << secondDaughterRECO.flightLength[j] << endl;}
    
    if ((firstDaughterRECO.omega[i] == 0) || (secondDaughterRECO.omega[j] == 0) )    // algorithm never converged
    {
        noDisplacedVertexConvergence++;
        decayLength1 = dummy;   // return a negative decay length
        firstDaughterRECO.omega[i] = dummy;
        secondDaughterRECO.omega[j] = dummy;
    }
    hists_["t21"]->Fill(firstDaughterRECO.PrimaryVertexToHit);   // firstDaughterRECO PrimaryVertexToHit
    hists_["t22"]->Fill(secondDaughterRECO.PrimaryVertexToHit);   // secondDaughterRECO PrimaryVertexToHit
    
    // Simply take half the average distance for the two electrons    
    // hists_["t20"]->Fill((firstDaughterRECO.PrimaryVertexToHit + secondDaughterRECO.PrimaryVertexToHit)/4.0);   // Initial flight length estimate
    // return (firstDaughterRECO.PrimaryVertexToHit + secondDaughterRECO.PrimaryVertexToHit)/4.0;
    hists_["t20"]->Fill(decayLength1);   // Initial flight length estimate
    return decayLength1;
    
 }

// Get daughter particle energy, ECAL hit time, hit position and distance of hit from primary vertex
// The daughter particles are 'passed by reference' so their parameters are updated
template <class photonType> void 
RecoFirst::getDaughterParameters (myParticle &daughterParticle, photonType photon)
{
    if(Debug2){cout << endl << "Get new daughter" << endl;}
    // Reset values to ensure previous values are not retained
    daughterParticle.energy = dummy;
    daughterParticle.pt = dummy;
    daughterParticle.et = dummy;
    daughterParticle.eta = dummy;
    daughterParticle.phi = dummy;
    daughterParticle.measuredHitTime[medium] = dummy;
    daughterParticle.timedDistanceToHit[medium] = dummy;
    daughterParticle.PrimaryVertexToHit = dummy;

    if (photon->superCluster().isNull() )
        {
                return;
        }

    daughterParticle.energy = photon->energy();
    daughterParticle.pt = photon->pt();
    daughterParticle.et = photon->et();
    daughterParticle.eta = photon->eta();
    daughterParticle.phi = photon->phi();
    if(Debug2){cout << "Lepton energy: " << daughterParticle.energy << endl;}

    // Use the DetId of the photon to find the time of the ECAL hit
        // Then get the time of the crystal with the maximum hit energy

        EBRecHitCollection::const_iterator EBRecHitItr;
        EERecHitCollection::const_iterator EERecHitItr;

    EcalRecHit EcalHit;
    // reducedEcalRecHitsEB EcalHitEB; // AOD modification
    // reducedEcalRecHitsEE EcalHitEE;

        float MaxEcalHitEnergy = 0.0;   // Initialisation value
    
    // The photon's superCluster returns a collection of crystals' DetId and energy as a fraction of the total energy
        const vector< pair < DetId, float > > PhotonHitsAndFraction = photon->superCluster()->hitsAndFractions();
        // Loop over the crystals in the cluster to get a DetId for each crystal
        for (unsigned HitIndex = 0; HitIndex < PhotonHitsAndFraction.size(); ++HitIndex)
        {       
                const DetId PhotonHit = PhotonHitsAndFraction[HitIndex].first;  // 'first' here defeines the CMSSW parameter 'Hits'
                // Look for a hit in the Ecal barrel or endcap that matches the photon hit from the cluster
                EBRecHitItr = EBhits_->find(PhotonHit);

                if (EBRecHitItr != EBhits_->end()) // We have a hit in the barrel
                {
                        EcalHit = (*EBRecHitItr);
            // EcalHitEB = (*EBRecHitItr);   // AOD modification
                }
                else    // Look for a hit in the  endcap that matches the photon hit
                {       
                        EERecHitItr = EEhits_->find(PhotonHit);
                        if (EERecHitItr != EEhits_->end()) // We have a hit in the endcap
                        {
                                EcalHit = (*EERecHitItr);
                // EcalHitEE = (*EBRecHitItr);   // AOD modification

                        }
                }
                // If we have a hit in either the Ecal barrel or endcap
                // Get the hit time for the crystal with the largest energy
                if ((EBRecHitItr != EBhits_->end()) || (EERecHitItr != EEhits_->end()))
                {
            if ((EcalHit.energy() > MaxEcalHitEnergy) && (EcalHit.energy() > crystalMinEnergy)) // Get the hit time for the crystal with the greatest energy
                        {
                MaxEcalHitEnergy = EcalHit.energy();
                daughterParticle.measuredHitTime[medium] = EcalHit.time();
                        }
                }
        }

    // Now populate the measuredHitTime error values
    // [low] is minimim, [medium] is measured, [high] is maximum
    // In barrel if |eta| < 1.44
    // If hit in barrel error +- 0.27ns. If in endcap error +- 0.18ns
    if (fabs(daughterParticle.eta) <= barrelLimit)
    {
        daughterParticle.measuredHitTime[low] = daughterParticle.measuredHitTime[medium] - barrelTimingError;
        daughterParticle.measuredHitTime[high] = daughterParticle.measuredHitTime[medium] + barrelTimingError;
    }
    else // Must be in endcap
    {
        daughterParticle.measuredHitTime[low] = daughterParticle.measuredHitTime[medium] - endcapTimingError;
        daughterParticle.measuredHitTime[high] = daughterParticle.measuredHitTime[medium] + endcapTimingError;
    }

    if(Debug2)
    {
        if (daughterParticle.measuredHitTime[medium] == dummy)  // Show particles with no times set
        {
            cout << "particle eta: "                << daughterParticle.eta << endl
                 << "daughterParticle.measuredHitTime[low]: "   << daughterParticle.measuredHitTime[low] << endl
                 << "daughterParticle.measuredHitTime[medium]: "   << daughterParticle.measuredHitTime[medium] << endl
                 << "daughterParticle.measuredHitTime[high]: "   << daughterParticle.measuredHitTime[high] << endl;
        }
    }

    int a = 0;
    if (measured) {a=1;}
    else if (minimum) {a=0;}
    else if (maximum) {a=2;}

        // Distance from primary vertex to photon hit
    daughterParticle.endPoint = photon->caloPosition(); // Where the photon hits the ECAL
        math::XYZPoint PrimaryVertexToHit;
        PrimaryVertexToHit = daughterParticle.endPoint - PrimaryVertex_;
        float x = PrimaryVertexToHit.x();
        float y = PrimaryVertexToHit.y();
        float z = PrimaryVertexToHit.z();
        daughterParticle.PrimaryVertexToHit = sqrt(x*x + y*y + z*z);
        // daughterParticle.timedDistanceToHit[medium] = daughterParticle.PrimaryVertexToHit + (daughterParticle.measuredHitTime[medium] * speedOfLight_);
        // CMSSW time = 0.0 is for flight from the IP, NOT from primary vertex. 
    // So must account for IP to primary vertex displacement
    // Find distance (from IP to hit) + distance for hit time
    // Note: IP is at (0,0,0)
    // TimedDistance is time * speed of light (30 cm/ns)
    x = daughterParticle.endPoint.x();
        y = daughterParticle.endPoint.y();
        z = daughterParticle.endPoint.z();
    float IpToHitPlusTimeDistance [3];
    for (int i=0; i<3; i++)
    {
        IpToHitPlusTimeDistance [i]= sqrt(x*x + y*y + z*z) + (daughterParticle.measuredHitTime[i] * speedOfLight_);
    }
    // Distance from IP to primary vertex
    x = PrimaryVertex_.x();
        y = PrimaryVertex_.y();
        z = PrimaryVertex_.z();
        float IpToPrimaryVertex = sqrt(x*x + y*y + z*z);
    // Angle between hit->IP->primary vertex
    TVector3 x1;
    TVector3 x0;
    x0.SetXYZ(daughterParticle.endPoint.x(), daughterParticle.endPoint.y(), daughterParticle.endPoint.z() );
    x1.SetXYZ(PrimaryVertex_.x(), PrimaryVertex_.y(), PrimaryVertex_.z());
    float hitToIpToPrimaryVertexAngle = x0.Angle(x1);
    // Using the cosine rule and knowing two sides of a triangle and the included angle
    // c = sqrt(a*a + b*b - 2*a*b*cos(theta))
    for (int i=0; i<3; i++)
    {
        float a = IpToHitPlusTimeDistance[i];
        float b = IpToPrimaryVertex;
        // float c = sqrt(a*a + b*b - 2*a*b*cos(hitToIpToPrimaryVertexAngle));
        daughterParticle.timedDistanceToHit[i] = sqrt(a*a + b*b - 2*a*b*cos(hitToIpToPrimaryVertexAngle));
        if (Debug2){cout << "timedDistanceToHit[" << i << "]: " << daughterParticle.timedDistanceToHit[i] << ", PrimaryVertexToHit: " << daughterParticle.PrimaryVertexToHit << endl;}
        }
        return;
}

// Get daughter particle statistics: energy, ECAL hit time, hit position and distance of hit from primary vertex
// The daughter particles are 'passed by reference' so their parameters are updated
template <class photonType> void 
RecoFirst::getDaughterStats (myParticle &daughterParticle, photonType photon)
{
    numberOfDaughters++;

    if (photon->superCluster().isNull() )
        {
        nullSuperClusters++;
                return;
        }

    // Use the DetId of the photon to find the time of the ECAL hit
        // Then get the time of the crystal with the maximum hit energy

        EBRecHitCollection::const_iterator EBRecHitItr;
        EERecHitCollection::const_iterator EERecHitItr;
        EcalRecHit EcalHit;
        float MaxEcalHitEnergy = 0.0;   // Initialisation value
    
    // The photon's superCluster returns a collection of crystals' DetId and energy as a fraction of the total energy
        const vector< pair < DetId, float > > PhotonHitsAndFraction = photon->superCluster()->hitsAndFractions();
        // Loop over the crystals in the cluster to get a DetId for each crystal
        for (unsigned HitIndex = 0; HitIndex < PhotonHitsAndFraction.size(); ++HitIndex)
        {       
                const DetId PhotonHit = PhotonHitsAndFraction[HitIndex].first;  // 'first' here defines the CMSSW parameter 'Hits'
                // Look for a hit in the Ecal barrel or endcap that matches the photon hit from the cluster
                EBRecHitItr = EBhits_->find(PhotonHit);

                if (EBRecHitItr != EBhits_->end()) // We have a hit in the barrel
                {
                        EcalHit = (*EBRecHitItr);
                }
                else    // Look for a hit in the  endcap that matches the photon hit
                {       
                        EERecHitItr = EEhits_->find(PhotonHit);
                        if (EERecHitItr != EEhits_->end()) // We have a hit in the endcap
                        {
                                EcalHit = (*EERecHitItr);
                        }
                }
                // If we have a hit in either the Ecal barrel or endcap
                // Get the hit time for the crystal with the largest energy
                if ( (EBRecHitItr != EBhits_->end()) || (EERecHitItr != EEhits_->end()) )
                {
                        if (EcalHit.energy() > MaxEcalHitEnergy)        // Get the hit time for the crystal with the greatest energy
                        {
                                MaxEcalHitEnergy = EcalHit.energy();
                if (EcalHit.energy() > crystalMinEnergy)    // crystal energy below this minimum has very bad timing resolution
                {
                    daughterParticle.measuredHitTime[medium] = EcalHit.time();
                }
                else
                {
                    daughterParticle.measuredHitTime[medium] = dummy;    // Make this a negative value so don't retain the previous value
                }
                        }
                }
        }
    if(Debug2){cout << "daughterParticle.measuredHitTime[medium]: " << daughterParticle.measuredHitTime[medium] << endl;}
    if (daughterParticle.measuredHitTime[medium] == dummy)
    {
        lowCrystalHitEnergy++;
    }
    if ((daughterParticle.measuredHitTime[medium] < 0.0) && (daughterParticle.measuredHitTime[medium] != dummy))
    {
        negativeHitTime++;
    }

    if (daughterParticle.measuredHitTime[medium] != dummy)   // Don't plot dummy values when crystal energy below minimum
    {
        hists_["t8"]->Fill(daughterParticle.measuredHitTime[medium]);
    }
        return;
}


float 
RecoFirst::calcDeltaPhi(float phi1, float phi2)
{
    // Calculate the difference between phi1 & phi2
    // Phi ranges from 0 to +pi and 0 to -pi
        // deltaPhi lies between 0 to pi
    float deltaPhi;
        if (Debug2) {cout << "phi1: " <<  phi1 << ", phi2: " << phi2 << endl;}
    if (phi1 < 0) {phi1 = (2.0 * pi) + phi1;}
        if (phi2 < 0) {phi2 = (2.0 * pi) + phi2;}
        deltaPhi = fabs(phi1 - phi2);
        if (deltaPhi > pi) {deltaPhi = (2.0 * pi) - deltaPhi;}
        if (Debug2) { cout << "phi1: " << phi1 << ", phi2: " << phi2 << ", deltaPhi: " << deltaPhi << endl;}
        return deltaPhi;
}

float 
RecoFirst::calcDeltaEta(float eta1, float eta2)
{
    // Calculates the difference between eta1 & eta2
    // eta ranges from -10 to + 10 (or anyway - a large value and we don't care beyond abs(2.5))
    float deltaEta;
    if ( ((eta1 <= 0.0) && (eta2 <= 0.0)) || ((eta1 >= 0.0) && (eta2 >= 0.0)) )
    {
        return deltaEta = fabs(eta1 - eta2);
    }
    
    else if ( ((eta1 >= 0.0) && (eta2 <= 0.0)) || ((eta1 <= 0.0) && (eta2 >= 0.0)) )
    {
        return deltaEta = fabs(eta1) + fabs(eta2);
    }
    else
        return 0;
}

float 
RecoFirst::getOppositeLength (float a, float b, float theta)
{
        // Using the cosine rule we can find the length of one side of a triangle
        // c = sqrt(a*a + b*b - 2*a*b*cos(theta))
        // where
        // a = IP To ECAL Hit Length
        // b = IP to Primary Vertex length
    // c = Primary Vertex To Hit Lenght
    // theta = angle between ECAL hit, IP and Primary vertex
    if(Debug2){cout << "a: " << a << ", b: " << b << ", theta: " << theta << endl;}
    if(Debug2){cout << "(a*a + b*b - 2*a*b*cos(theta): " << (a*a + b*b - 2*a*b*cos(theta)) << endl;}
    float parameter = (a*a + b*b - 2*a*b*cos(theta)) ;
    // Check for sqrt of a negative value
    if (parameter < 0)
    {
        if(Debug2)
        {
            cout << "Can't take sqrt of " << parameter << ")" << endl;
            cout << "a: " << a << ", b: " << b << ", theta: " << theta << endl;
        }
        return 0.0;
    }
    else
    {
            float c = sqrt(a*a + b*b - 2*a*b*cos(theta));
        if(Debug2){cout << "c: " << c << endl;}
            return c;
    }
}

float 
RecoFirst::findIncludedAngleBetweenHits(float displacedVertex, myParticle firstDaughterRECO, myParticle secondDaughterRECO)
{
    // For both electrons get the internal angle between the DisplacedVertexLength_ and the (TimedDistanceToHit - DisplacedVertexLength_)
    // We now know 3 sides of a triangle - the DistanceFromPrimaryVertexToHit, the DisplacedVertexLength_, and the (TimedDistanceToHit - DisplacedVertexLength_)
    // Using the cosine rule we can find the internal angles.
    if (Debug2) {cout << "firsttimedDistanceToHit[medium] : " << firstDaughterRECO.timedDistanceToHit[medium] << endl;}
    if (Debug2) {cout << "secondtimedDistanceToHit[medium] : " << secondDaughterRECO.timedDistanceToHit[medium] << endl;}

    float AngleBetweenHits;

    float firstAlpha = getAngle(displacedVertex, firstDaughterRECO.flightLength[medium], firstDaughterRECO.PrimaryVertexToHit);
    if (Debug2) {cout << "firstAlpha: " << firstAlpha << ", displacedVertex: " << displacedVertex << ", firstDaughterRECO.flightLength[medium] : " << firstDaughterRECO.flightLength[medium] << ", firstDaughterRECO.PrimaryVertexToHit: " << firstDaughterRECO.PrimaryVertexToHit << endl;}
  
    float secondAlpha = getAngle(displacedVertex, secondDaughterRECO.flightLength[medium], secondDaughterRECO.PrimaryVertexToHit);
    if (Debug2) {cout << "secondAlpha: " << secondAlpha << ", displacedVertex: " << displacedVertex << ", secondDaughterRECO.flightLength[medium] : " << secondDaughterRECO.flightLength[medium] << ", secondDaughterRECO.PrimaryVertexToHit: " << secondDaughterRECO.PrimaryVertexToHit << endl;}

    if ((firstAlpha + secondAlpha) > pi)
    {
        // The total included angle between the two decay paths = (2 * pi - Alpha0 - Alpha1)
        AngleBetweenHits = (2.0 * pi) - firstAlpha - secondAlpha;
        if (Debug2) {cout << "firstAlpha: " << firstAlpha << ", secondAlpha: " << secondAlpha << ", AngleBetweenHits: " << AngleBetweenHits << endl;}
    }
    else
    {
        AngleBetweenHits = firstAlpha + secondAlpha;
    }
    return AngleBetweenHits;
}

float 
RecoFirst::getAngle (float a, float b, float c)
{
        // Using the cosine rule we can find the internal angle, alpha.
        // Alpha = arc-cos(a*a + b*b - c*c)/2*a*b
        // where
        // a = DisplacedVertexLength
        // c = DistanceFromPrimaryVertexToHit
        // b = timedDistanceToHit[medium] - DisplacedVertexLength
    if(Debug2){cout << "a: " << a << ", b: " << b << ", c: " << c << endl;}
    if(Debug2){cout << "(a*a + b*b - c*c)/(2*a*b): " << (a*a + b*b - c*c)/(2*a*b) << endl;}
    float parameter = (a*a + b*b - c*c)/(2*a*b);
    // Check for bad parameters for acos() i.e. > 1 or < -1
    float Alpha = 0.0;
    if (parameter > 1)
    {
        Alpha = acos(1.0);
        if(Debug2)
        {
            cout << "Bad angle: Acos(" << parameter << ")" << endl;
            cout << "a: " << a << ", b: " << b << ", c: " << c << endl;
            cout << "Angle: " << Alpha << endl;
        }
    }
    if (parameter < -1)
    {
        Alpha = acos(-1.0);
        if(Debug2)
        {
            cout << "Bad angle: Acos(" << parameter << ")" << endl;
            cout << "a: " << a << ", b: " << b << ", c: " << c << endl;
            cout << "Angle: " << Alpha << endl;
        }
    }
    else
    {
            Alpha = acos(parameter);
        if(Debug2){cout << "Angle: " << Alpha << endl;}
    }
    return Alpha;
}

// The angle between the firstDaughterRECO particle's ECAL hit, the IP and the secondDaughterRECO particle's ECAL hit
float 
RecoFirst::getInitialAngleBetweenLeptons(myParticle &firstDaughterRECO, myParticle &secondDaughterRECO)
{       TVector3 x1;
        TVector3 x0;  
    if(Debug2){cout << "firstDaughterRECO.endPoint.x(): " << firstDaughterRECO.endPoint.x() << endl;}               
    if(Debug2){cout << "secondDaughterRECO.endPoint.x(): " << secondDaughterRECO.endPoint.x() << endl;}
    if(Debug2){cout << "firstDaughterRECO.endPoint.y(): " << firstDaughterRECO.endPoint.y() << endl;}               
    if(Debug2){cout << "secondDaughterRECO.endPoint.y(): " << secondDaughterRECO.endPoint.y() << endl;}
    if(Debug2){cout << "firstDaughterRECO.endPoint.z(): " << firstDaughterRECO.endPoint.z() << endl;}               
    if(Debug2){cout << "secondDaughterRECO.endPoint.z(): " << secondDaughterRECO.endPoint.z() << endl;}

        x0.SetXYZ(firstDaughterRECO.endPoint.x(), firstDaughterRECO.endPoint.y(), firstDaughterRECO.endPoint.z() );
        x1.SetXYZ(secondDaughterRECO.endPoint.x(), secondDaughterRECO.endPoint.y(), secondDaughterRECO.endPoint.z() );
        if (Debug2) {cout << "Initial angle between hits: " << x0.Angle(x1) << endl;}
        return x0.Angle(x1);
}

// ------------ method called once each job just before starting event loop  ------------
void 
RecoFirst::beginJob()
{
}

void 
RecoFirst::book(Service<TFileService>& fs)
{
    if (Zee)
    {
            hists_["t1"] = fs->make<TH1F> ("Z->ee Ecal Timing Reconstructed Mass", "Reconstructed Mass (GeV)", 100, dummy, 200);
            hists_["t1"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");
    
    hists_["t2"] = fs->make<TH1F> ("Z->ee Calculated Flight Length (cm)", "Calculated Flight Length (cm)", 100, -5, 200);
        hists_["t2"]->GetXaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t3"] = fs->make<TH1F> ("MC (true) Flight Length (cm)", "MC (true) Flight Length (cm)", 100, dummy, 200);
        hists_["t3"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");

    hists_["t4"] = fs->make<TH2F> ("Z->ee MC (true) vs Calculated Flight Length (cm)", "MC (true) vs Calculated Flight Length (cm)", 100, 0, 200, 100, 0, 200);
        hists_["t4"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");
        hists_["t4"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t5"] = fs->make<TH2F> ("Z->ee Reconstructed Mass Vs Calculated Flight Length", "Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 100, 0, 2000, 100, 0, 200);
        hists_["t5"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");
        hists_["t5"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t6"] = fs->make<TH1F> ("Z->ee MC (true) Mass (GeV)", "MC (true) Mass (GeV)", 100, 0, 2000);
        hists_["t6"]->GetXaxis()->SetTitle("True MC Mass (GeV)");

    hists_["t7"] = fs->make<TH2F> ("Z->ee Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", "Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", 100, -0.5, 2, 100, 0, 20);
        hists_["t7"]->GetXaxis()->SetTitle("Decay Daughter's Hit Time (ns)");
        hists_["t7"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");

    hists_["t8"] = fs->make<TH1F> ("Z->ee Measured ECAL Hit Time", "Measured ECAL Hit Time", 100, -4, 3);
        hists_["t8"]->GetXaxis()->SetTitle("Measured ECAL Hit Time (ns)");
    
    hists_["t9"] = fs->make<TProfile> ("Z->ee Profile Reconstructed Mass Vs Calculated Flight Length", "Profile Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 200, 50, 1200, 0, 20);
 
    TGraphAsymmErrors_["t10"] = fs->make<TGraphAsymmErrors> ();

    hists_["t11"] = fs->make<TProfile> ("Z->ee Profile Reconstructed Mass Vs True Mass", "Profile Reconstructed Mass (GeV) vs MC Mass - Reconstructed Mass (GeV)", 200, 0, 1200, -50, 50);
    
    hists_["t12"] = fs->make<TH2F> ("Z->ee True MC Mass Vs Reconstructed Mass", "True MC Mass (GeV) vs Reconstructed Mass (GeV)", 100, dummy, 200, 100, dummy, 1000);
        hists_["t12"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t12"]->GetYaxis()->SetTitle("Reconstructed Mass (GeV)");
    
    hists_["t13"] = fs->make<TH2F> ("Z->ee Reconstructed Flight Length Residual vs Reconstructed Mass Residual", "Reconstructed Flight Length Residual vs Reconstructed Mass Residual", 100, dummy, 20, 100, -300, 300);
        hists_["t13"]->GetXaxis()->SetTitle("Reconstructed Flight Length - True MC Flight Length (cm)");
        hists_["t13"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");
    
    hists_["t14"] = fs->make<TH1F> ("Z->ee Valid Electron Pair Calculated Flight Length", "Valid Electron Pair Calculated Flight Length", 100, 0, 200);
        hists_["t14"]->GetXaxis()->SetTitle("Valid Electron Pair Calculated Flight Length (cm)");

    hists_["t15"] = fs->make<TH1F> ("Z->ee Valid Electrons Measured ECAL Hit Time", "Valid Electrons Measured ECAL Hit Time", 100, -0.5, 2);
        hists_["t15"]->GetXaxis()->SetTitle("Valid Electron Measured ECAL Hit Time (ns)");

    hists_["t16"] = fs->make<TH2F> ("Z->ee True MC Mass vs Reconstructed Mass Residual", "True MC Mass vs Reconstructed Mass Residual", 100, 0, 2000, 100, -1000, 1000);
        hists_["t16"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t16"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");
    }
    
    if (ZPrimeee)
    {
            hists_["t1"] = fs->make<TH1F> ("Z'->ee Ecal Timing Reconstructed Mass", "Z'->ee Reconstructed Mass (GeV)", 100, dummy, 1400);
            hists_["t1"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");

    hists_["t2"] = fs->make<TH1F> ("Z'->ee Calculated Flight Length (cm)", "Z'->ee Calculated Flight Length (cm)", 100, -5, 20);
        hists_["t2"]->GetXaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t3"] = fs->make<TH1F> ("Z'->ee MC (true) Flight Length (cm)", "Z'->ee MC (true) Flight Length (cm)", 100, dummy, 20);
        hists_["t3"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");

    hists_["t4"] = fs->make<TH2F> ("Z'->ee MC (true) vs Calculated Flight Length (cm)", "Z'->ee MC (true) vs Calculated Flight Length (cm)", 100, 0, 20, 100, 0, 20);
        hists_["t4"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");
        hists_["t4"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t5"] = fs->make<TH2F> ("Z'->ee Reconstructed Mass Vs Calculated Flight Length", "Z'->ee Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 100, 0, 1400, 100, 0, 20);
        hists_["t5"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");
        hists_["t5"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t6"] = fs->make<TH1F> ("Z'->ee MC (true) Mass (GeV)", "Z'->ee MC (true) Mass (GeV)", 100, 0, 1400);
        hists_["t6"]->GetXaxis()->SetTitle("True MC Mass (GeV)");

    hists_["t7"] = fs->make<TH2F> ("Z'->ee Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", "Z'->ee Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", 100, -0.5, 2, 100, 0, 20);
        hists_["t7"]->GetXaxis()->SetTitle("Decay Daughter's Hit Time (ns)");
        hists_["t7"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");

    hists_["t8"] = fs->make<TH1F> ("Z'->ee Measured ECAL Hit Time", "Z'->ee Measured ECAL Hit Time", 100, -4, 3);
        hists_["t8"]->GetXaxis()->SetTitle("Measured ECAL Hit Time (ns)");
    
    hists_["t9"] = fs->make<TProfile> ("Z'->ee Profile Reconstructed Mass Vs Calculated Flight Length", "Z'->ee Profile Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 200, 50, 1400, 0, 20);
 
    TGraphAsymmErrors_["t10"] = fs->make<TGraphAsymmErrors> ();

    hists_["t11"] = fs->make<TProfile> ("Z'->ee Profile Reconstructed Mass Vs True Mass", "Z'->ee Profile Reconstructed Mass (GeV) vs MC Mass - Reconstructed Mass (GeV)", 200, 0, 1400, -50, 50);
    
    hists_["t12"] = fs->make<TH2F> ("Z'->ee True MC Mass Vs Reconstructed Mass", "Z'->ee True MC Mass (GeV) vs Reconstructed Mass (GeV)", 100, dummy, 1400, 100, dummy, 1400);
        hists_["t12"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t12"]->GetYaxis()->SetTitle("Reconstructed Mass (GeV)");
    
    hists_["t13"] = fs->make<TH2F> ("Z'->ee Reconstructed Flight Length Residual vs Reconstructed Mass Residual", "Z'->ee Reconstructed Flight Length Residual vs Reconstructed Mass Residual", 100, -5, 20, 100, -300, 300);
        hists_["t13"]->GetXaxis()->SetTitle("Reconstructed Flight Length - True MC Flight Length (cm)");
        hists_["t13"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");
    
    hists_["t14"] = fs->make<TH1F> ("Z'->ee Valid Electron Pair Calculated Flight Length", "Z'->ee Valid Electron Pair Calculated Flight Length", 100, 0, 20);
        hists_["t14"]->GetXaxis()->SetTitle("Valid Electron Pair Calculated Flight Length (cm)");

    hists_["t15"] = fs->make<TH1F> ("Z'->ee Valid Electrons Measured ECAL Hit Time", "Z'->ee Valid Electrons Measured ECAL Hit Time", 100, -0.5, 2);
        hists_["t15"]->GetXaxis()->SetTitle("Valid Electron Measured ECAL Hit Time (ns)");

    hists_["t16"] = fs->make<TH2F> ("Z'->ee True MC Mass vs Reconstructed Mass Residual", "Z'->ee True MC Mass vs Reconstructed Mass Residual", 100, 0, 1400, 100, -400, 400);
        hists_["t16"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t16"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");
    }

    if (LLee)
    {
            hists_["t1"] = fs->make<TH1F> ("LL->ee Ecal Timing Reconstructed Mass", "LL->ee Reconstructed Mass (GeV)", 100, dummy, 1000);
            hists_["t1"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");

     hists_["t2"] = fs->make<TH1F> ("LL->ee Calculated Flight Length (cm)", "LL->ee Calculated Flight Length (cm)", 100, -5, 200);
        hists_["t2"]->GetXaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t3"] = fs->make<TH1F> ("LL->ee MC (true) Flight Length (cm)", "LL->ee MC (true) Flight Length (cm)", 100, dummy, 200);
        hists_["t3"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");

    hists_["t4"] = fs->make<TH2F> ("LL->ee MC (true) vs Calculated Flight Length (cm)", "LL->ee MC (true) vs Calculated Flight Length (cm)", 100, 0, 200, 100, 0, 200);
        hists_["t4"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");
        hists_["t4"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t5"] = fs->make<TH2F> ("LL->ee Reconstructed Mass Vs Calculated Flight Length", "LL->ee Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 100, 0, 200, 100, 0, 200);
        hists_["t5"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");
        hists_["t5"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");
    
    hists_["t6"] = fs->make<TH1F> ("LL->ee MC (true) Mass (GeV)", "LL->ee MC (true) Mass (GeV)", 100, 0, 200);
        hists_["t6"]->GetXaxis()->SetTitle("True MC Mass (GeV)");

    hists_["t7"] = fs->make<TH2F> ("LL->ee Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", "LL->ee Decay Daughter's Hit Time (ns) vs Calculated Flight Length (cm)", 100, -0.5, 2, 100, 0, 200);
        hists_["t7"]->GetXaxis()->SetTitle("Decay Daughter's Hit Time (ns)");
        hists_["t7"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");

    hists_["t8"] = fs->make<TH1F> ("LL->ee Measured ECAL Hit Time", "LL->ee Measured ECAL Hit Time", 100, -4, 3);
        hists_["t8"]->GetXaxis()->SetTitle("Measured ECAL Hit Time (ns)");
    
    hists_["t9"] = fs->make<TProfile> ("LL->ee Profile Reconstructed Mass Vs Calculated Flight Length", "LL->ee Profile Reconstructed Mass (GeV) vs Calculated Flight Length (cm)", 200, 50, 1200, 0, 200);
 
    TGraphAsymmErrors_["t10"] = fs->make<TGraphAsymmErrors> ();

    hists_["t11"] = fs->make<TProfile> ("LL->ee Profile Reconstructed Mass Vs True Mass", "LL->ee Profile Reconstructed Mass (GeV) vs MC Mass - Reconstructed Mass (GeV)", 200, 0, 1200, -50, 50);
    
    hists_["t12"] = fs->make<TH2F> ("LL->ee True MC Mass Vs Reconstructed Mass", "LL->ee True MC Mass (GeV) vs Reconstructed Mass (GeV)", 100, 140, 160, 100, 100, 200);
        hists_["t12"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t12"]->GetYaxis()->SetTitle("Reconstructed Mass (GeV)");
    
    hists_["t13"] = fs->make<TH2F> ("LL->ee Reconstructed Flight Length Residual vs Reconstructed Mass Residual", "LL->ee Reconstructed Flight Length Residual vs Reconstructed Mass Residual", 100, -150, 200, 100, -300, 1000);
        hists_["t13"]->GetXaxis()->SetTitle("Reconstructed Flight Length - True MC Flight Length (cm)");
        hists_["t13"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");
    
    hists_["t14"] = fs->make<TH1F> ("LL->ee Valid Electron Pair Calculated Flight Length", "LL->ee Valid Electron Pair Calculated Flight Length", 100, 0, 200);
        hists_["t14"]->GetXaxis()->SetTitle("Valid Electron Pair Calculated Flight Length (cm)");

    hists_["t15"] = fs->make<TH1F> ("LL->ee Valid Electrons Measured ECAL Hit Time", "LL->ee Valid Electrons Measured ECAL Hit Time", 100, -0.5, 2);
        hists_["t15"]->GetXaxis()->SetTitle("Valid Electron Measured ECAL Hit Time (ns)");

    hists_["t16"] = fs->make<TH2F> ("LL->ee True MC Mass vs Reconstructed Mass Residual", "LL->ee True MC Mass vs Reconstructed Mass Residual", 100, 140, 160, 100, -50, 50);
        hists_["t16"]->GetXaxis()->SetTitle("True MC Mass (GeV)");
        hists_["t16"]->GetYaxis()->SetTitle("Reconstructed Mass - True MC Mass (GeV)");

    hists_["t17"] = fs->make<TH1F> ("LL->ee Valid Reconstructed Mass (GeV)", "LL->ee Valid Reconstructed Mass (GeV)", 100, 0, 500);
        hists_["t17"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");

    hists_["t18"] = fs->make<TH1F> ("LL->ee Valid Reconstructed Flight Length (cm)", "LL->ee Valid Reconstructed Flight Length (cm)", 100, dummy, 100);
        hists_["t18"]->GetXaxis()->SetTitle("Reconstructed Flight Length (cm)");

    hists_["t19"] = fs->make<TH2F> ("LL->ee True MC Flight Length (cm) vs Reconstructed Flight Length (cm)", "LL->ee True MC Flight Length (cm) vs Reconstructed Flight Length (cm)", 100, 0, 200, 100, 0, 200);
        hists_["t19"]->GetXaxis()->SetTitle("True MC Flight Length (cm)");
        hists_["t19"]->GetYaxis()->SetTitle("Reconstructed Flight Length (cm)");
    
    hists_["t20"] = fs->make<TH1F> ("LL->ee Initial Decay Length Estimates (cm)", "LL->ee Initial Decay Length Estimates (cm)", 100, dummy, 500);
        hists_["t20"]->GetXaxis()->SetTitle("Initial Decay Length Estimates (cm)");

    hists_["t21"] = fs->make<TH1F> ("LL->ee firstDaughterRECO PrimaryVertexToHit (cm)", "LL->ee firstDaughterRECO PrimaryVertexToHit (cm)", 100, dummy, 500);
        hists_["t21"]->GetXaxis()->SetTitle("firstDaughterRECO PrimaryVertexToHit (cm)");

    hists_["t22"] = fs->make<TH1F> ("LL->ee secondDaughterRECO PrimaryVertexToHit (cm)", "LL->ee secondDaughterRECO PrimaryVertexToHit (cm)", 100, dummy, 500);
        hists_["t22"]->GetXaxis()->SetTitle("secondDaughterRECO PrimaryVertexToHit (cm)");

    hists_["t23"] = fs->make<TH2F> ("LL->ee Constructed Flight Length (cm) vs Included Angle (rad) vs ", "LL->ee Constructed Flight Length (cm) vs Included Angle (rad)", 100, 0, 20, 100, -pi, pi);
        hists_["t23"]->GetXaxis()->SetTitle("Constructed Flight Length (cm)");
        hists_["t23"]->GetYaxis()->SetTitle("Included Angle (rad)");

    hists_["t24"] = fs->make<TH2F> ("LL->ee True MC Flight Length (cm) vs Constructed Flight Length Residue (cm)", "LL->ee True MC Flight Length (cm) vs Constructed Flight Length Residue (cm)", 100, 0, 200, 100, -200, 20);
        hists_["t24"]->GetXaxis()->SetTitle("True MC Flight Length (cm)");
        hists_["t24"]->GetYaxis()->SetTitle("Constructed Flight Length Residue (cm)");

    hists_["t25"] = fs->make<TH2F> ("LL->ee Calculated Included Angle (rad) vs Valid Reconstructed Mass (GeV)", "Calculated Included Angle (rad) vs Valid Reconstructed Mass (GeV)", 100, -pi, pi, 100, 0, 200);
        hists_["t25"]->GetXaxis()->SetTitle("Calculated Included Angle (rad)");
        hists_["t25"]->GetYaxis()->SetTitle("Valid Reconstructed Mass (GeV)");

    hists_["t26"] = fs->make<TH2F> ("LL->ee ECAL Energy vs Valid Reconstructed Mass (GeV)", "ECAL Energy (GeV) vs Valid Reconstructed Mass (GeV)", 100, 0, 500, 100, 0, 200);
        hists_["t26"]->GetXaxis()->SetTitle("sqrt(Energy1 * Energy2) (GeV)");
        hists_["t26"]->GetYaxis()->SetTitle("Valid Reconstructed Mass (GeV)");

    hists_["t27"] = fs->make<TH2F> ("LL->ee ECAL Energy vs 1 - cos(Included Angle)", "ECAL Energy vs 1 - cos(Included Angle)", 100, 0, 500000, 100, 0, 2);
        hists_["t27"]->GetXaxis()->SetTitle("Energy1 * Energy2 (GeV)");
        hists_["t27"]->GetYaxis()->SetTitle("1 - cos(Included Angle)");

    hists_["t28"] = fs->make<TH2F> ("LL->ee MC (true) included angle vs MC (true) mass (GeV)", "MC (true) included angle vs MC (true) mass (GeV)", 100, 0, pi, 100, 140, 160);
        hists_["t28"]->GetXaxis()->SetTitle("MC (true) Included Angle (rad)");
        hists_["t28"]->GetYaxis()->SetTitle("MC (true) mass (GeV)");

    hists_["t29"] = fs->make<TH2F> ("LL->ee MC (true) included Angle vs Calculated included Angle", "MC (true) included angle vs Calculated included Angle", 100, 0, pi, 100, 0, pi);
        hists_["t29"]->GetXaxis()->SetTitle("MC (true) Included Angle (rad)");
        hists_["t29"]->GetYaxis()->SetTitle("Calculated included Angle (rad)");

    hists_["t30"] = fs->make<TH2F> ("LL->ee MC (true) Flight Length vs MC (true) included Angle", "MC (true) Flight Length vs MC (true) included Angle", 100, 0, 200, 100, 0, pi);
        hists_["t30"]->GetXaxis()->SetTitle("MC (true) Flight Length (cm)");
        hists_["t30"]->GetYaxis()->SetTitle("MC (true) included Angle (rad)");

    hists_["t31"] = fs->make<TH1F> ("LL->ee MC (true) included Angle (rad)", "LL->ee MC (true) included Angle (rad)", 100, 0, pi);
        hists_["t31"]->GetXaxis()->SetTitle("MC (true) included Angle (rad)");

    hists_["t32"] = fs->make<TH1F> ("LL->ee Calculated included Angle (rad)", "LL->ee Calculated included Angle (rad)", 100, 0, pi);
        hists_["t32"]->GetXaxis()->SetTitle("Calculated included Angle (rad)");

    hists_["t33"] = fs->make<TH2F> ("LL->ee MC (true) ECAL Hit Time vs Measured Hit Time (ns)", "MC (true) ECAL Hit Time vs Measured Hit Time (ns)", 100, 0, 10, 100, -2, 2);
        hists_["t33"]->GetXaxis()->SetTitle("MC (true) Calculated ECAL Hit Time (ns)");
        hists_["t33"]->GetYaxis()->SetTitle("Measured ECAL Hit Time (ns)");

    hists_["t34"] = fs->make<TH2F> ("LL->ee MC (true) Absolute ECAL Hit Time vs Absolute Measured Hit Time (ns)", "MC (true) Absolute ECAL Hit Time vs Absolute Measured Hit Time (ns)", 100, 0, 20, 100, 0, 20);
        hists_["t34"]->GetXaxis()->SetTitle("MC (true) Calculated Absolute ECAL Hit Time (ns)");
        hists_["t34"]->GetYaxis()->SetTitle("Measured Absolute ECAL Hit Time (ns)");

    hists_["t35"] = fs->make<TH1F> ("LL->ee MC LL Flight Time (ns)", "LL->ee MC LL Flight Time (ns)", 100, 0, 20);
        hists_["t32"]->GetXaxis()->SetTitle("MC LL Flight Time (ns)");

    hists_["t36"] = fs->make<TH1F> ("LL->ee MC Electron Flight Time (ns)", "LL->ee MC Electron Flight Time (ns)", 100, 0, 20);
        hists_["t36"]->GetXaxis()->SetTitle("MC Electron Flight Time (ns)");

    hists_["t37"] = fs->make<TH1F> ("LL->ee MC relatavistic Flight Time (ns)", "LL->ee MC relatavistic Flight Time (ns)", 100, 0, 20);
        hists_["t37"]->GetXaxis()->SetTitle("MC relatavistic Flight Time (ns)");

    }
}

// ------------ method called once each job just after ending the event loop  ------------
void 
RecoFirst::endJob() 
{
        Service<TFileService> fs;
        if (!fs)
        { throw Exception(errors::Configuration, "TFile Service is not registered in cfg file");
        }
    TGraphAsymmErrors_["t10"] = fs->make<TGraphAsymmErrors> ( n, bestMass, bestFlightLength, lowMass, lowFlightLength, highMass, highFlightLength);
    TGraphAsymmErrors_["t10"]->SetTitle("Reconstructed Mass (GeV) vs Calculated Flight Length (cm)");
    TGraphAsymmErrors_["t10"]->GetXaxis()->SetTitle("Reconstructed Mass (GeV)");
    TGraphAsymmErrors_["t10"]->GetYaxis()->SetTitle("Calculated Flight Length (cm)");

    if (Debug2)
    {
        cout << endl << endl << "Number of electrons: " << numberOfElectrons << endl;
        cout << "Number of photons: " << numberOfPhotons << endl;
        cout << "Number of number Of Daughters: " << numberOfDaughters << endl;
        cout << "Number of electron/photon pairs: " << numberOfElectronPhotonPairs << endl;
        cout << "Number of ECAL hits of energy < "<< crystalMinEnergy << " GeV: " << lowCrystalHitEnergy << endl;
        cout << "Number of negative hit times: " << negativeHitTime << endl;
        cout << "Number of Zs without two positive times: " << numberOfZsWithoutTwoPositiveTimes << endl;
        cout << "Number of Zs with a displaced vertex: " << numberOfDisplacedVertices << endl;
        cout << "Number of Zs without a displaced vertex convergence: " << noDisplacedVertexConvergence << endl;
        cout << "Number of Zs with a negative decay length: " << nonPhysicalDecayLength << endl;
        cout << "Number of Zs mass reconstructions: " << massReconstructions << endl;
        cout << "Number of validations: " << countValidations << endl;
        cout << "Number of validated reconstructions where both RECO & MC electrons are close : " << ValidatedReconstructions << endl;
        cout << "Number of close electrons with negative hit time : " << closeElectronsWithNegativeHitTime << endl;
        cout << "Number of negative flight lengths : " << negativeFlightLength << endl;
        cout << endl;

        cout << "Number of events with a primary vertex: " << numberOfPrimaryVertexEvents << endl;
        cout << "Number of validation events: " << ValidatedEvents << endl;
        cout << "Number of events with no particles: " << noParticles << endl;
        cout << "Number of Z particles: " << numberOfZParticles << endl;
        cout << "Number of events with Z above minimum mass of " << minimumZMass << ": " << numberOfZeventsAboveMinimumMass << endl;
        cout << "Number of events with minimum mass Z & 2 MC daughters: " << numberOfZeventsWithTwoMcDaughters << endl;
        cout << "Number of events with MC daughter not in final state: " << numberMcDaughterElectronsNotFinalState << endl;
        cout << "Number of Z decays with badRecoDeltaRMatch: " << badRecoDeltaRMatch << endl;
        cout << "Number of Z decays with goodRecoElectron: " << goodRecoElectron << endl;
        cout << "Number of Z decays with goodRecoPhoton: " << goodRecoPhoton << endl;
        cout << "Number of Z decays with noRecoMatch: " << noRecoMatch << endl;
        cout << "Number of leptons with nullSuperClusters: " << nullSuperClusters << endl;
        cout << "Number of leptons with zero energy: " << zeroEnergyLepton << endl;
        cout << "Number of Z events & found MC daughters: " << foundDaughters << endl;
        cout << "Number of Z decays with matched RECO daughters: " << zDecaysWithEnergeticDaughters << endl;
        cout << "Number of Zs with two RECO daughters: " << zWithTwoRecoDaughters << endl;
        cout << "Number of Zs entering veto: " << recoPairsGoingIntoVeto << endl;
        cout << "Number of Zs failing isolation: " << failedIsolation << endl;
        cout << "Number of Zs failing barrel-endcap gap: " << failedBarrelEndcapGap << endl;
        cout << "Number of Zs failing eta cut: " << failedEtaCut << endl;
        cout << "Number of Zs failing minimum Et cut: " << failedMinEt << endl;
        cout << "Number of Zs passing veto: " << recoPairsPassedVeto << endl;
        cout << "Number of reconstructions where electrons are too close: " << minimumAngleBetweenElectronsRejects << endl;
        cout << "Number of Zs with middle reconstructed mass: " << middleReconstructedMass << endl;
        cout << "Number of Zs reconstructed mass plotted: " << massCalcCounter << endl;
        cout << "Number of electron decays: " << ElectronDecays << endl;
        cout << "Number of photon decays: " << PhotonDecays << endl;
        cout << "Number of brem electrons: " << bremElectrons << endl;
        cout << "Number of muon decays: " << MuonDecays << endl;
        cout << "Number of tau decays: " << TauDecays << endl;
    }
}

// ------------ method called when starting to processes a run  ------------
void 
RecoFirst::beginRun(edm::Run const&, edm::EventSetup const&)
{
        Service<TFileService> fs;
        if (!fs)
        { throw Exception(errors::Configuration, "TFile Service is not registered in cfg file");
        }

        // book histogram
        book (fs);
}

// ------------ method called when ending the processing of a run  ------------
void 
RecoFirst::endRun(edm::Run const&, edm::EventSetup const&)
{
}

// ------------ method called when starting to processes a luminosity block  ------------
void 
RecoFirst::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{
}

// ------------ method called when ending the processing of a luminosity block  ------------
void 
RecoFirst::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void
RecoFirst::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(RecoFirst);
Revision:	1.1
Committed:	Mon Sep 19 15:48:29 2011 UTC (13 years, 7 months ago) by tonypoll
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	Electron Timing