#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)
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::InputTag	EBRecHitCollection_
edm::ESHandle< DetIdAssociator >	ecalDetIdAssociator_
edm::InputTag	EERecHitCollection_
TrackAssociatorParameters	parameters_
edm::ESHandle< CaloGeometry >	theCaloGeometry_
TrackDetectorAssociator	trackAssociator_

Detailed Description

Definition at line 31 of file BetaCalculatorECAL.h.

Constructor & Destructor Documentation

BetaCalculatorECAL::BetaCalculatorECAL ( const edm::ParameterSet & iConfig )

Definition at line 19 of file BetaCalculatorECAL.cc.

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

                                                                     :
  EBRecHitCollection_ (iConfig.getParameter<edm::InputTag>("EBRecHitCollection")),
  EERecHitCollection_ (iConfig.getParameter<edm::InputTag>("EERecHitCollection"))
{
   edm::ParameterSet trkParameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
   parameters_.loadParameters( trkParameters ); 
   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 30 of file BetaCalculatorECAL.cc.

                                                                                                                                                                                  {
   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.getByLabel(EBRecHitCollection_,ebRecHits);
   edm::Handle<EERecHitCollection> eeRecHits;
   iEvent.getByLabel(EERecHitCollection_,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(hit.outOfTimeEnergy());
       result.ecalOutOfTimeChi2s.push_back(hit.outOfTimeChi2());
       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.size() > 0)
   {
     setCalo = true;
     sort(crossedRecHits.begin(),crossedRecHits.end(),EcalRecHitLess());
     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.size() > 0)
   {
     // 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 316 of file BetaCalculatorECAL.cc.

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

Referenced by getDetailedTrackLengthInXtals().

{

  const CaloCellGeometry *cell_p = theSubdetGeometry->getGeometry(aDetId);
  GlobalPoint p = (dynamic_cast <const TruncatedPyramid *> (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(), CosmicsPD_Skims::maxZ, FiducialVolume::maxZ(), FiducialVolume::minR(), cmsRun_displayProdMFGeom_cfg::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, CaloGeometry::getSubdetectorGeometry(), CaloTopology::getSubdetectorTopology(), CaloSubdetectorTopology::getWindow(), gen::k, mag(), point, launcher::step, 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;
}