#include <PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h>

Public Member Functions
void	analyze (const edm::Event &event, const edm::EventSetup &setup)
void	clear ()
IdealHelixParameters *	getIdealHelixParameters ()
TrajectorySeedCollection *	getTrajectorySeedCollection ()
TrajectorySeedCollection *	getTrajectorySeedCollectionOfSourceTracks ()
void	init ()
	PhotonConversionTrajectorySeedProducerFromSingleLegAlgo (const edm::ParameterSet &)
	~PhotonConversionTrajectorySeedProducerFromSingleLegAlgo ()
Private Types
typedef Regions::const_iterator	IR
typedef std::vector < TrackingRegion * >	Regions
Private Member Functions
bool	inspectTrack (const reco::Track *track, const TrackingRegion &region, math::XYZPoint &primaryVertexPoint)
void	loopOnPriVtx (const reco::Track &tk, const std::vector< reco::Vertex > &selectedPriVtxCompatibleWithTrack)
void	loopOnTracks ()
bool	rejectTrack (const reco::Track &track)
bool	selectPriVtxCompatibleWithTrack (const reco::Track &tk, std::vector< reco::Vertex > &selectedPriVtxCompatibleWithTrack)
Private Attributes
bool	_applyTkVtxConstraint
edm::InputTag	_beamSpotInputTag
const edm::ParameterSet	_conf
size_t	_countSeedTracks
IdealHelixParameters	_IdealHelixParameters
double	_maxDZSigmas
size_t	_maxNumSelVtx
edm::InputTag	_primaryVtxInputTag
double	_vtxMinDoF
edm::ParameterSet	creatorPSet
edm::ParameterSet	hitsfactoryPSet
const MagneticField *	magField
const edm::EventSetup *	myEsetup
const edm::Event *	myEvent
PrintRecoObjects	po
reco::Vertex	primaryVertex
edm::Handle< reco::BeamSpot >	recoBeamSpotHandle
edm::ParameterSet	regfactoryPSet
Regions	regions
TrajectorySeedCollection *	seedCollection
TrajectorySeedCollection *	seedCollectionOfSourceTracks
std::stringstream	ss
const reco::BeamSpot *	theBeamSpot
ClusterChecker	theClusterCheck
CombinedHitPairGeneratorForPhotonConversion *	theHitsGenerator
GlobalTrackingRegionProducerFromBeamSpot *	theRegionProducer
SeedForPhotonConversion1Leg *	theSeedCreator
bool	theSilentOnClusterCheck
edm::Handle < reco::TrackCollection >	trackCollectionH
edm::Handle < TrajTrackAssociationCollection >	trajTrackAssociations
reco::VertexCollection	vertexCollection
edm::Handle < reco::VertexCollection >	vertexHandle

Detailed Description

Definition at line 32 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

Member Typedef Documentation

typedef Regions::const_iterator PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::IR [private]

Definition at line 80 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

typedef std::vector<TrackingRegion* > PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::Regions [private]

Definition at line 79 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

Constructor & Destructor Documentation

PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::PhotonConversionTrajectorySeedProducerFromSingleLegAlgo ( const edm::ParameterSet & conf )

Definition at line 9 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References init().

  :_conf(conf),seedCollection(0),seedCollectionOfSourceTracks(0),
   hitsfactoryPSet(conf.getParameter<edm::ParameterSet>("OrderedHitsFactoryPSet")),   
   creatorPSet(conf.getParameter<edm::ParameterSet>("SeedCreatorPSet")),
   regfactoryPSet(conf.getParameter<edm::ParameterSet>("RegionFactoryPSet")),
   theClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")),
   theSilentOnClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet").getUntrackedParameter<bool>("silentClusterCheck",false)),
   _vtxMinDoF(conf.getParameter<double>("vtxMinDoF")),
   _maxDZSigmas(conf.getParameter<double>("maxDZSigmas")),
   _maxNumSelVtx(conf.getParameter<uint32_t>("maxNumSelVtx")),
   _applyTkVtxConstraint(conf.getParameter<bool>("applyTkVtxConstraint")),
   _countSeedTracks(0),
   _primaryVtxInputTag(conf.getParameter<edm::InputTag>("primaryVerticesTag")),
   _beamSpotInputTag(conf.getParameter<edm::InputTag>("beamSpotInputTag"))
{
   init();  
}

PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::~PhotonConversionTrajectorySeedProducerFromSingleLegAlgo ( ) [inline]

