ConversionProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    ConversionProducer
// Class:      ConversionProducer
//
//
// Original Authors:  Hongliang Liu
//         Created:  Thu Mar 13 17:40:48 CDT 2008
//
//


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


#include "RecoEgamma/EgammaPhotonProducers/interface/ConversionProducer.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"

#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

#include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
#include "DataFormats/EgammaTrackReco/interface/ConversionTrack.h"


#include "DataFormats/VertexReco/interface/Vertex.h"

#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"

#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"

#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"

#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
#include "DataFormats/Math/interface/deltaPhi.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionVertexFinder.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/TangentApproachInRPhi.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionHitChecker.h"


//Kinematic constraint vertex fitter
#include "RecoVertex/KinematicFitPrimitives/interface/ParticleMass.h"
#include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
#include <RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h>
#include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
#include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
#include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
#include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
#include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
#include "RecoVertex/KinematicFit/interface/ColinearityKinematicConstraint.h"



ConversionProducer::ConversionProducer(const edm::ParameterSet& iConfig):
  theVertexFinder_(nullptr)

{
  algoName_ = iConfig.getParameter<std::string>( "AlgorithmName" );

  src_ = 
    consumes<edm::View<reco::ConversionTrack> >(iConfig.getParameter<edm::InputTag>("src"));

  maxNumOfTrackInPU_ = iConfig.getParameter<int>("maxNumOfTrackInPU");
  maxTrackRho_ = iConfig.getParameter<double>("maxTrackRho");
  maxTrackZ_ = iConfig.getParameter<double>("maxTrackZ");

  allowTrackBC_ = iConfig.getParameter<bool>("AllowTrackBC");
  allowD0_ = iConfig.getParameter<bool>("AllowD0");
  allowDeltaPhi_ = iConfig.getParameter<bool>("AllowDeltaPhi");
  allowDeltaCot_ = iConfig.getParameter<bool>("AllowDeltaCot");
  allowMinApproach_ = iConfig.getParameter<bool>("AllowMinApproach");
  allowOppCharge_ = iConfig.getParameter<bool>("AllowOppCharge");

  allowVertex_ = iConfig.getParameter<bool>("AllowVertex");

  bypassPreselGsf_ = iConfig.getParameter<bool>("bypassPreselGsf");
  bypassPreselEcal_ = iConfig.getParameter<bool>("bypassPreselEcal");
  bypassPreselEcalEcal_ = iConfig.getParameter<bool>("bypassPreselEcalEcal");
  
  deltaEta_ = iConfig.getParameter<double>("deltaEta");
  
  halfWayEta_ = iConfig.getParameter<double>("HalfwayEta");//open angle to search track matches with BC

  d0Cut_ = iConfig.getParameter<double>("d0");
    
  usePvtx_ = iConfig.getParameter<bool>("UsePvtx");//if use primary vertices

  vertexProducer_   = 
    consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertexProducer"));
  

  //Track-cluster matching eta and phi cuts
  dEtaTkBC_ = iConfig.getParameter<double>("dEtaTrackBC");//TODO research on cut endcap/barrel
  dPhiTkBC_ = iConfig.getParameter<double>("dPhiTrackBC");
  
  bcBarrelCollection_ = 
    consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("bcBarrelCollection"));
  bcEndcapCollection_ = 
    consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("bcEndcapCollection"));
  
  scBarrelProducer_   = 
    consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("scBarrelProducer"));
  scEndcapProducer_   = 
    consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("scEndcapProducer"));
  
  energyBC_               = iConfig.getParameter<double>("EnergyBC");//BC energy threshold
  energyTotalBC_          = iConfig.getParameter<double>("EnergyTotalBC");//BC energy threshold
  minSCEt_                = iConfig.getParameter<double>("minSCEt");//super cluster energy threshold
  dEtacutForSCmatching_     = iConfig.getParameter<double>("dEtacutForSCmatching");// dEta between conversion momentum direction and SC position
  dPhicutForSCmatching_     = iConfig.getParameter<double>("dPhicutForSCmatching");// dPhi between conversion momentum direction and SC position

  
   

  //Track cuts on left right track: at least one leg reaches ECAL
  //Left track: must exist, must reach Ecal and match BC, so loose cut on Chi2 and tight on hits
  //Right track: not necessary to exist (if allowSingleLeg_), not necessary to reach ECAL or match BC, so tight cut on Chi2 and loose on hits
  maxChi2Left_ =  iConfig.getParameter<double>("MaxChi2Left");
  maxChi2Right_ =  iConfig.getParameter<double>("MaxChi2Right");
  minHitsLeft_ = iConfig.getParameter<int>("MinHitsLeft");
  minHitsRight_ = iConfig.getParameter<int>("MinHitsRight");

  //Track Open angle cut on delta cot(theta) and delta phi
  deltaCotTheta_ = iConfig.getParameter<double>("DeltaCotTheta");
  deltaPhi_ = iConfig.getParameter<double>("DeltaPhi");
  minApproachLow_ = iConfig.getParameter<double>("MinApproachLow");
  minApproachHigh_ = iConfig.getParameter<double>("MinApproachHigh");
  

  // if allow single track collection, by default False
  allowSingleLeg_ = iConfig.getParameter<bool>("AllowSingleLeg");
  rightBC_ = iConfig.getParameter<bool>("AllowRightBC");

  //track inner position dz cut, need RECO
  dzCut_ = iConfig.getParameter<double>("dz");
  //track analytical cross cut
  r_cut = iConfig.getParameter<double>("rCut");
  vtxChi2_ = iConfig.getParameter<double>("vtxChi2");


  theVertexFinder_ = new ConversionVertexFinder ( iConfig );

  thettbuilder_ = nullptr;

  //output
  ConvertedPhotonCollection_     = iConfig.getParameter<std::string>("convertedPhotonCollection");

  produces< reco::ConversionCollection >(ConvertedPhotonCollection_);

}


