CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PATPackedCandidateProducer.cc
Go to the documentation of this file.
1 #include <string>
2 
3 
22 
23 /*#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
24 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
25 #include "MagneticField/Engine/interface/MagneticField.h"
26 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
27 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
28 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
29 */
30 //#define CRAZYSORT
31 
32 namespace pat {
34  const static int qualityMap[8] = {1,0,1,1,4,4,5,6};
35 
37  public:
40 
41  virtual void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
42 
43  //sorting of cands to maximize the zlib compression
45  if (std::abs(i.charge()) == std::abs(j.charge())) {
46  if(i.charge()!=0){
47  if(i.pt() > minPtForTrackProperties_ and j.pt() <= minPtForTrackProperties_ ) return true;
48  if(i.pt() <= minPtForTrackProperties_ and j.pt() > minPtForTrackProperties_ ) return false;
49  }
50  if(i.vertexRef() == j.vertexRef())
51  return i.eta() > j.eta();
52  else
53  return i.vertexRef().key() < j.vertexRef().key();
54  }
55  return std::abs(i.charge()) > std::abs(j.charge());
56  }
57  template <typename T>
58  std::vector<size_t> sort_indexes(const std::vector<T> &v ) const {
59  std::vector<size_t> idx(v.size());
60  for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;
61  std::sort(idx.begin(), idx.end(),[&v,this](size_t i1, size_t i2) { return candsOrdering(v[i1],v[i2]);});
62  return idx;
63  }
64 
65  private:
77  std::vector< edm::EDGetTokenT<edm::View<reco::CompositePtrCandidate> > > SVWhiteLists_;
78 
80  // for debugging
81  float calcDxy(float dx, float dy, float phi) const {
82  return - dx * std::sin(phi) + dy * std::cos(phi);
83  }
85  return p.Z()-v.Z() - ((p.X()-v.X()) * c.px() + (p.Y()-v.Y())*c.py()) * c.pz()/(c.pt()*c.pt());
86  }
87  };
88 }
89 
91  Cands_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("inputCollection"))),
92  PVs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("inputVertices"))),
93  PVAsso_(consumes<edm::Association<reco::VertexCollection> >(iConfig.getParameter<edm::InputTag>("vertexAssociator"))),
94  PVAssoQuality_(consumes<edm::ValueMap<int> >(iConfig.getParameter<edm::InputTag>("vertexAssociator"))),
95  PVOrigs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("originalVertices"))),
96  TKOrigs_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("originalTracks"))),
97  PuppiWeight_(consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("PuppiSrc"))),
98  PuppiWeightNoLep_(consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("PuppiNoLepSrc"))),
99  PuppiCandsMap_(consumes<edm::ValueMap<reco::CandidatePtr> >(iConfig.getParameter<edm::InputTag>("PuppiSrc"))),
100  PuppiCands_(consumes<std::vector< reco::PFCandidate > >(iConfig.getParameter<edm::InputTag>("PuppiSrc"))),
101  PuppiCandsNoLep_(consumes<std::vector< reco::PFCandidate > >(iConfig.getParameter<edm::InputTag>("PuppiNoLepSrc"))),
102  minPtForTrackProperties_(iConfig.getParameter<double>("minPtForTrackProperties"))
103 {
104  std::vector<edm::InputTag> sv_tags = iConfig.getParameter<std::vector<edm::InputTag> >("secondaryVerticesForWhiteList");
105  for(auto itag : sv_tags){
106  SVWhiteLists_.push_back(
108  );
109  }
110 
111  produces< std::vector<pat::PackedCandidate> > ();
112  produces< edm::Association<pat::PackedCandidateCollection> > ();
113  produces< edm::Association<reco::PFCandidateCollection> > ();
114 }
115 
117 
118 
119 
121 
123  iEvent.getByToken( Cands_, cands );
124  std::vector<reco::Candidate>::const_iterator cand;
125 
126  edm::Handle< edm::ValueMap<float> > puppiWeight;
127  iEvent.getByToken( PuppiWeight_, puppiWeight );
129  iEvent.getByToken( PuppiCandsMap_, puppiCandsMap );
131  iEvent.getByToken( PuppiCands_, puppiCands );
132  std::vector<int> mappingPuppi(puppiCands->size());
133 
134  edm::Handle< edm::ValueMap<float> > puppiWeightNoLep;
135  iEvent.getByToken( PuppiWeightNoLep_, puppiWeightNoLep );
137  iEvent.getByToken( PuppiCandsNoLep_, puppiCandsNoLep );
138 
139  std::vector<reco::CandidatePtr> puppiCandsNoLepPtrs;
140  if (puppiCandsNoLep.isValid()){
141  for (auto pup : *puppiCandsNoLep){
142  puppiCandsNoLepPtrs.push_back(pup.sourceCandidatePtr(0));
143  }
144  }
145  auto const& puppiCandsNoLepV = puppiCandsNoLep.product();
146 
148  iEvent.getByToken( PVOrigs_, PVOrigs );
149 
151  iEvent.getByToken(PVAsso_,assoHandle);
152  edm::Handle<edm::ValueMap<int> > assoQualityHandle;
153  iEvent.getByToken(PVAssoQuality_,assoQualityHandle);
154  const edm::Association<reco::VertexCollection> & associatedPV=*(assoHandle.product());
155  const edm::ValueMap<int> & associationQuality=*(assoQualityHandle.product());
156 
157 
158  std::set<unsigned int> whiteList;
159  for(auto itoken : SVWhiteLists_) {
161  iEvent.getByToken(itoken, svWhiteListHandle);
162  const edm::View<reco::CompositePtrCandidate > & svWhiteList=*(svWhiteListHandle.product());
163  for(unsigned int i=0; i<svWhiteList.size();i++) {
164  for(unsigned int j=0; j< svWhiteList[i].numberOfSourceCandidatePtrs(); j++) {
165  const edm::Ptr<reco::Candidate> & c = svWhiteList[i].sourceCandidatePtr(j);
166  if(c.id() == cands.id()) whiteList.insert(c.key());
167  }
168  }
169  }
170 
171 
173  iEvent.getByToken( PVs_, PVs );
174  reco::VertexRef PV(PVs.id());
175  reco::VertexRefProd PVRefProd(PVs);
176  math::XYZPoint PVpos;
177 
178 
180  iEvent.getByToken( TKOrigs_, TKOrigs );
181  std::auto_ptr< std::vector<pat::PackedCandidate> > outPtrP( new std::vector<pat::PackedCandidate> );
182  std::vector<int> mapping(cands->size());
183  std::vector<int> mappingReverse(cands->size());
184  std::vector<int> mappingTk(TKOrigs->size(), -1);
185 
186  for(unsigned int ic=0, nc = cands->size(); ic < nc; ++ic) {
187  const reco::PFCandidate &cand=(*cands)[ic];
188  float phiAtVtx = cand.phi();
189  const reco::Track *ctrack = 0;
190  if ((abs(cand.pdgId()) == 11 || cand.pdgId() == 22) && cand.gsfTrackRef().isNonnull()) {
191  ctrack = &*cand.gsfTrackRef();
192  } else if (cand.trackRef().isNonnull()) {
193  ctrack = &*cand.trackRef();
194  }
195  if (ctrack) {
196  float dist=1e99;
197  int pvi=-1;
198  for(size_t ii=0;ii<PVs->size();ii++){
199  float dz=std::abs(ctrack->dz( ((*PVs)[ii]).position()));
200  if(dz<dist) {pvi=ii;dist=dz; }
201  }
202  PV = reco::VertexRef(PVs, pvi);
203  math::XYZPoint vtx = cand.vertex();
205  const reco::VertexRef & PVOrig = associatedPV[reco::CandidatePtr(cands,ic)];
206  if(PVOrig.isNonnull()) PV = reco::VertexRef(PVs, PVOrig.key()); // WARNING: assume the PV slimmer is keeping same order
207  int quality=associationQuality[reco::CandidatePtr(cands,ic)];
208 // if ((size_t)pvi!=PVOrig.key()) std::cout << "not closest in Z" << pvi << " " << PVOrig.key() << " " << cand.pt() << " " << quality << std::endl;
209  // TrajectoryStateOnSurface tsos = extrapolator.extrapolate(trajectoryStateTransform::initialFreeState(*ctrack,&*magneticField), RecoVertex::convertPos(PV->position()));
210  // vtx = tsos.globalPosition();
211  // phiAtVtx = tsos.globalDirection().phi();
212  vtx = ctrack->referencePoint();
213  phiAtVtx = ctrack->phi();
215  if (nlost == 0) {
216  if (ctrack->hitPattern().hasValidHitInFirstPixelBarrel()) {
218  }
219  } else {
221  }
222 
223 
224  outPtrP->push_back( pat::PackedCandidate(cand.polarP4(), vtx, phiAtVtx, cand.pdgId(), PVRefProd, PV.key()));
225  outPtrP->back().setAssociationQuality(pat::PackedCandidate::PVAssociationQuality(qualityMap[quality]));
226  if(cand.trackRef().isNonnull() && PVOrig->trackWeight(cand.trackRef()) > 0.5 && quality == 7) {
227  outPtrP->back().setAssociationQuality(pat::PackedCandidate::UsedInFitTight);
228  }
229  // properties of the best track
230  outPtrP->back().setLostInnerHits( lostHits );
231  if(outPtrP->back().pt() > minPtForTrackProperties_ || whiteList.find(ic)!=whiteList.end()) {
232  outPtrP->back().setTrackProperties(*ctrack);
233  //outPtrP->back().setTrackProperties(*ctrack,tsos.curvilinearError());
234  }
235 
236  // these things are always for the CKF track
237  outPtrP->back().setTrackHighPurity( cand.trackRef().isNonnull() && cand.trackRef()->quality(reco::Track::highPurity) );
238  if (cand.muonRef().isNonnull()) {
239  outPtrP->back().setMuonID(cand.muonRef()->isStandAloneMuon(), cand.muonRef()->isGlobalMuon());
240  }
241  } else {
242 
243  if (!PVs->empty()) {
244  PV = reco::VertexRef(PVs, 0);
245  PVpos = PV->position();
246  }
247 
248  outPtrP->push_back( pat::PackedCandidate(cand.polarP4(), PVpos, cand.phi(), cand.pdgId(), PVRefProd, PV.key()));
250  }
251 
252  // neutrals
253 
254  if(abs(cand.pdgId()) == 1 || abs(cand.pdgId()) == 130) {
255  outPtrP->back().setHcalFraction(cand.hcalEnergy()/(cand.ecalEnergy()+cand.hcalEnergy()));
256  } else {
257  outPtrP->back().setHcalFraction(0);
258  }
259 
260 
261  if (puppiWeight.isValid()){
262  reco::PFCandidateRef pkref( cands, ic );
263  // outPtrP->back().setPuppiWeight( (*puppiWeight)[pkref]);
264 
265  float puppiWeightVal = (*puppiWeight)[pkref];
266  float puppiWeightNoLepVal = 0.0;
267 
268  // Check the "no lepton" puppi weights.
269  // If present, then it is not a lepton, use stored weight
270  // If absent, it is a lepton, so set the weight to 1.0
271  if ( puppiWeightNoLep.isValid() ) {
272  // Look for the pointer inside the "no lepton" candidate collection.
273  auto pkrefPtr = pkref->sourceCandidatePtr(0);
274 
275  bool foundNoLep = false;
276  for ( size_t ipcnl = 0; ipcnl < puppiCandsNoLepPtrs.size(); ipcnl++){
277  if (puppiCandsNoLepPtrs[ipcnl] == pkrefPtr){
278  foundNoLep = true;
279  puppiWeightNoLepVal = puppiCandsNoLepV->at(ipcnl).pt()/cand.pt(); // a hack for now, should use the value map
280  break;
281  }
282  }
283  if ( !foundNoLep || puppiWeightNoLepVal > 1 ) {
284  puppiWeightNoLepVal = 1.0;
285  }
286  }
287  outPtrP->back().setPuppiWeight( puppiWeightVal, puppiWeightNoLepVal );
288 
289  mappingPuppi[((*puppiCandsMap)[pkref]).key()]=ic;
290  }
291 
292  mapping[ic] = ic; // trivial at the moment!
293  if (cand.trackRef().isNonnull() && cand.trackRef().id() == TKOrigs.id()) {
294  mappingTk[cand.trackRef().key()] = ic;
295  }
296 
297  }
298 
299  std::auto_ptr< std::vector<pat::PackedCandidate> > outPtrPSorted( new std::vector<pat::PackedCandidate> );
300  std::vector<size_t> order=sort_indexes(*outPtrP);
301  std::vector<size_t> reverseOrder(order.size());
302  for(size_t i=0,nc=cands->size();i<nc;i++) {
303  outPtrPSorted->push_back((*outPtrP)[order[i]]);
304  reverseOrder[order[i]] = i;
305  mappingReverse[order[i]]=i;
306  }
307 
308  // Fix track association for sorted candidates
309  for(size_t i=0,ntk=mappingTk.size();i<ntk;i++){
310  if(mappingTk[i] >= 0)
311  mappingTk[i]=reverseOrder[mappingTk[i]];
312  }
313 
314  for(size_t i=0,ntk=mappingPuppi.size();i<ntk;i++){
315  mappingPuppi[i]=reverseOrder[mappingPuppi[i]];
316  }
317 
318  edm::OrphanHandle<pat::PackedCandidateCollection> oh = iEvent.put( outPtrPSorted );
319 
320  // now build the two maps
321  std::auto_ptr<edm::Association<pat::PackedCandidateCollection> > pf2pc(new edm::Association<pat::PackedCandidateCollection>(oh ));
322  std::auto_ptr<edm::Association<reco::PFCandidateCollection > > pc2pf(new edm::Association<reco::PFCandidateCollection >(cands));
325  pf2pcFiller.insert(cands, mappingReverse.begin(), mappingReverse.end());
326  pc2pfFiller.insert(oh , order.begin(), order.end());
327  // include also the mapping track -> packed PFCand
328  pf2pcFiller.insert(TKOrigs, mappingTk.begin(), mappingTk.end());
329  pf2pcFiller.insert(puppiCands, mappingPuppi.begin(), mappingPuppi.end());
330 
331  pf2pcFiller.fill();
332  pc2pfFiller.fill();
333  iEvent.put(pf2pc);
334  iEvent.put(pc2pf);
335 
336 }
337 
338 
T getParameter(std::string const &) const
double ecalEnergy() const
return corrected Ecal energy
Definition: PFCandidate.h:221
const Point & referencePoint() const
Reference point on the track.
Definition: TrackBase.h:676
int i
Definition: DBlmapReader.cc:9
const edm::EDGetTokenT< edm::ValueMap< float > > PuppiWeight_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
key_type key() const
Definition: Ptr.h:186
ProductID id() const
Definition: HandleBase.cc:15
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
virtual double pt() const =0
transverse momentum
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
virtual double pz() const =0
z coordinate of momentum vector
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:52
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
virtual double phi() const final
momentum azimuthal angle
size_type size() const
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:640
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
static const int qualityMap[8]
conversion map from quality flags used in PV association and miniAOD one
const edm::EDGetTokenT< reco::PFCandidateCollection > Cands_
virtual double pt() const
transverse momentum
key_type key() const
Accessor for product key.
Definition: Ref.h:264
int ii
Definition: cuy.py:588
virtual double eta() const
momentum pseudorapidity
const reco::VertexRef vertexRef() const
ProductID id() const
Accessor for product ID.
Definition: Ref.h:258
const edm::EDGetTokenT< edm::ValueMap< int > > PVAssoQuality_
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:433
int iEvent
Definition: GenABIO.cc:230
PATPackedCandidateProducer(const edm::ParameterSet &)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
float calcDxy(float dx, float dy, float phi) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
virtual int charge() const
electric charge
virtual const Point & vertex() const
vertex position (overwritten by PF...)
Definition: PFCandidate.cc:647
const edm::EDGetTokenT< reco::TrackCollection > TKOrigs_
const edm::EDGetTokenT< std::vector< reco::PFCandidate > > PuppiCands_
std::vector< size_t > sort_indexes(const std::vector< T > &v) const
virtual double py() const =0
y coordinate of momentum vector
const edm::EDGetTokenT< edm::ValueMap< reco::CandidatePtr > > PuppiCandsMap_
bool isValid() const
Definition: HandleBase.h:75
const edm::EDGetTokenT< edm::ValueMap< float > > PuppiWeightNoLep_
edm::Ref< VertexCollection > VertexRef
persistent reference to a Vertex
Definition: VertexFwd.h:13
string key
FastSim: produces sample of signal events, overlayed with premixed minbias events.
virtual void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
reco::MuonRef muonRef() const
Definition: PFCandidate.cc:450
virtual double px() const =0
x coordinate of momentum vector
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:604
LostInnerHits
Enumerator specifying the.
std::vector< edm::EDGetTokenT< edm::View< reco::CompositePtrCandidate > > > SVWhiteLists_
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
float calcDz(reco::Candidate::Point p, reco::Candidate::Point v, const reco::Candidate &c) const
edm::Ptr< Candidate > CandidatePtr
persistent reference to an object in a collection of Candidate objects
Definition: CandidateFwd.h:25
T const * product() const
Definition: Handle.h:81
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:445
const edm::EDGetTokenT< std::vector< reco::PFCandidate > > PuppiCandsNoLep_
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
ProductID id() const
Accessor for product ID.
Definition: Ptr.h:181
const edm::EDGetTokenT< reco::VertexCollection > PVOrigs_
virtual const PolarLorentzVector & polarP4() const final
four-momentum Lorentz vector
int numberOfLostHits(HitCategory category) const
Definition: HitPattern.h:902
virtual int pdgId() const final
PDG identifier.
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:39
reco::GsfTrackRef gsfTrackRef() const
Definition: PFCandidate.cc:471
double hcalEnergy() const
return corrected Hcal energy
Definition: PFCandidate.h:231
math::XYZPoint Point
point in the space
Definition: Candidate.h:41
const edm::EDGetTokenT< edm::Association< reco::VertexCollection > > PVAsso_
bool hasValidHitInFirstPixelBarrel() const
Definition: HitPattern.cc:326
const edm::EDGetTokenT< reco::VertexCollection > PVs_
bool candsOrdering(pat::PackedCandidate i, pat::PackedCandidate j) const
virtual double pt() const final
transverse momentum