/data/git/CMSSW_5_3_11_patch5/src/SimGeneral/DataMixingModule/plugins/DataMixingEMDigiWorker.cc

Go to the documentation of this file.

// File: DataMixingEMDigiWorker.cc
// Description:  see DataMixingEMDigiWorker.h
// Author:  Mike Hildreth, University of Notre Dame
//
//--------------------------------------------

#include <map>
#include <cmath>
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/Framework/interface/ConstProductRegistry.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Provenance/interface/Provenance.h"
#include "DataFormats/Provenance/interface/BranchDescription.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
#include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
#include "CondFormats/EcalObjects/interface/EcalPedestals.h"
#include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
#include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"

//
//
#include "DataMixingEMDigiWorker.h"


using namespace std;

namespace edm
{

  // Virtual constructor

  DataMixingEMDigiWorker::DataMixingEMDigiWorker() { sel_=0;}

  // Constructor 
  DataMixingEMDigiWorker::DataMixingEMDigiWorker(const edm::ParameterSet& ps) : 
                                                            label_(ps.getParameter<std::string>("Label"))

  {                                                         

    // get the subdetector names
    //    this->getSubdetectorNames();  //something like this may be useful to check what we are supposed to do...

    // create input selector
    if (label_.size()>0){
      sel_=new Selector( ModuleLabelSelector(label_));
    }
    else {
      sel_=new Selector( MatchAllSelector());
    }

    // declare the products to produce, retrieve

    EBProducerSig_ = ps.getParameter<edm::InputTag>("EBdigiProducerSig");
    EEProducerSig_ = ps.getParameter<edm::InputTag>("EEdigiProducerSig");
    ESProducerSig_ = ps.getParameter<edm::InputTag>("ESdigiProducerSig");

    EBPileInputTag_ = ps.getParameter<edm::InputTag>("EBPileInputTag");
    EEPileInputTag_ = ps.getParameter<edm::InputTag>("EEPileInputTag");
    ESPileInputTag_ = ps.getParameter<edm::InputTag>("ESPileInputTag");

    EBDigiCollectionDM_        = ps.getParameter<std::string>("EBDigiCollectionDM");
    EEDigiCollectionDM_        = ps.getParameter<std::string>("EEDigiCollectionDM");
    ESDigiCollectionDM_        = ps.getParameter<std::string>("ESDigiCollectionDM");
   //   nMaxPrintout_            = ps.getUntrackedParameter<int>("nMaxPrintout",10);

   //EBalgo_ = new EcalDigiSimpleAlgo();
   //EEalgo_ = new EcalDigiSimpleAlgo();

   // don't think I can "produce" in a sub-class...

   //produces< EBDigiCollection >(EBDigiCollectionDM_);
   //produces< EEDigiCollection >(EEDigiCollectionDM_);

  }
               

  // Virtual destructor needed.
  DataMixingEMDigiWorker::~DataMixingEMDigiWorker() { 
    delete sel_;
    sel_=0;
  }  

