BetaCalculatorECAL Class Reference

#include <BetaCalculatorECAL.h>

Public Member Functions
void	addInfoToCandidate (susybsm::HSCParticle &candidate, edm::Handle< reco::TrackCollection > &tracks, edm::Event &iEvent, const edm::EventSetup &iSetup, susybsm::HSCPCaloInfo &caloInfo)
	BetaCalculatorECAL (const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)

Private Member Functions
void	addStepToXtal (std::map< int, GlobalPoint > &trackExitPositionMap, std::map< int, float > &trackCrossedXtalMap, DetId aDetId, float step, GlobalPoint point, const CaloSubdetectorGeometry *theSubdetGeometry)
std::vector< SteppingHelixStateInfo >	calcEcalDeposit (const FreeTrajectoryState *tkInnerState, const DetIdAssociator &associator)
int	getDetailedTrackLengthInXtals (std::map< int, GlobalPoint > &trackExitPositionMap, std::map< int, float > &trackCrossedXtalMap, double &totalLengthCurved, GlobalPoint &internalPointCurved, GlobalPoint &externalPointCurved, const CaloGeometry *theGeometry, const CaloTopology *theTopology, const std::vector< SteppingHelixStateInfo > &neckLace)

Private Attributes
edm::ESHandle< MagneticField >	bField_
edm::EDGetTokenT< EBRecHitCollection >	EBRecHitCollectionToken_
edm::ESHandle< DetIdAssociator >	ecalDetIdAssociator_
edm::EDGetTokenT< EERecHitCollection >	EERecHitCollectionToken_
TrackAssociatorParameters	parameters_
edm::ESHandle< CaloGeometry >	theCaloGeometry_
TrackDetectorAssociator	trackAssociator_

Detailed Description

Definition at line 33 of file BetaCalculatorECAL.h.

Constructor & Destructor Documentation

BetaCalculatorECAL::BetaCalculatorECAL	(	const edm::ParameterSet &	iConfig,
		edm::ConsumesCollector &&	iC
	)

Definition at line 18 of file BetaCalculatorECAL.cc.

References edm::ParameterSet::getParameter(), TrackAssociatorParameters::loadParameters(), parameters_, trackAssociator_, and TrackDetectorAssociator::useDefaultPropagator().

                                                                                                :
  EBRecHitCollectionToken_(iC.consumes<EBRecHitCollection>(iConfig.getParameter<edm::InputTag>("EBRecHitCollection"))),
  EERecHitCollectionToken_(iC.consumes<EERecHitCollection>(iConfig.getParameter<edm::InputTag>("EERecHitCollection")))
{
   edm::ParameterSet trkParameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
   parameters_.loadParameters( trkParameters, iC );
   trackAssociator_.useDefaultPropagator();

}

Member Function Documentation

void BetaCalculatorECAL::addInfoToCandidate	(	susybsm::HSCParticle &	candidate,
		edm::Handle< reco::TrackCollection > &	tracks,
		edm::Event &	iEvent,
		const edm::EventSetup &	iSetup,
		susybsm::HSCPCaloInfo &	caloInfo
	)

Definition at line 29 of file BetaCalculatorECAL.cc.

