ConvertedPhotonProducer.cc

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include <memory>

// Framework
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
//
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/CaloRecHit/interface/CaloCluster.h"

//
#include "DataFormats/EgammaTrackReco/interface/TrackCaloClusterAssociation.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
#include "DataFormats/EgammaCandidates/interface/Conversion.h"
#include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
//
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
//
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
//
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackEcalImpactPoint.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackPairFinder.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionVertexFinder.h"
#include "RecoEgamma/EgammaPhotonProducers/interface/ConvertedPhotonProducer.h"
//
#include "RecoTracker/Record/interface/CkfComponentsRecord.h"
#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
//
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"

ConvertedPhotonProducer::ConvertedPhotonProducer(const edm::ParameterSet& config) : 
  conf_(config), 
  theTrackPairFinder_(0), 
  theVertexFinder_(0), 
  theLikelihoodCalc_(0)
{


  
  //cout<< " ConvertedPhotonProducer CTOR " << "\n";
  
  
  
  // use onfiguration file to setup input collection names
  bcBarrelCollection_     = conf_.getParameter<edm::InputTag>("bcBarrelCollection");
  bcEndcapCollection_     = conf_.getParameter<edm::InputTag>("bcEndcapCollection");
  
  scHybridBarrelProducer_       = conf_.getParameter<edm::InputTag>("scHybridBarrelProducer");
  scIslandEndcapProducer_       = conf_.getParameter<edm::InputTag>("scIslandEndcapProducer");

  conversionOITrackProducer_ = conf_.getParameter<std::string>("conversionOITrackProducer");
  conversionIOTrackProducer_ = conf_.getParameter<std::string>("conversionIOTrackProducer");
  
  outInTrackSCAssociationCollection_ = conf_.getParameter<std::string>("outInTrackSCAssociation");
  inOutTrackSCAssociationCollection_ = conf_.getParameter<std::string>("inOutTrackSCAssociation");
  
  algoName_ = conf_.getParameter<std::string>( "AlgorithmName" );  

  hcalTowers_ = conf_.getParameter<edm::InputTag>("hcalTowers");
  hOverEConeSize_   = conf_.getParameter<double>("hOverEConeSize");
  maxHOverE_        = conf_.getParameter<double>("maxHOverE");
  minSCEt_        = conf_.getParameter<double>("minSCEt");
  recoverOneTrackCase_ = conf_.getParameter<bool>( "recoverOneTrackCase" );  
  dRForConversionRecovery_ = conf_.getParameter<double>("dRForConversionRecovery");
  deltaCotCut_ = conf_.getParameter<double>("deltaCotCut");
  minApproachDisCut_ = conf_.getParameter<double>("minApproachDisCut");

  maxNumOfCandidates_        = conf_.getParameter<int>("maxNumOfCandidates");
  risolveAmbiguity_ = conf_.getParameter<bool>("risolveConversionAmbiguity");  
  likelihoodWeights_= conf_.getParameter<std::string>("MVA_weights_location");
  
 
  // use onfiguration file to setup output collection names
  ConvertedPhotonCollection_     = conf_.getParameter<std::string>("convertedPhotonCollection");
  CleanedConvertedPhotonCollection_     = conf_.getParameter<std::string>("cleanedConvertedPhotonCollection");
  
  
  // Register the product
  produces< reco::ConversionCollection >(ConvertedPhotonCollection_);
  produces< reco::ConversionCollection >(CleanedConvertedPhotonCollection_);
  
  // instantiate the Track Pair Finder algorithm
  theTrackPairFinder_ = new ConversionTrackPairFinder ();
  edm::FileInPath path_mvaWeightFile(likelihoodWeights_.c_str() );      
  theLikelihoodCalc_ = new ConversionLikelihoodCalculator();
  theLikelihoodCalc_->setWeightsFile(path_mvaWeightFile.fullPath().c_str());
  // instantiate the Vertex Finder algorithm
  theVertexFinder_ = new ConversionVertexFinder ( conf_);


  // Inizilize my global event counter
  nEvt_=0;
  
}

