From 674847dcdaeb0aab2c647205d807c2da2f685412 Mon Sep 17 00:00:00 2001
From: moukaddam <mhd.moukaddam@gmail.com>
Date: Mon, 13 Feb 2017 10:40:47 +0000
Subject: [PATCH] coding new add-back method based on segment-core pixel

---
 NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx | 132 ++++++++++++++--------
 NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h   |  17 +--
 2 files changed, 93 insertions(+), 56 deletions(-)

diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
index 8618c4293..9047579bf 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
@@ -53,6 +53,7 @@ ClassImp(TGeTAMUPhysics)
 void TGeTAMUPhysics::BuildPhysicalEvent(){
   PreTreat();
 
+/*
  //Treat singles
   for(unsigned int iPixel = 0 ; iPixel < Singles_E.size() ; iPixel++){
     int clv = Singles_Clover[iPixel];
@@ -82,46 +83,44 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
   map<int,TVector3> max_pos; // first hit
   map<int,int> cry_segment;    
 
-/*
-
-
-  for(unsigned int i = 0 ; i < c_size_e ; i++){
-    int clv = m_PreTreatedData->GetCoreCloverNbrE(i);
-    int cry = m_PreTreatedData->GetCoreCrystalNbrE(i);
-    double energy = m_PreTreatedData->GetCoreEnergy(i);
-    // Add back energy
-    clv_energy[clv] += energy;
-    // Pick up the crystal with the maximum energy in every clover 
-    if(energy > max_cry[clv]){
-      max_cry[clv] = energy;
-      clv_crystal[clv] = cry;
-    }
-    // Pick up the segment with the maximum energy in every clover
-    for(unsigned int j = 0 ; j < s_size_e ; j++){
-      double s_energy = m_PreTreatedData->GetSegmentEnergy(j); 
-      if(s_energy > max_segment[clv]){
-        max_segment[clv] = s_energy;
-        clv_segment[clv] = m_PreTreatedData->GetSegmentSegmentNbrE(j);
-      }
-    }
-  }
-
-    // Pick up the segment with the maximum energy in every clover
-  for(unsigned int i = 0 ; i < c_size_e ; i++){
-    int clv = m_PreTreatedData->GetCoreCloverNbrE(i);
-    int cry = m_PreTreatedData->GetCoreCrystalNbrE(i);
-    double energy = m_PreTreatedData->GetCoreEnergy(i);
-    }
+  
 
-    for(unsigned int j = 0 ; j < s_size_e ; j++){
-      double s_energy = m_PreTreatedData->GetSegmentEnergy(j); 
-      if(s_energy > max_segment[clv]){
-        max_segment[clv] = s_energy;
-        clv_segment[clv] = m_PreTreatedData->GetSegmentSegmentNbrE(j);
-      }
-    }
 
-*/
+                    for(unsigned int i = 0 ; i < c_size_e ; i++){
+                      int clv = m_PreTreatedData->GetCoreCloverNbrE(i);
+                      int cry = m_PreTreatedData->GetCoreCrystalNbrE(i);
+                      double energy = m_PreTreatedData->GetCoreEnergy(i);
+                      // Add back energy
+                      clv_energy[clv] += energy;
+                      // Pick up the crystal with the maximum energy in every clover 
+                      if(energy > max_cry[clv]){
+                        max_cry[clv] = energy;
+                        clv_crystal[clv] = cry;
+                      }
+                      // Pick up the segment with the maximum energy in every clover
+                      for(unsigned int j = 0 ; j < s_size_e ; j++){
+                        double s_energy = m_PreTreatedData->GetSegmentEnergy(j); 
+                        if(s_energy > max_segment[clv]){
+                          max_segment[clv] = s_energy;
+                          clv_segment[clv] = m_PreTreatedData->GetSegmentSegmentNbrE(j);
+                        }
+                      }
+                    }
+
+                      // Pick up the segment with the maximum energy in every clover
+                    for(unsigned int i = 0 ; i < c_size_e ; i++){
+                      int clv = m_PreTreatedData->GetCoreCloverNbrE(i);
+                      int cry = m_PreTreatedData->GetCoreCrystalNbrE(i);
+                      double energy = m_PreTreatedData->GetCoreEnergy(i);
+                      }
+
+                      for(unsigned int j = 0 ; j < s_size_e ; j++){
+                        double s_energy = m_PreTreatedData->GetSegmentEnergy(j); 
+                        if(s_energy > max_segment[clv]){
+                          max_segment[clv] = s_energy;
+                          clv_segment[clv] = m_PreTreatedData->GetSegmentSegmentNbrE(j);
+                        }
+                      }
 
 
   for(unsigned int i = 0 ; i < c_size_e ; i++){
@@ -170,6 +169,7 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
     AddBack_Z.push_back(-1000);
   }
 
+*/
 //Fill the time OR
 for (unsigned i = 0 ; i < m_PreTreatedData->GetMultiplicityCoreT(); i++)
   GeTime.push_back(m_PreTreatedData->GetCoreTime(i));
@@ -388,6 +388,52 @@ double TGeTAMUPhysics::GetDopplerCorrectedEnergy(double& energy , TVector3 posit
   return m_GammaLV.Energy();
 }
 
+void TGeTAMUPhysics::FillAddBack(int scheme, TVector3& beta){
+
+    vector<int>::iterator itClover;
+
+  if (scheme==1){
+   //cout << " Clover Add-Back, singles: " << Singles_E.size()<< endl; 
+   //Treat singles
+    for(unsigned int iPixel = 0 ; iPixel < Singles_E.size() ; iPixel++){
+      int clv = Singles_Clover[iPixel];
+      int cry = Singles_Crystal[iPixel];
+      int seg = Singles_Segment[iPixel];
+      double energy = Singles_E[iPixel];
+      //cout << clv << " " << cry << " " << seg << " "<< energy << endl; 
+      itClover = find (AddBack_Clover.begin(), AddBack_Clover.end(), clv); 
+      bool end = (itClover == AddBack_Clover.end());
+      if ( end ){ // Clover is not found
+        // Fill these values only for the first hit
+        AddBack_Clover.push_back(clv);
+        AddBack_Crystal.push_back(cry);
+        AddBack_Segment.push_back(seg);
+        TVector3 position = GetSegmentPosition(clv,cry,seg);
+        AddBack_X.push_back(position.X());
+        AddBack_Y.push_back(position.Y());
+        AddBack_Z.push_back(position.Z());
+        AddBack_Theta.push_back(position.Theta()); 
+        AddBack_E.push_back(energy);      
+        // Define reference axis as the beam direction
+        //beta.Unit().Dump();
+        //position.Dump();
+        position.RotateUz(beta.Unit()); // CHECK
+        //position.Dump();
+        //AddBack_DC.push_back(energy); // Doppler Corrected for highest energy
+        AddBack_DC.push_back(GetDopplerCorrectedEnergy(energy, position, beta)); // Doppler Corrected for highest energy
+        }
+      else{
+        AddBack_E.back()+=energy;      // E1+E2+E3...
+        AddBack_DC.back()+=energy;     // DC(E1)+E2+E3...
+        }
+      }
+    } 
+    else 
+      cout << " Addback scheme " << scheme << " is not supported " << endl;
+   
+
+} // end of add back
+
 /////////////////////////////////////////////////
 void TGeTAMUPhysics::AddClover(unsigned int ID ,double R, double Theta, double Phi){
   TVector3 Pos(0,0,1);
@@ -454,16 +500,6 @@ TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& Se
     cout << "Warning: GeTAMU segment number " << SegmentNbr 
     << " is out of range\n accepted values: 0 (reserved for core) or 1-3 (Left, Middle, Right segment) " << endl;
 
-/* // Not need in case of geTAMU, CHECK
-  // Each crystal is a rotation of the previous one
-  if (CoreNbr == 2 )
-    Pos.RotateZ(90*deg);
-  else if (CoreNbr == 3 )
-    Pos.RotateZ(180*deg);
-  else if (CoreNbr == 4)
-    Pos.RotateZ(270*deg);
-*/
-
   // Define reference axis as the Clover position, 
   // This is a special case to GeTAMU where segmentation is with respect to clover
   Pos.RotateUz(CloverPos.Unit());
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
index 05ea52d91..230f2f913 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
@@ -88,14 +88,14 @@ class TGeTAMUPhysics :  public TObject, public NPL::VDetector{
 
   public: // Data Member
     //singles sorting tools
-    map<int, vector <int> > Singles_CloverMap_CryEN; // cry number energy
-    map<int, vector <int> > Singles_CloverMap_SegEN; // seg number
-    map<int, vector <double> > Singles_CloverMap_CryE; // cry energy
-    map<int, vector <double> > Singles_CloverMap_SegE; // seg energy 
-    map<int, vector <int> > Singles_CloverMap_CryTN; // cry number time
-    map<int, vector <int> > Singles_CloverMap_SegTN; // seg number
-    map<int, vector <double> > Singles_CloverMap_CryT; // cry energy
-    map<int, vector <double> > Singles_CloverMap_SegT; // seg energy 
+    map<int, vector <int> > Singles_CloverMap_CryEN; //! cry number energy
+    map<int, vector <int> > Singles_CloverMap_SegEN; //1 seg number
+    map<int, vector <double> > Singles_CloverMap_CryE; //! cry energy
+    map<int, vector <double> > Singles_CloverMap_SegE; //! seg energy 
+    map<int, vector <int> > Singles_CloverMap_CryTN; //! cry number time
+    map<int, vector <int> > Singles_CloverMap_SegTN; //! seg number
+    map<int, vector <double> > Singles_CloverMap_CryT; //! cry energy
+    map<int, vector <double> > Singles_CloverMap_SegT; //! seg energy 
     //sorting parameters
     vector<double> Singles_E;    
     vector<double> Singles_T;    
@@ -132,6 +132,7 @@ class TGeTAMUPhysics :  public TObject, public NPL::VDetector{
     TVector3 GetCloverPosition(int& CloverNbr);
     TVector3 GetCorePosition(int& CloverNbr, int& CoreNbr);
     TVector3 GetSegmentPosition(int& CloverNbr, int& CoreNbr, int& SegmentNbr);
+    void FillAddBack(int scheme, TVector3& beta);
     inline TVector3 GetCrystalPosition(int& CloverNbr, int& CoreNbr){return GetCorePosition(CloverNbr,CoreNbr);};
 
   private:
-- 
GitLab