Skip to content
Snippets Groups Projects
Commit 4697c1b1 authored by Eric Torrence's avatar Eric Torrence
Browse files

Merge remote-tracking branch 'upstream/master'

parents 787d1a8b 1c5608ed
No related branches found
No related tags found
16 merge requests!197Reco Script update for Segment Finder,!194Reco script update,!176Waveform Reco updates,!175Waveform reco update,!173Update production to latest master,!167Catch raw data format errors,!166Waveform raw integral bug,!165Reco scripts,!162Production Reco Scripts,!161Add waveform identifiers, configure digitizer channels,!159Move Waveform code to its own major area,!157Catch error on truncated events,!155Cluster Limit,!138Made NoisyStripFinder package,!134Added raw data reconstruction tests to pipeline,!123Waveform reconstruction bugfix
Hide whitespace changes
Inline Side-by-side
Showing
with 122 additions and 31 deletions
Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx
+ 113
− 25
View file @ 4697c1b1
......@@ -40,6 +40,7 @@ ClusterFitAlg::ClusterFitAlg(const std::string& name, ISvcLocator* pSvcLocator)
// Initialize method:
StatusCode ClusterFitAlg::initialize() {
ATH_MSG_INFO("ClusterFitAlg::initialize()");
ATH_MSG_ALWAYS("zCenter: " << m_zCenter[0] << " / " << m_zCenter[1] << " / " << m_zCenter[2] << " / " << m_zCenter[3]);
ATH_CHECK(m_faserTriggerDataKey.initialize());
ATH_CHECK(m_clusterContainerKey.initialize());
......@@ -76,9 +77,11 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
SG::ReadHandle<xAOD::FaserTriggerData> triggerContainer{m_faserTriggerDataKey, ctx};
ATH_CHECK(triggerContainer.isValid());
uint32_t eventNumber = triggerContainer->eventId();
uint8_t triggerBits = triggerContainer->tap();
// if ((triggerBits & 0x4) == 0)
if (m_triggerMask != 0 && ((triggerBits & m_triggerMask) == 0))
// if (eventNumber != 10841)
{
ATH_CHECK(trackContainer.record(std::move(outputTracks)));
return StatusCode::SUCCESS;
......@@ -111,22 +114,65 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
const TrackerDD::SiDetectorElement* elem = (*clusters)->detectorElement();
if (elem != nullptr)
{
ATH_MSG_VERBOSE("( " << m_idHelper->wafer_hash(elem->identify()) << " = " << m_idHelper->layer(elem->identify()) << "/" << m_idHelper->eta_module(elem->identify()) << "/" << m_idHelper->phi_module(elem->identify()) << "|" << m_idHelper->side(elem->identify()) <<
ATH_MSG_VERBOSE("( " << m_idHelper->wafer_hash(elem->identify()) << " = " << m_idHelper->station(elem->identify()) << "/" << m_idHelper->layer(elem->identify()) << "/" << m_idHelper->eta_module(elem->identify()) << "/" << m_idHelper->phi_module(elem->identify()) << "|" << m_idHelper->side(elem->identify()) <<
" ): isStereo? " << elem->isStereo() << " sinStereo: " << elem->sinStereo() << " Swap phi? " << elem->swapPhiReadoutDirection() << " Depth, Eta, Phi directions: " <<
elem->hitDepthDirection() << " : " << elem->hitEtaDirection() << " : " << elem->hitPhiDirection() << " / zGlobal = " << (*clusters)->globalPosition().z());
stations.emplace(m_idHelper->station(elem->identify()));
Identifier id = elem->identify();
stations.emplace(m_idHelper->station(id));
clusterInfo* info = new clusterInfo { };
info->cluster = *clusters;
info->sinAlpha = elem->sinStereo();
info->cosAlpha = sqrt(1.0 - info->sinAlpha * info->sinAlpha);
double alpha = std::abs(asin(elem->sinStereo()));
int layer = m_idHelper->layer(id);
int eta = m_idHelper->eta_module(id); // -1 or +1
int etaIndex = (eta + 1)/2; // 0 or 1
int phi = m_idHelper->phi_module(id); // 0 to 3
int moduleIndex = ((etaIndex + phi)%2 == 1 ? phi : phi + 4); // 0 to 7
switch (moduleIndex)
{
case 0:
case 2:
{
// info->sinAlpha = sin(alpha);
// info->cosAlpha = -cos(alpha);
info->sinAlpha = -sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 1:
case 3:
{
info->sinAlpha = -sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 4:
case 6:
{
// info->sinAlpha = -sin(alpha);
// info->cosAlpha = -cos(alpha);
info->sinAlpha = sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 5:
case 7:
{
info->sinAlpha = sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
}
int side = m_idHelper->side(id);
if (side > 0) info->sinAlpha = - info->sinAlpha; // always reverse angle for back side wafer
info->u = (*clusters)->globalPosition().y() * info->cosAlpha + (*clusters)->globalPosition().x() * info->sinAlpha;
info->z = (*clusters)->globalPosition().z();
ATH_MSG_VERBOSE("Event: " << eventNumber << " Adding cluster " << layer << "/" << eta << "/" << phi << "/" << side << ", u= " << info->u << ", z= " << info->z << ", sin= " << info->sinAlpha << ", cos= " << info->cosAlpha );
info->sigmaSq = (*clusters)->localCovariance()(0,0);
if (info->sinAlpha < 0)
if (info->sinAlpha * info->cosAlpha < 0)
{
clusterMap[hash/2].first.push_back(info);
}
else if (info->sinAlpha > 0)
else if (info->sinAlpha * info->cosAlpha > 0)
{
clusterMap[hash/2].second.push_back(info);
}
......@@ -162,18 +208,21 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
std::vector<clusterInfo*>& stereoMinus = entry.second.first;
std::vector<clusterInfo*>& stereoPlus = entry.second.second;
ATH_MSG_VERBOSE("Module entry in layer " << layer << " has " << stereoMinus.size() << " negative and " << stereoPlus.size() << " positive stereo clusters");
if (stereoMinus.size() > 0 && stereoPlus.size() > 0 && !layerGood[layer])
if (stereoMinus.size() > 0 && stereoPlus.size() > 0)
{
nGoodLayers++;
layerGood[layer] = true;
if (!layerGood[layer])
{
nGoodLayers++;
layerGood[layer] = true;
}
for (clusterInfo* minus : stereoMinus)
{
for (clusterInfo* plus : stereoPlus)
{
layerCombo* combo = new layerCombo { };
combo->initialize();
combo->addCluster(plus, m_zCenter);
combo->addCluster(minus, m_zCenter);
combo->addCluster(plus, m_zCenter[thisStation]);
combo->addCluster(minus, m_zCenter[thisStation]);
layerCombos[layer].push_back( combo );
}
}
......@@ -218,11 +267,40 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
}
}
}
ATH_MSG_VERBOSE("best fit: (" << bestFit(0) << "+/-" << sqrt(bestCov(0,0)) << "," <<
ATH_MSG_VERBOSE("Event #" << eventNumber << " best fit: (" << bestFit(0) << "+/-" << sqrt(bestCov(0,0)) << "," <<
bestFit(1) << "+/-" << sqrt(bestCov(1,1)) << "," <<
bestFit(2) << "+/-" << sqrt(bestCov(2,2)) << "," <<
bestFit(3) << "+/-" << sqrt(bestCov(3,3)) << "); chi2 = " << bestChi2 );
Residuals(bestClusters);
if (eventNumber == 10841) std::cout << "fit = ParametricPlot3D[{" << bestFit(0) << " + (" << bestFit(2) << ") z, " << bestFit(1) << " + ( " << bestFit(3) << ") z, z}, {z, -50, 50}, Boxed->False];" << std::endl;
int iC = 0;
for (auto info : bestClusters)
{
if (info->sinAlpha * info->cosAlpha < 0) continue;
double dz = info->z - m_zCenter[thisStation];
double xFit = bestFit(0) + bestFit(2) * dz;
double yFit = bestFit(1) + bestFit(3) * dz;
double uFit = yFit * info->cosAlpha + xFit * info->sinAlpha;
ATH_MSG_VERBOSE("Event: " << eventNumber << " z: " << info->z << " z-zC: " << dz<< " uFit: " << uFit << " u: " << info->u);
auto stripEnds = info->cluster->detectorElement()->endsOfStrip(info->cluster->localPosition());
if (eventNumber == 10841) std::cout << "l" << iC << " = Line[{{" << stripEnds.first.x() << ", " << stripEnds.first.y() << ", " << stripEnds.first.z() - m_zCenter[thisStation] << "}, {" <<
stripEnds.second.x() << ", " << stripEnds.second.y() << ", " << stripEnds.second.z() - m_zCenter[thisStation] << "}}];" << std::endl;
iC++;
}
for (auto info : bestClusters)
{
if (info->sinAlpha * info->cosAlpha > 0) continue;
double dz = info->z - m_zCenter[thisStation];
double xFit = bestFit(0) + bestFit(2) * dz;
double yFit = bestFit(1) + bestFit(3) * dz;
double uFit = yFit * info->cosAlpha + xFit * info->sinAlpha;
ATH_MSG_VERBOSE("Event: " << eventNumber << " z: " << info->z << " z-zC: " << dz<< " uFit: " << uFit << " u: " << info->u);
auto stripEnds = info->cluster->detectorElement()->endsOfStrip(info->cluster->localPosition());
if (eventNumber == 10841) std::cout << "l" << iC << " = Line[{{" << stripEnds.first.x() << ", " << stripEnds.first.y() << ", " << stripEnds.first.z() - m_zCenter[thisStation] << "}, {" <<
stripEnds.second.x() << ", " << stripEnds.second.y() << ", " << stripEnds.second.z() - m_zCenter[thisStation] << "}}];" << std::endl;
iC++;
}
if (eventNumber == 10841) std::cout << "Show[Graphics3D[{l1, l2, l3, l4, l5, l6}, Boxed -> False], fit]" << std::endl;
//Residuals(bestClusters);
ATH_CHECK(AddTrack(outputTracks, bestFit, bestCov, bestClusters, bestChi2, bestClusters.size()-4));
m_chi2->Fill(bestChi2);
setFilterPassed(true, ctx);
......@@ -293,13 +371,15 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
Trk::FitQuality* q = new Trk::FitQuality {chi2, ndof};
DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
int station = m_idHelper->station(fitClusters[0]->cluster->detectorElement()->identify());
// translate parameters to nominal fit point
s->push_back( GetState(fitResult, fitCovariance, nullptr) );
s->push_back( GetState(fitResult, fitCovariance, nullptr, station) );
for (const clusterInfo* cInfo : fitClusters)
{
s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster) );
}
// for (const clusterInfo* cInfo : fitClusters)
// {
// s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
// }
// Create and store track
tracks->push_back(new Trk::Track(i, s , q));
......@@ -309,12 +389,13 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
Trk::TrackStateOnSurface*
ClusterFitAlg::GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
const Eigen::Matrix< double, 4, 4 >& fitCovariance,
const FaserSCT_Cluster* fitCluster ) const
const FaserSCT_Cluster* fitCluster, int station ) const
{
// position of fit point:
double zFit = m_zCenter;
// int station = m_idHelper->station(fitCluster->detectorElement()->identify());
double zFit = m_zCenter[station];
if (fitCluster != nullptr) zFit = fitCluster->globalPosition()[2];
Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - m_zCenter), fitResult[1] + fitResult[3] * (zFit - m_zCenter), zFit };
Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - m_zCenter[station]), fitResult[1] + fitResult[3] * (zFit - m_zCenter[station]), zFit };
double phi = atan2( fitResult[3], fitResult[2] );
double theta = atan( sqrt(fitResult[2]*fitResult[2] + fitResult[3]*fitResult[3]) );
double qoverp = 1.0/100000.0;
......@@ -352,17 +433,24 @@ ClusterFitAlg::ClusterFit(layerCombo* c0, layerCombo* c1, layerCombo* c2) const
{
sums[i] = c0->sums[i] + c1->sums[i] + c2->sums[i];
}
Amg::MatrixX s(5,5);
Eigen::Matrix< double, 5, 5 > s;
s << sums[1] , sums[2] , sums[4] , sums[5] , sums[10] ,
sums[2] , sums[3] , sums[5] , sums[6] , sums[11] ,
sums[4] , sums[5] , sums[7] , sums[8] , sums[12] ,
sums[5] , sums[6] , sums[8] , sums[9] , sums[13] ,
sums[10] , sums[11] , sums[12] , sums[13] , sums[14];
Amg::MatrixX s4(4,4);
Eigen::Matrix< double, 4, 4 > s4;
s4 = s.block<4,4>(0,0);
Eigen::Matrix< double, 4, 1 > v;
v << sums[10] , sums[11] , sums[12] , sums[13];
Eigen::Matrix< double, 4, 1 > x = s4.colPivHouseholderQr().solve(v);
ATH_MSG_VERBOSE(s4);
ATH_MSG_VERBOSE(v);
ATH_MSG_VERBOSE(s4.inverse());
// Eigen::Matrix< double, 4, 1 > x = s4.colPivHouseholderQr().solve(v);
Eigen::Matrix< double, 4, 1 > x = s4.fullPivLu().solve(v);
ATH_MSG_VERBOSE("Check solution");
ATH_MSG_VERBOSE(s4*x);
ATH_MSG_VERBOSE(v);
Eigen::Matrix< double, 5, 1 > x5;
x5 << x(0), x(1), x(2), x(3), -1;
double chi2 = x5.transpose()*s*x5;
......
Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h
+ 3
− 2
View file @ 4697c1b1
......@@ -173,7 +173,8 @@ class ClusterFitAlg : public AthReentrantAlgorithm, AthHistogramming
Trk::TrackStateOnSurface* GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
const Eigen::Matrix< double, 4, 4 >& fitCovariance,
const FaserSCT_Cluster* fitCluster ) const;
const FaserSCT_Cluster* fitCluster,
int station ) const;
std::tuple<Eigen::Matrix<double, 4, 1>,
Eigen::Matrix<double, 4, 4>,
......@@ -189,7 +190,7 @@ class ClusterFitAlg : public AthReentrantAlgorithm, AthHistogramming
SG::ReadHandleKey<FaserSCT_ClusterContainer> m_clusterContainerKey { this, "ClustersName", "SCT_ClusterContainer", "FaserSCT cluster container" };
SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "ClusterFit", "Output track collection name" };
DoubleProperty m_zCenter { this, "ZCenter", 2300.0, "Global z position at which to reconstruct track parameters"};
DoubleArrayProperty m_zCenter { this, "ZCenter", { -1452.2925, 47.7075 , 1237.7075, 2427.7075 }, "Global z position at which to reconstruct track parameters"};
UnsignedIntegerProperty m_triggerMask { this, "TriggerMask", 0x0, "Trigger mask to analyze (0 = pass all)" };
/// a handle on the Hist/TTree registration service
ServiceHandle<ITHistSvc> m_histSvc;
......
Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py
+ 2
− 2
View file @ 4697c1b1
......@@ -26,7 +26,7 @@ Configurable.configurableRun3Behavior = True
# Configure
ConfigFlags.Input.Files = [
'/eos/project-f/faser-commissioning/TI12Data/Run-001261/Faser-Physics-001261-00000.raw',
'/eos/project-f/faser-commissioning/TI12Data/Run-001332/Faser-Physics-001332-00000.raw',
#'/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00000.raw',
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00001.raw',
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00002.raw',
......@@ -92,7 +92,7 @@ ConfigFlags.Input.Files = [
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00062.raw'
]
#ConfigFlags.Output.ESDFileName = "run608.ESD.pool.root"
ConfigFlags.Output.ESDFileName = "run001261.ESD.pool.root"
ConfigFlags.Output.ESDFileName = "run001332.ESD.pool.root"
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-XXXX-XXX-XX" # Needed to bypass autoconfig, only the "OFLCOND" matters at the moment
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
......
Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx
+ 2
− 1
View file @ 4697c1b1
......@@ -593,7 +593,8 @@ namespace Tracker
// when this is the case here, and set hitsInThirdTimeBin to zero later on
//
double iphipitch = 1./element->phiPitch();
double shift = m_lorentzAngleTool->getLorentzShift(element->identifyHash());
// double shift = m_lorentzAngleTool->getLorentzShift(element->identifyHash());
double shift = 0;
double stripPitch = design->stripPitch();
bool badStripInClusterOnThisModuleSide = (idGroups.size() != tbinGroups.size());
bool rotate = element->design().shape() == TrackerDD::Trapezoid || element->design().shape() == TrackerDD::Annulus;
......
Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx
+ 2
− 1
View file @ 4697c1b1
......@@ -81,7 +81,8 @@ Tracker::FaserSCT_Cluster* TrackerClusterMakerTool::sctCluster(
const TrackerDD::SiDetectorElement* element,
int errorStrategy) const{
double shift = m_sctLorentzAngleTool->getLorentzShift(element->identifyHash());
// double shift = m_sctLorentzAngleTool->getLorentzShift(element->identifyHash());
double shift = 0;
// const InDetDD::SiLocalPosition& localPosition =
// InDetDD::SiLocalPosition(localPos[Trk::locY),
// localPos[Trk::locX)+shift,0);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment