CMS 3D CMS Logo

SiPixelPhase1TrackClusters.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: SiPixelPhase1TrackClusters
4 // Class: SiPixelPhase1TrackClusters
5 //
6 
7 // Original Author: Marcel Schneider
8 
11 
24 
25 
27  SiPixelPhase1Base(iConfig)
28 {
29  clustersToken_ = consumes<edmNew::DetSetVector<SiPixelCluster>>(iConfig.getParameter<edm::InputTag>("clusters"));
30  tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
31 }
32 
34 
35  // get geometry
37  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
38  assert(tracker.isValid());
39 
40  //get the map
42  iEvent.getByToken( tracksToken_, tracks);
43  if ( !tracks.isValid() ) {
44  edm::LogWarning("SiPixelPhase1TrackClusters") << "track collection is not valid";
45  return;
46  }
47 
48  // get clusters
50  iEvent.getByToken( clustersToken_, clusterColl );
51  if ( !clusterColl.isValid() ) {
52  edm::LogWarning("SiPixelPhase1TrackClusters") << "pixel cluster collection is not valid";
53  return;
54  }
55 
56  // we need to store some per-cluster data. Instead of a map, we use a vector,
57  // exploiting the fact that all custers live in the DetSetVector and we can
58  // use the same indices to refer to them.
59  // corr_charge is not strictly needed but cleaner to have it.
60  std::vector<bool> ontrack (clusterColl->data().size(), false);
61  std::vector<float> corr_charge(clusterColl->data().size(), -1.0f);
62 
63  for (auto const & track : *tracks) {
64 
65  bool isBpixtrack = false, isFpixtrack = false, crossesPixVol=false;
66 
67  // find out whether track crosses pixel fiducial volume (for cosmic tracks)
68  double d0 = track.d0(), dz = track.dz();
69  if(std::abs(d0)<15 && std::abs(dz)<50) crossesPixVol = true;
70 
71  auto const & trajParams = track.extra()->trajParams();
72  assert(trajParams.size()==track.recHitsSize());
73  auto hb = track.recHitsBegin();
74  for(unsigned int h=0;h<track.recHitsSize();h++){
75  auto hit = *(hb+h);
76  if (!hit->isValid()) continue;
77  DetId id = hit->geographicalId();
78 
79  // check that we are in the pixel
80  uint32_t subdetid = (id.subdetId());
81  if (subdetid == PixelSubdetector::PixelBarrel) isBpixtrack = true;
82  if (subdetid == PixelSubdetector::PixelEndcap) isFpixtrack = true;
83  if (subdetid != PixelSubdetector::PixelBarrel && subdetid != PixelSubdetector::PixelEndcap) continue;
84  auto pixhit = dynamic_cast<const SiPixelRecHit*>(hit->hit());
85  if (!pixhit) continue;
86 
87  // get the cluster
88  auto clust = pixhit->cluster();
89  if (clust.isNull()) continue;
90  ontrack[clust.key()] = true; // mark cluster as ontrack
91 
92 
93  // correct charge for track impact angle
94  auto const & ltp = trajParams[h];
95  LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
96 
97  float clust_alpha = atan2(localDir.z(), localDir.x());
98  float clust_beta = atan2(localDir.z(), localDir.y());
99  double corrCharge = clust->charge() * sqrt( 1.0 / ( 1.0/pow( tan(clust_alpha), 2 ) +
100  1.0/pow( tan(clust_beta ), 2 ) +
101  1.0 ));
102  corr_charge[clust.key()] = (float) corrCharge;
103  }
104 
105  // statistics on tracks
106  histo[NTRACKS].fill(1, DetId(0), &iEvent);
107  if (isBpixtrack || isFpixtrack)
108  histo[NTRACKS].fill(2, DetId(0), &iEvent);
109  if (isBpixtrack)
110  histo[NTRACKS].fill(3, DetId(0), &iEvent);
111  if (isFpixtrack)
112  histo[NTRACKS].fill(4, DetId(0), &iEvent);
113 
114  if (crossesPixVol) {
115  if (isBpixtrack || isFpixtrack)
116  histo[NTRACKS_VOLUME].fill(1, DetId(0), &iEvent);
117  else
118  histo[NTRACKS_VOLUME].fill(0, DetId(0), &iEvent);
119  }
120  }
121 
123  for (it = clusterColl->begin(); it != clusterColl->end(); ++it) {
124  auto id = DetId(it->detId());
125 
126  const PixelGeomDetUnit* geomdetunit = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(id) );
127  const PixelTopology& topol = geomdetunit->specificTopology();
128 
129  for(auto subit = it->begin(); subit != it->end(); ++subit) {
130  // we could do subit-...->data().front() as well, but this seems cleaner.
131  auto key = edmNew::makeRefTo(clusterColl, subit).key();
132  bool is_ontrack = ontrack[key];
133  float corrected_charge = corr_charge[key];
134  SiPixelCluster const& cluster = *subit;
135 
136  LocalPoint clustlp = topol.localPosition(MeasurementPoint(cluster.x(), cluster.y()));
137  GlobalPoint clustgp = geomdetunit->surface().toGlobal(clustlp);
138 
139  if (is_ontrack) {
140  histo[ONTRACK_NCLUSTERS ].fill(id, &iEvent);
141  histo[ONTRACK_CHARGE ].fill(double(corrected_charge), id, &iEvent);
142  histo[ONTRACK_SIZE ].fill(double(cluster.size() ), id, &iEvent);
143  histo[ONTRACK_POSITION_B].fill(clustgp.z(), clustgp.phi(), id, &iEvent);
144  histo[ONTRACK_POSITION_F].fill(clustgp.x(), clustgp.y(), id, &iEvent);
145  } else {
146  histo[OFFTRACK_NCLUSTERS ].fill(id, &iEvent);
147  histo[OFFTRACK_CHARGE ].fill(double(cluster.charge()), id, &iEvent);
148  histo[OFFTRACK_SIZE ].fill(double(cluster.size() ), id, &iEvent);
149  histo[OFFTRACK_POSITION_B].fill(clustgp.z(), clustgp.phi(), id, &iEvent);
150  histo[OFFTRACK_POSITION_F].fill(clustgp.x(), clustgp.y(), id, &iEvent);
151  }
152  }
153  }
154 
155  histo[ONTRACK_NCLUSTERS].executePerEventHarvesting(&iEvent);
156  histo[OFFTRACK_NCLUSTERS].executePerEventHarvesting(&iEvent);
157 }
158 
160 
T getParameter(std::string const &) const
edm::Ref< typename HandleT::element_type, typename HandleT::element_type::value_type::value_type > makeRefTo(const HandleT &iHandle, typename HandleT::element_type::value_type::const_iterator itIter)
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
edm::EDGetTokenT< reco::TrackCollection > tracksToken_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
int charge() const
SiPixelPhase1TrackClusters(const edm::ParameterSet &conf)
Measurement2DPoint MeasurementPoint
Measurement points are two-dimensional by default.
int iEvent
Definition: GenABIO.cc:230
T mag() const
Definition: PV3DBase.h:67
T sqrt(T t)
Definition: SSEVec.h:18
void analyze(const edm::Event &, const edm::EventSetup &)
T z() const
Definition: PV3DBase.h:64
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool isValid() const
Definition: HandleBase.h:74
Definition: DetId.h:18
ClusterRef cluster() const
Definition: SiPixelRecHit.h:49
const T & get() const
Definition: EventSetup.h:56
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
std::vector< HistogramManager > histo
Pixel cluster – collection of neighboring pixels above threshold.
float y() const
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > clustersToken_
bool isValid() const
Definition: ESHandle.h:47
T x() const
Definition: PV3DBase.h:62
float x() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
const_iterator begin(bool update=false) const
const TrackerGeomDet * idToDet(DetId) const
Our base class.
Definition: SiPixelRecHit.h:23
int size() const