| 1 |
// -*- C++ -*-
|
| 2 |
//
|
| 3 |
// Package: TrackingAnalysis/Cosmics
|
| 4 |
// Class: HitRes
|
| 5 |
//
|
| 6 |
/**\class TrackingAnalysis/Cosmics HitRes.cc TrackingAnalysis/Cosmics/src/HitRes.cc
|
| 7 |
|
| 8 |
Description: <one line class summary>
|
| 9 |
Use overlaps in TIF cosmics data to evaluate hit resolution
|
| 10 |
|
| 11 |
Implementation:
|
| 12 |
See sample cfg files in TrackingAnalysis/Cosmics/test/hitRes*cfg
|
| 13 |
*/
|
| 14 |
//
|
| 15 |
// Original Authors: Wolfgang Adam, Keith Ulmer
|
| 16 |
// Created: Thu Oct 11 14:53:32 CEST 2007
|
| 17 |
// $Id: HitRes.cc,v 1.18 2010/05/19 15:31:58 speer Exp $
|
| 18 |
//
|
| 19 |
//
|
| 20 |
|
| 21 |
|
| 22 |
// system include files
|
| 23 |
#include <memory>
|
| 24 |
|
| 25 |
// user include files
|
| 26 |
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
| 27 |
#include "FWCore/Framework/interface/EDAnalyzer.h"
|
| 28 |
|
| 29 |
#include "FWCore/Framework/interface/Event.h"
|
| 30 |
#include "FWCore/Framework/interface/EventSetup.h"
|
| 31 |
#include "FWCore/Framework/interface/MakerMacros.h"
|
| 32 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 33 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 34 |
|
| 35 |
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
| 36 |
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
|
| 37 |
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
|
| 38 |
#include "MagneticField/Engine/interface/MagneticField.h"
|
| 39 |
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
| 40 |
#include "FWCore/ServiceRegistry/interface/Service.h"
|
| 41 |
#include "CommonTools/UtilAlgos/interface/TFileService.h"
|
| 42 |
|
| 43 |
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
|
| 44 |
#include "FWCore/ParameterSet/interface/FileInPath.h"
|
| 45 |
#include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
|
| 46 |
|
| 47 |
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
|
| 48 |
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
|
| 49 |
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
|
| 50 |
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
|
| 51 |
#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h"
|
| 52 |
|
| 53 |
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
| 54 |
#include "DataFormats/TrackReco/interface/Track.h"
|
| 55 |
#include "DataFormats/TrackReco/interface/TrackExtra.h"
|
| 56 |
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
|
| 57 |
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
|
| 58 |
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
|
| 59 |
#include "DataFormats/DetId/interface/DetId.h"
|
| 60 |
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
|
| 61 |
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
|
| 62 |
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
|
| 63 |
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
|
| 64 |
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
|
| 65 |
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
|
| 66 |
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
|
| 67 |
#include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
|
| 68 |
#include "TrackingTools/PatternTools/interface/Trajectory.h"
|
| 69 |
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
|
| 70 |
#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripMatchedRecHit.h"
|
| 71 |
#include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
|
| 72 |
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
|
| 73 |
#include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
|
| 74 |
|
| 75 |
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
|
| 76 |
#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
|
| 77 |
#include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
|
| 78 |
#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
|
| 79 |
#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
|
| 80 |
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
|
| 81 |
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
|
| 82 |
|
| 83 |
#include "TFile.h"
|
| 84 |
#include "TTree.h"
|
| 85 |
|
| 86 |
#include <vector>
|
| 87 |
#include <utility>
|
| 88 |
|
| 89 |
using namespace std;
|
| 90 |
//
|
| 91 |
// class decleration
|
| 92 |
//
|
| 93 |
|
| 94 |
|
| 95 |
class HitRes : public edm::EDAnalyzer {
|
| 96 |
public:
|
| 97 |
explicit HitRes(const edm::ParameterSet&);
|
| 98 |
~HitRes();
|
| 99 |
|
| 100 |
private:
|
| 101 |
typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
|
| 102 |
virtual void analyze(const edm::Event&, const edm::EventSetup&);
|
| 103 |
virtual void endJob() ;
|
| 104 |
|
| 105 |
virtual void analyze(const Trajectory&, const Propagator&, TrackerHitAssociator&);
|
| 106 |
int layerFromId (const DetId&) const;
|
| 107 |
|
| 108 |
// ----------member data ---------------------------
|
| 109 |
edm::ParameterSet config_;
|
| 110 |
edm::InputTag trajectoryTag_;
|
| 111 |
SiStripDetInfoFileReader* reader;
|
| 112 |
bool doSimHit_;
|
| 113 |
const TrackerGeometry* trackerGeometry_;
|
| 114 |
const MagneticField* magField_;
|
| 115 |
|
| 116 |
TrajectoryStateCombiner combiner_;
|
| 117 |
int overlapCounts_[3];
|
| 118 |
|
| 119 |
TTree* rootTree_;
|
| 120 |
edm::FileInPath FileInPath_;
|
| 121 |
float overlapPath_;
|
| 122 |
uint layer_;
|
| 123 |
unsigned short hitCounts_[2];
|
| 124 |
float chi2_[2];
|
| 125 |
unsigned int overlapIds_[2];
|
| 126 |
float predictedPositions_[3][2];
|
| 127 |
float predictedLocalParameters_[5][2];
|
| 128 |
float predictedLocalErrors_[5][2];
|
| 129 |
float predictedDeltaXError_;
|
| 130 |
float predictedDeltaYError_;
|
| 131 |
char relativeXSign_;
|
| 132 |
char relativeYSign_;
|
| 133 |
float hitPositions_[2];
|
| 134 |
float hitErrors_[2];
|
| 135 |
float hitPositionsY_[2];
|
| 136 |
float hitErrorsY_[2];
|
| 137 |
float simHitPositions_[2];
|
| 138 |
float simHitPositionsY_[2];
|
| 139 |
float clusterWidthX_[2];
|
| 140 |
float clusterWidthY_[2];
|
| 141 |
float clusterSize_[2];
|
| 142 |
uint clusterCharge_[2];
|
| 143 |
int edge_[2];
|
| 144 |
|
| 145 |
vector<bool> acceptLayer;
|
| 146 |
float momentum_;
|
| 147 |
uint run_, event_;
|
| 148 |
bool barrelOnly_;
|
| 149 |
};
|
| 150 |
|
| 151 |
//
|
| 152 |
// constants, enums and typedefs
|
| 153 |
//
|
| 154 |
|
| 155 |
//
|
| 156 |
// static data member definitions
|
| 157 |
//
|
| 158 |
|
| 159 |
using std::vector;
|
| 160 |
using std::cout;
|
| 161 |
using std::endl;
|
| 162 |
//
|
| 163 |
// constructors and destructor
|
| 164 |
//
|
| 165 |
HitRes::HitRes(const edm::ParameterSet& iConfig) :
|
| 166 |
config_(iConfig), rootTree_(0),
|
| 167 |
FileInPath_("CalibTracker/SiStripCommon/data/SiStripDetInfo.dat")
|
| 168 |
{
|
| 169 |
//now do what ever initialization is needed
|
| 170 |
trajectoryTag_ = iConfig.getParameter<edm::InputTag>("trajectories");
|
| 171 |
doSimHit_ = iConfig.getParameter<bool>("associateStrip");
|
| 172 |
reader=new SiStripDetInfoFileReader(FileInPath_.fullPath());
|
| 173 |
|
| 174 |
overlapCounts_[0] = 0; // #trajectories
|
| 175 |
overlapCounts_[1] = 0; // #hits
|
| 176 |
overlapCounts_[2] = 0; // #overlap hits
|
| 177 |
acceptLayer.resize(7,false);
|
| 178 |
acceptLayer[PixelSubdetector::PixelBarrel] = iConfig.getParameter<bool>("usePXB") ;
|
| 179 |
acceptLayer[PixelSubdetector::PixelEndcap] = iConfig.getParameter<bool>("usePXF") ;
|
| 180 |
acceptLayer[StripSubdetector::TIB] = iConfig.getParameter<bool>("useTIB") ;
|
| 181 |
acceptLayer[StripSubdetector::TOB] = iConfig.getParameter<bool>("useTOB") ;
|
| 182 |
acceptLayer[StripSubdetector::TID] = iConfig.getParameter<bool>("useTID") ;
|
| 183 |
acceptLayer[StripSubdetector::TEC] = iConfig.getParameter<bool>("useTEC") ;
|
| 184 |
barrelOnly_ = iConfig.getParameter<bool>("barrelOnly");
|
| 185 |
|
| 186 |
edm::Service<TFileService> fs;
|
| 187 |
//
|
| 188 |
// root output
|
| 189 |
//
|
| 190 |
rootTree_ = fs->make<TTree>("Overlaps","Overlaps");
|
| 191 |
rootTree_->Branch("hitCounts",hitCounts_,"found/s:lost/s");
|
| 192 |
rootTree_->Branch("chi2",chi2_,"chi2/F:ndf/F");
|
| 193 |
rootTree_->Branch("path",&overlapPath_,"path/F");
|
| 194 |
rootTree_->Branch("layer",&layer_,"layer/i");
|
| 195 |
rootTree_->Branch("detids",overlapIds_,"id[2]/i");
|
| 196 |
rootTree_->Branch("predPos",predictedPositions_,"gX[2]/F:gY[2]/F:gZ[2]/F");
|
| 197 |
rootTree_->Branch("predPar",predictedLocalParameters_,
|
| 198 |
"predQP[2]/F:predDX[2]/F:predDY[2]/F:predX[2]/F:predY[2]/F");
|
| 199 |
rootTree_->Branch("predErr",predictedLocalErrors_,
|
| 200 |
"predEQP[2]/F:predEDX[2]/F:predEDY[2]/F:predEX[2]/F:predEY[2]/F");
|
| 201 |
rootTree_->Branch("predEDeltaX",&predictedDeltaXError_,"sigDeltaX/F");
|
| 202 |
rootTree_->Branch("predEDeltaY",&predictedDeltaYError_,"sigDeltaY/F");
|
| 203 |
rootTree_->Branch("relSignX",&relativeXSign_,"relSignX/B");
|
| 204 |
rootTree_->Branch("relSignY",&relativeYSign_,"relSignY/B");
|
| 205 |
rootTree_->Branch("hitX",hitPositions_,"hitX[2]/F");
|
| 206 |
rootTree_->Branch("hitEX",hitErrors_,"hitEX[2]/F");
|
| 207 |
rootTree_->Branch("hitY",hitPositionsY_,"hitY[2]/F");
|
| 208 |
rootTree_->Branch("hitEY",hitErrorsY_,"hitEY[2]/F");
|
| 209 |
rootTree_->Branch("simX",simHitPositions_,"simX[2]/F");
|
| 210 |
rootTree_->Branch("simY",simHitPositionsY_,"simY[2]/F");
|
| 211 |
rootTree_->Branch("clusterSize",clusterSize_,"clusterSize[2]/F");
|
| 212 |
rootTree_->Branch("clusterWidthX",clusterWidthX_,"clusterWidthX[2]/F");
|
| 213 |
rootTree_->Branch("clusterWidthY",clusterWidthY_,"clusterWidthY[2]/F");
|
| 214 |
rootTree_->Branch("clusterCharge",clusterCharge_,"clusterCharge[2]/i");
|
| 215 |
rootTree_->Branch("edge",edge_,"edge[2]/I");
|
| 216 |
rootTree_->Branch("momentum",&momentum_,"momentum/F");
|
| 217 |
rootTree_->Branch("run",&run_,"run/i");
|
| 218 |
rootTree_->Branch("event",&event_,"event/i");
|
| 219 |
|
| 220 |
}
|
| 221 |
|
| 222 |
|
| 223 |
HitRes::~HitRes()
|
| 224 |
{
|
| 225 |
|
| 226 |
// do anything here that needs to be done at desctruction time
|
| 227 |
// (e.g. close files, deallocate resources etc.)
|
| 228 |
|
| 229 |
cout << "Counters =" ;
|
| 230 |
cout << " Number of tracks: " << overlapCounts_[0]<<endl;
|
| 231 |
cout << " Number of valid hits: " << overlapCounts_[1]<<endl;
|
| 232 |
cout << " Number of overlaps: " << overlapCounts_[2]<<endl;
|
| 233 |
}
|
| 234 |
|
| 235 |
|
| 236 |
//
|
| 237 |
// member functions
|
| 238 |
//
|
| 239 |
|
| 240 |
|
| 241 |
// ------------ method called to for each event ------------
|
| 242 |
void
|
| 243 |
HitRes::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 244 |
{
|
| 245 |
using namespace edm;
|
| 246 |
//
|
| 247 |
// mag field & search tracker
|
| 248 |
//
|
| 249 |
edm::ESHandle<MagneticField> magFieldHandle;
|
| 250 |
iSetup.get<IdealMagneticFieldRecord>().get(magFieldHandle);
|
| 251 |
magField_ = magFieldHandle.product();
|
| 252 |
//
|
| 253 |
// propagator
|
| 254 |
//
|
| 255 |
AnalyticalPropagator propagator(magField_,anyDirection);
|
| 256 |
//
|
| 257 |
// geometry
|
| 258 |
//
|
| 259 |
edm::ESHandle<TrackerGeometry> geometryHandle;
|
| 260 |
iSetup.get<TrackerDigiGeometryRecord>().get(geometryHandle);
|
| 261 |
trackerGeometry_ = geometryHandle.product();
|
| 262 |
//
|
| 263 |
// make associator for SimHits
|
| 264 |
//
|
| 265 |
TrackerHitAssociator* associator;
|
| 266 |
if(doSimHit_) associator = new TrackerHitAssociator(iEvent, config_); else associator = 0;
|
| 267 |
|
| 268 |
//
|
| 269 |
// trajectories (from refit)
|
| 270 |
//
|
| 271 |
typedef vector<Trajectory> TrajectoryCollection;
|
| 272 |
edm::Handle<TrajectoryCollection> trajectoryCollectionHandle;
|
| 273 |
iEvent.getByLabel(trajectoryTag_,trajectoryCollectionHandle);
|
| 274 |
const TrajectoryCollection* trajectoryCollection = trajectoryCollectionHandle.product();
|
| 275 |
|
| 276 |
//
|
| 277 |
// loop over trajectories from refit
|
| 278 |
cout << "Tracks: " << trajectoryCollection->size()<<endl;
|
| 279 |
for ( TrajectoryCollection::const_iterator it=trajectoryCollection->begin();
|
| 280 |
it!=trajectoryCollection->end(); ++it ) analyze(*it,propagator,*associator);
|
| 281 |
|
| 282 |
run_ = iEvent.id().run();
|
| 283 |
event_ = iEvent.id().event();
|
| 284 |
|
| 285 |
}
|
| 286 |
|
| 287 |
|
| 288 |
void
|
| 289 |
HitRes::analyze (const Trajectory& trajectory,
|
| 290 |
const Propagator& propagator,
|
| 291 |
TrackerHitAssociator& associator)
|
| 292 |
{
|
| 293 |
typedef std::pair<const TrajectoryMeasurement*, const TrajectoryMeasurement*> Overlap;
|
| 294 |
typedef vector<Overlap> OverlapContainer;
|
| 295 |
++overlapCounts_[0];
|
| 296 |
|
| 297 |
OverlapContainer overlapHits;
|
| 298 |
|
| 299 |
// quality cuts on trajectory
|
| 300 |
// min. # hits / matched hits
|
| 301 |
|
| 302 |
if ( trajectory.foundHits()<6 ) return;
|
| 303 |
if ( ChiSquaredProbability((double)( trajectory.chiSquared() ),(double)( trajectory.ndof(false) )) < 0.001 ) return;
|
| 304 |
//
|
| 305 |
// loop over measurements in the trajectory and calculate residuals
|
| 306 |
//
|
| 307 |
|
| 308 |
vector<TrajectoryMeasurement> measurements(trajectory.measurements());
|
| 309 |
for ( vector<TrajectoryMeasurement>::const_iterator itm=measurements.begin();
|
| 310 |
itm!=measurements.end(); ++itm ) {
|
| 311 |
//
|
| 312 |
// skip "invalid" (i.e. missing) hits
|
| 313 |
//
|
| 314 |
ConstRecHitPointer hit = itm->recHit();
|
| 315 |
DetId id = hit->geographicalId();
|
| 316 |
int layer(layerFromId(id));
|
| 317 |
int subDet = id.subdetId();
|
| 318 |
|
| 319 |
if ( !hit->isValid() ) {
|
| 320 |
edm::LogVerbatim("HitRes") << "Invalid";
|
| 321 |
continue;
|
| 322 |
}
|
| 323 |
if (barrelOnly_ && (subDet==4 || subDet==6) ) return;
|
| 324 |
|
| 325 |
//edm::LogVerbatim("HitRes") << "Check " <<subDet << ", layer = " << layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(id).stereo()):2);
|
| 326 |
//cout << "Check SubID " <<subDet << ", layer = " << layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(id).stereo()):2) << endl;
|
| 327 |
|
| 328 |
//
|
| 329 |
// check for overlap: same subdet-id && layer number for
|
| 330 |
// two consecutive hits
|
| 331 |
//
|
| 332 |
++overlapCounts_[1];
|
| 333 |
if ( (layer!=-1 )&&(acceptLayer[subDet])) {
|
| 334 |
for (vector<TrajectoryMeasurement>::const_iterator itmCompare = itm-1;
|
| 335 |
itmCompare >= measurements.begin() && itmCompare > itm - 4;
|
| 336 |
--itmCompare){
|
| 337 |
|
| 338 |
DetId compareId = itmCompare->recHit()->geographicalId();
|
| 339 |
|
| 340 |
if ( subDet != compareId.subdetId() || layer != layerFromId(compareId)) break;
|
| 341 |
if (!itmCompare->recHit()->isValid()) continue;
|
| 342 |
if ( (subDet<=2) ||
|
| 343 |
(subDet > 2 && SiStripDetId(id).stereo()==SiStripDetId(compareId).stereo() )) {
|
| 344 |
overlapHits.push_back(std::make_pair(&(*itmCompare),&(*itm)));
|
| 345 |
//edm::LogVerbatim("HitRes") << "adding pair "<< ((subDet > 2)?(SiStripDetId(id).stereo()) : 2)
|
| 346 |
// << " from layer = " << layer;
|
| 347 |
//cout << "adding pair "<< ((subDet > 2)?(SiStripDetId(id).stereo()) : 2) << " from subDet = " << subDet << " and layer = " << layer;
|
| 348 |
//cout << " \t"<<run_<< "\t"<<event_<<"\t";
|
| 349 |
//cout << min(id.rawId(),compareId.rawId())<<"\t"<<max(id.rawId(),compareId.rawId())<<endl;
|
| 350 |
if ( SiStripDetId(id).glued() == id.rawId() ) cout << "BAD GLUED: Have glued layer with id = " << id.rawId() << " and glued id = " << SiStripDetId(id).glued() << " and stereo = " << SiStripDetId(id).stereo() << endl;
|
| 351 |
if ( SiStripDetId(compareId).glued() == compareId.rawId() ) cout << "BAD GLUED: Have glued layer with id = " << compareId.rawId() << " and glued id = " << SiStripDetId(compareId).glued() << " and stereo = " << SiStripDetId(compareId).stereo() << endl;
|
| 352 |
break;
|
| 353 |
}
|
| 354 |
}
|
| 355 |
}
|
| 356 |
}
|
| 357 |
//
|
| 358 |
// Loop over all overlap pairs.
|
| 359 |
//
|
| 360 |
hitCounts_[0] = trajectory.foundHits();
|
| 361 |
hitCounts_[1] = trajectory.lostHits();
|
| 362 |
chi2_[0] = trajectory.chiSquared();
|
| 363 |
chi2_[1] = trajectory.ndof(false);
|
| 364 |
|
| 365 |
for ( OverlapContainer::const_iterator iol=overlapHits.begin();
|
| 366 |
iol!=overlapHits.end(); ++iol ) {
|
| 367 |
//
|
| 368 |
// create reference state @ module 1 (no info from overlap hits)
|
| 369 |
//
|
| 370 |
++overlapCounts_[2];
|
| 371 |
// backward predicted state at module 1
|
| 372 |
TrajectoryStateOnSurface bwdPred1 = (*iol).first->backwardPredictedState();
|
| 373 |
if ( !bwdPred1.isValid() ) continue;
|
| 374 |
//cout << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl;
|
| 375 |
// forward predicted state at module 2
|
| 376 |
TrajectoryStateOnSurface fwdPred2 = (*iol).second->forwardPredictedState();
|
| 377 |
//cout << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl;
|
| 378 |
if ( !fwdPred2.isValid() ) continue;
|
| 379 |
// extrapolate fwdPred2 to module 1
|
| 380 |
TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2,bwdPred1.surface());
|
| 381 |
if ( !fwdPred2At1.isValid() ) continue;
|
| 382 |
// combine fwdPred2At1 with bwdPred1 (ref. state, best estimate without hits 1 and 2)
|
| 383 |
TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1,fwdPred2At1);
|
| 384 |
if ( !comb1.isValid() ) continue;
|
| 385 |
//
|
| 386 |
// propagation of reference parameters to module 2
|
| 387 |
//
|
| 388 |
std::pair<TrajectoryStateOnSurface,double> tsosWithS =
|
| 389 |
propagator.propagateWithPath(comb1,fwdPred2.surface());
|
| 390 |
TrajectoryStateOnSurface comb1At2 = tsosWithS.first;
|
| 391 |
if ( !comb1At2.isValid() ) continue;
|
| 392 |
//distance of propagation from one surface to the next==could cut here
|
| 393 |
overlapPath_ = tsosWithS.second;
|
| 394 |
if (abs(overlapPath_) > 15 ) continue; //cut to remove hit pairs > 15 cm apart
|
| 395 |
// global position on module 1
|
| 396 |
GlobalPoint position = comb1.globalPosition();
|
| 397 |
predictedPositions_[0][0] = position.x();
|
| 398 |
predictedPositions_[1][0] = position.y();
|
| 399 |
predictedPositions_[2][0] = position.z();
|
| 400 |
momentum_ = comb1.globalMomentum().mag();
|
| 401 |
//cout << "momentum from combination = " << momentum_ << endl;
|
| 402 |
//cout << "magnetic field from TSOS = " << comb1.magneticField()->inTesla(position).mag() << endl;
|
| 403 |
// local parameters and errors on module 1
|
| 404 |
AlgebraicVector5 pars = comb1.localParameters().vector();
|
| 405 |
AlgebraicSymMatrix55 errs = comb1.localError().matrix();
|
| 406 |
for ( int i=0; i<5; ++i ) {
|
| 407 |
predictedLocalParameters_[i][0] = pars[i];
|
| 408 |
predictedLocalErrors_[i][0] = sqrt(errs(i,i));
|
| 409 |
}
|
| 410 |
// global position on module 2
|
| 411 |
position = comb1At2.globalPosition();
|
| 412 |
predictedPositions_[0][1] = position.x();
|
| 413 |
predictedPositions_[1][1] = position.y();
|
| 414 |
predictedPositions_[2][1] = position.z();
|
| 415 |
// local parameters and errors on module 2
|
| 416 |
pars = comb1At2.localParameters().vector();
|
| 417 |
errs = comb1At2.localError().matrix();
|
| 418 |
for ( int i=0; i<5; ++i ) {
|
| 419 |
predictedLocalParameters_[i][1] = pars[i];
|
| 420 |
predictedLocalErrors_[i][1] = sqrt(errs(i,i));
|
| 421 |
}
|
| 422 |
|
| 423 |
//print out local errors in X to check
|
| 424 |
//cout << "Predicted local error in X at 1 = " << predictedLocalErrors_[3][0] << " and predicted local error in X at 2 is = " << predictedLocalErrors_[3][1] << endl;
|
| 425 |
//cout << "Predicted local error in Y at 1 = " << predictedLocalErrors_[4][0] << " and predicted local error in Y at 2 is = " << predictedLocalErrors_[4][1] << endl;
|
| 426 |
|
| 427 |
//
|
| 428 |
// jacobians (local-to-global@1,global 1-2,global-to-local@2)
|
| 429 |
//
|
| 430 |
JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(),
|
| 431 |
comb1.localParameters(),
|
| 432 |
*magField_);
|
| 433 |
AnalyticalCurvilinearJacobian jacCurvToCurv(comb1.globalParameters(),
|
| 434 |
comb1At2.globalPosition(),
|
| 435 |
comb1At2.globalMomentum(),
|
| 436 |
tsosWithS.second);
|
| 437 |
JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(),
|
| 438 |
comb1At2.localParameters(),
|
| 439 |
*magField_);
|
| 440 |
// combined jacobian local-1-to-local-2
|
| 441 |
AlgebraicMatrix55 jacobian =
|
| 442 |
jacLocToCurv.jacobian()*jacCurvToCurv.jacobian()*jacCurvToLoc.jacobian();
|
| 443 |
// covariance on module 1
|
| 444 |
AlgebraicSymMatrix55 covComb1 = comb1.localError().matrix();
|
| 445 |
// variance and correlations for predicted local_x on modules 1 and 2
|
| 446 |
double c00 = covComb1(3,3);
|
| 447 |
double c10(0.);
|
| 448 |
double c11(0.);
|
| 449 |
for ( int i=1; i<5; ++i ) {
|
| 450 |
c10 += jacobian(3,i)*covComb1(i,3);
|
| 451 |
for ( int j=1; j<5; ++j ) c11 += jacobian(3,i)*covComb1(i,j)*jacobian(3,j);
|
| 452 |
}
|
| 453 |
// choose relative sign in order to minimize error on difference
|
| 454 |
double diff = c00 - 2*fabs(c10) + c11;
|
| 455 |
diff = diff>0 ? sqrt(diff) : -sqrt(-diff);
|
| 456 |
predictedDeltaXError_ = diff;
|
| 457 |
relativeXSign_ = c10>0 ? -1 : 1;
|
| 458 |
//
|
| 459 |
// now find variance and correlations for predicted local_y
|
| 460 |
double c00Y = covComb1(4,4);
|
| 461 |
double c10Y(0.);
|
| 462 |
double c11Y(0.);
|
| 463 |
for ( int i=1; i<5; ++i ) {
|
| 464 |
c10Y += jacobian(4,i)*covComb1(i,4);
|
| 465 |
for ( int j=1; j<5; ++j ) c11Y += jacobian(4,i)*covComb1(i,j)*jacobian(4,j);
|
| 466 |
}
|
| 467 |
double diffY = c00Y - 2*fabs(c10Y) + c11Y;
|
| 468 |
diffY = diffY>0 ? sqrt(diffY) : -sqrt(-diffY);
|
| 469 |
predictedDeltaYError_ = diffY;
|
| 470 |
relativeYSign_ = c10Y>0 ? -1 : 1;
|
| 471 |
|
| 472 |
// information on modules and hits
|
| 473 |
overlapIds_[0] = (*iol).first->recHit()->geographicalId().rawId();
|
| 474 |
overlapIds_[1] = (*iol).second->recHit()->geographicalId().rawId();
|
| 475 |
|
| 476 |
if ( (*iol).first->recHit()->geographicalId().subdetId()==3 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId());
|
| 477 |
else if ( (*iol).first->recHit()->geographicalId().subdetId()==5 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId())+4;
|
| 478 |
else if ( (*iol).first->recHit()->geographicalId().subdetId()==4 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId())+10;
|
| 479 |
else if ( (*iol).first->recHit()->geographicalId().subdetId()==6 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId())+13;
|
| 480 |
else if ( (*iol).first->recHit()->geographicalId().subdetId()==1 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId())+20;
|
| 481 |
else if ( (*iol).first->recHit()->geographicalId().subdetId()==2 ) layer_ = layerFromId((*iol).first->recHit()->geographicalId().rawId())+30;
|
| 482 |
else layer_ = 99;
|
| 483 |
|
| 484 |
if ( overlapIds_[0] == SiStripDetId((*iol).first->recHit()->geographicalId()).glued() )
|
| 485 |
cout << "BAD GLUED: First Id = " << overlapIds_[0] << " has glued = " << SiStripDetId((*iol).first->recHit()->geographicalId()).glued() << " and stereo = " << SiStripDetId((*iol).first->recHit()->geographicalId()).stereo() << endl;
|
| 486 |
if ( overlapIds_[1] == SiStripDetId((*iol).second->recHit()->geographicalId()).glued() )
|
| 487 |
cout << "BAD GLUED: Second Id = " <<overlapIds_[1] << " has glued = " << SiStripDetId((*iol).second->recHit()->geographicalId()).glued() << " and stereo = " << SiStripDetId((*iol).second->recHit()->geographicalId()).stereo() << endl;
|
| 488 |
|
| 489 |
const TransientTrackingRecHit::ConstRecHitPointer firstRecHit = &(*(*iol).first->recHit());
|
| 490 |
const TransientTrackingRecHit::ConstRecHitPointer secondRecHit = &(*(*iol).second->recHit());
|
| 491 |
|
| 492 |
hitPositions_[0] = firstRecHit->localPosition().x();
|
| 493 |
hitErrors_[0] = sqrt(firstRecHit->localPositionError().xx());
|
| 494 |
hitPositions_[1] = secondRecHit->localPosition().x();
|
| 495 |
hitErrors_[1] = sqrt(secondRecHit->localPositionError().xx());
|
| 496 |
|
| 497 |
hitPositionsY_[0] = firstRecHit->localPosition().y();
|
| 498 |
hitErrorsY_[0] = sqrt(firstRecHit->localPositionError().yy());
|
| 499 |
hitPositionsY_[1] = secondRecHit->localPosition().y();
|
| 500 |
hitErrorsY_[1] = sqrt(secondRecHit->localPositionError().yy());
|
| 501 |
|
| 502 |
//cout << "printing local X hit position and error for the overlap hits. Hit 1 = " << hitPositions_[0] << "+-" << hitErrors_[0] << " and hit 2 is " << hitPositions_[1] << "+-" << hitErrors_[1] << endl;
|
| 503 |
|
| 504 |
DetId id1 = (*iol).first->recHit()->geographicalId();
|
| 505 |
DetId id2 = (*iol).second->recHit()->geographicalId();
|
| 506 |
int layer1 = layerFromId(id1);
|
| 507 |
int subDet1 = id1.subdetId();
|
| 508 |
int layer2 = layerFromId(id2);
|
| 509 |
int subDet2 = id2.subdetId();
|
| 510 |
if (abs(hitPositions_[0])>5) cout << "BAD: Bad hit position: Id = " << id1.rawId() << " stereo = " << SiStripDetId(id1).stereo() << " glued = " << SiStripDetId(id1).glued() << " from subdet = " << subDet1 << " and layer = " << layer1 << endl;
|
| 511 |
if (abs(hitPositions_[1])>5) cout << "BAD: Bad hit position: Id = " << id2.rawId() << " stereo = " << SiStripDetId(id2).stereo() << " glued = " << SiStripDetId(id2).glued() << " from subdet = " << subDet2 << " and layer = " << layer2 << endl;
|
| 512 |
|
| 513 |
// get track momentum
|
| 514 |
momentum_ = comb1.globalMomentum().mag();
|
| 515 |
|
| 516 |
// get cluster size
|
| 517 |
if (subDet1>2) { //strip
|
| 518 |
const TransientTrackingRecHit::ConstRecHitPointer thit1=(*iol).first->recHit();
|
| 519 |
const SiStripRecHit1D* hit1=dynamic_cast<const SiStripRecHit1D*>((*thit1).hit());
|
| 520 |
if (hit1) {
|
| 521 |
//check cluster width
|
| 522 |
const SiStripRecHit1D::ClusterRef & cluster1=hit1->cluster();
|
| 523 |
clusterSize_[0] = (cluster1->amplitudes()).size();
|
| 524 |
clusterWidthX_[0] = (cluster1->amplitudes()).size();
|
| 525 |
clusterWidthY_[0] = -1;
|
| 526 |
|
| 527 |
//check if cluster at edge of sensor
|
| 528 |
uint16_t firstStrip = cluster1->firstStrip();
|
| 529 |
uint16_t lastStrip = firstStrip + (cluster1->amplitudes()).size() -1;
|
| 530 |
unsigned short Nstrips;
|
| 531 |
Nstrips=reader->getNumberOfApvsAndStripLength(id1).first*128;
|
| 532 |
bool atEdge = false;
|
| 533 |
if (firstStrip == 0 || lastStrip == (Nstrips-1) ) atEdge = true;
|
| 534 |
if (atEdge) edge_[0] = 1; else edge_[0] = -1;
|
| 535 |
|
| 536 |
// get cluster total charge
|
| 537 |
const std::vector<uint8_t>& stripCharges = cluster1->amplitudes();
|
| 538 |
uint16_t charge = 0;
|
| 539 |
for (uint i = 0; i < stripCharges.size(); i++) {
|
| 540 |
charge += stripCharges.at(i);
|
| 541 |
}
|
| 542 |
clusterCharge_[0] = charge;
|
| 543 |
} else
|
| 544 |
cout << "Couldn't find SiStripRecHit1D first" << endl;
|
| 545 |
const TransientTrackingRecHit::ConstRecHitPointer thit2=(*iol).second->recHit();
|
| 546 |
const SiStripRecHit1D* hit2=dynamic_cast<const SiStripRecHit1D*>((*thit2).hit());
|
| 547 |
if (hit2) {
|
| 548 |
const SiStripRecHit1D::ClusterRef & cluster2=hit2->cluster();
|
| 549 |
clusterSize_[1] = (cluster2->amplitudes()).size();
|
| 550 |
clusterWidthX_[1] = (cluster2->amplitudes()).size();
|
| 551 |
clusterWidthY_[1] = -1;
|
| 552 |
|
| 553 |
uint16_t firstStrip = cluster2->firstStrip();
|
| 554 |
uint16_t lastStrip = firstStrip + (cluster2->amplitudes()).size() -1;
|
| 555 |
unsigned short Nstrips;
|
| 556 |
Nstrips=reader->getNumberOfApvsAndStripLength(id2).first*128;
|
| 557 |
bool atEdge = false;
|
| 558 |
if (firstStrip == 0 || lastStrip == (Nstrips-1) ) atEdge = true;
|
| 559 |
if (atEdge) edge_[1] = 1; else edge_[1] = -1;
|
| 560 |
|
| 561 |
// get cluster total charge
|
| 562 |
const std::vector<uint8_t>& stripCharges = cluster2->amplitudes();
|
| 563 |
uint16_t charge = 0;
|
| 564 |
for (uint i = 0; i < stripCharges.size(); i++) {
|
| 565 |
charge += stripCharges.at(i);
|
| 566 |
}
|
| 567 |
clusterCharge_[1] = charge;
|
| 568 |
|
| 569 |
} else
|
| 570 |
cout << "Couldn't find SiStripRecHit1D second" << endl;
|
| 571 |
//cout << "strip cluster size2 = " << clusterWidthX_[0] << " and size 2 = " << clusterWidthX_[1] << endl;
|
| 572 |
}
|
| 573 |
|
| 574 |
if (subDet2<3) { //pixel
|
| 575 |
|
| 576 |
const TransientTrackingRecHit::ConstRecHitPointer thit1=(*iol).first->recHit();
|
| 577 |
const SiPixelRecHit * recHitPix1 = dynamic_cast<const SiPixelRecHit *>((*thit1).hit());
|
| 578 |
if(recHitPix1) {
|
| 579 |
// check for cluster size and width
|
| 580 |
SiPixelRecHit::ClusterRef const& cluster1 = recHitPix1->cluster();
|
| 581 |
|
| 582 |
clusterSize_[0] = cluster1->size();
|
| 583 |
clusterWidthX_[0] = cluster1->sizeX();
|
| 584 |
clusterWidthY_[0] = cluster1->sizeY();
|
| 585 |
|
| 586 |
// check for cluster at edge
|
| 587 |
const PixelGeomDetUnit * theGeomDet =
|
| 588 |
dynamic_cast<const PixelGeomDetUnit*> ((*trackerGeometry_).idToDet(id1) );
|
| 589 |
const RectangularPixelTopology * topol =
|
| 590 |
dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
|
| 591 |
|
| 592 |
int minPixelRow = cluster1->minPixelRow(); //x
|
| 593 |
int maxPixelRow = cluster1->maxPixelRow();
|
| 594 |
int minPixelCol = cluster1->minPixelCol(); //y
|
| 595 |
int maxPixelCol = cluster1->maxPixelCol();
|
| 596 |
|
| 597 |
bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) ||
|
| 598 |
(topol->isItEdgePixelInX(maxPixelRow));
|
| 599 |
bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) ||
|
| 600 |
(topol->isItEdgePixelInY(maxPixelCol));
|
| 601 |
if (edgeHitX||edgeHitY) edge_[0] = 1; else edge_[0] = -1;
|
| 602 |
|
| 603 |
clusterCharge_[0] = (uint)cluster1->charge();
|
| 604 |
|
| 605 |
} else {
|
| 606 |
cout << "didn't find pixel cluster" << endl;
|
| 607 |
continue;
|
| 608 |
}
|
| 609 |
|
| 610 |
const TransientTrackingRecHit::ConstRecHitPointer thit2=(*iol).second->recHit();
|
| 611 |
const SiPixelRecHit * recHitPix2 = dynamic_cast<const SiPixelRecHit *>((*thit2).hit());
|
| 612 |
if(recHitPix2) {
|
| 613 |
SiPixelRecHit::ClusterRef const& cluster2 = recHitPix2->cluster();
|
| 614 |
|
| 615 |
clusterSize_[1] = cluster2->size();
|
| 616 |
clusterWidthX_[1] = cluster2->sizeX();
|
| 617 |
clusterWidthY_[1] = cluster2->sizeY();
|
| 618 |
//cout << "second pixel cluster is " << clusterSize_[1] << " pixels with x width = " << clusterWidthX_[1] << " and y width = " << clusterWidthY_[1] << endl;
|
| 619 |
|
| 620 |
const PixelGeomDetUnit * theGeomDet =
|
| 621 |
dynamic_cast<const PixelGeomDetUnit*> ((*trackerGeometry_).idToDet(id2) );
|
| 622 |
const RectangularPixelTopology * topol =
|
| 623 |
dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
|
| 624 |
|
| 625 |
int minPixelRow = cluster2->minPixelRow(); //x
|
| 626 |
int maxPixelRow = cluster2->maxPixelRow();
|
| 627 |
int minPixelCol = cluster2->minPixelCol(); //y
|
| 628 |
int maxPixelCol = cluster2->maxPixelCol();
|
| 629 |
|
| 630 |
bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) ||
|
| 631 |
(topol->isItEdgePixelInX(maxPixelRow));
|
| 632 |
bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) ||
|
| 633 |
(topol->isItEdgePixelInY(maxPixelCol));
|
| 634 |
if (edgeHitX||edgeHitY) edge_[1] = 1; else edge_[1] = -1;
|
| 635 |
|
| 636 |
clusterCharge_[1] = (uint)cluster2->charge();
|
| 637 |
|
| 638 |
} else {
|
| 639 |
cout << "didn't find pixel cluster" << endl;
|
| 640 |
continue;
|
| 641 |
}
|
| 642 |
|
| 643 |
}
|
| 644 |
|
| 645 |
|
| 646 |
//also check for edge pixels
|
| 647 |
|
| 648 |
//try writing out the SimHit info (for MC only)
|
| 649 |
if(doSimHit_){
|
| 650 |
std::vector<PSimHit> psimHits1;
|
| 651 |
std::vector<PSimHit> psimHits2;
|
| 652 |
//calculate layer
|
| 653 |
DetId id = (*iol).first->recHit()->geographicalId();
|
| 654 |
int layer(-1);
|
| 655 |
layer = layerFromId(id);
|
| 656 |
int subDet = id.subdetId();
|
| 657 |
edm::LogVerbatim("HitRes") << "Subdet = " << subDet << " ; layer = " << layer;
|
| 658 |
|
| 659 |
psimHits1 = associator.associateHit( *(firstRecHit->hit()) );
|
| 660 |
edm::LogVerbatim("HitRes") << "single hit ";
|
| 661 |
edm::LogVerbatim("HitRes") << "length of psimHits1: " << psimHits1.size();
|
| 662 |
if ( !psimHits1.empty() ) {
|
| 663 |
float closest_dist = 99999.9;
|
| 664 |
std::vector<PSimHit>::const_iterator closest_simhit = psimHits1.begin();
|
| 665 |
for (std::vector<PSimHit>::const_iterator m = psimHits1.begin(); m < psimHits1.end(); m++) {
|
| 666 |
//find closest simHit to the recHit
|
| 667 |
float simX = (*m).localPosition().x();
|
| 668 |
float dist = fabs( simX - ((*iol).first->recHit()->localPosition().x()) );
|
| 669 |
edm::LogVerbatim("HitRes") << "simHit1 simX = " << simX << " hitX = " << (*iol).first->recHit()->localPosition().x() << " distX = " << dist << " layer = " << layer;
|
| 670 |
if (dist<closest_dist) {
|
| 671 |
//cout << "found newest closest dist for simhit1" << endl;
|
| 672 |
closest_dist = dist;
|
| 673 |
closest_simhit = m;
|
| 674 |
}
|
| 675 |
}
|
| 676 |
//if glued layer, convert sim hit position to matchedhit surface
|
| 677 |
//layer index from 1-4 for TIB, 1-6 for TOB
|
| 678 |
// Are the sim hits on the glued layers or are they split???
|
| 679 |
if ( subDet > 2 && !SiStripDetId(id).glued() ) {
|
| 680 |
const GluedGeomDet* gluedDet = (const GluedGeomDet*)(*trackerGeometry_).idToDet((*firstRecHit).hit()->geographicalId());
|
| 681 |
const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) gluedDet->monoDet();
|
| 682 |
GlobalPoint gp = stripDet->surface().toGlobal( (*closest_simhit).localPosition() );
|
| 683 |
LocalPoint lp = gluedDet->surface().toLocal( gp );
|
| 684 |
LocalVector localdirection = (*closest_simhit).localDirection();
|
| 685 |
GlobalVector globaldirection = stripDet->surface().toGlobal(localdirection);
|
| 686 |
LocalVector direction = gluedDet->surface().toLocal(globaldirection);
|
| 687 |
float scale = -lp.z() / direction.z();
|
| 688 |
LocalPoint projectedPos = lp + scale*direction;
|
| 689 |
simHitPositions_[0] = projectedPos.x();
|
| 690 |
edm::LogVerbatim("HitRes") << "simhit position from matched layer = " << simHitPositions_[0];
|
| 691 |
simHitPositionsY_[0] = projectedPos.y();
|
| 692 |
} else {
|
| 693 |
simHitPositions_[0] = (*closest_simhit).localPosition().x();
|
| 694 |
simHitPositionsY_[0] = (*closest_simhit).localPosition().y();
|
| 695 |
edm::LogVerbatim("HitRes") << "simhit position from non-matched layer = " << simHitPositions_[0];
|
| 696 |
}
|
| 697 |
edm::LogVerbatim("HitRes") << "hit position = " << hitPositions_[0];
|
| 698 |
} else {
|
| 699 |
simHitPositions_[0] = -99.;
|
| 700 |
simHitPositionsY_[0] = -99.;
|
| 701 |
//cout << " filling simHitX: " << -99 << endl;
|
| 702 |
}
|
| 703 |
|
| 704 |
psimHits2 = associator.associateHit( *(secondRecHit->hit()) );
|
| 705 |
if ( !psimHits2.empty() ) {
|
| 706 |
float closest_dist = 99999.9;
|
| 707 |
std::vector<PSimHit>::const_iterator closest_simhit = psimHits2.begin();
|
| 708 |
for (std::vector<PSimHit>::const_iterator m = psimHits2.begin(); m < psimHits2.end(); m++) {
|
| 709 |
float simX = (*m).localPosition().x();
|
| 710 |
float dist = fabs( simX - ((*iol).second->recHit()->localPosition().x()) );
|
| 711 |
if (dist<closest_dist) {
|
| 712 |
closest_dist = dist;
|
| 713 |
closest_simhit = m;
|
| 714 |
}
|
| 715 |
}
|
| 716 |
//if glued layer, convert sim hit position to matchedhit surface
|
| 717 |
// if no sim hits on matched layers then this section can be removed
|
| 718 |
if ( subDet > 2 && !SiStripDetId(id).glued() ) {
|
| 719 |
const GluedGeomDet* gluedDet = (const GluedGeomDet*)(*trackerGeometry_).idToDet((*secondRecHit).hit()->geographicalId());
|
| 720 |
const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) gluedDet->monoDet();
|
| 721 |
GlobalPoint gp = stripDet->surface().toGlobal( (*closest_simhit).localPosition() );
|
| 722 |
LocalPoint lp = gluedDet->surface().toLocal( gp );
|
| 723 |
LocalVector localdirection = (*closest_simhit).localDirection();
|
| 724 |
GlobalVector globaldirection = stripDet->surface().toGlobal(localdirection);
|
| 725 |
LocalVector direction = gluedDet->surface().toLocal(globaldirection);
|
| 726 |
float scale = -lp.z() / direction.z();
|
| 727 |
LocalPoint projectedPos = lp + scale*direction;
|
| 728 |
simHitPositions_[1] = projectedPos.x();
|
| 729 |
simHitPositionsY_[1] = projectedPos.y();
|
| 730 |
} else {
|
| 731 |
simHitPositions_[1] = (*closest_simhit).localPosition().x();
|
| 732 |
simHitPositionsY_[1] = (*closest_simhit).localPosition().y();
|
| 733 |
}
|
| 734 |
} else {
|
| 735 |
simHitPositions_[1] = -99.;
|
| 736 |
simHitPositionsY_[1] = -99.;
|
| 737 |
}
|
| 738 |
}
|
| 739 |
rootTree_->Fill();
|
| 740 |
}
|
| 741 |
}
|
| 742 |
|
| 743 |
int
|
| 744 |
HitRes::layerFromId (const DetId& id) const
|
| 745 |
{
|
| 746 |
if ( id.subdetId()==PixelSubdetector::PixelBarrel ) {
|
| 747 |
PXBDetId tobId(id);
|
| 748 |
return tobId.layer();
|
| 749 |
}
|
| 750 |
else if ( id.subdetId()==PixelSubdetector::PixelEndcap ) {
|
| 751 |
PXFDetId tobId(id);
|
| 752 |
return tobId.disk() + (3*(tobId.side()-1));
|
| 753 |
}
|
| 754 |
else if ( id.subdetId()==StripSubdetector::TIB ) {
|
| 755 |
TIBDetId tibId(id);
|
| 756 |
return tibId.layer();
|
| 757 |
}
|
| 758 |
else if ( id.subdetId()==StripSubdetector::TOB ) {
|
| 759 |
TOBDetId tobId(id);
|
| 760 |
return tobId.layer();
|
| 761 |
}
|
| 762 |
else if ( id.subdetId()==StripSubdetector::TEC ) {
|
| 763 |
TECDetId tobId(id);
|
| 764 |
return tobId.wheel() + (9*(tobId.side()-1));
|
| 765 |
}
|
| 766 |
else if ( id.subdetId()==StripSubdetector::TID ) {
|
| 767 |
TIDDetId tobId(id);
|
| 768 |
return tobId.wheel() + (3*(tobId.side()-1));
|
| 769 |
}
|
| 770 |
return -1;
|
| 771 |
}
|
| 772 |
|
| 773 |
void
|
| 774 |
HitRes::endJob() {
|
| 775 |
if ( rootTree_ ) {
|
| 776 |
rootTree_->GetDirectory()->cd();
|
| 777 |
rootTree_->Write();
|
| 778 |
delete rootTree_;
|
| 779 |
}
|
| 780 |
}
|
| 781 |
|
| 782 |
//define this as a plug-in
|
| 783 |
DEFINE_FWK_MODULE(HitRes);
|