ConvertedPhotonProducer::~ConvertedPhotonProducer() {
  delete theTrackPairFinder_;
  delete theLikelihoodCalc_;
  delete theVertexFinder_;
}



void  ConvertedPhotonProducer::beginRun (edm::Run const& r, edm::EventSetup const & theEventSetup) {
 

    //get magnetic field
  //edm::LogInfo("ConvertedPhotonProducer") << " get magnetic field" << "\n";
  theEventSetup.get<IdealMagneticFieldRecord>().get(theMF_);  

  // Transform Track into TransientTrack (needed by the Vertex fitter)
  theEventSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theTransientTrackBuilder_);

  
}



void ConvertedPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
  
  using namespace edm;
  nEvt_++;  


  //  LogDebug("ConvertedPhotonProducer")   << "ConvertedPhotonProduce::produce event number " <<   theEvent.id() << " Global counter " << nEvt_ << "\n";
  //  std::cout    << "ConvertedPhotonProduce::produce event number " <<   theEvent.id() << " Global counter " << nEvt_ << "\n";
  
  //
  // create empty output collections
  //
  // Converted photon candidates
  reco::ConversionCollection outputConvPhotonCollection;
  std::auto_ptr<reco::ConversionCollection> outputConvPhotonCollection_p(new reco::ConversionCollection);
  // Converted photon candidates
  reco::ConversionCollection cleanedConversionCollection;
  std::auto_ptr<reco::ConversionCollection> cleanedConversionCollection_p(new reco::ConversionCollection);

  
  // Get the Super Cluster collection in the Barrel
  bool validBarrelSCHandle=true;
  edm::Handle<edm::View<reco::CaloCluster> > scBarrelHandle;
  theEvent.getByLabel(scHybridBarrelProducer_,scBarrelHandle);
  if (!scBarrelHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product "<<scHybridBarrelProducer_.label();
    validBarrelSCHandle=false;
  }
   
  // Get the Super Cluster collection in the Endcap
  bool validEndcapSCHandle=true;
  edm::Handle<edm::View<reco::CaloCluster> > scEndcapHandle;
  theEvent.getByLabel(scIslandEndcapProducer_,scEndcapHandle);
  if (!scEndcapHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product "<<scIslandEndcapProducer_.label();
    validEndcapSCHandle=false;
  }
  
    
  bool validTrackInputs=true;
  Handle<reco::TrackCollection> outInTrkHandle;
  theEvent.getByLabel(conversionOITrackProducer_,  outInTrkHandle);
  if (!outInTrkHandle.isValid()) {
    //std::cout << "Error! Can't get the conversionOITrack " << "\n";
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionOITrack " << "\n";
    validTrackInputs=false;
  }
  //  LogDebug("ConvertedPhotonProducer")<< "ConvertedPhotonProducer  outInTrack collection size " << (*outInTrkHandle).size() << "\n";
  
   
  Handle<reco::TrackCaloClusterPtrAssociation> outInTrkSCAssocHandle;
  theEvent.getByLabel( conversionOITrackProducer_, outInTrackSCAssociationCollection_, outInTrkSCAssocHandle);
  if (!outInTrkSCAssocHandle.isValid()) {
    //  std::cout << "Error! Can't get the product " <<  outInTrackSCAssociationCollection_.c_str() <<"\n";
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product " <<  outInTrackSCAssociationCollection_.c_str() <<"\n";
    validTrackInputs=false;
  }

  Handle<reco::TrackCollection> inOutTrkHandle;
  theEvent.getByLabel(conversionIOTrackProducer_, inOutTrkHandle);
  if (!inOutTrkHandle.isValid()) {
    // std::cout << "Error! Can't get the conversionIOTrack " << "\n";
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionIOTrack " << "\n";
    validTrackInputs=false;
  }
  //  LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer inOutTrack collection size " << (*inOutTrkHandle).size() << "\n";



  Handle<reco::TrackCollection> generalTrkHandle;
  if (  recoverOneTrackCase_ ) {
    theEvent.getByLabel("generalTracks", generalTrkHandle);
    if (!generalTrkHandle.isValid()) {
      //std::cout << "Error! Can't get the genralTracks " << "\n";
      edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the genralTracks " << "\n";
    }
  }  
  
  Handle<reco::TrackCaloClusterPtrAssociation> inOutTrkSCAssocHandle;
  theEvent.getByLabel( conversionIOTrackProducer_, inOutTrackSCAssociationCollection_, inOutTrkSCAssocHandle);
  if (!inOutTrkSCAssocHandle.isValid()) {
    //std::cout << "Error! Can't get the product " <<  inOutTrackSCAssociationCollection_.c_str() <<"\n";
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product " <<  inOutTrackSCAssociationCollection_.c_str() <<"\n";
    validTrackInputs=false;
  }


  

  // Get the basic cluster collection in the Barrel 
  edm::Handle<edm::View<reco::CaloCluster> > bcBarrelHandle;
  theEvent.getByLabel( bcBarrelCollection_, bcBarrelHandle);
  if (!bcBarrelHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product "<<bcBarrelCollection_.label();
  }

    
  // Get the basic cluster collection in the Endcap 
  edm::Handle<edm::View<reco::CaloCluster> > bcEndcapHandle;
  theEvent.getByLabel( bcEndcapCollection_, bcEndcapHandle);
  if (!bcEndcapHandle.isValid()) {
    edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the product "<<bcEndcapCollection_.label();
  }
 

