CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/SimCalorimetry/HcalSimProducers/src/HcalDigitizer.cc

Go to the documentation of this file.
00001 #include "SimCalorimetry/HcalSimProducers/interface/HcalDigitizer.h"
00002 #include "SimCalorimetry/HcalSimProducers/src/HcalTestHitGenerator.h"
00003 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
00004 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
00005 #include "SimCalorimetry/HcalSimAlgos/interface/HcalShapes.h"
00006 #include "SimCalorimetry/HcalSimAlgos/interface/HcalElectronicsSim.h"
00007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
00008 #include "SimCalorimetry/HcalSimAlgos/interface/HcalAmplifier.h"
00009 #include "SimCalorimetry/HcalSimAlgos/interface/HcalCoderFactory.h"
00010 #include "SimCalorimetry/HcalSimAlgos/interface/HcalHitCorrection.h"
00011 #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
00012 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
00013 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h"
00014 #include "SimCalorimetry/HcalSimAlgos/interface/HPDIonFeedbackSim.h"
00015 #include "DataFormats/Common/interface/Handle.h"
00016 #include "DataFormats/Common/interface/Wrapper.h"
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include "FWCore/Framework/interface/Event.h"
00019 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
00020 #include "FWCore/Framework/interface/EventSetup.h"
00021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 #include "SimCalorimetry/CaloSimAlgos/interface/CaloTDigitizer.h"
00024 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00025 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00026 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00027 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00028 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
00029 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
00030 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00031 #include "FWCore/ServiceRegistry/interface/Service.h"
00032 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
00033 #include "SimCalorimetry/HcalSimAlgos/interface/HPDNoiseGenerator.h"
00034 #include "CondFormats/HcalObjects/interface/HcalCholeskyMatrix.h"
00035 #include "CondFormats/HcalObjects/interface/HcalCholeskyMatrices.h"
00036 #include <boost/foreach.hpp>
00037 #include "Geometry/CaloTopology/interface/HcalTopology.h"
00038 
00039 namespace HcalDigitizerImpl {
00040 
00041   template<typename SIPMDIGITIZER>
00042   void fillSiPMCells(const std::vector<int> & siPMCells, SIPMDIGITIZER * siPMDigitizer)
00043   {
00044     std::vector<DetId> siPMDetIds;
00045     siPMDetIds.reserve(siPMCells.size());
00046     for(std::vector<int>::const_iterator idItr = siPMCells.begin();
00047         idItr != siPMCells.end(); ++idItr)
00048     {
00049       siPMDetIds.emplace_back(*idItr);
00050     }
00051     siPMDigitizer->setDetIds(siPMDetIds);
00052   }
00053 
00054   // if both exist, assume the SiPM one has cells filled, and
00055   // assign the rest to the HPD
00056   template<typename HPDDIGITIZER, typename SIPMDIGITIZER>
00057   void fillCells(const std::vector<DetId>& allCells,
00058                  HPDDIGITIZER * hpdDigitizer,
00059                  SIPMDIGITIZER * siPMDigitizer)
00060   {
00061     // if both digitizers exist, split up the cells
00062     if(siPMDigitizer && hpdDigitizer)
00063     {
00064       std::vector<DetId> siPMDetIds = siPMDigitizer->detIds();
00065       std::sort(siPMDetIds.begin(), siPMDetIds.end());
00066       std::vector<DetId> sortedCells = allCells;
00067       std::sort(sortedCells.begin(), sortedCells.end());
00068       std::vector<DetId> hpdCells;
00069       std::set_difference(sortedCells.begin(), sortedCells.end(),
00070                           siPMDetIds.begin(), siPMDetIds.end(),
00071                           std::back_inserter(hpdCells) );
00072       hpdDigitizer->setDetIds(hpdCells);
00073     }
00074     else
00075     {
00076       if(siPMDigitizer) siPMDigitizer->setDetIds(allCells);
00077       if(hpdDigitizer) hpdDigitizer->setDetIds(allCells);
00078     }
00079   }
00080 } // namespace HcaiDigitizerImpl
00081 
00082 
00083 HcalDigitizer::HcalDigitizer(const edm::ParameterSet& ps) 
00084 : theGeometry(0),
00085   theParameterMap(new HcalSimParameterMap(ps)),
00086   theShapes(new HcalShapes()),
00087   theHBHEResponse(0),
00088   theHBHESiPMResponse(0),
00089   theHOResponse(0),   
00090   theHOSiPMResponse(0),
00091   theHFResponse(new CaloHitResponse(theParameterMap, theShapes)),
00092   theZDCResponse(new CaloHitResponse(theParameterMap, theShapes)),
00093   theHBHEAmplifier(0),
00094   theHFAmplifier(0),
00095   theHOAmplifier(0),
00096   theZDCAmplifier(0),
00097   theIonFeedback(0),
00098   theCoderFactory(0),
00099   theHBHEElectronicsSim(0),
00100   theHFElectronicsSim(0),
00101   theHOElectronicsSim(0),
00102   theZDCElectronicsSim(0),
00103   theHBHEHitFilter(),
00104   theHFHitFilter(ps.getParameter<bool>("doHFWindow")),
00105   theHOHitFilter(),
00106   theHOSiPMHitFilter(HcalOuter),
00107   theZDCHitFilter(),
00108   theHitCorrection(0),
00109   theNoiseGenerator(0),
00110   theNoiseHitGenerator(0),
00111   theHBHEDigitizer(0),
00112   theHBHESiPMDigitizer(0),
00113   theHODigitizer(0),
00114   theHOSiPMDigitizer(0),
00115   theHFDigitizer(0),
00116   theZDCDigitizer(0),
00117   theHBHEDetIds(),
00118   theHOHPDDetIds(),
00119   theHOSiPMDetIds(),
00120   isZDC(true),
00121   isHCAL(true),
00122   zdcgeo(true),
00123   hbhegeo(true),
00124   hogeo(true),
00125   hfgeo(true),
00126   theHOSiPMCode(ps.getParameter<edm::ParameterSet>("ho").getParameter<int>("siPMCode"))
00127 {
00128   bool doNoise = ps.getParameter<bool>("doNoise");
00129   bool useOldNoiseHB = ps.getParameter<bool>("useOldHB");
00130   bool useOldNoiseHE = ps.getParameter<bool>("useOldHE");
00131   bool useOldNoiseHF = ps.getParameter<bool>("useOldHF");
00132   bool useOldNoiseHO = ps.getParameter<bool>("useOldHO");
00133   bool doEmpty = ps.getParameter<bool>("doEmpty");
00134   double HBtp = ps.getParameter<double>("HBTuningParameter");
00135   double HEtp = ps.getParameter<double>("HETuningParameter");
00136   double HFtp = ps.getParameter<double>("HFTuningParameter");
00137   double HOtp = ps.getParameter<double>("HOTuningParameter");
00138 
00139   // need to make copies, because they might get different noise generators
00140   theHBHEAmplifier = new HcalAmplifier(theParameterMap, doNoise);
00141   theHFAmplifier = new HcalAmplifier(theParameterMap, doNoise);
00142   theHOAmplifier = new HcalAmplifier(theParameterMap, doNoise);
00143   theZDCAmplifier = new HcalAmplifier(theParameterMap, doNoise);
00144   theHBHEAmplifier->setHBtuningParameter(HBtp);
00145   theHBHEAmplifier->setHEtuningParameter(HEtp);
00146   theHFAmplifier->setHFtuningParameter(HFtp);
00147   theHOAmplifier->setHOtuningParameter(HOtp);
00148   theHBHEAmplifier->setUseOldHB(useOldNoiseHB);
00149   theHBHEAmplifier->setUseOldHE(useOldNoiseHE);
00150   theHFAmplifier->setUseOldHF(useOldNoiseHF);
00151   theHOAmplifier->setUseOldHO(useOldNoiseHO);
00152 
00153   theCoderFactory = new HcalCoderFactory(HcalCoderFactory::DB);
00154   theHBHEElectronicsSim = new HcalElectronicsSim(theHBHEAmplifier, theCoderFactory);
00155   theHFElectronicsSim = new HcalElectronicsSim(theHFAmplifier, theCoderFactory);
00156   theHOElectronicsSim = new HcalElectronicsSim(theHOAmplifier, theCoderFactory);
00157   theZDCElectronicsSim = new HcalElectronicsSim(theZDCAmplifier, theCoderFactory);
00158 
00159   // a code of 1 means make all cells SiPM
00160   std::vector<int> hbSiPMCells(ps.getParameter<edm::ParameterSet>("hb").getParameter<std::vector<int> >("siPMCells"));
00161   //std::vector<int> hoSiPMCells(ps.getParameter<edm::ParameterSet>("ho").getParameter<std::vector<int> >("siPMCells"));
00162   // 0 means none, 1 means all, and 2 means use hardcoded
00163 
00164   bool doHBHEHPD = hbSiPMCells.empty() || (hbSiPMCells[0] != 1);
00165   bool doHOHPD = (theHOSiPMCode != 1);
00166   bool doHBHESiPM = !hbSiPMCells.empty();
00167   bool doHOSiPM = (theHOSiPMCode != 0);
00168   if(doHBHEHPD)
00169   {
00170     theHBHEResponse = new CaloHitResponse(theParameterMap, theShapes);
00171     edm::LogInfo("HcalDigitizer") <<"Set scale for HB towers";
00172     theHBHEResponse->initHBHEScale(); //GMA
00173 
00174     theHBHEResponse->setHitFilter(&theHBHEHitFilter);
00175     theHBHEDigitizer = new HBHEDigitizer(theHBHEResponse, theHBHEElectronicsSim, doEmpty);
00176     bool    changeResponse = ps.getParameter<bool>("ChangeResponse");
00177     edm::FileInPath fname  = ps.getParameter<edm::FileInPath>("CorrFactorFile");
00178     if (changeResponse) {
00179       std::string corrFileName = fname.fullPath();
00180       edm::LogInfo("HcalDigitizer") << "Set scale for HB towers from " << corrFileName;
00181       theHBHEResponse->setHBHEScale(corrFileName); //GMA
00182     }
00183   }
00184   if(doHOHPD) 
00185   {
00186     theHOResponse = new CaloHitResponse(theParameterMap, theShapes);
00187     theHOResponse->setHitFilter(&theHOHitFilter);
00188     theHODigitizer = new HODigitizer(theHOResponse, theHOElectronicsSim, doEmpty);
00189   }
00190 
00191   if(doHBHESiPM)
00192   {
00193     theHBHESiPMResponse = new HcalSiPMHitResponse(theParameterMap, theShapes);
00194     theHBHESiPMResponse->setHitFilter(&theHBHEHitFilter);
00195     theHBHESiPMDigitizer = new HBHEDigitizer(theHBHESiPMResponse, theHBHEElectronicsSim, doEmpty);
00196   }
00197   if(doHOSiPM)
00198   {
00199     theHOSiPMResponse = new HcalSiPMHitResponse(theParameterMap, theShapes);
00200     theHOSiPMResponse->setHitFilter(&theHOSiPMHitFilter);
00201     theHOSiPMDigitizer = new HODigitizer(theHOSiPMResponse, theHOElectronicsSim, doEmpty);
00202   }
00203 
00204   // if both are present, fill the SiPM cells now
00205   if(doHBHEHPD && doHBHESiPM)
00206   {
00207     HcalDigitizerImpl::fillSiPMCells(hbSiPMCells, theHBHESiPMDigitizer);
00208   }
00209 
00210   theHFResponse->setHitFilter(&theHFHitFilter);
00211   theZDCResponse->setHitFilter(&theZDCHitFilter);
00212 
00213   bool doTimeSlew = ps.getParameter<bool>("doTimeSlew");
00214   if(doTimeSlew) {
00215     // no time slewing for HF
00216     /*  
00217     theHitCorrection = new HcalHitCorrection(theParameterMap);
00218     if(theHBHEResponse) theHBHEResponse->setHitCorrection(theHitCorrection);
00219     if(theHBHESiPMResponse) theHBHESiPMResponse->setHitCorrection(theHitCorrection);
00220     if(theHOResponse) theHOResponse->setHitCorrection(theHitCorrection);
00221     if(theHOSiPMResponse) theHOSiPMResponse->setHitCorrection(theHitCorrection);
00222     theZDCResponse->setHitCorrection(theHitCorrection);
00223     */ 
00224     theTimeSlewSim = new HcalTimeSlewSim(theParameterMap);
00225     theHBHEAmplifier->setTimeSlewSim(theTimeSlewSim);
00226     theHOAmplifier->setTimeSlewSim(theTimeSlewSim);
00227     theZDCAmplifier->setTimeSlewSim(theTimeSlewSim);
00228   }
00229 
00230   theHFDigitizer = new HFDigitizer(theHFResponse, theHFElectronicsSim, doEmpty);
00231   theZDCDigitizer = new ZDCDigitizer(theZDCResponse, theZDCElectronicsSim, doEmpty);
00232 
00233   bool doHPDNoise = ps.getParameter<bool>("doHPDNoise");
00234   if(doHPDNoise) {
00235     //edm::ParameterSet hpdNoisePset = ps.getParameter<edm::ParameterSet>("HPDNoiseLibrary");
00236     theNoiseGenerator = new HPDNoiseGenerator(ps); 
00237     if(theHBHEDigitizer) theHBHEDigitizer->setNoiseSignalGenerator(theNoiseGenerator);
00238     if(theHBHESiPMDigitizer) theHBHESiPMDigitizer->setNoiseSignalGenerator(theNoiseGenerator);
00239   }
00240 
00241   if(ps.getParameter<bool>("doIonFeedback") && theHBHEResponse)
00242   {
00243     theIonFeedback = new HPDIonFeedbackSim(ps, theShapes);
00244     theHBHEResponse->setPECorrection(theIonFeedback);
00245     if(ps.getParameter<bool>("doThermalNoise"))
00246     {
00247       theHBHEAmplifier->setIonFeedbackSim(theIonFeedback);
00248     }
00249   }
00250 
00251   if(ps.getParameter<bool>("injectTestHits") ){
00252     theNoiseHitGenerator = new HcalTestHitGenerator(ps);
00253     if(theHBHEDigitizer) theHBHEDigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00254     if(theHBHESiPMDigitizer) theHBHESiPMDigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00255     if(theHODigitizer) theHODigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00256     if(theHOSiPMDigitizer) theHOSiPMDigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00257     theHFDigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00258     theZDCDigitizer->setNoiseHitGenerator(theNoiseHitGenerator);
00259   }
00260 
00261   edm::Service<edm::RandomNumberGenerator> rng;
00262   if ( ! rng.isAvailable()) {
00263     throw cms::Exception("Configuration")
00264       << "HcalDigitizer requires the RandomNumberGeneratorService\n"
00265          "which is not present in the configuration file.  You must add the service\n"
00266          "in the configuration file or remove the modules that require it.";
00267   }
00268 
00269   CLHEP::HepRandomEngine& engine = rng->getEngine();
00270   if(theHBHEDigitizer) theHBHEDigitizer->setRandomEngine(engine);
00271   if(theHBHESiPMDigitizer) theHBHESiPMDigitizer->setRandomEngine(engine);
00272   if(theHODigitizer) theHODigitizer->setRandomEngine(engine);
00273   if(theHOSiPMDigitizer) theHOSiPMDigitizer->setRandomEngine(engine);
00274   if(theIonFeedback) theIonFeedback->setRandomEngine(engine);
00275   if(theTimeSlewSim) theTimeSlewSim->setRandomEngine(engine);
00276   theHFDigitizer->setRandomEngine(engine);
00277   theZDCDigitizer->setRandomEngine(engine);
00278 
00279   if (theHitCorrection!=0) theHitCorrection->setRandomEngine(engine);
00280 
00281   hitsProducer_ = ps.getParameter<std::string>("hitsProducer");
00282 }
00283 
00284 
00285 HcalDigitizer::~HcalDigitizer() {
00286   delete theHBHEDigitizer;
00287   delete theHBHESiPMDigitizer;
00288   delete theHODigitizer;
00289   delete theHOSiPMDigitizer;
00290   delete theHFDigitizer;
00291   delete theZDCDigitizer;
00292   delete theParameterMap;
00293   delete theHBHEResponse;
00294   delete theHBHESiPMResponse;
00295   delete theHOResponse;
00296   delete theHOSiPMResponse;
00297   delete theHFResponse;
00298   delete theZDCResponse;
00299   delete theHBHEElectronicsSim;
00300   delete theHFElectronicsSim;
00301   delete theHOElectronicsSim;
00302   delete theZDCElectronicsSim;
00303   delete theHBHEAmplifier;
00304   delete theHFAmplifier;
00305   delete theHOAmplifier;
00306   delete theZDCAmplifier;
00307   delete theCoderFactory;
00308   delete theHitCorrection;
00309   delete theNoiseGenerator;
00310 }
00311 
00312 
00313 void HcalDigitizer::setHBHENoiseSignalGenerator(HcalBaseSignalGenerator * noiseGenerator)
00314 {
00315   noiseGenerator->setParameterMap(theParameterMap);
00316   noiseGenerator->setElectronicsSim(theHBHEElectronicsSim);
00317   theHBHEDigitizer->setNoiseSignalGenerator(noiseGenerator);
00318   theHBHEAmplifier->setNoiseSignalGenerator(noiseGenerator);
00319 }
00320 
00321 void HcalDigitizer::setHFNoiseSignalGenerator(HcalBaseSignalGenerator * noiseGenerator)
00322 {
00323   noiseGenerator->setParameterMap(theParameterMap);
00324   noiseGenerator->setElectronicsSim(theHFElectronicsSim);
00325   theHFDigitizer->setNoiseSignalGenerator(noiseGenerator);
00326   theHFAmplifier->setNoiseSignalGenerator(noiseGenerator);
00327 }
00328 
00329 void HcalDigitizer::setHONoiseSignalGenerator(HcalBaseSignalGenerator * noiseGenerator)
00330 {
00331   noiseGenerator->setParameterMap(theParameterMap);
00332   noiseGenerator->setElectronicsSim(theHOElectronicsSim);
00333   theHODigitizer->setNoiseSignalGenerator(noiseGenerator);
00334   theHOAmplifier->setNoiseSignalGenerator(noiseGenerator);
00335 }
00336 
00337 void HcalDigitizer::setZDCNoiseSignalGenerator(HcalBaseSignalGenerator * noiseGenerator)
00338 {
00339   noiseGenerator->setParameterMap(theParameterMap);
00340   noiseGenerator->setElectronicsSim(theZDCElectronicsSim);
00341   theZDCDigitizer->setNoiseSignalGenerator(noiseGenerator);
00342   theZDCAmplifier->setNoiseSignalGenerator(noiseGenerator);
00343 }
00344 
00345 void HcalDigitizer::initializeEvent(edm::Event const& e, edm::EventSetup const& eventSetup) {
00346   // get the appropriate gains, noises, & widths for this event
00347   edm::ESHandle<HcalDbService> conditions;
00348   eventSetup.get<HcalDbRecord>().get(conditions);
00349   theHBHEAmplifier->setDbService(conditions.product());
00350   theHFAmplifier->setDbService(conditions.product());
00351   theHOAmplifier->setDbService(conditions.product());
00352   theZDCAmplifier->setDbService(conditions.product());
00353 
00354   theCoderFactory->setDbService(conditions.product());
00355   theParameterMap->setDbService(conditions.product());
00356 
00357    edm::ESHandle<HcalCholeskyMatrices> refCholesky;
00358    eventSetup.get<HcalCholeskyMatricesRcd>().get(refCholesky);
00359    const HcalCholeskyMatrices * myCholesky = refCholesky.product();
00360 
00361    edm::ESHandle<HcalPedestals> pedshandle;
00362    eventSetup.get<HcalPedestalsRcd>().get(pedshandle);
00363    const HcalPedestals *  myADCPedestals = pedshandle.product();
00364 
00365   theHBHEAmplifier->setCholesky(myCholesky);
00366   theHFAmplifier->setCholesky(myCholesky);
00367   theHOAmplifier->setCholesky(myCholesky);
00368 
00369   theHBHEAmplifier->setADCPeds(myADCPedestals);
00370   theHFAmplifier->setADCPeds(myADCPedestals);
00371   theHOAmplifier->setADCPeds(myADCPedestals);
00372 
00373   if(theHitCorrection != 0) {
00374     theHitCorrection->clear();
00375   }
00376 }
00377 
00378 void HcalDigitizer::accumulateCaloHits(edm::Handle<std::vector<PCaloHit> > const& hcalHandle, edm::Handle<std::vector<PCaloHit> > const& zdcHandle, int bunchCrossing) {
00379   // Step A: pass in inputs, and accumulate digirs
00380   if(isHCAL) {
00381     std::vector<PCaloHit> const& hcalHits = *hcalHandle.product();
00382     if(theHitCorrection != 0) {
00383       theHitCorrection->fillChargeSums(hcalHits);
00384     }
00385 
00386     if(hbhegeo) {
00387       if(theHBHEDigitizer) theHBHEDigitizer->add(hcalHits, bunchCrossing);
00388       if(theHBHESiPMDigitizer) theHBHESiPMDigitizer->add(hcalHits, bunchCrossing);
00389     }
00390 
00391     if(hogeo) {
00392       if(theHODigitizer) theHODigitizer->add(hcalHits, bunchCrossing);
00393       if(theHOSiPMDigitizer) theHOSiPMDigitizer->add(hcalHits, bunchCrossing);
00394     }
00395 
00396     if(hfgeo) {
00397       theHFDigitizer->add(hcalHits, bunchCrossing);
00398     } 
00399   } else {
00400     edm::LogInfo("HcalDigitizer") << "We don't have HCAL hit collection available ";
00401   }
00402 
00403   if(isZDC) {
00404     if(zdcgeo) {
00405       theZDCDigitizer->add(*zdcHandle.product(), bunchCrossing);
00406     } 
00407   } else {
00408     edm::LogInfo("HcalDigitizer") << "We don't have ZDC hit collection available ";
00409   }
00410 }
00411 
00412 void HcalDigitizer::accumulate(edm::Event const& e, edm::EventSetup const& eventSetup) {
00413   // Step A: Get Inputs
00414   edm::InputTag zdcTag(hitsProducer_, "ZDCHITS");
00415   edm::Handle<std::vector<PCaloHit> > zdcHandle;
00416   e.getByLabel(zdcTag, zdcHandle);
00417   isZDC = zdcHandle.isValid();
00418 
00419   edm::InputTag hcalTag(hitsProducer_, "HcalHits");
00420   edm::Handle<std::vector<PCaloHit> > hcalHandle;
00421   e.getByLabel(hcalTag, hcalHandle);
00422   isHCAL = hcalHandle.isValid();
00423 
00424   accumulateCaloHits(hcalHandle, zdcHandle, 0);
00425 }
00426 
00427 void HcalDigitizer::accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& eventSetup) {
00428   // Step A: Get Inputs
00429   edm::InputTag zdcTag(hitsProducer_, "ZDCHITS");
00430   edm::Handle<std::vector<PCaloHit> > zdcHandle;
00431   e.getByLabel(zdcTag, zdcHandle);
00432   isZDC = zdcHandle.isValid();
00433 
00434   edm::InputTag hcalTag(hitsProducer_, "HcalHits");
00435   edm::Handle<std::vector<PCaloHit> > hcalHandle;
00436   e.getByLabel(hcalTag, hcalHandle);
00437   isHCAL = hcalHandle.isValid();
00438 
00439   accumulateCaloHits(hcalHandle, zdcHandle, e.bunchCrossing());
00440 }
00441 
00442 void HcalDigitizer::finalizeEvent(edm::Event& e, const edm::EventSetup& eventSetup) {
00443 
00444   // Step B: Create empty output
00445   std::auto_ptr<HBHEDigiCollection> hbheResult(new HBHEDigiCollection());
00446   std::auto_ptr<HODigiCollection> hoResult(new HODigiCollection());
00447   std::auto_ptr<HFDigiCollection> hfResult(new HFDigiCollection());
00448   std::auto_ptr<ZDCDigiCollection> zdcResult(new ZDCDigiCollection());
00449 
00450   // Step C: Invoke the algorithm, getting back outputs.
00451   if(isHCAL&&hbhegeo)
00452   {
00453     if(theHBHEDigitizer) theHBHEDigitizer->run(*hbheResult);
00454     if(theHBHESiPMDigitizer) theHBHESiPMDigitizer->run(*hbheResult);
00455   }
00456   if(isHCAL&&hogeo)
00457   {
00458     if(theHODigitizer) theHODigitizer->run(*hoResult);
00459     if(theHOSiPMDigitizer) theHOSiPMDigitizer->run(*hoResult);
00460   }
00461   if(isHCAL&&hfgeo)
00462     theHFDigitizer->run(*hfResult);  
00463   if(isZDC&&zdcgeo) 
00464     theZDCDigitizer->run(*zdcResult);
00465   
00466   edm::LogInfo("HcalDigitizer") << "HCAL HBHE digis : " << hbheResult->size();
00467   edm::LogInfo("HcalDigitizer") << "HCAL HO digis   : " << hoResult->size();
00468   edm::LogInfo("HcalDigitizer") << "HCAL HF digis   : " << hfResult->size();
00469   edm::LogInfo("HcalDigitizer") << "HCAL ZDC digis   : " << zdcResult->size();
00470   // Step D: Put outputs into event
00471   e.put(hbheResult);
00472   e.put(hoResult);
00473   e.put(hfResult);
00474   e.put(zdcResult);
00475 
00476   if(theHitCorrection) {
00477     theHitCorrection->clear();
00478   }
00479 }
00480 
00481 
00482 void HcalDigitizer::beginRun(const edm::EventSetup & es)
00483 {
00484   checkGeometry(es);
00485   theShapes->beginRun(es);
00486 }
00487 
00488 
00489 void HcalDigitizer::endRun()
00490 {
00491   theShapes->endRun();
00492 }
00493 
00494 
00495 void HcalDigitizer::checkGeometry(const edm::EventSetup & eventSetup) {
00496   // TODO find a way to avoid doing this every event
00497   edm::ESHandle<CaloGeometry> geometry;
00498   eventSetup.get<CaloGeometryRecord>().get(geometry);
00499   // See if it's been updated
00500   if(&*geometry != theGeometry)
00501   {
00502     theGeometry = &*geometry;
00503     updateGeometry(eventSetup);
00504   }
00505 }
00506 
00507 
00508 void  HcalDigitizer::updateGeometry(const edm::EventSetup & eventSetup)
00509 {
00510   if(theHBHEResponse) theHBHEResponse->setGeometry(theGeometry);
00511   if(theHBHESiPMResponse) theHBHESiPMResponse->setGeometry(theGeometry);
00512   if(theHOResponse) theHOResponse->setGeometry(theGeometry);
00513   if(theHOSiPMResponse) theHOSiPMResponse->setGeometry(theGeometry);
00514   theHFResponse->setGeometry(theGeometry);
00515   theZDCResponse->setGeometry(theGeometry);
00516 
00517   const std::vector<DetId>& hbCells = theGeometry->getValidDetIds(DetId::Hcal, HcalBarrel);
00518   const std::vector<DetId>& heCells = theGeometry->getValidDetIds(DetId::Hcal, HcalEndcap);
00519   const std::vector<DetId>& hoCells = theGeometry->getValidDetIds(DetId::Hcal, HcalOuter);
00520   const std::vector<DetId>& hfCells = theGeometry->getValidDetIds(DetId::Hcal, HcalForward);
00521   const std::vector<DetId>& zdcCells = theGeometry->getValidDetIds(DetId::Calo, HcalZDCDetId::SubdetectorId);
00522   //const std::vector<DetId>& hcalTrigCells = geometry->getValidDetIds(DetId::Hcal, HcalTriggerTower);
00523   //const std::vector<DetId>& hcalCalib = geometry->getValidDetIds(DetId::Calo, HcalCastorDetId::SubdetectorId);
00524   //std::cout<<"HcalDigitizer::CheckGeometry number of cells: "<<zdcCells.size()<<std::endl;
00525   if(zdcCells.empty()) zdcgeo = false;
00526   if(hbCells.empty() && heCells.empty()) hbhegeo = false;
00527   if(hoCells.empty()) hogeo = false;
00528   if(hfCells.empty()) hfgeo = false;
00529   // combine HB & HE
00530 
00531 
00532   theHBHEDetIds = hbCells;
00533   theHBHEDetIds.insert(theHBHEDetIds.end(), heCells.begin(), heCells.end());
00534 
00535   HcalDigitizerImpl::fillCells(theHBHEDetIds, theHBHEDigitizer, theHBHESiPMDigitizer);
00536   //HcalDigitizerImpl::fillCells(hoCells, theHODigitizer, theHOSiPMDigitizer);
00537   buildHOSiPMCells(hoCells, eventSetup);
00538   theHFDigitizer->setDetIds(hfCells);
00539   theZDCDigitizer->setDetIds(zdcCells); 
00540 }
00541 
00542 
00543 void HcalDigitizer::buildHOSiPMCells(const std::vector<DetId>& allCells, const edm::EventSetup & eventSetup)
00544 {
00545   // all HPD
00546   if(theHOSiPMCode == 0)
00547   {
00548     theHODigitizer->setDetIds(allCells);
00549   }
00550   else if(theHOSiPMCode == 1)
00551   {
00552     theHOSiPMDigitizer->setDetIds(allCells);
00553     // FIXME pick Zecotek or hamamatsu?
00554   } 
00555   else if(theHOSiPMCode == 2)
00556   {
00557     std::vector<HcalDetId> zecotekDetIds, hamamatsuDetIds;
00558     edm::ESHandle<HcalMCParams> p;
00559     eventSetup.get<HcalMCParamsRcd>().get(p);
00560     edm::ESHandle<HcalTopology> htopo;
00561     eventSetup.get<IdealGeometryRecord>().get(htopo);
00562    
00563     HcalMCParams mcParams(*p.product());
00564     if (mcParams.topo()==0) {
00565       mcParams.setTopo(htopo.product());
00566     }
00567 
00568     for(std::vector<DetId>::const_iterator detItr = allCells.begin();
00569         detItr != allCells.end(); ++detItr)
00570     {
00571       int shapeType = mcParams.getValues(*detItr)->signalShape();
00572       if(shapeType == HcalShapes::ZECOTEK) {
00573         zecotekDetIds.emplace_back(*detItr);
00574         theHOSiPMDetIds.push_back(*detItr);
00575       }
00576       else if(shapeType == HcalShapes::HAMAMATSU) {
00577         hamamatsuDetIds.emplace_back(*detItr);
00578         theHOSiPMDetIds.push_back(*detItr);
00579       }
00580       else {
00581         theHOHPDDetIds.push_back(*detItr);
00582       }
00583     }
00584 
00585     assert(theHODigitizer);
00586     assert(theHOSiPMDigitizer);
00587     theHODigitizer->setDetIds(theHOHPDDetIds);
00588     theHOSiPMDigitizer->setDetIds(theHOSiPMDetIds);
00589     theHOSiPMHitFilter.setDetIds(theHOSiPMDetIds);
00590     // FIXME not applying a HitFilter to the HPDs, for now
00591     theParameterMap->setHOZecotekDetIds(zecotekDetIds);
00592     theParameterMap->setHOHamamatsuDetIds(hamamatsuDetIds);
00593 
00594     // make sure we don't got through this exercise again
00595     theHOSiPMCode = -2;
00596   }
00597 }
00598 
00599       
00600     
00601