| 1 |
msegala |
1.1 |
////////////////////////////////////////////////////////////////////////////////
|
| 2 |
|
|
// Package: TrackingAnalysis/Cosmics
|
| 3 |
|
|
// Class: HitEff
|
| 4 |
|
|
// Original Author: Daniele Benedetti-INFN perugia
|
| 5 |
|
|
// Update for P5: Keith Ulmer--University of Colorado
|
| 6 |
|
|
//
|
| 7 |
|
|
///////////////////////////////////////////////////////////////////////////////
|
| 8 |
|
|
|
| 9 |
|
|
// system include files
|
| 10 |
|
|
#include <memory>
|
| 11 |
|
|
#include <string>
|
| 12 |
|
|
#include <iostream>
|
| 13 |
|
|
|
| 14 |
|
|
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
| 15 |
|
|
#include "FWCore/Framework/interface/EDAnalyzer.h"
|
| 16 |
|
|
#include "FWCore/Framework/interface/Event.h"
|
| 17 |
|
|
#include "FWCore/Framework/interface/MakerMacros.h"
|
| 18 |
|
|
|
| 19 |
|
|
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
| 20 |
|
|
#include "TrackingAnalysis/Cosmics/interface/HitEff.h"
|
| 21 |
|
|
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
|
| 22 |
|
|
#include "DataFormats/Common/interface/Handle.h"
|
| 23 |
|
|
#include "FWCore/Framework/interface/ESHandle.h"
|
| 24 |
|
|
#include "FWCore/Framework/interface/EventSetup.h"
|
| 25 |
|
|
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
|
| 26 |
|
|
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
|
| 27 |
|
|
#include "DataFormats/GeometryVector/interface/LocalVector.h"
|
| 28 |
|
|
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
|
| 29 |
|
|
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
|
| 30 |
|
|
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
|
| 31 |
|
|
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
|
| 32 |
|
|
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
|
| 33 |
|
|
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
|
| 34 |
|
|
#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
|
| 35 |
|
|
#include "DataFormats/TrackReco/interface/Track.h"
|
| 36 |
|
|
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
| 37 |
|
|
#include "DataFormats/TrackReco/interface/TrackExtra.h"
|
| 38 |
|
|
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
| 39 |
|
|
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
|
| 40 |
|
|
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
|
| 41 |
|
|
#include "TrackingAnalysis/Cosmics/interface/TrajectoryAtValidHit.h"
|
| 42 |
|
|
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
|
| 43 |
|
|
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
|
| 44 |
|
|
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
|
| 45 |
|
|
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
|
| 46 |
|
|
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
|
| 47 |
|
|
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
|
| 48 |
|
|
#include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
|
| 49 |
|
|
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
|
| 50 |
|
|
#include "DataFormats/TrackReco/interface/DeDxData.h"
|
| 51 |
|
|
|
| 52 |
|
|
#include "AnalysisDataFormats/SiStripClusterInfo/interface/SiStripClusterInfo.h"
|
| 53 |
|
|
#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
|
| 54 |
|
|
#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
|
| 55 |
|
|
#include "CalibTracker/Records/interface/SiStripQualityRcd.h"
|
| 56 |
|
|
#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
|
| 57 |
|
|
#include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
|
| 58 |
|
|
#include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
|
| 59 |
|
|
#include "DataFormats/Common/interface/DetSetVector.h"
|
| 60 |
|
|
#include "DataFormats/Common/interface/DetSetVectorNew.h"
|
| 61 |
|
|
#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
|
| 62 |
|
|
#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
|
| 63 |
|
|
|
| 64 |
|
|
#include "DataFormats/MuonReco/interface/Muon.h"
|
| 65 |
|
|
#include "DataFormats/MuonReco/interface/MuonFwd.h"
|
| 66 |
|
|
|
| 67 |
|
|
#include "CommonTools/UtilAlgos/interface/TFileService.h"
|
| 68 |
|
|
|
| 69 |
|
|
|
| 70 |
|
|
#include "TMath.h"
|
| 71 |
|
|
#include "TH1F.h"
|
| 72 |
|
|
|
| 73 |
|
|
//
|
| 74 |
|
|
// constructors and destructor
|
| 75 |
|
|
//
|
| 76 |
|
|
|
| 77 |
|
|
using namespace std;
|
| 78 |
|
|
HitEff::HitEff(const edm::ParameterSet& conf) :
|
| 79 |
|
|
conf_(conf)
|
| 80 |
|
|
{
|
| 81 |
|
|
layers =conf_.getParameter<int>("Layer");
|
| 82 |
|
|
DEBUG = conf_.getParameter<bool>("Debug");
|
| 83 |
|
|
|
| 84 |
|
|
edm::Service<TFileService> fs;
|
| 85 |
|
|
|
| 86 |
|
|
traj = fs->make<TTree>("traj","tree of trajectory positions");
|
| 87 |
|
|
traj->Branch("TrajGlbX",&TrajGlbX,"TrajGlbX/F");
|
| 88 |
|
|
traj->Branch("TrajGlbY",&TrajGlbY,"TrajGlbY/F");
|
| 89 |
|
|
traj->Branch("TrajGlbZ",&TrajGlbZ,"TrajGlbZ/F");
|
| 90 |
|
|
traj->Branch("TrajLocX",&TrajLocX,"TrajLocX/F");
|
| 91 |
|
|
traj->Branch("TrajLocY",&TrajLocY,"TrajLocY/F");
|
| 92 |
|
|
traj->Branch("TrajLocErrX",&TrajLocErrX,"TrajLocErrX/F");
|
| 93 |
|
|
traj->Branch("TrajLocErrY",&TrajLocErrY,"TrajLocErrY/F");
|
| 94 |
|
|
traj->Branch("TrajLocAngleX",&TrajLocAngleX,"TrajLocAngleX/F");
|
| 95 |
|
|
traj->Branch("TrajLocAngleY",&TrajLocAngleY,"TrajLocAngleY/F");
|
| 96 |
|
|
traj->Branch("ClusterLocX",&ClusterLocX,"ClusterLocX/F");
|
| 97 |
|
|
traj->Branch("ClusterLocY",&ClusterLocY,"ClusterLocY/F");
|
| 98 |
|
|
traj->Branch("ClusterLocErrX",&ClusterLocErrX,"ClusterLocErrX/F");
|
| 99 |
|
|
traj->Branch("ClusterLocErrY",&ClusterLocErrY,"ClusterLocErrY/F");
|
| 100 |
|
|
traj->Branch("ClusterStoN",&ClusterStoN,"ClusterStoN/F");
|
| 101 |
|
|
traj->Branch("ResX",&ResX,"ResX/F");
|
| 102 |
|
|
traj->Branch("ResXSig",&ResXSig,"ResXSig/F");
|
| 103 |
|
|
traj->Branch("ModIsBad",&ModIsBad,"ModIsBad/i");
|
| 104 |
|
|
traj->Branch("SiStripQualBad",&SiStripQualBad,"SiStripQualBad/i");
|
| 105 |
|
|
traj->Branch("withinAcceptance",&withinAcceptance, "withinAcceptance/O");
|
| 106 |
|
|
traj->Branch("Id",&Id,"Id/i");
|
| 107 |
|
|
traj->Branch("run",&run,"run/i");
|
| 108 |
|
|
traj->Branch("event",&event,"event/i");
|
| 109 |
|
|
traj->Branch("layer",&whatlayer,"layer/i");
|
| 110 |
|
|
traj->Branch("timeDT",&timeDT,"timeDT/F");
|
| 111 |
|
|
traj->Branch("timeDTErr",&timeDTErr,"timeDTErr/F");
|
| 112 |
|
|
traj->Branch("timeDTDOF",&timeDTDOF,"timeDTDOF/I");
|
| 113 |
|
|
traj->Branch("timeECAL",&timeECAL,"timeECAL/F");
|
| 114 |
|
|
traj->Branch("dedx",&dedx,"dedx/F");
|
| 115 |
|
|
traj->Branch("dedxNOM",&dedxNOM,"dedxNOM/I");
|
| 116 |
|
|
|
| 117 |
|
|
events = 0;
|
| 118 |
|
|
EventTrackCKF = 0;
|
| 119 |
|
|
|
| 120 |
|
|
|
| 121 |
|
|
}
|
| 122 |
|
|
|
| 123 |
|
|
// Virtual destructor needed.
|
| 124 |
|
|
HitEff::~HitEff() { }
|
| 125 |
|
|
/*
|
| 126 |
|
|
void HitEff::beginJob(const edm::EventSetup& c){
|
| 127 |
|
|
|
| 128 |
|
|
edm::Service<TFileService> fs;
|
| 129 |
|
|
|
| 130 |
|
|
traj = fs->make<TTree>("traj","tree of trajectory positions");
|
| 131 |
|
|
traj->Branch("TrajGlbX",&TrajGlbX,"TrajGlbX/F");
|
| 132 |
|
|
traj->Branch("TrajGlbY",&TrajGlbY,"TrajGlbY/F");
|
| 133 |
|
|
traj->Branch("TrajGlbZ",&TrajGlbZ,"TrajGlbZ/F");
|
| 134 |
|
|
traj->Branch("TrajLocX",&TrajLocX,"TrajLocX/F");
|
| 135 |
|
|
traj->Branch("TrajLocY",&TrajLocY,"TrajLocY/F");
|
| 136 |
|
|
traj->Branch("TrajLocErrX",&TrajLocErrX,"TrajLocErrX/F");
|
| 137 |
|
|
traj->Branch("TrajLocErrY",&TrajLocErrY,"TrajLocErrY/F");
|
| 138 |
|
|
traj->Branch("TrajLocAngleX",&TrajLocAngleX,"TrajLocAngleX/F");
|
| 139 |
|
|
traj->Branch("TrajLocAngleY",&TrajLocAngleY,"TrajLocAngleY/F");
|
| 140 |
|
|
traj->Branch("ClusterLocX",&ClusterLocX,"ClusterLocX/F");
|
| 141 |
|
|
traj->Branch("ClusterLocY",&ClusterLocY,"ClusterLocY/F");
|
| 142 |
|
|
traj->Branch("ClusterLocErrX",&ClusterLocErrX,"ClusterLocErrX/F");
|
| 143 |
|
|
traj->Branch("ClusterLocErrY",&ClusterLocErrY,"ClusterLocErrY/F");
|
| 144 |
|
|
traj->Branch("ClusterStoN",&ClusterStoN,"ClusterStoN/F");
|
| 145 |
|
|
traj->Branch("ResX",&ResX,"ResX/F");
|
| 146 |
|
|
traj->Branch("ResXSig",&ResXSig,"ResXSig/F");
|
| 147 |
|
|
traj->Branch("ModIsBad",&ModIsBad,"ModIsBad/i");
|
| 148 |
|
|
traj->Branch("SiStripQualBad",&SiStripQualBad,"SiStripQualBad/i");
|
| 149 |
|
|
traj->Branch("withinAcceptance",&withinAcceptance, "withinAcceptance/O");
|
| 150 |
|
|
traj->Branch("Id",&Id,"Id/i");
|
| 151 |
|
|
traj->Branch("run",&run,"run/i");
|
| 152 |
|
|
traj->Branch("event",&event,"event/i");
|
| 153 |
|
|
traj->Branch("layer",&whatlayer,"layer/i");
|
| 154 |
|
|
traj->Branch("timeDT",&timeDT,"timeDT/F");
|
| 155 |
|
|
traj->Branch("timeDTErr",&timeDTErr,"timeDTErr/F");
|
| 156 |
|
|
traj->Branch("timeDTDOF",&timeDTDOF,"timeDTDOF/I");
|
| 157 |
|
|
traj->Branch("timeECAL",&timeECAL,"timeECAL/F");
|
| 158 |
|
|
traj->Branch("dedx",&dedx,"dedx/F");
|
| 159 |
|
|
traj->Branch("dedxNOM",&dedxNOM,"dedxNOM/I");
|
| 160 |
|
|
|
| 161 |
|
|
events = 0;
|
| 162 |
|
|
EventTrackCKF = 0;
|
| 163 |
|
|
|
| 164 |
|
|
}
|
| 165 |
|
|
*/
|
| 166 |
|
|
|
| 167 |
|
|
void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es){
|
| 168 |
|
|
|
| 169 |
|
|
// bool DEBUG = false;
|
| 170 |
|
|
|
| 171 |
|
|
if (DEBUG) cout << "beginning analyze from HitEff" << endl;
|
| 172 |
|
|
|
| 173 |
|
|
using namespace edm;
|
| 174 |
|
|
using namespace reco;
|
| 175 |
|
|
// Step A: Get Inputs
|
| 176 |
|
|
|
| 177 |
|
|
int run_nr = e.id().run();
|
| 178 |
|
|
int ev_nr = e.id().event();
|
| 179 |
|
|
|
| 180 |
|
|
//CombinatoriaTrack
|
| 181 |
|
|
edm::Handle<reco::TrackCollection> trackCollectionCKF;
|
| 182 |
|
|
edm::InputTag TkTagCKF = conf_.getParameter<edm::InputTag>("combinatorialTracks");
|
| 183 |
|
|
e.getByLabel(TkTagCKF,trackCollectionCKF);
|
| 184 |
|
|
|
| 185 |
|
|
edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
|
| 186 |
|
|
edm::InputTag TkTrajCKF = conf_.getParameter<edm::InputTag>("trajectories");
|
| 187 |
|
|
e.getByLabel(TkTrajCKF,TrajectoryCollectionCKF);
|
| 188 |
|
|
|
| 189 |
|
|
// Clusters
|
| 190 |
|
|
// get the SiStripClusters from the event
|
| 191 |
|
|
edm::Handle< edmNew::DetSetVector<SiStripCluster> > theStripClusters;
|
| 192 |
|
|
e.getByLabel("siStripClusters", theStripClusters);
|
| 193 |
|
|
|
| 194 |
|
|
//should try getting the pixel clusters here, too.
|
| 195 |
|
|
edm::Handle<edmNew::DetSetVector<SiPixelCluster> > thePixelClusters;
|
| 196 |
|
|
e.getByLabel("siPixelClusters", thePixelClusters);
|
| 197 |
|
|
|
| 198 |
|
|
//get tracker geometry
|
| 199 |
|
|
edm::ESHandle<TrackerGeometry> tracker;
|
| 200 |
|
|
es.get<TrackerDigiGeometryRecord>().get(tracker);
|
| 201 |
|
|
const TrackerGeometry * tkgeom=&(* tracker);
|
| 202 |
|
|
|
| 203 |
|
|
//get strip Cluster Parameter Estimator
|
| 204 |
|
|
edm::ESHandle<StripClusterParameterEstimator> stripCPE;
|
| 205 |
|
|
es.get<TkStripCPERecord>().get("StripCPEfromTrackAngle", stripCPE);
|
| 206 |
|
|
const StripClusterParameterEstimator &stripcpe(*stripCPE);
|
| 207 |
|
|
|
| 208 |
|
|
//get pixel Cluster Parameter Estimator
|
| 209 |
|
|
edm::ESHandle<PixelClusterParameterEstimator> pixelCPE;
|
| 210 |
|
|
//es.get<TkPixelCPERecord>().get("PixelCPEfromTrackAngle", pixelCPE);
|
| 211 |
|
|
es.get<TkPixelCPERecord>().get("PixelCPETemplateReco", pixelCPE);
|
| 212 |
|
|
const PixelClusterParameterEstimator &pixelcpe(*pixelCPE);
|
| 213 |
|
|
|
| 214 |
|
|
// get the SiStripQuality records
|
| 215 |
|
|
edm::ESHandle<SiStripQuality> SiStripQuality_;
|
| 216 |
|
|
try {
|
| 217 |
|
|
es.get<SiStripQualityRcd>().get("forCluster",SiStripQuality_);
|
| 218 |
|
|
}
|
| 219 |
|
|
catch (...) {
|
| 220 |
|
|
es.get<SiStripQualityRcd>().get(SiStripQuality_);
|
| 221 |
|
|
}
|
| 222 |
|
|
|
| 223 |
|
|
edm::ESHandle<MagneticField> magFieldHandle;
|
| 224 |
|
|
es.get<IdealMagneticFieldRecord>().get(magFieldHandle);
|
| 225 |
|
|
const MagneticField* magField_ = magFieldHandle.product();
|
| 226 |
|
|
|
| 227 |
|
|
events++;
|
| 228 |
|
|
|
| 229 |
|
|
// *************** SiStripCluster Collection
|
| 230 |
|
|
const edmNew::DetSetVector<SiStripCluster>& stripClusters = *theStripClusters;
|
| 231 |
|
|
|
| 232 |
|
|
//go through clusters to write out global position of good clusters for the layer understudy for comparison
|
| 233 |
|
|
// Loop through clusters just to print out locations
|
| 234 |
|
|
|
| 235 |
|
|
for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = stripClusters.begin(); DSViter != stripClusters.end(); DSViter++) {
|
| 236 |
|
|
// DSViter is a vector of SiStripClusters located on a single module
|
| 237 |
|
|
unsigned int ClusterId = DSViter->id();
|
| 238 |
|
|
DetId ClusterDetId(ClusterId);
|
| 239 |
|
|
const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
|
| 240 |
|
|
|
| 241 |
|
|
edmNew::DetSet<SiStripCluster>::const_iterator begin=DSViter->begin();
|
| 242 |
|
|
edmNew::DetSet<SiStripCluster>::const_iterator end =DSViter->end();
|
| 243 |
|
|
for(edmNew::DetSet<SiStripCluster>::const_iterator iter=begin;iter!=end;++iter) {
|
| 244 |
|
|
//iter is a single SiStripCluster
|
| 245 |
|
|
StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
|
| 246 |
|
|
|
| 247 |
|
|
const Surface* surface;
|
| 248 |
|
|
surface = &(tracker->idToDet(ClusterDetId)->surface());
|
| 249 |
|
|
LocalPoint lp = parameters.first;
|
| 250 |
|
|
GlobalPoint gp = surface->toGlobal(lp);
|
| 251 |
|
|
|
| 252 |
|
|
uint layer=0;
|
| 253 |
|
|
unsigned int subid=ClusterDetId.subdetId();
|
| 254 |
|
|
|
| 255 |
|
|
if (subid== StripSubdetector::TIB) {
|
| 256 |
|
|
layer = TIBDetId(ClusterDetId).layer();
|
| 257 |
|
|
}
|
| 258 |
|
|
if (subid== StripSubdetector::TOB) {
|
| 259 |
|
|
layer = TOBDetId(ClusterDetId).layer() + 4;
|
| 260 |
|
|
}
|
| 261 |
|
|
if (subid== StripSubdetector::TID) {
|
| 262 |
|
|
layer = TIDDetId(ClusterDetId).wheel() + 10;
|
| 263 |
|
|
}
|
| 264 |
|
|
if (subid== StripSubdetector::TEC) {
|
| 265 |
|
|
layer = TECDetId(ClusterDetId).wheel() + 13;
|
| 266 |
|
|
}
|
| 267 |
|
|
if ( subid==PixelSubdetector::PixelBarrel ) {
|
| 268 |
|
|
layer = PXBDetId(ClusterDetId).layer() + 22;
|
| 269 |
|
|
}
|
| 270 |
|
|
if ( subid==PixelSubdetector::PixelEndcap ) {
|
| 271 |
|
|
layer = PXFDetId(ClusterDetId).disk() + 25;
|
| 272 |
|
|
}
|
| 273 |
|
|
|
| 274 |
|
|
if(DEBUG) cout << "Found hit in cluster collection layer = " << layer << " with id = " << ClusterId << " local X position = " << lp.x() << " +- " << sqrt(parameters.second.xx()) << " matched/stereo/rphi = " << ((ClusterId & 0x3)==0) << "/" << ((ClusterId & 0x3)==1) << "/" << ((ClusterId & 0x3)==2) << endl;
|
| 275 |
|
|
}
|
| 276 |
|
|
}
|
| 277 |
|
|
|
| 278 |
|
|
// *************** SiPixelCluster Collection
|
| 279 |
|
|
const edmNew::DetSetVector<SiPixelCluster>& pixelClusters = *thePixelClusters;
|
| 280 |
|
|
|
| 281 |
|
|
//go through clusters to write out global position of good clusters for the layer understudy for comparison
|
| 282 |
|
|
// Loop through clusters just to print out locations
|
| 283 |
|
|
|
| 284 |
|
|
for (edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = pixelClusters.begin(); DSViter != pixelClusters.end(); DSViter++) {
|
| 285 |
|
|
// DSViter is a vector of SiStripClusters located on a single module
|
| 286 |
|
|
unsigned int pixelClusterId = DSViter->id();
|
| 287 |
|
|
|
| 288 |
|
|
cout << "pixel cluster ID = " << pixelClusterId << endl;
|
| 289 |
|
|
|
| 290 |
|
|
}
|
| 291 |
|
|
|
| 292 |
|
|
// Tracking
|
| 293 |
|
|
const reco::TrackCollection *tracksCKF=trackCollectionCKF.product();
|
| 294 |
|
|
if (DEBUG) cout << "number ckf tracks found = " << tracksCKF->size() << endl;
|
| 295 |
|
|
if (tracksCKF->size() == 1 ){
|
| 296 |
|
|
if (DEBUG) cout << "starting checking good single track event" << endl;
|
| 297 |
|
|
reco::TrackCollection::const_iterator iCKF=trackCollectionCKF.product()->begin();
|
| 298 |
|
|
EventTrackCKF++;
|
| 299 |
|
|
|
| 300 |
|
|
//get dEdx info if available
|
| 301 |
|
|
Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
|
| 302 |
|
|
if (e.getByLabel("dedxMedianCTF", dEdxUncalibHandle)) {
|
| 303 |
|
|
const ValueMap<DeDxData> dEdxTrackUncalib = *dEdxUncalibHandle.product();
|
| 304 |
|
|
|
| 305 |
|
|
reco::TrackRef itTrack = reco::TrackRef( trackCollectionCKF, 0 );
|
| 306 |
|
|
dedx = dEdxTrackUncalib[itTrack].dEdx();
|
| 307 |
|
|
dedxNOM = dEdxTrackUncalib[itTrack].numberOfMeasurements();
|
| 308 |
|
|
} else {
|
| 309 |
|
|
dedx = -999.0; dedxNOM = -999;
|
| 310 |
|
|
}
|
| 311 |
|
|
|
| 312 |
|
|
//get muon and ecal timing info if available
|
| 313 |
|
|
Handle<MuonCollection> muH;
|
| 314 |
|
|
if(e.getByLabel("muonsWitht0Correction",muH)){
|
| 315 |
|
|
const MuonCollection & muonsT0 = *muH.product();
|
| 316 |
|
|
if(muonsT0.size()!=0) {
|
| 317 |
|
|
MuonTime mt0 = muonsT0[0].time();
|
| 318 |
|
|
timeDT = mt0.timeAtIpInOut;
|
| 319 |
|
|
timeDTErr = mt0.timeAtIpInOutErr;
|
| 320 |
|
|
timeDTDOF = mt0.nDof;
|
| 321 |
|
|
|
| 322 |
|
|
bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
|
| 323 |
|
|
if (hasCaloEnergyInfo) timeECAL = muonsT0[0].calEnergy().ecal_time;
|
| 324 |
|
|
}
|
| 325 |
|
|
} else {
|
| 326 |
|
|
timeDT = -999.0; timeDTErr = -999.0; timeDTDOF = -999; timeECAL = -999.0;
|
| 327 |
|
|
}
|
| 328 |
|
|
|
| 329 |
|
|
const Trajectory traject = *(TrajectoryCollectionCKF.product()->begin());
|
| 330 |
|
|
std::vector<TrajectoryMeasurement> TMeas=traject.measurements();
|
| 331 |
|
|
vector<TrajectoryMeasurement>::iterator itm;
|
| 332 |
|
|
double xloc = 0.;
|
| 333 |
|
|
double yloc = 0.;
|
| 334 |
|
|
double xErr = 0.;
|
| 335 |
|
|
double yErr = 0.;
|
| 336 |
|
|
double angleX = -999.;
|
| 337 |
|
|
double angleY = -999.;
|
| 338 |
|
|
double xglob,yglob,zglob;
|
| 339 |
|
|
|
| 340 |
|
|
for (itm=TMeas.begin();itm!=TMeas.end();itm++){
|
| 341 |
|
|
ConstReferenceCountingPointer<TransientTrackingRecHit> theInHit;
|
| 342 |
|
|
theInHit = (*itm).recHit();
|
| 343 |
|
|
|
| 344 |
|
|
if(DEBUG) cout << "theInHit is valid = " << theInHit->isValid() << endl;
|
| 345 |
|
|
|
| 346 |
|
|
uint iidd = theInHit->geographicalId().rawId();
|
| 347 |
|
|
|
| 348 |
|
|
StripSubdetector strip=StripSubdetector(iidd);
|
| 349 |
|
|
unsigned int subid=strip.subdetId();
|
| 350 |
|
|
uint TKlayers = 0;
|
| 351 |
|
|
if (subid == StripSubdetector::TIB) {
|
| 352 |
|
|
TIBDetId tibid(iidd);
|
| 353 |
|
|
TKlayers = tibid.layer();
|
| 354 |
|
|
}
|
| 355 |
|
|
if (subid == StripSubdetector::TOB) {
|
| 356 |
|
|
TOBDetId tobid(iidd);
|
| 357 |
|
|
TKlayers = tobid.layer() + 4 ;
|
| 358 |
|
|
}
|
| 359 |
|
|
if (subid == StripSubdetector::TID) {
|
| 360 |
|
|
TIDDetId tidid(iidd);
|
| 361 |
|
|
TKlayers = tidid.wheel() + 10;
|
| 362 |
|
|
}
|
| 363 |
|
|
if (subid == StripSubdetector::TEC) {
|
| 364 |
|
|
TECDetId tecid(iidd);
|
| 365 |
|
|
TKlayers = tecid.wheel() + 13 ;
|
| 366 |
|
|
}
|
| 367 |
|
|
if ( subid==PixelSubdetector::PixelBarrel ) {
|
| 368 |
|
|
TKlayers = PXBDetId(iidd).layer() + 22;
|
| 369 |
|
|
}
|
| 370 |
|
|
if ( subid==PixelSubdetector::PixelEndcap ) {
|
| 371 |
|
|
TKlayers = PXFDetId(iidd).disk() + 25;
|
| 372 |
|
|
}
|
| 373 |
|
|
if (DEBUG) cout << "TKlayer from trajectory: " << TKlayers << " from module = " << iidd << " matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
|
| 374 |
|
|
|
| 375 |
|
|
// Make vector of TrajectoryAtValidHits to hold the trajectories
|
| 376 |
|
|
std::vector<TrajectoryAtValidHit> TMs;
|
| 377 |
|
|
|
| 378 |
|
|
// Make AnalyticalPropagator to use in TAVH constructor
|
| 379 |
|
|
AnalyticalPropagator propagator(magField_,anyDirection);
|
| 380 |
|
|
|
| 381 |
|
|
// for double sided layers check both sensors--if no hit was found on either sensor surface,
|
| 382 |
|
|
// the trajectory measurements only have one invalid hit entry on the matched surface
|
| 383 |
|
|
// so get the TrajectoryAtValidHit for both surfaces and include them in the study
|
| 384 |
|
|
if (isDoubleSided(iidd) && ((iidd & 0x3)==0) ) {
|
| 385 |
|
|
// do hit eff check twice--once for each sensor
|
| 386 |
|
|
//add a TM for each surface
|
| 387 |
|
|
TMs.push_back(TrajectoryAtValidHit(*itm,tkgeom, propagator, 1));
|
| 388 |
|
|
TMs.push_back(TrajectoryAtValidHit(*itm,tkgeom, propagator, 2));
|
| 389 |
|
|
} else if ( isDoubleSided(iidd) && (!check2DPartner(iidd, TMeas)) ) {
|
| 390 |
|
|
// if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
|
| 391 |
|
|
// the matched layer should be added to the study as well
|
| 392 |
|
|
TMs.push_back(TrajectoryAtValidHit(*itm,tkgeom, propagator, 1));
|
| 393 |
|
|
TMs.push_back(TrajectoryAtValidHit(*itm,tkgeom, propagator, 2));
|
| 394 |
|
|
if (DEBUG) cout << " found a hit with a missing partner" << endl;
|
| 395 |
|
|
} else {
|
| 396 |
|
|
//only add one TM for the single surface and the other will be added in the next iteration
|
| 397 |
|
|
TMs.push_back(TrajectoryAtValidHit(*itm,tkgeom, propagator));
|
| 398 |
|
|
}
|
| 399 |
|
|
|
| 400 |
|
|
// Modules Constraints
|
| 401 |
|
|
|
| 402 |
|
|
for(std::vector<TrajectoryAtValidHit>::const_iterator TM=TMs.begin();TM!=TMs.end();++TM) {
|
| 403 |
|
|
|
| 404 |
|
|
// --> Get trajectory from combinatedState
|
| 405 |
|
|
iidd = TM->monodet_id();
|
| 406 |
|
|
if (DEBUG) cout << "setting iidd = " << iidd << " before checking efficiency and ";
|
| 407 |
|
|
|
| 408 |
|
|
xloc = TM->localX();
|
| 409 |
|
|
yloc = TM->localY();
|
| 410 |
|
|
|
| 411 |
|
|
xErr = TM->localErrorX();
|
| 412 |
|
|
yErr = TM->localErrorY();
|
| 413 |
|
|
|
| 414 |
|
|
angleX = atan( TM->localDxDz() );
|
| 415 |
|
|
angleY = atan( TM->localDyDz() );
|
| 416 |
|
|
|
| 417 |
|
|
xglob = TM->globalX();
|
| 418 |
|
|
yglob = TM->globalY();
|
| 419 |
|
|
zglob = TM->globalZ();
|
| 420 |
|
|
withinAcceptance = TM->withinAcceptance();
|
| 421 |
|
|
|
| 422 |
|
|
if ( (layers == TKlayers) || (layers == 0) ) { // Look at the layer not used to reconstruct the track
|
| 423 |
|
|
whatlayer = TKlayers;
|
| 424 |
|
|
if (DEBUG) cout << "Looking at layer under study" << endl;
|
| 425 |
|
|
TrajGlbX = 0.0; TrajGlbY = 0.0; TrajGlbZ = 0.0; ModIsBad = 0; Id = 0; SiStripQualBad = 0;
|
| 426 |
|
|
run = 0; event = 0; TrajLocX = 0.0; TrajLocY = 0.0; TrajLocErrX = 0.0; TrajLocErrY = 0.0;
|
| 427 |
|
|
TrajLocAngleX = -999.0; TrajLocAngleY = -999.0; ResX = 0.0; ResXSig = 0.0;
|
| 428 |
|
|
ClusterLocX = 0.0; ClusterLocY = 0.0; ClusterLocErrX = 0.0; ClusterLocErrY = 0.0; ClusterStoN = 0.0;
|
| 429 |
|
|
|
| 430 |
|
|
// RPhi RecHit Efficiency
|
| 431 |
|
|
|
| 432 |
|
|
if (stripClusters.size() > 0 ) {
|
| 433 |
|
|
if (DEBUG) cout << "Checking clusters with size = " << stripClusters.size() << endl;
|
| 434 |
|
|
int nClusters = 0;
|
| 435 |
|
|
std::vector< std::vector<float> > VCluster_info; //fill with X residual, X residual pull, local X, sig(X), local Y, sig(Y), StoN
|
| 436 |
|
|
// check sistripclusters for a match
|
| 437 |
|
|
for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = stripClusters.begin(); DSViter != stripClusters.end(); DSViter++) {
|
| 438 |
|
|
// DSViter is a vector of SiStripClusters located on a single module
|
| 439 |
|
|
if (DEBUG) cout << "the ID from the DSViter over strip clusters = " << DSViter->id() << endl;
|
| 440 |
|
|
unsigned int ClusterId = DSViter->id();
|
| 441 |
|
|
if (ClusterId == iidd) {
|
| 442 |
|
|
if (DEBUG) cout << "found strip (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
|
| 443 |
|
|
DetId ClusterDetId(ClusterId);
|
| 444 |
|
|
const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
|
| 445 |
|
|
|
| 446 |
|
|
for(edmNew::DetSet<SiStripCluster>::const_iterator iter=DSViter->begin();iter!=DSViter->end();++iter) {
|
| 447 |
|
|
//iter is a single SiStripCluster
|
| 448 |
|
|
StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
|
| 449 |
|
|
float res = (parameters.first.x() - xloc);
|
| 450 |
|
|
float sigma = checkConsistency(parameters , xloc, xErr);
|
| 451 |
|
|
// The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
|
| 452 |
|
|
// you need a TransientTrackingRecHit instead of the cluster
|
| 453 |
|
|
//theEstimator= new Chi2MeasurementEstimator(30);
|
| 454 |
|
|
//const Chi2MeasurementEstimator *theEstimator(100);
|
| 455 |
|
|
//theEstimator->estimate(TM->tsos(), TransientTrackingRecHit);
|
| 456 |
|
|
|
| 457 |
|
|
//SiStripClusterInfo clusterInfo = SiStripClusterInfo(*iter, es);
|
| 458 |
|
|
// signal to noise from SiStripClusterInfo not working in 225. I'll fix this after the interface
|
| 459 |
|
|
// redesign in 300 -ku
|
| 460 |
|
|
//float cluster_info[7] = {res, sigma, parameters.first.x(), sqrt(parameters.second.xx()), parameters.first.y(), sqrt(parameters.second.yy()), signal_to_noise};
|
| 461 |
|
|
std::vector< float > cluster_info;
|
| 462 |
|
|
cluster_info.push_back(res);
|
| 463 |
|
|
cluster_info.push_back(sigma);
|
| 464 |
|
|
cluster_info.push_back(parameters.first.x());
|
| 465 |
|
|
cluster_info.push_back(sqrt(parameters.second.xx()));
|
| 466 |
|
|
cluster_info.push_back(parameters.first.y());
|
| 467 |
|
|
cluster_info.push_back(sqrt(parameters.second.yy()));
|
| 468 |
|
|
//cout << "before getting signal over noise" << endl;
|
| 469 |
|
|
//cluster_info.push_back( clusterInfo.signalOverNoise() );
|
| 470 |
|
|
//cluster_info.push_back( clusterInfo.getSignalOverNoise() );
|
| 471 |
|
|
//cout << "after getting signal over noise" << endl;
|
| 472 |
|
|
VCluster_info.push_back(cluster_info);
|
| 473 |
|
|
nClusters++;
|
| 474 |
|
|
if (DEBUG) cout << "Have ID match. residual = " << VCluster_info.back()[0] << " res sigma = " << VCluster_info.back()[1] << endl;
|
| 475 |
|
|
if (DEBUG) cout << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
|
| 476 |
|
|
if (DEBUG) cout << "hit position = " << parameters.first.x() << " hit error = " << sqrt(parameters.second.xx()) << " trajectory position = " << xloc << " traj error = " << xErr << endl;
|
| 477 |
|
|
}
|
| 478 |
|
|
}
|
| 479 |
|
|
}
|
| 480 |
|
|
|
| 481 |
|
|
// check sipixelclusters for a match
|
| 482 |
|
|
for (edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = pixelClusters.begin(); DSViter != pixelClusters.end(); DSViter++) {
|
| 483 |
|
|
// DSViter is a vector of SiPixelClusters located on a single module
|
| 484 |
|
|
if (DEBUG) cout << "the ID from the DSViter over pixel clusters = " << DSViter->id() << endl;
|
| 485 |
|
|
unsigned int ClusterId = DSViter->id();
|
| 486 |
|
|
if (ClusterId == iidd) {
|
| 487 |
|
|
|
| 488 |
|
|
if (DEBUG) cout << "found pixel (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
|
| 489 |
|
|
|
| 490 |
|
|
DetId ClusterDetId(ClusterId);
|
| 491 |
|
|
const PixelGeomDetUnit * pixeldet=(const PixelGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
|
| 492 |
|
|
for(edmNew::DetSet<SiPixelCluster>::const_iterator iter=DSViter->begin();iter!=DSViter->end();++iter) {
|
| 493 |
|
|
//iter is a single SiPixelCluster
|
| 494 |
|
|
PixelClusterParameterEstimator::LocalValues parameters=pixelcpe.localParameters(*iter,*pixeldet);
|
| 495 |
|
|
float res = (parameters.first.x() - xloc);
|
| 496 |
|
|
float sigma = checkConsistency(parameters , xloc, xErr);
|
| 497 |
|
|
|
| 498 |
|
|
std::vector< float > cluster_info;
|
| 499 |
|
|
cluster_info.push_back(res);
|
| 500 |
|
|
cluster_info.push_back(sigma);
|
| 501 |
|
|
cluster_info.push_back(parameters.first.x());
|
| 502 |
|
|
cluster_info.push_back(sqrt(parameters.second.xx()));
|
| 503 |
|
|
cluster_info.push_back(parameters.first.y());
|
| 504 |
|
|
cluster_info.push_back(sqrt(parameters.second.yy()));
|
| 505 |
|
|
//cout << "before getting signal over noise" << endl;
|
| 506 |
|
|
//cluster_info.push_back( clusterInfo.signalOverNoise() );
|
| 507 |
|
|
//cluster_info.push_back( clusterInfo.getSignalOverNoise() );
|
| 508 |
|
|
//cout << "after getting signal over noise" << endl;
|
| 509 |
|
|
VCluster_info.push_back(cluster_info);
|
| 510 |
|
|
nClusters++;
|
| 511 |
|
|
if (DEBUG) cout << "Have ID match. residual = " << VCluster_info.back()[0] << " res sigma = " << VCluster_info.back()[1] << endl;
|
| 512 |
|
|
if (DEBUG) cout << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
|
| 513 |
|
|
if (DEBUG) cout << "hit position = " << parameters.first.x() << " hit error = " << sqrt(parameters.second.xx()) << " trajectory position = " << xloc << " traj error = " << xErr << endl;
|
| 514 |
|
|
}
|
| 515 |
|
|
}
|
| 516 |
|
|
}
|
| 517 |
|
|
|
| 518 |
|
|
float FinalResSig = 1000.0;
|
| 519 |
|
|
float FinalCluster[7]= {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
|
| 520 |
|
|
if (nClusters > 0) {
|
| 521 |
|
|
if (DEBUG) cout << "found clusters > 0" << endl;
|
| 522 |
|
|
if (nClusters > 1) {
|
| 523 |
|
|
//get the smallest one
|
| 524 |
|
|
vector< vector<float> >::iterator ires;
|
| 525 |
|
|
for (ires=VCluster_info.begin(); ires!=VCluster_info.end(); ires++){
|
| 526 |
|
|
if ( abs((*ires)[1]) < abs(FinalResSig)) {
|
| 527 |
|
|
FinalResSig = (*ires)[1];
|
| 528 |
|
|
for (uint i = 0; i<ires->size(); i++) {
|
| 529 |
|
|
if (DEBUG) cout << "filling final cluster. i = " << i << " before fill FinalCluster[i]=" << FinalCluster[i] << " and (*ires)[i] =" << (*ires)[i] << endl;
|
| 530 |
|
|
if (DEBUG) FinalCluster[i] = (*ires)[i];
|
| 531 |
|
|
if (DEBUG) cout << "filling final cluster. i = " << i << " after fill FinalCluster[i]=" << FinalCluster[i] << " and (*ires)[i] =" << (*ires)[i] << endl;
|
| 532 |
|
|
}
|
| 533 |
|
|
}
|
| 534 |
|
|
if (DEBUG) cout << "iresidual = " << (*ires)[0] << " isigma = " << (*ires)[1] << " and FinalRes = " << FinalCluster[0] << endl;
|
| 535 |
|
|
}
|
| 536 |
|
|
}
|
| 537 |
|
|
else {
|
| 538 |
|
|
FinalResSig = VCluster_info.at(0)[1];
|
| 539 |
|
|
for (uint i = 0; i<VCluster_info.at(0).size(); i++) {
|
| 540 |
|
|
FinalCluster[i] = VCluster_info.at(0)[i];
|
| 541 |
|
|
}
|
| 542 |
|
|
}
|
| 543 |
|
|
nClusters=0;
|
| 544 |
|
|
VCluster_info.clear();
|
| 545 |
|
|
}
|
| 546 |
|
|
|
| 547 |
|
|
if (DEBUG) cout << "Final residual in X = " << FinalCluster[0] << "+-" << (FinalCluster[0]/FinalResSig) << endl;
|
| 548 |
|
|
if (DEBUG) cout << "Checking location of trajectory: abs(yloc) = " << abs(yloc) << " abs(xloc) = " << abs(xloc) << endl;
|
| 549 |
|
|
if (DEBUG) cout << "Checking location of cluster hit: yloc = " << FinalCluster[4] << "+-" << FinalCluster[5] << " xloc = " << FinalCluster[2] << "+-" << FinalCluster[3] << endl;
|
| 550 |
|
|
if (DEBUG) cout << "Final cluster signal to noise = " << FinalCluster[6] << endl;
|
| 551 |
|
|
|
| 552 |
|
|
float exclusionWidth = 0.4;
|
| 553 |
|
|
float TOBexclusion = 0.0;
|
| 554 |
|
|
float TECexRing5 = -0.89;
|
| 555 |
|
|
float TECexRing6 = -0.56;
|
| 556 |
|
|
float TECexRing7 = 0.60;
|
| 557 |
|
|
//Added by Chris Edelmaier to do TEC bonding exclusion
|
| 558 |
|
|
int subdetector = ((iidd>>25) & 0x7);
|
| 559 |
|
|
int ringnumber = ((iidd>>5) & 0x7);
|
| 560 |
|
|
|
| 561 |
|
|
//New TOB and TEC bonding region exclusion zone
|
| 562 |
|
|
if((TKlayers >= 5 && TKlayers < 11)||((subdetector == 6)&&( (ringnumber >= 5)&&(ringnumber <=7) ))) {
|
| 563 |
|
|
//There are only 2 cases that we need to exclude for
|
| 564 |
|
|
float highzone = 0.0;
|
| 565 |
|
|
float lowzone = 0.0;
|
| 566 |
|
|
float higherr = yloc + 5.0*yErr;
|
| 567 |
|
|
float lowerr = yloc - 5.0*yErr;
|
| 568 |
|
|
if(TKlayers >= 5 && TKlayers < 11) {
|
| 569 |
|
|
//TOB zone
|
| 570 |
|
|
highzone = TOBexclusion + exclusionWidth;
|
| 571 |
|
|
lowzone = TOBexclusion - exclusionWidth;
|
| 572 |
|
|
} else if (ringnumber == 5) {
|
| 573 |
|
|
//TEC ring 5
|
| 574 |
|
|
highzone = TECexRing5 + exclusionWidth;
|
| 575 |
|
|
lowzone = TECexRing5 - exclusionWidth;
|
| 576 |
|
|
} else if (ringnumber == 6) {
|
| 577 |
|
|
//TEC ring 6
|
| 578 |
|
|
highzone = TECexRing6 + exclusionWidth;
|
| 579 |
|
|
lowzone = TECexRing6 - exclusionWidth;
|
| 580 |
|
|
} else if (ringnumber == 7) {
|
| 581 |
|
|
//TEC ring 7
|
| 582 |
|
|
highzone = TECexRing7 + exclusionWidth;
|
| 583 |
|
|
lowzone = TECexRing7 - exclusionWidth;
|
| 584 |
|
|
}
|
| 585 |
|
|
//Now that we have our exclusion region, we just have to properly identify it
|
| 586 |
|
|
if((highzone <= higherr)&&(highzone >= lowerr)) withinAcceptance = false;
|
| 587 |
|
|
if((lowzone >= lowerr)&&(lowzone <= higherr)) withinAcceptance = false;
|
| 588 |
|
|
if((higherr <= highzone)&&(higherr >= lowzone)) withinAcceptance = false;
|
| 589 |
|
|
if((lowerr >= lowzone)&&(lowerr <= highzone)) withinAcceptance = false;
|
| 590 |
|
|
}
|
| 591 |
|
|
|
| 592 |
|
|
// fill ntuple varibles
|
| 593 |
|
|
//get global position from module id number iidd
|
| 594 |
|
|
TrajGlbX = xglob;
|
| 595 |
|
|
TrajGlbY = yglob;
|
| 596 |
|
|
TrajGlbZ = zglob;
|
| 597 |
|
|
Id = iidd;
|
| 598 |
|
|
run = run_nr;
|
| 599 |
|
|
event = ev_nr;
|
| 600 |
|
|
//if ( SiStripQuality_->IsModuleBad(iidd) ) {
|
| 601 |
|
|
if ( SiStripQuality_->getBadApvs(iidd)!=0 ) {
|
| 602 |
|
|
SiStripQualBad = 1;
|
| 603 |
|
|
if(DEBUG) cout << "strip is bad from SiStripQuality" << endl;
|
| 604 |
|
|
} else {
|
| 605 |
|
|
SiStripQualBad = 0;
|
| 606 |
|
|
if(DEBUG) cout << "strip is good from SiStripQuality" << endl;
|
| 607 |
|
|
}
|
| 608 |
|
|
|
| 609 |
|
|
TrajLocX = xloc;
|
| 610 |
|
|
TrajLocY = yloc;
|
| 611 |
|
|
TrajLocErrX = xErr;
|
| 612 |
|
|
TrajLocErrY = yErr;
|
| 613 |
|
|
TrajLocAngleX = angleX;
|
| 614 |
|
|
TrajLocAngleY = angleY;
|
| 615 |
|
|
ResX = FinalCluster[0];
|
| 616 |
|
|
ResXSig = FinalResSig;
|
| 617 |
|
|
if (FinalResSig != FinalCluster[1]) if (DEBUG) cout << "Problem with best cluster selection because FinalResSig = " << FinalResSig << " and FinalCluster[1] = " << FinalCluster[1] << endl;
|
| 618 |
|
|
ClusterLocX = FinalCluster[2];
|
| 619 |
|
|
ClusterLocY = FinalCluster[4];
|
| 620 |
|
|
ClusterLocErrX = FinalCluster[3];
|
| 621 |
|
|
ClusterLocErrY = FinalCluster[5];
|
| 622 |
|
|
ClusterStoN = FinalCluster[6];
|
| 623 |
|
|
|
| 624 |
|
|
if (DEBUG) cout << "before check good" << endl;
|
| 625 |
|
|
|
| 626 |
|
|
if ( FinalResSig < 999.0) { //could make requirement on track/hit consistency, but for
|
| 627 |
|
|
//now take anything with a hit on the module
|
| 628 |
|
|
if (DEBUG) cout << "hit being counted as good" << endl;
|
| 629 |
|
|
ModIsBad = 0;
|
| 630 |
|
|
traj->Fill();
|
| 631 |
|
|
}
|
| 632 |
|
|
else {
|
| 633 |
|
|
if (DEBUG) cout << "hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0] << " FinalRecHit " <<
|
| 634 |
|
|
iidd << " TKlayers " << TKlayers << " xloc " << xloc << " yloc " << yloc << " module " << iidd <<
|
| 635 |
|
|
" matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
|
| 636 |
|
|
ModIsBad = 1;
|
| 637 |
|
|
traj->Fill();
|
| 638 |
|
|
|
| 639 |
|
|
if (DEBUG) cout << " RPhi Error " << sqrt(xErr*xErr + yErr*yErr) << " ErrorX " << xErr << " yErr " << yErr << endl;
|
| 640 |
|
|
} if (DEBUG) cout << "after good location check" << endl;
|
| 641 |
|
|
} if (DEBUG) cout << "after list of clusters" << endl;
|
| 642 |
|
|
} if (DEBUG) cout << "After layers=TKLayers if" << endl;
|
| 643 |
|
|
} if (DEBUG) cout << "After looping over TrajAtValidHit list" << endl;
|
| 644 |
|
|
} if (DEBUG) cout << "end TMeasurement loop" << endl;
|
| 645 |
|
|
}
|
| 646 |
|
|
}
|
| 647 |
|
|
|
| 648 |
|
|
void HitEff::endJob(){
|
| 649 |
|
|
traj->GetDirectory()->cd();
|
| 650 |
|
|
traj->Write();
|
| 651 |
|
|
|
| 652 |
|
|
cout << " Events Analysed " << events << endl;
|
| 653 |
|
|
cout << " Number Of Tracked events " << EventTrackCKF << endl;
|
| 654 |
|
|
}
|
| 655 |
|
|
|
| 656 |
|
|
double HitEff::checkConsistency(StripClusterParameterEstimator::LocalValues parameters, double xx, double xerr) {
|
| 657 |
|
|
double error = sqrt(parameters.second.xx() + xerr*xerr);
|
| 658 |
|
|
double separation = abs(parameters.first.x() - xx);
|
| 659 |
|
|
double consistency = separation/error;
|
| 660 |
|
|
return consistency;
|
| 661 |
|
|
}
|
| 662 |
|
|
|
| 663 |
|
|
double HitEff::checkConsistency(const SiStripRecHit2D* rechit, double xx, double xerr)
|
| 664 |
|
|
{
|
| 665 |
|
|
double error = sqrt(rechit->localPositionError().xx() + xerr*xerr);
|
| 666 |
|
|
double separation = abs(rechit->localPosition().x() - xx);
|
| 667 |
|
|
double consistency = separation/error;
|
| 668 |
|
|
return consistency;
|
| 669 |
|
|
}
|
| 670 |
|
|
|
| 671 |
|
|
bool HitEff::isDoubleSided(uint iidd) const {
|
| 672 |
|
|
StripSubdetector strip=StripSubdetector(iidd);
|
| 673 |
|
|
unsigned int subid=strip.subdetId();
|
| 674 |
|
|
uint layer = 0;
|
| 675 |
|
|
if (subid == StripSubdetector::TIB) {
|
| 676 |
|
|
TIBDetId tibid(iidd);
|
| 677 |
|
|
layer = tibid.layer();
|
| 678 |
|
|
if (layer == 1 || layer == 2) return true;
|
| 679 |
|
|
else return false;
|
| 680 |
|
|
}
|
| 681 |
|
|
else if (subid == StripSubdetector::TOB) {
|
| 682 |
|
|
TOBDetId tobid(iidd);
|
| 683 |
|
|
layer = tobid.layer() + 4 ;
|
| 684 |
|
|
if (layer == 5 || layer == 6) return true;
|
| 685 |
|
|
else return false;
|
| 686 |
|
|
}
|
| 687 |
|
|
else if (subid == StripSubdetector::TID) {
|
| 688 |
|
|
TIDDetId tidid(iidd);
|
| 689 |
|
|
layer = tidid.ring() + 10;
|
| 690 |
|
|
if (layer == 11 || layer == 12) return true;
|
| 691 |
|
|
else return false;
|
| 692 |
|
|
}
|
| 693 |
|
|
else if (subid == StripSubdetector::TEC) {
|
| 694 |
|
|
TECDetId tecid(iidd);
|
| 695 |
|
|
layer = tecid.ring() + 13 ;
|
| 696 |
|
|
if (layer == 14 || layer == 15 || layer == 18) return true;
|
| 697 |
|
|
else return false;
|
| 698 |
|
|
}
|
| 699 |
|
|
else
|
| 700 |
|
|
return false;
|
| 701 |
|
|
}
|
| 702 |
|
|
|
| 703 |
|
|
bool HitEff::check2DPartner(uint iidd, std::vector<TrajectoryMeasurement> traj) {
|
| 704 |
|
|
uint partner_iidd = 0;
|
| 705 |
|
|
bool found2DPartner = false;
|
| 706 |
|
|
// first get the id of the other detector
|
| 707 |
|
|
if ((iidd & 0x3)==1) partner_iidd = iidd+1;
|
| 708 |
|
|
if ((iidd & 0x3)==2) partner_iidd = iidd-1;
|
| 709 |
|
|
// next look in the trajectory measurements for a measurement from that detector
|
| 710 |
|
|
// loop through trajectory measurements to find the partner_iidd
|
| 711 |
|
|
for (std::vector<TrajectoryMeasurement>::const_iterator iTM=traj.begin(); iTM!=traj.end(); ++iTM) {
|
| 712 |
|
|
if (iTM->recHit()->geographicalId().rawId()==partner_iidd) {
|
| 713 |
|
|
found2DPartner = true;
|
| 714 |
|
|
}
|
| 715 |
|
|
}
|
| 716 |
|
|
return found2DPartner;
|
| 717 |
|
|
}
|
| 718 |
|
|
|
| 719 |
|
|
//define this as a plug-in
|
| 720 |
|
|
DEFINE_FWK_MODULE(HitEff);
|