Definition at line 38 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

{};

Member Function Documentation

void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	setup
	)

Definition at line 45 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _beamSpotInputTag, _IdealHelixParameters, _primaryVtxInputTag, CombinedHitPairGeneratorForPhotonConversion::clearLayerCache(), event(), edm::EventSetup::get(), edm::HandleBase::isValid(), loopOnTracks(), magField, myEsetup, myEvent, po, PrintRecoObjects::print(), edm::ESHandle< T >::product(), recoBeamSpotHandle, regions, GlobalTrackingRegionProducerFromBeamSpot::regions(), seedCollection, seedCollectionOfSourceTracks, IdealHelixParameters::setMagnField(), HcalObjRepresent::setup(), ss, theClusterCheck, theHitsGenerator, theRegionProducer, theSilentOnClusterCheck, ClusterChecker::tooManyClusters(), and vertexHandle.

Referenced by PhotonConversionTrajectorySeedProducerFromSingleLeg::produce().

                                                           {

  myEsetup = &setup;
  myEvent = &event;

  if(seedCollection!=0)
    delete seedCollection;

  if(seedCollectionOfSourceTracks!=0)
    delete seedCollectionOfSourceTracks;

  seedCollection= new TrajectorySeedCollection();
  seedCollectionOfSourceTracks= new TrajectorySeedCollection();

  size_t clustsOrZero = theClusterCheck.tooManyClusters(event);
  if (clustsOrZero){
    if (!theSilentOnClusterCheck)
      edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
    return ;
  }


  edm::ESHandle<MagneticField> handleMagField;
  setup.get<IdealMagneticFieldRecord>().get(handleMagField);
  magField = handleMagField.product();

  _IdealHelixParameters.setMagnField(magField);


  event.getByLabel(_primaryVtxInputTag, vertexHandle);
  if (!vertexHandle.isValid() || vertexHandle->empty()){
      edm::LogError("PhotonConversionFinderFromTracks") << "Error! Can't get the product primary Vertex Collection "<< _primaryVtxInputTag <<  "\n";
      return;
  }

  event.getByLabel(_beamSpotInputTag,recoBeamSpotHandle);
  

  regions = theRegionProducer->regions(event,setup);
  
  //Do the analysis
  loopOnTracks();

 
#ifdef debugTSPFSLA 
  std::stringstream ss;
  ss.str("");
  ss << "\n++++++++++++++++++\n";
  ss << "seed collection size " << seedCollection->size();
  BOOST_FOREACH(TrajectorySeed tjS,*seedCollection){
    po.print(ss, tjS);
  }
  edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
  //-------------------------------------------------
#endif

   // clear memory
  theHitsGenerator->clearLayerCache();
  for (IR ir=regions.begin(), irEnd=regions.end(); ir < irEnd; ++ir) delete (*ir);

}

void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::clear ( void )

Definition at line 28 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References NULL, theHitsGenerator, theRegionProducer, and theSeedCreator.

Referenced by PhotonConversionTrajectorySeedProducerFromSingleLeg::endRun().

       {
  if(theHitsGenerator!=NULL)
    delete theHitsGenerator;
  if(theSeedCreator!=NULL)
    delete theSeedCreator;
  if(theRegionProducer!=NULL)
    delete theRegionProducer;
}

IdealHelixParameters* PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::getIdealHelixParameters ( ) [inline]

Definition at line 44 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

References _IdealHelixParameters.

{return &_IdealHelixParameters;}

TrajectorySeedCollection* PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::getTrajectorySeedCollection ( ) [inline]

Definition at line 45 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

References seedCollection.

Referenced by PhotonConversionTrajectorySeedProducerFromSingleLeg::produce().

{return seedCollection;}

TrajectorySeedCollection* PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::getTrajectorySeedCollectionOfSourceTracks ( ) [inline]

Definition at line 46 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h.

References seedCollectionOfSourceTracks.

Referenced by PhotonConversionTrajectorySeedProducerFromSingleLeg::produce().

{return seedCollectionOfSourceTracks;}

void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::init ( void )

Definition at line 38 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References creatorPSet, hitsfactoryPSet, regfactoryPSet, theHitsGenerator, theRegionProducer, and theSeedCreator.

