CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/SimGeneral/DataMixingModule/plugins/DataMixingEMWorker.cc

Go to the documentation of this file.
00001 // File: DataMixingEMWorker.cc
00002 // Description:  see DataMixingEMWorker.h
00003 // Author:  Mike Hildreth, University of Notre Dame
00004 //
00005 //--------------------------------------------
00006 
00007 #include <map>
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "FWCore/Utilities/interface/EDMException.h"
00010 #include "FWCore/Framework/interface/ConstProductRegistry.h"
00011 #include "FWCore/ServiceRegistry/interface/Service.h"
00012 #include "DataFormats/Common/interface/Handle.h"
00013 #include "DataFormats/Provenance/interface/Provenance.h"
00014 #include "DataFormats/Provenance/interface/BranchDescription.h"
00015 //
00016 //
00017 #include "DataMixingEMWorker.h"
00018 
00019 
00020 using namespace std;
00021 
00022 namespace edm
00023 {
00024 
00025   // Virtual constructor
00026 
00027   DataMixingEMWorker::DataMixingEMWorker() { }
00028 
00029   // Constructor 
00030   DataMixingEMWorker::DataMixingEMWorker(const edm::ParameterSet& ps) : 
00031                                                             label_(ps.getParameter<std::string>("Label"))
00032 
00033   {                                                         
00034 
00035     // get the subdetector names
00036     //    this->getSubdetectorNames();  //something like this may be useful to check what we are supposed to do...
00037 
00038     // declare the products to produce, retrieve
00039 
00040     EBProducerSig_ = ps.getParameter<edm::InputTag>("EBProducerSig");
00041     EEProducerSig_ = ps.getParameter<edm::InputTag>("EEProducerSig");
00042     ESProducerSig_ = ps.getParameter<edm::InputTag>("ESProducerSig");
00043 
00044 
00045     EBrechitCollectionSig_ = ps.getParameter<edm::InputTag>("EBrechitCollectionSig");
00046     EErechitCollectionSig_ = ps.getParameter<edm::InputTag>("EErechitCollectionSig");
00047     ESrechitCollectionSig_ = ps.getParameter<edm::InputTag>("ESrechitCollectionSig");
00048 
00049     EBPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("EBPileRecHitInputTag");
00050     EEPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("EEPileRecHitInputTag");
00051     ESPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("ESPileRecHitInputTag");
00052 
00053 
00054     EBRecHitCollectionDM_        = ps.getParameter<std::string>("EBRecHitCollectionDM");
00055     EERecHitCollectionDM_        = ps.getParameter<std::string>("EERecHitCollectionDM");
00056     ESRecHitCollectionDM_        = ps.getParameter<std::string>("ESRecHitCollectionDM");
00057    //   nMaxPrintout_            = ps.getUntrackedParameter<int>("nMaxPrintout",10);
00058 
00059    //EBalgo_ = new EcalRecHitSimpleAlgo();
00060    //EEalgo_ = new EcalRecHitSimpleAlgo();
00061 
00062    // don't think I can "produce" in a sub-class...
00063 
00064    //produces< EBRecHitCollection >(EBRecHitCollectionDM_);
00065    //produces< EERecHitCollection >(EERecHitCollectionDM_);
00066 
00067   }
00068                
00069 
00070   // Virtual destructor needed.
00071   DataMixingEMWorker::~DataMixingEMWorker() { 
00072   }  
00073 
00074   void DataMixingEMWorker::addEMSignals(const edm::Event &e) { 
00075     // fill in maps of hits
00076 
00077     LogInfo("DataMixingEMWorker")<<"===============> adding MC signals for "<<e.id();
00078 
00079     // EB first
00080 
00081    Handle< EBRecHitCollection > pEBRecHits;
00082 
00083    const EBRecHitCollection*  EBRecHits = 0;
00084 
00085    if(e.getByLabel(EBProducerSig_, pEBRecHits) ){
00086      EBRecHits = pEBRecHits.product(); // get a ptr to the product
00087      LogDebug("DataMixingEMWorker") << "total # EB rechits: " << EBRecHits->size();
00088    }
00089    
00090  
00091    if (EBRecHits)
00092      {
00093        // loop over rechits, storing them in a map so we can add pileup later
00094        for(EBRecHitCollection::const_iterator it  = EBRecHits->begin(); 
00095            it != EBRecHits->end(); ++it) {
00096 
00097          EBRecHitStorage_.insert(EBRecHitMap::value_type( ( it->id() ), *it ));
00098          
00099          LogDebug("DataMixingEMWorker") << "processed EBRecHit with rawId: "
00100                                       << it->id().rawId() << "\n"
00101                                       << " rechit energy: " << it->energy();
00102 
00103        }
00104      }
00105 
00106    // EE next
00107 
00108    Handle< EERecHitCollection > pEERecHits;
00109 
00110    const EERecHitCollection*  EERecHits = 0;
00111 
00112    
00113    if(e.getByLabel(EEProducerSig_, pEERecHits) ){
00114      EERecHits = pEERecHits.product(); // get a ptr to the product
00115      LogDebug("DataMixingEMWorker") << "total # EE rechits: " << EERecHits->size();
00116    } 
00117    
00118  
00119    if (EERecHits)
00120      {
00121        // loop over rechits, storing them in a map so we can add pileup later
00122        for(EERecHitCollection::const_iterator it  = EERecHits->begin(); 
00123            it != EERecHits->end(); ++it) {
00124 
00125          EERecHitStorage_.insert(EERecHitMap::value_type( ( it->id() ), *it ));
00126 #ifdef DEBUG     
00127          LogDebug("DataMixingEMWorker") << "processed EERecHit with rawId: "
00128                                       << it->id().rawId() << "\n"
00129                                       << " rechit energy: " << it->energy();
00130 #endif
00131 
00132        }
00133      }
00134    // ES next
00135 
00136    Handle< ESRecHitCollection > pESRecHits;
00137 
00138    const ESRecHitCollection*  ESRecHits = 0;
00139 
00140    
00141    if(e.getByLabel( ESProducerSig_, pESRecHits) ){
00142      ESRecHits = pESRecHits.product(); // get a ptr to the product
00143 #ifdef DEBUG
00144      LogDebug("DataMixingEMWorker") << "total # ES rechits: " << ESRecHits->size();
00145 #endif
00146    } 
00147    
00148  
00149    if (ESRecHits)
00150      {
00151        // loop over rechits, storing them in a map so we can add pileup later
00152        for(ESRecHitCollection::const_iterator it  = ESRecHits->begin(); 
00153            it != ESRecHits->end(); ++it) {
00154 
00155          ESRecHitStorage_.insert(ESRecHitMap::value_type( ( it->id() ), *it ));
00156          
00157 #ifdef DEBUG     
00158          LogDebug("DataMixingEMWorker") << "processed ESRecHit with rawId: "
00159                                       << it->id().rawId() << "\n"
00160                                       << " rechit energy: " << it->energy();
00161 #endif
00162 
00163        }
00164      }
00165     
00166   } // end of addEMSignals
00167 
00168   void DataMixingEMWorker::addEMPileups(const int bcr, const EventPrincipal *ep, unsigned int eventNr) {
00169   
00170     LogInfo("DataMixingEMWorker") <<"\n===============> adding pileups from event  "<<ep->id()<<" for bunchcrossing "<<bcr;
00171 
00172     // fill in maps of hits; same code as addSignals, except now applied to the pileup events
00173 
00174     // EB first
00175 
00176     boost::shared_ptr<Wrapper<EBRecHitCollection>  const> EBRecHitsPTR =
00177       getProductByTag<EBRecHitCollection>(*ep, EBPileRecHitInputTag_ );
00178 
00179     if(EBRecHitsPTR ) {
00180 
00181       const EBRecHitCollection*  EBRecHits = const_cast< EBRecHitCollection * >(EBRecHitsPTR->product());
00182 
00183       LogDebug("DataMixingEMWorker") << "total # EB rechits: " << EBRecHits->size();
00184 
00185       // loop over digis, adding these to the existing maps                                         
00186       for(EBRecHitCollection::const_iterator it  = EBRecHits->begin();
00187           it != EBRecHits->end(); ++it) {
00188 
00189         EBRecHitStorage_.insert(EBRecHitMap::value_type( (it->id()), *it ));
00190 
00191 #ifdef DEBUG
00192         LogDebug("DataMixingEMWorker") << "processed EBRecHit with rawId: "
00193                                            << it->id().rawId() << "\n"
00194                                            << " rechit energy: " << it->energy();
00195 #endif
00196       }
00197     }
00198 
00199     // EE Next
00200 
00201     boost::shared_ptr<Wrapper<EERecHitCollection>  const> EERecHitsPTR =
00202       getProductByTag<EERecHitCollection>(*ep, EEPileRecHitInputTag_ );
00203 
00204     if(EERecHitsPTR ) {
00205 
00206       const EERecHitCollection*  EERecHits = const_cast< EERecHitCollection * >(EERecHitsPTR->product());
00207 
00208       LogDebug("DataMixingEMWorker") << "total # EE rechits: " << EERecHits->size();
00209 
00210       // loop over digis, adding these to the existing maps                                         
00211       for(EERecHitCollection::const_iterator it  = EERecHits->begin();
00212           it != EERecHits->end(); ++it) {
00213 
00214         EERecHitStorage_.insert(EERecHitMap::value_type( (it->id()), *it ));
00215 
00216 #ifdef DEBUG
00217         LogDebug("DataMixingEMWorker") << "processed EERecHit with rawId: "
00218                                        << it->id().rawId() << "\n"
00219                                        << " rechit energy: " << it->energy();
00220 #endif
00221       }
00222     }
00223 
00224     // ES Next
00225 
00226     boost::shared_ptr<Wrapper<ESRecHitCollection>  const> ESRecHitsPTR =
00227       getProductByTag<ESRecHitCollection>(*ep, ESPileRecHitInputTag_ );
00228 
00229     if(ESRecHitsPTR ) {
00230 
00231       const ESRecHitCollection*  ESRecHits = const_cast< ESRecHitCollection * >(ESRecHitsPTR->product());
00232 
00233       LogDebug("DataMixingEMWorker") << "total # ES rechits: " << ESRecHits->size();
00234 
00235       // loop over digis, adding these to the existing maps                                         
00236       for(ESRecHitCollection::const_iterator it  = ESRecHits->begin();
00237           it != ESRecHits->end(); ++it) {
00238 
00239         ESRecHitStorage_.insert(ESRecHitMap::value_type( (it->id()), *it ));
00240 
00241 #ifdef DEBUG
00242         LogDebug("DataMixingEMWorker") << "processed ESRecHit with rawId: "
00243                                        << it->id().rawId() << "\n"
00244                                        << " rechit energy: " << it->energy();
00245 #endif
00246       }
00247     }
00248 
00249 
00250   }
00251  
00252   void DataMixingEMWorker::putEM(edm::Event &e) {
00253 
00254     // collection of rechits to put in the event
00255     std::auto_ptr< EBRecHitCollection > EBrechits( new EBRecHitCollection );
00256     std::auto_ptr< EERecHitCollection > EErechits( new EERecHitCollection );
00257     std::auto_ptr< ESRecHitCollection > ESrechits( new ESRecHitCollection );
00258 
00259     // loop over the maps we have, re-making individual hits or digis if necessary.
00260     DetId formerID = 0;
00261     DetId currentID;
00262     float ESum = 0.;
00263     float EBTime = 0.;
00264 
00265     // EB first...
00266 
00267     EBRecHitMap::const_iterator iEBchk;
00268 
00269     for(EBRecHitMap::const_iterator iEB  = EBRecHitStorage_.begin();
00270         iEB != EBRecHitStorage_.end(); ++iEB) {
00271 
00272       currentID = iEB->first; 
00273 
00274       if (currentID == formerID) { // we have to add these rechits together
00275 
00276         ESum+=(iEB->second).energy(); 
00277       }
00278       else {
00279           if(formerID>0) {
00280             // cutoff for ESum?
00281             EcalRecHit aHit(formerID, ESum, EBTime);
00282             EBrechits->push_back( aHit );
00283           }
00284           //save pointers for next iteration
00285           formerID = currentID;
00286           ESum = (iEB->second).energy();
00287           EBTime = (iEB->second).time();  // take time of first hit in sequence - is this ok?
00288       }
00289 
00290       iEBchk = iEB;
00291       if((++iEBchk) == EBRecHitStorage_.end()) {  //make sure not to lose the last one
00292         EcalRecHit aHit(formerID, ESum, EBTime);
00293         EBrechits->push_back( aHit );     
00294       }
00295     }
00296 
00297     // EE next...
00298 
00299     // loop over the maps we have, re-making individual hits or digis if necessary.
00300     formerID = 0;
00301     ESum = 0.;
00302     float EETime = 0.;
00303     
00304     EERecHitMap::const_iterator iEEchk;
00305 
00306     for(EERecHitMap::const_iterator iEE  = EERecHitStorage_.begin();
00307         iEE != EERecHitStorage_.end(); ++iEE) {
00308 
00309       currentID = iEE->first; 
00310 
00311       if (currentID == formerID) { // we have to add these rechits together
00312 
00313         ESum+=(iEE->second).energy(); 
00314       }
00315       else {
00316           if(formerID>0) {
00317             // cutoff for ESum?
00318             EcalRecHit aHit(formerID, ESum, EETime);
00319             EErechits->push_back( aHit );
00320           }
00321           //save pointers for next iteration
00322           formerID = currentID;
00323           ESum = (iEE->second).energy();
00324           EETime = (iEE->second).time();  // take time of first hit in sequence - is this ok?
00325       }
00326 
00327       iEEchk = iEE;
00328       if((++iEEchk) == EERecHitStorage_.end()) {  //make sure not to lose the last one
00329         EcalRecHit aHit(formerID, ESum, EETime);
00330         EErechits->push_back( aHit );     
00331       }
00332     }
00333 
00334     // ES next...
00335 
00336     // loop over the maps we have, re-making individual hits or digis if necessary.
00337     formerID = 0;
00338     ESum = 0.;
00339     float ESTime = 0.;
00340 
00341     ESRecHitMap::const_iterator iESchk;
00342 
00343     for(ESRecHitMap::const_iterator iES  = ESRecHitStorage_.begin();
00344         iES != ESRecHitStorage_.end(); ++iES) {
00345 
00346       currentID = iES->first; 
00347 
00348       if (currentID == formerID) { // we have to add these rechits together
00349 
00350         ESum+=(iES->second).energy(); 
00351       }
00352       else {
00353           if(formerID>0) {
00354             // cutoff for ESum?
00355             EcalRecHit aHit(formerID, ESum, ESTime);
00356             ESrechits->push_back( aHit );
00357           }
00358           //save pointers for next iteration
00359           formerID = currentID;
00360           ESum = (iES->second).energy();
00361           ESTime = (iES->second).time();  // take time of first hit in sequence - is this ok?
00362       }
00363 
00364       iESchk = iES;
00365       if((++iESchk) == ESRecHitStorage_.end()) {  //make sure not to lose the last one
00366         EcalRecHit aHit(formerID, ESum, ESTime);
00367         ESrechits->push_back( aHit );     
00368       }
00369     }
00370 
00371     // done merging
00372 
00373     // put the collection of reconstructed hits in the event   
00374     LogInfo("DataMixingEMWorker") << "total # EB Merged rechits: " << EBrechits->size() ;
00375     LogInfo("DataMixingEMWorker") << "total # EE Merged rechits: " << EErechits->size() ;
00376     LogInfo("DataMixingEMWorker") << "total # ES Merged rechits: " << ESrechits->size() ;
00377 
00378     e.put( EBrechits, EBRecHitCollectionDM_ );
00379     e.put( EErechits, EERecHitCollectionDM_ );
00380     e.put( ESrechits, ESRecHitCollectionDM_ );
00381     
00382     // clear local storage after this event
00383 
00384     EBRecHitStorage_.clear();
00385     EERecHitStorage_.clear();
00386     ESRecHitStorage_.clear();
00387 
00388   }
00389 
00390 } //edm