Reference Manual

Go to the documentation of this file.
00001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00002 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00003 
00004 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00005 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
00006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00007 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
00008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00009 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixStateInfo.h"
00010 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00011 
00012 #include "RecoEcal/EgammaCoreTools/interface/EcalRecHitLess.h"
00013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00014 
00015 #include "SUSYBSMAnalysis/HSCP/interface/BetaCalculatorECAL.h"
00016 
00017 using namespace susybsm;
00018 
00019 BetaCalculatorECAL::BetaCalculatorECAL(const edm::ParameterSet& iConfig) :
00020   EBRecHitCollection_ (iConfig.getParameter<edm::InputTag>("EBRecHitCollection")),
00021   EERecHitCollection_ (iConfig.getParameter<edm::InputTag>("EERecHitCollection"))
00022 {
00023    edm::ParameterSet trkParameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00024    parameters_.loadParameters( trkParameters ); 
00025    trackAssociator_.useDefaultPropagator();
00026 
00027 }
00028 
00029 
00030 void BetaCalculatorECAL::addInfoToCandidate(HSCParticle& candidate, edm::Handle<reco::TrackCollection>& tracks, edm::Event& iEvent, const edm::EventSetup& iSetup,  HSCPCaloInfo& caloInfo) {
00031    bool setCalo = false;
00032    HSCPCaloInfo result;
00033 
00034    // EcalDetIdAssociator
00035    iSetup.get<DetIdAssociatorRecord>().get("EcalDetIdAssociator", ecalDetIdAssociator_);
00036    // Get the Bfield
00037    iSetup.get<IdealMagneticFieldRecord>().get(bField_);
00038    // Geometry
00039    iSetup.get<CaloGeometryRecord>().get(theCaloGeometry_);
00040    const CaloGeometry* theGeometry = theCaloGeometry_.product();
00041    // Topology
00042    edm::ESHandle<CaloTopology> pCaloTopology;
00043    iSetup.get<CaloTopologyRecord>().get(pCaloTopology);
00044    const CaloTopology* theCaloTopology = pCaloTopology.product();
00045    // EcalRecHits
00046    edm::Handle<EBRecHitCollection> ebRecHits;
00047    iEvent.getByLabel(EBRecHitCollection_,ebRecHits);
00048    edm::Handle<EERecHitCollection> eeRecHits;
00049    iEvent.getByLabel(EERecHitCollection_,eeRecHits);
00050    
00051    // select the track
00052    reco::Track track;
00053    if(candidate.hasTrackRef())
00054      track = *(candidate.trackRef());
00055    else
00056      return; // in case there is no track ref, we can't do much
00057 
00058    // compute the track isolation
00059    result.trkIsoDr=100;
00060    for(reco::TrackCollection::const_iterator ndTrack = tracks->begin(); ndTrack != tracks->end(); ++ndTrack) {
00061        double dr=sqrt(pow((track.outerEta()-ndTrack->outerEta()),2)+pow((track.outerPhi()-ndTrack->outerPhi()),2));
00062        if(dr>0.00001 && dr<result.trkIsoDr) result.trkIsoDr=dr;
00063    }
00064 
00065    // use the track associator to propagate to the calo
00066    TrackDetMatchInfo info = trackAssociator_.associate( iEvent, iSetup, 
00067                                                         trackAssociator_.getFreeTrajectoryState(iSetup, track),
00068                                                         parameters_ );
00069 
00070    // do a custom propagation through Ecal
00071    std::map<int,GlobalPoint> trackExitPositionMap; // rawId to exit position (subtracting cry center)
00072    std::map<int,float> trackCrossedXtalCurvedMap; // rawId to trackLength
00073    TrajectoryStateTransform tsTransform;
00074    FreeTrajectoryState tkInnerState = tsTransform.innerFreeState(track, &*bField_);
00075    // Build set of points in Ecal (necklace) using the propagator
00076    std::vector<SteppingHelixStateInfo> neckLace;
00077    neckLace = calcEcalDeposit(&tkInnerState,*ecalDetIdAssociator_);
00078    // Initialize variables to be filled by the track-length function
00079    double totalLengthCurved = 0.;
00080    GlobalPoint internalPointCurved(0., 0., 0.);
00081    GlobalPoint externalPointCurved(0., 0., 0.);
00082    if(neckLace.size() > 1)
00083    {
00084      getDetailedTrackLengthInXtals(trackExitPositionMap,
00085          trackCrossedXtalCurvedMap,
00086          totalLengthCurved,
00087          internalPointCurved,
00088          externalPointCurved,
00089          & (*theGeometry),
00090          & (*theCaloTopology),
00091          neckLace);
00092    }
00093    int numCrysCrossed = -1;
00094    numCrysCrossed = trackCrossedXtalCurvedMap.size();
00095 
00096    // Make weighted sum of times
00097    float sumWeightedTime = 0;
00098    float sumTimeErrorSqr = 0;
00099    float sumEnergy = 0;
00100    float sumTrackLength = 0;
00101    std::vector<EcalRecHit> crossedRecHits;
00102    EcalRecHitCollection::const_iterator thisHit;
00103 
00104    std::map<int,GlobalPoint>::const_iterator trackExitMapIt = trackExitPositionMap.begin();
00105    for(std::map<int,float>::const_iterator mapIt = trackCrossedXtalCurvedMap.begin();
00106        mapIt != trackCrossedXtalCurvedMap.end(); ++mapIt)
00107    {
00108      if(DetId(mapIt->first).subdetId()==EcalBarrel)
00109      {
00110        EBDetId ebDetId(mapIt->first);
00111        thisHit = ebRecHits->find(ebDetId);
00112        if(thisHit == ebRecHits->end())
00113        {
00114          //std::cout << "\t Could not find crossedEcal detId: " << ebDetId << " in EBRecHitCollection!" << std::endl;
00115          continue;
00116        }
00117        const EcalRecHit hit = *thisHit;
00118        // Cut out badly-reconstructed hits
00119        if(!hit.isTimeValid())
00120          continue;
00121        uint32_t rhFlag = hit.recoFlag();
00122        if((rhFlag != EcalRecHit::kGood) && (rhFlag != EcalRecHit::kOutOfTime) && (rhFlag != EcalRecHit::kPoorCalib))
00123          continue;
00124 
00125        float errorOnThis = hit.timeError();
00126        sumTrackLength+=mapIt->second;
00127        sumEnergy+=hit.energy();
00128        crossedRecHits.push_back(hit);
00129 //       result.ecalSwissCrossKs.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kSwissCross));
00130 //       result.ecalE1OverE9s.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kE1OverE9));
00131        result.ecalTrackLengths.push_back(mapIt->second);
00132        result.ecalTrackExitPositions.push_back(trackExitMapIt->second);
00133        result.ecalEnergies.push_back(hit.energy());
00134        result.ecalTimes.push_back(hit.time());
00135        result.ecalTimeErrors.push_back(hit.timeError());
00136        result.ecalOutOfTimeEnergies.push_back(hit.outOfTimeEnergy());
00137        result.ecalOutOfTimeChi2s.push_back(hit.outOfTimeChi2());
00138        result.ecalChi2s.push_back(hit.chi2());
00139        result.ecalDetIds.push_back(ebDetId);
00140        // SIC DEBUG
00141        //std::cout << " SIC DEBUG: time error on this crossed RecHit: " << errorOnThis << " energy of hit: "
00142        //  << hit.energy() << " time of hit: " << hit.time() << " trackLength: " << mapIt->second << std::endl;
00143 
00144        if(hit.isTimeErrorValid()) // use hit time for weighted time average
00145        {
00146          sumWeightedTime+=hit.time()/(errorOnThis*errorOnThis);
00147          sumTimeErrorSqr+=1/(errorOnThis*errorOnThis);
00148        }
00149      }
00150      trackExitMapIt++;
00151    }
00152 
00153    if(crossedRecHits.size() > 0)
00154    {
00155      setCalo = true;
00156      sort(crossedRecHits.begin(),crossedRecHits.end(),EcalRecHitLess());
00157      result.ecalCrossedEnergy = sumEnergy;
00158      result.ecalCrysCrossed = crossedRecHits.size();
00159      result.ecalDeDx = sumEnergy/sumTrackLength;
00160      // replace the below w/o trackassociator quantities?
00161      result.ecal3by3dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
00162      result.ecal5by5dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
00163 
00164      if(sumTimeErrorSqr > 0)
00165      {
00166        result.ecalTime = sumWeightedTime/sumTimeErrorSqr;
00167        result.ecalTimeError = sqrt(1/sumTimeErrorSqr);
00168        DetId maxEnergyId = crossedRecHits.begin()->id();
00169 
00170        if(maxEnergyId != DetId()) // double check
00171        {
00172          // To get beta, we assume photon propagation time is about the same for muons and e/gamma
00173          // Since the typical path length is >> crystal length, this shouldn't be too bad
00174          GlobalPoint position = info.getPosition(maxEnergyId); // position of crystal center on front face
00175          double frontFaceR = sqrt(pow(position.x(),2)+pow(position.y(),2)+pow(position.z(),2));
00176          double muonShowerMax = frontFaceR+11.5; // assume muon "showerMax" is halfway into the crystal
00177          double gammaShowerMax = frontFaceR+6.23; // 7 X0 for e/gamma showerMax
00178          double speedOfLight = 29.979; // cm/ns
00179          result.ecalBeta = (muonShowerMax)/(result.ecalTime*speedOfLight+gammaShowerMax);
00180          result.ecalBetaError = (speedOfLight*muonShowerMax*result.ecalTimeError)/pow(speedOfLight*result.ecalTime+gammaShowerMax,2);
00181          result.ecalInvBetaError = speedOfLight*result.ecalTimeError/muonShowerMax;
00182        }
00183        // SIC debug
00184        //std::cout << "BetaCalcEcal: CrossedRecHits: " << crossedRecHits.size()
00185        //  << " ecalTime: " << result.ecalTime << " timeError: " << result.ecalTimeError
00186        //  << " ecalCrossedEnergy: " << result.ecalCrossedEnergy << " ecalBeta: " << result.ecalBeta 
00187        //  << " ecalBetaError: " << result.ecalBetaError <<  " ecalDeDx (MeV/cm): " << 1000*result.ecalDeDx << std::endl;
00188      }
00189    }
00190 
00191    if(info.crossedHcalRecHits.size() > 0)
00192    {
00193      // HCAL (not ECAL) info
00194      result.hcalCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
00195      result.hoCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
00196      //maxEnergyId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
00197      result.hcal3by3dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
00198      result.hcal5by5dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2);
00199    }
00200 
00201    if(setCalo) 
00202      caloInfo = result;
00203 }
00204 
00205 std::vector<SteppingHelixStateInfo> BetaCalculatorECAL::calcEcalDeposit(const FreeTrajectoryState* tkInnerState,
00206             const DetIdAssociator& associator)
00207 {   
00208    // Set some parameters
00209    double minR = associator.volume().minR () ;
00210    double minZ = associator.volume().minZ () ;
00211    double maxR = associator.volume().maxR () ;
00212    double maxZ = associator.volume().maxZ () ;
00213 
00214    // Define the TrackOrigin (where the propagation starts)
00215    SteppingHelixStateInfo trackOrigin(*tkInnerState);
00216 
00217    // Define Propagator
00218    SteppingHelixPropagator* prop = new SteppingHelixPropagator (&*bField_, alongMomentum);
00219    prop -> setMaterialMode(false); 
00220    prop -> applyRadX0Correction(true);
00221 
00222    return propagateThoughFromIP(trackOrigin,prop,associator.volume(), 500,0.1,minR,minZ,maxR,maxZ);
00223 }
00224 
00225 int BetaCalculatorECAL::getDetailedTrackLengthInXtals(std::map<int,GlobalPoint>& trackExitPositionMap,
00226    std::map<int,float>& trackCrossedXtalMap,
00227    double& totalLengthCurved,
00228    GlobalPoint& internalPointCurved,
00229    GlobalPoint& externalPointCurved,
00230    const CaloGeometry* theGeometry,
00231    const CaloTopology* theTopology,
00232    const std::vector<SteppingHelixStateInfo>& neckLace)
00233 {
00234    GlobalPoint origin (0., 0., 0.);
00235    internalPointCurved = origin ;
00236    externalPointCurved = origin ;
00237 
00238    bool firstPoint = false;
00239    trackCrossedXtalMap.clear();
00240 
00241    const CaloSubdetectorGeometry* theBarrelSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,1);
00242    const CaloSubdetectorGeometry* theEndcapSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,2);
00243 
00244    for(std::vector<SteppingHelixStateInfo>::const_iterator itr = (neckLace.begin() + 1); itr != neckLace.end(); ++itr)
00245    {
00246      GlobalPoint probe_gp = (*itr).position();
00247      GlobalVector probe_gv = (*itr).momentum();
00248      std::vector<DetId> surroundingMatrix;
00249 
00250      EBDetId closestBarrelDetIdToProbe = ((theBarrelSubdetGeometry -> getClosestCell(probe_gp)).rawId());
00251      EEDetId closestEndcapDetIdToProbe = ((theEndcapSubdetGeometry -> getClosestCell(probe_gp)).rawId());
00252 
00253      // check if the probe is inside the xtal
00254      if( (closestEndcapDetIdToProbe) && (theGeometry->getSubdetectorGeometry(closestEndcapDetIdToProbe)->
00255            getGeometry(closestEndcapDetIdToProbe)->inside(probe_gp)) )
00256      {
00257        double step = ((*itr).position() - (*(itr-1)).position()).mag();
00258        GlobalPoint point = itr->position();
00259        addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestEndcapDetIdToProbe, step, point, theEndcapSubdetGeometry);
00260        totalLengthCurved += step;
00261 
00262        if (firstPoint == false)
00263        {
00264          internalPointCurved = probe_gp ;
00265          firstPoint = true ;
00266        }
00267 
00268        externalPointCurved = probe_gp ;
00269      }
00270 
00271      if( (closestBarrelDetIdToProbe) && (theGeometry->getSubdetectorGeometry(closestBarrelDetIdToProbe)->
00272            getGeometry(closestBarrelDetIdToProbe)->inside(probe_gp)) )
00273      {
00274        double step = ((*itr).position() - (*(itr-1)).position()).mag();
00275        GlobalPoint point = itr->position();
00276        addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestBarrelDetIdToProbe, step, point, theBarrelSubdetGeometry);
00277        totalLengthCurved += step;
00278 
00279        if (firstPoint == false)
00280        {
00281          internalPointCurved = probe_gp ;
00282          firstPoint = true ;
00283        }
00284 
00285        externalPointCurved = probe_gp ;
00286      }
00287      else
00288      {
00289        // 3x3 matrix surrounding the probe
00290        surroundingMatrix = theTopology->getSubdetectorTopology(closestBarrelDetIdToProbe)->getWindow(closestBarrelDetIdToProbe,3,3);
00291 
00292        for( unsigned int k=0; k<surroundingMatrix.size(); ++k ) {
00293          if(theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k))->getGeometry(surroundingMatrix.at(k))->inside(probe_gp))
00294          {
00295            double step = ((*itr).position() - (*(itr-1)).position()).mag();
00296            GlobalPoint point = itr->position();
00297            addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, surroundingMatrix[k], step,
00298                point, theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k)));
00299            totalLengthCurved += step;
00300 
00301            if (firstPoint == false)
00302            {
00303              internalPointCurved = probe_gp ;
00304              firstPoint = true ;
00305            }
00306 
00307            externalPointCurved = probe_gp ;
00308          }
00309        }
00310 
00311        // clear neighborhood matrix
00312        surroundingMatrix.clear();
00313      }
00314    }
00315 
00316    return 0;
00317 }
00318 
00319 void BetaCalculatorECAL::addStepToXtal(std::map<int,GlobalPoint>& trackExitPositionMap,
00320     std::map<int,float>& trackCrossedXtalMap,
00321     DetId aDetId,
00322     float step,
00323     GlobalPoint point,
00324     const CaloSubdetectorGeometry* theSubdetGeometry)
00325 {
00326 
00327   const CaloCellGeometry *cell_p = theSubdetGeometry->getGeometry(aDetId);
00328   GlobalPoint p = (dynamic_cast <const TruncatedPyramid *> (cell_p))->getPosition(23);
00329   GlobalPoint diff(point.x()-p.x(),point.y()-p.y(),point.z()-p.z());
00330 
00331   std::map<int,GlobalPoint>::iterator xtal = trackExitPositionMap.find(aDetId.rawId());
00332   if (xtal!=trackExitPositionMap.end())
00333     ((*xtal).second)=diff;
00334   else
00335     trackExitPositionMap.insert(std::pair<int,GlobalPoint>(aDetId.rawId(),diff));
00336 
00337   std::map<int,float>::iterator xtal2 = trackCrossedXtalMap.find(aDetId.rawId());
00338   if (xtal2!= trackCrossedXtalMap.end())
00339     ((*xtal2).second)+=step;
00340   else
00341     trackCrossedXtalMap.insert(std::pair<int,float>(aDetId.rawId(),step));
00342 }
00343 
00344 
00345
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/SUSYBSMAnalysis/HSCP/src/BetaCalculatorECAL.cc
