runMmGeoComparison.cxx File Reference

#include <GeoPrimitives/GeoPrimitives.h>
#include <GeoPrimitives/GeoPrimitivesHelpers.h>
#include <GeoPrimitives/GeoPrimitivesToStringConverter.h>
#include <GaudiKernel/SystemOfUnits.h>
#include <MuonReadoutGeometryR4/MuonDetectorDefs.h>
#include <string>
#include <set>
#include <vector>
#include <map>
#include <iostream>
#include <cmath>
#include <PathResolver/PathResolver.h>
#include <TFile.h>
#include <TTreeReader.h>

Go to the source code of this file.

Classes
struct	MmChamber
	Helper struct to represent a full MicroMegas chamber. More...
struct	MmChamber::MmChannel
struct	MmChamber::MmLayer
	Helper struct to assess that the layers are properly oriented. More...

Macros
#define	TEST_BASICPROP(attribute, propName)

Functions
std::ostream &	operator<< (std::ostream &ostr, const MmChamber &chamb)
	Translation of the station Index -> station Name. More...
std::ostream &	operator<< (std::ostream &ostr, const MmChamber::MmChannel &channel)
std::ostream &	operator<< (std::ostream &ostr, const MmChamber::MmLayer &layer)
std::set< MmChamber >	readTreeDump (const std::string &inputFile)
int	main (int argc, char **argv)

Variables
constexpr double	tolerance = 0.003*Gaudi::Units::millimeter

Macro Definition Documentation

TEST_BASICPROP

#define TEST_BASICPROP	(		attribute,
			propName
	)

Value:

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

Definition at line 278 of file runMmGeoComparison.cxx.

Function Documentation

main()

