Skip to content
Snippets Groups Projects
Commit 5e45f700 authored by Alexandre Laurier's avatar Alexandre Laurier
Browse files

Add funtions to MM Digitization to add only read strips

parent 37f90987
No related branches found
No related tags found
No related merge requests found
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
......@@ -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; };
......
......@@ -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);
......
......@@ -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);
}
}
}
......
......@@ -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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment