DigitizerFP420.cc

Go to the documentation of this file.

// File: DigitizerFP420.cc
// Date: 12.2006
// Description: DigitizerFP420 for FP420
// Modifications:
//
// system include files
#include <memory>
 
#include "DataFormats/Common/interface/DetSetVector.h"
 
// user include files
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
 
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
 
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
#include "SimRomanPot/SimFP420/interface/DigitizerFP420.h"
 
#include "DataFormats/FP420Digi/interface/DigiCollectionFP420.h"
#include "DataFormats/FP420Digi/interface/HDigiFP420.h"
 
// needed for the geometry:
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
 
// Random Number
#include "CLHEP/Random/RandomEngine.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
 
// G4 stuff
#include "G4HCofThisEvent.hh"
#include "G4ProcessManager.hh"
#include "G4SDManager.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4TransportationManager.hh"
#include "G4UserEventAction.hh"
#include "G4VProcess.hh"
 
#include <cstdlib>
#include <vector>
 
using namespace std;
 
namespace cms {
  DigitizerFP420::DigitizerFP420(const edm::ParameterSet &conf)
      : conf_(conf), stripDigitizer_(new FP420DigiMain(conf)) {
    std::string alias(conf.getParameter<std::string>("@module_label"));
 
    //  produces<edm::DetSetVector<HDigiFP420> >().setBranchAlias( alias );
    //  produces<edm::DetSetVector<HDigiFP420SimLink> >().setBranchAlias ( alias +
    //  "hDigiFP420SimLink");
    produces<DigiCollectionFP420>().setBranchAlias(alias);
 
    trackerContainers.clear();
    trackerContainers = conf.getParameter<std::vector<std::string>>("ROUList");
 
    verbosity = conf_.getUntrackedParameter<int>("VerbosityLevel");
    dn0 = conf_.getParameter<int>("NumberFP420Detectors");
    sn0 = conf_.getParameter<int>("NumberFP420Stations");
    pn0 = conf_.getParameter<int>("NumberFP420SPlanes");
    rn0 = 7;
    // rn0 = 3;
 
    //  produces<DigiCollectionFP420>();
 
    //  produces<StripDigiCollection>();
    //   produces<HDigiFP420>();
    //  produces<edm::DetSetVector<HDigiFP420> >().setBranchAlias( alias );
 
    //  produces<DigiCollectionFP420>();
    // produces<DigiCollectionFP420>("HDigiFP420");
 
    //  produces<edm::DigiCollectionFP420>();
 
    //   produces<edm::DetSetVector<DigiCollectionFP420> >();
 
    if (verbosity > 0) {
      std::cout << "Creating a DigitizerFP420" << std::endl;
      std::cout << "DigitizerFP420: dn0=" << dn0 << " sn0=" << sn0 << " pn0=" << pn0 << " rn0=" << rn0 << std::endl;
      std::cout << "DigitizerFP420:trackerContainers.size()=" << trackerContainers.size() << std::endl;
    }
  }
 
  // Virtual destructor needed.
  DigitizerFP420::~DigitizerFP420() {
    if (verbosity > 0) {
      std::cout << "Destroying a DigitizerFP420" << std::endl;
    }
    delete stripDigitizer_;
  }
  //  void DigitizerFP420::produce(PSimHitCollection *   theCAFI,
  //  DigiCollectionFP420 & output) {
  void DigitizerFP420::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
    // be lazy and include the appropriate namespaces
    using namespace edm;
    using namespace std;
 
    if (verbosity > 0) {
      std::cout << " ===" << std::endl;
      std::cout << " ============== DigitizerFP420: start   produce= " << std::endl;
      std::cout << " ===" << std::endl;
    }
    // Get input
    //    std::cout << "DigitizerFP420 start produce" << std::endl;
    //  edm::ESHandle < ParticleDataTable > pdt;
    //  iSetup.getData( pdt );
 
    // Step A: Get Inputs for allTrackerHits
 
    edm::Handle<CrossingFrame<PSimHit>> cf_simhit;
    std::vector<const CrossingFrame<PSimHit> *> cf_simhitvec;
    for (uint32_t i = 0; i < trackerContainers.size(); i++) {
      iEvent.getByLabel("mix", trackerContainers[i], cf_simhit);
      cf_simhitvec.push_back(cf_simhit.product());
    }
    std::unique_ptr<MixCollection<PSimHit>> allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
 
