Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:56 

0001 // -*- C++ -*-
0002 //
0003 // Package:    CalibTracker/SiStripHitResolution
0004 // Class:      SiStripCPEAnalyzer
0005 //
0006 /**\class SiStripCPEAnalyzer SiStripCPEAnalyzer.cc CalibTracker/SiStripHitResolution/plugins/SiStripCPEAnalyzer.cc
0007 
0008  Description: CPE analysis (resolution param, etc.)
0009 
0010  Implementation:
0011      Used for Strip Tracker DPG studies
0012 */
0013 //
0014 // Original Author:  Christophe Delaere
0015 //         Created:  Thu, 12 Sep 2019 13:51:00 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022 #include <algorithm>
0023 
0024 // user include files
0025 #include "CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027 #include "DataFormats/Common/interface/DetSetVector.h"
0028 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0031 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0040 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0043 #include "RecoLocalTracker/Records/interface/TkStripCPERecord.h"
0044 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0045 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0046 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0047 
0048 #include "TH1.h"
0049 #include "TTree.h"
0050 
0051 //
0052 // class declaration
0053 //
0054 
0055 using reco::TrackCollection;
0056 typedef std::vector<Trajectory> TrajectoryCollection;
0057 
0058 struct CPEBranch {
0059   float x, y, z;
0060   float distance, mdistance, shift, offsetA, offsetB, angle;
0061   float x1r, x2r, x1m, x2m, y1m, y2m;
0062 };
0063 
0064 struct TrackBranch {
0065   float px, py, pz, pt, eta, phi, charge;
0066 };
0067 
0068 struct GeometryBranch {
0069   int subdet, moduleGeometry, stereo, layer, side, ring;
0070   float pitch;
0071   int detid;
0072 };
0073 
0074 struct ClusterBranch {
0075   int strips[11];
0076   int size, firstStrip;
0077   float barycenter;
0078 };
0079 
0080 class SiStripCPEAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0081 public:
0082   explicit SiStripCPEAnalyzer(const edm::ParameterSet&);
0083   ~SiStripCPEAnalyzer() override = default;
0084 
0085   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0086 
0087 private:
0088   void analyze(const edm::Event&, const edm::EventSetup&) override;
0089   static bool goodMeasurement(TrajectoryMeasurement const& m);
0090 
0091   // ----------member data ---------------------------
0092 
0093   // tokens for event data
0094   const edm::EDGetTokenT<TrackCollection> tracksToken_;  //used to select what tracks to read from configuration file
0095   const edm::EDGetTokenT<TrajectoryCollection>
0096       trajsToken_;  //used to select what trajectories to read from configuration file
0097   const edm::EDGetTokenT<TrajTrackAssociationCollection> tjTagToken_;  //association map between tracks and trajectories
0098   const edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > clusToken_;  //clusters
0099 
0100   // tokens for ES data
0101   const edm::ESInputTag cpeTag_;
0102   const edm::ESGetToken<StripClusterParameterEstimator, TkStripCPERecord> cpeToken_;
0103   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0104   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0105 
0106   const StripClusterParameterEstimator* parameterestimator_;  //CPE
0107 
0108   // to fill the tree
0109   TTree* tree_;
0110   CPEBranch treeBranches_;
0111   TrackBranch trackBranches_;
0112   GeometryBranch geom1Branches_;
0113   GeometryBranch geom2Branches_;
0114   ClusterBranch cluster1Branches_;
0115   ClusterBranch cluster2Branches_;
0116 };
0117 
0118 //
0119 // constructors and destructor
0120 //
0121 SiStripCPEAnalyzer::SiStripCPEAnalyzer(const edm::ParameterSet& iConfig)
0122     : tracksToken_(consumes<TrackCollection>(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
0123       trajsToken_(consumes<TrajectoryCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trajectories"))),
0124       tjTagToken_(
0125           consumes<TrajTrackAssociationCollection>(iConfig.getUntrackedParameter<edm::InputTag>("association"))),
0126       clusToken_(
0127           consumes<edmNew::DetSetVector<SiStripCluster> >(iConfig.getUntrackedParameter<edm::InputTag>("clusters"))),
0128       cpeTag_(iConfig.getParameter<edm::ESInputTag>("StripCPE")),
0129       cpeToken_(esConsumes(cpeTag_)),
0130       topoToken_(esConsumes()),
0131       geomToken_(esConsumes()) {
0132   //now do what ever initialization is needed
0133   usesResource("TFileService");
0134   edm::Service<TFileService> fs;
0135   //histo = fs->make<TH1I>("charge" , "Charges" , 3 , -1 , 2 );
0136   tree_ = fs->make<TTree>("CPEanalysis", "CPE analysis tree");
0137   tree_->Branch(
0138       "Overlaps", &treeBranches_, "x:y:z:distance:mdistance:shift:offsetA:offsetB:angle:x1r:x2r:x1m:x2m:y1m:y2m");
0139   tree_->Branch("Tracks", &trackBranches_, "px:py:pz:pt:eta:phi:charge");
0140   tree_->Branch("Cluster1", &cluster1Branches_, "strips[11]/I:size/I:firstStrip/I:barycenter/F");
0141   tree_->Branch("Cluster2", &cluster2Branches_, "strips[11]/I:size/I:firstStrip/I:barycenter/F");
0142   tree_->Branch("Geom1", &geom1Branches_, "subdet/I:moduleGeometry/I:stereo/I:layer/I:side/I:ring/I:pitch/F:detid/I");
0143   tree_->Branch("Geom2", &geom2Branches_, "subdet/I:moduleGeometry/I:stereo/I:layer/I:side/I:ring/I:pitch/F:detid/I");
0144 }
0145 
0146 //
0147 // member functions
0148 //
0149 // ------------ method called for each event  ------------
0150 void SiStripCPEAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0151   using namespace edm;
0152 
0153   // load the CPE, geometry and topology
0154   parameterestimator_ = &iSetup.getData(cpeToken_);  //CPE
0155   const TrackerGeometry* tracker = &iSetup.getData(geomToken_);
0156   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0157 
0158   // prepare the output
0159   std::vector<SiStripOverlapHit> hitpairs;
0160 
0161   // loop on trajectories and tracks
0162   //for(const auto& tt : iEvent.get(tjTagToken_) ) {
0163   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociations;
0164   iEvent.getByToken(tjTagToken_, trajTrackAssociations);
0165   for (const auto& tt : *trajTrackAssociations) {
0166     auto& traj = *(tt.key);
0167     auto& track = *(tt.val);
0168     // track quantities
0169     trackBranches_.px = track.px();
0170     trackBranches_.py = track.py();
0171     trackBranches_.pz = track.pz();
0172     trackBranches_.pt = track.pt();
0173     trackBranches_.eta = track.eta();
0174     trackBranches_.phi = track.phi();
0175     trackBranches_.charge = track.charge();
0176     // loop on measurements
0177     for (auto it = traj.measurements().begin(); it != traj.measurements().end(); ++it) {
0178       auto& meas = *it;
0179 
0180       // only OT measurements on valid hits (not glued det)
0181       if (!goodMeasurement(meas))
0182         continue;
0183 
0184       // restrict the search for compatible hits in the same layer (measurements are sorted by layer)
0185       auto layerRangeEnd = it + 1;
0186       for (; layerRangeEnd < traj.measurements().end(); ++layerRangeEnd) {
0187         if (layerRangeEnd->layer()->seqNum() != meas.layer()->seqNum())
0188           break;
0189       }
0190 
0191       // now look for a matching hit in that range.
0192       auto meas2it = std::find_if(it + 1, layerRangeEnd, [meas](const TrajectoryMeasurement& m) -> bool {
0193         return goodMeasurement(m) && (m.recHit()->rawId() & 0x3) == (meas.recHit()->rawId() & 0x3);
0194       });
0195 
0196       // check if we found something, build a pair object and add to vector
0197       if (meas2it != layerRangeEnd) {
0198         auto& meas2 = *meas2it;
0199         hitpairs.push_back(SiStripOverlapHit(meas, meas2));
0200       }
0201     }
0202   }
0203 
0204   // load clusters
0205   Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
0206   iEvent.getByToken(clusToken_, clusters);
0207 
0208   // At this stage we will have a vector of pairs of measurement. Fill a ntuple and some histograms.
0209   for (const auto& pair : hitpairs) {
0210     //TODO basically everything below is done twice. -> can be factorized.
0211 
0212     // store generic information about the pair
0213     treeBranches_.x = pair.position().x();
0214     treeBranches_.y = pair.position().y();
0215     treeBranches_.z = pair.position().z();
0216     treeBranches_.distance = pair.distance(false);
0217     treeBranches_.mdistance = pair.distance(true);
0218     treeBranches_.shift = pair.shift();
0219     treeBranches_.offsetA = pair.offset(0);
0220     treeBranches_.offsetB = pair.offset(1);
0221     treeBranches_.angle = pair.getTrackLocalAngle(0);
0222     treeBranches_.x1r = pair.hitA()->localPosition().x();
0223     treeBranches_.x1m = pair.trajectoryStateOnSurface(0, false).localPosition().x();
0224     treeBranches_.y1m = pair.trajectoryStateOnSurface(0, false).localPosition().y();
0225     treeBranches_.x2r = pair.hitB()->localPosition().x();
0226     treeBranches_.x2m = pair.trajectoryStateOnSurface(1, false).localPosition().x();
0227     treeBranches_.y2m = pair.trajectoryStateOnSurface(1, false).localPosition().y();
0228 
0229     // load the clusters for the detectors
0230     auto detset1 = (*clusters)[pair.hitA()->rawId()];
0231     auto detset2 = (*clusters)[pair.hitB()->rawId()];
0232 
0233     // look for the proper cluster
0234     const GeomDetUnit* du;
0235     du = tracker->idToDetUnit(pair.hitA()->rawId());
0236     auto cluster1 = std::min_element(
0237         detset1.begin(), detset1.end(), [du, this](SiStripCluster const& cl1, SiStripCluster const& cl2) {
0238           return (fabs(parameterestimator_->localParameters(cl1, *du).first.x() - treeBranches_.x1r) <
0239                   fabs(parameterestimator_->localParameters(cl2, *du).first.x() - treeBranches_.x1r));
0240         });
0241     du = tracker->idToDetUnit(pair.hitB()->rawId());
0242     auto cluster2 = std::min_element(
0243         detset2.begin(), detset2.end(), [du, this](SiStripCluster const& cl1, SiStripCluster const& cl2) {
0244           return (fabs(parameterestimator_->localParameters(cl1, *du).first.x() - treeBranches_.x2r) <
0245                   fabs(parameterestimator_->localParameters(cl2, *du).first.x() - treeBranches_.x2r));
0246         });
0247 
0248     // store the amplitudes centered around the maximum
0249     auto amplitudes1 = cluster1->amplitudes();
0250     auto amplitudes2 = cluster2->amplitudes();
0251     auto max1 = std::max_element(amplitudes1.begin(), amplitudes1.end());
0252     auto max2 = std::max_element(amplitudes2.begin(), amplitudes2.end());
0253     for (unsigned int i = 0; i < 11; ++i) {
0254       cluster1Branches_.strips[i] = 0.;
0255       cluster2Branches_.strips[i] = 0.;
0256     }
0257     cluster1Branches_.size = amplitudes1.size();
0258     cluster2Branches_.size = amplitudes2.size();
0259     cluster1Branches_.firstStrip = cluster1->firstStrip();
0260     cluster1Branches_.barycenter = cluster1->barycenter();
0261     cluster2Branches_.firstStrip = cluster2->firstStrip();
0262     cluster2Branches_.barycenter = cluster2->barycenter();
0263 
0264     int cnt = 0;
0265     int offset = 5 - (max1 - amplitudes1.begin());
0266     for (auto& s : amplitudes1) {
0267       if ((offset + cnt) >= 0 && (offset + cnt) < 11) {
0268         cluster1Branches_.strips[offset + cnt] = s;
0269       }
0270       cnt++;
0271     }
0272     cnt = 0;
0273     offset = 5 - (max2 - amplitudes2.begin());
0274     for (auto& s : amplitudes2) {
0275       if ((offset + cnt) >= 0 && (offset + cnt) < 11) {
0276         cluster2Branches_.strips[offset + cnt] = s;
0277       }
0278       cnt++;
0279     }
0280 
0281     // store information about the geometry (for both sensors)
0282     DetId detid1 = pair.hitA()->geographicalId();
0283     DetId detid2 = pair.hitB()->geographicalId();
0284     geom1Branches_.detid = detid1.rawId();
0285     geom2Branches_.detid = detid2.rawId();
0286     geom1Branches_.subdet = detid1.subdetId();
0287     geom2Branches_.subdet = detid2.subdetId();
0288     //geom1Branches_.moduleGeometry = tTopo->moduleGeometry(detid1);
0289     //geom2Branches_.moduleGeometry = tTopo->moduleGeometry(detid2);
0290     geom1Branches_.stereo = tTopo->isStereo(detid1);
0291     geom2Branches_.stereo = tTopo->isStereo(detid2);
0292     geom1Branches_.layer = tTopo->layer(detid1);
0293     geom2Branches_.layer = tTopo->layer(detid2);
0294     geom1Branches_.side = tTopo->side(detid1);
0295     geom2Branches_.side = tTopo->side(detid2);
0296     if (geom1Branches_.subdet == StripSubdetector::TID) {
0297       geom1Branches_.ring = tTopo->tidRing(detid1);
0298     } else if (geom1Branches_.subdet == StripSubdetector::TEC) {
0299       geom1Branches_.ring = tTopo->tecRing(detid1);
0300     } else {
0301       geom1Branches_.ring = 0;
0302     }
0303     if (geom2Branches_.subdet == StripSubdetector::TID) {
0304       geom2Branches_.ring = tTopo->tidRing(detid2);
0305     } else if (geom2Branches_.subdet == StripSubdetector::TEC) {
0306       geom2Branches_.ring = tTopo->tecRing(detid2);
0307     } else {
0308       geom2Branches_.ring = 0;
0309     }
0310 
0311     const StripGeomDetUnit* theStripDet;
0312     theStripDet = dynamic_cast<const StripGeomDetUnit*>(tracker->idToDetUnit(pair.hitA()->rawId()));
0313     geom1Branches_.pitch =
0314         theStripDet->specificTopology().localPitch(pair.trajectoryStateOnSurface(0, false).localPosition());
0315     theStripDet = dynamic_cast<const StripGeomDetUnit*>(tracker->idToDetUnit(pair.hitB()->rawId()));
0316     geom2Branches_.pitch =
0317         theStripDet->specificTopology().localPitch(pair.trajectoryStateOnSurface(1, false).localPosition());
0318 
0319     //fill the tree.
0320     tree_->Fill();  // we fill one entry per overlap, to track info is multiplicated
0321   }
0322 }
0323 
0324 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0325 void SiStripCPEAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0326   edm::ParameterSetDescription desc;
0327   desc.addUntracked<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0328   desc.addUntracked<edm::InputTag>("trajectories", edm::InputTag("generalTracks"));
0329   desc.addUntracked<edm::InputTag>("association", edm::InputTag("generalTracks"));
0330   desc.addUntracked<edm::InputTag>("clusters", edm::InputTag("siStripClusters"));
0331   desc.add<edm::ESInputTag>("StripCPE", edm::ESInputTag("StripCPEfromTrackAngleESProducer", "StripCPEfromTrackAngle"));
0332   descriptions.addWithDefaultLabel(desc);
0333 }
0334 
0335 bool SiStripCPEAnalyzer::goodMeasurement(TrajectoryMeasurement const& m) {
0336   return m.recHit()->isValid() &&                        // valid
0337          m.recHit()->geographicalId().subdetId() > 2 &&  // not IT
0338          (m.recHit()->rawId() & 0x3) != 0 &&             // not glued DetLayer
0339          m.recHit()->getType() == 0;                     // only valid hits (redundant?)
0340 }
0341 
0342 //define this as a plug-in
0343 DEFINE_FWK_MODULE(SiStripCPEAnalyzer);
[ 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!