CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    HTrackAssociator
00004 // Class:      HTrackAssociator
00005 // 
00006 /*
00007 
00008  Description: <one line class summary>
00009 
00010  Implementation:
00011      <Notes on implementation>
00012 */
00013 //
00014 // Original Author:  Dmytro Kovalskyi
00015 // Modified for ECAL+HCAL by:  Michal Szleper
00016 //         Created:  Fri Apr 21 10:59:41 PDT 2006
00017 // $Id: TrackAssociator.cc,v 1.5 2008/05/27 16:36:33 heltsley Exp $
00018 //
00019 //
00020 
00021 #include "Calibration/Tools/interface/TrackAssociator.h"
00022 
00023 // user include files
00024 
00025 
00026 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00027 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00028 #include "DataFormats/DetId/interface/DetId.h"
00029 
00030 // calorimeter info
00031 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00032 
00033 
00034 
00035 
00036 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00037 
00038 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00039 #include "Utilities/Timing/interface/TimingReport.h"
00040 #include <stack>
00041 #include <set>
00042 
00043 
00044 #include "Calibration/Tools/interface/TimerStack.h"
00045 
00046 
00047 //
00048 // class declaration
00049 //
00050 
00051 using namespace reco;
00052 
00053 HTrackAssociator::HTrackAssociator() 
00054 {
00055    ivProp_ = 0;
00056    defProp_ = 0;
00057    debug_ = 0;
00058    caloTowerMap_ = 0;
00059    useDefaultPropagator_ = false;
00060 }
00061 
00062 HTrackAssociator::~HTrackAssociator()
00063 {
00064    if (defProp_) delete defProp_;
00065 }
00066 
00067 //-----------------------------------------------------------------------------
00068 void HTrackAssociator::addDataLabels( const std::string className,
00069                                      const std::string moduleLabel,
00070                                      const std::string productInstanceLabel)
00071 {
00072    if (className == "EBRecHitCollection")
00073      {
00074         EBRecHitCollectionLabels.clear();
00075         EBRecHitCollectionLabels.push_back(moduleLabel);
00076         EBRecHitCollectionLabels.push_back(productInstanceLabel);
00077      }
00078    if (className == "EERecHitCollection")
00079      {
00080         EERecHitCollectionLabels.clear();
00081         EERecHitCollectionLabels.push_back(moduleLabel);
00082         EERecHitCollectionLabels.push_back(productInstanceLabel);
00083      }
00084    if (className == "HBHERecHitCollection")
00085      {
00086         HBHERecHitCollectionLabels.clear();
00087         HBHERecHitCollectionLabels.push_back(moduleLabel);
00088         HBHERecHitCollectionLabels.push_back(productInstanceLabel);
00089      }
00090    if (className == "CaloTowerCollection")
00091      {
00092         CaloTowerCollectionLabels.clear();
00093         CaloTowerCollectionLabels.push_back(moduleLabel);
00094         CaloTowerCollectionLabels.push_back(productInstanceLabel);
00095      }
00096 }
00097 
00098 
00099 //-----------------------------------------------------------------------------
00100 void HTrackAssociator::setPropagator( Propagator* ptr)
00101 {
00102    ivProp_ = ptr; 
00103    caloDetIdAssociator_.setPropagator(ivProp_);
00104    ecalDetIdAssociator_.setPropagator(ivProp_);
00105    hcalDetIdAssociator_.setPropagator(ivProp_);
00106 }
00107 
00108 //-----------------------------------------------------------------------------
00109 void HTrackAssociator::useDefaultPropagator()
00110 {
00111    useDefaultPropagator_ = true;
00112 }
00113 
00114 
00115 //-----------------------------------------------------------------------------
00116 void HTrackAssociator::init( const edm::EventSetup& iSetup )
00117 {
00118    // access the calorimeter geometry
00119    iSetup.get<CaloGeometryRecord>().get(theCaloGeometry_);
00120    if (!theCaloGeometry_.isValid()) 
00121      throw cms::Exception("FatalError") << "Unable to find IdealGeometryRecord in event!\n";
00122    
00123    if (useDefaultPropagator_ && ! defProp_ ) {
00124       // setup propagator
00125       edm::ESHandle<MagneticField> bField;
00126       iSetup.get<IdealMagneticFieldRecord>().get(bField);
00127       
00128       SteppingHelixPropagator* prop  = new SteppingHelixPropagator(&*bField,anyDirection);
00129       prop->setMaterialMode(false);
00130       prop->applyRadX0Correction(true);
00131       // prop->setDebug(true); // tmp
00132       defProp_ = prop;
00133       setPropagator(defProp_);
00134    }
00135    
00136         
00137 }
00138 
00139 //-----------------------------------------------------------------------------
00140 HTrackDetMatchInfo HTrackAssociator::associate( const edm::Event& iEvent,
00141                                               const edm::EventSetup& iSetup,
00142                                               const FreeTrajectoryState& trackOrigin,
00143                                               const HAssociatorParameters& parameters )
00144 {
00145    HTrackDetMatchInfo info;
00146    using namespace edm;
00147    HTimerStack timers;
00148    
00149    init( iSetup );
00150    
00151    FreeTrajectoryState currentPosition(trackOrigin);
00152 
00153    if (parameters.useEcal) fillEcal( iEvent, iSetup, info, currentPosition,parameters.idREcal, parameters.dREcal);
00154    if (parameters.useHcal) fillHcal( iEvent, iSetup, info, currentPosition,parameters.idRHcal,parameters.dRHcal);
00155    if (parameters.useCalo) fillCaloTowers( iEvent, iSetup, info, currentPosition,parameters.idRCalo,parameters.dRCalo);
00156 
00157    return info;
00158 }
00159 
00160 //-----------------------------------------------------------------------------
00161 std::vector<EcalRecHit> HTrackAssociator::associateEcal( const edm::Event& iEvent,
00162                                                         const edm::EventSetup& iSetup,
00163                                                         const FreeTrajectoryState& trackOrigin,
00164                                                         const double dR )
00165 {
00166    HAssociatorParameters parameters;
00167    parameters.useHcal = false;
00168    parameters.dREcal = dR;
00169    HTrackDetMatchInfo info( associate(iEvent, iSetup, trackOrigin, parameters ));
00170    if (dR>0) 
00171      return info.coneEcalRecHits;
00172    else
00173      return info.crossedEcalRecHits;
00174 }
00175 
00176 //-----------------------------------------------------------------------------
00177 double HTrackAssociator::getEcalEnergy( const edm::Event& iEvent,
00178                                        const edm::EventSetup& iSetup,
00179                                        const FreeTrajectoryState& trackOrigin,
00180                                        const double dR )
00181 {
00182    HAssociatorParameters parameters;
00183    parameters.useHcal = false;
00184    parameters.dREcal = dR;
00185    HTrackDetMatchInfo info = associate(iEvent, iSetup, trackOrigin, parameters );
00186    if(dR>0) 
00187      return info.ecalConeEnergyFromRecHits();
00188    else
00189      return info.ecalEnergyFromRecHits();
00190 }
00191 
00192 //-----------------------------------------------------------------------------
00193 std::vector<CaloTower> HTrackAssociator::associateHcal( const edm::Event& iEvent,
00194                                                        const edm::EventSetup& iSetup,
00195                                                        const FreeTrajectoryState& trackOrigin,
00196                                                        const double dR )
00197 {
00198    HAssociatorParameters parameters;
00199    parameters.useEcal = false;
00200    parameters.dRHcal = dR;
00201    HTrackDetMatchInfo info( associate(iEvent, iSetup, trackOrigin, parameters ));
00202    if (dR>0) 
00203      return info.coneTowers;
00204    else
00205      return info.crossedTowers;
00206    
00207 }
00208 
00209 //-----------------------------------------------------------------------------
00210 double HTrackAssociator::getHcalEnergy( const edm::Event& iEvent,
00211                                        const edm::EventSetup& iSetup,
00212                                        const FreeTrajectoryState& trackOrigin,
00213                                        const double dR )
00214 {
00215    HAssociatorParameters parameters;
00216    parameters.useEcal = false;
00217    parameters.dRHcal = dR;
00218    HTrackDetMatchInfo info( associate(iEvent, iSetup, trackOrigin, parameters ));
00219    if (dR>0) 
00220      return info.hcalConeEnergyFromRecHits();
00221    else
00222      return info.hcalEnergyFromRecHits();
00223 }
00224 
00225 //-----------------------------------------------------------------------------
00226 void HTrackAssociator::fillEcal( const edm::Event& iEvent,
00227                                 const edm::EventSetup& iSetup,
00228                                 HTrackDetMatchInfo& info,
00229                                 const FreeTrajectoryState& trajectoryPoint,
00230                                 const int idR,
00231                                 const double dR)
00232 {
00233    HTimerStack timers;
00234    timers.push("HTrackAssociator::fillEcal");
00235    
00236    ecalDetIdAssociator_.setGeometry(&*theCaloGeometry_);
00237    
00238    timers.push("HTrackAssociator::fillEcal::propagation");
00239    // ECAL points (EB+EE)
00240    std::vector<GlobalPoint> ecalPoints;
00241    ecalPoints.push_back(GlobalPoint(135.,0,310.));
00242    ecalPoints.push_back(GlobalPoint(150.,0,340.));
00243    ecalPoints.push_back(GlobalPoint(170.,0,370.));
00244    
00245    std::vector<GlobalPoint> ecalTrajectory = ecalDetIdAssociator_.getTrajectory(trajectoryPoint, ecalPoints);
00246 //   if(ecalTrajectory.empty()) throw cms::Exception("FatalError") << "Failed to propagate a track to ECAL\n";
00247 
00248    if(ecalTrajectory.empty()) {
00249       LogTrace("HTrackAssociator::fillEcal") << "Failed to propagate a track to ECAL; moving on\n";
00250       info.isGoodEcal = 0;
00251       return;
00252    }
00253    
00254    info.isGoodEcal = 1;
00255 
00256    info.trkGlobPosAtEcal = getPoint(ecalTrajectory[0]);
00257 
00258    // Find ECAL crystals
00259    timers.pop_and_push("HTrackAssociator::fillEcal::access::EcalBarrel");
00260    edm::Handle<EBRecHitCollection> EBRecHits;
00261    edm::Handle<EERecHitCollection> EERecHits;
00262 //   if (EBRecHitCollectionLabels.empty())
00263 //     throw cms::Exception("FatalError") << "Module lable is not set for EBRecHitCollection.\n";
00264 //   else
00265      iEvent.getByLabel (EBRecHitCollectionLabels[0], EBRecHitCollectionLabels[1], EBRecHits);
00266 //   if (!EBRecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in event!\n";
00267      if (EERecHitCollectionLabels[1]!=EBRecHitCollectionLabels[1]) iEvent.getByLabel (EERecHitCollectionLabels[0], EERecHitCollectionLabels[1], EERecHits);
00268 //   if (!EERecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n";
00269 
00270    timers.pop_and_push("HTrackAssociator::fillEcal::matching");
00271    std::set<DetId> ecalIdsInRegion = ecalDetIdAssociator_.getDetIdsCloseToAPoint(ecalTrajectory[0],idR);
00272    std::set<DetId> ecalIdsInACone =  ecalDetIdAssociator_.getDetIdsInACone(ecalIdsInRegion, ecalTrajectory, dR);
00273    std::set<DetId> crossedEcalIds =  ecalDetIdAssociator_.getCrossedDetIds(ecalIdsInRegion, ecalTrajectory);
00274    
00275    // add EcalRecHits
00276    timers.pop_and_push("HTrackAssociator::fillEcal::addEcalRecHits");
00277    for(std::set<DetId>::const_iterator itr=crossedEcalIds.begin(); itr!=crossedEcalIds.end();itr++) {
00278      std::vector<EcalRecHit>::const_iterator hit = (*EBRecHits).find(*itr);
00279      if(hit != (*EBRecHits).end()) 
00280        info.crossedEcalRecHits.push_back(*hit);
00281      else  
00282        LogTrace("HTrackAssociator::fillEcal") << "EcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00283    }
00284    for(std::set<DetId>::const_iterator itr=ecalIdsInACone.begin(); itr!=ecalIdsInACone.end();itr++) {
00285      std::vector<EcalRecHit>::const_iterator hit = (*EBRecHits).find(*itr);
00286      if(hit != (*EBRecHits).end()) {
00287        info.coneEcalRecHits.push_back(*hit);
00288      }
00289      else 
00290        LogTrace("HTrackAssociator::fillEcal") << "EcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00291    }
00292    if (EERecHitCollectionLabels[1]==EBRecHitCollectionLabels[1])return;
00293    for(std::set<DetId>::const_iterator itr=crossedEcalIds.begin(); itr!=crossedEcalIds.end();itr++) {
00294      std::vector<EcalRecHit>::const_iterator hit = (*EERecHits).find(*itr);
00295      if(hit != (*EERecHits).end()) 
00296        info.crossedEcalRecHits.push_back(*hit);
00297      else  
00298        LogTrace("HTrackAssociator::fillEcal") << "EcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00299    }
00300    for(std::set<DetId>::const_iterator itr=ecalIdsInACone.begin(); itr!=ecalIdsInACone.end();itr++) {
00301      std::vector<EcalRecHit>::const_iterator hit = (*EERecHits).find(*itr);
00302      if(hit != (*EERecHits).end()) {
00303        info.coneEcalRecHits.push_back(*hit);
00304      }
00305      else 
00306        LogTrace("HTrackAssociator::fillEcal") << "EcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00307    }
00308 }
00309 
00310 //-----------------------------------------------------------------------------
00311 void HTrackAssociator::fillCaloTowers( const edm::Event& iEvent,
00312                                       const edm::EventSetup& iSetup,
00313                                       HTrackDetMatchInfo& info,
00314                                       const FreeTrajectoryState& trajectoryPoint,
00315                                       const int idR,
00316                                       const double dR)
00317 {
00318    // ECAL hits are not used for the CaloTower identification
00319    HTimerStack timers;
00320    timers.push("HTrackAssociator::fillCaloTowers");
00321 
00322    caloDetIdAssociator_.setGeometry(&*theCaloGeometry_);
00323    
00324    // HCAL points (HB+HE)
00325    timers.push("HTrackAssociator::fillCaloTowers::propagation");
00326    std::vector<GlobalPoint> hcalPoints;   
00327    hcalPoints.push_back(GlobalPoint(135.,0,310.));
00328    hcalPoints.push_back(GlobalPoint(150.,0,340.));
00329    hcalPoints.push_back(GlobalPoint(170.,0,370.));
00330    hcalPoints.push_back(GlobalPoint(190.,0,400.));
00331    hcalPoints.push_back(GlobalPoint(240.,0,500.));
00332    hcalPoints.push_back(GlobalPoint(280.,0,550.));
00333    
00334    std::vector<GlobalPoint> hcalTrajectory = caloDetIdAssociator_.getTrajectory(trajectoryPoint, hcalPoints);
00335 //   if(hcalTrajectory.empty()) throw cms::Exception("FatalError") << "Failed to propagate the track to HCAL\n";
00336 
00337    if(hcalTrajectory.empty()) {
00338       LogTrace("HTrackAssociator::fillEcal") << "Failed to propagate a track to ECAL; moving on\n";
00339       info.isGoodCalo = 0;
00340       info.isGoodEcal = 0;
00341       std::cout<<" HTrackAssociator::fillCaloTowers::Failed to propagate a track to ECAL "<<std::endl;
00342       return;
00343    }
00344    
00345    info.isGoodCalo = 1;
00346    info.isGoodEcal = 1;
00347    info.trkGlobPosAtEcal = getPoint(hcalTrajectory[0]);
00348    
00349    if(hcalTrajectory.size()<4) {
00350       LogTrace("HTrackAssociator::fillEcal") << "Failed to propagate a track to HCAL; moving on\n";
00351       info.isGoodHcal = 0;
00352    }
00353    
00354    info.isGoodHcal = 1;
00355    
00356    info.trkGlobPosAtHcal = getPoint(hcalTrajectory[4]);
00357 
00358    // find crossed CaloTowers
00359    timers.pop_and_push("HTrackAssociator::fillCaloTowers::access::CaloTowers");
00360    edm::Handle<CaloTowerCollection> caloTowers;
00361 
00362    if (CaloTowerCollectionLabels.empty())
00363      throw cms::Exception("FatalError") << "Module lable is not set for CaloTowers.\n";
00364    else
00365      iEvent.getByLabel (CaloTowerCollectionLabels[0], CaloTowerCollectionLabels[1], caloTowers);
00366    if (!caloTowers.isValid())  throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
00367    
00368    timers.push("HTrackAssociator::fillCaloTowers::matching");
00369 
00370 // first get DetIds in a predefined NxN region
00371 //   std::set<DetId> caloTowerIdsInBigRegion = caloDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR+1);
00372    std::set<DetId> caloTowerIdsInRegion = caloDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR);
00373 
00374    std::set<DetId> caloTowerIdsInACone;
00375    std::set<DetId> crossedCaloTowerIds;
00376    std::set<DetId> caloTowerIdsInBox;
00377    caloTowerIdsInACone = caloDetIdAssociator_.getDetIdsInACone(caloTowerIdsInRegion, hcalTrajectory, dR);
00378 // get DetId of the most energetic tower in that region
00379    crossedCaloTowerIds = caloDetIdAssociator_.getMaxEDetId(caloTowerIdsInRegion, caloTowers);
00380 // get DetIds of the towers surrounding the most energetic one
00381    caloTowerIdsInBox = caloDetIdAssociator_.getDetIdsInACone(crossedCaloTowerIds, hcalTrajectory, -1.);
00382 
00383 //
00384 //  Debug prints
00385 //
00386 //   std::cout <<" Debug printout in CaloTowers "<<std::endl;
00387 //   std::cout <<" with position at outer layer:r,z,phi "<<trajectoryPoint.position().eta()<<
00388 //               " "<<trajectoryPoint.position().phi()<<
00389 //             " "<<trajectoryPoint.position().perp()<<
00390 //               " "<<trajectoryPoint.position().z()<<
00391 //             " "<<trajectoryPoint.charge()<<std::endl;
00392 //   std::cout <<" Trajectory point at ECAL surface:eta:phi:radius:z "<<(hcalTrajectory[0]).eta()<<
00393 //               " "<<(hcalTrajectory[0]).phi()<<
00394 //               " "<<(hcalTrajectory[0]).perp()<<
00395 //             " "<<(hcalTrajectory[0]).z()<<
00396 //             " momentum "<<trajectoryPoint.momentum().perp()<<std::endl; 
00397 //
00398 //   std::cout<<" Number of towers in the region "<<caloTowerIdsInRegion.size()<<" idR= "<<idR<<std::endl;
00399    
00400    // add CaloTowers
00401    timers.push("HTrackAssociator::fillCaloTowers::addCaloTowers");
00402    for(std::set<DetId>::const_iterator itr=crossedCaloTowerIds.begin(); itr!=crossedCaloTowerIds.end();itr++)
00403      {
00404         DetId id(*itr);
00405         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00406         if(tower != (*caloTowers).end()) 
00407           info.crossedTowers.push_back(*tower);
00408         else
00409           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00410      }
00411 
00412    for(std::set<DetId>::const_iterator itr=caloTowerIdsInACone.begin(); itr!=caloTowerIdsInACone.end();itr++)
00413      {
00414         DetId id(*itr);
00415         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00416         if(tower != (*caloTowers).end()) {
00417           info.coneTowers.push_back(*tower);
00418         }
00419         else 
00420           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00421      }
00422 
00423    for(std::set<DetId>::const_iterator itr=caloTowerIdsInBox.begin(); itr!=caloTowerIdsInBox.end();itr++)
00424      {
00425         DetId id(*itr);
00426         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00427         if(tower != (*caloTowers).end()) {
00428           info.boxTowers.push_back(*tower);
00429         }
00430         else
00431           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00432      }
00433    
00434    for(std::set<DetId>::const_iterator itr=caloTowerIdsInRegion.begin(); itr!=caloTowerIdsInRegion.end();itr++)
00435      {
00436         DetId id(*itr);
00437         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00438         if(tower != (*caloTowers).end()) {
00439           info.regionTowers.push_back(*tower);
00440         }
00441         else 
00442           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00443      }
00444 
00445 }
00446 
00447 //-----------------------------------------------------------------------------
00448 void HTrackAssociator::fillHcal( const edm::Event& iEvent,
00449                                       const edm::EventSetup& iSetup,
00450                                       HTrackDetMatchInfo& info,
00451                                       const FreeTrajectoryState& trajectoryPoint,
00452                                       const int idR,
00453                                       const double dR) {
00454   HTimerStack timers;
00455   timers.push("HTrackAssociator::fillHcal");
00456 
00457   hcalDetIdAssociator_.setGeometry(&*theCaloGeometry_);
00458 
00459 // HCAL points (HB+HE)
00460   timers.push("HTrackAssociator::fillHcal::propagation");
00461   std::vector<GlobalPoint> hcalPoints;
00462   hcalPoints.push_back(GlobalPoint(190.,0,400.));
00463   hcalPoints.push_back(GlobalPoint(240.,0,500.));
00464   hcalPoints.push_back(GlobalPoint(280.,0,550.));
00465 
00466   std::vector<GlobalPoint> hcalTrajectory = hcalDetIdAssociator_.getTrajectory(trajectoryPoint, hcalPoints);
00467 //   if(hcalTrajectory.empty()) throw cms::Exception("FatalError") << "Failed to propagate the track to HCAL\n";
00468 
00469   if(hcalTrajectory.empty()) {
00470     LogTrace("HTrackAssociator::fillHcal") << "Failed to propagate a track to HCAL; moving on\n";
00471     info.isGoodHcal = 0;
00472 //    std::cout<<" HTrackAssociator::fillHcal::Failed to propagate a track to HCAL "<<std::endl;
00473     return;
00474   }
00475 
00476   info.isGoodHcal = 1;
00477 
00478   info.trkGlobPosAtHcal = getPoint(hcalTrajectory[0]);
00479   timers.pop_and_push("HTrackAssociator::fillHcal::access::Hcal");
00480 
00481   edm::Handle<HBHERecHitCollection> HBHERecHits;
00482 //  if (HBHERecHitCollectionLabels.empty())
00483 //    throw cms::Exception("FatalError") << "Module label is not set for HBHERecHitCollection.\n";
00484 //  else
00485     iEvent.getByLabel (HBHERecHitCollectionLabels[0], HBHERecHitCollectionLabels[1], HBHERecHits);
00486   if (!HBHERecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find HBHERecHitCollection in event!\n";
00487 
00488   timers.pop_and_push("HTrackAssociator::fillHcal::matching");
00489 
00490 // first get DetIds in a predefined NxN region
00491 //  std::set<DetId> hcalIdsInBigRegion = hcalDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR+1);
00492   std::set<DetId> hcalIdsInRegion = hcalDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR);
00493 
00494   std::set<DetId> hcalIdsInACone;
00495   std::set<DetId> crossedHcalIds;
00496   std::set<DetId> hcalIdsInBox;
00497   hcalIdsInACone = hcalDetIdAssociator_.getDetIdsInACone(hcalIdsInRegion, hcalTrajectory, dR);
00498 // get DetId of the most energetic tower in that region
00499   crossedHcalIds = hcalDetIdAssociator_.getMaxEDetId(hcalIdsInRegion, HBHERecHits);
00500 // get DetIds of the towers surrounding the most energetic one
00501   hcalIdsInBox = hcalDetIdAssociator_.getDetIdsInACone(crossedHcalIds, hcalTrajectory, -1.);
00502 
00503 // add HcalRecHits
00504   timers.pop_and_push("HTrackAssociator::fillHcal::addHcalRecHits");
00505   for(std::set<DetId>::const_iterator itr=crossedHcalIds.begin(); itr!=crossedHcalIds.end();itr++) {
00506     std::vector<HBHERecHit>::const_iterator hit = (*HBHERecHits).find(*itr);
00507     if(hit != (*HBHERecHits).end())
00508       info.crossedHcalRecHits.push_back(*hit);
00509     else
00510       LogTrace("HTrackAssociator::fillHcal") << "HcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00511   }
00512   for(std::set<DetId>::const_iterator itr=hcalIdsInACone.begin(); itr!=hcalIdsInACone.end();itr++) {
00513     std::vector<HBHERecHit>::const_iterator hit = (*HBHERecHits).find(*itr);
00514     if(hit != (*HBHERecHits).end())
00515       info.coneHcalRecHits.push_back(*hit);
00516     else
00517       LogTrace("HTrackAssociator::fillHcal") << "HcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00518   }
00519   for(std::set<DetId>::const_iterator itr=hcalIdsInBox.begin(); itr!=hcalIdsInBox.end();itr++) {
00520     std::vector<HBHERecHit>::const_iterator hit = (*HBHERecHits).find(*itr);
00521     if(hit != (*HBHERecHits).end())
00522       info.boxHcalRecHits.push_back(*hit);
00523     else
00524       LogTrace("HTrackAssociator::fillHcal") << "HcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00525   }
00526   for(std::set<DetId>::const_iterator itr=hcalIdsInRegion.begin(); itr!=hcalIdsInRegion.end();itr++) {
00527     std::vector<HBHERecHit>::const_iterator hit = (*HBHERecHits).find(*itr);
00528     if(hit != (*HBHERecHits).end())
00529       info.regionHcalRecHits.push_back(*hit);
00530     else
00531       LogTrace("HTrackAssociator::fillHcal") << "HcalRecHit is not found for DetId: " << itr->rawId() <<"\n";
00532   }
00533 }
00534 
00535 //-----------------------------------------------------------------------------
00536 void HTrackAssociator::fillHcalTowers( const edm::Event& iEvent,
00537                                       const edm::EventSetup& iSetup,
00538                                       HTrackDetMatchInfo& info,
00539                                       const FreeTrajectoryState& trajectoryPoint,
00540                                       const int idR,
00541                                       const double dR)
00542 {
00543    // ECAL hits are not used for the CaloTower identification
00544    HTimerStack timers;
00545    timers.push("HTrackAssociator::fillCaloTowers");
00546 
00547    caloDetIdAssociator_.setGeometry(&*theCaloGeometry_);
00548    
00549    // HCAL points (HB+HE)
00550    timers.push("HTrackAssociator::fillCaloTowers::propagation");
00551    std::vector<GlobalPoint> hcalPoints;
00552    hcalPoints.push_back(GlobalPoint(190.,0,400.));
00553    hcalPoints.push_back(GlobalPoint(240.,0,500.));
00554    hcalPoints.push_back(GlobalPoint(280.,0,550.));
00555    
00556    std::vector<GlobalPoint> hcalTrajectory = caloDetIdAssociator_.getTrajectory(trajectoryPoint, hcalPoints);
00557 //   if(hcalTrajectory.empty()) throw cms::Exception("FatalError") << "Failed to propagate the track to HCAL\n";
00558 
00559    if(hcalTrajectory.empty()) {
00560       LogTrace("HTrackAssociator::fillEcal") << "Failed to propagate a track to HCAL; moving on\n";
00561       info.isGoodCalo = 0;
00562       std::cout<<" HTrackAssociator::fillCaloTowers::Failed to propagate a track to HCAL "<<std::endl;
00563       return;
00564    }
00565    
00566    info.isGoodCalo = 1;
00567 
00568    info.trkGlobPosAtHcal = getPoint(hcalTrajectory[0]);
00569    
00570    // find crossed CaloTowers
00571    timers.pop_and_push("HTrackAssociator::fillCaloTowers::access::CaloTowers");
00572    edm::Handle<CaloTowerCollection> caloTowers;
00573 
00574    if (CaloTowerCollectionLabels.empty())
00575      throw cms::Exception("FatalError") << "Module lable is not set for CaloTowers.\n";
00576    else
00577      iEvent.getByLabel (CaloTowerCollectionLabels[0], CaloTowerCollectionLabels[1], caloTowers);
00578    if (!caloTowers.isValid())  throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
00579    
00580    timers.push("HTrackAssociator::fillCaloTowers::matching");
00581 
00582 // first get DetIds in a predefined NxN region
00583    std::set<DetId> caloTowerIdsInBigRegion = caloDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR+1);
00584    std::set<DetId> caloTowerIdsInRegion = caloDetIdAssociator_.getDetIdsCloseToAPoint(hcalTrajectory[0],idR);
00585 
00586    std::set<DetId> caloTowerIdsInACone;
00587    std::set<DetId> crossedCaloTowerIds;
00588    std::set<DetId> caloTowerIdsInBox;
00589    caloTowerIdsInACone = caloDetIdAssociator_.getDetIdsInACone(caloTowerIdsInBigRegion, hcalTrajectory, dR);
00590 // get DetId of the most energetic tower in that region
00591    crossedCaloTowerIds = caloDetIdAssociator_.getMaxEDetId(caloTowerIdsInRegion, caloTowers);
00592 // get DetIds of the towers surrounding the most energetic one
00593    caloTowerIdsInBox = caloDetIdAssociator_.getDetIdsInACone(crossedCaloTowerIds, hcalTrajectory, -1.);
00594    
00595    // add CaloTowers
00596    timers.push("HTrackAssociator::fillCaloTowers::addCaloTowers");
00597    for(std::set<DetId>::const_iterator itr=crossedCaloTowerIds.begin(); itr!=crossedCaloTowerIds.end();itr++)
00598      {
00599         DetId id(*itr);
00600         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00601         if(tower != (*caloTowers).end()) 
00602           info.crossedTowers.push_back(*tower);
00603         else
00604           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00605      }
00606 
00607    for(std::set<DetId>::const_iterator itr=caloTowerIdsInACone.begin(); itr!=caloTowerIdsInACone.end();itr++)
00608      {
00609         DetId id(*itr);
00610         CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00611         if(tower != (*caloTowers).end()) 
00612           info.coneTowers.push_back(*tower);
00613         else 
00614           LogTrace("HTrackAssociator::fillEcal") << "CaloTower is not found for DetId: " << id.rawId() << "\n";
00615      }
00616    
00617 }
00618 
00619 //-----------------------------------------------------------------------------
00620 FreeTrajectoryState HTrackAssociator::getFreeTrajectoryState( const edm::EventSetup& iSetup, 
00621                                                              const SimTrack& track, 
00622                                                              const SimVertex& vertex )
00623 {
00624    edm::ESHandle<MagneticField> bField;
00625    iSetup.get<IdealMagneticFieldRecord>().get(bField);
00626    
00627    GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );
00628    // convert mm to cm
00629    GlobalPoint point( vertex.position().x()*.1, vertex.position().y()*.1, vertex.position().z()*.1 );
00630    int charge = track.type( )> 0 ? -1 : 1;
00631    GlobalTrajectoryParameters tPars(point, vector, charge, &*bField);
00632    
00633    HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
00634    CartesianTrajectoryError tCov(covT);
00635    
00636    return FreeTrajectoryState(tPars, tCov);
00637 }
00638 
00639 //-----------------------------------------------------------------------------
00640 FreeTrajectoryState HTrackAssociator::getFreeTrajectoryState( const edm::EventSetup& iSetup,
00641                                                              const reco::Track& track )
00642 {
00643    edm::ESHandle<MagneticField> bField;
00644    iSetup.get<IdealMagneticFieldRecord>().get(bField);
00645    
00646    GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );
00647 
00648    GlobalPoint point( track.vertex().x(), track.vertex().y(),  track.vertex().z() );
00649 
00650    GlobalTrajectoryParameters tPars(point, vector, track.charge(), &*bField);
00651    
00652    // FIX THIS !!!
00653    // need to convert from perigee to global or helix (curvilinear) frame
00654    // for now just an arbitrary matrix.
00655    HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
00656    CartesianTrajectoryError tCov(covT);
00657    
00658    return FreeTrajectoryState(tPars, tCov);
00659 }
00660