algorithms/fardetectors/MatrixTransferStatic.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022, 2023 Alex Jentsch, Wouter Deconinck, Sylvester Joosten, David Lawrence, Simon Gardner
0003 //
0004 // This converted from: https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/FarForwardParticles.cpp
0005
0006 #include "MatrixTransferStatic.h"
0007
0008 #include <DD4hep/Alignments.h>
0009 #include <DD4hep/DetElement.h>
0010 #include <DD4hep/Objects.h>
0011 #include <DD4hep/VolumeManager.h>
0012 #include <Evaluator/DD4hepUnits.h>
0013 #include <Math/GenVector/Cartesian3D.h>
0014 #include <Math/GenVector/DisplacementVector3D.h>
0015 #include <edm4hep/Vector3f.h>
0016 #include <edm4hep/utils/vector_utils.h>
0017 #include <cmath>
0018 #include <gsl/pointers>
0019 #include <vector>
0020
0021 #include "algorithms/fardetectors/MatrixTransferStaticConfig.h"
0022
0023 void eicrecon::MatrixTransferStatic::init() {
0024
0025 }
0026
0027 void eicrecon::MatrixTransferStatic::process(
0028     const MatrixTransferStatic::Input& input,
0029     const MatrixTransferStatic::Output& output) const {
0030
0031   const auto [mcparts, rechits] = input;
0032   auto [outputParticles] = output;
0033
0034   std::vector<std::vector<double>> aX(m_cfg.aX);
0035   std::vector<std::vector<double>> aY(m_cfg.aY);
0036
0037   //----- Define constants here ------
0038   double aXinv[2][2] = {{0.0, 0.0},
0039                         {0.0, 0.0}};
0040   double aYinv[2][2] = {{0.0, 0.0},
0041                         {0.0, 0.0}};
0042
0043   double nomMomentum     = m_cfg.nomMomentum; //extract the nominal value first -- will be overwritten by MCParticle
0044   double local_x_offset  = m_cfg.local_x_offset;
0045   double local_y_offset  = m_cfg.local_y_offset;
0046   double local_x_slope_offset  = m_cfg.local_x_slope_offset;
0047   double local_y_slope_offset  = m_cfg.local_y_slope_offset;
0048
0049   double numBeamProtons = 0;
0050   double runningMomentum = 0.0;
0051
0052   for (const auto& p: *mcparts) {
0053           if(mcparts->size() == 1 && p.getPDG() == 2212){
0054                 runningMomentum = p.getMomentum().z;
0055                 numBeamProtons++;
0056           }
0057         if (p.getGeneratorStatus() == 4 && p.getPDG() == 2212) { //look for "beam" proton
0058                 runningMomentum += p.getMomentum().z;
0059                 numBeamProtons++;
0060         }
0061   }
0062
0063   if(numBeamProtons == 0) {error("No beam protons to choose matrix!! Skipping!!"); return;}
0064
0065   nomMomentum = runningMomentum/numBeamProtons;
0066
0067   double nomMomentumError = 0.05;
0068
0069   //This is a temporary solution to get the beam energy information
0070   //needed to select the correct matrix
0071
0072   if(abs(275.0 - nomMomentum)/275.0 < nomMomentumError){
0073
0074      aX[0][0] = 3.251116; //a
0075      aX[0][1] = 30.285734; //b
0076      aX[1][0] = 0.186036375; //c
0077      aX[1][1] = 0.196439472; //d
0078
0079      aY[0][0] = 0.4730500000; //a
0080      aY[0][1] = 3.062999454; //b
0081      aY[1][0] = 0.0204108951; //c
0082      aY[1][1] = -0.139318692; //d
0083
0084      local_x_offset       = -0.339334;
0085      local_y_offset       = -0.000299454;
0086      local_x_slope_offset = -0.219603248;
0087      local_y_slope_offset = -0.000176128;
0088
0089   }
0090   else if(abs(100.0 - nomMomentum)/100.0 < nomMomentumError){
0091
0092      aX[0][0] = 3.152158; //a
0093      aX[0][1] = 20.852072; //b
0094      aX[1][0] = 0.181649517; //c
0095      aX[1][1] = -0.303998487; //d
0096
0097      aY[0][0] = 0.5306100000; //a
0098      aY[0][1] = 3.19623343; //b
0099      aY[1][0] = 0.0226283320; //c
0100      aY[1][1] = -0.082666019; //d
0101
0102      local_x_offset       = -0.329072;
0103      local_y_offset       = -0.00028343;
0104      local_x_slope_offset = -0.218525084;
0105      local_y_slope_offset = -0.00015321;
0106
0107   }
0108   else if(abs(41.0 - nomMomentum)/41.0 < nomMomentumError){
0109
0110          aX[0][0] = 3.135997; //a
0111          aX[0][1] = 18.482273; //b
0112          aX[1][0] = 0.176479921; //c
0113          aX[1][1] = -0.497839483; //d
0114
0115          aY[0][0] = 0.4914400000; //a
0116          aY[0][1] = 4.53857451; //b
0117          aY[1][0] = 0.0179664765; //c
0118          aY[1][1] = 0.004160679; //d
0119
0120          local_x_offset       = -0.283273;
0121          local_y_offset       = -0.00552451;
0122          local_x_slope_offset = -0.21174031;
0123          local_y_slope_offset = -0.003212011;
0124
0125   }
0126   else if(abs(135.0 - nomMomentum)/135.0 < nomMomentumError){ //135 GeV deuterons
0127
0128       aX[0][0] = 1.6248;
0129       aX[0][1] = 12.966293;
0130       aX[1][0] = 0.1832;
0131       aX[1][1] = -2.8636535;
0132
0133       aY[0][0] = 0.0001674; //a
0134       aY[0][1] = -28.6003; //b
0135       aY[1][0] = 0.0000837; //c
0136       aY[1][1] = -2.87985; //d
0137
0138       local_x_offset       = -11.9872;
0139       local_y_offset       = -0.0146;
0140       local_x_slope_offset = -14.75315;
0141       local_y_slope_offset = -0.0073;
0142
0143   }
0144   else {
0145     error("MatrixTransferStatic:: No valid matrix found to match beam momentum!! Skipping!!");
0146     return;
0147   }
0148
0149   double determinant = aX[0][0] * aX[1][1] - aX[0][1] * aX[1][0];
0150
0151   if (determinant == 0) {
0152     error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0153     return;
0154   }
0155
0156   aXinv[0][0] =  aX[1][1] / determinant;
0157   aXinv[0][1] = -aX[0][1] / determinant;
0158   aXinv[1][0] = -aX[1][0] / determinant;
0159   aXinv[1][1] =  aX[0][0] / determinant;
0160
0161
0162   determinant = aY[0][0] * aY[1][1] - aY[0][1] * aY[1][0];
0163
0164   if (determinant == 0) {
0165     error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0166     return;
0167   }
0168
0169   aYinv[0][0] =  aY[1][1] / determinant;
0170   aYinv[0][1] = -aY[0][1] / determinant;
0171   aYinv[1][0] = -aY[1][0] / determinant;
0172   aYinv[1][1] =  aY[0][0] / determinant;
0173
0174   //---- begin Reconstruction code ----
0175
0176   edm4hep::Vector3f goodHit[2] = {{0.0,0.0,0.0},{0.0,0.0,0.0}};
0177
0178   double goodHitX[2] = {0.0, 0.0};
0179   double goodHitY[2] = {0.0, 0.0};
0180   double goodHitZ[2] = {0.0, 0.0};
0181
0182   bool goodHit1 = false;
0183   bool goodHit2 = false;
0184
0185   for (const auto &h: *rechits) {
0186
0187     auto cellID = h.getCellID();
0188     // The actual hit position in Global Coordinates
0189     auto gpos = m_converter->position(cellID);
0190     // local positions
0191     auto volman = m_detector->volumeManager();
0192     auto local = volman.lookupDetElement(cellID);
0193
0194     auto pos0 = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); // hit position in local coordinates
0195
0196     // convert into mm
0197     gpos = gpos/dd4hep::mm;
0198     pos0 = pos0/dd4hep::mm;
0199
0200     //std::cout << "gpos.z() = " << gpos.z() << " pos0.z() = " << pos0.z() << std::endl;
0201     //std::cout << "[gpos.x(), gpos.y()] = " << gpos.x() <<", "<< gpos.y() << "  and [pos0.x(), pos0.y()] = "<< pos0.x()<< ", " << pos0.y() << std::endl;
0202
0203     if(!goodHit2 && gpos.z() > m_cfg.hit2minZ && gpos.z() < m_cfg.hit2maxZ){
0204
0205       goodHit[1].x = pos0.x();
0206       goodHit[1].y = pos0.y();
0207       goodHit[1].z = gpos.z();
0208       goodHit2 = true;
0209
0210     }
0211     if(!goodHit1 && gpos.z() > m_cfg.hit1minZ && gpos.z() < m_cfg.hit1maxZ){
0212
0213       goodHit[0].x = pos0.x();
0214       goodHit[0].y = pos0.y();
0215       goodHit[0].z = gpos.z();
0216       goodHit1 = true;
0217
0218     }
0219
0220   }
0221
0222   // NB:
0223   // This is a "dumb" algorithm - I am just checking the basic thing works with a simple single-proton test.
0224   // Will need to update and modify for generic # of hits for more complicated final-states.
0225
0226   if (goodHit1 && goodHit2) {
0227
0228     // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit
0229     // trajectory
0230     double XL[2] = {goodHit[0].x - local_x_offset, goodHit[1].x - local_x_offset};
0231     double YL[2] = {goodHit[0].y - local_y_offset, goodHit[1].y - local_y_offset};
0232
0233     double base = goodHit[1].z - goodHit[0].z;
0234
0235     if (base == 0) {
0236       info("Detector separation = 0! Cannot calculate slope!");
0237     }
0238     else{
0239
0240       double Xip[2] = {0.0, 0.0};
0241       double Xrp[2] = {XL[1], ((XL[1] - XL[0]) / (base))/dd4hep::mrad - local_x_slope_offset};
0242       double Yip[2] = {0.0, 0.0};
0243       double Yrp[2] = {YL[1], ((YL[1] - YL[0]) / (base))/dd4hep::mrad - local_y_slope_offset};
0244
0245       // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the
0246       // Polar Angle and deltaP at the IP
0247
0248       for (unsigned i0 = 0; i0 < 2; i0++) {
0249         for (unsigned i1 = 0; i1 < 2; i1++) {
0250           Xip[i0] += aXinv[i0][i1] * Xrp[i1];
0251           Yip[i0] += aYinv[i0][i1] * Yrp[i1];
0252         }
0253       }
0254
0255       // convert polar angles to radians
0256       double rsx = Xip[1] * dd4hep::mrad;
0257       double rsy = Yip[1] * dd4hep::mrad;
0258
0259       // calculate momentum magnitude from measured deltaP – using thin lens optics.
0260       double p = nomMomentum * (1 + 0.01 * Xip[0]);
0261       double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
0262
0263       edm4hep::Vector3f prec = {static_cast<float>(p * rsx / norm), static_cast<float>(p * rsy / norm),
0264                                 static_cast<float>(p / norm)};
0265       auto refPoint = goodHit[0];
0266
0267       //----- end reconstruction code ------
0268
0269       edm4eic::MutableReconstructedParticle reconTrack;
0270       reconTrack.setType(0);
0271       reconTrack.setMomentum(prec);
0272       reconTrack.setEnergy(std::hypot(edm4hep::utils::magnitude(reconTrack.getMomentum()), m_cfg.partMass));
0273       reconTrack.setReferencePoint(refPoint);
0274       reconTrack.setCharge(m_cfg.partCharge);
0275       reconTrack.setMass(m_cfg.partMass);
0276       reconTrack.setGoodnessOfPID(1.);
0277       reconTrack.setPDG(m_cfg.partPDG);
0278       //reconTrack.covMatrix(); // @TODO: Errors
0279       outputParticles->push_back(reconTrack);
0280     }
0281   } // end enough hits for matrix reco
0282
0283 }