ConversionProducer::~ConversionProducer()
{

  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  delete theVertexFinder_;
}


// ------------ method called to produce the data  ------------
void
ConversionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;

  reco::ConversionCollection outputConvPhotonCollection;
  auto outputConvPhotonCollection_p = std::make_unique<reco::ConversionCollection>();

  //std::cout << " ConversionProducer::produce " << std::endl;
  //Read multiple track input collections

  edm::Handle<edm::View<reco::ConversionTrack> > trackCollectionHandle;
  iEvent.getByToken(src_,trackCollectionHandle);    

  //build map of ConversionTracks ordered in eta
  std::multimap<float, edm::Ptr<reco::ConversionTrack> > convTrackMap;
  for(auto const& t : trackCollectionHandle->ptrs()) convTrackMap.emplace(t->track()->eta(),t);

  edm::Handle<reco::VertexCollection> vertexHandle;
  reco::VertexCollection vertexCollection;
  if (usePvtx_){
    iEvent.getByToken(vertexProducer_, vertexHandle);
    if (!vertexHandle.isValid()) {
      edm::LogError("ConversionProducer") 
    << "Error! Can't get the product primary Vertex Collection "<< "\n";
      usePvtx_ = false;
    }
    if (usePvtx_)
      vertexCollection = *(vertexHandle.product());
  }
    
  edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
  thettbuilder_ = hTransientTrackBuilder.product();
    
  reco::Vertex the_pvtx;
  //because the priamry vertex is sorted by quality, the first one is the best
  if (!vertexCollection.empty())
    the_pvtx = *(vertexCollection.begin());
    
  if (trackCollectionHandle->size()> maxNumOfTrackInPU_){
    iEvent.put(std::move(outputConvPhotonCollection_p), ConvertedPhotonCollection_);
    return;
  }
    
    
  // build Super and Basic cluster geometry map to search in eta bounds for clusters
  std::multimap<double, reco::CaloClusterPtr> basicClusterPtrs;
  std::multimap<double, reco::CaloClusterPtr> superClusterPtrs;

  
  buildSuperAndBasicClusterGeoMap(iEvent,basicClusterPtrs,superClusterPtrs);
      
  buildCollection( iEvent, iSetup, convTrackMap,  superClusterPtrs, basicClusterPtrs, the_pvtx, outputConvPhotonCollection);//allow empty basicClusterPtrs
    
  outputConvPhotonCollection_p->assign(outputConvPhotonCollection.begin(), outputConvPhotonCollection.end());
  iEvent.put(std::move(outputConvPhotonCollection_p), ConvertedPhotonCollection_);
    
}


void ConversionProducer::buildSuperAndBasicClusterGeoMap(const edm::Event& iEvent,  
                             std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
                             std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs){

  // Get the Super Cluster collection in the Barrel
  edm::Handle<edm::View<reco::CaloCluster> > scBarrelHandle;
  iEvent.getByToken(scBarrelProducer_,scBarrelHandle);
  if (!scBarrelHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") 
      << "Error! Can't get the barrel superclusters!";
  }
    
  // Get the Super Cluster collection in the Endcap
  edm::Handle<edm::View<reco::CaloCluster> > scEndcapHandle;
  iEvent.getByToken(scEndcapProducer_,scEndcapHandle);
  if (!scEndcapHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") 
      << "Error! Can't get the endcap superclusters!";
  }
    
    
  edm::Handle<edm::View<reco::CaloCluster> > bcBarrelHandle;
  edm::Handle<edm::View<reco::CaloCluster> > bcEndcapHandle;//TODO check cluster type if BasicCluster or PFCluster

  iEvent.getByToken( bcBarrelCollection_, bcBarrelHandle);
  if (!bcBarrelHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") 
      << "Error! Can't get the barrel basic clusters!";
  }

  iEvent.getByToken( bcEndcapCollection_, bcEndcapHandle);
  if (! bcEndcapHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") 
      << "Error! Can't get the endcap basic clusters!";
  }

  if( bcBarrelHandle.isValid()  ) {    
    for(auto const& handle : {bcBarrelHandle, bcEndcapHandle} ) {
      for(auto const& bc : handle->ptrs()) {
        if(bc->energy() > energyBC_) basicClusterPtrs.emplace(bc->position().eta(), bc);
      }
    }
  }

  if( scBarrelHandle.isValid()  ) {    
    for(auto const& handle : {scBarrelHandle, scEndcapHandle} ) {
      for(auto const& sc : handle->ptrs()) {
        if(sc->energy() > minSCEt_) superClusterPtrs.emplace(sc->position().eta(), sc);
      }
    } 
  }

}


