File indexing completed on 2024-04-06 12:28:03
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <vector>
0022
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 #include "DataFormats/TrackReco/interface/Track.h"
0033 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0034 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0035 #include <TH1.h>
0036 #include <TH2.h>
0037 #include <TFile.h>
0038 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0039 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0040 #include "MagneticField/Engine/interface/MagneticField.h"
0041 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0043 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0046 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0047 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0048 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0049
0050
0051
0052
0053
0054 class TestOutliers : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0055 public:
0056 explicit TestOutliers(const edm::ParameterSet &);
0057 ~TestOutliers() override;
0058
0059 private:
0060 void beginRun(edm::Run const &run, const edm::EventSetup &) override;
0061 void analyze(const edm::Event &, const edm::EventSetup &) override;
0062 void endRun(edm::Run const &run, const edm::EventSetup &) override {}
0063 void endJob() override;
0064
0065
0066 edm::InputTag trackTagsOut_;
0067 edm::InputTag trackTagsOld_;
0068 edm::InputTag tpTags_;
0069 TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0070 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOldToken;
0071 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOutToken;
0072 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> theGToken;
0073 edm::ESHandle<TrackerGeometry> theG;
0074 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> theTopoToken;
0075 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken;
0076 std::string out;
0077 TFile *file;
0078 TH1F *histoPtOut, *histoPtOld;
0079 TH1F *histoQoverpOut, *histoQoverpOld;
0080 TH1F *histoPhiOut, *histoPhiOld;
0081 TH1F *histoD0Out, *histoD0Old;
0082 TH1F *histoDzOut, *histoDzOld;
0083 TH1F *histoLambdaOut, *histoLambdaOld;
0084 TH1F *deltahits, *deltahitsAssocGained, *deltahitsAssocLost, *hitsPerTrackLost, *hitsPerTrackAssocLost,
0085 *hitsPerTrackGained, *hitsPerTrackAssocGained, *deltahitsOK, *deltahitsNO;
0086 TH1F *okcutsOut, *okcutsOld;
0087 TH1F *tracks, *goodbadhits, *process, *goodcluster, *goodprocess, *badcluster, *badprocess;
0088 TH1F *goodhittype, *goodlayer, *goodhittype_clgteq4, *goodlayer_clgteq4, *goodhittype_cllt4, *goodlayer_cllt4,
0089 *badhittype, *badlayer;
0090 TH2F *posxy, *poszr;
0091 TH1F *goodpixgteq4_simvecsize, *goodpixlt4_simvecsize;
0092 TH1F *goodst1gteq4_simvecsize, *goodst1lt4_simvecsize;
0093 TH1F *goodst2gteq4_simvecsize, *goodst2lt4_simvecsize;
0094 TH1F *goodprjgteq4_simvecsize, *goodprjlt4_simvecsize;
0095 TH1F *goodpix_clustersize;
0096 TH1F *goodst1_clustersize;
0097 TH1F *goodst2_clustersize;
0098 TH1F *goodprj_clustersize;
0099 TH1F *goodpix_simvecsize;
0100 TH1F *goodst1_simvecsize;
0101 TH1F *goodst2_simvecsize;
0102 TH1F *goodprj_simvecsize;
0103 TH1F *goodhittype_simvecsmall, *goodlayer_simvecsmall, *goodhittype_simvecbig, *goodlayer_simvecbig, *goodbadmerged,
0104 *goodbadmergedLost, *goodbadmergedGained;
0105 TH1F *energyLoss, *energyLossMax, *energyLossRatio, *nOfTrackIds, *mergedPull, *mergedlayer, *mergedcluster,
0106 *mergedhittype;
0107 TH1F *sizeOut, *sizeOld, *sizeOutT, *sizeOldT;
0108 TH1F *countOutA, *countOutT, *countOldT;
0109 TH1F *gainedhits, *gainedhits2;
0110 TH1F *probXgood, *probXbad, *probXdelta, *probXshared, *probXnoshare;
0111 TH1F *probYgood, *probYbad, *probYdelta, *probYshared, *probYnoshare;
0112 };
0113
0114
0115
0116
0117
0118
0119
0120
0121 using namespace edm;
0122
0123
0124
0125
0126 TestOutliers::TestOutliers(const edm::ParameterSet &iConfig)
0127 : trackTagsOut_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOut")),
0128 trackTagsOld_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOld")),
0129 tpTags_(iConfig.getUntrackedParameter<edm::InputTag>("tp")),
0130 trackerHitAssociatorConfig_(consumesCollector()),
0131 theAssociatorOldToken(consumes<reco::TrackToTrackingParticleAssociator>(
0132 iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOld"))),
0133 theAssociatorOutToken(consumes<reco::TrackToTrackingParticleAssociator>(
0134 iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOut"))),
0135 theGToken(esConsumes<edm::Transition::BeginRun>()),
0136 theTopoToken(esConsumes()),
0137 theMFToken(esConsumes()),
0138 out(iConfig.getParameter<std::string>("out")) {
0139 LogTrace("TestOutliers") << "constructor";
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149 LogTrace("TestOutliers") << "end constructor";
0150 }
0151
0152 TestOutliers::~TestOutliers() {
0153
0154
0155 }
0156
0157
0158
0159
0160
0161
0162 void TestOutliers::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0163
0164 edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(theTopoToken);
0165
0166 theG = iSetup.getHandle(theGToken);
0167 edm::ESHandle<MagneticField> theMF = iSetup.getHandle(theMFToken);
0168
0169 using namespace edm;
0170 using namespace std;
0171 using reco::TrackCollection;
0172
0173 LogTrace("TestOutliers") << "analyze event #" << iEvent.id();
0174
0175 edm::Handle<edm::View<reco::Track> > tracksOut;
0176 iEvent.getByLabel(trackTagsOut_, tracksOut);
0177 edm::Handle<edm::View<reco::Track> > tracksOld;
0178 iEvent.getByLabel(trackTagsOld_, tracksOld);
0179 Handle<TrackingParticleCollection> tps;
0180 iEvent.getByLabel(tpTags_, tps);
0181 edm::Handle<reco::BeamSpot> beamSpot;
0182 iEvent.getByLabel("offlineBeamSpot", beamSpot);
0183
0184 TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);
0185
0186 edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOld;
0187 iEvent.getByToken(theAssociatorOldToken, hAssociatorOld);
0188 const reco::TrackToTrackingParticleAssociator *theAssociatorOld = hAssociatorOld.product();
0189
0190 edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOut;
0191 iEvent.getByToken(theAssociatorOutToken, hAssociatorOut);
0192 const reco::TrackToTrackingParticleAssociator *theAssociatorOut = hAssociatorOut.product();
0193
0194 reco::RecoToSimCollection recSimCollOut = theAssociatorOut->associateRecoToSim(tracksOut, tps);
0195 reco::RecoToSimCollection recSimCollOld = theAssociatorOld->associateRecoToSim(tracksOld, tps);
0196 sizeOut->Fill(recSimCollOut.size());
0197 sizeOld->Fill(recSimCollOld.size());
0198 sizeOutT->Fill(tracksOut->size());
0199 sizeOldT->Fill(tracksOld->size());
0200
0201 LogTrace("TestOutliers") << "tOld size=" << tracksOld->size() << " tOut size=" << tracksOut->size()
0202 << " aOld size=" << recSimCollOld.size() << " aOut size=" << recSimCollOut.size();
0203
0204 #if 0
0205 LogTrace("TestOutliers") << "recSimCollOld.size()=" << recSimCollOld.size() ;
0206 for(reco::TrackCollection::size_type j=0; j<tracksOld->size(); ++j){
0207 reco::TrackRef trackOld(tracksOld, j);
0208
0209 LogTrace("TestOutliers") << "trackOld->pt()=" << trackOld->pt() << " trackOld->numberOfValidHits()=" << trackOld->numberOfValidHits();
0210 std::vector<pair<TrackingParticleRef, double> > tpOld;
0211 if(recSimCollOld.find(trackOld) != recSimCollOld.end()){
0212 tpOld = recSimCollOld[trackOld];
0213 if (tpOld.size()!=0) LogTrace("TestOutliers") << " associated" ;
0214 else LogTrace("TestOutliers") << " NOT associated" ;
0215 } else LogTrace("TestOutliers") << " NOT associated" ;
0216 }
0217 LogTrace("TestOutliers") << "recSimCollOut.size()=" << recSimCollOut.size() ;
0218 for(reco::TrackCollection::size_type j=0; j<tracksOut->size(); ++j){
0219 reco::TrackRef trackOut(tracksOut, j);
0220
0221 LogTrace("TestOutliers") << "trackOut->pt()=" << trackOut->pt() << " trackOut->numberOfValidHits()=" << trackOut->numberOfValidHits();
0222 std::vector<pair<TrackingParticleRef, double> > tpOut;
0223 if(recSimCollOut.find(trackOut) != recSimCollOut.end()){
0224 tpOut = recSimCollOut[trackOut];
0225 if (tpOut.size()!=0) LogTrace("TestOutliers") << " associated" ;
0226 else LogTrace("TestOutliers") << " NOT associated" ;
0227 } else LogTrace("TestOutliers") << " NOT associated" ;
0228 }
0229 #endif
0230
0231 LogTrace("TestOutliers") << "tracksOut->size()=" << tracksOut->size();
0232 LogTrace("TestOutliers") << "tracksOld->size()=" << tracksOld->size();
0233
0234 std::vector<unsigned int> outused;
0235 for (unsigned int j = 0; j < tracksOld->size(); ++j) {
0236 countOldT->Fill(1);
0237 edm::RefToBase<reco::Track> trackOld(tracksOld, j);
0238 LogTrace("TestOutliers") << "now track old with id=" << j << " seed ref=" << trackOld->seedRef().get()
0239 << " pt=" << trackOld->pt() << " eta=" << trackOld->eta()
0240 << " chi2=" << trackOld->normalizedChi2()
0241 << " tip= " << fabs(trackOld->dxy(beamSpot->position()))
0242 << " lip= " << fabs(trackOld->dsz(beamSpot->position()));
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257 std::vector<unsigned int> outtest;
0258 for (unsigned int k = 0; k < tracksOut->size(); ++k) {
0259 edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, k);
0260 if (tmpOut->seedRef() == trackOld->seedRef()) {
0261 outtest.push_back(k);
0262 }
0263 }
0264
0265 edm::RefToBase<reco::Track> trackOut;
0266 if (outtest.size() == 1) {
0267 trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[0]);
0268 LogTrace("TestOutliers") << "now track out with id=" << outtest[0] << " seed ref=" << trackOut->seedRef().get()
0269 << " pt=" << trackOut->pt();
0270 outused.push_back(outtest[0]);
0271 } else if (outtest.size() > 1) {
0272 for (unsigned int p = 0; p < outtest.size(); ++p) {
0273 edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
0274 bool allhits = true;
0275 for (trackingRecHit_iterator itOut = tmpOut->recHitsBegin(); itOut != tmpOut->recHitsEnd(); itOut++) {
0276 if ((*itOut)->isValid()) {
0277 bool thishit = false;
0278 for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0279 if ((*itOld)->isValid()) {
0280 const TrackingRecHit *kt = &(**itOld);
0281 if ((*itOut)->sharesInput(kt, TrackingRecHit::some)) {
0282 thishit = true;
0283 break;
0284 }
0285 }
0286 }
0287 if (!thishit)
0288 allhits = false;
0289 }
0290 }
0291 if (allhits) {
0292 trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
0293 LogTrace("TestOutliers") << "now track out with id=" << outtest[p]
0294 << " seed ref=" << trackOut->seedRef().get() << " pt=" << trackOut->pt();
0295 outused.push_back(outtest[p]);
0296 }
0297 }
0298 }
0299
0300 if (outtest.empty() || trackOut.get() == nullptr) {
0301 if (recSimCollOld.find(trackOld) != recSimCollOld.end()) {
0302 vector<pair<TrackingParticleRef, double> > tpOld;
0303 tpOld = recSimCollOld[trackOld];
0304 if (!tpOld.empty()) {
0305 LogTrace("TestOutliers") << "no match: old associated and out lost! old has #hits="
0306 << trackOld->numberOfValidHits() << " and fraction=" << tpOld.begin()->second;
0307 if (tpOld.begin()->second > 0.5)
0308 hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
0309 }
0310 }
0311 LogTrace("TestOutliers") << "...skip to next old track";
0312 continue;
0313 }
0314
0315
0316 LogTrace("TestOutliers") << "trackOut->seedRef()=" << trackOut->seedRef().get()
0317 << " trackOld->seedRef()=" << trackOld->seedRef().get();
0318 bool oldAssoc = recSimCollOld.find(trackOld) != recSimCollOld.end();
0319 bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
0320 LogTrace("TestOutliers") << "outAssoc=" << outAssoc << " oldAssoc=" << oldAssoc;
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338 countOutT->Fill(1);
0339
0340
0341
0342 tracks->Fill(0);
0343
0344 TrackingParticleRef tpr;
0345 TrackingParticleRef tprOut;
0346 TrackingParticleRef tprOld;
0347 double fracOut;
0348 std::vector<unsigned int> tpids;
0349 std::vector<std::pair<TrackingParticleRef, double> > tpOut;
0350 std::vector<pair<TrackingParticleRef, double> > tpOld;
0351
0352
0353 if (outAssoc) {
0354 tpOut = recSimCollOut[trackOut];
0355 if (!tpOut.empty()) {
0356 countOutA->Fill(1);
0357 tprOut = tpOut.begin()->first;
0358 fracOut = tpOut.begin()->second;
0359 for (TrackingParticle::g4t_iterator g4T = tprOut->g4Track_begin(); g4T != tprOut->g4Track_end(); ++g4T) {
0360 tpids.push_back(g4T->trackId());
0361 }
0362 }
0363 }
0364
0365 if (oldAssoc) {
0366 tpOld = recSimCollOld[trackOld];
0367 if (!tpOld.empty()) {
0368 tprOld = tpOld.begin()->first;
0369
0370
0371
0372 if (tpOut.empty()) {
0373 for (TrackingParticle::g4t_iterator g4T = tprOld->g4Track_begin(); g4T != tprOld->g4Track_end(); ++g4T) {
0374 tpids.push_back(g4T->trackId());
0375 }
0376 }
0377 }
0378 }
0379
0380 if (tprOut.get() != nullptr || tprOld.get() != nullptr) {
0381
0382 tpr = tprOut.get() != nullptr ? tprOut : tprOld;
0383
0384 const SimTrack *assocTrack = &(*tpr->g4Track_begin());
0385
0386
0387 if (trackOut->numberOfValidHits() != trackOld->numberOfValidHits() ||
0388 !(*trackOut->recHitsBegin())->sharesInput((*trackOld->recHitsBegin()), TrackingRecHit::some) ||
0389 !(*(trackOut->recHitsEnd() - 1))
0390 ->sharesInput(
0391 (*(trackOld->recHitsEnd() - 1)),
0392 TrackingRecHit::
0393 some)) {
0394 LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
0395 tracks->Fill(1);
0396 LogTrace("TestOutliers")
0397 << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
0398 << sqrt(tpr->momentum().perp2())
0399
0400 << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
0401 << trackOut->numberOfValidHits()
0402
0403 << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0404
0405
0406 double PtPullOut = (trackOut->pt() - sqrt(tpr->momentum().perp2())) / trackOut->ptError();
0407 double PtPullOld = (trackOld->pt() - sqrt(tpr->momentum().perp2())) / trackOld->ptError();
0408 histoPtOut->Fill(PtPullOut);
0409 histoPtOld->Fill(PtPullOld);
0410
0411
0412 FreeTrajectoryState ftsAtProduction(
0413 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0414 GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
0415 TrackCharge(trackOld->charge()),
0416 theMF.product());
0417 TSCPBuilderNoMaterial tscpBuilder;
0418 TrajectoryStateClosestToPoint tsAtClosestApproach =
0419 tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));
0420 GlobalPoint v = tsAtClosestApproach.theState().position();
0421 GlobalVector p = tsAtClosestApproach.theState().momentum();
0422
0423
0424 double qoverpSim = tsAtClosestApproach.charge() / p.mag();
0425 double lambdaSim = M_PI / 2 - p.theta();
0426 double phiSim = p.phi();
0427 double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0428 double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0429 double d0Sim = -dxySim;
0430 double dzSim = dszSim * p.mag() / p.perp();
0431
0432
0433 double qoverpPullOut = (trackOut->qoverp() - qoverpSim) / trackOut->qoverpError();
0434 double qoverpPullOld = (trackOld->qoverp() - qoverpSim) / trackOld->qoverpError();
0435 double lambdaPullOut = (trackOut->lambda() - lambdaSim) / trackOut->thetaError();
0436 double lambdaPullOld = (trackOld->lambda() - lambdaSim) / trackOld->thetaError();
0437 double phi0PullOut = (trackOut->phi() - phiSim) / trackOut->phiError();
0438 double phi0PullOld = (trackOld->phi() - phiSim) / trackOld->phiError();
0439 double d0PullOut = (trackOut->d0() - d0Sim) / trackOut->d0Error();
0440 double d0PullOld = (trackOld->d0() - d0Sim) / trackOld->d0Error();
0441 double dzPullOut = (trackOut->dz() - dzSim) / trackOut->dzError();
0442 double dzPullOld = (trackOld->dz() - dzSim) / trackOld->dzError();
0443
0444
0445 histoQoverpOut->Fill(qoverpPullOut);
0446 histoQoverpOld->Fill(qoverpPullOld);
0447 histoPhiOut->Fill(phi0PullOut);
0448 histoPhiOld->Fill(phi0PullOld);
0449 histoD0Out->Fill(d0PullOut);
0450 histoD0Old->Fill(d0PullOld);
0451 histoDzOut->Fill(dzPullOut);
0452 histoDzOld->Fill(dzPullOld);
0453 histoLambdaOut->Fill(lambdaPullOut);
0454 histoLambdaOld->Fill(lambdaPullOld);
0455
0456
0457 LogTrace("TestOutliers") << "deltahits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0458 deltahits->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0459
0460 if (tprOut.get() != nullptr && tprOld.get() == nullptr) {
0461 if (!tpOld.empty() && tpOld.begin()->second <= 0.5) {
0462 deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0463 hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
0464 LogTrace("TestOutliers") << "a) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
0465 << " old #hits=" << trackOld->numberOfValidHits();
0466 } else {
0467 deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0468 hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
0469 LogTrace("TestOutliers") << "b) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
0470 << " old #hits=" << trackOld->numberOfValidHits();
0471 }
0472 } else if (tprOut.get() == nullptr && tprOld.get() != nullptr) {
0473 LogTrace("TestOutliers") << "old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
0474 << " and fraction=" << tpOld.begin()->second;
0475 if (tpOld.begin()->second > 0.5) {
0476 hitsPerTrackAssocLost->Fill(trackOld->numberOfValidHits());
0477 deltahitsAssocLost->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0478 }
0479 }
0480
0481 if (fabs(PtPullOut) < fabs(PtPullOld))
0482 deltahitsOK->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0483 else
0484 deltahitsNO->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0485
0486
0487
0488
0489
0490 LogTrace("TestOutliers") << "track old";
0491 for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0492 LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOld)->geographicalId().rawId();
0493 }
0494 LogTrace("TestOutliers") << "track out";
0495 for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0496 LogTrace("TestOutliers") << (*itOut)->isValid() << " " << (*itOut)->geographicalId().rawId();
0497 }
0498
0499
0500 vector<pair<int, trackingRecHit_iterator> > gainedlostoutliers;
0501
0502 for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0503 bool gained = true;
0504 if ((*itOut)->isValid()) {
0505 for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0506 if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId())
0507 gained = false;
0508 }
0509 if (gained) {
0510 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(1, itOut));
0511 LogTrace("TestOutliers") << "broken trajectory during old fit... gained hit "
0512 << (*itOut)->geographicalId().rawId();
0513 gainedhits->Fill(1);
0514 }
0515 }
0516 }
0517
0518
0519 for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0520 bool outlier = false;
0521 bool lost = true;
0522
0523 for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0524 if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
0525 lost = false;
0526 if ((*itOld)->isValid() && !(*itOut)->isValid() &&
0527 (*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
0528 LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOut)->isValid() << " "
0529 << (*itOld)->geographicalId().rawId() << " "
0530 << (*itOut)->geographicalId().rawId();
0531 outlier = true;
0532 }
0533 }
0534 }
0535 if (lost)
0536 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(2, itOld));
0537 if (lost)
0538 LogTrace("TestOutliers") << "lost";
0539 else if (outlier)
0540 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(3, itOld));
0541 }
0542
0543 for (std::vector<pair<int, trackingRecHit_iterator> >::iterator it = gainedlostoutliers.begin();
0544 it != gainedlostoutliers.end();
0545 ++it) {
0546 LogTrace("TestOutliers") << "type of processed hit:" << it->first;
0547 trackingRecHit_iterator itHit = it->second;
0548 bool gained = false;
0549 bool lost = false;
0550 bool outlier = false;
0551 if (it->first == 1)
0552 gained = true;
0553 else if (it->first == 2)
0554 lost = true;
0555 else if (it->first == 3)
0556 outlier = true;
0557
0558 if (outlier || lost || gained) {
0559
0560
0561 if (lost && (*itHit)->isValid() == false) {
0562 goodbadmergedLost->Fill(0);
0563 LogTrace("TestOutliers") << "lost invalid";
0564 continue;
0565 } else if (gained && (*itHit)->isValid() == false) {
0566 goodbadmergedGained->Fill(0);
0567 LogTrace("TestOutliers") << "gained invalid";
0568 continue;
0569 }
0570
0571
0572
0573 std::vector<SimHitIdpr> simTrackIds = hitAssociator.associateHitId(**itHit);
0574 bool goodhit = false;
0575 for (size_t j = 0; j < simTrackIds.size(); j++) {
0576 for (size_t jj = 0; jj < tpids.size(); jj++) {
0577 if (simTrackIds[j].first == tpids[jj])
0578 goodhit = true;
0579 break;
0580 }
0581 }
0582
0583
0584 int clustersize = 0;
0585 int hittypeval = 0;
0586 int layerval = 0;
0587 if (dynamic_cast<const SiPixelRecHit *>(&**itHit)) {
0588 LogTrace("TestOutliers") << "SiPixelRecHit";
0589 clustersize = ((const SiPixelRecHit *)(&**itHit))->cluster()->size();
0590 hittypeval = 1;
0591 } else if (dynamic_cast<const SiStripRecHit2D *>(&**itHit)) {
0592 LogTrace("TestOutliers") << "SiStripRecHit2D";
0593 clustersize = ((const SiStripRecHit2D *)(&**itHit))->cluster()->amplitudes().size();
0594 hittypeval = 2;
0595 } else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit)) {
0596 LogTrace("TestOutliers") << "SiStripMatchedRecHit2D";
0597 int clsize1 = ((const SiStripMatchedRecHit2D *)(&**itHit))->monoCluster().amplitudes().size();
0598 int clsize2 = ((const SiStripMatchedRecHit2D *)(&**itHit))->stereoCluster().amplitudes().size();
0599 if (clsize1 > clsize2)
0600 clustersize = clsize1;
0601 else
0602 clustersize = clsize2;
0603 hittypeval = 3;
0604 } else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&**itHit)) {
0605 LogTrace("TestOutliers") << "ProjectedSiStripRecHit2D";
0606 clustersize =
0607 ((const ProjectedSiStripRecHit2D *)(&**itHit))->originalHit().cluster()->amplitudes().size();
0608 hittypeval = 4;
0609 }
0610
0611
0612 int subdetId = (*itHit)->geographicalId().subdetId();
0613 DetId id = (*itHit)->geographicalId();
0614 int layerId = tTopo->layer(id);
0615 layerval = subdetId * 10 + layerId;
0616
0617
0618 GlobalPoint gpos = theG->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition());
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633 double energyLoss_ = 0.;
0634 unsigned int monoId = 0;
0635 std::vector<double> energyLossM;
0636 std::vector<double> energyLossS;
0637 std::vector<PSimHit> assSimHits = hitAssociator.associateHit(**itHit);
0638 if (assSimHits.empty())
0639 continue;
0640 PSimHit shit;
0641 std::vector<unsigned int> trackIds;
0642 energyLossS.clear();
0643 energyLossM.clear();
0644
0645 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0646 if (outlier)
0647 energyLoss->Fill(m->energyLoss());
0648 unsigned int tId = m->trackId();
0649 if (find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
0650 trackIds.push_back(tId);
0651 LogTrace("TestOutliers") << "id=" << tId;
0652 if (m->energyLoss() > energyLoss_) {
0653 shit = *m;
0654 energyLoss_ = m->energyLoss();
0655 }
0656 if (hittypeval == 3) {
0657 if (monoId == 0)
0658 monoId = m->detUnitId();
0659 if (monoId == m->detUnitId()) {
0660 energyLossM.push_back(m->energyLoss());
0661 } else {
0662 energyLossS.push_back(m->energyLoss());
0663 }
0664
0665 } else {
0666 energyLossM.push_back(m->energyLoss());
0667 }
0668 }
0669 unsigned int nIds = trackIds.size();
0670
0671 if (outlier) {
0672 goodbadhits->Fill(goodhit);
0673 posxy->Fill(fabs(gpos.x()), fabs(gpos.y()));
0674 poszr->Fill(fabs(gpos.z()), sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y()));
0675 process->Fill(shit.processType());
0676 energyLossMax->Fill(energyLoss_);
0677 }
0678
0679
0680 bool shared = true;
0681 bool ioniOnly = true;
0682 unsigned int idc = 0;
0683 for (size_t jj = 0; jj < tpids.size(); jj++) {
0684 idc += std::count(trackIds.begin(), trackIds.end(), tpids[jj]);
0685 }
0686 if (idc == trackIds.size()) {
0687 shared = false;
0688 }
0689 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin() + 1; m < assSimHits.end(); m++) {
0690 if ((m->processType() != 7 && m->processType() != 8 && m->processType() != 9) &&
0691 abs(m->particleType()) != 11) {
0692 ioniOnly = false;
0693 break;
0694 }
0695 }
0696 if (ioniOnly && !shared) {
0697 LogTrace("TestOutliers") << "delta";
0698 }
0699
0700 if (goodhit) {
0701 if (outlier) {
0702 goodprocess->Fill(shit.processType());
0703 if (clustersize >= 4) {
0704 goodhittype_clgteq4->Fill(hittypeval);
0705 goodlayer_clgteq4->Fill(layerval);
0706 } else {
0707 goodhittype_cllt4->Fill(hittypeval);
0708 goodlayer_cllt4->Fill(layerval);
0709 }
0710
0711 if (hittypeval == 1 && clustersize >= 4)
0712 goodpixgteq4_simvecsize->Fill(assSimHits.size());
0713 if (hittypeval == 1 && clustersize < 4)
0714 goodpixlt4_simvecsize->Fill(assSimHits.size());
0715 if (hittypeval == 2 && clustersize >= 4)
0716 goodst1gteq4_simvecsize->Fill(assSimHits.size());
0717 if (hittypeval == 2 && clustersize < 4)
0718 goodst1lt4_simvecsize->Fill(assSimHits.size());
0719 if (hittypeval == 3 && clustersize >= 4)
0720 goodst2gteq4_simvecsize->Fill(assSimHits.size());
0721 if (hittypeval == 3 && clustersize < 4)
0722 goodst2lt4_simvecsize->Fill(assSimHits.size());
0723 if (hittypeval == 4 && clustersize >= 4)
0724 goodprjgteq4_simvecsize->Fill(assSimHits.size());
0725 if (hittypeval == 4 && clustersize < 4)
0726 goodprjlt4_simvecsize->Fill(assSimHits.size());
0727
0728
0729 if (hittypeval == 1)
0730 goodpix_clustersize->Fill(clustersize);
0731 if (hittypeval == 2)
0732 goodst1_clustersize->Fill(clustersize);
0733 if (hittypeval == 3)
0734 goodst2_clustersize->Fill(clustersize);
0735 if (hittypeval == 4)
0736 goodprj_clustersize->Fill(clustersize);
0737 if (hittypeval == 1)
0738 goodpix_simvecsize->Fill(assSimHits.size());
0739 if (hittypeval == 2)
0740 goodst1_simvecsize->Fill(assSimHits.size());
0741 if (hittypeval == 3)
0742 goodst2_simvecsize->Fill(assSimHits.size());
0743 if (hittypeval == 4)
0744 goodprj_simvecsize->Fill(assSimHits.size());
0745
0746
0747 nOfTrackIds->Fill(nIds);
0748 if (hittypeval != 3) {
0749 if (energyLossM.size() > 1) {
0750 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0751 energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
0752 }
0753 } else {
0754
0755 if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
0756 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0757 energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
0758 } else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
0759 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0760 energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0761 } else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
0762 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0763 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0764 energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
0765 ? energyLossRatio->Fill(energyLossM[1] / energyLossM[0])
0766 : energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0767 }
0768 }
0769
0770 LogTrace("TestOutliers") << "before merged";
0771 const SiStripMatchedRecHit2D *tmp = dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit);
0772 LogTrace("TestOutliers") << "tmp=" << tmp;
0773 LogTrace("TestOutliers") << "assSimHits.size()=" << assSimHits.size();
0774 if ((assSimHits.size() > 1 && tmp == nullptr) || (assSimHits.size() > 2 && tmp != nullptr)) {
0775
0776
0777 mergedlayer->Fill(layerval);
0778 mergedcluster->Fill(clustersize);
0779 mergedhittype->Fill(hittypeval);
0780
0781 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0782 unsigned int tId = m->trackId();
0783 LogTrace("TestOutliers") << "component with id=" << tId << " eLoss=" << m->energyLoss()
0784 << " proc=" << m->processType() << " part=" << m->particleType();
0785 if (find(tpids.begin(), tpids.end(), tId) == tpids.end())
0786 continue;
0787 if (m->processType() == 2) {
0788
0789
0790 AlgebraicSymMatrix33 ger = ErrorFrameTransformer()
0791 .transform((*itHit)->localPositionError(),
0792 theG->idToDet((*itHit)->geographicalId())->surface())
0793 .matrix();
0794
0795 GlobalPoint gps =
0796 theG->idToDet((*itHit)->geographicalId())->surface().toGlobal(m->localPosition());
0797
0798 ROOT::Math::SVector<double, 3> delta;
0799 delta[0] = gpos.x() - gps.x();
0800 delta[1] = gpos.y() - gps.y();
0801 delta[2] = gpos.z() - gps.z();
0802 LogTrace("TestOutliers") << delta << " " << ger;
0803 double mpull = sqrt(delta[0] * delta[0] / ger[0][0] + delta[1] * delta[1] / ger[1][1] +
0804 delta[2] * delta[2] / ger[2][2]);
0805 LogTrace("TestOutliers") << "hit pull=" << mpull;
0806 mergedPull->Fill(mpull);
0807 }
0808 }
0809 LogTrace("TestOutliers") << "end merged";
0810 } else {
0811
0812 goodlayer->Fill(layerval);
0813 goodcluster->Fill(clustersize);
0814 goodhittype->Fill(hittypeval);
0815 }
0816 }
0817
0818 const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
0819 if ((hittypeval != 3 && assSimHits.size() < 2) || (hittypeval == 3 && assSimHits.size() < 3)) {
0820 if (outlier) {
0821 goodhittype_simvecsmall->Fill(hittypeval);
0822 goodlayer_simvecsmall->Fill(layerval);
0823 goodbadmerged->Fill(1);
0824 if (pix) {
0825 probXgood->Fill(pix->probabilityX());
0826 probYgood->Fill(pix->probabilityY());
0827 }
0828 LogTrace("TestOutliers") << "out good";
0829 } else if (lost) {
0830 goodbadmergedLost->Fill(1);
0831 LogTrace("TestOutliers") << "lost good";
0832 } else if (gained) {
0833 goodbadmergedGained->Fill(1);
0834 LogTrace("TestOutliers") << "gained good";
0835 }
0836 LogTrace("TestOutliers") << "good";
0837 } else {
0838 if (outlier) {
0839 goodhittype_simvecbig->Fill(hittypeval);
0840 goodlayer_simvecbig->Fill(layerval);
0841 if (ioniOnly && !shared) {
0842 goodbadmerged->Fill(3);
0843 if (pix) {
0844 probXdelta->Fill(pix->probabilityX());
0845 probYdelta->Fill(pix->probabilityY());
0846 }
0847 } else if (!ioniOnly && !shared) {
0848 goodbadmerged->Fill(4);
0849 if (pix) {
0850 probXnoshare->Fill(pix->probabilityX());
0851 probYnoshare->Fill(pix->probabilityY());
0852 }
0853 } else {
0854 goodbadmerged->Fill(5);
0855 if (pix) {
0856 probXshared->Fill(pix->probabilityX());
0857 probYshared->Fill(pix->probabilityY());
0858 }
0859 }
0860 LogTrace("TestOutliers") << "out merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0861 } else if (lost) {
0862 if (ioniOnly && !shared)
0863 goodbadmergedLost->Fill(3);
0864 else if (!ioniOnly && !shared)
0865 goodbadmergedLost->Fill(4);
0866 else
0867 goodbadmergedLost->Fill(5);
0868 LogTrace("TestOutliers") << "lost merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0869 } else if (gained) {
0870 if (ioniOnly && !shared)
0871 goodbadmergedGained->Fill(3);
0872 else if (!ioniOnly && !shared)
0873 goodbadmergedGained->Fill(4);
0874 else
0875 goodbadmergedGained->Fill(5);
0876 LogTrace("TestOutliers") << "gained merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0877 }
0878 LogTrace("TestOutliers") << "merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0879 }
0880 }
0881 else {
0882
0883 if (outlier) {
0884 badcluster->Fill(clustersize);
0885 badhittype->Fill(hittypeval);
0886 badlayer->Fill(layerval);
0887 badprocess->Fill(shit.processType());
0888 goodbadmerged->Fill(2);
0889 const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
0890 if (pix) {
0891 probXbad->Fill(pix->probabilityX());
0892 probYbad->Fill(pix->probabilityY());
0893 }
0894 LogTrace("TestOutliers") << "out bad";
0895 } else if (lost) {
0896 goodbadmergedLost->Fill(2);
0897 LogTrace("TestOutliers") << "lost bad";
0898 } else if (gained) {
0899 goodbadmergedGained->Fill(2);
0900 LogTrace("TestOutliers") << "gained bad";
0901 }
0902 LogTrace("TestOutliers") << "bad";
0903 }
0904 }
0905 }
0906 }
0907
0908 else if (false) {
0909 LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
0910 tracks->Fill(1);
0911 LogTrace("TestOutliers")
0912 << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
0913 << sqrt(tpr->momentum().perp2())
0914
0915 << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
0916 << trackOut->numberOfValidHits()
0917
0918 << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0919 LogTrace("TestOutliers") << "track with gained hits";
0920 gainedhits2->Fill(trackOut->numberOfValidHits() - trackOld->numberOfValidHits());
0921 } else {
0922 LogTrace("TestOutliers") << "no outliers for track with #hits=" << trackOut->numberOfValidHits();
0923 }
0924 }
0925 LogTrace("TestOutliers") << "end track old #" << j;
0926 }
0927
0928 for (unsigned int k = 0; k < tracksOut->size(); ++k) {
0929 if (find(outused.begin(), outused.end(), k) == outused.end()) {
0930 edm::RefToBase<reco::Track> trackOut(tracksOut, k);
0931 bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
0932 if (outAssoc) {
0933 hitsPerTrackGained->Fill(trackOut->numberOfValidHits());
0934 LogTrace("TestOutliers") << "gained track out id=" << k << " #hits==" << trackOut->numberOfValidHits();
0935 }
0936 }
0937 }
0938 }
0939
0940
0941 void TestOutliers::beginRun(edm::Run const &run, const edm::EventSetup &es) {
0942 const bool oldAddDir = TH1::AddDirectoryStatus();
0943 TH1::AddDirectory(true);
0944 file = new TFile(out.c_str(), "recreate");
0945 histoPtOut = new TH1F("histoPtOut", "histoPtOut", 100, -10, 10);
0946 histoPtOld = new TH1F("histoPtOld", "histoPtOld", 100, -10, 10);
0947 histoQoverpOut = new TH1F("histoQoverpOut", "histoQoverpOut", 250, -25, 25);
0948 histoQoverpOld = new TH1F("histoQoverpOld", "histoQoverpOld", 250, -25, 25);
0949 histoPhiOut = new TH1F("histoPhiOut", "histoPhiOut", 250, -25, 25);
0950 histoPhiOld = new TH1F("histoPhiOld", "histoPhiOld", 250, -25, 25);
0951 histoD0Out = new TH1F("histoD0Out", "histoD0Out", 250, -25, 25);
0952 histoD0Old = new TH1F("histoD0Old", "histoD0Old", 250, -25, 25);
0953 histoDzOut = new TH1F("histoDzOut", "histoDzOut", 250, -25, 25);
0954 histoDzOld = new TH1F("histoDzOld", "histoDzOld", 250, -25, 25);
0955 histoLambdaOut = new TH1F("histoLambdaOut", "histoLambdaOut", 250, -25, 25);
0956 histoLambdaOld = new TH1F("histoLambdaOld", "histoLambdaOld", 250, -25, 25);
0957 deltahits = new TH1F("deltahits", "deltahits", 80, -40, 40);
0958 deltahitsOK = new TH1F("deltahitsOK", "deltahitsOK", 20, 0, 20);
0959 deltahitsNO = new TH1F("deltahitsNO", "deltahitsNO", 20, 0, 20);
0960 okcutsOut = new TH1F("okcutsOut", "okcutsOut", 2, -0.5, 1.5);
0961 okcutsOld = new TH1F("okcutsOld", "okcutsOld", 2, -0.5, 1.5);
0962 tracks = new TH1F("tracks_", "tracks_", 2, -0.5, 1.5);
0963 goodbadhits = new TH1F("goodbadhits", "goodbadhits", 2, -0.5, 1.5);
0964 process = new TH1F("process", "process", 20, -0.5, 19.5);
0965 goodcluster = new TH1F("goodcluster", "goodcluster", 40, -0.5, 39.5);
0966 goodprocess = new TH1F("goodprocess", "goodprocess", 20, -0.5, 19.5);
0967 badcluster = new TH1F("badcluster", "badcluster", 40, -0.5, 39.5);
0968 badprocess = new TH1F("badprocess", "badprocess", 20, -0.5, 19.5);
0969 goodhittype = new TH1F("goodhittype", "goodhittype", 5, -0.5, 4.5);
0970 goodlayer = new TH1F("goodlayer", "goodlayer", 70, -0.5, 69.5);
0971 goodhittype_clgteq4 = new TH1F("goodhittype_clgteq4", "goodhittype_clgteq4", 5, -0.5, 4.5);
0972 goodlayer_clgteq4 = new TH1F("goodlayer_clgteq4", "goodlayer_clgteq4", 70, -0.5, 69.5);
0973 goodhittype_cllt4 = new TH1F("goodhittype_cllt4", "goodhittype_cllt4", 5, -0.5, 4.5);
0974 goodlayer_cllt4 = new TH1F("goodlayer_cllt4", "goodlayer_cllt4", 70, -0.5, 69.5);
0975 badhittype = new TH1F("badhittype", "badhittype", 5, -0.5, 4.5);
0976 badlayer = new TH1F("badlayer", "badlayer", 70, -0.5, 69.5);
0977 posxy = new TH2F("posxy", "posxy", 1200, 0, 120, 1200, 0, 120);
0978 poszr = new TH2F("poszr", "poszr", 3000, 0, 300, 1200, 0, 120);
0979 goodpixgteq4_simvecsize = new TH1F("goodpixgteq4_simvecsize", "goodpixgteq4_simvecsize", 40, -0.5, 39.5);
0980 goodpixlt4_simvecsize = new TH1F("goodpixlt4_simvecsize", "goodpixlt4_simvecsize", 40, -0.5, 39.5);
0981 goodst1gteq4_simvecsize = new TH1F("goodst1gteq4_simvecsize", "goodst1gteq4_simvecsize", 40, -0.5, 39.5);
0982 goodst1lt4_simvecsize = new TH1F("goodst1lt4_simvecsize", "goodst1lt4_simvecsize", 40, -0.5, 39.5);
0983 goodst2gteq4_simvecsize = new TH1F("goodst2gteq4_simvecsize", "goodst2gteq4_simvecsize", 40, -0.5, 39.5);
0984 goodst2lt4_simvecsize = new TH1F("goodst2lt4_simvecsize", "goodst2lt4_simvecsize", 40, -0.5, 39.5);
0985 goodprjgteq4_simvecsize = new TH1F("goodprjgteq4_simvecsize", "goodprjgteq4_simvecsize", 40, -0.5, 39.5);
0986 goodprjlt4_simvecsize = new TH1F("goodprjlt4_simvecsize", "goodprjlt4_simvecsize", 40, -0.5, 39.5);
0987 goodpix_clustersize = new TH1F("goodpix_clustersize", "goodpix_clustersize", 40, -0.5, 39.5);
0988 goodst1_clustersize = new TH1F("goodst1_clustersize", "goodst1_clustersize", 40, -0.5, 39.5);
0989 goodst2_clustersize = new TH1F("goodst2_clustersize", "goodst2_clustersize", 40, -0.5, 39.5);
0990 goodprj_clustersize = new TH1F("goodprj_clustersize", "goodprj_clustersize", 40, -0.5, 39.5);
0991 goodpix_simvecsize = new TH1F("goodpix_simvecsize", "goodpix_simvecsize", 40, -0.5, 39.5);
0992 goodst1_simvecsize = new TH1F("goodst1_simvecsize", "goodst1_simvecsize", 40, -0.5, 39.5);
0993 goodst2_simvecsize = new TH1F("goodst2_simvecsize", "goodst2_simvecsize", 40, -0.5, 39.5);
0994 goodprj_simvecsize = new TH1F("goodprj_simvecsize", "goodprj_simvecsize", 40, -0.5, 39.5);
0995 goodhittype_simvecsmall = new TH1F("goodhittype_simvecsmall", "goodhittype_simvecsmall", 5, -0.5, 4.5);
0996 goodlayer_simvecsmall = new TH1F("goodlayer_simvecsmall", "goodlayer_simvecsmall", 70, -0.5, 69.5);
0997 goodhittype_simvecbig = new TH1F("goodhittype_simvecbig", "goodhittype_simvecbig", 5, -0.5, 4.5);
0998 goodlayer_simvecbig = new TH1F("goodlayer_simvecbig", "goodlayer_simvecbig", 70, -0.5, 69.5);
0999 goodbadmerged = new TH1F("goodbadmerged", "goodbadmerged", 5, 0.5, 5.5);
1000 goodbadmergedLost = new TH1F("goodbadmergedLost", "goodbadmergedLost", 5, 0.5, 5.5);
1001 goodbadmergedGained = new TH1F("goodbadmergedGained", "goodbadmergedGained", 5, 0.5, 5.5);
1002 energyLoss = new TH1F("energyLoss", "energyLoss", 1000, 0, 0.1);
1003 energyLossMax = new TH1F("energyLossMax", "energyLossMax", 1000, 0, 0.1);
1004 energyLossRatio = new TH1F("energyLossRatio", "energyLossRatio", 100, 0, 1);
1005 nOfTrackIds = new TH1F("nOfTrackIds", "nOfTrackIds", 10, 0, 10);
1006 mergedPull = new TH1F("mergedPull", "mergedPull", 100, 0, 10);
1007 mergedlayer = new TH1F("mergedlayer", "mergedlayer", 70, -0.5, 69.5);
1008 mergedhittype = new TH1F("mergedhittype", "mergedhittype", 5, -0.5, 4.5);
1009 mergedcluster = new TH1F("mergedcluster", "mergedcluster", 40, -0.5, 39.5);
1010 deltahitsAssocGained = new TH1F("deltahitsAssocGained", "deltahitsAssocGained", 80, -40, 40);
1011 deltahitsAssocLost = new TH1F("deltahitsAssocLost", "deltahitsAssocLost", 80, -40, 40);
1012 hitsPerTrackLost = new TH1F("hitsPerTrackLost", "hitsPerTrackLost", 40, -0.5, 39.5);
1013 hitsPerTrackAssocLost = new TH1F("hitsPerTrackAssocLost", "hitsPerTrackAssocLost", 40, -0.5, 39.5);
1014 hitsPerTrackGained = new TH1F("hitsPerTrackGained", "hitsPerTrackGained", 40, -0.5, 39.5);
1015 hitsPerTrackAssocGained = new TH1F("hitsPerTrackAssocGained", "hitsPerTrackAssocGained", 40, -0.5, 39.5);
1016 sizeOut = new TH1F("sizeOut", "sizeOut", 900, -0.5, 899.5);
1017 sizeOld = new TH1F("sizeOld", "sizeOld", 900, -0.5, 899.5);
1018 sizeOutT = new TH1F("sizeOutT", "sizeOutT", 900, -0.5, 899.5);
1019 sizeOldT = new TH1F("sizeOldT", "sizeOldT", 900, -0.5, 899.5);
1020 countOutA = new TH1F("countOutA", "countOutA", 2, 0, 2);
1021 countOutT = new TH1F("countOutT", "countOutT", 2, 0, 2);
1022 countOldT = new TH1F("countOldT", "countOldT", 2, 0, 2);
1023 gainedhits = new TH1F("gainedhits", "gainedhits", 2, 0, 2);
1024 gainedhits2 = new TH1F("gainedhits2", "gainedhits2", 30, -0.5, 29.5);
1025 probXgood = new TH1F("probXgood", "probXgood", 110, 0, 1.1);
1026 probXbad = new TH1F("probXbad", "probXbad", 110, 0, 1.1);
1027 probXdelta = new TH1F("probXdelta", "probXdelta", 110, 0, 1.1);
1028 probXshared = new TH1F("probXshared", "probXshared", 110, 0, 1.1);
1029 probXnoshare = new TH1F("probXnoshare", "probXnoshare", 110, 0, 1.1);
1030 probYgood = new TH1F("probYgood", "probYgood", 110, 0, 1.1);
1031 probYbad = new TH1F("probYbad", "probYbad", 110, 0, 1.1);
1032 probYdelta = new TH1F("probYdelta", "probYdelta", 110, 0, 1.1);
1033 probYshared = new TH1F("probYshared", "probYshared", 110, 0, 1.1);
1034 probYnoshare = new TH1F("probYnoshare", "probYnoshare", 110, 0, 1.1);
1035 TH1::AddDirectory(oldAddDir);
1036 }
1037
1038 void TestOutliers::endJob() {
1039 LogTrace("TestOutliers") << "TestOutliers::endJob";
1040 file->Write();
1041 LogTrace("TestOutliers") << "outfile written";
1042 file->Close();
1043 LogTrace("TestOutliers") << "oufile closed";
1044 LogTrace("TestOutliers") << "exiting TestOutliers::endJob";
1045 }
1046
1047
1048 DEFINE_FWK_MODULE(TestOutliers);
