| 1 |
#include "UserCode/L1RpcTriggerAnalysis/interface/DetHitCompatibleCollector.h"
|
| 2 |
|
| 3 |
#include "FWCore/Framework/interface/Event.h"
|
| 4 |
#include "FWCore/Framework/interface/EventSetup.h"
|
| 5 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 6 |
#include "DataFormats/Common/interface/Handle.h"
|
| 7 |
|
| 8 |
#include "DataFormats/MuonReco/interface/Muon.h"
|
| 9 |
#include "DataFormats/MuonReco/interface/MuonFwd.h"
|
| 10 |
|
| 11 |
//#include "Geometry/RPCGeometry/interface/RPCRoll.h"
|
| 12 |
//#include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
|
| 13 |
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
|
| 14 |
#include "DataFormats/MuonDetId/interface/RPCDetId.h"
|
| 15 |
#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
|
| 16 |
#include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
|
| 17 |
#include "DataFormats/DetId/interface/DetIdCollection.h"
|
| 18 |
#include "DataFormats/TrackReco/interface/HitPattern.h"
|
| 19 |
|
| 20 |
#include "Geometry/RPCGeometry/interface/RPCGeometry.h"
|
| 21 |
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
|
| 22 |
//#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
|
| 23 |
//#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
|
| 24 |
//#include "MagneticField/Engine/interface/MagneticField.h"
|
| 25 |
//#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
| 26 |
|
| 27 |
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
|
| 28 |
#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
|
| 29 |
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
|
| 30 |
|
| 31 |
#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
|
| 32 |
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
|
| 33 |
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
|
| 34 |
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
|
| 35 |
|
| 36 |
#include "DataFormats/Math/interface/deltaR.h"
|
| 37 |
|
| 38 |
#include "TObjArray.h"
|
| 39 |
#include "TH1F.h"
|
| 40 |
#include <sstream>
|
| 41 |
|
| 42 |
#include "UserCode/L1RpcTriggerAnalysis/interface/TrackAtSurface.h"
|
| 43 |
#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
|
| 44 |
|
| 45 |
|
| 46 |
DetHitCompatibleCollector::DetHitCompatibleCollector(const edm::ParameterSet& cfg)
|
| 47 |
: theNoDigiWarning(false)
|
| 48 |
{ }
|
| 49 |
|
| 50 |
DetHitCompatibleCollector::~ DetHitCompatibleCollector()
|
| 51 |
{
|
| 52 |
if (theNoDigiWarning) std::cout <<" **** DetHitCompatibleCollector **** WARNING - NoDigiWarning was set!" << std::endl;
|
| 53 |
}
|
| 54 |
|
| 55 |
std::vector<uint32_t> DetHitCompatibleCollector::compatibleSIMU( const reco::Muon* muon, const edm::Event &ev, const edm::EventSetup &es)
|
| 56 |
{
|
| 57 |
std::vector<uint32_t> detsSIMU;
|
| 58 |
edm::Handle< edm::PSimHitContainer > simHits;
|
| 59 |
ev.getByLabel("g4SimHits","MuonRPCHits",simHits);
|
| 60 |
std::cout <<" SIZE OF SIM HITS IS: "<<simHits->size() << std::endl;
|
| 61 |
for (edm::PSimHitContainer::const_iterator hitItr = simHits->begin(); hitItr != simHits->end(); ++hitItr) {
|
| 62 |
uint32_t rpcRawId = hitItr->detUnitId();
|
| 63 |
if ( !(13 == abs(hitItr->particleType()) ) ) continue;
|
| 64 |
if (detsSIMU.end() == find(detsSIMU.begin(),detsSIMU.end(),rpcRawId) ) detsSIMU.push_back(rpcRawId);
|
| 65 |
}
|
| 66 |
return detsSIMU;
|
| 67 |
}
|
| 68 |
|
| 69 |
std::vector<DetCluDigiObj> DetHitCompatibleCollector::compatibleHits( const reco::Muon* muon, const edm::Event &ev, const edm::EventSetup &es)
|
| 70 |
{
|
| 71 |
std::vector<DetCluDigiObj> detsHitsCompatibleWithMuon;
|
| 72 |
|
| 73 |
edm::Handle<RPCRecHitCollection> recHits;
|
| 74 |
ev.getByLabel("rpcRecHits", recHits);
|
| 75 |
|
| 76 |
edm::Handle<RPCDigiCollection> rpcDigis;
|
| 77 |
ev.getByLabel("muonRPCDigis", rpcDigis);
|
| 78 |
|
| 79 |
edm::ESHandle<RPCGeometry> rpcGeometry;
|
| 80 |
es.get<MuonGeometryRecord>().get(rpcGeometry);
|
| 81 |
|
| 82 |
TrackAtSurface trackAtSurface(muon, ev,es);
|
| 83 |
|
| 84 |
std::map< uint32_t, std::pair<uint32_t,uint32_t> > aMap;
|
| 85 |
typedef RPCRecHitCollection::const_iterator IH;
|
| 86 |
for (IH ih=recHits->begin(); ih != recHits->end(); ++ih) {
|
| 87 |
|
| 88 |
if (ih->BunchX() != 0) continue;
|
| 89 |
if (! ih->isValid() ) continue;
|
| 90 |
RPCDetId rpcDet = ih->rpcId();
|
| 91 |
GlobalPoint detPosition = rpcGeometry->idToDet(rpcDet)->position();
|
| 92 |
GlobalPoint hitPosition = rpcGeometry->idToDet(rpcDet)->toGlobal(ih->localPosition());
|
| 93 |
if (fabs(reco::deltaPhi(detPosition.phi(),muon->phi())) > M_PI/2.) continue;
|
| 94 |
|
| 95 |
/*
|
| 96 |
TrajectoryStateOnSurface trackAtHit= trackAtSurface.atPoint(hitPosition);
|
| 97 |
if(!trackAtHit.isValid()) std::cout <<" TRAJ NOT VALID! "<< std::endl;
|
| 98 |
if (!trackAtHit.isValid()) continue;
|
| 99 |
GlobalPoint hitPoint = rpcGeometry->idToDet(rpcDet)->toGlobal(ih->localPosition());
|
| 100 |
LocalError hitError = ih->localPositionError();
|
| 101 |
GlobalPoint trackAtHitPoint = trackAtHit.globalPosition();
|
| 102 |
LocalError trackAtHitError = trackAtHit.localError().positionError();
|
| 103 |
float distX = hitPoint.perp()*deltaPhi((float) hitPoint.phi(), (float) trackAtHitPoint.phi());
|
| 104 |
float distY = RPCDetIdUtil(rpcDet).isBarrel() ? hitPoint.z()-trackAtHitPoint.z() : hitPoint.perp()-trackAtHitPoint.perp() ;
|
| 105 |
*/
|
| 106 |
//
|
| 107 |
TrajectoryStateOnSurface trackAtHit= trackAtSurface.atDetFromClose(rpcDet,hitPosition);
|
| 108 |
// if(!trackAtHit.isValid()) std::cout <<" TRAJ NOT VALID! "<< std::endl;
|
| 109 |
if (!trackAtHit.isValid()) continue;
|
| 110 |
LocalPoint hitPoint = ih->localPosition();
|
| 111 |
LocalError hitError = ih->localPositionError();
|
| 112 |
LocalPoint trackAtHitPoint = trackAtHit.localPosition() ;
|
| 113 |
LocalError trackAtHitError = trackAtHit.localError().positionError();
|
| 114 |
float distX = hitPoint.x()-trackAtHitPoint.x();
|
| 115 |
float distY = hitPoint.y()-trackAtHitPoint.y();
|
| 116 |
//
|
| 117 |
|
| 118 |
float pullX = distX/ sqrt( trackAtHitError.xx()+hitError.xx());
|
| 119 |
float pullY = distY/ sqrt( trackAtHitError.yy()+hitError.yy());
|
| 120 |
|
| 121 |
bool hitCompatible = (fabs(pullX) < 3.5 || fabs(distX) < 10.) && fabs(pullY) < 3.5 ;
|
| 122 |
if (hitCompatible) {
|
| 123 |
if (aMap.find(rpcDet.rawId()) == aMap.end() ) aMap[rpcDet.rawId()] = std::make_pair(0,0);
|
| 124 |
|
| 125 |
unsigned int clusterSize = ih->clusterSize();
|
| 126 |
if ( aMap[rpcDet.rawId()].first < clusterSize) aMap[rpcDet.rawId()].first = clusterSize;
|
| 127 |
|
| 128 |
if (rpcDigis.isValid()) {
|
| 129 |
const RPCDigiCollection::Range range = rpcDigis->get( rpcDet.rawId() );
|
| 130 |
std::map<int, bool> strips;
|
| 131 |
for (RPCDigiCollection::const_iterator id = range.first; id != range.second; ++id) if (id->bx() == 0) strips[id->strip()] = true;
|
| 132 |
if ( strips.size() == 0 ) std::cout <<"WARNING ***************"<<std::endl;
|
| 133 |
aMap[rpcDet.rawId()].second = strips.size();
|
| 134 |
} else theNoDigiWarning = true;
|
| 135 |
}
|
| 136 |
|
| 137 |
RPCDetIdUtil place(rpcDet);
|
| 138 |
if (place.isBarrel()) {
|
| 139 |
hPullX_B[place.layer()-1]->Fill(pullX);
|
| 140 |
hDistX_B[place.layer()-1]->Fill(distX);
|
| 141 |
hPullY_B[place.layer()-1]->Fill(pullY);
|
| 142 |
hDistY_B[place.layer()-1]->Fill(distY);
|
| 143 |
} else {
|
| 144 |
hPullX_E[place.layer()-1]->Fill(pullX);
|
| 145 |
hDistX_E[place.layer()-1]->Fill(distX);
|
| 146 |
hPullY_E[place.layer()-1]->Fill(pullY);
|
| 147 |
hDistY_E[place.layer()-1]->Fill(distY);
|
| 148 |
}
|
| 149 |
hPullX->Fill(pullX);
|
| 150 |
hPullY->Fill(pullY);
|
| 151 |
}
|
| 152 |
|
| 153 |
|
| 154 |
for ( std::map< uint32_t, std::pair<uint32_t,uint32_t> >::const_iterator it = aMap.begin(); it != aMap.end(); ++it) {
|
| 155 |
DetCluDigiObj obj;
|
| 156 |
obj.det = (*it).first;
|
| 157 |
obj.clusterSize = (*it).second.first;
|
| 158 |
obj.nDigis = (*it).second.second;
|
| 159 |
detsHitsCompatibleWithMuon.push_back(obj);
|
| 160 |
}
|
| 161 |
return detsHitsCompatibleWithMuon;
|
| 162 |
}
|
| 163 |
|
| 164 |
/*
|
| 165 |
std::vector<uint32_t> DetHitCompatibleCollector::clSizeCompDets(const std::vector<uint32_t> & detIds, const edm::Event &ev, const edm::EventSetup &es)
|
| 166 |
{
|
| 167 |
std::vector<uint32_t> cluInDets;
|
| 168 |
edm::Handle<RPCRecHitCollection> recHits;
|
| 169 |
ev.getByLabel("rpcRecHits", recHits);
|
| 170 |
for (std::vector<uint32_t>::const_iterator it = detIds.begin(); it != detIds.end(); ++it) {
|
| 171 |
RPCRecHitCollection::range range = recHits->get(*it);
|
| 172 |
// std::cout <<"DET: "<<*it << std::endl;
|
| 173 |
for(RPCRecHitCollection::const_iterator ih= range.first; ih != range.second; ++ih) {
|
| 174 |
if (ih->BunchX() == 0) { cluInDets.push_back(ih->clusterSize()); break; }
|
| 175 |
}
|
| 176 |
// for(RPCRecHitCollection::const_iterator ih= range.first; ih != range.second; ++ih)
|
| 177 |
// std::cout <<" Hit: cluster:"<<ih->clusterSize() << "first strip: "<<ih->firstClusterStrip()<<std::endl;
|
| 178 |
// if (range.second-range.first) cluInDets.push_back(range.first->clusterSize());
|
| 179 |
}
|
| 180 |
return cluInDets;
|
| 181 |
}
|
| 182 |
|
| 183 |
std::vector<uint32_t> DetHitCompatibleCollector::nDigisCompDets(const std::vector<uint32_t> & detIds, const edm::Event &ev, const edm::EventSetup &es)
|
| 184 |
{
|
| 185 |
std::vector<uint32_t> digisInDets;
|
| 186 |
edm::Handle<RPCDigiCollection> rpcDigis;
|
| 187 |
ev.getByLabel("muonRPCDigis", rpcDigis);
|
| 188 |
for (std::vector<uint32_t>::const_iterator it = detIds.begin(); it != detIds.end(); ++it) {
|
| 189 |
// std::cout <<"DET: "<<*it << std::endl;
|
| 190 |
const RPCDigiCollection::Range range = rpcDigis->get(*it);
|
| 191 |
std::map<int, bool> strips;
|
| 192 |
// for (RPCDigiCollection::const_iterator id = range.first; id != range.second; ++id) id->print() ;
|
| 193 |
for (RPCDigiCollection::const_iterator id = range.first; id != range.second; ++id) if (id->bx() == 0) strips[id->strip()] = true;
|
| 194 |
if ( strips.size() == 0 ) std::cout <<"WARNING ***************"<<std::endl;
|
| 195 |
digisInDets.push_back( strips.size());
|
| 196 |
// std::cout <<"SIZE is : "<<strips.size() << std::endl;
|
| 197 |
}
|
| 198 |
return digisInDets;
|
| 199 |
}
|
| 200 |
*/
|
| 201 |
|
| 202 |
std::vector<uint32_t> DetHitCompatibleCollector::compatibleDets( const reco::Muon* muon, const edm::Event &ev, const edm::EventSetup &es, bool deepInside)
|
| 203 |
{
|
| 204 |
std::vector<uint32_t> detsCrossedByMuon;
|
| 205 |
TrackAtSurface trackAtSurface(muon, ev,es);
|
| 206 |
edm::ESHandle<RPCGeometry> rpcGeometry;
|
| 207 |
es.get<MuonGeometryRecord>().get(rpcGeometry);
|
| 208 |
const TrackingGeometry::DetIdContainer & detIds = rpcGeometry->detIds();
|
| 209 |
for (TrackingGeometry::DetIdContainer::const_iterator it = detIds.begin(); it != detIds.end(); ++it) {
|
| 210 |
|
| 211 |
const GeomDet * det = rpcGeometry->idToDet(*it);
|
| 212 |
GlobalPoint detPosition = det->position();
|
| 213 |
RPCDetId rpcDet(*it);
|
| 214 |
if (rpcDet.roll()==0) continue;
|
| 215 |
if (deltaR(muon->eta(), muon->phi(), detPosition.eta(), detPosition.phi()) > 0.5) continue;
|
| 216 |
|
| 217 |
// TrajectoryStateOnSurface stateAtDet = deepInside ? trackAtSurface.atDetFromClose(rpcDet,detPosition) : trackAtSurface.atDetFromTrack(rpcDet);
|
| 218 |
TrajectoryStateOnSurface stateAtDet = trackAtSurface.atDetFromClose(rpcDet,detPosition);
|
| 219 |
// : trackAtSurface.atDetFromTrack(rpcDet);
|
| 220 |
if (!stateAtDet.isValid()) continue;
|
| 221 |
// std::cout <<" det:width="<< det->surface().bounds().width()/2.<<" det:lengths="<<det->surface().bounds().length()/2.<<std::endl;
|
| 222 |
if (! (det->surface().bounds().inside(stateAtDet.localPosition()))) continue;
|
| 223 |
if (deepInside) {
|
| 224 |
LocalError stateAtDetError = stateAtDet.localError().positionError();
|
| 225 |
LocalPoint farEdge( fabs(stateAtDet.localPosition().x())+3.5*sqrt(stateAtDetError.xx())+5., fabs(stateAtDet.localPosition().y())+3.5*sqrt(stateAtDetError.yy())+5.);
|
| 226 |
if (! det->surface().bounds().inside(farEdge) ) continue;
|
| 227 |
}
|
| 228 |
if (detsCrossedByMuon.end() == find(detsCrossedByMuon.begin(),detsCrossedByMuon.end(),rpcDet.rawId()) )
|
| 229 |
detsCrossedByMuon.push_back(rpcDet.rawId());
|
| 230 |
|
| 231 |
hPropToDetDeltaR->Fill(deltaR(muon->eta(), muon->phi(), detPosition.eta(), detPosition.phi()));
|
| 232 |
}
|
| 233 |
return detsCrossedByMuon;
|
| 234 |
}
|
| 235 |
|
| 236 |
|
| 237 |
void DetHitCompatibleCollector::initHistos( TObjArray & histos)
|
| 238 |
{
|
| 239 |
for (unsigned int i=1; i<=6; ++i) {
|
| 240 |
std::stringstream namePullX_B; namePullX_B <<"hPullX_B"<< i;
|
| 241 |
std::stringstream nameDistX_B; nameDistX_B <<"hDistX_B"<< i;
|
| 242 |
hPullX_B[i-1]= new TH1F(namePullX_B.str().c_str(),namePullX_B.str().c_str(),50,-5.,5.); histos.Add( hPullX_B[i-1]);
|
| 243 |
hDistX_B[i-1]= new TH1F(nameDistX_B.str().c_str(),nameDistX_B.str().c_str(),100,-20.,20.); histos.Add( hDistX_B[i-1]);
|
| 244 |
std::stringstream namePullY_B; namePullY_B <<"hPullY_B"<< i;
|
| 245 |
std::stringstream nameDistY_B; nameDistY_B <<"hDistY_B"<< i;
|
| 246 |
hPullY_B[i-1]= new TH1F(namePullY_B.str().c_str(),namePullY_B.str().c_str(),50,-5.,5.); histos.Add( hPullY_B[i-1]);
|
| 247 |
hDistY_B[i-1]= new TH1F(nameDistY_B.str().c_str(),nameDistY_B.str().c_str(),100,-100.,100.); histos.Add( hDistY_B[i-1]);
|
| 248 |
}
|
| 249 |
for (unsigned int i=1; i<=3; ++i) {
|
| 250 |
std::stringstream namePullX_E; namePullX_E <<"hPullX_E" << i;
|
| 251 |
std::stringstream nameDistX_E; nameDistX_E <<"hDistX_E"<< i;
|
| 252 |
hPullX_E[i-1]= new TH1F(namePullX_E.str().c_str(),namePullX_E.str().c_str(),50,-5.,5.); histos.Add( hPullX_E[i-1]);
|
| 253 |
hDistX_E[i-1]= new TH1F(nameDistX_E.str().c_str(),nameDistX_E.str().c_str(),100,-20.,20.); histos.Add( hDistX_E[i-1]);
|
| 254 |
std::stringstream namePullY_E; namePullY_E <<"hPullY_E" << i;
|
| 255 |
std::stringstream nameDistY_E; nameDistY_E <<"hDistY_E"<< i;
|
| 256 |
hPullY_E[i-1]= new TH1F(namePullY_E.str().c_str(),namePullY_E.str().c_str(),50,-5.,5.); histos.Add( hPullY_E[i-1]);
|
| 257 |
hDistY_E[i-1]= new TH1F(nameDistY_E.str().c_str(),nameDistY_E.str().c_str(),100,-100.,100.); histos.Add( hDistY_E[i-1]);
|
| 258 |
}
|
| 259 |
hPullX = new TH1F("hPullX","hPullX",50,-5.,5.); histos.Add( hPullX);
|
| 260 |
hPullY = new TH1F("hPullY","hPullY",50,-5.,5.); histos.Add( hPullY);
|
| 261 |
|
| 262 |
hPropToDetDeltaR = new TH1F("hPropToDetDeltaR","DR good Mu-DU before propagation",100,0.,3.); histos.Add(hPropToDetDeltaR);
|
| 263 |
}
|
