File indexing completed on 2025-04-04 01:27:01
0001
0002
0003
0004
0005
0006
0007 #include <memory>
0008
0009
0010 #include <cmath>
0011
0012
0013 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include <CLHEP/Random/RandFlat.h>
0016 #include <CLHEP/Random/RandGaussQ.h>
0017
0018
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027
0028
0029 #include "Geometry/DTGeometry/interface/DTLayer.h"
0030 #include "Geometry/DTGeometry/interface/DTTopology.h"
0031
0032 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0033
0034
0035 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0036 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0037 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0038
0039
0040 #include "SimMuon/DTDigitizer/interface/DTDigiSyncBase.h"
0041 #include "SimMuon/DTDigitizer/interface/DTDigiSyncFactory.h"
0042
0043 #include "SimMuon/DTDigitizer/src/DTDigitizer.h"
0044 #include "SimMuon/DTDigitizer/src/DTDriftTimeParametrization.h"
0045
0046
0047 using namespace edm;
0048 using namespace std;
0049
0050
0051 DTDigitizer::DTDigitizer(const ParameterSet &conf_)
0052 :
0053 theSync{DTDigiSyncFactory::get()->create(conf_.getParameter<string>("SyncName"),
0054 conf_.getParameter<ParameterSet>("pset"))} {
0055
0056 debug = conf_.getUntrackedParameter<bool>("debug");
0057
0058 if (debug)
0059 LogPrint("DTDigitizer") << "Creating a DTDigitizer" << endl;
0060
0061
0062
0063
0064 produces<DTDigiCollection>();
0065
0066 produces<DTDigiSimLinkCollection>();
0067
0068
0069
0070
0071 onlyMuHits = conf_.getParameter<bool>("onlyMuHits");
0072
0073
0074 interpolate = conf_.getParameter<bool>("interpolate");
0075
0076
0077
0078
0079
0080 vPropWire = conf_.getParameter<double>("vPropWire");
0081
0082
0083 deadTime = conf_.getParameter<double>("deadTime");
0084
0085
0086 smearing = conf_.getParameter<double>("Smearing");
0087
0088
0089
0090 IdealModel = conf_.getParameter<bool>("IdealModel");
0091
0092
0093
0094 if (conf_.getParameter<bool>("phase2Digis"))
0095 base = 30;
0096 else
0097 base = 32;
0098
0099
0100 if (IdealModel)
0101 theConstVDrift = conf_.getParameter<double>("IdealModelConstantDriftVelocity");
0102 else
0103 theConstVDrift = 55.;
0104
0105
0106 edm::Service<edm::RandomNumberGenerator> rng;
0107 if (!rng.isAvailable()) {
0108 throw cms::Exception("Configuration") << "RandomNumberGeneratorService for DTDigitizer missing in cfg file";
0109 }
0110
0111
0112 MultipleLinks = conf_.getParameter<bool>("MultipleLinks");
0113
0114
0115 LinksTimeWindow = conf_.getParameter<double>("LinksTimeWindow");
0116
0117
0118 mix_ = conf_.getParameter<std::string>("mixLabel");
0119 collection_for_XF = conf_.getParameter<std::string>("InputCollection");
0120 cf_token = consumes<CrossingFrame<PSimHit>>(edm::InputTag(mix_, collection_for_XF));
0121 magnField_token = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0122
0123
0124
0125 geometryType = conf_.getParameter<std::string>("GeometryType");
0126 muonGeom_token = esConsumes<DTGeometry, MuonGeometryRecord>(edm::ESInputTag("", geometryType));
0127 }
0128
0129
0130 void DTDigitizer::produce(Event &iEvent, const EventSetup &iSetup) {
0131 edm::Service<edm::RandomNumberGenerator> rng;
0132 CLHEP::HepRandomEngine *engine = &rng->getEngine(iEvent.streamID());
0133
0134 if (debug)
0135 LogPrint("DTDigitizer") << "--- Run: " << iEvent.id().run() << " Event: " << iEvent.id().event() << endl;
0136
0137
0138
0139
0140
0141
0142
0143 Handle<CrossingFrame<PSimHit>> xFrame;
0144 iEvent.getByToken(cf_token, xFrame);
0145
0146 unique_ptr<MixCollection<PSimHit>> simHits(new MixCollection<PSimHit>(xFrame.product()));
0147
0148
0149 unique_ptr<DTDigiCollection> output(new DTDigiCollection());
0150
0151 unique_ptr<DTDigiSimLinkCollection> outputLinks(new DTDigiSimLinkCollection());
0152
0153
0154 ESHandle<DTGeometry> muonGeom = iSetup.getHandle(muonGeom_token);
0155
0156
0157 ESHandle<MagneticField> magnField = iSetup.getHandle(magnField_token);
0158
0159
0160
0161
0162
0163 DTWireIdMap wireMap;
0164
0165 for (MixCollection<PSimHit>::MixItr simHit = simHits->begin(); simHit != simHits->end(); simHit++) {
0166
0167
0168 DTWireId wireId((*simHit).detUnitId());
0169
0170 wireMap[wireId].push_back(&(*simHit));
0171 }
0172
0173 pair<float, bool> time(0., false);
0174
0175
0176
0177 for (DTWireIdMapConstIter wire = wireMap.begin(); wire != wireMap.end(); wire++) {
0178
0179 const vector<const PSimHit *> &vhit = (*wire).second;
0180 if (!vhit.empty()) {
0181 TDContainer tdCont;
0182
0183
0184 DTWireId wireId = (*wire).first;
0185
0186
0187
0188 const DTLayer *layer = muonGeom->layer(wireId.layerId());
0189
0190
0191 for (vector<const PSimHit *>::const_iterator hit = vhit.begin(); hit != vhit.end(); hit++) {
0192
0193 LocalPoint locPos = (*hit)->localPosition();
0194
0195 const LocalVector BLoc = layer->surface().toLocal(magnField->inTesla(layer->surface().toGlobal(locPos)));
0196
0197 time = computeTime(layer, wireId, *hit, BLoc, engine);
0198
0199
0200 if (time.second) {
0201 tdCont.push_back(make_pair((*hit), time.first));
0202 } else {
0203 if (debug)
0204 LogPrint("DTDigitizer") << "hit discarded" << endl;
0205 }
0206 }
0207
0208
0209
0210
0211
0212
0213
0214 storeDigis(wireId, tdCont, *output, *outputLinks);
0215 }
0216 }
0217
0218
0219
0220
0221 iEvent.put(std::move(output));
0222 iEvent.put(std::move(outputLinks));
0223 }
0224
0225 pair<float, bool> DTDigitizer::computeTime(const DTLayer *layer,
0226 const DTWireId &wireId,
0227 const PSimHit *hit,
0228 const LocalVector &BLoc,
0229 CLHEP::HepRandomEngine *engine) {
0230 LocalPoint entryP = hit->entryPoint();
0231 LocalPoint exitP = hit->exitPoint();
0232 int partType = hit->particleType();
0233
0234 const DTTopology &topo = layer->specificTopology();
0235
0236
0237
0238 if (debug)
0239 LogPrint("DTDigitizer") << "Hit local entry point: " << entryP << endl << "Hit local exit point: " << exitP << endl;
0240
0241 float xwire = topo.wirePosition(wireId.wire());
0242 float xEntry = entryP.x() - xwire;
0243 float xExit = exitP.x() - xwire;
0244
0245 if (debug)
0246 LogPrint("DTDigitizer") << "wire position: " << xwire << " x entry in cell rf: " << xEntry
0247 << " x exit in cell rf: " << xExit << endl;
0248
0249 DTTopology::Side entrySide = topo.onWhichBorder(xEntry, entryP.y(), entryP.z());
0250 DTTopology::Side exitSide = topo.onWhichBorder(xExit, exitP.y(), exitP.z());
0251
0252 if (debug)
0253 dumpHit(hit, xEntry, xExit, topo);
0254
0255
0256 pair<float, bool> driftTime(0., false);
0257
0258
0259
0260
0261
0262 if (partType == 11 && entrySide == DTTopology::none) {
0263 if (debug)
0264 LogPrint("DTDigitizer") << " e- hit in gas; discarding " << endl;
0265 return driftTime;
0266 }
0267
0268 float By = BLoc.y();
0269 float Bz = BLoc.z();
0270
0271
0272
0273
0274 LocalVector d = (exitP - entryP);
0275 LocalVector pHat = hit->localDirection().unit();
0276 LocalVector hHat = (d.cross(pHat.cross(d))).unit();
0277 float cosAlpha = hHat.dot(pHat);
0278 float sinAlpha = sqrt(1. - cosAlpha * cosAlpha);
0279 float radius_P = (d.mag()) / (2. * cosAlpha);
0280 float sagitta_P = radius_P * (1. - sinAlpha);
0281
0282
0283
0284 float halfd = d.mag() / 2.;
0285 float BMag = BLoc.mag();
0286 LocalVector pT = (pHat - (BLoc.unit() * pHat.dot(BLoc.unit()))) * (hit->pabs());
0287 float radius_B = (pT.mag() / (0.3 * BMag)) * 100.;
0288 float sagitta_B;
0289 if (radius_B > halfd) {
0290 sagitta_B = radius_B - sqrt(radius_B * radius_B - halfd * halfd);
0291 } else {
0292 sagitta_B = radius_B;
0293 }
0294
0295
0296
0297 float delta = pHat.dot(d.unit());
0298 if (debug)
0299 LogPrint("DTDigitizer") << " delta = " << delta << endl
0300 << " cosAlpha = " << cosAlpha << endl
0301 << " sinAlpha = " << sinAlpha << endl
0302 << " pMag = " << pT.mag() << endl
0303 << " bMag = " << BMag << endl
0304 << " pT = " << pT << endl
0305 << " halfd = " << halfd << endl
0306 << " radius_P (cm) = " << radius_P << endl
0307 << " sagitta_P (um) = " << sagitta_P * 10000. << endl
0308 << " radius_B (cm) = " << radius_B << endl
0309 << " sagitta_B (um) = " << sagitta_B * 10000. << endl;
0310
0311
0312 bool noParametrisation = ((entrySide == DTTopology::none || exitSide == DTTopology::none)
0313 || (entrySide == exitSide)
0314 || ((entrySide == DTTopology::xMin && exitSide == DTTopology::xMax) ||
0315 (entrySide == DTTopology::xMax && exitSide == DTTopology::xMin))
0316 );
0317
0318
0319
0320 if (delta < 0.99996
0321 && (noParametrisation == false)) {
0322 if (debug)
0323 LogPrint("DTDigitizer") << "*** WARNING: hit is not straight, type = " << partType << endl;
0324 }
0325
0326
0327
0328 if (!noParametrisation) {
0329 LocalVector dir = hit->momentumAtEntry();
0330
0331 float theta = atan(dir.x() / -dir.z()) * 180 / M_PI;
0332
0333
0334
0335
0336
0337 float x;
0338
0339 Local3DPoint pt = hit->localPosition();
0340
0341
0342 if (fabs(pt.z()) < 0.002) {
0343
0344 x = pt.x() - xwire;
0345 } else {
0346 x = xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z());
0347 }
0348
0349 if (IdealModel)
0350 return make_pair(fabs(x) / theConstVDrift, true);
0351 else
0352 driftTime = driftTimeFromParametrization(x, theta, By, Bz, engine);
0353 }
0354
0355 if ((driftTime.second) == false) {
0356
0357 driftTime = driftTimeFromTimeMap();
0358 }
0359
0360
0361
0362
0363 if (driftTime.second) {
0364 driftTime.first += externalDelays(layer, wireId, hit);
0365 }
0366 return driftTime;
0367 }
0368
0369
0370
0371 pair<float, bool> DTDigitizer::driftTimeFromParametrization(
0372 float x, float theta, float By, float Bz, CLHEP::HepRandomEngine *engine) const {
0373
0374 x *= 10.;
0375
0376
0377
0378 if (fabs(x) > 21.) {
0379 if (debug)
0380 LogPrint("DTDigitizer") << "*** WARNING: parametrisation: x out of range = " << x << ", skipping" << endl;
0381 return pair<float, bool>(0.f, false);
0382 }
0383
0384
0385
0386
0387 float By_par = Bz;
0388 float Bz_par = -By;
0389 float theta_par = theta;
0390
0391
0392
0393 if (fabs(theta_par) > 45.) {
0394 if (debug)
0395 LogPrint("DTDigitizer") << "*** WARNING: extrapolating theta > 45: " << theta << endl;
0396
0397 }
0398 if (fabs(By_par) > 0.75) {
0399 if (debug)
0400 LogPrint("DTDigitizer") << "*** WARNING: extrapolating Bnorm > 0.75: " << By_par << endl;
0401
0402 }
0403 if (fabs(Bz_par) > 0.4) {
0404 if (debug)
0405 LogPrint("DTDigitizer") << "*** WARNING: extrapolating Bwire >0.4: " << Bz_par << endl;
0406
0407 }
0408
0409 DTDriftTimeParametrization::drift_time DT;
0410 static const DTDriftTimeParametrization par;
0411 unsigned short flag = par.MB_DT_drift_time(x, theta_par, By_par, Bz_par, 0, &DT, interpolate);
0412
0413 if (debug) {
0414 LogPrint("DTDigitizer") << " Parametrisation: x, theta, Bnorm, Bwire = " << x << " " << theta_par << " "
0415 << By_par << " " << Bz_par << endl
0416 << " time=" << DT.t_drift << " sigma_m=" << DT.t_width_m << " sigma_p=" << DT.t_width_p
0417 << endl;
0418 if (flag != 1) {
0419 LogPrint("DTDigitizer") << "*** WARNING: call to parametrisation failed" << endl;
0420 return pair<float, bool>(0.f, false);
0421 }
0422 }
0423
0424
0425 float time = asymGausSmear(DT.t_drift, DT.t_width_m, DT.t_width_p, engine);
0426
0427
0428 time = max(time, 0.f);
0429
0430
0431
0432
0433
0434 double u = CLHEP::RandGaussQ::shoot(engine, 0., smearing);
0435 time += u;
0436
0437 if (debug)
0438 LogPrint("DTDigitizer") << " drift time = " << time << endl;
0439
0440 return pair<float, bool>(time, true);
0441 }
0442
0443 float DTDigitizer::asymGausSmear(double mean,
0444 double sigmaLeft,
0445 double sigmaRight,
0446 CLHEP::HepRandomEngine *engine) const {
0447 double f = sigmaLeft / (sigmaLeft + sigmaRight);
0448 double t;
0449
0450 if (CLHEP::RandFlat::shoot(engine) <= f) {
0451 t = CLHEP::RandGaussQ::shoot(engine, mean, sigmaLeft);
0452 t = mean - fabs(t - mean);
0453 } else {
0454 t = CLHEP::RandGaussQ::shoot(engine, mean, sigmaRight);
0455 t = mean + fabs(t - mean);
0456 }
0457 return static_cast<float>(t);
0458 }
0459
0460 pair<float, bool> DTDigitizer::driftTimeFromTimeMap() const {
0461
0462 if (debug)
0463 LogPrint("DTDigitizer") << " TimeMap " << endl;
0464 return pair<float, bool>(0., false);
0465 }
0466
0467
0468
0469 float DTDigitizer::externalDelays(const DTLayer *layer, const DTWireId &wireId, const PSimHit *hit) const {
0470
0471
0472 float wireCoord = hit->localPosition().y();
0473 float halfL = (layer->specificTopology().cellLenght()) / 2.;
0474 float propgL = halfL - wireCoord;
0475
0476 float propDelay = propgL / vPropWire;
0477
0478
0479 float tof = hit->tof();
0480
0481
0482
0483 double sync = theSync->digitizerOffset(&wireId, layer);
0484
0485 if (debug) {
0486 LogPrint("DTDigitizer") << " propDelay =" << propDelay << "; TOF=" << tof << "; sync= " << sync << endl;
0487 }
0488
0489 return propDelay + tof + sync;
0490 }
0491
0492
0493
0494 void DTDigitizer::storeDigis(DTWireId &wireId,
0495 TDContainer &hits,
0496 DTDigiCollection &output,
0497 DTDigiSimLinkCollection &outputLinks) {
0498
0499
0500
0501 sort(hits.begin(), hits.end(), hitLessT());
0502
0503
0504
0505 float wakeTime = -999999.0;
0506 float resolTime = -999999.0;
0507 int digiN = -1;
0508 DTDigi digi;
0509
0510
0511 for (TDContainer::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0512 if (onlyMuHits && abs((*hit).first->particleType()) != 13)
0513 continue;
0514
0515
0516
0517 float time = (*hit).second;
0518 if (time > wakeTime) {
0519
0520 int wireN = wireId.wire();
0521 digiN++;
0522 digi = DTDigi(wireN, time, digiN, base);
0523
0524
0525 unsigned int SimTrackId = (*hit).first->trackId();
0526 EncodedEventId evId = (*hit).first->eventId();
0527 DTDigiSimLink digisimLink(wireN, digiN, time, SimTrackId, evId, base);
0528
0529 if (debug) {
0530 LogPrint("DTDigitizer") << endl << "---- DTDigitizer ----" << endl;
0531 LogPrint("DTDigitizer") << "wireId: " << wireId << endl;
0532 LogPrint("DTDigitizer") << "sim. time = " << time << endl;
0533 LogPrint("DTDigitizer") << "digi number = " << digi.number() << ", digi time = " << digi.time()
0534 << ", linked to SimTrack Id = " << SimTrackId << endl;
0535 }
0536
0537
0538 DTLayerId layerID = wireId.layerId();
0539 output.insertDigi(layerID, digi);
0540 outputLinks.insertDigi(layerID, digisimLink);
0541 wakeTime = time + deadTime;
0542 resolTime = time + LinksTimeWindow;
0543 } else if (MultipleLinks && time < resolTime) {
0544 int wireN = wireId.wire();
0545 unsigned int SimTrackId = (*hit).first->trackId();
0546 EncodedEventId evId = (*hit).first->eventId();
0547 DTDigiSimLink digisimLink(wireN, digiN, time, SimTrackId, evId, base);
0548 DTLayerId layerID = wireId.layerId();
0549 outputLinks.insertDigi(layerID, digisimLink);
0550
0551 if (debug) {
0552 LogPrint("DTDigitizer") << "\nAdded multiple link: \n"
0553 << "digi number = " << digi.number() << ", digi time = " << digi.time()
0554 << " (sim. time = " << time << ")"
0555 << ", linked to SimTrack Id = " << SimTrackId << endl;
0556 }
0557 }
0558 }
0559 }
0560
0561 void DTDigitizer::dumpHit(const PSimHit *hit, float xEntry, float xExit, const DTTopology &topo) {
0562 LocalPoint entryP = hit->entryPoint();
0563 LocalPoint exitP = hit->exitPoint();
0564
0565 DTTopology::Side entrySide = topo.onWhichBorder(xEntry, entryP.y(), entryP.z());
0566 DTTopology::Side exitSide = topo.onWhichBorder(xExit, exitP.y(), exitP.z());
0567
0568
0569 LogPrint("DTDigitizer") << endl
0570 << "------- SimHit: " << endl
0571 << " Particle type = " << hit->particleType() << endl
0572 << " process type = " << hit->processType() << endl
0573 << " process type = " << hit->processType()
0574 << endl
0575
0576 << " trackId = " << hit->trackId()
0577 << endl
0578
0579 << " |p| = " << hit->pabs() << endl
0580 << " Energy loss = " << hit->energyLoss()
0581 << endl
0582
0583
0584
0585
0586 << " localDirection = " << hit->momentumAtEntry().unit()
0587 << endl
0588 << " Entry point = " << entryP << " cell x = " << xEntry << endl
0589 << " Exit point = " << exitP << " cell x = " << xExit << endl
0590 << " DR = = " << (exitP - entryP).mag() << endl
0591 << " Dx = = " << (exitP - entryP).x() << endl
0592 << " Cell w,h,l = (" << topo.cellWidth() << " , " << topo.cellHeight() << " , "
0593 << topo.cellLenght() << ") cm" << endl
0594 << " DY entry from edge = " << topo.cellLenght() / 2. - fabs(entryP.y())
0595 << " DY exit from edge = " << topo.cellLenght() / 2. - fabs(exitP.y())
0596 << " entrySide = " << (int)entrySide << " ; exitSide = " << (int)exitSide << endl;
0597 }