  void DataMixingEMDigiWorker::addEMSignals(const edm::Event &e,const edm::EventSetup& ES) { 
    // fill in maps of hits

    LogInfo("DataMixingEMDigiWorker")<<"===============> adding MC signals for "<<e.id();

    // EB first

   Handle< EBDigiCollection > pEBDigis;

   const EBDigiCollection*  EBDigis = 0;

   if(e.getByLabel(EBProducerSig_, pEBDigis) ){
     EBDigis = pEBDigis.product(); // get a ptr to the product
     LogDebug("DataMixingEMDigiWorker") << "total # EB digis: " << EBDigis->size();
   }
   else { std::cout << "NO EBDigis! " << EBProducerSig_.label() << " " << EBdigiCollectionSig_.label() << std::endl;}
 
   if (EBDigis)
     {
       // loop over digis, storing them in a map so we can add pileup later


       for(EBDigiCollection::const_iterator it  = EBDigis->begin();     
           it != EBDigis->end(); ++it) {

         EBDigiStorage_.insert(EBDigiMap::value_type( ( it->id() ), *it ));
#ifdef DEBUG     
         LogDebug("DataMixingEMDigiWorker") << "processed EBDigi with rawId: "
                                      << it->id().rawId() << "\n"
                                      << " digi energy: " << it->energy();
#endif
       }
     }

   // EE next

   Handle< EEDigiCollection > pEEDigis;

   const EEDigiCollection*  EEDigis = 0;

   
   if(e.getByLabel(EEProducerSig_, pEEDigis) ){
     EEDigis = pEEDigis.product(); // get a ptr to the product
     LogDebug("DataMixingEMDigiWorker") << "total # EE digis: " << EEDigis->size();
   } 
   
 
   if (EEDigis)
     {
       // loop over digis, storing them in a map so we can add pileup later
       for(EEDigiCollection::const_iterator it  = EEDigis->begin();     
           it != EEDigis->end(); ++it) {

         EEDigiStorage_.insert(EEDigiMap::value_type( ( it->id() ), *it ));
#ifdef DEBUG     
         LogDebug("DataMixingEMDigiWorker") << "processed EEDigi with rawId: "
                                      << it->id().rawId() << "\n"
                                      << " digi energy: " << it->energy();
#endif

       }
     }
   // ES next

   Handle< ESDigiCollection > pESDigis;

   const ESDigiCollection*  ESDigis = 0;

   
   if(e.getByLabel( ESProducerSig_, pESDigis) ){
     ESDigis = pESDigis.product(); // get a ptr to the product
#ifdef DEBUG
     LogDebug("DataMixingEMDigiWorker") << "total # ES digis: " << ESDigis->size();
#endif
   } 

 
   if (ESDigis)
     {

       // loop over digis, storing them in a map so we can add pileup later
       for(ESDigiCollection::const_iterator it  = ESDigis->begin();     
           it != ESDigis->end(); ++it) {

         ESDigiStorage_.insert(ESDigiMap::value_type( ( it->id() ), *it ));
         
#ifdef DEBUG     
         LogDebug("DataMixingEMDigiWorker") << "processed ESDigi with rawId: "
                                      << it->id().rawId() << "\n"
                                      << " digi energy: " << it->energy();
#endif

       }
     }
    
  } // end of addEMSignals

  void DataMixingEMDigiWorker::addEMPileups(const int bcr, const EventPrincipal *ep, unsigned int eventNr, const edm::EventSetup& ES) {
  
    LogInfo("DataMixingEMDigiWorker") <<"\n===============> adding pileups from event  "<<ep->id()<<" for bunchcrossing "<<bcr;

    // fill in maps of hits; same code as addSignals, except now applied to the pileup events

    // EB first

    boost::shared_ptr<Wrapper<EBDigiCollection>  const> EBDigisPTR = 
          getProductByTag<EBDigiCollection>(*ep, EBPileInputTag_ );
 
   if(EBDigisPTR ) {

     const EBDigiCollection*  EBDigis = const_cast< EBDigiCollection * >(EBDigisPTR->product());

     LogDebug("DataMixingEMDigiWorker") << "total # EB digis: " << EBDigis->size();

       // loop over digis, adding these to the existing maps
     for(EBDigiCollection::const_iterator it  = EBDigis->begin();
         it != EBDigis->end(); ++it) {

         EBDigiStorage_.insert(EBDigiMap::value_type( (it->id()), *it ));
         
#ifdef DEBUG     
         LogDebug("DataMixingEMDigiWorker") << "processed EBDigi with rawId: "
                      << it->id().rawId() << "\n"
                      << " digi energy: " << it->energy();
#endif
      }
    }

    // EE Next

    boost::shared_ptr<Wrapper<EEDigiCollection>  const> EEDigisPTR =
          getProductByTag<EEDigiCollection>(*ep, EEPileInputTag_ );

    if(EEDigisPTR ) {

     const EEDigiCollection*  EEDigis = const_cast< EEDigiCollection * >(EEDigisPTR->product()); 

     LogDebug("DataMixingEMDigiWorker") << "total # EE digis: " << EEDigis->size();

     for(EEDigiCollection::const_iterator it  = EEDigis->begin();
         it != EEDigis->end(); ++it) {

       EEDigiStorage_.insert(EEDigiMap::value_type( (it->id()), *it ));
         
#ifdef DEBUG     
       LogDebug("DataMixingEMDigiWorker") << "processed EEDigi with rawId: "
                                      << it->id().rawId() << "\n"
                                      << " digi energy: " << it->energy();
#endif
     }
   }
    // ES Next

    boost::shared_ptr<Wrapper<ESDigiCollection>  const> ESDigisPTR =
      getProductByTag<ESDigiCollection>(*ep, ESPileInputTag_ );

    if(ESDigisPTR ) {

      const ESDigiCollection*  ESDigis = const_cast< ESDigiCollection * >(ESDigisPTR->product());

      LogDebug("DataMixingEMDigiWorker") << "total # ES digis: " << ESDigis->size();

      for(ESDigiCollection::const_iterator it  = ESDigis->begin();
          it != ESDigis->end(); ++it) {

        ESDigiStorage_.insert(ESDigiMap::value_type( (it->id()), *it ));

#ifdef DEBUG
        LogDebug("DataMixingEMDigiWorker") << "processed ESDigi with rawId: "
                                           << it->id().rawId() << "\n"
                                           << " digi energy: " << it->energy();
#endif
      }
    }
 
  
 
  }
 