int main	(	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

GasGap Dimensions for debugging

Eta strips have their centres at local Y = 0 --> Their positions must match in absolute terms w.r.t legacy geometry

The centres of the Stereo layers are defined as the bisector of the line between the two frame edges. However, thus far the parameter book deviates from the legacy Run-3 implementation. --> Cannot compare the absolute position of the stereo layers. Instead check that the left edge, right edge and center point in the new geometry are on the same line as defined by the reference system.

Definition at line 285 of file runMmGeoComparison.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<MmChamber> refChambers = readTreeDump(refFile);
    if (refChambers.empty()) {
        std::cerr<<"The file "<<refFile<<" should contain at least one chamber "<<std::endl;
        return EXIT_FAILURE;
    }
    std::set<MmChamber> 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;
    for (const MmChamber& reference : refChambers) {
        std::set<MmChamber>::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 MmChamber& test = {*test_itr};
 
        TEST_BASICPROP(ActiveWidthS, "GasGap length on the short side");
        TEST_BASICPROP(ActiveWidthL, "GasGap length on the long side");
        TEST_BASICPROP(ActiveHeightR, "GasGap Height");
        TEST_BASICPROP(stripPitch, "Strip pitch");
        // if (!chamberOkay) continue;
        using MmLayer = MmChamber::MmLayer;
        for (const MmLayer& refLayer : reference.layers) {
            break;
            std::set<MmLayer>::const_iterator lay_itr = test.layers.find(refLayer);
            if (lay_itr == test.layers.end()) {
                std::cerr<<"runMmGeoComparison() "<<__LINE__<<": in "<<test<<" "
                         <<refLayer<<" is not found. "<<std::endl;
                chamberOkay = false;
                continue;
            }
            const MmLayer& testLayer{*lay_itr};
            if ( (refLayer.firstStripPos- testLayer.firstStripPos).mag() > tolerance) {
                std::cerr<<"runMmGeoComparison() "<<__LINE__<<": in "<<test<<" "
                         <<testLayer.gasGap<<" has different starting position "
                         <<Amg::toString(refLayer.firstStripPos, 2) <<" vs. "
                         <<Amg::toString(testLayer.firstStripPos, 2)
                         <<"difference: "<<Amg::toString(refLayer.firstStripPos- testLayer.firstStripPos, 2)
                         <<" / "<<(refLayer.firstStripPos- testLayer.firstStripPos).mag()/reference.stripPitch
                         <<std::endl;
                chamberOkay = false;
            }
            if (refLayer.firstStrip != testLayer.firstStrip) {
                 std::cerr<<"runMmGeoComparison() "<<__LINE__<<": in "<<test<<" "
                         <<testLayer.gasGap<<" starts from different strip "<<refLayer.firstStrip<<" vs. "
                         <<testLayer.firstStrip<<std::endl;
                chamberOkay = false;
            }
            
            
            if (!Amg::doesNotDeform(refLayer.transform.inverse()* testLayer.transform)){
                std::cerr<<"runMmGeoComparison() "<<__LINE__<<": in "<<test<<" "
                         <<testLayer<<" differs w.r.t. reference "<<Amg::toString(refLayer.transform)<<". vs. "
                         <<Amg::toString(refLayer.transform.inverse()*testLayer.transform)<<std::endl;
                chamberOkay = false;
            }           
        }
        unsigned int failedEta{0}, failedStereo{0};
        for (const MmChamber::MmChannel& refStrip : reference.channels) {
            std::set<MmChamber::MmChannel>::const_iterator strip_itr = test.channels.find(refStrip);
            if (strip_itr == test.channels.end()) {
                std::cerr<<"runMmGeoComparison() "<<__LINE__<<": in "<<test<<" "
                          <<refStrip<<" is not found. "<<std::endl;
                chamberOkay = false;
                continue;
            }
            const MmChamber::MmChannel& testStrip{*strip_itr};
        
            if (!refStrip.isStereo && (failedEta <=10)) {
                const Amg::Vector3D diffStrip{testStrip.globCenter - refStrip.globCenter};
                if (diffStrip.mag() > tolerance) {
                    std::cerr<<"runMmGeoComparison() "<<__LINE__<<": In "
                             <<test<<" " <<testStrip <<"/local: "<<Amg::toString(testStrip.locCenter, 2) 
                             <<" should be located at "<<Amg::toString(refStrip.globCenter, 2)
                             <<"/local: "<<Amg::toString(refStrip.locCenter, 2)
                             <<" displacement: "<<Amg::toString(diffStrip, 2)<<"  "<<diffStrip.mag()<<std::endl;
                    chamberOkay = false;
                }
                ++failedEta;
            } 
            else if (failedStereo <= 10) {
                const Amg::Vector3D stripDir{Amg::getRotateZ3D(90*Gaudi::Units::deg)*
                                             (refStrip.leftEdge - refStrip.rightEdge).unit()};
                const Amg::Vector3D testDir{Amg::getRotateZ3D(90*Gaudi::Units::deg)*
                                            (testStrip.rightEdge - testStrip.leftEdge).unit()};
                const double centerDist = stripDir.dot(testStrip.globCenter - refStrip.globCenter);
                const double leftDist = stripDir.dot(testStrip.leftEdge -refStrip.globCenter);
                const double rightDist = stripDir.dot(testStrip.rightEdge - refStrip.globCenter);
                if ( std::abs(centerDist) > tolerance || std::abs(leftDist) > tolerance || std::abs(rightDist) > tolerance) {
                     std::cerr<<"runMmGeoComparison() "<<__LINE__<<": In "
                             <<test<<" " <<testStrip <<" + mu "<<Amg::toString(testDir,2) 
                             <<"/local: "<<Amg::toString(testStrip.locCenter, 2) 
                             <<" does not describe the same stereo strip as "
                             <<Amg::toString(refStrip.globCenter, 2)<<" + lambda "<<Amg::toString(stripDir,2)
                             <<". Distances to the left-edge/center/right-edge: "
                             <<leftDist<<"/"<<centerDist<<"/"<<rightDist<<", dot: "
                             <<std::acos(std::clamp(stripDir.dot(testDir),- 1., 1.)) / Gaudi::Units::deg<<std::endl;
                    chamberOkay = false;
                }
                ++failedStereo;
            }
        }
 
        if (!chamberOkay) {
            return_code = EXIT_FAILURE;
        }
    }
    return return_code;
 
}

operator<<() [1/3]

