Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:59 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ConversionProducer
0004 // Class:      ConversionProducer
0005 //
0006 /**\class ConversionProducer 
0007 
0008 Description: Produces converted photon objects using default track collections
0009 
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Authors:  Hongliang Liu, UC of Riverside US, Nancy Marinelli Univ of Notre Dame
0015 //         Created:  Thu Mar 13 17:40:48 CDT 2008
0016 //
0017 //
0018 
0019 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0020 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0021 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0022 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0023 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0024 #include "DataFormats/EgammaTrackReco/interface/ConversionTrack.h"
0025 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0026 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0027 #include "DataFormats/GeometrySurface/interface/ReferenceCounted.h"
0028 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
0029 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0030 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0031 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0032 #include "DataFormats/Math/interface/deltaPhi.h"
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0035 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0036 #include "DataFormats/VertexReco/interface/Vertex.h"
0037 #include "FWCore/Framework/interface/ESHandle.h"
0038 #include "FWCore/Framework/interface/Event.h"
0039 #include "FWCore/Framework/interface/stream/EDProducer.h"
0040 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/Utilities/interface/ESGetToken.h"
0043 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0046 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0047 #include "MagneticField/Engine/interface/MagneticField.h"
0048 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0049 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionHitChecker.h"
0050 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionVertexFinder.h"
0051 #include "RecoEgamma/EgammaPhotonAlgos/interface/TangentApproachInRPhi.h"
0052 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0053 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0054 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0055 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0056 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0057 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0058 #include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
0059 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0060 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0061 
0062 #include <map>
0063 #include <memory>
0064 
0065 class ConversionProducer : public edm::stream::EDProducer<> {
0066 public:
0067   explicit ConversionProducer(const edm::ParameterSet&);
0068 
0069 private:
0070   void produce(edm::Event&, const edm::EventSetup&) override;
0071 
0072   void buildSuperAndBasicClusterGeoMap(const edm::Event&,
0073                                        std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
0074                                        std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs);
0075 
0076   // ----------member data ---------------------------
0077   std::string algoName_;
0078 
0079   typedef math::XYZPointF Point;
0080   typedef std::vector<Point> PointCollection;
0081 
0082   edm::EDGetTokenT<edm::View<reco::ConversionTrack> > src_;
0083 
0084   edm::EDGetTokenT<edm::View<reco::CaloCluster> > scBarrelProducer_;
0085   edm::EDGetTokenT<edm::View<reco::CaloCluster> > scEndcapProducer_;
0086   edm::EDGetTokenT<edm::View<reco::CaloCluster> > bcBarrelCollection_;
0087   edm::EDGetTokenT<edm::View<reco::CaloCluster> > bcEndcapCollection_;
0088   std::string ConvertedPhotonCollection_;
0089 
0090   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackBuilder_;
0091   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometry_;
0092   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticField_;
0093 
0094   bool allowD0_, allowDeltaPhi_, allowTrackBC_, allowDeltaCot_, allowMinApproach_, allowOppCharge_, allowVertex_;
0095 
0096   bool bypassPreselGsf_, bypassPreselEcal_, bypassPreselEcalEcal_;
0097 
0098   bool usePvtx_;  //if use primary vertices
0099   edm::EDGetTokenT<reco::VertexCollection> vertexProducer_;
0100   ConversionVertexFinder vertexFinder_;
0101 
0102   const TransientTrackBuilder* thettbuilder_;
0103 
0104   double deltaEta_;
0105 
0106   double halfWayEta_, halfWayPhi_;  //halfway open angle to search in basic clusters
0107   unsigned int maxNumOfTrackInPU_;
0108   double maxTrackZ_;
0109   double maxTrackRho_;
0110   double minSCEt_;
0111   double dEtacutForSCmatching_;
0112   double dPhicutForSCmatching_;
0113   double energyBC_;             //1.5GeV for track BC selection
0114   double energyTotalBC_;        //5GeV for track pair BC selection
0115   double d0Cut_;                //0 for d0*charge cut
0116   double dzCut_;                //innerposition of z diff cut
0117   double dEtaTkBC_, dPhiTkBC_;  //0.06 0.6 for track and BC matching
0118 
0119   double maxChi2Left_, maxChi2Right_;  //5. 5. for track chi2 quality
0120   double minHitsLeft_, minHitsRight_;  //5 2 for track hits quality
0121 
0122   double deltaCotTheta_, deltaPhi_, minApproachLow_,
0123       minApproachHigh_;  //0.02 0.2 for track pair open angle and > -0.1 cm
0124 
0125   double r_cut;     //cross_r cut
0126   double vtxChi2_;  //vertex chi2 probablity cut
0127 
0128   bool allowSingleLeg_;  //if single track conversion ?
0129   bool rightBC_;         //if right leg requires matching BC?
0130 
0131   void buildCollection(edm::Event& iEvent,
0132                        const edm::EventSetup& iSetup,
0133                        const std::multimap<float, edm::Ptr<reco::ConversionTrack> >& allTracks,
0134                        const std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs,
0135                        const std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
0136                        const reco::Vertex& the_pvtx,
0137                        reco::ConversionCollection& outputConvPhotonCollection);
0138 
0139   //track quality cut, returns pass or no
0140   inline bool trackQualityFilter(const edm::RefToBase<reco::Track>& ref, bool isLeft);
0141   inline bool trackD0Cut(const edm::RefToBase<reco::Track>& ref);
0142   inline bool trackD0Cut(const edm::RefToBase<reco::Track>& ref, const reco::Vertex& the_pvtx);
0143 
0144   //track impact point at ECAL wall, returns validity to access position ew
0145   bool getTrackImpactPosition(const reco::Track* tk_ref,
0146                               TrackerGeometry const& trackerGeom,
0147                               MagneticField const& magField,
0148                               math::XYZPointF& ew);
0149 
0150   //distance at min approaching point, returns distance
0151   //      double getMinApproach(const edm::RefToBase<reco::Track>& ll, const edm::RefToBase<reco::Track>& rr,
0152   //          const MagneticField* magField);
0153 
0154   bool preselectTrackPair(const reco::TransientTrack& ttk_l, const reco::TransientTrack& ttk_r, double& appDist);
0155 
0156   //cut-based selection, TODO remove global cut variables
0157   bool checkTrackPair(const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& ll,
0158                       const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& rr);
0159 
0160   //kinematic vertex fitting, return true for valid vertex
0161   inline bool checkVertex(const reco::TransientTrack& ttk_l,
0162                           const reco::TransientTrack& ttk_r,
0163                           MagneticField const& magField,
0164                           reco::Vertex& the_vertex) {
0165     return vertexFinder_.run({ttk_l, ttk_r}, the_vertex);
0166   }
0167 
0168   bool checkPhi(const edm::RefToBase<reco::Track>& tk_l,
0169                 const edm::RefToBase<reco::Track>& tk_r,
0170                 TrackerGeometry const& trackerGeom,
0171                 MagneticField const& magField,
0172                 const reco::Vertex& the_vertex);
0173 
0174   //check the closest BC, returns true for found a BC
0175   bool getMatchedBC(const std::multimap<double, reco::CaloClusterPtr>& bcMap,
0176                     const math::XYZPointF& trackImpactPosition,
0177                     reco::CaloClusterPtr& closestBC);
0178 
0179   // finds the super cluster matching with at least one track in the pair
0180   bool matchingSC(const std::multimap<double, reco::CaloClusterPtr>& scMap,
0181                   reco::Conversion& conv,
0182                   reco::CaloClusterPtrVector& mSC);
0183 
0184   double etaTransformation(float EtaParticle, float Zvertex);
0185 
0186   math::XYZPointF toFConverterP(const math::XYZPoint& val) { return math::XYZPointF(val.x(), val.y(), val.z()); }
0187 
0188   math::XYZVectorF toFConverterV(const math::XYZVector& val) { return math::XYZVectorF(val.x(), val.y(), val.z()); }
0189 };
0190 
0191 #include "FWCore/Framework/interface/MakerMacros.h"
0192 DEFINE_FWK_MODULE(ConversionProducer);
0193 
0194 inline const GeomDet* recHitDet(const TrackingRecHit& hit, const TrackingGeometry* geom) {
0195   return geom->idToDet(hit.geographicalId());
0196 }
0197 
0198 inline const BoundPlane& recHitSurface(const TrackingRecHit& hit, const TrackingGeometry* geom) {
0199   return recHitDet(hit, geom)->surface();
0200 }
0201 
0202 inline LocalVector toLocal(const reco::Track::Vector& v, const Surface& s) {
0203   return s.toLocal(GlobalVector(v.x(), v.y(), v.z()));
0204 }
0205 
0206 ConversionProducer::ConversionProducer(const edm::ParameterSet& iConfig) : vertexFinder_{iConfig} {
0207   algoName_ = iConfig.getParameter<std::string>("AlgorithmName");
0208 
0209   src_ = consumes<edm::View<reco::ConversionTrack> >(iConfig.getParameter<edm::InputTag>("src"));
0210 
0211   maxNumOfTrackInPU_ = iConfig.getParameter<int>("maxNumOfTrackInPU");
0212   maxTrackRho_ = iConfig.getParameter<double>("maxTrackRho");
0213   maxTrackZ_ = iConfig.getParameter<double>("maxTrackZ");
0214 
0215   allowTrackBC_ = iConfig.getParameter<bool>("AllowTrackBC");
0216   allowD0_ = iConfig.getParameter<bool>("AllowD0");
0217   allowDeltaPhi_ = iConfig.getParameter<bool>("AllowDeltaPhi");
0218   allowDeltaCot_ = iConfig.getParameter<bool>("AllowDeltaCot");
0219   allowMinApproach_ = iConfig.getParameter<bool>("AllowMinApproach");
0220   allowOppCharge_ = iConfig.getParameter<bool>("AllowOppCharge");
0221 
0222   allowVertex_ = iConfig.getParameter<bool>("AllowVertex");
0223 
0224   bypassPreselGsf_ = iConfig.getParameter<bool>("bypassPreselGsf");
0225   bypassPreselEcal_ = iConfig.getParameter<bool>("bypassPreselEcal");
0226   bypassPreselEcalEcal_ = iConfig.getParameter<bool>("bypassPreselEcalEcal");
0227 
0228   deltaEta_ = iConfig.getParameter<double>("deltaEta");
0229 
0230   halfWayEta_ = iConfig.getParameter<double>("HalfwayEta");  //open angle to search track matches with BC
0231 
0232   d0Cut_ = iConfig.getParameter<double>("d0");
0233 
0234   usePvtx_ = iConfig.getParameter<bool>("UsePvtx");  //if use primary vertices
0235 
0236   vertexProducer_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertexProducer"));
0237 
0238   transientTrackBuilder_ = esConsumes(edm::ESInputTag("", "TransientTrackBuilder"));
0239   trackerGeometry_ = esConsumes();
0240   magneticField_ = esConsumes();
0241 
0242   //Track-cluster matching eta and phi cuts
0243   dEtaTkBC_ = iConfig.getParameter<double>("dEtaTrackBC");  //TODO research on cut endcap/barrel
0244   dPhiTkBC_ = iConfig.getParameter<double>("dPhiTrackBC");
0245 
0246   bcBarrelCollection_ =
0247       consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("bcBarrelCollection"));
0248   bcEndcapCollection_ =
0249       consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("bcEndcapCollection"));
0250 
0251   scBarrelProducer_ = consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("scBarrelProducer"));
0252   scEndcapProducer_ = consumes<edm::View<reco::CaloCluster> >(iConfig.getParameter<edm::InputTag>("scEndcapProducer"));
0253 
0254   energyBC_ = iConfig.getParameter<double>("EnergyBC");            //BC energy threshold
0255   energyTotalBC_ = iConfig.getParameter<double>("EnergyTotalBC");  //BC energy threshold
0256   minSCEt_ = iConfig.getParameter<double>("minSCEt");              //super cluster energy threshold
0257   dEtacutForSCmatching_ = iConfig.getParameter<double>(
0258       "dEtacutForSCmatching");  // dEta between conversion momentum direction and SC position
0259   dPhicutForSCmatching_ = iConfig.getParameter<double>(
0260       "dPhicutForSCmatching");  // dPhi between conversion momentum direction and SC position
0261 
0262   //Track cuts on left right track: at least one leg reaches ECAL
0263   //Left track: must exist, must reach Ecal and match BC, so loose cut on Chi2 and tight on hits
0264   //Right track: not necessary to exist (if allowSingleLeg_), not necessary to reach ECAL or match BC, so tight cut on Chi2 and loose on hits
0265   maxChi2Left_ = iConfig.getParameter<double>("MaxChi2Left");
0266   maxChi2Right_ = iConfig.getParameter<double>("MaxChi2Right");
0267   minHitsLeft_ = iConfig.getParameter<int>("MinHitsLeft");
0268   minHitsRight_ = iConfig.getParameter<int>("MinHitsRight");
0269 
0270   //Track Open angle cut on delta cot(theta) and delta phi
0271   deltaCotTheta_ = iConfig.getParameter<double>("DeltaCotTheta");
0272   deltaPhi_ = iConfig.getParameter<double>("DeltaPhi");
0273   minApproachLow_ = iConfig.getParameter<double>("MinApproachLow");
0274   minApproachHigh_ = iConfig.getParameter<double>("MinApproachHigh");
0275 
0276   // if allow single track collection, by default False
0277   allowSingleLeg_ = iConfig.getParameter<bool>("AllowSingleLeg");
0278   rightBC_ = iConfig.getParameter<bool>("AllowRightBC");
0279 
0280   //track inner position dz cut, need RECO
0281   dzCut_ = iConfig.getParameter<double>("dz");
0282   //track analytical cross cut
0283   r_cut = iConfig.getParameter<double>("rCut");
0284   vtxChi2_ = iConfig.getParameter<double>("vtxChi2");
0285 
0286   thettbuilder_ = nullptr;
0287 
0288   //output
0289   ConvertedPhotonCollection_ = iConfig.getParameter<std::string>("convertedPhotonCollection");
0290 
0291   produces<reco::ConversionCollection>(ConvertedPhotonCollection_);
0292 }
0293 
0294 // ------------ method called to produce the data  ------------
0295 void ConversionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0296   using namespace edm;
0297 
0298   reco::ConversionCollection outputConvPhotonCollection;
0299   auto outputConvPhotonCollection_p = std::make_unique<reco::ConversionCollection>();
0300 
0301   //std::cout << " ConversionProducer::produce " << std::endl;
0302   //Read multiple track input collections
0303 
0304   edm::Handle<edm::View<reco::ConversionTrack> > trackCollectionHandle;
0305   iEvent.getByToken(src_, trackCollectionHandle);
0306 
0307   //build map of ConversionTracks ordered in eta
0308   std::multimap<float, edm::Ptr<reco::ConversionTrack> > convTrackMap;
0309   for (auto const& t : trackCollectionHandle->ptrs())
0310     convTrackMap.emplace(t->track()->eta(), t);
0311 
0312   edm::Handle<reco::VertexCollection> vertexHandle;
0313   reco::VertexCollection vertexCollection;
0314   if (usePvtx_) {
0315     iEvent.getByToken(vertexProducer_, vertexHandle);
0316     if (!vertexHandle.isValid()) {
0317       edm::LogError("ConversionProducer") << "Error! Can't get the product primary Vertex Collection "
0318                                           << "\n";
0319       usePvtx_ = false;
0320     }
0321     if (usePvtx_)
0322       vertexCollection = *(vertexHandle.product());
0323   }
0324 
0325   thettbuilder_ = &iSetup.getData(transientTrackBuilder_);
0326 
0327   reco::Vertex the_pvtx;
0328   //because the priamry vertex is sorted by quality, the first one is the best
0329   if (!vertexCollection.empty())
0330     the_pvtx = *(vertexCollection.begin());
0331 
0332   if (trackCollectionHandle->size() > maxNumOfTrackInPU_) {
0333     iEvent.put(std::move(outputConvPhotonCollection_p), ConvertedPhotonCollection_);
0334     return;
0335   }
0336 
0337   // build Super and Basic cluster geometry map to search in eta bounds for clusters
0338   std::multimap<double, reco::CaloClusterPtr> basicClusterPtrs;
0339   std::multimap<double, reco::CaloClusterPtr> superClusterPtrs;
0340 
0341   buildSuperAndBasicClusterGeoMap(iEvent, basicClusterPtrs, superClusterPtrs);
0342 
0343   buildCollection(iEvent,
0344                   iSetup,
0345                   convTrackMap,
0346                   superClusterPtrs,
0347                   basicClusterPtrs,
0348                   the_pvtx,
0349                   outputConvPhotonCollection);  //allow empty basicClusterPtrs
0350 
0351   outputConvPhotonCollection_p->assign(outputConvPhotonCollection.begin(), outputConvPhotonCollection.end());
0352   iEvent.put(std::move(outputConvPhotonCollection_p), ConvertedPhotonCollection_);
0353 }
0354 
0355 void ConversionProducer::buildSuperAndBasicClusterGeoMap(const edm::Event& iEvent,
0356                                                          std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
0357                                                          std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs) {
0358   // Get the Super Cluster collection in the Barrel
0359   edm::Handle<edm::View<reco::CaloCluster> > scBarrelHandle;
0360   iEvent.getByToken(scBarrelProducer_, scBarrelHandle);
0361   if (!scBarrelHandle.isValid()) {
0362     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the barrel superclusters!";
0363   }
0364 
0365   // Get the Super Cluster collection in the Endcap
0366   edm::Handle<edm::View<reco::CaloCluster> > scEndcapHandle;
0367   iEvent.getByToken(scEndcapProducer_, scEndcapHandle);
0368   if (!scEndcapHandle.isValid()) {
0369     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the endcap superclusters!";
0370   }
0371 
0372   edm::Handle<edm::View<reco::CaloCluster> > bcBarrelHandle;
0373   edm::Handle<edm::View<reco::CaloCluster> > bcEndcapHandle;  //TODO check cluster type if BasicCluster or PFCluster
0374 
0375   iEvent.getByToken(bcBarrelCollection_, bcBarrelHandle);
0376   if (!bcBarrelHandle.isValid()) {
0377     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the barrel basic clusters!";
0378   }
0379 
0380   iEvent.getByToken(bcEndcapCollection_, bcEndcapHandle);
0381   if (!bcEndcapHandle.isValid()) {
0382     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the endcap basic clusters!";
0383   }
0384 
0385   if (bcBarrelHandle.isValid()) {
0386     for (auto const& handle : {bcBarrelHandle, bcEndcapHandle}) {
0387       for (auto const& bc : handle->ptrs()) {
0388         if (bc->energy() > energyBC_)
0389           basicClusterPtrs.emplace(bc->position().eta(), bc);
0390       }
0391     }
0392   }
0393 
0394   if (scBarrelHandle.isValid()) {
0395     for (auto const& handle : {scBarrelHandle, scEndcapHandle}) {
0396       for (auto const& sc : handle->ptrs()) {
0397         if (sc->energy() > minSCEt_)
0398           superClusterPtrs.emplace(sc->position().eta(), sc);
0399       }
0400     }
0401   }
0402 }
0403 
0404 void ConversionProducer::buildCollection(edm::Event& iEvent,
0405                                          const edm::EventSetup& iSetup,
0406                                          const std::multimap<float, edm::Ptr<reco::ConversionTrack> >& allTracks,
0407                                          const std::multimap<double, reco::CaloClusterPtr>& superClusterPtrs,
0408                                          const std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
0409                                          const reco::Vertex& the_pvtx,
0410                                          reco::ConversionCollection& outputConvPhotonCollection) {
0411   TrackerGeometry const& trackerGeom = iSetup.getData(trackerGeometry_);
0412   MagneticField const& magField = iSetup.getData(magneticField_);
0413 
0414   //   std::vector<math::XYZPointF> trackImpactPosition;
0415   //   trackImpactPosition.reserve(allTracks.size());//track impact position at ECAL
0416   //   std::vector<bool> trackValidECAL;//Does this track reach ECAL basic cluster (reach ECAL && match with BC)
0417   //   trackValidECAL.assign(allTracks.size(), false);
0418   //
0419   //   std::vector<reco::CaloClusterPtr> trackMatchedBC;
0420   //   reco::CaloClusterPtr empty_bc;
0421   //   trackMatchedBC.assign(allTracks.size(), empty_bc);//TODO find a better way to avoid copy constructor
0422   //
0423   //   std::vector<int> bcHandleId;//the associated BC handle id, -1 invalid, 0 barrel 1 endcap
0424   //   bcHandleId.assign(allTracks.size(), -1);
0425 
0426   // not used    std::multimap<double, int> trackInnerEta;//Track innermost state Eta map to TrackRef index, to be used in track pair sorting
0427 
0428   std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF> trackImpactPosition;
0429   std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr> trackMatchedBC;
0430 
0431   ConversionHitChecker hitChecker;
0432 
0433   //2 propagate all tracks into ECAL, record its eta and phi
0434 
0435   for (auto const& tk_ref : allTracks) {
0436     const reco::Track* tk = tk_ref.second->trackRef().get();
0437 
0438     //check impact position then match with BC
0439     math::XYZPointF ew;
0440     if (getTrackImpactPosition(tk, trackerGeom, magField, ew)) {
0441       trackImpactPosition[tk_ref.second] = ew;
0442 
0443       reco::CaloClusterPtr closest_bc;  //the closest matching BC to track
0444 
0445       if (getMatchedBC(basicClusterPtrs, ew, closest_bc)) {
0446         trackMatchedBC[tk_ref.second] = closest_bc;
0447       }
0448     }
0449   }
0450 
0451   //3. pair up tracks:
0452   //TODO it is k-Closest pair of point problem
0453   //std::cout << " allTracks.size() " <<  allTracks.size() << std::endl;
0454   for (auto ll = allTracks.begin(); ll != allTracks.end(); ++ll) {
0455     bool track1HighPurity = true;
0456     //std::cout << " Loop on allTracks " << std::endl;
0457     const edm::RefToBase<reco::Track>& left = ll->second->trackRef();
0458 
0459     //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
0460     //(Note that the TrackRef and GsfTrackRef versions of the constructor are needed
0461     // to properly get refit tracks in the output vertex)
0462     reco::TransientTrack ttk_l;
0463     if (dynamic_cast<const reco::GsfTrack*>(left.get())) {
0464       ttk_l = thettbuilder_->build(left.castTo<reco::GsfTrackRef>());
0465     } else {
0466       ttk_l = thettbuilder_->build(left.castTo<reco::TrackRef>());
0467     }
0468 
0469     ////  Remove BC matching from track selection
0470     //      if ((allowTrackBC_ && !trackValidECAL[ll-allTracks.begin()]) )//this Left leg should have valid BC
0471     // continue;
0472 
0473     if (the_pvtx.isValid()) {
0474       if (!(trackD0Cut(left, the_pvtx)))
0475         track1HighPurity = false;
0476     } else {
0477       if (!(trackD0Cut(left)))
0478         track1HighPurity = false;
0479     }
0480 
0481     std::vector<int> right_candidates;  //store all right legs passed the cut (theta/approach and ref pair)
0482     std::vector<double> right_candidate_theta, right_candidate_approach;
0483     std::vector<std::pair<bool, reco::Vertex> > vertex_candidates;
0484 
0485     //inner loop only over tracks between eta and eta + deltaEta of the first track
0486     float etasearch = ll->first + deltaEta_;
0487     std::multimap<float, edm::Ptr<reco::ConversionTrack> >::const_iterator rr = ll;
0488     ++rr;
0489     for (; rr != allTracks.lower_bound(etasearch); ++rr) {
0490       bool track2HighPurity = true;
0491       bool highPurityPair = true;
0492 
0493       const edm::RefToBase<reco::Track>& right = rr->second->trackRef();
0494 
0495       //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
0496       reco::TransientTrack ttk_r;
0497       if (dynamic_cast<const reco::GsfTrack*>(right.get())) {
0498         ttk_r = thettbuilder_->build(right.castTo<reco::GsfTrackRef>());
0499       } else {
0500         ttk_r = thettbuilder_->build(right.castTo<reco::TrackRef>());
0501       }
0502       //std::cout << " This track is " <<  right->algoName() << std::endl;
0503 
0504       //all vertexing preselection should go here
0505 
0506       //check for opposite charge
0507       if (allowOppCharge_ && (left->charge() * right->charge() > 0))
0508         continue;  //same sign, reject pair
0509 
0510       ////  Remove BC matching from track selection
0511       //if ( (allowTrackBC_ && !trackValidECAL[rr-allTracks.begin()] && rightBC_) )// if right track matches ECAL
0512       //  continue;
0513 
0514       double approachDist = -999.;
0515       //apply preselection to track pair, overriding preselection for gsf+X or ecalseeded+X pairs if so configured
0516       bool preselected = preselectTrackPair(ttk_l, ttk_r, approachDist);
0517       preselected = preselected || (bypassPreselGsf_ &&
0518                                     (left->algo() == reco::TrackBase::gsf || right->algo() == reco::TrackBase::gsf));
0519       preselected = preselected || (bypassPreselEcal_ && (left->algo() == reco::TrackBase::outInEcalSeededConv ||
0520                                                           right->algo() == reco::TrackBase::outInEcalSeededConv ||
0521                                                           left->algo() == reco::TrackBase::inOutEcalSeededConv ||
0522                                                           right->algo() == reco::TrackBase::inOutEcalSeededConv));
0523       preselected = preselected || (bypassPreselEcalEcal_ &&
0524                                     (left->algo() == reco::TrackBase::outInEcalSeededConv ||
0525                                      left->algo() == reco::TrackBase::inOutEcalSeededConv) &&
0526                                     (right->algo() == reco::TrackBase::outInEcalSeededConv ||
0527                                      right->algo() == reco::TrackBase::inOutEcalSeededConv));
0528 
0529       if (!preselected) {
0530         continue;
0531       }
0532 
0533       //do the actual vertex fit
0534       reco::Vertex theConversionVertex;  //by default it is invalid
0535       bool goodVertex = checkVertex(ttk_l, ttk_r, magField, theConversionVertex);
0536 
0537       //bail as early as possible in case the fit didn't return a good vertex
0538       if (!goodVertex) {
0539         continue;
0540       }
0541 
0542       //track pair pass the quality cut
0543       if (!((trackQualityFilter(left, true) && trackQualityFilter(right, false)) ||
0544             (trackQualityFilter(left, false) && trackQualityFilter(right, true)))) {
0545         highPurityPair = false;
0546       }
0547 
0548       if (the_pvtx.isValid()) {
0549         if (!(trackD0Cut(right, the_pvtx)))
0550           track2HighPurity = false;
0551       } else {
0552         if (!(trackD0Cut(right)))
0553           track2HighPurity = false;
0554       }
0555 
0556       //if all cuts passed, go ahead to make conversion candidates
0557       std::vector<edm::RefToBase<reco::Track> > trackPairRef;
0558       trackPairRef.push_back(left);   //left track
0559       trackPairRef.push_back(right);  //right track
0560 
0561       std::vector<math::XYZVectorF> trackPin;
0562       std::vector<math::XYZVectorF> trackPout;
0563       std::vector<math::XYZPointF> trackInnPos;
0564       std::vector<uint8_t> nHitsBeforeVtx;
0565       std::vector<Measurement1DFloat> dlClosestHitToVtx;
0566 
0567       if (left->extra().isNonnull() && right->extra().isNonnull()) {  //only available on TrackExtra
0568         trackInnPos.push_back(toFConverterP(left->innerPosition()));
0569         trackInnPos.push_back(toFConverterP(right->innerPosition()));
0570         trackPin.push_back(toFConverterV(left->innerMomentum()));
0571         trackPin.push_back(toFConverterV(right->innerMomentum()));
0572         trackPout.push_back(toFConverterV(left->outerMomentum()));
0573         trackPout.push_back(toFConverterV(right->outerMomentum()));
0574         auto leftWrongHits = hitChecker.nHitsBeforeVtx(*left->extra(), theConversionVertex);
0575         auto rightWrongHits = hitChecker.nHitsBeforeVtx(*right->extra(), theConversionVertex);
0576         nHitsBeforeVtx.push_back(leftWrongHits.first);
0577         nHitsBeforeVtx.push_back(rightWrongHits.first);
0578         dlClosestHitToVtx.push_back(leftWrongHits.second);
0579         dlClosestHitToVtx.push_back(rightWrongHits.second);
0580       }
0581 
0582       uint8_t nSharedHits = hitChecker.nSharedHits(*left.get(), *right.get());
0583 
0584       //if using kinematic fit, check with chi2 post cut
0585       if (theConversionVertex.isValid()) {
0586         const float chi2Prob = ChiSquaredProbability(theConversionVertex.chi2(), theConversionVertex.ndof());
0587         if (chi2Prob < vtxChi2_)
0588           highPurityPair = false;
0589       }
0590 
0591       //std::cout << "  highPurityPair after vertex cut " <<  highPurityPair << std::endl;
0592       std::vector<math::XYZPointF> trkPositionAtEcal;
0593       std::vector<reco::CaloClusterPtr> matchingBC;
0594 
0595       if (allowTrackBC_) {  //TODO find out the BC ptrs if not doing matching, otherwise, leave it empty
0596         //const int lbc_handle = bcHandleId[ll-allTracks.begin()],
0597         //        rbc_handle = bcHandleId[rr-allTracks.begin()];
0598 
0599         std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF>::const_iterator trackImpactPositionLeft =
0600             trackImpactPosition.find(ll->second);
0601         std::map<edm::Ptr<reco::ConversionTrack>, math::XYZPointF>::const_iterator trackImpactPositionRight =
0602             trackImpactPosition.find(rr->second);
0603         std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr>::const_iterator trackMatchedBCLeft =
0604             trackMatchedBC.find(ll->second);
0605         std::map<edm::Ptr<reco::ConversionTrack>, reco::CaloClusterPtr>::const_iterator trackMatchedBCRight =
0606             trackMatchedBC.find(rr->second);
0607 
0608         if (trackImpactPositionLeft != trackImpactPosition.end()) {
0609           trkPositionAtEcal.push_back(trackImpactPositionLeft->second);  //left track
0610         } else {
0611           trkPositionAtEcal.push_back(math::XYZPointF());  //left track
0612         }
0613         if (trackImpactPositionRight != trackImpactPosition.end()) {  //second track ECAL position may be invalid
0614           trkPositionAtEcal.push_back(trackImpactPositionRight->second);
0615         }
0616 
0617         double total_e_bc = 0.;
0618         if (trackMatchedBCLeft != trackMatchedBC.end()) {
0619           matchingBC.push_back(trackMatchedBCLeft->second);  //left track
0620           total_e_bc += trackMatchedBCLeft->second->energy();
0621         } else {
0622           matchingBC.push_back(reco::CaloClusterPtr());  //left track
0623         }
0624         if (trackMatchedBCRight != trackMatchedBC.end()) {  //second track ECAL position may be invalid
0625           matchingBC.push_back(trackMatchedBCRight->second);
0626           total_e_bc += trackMatchedBCRight->second->energy();
0627         }
0628 
0629         if (total_e_bc < energyTotalBC_) {
0630           highPurityPair = false;
0631         }
0632       }
0633       //signature cuts, then check if vertex, then post-selection cuts
0634       highPurityPair = highPurityPair && track1HighPurity && track2HighPurity && goodVertex &&
0635                        checkPhi(left, right, trackerGeom, magField, theConversionVertex);
0636 
0637       /// match the track pair with a SC. If at least one track matches, store the SC
0638       /*
0639         for ( std::vector<edm::RefToBase<reco::Track> >::iterator iTk=trackPairRef.begin();  iTk!=trackPairRef.end(); iTk++) {
0640         math::XYZPointF impPos;
0641         if ( getTrackImpactPosition(*iTk, trackerGeom, magField, impPos) ) {
0642            
0643         }
0644 
0645         }
0646       */
0647 
0648       const float minAppDist = approachDist;
0649       reco::Conversion::ConversionAlgorithm algo = reco::Conversion::algoByName(algoName_);
0650       float dummy = 0;
0651       reco::CaloClusterPtrVector scPtrVec;
0652       reco::Conversion newCandidate(scPtrVec,
0653                                     trackPairRef,
0654                                     trkPositionAtEcal,
0655                                     theConversionVertex,
0656                                     matchingBC,
0657                                     minAppDist,
0658                                     trackInnPos,
0659                                     trackPin,
0660                                     trackPout,
0661                                     nHitsBeforeVtx,
0662                                     dlClosestHitToVtx,
0663                                     nSharedHits,
0664                                     dummy,
0665                                     algo);
0666       // Fill in scPtrVec with the macthing SC
0667       if (matchingSC(superClusterPtrs, newCandidate, scPtrVec))
0668         newCandidate.setMatchingSuperCluster(scPtrVec);
0669 
0670       newCandidate.setQuality(reco::Conversion::highPurity, highPurityPair);
0671       bool generalTracksOnly = ll->second->isTrackerOnly() && rr->second->isTrackerOnly() &&
0672                                !dynamic_cast<const reco::GsfTrack*>(ll->second->trackRef().get()) &&
0673                                !dynamic_cast<const reco::GsfTrack*>(rr->second->trackRef().get());
0674       bool gsfTracksOpenOnly = ll->second->isGsfTrackOpen() && rr->second->isGsfTrackOpen();
0675       bool arbitratedEcalSeeded = ll->second->isArbitratedEcalSeeded() && rr->second->isArbitratedEcalSeeded();
0676       bool arbitratedMerged = ll->second->isArbitratedMerged() && rr->second->isArbitratedMerged();
0677       bool arbitratedMergedEcalGeneral =
0678           ll->second->isArbitratedMergedEcalGeneral() && rr->second->isArbitratedMergedEcalGeneral();
0679 
0680       newCandidate.setQuality(reco::Conversion::generalTracksOnly, generalTracksOnly);
0681       newCandidate.setQuality(reco::Conversion::gsfTracksOpenOnly, gsfTracksOpenOnly);
0682       newCandidate.setQuality(reco::Conversion::arbitratedEcalSeeded, arbitratedEcalSeeded);
0683       newCandidate.setQuality(reco::Conversion::arbitratedMerged, arbitratedMerged);
0684       newCandidate.setQuality(reco::Conversion::arbitratedMergedEcalGeneral, arbitratedMergedEcalGeneral);
0685 
0686       outputConvPhotonCollection.push_back(newCandidate);
0687     }
0688   }
0689 }
0690 
0691 //
0692 // member functions
0693 //
0694 
0695 inline bool ConversionProducer::trackQualityFilter(const edm::RefToBase<reco::Track>& ref, bool isLeft) {
0696   bool pass = true;
0697   if (isLeft) {
0698     pass = (ref->normalizedChi2() < maxChi2Left_ && ref->found() >= minHitsLeft_);
0699   } else {
0700     pass = (ref->normalizedChi2() < maxChi2Right_ && ref->found() >= minHitsRight_);
0701   }
0702 
0703   return pass;
0704 }
0705 
0706 inline bool ConversionProducer::trackD0Cut(const edm::RefToBase<reco::Track>& ref) {
0707   //NOTE if not allow d0 cut, always true
0708   return ((!allowD0_) || !(ref->d0() * ref->charge() / ref->d0Error() < d0Cut_));
0709 }
0710 
0711 inline bool ConversionProducer::trackD0Cut(const edm::RefToBase<reco::Track>& ref, const reco::Vertex& the_pvtx) {
0712   //
0713   return ((!allowD0_) || !(-ref->dxy(the_pvtx.position()) * ref->charge() / ref->dxyError() < d0Cut_));
0714 }
0715 
0716 bool ConversionProducer::getTrackImpactPosition(const reco::Track* tk_ref,
0717                                                 TrackerGeometry const& trackerGeom,
0718                                                 MagneticField const& magField,
0719                                                 math::XYZPointF& ew) {
0720   PropagatorWithMaterial propag(alongMomentum, 0.000511, &magField);
0721 
0722   ReferenceCountingPointer<Surface> ecalWall(new BoundCylinder(
0723       129.f, GlobalPoint(0., 0., 0.), TkRotation<float>(), new SimpleCylinderBounds(129, 129, -320.5, 320.5)));
0724   constexpr float epsilon = 0.001;
0725   Surface::RotationType rot;  // unit rotation matrix
0726   constexpr float barrelRadius = 129.f;
0727   constexpr float barrelHalfLength = 270.9f;
0728   constexpr float endcapRadius = 171.1f;
0729   constexpr float endcapZ = 320.5f;
0730   ReferenceCountingPointer<BoundCylinder> theBarrel_(new BoundCylinder(
0731       barrelRadius,
0732       Surface::PositionType(0, 0, 0),
0733       rot,
0734       new SimpleCylinderBounds(barrelRadius - epsilon, barrelRadius + epsilon, -barrelHalfLength, barrelHalfLength)));
0735   ReferenceCountingPointer<BoundDisk> theNegativeEtaEndcap_(new BoundDisk(
0736       Surface::PositionType(0, 0, -endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon)));
0737   ReferenceCountingPointer<BoundDisk> thePositiveEtaEndcap_(new BoundDisk(
0738       Surface::PositionType(0, 0, endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon)));
0739 
0740   //const TrajectoryStateOnSurface myTSOS = trajectoryStateTransform::innerStateOnSurface(*(*ref), *trackerGeom, magField);
0741   const auto myTSOS = trajectoryStateTransform::outerStateOnSurface(*tk_ref, trackerGeom, &magField);
0742   TrajectoryStateOnSurface stateAtECAL;
0743   stateAtECAL = propag.propagate(myTSOS, *theBarrel_);
0744   if (!stateAtECAL.isValid() || (stateAtECAL.isValid() && fabs(stateAtECAL.globalPosition().eta()) > 1.479f)) {
0745     //endcap propagator
0746     if (myTSOS.globalPosition().z() > 0.) {
0747       stateAtECAL = propag.propagate(myTSOS, *thePositiveEtaEndcap_);
0748     } else {
0749       stateAtECAL = propag.propagate(myTSOS, *theNegativeEtaEndcap_);
0750     }
0751   }
0752   if (stateAtECAL.isValid()) {
0753     ew = stateAtECAL.globalPosition();
0754     return true;
0755   } else
0756     return false;
0757 }
0758 
0759 bool ConversionProducer::matchingSC(const std::multimap<double, reco::CaloClusterPtr>& scMap,
0760                                     reco::Conversion& aConv,
0761                                     // reco::CaloClusterPtr& mSC){
0762                                     reco::CaloClusterPtrVector& mSC) {
0763   //  double dRMin=999.;
0764   double detaMin = 999.;
0765   double dphiMin = 999.;
0766   reco::CaloClusterPtr match;
0767   for (std::multimap<double, reco::CaloClusterPtr>::const_iterator scItr = scMap.begin(); scItr != scMap.end();
0768        scItr++) {
0769     const reco::CaloClusterPtr& sc = scItr->second;
0770     const double delta_phi = reco::deltaPhi(aConv.refittedPairMomentum().phi(), sc->phi());
0771     double sceta = sc->eta();
0772     double conveta = etaTransformation(aConv.refittedPairMomentum().eta(), aConv.zOfPrimaryVertexFromTracks());
0773     const double delta_eta = fabs(conveta - sceta);
0774     if (fabs(delta_eta) < fabs(detaMin) && fabs(delta_phi) < fabs(dphiMin)) {
0775       detaMin = fabs(delta_eta);
0776       dphiMin = fabs(delta_phi);
0777       match = sc;
0778     }
0779   }
0780 
0781   if (fabs(detaMin) < dEtacutForSCmatching_ && fabs(dphiMin) < dPhicutForSCmatching_) {
0782     mSC.push_back(match);
0783     return true;
0784   } else
0785     return false;
0786 }
0787 
0788 bool ConversionProducer::getMatchedBC(const std::multimap<double, reco::CaloClusterPtr>& bcMap,
0789                                       const math::XYZPointF& trackImpactPosition,
0790                                       reco::CaloClusterPtr& closestBC) {
0791   const double track_eta = trackImpactPosition.eta();
0792   const double track_phi = trackImpactPosition.phi();
0793 
0794   double min_eta = 999., min_phi = 999.;
0795   reco::CaloClusterPtr closest_bc;
0796   for (std::multimap<double, reco::CaloClusterPtr>::const_iterator bc = bcMap.lower_bound(track_eta - halfWayEta_);
0797        bc != bcMap.upper_bound(track_eta + halfWayEta_);
0798        ++bc) {  //use eta map to select possible BC collection then loop in
0799     const reco::CaloClusterPtr& ebc = bc->second;
0800     const double delta_eta = track_eta - (ebc->position().eta());
0801     const double delta_phi = reco::deltaPhi(track_phi, (ebc->position().phi()));
0802     if (fabs(delta_eta) < dEtaTkBC_ && fabs(delta_phi) < dPhiTkBC_) {
0803       if (fabs(min_eta) > fabs(delta_eta) && fabs(min_phi) > fabs(delta_phi)) {  //take the closest to track BC
0804         min_eta = delta_eta;
0805         min_phi = delta_phi;
0806         closest_bc = bc->second;
0807         //TODO check if min_eta>delta_eta but min_phi<delta_phi
0808       }
0809     }
0810   }
0811 
0812   if (min_eta < 999.) {
0813     closestBC = closest_bc;
0814     return true;
0815   } else
0816     return false;
0817 }
0818 
0819 //check track open angle of phi at vertex
0820 bool ConversionProducer::checkPhi(const edm::RefToBase<reco::Track>& tk_l,
0821                                   const edm::RefToBase<reco::Track>& tk_r,
0822                                   TrackerGeometry const& trackerGeom,
0823                                   MagneticField const& magField,
0824                                   const reco::Vertex& vtx) {
0825   if (!allowDeltaPhi_)
0826     return true;
0827   //if track has innermost momentum, check with innermost phi
0828   //if track also has valid vertex, propagate to vertex then calculate phi there
0829   //if track has no innermost momentum, just return true, because track->phi() makes no sense
0830   if (tk_l->extra().isNonnull() && tk_r->extra().isNonnull()) {
0831     double iphi1 = tk_l->innerMomentum().phi(), iphi2 = tk_r->innerMomentum().phi();
0832     if (vtx.isValid()) {
0833       PropagatorWithMaterial propag(anyDirection, 0.000511, &magField);
0834 
0835       double recoPhoR = vtx.position().Rho();
0836       Surface::RotationType rot;
0837       ReferenceCountingPointer<BoundCylinder> theBarrel_(new BoundCylinder(
0838           recoPhoR,
0839           Surface::PositionType(0, 0, 0),
0840           rot,
0841           new SimpleCylinderBounds(
0842               recoPhoR - 0.001, recoPhoR + 0.001, -fabs(vtx.position().z()), fabs(vtx.position().z()))));
0843       ReferenceCountingPointer<BoundDisk> theDisk_(new BoundDisk(
0844           Surface::PositionType(0, 0, vtx.position().z()), rot, new SimpleDiskBounds(0, recoPhoR, -0.001, 0.001)));
0845 
0846       const auto myTSOS1 = trajectoryStateTransform::innerStateOnSurface(*tk_l, trackerGeom, &magField);
0847       const auto myTSOS2 = trajectoryStateTransform::innerStateOnSurface(*tk_r, trackerGeom, &magField);
0848       TrajectoryStateOnSurface stateAtVtx1, stateAtVtx2;
0849       stateAtVtx1 = propag.propagate(myTSOS1, *theBarrel_);
0850       if (!stateAtVtx1.isValid()) {
0851         stateAtVtx1 = propag.propagate(myTSOS1, *theDisk_);
0852       }
0853       if (stateAtVtx1.isValid()) {
0854         iphi1 = stateAtVtx1.globalDirection().phi();
0855       }
0856       stateAtVtx2 = propag.propagate(myTSOS2, *theBarrel_);
0857       if (!stateAtVtx2.isValid()) {
0858         stateAtVtx2 = propag.propagate(myTSOS2, *theDisk_);
0859       }
0860       if (stateAtVtx2.isValid()) {
0861         iphi2 = stateAtVtx2.globalDirection().phi();
0862       }
0863     }
0864     const double dPhi = reco::deltaPhi(iphi1, iphi2);
0865     return (fabs(dPhi) < deltaPhi_);
0866   } else {
0867     return true;
0868   }
0869 }
0870 
0871 bool ConversionProducer::preselectTrackPair(const reco::TransientTrack& ttk_l,
0872                                             const reco::TransientTrack& ttk_r,
0873                                             double& appDist) {
0874   double dCotTheta = 1. / tan(ttk_l.track().innerMomentum().theta()) - 1. / tan(ttk_r.track().innerMomentum().theta());
0875   if (allowDeltaCot_ && (std::abs(dCotTheta) > deltaCotTheta_)) {
0876     return false;
0877   }
0878 
0879   //non-conversion hypothesis, reject prompt track pairs
0880   ClosestApproachInRPhi closest;
0881   closest.calculate(ttk_l.innermostMeasurementState(), ttk_r.innermostMeasurementState());
0882   if (!closest.status()) {
0883     return false;
0884   }
0885 
0886   if (closest.crossingPoint().perp() < r_cut) {
0887     return false;
0888   }
0889 
0890   //compute tangent point btw tracks (conversion hypothesis)
0891   TangentApproachInRPhi tangent;
0892   tangent.calculate(ttk_l.innermostMeasurementState(), ttk_r.innermostMeasurementState());
0893   if (!tangent.status()) {
0894     return false;
0895   }
0896 
0897   GlobalPoint tangentPoint = tangent.crossingPoint();
0898   double rho = tangentPoint.perp();
0899 
0900   //reject candidates well outside of tracker bounds
0901   if (rho > maxTrackRho_) {
0902     return false;
0903   }
0904 
0905   if (std::abs(tangentPoint.z()) > maxTrackZ_) {
0906     return false;
0907   }
0908 
0909   std::pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters> trajs = tangent.trajectoryParameters();
0910 
0911   //very large separation in z, no hope
0912   if (std::abs(trajs.first.position().z() - trajs.second.position().z()) > dzCut_) {
0913     return false;
0914   }
0915 
0916   float minApproach = tangent.perpdist();
0917   appDist = minApproach;
0918 
0919   if (allowMinApproach_ && (minApproach < minApproachLow_ || minApproach > minApproachHigh_)) {
0920     return false;
0921   }
0922 
0923   return true;
0924 }
0925 
0926 bool ConversionProducer::checkTrackPair(const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& ll,
0927                                         const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& rr) {
0928   const reco::CaloClusterPtr& bc_l = ll.second;  //can be null, so check isNonnull()
0929   const reco::CaloClusterPtr& bc_r = rr.second;
0930 
0931   //The cuts should be ordered by considering if takes time and if cuts off many fakes
0932   if (allowTrackBC_) {
0933     //check energy of BC
0934     double total_e_bc = 0;
0935     if (bc_l.isNonnull())
0936       total_e_bc += bc_l->energy();
0937     if (rightBC_)
0938       if (bc_r.isNonnull())
0939         total_e_bc += bc_r->energy();
0940 
0941     if (total_e_bc < energyTotalBC_)
0942       return false;
0943   }
0944 
0945   return true;
0946 }
0947 
0948 double ConversionProducer::etaTransformation(float EtaParticle, float Zvertex) {
0949   //---Definitions
0950   const float PI = 3.1415927;
0951 
0952   //---Definitions for ECAL
0953   const float R_ECAL = 136.5;
0954   const float Z_Endcap = 328.0;
0955   const float etaBarrelEndcap = 1.479;
0956 
0957   //---ETA correction
0958 
0959   float Theta = 0.0;
0960   float ZEcal = R_ECAL * sinh(EtaParticle) + Zvertex;
0961 
0962   if (ZEcal != 0.0)
0963     Theta = atan(R_ECAL / ZEcal);
0964   if (Theta < 0.0)
0965     Theta = Theta + PI;
0966   double ETA = -log(tan(0.5 * Theta));
0967 
0968   if (fabs(ETA) > etaBarrelEndcap) {
0969     float Zend = Z_Endcap;
0970     if (EtaParticle < 0.0)
0971       Zend = -Zend;
0972     float Zlen = Zend - Zvertex;
0973     float RR = Zlen / sinh(EtaParticle);
0974     Theta = atan(RR / Zend);
0975     if (Theta < 0.0)
0976       Theta = Theta + PI;
0977     ETA = -log(tan(0.5 * Theta));
0978   }
0979   //---Return the result
0980   return ETA;
0981   //---end
0982 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!