void ConversionProducer::buildCollection(edm::Event& iEvent, const edm::EventSetup& iSetup,
                                         const std::multimap<float, edm::Ptr<reco::ConversionTrack> >& allTracks,
                                         const std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs,
                                         const std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
                                         const reco::Vertex& the_pvtx,
                                         reco::ConversionCollection & outputConvPhotonCollection){
  
  
  edm::ESHandle<TrackerGeometry> trackerGeomHandle;
  edm::ESHandle<MagneticField> magFieldHandle;
  iSetup.get<TrackerDigiGeometryRecord>().get( trackerGeomHandle );
  iSetup.get<IdealMagneticFieldRecord>().get( magFieldHandle );
  
  const TrackerGeometry* trackerGeom = trackerGeomHandle.product();
  const MagneticField* magField = magFieldHandle.product();
  
//   std::vector<math::XYZPointF> trackImpactPosition;
//   trackImpactPosition.reserve(allTracks.size());//track impact position at ECAL
//   std::vector<bool> trackValidECAL;//Does this track reach ECAL basic cluster (reach ECAL && match with BC)
//   trackValidECAL.assign(allTracks.size(), false);
//   
//   std::vector<reco::CaloClusterPtr> trackMatchedBC;
//   reco::CaloClusterPtr empty_bc;
//   trackMatchedBC.assign(allTracks.size(), empty_bc);//TODO find a better way to avoid copy constructor
//   
//   std::vector<int> bcHandleId;//the associated BC handle id, -1 invalid, 0 barrel 1 endcap
//   bcHandleId.assign(allTracks.size(), -1);
  
  // not used    std::multimap<double, int> trackInnerEta;//Track innermost state Eta map to TrackRef index, to be used in track pair sorting
  
  std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF> trackImpactPosition;
  std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr> trackMatchedBC;
  
  ConversionHitChecker hitChecker;
  
  
  //2 propagate all tracks into ECAL, record its eta and phi
 
  for (std::multimap<float, edm::Ptr<reco::ConversionTrack> >::const_iterator tk_ref = allTracks.begin(); tk_ref != allTracks.end(); ++tk_ref ){
    const reco::Track* tk = tk_ref->second->trackRef().get()  ;
    
    
    //check impact position then match with BC
    math::XYZPointF ew;
    if ( getTrackImpactPosition(tk, trackerGeom, magField, ew) ){
      trackImpactPosition[tk_ref->second] = ew;
      
      reco::CaloClusterPtr closest_bc;//the closest matching BC to track
      
      if ( getMatchedBC(basicClusterPtrs, ew, closest_bc) ){
    trackMatchedBC[tk_ref->second] = closest_bc;
      }
    }    
  }
  
  
  
  //3. pair up tracks: 
  //TODO it is k-Closest pair of point problem
  //std::cout << " allTracks.size() " <<  allTracks.size() << std::endl;
  for(std::multimap<float, edm::Ptr<reco::ConversionTrack> >::const_iterator ll = allTracks.begin(); ll != allTracks.end();  ++ll ) {
    bool track1HighPurity=true;
    //std::cout << " Loop on allTracks " << std::endl;
    const  edm::RefToBase<reco::Track> & left = ll->second->trackRef();
    

    //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
    //(Note that the TrackRef and GsfTrackRef versions of the constructor are needed
    // to properly get refit tracks in the output vertex)
    reco::TransientTrack ttk_l;
    if (dynamic_cast<const reco::GsfTrack*>(left.get())) {
      ttk_l = thettbuilder_->build(left.castTo<reco::GsfTrackRef>());
    }
    else {
      ttk_l = thettbuilder_->build(left.castTo<reco::TrackRef>());
    }
    
    //      if ((allowTrackBC_ && !trackValidECAL[ll-allTracks.begin()]) )//this Left leg should have valid BC
    // continue;
    
    
    if (the_pvtx.isValid()){
      if (!(trackD0Cut(left, the_pvtx)))   track1HighPurity=false;
    } else {
      if (!(trackD0Cut(left)))  track1HighPurity=false;
    }
    
    std::vector<int> right_candidates;//store all right legs passed the cut (theta/approach and ref pair)
    std::vector<double> right_candidate_theta, right_candidate_approach;
    std::vector<std::pair<bool, reco::Vertex> > vertex_candidates;
    
    //inner loop only over tracks between eta and eta + deltaEta of the first track
    float etasearch = ll->first + deltaEta_;
    std::multimap<float, edm::Ptr<reco::ConversionTrack> >::const_iterator rr = ll;
    ++rr;    
    for (; rr != allTracks.lower_bound(etasearch); ++rr ) {
      bool track2HighPurity = true;
      bool highPurityPair = true;
      
      const  edm::RefToBase<reco::Track> & right = rr->second->trackRef();
      
      
      //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
      reco::TransientTrack ttk_r;
      if (dynamic_cast<const reco::GsfTrack*>(right.get())) {
        ttk_r = thettbuilder_->build(right.castTo<reco::GsfTrackRef>());
      }
      else {
        ttk_r = thettbuilder_->build(right.castTo<reco::TrackRef>());
      }
      //std::cout << " This track is " <<  right->algoName() << std::endl;
      
      
      //all vertexing preselection should go here
      
      //check for opposite charge
      if (  allowOppCharge_ && (left->charge()*right->charge() > 0) )  
        continue; //same sign, reject pair
          
      //if ( (allowTrackBC_ && !trackValidECAL[rr-allTracks.begin()] && rightBC_) )// if right track matches ECAL
      //  continue;
          
          
      double approachDist = -999.;
      //apply preselection to track pair, overriding preselection for gsf+X or ecalseeded+X pairs if so configured
      bool preselected = preselectTrackPair(ttk_l,ttk_r, approachDist);
      preselected = preselected || (bypassPreselGsf_ && (left->algo()==reco::TrackBase::gsf || right->algo()==reco::TrackBase::gsf));
      preselected = preselected || (bypassPreselEcal_ && (left->algo()==reco::TrackBase::outInEcalSeededConv || right->algo()==reco::TrackBase::outInEcalSeededConv || left->algo()==reco::TrackBase::inOutEcalSeededConv || right->algo()==reco::TrackBase::inOutEcalSeededConv));
      preselected = preselected || (bypassPreselEcalEcal_ && (left->algo()==reco::TrackBase::outInEcalSeededConv || left->algo()==reco::TrackBase::inOutEcalSeededConv) && (right->algo()==reco::TrackBase::outInEcalSeededConv || right->algo()==reco::TrackBase::inOutEcalSeededConv));
      
      if (!preselected) {
        continue;
      }
          
      //do the actual vertex fit
      reco::Vertex theConversionVertex;//by default it is invalid          
      bool goodVertex = checkVertex(ttk_l, ttk_r, magField, theConversionVertex);
          
      //bail as early as possible in case the fit didn't return a good vertex
      if (!goodVertex) {
        continue;
      }

          
          
      //track pair pass the quality cut
      if (   !( (trackQualityFilter(left, true) && trackQualityFilter(right, false))
                || (trackQualityFilter(left, false) && trackQualityFilter(right, true)) ) ) {
        highPurityPair=false;
      }

      if (the_pvtx.isValid()){
        if (!(trackD0Cut(right, the_pvtx))) track2HighPurity=false; 
      } else {
        if (!(trackD0Cut(right))) track2HighPurity=false; 
      }
        

      //if all cuts passed, go ahead to make conversion candidates
      std::vector<edm::RefToBase<reco::Track> > trackPairRef;
      trackPairRef.push_back(left);//left track
      trackPairRef.push_back(right);//right track

      std::vector<math::XYZVectorF> trackPin;
      std::vector<math::XYZVectorF> trackPout;
      std::vector<math::XYZPointF> trackInnPos;
      std::vector<uint8_t> nHitsBeforeVtx;
      std::vector<Measurement1DFloat> dlClosestHitToVtx;

      if (left->extra().isNonnull() && right->extra().isNonnull()){//only available on TrackExtra
        trackInnPos.push_back(  toFConverterP(left->innerPosition()));
        trackInnPos.push_back(  toFConverterP(right->innerPosition()));
        trackPin.push_back(toFConverterV(left->innerMomentum()));
        trackPin.push_back(toFConverterV(right->innerMomentum()));
        trackPout.push_back(toFConverterV(left->outerMomentum()));
    trackPout.push_back(toFConverterV(right->outerMomentum()));
        auto leftWrongHits = hitChecker.nHitsBeforeVtx(*left->extra(),theConversionVertex);
        auto rightWrongHits = hitChecker.nHitsBeforeVtx(*right->extra(),theConversionVertex);
        nHitsBeforeVtx.push_back(leftWrongHits.first);
        nHitsBeforeVtx.push_back(rightWrongHits.first);
        dlClosestHitToVtx.push_back(leftWrongHits.second);
        dlClosestHitToVtx.push_back(rightWrongHits.second);            
      }
          
      uint8_t nSharedHits = hitChecker.nSharedHits(*left.get(),*right.get());


      //if using kinematic fit, check with chi2 post cut
      if (theConversionVertex.isValid()){
        const float chi2Prob = ChiSquaredProbability(theConversionVertex.chi2(), theConversionVertex.ndof());
        if (chi2Prob<vtxChi2_)  highPurityPair=false;
      }

      //std::cout << "  highPurityPair after vertex cut " <<  highPurityPair << std::endl;
      std::vector<math::XYZPointF> trkPositionAtEcal;
      std::vector<reco::CaloClusterPtr> matchingBC;

      if (allowTrackBC_){//TODO find out the BC ptrs if not doing matching, otherwise, leave it empty
        //const int lbc_handle = bcHandleId[ll-allTracks.begin()],
        //        rbc_handle = bcHandleId[rr-allTracks.begin()];

        std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF>::const_iterator trackImpactPositionLeft = trackImpactPosition.find(ll->second);
        std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF>::const_iterator trackImpactPositionRight = trackImpactPosition.find(rr->second);
        std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr>::const_iterator trackMatchedBCLeft = trackMatchedBC.find(ll->second);        
        std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr>::const_iterator trackMatchedBCRight = trackMatchedBC.find(rr->second);        
        
        if (trackImpactPositionLeft!=trackImpactPosition.end()) {
          trkPositionAtEcal.push_back(trackImpactPositionLeft->second);//left track
        }
        else {
          trkPositionAtEcal.push_back(math::XYZPointF());//left track
        }
        if (trackImpactPositionRight!=trackImpactPosition.end()) {//second track ECAL position may be invalid
          trkPositionAtEcal.push_back(trackImpactPositionRight->second);
        }

        double total_e_bc = 0.;
        if (trackMatchedBCLeft!=trackMatchedBC.end()) {
          matchingBC.push_back(trackMatchedBCLeft->second);//left track
          total_e_bc += trackMatchedBCLeft->second->energy();
        }
        else {
          matchingBC.push_back( reco::CaloClusterPtr() );//left track
        }
        if (trackMatchedBCRight!=trackMatchedBC.end()) {//second track ECAL position may be invalid
          matchingBC.push_back(trackMatchedBCRight->second);
          total_e_bc += trackMatchedBCRight->second->energy();
        }
        
        if (total_e_bc<energyTotalBC_) {
          highPurityPair = false;
        }


      }
      //signature cuts, then check if vertex, then post-selection cuts
      highPurityPair = highPurityPair && track1HighPurity && track2HighPurity && goodVertex && checkPhi(left, right, trackerGeom, magField, theConversionVertex) ;


      /*
        for ( std::vector<edm::RefToBase<reco::Track> >::iterator iTk=trackPairRef.begin();  iTk!=trackPairRef.end(); iTk++) {
        math::XYZPointF impPos;
        if ( getTrackImpactPosition(*iTk, trackerGeom, magField, impPos) ) {
           
        }

        }
      */

      const float minAppDist = approachDist;
      reco::Conversion::ConversionAlgorithm algo = reco::Conversion::algoByName(algoName_);
      float dummy=0;
      reco::CaloClusterPtrVector scPtrVec;
      reco::Conversion  newCandidate(scPtrVec,  trackPairRef, trkPositionAtEcal, theConversionVertex, matchingBC, minAppDist,  trackInnPos, trackPin, trackPout, nHitsBeforeVtx, dlClosestHitToVtx, nSharedHits, dummy, algo );
      // Fill in scPtrVec with the macthing SC
      if ( matchingSC ( superClusterPtrs, newCandidate, scPtrVec) ) 
        newCandidate.setMatchingSuperCluster( scPtrVec);
          
          
      newCandidate.setQuality(reco::Conversion::highPurity,  highPurityPair);
      bool generalTracksOnly = ll->second->isTrackerOnly() && rr->second->isTrackerOnly() && !dynamic_cast<const reco::GsfTrack*>(ll->second->trackRef().get()) && !dynamic_cast<const reco::GsfTrack*>(rr->second->trackRef().get());
      bool gsfTracksOpenOnly = ll->second->isGsfTrackOpen() && rr->second->isGsfTrackOpen();
      bool arbitratedEcalSeeded = ll->second->isArbitratedEcalSeeded() && rr->second->isArbitratedEcalSeeded();
      bool arbitratedMerged = ll->second->isArbitratedMerged() && rr->second->isArbitratedMerged();
      bool arbitratedMergedEcalGeneral = ll->second->isArbitratedMergedEcalGeneral() && rr->second->isArbitratedMergedEcalGeneral();          
          
      newCandidate.setQuality(reco::Conversion::generalTracksOnly,  generalTracksOnly);
      newCandidate.setQuality(reco::Conversion::gsfTracksOpenOnly,  gsfTracksOpenOnly);
      newCandidate.setQuality(reco::Conversion::arbitratedEcalSeeded,  arbitratedEcalSeeded);
      newCandidate.setQuality(reco::Conversion::arbitratedMerged,  arbitratedMerged);
      newCandidate.setQuality(reco::Conversion::arbitratedMergedEcalGeneral,  arbitratedMergedEcalGeneral);          
          
      outputConvPhotonCollection.push_back(newCandidate);

    }
      
  }






}





