algorithms/reco/FarForwardLambdaReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Sebouh Paul
0003
0004 #include <Evaluator/DD4hepUnits.h>
0005 #include <TVector3.h>
0006 #include <edm4eic/ReconstructedParticleCollection.h>
0007 #include <edm4hep/Vector3f.h>
0008 #include <edm4hep/utils/vector_utils.h>
0009 #include <fmt/core.h>
0010 #include <math.h>
0011 #include <cstddef>
0012 #include <gsl/pointers>
0013 #include <vector>
0014
0015 #include "FarForwardLambdaReconstruction.h"
0016 #include "TLorentzVector.h"
0017 /**
0018 Creates Lambda candidates from a neutron and two photons from a pi0 decay
0019  */
0020
0021 namespace eicrecon {
0022
0023     void FarForwardLambdaReconstruction::init() {
0024
0025     }
0026
0027   /* converts one type of vector format to another */
0028   void toTVector3(TVector3& v1, const edm4hep::Vector3f& v2){
0029     v1.SetXYZ(v2.x,v2.y,v2.z);
0030   }
0031
0032     void FarForwardLambdaReconstruction::process(const FarForwardLambdaReconstruction::Input& input,
0033                       const FarForwardLambdaReconstruction::Output& output) const {
0034       const auto [neutrals] = input;
0035       auto [out_lambdas, out_decay_products] = output;
0036       std::vector<edm4eic::ReconstructedParticle> neutrons={};
0037       std::vector<edm4eic::ReconstructedParticle> gammas={};
0038       for (auto part: *neutrals){
0039         if (part.getPDG()==2112){
0040           neutrons.push_back(part);
0041         }
0042         if (part.getPDG()==22){
0043           gammas.push_back(part);
0044         }
0045       }
0046
0047       if (neutrons.size()<1 || gammas.size()<2)
0048         return;
0049
0050       static const double m_neutron = m_particleSvc.particle(2112).mass;
0051       static const double m_pi0 = m_particleSvc.particle(111).mass;
0052       static const double m_lambda = m_particleSvc.particle(3122).mass;
0053
0054       for (std::size_t i_n=0; i_n<neutrons.size(); i_n++){
0055         for (std::size_t i_1=0; i_1<gammas.size()-1; i_1++){
0056           for (std::size_t i_2=i_1+1; i_2<gammas.size(); i_2++){
0057
0058             double En=neutrons[i_n].getEnergy();
0059             double pn=sqrt(En*En-m_neutron*m_neutron);
0060             double E1=gammas[i_1].getEnergy();
0061             double E2=gammas[i_2].getEnergy();
0062             TVector3 xn;
0063             TVector3 x1;
0064             TVector3 x2;
0065
0066             toTVector3(xn,neutrons[0].getReferencePoint()*dd4hep::mm);
0067             toTVector3(x1,gammas[0].getReferencePoint()*dd4hep::mm);
0068             toTVector3(x2,gammas[1].getReferencePoint()*dd4hep::mm);
0069
0070             xn.RotateY(-m_cfg.globalToProtonRotation);
0071             x1.RotateY(-m_cfg.globalToProtonRotation);
0072             x2.RotateY(-m_cfg.globalToProtonRotation);
0073
0074             debug("nx recon = {}, g1x recon = {}, g2x recon = {}", xn.X(), x1.X(), x2.X());
0075             debug("nz recon = {}, g1z recon = {}, g2z recon = {}, z face = {}", xn.Z(), x1.Z(), x2.Z(), m_cfg.zMax);
0076
0077             TVector3 vtx(0,0,0);
0078             double f=0;
0079             double df=0.5;
0080             double theta_open_expected=2*asin(m_pi0/(2*sqrt(E1*E2)));
0081             TLorentzVector n, g1, g2, lambda;
0082             for(int i = 0; i<m_cfg.iterations; i++){
0083               n={pn*(xn-vtx).Unit(), En};
0084               g1={E1*(x1-vtx).Unit(), E1};
0085               g2={E2*(x2-vtx).Unit(), E2};
0086               double theta_open=g1.Angle(g2.Vect());
0087               lambda=n+g1+g2;
0088               if (theta_open>theta_open_expected) {
0089                 f-=df;
0090               } else if (theta_open<theta_open_expected) {
0091                 f+=df;
0092               }
0093
0094               vtx=lambda.Vect()*(f*m_cfg.zMax/lambda.Z());
0095               df/=2;
0096             }
0097
0098             double mass_rec=lambda.M();
0099             if (abs(mass_rec-m_lambda)>m_cfg.lambdaMaxMassDev)
0100               continue;
0101
0102             // rotate everything back to the lab coordinates.
0103             vtx.RotateY(m_cfg.globalToProtonRotation);
0104             lambda.RotateY(m_cfg.globalToProtonRotation);
0105             n.RotateY(m_cfg.globalToProtonRotation);
0106             g1.RotateY(m_cfg.globalToProtonRotation);
0107             g2.RotateY(m_cfg.globalToProtonRotation);
0108
0109             auto b=-lambda.BoostVector();
0110             n.Boost(b);
0111             g1.Boost(b);
0112             g2.Boost(b);
0113
0114             //convert vertex to mm:
0115             vtx=vtx*(1/dd4hep::mm);
0116
0117             auto rec_lambda = out_lambdas->create();
0118             rec_lambda.setPDG(3122);
0119
0120             rec_lambda.setEnergy(lambda.E());
0121             rec_lambda.setMomentum({static_cast<float>(lambda.X()), static_cast<float>(lambda.Y()), static_cast<float>(lambda.Z())});
0122             rec_lambda.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
0123             rec_lambda.setCharge(0);
0124             rec_lambda.setMass(mass_rec);
0125
0126             auto neutron_cm = out_decay_products->create();
0127             neutron_cm.setPDG(2112);
0128             neutron_cm.setEnergy(n.E());
0129             neutron_cm.setMomentum({static_cast<float>(n.X()), static_cast<float>(n.Y()), static_cast<float>(n.Z())});
0130             neutron_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
0131             neutron_cm.setCharge(0);
0132             neutron_cm.setMass(m_neutron);
0133             //link the reconstructed lambda to the input neutron,
0134             // and the cm neutron to the input neutron
0135             rec_lambda.addToParticles(neutrons[i_n]);
0136             neutron_cm.addToParticles(neutrons[i_n]);
0137
0138             auto gamma1_cm = out_decay_products->create();
0139             gamma1_cm.setPDG(22);
0140             gamma1_cm.setEnergy(g1.E());
0141             gamma1_cm.setMomentum({static_cast<float>(g1.X()), static_cast<float>(g1.Y()), static_cast<float>(g1.Z())});
0142             gamma1_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
0143             gamma1_cm.setCharge(0);
0144             gamma1_cm.setMass(0);
0145             rec_lambda.addToParticles(gammas[i_1]);
0146             gamma1_cm.addToParticles(gammas[i_1]);
0147
0148             auto gamma2_cm = out_decay_products->create();
0149             gamma2_cm.setPDG(22);
0150             gamma2_cm.setEnergy(g2.E());
0151             gamma2_cm.setMomentum({static_cast<float>(g2.X()), static_cast<float>(g2.Y()), static_cast<float>(g2.Z())});
0152             gamma2_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
0153             gamma2_cm.setCharge(0);
0154             gamma2_cm.setMass(0);
0155             rec_lambda.addToParticles(gammas[i_2]);
0156             gamma2_cm.addToParticles(gammas[i_2]);
0157             continue;
0158           }
0159         }
0160       }
0161     }
0162 }