Wide angle scattering events degree project

Recon/include/Recon/ParticleFlow.h

@@ -61,6 +61,11 @@ class ParticleFlow : public framework::Producer {
   std::string outputCollName_;
   // configuration
   bool singleParticle_;
+
+  // configurable matching criteria for Track + (side) HCal cluster matching
+  double tkHadCaloMatchDist_;
+  double tkHadCaloMinEnergyRatio_;
+  double tkHadCaloMaxEnergyRatio_;
 };
 }  // namespace recon
 

Recon/python/pfReco.py

@@ -56,6 +56,11 @@ def __init__(self, name='PFlow') :
         self.inputTrackCollName = 'PFTracks'
         self.outputCollName     = 'PFCandidates'
         self.singleParticle     = False
+
+        # Matching criteria for Track + (side) HCal cluster matching (need to be configured)
+        self.tkHadCaloMatchDist         = 2.0
+        self.tkHadCaloMinEnergyRatio    = 0.3
+        self.tkHadCaloMaxEnergyRatio    = 2.0
 
         self.input_ecal_passname  = ''
         self.input_hcal_passname  = ''
@@ -79,4 +84,4 @@ def __init__(self, name='PFTruth') :
         self.sim_particles_event_passname = ''
         self.ecal_sp_hits_event_passname = ''
         self.target_sp_hits_event_passname = ''
-        self.target_sp_passname = ''
+        self.target_sp_passname = ''

Recon/src/Recon/ParticleFlow.cxx