    // use instead of the previous
    /*
        std::cout <<" ============== DigitizerFP420: start loop           1   "
     << std::endl; edm::Handle<CrossingFrame<PSimHit> > xFrame; std::cout <<"
     ============== DigitizerFP420: start loop           2   " << std::endl;
        iEvent.getByLabel("mix","FP420SI",xFrame);
        std::cout <<" ============== DigitizerFP420: start loop           3   "
     << std::endl; std::unique_ptr<MixCollection<PSimHit> > allTrackerHits( new
     MixCollection<PSimHit>(xFrame.product()) ); std::cout <<" ==============
     DigitizerFP420: start loop           4   " << std::endl;
  */
 
    // Loop on PSimHit
 
    // Step C: create empty output collection
    std::unique_ptr<DigiCollectionFP420> output(new DigiCollectionFP420);
    //  std::unique_ptr<edm::DetSetVector<HDigiFP420> > outputfinal(new
    //  edm::DetSetVector<HDigiFP420>(output) );
    //  std::unique_ptr<edm::DetSetVector<HDigiFP420> > outputfinal(new
    //  edm::DetSetVector<HDigiFP420>(output) );
    //  std::unique_ptr<edm::DetSetVector<HDigiFP420SimLink> > outputlink(new
    //  edm::DetSetVector<HDigiFP420SimLink>(output) );
 
    SimHitMap.clear();
 
    // ==================================
    if (verbosity > 0) {
      std::cout << " ===" << std::endl;
      std::cout << " ============== DigitizerFP420: MixCollection treatment= " << std::endl;
      std::cout << " ===" << std::endl;
    }
 
    MixCollection<PSimHit>::iterator isim;
    for (isim = allTrackerHits->begin(); isim != allTrackerHits->end(); isim++) {
      unsigned int unitID = (*isim).detUnitId();
      int det, zside, sector, zmodule;
      FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
      // below, the continues plane index should be (for even different sensor
      // index zside)
      //    unsigned int intindex = packMYIndex(rn0, pn0, sn0, det, zside, sector,
      // zmodule);
      unsigned int intindex = FP420NumberingScheme::packMYIndex(rn0, pn0, sn0, det, zside, sector, zmodule);
      //    int zScale=(rn0-1), sScale = (rn0-1)*(pn0-1), dScale =
      //(rn0-1)*(pn0-1)*(sn0-1);    unsigned int intindex = dScale*(det -
      // 1)+sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
 
      if (verbosity > 0) {
        double losenergy = (*isim).energyLoss();
        std::cout << " ===" << std::endl;
        std::cout << " ============== DigitizerFP420:  losenergy= " << losenergy << std::endl;
        std::cout << " === for intindex = " << intindex << std::endl;
      }
      // does not matter which index is used: intindex or unitID - mainly to
      // collect hits under every index
      SimHitMap[intindex].push_back((*isim));
      // for development later one( cal be used another index):
      //    SimHitMap[unitID].push_back((*isim));
    }
    //============================================================================================================================
 
    if (verbosity > 0) {
      std::cout << " ===" << std::endl;
      std::cout << " ============== DigitizerFP420: put zero to container " << std::endl;
      std::cout << " ===" << std::endl;
    }
    //    put zero to container info from the beginning (important! because not
    //    any detID is updated with coming of new event     !!!!!!
    // clean info of container from previous event
    for (int det = 1; det < dn0; det++) {
      for (int sector = 1; sector < sn0; sector++) {
        for (int zmodule = 1; zmodule < pn0; zmodule++) {
          for (int zside = 1; zside < rn0; zside++) {
            // intindex is a continues numbering of FP420
            unsigned int detID = FP420NumberingScheme::packMYIndex(rn0, pn0, sn0, det, zside, sector, zmodule);
            // int zScale=(rn0-1), sScale = (rn0-1)*(pn0-1), dScale =
            // (rn0-1)*(pn0-1)*(sn0-1); unsigned int detID = dScale*(det -
            // 1)+sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
            std::vector<HDigiFP420> collector;
            collector.clear();
            DigiCollectionFP420::Range inputRange;
            inputRange.first = collector.begin();
            inputRange.second = collector.end();
            output->putclear(inputRange, detID);
          }  // for
        }    // for
      }      // for
    }        // for
    //                                                                                                                      !!!!!!
    // if we want to keep Digi container/Collection for one event uncomment the
    // line below and vice versa
    output->clear();  // container_.clear() --> start from the beginning of the
                      // container
 
    //============================================================================================================================================
 
    if (verbosity > 0) {
      std::cout << " ===" << std::endl;
      std::cout << " ============== DigitizerFP420: start loop over det iu " << std::endl;
      std::cout << " ============== DigitizerFP420: SimHitMap.size()= " << SimHitMap.size() << std::endl;
      std::cout << " ===" << std::endl;
    }
    bool first = true;
 