//
// member functions
//

inline bool ConversionProducer::trackQualityFilter(const  edm::RefToBase<reco::Track>&   ref, bool isLeft){
  bool pass = true;
  if (isLeft){
    pass = (ref->normalizedChi2() < maxChi2Left_ && ref->found() >= minHitsLeft_);
  } else {
    pass = (ref->normalizedChi2() < maxChi2Right_ && ref->found() >= minHitsRight_);
  }

  return pass;
}

inline bool ConversionProducer::trackD0Cut(const  edm::RefToBase<reco::Track>&  ref){
  //NOTE if not allow d0 cut, always true
  return ((!allowD0_) || !(ref->d0()*ref->charge()/ref->d0Error()<d0Cut_));
}

inline bool ConversionProducer::trackD0Cut(const edm::RefToBase<reco::Track>&  ref, const reco::Vertex& the_pvtx){
  //
  return ((!allowD0_) || !(-ref->dxy(the_pvtx.position())*ref->charge()/ref->dxyError()<d0Cut_));
}


bool ConversionProducer::getTrackImpactPosition(const reco::Track* tk_ref,
                                                const TrackerGeometry* trackerGeom, const MagneticField* magField,
                                                math::XYZPointF& ew){

  PropagatorWithMaterial propag( alongMomentum, 0.000511, magField );
  
  ReferenceCountingPointer<Surface> ecalWall(
                                             new  BoundCylinder(129.f, GlobalPoint(0.,0.,0.), TkRotation<float>(),
                                                                 new SimpleCylinderBounds( 129, 129, -320.5, 320.5 ) ) );
  const float epsilon = 0.001;
  Surface::RotationType rot; // unit rotation matrix
  const float barrelRadius = 129.f;
  const float barrelHalfLength = 270.9f;
  const float endcapRadius = 171.1f;
  const float endcapZ = 320.5f;
  ReferenceCountingPointer<BoundCylinder>  theBarrel_(new BoundCylinder(barrelRadius, Surface::PositionType(0,0,0), rot,
                                                                         new SimpleCylinderBounds( barrelRadius-epsilon, barrelRadius+epsilon, 
-barrelHalfLength, barrelHalfLength)));
  ReferenceCountingPointer<BoundDisk>      theNegativeEtaEndcap_(
                                                                 new BoundDisk( Surface::PositionType( 0, 0, -endcapZ), rot,
                                                                                new SimpleDiskBounds( 0, endcapRadius, -epsilon, epsilon)));
  ReferenceCountingPointer<BoundDisk>      thePositiveEtaEndcap_(
                                                                 new BoundDisk( Surface::PositionType( 0, 0, endcapZ), rot,
                                                                                new SimpleDiskBounds( 0, endcapRadius, -epsilon, epsilon)));

  //const TrajectoryStateOnSurface myTSOS = trajectoryStateTransform::innerStateOnSurface(*(*ref), *trackerGeom, magField);
  const TrajectoryStateOnSurface myTSOS = trajectoryStateTransform::outerStateOnSurface(*tk_ref, *trackerGeom, magField);
  TrajectoryStateOnSurface  stateAtECAL;
  stateAtECAL = propag.propagate(myTSOS, *theBarrel_);
  if (!stateAtECAL.isValid() || ( stateAtECAL.isValid() && fabs(stateAtECAL.globalPosition().eta() ) >1.479f )  ) {
    //endcap propagator
    if (myTSOS.globalPosition().z() > 0.) {
        stateAtECAL = propag.propagate(myTSOS, *thePositiveEtaEndcap_);
    } else {
        stateAtECAL = propag.propagate(myTSOS, *theNegativeEtaEndcap_);
    }
  }       
  if (stateAtECAL.isValid()){
    ew = stateAtECAL.globalPosition();
    return true;
  }
  else
    return false;
}




