runMdtGeoComparison.cxx File Reference

#include <GeoPrimitives/GeoPrimitivesHelpers.h>
#include <GeoPrimitives/GeoPrimitivesToStringConverter.h>
#include "GeoModelHelpers/TransformToStringConverter.h"
#include <MuonCablingData/MdtCablingData.h>
#include <MuonReadoutGeometryR4/MuonDetectorDefs.h>
#include <GaudiKernel/SystemOfUnits.h>
#include <iostream>
#include "Acts/Utilities/UnitVectors.hpp"
#include "Acts/Definitions/Units.hpp"
#include <PathResolver/PathResolver.h>
#include <TFile.h>
#include <TTreeReader.h>

Go to the source code of this file.

Classes
struct	MdtChamber
	Helper struct to represent a full Mdt chamber. More...
struct	MdtChamber::TubePositioning

Macros
#define	TEST_BASICPROP(attribute, propName)
#define	TEST_TUBEPROP(attribute, propName)

Functions
Amg::Vector3D	makeDir (const double theta, const double phi)
std::ostream &	operator<< (std::ostream &ostr, const MdtChamber &chamb)
	Translation of the station Index -> station Name.
std::set< MdtChamber >	readTreeDump (const std::string &inputFile)
int	main1 (int argc, char **argv)
int	main (int argc, char **argv)

Variables
constexpr double	tolerance = 10 * Gaudi::Units::micrometer

Macro Definition Documentation

TEST_BASICPROP

#define TEST_BASICPROP	(		attribute,
			propName )

Value:

if (std::abs(1.*test.attribute - 1.*reference.attribute) > tolerance) {           \
    std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": "<<test                    \
             <<" differs w.r.t "<<propName<<" "<< reference.attribute             \
             <<" (ref) vs. " <<test.attribute << " (test)" << std::endl;          \
    chamberOkay = false;                                                          \
}

Definition at line 113 of file runMdtGeoComparison.cxx.

#define TEST_BASICPROP(attribute, propName) \
    if (std::abs(1.*test.attribute - 1.*reference.attribute) > tolerance) {           \
        std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": "<<test                    \
                 <<" differs w.r.t "<<propName<<" "<< reference.attribute             \
                 <<" (ref) vs. " <<test.attribute << " (test)" << std::endl;          \
        chamberOkay = false;                                                          \
    }

TEST_TUBEPROP

#define TEST_TUBEPROP	(		attribute,
			propName )

Value:

if (std::abs(1.*refTube.attribute - 1.*testTube.attribute) > tolerance) {         \
    std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": In "<<test<<" the tubes (" \
             <<layer<<","<<std::setfill('0')<<std::setw(3)<<tube<<")"             \
             <<" differ w.r.t "<<propName<<". "                                   \
             << refTube.attribute <<" (ref) vs. " <<testTube.attribute            \
             << " (test)" << std::endl;                                           \
    chamberOkay = false;                                                          \
}

Definition at line 121 of file runMdtGeoComparison.cxx.

#define TEST_TUBEPROP(attribute, propName) \
    if (std::abs(1.*refTube.attribute - 1.*testTube.attribute) > tolerance) {         \
        std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": In "<<test<<" the tubes (" \
                 <<layer<<","<<std::setfill('0')<<std::setw(3)<<tube<<")"             \
                 <<" differ w.r.t "<<propName<<". "                                   \
                 << refTube.attribute <<" (ref) vs. " <<testTube.attribute            \
                 << " (test)" << std::endl;                                           \
        chamberOkay = false;                                                          \
    }

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 395 of file runMdtGeoComparison.cxx.

{
  int ret = 1;
  try {
    ret = main1 (argc, argv);
  }
  catch (const std::exception& e) {
    std::cerr << e.what() << "\n";
  }
  return ret;
}

main1()

int main1	(	int	argc,
		char **	argv )

check whether the files are xroot d -> otherwise call path resovler

Parse the tree dump

Start to loop over the chambers

We do not care whether the orientation of the coordinate system along the wire flips for negative chambers or not

The ultimate goal is to have the tube positioned at the same place. We maybe need the origin position later when we are adding the alignable transforms...

Check the tube transformations

Remember the tube staggering is in the (x-y) plane. Allow for deviations in the z-axis due to different cutouts

In cases where the tube coordinate systems are not aligned, there's no point in checking the position of the reaodout

