Neutron CrossTalk : first complete untested version

NPLib/Neutron/NPCrossTalk.cxx

@@ -58,19 +58,10 @@ void CrossTalk::AddHitVector(const vector<double>& X, const vector<double>& Y, c
 vector<double> CrossTalk::ComputeCrossTalk(){
 
   FirstHit = -1;
-  // take back the first hit, automaticaly identified as the first Neutron -> probably useless with this new method
-  for(int i = 0; i < sizeHit; i++){
-    static double minT;
-    minT = 1000000;
-    if((*HitT)[i] < minT){
-      minT = (*HitT)[i];
-      FirstHit = i;
-    }
-  }
 
   static double x1,y1,z1,dx1,dy1,dz1,t1;
   static double x2,y2,z2,dx2,dy2,dz2,t2;
-  static double Dist, dR1, dR2, v1;
+  static double Dist, dR1, dR2;
 
   // A different Cluster number (starting at 1) is assigned to each hit
   static vector<int> ID_ClustHit;
@@ -79,28 +70,6 @@ vector<double> CrossTalk::ComputeCrossTalk(){
     ID_ClustHit.push_back(i+1);
   }
   
-  /*
-  // Identify Hits around the 1st neutron, 1st cluster
-  x1 = (*HitX)[FirstHit], y1 = (*HitY)[FirstHit], z1 = (*HitZ)[FirstHit], dx1 = (*HitdX)[FirstHit], dy1 = (*HitdY)[FirstHit], dz1 = (*HitdZ)[FirstHit], t1 = (*HitT)[FirstHit];   
-  dR1 = sqrt(x1*dx1 + dy1*dy1 + dz1*dz1);
-  v1 = sqrt(x1*x1 + y1*y1 + z1*z1)/minT; 
-
-  for(int i = 0; i < sizeHit; i++){
-    if(i == FirstHit){
-      ID_ClustHit.push_back(1);
-      continue;
-    }
-    x2 = (*HitX)[i], y2 = (*HitY)[i], z2 = (*HitZ)[i], dx2 = (*HitdX)[i], dy2 = (*HitdX)[i], dz2 = (*HitdZ)[i], t2 = (*HitT)[i];   
-    Dist = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
-    if(Dist < coef*dR1){ //coef*dR1 : Radius of the cluster
-      ID_ClustHit.push_back(1); // this hit is part of the cluster 1
-    }
-    else(){
-      ID_ClustHit.push_back(-1);
-    }
-  }
-  */
-  
   //Test each Dist(n-n) to find clusters
   //When 2 neutrons are part of a the same cluster, change their ID_ClustHit for the lowest
   for(int j = 0; j < sizeHit-1; j++){
@@ -125,11 +94,12 @@ vector<double> CrossTalk::ComputeCrossTalk(){
   for(unsigned int i = 0; i < sizeHit; i++){
     mapOfClust[ID_ClustHit[i]].push_back(i);
   }
-  
+
   //Put first hit in time of each cluster in a new map mapOfFirstN
   //And also find the first hit of all : FirstN
   static map< unsigned int, unsigned int > mapOfFirstN;
   static unsigned int NbrOfClust, FirstN;
+  static int FirstClust;
   NbrOfClust = 0;
   NbrOfClust = mapOfClust.size();
   static double GeneralTmin; 
@@ -147,15 +117,50 @@ vector<double> CrossTalk::ComputeCrossTalk(){
       if((*HitT)[mapOfClust[i][j]] < GeneralTmin){
         GeneralTmin = (*HitT)[mapOfClust[i][j]];
         FirstN = mapOfClust[i][j];
+        FirstClust = i;
       }
     }
   }
    
   FirstHit = FirstN;
+  static double v_n, Dmax;
+  static vector<double> CrossTalk;
+  CrossTalk.clear();
+  //Test now the "Friend" clusters
+  for(int i = 1; i < NbrOfClust-1; i++){ 
+    x1 = (*HitX)[mapOfFirstN[i]], y1 = (*HitY)[mapOfFirstN[i]], z1 = (*HitZ)[mapOfFirstN[i]], t1 = (*HitT)[mapOfFirstN[i]]; 
+    v_n = sqrt(x1*x1+y1*y1+z1*z1)/t1;
+    for(int j = i+1; i < NbrOfClust; i++){
+      x2 = (*HitX)[mapOfFirstN[i]], y2 = (*HitY)[mapOfFirstN[i]], z2 = (*HitZ)[mapOfFirstN[i]], t2 = (*HitT)[mapOfFirstN[i]]; 
+      Dist = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
+      Dmax = (t2-t1)*v_n;
+      if(Dist < Dmax){
+        if(t1 < t2){
+          CrossTalk.push_back(j);
+        }
+        else{
+          CrossTalk.push_back(i);
+        }
+      }
+    }
+  }
 
-  static vector<double> m_Neutrons;
+  //elimate potential CrossTalk in mapOfFirstN
+  static unsigned int sizeCT;
+  sizeCT = CrossTalk.size();
+  for(unsigned int i = 0; i < sizeCT; i++  ){
+    mapOfFirstN.erase(i);
+  }
+  
+  //return vector of real neutrons IDs 
   m_Neutrons.clear();
+  for(auto itr = mapOfFirstN.begin(); itr != mapOfFirstN.end(); itr++){
+    m_Neutrons.push_back(itr->second);
+  }
+
   return m_Neutrons;
+
+
 }
 
 int CrossTalk::GetFirstN(){

NPLib/Neutron/NPCrossTalk.h

@@ -51,9 +51,11 @@ namespace NPL{
       const std::vector<double>* HitT;
 
       std::vector<int> ClustHit;
+      std::vector<double> m_Neutrons;
       double coef;
       unsigned int sizeHit ;
-      int FirstHit ;     
+      int FirstHit;
+      
   };
 }