File indexing completed on 2024-09-27 07:03:27
0001
0002
0003 #ifndef EXCLUDE_DELPHES
0004
0005
0006
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 {};
0014
0015
0016
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
0032
0033 quarkpT = partTrue->PT;
0034 }
0035 }
0036
0037 std::vector<fastjet::PseudoJet> particles;
0038 std::vector<fastjet::PseudoJet> particlesTrue;
0039 jetConstituents.clear();
0040
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
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
0096 fastjet::JetDefinition jet_def(fastjet::kt_algorithm, jetRadius);
0097
0098
0099
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
0113 TLorentzVector pjet(jet.px(), jet.py(), jet.pz(), jet.E());
0114 TVector3 qTjetVect( vecElectron.Px()+pjet.Px(), vecElectron.Py()+pjet.Py(), 0);
0115 qTjet = qTjetVect.Mag();
0116
0117 zjet = (vecIonBeam.Dot(pjet))/((vecIonBeam).Dot(vecQ));
0118 pTjet = jet.pt();
0119 mTjet = jet.mt();
0120 etajet = jet.eta();
0121 phijet = jet.phi();
0122 mjet = jet.m();
0123 ejet = jet.e();
0124
0125 jperp.clear();
0126 zhad_jet.clear();
0127 std::vector<fastjet::PseudoJet> constituents = jet.constituents();
0128 int constituentPID = 0;
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
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
0174 if(matchIndex == -1)
0175 deltaRjet = -1.0;
0176 else
0177 deltaRjet = minDeltaR;
0178
0179
0180 if(minDeltaR < deltaRCut && matchIndex != -1) matchStatusjet = 0;
0181 if(minDeltaR > deltaRCut && matchIndex != -1) matchStatusjet = 1;
0182 if(matchIndex == -1) matchStatusjet = 2;
0183 if(minDeltaR < deltaRCut && matchIndex != -1)
0184 {
0185 pTmtjet = jetsTrue[matchIndex].pt();
0186 mTmtjet = jetsTrue[matchIndex].mt();
0187 etamtjet = jetsTrue[matchIndex].eta();
0188 phimtjet = jetsTrue[matchIndex].phi();
0189 mmtjet = jetsTrue[matchIndex].m();
0190 emtjet= jetsTrue[matchIndex].e();
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(¢auroPlugin);
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;
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
0348
0349
0350
0351 KinematicsJets::~KinematicsJets() {};
0352 #endif