File indexing completed on 2024-04-06 12:11:18
0001
0002 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008
0009 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0014 #include "FastSimulation/Event/interface/FSimEvent.h"
0015 #include "FastSimulation/Event/interface/FSimTrack.h"
0016 #include "FastSimulation/Event/interface/FSimVertex.h"
0017 #include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"
0018
0019 #include <memory>
0020 #include <string>
0021 #include <vector>
0022 #include "TFile.h"
0023 #include "TTree.h"
0024 #include "TProcessID.h"
0025
0026 class testNuclearInteractions : public DQMEDAnalyzer {
0027 public:
0028 explicit testNuclearInteractions(const edm::ParameterSet&);
0029 ~testNuclearInteractions();
0030
0031 virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0032 virtual void dqmBeginRun(edm::Run const&, const edm::EventSetup&) override;
0033 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0034
0035 private:
0036
0037 const edm::ParameterSet particleFilter_;
0038 const bool saveNU;
0039 const int maxNU;
0040 const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> tok_pdt_;
0041 const edm::EDGetTokenT<std::vector<SimTrack>> tok_fullSimTk_;
0042 const edm::EDGetTokenT<std::vector<SimVertex>> tok_fullSimVx_;
0043 const edm::EDGetTokenT<std::vector<SimTrack>> tok_fastSimTk_;
0044 const edm::EDGetTokenT<std::vector<SimVertex>> tok_fastSimVx_;
0045 std::vector<std::unique_ptr<FSimEvent>> mySimEvent;
0046 NUEvent* nuEvent;
0047 TTree* nuTree;
0048 TFile* outFile;
0049 int ObjectNumber;
0050 std::string simModuleLabel_;
0051
0052
0053 std::vector<MonitorElement*> h0;
0054 std::vector<MonitorElement*> h1;
0055 std::vector<MonitorElement*> h2;
0056 std::vector<MonitorElement*> h3;
0057 std::vector<MonitorElement*> h4;
0058 std::vector<MonitorElement*> h5;
0059 std::vector<MonitorElement*> h6;
0060 std::vector<MonitorElement*> h7;
0061 std::vector<MonitorElement*> h8;
0062 std::vector<MonitorElement*> h9;
0063 std::vector<MonitorElement*> h10;
0064 std::vector<MonitorElement*> htmp;
0065 std::vector<MonitorElement*> totalCharge;
0066
0067 std::vector<std::vector<MonitorElement*>> h100;
0068 std::vector<std::vector<MonitorElement*>> h200;
0069 std::vector<std::vector<MonitorElement*>> h300;
0070 std::vector<std::vector<double>> trackerRadius;
0071 std::vector<std::vector<double>> trackerLength;
0072 std::vector<std::vector<double>> blockTrackerRadius;
0073 std::vector<std::vector<double>> blockTrackerLength;
0074 std::vector<std::vector<double>> subTrackerRadius;
0075 std::vector<std::vector<double>> subTrackerLength;
0076 std::vector<double> tmpRadius;
0077 std::vector<double> tmpLength;
0078
0079 std::vector<unsigned> stoppedPions;
0080 std::vector<unsigned> interactingPions;
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 int intfull;
0091 int intfast;
0092
0093 std::string NUEventFileName;
0094
0095 int totalNEvt;
0096 int totalNU;
0097 };
0098
0099 testNuclearInteractions::testNuclearInteractions(const edm::ParameterSet& p)
0100 : particleFilter_(p.getParameter<edm::ParameterSet>("TestParticleFilter")),
0101 saveNU(p.getParameter<bool>("SaveNuclearInteractions")),
0102 maxNU(p.getParameter<unsigned>("MaxNumberOfNuclearInteractions")),
0103 tok_pdt_(esConsumes<HepPDT::ParticleDataTable, PDTRecord>()),
0104 tok_fullSimTk_(consumes<std::vector<SimTrack>>(edm::InputTag("g4SimHits"))),
0105 tok_fullSimVx_(consumes<std::vector<SimVertex>>(edm::InputTag("g4SimHits"))),
0106 tok_fastSimTk_(consumes<std::vector<SimTrack>>(edm::InputTag("famosSimHits"))),
0107 tok_fastSimVx_(consumes<std::vector<SimVertex>>(edm::InputTag("famosSimHits"))),
0108 h0(2, static_cast<MonitorElement*>(0)),
0109 h1(2, static_cast<MonitorElement*>(0)),
0110 h2(2, static_cast<MonitorElement*>(0)),
0111 h3(2, static_cast<MonitorElement*>(0)),
0112 h4(2, static_cast<MonitorElement*>(0)),
0113 h5(2, static_cast<MonitorElement*>(0)),
0114 h6(2, static_cast<MonitorElement*>(0)),
0115 h7(2, static_cast<MonitorElement*>(0)),
0116 h8(2, static_cast<MonitorElement*>(0)),
0117 h9(2, static_cast<MonitorElement*>(0)),
0118 h10(2, static_cast<MonitorElement*>(0)),
0119 htmp(2, static_cast<MonitorElement*>(0)),
0120 totalCharge(2, static_cast<MonitorElement*>(0)),
0121 tmpRadius(2, static_cast<double>(0.)),
0122 tmpLength(2, static_cast<double>(0.)),
0123 stoppedPions(2, static_cast<unsigned>(0)),
0124 interactingPions(2, static_cast<unsigned>(0)),
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134 intfull(0),
0135 intfast(0),
0136 totalNEvt(0),
0137 totalNU(0) {
0138
0139
0140 mySimEvent.emplace_back(std::make_unique<FSimEvent>(particleFilter_));
0141
0142 mySimEvent.emplace_back(std::make_unique<FSimEvent>(particleFilter_));
0143
0144
0145 if (saveNU)
0146 std::cout << "Nuclear Interactions will be saved ! " << std::endl;
0147
0148
0149 NUEventFileName = "none";
0150 if (saveNU) {
0151 nuEvent = new NUEvent();
0152
0153 NUEventFileName = p.getUntrackedParameter<std::string>("NUEventFile", "NuclearInteractionsTest.root");
0154
0155 outFile = new TFile(NUEventFileName.c_str(), "recreate");
0156
0157
0158 nuTree = new TTree("NuclearInteractions", "");
0159 nuTree->Branch("nuEvent", "NUEvent", &nuEvent, 32000, 99);
0160 }
0161
0162
0163 ObjectNumber = -1;
0164
0165
0166 std::vector<double> tmpRadius = p.getUntrackedParameter<std::vector<double>>("BPCylinderRadius");
0167 std::vector<double> tmpLength = p.getUntrackedParameter<std::vector<double>>("BPCylinderLength");
0168 trackerRadius.push_back(tmpRadius);
0169 trackerLength.push_back(tmpLength);
0170
0171
0172 subTrackerRadius.push_back(tmpRadius);
0173 subTrackerLength.push_back(tmpLength);
0174
0175
0176 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB1CylinderRadius");
0177 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB1CylinderLength");
0178 trackerRadius.push_back(tmpRadius);
0179 trackerLength.push_back(tmpLength);
0180
0181
0182 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB2CylinderRadius");
0183 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB2CylinderLength");
0184 trackerRadius.push_back(tmpRadius);
0185 trackerLength.push_back(tmpLength);
0186
0187
0188 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB3CylinderRadius");
0189 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB3CylinderLength");
0190 trackerRadius.push_back(tmpRadius);
0191 trackerLength.push_back(tmpLength);
0192
0193
0194 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXBCablesCylinderRadius");
0195 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXBCablesCylinderLength");
0196 trackerRadius.push_back(tmpRadius);
0197 trackerLength.push_back(tmpLength);
0198
0199
0200 blockTrackerRadius.push_back(tmpRadius);
0201 blockTrackerLength.push_back(tmpLength);
0202
0203
0204 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXD1CylinderRadius");
0205 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXD1CylinderLength");
0206 trackerRadius.push_back(tmpRadius);
0207 trackerLength.push_back(tmpLength);
0208
0209
0210 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXD2CylinderRadius");
0211 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXD2CylinderLength");
0212 trackerRadius.push_back(tmpRadius);
0213 trackerLength.push_back(tmpLength);
0214
0215
0216 tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXDCablesCylinderRadius");
0217 tmpLength = p.getUntrackedParameter<std::vector<double>>("PXDCablesCylinderLength");
0218 trackerRadius.push_back(tmpRadius);
0219 trackerLength.push_back(tmpLength);
0220
0221
0222 blockTrackerRadius.push_back(tmpRadius);
0223 blockTrackerLength.push_back(tmpLength);
0224
0225
0226 subTrackerRadius.push_back(tmpRadius);
0227 subTrackerLength.push_back(tmpLength);
0228
0229
0230 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB1CylinderRadius");
0231 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB1CylinderLength");
0232 trackerRadius.push_back(tmpRadius);
0233 trackerLength.push_back(tmpLength);
0234
0235
0236 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB2CylinderRadius");
0237 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB2CylinderLength");
0238 trackerRadius.push_back(tmpRadius);
0239 trackerLength.push_back(tmpLength);
0240
0241
0242 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB3CylinderRadius");
0243 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB3CylinderLength");
0244 trackerRadius.push_back(tmpRadius);
0245 trackerLength.push_back(tmpLength);
0246
0247
0248 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB4CylinderRadius");
0249 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB4CylinderLength");
0250 trackerRadius.push_back(tmpRadius);
0251 trackerLength.push_back(tmpLength);
0252
0253
0254 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIBCablesCylinderRadius");
0255 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIBCablesCylinderLength");
0256 trackerRadius.push_back(tmpRadius);
0257 trackerLength.push_back(tmpLength);
0258
0259
0260 blockTrackerRadius.push_back(tmpRadius);
0261 blockTrackerLength.push_back(tmpLength);
0262
0263
0264 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID1CylinderRadius");
0265 tmpLength = p.getUntrackedParameter<std::vector<double>>("TID1CylinderLength");
0266 trackerRadius.push_back(tmpRadius);
0267 trackerLength.push_back(tmpLength);
0268
0269
0270 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID2CylinderRadius");
0271 tmpLength = p.getUntrackedParameter<std::vector<double>>("TID2CylinderLength");
0272 trackerRadius.push_back(tmpRadius);
0273 trackerLength.push_back(tmpLength);
0274
0275
0276 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID3CylinderRadius");
0277 tmpLength = p.getUntrackedParameter<std::vector<double>>("TID3CylinderLength");
0278 trackerRadius.push_back(tmpRadius);
0279 trackerLength.push_back(tmpLength);
0280
0281
0282 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIDCablesCylinderRadius");
0283 tmpLength = p.getUntrackedParameter<std::vector<double>>("TIDCablesCylinderLength");
0284 trackerRadius.push_back(tmpRadius);
0285 trackerLength.push_back(tmpLength);
0286
0287
0288 blockTrackerRadius.push_back(tmpRadius);
0289 blockTrackerLength.push_back(tmpLength);
0290
0291
0292 subTrackerRadius.push_back(tmpRadius);
0293 subTrackerLength.push_back(tmpLength);
0294
0295
0296 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB1CylinderRadius");
0297 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB1CylinderLength");
0298 trackerRadius.push_back(tmpRadius);
0299 trackerLength.push_back(tmpLength);
0300
0301
0302 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB2CylinderRadius");
0303 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB2CylinderLength");
0304 trackerRadius.push_back(tmpRadius);
0305 trackerLength.push_back(tmpLength);
0306
0307
0308 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB3CylinderRadius");
0309 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB3CylinderLength");
0310 trackerRadius.push_back(tmpRadius);
0311 trackerLength.push_back(tmpLength);
0312
0313
0314 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB4CylinderRadius");
0315 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB4CylinderLength");
0316 trackerRadius.push_back(tmpRadius);
0317 trackerLength.push_back(tmpLength);
0318
0319
0320 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB5CylinderRadius");
0321 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB5CylinderLength");
0322 trackerRadius.push_back(tmpRadius);
0323 trackerLength.push_back(tmpLength);
0324
0325
0326 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB6CylinderRadius");
0327 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB6CylinderLength");
0328 trackerRadius.push_back(tmpRadius);
0329 trackerLength.push_back(tmpLength);
0330
0331
0332 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOBCablesCylinderRadius");
0333 tmpLength = p.getUntrackedParameter<std::vector<double>>("TOBCablesCylinderLength");
0334 trackerRadius.push_back(tmpRadius);
0335 trackerLength.push_back(tmpLength);
0336
0337
0338 blockTrackerRadius.push_back(tmpRadius);
0339 blockTrackerLength.push_back(tmpLength);
0340
0341
0342 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC1CylinderRadius");
0343 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC1CylinderLength");
0344 trackerRadius.push_back(tmpRadius);
0345 trackerLength.push_back(tmpLength);
0346
0347
0348 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC2CylinderRadius");
0349 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC2CylinderLength");
0350 trackerRadius.push_back(tmpRadius);
0351 trackerLength.push_back(tmpLength);
0352
0353
0354 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC3CylinderRadius");
0355 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC3CylinderLength");
0356 trackerRadius.push_back(tmpRadius);
0357 trackerLength.push_back(tmpLength);
0358
0359
0360 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC4CylinderRadius");
0361 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC4CylinderLength");
0362 trackerRadius.push_back(tmpRadius);
0363 trackerLength.push_back(tmpLength);
0364
0365
0366 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC5CylinderRadius");
0367 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC5CylinderLength");
0368 trackerRadius.push_back(tmpRadius);
0369 trackerLength.push_back(tmpLength);
0370
0371
0372 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC6CylinderRadius");
0373 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC6CylinderLength");
0374 trackerRadius.push_back(tmpRadius);
0375 trackerLength.push_back(tmpLength);
0376
0377
0378 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC7CylinderRadius");
0379 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC7CylinderLength");
0380 trackerRadius.push_back(tmpRadius);
0381 trackerLength.push_back(tmpLength);
0382
0383
0384 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC8CylinderRadius");
0385 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC8CylinderLength");
0386 trackerRadius.push_back(tmpRadius);
0387 trackerLength.push_back(tmpLength);
0388
0389
0390 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC9CylinderRadius");
0391 tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC9CylinderLength");
0392 trackerRadius.push_back(tmpRadius);
0393 trackerLength.push_back(tmpLength);
0394
0395
0396 blockTrackerRadius.push_back(tmpRadius);
0397 blockTrackerLength.push_back(tmpLength);
0398
0399
0400 subTrackerRadius.push_back(tmpRadius);
0401 subTrackerLength.push_back(tmpLength);
0402
0403
0404 tmpRadius = p.getUntrackedParameter<std::vector<double>>("TrackerCablesCylinderRadius");
0405 tmpLength = p.getUntrackedParameter<std::vector<double>>("TrackerCablesCylinderLength");
0406 trackerRadius.push_back(tmpRadius);
0407 trackerLength.push_back(tmpLength);
0408
0409
0410 h200.push_back(htmp);
0411 blockTrackerRadius.push_back(tmpRadius);
0412 blockTrackerLength.push_back(tmpLength);
0413
0414
0415 h300.push_back(htmp);
0416 subTrackerRadius.push_back(tmpRadius);
0417 subTrackerLength.push_back(tmpLength);
0418 }
0419
0420 void testNuclearInteractions::bookHistograms(DQMStore::IBooker& ibooker,
0421 edm::Run const& iRun,
0422 edm::EventSetup const& iSetup) {
0423 ibooker.setCurrentFolder("testNuclearInteractions");
0424
0425 h0[0] = ibooker.book2D("radioFull", "Full Tracker radiography", 1000, 0., 320., 1000, 0., 150.);
0426 h0[1] = ibooker.book2D("radioFast", "Fast Tracker radiography", 1000, 0., 320., 1000, 0., 150.);
0427 h1[0] = ibooker.book1D("vertexFull", "Full Nb of Vertices", 20, -0.5, 19.5);
0428 h1[1] = ibooker.book1D("vertexFast", "Fast Nb of Vertices", 20, -0.5, 19.5);
0429 h2[0] = ibooker.book1D("daughterFull", "Full Nb of daughters", 20, -0.5, 19.5);
0430 h2[1] = ibooker.book1D("daughterFast", "Fast Nb of daughters", 20, -0.5, 19.5);
0431 h3[0] = ibooker.book1D("ecmFull", "Full centre-of-mass energy", 100, 0., 10.);
0432 h3[1] = ibooker.book1D("ecmFast", "Fast centre-of-mass energy", 100, 0., 10.);
0433 h4[0] = ibooker.book1D("FecmFull", "Full c.m. energy fraction", 100, 0., 2.);
0434 h4[1] = ibooker.book1D("FecmFast", "Fast c.m. energy fraction", 100, 0., 2.);
0435 h5[0] = ibooker.book1D("FmomFull", "Full momemtum", 100, 0., 10.);
0436 h5[1] = ibooker.book1D("FmomFast", "Fast momemtum", 100, 0., 10.);
0437 h6[0] = ibooker.book1D("DeltaEFull4", "Full DeltaE", 2000, -1., 4.);
0438 h6[1] = ibooker.book1D("DeltaEFast4", "Fast DetlaE", 2000, -1., 4.);
0439 h7[0] = ibooker.book1D("DeltaEFull3", "Full DeltaE 3 daugh", 2000, -1., 4.);
0440 h7[1] = ibooker.book1D("DeltaEFast3", "Fast DetlaE 3 daugh", 2000, -1., 4.);
0441 h8[0] = ibooker.book1D("DeltaMFull4", "Full DeltaE", 2000, -10., 40.);
0442 h8[1] = ibooker.book1D("DeltaMFast4", "Fast DetlaE", 2000, -10., 40.);
0443 h9[0] = ibooker.book1D("DeltaMFull3", "Full DeltaE 3 daugh", 2000, -10., 40.);
0444 h9[1] = ibooker.book1D("DeltaMFast3", "Fast DetlaE 3 daugh", 2000, -10., 40.);
0445 h10[0] = ibooker.book1D("EafterFull", "E(after)/E(before) full", 200, 0., 4.);
0446 h10[1] = ibooker.book1D("EafterFast", "E(after)/E(before) fast", 200, 0., 4.);
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464 totalCharge[0] = ibooker.book1D("ChargeFull", "Total Charge (full)", 19, -9.5, 9.5);
0465 totalCharge[1] = ibooker.book1D("ChargeFast", "Total Charge (fast)", 19, -9.5, 9.5);
0466
0467
0468 htmp[0] = ibooker.book1D("BeamPipeFull", "Full Beam Pipe", 120, 0., 3.);
0469 htmp[1] = ibooker.book1D("BeamPipeFast", "Fast Beam Pipe", 120, 0., 3.);
0470 h100.push_back(htmp);
0471
0472
0473 htmp[0] = ibooker.book1D("BPFull", "Full Beam Pipe", 120, 0., 3.);
0474 htmp[1] = ibooker.book1D("BPFast", "Fast Beam Pipe", 120, 0., 3.);
0475 h300.push_back(htmp);
0476
0477
0478 htmp[0] = ibooker.book1D("PXB1Full", "Full Pixel Barrel 1", 120, 0., 3.);
0479 htmp[1] = ibooker.book1D("PXB1Fast", "Fast Pixel Barrel 1", 120, 0., 3.);
0480 h100.push_back(htmp);
0481
0482
0483 htmp[0] = ibooker.book1D("PXB2Full", "Full Pixel Barrel 2", 120, 0., 3.);
0484 htmp[1] = ibooker.book1D("PXB2Fast", "Fast Pixel Barrel 2", 120, 0., 3.);
0485 h100.push_back(htmp);
0486
0487
0488 htmp[0] = ibooker.book1D("PXB3Full", "Full Pixel Barrel 3", 120, 0., 3.);
0489 htmp[1] = ibooker.book1D("PXB3Fast", "Fast Pixel Barrel 3", 120, 0., 3.);
0490 h100.push_back(htmp);
0491
0492
0493 htmp[0] = ibooker.book1D("PXBCFull", "Full Pixel Barrel Cables", 120, 0., 3.);
0494 htmp[1] = ibooker.book1D("PXBCFast", "Fast Pixel Barrel Cables", 120, 0., 3.);
0495 h100.push_back(htmp);
0496
0497
0498 htmp[0] = ibooker.book1D("PXBFull", "Full Pixel Barrel", 120, 0., 3.);
0499 htmp[1] = ibooker.book1D("PXBFast", "Fast Pixel Barrel", 120, 0., 3.);
0500 h200.push_back(htmp);
0501
0502
0503 htmp[0] = ibooker.book1D("PXD1Full", "Full Pixel Disk 1", 120, 0., 3.);
0504 htmp[1] = ibooker.book1D("PXD1Fast", "Fast Pixel Disk 1", 120, 0., 3.);
0505 h100.push_back(htmp);
0506
0507
0508 htmp[0] = ibooker.book1D("PXD2Full", "Full Pixel Disk 2", 120, 0., 3.);
0509 htmp[1] = ibooker.book1D("PXD2Fast", "Fast Pixel Disk 2", 120, 0., 3.);
0510 h100.push_back(htmp);
0511
0512
0513 htmp[0] = ibooker.book1D("PXDCFull", "Full Pixel Disk Cables", 120, 0., 3.);
0514 htmp[1] = ibooker.book1D("PXDCFast", "Fast Pixel Disk Cables", 120, 0., 3.);
0515 h100.push_back(htmp);
0516
0517
0518 htmp[0] = ibooker.book1D("PXDFull", "Full Pixel Disk", 120, 0., 3.);
0519 htmp[1] = ibooker.book1D("PXDFast", "Fast Pixel Disk", 120, 0., 3.);
0520 h200.push_back(htmp);
0521
0522
0523 htmp[0] = ibooker.book1D("PixelFull", "Full Pixel", 120, 0., 3.);
0524 htmp[1] = ibooker.book1D("PixelFast", "Fast Pixel", 120, 0., 3.);
0525 h300.push_back(htmp);
0526
0527
0528 htmp[0] = ibooker.book1D("TIB1Full", "Full Tracker Inner Barrel 1", 120, 0., 3.);
0529 htmp[1] = ibooker.book1D("TIB1Fast", "Fast Tracker Inner Barrel 1", 120, 0., 3.);
0530 h100.push_back(htmp);
0531
0532
0533 htmp[0] = ibooker.book1D("TIB2Full", "Full Tracker Inner Barrel 2", 120, 0., 3.);
0534 htmp[1] = ibooker.book1D("TIB2Fast", "Fast Tracker Inner Barrel 2", 120, 0., 3.);
0535 h100.push_back(htmp);
0536
0537
0538 htmp[0] = ibooker.book1D("TIB3Full", "Full Tracker Inner Barrel 3", 120, 0., 3.);
0539 htmp[1] = ibooker.book1D("TIB3Fast", "Fast Tracker Inner Barrel 3", 120, 0., 3.);
0540 h100.push_back(htmp);
0541
0542
0543 htmp[0] = ibooker.book1D("TIB4Full", "Full Tracker Inner Barrel 4", 120, 0., 3.);
0544 htmp[1] = ibooker.book1D("TIB4Fast", "Fast Tracker Inner Barrel 4", 120, 0., 3.);
0545 h100.push_back(htmp);
0546
0547
0548 htmp[0] = ibooker.book1D("TIBCFull", "Full Tracker Inner Barrel Cables", 120, 0., 3.);
0549 htmp[1] = ibooker.book1D("TIBCFast", "Fast Tracker Inner Barrel Cables", 120, 0., 3.);
0550 h100.push_back(htmp);
0551
0552
0553 htmp[0] = ibooker.book1D("TIBFull", "Full Tracker Inner Barrel", 120, 0., 3.);
0554 htmp[1] = ibooker.book1D("TIBFast", "Fast Tracker Inner Barrel", 120, 0., 3.);
0555 h200.push_back(htmp);
0556
0557
0558 htmp[0] = ibooker.book1D("TID1Full", "Full Tracker Inner Disk 1", 120, 0., 3.);
0559 htmp[1] = ibooker.book1D("TID1Fast", "Fast Tracker Inner Disk 1", 120, 0., 3.);
0560 h100.push_back(htmp);
0561
0562
0563 htmp[0] = ibooker.book1D("TID2Full", "Full Tracker Inner Disk 2", 120, 0., 3.);
0564 htmp[1] = ibooker.book1D("TID2Fast", "Fast Tracker Inner Disk 2", 120, 0., 3.);
0565 h100.push_back(htmp);
0566
0567
0568 htmp[0] = ibooker.book1D("TID3Full", "Full Tracker Inner Disk 3", 120, 0., 3.);
0569 htmp[1] = ibooker.book1D("TID3Fast", "Fast Tracker Inner Disk 3", 120, 0., 3.);
0570 h100.push_back(htmp);
0571
0572
0573 htmp[0] = ibooker.book1D("TIDCFull", "Full Tracker Inner Disk Cables", 120, 0., 3.);
0574 htmp[1] = ibooker.book1D("TIDCFast", "Fast Tracker Inner Disk Cables", 120, 0., 3.);
0575 h100.push_back(htmp);
0576
0577
0578 htmp[0] = ibooker.book1D("TIDFull", "Full Tracker Inner Disk", 120, 0., 3.);
0579 htmp[1] = ibooker.book1D("TIDFast", "Fast Tracker Inner Disk", 120, 0., 3.);
0580 h200.push_back(htmp);
0581
0582
0583 htmp[0] = ibooker.book1D("InnerFull", "Full Inner Tracker", 120, 0., 3.);
0584 htmp[1] = ibooker.book1D("InnerFast", "Fast Inner Tracker", 120, 0., 3.);
0585 h300.push_back(htmp);
0586
0587
0588 htmp[0] = ibooker.book1D("TOB1Full", "Full Tracker Outer Barrel 1", 120, 0., 3.);
0589 htmp[1] = ibooker.book1D("TOB1Fast", "Fast Tracker Outer Barrel 1", 120, 0., 3.);
0590 h100.push_back(htmp);
0591
0592
0593 htmp[0] = ibooker.book1D("TOB2Full", "Full Tracker Outer Barrel 2", 120, 0., 3.);
0594 htmp[1] = ibooker.book1D("TOB2Fast", "Fast Tracker Outer Barrel 2", 120, 0., 3.);
0595 h100.push_back(htmp);
0596
0597
0598 htmp[0] = ibooker.book1D("TOB3Full", "Full Tracker Outer Barrel 3", 120, 0., 3.);
0599 htmp[1] = ibooker.book1D("TOB3Fast", "Fast Tracker Outer Barrel 3", 120, 0., 3.);
0600 h100.push_back(htmp);
0601
0602
0603 htmp[0] = ibooker.book1D("TOB4Full", "Full Tracker Outer Barrel 4", 120, 0., 3.);
0604 htmp[1] = ibooker.book1D("TOB4Fast", "Fast Tracker Outer Barrel 4", 120, 0., 3.);
0605 h100.push_back(htmp);
0606
0607
0608 htmp[0] = ibooker.book1D("TOB5Full", "Full Tracker Outer Barrel 5", 120, 0., 3.);
0609 htmp[1] = ibooker.book1D("TOB5Fast", "Fast Tracker Outer Barrel 5", 120, 0., 3.);
0610 h100.push_back(htmp);
0611
0612
0613 htmp[0] = ibooker.book1D("TOB6Full", "Full Tracker Outer Barrel 6", 120, 0., 3.);
0614 htmp[1] = ibooker.book1D("TOB6Fast", "Fast Tracker Outer Barrel 6", 120, 0., 3.);
0615 h100.push_back(htmp);
0616
0617
0618 htmp[0] = ibooker.book1D("TOBCFull", "Full Tracker Outer Barrel Cables", 120, 0., 3.);
0619 htmp[1] = ibooker.book1D("TOBCFast", "Fast Tracker Outer Barrel Cables", 120, 0., 3.);
0620 h100.push_back(htmp);
0621
0622
0623 htmp[0] = ibooker.book1D("TOBFull", "Full Tracker Outer Barrel", 120, 0., 3.);
0624 htmp[1] = ibooker.book1D("TOBFast", "Fast Tracker Outer Barrel", 120, 0., 3.);
0625 h200.push_back(htmp);
0626
0627
0628 htmp[0] = ibooker.book1D("TEC1Full", "Full Tracker EndCap 1", 120, 0., 3.);
0629 htmp[1] = ibooker.book1D("TEC1Fast", "Fast Tracker Endcap 1", 120, 0., 3.);
0630 h100.push_back(htmp);
0631
0632
0633 htmp[0] = ibooker.book1D("TEC2Full", "Full Tracker EndCap 2", 120, 0., 3.);
0634 htmp[1] = ibooker.book1D("TEC2Fast", "Fast Tracker Endcap 2", 120, 0., 3.);
0635 h100.push_back(htmp);
0636
0637
0638 htmp[0] = ibooker.book1D("TEC3Full", "Full Tracker EndCap 3", 120, 0., 3.);
0639 htmp[1] = ibooker.book1D("TEC3Fast", "Fast Tracker Endcap 3", 120, 0., 3.);
0640 h100.push_back(htmp);
0641
0642
0643 htmp[0] = ibooker.book1D("TEC4Full", "Full Tracker EndCap 4", 120, 0., 3.);
0644 htmp[1] = ibooker.book1D("TEC4Fast", "Fast Tracker Endcap 4", 120, 0., 3.);
0645 h100.push_back(htmp);
0646
0647
0648 htmp[0] = ibooker.book1D("TEC5Full", "Full Tracker EndCap 5", 120, 0., 3.);
0649 htmp[1] = ibooker.book1D("TEC5Fast", "Fast Tracker Endcap 5", 120, 0., 3.);
0650 h100.push_back(htmp);
0651
0652
0653 htmp[0] = ibooker.book1D("TEC6Full", "Full Tracker EndCap 6", 120, 0., 3.);
0654 htmp[1] = ibooker.book1D("TEC6Fast", "Fast Tracker Endcap 6", 120, 0., 3.);
0655 h100.push_back(htmp);
0656
0657
0658 htmp[0] = ibooker.book1D("TEC7Full", "Full Tracker EndCap 7", 120, 0., 3.);
0659 htmp[1] = ibooker.book1D("TEC7Fast", "Fast Tracker Endcap 7", 120, 0., 3.);
0660 h100.push_back(htmp);
0661
0662
0663 htmp[0] = ibooker.book1D("TEC8Full", "Full Tracker EndCap 8", 120, 0., 3.);
0664 htmp[1] = ibooker.book1D("TEC8Fast", "Fast Tracker Endcap 8", 120, 0., 3.);
0665 h100.push_back(htmp);
0666
0667
0668 htmp[0] = ibooker.book1D("TEC9Full", "Full Tracker EndCap 9", 120, 0., 3.);
0669 htmp[1] = ibooker.book1D("TEC9Fast", "Fast Tracker Endcap 9", 120, 0., 3.);
0670 h100.push_back(htmp);
0671
0672
0673 htmp[0] = ibooker.book1D("TECFull", "Full Tracker EndCap", 120, 0., 3.);
0674 htmp[1] = ibooker.book1D("TECFast", "Fast Tracker EndCap", 120, 0., 3.);
0675 h200.push_back(htmp);
0676
0677
0678 htmp[0] = ibooker.book1D("OuterFull", "Full Outer Tracker", 120, 0., 3.);
0679 htmp[1] = ibooker.book1D("OuterFast", "Fast Outer Tracker", 120, 0., 3.);
0680 h300.push_back(htmp);
0681
0682
0683 htmp[0] = ibooker.book1D("TECCFull", "Full Tracker Outer Cables", 120, 0., 3.);
0684 htmp[1] = ibooker.book1D("TECCFast", "Fast Tracker Outer Cables", 120, 0., 3.);
0685 h100.push_back(htmp);
0686
0687
0688 htmp[0] = ibooker.book1D("CablesFull", "Full Tracker Cables", 120, 0., 3.);
0689 htmp[1] = ibooker.book1D("CablesFast", "Fast Tracker Cables", 120, 0., 3.);
0690 h200.push_back(htmp);
0691
0692
0693 htmp[0] = ibooker.book1D("TrackerFull", "Full Tracker", 120, 0., 3.);
0694 htmp[1] = ibooker.book1D("TrackerFast", "Fast Tracker", 120, 0., 3.);
0695 h300.push_back(htmp);
0696 }
0697
0698 testNuclearInteractions::~testNuclearInteractions() {
0699 if (saveNU) {
0700 outFile->cd();
0701
0702 nuTree->Fill();
0703
0704 nuTree->Write();
0705
0706 nuTree->Print();
0707
0708
0709 delete nuEvent;
0710 delete nuTree;
0711 delete outFile;
0712 }
0713
0714
0715 }
0716
0717 void testNuclearInteractions::dqmBeginRun(edm::Run const&, const edm::EventSetup& es) {
0718
0719 const edm::ESHandle<HepPDT::ParticleDataTable> pdt = es.getHandle(tok_pdt_);
0720
0721 mySimEvent[0]->initializePdt(&(*pdt));
0722 mySimEvent[1]->initializePdt(&(*pdt));
0723 }
0724
0725 void testNuclearInteractions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0726 ++totalNEvt;
0727 if (totalNEvt / 1000 * 1000 == totalNEvt)
0728 std::cout << "Number of event analysed/NU " << totalNEvt << " / " << totalNU << std::endl;
0729
0730 std::unique_ptr<edm::SimTrackContainer> nuclSimTracks(new edm::SimTrackContainer);
0731
0732
0733 const edm::Handle<std::vector<SimTrack>>& fullSimTracks = iEvent.getHandle(tok_fullSimTk_);
0734 const edm::Handle<std::vector<SimVertex>>& fullSimVertices = iEvent.getHandle(tok_fullSimVx_);
0735 mySimEvent[0]->fill(*fullSimTracks, *fullSimVertices);
0736
0737
0738
0739
0740 const edm::Handle<std::vector<SimTrack>>& fastSimTracks = iEvent.getHandle(tok_fastSimTk_);
0741 const edm::Handle<std::vector<SimVertex>>& fastSimVertices = iEvent.getHandle(tok_fastSimVx_);
0742 mySimEvent[1]->fill(*fastSimTracks, *fastSimVertices);
0743
0744
0745
0746
0747 XYZTLorentzVector theProtonMomentum(0., 0., 0., 0.939);
0748
0749
0750 if (saveNU) {
0751 if (ObjectNumber == -1 || nuEvent->nInteractions() == 1000) {
0752 ObjectNumber = TProcessID::GetObjectCount();
0753 nuEvent->reset();
0754 }
0755 }
0756
0757 for (unsigned ievt = 0; ievt < 2; ++ievt) {
0758
0759
0760
0761
0762
0763 if (!mySimEvent[ievt]->nVertices())
0764 continue;
0765 const FSimTrack& thePion = mySimEvent[ievt]->track(0);
0766
0767
0768
0769 h1[ievt]->Fill(mySimEvent[ievt]->nVertices());
0770
0771 if (mySimEvent[ievt]->nVertices() == 1) {
0772 if (thePion.trackerSurfaceMomentum().e() < 1E-10)
0773 ++stoppedPions[ievt];
0774 }
0775 if (mySimEvent[ievt]->nVertices() == 1)
0776 continue;
0777
0778 const FSimVertex& thePionVertex = mySimEvent[ievt]->vertex(1);
0779
0780 double zed = thePionVertex.position().z();
0781 double radius = thePionVertex.position().pt();
0782 double eta = thePionVertex.position().eta();
0783
0784 h0[ievt]->Fill(fabs(thePionVertex.position().z()), thePionVertex.position().pt());
0785
0786
0787
0788
0789 unsigned ndaugh = thePionVertex.nDaughters();
0790 h2[ievt]->Fill(ndaugh);
0791
0792
0793
0794
0795
0796
0797
0798 int firstDaughter = -1;
0799 int lastDaughter = -1;
0800 if (thePionVertex.nDaughters()) {
0801 lastDaughter = thePionVertex.daughters()[thePionVertex.nDaughters() - 1];
0802 firstDaughter = thePionVertex.daughters()[0];
0803 }
0804
0805
0806
0807
0808
0809
0810 XYZTLorentzVector totMoth = thePion.momentum();
0811 XYZTLorentzVector totDaugh(0., 0., 0., 0.);
0812
0813
0814 unsigned nleptons = 0;
0815 unsigned nothers = 0;
0816 if (!(firstDaughter < 0 || lastDaughter < 0)) {
0817 for (int idaugh = firstDaughter; idaugh <= lastDaughter; ++idaugh) {
0818 const FSimTrack& myDaugh = mySimEvent[ievt]->track(idaugh);
0819 totDaugh += myDaugh.momentum();
0820
0821
0822 if (abs(myDaugh.type()) == 11 || abs(myDaugh.type()) == 13)
0823 ++nleptons;
0824
0825 if (abs(myDaugh.type()) != 111 && abs(myDaugh.type()) != 211)
0826 ++nothers;
0827 }
0828 }
0829
0830
0831 if ((abs(thePion.type()) == 211 && ndaugh == 1) || (abs(thePion.type()) == 130 && ndaugh < 4) ||
0832 (abs(thePion.type()) == 321 && ndaugh < 4)) {
0833 double diffE = (totMoth - totDaugh).E();
0834 double diffP = std::sqrt((totMoth - totDaugh).Vect().Mag2());
0835 double diffm = totMoth.M2() - totDaugh.M2();
0836
0837
0838
0839 if (abs(thePion.type()) == 130 && fabs(diffE) < 1E-5 && diffP < 1E-5)
0840 return;
0841
0842
0843
0844 if (nleptons == 1)
0845 return;
0846
0847
0848
0849
0850
0851
0852
0853
0854
0855 h7[ievt]->Fill(diffE);
0856 h9[ievt]->Fill(diffm);
0857 if (abs(thePion.type()) != 211 &&
0858 diffE > -1E-5 && diffE < 0.1 &&
0859 nothers == 0)
0860 return;
0861 h6[ievt]->Fill(diffE);
0862 h8[ievt]->Fill(diffm);
0863 }
0864
0865 if (ndaugh)
0866 ++interactingPions[ievt];
0867 else
0868 ++stoppedPions[ievt];
0869
0870
0871 if (!(firstDaughter < 0 || lastDaughter < 0)) {
0872
0873 XYZTLorentzVector theBoost = thePion.momentum() + theProtonMomentum;
0874 double ecm = theBoost.mag();
0875 theBoost /= theBoost.e();
0876 if (ievt == 0 && saveNU && totalNEvt < maxNU) {
0877 NUEvent::NUInteraction interaction;
0878 interaction.first = nuEvent->nParticles();
0879 interaction.last = interaction.first + lastDaughter - firstDaughter;
0880 nuEvent->addNUInteraction(interaction);
0881 ++totalNU;
0882 }
0883
0884
0885 double qTot = 0.;
0886 double eBefore = thePion.momentum().E();
0887 double eAfter = 0.;
0888 XYZTLorentzVector theTotal(0., 0., 0., 0.);
0889
0890
0891 XYZVector zAxis(0., 0., 1.);
0892 XYZVector rotationAxis = (theBoost.Vect().Cross(zAxis)).Unit();
0893 double rotationAngle = std::acos(theBoost.Vect().Unit().Z());
0894 RawParticle::Rotation rotation(rotationAxis, rotationAngle);
0895
0896 for (int idaugh = firstDaughter; idaugh <= lastDaughter; ++idaugh) {
0897
0898 const FSimTrack& myDaugh = mySimEvent[ievt]->track(idaugh);
0899 qTot += myDaugh.charge();
0900
0901 RawParticle theMom(myDaugh.momentum());
0902 eAfter += theMom.E();
0903
0904
0905 theMom.boost(-theBoost.x(), -theBoost.y(), -theBoost.z());
0906 theTotal = theTotal + theMom.momentum();
0907
0908
0909 theMom.rotate(rotation);
0910
0911
0912 if (ievt == 0 && saveNU && totalNEvt <= maxNU) {
0913 NUEvent::NUParticle particle;
0914 particle.px = theMom.px() / ecm;
0915 particle.py = theMom.py() / ecm;
0916 particle.pz = theMom.pz() / ecm;
0917 particle.mass = theMom.mag();
0918 particle.id = myDaugh.type();
0919 nuEvent->addNUParticle(particle);
0920
0921
0922 }
0923 }
0924
0925
0926 h3[ievt]->Fill(ecm);
0927 h4[ievt]->Fill(theTotal.mag() / ecm);
0928 h5[ievt]->Fill(sqrt(theTotal.Vect().mag2()));
0929 h10[ievt]->Fill(eAfter / eBefore);
0930
0931
0932 totalCharge[ievt]->Fill(qTot);
0933
0934
0935 bool filled = false;
0936 for (unsigned hist = 0; hist < h100.size() && !filled; ++hist) {
0937 if (radius < trackerRadius[hist][ievt] && fabs(zed) < trackerLength[hist][ievt]) {
0938 h100[hist][ievt]->Fill(eta);
0939 filled = true;
0940 }
0941 }
0942
0943
0944 filled = false;
0945 for (unsigned hist = 0; hist < h200.size() && !filled; ++hist) {
0946 if (radius < blockTrackerRadius[hist][ievt] && fabs(zed) < blockTrackerLength[hist][ievt]) {
0947 h200[hist][ievt]->Fill(eta);
0948 filled = true;
0949 }
0950 }
0951
0952
0953 for (unsigned hist = 0; hist < h300.size(); ++hist) {
0954 if ((radius < subTrackerRadius[hist][ievt] && fabs(zed) < subTrackerLength[hist][ievt]) ||
0955 (hist == 2 && radius < subTrackerRadius[1][ievt] && fabs(zed) < subTrackerLength[1][ievt])) {
0956 h300[hist][ievt]->Fill(eta);
0957 }
0958 }
0959 }
0960
0961
0962 if (ievt == 0 && saveNU && totalNEvt <= maxNU) {
0963 if (nuEvent->nInteractions() == 1000) {
0964
0965 TProcessID::SetObjectCount(ObjectNumber);
0966
0967 std::cout << "Saved " << nuEvent->nInteractions() << " Interaction(s) with " << nuEvent->nParticles()
0968 << " Particles in the NUEvent " << std::endl;
0969 outFile->cd();
0970 nuTree->Fill();
0971
0972 }
0973 }
0974 }
0975 }
0976
0977
0978
0979 DEFINE_FWK_MODULE(testNuclearInteractions);