packages/trackbase/AlignmentTransformation.cc

0001 #include "AlignmentTransformation.h"
0002
0003 #include "ActsGeometry.h"
0004 #include "TpcDefs.h"
0005 #include "TrkrDefs.h"
0006
0007 #include <ffamodules/CDBInterface.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/PHDataNode.h>
0013 #include <phool/PHNode.h>
0014 #include <phool/PHNodeIterator.h>
0015 #include <phool/PHObject.h>
0016 #include <phool/PHTimer.h>
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>
0019
0020 #include <Acts/Surfaces/PerigeeSurface.hpp>
0021 #include <Acts/Surfaces/PlaneSurface.hpp>
0022 #include <Acts/Surfaces/Surface.hpp>
0023
0024 #include <Eigen/Dense>
0025 #include <Eigen/Geometry>
0026 #include <Eigen/LU>
0027
0028 #include <cmath>
0029 #include <fstream>
0030 #include <sstream>
0031
0032 void AlignmentTransformation::createMap(PHCompositeNode* topNode)
0033 {
0034   localVerbosity = 0;
0035   // The default is to use translation parameters that are in global coordinates
0036   std::cout << "AlignmentTransformation: use INTT survey geometry = " << use_intt_survey_geometry << std::endl;
0037   std::cout << "AlignmentTransformation: localVerbosity = " << localVerbosity << std::endl;
0038
0039   getNodes(topNode);
0040
0041   // Use construction transforms as a reference for making the map
0042   if (alignmentTransformationContainer::use_alignment)
0043   {
0044     alignmentTransformationContainer::use_alignment = false;
0045   }
0046
0047   // Define Parsing Variables
0048   TrkrDefs::hitsetkey hitsetkey = 0;
0049   float alpha = 0.0, beta = 0.0, gamma = 0.0, dx = 0.0, dy = 0.0, dz = 0.0, dgrx = 0.0, dgry = 0.0, dgrz = 0.0;
0050
0051   // load alignment constants file
0052   std::ifstream datafile;
0053   datafile.open(alignmentParamsFile);  //  looks for default file name on disk
0054   if (datafile.is_open())
0055   {
0056     std::cout << "AlignmentTransformation: Reading alignment parameters from disk file: "
0057               << alignmentParamsFile << " localVerbosity = " << localVerbosity << std::endl;
0058   }
0059   else
0060   {
0061     datafile.clear();
0062     // load alignment constants file from database
0063     alignmentParamsFile = CDBInterface::instance()->getUrl("TRACKINGALIGNMENT");
0064     std::cout << "AlignmentTransformation: Reading alignment parameters from database file: " << alignmentParamsFile << std::endl;
0065     datafile.open(alignmentParamsFile);
0066   }
0067
0068   ActsSurfaceMaps surfMaps = m_tGeometry->maps();
0069   Surface surf;
0070
0071   int linecount = 0;
0072   std::string str;
0073   while( std::getline(datafile, str) )
0074   {
0075     // trim leading space characters
0076     str.erase(str.begin(), std::find_if(str.begin(), str.end(), [](unsigned char ch) {return !std::isspace(ch);}));
0077
0078     // skip empty lines, or commented lines
0079     if( str.empty() ) { continue; }
0080     if( str.substr(0, 2) == "//" ) {continue;}
0081     if( str.substr(0, 1) == "#" ) {continue;}
0082
0083     // try read
0084     std::stringstream ss(str);
0085
0086   // check to see how many parameters per line in the file
0087   // If it is old, there may be only six. In that case, set the global rotation pars to zero, print a message.
0088     std::string dummy;
0089     int count = 0;
0090     while(ss >> dummy)
0091       {
0092     count ++;
0093       }
0094     if(count < 9)
0095       {
0096     std::stringstream str6(str);
0097     str6 >>  hitsetkey >> alpha >> beta >> gamma >> dx >> dy >> dz;
0098     if( str6.rdstate()&std::ios::failbit )
0099       {
0100         std::cout << "AlignmentTransformation::createMap - invalid line: " << str << " -------- Exiting" << std::endl;
0101         exit(1);
0102       }
0103     dgrx=0; dgry = 0; dgrz = 0;
0104
0105     if(linecount == 1 && localVerbosity > 0)
0106       {
0107         std::cout << PHWHERE << "The  alignment parameters file has only 6 parameters" << std::endl
0108               << "     --- setting global rotation parameters to zero!" << std::endl;
0109       }
0110       }
0111     else
0112       {
0113     std::stringstream str9(str);
0114     str9 >> hitsetkey >> alpha >> beta >> gamma >> dx >> dy >> dz >> dgrx >> dgry >> dgrz;
0115     if( str9.rdstate()&std::ios::failbit )
0116       {
0117         std::cout << "AlignmentTransformation::createMap - invalid line: " << str << " -------- Exiting" << std::endl;
0118         exit(1);
0119       }
0120       }
0121
0122     linecount ++;
0123
0124     if(localVerbosity > 0)
0125       {
0126     std::cout  <<  hitsetkey << "  " << alpha  << "  " << beta  << "  " << gamma  << " dx " << dx  << " dy " << dy << " dz "  << dz
0127            << " dgrx " << dgrx << " dgry " << dgry << " dgrz " << dgrz << std::endl;
0128       }
0129
0130     // Perturbation translations and angles for stave and sensor
0131     Eigen::Vector3d sensorAngles(alpha, beta, gamma);
0132     Eigen::Vector3d millepedeTranslation(dx, dy, dz);
0133     Eigen::Vector3d sensorAnglesGlobal(dgrx, dgry, dgrz);
0134
0135     unsigned int trkrId = TrkrDefs::getTrkrId(hitsetkey);  // specify between detectors
0136
0137     perturbationAngles = Eigen::Vector3d(0.0, 0.0, 0.0);
0138     perturbationAnglesGlobal = Eigen::Vector3d(0.0, 0.0, 0.0);
0139     perturbationTranslation = Eigen::Vector3d(0.0, 0.0, 0.0);
0140
0141     switch( trkrId )
0142     {
0143
0144       case TrkrDefs::mvtxId:
0145       {
0146         if (perturbMVTX)
0147         {
0148           generateRandomPerturbations(mvtxAngleDev, mvtxTransDev);
0149           sensorAngles = sensorAngles + perturbationAngles;
0150           millepedeTranslation = millepedeTranslation + perturbationTranslation;
0151         }
0152
0153         surf = surfMaps.getSiliconSurface(hitsetkey);
0154
0155         Acts::Transform3 transform;
0156         transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, sensorAnglesGlobal, false);
0157
0158         Acts::GeometryIdentifier id = surf->geometryId();
0159
0160         if (localVerbosity)
0161         {
0162           std::cout << " Add transform for MVTX with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl;
0163           std::cout << " final mvtx transform:" << std::endl
0164             << transform.matrix() << std::endl;
0165         }
0166         transformMap->addTransform(id, transform);
0167         transformMapTransient->addTransform(id, transform);
0168
0169         break;
0170       }
0171
0172       case TrkrDefs::inttId:
0173       {
0174         if (perturbINTT)
0175         {
0176           generateRandomPerturbations(inttAngleDev, inttTransDev);
0177           sensorAngles = sensorAngles + perturbationAngles;
0178           millepedeTranslation = millepedeTranslation + perturbationTranslation;
0179         }
0180
0181         surf = surfMaps.getSiliconSurface(hitsetkey);
0182
0183         Acts::Transform3 transform;
0184         transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, sensorAnglesGlobal, use_intt_survey_geometry);
0185         Acts::GeometryIdentifier id = surf->geometryId();
0186
0187         if (localVerbosity)
0188         {
0189           std::cout << " Add transform for INTT with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl;
0190         }
0191
0192         transformMap->addTransform(id, transform);
0193         transformMapTransient->addTransform(id, transform);
0194         break;
0195       }
0196
0197       case TrkrDefs::tpcId:
0198       {
0199         if (perturbTPC)
0200         {
0201           generateRandomPerturbations(tpcAngleDev, tpcTransDev);
0202           sensorAngles = sensorAngles + perturbationAngles;
0203           millepedeTranslation = millepedeTranslation + perturbationTranslation;
0204         }
0205
0206         unsigned int side = TpcDefs::getSide(hitsetkey);
0207         unsigned int sector = TpcDefs::getSectorId(hitsetkey);
0208         int subsurfkey_min = (1 - side) * 144 + (144 - sector * 12) - 12 - 6;
0209         int subsurfkey_max = subsurfkey_min + 12;
0210         // std::cout << " sector " << sector << " side " << side << " subsurfkey_min " << subsurfkey_min << " subsurfkey_max " << subsurfkey_max << std::endl;
0211
0212         for (int subsurfkey = subsurfkey_min; subsurfkey < subsurfkey_max; subsurfkey++)
0213         {
0214           int sskey = subsurfkey;
0215           if (sskey < 0)
0216           {
0217             sskey += 288;
0218           }
0219
0220           surf = surfMaps.getTpcSurface(hitsetkey, (unsigned int) sskey);
0221
0222           Acts::Transform3 transform;
0223           transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, sensorAnglesGlobal, false);
0224           Acts::GeometryIdentifier id = surf->geometryId();
0225
0226           if (localVerbosity)
0227           {
0228             unsigned int layer = TrkrDefs::getLayer(hitsetkey);
0229             std::cout << " Add transform for TPC with surface GeometryIdentifier " << id << std::endl
0230               << " trkrid " << trkrId << " hitsetkey " << hitsetkey << " layer " << layer << " sector " << sector << " side " << side
0231               << " subsurfkey " << subsurfkey << std::endl;
0232             Acts::Vector3 center = surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1;  // convert to cm
0233             std::cout << "Ideal surface center: " << std::endl
0234               << center << std::endl;
0235             std::cout << "transform matrix: " << std::endl
0236               << transform.matrix() << std::endl;
0237           }
0238           transformMap->addTransform(id, transform);
0239           transformMapTransient->addTransform(id, transform);
0240         }
0241
0242         break;
0243       }
0244
0245       case TrkrDefs::micromegasId:
0246       {
0247         if (perturbMM)
0248         {
0249           generateRandomPerturbations(mmAngleDev, mmTransDev);
0250
0251           sensorAngles = sensorAngles + perturbationAngles;
0252           millepedeTranslation = millepedeTranslation + perturbationTranslation;
0253         }
0254         surf = surfMaps.getMMSurface(hitsetkey);
0255
0256         Acts::Transform3 transform;
0257         transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, sensorAnglesGlobal, false);
0258         Acts::GeometryIdentifier id = surf->geometryId();
0259
0260         if (localVerbosity)
0261         {
0262           std::cout << " Add transform for Micromegas with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl;
0263         }
0264
0265         transformMap->addTransform(id, transform);
0266         transformMapTransient->addTransform(id, transform);
0267         break;
0268       }
0269
0270       default:
0271       {
0272         std::cout << "AlignmentTransformation::createMap - Invalid Hitsetkey: " << hitsetkey << std::endl;
0273         break;
0274       }
0275
0276     }
0277
0278   }
0279
0280   // copy map into geoContext
0281   m_tGeometry->geometry().geoContext = transformMap;
0282
0283   std::cout << " AlignmentTransformation processed " << linecount << " input lines " << std::endl;
0284
0285   // map is created, now we can use the transforms
0286   alignmentTransformationContainer::use_alignment = true;
0287 }
0288
0289 // currently used as the transform maker
0290 Acts::Transform3 AlignmentTransformation::newMakeTransform(const Surface& surf, Eigen::Vector3d& millepedeTranslation, Eigen::Vector3d& sensorAngles, Eigen::Vector3d& sensorAnglesGlobal, bool survey)
0291 {
0292   // define null matrices
0293   Eigen::Vector3d nullTranslation(0, 0, 0);
0294   Eigen::AngleAxisd a(0, Eigen::Vector3d::UnitX());
0295   Eigen::AngleAxisd b(0, Eigen::Vector3d::UnitY());
0296   Eigen::AngleAxisd g(0, Eigen::Vector3d::UnitZ());
0297   Eigen::Quaternion<double> qnull = g * b * a;
0298   Eigen::Matrix3d nullRotation = qnull.matrix();
0299
0300   // get the acts transform components
0301   // Note that Acts transforms local coordinates of (x,z,y) to global (x,y,z)
0302   Acts::Transform3 actsTransform = surf->transform(m_tGeometry->geometry().getGeoContext());
0303   Eigen::Matrix3d actsRotationPart = actsTransform.rotation();
0304   Eigen::Vector3d actsTranslationPart = actsTransform.translation();
0305
0306   // Create  alignment local coordinates rotation matrix
0307   Eigen::AngleAxisd alpha(sensorAngles(0), Eigen::Vector3d::UnitX());
0308   Eigen::AngleAxisd beta(sensorAngles(1), Eigen::Vector3d::UnitY());
0309   Eigen::AngleAxisd gamma(sensorAngles(2), Eigen::Vector3d::UnitZ());
0310   Eigen::Quaternion<double> q = gamma * beta * alpha;
0311   Eigen::Matrix3d millepedeRotation = q.matrix();
0312
0313  // Create alignment global coordinates rotation matrix
0314   Eigen::AngleAxisd grx(sensorAnglesGlobal(0), Eigen::Vector3d::UnitX());
0315   Eigen::AngleAxisd gry(sensorAnglesGlobal(1), Eigen::Vector3d::UnitY());
0316   Eigen::AngleAxisd grz(sensorAnglesGlobal(2), Eigen::Vector3d::UnitZ());
0317   Eigen::Quaternion<double> gqr = grz * gry * grx;
0318   Eigen::Matrix3d millepedeRotationGlobal = gqr.matrix();
0319
0320   // and make affine matrices from each
0321
0322   Acts::Transform3 mpRotationAffine;
0323   mpRotationAffine.linear() = millepedeRotation;
0324   mpRotationAffine.translation() = nullTranslation;
0325
0326   Acts::Transform3 mpRotationGlobalAffine;
0327   mpRotationGlobalAffine.linear() = millepedeRotationGlobal;
0328   mpRotationGlobalAffine.translation() = nullTranslation;
0329
0330   Acts::Transform3 mpTranslationAffine;
0331   mpTranslationAffine.linear() = nullRotation;
0332   mpTranslationAffine.translation() = millepedeTranslation;
0333
0334   Acts::Transform3 actsRotationAffine;
0335   actsRotationAffine.linear() = actsRotationPart;
0336   actsRotationAffine.translation() = nullTranslation;
0337   Acts::Transform3 actsTranslationAffine;
0338   actsTranslationAffine.linear() = nullRotation;
0339   actsTranslationAffine.translation() = actsTranslationPart;
0340
0341   // Put them together into a combined transform
0342   Acts::Transform3 transform;
0343   //! If we read the survey parameters directly, that is the full transform
0344   if (survey)
0345   {
0346     //! The millepede affines will just be what was read in, which was the
0347     //! survey information. This should (in principle) be equivalent to
0348     //! the ideal position + any misalignment
0349     transform = mpTranslationAffine  * mpRotationGlobalAffine * mpRotationAffine;
0350   }
0351   else
0352   {
0353     transform = mpTranslationAffine * mpRotationGlobalAffine * actsTranslationAffine * mpRotationAffine * actsRotationAffine;
0354   }
0355
0356   if (localVerbosity)
0357   {
0358     Acts::Transform3 actstransform = actsTranslationAffine * actsRotationAffine;
0359
0360     std::cout << "newMakeTransform" << std::endl;
0361     std::cout << "Input translation: " << std::endl << millepedeTranslation << std::endl;
0362     std::cout << "Input sensorAngles: " << std::endl << sensorAngles << std::endl;
0363     std::cout << "Input sensorAnglesGlobal: " << std::endl << sensorAnglesGlobal << std::endl;
0364     std::cout << "mpRotationAffine: " << std::endl
0365               << mpRotationAffine.matrix() << std::endl;
0366       std::cout << "millepederotation * acts " << std::endl
0367               << millepedeRotation * actsRotationPart << std::endl;
0368     std::cout << "actsRotationAffine: " << std::endl
0369               << actsRotationAffine.matrix() << std::endl;
0370     std::cout << "actsTranslationAffine: " << std::endl
0371               << actsTranslationAffine.matrix() << std::endl;
0372     std::cout << "full acts transform " << std::endl
0373               << actstransform.matrix() << std::endl;
0374     std::cout << "mpRotationGlobalAffine: " << std::endl
0375           << mpRotationGlobalAffine.matrix() << std::endl;
0376     std::cout << "mpTranslationAffine: " << std::endl
0377           << mpTranslationAffine.matrix() << std::endl;
0378     std::cout << "Overall transform: " << std::endl
0379               << transform.matrix() << std::endl;
0380     std::cout << "overall * idealinv " << std::endl
0381               << (transform * actstransform.inverse()).matrix() << std::endl;
0382     std::cout << "overall - ideal " << std::endl;
0383     for (int test = 0; test < transform.matrix().rows(); test++)
0384     {
0385       for (int test2 = 0; test2 < transform.matrix().cols(); test2++)
0386       {
0387         std::cout << transform(test, test2) - actstransform(test, test2) << ", ";
0388       }
0389       std::cout << std::endl;
0390     }
0391   }
0392
0393   return transform;
0394 }
0395
0396 int AlignmentTransformation::getNodes(PHCompositeNode* topNode)
0397 {
0398   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0399   if (!m_tGeometry)
0400   {
0401     std::cout << "ActsGeometry not on node tree. Exiting."
0402               << std::endl;
0403
0404     return Fun4AllReturnCodes::ABORTEVENT;
0405   }
0406
0407   return 0;
0408 }
0409
0410 void AlignmentTransformation::misalignmentFactor(uint8_t layer, const double factor)
0411 {
0412   transformMap->setMisalignmentFactor(layer, factor);
0413   transformMapTransient->setMisalignmentFactor(layer, factor);
0414 }
0415 void AlignmentTransformation::createAlignmentTransformContainer(PHCompositeNode* topNode)
0416 {
0417   // Get a pointer to the top of the node tree
0418   PHNodeIterator iter(topNode);
0419
0420   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0421   if (!dstNode)
0422   {
0423     std::cerr << "DST node is missing, quitting" << std::endl;
0424     throw std::runtime_error("Failed to find DST node in AlignmentTransformation::createNodes");
0425   }
0426
0427   transformMap = findNode::getClass<alignmentTransformationContainer>(topNode, "alignmentTransformationContainer");
0428   if (!transformMap)
0429   {
0430     transformMap = new alignmentTransformationContainer;
0431     auto node = new PHDataNode<alignmentTransformationContainer>(transformMap, "alignmentTransformationContainer");
0432     dstNode->addNode(node);
0433   }
0434
0435   transformMapTransient = findNode::getClass<alignmentTransformationContainer>(topNode, "alignmentTransformationContainerTransient");
0436   if (!transformMapTransient)
0437   {
0438     transformMapTransient = new alignmentTransformationContainer;
0439     auto node = new PHDataNode<alignmentTransformationContainer>(transformMapTransient, "alignmentTransformationContainerTransient");
0440     dstNode->addNode(node);
0441   }
0442 }
0443
0444 void AlignmentTransformation::generateRandomPerturbations(Eigen::Vector3d angleDev, Eigen::Vector3d transformDev)
0445 {
0446   /*Creates random perturbations for the correctional parameters with a given standard deviation and mean of zero*/
0447
0448   std::cout << "Generating Random Perturbations..." << std::endl;
0449
0450   if (angleDev(0) != 0)
0451   {
0452     std::normal_distribution<double> distribution(0, angleDev(0));
0453     perturbationAngles(0) = distribution(generator);
0454   }
0455   if (angleDev(1) != 0)
0456   {
0457     std::normal_distribution<double> distribution(0, angleDev(1));
0458     perturbationAngles(1) = distribution(generator);
0459   }
0460   if (angleDev(2) != 0)
0461   {
0462     std::normal_distribution<double> distribution(0, angleDev(2));
0463     perturbationAngles(2) = distribution(generator);
0464   }
0465   if (transformDev(0) != 0)
0466   {
0467     std::normal_distribution<double> distribution(0, transformDev(0));
0468     perturbationTranslation(0) = distribution(generator);
0469   }
0470   if (transformDev(1) != 0)
0471   {
0472     std::normal_distribution<double> distribution(0, transformDev(1));
0473     perturbationTranslation(1) = distribution(generator);
0474   }
0475   if (transformDev(2) != 0)
0476   {
0477     std::normal_distribution<double> distribution(0, transformDev(2));
0478     perturbationTranslation(2) = distribution(generator);
0479   }
0480   if (localVerbosity)
0481   {
0482     std::cout << "randomperturbationAngles" << perturbationAngles << " randomperturbationTrans:" << perturbationTranslation << std::endl;
0483   }
0484 }