bool ConversionProducer::matchingSC(const std::multimap<double, reco::CaloClusterPtr>& scMap, 
                                    reco::Conversion& aConv,
                                    // reco::CaloClusterPtr& mSC){
                                    reco::CaloClusterPtrVector& mSC) {

  //  double dRMin=999.;
  double detaMin=999.;
  double dphiMin=999.;                 
  reco::CaloClusterPtr match;
  for (std::multimap<double, reco::CaloClusterPtr>::const_iterator scItr = scMap.begin();  scItr != scMap.end(); scItr++) {
    const reco::CaloClusterPtr& sc = scItr->second; 
    const double delta_phi = reco::deltaPhi( aConv.refittedPairMomentum().phi(), sc->phi());
    double sceta = sc->eta();
    double conveta = etaTransformation(aConv.refittedPairMomentum().eta(), aConv.zOfPrimaryVertexFromTracks() );
    const double delta_eta = fabs(conveta - sceta);
    if ( fabs(delta_eta) < fabs(detaMin) && fabs(delta_phi) < fabs(dphiMin) ) {
      detaMin=  fabs(delta_eta);
      dphiMin=  fabs(delta_phi);
      match=sc;
    }
  }
  
  if ( fabs(detaMin) < dEtacutForSCmatching_ && fabs(dphiMin) < dPhicutForSCmatching_ ) {
    mSC.push_back(match);
    return true;
  } else 
    return false;
}