References TrackDetectorAssociator::associate(), bField_, calcEcalDeposit(), EcalRecHit::chi2(), TrackDetMatchInfo::crossedEnergy(), TrackDetMatchInfo::crossedHcalRecHits, runTauDisplay::dr, EBRecHitCollectionToken_, susybsm::HSCPCaloInfo::ecal3by3dir, susybsm::HSCPCaloInfo::ecal5by5dir, EcalBarrel, susybsm::HSCPCaloInfo::ecalBeta, susybsm::HSCPCaloInfo::ecalBetaError, susybsm::HSCPCaloInfo::ecalChi2s, susybsm::HSCPCaloInfo::ecalCrossedEnergy, susybsm::HSCPCaloInfo::ecalCrysCrossed, susybsm::HSCPCaloInfo::ecalDeDx, ecalDetIdAssociator_, susybsm::HSCPCaloInfo::ecalDetIds, susybsm::HSCPCaloInfo::ecalEnergies, susybsm::HSCPCaloInfo::ecalInvBetaError, susybsm::HSCPCaloInfo::ecalOutOfTimeChi2s, susybsm::HSCPCaloInfo::ecalOutOfTimeEnergies, TrackDetMatchInfo::EcalRecHits, susybsm::HSCPCaloInfo::ecalTime, susybsm::HSCPCaloInfo::ecalTimeError, susybsm::HSCPCaloInfo::ecalTimeErrors, susybsm::HSCPCaloInfo::ecalTimes, susybsm::HSCPCaloInfo::ecalTrackExitPositions, susybsm::HSCPCaloInfo::ecalTrackLengths, EERecHitCollectionToken_, edm::SortedCollection< T, SORT >::end(), EcalRecHit::energy(), edm::SortedCollection< T, SORT >::find(), edm::EventSetup::get(), edm::Event::getByToken(), getDetailedTrackLengthInXtals(), TrackDetectorAssociator::getFreeTrajectoryState(), TrackDetMatchInfo::getPosition(), susybsm::HSCParticle::hasTrackRef(), susybsm::HSCPCaloInfo::hcal3by3dir, susybsm::HSCPCaloInfo::hcal5by5dir, susybsm::HSCPCaloInfo::hcalCrossedEnergy, TrackDetMatchInfo::HcalRecHits, susybsm::HSCPCaloInfo::hoCrossedEnergy, TrackDetMatchInfo::HORecHits, info(), trajectoryStateTransform::innerFreeState(), EcalRecHit::isTimeErrorValid(), EcalRecHit::isTimeValid(), EcalRecHit::kGood, EcalRecHit::kOutOfTime, EcalRecHit::kPoorCalib, TrackDetMatchInfo::nXnEnergy(), reco::Track::outerEta(), reco::Track::outerPhi(), parameters_, position, funct::pow(), edm::ESHandle< T >::product(), EcalRecHit::recoFlag(), mps_fire::result, fastsim::Constants::speedOfLight, mathSSE::sqrt(), DetId::subdetId(), theCaloGeometry_, EcalRecHit::time(), EcalRecHit::timeError(), HiIsolationCommonParameters_cff::track, trackAssociator_, susybsm::HSCParticle::trackRef(), susybsm::HSCPCaloInfo::trkIsoDr, x, PV3DBase< T, PVType, FrameType >::x(), y, PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                                                                                                                  {
   bool setCalo = false;
   HSCPCaloInfo result;

   // EcalDetIdAssociator
   iSetup.get<DetIdAssociatorRecord>().get("EcalDetIdAssociator", ecalDetIdAssociator_);
   // Get the Bfield
   iSetup.get<IdealMagneticFieldRecord>().get(bField_);
   // Geometry
   iSetup.get<CaloGeometryRecord>().get(theCaloGeometry_);
   const CaloGeometry* theGeometry = theCaloGeometry_.product();
   // Topology
   edm::ESHandle<CaloTopology> pCaloTopology;
   iSetup.get<CaloTopologyRecord>().get(pCaloTopology);
   const CaloTopology* theCaloTopology = pCaloTopology.product();
   // EcalRecHits
   edm::Handle<EBRecHitCollection> ebRecHits;
   iEvent.getByToken(EBRecHitCollectionToken_,ebRecHits);
   edm::Handle<EERecHitCollection> eeRecHits;
   iEvent.getByToken(EERecHitCollectionToken_,eeRecHits);

   // select the track
   reco::Track track;
   if(candidate.hasTrackRef())
     track = *(candidate.trackRef());
   else
     return; // in case there is no track ref, we can't do much

   // compute the track isolation
   result.trkIsoDr=100;
   for(reco::TrackCollection::const_iterator ndTrack = tracks->begin(); ndTrack != tracks->end(); ++ndTrack) {
       double dr=sqrt(pow((track.outerEta()-ndTrack->outerEta()),2)+pow((track.outerPhi()-ndTrack->outerPhi()),2));
       if(dr>0.00001 && dr<result.trkIsoDr) result.trkIsoDr=dr;
   }

   // use the track associator to propagate to the calo
   TrackDetMatchInfo info = trackAssociator_.associate( iEvent, iSetup,
                                                        trackAssociator_.getFreeTrajectoryState(iSetup, track),
                                                        parameters_ );

   // do a custom propagation through Ecal
   std::map<int,GlobalPoint> trackExitPositionMap; // rawId to exit position (subtracting cry center)
   std::map<int,float> trackCrossedXtalCurvedMap; // rawId to trackLength

   FreeTrajectoryState tkInnerState = trajectoryStateTransform::innerFreeState(track, &*bField_);
   // Build set of points in Ecal (necklace) using the propagator
   std::vector<SteppingHelixStateInfo> neckLace;
   neckLace = calcEcalDeposit(&tkInnerState,*ecalDetIdAssociator_);
   // Initialize variables to be filled by the track-length function
   double totalLengthCurved = 0.;
   GlobalPoint internalPointCurved(0., 0., 0.);
   GlobalPoint externalPointCurved(0., 0., 0.);
   if(neckLace.size() > 1)
   {
     getDetailedTrackLengthInXtals(trackExitPositionMap,
         trackCrossedXtalCurvedMap,
         totalLengthCurved,
         internalPointCurved,
         externalPointCurved,
         & (*theGeometry),
         & (*theCaloTopology),
         neckLace);
   }

   // Make weighted sum of times
   float sumWeightedTime = 0;
   float sumTimeErrorSqr = 0;
   float sumEnergy = 0;
   float sumTrackLength = 0;
   std::vector<EcalRecHit> crossedRecHits;
   EcalRecHitCollection::const_iterator thisHit;

   std::map<int,GlobalPoint>::const_iterator trackExitMapIt = trackExitPositionMap.begin();
   for(std::map<int,float>::const_iterator mapIt = trackCrossedXtalCurvedMap.begin();
       mapIt != trackCrossedXtalCurvedMap.end(); ++mapIt)
   {
     if(DetId(mapIt->first).subdetId()==EcalBarrel)
     {
       EBDetId ebDetId(mapIt->first);
       thisHit = ebRecHits->find(ebDetId);
       if(thisHit == ebRecHits->end())
       {
         //std::cout << "\t Could not find crossedEcal detId: " << ebDetId << " in EBRecHitCollection!" << std::endl;
         continue;
       }
       const EcalRecHit hit = *thisHit;
       // Cut out badly-reconstructed hits
       if(!hit.isTimeValid())
         continue;
       uint32_t rhFlag = hit.recoFlag();
       if((rhFlag != EcalRecHit::kGood) && (rhFlag != EcalRecHit::kOutOfTime) && (rhFlag != EcalRecHit::kPoorCalib))
         continue;

       float errorOnThis = hit.timeError();
       sumTrackLength+=mapIt->second;
       sumEnergy+=hit.energy();
       crossedRecHits.push_back(hit);
//       result.ecalSwissCrossKs.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kSwissCross));
//       result.ecalE1OverE9s.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kE1OverE9));
       result.ecalTrackLengths.push_back(mapIt->second);
       result.ecalTrackExitPositions.push_back(trackExitMapIt->second);
       result.ecalEnergies.push_back(hit.energy());
       result.ecalTimes.push_back(hit.time());
       result.ecalTimeErrors.push_back(hit.timeError());
       result.ecalOutOfTimeEnergies.push_back(0.);
       result.ecalOutOfTimeChi2s.push_back(0.);
       result.ecalChi2s.push_back(hit.chi2());
       result.ecalDetIds.push_back(ebDetId);
       // SIC DEBUG
       //std::cout << " SIC DEBUG: time error on this crossed RecHit: " << errorOnThis << " energy of hit: "
       //  << hit.energy() << " time of hit: " << hit.time() << " trackLength: " << mapIt->second << std::endl;

       if(hit.isTimeErrorValid()) // use hit time for weighted time average
       {
         sumWeightedTime+=hit.time()/(errorOnThis*errorOnThis);
         sumTimeErrorSqr+=1/(errorOnThis*errorOnThis);
       }
     }
     trackExitMapIt++;
   }

   if(!crossedRecHits.empty())
   {
     setCalo = true;
     sort(crossedRecHits.begin(),crossedRecHits.end(),
             [](auto& x, auto& y){return (x.energy() > y.energy());});
     result.ecalCrossedEnergy = sumEnergy;
     result.ecalCrysCrossed = crossedRecHits.size();
     result.ecalDeDx = sumEnergy/sumTrackLength;
     // replace the below w/o trackassociator quantities?
     result.ecal3by3dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
     result.ecal5by5dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);

     if(sumTimeErrorSqr > 0)
     {
       result.ecalTime = sumWeightedTime/sumTimeErrorSqr;
       result.ecalTimeError = sqrt(1/sumTimeErrorSqr);
       DetId maxEnergyId = crossedRecHits.begin()->id();

       if(maxEnergyId != DetId()) // double check
       {
         // To get beta, we assume photon propagation time is about the same for muons and e/gamma
         // Since the typical path length is >> crystal length, this shouldn't be too bad
         GlobalPoint position = info.getPosition(maxEnergyId); // position of crystal center on front face
         double frontFaceR = sqrt(pow(position.x(),2)+pow(position.y(),2)+pow(position.z(),2));
         double muonShowerMax = frontFaceR+11.5; // assume muon "showerMax" is halfway into the crystal
         double gammaShowerMax = frontFaceR+6.23; // 7 X0 for e/gamma showerMax
         double speedOfLight = 29.979; // cm/ns
         result.ecalBeta = (muonShowerMax)/(result.ecalTime*speedOfLight+gammaShowerMax);
         result.ecalBetaError = (speedOfLight*muonShowerMax*result.ecalTimeError)/pow(speedOfLight*result.ecalTime+gammaShowerMax,2);
         result.ecalInvBetaError = speedOfLight*result.ecalTimeError/muonShowerMax;
       }
       // SIC debug
       //std::cout << "BetaCalcEcal: CrossedRecHits: " << crossedRecHits.size()
       //  << " ecalTime: " << result.ecalTime << " timeError: " << result.ecalTimeError
       //  << " ecalCrossedEnergy: " << result.ecalCrossedEnergy << " ecalBeta: " << result.ecalBeta
       //  << " ecalBetaError: " << result.ecalBetaError <<  " ecalDeDx (MeV/cm): " << 1000*result.ecalDeDx << std::endl;
     }
   }

   if(!info.crossedHcalRecHits.empty())
   {
     // HCAL (not ECAL) info
     result.hcalCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
     result.hoCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
     //maxEnergyId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
     result.hcal3by3dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
     result.hcal5by5dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2);
   }

   if(setCalo)
     caloInfo = result;
}