  void DataMixingEMDigiWorker::putEM(edm::Event &e, const edm::EventSetup& ES) {

    // collection of digis to put in the event
    std::auto_ptr< EBDigiCollection > EBdigis( new EBDigiCollection );
    std::auto_ptr< EEDigiCollection > EEdigis( new EEDigiCollection );
    std::auto_ptr< ESDigiCollection > ESdigis( new ESDigiCollection );


    // loop over the maps we have, re-making individual hits or digis if necessary.
    DetId formerID = 0;
    DetId currentID;

    EBDataFrame EB_old;

    int gain_new = 0;
    int gain_old = 0;
    int gain_consensus = 0;
    int adc_new;
    int adc_old;
    int adc_sum;
    uint16_t data;

    // EB first...

    EBDigiMap::const_iterator iEBchk;

    

    for(EBDigiMap::const_iterator iEB  = EBDigiStorage_.begin();
        iEB != EBDigiStorage_.end(); iEB++) {

      currentID = iEB->first; 

      if (currentID == formerID) { // we have to add these digis together
        /*      
        cout<< " Adding signals " << EBDetId(currentID).ieta() << " " 
                                  << EBDetId(currentID).iphi() << std::endl;

        cout << 1 << " " ; 
        for (int i=0; i<10;++i)  std::cout << EB_old[i].adc()<< "["<<EB_old[i].gainId()<< "] " ; std::cout << std::endl;
 
        cout << 2 << " " ; 
        for (int i=0; i<10;++i)  std::cout << (iEB->second)[i].adc()<< "["<<(iEB->second)[i].gainId()<< "] " ; std::cout << std::endl;
        */
        //loop over digi samples in each DataFrame
        unsigned int sizenew = (iEB->second).size();
        unsigned int sizeold = EB_old.size();

        unsigned int max_samp = std::max(sizenew, sizeold);

        // samples from different events can be of different lengths - sum all
        // that overlap.
        // check to see if gains match - if not, scale smaller cell down.

        int sw_gain_consensus=0;


        for(unsigned int isamp = 0; isamp<max_samp; isamp++) {
          if(isamp < sizenew) {
            gain_new = (iEB->second)[isamp].gainId();
            adc_new = (iEB->second)[isamp].adc();
          }
          else { adc_new = 0;}

          if(isamp < sizeold) {
              gain_old = EB_old[isamp].gainId();
              adc_old = EB_old[isamp].adc();
          }
          else { adc_old = 0;}

          const std::vector<float> pedeStals = GetPedestals(ES,currentID);
          const std::vector<float> gainRatios = GetGainRatios(ES,currentID);

          if(adc_new>0 && adc_old>0) {
            if(gain_old == gain_new) {  // we're happy - easy case
              gain_consensus = gain_old;
            }
            else {  // lower gain sample has more energy
                              

              if(gain_old < gain_new) { // old has higher gain than new, scale to lower gain
                
        
                float ratio = gainRatios[gain_new-1]/gainRatios[gain_old-1];
                adc_old = (int) round ((adc_old - pedeStals[gain_old-1]) / ratio + pedeStals[gain_new-1] );  
                gain_consensus = gain_new;
              }
              else { // scale to old (lower) gain
                float ratio = gainRatios[gain_old-1]/gainRatios[gain_new-1];
                adc_new = (int) round ( (adc_new - pedeStals[gain_new-1]) / ratio+ pedeStals[gain_old-1] );
                gain_consensus = gain_old;
              } 
            }
          }

         
          // add values, but don't count pedestals twice
          adc_sum = adc_new + adc_old - (int) round (pedeStals[gain_consensus-1]);


          // if we are now saturating that gain, switch to the next
          if (adc_sum> 4096) {
            if (gain_consensus<3){

              double ratio = gainRatios[gain_consensus]/gainRatios[gain_consensus-1];
              adc_sum = (int) round ((adc_sum - pedeStals[gain_consensus-1])/ ratio + pedeStals[gain_consensus]  )  ;
              sw_gain_consensus=++gain_consensus;             
            }
            else adc_sum = 4096;
                
          } 

          // furthermore, make sure we don't decrease our gain once we've switched up
          // in case go back 
          if (gain_consensus<sw_gain_consensus){

            double ratio = gainRatios[sw_gain_consensus-1]/gainRatios[gain_consensus-1];
            adc_sum = (int) round((adc_sum - pedeStals[gain_consensus-1] )/ratio + pedeStals[sw_gain_consensus-1]);
            gain_consensus = sw_gain_consensus;
          }

          EcalMGPASample sample(adc_sum, gain_consensus);
          EB_old.setSample(isamp,sample);  // overwrite old sample, adding new info
        } // for sample


      } // if current = former
      else {
          if(formerID>0) {
            EBdigis->push_back( formerID, EB_old.frame().begin() );
          }
          //save pointers for next iteration
          formerID = currentID;
          EB_old = iEB->second;

      }


      iEBchk = iEB;
      if((++iEBchk) == EBDigiStorage_.end()) {  //make sure not to lose the last one
            EBdigis->push_back( currentID, (iEB->second).frame().begin()) ;
      }
    }

    // EE next...

    formerID = 0;
    EEDataFrame EE_old;

    EEDigiMap::const_iterator iEEchk;

    for(EEDigiMap::const_iterator iEE  = EEDigiStorage_.begin();
        iEE != EEDigiStorage_.end(); iEE++) {

      currentID = iEE->first; 

      if (currentID == formerID) { // we have to add these digis together

        //loop over digi samples in each DataFrame
        unsigned int sizenew = (iEE->second).size();
        unsigned int sizeold = EE_old.size();

        unsigned int max_samp = std::max(sizenew, sizeold);

        // samples from different events can be of different lengths - sum all
        // that overlap.
        // check to see if gains match - if not, scale smaller cell down.

        for(unsigned int isamp = 0; isamp<max_samp; isamp++) {
          if(isamp < sizenew) {
            gain_new = (iEE->second)[isamp].gainId();
            adc_new = (iEE->second)[isamp].adc();
          }
          else { adc_new = 0;}

          if(isamp < sizeold) {
              gain_old = EE_old[isamp].gainId();
              adc_old = EE_old[isamp].adc();
          }
          else { adc_old = 0;}

          const std::vector<float> pedeStals = GetPedestals(ES,currentID);
          const std::vector<float> gainRatios = GetGainRatios(ES,currentID);

          if(adc_new>0 && adc_old>0) {
            if(gain_old == gain_new) {  // we're happy - easy case
              gain_consensus = gain_old;
            }
            else {  // lower gain sample has more energy

              if(gain_old < gain_new) { // old has higher gain than new, scale to lower gain
                
                
                float ratio = gainRatios[gain_new-1]/gainRatios[gain_old-1];
                adc_old = (int) round ((adc_old - pedeStals[gain_old-1]) / ratio + pedeStals[gain_new-1] );  
                gain_consensus = gain_new;
              }
              else { // scale to old (lower) gain
                float ratio = gainRatios[gain_old-1]/gainRatios[gain_new-1];
                adc_new = (int) round ( (adc_new - pedeStals[gain_new-1]) / ratio+ pedeStals[gain_old-1] );
                gain_consensus = gain_old;
              } 
            }
            
          }
             

          // add values
          adc_sum = adc_new + adc_old;
          
          // if the sum saturates this gain, switch
          if (adc_sum> 4096) {
            if (gain_consensus<3){
              
              double ratio = gainRatios[gain_consensus]/gainRatios[gain_consensus-1];
              adc_sum = (int) round ((adc_sum - pedeStals[gain_consensus-1])/ ratio + pedeStals[gain_consensus]  )  ;
              ++gain_consensus;
            }
            else adc_sum = 4096;
            
          } 
          
          EcalMGPASample sample(adc_sum, gain_consensus);
          EE_old.setSample(isamp,sample);
        }

      }
      else {
          if(formerID>0) {
            EEdigis->push_back(formerID, EE_old.frame().begin() );
             
          }
          //save pointers for next iteration
          formerID = currentID;
          EE_old = iEE->second;
      }
     
      iEEchk = iEE;
      if((++iEEchk) == EEDigiStorage_.end()) {  //make sure not to lose the last one
            EEdigis->push_back(currentID, (iEE->second).frame().begin());
      }
    }
   

    // ES next...

    formerID = 0;
    ESDataFrame ES_old;

    ESDigiMap::const_iterator iESchk;

    for(ESDigiMap::const_iterator iES  = ESDigiStorage_.begin();
        iES != ESDigiStorage_.end(); iES++) {

      currentID = iES->first; 

      if (currentID == formerID) { // we have to add these digis together

        //loop over digi samples in each DataFrame
        unsigned int sizenew = (iES->second).size();
        unsigned int sizeold = ES_old.size();
        uint16_t rawdat = 0;
        unsigned int max_samp = std::max(sizenew, sizeold);

        // samples from different events can be of different lengths - sum all
        // that overlap.
        // check to see if gains match - if not, scale smaller cell down.

        for(unsigned int isamp = 0; isamp<max_samp; isamp++) {
          if(isamp < sizenew) {
            adc_new = (iES->second)[isamp].adc();
            rawdat = (iES->second)[isamp].raw();
          }
          else { adc_new = 0;}

          if(isamp < sizeold) {
              adc_old = ES_old[isamp].adc();
              rawdat = ES_old[isamp].raw();
          }
          else { adc_old = 0;}

          // add values
          adc_sum = adc_new + adc_old;
          // make data word of gain, rawdata
          adc_sum = std::min(adc_sum,4095); //first 12 bits of (uint)
          data = adc_sum+(rawdat&0xF000); // data is 14 bit word with gain as MSBs
          ES_old.setSample(isamp,data);
        }

      }
      else {
          if(formerID>0) {
            ESdigis->push_back(ES_old);
          }
          //save pointers for next iteration
          formerID = currentID;
          ES_old = iES->second;
      }

      iESchk = iES;
      if((++iESchk) == ESDigiStorage_.end()) {  //make sure not to lose the last one
        ESdigis->push_back(iES->second);
        //      ESDataFrame df( (*ESdigis)->back() );
        //for(int isamp=0; isamp<(iES->second).size(); isamp++) {
        //  df.setSample(isamp,(iES->second).data[isamp]);
        //      }
      }
    }


    // done merging

    // put the collection of reconstructed hits in the event   
    LogInfo("DataMixingEMDigiWorker") << "total # EB Merged digis: " << EBdigis->size() ;
    LogInfo("DataMixingEMDigiWorker") << "total # EE Merged digis: " << EEdigis->size() ;
    LogInfo("DataMixingEMDigiWorker") << "total # ES Merged digis: " << ESdigis->size() ;

    e.put( EBdigis, EBDigiCollectionDM_ );
    e.put( EEdigis, EEDigiCollectionDM_ );
    e.put( ESdigis, ESDigiCollectionDM_ );
    
    // clear local storage after this event

    EBDigiStorage_.clear();
    EEDigiStorage_.clear();
    ESDigiStorage_.clear();

  }
  const std::vector<float>  DataMixingEMDigiWorker::GetPedestals (const edm::EventSetup& ES, const DetId& detid) {
    
    std::vector<float> pedeStals(3);

    // get pedestals
    edm::ESHandle<EcalPedestals> pedHandle;
    ES.get<EcalPedestalsRcd>().get( pedHandle );
    
    
    const EcalPedestalsMap & pedMap = pedHandle.product()->getMap(); // map of pedestals
    EcalPedestalsMapIterator pedIter; // pedestal iterator
    EcalPedestals::Item aped; // pedestal object for a single xtal
    
    
    pedIter = pedMap.find(detid);
    if( pedIter != pedMap.end() ) {
      aped = (*pedIter);
      pedeStals[0] = aped.mean_x12;
      pedeStals[1] = aped.mean_x6;
      pedeStals[2] = aped.mean_x1;
    } else {
      edm::LogError("DataMixingMissingInput") << "Cannot find pedestals";  
      pedeStals[0] = 0;
      pedeStals[1] = 0;
      pedeStals[2] = 0;
    }
    
    
    return pedeStals;
  }

  const std::vector<float>  DataMixingEMDigiWorker::GetGainRatios(const edm::EventSetup& ES, const DetId& detid) {

    std::vector<float> gainRatios(3);
    // get gain ratios  
    edm::ESHandle<EcalGainRatios> grHandle;
    ES.get<EcalGainRatiosRcd>().get(grHandle);
    EcalMGPAGainRatio theRatio= (*grHandle)[detid];
    
    
    gainRatios[0] = 1.;
    gainRatios[1] = theRatio.gain12Over6();
    gainRatios[2] = theRatio.gain6Over1()  * theRatio.gain12Over6();

    return gainRatios;
  }


} //edm