Referenced by PhotonConversionTrajectorySeedProducerFromSingleLeg::beginRun(), and PhotonConversionTrajectorySeedProducerFromSingleLegAlgo().

      {
  theHitsGenerator  = new CombinedHitPairGeneratorForPhotonConversion(hitsfactoryPSet);
  theSeedCreator    = new SeedForPhotonConversion1Leg(creatorPSet);
  theRegionProducer = new GlobalTrackingRegionProducerFromBeamSpot(regfactoryPSet);
}

bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::inspectTrack	(	const reco::Track *	track,
		const TrackingRegion &	region,
		math::XYZPoint &	primaryVertexPoint
	)		`[private]`

Definition at line 298 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _IdealHelixParameters, abs, reco::TrackBase::charge(), f, IdealHelixParameters::GetMomentumAtTangentPoint(), IdealHelixParameters::GetTangentPoint(), myEsetup, myEvent, PV3DBase< T, PVType, FrameType >::perp(), ptmin, CombinedHitPairGeneratorForPhotonConversion::run(), seedCollection, IdealHelixParameters::setData(), OrderedSeedingHits::size(), ss, theHitsGenerator, theSeedCreator, reco::TrackBase::thetaError(), and SeedForPhotonConversion1Leg::trajectorySeed().

Referenced by loopOnPriVtx().

                                                                                                     {

  _IdealHelixParameters.setData(track,primaryVertexPoint);   
    
  if(_IdealHelixParameters.GetTangentPoint().r()==0){
    //this case means a null results on the _IdealHelixParameters side
    return false;
  }

  float rMin=3.; //cm
  if(_IdealHelixParameters.GetTangentPoint().rho()<rMin){
    //this case means a track that has the tangent point nearby the primary vertex
    // if the track is primary, this number tends to be the primary vertex itself
    //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
    //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
    return false;
  }

  float ptmin = 0.5;
  float originRBound = 3;
  float originZBound  = 3.;

  GlobalPoint originPos;
  originPos = GlobalPoint(_IdealHelixParameters.GetTangentPoint().x(),
                          _IdealHelixParameters.GetTangentPoint().y(),
                          _IdealHelixParameters.GetTangentPoint().z()
                          );
  float cotTheta;
  if( std::abs(_IdealHelixParameters.GetMomentumAtTangentPoint().rho()) > 1.e-4f ){
    cotTheta=_IdealHelixParameters.GetMomentumAtTangentPoint().z()/_IdealHelixParameters.GetMomentumAtTangentPoint().rho();
  }else{
    if(_IdealHelixParameters.GetMomentumAtTangentPoint().z()>0)
      cotTheta=99999.f; 
    else
      cotTheta=-99999.f; 
  }
  GlobalVector originBounds(originRBound,originRBound,originZBound);

  GlobalPoint pvtxPoint(primaryVertexPoint.x(),
                        primaryVertexPoint.y(),
                        primaryVertexPoint.z()
                        );
  ConversionRegion convRegion(originPos, pvtxPoint, cotTheta, track->thetaError(), -1*track->charge());

#ifdef debugTSPFSLA 
  ss << "\nConversion Point " << originPos << " " << originPos.perp() << "\n";
#endif

  const OrderedSeedingHits & hitss = theHitsGenerator->run(convRegion, region, *myEvent, *myEsetup);
  
  unsigned int nHitss =  hitss.size();

  if(nHitss==0)
    return false;

#ifdef debugTSPFSLA 
  ss << "\n nHitss " << nHitss << "\n";
#endif

  if (seedCollection->empty()) seedCollection->reserve(nHitss); // don't do multiple reserves in the case of multiple regions: it would make things even worse
                                                               // as it will cause N re-allocations instead of the normal log(N)/log(2)
  for (unsigned int iHits = 0; iHits < nHitss; ++iHits) { 

#ifdef debugTSPFSLA 
    ss << "\n iHits " << iHits << "\n";
#endif
    const SeedingHitSet & hits =  hitss[iHits];
    //if (!theComparitor || theComparitor->compatible( hits, es) ) {
    //try{
    theSeedCreator->trajectorySeed(*seedCollection,hits, originPos, originBounds, ptmin, *myEsetup,convRegion.cotTheta(),ss);
    //}catch(cms::Exception& er){
    //  edm::LogError("SeedingConversion") << " Problem in the Single Leg Seed creator " <<er.what()<<std::endl;
    //}catch(std::exception& er){
    //  edm::LogError("SeedingConversion") << " Problem in the Single Leg Seed creator " << er.what()<<std::endl;
    //}
  }
  return true;
}

