TrackAssociator/test/CaloMatchingExample.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackAssociator
0004 // Class:      CaloMatchingExample
0005 //
0006 /*
0007
0008  Description: Example shows how to access various forms of energy deposition and store them in an ntuple
0009
0010 */
0011 //
0012 // Original Author:  Dmytro Kovalskyi
0013 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0014 //
0015 //
0016
0017 // system include files
0018 #include <memory>
0019
0020 // user include files
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "DataFormats/Common/interface/OrphanHandle.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0035 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0036 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0037 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0038 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0039 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0040
0041 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0042 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0043 #include "SimDataFormats/Track/interface/SimTrack.h"
0044 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0045 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0046 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0047 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0049
0050 // calorimeter info
0051 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0052 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0053 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0054 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0055 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0056 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0057
0058 #include "Geometry/DTGeometry/interface/DTLayer.h"
0059 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0060 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0061
0062 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0063 #include "DataFormats/GeometrySurface/interface/Plane.h"
0064
0065 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0066
0067 #include "MagneticField/Engine/interface/MagneticField.h"
0068 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0069
0070 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0071
0072 #include <boost/regex.hpp>
0073
0074 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0075 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0076
0077 #include "TFile.h"
0078 #include "TTree.h"
0079
0080 #include "CLHEP/Random/Random.h"
0081
0082 class CaloMatchingExample : public edm::one::EDAnalyzer<> {
0083 public:
0084   explicit CaloMatchingExample(const edm::ParameterSet&);
0085   virtual ~CaloMatchingExample() {
0086     file_->cd();
0087     tree_->Write();
0088     file_->Close();
0089   }
0090
0091   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0092
0093 private:
0094   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0095
0096   TrackDetectorAssociator trackAssociator_;
0097   bool useEcal_;
0098   bool useHcal_;
0099   bool useMuon_;
0100   bool useOldMuonMatching_;
0101   TFile* file_;
0102   TTree* tree_;
0103   int nTracks_;
0104
0105   float ecalCrossedEnergy_[1000];
0106   float ecal3x3Energy_[1000];
0107   float ecal5x5Energy_[1000];
0108   float ecal3x3EnergyMax_[1000];
0109   float ecal5x5EnergyMax_[1000];
0110   float ecalTrueEnergy_[1000];
0111   float trkPosAtEcal_[1000][2];
0112   float ecalMaxPos_[1000][2];
0113
0114   float ecalCrossedEnergyRandom_[1000];
0115   float ecal3x3EnergyRandom_[1000];
0116   float ecal5x5EnergyRandom_[1000];
0117   float ecal3x3EnergyMaxRandom_[1000];
0118   float ecal5x5EnergyMaxRandom_[1000];
0119   float trkPosAtEcalRandom_[1000][2];
0120   float ecalMaxPosRandom_[1000][2];
0121
0122   float hcalCrossedEnergy_[1000];
0123   float hcal3x3Energy_[1000];
0124   float hcal5x5Energy_[1000];
0125   float hcal3x3EnergyMax_[1000];
0126   float hcal5x5EnergyMax_[1000];
0127   float hcalTrueEnergy_[1000];
0128   float hcalTrueEnergyCorrected_[1000];
0129   float trkPosAtHcal_[1000][2];
0130   float hcalMaxPos_[1000][2];
0131   float trackPt_[1000];
0132
0133   float hcalCrossedEnergyRandom_[1000];
0134   float hcal3x3EnergyRandom_[1000];
0135   float hcal5x5EnergyRandom_[1000];
0136   float hcal3x3EnergyMaxRandom_[1000];
0137   float hcal5x5EnergyMaxRandom_[1000];
0138   float trkPosAtHcalRandom_[1000][2];
0139   float hcalMaxPosRandom_[1000][2];
0140
0141   edm::InputTag inputRecoTrackColl_;
0142   TrackAssociatorParameters parameters_;
0143 };
0144
0145 CaloMatchingExample::CaloMatchingExample(const edm::ParameterSet& iConfig) : caloGeomToken_(esConsumes()) {
0146   CLHEP::HepRandom::createInstance();
0147   file_ = new TFile(iConfig.getParameter<std::string>("outputfile").c_str(), "recreate");
0148   tree_ = new TTree("calomatch", "calomatch");
0149   tree_->Branch("nTracks", &nTracks_, "nTracks/I");
0150
0151   tree_->Branch("ecalCrossedEnergy", ecalCrossedEnergy_, "ecalCrossedEnergy[nTracks]/F");
0152   tree_->Branch("ecal3x3Energy", ecal3x3Energy_, "ecal3x3Energy[nTracks]/F");
0153   tree_->Branch("ecal5x5Energy", ecal5x5Energy_, "ecal5x5Energy[nTracks]/F");
0154   tree_->Branch("ecal3x3EnergyMax", ecal3x3EnergyMax_, "ecal3x3EnergyMax[nTracks]/F");
0155   tree_->Branch("ecal5x5EnergyMax", ecal5x5EnergyMax_, "ecal5x5EnergyMax[nTracks]/F");
0156   tree_->Branch("ecalTrueEnergy", ecalTrueEnergy_, "ecalTrueEnergy[nTracks]/F");
0157   tree_->Branch("trkPosAtEcal", trkPosAtEcal_, "trkPosAtEcal[nTracks][2]/F");
0158   tree_->Branch("ecalMaxPos", ecalMaxPos_, "ecalMaxPos[nTracks][2]/F");
0159
0160   tree_->Branch("ecalCrossedEnergyRandom", ecalCrossedEnergyRandom_, "ecalCrossedEnergyRandom[nTracks]/F");
0161   tree_->Branch("ecal3x3EnergyRandom", ecal3x3EnergyRandom_, "ecal3x3EnergyRandom[nTracks]/F");
0162   tree_->Branch("ecal5x5EnergyRandom", ecal5x5EnergyRandom_, "ecal5x5EnergyRandom[nTracks]/F");
0163   tree_->Branch("ecal3x3EnergyMaxRandom", ecal3x3EnergyMaxRandom_, "ecal3x3EnergyMaxRandom[nTracks]/F");
0164   tree_->Branch("ecal5x5EnergyMaxRandom", ecal5x5EnergyMaxRandom_, "ecal5x5EnergyMaxRandom[nTracks]/F");
0165   tree_->Branch("trkPosAtEcalRandom", trkPosAtEcalRandom_, "trkPosAtEcalRandom[nTracks][2]/F");
0166   tree_->Branch("ecalMaxPosRandom", ecalMaxPosRandom_, "ecalMaxPosRandom[nTracks][2]/F");
0167
0168   tree_->Branch("hcalCrossedEnergy", hcalCrossedEnergy_, "hcalCrossedEnergy[nTracks]/F");
0169   tree_->Branch("hcal3x3Energy", hcal3x3Energy_, "hcal3x3Energy[nTracks]/F");
0170   tree_->Branch("hcal5x5Energy", hcal5x5Energy_, "hcal5x5Energy[nTracks]/F");
0171   tree_->Branch("hcal3x3EnergyMax", hcal3x3EnergyMax_, "hcal3x3EnergyMax[nTracks]/F");
0172   tree_->Branch("hcal5x5EnergyMax", hcal5x5EnergyMax_, "hcal5x5EnergyMax[nTracks]/F");
0173   tree_->Branch("hcalTrueEnergy", hcalTrueEnergy_, "hcalTrueEnergy[nTracks]/F");
0174   tree_->Branch("hcalTrueEnergyCorrected", hcalTrueEnergyCorrected_, "hcalTrueEnergyCorrected[nTracks]/F");
0175   tree_->Branch("trkPosAtHcal", trkPosAtHcal_, "trkPosAtHcal[nTracks][2]/F");
0176   tree_->Branch("hcalMaxPos", hcalMaxPos_, "hcalMaxPos[nTracks][2]/F");
0177
0178   tree_->Branch("hcalCrossedEnergyRandom", hcalCrossedEnergyRandom_, "hcalCrossedEnergyRandom[nTracks]/F");
0179   tree_->Branch("hcal3x3EnergyRandom", hcal3x3EnergyRandom_, "hcal3x3EnergyRandom[nTracks]/F");
0180   tree_->Branch("hcal5x5EnergyRandom", hcal5x5EnergyRandom_, "hcal5x5EnergyRandom[nTracks]/F");
0181   tree_->Branch("hcal3x3EnergyMaxRandom", hcal3x3EnergyMaxRandom_, "hcal3x3EnergyMaxRandom[nTracks]/F");
0182   tree_->Branch("hcal5x5EnergyMaxRandom", hcal5x5EnergyMaxRandom_, "hcal5x5EnergyMaxRandom[nTracks]/F");
0183   tree_->Branch("trkPosAtHcalRandom", trkPosAtHcalRandom_, "trkPosAtHcalRandom[nTracks][2]/F");
0184   tree_->Branch("hcalMaxPosRandom", hcalMaxPosRandom_, "hcalMaxPosRandom[nTracks][2]/F");
0185
0186   tree_->Branch("trackPt", trackPt_, "trackPt_[nTracks]/F");
0187
0188   inputRecoTrackColl_ = iConfig.getParameter<edm::InputTag>("inputRecoTrackColl");
0189
0190   // TrackAssociator parameters
0191   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0192   edm::ConsumesCollector iC = consumesCollector();
0193   parameters_.loadParameters(parameters, iC);
0194
0195   trackAssociator_.useDefaultPropagator();
0196 }
0197
0198 void CaloMatchingExample::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0199   using namespace edm;
0200
0201   // calo geometry
0202   const auto& geometry = iSetup.getHandle(caloGeomToken_);
0203   if (!geometry.isValid())
0204     throw cms::Exception("FatalError") << "Unable to find CaloGeometryRecord in event!\n";
0205
0206   nTracks_ = 0;
0207
0208   // get reco tracks
0209   Handle<reco::TrackCollection> recoTracks;
0210   iEvent.getByLabel(inputRecoTrackColl_, recoTracks);
0211   if (!recoTracks.isValid())
0212     throw cms::Exception("FatalError") << "No reco tracks were found\n";
0213
0214   // loop over reconstructed tracks
0215   LogTrace("TrackAssociator") << "Number of reco tracks found in the event: " << recoTracks->size();
0216   // for(SimTrackContainer::const_iterator tracksCI = simTracks->begin();
0217   //    tracksCI != simTracks->end(); tracksCI++){
0218   for (reco::TrackCollection::const_iterator recoTrack = recoTracks->begin(); recoTrack != recoTracks->end();
0219        ++recoTrack) {
0220     // skip low Pt tracks
0221     if (recoTrack->pt() < 2) {
0222       LogTrace("TrackAssociator") << "Skipped low Pt track (Pt: " << recoTrack->pt() << ")";
0223       continue;
0224     }
0225
0226     LogTrace("TrackAssociator") << "\n-------------------------------------------------------\n Track (pt,eta,phi): "
0227                                 << recoTrack->pt() << " , " << recoTrack->eta() << " , " << recoTrack->phi();
0228
0229     // Get track matching info
0230
0231     LogTrace("TrackAssociator")
0232         << "===========================================================================\nDetails:\n";
0233     TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, *recoTrack, parameters_);
0234     // get some noise info (random direction)
0235     ROOT::Math::RhoEtaPhiVector randomVector(10,
0236                                              (CLHEP::HepRandom::getTheEngine()->flat() - 0.5) * 6,
0237                                              (CLHEP::HepRandom::getTheEngine()->flat() - 0.5) * 2 * M_PI);
0238     TrackDetMatchInfo infoRandom =
0239         trackAssociator_.associate(iEvent,
0240                                    iSetup,
0241                                    GlobalVector(randomVector.x(), randomVector.y(), randomVector.z()),
0242                                    GlobalPoint(0, 0, 0),
0243                                    +1,
0244                                    parameters_);
0245
0246     ///////////////////////////////////////////////////
0247     //
0248     //   Fill ntuple
0249     //
0250     ///////////////////////////////////////////////////
0251
0252     DetId centerId;
0253     DetId centerIdRandom;
0254
0255     trackPt_[nTracks_] = recoTrack->pt();
0256     ecalCrossedEnergy_[nTracks_] = info.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
0257     centerId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits);
0258     ecal3x3Energy_[nTracks_] = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
0259     ecal5x5Energy_[nTracks_] = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
0260     ecal3x3EnergyMax_[nTracks_] = info.nXnEnergy(centerId, TrackDetMatchInfo::EcalRecHits, 1);
0261     ecal5x5EnergyMax_[nTracks_] = info.nXnEnergy(centerId, TrackDetMatchInfo::EcalRecHits, 2);
0262     ecalTrueEnergy_[nTracks_] = info.ecalTrueEnergy;
0263     trkPosAtEcal_[nTracks_][0] = info.trkGlobPosAtEcal.eta();
0264     trkPosAtEcal_[nTracks_][1] = info.trkGlobPosAtEcal.phi();
0265     if (geometry->getSubdetectorGeometry(centerId) &&
0266         geometry->getSubdetectorGeometry(centerId)->getGeometry(centerId)) {
0267       GlobalPoint position = geometry->getSubdetectorGeometry(centerId)->getGeometry(centerId)->getPosition();
0268       ecalMaxPos_[nTracks_][0] = position.eta();
0269       ecalMaxPos_[nTracks_][1] = position.phi();
0270     } else {
0271       ecalMaxPos_[nTracks_][0] = -999;
0272       ecalMaxPos_[nTracks_][1] = -999;
0273     }
0274     ecalCrossedEnergyRandom_[nTracks_] = infoRandom.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
0275     centerIdRandom = infoRandom.findMaxDeposition(TrackDetMatchInfo::EcalRecHits);
0276     ecal3x3EnergyRandom_[nTracks_] = infoRandom.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
0277     ecal5x5EnergyRandom_[nTracks_] = infoRandom.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
0278     ecal3x3EnergyMaxRandom_[nTracks_] = infoRandom.nXnEnergy(centerIdRandom, TrackDetMatchInfo::EcalRecHits, 1);
0279     ecal5x5EnergyMaxRandom_[nTracks_] = infoRandom.nXnEnergy(centerIdRandom, TrackDetMatchInfo::EcalRecHits, 2);
0280     trkPosAtEcalRandom_[nTracks_][0] = infoRandom.trkGlobPosAtEcal.eta();
0281     trkPosAtEcalRandom_[nTracks_][1] = infoRandom.trkGlobPosAtEcal.phi();
0282
0283     hcalCrossedEnergy_[nTracks_] = info.hcalEnergy();
0284     centerId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
0285     hcal3x3Energy_[nTracks_] = info.nXnEnergy(centerId, TrackDetMatchInfo::HcalRecHits, 1);
0286     hcal5x5Energy_[nTracks_] = info.nXnEnergy(centerId, TrackDetMatchInfo::HcalRecHits, 2);
0287     hcalTrueEnergy_[nTracks_] = info.hcalTrueEnergy;
0288     hcalTrueEnergyCorrected_[nTracks_] = info.hcalTrueEnergyCorrected;
0289     trkPosAtHcal_[nTracks_][0] = info.trkGlobPosAtHcal.eta();
0290     trkPosAtHcal_[nTracks_][1] = info.trkGlobPosAtHcal.phi();
0291     if (geometry->getSubdetectorGeometry(centerId) &&
0292         geometry->getSubdetectorGeometry(centerId)->getGeometry(centerId)) {
0293       GlobalPoint position = geometry->getSubdetectorGeometry(centerId)->getGeometry(centerId)->getPosition();
0294       hcalMaxPos_[nTracks_][0] = position.eta();
0295       hcalMaxPos_[nTracks_][1] = position.phi();
0296     } else {
0297       hcalMaxPos_[nTracks_][0] = -999;
0298       hcalMaxPos_[nTracks_][1] = -999;
0299     }
0300     hcalCrossedEnergyRandom_[nTracks_] = infoRandom.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
0301     centerIdRandom = infoRandom.findMaxDeposition(TrackDetMatchInfo::HcalRecHits);
0302     hcal3x3EnergyRandom_[nTracks_] = infoRandom.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
0303     hcal5x5EnergyRandom_[nTracks_] = infoRandom.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2);
0304     hcal3x3EnergyMaxRandom_[nTracks_] = infoRandom.nXnEnergy(centerIdRandom, TrackDetMatchInfo::HcalRecHits, 1);
0305     hcal5x5EnergyMaxRandom_[nTracks_] = infoRandom.nXnEnergy(centerIdRandom, TrackDetMatchInfo::HcalRecHits, 2);
0306     trkPosAtHcalRandom_[nTracks_][0] = infoRandom.trkGlobPosAtHcal.eta();
0307     trkPosAtHcalRandom_[nTracks_][1] = infoRandom.trkGlobPosAtHcal.phi();
0308
0309     // Debugging information
0310
0311     LogTrace("TrackAssociator") << "ECAL, number of crossed cells: " << info.crossedEcalRecHits.size();
0312     LogTrace("TrackAssociator") << "ECAL, energy of crossed cells: " << info.ecalEnergy() << " GeV";
0313     LogTrace("TrackAssociator") << "ECAL, number of cells in the cone: " << info.ecalRecHits.size();
0314     LogTrace("TrackAssociator") << "ECAL, energy in the cone: " << info.ecalConeEnergy() << " GeV";
0315     LogTrace("TrackAssociator") << "ECAL, true energy: " << info.ecalTrueEnergy << " GeV";
0316     LogTrace("TrackAssociator") << "ECAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtEcal.z() << ", "
0317                                 << info.trkGlobPosAtEcal.R() << " , " << info.trkGlobPosAtEcal.eta() << " , "
0318                                 << info.trkGlobPosAtEcal.phi();
0319
0320     LogTrace("TrackAssociator") << "HCAL, number of crossed elements (towers): " << info.crossedTowers.size();
0321     LogTrace("TrackAssociator") << "HCAL, energy of crossed elements (towers): " << info.hcalTowerEnergy() << " GeV";
0322     LogTrace("TrackAssociator") << "HCAL, number of crossed elements (hits): " << info.crossedHcalRecHits.size();
0323     LogTrace("TrackAssociator") << "HCAL, energy of crossed elements (hits): " << info.hcalEnergy() << " GeV";
0324     LogTrace("TrackAssociator") << "HCAL, number of elements in the cone (towers): " << info.towers.size();
0325     LogTrace("TrackAssociator") << "HCAL, energy in the cone (towers): " << info.hcalTowerConeEnergy() << " GeV";
0326     LogTrace("TrackAssociator") << "HCAL, number of elements in the cone (hits): " << info.hcalRecHits.size();
0327     LogTrace("TrackAssociator") << "HCAL, energy in the cone (hits): " << info.hcalConeEnergy() << " GeV";
0328     LogTrace("TrackAssociator") << "HCAL, true energy: " << info.hcalTrueEnergy << " GeV";
0329     LogTrace("TrackAssociator") << "HCAL, true energy corrected: " << info.hcalTrueEnergyCorrected << " GeV";
0330     LogTrace("TrackAssociator") << "HCAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtHcal.z() << ", "
0331                                 << info.trkGlobPosAtHcal.R() << " , " << info.trkGlobPosAtHcal.eta() << " , "
0332                                 << info.trkGlobPosAtHcal.phi();
0333
0334     LogTrace("TrackAssociator") << "HO, number of crossed elements (hits): " << info.crossedHORecHits.size();
0335     LogTrace("TrackAssociator") << "HO, energy of crossed elements (hits): " << info.hoEnergy() << " GeV";
0336     LogTrace("TrackAssociator") << "HO, number of elements in the cone: " << info.hoRecHits.size();
0337     LogTrace("TrackAssociator") << "HO, energy in the cone: " << info.hoConeEnergy() << " GeV";
0338     LogTrace("TrackAssociator") << "HCAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtHO.z() << ", "
0339                                 << info.trkGlobPosAtHO.R() << " , " << info.trkGlobPosAtHO.eta() << " , "
0340                                 << info.trkGlobPosAtHO.phi();
0341     nTracks_++;
0342   }
0343   tree_->Fill();
0344 }
0345
0346 //define this as a plug-in
0347 DEFINE_FWK_MODULE(CaloMatchingExample);