bool ConversionProducer::getMatchedBC(const std::multimap<double, reco::CaloClusterPtr>& bcMap, 
                                      const math::XYZPointF& trackImpactPosition,
                                      reco::CaloClusterPtr& closestBC){
  const double track_eta = trackImpactPosition.eta();
  const double track_phi = trackImpactPosition.phi();

  double min_eta = 999., min_phi = 999.;
  reco::CaloClusterPtr closest_bc;
  for (std::multimap<double, reco::CaloClusterPtr>::const_iterator bc = bcMap.lower_bound(track_eta - halfWayEta_);
       bc != bcMap.upper_bound(track_eta + halfWayEta_); ++bc){//use eta map to select possible BC collection then loop in
    const reco::CaloClusterPtr& ebc = bc->second;
    const double delta_eta = track_eta-(ebc->position().eta());
    const double delta_phi = reco::deltaPhi(track_phi, (ebc->position().phi()));
    if (fabs(delta_eta)<dEtaTkBC_ && fabs(delta_phi)<dPhiTkBC_){
        if (fabs(min_eta)>fabs(delta_eta) && fabs(min_phi)>fabs(delta_phi)){//take the closest to track BC
        min_eta = delta_eta;
        min_phi = delta_phi;
        closest_bc = bc->second;
        //TODO check if min_eta>delta_eta but min_phi<delta_phi
        }
    }
  }

  if (min_eta < 999.){
    closestBC = closest_bc;
    return true;
  } else
    return false;
}