void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnPriVtx	(	const reco::Track &	tk,
		const std::vector< reco::Vertex > &	selectedPriVtxCompatibleWithTrack
	)		`[private]`

Definition at line 205 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _countSeedTracks, inspectTrack(), reco::Vertex::position(), TrackingRegion::print(), regions, and ss.

Referenced by loopOnTracks().

                                                                                                 {

  bool foundAtLeastASeedCand=false;
  BOOST_FOREACH(const reco::Vertex vtx, selectedPriVtxCompatibleWithTrack){

    math::XYZPoint primaryVertexPoint=math::XYZPoint(vtx.position());
    
    for (IR ir=regions.begin(), irEnd=regions.end(); ir < irEnd; ++ir) {
      const TrackingRegion & region = **ir;

#ifdef debugTSPFSLA 
      ss << "[PrintRegion] " << region.print() << std::endl;
#endif
      
      //This if is just for the _countSeedTracks. otherwise 
      //inspectTrack(&tk,region, primaryVertexPoint);
      //would be enough

      if(
         inspectTrack(&tk,region, primaryVertexPoint)
         and
         !foundAtLeastASeedCand
         ){
        foundAtLeastASeedCand=true;
        _countSeedTracks++; 
      }

    }
  }
}

void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnTracks ( ) [private]

Definition at line 109 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _applyTkVtxConstraint, _conf, _countSeedTracks, edm::Event::getByLabel(), edm::ParameterSet::getParameter(), edm::HandleBase::isValid(), loopOnPriVtx(), myEvent, rejectTrack(), runTheMatrix::sel, selectPriVtxCompatibleWithTrack(), ss, trackCollectionH, and vertexHandle.

Referenced by analyze().

              {

  //--- Get Tracks
  myEvent->getByLabel(_conf.getParameter<edm::InputTag>("TrackRefitter"),trackCollectionH);

  if(trackCollectionH.isValid()==0){
    edm::LogError("MissingInput")<<" could not find track collecion:"<<_conf.getParameter<edm::InputTag>("TrackRefitter");
    return;
  }
  size_t idx=0, sel=0;
  _countSeedTracks=0;

  ss.str("");
  
  for( reco::TrackCollection::const_iterator tr = trackCollectionH->begin(); 
       tr != trackCollectionH->end(); tr++, idx++) {
    
    // #ifdef debugTSPFSLA 
    //     ss << "\nStuding track Nb " << idx;
    // #endif

    if(rejectTrack(*tr))  continue;
    std::vector<reco::Vertex> selectedPriVtxCompatibleWithTrack;  
    if(!_applyTkVtxConstraint){
      selectedPriVtxCompatibleWithTrack.push_back(*(vertexHandle->begin())); //Same approach as before
    }else{
      if(!selectPriVtxCompatibleWithTrack(*tr,selectedPriVtxCompatibleWithTrack)) continue;
    }

    sel++;
    loopOnPriVtx(*tr,selectedPriVtxCompatibleWithTrack);
  }
#ifdef debugTSPFSLA 
  edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
  edm::LogInfo("debugTrajSeedFromSingleLeg") << "Inspected " << sel << " tracks over " << idx << " tracks. \t # tracks providing at least one seed " << _countSeedTracks ;
#endif
}

bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::rejectTrack ( const reco::Track & track ) [private]

Definition at line 237 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _IdealHelixParameters, IdealHelixParameters::GetTangentPoint(), edm::HandleBase::isValid(), recoBeamSpotHandle, and IdealHelixParameters::setData().