// get Hcal towers collection 
  Handle<CaloTowerCollection> hcalTowersHandle;
  theEvent.getByLabel(hcalTowers_, hcalTowersHandle);

  // get the geometry from the event setup:
  theEventSetup.get<CaloGeometryRecord>().get(theCaloGeom_);


  if (  validTrackInputs ) {
    //do the conversion:
    std::vector<reco::TransientTrack> t_outInTrk = ( *theTransientTrackBuilder_ ).build(outInTrkHandle );
    std::vector<reco::TransientTrack> t_inOutTrk = ( *theTransientTrackBuilder_ ).build(inOutTrkHandle );
    
    
    std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairs;
    allPairs = theTrackPairFinder_->run(t_outInTrk, outInTrkHandle, outInTrkSCAssocHandle, t_inOutTrk, inOutTrkHandle, inOutTrkSCAssocHandle  );
    //LogDebug("ConvertedPhotonProducer")  << "ConvertedPhotonProducer  allPairs.size " << allPairs.size() << "\n";      

    buildCollections(theEventSetup, scBarrelHandle, bcBarrelHandle, hcalTowersHandle, generalTrkHandle, allPairs, outputConvPhotonCollection);
    buildCollections(theEventSetup, scEndcapHandle, bcEndcapHandle, hcalTowersHandle, generalTrkHandle, allPairs, outputConvPhotonCollection);
  }
  
  // put the product in the event
  outputConvPhotonCollection_p->assign(outputConvPhotonCollection.begin(),outputConvPhotonCollection.end());
  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Putting in the event    converted photon candidates " << (*outputConvPhotonCollection_p).size() << "\n";  
  const edm::OrphanHandle<reco::ConversionCollection> conversionHandle= theEvent.put( outputConvPhotonCollection_p, ConvertedPhotonCollection_);


  // Loop over barrel and endcap SC collections and fill the  photon collection
  if ( validBarrelSCHandle) cleanCollections(scBarrelHandle,
                         conversionHandle,
                         cleanedConversionCollection);
  if ( validEndcapSCHandle) cleanCollections(scEndcapHandle,
                         conversionHandle,
                         cleanedConversionCollection);
                        

  cleanedConversionCollection_p->assign(cleanedConversionCollection.begin(),cleanedConversionCollection.end());   
  theEvent.put( cleanedConversionCollection_p, CleanedConvertedPhotonCollection_);

  
}