//check track open angle of phi at vertex
bool ConversionProducer::checkPhi(const edm::RefToBase<reco::Track>& tk_l, const edm::RefToBase<reco::Track>& tk_r,
                                  const TrackerGeometry* trackerGeom, const MagneticField* magField,
                                  const reco::Vertex& vtx){
  if (!allowDeltaPhi_)
    return true;
  //if track has innermost momentum, check with innermost phi
  //if track also has valid vertex, propagate to vertex then calculate phi there
  //if track has no innermost momentum, just return true, because track->phi() makes no sense
  if (tk_l->extra().isNonnull() && tk_r->extra().isNonnull()){
    double iphi1 = tk_l->innerMomentum().phi(), iphi2 = tk_r->innerMomentum().phi();
    if (vtx.isValid()){
        PropagatorWithMaterial propag( anyDirection, 0.000511, magField );
        
        double recoPhoR = vtx.position().Rho();
        Surface::RotationType rot;
        ReferenceCountingPointer<BoundCylinder>  theBarrel_(new BoundCylinder(recoPhoR, Surface::PositionType(0,0,0), rot,
                                                                 new SimpleCylinderBounds( recoPhoR-0.001, recoPhoR+0.001, 
                                                                 -fabs(vtx.position().z()), fabs(vtx.position().z()))));
        ReferenceCountingPointer<BoundDisk>      theDisk_(
                                                        new BoundDisk( Surface::PositionType( 0, 0, vtx.position().z()), rot,
                                                                       new SimpleDiskBounds( 0, recoPhoR, -0.001, 0.001)));

        const TrajectoryStateOnSurface myTSOS1 = trajectoryStateTransform::innerStateOnSurface(*tk_l, *trackerGeom, magField);
        const TrajectoryStateOnSurface myTSOS2 = trajectoryStateTransform::innerStateOnSurface(*tk_r, *trackerGeom, magField);
        TrajectoryStateOnSurface  stateAtVtx1, stateAtVtx2;
        stateAtVtx1 = propag.propagate(myTSOS1, *theBarrel_);
        if (!stateAtVtx1.isValid() ) {
        stateAtVtx1 = propag.propagate(myTSOS1, *theDisk_);
        }
        if (stateAtVtx1.isValid()){
        iphi1 = stateAtVtx1.globalDirection().phi();
        }
        stateAtVtx2 = propag.propagate(myTSOS2, *theBarrel_);
        if (!stateAtVtx2.isValid() ) {
        stateAtVtx2 = propag.propagate(myTSOS2, *theDisk_);
        }
        if (stateAtVtx2.isValid()){
        iphi2 = stateAtVtx2.globalDirection().phi();
        }
    }
    const double dPhi = reco::deltaPhi(iphi1, iphi2);
    return (fabs(dPhi) < deltaPhi_);
  } else {
    return true;
  }
}