void BetaCalculatorECAL::addStepToXtal ( std::map< int, GlobalPoint > &  trackExitPositionMap, std::map< int, float > &  trackCrossedXtalMap, DetId  aDetId, float  step, GlobalPoint  point, const CaloSubdetectorGeometry *  theSubdetGeometry  )

private

Definition at line 318 of file BetaCalculatorECAL.cc.

References diffTreeTool::diff, CaloSubdetectorGeometry::getGeometry(), AlCaHLTBitMon_ParallelJobs::p, DetId::rawId(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by getDetailedTrackLengthInXtals().

{

  auto   cell_p = theSubdetGeometry->getGeometry(aDetId);
  GlobalPoint p = cell_p->getPosition(23);
  GlobalPoint diff(point.x()-p.x(),point.y()-p.y(),point.z()-p.z());

  std::map<int,GlobalPoint>::iterator xtal = trackExitPositionMap.find(aDetId.rawId());
  if (xtal!=trackExitPositionMap.end())
    ((*xtal).second)=diff;
  else
    trackExitPositionMap.insert(std::pair<int,GlobalPoint>(aDetId.rawId(),diff));

  std::map<int,float>::iterator xtal2 = trackCrossedXtalMap.find(aDetId.rawId());
  if (xtal2!= trackCrossedXtalMap.end())
    ((*xtal2).second)+=step;
  else
    trackCrossedXtalMap.insert(std::pair<int,float>(aDetId.rawId(),step));
}

std::vector< SteppingHelixStateInfo > BetaCalculatorECAL::calcEcalDeposit ( const FreeTrajectoryState *  tkInnerState, const DetIdAssociator &  associator  )

private

Definition at line 203 of file BetaCalculatorECAL.cc.

References alongMomentum, bField_, FiducialVolume::maxR(), pvSelector_cfi::maxZ, FiducialVolume::maxZ(), FiducialVolume::minR(), HiClusterCompatibility_cfi::minZ, FiducialVolume::minZ(), propagateThoughFromIP(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, and DetIdAssociator::volume().

Referenced by addInfoToCandidate().

{
   // Set some parameters
   double minR = associator.volume().minR () ;
   double minZ = associator.volume().minZ () ;
   double maxR = associator.volume().maxR () ;
   double maxZ = associator.volume().maxZ () ;

   // Define the TrackOrigin (where the propagation starts)
   SteppingHelixStateInfo trackOrigin(*tkInnerState);

   // Define Propagator
   SteppingHelixPropagator* prop = new SteppingHelixPropagator (&*bField_, alongMomentum);
   prop -> setMaterialMode(false);
   prop -> applyRadX0Correction(true);

   return propagateThoughFromIP(trackOrigin,prop,associator.volume(), 500,0.1,minR,minZ,maxR,maxZ);
}

int BetaCalculatorECAL::getDetailedTrackLengthInXtals ( std::map< int, GlobalPoint > &  trackExitPositionMap, std::map< int, float > &  trackCrossedXtalMap, double &  totalLengthCurved, GlobalPoint &  internalPointCurved, GlobalPoint &  externalPointCurved, const CaloGeometry *  theGeometry, const CaloTopology *  theTopology, const std::vector< SteppingHelixStateInfo > &  neckLace  )

private

Definition at line 223 of file BetaCalculatorECAL.cc.

References addStepToXtal(), DetId::Ecal, ecaldqm::getGeometry(), CaloGeometry::getSubdetectorGeometry(), CaloTopology::getSubdetectorTopology(), CaloSubdetectorTopology::getWindow(), gen::k, mag(), point, and funct::true.

Referenced by addInfoToCandidate().

{
   GlobalPoint origin (0., 0., 0.);
   internalPointCurved = origin ;
   externalPointCurved = origin ;

   bool firstPoint = false;
   trackCrossedXtalMap.clear();

   const CaloSubdetectorGeometry* theBarrelSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,1);
   const CaloSubdetectorGeometry* theEndcapSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,2);

   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = (neckLace.begin() + 1); itr != neckLace.end(); ++itr)
   {
     GlobalPoint probe_gp = (*itr).position();
     std::vector<DetId> surroundingMatrix;

     EBDetId closestBarrelDetIdToProbe = ((theBarrelSubdetGeometry -> getClosestCell(probe_gp)).rawId());
     EEDetId closestEndcapDetIdToProbe = ((theEndcapSubdetGeometry -> getClosestCell(probe_gp)).rawId());

     // check if the probe is inside the xtal
     if( (closestEndcapDetIdToProbe) && 
     (theGeometry->getSubdetectorGeometry(closestEndcapDetIdToProbe)->
      getGeometry(closestEndcapDetIdToProbe)->inside(probe_gp)) ) 
     {
       double step = ((*itr).position() - (*(itr-1)).position()).mag();
       GlobalPoint point = itr->position();
       addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestEndcapDetIdToProbe, step, point, theEndcapSubdetGeometry);
       totalLengthCurved += step;

       if (firstPoint == false)
       {
         internalPointCurved = probe_gp ;
         firstPoint = true ;
       }

       externalPointCurved = probe_gp ;
     }

     if( (closestBarrelDetIdToProbe) && 
     (theGeometry->getSubdetectorGeometry(closestBarrelDetIdToProbe)->
      getGeometry(closestBarrelDetIdToProbe)->inside(probe_gp)) )
     {
       double step = ((*itr).position() - (*(itr-1)).position()).mag();
       GlobalPoint point = itr->position();
       addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestBarrelDetIdToProbe, step, point, theBarrelSubdetGeometry);
       totalLengthCurved += step;

       if (firstPoint == false)
       {
         internalPointCurved = probe_gp ;
         firstPoint = true ;
       }

       externalPointCurved = probe_gp ;
     }
     else
     {
       // 3x3 matrix surrounding the probe
       surroundingMatrix = theTopology->getSubdetectorTopology(closestBarrelDetIdToProbe)->getWindow(closestBarrelDetIdToProbe,3,3);

       for( unsigned int k=0; k<surroundingMatrix.size(); ++k ) {
         if(theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k))->getGeometry(surroundingMatrix.at(k))->inside(probe_gp))
         {
           double step = ((*itr).position() - (*(itr-1)).position()).mag();
           GlobalPoint point = itr->position();
           addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, surroundingMatrix[k], step,
               point, theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k)));
           totalLengthCurved += step;

           if (firstPoint == false)
           {
             internalPointCurved = probe_gp ;
             firstPoint = true ;
           }

           externalPointCurved = probe_gp ;
         }
       }

       // clear neighborhood matrix
       surroundingMatrix.clear();
     }
   }

   return 0;
}

Member Data Documentation

edm::ESHandle<MagneticField> BetaCalculatorECAL::bField_

private

Definition at line 64 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate(), and calcEcalDeposit().

edm::EDGetTokenT<EBRecHitCollection> BetaCalculatorECAL::EBRecHitCollectionToken_

private

Definition at line 60 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate().

edm::ESHandle<DetIdAssociator> BetaCalculatorECAL::ecalDetIdAssociator_

private

Definition at line 63 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate().

edm::EDGetTokenT<EERecHitCollection> BetaCalculatorECAL::EERecHitCollectionToken_

private

Definition at line 61 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate().

TrackAssociatorParameters BetaCalculatorECAL::parameters_

private

Definition at line 59 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate(), BetaCalculatorECAL(), Mixins._TypedParameterizable::clone(), Types.PSet::clone(), Types.EDAlias::clone(), and Mixins._TypedParameterizable::copy().

edm::ESHandle<CaloGeometry> BetaCalculatorECAL::theCaloGeometry_

private

Definition at line 65 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate().

TrackDetectorAssociator BetaCalculatorECAL::trackAssociator_

private

Definition at line 58 of file BetaCalculatorECAL.h.

Referenced by addInfoToCandidate(), and BetaCalculatorECAL().