SiPixelDigitizer/plugins/SiPixelDigitizer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    SiPixelDigitizer
0004 // Class:      SiPixelDigitizer
0005 //
0006 /**\class SiPixelDigitizer SiPixelDigitizer.cc SimTracker/SiPixelDigitizer/src/SiPixelDigitizer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Michele Pioppi-INFN perugia
0015 //   Modifications: Freya Blekman - Cornell University
0016 //         Created:  Mon Sep 26 11:08:32 CEST 2005
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <set>
0023
0024 // user include files
0025 #include "SiPixelDigitizer.h"
0026 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelSimHitExtraInfo.h"
0027 #include "PixelDigiAddTempInfo.h"
0028 #include "SiPixelDigitizerAlgorithm.h"
0029
0030 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0031 #include "DataFormats/Common/interface/Handle.h"
0032 #include "FWCore/Framework/interface/ConsumesCollector.h"
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0036 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0037 #include "DataFormats/Common/interface/DetSet.h"
0038 #include "DataFormats/Common/interface/DetSetVector.h"
0039 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0040 #include "DataFormats/SiPixelDigi/interface/PixelDigiCollection.h"
0041 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0042 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0043 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0046 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0047
0048 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0049 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0050
0051 // user include files
0052 #include "FWCore/Framework/interface/Frameworkfwd.h"
0053
0054 #include "FWCore/Framework/interface/Event.h"
0055 #include "FWCore/Framework/interface/MakerMacros.h"
0056
0057 #include "MagneticField/Engine/interface/MagneticField.h"
0058
0059 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0060 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0061 #include "DataFormats/SiPixelDetId/interface/PixelFEDChannel.h"
0062 //Random Number
0063 #include "FWCore/ServiceRegistry/interface/Service.h"
0064 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0065 #include "FWCore/Utilities/interface/StreamID.h"
0066 #include "FWCore/Utilities/interface/Exception.h"
0067
0068 //
0069 // constants, enums and typedefs
0070 //
0071
0072 //
0073 // static data member definitions
0074 //
0075
0076 //
0077 // constructors and destructor
0078 //
0079 //using namespace std;
0080
0081 namespace cms {
0082   SiPixelDigitizer::SiPixelDigitizer(const edm::ParameterSet& iConfig,
0083                                      edm::ProducesCollector producesCollector,
0084                                      edm::ConsumesCollector& iC)
0085       : firstInitializeEvent_(true),
0086         firstFinalizeEvent_(true),
0087         applyLateReweighting_(
0088             iConfig.exists("applyLateReweighting") ? iConfig.getParameter<bool>("applyLateReweighting") : false),
0089         store_SimHitEntryExitPoints_(iConfig.exists("store_SimHitEntryExitPoints")
0090                                          ? iConfig.getParameter<bool>("store_SimHitEntryExitPoints")
0091                                          : false),
0092         _pixeldigialgo(),
0093         hitsProducer(iConfig.getParameter<std::string>("hitsProducer")),
0094         trackerContainers(iConfig.getParameter<std::vector<std::string> >("RoutList")),
0095         pilotBlades(iConfig.exists("enablePilotBlades") ? iConfig.getParameter<bool>("enablePilotBlades") : false),
0096         NumberOfEndcapDisks(iConfig.exists("NumPixelEndcap") ? iConfig.getParameter<int>("NumPixelEndcap") : 2),
0097         tTopoToken_(iC.esConsumes()),
0098         pDDToken_(iC.esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("PixGeometryType")))),
0099         pSetupToken_(iC.esConsumes()) {
0100     edm::LogInfo("PixelDigitizer ") << "Enter the Pixel Digitizer";
0101
0102     edm::LogInfo("PixelDigitizer ") << " applyLateReweighting_ " << applyLateReweighting_;
0103
0104     const std::string alias("simSiPixelDigis");
0105
0106     producesCollector.produces<edm::DetSetVector<PixelDigi> >().setBranchAlias(alias);
0107     producesCollector.produces<edm::DetSetVector<PixelDigiSimLink> >().setBranchAlias(alias + "siPixelDigiSimLink");
0108     if (store_SimHitEntryExitPoints_)
0109       producesCollector.produces<edm::DetSetVector<PixelSimHitExtraInfo> >().setBranchAlias(alias +
0110                                                                                             "siPixelExtraSimHit");
0111
0112     for (auto const& trackerContainer : trackerContainers) {
0113       edm::InputTag tag(hitsProducer, trackerContainer);
0114       iC.consumes<std::vector<PSimHit> >(edm::InputTag(hitsProducer, trackerContainer));
0115     }
0116     edm::Service<edm::RandomNumberGenerator> rng;
0117     if (!rng.isAvailable()) {
0118       throw cms::Exception("Configuration")
0119           << "SiPixelDigitizer requires the RandomNumberGeneratorService\n"
0120              "which is not present in the configuration file.  You must add the service\n"
0121              "in the configuration file or remove the modules that require it.";
0122     }
0123
0124     _pixeldigialgo = std::make_unique<SiPixelDigitizerAlgorithm>(iConfig, iC);
0125     if (NumberOfEndcapDisks != 2)
0126       producesCollector.produces<PixelFEDChannelCollection>();
0127   }
0128
0129   SiPixelDigitizer::~SiPixelDigitizer() { edm::LogInfo("PixelDigitizer ") << "Destruct the Pixel Digitizer"; }
0130
0131   //
0132   // member functions
0133   //
0134
0135   void SiPixelDigitizer::accumulatePixelHits(edm::Handle<std::vector<PSimHit> > hSimHits,
0136                                              size_t globalSimHitIndex,
0137                                              const unsigned int tofBin,
0138                                              edm::EventSetup const& iSetup) {
0139     if (hSimHits.isValid()) {
0140       std::set<unsigned int> detIds;
0141       std::vector<PSimHit> const& simHits = *hSimHits.product();
0142       const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0143       _pixeldigialgo->fillSimHitMaps(simHits, tofBin);
0144       for (std::vector<PSimHit>::const_iterator it = simHits.begin(), itEnd = simHits.end(); it != itEnd;
0145            ++it, ++globalSimHitIndex) {
0146         unsigned int detId = (*it).detUnitId();
0147         if (detIds.insert(detId).second) {
0148           // The insert succeeded, so this detector element has not yet been processed.
0149           assert(detectorUnits[detId]);
0150           if (detectorUnits[detId] &&
0151               detectorUnits[detId]
0152                   ->type()
0153                   .isTrackerPixel()) {  // this test could be avoided and changed into a check of pixdet!=0
0154             std::map<unsigned int, PixelGeomDetUnit const*>::iterator itDet = detectorUnits.find(detId);
0155             if (itDet == detectorUnits.end())
0156               continue;
0157             auto pixdet = itDet->second;
0158             assert(pixdet != nullptr);
0159             //access to magnetic field in global coordinates
0160             GlobalVector bfield = pSetup->inTesla(pixdet->surface().position());
0161             LogDebug("PixelDigitizer ") << "B-field(T) at " << pixdet->surface().position()
0162                                         << "(cm): " << pSetup->inTesla(pixdet->surface().position());
0163             _pixeldigialgo->accumulateSimHits(
0164                 it, itEnd, globalSimHitIndex, tofBin, pixdet, bfield, tTopo, randomEngine_);
0165           }
0166         }
0167       }
0168     }
0169   }
0170
0171   void SiPixelDigitizer::initializeEvent(edm::Event const& e, edm::EventSetup const& iSetup) {
0172     if (firstInitializeEvent_) {
0173       _pixeldigialgo->init(iSetup);
0174       firstInitializeEvent_ = false;
0175     }
0176
0177     // Make sure that the first crossing processed starts indexing the sim hits from zero.
0178     // This variable is used so that the sim hits from all crossing frames have sequential
0179     // indices used to create the digi-sim link (if configured to do so) rather than starting
0180     // from zero for each crossing.
0181     crossingSimHitIndexOffset_.clear();
0182
0183     // Cache random number engine
0184     edm::Service<edm::RandomNumberGenerator> rng;
0185     randomEngine_ = &rng->getEngine(e.streamID());
0186
0187     _pixeldigialgo->initializeEvent();
0188     pDD = &iSetup.getData(pDDToken_);
0189     pSetup = &iSetup.getData(pSetupToken_);
0190     const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0191
0192     // FIX THIS! We only need to clear and (re)fill this map when the geometry type IOV changes.  Use ESWatcher to determine this.
0193     if (true) {  // Replace with ESWatcher
0194       detectorUnits.clear();
0195       for (const auto& iu : pDD->detUnits()) {
0196         unsigned int detId = iu->geographicalId().rawId();
0197         if (iu->type().isTrackerPixel()) {
0198           auto pixdet = dynamic_cast<const PixelGeomDetUnit*>(iu);
0199           assert(pixdet != nullptr);
0200           if (iu->subDetector() ==
0201               GeomDetEnumerators::SubDetector::PixelEndcap) {  // true ONLY for the phase 0 pixel deetctor
0202             unsigned int disk = tTopo->layer(detId);           // using the generic layer method
0203             //if using pilot blades, then allowing it for current detector only
0204             if ((disk == 3) && ((!pilotBlades) && (NumberOfEndcapDisks == 2)))
0205               continue;
0206           }
0207           detectorUnits.insert(std::make_pair(detId, pixdet));
0208         }
0209       }
0210     }
0211   }
0212
0213   void SiPixelDigitizer::accumulate(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0214     // Step A: Get Inputs
0215     for (vstring::const_iterator i = trackerContainers.begin(), iEnd = trackerContainers.end(); i != iEnd; ++i) {
0216       edm::Handle<std::vector<PSimHit> > simHits;
0217       edm::InputTag tag(hitsProducer, *i);
0218
0219       iEvent.getByLabel(tag, simHits);
0220       unsigned int tofBin = PixelDigiSimLink::LowTof;
0221       if ((*i).find(std::string("HighTof")) != std::string::npos)
0222         tofBin = PixelDigiSimLink::HighTof;
0223       accumulatePixelHits(simHits, crossingSimHitIndexOffset_[tag.encode()], tofBin, iSetup);
0224       // Now that the hits have been processed, I'll add the amount of hits in this crossing on to
0225       // the global counter. Next time accumulateStripHits() is called it will count the sim hits
0226       // as though they were on the end of this collection.
0227       // Note that this is only used for creating digi-sim links (if configured to do so).
0228       //       std::cout << "index offset, current hit count = " << crossingSimHitIndexOffset_[tag.encode()] << ", " << simHits->size() << std::endl;
0229       if (simHits.isValid())
0230         crossingSimHitIndexOffset_[tag.encode()] += simHits->size();
0231     }
0232   }
0233
0234   void SiPixelDigitizer::accumulate(PileUpEventPrincipal const& iEvent,
0235                                     edm::EventSetup const& iSetup,
0236                                     edm::StreamID const& streamID) {
0237     // Step A: Get Inputs
0238     for (vstring::const_iterator i = trackerContainers.begin(), iEnd = trackerContainers.end(); i != iEnd; ++i) {
0239       edm::Handle<std::vector<PSimHit> > simHits;
0240       edm::InputTag tag(hitsProducer, *i);
0241
0242       iEvent.getByLabel(tag, simHits);
0243       unsigned int tofBin = PixelDigiSimLink::LowTof;
0244       if ((*i).find(std::string("HighTof")) != std::string::npos)
0245         tofBin = PixelDigiSimLink::HighTof;
0246       accumulatePixelHits(simHits, crossingSimHitIndexOffset_[tag.encode()], tofBin, iSetup);
0247       // Now that the hits have been processed, I'll add the amount of hits in this crossing on to
0248       // the global counter. Next time accumulateStripHits() is called it will count the sim hits
0249       // as though they were on the end of this collection.
0250       // Note that this is only used for creating digi-sim links (if configured to do so).
0251       //       std::cout << "index offset, current hit count = " << crossingSimHitIndexOffset_[tag.encode()] << ", " << simHits->size() << std::endl;
0252       if (simHits.isValid())
0253         crossingSimHitIndexOffset_[tag.encode()] += simHits->size();
0254     }
0255   }
0256
0257   // ------------ method called to produce the data  ------------
0258   void SiPixelDigitizer::finalizeEvent(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0259     const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0260
0261     std::vector<edm::DetSet<PixelDigi> > theDigiVector;
0262     std::vector<edm::DetSet<PixelDigiSimLink> > theDigiLinkVector;
0263     std::vector<edm::DetSet<PixelSimHitExtraInfo> > theExtraSimHitInfoVector;
0264
0265     if (firstFinalizeEvent_) {
0266       _pixeldigialgo->init_DynIneffDB(iSetup);
0267       firstFinalizeEvent_ = false;
0268     }
0269     _pixeldigialgo->calculateInstlumiFactor(PileupInfo_.get());
0270
0271     if (_pixeldigialgo->killBadFEDChannels()) {
0272       std::unique_ptr<PixelFEDChannelCollection> PixelFEDChannelCollection_ =
0273           _pixeldigialgo->chooseScenario(PileupInfo_.get(), randomEngine_);
0274       if (PixelFEDChannelCollection_ == nullptr) {
0275         throw cms::Exception("NullPointerError") << "PixelFEDChannelCollection not set in chooseScenario function.\n";
0276       }
0277       iEvent.put(std::move(PixelFEDChannelCollection_));
0278     }
0279
0280     for (const auto& iu : pDD->detUnits()) {
0281       if (iu->type().isTrackerPixel()) {
0282         //
0283
0284         edm::DetSet<PixelDigi> collector(iu->geographicalId().rawId());
0285         edm::DetSet<PixelDigiSimLink> linkcollector(iu->geographicalId().rawId());
0286         std::vector<PixelDigiAddTempInfo> tempcollector;
0287         edm::DetSet<PixelSimHitExtraInfo> tempSHcollector(iu->geographicalId().rawId());
0288
0289         _pixeldigialgo->digitize(dynamic_cast<const PixelGeomDetUnit*>(iu),
0290                                  collector.data,
0291                                  linkcollector.data,
0292                                  tempcollector,
0293                                  tTopo,
0294                                  randomEngine_);
0295
0296         // transformation of the tempcollector  (October 15, 2021)
0297         if (!tempcollector.empty()) {
0298           std::vector<PixelDigiAddTempInfo>::const_iterator loopNewClass;
0299           unsigned int channelPrevious2 = -1;
0300           size_t hitFirstOne2 = -1;
0301           for (loopNewClass = tempcollector.begin(); loopNewClass != tempcollector.end(); ++loopNewClass) {
0302             //   check if the new SimHit already exists in the class
0303             //      if yes : add only the Digi info to the existing entry
0304             //      if not : create a new entry
0305             //          PixelSimHitExtraInfo(
0306             //          size_t Hindex, Local3DPoint entryP , Local3DPoint exitP, uint32_t detID, std::vector<unsigned int> ch, std::vector<float> InitCharge)
0307             //   what about the duplicates : 2 SimHits associated to the same Digis ?
0308             //
0309             //
0310
0311             bool checkTwoSimHits = false;
0312             if (channelPrevious2 == loopNewClass->channel() && hitFirstOne2 != loopNewClass->hitIndex()) {
0313               // case of one Digi associated to a second SimHit
0314               checkTwoSimHits = true;
0315             } else {
0316               channelPrevious2 = loopNewClass->channel();
0317               hitFirstOne2 = loopNewClass->hitIndex();
0318             }
0319
0320             bool checkInTheList = false;
0321             if (!checkTwoSimHits) {
0322               std::vector<PixelSimHitExtraInfo>::iterator loopTempSH;
0323               for (loopTempSH = tempSHcollector.begin(); loopTempSH != tempSHcollector.end(); ++loopTempSH) {
0324                 if (loopNewClass->hitIndex() == loopTempSH->hitIndex()) {
0325                   checkInTheList = true;
0326                   loopTempSH->addDigiInfo(loopNewClass->channel());
0327                 }
0328               }
0329               if (!checkInTheList) {
0330                 PixelSimHitExtraInfo newSHEntry(loopNewClass->hitIndex(),
0331                                                 loopNewClass->entryPoint(),
0332                                                 loopNewClass->exitPoint(),
0333                                                 loopNewClass->channel());
0334                 tempSHcollector.push_back(newSHEntry);
0335               }
0336             }
0337           }
0338         }
0339
0340         if (applyLateReweighting_) {
0341           // if applyLateReweighting_  is true, the charge reweighting has to be applied on top of the digis
0342           _pixeldigialgo->lateSignalReweight(
0343               dynamic_cast<const PixelGeomDetUnit*>(iu), collector.data, tempSHcollector.data, tTopo, randomEngine_);
0344         }
0345
0346         if (!collector.data.empty()) {
0347           theDigiVector.push_back(std::move(collector));
0348         }
0349         if (!linkcollector.data.empty()) {
0350           theDigiLinkVector.push_back(std::move(linkcollector));
0351         }
0352         if (!tempSHcollector.data.empty()) {
0353           theExtraSimHitInfoVector.push_back(std::move(tempSHcollector));
0354         }
0355       }
0356     }
0357     _pixeldigialgo->resetSimHitMaps();
0358
0359     // Step C: create collection with the cache vector of DetSet
0360     std::unique_ptr<edm::DetSetVector<PixelDigi> > output(new edm::DetSetVector<PixelDigi>(theDigiVector));
0361     std::unique_ptr<edm::DetSetVector<PixelDigiSimLink> > outputlink(
0362         new edm::DetSetVector<PixelDigiSimLink>(theDigiLinkVector));
0363     std::unique_ptr<edm::DetSetVector<PixelSimHitExtraInfo> > outputExtraSim(
0364         new edm::DetSetVector<PixelSimHitExtraInfo>(theExtraSimHitInfoVector));
0365
0366     // Step D: write output to file
0367     iEvent.put(std::move(output));
0368     iEvent.put(std::move(outputlink));
0369     if (store_SimHitEntryExitPoints_)
0370       iEvent.put(std::move(outputExtraSim));
0371
0372     randomEngine_ = nullptr;  // to prevent access outside event
0373   }
0374 }  // namespace cms