    /*
  if(verbosity>0) std::cout <<"=======  DigitizerFP420:  SimHitMap size = " <<
  SimHitMap.size() << std::endl; for(unsigned int i = 0; i < SimHitMap.size();
  i++ ) {
    //    std::cout <<" ====== DigitizerFP420: i= " << i << std::endl;
    vector<PSimHit>::const_iterator simHitIter = SimHitMap[i].begin();
    vector<PSimHit>::const_iterator simHitIterEnd = SimHitMap[i].end();
    for (;simHitIter != simHitIterEnd; ++simHitIter) {
      const PSimHit ihit = *simHitIter;
      unsigned int unitID = ihit.detUnitId();
      if(verbosity>0) std::cout <<" ====== DigitizerFP420: unitID= " << unitID
  << " i=  " << i << std::endl; int det, zside, sector, zmodule;
      FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector,
  zmodule); int zScale=(rn0-1), sScale = (rn0-1)*(pn0-1), dScale =
  (rn0-1)*(pn0-1)*(sn0-1); unsigned int iu = dScale*(det - 1)+sScale*(sector -
  1)+zScale*(zmodule - 1)+zside;
    }
  }
*/
    //============================================================================================================================================
    // new:   <------
    for (unsigned int i = 0; i < SimHitMap.size(); i++) {
      vector<PSimHit>::const_iterator simHitIter = SimHitMap[i].begin();
      vector<PSimHit>::const_iterator simHitIterEnd = SimHitMap[i].end();
      for (; simHitIter != simHitIterEnd; ++simHitIter) {
        const PSimHit ihit = *simHitIter;
        unsigned int unitID = ihit.detUnitId();
        if (verbosity > 0 || verbosity == -50)
          std::cout << " ====== DigitizerFP420: unitID= " << unitID << "Hit number i=  " << i << std::endl;
        int det, zside, sector, zmodule;
        FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
 
        unsigned int iu = FP420NumberingScheme::packMYIndex(rn0, pn0, sn0, det, zside, sector, zmodule);
        if (verbosity > 0 || verbosity == -50)
          std::cout << "for Hits iu = " << iu << " sector = " << sector << " zmodule = " << zmodule
                    << " zside = " << zside << "  det=" << det << std::endl;
        //   int zScale=(rn0-1), sScale = (rn0-1)*(pn0-1), dScale =
        //   (rn0-1)*(pn0-1)*(sn0-1);
        //  unsigned int iu = dScale*(det - 1)+sScale*(sector - 1)+zScale*(zmodule
        //  - 1)+zside;
 
        if (verbosity > 0) {
          unsigned int index = FP420NumberingScheme::packFP420Index(det, zside, sector, zmodule);
          std::cout << " DigitizerFP420:       index = " << index << " iu = " << iu << std::endl;
        }
 
        //    GlobalVector bfield=pSetup->inTesla((*iu)->surface().position());
        //      CLHEP::Hep3Vector Bfieldloc=bfield();
        G4ThreeVector bfield(0., 0., 0.0);
        //  G4ThreeVector bfield(  0.5,  0.5,  1.0  );
 
        if (verbosity > 0) {
          std::cout << " ===" << std::endl;
          std::cout << " ============== DigitizerFP420:  call run for iu= " << iu << std::endl;
          std::cout << " ===" << std::endl;
        }
        collector.clear();
 
        collector = stripDigitizer_->run(SimHitMap[iu], bfield,
                                         iu);  // stripDigitizer_.run...  return
        //                       ,sScale
 
        if (verbosity > 0) {
          std::cout << " ===" << std::endl;
          std::cout << " ===" << std::endl;
          std::cout << "=======  DigitizerFP420:  collector size = " << collector.size() << std::endl;
          std::cout << " ===" << std::endl;
          std::cout << " ===" << std::endl;
        }
        /*
 
      std::vector<HDigiFP420> collector;
      collector.clear();
      DigiCollectionFP420::Range inputRange;
      inputRange.first = collector.begin();
      inputRange.second = collector.end();
      output->putclear(inputRange,detID);
      */
        if (!collector.empty()) {
          if (verbosity > 0) {
            std::cout << "         ============= DigitizerFP420:collector "
                         "start!!!!!!!!!!!!!!"
                      << std::endl;
          }
          DigiCollectionFP420::Range outputRange;
          outputRange.first = collector.begin();
          outputRange.second = collector.end();
 
          if (first) {
            // use it only if ClusterCollectionFP420 is the ClusterCollection of
            // one event, otherwise, do not use (loose 1st cl. of 1st event only)
            first = false;
            unsigned int detID0 = 0;
            output->put(outputRange, detID0);  // !!! put into adress 0 for detID
                                               // which will not be used never
          }                                    // if ( first )
 
          // put !!!
          output->put(outputRange, iu);
 
        }  // if(collector.size()>0
 
        //   }   // for
        // }   // for
        //  }   // for
        // }   // for
 
      }  // for
    }    // for
 