@@ -25,6 +25,9 @@ void ParticleFlow::configure(framework::config::Parameters& ps) {
 
   // Algorithm configuration
   singleParticle_ = ps.getParameter<bool>("singleParticle");
+  tkHadCaloMatchDist_ = ps.getParameter<double>("tkHadCaloMatchDist");
+  tkHadCaloMinEnergyRatio_ = ps.getParameter<double>("tkHadCaloMinEnergyRatio");
+  tkHadCaloMaxEnergyRatio_ = ps.getParameter<double>("tkHadCaloMaxEnergyRatio");
 
   // Calibration factors, from jason, temperary
   std::vector<float> em1{250.0,  750.0,  1250.0, 1750.0, 2250.0, 2750.0,
@@ -72,7 +75,7 @@ void ParticleFlow::fillCandEMCalo(ldmx::PFCandidate& cand,
     corr = eCorr_->GetY()[0];
   } else if (e > eCorr_->GetX()[eCorr_->GetN() - 1]) {
     corr = eCorr_->GetY()[eCorr_->GetN() - 1];
-  } else {  // else look up calibration factor
+  } else {
     corr = eCorr_->Eval(e);
   }
   cand.setEcalEnergy(e * corr);
@@ -135,7 +138,7 @@ void ParticleFlow::produce(framework::Event& event) {
       4c. (Upstream?) Categorize hcal clusters as: EM/Had-like
       5. Build candidates by category, moving from Tk-Ecal-Hcal
     */
-
+    std::cout << "PF Multiple Particles for Event: "<< event.getEventNumber() << std::endl;
     //
     // track-calo linking
     //
@@ -209,21 +212,43 @@ void ParticleFlow::produce(framework::Event& event) {
       }
     }
 
-    // NOT YET IMPLEMENTED...
-    // tk-hadcalo linking (Side HCal)
+    // tk-hadcalo linking (Side HCal) BY TILLRUE
     std::map<int, std::vector<int> > tkHadCaloMap;
-    // for(int i=0; i<tracks.size(); i++){
-    //   const auto& tk = tracks[i];
-    //   for(int j=0; j<hcalClusters.size(); j++){
-    // 	const auto& hcal = hcalClusters[j];
-    // 	// TODO: add the matching logic here...
-    // 	bool isMatch = true;
-    // 	if(isMatch){
-    // 	  if (tkHadCaloMap.count(i)) tkHadCaloMap[i].push_back(j);
-    // 	  else tkHadCaloMap[i] = {j};
-    // 	}
-    //   }
-    // }
+    std::map<int, std::vector<int> > HadCaloTkMap;
+    std::map<std::pair<int, int>, float> tkHadCaloDist;
+    for(int i=0; i<tracks.size(); i++){
+      const auto& tk = tracks[i];
+      const std::vector<float> xyz = tk.getPosition();
+      const std::vector<double> pxyz = tk.getMomentum();
+      const float p = sqrt(pow(pxyz[0], 2) + pow(pxyz[1], 2) + pow(pxyz[2], 2));
+
+      for(int j=0; j<hcalClusters.size(); j++){
+       	const auto& hcal = hcalClusters[j];
+        // Matching logic same as for track ecalCluster matching
+        const float hcalClusZ = hcal.getCentroidZ();
+        const float tkXAtClus =
+            xyz[0] + pxyz[0] / pxyz[2] * (hcalClusZ - xyz[2]);  // extrapolation
+        const float tkYAtClus =
+            xyz[1] + pxyz[1] / pxyz[2] * (hcalClusZ - xyz[2]);
+        float dist = hypot(
+            (tkXAtClus - hcal.getCentroidX()) / std::max(1.0, hcal.getRMSX()),
+            (tkYAtClus - hcal.getCentroidY()) / std::max(1.0, hcal.getRMSY()));
+        tkHadCaloDist[{i, j}] = dist;
+        bool isMatch =
+            (dist < tkHadCaloMatchDist_) && (hcal.getEnergy() > tkHadCaloMinEnergyRatio_ * p &&
+                           hcal.getEnergy() < tkHadCaloMaxEnergyRatio_ * p);  // matching criteria, need to be configured
+         if (isMatch) {
+          if (tkHadCaloMap.count(i))
+            tkHadCaloMap[i].push_back(j);
+          else
+            tkHadCaloMap[i] = {j};
+          if (HadCaloTkMap.count(j))
+            HadCaloTkMap[j].push_back(i);
+          else
+            HadCaloTkMap[j] = {i};
+        }
+      }
+    }
 
     //
     // track / ecal cluster arbitration
@@ -245,10 +270,12 @@ void ParticleFlow::produce(framework::Event& event) {
       }
     }
 
-    // track / hcal cluster arbitration
+    // ecal / hcal cluster arbitration
     std::vector<bool> EMIsHadLinked(ecalClusters.size(), false);
     std::vector<bool> HadIsEMLinked(hcalClusters.size(), false);
     std::map<int, int> EMHadPairs{};
+    std::map<int, int> HadEMPairs{};
+
     for (int i = 0; i < ecalClusters.size(); i++) {
       if (emHadCaloMap.count(i)) {
         // pick first (highest-energy) unused matching cluster
@@ -257,6 +284,27 @@ void ParticleFlow::produce(framework::Event& event) {
             HadIsEMLinked[had_idx] = true;
             EMIsHadLinked[i] = true;
             EMHadPairs[i] = had_idx;
+            HadEMPairs[had_idx] = i;
+            break;
+          }
+        }
+      }
+    }
+
+    //
+    // track / hcal cluster arbitration BY TILLRUE
+    //
+    std::vector<bool> tkIsHadLinked(tracks.size(), false);
+    std::vector<bool> HadIsTkLinked(hcalClusters.size(), false);
+    std::map<int, int> tkHadPairs{};
+    for (int i = 0; i < tracks.size(); i++) {
+      if (tkHadCaloMap.count(i)) {
+        // pick first (highest-energy) unused matching cluster
+        for (int had_idx : tkHadCaloMap[i]) {
+          if (!HadIsTkLinked[had_idx]) {
+            HadIsTkLinked[had_idx] = true;
+            tkIsHadLinked[i] = true;
+            tkHadPairs[i] = had_idx;
             break;
           }
         }
@@ -271,49 +319,62 @@ void ParticleFlow::produce(framework::Event& event) {
     // Begin building pf candidates from tracks
     //
 
-    // std::vector<ldmx::PFCandidate> chargedMatch;
-    // std::vector<ldmx::PFCandidate> chargedUnmatch;
+    // loop through tracks
     for (int i = 0; i < tracks.size(); i++) {
       ldmx::PFCandidate cand;
       fillCandTrack(cand, tracks[i]);  // append track info to candidate
-
-      if (!tkIsEMLinked[i]) {
-        // chargedUnmatch.push_back(cand);
-      } else {  // if track is linked with ECal cluster
+      if (tkIsEMLinked[i]){ // if track is linked with ECal cluster
         fillCandEMCalo(cand, ecalClusters[tkEMPairs[i]]);
-        if (EMIsHadLinked[tkEMPairs[i]]) {  // if ECal is linked with HCal
-                                            // cluster
-          fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]);
+        if (EMIsHadLinked[tkEMPairs[i]]) {  // if ECal is linked with HCal cluster
+          if (!tkIsHadLinked[i]){   // if track is NOT also linked with Hcal
+          fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]);}
+        }
+      }
+
+      if(tkIsHadLinked[i]){ // if track is linked with Hcal
+        fillCandHadCalo(cand, hcalClusters[tkHadPairs[i]]);
+        if(HadIsEMLinked[tkHadPairs[i]]){ // if hcal is linked with ecal
+          if(!tkIsEMLinked[i]){ // if ecal is NOT also linked with track
+          fillCandEMCalo(cand,ecalClusters[HadEMPairs[tkHadPairs[i]]]);
+          }
+        } else if(!tkIsEMLinked[i]){
+          cand.setDistTkHcalMatch(tkHadCaloDist[{i,tkHadPairs[i]}]);
+          if(hcalClusters.size() > 1){
+            cand.setHcalSecondEnergy(hcalClusters[tkHadPairs[i]+1].getEnergy());
+            cand.setDistTkHcalSecondCluster(tkHadCaloDist[{i,tkHadPairs[i]+1}]);
+          }
         }
-        // chargedMatch.push_back(cand);
       }
       pfCands.push_back(cand);
     }
 
     // std::vector<ldmx::PFCandidate> emMatch;
     // std::vector<ldmx::PFCandidate> emUnmatch;
     for (int i = 0; i < ecalClusters.size(); i++) {
-      // already linked with ECal in the previous step
-      if (EMIsTkLinked[i]) continue;
+      if (EMIsTkLinked[i]) continue; // already linked with ECal with Track in the previous step
+      if (EMIsHadLinked[i] && HadIsTkLinked[EMHadPairs[i]]) continue; // already linked with track through Hcal
 
       ldmx::PFCandidate cand;
       fillCandEMCalo(cand, ecalClusters[i]);
-      if (EMIsHadLinked[tkEMPairs[i]]) {
+      if (EMIsHadLinked[i]) { // is Ecal is linked with Hcal
         fillCandHadCalo(cand, hcalClusters[EMHadPairs[i]]);
         // emMatch.push_back(cand);
       } else {
         // emUnmatch.push_back(cand);
       }
       pfCands.push_back(cand);
     }
-    std::vector<ldmx::PFCandidate> hadOnly;
+
+    // std::vector<ldmx::PFCandidate> hadOnly;
     for (int i = 0; i < hcalClusters.size(); i++) {
-      if (HadIsEMLinked[i]) continue;
+      if (HadIsEMLinked[i]) continue; // already linked with ecal
+      if (HadIsTkLinked[i]) continue; // already linked with track
       ldmx::PFCandidate cand;
       fillCandHadCalo(cand, hcalClusters[i]);
       // hadOnly.push_back(cand);
       pfCands.push_back(cand);
     }
+
 
     // // track / ecal cluster arbitration
     // std::vector<ldmx::PFCandidate> caloMatchedTks;
@@ -386,6 +447,7 @@ void ParticleFlow::produce(framework::Event& event) {
     // }
 
   } else {
+    std::cout << "Matching Single Particle Signals" << std::endl;
     // Single-particle builder
     ldmx::PFCandidate pf;
     int pid = 0;  // initialize pid to add

SimCore/python/generators.py

@@ -359,3 +359,32 @@ def single_backwards_positron(energy: float):
     beam.energy = energy
     return beam
 
+def single_e_wide_angle_downstream_target(minTheta = 30, maxTheta = 70, minPhi = 0, maxPhi = 360):
+    """A general particle source configured to shoot electrons downstream from after the target at wide angles such that they hit the side hcal.
+
+    This generator is helpful to study the matching criteria between tracks and signals in the (side) hcal. Theta is the angle with respect to the downstream z axis. Note that the angular distribution and energy of the beam electrons can easily be modified. The particle source starts behind the target.
+
+    Returns
+    -------
+    gun:
+        configured general particle source to shoot electrons downstream the target at wide angles
+    """
+    myGPS = gps( 'myGPS' , [
+            "/gps/particle e-",
+            "/gps/number 1",
+            "/gps/pos/type Plane",
+            "/gps/pos/shape Rectangle",
+            "/gps/pos/centre 0 0 0 mm",
+            "/gps/pos/halfx 10 mm",
+            "/gps/pos/halfy 40 mm",
+            "/gps/ang/type iso",
+            f"/gps/ang/mintheta {minTheta} deg",
+            f"/gps/ang/maxtheta {maxTheta} deg",
+            f"/gps/ang/minphi {minPhi} deg",
+            f"/gps/ang/maxphi {maxPhi} deg",
+            f"/gps/ang/rot1 0 1 0", # These have been determined by trial and error so that theta is the angle from the ldmx z axis and
+            f"/gps/ang/rot2 1 0 0", # phi is rotates clockwise in the XY plane starting from the negative Y axis
+            "/gps/ene/mono 8 GeV",
+            ] )
+    return myGPS
+