CMSSW Reference Manual

00001 
00002 #include "EcalDigiProducer.h"
00003 #include "DataFormats/Common/interface/EDProduct.h"
00004 #include "FWCore/Framework/interface/ESHandle.h"
00005 #include "DataFormats/Common/interface/Handle.h"
00006 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 #include "DataFormats/Provenance/interface/Provenance.h"
00009 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
00010 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
00011 #include "SimCalorimetry/CaloSimAlgos/interface/CaloDigiCollectionSorter.h"
00012 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00013 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
00014 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
00015 
00016 EcalDigiProducer::EcalDigiProducer(const edm::ParameterSet& params) 
00017   :  theGeometry(0)
00018 {
00019 
00021 
00022   EBdigiCollection_ = params.getParameter<std::string>("EBdigiCollection");
00023   EEdigiCollection_ = params.getParameter<std::string>("EEdigiCollection");
00024   ESdigiCollection_ = params.getParameter<std::string>("ESdigiCollection");
00025 
00026   hitsProducer_ = params.getParameter<std::string>("hitsProducer");
00027   
00028   produces<EBDigiCollection>(EBdigiCollection_);
00029   produces<EEDigiCollection>(EEdigiCollection_);
00030   produces<ESDigiCollection>(ESdigiCollection_);
00031 
00032 
00033   // initialize the default valuer for hardcoded parameters and the EB/EE shape
00034 
00035   bool addNoise = params.getParameter<bool>("doNoise"); 
00036   double simHitToPhotoelectronsBarrel = params.getParameter<double>("simHitToPhotoelectronsBarrel");
00037   double simHitToPhotoelectronsEndcap = params.getParameter<double>("simHitToPhotoelectronsEndcap");
00038   double photoelectronsToAnalogBarrel = params.getParameter<double>("photoelectronsToAnalogBarrel");
00039   double photoelectronsToAnalogEndcap = params.getParameter<double>("photoelectronsToAnalogEndcap");
00040   double samplingFactor = params.getParameter<double>("samplingFactor");
00041   double timePhase = params.getParameter<double>("timePhase");
00042   int readoutFrameSize = params.getParameter<int>("readoutFrameSize");
00043   int binOfMaximum = params.getParameter<int>("binOfMaximum");
00044   bool doPhotostatistics = params.getParameter<bool>("doPhotostatistics");
00045   bool syncPhase = params.getParameter<bool>("syncPhase");
00046 
00047   theParameterMap = new EcalSimParameterMap(simHitToPhotoelectronsBarrel, simHitToPhotoelectronsEndcap, 
00048                                             photoelectronsToAnalogBarrel, photoelectronsToAnalogEndcap, 
00049                                             samplingFactor, timePhase, readoutFrameSize, binOfMaximum,
00050                                             doPhotostatistics, syncPhase);
00051 
00052   theEcalShape = new EcalShape(timePhase);
00053   
00054   int ESGain = params.getParameter<int>("ESGain");
00055   theESShape = new ESShape(ESGain);
00056 
00057   theEcalResponse = new CaloHitResponse(theParameterMap, theEcalShape);
00058   theESResponse = new CaloHitResponse(theParameterMap, theESShape);
00059 
00060   // further phase for cosmics studies
00061   cosmicsPhase = params.getParameter<bool>("cosmicsPhase");
00062   cosmicsShift = params.getParameter<double>("cosmicsShift");
00063   if (cosmicsPhase) {
00064     theEcalResponse->setPhaseShift(1.+cosmicsShift);
00065   }
00066 
00067   EcalCorrMatrix thisMatrix;
00068 
00069   std::vector<double> corrNoiseMatrix = params.getParameter< std::vector<double> >("CorrelatedNoiseMatrix");
00070   if ( corrNoiseMatrix.size() == (unsigned int)(readoutFrameSize*readoutFrameSize) ) {
00071     for ( int row = 0 ; row < readoutFrameSize; ++row ) {
00072       for ( int column = 0 ; column < readoutFrameSize; ++column ) {
00073         int index = column + readoutFrameSize*row;
00074         thisMatrix(row,column) = corrNoiseMatrix[index];
00075       }
00076     }
00077   }
00078   theNoiseMatrix = new EcalCorrelatedNoiseMatrix(thisMatrix);
00079 
00080   theCorrNoise = new CorrelatedNoisifier<EcalCorrMatrix>(thisMatrix);
00081 
00082   theCoder = new EcalCoder(addNoise, theCorrNoise);
00083   bool applyConstantTerm = params.getParameter<bool>("applyConstantTerm");
00084   double rmsConstantTerm = params.getParameter<double> ("ConstantTerm");
00085   theElectronicsSim = new EcalElectronicsSim(theParameterMap, theCoder, applyConstantTerm, rmsConstantTerm);
00086 
00087   doFast = params.getParameter<bool>("doFast");
00088   bool addESNoise = params.getParameter<bool>("doESNoise");
00089   double ESNoiseSigma = params.getParameter<double> ("ESNoiseSigma");
00090   int ESBaseline  = params.getParameter<int>("ESBaseline");
00091   double ESMIPADC = params.getParameter<double>("ESMIPADC");
00092   double ESMIPkeV = params.getParameter<double>("ESMIPkeV");
00093   int numESdetId  = params.getParameter<int>("numESdetId");
00094   double zsThreshold = params.getParameter<double>("zsThreshold");
00095   std::string refFile = params.getParameter<std::string>("refHistosFile");
00096 
00097   theESElectronicsSim     = 0;
00098   theESElectronicsSimFast = 0;
00099   if (!doFast) { theESElectronicsSim     = new ESElectronicsSim(addESNoise, ESNoiseSigma, ESGain, ESBaseline, ESMIPADC, ESMIPkeV); }
00100   if ( doFast) { theESElectronicsSimFast = new ESElectronicsSimFast(addESNoise, ESNoiseSigma, ESGain, ESBaseline, ESMIPADC, ESMIPkeV);  }
00101 
00102   theBarrelDigitizer = new EBDigitizer(theEcalResponse, theElectronicsSim, addNoise);
00103   theEndcapDigitizer = new EEDigitizer(theEcalResponse, theElectronicsSim, addNoise);
00104 
00105   theESDigitizer = 0;
00106   theESDigitizerFast = 0;
00107   if (!doFast){ theESDigitizer = new ESDigitizer(theESResponse, theESElectronicsSim, addESNoise); }
00108   if ( doFast){ theESDigitizerFast = new ESFastTDigitizer(theESResponse, theESElectronicsSimFast, addESNoise, numESdetId, zsThreshold, refFile);}
00109 
00110   // not containment corrections
00111   EBs25notCont = params.getParameter<double>("EBs25notContainment");
00112   EEs25notCont = params.getParameter<double>("EEs25notContainment");
00113 }
00114 
00115 
00116 EcalDigiProducer::~EcalDigiProducer() 
00117 {
00118   if (theParameterMap)  { delete theParameterMap; }
00119   if (theEcalShape)     { delete theEcalShape; }
00120   if (theESShape)       { delete theESShape; }
00121   if (theEcalResponse)  { delete theEcalResponse; }
00122   if (theESResponse)    { delete theESResponse; }
00123   if (theCorrNoise)     { delete theCorrNoise; }
00124   if (theNoiseMatrix)   { delete theNoiseMatrix; }
00125   if (theCoder)         { delete theCoder; }
00126   if (theElectronicsSim){ delete theElectronicsSim; }
00127   if (theESElectronicsSim)    { delete theESElectronicsSim; }
00128   if (theESElectronicsSimFast){ delete theESElectronicsSimFast; }
00129   if (theBarrelDigitizer){ delete theBarrelDigitizer; }
00130   if (theEndcapDigitizer){ delete theEndcapDigitizer; }
00131   if (theESDigitizer)    { delete theESDigitizer; }
00132   if (theESDigitizerFast){ delete theESDigitizerFast; }
00133 }
00134 
00135 
00136 void EcalDigiProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) 
00137 {
00138 
00139   // Step A: Get Inputs
00140 
00141   checkGeometry(eventSetup);
00142   checkCalibrations(eventSetup);
00143 
00144   // Get input
00145   edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
00146 
00147   const std::vector<DetId>& theBarrelDets =  theGeometry->getValidDetIds(DetId::Ecal, EcalBarrel);
00148   const std::vector<DetId>& theEndcapDets =  theGeometry->getValidDetIds(DetId::Ecal, EcalEndcap);
00149   const std::vector<DetId>& theESDets     =  theGeometry->getValidDetIds(DetId::Ecal, EcalPreshower);
00150 
00151   // test access to SimHits
00152   const std::string barrelHitsName(hitsProducer_+"EcalHitsEB");
00153   const std::string endcapHitsName(hitsProducer_+"EcalHitsEE");
00154   const std::string preshowerHitsName(hitsProducer_+"EcalHitsES");
00155 
00156   bool isEB = true;
00157   event.getByLabel("mix",barrelHitsName,crossingFrame);
00158   MixCollection<PCaloHit> * EBHits = 0 ;
00159   if (crossingFrame.isValid()) { 
00160     EBHits = new MixCollection<PCaloHit>(crossingFrame.product());
00161   }
00162   else { 
00163     edm::LogError("EcalDigiProducer") << "Error! can't get the product " << barrelHitsName.c_str() ;
00164     isEB = false;
00165   }
00166   if ( ! EBHits->inRegistry() || theBarrelDets.size() == 0 ) isEB = false;
00167 
00168   bool isEE = true;
00169   event.getByLabel("mix",endcapHitsName,crossingFrame);
00170   MixCollection<PCaloHit> * EEHits = 0 ;
00171   if (crossingFrame.isValid()) {
00172     EEHits = new MixCollection<PCaloHit>(crossingFrame.product());
00173   }
00174   else {
00175     edm::LogError("EcalDigiProducer") << "Error! can't get the product " << endcapHitsName.c_str() ;
00176     isEE = false;
00177   }
00178   if ( ! EEHits->inRegistry() || theEndcapDets.size() == 0 ) isEE = false;
00179 
00180   bool isES = true;
00181   event.getByLabel("mix",preshowerHitsName,crossingFrame);
00182   MixCollection<PCaloHit> * ESHits = 0 ;
00183   if (crossingFrame.isValid()) {
00184     ESHits = new MixCollection<PCaloHit>(crossingFrame.product());
00185   }
00186   else { 
00187     edm::LogError("EcalDigiProducer") << "Error! can't get the product " << preshowerHitsName.c_str() ;
00188     isES = false; 
00189   }    
00190   if ( ! ESHits->inRegistry() || theESDets.size() == 0 ) isES = false;
00191 
00192   // Step B: Create empty output
00193   std::auto_ptr<EBDigiCollection> barrelResult(new EBDigiCollection());
00194   std::auto_ptr<EEDigiCollection> endcapResult(new EEDigiCollection());
00195   std::auto_ptr<ESDigiCollection> preshowerResult(new ESDigiCollection());
00196 
00197   // run the algorithm
00198 
00199   CaloDigiCollectionSorter sorter(5);
00200 
00201   if ( isEB ) {
00202 
00203     std::auto_ptr<MixCollection<PCaloHit> >  barrelHits( EBHits );
00204     theBarrelDigitizer->run(*barrelHits, *barrelResult);
00205     edm::LogInfo("DigiInfo") << "EB Digis: " << barrelResult->size();
00206 
00207     /*
00208     std::vector<EBDataFrame> sortedDigisEB = sorter.sortedVector(*barrelResult);
00209     LogDebug("EcalDigi") << "Top 10 EB digis";
00210     for(int i = 0; i < std::min(10,(int) sortedDigisEB.size()); ++i) 
00211       {
00212         LogDebug("EcalDigi") << sortedDigisEB[i];
00213       }
00214     */
00215   }
00216 
00217   if ( isEE ) {
00218 
00219     std::auto_ptr<MixCollection<PCaloHit> >  endcapHits( EEHits );
00220     theEndcapDigitizer->run(*endcapHits, *endcapResult);
00221     edm::LogInfo("EcalDigi") << "EE Digis: " << endcapResult->size();
00222 
00223     /*
00224     std::vector<EEDataFrame> sortedDigisEE = sorter.sortedVector(*endcapResult);
00225     LogDebug("EcalDigi")  << "Top 10 EE digis";
00226     for(int i = 0; i < std::min(10,(int) sortedDigisEE.size()); ++i) 
00227       {
00228         LogDebug("EcalDigi") << sortedDigisEE[i];
00229       }
00230     */
00231   }
00232 
00233   if ( isES ) {
00234 
00235     std::auto_ptr<MixCollection<PCaloHit> >  preshowerHits( ESHits );
00236     if (!doFast) { theESDigitizer->run(*preshowerHits, *preshowerResult); }
00237     if ( doFast) { theESDigitizerFast->run(*preshowerHits, *preshowerResult); }
00238     edm::LogInfo("EcalDigi") << "ES Digis: " << preshowerResult->size();
00239     
00240 //   CaloDigiCollectionSorter sorter_es(7);
00241 //   std::vector<ESDataFrame> sortedDigis_es = sorter_es.sortedVector(*preshowerResult);
00242 //   LogDebug("DigiDump") << "List all ES digis";
00243 //   for(int i = 0; i < sortedDigis_es.size(); ++i) 
00244 //     {
00245 //       LogDebug("DigiDump") << sortedDigis_es[i];
00246 //     }
00247   }
00248 
00249   // Step D: Put outputs into event
00250   event.put(barrelResult, EBdigiCollection_);
00251   event.put(endcapResult, EEdigiCollection_);
00252   event.put(preshowerResult, ESdigiCollection_);
00253 
00254 }
00255 
00256 
00257 void  EcalDigiProducer::checkCalibrations(const edm::EventSetup & eventSetup) 
00258 {
00259 
00260   // Pedestals from event setup
00261 
00262   edm::ESHandle<EcalPedestals> dbPed;
00263   eventSetup.get<EcalPedestalsRcd>().get( dbPed );
00264   const EcalPedestals* thePedestals=dbPed.product();
00265   
00266   theCoder->setPedestals( thePedestals );
00267 
00268   // ADC -> GeV Scale
00269   edm::ESHandle<EcalADCToGeVConstant> pAgc;
00270   eventSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
00271   const EcalADCToGeVConstant* agc = pAgc.product();
00272   
00273   // Gain Ratios
00274   edm::ESHandle<EcalGainRatios> pRatio;
00275   eventSetup.get<EcalGainRatiosRcd>().get(pRatio);
00276   const EcalGainRatios* gr = pRatio.product();
00277 
00278   theCoder->setGainRatios( gr );
00279 
00280   EcalMGPAGainRatio * defaultRatios = new EcalMGPAGainRatio();
00281 
00282   double theGains[theCoder->NGAINS+1];
00283   theGains[0] = 0.;
00284   theGains[3] = 1.;
00285   theGains[2] = defaultRatios->gain6Over1() ;
00286   theGains[1] = theGains[2]*(defaultRatios->gain12Over6()) ;
00287 
00288   LogDebug("EcalDigi") << " Gains: " << "\n" << " g1 = " << theGains[1] << "\n" << " g2 = " << theGains[2] << "\n" << " g3 = " << theGains[3];
00289 
00290   delete defaultRatios;
00291 
00292   const double EBscale = (agc->getEBValue())*theGains[1]*(theCoder->MAXADC)*EBs25notCont;
00293   LogDebug("EcalDigi") << " GeV/ADC = " << agc->getEBValue() << "\n" << " notCont = " << EBs25notCont << "\n" << " saturation for EB = " << EBscale << ", " << EBs25notCont;
00294   const double EEscale = (agc->getEEValue())*theGains[1]*(theCoder->MAXADC)*EEs25notCont;
00295   LogDebug("EcalDigi") << " GeV/ADC = " << agc->getEEValue() << "\n" << " notCont = " << EEs25notCont << "\n" << " saturation for EB = " << EEscale << ", " << EEs25notCont;
00296   theCoder->setFullScaleEnergy( EBscale , EEscale );
00297 
00298 }
00299 
00300 
00301 void EcalDigiProducer::checkGeometry(const edm::EventSetup & eventSetup) 
00302 {
00303   // TODO find a way to avoid doing this every event
00304   edm::ESHandle<CaloGeometry> hGeometry;
00305   eventSetup.get<CaloGeometryRecord>().get(hGeometry);
00306 
00307   const CaloGeometry * pGeometry = &*hGeometry;
00308   
00309   // see if we need to update
00310   if(pGeometry != theGeometry) {
00311     theGeometry = pGeometry;
00312     updateGeometry();
00313   }
00314 }
00315 
00316 
00317 void EcalDigiProducer::updateGeometry() {
00318   theEcalResponse->setGeometry(theGeometry);
00319   theESResponse->setGeometry(theGeometry);
00320 
00321   const std::vector<DetId>& theBarrelDets =  theGeometry->getValidDetIds(DetId::Ecal, EcalBarrel);
00322   const std::vector<DetId>& theEndcapDets =  theGeometry->getValidDetIds(DetId::Ecal, EcalEndcap);
00323   const std::vector<DetId>& theESDets     =  theGeometry->getValidDetIds(DetId::Ecal, EcalPreshower);
00324 
00325   edm::LogInfo("EcalDigi") << "deb geometry: " << "\n" 
00326                       << "\t barrel: " << theBarrelDets.size () << "\n"
00327                       << "\t endcap: " << theEndcapDets.size () << "\n"
00328                       << "\t preshower: " << theESDets.size();
00329   
00330   theBarrelDigitizer->setDetIds(theBarrelDets);
00331   theEndcapDigitizer->setDetIds(theEndcapDets);
00332   if (!doFast) { theESDigitizer->setDetIds(theESDets); }
00333   if ( doFast) { theESDigitizerFast->setDetIds(theESDets); }
00334 }
00335
EcalDigiProducer.cc