Definition at line 249 of file runMdtGeoComparison.cxx.

                                   {
    std::string refFile{}, testFile{};
    
    for (int arg = 1; arg < argc; ++arg) {
       std::string the_arg{argv[arg]};
       if (the_arg == "--refFile" && arg +1 < argc) {
          refFile = std::string{argv[arg+1]};
          ++arg;
       } else if (the_arg == "--testFile" && arg + 1 < argc) {
            testFile = std::string{argv[arg+1]};
            ++arg;
       }
    }
    if (refFile.empty()) {
        std::cerr<<"Please parse the path of the reference file via --refFile "<<std::endl;
        return EXIT_FAILURE;
    }
    if (testFile.empty()) {
        std::cerr<<"Please parse the path of the test file via --testFile "<<std::endl;
        return EXIT_FAILURE;
    }
    if (!refFile.starts_with( "root://")) refFile = PathResolver::FindCalibFile(refFile);
    if (!testFile.starts_with( "root://")) testFile = PathResolver::FindCalibFile(testFile);
    std::set<MdtChamber> refChambers = readTreeDump(refFile);
    if (refChambers.empty()) {
        std::cerr<<"The file "<<refFile<<" should contain at least one chamber "<<std::endl;
        return EXIT_FAILURE;
    }
    std::set<MdtChamber> testChambers = readTreeDump(testFile);
    if (testChambers.empty()) {
        std::cerr<<"The file "<<testFile<<" should contain at least one chamber "<<std::endl;
        return EXIT_FAILURE;
    }
    int return_code = EXIT_SUCCESS;
    unsigned int goodChamb{0};
    for (const MdtChamber& reference : refChambers) {
        std::set<MdtChamber>::const_iterator test_itr = testChambers.find(reference);
        
        if (test_itr == testChambers.end()) {
            std::cerr<<"The chamber "<<reference<<" is not part of the testing "<<std::endl;
            return_code = EXIT_FAILURE;
            continue;
        }
        bool chamberOkay = true;
        const MdtChamber& test = {*test_itr};
        TEST_BASICPROP(numLayers, "number of layers");
        TEST_BASICPROP(numTubes, "number of tubes");
        TEST_BASICPROP(tubePitch, "tube pitch");
        TEST_BASICPROP(tubeRadius, "tube radius");
        
        const Amg::Transform3D distortion = test.geoModelTransform.inverse() * reference.geoModelTransform;
        bool flippedChamb = {reference.id.eta < 0 && Amg::doesNotDeform(distortion * Amg::getRotateX3D(M_PI))};
        if (!Amg::doesNotDeform(distortion) && !flippedChamb) {   
            std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": The chamber coordinate systems rotate differently for  "
                     <<reference<<". Difference in the coordinate transformation: "<<Amg::toString(distortion)
                     <<" --- refTrf: "<<Amg::toString(reference.geoModelTransform)
                     <<" --- testTrf: "<<Amg::toString(test.geoModelTransform)<<std::endl;
            chamberOkay = false;            
        }
        if (distortion.translation().mag() > tolerance) {
            std::cout<<"The origins of the chamber coordinate systems are not exactly at the same point for "
                     <<reference<<". Translation shift: "<<Amg::toString(distortion.translation(), 2)<<std::endl;
        }
        bool stagFailure{false}, alignFailure{false}, readoutOrient{false};            
        for (unsigned int layer = 1; layer<= std::min(reference.numLayers, test.numLayers) ; ++layer) {
            for (unsigned int tube = 1; tube <= std::min(reference.numTubes, test.numTubes); ++tube) {
                using TubePositioning = MdtChamber::TubePositioning;
                const TubePositioning& refTube = reference.getTube(layer, tube);
                const TubePositioning& testTube = test.getTube(layer, tube);
                const Amg::Transform3D tubeDistortion = testTube.localToGlobal.inverse() * refTube.localToGlobal;
                bool flippedTube{reference.id.eta < 0 && Amg::doesNotDeform(tubeDistortion * Amg::getRotateX3D(M_PI))};
        
                if (!alignFailure && !(Amg::doesNotDeform(tubeDistortion)  || flippedTube)) {
                    std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": In chamber "<<reference<<" the tube reference systems for ("<<layer<<", "
                             <<std::setfill('0')<<std::setw(3)<<tube<<") are not exactly aligned. "<<GeoTrf::toString(tubeDistortion)<<std::endl;                   
                    alignFailure = true;
                }
                if (!stagFailure && tubeDistortion.translation().perp() > tolerance) {
                    std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": Misplaced staggering found in "<<reference<<" the tube ("<<layer
                             <<", "<<std::setfill('0')<<std::setw(3)<<tube<<") "
                             << Amg::toString(tubeDistortion.translation(), 3)<<", mag: "<<tubeDistortion.translation().mag()<<
                             ", perp: "<<tubeDistortion.translation().perp()<<std::endl;                    
                    if (tube > 1) stagFailure = true;
                }
                
                // TEST_TUBEPROP(tubeLength, "tube length");
                // TEST_TUBEPROP(wireLength, "wire length");
                TEST_TUBEPROP(activeLength, "active length");
                if (alignFailure || readoutOrient) continue;
                const Amg::Transform3D refSystem = refTube.localToGlobal.inverse();
 
                const Amg::Vector3D refRO = refSystem * refTube.readoutPos;
                const Amg::Vector3D testRO = refSystem* testTube.readoutPos;
                if (refRO.z()* testRO.z() < 0.){
                    std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": The readout is on different sites for chamber: "<<reference<<
                             ", layer "<<layer<<", tube: "<<std::setfill('0')<<std::setw(3)<<tube<<". " 
                             <<Amg::toString(refRO, 2)<<" vs. "<<Amg::toString(testRO)<<std::endl;
                    readoutOrient = true;
                }
            }
            if (stagFailure || alignFailure) {
                chamberOkay = false;
            }
        }
        if (!chamberOkay) {
            return_code = EXIT_FAILURE;
            continue;
        }
        const Amg::Transform3D alignableDistort = test.alignableTransform.inverse()*(reference.alignableTransform );
        if (!Amg::doesNotDeform(alignableDistort) || alignableDistort.translation().mag() > tolerance) {
            std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": The alignable nodes are at differnt places for  "
                     <<test<<". " <<GeoTrf::toString(alignableDistort, true)<<" chamber length: "<<
                      (reference.tubePitch * (1.*reference.numTubes + 0.5))<<std::endl;
            chamberOkay = false;
 
        } else {
            std::cout<<"runMdtGeoComparision() "<<__LINE__<<": Found perfect agreement between new & old geometry for "<<reference<<std::endl;
            ++goodChamb;
        }
    }
    for (const MdtChamber& test : testChambers) {
        if (!refChambers.count(test)) {
            std::cerr<<"runMdtGeoComparision() "<<__LINE__<<": "<<test<<" is only in the test set."<<std::endl;
            return_code = EXIT_FAILURE;
        }
    }
    std::cout<<"runMdtGeoComparision() "<<__LINE__<<": "<<goodChamb<<"/"<<refChambers.size()<<" chambers are in perfect agreement. "<<std::endl;
 
    return return_code;
 
}

