Inheritance diagram for TrackClusterSplitter:

Classes
struct	ClusterWithTracks

class	FindCluster

struct	TrackAndState

Public Member Functions
void	produce (edm::Event &iEvent, const edm::EventSetup &iSetup)

	TrackClusterSplitter (const edm::ParameterSet &iConfig)

	~TrackClusterSplitter ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Private Types
typedef std::pair< uint32_t, EncodedEventId >	SimHitIdpr

typedef boost::sub_range < std::vector < SiPixelClusterWithTracks > >	SiPixelClustersWithTracks

typedef ClusterWithTracks < SiPixelCluster >	SiPixelClusterWithTracks

typedef boost::sub_range < std::vector < SiStripClusterWithTracks > >	SiStripClustersWithTracks

typedef ClusterWithTracks < SiStripCluster >	SiStripClusterWithTracks

Private Member Functions
template<>
const SiPixelCluster *	getCluster (const TrackingRecHit *hit)

template<>
const SiStripCluster *	getCluster (const TrackingRecHit *hit)

template<>
const SiStripCluster *	getCluster (const TrackingRecHit *hit)

template<>
const SiPixelCluster *	getCluster (const TrackingRecHit *hit)

template<typename Cluster >
void	markClusters (std::map< uint32_t, boost::sub_range< std::vector< ClusterWithTracks< Cluster > > > > &map, const TrackingRecHit hit, const reco::Track track, const TrajectoryStateOnSurface &tsos) const

template<typename Cluster >
void	splitCluster (const ClusterWithTracks< Cluster > &cluster, const GlobalVector &dir, typename edmNew::DetSetVector< Cluster >::FastFiller &output, DetId detId) const

template<>
void	splitCluster (const SiPixelClusterWithTracks &cluster, const GlobalVector &dir, edmNew::DetSetVector< SiPixelCluster >::FastFiller &output, DetId detId) const

template<>
void	splitCluster (const SiStripClusterWithTracks &cluster, const GlobalVector &dir, edmNew::DetSetVector< SiStripCluster >::FastFiller &output, DetId detId) const

template<>
void	splitCluster (const SiStripClusterWithTracks &c, const GlobalVector &vtx, edmNew::DetSetVector< SiStripCluster >::FastFiller &output, DetId detId) const

template<>
void	splitCluster (const SiPixelClusterWithTracks &c, const GlobalVector &vtx, edmNew::DetSetVector< SiPixelCluster >::FastFiller &output, DetId detId) const

template<typename Cluster >
std::auto_ptr < edmNew::DetSetVector < Cluster > >	splitClusters (const std::map< uint32_t, boost::sub_range< std::vector< ClusterWithTracks< Cluster > > > > &input, const reco::Vertex &vtx) const

Static Private Member Functions
template<typename C >
static const C *	equalClusters (const C &c1, const C &c2)

template<typename C >
static const C *	getCluster (const TrackingRecHit *hit)

Private Attributes
std::vector < SiPixelClusterWithTracks >	allSiPixelClusters
	working data More...

std::vector < SiStripClusterWithTracks >	allSiStripClusters

edm::ESHandle < GlobalTrackingGeometry >	geometry_

TrackerHitAssociator *	hitAssociator

edm::ESHandle< MagneticField >	magfield_

edm::InputTag	pixelClusters_

edm::Handle< edm::DetSetVector < PixelDigiSimLink > >	pixeldigisimlink

edm::ESHandle< Propagator >	propagator_

std::string	propagatorName_

bool	simSplitPixel_

bool	simSplitStrip_

std::map< uint32_t, SiPixelClustersWithTracks >	siPixelDetsWithClusters

std::map< uint32_t, SiStripClustersWithTracks >	siStripDetsWithClusters

SiStripTemplate	strip_templ_

edm::InputTag	stripClusters_

edm::Handle< edm::DetSetVector < StripDigiSimLink > >	stripdigisimlink

SiPixelTemplate2D	templ2D_

SiPixelTemplate	templ_

bool	tmpSplitPixel_

bool	tmpSplitStrip_

edm::InputTag	trajectories_

bool	useStraightTracks_

bool	useTrajectories_

edm::InputTag	vertices_

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

Detailed Description

Definition at line 64 of file TrackClusterSplitter.cc.

Member Typedef Documentation

typedef std::pair<uint32_t, EncodedEventId> TrackClusterSplitter::SimHitIdpr

private

Definition at line 149 of file TrackClusterSplitter.cc.

typedef boost::sub_range<std::vector<SiPixelClusterWithTracks> > TrackClusterSplitter::SiPixelClustersWithTracks

private

Definition at line 144 of file TrackClusterSplitter.cc.

typedef ClusterWithTracks<SiPixelCluster> TrackClusterSplitter::SiPixelClusterWithTracks

private

Definition at line 139 of file TrackClusterSplitter.cc.

typedef boost::sub_range<std::vector<SiStripClusterWithTracks> > TrackClusterSplitter::SiStripClustersWithTracks

private

Definition at line 145 of file TrackClusterSplitter.cc.

typedef ClusterWithTracks<SiStripCluster> TrackClusterSplitter::SiStripClusterWithTracks

private

Definition at line 140 of file TrackClusterSplitter.cc.

Constructor & Destructor Documentation

TrackClusterSplitter::TrackClusterSplitter ( const edm::ParameterSet & iConfig )

Definition at line 231 of file TrackClusterSplitter.cc.

References edm::ParameterSet::getParameter(), propagatorName_, SiPixelTemplate2D::pushfile(), SiStripTemplate::pushfile(), SiPixelTemplate::pushfile(), simSplitPixel_, simSplitStrip_, strip_templ_, templ2D_, templ_, tmpSplitPixel_, tmpSplitStrip_, useStraightTracks_, and useTrajectories_.

                                                                         :
   stripClusters_(iConfig.getParameter<edm::InputTag>("stripClusters")),
   pixelClusters_(iConfig.getParameter<edm::InputTag>("pixelClusters")),
   useTrajectories_(iConfig.getParameter<bool>("useTrajectories")),
   trajectories_(iConfig.getParameter<edm::InputTag>(useTrajectories_ ? "trajTrackAssociations" : "tracks")),
   vertices_(iConfig.getParameter<edm::InputTag>("vertices"))
 {
   if ( !useTrajectories_ ) 
     propagatorName_ = iConfig.getParameter<std::string>("propagator");
   
   produces< edmNew::DetSetVector<SiPixelCluster> >();
 
   produces< edmNew::DetSetVector<SiStripCluster> >();
   
   simSplitPixel_ = (iConfig.getParameter<bool>("simSplitPixel"));
   simSplitStrip_ = (iConfig.getParameter<bool>("simSplitStrip"));
   tmpSplitPixel_ = (iConfig.getParameter<bool>("tmpSplitPixel")); // not so nice... you don't want two bool but some switch
   tmpSplitStrip_ = (iConfig.getParameter<bool>("tmpSplitStrip"));
 
   useStraightTracks_ = (iConfig.getParameter<bool>("useStraightTracks"));
 
   /*
     cout << "TrackClusterSplitter : " << endl;
     cout << endl << endl << endl;
     cout << "(int)simSplitPixel_   = " << (int)simSplitPixel_  << endl;
     cout << "(int)simSplitStrip_   = " << (int)simSplitStrip_  << endl;
     cout << "(int)tmpSplitPixel_   = " << (int)tmpSplitPixel_  << endl;
     cout << "(int)tmpSplitStrip_   = " << (int)tmpSplitStrip_  << endl;
     cout << "stripClusters_        = " << stripClusters_        << endl;
     cout << "pixelClusters_        = " << pixelClusters_        << endl;
     cout << "(int)useTrajectories_ = " << (int)useTrajectories_ << endl;
     cout << "trajectories_         = " << trajectories_         << endl;
     cout << "propagatorName_       = " << propagatorName_       << endl;
     cout << "vertices_             = " << vertices_             << endl;
     cout << "useStraightTracks_    = " << useStraightTracks_    << endl;
     cout << endl << endl << endl;
   */
 
   // Load template; 40 for barrel and 41 for endcaps
   templ_.pushfile( 40 );
   templ_.pushfile( 41 );
   templ2D_.pushfile( 40 );
   templ2D_.pushfile( 41 );
 
   // Load strip templates
   strip_templ_.pushfile( 11 );
   strip_templ_.pushfile( 12 );
   strip_templ_.pushfile( 13 );
   strip_templ_.pushfile( 14 );
   strip_templ_.pushfile( 15 );
   strip_templ_.pushfile( 16 );
 
 }

TrackClusterSplitter::~TrackClusterSplitter ( )

Definition at line 314 of file TrackClusterSplitter.cc.

315 {

316 }

Member Function Documentation

template<typename C >

static const C* TrackClusterSplitter::equalClusters	(	const C &	c1,
		const C &	c2
	)

inlinestaticprivate

Definition at line 156 of file TrackClusterSplitter.cc.

   { 
     return nullptr; 
   }

template<typename C >

static const C* TrackClusterSplitter::getCluster ( const TrackingRecHit * hit )

staticprivate

template<>

