#include <BetaCalculatorTK.h>

Public Member Functions
void	addInfoToCandidate (HSCParticle &candidate, edm::Event &iEvent, const edm::EventSetup &iSetup)
	BetaCalculatorTK (const edm::ParameterSet &iConfig)
Public Attributes
edm::InputTag	m_dedxDiscriminator1Tag
edm::InputTag	m_dedxDiscriminator2Tag
edm::InputTag	m_dedxDiscriminator3Tag
edm::InputTag	m_dedxDiscriminator4Tag
edm::InputTag	m_dedxDiscriminator5Tag
edm::InputTag	m_dedxDiscriminator6Tag
edm::InputTag	m_dedxEstimator1Tag
edm::InputTag	m_dedxEstimator2Tag
edm::InputTag	m_dedxEstimator3Tag
edm::InputTag	m_dedxEstimator4Tag
edm::InputTag	m_dedxEstimator5Tag
edm::InputTag	m_dedxEstimator6Tag

Detailed Description

Definition at line 27 of file BetaCalculatorTK.h.

Constructor & Destructor Documentation

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

Definition at line 3 of file BetaCalculatorTK.cc.

                                                                {
/*
  m_dedxEstimator1Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator1");
  m_dedxEstimator2Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator2");
  m_dedxEstimator3Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator3");
  m_dedxEstimator4Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator4");
  m_dedxEstimator5Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator5");
  m_dedxEstimator6Tag     = iConfig.getParameter<edm::InputTag>("dedxEstimator6");
  m_dedxDiscriminator1Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator1");
  m_dedxDiscriminator2Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator2");
  m_dedxDiscriminator3Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator3");
  m_dedxDiscriminator4Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator4");
  m_dedxDiscriminator5Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator5");
  m_dedxDiscriminator6Tag = iConfig.getParameter<edm::InputTag>("dedxDiscriminator6");
*/
}

Member Function Documentation

void BetaCalculatorTK::addInfoToCandidate	(	HSCParticle &	candidate,
		edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

Definition at line 27 of file BetaCalculatorECAL.cc.

Referenced by HSCParticleProducer::filter().

                                                                                                                                                                                  {
   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
   TrajectoryStateTransform tsTransform;
   FreeTrajectoryState tkInnerState = tsTransform.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);
   }
   int numCrysCrossed = -1;
   numCrysCrossed = trackCrossedXtalCurvedMap.size();

   // 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;
}