00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004
00005
00006 #include "FWCore/Framework/interface/Event.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "FWCore/Framework/interface/ESHandle.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012
00013 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00014 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00015 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00016 #include "DataFormats/EgammaTrackReco/interface/TrackCandidateCaloClusterAssociation.h"
00017
00018 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00019
00020 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00021
00022 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionSeedFinder.h"
00023 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackFinder.h"
00024
00025 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00026 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00027 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
00028
00029 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionSeedFinder.h"
00030 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionSeedFinder.h"
00031 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionTrackFinder.h"
00032 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionTrackFinder.h"
00033 #include "Math/GenVector/VectorUtil.h"
00034
00035 #include "RecoEgamma/EgammaPhotonProducers/interface/SoftConversionTrackCandidateProducer.h"
00036
00037 bool IsGoodSeed(const TrajectorySeedCollection& seeds, const TrajectorySeed& seed){
00038
00039
00040
00041
00042
00043 bool found = false;
00044 for(TrajectorySeedCollection::const_iterator it = seeds.begin(); it != seeds.end(); it++){
00045 if(it->nHits() != seed.nHits()) continue;
00046 if(it->startingState().detId() != seed.startingState().detId()) continue;
00047 if(it->startingState().surfaceSide() != seed.startingState().surfaceSide()) continue;
00048 if((it->startingState().parameters().position() - seed.startingState().parameters().position()).mag() > 1.0e-6) continue;
00049 if((it->startingState().parameters().momentum() - seed.startingState().parameters().momentum()).mag() > 1.0e-6) continue;
00050 found = true;
00051 break;
00052 }
00053
00054 return found;
00055 }
00056
00057 SoftConversionTrackCandidateProducer::SoftConversionTrackCandidateProducer(const edm::ParameterSet& config) :
00058 conf_(config),
00059 theNavigationSchool_(0),
00060 theOutInSeedFinder_(0),
00061 theOutInTrackFinder_(0),
00062 theInOutSeedFinder_(0),
00063 theInOutTrackFinder_(0)
00064 {
00065 LogDebug("SoftConversionTrackCandidateProducer") << "SoftConversionTrackCandidateProducer CTOR " << "\n";
00066
00067 clusterType_ = conf_.getParameter<std::string>("clusterType");
00068 clusterBarrelCollection_ = conf_.getParameter<edm::InputTag>("clusterBarrelCollection");
00069 clusterEndcapCollection_ = conf_.getParameter<edm::InputTag>("clusterEndcapCollection");
00070
00071 OutInTrackCandidateCollection_ = conf_.getParameter<std::string>("outInTrackCandidateCollection");
00072 InOutTrackCandidateCollection_ = conf_.getParameter<std::string>("inOutTrackCandidateCollection");
00073
00074 OutInTrackClusterAssociationCollection_ = conf_.getParameter<std::string>("outInTrackCandidateClusterAssociationCollection");
00075 InOutTrackClusterAssociationCollection_ = conf_.getParameter<std::string>("inOutTrackCandidateClusterAssociationCollection");
00076
00077
00078 produces< TrackCandidateCollection > (OutInTrackCandidateCollection_);
00079 produces< TrackCandidateCollection > (InOutTrackCandidateCollection_);
00080 produces< reco::TrackCandidateCaloClusterPtrAssociation > ( OutInTrackClusterAssociationCollection_);
00081 produces< reco::TrackCandidateCaloClusterPtrAssociation > ( InOutTrackClusterAssociationCollection_);
00082
00083 }
00084
00085 SoftConversionTrackCandidateProducer::~SoftConversionTrackCandidateProducer() {}
00086
00087 void SoftConversionTrackCandidateProducer::setEventSetup (const edm::EventSetup & theEventSetup) {
00088 theOutInSeedFinder_->setEventSetup(theEventSetup);
00089 theInOutSeedFinder_->setEventSetup(theEventSetup);
00090 theOutInTrackFinder_->setEventSetup(theEventSetup);
00091 theInOutTrackFinder_->setEventSetup(theEventSetup);
00092 }
00093
00094
00095
00096 void SoftConversionTrackCandidateProducer::beginRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
00097 nEvt_=0;
00098
00099 edm::LogInfo("SoftConversionTrackCandidateProducer") << " get magnetic field" << "\n";
00100
00101 edm::ESHandle<NavigationSchool> nav;
00102 theEventSetup.get<NavigationSchoolRecord>().get("SimpleNavigationSchool", nav);
00103 theNavigationSchool_ = nav.product();
00104
00105
00106 edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the OutInSeedFinder" << "\n";
00107 theOutInSeedFinder_ = new OutInConversionSeedFinder ( conf_ );
00108
00109
00110 edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the OutInTrackFinder" << "\n";
00111 theOutInTrackFinder_ = new OutInConversionTrackFinder ( theEventSetup, conf_ );
00112
00113
00114 edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the InOutSeedFinder" << "\n";
00115 theInOutSeedFinder_ = new InOutConversionSeedFinder ( conf_ );
00116
00117
00118 edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the InOutTrackFinder" << "\n";
00119 theInOutTrackFinder_ = new InOutConversionTrackFinder ( theEventSetup, conf_ );
00120 }
00121
00122 void SoftConversionTrackCandidateProducer::endRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
00123 delete theOutInSeedFinder_;
00124 delete theOutInTrackFinder_;
00125 delete theInOutSeedFinder_;
00126 delete theInOutTrackFinder_;
00127 }
00128
00129
00130 void SoftConversionTrackCandidateProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
00131
00132 using namespace edm;
00133 nEvt_++;
00134 edm::LogInfo("SoftConversionTrackCandidateProducer") << "SoftConversionTrackCandidateProducer Analyzing event number: " << theEvent.id() << " Global Counter " << nEvt_ << "\n";
00135
00136 std::cout << "SoftConversionTrackCandidateProducer Analyzing event number: " << theEvent.id() << " Global Counter " << nEvt_ << "\n";
00137
00138 setEventSetup( theEventSetup );
00139 theOutInSeedFinder_->setEvent(theEvent);
00140 theInOutSeedFinder_->setEvent(theEvent);
00141 theOutInTrackFinder_->setEvent(theEvent);
00142 theInOutTrackFinder_->setEvent(theEvent);
00143
00144
00145 NavigationSetter setter(*theNavigationSchool_);
00146
00147
00148 std::auto_ptr<TrackCandidateCollection> outInTrackCandidate_p(new TrackCandidateCollection);
00149 std::auto_ptr<TrackCandidateCollection> inOutTrackCandidate_p(new TrackCandidateCollection);
00150 std::auto_ptr<reco::TrackCandidateCaloClusterPtrAssociation> outInAssoc_p(new reco::TrackCandidateCaloClusterPtrAssociation);
00151 std::auto_ptr<reco::TrackCandidateCaloClusterPtrAssociation> inOutAssoc_p(new reco::TrackCandidateCaloClusterPtrAssociation);
00152
00153 std::vector<edm::Ptr<reco::CaloCluster> > vecOfClusterRefForOutIn;
00154 std::vector<edm::Ptr<reco::CaloCluster> > vecOfClusterRefForInOut;
00155
00156
00157 edm::Handle<edm::View<reco::CaloCluster> > clusterBarrelHandle;
00158 theEvent.getByLabel(clusterBarrelCollection_, clusterBarrelHandle);
00159
00160 if (!clusterBarrelHandle.isValid()) {
00161 edm::LogError("SoftConverionTrackCandidateProducer") << "Error! Can't get the product "<<clusterBarrelCollection_.label();
00162 return;
00163 }
00164
00165 buildCollections(clusterBarrelHandle, *outInTrackCandidate_p, *inOutTrackCandidate_p, vecOfClusterRefForOutIn, vecOfClusterRefForInOut);
00166
00167 if(clusterType_ == "BasicCluster" ) {
00168
00169 edm::Handle<edm::View<reco::CaloCluster> > clusterEndcapHandle;
00170 theEvent.getByLabel(clusterEndcapCollection_, clusterEndcapHandle);
00171
00172 if (!clusterEndcapHandle.isValid()) {
00173 edm::LogError("SoftConversionTrackCandidateProducer") << "Error! Can't get the product "<<clusterEndcapCollection_.label();
00174 return;
00175 }
00176
00177 buildCollections(clusterEndcapHandle, *outInTrackCandidate_p, *inOutTrackCandidate_p, vecOfClusterRefForOutIn, vecOfClusterRefForInOut);
00178
00179 }
00180
00181
00182 const edm::OrphanHandle<TrackCandidateCollection> refprodOutInTrackC = theEvent.put( outInTrackCandidate_p, OutInTrackCandidateCollection_ );
00183 const edm::OrphanHandle<TrackCandidateCollection> refprodInOutTrackC = theEvent.put( inOutTrackCandidate_p, InOutTrackCandidateCollection_ );
00184
00185 edm::ValueMap<reco::CaloClusterPtr>::Filler fillerOI(*outInAssoc_p);
00186 fillerOI.insert(refprodOutInTrackC, vecOfClusterRefForOutIn.begin(), vecOfClusterRefForOutIn.end());
00187 fillerOI.fill();
00188 edm::ValueMap<reco::CaloClusterPtr>::Filler fillerIO(*inOutAssoc_p);
00189 fillerIO.insert(refprodInOutTrackC, vecOfClusterRefForInOut.begin(), vecOfClusterRefForInOut.end());
00190 fillerIO.fill();
00191
00192 theEvent.put( outInAssoc_p, OutInTrackClusterAssociationCollection_);
00193 theEvent.put( inOutAssoc_p, InOutTrackClusterAssociationCollection_);
00194 }
00195
00196
00197 void SoftConversionTrackCandidateProducer::buildCollections(const edm::Handle<edm::View<reco::CaloCluster> > & clusterHandle,
00198 TrackCandidateCollection& outInTrackCandidates,
00199 TrackCandidateCollection& inOutTrackCandidates,
00200 std::vector<edm::Ptr<reco::CaloCluster> >& vecRecOI,
00201 std::vector<edm::Ptr<reco::CaloCluster> >& vecRecIO) {
00202
00203
00204 TrackCandidateCollection tempTCC;
00205 TrajectorySeedCollection totalOISeeds;
00206 TrajectorySeedCollection totalIOSeeds;
00207
00208 int nClusters = (int) clusterHandle->size();
00209
00210
00211
00212 for(int iCluster=0; iCluster<nClusters; iCluster++){
00213 reco::CaloClusterPtr clusterRefOutIn = clusterHandle->ptrAt(iCluster);
00214 math::XYZPoint position = clusterRefOutIn->position();
00215 GlobalPoint gp(position.x(),position.y(),position.z());
00216 theOutInSeedFinder_->setCandidate(clusterRefOutIn->energy(),gp);
00217 theOutInSeedFinder_->makeSeeds(clusterRefOutIn);
00218
00219
00220 TrajectorySeedCollection oISeeds = theOutInSeedFinder_->seeds();
00221 for(TrajectorySeedCollection::const_iterator it = oISeeds.begin(); it != oISeeds.end(); it++){
00222 totalOISeeds.push_back(*it);
00223 }
00224
00225 std::vector<Trajectory> theOutInTracks= theOutInTrackFinder_->tracks(oISeeds, tempTCC);
00226 tempTCC.clear();
00227
00228 for(int jCluster=iCluster; jCluster<nClusters; jCluster++){
00229 reco::CaloClusterPtr clusterRefInOut = clusterHandle->ptrAt(jCluster);
00230
00231 math::XYZPoint position2 = clusterRefInOut->position();
00232 GlobalPoint gp2(position2.x(),position2.y(),position2.z());
00233 double dEta = std::abs(position.Eta() - position2.Eta());
00234 if(dEta > 0.1) continue;
00235
00236 double dPhi = std::abs(ROOT::Math::VectorUtil::DeltaPhi(position, position2));
00237 if(dPhi > 0.5) continue;
00238
00239 theInOutSeedFinder_->setCandidate(clusterRefInOut->energy(),gp2);
00240 theInOutSeedFinder_->setTracks(theOutInTracks);
00241 theInOutSeedFinder_->makeSeeds(clusterHandle);
00242
00243 TrajectorySeedCollection iOSeeds = theInOutSeedFinder_->seeds();
00244
00245 for(TrajectorySeedCollection::const_iterator it = iOSeeds.begin(); it != iOSeeds.end(); it++){
00246 totalIOSeeds.push_back(*it);
00247 }
00248
00249 }
00250 }
00251
00252
00253
00254 TrajectorySeedCollection oIFilteredSeeds;
00255 TrajectorySeedCollection iOFilteredSeeds;
00256
00257 tempTCC.clear();
00258 std::vector<Trajectory> tempTrj = theOutInTrackFinder_->tracks(totalOISeeds,tempTCC);
00259 for(std::vector<Trajectory>::iterator it = tempTrj.begin(); it!= tempTrj.end(); it++){
00260 oIFilteredSeeds.push_back(it->seed());
00261 }
00262
00263 tempTrj.clear();
00264 tempTCC.clear();
00265 tempTrj = theInOutTrackFinder_->tracks(totalIOSeeds,tempTCC);
00266 for(std::vector<Trajectory>::iterator it = tempTrj.begin(); it!= tempTrj.end(); it++){
00267 iOFilteredSeeds.push_back(it->seed());
00268 }
00269
00270 tempTCC.clear();
00271 tempTrj.clear();
00272 totalOISeeds.clear();
00273 totalIOSeeds.clear();
00274
00275
00276
00277
00278 for(int iCluster=0; iCluster<nClusters; iCluster++){
00279 reco::CaloClusterPtr clusterRefOutIn = clusterHandle->ptrAt(iCluster);
00280 math::XYZPoint position = clusterRefOutIn->position();
00281 GlobalPoint gp(position.x(),position.y(),position.z());
00282 theOutInSeedFinder_->setCandidate(clusterRefOutIn->energy(),gp);
00283 theOutInSeedFinder_->makeSeeds(clusterRefOutIn);
00284
00285 TrajectorySeedCollection oISeeds_all = theOutInSeedFinder_->seeds();
00286 TrajectorySeedCollection oISeeds;
00287 for(TrajectorySeedCollection::iterator it = oISeeds_all.begin(); it != oISeeds_all.end(); it++){
00288 if(IsGoodSeed(oIFilteredSeeds,*it)) oISeeds.push_back(*it);
00289 }
00290
00291 if(oISeeds.size() == 0) continue;
00292
00293 std::vector<Trajectory> theOutInTracks= theOutInTrackFinder_->tracks(oISeeds, outInTrackCandidates);
00294
00295 int nOITrj = (int) theOutInTracks.size();
00296 for(int itrj=0; itrj < nOITrj; itrj++) vecRecOI.push_back( clusterRefOutIn );
00297
00298 for(int jCluster=iCluster; jCluster<nClusters; jCluster++){
00299 reco::CaloClusterPtr clusterRefInOut = clusterHandle->ptrAt(jCluster);
00300
00301 math::XYZPoint position2 = clusterRefInOut->position();
00302 GlobalPoint gp2(position2.x(),position2.y(),position2.z());
00303 double dEta = std::abs(position.Eta() - position2.Eta());
00304 if(dEta > 0.1) continue;
00305
00306 double dPhi = std::abs(ROOT::Math::VectorUtil::DeltaPhi(position, position2));
00307 if(dPhi > 0.5) continue;
00308
00309 theInOutSeedFinder_->setCandidate(clusterRefInOut->energy(),gp2);
00310 theInOutSeedFinder_->setTracks(theOutInTracks);
00311 theInOutSeedFinder_->makeSeeds(clusterHandle);
00312
00313 TrajectorySeedCollection iOSeeds_all = theInOutSeedFinder_->seeds();
00314 TrajectorySeedCollection iOSeeds;
00315 for(TrajectorySeedCollection::iterator it = iOSeeds_all.begin(); it != iOSeeds_all.end(); it++){
00316 if(IsGoodSeed(iOFilteredSeeds,*it)) iOSeeds.push_back(*it);
00317 }
00318
00319 if(iOSeeds.size() == 0) continue;
00320
00321 std::vector<Trajectory> theInOutTracks= theInOutTrackFinder_->tracks(iOSeeds, inOutTrackCandidates);
00322
00323 int nIOTrj = (int) theInOutTracks.size();
00324 for(int itrj=0; itrj < nIOTrj; itrj++) vecRecIO.push_back( clusterRefInOut );
00325
00326 }
00327 }
00328
00329
00330 }
00331