    // END
 
    /*
  if(verbosity>0) {
    std::vector<HDigiFP420> theAllDigis;
    theAllDigis.clear();
    DigiCollectionFP420::Range outputRange;
    DigiCollectionFP420::ContainerIterator sort_begin = outputRange.first;
    DigiCollectionFP420::ContainerIterator sort_end = outputRange.second;
    theAllDigis.insert(theAllDigis.end(), sort_begin, sort_end);
    std::cout <<"======  theAllDigis size = " << theAllDigis.size() <<
  std::endl; for (std::vector<HDigiFP420>::iterator isim = theAllDigis.begin();
         isim != theAllDigis.end(); ++isim){
      const HDigiFP420 istrip = *isim;
      std::cout <<
  "*******************************************DigitizerFP420:check1" <<
  std::endl; std::cout << " strip number=" << istrip.strip() << "  adc=" <<
  istrip.adc() << std::endl; std::cout <<" channel =" << istrip.channel() <<" V
  " << istrip.stripV() <<" VW " << istrip.stripVW() << std::endl; std::cout <<"
  ===" << std::endl; std::cout <<" ===" << std::endl; std::cout <<"
  =======================" << std::endl;
    }// for
  }
*/
    if (verbosity == -50) {
      //     check of access to the collector:
      for (int det = 1; det < dn0; det++) {
        for (int sector = 1; sector < sn0; sector++) {
          for (int zmodule = 1; zmodule < pn0; zmodule++) {
            for (int zside = 1; zside < rn0; zside++) {
              unsigned int iu = FP420NumberingScheme::packMYIndex(rn0, pn0, sn0, det, zside, sector, zmodule);
              int layer = FP420NumberingScheme::unpackLayerIndex(rn0, zside);
              int orient = FP420NumberingScheme::unpackOrientation(rn0, zside);
              std::cout << "****DigitizerFP420:check2" << std::endl;
              //    std::cout <<" iu = " << iu <<" sector = " << sector <<" zmodule
              //= " << zmodule <<" zside = " << zside << "  det=" << det <<
              // std::endl;     std::cout <<" layer = " << layer <<" orient = "
              // << orient << std::endl;
              int newdet, newzside, newsector, newzmodule;
              FP420NumberingScheme::unpackMYIndex(iu, rn0, pn0, sn0, newdet, newzside, newsector, newzmodule);
              std::cout << " newdet = " << newdet << " newsector = " << newsector << " newzmodule = " << newzmodule
                        << " newzside = " << newzside << std::endl;
 
              collector.clear();
              DigiCollectionFP420::Range outputRange;
              //    outputRange = output->get(iu);
              outputRange = output->get(iu);
 
              // fill output in collector vector (for may be sorting? or other
              // checks)
              std::vector<HDigiFP420> collector;
              //  collector.clear();
              DigiCollectionFP420::ContainerIterator sort_begin = outputRange.first;
              DigiCollectionFP420::ContainerIterator sort_end = outputRange.second;
              for (; sort_begin != sort_end; ++sort_begin) {
                collector.push_back(*sort_begin);
              }  // for
              //  std::sort(collector.begin(),collector.end());
              std::cout << " ===" << std::endl;
              std::cout << "======  collector size = " << collector.size() << std::endl;
              if (!collector.empty()) {
                std::cout << " iu = " << iu << " sector = " << sector << " zmodule = " << zmodule
                          << " zside = " << zside << "  det=" << det << " layer = " << layer << " orient = " << orient
                          << std::endl;
                std::cout << " ===" << std::endl;
              }
              vector<HDigiFP420>::const_iterator simHitIter = collector.begin();
              vector<HDigiFP420>::const_iterator simHitIterEnd = collector.end();
              for (; simHitIter != simHitIterEnd; ++simHitIter) {
                const HDigiFP420 istrip = *simHitIter;
                std::cout << " strip number=" << istrip.strip() << "  adc=" << istrip.adc() << std::endl;
                std::cout << " channel =" << istrip.channel() << " V " << istrip.stripV() << " VW " << istrip.stripVW()
                          << std::endl;
                std::cout << " ===" << std::endl;
                std::cout << " ===" << std::endl;
                std::cout << " ===================================================" << std::endl;
              }
 
              //==================================
 
            }  // for
          }    // for
        }      // for
      }        // for
 
      //     end of check of access to the strip collection
 
    }  // if(verbosity
 
    // Step D: write output to file
    if (verbosity > 0) {
      std::cout << "DigitizerFP420 recoutput" << std::endl;
    }
    // Step D: write output to file
    iEvent.put(std::move(output));
  }  // produce
 
}  // namespace cms