diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.cxx index 5ce28622b55d6885c51b6d9b1f82087b03abe7d9..7e7d40ba22227fde387894493502fe974c6015e9 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.cxx +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.cxx @@ -4,6 +4,11 @@ #include "ConstraintAngleMMClusterBuilderTool.h" +#include <numeric> +#include <cmath> +#include <algorithm> +#include <memory> + #include "GaudiKernel/PhysicalConstants.h" #include "GaudiKernel/SystemOfUnits.h" @@ -12,19 +17,8 @@ #include "Minuit2/Minuit2Minimizer.h" #include "Math/IFunction.h" -#include <cmath> namespace{ - - static constexpr double const& reciprocalSpeedOfLight = Gaudi::Units::c_light * 1e-6; // mm/ns - static constexpr double const& sqrt12PitchSquared = 0.425*0.425/12; // pitch/sqrt(12) * pitch/sqrt(12) to dertermine the transversal uncertainty of the drifted electron - - template <class T> - std::string printVec(std::vector<T> vec){ - std::ostringstream outstream; - std::copy(vec.begin(),vec.end(),std::ostream_iterator<T>(outstream,", ")); - return outstream.str(); - } //class to provide fit function to minuit class MyFCN : public ROOT::Math::IBaseFunctionMultiDim { public: @@ -80,14 +74,7 @@ Muon::ConstraintAngleMMClusterBuilderTool::ConstraintAngleMMClusterBuilderTool(c AthAlgTool(t,n,p) { declareInterface<IMMClusterBuilderTool>(this); - declareProperty("vDrift", m_vDrift = 0.047); - declareProperty("t0", m_t0 = -100); declareProperty("nSigmaSelection",m_nSigmaSelection = 3); - declareProperty("longDiff",m_longDiff=0.019); //mm/mm - declareProperty("transDiff",m_transDiff=0.036); //mm/mm - declareProperty("ionUncertainty",m_ionUncertainty=4.0); //ns - declareProperty("scaleXerror",m_scaleXerror=1.); - declareProperty("scaleYerror",m_scaleYerror=1.); declareProperty("sTheta",m_sigmaTheta = 3); declareProperty("fitAngleCut",m_fitAngleCut = 10); //degree } @@ -154,7 +141,6 @@ const{ ATH_MSG_DEBUG("Hit channel "<< m_idHelperSvc->mmIdHelper().channel(prd.identify()) << " local positionX " << prd.localPosition().x() << " time " << prd.time() - << " corrected time " << prd.time()-prd.globalPosition().norm()*reciprocalSpeedOfLight+m_t0 << " angle to IP" << std::atan(prd.globalPosition().perp()/std::abs(prd.globalPosition().z())) << " "<< std::atan(prd.globalPosition().perp()/std::abs(prd.globalPosition().z()))/Gaudi::Units::degree << " stationEta " <<m_idHelperSvc->mmIdHelper().stationEta(prd.identify()) << " stationPhi "<< m_idHelperSvc->mmIdHelper().stationEta(prd.identify()) @@ -174,7 +160,7 @@ const{ } ATH_MSG_DEBUG("Found "<< idxClusters.size() <<" clusters"); for(const auto idxCluster:idxClusters){ - ATH_MSG_DEBUG("cluster: "<< printVec(idxCluster)); + ATH_MSG_DEBUG("cluster: "<< idxCluster); } //determine mean IP pointing theta and mean position per cluster @@ -216,8 +202,8 @@ const{ for(const auto &idx:idxClusters.at(i_cluster)){ double x = sign * (mmPrdsPerLayer.at(idx).localPosition().x() - clusterMeanX.at(i_cluster)); - double y = (mmPrdsPerLayer.at(idx).time()-mmPrdsPerLayer.at(idx).globalPosition().norm()*reciprocalSpeedOfLight+m_t0)*m_vDrift; - double sy = std::sqrt(m_ionUncertainty*m_ionUncertainty*m_vDrift*m_vDrift + m_longDiff*m_longDiff*y*y); + double y = mmPrdsPerLayer.at(idx).driftDist(); + double sy = std::sqrt(mmPrdsPerLayer.at(idx).localCovariance()(1,1)); ATH_MSG_VERBOSE("selection: x "<< x <<" y "<< y << " fUpper(x) " << fUpper->Eval(x) <<" fLower(x) "<< fLower->Eval(x) << " ylower " << y-m_nSigmaSelection*sy << "yupper" << y+m_nSigmaSelection*sy); //select points which are within x sigma of the road if(x<=0) ATH_MSG_VERBOSE("x<=0: "<< fUpper->Eval(x)-y<<" "<< y-fLower->Eval(x) << " x*sigma"<< m_nSigmaSelection*sy); @@ -235,7 +221,7 @@ const{ } } // end of probe intercept idxClusters.at(i_cluster) = idxBest; - ATH_MSG_DEBUG(" cluster after filtering "<< printVec(idxClusters.at(i_cluster))); + ATH_MSG_DEBUG(" cluster after filtering "<< idxClusters.at(i_cluster)); } // end of i_cluster for(int i_cluster = idxClusters.size()-1; i_cluster>=0; i_cluster--){ if(idxClusters.at(i_cluster).size()<2){ // delete to small clusters after pattern reco @@ -278,19 +264,17 @@ const{ for(uint i_strip=0; i_strip < idxCluster.size(); i_strip++){ double x = prdPerLayer.at(idxCluster.at(i_strip)).localPosition().x()-meanX; - double y = (prdPerLayer.at(idxCluster.at(i_strip)).time()-prdPerLayer.at(i_strip).globalPosition().norm()*reciprocalSpeedOfLight+m_t0)*m_vDrift; - double xerror = std::sqrt(sqrt12PitchSquared + m_transDiff*m_transDiff*y*y); - double yerror = std::sqrt(m_ionUncertainty*m_ionUncertainty*m_vDrift*m_vDrift + m_longDiff*m_longDiff*y*y); + double y = prdPerLayer.at(idxCluster.at(i_strip)).driftDist(); + double xerror = std::sqrt(prdPerLayer.at(idxCluster.at(i_strip)).localCovariance()(0,0)); + double yerror = std::sqrt(prdPerLayer.at(idxCluster.at(i_strip)).localCovariance()(1,1)); fitGraph->SetPoint(i_strip,x,y); fitGraph->SetPointError(i_strip, xerror, yerror); fitFunc.addPoint(x,y,xerror,yerror); stripsOfCluster.push_back(prdPerLayer.at(idxCluster.at(i_strip)).identify()); stripsOfClusterCharges.push_back(prdPerLayer.at(idxCluster.at(i_strip)).charge()); - stripsOfClusterTimes.push_back(prdPerLayer.at(idxCluster.at(i_strip)).time()-prdPerLayer.at(i_strip).globalPosition().norm()*reciprocalSpeedOfLight+m_t0); stripsOfClusterChannels.push_back(m_idHelperSvc->mmIdHelper().channel(prdPerLayer.at(idxCluster.at(i_strip)).identify())); - - + stripsOfClusterTimes.push_back(prdPerLayer.at(idxCluster.at(i_strip)).time()); } if(std::abs(ffit->GetParameter(0)-clusterTheta)> m_fitAngleCut * Gaudi::Units::degree){return StatusCode::FAILURE;}; // very loose cut for now diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.h b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.h index ab13a98721e890111255604aae10aae2f4412b94..1c105ec983d061aad375b1859a120800a4638af5 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.h +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ConstraintAngleMMClusterBuilderTool.h @@ -4,15 +4,16 @@ #ifndef CONSTRIANTANGLEMMCLUSTERBUILDERTOOL_h #define CONSTRIANTANGLEMMCLUSTERBUILDERTOOL_h +#include <string> +#include <vector> -#include "GaudiKernel/ServiceHandle.h" #include "MMClusterization/IMMClusterBuilderTool.h" #include "MuonPrepRawData/MMPrepData.h" #include "AthenaBaseComps/AthAlgTool.h" +#include "GaudiKernel/ServiceHandle.h" #include "MuonIdHelpers/IMuonIdHelperSvc.h" -#include <numeric> namespace Muon{ class ConstraintAngleMMClusterBuilderTool: virtual public IMMClusterBuilderTool, public AthAlgTool { @@ -38,30 +39,11 @@ namespace Muon{ StatusCode scanLayer(const std::vector<Muon::MMPrepData> &mmPrdsPerLayer,std::vector<std::vector<uint>> &idxClusters, std::vector<double> &clusterTheta)const ; StatusCode fitCluster(const std::vector<Muon::MMPrepData> &prdPerLayer,const std::vector<uint>& idxCluster,const double &clusterTheta,std::vector<Muon::MMPrepData*>& clustersVec) const; - double m_t0; - double m_vDrift; int m_nSigmaSelection; double m_sigmaTheta; - double m_transDiff,m_longDiff, m_ionUncertainty, m_scaleXerror, m_scaleYerror; double m_fitAngleCut; + }; // class ConstraintAngleMMClusterBuilderTool - - - - - - - - - }; +} // namespace Muon #endif - - - - - - - - -} diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.cxx index 80a09166f9a9e404a7d4e7f85a20936d43447d07..53fda070640dbeb502b32e96f0d63edb23fe87c5 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.cxx +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.cxx @@ -2,26 +2,21 @@ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ #include "ProjectionMMClusterBuilderTool.h" + +#include <cmath> +#include <algorithm> + #include "MuonPrepRawData/MMPrepData.h" #include "MuonReadoutGeometry/MuonDetectorManager.h" #include "MuonIdHelpers/MmIdHelper.h" #include "GaudiKernel/MsgStream.h" -#include <cmath> +#include "GaudiKernel/SystemOfUnits.h" namespace{ -template <class T> -std::string printVec(std::vector<T> vec){ - std::ostringstream outstream; - std::copy(vec.begin(),vec.end(),std::ostream_iterator<T>(outstream,", ")); - return outstream.str(); -} - // Parametrization of the strip error after the projection constexpr double stripErrorSlope = 0.2; constexpr double stripErrorIntercept = 0.15; - - } @@ -35,10 +30,9 @@ Muon::ProjectionMMClusterBuilderTool::ProjectionMMClusterBuilderTool(const std:: AthAlgTool(t,n,p) { declareInterface<IMMClusterBuilderTool>(this); - declareProperty("vDrift",m_vDrift = 0.048); - declareProperty("tmin", m_tmin=30); - declareProperty("tmax", m_tmax=130); - declareProperty("tOffset", m_tOffset=10.); + declareProperty("tmin", m_tmin=0.0); + declareProperty("tmax", m_tmax=5.0); + declareProperty("tOffset", m_tOffset=0); declareProperty("corr_p0",m_p0=2.48741e-01); declareProperty("corr_p1",m_p1=-4.54356e-01); declareProperty("corr_p2",m_p2=3.77574e-01); @@ -51,7 +45,6 @@ Muon::ProjectionMMClusterBuilderTool::ProjectionMMClusterBuilderTool(const std:: StatusCode Muon::ProjectionMMClusterBuilderTool::initialize() { ATH_CHECK(m_idHelperSvc.retrieve()); - ATH_MSG_INFO("Set vDrift to "<< m_vDrift <<" um/ns"); return StatusCode::SUCCESS; } @@ -83,22 +76,17 @@ for(const auto& prdsOfLayer : prdsPerLayer){ if(prdsOfLayer.size()<2) continue; for(const auto& prd:prdsOfLayer){ Identifier id_prd = prd.identify(); - //determine time of flight from strip position - double distToIp=prd.globalPosition().norm(); - double tof=distToIp/299.792; // divide by speed of light in mm/ns //determine angle to IP for debug reasons - double globalR=std::sqrt(std::pow(prd.globalPosition().x(),2)+std::pow(prd.globalPosition().y(),2)); - double Tan=globalR/prd.globalPosition().z(); - double angleToIp=std::atan(Tan)/(180./M_PI); - ATH_MSG_DEBUG("Hit channel: "<< m_idHelperSvc->mmIdHelper().channel(id_prd) <<" time "<< prd.time() << " localPosX "<< prd.localPosition().x() << " tof "<<tof <<" angleToIp " << angleToIp<<" gas_gap "<< m_idHelperSvc->mmIdHelper().gasGap(id_prd) << " multiplet " << m_idHelperSvc->mmIdHelper().multilayer(id_prd) << " stationname " <<m_idHelperSvc->mmIdHelper().stationName(id_prd) << " stationPhi " <<m_idHelperSvc->mmIdHelper().stationPhi(id_prd) << " stationEta "<<m_idHelperSvc->mmIdHelper().stationEta(id_prd)); + double Tan=prd.globalPosition().perp() / prd.globalPosition().z(); + double angleToIp=std::atan(Tan)/Gaudi::Units::degree; + ATH_MSG_DEBUG("Hit channel: "<< m_idHelperSvc->mmIdHelper().channel(id_prd) <<" time "<< prd.time() << " drift dist " << prd.driftDist()<< " localPosX "<< prd.localPosition().x() <<" angleToIp " << angleToIp<<" gas_gap "<< m_idHelperSvc->mmIdHelper().gasGap(id_prd) << " multiplet " << m_idHelperSvc->mmIdHelper().multilayer(id_prd) << " stationname " <<m_idHelperSvc->mmIdHelper().stationName(id_prd) << " stationPhi " <<m_idHelperSvc->mmIdHelper().stationPhi(id_prd) << " stationEta "<<m_idHelperSvc->mmIdHelper().stationEta(id_prd)); } - //vector<double> xpos,time,corr,newXpos,charge,globalR; std::vector<double> v_posxc,v_cor; StatusCode sc=calculateCorrection(prdsOfLayer,v_posxc,v_cor); ATH_MSG_DEBUG("calculated correction"); - ATH_MSG_DEBUG("v_posxc: "<<printVec(v_posxc)); - ATH_MSG_DEBUG("v_cor: "<<printVec(v_cor)); + ATH_MSG_DEBUG("v_posxc: "<<v_posxc); + ATH_MSG_DEBUG("v_cor: "<<v_cor); if(sc.isFailure()) return sc; std::vector<int> flag(prdsOfLayer.size(),0); // flag for already used strips (1) or unused strips(0) @@ -109,8 +97,8 @@ for(const auto& prdsOfLayer : prdsPerLayer){ double xmean,xerr,qtot; sc = doFineScan(flag,v_posxc,v_cor,idx_selected); - ATH_MSG_DEBUG("flag"<<printVec(flag)); - ATH_MSG_DEBUG("idx_selected"<<printVec(idx_selected)); + ATH_MSG_DEBUG("flag"<<flag); + ATH_MSG_DEBUG("idx_selected"<<idx_selected); ATH_MSG_DEBUG("did fine scan and selected strips: "<< idx_selected.size() << " Status code: " << sc); if(sc.isFailure()){ATH_MSG_DEBUG("Fine scan failed"); break;} //returns FAILURE if less then three strips in a cluster @@ -140,16 +128,13 @@ StatusCode Muon::ProjectionMMClusterBuilderTool::calculateCorrection( const std::vector<Muon::MMPrepData> &prdsOfLayer,std::vector<double>& v_posxc,std::vector<double>& v_cor)const { for(const auto& prd : prdsOfLayer){ double posx = prd.localPosition().x(); - double posz = std::fabs(prd.globalPosition().z()); // use std::fabs of z to not break symetrie for A and C side - double tm = prd.time()-prd.globalPosition().norm()/299.792-m_t0; - //std::printf("Time after tof: %.2f tmin: %.1f tmax: %.1f \n",tm,m_tmin,m_tmax); + double tm = prd.driftDist(); if(tm<m_tmin) tm=m_tmin; if(tm>m_tmax) tm=m_tmax; - double globalR=std::sqrt(std::pow(prd.globalPosition().x(),2) + std::pow(prd.globalPosition().y(),2)); - double cor=m_vDrift*(tm-m_tmin-m_tOffset)*(1.0*globalR/posz); - //cor*=(m_idHelperSvc->mmIdHelper().gasGap(prd.identify())%2==1? 1:-1);//correct for the local drift time directions. Gas gap 2 and 4 have "negativ" time, therefore correction is applied negative + double cor=(tm-m_tmin-m_tOffset)*(1.0*prd.globalPosition().perp()/std::fabs(prd.globalPosition().z())); + cor*=(m_idHelperSvc->mmIdHelper().gasGap(prd.identify())%2==1? 1:-1);//correct for the local drift time directions. Gas gap 2 and 4 have "negativ" time, therefore correction is applied negative - ATH_MSG_DEBUG("globalR: "<< globalR << " Time: " << tm-m_tmin-m_tOffset << "globalZ: " << posz << " R/Z: "<<1.0*globalR/posz); + ATH_MSG_DEBUG("globalR: "<< prd.globalPosition().perp() << " drift dist: " << tm-m_tmin-m_tOffset << "globalZ: " << prd.globalPosition().z() << " R/Z: "<< 1.0*prd.globalPosition().perp()/std::fabs(prd.globalPosition().z())); v_cor.push_back(cor); // v_posxc.push_back(posx + cor); @@ -230,7 +215,7 @@ StatusCode Muon::ProjectionMMClusterBuilderTool::writeNewPrd(std::vector<Muon::M double meanTime=0; for(const auto& id_goodStrip:idx_selected){ stripsOfCluster.push_back(prdsOfLayer.at(id_goodStrip).identify()); - stripsOfClusterDriftTime.push_back(int(prdsOfLayer.at(id_goodStrip).time()-prdsOfLayer.at(id_goodStrip).globalPosition().norm()/299.792-m_t0)); + stripsOfClusterDriftTime.push_back(int(prdsOfLayer.at(id_goodStrip).time())); stripsOfClusterCharge.push_back(int(prdsOfLayer.at(id_goodStrip).charge())); stripsOfClusterStripNumber.push_back(m_idHelperSvc->mmIdHelper().channel(prdsOfLayer.at(id_goodStrip).identify())); diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.h b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.h index 8084bdef6bd9ba46250cf1a1abdcd8d4df880058..2085327f52233d732101eef9188cd29c833348ff 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.h +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/ProjectionMMClusterBuilderTool.h @@ -4,6 +4,10 @@ #ifndef ProjectionMMClusterBuilderTool_h #define ProjectionMMClusterBuilderTool_h +#include <numeric> +#include <string> +#include <vector> + #include "MMClusterization/IMMClusterBuilderTool.h" #include "MuonPrepRawData/MMPrepData.h" #include "AthenaBaseComps/AthAlgTool.h" @@ -11,10 +15,6 @@ #include "GaudiKernel/ServiceHandle.h" #include "MuonIdHelpers/IMuonIdHelperSvc.h" - -#include <numeric> - - class MmIdHelper; namespace MuonGM { @@ -47,7 +47,7 @@ namespace Muon ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}; - double m_vDrift,m_tmin,m_tmax,m_tOffset; + double m_tmin,m_tmax,m_tOffset; double m_p0,m_p1,m_p2; //correction factors for charge dependence int m_t0; diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.cxx index 7f71354cd957c5c8661be1e4bf1f56136fa9f156..192b462a07d5aeda0db53f3e41e2c90f7e462f43 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.cxx +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.cxx @@ -3,21 +3,15 @@ */ #include "UTPCMMClusterBuilderTool.h" -#include "GaudiKernel/MsgStream.h" #include <cmath> +#include <numeric> -namespace{ - -static constexpr double const& toRad = M_PI/180; +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/SystemOfUnits.h" -template <class T> -std::string printVec(std::vector<T> vec){ - std::ostringstream outstream; - std::copy(vec.begin(),vec.end(),std::ostream_iterator<T>(outstream,", ")); - return outstream.str(); -} +namespace { bool sortTracks(std::tuple<int,double> &a, std::tuple<int,double> &b){ if(std::get<0>(a)==std::get<0>(b)){ @@ -47,27 +41,21 @@ Muon::UTPCMMClusterBuilderTool::UTPCMMClusterBuilderTool(const std::string& t, declareProperty("HoughDMin",m_dMin=0); declareProperty("HoughExpectedDriftRange",m_driftRange=12); //mm declareProperty("HoughDResolution",m_dResolution=0.125); - declareProperty("HoughMinCounts",m_houghMinCounts=3); - declareProperty("timeOffset",m_timeOffset=0); //ns - declareProperty("uTPCDHalf",m_dHalf=2.5); //mm + declareProperty("HoughMinCounts",m_houghMinCounts=3); declareProperty("HoughSelectionCut",m_selectionCut=1.0); //mm - declareProperty("vDrift",m_vDrift=0.047); //mm/us - declareProperty("longDiff",m_longDiff=0.019); //mm/mm - declareProperty("transDiff",m_transDiff=0.036); //mm/mm - declareProperty("ionUncertainty",m_ionUncertainty=4.0); //ns - declareProperty("scaleXerror",m_scaleXerror=1.); - declareProperty("scaleYerror",m_scaleYerror=1.); declareProperty("outerChargeRatioCut",m_outerChargeRatioCut=0.); //charge ratio cut to supress cross talk declareProperty("maxStripRemove",m_maxStripsCut=4); //max number of strips cut by cross talk cut declareProperty("digiHasNegativeAngle",m_digiHasNegativeAngles=true); declareProperty("scaleClusterError",m_scaleClusterError=2); + ATH_MSG_INFO("outerChargeRatioCut: " << m_outerChargeRatioCut); + ATH_MSG_INFO("maxStripRemove " << m_maxStripsCut); + ATH_MSG_INFO("scale cluster error: "<< m_scaleClusterError); } StatusCode Muon::UTPCMMClusterBuilderTool::initialize(){ ATH_CHECK( m_idHelperSvc.retrieve() ); - ATH_MSG_INFO("Set vDrift to "<< m_vDrift <<" um/ns"); return StatusCode::SUCCESS; } @@ -91,29 +79,20 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD for(const auto& MMprdsOfLayer:MMprdsPerLayer){ - std::vector<double> stripsTime,stripsLocalPosX; - std::vector<Identifier> strips_id; - std::vector<int> stripsCharge,stripsChannel; + std::vector<double> stripsLocalPosX; + stripsLocalPosX.reserve(MMprdsOfLayer.size()); - std::vector<int> flag=std::vector<int>(MMprdsOfLayer.size(),0); // vector of 0s and ones to indicate if strip was already used (0=unused;1=used;2=rejected as cross talk) for( auto MMprd: MMprdsOfLayer){ Identifier id_prd=MMprd.identify(); - strips_id.push_back(id_prd); int gasGap=m_idHelperSvc->mmIdHelper().gasGap(id_prd); stripsLocalPosX.push_back(((gasGap%2==0 && m_digiHasNegativeAngles) ? -1:1)*MMprd.localPosition().x()); //flip local positions for even gas gaps to reduce hough space to positiv angles - //determine time of flight from strip position - double distToIp=MMprd.globalPosition().norm(); - double tof=distToIp/299.792; // divide by speed of light in mm/ns + //determine angle to IP for debug reasons double globalR=std::sqrt(std::pow(MMprd.globalPosition().x(),2)+std::pow(MMprd.globalPosition().y(),2)); double Tan=globalR/MMprd.globalPosition().z(); - double angleToIp=std::atan(Tan)/toRad; - - stripsTime.push_back(MMprd.time()-tof+m_timeOffset); // use time minus tof minus vmm integration time as actual time - stripsChannel.push_back(m_idHelperSvc->mmIdHelper().channel(id_prd)); - stripsCharge.push_back(MMprd.charge()); + double angleToIp=std::atan(Tan)/Gaudi::Units::degree; - ATH_MSG_DEBUG("Hit channel: "<< m_idHelperSvc->mmIdHelper().channel(id_prd) <<" time "<< MMprd.time()-tof+m_timeOffset << " localPosX "<< MMprd.localPosition().x() << " tof "<<tof <<" angleToIp " << angleToIp<<" gas_gap "<< m_idHelperSvc->mmIdHelper().gasGap(id_prd) << " multiplet " << m_idHelperSvc->mmIdHelper().multilayer(id_prd) << " stationname " <<m_idHelperSvc->mmIdHelper().stationName(id_prd) << " stationPhi " <<m_idHelperSvc->mmIdHelper().stationPhi(id_prd) << " stationEta "<<m_idHelperSvc->mmIdHelper().stationEta(id_prd)); + ATH_MSG_DEBUG("Hit channel: "<< m_idHelperSvc->mmIdHelper().channel(id_prd) <<" time "<< MMprd.time() << " drift dist " << MMprd.driftDist() << " localPosX "<< MMprd.localPosition().x() <<" angleToIp " << angleToIp<<" gas_gap "<< m_idHelperSvc->mmIdHelper().gasGap(id_prd) << " multiplet " << m_idHelperSvc->mmIdHelper().multilayer(id_prd) << " stationname " <<m_idHelperSvc->mmIdHelper().stationName(id_prd) << " stationPhi " <<m_idHelperSvc->mmIdHelper().stationPhi(id_prd) << " stationEta "<<m_idHelperSvc->mmIdHelper().stationEta(id_prd)); } if(MMprdsOfLayer.size()<2) continue; @@ -121,12 +100,15 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD double meanX = std::accumulate(stripsLocalPosX.begin(),stripsLocalPosX.end(),0.0)/stripsLocalPosX.size(); ATH_MSG_DEBUG("Got mean element "<<meanX ); std::for_each(stripsLocalPosX.begin(),stripsLocalPosX.end(),[meanX](double& d){d-=meanX;}); - ATH_MSG_VERBOSE(printVec(stripsLocalPosX)); - + ATH_MSG_VERBOSE(stripsLocalPosX); + + // vector of 0s, ones and twos to indicate if strip was already used (0=unused;1=used;2=rejected as cross talk) + std::vector<int> flag=std::vector<int>(MMprdsOfLayer.size(),0); + while(true){ ATH_MSG_DEBUG("while true"); std::vector<int> idx_goodStrips; - StatusCode sc=runHoughTrafo(flag,stripsLocalPosX,stripsTime,idx_goodStrips); + StatusCode sc=runHoughTrafo(MMprdsOfLayer,stripsLocalPosX,flag,idx_goodStrips); ATH_MSG_DEBUG("Hough done"); //if(sc.isFailure()) break; @@ -142,12 +124,12 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD //remove strips induced by crosstalk till no more strip was removed int nStripsCutByCrosstalk=0; while(!applyCrossTalkCut(idx_goodStrips,MMprdsOfLayer,flag,nStripsCutByCrosstalk).isFailure()){} - + ATH_MSG_DEBUG("After cutting CT idx_goddStrips " << idx_goodStrips<< " flag " << flag ); double localClusterPosition=-9999; double sigmaLocalClusterPosition=0; double finalFitAngle,finalFitChiSqProb; - sc=finalFit(MMprdsOfLayer,stripsTime,idx_goodStrips,localClusterPosition, sigmaLocalClusterPosition,finalFitAngle,finalFitChiSqProb); + sc=finalFit(MMprdsOfLayer,idx_goodStrips,localClusterPosition, sigmaLocalClusterPosition,finalFitAngle,finalFitChiSqProb); ATH_MSG_DEBUG("final fit done"); if(sc.isFailure()) break; @@ -157,6 +139,11 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD std::vector<short int> stripsOfClusterTimes; std::vector<int> stripsOfClusterCharges; + stripsOfCluster.reserve(idx_goodStrips.size()); + stripsOfClusterChannels.reserve(idx_goodStrips.size()); + stripsOfClusterTimes.reserve(idx_goodStrips.size()); + stripsOfClusterCharges.reserve(idx_goodStrips.size()); + ATH_MSG_DEBUG("Found good Strips: "<< idx_goodStrips.size()); @@ -164,10 +151,10 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD ATH_MSG_DEBUG("Setting good strip: "<<idx<< " size of strips flag: "<<flag.size()); flag.at(idx)=1; ATH_MSG_DEBUG("Set Good strips"); - stripsOfCluster.push_back(strips_id.at(idx)); - stripsOfClusterChannels.push_back(stripsChannel.at(idx)); - stripsOfClusterTimes.push_back(stripsTime.at(idx)); - stripsOfClusterCharges.push_back(stripsCharge.at(idx)); + stripsOfCluster.push_back(MMprdsOfLayer.at(idx).identify()); + stripsOfClusterChannels.push_back(m_idHelperSvc->mmIdHelper().channel(MMprdsOfLayer.at(idx).identify())); + stripsOfClusterTimes.push_back(MMprdsOfLayer.at(idx).time()); + stripsOfClusterCharges.push_back(MMprdsOfLayer.at(idx).charge()); } Amg::MatrixX* covN = new Amg::MatrixX(1,1); covN->coeffRef(0,0)=sigmaLocalClusterPosition; @@ -208,7 +195,8 @@ StatusCode Muon::UTPCMMClusterBuilderTool::getClusters(std::vector<Muon::MMPrepD return StatusCode::SUCCESS; } -StatusCode Muon::UTPCMMClusterBuilderTool::runHoughTrafo(std::vector<int>& flag,std::vector<double>& xpos, std::vector<double>& time,std::vector<int>& idx_selected)const{ +StatusCode Muon::UTPCMMClusterBuilderTool::runHoughTrafo(const std::vector<Muon::MMPrepData> &mmPrd, + std::vector<double>& xpos, std::vector<int>& flag, std::vector<int>& idx_selected)const{ ATH_MSG_DEBUG("beginning of runHoughTrafo() with: " << xpos.size() <<" hits" ); double maxX = *std::max_element(xpos.begin(),xpos.end()); double minX = *std::min_element(xpos.begin(),xpos.end()); @@ -216,7 +204,7 @@ StatusCode Muon::UTPCMMClusterBuilderTool::runHoughTrafo(std::vector<int>& flag, StatusCode sc = houghInitCummulator(h_hough,maxX,minX); if(sc.isFailure()) return sc; - sc = fillHoughTrafo(h_hough,flag,xpos,time); + sc = fillHoughTrafo(mmPrd, xpos, flag, h_hough); ATH_MSG_DEBUG("filled Hough"); if(sc.isFailure()) return sc; @@ -225,13 +213,12 @@ StatusCode Muon::UTPCMMClusterBuilderTool::runHoughTrafo(std::vector<int>& flag, if(sc.isFailure()) return sc; ATH_MSG_DEBUG("Found HoughPeaks"); - sc = selectTrack(houghPeaks,xpos,time,flag,idx_selected); - ATH_MSG_DEBUG("Selected track"); + sc = selectTrack(mmPrd, xpos, flag, houghPeaks, idx_selected); if(sc.isFailure()) return sc; return StatusCode::SUCCESS; } StatusCode Muon::UTPCMMClusterBuilderTool::houghInitCummulator(std::unique_ptr<TH2D>& h_hough,double xmax,double xmin)const{ - ATH_MSG_VERBOSE("xmax: "<< xmax <<" xmin: "<< xmin <<" m_dResolution "<< m_dResolution <<" m_alphaMin "<< m_alphaMin <<" m_alphaMax: "<< m_alphaMax <<" toRad: "<< toRad <<" m_alphaResolution: "<<m_alphaResolution); + ATH_MSG_VERBOSE("xmax: "<< xmax <<" xmin: "<< xmin <<" m_dResolution "<< m_dResolution <<" m_alphaMin "<< m_alphaMin <<" m_alphaMax: "<< m_alphaMax <<" m_alphaResolution: "<<m_alphaResolution); double dmax=std::max(std::fabs(xmin),std::fabs(xmax)); dmax=std::sqrt(std::pow(dmax,2)+std::pow(m_driftRange,2)); // rspace =sqrt(xmax*xmax+driftrange*driftrange) where driftrange is assumed to be 6mm int nbinsd = static_cast<int>(1.0*dmax*2/m_dResolution); @@ -239,24 +226,24 @@ StatusCode Muon::UTPCMMClusterBuilderTool::houghInitCummulator(std::unique_ptr<T ATH_MSG_VERBOSE("Hough Using nBinsA "<< nbinsa <<" nBinsd "<< nbinsd<<" dmax "<< dmax); - h_hough = std::unique_ptr<TH2D>(new TH2D("h_hough","h_hough",nbinsa,m_alphaMin*toRad,m_alphaMax*toRad,nbinsd,-dmax,dmax)); + h_hough = std::unique_ptr<TH2D>(new TH2D("h_hough","h_hough",nbinsa,m_alphaMin*Gaudi::Units::degree,m_alphaMax*Gaudi::Units::degree,nbinsd,-dmax,dmax)); return StatusCode::SUCCESS; } -StatusCode Muon::UTPCMMClusterBuilderTool::fillHoughTrafo(std::unique_ptr<TH2D>& h_hough,std::vector<int>& flag,std::vector<double>& xpos, std::vector<double>& time)const{ - for(size_t i_hit=0; i_hit<xpos.size(); i_hit++){ - if(flag.at(i_hit)!=0) continue; //skip hits which have been already used or been rejected as cross talk - double x=xpos.at(i_hit); - double y=time.at(i_hit)*m_vDrift; - for(int i_alpha=1;i_alpha<h_hough->GetNbinsX();i_alpha++){ - double alpha=h_hough->GetXaxis()->GetBinCenter(i_alpha); - double d= x*std::cos(alpha) - y*std::sin(alpha); - ATH_MSG_VERBOSE("Fill Hough alpha: "<< alpha <<" d "<< d << " x "<< x << " y "<< y); - h_hough->Fill(alpha,d); - } - } - return StatusCode::SUCCESS; +StatusCode Muon::UTPCMMClusterBuilderTool::fillHoughTrafo( + const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& xpos, std::vector<int>& flag, + std::unique_ptr<TH2D>& h_hough)const{ + for(size_t i_hit=0; i_hit<mmPrd.size(); i_hit++){ + if(flag.at(i_hit)!=0) continue; //skip hits which have been already used or been rejected as cross talk + for(int i_alpha=1;i_alpha<h_hough->GetNbinsX();i_alpha++){ + double alpha = h_hough->GetXaxis()->GetBinCenter(i_alpha); + double d = xpos.at(i_hit)*std::cos(alpha) - mmPrd.at(i_hit).driftDist()*std::sin(alpha); + ATH_MSG_VERBOSE("Fill Hough alpha: "<< alpha <<" d "<< d << " x "<< xpos.at(i_hit) << " y "<< mmPrd.at(i_hit).driftDist() ); + h_hough->Fill(alpha,d); + } + } + return StatusCode::SUCCESS; } StatusCode Muon::UTPCMMClusterBuilderTool::findAlphaMax(std::unique_ptr<TH2D>& h_hough, std::vector<std::tuple<double,double>> &maxPos)const{ @@ -269,7 +256,7 @@ StatusCode Muon::UTPCMMClusterBuilderTool::findAlphaMax(std::unique_ptr<TH2D>& h for(int i_binX=0; i_binX<=h_hough->GetNbinsX();i_binX++){ for(int i_binY=0; i_binY<=h_hough->GetNbinsY();i_binY++){ if(h_hough->GetBinContent(i_binX,i_binY)==cmax){ - ATH_MSG_DEBUG("Find Max Alpha: BinX "<< i_binX << " Alpha: "<< h_hough->GetXaxis()->GetBinCenter(i_binX)/toRad <<" BinY: "<< i_binY <<" over threshold: "<< h_hough->GetBinContent(i_binX,i_binY)); + ATH_MSG_DEBUG("Find Max Alpha: BinX "<< i_binX << " Alpha: "<< h_hough->GetXaxis()->GetBinCenter(i_binX)/Gaudi::Units::degree <<" BinY: "<< i_binY <<" over threshold: "<< h_hough->GetBinContent(i_binX,i_binY)); maxPos.push_back(std::make_tuple(h_hough->GetXaxis()->GetBinCenter(i_binX),h_hough->GetYaxis()->GetBinCenter(i_binY))); } } @@ -277,8 +264,8 @@ StatusCode Muon::UTPCMMClusterBuilderTool::findAlphaMax(std::unique_ptr<TH2D>& h return StatusCode::SUCCESS; } -StatusCode Muon::UTPCMMClusterBuilderTool::selectTrack(std::vector<std::tuple<double,double>> &tracks,std::vector<double>& xpos, std::vector<double>& time, - std::vector<int>& flag, std::vector<int> &idxGoodStrips)const{ +StatusCode Muon::UTPCMMClusterBuilderTool::selectTrack(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& xpos, std::vector<int>& flag, + std::vector<std::tuple<double,double>> &tracks, std::vector<int> &idxGoodStrips)const{ std::vector<double> chi2; std::vector<std::vector<int>> allGoodStrips; std::vector<int> ngoodStrips; @@ -295,27 +282,27 @@ StatusCode Muon::UTPCMMClusterBuilderTool::selectTrack(std::vector<std::tuple<do std::vector<int> indexUsedForDist,goodStrips; for(size_t i_hit=0;i_hit<xpos.size();i_hit++){ if (flag.at(i_hit) != 1){ - dist.push_back(time.at(i_hit)*m_vDrift-slope*xpos.at(i_hit)-intercept); + dist.push_back(mmPrd.at(i_hit).driftDist()-slope*xpos.at(i_hit)-intercept); indexUsedForDist.push_back(i_hit); } } - ATH_MSG_VERBOSE("indexUsedForDist "<< printVec(indexUsedForDist)); - ATH_MSG_VERBOSE("distVec: "<<printVec(dist)); + ATH_MSG_VERBOSE("indexUsedForDist "<< indexUsedForDist); + ATH_MSG_VERBOSE("distVec: "<< dist); dist.clear(); for(size_t i_hit: indexUsedForDist){ - double d=time.at(i_hit)*m_vDrift-slope*xpos.at(i_hit)-intercept; + double d = mmPrd.at(i_hit).driftDist() - slope * xpos.at(i_hit) - intercept; ATH_MSG_VERBOSE("dist: "<<d<<" for index "<<i_hit); if (std::fabs(d)<m_selectionCut){ - dist.push_back(std::pow(d/m_vDrift/10,2)); //determine dummy chi2 + dist.push_back(d*d/mmPrd.at(i_hit).localCovariance()(1,1)); //determine chi2 goodStrips.push_back(i_hit); } } if(true){ - ATH_MSG_DEBUG("Angle estimate ="<< std::get<0>(track)<<" "<<std::get<0>(track)/M_PI*180.); + ATH_MSG_DEBUG("Angle estimate ="<< std::get<0>(track)<<" "<<std::get<0>(track)/Gaudi::Units::degree); ATH_MSG_DEBUG("restimate ="<< std::get<1>(track)); ATH_MSG_DEBUG("slope estimate ="<< slope); ATH_MSG_DEBUG("intercept estimate ="<< intercept); - ATH_MSG_DEBUG("good strips ="<< printVec(goodStrips)); + ATH_MSG_DEBUG("good strips ="<< goodStrips); ATH_MSG_DEBUG("n good points ="<< goodStrips.size()); ATH_MSG_DEBUG("chi2 ="<< std::accumulate(dist.begin(),dist.end(),0.0)/(dist.size()-2)); } @@ -325,7 +312,6 @@ StatusCode Muon::UTPCMMClusterBuilderTool::selectTrack(std::vector<std::tuple<do trackQuality.push_back(std::make_tuple(goodStrips.size(),std::accumulate(dist.begin(),dist.end(),0.0)/(dist.size()-2))); } } - if(chi2.size()==0) return StatusCode::FAILURE; int trackIndex=std::min_element(trackQuality.begin(),trackQuality.end(),sortTracks)-trackQuality.begin(); idxGoodStrips=allGoodStrips.at(trackIndex); @@ -345,28 +331,35 @@ StatusCode Muon::UTPCMMClusterBuilderTool::transformParameters(double alpha, dou StatusCode Muon::UTPCMMClusterBuilderTool::applyCrossTalkCut(std::vector<int> &idxSelected,const std::vector<MMPrepData> &MMPrdsOfLayer,std::vector<int> &flag, int &nStripsCut)const{ int N=idxSelected.size(); bool removedStrip=false; - if (N<3 || nStripsCut>=m_maxStripsCut ) return StatusCode::FAILURE; - //reject outer strip for fit if charge ratio to second outer strip indicates strip charge is created by crosstalk only - if(1.0*MMPrdsOfLayer.at(idxSelected.at(N-1)).charge()/MMPrdsOfLayer.at(idxSelected.at(N-2)).charge()<m_outerChargeRatioCut){ + ATH_MSG_DEBUG("starting to remove cross talk strips "<< idxSelected << " flag "<< flag ); + if (N<3 || nStripsCut>=m_maxStripsCut ) {ATH_MSG_DEBUG("first cut failed");return StatusCode::FAILURE;} + //reject outer strip for fit if charge ratio to second outer strip indicates strip charge is created by crosstalk only + double ratioLastStrip = 1.0*MMPrdsOfLayer.at(idxSelected.at(N-1)).charge()/MMPrdsOfLayer.at(idxSelected.at(N-2)).charge(); + double ratioFirstStrip = 1.0*MMPrdsOfLayer.at(idxSelected.at(0)).charge()/MMPrdsOfLayer.at(idxSelected.at(1)).charge(); + ATH_MSG_DEBUG("ratioLastStrip " << ratioLastStrip << " ratioFirstStrip " << ratioFirstStrip); + if(ratioLastStrip<m_outerChargeRatioCut){ flag.at(idxSelected.at(N-1))=2; idxSelected.erase(idxSelected.begin()+(N-1)); removedStrip=true; nStripsCut++; + ATH_MSG_DEBUG("cutting last strip nStripCuts: "<< nStripsCut << " flag "<< flag << " idxSelected "<< idxSelected); } - if(idxSelected.size()<3 || nStripsCut>=m_maxStripsCut) return StatusCode::FAILURE; + if(idxSelected.size()<3 || nStripsCut>=m_maxStripsCut) {ATH_MSG_DEBUG("first cut failed");return StatusCode::FAILURE;} - if(1.0*MMPrdsOfLayer.at(idxSelected.at(0)).charge()/MMPrdsOfLayer.at(idxSelected.at(1)).charge()<m_outerChargeRatioCut){ + if(ratioFirstStrip<m_outerChargeRatioCut){ flag.at(idxSelected.at(0))=2; idxSelected.erase(idxSelected.begin()+0); removedStrip=true; nStripsCut++; + ATH_MSG_DEBUG("cutting first strip nStripCuts: "<< nStripsCut << " flag "<< flag << " idxSelected "<< idxSelected); } - if(nStripsCut>=m_maxStripsCut) return StatusCode::FAILURE; + if(nStripsCut>=m_maxStripsCut) {ATH_MSG_DEBUG("first cut failed");return StatusCode::FAILURE;} + ATH_MSG_DEBUG("return value "<< (!removedStrip ? "FAILURE":"SUCCESS")); return (!removedStrip ? StatusCode::FAILURE:StatusCode::SUCCESS); //return success if at least one strip was removed } -StatusCode Muon::UTPCMMClusterBuilderTool::finalFit(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& time, std::vector<int>& idxSelected, double& x0, double &sigmaX0 ,double &fitAngle, double &chiSqProb)const{ +StatusCode Muon::UTPCMMClusterBuilderTool::finalFit(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<int>& idxSelected,double& x0, double &sigmaX0, double &fitAngle, double &chiSqProb)const{ std::unique_ptr<TGraphErrors> fitGraph=std::unique_ptr<TGraphErrors>(new TGraphErrors()); std::unique_ptr<TF1> ffit=std::unique_ptr<TF1>(new TF1("ffit","pol1")); @@ -380,9 +373,9 @@ StatusCode Muon::UTPCMMClusterBuilderTool::finalFit(const std::vector<Muon::MMPr double xpos=mmPrd.at(idx).localPosition().x()-xmean; // substract mean value from x to center fit around 0 if(xmin>xpos) xmin = xpos; else if(xmax<xpos) xmax = xpos; - lf->AddPoint(&xpos,time.at(idx)*m_vDrift); - fitGraph->SetPoint(fitGraph->GetN(),xpos,time.at(idx)*m_vDrift); - fitGraph->SetPointError(fitGraph->GetN()-1,m_scaleXerror*std::pow(pow(0.425/3.464101615/*sqrt(12)*/,2)+std::pow(m_transDiff*time.at(idx)*m_vDrift,2),0.5), m_scaleYerror*std::pow(std::pow(m_ionUncertainty*m_vDrift,2)+std::pow(m_longDiff*time.at(idx)*m_vDrift,2),0.5)); + lf->AddPoint(&xpos, mmPrd.at(idx).driftDist()); + fitGraph->SetPoint(fitGraph->GetN(),xpos,mmPrd.at(idx).driftDist()); + fitGraph->SetPointError(fitGraph->GetN()-1, std::sqrt(mmPrd.at(idx).localCovariance()(0,0)),std::sqrt(mmPrd.at(idx).localCovariance()(1,1))); } lf->Eval(); @@ -399,7 +392,8 @@ StatusCode Muon::UTPCMMClusterBuilderTool::finalFit(const std::vector<Muon::MMPr double cov = (s->GetCovarianceMatrix())(0,1); ATH_MSG_DEBUG("covariance is: "<<cov); double fitSlope = ffit->GetParameter(1); - double interceptCorr = m_dHalf - ffit->GetParameter(0); + double dHalfGap = 2.5; // half the thickness of the gas gap + double interceptCorr = dHalfGap - ffit->GetParameter(0); x0 = interceptCorr/fitSlope; sigmaX0 = pow(m_scaleClusterError,2)*(pow(ffit->GetParError(0)/interceptCorr,2)+pow(ffit->GetParError(1)/fitSlope,2) - 2./(fitSlope*interceptCorr)*cov)*pow(x0,2); x0 += xmean; //adding back mean value after error was calculated diff --git a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.h b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.h index 81073b0dce2029315d8702cd84da97f7b77fb1d5..ef6bcf1c18c1c881ddb3522981bb5f1225dd429c 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.h +++ b/MuonSpectrometer/MuonReconstruction/MuonDataPrep/MMClusterization/src/UTPCMMClusterBuilderTool.h @@ -4,14 +4,17 @@ #ifndef UTPCMMClusterBuilderTool_h #define UTPCMMClusterBuilderTool_h -#include "GaudiKernel/ServiceHandle.h" +#include <tuple> +#include <vector> +#include <memory> +#include <string> + #include "MMClusterization/IMMClusterBuilderTool.h" #include "MuonPrepRawData/MMPrepData.h" #include "AthenaBaseComps/AthAlgTool.h" -#include "MuonIdHelpers/IMuonIdHelperSvc.h" -#include <numeric> -#include <cmath> +#include "GaudiKernel/ServiceHandle.h" +#include "MuonIdHelpers/IMuonIdHelperSvc.h" #include "TH2D.h" #include "TF1.h" @@ -32,7 +35,7 @@ namespace Muon public: /** Default constructor */ UTPCMMClusterBuilderTool(const std::string&, const std::string&, const IInterface*); - + /** Default destructor */ virtual ~UTPCMMClusterBuilderTool() = default; @@ -53,25 +56,24 @@ namespace Muon double m_dMin,m_dMax,m_dResolution,m_driftRange; int m_houghMinCounts; - double m_timeOffset,m_dHalf,m_vDrift,m_transDiff,m_longDiff, m_ionUncertainty, m_scaleXerror, m_scaleYerror; double m_outerChargeRatioCut; int m_maxStripsCut; bool m_digiHasNegativeAngles; float m_scaleClusterError; - StatusCode runHoughTrafo(std::vector<int>& flag,std::vector<double>& xpos, std::vector<double>& time,std::vector<int>& idx_selected)const; - StatusCode fillHoughTrafo(std::unique_ptr<TH2D>& cummulator,std::vector<int>& flag, std::vector<double>& xpos, std::vector<double>& time)const; - StatusCode houghInitCummulator(std::unique_ptr<TH2D>& cummulator,double xmax,double xmin)const; + StatusCode runHoughTrafo(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& xpos, std::vector<int>& flag, std::vector<int>& idx_selected)const; + StatusCode fillHoughTrafo(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& xpos, std::vector<int>& flag, std::unique_ptr<TH2D>& h_hough)const; + StatusCode houghInitCummulator(std::unique_ptr<TH2D>& cummulator, double xmax, double xmin)const; StatusCode findAlphaMax(std::unique_ptr<TH2D>& h_hough, std::vector<std::tuple<double,double>> &maxPos)const; - StatusCode selectTrack(std::vector<std::tuple<double,double>> &tracks,std::vector<double>& xpos, std::vector<double>& time,std::vector<int>& flag,std::vector<int>& idx_selected)const; + StatusCode selectTrack(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& xpos, std::vector<int>& flag, std::vector<std::tuple<double,double>> &tracks, std::vector<int> &idxGoodStrips)const; - StatusCode transformParameters(double alpha, double d, double dRMS, double& slope,double& intercept, double& interceptRMS)const; + StatusCode transformParameters(double alpha, double d, double dRMS, double& slope, double& intercept, double& interceptRMS)const; StatusCode applyCrossTalkCut(std::vector<int> &idxSelected,const std::vector<MMPrepData> &MMPrdsOfLayer,std::vector<int> &flag,int &nStripsCut)const; - StatusCode finalFit(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<double>& time, std::vector<int>& idxSelected,double& x0, double &sigmaX0, double &fitAngle, double &chiSqProb)const; -}; + StatusCode finalFit(const std::vector<Muon::MMPrepData> &mmPrd, std::vector<int>& idxSelected,double& x0, double &sigmaX0, double &fitAngle, double &chiSqProb)const; + }; -} +} // namespace Muon #endif