Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-06-17 07:06:50 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Kevin Adkins, Brian Page, Connor Pecar
0003 #ifndef EXCLUDE_DELPHES
0004 
0005 /* NOTE:
0006  * if you make changes, MAINTAIN DOCUMENTATION IN ../doc/kinematicsJets.md
0007  */
0008 
0009 #include "KinematicsJets.h"
0010 ClassImp(KinematicsJets)
0011 
0012 KinematicsJets::KinematicsJets(Double_t enEleBeam_, Double_t enIonBeam_, Double_t crossAng_) : Kinematics(enEleBeam_,enIonBeam_,crossAng_)
0013 {/*Note: Energy units are GeV, angle units are mrad*/};
0014 
0015 
0016 // JET METHODS ///////////////////////
0017 
0018 void KinematicsJets::GetJets(
0019     TObjArrayIter itEFlowTrack, TObjArrayIter itEFlowPhoton,
0020     TObjArrayIter itEFlowNeutralHadron, TObjArrayIter itParticle,
0021     int jetAlgo, double jetRadius, double jetMinPt
0022     )
0023 {
0024   itEFlowTrack.Reset();
0025   itEFlowPhoton.Reset();
0026   itEFlowNeutralHadron.Reset();
0027   itParticle.Reset();
0028 
0029   while(GenParticle *partTrue = (GenParticle*)itParticle() ){
0030     if( (partTrue->PID == 1 || partTrue->PID == 2) && (partTrue->Status == 23) ){
0031       // Status: 23->outgoing, but there's also 63->outgoing beam remnant. TODO: Which do we want?
0032       // from pythia 8 documentation
0033       quarkpT = partTrue->PT;
0034     }
0035   }
0036 
0037   std::vector<fastjet::PseudoJet> particles;
0038   std::vector<fastjet::PseudoJet> particlesTrue;
0039   jetConstituents.clear();
0040   // looping over final state particles, adding to particles vector
0041   while(Track *eflowTrack = (Track*)itEFlowTrack() ){
0042     TLorentzVector eflowTrackp4 = eflowTrack->P4();
0043     GenParticle *trackParticle = (GenParticle*)eflowTrack->Particle.GetObject();
0044     int partPID = std::abs(trackParticle->PID);
0045     if(!isnan(eflowTrackp4.E())){
0046       if(std::abs(eflowTrack->Eta) < 4.0 && eflowTrack->PT > 0.1){
0047         this->TransformToHeadOnFrame(eflowTrackp4,eflowTrackp4);
0048         if(partPID != 11) particles.push_back(fastjet::PseudoJet(eflowTrackp4.Px(),eflowTrackp4.Py(),eflowTrackp4.Pz(),eflowTrackp4.E()));
0049 
0050         //GenParticle *trackParticle = (GenParticle*)eflowTrack->Particle.GetObject();
0051         TLorentzVector partp4 = trackParticle->P4();
0052         this->TransformToHeadOnFrame(partp4,partp4);
0053         if(partPID != 11) particlesTrue.push_back(fastjet::PseudoJet(partp4.Px(),partp4.Py(),partp4.Pz(),partp4.E()));
0054 
0055         jetConstituents.insert({eflowTrackp4.Px(), eflowTrack->PID});
0056       }
0057     }
0058   }
0059   while(Tower* towerPhoton = (Tower*)itEFlowPhoton() ){
0060     TLorentzVector  towerPhotonp4 = towerPhoton->P4();
0061     if(!isnan(towerPhotonp4.E())){
0062       if( std::abs(towerPhoton->Eta) < 4.0){
0063         this->TransformToHeadOnFrame(towerPhotonp4,towerPhotonp4);
0064         particles.push_back(fastjet::PseudoJet(towerPhotonp4.Px(),towerPhotonp4.Py(),towerPhotonp4.Pz(),towerPhotonp4.E()));
0065 
0066         for(int i = 0; i < towerPhoton->Particles.GetEntries(); i++){
0067           GenParticle *photonPart = (GenParticle*)towerPhoton->Particles.At(i);
0068           TLorentzVector photonp4 = photonPart->P4();
0069           this->TransformToHeadOnFrame(photonp4,photonp4);
0070           particlesTrue.push_back(fastjet::PseudoJet(photonp4.Px(),photonp4.Py(),photonp4.Pz(),photonp4.E()));
0071         }
0072       }
0073     }
0074   }
0075 
0076   while(Tower* towerNeutralHadron = (Tower*)itEFlowNeutralHadron() ){
0077     TLorentzVector  towerNeutralHadronp4 = towerNeutralHadron->P4();
0078     if(!isnan(towerNeutralHadronp4.E())){
0079       if( std::abs(towerNeutralHadron->Eta) < 4.0){
0080         this->TransformToHeadOnFrame(towerNeutralHadronp4,towerNeutralHadronp4);
0081         particles.push_back(
0082           fastjet::PseudoJet(towerNeutralHadronp4.Px(),towerNeutralHadronp4.Py(),towerNeutralHadronp4.Pz(),towerNeutralHadronp4.E())
0083           );
0084 
0085         for(int i = 0; i < towerNeutralHadron->Particles.GetEntries(); i++){
0086           GenParticle *nhadPart = (GenParticle*)towerNeutralHadron->Particles.At(i);
0087           TLorentzVector nhadp4 = nhadPart->P4();
0088           this->TransformToHeadOnFrame(nhadp4,nhadp4);
0089           particlesTrue.push_back(fastjet::PseudoJet(nhadp4.Px(),nhadp4.Py(),nhadp4.Pz(),nhadp4.E()));
0090         }
0091       }
0092     }
0093   }
0094 
0095   // Set Jet Definition
0096   fastjet::JetDefinition jet_def(fastjet::kt_algorithm, jetRadius);
0097 
0098   // Redefine Jet Algo Based on Macro Input
0099   // Default is jetAlgo = 0 which is kt_algorithm
0100   if(jetAlgo == 1) jet_def.set_jet_algorithm(fastjet::cambridge_algorithm);
0101   if(jetAlgo == 2) jet_def.set_jet_algorithm(fastjet::antikt_algorithm);
0102 
0103   csRec = fastjet::ClusterSequence(particles, jet_def);
0104   csTrue = fastjet::ClusterSequence(particlesTrue, jet_def);
0105   jetsRec = sorted_by_pt(csRec.inclusive_jets(jetMinPt));
0106   jetsTrue = sorted_by_pt(csTrue.inclusive_jets(jetMinPt));
0107 
0108 };
0109 
0110 
0111 void KinematicsJets::CalculateJetKinematics(fastjet::PseudoJet jet){
0112   // `jet` is already in the head-on frame, since `jetsRec` was filled with head-on frame momenta
0113   TLorentzVector pjet(jet.px(), jet.py(), jet.pz(), jet.E());
0114   TVector3 qTjetVect( vecElectron.Px()+pjet.Px(), vecElectron.Py()+pjet.Py(), 0); // (used only in Lorentz invariant calculations)
0115   qTjet = qTjetVect.Mag();
0116 
0117   zjet = (vecIonBeam.Dot(pjet))/((vecIonBeam).Dot(vecQ));
0118   pTjet = jet.pt(); // lab frame pT
0119   mTjet = jet.mt(); // Transverse Mass
0120   etajet = jet.eta(); // Jet Pseudorapidity
0121   phijet = jet.phi(); // Jet Phi
0122   mjet = jet.m(); // Jet Mass
0123   ejet = jet.e(); // Jet Energy
0124 
0125   jperp.clear();
0126   zhad_jet.clear();
0127   std::vector<fastjet::PseudoJet> constituents = jet.constituents();
0128   int constituentPID = 0; // if we only want zh/jperp for pi+, other tracks
0129 
0130   if(constituentPID == 0){
0131     for(int i = 0; i < constituents.size(); i++){
0132       TLorentzVector partVec(constituents[i].px(), constituents[i].py(), constituents[i].pz(), constituents[i].E());
0133       TVector3 jperpVec = Reject(partVec.Vect(),pjet.Vect());
0134       jperp.push_back(jperpVec.Mag());
0135       zhad_jet.push_back( (partVec.Vect()).Mag()/((pjet.Vect()).Mag()) );
0136     }
0137   }
0138   else {
0139     for(int i = 0; i < constituents.size(); i++){
0140       TLorentzVector partVec(constituents[i].px(), constituents[i].py(), constituents[i].pz(), constituents[i].E());
0141       std::map<double,int>::iterator it;
0142       it = jetConstituents.find(partVec.Px());
0143       if(it != jetConstituents.end()){
0144         int pidTrack = it->second;
0145         if(pidTrack == constituentPID){
0146           TVector3 jperpVec = Reject(partVec.Vect(),pjet.Vect());
0147           jperp.push_back(jperpVec.Mag());
0148           zhad_jet.push_back( (partVec.Vect()).Mag()/((pjet.Vect()).Mag()) );
0149         }
0150       }
0151     }
0152   }
0153 };
0154 
0155 
0156 void KinematicsJets::CalculateJetResolution(double deltaRCut){
0157   // Find Matching Truth-level Jet for each Reco Jet
0158   double minDeltaR = 100000.0;
0159   double matchIndex = -1;
0160   for(unsigned int jt=0; jt<jetsTrue.size(); jt++)
0161     {
0162       double deltaEta = etajet - jetsTrue[jt].eta();
0163       double deltaPhi = TVector2::Phi_mpi_pi(phijet - jetsTrue[jt].phi());
0164       double deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
0165 
0166       if(deltaR < minDeltaR)
0167     {
0168       minDeltaR = deltaR;
0169       matchIndex = jt;
0170     }
0171     }
0172 
0173   // Pass DeltaR to Saved Variables
0174   if(matchIndex == -1)
0175     deltaRjet = -1.0;
0176   else
0177     deltaRjet = minDeltaR;
0178 
0179   // Calculate Resolution for Matched Jets
0180   if(minDeltaR < deltaRCut && matchIndex != -1) matchStatusjet = 0; // Matched Jet Found
0181   if(minDeltaR > deltaRCut && matchIndex != -1) matchStatusjet = 1; // True Jet Found, but Outside Match Criteria
0182   if(matchIndex == -1) matchStatusjet = 2; // No True Jet Found
0183   if(minDeltaR < deltaRCut && matchIndex != -1)
0184     {
0185       pTmtjet = jetsTrue[matchIndex].pt(); // lab frame matched true jet pT
0186       mTmtjet = jetsTrue[matchIndex].mt(); // Matched true Jet Transverse Mass
0187       etamtjet = jetsTrue[matchIndex].eta(); // Matched true Jet Pseudorapidity
0188       phimtjet = jetsTrue[matchIndex].phi(); // Matched true Jet Phi
0189       mmtjet = jetsTrue[matchIndex].m(); // Matched True Jet Mass
0190       emtjet= jetsTrue[matchIndex].e(); // Matched true Jet Energy
0191     }
0192 };
0193 
0194 
0195 #ifdef INCLUDE_CENTAURO
0196 void KinematicsJets::GetBreitFrameJets(
0197   TObjArrayIter itEFlowTrack, TObjArrayIter itEFlowPhoton,
0198   TObjArrayIter itEFlowNeutralHadron, TObjArrayIter itParticle
0199   )
0200 {
0201   itEFlowTrack.Reset();
0202   itEFlowPhoton.Reset();
0203   itEFlowNeutralHadron.Reset();
0204   itParticle.Reset();
0205   std::vector<fastjet::PseudoJet> particles;
0206   std::vector<fastjet::PseudoJet> particlesTrue;
0207 
0208   jetConstituents.clear();
0209 
0210   double highPT = -1;
0211   TLorentzVector eleTrue;
0212   while(GenParticle* part = (GenParticle*) itParticle()){
0213     if(part->PID == 11){
0214       if(part->PT > highPT){
0215         highPT = part->PT;
0216         eleTrue = part->P4();
0217       }
0218     }
0219   }
0220   TLorentzVector vecQTrue = vecEleBeam - eleTrue;
0221   double Q2true = -1*vecQTrue.M2();
0222   double xtrue = Q2true / ( 2 * vecQTrue.Dot(vecIonBeam) );
0223 
0224   TLorentzVector breitVecTrue = vecQTrue + 2*xtrue*vecIonBeam;
0225   TLorentzVector breitVec = vecQ + 2*x*vecIonBeam;
0226   TVector3 breitBoostTrue = -1*breitVecTrue.BoostVector();
0227   TVector3 breitBoost = -1*breitVec.BoostVector();
0228   itParticle.Reset();
0229 
0230   while(Track *eflowTrack = (Track*)itEFlowTrack() ){
0231     TLorentzVector eflowTrackp4 = eflowTrack->P4();
0232     if(!isnan(eflowTrackp4.E()) && eflowTrackp4 != vecElectron){
0233       if(std::abs(eflowTrack->Eta) < 4.0 && eflowTrack->PT > 0.2){
0234         eflowTrackp4.Boost(breitBoost);
0235         particles.push_back(fastjet::PseudoJet(eflowTrackp4.Px(),eflowTrackp4.Py(),eflowTrackp4.Pz(),eflowTrackp4.E()));
0236 
0237         GenParticle *trackParticle = (GenParticle*)eflowTrack->Particle.GetObject();
0238         TLorentzVector partp4 = trackParticle->P4();
0239         partp4.Boost(breitBoostTrue);
0240         particlesTrue.push_back(fastjet::PseudoJet(partp4.Px(),partp4.Py(),partp4.Pz(),partp4.E()));
0241 
0242         jetConstituents.insert({eflowTrackp4.Px(), eflowTrack->PID});
0243 
0244       }
0245     }
0246   }
0247   while(Tower* towerPhoton = (Tower*)itEFlowPhoton() ){
0248     TLorentzVector  towerPhotonp4 = towerPhoton->P4();
0249     if(!isnan(towerPhotonp4.E())){
0250       if( std::abs(towerPhoton->Eta) < 4.0 &&
0251           sqrt(towerPhotonp4.Px()*towerPhotonp4.Px()+towerPhotonp4.Py()*towerPhotonp4.Py()) > 0.2
0252           )
0253       {
0254         towerPhotonp4.Boost(breitBoost);
0255         particles.push_back(fastjet::PseudoJet(towerPhotonp4.Px(),towerPhotonp4.Py(),towerPhotonp4.Pz(),towerPhotonp4.E()));
0256 
0257         for(int i = 0; i < towerPhoton->Particles.GetEntries(); i++){
0258           GenParticle *photonPart = (GenParticle*)towerPhoton->Particles.At(i);
0259           TLorentzVector photonp4 = photonPart->P4();
0260           photonp4.Boost(breitBoostTrue);
0261           particlesTrue.push_back(fastjet::PseudoJet(photonp4.Px(),photonp4.Py(),photonp4.Pz(),photonp4.E()));
0262         }
0263       }
0264     }
0265   }
0266   while(Tower* towerNeutralHadron = (Tower*)itEFlowNeutralHadron() ){
0267     TLorentzVector  towerNeutralHadronp4 = towerNeutralHadron->P4();
0268     if( !isnan(towerNeutralHadronp4.E()) &&
0269         sqrt(towerNeutralHadronp4.Px()*towerNeutralHadronp4.Px()+towerNeutralHadronp4.Py()*towerNeutralHadronp4.Py()) > 0.2
0270         )
0271     {
0272       if( std::abs(towerNeutralHadron->Eta) < 4.0 ){
0273         towerNeutralHadronp4.Boost(breitBoost);
0274         particles.push_back(
0275             fastjet::PseudoJet(towerNeutralHadronp4.Px(),towerNeutralHadronp4.Py(),towerNeutralHadronp4.Pz(),towerNeutralHadronp4.E())
0276           );
0277 
0278         for(int i = 0; i < towerNeutralHadron->Particles.GetEntries(); i++){
0279           GenParticle *nhadPart = (GenParticle*)towerNeutralHadron->Particles.At(i);
0280           TLorentzVector nhadp4 = nhadPart->P4();
0281           nhadp4.Boost(breitBoostTrue);
0282           particlesTrue.push_back(fastjet::PseudoJet(nhadp4.Px(),nhadp4.Py(),nhadp4.Pz(),nhadp4.E()));
0283         }
0284       }
0285     }
0286   }
0287 
0288   double R = 0.8;
0289   contrib::CentauroPlugin centauroPlugin(R);
0290   fastjet::JetDefinition jet_def(&centauroPlugin);
0291 
0292   csRec = fastjet::ClusterSequence(particles, jet_def);
0293   csTrue = fastjet::ClusterSequence(particlesTrue, jet_def);
0294   breitJetsRec = sorted_by_pt(csRec.inclusive_jets());
0295   breitJetsTrue = sorted_by_pt(csTrue.inclusive_jets());
0296 };
0297 
0298 
0299 void KinematicsJets::CalculateBreitJetKinematics(fastjet::PseudoJet jet){
0300   TLorentzVector pjet(jet.px(), jet.py(), jet.pz(), jet.E());
0301   TLorentzVector pjetLab = pjet;
0302 
0303   TLorentzVector breitVec = vecQ + 2*x*vecIonBeam;
0304   TVector3 breitBoost = -1*breitVec.BoostVector();
0305 
0306   pjetLab.Boost(-1*breitBoost);
0307   pTjet = sqrt(pjetLab.Px()*pjetLab.Px() + pjetLab.Py()*pjetLab.Py());
0308 
0309   TLorentzVector vecElectronBreit = vecElectron;
0310   vecElectronBreit.Boost(breitBoost);
0311   TVector3 qTjetVect(vecElectronBreit.Px()+pjet.Px(), vecElectronBreit.Py()+pjet.Py(), 0);
0312   qTjet = qTjetVect.Mag();
0313 
0314   TLorentzVector nbreit(0,0,1/sqrt(Q2),1/sqrt(Q2));
0315   double zjet = nbreit*pjet;
0316 
0317   jperp.clear();
0318   zhad_jet.clear();
0319   std::vector<fastjet::PseudoJet> constituents = jet.constituents();
0320   int constituentPID = 0; // if we only want zh/jperp for pi+, other tracks
0321 
0322   if(constituentPID == 0){
0323     for(int i = 0; i < constituents.size(); i++){
0324       TLorentzVector partVec(constituents[i].px(), constituents[i].py(), constituents[i].pz(), constituents[i].E());
0325       TVector3 jperpVec = Reject(partVec.Vect(),pjet.Vect());
0326       jperp.push_back(jperpVec.Mag());
0327       zhad_jet.push_back( (partVec.Vect()).Mag()/((pjet.Vect()).Mag()) );
0328     }
0329   }
0330   else{
0331     for(int i = 0; i < constituents.size(); i++){
0332       TLorentzVector partVec(constituents[i].px(), constituents[i].py(), constituents[i].pz(), constituents[i].E());
0333       std::map<double,int>::iterator it;
0334       it = jetConstituents.find(partVec.Px());
0335       if(it != jetConstituents.end()){
0336         int pidTrack = it->second;
0337         if( pidTrack == constituentPID){
0338           TVector3 jperpVec = Reject(partVec.Vect(),pjet.Vect());
0339           jperp.push_back(jperpVec.Mag());
0340           zhad_jet.push_back( (partVec.Vect()).Mag()/((pjet.Vect()).Mag()) );
0341         }
0342       }
0343     }
0344   }
0345 
0346 };
0347 #endif // ifdef INCLUDE_CENTAURO
0348 // end DELPHES-only methods //////////////////////////////////////////////////////
0349 
0350 
0351 KinematicsJets::~KinematicsJets() {};
0352 #endif // ifndef EXCLUDE_DELPHES
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.5 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!