00001 #include "CalibTracker/SiStripCommon/interface/ShallowSimhitClustersProducer.h"
00002
00003 #include "CalibTracker/SiStripCommon/interface/ShallowTools.h"
00004 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
00005 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
00006
00007 #include "MagneticField/Engine/interface/MagneticField.h"
00008 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00009 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
00010 #include "CondFormats/DataRecord/interface/SiStripLorentzAngleRcd.h"
00011 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00012 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00013 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00014 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "FWCore/Framework/interface/Event.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 #include "boost/foreach.hpp"
00019
00020 ShallowSimhitClustersProducer::ShallowSimhitClustersProducer(const edm::ParameterSet& iConfig)
00021 : inputTags( iConfig.getParameter<std::vector<edm::InputTag> >("InputTags") ),
00022 theClustersLabel( iConfig.getParameter<edm::InputTag>("Clusters")),
00023 Prefix( iConfig.getParameter<std::string>("Prefix") )
00024 {
00025 produces <std::vector<unsigned> > ( Prefix + "hits" );
00026 produces <std::vector<float> > ( Prefix + "strip" );
00027 produces <std::vector<float> > ( Prefix + "localtheta" );
00028 produces <std::vector<float> > ( Prefix + "localphi" );
00029 produces <std::vector<float> > ( Prefix + "localx" );
00030 produces <std::vector<float> > ( Prefix + "localy" );
00031 produces <std::vector<float> > ( Prefix + "localz" );
00032 produces <std::vector<float> > ( Prefix + "momentum" );
00033 produces <std::vector<float> > ( Prefix + "energyloss" );
00034 produces <std::vector<float> > ( Prefix + "time" );
00035 produces <std::vector<int> > ( Prefix + "particle" );
00036 produces <std::vector<unsigned short> > ( Prefix + "process" );
00037 }
00038
00039 void ShallowSimhitClustersProducer::
00040 produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
00041 shallow::CLUSTERMAP clustermap = shallow::make_cluster_map(iEvent,theClustersLabel);
00042
00043 int size = clustermap.size();
00044 std::auto_ptr<std::vector<unsigned> > hits ( new std::vector<unsigned> (size, 0) );
00045 std::auto_ptr<std::vector<float> > strip ( new std::vector<float> (size, -100) );
00046 std::auto_ptr<std::vector<float> > localtheta ( new std::vector<float> (size, -100) );
00047 std::auto_ptr<std::vector<float> > localphi ( new std::vector<float> (size, -100) );
00048 std::auto_ptr<std::vector<float> > localx ( new std::vector<float> (size, -100) );
00049 std::auto_ptr<std::vector<float> > localy ( new std::vector<float> (size, -100) );
00050 std::auto_ptr<std::vector<float> > localz ( new std::vector<float> (size, -100) );
00051 std::auto_ptr<std::vector<float> > momentum ( new std::vector<float> (size, 0) );
00052 std::auto_ptr<std::vector<float> > energyloss ( new std::vector<float> (size, -1) );
00053 std::auto_ptr<std::vector<float> > time ( new std::vector<float> (size, -1) );
00054 std::auto_ptr<std::vector<int> > particle ( new std::vector<int> (size,-500) );
00055 std::auto_ptr<std::vector<unsigned short> > process ( new std::vector<unsigned short> (size,0) );
00056
00057 edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
00058 edm::ESHandle<MagneticField> magfield; iSetup.get<IdealMagneticFieldRecord>().get(magfield);
00059 edm::ESHandle<SiStripLorentzAngle> SiStripLorentzAngle; iSetup.get<SiStripLorentzAngleRcd>().get(SiStripLorentzAngle);
00060 edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters; iEvent.getByLabel("siStripClusters", "", clusters);
00061
00062 BOOST_FOREACH( const edm::InputTag inputTag, inputTags ) { edm::Handle<std::vector<PSimHit> > simhits; iEvent.getByLabel(inputTag, simhits);
00063 BOOST_FOREACH( const PSimHit hit, *simhits ) {
00064
00065 const uint32_t id = hit.detUnitId();
00066 const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTrackerGeometry->idToDet( id ) );
00067 const LocalVector drift = shallow::drift(theStripDet, *magfield, *SiStripLorentzAngle);
00068
00069 const float driftedstrip_ = theStripDet->specificTopology().strip( hit.localPosition()+0.5*drift );
00070 const float hitstrip_ = theStripDet->specificTopology().strip( hit.localPosition() );
00071
00072 shallow::CLUSTERMAP::const_iterator cluster = match_cluster( id, driftedstrip_, clustermap, *clusters);
00073 if(cluster != clustermap.end()) {
00074 unsigned i = cluster->second;
00075 hits->at(i)+=1;
00076 if(hits->at(i) == 1) {
00077 strip->at(i) = hitstrip_;
00078 localtheta->at(i) = hit.thetaAtEntry();
00079 localphi->at(i) = hit.phiAtEntry();
00080 localx->at(i) = hit.localPosition().x();
00081 localy->at(i) = hit.localPosition().y();
00082 localz->at(i) = hit.localPosition().z();
00083 momentum->at(i) = hit.pabs();
00084 energyloss->at(i) = hit.energyLoss();
00085 time->at(i) = hit.timeOfFlight();
00086 particle->at(i) = hit.particleType();
00087 process->at(i) = hit.processType();
00088 }
00089 }
00090 }
00091 }
00092
00093 iEvent.put( hits, Prefix + "hits" );
00094 iEvent.put( strip, Prefix + "strip" );
00095 iEvent.put( localtheta, Prefix + "localtheta" );
00096 iEvent.put( localphi, Prefix + "localphi" );
00097 iEvent.put( localx, Prefix + "localx" );
00098 iEvent.put( localy, Prefix + "localy" );
00099 iEvent.put( localz, Prefix + "localz" );
00100 iEvent.put( momentum, Prefix + "momentum" );
00101 iEvent.put( energyloss, Prefix + "energyloss" );
00102 iEvent.put( time, Prefix + "time" );
00103 iEvent.put( particle, Prefix + "particle" );
00104 iEvent.put( process, Prefix + "process" );
00105 }
00106
00107 shallow::CLUSTERMAP::const_iterator ShallowSimhitClustersProducer::
00108 match_cluster( const unsigned& id, const float& strip_, const shallow::CLUSTERMAP& clustermap, const edmNew::DetSetVector<SiStripCluster>& clusters) const {
00109 shallow::CLUSTERMAP::const_iterator cluster = clustermap.end();
00110 edmNew::DetSetVector<SiStripCluster>::const_iterator clustersDetSet = clusters.find(id);
00111 if( clustersDetSet != clusters.end() ) {
00112 edmNew::DetSet<SiStripCluster>::const_iterator left, right=clustersDetSet->begin();
00113 while( right != clustersDetSet->end() && strip_ > right->barycenter() )
00114 right++;
00115 left = right-1;
00116 if(right!=clustersDetSet->end() && right!=clustersDetSet->begin()) {
00117 unsigned firstStrip = (right->barycenter()-strip_) < (strip_-left->barycenter()) ? right->firstStrip() : left->firstStrip();
00118 cluster = clustermap.find( std::make_pair( id, firstStrip));
00119 }
00120 else if(right != clustersDetSet->begin())
00121 cluster = clustermap.find( std::make_pair( id, left->firstStrip()));
00122 else
00123 cluster = clustermap.find( std::make_pair( id, right->firstStrip()));
00124 }
00125 return cluster;
00126 }
00127
1.7.3