Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:38

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackAssociator
0004 // Class:      TestTrackAssociator
0005 //
0006 /*
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Dmytro Kovalskyi
0015 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025 
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "DataFormats/Common/interface/Handle.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "DataFormats/Common/interface/OrphanHandle.h"
0032 
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 
0035 #include "DataFormats/TrackReco/interface/Track.h"
0036 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0037 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0038 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0039 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0040 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0041 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0042 
0043 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0044 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0045 #include "SimDataFormats/Track/interface/SimTrack.h"
0046 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0047 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0048 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0049 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0050 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0051 
0052 // calorimeter info
0053 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0054 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0055 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0056 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0057 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0058 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0059 
0060 #include "Geometry/DTGeometry/interface/DTLayer.h"
0061 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0062 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0063 
0064 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0065 #include "DataFormats/GeometrySurface/interface/Plane.h"
0066 
0067 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0068 
0069 #include "MagneticField/Engine/interface/MagneticField.h"
0070 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0071 
0072 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0073 
0074 #include <boost/regex.hpp>
0075 
0076 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0077 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0078 
0079 #include "DataFormats/MuonReco/interface/Muon.h"
0080 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0081 
0082 class TestTrackAssociator : public edm::one::EDAnalyzer<> {
0083 public:
0084   explicit TestTrackAssociator(const edm::ParameterSet&);
0085   virtual ~TestTrackAssociator() {}
0086 
0087   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0088 
0089 private:
0090   TrackDetectorAssociator trackAssociator_;
0091   TrackAssociatorParameters parameters_;
0092 };
0093 
0094 TestTrackAssociator::TestTrackAssociator(const edm::ParameterSet& iConfig) {
0095   // TrackAssociator parameters
0096   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0097   edm::ConsumesCollector iC = consumesCollector();
0098   parameters_.loadParameters(parameters, iC);
0099 
0100   trackAssociator_.useDefaultPropagator();
0101 }
0102 
0103 void TestTrackAssociator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0104   using namespace edm;
0105 
0106   // get list of tracks and their vertices
0107   Handle<reco::MuonCollection> muons;
0108   iEvent.getByLabel("muons", muons);
0109 
0110   // loop
0111   LogVerbatim("TrackAssociator") << "Number of muons found in the event: " << muons->size();
0112   for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0113     // skip low Pt tracks
0114     if (muon->pt() < 2) {
0115       LogVerbatim("TrackAssociator") << "Skipped low Pt muon (Pt: " << muon->pt() << ")";
0116       continue;
0117     }
0118 
0119     // skip tracks originated away from the IP
0120     if (fabs(muon->vertex().rho()) > 50) {
0121       LogVerbatim("TrackAssociator") << "Skipped track with large impact parameter: " << muon->vertex().rho();
0122       continue;
0123     }
0124 
0125     LogVerbatim("TrackAssociator") << "\n-------------------------------------------------------\n Track (pt,eta,phi): "
0126                                    << muon->pt() << " , " << muon->eta() << " , " << muon->phi();
0127 
0128     TrackDetMatchInfo info;
0129     if (muon->innerTrack().isAvailable())
0130       info = trackAssociator_.associate(iEvent, iSetup, *muon->innerTrack(), parameters_);
0131     else {
0132       if (!muon->outerTrack().isAvailable()) {
0133         LogVerbatim("TrackAssociator") << "No refernced tracks are available, skim the muon";
0134         continue;
0135       }
0136       info = trackAssociator_.associate(iEvent, iSetup, *muon->outerTrack(), parameters_);
0137     }
0138 
0139     LogVerbatim("TrackAssociator") << "===========================================================================";
0140     LogVerbatim("TrackAssociator")
0141         << "ECAL RecHit energy: crossed, 3x3(max), 5x5(max), 3x3(direction), 5x5(direction), cone R0.5, generator";
0142     DetId centerId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits);
0143     LogVerbatim("TrackAssociator") << "     " << info.crossedEnergy(TrackDetMatchInfo::EcalRecHits) << ", \t"
0144                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::EcalRecHits, 1) << ", \t"
0145                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::EcalRecHits, 2) << ", \t"
0146                                    << info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1) << ", \t"
0147                                    << info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2) << ", \t"
0148                                    << info.coneEnergy(0.5, TrackDetMatchInfo::EcalRecHits) << ", \t"
0149                                    << info.ecalTrueEnergy;
0150     LogVerbatim("TrackAssociator") << "ECAL trajectory point (z,Rho,eta,phi), max deposit DetId";
0151     LogVerbatim("TrackAssociator") << "     "
0152                                    << "(" << info.trkGlobPosAtEcal.z() << ", " << info.trkGlobPosAtEcal.Rho() << ", "
0153                                    << info.trkGlobPosAtEcal.eta() << ", " << info.trkGlobPosAtEcal.phi() << "), "
0154                                    << centerId.rawId();
0155     LogVerbatim("TrackAssociator") << "ECAL crossed DetIds with associated hits: (id, energy, z, perp, eta, phi)";
0156     for (std::vector<const EcalRecHit*>::const_iterator hit = info.crossedEcalRecHits.begin();
0157          hit != info.crossedEcalRecHits.end();
0158          ++hit) {
0159       GlobalPoint point = info.getPosition((*hit)->detid());
0160       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0161                                      << point.z() << ", \t" << point.perp() << ", \t" << point.eta() << ", \t"
0162                                      << point.phi() << ")";
0163     }
0164     LogVerbatim("TrackAssociator") << "ECAL crossed DetIds: (id, z, perp, eta, phi)";
0165     for (std::vector<DetId>::const_iterator id = info.crossedEcalIds.begin(); id != info.crossedEcalIds.end(); ++id) {
0166       GlobalPoint point = info.getPosition(*id);
0167       LogVerbatim("TrackAssociator") << "\t" << id->rawId() << " \t(" << point.z() << ", \t" << point.perp() << ", \t"
0168                                      << point.eta() << ", \t" << point.phi() << ")";
0169     }
0170     LogVerbatim("TrackAssociator") << "ECAL associated DetIds: (id, energy, z, perp, eta, phi)";
0171     for (std::vector<const EcalRecHit*>::const_iterator hit = info.ecalRecHits.begin(); hit != info.ecalRecHits.end();
0172          ++hit) {
0173       GlobalPoint point = info.getPosition((*hit)->detid());
0174       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0175                                      << point.z() << ", \t" << point.perp() << ", \t" << point.eta() << ", \t"
0176                                      << point.phi() << ")";
0177     }
0178     LogVerbatim("TrackAssociator") << "Preshower crossed DetIds: (id, z, perp, eta, phi)";
0179     for (std::vector<DetId>::const_iterator id = info.crossedPreshowerIds.begin(); id != info.crossedPreshowerIds.end();
0180          ++id) {
0181       GlobalPoint point = info.getPosition(*id);
0182       LogVerbatim("TrackAssociator") << "\t" << id->rawId() << " \t(" << point.z() << ", \t" << point.perp() << ", \t"
0183                                      << point.eta() << ", \t" << point.phi() << ")";
0184     }
0185 
0186     LogVerbatim("TrackAssociator") << "---------------------------------------------------------------------------";
0187     LogVerbatim("TrackAssociator")
0188         << "HCAL RecHit energy: crossed, 3x3(max), 5x5(max), 3x3(direction), 5x5(direction), cone R0.5, generator";
0189     centerId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
0190     LogVerbatim("TrackAssociator") << "     " << info.crossedEnergy(TrackDetMatchInfo::HcalRecHits) << ", \t"
0191                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::HcalRecHits, 1) << ", \t"
0192                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::HcalRecHits, 2) << ", \t"
0193                                    << info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1) << ", \t"
0194                                    << info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2) << ", \t"
0195                                    << info.coneEnergy(0.5, TrackDetMatchInfo::HcalRecHits) << ", \t"
0196                                    << info.hcalTrueEnergyCorrected;
0197     LogVerbatim("TrackAssociator") << "HCAL trajectory point (z,Rho,eta,phi), max deposit DetId";
0198     LogVerbatim("TrackAssociator") << "     "
0199                                    << "(" << info.trkGlobPosAtHcal.z() << ", " << info.trkGlobPosAtHcal.Rho() << ", "
0200                                    << info.trkGlobPosAtHcal.eta() << ", " << info.trkGlobPosAtHcal.phi() << "), "
0201                                    << centerId.rawId();
0202     LogVerbatim("TrackAssociator") << "HCAL crossed DetIds with hits:";
0203     for (std::vector<const HBHERecHit*>::const_iterator hit = info.crossedHcalRecHits.begin();
0204          hit != info.crossedHcalRecHits.end();
0205          ++hit) {
0206       GlobalPoint point = info.getPosition((*hit)->detid());
0207       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0208                                      << point.z() << ", \t" << point.perp() << ", \t" << (*hit)->id().depth() << ", \t"
0209                                      << point.eta() << ", \t" << point.phi() << ")";
0210     }
0211     LogVerbatim("TrackAssociator") << "HCAL crossed DetIds: (id, z, perp, eta, phi)";
0212     for (std::vector<DetId>::const_iterator id = info.crossedHcalIds.begin(); id != info.crossedHcalIds.end(); ++id) {
0213       GlobalPoint point = info.getPosition(*id);
0214       LogVerbatim("TrackAssociator") << "\t" << id->rawId() << " \t(" << point.z() << ", \t" << point.perp() << ", \t"
0215                                      << point.eta() << ", \t" << point.phi() << ")";
0216     }
0217     LogVerbatim("TrackAssociator") << "HCAL associated DetIds: id, (energy, z, perp, depth, eta, phi)";
0218     for (std::vector<const HBHERecHit*>::const_iterator hit = info.hcalRecHits.begin(); hit != info.hcalRecHits.end();
0219          ++hit) {
0220       GlobalPoint point = info.getPosition((*hit)->detid());
0221       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0222                                      << point.z() << ", \t" << point.perp() << ", \t" << (*hit)->id().depth() << ", \t"
0223                                      << point.eta() << ", \t" << point.phi() << ")";
0224     }
0225 
0226     LogVerbatim("TrackAssociator") << "---------------------------------------------------------------------------";
0227     LogVerbatim("TrackAssociator")
0228         << "HO RecHit energy: crossed, 3x3(max), 5x5(max), 3x3(direction), 5x5(direction), cone R0.5";
0229     centerId = info.findMaxDeposition(TrackDetMatchInfo::HORecHits);
0230     LogVerbatim("TrackAssociator") << "     " << info.crossedEnergy(TrackDetMatchInfo::HORecHits) << ", \t"
0231                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::HORecHits, 1) << ", \t"
0232                                    << info.nXnEnergy(centerId, TrackDetMatchInfo::HORecHits, 2) << ", \t"
0233                                    << info.nXnEnergy(TrackDetMatchInfo::HORecHits, 1) << ", \t"
0234                                    << info.nXnEnergy(TrackDetMatchInfo::HORecHits, 2) << ", \t"
0235                                    << info.coneEnergy(0.5, TrackDetMatchInfo::HORecHits);
0236     LogVerbatim("TrackAssociator") << "HO trajectory point (z,Rho,eta,phi), max deposit DetId";
0237     LogVerbatim("TrackAssociator") << "     "
0238                                    << "(" << info.trkGlobPosAtHO.z() << ", " << info.trkGlobPosAtHO.Rho() << ", "
0239                                    << info.trkGlobPosAtHO.eta() << ", " << info.trkGlobPosAtHO.phi() << "), "
0240                                    << centerId.rawId();
0241     LogVerbatim("TrackAssociator") << "HO crossed DetIds with hits:";
0242     for (std::vector<const HORecHit*>::const_iterator hit = info.crossedHORecHits.begin();
0243          hit != info.crossedHORecHits.end();
0244          ++hit) {
0245       GlobalPoint point = info.getPosition((*hit)->detid());
0246       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0247                                      << point.z() << ", \t" << point.perp() << ", \t" << point.eta() << ", \t"
0248                                      << point.phi() << ")";
0249     }
0250     LogVerbatim("TrackAssociator") << "HO crossed DetIds: (id, z, perp, eta, phi)";
0251     for (std::vector<DetId>::const_iterator id = info.crossedHOIds.begin(); id != info.crossedHOIds.end(); ++id) {
0252       GlobalPoint point = info.getPosition(*id);
0253       LogVerbatim("TrackAssociator") << "\t" << id->rawId() << " \t(" << point.z() << ", \t" << point.perp() << ", \t"
0254                                      << point.eta() << ", \t" << point.phi() << ")";
0255     }
0256     LogVerbatim("TrackAssociator") << "HO associated DetIds: (id, energy,position)";
0257     for (std::vector<const HORecHit*>::const_iterator hit = info.hoRecHits.begin(); hit != info.hoRecHits.end();
0258          ++hit) {
0259       GlobalPoint point = info.getPosition((*hit)->detid());
0260       LogVerbatim("TrackAssociator") << "\t" << (*hit)->detid().rawId() << ", " << (*hit)->energy() << " \t("
0261                                      << point.z() << ", \t" << point.perp() << ", \t" << point.eta() << ", \t"
0262                                      << point.phi();
0263       // << ") ## " << info.dumpGeometry((*hit)->detid());
0264     }
0265 
0266     if (parameters_.useMuon) {
0267       LogVerbatim("TrackAssociator") << "Muon detector matching details: ";
0268       for (std::vector<TAMuonChamberMatch>::const_iterator chamber = info.chambers.begin();
0269            chamber != info.chambers.end();
0270            chamber++) {
0271         LogVerbatim("TrackAssociator") << chamber->info()
0272                                        << "\n\t(DetId, station, edgeX, edgeY): " << chamber->id.rawId() << ", "
0273                                        << chamber->station() << ", " << chamber->localDistanceX << ", "
0274                                        << chamber->localDistanceY << ", ";
0275         LogVerbatim("TrackAssociator") << "\t trajectory global point (z,perp,eta,phi): "
0276                                        << chamber->tState.globalPosition().z() << ", "
0277                                        << chamber->tState.globalPosition().perp() << ", "
0278                                        << chamber->tState.globalPosition().eta() << ", "
0279                                        << chamber->tState.globalPosition().phi();
0280         LogVerbatim("TrackAssociator") << "\t trajectory local point (x,y): " << chamber->tState.localPosition().x()
0281                                        << ", " << chamber->tState.localPosition().y();
0282 
0283         for (std::vector<TAMuonSegmentMatch>::const_iterator segment = chamber->segments.begin();
0284              segment != chamber->segments.end();
0285              segment++) {
0286           LogVerbatim("TrackAssociator") << "\t segment position (z,Rho,eta,phi,DetId): "
0287                                          << segment->segmentGlobalPosition.z() << ", "
0288                                          << segment->segmentGlobalPosition.Rho() << ", "
0289                                          << segment->segmentGlobalPosition.eta() << ", "
0290                                          << segment->segmentGlobalPosition.phi() << ", " << chamber->id.rawId();
0291           LogVerbatim("TrackAssociator") << "\t segment local position (x,y): " << segment->segmentLocalPosition.x()
0292                                          << ", " << segment->segmentLocalPosition.y();
0293         }
0294       }
0295     }
0296   }
0297 }
0298 
0299 //define this as a plug-in
0300 DEFINE_FWK_MODULE(TestTrackAssociator);