void ConvertedPhotonProducer::buildCollections ( edm::EventSetup const & es,  
                         const edm::Handle<edm::View<reco::CaloCluster> > & scHandle,
                         const edm::Handle<edm::View<reco::CaloCluster> > & bcHandle,
                         const edm::Handle<CaloTowerCollection> & hcalTowersHandle, 
                         const edm::Handle<reco::TrackCollection>  & generalTrkHandle,
                         std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
                         reco::ConversionCollection & outputConvPhotonCollection)

{

 // instantiate the algorithm for finding the position of the track extrapolation at the Ecal front face
  ConversionTrackEcalImpactPoint theEcalImpactPositionFinder( &(*theMF_) );
  

  reco::Conversion::ConversionAlgorithm algo = reco::Conversion::algoByName(algoName_);
 
  std::vector<reco::TransientTrack> t_generalTrk;
  if (  recoverOneTrackCase_ )  t_generalTrk = ( *theTransientTrackBuilder_ ).build(generalTrkHandle );
  //const CaloGeometry* geometry = theCaloGeom_.product(); 

  //  Loop over SC in the barrel and reconstruct converted photons
  int myCands=0;
  reco::CaloClusterPtrVector scPtrVec;
  for (unsigned i = 0; i < scHandle->size(); ++i ) {
    reco::CaloClusterPtr aClus= scHandle->ptrAt(i);

    // preselection based in Et and H/E cut
    if (aClus->energy()/cosh(aClus->eta()) <= minSCEt_) continue;
    const reco::CaloCluster* pClus=&(*aClus);
    const reco::SuperCluster*  sc=dynamic_cast<const reco::SuperCluster*>(pClus);
    const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
    EgammaTowerIsolation towerIso(hOverEConeSize_,0.,0.,-1,hcalTowersColl) ;
    double HoE=towerIso.getTowerESum(sc)/sc->energy();
    if (HoE>=maxHOverE_)  continue;


    std::vector<edm::Ref<reco::TrackCollection> > trackPairRef;
    std::vector<math::XYZPointF> trackInnPos;
    std::vector<math::XYZVectorF> trackPin;
    std::vector<math::XYZVectorF> trackPout;
    float minAppDist=-99;

    //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer SC energy " << aClus->energy() << " eta " <<  aClus->eta() << " phi " <<  aClus->phi() << "\n";

    
    const reco::Particle::Point  vtx( 0, 0, 0 );

   
    math::XYZVector direction =aClus->position() - vtx;
    math::XYZVector momentum = direction.unit() * aClus->energy();
    const reco::Particle::LorentzVector  p4(momentum.x(), momentum.y(), momentum.z(), aClus->energy() );
    
    int nFound=0;    
    if ( allPairs.size() ) {

      nFound=0;


      for (  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairs.begin(); iPair!= allPairs.end(); ++iPair ) {
    scPtrVec.clear();

    reco::Vertex  theConversionVertex;       
    reco::CaloClusterPtr caloPtr=iPair->second;
    if ( !( aClus == caloPtr ) ) continue;
            
        scPtrVec.push_back(aClus);     
    nFound++;

    std::vector<math::XYZPointF> trkPositionAtEcal = theEcalImpactPositionFinder.find(  iPair->first, bcHandle );
    std::vector<reco::CaloClusterPtr>  matchingBC = theEcalImpactPositionFinder.matchingBC();
    

        minAppDist=-99;
    const std::string metname = "ConvertedPhotons|ConvertedPhotonProducer";
    if ( (iPair->first).size()  > 1 ) {
      try{
        
        theVertexFinder_->run(iPair->first, theConversionVertex ); 
        
        
      }
      catch ( cms::Exception& e ) {
        //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
        edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
                     << e.explainSelf();
        
      }
     
      // Old TwoTrackMinimumDistance md;
      // Old md.calculate  (  (iPair->first)[0].initialFreeState(),  (iPair->first)[1].initialFreeState() );
          // Old minAppDist = md.distance(); 
 
    
    

        
    

    /*
    for ( unsigned int i=0; i< matchingBC.size(); ++i) {
          if (  matchingBC[i].isNull() )  std::cout << " This ref to BC is null: skipping " <<  "\n";
          else 
        std::cout << " BC energy " << matchingBC[i]->energy() <<  "\n";
    }
    */


    trackPairRef.clear();
        trackInnPos.clear();
    trackPin.clear();
    trackPout.clear();
      
    
    for ( std::vector<reco::TransientTrack>::const_iterator iTk=(iPair->first).begin(); iTk!= (iPair->first).end(); ++iTk) {
      //LogDebug("ConvertedPhotonProducer")  << "  ConvertedPhotonProducer Transient Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
      
      const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
      reco::TrackRef myTkRef= ttt->persistentTrackRef(); 
      
      //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer Ref to Rec Tracks in the pair  charge " << myTkRef->charge() << " Num of RecHits " << myTkRef->recHitsSize() << " inner momentum " << myTkRef->innerMomentum() << "\n";  
      if ( myTkRef->extra().isNonnull() ) {
        trackInnPos.push_back(  toFConverterP(myTkRef->innerPosition()));
        trackPin.push_back(  toFConverterV( myTkRef->innerMomentum()));
        trackPout.push_back(  toFConverterV(myTkRef->outerMomentum()));
      }
      trackPairRef.push_back(myTkRef);
      
    }
    
    //  std::cout << " ConvertedPhotonProducer trackPin size " << trackPin.size() << std::endl;
    //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer SC energy " <<  aClus->energy() << "\n";
    //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer photon p4 " << p4  << "\n";
    //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer vtx " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
        if( theConversionVertex.isValid() ) {
      
      //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer theConversionVertex " << theConversionVertex.position().x() << " " << theConversionVertex.position().y() << " " << theConversionVertex.position().z() << "\n";
      
    }
    //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer trackPairRef  " << trackPairRef.size() <<  "\n";

        
    minAppDist=calculateMinApproachDistance( trackPairRef[0],  trackPairRef[1]);

    double like = -999.;
    reco::Conversion  newCandidate(scPtrVec,  trackPairRef,  trkPositionAtEcal, theConversionVertex, matchingBC,  minAppDist, trackInnPos, trackPin, trackPout, like, algo);
//    like = theLikelihoodCalc_->calculateLikelihood(newCandidate, es );
    like = theLikelihoodCalc_->calculateLikelihood( newCandidate );
//    std::cout << "like = " << like << std::endl;
    newCandidate.setMVAout(like);
    outputConvPhotonCollection.push_back(newCandidate);
    
    
    myCands++;
    //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Put the ConvertedPhotonCollection a candidate in the Barrel " << "\n";
    
    } else {
      
      
      //      std::cout << "   ConvertedPhotonProducer case with only one track found " <<  "\n";
 
        //std::cout << "   ConvertedPhotonProducer recovering one track " <<  "\n";
        trackPairRef.clear();
        trackInnPos.clear();
        trackPin.clear();
        trackPout.clear();
        std::vector<reco::TransientTrack>::const_iterator iTk=(iPair->first).begin();
        //std::cout  << "  ConvertedPhotonProducer Transient Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << " pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";     
        const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
        reco::TrackRef myTk= ttt->persistentTrackRef(); 
        if ( myTk->extra().isNonnull() ) {
          trackInnPos.push_back(  toFConverterP(myTk->innerPosition()));
          trackPin.push_back(  toFConverterV(myTk->innerMomentum()));
          trackPout.push_back(  toFConverterV(myTk->outerMomentum()));
        }
        trackPairRef.push_back(myTk);
        //std::cout << " Provenance " << myTk->algoName() << std::endl;
    
        if (  recoverOneTrackCase_ ) {    
          float theta1 = myTk->innerMomentum().Theta();
          float dCot=999.;
          float dCotTheta=-999.;
          reco::TrackRef goodRef;
          std::vector<reco::TransientTrack>::const_iterator iGoodGenTran;
          for ( std::vector<reco::TransientTrack>::const_iterator iTran= t_generalTrk.begin(); iTran != t_generalTrk.end(); ++iTran) {
        const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTran->basicTransientTrack());
        reco::TrackRef trRef= ttt->persistentTrackRef(); 
        if ( trRef->charge()*myTk->charge() > 0 ) continue;
        float dEta =  trRef->eta() - myTk->eta();
        float dPhi =  trRef->phi() - myTk->phi();
        if ( sqrt (dEta*dEta + dPhi*dPhi) > dRForConversionRecovery_ ) continue; 
        float theta2 = trRef->innerMomentum().Theta();
        dCotTheta =  1./tan(theta1) - 1./tan(theta2) ;
        //  std::cout << "  ConvertedPhotonProducer recovering general transient track charge " << trRef->charge() << " momentum " << trRef->innerMomentum() << " dcotTheta " << fabs(dCotTheta) << std::endl;
        if ( fabs(dCotTheta) < dCot ) {
          dCot = fabs(dCotTheta);
          goodRef = trRef;
          iGoodGenTran=iTran;
        }
          }
          
          if ( goodRef.isNonnull() ) {
        
        minAppDist=calculateMinApproachDistance( myTk, goodRef);
        
        // std::cout << "  ConvertedPhotonProducer chosen dCotTheta " <<  fabs(dCotTheta) << std::endl;
        if ( fabs(dCotTheta) < deltaCotCut_ && minAppDist > minApproachDisCut_ ) {
          trackInnPos.push_back(  toFConverterP(goodRef->innerPosition())); 
          trackPin.push_back(  toFConverterV(goodRef->innerMomentum()));
          trackPout.push_back(  toFConverterV(goodRef->outerMomentum()));
          trackPairRef.push_back( goodRef );
          //        std::cout << " ConvertedPhotonProducer adding opposite charge track from generalTrackCollection charge " <<  goodRef ->charge() << " pt " << sqrt(goodRef->innerMomentum().perp2())  << " trackPairRef size " << trackPairRef.size() << std::endl;            
          //std::cout << " Track Provenenance " << goodRef->algoName() << std::endl; 
          std::vector<reco::TransientTrack> mypair;
          mypair.push_back(*iTk); 
          mypair.push_back(*iGoodGenTran); 
          
          try{
            theVertexFinder_->run(iPair->first, theConversionVertex ); 
            
          }
          catch ( cms::Exception& e ) {
            //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
            edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
                         << e.explainSelf();
            
          }
        } 
        
          }     
          
        } // bool On/Off one track case recovery using generalTracks  
        double like = -999.;
        reco::Conversion  newCandidate(scPtrVec,  trackPairRef,  trkPositionAtEcal, theConversionVertex, matchingBC,  minAppDist, trackInnPos, trackPin, trackPout, like, algo);
        like = theLikelihoodCalc_->calculateLikelihood(newCandidate);
        newCandidate.setMVAout(like);
        outputConvPhotonCollection.push_back(newCandidate);
          
          
          
        
    } // case with only on track: looking in general tracks
    
    
    
    
      } 
      
    }
    


    
  }
  




}


