DataMixingModule/plugins/DataMixingHcalWorker.cc

0001 // File: DataMixingHcalWorker.cc
0002 // Description:  see DataMixingHcalWorker.h
0003 // Author:  Mike Hildreth, University of Notre Dame
0004 //
0005 //--------------------------------------------
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include <map>
0010 #include <memory>
0011 //
0012 //
0013 #include "DataMixingHcalWorker.h"
0014
0015 using namespace std;
0016
0017 namespace edm {
0018
0019   // Virtual constructor
0020
0021   DataMixingHcalWorker::DataMixingHcalWorker() {}
0022
0023   // Constructor
0024   DataMixingHcalWorker::DataMixingHcalWorker(const edm::ParameterSet &ps, edm::ConsumesCollector &&iC)
0025       : label_(ps.getParameter<std::string>("Label"))
0026
0027   {
0028     // get the subdetector names
0029     //    this->getSubdetectorNames();  //something like this may be useful to
0030     //    check what we are supposed to do...
0031
0032     // declare the products to produce
0033
0034     // Hcal
0035
0036     HBHErechitCollectionSig_ = ps.getParameter<edm::InputTag>("HBHEProducerSig");
0037     HOrechitCollectionSig_ = ps.getParameter<edm::InputTag>("HOProducerSig");
0038     HFrechitCollectionSig_ = ps.getParameter<edm::InputTag>("HFProducerSig");
0039     ZDCrechitCollectionSig_ = ps.getParameter<edm::InputTag>("ZDCrechitCollectionSig");
0040
0041     HBHEPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HBHEPileRecHitInputTag");
0042     HOPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HOPileRecHitInputTag");
0043     HFPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HFPileRecHitInputTag");
0044     ZDCPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("ZDCPileRecHitInputTag");
0045
0046     HBHERecHitCollectionDM_ = ps.getParameter<std::string>("HBHERecHitCollectionDM");
0047     HORecHitCollectionDM_ = ps.getParameter<std::string>("HORecHitCollectionDM");
0048     HFRecHitCollectionDM_ = ps.getParameter<std::string>("HFRecHitCollectionDM");
0049     ZDCRecHitCollectionDM_ = ps.getParameter<std::string>("ZDCRecHitCollectionDM");
0050
0051     HBHERecHitToken_ = iC.consumes<HBHERecHitCollection>(HBHErechitCollectionSig_);
0052     HORecHitToken_ = iC.consumes<HORecHitCollection>(HOrechitCollectionSig_);
0053     HFRecHitToken_ = iC.consumes<HFRecHitCollection>(HFrechitCollectionSig_);
0054
0055     HBHERecHitPToken_ = iC.consumes<HBHERecHitCollection>(HBHEPileRecHitInputTag_);
0056     HORecHitPToken_ = iC.consumes<HORecHitCollection>(HOPileRecHitInputTag_);
0057     HFRecHitPToken_ = iC.consumes<HFRecHitCollection>(HFPileRecHitInputTag_);
0058
0059     ZDCRecHitToken_ = iC.consumes<ZDCRecHitCollection>(ZDCrechitCollectionSig_);
0060     ZDCRecHitPToken_ = iC.consumes<ZDCRecHitCollection>(ZDCPileRecHitInputTag_);
0061   }
0062
0063   // Virtual destructor needed.
0064   DataMixingHcalWorker::~DataMixingHcalWorker() {}
0065
0066   void DataMixingHcalWorker::addHcalSignals(const edm::Event &e) {
0067     // fill in maps of hits
0068
0069     LogInfo("DataMixingHcalWorker") << "===============> adding MC signals for " << e.id();
0070
0071     // HBHE first
0072
0073     Handle<HBHERecHitCollection> pHBHERecHits;
0074
0075     const HBHERecHitCollection *HBHERecHits = nullptr;
0076
0077     if (e.getByToken(HBHERecHitToken_, pHBHERecHits)) {
0078       HBHERecHits = pHBHERecHits.product();  // get a ptr to the product
0079       LogDebug("DataMixingHcalWorker") << "total # HBHE rechits: " << HBHERecHits->size();
0080     }
0081
0082     if (HBHERecHits) {
0083       // loop over rechits, storing them in a map so we can add pileup later
0084       for (HBHERecHitCollection::const_iterator it = HBHERecHits->begin(); it != HBHERecHits->end(); ++it) {
0085         HBHERecHitStorage_.insert(HBHERecHitMap::value_type((it->id()), *it));
0086
0087 #ifdef DEBUG
0088         LogDebug("DataMixingHcalWorker") << "processed HBHERecHit with rawId: " << it->id() << "\n"
0089                                          << " rechit energy: " << it->energy();
0090 #endif
0091       }
0092     }
0093
0094     // HO next
0095
0096     Handle<HORecHitCollection> pHORecHits;
0097
0098     const HORecHitCollection *HORecHits = nullptr;
0099
0100     if (e.getByToken(HORecHitToken_, pHORecHits)) {
0101       HORecHits = pHORecHits.product();  // get a ptr to the product
0102 #ifdef DEBUG
0103       LogDebug("DataMixingHcalWorker") << "total # HO rechits: " << HORecHits->size();
0104 #endif
0105     }
0106
0107     if (HORecHits) {
0108       // loop over rechits, storing them in a map so we can add pileup later
0109       for (HORecHitCollection::const_iterator it = HORecHits->begin(); it != HORecHits->end(); ++it) {
0110         HORecHitStorage_.insert(HORecHitMap::value_type((it->id()), *it));
0111
0112 #ifdef DEBUG
0113         LogDebug("DataMixingHcalWorker") << "processed HORecHit with rawId: " << it->id() << "\n"
0114                                          << " rechit energy: " << it->energy();
0115 #endif
0116       }
0117     }
0118
0119     // HF next
0120
0121     Handle<HFRecHitCollection> pHFRecHits;
0122
0123     const HFRecHitCollection *HFRecHits = nullptr;
0124
0125     if (e.getByToken(HFRecHitToken_, pHFRecHits)) {
0126       HFRecHits = pHFRecHits.product();  // get a ptr to the product
0127 #ifdef DEBUG
0128       LogDebug("DataMixingHcalWorker") << "total # HF rechits: " << HFRecHits->size();
0129 #endif
0130     }
0131
0132     if (HFRecHits) {
0133       // loop over rechits, storing them in a map so we can add pileup later
0134       for (HFRecHitCollection::const_iterator it = HFRecHits->begin(); it != HFRecHits->end(); ++it) {
0135         HFRecHitStorage_.insert(HFRecHitMap::value_type((it->id()), *it));
0136
0137 #ifdef DEBUG
0138         LogDebug("DataMixingHcalWorker") << "processed HFRecHit with rawId: " << it->id() << "\n"
0139                                          << " rechit energy: " << it->energy();
0140 #endif
0141       }
0142     }
0143
0144     // ZDC next
0145
0146     Handle<ZDCRecHitCollection> pZDCRecHits;
0147
0148     const ZDCRecHitCollection *ZDCRecHits = nullptr;
0149
0150     if (e.getByToken(ZDCRecHitToken_, pZDCRecHits)) {
0151       ZDCRecHits = pZDCRecHits.product();  // get a ptr to the product
0152 #ifdef DEBUG
0153       LogDebug("DataMixingHcalWorker") << "total # ZDC rechits: " << ZDCRecHits->size();
0154 #endif
0155     }
0156
0157     if (ZDCRecHits) {
0158       // loop over rechits, storing them in a map so we can add pileup later
0159       for (ZDCRecHitCollection::const_iterator it = ZDCRecHits->begin(); it != ZDCRecHits->end(); ++it) {
0160         ZDCRecHitStorage_.insert(ZDCRecHitMap::value_type((it->id()), *it));
0161
0162 #ifdef DEBUG
0163         LogDebug("DataMixingHcalWorker") << "processed ZDCRecHit with rawId: " << it->id() << "\n"
0164                                          << " rechit energy: " << it->energy();
0165 #endif
0166       }
0167     }
0168
0169   }  // end of addEMSignals
0170
0171   void DataMixingHcalWorker::addHcalPileups(const int bcr,
0172                                             const EventPrincipal *ep,
0173                                             unsigned int eventNr,
0174                                             ModuleCallingContext const *mcc) {
0175     LogDebug("DataMixingHcalWorker") << "\n===============> adding pileups from event  " << ep->id()
0176                                      << " for bunchcrossing " << bcr;
0177
0178     // fill in maps of hits; same code as addSignals, except now applied to the
0179     // pileup events
0180
0181     // HBHE first
0182
0183     std::shared_ptr<Wrapper<HBHERecHitCollection> const> HBHERecHitsPTR =
0184         getProductByTag<HBHERecHitCollection>(*ep, HBHEPileRecHitInputTag_, mcc);
0185
0186     if (HBHERecHitsPTR) {
0187       const HBHERecHitCollection *HBHERecHits = const_cast<HBHERecHitCollection *>(HBHERecHitsPTR->product());
0188
0189       LogDebug("DataMixingEMWorker") << "total # HBHE rechits: " << HBHERecHits->size();
0190
0191       // loop over rechits, adding these to the existing maps
0192       for (HBHERecHitCollection::const_iterator it = HBHERecHits->begin(); it != HBHERecHits->end(); ++it) {
0193         HBHERecHitStorage_.insert(HBHERecHitMap::value_type((it->id()), *it));
0194
0195 #ifdef DEBUG
0196         LogDebug("DataMixingEMWorker") << "processed HBHERecHit with rawId: " << it->id().rawId() << "\n"
0197                                        << " rechit energy: " << it->energy();
0198 #endif
0199       }
0200     }
0201
0202     // HO Next
0203
0204     std::shared_ptr<Wrapper<HORecHitCollection> const> HORecHitsPTR =
0205         getProductByTag<HORecHitCollection>(*ep, HOPileRecHitInputTag_, mcc);
0206
0207     if (HORecHitsPTR) {
0208       const HORecHitCollection *HORecHits = const_cast<HORecHitCollection *>(HORecHitsPTR->product());
0209
0210       LogDebug("DataMixingEMWorker") << "total # HO rechits: " << HORecHits->size();
0211
0212       // loop over rechits, adding these to the existing maps
0213       for (HORecHitCollection::const_iterator it = HORecHits->begin(); it != HORecHits->end(); ++it) {
0214         HORecHitStorage_.insert(HORecHitMap::value_type((it->id()), *it));
0215
0216 #ifdef DEBUG
0217         LogDebug("DataMixingEMWorker") << "processed HORecHit with rawId: " << it->id().rawId() << "\n"
0218                                        << " rechit energy: " << it->energy();
0219 #endif
0220       }
0221     }
0222
0223     // HF Next
0224
0225     std::shared_ptr<Wrapper<HFRecHitCollection> const> HFRecHitsPTR =
0226         getProductByTag<HFRecHitCollection>(*ep, HFPileRecHitInputTag_, mcc);
0227
0228     if (HFRecHitsPTR) {
0229       const HFRecHitCollection *HFRecHits = const_cast<HFRecHitCollection *>(HFRecHitsPTR->product());
0230
0231       LogDebug("DataMixingEMWorker") << "total # HF rechits: " << HFRecHits->size();
0232
0233       // loop over rechits, adding these to the existing maps
0234       for (HFRecHitCollection::const_iterator it = HFRecHits->begin(); it != HFRecHits->end(); ++it) {
0235         HFRecHitStorage_.insert(HFRecHitMap::value_type((it->id()), *it));
0236
0237 #ifdef DEBUG
0238         LogDebug("DataMixingEMWorker") << "processed HFRecHit with rawId: " << it->id().rawId() << "\n"
0239                                        << " rechit energy: " << it->energy();
0240 #endif
0241       }
0242     }
0243
0244     // ZDC Next
0245
0246     std::shared_ptr<Wrapper<ZDCRecHitCollection> const> ZDCRecHitsPTR =
0247         getProductByTag<ZDCRecHitCollection>(*ep, ZDCPileRecHitInputTag_, mcc);
0248
0249     if (ZDCRecHitsPTR) {
0250       const ZDCRecHitCollection *ZDCRecHits = const_cast<ZDCRecHitCollection *>(ZDCRecHitsPTR->product());
0251
0252       LogDebug("DataMixingEMWorker") << "total # ZDC rechits: " << ZDCRecHits->size();
0253
0254       // loop over rechits, adding these to the existing maps
0255       for (ZDCRecHitCollection::const_iterator it = ZDCRecHits->begin(); it != ZDCRecHits->end(); ++it) {
0256         ZDCRecHitStorage_.insert(ZDCRecHitMap::value_type((it->id()), *it));
0257
0258 #ifdef DEBUG
0259         LogDebug("DataMixingEMWorker") << "processed ZDCRecHit with rawId: " << it->id().rawId() << "\n"
0260                                        << " rechit energy: " << it->energy();
0261 #endif
0262       }
0263     }
0264   }
0265
0266   void DataMixingHcalWorker::putHcal(edm::Event &e) {
0267     // collection of rechits to put in the event
0268     std::unique_ptr<HBHERecHitCollection> HBHErechits(new HBHERecHitCollection);
0269     std::unique_ptr<HORecHitCollection> HOrechits(new HORecHitCollection);
0270     std::unique_ptr<HFRecHitCollection> HFrechits(new HFRecHitCollection);
0271     std::unique_ptr<ZDCRecHitCollection> ZDCrechits(new ZDCRecHitCollection);
0272
0273     // loop over the maps we have, re-making individual hits or rechits if
0274     // necessary.
0275     DetId formerID = 0;
0276     DetId currentID;
0277     float ESum = 0.;
0278     float HBTime = 0.;
0279
0280     // HB first...
0281
0282     HBHERecHitMap::const_iterator iHBchk;
0283
0284     for (HBHERecHitMap::const_iterator iHB = HBHERecHitStorage_.begin(); iHB != HBHERecHitStorage_.end(); ++iHB) {
0285       currentID = iHB->first;
0286
0287       if (currentID == formerID) {  // we have to add these rechits together
0288
0289         ESum += (iHB->second).energy();
0290
0291       } else {
0292         if (formerID > 0) {
0293           // cutoff for ESum?
0294           HBHERecHit aHit(formerID, ESum, HBTime);
0295           HBHErechits->push_back(aHit);
0296         }
0297         // save pointers for next iteration
0298         formerID = currentID;
0299         ESum = (iHB->second).energy();
0300         HBTime = (iHB->second).time();  // take time of first hit in sequence - is this ok?
0301       }
0302
0303       iHBchk = iHB;
0304       if ((++iHBchk) == HBHERecHitStorage_.end()) {  // make sure not to lose the last one
0305         HBHERecHit aHit(formerID, ESum, HBTime);
0306         HBHErechits->push_back(aHit);
0307       }
0308     }
0309
0310     // HO next...
0311
0312     // loop over the maps we have, re-making individual hits or rechits if
0313     // necessary.
0314     formerID = 0;
0315     ESum = 0.;
0316     float HOTime = 0.;
0317
0318     HORecHitMap::const_iterator iHOchk;
0319
0320     for (HORecHitMap::const_iterator iHO = HORecHitStorage_.begin(); iHO != HORecHitStorage_.end(); ++iHO) {
0321       currentID = iHO->first;
0322
0323       if (currentID == formerID) {  // we have to add these rechits together
0324
0325         ESum += (iHO->second).energy();
0326
0327       } else {
0328         if (formerID > 0) {
0329           // cutoff for ESum?
0330           HORecHit aHit(formerID, ESum, HOTime);
0331           HOrechits->push_back(aHit);
0332         }
0333         // save pointers for next iteration
0334         formerID = currentID;
0335         ESum = (iHO->second).energy();
0336         HOTime = (iHO->second).time();  // take time of first hit in sequence - is this ok?
0337       }
0338
0339       iHOchk = iHO;
0340       if ((++iHOchk) == HORecHitStorage_.end()) {  // make sure not to lose the last one
0341         HORecHit aHit(formerID, ESum, HOTime);
0342         HOrechits->push_back(aHit);
0343       }
0344     }
0345
0346     // HF next...
0347
0348     // loop over the maps we have, re-making individual hits or rechits if
0349     // necessary.
0350     formerID = 0;
0351     ESum = 0.;
0352     float HFTime = 0.;
0353
0354     HFRecHitMap::const_iterator iHFchk;
0355
0356     for (HFRecHitMap::const_iterator iHF = HFRecHitStorage_.begin(); iHF != HFRecHitStorage_.end(); ++iHF) {
0357       currentID = iHF->first;
0358
0359       if (currentID == formerID) {  // we have to add these rechits together
0360
0361         ESum += (iHF->second).energy();
0362
0363       } else {
0364         if (formerID > 0) {
0365           // cutoff for ESum?
0366           HFRecHit aHit(formerID, ESum, HFTime);
0367           HFrechits->push_back(aHit);
0368         }
0369         // save pointers for next iteration
0370         formerID = currentID;
0371         ESum = (iHF->second).energy();
0372         HFTime = (iHF->second).time();  // take time of first hit in sequence - is this ok?
0373       }
0374
0375       iHFchk = iHF;
0376       if ((++iHFchk) == HFRecHitStorage_.end()) {  // make sure not to lose the last one
0377         HFRecHit aHit(formerID, ESum, HBTime);
0378         HFrechits->push_back(aHit);
0379       }
0380     }
0381
0382     // ZDC next...
0383
0384     // loop over the maps we have, re-making individual hits or rechits if
0385     // necessary.
0386     formerID = 0;
0387     ESum = 0.;
0388     float ZDCTime = 0.;
0389     float lowGainEnergy = 0;
0390     ZDCRecHit ZOldHit;
0391
0392     ZDCRecHitMap::const_iterator iZDCchk;
0393
0394     for (ZDCRecHitMap::const_iterator iZDC = ZDCRecHitStorage_.begin(); iZDC != ZDCRecHitStorage_.end(); ++iZDC) {
0395       currentID = iZDC->first;
0396
0397       if (currentID == formerID) {  // we have to add these rechits together
0398
0399         ESum += (iZDC->second).energy();
0400
0401       } else {
0402         if (formerID > 0) {
0403           // cutoff for ESum?
0404           ZDCRecHit aHit(formerID, ESum, ZDCTime, lowGainEnergy);
0405           ZDCrechits->push_back(aHit);
0406         }
0407         // save pointers for next iteration
0408         formerID = currentID;
0409         ESum = (iZDC->second).energy();
0410         lowGainEnergy = (iZDC->second).lowGainEnergy();
0411         ZDCTime = (iZDC->second).time();  // take time of first hit in sequence - is this ok?
0412       }
0413
0414       iZDCchk = iZDC;
0415       if ((++iZDCchk) == ZDCRecHitStorage_.end()) {  // make sure not to lose the last one
0416         ZDCRecHit aHit(formerID, ESum, HBTime, lowGainEnergy);
0417         ZDCrechits->push_back(aHit);
0418       }
0419     }
0420
0421     // done merging
0422
0423     // put the collection of recunstructed hits in the event
0424     LogInfo("DataMixingHcalWorker") << "total # HBHE Merged rechits: " << HBHErechits->size();
0425     LogInfo("DataMixingHcalWorker") << "total # HO Merged rechits: " << HOrechits->size();
0426     LogInfo("DataMixingHcalWorker") << "total # HF Merged rechits: " << HFrechits->size();
0427     LogInfo("DataMixingHcalWorker") << "total # ZDC Merged rechits: " << ZDCrechits->size();
0428
0429     e.put(std::move(HBHErechits), HBHERecHitCollectionDM_);
0430     e.put(std::move(HOrechits), HORecHitCollectionDM_);
0431     e.put(std::move(HFrechits), HFRecHitCollectionDM_);
0432     e.put(std::move(ZDCrechits), ZDCRecHitCollectionDM_);
0433
0434     // clear local storage after this event
0435     HBHERecHitStorage_.clear();
0436     HORecHitStorage_.clear();
0437     HFRecHitStorage_.clear();
0438     ZDCRecHitStorage_.clear();
0439   }
0440
0441 }  // namespace edm