diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/cmt/requirements b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..ed274ba39c4f84c80205e4f1bf3e6cae5a730336
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/cmt/requirements
@@ -0,0 +1,35 @@
+package MuonFastDigitization
+
+public
+use AtlasPolicy AtlasPolicy-*
+
+private
+use AtlasCLHEP AtlasCLHEP-* External
+use AtlasROOT AtlasROOT-* External
+use GaudiInterface GaudiInterface-* External
+
+use StoreGate StoreGate-* Control
+use AthenaKernel AthenaKernel-* Control
+use AthenaBaseComps AthenaBaseComps-* Control
+
+use Identifier Identifier-* DetectorDescription
+
+use MuonIdHelpers MuonIdHelpers-* MuonSpectrometer
+use MuonSimEvent MuonSimEvent-* MuonSpectrometer
+use MuonSimData MuonSimData-* MuonSpectrometer
+use MuonReadoutGeometry MuonReadoutGeometry-* MuonSpectrometer/MuonDetDescr
+use MuonPrepRawData MuonPrepRawData-* MuonSpectrometer/MuonReconstruction/MuonRecEvent
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+use MuonRecToolInterfaces MuonRecToolInterfaces-* MuonSpectrometer/MuonReconstruction/MuonRecTools
+use PathResolver PathResolver-* Tools
+
+
+public
+library MuonFastDigitization *.cxx components/*.cxx
+apply_pattern component_library
+
+#private
+#macro cppdebugflags '$(cppdebugflags_s)'
+#macro_remove componentshr_linkopts "-Wl,-s"
+
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4b63f026026fc0bd44a849b816e01e49a06ab179
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx
@@ -0,0 +1,552 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//Gaudi - Core
+#include "StoreGate/StoreGateSvc.h"
+
+//Inputs
+#include "MuonSimData/MuonSimDataCollection.h"
+#include "MuonSimData/MuonSimData.h"
+
+//Outputs
+#include "MuonPrepRawData/MMPrepDataContainer.h"
+#include "MuonPrepRawData/MMPrepData.h"
+
+#include "MM_FastDigitizer.h"
+#include "MuonSimEvent/MM_SimIdToOfflineId.h"
+#include "MuonSimEvent/GenericMuonSimHitCollection.h"
+#include "MuonSimEvent/MicromegasHitIdHelper.h"
+
+#include "MuonReadoutGeometry/NSWenumeration.h"
+#include "MuonReadoutGeometry/NSWgeometry.h"
+#include "MuonReadoutGeometry/MuonDetectorManager.h"
+#include "MuonReadoutGeometry/MMReadoutElement.h"
+#include "MuonIdHelpers/MmIdHelper.h"
+#include "TrkEventPrimitives/LocalDirection.h"
+
+#include "Identifier/Identifier.h"
+#include "TH1.h"
+#include "TTree.h"
+#include "TFile.h"
+
+//Random Numbers
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandGauss.h"
+
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+
+using namespace Muon;
+/*******************************************************************************/
+ MM_FastDigitizer::MM_FastDigitizer(const std::string& name, ISvcLocator* pSvcLocator)
+: AthAlgorithm(name, pSvcLocator)
+ , m_idHelperTool("Muon::MuonIdHelperTool/MuonIdHelperTool" )
+ , m_muonClusterCreator("Muon::MuonClusterOnTrackCreator/MuonClusterOnTrackCreator")
+ , m_rndmSvc("AtRndmGenSvc", name )
+ , m_inputObjectName("MicromegasSensitiveDetector")
+{
+ declareProperty("InputObjectName", m_inputObjectName = "MicromegasSensitiveDetector", "name of the input object");
+ declareProperty("RndmEngine", m_rndmEngineName, "Random engine name");
+ declareProperty("RndmSvc", m_rndmSvc, "Random Number Service used in Muon digitization");
+ declareProperty("UseTimeShift", m_useTimeShift = true, "Use time shift");
+ declareProperty("EnergyThreshold", m_energyThreshold = 50, "Minimal energy to produce a PRD" );
+ declareProperty("CheckIds", m_checkIds = false, "Turn on validity checking of Identifiers" );
+ declareProperty("MaskMultiplet", m_maskMultiplet = 0, "0: all, 1: first, 2: second, 3: both" );
+}
+/*******************************************************************************/
+MM_FastDigitizer::~MM_FastDigitizer() {
+}
+/*******************************************************************************/
+StatusCode MM_FastDigitizer::initialize() {
+ StatusCode status(StatusCode::SUCCESS);
+
+ // initialize transient event store
+ ATH_MSG_DEBUG ( "Retrieved StoreGateSvc." );
+
+ status = service("ActiveStoreSvc", m_activeStore);
+ if ( !status.isSuccess() ) {
+ ATH_MSG_FATAL ( "Could not get active store service" );
+ return status;
+ }
+
+ if( m_muonClusterCreator.retrieve().isFailure() ){
+ ATH_MSG_WARNING("Failed to retrieve " << m_muonClusterCreator);
+ return StatusCode::FAILURE;
+ }
+
+
+ // initialize transient detector store and MuonGeoModel OR MuonDetDescrManager
+ StoreGateSvc* detStore=0;
+ m_detManager=0;
+ status = serviceLocator()->service("DetectorStore", detStore);
+ if (status.isSuccess()) {
+ if(detStore->contains<MuonGM::MuonDetectorManager>( "Muon" )){
+
+ status = detStore->retrieve(m_detManager);
+ if (status.isFailure()) {
+ ATH_MSG_FATAL ( "Could not retrieve MuonGeoModelDetectorManager!" );
+ return status;
+ }
+ else {
+ ATH_MSG_DEBUG ( "Retrieved MuonGeoModelDetectorManager from StoreGate" );
+ //initialize the MmIdHelper
+ m_idHelper = m_detManager->mmIdHelper();
+ if(!m_idHelper) return status;
+ ATH_MSG_DEBUG ( "Retrieved MmIdHelper " << m_idHelper );
+ }
+ }
+ }
+ else {
+ ATH_MSG_FATAL ( "Could not retrieve DetectorStore!" );
+ return status;
+ }
+
+ if( m_idHelperTool.retrieve().isFailure() ){
+ ATH_MSG_WARNING("Failed to retrieve " << m_idHelperTool);
+ return StatusCode::FAILURE;
+ }
+
+ // check the input object name
+ if (m_inputObjectName=="") {
+ ATH_MSG_FATAL ( "Property InputObjectName not set !" );
+ return StatusCode::FAILURE;
+ }
+ else {
+ ATH_MSG_DEBUG ( "Input objects: '" << m_inputObjectName << "'" );
+ }
+
+ // getting our random numbers stream
+ ATH_MSG_DEBUG ( "Getting random number engine : <" << m_rndmEngineName << ">" );
+ if (!m_rndmSvc.retrieve().isSuccess()){
+ ATH_MSG_ERROR("Could not initialize Random Number Service");
+ }
+ m_rndmEngine = m_rndmSvc->GetEngine(m_rndmEngineName);
+ if (m_rndmEngine==0) {
+ ATH_MSG_ERROR("Could not find RndmEngine : " << m_rndmEngineName );
+ return StatusCode::FAILURE;
+ }
+
+
+ m_file = new TFile("MM_plots.root","RECREATE");
+ m_ntuple = new TTree("a","a");
+ m_ntuple->Branch("slx",&slx);
+ m_ntuple->Branch("sly",&sly);
+ m_ntuple->Branch("slz",&slz);
+ m_ntuple->Branch("dlx",&dlx);
+ m_ntuple->Branch("dly",&dly);
+ m_ntuple->Branch("dlz",&dlz);
+ m_ntuple->Branch("sulx",&sulx);
+ m_ntuple->Branch("suly",&suly);
+ m_ntuple->Branch("tsulx",&tsulx);
+ m_ntuple->Branch("tsuly",&tsuly);
+ m_ntuple->Branch("tsulz",&tsulz);
+ m_ntuple->Branch("stsulx",&stsulx);
+ m_ntuple->Branch("stsuly",&stsuly);
+ m_ntuple->Branch("stsulz",&stsulz);
+ m_ntuple->Branch("ang",&ang);
+ m_ntuple->Branch("shift",&shift);
+ m_ntuple->Branch("resx",&resx);
+ m_ntuple->Branch("resy",&resy);
+ m_ntuple->Branch("suresx",&suresx);
+ m_ntuple->Branch("suresy",&suresy);
+ m_ntuple->Branch("resz",&resz);
+ m_ntuple->Branch("err",&err);
+ m_ntuple->Branch("res",&res);
+ m_ntuple->Branch("pull",&pull);
+ m_ntuple->Branch("is",&is);
+ m_ntuple->Branch("seta",&seta);
+ m_ntuple->Branch("sphi",&sphi);
+ m_ntuple->Branch("sml",&sml);
+ m_ntuple->Branch("sl",&sl);
+ m_ntuple->Branch("ss",&ss);
+ m_ntuple->Branch("ieta",&ieta);
+ m_ntuple->Branch("iphi",&iphi);
+ m_ntuple->Branch("iml",&iml);
+ m_ntuple->Branch("il",&il);
+ m_ntuple->Branch("ich",&ich);
+ m_ntuple->Branch("istr",&istr);
+ m_ntuple->Branch("exitcode",&exitcode);
+ m_ntuple->Branch("mode",&mode);
+ m_ntuple->Branch("pdg",&pdg);
+ m_ntuple->Branch("trkid",&trkid);
+ m_ntuple->Branch("gpx",&gpx);
+ m_ntuple->Branch("gpy",&gpy);
+ m_ntuple->Branch("gpz",&gpz);
+ m_ntuple->Branch("gpr",&gpr);
+ m_ntuple->Branch("gpp",&gpp);
+ m_ntuple->Branch("dgpx",&dgpx);
+ m_ntuple->Branch("dgpy",&dgpy);
+ m_ntuple->Branch("dgpz",&dgpz);
+ m_ntuple->Branch("dgpr",&dgpr);
+ m_ntuple->Branch("dgpp",&dgpp);
+ m_ntuple->Branch("globalHitTime", &globalHitTime);
+ m_ntuple->Branch("tofCorrection", &tofCorrection);
+ m_ntuple->Branch("bunchTime", &bunchTime);
+ m_ntuple->Branch("e", &e);
+ m_ntuple->Branch("edep", &edep);
+ return StatusCode::SUCCESS;
+
+}
+/*******************************************************************************/
+StatusCode MM_FastDigitizer::execute() {
+
+
+
+
+ // Create and record the SDO container in StoreGate
+ MuonSimDataCollection* sdoContainer = new MuonSimDataCollection();
+ if (evtStore()->record(sdoContainer,"MM_SDO").isFailure()) {
+ ATH_MSG_ERROR ( "Unable to record MM SDO collection in StoreGate" );
+ return StatusCode::FAILURE;
+ }
+
+ MMPrepDataContainer* prdContainer = new MMPrepDataContainer(m_idHelper->detectorElement_hash_max());
+ std::string key = "MM_Measurements";
+ ATH_MSG_DEBUG(" Done! Total number of of MM chambers with PRDS: " << prdContainer->numberOfCollections() << " key " << key);
+ if (evtStore()->record(prdContainer,key).isFailure()) {
+ ATH_MSG_ERROR ( "Unable to record MMPrepData container in StoreGate" );
+ return StatusCode::FAILURE;
+ }
+
+ if( m_maskMultiplet == 3 ) {
+
+ return StatusCode::SUCCESS;
+ }
+ // as the MMPrepDataContainer only allows const accesss, need a local vector as well.
+ std::vector<MMPrepDataCollection*> localMMVec(m_idHelper->detectorElement_hash_max());
+
+ const DataHandle< GenericMuonSimHitCollection > collGMSH;
+ if ( evtStore()->retrieve( collGMSH,m_inputObjectName ).isFailure()) {
+ ATH_MSG_WARNING("No MM hits found in SG");
+ return StatusCode::FAILURE;
+ }
+
+ ATH_MSG_DEBUG( "Retrieved " << collGMSH->size() << " MM hits!");
+
+ // Get the MicroMegas Id hit helper
+ MM_SimIdToOfflineId simToOffline(*m_idHelper);
+
+ std::map<Identifier,int> hitsPerChannel;
+ int nhits = 0;
+
+ GenericMuonSimHitCollection::const_iterator iterMM;
+
+ const GenericMuonSimHit* previousHit = 0;
+
+ for (iterMM=collGMSH->begin();iterMM!=collGMSH->end();++iterMM) {
+ const GenericMuonSimHit& hit = *iterMM;
+
+ // SimHits without energy loss are not recorded.
+ // not needed because of already done in sensitive detector
+ // https://svnweb.cern.ch/trac/atlasoff/browser/MuonSpectrometer/MuonG4/MuonG4SD/trunk/src/MicromegasSensitiveDetector.cxx?rev=542333#L65
+ // if(hit.depositEnergy()==0.) continue;
+
+ if( previousHit && abs(hit.particleEncoding())==13 && abs(previousHit->particleEncoding())==13 ) {
+ Amg::Vector3D diff = previousHit->localPosition() - hit.localPrePosition();
+ ATH_MSG_VERBOSE("second hit from a muon: prev " << previousHit->localPosition() << " current " << hit.localPrePosition()
+ << " diff " << diff );
+ if( diff.mag() < 0.1 ) continue;
+ }
+
+ globalHitTime = hit.globalpreTime();
+ tofCorrection = hit.globalPosition().mag()/CLHEP::c_light;
+ bunchTime = globalHitTime - tofCorrection;
+ const float stripPropagationTime = 0.;
+ static const float jitter = 0.;
+ float MMDigitTime = bunchTime + jitter + stripPropagationTime;
+
+ float timeWindowStrip = 120.; //driftvelocity gap;
+ if (MMDigitTime < -timeWindowStrip || MMDigitTime > timeWindowStrip) {
+ exitcode = 4;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // Read the information about the Micro Megas hit
+ ATH_MSG_VERBOSE("MM hit: r " << hit.globalPosition().perp() << " z " << hit.globalPosition().z() << " mclink " << hit.particleLink() );
+
+ // convert simHit id to offline layer id; make sanity checks; retrieve the associated detector element.
+ int simId = hit.GenericId();
+ Identifier layid = simToOffline.convert(simId);
+
+ // sanity checks
+ if( !m_idHelper->is_mm(layid) ){
+ ATH_MSG_WARNING("MM id is not a mm id! " << m_idHelperTool->toString(layid));
+ }
+ if( !m_idHelper->is_muon(layid) ){
+ ATH_MSG_WARNING("MM id is not a muon id! " << m_idHelperTool->toString(layid));
+ }
+ if( m_idHelper->is_mdt(layid)||m_idHelper->is_rpc(layid)||m_idHelper->is_tgc(layid)||m_idHelper->is_csc(layid)||m_idHelper->is_stgc(layid) ){
+ ATH_MSG_WARNING("MM id has wrong technology type! " << m_idHelper->is_mdt(layid) << " " << m_idHelper->is_rpc(layid)
+ << " " << m_idHelper->is_tgc(layid) << " " << m_idHelper->is_csc(layid) << " " << m_idHelper->is_stgc(layid) );
+ }
+
+ // remove hits in masked multiplet
+ if( m_maskMultiplet == m_idHelperTool->stgcIdHelper().multilayer(layid) ) continue;
+
+ // get readout element
+ const MuonGM::MMReadoutElement* detEl = m_detManager->getMMReadoutElement(layid);
+ if( !detEl ){
+ ATH_MSG_WARNING("Failed to retrieve detector element for: " << m_idHelperTool->toString(layid) );
+ continue;
+ }
+
+ // collections stored per readout element
+ IdentifierHash hash;
+ m_idHelper->get_detectorElement_hash(layid, hash);
+ MuonPrepDataCollection<Muon::MMPrepData>* col = localMMVec[hash];
+ if( !col ){
+ col = new MMPrepDataCollection(hash);
+ col->setIdentifier(m_idHelper->channelID(m_idHelper->parentID(layid), m_idHelper->multilayer(layid),1,1) );
+ if( prdContainer->addCollection(col,hash).isFailure() ){
+ ATH_MSG_WARNING("Failed to add collection with hash " << (int)hash );
+ delete col;col=0;
+ continue;
+ }
+ localMMVec[hash] = col;
+ }
+
+ // surface
+ const Trk::PlaneSurface& surf = detEl->surface(layid);
+
+ //Angle
+ const Amg::Vector3D GloDire(hit.globalDirection().x(), hit.globalDirection().y(), hit.globalDirection().z());
+ Trk::LocalDirection locDire;
+ surf.globalToLocalDirection(GloDire, locDire);
+ float inAngle_XZ = fabs( locDire.angleXZ() / CLHEP::degree);
+ inAngle_XZ = 90. - inAngle_XZ ;
+
+ // compute hit position within the detector element/surfaces
+ Amg::Transform3D gToL = detEl->absTransform().inverse();
+ Amg::Vector3D hpos(hit.globalPosition().x(),hit.globalPosition().y(),hit.globalPosition().z());
+ Amg::Vector3D lpos = gToL*hpos;
+ Amg::Vector3D slpos = surf.transform().inverse()*hpos;
+
+ Amg::Vector3D ldir;
+ if (MM_READOUT[m_idHelper->gasGap(layid)-1]==1)
+ ldir = surf.transform().inverse().linear()*Amg::Vector3D(hit.globalDirection().x(), hit.globalDirection().y(), hit.globalDirection().z());
+ else
+ ldir = surf.transform().inverse().linear()*Amg::Vector3D(hit.globalDirection().x(), hit.globalDirection().y(), -hit.globalDirection().z());
+
+ double scale;//, scaletop;
+// double gasgap = 5.;
+
+ scale = -slpos.z()/ldir.z();
+// scaletop = (gasgap+slpos.z())/ldir.z();
+
+ Amg::Vector3D hitOnSurface = slpos + scale*ldir;
+// Amg::Vector3D hitOnTopSurface = slpos + scaletop*ldir;
+
+ // account for time offset
+ const double vdrift = 0.047;
+ double shiftTimeOffset = MMDigitTime*vdrift;
+ Amg::Vector3D hitAfterTimeShift(hitOnSurface.x(),hitOnSurface.y(),shiftTimeOffset);
+ Amg::Vector3D hitAfterTimeShiftOnSurface = hitAfterTimeShift - (shiftTimeOffset/ldir.z())*ldir;
+
+ // resolution = -.01/3 * angle + .64/3.
+ double resolution;
+ if (fabs(inAngle_XZ)<3)
+ resolution = .07;
+ else
+ resolution = ( -.001/3.*fabs(inAngle_XZ) ) + .28/3.;
+ double sp = CLHEP::RandGauss::shoot(m_rndmEngine, 0, resolution);
+
+ ATH_MSG_VERBOSE("slpos.z " << slpos.z() << ", ldir " << ldir.z() << ", scale " << scale << ", hitOnSurface.z " << hitOnSurface.z() );
+
+ if( fabs(hitOnSurface.z()) > 0.1 ) ATH_MSG_WARNING("bad propagation to surface " << hitOnSurface );
+ if( fabs(hitAfterTimeShiftOnSurface.z()) > 0.1 ) ATH_MSG_WARNING("bad propagation to surface after time shift " << hitAfterTimeShiftOnSurface );
+ // smeared local position
+ Amg::Vector2D posOnSurf(hitOnSurface.x()+sp,hitOnSurface.y());
+
+// int digiMode = 0;
+ // for large angles project perpendicular to surface
+ if( fabs(inAngle_XZ) > 70 ){
+ posOnSurf[0]=(slpos.x()+sp);
+// digiMode = 1;
+ // if using timing information use hit position after shift
+ }else if( m_useTimeShift ){
+ posOnSurf[0]=(hitAfterTimeShiftOnSurface.x()+sp);
+// digiMode = 2;
+ }
+
+ ////// fill first part of ntuple
+ Amg::Vector3D repos = detEl->globalPosition();
+ MicromegasHitIdHelper* hitHelper = MicromegasHitIdHelper::GetHelper();
+
+ std::string stName = m_idHelper->stationNameString(m_idHelper->stationName(layid));
+ int isSmall = stName[2] == 'S';
+ slx = hit.localPosition().x();
+ sly = hit.localPosition().y();
+ slz = hit.localPosition().z();
+ dlx = lpos.x();
+ dly = lpos.y();
+ sulx = posOnSurf.x();
+ suly = posOnSurf.y();
+ tsulx = hitOnSurface.x();
+ tsuly = hitOnSurface.y();
+ tsulz = hitOnSurface.z();
+ stsulx = hitAfterTimeShiftOnSurface.x();
+ stsuly = hitAfterTimeShiftOnSurface.y();
+ stsulz = hitAfterTimeShiftOnSurface.z();
+ ang = inAngle_XZ;
+ shift = shiftTimeOffset;
+ resx = hit.localPosition().x() - lpos.x();
+ resy = hit.localPosition().y() - lpos.y();
+ resz = hit.localPosition().z() - lpos.z();
+ suresx = posOnSurf.x()-hitOnSurface.x();
+ suresy = posOnSurf.y()-hitOnSurface.y();
+ err = -99999.;
+ res = -99999.;
+ pull = -99999.;
+ is = isSmall;
+ seta = hitHelper->GetZSector(simId);
+ sphi = hitHelper->GetPhiSector(simId);
+ sml = hitHelper->GetMultiLayer(simId);
+ sl = hitHelper->GetLayer(simId);
+ ss = hitHelper->GetSide(simId);
+ ieta = m_idHelper->stationEta(layid);
+ iphi = m_idHelper->stationPhi(layid);
+ iml = m_idHelper->multilayer(layid);
+ il = m_idHelper->gasGap(layid);
+ ich = -99999;
+ istr = -99999;
+ exitcode = 0;
+ mode = 0;
+ pdg = hit.particleEncoding();
+ trkid = hit.trackNumber();
+ gpx = hit.globalPosition().x();
+ gpy = hit.globalPosition().y();
+ gpz = hit.globalPosition().z();
+ gpr = hit.globalPosition().perp();
+ gpp = hit.globalPosition().phi();
+ dgpx = repos.x();
+ dgpy = repos.y();
+ dgpz = repos.z();
+ dgpr = repos.perp();
+ dgpp = repos.phi();
+ edep = hit.depositEnergy();
+ e = hit.kineticEnergy();
+
+
+ /////
+ if(hit.kineticEnergy()<m_energyThreshold && abs(hit.particleEncoding())==11) {
+ exitcode = 5;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // perform bound check
+ if( !surf.insideBounds(posOnSurf) ){
+ exitcode = 1;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // recalculate strip number after smearing
+ int stripNumber = detEl->stripNumber(posOnSurf,layid);
+ //int LastStripNumber = detEl->stripNumber(posOnTopSurf, layid);
+
+ if( stripNumber == -1 ){
+ ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperTool->toString(layid) );
+ ATH_MSG_WARNING(" pos " << posOnSurf << " z " << slpos.z() );
+ exitcode = 2;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // Recalculate the Identifier using the strip number
+ // here need to use parent ID to avoid creating corrupted identifiers
+ Identifier parentID = m_idHelper->parentID(layid);
+ Identifier id = m_idHelper->channelID(parentID, m_idHelper->multilayer(layid), m_idHelper->gasGap(layid),stripNumber,m_checkIds);
+
+ // only one hit per channel
+ int& counts = hitsPerChannel[id];
+ ++counts;
+ if( counts > 1 ) continue;
+ ++nhits;
+
+
+ if( stripNumber >= detEl->numberOfStrips(layid) ){
+ ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperTool->toString(layid) );
+ ATH_MSG_WARNING(" pos " << posOnSurf << " z " << slpos.z() );
+ exitcode = 3;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // Recalculate the Identifier using the strip number
+ // here need to use parent ID to avoid creating corrupted identifiers
+ id = m_idHelper->channelID(parentID, m_idHelper->multilayer(layid), m_idHelper->gasGap(layid),stripNumber,m_checkIds);
+
+ if( stripNumber != m_idHelper->channel(id) ) {
+ ATH_MSG_WARNING(" bad stripNumber: in " << stripNumber << " from id " << m_idHelper->channel(id));
+ exitcode = 4;
+ continue;
+ }
+
+ std::vector<Identifier> rdoList;
+ rdoList.push_back(id);
+ Amg::MatrixX cov(1,1); cov.setIdentity();
+ cov *= resolution*resolution;
+ MMPrepData* prd = new MMPrepData( id,hash,posOnSurf,rdoList,&cov,detEl);
+ prd->setHashAndIndex(col->identifyHash(), col->size()); // <<< add this line to the MM code as well
+ col->push_back(prd);
+ // ATH_MSG_VERBOSE("Global hit: r " << hit.globalPosition().perp() << " phi " << hit.globalPosition().phi() << " z " << hit.globalPosition().z()
+ // << " detEl: r " << repos.perp() << " phi " << repos.phi() << " z " << repos.z() << " " << m_idHelper->print_to_string(id) );
+
+ // ATH_MSG_VERBOSE("Local hit: x " << hit.localPosition().x() << " y " << hit.localPosition().y() << " z " << hit.localPosition().z()
+ // << " detEl: x " << lpos.x() << " y " << lpos.y() << " z " << lpos.z() << " strip " << stripNumber);
+ ATH_MSG_DEBUG(" hit: " << m_idHelperTool->toString(id) << " hitx " << posOnSurf.x() << " residual " << posOnSurf.x() - hitOnSurface.x()
+ << " pull " << (posOnSurf.x() - hitOnSurface.x())/resolution );
+ err = resolution;
+ ich = m_idHelper->channel(id);
+ istr = stripNumber;
+
+// const MuonClusterOnTrack* rot = m_muonClusterCreator->createRIO_OnTrack( *prd, hit.globalPosition() );
+// if( rot ){
+// res = rot->localParameters().get(Trk::locX)-hitOnSurface.x();
+// pull = res/rot->localErrorMatrix().error(Trk::locX);
+// delete rot;
+// }
+
+ m_ntuple->Fill();
+ // create SDO
+ MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(hitOnSurface.x(),hitOnSurface.y()));
+ //Record the SDO collection in StoreGate
+ std::vector<MuonSimData::Deposit> deposits;
+ deposits.push_back(deposit);
+ sdoContainer->insert ( std::make_pair ( id, MuonSimData(deposits,0) ) );
+ // OLD CODE ENDS HERE
+
+ previousHit = &hit;
+ }
+
+ if( msgLvl(MSG::DEBUG) ){
+ std::map<Identifier,int>::const_iterator hit = hitsPerChannel.begin();
+ std::map<Identifier,int>::const_iterator hit_end = hitsPerChannel.end();
+ msg(MSG::DEBUG) << " number of channels with hit " << hitsPerChannel.size() << " nhits " << nhits << std::endl;
+ for( ;hit!=hit_end;++hit ){
+ msg(MSG::DEBUG) << " " << m_idHelperTool->toString(hit->first) << " -> " << hit->second << std::endl;
+ }
+ msg(MSG::DEBUG) << endreq;
+ }
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+StatusCode MM_FastDigitizer::finalize() {
+ m_ntuple->Write();
+ m_file->Write();
+ m_file->Close();
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+/** Function to convert Radians into Degrees */
+float MM_FastDigitizer::RadsToDegrees(float Radians)
+{
+ float Degrees = Radians * (180.) / M_PI;
+ return Degrees;
+}
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..8bdfbdf7ee0eeefc99481bda6dfd6d17da316ab3
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h
@@ -0,0 +1,134 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONDIGITIZATION_MM_FASTDIGITIZER_H
+#define MUONDIGITIZATION_MM_FASTDIGITIZER_H
+
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "StoreGate/StoreGateSvc.h"
+
+//Random
+#include "CLHEP/Random/RandomEngine.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandGauss.h"
+
+#include "MuonIdHelpers/MuonIdHelperTool.h"
+#include "MuonRecToolInterfaces/IMuonClusterOnTrackCreator.h"
+
+class TTree;
+class TFile;
+
+class MmIdHelper;
+namespace MuonGM {
+ class MuonDetectorManager;
+}
+
+//Random
+namespace CLHEP{
+ class HepRandomEngine;
+}
+
+
+class IAtRndmGenSvc;
+class ActiveStoreSvc;
+
+class MM_FastDigitizer : public AthAlgorithm {
+
+ public:
+
+ MM_FastDigitizer(const std::string& name, ISvcLocator* pSvcLocator);
+ ~MM_FastDigitizer();
+
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+ ServiceHandle<IAtRndmGenSvc> getRndmSvc() const { return m_rndmSvc; } // Random number service
+ CLHEP::HepRandomEngine *getRndmEngine() const { return m_rndmEngine; } // Random number engine used
+
+ float RadsToDegrees(float Radians);
+
+ private:
+
+ ActiveStoreSvc* m_activeStore;
+
+ const MuonGM::MuonDetectorManager* m_detManager;
+ const MmIdHelper* m_idHelper;
+
+ TFile* m_file;
+ TTree* m_ntuple;
+ float slx; // local position simhit in G4
+ float sly;
+ float slz;
+ float dlx; // local position simhit in GeoModel frame
+ float dly;
+ float dlz;
+ float sulx; // local position from REin tracking surface frame
+ float suly;
+ float tsulx; // local position simhit in tracking surface frame
+ float tsuly;
+ float tsulz;
+ float stsulx; // local position simhit in tracking surface frame
+ float stsuly;
+ float stsulz;
+ float ang; // local angel
+ float shift; // shift due to bunch offset
+ float resx; // residuals in local G4 - GeoModel
+ float resy;
+ float resz;
+ float suresx; // residuals g4 - RE stripPosition in tracking frame
+ float suresy;
+ float err; // error
+ float res; // residual + pull
+ float pull;
+ int is; // simulation Identifier: issmall
+ int seta; // eta
+ int sphi; // phi
+ int sml; // multi layer
+ int sl; // layer
+ int ss; // strip
+ int ieta; // offline id fields
+ int iphi;
+ int iml;
+ int il;
+ int ich; // strip number of Id
+ int istr; // strip number from RE
+ int exitcode; // flag reason why prd not made
+ int mode; // flag digitization strategy
+ int pdg; // pdg ID
+ int trkid; // track number in G4
+ float gpx; // global position of the simhit
+ float gpy;
+ float gpz;
+ float gpr; // radial position pos.perp()
+ float gpp; // phi position pos.phi()
+ float dgpx; // same for readout element center
+ float dgpy;
+ float dgpz;
+ float dgpr;
+ float dgpp;
+ float tofCorrection;
+ float bunchTime;
+ float globalHitTime;
+ float e;
+ float edep;
+
+ protected:
+ ToolHandle <Muon::MuonIdHelperTool> m_idHelperTool;
+ ToolHandle <Muon::IMuonClusterOnTrackCreator> m_muonClusterCreator;
+ ServiceHandle <IAtRndmGenSvc> m_rndmSvc; // Random number service
+ CLHEP::HepRandomEngine *m_rndmEngine; // Random number engine used - not init in SiDigitization
+ std::string m_rndmEngineName;// name of random engine
+ std::string m_inputObjectName; // name of the input objects
+ bool m_useTimeShift;
+ double m_energyThreshold;
+ bool m_checkIds;
+ int m_maskMultiplet;
+};
+
+#endif // MUONDIGITIZATION_MM_DIGITIZER_H
+
+
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_entries.cxx b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2fba763438624bbb56fbdd46af7d15fa38b54f52
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_entries.cxx
@@ -0,0 +1,12 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "../sTgcFastDigitizer.h"
+#include "../MM_FastDigitizer.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+DECLARE_ALGORITHM_FACTORY( sTgcFastDigitizer )
+DECLARE_ALGORITHM_FACTORY( MM_FastDigitizer )
+
+DECLARE_FACTORY_ENTRIES(MuonFastDigitization) {
+ DECLARE_ALGORITHM( sTgcFastDigitizer )
+ DECLARE_ALGORITHM( MM_FastDigitizer )
+}
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_load.cxx b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_load.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c436585549612e5641d7c83944d562f4524e5c19
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/components/MuonFastDigitization_load.cxx
@@ -0,0 +1,3 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES(MuonFastDigitization)
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.cxx b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a6d330d5b6995ddd6386834a37e8c177339b8f55
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.cxx
@@ -0,0 +1,674 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "sTgcFastDigitizer.h"
+#include "MuonSimEvent/sTgcSimIdToOfflineId.h"
+
+#include "MuonSimEvent/GenericMuonSimHitCollection.h"
+
+#include "MuonSimEvent/MicromegasHitIdHelper.h"
+#include "MuonReadoutGeometry/MuonDetectorManager.h"
+#include "MuonReadoutGeometry/sTgcReadoutElement.h"
+#include "MuonIdHelpers/sTgcIdHelper.h"
+
+#include "MuonSimData/MuonSimDataCollection.h"
+#include "MuonSimData/MuonSimData.h"
+#include "MuonPrepRawData/sTgcPrepDataContainer.h"
+#include "MuonPrepRawData/sTgcPrepData.h"
+#include "TrkEventPrimitives/LocalDirection.h"
+#include "TrkSurfaces/Surface.h"
+
+#include "Identifier/Identifier.h"
+#include "TH1.h"
+#include "TTree.h"
+#include "TFile.h"
+
+//Random Numbers
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandGauss.h"
+#include "CLHEP/Random/RandFlat.h"
+#include "PathResolver/PathResolver.h"
+
+using namespace Muon;
+
+sTgcFastDigitizer::sTgcFastDigitizer(const std::string& name, ISvcLocator* pSvcLocator)
+ : AthAlgorithm(name, pSvcLocator)
+ , m_idHelperTool("Muon::MuonIdHelperTool/MuonIdHelperTool" )
+ , m_muonClusterCreator("Muon::MuonClusterOnTrackCreator/MuonClusterOnTrackCreator")
+ , m_rndmSvc("AtRndmGenSvc", name )
+ , m_timeWindowOffsetWire(0.)
+ , m_timeWindowOffsetStrip(0.)
+ , m_timeWindowWire(24.95) // TGC 29.32; // 29.32 ns = 26 ns + 4 * 0.83 ns
+ , m_timeWindowStrip(24.95) // TGC 40.94; // 40.94 ns = 26 ns + 18 * 0.83 ns
+ , m_bunchCrossingTime(24.95) // 24.95 ns =(40.08 MHz)^(-1)
+{
+ declareProperty("IdHelperTool", m_idHelperTool);
+ declareProperty("ClusterCreator", m_muonClusterCreator);
+ declareProperty("ChannelTypes", m_channelTypes = 3);
+ declareProperty("RndmEngine", m_rndmEngineName, "Random engine name");
+ declareProperty("RndmSvc", m_rndmSvc, "Random Number Service used in Muon digitization");
+ declareProperty("EnergyThreshold", m_energyThreshold = 50, "Minimal energy of incoming particle to produce a PRD" );
+ declareProperty("EnergyDepositThreshold", m_energyDepositThreshold = 0.00052, "Minimal energy deposit to produce a PRD" );
+ declareProperty("CheckIds", m_checkIds = false, "Turn on validity checking of Identifiers" );
+}
+
+sTgcFastDigitizer::~sTgcFastDigitizer() {
+
+}
+
+StatusCode sTgcFastDigitizer::initialize() {
+
+ if ( detStore()->retrieve( m_detManager ).isFailure()) {
+ ATH_MSG_WARNING("Failed to retrieve MuonDetectorManager");
+ return StatusCode::FAILURE;
+ }
+ if( detStore()->retrieve( m_idHelper ).isFailure() ){
+ ATH_MSG_WARNING("Failed to retrieve sTgcIdHelper");
+ return StatusCode::FAILURE;
+ }
+ m_rndmEngine = m_rndmSvc->GetEngine(m_rndmEngineName);
+ if (m_rndmEngine==0) {
+ ATH_MSG_ERROR("Could not find RndmEngine : " << m_rndmEngineName );
+ return StatusCode::FAILURE;
+ }
+ if( m_idHelperTool.retrieve().isFailure() ){
+ ATH_MSG_WARNING("Failed to retrieve " << m_idHelperTool);
+ return StatusCode::FAILURE;
+ }
+
+ if( m_muonClusterCreator.retrieve().isFailure() ){
+ ATH_MSG_WARNING("Failed to retrieve " << m_muonClusterCreator);
+ return StatusCode::FAILURE;
+ }
+
+ if( !readFileOfTimeJitter() ) return StatusCode::FAILURE;
+
+ m_file = new TFile("sTgcplots.root","RECREATE");
+
+ m_ntuple = new TTree("a","a");
+ m_ntuple->Branch("slx",&slx);
+ m_ntuple->Branch("sly",&sly);
+ m_ntuple->Branch("slz",&slz);
+ m_ntuple->Branch("dlx",&dlx);
+ m_ntuple->Branch("dly",&dly);
+ m_ntuple->Branch("dlz",&dlz);
+ m_ntuple->Branch("sulx",&sulx);
+ m_ntuple->Branch("suly",&suly);
+ m_ntuple->Branch("tsulx",&tsulx);
+ m_ntuple->Branch("tsuly",&tsuly);
+ m_ntuple->Branch("tsulz",&tsulz);
+ m_ntuple->Branch("resx",&resx);
+ m_ntuple->Branch("resy",&resy);
+ m_ntuple->Branch("suresx",&suresx);
+ m_ntuple->Branch("suresy",&suresy);
+ m_ntuple->Branch("resz",&resz);
+ m_ntuple->Branch("errx",&errx);
+ m_ntuple->Branch("erry",&erry);
+ m_ntuple->Branch("res",&res);
+ m_ntuple->Branch("pull",&pull);
+ m_ntuple->Branch("is",&is);
+ m_ntuple->Branch("seta",&seta);
+ m_ntuple->Branch("sphi",&sphi);
+ m_ntuple->Branch("sml",&sml);
+ m_ntuple->Branch("sl",&sl);
+ m_ntuple->Branch("ss",&ss);
+ m_ntuple->Branch("stype",&stype);
+ m_ntuple->Branch("ieta",&ieta);
+ m_ntuple->Branch("iphi",&iphi);
+ m_ntuple->Branch("iml",&iml);
+ m_ntuple->Branch("il",&il);
+ m_ntuple->Branch("ich",&ich);
+ m_ntuple->Branch("istr",&istr);
+ m_ntuple->Branch("itype",&itype);
+ m_ntuple->Branch("ipadeta",&ipadeta);
+ m_ntuple->Branch("ipadphi",&ipadphi);
+ m_ntuple->Branch("exitcode",&exitcode);
+ m_ntuple->Branch("mode",&mode);
+ m_ntuple->Branch("pdg",&pdg);
+ m_ntuple->Branch("trkid",&trkid);
+ m_ntuple->Branch("bct",&bct);
+ m_ntuple->Branch("tb",&tb);
+ m_ntuple->Branch("tj",&tj);
+ m_ntuple->Branch("tg4",&tg4);
+ m_ntuple->Branch("ttof",&ttof);
+ m_ntuple->Branch("gpx",&gpx);
+ m_ntuple->Branch("gpy",&gpy);
+ m_ntuple->Branch("gpz",&gpz);
+ m_ntuple->Branch("gpr",&gpr);
+ m_ntuple->Branch("gpp",&gpp);
+ m_ntuple->Branch("dgpx",&dgpx);
+ m_ntuple->Branch("dgpy",&dgpy);
+ m_ntuple->Branch("dgpz",&dgpz);
+ m_ntuple->Branch("dgpr",&dgpr);
+ m_ntuple->Branch("dgpp",&dgpp);
+ m_ntuple->Branch("edep",&edep);
+ m_ntuple->Branch("e",&e);
+ m_ntuple->Branch("at",&at);
+ m_ntuple->Branch("as",&as);
+
+ return StatusCode::SUCCESS;
+
+}
+
+StatusCode sTgcFastDigitizer::execute() {
+
+// Create and record the SDO container in StoreGate
+ MuonSimDataCollection* sdoContainer = new MuonSimDataCollection();
+ if (evtStore()->record(sdoContainer,"STGC_SDO").isFailure()) {
+ ATH_MSG_ERROR ( "Unable to record sTgc SDO collection in StoreGate" );
+ return StatusCode::FAILURE;
+ }
+
+ sTgcPrepDataContainer* prdContainer = new sTgcPrepDataContainer(m_idHelper->detectorElement_hash_max());
+
+ // as the sTgcPrepDataContainer only allows const accesss, need a local vector as well.
+ std::vector<sTgcPrepDataCollection*> localsTgcVec(m_idHelper->detectorElement_hash_max());
+
+ const DataHandle< GenericMuonSimHitCollection > collGMSH;
+ if ( evtStore()->retrieve( collGMSH,"sTGCSensitiveDetector").isFailure()) {
+ ATH_MSG_WARNING("No sTgc hits found in SG");
+ return StatusCode::FAILURE;
+ }
+
+ ATH_MSG_DEBUG( "Retrieved " << collGMSH->size() << " sTgc hits!");
+ sTgcHitIdHelper* hitHelper=sTgcHitIdHelper::GetHelper();
+ sTgcSimIdToOfflineId simToOffline(*m_idHelper);
+
+ const GenericMuonSimHit* previousHit = 0;
+
+ std::map<Identifier,int> hitsPerChannel;
+ int nhits = 0;
+ GenericMuonSimHitCollection::const_iterator itersTgc;
+ for (itersTgc=collGMSH->begin();itersTgc!=collGMSH->end();++itersTgc) {
+ const GenericMuonSimHit& hit = *itersTgc;
+
+ float globalHitTime = hit.globalpreTime();
+ float tofCorrection = hit.globalPosition().mag()/CLHEP::c_light;
+
+ // bunch time
+ float bunchTime = globalHitTime - tofCorrection;
+
+ int simId = hit.GenericId();
+ Identifier layid = simToOffline.convert(simId);
+ ATH_MSG_VERBOSE("sTgc hit: r " << hit.globalPosition().perp() << " z " << hit.globalPosition().z() << " mclink " << hit.particleLink()
+ << " -- " << m_idHelperTool->toString(layid));
+
+ std::string stName = m_idHelper->stationNameString(m_idHelper->stationName(layid));
+ int isSmall = stName[2] == 'S';
+ const MuonGM::sTgcReadoutElement* detEl = m_detManager->getsTgcReadoutElement(layid);
+ if( !detEl ){
+ ATH_MSG_WARNING("Failed to retrieve detector element for: " << m_idHelperTool->toString(layid) );
+ continue;
+ }
+
+ IdentifierHash hash;
+ m_idHelper->get_detectorElement_hash(layid, hash);
+ MuonPrepDataCollection<Muon::sTgcPrepData>* col = localsTgcVec[hash];
+ if( !col ){
+ col = new sTgcPrepDataCollection(hash);
+ col->setIdentifier(m_idHelper->channelID(m_idHelper->parentID(layid), m_idHelper->multilayer(layid),1,1,1) );
+ if( prdContainer->addCollection(col,hash).isFailure() ){
+ ATH_MSG_WARNING("Failed to add collection with hash " << (int)hash );
+ delete col;col=0;
+ continue;
+ }
+ localsTgcVec[hash] = col;
+ }
+
+ Identifier parentId = m_idHelper->parentID(layid);
+
+ // SimHits without energy loss are not recorded.
+ // not needed because of already done in sensitive detector
+ // https://svnweb.cern.ch/trac/atlasoff/browser/MuonSpectrometer/MuonG4/MuonG4SD/trunk/src/sTGCSensitiveDetector.cxx?rev=542333#L66
+ // if(hit.depositEnergy()==0.) continue;
+
+
+ if( previousHit && abs(hit.particleEncoding())==13 && abs(previousHit->particleEncoding())==13 ) {
+ Amg::Vector3D diff = previousHit->localPosition() - hit.localPrePosition();
+ if( diff.mag() < 0.1 ) {
+ ATH_MSG_VERBOSE("second hit from a muon: prev " << previousHit->localPosition() << " current " << hit.localPrePosition()
+ << " diff " << diff );
+ continue;
+ }
+ }
+
+ // select whether to produce only strips or strips + wires or strips + wires + pads
+ int ftype = m_channelTypes == 3 ? 0 : 1;
+ int ltype = m_channelTypes == 1 ? 1 : 2;
+
+ for( int type=ftype;type<=ltype;++type ){
+
+ // only produce wire hits for outer two stations
+ if( type == 2 && abs(m_idHelper->stationEta(layid)) < 3 ) continue;
+
+ // for debugging only treat pads
+ //if( type != 0 ) continue;
+
+ // first create surface identifier
+ Identifier id = m_idHelper->channelID(parentId, m_idHelper->multilayer(layid), m_idHelper->gasGap(layid),type,1,m_checkIds);
+ ATH_MSG_VERBOSE("handling layer " << m_idHelperTool->toString(id) << " type " << type );
+ const Trk::PlaneSurface& surf = detEl->surface(id);
+ Amg::Transform3D gToL = detEl->absTransform().inverse();
+ Amg::Vector3D hpos(hit.globalPosition().x(),hit.globalPosition().y(),hit.globalPosition().z());
+ Amg::Vector3D lpos = gToL*hpos;
+ Amg::Vector3D slpos = surf.transform().inverse()*hpos;
+
+ // propagate sim hit position to surface
+ Amg::Vector3D ldir = surf.transform().inverse().linear()*Amg::Vector3D(hit.globalDirection().x(), hit.globalDirection().y(), hit.globalDirection().z()); //sleontsi
+ double scale = -slpos.z()/ldir.z();
+ Amg::Vector3D hitOnSurface = slpos + scale*ldir;
+
+ // the position resolution is depending on the incident angle
+ //double resolution = .18;
+ const Amg::Vector3D GloDire(hit.globalDirection().x(), hit.globalDirection().y(), hit.globalDirection().z());
+ Trk::LocalDirection LocDirection;
+ surf.globalToLocalDirection(GloDire, LocDirection);
+
+ float inAngle_time = 90 - fabs( LocDirection.angleXZ() / CLHEP::degree);
+ float inAngle_space = 90 - fabs( LocDirection.angleYZ() / CLHEP::degree);
+ if(inAngle_time > 90) inAngle_time = inAngle_time -90.;
+ if(inAngle_space > 90) inAngle_space = inAngle_space -90.;
+
+
+
+ // bctagging
+// static const float jitterInitial = 9999.;
+ float jitter = 0;//jitterInitial; // calculated at central strip but also used in all the strips fired by the same hit
+// if(jitter > jitterInitial-0.1) {
+// jitter = timeJitter(inAngle_time);
+// }
+
+ //const float stripPropagationTime = 3.3*CLHEP::ns/CLHEP::m * detEl->distanceToReadout(posOnSurf_strip, elemId); // 8.5*ns/m was used until MC10.
+ const float stripPropagationTime = 0.; // 8.5*ns/m was used until MC10.
+
+ // shift by 0.5 to get all current bc in time
+ float sDigitTime = bunchTime + jitter + stripPropagationTime+0.5;
+
+ uint16_t bctag = bcTagging(sDigitTime, type);
+
+ // smear
+ // for the strips apply gaussian smearing, for the wire gangs or pads do nothing here
+ double sp = hitOnSurface.x();
+ double resolution = 0;
+ if( type == 1 ){
+ resolution = getResolution(inAngle_time);
+ sp = CLHEP::RandGauss::shoot(m_rndmEngine, hitOnSurface.x(), resolution);
+ }
+ //ATH_MSG_DEBUG("slpos.z " << slpos.z() << ", ldir " << ldir.z() << ", scale " << scale << ", hitOnSurface.z " << hitOnSurface.z());
+
+
+ Amg::Vector2D posOnSurf(sp,hitOnSurface.y());
+ Amg::Vector2D locHitPosOnSurf(hitOnSurface.x(),hitOnSurface.y());
+
+ /// now fill most of the ntuple
+ Amg::Vector3D repos = detEl->globalPosition();
+
+ slx = hit.localPosition().x();
+ sly = hit.localPosition().y();
+ slz = hit.localPosition().z();
+ dlx = lpos.x();
+ dly = lpos.y();
+ dlz = lpos.z();
+ sulx = posOnSurf.x();
+ suly = posOnSurf.y();
+ tsulx = hitOnSurface.x();
+ tsuly = hitOnSurface.y();
+ tsulz = hitOnSurface.z();
+ resx = hit.localPosition().x() - lpos.x();
+ resy = hit.localPosition().y() - lpos.y();
+ resz = hit.localPosition().z() - lpos.z();
+ suresx = posOnSurf.x()-hitOnSurface.x();
+ suresy = posOnSurf.y()-hitOnSurface.y();
+ errx = -99999.;
+ erry = -99999.;
+ res = -99999.;
+ pull = -99999.;
+ is = isSmall;
+ seta = hitHelper->GetZSector(simId);
+ sphi = hitHelper->GetPhiSector(simId);
+ sml = hitHelper->GetMultiLayer(simId);
+ sl = hitHelper->GetLayer(simId);
+ ss = hitHelper->GetSide(simId);
+ stype = type;
+ ieta = m_idHelper->stationEta(id);
+ iphi = m_idHelper->stationPhi(id);
+ iml = m_idHelper->multilayer(id);
+ il = m_idHelper->gasGap(id);
+ ich = -99999;
+ istr = -99999;
+ itype = -99999;
+ ipadeta = -99999;
+ ipadphi = -99999;
+ exitcode = 0;
+ mode = 0;
+ pdg = hit.particleEncoding();
+ trkid = hit.trackNumber();
+ bct = bctag;
+ tb = bunchTime;
+ tj = jitter;
+ tg4 = globalHitTime;
+ ttof = tofCorrection;
+ gpx = hit.globalPosition().x();
+ gpy = hit.globalPosition().y();
+ gpz = hit.globalPosition().z();
+ gpr = hit.globalPosition().perp();
+ gpp = hit.globalPosition().phi();
+ dgpx = repos.x();
+ dgpy = repos.y();
+ dgpz = repos.z();
+ dgpr = repos.perp();
+ dgpp = repos.phi();
+ edep = hit.depositEnergy();
+ e = hit.kineticEnergy();
+ at = inAngle_time;
+ as = inAngle_space;
+ // cut on the kineticEnergy = 50MeV
+ if(hit.kineticEnergy()< m_energyThreshold ) {
+ exitcode = 5;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // cut on depositEnergy(0.52KeV) to simulation the detector efficiency(95% for strips)
+ if( type ==1 && hit.depositEnergy()<m_energyDepositThreshold) {
+ ATH_MSG_VERBOSE("Drop SimHit with depositEnergy = " << edep << " in the strip response!");
+ exitcode = 6;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ /// done filling
+ if ( !surf.insideBounds(posOnSurf) ) {
+ exitcode = 1;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ if(sDigitTime < -m_bunchCrossingTime+m_timeWindowOffsetStrip || +m_bunchCrossingTime+m_timeWindowOffsetStrip+ m_timeWindowStrip < sDigitTime) {
+ exitcode = 4;
+ m_ntuple->Fill();
+ ATH_MSG_DEBUG( "Strip: Digitized time is outside of time window. " << sDigitTime
+ << " bunchTime: " << bunchTime << " time jitter: " << jitter
+ << " propagation time: " << stripPropagationTime );
+ continue;
+ }
+
+ bctag = bcTagging(sDigitTime, type);
+
+ int stripNumber = detEl->stripNumber(locHitPosOnSurf,id);
+ if( stripNumber == -1 ){
+ ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperTool->toString(id)
+ << " pos " << posOnSurf.x() << " - " << hitOnSurface.x() << " z " << slpos.z() << " type " << type );
+ exitcode = 2;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // create channel identifier
+ id = m_idHelper->channelID(parentId, m_idHelper->multilayer(id), m_idHelper->gasGap(id),type,stripNumber,m_checkIds);
+ ATH_MSG_VERBOSE(" Unsmeared hit id " << m_idHelperTool->toString(id) );
+
+ int& counts = hitsPerChannel[id];
+ ++counts;
+ if( counts > 1 ) continue;
+ ++nhits;
+
+
+ // recalculate using smeared position
+ stripNumber = detEl->stripNumber(posOnSurf,id);
+ if( stripNumber == -1 ){
+ ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperTool->toString(id)
+ << " pos " << posOnSurf.x() << " - " << hitOnSurface.x() << " z " << slpos.z() << " type " << type );
+ exitcode = 2;
+ m_ntuple->Fill();
+ continue;
+ }
+
+ // create channel identifier
+ id = m_idHelper->channelID(parentId, m_idHelper->multilayer(id), m_idHelper->gasGap(id),type,stripNumber,m_checkIds);
+
+ if( type != m_idHelper->channelType(id) ) ATH_MSG_WARNING(" bad type: in " << type << " from id " << m_idHelper->channelType(id)
+ << " strip " << stripNumber << " from id " << m_idHelper->channel(id)
+ << " eta " << m_idHelper->stationEta(id)
+ << " local pos " << slpos );
+
+ // assign strip position to PRD for wires and pads
+ if( type != 1 ){
+ Amg::Vector2D locpos(0,0);
+ if( !detEl->stripPosition(id,locpos ) ){
+ ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperTool->toString(id) );
+ exitcode = 3;
+ m_ntuple->Fill();
+ continue;
+ }
+ posOnSurf = locpos;
+ }
+
+ std::vector<Identifier> rdoList;
+ rdoList.push_back(id);
+
+ double errX = 500;
+
+ if( type == 2 ){
+ const MuonGM::MuonChannelDesign* design = detEl->getDesign(id);
+ if( !design ){
+ ATH_MSG_WARNING("Failed to get design for " << m_idHelperTool->toString(id) );
+ }else{
+ errX = fabs(design->inputPitch)/sqrt(12);
+ }
+ }else if( type == 1 ){
+ errX = resolution;
+ }
+ if( type == 0 ){
+
+ const MuonGM::MuonPadDesign* design = detEl->getPadDesign(id);
+ if( !design ){
+ ATH_MSG_WARNING("Failed to get design for " << m_idHelperTool->toString(id) );
+ }else{
+ errX = design->channelWidth(posOnSurf,true)/sqrt(12);
+ }
+ }
+
+
+ Amg::MatrixX cov(1,1); cov.setIdentity();
+ cov*=errX*errX;
+
+// ATH_MSG_DEBUG(" New hit " << m_idHelperTool->toString(id) << " chtype " << type << " lpos " << posOnSurf << " from truth "
+// << hitOnSurface << " error " << locErrMat->error(Trk::locX) << " pull " << (posOnSurf.x()-hitOnSurface.x())/locErrMat->error(Trk::locX) );
+
+ //sTgcPrepData* prd = new sTgcPrepData( id,hash,posOnSurf,rdoList,locErrMat,detEl);
+ sTgcPrepData* prd = new sTgcPrepData( id,hash,posOnSurf,rdoList,&cov,detEl, bctag);
+ prd->setHashAndIndex(col->identifyHash(), col->size());
+ col->push_back(prd);
+
+ /// fill final bits of the ntuple
+ suresx = posOnSurf.x()-hitOnSurface.x();
+ suresy = posOnSurf.y()-hitOnSurface.y();
+ sulx = posOnSurf.x();
+ suly = posOnSurf.y();
+ errx = errX;
+ ich = m_idHelper->channel(id);
+ istr = stripNumber;
+ itype = m_idHelper->channelType(id);
+ if( type == 0 ){
+ ipadeta = m_idHelper->padEta(id);
+ ipadphi = m_idHelper->padPhi(id);
+ }
+
+// const MuonClusterOnTrack* rot = m_muonClusterCreator->createRIO_OnTrack( *prd, hit.globalPosition() );
+// if( rot ){
+// res = rot->localParameters().get(Trk::locX)-hitOnSurface.x();
+// pull = res/rot->localErrorMatrix().error(Trk::locX);
+// delete rot;
+// }
+
+ m_ntuple->Fill();
+ // create SDO
+ MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(hitOnSurface.x(),hitOnSurface.y()));
+ //Record the SDO collection in StoreGate
+ std::vector<MuonSimData::Deposit> deposits;
+ deposits.push_back(deposit);
+ sdoContainer->insert ( std::make_pair ( id, MuonSimData(deposits,0) ) );
+
+ previousHit = &hit;
+
+ } // end of loop channelType
+ } // end of loop SimHits
+ if( msgLvl(MSG::DEBUG) ){
+ std::map<Identifier,int>::const_iterator hit = hitsPerChannel.begin();
+ std::map<Identifier,int>::const_iterator hit_end = hitsPerChannel.end();
+ msg(MSG::DEBUG) << " number of channels with hit " << hitsPerChannel.size() << " nhits " << nhits << std::endl;
+ for( ;hit!=hit_end;++hit ){
+ msg(MSG::DEBUG) << " " << m_idHelperTool->toString(hit->first) << " -> " << hit->second << std::endl;
+ }
+ msg(MSG::DEBUG) << endreq;
+ }
+ std::string key = "STGC_Measurements";
+ ATH_MSG_DEBUG(" Done! Total number of sTgc chambers with PRDS: " << prdContainer->numberOfCollections() << " key " << key);
+ if (evtStore()->record(prdContainer,key).isFailure()) {
+ ATH_MSG_ERROR ( "Unable to record sTgcPrepData container in StoreGate" );
+ return StatusCode::FAILURE;
+ }
+ return StatusCode::SUCCESS;
+}
+
+StatusCode sTgcFastDigitizer::finalize() {
+ //m_ntuple_SimHit->Write();
+ m_ntuple->Write();
+ m_file->Write();
+ m_file->Close();
+
+ return StatusCode::SUCCESS;
+}
+
+
+//+++++++++++++++++++++++++++++++++++++++++++++++
+double sTgcFastDigitizer::getResolution(float inAngle_space) const
+{
+ // 0 degree --> 0.18mm; 10 degree -->
+ const int NIndex = 2;
+ //double Reso[NIndex] = {0.18, 0.20, 0.24};
+ double Reso[NIndex] = {0.06, 0.30};
+ double Angle[NIndex] = {0., 60.};
+
+ int index = 1;
+ inAngle_space = fabs(inAngle_space);
+ while(index<NIndex && inAngle_space>Angle[index]){
+ index++;
+ }
+
+ double Resolution = 0.;
+
+ if(index==NIndex) Resolution = Reso[NIndex-1];
+ else Resolution = Reso[index-1] + (inAngle_space - Angle[index-1]) / (Angle[index]-Angle[index-1]) * (Reso[index] - Reso[index-1]);
+
+ ATH_MSG_VERBOSE("sTgcFastDigitizer::getResolution(" << inAngle_space << ") = " <<Resolution );
+
+ return Resolution;
+}
+
+
+//+++++++++++++++++++++++++++++++++++++++++++++++
+uint16_t sTgcFastDigitizer::bcTagging(float digitTime, int channelType) const {
+
+ double offset, window;
+ if(channelType == 1) { // strip
+ offset = m_timeWindowOffsetStrip;
+ window = m_timeWindowStrip;
+ }
+ else { // wire gangs or pad
+ offset = m_timeWindowOffsetWire;
+ window = m_timeWindowWire;
+ }
+
+ uint16_t bctag = 0;
+ if(-m_bunchCrossingTime+offset < digitTime && digitTime < -m_bunchCrossingTime+offset+window) {
+ bctag = (bctag | 0x1);
+ }
+ if( offset < digitTime && digitTime < offset+window) {
+ bctag = (bctag | 0x2);
+ }
+ if(+m_bunchCrossingTime+offset < digitTime && digitTime < +m_bunchCrossingTime+offset+window) {
+ bctag = (bctag | 0x4);
+ }
+ return bctag;
+}
+
+
+//+++++++++++++++++++++++++++++++++++++++++++++++
+float sTgcFastDigitizer::timeJitter(float inAngle_time) const
+{
+
+ int ithAngle = static_cast<int>(inAngle_time/5.);
+ float wAngle = inAngle_time/5. - static_cast<float>(ithAngle);
+ int jthAngle;
+ if (ithAngle > 11) {
+ ithAngle = 12;
+ jthAngle = 12;
+ }
+ else {
+ jthAngle = ithAngle+1;
+ }
+
+ float jitter;
+ float prob = 1.;
+ float probRef = 0.;
+
+ while (prob > probRef) {
+ prob = CLHEP::RandFlat::shoot(m_rndmEngine, 0.0, 1.0);
+ jitter = CLHEP::RandFlat::shoot(m_rndmEngine, 0.0, 1.0)*40.; // trial time jitter in nsec
+ int ithJitter = static_cast<int>(jitter);
+ // probability distribution calculated from weighted sum between neighboring bins of angles
+ probRef = (1.-wAngle)*m_vecAngle_Time[ithAngle][ithJitter]
+ + wAngle *m_vecAngle_Time[jthAngle][ithJitter];
+ }
+ return jitter;
+}
+
+
+//+++++++++++++++++++++++++++++++++++++++++++++++
+bool sTgcFastDigitizer::readFileOfTimeJitter()
+{
+
+ const char* const fileName = "sTGC_Digitization_timejitter.dat";
+ std::string fileWithPath = PathResolver::find_file (fileName, "DATAPATH");
+
+ std::ifstream ifs;
+ if (fileWithPath != "") {
+ ifs.open(fileWithPath.c_str(), std::ios::in);
+ }
+ else {
+ ATH_MSG_ERROR("readFileOfTimeJitter(): Could not find file " << fileName);
+ return false;
+ }
+
+ if(ifs.bad()){
+ ATH_MSG_ERROR("readFileOfTimeJitter(): Could not find file " << fileName);
+ return false;
+ }
+
+ int angle = 0;
+ int bins = 0;
+ int i = 0;
+ float prob = 0.;
+
+ while(ifs.good()){
+ ifs >> angle >> bins;
+ if (ifs.eof()) break;
+ if( msgLvl(MSG::VERBOSE) ) msg(MSG::VERBOSE) << "readFileOfTimeJitter(): Timejitter, angle, Number of bins, prob. dist.: " << angle << " " << bins << " ";
+ m_vecAngle_Time.resize(i + 1);
+ for (int j = 0; j < 41/*bins*/; j++) {
+ ifs >> prob;
+ m_vecAngle_Time[i].push_back(prob);
+ if( msgLvl(MSG::VERBOSE) ){
+ if( j == 0) msg(MSG::VERBOSE) << "readFileOfTimeJitter(): ";
+ msg(MSG::VERBOSE) << prob << " ";
+ }
+ }
+ if( msgLvl(MSG::VERBOSE) ) msg(MSG::VERBOSE) << endreq;
+ i++;
+ }
+ ifs.close();
+ return true;
+}
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.h b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..9fb4ad2cc2cd931032e528d3a36cc9ea3a7b6048
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/sTgcFastDigitizer.h
@@ -0,0 +1,146 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONDIGITIZATION_sTgcFASTDIGITIZER_H
+#define MUONDIGITIZATION_sTgcFASTDIGITIZER_H
+
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "StoreGate/StoreGateSvc.h"
+
+//Random
+#include "CLHEP/Random/RandomEngine.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandGauss.h"
+#include "MuonIdHelpers/MuonIdHelperTool.h"
+#include "MuonRecToolInterfaces/IMuonClusterOnTrackCreator.h"
+
+class TTree;
+class TFile;
+
+class sTgcIdHelper;
+namespace MuonGM {
+ class MuonDetectorManager;
+}
+
+class IAtRndmGenSvc;
+class ActiveStoreSvc;
+
+class sTgcFastDigitizer : public AthAlgorithm {
+
+ public:
+
+ sTgcFastDigitizer(const std::string& name, ISvcLocator* pSvcLocator);
+ ~sTgcFastDigitizer();
+
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+ ServiceHandle<IAtRndmGenSvc> getRndmSvc() const { return m_rndmSvc; } // Random number service
+ CLHEP::HepRandomEngine *getRndmEngine() const { return m_rndmEngine; } // Random number engine used
+
+ private:
+ const MuonGM::MuonDetectorManager* m_detManager;
+ const sTgcIdHelper* m_idHelper;
+
+ int m_channelTypes; // 1 -> strips, 2 -> strips+wires, 3 -> strips/wires/pads
+
+ TFile* m_file;
+ //TTree* m_ntuple_SimHit;
+ TTree* m_ntuple;
+ float slx; // local position simhit in G4
+ float sly;
+ float slz;
+ float dlx; // local position simhit in GeoModel frame
+ float dly;
+ float dlz;
+ float sulx; // local position from REin tracking surface frame
+ float suly;
+ float tsulx; // local position simhit in tracking surface frame
+ float tsuly;
+ float tsulz;
+ float resx; // residuals in local G4 - GeoModel
+ float resy;
+ float resz;
+ float suresx; // residuals g4 - RE stripPosition in tracking frame
+ float suresy;
+ float errx; // error on x and y
+ float erry;
+ float res; // residual + pull
+ float pull;
+ int is; // simulation Identifier: issmall
+ int seta; // eta
+ int sphi; // phi
+ int sml; // multi layer
+ int sl; // layer
+ int ss; // side
+ int stype; // channel type
+ int ieta; // offline id fields
+ int iphi;
+ int iml;
+ int il;
+ int ich; // strip number of Id
+ int istr; // strip number from RE
+ int itype; // channel type
+ int ipadeta; // channel type
+ int ipadphi; // channel type
+ int exitcode; // flag reason why prd not made
+ int mode; // flag digitization mode
+ int pdg; // pdg ID
+ int trkid; // track number in G4
+ int bct; // bc tag
+ float tb; // bunchTime
+ float tj; // time jitter
+ float tg4; // globalHitTime in G4
+ float ttof; // tofCorrection assuming IP;
+ float gpx; // global position of the simhit
+ float gpy;
+ float gpz;
+ float gpr; // radial position pos.perp()
+ float gpp; // phi position pos.phi()
+ float dgpx; // same for readout element center
+ float dgpy;
+ float dgpz;
+ float dgpr;
+ float dgpp;
+ float edep;
+ float e;
+ float as;
+ float at;
+
+
+ double getResolution(float inAngle_space) const;
+ uint16_t bcTagging(float digittime, int channelType) const;
+ float timeJitter(float inAngle_time) const;
+
+ /**
+ Reads parameters for intrinsic time response from timejitter.dat.
+ */
+ bool readFileOfTimeJitter();
+
+
+ protected:
+ ToolHandle <Muon::MuonIdHelperTool> m_idHelperTool; // IdHelperTool
+ ToolHandle <Muon::IMuonClusterOnTrackCreator> m_muonClusterCreator;
+ ServiceHandle <IAtRndmGenSvc> m_rndmSvc; // Random number service
+ CLHEP::HepRandomEngine *m_rndmEngine; // Random number engine used - not init in SiDigitization
+ std::string m_rndmEngineName;// name of random engine
+ std::string m_inputObjectName; // name of the input objects
+ double m_timeWindowOffsetWire;
+ double m_timeWindowOffsetStrip;
+ double m_timeWindowWire;
+ double m_timeWindowStrip;
+ double m_bunchCrossingTime;
+ double m_energyThreshold;
+ double m_energyDepositThreshold;
+ bool m_checkIds;
+ std::vector<std::vector<float> > m_vecAngle_Time;
+
+};
+
+#endif // MUONDIGITIZATION_sTgcDIGITIZER_H
+
+