CMS 3D CMS Logo

MultiTrackValidatorGenPs.cc
Go to the documentation of this file.
2 
5 
20 
24 
25 #include "TMath.h"
26 #include <TF1.h>
27 
28 //#include <iostream>
29 
30 using namespace std;
31 using namespace edm;
32 
33 static const std::string kTrackAssociatorByChi2("trackAssociatorByChi2");
34 
36  gpSelector = GenParticleCustomSelector(pset.getParameter<double>("ptMinGP"),
37  pset.getParameter<double>("minRapidityGP"),
38  pset.getParameter<double>("maxRapidityGP"),
39  pset.getParameter<double>("tipGP"),
40  pset.getParameter<double>("lipGP"),
41  pset.getParameter<bool>("chargedOnlyGP"),
42  pset.getParameter<int>("statusGP"),
43  pset.getParameter<std::vector<int> >("pdgIdGP"));
44 
45  if (useAssociators_) {
46  for (auto const& src : associators) {
47  if (src.label() == kTrackAssociatorByChi2) {
48  label_gen_associator = consumes<reco::TrackToGenParticleAssociator>(src);
49  break;
50  }
51  }
52  } else {
53  for (auto const& src : associators) {
54  associatormapGtR = consumes<reco::GenToRecoCollection>(src);
55  associatormapRtG = consumes<reco::RecoToGenCollection>(src);
56  break;
57  }
58  }
59 }
60 
62 
64  const edm::EventSetup& setup,
65  const Histograms& histograms) const {
66  using namespace reco;
67 
68  edm::LogInfo("TrackValidator") << "\n===================================================="
69  << "\n"
70  << "Analyzing new event"
71  << "\n"
72  << "====================================================\n"
73  << "\n";
74 
75  const TrackerTopology& ttopo = setup.getData(tTopoEsToken);
76 
77  edm::Handle<GenParticleCollection> TPCollectionHeff;
78  event.getByToken(label_tp_effic, TPCollectionHeff);
79  const GenParticleCollection tPCeff = *(TPCollectionHeff.product());
80 
81  edm::Handle<GenParticleCollection> TPCollectionHfake;
82  event.getByToken(label_tp_fake, TPCollectionHfake);
83  const GenParticleCollection tPCfake = *(TPCollectionHfake.product());
84 
85  //if (tPCeff.size()==0) {edm::LogInfo("TrackValidator")
86  //<< "TP Collection for efficiency studies has size = 0! Skipping Event." ; return;}
87  //if (tPCfake.size()==0) {edm::LogInfo("TrackValidator")
88  //<< "TP Collection for fake rate studies has size = 0! Skipping Event." ; return;}
89 
90  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
91  event.getByToken(bsSrc, recoBeamSpotHandle);
92  reco::BeamSpot bs = *recoBeamSpotHandle;
93 
95  event.getByToken(label_pileupinfo, puinfoH);
96  PileupSummaryInfo puinfo;
97 
98  for (unsigned int puinfo_ite = 0; puinfo_ite < (*puinfoH).size(); ++puinfo_ite) {
99  if ((*puinfoH)[puinfo_ite].getBunchCrossing() == 0) {
100  puinfo = (*puinfoH)[puinfo_ite];
101  break;
102  }
103  }
104 
105  const reco::TrackToGenParticleAssociator* trackGenAssociator = nullptr;
106  if (useAssociators_) {
108  return;
109  } else {
111  event.getByToken(label_gen_associator, trackGenAssociatorH);
112  trackGenAssociator = trackGenAssociatorH.product();
113  }
114  } else if (associatormapGtR.isUninitialized()) {
115  return;
116  }
117 
118  // dE/dx
119  // at some point this could be generalized, with a vector of tags and a corresponding vector of Handles
120  // I'm writing the interface such to take vectors of ValueMaps
121  std::vector<const edm::ValueMap<reco::DeDxData>*> v_dEdx;
122  if (dodEdxPlots_) {
125  event.getByToken(m_dEdx1Tag, dEdx1Handle);
126  event.getByToken(m_dEdx2Tag, dEdx2Handle);
127  v_dEdx.push_back(dEdx1Handle.product());
128  v_dEdx.push_back(dEdx2Handle.product());
129  }
130 
131  std::vector<float> mvaDummy;
132 
133  int w = 0; //counter counting the number of sets of histograms
134  for (unsigned int www = 0; www < label.size(); www++) {
135  //
136  //get collections from the event
137  //
140  continue;
141  //if (trackCollection->size()==0)
142  //edm::LogInfo("TrackValidator") << "TrackCollection size = 0!" ;
143  //continue;
144  //}
145  reco::RecoToGenCollection recGenColl;
146  reco::GenToRecoCollection genRecColl;
147  //associate tracks
148  if (useAssociators_) {
149  edm::LogVerbatim("TrackValidator") << "Analyzing " << label[www].process() << ":" << label[www].label() << ":"
150  << label[www].instance() << " with " << kTrackAssociatorByChi2 << "\n";
151 
152  LogTrace("TrackValidator") << "Calling associateRecoToGen method"
153  << "\n";
154  recGenColl = trackGenAssociator->associateRecoToGen(trackCollection, TPCollectionHfake);
155  LogTrace("TrackValidator") << "Calling associateGenToReco method"
156  << "\n";
157  genRecColl = trackGenAssociator->associateGenToReco(trackCollection, TPCollectionHeff);
158  } else {
159  edm::LogVerbatim("TrackValidator") << "Analyzing " << label[www].process() << ":" << label[www].label() << ":"
160  << label[www].instance() << " with " << associators[0] << "\n";
161 
162  Handle<reco::GenToRecoCollection> gentorecoCollectionH;
163  event.getByToken(associatormapGtR, gentorecoCollectionH);
164  genRecColl = *(gentorecoCollectionH.product());
165 
166  Handle<reco::RecoToGenCollection> recotogenCollectionH;
167  event.getByToken(associatormapRtG, recotogenCollectionH);
168  recGenColl = *(recotogenCollectionH.product());
169  }
170 
171  // ########################################################
172  // fill simulation histograms (LOOP OVER TRACKINGPARTICLES)
173  // ########################################################
174 
175  //compute number of tracks per eta interval
176  //
177  edm::LogVerbatim("TrackValidator") << "\n# of GenParticles: " << tPCeff.size() << "\n";
178  int ats(0); //This counter counts the number of simTracks that are "associated" to recoTracks
179  int st(0); //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
180  unsigned sts(0); //This counter counts the number of simTracks surviving the bunchcrossing cut
181  unsigned asts(
182  0); //This counter counts the number of simTracks that are "associated" to recoTracks surviving the bunchcrossing cut
183  for (GenParticleCollection::size_type i = 0; i < tPCeff.size();
184  i++) { //loop over TPs collection for tracking efficiency
185  GenParticleRef tpr(TPCollectionHeff, i);
186  GenParticle* tp = const_cast<GenParticle*>(tpr.get());
187  TrackingParticle::Vector momentumTP;
188  TrackingParticle::Point vertexTP;
189  double dxyGen(0);
190  double dzGen(0);
191 
192  //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
193  //If the GenParticle is collison like, get the momentum and vertex at production state
195  //fixme this one shold be implemented
196  if (!gpSelector(*tp))
197  continue;
198  momentumTP = tp->momentum();
199  vertexTP = tp->vertex();
200  //Calcualte the impact parameters w.r.t. PCA
203  dxyGen = (-vertex.x() * sin(momentum.phi()) + vertex.y() * cos(momentum.phi()));
204  dzGen = vertex.z() - (vertex.x() * momentum.x() + vertex.y() * momentum.y()) / sqrt(momentum.perp2()) *
205  momentum.z() / sqrt(momentum.perp2());
206  }
207  //If the GenParticle is comics, get the momentum and vertex at PCA
208  else {
209  //if(! cosmictpSelector(*tp,&bs,event,setup)) continue;
210  momentumTP = parametersDefinerTP_->momentum(event, setup, *tp);
211  vertexTP = parametersDefinerTP_->vertex(event, setup, *tp);
212  dxyGen = (-vertexTP.x() * sin(momentumTP.phi()) + vertexTP.y() * cos(momentumTP.phi()));
213  dzGen = vertexTP.z() - (vertexTP.x() * momentumTP.x() + vertexTP.y() * momentumTP.y()) /
214  sqrt(momentumTP.perp2()) * momentumTP.z() / sqrt(momentumTP.perp2());
215  }
216  //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
217 
218  st++; //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
219 
220  // in the coming lines, histos are filled using as input
221  // - momentumTP
222  // - vertexTP
223  // - dxyGen
224  // - dzGen
225 
226  if (doSimPlots_ && w == 0) {
227  histoProducerAlgo_->fill_generic_simTrack_histos(histograms.histoProducerAlgo,
228  momentumTP,
229  vertexTP,
230  tp->collisionId()); //fixme: check meaning of collisionId
231  }
232  if (!doSimTrackPlots_)
233  continue;
234 
235  // ##############################################
236  // fill RecoAssociated GenTracks' histograms
237  // ##############################################
238  // bool isRecoMatched(false); // UNUSED
239  const reco::Track* matchedTrackPointer = nullptr;
240  std::vector<std::pair<RefToBase<Track>, double> > rt;
241  if (genRecColl.find(tpr) != genRecColl.end()) {
242  rt = (std::vector<std::pair<RefToBase<Track>, double> >)genRecColl[tpr];
243  if (!rt.empty()) {
244  ats++; //This counter counts the number of simTracks that have a recoTrack associated
245  // isRecoMatched = true; // UNUSED
246  matchedTrackPointer = rt.begin()->first.get();
247  edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt=" << sqrt(momentumTP.perp2())
248  << " associated with quality:" << rt.begin()->second << "\n";
249  }
250  } else {
251  edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt,eta,phi: " << sqrt(momentumTP.perp2())
252  << " , " << momentumTP.eta() << " , " << momentumTP.phi() << " , "
253  << " NOT associated to any reco::Track"
254  << "\n";
255  }
256 
257  int nSimHits = 0;
258  histoProducerAlgo_->fill_recoAssociated_simTrack_histos(histograms.histoProducerAlgo,
259  w,
260  *tp,
261  momentumTP,
262  vertexTP,
263  dxyGen,
264  dzGen,
265  nSimHits,
266  matchedTrackPointer,
267  puinfo.getPU_NumInteractions());
268 
269  sts++;
270  if (matchedTrackPointer)
271  asts++;
272 
273  } // End for (GenParticleCollection::size_type i=0; i<tPCeff.size(); i++){
274 
275  if (doSimPlots_ && w == 0) {
276  histoProducerAlgo_->fill_simTrackBased_histos(histograms.histoProducerAlgo, st);
277  }
278 
279  // ##############################################
280  // fill recoTracks histograms (LOOP OVER TRACKS)
281  // ##############################################
282  if (!doRecoTrackPlots_)
283  continue;
284  edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with " << label[www].process() << ":"
285  << label[www].label() << ":" << label[www].instance() << ": "
286  << trackCollection->size() << "\n";
287 
288  //int sat(0); //This counter counts the number of recoTracks that are associated to GenTracks from Signal only
289  int at(0); //This counter counts the number of recoTracks that are associated to GenTracks
290  int rT(0); //This counter counts the number of recoTracks in general
291 
292  for (View<Track>::size_type i = 0; i < trackCollection->size(); ++i) {
294  rT++;
295 
296  bool isSigGenMatched(false);
297  bool isGenMatched(false);
298  bool isChargeMatched(true);
299  int numAssocRecoTracks = 0;
300  int nSimHits = 0;
301  double sharedFraction = 0.;
302  std::vector<std::pair<GenParticleRef, double> > tp;
303  if (recGenColl.find(track) != recGenColl.end()) {
304  tp = recGenColl[track];
305  if (!tp.empty()) {
306  /*
307  std::vector<PSimHit> simhits=tp[0].first->trackPSimHit(DetId::Tracker);
308  nSimHits = simhits.end()-simhits.begin();
309  */
310  sharedFraction = tp[0].second;
311  isGenMatched = true;
312  if (tp[0].first->charge() != track->charge())
313  isChargeMatched = false;
314  if (genRecColl.find(tp[0].first) != genRecColl.end())
315  numAssocRecoTracks = genRecColl[tp[0].first].size();
316  //std::cout << numAssocRecoTracks << std::endl;
317  at++;
318  for (unsigned int tp_ite = 0; tp_ite < tp.size(); ++tp_ite) {
319  GenParticle trackpart = *(tp[tp_ite].first);
320  /*
321  if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)){
322  isSigGenMatched = true;
323  sat++;
324  break;
325  }
326  */
327  }
328  edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
329  << " associated with quality:" << tp.begin()->second << "\n";
330  }
331  } else {
332  edm::LogVerbatim("TrackValidator")
333  << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any GenParticle"
334  << "\n";
335  }
336 
337  double dR = 0; //fixme: plots vs dR and vs dRjet not implemented for now
338  histoProducerAlgo_->fill_generic_recoTrack_histos(histograms.histoProducerAlgo,
339  w,
340  *track,
341  ttopo,
342  bs.position(),
343  nullptr,
344  nullptr,
345  isGenMatched,
346  isSigGenMatched,
347  isChargeMatched,
348  numAssocRecoTracks,
349  puinfo.getPU_NumInteractions(),
350  nSimHits,
351  sharedFraction,
352  dR,
353  dR,
354  mvaDummy,
355  0,
356  0);
357 
358  // dE/dx
359  if (dodEdxPlots_)
360  histoProducerAlgo_->fill_dedx_recoTrack_histos(histograms.histoProducerAlgo, w, track, v_dEdx);
361 
362  //Fill other histos
363  //try{ //Is this really necessary ????
364 
365  if (tp.empty())
366  continue;
367 
368  histoProducerAlgo_->fill_simAssociated_recoTrack_histos(histograms.histoProducerAlgo, w, *track);
369 
370  GenParticleRef tpr = tp.begin()->first;
371 
372  /* TO BE FIXED LATER
373  if (associators[ww]=="TrackAssociatorByChi2"){
374  //association chi2
375  double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
376  h_assochi2[www]->Fill(assocChi2);
377  h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
378  }
379  else if (associators[ww]=="quickTrackAssociatorByHits"){
380  double fraction = tp.begin()->second;
381  h_assocFraction[www]->Fill(fraction);
382  h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
383  }
384  */
385 
386  //Get tracking particle parameters at point of closest approach to the beamline
387  TrackingParticle::Vector momentumTP = parametersDefinerTP_->momentum(event, setup, *(tpr.get()));
388  TrackingParticle::Point vertexTP = parametersDefinerTP_->vertex(event, setup, *(tpr.get()));
389  int chargeTP = tpr->charge();
390 
391  histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(
392  histograms.histoProducerAlgo, w, momentumTP, vertexTP, chargeTP, *track, bs.position());
393 
394  //TO BE FIXED
395  //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
396  //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
397 
398  /*
399  } // End of try{
400  catch (cms::Exception e){
401  LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
402  }
403  */
404 
405  } // End of for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
406 
407  histoProducerAlgo_->fill_trackBased_histos(histograms.histoProducerAlgo, w, at, rT, rT, st);
408 
409  edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
410  << "Total Associated (genToReco): " << ats << "\n"
411  << "Total Reconstructed: " << rT << "\n"
412  << "Total Associated (recoToGen): " << at << "\n"
413  << "Total Fakes: " << rT - at << "\n";
414 
415  w++;
416  } // End of for (unsigned int www=0;www<label.size();www++){
417 }
size
Write out results.
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > tTopoEsToken
Log< level::Info, true > LogVerbatim
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
unsigned int size_type
Definition: View.h:90
edm::EDGetTokenT< reco::GenToRecoCollection > associatormapGtR
edm::EDGetTokenT< edm::ValueMap< reco::DeDxData > > m_dEdx2Tag
static const std::string kTrackAssociatorByChi2("trackAssociatorByChi2")
reco::RecoToGenCollection associateRecoToGen(const edm::RefToBaseVector< reco::Track > &tracks, const edm::RefVector< reco::GenParticleCollection > &gens) const
Association Sim To Reco with Collections (Gen Particle version)
T w() const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
T const * product() const
Definition: Handle.h:70
constexpr bool isUninitialized() const noexcept
Definition: EDGetToken.h:104
void dqmAnalyze(const edm::Event &, const edm::EventSetup &, const Histograms &) const override
Method called once per event.
uint16_t size_type
edm::EDGetTokenT< reco::BeamSpot > bsSrc
const_iterator find(const key_type &k) const
find element with specified reference key
MultiTrackValidatorGenPs(const edm::ParameterSet &pset)
Constructor.
#define LogTrace(id)
const_iterator end() const
last iterator over the map (read only)
math::XYZPointD Point
point in the space
std::vector< edm::InputTag > associators
edm::EDGetTokenT< reco::TrackToGenParticleAssociator > label_gen_associator
std::vector< edm::EDGetTokenT< edm::View< reco::Track > > > labelToken
edm::EDGetTokenT< edm::ValueMap< reco::DeDxData > > m_dEdx1Tag
GenParticleCustomSelector gpSelector
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > label_pileupinfo
std::vector< edm::InputTag > label
SingleObjectSelector< GenParticleCollection, ::GenParticleCustomSelector > GenParticleCustomSelector
std::unique_ptr< ParametersDefinerForTP > parametersDefinerTP_
trackCollection
Definition: JetHT_cfg.py:51
std::unique_ptr< MTVHistoProducerAlgoForTracker > histoProducerAlgo_
Log< level::Info, false > LogInfo
~MultiTrackValidatorGenPs() override
Destructor.
const bool ignoremissingtkcollection_
fixed size matrix
HLT enums.
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
edm::EDGetTokenT< TrackingParticleCollection > label_tp_fake
math::XYZVectorD Vector
point in the space
edm::EDGetTokenT< reco::RecoToGenCollection > associatormapRtG
edm::EDGetTokenT< TrackingParticleCollection > label_tp_effic
const int getPU_NumInteractions() const
reco::GenToRecoCollection associateGenToReco(const edm::RefToBaseVector< reco::Track > &tracks, const edm::RefVector< reco::GenParticleCollection > &gens) const
Association Sim To Reco with Collections (Gen Particle version)
Definition: event.py:1
const bool parametersDefinerIsCosmic_