DTDigitizer/src/DTDigitizer.cc

0001 /** \file
0002  *
0003  *  \authors: G. Bevilacqua, N. Amapane, G. Cerminara, R. Bellan
0004  */
0005
0006 // system include files
0007 #include <memory>
0008
0009 // C++ headers
0010 #include <cmath>
0011
0012 // Random generator
0013 #include "FWCore/Utilities/interface/Exception.h"
0014 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0015 #include <CLHEP/Random/RandFlat.h>
0016 #include <CLHEP/Random/RandGaussQ.h>
0017
0018 // Framework
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 // Geometry
0029 #include "Geometry/DTGeometry/interface/DTLayer.h"
0030 #include "Geometry/DTGeometry/interface/DTTopology.h"
0031
0032 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0033
0034 // Digis
0035 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0036 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0037 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0038
0039 // DTDigitizer
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 // namespaces
0047 using namespace edm;
0048 using namespace std;
0049
0050 // Constructor
0051 DTDigitizer::DTDigitizer(const ParameterSet &conf_)
0052     :  // Sync Algo
0053       theSync{DTDigiSyncFactory::get()->create(conf_.getParameter<string>("SyncName"),
0054                                                conf_.getParameter<ParameterSet>("pset"))} {
0055   // Set verbose output
0056   debug = conf_.getUntrackedParameter<bool>("debug");
0057
0058   if (debug)
0059     LogPrint("DTDigitizer") << "Creating a DTDigitizer" << endl;
0060
0061   // register the Producer with a label
0062   // produces<DTDigiCollection>("MuonDTDigis"); // FIXME: Do I pass it by
0063   // ParameterSet?
0064   produces<DTDigiCollection>();  // FIXME: Do I pass it by ParameterSet?
0065   //  produces<DTDigiSimLinkCollection>("MuonDTDigiSimLinks");
0066   produces<DTDigiSimLinkCollection>();
0067
0068   // Parameters:
0069
0070   // build digis only for mu hits (for debug purposes)
0071   onlyMuHits = conf_.getParameter<bool>("onlyMuHits");
0072
0073   // interpolate parametrization function
0074   interpolate = conf_.getParameter<bool>("interpolate");
0075
0076   // Velocity of signal propagation along the wire (cm/ns)
0077   // For the default value
0078   // cfr. CMS-IN 2000-021:   (2.56+-0.17)x1e8 m/s
0079   //      CMS NOTE 2003-17:  (0.244)  m/ns
0080   vPropWire = conf_.getParameter<double>("vPropWire");  // 24.4
0081
0082   // Dead time for signals on the same wire (number from M. Pegoraro)
0083   deadTime = conf_.getParameter<double>("deadTime");  // 150
0084
0085   // further configurable smearing
0086   smearing = conf_.getParameter<double>("Smearing");  // 3.
0087
0088   // Debug flag to switch to the Ideal model
0089   // it uses a constant drift velocity and doesn't set any external delay
0090   IdealModel = conf_.getParameter<bool>("IdealModel");
0091
0092   // Flag to specify that we want digis in phase-2 units (25./30. ns)
0093   // instead that the old units (25./32.)
0094   if (conf_.getParameter<bool>("phase2Digis"))
0095     base = 30;
0096   else
0097     base = 32;
0098
0099   // Constant drift velocity needed by the above flag
0100   if (IdealModel)
0101     theConstVDrift = conf_.getParameter<double>("IdealModelConstantDriftVelocity");  // 55 um/ns
0102   else
0103     theConstVDrift = 55.;
0104
0105   // get random engine
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   // MultipleLinks=false ==> one-to-one correspondence between digis and SimHits
0112   MultipleLinks = conf_.getParameter<bool>("MultipleLinks");
0113   // MultipleLinks=true ==> association of SimHits within a time window
0114   // LinksTimeWindow (of the order of the resolution)
0115   LinksTimeWindow = conf_.getParameter<double>("LinksTimeWindow");  // (10 ns)
0116
0117   // Name of Collection used for create the XF
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   // String to choose between ideal (the default) and (mis)aligned geometry
0124   // for the digitization step
0125   geometryType = conf_.getParameter<std::string>("GeometryType");
0126   muonGeom_token = esConsumes<DTGeometry, MuonGeometryRecord>(edm::ESInputTag("", geometryType));
0127 }
0128
0129 // method called to produce the data
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   //************ 1 ***************
0138   // create the container for the SimHits
0139   //  Handle<PSimHitContainer> simHits;
0140   //  iEvent.getByLabel("g4SimHits","MuonDTHits",simHits);
0141
0142   // use MixCollection instead of the previous
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   // create the pointer to the Digi container
0149   unique_ptr<DTDigiCollection> output(new DTDigiCollection());
0150   // pointer to the DigiSimLink container
0151   unique_ptr<DTDigiSimLinkCollection> outputLinks(new DTDigiSimLinkCollection());
0152
0153   // Muon Geometry
0154   ESHandle<DTGeometry> muonGeom = iSetup.getHandle(muonGeom_token);
0155
0156   // Magnetic Field
0157   ESHandle<MagneticField> magnField = iSetup.getHandle(magnField_token);
0158
0159   //************ 2 ***************
0160
0161   // These are sorted by DetId, i.e. by layer and then by wire #
0162   //  map<DTDetId, vector<const PSimHit*> > wireMap;
0163   DTWireIdMap wireMap;
0164
0165   for (MixCollection<PSimHit>::MixItr simHit = simHits->begin(); simHit != simHits->end(); simHit++) {
0166     // Create the id of the wire, the simHits in the DT known also the wireId
0167
0168     DTWireId wireId((*simHit).detUnitId());
0169     // Fill the map
0170     wireMap[wireId].push_back(&(*simHit));
0171   }
0172
0173   pair<float, bool> time(0., false);
0174
0175   //************ 3 ***************
0176   // Loop over the wires
0177   for (DTWireIdMapConstIter wire = wireMap.begin(); wire != wireMap.end(); wire++) {
0178     // SimHit Container associated to the wire
0179     const vector<const PSimHit *> &vhit = (*wire).second;
0180     if (!vhit.empty()) {
0181       TDContainer tdCont;  // It is a vector<pair<const PSimHit*,float> >;
0182
0183       //************ 4 ***************
0184       DTWireId wireId = (*wire).first;
0185
0186       // const DTLayer* layer = dynamic_cast< const DTLayer* >
0187       // (muonGeom->idToDet(wireId.layerId()));
0188       const DTLayer *layer = muonGeom->layer(wireId.layerId());
0189
0190       // Loop on the hits of this wire
0191       for (vector<const PSimHit *>::const_iterator hit = vhit.begin(); hit != vhit.end(); hit++) {
0192         //************ 5 ***************
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         //************ 6 ***************
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       //************ 7 ***************
0209
0210       // the loading must be done by layer but
0211       // the digitization must be done by wire (in order to take into account
0212       // the dead time)
0213
0214       storeDigis(wireId, tdCont, *output, *outputLinks);
0215     }
0216   }
0217
0218   //************ 8 ***************
0219   // Load the Digi Container in the Event
0220   // iEvent.put(std::move(output),"MuonDTDigis");
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   // Pay attention: in CMSSW the rf of the SimHit is in the layer's rf
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   // The bolean is used to flag the drift time computation
0256   pair<float, bool> driftTime(0., false);
0257
0258   // if delta in gas->ignore, since it is included in the parametrisation.
0259   // FIXME: should check that it is actually a delta ray produced by a nearby
0260   // muon hit.
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   // Radius and sagitta according to direction of momentum
0272   // (just for printing)
0273   // NOTE: in cmsim, d is always taken // pHat!
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   // Radius, sagitta according to field bending
0283   // (just for printing)
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   // cos(delta), delta= angle between direction at entry and hit segment
0296   // (just for printing)
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   // Select cases where parametrization can not be used.
0312   bool noParametrisation = ((entrySide == DTTopology::none || exitSide == DTTopology::none)  // case # 2,3,8,9 or 11
0313                             || (entrySide == exitSide)                                       // case # 4 or 10
0314                             || ((entrySide == DTTopology::xMin && exitSide == DTTopology::xMax) ||
0315                                 (entrySide == DTTopology::xMax && exitSide == DTTopology::xMin))  // Hit is case # 7
0316   );
0317
0318   // FIXME: now, debug warning only; consider treating those
0319   // with TM algo.
0320   if (delta < 0.99996  // Track is not straight. FIXME: use sagitta?
0321       && (noParametrisation == false)) {
0322     if (debug)
0323       LogPrint("DTDigitizer") << "*** WARNING: hit is not straight, type = " << partType << endl;
0324   }
0325
0326   //************ 5A ***************
0327
0328   if (!noParametrisation) {
0329     LocalVector dir = hit->momentumAtEntry();  // ex Measurement3DVector dir =
0330                                                // hit->measurementDirection(); //FIXME?
0331     float theta = atan(dir.x() / -dir.z()) * 180 / M_PI;
0332
0333     // FIXME: use dir if M.S. is included as GARFIELD option...
0334     //        otherwise use hit segment dirction.
0335     //    LocalVector dir0 = (exitP-entryP).unit();
0336     //    float theta = atan(dir0.x()/-dir0.z())*180/M_PI;
0337     float x;
0338
0339     Local3DPoint pt = hit->localPosition();  // ex Measurement3DPoint pt =
0340                                              // hit->measurementPosition(); // FIXME?
0341
0342     if (fabs(pt.z()) < 0.002) {
0343       // hit center within 20 um from z=0, no need to extrapolate.
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     // Parametrisation not applicable, or failed. Use time map.
0357     driftTime = driftTimeFromTimeMap();
0358   }
0359
0360   //************ 5B ***************
0361
0362   // Signal propagation, TOF etc.
0363   if (driftTime.second) {
0364     driftTime.first += externalDelays(layer, wireId, hit);
0365   }
0366   return driftTime;
0367 }
0368
0369 //************ 5A ***************
0370
0371 pair<float, bool> DTDigitizer::driftTimeFromParametrization(
0372     float x, float theta, float By, float Bz, CLHEP::HepRandomEngine *engine) const {
0373   // Convert from CMSSW frame/units r.f. to parametrization ones.
0374   x *= 10.;  // cm -> mm
0375
0376   // FIXME: Current parametrisation can extrapolate above 21 mm,
0377   // however a detailed study is needed before using this.
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   // Different r.f. of the parametrization:
0385   // X_par = X_ORCA; Y_par=Z_ORCA; Z_par = -Y_ORCA
0386
0387   float By_par = Bz;   // Bnorm
0388   float Bz_par = -By;  // Bwire
0389   float theta_par = theta;
0390
0391   // Parametrisation uses interpolation up to |theta|=45 deg,
0392   // |Bwire|=0.4, |Bnorm|=0.75; extrapolation above.
0393   if (fabs(theta_par) > 45.) {
0394     if (debug)
0395       LogPrint("DTDigitizer") << "*** WARNING: extrapolating theta > 45: " << theta << endl;
0396     // theta_par = min(fabs(theta_par),45.f)*((theta_par<0.)?-1.:1.);
0397   }
0398   if (fabs(By_par) > 0.75) {
0399     if (debug)
0400       LogPrint("DTDigitizer") << "*** WARNING: extrapolating Bnorm > 0.75: " << By_par << endl;
0401     // By_par = min(fabs(By_par),0.75f)*((By_par<0.)?-1.:1.);
0402   }
0403   if (fabs(Bz_par) > 0.4) {
0404     if (debug)
0405       LogPrint("DTDigitizer") << "*** WARNING: extrapolating Bwire >0.4: " << Bz_par << endl;
0406     // Bz_par = min(fabs(Bz_par),0.4)*((Bz_par<0.)?-1.:1.);
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   // Double half-gaussian smearing
0425   float time = asymGausSmear(DT.t_drift, DT.t_width_m, DT.t_width_p, engine);
0426
0427   // Do not allow the smearing to lead to negative values
0428   time = max(time, 0.f);
0429
0430   // Apply a Gaussian smearing to account for electronic effects (cf. 2004 TB
0431   // analysis) The width of the Gaussian can be configured with the "Smearing"
0432   // parameter
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   // FIXME: not yet implemented.
0462   if (debug)
0463     LogPrint("DTDigitizer") << "   TimeMap " << endl;
0464   return pair<float, bool>(0., false);
0465 }
0466
0467 //************ 5B ***************
0468
0469 float DTDigitizer::externalDelays(const DTLayer *layer, const DTWireId &wireId, const PSimHit *hit) const {
0470   // Time of signal propagation along wire.
0471
0472   float wireCoord = hit->localPosition().y();
0473   float halfL = (layer->specificTopology().cellLenght()) / 2.;
0474   float propgL = halfL - wireCoord;  // the FE is always located at the pos coord.
0475
0476   float propDelay = propgL / vPropWire;
0477
0478   // Real TOF.
0479   float tof = hit->tof();
0480
0481   // Delays and t0 according to theSync
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 // accumulate digis by layer
0493
0494 void DTDigitizer::storeDigis(DTWireId &wireId,
0495                              TDContainer &hits,
0496                              DTDigiCollection &output,
0497                              DTDigiSimLinkCollection &outputLinks) {
0498   //************ 7A ***************
0499
0500   // sort signal times
0501   sort(hits.begin(), hits.end(), hitLessT());
0502
0503   //************ 7B ***************
0504
0505   float wakeTime = -999999.0;
0506   float resolTime = -999999.0;
0507   int digiN = -1;  // Digi identifier within the cell (for multiple digis)
0508   DTDigi digi;
0509
0510   // loop over signal times and drop signals inside dead time
0511   for (TDContainer::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0512     if (onlyMuHits && abs((*hit).first->particleType()) != 13)
0513       continue;
0514
0515     //************ 7C ***************
0516
0517     float time = (*hit).second;
0518     if (time > wakeTime) {
0519       // Note that digi is constructed with a float value (in ns)
0520       int wireN = wireId.wire();
0521       digiN++;
0522       digi = DTDigi(wireN, time, digiN, base);
0523
0524       // Add association between THIS digi and the corresponding SimTrack
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       //************ 7D ***************
0538       DTLayerId layerID = wireId.layerId();  // taking the layer of the wire
0539       output.insertDigi(layerID, digi);      // ordering Digis by layer
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   //  ProcessTypeEnumerator pTypes;
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                           // << "   packedTrackId         = " << hit->packedTrackId() << endl
0576                           << "   trackId               = " << hit->trackId()
0577                           << endl  // new,is the same as the
0578                                    // previous?? FIXME-Check
0579                           << "   |p|                   = " << hit->pabs() << endl
0580                           << "   Energy loss           = " << hit->energyLoss()
0581                           << endl
0582                           // << "   timeOffset            = " << hit->timeOffset() << endl
0583                           // << "   measurementPosition   = " << hit->measurementPosition() << endl
0584                           // << "   measurementDirection  = " << hit->measurementDirection() << endl
0585                           // //FIXME
0586                           << "   localDirection        = " << hit->momentumAtEntry().unit()
0587                           << endl  // FIXME is it a versor?
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 }