const SiPixelCluster* TrackClusterSplitter::getCluster ( const TrackingRecHit * hit )

private

template<>

const SiStripCluster* TrackClusterSplitter::getCluster ( const TrackingRecHit * hit )

private

template<>

const SiStripCluster* TrackClusterSplitter::getCluster ( const TrackingRecHit * hit )

private

Definition at line 288 of file TrackClusterSplitter.cc.

References edm::hlt::Exception.

 {
   if ( typeid(*hit) == typeid(SiStripRecHit2D) ) 
     {
       return (static_cast<const SiStripRecHit2D &>(*hit)).cluster().get();
     } 
   else if ( typeid(*hit) == typeid(SiStripRecHit1D) ) 
     {
       return (static_cast<const SiStripRecHit1D &>(*hit)).cluster().get();
     } 
   else 
     throw cms::Exception("Unsupported") << "Detid of type " << typeid(*hit).name() << " not supported.\n";
 }

template<>

const SiPixelCluster* TrackClusterSplitter::getCluster ( const TrackingRecHit * hit )

private

Definition at line 304 of file TrackClusterSplitter.cc.

References edm::hlt::Exception.

 {
   if ( typeid(*hit) == typeid(SiPixelRecHit) ) 
     {
       return (static_cast<const SiPixelRecHit&>(*hit)).cluster().get();
     } 
   else 
     throw cms::Exception("Unsupported") << "Detid of type " << typeid(*hit).name() << " not supported.\n";
 }

template<typename Cluster >

void TrackClusterSplitter::markClusters	(	std::map< uint32_t, boost::sub_range< std::vector< ClusterWithTracks< Cluster > > > > &	map,
		const TrackingRecHit *	hit,
		const reco::Track *	track,
		const TrajectoryStateOnSurface &	tsos
	)		const

private

Definition at line 518 of file TrackClusterSplitter.cc.

References TrackingRecHit::geographicalId(), Association::map, match(), mergeVDriftHistosByStation::name, and DetId::rawId().

 {
   boost::sub_range<std::vector<ClusterWithTracks<Cluster> > >& range = map[hit->geographicalId().rawId()];
   
   typedef typename std::vector<ClusterWithTracks<Cluster> >::iterator IT;
   IT match = std::find_if(range.begin(), range.end(), FindCluster<Cluster>(hit));
   
   if ( match != range.end() ) 
     {
       match->tracks.push_back( TrackAndState(track,tsos) );
     } 
   else 
     {
       edm::LogWarning("ClusterNotFound") << "Cluster of type " << typeid(Cluster).name() << " on detid " 
                      << hit->geographicalId().rawId() << " from hit of type " << typeid(*hit).name();
     }
 }

void TrackClusterSplitter::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

virtual

Implements edm::EDProducer.

Definition at line 319 of file TrackClusterSplitter.cc.

