From 5e45f70056e4cb9d0ee6a2935381ab93659ac4a1 Mon Sep 17 00:00:00 2001
From: Alexandre Laurier <alexandre.laurier@cern.ch>
Date: Tue, 28 Apr 2020 21:17:54 +0200
Subject: [PATCH] Add funtions to MM Digitization to add only read strips

---
 .../MuonReadoutGeometry/MMReadoutElement.h    | 20 +++++-
 .../MM_Digitization/MM_DigitToolInput.h       |  5 +-
 .../MM_Digitization/MM_StripResponse.h        |  3 +-
 .../MM_StripsResponseSimulation.h             |  2 +-
 .../src/MM_DigitizationTool.cxx               | 69 +++++++++----------
 .../MM_Digitization/src/MM_StripResponse.cxx  | 42 ++++++-----
 .../src/MM_StripsResponseSimulation.cxx       |  4 +-
 7 files changed, 88 insertions(+), 57 deletions(-)

diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
index 191d84d3eed..5d7f743a66f 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
@@ -59,7 +59,7 @@ namespace MuonGM {
     bool stripPosition(       const Identifier& id, Amg::Vector2D& pos )  const;
     bool stripGlobalPosition( const Identifier& id, Amg::Vector3D& gpos ) const;
       
-      double stripLength( const Identifier& id) const;
+    double stripLength( const Identifier& id) const;
 
     /** number of layers in phi/eta projection */
     int numberOfLayers( bool ) const;
@@ -68,6 +68,10 @@ namespace MuonGM {
     int numberOfStrips( const Identifier& layerId )   const;
     int numberOfStrips( int , bool measuresPhi ) const;
 
+    /** Number of missing bottom and top strips (not read out) */
+    int numberOfMissingTopStrips( const Identifier& layerId )   const;
+    int numberOfMissingBottomStrips( const Identifier& layerId )   const;
+
     /** space point position for a given pair of phi and eta identifiers 
 	The LocalPosition is expressed in the reference frame of the phi surface.
 	If one of the identifiers is outside the valid range, the function will return false */
@@ -254,6 +258,20 @@ namespace MuonGM {
     else return -1;
   }
 
+  inline int MMReadoutElement::numberOfMissingTopStrips( const Identifier& id )   const {
+    const MuonChannelDesign* design = getDesign(id);
+    if( !design ) return -1;
+    int nStrips = design->sAngle == 0 ? design->nMissedTopEta : design->nMissedTopStereo;
+    return nStrips;
+  }
+
+  inline int MMReadoutElement::numberOfMissingBottomStrips( const Identifier& id )   const {
+    const MuonChannelDesign* design = getDesign(id);
+    if( !design ) return -1;
+    int nStrips = design->sAngle == 0 ? design->nMissedBottomEta : design->nMissedBottomStereo;
+    return nStrips;
+  }
+
   inline bool MMReadoutElement::spacePointPosition( const Identifier& phiId, const Identifier& etaId, Amg::Vector2D& pos ) const {
     Amg::Vector2D phiPos;
     Amg::Vector2D etaPos;
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitToolInput.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitToolInput.h
index f44273073e7..ed3bae9741f 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitToolInput.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitToolInput.h
@@ -25,13 +25,14 @@ Class to store input needed for the MM_Digitization tools:
 class MM_DigitToolInput {
  public:
 
- MM_DigitToolInput(int stripIdLocal, double posx, double incomingAngleXZ, double incomingAngleYZ, const Amg::Vector3D& magneticField, int stripMaxId, int gasgap, float eventTime)
+ MM_DigitToolInput(int stripIdLocal, double posx, double incomingAngleXZ, double incomingAngleYZ, const Amg::Vector3D& magneticField, int stripMinId, int stripMaxId, int gasgap, float eventTime)
 
    :  m_stripIDLocal(stripIdLocal),
       m_xpos(posx),
       m_incomingAngleXZ(incomingAngleXZ),
       m_incomingAngleYZ(incomingAngleYZ),
       m_magneticField(magneticField),
+      m_stripMinId(stripMinId),
       m_stripMaxId(stripMaxId),
       m_gasgap(gasgap),
       m_eventTime(eventTime)
@@ -45,6 +46,7 @@ class MM_DigitToolInput {
       double incomingAngleXZ()       const { return m_incomingAngleXZ; }
       double incomingAngleYZ()       const { return m_incomingAngleYZ; }
       const Amg::Vector3D& magneticField()       const { return m_magneticField; } // kT unit, local cordinate
+      int    stripMinID()          const { return m_stripMinId; }
       int    stripMaxID()          const { return m_stripMaxId; }
       int    gasgap()              const { return m_gasgap; }
       Identifier getHitID()        const { return m_hitID; }
@@ -56,6 +58,7 @@ class MM_DigitToolInput {
       double m_incomingAngleXZ;
       double m_incomingAngleYZ;
       Amg::Vector3D m_magneticField;
+      int    m_stripMinId;
       int    m_stripMaxId;
       int    m_gasgap;
       Identifier m_hitID;
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripResponse.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripResponse.h
index a2bddb14e5d..f3b5d1612d6 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripResponse.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripResponse.h
@@ -22,7 +22,7 @@ class MM_StripResponse {
  public:
 
   MM_StripResponse();
-  MM_StripResponse(std::vector<MM_IonizationCluster> IonizationClusters, float timeResolution, float stripPitch, int stripID, int maxstripID);
+  MM_StripResponse(std::vector<MM_IonizationCluster> IonizationClusters, float timeResolution, float stripPitch, int stripID, int minstripID, int maxstripID);
   void timeOrderElectrons();
   void calculateTimeSeries(float thetaD, int gasgap);
   //  void calculateTimeSeries();
@@ -47,6 +47,7 @@ class MM_StripResponse {
   float m_timeResolution;
   float m_stripPitch;
   int m_stripID;
+  int m_minstripID;
   int m_maxstripID;
 
   std::vector<MM_Electron*> m_Electrons;
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
index b65b7c24260..8f8928799f3 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
@@ -86,7 +86,7 @@ public :
   void initHistos ();
   void clearValues ();
   void initFunctions ();
-  void whichStrips(const float & hitx, const int & stripOffest, const float & incidentAngleXZ, const float & incidentAngleYZ, const int & stripMaxID, const MM_DigitToolInput & digiInput);
+  void whichStrips(const float & hitx, const int & stripOffest, const float & incidentAngleXZ, const float & incidentAngleYZ, const int & stripMinID, const int & stripMaxID, const MM_DigitToolInput & digiInput);
 
   inline void setQThreshold (float val) { m_qThreshold = val; };
   inline void setTransverseDiffusionSigma (float val) { m_transverseDiffusionSigma = val; };
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
index 0f934c51901..e4cbb08d767 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
@@ -929,23 +929,23 @@ StatusCode MM_DigitizationTool::doDigitization() {
       Amg::Vector2D tmp (stripLayerPosition.x(), stripLayerPosition.y());
       
       if( stripNumber == -1 ){
-	ATH_MSG_WARNING("!!! Failed to obtain strip number "
-			<< m_idHelper->print_to_string(layerID)
-			<<  "\n\t\t with pos "
-			<< positionOnSurface
-			<< " z "
-			<< stripLayerPosition.z()
-			<< " eKin: "
-			<< hit.kineticEnergy()
-			<< " eDep: "
-			<< hit.depositEnergy()
-			<< " unprojectedStrip: "
-			<< detectorReadoutElement->stripNumber(positionOnSurfaceUnprojected, layerID)
-			);
-	m_exitcode = 2;
-	if(m_writeOutputFile) m_ntuple->Fill();
-	ATH_MSG_DEBUG( "m_exitcode = 2 " );
-	continue;
+	      ATH_MSG_WARNING("!!! Failed to obtain strip number "
+			    << m_idHelper->print_to_string(layerID)
+			    <<  "\n\t\t with pos "
+			    << positionOnSurface
+			    << " z "
+			    << stripLayerPosition.z()
+			    << " eKin: "
+			    << hit.kineticEnergy()
+			    << " eDep: "
+			    << hit.depositEnergy()
+			    << " unprojectedStrip: "
+			    << detectorReadoutElement->stripNumber(positionOnSurfaceUnprojected, layerID)
+        );
+	      m_exitcode = 2;
+	      if(m_writeOutputFile) m_ntuple->Fill();
+	      ATH_MSG_DEBUG( "m_exitcode = 2 " );
+	      continue;
       }
       
       // Re-definition Of ID
@@ -1003,18 +1003,17 @@ StatusCode MM_DigitizationTool::doDigitization() {
       //
       ////////////////////////////////////////////////////////////////////
       
-      
       ////////////////////////////////////////////////////////////////////
       //
       // Strip Response Simulation For This Hit
       //
-      
       const MM_DigitToolInput stripDigitInput( stripNumber,
 					       distToChannel,
 					       inAngle_XZ,
 					       inAngle_YZ,
 					       localMagneticField,
-					       detectorReadoutElement->numberOfStrips(layerID),
+					       detectorReadoutElement->numberOfMissingBottomStrips(layerID),
+					       detectorReadoutElement->numberOfStrips(layerID)-detectorReadoutElement->numberOfMissingTopStrips(layerID),
 					       m_idHelper->gasGap(layerID),
 					       m_eventTime+m_globalHitTime
 					       );
@@ -1050,24 +1049,22 @@ StatusCode MM_DigitizationTool::doDigitization() {
       MM_ElectronicsToolInput stripDigitOutput( tmpStripOutput.NumberOfStripsPos(), tmpStripOutput.chipCharge(), tmpStripOutput.chipTime(), digitID , hit.kineticEnergy());
       
       // This block is purely validation
-      if (stripNumber!=1){ // Extra if statement from quick fix from deadstrip = #1
-	for(size_t i = 0; i<tmpStripOutput.NumberOfStripsPos().size(); i++){
-	  int tmpStripID = tmpStripOutput.NumberOfStripsPos().at(i);
-	  bool isValid;
-	  Identifier cr_id = m_idHelper->channelID(stName, m_idHelper->stationEta(layerID), m_idHelper->stationPhi(layerID), m_idHelper->multilayer(layerID), m_idHelper->gasGap(layerID), tmpStripID, true, &isValid);
-	  if (!isValid) {
-	    ATH_MSG_WARNING( "MicroMegas digitization: failed to create a valid ID for (chip response) strip n. " << tmpStripID << "; associated positions will be set to 0.0." );
-	  } else {
-	    Amg::Vector2D cr_strip_pos(0., 0.);
-	    if ( !detectorReadoutElement->stripPosition(cr_id,cr_strip_pos) ) {
-	      ATH_MSG_WARNING("MicroMegas digitization: failed to associate a valid local position for (chip response) strip n. "
-			      << tmpStripID
-			      << "; associated positions will be set to 0.0."
+	    for(size_t i = 0; i<tmpStripOutput.NumberOfStripsPos().size(); i++){
+	      int tmpStripID = tmpStripOutput.NumberOfStripsPos().at(i);
+	      bool isValid;
+	      Identifier cr_id = m_idHelper->channelID(stName, m_idHelper->stationEta(layerID), m_idHelper->stationPhi(layerID), m_idHelper->multilayer(layerID), m_idHelper->gasGap(layerID), tmpStripID, true, &isValid);
+	      if (!isValid) {
+	        ATH_MSG_WARNING( "MicroMegas digitization: failed to create a valid ID for (chip response) strip n. " << tmpStripID << "; associated positions will be set to 0.0." );
+	      } else {
+	        Amg::Vector2D cr_strip_pos(0., 0.);
+	        if ( !detectorReadoutElement->stripPosition(cr_id,cr_strip_pos) ) {
+	          ATH_MSG_WARNING("MicroMegas digitization: failed to associate a valid local position for (chip response) strip n. "
+			        << tmpStripID
+			        << "; associated positions will be set to 0.0."
 			      );
+	        }
+	      }
 	    }
-	  }
-	}
-      }
       
       
       v_stripDigitOutput.push_back(stripDigitOutput);
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
index 2f10fa07395..acf4e6c33b3 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
@@ -8,7 +8,7 @@ using namespace std;
 
 MM_StripResponse::MM_StripResponse() {}
 
-MM_StripResponse::MM_StripResponse(std::vector<MM_IonizationCluster> IonizationClusters, float timeResolution, float stripPitch, int stripID, int maxstripID) : m_timeResolution(timeResolution), m_stripPitch(stripPitch), m_stripID(stripID), m_maxstripID(maxstripID) {
+MM_StripResponse::MM_StripResponse(std::vector<MM_IonizationCluster> IonizationClusters, float timeResolution, float stripPitch, int stripID, int minstripID, int maxstripID) : m_timeResolution(timeResolution), m_stripPitch(stripPitch), m_stripID(stripID), m_minstripID(minstripID), m_maxstripID(maxstripID) {
 
 	for (auto& IonizationCluster : IonizationClusters)
 		for (auto& Electron : IonizationCluster.getElectrons())
@@ -45,8 +45,10 @@ void MM_StripResponse::calculateTimeSeries(float /*thetaD*/, int /*gasgap*/) {
 		}
 		else stripVal = m_stripID;
 
-		if (stripVal < 0 || stripVal > m_maxstripID) stripVal = -1;
-		(m_stripCharges[timeBin])[stripVal] += Electron->getCharge();
+		// Only add the strips that are either read out, or can cross talk to the read out strips
+    if (stripVal < m_minstripID-2 || stripVal > m_maxstripID+1) stripVal = -1;
+    if (stripVal > 0)
+      (m_stripCharges[timeBin])[stripVal] += Electron->getCharge();
 
 	}
 }
@@ -70,14 +72,16 @@ void MM_StripResponse::simulateCrossTalk(float crossTalk1, float crossTalk2) {
 
 				if (stripChargeVal==0.) continue;
 
-				if (stripVal-1 > -1) (m_stripCharges[timeBin])[stripVal-1] += stripChargeVal * crossTalk1;
-				if (stripVal+1 > -1) (m_stripCharges[timeBin])[stripVal+1] += stripChargeVal * crossTalk1;
-				(m_stripCharges[timeBin])[stripVal] -= stripChargeVal * crossTalk1 * ( (stripVal-1 > -1) + (stripVal+1 > -1) );
+        // Allow crosstalk between strips that exist.
+        // Will check for read out strips in calculateSummaries function
+				if (stripVal-1 > 0) (m_stripCharges[timeBin])[stripVal-1] += stripChargeVal * crossTalk1;
+				if (stripVal+1 < m_maxstripID) (m_stripCharges[timeBin])[stripVal+1] += stripChargeVal * crossTalk1;
+				(m_stripCharges[timeBin])[stripVal] -= stripChargeVal * crossTalk1 * ( (stripVal-1 > 0) + (stripVal+1 < m_maxstripID) );
 
 				if (crossTalk2 > 0.){
-					if (stripVal-2 > -1) (m_stripCharges[timeBin])[stripVal-2] += stripChargeVal * crossTalk2;
-					if (stripVal+2 > -1) (m_stripCharges[timeBin])[stripVal+2] += stripChargeVal * crossTalk2;
-					(m_stripCharges[timeBin])[stripVal] -= stripChargeVal * crossTalk2 * ( (stripVal-2 > -1) + (stripVal+2 > -1) );
+					if (stripVal-2 > 0) (m_stripCharges[timeBin])[stripVal-2] += stripChargeVal * crossTalk2;
+					if (stripVal+2 < m_maxstripID) (m_stripCharges[timeBin])[stripVal+2] += stripChargeVal * crossTalk2;
+					(m_stripCharges[timeBin])[stripVal] -= stripChargeVal * crossTalk2 * ( (stripVal-2 > 0) + (stripVal+2 < m_maxstripID) );
 				}
 			}
 		}
@@ -92,6 +96,12 @@ void MM_StripResponse::calculateSummaries(float chargeThreshold) {
 		int timeBin = stripTimeSeries.first;
 		for (auto & stripCharge : stripTimeSeries.second ){
 			int stripVal = stripCharge.first;
+      // remove dead (missing) strips
+      // First active strip starts at m_minstripID
+      // Last active strip numbrer is maxStripID-1
+      if (stripVal < m_minstripID || stripVal > m_maxstripID-1) continue;
+      // remove PCB gap strips
+      if (stripVal == 1023 || stripVal == 1024 || stripVal == 2047 || stripVal == 2048 || stripVal == 3071 || stripVal == 3072 || stripVal == 4095 || stripVal == 4096) continue;
 			float stripChargeVal = stripCharge.second;
 			if(stripChargeVal < chargeThreshold) continue;
 			
@@ -105,13 +115,13 @@ void MM_StripResponse::calculateSummaries(float chargeThreshold) {
 				}
 			}
 			if(!found){ // 	// strip not in vector, add new entry
-				m_v_strip.push_back(stripVal);
-				vector<float> qTemp;
-				qTemp.push_back(stripChargeVal);
-				m_v_stripTotalCharge.push_back(qTemp);
-				vector<float> tTemp;
-				tTemp.push_back(timeBin*m_timeResolution);
-				m_v_stripTimeThreshold.push_back(tTemp);
+        m_v_strip.push_back(stripVal);
+        vector<float> qTemp;
+        qTemp.push_back(stripChargeVal);
+        m_v_stripTotalCharge.push_back(qTemp);
+        vector<float> tTemp;
+        tTemp.push_back(timeBin*m_timeResolution);
+        m_v_stripTimeThreshold.push_back(tTemp);
 			}
 		}
 	}
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
index 1e7d8e4657c..b1748effff9 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
@@ -149,6 +149,7 @@ MM_StripToolOutput MM_StripsResponseSimulation::GetResponseFrom(const MM_DigitTo
 		digiInput.stripIDLocal(),
 		digiInput.incomingAngleXZ(), //degrees
 		digiInput.incomingAngleYZ(), //degrees
+		digiInput.stripMinID(),
 		digiInput.stripMaxID(),
 		digiInput
 		);
@@ -169,6 +170,7 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 											const int & stripID,
 											const float & incidentAngleXZ,
 											const float & incidentAngleYZ,
+											const int & stripMinID,
 											const int & stripMaxID,
 											const MM_DigitToolInput & digiInput)
 {
@@ -281,7 +283,7 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 
 	float timeresolution = 0.01; //ns
 
-	MM_StripResponse stripResponseObject(m_IonizationClusters, timeresolution, m_pitch, stripID, stripMaxID);
+	MM_StripResponse stripResponseObject(m_IonizationClusters, timeresolution, m_pitch, stripID, stripMinID, stripMaxID);
 	stripResponseObject.timeOrderElectrons();
 	stripResponseObject.calculateTimeSeries(incidentAngleXZ, digiInput.gasgap());
 	stripResponseObject.simulateCrossTalk( m_crossTalk1,  m_crossTalk2);
-- 
GitLab