Referenced by loopOnTracks().

                                   {
  
  math::XYZVector beamSpot;
  if(recoBeamSpotHandle.isValid()) {
    beamSpot =  math::XYZVector(recoBeamSpotHandle->position());

    _IdealHelixParameters.setData(&track,beamSpot);   
    if(_IdealHelixParameters.GetTangentPoint().r()==0){
      //this case means a null results on the _IdealHelixParameters side
      return true;
      }
      
    float rMin=2.; //cm
    if(_IdealHelixParameters.GetTangentPoint().rho()<rMin){
      //this case means a track that has the tangent point nearby the primary vertex
      // if the track is primary, this number tends to be the primary vertex itself
      //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
      //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
      return true;
    }
  }

  //-------------------------------------------------------
  /*
  float maxPt2=64.; //Cut on pt^2 Indeed doesn't do almost nothing
  if(track.momentum().Perp2() > maxPt2)
    return true;
  */
  //-------------------------------------------------------
  //Cut in the barrel eta region FIXME: to be extended to endcaps
  /*
  float maxEta=1.3; 
  if(fabs(track.eta()) > maxEta)
    return true;
  */
  //-------------------------------------------------------
  //Reject tracks that have a first valid hit in the pixel barrel/endcap layer/disk 1
  //assume that the hits are aligned along momentum
  /*
  const reco::HitPattern& p=track.hitPattern();
  for (int i=0; i<p.numberOfHits(); i++) {
    uint32_t hit = p.getHitPattern(i);
    // if the hit is valid and in pixel barrel, print out the layer
    if (! p.validHitFilter(hit) ) continue;
    if( (p.pixelBarrelHitFilter(hit) || p.pixelEndcapHitFilter(hit))
        &&
        p.getLayer(hit) == 1
        )
      return true;
    else
      break; //because the first valid hit is in a different layer
  }
  */
  //-------------------------------------------------------


  return false;
}

bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::selectPriVtxCompatibleWithTrack	(	const reco::Track &	tk,
		std::vector< reco::Vertex > &	selectedPriVtxCompatibleWithTrack
	)		`[private]`

Definition at line 148 of file PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc.

References _maxDZSigmas, _maxNumSelVtx, _vtxMinDoF, prof2calltree::count, reco::Vertex::covariance(), reco::TrackBase::covariance(), reco::TrackBase::dz(), reco::TrackBase::dzError(), i, lt_(), reco::Vertex::ndof(), reco::Vertex::position(), reco::TrackBase::pt(), reco::TrackBase::px(), reco::TrackBase::py(), reco::TrackBase::pz(), edm::second(), funct::sqr(), mathSSE::sqrt(), ss, vertexHandle, reco::TrackBase::vx(), reco::TrackBase::vy(), and reco::TrackBase::vz().

Referenced by loopOnTracks().

                                                                                                              {
  
  std::vector< std::pair< double, short> > idx;
  short count=-1;

  double cosPhi=tk.px()/tk.pt();
  double sinPhi=tk.py()/tk.pt();
  double sphi2=tk.covariance(2,2);
  double stheta2=tk.covariance(1,1);

  BOOST_FOREACH(const reco::Vertex& vtx,  *vertexHandle){
    count++;
    if(vtx.ndof()<= _vtxMinDoF) continue;

    double _dz= tk.dz(vtx.position());
    double _dzError=tk.dzError();

    double cotTheta=tk.pz()/tk.pt();
    double dx = vtx.position().x();
    double dy = vtx.position().y();
    double sx2=vtx.covariance(0,0);
    double sy2=vtx.covariance(1,1);

    double sxy2= sqr(cosPhi*cotTheta)*sx2+
      sqr(sinPhi*cotTheta)*sy2+
      sqr(cotTheta*(-dx*sinPhi+dy*cosPhi))*sphi2+
      sqr((1+cotTheta*cotTheta)*(dx*cosPhi+dy*sinPhi))*stheta2;
      
    _dzError=sqrt(_dzError*_dzError+vtx.covariance(2,2)+sxy2); //there is a missing component, related to the element (vtx.x*px+vtx.y*py)/pt * pz/pt. since the tk ref point is at the point of closest approach, this scalar product should be almost zero.

#ifdef debugTSPFSLA
    ss << " primary vtx " << vtx.position()  << " \tk vz " << tk.vz() << " vx " << tk.vx() << " vy " << tk.vy() << " pz/pt " << tk.pz()/tk.pt() << " \t dz " << _dz << " \t " << _dzError << " sxy2 "<< sxy2<< " \t dz/dzErr " << _dz/_dzError<< std::endl;
#endif

    if(fabs(_dz)/_dzError > _maxDZSigmas) continue; 
      
    idx.push_back(std::pair<double,short>(fabs(_dz),count));
  }
  if(idx.size()==0) {
#ifdef debugTSPFSLA
    ss << "no vertex selected " << std::endl; 
#endif
    return false;
}
  
  std::stable_sort(idx.begin(),idx.end(),lt_);
  for(size_t i=0;i<_maxNumSelVtx && i<idx.size();++i){
    selectedPriVtxCompatibleWithTrack.push_back((*vertexHandle)[idx[i].second]);
#ifdef debugTSPFSLA
    ss << "selected vtx dz " << idx[0].first << "  position" << idx[0].second << std::endl;
#endif
  }

  return true;
}