CMS 3D CMS Logo

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

Go to the documentation of this file.
00001 // File: DataMixingHcalWorker.cc
00002 // Description:  see DataMixingHcalWorker.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 "DataMixingHcalWorker.h"
00018 
00019 
00020 using namespace std;
00021 
00022 namespace edm
00023 {
00024 
00025   // Virtual constructor
00026 
00027   DataMixingHcalWorker::DataMixingHcalWorker() { } 
00028 
00029   // Constructor 
00030   DataMixingHcalWorker::DataMixingHcalWorker(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
00039 
00040     // Hcal 
00041 
00042     HBHErechitCollectionSig_  = ps.getParameter<edm::InputTag>("HBHEProducerSig");
00043     HOrechitCollectionSig_    = ps.getParameter<edm::InputTag>("HOProducerSig");
00044     HFrechitCollectionSig_    = ps.getParameter<edm::InputTag>("HFProducerSig");
00045     ZDCrechitCollectionSig_   = ps.getParameter<edm::InputTag>("ZDCrechitCollectionSig");
00046 
00047     HBHEPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HBHEPileRecHitInputTag");
00048     HOPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HOPileRecHitInputTag");
00049     HFPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("HFPileRecHitInputTag");
00050     ZDCPileRecHitInputTag_ = ps.getParameter<edm::InputTag>("ZDCPileRecHitInputTag");
00051 
00052     HBHERecHitCollectionDM_ = ps.getParameter<std::string>("HBHERecHitCollectionDM");
00053     HORecHitCollectionDM_   = ps.getParameter<std::string>("HORecHitCollectionDM");
00054     HFRecHitCollectionDM_   = ps.getParameter<std::string>("HFRecHitCollectionDM");
00055     ZDCRecHitCollectionDM_  = ps.getParameter<std::string>("ZDCRecHitCollectionDM");
00056 
00057 
00058   }
00059                
00060   // Virtual destructor needed.
00061   DataMixingHcalWorker::~DataMixingHcalWorker() { 
00062   }  
00063 
00064   void DataMixingHcalWorker::addHcalSignals(const edm::Event &e) { 
00065     // fill in maps of hits
00066 
00067     LogInfo("DataMixingHcalWorker")<<"===============> adding MC signals for "<<e.id();
00068 
00069     // HBHE first
00070 
00071    Handle< HBHERecHitCollection > pHBHERecHits;
00072 
00073    const HBHERecHitCollection*  HBHERecHits = 0;
00074 
00075    if( e.getByLabel( HBHErechitCollectionSig_, pHBHERecHits) ) {
00076      HBHERecHits = pHBHERecHits.product(); // get a ptr to the product
00077      LogDebug("DataMixingHcalWorker") << "total # HBHE rechits: " << HBHERecHits->size();
00078    } 
00079    
00080  
00081    if (HBHERecHits)
00082      {
00083        // loop over rechits, storing them in a map so we can add pileup later
00084        for(HBHERecHitCollection::const_iterator it  = HBHERecHits->begin();     
00085            it != HBHERecHits->end(); ++it) {
00086 
00087          HBHERecHitStorage_.insert(HBHERecHitMap::value_type( ( it->id() ), *it ));
00088          
00089 #ifdef DEBUG     
00090          LogDebug("DataMixingHcalWorker") << "processed HBHERecHit with rawId: "
00091                                       << it->id() << "\n"
00092                                       << " rechit energy: " << it->energy();
00093 #endif
00094 
00095        }
00096      }
00097 
00098    // HO next
00099 
00100    Handle< HORecHitCollection > pHORecHits;
00101 
00102    const HORecHitCollection*  HORecHits = 0;
00103 
00104    if( e.getByLabel( HOrechitCollectionSig_, pHORecHits) ){
00105      HORecHits = pHORecHits.product(); // get a ptr to the product
00106 #ifdef DEBUG
00107      LogDebug("DataMixingHcalWorker") << "total # HO rechits: " << HORecHits->size();
00108 #endif
00109    } 
00110    
00111  
00112    if (HORecHits)
00113      {
00114        // loop over rechits, storing them in a map so we can add pileup later
00115        for(HORecHitCollection::const_iterator it  = HORecHits->begin(); 
00116            it != HORecHits->end(); ++it) {
00117 
00118          HORecHitStorage_.insert(HORecHitMap::value_type( ( it->id() ), *it ));
00119          
00120 #ifdef DEBUG     
00121          LogDebug("DataMixingHcalWorker") << "processed HORecHit with rawId: "
00122                                       << it->id() << "\n"
00123                                       << " rechit energy: " << it->energy();
00124 #endif
00125 
00126        }
00127      }
00128 
00129    // HF next
00130 
00131    Handle< HFRecHitCollection > pHFRecHits;
00132 
00133    const HFRecHitCollection*  HFRecHits = 0;
00134 
00135    if( e.getByLabel( HFrechitCollectionSig_, pHFRecHits) ) {
00136      HFRecHits = pHFRecHits.product(); // get a ptr to the product
00137 #ifdef DEBUG
00138      LogDebug("DataMixingHcalWorker") << "total # HF rechits: " << HFRecHits->size();
00139 #endif
00140    } 
00141    
00142  
00143    if (HFRecHits)
00144      {
00145        // loop over rechits, storing them in a map so we can add pileup later
00146        for(HFRecHitCollection::const_iterator it  = HFRecHits->begin(); 
00147            it != HFRecHits->end(); ++it) {
00148 
00149          HFRecHitStorage_.insert(HFRecHitMap::value_type( ( it->id() ), *it ));
00150          
00151 #ifdef DEBUG     
00152          LogDebug("DataMixingHcalWorker") << "processed HFRecHit with rawId: "
00153                                       << it->id() << "\n"
00154                                       << " rechit energy: " << it->energy();
00155 #endif
00156 
00157        }
00158      }
00159 
00160    // ZDC next
00161 
00162    Handle< ZDCRecHitCollection > pZDCRecHits;
00163 
00164    const ZDCRecHitCollection*  ZDCRecHits = 0;
00165 
00166    if( e.getByLabel( ZDCrechitCollectionSig_, pZDCRecHits) ) {
00167      ZDCRecHits = pZDCRecHits.product(); // get a ptr to the product
00168 #ifdef DEBUG
00169      LogDebug("DataMixingHcalWorker") << "total # ZDC rechits: " << ZDCRecHits->size();
00170 #endif
00171    } 
00172    
00173  
00174    if (ZDCRecHits)
00175      {
00176        // loop over rechits, storing them in a map so we can add pileup later
00177        for(ZDCRecHitCollection::const_iterator it  = ZDCRecHits->begin();       
00178            it != ZDCRecHits->end(); ++it) {
00179 
00180          ZDCRecHitStorage_.insert(ZDCRecHitMap::value_type( ( it->id() ), *it ));
00181          
00182 #ifdef DEBUG     
00183          LogDebug("DataMixingHcalWorker") << "processed ZDCRecHit with rawId: "
00184                                       << it->id() << "\n"
00185                                       << " rechit energy: " << it->energy();
00186 #endif
00187 
00188        }
00189      }
00190     
00191   } // end of addEMSignals
00192 
00193   void DataMixingHcalWorker::addHcalPileups(const int bcr, const EventPrincipal *ep, unsigned int eventNr) {
00194   
00195     LogDebug("DataMixingHcalWorker") <<"\n===============> adding pileups from event  "<<ep->id()<<" for bunchcrossing "<<bcr;
00196 
00197     // fill in maps of hits; same code as addSignals, except now applied to the pileup events
00198 
00199     // HBHE first
00200 
00201     boost::shared_ptr<Wrapper<HBHERecHitCollection>  const> HBHERecHitsPTR =
00202       getProductByTag<HBHERecHitCollection>(*ep, HBHEPileRecHitInputTag_ );
00203 
00204     if(HBHERecHitsPTR ) {
00205 
00206       const HBHERecHitCollection*  HBHERecHits = const_cast< HBHERecHitCollection * >(HBHERecHitsPTR->product());
00207 
00208       LogDebug("DataMixingEMWorker") << "total # HBHE rechits: " << HBHERecHits->size();
00209 
00210       // loop over digis, adding these to the existing maps                                                     
00211       for(HBHERecHitCollection::const_iterator it  = HBHERecHits->begin();
00212           it != HBHERecHits->end(); ++it) {
00213 
00214         HBHERecHitStorage_.insert(HBHERecHitMap::value_type( (it->id()), *it ));
00215 
00216 #ifdef DEBUG
00217         LogDebug("DataMixingEMWorker") << "processed HBHERecHit with rawId: "
00218                                        << it->id().rawId() << "\n"
00219                                        << " rechit energy: " << it->energy();
00220 #endif
00221       }
00222     }
00223 
00224     // HO Next
00225 
00226     boost::shared_ptr<Wrapper<HORecHitCollection>  const> HORecHitsPTR =
00227       getProductByTag<HORecHitCollection>(*ep, HOPileRecHitInputTag_ );
00228 
00229     if(HORecHitsPTR ) {
00230 
00231       const HORecHitCollection*  HORecHits = const_cast< HORecHitCollection * >(HORecHitsPTR->product());
00232 
00233       LogDebug("DataMixingEMWorker") << "total # HO rechits: " << HORecHits->size();
00234 
00235       // loop over digis, adding these to the existing maps                                                     
00236       for(HORecHitCollection::const_iterator it  = HORecHits->begin();
00237           it != HORecHits->end(); ++it) {
00238 
00239         HORecHitStorage_.insert(HORecHitMap::value_type( (it->id()), *it ));
00240 
00241 #ifdef DEBUG
00242         LogDebug("DataMixingEMWorker") << "processed HORecHit with rawId: "
00243                                        << it->id().rawId() << "\n"
00244                                        << " rechit energy: " << it->energy();
00245 #endif
00246       }
00247     }
00248 
00249     // HF Next
00250 
00251     boost::shared_ptr<Wrapper<HFRecHitCollection>  const> HFRecHitsPTR =
00252       getProductByTag<HFRecHitCollection>(*ep, HFPileRecHitInputTag_ );
00253 
00254     if(HFRecHitsPTR ) {
00255 
00256       const HFRecHitCollection*  HFRecHits = const_cast< HFRecHitCollection * >(HFRecHitsPTR->product());
00257 
00258       LogDebug("DataMixingEMWorker") << "total # HF rechits: " << HFRecHits->size();
00259 
00260       // loop over digis, adding these to the existing maps                                                     
00261       for(HFRecHitCollection::const_iterator it  = HFRecHits->begin();
00262           it != HFRecHits->end(); ++it) {
00263 
00264         HFRecHitStorage_.insert(HFRecHitMap::value_type( (it->id()), *it ));
00265 
00266 #ifdef DEBUG
00267         LogDebug("DataMixingEMWorker") << "processed HFRecHit with rawId: "
00268                                        << it->id().rawId() << "\n"
00269                                        << " rechit energy: " << it->energy();
00270 #endif
00271       }
00272     }
00273 
00274     // ZDC Next
00275 
00276     boost::shared_ptr<Wrapper<ZDCRecHitCollection>  const> ZDCRecHitsPTR =
00277       getProductByTag<ZDCRecHitCollection>(*ep, ZDCPileRecHitInputTag_ );
00278 
00279     if(ZDCRecHitsPTR ) {
00280 
00281       const ZDCRecHitCollection*  ZDCRecHits = const_cast< ZDCRecHitCollection * >(ZDCRecHitsPTR->product());
00282 
00283       LogDebug("DataMixingEMWorker") << "total # ZDC rechits: " << ZDCRecHits->size();
00284 
00285       // loop over digis, adding these to the existing maps                                                     
00286       for(ZDCRecHitCollection::const_iterator it  = ZDCRecHits->begin();
00287           it != ZDCRecHits->end(); ++it) {
00288 
00289         ZDCRecHitStorage_.insert(ZDCRecHitMap::value_type( (it->id()), *it ));
00290 
00291 #ifdef DEBUG
00292         LogDebug("DataMixingEMWorker") << "processed ZDCRecHit with rawId: "
00293                                        << it->id().rawId() << "\n"
00294                                        << " rechit energy: " << it->energy();
00295 #endif
00296       }
00297     }
00298 
00299 
00300   }
00301  
00302   void DataMixingHcalWorker::putHcal(edm::Event &e) {
00303 
00304     // collection of rechits to put in the event
00305     std::auto_ptr< HBHERecHitCollection > HBHErechits( new HBHERecHitCollection );
00306     std::auto_ptr< HORecHitCollection > HOrechits( new HORecHitCollection );
00307     std::auto_ptr< HFRecHitCollection > HFrechits( new HFRecHitCollection );
00308     std::auto_ptr< ZDCRecHitCollection > ZDCrechits( new ZDCRecHitCollection );
00309 
00310     // loop over the maps we have, re-making individual hits or digis if necessary.
00311     DetId formerID = 0;
00312     DetId currentID;
00313     float ESum = 0.;
00314     float HBTime = 0.;
00315 
00316     // HB first...
00317 
00318     HBHERecHitMap::const_iterator iHBchk;
00319 
00320     for(HBHERecHitMap::const_iterator iHB  = HBHERecHitStorage_.begin();
00321         iHB != HBHERecHitStorage_.end(); ++iHB) {
00322 
00323       currentID = iHB->first; 
00324 
00325       if (currentID == formerID) { // we have to add these rechits together
00326 
00327         ESum+=(iHB->second).energy();  
00328 
00329       }
00330       else {
00331         if(formerID>0) {
00332           // cutoff for ESum?                                                                                 
00333           HBHERecHit aHit(formerID, ESum, HBTime);
00334           HBHErechits->push_back( aHit );
00335         }
00336         //save pointers for next iteration                                                                    
00337         formerID = currentID;
00338         ESum = (iHB->second).energy();
00339         HBTime = (iHB->second).time();  // take time of first hit in sequence - is this ok?
00340       }
00341 
00342       iHBchk = iHB;
00343       if((++iHBchk) == HBHERecHitStorage_.end()) {  //make sure not to lose the last one  
00344         HBHERecHit aHit(formerID, ESum, HBTime);
00345         HBHErechits->push_back( aHit );
00346       }
00347     }
00348 
00349     // HO next...
00350 
00351     // loop over the maps we have, re-making individual hits or digis if necessary.
00352     formerID = 0;
00353     ESum = 0.;
00354     float HOTime = 0.;
00355 
00356     HORecHitMap::const_iterator iHOchk;
00357 
00358     for(HORecHitMap::const_iterator iHO  = HORecHitStorage_.begin();
00359         iHO != HORecHitStorage_.end(); ++iHO) {
00360 
00361       currentID = iHO->first; 
00362 
00363       if (currentID == formerID) { // we have to add these rechits together
00364 
00365         ESum+=(iHO->second).energy();  
00366 
00367       }
00368       else {
00369         if(formerID>0) {
00370           // cutoff for ESum?                                                                                 
00371           HORecHit aHit(formerID, ESum, HOTime);
00372           HOrechits->push_back( aHit );
00373         }
00374         //save pointers for next iteration                                                                    
00375         formerID = currentID;
00376         ESum = (iHO->second).energy();
00377         HOTime = (iHO->second).time();  // take time of first hit in sequence - is this ok?
00378       }
00379 
00380       iHOchk = iHO;
00381       if((++iHOchk) == HORecHitStorage_.end()) {  //make sure not to lose the last one  
00382         HORecHit aHit(formerID, ESum, HOTime);
00383         HOrechits->push_back( aHit );
00384       }
00385     }
00386 
00387 
00388     // HF next...
00389 
00390     // loop over the maps we have, re-making individual hits or digis if necessary.
00391     formerID = 0;
00392     ESum = 0.;
00393     float HFTime = 0.;
00394     HFRecHit HFOldHit;
00395 
00396     HFRecHitMap::const_iterator iHFchk;
00397 
00398     for(HFRecHitMap::const_iterator iHF  = HFRecHitStorage_.begin();
00399         iHF != HFRecHitStorage_.end(); ++iHF) {
00400 
00401       currentID = iHF->first; 
00402 
00403       if (currentID == formerID) { // we have to add these rechits together
00404 
00405         ESum+=(iHF->second).energy();  
00406 
00407       }
00408       else {
00409         if(formerID>0) {
00410           // cutoff for ESum?                                                                                 
00411           HFRecHit aHit(formerID, ESum, HFTime);
00412           HFrechits->push_back( aHit );
00413         }
00414         //save pointers for next iteration                                                                    
00415         formerID = currentID;
00416         ESum = (iHF->second).energy();
00417         HFTime = (iHF->second).time();  // take time of first hit in sequence - is this ok?
00418       }
00419 
00420       iHFchk = iHF;
00421       if((++iHFchk) == HFRecHitStorage_.end()) {  //make sure not to lose the last one  
00422         HFRecHit aHit(formerID, ESum, HBTime);
00423         HFrechits->push_back( aHit );
00424       }
00425     }
00426 
00427     // ZDC next...
00428 
00429     // loop over the maps we have, re-making individual hits or digis if necessary.
00430     formerID = 0;
00431     ESum = 0.;
00432     float ZDCTime = 0.;
00433     float lowGainEnergy = 0;
00434     ZDCRecHit ZOldHit;
00435 
00436     ZDCRecHitMap::const_iterator iZDCchk;
00437 
00438     for(ZDCRecHitMap::const_iterator iZDC  = ZDCRecHitStorage_.begin();
00439         iZDC != ZDCRecHitStorage_.end(); ++iZDC) {
00440 
00441       currentID = iZDC->first; 
00442 
00443       if (currentID == formerID) { // we have to add these rechits together
00444 
00445         ESum+=(iZDC->second).energy();  
00446         
00447       }
00448       else {
00449         if(formerID>0) {
00450           // cutoff for ESum?                                                                                 
00451           ZDCRecHit aHit(formerID, ESum, ZDCTime, lowGainEnergy);
00452           ZDCrechits->push_back( aHit );
00453         }
00454         //save pointers for next iteration                                                                    
00455         formerID = currentID;
00456         ESum = (iZDC->second).energy();
00457         lowGainEnergy = (iZDC->second).lowGainEnergy();
00458         ZDCTime = (iZDC->second).time();  // take time of first hit in sequence - is this ok?
00459       }
00460       
00461       iZDCchk = iZDC;
00462       if((++iZDCchk) == ZDCRecHitStorage_.end()) {  //make sure not to lose the last one  
00463         ZDCRecHit aHit(formerID, ESum, HBTime, lowGainEnergy);
00464         ZDCrechits->push_back( aHit );
00465       }
00466     } 
00467   
00468    //done merging
00469 
00470     // put the collection of recunstructed hits in the event   
00471     LogInfo("DataMixingHcalWorker") << "total # HBHE Merged rechits: " << HBHErechits->size() ;
00472     LogInfo("DataMixingHcalWorker") << "total # HO Merged rechits: " << HOrechits->size() ;
00473     LogInfo("DataMixingHcalWorker") << "total # HF Merged rechits: " << HFrechits->size() ;
00474     LogInfo("DataMixingHcalWorker") << "total # ZDC Merged rechits: " << ZDCrechits->size() ;
00475 
00476     e.put( HBHErechits, HBHERecHitCollectionDM_ );
00477     e.put( HOrechits, HORecHitCollectionDM_ );
00478     e.put( HFrechits, HFRecHitCollectionDM_ );
00479     e.put( ZDCrechits, ZDCRecHitCollectionDM_ );
00480 
00481     // clear local storage after this event
00482     HBHERecHitStorage_.clear();
00483     HORecHitStorage_.clear();
00484     HFRecHitStorage_.clear();
00485     ZDCRecHitStorage_.clear();
00486 
00487   }
00488 
00489 } //edm