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    
00074    FreeTrajectoryState tkInnerState = trajectoryStateTransform::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 
00094    // Make weighted sum of times
00095    float sumWeightedTime = 0;
00096    float sumTimeErrorSqr = 0;
00097    float sumEnergy = 0;
00098    float sumTrackLength = 0;
00099    std::vector<EcalRecHit> crossedRecHits;
00100    EcalRecHitCollection::const_iterator thisHit;
00101 
00102    std::map<int,GlobalPoint>::const_iterator trackExitMapIt = trackExitPositionMap.begin();
00103    for(std::map<int,float>::const_iterator mapIt = trackCrossedXtalCurvedMap.begin();
00104        mapIt != trackCrossedXtalCurvedMap.end(); ++mapIt)
00105    {
00106      if(DetId(mapIt->first).subdetId()==EcalBarrel)
00107      {
00108        EBDetId ebDetId(mapIt->first);
00109        thisHit = ebRecHits->find(ebDetId);
00110        if(thisHit == ebRecHits->end())
00111        {
00112          //std::cout << "\t Could not find crossedEcal detId: " << ebDetId << " in EBRecHitCollection!" << std::endl;
00113          continue;
00114        }
00115        const EcalRecHit hit = *thisHit;
00116        // Cut out badly-reconstructed hits
00117        if(!hit.isTimeValid())
00118          continue;
00119        uint32_t rhFlag = hit.recoFlag();
00120        if((rhFlag != EcalRecHit::kGood) && (rhFlag != EcalRecHit::kOutOfTime) && (rhFlag != EcalRecHit::kPoorCalib))
00121          continue;
00122 
00123        float errorOnThis = hit.timeError();
00124        sumTrackLength+=mapIt->second;
00125        sumEnergy+=hit.energy();
00126        crossedRecHits.push_back(hit);
00127 //       result.ecalSwissCrossKs.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kSwissCross));
00128 //       result.ecalE1OverE9s.push_back(EcalSeverityLevelAlgo::spikeFromNeighbours(ebDetId,(*ebRecHits),0.2,EcalSeverityLevelAlgo::kE1OverE9));
00129        result.ecalTrackLengths.push_back(mapIt->second);
00130        result.ecalTrackExitPositions.push_back(trackExitMapIt->second);
00131        result.ecalEnergies.push_back(hit.energy());
00132        result.ecalTimes.push_back(hit.time());
00133        result.ecalTimeErrors.push_back(hit.timeError());
00134        result.ecalOutOfTimeEnergies.push_back(hit.outOfTimeEnergy());
00135        result.ecalOutOfTimeChi2s.push_back(hit.outOfTimeChi2());
00136        result.ecalChi2s.push_back(hit.chi2());
00137        result.ecalDetIds.push_back(ebDetId);
00138        // SIC DEBUG
00139        //std::cout << " SIC DEBUG: time error on this crossed RecHit: " << errorOnThis << " energy of hit: "
00140        //  << hit.energy() << " time of hit: " << hit.time() << " trackLength: " << mapIt->second << std::endl;
00141 
00142        if(hit.isTimeErrorValid()) // use hit time for weighted time average
00143        {
00144          sumWeightedTime+=hit.time()/(errorOnThis*errorOnThis);
00145          sumTimeErrorSqr+=1/(errorOnThis*errorOnThis);
00146        }
00147      }
00148      trackExitMapIt++;
00149    }
00150 
00151    if(crossedRecHits.size() > 0)
00152    {
00153      setCalo = true;
00154      sort(crossedRecHits.begin(),crossedRecHits.end(),EcalRecHitLess());
00155      result.ecalCrossedEnergy = sumEnergy;
00156      result.ecalCrysCrossed = crossedRecHits.size();
00157      result.ecalDeDx = sumEnergy/sumTrackLength;
00158      // replace the below w/o trackassociator quantities?
00159      result.ecal3by3dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
00160      result.ecal5by5dir = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
00161 
00162      if(sumTimeErrorSqr > 0)
00163      {
00164        result.ecalTime = sumWeightedTime/sumTimeErrorSqr;
00165        result.ecalTimeError = sqrt(1/sumTimeErrorSqr);
00166        DetId maxEnergyId = crossedRecHits.begin()->id();
00167 
00168        if(maxEnergyId != DetId()) // double check
00169        {
00170          // To get beta, we assume photon propagation time is about the same for muons and e/gamma
00171          // Since the typical path length is >> crystal length, this shouldn't be too bad
00172          GlobalPoint position = info.getPosition(maxEnergyId); // position of crystal center on front face
00173          double frontFaceR = sqrt(pow(position.x(),2)+pow(position.y(),2)+pow(position.z(),2));
00174          double muonShowerMax = frontFaceR+11.5; // assume muon "showerMax" is halfway into the crystal
00175          double gammaShowerMax = frontFaceR+6.23; // 7 X0 for e/gamma showerMax
00176          double speedOfLight = 29.979; // cm/ns
00177          result.ecalBeta = (muonShowerMax)/(result.ecalTime*speedOfLight+gammaShowerMax);
00178          result.ecalBetaError = (speedOfLight*muonShowerMax*result.ecalTimeError)/pow(speedOfLight*result.ecalTime+gammaShowerMax,2);
00179          result.ecalInvBetaError = speedOfLight*result.ecalTimeError/muonShowerMax;
00180        }
00181        // SIC debug
00182        //std::cout << "BetaCalcEcal: CrossedRecHits: " << crossedRecHits.size()
00183        //  << " ecalTime: " << result.ecalTime << " timeError: " << result.ecalTimeError
00184        //  << " ecalCrossedEnergy: " << result.ecalCrossedEnergy << " ecalBeta: " << result.ecalBeta 
00185        //  << " ecalBetaError: " << result.ecalBetaError <<  " ecalDeDx (MeV/cm): " << 1000*result.ecalDeDx << std::endl;
00186      }
00187    }
00188 
00189    if(info.crossedHcalRecHits.size() > 0)
00190    {
00191      // HCAL (not ECAL) info
00192      result.hcalCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
00193      result.hoCrossedEnergy = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
00194      //maxEnergyId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
00195      result.hcal3by3dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
00196      result.hcal5by5dir = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2);
00197    }
00198 
00199    if(setCalo) 
00200      caloInfo = result;
00201 }
00202 
00203 std::vector<SteppingHelixStateInfo> BetaCalculatorECAL::calcEcalDeposit(const FreeTrajectoryState* tkInnerState,
00204             const DetIdAssociator& associator)
00205 {   
00206    // Set some parameters
00207    double minR = associator.volume().minR () ;
00208    double minZ = associator.volume().minZ () ;
00209    double maxR = associator.volume().maxR () ;
00210    double maxZ = associator.volume().maxZ () ;
00211 
00212    // Define the TrackOrigin (where the propagation starts)
00213    SteppingHelixStateInfo trackOrigin(*tkInnerState);
00214 
00215    // Define Propagator
00216    SteppingHelixPropagator* prop = new SteppingHelixPropagator (&*bField_, alongMomentum);
00217    prop -> setMaterialMode(false); 
00218    prop -> applyRadX0Correction(true);
00219 
00220    return propagateThoughFromIP(trackOrigin,prop,associator.volume(), 500,0.1,minR,minZ,maxR,maxZ);
00221 }
00222 
00223 int BetaCalculatorECAL::getDetailedTrackLengthInXtals(std::map<int,GlobalPoint>& trackExitPositionMap,
00224    std::map<int,float>& trackCrossedXtalMap,
00225    double& totalLengthCurved,
00226    GlobalPoint& internalPointCurved,
00227    GlobalPoint& externalPointCurved,
00228    const CaloGeometry* theGeometry,
00229    const CaloTopology* theTopology,
00230    const std::vector<SteppingHelixStateInfo>& neckLace)
00231 {
00232    GlobalPoint origin (0., 0., 0.);
00233    internalPointCurved = origin ;
00234    externalPointCurved = origin ;
00235 
00236    bool firstPoint = false;
00237    trackCrossedXtalMap.clear();
00238 
00239    const CaloSubdetectorGeometry* theBarrelSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,1);
00240    const CaloSubdetectorGeometry* theEndcapSubdetGeometry = theGeometry->getSubdetectorGeometry(DetId::Ecal,2);
00241 
00242    for(std::vector<SteppingHelixStateInfo>::const_iterator itr = (neckLace.begin() + 1); itr != neckLace.end(); ++itr)
00243    {
00244      GlobalPoint probe_gp = (*itr).position();
00245      std::vector<DetId> surroundingMatrix;
00246 
00247      EBDetId closestBarrelDetIdToProbe = ((theBarrelSubdetGeometry -> getClosestCell(probe_gp)).rawId());
00248      EEDetId closestEndcapDetIdToProbe = ((theEndcapSubdetGeometry -> getClosestCell(probe_gp)).rawId());
00249 
00250      // check if the probe is inside the xtal
00251      if( (closestEndcapDetIdToProbe) && (theGeometry->getSubdetectorGeometry(closestEndcapDetIdToProbe)->
00252            getGeometry(closestEndcapDetIdToProbe)->inside(probe_gp)) )
00253      {
00254        double step = ((*itr).position() - (*(itr-1)).position()).mag();
00255        GlobalPoint point = itr->position();
00256        addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestEndcapDetIdToProbe, step, point, theEndcapSubdetGeometry);
00257        totalLengthCurved += step;
00258 
00259        if (firstPoint == false)
00260        {
00261          internalPointCurved = probe_gp ;
00262          firstPoint = true ;
00263        }
00264 
00265        externalPointCurved = probe_gp ;
00266      }
00267 
00268      if( (closestBarrelDetIdToProbe) && (theGeometry->getSubdetectorGeometry(closestBarrelDetIdToProbe)->
00269            getGeometry(closestBarrelDetIdToProbe)->inside(probe_gp)) )
00270      {
00271        double step = ((*itr).position() - (*(itr-1)).position()).mag();
00272        GlobalPoint point = itr->position();
00273        addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, closestBarrelDetIdToProbe, step, point, theBarrelSubdetGeometry);
00274        totalLengthCurved += step;
00275 
00276        if (firstPoint == false)
00277        {
00278          internalPointCurved = probe_gp ;
00279          firstPoint = true ;
00280        }
00281 
00282        externalPointCurved = probe_gp ;
00283      }
00284      else
00285      {
00286        // 3x3 matrix surrounding the probe
00287        surroundingMatrix = theTopology->getSubdetectorTopology(closestBarrelDetIdToProbe)->getWindow(closestBarrelDetIdToProbe,3,3);
00288 
00289        for( unsigned int k=0; k<surroundingMatrix.size(); ++k ) {
00290          if(theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k))->getGeometry(surroundingMatrix.at(k))->inside(probe_gp))
00291          {
00292            double step = ((*itr).position() - (*(itr-1)).position()).mag();
00293            GlobalPoint point = itr->position();
00294            addStepToXtal(trackExitPositionMap, trackCrossedXtalMap, surroundingMatrix[k], step,
00295                point, theGeometry->getSubdetectorGeometry(surroundingMatrix.at(k)));
00296            totalLengthCurved += step;
00297 
00298            if (firstPoint == false)
00299            {
00300              internalPointCurved = probe_gp ;
00301              firstPoint = true ;
00302            }
00303 
00304            externalPointCurved = probe_gp ;
00305          }
00306        }
00307 
00308        // clear neighborhood matrix
00309        surroundingMatrix.clear();
00310      }
00311    }
00312 
00313    return 0;
00314 }
00315 
00316 void BetaCalculatorECAL::addStepToXtal(std::map<int,GlobalPoint>& trackExitPositionMap,
00317     std::map<int,float>& trackCrossedXtalMap,
00318     DetId aDetId,
00319     float step,
00320     GlobalPoint point,
00321     const CaloSubdetectorGeometry* theSubdetGeometry)
00322 {
00323 
00324   const CaloCellGeometry *cell_p = theSubdetGeometry->getGeometry(aDetId);
00325   GlobalPoint p = (dynamic_cast <const TruncatedPyramid *> (cell_p))->getPosition(23);
00326   GlobalPoint diff(point.x()-p.x(),point.y()-p.y(),point.z()-p.z());
00327 
00328   std::map<int,GlobalPoint>::iterator xtal = trackExitPositionMap.find(aDetId.rawId());
00329   if (xtal!=trackExitPositionMap.end())
00330     ((*xtal).second)=diff;
00331   else
00332     trackExitPositionMap.insert(std::pair<int,GlobalPoint>(aDetId.rawId(),diff));
00333 
00334   std::map<int,float>::iterator xtal2 = trackCrossedXtalMap.find(aDetId.rawId());
00335   if (xtal2!= trackCrossedXtalMap.end())
00336     ((*xtal2).second)+=step;
00337   else
00338     trackCrossedXtalMap.insert(std::pair<int,float>(aDetId.rawId(),step));
00339 }
00340 
00341 
00342
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/SUSYBSMAnalysis/HSCP/src/BetaCalculatorECAL.cc
