Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:06:37 

0001 // -*- C++ -*-
0002 //
0003 // Package:    DPGAnalysis
0004 // Class:      TrackerDpgAnalysis
0005 //
0006 /**\class TrackerDpgAnalysis TrackerDpgAnalysis.cc DPGAnalysis/SiStripTools/plugins/TrackerDpgAnalysis.cc
0007 
0008  Description: analysis of the clusters and digis in the tracker
0009 
0010  Implementation:
0011       analysis of the clusters and digis in the tracker
0012 */
0013 //
0014 // Original Author:  Christophe DELAERE
0015 //         Created:  Tue Sep 23 02:11:44 CEST 2008
0016 //         Revised:  Thu Nov 26 10:00:00 CEST 2009
0017 // part of the code was inspired by http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/YGao/LhcTrackAnalyzer/
0018 // part of the code was inspired by
0019 // other inputs from Andrea Giammanco, Gaelle Boudoul, Andrea Venturi, Steven Lowette, Gavril Giurgiu
0020 // $Id: TrackerDpgAnalysis.cc,v 1.14 2013/02/27 19:49:47 wmtan Exp $
0021 //
0022 //
0023 
0024 // system include files
0025 #include <memory>
0026 #include <iostream>
0027 #include <limits>
0028 #include <utility>
0029 #include <vector>
0030 #include <algorithm>
0031 #include <functional>
0032 #include <cstring>
0033 #include <sstream>
0034 #include <fstream>
0035 
0036 // root include files
0037 #include "TTree.h"
0038 #include "TFile.h"
0039 
0040 // user include files
0041 #include "FWCore/Framework/interface/Frameworkfwd.h"
0042 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0043 #include "FWCore/Framework/interface/Event.h"
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 #include "FWCore/Utilities/interface/InputTag.h"
0046 #include "FWCore/Utilities/interface/transform.h"
0047 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0050 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0051 #include <CondFormats/SiStripObjects/interface/FedChannelConnection.h>
0052 #include <CondFormats/SiStripObjects/interface/SiStripFedCabling.h>
0053 #include <CondFormats/DataRecord/interface/SiStripFedCablingRcd.h>
0054 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0055 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0056 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0057 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0058 #include <Geometry/CommonTopologies/interface/Topology.h>
0059 #include <Geometry/CommonTopologies/interface/StripTopology.h>
0060 #include <Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h>
0061 #include <Geometry/CommonTopologies/interface/PixelTopology.h>
0062 #include "DataFormats/Common/interface/Ref.h"
0063 #include "DataFormats/DetId/interface/DetId.h"
0064 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0065 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0066 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0067 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0068 #include "DataFormats/SiStripCommon/interface/SiStripEventSummary.h"
0069 #include <DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h>
0070 #include <DataFormats/SiStripDetId/interface/SiStripDetId.h>
0071 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0072 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0073 #include <DataFormats/SiPixelCluster/interface/SiPixelCluster.h>
0074 #include <DataFormats/TrackReco/interface/TrackFwd.h>
0075 #include <DataFormats/TrackReco/interface/Track.h>
0076 #include "DataFormats/TrackReco/interface/DeDxData.h"
0077 #include <DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h>
0078 #include <DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h>
0079 #include <DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h>
0080 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0081 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0082 #include <DataFormats/L1GlobalTrigger/interface/L1GtFdlWord.h>
0083 #include <DataFormats/Common/interface/TriggerResults.h>
0084 #include "DataFormats/VertexReco/interface/Vertex.h"
0085 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0086 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0087 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0088 #include <MagneticField/Engine/interface/MagneticField.h>
0089 #include <MagneticField/Records/interface/IdealMagneticFieldRecord.h>
0090 #include <RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h>
0091 #include <RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h>
0092 #include <RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit1D.h>
0093 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0094 #include <TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h>
0095 #include <TrackingTools/PatternTools/interface/Trajectory.h>
0096 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0097 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0098 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
0099 
0100 // topology
0101 #include "DPGAnalysis/SiStripTools/interface/EventShape.h"
0102 
0103 //
0104 // class decleration
0105 //
0106 typedef math::XYZPoint Point;
0107 
0108 class TrackerDpgAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0109 public:
0110   explicit TrackerDpgAnalysis(const edm::ParameterSet&);
0111   ~TrackerDpgAnalysis() override;
0112 
0113 protected:
0114   std::vector<double> onTrackAngles(edm::Handle<edmNew::DetSetVector<SiStripCluster> >&,
0115                                     const std::vector<Trajectory>&);
0116   void insertMeasurement(std::multimap<const uint32_t, std::pair<LocalPoint, double> >&,
0117                          const TransientTrackingRecHit*,
0118                          double);
0119   std::vector<int> onTrack(edm::Handle<edmNew::DetSetVector<SiStripCluster> >&, const reco::TrackCollection&, uint32_t);
0120   void insertMeasurement(std::multimap<const uint32_t, std::pair<int, int> >&, const TrackingRecHit*, int);
0121   std::map<size_t, int> inVertex(const reco::TrackCollection&, const reco::VertexCollection&, uint32_t);
0122   std::vector<std::pair<double, double> > onTrackAngles(edm::Handle<edmNew::DetSetVector<SiPixelCluster> >&,
0123                                                         const std::vector<Trajectory>&);
0124   void insertMeasurement(std::multimap<const uint32_t, std::pair<LocalPoint, std::pair<double, double> > >&,
0125                          const TransientTrackingRecHit*,
0126                          double,
0127                          double);
0128   std::vector<int> onTrack(edm::Handle<edmNew::DetSetVector<SiPixelCluster> >&, const reco::TrackCollection&, uint32_t);
0129   void insertMeasurement(std::multimap<const uint32_t, std::pair<std::pair<float, float>, int> >&,
0130                          const TrackingRecHit*,
0131                          int);
0132   std::string toStringName(uint32_t, const TrackerTopology&);
0133   std::string toStringId(uint32_t);
0134   double sumPtSquared(const reco::Vertex&);
0135   float delay(const SiStripEventSummary&);
0136   std::map<uint32_t, float> delay(const std::vector<std::string>&);
0137 
0138 private:
0139   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0140   void endRun(const edm::Run&, const edm::EventSetup&) override {}
0141   void analyze(const edm::Event&, const edm::EventSetup&) override;
0142   void endJob() override;
0143 
0144   // ----------member data ---------------------------
0145   static const int nMaxPVs_ = 50;
0146   SiStripClusterInfo siStripClusterInfo_;
0147   edm::EDGetTokenT<SiStripEventSummary> summaryToken_;
0148   edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > clusterToken_;
0149   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > pixelclusterToken_;
0150   edm::EDGetTokenT<edm::ValueMap<reco::DeDxData> > dedx1Token_;
0151   edm::EDGetTokenT<edm::ValueMap<reco::DeDxData> > dedx2Token_;
0152   edm::EDGetTokenT<edm::ValueMap<reco::DeDxData> > dedx3Token_;
0153   edm::EDGetTokenT<reco::VertexCollection> pixelVertexToken_;
0154   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0155   edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0156   edm::EDGetTokenT<L1GlobalTriggerReadoutRecord> L1Token_;
0157   edm::InputTag HLTTag_;
0158   edm::EDGetTokenT<edm::TriggerResults> HLTToken_;
0159   std::vector<edm::InputTag> trackLabels_;
0160   std::vector<edm::EDGetTokenT<reco::TrackCollection> > trackTokens_;
0161   std::vector<edm::EDGetTokenT<std::vector<Trajectory> > > trajectoryTokens_;
0162   std::vector<edm::EDGetTokenT<TrajTrackAssociationCollection> > trajTrackAssoTokens_;
0163   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0164   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0165   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0166   edm::ESGetToken<SiStripFedCabling, SiStripFedCablingRcd> fedCablingToken_;
0167   const TrackerGeometry* tracker_;
0168   std::multimap<const uint32_t, const FedChannelConnection*> connections_;
0169   bool functionality_offtrackClusters_, functionality_ontrackClusters_, functionality_pixclusters_,
0170       functionality_pixvertices_, functionality_missingHits_, functionality_tracks_, functionality_vertices_,
0171       functionality_events_;
0172   TTree* clusters_;
0173   TTree* pixclusters_;
0174   std::vector<TTree*> tracks_;
0175   std::vector<TTree*> missingHits_;
0176   TTree* vertices_;
0177   TTree* pixelVertices_;
0178   TTree* event_;
0179   TTree* psumap_;
0180   TTree* readoutmap_;
0181   bool onTrack_;
0182   uint32_t vertexid_;
0183   edm::EventNumber_t eventid_;
0184   uint32_t runid_;
0185   uint32_t globalvertexid_;
0186   uint32_t *globaltrackid_, *trackid_;
0187   float globalX_, globalY_, globalZ_;
0188   float measX_, measY_, errorX_, errorY_;
0189   float angle_, maxCharge_;
0190   float clCorrectedCharge_, clCorrectedSignalOverNoise_;
0191   float clNormalizedCharge_, clNormalizedNoise_, clSignalOverNoise_;
0192   float clBareNoise_, clBareCharge_;
0193   float clWidth_, clPosition_, thickness_, stripLength_, distance_;
0194   float eta_, phi_, chi2_;
0195   float dedx1_, dedx2_, dedx3_;
0196   uint32_t detid_, dcuId_, type_;
0197   uint16_t fecCrate_, fecSlot_, fecRing_, ccuAdd_, ccuChan_, lldChannel_, fedId_, fedCh_, fiberLength_;
0198   uint32_t nclusters_, npixClusters_, nclustersOntrack_, npixClustersOntrack_, dedxNoM_, quality_, foundhits_,
0199       lostHits_, ndof_;
0200   uint32_t *ntracks_, *ntrajs_;
0201   float* lowPixelProbabilityFraction_;
0202   uint32_t nVertices_, nPixelVertices_, nLayers_, foundhitsStrips_, foundhitsPixels_, losthitsStrips_, losthitsPixels_;
0203   uint32_t nTracks_pvtx_;
0204   uint32_t clSize_, clSizeX_, clSizeY_;
0205   float fBz_, clPositionX_, clPositionY_, alpha_, beta_, chargeCorr_;
0206   float recx_pvtx_, recy_pvtx_, recz_pvtx_, recx_err_pvtx_, recy_err_pvtx_, recz_err_pvtx_, sumptsq_pvtx_;
0207   float pterr_, etaerr_, phierr_;
0208   float dz_, dzerr_, dzCorr_, dxy_, dxyerr_, dxyCorr_;
0209   float qoverp_, xPCA_, yPCA_, zPCA_, trkWeightpvtx_;
0210   bool isValid_pvtx_, isFake_pvtx_;
0211   float charge_, p_, pt_;
0212   float bsX0_, bsY0_, bsZ0_, bsSigmaZ_, bsDxdz_, bsDydz_;
0213   float thrustValue_, thrustX_, thrustY_, thrustZ_, sphericity_, planarity_, aplanarity_, delay_;
0214   bool L1DecisionBits_[192], L1TechnicalBits_[64], HLTDecisionBits_[256];
0215   uint32_t orbit_, orbitL1_, bx_, store_, time_;
0216   uint16_t lumiSegment_, physicsDeclared_;
0217   char *moduleName_, *moduleId_, *PSUname_;
0218   std::string cablingFileName_;
0219   std::vector<std::string> delayFileNames_;
0220   edm::ParameterSet pset_;
0221   std::vector<std::string> hlNames_;  // name of each HLT algorithm
0222   HLTConfigProvider hltConfig_;       // to get configuration for L1s/Pre
0223 };
0224 
0225 //
0226 // constructors and destructor
0227 //
0228 TrackerDpgAnalysis::TrackerDpgAnalysis(const edm::ParameterSet& iConfig)
0229     : siStripClusterInfo_(consumesCollector(), std::string("")),
0230       magFieldToken_(esConsumes()),
0231       tTopoToken_(esConsumes<edm::Transition::BeginRun>()),
0232       tkGeomToken_(esConsumes<edm::Transition::BeginRun>()),
0233       fedCablingToken_(esConsumes<edm::Transition::BeginRun>()),
0234       hltConfig_() {
0235   // members
0236   moduleName_ = new char[256];
0237   moduleId_ = new char[256];
0238   PSUname_ = new char[256];
0239   pset_ = iConfig;
0240 
0241   usesResource(TFileService::kSharedResource);
0242 
0243   // enable/disable functionalities
0244   functionality_offtrackClusters_ = iConfig.getUntrackedParameter<bool>("keepOfftrackClusters", true);
0245   functionality_ontrackClusters_ = iConfig.getUntrackedParameter<bool>("keepOntrackClusters", true);
0246   functionality_pixclusters_ = iConfig.getUntrackedParameter<bool>("keepPixelClusters", true);
0247   functionality_pixvertices_ = iConfig.getUntrackedParameter<bool>("keepPixelVertices", true);
0248   functionality_missingHits_ = iConfig.getUntrackedParameter<bool>("keepMissingHits", true);
0249   functionality_tracks_ = iConfig.getUntrackedParameter<bool>("keepTracks", true);
0250   functionality_vertices_ = iConfig.getUntrackedParameter<bool>("keepVertices", true);
0251   functionality_events_ = iConfig.getUntrackedParameter<bool>("keepEvents", true);
0252 
0253   // parameters
0254   summaryToken_ = consumes<SiStripEventSummary>(edm::InputTag("siStripDigis"));
0255   clusterToken_ = consumes<edmNew::DetSetVector<SiStripCluster> >(iConfig.getParameter<edm::InputTag>("ClustersLabel"));
0256   pixelclusterToken_ =
0257       consumes<edmNew::DetSetVector<SiPixelCluster> >(iConfig.getParameter<edm::InputTag>("PixelClustersLabel"));
0258   trackLabels_ = iConfig.getParameter<std::vector<edm::InputTag> >("TracksLabel");
0259   trackTokens_ = edm::vector_transform(
0260       trackLabels_, [this](edm::InputTag const& tag) { return consumes<reco::TrackCollection>(tag); });
0261   trajectoryTokens_ = edm::vector_transform(
0262       trackLabels_, [this](edm::InputTag const& tag) { return consumes<std::vector<Trajectory> >(tag); });
0263   trajTrackAssoTokens_ = edm::vector_transform(
0264       trackLabels_, [this](edm::InputTag const& tag) { return consumes<TrajTrackAssociationCollection>(tag); });
0265   dedx1Token_ = consumes<edm::ValueMap<reco::DeDxData> >(iConfig.getParameter<edm::InputTag>("DeDx1Label"));
0266   dedx2Token_ = consumes<edm::ValueMap<reco::DeDxData> >(iConfig.getParameter<edm::InputTag>("DeDx2Label"));
0267   dedx3Token_ = consumes<edm::ValueMap<reco::DeDxData> >(iConfig.getParameter<edm::InputTag>("DeDx3Label"));
0268   vertexToken_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexLabel"));
0269   pixelVertexToken_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("pixelVertexLabel"));
0270   bsToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotLabel"));
0271   L1Token_ = consumes<L1GlobalTriggerReadoutRecord>(iConfig.getParameter<edm::InputTag>("L1Label"));
0272   HLTTag_ = iConfig.getParameter<edm::InputTag>("HLTLabel");
0273   HLTToken_ = consumes<edm::TriggerResults>(HLTTag_);
0274 
0275   // initialize arrays
0276   size_t trackSize(trackLabels_.size());
0277   ntracks_ = new uint32_t[trackSize];
0278   ntrajs_ = new uint32_t[trackSize];
0279   globaltrackid_ = new uint32_t[trackSize];
0280   trackid_ = new uint32_t[trackSize];
0281   lowPixelProbabilityFraction_ = new float[trackSize];
0282   globalvertexid_ = iConfig.getParameter<uint32_t>("InitalCounter");
0283   for (size_t i = 0; i < trackSize; ++i) {
0284     ntracks_[i] = 0;
0285     ntrajs_[i] = 0;
0286     globaltrackid_[i] = iConfig.getParameter<uint32_t>("InitalCounter");
0287     trackid_[i] = 0;
0288     lowPixelProbabilityFraction_[i] = 0;
0289   }
0290 
0291   // create output
0292   edm::Service<TFileService> fileService;
0293   TFileDirectory* dir = new TFileDirectory(fileService->mkdir("trackerDPG"));
0294 
0295   // create a TTree for clusters
0296   clusters_ = dir->make<TTree>("clusters", "cluster information");
0297   clusters_->Branch("eventid", &eventid_, "eventid/i");
0298   clusters_->Branch("runid", &runid_, "runid/i");
0299   for (size_t i = 0; i < trackSize; ++i) {
0300     char buffer1[256];
0301     char buffer2[256];
0302     sprintf(buffer1, "trackid%lu", (unsigned long)i);
0303     sprintf(buffer2, "trackid%lu/i", (unsigned long)i);
0304     clusters_->Branch(buffer1, trackid_ + i, buffer2);
0305   }
0306   clusters_->Branch("onTrack", &onTrack_, "onTrack/O");
0307   clusters_->Branch("clWidth", &clWidth_, "clWidth/F");
0308   clusters_->Branch("clPosition", &clPosition_, "clPosition/F");
0309   clusters_->Branch("clglobalX", &globalX_, "clglobalX/F");
0310   clusters_->Branch("clglobalY", &globalY_, "clglobalY/F");
0311   clusters_->Branch("clglobalZ", &globalZ_, "clglobalZ/F");
0312   clusters_->Branch("angle", &angle_, "angle/F");
0313   clusters_->Branch("thickness", &thickness_, "thickness/F");
0314   clusters_->Branch("maxCharge", &maxCharge_, "maxCharge/F");
0315   clusters_->Branch("clNormalizedCharge", &clNormalizedCharge_, "clNormalizedCharge/F");
0316   clusters_->Branch("clNormalizedNoise", &clNormalizedNoise_, "clNormalizedNoise/F");
0317   clusters_->Branch("clSignalOverNoise", &clSignalOverNoise_, "clSignalOverNoise/F");
0318   clusters_->Branch("clCorrectedCharge", &clCorrectedCharge_, "clCorrectedCharge/F");
0319   clusters_->Branch("clCorrectedSignalOverNoise", &clCorrectedSignalOverNoise_, "clCorrectedSignalOverNoise/F");
0320   clusters_->Branch("clBareCharge", &clBareCharge_, "clBareCharge/F");
0321   clusters_->Branch("clBareNoise", &clBareNoise_, "clBareNoise/F");
0322   clusters_->Branch("stripLength", &stripLength_, "stripLength/F");
0323   clusters_->Branch("detid", &detid_, "detid/i");
0324   clusters_->Branch("lldChannel", &lldChannel_, "lldChannel/s");
0325 
0326   // create a TTree for pixel clusters
0327   pixclusters_ = dir->make<TTree>("pixclusters", "pixel cluster information");
0328   pixclusters_->Branch("eventid", &eventid_, "eventid/i");
0329   pixclusters_->Branch("runid", &runid_, "runid/i");
0330   for (size_t i = 0; i < trackSize; ++i) {
0331     char buffer1[256];
0332     char buffer2[256];
0333     sprintf(buffer1, "trackid%lu", (unsigned long)i);
0334     sprintf(buffer2, "trackid%lu/i", (unsigned long)i);
0335     pixclusters_->Branch(buffer1, trackid_ + i, buffer2);
0336   }
0337   pixclusters_->Branch("onTrack", &onTrack_, "onTrack/O");
0338   pixclusters_->Branch("clPositionX", &clPositionX_, "clPositionX/F");
0339   pixclusters_->Branch("clPositionY", &clPositionY_, "clPositionY/F");
0340   pixclusters_->Branch("clSize", &clSize_, "clSize/i");
0341   pixclusters_->Branch("clSizeX", &clSizeX_, "clSizeX/i");
0342   pixclusters_->Branch("clSizeY", &clSizeY_, "clSizeY/i");
0343   pixclusters_->Branch("alpha", &alpha_, "alpha/F");
0344   pixclusters_->Branch("beta", &beta_, "beta/F");
0345   pixclusters_->Branch("charge", &charge_, "charge/F");
0346   pixclusters_->Branch("chargeCorr", &chargeCorr_, "chargeCorr/F");
0347   pixclusters_->Branch("clglobalX", &globalX_, "clglobalX/F");
0348   pixclusters_->Branch("clglobalY", &globalY_, "clglobalY/F");
0349   pixclusters_->Branch("clglobalZ", &globalZ_, "clglobalZ/F");
0350   pixclusters_->Branch("detid", &detid_, "detid/i");
0351 
0352   // create a tree for tracks
0353   for (size_t i = 0; i < trackSize; ++i) {
0354     char buffer1[256];
0355     char buffer2[256];
0356     sprintf(buffer1, "tracks%lu", (unsigned long)i);
0357     sprintf(buffer2, "track%lu information", (unsigned long)i);
0358     TTree* thetracks_ = dir->make<TTree>(buffer1, buffer2);
0359     sprintf(buffer1, "trackid%lu", (unsigned long)i);
0360     sprintf(buffer2, "trackid%lu/i", (unsigned long)i);
0361     thetracks_->Branch(buffer1, globaltrackid_ + i, buffer2);
0362     thetracks_->Branch("eventid", &eventid_, "eventid/i");
0363     thetracks_->Branch("runid", &runid_, "runid/i");
0364     thetracks_->Branch("chi2", &chi2_, "chi2/F");
0365     thetracks_->Branch("eta", &eta_, "eta/F");
0366     thetracks_->Branch("etaerr", &etaerr_, "etaerr/F");
0367     thetracks_->Branch("phi", &phi_, "phi/F");
0368     thetracks_->Branch("phierr", &phierr_, "phierr/F");
0369     thetracks_->Branch("dedx1", &dedx1_, "dedx1/F");
0370     thetracks_->Branch("dedx2", &dedx2_, "dedx2/F");
0371     thetracks_->Branch("dedx3", &dedx3_, "dedx3/F");
0372     thetracks_->Branch("dedxNoM", &dedxNoM_, "dedxNoM/i");
0373     thetracks_->Branch("charge", &charge_, "charge/F");
0374     thetracks_->Branch("quality", &quality_, "quality/i");
0375     thetracks_->Branch("foundhits", &foundhits_, "foundhits/i");
0376     thetracks_->Branch("lostHits", &lostHits_, "lostHits/i");
0377     thetracks_->Branch("foundhitsStrips", &foundhitsStrips_, "foundhitsStrips/i");
0378     thetracks_->Branch("foundhitsPixels", &foundhitsPixels_, "foundhitsPixels/i");
0379     thetracks_->Branch("losthitsStrips", &losthitsStrips_, "losthitsStrips/i");
0380     thetracks_->Branch("losthitsPixels", &losthitsPixels_, "losthitsPixels/i");
0381     thetracks_->Branch("p", &p_, "p/F");
0382     thetracks_->Branch("pt", &pt_, "pt/F");
0383     thetracks_->Branch("pterr", &pterr_, "pterr/F");
0384     thetracks_->Branch("ndof", &ndof_, "ndof/i");
0385     thetracks_->Branch("dz", &dz_, "dz/F");
0386     thetracks_->Branch("dzerr", &dzerr_, "dzerr/F");
0387     thetracks_->Branch("dzCorr", &dzCorr_, "dzCorr/F");
0388     thetracks_->Branch("dxy", &dxy_, "dxy/F");
0389     thetracks_->Branch("dxyerr", &dxyerr_, "dxyerr/F");
0390     thetracks_->Branch("dxyCorr", &dxyCorr_, "dxyCorr/F");
0391     thetracks_->Branch("qoverp", &qoverp_, "qoverp/F");
0392     thetracks_->Branch("xPCA", &xPCA_, "xPCA/F");
0393     thetracks_->Branch("yPCA", &yPCA_, "yPCA/F");
0394     thetracks_->Branch("zPCA", &zPCA_, "zPCA/F");
0395     thetracks_->Branch("nLayers", &nLayers_, "nLayers/i");
0396     thetracks_->Branch("trkWeightpvtx", &trkWeightpvtx_, "trkWeightpvtx/F");
0397     thetracks_->Branch("vertexid", &vertexid_, "vertexid/i");
0398     tracks_.push_back(thetracks_);
0399   }
0400 
0401   // create a tree for missing hits
0402   for (size_t i = 0; i < trackSize; ++i) {
0403     char buffer1[256];
0404     char buffer2[256];
0405     sprintf(buffer1, "misingHits%lu", (unsigned long)i);
0406     sprintf(buffer2, "missing hits from track collection %lu", (unsigned long)i);
0407     TTree* themissingHits_ = dir->make<TTree>(buffer1, buffer2);
0408     sprintf(buffer1, "trackid%lu", (unsigned long)i);
0409     sprintf(buffer2, "trackid%lu/i", (unsigned long)i);
0410     themissingHits_->Branch(buffer1, globaltrackid_ + i, buffer2);
0411     themissingHits_->Branch("eventid", &eventid_, "eventid/i");
0412     themissingHits_->Branch("runid", &runid_, "runid/i");
0413     themissingHits_->Branch("detid", &detid_, "detid/i");
0414     themissingHits_->Branch("type", &type_, "type/i");
0415     themissingHits_->Branch("localX", &clPositionX_, "localX/F");
0416     themissingHits_->Branch("localY", &clPositionY_, "localY/F");
0417     themissingHits_->Branch("globalX", &globalX_, "globalX/F");
0418     themissingHits_->Branch("globalY", &globalY_, "globalY/F");
0419     themissingHits_->Branch("globalZ", &globalZ_, "globalZ/F");
0420     themissingHits_->Branch("measX", &measX_, "measX/F");
0421     themissingHits_->Branch("measY", &measY_, "measY/F");
0422     themissingHits_->Branch("errorX", &errorX_, "errorX/F");
0423     themissingHits_->Branch("errorY", &errorY_, "errorY/F");
0424     missingHits_.push_back(themissingHits_);
0425   }
0426 
0427   // create a tree for the vertices
0428   vertices_ = dir->make<TTree>("vertices", "vertex information");
0429   vertices_->Branch("vertexid", &globalvertexid_, "vertexid/i");
0430   vertices_->Branch("eventid", &eventid_, "eventid/i");
0431   vertices_->Branch("runid", &runid_, "runid/i");
0432   vertices_->Branch("nTracks", &nTracks_pvtx_, "nTracks/i");
0433   vertices_->Branch("sumptsq", &sumptsq_pvtx_, "sumptsq/F");
0434   vertices_->Branch("isValid", &isValid_pvtx_, "isValid/O");
0435   vertices_->Branch("isFake", &isFake_pvtx_, "isFake/O");
0436   vertices_->Branch("recx", &recx_pvtx_, "recx/F");
0437   vertices_->Branch("recy", &recy_pvtx_, "recy/F");
0438   vertices_->Branch("recz", &recz_pvtx_, "recz/F");
0439   vertices_->Branch("recx_err", &recx_err_pvtx_, "recx_err/F");
0440   vertices_->Branch("recy_err", &recy_err_pvtx_, "recy_err/F");
0441   vertices_->Branch("recz_err", &recz_err_pvtx_, "recz_err/F");
0442 
0443   // create a tree for the vertices
0444   pixelVertices_ = dir->make<TTree>("pixelVertices", "pixel vertex information");
0445   pixelVertices_->Branch("eventid", &eventid_, "eventid/i");
0446   pixelVertices_->Branch("runid", &runid_, "runid/i");
0447   pixelVertices_->Branch("nTracks", &nTracks_pvtx_, "nTracks/i");
0448   pixelVertices_->Branch("sumptsq", &sumptsq_pvtx_, "sumptsq/F");
0449   pixelVertices_->Branch("isValid", &isValid_pvtx_, "isValid/O");
0450   pixelVertices_->Branch("isFake", &isFake_pvtx_, "isFake/O");
0451   pixelVertices_->Branch("recx", &recx_pvtx_, "recx/F");
0452   pixelVertices_->Branch("recy", &recy_pvtx_, "recy/F");
0453   pixelVertices_->Branch("recz", &recz_pvtx_, "recz/F");
0454   pixelVertices_->Branch("recx_err", &recx_err_pvtx_, "recx_err/F");
0455   pixelVertices_->Branch("recy_err", &recy_err_pvtx_, "recy_err/F");
0456   pixelVertices_->Branch("recz_err", &recz_err_pvtx_, "recz_err/F");
0457 
0458   // create a tree for the events
0459   event_ = dir->make<TTree>("events", "event information");
0460   event_->Branch("eventid", &eventid_, "eventid/i");
0461   event_->Branch("runid", &runid_, "runid/i");
0462   event_->Branch("L1DecisionBits", L1DecisionBits_, "L1DecisionBits[192]/O");
0463   event_->Branch("L1TechnicalBits", L1TechnicalBits_, "L1TechnicalBits[64]/O");
0464   event_->Branch("orbit", &orbit_, "orbit/i");
0465   event_->Branch("orbitL1", &orbitL1_, "orbitL1/i");
0466   event_->Branch("bx", &bx_, "bx/i");
0467   event_->Branch("store", &store_, "store/i");
0468   event_->Branch("time", &time_, "time/i");
0469   event_->Branch("delay", &delay_, "delay/F");
0470   event_->Branch("lumiSegment", &lumiSegment_, "lumiSegment/s");
0471   event_->Branch("physicsDeclared", &physicsDeclared_, "physicsDeclared/s");
0472   event_->Branch("HLTDecisionBits", HLTDecisionBits_, "HLTDecisionBits[256]/O");
0473   char buffer[256];
0474   sprintf(buffer, "ntracks[%lu]/i", (unsigned long)trackSize);
0475   event_->Branch("ntracks", ntracks_, buffer);
0476   sprintf(buffer, "ntrajs[%lu]/i", (unsigned long)trackSize);
0477   event_->Branch("ntrajs", ntrajs_, buffer);
0478   sprintf(buffer, "lowPixelProbabilityFraction[%lu]/F", (unsigned long)trackSize);
0479   event_->Branch("lowPixelProbabilityFraction", lowPixelProbabilityFraction_, buffer);
0480   event_->Branch("nclusters", &nclusters_, "nclusters/i");
0481   event_->Branch("npixClusters", &npixClusters_, "npixClusters/i");
0482   event_->Branch("nclustersOntrack", &nclustersOntrack_, "nclustersOntrack/i");
0483   event_->Branch("npixClustersOntrack", &npixClustersOntrack_, "npixClustersOntrack/i");
0484   event_->Branch("bsX0", &bsX0_, "bsX0/F");
0485   event_->Branch("bsY0", &bsY0_, "bsY0/F");
0486   event_->Branch("bsZ0", &bsZ0_, "bsZ0/F");
0487   event_->Branch("bsSigmaZ", &bsSigmaZ_, "bsSigmaZ/F");
0488   event_->Branch("bsDxdz", &bsDxdz_, "bsDxdz/F");
0489   event_->Branch("bsDydz", &bsDydz_, "bsDydz/F");
0490   event_->Branch("nVertices", &nVertices_, "nVertices/i");
0491   event_->Branch("thrustValue", &thrustValue_, "thrustValue/F");
0492   event_->Branch("thrustX", &thrustX_, "thrustX/F");
0493   event_->Branch("thrustY", &thrustY_, "thrustY/F");
0494   event_->Branch("thrustZ", &thrustZ_, "thrustZ/F");
0495   event_->Branch("sphericity", &sphericity_, "sphericity/F");
0496   event_->Branch("planarity", &planarity_, "planarity/F");
0497   event_->Branch("aplanarity", &aplanarity_, "aplanarity/F");
0498   event_->Branch("MagneticField", &fBz_, "MagneticField/F");
0499 
0500   // cabling
0501   cablingFileName_ = iConfig.getUntrackedParameter<std::string>("PSUFileName", "PSUmapping.csv");
0502   delayFileNames_ =
0503       iConfig.getUntrackedParameter<std::vector<std::string> >("DelayFileNames", std::vector<std::string>(0));
0504   psumap_ = dir->make<TTree>("psumap", "PSU map");
0505   psumap_->Branch("PSUname", PSUname_, "PSUname/C");
0506   psumap_->Branch("dcuId", &dcuId_, "dcuId/i");
0507   readoutmap_ = dir->make<TTree>("readoutMap", "cabling map");
0508   readoutmap_->Branch("detid", &detid_, "detid/i");
0509   readoutmap_->Branch("dcuId", &dcuId_, "dcuId/i");
0510   readoutmap_->Branch("fecCrate", &fecCrate_, "fecCrate/s");
0511   readoutmap_->Branch("fecSlot", &fecSlot_, "fecSlot/s");
0512   readoutmap_->Branch("fecRing", &fecRing_, "fecRing/s");
0513   readoutmap_->Branch("ccuAdd", &ccuAdd_, "ccuAdd/s");
0514   readoutmap_->Branch("ccuChan", &ccuChan_, "ccuChan/s");
0515   readoutmap_->Branch("lldChannel", &lldChannel_, "lldChannel/s");
0516   readoutmap_->Branch("fedId", &fedId_, "fedId/s");
0517   readoutmap_->Branch("fedCh", &fedCh_, "fedCh/s");
0518   readoutmap_->Branch("fiberLength", &fiberLength_, "fiberLength/s");
0519   readoutmap_->Branch("moduleName", moduleName_, "moduleName/C");
0520   readoutmap_->Branch("moduleId", moduleId_, "moduleId/C");
0521   readoutmap_->Branch("delay", &delay_, "delay/F");
0522   readoutmap_->Branch("globalX", &globalX_, "globalX/F");
0523   readoutmap_->Branch("globalY", &globalY_, "globalY/F");
0524   readoutmap_->Branch("globalZ", &globalZ_, "globalZ/F");
0525 }
0526 
0527 TrackerDpgAnalysis::~TrackerDpgAnalysis() {
0528   delete[] moduleName_;
0529   delete[] moduleId_;
0530 }
0531 
0532 //
0533 // member functions
0534 //
0535 
0536 // ------------ method called to for each event  ------------
0537 void TrackerDpgAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0538   using namespace edm;
0539   using namespace reco;
0540   using namespace std;
0541   using reco::TrackCollection;
0542 
0543   // load event info
0544   eventid_ = iEvent.id().event();
0545   runid_ = iEvent.id().run();
0546   bx_ = iEvent.eventAuxiliary().bunchCrossing();
0547   orbit_ = iEvent.eventAuxiliary().orbitNumber();
0548   store_ = iEvent.eventAuxiliary().storeNumber();
0549   time_ = iEvent.eventAuxiliary().time().value();
0550   lumiSegment_ = iEvent.eventAuxiliary().luminosityBlock();
0551 
0552   // Retrieve commissioning information from "event summary", when available (for standard fine delay)
0553   edm::Handle<SiStripEventSummary> summary;
0554   iEvent.getByToken(summaryToken_, summary);
0555   if (summary.isValid())
0556     delay_ = delay(*summary.product());
0557   else
0558     delay_ = 0.;
0559 
0560   // -- Magnetic field
0561   fBz_ = fabs(iSetup.getData(magFieldToken_).inTesla(GlobalPoint(0, 0, 0)).z());
0562 
0563   siStripClusterInfo_.initEvent(iSetup);
0564 
0565   // load trigger info
0566   edm::Handle<L1GlobalTriggerReadoutRecord> gtrr_handle;
0567   iEvent.getByToken(L1Token_, gtrr_handle);
0568   L1GlobalTriggerReadoutRecord const* gtrr = gtrr_handle.product();
0569   L1GtFdlWord fdlWord = gtrr->gtFdlWord();
0570   DecisionWord L1decision = fdlWord.gtDecisionWord();
0571   for (int bit = 0; bit < 128; ++bit) {
0572     L1DecisionBits_[bit] = L1decision[bit];
0573   }
0574   DecisionWordExtended L1decisionE = fdlWord.gtDecisionWordExtended();
0575   for (int bit = 0; bit < 64; ++bit) {
0576     L1DecisionBits_[bit + 128] = L1decisionE[bit];
0577   }
0578   TechnicalTriggerWord L1technical = fdlWord.gtTechnicalTriggerWord();
0579   for (int bit = 0; bit < 64; ++bit) {
0580     L1TechnicalBits_[bit] = L1technical[bit];
0581   }
0582   orbitL1_ = fdlWord.orbitNr();
0583   physicsDeclared_ = fdlWord.physicsDeclared();
0584   edm::Handle<edm::TriggerResults> trh;
0585   iEvent.getByToken(HLTToken_, trh);
0586   size_t ntrh = trh->size();
0587   for (size_t bit = 0; bit < 256; ++bit)
0588     HLTDecisionBits_[bit] = bit < ntrh ? (bool)(trh->accept(bit)) : false;
0589 
0590   // load beamspot
0591   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0592   iEvent.getByToken(bsToken_, recoBeamSpotHandle);
0593   reco::BeamSpot bs = *recoBeamSpotHandle;
0594   const Point beamSpot = recoBeamSpotHandle.isValid()
0595                              ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0())
0596                              : Point(0, 0, 0);
0597   if (recoBeamSpotHandle.isValid()) {
0598     bsX0_ = bs.x0();
0599     bsY0_ = bs.y0();
0600     bsZ0_ = bs.z0();
0601     bsSigmaZ_ = bs.sigmaZ();
0602     bsDxdz_ = bs.dxdz();
0603     bsDydz_ = bs.dydz();
0604   } else {
0605     bsX0_ = 0.;
0606     bsY0_ = 0.;
0607     bsZ0_ = 0.;
0608     bsSigmaZ_ = 0.;
0609     bsDxdz_ = 0.;
0610     bsDydz_ = 0.;
0611   }
0612 
0613   // load primary vertex
0614   static const reco::VertexCollection s_empty_vertexColl;
0615   edm::Handle<reco::VertexCollection> vertexCollectionHandle;
0616   iEvent.getByToken(vertexToken_, vertexCollectionHandle);
0617   const reco::VertexCollection vertexColl = *(vertexCollectionHandle.product());
0618   nVertices_ = 0;
0619   for (reco::VertexCollection::const_iterator v = vertexColl.begin(); v != vertexColl.end(); ++v) {
0620     if (v->isValid() && !v->isFake())
0621       ++nVertices_;
0622   }
0623 
0624   // load pixel vertices
0625   // Pixel vertices are handled as primary vertices, but not linked to tracks.
0626   edm::Handle<reco::VertexCollection> pixelVertexCollectionHandle;
0627   iEvent.getByToken(pixelVertexToken_, pixelVertexCollectionHandle);
0628   const reco::VertexCollection pixelVertexColl = *(pixelVertexCollectionHandle.product());
0629   nPixelVertices_ = pixelVertexColl.size();
0630 
0631   // load the clusters
0632   edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
0633   iEvent.getByToken(clusterToken_, clusters);
0634   edm::Handle<edmNew::DetSetVector<SiPixelCluster> > pixelclusters;
0635   iEvent.getByToken(pixelclusterToken_, pixelclusters);
0636 
0637   // load dedx info
0638   Handle<ValueMap<DeDxData> > dEdx1Handle;
0639   Handle<ValueMap<DeDxData> > dEdx2Handle;
0640   Handle<ValueMap<DeDxData> > dEdx3Handle;
0641   try {
0642     iEvent.getByToken(dedx1Token_, dEdx1Handle);
0643   } catch (cms::Exception&) {
0644     ;
0645   }
0646   try {
0647     iEvent.getByToken(dedx2Token_, dEdx2Handle);
0648   } catch (cms::Exception&) {
0649     ;
0650   }
0651   try {
0652     iEvent.getByToken(dedx3Token_, dEdx3Handle);
0653   } catch (cms::Exception&) {
0654     ;
0655   }
0656   const ValueMap<DeDxData> dEdxTrack1 = *dEdx1Handle.product();
0657   const ValueMap<DeDxData> dEdxTrack2 = *dEdx2Handle.product();
0658   const ValueMap<DeDxData> dEdxTrack3 = *dEdx3Handle.product();
0659 
0660   // load track collections
0661   const size_t trackSize(trackLabels_.size());
0662   std::vector<reco::TrackCollection> trackCollection;
0663   std::vector<edm::Handle<reco::TrackCollection> > trackCollectionHandle;
0664   trackCollectionHandle.resize(trackSize);
0665   size_t index = 0;
0666   for (std::vector<edm::EDGetTokenT<reco::TrackCollection> >::const_iterator token = trackTokens_.begin();
0667        token != trackTokens_.end();
0668        ++token, ++index) {
0669     try {
0670       iEvent.getByToken(*token, trackCollectionHandle[index]);
0671     } catch (cms::Exception&) {
0672       ;
0673     }
0674     trackCollection.push_back(*trackCollectionHandle[index].product());
0675     ntracks_[index] = trackCollection[index].size();
0676   }
0677 
0678   // load the trajectory collections
0679   std::vector<std::vector<Trajectory> > trajectoryCollection;
0680   std::vector<edm::Handle<std::vector<Trajectory> > > trajectoryCollectionHandle;
0681   trajectoryCollectionHandle.resize(trackSize);
0682   index = 0;
0683   for (std::vector<edm::EDGetTokenT<std::vector<Trajectory> > >::const_iterator token = trajectoryTokens_.begin();
0684        token != trajectoryTokens_.end();
0685        ++token, ++index) {
0686     try {
0687       iEvent.getByToken(*token, trajectoryCollectionHandle[index]);
0688     } catch (cms::Exception&) {
0689       ;
0690     }
0691     trajectoryCollection.push_back(*trajectoryCollectionHandle[index].product());
0692     ntrajs_[index] = trajectoryCollection[index].size();
0693   }
0694 
0695   // load the tracks/traj association maps
0696   std::vector<TrajTrackAssociationCollection> TrajToTrackMap;
0697   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0698   for (std::vector<edm::EDGetTokenT<TrajTrackAssociationCollection> >::const_iterator token =
0699            trajTrackAssoTokens_.begin();
0700        token != trajTrackAssoTokens_.end();
0701        ++token) {
0702     try {
0703       iEvent.getByToken(*token, trajTrackAssociationHandle);
0704     } catch (cms::Exception&) {
0705       ;
0706     }
0707     TrajToTrackMap.push_back(*trajTrackAssociationHandle.product());
0708   }
0709 
0710   // sanity check
0711   if (!(!trackCollection.empty() && !trajectoryCollection.empty()))
0712     return;
0713 
0714   // build the reverse map tracks -> vertex
0715   std::vector<std::map<size_t, int> > trackVertices;
0716   for (size_t i = 0; i < trackSize; ++i) {
0717     trackVertices.push_back(inVertex(trackCollection[0], vertexColl, globalvertexid_ + 1));
0718   }
0719 
0720   // iterate over vertices
0721   if (functionality_vertices_) {
0722     for (reco::VertexCollection::const_iterator v = vertexColl.begin(); v != vertexColl.end(); ++v) {
0723       nTracks_pvtx_ = v->tracksSize();
0724       sumptsq_pvtx_ = sumPtSquared(*v);
0725       isValid_pvtx_ = int(v->isValid());
0726       isFake_pvtx_ = int(v->isFake());
0727       recx_pvtx_ = v->x();
0728       recy_pvtx_ = v->y();
0729       recz_pvtx_ = v->z();
0730       recx_err_pvtx_ = v->xError();
0731       recy_err_pvtx_ = v->yError();
0732       recz_err_pvtx_ = v->zError();
0733       globalvertexid_++;
0734       vertices_->Fill();
0735     }
0736   }
0737 
0738   // iterate over pixel vertices
0739   if (functionality_pixvertices_) {
0740     for (reco::VertexCollection::const_iterator v = pixelVertexColl.begin(); v != pixelVertexColl.end(); ++v) {
0741       nTracks_pvtx_ = v->tracksSize();
0742       sumptsq_pvtx_ = sumPtSquared(*v);
0743       isValid_pvtx_ = int(v->isValid());
0744       isFake_pvtx_ = int(v->isFake());
0745       recx_pvtx_ = v->x();
0746       recy_pvtx_ = v->y();
0747       recz_pvtx_ = v->z();
0748       recx_err_pvtx_ = v->xError();
0749       recy_err_pvtx_ = v->yError();
0750       recz_err_pvtx_ = v->zError();
0751       pixelVertices_->Fill();
0752     }
0753   }
0754 
0755   // determine if each cluster is on a track or not, and record the local angle
0756   // to do this, we use the first track/traj collection
0757   std::vector<double> clusterOntrackAngles = onTrackAngles(clusters, trajectoryCollection[0]);
0758   std::vector<std::pair<double, double> > pixclusterOntrackAngles =
0759       onTrackAngles(pixelclusters, trajectoryCollection[0]);
0760 
0761   /*
0762   // iterate over trajectories
0763   // note: when iterating over trajectories, it might be simpler to use the tracks/trajectories association map
0764   for(std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj< trajVec.end(); ++traj) {
0765   }
0766   // loop over all rechits from trajectories
0767   //iterate over trajectories
0768   for(std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj< trajVec.end(); ++traj) {
0769     Trajectory::DataContainer measurements = traj->measurements();
0770     // iterate over measurements
0771     for(Trajectory::DataContainer::iterator meas = measurements.begin(); meas!= measurements.end(); ++meas) {
0772     }
0773   }
0774 */
0775 
0776   // determine if each cluster is on a track or not, and record the trackid
0777   std::vector<std::vector<int> > stripClusterOntrackIndices;
0778   for (size_t i = 0; i < trackSize; ++i) {
0779     stripClusterOntrackIndices.push_back(onTrack(clusters, trackCollection[i], globaltrackid_[i] + 1));
0780   }
0781   std::vector<std::vector<int> > pixelClusterOntrackIndices;
0782   for (size_t i = 0; i < trackSize; ++i) {
0783     pixelClusterOntrackIndices.push_back(onTrack(pixelclusters, trackCollection[i], globaltrackid_[i] + 1));
0784   }
0785   nclustersOntrack_ = count_if(
0786       stripClusterOntrackIndices[0].begin(), stripClusterOntrackIndices[0].end(), [](auto c) { return c != -1; });
0787   npixClustersOntrack_ = count_if(
0788       pixelClusterOntrackIndices[0].begin(), pixelClusterOntrackIndices[0].end(), [](auto c) { return c != -1; });
0789 
0790   // iterate over tracks
0791   for (size_t coll = 0; coll < trackCollection.size(); ++coll) {
0792     uint32_t n_hits_barrel = 0;
0793     uint32_t n_hits_lowprob = 0;
0794     for (TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap[coll].begin();
0795          it != TrajToTrackMap[coll].end();
0796          ++it) {
0797       reco::TrackRef itTrack = it->val;
0798       edm::Ref<std::vector<Trajectory> > traj = it->key;  // bug to find type of the key
0799       eta_ = itTrack->eta();
0800       phi_ = itTrack->phi();
0801       try {  // not all track collections have the dedx info... indeed at best one.
0802         dedxNoM_ = dEdxTrack1[itTrack].numberOfMeasurements();
0803         dedx1_ = dEdxTrack1[itTrack].dEdx();
0804         dedx2_ = dEdxTrack2[itTrack].dEdx();
0805         dedx3_ = dEdxTrack3[itTrack].dEdx();
0806       } catch (cms::Exception&) {
0807         dedxNoM_ = 0;
0808         dedx1_ = 0.;
0809         dedx2_ = 0.;
0810         dedx3_ = 0.;
0811       }
0812       charge_ = itTrack->charge();
0813       quality_ = itTrack->qualityMask();
0814       foundhits_ = itTrack->found();
0815       lostHits_ = itTrack->lost();
0816       foundhitsStrips_ = itTrack->hitPattern().numberOfValidStripHits();
0817       foundhitsPixels_ = itTrack->hitPattern().numberOfValidPixelHits();
0818       losthitsStrips_ = itTrack->hitPattern().numberOfLostStripHits(reco::HitPattern::TRACK_HITS);
0819       losthitsPixels_ = itTrack->hitPattern().numberOfLostPixelHits(reco::HitPattern::TRACK_HITS);
0820       nLayers_ = uint32_t(itTrack->hitPattern().trackerLayersWithMeasurement());
0821       p_ = itTrack->p();
0822       pt_ = itTrack->pt();
0823       chi2_ = itTrack->chi2();
0824       ndof_ = (uint32_t)itTrack->ndof();
0825       dz_ = itTrack->dz();
0826       dzerr_ = itTrack->dzError();
0827       dzCorr_ = itTrack->dz(beamSpot);
0828       dxy_ = itTrack->dxy();
0829       dxyerr_ = itTrack->dxyError();
0830       dxyCorr_ = itTrack->dxy(beamSpot);
0831       pterr_ = itTrack->ptError();
0832       etaerr_ = itTrack->etaError();
0833       phierr_ = itTrack->phiError();
0834       qoverp_ = itTrack->qoverp();
0835       xPCA_ = itTrack->vertex().x();
0836       yPCA_ = itTrack->vertex().y();
0837       zPCA_ = itTrack->vertex().z();
0838       try {  // only one track collection (at best) is connected to the main vertex
0839         if (!vertexColl.empty() && !vertexColl.begin()->isFake()) {
0840           trkWeightpvtx_ = vertexColl.begin()->trackWeight(itTrack);
0841         } else
0842           trkWeightpvtx_ = 0.;
0843       } catch (cms::Exception&) {
0844         trkWeightpvtx_ = 0.;
0845       }
0846       globaltrackid_[coll]++;
0847       std::map<size_t, int>::const_iterator theV = trackVertices[coll].find(itTrack.key());
0848       vertexid_ = (theV != trackVertices[coll].end()) ? theV->second : 0;
0849       // add missing hits (separate tree, common strip + pixel)
0850       Trajectory::DataContainer const& measurements = traj->measurements();
0851       if (functionality_missingHits_) {
0852         for (Trajectory::DataContainer::const_iterator it = measurements.begin(); it != measurements.end(); ++it) {
0853           TrajectoryMeasurement::ConstRecHitPointer rechit = it->recHit();
0854           if (!rechit->isValid()) {
0855             // detid
0856             detid_ = rechit->geographicalId();
0857             // status
0858             type_ = rechit->getType();
0859             // position
0860             LocalPoint local = it->predictedState().localPosition();
0861             clPositionX_ = local.x();
0862             clPositionY_ = local.y();
0863             // global position
0864             GlobalPoint global = it->predictedState().globalPosition();
0865             globalX_ = global.x();
0866             globalY_ = global.y();
0867             globalZ_ = global.z();
0868             // position in the measurement frame
0869             measX_ = 0;
0870             measY_ = 0;
0871             if (type_ != TrackingRecHit::inactive) {
0872               const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker_->idToDetUnit(detid_));
0873               if (gdu && gdu->type().isTracker()) {
0874                 const Topology& topo = gdu->topology();
0875                 MeasurementPoint meas = topo.measurementPosition(local);
0876                 measX_ = meas.x();
0877                 measY_ = meas.y();
0878               }
0879             }
0880             // local error
0881             LocalError error = it->predictedState().localError().positionError();
0882             errorX_ = error.xx();
0883             errorY_ = error.yy();
0884             // fill
0885             missingHits_[coll]->Fill();
0886           }
0887         }
0888       }
0889       // compute the fraction of low probability pixels... will be added to the event tree
0890       for (trackingRecHit_iterator it = itTrack->recHitsBegin(); it != itTrack->recHitsEnd(); ++it) {
0891         const TrackingRecHit* hit = &(**it);
0892         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
0893         if (pixhit) {
0894           DetId detId = pixhit->geographicalId();
0895           if (detId.subdetId() == PixelSubdetector::PixelBarrel) {
0896             ++n_hits_barrel;
0897             double proba = pixhit->clusterProbability(0);
0898             if (proba <= 0.0)
0899               ++n_hits_lowprob;
0900           }
0901         }
0902       }
0903       // fill the track tree
0904       if (functionality_tracks_)
0905         tracks_[coll]->Fill();
0906     }
0907     lowPixelProbabilityFraction_[coll] = n_hits_barrel > 0 ? (float)n_hits_lowprob / n_hits_barrel : -1.;
0908   }
0909 
0910   // iterate over clusters
0911   nclusters_ = 0;
0912   std::vector<double>::const_iterator angleIt = clusterOntrackAngles.begin();
0913   uint32_t localCounter = 0;
0914   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
0915        DSViter++) {
0916     edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
0917     edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
0918     uint32_t detid = DSViter->id();
0919     nclusters_ += DSViter->size();
0920     if (functionality_offtrackClusters_ || functionality_ontrackClusters_) {
0921       for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end;
0922            ++iter, ++angleIt, ++localCounter) {
0923         siStripClusterInfo_.setCluster(*iter, detid);
0924         // general quantities
0925         for (size_t i = 0; i < trackSize; ++i) {
0926           trackid_[i] = stripClusterOntrackIndices[i][localCounter];
0927         }
0928         onTrack_ = (trackid_[0] != (uint32_t)-1);
0929         clWidth_ = siStripClusterInfo_.width();
0930         clPosition_ = siStripClusterInfo_.baryStrip();
0931         angle_ = *angleIt;
0932         thickness_ = ((((DSViter->id() >> 25) & 0x7f) == 0xd) ||
0933                       ((((DSViter->id() >> 25) & 0x7f) == 0xe) && (((DSViter->id() >> 5) & 0x7) > 4)))
0934                          ? 500
0935                          : 300;
0936         stripLength_ = static_cast<const StripGeomDetUnit*>(tracker_->idToDet(detid))->specificTopology().stripLength();
0937         int nstrips = static_cast<const StripGeomDetUnit*>(tracker_->idToDet(detid))->specificTopology().nstrips();
0938         maxCharge_ = siStripClusterInfo_.maxCharge();
0939         // signal and noise with gain corrections
0940         clNormalizedCharge_ = siStripClusterInfo_.charge();
0941         clNormalizedNoise_ = siStripClusterInfo_.noiseRescaledByGain();
0942         clSignalOverNoise_ = siStripClusterInfo_.signalOverNoise();
0943         // signal and noise with gain corrections and angle corrections
0944         clCorrectedCharge_ = clNormalizedCharge_ * fabs(cos(angle_));          // corrected for track angle
0945         clCorrectedSignalOverNoise_ = clSignalOverNoise_ * fabs(cos(angle_));  // corrected for track angle
0946         // signal and noise without gain corrections
0947         clBareNoise_ = siStripClusterInfo_.noise();
0948         clBareCharge_ = clSignalOverNoise_ * clBareNoise_;
0949         // global position
0950         const StripGeomDetUnit* sgdu = static_cast<const StripGeomDetUnit*>(tracker_->idToDet(detid));
0951         Surface::GlobalPoint gp =
0952             sgdu->surface().toGlobal(sgdu->specificTopology().localPosition(MeasurementPoint(clPosition_, 0)));
0953         globalX_ = gp.x();
0954         globalY_ = gp.y();
0955         globalZ_ = gp.z();
0956         // cabling
0957         detid_ = detid;
0958         lldChannel_ = 1 + (int(floor(iter->barycenter())) / 256);
0959         if (lldChannel_ == 2 && nstrips == 512)
0960           lldChannel_ = 3;
0961         if ((functionality_offtrackClusters_ && !onTrack_) || (functionality_ontrackClusters_ && onTrack_))
0962           clusters_->Fill();
0963       }
0964     }
0965   }
0966 
0967   // iterate over pixel clusters
0968   npixClusters_ = 0;
0969   std::vector<std::pair<double, double> >::const_iterator pixAngleIt = pixclusterOntrackAngles.begin();
0970   localCounter = 0;
0971   for (edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = pixelclusters->begin();
0972        DSViter != pixelclusters->end();
0973        DSViter++) {
0974     edmNew::DetSet<SiPixelCluster>::const_iterator begin = DSViter->begin();
0975     edmNew::DetSet<SiPixelCluster>::const_iterator end = DSViter->end();
0976     uint32_t detid = DSViter->id();
0977     npixClusters_ += DSViter->size();
0978     if (functionality_pixclusters_) {
0979       for (edmNew::DetSet<SiPixelCluster>::const_iterator iter = begin; iter != end;
0980            ++iter, ++pixAngleIt, ++localCounter) {
0981         // general quantities
0982         for (size_t i = 0; i < trackSize; ++i) {
0983           trackid_[i] = pixelClusterOntrackIndices[i][localCounter];
0984         }
0985         onTrack_ = (trackid_[0] != (uint32_t)-1);
0986         clPositionX_ = iter->x();
0987         clPositionY_ = iter->y();
0988         clSize_ = iter->size();
0989         clSizeX_ = iter->sizeX();
0990         clSizeY_ = iter->sizeY();
0991         alpha_ = pixAngleIt->first;
0992         beta_ = pixAngleIt->second;
0993         charge_ = (iter->charge()) / 1000.;
0994         chargeCorr_ = charge_ * sqrt(1.0 / (1.0 / pow(tan(alpha_), 2) + 1.0 / pow(tan(beta_), 2) + 1.0)) / 1000.;
0995         // global position
0996         const PixelGeomDetUnit* pgdu = static_cast<const PixelGeomDetUnit*>(tracker_->idToDet(detid));
0997         Surface::GlobalPoint gp = pgdu->surface().toGlobal(
0998             pgdu->specificTopology().localPosition(MeasurementPoint(clPositionX_, clPositionY_)));
0999         globalX_ = gp.x();
1000         globalY_ = gp.y();
1001         globalZ_ = gp.z();
1002         // cabling
1003         detid_ = detid;
1004         // fill
1005         pixclusters_->Fill();
1006       }
1007     }
1008   }
1009 
1010   // topological quantities - uses the first track collection
1011   EventShape shape(trackCollection[0]);
1012   math::XYZTLorentzVectorF thrust = shape.thrust();
1013   thrustValue_ = thrust.t();
1014   thrustX_ = thrust.x();
1015   thrustY_ = thrust.y();
1016   thrustZ_ = thrust.z();
1017   sphericity_ = shape.sphericity();
1018   planarity_ = shape.planarity();
1019   aplanarity_ = shape.aplanarity();
1020 
1021   // fill event tree
1022   if (functionality_events_)
1023     event_->Fill();
1024 }
1025 
1026 // ------------ method called once each job just before starting event loop  ------------
1027 void TrackerDpgAnalysis::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
1028   const auto& tTopo = iSetup.getData(tTopoToken_);
1029   tracker_ = &iSetup.getData(tkGeomToken_);
1030 
1031   //HLT names
1032   bool changed(true);
1033   if (hltConfig_.init(iRun, iSetup, HLTTag_.process(), changed)) {
1034     if (changed) {
1035       hlNames_ = hltConfig_.triggerNames();
1036     }
1037   }
1038   int i = 0;
1039   for (std::vector<std::string>::const_iterator it = hlNames_.begin(); it < hlNames_.end(); ++it) {
1040     std::cout << (i++) << " = " << (*it) << std::endl;
1041   }
1042 
1043   // read the delay offsets for each device from input files
1044   // this is only for the so-called "random delay" run
1045   std::map<uint32_t, float> delayMap = delay(delayFileNames_);
1046   TrackerMap tmap("Delays");
1047 
1048   // cabling I (readout)
1049   const auto& cabling = iSetup.getData(fedCablingToken_);
1050   auto feds = cabling.fedIds();
1051   for (auto fedid = feds.begin(); fedid < feds.end(); ++fedid) {
1052     auto connections = cabling.fedConnections(*fedid);
1053     for (auto conn = connections.begin(); conn < connections.end(); ++conn) {
1054       // Fill the "old" map to be used for lookup during analysis
1055       if (conn->isConnected())
1056         connections_.insert(std::make_pair(conn->detId(), new FedChannelConnection(*conn)));
1057       // Fill the standalone tree (once for all)
1058       if (conn->isConnected()) {
1059         detid_ = conn->detId();
1060         strncpy(moduleName_, toStringName(detid_, tTopo).c_str(), 256);
1061         strncpy(moduleId_, toStringId(detid_).c_str(), 256);
1062         lldChannel_ = conn->lldChannel();
1063         dcuId_ = conn->dcuId();
1064         fecCrate_ = conn->fecCrate();
1065         fecSlot_ = conn->fecSlot();
1066         fecRing_ = conn->fecRing();
1067         ccuAdd_ = conn->ccuAddr();
1068         ccuChan_ = conn->ccuChan();
1069         fedId_ = conn->fedId();
1070         fedCh_ = conn->fedCh();
1071         fiberLength_ = conn->fiberLength();
1072         delay_ = delayMap[dcuId_];
1073         const StripGeomDetUnit* sgdu = static_cast<const StripGeomDetUnit*>(tracker_->idToDet(detid_));
1074         Surface::GlobalPoint gp = sgdu->surface().toGlobal(LocalPoint(0, 0));
1075         globalX_ = gp.x();
1076         globalY_ = gp.y();
1077         globalZ_ = gp.z();
1078         readoutmap_->Fill();
1079         tmap.fill_current_val(detid_, delay_);
1080       }
1081     }
1082   }
1083   if (!delayMap.empty())
1084     tmap.save(true, 0, 0, "delaymap.png");
1085 
1086   // cabling II (DCU map)
1087   std::ifstream cablingFile(cablingFileName_.c_str());
1088   if (cablingFile.is_open()) {
1089     char buffer[1024];
1090     cablingFile.getline(buffer, 1024);
1091     while (!cablingFile.eof()) {
1092       std::istringstream line(buffer);
1093       std::string name;
1094       // one line contains the PSU name + all dcuids connected to it.
1095       line >> name;
1096       strncpy(PSUname_, name.c_str(), 256);
1097       while (!line.eof()) {
1098         line >> dcuId_;
1099         psumap_->Fill();
1100       }
1101       cablingFile.getline(buffer, 1024);
1102     }
1103   } else {
1104     edm::LogWarning("BadConfig") << " The PSU file does not exist. The psumap tree will not be filled." << std::endl
1105                                  << " Looking for " << cablingFileName_.c_str() << "." << std::endl
1106                                  << " Please specify a valid filename through the PSUFileName untracked parameter.";
1107   }
1108 }
1109 
1110 // ------------ method called once each job just after ending the event loop  ------------
1111 void TrackerDpgAnalysis::endJob() {
1112   for (size_t i = 0; i < tracks_.size(); ++i) {
1113     char buffer[256];
1114     sprintf(buffer, "trackid%lu", (unsigned long)i);
1115     if (tracks_[i]->GetEntries())
1116       tracks_[i]->BuildIndex(buffer, "eventid");
1117   }
1118   /* not needed: missing hits is a high-level quantity
1119   for(size_t i = 0; i<missingHits_.size();++i) {
1120     char buffer[256];
1121     sprintf(buffer,"trackid%lu",(unsigned long)i);
1122     if(missingHits_[i]->GetEntries()) missingHits_[i]->BuildIndex(buffer);
1123   }
1124   */
1125   if (vertices_->GetEntries())
1126     vertices_->BuildIndex("vertexid", "eventid");
1127   if (event_->GetEntries())
1128     event_->BuildIndex("runid", "eventid");
1129   if (psumap_->GetEntries())
1130     psumap_->BuildIndex("dcuId");
1131   if (readoutmap_->GetEntries())
1132     readoutmap_->BuildIndex("detid", "lldChannel");
1133 }
1134 
1135 std::vector<double> TrackerDpgAnalysis::onTrackAngles(edm::Handle<edmNew::DetSetVector<SiStripCluster> >& clusters,
1136                                                       const std::vector<Trajectory>& trajVec) {
1137   std::vector<double> result;
1138   // first, build a list of positions and angles on trajectories
1139   std::multimap<const uint32_t, std::pair<LocalPoint, double> > onTrackPositions;
1140   for (std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj < trajVec.end(); ++traj) {
1141     Trajectory::DataContainer measurements = traj->measurements();
1142     for (Trajectory::DataContainer::iterator meas = measurements.begin(); meas != measurements.end(); ++meas) {
1143       double tla = meas->updatedState().localDirection().theta();
1144       insertMeasurement(onTrackPositions, &(*(meas->recHit())), tla);
1145     }
1146   }
1147   // then loop over the clusters to check
1148   double angle = 0.;
1149   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
1150        DSViter++) {
1151     edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
1152     edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
1153     std::pair<std::multimap<uint32_t, std::pair<LocalPoint, double> >::const_iterator,
1154               std::multimap<uint32_t, std::pair<LocalPoint, double> >::const_iterator>
1155         range = onTrackPositions.equal_range(DSViter->id());
1156     const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker_->idToDet(DSViter->id()));
1157     for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
1158       angle = 0.;
1159       for (std::multimap<uint32_t, std::pair<LocalPoint, double> >::const_iterator cl = range.first; cl != range.second;
1160            ++cl) {
1161         if (fabs(gdu->topology().measurementPosition(cl->second.first).x() - iter->barycenter()) < 2) {
1162           angle = cl->second.second;
1163         }
1164       }
1165       result.push_back(angle);
1166     }
1167   }
1168   return result;
1169 }
1170 
1171 void TrackerDpgAnalysis::insertMeasurement(std::multimap<const uint32_t, std::pair<LocalPoint, double> >& collection,
1172                                            const TransientTrackingRecHit* hit,
1173                                            double tla) {
1174   if (!hit)
1175     return;
1176   const SiTrackerMultiRecHit* multihit = dynamic_cast<const SiTrackerMultiRecHit*>(hit);
1177   const SiStripRecHit2D* singlehit = dynamic_cast<const SiStripRecHit2D*>(hit);
1178   const SiStripRecHit1D* hit1d = dynamic_cast<const SiStripRecHit1D*>(hit);
1179   if (hit1d) {  //...->33X
1180     collection.insert(std::make_pair(hit1d->geographicalId().rawId(), std::make_pair(hit1d->localPosition(), tla)));
1181   } else if (singlehit) {  // 41X->...
1182     collection.insert(
1183         std::make_pair(singlehit->geographicalId().rawId(), std::make_pair(singlehit->localPosition(), tla)));
1184   } else if (multihit) {
1185     std::vector<const TrackingRecHit*> childs = multihit->recHits();
1186     for (std::vector<const TrackingRecHit*>::const_iterator it = childs.begin(); it != childs.end(); ++it) {
1187       insertMeasurement(collection, dynamic_cast<const TrackingRecHit*>(*it), tla);
1188     }
1189   }
1190 }
1191 
1192 std::vector<int> TrackerDpgAnalysis::onTrack(edm::Handle<edmNew::DetSetVector<SiStripCluster> >& clusters,
1193                                              const reco::TrackCollection& trackVec,
1194                                              uint32_t firstTrack) {
1195   std::vector<int> result;
1196   // first, build a list of positions and trackid on tracks
1197   std::multimap<const uint32_t, std::pair<int, int> > onTrackPositions;
1198   uint32_t trackid = firstTrack;
1199   for (reco::TrackCollection::const_iterator itTrack = trackVec.begin(); itTrack != trackVec.end();
1200        ++itTrack, ++trackid) {
1201     for (trackingRecHit_iterator it = itTrack->recHitsBegin(); it != itTrack->recHitsEnd(); ++it) {
1202       const TrackingRecHit* hit = &(**it);
1203       insertMeasurement(onTrackPositions, hit, trackid);
1204     }
1205   }
1206   // then loop over the clusters to check
1207   int thetrackid = -1;
1208   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
1209        DSViter++) {
1210     edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
1211     edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
1212     std::pair<std::multimap<uint32_t, std::pair<int, int> >::const_iterator,
1213               std::multimap<uint32_t, std::pair<int, int> >::const_iterator>
1214         range = onTrackPositions.equal_range(DSViter->id());
1215     for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
1216       thetrackid = -1;
1217       for (std::multimap<uint32_t, std::pair<int, int> >::const_iterator cl = range.first; cl != range.second; ++cl) {
1218         if (fabs(cl->second.first - iter->barycenter()) < 2) {
1219           thetrackid = cl->second.second;
1220         }
1221       }
1222       result.push_back(thetrackid);
1223     }
1224   }
1225   return result;
1226 }
1227 
1228 void TrackerDpgAnalysis::insertMeasurement(std::multimap<const uint32_t, std::pair<int, int> >& collection,
1229                                            const TrackingRecHit* hit,
1230                                            int trackid) {
1231   if (!hit)
1232     return;
1233   const SiTrackerMultiRecHit* multihit = dynamic_cast<const SiTrackerMultiRecHit*>(hit);
1234   const SiStripRecHit2D* singlehit = dynamic_cast<const SiStripRecHit2D*>(hit);
1235   const SiStripRecHit1D* hit1d = dynamic_cast<const SiStripRecHit1D*>(hit);
1236   if (hit1d) {  // 41X->...
1237     collection.insert(
1238         std::make_pair(hit1d->geographicalId().rawId(), std::make_pair(int(hit1d->cluster()->barycenter()), trackid)));
1239   } else if (singlehit) {  //...->33X
1240     collection.insert(std::make_pair(singlehit->geographicalId().rawId(),
1241                                      std::make_pair(int(singlehit->cluster()->barycenter()), trackid)));
1242   } else if (multihit) {
1243     std::vector<const TrackingRecHit*> childs = multihit->recHits();
1244     for (std::vector<const TrackingRecHit*>::const_iterator it = childs.begin(); it != childs.end(); ++it) {
1245       insertMeasurement(collection, *it, trackid);
1246     }
1247   }
1248 }
1249 
1250 std::map<size_t, int> TrackerDpgAnalysis::inVertex(const reco::TrackCollection& tracks,
1251                                                    const reco::VertexCollection& vertices,
1252                                                    uint32_t firstVertex) {
1253   // build reverse map track -> vertex
1254   std::map<size_t, int> output;
1255   uint32_t vertexid = firstVertex;
1256   for (reco::VertexCollection::const_iterator v = vertices.begin(); v != vertices.end(); ++v, ++vertexid) {
1257     reco::Vertex::trackRef_iterator it = v->tracks_begin();
1258     reco::Vertex::trackRef_iterator lastTrack = v->tracks_end();
1259     for (; it != lastTrack; ++it) {
1260       output[it->key()] = vertexid;
1261     }
1262   }
1263   return output;
1264 }
1265 
1266 std::vector<std::pair<double, double> > TrackerDpgAnalysis::onTrackAngles(
1267     edm::Handle<edmNew::DetSetVector<SiPixelCluster> >& clusters, const std::vector<Trajectory>& trajVec) {
1268   std::vector<std::pair<double, double> > result;
1269   // first, build a list of positions and angles on trajectories
1270   std::multimap<const uint32_t, std::pair<LocalPoint, std::pair<double, double> > > onTrackPositions;
1271   for (std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj < trajVec.end(); ++traj) {
1272     Trajectory::DataContainer measurements = traj->measurements();
1273     for (Trajectory::DataContainer::iterator meas = measurements.begin(); meas != measurements.end(); ++meas) {
1274       LocalVector localDir = meas->updatedState().localDirection();
1275       double alpha = atan2(localDir.z(), localDir.x());
1276       double beta = atan2(localDir.z(), localDir.y());
1277       insertMeasurement(onTrackPositions, &(*(meas->recHit())), alpha, beta);
1278     }
1279   }
1280   // then loop over the clusters to check
1281   double alpha = 0.;
1282   double beta = 0.;
1283   for (edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
1284        DSViter++) {
1285     edmNew::DetSet<SiPixelCluster>::const_iterator begin = DSViter->begin();
1286     edmNew::DetSet<SiPixelCluster>::const_iterator end = DSViter->end();
1287     for (edmNew::DetSet<SiPixelCluster>::const_iterator iter = begin; iter != end; ++iter) {
1288       alpha = 0.;
1289       beta = 0.;
1290       std::pair<std::multimap<uint32_t, std::pair<LocalPoint, std::pair<double, double> > >::const_iterator,
1291                 std::multimap<uint32_t, std::pair<LocalPoint, std::pair<double, double> > >::const_iterator>
1292           range = onTrackPositions.equal_range(DSViter->id());
1293       const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker_->idToDet(DSViter->id()));
1294       for (std::multimap<uint32_t, std::pair<LocalPoint, std::pair<double, double> > >::const_iterator cl = range.first;
1295            cl != range.second;
1296            ++cl) {
1297         if (fabs(gdu->topology().measurementPosition(cl->second.first).x() - iter->x()) < 2 &&
1298             fabs(gdu->topology().measurementPosition(cl->second.first).y() - iter->y()) < 2) {
1299           alpha = cl->second.second.first;
1300           beta = cl->second.second.second;
1301         }
1302       }
1303       result.push_back(std::make_pair(alpha, beta));
1304     }
1305   }
1306   return result;
1307 }
1308 
1309 void TrackerDpgAnalysis::insertMeasurement(
1310     std::multimap<const uint32_t, std::pair<LocalPoint, std::pair<double, double> > >& collection,
1311     const TransientTrackingRecHit* hit,
1312     double alpha,
1313     double beta) {
1314   if (!hit)
1315     return;
1316   const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
1317   if (pixhit) {
1318     collection.insert(std::make_pair(pixhit->geographicalId().rawId(),
1319                                      std::make_pair(pixhit->localPosition(), std::make_pair(alpha, beta))));
1320   }
1321 }
1322 
1323 std::vector<int> TrackerDpgAnalysis::onTrack(edm::Handle<edmNew::DetSetVector<SiPixelCluster> >& clusters,
1324                                              const reco::TrackCollection& trackVec,
1325                                              uint32_t firstTrack) {
1326   std::vector<int> result;
1327   // first, build a list of positions and trackid on tracks
1328   std::multimap<const uint32_t, std::pair<std::pair<float, float>, int> > onTrackPositions;
1329   uint32_t trackid = firstTrack;
1330   for (reco::TrackCollection::const_iterator itTrack = trackVec.begin(); itTrack != trackVec.end();
1331        ++itTrack, ++trackid) {
1332     for (trackingRecHit_iterator it = itTrack->recHitsBegin(); it != itTrack->recHitsEnd(); ++it) {
1333       const TrackingRecHit* hit = &(**it);
1334       insertMeasurement(onTrackPositions, hit, trackid);
1335     }
1336   }
1337   // then loop over the clusters to check
1338   int thetrackid = -1;
1339   for (edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
1340        DSViter++) {
1341     edmNew::DetSet<SiPixelCluster>::const_iterator begin = DSViter->begin();
1342     edmNew::DetSet<SiPixelCluster>::const_iterator end = DSViter->end();
1343     for (edmNew::DetSet<SiPixelCluster>::const_iterator iter = begin; iter != end; ++iter) {
1344       thetrackid = -1;
1345       std::pair<std::multimap<uint32_t, std::pair<std::pair<float, float>, int> >::const_iterator,
1346                 std::multimap<uint32_t, std::pair<std::pair<float, float>, int> >::const_iterator>
1347           range = onTrackPositions.equal_range(DSViter->id());
1348       for (std::multimap<uint32_t, std::pair<std::pair<float, float>, int> >::const_iterator cl = range.first;
1349            cl != range.second;
1350            ++cl) {
1351         if ((fabs(cl->second.first.first - iter->x()) < 2) && (fabs(cl->second.first.second - iter->y()) < 2)) {
1352           thetrackid = cl->second.second;
1353         }
1354       }
1355       result.push_back(thetrackid);
1356     }
1357   }
1358   return result;
1359 }
1360 
1361 void TrackerDpgAnalysis::insertMeasurement(
1362     std::multimap<const uint32_t, std::pair<std::pair<float, float>, int> >& collection,
1363     const TrackingRecHit* hit,
1364     int trackid) {
1365   if (!hit)
1366     return;
1367   const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
1368   if (pixhit) {
1369     collection.insert(
1370         std::make_pair(pixhit->geographicalId().rawId(),
1371                        std::make_pair(std::make_pair(pixhit->cluster()->x(), pixhit->cluster()->y()), trackid)));
1372   }
1373 }
1374 
1375 std::string TrackerDpgAnalysis::toStringName(uint32_t rawid, const TrackerTopology& tTopo) {
1376   SiStripDetId detid(rawid);
1377   std::string out;
1378   std::stringstream output;
1379   switch (detid.subDetector()) {
1380     case 3: {
1381       output << "TIB";
1382 
1383       output << (tTopo.tibIsZPlusSide(rawid) ? "+" : "-");
1384       output << " layer ";
1385       output << tTopo.tibLayer(rawid);
1386       output << ", string ";
1387       output << tTopo.tibString(rawid);
1388       output << (tTopo.tibIsExternalString(rawid) ? " external" : " internal");
1389       output << ", module ";
1390       output << tTopo.tibModule(rawid);
1391       if (tTopo.tibIsDoubleSide(rawid)) {
1392         output << " (double)";
1393       } else {
1394         output << (tTopo.tibIsRPhi(rawid) ? " (rphi)" : " (stereo)");
1395       }
1396       break;
1397     }
1398     case 4: {
1399       output << "TID";
1400 
1401       output << (tTopo.tidIsZPlusSide(rawid) ? "+" : "-");
1402       output << " disk ";
1403       output << tTopo.tidWheel(rawid);
1404       output << ", ring ";
1405       output << tTopo.tidRing(rawid);
1406       output << (tTopo.tidIsFrontRing(rawid) ? " front" : " back");
1407       output << ", module ";
1408       output << tTopo.tidModule(rawid);
1409       if (tTopo.tidIsDoubleSide(rawid)) {
1410         output << " (double)";
1411       } else {
1412         output << (tTopo.tidIsRPhi(rawid) ? " (rphi)" : " (stereo)");
1413       }
1414       break;
1415     }
1416     case 5: {
1417       output << "TOB";
1418 
1419       output << (tTopo.tobIsZPlusSide(rawid) ? "+" : "-");
1420       output << " layer ";
1421       output << tTopo.tobLayer(rawid);
1422       output << ", rod ";
1423       output << tTopo.tobRod(rawid);
1424       output << ", module ";
1425       output << tTopo.tobModule(rawid);
1426       if (tTopo.tobIsDoubleSide(rawid)) {
1427         output << " (double)";
1428       } else {
1429         output << (tTopo.tobIsRPhi(rawid) ? " (rphi)" : " (stereo)");
1430       }
1431       break;
1432     }
1433     case 6: {
1434       output << "TEC";
1435 
1436       output << (tTopo.tecIsZPlusSide(rawid) ? "+" : "-");
1437       output << " disk ";
1438       output << tTopo.tecWheel(rawid);
1439       output << " sector ";
1440       output << tTopo.tecPetalNumber(rawid);
1441       output << (tTopo.tecIsFrontPetal(rawid) ? " Front Petal" : " Back Petal");
1442       output << ", module ";
1443       output << tTopo.tecRing(rawid);
1444       output << tTopo.tecModule(rawid);
1445       if (tTopo.tecIsDoubleSide(rawid)) {
1446         output << " (double)";
1447       } else {
1448         output << (tTopo.tecIsRPhi(rawid) ? " (rphi)" : " (stereo)");
1449       }
1450       break;
1451     }
1452     default: {
1453       output << "UNKNOWN";
1454     }
1455   }
1456   out = output.str();
1457   return out;
1458 }
1459 
1460 std::string TrackerDpgAnalysis::toStringId(uint32_t rawid) {
1461   std::string out;
1462   std::stringstream output;
1463   output << rawid << " (0x" << std::hex << rawid << std::dec << ")";
1464   out = output.str();
1465   return out;
1466 }
1467 
1468 double TrackerDpgAnalysis::sumPtSquared(const reco::Vertex& v) {
1469   double sum = 0.;
1470   double pT;
1471   for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1472     pT = (**it).pt();
1473     sum += pT * pT;
1474   }
1475   return sum;
1476 }
1477 
1478 float TrackerDpgAnalysis::delay(const SiStripEventSummary& summary) {
1479   float delay = const_cast<SiStripEventSummary&>(summary).ttcrx();
1480   uint32_t latencyCode = (const_cast<SiStripEventSummary&>(summary).layerScanned() >> 24) & 0xff;
1481   int latencyShift =
1482       latencyCode & 0x3f;  // number of bunch crossings between current value and start of scan... must be positive
1483   if (latencyShift > 32)
1484     latencyShift -= 64;  // allow negative values: we cover [-32,32].. should not be needed.
1485   if ((latencyCode >> 6) == 2)
1486     latencyShift -= 3;  // layer in deconv, rest in peak
1487   if ((latencyCode >> 6) == 1)
1488     latencyShift += 3;  // layer in peak, rest in deconv
1489   float correctedDelay =
1490       delay - (latencyShift * 25.);  // shifts the delay so that 0 corresponds to the current settings.
1491   return correctedDelay;
1492 }
1493 
1494 std::map<uint32_t, float> TrackerDpgAnalysis::delay(const std::vector<std::string>& files) {
1495   // prepare output
1496   uint32_t dcuid;
1497   float delay;
1498   std::map<uint32_t, float> delayMap;
1499   //iterator over input files
1500   for (std::vector<std::string>::const_iterator file = files.begin(); file < files.end(); ++file) {
1501     // open the file
1502     std::ifstream cablingFile(file->c_str());
1503     if (cablingFile.is_open()) {
1504       char buffer[1024];
1505       // read one line
1506       cablingFile.getline(buffer, 1024);
1507       while (!cablingFile.eof()) {
1508         std::string line(buffer);
1509         size_t pos = line.find("dcuid");
1510         // one line containing dcuid
1511         if (pos != std::string::npos) {
1512           // decode dcuid
1513           std::string dcuids = line.substr(pos + 7, line.find(' ', pos) - pos - 8);
1514           std::istringstream dcuidstr(dcuids);
1515           dcuidstr >> std::hex >> dcuid;
1516           // decode delay
1517           pos = line.find("difpll");
1518           std::string diffs = line.substr(pos + 8, line.find(' ', pos) - pos - 9);
1519           std::istringstream diffstr(diffs);
1520           diffstr >> delay;
1521           // fill the map
1522           delayMap[dcuid] = delay;
1523         }
1524         // iterate
1525         cablingFile.getline(buffer, 1024);
1526       }
1527     } else {
1528       edm::LogWarning("BadConfig") << " The delay file does not exist. The delay map will not be filled properly."
1529                                    << std::endl
1530                                    << " Looking for " << file->c_str() << "." << std::endl
1531                                    << " Please specify valid filenames through the DelayFileNames untracked parameter.";
1532     }
1533   }
1534   return delayMap;
1535 }
1536 
1537 //define this as a plug-in
1538 DEFINE_FWK_MODULE(TrackerDpgAnalysis);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!