bool ConversionProducer::preselectTrackPair(const reco::TransientTrack &ttk_l, const reco::TransientTrack &ttk_r,
                                            double& appDist) {
  

  double dCotTheta =  1./tan(ttk_l.track().innerMomentum().theta()) - 1./tan(ttk_r.track().innerMomentum().theta());
  if (allowDeltaCot_ && (std::abs(dCotTheta) > deltaCotTheta_)) {
    return false;
  }
  
  //non-conversion hypothesis, reject prompt track pairs
  ClosestApproachInRPhi closest;
  closest.calculate(ttk_l.innermostMeasurementState(),ttk_r.innermostMeasurementState());
  if (!closest.status()) {
    return false;
  }
  
  if (closest.crossingPoint().perp() < r_cut) {
    return false;
  }

  
  //compute tangent point btw tracks (conversion hypothesis)
  TangentApproachInRPhi tangent;
  tangent.calculate(ttk_l.innermostMeasurementState(),ttk_r.innermostMeasurementState());
  if (!tangent.status()) {
    return false;
  }
  
  GlobalPoint tangentPoint = tangent.crossingPoint();
  double rho = tangentPoint.perp();
  
  //reject candidates well outside of tracker bounds
  if (rho > maxTrackRho_) {
    return false;
  }
  
  if (std::abs(tangentPoint.z()) > maxTrackZ_) {
    return false;
  }
  
  std::pair<GlobalTrajectoryParameters,GlobalTrajectoryParameters> trajs = tangent.trajectoryParameters();
  
  //very large separation in z, no hope
  if (std::abs(trajs.first.position().z() - trajs.second.position().z()) > dzCut_) {
    return false;
  }
  
  
  float minApproach = tangent.perpdist();
  appDist = minApproach;
  
  if (allowMinApproach_ && (minApproach < minApproachLow_ || minApproach > minApproachHigh_) ) {
    return false;
  }
  
  return true;
  
  
}

bool ConversionProducer::checkTrackPair(const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& ll, 
                                        const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& rr){

  const reco::CaloClusterPtr& bc_l = ll.second;//can be null, so check isNonnull()
  const reco::CaloClusterPtr& bc_r = rr.second;
    
  //The cuts should be ordered by considering if takes time and if cuts off many fakes
  if (allowTrackBC_){
    //check energy of BC
    double total_e_bc = 0;
    if (bc_l.isNonnull()) total_e_bc += bc_l->energy();
    if (rightBC_) 
        if (bc_r.isNonnull())
        total_e_bc += bc_r->energy();

    if (total_e_bc < energyTotalBC_) return false;
  }

  return true;
}



//because reco::vertex uses track ref, so have to keep them
bool ConversionProducer::checkVertex(const reco::TransientTrack &ttk_l, const reco::TransientTrack &ttk_r, 
                                     const MagneticField* magField,
                                     reco::Vertex& the_vertex){
  bool found = false;

  std::vector<reco::TransientTrack>  pair;
  pair.push_back(ttk_l);
  pair.push_back(ttk_r);
   
  found = theVertexFinder_->run(pair, the_vertex);



  return found;
}



double ConversionProducer::etaTransformation(  float EtaParticle , float Zvertex)  {

  //---Definitions
  const float PI    = 3.1415927;

  //---Definitions for ECAL
  const float R_ECAL           = 136.5;
  const float Z_Endcap         = 328.0;
  const float etaBarrelEndcap  = 1.479; 
   
  //---ETA correction

  float Theta = 0.0  ; 
  float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;

  if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
  if(Theta<0.0) Theta = Theta+PI ;
  double ETA = - log(tan(0.5*Theta));
         
  if( fabs(ETA) > etaBarrelEndcap )
    {
      float Zend = Z_Endcap ;
      if(EtaParticle<0.0 )  Zend = -Zend ;
      float Zlen = Zend - Zvertex ;
      float RR = Zlen/sinh(EtaParticle); 
      Theta = atan(RR/Zend);
      if(Theta<0.0) Theta = Theta+PI ;
      ETA = - log(tan(0.5*Theta));            
    } 
  //---Return the result
  return ETA;
  //---end
}