void ConvertedPhotonProducer::cleanCollections(const edm::Handle<edm::View<reco::CaloCluster> > & scHandle,
                          const edm::OrphanHandle<reco::ConversionCollection> & conversionHandle,
                          reco::ConversionCollection & outputConversionCollection) {


  reco::Conversion* newCandidate=0;
  for(unsigned int lSC=0; lSC < scHandle->size(); lSC++) {
    
    // get pointer to SC
    reco::CaloClusterPtr aClus= scHandle->ptrAt(lSC);    
        
    // SC energy preselection
    if (aClus->energy()/cosh(aClus->eta()) <= minSCEt_) continue;


    if (  conversionHandle.isValid() ) {

      if ( risolveAmbiguity_ ) {
    std::vector<reco::ConversionRef> bestRef=solveAmbiguity( conversionHandle , aClus);
   
        for ( std::vector<reco::ConversionRef>::iterator iRef=bestRef.begin(); iRef!=bestRef.end(); iRef++ ) {
          if ( iRef->isNonnull() )  {
        newCandidate= (*iRef)->clone();
        outputConversionCollection.push_back(*newCandidate);
            delete newCandidate;
      }
    }
    
      } else {
    

    for( unsigned int icp = 0;  icp < conversionHandle->size(); icp++) {
      reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle,icp));
      if (!( aClus.id() == cpRef->caloCluster()[0].id() && aClus.key() == cpRef->caloCluster()[0].key() )) continue; 
      if ( !cpRef->isConverted() ) continue;  
      if ( cpRef->nTracks() <2 ) continue;
      newCandidate= (&(*cpRef))->clone();    
      outputConversionCollection.push_back(*newCandidate);
      delete newCandidate;

    }     
    
      } // solve or not the ambiguity of many conversion candidates 
    
    }     

   
  }    
}




