File indexing completed on 2021-11-22 06:26:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/EDMException.h"
0020 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0023 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0024 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0025
0026
0027 #include "TH1.h"
0028 #include "TH2.h"
0029 #include "TFile.h"
0030
0031
0032 #include <iostream>
0033
0034 class SimAnalyzer : public edm::one::EDAnalyzer<> {
0035 public:
0036
0037 explicit SimAnalyzer(edm::ParameterSet const &theConf);
0038
0039 ~SimAnalyzer();
0040
0041
0042 virtual void analyze(edm::Event const &theEvent, edm::EventSetup const &theSetup);
0043
0044 virtual void beginJob();
0045
0046 private:
0047
0048 double angle(double theAngle);
0049
0050 void closeRootFile();
0051
0052 void initHistograms();
0053
0054
0055 void trackerStatistics(edm::Event const &theEvent, edm::EventSetup const &theSetup);
0056
0057 private:
0058 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0059 int theEvents;
0060 int theDebugLevel;
0061 double theSearchPhiTIB;
0062 double theSearchPhiTOB;
0063 double theSearchPhiTEC;
0064 double theSearchZTIB;
0065 double theSearchZTOB;
0066
0067
0068 TFile *theFile;
0069 int theCompression;
0070 std::string theFileName;
0071
0072
0073 TH1D *theBarrelSimHitsX;
0074 TH1D *theBarrelSimHitsY;
0075 TH1D *theBarrelSimHitsZ;
0076 TH2D *theBarrelSimHitsYvsX;
0077 TH2D *theBarrelSimHitsXvsZ;
0078 TH2D *theBarrelSimHitsYvsZ;
0079 TH2D *theBarrelSimHitsRvsZ;
0080 TH2D *theBarrelSimHitsPhivsX;
0081 TH2D *theBarrelSimHitsPhivsY;
0082 TH2D *theBarrelSimHitsPhivsZ;
0083
0084
0085 TH1D *theEndcapSimHitsX;
0086 TH1D *theEndcapSimHitsY;
0087 TH1D *theEndcapSimHitsZ;
0088 TH2D *theEndcapSimHitsYvsX;
0089 TH2D *theEndcapSimHitsXvsZ;
0090 TH2D *theEndcapSimHitsYvsZ;
0091 TH2D *theEndcapSimHitsRvsZ;
0092 TH2D *theEndcapSimHitsPhivsX;
0093 TH2D *theEndcapSimHitsPhivsY;
0094 TH2D *theEndcapSimHitsPhivsZ;
0095
0096
0097 TH2D *theSimHitsRvsZ;
0098 TH2D *theSimHitsPhivsZ;
0099 };
0100
0101 SimAnalyzer::SimAnalyzer(edm::ParameterSet const &theConf)
0102 : tkGeomToken_(esConsumes()),
0103 theEvents(0),
0104 theDebugLevel(theConf.getUntrackedParameter<int>("DebugLevel", 0)),
0105 theSearchPhiTIB(theConf.getUntrackedParameter<double>("SearchWindowPhiTIB", 0.05)),
0106 theSearchPhiTOB(theConf.getUntrackedParameter<double>("SearchWindowPhiTOB", 0.05)),
0107 theSearchPhiTEC(theConf.getUntrackedParameter<double>("SearchWindowPhiTEC", 0.05)),
0108 theSearchZTIB(theConf.getUntrackedParameter<double>("SearchWindowZTIB", 1.0)),
0109 theSearchZTOB(theConf.getUntrackedParameter<double>("SearchWindowZTOB", 1.0)),
0110 theFile(),
0111 theCompression(theConf.getUntrackedParameter<int>("ROOTFileCompression", 1)),
0112 theFileName(theConf.getUntrackedParameter<std::string>("ROOTFileName", "test.root")),
0113 theBarrelSimHitsX(0),
0114 theBarrelSimHitsY(0),
0115 theBarrelSimHitsZ(0),
0116 theBarrelSimHitsYvsX(0),
0117 theBarrelSimHitsXvsZ(0),
0118 theBarrelSimHitsYvsZ(0),
0119 theBarrelSimHitsRvsZ(0),
0120 theBarrelSimHitsPhivsX(0),
0121 theBarrelSimHitsPhivsY(0),
0122 theBarrelSimHitsPhivsZ(0),
0123
0124 theEndcapSimHitsX(0),
0125 theEndcapSimHitsY(0),
0126 theEndcapSimHitsZ(0),
0127 theEndcapSimHitsYvsX(0),
0128 theEndcapSimHitsXvsZ(0),
0129 theEndcapSimHitsYvsZ(0),
0130 theEndcapSimHitsRvsZ(0),
0131 theEndcapSimHitsPhivsX(0),
0132 theEndcapSimHitsPhivsY(0),
0133 theEndcapSimHitsPhivsZ(0),
0134
0135 theSimHitsRvsZ(0),
0136 theSimHitsPhivsZ(0) {
0137
0138 edm::LogInfo("SimAnalyzer") << "==========================================================="
0139 << "=== Start configuration ==="
0140 << " theDebugLevel = " << theDebugLevel
0141 << " theSearchPhiTIB = " << theSearchPhiTIB
0142 << " theSearchPhiTOB = " << theSearchPhiTOB
0143 << " theSearchPhiTEC = " << theSearchPhiTEC
0144 << " theSearchZTIB = " << theSearchZTIB
0145 << " theSearchZTOB = " << theSearchZTOB
0146 << " ROOT filename = " << theFileName
0147 << " compression = " << theCompression
0148 << "===========================================================";
0149 }
0150
0151 SimAnalyzer::~SimAnalyzer() {
0152 if (theFile != 0) {
0153
0154 closeRootFile();
0155
0156 delete theFile;
0157 }
0158 }
0159
0160 void SimAnalyzer::analyze(edm::Event const &theEvent, edm::EventSetup const &theSetup) {
0161 LogDebug("SimAnalyzer") << "==========================================================="
0162 << " Private analysis of event #" << theEvent.id().event() << " in run #"
0163 << theEvent.id().run();
0164
0165
0166
0167
0168 theEvents++;
0169
0170 LogDebug("SimAnalyzer") << " Total Event number = " << theEvents;
0171
0172
0173 trackerStatistics(theEvent, theSetup);
0174
0175 LogDebug("SimAnalyzer") << "===========================================================";
0176 }
0177
0178 void SimAnalyzer::beginJob() {
0179 LogDebug("SimAnalyzer") << "==========================================================="
0180 << "=== Start beginJob() ==="
0181 << " creating a CMS Root Tree ...";
0182
0183 theFile = new TFile(theFileName.c_str(), "RECREATE", "CMS ROOT file");
0184
0185
0186 if (theFile) {
0187 theFile->SetCompressionLevel(theCompression);
0188 this->initHistograms();
0189 } else {
0190 throw edm::Exception(edm::errors::LogicError,
0191 "<SimAnalyzer::initSetup()>: ERROR!!! something wrong "
0192 "with the RootFile");
0193 }
0194
0195 LogDebug("SimAnalyzer") << "=== Done beginJob() ==="
0196 << "===========================================================";
0197 }
0198
0199 double SimAnalyzer::angle(double theAngle) { return theAngle += (theAngle >= 0.0) ? 0.0 : 2.0 * M_PI; }
0200
0201 void SimAnalyzer::closeRootFile() {
0202 LogDebug("SimAnalyzer") << " writing all information into a file ";
0203
0204 theFile->Write();
0205 }
0206
0207 void SimAnalyzer::initHistograms() {
0208 LogDebug("SimAnalyzer") << " Start of initialisation monitor histograms ...";
0209
0210
0211 TDirectory *SimHitDir = theFile->mkdir("SimHits");
0212 TDirectory *BarrelDir = SimHitDir->mkdir("Barrel");
0213 TDirectory *EndcapDir = SimHitDir->mkdir("Endcap");
0214
0215
0216 theBarrelSimHitsX = new TH1D("BarrelSimHitsX", "X Position of the SimHits", 200, -100.0, 100.0);
0217 theBarrelSimHitsX->SetDirectory(BarrelDir);
0218 theBarrelSimHitsX->Sumw2();
0219
0220 theBarrelSimHitsY = new TH1D("BarrelSimHitsY", "Y Position of the SimHits", 200, -100.0, 100.0);
0221 theBarrelSimHitsY->SetDirectory(BarrelDir);
0222 theBarrelSimHitsY->Sumw2();
0223
0224 theBarrelSimHitsZ = new TH1D("BarrelSimHitsZ", "Z Position of the SimHits", 240, -120.0, 120.0);
0225 theBarrelSimHitsZ->SetDirectory(BarrelDir);
0226 theBarrelSimHitsZ->Sumw2();
0227
0228 theBarrelSimHitsYvsX =
0229 new TH2D("BarrelSimHitsYvsX", "Y vs X Position of the SimHits", 200, -100.0, 100.0, 200, -100.0, 100.0);
0230 theBarrelSimHitsYvsX->SetDirectory(BarrelDir);
0231 theBarrelSimHitsYvsX->Sumw2();
0232
0233 theBarrelSimHitsXvsZ =
0234 new TH2D("BarrelSimHitsXvsZ", "X vs Z Position of the SimHits", 240, -120.0, 120.0, 200, -100.0, 100.0);
0235 theBarrelSimHitsXvsZ->SetDirectory(BarrelDir);
0236 theBarrelSimHitsXvsZ->Sumw2();
0237
0238 theBarrelSimHitsYvsZ =
0239 new TH2D("BarrelSimHitsYvsZ", "Y vs Z Position of the SimHits", 240, -120.0, 120.0, 200, -100.0, 100.0);
0240 theBarrelSimHitsYvsZ->SetDirectory(BarrelDir);
0241 theBarrelSimHitsYvsZ->Sumw2();
0242
0243 theBarrelSimHitsRvsZ =
0244 new TH2D("BarrelSimHitsRvsZ", "R vs Z Position of the SimHits", 240, -120.0, 120.0, 130, 0.0, 65.0);
0245 theBarrelSimHitsRvsZ->SetDirectory(BarrelDir);
0246 theBarrelSimHitsRvsZ->Sumw2();
0247
0248 theBarrelSimHitsPhivsX =
0249 new TH2D("BarrelSimHitsPhivsX", "Phi [rad] vs X Position of the SimHits", 200, -100.0, 100.0, 70, 0.0, 7.0);
0250 theBarrelSimHitsPhivsX->SetDirectory(BarrelDir);
0251 theBarrelSimHitsPhivsX->Sumw2();
0252
0253 theBarrelSimHitsPhivsY =
0254 new TH2D("BarrelSimHitsPhivsY", "Phi [rad] vs Y Position of the SimHits", 200, -100.0, 100.0, 70, 0.0, 7.0);
0255 theBarrelSimHitsPhivsY->SetDirectory(BarrelDir);
0256 theBarrelSimHitsPhivsY->Sumw2();
0257
0258 theBarrelSimHitsPhivsZ =
0259 new TH2D("BarrelSimHitsPhivsZ", "Phi [rad] vs Z Position of the SimHits", 240, -120.0, 120.0, 70, 0.0, 7.0);
0260 theBarrelSimHitsPhivsZ->SetDirectory(BarrelDir);
0261 theBarrelSimHitsPhivsZ->Sumw2();
0262
0263
0264 theEndcapSimHitsX = new TH1D("EndcapSimHitsX", "X Position of the SimHits", 200, -100.0, 100.0);
0265 theEndcapSimHitsX->SetDirectory(EndcapDir);
0266 theEndcapSimHitsX->Sumw2();
0267
0268 theEndcapSimHitsY = new TH1D("EndcapSimHitsY", "Y Position of the SimHits", 200, -100.0, 100.0);
0269 theEndcapSimHitsY->SetDirectory(EndcapDir);
0270 theEndcapSimHitsY->Sumw2();
0271
0272 theEndcapSimHitsZ = new TH1D("EndcapSimHitsZ", "Z Position of the SimHits", 600, -300.0, 300.0);
0273 theEndcapSimHitsZ->SetDirectory(EndcapDir);
0274 theEndcapSimHitsZ->Sumw2();
0275
0276 theEndcapSimHitsYvsX =
0277 new TH2D("EndcapSimHitsYvsX", "Y vs X Position of the SimHits", 200, -100.0, 100.0, 200, -100.0, 100.0);
0278 theEndcapSimHitsYvsX->SetDirectory(EndcapDir);
0279 theEndcapSimHitsYvsX->Sumw2();
0280
0281 theEndcapSimHitsXvsZ =
0282 new TH2D("EndcapSimHitsXvsZ", "X vs Z Position of the SimHits", 600, -300.0, 300.0, 200, -100.0, 100.0);
0283 theEndcapSimHitsXvsZ->SetDirectory(EndcapDir);
0284 theEndcapSimHitsXvsZ->Sumw2();
0285
0286 theEndcapSimHitsYvsZ =
0287 new TH2D("EndcapSimHitsYvsZ", "Y vs Z Position of the SimHits", 600, -300.0, 300.0, 200, -100.0, 100.0);
0288 theEndcapSimHitsYvsZ->SetDirectory(EndcapDir);
0289 theEndcapSimHitsYvsZ->Sumw2();
0290
0291 theEndcapSimHitsRvsZ =
0292 new TH2D("EndcapSimHitsRvsZ", "R vs Z Position of the SimHits", 600, -300.0, 300.0, 1000, 0.0, 100.0);
0293 theEndcapSimHitsRvsZ->SetDirectory(EndcapDir);
0294 theEndcapSimHitsRvsZ->Sumw2();
0295
0296 theEndcapSimHitsPhivsX =
0297 new TH2D("EndcapSimHitsPhivsX", "Phi [rad] vs X Positon of the SimHits", 200, -100.0, 100.0, 70, 0.0, 7.0);
0298 theEndcapSimHitsPhivsX->SetDirectory(EndcapDir);
0299 theEndcapSimHitsPhivsX->Sumw2();
0300
0301 theEndcapSimHitsPhivsY =
0302 new TH2D("EndcapSimHitsPhivsY", "Phi [rad] vs Y Position of the SimHits", 200, -100.0, 100.0, 70, 0.0, 7.0);
0303 theEndcapSimHitsPhivsY->SetDirectory(EndcapDir);
0304 theEndcapSimHitsPhivsY->Sumw2();
0305
0306 theEndcapSimHitsPhivsZ =
0307 new TH2D("EndcapSimHitsPhivsZ", "Phi [rad] vs Z Position of the SimHits", 600, -300.0, 300.0, 70, 0.0, 7.0);
0308 theEndcapSimHitsPhivsZ->SetDirectory(EndcapDir);
0309 theEndcapSimHitsPhivsZ->Sumw2();
0310
0311
0312 theSimHitsRvsZ = new TH2D("SimHitsRvsZ", "R vs Z Position of the SimHits", 600, -300.0, 300.0, 1000, 0.0, 100.0);
0313 theSimHitsRvsZ->SetDirectory(SimHitDir);
0314 theSimHitsRvsZ->Sumw2();
0315
0316 theSimHitsPhivsZ =
0317 new TH2D("SimHitsPhivsZ", "Phi [rad] vs Z Position of the SimHits", 600, -300.0, 300.0, 700, 0.0, 7.0);
0318 theSimHitsPhivsZ->SetDirectory(SimHitDir);
0319 theSimHitsPhivsZ->Sumw2();
0320 }
0321
0322 void SimAnalyzer::trackerStatistics(edm::Event const &theEvent, edm::EventSetup const &theSetup) {
0323 LogDebug("SimAnalyzer") << "<SimAnalyzer:t:rackerStatistics(edm::Event "
0324 "const& theEvent): filling the histograms ... ";
0325
0326 int theBarrelHits = 0;
0327 int theEndcapHits = 0;
0328
0329
0330 const TrackerGeometry *theTracker = &theSetup.getData(tkGeomToken_);
0331
0332
0333 TrackingGeometry::DetContainer theDetUnits = theTracker->dets();
0334
0335
0336 std::vector<edm::Handle<edm::PSimHitContainer>> theSimHitContainers;
0337 theEvent.getManyByType(theSimHitContainers);
0338
0339 LogDebug("SimAnalyzer") << " Geometry node for TrackingGeometry is " << theTracker << "\n I have "
0340 << theTracker->dets().size() << " detectors "
0341 << "\n I have " << theTracker->detTypes().size() << " types "
0342 << "\n theDetUnits has " << theDetUnits.size() << " dets ";
0343
0344
0345 std::vector<PSimHit> theSimHits;
0346 for (int i = 0; i < int(theSimHitContainers.size()); i++) {
0347 theSimHits.insert(theSimHits.end(), theSimHitContainers.at(i)->begin(), theSimHitContainers.at(i)->end());
0348 }
0349
0350
0351 for (std::vector<PSimHit>::const_iterator iHit = theSimHits.begin(); iHit != theSimHits.end(); iHit++) {
0352
0353 DetId theDetUnitId((*iHit).detUnitId());
0354
0355
0356 const GeomDet *theDet = theTracker->idToDet(theDetUnitId);
0357
0358
0359 std::string thePart = "";
0360
0361 LogDebug("SimAnalyzer:SimHitInfo") << " ============================================================ "
0362 << "\n Some Information for this SimHit: "
0363 << "\n"
0364 << *iHit << "\n GlobalPosition = ("
0365 << theDet->surface().toGlobal((*iHit).localPosition()).x() << ","
0366 << theDet->surface().toGlobal((*iHit).localPosition()).y() << ","
0367 << theDet->surface().toGlobal((*iHit).localPosition()).z() << ") "
0368 << "\n R = " << theDet->surface().toGlobal((*iHit).localPosition()).perp()
0369 << "\n phi = " << theDet->surface().toGlobal((*iHit).localPosition()).phi()
0370 << "\n angle(phi) = "
0371 << angle(theDet->surface().toGlobal((*iHit).localPosition()).phi())
0372 << "\n GlobalDirection = ("
0373 << theDet->surface().toGlobal((*iHit).localDirection()).x() << ","
0374 << theDet->surface().toGlobal((*iHit).localDirection()).y() << ","
0375 << theDet->surface().toGlobal((*iHit).localDirection()).z() << ") "
0376 << "\n Total momentum = " << (*iHit).pabs()
0377 << "\n ParticleType = " << (*iHit).particleType()
0378 << "\n detUnitId = " << (*iHit).detUnitId();
0379
0380
0381 switch (theDetUnitId.subdetId()) {
0382 case StripSubdetector::TIB: {
0383 thePart = "TIB";
0384 break;
0385 }
0386 case StripSubdetector::TOB: {
0387 thePart = "TOB";
0388 break;
0389 }
0390 case StripSubdetector::TEC: {
0391 thePart = "TEC";
0392 break;
0393 }
0394 }
0395
0396 if ((thePart == "TIB") || (thePart == "TOB")) {
0397 theBarrelHits++;
0398
0399
0400 theBarrelSimHitsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x());
0401 theBarrelSimHitsY->Fill(theDet->surface().toGlobal((*iHit).localPosition()).y());
0402 theBarrelSimHitsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z());
0403 theBarrelSimHitsYvsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x(),
0404 theDet->surface().toGlobal((*iHit).localPosition()).y());
0405 theBarrelSimHitsXvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0406 theDet->surface().toGlobal((*iHit).localPosition()).x());
0407 theBarrelSimHitsYvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0408 theDet->surface().toGlobal((*iHit).localPosition()).y());
0409 theBarrelSimHitsRvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0410 theDet->surface().toGlobal((*iHit).localPosition()).perp());
0411 theBarrelSimHitsPhivsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x(),
0412 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0413 theBarrelSimHitsPhivsY->Fill(theDet->surface().toGlobal((*iHit).localPosition()).y(),
0414 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0415 theBarrelSimHitsPhivsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0416 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0417 } else if (thePart == "TEC") {
0418 theEndcapHits++;
0419
0420
0421 theEndcapSimHitsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x());
0422 theEndcapSimHitsY->Fill(theDet->surface().toGlobal((*iHit).localPosition()).y());
0423 theEndcapSimHitsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z());
0424 theEndcapSimHitsYvsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x(),
0425 theDet->surface().toGlobal((*iHit).localPosition()).y());
0426 theEndcapSimHitsXvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0427 theDet->surface().toGlobal((*iHit).localPosition()).x());
0428 theEndcapSimHitsYvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0429 theDet->surface().toGlobal((*iHit).localPosition()).y());
0430 theEndcapSimHitsRvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0431 theDet->surface().toGlobal((*iHit).localPosition()).perp());
0432 theEndcapSimHitsPhivsX->Fill(theDet->surface().toGlobal((*iHit).localPosition()).x(),
0433 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0434 theEndcapSimHitsPhivsY->Fill(theDet->surface().toGlobal((*iHit).localPosition()).y(),
0435 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0436 theEndcapSimHitsPhivsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0437 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0438 }
0439
0440
0441 if ((thePart == "TIB") || (thePart == "TOB") || (thePart == "TEC")) {
0442 theSimHitsRvsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0443 theDet->surface().toGlobal((*iHit).localPosition()).perp());
0444 theSimHitsPhivsZ->Fill(theDet->surface().toGlobal((*iHit).localPosition()).z(),
0445 angle(theDet->surface().toGlobal((*iHit).localPosition()).phi()));
0446 }
0447
0448 }
0449
0450
0451 edm::LogInfo("SimAnalyzer") << " number of SimHits = " << theBarrelHits << "/" << theEndcapHits
0452 << " (Barrel/Endcap) ";
0453 }
0454
0455 #include "FWCore/Framework/interface/MakerMacros.h"
0456 DEFINE_FWK_MODULE(SimAnalyzer);