Elnura code updates for ALPs

.gitignore

@@ -112,3 +112,8 @@ benchmark*json
 # Graphviz graphs
 *.dot
 *.png
+
+# Condor logs
+*.err
+*.out
+*.log

analyzers/dataframe/FCCAnalyses/MCParticle.h

@@ -39,6 +39,13 @@ namespace MCParticle{
     ROOT::VecOps::RVec<edm4hep::MCParticleData>  operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
   };
 
+  /// select MCParticles with absolute pseudorapidity less than a max value
+  struct sel_eta {
+    sel_eta(float arg_min_eta);
+    float m_min_eta = 2.5; //> pseudorapidity threshold
+    ROOT::VecOps::RVec<edm4hep::MCParticleData>  operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
+  };
+
   /// select MCParticles with their status
   struct sel_genStatus {
     sel_genStatus(int arg_status);
@@ -188,9 +195,18 @@ namespace MCParticle{
   /// return the phi of the input MCParticles
   ROOT::VecOps::RVec<float> get_phi(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
 
-  /// return the delta r of the input ReconstructedParticles
+  /// return the delta eta of the input MCParticles
+  ROOT::VecOps::RVec<float> get_delta_eta(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
+
+  /// return the delta phi of the input MCParticles
+  ROOT::VecOps::RVec<float> get_delta_phi(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
+
+  /// return the delta r of the input MCParticles
   ROOT::VecOps::RVec<float> get_delta_r(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
 
+   /// return the min delta r of the input MCParticles
+  ROOT::VecOps::RVec<float> get_min_delta_r(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
+
   /// return the energy of the input MCParticles
   ROOT::VecOps::RVec<float> get_e(ROOT::VecOps::RVec<edm4hep::MCParticleData> in);
 

analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h

@@ -110,6 +110,8 @@ namespace ReconstructedParticle{
   /// return the momenta of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_p(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
+  float get_diphoton_deltaR(edm4hep::ReconstructedParticleData p, edm4hep::ReconstructedParticleData diphoton1 , edm4hep::ReconstructedParticleData diphoton2);
+
   /// return the momenta of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_px(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
@@ -128,6 +130,14 @@ namespace ReconstructedParticle{
   /// return the theta of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_theta(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
+  /// return the distance of particle to the calorimeter barrel or endcap based on its current trajectory
+  //ROOT::VecOps::RVec<float> get_L2calorimeter(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> photon, float ALP_Lx, float ALP_Ly, float ALP_Lz);
+
+  /// return coordiantes of calorimeter hit of a prompt photon coming from the interaction point
+  ROOT::VecOps::RVec<std::vector<float>> get_prompt_photon_calorimeter_hits(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> photon);
+
+  ROOT::VecOps::RVec<std::vector<float>> get_displaced_photon_calorimeter_hits(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> photon,float ALP_Lx, float ALP_Ly, float ALP_Lz);
+
   /// return the phi of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_phi(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
@@ -140,9 +150,15 @@ namespace ReconstructedParticle{
   /// return the delta r of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
+    /// return the delta r of the input ReconstructedParticles
+  float get_delta_r(edm4hep::ReconstructedParticleData p1, edm4hep::ReconstructedParticleData p2);
+
   /// return the minimum delta r of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_min_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 
+  /// return the photon indices of the minimum delta r of the input ReconstructedParticles
+  ROOT::VecOps::RVec<float> get_pidx_min_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
+
   /// return the energy of the input ReconstructedParticles
   ROOT::VecOps::RVec<float> get_reference_point_x(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in);
 

analyzers/dataframe/src/MCParticle.cc

@@ -72,6 +72,21 @@ ROOT::VecOps::RVec<edm4hep::MCParticleData>  sel_pt::operator() (ROOT::VecOps::R
   return result;
 }
 
+sel_eta::sel_eta(float arg_min_eta) : m_min_eta(arg_min_eta) {};
+ROOT::VecOps::RVec<edm4hep::MCParticleData>  sel_eta::operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
+  ROOT::VecOps::RVec<edm4hep::MCParticleData> result;
+  result.reserve(in.size());
+  for (size_t i = 0; i < in.size(); ++i) {
+    auto & p = in[i];
+    TLorentzVector tv1;
+    tv1.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass);
+    if (abs(tv1.Eta()) < abs(m_min_eta)){
+      result.emplace_back(p);
+    }
+  }
+  return result;
+}
+
 
 filter_pdgID::filter_pdgID(int arg_pdgid, bool arg_abs){m_pdgid = arg_pdgid; m_abs = arg_abs;};
 bool  filter_pdgID::operator() (ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
@@ -337,6 +352,38 @@ ROOT::VecOps::RVec<float> get_phi(ROOT::VecOps::RVec<edm4hep::MCParticleData> in
   return result;
 }
 
+
+ROOT::VecOps::RVec<float> get_delta_eta(ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  for (int i = 0; i < in.size(); i++) {
+    TLorentzVector tlv1;
+    tlv1.SetXYZM(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+    for (auto j = i + 1; j < in.size(); j++) {
+      TLorentzVector tlv2;
+      tlv2.SetXYZM(in[j].momentum.x, in[j].momentum.y, in[j].momentum.z, in[j].mass);
+      float delta_eta = abs(tlv1.Eta() - tlv2.Eta());
+      result.push_back(delta_eta);
+    }
+  }
+  return result;
+}
+
+ROOT::VecOps::RVec<float> get_delta_phi(ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  for (int i = 0; i < in.size(); i++) {
+    TLorentzVector tlv1;
+    tlv1.SetXYZM(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+    for (auto j = i + 1; j < in.size(); j++) {
+      TLorentzVector tlv2;
+      tlv2.SetXYZM(in[j].momentum.x, in[j].momentum.y, in[j].momentum.z, in[j].mass);
+      float delta_phi = tlv1.DeltaPhi(tlv2);
+      result.push_back(delta_phi);
+    }
+  }
+  return result;
+}
+
+
 ROOT::VecOps::RVec<float> get_delta_r(ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   for (int i = 0; i < in.size(); i++) {
@@ -352,6 +399,32 @@ ROOT::VecOps::RVec<float> get_delta_r(ROOT::VecOps::RVec<edm4hep::MCParticleData
   return result;
 }
 
+ROOT::VecOps::RVec<float> get_min_delta_r(ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  float min_delta_r = 999;
+  // p1idx = -1
+  // p2idx = -1
+  for (int i = 0; i < in.size(); i++) {
+    TLorentzVector tlv1;
+    tlv1.SetXYZM(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+    for (int j = i + 1; j < in.size(); j++) {
+      TLorentzVector tlv2;
+      tlv2.SetXYZM(in[j].momentum.x, in[j].momentum.y, in[j].momentum.z, in[j].mass);
+      float delta_r = tlv1.DeltaR(tlv2);
+      if (delta_r < min_delta_r) {
+        min_delta_r = delta_r;
+        // p1idx = i
+        // p2idx = j
+      }
+    }
+  }
+  result.push_back(min_delta_r);
+  // if(p1idx > -1)
+    // result.emplace_back(in[p1idx]);
+    // result.emplace_back(in[p2idx]);
+  return result;
+}
+
 ROOT::VecOps::RVec<float> get_e(ROOT::VecOps::RVec<edm4hep::MCParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   for (auto & p: in) {

analyzers/dataframe/src/ReconstructedParticle.cc

@@ -387,6 +387,15 @@ ROOT::VecOps::RVec<float> get_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedP
   return result;
 }
 
+float get_delta_r(edm4hep::ReconstructedParticleData p1, edm4hep::ReconstructedParticleData p2) {
+  TLorentzVector tlv1;
+  tlv1.SetXYZM(p1.momentum.x, p1.momentum.y, p1.momentum.z, p1.mass);
+  TLorentzVector tlv2;
+  tlv2.SetXYZM(p2.momentum.x, p2.momentum.y, p2.momentum.z, p2.mass);
+  float delta_r = tlv1.DeltaR(tlv2);
+  return delta_r;
+}
+
 ROOT::VecOps::RVec<float> get_min_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   float min_delta_r = 999;
@@ -406,6 +415,35 @@ ROOT::VecOps::RVec<float> get_min_delta_r(ROOT::VecOps::RVec<edm4hep::Reconstruc
   return result;
 }
 
+ROOT::VecOps::RVec<float> get_pidx_min_delta_r(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  float min_delta_r = 999;
+  int p1idx = -1;
+  int p2idx = -1;
+  for (int i = 0; i < in.size(); i++) {
+    TLorentzVector tlv1;
+    tlv1.SetXYZM(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+    for (int j = i + 1; j < in.size(); j++) {
+      TLorentzVector tlv2;
+      tlv2.SetXYZM(in[j].momentum.x, in[j].momentum.y, in[j].momentum.z, in[j].mass);
+      float delta_r = tlv1.DeltaR(tlv2);
+      if (delta_r < min_delta_r) {
+        min_delta_r = delta_r;
+        p1idx = i;
+        p2idx = j;
+      }
+    }
+  }
+
+  if(p1idx > -1 && p2idx > -1) 
+    result.emplace_back(p1idx);
+    result.emplace_back(p2idx);
+    // std::cout << p1idx << " " << p2idx << std::endl;
+  return result;
+}
+
+
+
 ROOT::VecOps::RVec<float> get_reference_point_x(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   for (auto & p: in) {
@@ -414,6 +452,22 @@ ROOT::VecOps::RVec<float> get_reference_point_x(ROOT::VecOps::RVec<edm4hep::Reco
   return result;
 }
 
+ROOT::VecOps::RVec<float> get_reference_point_y(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  for (auto & p: in) {
+    result.push_back(p.referencePoint.y);
+  }
+  return result;
+}
+
+ROOT::VecOps::RVec<float> get_reference_point_z(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
+  ROOT::VecOps::RVec<float> result;
+  for (auto & p: in) {
+    result.push_back(p.referencePoint.z);
+  }
+  return result;
+}
+
 ROOT::VecOps::RVec<float> get_e(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   for (auto & p: in) {
@@ -440,6 +494,15 @@ ROOT::VecOps::RVec<float> get_p(ROOT::VecOps::RVec<edm4hep::ReconstructedParticl
   return result;
 }
 
+float get_diphoton_deltaR(edm4hep::ReconstructedParticleData p, edm4hep::ReconstructedParticleData diphoton1 , edm4hep::ReconstructedParticleData diphoton2) {
+  ROOT::VecOps::RVec<float> result;
+  TLorentzVector diphoton_tlv, p_tlv;
+  diphoton_tlv.SetXYZM(diphoton1.momentum.x+diphoton2.momentum.x, diphoton1.momentum.y+diphoton2.momentum.y, diphoton1.momentum.z+diphoton2.momentum.z, diphoton1.mass+diphoton2.mass);
+  p_tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass);
+  float delta_r = diphoton_tlv.DeltaR(p_tlv);
+  return delta_r;
+}
+
 ROOT::VecOps::RVec<float> get_px(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
   ROOT::VecOps::RVec<float> result;
   for (auto & p: in) {
@@ -493,6 +556,82 @@ ROOT::VecOps::RVec<float> get_theta(ROOT::VecOps::RVec<edm4hep::ReconstructedPar
   return result;
 }
 
+
+ROOT::VecOps::RVec<std::vector<float>> get_prompt_photon_calorimeter_hits(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> photon) {
+  const float R = 2500.0; //calorimeter radius and length is 2500mm
+  const float Zmax = 2500.0;
+  ROOT::VecOps::RVec<std::vector<float>> result;
+  for (auto & p : photon) {
+    TLorentzVector tlv;
+    tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass);
+    TVector3 direction = tlv.Vect().Unit(); 
+    float ux = static_cast<float>(direction.X());
+    float uy = static_cast<float>(direction.Y());  
+    float uz = static_cast<float>(direction.Z());  
+    float t_r = R / sqrt(ux*ux + uy*uy);  
+    float t_z = Zmax / fabs(uz);             
+    float t_hit = std::min(t_r, t_z); 
+
+    // float ux2 = sin(tlv.Theta())*cos(tlv.Phi());
+    // float uy2 = sin(tlv.Theta())*sin(tlv.Phi());
+    // float uz2 = cos(tlv.Theta());
+
+    std::vector<float> hit_pos = {t_hit * ux, t_hit * uy, t_hit * uz};
+
+
+
+    result.push_back(hit_pos);
+  }
+  return result;
+}
+
+ROOT::VecOps::RVec<std::vector<float>> get_displaced_photon_calorimeter_hits(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> photon,float ALP_Lx, float ALP_Ly, float ALP_Lz) {
+  const float R = 2500.0; //calorimeter radius and length is 2500mm
+  const float Zmax = 2500.0;
+  ROOT::VecOps::RVec<std::vector<float>> result;
+  for (auto & p : photon) {
+    TLorentzVector tlv;
+    tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass);
+    TVector3 direction = tlv.Vect().Unit(); 
+    float ux = static_cast<float>(direction.X());
+    float uy = static_cast<float>(direction.Y());  
+    float uz = static_cast<float>(direction.Z());  
+
+    ///solving quadratic equation to find intersection with calorimeter
+    /// (ALP_Lx + t * ux)^2 + (ALP_Ly + t * uy)^2 = R^2 and solve for t
+    /// (ux^2 + uy^2) * t^2 + 2 * (ALP_Lx * ux + ALP_Ly * uy) * t + (ALP_Lx^2 + ALP_Ly^2 - R^2) = 0
+    float a = 2.0f * (ALP_Lx * ux + ALP_Ly * uy);
+    float b = a * a - 4 * (ux * ux + uy * uy) * (ALP_Lx * ALP_Lx + ALP_Ly * ALP_Ly - R * R);
+
+    float t_r = 99999999;
+
+    if (b >= 0) {
+      float t1 = (-a + sqrt(b)) / (2 * (ux * ux + uy * uy));
+      float t2 = (-a - sqrt(b)) / (2 * (ux * ux + uy * uy));
+      if (t1 > 0) t_r = t1;
+      if (t2 > 0) t_r = std::min(t_r, t2);
+    }
+
+    float t_z_pos = (Zmax - ALP_Lz) / uz;
+    float t_z_neg = (-Zmax - ALP_Lz) / uz;
+    float t_z = 99999999;
+    if (t_z_pos > 0 && fabs(ALP_Lz + t_z_pos * uz) <= Zmax) t_z = t_z_pos;
+    if (t_z_neg > 0 && fabs(ALP_Lz + t_z_neg * uz) <= Zmax) t_z = std::min(t_z, t_z_neg);
+
+    if (t_r < 0) t_r = 99999999; 
+    if (t_z < 0) t_z = 99999999; 
+
+    float t_hit = std::min(t_r, t_z); 
+
+    std::vector<float> hit_pos = {ALP_Lx + t_hit * ux, ALP_Ly + t_hit * uy, ALP_Lz + t_hit * uz};
+    result.push_back(hit_pos);
+  }
+  return result;
+}
+
+
+
+
 ROOT::VecOps::RVec<TLorentzVector> get_tlv(ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> in) {
   ROOT::VecOps::RVec<TLorentzVector> result;
   for (auto & p: in) {

examples/FCCee/bsm/LLPs/ALPs/ALP_sample_creation/ALP_pythia.cmnd

@@ -1,17 +1,19 @@
 ! File: ALP_pythia.cmd
 Random:setSeed = on
 Main:timesAllowErrors = 10         ! how many aborts before run stops
-Main:numberOfEvents = 100000         ! number of events to generate
+Main:numberOfEvents = 1000000         ! number of events to generate
 
 
 ! 2) Settings related to output in init(), next() and stat().
-Next:numberCount = 100             ! print message every n events
+Next:numberCount = 1000             ! print message every n events
 !Beams:idA = 11                     ! first beam, e+ = 11
 !Beams:idB = -11                    ! second beam, e- = -11
 
 Beams:frameType = 4                ! read info from a LHEF
 ! Change the LHE file here
-Beams:LHEF = ALP_Z_aa_1GeV_cYY_0p5.lhe
+!Beams:LHEF = /eos/experiment/fcc/ee/analyses_storage/BSM/LLPs/ALPs_3photons/mg5amcnlo/madgraph_5_1.0/Events/run_01/unweighted_events.lhe
+Beams:LHEF = /eos/user/e/ebakhish/MG/MGOutput/Backgrounds/background_ee_aa/Events/run_01/unweighted_events.lhe.gz
+
 
 ! 3) Settings for the event generation process in the Pythia8 library.
 PartonLevel:ISR = on               ! initial-state radiation