std::vector<reco::ConversionRef>  ConvertedPhotonProducer::solveAmbiguity(const edm::OrphanHandle<reco::ConversionCollection> & conversionHandle, reco::CaloClusterPtr& scRef) {
  std::multimap<double, reco::ConversionRef, std::greater<double> >   convMap;

  for ( unsigned int icp=0; icp< conversionHandle->size(); icp++) {
    reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle,icp));
    
    //std::cout << " cpRef " << cpRef->nTracks() << " " <<  cpRef ->caloCluster()[0]->energy() << std::endl;    
    if (!( scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key() )) continue;    
    if ( !cpRef->isConverted() ) continue;  
    double like = cpRef->MVAout();
    if ( cpRef->nTracks() <2 ) continue;
    //    std::cout << " Like " << like << std::endl;
    convMap.insert ( std::make_pair(like,cpRef) );
    
  }          
  
  //  std::cout << " convMap size " << convMap.size() << std::endl;

  std::multimap<double, reco::ConversionRef>::iterator  iMap; 
  std::vector<reco::ConversionRef> bestRefs;
  for (iMap=convMap.begin();  iMap!=convMap.end(); iMap++) {
    //    std::cout << " Like list in the map " <<  iMap->first << " " << (iMap->second)->EoverP() << std::endl;
    bestRefs.push_back(  iMap->second );
    if (  int(bestRefs.size()) ==  maxNumOfCandidates_ ) break;   
  }


   return bestRefs;
  
  
} 





