| 1 |
#include "UserCode/ShallowTools/interface/ShallowSimhitClustersProducer.h"
|
| 2 |
|
| 3 |
#include "UserCode/ShallowTools/interface/ShallowTools.h"
|
| 4 |
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
|
| 5 |
#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
|
| 6 |
|
| 7 |
#include "MagneticField/Engine/interface/MagneticField.h"
|
| 8 |
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
| 9 |
#include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
|
| 10 |
#include "CondFormats/DataRecord/interface/SiStripLorentzAngleRcd.h"
|
| 11 |
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
|
| 12 |
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
|
| 13 |
#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
|
| 14 |
#include "Geometry/CommonTopologies/interface/StripTopology.h"
|
| 15 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 16 |
#include "FWCore/Framework/interface/Event.h"
|
| 17 |
#include "boost/foreach.hpp"
|
| 18 |
|
| 19 |
ShallowSimhitClustersProducer::ShallowSimhitClustersProducer(const edm::ParameterSet& iConfig)
|
| 20 |
: inputTags( iConfig.getParameter<std::vector<edm::InputTag> >("InputTags") ),
|
| 21 |
theClustersLabel( iConfig.getParameter<edm::InputTag>("Clusters")),
|
| 22 |
Prefix( iConfig.getParameter<std::string>("Prefix") )
|
| 23 |
{
|
| 24 |
produces <std::vector<unsigned> > ( Prefix + "hits" );
|
| 25 |
produces <std::vector<float> > ( Prefix + "strip" );
|
| 26 |
produces <std::vector<float> > ( Prefix + "localtheta" );
|
| 27 |
produces <std::vector<float> > ( Prefix + "localphi" );
|
| 28 |
produces <std::vector<float> > ( Prefix + "localx" );
|
| 29 |
produces <std::vector<float> > ( Prefix + "localy" );
|
| 30 |
produces <std::vector<float> > ( Prefix + "localz" );
|
| 31 |
produces <std::vector<float> > ( Prefix + "momentum" );
|
| 32 |
produces <std::vector<float> > ( Prefix + "energyloss" );
|
| 33 |
produces <std::vector<float> > ( Prefix + "time" );
|
| 34 |
produces <std::vector<int> > ( Prefix + "particle" );
|
| 35 |
produces <std::vector<unsigned short> > ( Prefix + "process" );
|
| 36 |
}
|
| 37 |
|
| 38 |
void ShallowSimhitClustersProducer::
|
| 39 |
produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
|
| 40 |
shallow::CLUSTERMAP clustermap = shallow::make_cluster_map(iEvent,theClustersLabel);
|
| 41 |
|
| 42 |
int size = clustermap.size();
|
| 43 |
std::auto_ptr<std::vector<unsigned> > hits ( new std::vector<unsigned> (size, 0) );
|
| 44 |
std::auto_ptr<std::vector<float> > strip ( new std::vector<float> (size, -100) );
|
| 45 |
std::auto_ptr<std::vector<float> > localtheta ( new std::vector<float> (size, -100) );
|
| 46 |
std::auto_ptr<std::vector<float> > localphi ( new std::vector<float> (size, -100) );
|
| 47 |
std::auto_ptr<std::vector<float> > localx ( new std::vector<float> (size, -100) );
|
| 48 |
std::auto_ptr<std::vector<float> > localy ( new std::vector<float> (size, -100) );
|
| 49 |
std::auto_ptr<std::vector<float> > localz ( new std::vector<float> (size, -100) );
|
| 50 |
std::auto_ptr<std::vector<float> > momentum ( new std::vector<float> (size, 0) );
|
| 51 |
std::auto_ptr<std::vector<float> > energyloss ( new std::vector<float> (size, -1) );
|
| 52 |
std::auto_ptr<std::vector<float> > time ( new std::vector<float> (size, -1) );
|
| 53 |
std::auto_ptr<std::vector<int> > particle ( new std::vector<int> (size,-500) );
|
| 54 |
std::auto_ptr<std::vector<unsigned short> > process ( new std::vector<unsigned short> (size,0) );
|
| 55 |
|
| 56 |
edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
|
| 57 |
edm::ESHandle<MagneticField> magfield; iSetup.get<IdealMagneticFieldRecord>().get(magfield);
|
| 58 |
edm::ESHandle<SiStripLorentzAngle> SiStripLorentzAngle; iSetup.get<SiStripLorentzAngleRcd>().get(SiStripLorentzAngle);
|
| 59 |
edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters; iEvent.getByLabel("siStripClusters", "", clusters);
|
| 60 |
|
| 61 |
BOOST_FOREACH( const edm::InputTag inputTag, inputTags ) { edm::Handle<std::vector<PSimHit> > simhits; iEvent.getByLabel(inputTag, simhits);
|
| 62 |
BOOST_FOREACH( const PSimHit hit, *simhits ) {
|
| 63 |
|
| 64 |
const uint32_t id = hit.detUnitId();
|
| 65 |
const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTrackerGeometry->idToDet( id ) );
|
| 66 |
const LocalVector drift = shallow::drift(theStripDet, *magfield, *SiStripLorentzAngle);
|
| 67 |
|
| 68 |
const float driftedstrip_ = theStripDet->specificTopology().strip( hit.localPosition()+0.5*drift );
|
| 69 |
const float hitstrip_ = theStripDet->specificTopology().strip( hit.localPosition() );
|
| 70 |
|
| 71 |
shallow::CLUSTERMAP::const_iterator cluster = match_cluster( id, driftedstrip_, clustermap, *clusters);
|
| 72 |
if(cluster != clustermap.end()) {
|
| 73 |
unsigned i = cluster->second;
|
| 74 |
hits->at(i)+=1;
|
| 75 |
if(hits->at(i) == 1) {
|
| 76 |
strip->at(i) = hitstrip_;
|
| 77 |
localtheta->at(i) = hit.thetaAtEntry();
|
| 78 |
localphi->at(i) = hit.phiAtEntry();
|
| 79 |
localx->at(i) = hit.localPosition().x();
|
| 80 |
localy->at(i) = hit.localPosition().y();
|
| 81 |
localz->at(i) = hit.localPosition().z();
|
| 82 |
momentum->at(i) = hit.pabs();
|
| 83 |
energyloss->at(i) = hit.energyLoss();
|
| 84 |
time->at(i) = hit.timeOfFlight();
|
| 85 |
particle->at(i) = hit.particleType();
|
| 86 |
process->at(i) = hit.processType();
|
| 87 |
}
|
| 88 |
}
|
| 89 |
}
|
| 90 |
}
|
| 91 |
|
| 92 |
iEvent.put( hits, Prefix + "hits" );
|
| 93 |
iEvent.put( strip, Prefix + "strip" );
|
| 94 |
iEvent.put( localtheta, Prefix + "localtheta" );
|
| 95 |
iEvent.put( localphi, Prefix + "localphi" );
|
| 96 |
iEvent.put( localx, Prefix + "localx" );
|
| 97 |
iEvent.put( localy, Prefix + "localy" );
|
| 98 |
iEvent.put( localz, Prefix + "localz" );
|
| 99 |
iEvent.put( momentum, Prefix + "momentum" );
|
| 100 |
iEvent.put( energyloss, Prefix + "energyloss" );
|
| 101 |
iEvent.put( time, Prefix + "time" );
|
| 102 |
iEvent.put( particle, Prefix + "particle" );
|
| 103 |
iEvent.put( process, Prefix + "process" );
|
| 104 |
}
|
| 105 |
|
| 106 |
shallow::CLUSTERMAP::const_iterator ShallowSimhitClustersProducer::
|
| 107 |
match_cluster( const unsigned& id, const float& strip_, const shallow::CLUSTERMAP& clustermap, const edmNew::DetSetVector<SiStripCluster>& clusters) const {
|
| 108 |
shallow::CLUSTERMAP::const_iterator cluster = clustermap.end();
|
| 109 |
edmNew::DetSetVector<SiStripCluster>::const_iterator clustersDetSet = clusters.find(id);
|
| 110 |
if( clustersDetSet != clusters.end() ) {
|
| 111 |
edmNew::DetSet<SiStripCluster>::const_iterator left, right=clustersDetSet->begin();
|
| 112 |
while( right != clustersDetSet->end() && strip_ > right->barycenter() )
|
| 113 |
right++;
|
| 114 |
left = right-1;
|
| 115 |
if(right!=clustersDetSet->end() && right!=clustersDetSet->begin()) {
|
| 116 |
unsigned firstStrip = (right->barycenter()-strip_) < (strip_-left->barycenter()) ? right->firstStrip() : left->firstStrip();
|
| 117 |
cluster = clustermap.find( std::make_pair( id, firstStrip));
|
| 118 |
}
|
| 119 |
else if(right != clustersDetSet->begin())
|
| 120 |
cluster = clustermap.find( std::make_pair( id, left->firstStrip()));
|
| 121 |
else
|
| 122 |
cluster = clustermap.find( std::make_pair( id, right->firstStrip()));
|
| 123 |
}
|
| 124 |
return cluster;
|
| 125 |
}
|
| 126 |
|