std::ostream& operator<<	(	std::ostream &	ostr,
		const MmChamber &	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 114 of file runMmGeoComparison.cxx.

                                                                 {
    static const std::map<int, std::string> stationDict{
        {55, "MMS"}, {56, "MML"}
    };
    ostr<<"MicroMegas chamber "<<stationDict.at(chamb.stationIndex)<<std::abs(chamb.stationEta)<<
        (chamb.stationEta>0 ? "A" : "C")<<chamb.stationPhi<<"-"<<chamb.multilayer<<" ";
    return ostr;
}

operator<<() [2/3]

std::ostream& operator<<	(	std::ostream &	ostr,
		const MmChamber::MmChannel &	channel
	)

Definition at line 123 of file runMmGeoComparison.cxx.

                                                                            {
    ostr<<"channel (gasGap/number): ";
    ostr<<channel.gasGap<<"/";
    ostr<<std::setfill('0')<<std::setw(4)<<channel.channel<<", ";
    ostr<<"center: "<<Amg::toString(channel.globCenter, 2);
    return ostr;
}

operator<<() [3/3]

std::ostream& operator<<	(	std::ostream &	ostr,
		const MmChamber::MmLayer &	layer
	)

Definition at line 132 of file runMmGeoComparison.cxx.

                                                                        {
    ostr<<"Mmlayer (gasGap): ";
    ostr<<layer.gasGap<<", ";
    ostr<<"transform: "<<Amg::toString(layer.transform);
    return ostr;
}

readTreeDump()

std::set<MmChamber> readTreeDump

(

const std::string &

inputFile

)

Identifier of the readout element

Strip Length and Positions

GasGap Dimensions for debugging

Geo Model transformation

Identifier of the readout element

Gas Gap lengths for debug

Definition at line 139 of file runMmGeoComparison.cxx.

                                                           {
    std::set<MmChamber> to_ret{};
    std::cout<<"Read the MicroMegas 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("MmGeoModelTree", inFile.get());
    if (treeReader.IsInvalid()){
        std::cerr<<__FILE__<<":"<<__LINE__<<" The file "<<inputFile<<" does not contain the 'MmGeoModelTree'"<<std::endl;
        return to_ret;
    }
    
    TTreeReaderValue<unsigned short> stationIndex{treeReader, "stationIndex"};
    TTreeReaderValue<short> stationEta{treeReader, "stationEta"};
    TTreeReaderValue<short> stationPhi{treeReader, "stationPhi"};
    TTreeReaderValue<short> multilayer{treeReader, "multilayer"};
 
    TTreeReaderValue<std::vector<uint>> channel{treeReader, "channel"};
    TTreeReaderValue<std::vector<short>> gasGap{treeReader, "gasGap"};
    TTreeReaderValue<std::vector<float>> stripLength{treeReader, "stripLength"};
 
    TTreeReaderValue<std::vector<bool>> isStereo{treeReader, "isStereo"};
    TTreeReaderValue<std::vector<float>> locStripCenterX{treeReader, "locStripCenterX"};
    TTreeReaderValue<std::vector<float>> locStripCenterY{treeReader, "locStripCenterY"};
 
    TTreeReaderValue<std::vector<float>> stripCenterX{treeReader, "stripCenterX"};
    TTreeReaderValue<std::vector<float>> stripCenterY{treeReader, "stripCenterY"};
    TTreeReaderValue<std::vector<float>> stripCenterZ{treeReader, "stripCenterZ"};
    TTreeReaderValue<std::vector<float>> stripLeftEdgeX{treeReader, "stripLeftEdgeX"};
    TTreeReaderValue<std::vector<float>> stripLeftEdgeY{treeReader, "stripLeftEdgeY"};
    TTreeReaderValue<std::vector<float>> stripLeftEdgeZ{treeReader, "stripLeftEdgeZ"};
    TTreeReaderValue<std::vector<float>> stripRightEdgeX{treeReader, "stripRightEdgeX"};
    TTreeReaderValue<std::vector<float>> stripRightEdgeY{treeReader, "stripRightEdgeY"};
    TTreeReaderValue<std::vector<float>> stripRightEdgeZ{treeReader, "stripRightEdgeZ"};
 
    TTreeReaderValue<float> ActiveHeightR{treeReader, "ActiveHeightR"};
    TTreeReaderValue<float> ActiveWidthS{treeReader, "ActiveWidthS"};
    TTreeReaderValue<float> ActiveWidthL{treeReader, "ActiveWidthL"};
    TTreeReaderValue<float> stripPitch{treeReader, "stripPitch"};
 
    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>> stripRotCol1X{treeReader, "stripRotLinearCol1X"};
    TTreeReaderValue<std::vector<float>> stripRotCol1Y{treeReader, "stripRotLinearCol1Y"};
    TTreeReaderValue<std::vector<float>> stripRotCol1Z{treeReader, "stripRotLinearCol1Z"};
 
    TTreeReaderValue<std::vector<float>> stripRotCol2X{treeReader, "stripRotLinearCol2X"};
    TTreeReaderValue<std::vector<float>> stripRotCol2Y{treeReader, "stripRotLinearCol2Y"};
    TTreeReaderValue<std::vector<float>> stripRotCol2Z{treeReader, "stripRotLinearCol2Z"};
 
    TTreeReaderValue<std::vector<float>> stripRotCol3X{treeReader, "stripRotLinearCol3X"};
    TTreeReaderValue<std::vector<float>> stripRotCol3Y{treeReader, "stripRotLinearCol3Y"};
    TTreeReaderValue<std::vector<float>> stripRotCol3Z{treeReader, "stripRotLinearCol3Z"};
 
    TTreeReaderValue<std::vector<float>> stripRotTransX{treeReader, "stripRotTranslationX"};
    TTreeReaderValue<std::vector<float>> stripRotTransY{treeReader, "stripRotTranslationY"};
    TTreeReaderValue<std::vector<float>> stripRotTransZ{treeReader, "stripRotTranslationZ"};
 
    TTreeReaderValue<std::vector<uint8_t>> stripRotGasGap{treeReader, "stripRotGasGap"};
 
    TTreeReaderValue<std::vector<float>> firstStripPosX{treeReader, "firstStripPosX"};
    TTreeReaderValue<std::vector<float>> firstStripPosY{treeReader, "firstStripPosY"};
 
    TTreeReaderValue<std::vector<int>> readoutSide{treeReader, "stripReadoutSide"};
    TTreeReaderValue<std::vector<unsigned int>> firstStripNum{treeReader, "stripFirstStrip"};
 
 
    while (treeReader.Next()) {
        MmChamber newchamber{};
 
        newchamber.stationIndex = (*stationIndex);
        newchamber.stationEta = (*stationEta);
        newchamber.stationPhi = (*stationPhi);
        newchamber.multilayer = (*multilayer);
 
        newchamber.ActiveHeightR = (*ActiveHeightR);
        newchamber.ActiveWidthS = (*ActiveWidthS);
        newchamber.ActiveWidthL = (*ActiveWidthL);
        newchamber.stripPitch = (*stripPitch);
 
        Amg::Vector3D geoTrans{(*geoModelTransformX)[0], (*geoModelTransformY)[0], (*geoModelTransformZ)[0]};
        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]);       
        newchamber.geoModelTransform = Amg::getTransformFromRotTransl(std::move(geoRot), std::move(geoTrans)); 
        
        //Strips
        for (size_t s = 0; s < stripCenterX->size(); ++s){
            MmChamber::MmChannel newStrip{};
            newStrip.globCenter = Amg::Vector3D{(*stripCenterX)[s], (*stripCenterY)[s], (*stripCenterZ)[s]}; 
            newStrip.locCenter = Amg::Vector2D{(*locStripCenterX)[s], (*locStripCenterY)[s]};
            newStrip.leftEdge = Amg::Vector3D{(*stripLeftEdgeX)[s], (*stripLeftEdgeY)[s], (*stripLeftEdgeZ)[s]};
            newStrip.rightEdge = Amg::Vector3D{(*stripRightEdgeX)[s], (*stripRightEdgeY)[s], (*stripRightEdgeZ)[s]};
            
            newStrip.gasGap = (*gasGap)[s];
            newStrip.channel = (*channel)[s];
            newStrip.isStereo = (*isStereo)[s];
            newchamber.channels.insert(std::move(newStrip));
        }
 
        for (size_t l = 0; l < stripRotGasGap->size(); ++l){
            MmChamber::MmLayer newLayer{};
            newLayer.gasGap = (*stripRotGasGap)[l];
            Amg::RotationMatrix3D stripRot{Amg::RotationMatrix3D::Identity()};
            stripRot.col(0) = Amg::Vector3D((*stripRotCol1X)[l],(*stripRotCol1Y)[l], (*stripRotCol1Z)[l]);
            stripRot.col(1) = Amg::Vector3D((*stripRotCol2X)[l],(*stripRotCol2Y)[l], (*stripRotCol2Z)[l]);
            stripRot.col(2) = Amg::Vector3D((*stripRotCol3X)[l],(*stripRotCol3Y)[l], (*stripRotCol3Z)[l]);
            Amg::Vector3D layTrans{(*stripRotTransX)[l], (*stripRotTransY)[l], (*stripRotTransZ)[l]};
            newLayer.transform = Amg::getTransformFromRotTransl(std::move(stripRot), std::move(layTrans));
            newLayer.firstStripPos = Amg::Vector2D{(*firstStripPosX)[l], (*firstStripPosY)[l]};
            newLayer.readoutSide = (*readoutSide)[l];
            newLayer.firstStrip = (*firstStripNum)[l];
            newchamber.layers.insert(std::move(newLayer));
        }
        
        auto insert_itr = to_ret.insert(std::move(newchamber));
        if (!insert_itr.second) {
            std::stringstream err{};
            err<<__FILE__<<":"<<__LINE__<<" The chamber "<<(*insert_itr.first).stationIndex
               <<" 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

constexpr double tolerance = 0.003*Gaudi::Units::millimeter

constexpr

Definition at line 30 of file runMmGeoComparison.cxx.