makeDir()

Amg::Vector3D makeDir	(	const double	theta,
		const double	phi )

Definition at line 30 of file runMdtGeoComparison.cxx.

                                                          {
    using namespace Acts::UnitLiterals;
    return Acts::makeDirectionFromPhiTheta(phi *1._degree, theta* 1._degree);
}

operator<<()

std::ostream & operator<<	(	std::ostream &	ostr,
		const MdtChamber &	chamb )

Translation of the station Index -> station Name.

Dictionary taken from https://gitlab.cern.ch/atlas/athena/-/blob/main/DetectorDescription/IdDictParser/data/IdDictMuonSpectrometer_R.09.03.xml

Definition at line 97 of file runMdtGeoComparison.cxx.

                                                                  {
    static const std::map<int, std::string> stationDict{
        {0, "BIL"}, {1, "BIS"}, {7, "BIR"},
        {2, "BML"}, {3, "BMS"}, {8, "BMF"}, {53, "BME"}, {54, "BMG"}, {52, "BIM"},
        {4, "BOL"}, {5, "BOS"}, {9, "BOF"}, {10, "BOG"},
        {6, "BEE"}, {14, "EEL"}, {15, "EES"}, 
        {13, "EIL"}, {49, "EIS"},
        {17, "EML"}, {18, "EMS"}, 
        {20, "EOL"}, {21, "EOS"}
    };
    ostr<<stationDict.at(chamb.id.stationIndex)<<" "<<chamb.id<<" "<<chamb.design;    
    return ostr;
} 

readTreeDump()

std::set< MdtChamber > readTreeDump

