|
|
|
Definition at line 64 of file TrackClusterSplitter.cc.
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.
| TrackClusterSplitter::TrackClusterSplitter | ( | const edm::ParameterSet & | iConfig | ) |
Definition at line 231 of file TrackClusterSplitter.cc.
References edm::ParameterSet::getParameter(), propagatorName_, SiPixelTemplate::pushfile(), SiPixelTemplate2D::pushfile(), SiStripTemplate::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.
{
}
| static const C* TrackClusterSplitter::equalClusters | ( | const C & | c1, |
| const C & | c2 | ||
| ) | [inline, static, private] |
Definition at line 156 of file TrackClusterSplitter.cc.
{
return false;
}
| const SiPixelCluster * TrackClusterSplitter::getCluster< SiPixelCluster > | ( | const TrackingRecHit * | hit | ) | [static, private] |
Definition at line 288 of file TrackClusterSplitter.cc.
References 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";
}
| 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(), 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(), edmNew::DetSet< T >::end(), end, 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(), 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();
}
| void TrackClusterSplitter::splitCluster | ( | const SiPixelClusterWithTracks & | cluster, |
| const GlobalVector & | dir, | ||
| edmNew::DetSetVector< SiPixelCluster >::FastFiller & | output, | ||
| DetId | detId | ||
| ) | const [private] |
| void TrackClusterSplitter::splitCluster | ( | const SiStripClusterWithTracks & | cluster, |
| const GlobalVector & | dir, | ||
| edmNew::DetSetVector< SiStripCluster >::FastFiller & | output, | ||
| DetId | detId | ||
| ) | const [private] |
| 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, PixelDigi::channelToPixel(), SiPixelCluster::charge(), cluster_matrix_size_x, cluster_matrix_size_y, edm::DetSet< T >::data, Exception, geometry_, i, PixelTopology::isItBigPixelInX(), PixelTopology::isItBigPixelInY(), j, TrajectoryStateOnSurface::localParameters(), Topology::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), SiPixelCluster::minPixelCol(), SiPixelCluster::minPixelRow(), LocalTrajectoryParameters::momentum(), PixelSubdetector::PixelBarrel, pixeldigisimlink, PixelSubdetector::PixelEndcap, SiPixelCluster::pixels(), SiPixelTemplateSplit::PixelTempSplit(), reco::TrackBase::pt(), SiPixelTemplate::qbin(), SiPixelTemplate::s50(), SiPixelCluster::setSplitClusterErrorX(), SiPixelCluster::setSplitClusterErrorY(), SiPixelTemplate::simpletemplate2D(), simSplitPixel_, findQualityFiles::size, SiPixelCluster::size(), PixelGeomDetUnit::specificTopology(), mathSSE::sqrt(), GeomDet::surface(), templ2D_, templ_, tmpSplitPixel_, Surface::toGlobal(), TXSIZE, TYSIZE, useStraightTracks_, SiPixelCluster::x(), PV3DBase< T, PVType, FrameType >::x(), SiPixelTemplate::xsize(), SiPixelTemplate2D::xytemp(), SiPixelCluster::y(), PV3DBase< T, PVType, FrameType >::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 );
}
}
| 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, edm::DetSet< T >::data, end, 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(), 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 );
}
}
| 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 Exception.
Referenced by splitClusters().
{
//cout << "Should never be here: TrackClusterSplitter, TrackClusterSplitter::splitCluster(...) !!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
throw cms::Exception("LogicError", "This should not be called");
}
| 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_, collect_tpl::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;
}
std::vector<SiPixelClusterWithTracks> TrackClusterSplitter::allSiPixelClusters [private] |
std::vector<SiStripClusterWithTracks> TrackClusterSplitter::allSiStripClusters [private] |
Definition at line 204 of file TrackClusterSplitter.cc.
Referenced by produce().
Definition at line 105 of file TrackClusterSplitter.cc.
Referenced by produce(), splitCluster(), and splitClusters().
Definition at line 150 of file TrackClusterSplitter.cc.
Referenced by produce(), and splitCluster().
Definition at line 103 of file TrackClusterSplitter.cc.
Referenced by produce(), and splitCluster().
Definition at line 74 of file TrackClusterSplitter.cc.
Referenced by produce().
Definition at line 108 of file TrackClusterSplitter.cc.
Referenced by produce(), and splitCluster().
Definition at line 104 of file TrackClusterSplitter.cc.
Referenced by produce().
std::string TrackClusterSplitter::propagatorName_ [private] |
Definition at line 102 of file TrackClusterSplitter.cc.
Referenced by TrackClusterSplitter().
bool TrackClusterSplitter::simSplitPixel_ [private] |
Definition at line 76 of file TrackClusterSplitter.cc.
Referenced by produce(), splitCluster(), and TrackClusterSplitter().
bool TrackClusterSplitter::simSplitStrip_ [private] |
Definition at line 77 of file TrackClusterSplitter.cc.
Referenced by produce(), splitCluster(), and TrackClusterSplitter().
std::map<uint32_t, SiPixelClustersWithTracks> TrackClusterSplitter::siPixelDetsWithClusters [private] |
Definition at line 202 of file TrackClusterSplitter.cc.
Referenced by produce().
std::map<uint32_t, SiStripClustersWithTracks> TrackClusterSplitter::siStripDetsWithClusters [private] |
Definition at line 205 of file TrackClusterSplitter.cc.
Referenced by produce().
SiStripTemplate TrackClusterSplitter::strip_templ_ [mutable, private] |
Definition at line 119 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
Definition at line 73 of file TrackClusterSplitter.cc.
Referenced by produce().
Definition at line 111 of file TrackClusterSplitter.cc.
Referenced by produce(), and splitCluster().
SiPixelTemplate2D TrackClusterSplitter::templ2D_ [mutable, private] |
Definition at line 117 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
SiPixelTemplate TrackClusterSplitter::templ_ [mutable, private] |
Definition at line 116 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
bool TrackClusterSplitter::tmpSplitPixel_ [private] |
Definition at line 78 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
bool TrackClusterSplitter::tmpSplitStrip_ [private] |
Definition at line 79 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
Definition at line 96 of file TrackClusterSplitter.cc.
Referenced by produce().
bool TrackClusterSplitter::useStraightTracks_ [private] |
Definition at line 87 of file TrackClusterSplitter.cc.
Referenced by splitCluster(), and TrackClusterSplitter().
bool TrackClusterSplitter::useTrajectories_ [private] |
Definition at line 93 of file TrackClusterSplitter.cc.
Referenced by produce(), and TrackClusterSplitter().
edm::InputTag TrackClusterSplitter::vertices_ [private] |
Definition at line 99 of file TrackClusterSplitter.cc.
Referenced by produce().
1.7.3