00001 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
00002 #include "DataFormats/TrackReco/interface/Track.h"
00003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00004 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00005
00006 using namespace reco;
00007
00008
00009 Conversion::Conversion( const reco::CaloClusterPtrVector sc,
00010 const std::vector<reco::TrackRef>& tr,
00011 const std::vector<math::XYZPointF>& trackPositionAtEcal,
00012 const reco::Vertex & convVtx,
00013 const std::vector<reco::CaloClusterPtr> & matchingBC,
00014 const float DCA,
00015 const std::vector<math::XYZPointF> & innPoint,
00016 const std::vector<math::XYZVectorF> & trackPin,
00017 const std::vector<math::XYZVectorF> & trackPout,
00018 const float mva,
00019 ConversionAlgorithm algo):
00020
00021
00022 caloCluster_(sc), tracks_(tr),
00023 thePositionAtEcal_(trackPositionAtEcal),
00024 theConversionVertex_(convVtx),
00025 theMatchingBCs_(matchingBC),
00026 theMinDistOfApproach_(DCA),
00027 theTrackInnerPosition_(innPoint),
00028 theTrackPin_(trackPin),
00029 theTrackPout_(trackPout),
00030 nSharedHits_(0),
00031 theMVAout_(mva),
00032 algorithm_(algo),
00033 qualityMask_(0)
00034 {
00035
00036 }
00037
00038
00039
00040
00041 Conversion::Conversion( const reco::CaloClusterPtrVector sc,
00042 const std::vector<edm::RefToBase<reco::Track> >& tr,
00043 const std::vector<math::XYZPointF>& trackPositionAtEcal,
00044 const reco::Vertex & convVtx,
00045 const std::vector<reco::CaloClusterPtr> & matchingBC,
00046 const float DCA,
00047 const std::vector<math::XYZPointF> & innPoint,
00048 const std::vector<math::XYZVectorF> & trackPin,
00049 const std::vector<math::XYZVectorF> & trackPout,
00050 const std::vector<uint8_t>& nHitsBeforeVtx,
00051 const std::vector<Measurement1DFloat> & dlClosestHitToVtx,
00052 uint8_t nSharedHits,
00053 const float mva,
00054 ConversionAlgorithm algo):
00055
00056
00057 caloCluster_(sc), trackToBaseRefs_(tr),
00058 thePositionAtEcal_(trackPositionAtEcal),
00059 theConversionVertex_(convVtx),
00060 theMatchingBCs_(matchingBC),
00061 theMinDistOfApproach_(DCA),
00062 theTrackInnerPosition_(innPoint),
00063 theTrackPin_(trackPin),
00064 theTrackPout_(trackPout),
00065 nHitsBeforeVtx_(nHitsBeforeVtx),
00066 dlClosestHitToVtx_(dlClosestHitToVtx),
00067 nSharedHits_(nSharedHits),
00068 theMVAout_(mva),
00069 algorithm_(algo),
00070 qualityMask_(0)
00071 {
00072
00073 }
00074
00075
00076
00077
00078 Conversion::Conversion( const reco::CaloClusterPtrVector sc,
00079 const std::vector<reco::TrackRef>& tr,
00080 const reco::Vertex & convVtx,
00081 ConversionAlgorithm algo):
00082 caloCluster_(sc), tracks_(tr),
00083 theConversionVertex_(convVtx),
00084 nSharedHits_(0),
00085 algorithm_(algo),
00086 qualityMask_(0)
00087 {
00088
00089
00090 theMinDistOfApproach_ = 9999.;
00091 theMVAout_ = 9999.;
00092 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00093 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00094 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00095 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00096 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00097 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00098 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00099 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00100
00101
00102 }
00103
00104
00105 Conversion::Conversion( const reco::CaloClusterPtrVector sc,
00106 const std::vector<edm::RefToBase<reco::Track> >& tr,
00107 const reco::Vertex & convVtx,
00108 ConversionAlgorithm algo):
00109 caloCluster_(sc), trackToBaseRefs_(tr),
00110 theConversionVertex_(convVtx),
00111 nSharedHits_(0),
00112 algorithm_(algo),
00113 qualityMask_(0)
00114 {
00115
00116
00117 theMinDistOfApproach_ = 9999.;
00118 theMVAout_ = 9999.;
00119 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00120 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00121 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00122 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00123 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00124 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00125 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00126 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00127
00128
00129 }
00130
00131
00132
00133 Conversion::Conversion() {
00134
00135 algorithm_=0;
00136 qualityMask_=0;
00137 theMinDistOfApproach_ = 9999.;
00138 nSharedHits_ = 0;
00139 theMVAout_ = 9999.;
00140 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00141 thePositionAtEcal_.push_back(math::XYZPointF(0.,0.,0.));
00142 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00143 theTrackInnerPosition_.push_back(math::XYZPointF(0.,0.,0.));
00144 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00145 theTrackPin_.push_back(math::XYZVectorF(0.,0.,0.));
00146 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00147 theTrackPout_.push_back(math::XYZVectorF(0.,0.,0.));
00148
00149 }
00150
00151
00152 Conversion::~Conversion() { }
00153
00154
00155 std::string const Conversion::algoNames[] = { "undefined","ecalSeeded","trackerOnly","mixed","pflow"};
00156
00157 Conversion::ConversionAlgorithm Conversion::algoByName(const std::string &name){
00158 ConversionAlgorithm size = algoSize;
00159 int index = std::find(algoNames, algoNames+size, name)-algoNames;
00160 if(index == size) return undefined;
00161
00162 return ConversionAlgorithm(index);
00163 }
00164
00165 Conversion * Conversion::clone() const {
00166 return new Conversion( * this );
00167 }
00168
00169
00170 std::vector<edm::RefToBase<reco::Track> > Conversion::tracks() const {
00171 if (trackToBaseRefs_.size() ==0 ) {
00172
00173 for (std::vector<reco::TrackRef>::const_iterator ref=tracks_.begin(); ref!=tracks_.end(); ref++ )
00174 {
00175 edm::RefToBase<reco::Track> tt(*ref);
00176 trackToBaseRefs_.push_back(tt);
00177
00178 }
00179 }
00180
00181 return trackToBaseRefs_;
00182 }
00183
00184
00185
00186 bool Conversion::isConverted() const {
00187
00188 if ( this->nTracks() == 2 )
00189 return true;
00190 else
00191 return false;
00192 }
00193
00194 double Conversion::pairInvariantMass() const{
00195 double invMass=-99.;
00196 const float mElec= 0.000511;
00197 if ( nTracks()==2 ) {
00198 double px= tracksPin()[0].x() + tracksPin()[1].x();
00199 double py= tracksPin()[0].y() + tracksPin()[1].y();
00200 double pz= tracksPin()[0].z() + tracksPin()[1].z();
00201 double mom1= tracksPin()[0].Mag2();
00202 double mom2= tracksPin()[1].Mag2();
00203 double e = sqrt( mom1+ mElec*mElec ) + sqrt( mom2 + mElec*mElec );
00204 invMass= ( e*e - px*px -py*py - pz*pz);
00205 if ( invMass>0) invMass = sqrt(invMass);
00206 else
00207 invMass = -1;
00208 }
00209
00210 return invMass;
00211 }
00212
00213 double Conversion::pairCotThetaSeparation() const {
00214 double dCotTheta=-99.;
00215
00216 if ( nTracks()==2 ) {
00217 double theta1=this->tracksPin()[0].Theta();
00218 double theta2=this->tracksPin()[1].Theta();
00219 dCotTheta = 1./tan(theta1) - 1./tan(theta2) ;
00220 }
00221
00222 return dCotTheta;
00223
00224 }
00225
00226
00227 math::XYZVectorF Conversion::pairMomentum() const {
00228
00229 if ( nTracks()==2 ) {
00230 return this->tracksPin()[0] + this->tracksPin()[1];
00231 }
00232 return math::XYZVectorF(0.,0.,0.);
00233
00234
00235
00236 }
00237
00238
00239 math::XYZTLorentzVectorF Conversion::refittedPair4Momentum() const {
00240
00241 math::XYZTLorentzVectorF p4;
00242 if ( this->conversionVertex().isValid() )
00243 p4 = this->conversionVertex().p4( 0.000511, 0.5);
00244
00245 return p4;
00246
00247
00248 }
00249
00250
00251
00252 math::XYZVectorF Conversion::refittedPairMomentum() const {
00253
00254 if ( this->conversionVertex().isValid() ) {
00255 return this->refittedPair4Momentum().Vect();
00256 }
00257 return math::XYZVectorF(0.,0.,0.);
00258
00259 }
00260
00261
00262
00263 double Conversion::EoverP() const {
00264
00265
00266 double ep=-99.;
00267
00268 if ( nTracks() > 0 ) {
00269 unsigned int size= this->caloCluster().size();
00270 float etot=0.;
00271 for ( unsigned int i=0; i<size; i++) {
00272 etot+= caloCluster()[i]->energy();
00273 }
00274 if (this->pairMomentum().Mag2() !=0) ep= etot/sqrt(this->pairMomentum().Mag2());
00275 }
00276
00277
00278
00279 return ep;
00280
00281 }
00282
00283
00284
00285 double Conversion::EoverPrefittedTracks() const {
00286
00287
00288 double ep=-99.;
00289
00290 if ( nTracks() > 0 ) {
00291 unsigned int size= this->caloCluster().size();
00292 float etot=0.;
00293 for ( unsigned int i=0; i<size; i++) {
00294 etot+= caloCluster()[i]->energy();
00295 }
00296 if (this->refittedPairMomentum().Mag2() !=0) ep= etot/sqrt(this->refittedPairMomentum().Mag2());
00297 }
00298
00299
00300
00301 return ep;
00302
00303 }
00304
00305
00306
00307 std::vector<double> Conversion::tracksSigned_d0() const {
00308 std::vector<double> result;
00309
00310 for (unsigned int i=0; i< nTracks(); i++)
00311 result.push_back(tracks()[i]->d0()* tracks()[i]->charge()) ;
00312
00313 return result;
00314
00315
00316 }
00317
00318 double Conversion::dPhiTracksAtVtx() const {
00319 double result=-99;
00320 if ( nTracks()==2 ) {
00321 result = tracksPin()[0].phi() - tracksPin()[1].phi();
00322 if( result > pi) { result = result - twopi;}
00323 if( result < -pi) { result = result + twopi;}
00324 }
00325
00326 return result;
00327
00328
00329 }
00330
00331 double Conversion::dPhiTracksAtEcal() const {
00332 double result =-99.;
00333
00334 if ( nTracks()==2 && bcMatchingWithTracks()[0].isNonnull() && bcMatchingWithTracks()[1].isNonnull() ) {
00335
00336 float recoPhi1 = ecalImpactPosition()[0].phi();
00337 if( recoPhi1 > pi) { recoPhi1 = recoPhi1 - twopi;}
00338 if( recoPhi1 < -pi) { recoPhi1 = recoPhi1 + twopi;}
00339
00340 float recoPhi2 = ecalImpactPosition()[1].phi();
00341 if( recoPhi2 > pi) { recoPhi2 = recoPhi2 - twopi;}
00342 if( recoPhi2 < -pi) { recoPhi2 = recoPhi2 + twopi;}
00343
00344 result = recoPhi1 -recoPhi2;
00345
00346 if( result > pi) { result = result - twopi;}
00347 if( result < -pi) { result = result + twopi;}
00348
00349 }
00350
00351 return result;
00352
00353
00354 }
00355
00356 double Conversion::dEtaTracksAtEcal() const {
00357 double result=-99.;
00358
00359
00360 if ( nTracks()==2 && bcMatchingWithTracks()[0].isNonnull() && bcMatchingWithTracks()[1].isNonnull() ) {
00361
00362 result =ecalImpactPosition()[0].eta() - ecalImpactPosition()[1].eta();
00363
00364 }
00365
00366
00367
00368 return result;
00369
00370
00371 }
00372
00373 double Conversion::dxy(const math::XYZPoint& myBeamSpot) const {
00374
00375 const reco::Vertex &vtx = conversionVertex();
00376 if (!vtx.isValid()) return -9999.;
00377
00378 math::XYZVectorF mom = refittedPairMomentum();
00379
00380 double dxy = (-(vtx.x() - myBeamSpot.x())*mom.y() + (vtx.y() - myBeamSpot.y())*mom.x())/mom.rho();
00381 return dxy;
00382
00383 }
00384
00385 double Conversion::dz(const math::XYZPoint& myBeamSpot) const {
00386
00387 const reco::Vertex &vtx = conversionVertex();
00388 if (!vtx.isValid()) return -9999.;
00389
00390 math::XYZVectorF mom = refittedPairMomentum();
00391
00392 double dz = (vtx.z()-myBeamSpot.z()) - ((vtx.x()-myBeamSpot.x())*mom.x()+(vtx.y()-myBeamSpot.y())*mom.y())/mom.rho() * mom.z()/mom.rho();
00393 return dz;
00394
00395 }
00396
00397 double Conversion::lxy(const math::XYZPoint& myBeamSpot) const {
00398
00399 const reco::Vertex &vtx = conversionVertex();
00400 if (!vtx.isValid()) return -9999.;
00401
00402 math::XYZVectorF mom = refittedPairMomentum();
00403
00404 double dbsx = vtx.x() - myBeamSpot.x();
00405 double dbsy = vtx.y() - myBeamSpot.y();
00406 double lxy = (mom.x()*dbsx + mom.y()*dbsy)/mom.rho();
00407 return lxy;
00408
00409 }
00410
00411 double Conversion::lz(const math::XYZPoint& myBeamSpot) const {
00412
00413 const reco::Vertex &vtx = conversionVertex();
00414 if (!vtx.isValid()) return -9999.;
00415
00416 math::XYZVectorF mom = refittedPairMomentum();
00417
00418 double lz = (vtx.z() - myBeamSpot.z())*mom.z()/std::abs(mom.z());
00419 return lz;
00420
00421 }
00422
1.7.3