(

const std::string &

inputFile

)

Information to access each tube individually

Readout the information of each tube specifically

Definition at line 131 of file runMdtGeoComparison.cxx.

                                                            {
    std::set<MdtChamber> to_ret{};
    std::cout<<"Read the Mdt geometry tree dump from "<<inputFile<<std::endl;
    std::unique_ptr<TFile> inFile{TFile::Open(inputFile.c_str())};
    if (!inFile || !inFile->IsOpen()) {
        std::cerr<<__FILE__<<":"<<__LINE__<<" Failed to open "<<inputFile<<std::endl;
        return to_ret;
    }
    TTreeReader treeReader("MdtGeoModelTree", inFile.get());
    if (treeReader.IsInvalid()){
        std::cerr<<__FILE__<<":"<<__LINE__<<" The file "<<inputFile<<" does not contain the 'MdtGeoModelTree'"<<std::endl;
        return to_ret;
    }
 
    TTreeReaderValue<unsigned short> stationIndex{treeReader, "stationIndex"};
    TTreeReaderValue<short> stationEta{treeReader, "stationEta"};
    TTreeReaderValue<short> stationPhi{treeReader, "stationPhi"};
    TTreeReaderValue<short> stationML{treeReader, "stationMultiLayer"};
    TTreeReaderValue<std::string> chamberDesign{treeReader,"chamberDesign"};
    
    TTreeReaderValue<double> tubeRadius{treeReader, "tubeRadius"};
    TTreeReaderValue<double> tubePitch{treeReader, "tubePitch"}; 
    TTreeReaderValue<unsigned short> numTubes{treeReader, "numTubes"};
    TTreeReaderValue<unsigned short> numLayers{treeReader, "numLayers"};
    
    TTreeReaderValue<std::vector<float>> geoModelTransformX{treeReader, "GeoModelTransformX"};
    TTreeReaderValue<std::vector<float>> geoModelTransformY{treeReader, "GeoModelTransformY"};
    TTreeReaderValue<std::vector<float>> geoModelTransformZ{treeReader, "GeoModelTransformZ"};
 
    TTreeReaderValue<std::vector<float>> alignableNodeX{treeReader, "AlignableNodeX"};
    TTreeReaderValue<std::vector<float>> alignableNodeY{treeReader, "AlignableNodeY"};
    TTreeReaderValue<std::vector<float>> alignableNodeZ{treeReader, "AlignableNodeZ"};
 

    TTreeReaderValue<std::vector<unsigned short>> tubeLayer{treeReader, "tubeLayer"}; 
    TTreeReaderValue<std::vector<unsigned short>> tubeNumber{treeReader, "tubeNumber"};
 
    TTreeReaderValue<std::vector<double>> tubeLength{treeReader, "tubeLength"};
    TTreeReaderValue<std::vector<double>> activeTubeLength{treeReader, "activeTubeLength"};
    TTreeReaderValue<std::vector<double>> wireLength{treeReader, "wireLength"};
 
    TTreeReaderValue<std::vector<float>> tubeTransformTransX{treeReader, "tubeTransformTranslationX"};
    TTreeReaderValue<std::vector<float>> tubeTransformTransY{treeReader, "tubeTransformTranslationY"};
    TTreeReaderValue<std::vector<float>> tubeTransformTransZ{treeReader, "tubeTransformTranslationZ"};
    
    TTreeReaderValue<std::vector<float>> tubeTransformCol0Theta{treeReader, "tubeTransformLinearCol0Theta"};
    TTreeReaderValue<std::vector<float>> tubeTransformCol0Phi{treeReader, "tubeTransformLinearCol0Phi"};
 
    TTreeReaderValue<std::vector<float>> tubeTransformCol1Theta{treeReader, "tubeTransformLinearCol1Theta"};
    TTreeReaderValue<std::vector<float>> tubeTransformCol1Phi{treeReader, "tubeTransformLinearCol1Phi"};
 
    TTreeReaderValue<std::vector<float>> tubeTransformCol2Theta{treeReader, "tubeTransformLinearCol2Theta"};
    TTreeReaderValue<std::vector<float>> tubeTransformCol2Phi{treeReader, "tubeTransformLinearCol2Phi"};
 
    TTreeReaderValue<std::vector<float>> readOutPosX{treeReader, "readOutPosX"};
    TTreeReaderValue<std::vector<float>> readOutPosY{treeReader, "readOutPosY"};
    TTreeReaderValue<std::vector<float>> readOutPosZ{treeReader, "readOutPosZ"};
 
    while (treeReader.Next()) {
        MdtChamber newchamber{};
 
        newchamber.id.stationIndex = (*stationIndex);
        newchamber.id.eta = (*stationEta);
        newchamber.id.phi = (*stationPhi);
        newchamber.id.multilayer = (*stationML);
        newchamber.design = (*chamberDesign);
 
        newchamber.tubeRadius = (*tubeRadius);
        newchamber.tubePitch = (*tubePitch);
        
        newchamber.numTubes = (*numTubes);
        newchamber.numLayers = (*numLayers);
 
        Amg::RotationMatrix3D geoRot{Amg::RotationMatrix3D::Identity()};
        geoRot.col(0) = Amg::Vector3D((*geoModelTransformX)[1], (*geoModelTransformY)[1], (*geoModelTransformZ)[1]);
        geoRot.col(1) = Amg::Vector3D((*geoModelTransformX)[2], (*geoModelTransformY)[2], (*geoModelTransformZ)[2]);
        geoRot.col(2) = Amg::Vector3D((*geoModelTransformX)[3], (*geoModelTransformY)[3], (*geoModelTransformZ)[3]);       
        Amg::Vector3D geoTrans{(*geoModelTransformX)[0], (*geoModelTransformY)[0], (*geoModelTransformZ)[0]};
        newchamber.geoModelTransform = Amg::getTransformFromRotTransl(std::move(geoRot), std::move(geoTrans));       
 
        geoRot.col(0) = Amg::Vector3D((*alignableNodeX)[1], (*alignableNodeY)[1], (*alignableNodeZ)[1]);
        geoRot.col(1) = Amg::Vector3D((*alignableNodeX)[2], (*alignableNodeY)[2], (*alignableNodeZ)[2]);
        geoRot.col(2) = Amg::Vector3D((*alignableNodeX)[3], (*alignableNodeY)[3], (*alignableNodeZ)[3]);       
        geoTrans = Amg::Vector3D{(*alignableNodeX)[0], (*alignableNodeY)[0], (*alignableNodeZ)[0]};
        newchamber.alignableTransform = Amg::getTransformFromRotTransl(std::move(geoRot), std::move(geoTrans));
 
       
        for (size_t t = 0; t < tubeLayer->size(); ++t){
            using TubePositioning = MdtChamber::TubePositioning;
            TubePositioning newTube{};
            newTube.layerNum = (*tubeLayer)[t];
            newTube.tubeNum = (*tubeNumber)[t];
            newTube.tubeLength = (*tubeLength)[t];
            newTube.wireLength = (*wireLength)[t];
            Amg::RotationMatrix3D tubeRot{Amg::RotationMatrix3D::Identity()};
            tubeRot.col(0) = makeDir((*tubeTransformCol0Theta)[t], (*tubeTransformCol0Phi)[t]);
            tubeRot.col(1) = makeDir((*tubeTransformCol1Theta)[t], (*tubeTransformCol1Phi)[t]);
            tubeRot.col(2) = makeDir((*tubeTransformCol2Theta)[t], (*tubeTransformCol2Phi)[t]);
            Amg::Vector3D tubeTrans{(*tubeTransformTransX)[t],(*tubeTransformTransY)[t], (*tubeTransformTransZ)[t]};
            newTube.localToGlobal = Amg::getTransformFromRotTransl(std::move(tubeRot), std::move(tubeTrans));
            newTube.readoutPos = Amg::Vector3D{(*readOutPosX)[t],(*readOutPosY)[t],(*readOutPosZ)[t]};
            newchamber.insertTube(std::move(newTube));
        }
 
        auto insert_itr = to_ret.insert(std::move(newchamber));
        if (!insert_itr.second) {
            std::stringstream err{};
            err<<__FILE__<<":"<<__LINE__<<" The chamber "<<(*insert_itr.first).id
               <<" has already been inserted. "<<std::endl;
            throw std::runtime_error(err.str());
        }
    }
    std::cout<<"File parsing is finished. Found in total "<<to_ret.size()<<" readout element dumps "<<std::endl;
    return to_ret;
}

Variable Documentation

tolerance

double tolerance = 10 * Gaudi::Units::micrometer

constexpr

Definition at line 111 of file runMdtGeoComparison.cxx.