CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoEgamma/EgammaPhotonProducers/src/SoftConversionTrackCandidateProducer.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004 
00005 // Framework
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 //  Abstract classes for the conversion tracking components
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 // Class header file
00035 #include "RecoEgamma/EgammaPhotonProducers/interface/SoftConversionTrackCandidateProducer.h"
00036 
00037 bool IsGoodSeed(const TrajectorySeedCollection& seeds, const TrajectorySeed& seed){
00038 
00039   // This function is not satisfactory. I don't know how to check equality of TrajectorySeed
00040   // So I compare all possible quantities in them.
00041   // This function can be dropped when I find to check equality of TrajectorySeed.
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   // Register the product
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   //get magnetic field
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   // get the Out In Seed Finder  
00106   edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the OutInSeedFinder" << "\n";
00107   theOutInSeedFinder_ = new OutInConversionSeedFinder (  conf_ );
00108   
00109   // get the Out In Track Finder
00110   edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the OutInTrackFinder" << "\n";
00111   theOutInTrackFinder_ = new OutInConversionTrackFinder ( theEventSetup, conf_  );
00112   
00113   // get the In Out Seed Finder  
00114   edm::LogInfo("SoftConversionTrackCandidateProducer") << " get the InOutSeedFinder" << "\n";
00115   theInOutSeedFinder_ = new InOutConversionSeedFinder ( conf_ );
00116   
00117   // get the In Out Track Finder
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   // Set the navigation school  
00145   NavigationSetter setter(*theNavigationSchool_);  
00146 
00147   // collections to be stored in events
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   // Get the basic cluster collection in the Barrel 
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     // Get the basic cluster collection in the Endcap 
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   // put all products in the event
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   // temporary collection
00204   TrackCandidateCollection tempTCC;
00205   TrajectorySeedCollection totalOISeeds; // total number of out-in trajectory seeds through entire cluster collection loop
00206   TrajectorySeedCollection totalIOSeeds; // total number of in-out trajectory seeds through entire cluster collection loop
00207 
00208   int nClusters = (int) clusterHandle->size();
00209 
00210   // first loop to fill totalOISeeds and totalIOSeeds
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     }// for jCluster
00250   }// for iCluster
00251 
00252 
00253   // Now we have total OI/IO seeds. Let's clean them up and save them with only giving good trajectories
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   // Now start normal procedure and consider seeds that belong to filtered ones.
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     }// for jCluster
00327   }// for iCluster
00328 
00329 
00330 }
00331