CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SoftConversionProducer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <memory>
4 
5 // Framework
10 //
15 //
24 
30 
33 
34 #include "Math/GenVector/VectorUtil.h"
35 
36 
38 
39  LogDebug("SoftConversionProducer") << " SoftConversionProducer CTOR " << "\n";
40 
41  conversionOITrackProducer_ = conf_.getParameter<std::string>("conversionOITrackProducer");
42  conversionIOTrackProducer_ = conf_.getParameter<std::string>("conversionIOTrackProducer");
43  outInTrackClusterAssociationCollection_ = conf_.getParameter<std::string>("outInTrackClusterAssociationCollection");
44  inOutTrackClusterAssociationCollection_ = conf_.getParameter<std::string>("inOutTrackClusterAssociationCollection");
45  clusterType_ = conf_.getParameter<std::string>("clusterType");
46  clusterBarrelCollection_ = conf_.getParameter<edm::InputTag>("clusterBarrelCollection");
47  clusterEndcapCollection_ = conf_.getParameter<edm::InputTag>("clusterEndcapCollection");
48  softConversionCollection_ = conf_.getParameter<std::string>("softConversionCollection");
49  trackMaxChi2_ = conf_.getParameter<double>("trackMaxChi2");
50  trackMinHits_ = conf_.getParameter<double>("trackMinHits");
51  clustersMaxDeltaEta_ = conf_.getParameter<double>("clustersMaxDeltaEta");
52  clustersMaxDeltaPhi_ = conf_.getParameter<double>("clustersMaxDeltaPhi");
53 
55  theVertexFinder_ = 0;
57 
58  // Register the product
59  produces< reco::ConversionCollection >(softConversionCollection_);
60 
61 }
62 
63 
65 
66 
68 
69  //get magnetic field
70  theEventSetup.get<IdealMagneticFieldRecord>().get(theMF_);
71 
72  // instantiate the Track Pair Finder algorithm
74 
75  // instantiate the Vertex Finder algorithm
77 
78  // instantiate the algorithm for finding the position of the track extrapolation at the Ecal front face
80 
81 }
82 
83 
84 void SoftConversionProducer::endRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
85 
89 
90 }
91 
92 
93 
94 void SoftConversionProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
95 
97  theEvent.getByLabel(conversionOITrackProducer_, outInTrkHandle);
98 
99  edm::Handle<reco::TrackCaloClusterPtrAssociation> outInTrkClusterAssocHandle;
100  theEvent.getByLabel( conversionOITrackProducer_, outInTrackClusterAssociationCollection_, outInTrkClusterAssocHandle);
101 
102  edm::Handle<reco::TrackCollection> inOutTrkHandle;
103  theEvent.getByLabel(conversionIOTrackProducer_, inOutTrkHandle);
104 
105  edm::Handle<reco::TrackCaloClusterPtrAssociation> inOutTrkClusterAssocHandle;
106  theEvent.getByLabel( conversionIOTrackProducer_, inOutTrackClusterAssociationCollection_, inOutTrkClusterAssocHandle);
107 
108  edm::Handle<edm::View<reco::CaloCluster> > clusterBarrelHandle;
109  theEvent.getByLabel(clusterBarrelCollection_, clusterBarrelHandle);
110 
111  edm::Handle<edm::View<reco::CaloCluster> > clusterEndcapHandle;
112  if(clusterType_ == "BasicCluster"){
113  //theEvent.getByLabel(clusterEndcapProducer_, clusterEndcapCollection_, clusterEndcapHandle);
114  theEvent.getByLabel(clusterEndcapCollection_, clusterEndcapHandle);
115  }
116 
117  // create temporary map to loop over tracks conveniently
118  TrackClusterMap trackClusterMap;
119 
120  int nTracksOI = (int) outInTrkHandle->size();
121  // std::cout << " nTracksOI " << nTracksOI << std::endl;
122  for(int itrk=0; itrk<nTracksOI; itrk++){
123  reco::TrackRef tRef(outInTrkHandle,itrk);
124  if(!trackQualityCut(tRef)) continue;
125  reco::CaloClusterPtr cRef = (*outInTrkClusterAssocHandle)[tRef];
126  trackClusterMap.push_back(make_pair(tRef,cRef));
127  }
128 
129  int nTracksIO = (int) inOutTrkHandle->size();
130  // std::cout << " nTracksIO " << nTracksIO << std::endl;
131  for(int itrk=0; itrk<nTracksIO; itrk++){
132  reco::TrackRef tRef(inOutTrkHandle,itrk);
133  if(!trackQualityCut(tRef)) continue;
134  reco::CaloClusterPtr cRef = (*inOutTrkClusterAssocHandle)[tRef];
135  trackClusterMap.push_back(make_pair(tRef,cRef));
136  }
137 
138  // Transform Track into TransientTrack (needed by the Vertex fitter)
139  edm::ESHandle<TransientTrackBuilder> theTransientTrackBuilder;
140  theEventSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theTransientTrackBuilder);
141 
142  // the output collection to be produced from this producer
143  std::auto_ptr< reco::ConversionCollection > outputColl(new reco::ConversionCollection);
144 
145  // prepare iterator
146  TrackClusterMap::iterator iter1 = trackClusterMap.begin();
147  TrackClusterMap::iterator iter2 = trackClusterMap.begin();
148  TrackClusterMap::iterator iter_end = trackClusterMap.end();
149 
150 
151 
152  // double-loop to make pairs
153  for( ; iter1 != iter_end; iter1++) {
154  // std::cout << " back to start of loop 1 " << std::endl;
155  const reco::TrackRef trk1 = iter1->first;
156 
157  for(iter2 = iter1+1; iter2 != iter_end; iter2++) {
158  // std::cout << " back to start of loop 2 " << std::endl;
159  const reco::TrackRef trk2 = iter2->first;
160  if(trk1 == trk2) continue;
161 
162 
163  const reco::CaloClusterPtr cls1 = iter1->second;
164  reco::TransientTrack tsk1 = theTransientTrackBuilder->build(*trk1);
165 
166  const reco::CaloClusterPtr cls2 = iter2->second;
167  reco::TransientTrack tsk2 = theTransientTrackBuilder->build(*trk2);
168 
169  // std::cout << " after transient " << std::endl;
170  // std::cout << " eta1 " << cls1->position().Eta() << " eta2 " << cls2->position().Eta() << std::endl;
171 
174 
175 
176  double dEta = std::abs(cls1->position().Eta() - cls2->position().Eta());
177  if(dEta > clustersMaxDeltaEta_) continue;
178  double dPhi = std::abs(ROOT::Math::VectorUtil::DeltaPhi(cls1->position(),cls2->position()));
179  if(dPhi > clustersMaxDeltaPhi_) continue;
180 
181  std::vector<reco::TransientTrack> toBeFitted;
182  toBeFitted.push_back(tsk1);
183  toBeFitted.push_back(tsk2);
184 
185  //
187  md.calculate ( tsk1.initialFreeState(), tsk2.initialFreeState() );
188  float minAppDist = md.distance();
189 
190 
191 
192  // std::cout << " Before vertex " << std::endl;
193  reco::Vertex theConversionVertex = (reco::Vertex) theVertexFinder_->run(toBeFitted);
194  //std::cout << " After vertex " << std::endl;
195 
196  if(theConversionVertex.isValid()){
197  // std::cout << " valid vertex " << std::endl;
199  scRefs.push_back(cls1);
200  scRefs.push_back(cls2);
201  //std::cout << " after push back scref " << std::endl;
202  std::vector<reco::CaloClusterPtr> clusterRefs;
203  clusterRefs.push_back(cls1);
204  clusterRefs.push_back(cls2);
205  //std::cout << " after push back slusterref " << std::endl;
206  std::vector<reco::TrackRef> trkRefs;
207  trkRefs.push_back(trk1);
208  trkRefs.push_back(trk2);
209 
210  std::vector<math::XYZVector> trackPin;
211  std::vector<math::XYZVector> trackPout;
212  trackPin.push_back( trk1->innerMomentum());
213  trackPin.push_back( trk2->innerMomentum());
214  trackPout.push_back( trk1->outerMomentum());
215  trackPout.push_back( trk2->outerMomentum());
216 
217  // std::cout << " Before impact finder " << std::endl;
218  std::vector<math::XYZPoint> trkPositionAtEcal = theEcalImpactPositionFinder_->find(toBeFitted,clusterBarrelHandle);
219  // std::cout << " After impact finder " << std::endl;
220  if((clusterType_ == "BasicCluster") && (std::abs(cls2->position().Eta()) > 1.5)){
221  trkPositionAtEcal.clear();
222  trkPositionAtEcal = theEcalImpactPositionFinder_->find(toBeFitted,clusterEndcapHandle);
223  }
224 
225  double dummy = -9999.;
226  std::vector<math::XYZPoint> dummyVec;
227  reco::Conversion newCandidate(scRefs, trkRefs, trkPositionAtEcal, theConversionVertex, clusterRefs, minAppDist,dummyVec, trackPin, trackPout, dummy );
228 
229  // Check this candidate is already in the collection.
230  // This is checking that two tracks in a conversion candidate are identicial.
231 
232  if(NotAlreadyIn(newCandidate,outputColl)) outputColl->push_back(newCandidate);
233 
234  // printf("=====> run(%d), event(%d) <=====\n",theEvent.id().run(),theEvent.id().event());
235  // printf("Found a softConverion with vtxR(%f), vtxEta(%f), pt(%f), pt1(%f), pt2(%f)\n",
236  // newCandidate.conversionVertex().position().rho(),newCandidate.conversionVertex().position().eta(),
237  // newCandidate.pairMomentum().perp(),trk1->momentum().rho(),trk2->momentum().rho());
238 
239  clusterRefs.clear();
240  trkRefs.clear();
241  trkPositionAtEcal.clear();
242  }// if(theConversionVertex.isValid()
243 
244  toBeFitted.clear();
245 
246  }// end of iter2
247  }// end of iter1
248 
249  // std::cout << " Putting stuff in the event " << std::endl;
250  // put the collection into the event
251  theEvent.put( outputColl, softConversionCollection_);
252 
253 }
254 
255 
257  return (trk->numberOfValidHits() >= trackMinHits_ && trk->normalizedChi2() < trackMaxChi2_);
258 }
259 
260 
262  const std::auto_ptr<reco::ConversionCollection>& outputColl) const {
263 
264  if(outputColl->size() == 0) return true;
265 
266  reco::ConversionCollection::const_iterator it = outputColl->begin();
267  reco::ConversionCollection::const_iterator it_end = outputColl->end();
268  for( ; it != it_end; it++){
269  const reco::Conversion& conv = *it;
270  if((thisConv.tracks()[0] == conv.tracks()[0] && thisConv.tracks()[1] == conv.tracks()[1]) ||
271  (thisConv.tracks()[0] == conv.tracks()[1] && thisConv.tracks()[1] == conv.tracks()[0]))
272  return false;
273  }// for
274 
275  return true;
276 }
#define LogDebug(id)
virtual float distance() const
virtual void beginRun(edm::Run &r, edm::EventSetup const &es)
T getParameter(std::string const &) const
virtual bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb)
ConversionTrackPairFinder * theTrackPairFinder_
std::string inOutTrackClusterAssociationCollection_
edm::InputTag clusterEndcapCollection_
edm::ESHandle< MagneticField > theMF_
static HepMC::IO_HEPEVT conv
ConversionTrackEcalImpactPoint * theEcalImpactPositionFinder_
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:135
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:61
#define abs(x)
Definition: mlp_lapack.h:159
ConversionVertexFinder * theVertexFinder_
std::vector< Conversion > ConversionCollection
collectin of Conversion objects
Definition: ConversionFwd.h:9
std::string outInTrackClusterAssociationCollection_
TrajectoryStateOnSurface innermostMeasurementState() const
std::vector< edm::RefToBase< reco::Track > > tracks() const
vector of track to base references
Definition: Conversion.cc:180
std::vector< math::XYZPoint > find(const std::vector< reco::TransientTrack > &tracks, const edm::Handle< edm::View< reco::CaloCluster > > &bcHandle)
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
TransientVertex run(std::vector< reco::TransientTrack > pair)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
edm::InputTag clusterBarrelCollection_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
virtual void produce(edm::Event &evt, const edm::EventSetup &es)
FreeTrajectoryState initialFreeState() const
const GlobalTrajectoryParameters & globalParameters() const
bool NotAlreadyIn(const reco::Conversion &thisConv, const std::auto_ptr< reco::ConversionCollection > &outputColl) const
const T & get() const
Definition: EventSetup.h:55
SoftConversionProducer(const edm::ParameterSet &ps)
std::vector< std::pair< reco::TrackRef, reco::CaloClusterPtr > > TrackClusterMap
tuple config
Definition: cmsDriver.py:17
bool trackQualityCut(const reco::TrackRef &trk)
virtual void endRun(edm::Run &, edm::EventSetup const &)
Definition: Run.h:31