float ConvertedPhotonProducer::calculateMinApproachDistance ( const reco::TrackRef& track1, const reco::TrackRef& track2) {
  float dist=9999.;

  double x1, x2, y1, y2;
  double xx_1 = track1->innerPosition().x(), yy_1 = track1->innerPosition().y(), zz_1 = track1->innerPosition().z();
  double xx_2 = track2->innerPosition().x(), yy_2 = track2->innerPosition().y(), zz_2 = track2->innerPosition().z();
  double radius1 = track1->innerMomentum().Rho()/(.3*(theMF_->inTesla(GlobalPoint(xx_1, yy_1, zz_1)).z()))*100;
  double radius2 = track2->innerMomentum().Rho()/(.3*(theMF_->inTesla(GlobalPoint(xx_2, yy_2, zz_2)).z()))*100;
  getCircleCenter(track1, radius1, x1, y1);
  getCircleCenter(track2, radius2, x2, y2);
  dist = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)) - radius1 - radius2;
  
 return dist;

} 


void ConvertedPhotonProducer::getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0){
  double x1, y1, phi;
  x1 = tk->innerPosition().x();//inner position and inner momentum need track Extra!
  y1 = tk->innerPosition().y();
  phi = tk->innerMomentum().phi();
  const int charge = tk->charge();
  x0 = x1 + r*sin(phi)*charge;
  y0 = y1 - r*cos(phi)*charge;

}