References allSiPixelClusters, allSiStripClusters, edmNew::DetSet< T >::begin(), edmNew::DetSet< T >::detId(), end, edmNew::DetSet< T >::end(), edm::hlt::Exception, first, TrackingRecHit::geographicalId(), geometry_, edm::EventSetup::get(), edm::Event::getByLabel(), hitAssociator, trajectoryStateTransform::innerFreeState(), TrackingRecHit::isValid(), magfield_, Trajectory::measurements(), pixelClusters_, pixeldigisimlink, propagator_, edm::Event::put(), DetId::rawId(), reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), reco::Track::recHitsSize(), simSplitPixel_, simSplitStrip_, siPixelDetsWithClusters, siStripDetsWithClusters, splitClusters(), errorMatrix2Lands_multiChannel::start, stripClusters_, stripdigisimlink, DetId::subdetId(), GeomDet::surface(), testEve_cfg::tracks, trajectories_, useTrajectories_, and vertices_.

 {
   using namespace edm;
 
   iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);  
   
   if ( !useTrajectories_ ) 
     {
       iSetup.get<IdealMagneticFieldRecord>().get( magfield_ );
       iSetup.get<TrackingComponentsRecord>().get( "AnalyticalPropagator", propagator_ );
     }
  
   Handle<edmNew::DetSetVector<SiPixelCluster> > inputPixelClusters;
   Handle<edmNew::DetSetVector<SiStripCluster> > inputStripClusters;
   
   iEvent.getByLabel(pixelClusters_, inputPixelClusters);
   
   iEvent.getByLabel(stripClusters_, inputStripClusters);
   
   if(simSplitStrip_)
     hitAssociator = new TrackerHitAssociator(iEvent);
 
     
   allSiPixelClusters.clear(); siPixelDetsWithClusters.clear();
   allSiStripClusters.clear(); siStripDetsWithClusters.clear();
 
   
   allSiPixelClusters.reserve(inputPixelClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators
   allSiStripClusters.reserve(inputStripClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators
 
   
   // fill in the list of all tracks
   foreach(const edmNew::DetSet<SiPixelCluster> &ds, *inputPixelClusters) 
     {
       std::vector<SiPixelClusterWithTracks>::iterator start = allSiPixelClusters.end();
       allSiPixelClusters.insert(start, ds.begin(), ds.end());
     
       std::vector<SiPixelClusterWithTracks>::iterator end   = allSiPixelClusters.end();
       siPixelDetsWithClusters[ds.detId()] = SiPixelClustersWithTracks(start,end);
     }
  
   foreach(const edmNew::DetSet<SiStripCluster> &ds, *inputStripClusters) 
     {
       std::vector<SiStripClusterWithTracks>::iterator start = allSiStripClusters.end();
       allSiStripClusters.insert(start, ds.begin(), ds.end());
      
       std::vector<SiStripClusterWithTracks>::iterator end   = allSiStripClusters.end();
       siStripDetsWithClusters[ds.detId()] = SiStripClustersWithTracks(start,end);
     }
   
   if ( useTrajectories_ ) 
     { 
       // Here use the fully reconstructed tracks to get the track angle
 
       Handle<TrajTrackAssociationCollection> trajectories; 
       iEvent.getByLabel(trajectories_, trajectories);
       for ( TrajTrackAssociationCollection::const_iterator it = trajectories->begin(), 
           ed = trajectories->end(); it != ed; ++it ) 
     { 
       const Trajectory  & traj =  *it->key;
       const reco::Track * tk   = &*it->val;
       
       if ( traj.measurements().size() != tk->recHitsSize() ) 
         throw cms::Exception("Aargh") << "Sizes don't match: traj " << traj.measurements().size()  
                       << ", tk " << tk->recHitsSize() << "\n";
       
       trackingRecHit_iterator it_hit = tk->recHitsBegin(), ed_hit = tk->recHitsEnd();
       
       const Trajectory::DataContainer & tms = traj.measurements();
       
       size_t i_hit = 0, last_hit = tms.size()-1;
       
       bool first = true, reversed = false;
       
       for (; it_hit != ed_hit; ++it_hit, ++i_hit) 
         {
           // ignore hits with no detid
 
           if ((*it_hit)->geographicalId().rawId() == 0)
         {
           //cout << "It should never happen that a trackingRecHit has no detector ID !!!!!!!!!!!!!!!!! " << endl;
           continue;
         }
 
           // if it's the first hit, check the ordering of track vs trajectory
           if (first) 
         {
           
           if ((*it_hit)->geographicalId() == tms[i_hit].recHit()->hit()->geographicalId()) 
             {
               reversed = false;
                     } 
           else if ((*it_hit)->geographicalId() == tms[last_hit-i_hit].recHit()->hit()->geographicalId()) 
             {
               reversed = true;
                     } 
           else 
             {
               throw cms::Exception("Aargh") << "DetIDs don't match either way :-( \n";
                     }
                 }  
 
           const TrackingRecHit *hit = it_hit->get();
           if ( hit == 0 || !hit->isValid() )
         continue;
           
           int subdet = hit->geographicalId().subdetId();
           
           if (subdet >= 3) 
         { // strip
           markClusters<SiStripCluster>(siStripDetsWithClusters, hit, tk, tms[reversed ? last_hit-i_hit : i_hit].updatedState());
                 } 
           else if (subdet >= 1) 
         { // pixel
           markClusters<SiPixelCluster>(siPixelDetsWithClusters, hit, tk, tms[reversed ? last_hit-i_hit : i_hit].updatedState());
                 } 
           else 
         {
           edm::LogWarning("HitNotFound") << "Hit of type " << typeid(*hit).name() << ",  detid " 
                          << hit->geographicalId().rawId() << ", subdet " << subdet;
                 }
             }
         }
     } 
   else 
     {
       // Here use the pixel tracks to get the track angles
 
       Handle<std::vector<reco::Track> > tracks; 
       iEvent.getByLabel(trajectories_, tracks);
       //TrajectoryStateTransform transform;
       foreach (const reco::Track &track, *tracks) 
     {
       FreeTrajectoryState atVtx   =  trajectoryStateTransform::innerFreeState(track, &*magfield_);
       trackingRecHit_iterator it_hit = track.recHitsBegin(), ed_hit = track.recHitsEnd();
       for (; it_hit != ed_hit; ++it_hit) 
         {
           const TrackingRecHit *hit = it_hit->get();
           if ( hit == 0 || !hit->isValid() ) 
         continue;
           
           int subdet = hit->geographicalId().subdetId();
 
           if ( subdet == 0 ) 
         continue;
           
           const GeomDet *det = geometry_->idToDet( hit->geographicalId() );
           
           if ( det == 0 ) 
         {
           edm::LogError("MissingDetId") << "DetIDs " << (int)(hit->geographicalId()) << " is not in geometry.\n";
           continue;
                 }
           
           TrajectoryStateOnSurface prop = propagator_->propagate(atVtx, det->surface());
           if ( subdet >= 3 ) 
         { // strip
           markClusters<SiStripCluster>(siStripDetsWithClusters, hit, &track, prop);
                 } 
           else if (subdet >= 1) 
         { // pixel
           markClusters<SiPixelCluster>(siPixelDetsWithClusters, hit, &track, prop);
                 } 
           else 
         {
           edm::LogWarning("HitNotFound") << "Hit of type " << typeid(*hit).name() << ",  detid " 
                          << hit->geographicalId().rawId() << ", subdet " << subdet;
         }
             }
         }
     }
 
   Handle<std::vector<reco::Vertex> > vertices; 
   iEvent.getByLabel(vertices_, vertices);
   
   // Needed in case of simsplit
   if ( simSplitPixel_ )
     iEvent.getByLabel("simSiPixelDigis", pixeldigisimlink);
     
   // Needed in case of strip simsplit
   if ( simSplitStrip_ )
     iEvent.getByLabel("simSiStripDigis", stripdigisimlink);
 
   // gavril : to do: choose the best vertex here instead of just choosing the first one ? 
   std::auto_ptr<edmNew::DetSetVector<SiPixelCluster> > newPixelClusters( splitClusters( siPixelDetsWithClusters, vertices->front() ) );
   std::auto_ptr<edmNew::DetSetVector<SiStripCluster> > newStripClusters( splitClusters( siStripDetsWithClusters, vertices->front() ) );
   
   if ( simSplitStrip_ )
     delete hitAssociator;
 
   iEvent.put(newPixelClusters);
   iEvent.put(newStripClusters);
     
   allSiPixelClusters.clear(); siPixelDetsWithClusters.clear();
   allSiStripClusters.clear(); siStripDetsWithClusters.clear();
 
 }

template<typename Cluster >

void TrackClusterSplitter::splitCluster	(	const ClusterWithTracks< Cluster > &	cluster,
		const GlobalVector &	dir,
		typename edmNew::DetSetVector< Cluster >::FastFiller &	output,
		DetId	detId
	)		const

private

Definition at line 575 of file TrackClusterSplitter.cc.

References edm::hlt::Exception.

Referenced by splitClusters().

 {
   //cout << "Should never be here: TrackClusterSplitter, TrackClusterSplitter::splitCluster(...)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
   throw cms::Exception("LogicError", "This should not be called");
 }

template<>

void TrackClusterSplitter::splitCluster	(	const SiPixelClusterWithTracks &	cluster,
		const GlobalVector &	dir,
		edmNew::DetSetVector< SiPixelCluster >::FastFiller &	output,
		DetId	detId
	)		const

private

template<>

void TrackClusterSplitter::splitCluster	(	const SiStripClusterWithTracks &	cluster,
		const GlobalVector &	dir,
		edmNew::DetSetVector< SiStripCluster >::FastFiller &	output,
		DetId	detId
	)		const

private

template<>

void TrackClusterSplitter::splitCluster	(	const SiStripClusterWithTracks &	c,
		const GlobalVector &	vtx,
		edmNew::DetSetVector< SiStripCluster >::FastFiller &	output,
		DetId	detId
	)		const

private

Definition at line 585 of file TrackClusterSplitter.cc.

References SiStripCluster::amplitudes(), TrackerHitAssociator::associateSimpleRecHitCluster(), SiStripCluster::barycenter(), begin, BXSIZE, trackerHits::c, GetRecoTauVFromDQM_MC_cff::cl2, edm::DetSet< T >::data, end, edm::hlt::Exception, first, SiStripCluster::firstStrip(), geometry_, hitAssociator, i, SiStripTemplate::interpolate(), j, prof2calltree::last, TrajectoryStateOnSurface::localParameters(), StripTopology::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), magfield_, max(), min, SiStripDetId::moduleGeometry(), LocalTrajectoryParameters::momentum(), convertSQLitetoXML_cfg::output, GloballyPositioned< T >::position(), reco::TrackBase::pt(), SiStripTemplate::qbin(), SiStripCluster::setSplitClusterError(), simSplitStrip_, findQualityFiles::size, StripGeomDetUnit::specificTopology(), mathSSE::sqrt(), strip_templ_, stripdigisimlink, SiStripTemplateSplit::StripTempSplit(), GeomDet::surface(), SiStripTemplate::sxtemp(), tmpSplitStrip_, Surface::toGlobal(), toLocal(), useStraightTracks_, PV3DBase< T, PVType, FrameType >::x(), SiStripTemplate::xsize(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

 { 
   if ( simSplitStrip_ )
     {      
       bool cluster_was_successfully_split = false;
             
       const SiStripCluster* clust = static_cast<const SiStripCluster*>(c.cluster);
       
       std::vector<SimHitIdpr> associatedIdpr;
       
       hitAssociator->associateSimpleRecHitCluster(clust, detId, associatedIdpr);
 
       size_t splittableClusterSize = 0;
       splittableClusterSize = associatedIdpr.size();
       std::vector<uint8_t> amp = clust->amplitudes();
       int clusiz = amp.size();
       associatedIdpr.clear();
 
       SiStripDetId ssdid( detId );
 
       // gavril : sim splitting can be applied to the forward detectors as well...
 
       if ( ( splittableClusterSize > 1 && amp.size() > 2 ) && 
        ( (int)ssdid.moduleGeometry() == 1 || 
          (int)ssdid.moduleGeometry() == 2 || 
          (int)ssdid.moduleGeometry() == 3 || 
          (int)ssdid.moduleGeometry() == 4 ) )
     {
       
       edm::DetSetVector<StripDigiSimLink>::const_iterator isearch = stripdigisimlink->find(detId);
       
       int first  = clust->firstStrip();
       int last   = first + clusiz;
       uint16_t rawAmpl = 0, currentAmpl = 0;
       
       std::vector<uint16_t> tmp1, tmp2;
       
       std::vector<int> firstStrip;
       std::vector<bool> trackInStrip;
       std::vector<unsigned int> trackID;
       std::vector<float> trackFraction;
       std::vector< std::vector<uint16_t> > trackAmp;
       unsigned int currentChannel( 9999 );
       unsigned int thisTrackID = 0;
 
       if ( isearch != stripdigisimlink->end() ) 
         {
           edm::DetSet<StripDigiSimLink> link_detset = (*isearch);
           
           for ( edm::DetSet<StripDigiSimLink>::const_iterator linkiter = link_detset.data.begin();
           linkiter != link_detset.data.end(); linkiter++)
         {
           if ( (int)(linkiter->channel()) >= first  && (int)(linkiter->channel()) < last  )
             {
               int stripIdx = (int)linkiter->channel()-first;
               rawAmpl = (uint16_t)(amp[stripIdx]);  
               
               // DigiSimLinks are ordered first by channel; there can be > 1 track, and > 1 simHit for each track
               
               if ( linkiter->channel() != currentChannel ) 
             {
               // New strip; store amplitudes for the previous one
               uint16_t thisAmpl;
               
               for (size_t i=0; i<trackID.size(); ++i)  
                 {
                   if ( trackInStrip[i] ) 
                 {
                   if ( ( thisAmpl=currentAmpl ) < 254 )
                     thisAmpl = min( uint16_t(253), max(uint16_t(0), (uint16_t)(currentAmpl*trackFraction[i]+0.5)) );
                   trackAmp[i].push_back( thisAmpl );
                 }
                   
                   trackFraction[i] = 0;
                   trackInStrip[i] = false;
                 }
               
               currentChannel = linkiter->channel();
               currentAmpl = rawAmpl;
             }
               
               // Now deal with this new DigiSimLink
               thisTrackID = linkiter->SimTrackId();
               
               // Have we seen this track yet?
               bool newTrack = true;
               int thisTrackIdx = 9999;
               
               for (size_t i=0; i<trackID.size(); ++i) 
             {
               if ( trackID[i] == thisTrackID ) 
                 {
                   thisTrackIdx = i;
                   newTrack = false;
                 }
             }
              
               if ( newTrack ) 
             {
               trackInStrip.push_back(false);  // We'll set it true below
               trackID.push_back(thisTrackID);
               firstStrip.push_back(currentChannel);
               std::vector<uint16_t> ampTmp;
               trackAmp.push_back(ampTmp);
               trackFraction.push_back(0);
               thisTrackIdx = trackID.size()-1;
             }
               
               trackInStrip[thisTrackIdx] = true;
               trackFraction[thisTrackIdx] += linkiter->fraction();
               currentAmpl = rawAmpl;
                
             }
           
         }// end of loop over DigiSimLinks
           
           // we want to continue here!!!!
           
           std::vector<SiStripCluster> newCluster;
           // Fill amplitudes for the last strip and create a cluster for each track
           uint16_t thisAmpl;
           
           for (size_t i=0; i < trackID.size(); ++i) 
         { 
           if ( trackInStrip[i] ) 
             {
               if ( ( thisAmpl=rawAmpl ) < 254 ) 
             thisAmpl = min(uint16_t(253), max(uint16_t(0), (uint16_t)(rawAmpl*trackFraction[i]+0.5)));
               
               if ( thisAmpl > 0 )
             trackAmp[i].push_back( thisAmpl );
               //else
               //cout << "thisAmpl = " << (int)thisAmpl << endl;
             }
           
           newCluster.push_back( SiStripCluster( detId, 
                             firstStrip[i],
                             trackAmp[i].begin(), 
                             trackAmp[i].end() ) );
           
         }
           
           
           for ( size_t i=0; i<newCluster.size(); ++i )
         {
           
           // gavril : This is never used. Will use it below 
           float clusterAmp = 0.0;
           for (size_t j=0; j<trackAmp[i].size(); ++j ) 
             clusterAmp += (float)(trackAmp[i])[j];
           
           if ( clusterAmp > 0.0 && firstStrip[i] != 9999 && trackAmp[i].size() > 0  ) 
             { 
               // gavril :  I think this should work
               output.push_back( newCluster[i] );
               
               //cout << endl << endl << endl; 
               //cout << "(int)(newCluster[i].amplitudes().size()) = " << (int)(newCluster[i].amplitudes().size()) << endl;
               //for ( int j=0; j<(int)(newCluster[i].amplitudes().size()); ++j )
               //cout << "(newCluster[i].amplitudes())[j] = " << (int)(newCluster[i].amplitudes())[j] << endl;
               
               cluster_was_successfully_split = true;
             } 
           else 
             {
               //std::cout << "\t\t Rejecting new cluster" << std::endl;
               
               // gavril : I think this pointer should be deleted above
               //delete newCluster[i];
             }
         }
                 
         } // if ( isearch != stripdigisimlink->end() ) ... 
       else 
         {
           // Do nothing...
         }
     }
       
       
       if ( !cluster_was_successfully_split )
     output.push_back( *c.cluster ); 
       
     } // end of   if ( strip_simSplit_ )...
   
   else if ( tmpSplitStrip_ )
    {
       bool cluster_was_successfully_split = false;
       
       const SiStripCluster* theStripCluster = static_cast<const SiStripCluster*>(c.cluster);
       
       if ( theStripCluster )
     {
       //SiStripDetId ssdid( theStripCluster->geographicalId() );
       SiStripDetId ssdid( detId.rawId() );
 
       // Do not attempt to split clusters of size less than or equal to one.
       // Only split clusters in IB1, IB2, OB1, OB2 (TIB and TOB).
 
       if ( (int)theStripCluster->amplitudes().size() <= 1 || 
            ( (int)ssdid.moduleGeometry() != 1 && 
          (int)ssdid.moduleGeometry() != 2 && 
          (int)ssdid.moduleGeometry() != 3 && 
          (int)ssdid.moduleGeometry() != 4 ) )
         {    
           // Do nothing.
           //cout << endl;
           //cout << "Will NOT attempt to split this clusters: " << endl;
           //cout << "(int)theStripCluster->amplitudes().size() = " << (int)theStripCluster->amplitudes().size() << endl;
           //cout << "(int)ssdid.moduleGeometry() = " << (int)ssdid.moduleGeometry() << endl;
           
         }
       else // if ( (int)theStripCluster->amplitudes().size() <= 1 )
         {
          
           //cout << endl;
           //cout << "Will attempt to split this clusters: " << endl;
 
           int ID = -99999;
                  
           int is_stereo = (int)( ssdid.stereo() ); 
           
           if      ( ssdid.moduleGeometry() == 1 ) // IB1 
         {
           if ( !is_stereo ) 
             ID = 11;
           else
             ID = 12;
         }
           else if ( ssdid.moduleGeometry() == 2 ) // IB2
         {
           ID = 13;
         }
           else if ( ssdid.moduleGeometry() == 3 ) // OB1
         {
           ID = 16; 
         }
           else if ( ssdid.moduleGeometry() == 4 ) // OB2
         {
           if ( !is_stereo )
             ID = 14;
           else
             ID = 15;
         }
           else 
         {
           throw cms::Exception("TrackClusterSplitter::splitCluster") 
             << "\nERROR: Wrong strip teplate ID. Should only use templates for IB1, IB2, OB1 and OB2 !!!" << "\n\n";
         }
           
 
 
 
           // Begin: determine incident angles ============================================================
           
           float cotalpha_ = -99999.9;
           float cotbeta_  = -99999.9;
           
           // First, determine track angles from straight track approximation
           
           // Find crude cluster center
           // gavril : This is in local coordinates ? 
           float xcenter = theStripCluster->barycenter();
          
           const GeomDetUnit* theDet = geometry_->idToDetUnit( detId );
           const StripGeomDetUnit* stripDet = dynamic_cast<const StripGeomDetUnit*>( theDet );
               
           if ( !stripDet )
         {
           throw cms::Exception("TrackClusterSplitter : ") 
             << "\nERROR: Wrong stripDet !!! " << "\n\n";
         }
 
 
           const StripTopology* theTopol = &( stripDet->specificTopology() );
 
           // Transform from measurement to local coordinates (in cm)
           // gavril: may have to differently if kicks/bows are introduced. However, at this point there are no tracks...:
           // LocalPoint lp = theTopol->localPosition( xcenter, /*const Topology::LocalTrackPred & */ trkPred );
 
           // gavril : What is lp.y() for strips ? It is zero, but is that the strip center or one of the ends ? 
           LocalPoint lp = theTopol->localPosition( xcenter );
           
           // Transform from local to global coordinates
           GlobalPoint gp0 = theDet->surface().toGlobal( lp );
           
            // Make a vector pointing from the PV to the cluster center 
           GlobalPoint gp(gp0.x()-vtx.x(), gp0.y()-vtx.y(), gp0.z()-vtx.z() );
           
           // Make gp a unit vector, gv, pointing from the PV to the cluster center
           float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
           float gpx = gp.x()/gp_mod;
           float gpy = gp.y()/gp_mod;
           float gpz = gp.z()/gp_mod;
           GlobalVector gv(gpx, gpy, gpz);
 
           // Make unit vectors in local coordinates and then transform them in global coordinates
           const Local3DVector lvx(1.0, 0.0, 0.0);
           GlobalVector gvx = theDet->surface().toGlobal( lvx );
           const Local3DVector lvy(0.0, 1.0, 0.0);
           GlobalVector gvy = theDet->surface().toGlobal( lvy );
           const Local3DVector lvz(0.0, 0.0, 1.0);
           GlobalVector gvz = theDet->surface().toGlobal( lvz );
           
           // Calculate angles alpha and beta
           float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
           float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
           float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
           
           // gavril : check beta !!!!!!!!!!!!
           //float alpha_ = atan2( gv_dot_gvz, gv_dot_gvx );
           //float beta_  = atan2( gv_dot_gvz, gv_dot_gvy );
           
           cotalpha_ = gv_dot_gvx / gv_dot_gvz;
           cotbeta_  = gv_dot_gvy / gv_dot_gvz;
            
           // Attempt to get a better angle from tracks (either pixel tracks or full tracks)
           if ( !useStraightTracks_ )
         {
           //cout << "TrackClusterSplitter.cc : " << endl;
           //cout << "Should not be here for now !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
 
           // Use either pixel tracks (useTrajectories_ = False) or fully reconstructed tracks (useTrajectories_ = True)
           // When pixel/full tracks are not associated with the current cluster, will use angles from straight tracks
           
           // These are the tracks associated with this cluster
           std::vector<TrackAndState> vec_tracks_states = c.tracks;
 
           if ( (int)vec_tracks_states.size() > 0 )
             {
               //cout << "There is at least one track associated with this cluster. Pick the one with largest Pt." << endl;
               //cout << "(int)vec_tracks_states.size() = " << (int)vec_tracks_states.size() << endl;
               
               int index_max_pt = -99999;  // index of the track with the highest Pt
               float max_pt = -99999.9;
               
               for (int i=0; i<(int)vec_tracks_states.size(); ++i )
             {
               const reco::Track* one_track = vec_tracks_states[i].track;
               
               if ( one_track->pt() > max_pt )
                 {
                   index_max_pt = i;
                   max_pt = one_track->pt();
                 }  
             }
               
               // Pick the tsos from the track with highest Pt
               // gavril: Should we use highest Pt or best Chi2 ?
               TrajectoryStateOnSurface one_tsos = vec_tracks_states[index_max_pt].state;
               
               LocalTrajectoryParameters ltp = one_tsos.localParameters();
               
               LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
               
               float locx = localDir.x();
               float locy = localDir.y();
               float locz = localDir.z();
               
               //alpha_ = atan2( locz, locx );
               //beta_  = atan2( locz, locy );
               
               cotalpha_ = locx/locz;
               cotbeta_  = locy/locz;
               
             } // if ( (int)vec_tracks_states.size() > 0 )
               
         } // if ( !useStraightTracks_ )
 
 
           // End: determine incident angles ============================================================
 
 
           // Calculate strip cluster charge and store amplitudes in vector for later use
           
           //cout << endl;
           //cout << "Calculate strip cluster charge and store amplitudes in vector for later use" << endl;
 
           float strip_cluster_charge = 0.0;
           std::vector<float> vec_cluster_charge; 
           vec_cluster_charge.clear();
           int cluster_size = (int)( (theStripCluster->amplitudes()).size() );
           
           int cluster_charge = 0;
           for (int i=0; i<cluster_size; ++i)
         {
           float current_strip_charge = (float)( (theStripCluster->amplitudes())[i] );
 
           strip_cluster_charge += current_strip_charge;
           vec_cluster_charge.push_back( current_strip_charge );  
         
           cluster_charge +=current_strip_charge; 
         }
          
        
           //cout << endl;
           //cout << "Calling strip qbin to see if the strip cluster has to be split..." << endl;          
           int strip_templQbin_ = strip_templ_.qbin( ID, cotalpha_, cotbeta_, strip_cluster_charge );
 
           if ( strip_templQbin_ < 0 || strip_templQbin_ > 5 )
         {
           // Do nothing...
           // cout << "Wrong strip strip_templQbin_ = " << strip_templQbin_ << endl;
         } // if ( strip_templQbin_ < 0 || strip_templQbin_ > 5 )
           else // if ( strip_templQbin_ < 0 || strip_templQbin_ > 5 ) {...} else
         {
           if ( strip_templQbin_ != 0 ) 
             {
               // Do not split this cluster. Do nothing.
               //cout << endl;
               //cout << "Do NOT split this cluster" << endl;
 
             } // if ( strip_templQbin_ != 0 ) 
           else // if ( templQbin_ != 0 ) {...} else
             {
               //cout << endl;
               //cout << "Split this cluster" << endl;
 
               // gavril : Is this OK ?
               uint16_t first_strip = theStripCluster->firstStrip();
  
               LocalVector lbfield = ( stripDet->surface() ).toLocal( magfield_->inTesla( stripDet->surface().position() ) ); 
               float locBy = lbfield.y();
               
               // Initialize values coming back from SiStripTemplateSplit::StripTempSplit
               float stripTemplXrec1_  = -99999.9;
               float stripTemplXrec2_  = -99999.9;
               float stripTemplSigmaX_ = -99999.9;
               float stripTemplProbX_  = -99999.9;
               int   stripTemplQbin_   = -99999;
               float stripTemplProbQ_  = -99999.9;
 
 
               /*
             cout << endl;
             cout << "About to call SiStripTemplateSplit::StripTempSplit(...)" << endl;
             
             cout << endl;
             cout << "ID        = " << ID        << endl;
             cout << "cotalpha_ = " << cotalpha_ << endl;
             cout << "cotbeta_  = " << cotbeta_  << endl;
             cout << "locBy     = " << locBy     << endl;
             cout << "Amplitudes: "; 
             for (int i=0; i<(int)vec_cluster_charge.size(); ++i)
             cout << vec_cluster_charge[i] << ", ";
             cout << endl;
               */
 
               int ierr =
             SiStripTemplateSplit::StripTempSplit(ID, 
                                  cotalpha_, cotbeta_, 
                                  locBy, 
                                  vec_cluster_charge, 
                                  strip_templ_, 
                                  stripTemplXrec1_, 
                                  stripTemplXrec2_, 
                                  stripTemplSigmaX_, 
                                  stripTemplProbX_,
                                  stripTemplQbin_,
                                  stripTemplProbQ_ );
               
              
 
               stripTemplXrec1_ += 2*strip_templ_.xsize();
               stripTemplXrec2_ += 2*strip_templ_.xsize();
               
             
 
               if ( ierr != 0 )
             {
               //cout << endl;
               //cout << "Strip cluster splitting failed: ierr = " << ierr << endl;
             }
               else // if ( ierr != 0 )
             {
               // Cluster was successfully split. 
               // Make the two split clusters and put them in the split cluster collection
               
               //cout << endl; 
               //cout << "Cluster was successfully split" << endl;
 
               cluster_was_successfully_split = true;
 
               std::vector<float> strip_cluster1;
               std::vector<float> strip_cluster2;
 
               strip_cluster1.clear();
               strip_cluster2.clear();
               
               // gavril : Is this OK ?!?!?!?! 
               for (int i=0; i<BXSIZE; ++i)
                 {
                   strip_cluster1.push_back(0.0);
                   strip_cluster2.push_back(0.0);
                 }
 
               //cout << endl;
               //cout << "About to call interpolate and sxtemp" << endl;
 
               strip_templ_.interpolate(ID, cotalpha_, cotbeta_, locBy);
               strip_templ_.sxtemp(stripTemplXrec1_, strip_cluster1);
               strip_templ_.sxtemp(stripTemplXrec2_, strip_cluster2);
              
              
               
               vector<SiStripDigi> vecSiStripDigi1;
               vecSiStripDigi1.clear();  
               int strip_cluster1_size = (int)strip_cluster1.size();
               for (int i=2; i<strip_cluster1_size; ++i)
                 {
                   if ( strip_cluster1[i] > 0.0 )
                 {
                   SiStripDigi current_digi1( first_strip + i-2, strip_cluster1[i] ); 
                   
                   vecSiStripDigi1.push_back( current_digi1 );
                 }
                 }
                
               
               
               vector<SiStripDigi> vecSiStripDigi2;
               vecSiStripDigi2.clear();  
               int strip_cluster2_size = (int)strip_cluster2.size();
               for (int i=2; i<strip_cluster2_size; ++i)
                 {
                   if ( strip_cluster2[i] > 0.0 )
                 {
                   SiStripDigi current_digi2( first_strip + i-2, strip_cluster2[i] ); 
                   
                   vecSiStripDigi2.push_back( current_digi2 );
                 }
                 }
               
               
             
               
               std::vector<SiStripDigi>::const_iterator SiStripDigiIterBegin1 = vecSiStripDigi1.begin();
               std::vector<SiStripDigi>::const_iterator SiStripDigiIterEnd1   = vecSiStripDigi1.end();
               std::pair<std::vector<SiStripDigi>::const_iterator, 
                 std::vector<SiStripDigi>::const_iterator> SiStripDigiRange1 
                 = make_pair(SiStripDigiIterBegin1, SiStripDigiIterEnd1);
 
               //if ( SiStripDigiIterBegin1 == SiStripDigiIterEnd1 )
               //{
               //  throw cms::Exception("TrackClusterSplitter : ") 
               //<< "\nERROR: SiStripDigiIterBegin1 = SiStripDigiIterEnd1 !!!" << "\n\n";
               //}
               
               std::vector<SiStripDigi>::const_iterator SiStripDigiIterBegin2 = vecSiStripDigi2.begin();
               std::vector<SiStripDigi>::const_iterator SiStripDigiIterEnd2   = vecSiStripDigi2.end();
               std::pair<std::vector<SiStripDigi>::const_iterator, 
                 std::vector<SiStripDigi>::const_iterator> SiStripDigiRange2 
                 = make_pair(SiStripDigiIterBegin2, SiStripDigiIterEnd2);
 
               // Sanity check...
               //if ( SiStripDigiIterBegin2 == SiStripDigiIterEnd2 )
               //{
               //  cout << endl;
               //  cout << "SiStripDigiIterBegin2 = SiStripDigiIterEnd2 !!!!!!!!!!!!!!!" << endl;
               //}
               
               
               // Save the split clusters
 
               if ( SiStripDigiIterBegin1 != SiStripDigiIterEnd1 )
                 {    
                   // gavril : Raw id ?
                   SiStripCluster cl1( detId.rawId(), SiStripDigiRange1 );
 
                   cl1.setSplitClusterError( stripTemplSigmaX_ );
 
                   output.push_back( cl1 );
                 
                   if ( (int)cl1.amplitudes().size() <= 0 )
                 {
                   throw cms::Exception("TrackClusterSplitter : ") 
                     << "\nERROR: (1) Wrong split cluster of size = " << (int)cl1.amplitudes().size() << "\n\n";
                 }
                   
                 } // if ( SiStripDigiIterBegin1 != SiStripDigiIterEnd1 )
 
               if ( SiStripDigiIterBegin2 != SiStripDigiIterEnd2 )
                 {
                   // gavril : Raw id ?
                   SiStripCluster cl2( detId.rawId(), SiStripDigiRange2 );
                   cl2.setSplitClusterError( stripTemplSigmaX_ );
                   output.push_back( cl2 ); 
         
                   if ( (int)cl2.amplitudes().size() <= 0 )
                 {
                   throw cms::Exception("TrackClusterSplitter : ") 
                     << "\nERROR: (2) Wrong split cluster of size = " << (int)cl2.amplitudes().size() << "\n\n";
                 }
                   
                 } // if ( SiStripDigiIterBegin2 != SiStripDigiIterEnd2 )
  
             
                
             } // else // if ( ierr != 0 )
 
             } // else // if ( strip_templQbin_ != 0 ) {...} else 
           
         } // else // if ( strip_templQbin_ < 0 || strip_templQbin_ > 5 ) {...} else
 
         } // else // if ( (int)theStripCluster->amplitudes().size() <= 1 )
 
       
       if ( !cluster_was_successfully_split )
         output.push_back( *c.cluster );
 
     } // if ( theStripCluster )
       else
     {
       throw cms::Exception("TrackClusterSplitter : ") 
         << "\nERROR: This is not a SiStripCluster  !!!" << "\n\n";
     }
 
     } // else if ( strip_tmpSplit_ )
   else
     {
       // gavril : Neither sim nor template splitter. Just add the cluster as it was. 
       output.push_back( *c.cluster );
     }
 }

template<>

void TrackClusterSplitter::splitCluster	(	const SiPixelClusterWithTracks &	c,
		const GlobalVector &	vtx,
		edmNew::DetSetVector< SiPixelCluster >::FastFiller &	output,
		DetId	detId
	)		const

private

Definition at line 1215 of file TrackClusterSplitter.cc.

References SiPixelCluster::add(), BXM2, BYM2, trackerHits::c, PixelDigi::channelToPixel(), SiPixelCluster::charge(), GetRecoTauVFromDQM_MC_cff::cl2, cluster_matrix_size_x, cluster_matrix_size_y, edm::DetSet< T >::data, edm::hlt::Exception, geometry_, i, PixelTopology::isItBigPixelInX(), PixelTopology::isItBigPixelInY(), j, TrajectoryStateOnSurface::localParameters(), Topology::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), SiPixelCluster::minPixelCol(), SiPixelCluster::minPixelRow(), LocalTrajectoryParameters::momentum(), convertSQLitetoXML_cfg::output, PixelSubdetector::PixelBarrel, pixeldigisimlink, PixelSubdetector::PixelEndcap, SiPixelCluster::pixels(), SiPixelTemplateSplit::PixelTempSplit(), reco::TrackBase::pt(), SiPixelTemplate::qbin(), SiPixelTemplate::s50(), SiPixelCluster::setSplitClusterErrorX(), SiPixelCluster::setSplitClusterErrorY(), SiPixelTemplate::simpletemplate2D(), simSplitPixel_, SiPixelCluster::size(), findQualityFiles::size, PixelGeomDetUnit::specificTopology(), mathSSE::sqrt(), GeomDet::surface(), templ2D_, templ_, tmpSplitPixel_, Surface::toGlobal(), TXSIZE, TYSIZE, useStraightTracks_, PV3DBase< T, PVType, FrameType >::x(), SiPixelCluster::x(), SiPixelTemplate::xsize(), SiPixelTemplate2D::xytemp(), PV3DBase< T, PVType, FrameType >::y(), SiPixelCluster::y(), SiPixelTemplate::ysize(), and PV3DBase< T, PVType, FrameType >::z().

 { 
   // The sim splitter:
   if ( simSplitPixel_ ) 
     {
       // cout << "Cluster splitting using simhits " << endl;
  
       int minPixelRow = (*c.cluster).minPixelRow();
       int maxPixelRow = (*c.cluster).maxPixelRow();
       int minPixelCol = (*c.cluster).minPixelCol();
       int maxPixelCol = (*c.cluster).maxPixelCol();    
       int dsl = 0; // number of digisimlinks
       
       edm::DetSetVector<PixelDigiSimLink>::const_iterator isearch = pixeldigisimlink->find(output.id());
       edm::DetSet<PixelDigiSimLink> digiLink = (*isearch);
       
       edm::DetSet<PixelDigiSimLink>::const_iterator linkiter = digiLink.data.begin();
       //create a vector for the track ids in the digisimlinks
       std::vector<int> simTrackIdV;  
       simTrackIdV.clear();
       //create a vector for the new splittedClusters 
       std::vector<SiPixelCluster> splittedCluster;
       splittedCluster.clear();
       
       for ( ; linkiter != digiLink.data.end(); linkiter++) 
     { // loop over all digisimlinks 
       dsl++;
       std::pair<int,int> pixel_coord = PixelDigi::channelToPixel(linkiter->channel());
       
       // is the digisimlink inside the cluster boundaries?
       if ( pixel_coord.first  <= maxPixelRow && 
            pixel_coord.first  >= minPixelRow &&
            pixel_coord.second <= maxPixelCol &&
            pixel_coord.second >= minPixelCol ) 
         {
           bool inStock(false); // did we see this simTrackId before?
          
           SiPixelCluster::PixelPos newPixelPos(pixel_coord.first, pixel_coord.second); // coordinates to the pixel
           
           //loop over the pixels from the cluster to get the charge in this pixel
           int newPixelCharge(0); //fraction times charge in the original cluster pixel
 
           const std::vector<SiPixelCluster::Pixel>& pixvector = (*c.cluster).pixels();
           
           for(std::vector<SiPixelCluster::Pixel>::const_iterator itPix = pixvector.begin(); itPix != pixvector.end(); itPix++)
         {
           if (((int) itPix->x) == ((int) pixel_coord.first)&&(((int) itPix->y) == ((int) pixel_coord.second)))
             {
               newPixelCharge = (int) (linkiter->fraction()*itPix->adc); 
             }
         }
           
           if ( newPixelCharge < 2500 ) 
         continue; 
           
           //add the pixel to an already existing cluster if the charge is above the threshold
           int clusVecPos = 0;
           std::vector<int>::const_iterator sTIter =  simTrackIdV.begin();
           
           for ( ; sTIter < simTrackIdV.end(); sTIter++) 
         {
           if (((*sTIter)== (int) linkiter->SimTrackId())) 
             {
               inStock=true; // now we saw this id before
               //      //          std::cout << " adding a pixel to the cluster " << (int) (clusVecPos) <<std::endl;
               //      //            std::cout << "newPixelCharge " << newPixelCharge << std::endl;
               splittedCluster.at(clusVecPos).add(newPixelPos,newPixelCharge); // add the pixel to the cluster
             }
           clusVecPos++;
         }
           
           //look if the splitted cluster was already made before, if not create one
           
           if ( !inStock ) 
         {
           //        std::cout << "creating a new cluster " << std::endl;
           simTrackIdV.push_back(linkiter->SimTrackId()); // add the track id to the vector
           splittedCluster.push_back(SiPixelCluster(newPixelPos,newPixelCharge)); // add the cluster to the vector
         }
         }
     }
       
       //    std::cout << "will add clusters : simTrackIdV.size() " << simTrackIdV.size() << std::endl;
       
       if ( ( ( (int)simTrackIdV.size() ) == 1 ) || ( *c.cluster).size()==1 ) 
     { 
       //        cout << "putting in this cluster" << endl;
       output.push_back(*c.cluster );
       //      std::cout << "cluster added " << output.size() << std::endl;
     }
       else 
     {     
       for (std::vector<SiPixelCluster>::const_iterator cIter = splittedCluster.begin(); cIter != splittedCluster.end(); cIter++ )
         {
           output.push_back( (*cIter) );   
         }
     }
   
       simTrackIdV.clear();  
       splittedCluster.clear();
     }
   else if ( tmpSplitPixel_ )
     { 
       bool cluster_was_successfully_split = false;
       
       const SiPixelCluster* thePixelCluster = static_cast<const SiPixelCluster*>(c.cluster);
 
       if ( thePixelCluster )
     { 
       // Do not attempt to split clusters of size one
       if ( (int)thePixelCluster->size() <= 1 )
         {    
           // Do nothing.
           //cout << "Will not attempt to split this clusters: " << endl;
           //cout << "(int)thePixelCluster->size() = " << (int)thePixelCluster->size() << endl;
         }
       else
         {
           // For barrel use template id 40 and for endcaps use template id 41
           int ID = -99999;
           if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel  )
         {
           //          cout << "We are in the barrel : " << (int)PixelSubdetector::PixelBarrel << endl;
           ID = 40;
         }
           else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
         {
           //          cout << "We are in the endcap : " << (int)PixelSubdetector::PixelEndcap << endl;
           ID = 41;
         }
           else 
         {
           // cout << "Not a pixel detector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
         }
           
 
           // Begin: determine incident angles ============================================================
 
           float cotalpha_ = -99999.9;
           float cotbeta_  = -99999.9;
 
           // First, determine track angles from straight track approximation
  
           // Find crude cluster center. 
           float xcenter = thePixelCluster->x();
           float ycenter = thePixelCluster->y();
           
           const GeomDetUnit* theDet = geometry_->idToDetUnit( detId );
           const PixelGeomDetUnit* pixDet = dynamic_cast<const PixelGeomDetUnit*>( theDet );
           
           const PixelTopology* theTopol = (&(pixDet->specificTopology()));
         
           // Transform from measurement to local coordinates (in cm)
           LocalPoint lp = theTopol->localPosition( MeasurementPoint(xcenter, ycenter) );
           
           // Transform from local to global coordinates
           GlobalPoint gp0 = theDet->surface().toGlobal( lp );
           
           // Make a vector pointing from the PV to the cluster center 
           GlobalPoint gp(gp0.x()-vtx.x(), gp0.y()-vtx.y(), gp0.z()-vtx.z() );
           
           // Make gp a unit vector, gv, pointing from the PV to the cluster center
           float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
           float gpx = gp.x()/gp_mod;
           float gpy = gp.y()/gp_mod;
           float gpz = gp.z()/gp_mod;
           GlobalVector gv(gpx, gpy, gpz);
 
           // Make unit vectors in local coordinates and then transform them in global coordinates
           const Local3DVector lvx(1.0, 0.0, 0.0);
           GlobalVector gvx = theDet->surface().toGlobal( lvx );
           const Local3DVector lvy(0.0, 1.0, 0.0);
           GlobalVector gvy = theDet->surface().toGlobal( lvy );
           const Local3DVector lvz(0.0, 0.0, 1.0);
           GlobalVector gvz = theDet->surface().toGlobal( lvz );
           
           // Calculate angles alpha and beta
           float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
           float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
           float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
           
           //float alpha_ = atan2( gv_dot_gvz, gv_dot_gvx );
           //float beta_  = atan2( gv_dot_gvz, gv_dot_gvy );
           
           cotalpha_ = gv_dot_gvx / gv_dot_gvz;
           cotbeta_  = gv_dot_gvy / gv_dot_gvz;
           
           
 
 
           // Attempt to get a better angle from tracks (either pixel tracks or full tracks)
           if ( !useStraightTracks_ )
         {
           // Use either pixel tracks (useTrajectories_ = False) or fully reconstructed tracks (useTrajectories_ = True)
           // When pixel/full tracks are not associated with the current cluster, will use angles from straight tracks
           
           // These are the tracks associated with this cluster
           std::vector<TrackAndState> vec_tracks_states = c.tracks;
           
           
 
           if ( (int)vec_tracks_states.size() > 0 )
             {
               //cout << "There is at least one track associated with this cluster. Pick the one with largest Pt." << endl;
               //cout << "(int)vec_tracks_states.size() = " << (int)vec_tracks_states.size() << endl;
               
               int index_max_pt = -99999;  // index of the track with the highest Pt
               float max_pt = -99999.9;
               
               for (int i=0; i<(int)vec_tracks_states.size(); ++i )
             {
               const reco::Track* one_track = vec_tracks_states[i].track;
               
               if ( one_track->pt() > max_pt )
                 {
                   index_max_pt = i;
                   max_pt = one_track->pt();
                 }  
             }
               
               // Pick the tsos from the track with highest Pt
               // gavril: Should we use highest Pt or best Chi2 ?
               TrajectoryStateOnSurface one_tsos = vec_tracks_states[index_max_pt].state;
               
               LocalTrajectoryParameters ltp = one_tsos.localParameters();
               
               LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
               
               float locx = localDir.x();
               float locy = localDir.y();
               float locz = localDir.z();
               
               //alpha_ = atan2( locz, locx );
               //beta_  = atan2( locz, locy );
               
               cotalpha_ = locx/locz;
               cotbeta_  = locy/locz;
               
              
               
             } // if ( (int)vec_tracks_states.size() > 0 )
               
         } // if ( !useStraightTracks_ )
 
           // End: determine incident angles ============================================================
 
 
           
           //cout << "Calling qbin to see if the cluster has to be split..." << endl;          
           int templQbin_ = templ_.qbin( ID, cotalpha_, cotbeta_, thePixelCluster->charge() );
 
           if ( templQbin_ < 0 || templQbin_ > 5  )
         {
           // gavril : check this....
           // cout << "Template call failed. Cannot decide if cluster should be split !!!!!!! " << endl;
           // cout << "Do nothing." << endl;
         }
           else // if ( templQbin_ < 0 || templQbin_ > 5  ) {...} else
                 {
           //cout << " Returned OK from PixelTempReco2D..." << endl;
           
           // Check for split clusters: we split the clusters with larger than expected charge: templQbin_ == 0 
           if ( templQbin_ != 0 ) 
             {
               // gavril: do not split this cluster
               //cout << "TEMPLATE SPLITTER SAYS : NO SPLIT " << endl;
               //cout << "This cluster will note be split !!!!!!!!!! " << endl; 
             }
           else // if ( templQbin_ != 0 ) {...} else
             {
               //cout << "TEMPLATE SPLITTER SAYS : SPLIT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
               //cout << "Found a cluster that has to be split. templQbin_ = " << templQbin_ << endl;
               
               // gavril: Call the template splitter
               //cout << "Calling the splitter..." << endl;
              
               // Put the pixels of this clusters in a 2x2 array to be used by the template splitter
               
               const std::vector<SiPixelCluster::Pixel> & pixVec = thePixelCluster->pixels();
               std::vector<SiPixelCluster::Pixel>::const_iterator pixIter = pixVec.begin(), pixEnd = pixVec.end();
               
               const int cluster_matrix_size_x = 13;
               const int cluster_matrix_size_y = 21;
               
               boost::multi_array<float, 2> clust_array_2d(boost::extents[cluster_matrix_size_x][cluster_matrix_size_y]);
               
               int row_offset = thePixelCluster->minPixelRow();
               int col_offset = thePixelCluster->minPixelCol();
               
               // Copy clust's pixels (calibrated in electrons) into clust_array_2d;
           
               for ( ; pixIter != pixEnd; ++pixIter ) 
             {
               int irow = int(pixIter->x) - row_offset;   // do we need +0.5 ???
               int icol = int(pixIter->y) - col_offset;   // do we need +0.5 ???
               
               if ( irow<cluster_matrix_size_x && icol<cluster_matrix_size_y )
                 {
                   clust_array_2d[irow][icol] = (float)pixIter->adc;
                 }
             }
               
               // Make and fill the bool arrays flagging double pixels
               std::vector<bool> ydouble(cluster_matrix_size_y), xdouble(cluster_matrix_size_x);
               
               // x directions (shorter), rows
               for (int irow = 0; irow < cluster_matrix_size_x; ++irow)
             {
               xdouble[irow] = theTopol->isItBigPixelInX( irow+row_offset );
             }
               
               // y directions (longer), columns
               for (int icol = 0; icol < cluster_matrix_size_y; ++icol) 
             {
               ydouble[icol] = theTopol->isItBigPixelInY( icol+col_offset );
             }
 
               // gavril: Initialize values coming back from SiPixelTemplateSplit::PixelTempSplit
               float templYrec1_  = -99999.9;
               float templYrec2_  = -99999.9;
               float templSigmaY_ = -99999.9;
               float templProbY_  = -99999.9;
               float templXrec1_  = -99999.9;
               float templXrec2_  = -99999.9;
               float templSigmaX_ = -99999.9;
               float templProbX_  = -99999.9;
               float dchisq       = -99999.9;
               float templProbQ_  = -99999.9;
 
               int ierr =
             SiPixelTemplateSplit::PixelTempSplit( ID, 
                                   cotalpha_, cotbeta_,
                                   clust_array_2d, 
                                   ydouble, xdouble,
                                   templ_,
                                   templYrec1_, templYrec2_, templSigmaY_, templProbY_,
                                   templXrec1_, templXrec2_, templSigmaX_, templProbX_,
                                   templQbin_, 
                                   templProbQ_, 
                                   true, 
                                   dchisq, 
                                   templ2D_ );
                 
               if ( ierr != 0 )
             {
               // cout << "Cluster splitting failed: ierr = " << ierr << endl;
             }
               else
             {
               // gavril: Cluster was successfully split. 
               // gavril: Make the two split clusters and put them in the split cluster collection
               //cout << "Cluster splitting OK: ierr = " << ierr << endl;
               
               // 2D templates have origin at the lower left corner of template2d[1][1] which is 
               // also 2 pixels larger than the cluster container
               
               float xsize = templ_.xsize();  // this is the pixel x-size in microns
               float ysize = templ_.ysize();  // this is the pixel y-size in microns
               
               // Shift the coordinates to the 2-D template system             
               float yrecp1 = -99999.9;
               float yrecp2 = -99999.9;
               float xrecp1 = -99999.9;
               float xrecp2 = -99999.9;
               
               if ( ydouble[0] ) 
                 {               
                   yrecp1 = templYrec1_ + ysize;
                   yrecp2 = templYrec2_ + ysize;
                 } 
               else 
                 {
                   yrecp1 = templYrec1_ + ysize/2.0;
                   yrecp2 = templYrec2_ + ysize/2.0;
                 }
               
               if ( xdouble[0] ) 
                 {
                   xrecp1 = templXrec1_ + xsize;
                   xrecp2 = templXrec2_ + xsize;
                 } 
               else 
                 {
                   xrecp1 = templXrec1_ + xsize/2.0;
                   xrecp2 = templXrec2_ + xsize/2.0;
                 }
               
               //  The xytemp method adds charge to a zeroed buffer
               
               float template2d1[BXM2][BYM2];
               float template2d2[BXM2][BYM2];
               
               for ( int j=0; j<BXM2; ++j ) 
                 {
                   for ( int i=0; i<BYM2; ++i ) 
                 {
                   template2d1[j][i] = 0.0;
                   template2d2[j][i] = 0.0;
                 }
                 }
                
    
               bool template_OK 
                 = templ2D_.xytemp(ID, cotalpha_, cotbeta_, 
                           xrecp1, yrecp1, 
                           ydouble, xdouble, 
                           template2d1);
               
               template_OK 
                 = template_OK && 
                 templ2D_.xytemp(ID, cotalpha_, cotbeta_, 
                         xrecp2, yrecp2, 
                         ydouble, xdouble, 
                         template2d2);
               
               if ( !template_OK ) 
                 {
                   // gavril: check this
                   //cout << "Template is not OK. Fill out with zeros !!!!!!!!!!!!!!! " << endl; 
                   
                   for ( int j=0; j<BXM2; ++j ) 
                 {
                   for ( int i=0; i<BYM2; ++i ) 
                     {
                       template2d1[j][i] = 0.0;
                       template2d2[j][i] = 0.0;
                     }
                 }
                   
                   if ( !templ_.simpletemplate2D(xrecp1, yrecp1, 
                                 xdouble, ydouble, 
                                 template2d1) ) 
                 {
                   cluster_was_successfully_split = false;
                 }
                   
                   if ( !templ_.simpletemplate2D(xrecp2, yrecp2, 
                                 xdouble, ydouble, 
                                 template2d2) ) 
                 {
                   cluster_was_successfully_split = false;
                 }
                   
                 } // if ( !template_OK ) 
               else
                 {
                   cluster_was_successfully_split = true;
                   
                   // Next, copy the 2-d templates into cluster containers, replace small signals with zero.
                   // Cluster 1 and cluster 2 should line up with clust_array_2d (same origin and pixel indexing)
                   
                   float q50 = templ_.s50();
                   
                   
                   float cluster1[TXSIZE][TYSIZE];
                   float cluster2[TXSIZE][TYSIZE];  //Note that TXSIZE = BXM2 - 2, TYSIZE = BYM2 - 2
                   
                   for ( int j=0; j<TXSIZE; ++j ) 
                 {
                   for ( int i=0; i<TYSIZE; ++i ) 
                     {
                       cluster1[j][i] = template2d1[j+1][i+1];
                       
                       if ( cluster1[j][i] < q50 ) 
                     cluster1[j][i] = 0.0;
                       
                       cluster2[j][i] = template2d2[j+1][i+1];
                       
                       if ( cluster2[j][i] < q50 ) 
                     cluster2[j][i] = 0.0;
 
                       //cout << "cluster1[j][i] = " << cluster1[j][i] << endl;
                       //cout << "cluster2[j][i] = " << cluster2[j][i] << endl;
                     }
                 }
                 
               
                   // Find the coordinates of first pixel in each of the two split clusters 
                   int i1 = -999;
                   int j1 = -999;
                   int i2 = -999;
                   int j2 = -999;
                                  
                   bool done_searching = false; 
                   for ( int i=0; i<13 && !done_searching; ++i)
                 {
                   for (int j=0; j<21 && !done_searching; ++j)
                     {
                       if ( cluster1[i][j] > 0 )
                     {
                       i1 = i;
                       j1 = j;
                       done_searching = true;
                     }    
                     }
                 }
 
                   done_searching = false; 
                   for ( int i=0; i<13 && !done_searching; ++i)
                 {
                   for (int j=0; j<21 && !done_searching; ++j)
                     {
                       if ( cluster2[i][j] > 0 )
                     {
                       i2 = i;
                       j2 = j;
                       done_searching = true;
                     }
                     }
                 }
                   
                   
                   // Make clusters out of the first pixels in each of the two split clsuters
 
                   SiPixelCluster cl1( SiPixelCluster::PixelPos( i1 + row_offset, j1 + col_offset), 
                           cluster1[i1][j1] );
                   
                   SiPixelCluster cl2( SiPixelCluster::PixelPos( i2 + row_offset, j2 + col_offset), 
                           cluster2[i2][j2] );
                   
 
                   // Now add the rest of the pixels to the split clusters
 
                   for ( int i=0; i<13; ++i)
                 {
                   for (int j=0; j<21; ++j)
                     {
                       
                       if ( cluster1[i][j] > 0.0 && (i!=i1 || j!=j1 ) )
                     {
                       cl1.add( SiPixelCluster::PixelPos( i + row_offset, j + col_offset), 
                            cluster1[i][j] ); 
 
                       //cout << "cluster1[i][j] = " << cluster1[i][j] << endl;
                     }
                       
                       
                       if ( cluster2[i][j] > 0.0 && (i!=i2 || j!=j2 ) )
                     {                             
                       cl2.add( SiPixelCluster::PixelPos( i + row_offset, j + col_offset), 
                            cluster2[i][j] );
                       
                       //cout << "cluster2[i][j] = " << cluster2[i][j] << endl;
                     }
                     }
                 }
                   
                   // Attach errors and probabilities to the split Clusters
                   // The errors will be later attahed to the SiPixelRecHit
 
 
                   cl1.setSplitClusterErrorX( templSigmaX_ );
                   cl1.setSplitClusterErrorY( templSigmaY_ );
                   //cl1.prob_q = templProbQ_;
                   
                   cl2.setSplitClusterErrorX( templSigmaX_ );
                   cl2.setSplitClusterErrorY( templSigmaY_ );
                   //cl2.prob_q = templProbQ_;
 
 
                   // Save the split clusters
                   output.push_back( cl1 );
                   output.push_back( cl2 );     
                  
                   // Some sanity checks
 
                   if ( (int)cl1.size() <= 0 )
                 {
                   edm::LogError("TrackClusterSplitter : ") 
                     << "1) Cluster of size = " << (int)cl1.size() << " !!! " << endl;
                 }
 
                   if ( (int)cl2.size() <= 0 )
                 {
                   edm::LogError("TrackClusterSplitter : ") 
                     << "2) Cluster of size = " << (int)cl2.size() << " !!! " << endl;
                 }
 
                   if ( cl1.charge() <= 0 )
                 {
                   edm::LogError("TrackClusterSplitter : ") 
                     << "1) Cluster of charge = " << (int)cl1.charge() << " !!! " << endl;
                   
                 }
 
                   if ( cl2.charge() <= 0 )
                 {
                   edm::LogError("TrackClusterSplitter : ") 
                     << "2) Cluster of charge = " << (int)cl2.charge() << " !!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
                 }
                   
                  
                 } // if ( !template_OK ) ... else 
               
             } // if ( ierr2 != 0 ) ... else
               
             } // if ( templQbin_ != 0 ) ... else
           
         } // else // if ( templQbin_ < 0 || templQbin_ > 5  ) {...} else
           
         } // if ( (int)thePixelCluster->size() <= 1 ) {... } else
       
       if ( !cluster_was_successfully_split )
         output.push_back(*c.cluster);
       
     } // if ( theSiPixelCluster )
       else 
     {
       throw cms::Exception("TrackClusterSplitter :") 
         << "This is not a SiPixelCluster !!! " << "\n"; 
     }
     }
   else 
     {
       // gavril : Neither sim nor template splitter. Just add the cluster as it was.
       // original from G.P.
       output.push_back( *c.cluster );
     }
 }

