analysis/EICAnalysis/DVMPHelper.C

0001 #include "DVMPHelper.h"
0002
0003 #include <cmath>
0004 #include <iostream>
0005 DVMPHelper::DVMPHelper(std::vector<float> reco_eta,
0006                std::vector<float> reco_phi,
0007                std::vector<float> reco_ptotal,
0008                std::vector<int> reco_charge,
0009                std::vector<float> reco_cluster_e,
0010                std::vector<bool> reco_is_scattered_lepton,
0011                std::vector<float> true_eta,
0012                std::vector<float> true_phi,
0013                std::vector<float> true_ptotal,
0014                std::vector<int> pid,
0015                std::vector<bool> is_scattered_lepton)
0016 {
0017
0018   _size_reco = reco_eta.size();
0019   _size_truth = true_eta.size();
0020
0021   rparticles = (DVMPHelper::particle_reco*)malloc(_size_reco * sizeof(DVMPHelper::particle_reco));
0022   tparticles = (DVMPHelper::particle_truth*)malloc(_size_truth * sizeof(DVMPHelper::particle_truth));
0023
0024   for(int i = 0 ; i < _size_reco ; i++)
0025     {
0026       rparticles[i]=GetParticleReco(reco_eta.at(i),reco_phi.at(i),reco_ptotal.at(i),reco_charge.at(i),reco_cluster_e.at(i),reco_is_scattered_lepton.at(i));
0027     }
0028   for(int j = 0 ; j <_size_truth ; j++)
0029     {
0030       tparticles[j]=GetParticleTruth(true_eta.at(j),true_phi.at(j),true_ptotal.at(j),pid.at(j),is_scattered_lepton.at(j));
0031     }
0032 }
0033
0034 DVMPHelper::particle_reco
0035 DVMPHelper::GetParticleReco(float eta, float phi, float ptotal, int charge,float e,bool is_scattered_lepton)
0036 {
0037   DVMPHelper::particle_reco particle = {eta,phi,ptotal,charge,e,is_scattered_lepton};
0038   return particle;
0039 }
0040
0041 DVMPHelper::particle_truth
0042 DVMPHelper::GetParticleTruth(float eta, float phi, float ptotal, int pid, bool is_scattered_lepton)
0043 {
0044   DVMPHelper::particle_truth particle = {eta,phi,ptotal,pid,is_scattered_lepton};
0045   return particle;
0046 }
0047
0048 float
0049 DVMPHelper::get_invariant_mass(DVMPHelper::particle_reco p1,
0050                    DVMPHelper::particle_reco p2)
0051 {
0052   float pt1 = DVMPHelper::get_pt(p1.eta,p1.ptotal);
0053   float pt2 = DVMPHelper::get_pt(p2.eta,p2.ptotal);
0054
0055   float eta1= p1.eta;
0056   float eta2= p2.eta;
0057
0058   float phi1= p1.phi;
0059   float phi2= p2.phi;
0060
0061   float Msquared = 2*pt1*pt2*(cosh(eta1-eta2)-cos(phi1-phi2));
0062   return sqrt(Msquared);
0063 }
0064
0065 float
0066 DVMPHelper::get_invariant_mass(DVMPHelper::particle_truth p1,
0067                    DVMPHelper::particle_truth p2)
0068 {
0069   float pt1 = DVMPHelper::get_pt(p1.eta,p1.ptotal);
0070   float pt2 = DVMPHelper::get_pt(p2.eta,p2.ptotal);
0071
0072   float eta1= p1.eta;
0073   float eta2= p2.eta;
0074
0075   float phi1= p1.phi;
0076   float phi2= p2.phi;
0077
0078   float Msquared = 2*pt1*pt2*(cosh(eta1-eta2)-cos(phi1-phi2));
0079   return sqrt(Msquared);
0080 }
0081
0082 float
0083 DVMPHelper::get_pt(float eta, float ptotal)
0084 {
0085   float theta = 2*atan(exp(-1*eta));
0086   float pt = ptotal * sin(theta);
0087   return pt;
0088 }
0089
0090 bool
0091 DVMPHelper::find_positron()
0092 {
0093   for(int i = 0 ; i < _size_reco ; i++)
0094     {
0095       if(rparticles[i].charge==1)
0096     return true;
0097     }
0098   return false;
0099 }
0100
0101 bool
0102 DVMPHelper::pass_cut(int index)
0103 {
0104   const float ep_cut = 0.7;
0105   if( rparticles[index].e / rparticles[index].ptotal > ep_cut )
0106     return true;
0107   else
0108     return false;
0109 }
0110
0111 std::vector<float>
0112 DVMPHelper::calculateInvariantMass_1()
0113 {
0114   std::vector<float> inv_mass;
0115
0116   // Assert number of particles in event is greater than 1
0117   if(_size_reco<=1)
0118     {
0119       inv_mass.push_back(NAN);
0120       return inv_mass;
0121     }
0122
0123   // Assert a positron exists
0124   if(!DVMPHelper::find_positron())
0125     {
0126       inv_mass.push_back(NAN);
0127       return inv_mass;
0128     }
0129
0130   // Record index of positron and electrons
0131   int idx_positron=-1;
0132   std::vector<int> idx_electron;
0133   for(int i = 0 ; i < _size_reco ; i++)
0134     {
0135       if( rparticles[i].charge == 1 )
0136     idx_positron=i;
0137       else if( rparticles[i].charge == -1 )
0138     idx_electron.push_back(i);
0139     }
0140   // Ensure positron and electron(s) pass ep_cut
0141   if(!(pass_cut(idx_positron)))
0142     {
0143       inv_mass.push_back(NAN);
0144       return inv_mass;
0145     }
0146
0147   for(unsigned i = 0 ; i < idx_electron.size() ; i++)
0148     {
0149       if( !pass_cut(idx_electron.at(i)) )
0150     idx_electron.erase(idx_electron.begin()+i);
0151     }
0152
0153   if(idx_electron.size()==0) // all found electrons fail cut
0154     {
0155       inv_mass.push_back(NAN);
0156       return inv_mass;
0157     }
0158   // Calculate invariant mass for all electron and positron pairs
0159   for(unsigned i = 0; i < idx_electron.size() ; i++)
0160     {
0161       inv_mass.push_back(DVMPHelper::get_invariant_mass( rparticles[idx_electron.at(i)], rparticles[idx_positron]));
0162     }
0163   return inv_mass;
0164 }
0165
0166 std::vector<float>
0167 DVMPHelper::calculateInvariantMass_2()
0168 {
0169   std::vector<float> inv_mass;
0170
0171   // Record index of positron and electrons
0172   int idx_positron=-1;
0173   std::vector<int> idx_electron;
0174   for(int i = 0 ; i < _size_truth ; i++)
0175     {
0176       if( tparticles[i].pid == -11 )
0177     idx_positron=i;
0178       else if( tparticles[i].pid == 11 )
0179     idx_electron.push_back(i);
0180     }
0181
0182   // Calculate invariant mass for all electron positron pairs
0183   for(unsigned i = 0; i < idx_electron.size() ; i++)
0184     {
0185       inv_mass.push_back(DVMPHelper::get_invariant_mass( tparticles[idx_electron.at(i)], tparticles[idx_positron]));
0186     }
0187   return inv_mass;
0188 }
0189
0190 std::vector<float>
0191 DVMPHelper::calculateInvariantMass_3()
0192 {
0193   std::vector<float> inv_mass;
0194
0195   // Assert number of particles in event is greater than 1
0196   if(_size_reco<=1)
0197     {
0198       inv_mass.push_back(NAN);
0199       return inv_mass;
0200     }
0201
0202   // Assert a positron exists
0203   if(!DVMPHelper::find_positron())
0204     {
0205       inv_mass.push_back(NAN);
0206       return inv_mass;
0207     }
0208
0209   // Record index of positron and electrons
0210   int idx_positron=-1;
0211   std::vector<int> idx_electron;
0212   for(int i = 0 ; i < _size_reco ; i++)
0213     {
0214       if( rparticles[i].charge == 1 )
0215     idx_positron=i;
0216       else if( rparticles[i].charge == -1 && rparticles[i].is_scattered_lepton==false )
0217     idx_electron.push_back(i);
0218     }
0219   // Ensure positron and electron(s) pass ep_cut
0220   if(!(pass_cut(idx_positron)))
0221     {
0222       inv_mass.push_back(NAN);
0223       return inv_mass;
0224     }
0225
0226   for(unsigned i = 0 ; i < idx_electron.size() ; i++)
0227     {
0228       if( !pass_cut(idx_electron.at(i)) )
0229     idx_electron.erase(idx_electron.begin()+i);
0230     }
0231
0232   if(idx_electron.size()==0) // all found electrons fail cut
0233     {
0234       inv_mass.push_back(NAN);
0235       return inv_mass;
0236     }
0237   // Calculate invariant mass for all electron and positron pairs
0238   for(unsigned i = 0; i < idx_electron.size() ; i++)
0239     {
0240       inv_mass.push_back(DVMPHelper::get_invariant_mass( rparticles[idx_electron.at(i)], rparticles[idx_positron]));
0241     }
0242   return inv_mass;
0243 }
0244
0245 std::vector<float>
0246 DVMPHelper::calculateInvariantMass_4()
0247 {
0248   std::vector<float> inv_mass;
0249
0250   // Assert number of particles in event is greater than 1
0251   if(_size_reco<=1)
0252     {
0253       inv_mass.push_back(NAN);
0254       return inv_mass;
0255     }
0256
0257   // Assert a positron exists
0258   if(!DVMPHelper::find_positron())
0259     {
0260       inv_mass.push_back(NAN);
0261       return inv_mass;
0262     }
0263
0264   // Record index of positron and electrons
0265   int idx_positron=-1;
0266   std::vector<int> idx_electron;
0267   for(int i = 0 ; i < _size_reco ; i++)
0268     {
0269       if( rparticles[i].charge == 1 )
0270     idx_positron=i;
0271       else if( rparticles[i].charge == -1 && rparticles[i].is_scattered_lepton==true )
0272     idx_electron.push_back(i);
0273     }
0274   // Ensure positron and electron(s) pass ep_cut
0275   if(!(pass_cut(idx_positron)))
0276     {
0277       inv_mass.push_back(NAN);
0278       return inv_mass;
0279     }
0280
0281   for(unsigned i = 0 ; i < idx_electron.size() ; i++)
0282     {
0283       if( !pass_cut(idx_electron.at(i)) )
0284     idx_electron.erase(idx_electron.begin()+i);
0285     }
0286
0287   if(idx_electron.size()==0) // all found electrons fail cut
0288     {
0289       inv_mass.push_back(NAN);
0290       return inv_mass;
0291     }
0292   // Calculate invariant mass for all electron and positron pairs
0293   for(unsigned i = 0; i < idx_electron.size() ; i++)
0294     {
0295       inv_mass.push_back(DVMPHelper::get_invariant_mass( rparticles[idx_electron.at(i)], rparticles[idx_positron]));
0296     }
0297   return inv_mass;
0298 }
0299
0300 std::vector<float>
0301 DVMPHelper::calculateInvariantMass_5()
0302 {
0303   std::vector<float> inv_mass;
0304
0305   // Record index of positron and electrons
0306   int idx_positron=-1;
0307   std::vector<int> idx_electron;
0308   for(int i = 0 ; i < _size_truth ; i++)
0309     {
0310       if( tparticles[i].pid == -11 )
0311     idx_positron=i;
0312       else if( tparticles[i].pid == 11 && tparticles[i].is_scattered_lepton==false)
0313     idx_electron.push_back(i);
0314     }
0315
0316   // Calculate invariant mass for all electron positron pairs
0317   for(unsigned i = 0; i < idx_electron.size() ; i++)
0318     {
0319       inv_mass.push_back(DVMPHelper::get_invariant_mass( tparticles[idx_electron.at(i)], tparticles[idx_positron]));
0320     }
0321
0322   return inv_mass;
0323 }
0324
0325 std::vector<float>
0326 DVMPHelper::calculateInvariantMass_6()
0327 {
0328   std::vector<float> inv_mass;
0329   // Record index of positron and electrons
0330   int idx_positron=-1;
0331   std::vector<int> idx_electron;
0332   for(int i = 0 ; i < _size_truth ; i++)
0333     {
0334       if( tparticles[i].pid == -11 )
0335     idx_positron=i;
0336       else if( tparticles[i].pid == 11 && tparticles[i].is_scattered_lepton==true)
0337     idx_electron.push_back(i);
0338     }
0339
0340   // Calculate invariant mass for all electron positron pairs
0341   for(unsigned i = 0; i < idx_electron.size() ; i++)
0342     {
0343       inv_mass.push_back(DVMPHelper::get_invariant_mass( tparticles[idx_electron.at(i)], tparticles[idx_positron]));
0344     }
0345   return inv_mass;
0346 }
0347