diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/ElectronicsResponse.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/ElectronicsResponse.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd4cc0f53f446f19710e93ce7db3ca9593631c31
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/ElectronicsResponse.h
@@ -0,0 +1,108 @@
+/* -*- C++ -*- */
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef ElectronicsResponse_h
+#define ElectronicsResponse_h
+/** @class ElectronicsResponse
+
+// ------------
+// Authors:
+// Iakovidis George <george.iakovidis@cern.ch>
+// Karakostas Konstantinos <konstantinos.karakostas@cern.ch>
+// Leontsinis Stefanos <stefanos.leontsinis@cern.ch>
+// Nektarios Chr. Benekos <nbenekos@cern.ch>
+// Jessica Metcalfe <jessica.metcalfe@gmail.com>
+//////////////////////////////////////////////////////////////////////////////
+
+Comments to be added here....
+
+
+*/
+
+/// ROOT
+#include <TROOT.h>
+#include <TH1.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TMath.h>
+
+/// STD'S
+#include <algorithm>
+#include <cmath>
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+#include <utility>
+#include <string>
+#include <sstream>
+#include <sys/stat.h>
+
+/// Projects
+#include "MM_Digitization/MmElectronicsToolInput.h"
+#include "MM_Digitization/MmDigitToolOutput.h"
+#include "MM_Digitization/StripsResponse.h"
+
+using std::vector;
+using std::cout;
+using std::endl;
+
+/// ROOT Classed
+class TF1;
+class TH1;
+class TH1F;
+class TH2F;
+
+class ElectronicsResponse {
+
+private:
+ /** power of responce function */
+ float alpha;
+ /** */
+ float RC ;
+ /** hreshold "Voltage" for histoBNL */
+ float electronicsThreshold;
+ float m_StripResponse_qThreshold;
+ float m_StripResponse_driftGap;
+ float m_StripResponse_driftVelocity;
+
+ TF1 *intFn;
+ StripsResponse* stripObject ;
+public :
+
+ ElectronicsResponse();
+ virtual ~ElectronicsResponse();
+ void clearValues ();
+ void bnlResponceFunction(const vector <int> & numberofStrip, const vector <float> & qStrip, const vector <float> & tStrip);
+ MmDigitToolOutput GetResponceFrom(const MmElectronicsToolInput & digiInput);
+
+ float tMinFromIntegration;
+ float tminFromIntegrationAboveThreshold;
+ float tMinFromFirstPeak;
+ vector <float> tStripElectronicsAbThr;
+ vector <float> qStripElectronics;
+ vector <int> nStripElectronics;
+
+ inline void set_alpha(float val) {alpha = val;};
+ inline void set_RC (float val) {RC = val;};
+ inline void set_electronicsThreshold(float val) {electronicsThreshold = val;};
+
+ float get_alpha() const { return alpha;};
+ float get_RC () const { return RC;};
+ float get_electronicsThreshold() const { return electronicsThreshold;};
+
+ float get_tMinFromIntegration () const { return tMinFromIntegration;};
+ float get_tminFromIntegrationAboveThreshold () const { return tminFromIntegrationAboveThreshold;};
+ float get_tMinFromFirstPeak () const { return tMinFromFirstPeak;};
+
+ vector <float> get_tStripElectronicsAbThr () const { return tStripElectronicsAbThr;};
+ vector <float> get_qStripElectronics () const { return qStripElectronics;};
+ vector <int> get_nStripElectronics () const { return nStripElectronics;};
+
+
+};
+
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/IMM_DigitizationTool.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/IMM_DigitizationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d24f75a54e85d38f6e0b03ebecc93890779a434
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/IMM_DigitizationTool.h
@@ -0,0 +1,30 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_IMM_DIGITIZATIONTOOL_H
+#define MM_DIGITIZATION_IMM_DIGITIZATIONTOOL_H
+
+#include "MM_Digitization/MmDigitToolOutput.h"
+#include "GaudiKernel/IAlgTool.h"
+/*-----------------------------------------------
+
+Created March 2013 by Nektarios Chr. Benekos
+
+Interface for tools which convert MicroMegas digitization input quantities into the signal
+-----------------------------------------------*/
+
+class MmDigitToolInput;
+
+static const InterfaceID IID_IMM_DigitizationTool("IMM_DigitizationTool",1,0);
+
+class IMM_DigitizationTool : virtual public IAlgTool {
+ public:
+
+ virtual MmDigitToolOutput digitize(/*const MmDigitToolInput& input*/) = 0;
+
+ static const InterfaceID& interfaceID()
+ { return IID_IMM_DigitizationTool; }
+};
+
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Digitizer.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Digitizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..5b9677399c7e56a0e8616cf54e65b0239e09a087
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Digitizer.h
@@ -0,0 +1,34 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONDIGITIZATION_MM_DIGITIZER_H
+#define MUONDIGITIZATION_MM_DIGITIZER_H
+
+/// Gaudi External
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+
+class IMuonDigitizationTool;
+
+/*******************************************************************************/
+class MM_Digitizer : public AthAlgorithm {
+
+ public:
+
+ MM_Digitizer(const std::string& name, ISvcLocator* pSvcLocator);
+ ~MM_Digitizer();
+
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+ private:
+
+ ToolHandle<IMuonDigitizationTool> m_digTool;
+
+};
+/*******************************************************************************/
+#endif // MUONDIGITIZATION_MM_DIGITIZER_H
+
+
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Response_DigitTool.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Response_DigitTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..59fe1263cbaa344f07fcbc59505d3be8b9d585d4
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_Response_DigitTool.h
@@ -0,0 +1,54 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_MM_RESPONSE_DIGITTOOL_H
+#define MM_DIGITIZATION_MM_RESPONSE_DIGITTOOL_H
+
+#include "GaudiKernel/AlgTool.h"
+#include "MM_Digitization/MmDigitToolOutput.h"
+#include "MM_Digitization/IMM_DigitizationTool.h"
+
+#include "GaudiKernel/ServiceHandle.h"
+
+#include "CLHEP/Random/RandomEngine.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+/*-----------------------------------------------
+
+Created March 2013 by Nektarios Chr. Benekos
+
+Digitization tool which uses MM_Response to convert MM digitization
+input quantities into the output
+
+-----------------------------------------------*/
+/*******************************************************************************/
+
+namespace MuonGM{
+ class MuonDetectorManager;
+}
+class MmIdHelper;
+class IAtRndmGenSvc;
+
+class MM_Response_DigitTool : public AlgTool, virtual public IMM_DigitizationTool {
+ public:
+ MM_Response_DigitTool( const std::string& type,
+ const std::string& name,
+ const IInterface* parent );
+
+ MmDigitToolOutput digitize(/* const MmDigitToolInput& input */ );
+ StatusCode initialize();
+ bool initializeStrip();
+
+ private:
+
+ const MuonGM::MuonDetectorManager* m_muonGeoMgr;
+ const MmIdHelper* m_idHelper;
+
+ protected:
+ CLHEP::HepRandomEngine *m_rndmEngine; // Random number engine used - not init in SiDigitization
+ std::string m_rndmEngineName;// name of random engine
+ ServiceHandle <IAtRndmGenSvc> m_rndmSvc; // Random number service
+
+};
+
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolInput.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolInput.h
new file mode 100644
index 0000000000000000000000000000000000000000..6fbb06e60d1a933602f89ca5e67f23c25998fdf0
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolInput.h
@@ -0,0 +1,54 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_MMDIGITTOOLINPUT_H
+#define MM_DIGITIZATION_MMDIGITTOOLINPUT_H
+#include "Identifier/Identifier.h"
+/*-----------------------------------------------
+
+Created March 2013 by Nektarios Chr. Benekos
+ Karakostas Konstantinos <Konstantinos.Karakostas@cern.ch>
+
+Class to store input needed for the MM_Digitization tools:
+- G4 driftradius
+- position along tube
+- magnetic field strength at hit position <-- ToDo: get the lorentz angle from the mag field
+- local temperature
+- electric charge
+- gamma factor
+- hit Identifier
+
+-----------------------------------------------*/
+/*******************************************************************************/
+class MmDigitToolInput {
+ public:
+
+ MmDigitToolInput(int stripIdLocal, double posx, double incomingAngle, double field, int stripMaxId)
+ : m_stripIDLocal(stripIdLocal),
+ m_xpos(posx),
+ m_incomingAngle(incomingAngle),
+ m_field(field),
+ m_stripMaxId(stripMaxId)
+ { }
+
+
+ ~MmDigitToolInput() {}
+
+ int stripIDLocal() const { return m_stripIDLocal; }
+ double positionWithinStrip() const { return m_xpos; }
+ double incomingAngle() const { return m_incomingAngle; }
+ double magneticField() const { return m_field; }
+ int stripMaxID() const { return m_stripMaxId; }
+ Identifier getHitID() const { return m_hitID; }
+
+ private:
+ int m_stripIDLocal;
+ double m_xpos;
+ double m_incomingAngle;
+ double m_field;
+ int m_stripMaxId;
+ Identifier m_hitID;
+};
+/*******************************************************************************/
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolOutput.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolOutput.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a71e86050cb51f312c6de771bb7494a56c52bf6
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitToolOutput.h
@@ -0,0 +1,59 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_MMDIGITTOOLOUTPUT_H
+#define MM_DIGITIZATION_MMDIGITTOOLOUTPUT_H
+
+/*-----------------------------------------------
+
+Created March 2013 by Nektarios Chr. Benekos
+ Karakostas Konstantinos <Konstantinos.Karakostas@cern.ch>
+
+Class to store output produced by MDT_Digitization tools:
+- was hit efficient
+- strip position
+- strip charge
+- strip time
+- first strip position [for Trigger study]
+- first strip time [for Trigger study]
+
+-----------------------------------------------*/
+
+class MmDigitToolOutput {
+ public:
+ MmDigitToolOutput(bool hitWasEff, std::vector <int> strpos, std::vector<float> time, std::vector<float> charge, int strTrig, float strTimeTrig )
+ : m_hitWasEff(hitWasEff),
+ m_strpos(strpos),
+ m_time(time),
+ m_charge(charge),
+ m_stripForTrigger(strTrig),
+ m_stripTimeForTrigger(strTimeTrig),
+ m_isValid(false)
+ {
+ if(m_strpos.size() > 0 && m_time.size() > 0 && m_charge.size() > 0) m_isValid = true;
+ }
+
+ ~MmDigitToolOutput() {}
+
+ bool hitWasEfficient() const { return m_hitWasEff; }
+ std::vector<int> stripPos() const { return m_strpos; }
+ std::vector<float> stripTime() const { return m_time; }
+ std::vector<float> stripCharge() const { return m_charge; }
+ int stripForTrigger() const { return m_stripForTrigger; }
+ float stripTimeForTrigger() const { return m_stripTimeForTrigger; }
+ void set_StripForTrigger(int val) { m_stripForTrigger = val; }
+ void set_StripTimeForTrigger(float val) { m_stripTimeForTrigger = val; }
+ bool isValid() { return m_isValid; }
+
+ private:
+ bool m_hitWasEff;
+ std::vector<int> m_strpos;
+ std::vector<float> m_time;
+ std::vector<float> m_charge;
+ int m_stripForTrigger;
+ float m_stripTimeForTrigger;
+ bool m_isValid;
+};
+
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitizationTool.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitizationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..b5515167e3c3a7fd8d1b82914e8e4355e1e52708
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmDigitizationTool.h
@@ -0,0 +1,275 @@
+/* -*- C++ -*- */
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATIONTOOL_H
+#define MM_DIGITIZATIONTOOL_H
+/** @class MmDigitizationTool
+
+ @section MM_DigitizerDetails Class methods and properties
+ // ------------
+ // Authors:
+ // Nektarios Chr. Benekos <nectarios.benekos@cern.ch>
+ // Konstantinos Karakostas <Konstantinos.Karakostas@cern.ch>
+ ////////////////////////////////////////////////////////////////////////////////
+
+ In the initialize() method, the PileUpMerge and StoreGate services are initialized, and a pointer to an instance of the class
+ MuonDetectorManager is retrieved from the detector store and used to obtain a MmIdHelper.
+ The MMDigitContainer is initialized and the simulation identifier helper retrieved, together with the pointer to the digitization tool.
+ Random numbers are obtained in the code from a dedicated stream via AtRndmSvc,
+ which is also initialized in the initialize() method.
+ In the execute() method, the digits and the SDOs (Simulation Data Object,
+ container for simulation data to be preserved after the digitization procedue, and persistified together with the RDOs)
+ containers are created and recorded on StoreGate; the GenericMuonSimHit collection are merged using the
+ TimedHitCollection sorted container (done in handleMicroMegasSimhit(TimedHitPtr<GenericMuonSimHit>& hit)) method);
+ into a loop over the TimedHitCollection for the given DetectorElement,
+ the handleMicroMegasSimhit() method converts the SimID into the Offline ID to be associated to the Digit and pass to the
+ digitization tool the drift radius and the distance to the chamber RO side (for the propagation delay computation).
+ The digitization tool returns a drift time, charge and strip position which are used together with the Offline ID,
+ to create the digit object (in doDigitization() method).
+ The finalize() method returns a SUCCESS StatusCode if the digitization procedure ends succesfully.
+
+ In the initialize() method...
+ In the execute() method...
+*/
+
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "HitManagement/TimedHitCollection.h"
+#include "MuonSimEvent/GenericMuonSimHitCollection.h"
+#include "MuonSimEvent/GenericMuonSimHit.h"
+#include "PileUpTools/PileUpToolBase.h"
+#include "Identifier/Identifier.h"
+
+#include "MM_Digitization/MmSortedHitVector.h"
+
+#include "CLHEP/Random/RandomEngine.h"
+#include "CLHEP/Geometry/Point3D.h"
+#include "CLHEP/Vector/ThreeVector.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include "MuonDigToolInterfaces/IMuonDigitizationTool.h"
+#include "MM_Digitization/StripsResponse.h"
+#include "MM_Digitization/ElectronicsResponse.h"
+
+#include <string>
+#include <sstream>
+#include <vector>
+#include <map>
+/*******************************************************************************/
+namespace MuonGM{
+ class MuonDetectorManager;
+ class MMReadoutElement;
+ class MuonChannelDesign;
+}
+namespace CLHEP{
+ class HepRandomEngine;
+}
+
+class StoreGateSvc;
+class ActiveStoreSvc;
+class PileUpMergeSvc;
+
+class MmDigitContainer;
+class MmDigitCollection;
+class MmIdHelper;
+class MicromegasHitIdHelper;
+class IAtRndmGenSvc;
+class MsgStream;
+
+class IMM_DigitizationTool;
+class MuonSimDataCollection;
+
+class StripsResponse;
+class ElectronicsResponse;
+
+class TTree;
+class TFile;
+
+/*******************************************************************************/
+
+class MmDigitizationTool : virtual public IMuonDigitizationTool, public PileUpToolBase {
+
+public:
+ MmDigitizationTool(const std::string& type, const std::string& name, const IInterface* parent);
+
+ /** Initialize */
+ virtual StatusCode initialize();
+
+ /** When being run from PileUpToolsAlgs, this method is called at the start of the subevts loop. Not able to access SubEvents */
+ StatusCode prepareEvent(const unsigned int /*nInputEvents*/);
+
+ /** When being run from PileUpToolsAlgs, this method is called for each active bunch-crossing to process current SubEvents bunchXing is in ns */
+ StatusCode processBunchXing(int bunchXing,
+ PileUpEventInfo::SubEvent::const_iterator bSubEvents,
+ PileUpEventInfo::SubEvent::const_iterator eSubEvents);
+
+ /** When being run from PileUpToolsAlgs, this method is called at the end of the subevts loop. Not (necessarily) able to access SubEvents */
+ StatusCode mergeEvent();
+
+ /** When being run from MM_Digitizer, this method is called during the event loop */
+
+ /** alternative interface which uses the PileUpMergeSvc to obtain
+ all the required SubEvents. */
+ virtual StatusCode processAllSubEvents();
+
+ /** Just calls processAllSubEvents - leaving for back-compatibility
+ (IMuonDigitizationTool) */
+
+ StatusCode digitize();
+
+ /** Finalize */
+ StatusCode finalize();
+
+ /** accessors */
+ ServiceHandle<IAtRndmGenSvc> getRndmSvc() const { return m_rndmSvc; } // Random number service
+ CLHEP::HepRandomEngine *getRndmEngine() const { return m_rndmEngine; } // Random number engine used
+
+ void setAmplification(const double amplification) {
+ m_amplification = amplification;
+ }
+
+ void set (const double bunchTime);
+
+private:
+
+ ServiceHandle<StoreGateSvc> m_sgSvc;
+ ActiveStoreSvc* m_activeStore;
+
+ /** Record MmDigitContainer and MuonSimDataCollection */
+ StatusCode recordDigitAndSdoContainers();
+
+ MmDigitContainer* m_digitContainer;
+ MuonSimDataCollection* m_sdoContainer;
+ MmSortedHitVector m_hits;
+
+ const MmIdHelper* m_idHelper;
+ MicromegasHitIdHelper* muonHelper;
+
+ const MuonGM::MuonDetectorManager* m_MuonGeoMgr;
+
+ ToolHandle <IMM_DigitizationTool> m_digitTool;
+ std::list<GenericMuonSimHitCollection*> m_MMHitCollList;
+ bool checkMMSimHit(const GenericMuonSimHit& /* hit */ ) const;
+
+ //CONFIGURATION
+ bool m_validationSetup;
+ double m_energyThreshold;
+ bool m_saveInternalHistos;
+ bool m_checkMMSimHits;
+ bool m_useTof;
+ bool m_useAttenuation;
+ bool m_useProp;
+
+ //pile-up
+ TimedHitCollection<GenericMuonSimHit>* m_thpcMM; // the hits
+
+ ///////////////////////////////////////////////////////////////////
+ // Access to the event methods:
+ ///////////////////////////////////////////////////////////////////
+ // Get next event and extract collection of hit collections:
+ StatusCode getNextEvent();
+ StatusCode doDigitization();
+
+ void fillMaps(const GenericMuonSimHit * mmHit, const Identifier digitId, const double driftR);
+ int digitizeTime(double time) const;
+ bool outsideWindow(double time) const; // default +-50...
+
+ //TIMING SCHEME
+ bool m_useTimeWindow;
+ double m_inv_c_light;
+ double m_DiffMagSecondMuonHit;
+
+ //TDC ELECTRONICS
+ double m_ns2TDC;
+ double m_resTDC;
+
+ // RandomGenerators
+ // double m_FlatDist;
+ // double m_GaussDist;
+ // double m_PoissonDist;
+ // double m_GammaDist;
+ double m_Polya;
+ double m_timeWindowLowerOffset;
+ double m_timeWindowUpperOffset;
+ double m_bunchTime;
+ double * m_sprob;
+ double m_amplification;
+
+ // StripsResponse stuff...
+ StripsResponse *m_StripsResponse;
+ float m_qThreshold, m_diffusSigma, m_LogitundinalDiffusSigma, m_driftGap, m_driftVelocity, m_crossTalk1, m_crossTalk2;
+ float m_qThresholdForTrigger;
+
+ // ElectronicsResponse stuff...
+ ElectronicsResponse *m_ElectronicsResponse;
+ float m_alpha;// power of responce function
+ float m_RC ;// time constant of responce function
+ float m_electronicsThreshold; // threshold "Voltage" for histoBNL
+ TFile *m_file;
+ TTree *m_ntuple;
+ TH1I *m_AngleDistr, *m_AbsAngleDistr, *m_ClusterLength2D, *m_ClusterLength, *m_gasGap, *m_gasGapDir ;
+
+ int m_n_Station_side, m_n_Station_eta, m_n_Station_phi, m_n_Station_multilayer, m_n_Station_layer, m_n_hitStripID, m_n_StrRespTrg_ID, m_n_strip_multiplicity, m_n_strip_multiplicity_2;
+ int exitcode, m_n_hitPDGId;
+ double m_n_hitOnSurface_x, m_n_hitOnSurface_y, m_n_hitDistToChannel, m_n_hitIncomingAngle,m_n_StrRespTrg_Time, m_n_hitIncomingAngleRads, m_n_hitKineticEnergy, m_n_hitDepositEnergy;
+ float tofCorrection, bunchTime, globalHitTime;
+ std::vector<int> m_n_StrRespID;
+ std::vector<float> m_n_StrRespCharge, m_n_StrRespTime;
+
+ // vector <TH1F*> m_histoBNL;
+
+ // private methods need to compute the induce charge on the strip
+ double qStrip (const int & nElectrons, const double & gammaDist) const;
+ double qStripR (const double x) const;
+ double qStripR (const double x, const std::string stationType) const;
+ double qStripPhi (const double x, const std::string stationType) const;
+ double fparamPhi (const double x, const int N, const double * p) const;
+ // double getDriftTime(const MMReadoutElement* descriptor, const Amg::Vector3D& pos) const;
+
+
+protected:
+
+ PileUpMergeSvc *m_mergeSvc; // Pile up service
+ std::string m_inputObjectName; // name of the input objects
+ std::string m_outputObjectName; // name of the output digits
+ std::string m_outputSDOName; // name of the output SDOs
+
+ 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
+ int m_maskMultiplet;
+
+};
+/*******************************************************************************/
+inline bool MmDigitizationTool::outsideWindow(double driftTime) const {
+ double time = driftTime; //m_bunchTime is included already....updated one..
+ // double time = driftTime+m_bunchTime;
+ return time < m_timeWindowLowerOffset || time > m_timeWindowUpperOffset;
+}
+/*******************************************************************************/
+inline double MmDigitizationTool::qStrip (const int & nElectrons, const double& gammaDist) const {
+
+ // find the charge on the wire
+
+ // double amplification = 0.58e5;
+ double amplification = m_amplification;
+ double stripCharge = 0.0;
+
+ for (int i=0; i<nElectrons; i++) {
+ double RNDpol = 0.0;
+ if (m_Polya > -1.0) {
+ RNDpol = gammaDist/(1.0+m_Polya);
+ }
+ stripCharge += amplification*RNDpol;
+ }
+ return stripCharge;
+}
+/*******************************************************************************/
+inline void MmDigitizationTool::set (const double bunchTime) {m_bunchTime = bunchTime;}
+
+#endif // MmDigitizationTool
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmElectronicsToolInput.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmElectronicsToolInput.h
new file mode 100644
index 0000000000000000000000000000000000000000..d9f3e2b9fe9a726f994c19a5cb641cab90aa84e2
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmElectronicsToolInput.h
@@ -0,0 +1,28 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_MMELECTRONICSTOOLINPUT_H
+#define MM_DIGITIZATION_MMELECTRONICSTOOLINPUT_H
+
+/*******************************************************************************/
+class MmElectronicsToolInput {
+
+ public:
+
+ MmElectronicsToolInput(std::vector<int> NumberOfStripsPos, std::vector<float> chipCharge, std::vector<float> chipTime ) :
+ m_NumberOfStripsPos(NumberOfStripsPos),m_chipCharge(chipCharge),m_chipTime(chipTime) {}
+
+ ~MmElectronicsToolInput() {}
+
+ std::vector<int> NumberOfStripsPos() const { return m_NumberOfStripsPos; }
+ std::vector<float> chipCharge() const { return m_chipCharge; }
+ std::vector<float> chipTime() const { return m_chipTime; }
+
+ private:
+ std::vector<int> m_NumberOfStripsPos;
+ std::vector<float> m_chipCharge;
+ std::vector<float> m_chipTime;
+};
+/*******************************************************************************/
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmSortedHitVector.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmSortedHitVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbfb64af3c9c6275cec04c25f6650adb027b986f
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MmSortedHitVector.h
@@ -0,0 +1,52 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MM_DIGITIZATION_MMSORTEDHITVECTOR_H
+#define MM_DIGITIZATION_MMSORTEDHITVECTOR_H
+/*******************************************************************************/
+class micromegas_hit_info {
+ public:
+ micromegas_hit_info(Identifier i,double t,int ch,double strP,const TimedHitPtr<GenericMuonSimHit>* aHit) :
+ id(i),time(t),charge(ch),stripPos(strP),simhit(aHit) {}
+ micromegas_hit_info(Identifier i,double t,int ch) : id(i),time(t),charge(ch) {}
+ micromegas_hit_info() :time(0.),charge(0) {}
+ Identifier id;
+ double time;
+ int charge;
+ double stripPos;
+ const TimedHitPtr<GenericMuonSimHit>* simhit;
+ bool operator < (const micromegas_hit_info& aInfo) const
+ {
+ if(id < aInfo.id)
+ return true;
+ else if(id == aInfo.id)
+ return time < aInfo.time;
+ return false;
+ }
+};
+/*******************************************************************************/
+// instead of using a map we use a sorted vector to temporary store the hits.
+typedef std::vector<micromegas_hit_info> HitVector;
+typedef HitVector::iterator HitIt;
+/*******************************************************************************/
+class MmSortedHitVector : public HitVector {
+ public:
+ void insert(const micromegas_hit_info& hit);
+ void sort();
+ bool isSorted() { return m_isSorted; }
+ private:
+ bool m_isSorted;
+};
+/*******************************************************************************/
+inline void MmSortedHitVector::sort() {
+ std::sort(HitVector::begin(),HitVector::end());
+ m_isSorted = true;
+}
+/*******************************************************************************/
+inline void MmSortedHitVector::insert(const micromegas_hit_info& hit) {
+ push_back(hit);
+ m_isSorted = false;
+}
+/*******************************************************************************/
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/StripsResponse.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/StripsResponse.h
new file mode 100644
index 0000000000000000000000000000000000000000..23e75a13a981ebb527cc15a8dd51d738d5445cea
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/StripsResponse.h
@@ -0,0 +1,179 @@
+/* -*- C++ -*- */
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef StripsResponse_H
+#define StripsResponse_H
+/** @class StripsResponse
+
+// ------------
+// Authors:
+// Iakovidis George <george.iakovidis@cern.ch>
+// Karakostas Konstantinos <konstantinos.karakostas@cern.ch>
+// Leontsinis Stefanos <stefanos.leontsinis@cern.ch>
+// Nektarios Chr. Benekos <nbenekos@cern.ch>
+//////////////////////////////////////////////////////////////////////////////
+
+Comments to be added here...
+
+*/
+
+#include "GaudiKernel/AlgFactory.h"
+#include "GaudiKernel/IToolSvc.h"
+#include "GaudiKernel/Service.h"
+#include "AthenaKernel/MsgStreamMember.h"
+#include "GaudiKernel/StatusCode.h"
+
+/// ROOT
+#include <TROOT.h>
+#include <TFile.h>
+#include <TF1.h>
+#include <TRandom.h>
+#include <TRandom3.h>
+#include <TCanvas.h>
+#include <TMath.h>
+
+/// Projects
+#include "MM_Digitization/MmDigitToolInput.h"
+#include "MM_Digitization/MmDigitToolOutput.h"
+
+/// STD'S
+#include <algorithm>
+#include <cmath>
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+#include <utility>
+#include <string>
+#include <sstream>
+#include <sys/stat.h>
+
+using std::vector;
+using std::cout;
+using std::endl;
+
+/// ROOT Classed
+class TF1;
+class TRandom;
+class TRandom3;
+
+class ElectronicsResponse;
+class MmDigitToolInput;
+class MmDigitToolOutput;
+
+class StripsResponse {
+
+private:
+
+ /** qThreshold=2e, we accept a good strip if the charge is >=2e */
+ float qThreshold;
+ /** // 0.350/10 diffusSigma=transverse diffusion (350 μm per 1cm ) for 93:7 @ 600 V/cm, according to garfield */
+ float diffusSigma;
+ float LogitundinalDiffusSigma;
+ float pitch;
+ /** //pitch=0.500 properties of the micromegas ToDo: to be reviewed */
+ float stripwidth;
+ /** crosstalk of neighbor strips, it's 15% */
+ float crossTalk1;//0.10; //
+ /** // crosstalk of second neighbor strips, it's 6% */
+ float crossTalk2;//0.03;
+ /** // (degrees) Magnetic Field 0.5 T */
+ float Lorentz_Angle;
+
+ /// ToDo: random number from custom functions
+ TF1 *polya, *conv_gaus;
+
+public :
+
+ float driftGap;
+ /** //0.050 drift velocity in [mm/ns], driftGap=5 mm +0.128 mm (the amplification gap) */
+ float driftVelocity;
+
+ StripsResponse();
+
+ virtual ~StripsResponse();
+ MmDigitToolOutput GetResponceFrom(const MmDigitToolInput & digiInput);
+
+ void calculateResponceFrom (const float & hitx, const int & stripOffest, const float & thetaD, const int & stripMaxID);
+ void initializationFrom ();
+ void testRandomNumGens (UInt_t sd);
+ void writeHistos();
+ void initHistos ();
+ void clearValues ();
+ void initFunctions ();
+ void whichStrips(const float & hitx, const int & stripOffest, const float & thetaD, const int & stripMaxID);
+
+ inline void set_qThreshold (float val) { qThreshold = val; };
+ inline void set_diffusSigma (float val) { diffusSigma = val; };
+ inline void set_LogitundinalDiffusSigma (float val) { LogitundinalDiffusSigma = val; };
+ inline void set_driftVelocity (float val) { driftVelocity = val; };
+ inline void set_crossTalk1 (float val) { crossTalk1 = val; };
+ inline void set_crossTalk2 (float val) { crossTalk2 = val; };
+ inline void set_driftGap (float val) {driftGap = val;};
+ inline void set_stripWidth (float val) {stripwidth = val;};
+
+ float get_stripWidth () const { return stripwidth ;};
+ float get_qThreshold () const { return qThreshold ;};
+ float get_driftGap () const { return driftGap ;};
+ float get_driftVelocity () const { return driftVelocity;};
+
+/*
+ Coverity defects
+ float get_tMinFirstPeak () const { return tMinFirstPeak;};
+ float get_tMinIntegration () const { return tMinIntegration;};
+ float get_tminIntegrationAboveThreshold () const { return tminIntegrationAboveThreshold;};
+*/
+
+ vector <float> get_tStripElectronicsAbThr() const { return tStripElectronicsAbThr;};
+ vector <float> get_qStripElectronics() const { return qStripElectronics;};
+ vector <float> get_finaltStripNoSlewing() const { return finaltStripNoSlewing;};
+ vector <float> get_finalqStrip() const { return finalqStrip;};
+ vector <float> get_finaltStrip() const { return finaltStrip;};
+ vector <int> get_nStripElectronics() const { return nStripElectronics;};
+ vector <int> get_finalNumberofStrip() const { return finalNumberofStrip;};
+
+ vector <int> finalNumberofStrip;
+ vector <int> nStripElectronics;
+ vector <float> finalqStrip;
+ vector <float> finaltStrip;
+ vector <float> finaltStripNoSlewing;
+ vector <float> tStripElectronicsAbThr;
+ vector <float> qStripElectronics;
+
+ int nstrip;
+ float temp_polya;
+ float polya_theta;
+ float numberOfElectrons;
+/*
+ Coverity defects
+ float tMinFirstPeak;
+ float tMinIntegration;
+ float tminIntegrationAboveThreshold;
+*/
+
+// whichStrips()
+ int dimClusters; //dimClusters=total number of collisions
+ int MaxPrimaryIons;
+ vector <int> stripNumber;
+ double pt,xx, xxDiffussion, yy, yyDiffussion ;
+ int primaryion;
+ vector <int> firstq;
+ float lmean;
+ vector <float> qstrip;
+ vector <float> cntTimes;
+ vector <float> tStrip;
+ vector <float> qStrip;
+ vector <float> time; //Drift velocity [mm/ns]
+ vector <int> numberofStrip;
+
+ vector <float> clusterelectrons;
+ vector <float> l;
+ float totalelectrons;
+ float ll;
+
+ TRandom * gRandom_loc;
+ TRandom3 *randomNum;
+};
+#endif
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/Makefile b/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..c2c9dbbed79a7894e694cf7a0c7d1bf0756334c3
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/Makefile
@@ -0,0 +1,4 @@
+include $(CMTROOT)/src/Makefile.header
+
+include $(CMTROOT)/src/constituents.make
+
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/requirements b/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..0bdcdf6ecec08012286acb0e4a932452a7d27335
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/cmt/requirements
@@ -0,0 +1,52 @@
+package MM_Digitization
+
+author Nektarios Benekos <nectarios.benekos@cern.ch>
+manager Nektarios Benekos <nectarios.benekos@cern.ch>
+
+public
+use AtlasPolicy AtlasPolicy-*
+use AthenaKernel AthenaKernel-* Control
+use AthenaBaseComps AthenaBaseComps-* Control
+use PileUpTools PileUpTools-* Control
+use AtlasCLHEP AtlasCLHEP-* External
+use GaudiInterface GaudiInterface-* External
+use AtlasROOT AtlasROOT-* External
+use Identifier Identifier-* DetectorDescription
+use HitManagement HitManagement-* Simulation
+use MuonDigToolInterfaces MuonDigToolInterfaces-* MuonSpectrometer/MuonDigitization
+use MuonSimEvent MuonSimEvent-* MuonSpectrometer
+
+private
+use StoreGate StoreGate-* Control
+use AtlasAIDA AtlasAIDA-* External
+use AtlasHepMC AtlasHepMC-* External
+use GeneratorObjects GeneratorObjects-* Generators
+use MuonDigitContainer MuonDigitContainer-* MuonSpectrometer
+use MuonIdHelpers MuonIdHelpers-* MuonSpectrometer
+use MuonSimData MuonSimData-* MuonSpectrometer
+use PathResolver PathResolver-* Tools
+use TrkDetDescrUtils TrkDetDescrUtils-* Tracking/TrkDetDescr
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
+use EventInfo EventInfo-* Event
+use EventInfoMgt EventInfoMgt-* Event
+###use MuonCondInterface MuonCondInterface-* MuonSpectrometer/MuonConditions/MuonCondGeneral
+use MuonReadoutGeometry MuonReadoutGeometry-* MuonSpectrometer/MuonDetDescr
+#use MuonDigiExample MuonDigiExample-* MuonSpectrometer/MuonDigitization
+use AtlasCLHEP_RandomGenerators AtlasCLHEP_RandomGenerators-* Simulation/Tools
+
+
+
+
+public
+library MM_Digitization *.cxx components/*.cxx
+apply_pattern component_library
+apply_pattern declare_joboptions files="-s=../share *.txt *.py"
+apply_pattern declare_python_modules files="*.py"
+#apply_pattern declare_runtime_extras extras="*.rt"
+
+apply_tag ROOTBasicLibs
+apply_tag ROOTMathLibs
+apply_tag ROOTSTLDictLibs
+apply_tag ROOTGraphicsLibs
+apply_tag ROOTTableLibs
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfig.py b/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..bfb13c2cbc15ccb107fd5a1b6b81f5a7c3f6377b
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfig.py
@@ -0,0 +1,65 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+#
+# Import MM_Digitization job properties
+#
+from Digitization.DigitizationFlags import jobproperties
+from AthenaCommon.BeamFlags import jobproperties
+from AthenaCommon import CfgMgr
+from AthenaCommon.AppMgr import ToolSvc, ServiceMgr
+
+# The earliest bunch crossing time for which interactions will be sent
+# to the MdtDigitizationTool.
+def MM_FirstXing():
+ return -250
+
+# The latest bunch crossing time for which interactions will be sent
+# to the MdtDigitizationTool.
+def MM_LastXing():
+ return 200
+
+def MmDigitizationTool(name="MmDigitizationTool",**kwargs):
+ kwargs.setdefault("RndmSvc", jobproperties.Digitization.rndmSvc() )
+ # set rndm seeds
+ mmRndm = kwargs.setdefault("RndmEngine","MM_Digitization")
+ jobproperties.Digitization.rndmSeedList.addSeed(mmRndm, 49261510, 105132394 )
+ if jobproperties.Digitization.doXingByXingPileUp():
+ kwargs.setdefault("FirstXing", MM_FirstXing() ) # this should match the range for the MM in Digitization/share/MuonDigitization.py
+ kwargs.setdefault("LastXing", MM_LastXing() ) # this should match the range for the MM in Digitization/share/MuonDigitization.py
+ kwargs.setdefault("CheckSimHits", True)
+ kwargs.setdefault("InputObjectName", "MicromegasSensitiveDetector")
+ kwargs.setdefault("OutputObjectName", "MM_DIGITS")
+ kwargs.setdefault("OutputSDOName", "MM_SDO")
+ kwargs.setdefault("ValidationSetup", False)
+ ## Strip response
+ kwargs.setdefault("qThreshold", 0.001)
+ kwargs.setdefault("DiffusSigma", 0.036)
+ kwargs.setdefault("LogitundinalDiffusSigma", 0.019)
+ kwargs.setdefault("DriftVelocity", 0.047)
+ kwargs.setdefault("crossTalk1", 0.0)
+ kwargs.setdefault("crossTalk2", 0.0)
+ ## Chip/electronics response
+ kwargs.setdefault("alpha", 2.5)
+ kwargs.setdefault("RC", 20.)
+
+ return CfgMgr.MmDigitizationTool(name,**kwargs)
+ #return CfgMgr.MM_PileUpTool(name,**kwargs)
+ #else:
+ #return CfgMgr.MdtDigitizationTool(name,**kwargs)
+
+def getMMRange(name="MMRange", **kwargs):
+ # bunch crossing range in ns
+ kwargs.setdefault('FirstXing', MM_FirstXing() )
+ kwargs.setdefault('LastXing', MM_LastXing() )
+ kwargs.setdefault('CacheRefreshFrequency', 1.0 ) #default 0 no dataproxy reset
+ kwargs.setdefault('ItemList', ["GenericMuonSimHitCollection#MicromegasSensitiveDetector"] )
+ return CfgMgr.PileUpXingFolder(name, **kwargs)
+
+
+def MM_Response_DigitTool(name="MM_Response_DigitTool",**kwargs):
+ kwargs.setdefault("RndmSvc", jobproperties.Digitization.rndmSvc())
+ mmRndm = kwargs.setdefault("RndmEngine", "MMResponse")
+ jobproperties.Digitization.rndmSeedList.addSeed(mmRndm, 49261510,105132394 )
+ return CfgMgr.MM_Response_DigitTool(name,**kwargs)
+
+
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfigDb.py b/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfigDb.py
new file mode 100644
index 0000000000000000000000000000000000000000..b99cf5c4996754bd3e7e1a4e6d40d8d3c14af8ab
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/python/MM_DigitizationConfigDb.py
@@ -0,0 +1,6 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon.CfgGetter import addTool
+addTool("MM_Digitization.MM_DigitizationConfig.MmDigitizationTool","MmDigitizationTool")
+addTool("MM_Digitization.MM_DigitizationConfig.MM_Response_DigitTool","MM_Response_DigitTool")
+addTool("MM_Digitization.MM_DigitizationConfig.getMMRange", "MMRange")
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/share/MicroMegas_Digitization_jobOptions.py b/MuonSpectrometer/MuonDigitization/MM_Digitization/share/MicroMegas_Digitization_jobOptions.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ee311ee3a86a270acd2f1fb7d1a55ffbcfa6704
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/share/MicroMegas_Digitization_jobOptions.py
@@ -0,0 +1,31 @@
+include.block ("MM_Digitization/MicroMegas_Digitization_jobOptions.py")
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+
+from AthenaCommon import CfgGetter
+from Digitization.DigitizationFlags import jobproperties
+
+job += CfgGetter.getAlgorithm("MM_Digitizer/MM_Digitizer", tryDefaultConfigurable=True)
+from MM_Digitization.MM_DigitizationConf import MmDigitizationTool
+MmDigitizationTool = MmDigitizationTool("MmDigitizationTool",
+ RndmSvc = jobproperties.Digitization.rndmSvc.get_Value(),
+ RndmEngine = "MM_Digitization",
+ InputObjectName = "MicromegasSensitiveDetector",
+ OutputObjectName = "MM_DIGITS",
+ OutputSDOName = "MM_SDO",
+ CheckSimHits = True,
+ ## Strip response
+ qThreshold = 0.00000001,
+ DiffusSigma = 0.036,
+ LogitundinalDiffusSigma = 0.019,
+ DriftVelocity = 0.047,
+ crossTalk1 = 0.0,
+ crossTalk2 = 0.0,
+ ## Chip/electronics response
+ alpha = 2.5,
+ RC = 20.,
+ SaveInternalHistos = True,
+ qThresholdForTrigger = 2,
+ ValidationSetup = False,
+ EnergyThreshold = 50,
+ OutputLevel = DEBUG)
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/share/postInclude.MicroMegasDigitization.py b/MuonSpectrometer/MuonDigitization/MM_Digitization/share/postInclude.MicroMegasDigitization.py
new file mode 100644
index 0000000000000000000000000000000000000000..5c64715f70c7d4a56f7e4d4644169e94e070f2cc
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/share/postInclude.MicroMegasDigitization.py
@@ -0,0 +1,36 @@
+#########################################################
+#
+# MM_Digitization/postInclude.MicroMegasDigitization.py
+# Nektarios Chr. Benekos
+# This job option should be added via the postInclude
+# command line argument.
+#########################################################
+include.block ("MM_Digitization/postInclude.MicroMegasDigitization.py")
+from Digitization.DigitizationFlags import digitizationFlags
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+if digitizationFlags.doXingByXingPileUp():
+ MMdigitization = job.PileUpToolsAlg.PileUpTools[ "MmDigitizationTool" ]
+else:
+ from MM_Digitization.MM_DigitizationConf import MmDigitizationTool as MMdigitization
+ MMdigitization.RndmSvc = jobproperties.Digitization.rndmSvc.get_Value()
+ MMdigitization.RndmEngine = "MM_Digitization"
+ MMdigitization.InputObjectName = "MicromegasSensitiveDetector"
+ MMdigitization.OutputObjectName = "MM_DIGITS"
+ MMdigitization.OutputSDOName = "MM_SDO"
+ MMdigitization.CheckSimHits = True
+ MMdigitization.ValidationSetup = False
+ ## Strip response
+ MMdigitization.qThreshold = 0.00000000001
+ MMdigitization.qThresholdForTrigger = 2.
+ MMdigitization.DiffusSigma = 0.036
+ MMdigitization.LogitundinalDiffusSigma = 0.019
+ MMdigitization.DriftVelocity = 0.047
+ MMdigitization.crossTalk1 = 0.0
+ MMdigitization.crossTalk2 = 0.0
+ ## Chip/electronics response
+ MMdigitization.alpha = 2.5
+ MMdigitization.RC = 20.
+ MMdigitization.SaveInternalHistos = True
+ MMdigitization.EnergyThreshold = 50
+ MMdigitization.OutputLevel = DEBUG
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/ElectronicsResponse.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/ElectronicsResponse.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2ede48af2e841c892df15fdd2a2f62d723722b04
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/ElectronicsResponse.cxx
@@ -0,0 +1,187 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * ElectronicsResponse.cxx
+ * MC for micromegas athena integration
+ *
+ **/
+//~
+/// PROJECTS
+#include "GaudiKernel/MsgStream.h"
+#include "MM_Digitization/ElectronicsResponse.h"
+
+/*******************************************************************************/
+ElectronicsResponse::ElectronicsResponse()
+{
+ //stripObject = 0;
+ stripObject = new StripsResponse();
+ m_StripResponse_qThreshold = stripObject->get_qThreshold();
+ m_StripResponse_driftGap = stripObject->get_driftGap();
+ m_StripResponse_driftVelocity = stripObject->get_driftVelocity();
+// if(stripObject) delete stripObject;
+// clearValues();
+ alpha = 2.5;
+ RC = 20.0;
+ electronicsThreshold = (m_StripResponse_qThreshold * ( TMath::Power(alpha,alpha)*TMath::Exp(-alpha)) ) ;
+
+ //---------------------------------------------------------
+
+ intFn = new TF1("intFn","[3]*TMath::Power((x-[2])/[1],[0])*TMath::Exp(-(x-[2])/[1])", 0, 3*( m_StripResponse_driftGap/m_StripResponse_driftVelocity ) );
+ intFn->SetParameter( 0, alpha);
+ intFn->SetParameter( 1, RC);
+}
+/*******************************************************************************/
+void ElectronicsResponse::clearValues()
+{
+ tMinFromIntegration = -1;
+ tminFromIntegrationAboveThreshold = -1;
+ tMinFromFirstPeak = -1;
+ tStripElectronicsAbThr.clear();
+ qStripElectronics.clear();
+ nStripElectronics.clear();
+}
+/*******************************************************************************/
+MmDigitToolOutput ElectronicsResponse::GetResponceFrom(const MmElectronicsToolInput & digiInput)
+{
+ clearValues();
+ bnlResponceFunction(digiInput.NumberOfStripsPos(), digiInput.chipCharge(), digiInput.chipTime() );
+ /// ToDo: include loop for calculating Trigger study vars
+ // MmDigitToolOutput(bool hitWasEff, std::vector <int> strpos, std::vector<float> time, std::vector<int> charge, int strTrig, float strTimeTrig ):
+ MmDigitToolOutput tmp(true, nStripElectronics, tStripElectronicsAbThr, qStripElectronics, 5, 0.3);
+ return tmp;
+}
+/*******************************************************************************/
+void ElectronicsResponse::bnlResponceFunction(const vector <int> & numberofStrip, const vector <float> & qStrip, const vector <float> & tStrip){
+
+ //int minStrip=1000000, maxStrip=-10000000;
+ //tStripElectronicsAbThr.clear();
+ //nStripElectronics.clear();
+ //qStripElectronics.clear();
+
+ //stripObject = new StripsResponse();
+ //cout << "stripObject->get_qThreshold()" << stripObject->get_qThreshold() << endl;
+ //cout << "m_StripResponse_driftGap" << m_StripResponse_driftGap << endl;
+ //cout << "m_StripResponse_driftVelocity" << m_StripResponse_driftVelocity << endl;
+
+
+ for (unsigned int ii = 0; ii < numberofStrip.size(); ii++) {
+
+ // cout << "----------for strip number " << numberofStrip.at(ii) << "---------" << endl;
+ //cout << "time max: " << 3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity) << endl;
+ //cout << "charge at (0,50,100,200): " << intFn->Eval(0,0,0) << ", " << intFn->Eval(50,0,0) << ", " << intFn->Eval(100,0,0) << ", " << intFn->Eval(200,0,0) << endl;
+
+ //find min and max times for each strip:
+
+ //maxChargeThisStrip = intFn->GetMaximum(0, 3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity), 1.0E-10, 100, false);
+ //maxChargeThisStrip = intFn->GetMaximum(50,300);
+ //maxChargeThisStrip = intFn->GetMaximum(tStrip.at(ii),3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity) );
+
+
+ //cout << "maxChargeThisStrip(with tmin = tstrip): " << maxChargeThisStrip << endl;
+
+ double maxChargeThisStrip = 0;
+ double maxChargeLeftNeighbor = 0;
+ double maxChargeRightNeighbor = 0;
+
+ Athena::MsgStreamMember log("ElectronicsResponse::bnlResponceFunction");
+ // find the maximum charge:
+ if ( ii > 0 ) {
+ log << MSG::DEBUG << "for Left neighbor: tStrip.at(ii-1) "<< tStrip.at(ii-1) << " qStrip.at(ii-1) " << qStrip.at(ii-1) << endreq;
+ intFn->SetParameter( 2, tStrip.at(ii-1));
+ intFn->SetParameter( 3, qStrip.at(ii-1));
+ maxChargeLeftNeighbor = intFn->GetMaximum(tStrip.at(ii-1),3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity));
+ }
+
+ if ( ii+1 < numberofStrip.size() ) {
+ log << MSG::DEBUG << "for Right neighbor: tStrip.at(ii+1) "<< tStrip.at(ii+1) << " qStrip.at(ii+1) " << qStrip.at(ii+1) << endreq;
+ intFn->SetParameter( 2, tStrip.at(ii+1));
+ intFn->SetParameter( 3, qStrip.at(ii+1));
+ maxChargeRightNeighbor = intFn->GetMaximum(tStrip.at(ii+1),3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity));
+ }
+
+ if ( ii > 0 ) {
+ log << MSG::DEBUG << "for This strip: tStrip.at(ii) "<< tStrip.at(ii) << " qStrip.at(ii) " << qStrip.at(ii) << endreq;
+ intFn->SetParameter( 2, tStrip.at(ii));
+ intFn->SetParameter( 3, qStrip.at(ii));
+ maxChargeThisStrip = intFn->GetMaximum(tStrip.at(ii),3*(m_StripResponse_driftGap/m_StripResponse_driftVelocity));
+ }
+
+ log << MSG::DEBUG << "Check if a strip or its neighbors were above electronics threshold (" << electronicsThreshold << "): maxChargeLeftNeighbor: " << maxChargeLeftNeighbor << ", maxChargeRightNeighbor: " << maxChargeRightNeighbor << ", maxChargeThisStrip: " << maxChargeThisStrip << endreq;
+
+ // Look at strip if it or its neighbor was above threshold:
+ if ( maxChargeLeftNeighbor > electronicsThreshold || maxChargeRightNeighbor > electronicsThreshold || maxChargeThisStrip > electronicsThreshold ) {
+
+ intFn->SetParameter( 2, tStrip.at(ii));
+ intFn->SetParameter( 3, qStrip.at(ii));
+
+ float localPeak = 0;
+ float localPeakt = 0;
+ float localPeakq = 0;
+
+ float stepSize = 0.1; //ns(?) step size corresponding to VMM discriminator
+
+ for (int jj = tStrip.at(ii); jj < (3*( m_StripResponse_driftGap/m_StripResponse_driftVelocity))/stepSize; jj++) {
+
+ //cout << "start looking for local peak: " << endl;
+
+
+ float thisStep = jj*stepSize;
+ float nextStep = (jj+1)*stepSize;
+ float oneAfterStep = (jj+2)*stepSize;
+
+ //cout << "intFn->Eval(thisStep,0,0): " << intFn->Eval(thisStep,0,0) << endl;
+ //cout << "intFn->Derivative(thisStep,0,0.001): " << intFn->Derivative(thisStep,0,0.001) << endl;
+
+ if (localPeak == 0 ) { //haven't found a local peak yet
+
+ //cout << "intFn->Eval(nextStep,0,0): " << intFn->Eval(nextStep,0,0) << endl;
+ //cout << "intFn->Derivative(nextStep,0,0.001): " << intFn->Derivative(nextStep,0,0.001) << endl;
+
+ //check if the charge for the next points is less than the current step and the derivative of the first point is positive or 0 and the next point is negative:
+ if ( ( intFn->Eval(thisStep,0,0) < intFn->Eval(nextStep,0,0) ) && ( intFn->Eval(nextStep,0,0) > intFn->Eval(oneAfterStep,0,0) ) &&
+ (intFn->Derivative(thisStep,0,0.001) > 0 && intFn->Derivative(oneAfterStep,0,0.001) < 0 ) ) {
+ // && (intFn->Derivative(thisStep,0,0.001)>0 && intFn->Derivative(nextStep,0,0.001)<=0 )
+ // && (intFn->Eval(oneAfterStep,0,0) < intFn->Eval(nextStep,0,0))
+
+
+ localPeak = 1;
+ localPeakt = nextStep;
+ localPeakq = intFn->Eval(nextStep,0,0);
+
+ log << MSG::DEBUG << "found a peak! for strip number: " << numberofStrip.at(ii) << " at time: " << nextStep << " with charge: " << intFn->Eval(nextStep,0,0) << endreq;
+
+ nStripElectronics.push_back(numberofStrip.at(ii));
+ tStripElectronicsAbThr.push_back(localPeakt);
+ qStripElectronics.push_back(localPeakq);
+
+ }
+ }
+ }
+ }
+ }
+
+ //---------------------------------------------------------
+ // /// DUMMY:
+ // tMinFromFirstPeak = tMinFromIntegration;
+
+ // nStripElectronics.push_back(ii);
+ // tStripElectronicsAbThr.push_back(tMinFromFirstPeak);
+ // if (maxOfIntegration>=electronicsThreshold) qStripElectronics.push_back(maxOfIntegration);
+ // else qStripElectronics.push_back(-50000.);
+ // }
+ // tMinFromFirstPeak = 50000.;
+ // for (size_t j = 0;j<tStripElectronicsAbThr.size();j++){
+ // if (tMinFromFirstPeak>tStripElectronicsAbThr.at(j)) tMinFromFirstPeak = tStripElectronicsAbThr.at(j);
+ // }
+ //delete spectrum;
+}///end of bnl responce function
+/*******************************************************************************/
+ElectronicsResponse::~ElectronicsResponse()
+{
+ if (stripObject) delete stripObject;
+ if (intFn) delete intFn;
+ clearValues();
+}
+/*******************************************************************************/
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Digitizer.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Digitizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7148332e1ac38d521b79b4c761989a1d5830cd55
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Digitizer.cxx
@@ -0,0 +1,44 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/// Gaudi External
+#include "AthenaBaseComps/AthMsgStreamMacros.h"
+#include "AthenaKernel/errorcheck.h"
+
+#include "MuonDigToolInterfaces/IMuonDigitizationTool.h"
+#include "MM_Digitization/MM_Digitizer.h"
+
+/*******************************************************************************/
+MM_Digitizer::MM_Digitizer(const std::string& name, ISvcLocator* pSvcLocator)
+ : AthAlgorithm(name, pSvcLocator),
+ m_digTool("MmDigitizationTool", this )
+{
+ declareProperty("MM_DigitizationTool", m_digTool , "AthAlgTool which performs the MicroMegas digitization");
+}
+
+/*******************************************************************************/
+MM_Digitizer::~MM_Digitizer() {}
+/*******************************************************************************/
+StatusCode MM_Digitizer::initialize() {
+ // intitialize store gate active store
+ if (m_digTool.retrieve().isFailure()) {
+ ATH_MSG_FATAL ("MM_Digitizer::Could not retrieve MM Digitization Tool!");
+ return StatusCode::FAILURE;
+ }
+ ATH_MSG_DEBUG ("MM_Digitizer::Retrieved MM Digitization Tool.");
+ return StatusCode::SUCCESS;
+
+}
+/*******************************************************************************/
+StatusCode MM_Digitizer::execute() {
+ ATH_MSG_DEBUG ("MM_Digitizer::in execute()");
+ return m_digTool->digitize();
+}
+/*******************************************************************************/
+StatusCode MM_Digitizer::finalize() {
+ ATH_MSG_INFO ("MM_Digitizer::finalize");
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Response_DigitTool.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Response_DigitTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..61d7597aa0307be402f2b772fd753687ab051a39
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_Response_DigitTool.cxx
@@ -0,0 +1,94 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GaudiKernel/ToolFactory.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "StoreGate/StoreGateSvc.h"
+
+#include "MM_Digitization/MmDigitToolInput.h"
+#include "MM_Digitization/MM_Response_DigitTool.h"
+
+
+#include "MuonIdHelpers/MmIdHelper.h"
+#include "MuonReadoutGeometry/MuonDetectorManager.h"
+
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include <iostream>
+#include <vector>
+
+using namespace MuonGM;
+using namespace std;
+/*******************************************************************************/
+MM_Response_DigitTool::MM_Response_DigitTool( const std::string& type,
+ const std::string& name,
+ const IInterface* parent )
+ : AlgTool(type,name,parent)
+ , m_rndmEngine(0)
+ , m_rndmEngineName("MuonDigitization")
+ , m_rndmSvc("AtRndmGenSvc", name )
+{
+ declareInterface<IMM_DigitizationTool>(this);
+ declareProperty("RndmSvc", m_rndmSvc, "Random Number Service used in Muon digitization" );
+ declareProperty("RndmEngine", m_rndmEngineName, "Random engine name");
+}
+/*******************************************************************************/
+MmDigitToolOutput MM_Response_DigitTool::digitize( /*const MmDigitToolInput& input*/ )
+{
+ vector<float> a, b;
+ vector<int> c;
+ MmDigitToolOutput output(false, c, b, a, 1, 1);
+ return output;
+}
+/*******************************************************************************/
+StatusCode MM_Response_DigitTool::initialize()
+{
+ MsgStream log(msgSvc(),name());
+
+ StoreGateSvc* detStore=0;
+ StatusCode status = serviceLocator()->service("DetectorStore", detStore);
+
+ if (status.isSuccess()) {
+ if(detStore->contains<MuonDetectorManager>( "Muon" )){
+ status = detStore->retrieve(m_muonGeoMgr);
+ if (status.isFailure()) {
+ if (log.level()<=MSG::FATAL) log << MSG::FATAL << "Could not retrieve MuonGeoModelDetectorManager!" << endreq;
+ return status;
+ }
+ else {
+ if (log.level()<=MSG::DEBUG) log << MSG::DEBUG << "MuonGeoModelDetectorManager retrieved from StoreGate."<< endreq;
+ //initialize the MdtIdHelper
+ m_idHelper = m_muonGeoMgr->mmIdHelper();
+ if (log.level()<=MSG::DEBUG) log << MSG::DEBUG << "MdtIdHelper: " << m_idHelper << endreq;
+ if(!m_idHelper) return status;
+ }
+ }
+ }
+
+ if (!m_rndmSvc.retrieve().isSuccess())
+ {
+ if (log.level()<=MSG::FATAL) log << MSG::FATAL << " Could not initialize Random Number Service" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ // getting our random numbers stream
+ if (log.level()<=MSG::DEBUG) log << MSG::DEBUG << "Getting random number engine : <" << m_rndmEngineName << ">" << endreq;
+ m_rndmEngine = m_rndmSvc->GetEngine(m_rndmEngineName);
+ if (m_rndmEngine==0) {
+ if (log.level()<=MSG::FATAL) log << MSG::FATAL << "Could not find RndmEngine : " << m_rndmEngineName << endreq;
+ return StatusCode::FAILURE;
+ }
+
+
+ initializeStrip();
+
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+bool MM_Response_DigitTool::initializeStrip(){
+ return true;
+}
+/*******************************************************************************/
+
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MmDigitizationTool.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MmDigitizationTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6c8bac9b6b3485d1ed1db1c0d4d1dde3d17cff3b
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MmDigitizationTool.cxx
@@ -0,0 +1,1069 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// MmDigitizationTool
+// ------------
+// Authors: Nektarios Chr. Benekos <nectarios.benekos@cern.ch>
+// Konstantinos Karakostas <Konstantinos.Karakostas@cern.ch>
+////////////////////////////////////////////////////////////////////////////////
+
+
+//Inputs
+#include "MuonSimData/MuonSimDataCollection.h"
+#include "MuonSimData/MuonSimData.h"
+
+//Outputs
+#include "MuonDigitContainer/MmDigitContainer.h"
+
+//MM digitization includes
+#include "MM_Digitization/MmDigitizationTool.h"
+#include "MM_Digitization/IMM_DigitizationTool.h"
+#include "MM_Digitization/MmDigitToolInput.h"
+#include "MuonSimEvent/MM_SimIdToOfflineId.h"
+
+#include "MuonReadoutGeometry/NSWenumeration.h"
+#include "MuonReadoutGeometry/NSWgeometry.h"
+
+//Gaudi - Core
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/AlgFactory.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "PathResolver/PathResolver.h"
+#include "AIDA/IHistogram1D.h"
+#include "EventInfo/TagInfo.h"
+#include "EventInfoMgt/ITagInfoMgr.h"
+
+//Geometry
+#include "MuonReadoutGeometry/MuonDetectorManager.h"
+#include "MuonReadoutGeometry/MMReadoutElement.h"
+#include "MuonReadoutGeometry/MuonChannelDesign.h"
+#include "MuonIdHelpers/MmIdHelper.h"
+#include "MuonSimEvent/MicromegasHitIdHelper.h"
+#include "TrkDetDescrUtils/GeometryStatics.h"
+#include "TrkEventPrimitives/LocalDirection.h"
+#include "TrkSurfaces/Surface.h"
+
+//Pile-up
+#include "PileUpTools/PileUpMergeSvc.h"
+
+//Truth
+#include "CLHEP/Units/PhysicalConstants.h"
+#include "GeneratorObjects/HepMcParticleLink.h"
+#include "HepMC/GenParticle.h"
+
+//Random Numbers
+#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "CLHEP/Random/RandGauss.h"
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGamma.h"
+#include "CLHEP/Random/RandPoisson.h"
+#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h"
+#include "CLHEP/Random/RandExponential.h"
+
+//ROOT
+#include "TH1.h"
+#include "TTree.h"
+#include "TFile.h"
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+
+using namespace MuonGM;
+
+/*******************************************************************************/
+MmDigitizationTool::MmDigitizationTool(const std::string& type, const std::string& name, const IInterface* parent)
+ : PileUpToolBase(type, name, parent)
+ , m_sgSvc("StoreGateSvc", name)
+ , m_digitContainer(NULL)
+ , m_sdoContainer(NULL)
+ , m_digitTool("MM_Response_DigitTool", this)
+ , m_validationSetup(false)
+ , m_energyThreshold(50.)
+ , m_saveInternalHistos(false)
+ , m_thpcMM(0)
+ , m_Polya(0.0)
+ , m_sprob(0)
+ , m_amplification(0.0)
+ , m_StripsResponse(0)
+ , m_qThreshold(0.001)
+ , m_diffusSigma(0.036)
+ , m_LogitundinalDiffusSigma(0.019)
+ , m_driftVelocity(0.047)
+ , m_crossTalk1(0.)
+ , m_crossTalk2(0.)
+ , m_qThresholdForTrigger(1.)
+ , m_ElectronicsResponse(0)
+ , m_file(0)
+ , m_ntuple(0)
+ , m_AngleDistr(0)
+ , m_AbsAngleDistr(0), m_ClusterLength2D(0), m_ClusterLength(0)
+ , m_mergeSvc(0)
+ , m_inputObjectName("")
+ , m_outputObjectName("")
+ , m_rndmSvc("AtRndmGenSvc", name )
+ , m_rndmEngine(0)
+ , m_rndmEngineName("MuonDigitization")
+{
+
+ declareInterface<IMuonDigitizationTool>(this);
+
+ declareProperty("MaskMultiplet", m_maskMultiplet = 0, "0: all, 1: first, 2: second, 3: both" );
+
+ //Random numbers service
+ declareProperty("RndmSvc", m_rndmSvc, "Random Number Service used in Muon digitization");
+ declareProperty("RndmEngine", m_rndmEngineName, "Random engine name");
+ declareProperty("DigitizationTool", m_digitTool, "Tool which handle the digitization process");
+ declareProperty("MCStore", m_sgSvc, "help");
+
+ //Object names
+ declareProperty("InputObjectName", m_inputObjectName = "MicromegasSensitiveDetector");
+ declareProperty("OutputObjectName", m_outputObjectName = "MM_DIGITS");
+ declareProperty("OutputSDOName", m_outputSDOName = "MM_SDO");
+
+ //Configurations
+ declareProperty("CheckSimHits", m_checkMMSimHits = true, "Control on the hit validity");
+
+ //Timing scheme
+ declareProperty("UseTimeWindow", m_useTimeWindow = true);
+ declareProperty("WindowLowerOffset", m_timeWindowLowerOffset = -50.);
+ declareProperty("WindowUpperOffset", m_timeWindowUpperOffset = +50.);
+ declareProperty("DiffMagSecondMuonHit",m_DiffMagSecondMuonHit = 0.1);
+
+ // electronics
+ declareProperty("ns2TDC", m_ns2TDC = 0.78125, "Conversion factor TDC/ns");
+ declareProperty("ResolutionTDC", m_resTDC = 0.5, "TDC resolution");
+
+ // Constants vars for the StripsResponse class
+ // qThreshold=2e, we accept a good strip if the charge is >=2e
+ declareProperty("qThreshold", m_qThreshold = 0.001);
+ // transverse diffusion (350 μm per 1cm ) for 93:7 @ 600 V/cm, according to garfield
+ declareProperty("DiffusSigma", m_diffusSigma = 0.036);
+ declareProperty("LogitundinalDiffusSigma", m_LogitundinalDiffusSigma = 0.019);
+ // 0.050 drift velocity in [mm/ns], driftGap=5 mm +0.128 mm (the amplification gap)
+ declareProperty("driftGap", m_driftGap = 5.128);
+ declareProperty("DriftVelocity", m_driftVelocity = 0.047);
+ // crosstalk of neighbor strips, it's 15%
+ declareProperty("crossTalk1", m_crossTalk1 = 0.0);
+ // crosstalk of neighbor strips, it's 6%
+ declareProperty("crossTalk2", m_crossTalk2 = 0.0);
+ declareProperty("qThresholdForTrigger", m_qThresholdForTrigger = 1.0);
+
+ // Constants vars for the ElectronicsResponse
+ declareProperty("alpha", m_alpha = 2.5);
+ declareProperty("RC", m_RC = 20.);
+ declareProperty("electronicsThreshold", m_electronicsThreshold = 0.000811174);
+
+ declareProperty("SaveInternalHistos", m_saveInternalHistos = false );
+ declareProperty("ValidationSetup", m_validationSetup = false );
+ declareProperty("EnergyThreshold", m_energyThreshold = 50., "Minimal energy to produce a PRD" );
+
+}
+/*******************************************************************************/
+// member function implementation
+//--------------------------------------------
+StatusCode MmDigitizationTool::initialize() {
+
+ StatusCode status(StatusCode::SUCCESS);
+ ATH_MSG_DEBUG ("MmDigitizationTool:: in initialize()") ;
+ ATH_MSG_DEBUG ( "Configuration MmDigitizationTool " );
+ ATH_MSG_DEBUG ( "RndmSvc " << m_rndmSvc );
+ ATH_MSG_DEBUG ( "RndmEngine " << m_rndmEngineName );
+ ATH_MSG_DEBUG ( "MCStore " << m_sgSvc );
+ ATH_MSG_DEBUG ( "DigitizationTool " << m_digitTool );
+ ATH_MSG_DEBUG ( "InputObjectName " << m_inputObjectName );
+ ATH_MSG_DEBUG ( "OutputObjectName " << m_outputObjectName );
+ ATH_MSG_DEBUG ( "OutputSDOName " << m_outputSDOName );
+ ATH_MSG_DEBUG ( "UseTimeWindow " << m_useTimeWindow );
+ ATH_MSG_DEBUG ( "CheckSimHits " << m_checkMMSimHits );
+ ATH_MSG_DEBUG ( "ns2TDC " << m_ns2TDC );
+ ATH_MSG_DEBUG ( "ResolutionTDC " << m_resTDC );
+ ATH_MSG_DEBUG ( "Threshold " << m_qThreshold );
+ ATH_MSG_DEBUG ( "DiffusSigma " << m_diffusSigma );
+ ATH_MSG_DEBUG ( "LogitundinalDiffusSigma" << m_LogitundinalDiffusSigma );
+ ATH_MSG_DEBUG ( "DriftVelocity " << m_driftVelocity );
+ ATH_MSG_DEBUG ( "crossTalk1 " << m_crossTalk1 );
+ ATH_MSG_DEBUG ( "crossTalk2 " << m_crossTalk2 );
+ ATH_MSG_INFO ( "ValidationSetup " << m_validationSetup );
+ ATH_MSG_INFO ( "EnergyThreshold " << m_energyThreshold );
+
+ // initialize random number generators
+ // then initialize the CSC identifier helper
+ // This method must be called before looping over the hits
+ // to digitize them using the method digitize_hit below
+ m_bunchTime = 0.0;
+
+ // initialize random number generators
+ // double average_int = 30; // average interactions per cm
+ // m_FlatDist = CLHEP::RandFlat::shoot(m_rndmEngine, 0.0,1.0);
+ // m_GaussDist = CLHEP::RandGauss::shoot(m_rndmEngine,0.0,1.0);
+ // m_GammaDist = CLHEP::RandGamma::shoot(m_rndmEngine, (1.0+m_Polya), 1.0);
+ // m_PoissonDist = CLHEP::RandPoisson::shoot(m_rndmEngine, average_int);
+
+
+ // initialize transient event store
+ if (m_sgSvc.retrieve().isFailure()) {
+ ATH_MSG_FATAL ( "Could not retrieve StoreGateSvc!" );
+ return StatusCode::FAILURE;
+ }
+ ATH_MSG_DEBUG ( "Retrieved StoreGateSvc." );
+
+ status = service("ActiveStoreSvc", m_activeStore);
+ if ( !status.isSuccess() ) {
+ ATH_MSG_FATAL ( "Could not get active store service" );
+ return status;
+ }
+
+ // initialize transient detector store and MuonGeoModel OR MuonDetDescrManager
+ StoreGateSvc* detStore=0;
+ m_MuonGeoMgr=0;
+ status = serviceLocator()->service("DetectorStore", detStore);
+ if (status.isSuccess()) {
+ if(detStore->contains<MuonGM::MuonDetectorManager>( "Muon" )){
+
+ status = detStore->retrieve(m_MuonGeoMgr);
+ 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_MuonGeoMgr->mmIdHelper();
+ if(!m_idHelper) return status;
+ ATH_MSG_DEBUG ( "Retrieved MmIdHelper " << m_idHelper );
+ }
+ }
+ }
+ else {
+ ATH_MSG_FATAL ( "Could not retrieve DetectorStore!" );
+ return status;
+ }
+
+ // 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 << "'" );
+ }
+
+ // check the output object name
+ if (m_outputObjectName=="") {
+ ATH_MSG_FATAL ( "Property OutputObjectName not set !" );
+ return StatusCode::FAILURE;
+ }
+ else {
+ ATH_MSG_DEBUG ( "Output digits: '" << m_outputObjectName << "'" );
+ }
+
+ //initialize the digit container
+ try{
+ m_digitContainer = new MmDigitContainer(m_idHelper->detectorElement_hash_max());
+ }
+ catch(std::bad_alloc){
+ ATH_MSG_FATAL ( "Could not create a new MicroMegas DigitContainer!" );
+ return StatusCode::FAILURE;
+ }
+ m_digitContainer->addRef();
+
+ //simulation identifier helper
+ muonHelper = MicromegasHitIdHelper::GetHelper();
+
+ //get the r->t conversion tool
+ status = m_digitTool.retrieve();
+ if( status.isFailure() ) {
+ ATH_MSG_FATAL("Could not retrieve digitization tool! " << m_digitTool);
+ return StatusCode::FAILURE;
+ }
+ else {
+ ATH_MSG_DEBUG("Retrieved digitization tool!" << m_digitTool);
+ }
+
+ // initialize inverse lightSpeed (c_light in m/s)
+ m_inv_c_light = 1./(CLHEP::c_light);
+
+ if (!m_rndmSvc.retrieve().isSuccess()) {
+ ATH_MSG_ERROR("Could not initialize Random Number Service");
+ }
+
+ // 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;
+ }
+
+ //locate the PileUpMergeSvc and initialize our local ptr
+ const bool CREATEIF(true);
+ if (!(service("PileUpMergeSvc", m_mergeSvc, CREATEIF)).isSuccess() || 0 == m_mergeSvc) {
+ ATH_MSG_ERROR ( "Could not find PileUpMergeSvc" );
+ return StatusCode::FAILURE;
+ }
+
+ // StripsResponse Creation
+ m_StripsResponse = new StripsResponse();
+ m_StripsResponse->set_qThreshold(m_qThreshold);
+ m_StripsResponse->set_diffusSigma(m_diffusSigma);
+ m_StripsResponse->set_LogitundinalDiffusSigma(m_LogitundinalDiffusSigma);
+ m_StripsResponse->set_driftGap(m_driftGap);
+ m_StripsResponse->set_driftVelocity(m_driftVelocity);
+ m_StripsResponse->set_crossTalk1(m_crossTalk1);
+ m_StripsResponse->set_crossTalk2(m_crossTalk2);
+
+ m_ElectronicsResponse = new ElectronicsResponse();
+ m_ElectronicsResponse->set_alpha(m_alpha);
+ m_ElectronicsResponse->set_RC(m_RC);
+ m_ElectronicsResponse->set_electronicsThreshold(m_electronicsThreshold);
+
+ //ROOT Staff for internal validation
+ if (m_saveInternalHistos) m_file = new TFile("MM_fullSimu_plots.root","RECREATE");
+ m_ntuple = new TTree("fullSim","fullSim");
+
+ m_ntuple->Branch("exitcode",&exitcode);
+ m_ntuple->Branch("Station_side",&m_n_Station_side);
+ m_ntuple->Branch("Station_eta",&m_n_Station_eta);
+ m_ntuple->Branch("Station_phi",&m_n_Station_phi);
+ m_ntuple->Branch("Station_multilayer",&m_n_Station_multilayer);
+ m_ntuple->Branch("Station_layer",&m_n_Station_layer);
+
+ m_ntuple->Branch("hitPDGId",&m_n_hitPDGId);
+ m_ntuple->Branch("hitKineticEnergy",&m_n_hitKineticEnergy);
+ m_ntuple->Branch("hitDepositEnergy",&m_n_hitDepositEnergy);
+ m_ntuple->Branch("hitOnSurface_x",&m_n_hitOnSurface_x);
+ m_ntuple->Branch("hitOnSurface_y",&m_n_hitOnSurface_y);
+ m_ntuple->Branch("hitStripID",&m_n_hitStripID);
+ m_ntuple->Branch("hitDistToChannel",&m_n_hitDistToChannel);
+ m_ntuple->Branch("hitIncomingAngle",&m_n_hitIncomingAngle);
+ m_ntuple->Branch("hitIncomingAngleRads",&m_n_hitIncomingAngleRads);
+
+ m_ntuple->Branch("StrRespID",&m_n_StrRespID);
+ m_ntuple->Branch("StrRespCharge",&m_n_StrRespCharge);
+ m_ntuple->Branch("StrRespTime",&m_n_StrRespTime);
+ m_ntuple->Branch("StrRespTrg_ID",&m_n_StrRespTrg_ID);
+ m_ntuple->Branch("StrRespTrg_Time",&m_n_StrRespTrg_Time);
+ m_ntuple->Branch("Strip_Multiplicity_byDiffer",&m_n_strip_multiplicity);
+ m_ntuple->Branch("Strip_Multiplicity_2",&m_n_strip_multiplicity_2);
+ m_ntuple->Branch("tofCorrection",&tofCorrection);
+ m_ntuple->Branch("bunchTime",&bunchTime);
+ m_ntuple->Branch("globalHitTime",&globalHitTime);
+
+ m_AngleDistr = new TH1I("m_AngleDistr", "m_AngleDistr", 360, -180., 180.);
+ m_AbsAngleDistr = new TH1I("m_AbsAngleDistr", "m_AbsAngleDistr", 180, 0., 180.);
+ m_ClusterLength2D = new TH1I("m_ClusterLength2D", "m_ClusterLength2D", 100, 0., 100.);
+ m_ClusterLength = new TH1I("m_ClusterLength", "m_ClusterLength", 100, 0., 100.);
+ m_gasGap = new TH1I("m_gasGap", "m_gasGap", 100, 0., 100.);
+ m_gasGapDir = new TH1I("m_gasGapDir", "m_gasGapDir", 20, -10., 10.);
+
+ // m_thpcMM->reserve(20000);
+
+
+ return status;
+}
+/*******************************************************************************/
+//----------------------------------------------------------------------
+// PrepareEvent method:
+//----------------------------------------------------------------------
+StatusCode MmDigitizationTool::prepareEvent(unsigned int nInputEvents) {
+
+ ATH_MSG_DEBUG("MmDigitizationTool::prepareEvent() called for " << nInputEvents << " input events" );
+ //m_digitContainer->cleanup();
+ m_MMHitCollList.clear();
+ // if (0 == m_thpcMM) { /*m_thpcMM->reserve(20000);*/ m_thpcMM = new TimedHitCollection<GenericMuonSimHit>(); }
+ //m_thpcMM = new TimedHitCollection<GenericMuonSimHit>();
+
+ if(!m_thpcMM) {
+ m_thpcMM = new TimedHitCollection<GenericMuonSimHit>();
+ }else{
+ ATH_MSG_ERROR ( "m_thpcMM is not null" );
+ return StatusCode::FAILURE;
+ }
+
+ //m_MMHitCollList.push_back(NULL);
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::processBunchXing(int bunchXing,
+ PileUpEventInfo::SubEvent::const_iterator bSubEvents,
+ PileUpEventInfo::SubEvent::const_iterator eSubEvents) {
+
+ ATH_MSG_DEBUG ( "MmDigitizationTool::in processBunchXing()" << bunchXing );
+
+ PileUpEventInfo::SubEvent::const_iterator iEvt = bSubEvents;
+
+ //loop on event and sub-events for the current bunch Xing
+ for (; iEvt!=eSubEvents; ++iEvt) {
+ StoreGateSvc& seStore = *iEvt->pSubEvtSG;
+ const EventInfo* pEI(0);
+ if (seStore.retrieve(pEI).isSuccess()) {
+ ATH_MSG_INFO( "SubEvt EventInfo from StoreGate " << seStore.name() << " :"
+ << " bunch crossing : " << bunchXing
+ << " time offset : " << iEvt->time()
+ << " event number : " << pEI->event_ID()->event_number()
+ << " run number : " << pEI->event_ID()->run_number() );
+ }
+
+ PileUpTimeEventIndex thisEventIndex = PileUpTimeEventIndex(static_cast<int>(iEvt->time()),iEvt->index());
+
+ const GenericMuonSimHitCollection* seHitColl = 0;
+ if (!seStore.retrieve(seHitColl,m_inputObjectName).isSuccess()) {
+ ATH_MSG_ERROR ( "SubEvent MicroMegas SimHitCollection not found in StoreGate " << seStore.name() );
+ return StatusCode::FAILURE;
+ }
+ ATH_MSG_VERBOSE ( "MicroMegas SimHitCollection found with " << seHitColl->size() << " hits" );
+
+ const double timeOfBCID(static_cast<double>(iEvt->time()));
+ ATH_MSG_DEBUG ( "timeOfBCID " << timeOfBCID );
+
+ //Copy hit Collection
+ GenericMuonSimHitCollection* MMHitColl = new GenericMuonSimHitCollection(m_inputObjectName);
+
+ GenericMuonSimHitCollection::const_iterator i = seHitColl->begin();
+ GenericMuonSimHitCollection::const_iterator e = seHitColl->end();
+
+ // Read hits from this collection
+ for (; i!=e; ++i){
+ GenericMuonSimHit mmhit(*i);
+ MMHitColl->Insert(mmhit);
+ }
+ m_thpcMM->insert(thisEventIndex, MMHitColl);
+
+ //store these for deletion at the end of mergeEvent
+ m_MMHitCollList.push_back(MMHitColl);
+
+ }// while (iEvt != eSubEvents) {
+ return StatusCode::SUCCESS;
+}
+
+/*******************************************************************************/
+StatusCode MmDigitizationTool::getNextEvent() {
+
+ // Get next event and extract collection of hit collections:
+ // This is applicable to non-PileUp Event...
+
+ ATH_MSG_DEBUG ( "MmDigitizationTool::getNextEvent()" );
+
+ if (!m_mergeSvc) {
+ const bool CREATEIF(true);
+ if (!(service("PileUpMergeSvc", m_mergeSvc, CREATEIF)).isSuccess() ||
+ 0 == m_mergeSvc) {
+ ATH_MSG_ERROR ("Could not find PileUpMergeSvc" );
+ return StatusCode::FAILURE;
+ }
+ }
+
+ // initialize pointer
+ //m_thpcMM = 0;
+
+ // get the container(s)
+ typedef PileUpMergeSvc::TimedList<GenericMuonSimHitCollection>::type TimedHitCollList;
+
+ //this is a list<info<time_t, DataLink<GenericMuonSimHitCollection> > >
+ TimedHitCollList hitCollList;
+
+ if (!(m_mergeSvc->retrieveSubEvtsData(m_inputObjectName, hitCollList).isSuccess()) ) {
+ ATH_MSG_ERROR ( "Could not fill TimedHitCollList" );
+ return StatusCode::FAILURE;
+ }
+ if (hitCollList.size()==0) {
+ ATH_MSG_ERROR ( "TimedHitCollList has size 0" );
+ return StatusCode::FAILURE;
+ }
+ else {
+ ATH_MSG_DEBUG ( hitCollList.size() << " MicroMegas SimHitCollections with key " << m_inputObjectName << " found" );
+ }
+
+ // create a new hits collection - Define Hit Collection
+ // m_thpcMM = new TimedHitCollection<GenericMuonSimHit>() ;
+ if(!m_thpcMM) {
+ m_thpcMM = new TimedHitCollection<GenericMuonSimHit>();
+ }else{
+ ATH_MSG_ERROR ( "m_thpcMM is not null" );
+ return StatusCode::FAILURE;
+ }
+
+ //now merge all collections into one
+ TimedHitCollList::iterator iColl(hitCollList.begin());
+ TimedHitCollList::iterator endColl(hitCollList.end());
+
+ // loop on the hit collections
+ while(iColl != endColl) {
+ const GenericMuonSimHitCollection* tmpColl(iColl->second);
+ m_thpcMM->insert(iColl->first, tmpColl);
+ ATH_MSG_DEBUG ( "MMSimHitCollection found with " << tmpColl->size() << " hits" );
+ ++iColl;
+ }
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::mergeEvent() {
+
+ StatusCode status = StatusCode::SUCCESS;
+
+ // ATH_MSG_DEBUG ( "MmDigitizationTool::in mergeEvent()" );
+
+ // Cleanup and record the Digit container in StoreGate
+ status = recordDigitAndSdoContainers();
+ if(!status.isSuccess()) {
+ ATH_MSG_FATAL("MMDigitizationTool::recordDigitAndSdoContainers failed.");
+ return StatusCode::FAILURE;
+ }
+
+ status = doDigitization();
+ if (status.isFailure()) { ATH_MSG_ERROR ( "doDigitization Failed" ); return StatusCode::FAILURE; }
+
+ // reset the pointer (delete null pointer should be safe)
+ if (m_thpcMM){
+ delete m_thpcMM;
+ m_thpcMM = 0;
+ }
+
+ // remove cloned one in processBunchXing......
+ std::list<GenericMuonSimHitCollection*>::iterator MMHitColl = m_MMHitCollList.begin();
+ std::list<GenericMuonSimHitCollection*>::iterator MMHitCollEnd = m_MMHitCollList.end();
+ while(MMHitColl!=MMHitCollEnd) {
+ delete (*MMHitColl);
+ ++MMHitColl;
+ }
+ m_MMHitCollList.clear();
+ return status;
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::digitize() {
+ return this->processAllSubEvents();
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::processAllSubEvents() {
+
+ ATH_MSG_DEBUG ("MmDigitizationTool::processAllSubEvents()");
+
+ StatusCode status = recordDigitAndSdoContainers();
+ //if(!status.isSuccess()) {
+ // ATH_MSG_FATAL("MMDigitizationTool::recordDigitAndSdoContainers failed.");
+ // return StatusCode::FAILURE; // there are no hits in this event
+ //}
+
+ //ATH_MSG_DEBUG ( "MicroMegas SDO collection (MmSDOCollection) recorded in StoreGate." );
+
+ //merging of the hit collection in getNextEvent method
+
+ if (0 == m_thpcMM ) {
+ status = getNextEvent();
+ if(status.isFailure()) {
+ ATH_MSG_FATAL( "There are no MicroMegas hits in this event" );
+ return StatusCode::FAILURE;
+ }
+ }
+
+ status = doDigitization();
+ if(!status.isSuccess()) {
+ ATH_MSG_ERROR( "MmDigitizationTool :: doDigitization() Failed" );
+ return StatusCode::FAILURE;
+ }
+
+ // reset the pointer (delete null pointer should be safe)
+ if (m_thpcMM){
+ delete m_thpcMM;
+ m_thpcMM = 0;
+ }
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::finalize() {
+ //ATH_MSG_DEBUG("MMDigitizationTool::finalize() ---- m_digitContainer->digit_size() = "<<m_digitContainer->digit_size() );
+ m_digitContainer->release();
+
+ if (m_saveInternalHistos) {
+ m_ntuple->Write();
+ m_AngleDistr->Write();
+ m_AbsAngleDistr->Write();
+ m_ClusterLength2D->Write();
+ m_ClusterLength->Write();
+ m_gasGap->Write();
+ m_gasGapDir->Write();
+ m_file->Close();
+ }
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+StatusCode MmDigitizationTool::recordDigitAndSdoContainers() {
+
+ // cleanup digit container
+ m_digitContainer->cleanup();
+
+ // record the digit container in StoreGate
+ m_activeStore->setStore(&*m_sgSvc);
+ StatusCode status = m_sgSvc->record(m_digitContainer, m_outputObjectName);
+ if(status.isFailure()) {
+ ATH_MSG_FATAL("Unable to record Micromegas digit container in StoreGate");
+ return status;
+ }else { ATH_MSG_DEBUG("MmDigitContainer recorded in StoreGate.");}
+
+ // create and record the SDO container in StoreGate
+ m_sdoContainer = new MuonSimDataCollection();
+ status = m_sgSvc->record(m_sdoContainer, m_outputSDOName);
+ if(status.isFailure()) {
+ ATH_MSG_FATAL("Unable to record MM SDO collection in StoreGate");
+ return status;
+ } else { ATH_MSG_DEBUG("MMSDOCollection recorded in StoreGate."); }
+
+ return status;
+}
+
+/*******************************************************************************/
+StatusCode MmDigitizationTool::doDigitization() {
+
+
+ GenericMuonSimHitCollection* inputSimHitColl=NULL;
+
+ if (m_validationSetup)
+ {
+ inputSimHitColl = new GenericMuonSimHitCollection("MicromegasSensitiveDetector");
+ StatusCode status = m_sgSvc->record(inputSimHitColl,"InputMicroMegasHits");
+ if (status.isFailure()) {
+ ATH_MSG_ERROR ( "Unable to record Input MicromegasSensitiveDetector HIT collection in StoreGate" );
+ return status;
+ }
+ }
+
+ if( m_maskMultiplet == 3 ) {
+
+ return StatusCode::SUCCESS;
+ }
+
+ // Perform null check on m_thpcCSC
+ // if(!m_thpcMM) {
+ // ATH_MSG_ERROR ( "m_thpcMM is null" );
+ // return StatusCode::FAILURE;
+ // }
+
+ //iterate over hits and fill id-keyed drift time map
+ TimedHitCollection< GenericMuonSimHit >::const_iterator i, e;
+ const GenericMuonSimHit* previousHit = 0;
+
+ // ATH_MSG_DEBUG("create PRD container of size " << m_idHelper->detectorElement_hash_max());
+ std::map<Identifier,int> hitsPerChannel;
+ int nhits = 0;
+
+ // nextDetectorElement-->sets an iterator range with the hits of current detector element , returns a bool when done
+ while( m_thpcMM->nextDetectorElement(i, e) ) {
+
+ // Loop over the hits:
+ while (i != e) {
+
+ TimedHitPtr<GenericMuonSimHit> phit = *i++;
+ const GenericMuonSimHit& hit(*phit);
+
+ // 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_DEBUG("second hit from a muon: prev " << previousHit->localPosition() << " current " << hit.localPrePosition() << " diff " << diff );
+ if( diff.mag() < m_DiffMagSecondMuonHit ) continue;
+ }
+ m_n_hitPDGId = hit.particleEncoding();
+ m_n_hitDepositEnergy = hit.depositEnergy();
+ m_n_hitKineticEnergy = hit.kineticEnergy();
+
+ 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;
+
+ // ATH_MSG_DEBUG ( "MMDigitTime " << timeOfBCID );
+
+ const float timeWindowStrip = 120.; //driftvelocity gap;
+ if (MMDigitTime < -timeWindowStrip || MMDigitTime > timeWindowStrip) {
+ exitcode = 4; m_ntuple->Fill();
+ continue;
+ }
+
+ // SimHits without energy loss are not recorded.
+ ///if(hit.depositEnergy()==0. && abs(hit.particleEncoding())!=13) continue;
+
+ m_n_Station_side=-111; m_n_Station_eta=-111; m_n_Station_phi=-111; m_n_Station_multilayer=-111; m_n_Station_layer=-111; m_n_hitStripID=-111; m_n_StrRespTrg_ID=-111;
+ m_n_hitOnSurface_x=-999.; m_n_hitOnSurface_y=-999.; m_n_hitDistToChannel=-999.; m_n_hitIncomingAngle=-999.; m_n_StrRespTrg_Time=-999.;
+ m_n_strip_multiplicity = -999.;
+ exitcode = 0;
+ m_n_StrRespID.clear();
+ m_n_StrRespCharge.clear(); m_n_StrRespTime.clear();
+
+ const int idHit = hit.GenericId();
+ // the G4 time or TOF from IP
+ double G4Time(hit.globalTime());
+ // see what are the members of GenericMuonSimHit
+ const Amg::Vector3D globPos = hit.globalPosition();
+
+ // convert sim id helper to offline id
+ muonHelper = MicromegasHitIdHelper::GetHelper();
+ MM_SimIdToOfflineId simToOffline(*m_idHelper);
+
+ //get the hit Identifier and info
+ int simId=hit.GenericId();
+ Identifier layid = simToOffline.convert(simId);
+
+ // Read the information about the Micro Megas hit
+ ATH_MSG_DEBUG ( "> idHit " << idHit << " Hit bunch time " << bunchTime << " tot " << globalHitTime << " tof/G4 time " << hit.globalTime() << " globalPosition " << globPos
+ << "hit: r " << hit.globalPosition().perp() << " z " << hit.globalPosition().z() << " mclink " << hit.particleLink() << " station eta " << m_idHelper->stationEta(layid) << " station phi " << m_idHelper->stationPhi(layid) << " multiplet " << m_idHelper->multilayer(layid) );
+
+ if (m_validationSetup){
+ GenericMuonSimHit* copyHit = new GenericMuonSimHit(idHit, hit.globalpreTime(), hit.globalTime(), hit.globalPosition(), hit.localPosition(), hit.globalPrePosition(), hit.localPrePosition(), hit.particleEncoding(), hit.kineticEnergy(), hit.globalDirection(), hit.depositEnergy(), hit.StepLength(), hit.trackNumber() );
+ ATH_MSG_INFO("Validation: globalHitTime, G4Time, BCtime = "<<globalHitTime<<" "<<G4Time<<" "<<bunchTime << "\n: " << copyHit->print() );
+ inputSimHitColl->Insert(*copyHit);
+ }
+ // Important checks for hits (global time, position along strip, charge, masked chambers etc..) DO NOT SET THIS CHECK TO FALSE IF YOU DON'T KNOW WHAT YOU'RE DOING !
+ if(m_checkMMSimHits) { if(checkMMSimHit(hit) == false) {exitcode = 8; m_ntuple->Fill();continue;} }
+
+
+ // sanity checks
+ if( !m_idHelper->is_mm(layid) ){
+ ATH_MSG_WARNING("MM id is not a mm id! " << m_idHelper->stationNameString(m_idHelper->stationName(layid)) );
+ }
+
+ std::string stName = m_idHelper->stationNameString(m_idHelper->stationName(layid));
+ int isSmall = stName[2] == 'S';
+
+ 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) );
+ exitcode = 9; m_ntuple->Fill();
+ }
+
+ if( m_idHelper->stationPhi(layid) == 0 ){
+ ATH_MSG_WARNING("unexpected phi range " << m_idHelper->stationPhi(layid) );
+ exitcode = 9; m_ntuple->Fill();
+ continue;
+ }
+
+
+ // remove hits in masked multiplet
+ if( m_maskMultiplet == m_idHelper->multilayer(layid) ) continue;
+
+
+ // get readout element
+ const MuonGM::MMReadoutElement* detEl = m_MuonGeoMgr->getMMReadoutElement(layid);
+ if( !detEl ){
+ ATH_MSG_WARNING( "Failed to retrieve detector element for: isSmall " << isSmall << " eta " << m_idHelper->stationEta(layid) << " phi " << m_idHelper->stationPhi(layid) << " ml " << m_idHelper->multilayer(layid) );
+ exitcode = 10; m_ntuple->Fill();
+ continue;
+ }
+
+ // surface
+ const Trk::PlaneSurface& surf = detEl->surface(layid);
+
+ // calculate the inclination angle
+ //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 ;
+ // inAngle_XZ = MM_READOUT [muonHelper->GetLayer(simId)-1]*inAngle_XZ ;
+ ATH_MSG_DEBUG(" At eta " << m_idHelper->stationEta(layid) << " phi " << m_idHelper->stationPhi(layid) << "\n IncomingAngle: " << locDire.angleXZ() / CLHEP::degree << "\n inAngle_XZ, " << inAngle_XZ << " , " << inAngle_XZ * CLHEP::degree << " .. " << CLHEP::degree );
+
+ // compute hit position within the detector element/surfaces
+ // Amg::Transform3D globalToLocal = detEl->absTransform().inverse();
+ Amg::Vector3D hpos(hit.globalPosition().x(),hit.globalPosition().y(),hit.globalPosition().z());
+ // Amg::Vector3D lpos = globalToLocal*hpos;
+ // compute the hit position on the readout plane (same as in MuonFastDigitization)
+ Amg::Vector3D slpos = surf.transform().inverse()*hpos;
+ Amg::Vector2D posOnSurfUnProjected(slpos.x(),slpos.y());
+
+ // double gasGapThickness = detEl->getDesign(layid)->gasGapThickness();
+
+ Amg::Vector3D locdir(0., 0., 0.);
+ if (MM_READOUT[m_idHelper->gasGap(layid)-1]==1) locdir = surf.transform().inverse().linear()*Amg::Vector3D(hit.globalDirection().x(), hit.globalDirection().y(), hit.globalDirection().z());
+ else locdir = surf.transform().inverse().linear()*Amg::Vector3D(hit.globalDirection().x(), hit.globalDirection().y(), -hit.globalDirection().z());
+
+ double scale, scaletop;
+ double gasgap = 5.;
+
+ scale = -slpos.z()/locdir.z();
+ scaletop = (gasgap+slpos.z())/locdir.z();
+ // scaletop = (fabs(gasGapThickness) + slpos.z())/locdir.z();
+
+ Amg::Vector3D hitOnSurface = slpos + scale*locdir;
+ Amg::Vector3D hitOnTopSurface = slpos + scaletop*locdir;
+ Amg::Vector2D posOnSurf (hitOnSurface.x(), hitOnSurface.y());
+ Amg::Vector2D posOnTopSurf (hitOnTopSurface.x(),hitOnTopSurface.y());
+
+ // account for time offset
+ double shiftTimeOffset = MMDigitTime* m_driftVelocity;
+ Amg::Vector3D hitAfterTimeShift(hitOnSurface.x(),hitOnSurface.y(),shiftTimeOffset);
+ Amg::Vector3D hitAfterTimeShiftOnSurface = hitAfterTimeShift - (shiftTimeOffset/locdir.z())*locdir;
+
+ ATH_MSG_DEBUG("slpos.z " << slpos.z() << ", locdir " << locdir.z() << ", scale " << scale << ", hitOnSurface.z " << hitOnSurface.z() << ", hitOnTopSurface.z " << hitOnTopSurface.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 );
+
+ m_n_Station_side = muonHelper->GetSide(simId);
+ m_n_Station_eta = muonHelper->GetZSector(simId);
+ m_n_Station_phi = muonHelper->GetPhiSector(simId);
+ m_n_Station_multilayer = muonHelper->GetMultiLayer(simId);
+ m_n_Station_layer = muonHelper->GetLayer(simId);
+
+ 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;
+ }
+
+ m_gasGap->Fill(m_n_Station_layer);
+ m_gasGapDir->Fill(MM_READOUT [m_n_Station_layer-1]);
+
+ int stripNumber = detEl->stripNumber(posOnSurf,layid);
+ int LastStripNumber = detEl->stripNumber(posOnTopSurf, layid);
+ if( stripNumber == -1 ){
+ ATH_MSG_WARNING("!!! Failed to obtain strip number " << m_idHelper->print_to_string(layid) << "\n\t\t with pos " << posOnSurf
+ << " z " << slpos.z() << " eKin: " << hit.kineticEnergy() << " eDep: " << hit.depositEnergy() << " unprojectedStrip: " << detEl->stripNumber(posOnSurfUnProjected, layid));
+ exitcode = 2; m_ntuple->Fill();
+ continue;
+ }
+ m_n_strip_multiplicity = LastStripNumber - stripNumber;
+
+
+ // re-definition of ID
+ Identifier parentID = m_idHelper->parentID(layid);
+ Identifier DigitId = m_idHelper->channelID(parentID, m_idHelper->multilayer(layid), m_idHelper->gasGap(layid), stripNumber);
+ // ATH_MSG_DEBUG(" re-definition of ID; " << m_idHelper->print_to_string(id) );
+
+ int& counts = hitsPerChannel[DigitId];
+ ++counts;
+ if( counts > 1 ) {exitcode = 11; m_ntuple->Fill(); continue;}
+ ++nhits;
+
+
+ IdentifierHash hash;
+ // contain (name, eta, phi, multiPlet)
+ m_idHelper->get_detectorElement_hash(layid, hash);
+ // ATH_MSG_DEBUG(" looking up collection using hash " << (int)hash << " " << m_idHelper->print_to_string(layid) );
+
+ const MuonGM::MuonChannelDesign* mmChannelDes = detEl->getDesign(DigitId);
+ double distToChannel_withStripID = mmChannelDes->distanceToChannel(posOnSurf, stripNumber);
+ double distToChannel = mmChannelDes->distanceToChannel(posOnSurf);
+ ATH_MSG_DEBUG(" looking up collection using hash " << (int)hash << " " << m_idHelper->print_to_string(layid) << " DigitId: " << m_idHelper->print_to_string(DigitId));
+
+ if ( fabs(distToChannel_withStripID - distToChannel) > mmChannelDes->channelWidth(posOnSurf)) {
+ ATH_MSG_WARNING( "Found: distToChannel_withStripID: " << distToChannel_withStripID << " != distToChannel: " << distToChannel );
+ exitcode = 12; m_ntuple->Fill();
+ continue;
+ }
+
+ // ATH_MSG_DEBUG( "Found: distToChannel_withStripID: " << distToChannel_withStripID << " ? distToChannel: " << distToChannel );
+
+ //starting the digitization
+ const Identifier elemId = m_idHelper -> elementID(DigitId); //
+ // ATH_MSG_DEBUG( "executeDigi() - element identifier is: " << m_idHelper->show_to_string(elemId) );
+
+ //store local hit position + sign
+ // ATH_MSG_DEBUG( " MmDigitToolInput create... " );
+ const MmDigitToolInput StripdigitInput(stripNumber, distToChannel, inAngle_XZ, 0., detEl->numberOfStrips(layid));
+ m_AngleDistr->Fill(inAngle_XZ);
+ m_AbsAngleDistr->Fill(fabs(inAngle_XZ));
+
+
+ // digitize input for strip response
+ m_StripsResponse->set_stripWidth(mmChannelDes->channelWidth(posOnSurf));
+ MmDigitToolOutput StripdigitOutput( m_StripsResponse->GetResponceFrom(StripdigitInput) );
+ if(!StripdigitOutput.isValid()) continue;
+
+ // digitize input for strip response
+ m_StripsResponse->set_stripWidth(mmChannelDes->channelWidth(posOnSurf));
+
+ // simulate strip response, check if strip fired
+ const std::vector<int> & StripsResponse_stripPos = StripdigitOutput.stripPos();
+ const std::vector<float> & StripsResponse_stripTime = StripdigitOutput.stripTime();
+ const std::vector<float> & StripsResponse_stripCharge = StripdigitOutput.stripCharge();
+
+ //ATH_MSG_DEBUG("posOnSurf " << posOnSurf.x() << " " << posOnSurf.y() );
+ //ATH_MSG_DEBUG( " Incoming at: " << stripNumber << " channelWidth: " << mmChannelDes->channelWidth(posOnSurf) << " ? " << m_StripsResponse->get_stripWidth());
+ //float maxChargeFound(0.);
+ float firstTimeFound(100000.);
+ float lastTimeFound(0.);
+ int posForTrigger(0), posForLast(0);
+ for (int i=0; i<(int)StripsResponse_stripPos.size(); i++) {
+ if ( (firstTimeFound > StripsResponse_stripTime[i]) && (StripsResponse_stripCharge[i] > m_qThresholdForTrigger)) {firstTimeFound = StripsResponse_stripTime[i]; posForTrigger = i;}
+ if ( (lastTimeFound < StripsResponse_stripTime[i]) && (StripsResponse_stripCharge[i] > m_qThresholdForTrigger)) {lastTimeFound = StripsResponse_stripTime[i]; posForLast = i;}
+ Amg::Vector2D tmpPos (posOnSurf.x()+ (StripsResponse_stripPos[i] - stripNumber )*0.425 , posOnSurf.y());
+ //ATH_MSG_DEBUG( " For: " << i<< " we have :" << StripsResponse_stripPos[i] << " ? "<< mmChannelDes->channelNumber(tmpPos) << " with charge: " << StripsResponse_stripCharge[i] << " time: " << StripsResponse_stripTime[i] << " timeForTrig: " << firstTimeFound);
+ }
+ m_n_strip_multiplicity_2 = MM_READOUT [muonHelper->GetLayer(simId)-1]*StripsResponse_stripPos[StripsResponse_stripPos.size()-1] - (-MM_READOUT [muonHelper->GetLayer(simId)-1])*StripsResponse_stripPos[0];
+
+ //ATH_MSG_DEBUG("Final: MM_READOUT ["<< m_n_Station_layer-1 <<"] "<< MM_READOUT [m_n_Station_layer-1] << " ?? " << m_n_strip_multiplicity << " from " << stripNumber << " to " << LastStripNumber << " while: " << StripsResponse_stripPos[0] << " to " << StripsResponse_stripPos[StripsResponse_stripPos.size()-1]);
+
+ StripdigitOutput.set_StripForTrigger(StripsResponse_stripPos[posForTrigger]);
+ StripdigitOutput.set_StripTimeForTrigger(StripsResponse_stripTime[posForTrigger]);
+ m_ClusterLength->Fill(StripsResponse_stripPos[StripsResponse_stripPos.size()-1] - StripsResponse_stripPos[0] );
+ m_ClusterLength2D->Fill(StripsResponse_stripPos[posForLast] - StripsResponse_stripPos[posForTrigger]);
+
+ m_n_hitStripID=stripNumber;
+ m_n_StrRespTrg_ID=StripsResponse_stripPos[posForTrigger];
+ m_n_hitDistToChannel=distToChannel;
+ m_n_hitIncomingAngle=inAngle_XZ;
+ m_n_hitIncomingAngleRads = inAngle_XZ * CLHEP::degree;
+ m_n_StrRespTrg_Time=StripsResponse_stripTime[posForTrigger];
+ m_n_hitOnSurface_x=posOnSurf.x();
+ m_n_hitOnSurface_y = posOnSurf.y();
+
+ m_n_StrRespID = StripsResponse_stripPos;
+ m_n_StrRespCharge = StripsResponse_stripCharge;
+ m_n_StrRespTime = StripsResponse_stripTime;
+
+ m_ntuple->Fill();
+
+ //ATH_MSG_DEBUG( " Cluster length: " << StripsResponse_stripPos[StripsResponse_stripPos.size()-1] - StripsResponse_stripPos[0] );
+
+ const MmElectronicsToolInput ElectronicDigitInput( StripdigitOutput.stripPos(), StripdigitOutput.stripCharge(), StripdigitOutput.stripTime() );
+ MmDigitToolOutput ElectronicOutput( m_ElectronicsResponse->GetResponceFrom(ElectronicDigitInput) );
+
+ if(!ElectronicOutput.isValid()) ATH_MSG_WARNING ( "MmDigitizationTool::doDigitization() -- there is no electronics response even though there is a strip response." );
+
+ /* unused variables to prevent compilation warnings */
+ /*
+
+ const std::vector<int> & chipResponsePos = ElectronicOutput.stripPos();
+ const std::vector<float> & chipResponseTime = ElectronicOutput.stripTime();
+ const std::vector<float> & chipResponseCharge = ElectronicOutput.stripCharge();
+ */
+
+ // The collections should use the detector element hashes not the module hashes to be consistent with the PRD granularity.
+ IdentifierHash detIdhash ;
+ // set RE hash id
+ int gethash_code = m_idHelper->get_detectorElement_hash(elemId, detIdhash);
+ if (gethash_code != 0) {
+ ATH_MSG_ERROR ( "MmDigitizationTool -- collection hash Id NOT computed for id = " << m_idHelper->show_to_string(elemId) );
+ // continue;
+ }
+
+ MmDigit* newDigit = new MmDigit(DigitId, StripdigitOutput.stripTime(), StripdigitOutput.stripPos(), StripdigitOutput.stripCharge(),
+ ElectronicOutput.stripTime(),ElectronicOutput.stripPos(),ElectronicOutput.stripCharge(),
+ StripdigitOutput.stripForTrigger(), StripdigitOutput.stripTimeForTrigger() );
+
+ //ATH_MSG_DEBUG ( "strip and chip response digit created" );
+
+ MmDigitCollection* digitCollection = 0;
+
+ // put new collection in storegate
+ // Get the messaging service, print where you are
+ m_activeStore->setStore( &*m_sgSvc );
+ MmDigitContainer::const_iterator it_coll = m_digitContainer->indexFind(detIdhash);
+ if (m_digitContainer->end() == it_coll) {
+ digitCollection = new MmDigitCollection(elemId, detIdhash);
+ digitCollection->push_back(newDigit);
+ m_activeStore->setStore( &*m_sgSvc );
+ StatusCode status = m_digitContainer->addCollection(digitCollection, detIdhash);
+ // ATH_MSG_DEBUG ( "MM digitCollection status: " << status );
+ if (status.isFailure()) {
+ ATH_MSG_ERROR ( "Couldn't record MicroMegas DigitCollection with key=" << detIdhash << " in StoreGate!" );
+ // delete digitCollection;
+ // digitCollection = 0 ;
+ // return StatusCode::RECOVERABLE; // consistent with ERROR message above.
+ }
+ //else ATH_MSG_DEBUG ( "New MicroMegas DigitCollection with key=" << detIdhash << " recorded in StoreGate." );
+ }
+ else {
+ digitCollection = const_cast<MmDigitCollection*>( it_coll->cptr() );
+ digitCollection->push_back(newDigit);
+ }
+
+ // fill the SDO collection in StoreGate
+ // create here deposit for MuonSimData, link and tof
+ ///MuonSimData::Deposit deposit(HepMcParticleLink(phit->trackNumber(), phit.eventId()), MuonMCData(G4Time, 0));
+ MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(hitOnSurface.x(), hitOnSurface.y()));
+ std::vector<MuonSimData::Deposit> deposits;
+ deposits.push_back(deposit);
+ m_sdoContainer->insert(std::make_pair(DigitId, MuonSimData(deposits, inAngle_XZ * CLHEP::degree*1000)));
+
+ // ATH_MSG_DEBUG( " MmDigitToolInput outsideWindow: " << outsideWindow(m_bunchTime) );
+ // it should be decided by bunchtime. -400 to +200
+ // if (outsideWindow(m_bunchTime)) continue;
+
+ //digitCollection = getDigitCollection(elemId);
+ //ATH_MSG_DEBUG ( "Finished MmDigitizationTool::MM digit and SDO stored in the DigitCollection" );
+ previousHit = &hit;
+ }//while(i != e)
+ }//while(m_thpcMM->nextDetectorElement(i, e))
+ // reset the pointer if it is not null
+
+ ATH_MSG_DEBUG ( "MmDigitization Done!" );
+
+ if (m_thpcMM){
+ delete m_thpcMM;
+ m_thpcMM = 0;
+ }
+
+ return StatusCode::SUCCESS;
+}
+/*******************************************************************************/
+bool MmDigitizationTool::checkMMSimHit( const GenericMuonSimHit& /*hit*/ ) const {
+ // ATH_MSG_DEBUG ("MmDigitizationTool::checkMMSimHit()");
+ /*
+ //get the hit Identifier and info
+ const int simId = hit.GenericId();
+ std::string stationName = muonHelper->GetStationName(simId);
+ int stationEta = muonHelper->GetZSector(simId);
+ int stationPhi = muonHelper->GetPhiSector(simId);
+ int multilayer = muonHelper->GetMultiLayer(simId); // mmMultiplet in MmIdHelpers
+ int layer = muonHelper->GetLayer(simId); // mmGasGap in MmIdHelpers
+ int side = muonHelper->GetSide(simId);
+
+ // ATH_MSG_VERBOSE("Micromegas hit: r " << hit.globalPosition().perp() << " z " << hit.globalPosition().z() << " mclink " << hit.particleLink()
+ // << " stName " << stationName << " eta " << stationEta << " phi " << stationPhi << " ml " << multilayer << " layer " << layer << " side " << side );
+
+ // the G4 time or TOF from IP
+ double MM_Hit_globalTime = hit.globalTime() ;
+ double MM_Hit_globalPreTime = hit.globalpreTime() ;
+
+ double timeOfFlight = 0.5* ( MM_Hit_globalPreTime - MM_Hit_globalTime );
+
+ // extract the distance to the origin of the module to Time of flight
+ // Use the coordinate of the center of the module to calculate the time of flight
+ Amg::Vector3D distance = 0.5* ( hit.globalPrePosition() + hit.globalPosition() );
+ double m_time = timeOfFlight/CLHEP::ns - distance.mag()/CLHEP::c_light/CLHEP::ns;
+
+ ATH_MSG_VERBOSE( "Micromegas hit "
+ << " MM_Hit_globalPreTime [" << MM_Hit_globalPreTime << "]"
+ << " MM_Hit_globalTime [" << MM_Hit_globalTime << "]"
+ << " timeOfFlight [" << timeOfFlight << "]"
+ << " point [" << distance << "]"
+ << " m_time [" << m_time << "]" );
+ */
+
+ return true;
+}
+/*******************************************************************************/
+int MmDigitizationTool::digitizeTime(double time) const {
+
+ int tdcCount;
+ double tmpCount = time/m_ns2TDC;
+ double rand = CLHEP::RandGaussZiggurat::shoot(m_rndmEngine, tmpCount, m_resTDC);
+ tdcCount = static_cast<long>(rand);
+
+ if (tdcCount < 0 || tdcCount > 4096){
+ //ATH_MSG_DEBUG( " Count outside TDC window: " << tdcCount );
+ }
+
+ return tdcCount;
+}
+/*******************************************************************************/
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/StripsResponse.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/StripsResponse.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0839204f7d1e23d92f5a3ea91a97e8b627a92fc0
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/StripsResponse.cxx
@@ -0,0 +1,319 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * StripsResponse.cxx
+ * MC for micromegas athena integration
+ *
+ **/
+
+/// PROJECTS
+#include "MM_Digitization/StripsResponse.h"
+#include "GaudiKernel/MsgStream.h"
+/*******************************************************************************/
+StripsResponse::StripsResponse():
+ qThreshold(0.001),
+ diffusSigma(0.360/10.),
+ LogitundinalDiffusSigma(0.190/10.),
+ pitch(0.500),
+ stripwidth(pitch - 0.0),
+ polya(0),
+ conv_gaus(0),
+ driftGap(5.128),
+ driftVelocity(0.047),
+ gRandom_loc(0),
+ randomNum(0)
+{
+ clearValues();
+ initializationFrom();
+}
+/*******************************************************************************/
+void StripsResponse::initHistos()
+{
+ // For validation purposes only...
+}
+/*******************************************************************************/
+void StripsResponse::writeHistos()
+{}
+/*******************************************************************************/
+void StripsResponse::initFunctions()
+{
+ polya = new TF1("polya","(1./[1])*(TMath::Power([0]+1,[0]+1)/TMath::Gamma([0]+1))*TMath::Power(x,[0])*TMath::Exp(-([0]+1)*x)", 0., 4.);
+ conv_gaus = new TF1("conv_gaus","1.*TMath::Exp(-TMath::Power(x,2.)/(2.*[0]*[0])) + 0.001*TMath::Exp(-TMath::Power(x,2)/(2.*[1]*[1]))", -1., 1.);
+ Athena::MsgStreamMember log("tripsResponse_uM::initFunctions");
+ log << MSG::DEBUG << "tripsResponse_uM::initFunctions DONE" << endreq;
+ // if(msgLevel(MSG::DEBUG)) msg(MSG::DEBUG) << "\t \t StripsResponse::initFunctions DONE " << endreq;
+}
+/*******************************************************************************/
+void StripsResponse::clearValues()
+{
+ /// Delete Gen:
+ if (conv_gaus !=0) delete conv_gaus;
+ if (polya !=0) delete polya;
+ if (gRandom_loc !=0) delete gRandom_loc;
+}
+/*******************************************************************************/
+void StripsResponse::initializationFrom()
+{
+ crossTalk1=0.0;
+ crossTalk2=0.0;
+ Lorentz_Angle = 0.0;
+
+ initHistos ();
+ initFunctions();
+ randomNum = new TRandom3(1);
+ randomNum->SetSeed(1);
+
+ Athena::MsgStreamMember log("StripsResponse::initializationFrom");
+ log << MSG::DEBUG << "StripsResponse::initializationFrom set values" << endreq;
+ polya->SetParameter(0, 2.3);
+ polya->SetParameter(1, 1.);
+ log << MSG::DEBUG << "StripsResponse::initializationFrom getRandom: " << polya->GetRandom() << endreq;
+}
+/*******************************************************************************/
+void StripsResponse::testRandomNumGens(UInt_t sd)
+{
+ TRandom3 *randomNum_ = new TRandom3(1);
+ TRandom *rsave = gRandom;
+ TRandom *r0 = new TRandom(1);
+ Athena::MsgStreamMember log("StripsResponse::testRandomNumGens");
+ rsave->Clear();
+ r0->Clear();
+ r0->SetSeed(1);
+ gRandom->Clear();
+ gRandom->SetSeed(sd);//always the same random numbers !!!
+ rsave->SetSeed(1);//always the same random numbers !!!
+ printf("Seeds a: %d | %d | %d\n", gRandom->GetSeed(), rsave->GetSeed(), r0->GetSeed());
+ float pt_ = gRandom->Uniform();
+ float pt = rsave->Uniform();
+ float pt0 = r0->Uniform();
+ for (int i=0; i< 6; i++) {
+ pt_ = gRandom->Uniform();
+ pt = rsave->Uniform();
+ pt0 = r0->Uniform();
+
+ log << MSG::DEBUG <<"\tpt <> " << randomNum_->Uniform() << " " << pt << " " << pt_ << " " << pt0 << " | " << gRandom->GetSeed() << endreq;
+ log << MSG::DEBUG << gRandom->GetSeed() << " " << randomNum_->GetSeed() << endreq;
+ }
+ log << MSG::DEBUG << gRandom->GetSeed() << " " << randomNum_->GetSeed() << endreq;
+}
+/*******************************************************************************/
+MmDigitToolOutput StripsResponse::GetResponceFrom(const MmDigitToolInput & digiInput)
+{
+ Athena::MsgStreamMember log("StripsResponse::GetResponceFrom");
+ log << MSG::DEBUG << "\t \t StripsResponse::GetResponceFrom start " << endreq;
+
+ /// Reduce internal calls of functions
+ calculateResponceFrom(digiInput.positionWithinStrip() , digiInput.stripIDLocal() , digiInput.incomingAngle() , digiInput.stripMaxID() );
+
+ log << MSG::DEBUG << "\t \t StripsResponse::GetResponceFrom MmDigitToolOutput create " << endreq;
+ /// ToDo: include loop for calculating Trigger study vars
+
+ // MmDigitToolOutput(bool hitWasEff, std::vector <int> strpos, std::vector<float> time, std::vector<int> charge, int strTrig, float strTimeTrig ):
+ MmDigitToolOutput tmp(true, finalNumberofStrip, finaltStrip, finalqStrip, 5, 0.3);
+ log << MSG::DEBUG << "StripsResponse::GetResponceFrom -> hitWasEfficient(): " << tmp.hitWasEfficient() << endreq;
+ return tmp;
+}
+/*******************************************************************************/
+void StripsResponse::calculateResponceFrom(const float & hitx, const int & stripID, const float & thetaD, const int & stripMaxID) /// <-------- stripID!!!
+{
+ Athena::MsgStreamMember log("StripsResponse::calculateResponceFrom");
+ log << MSG::DEBUG << "\t \t StripsResponse::calculateResponceFrom start " << endreq;
+ finalNumberofStrip.clear();
+ finalqStrip.clear();
+ finaltStrip.clear();
+ crossTalk1=0.0;
+ crossTalk2=0.0;
+ Lorentz_Angle = 0.0;
+ log << MSG::DEBUG << "\t \t StripsResponse::calculateResponceFrom call whichStrips " << endreq;
+ whichStrips(hitx, stripID, thetaD, stripMaxID);
+}
+//__________________________________________________________________________________________________
+//==================================================================================================
+// calculate the strips that got fired, charge and time. Assume any angle thetaD....
+//==================================================================================================
+void StripsResponse::whichStrips(const float & hitx, const int & stripID, const float & thetaD, const int & stripMaxID){
+
+ dimClusters=7000; //dimClusters=total number of collisions
+ MaxPrimaryIons = 300;
+ pt=0.;
+ xx=0.;
+ xxDiffussion=0.;
+ yy=0.;
+ yyDiffussion=0.;
+
+ // NelectronProb has the probability that a primary ionization cluster
+ float NelectronProb[23]={65.6,80.6,87.0,90.5,92.75,94.30,95.35,96.16,96.77,97.26,
+ 97.65,97.95,98.20,98.40,98.56,98.68,98.775,98.850,98.913,98.967,99.016,99.061,99.105};
+
+ lmean = 10./26.; // 24 e/8 mm or 23 e/10mm (as it has been used here) according to Sauli's paper, for Argon
+
+ float theta = thetaD*TMath::Pi()/180.0, xoff=0.;
+
+
+ //define the polya distribution to get the flactuations of the charges...
+ // Double_t polya_alpha=1., polya_beta=(0.02-1.)/4., temp_polya, polya_theta=0.43;
+ polya_theta=2.3;
+
+ //name, name of the above C++ function, min,max, number of parameters (a,b)
+ //Initialize parameters
+ polya->SetParameter(0, polya_theta);
+ polya->SetParameter(1, 1.);
+
+ // zero charge, time of primary ionization clusters etc
+ l.clear(); l.resize(dimClusters, 0.);
+ clusterelectrons.clear(); clusterelectrons.resize(dimClusters, 0.);
+ qstrip.clear(); qstrip.resize(dimClusters, 0.);
+ time.clear(); time.resize(dimClusters, 0.);
+ cntTimes.clear(); cntTimes.resize(dimClusters, 0.);
+ firstq.clear(); firstq.resize(dimClusters, 0.);
+ ll = 0;
+ totalelectrons = 0;
+ primaryion = 0;
+
+ // Now generate the primary ionization colisions, Poisson distributed
+ while (ll<driftGap/TMath::Cos(theta)){
+ pt = randomNum->Uniform();
+ ll = ll-lmean*log(pt);
+ if (ll>driftGap/TMath::Cos(theta)) break; // if you exceed the driftGap, go out...
+ l[primaryion]=ll;
+ pt = 100.0*randomNum->Uniform();
+ // Use the Nelectron probability to calculate number of electrons per collision
+
+ float limProb1=NelectronProb[22], limProb2;
+
+ for (int k = 0; k<dimClusters; k++){
+
+ if (k<23) {
+ if (pt <= NelectronProb[0]) {
+ clusterelectrons[primaryion] = clusterelectrons[primaryion] + k + 1;
+ totalelectrons = totalelectrons + clusterelectrons[primaryion];
+ break;
+ } else
+ if (pt <= NelectronProb[k] && pt > NelectronProb[k-1]){
+ clusterelectrons[primaryion] = clusterelectrons[primaryion] + k + 1;
+ totalelectrons = totalelectrons + clusterelectrons[primaryion];
+ break;
+ }
+ } else {
+ // more probabilities according to Fischle paper
+ limProb2 = limProb1 + 21.6/((k)*(k));
+ if (pt <= limProb2 && pt > limProb1){
+ clusterelectrons[primaryion] = clusterelectrons[primaryion] + k + 1;
+ totalelectrons = totalelectrons + clusterelectrons[primaryion];
+ break;
+ }
+
+ limProb1 = limProb2;
+ }
+ }// end of loop for total cluster ionization
+
+ primaryion = primaryion + 1;
+ if (primaryion >= MaxPrimaryIons) break; //don't create more than "MaxPrimaryIons" along a track....
+ }//end of clusters loop
+ for (int m=0; m<primaryion; m++){
+ yy = l[m]*cos(theta); // path of the drifted electrons
+ for(int j=1;j<=clusterelectrons[m];j++){
+ //
+ //===================================================
+ // calculate the new charge due to the Polya fluctiation of the avalanche
+ // do this for each electron separately rather than the full cluster (by George Iakovidis)
+ temp_polya = polya->GetRandom();
+ //===================================================
+ //
+ yyDiffussion = randomNum->Gaus(yy, LogitundinalDiffusSigma*yy); // insert longitundinal diffusion of electron
+ if(diffusSigma==0 || LogitundinalDiffusSigma==0) conv_gaus->SetParameters(diffusSigma*yy, 0.0);
+ else conv_gaus->SetParameters(diffusSigma*yy, 1.0);
+
+ float temp_gaus_conv = conv_gaus->GetRandom();
+ xxDiffussion = hitx + l[m]*sin(theta) + temp_gaus_conv;// insert transverese diffusion of electron
+
+ xx = xxDiffussion + TMath::Tan(Lorentz_Angle * TMath::Pi()/180.)*yyDiffussion;
+
+ nstrip = xx/(stripwidth+(pitch-stripwidth));
+ if (xx<0) nstrip = nstrip - 1;
+ nstrip = nstrip + stripID;
+
+ if (xx-hitx>l[m]*sin(theta)) {
+ xoff = xx-(hitx+l[m]*sin(theta));
+ } else {
+ xoff = (hitx+l[m]*sin(theta))-xx;
+ }
+ numberOfElectrons = 1.0*temp_polya;//single electron fluctuation
+ numberofStrip.push_back(nstrip); qStrip.push_back(numberOfElectrons); tStrip.push_back(TMath::Sqrt(yyDiffussion*yyDiffussion+xoff*xoff)/driftVelocity);
+ }
+
+ // for each charge, assigne a crosstalk to the two neighbor strips... crossTalk1 for the strip+-1 and crossTalk2 for the strip+-2
+ // do strip+-1 crosstalk only if it is greater than 0% (by George Iakovidis)
+ if (crossTalk1 > 0.0) {
+ if (nstrip-1>-1) {
+ numberofStrip.push_back(nstrip-1); qStrip.push_back(crossTalk1*numberOfElectrons);
+ tStrip.push_back(TMath::Sqrt(yyDiffussion*yyDiffussion+xoff*xoff)/driftVelocity);
+ }
+ if (nstrip+1>-1) {
+ numberofStrip.push_back(nstrip+1); qStrip.push_back(crossTalk1*numberOfElectrons);
+ tStrip.push_back(TMath::Sqrt(yyDiffussion*yyDiffussion+xoff*xoff)/driftVelocity);
+ }
+ }
+ // do strip+-2 crosstalk only if it is greater than 0% (by George Iakovidis)
+ if (crossTalk2 > 0.0) {
+ if (nstrip-2>-1) {
+ numberofStrip.push_back(nstrip-2); qStrip.push_back(crossTalk2*numberOfElectrons);
+ tStrip.push_back(TMath::Sqrt(yyDiffussion*yyDiffussion+xoff*xoff)/driftVelocity);
+ }
+ if (nstrip+2>-1) {
+ numberofStrip.push_back(nstrip+2); qStrip.push_back(crossTalk2*numberOfElectrons);
+ tStrip.push_back(TMath::Sqrt(yyDiffussion*yyDiffussion+xoff*xoff)/driftVelocity);
+ }
+ }
+ // end of crosstalk assignemnt
+
+ }
+ //calculate the time per strip, in case you cross the qThreshold=2e-
+ nstrip = 0;
+ for (size_t ii = 0; ii<numberofStrip.size(); ii++){
+ nstrip = numberofStrip.at(ii);
+ if (qstrip[nstrip] == 0) firstq[nstrip] = 0;
+ qstrip[nstrip] = qstrip[nstrip] + qStrip.at(ii);
+
+ // find the strip that receives charge above the qThreshold, assign as time the mean time of the arriving e
+ if ((qstrip[nstrip] > 0.0 && firstq[nstrip] == 0 && qstrip[nstrip]<qThreshold) ||
+ (qstrip[nstrip] > 0.0 && firstq[nstrip] == 0 && qstrip[nstrip]>=qThreshold && cntTimes[nstrip]>=0.)) {
+ time[nstrip] = time[nstrip] + tStrip.at(ii);
+ cntTimes[nstrip] = cntTimes[nstrip]+1;
+ }
+ if (qstrip[nstrip]>=qThreshold && firstq[nstrip] == 0) {
+ time[nstrip] = time[nstrip]/cntTimes[nstrip];
+ firstq[nstrip] = 1;
+ cntTimes[nstrip]=0;
+ stripNumber.push_back(nstrip); // record the strip number that got a hit for 1rst time...
+ }
+ }
+
+ int minStripEv=1000000, maxStripEv=-10000000;
+ for (size_t ii = 0; ii<numberofStrip.size(); ii++){
+ if (qStrip.at(ii)<qThreshold) {continue;}
+ if (minStripEv>=numberofStrip.at(ii)) { minStripEv=numberofStrip.at(ii);}
+ if (maxStripEv<=numberofStrip.at(ii)) { maxStripEv=numberofStrip.at(ii);}
+ }
+ for (int ii=minStripEv; ii<=maxStripEv; ii++) {
+ if ( qstrip[ii]==0 || qstrip[ii]<qThreshold) { continue;} // bad strip... ignore it...
+ if ( ii > stripMaxID || ii == 0 ) continue; // do not consider strips outside allowed range, i.e. outside chamber
+
+ finalNumberofStrip.push_back(ii); finalqStrip.push_back(qstrip[ii]); finaltStrip.push_back(time[ii]);
+ }
+
+ //
+ //****************************************************************************
+ // end of debug=10, plot the pulse of the strip... qStrip vs strip number
+ //****************************************************************************
+ //
+ numberofStrip.clear(); qStrip.clear(); tStrip.clear();
+} // end of whichStrips()
+/*******************************************************************************/
+StripsResponse::~StripsResponse()
+{
+ clearValues();
+}
+/*******************************************************************************/
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_entries.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9d856cfddfb44ea7ca6658f1fb4a1b910b71b534
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_entries.cxx
@@ -0,0 +1,16 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "MM_Digitization/MmDigitizationTool.h"
+#include "MM_Digitization/MM_Digitizer.h"
+#include "MM_Digitization/MM_Response_DigitTool.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+DECLARE_ALGORITHM_FACTORY( MM_Digitizer )
+
+DECLARE_TOOL_FACTORY( MM_Response_DigitTool )
+DECLARE_TOOL_FACTORY( MmDigitizationTool )
+
+DECLARE_FACTORY_ENTRIES(MM_Digitization) {
+ DECLARE_ALGORITHM( MM_Digitizer )
+ DECLARE_TOOL( MM_Response_DigitTool )
+ DECLARE_TOOL( MmDigitizationTool )
+}
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_load.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_load.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a4418a12d1ce5e43632e7b0ac7b5847cf626099a
--- /dev/null
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/components/MM_Digitization_load.cxx
@@ -0,0 +1,3 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES(MM_Digitization)