template<typename Cluster >

std::auto_ptr< edmNew::DetSetVector< Cluster > > TrackClusterSplitter::splitClusters	(	const std::map< uint32_t, boost::sub_range< std::vector< ClusterWithTracks< Cluster > > > > &	input,
		const reco::Vertex &	vtx
	)		const

private

Definition at line 541 of file TrackClusterSplitter.cc.

References trackerHits::c, geometry_, LaserDQM_cfg::input, convertSQLitetoXML_cfg::output, AlCaHLTBitMon_ParallelJobs::p, splitCluster(), reco::Vertex::x(), reco::Vertex::y(), and reco::Vertex::z().

Referenced by produce().

 {
   std::auto_ptr<edmNew::DetSetVector<Cluster> > output(new edmNew::DetSetVector<Cluster>());
   typedef std::pair<uint32_t, boost::sub_range<std::vector<ClusterWithTracks<Cluster> > > > pair;
   
   foreach(const pair &p, input) 
     {
       const GeomDet* det = geometry_->idToDet( DetId(p.first) );
       
       if ( det == 0 ) 
         { 
           edm::LogError("MissingDetId") << "DetIDs " << p.first << " is not in geometry.\n";
       continue;
     }
       
       // gavril: Pass the PV instead of direction
       // GlobalVector dir(det->position().x() - vtx.x(), det->position().y() - vtx.y(), det->position().z() - vtx.z());
       GlobalVector primary_vtx( vtx.x(), vtx.y(), vtx.z() );
 
       // Create output collection
       typename edmNew::DetSetVector<Cluster>::FastFiller detset(*output, p.first);
       
       // fill it
       foreach(const ClusterWithTracks<Cluster> &c, p.second) 
     {
       splitCluster(c, primary_vtx, detset, DetId(p.first) );
         }
     }
 
   return output;
 }