Reference Manual

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