diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx
index 191eb57f0851b85d6a9ffaa4dd2584dae87d956e..0f1eaa4c47a05fbf43fe5ff1ee11a8ac11b47b2d 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx
@@ -41,6 +41,10 @@ void TGeTAMUData::Clear(){
   fGeTAMU_Core_CloverNbr_E.clear();
   fGeTAMU_Core_CrystalNbr_E.clear();
   fGeTAMU_Core_Energy.clear();
+  
+  fGeTAMU_Core_CloverNbr_E_LG.clear();
+  fGeTAMU_Core_CrystalNbr_E_LG.clear();
+  fGeTAMU_Core_Energy_LG.clear();
 
   fGeTAMU_Core_CloverNbr_T.clear();
   fGeTAMU_Core_CrystalNbr_T.clear();
@@ -49,6 +53,10 @@ void TGeTAMUData::Clear(){
   fGeTAMU_Segment_CloverNbr_E.clear();
   fGeTAMU_Segment_SegmentNbr_E.clear();
   fGeTAMU_Segment_Energy.clear();
+  
+  fGeTAMU_Segment_CloverNbr_E_LG.clear();
+  fGeTAMU_Segment_SegmentNbr_E_LG.clear();
+  fGeTAMU_Segment_Energy_LG.clear();
 
   fGeTAMU_Segment_CloverNbr_T.clear();
   fGeTAMU_Segment_SegmentNbr_T.clear();
@@ -62,7 +70,14 @@ void TGeTAMUData::Dump() const{
   for (unsigned int i = 0; i <  fGeTAMU_Core_CloverNbr_E.size(); i++){
     cout << " Clover: " <<  fGeTAMU_Core_CloverNbr_E[i] 
          << " Crystal: " << fGeTAMU_Core_CrystalNbr_E[i]  
-         << " Energy: " <<  fGeTAMU_Core_Energy[i] << endl;  
+         << " Energy: " <<  fGeTAMU_Core_Energy[i]<< endl;  
+
+  }
+ cout << "Core Energy Multiplicity (LowGain) = " << fGeTAMU_Core_CloverNbr_E_LG.size() << endl;  
+  for (unsigned int i = 0; i <  fGeTAMU_Core_CloverNbr_E_LG.size(); i++){
+    cout << " Clover LG: " <<  fGeTAMU_Core_CloverNbr_E_LG[i] 
+         << " Crystal LG: " << fGeTAMU_Core_CrystalNbr_E_LG[i]  
+         << " Energy LG: " <<  fGeTAMU_Core_Energy_LG[i] << endl;  
   }
 
  cout << "Core Time Multiplicity = " << fGeTAMU_Core_CloverNbr_T.size() << endl;  
@@ -79,6 +94,13 @@ void TGeTAMUData::Dump() const{
          << " Energy: " <<  fGeTAMU_Segment_Energy[i] << endl;   
   }
 
+  cout << "Segment Energy Multiplicity (LowGain) = " <<  fGeTAMU_Segment_CloverNbr_E_LG.size() << endl;  
+  for (unsigned int i = 0; i < fGeTAMU_Segment_CloverNbr_E_LG.size(); i++){
+    cout << " Clover: " <<  fGeTAMU_Segment_CloverNbr_E_LG[i]  
+         << " Segment: " << fGeTAMU_Segment_SegmentNbr_E_LG[i]  
+         << " Energy: " <<  fGeTAMU_Segment_Energy_LG[i] << endl;   
+  }
+
   cout << "Segment Time Multiplicity = " <<  fGeTAMU_Segment_CloverNbr_T.size() << endl;  
   for (unsigned int i = 0; i < fGeTAMU_Segment_CloverNbr_T.size(); i++){
     cout << " Clover: " <<  fGeTAMU_Segment_CloverNbr_T[i]  
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUData.h b/NPLib/Detectors/GeTAMU/TGeTAMUData.h
index a97a2d422806752d92fbf2321023e5368cd3424e..cc7e657ce361af6d11692c18ae385ecae81b5864 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUData.h
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUData.h
@@ -11,7 +11,7 @@
  * Original Author: Adrien MATTA  contact address: a.matta@surrey.ac.uk      *
  *                                                                           *
  * Creation Date  : 2016                                                     *
- * Last update    : December 2016 m.moukaddam@suurey.ac.uk                   *
+ * Last update    : December 2016 m.moukaddam@surrey.ac.uk                   *
  *---------------------------------------------------------------------------*
  * Decription:                                                               *
  *  This class hold the GeTAMU  raw data                                     *
@@ -37,6 +37,9 @@ private:
   vector<unsigned short> fGeTAMU_Core_CloverNbr_E;
   vector<unsigned short> fGeTAMU_Core_CrystalNbr_E;
   vector<double> fGeTAMU_Core_Energy;
+  vector<unsigned short> fGeTAMU_Core_CloverNbr_E_LG;
+  vector<unsigned short> fGeTAMU_Core_CrystalNbr_E_LG;
+  vector<double> fGeTAMU_Core_Energy_LG;
   vector<unsigned short> fGeTAMU_Core_CloverNbr_T;
   vector<unsigned short> fGeTAMU_Core_CrystalNbr_T;
   vector<double> fGeTAMU_Core_Time;
@@ -44,6 +47,9 @@ private:
   vector<unsigned short> fGeTAMU_Segment_CloverNbr_E;
   vector<unsigned short> fGeTAMU_Segment_SegmentNbr_E;
   vector<double> fGeTAMU_Segment_Energy;
+  vector<unsigned short> fGeTAMU_Segment_CloverNbr_E_LG;
+  vector<unsigned short> fGeTAMU_Segment_SegmentNbr_E_LG;
+  vector<double> fGeTAMU_Segment_Energy_LG;
   vector<unsigned short> fGeTAMU_Segment_CloverNbr_T;
   vector<unsigned short> fGeTAMU_Segment_SegmentNbr_T;
   vector<double> fGeTAMU_Segment_Time;
@@ -62,6 +68,10 @@ public:
   inline void SetCoreCrystalNbrE(const unsigned short &CoreCrystalNbr){fGeTAMU_Core_CrystalNbr_E.push_back(CoreCrystalNbr);}
   inline void SetCoreEnergy(const double &CoreEnergy){fGeTAMU_Core_Energy.push_back(CoreEnergy);}
 
+  inline void SetCoreCloverNbrELowGain(const unsigned short &CoreCloverNbr){fGeTAMU_Core_CloverNbr_E_LG.push_back(CoreCloverNbr); }
+  inline void SetCoreCrystalNbrELowGain(const unsigned short &CoreCrystalNbr){fGeTAMU_Core_CrystalNbr_E_LG.push_back(CoreCrystalNbr);}
+  inline void SetCoreEnergyLowGain(const double &CoreEnergy){fGeTAMU_Core_Energy_LG.push_back(CoreEnergy);}
+
   inline void SetCoreCloverNbrT(const unsigned short &CoreCloverNbr){fGeTAMU_Core_CloverNbr_T.push_back(CoreCloverNbr); }
   inline void SetCoreCrystalNbrT(const unsigned short &CoreCrystalNbr){fGeTAMU_Core_CrystalNbr_T.push_back(CoreCrystalNbr);}
   inline void SetCoreTime(const double &CoreTime){fGeTAMU_Core_Time.push_back(CoreTime);}
@@ -73,6 +83,13 @@ public:
     fGeTAMU_Core_CrystalNbr_E.push_back(CoreCrystalNbr);
     fGeTAMU_Core_Energy.push_back(CoreEnergy);}
 
+  inline void SetCoreELowGain(const unsigned short &CoreCloverNbr,
+                      const unsigned short &CoreCrystalNbr,
+                      const double &CoreEnergy){
+    fGeTAMU_Core_CloverNbr_E_LG.push_back(CoreCloverNbr);
+    fGeTAMU_Core_CrystalNbr_E_LG.push_back(CoreCrystalNbr);
+    fGeTAMU_Core_Energy_LG.push_back(CoreEnergy);}
+
   inline void SetCoreT(const unsigned short &CoreCloverNbr,
                       const unsigned short &CoreCrystalNbr,
                       const double &CoreTime){
@@ -83,6 +100,10 @@ public:
   inline void SetSegmentCloverNbrE(const unsigned short &SegmentCloverNbr){fGeTAMU_Segment_CloverNbr_E.push_back(SegmentCloverNbr); }
   inline void SetSegmentSegmentNbrE(const unsigned short &SegmentNbr){fGeTAMU_Segment_SegmentNbr_E.push_back(SegmentNbr);}
   inline void SetSegmentEnergy(const double &SegmentEnergy){fGeTAMU_Segment_Energy.push_back(SegmentEnergy);}
+
+  inline void SetSegmentCloverNbrELowGain(const unsigned short &SegmentCloverNbr){fGeTAMU_Segment_CloverNbr_E_LG.push_back(SegmentCloverNbr); }
+  inline void SetSegmentSegmentNbrELowGain(const unsigned short &SegmentNbr){fGeTAMU_Segment_SegmentNbr_E_LG.push_back(SegmentNbr);}
+  inline void SetSegmentEnergyLowGain(const double &SegmentEnergy){fGeTAMU_Segment_Energy_LG.push_back(SegmentEnergy);}
   
   inline void SetSegmentCloverNbrT(const unsigned short &SegmentCloverNbr){fGeTAMU_Segment_CloverNbr_T.push_back(SegmentCloverNbr); }
   inline void SetSegmentSegmentNbrT(const unsigned short &SegmentNbr){fGeTAMU_Segment_SegmentNbr_T.push_back(SegmentNbr);}
@@ -95,6 +116,14 @@ public:
     fGeTAMU_Segment_SegmentNbr_E.push_back(SegmentNbr);
     fGeTAMU_Segment_Energy.push_back(SegmentEnergy);};
 
+  
+  inline void SetSegmentELowGain(const unsigned short &SegmentCloverNbr,
+                      const unsigned short &SegmentNbr,
+                      const double &SegmentEnergy){
+    fGeTAMU_Segment_CloverNbr_E_LG.push_back(SegmentCloverNbr);
+    fGeTAMU_Segment_SegmentNbr_E_LG.push_back(SegmentNbr);
+    fGeTAMU_Segment_Energy_LG.push_back(SegmentEnergy);};
+
   inline void SetSegmentT(const unsigned short &SegmentCloverNbr,
                       const unsigned short &SegmentNbr,
                       const double &SegmentTime){
@@ -108,6 +137,11 @@ public:
   inline unsigned short GetCoreCrystalNbrE(const unsigned int &i)  {return fGeTAMU_Core_CrystalNbr_E[i]; }
   inline double GetCoreEnergy(const unsigned int &i)      {return fGeTAMU_Core_Energy[i];}
 
+  inline unsigned int GetMultiplicityCoreELowGain()  {return fGeTAMU_Core_CloverNbr_E_LG.size();}
+  inline unsigned short GetCoreCloverNbrELowGain(const unsigned int &i)   {return fGeTAMU_Core_CloverNbr_E_LG[i]; }
+  inline unsigned short GetCoreCrystalNbrELowGain(const unsigned int &i)  {return fGeTAMU_Core_CrystalNbr_E_LG[i]; }
+  inline double GetCoreEnergyLowGain(const unsigned int &i)      {return fGeTAMU_Core_Energy_LG[i];}
+
   inline unsigned int GetMultiplicityCoreT()  {return fGeTAMU_Core_CloverNbr_T.size();}
   inline unsigned short GetCoreCloverNbrT(const unsigned int &i)   {return fGeTAMU_Core_CloverNbr_T[i]; }
   inline unsigned short GetCoreCrystalNbrT(const unsigned int &i)  {return fGeTAMU_Core_CrystalNbr_T[i]; }
@@ -117,6 +151,11 @@ public:
   inline unsigned short GetSegmentCloverNbrE(const unsigned int &i)   {return fGeTAMU_Segment_CloverNbr_E[i]; }
   inline unsigned short GetSegmentSegmentNbrE(const unsigned int &i)   {return fGeTAMU_Segment_SegmentNbr_E[i]; }
   inline double GetSegmentEnergy(const unsigned int &i)      {return fGeTAMU_Segment_Energy[i];}
+
+   inline unsigned int GetMultiplicitySegmentELowGain()  {return fGeTAMU_Segment_CloverNbr_E_LG.size();}
+  inline unsigned short GetSegmentCloverNbrELowGain(const unsigned int &i)   {return fGeTAMU_Segment_CloverNbr_E_LG[i]; }
+  inline unsigned short GetSegmentSegmentNbrELowGain(const unsigned int &i)   {return fGeTAMU_Segment_SegmentNbr_E_LG[i]; }
+  inline double GetSegmentEnergyLowGain(const unsigned int &i)      {return fGeTAMU_Segment_Energy_LG[i];}
   
   inline unsigned int GetMultiplicitySegmentT()  {return fGeTAMU_Segment_CloverNbr_T.size();}
   inline unsigned short GetSegmentCloverNbrT(const unsigned int &i)   {return fGeTAMU_Segment_CloverNbr_T[i]; }
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
index 5f0260365a90a39111d3bda8de1b36d0c7d47f28..a8c6e55eab1f8cb3d448a73f1ebaa31715864932 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
@@ -54,18 +54,18 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
   PreTreat();
 
  //Treat singles
-  for(unsigned int iSeg = 0 ; iSeg < Singles_SegE.size() ; iSeg++){
-    int clv = Singles_Clover[iSeg];
-    int cry = Singles_Crystal[iSeg];
-    int seg = Singles_Segment[iSeg];
-    double energy = Singles_SegE[iSeg];
-    TVector3 pos = GetPositionOfInteraction(iSeg);
-
-    Singles_Theta.push_back(pos.Theta()/deg);
-    Singles_X.push_back(pos.X());
-    Singles_Y.push_back(pos.Y());
-    Singles_Z.push_back(pos.Z());
-    }
+
+  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];
+    double X = Singles_X[iPixel];
+    double Y = Singles_Y[iPixel];
+    double Z = Singles_Z[iPixel];
+    double Theta = Singles_Theta[iPixel];
+    //cout << clv << " " << cry << " " << seg << " " << X << " " << Y << " " << Z << " "<< energy << endl;
+     }
 
   // Treat addback
   unsigned int c_size_e = m_PreTreatedData->GetMultiplicityCoreE();
@@ -77,7 +77,7 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
   map<int,double> clv_energy;
   map<int,int> clv_segment;
   map<int,int> clv_crystal;
-  map<int,double> max_core;
+  map<int,double> max_cry;
   map<int,double> max_segment;
 
   map<int,TVector3> max_pos; // first hit
@@ -93,8 +93,8 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
     // Add back energy
     clv_energy[clv] += energy;
     // Pick up the crystal with the maximum energy in every clover
-    if(energy > max_core[clv]){
-      max_core[clv] = energy;
+    if(energy > max_cry[clv]){
+      max_cry[clv] = energy;
       clv_crystal[clv] = cry;
     }
     // Pick up the segment with the maximum energy in every clover
@@ -132,8 +132,8 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
     // Add back energy
     clv_energy[clv] += energy;
     // Pick up the crystal with the maximum energy in every clover
-    if(energy > max_core[clv]){
-      max_core[clv] = energy;
+    if(energy > max_cry[clv]){
+      max_cry[clv] = energy;
       clv_crystal[clv] = cry;
     }
     // Pick up the segment with the maximum energy in every clover
@@ -175,7 +175,7 @@ void TGeTAMUPhysics::BuildPhysicalEvent(){
 for (unsigned i = 0 ; i < m_PreTreatedData->GetMultiplicityCoreT(); i++)
   GeTime.push_back(m_PreTreatedData->GetCoreTime(i));
 
-}
+} // end
 
 /////////////////////////////////////////////////
 void TGeTAMUPhysics::PreTreat(){
@@ -197,6 +197,8 @@ void TGeTAMUPhysics::PreTreat(){
       crystal = m_EventData->GetCoreCrystalNbrE(i);
       name = "GETAMU/D"+ NPL::itoa(clover)+"_CRY"+ NPL::itoa(crystal);
       Energy =  cal->ApplyCalibration(name+"_E", Eraw);
+      Singles_CloverMap_CryEN[clover].push_back(crystal);
+      Singles_CloverMap_CryE[clover].push_back(Energy);
       m_PreTreatedData->SetCoreE(clover,crystal,Energy);
     }
   }
@@ -208,6 +210,8 @@ void TGeTAMUPhysics::PreTreat(){
       crystal = m_EventData->GetCoreCrystalNbrT(i);
       name = "GETAMU/D"+ NPL::itoa(clover)+"_CRY"+ NPL::itoa(crystal);
       Time =  cal->ApplyCalibration(name+"_T", Traw);
+      Singles_CloverMap_CryTN[clover].push_back(crystal);
+      Singles_CloverMap_CryT[clover].push_back(Time);
       m_PreTreatedData->SetCoreT(clover,crystal,Time);
     }
   }
@@ -220,6 +224,8 @@ void TGeTAMUPhysics::PreTreat(){
       segment = m_EventData->GetSegmentSegmentNbrE(i);
       name = "GETAMU/D"+ NPL::itoa(clover)+"_SEG"+ NPL::itoa(segment);
       Energy =  cal->ApplyCalibration(name+"_E", Eraw);
+      Singles_CloverMap_SegEN[clover].push_back(segment);
+      Singles_CloverMap_SegE[clover].push_back(Energy);
       m_PreTreatedData->SetSegmentE(clover,segment,Energy);
     }
   }
@@ -232,34 +238,136 @@ void TGeTAMUPhysics::PreTreat(){
       segment = m_EventData->GetSegmentSegmentNbrT(i);
       name = "GETAMU/D"+ NPL::itoa(clover)+"_SEG"+ NPL::itoa(segment);
       Time =  cal->ApplyCalibration(name+"_T", Traw);
+      Singles_CloverMap_CryTN[clover].push_back(segment);
+      Singles_CloverMap_CryT[clover].push_back(Time);
       m_PreTreatedData->SetSegmentT(clover,segment,Time);
     }
   }
 
-//match data and fill singles
-  unsigned int sizec = m_PreTreatedData->GetMultiplicityCoreE();
-  unsigned int sizes = m_PreTreatedData->GetMultiplicitySegmentE();
-  //cores
-  for(unsigned int iCore = 0 ; iCore < sizec ; iCore++){
-    Energy = m_PreTreatedData->GetCoreEnergy(iCore);
-    Singles_CoreE.push_back(Energy);
+if(m_PreTreatedData->GetMultiplicityCoreE())   FillSingles();
+
+}
+
+void TGeTAMUPhysics::FillSingles(void){
+//match segments data to cores data,
+// Rule 0: If more than 1 core is fired, there are as many hits as there are fired cores
+// Rule 1: segments are auxilliary to cores, they only point out the reduced "pixel" if any, (pixel== overlap segment/core)
+//Ideas: try to kick one segment out by sum rule
+
+//cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
+//m_PreTreatedData->Dump();
+
+vector <int>    CryEN, SegEN;
+vector <double> CryE, SegE;
+
+for (unsigned iClover=0; iClover<4; iClover++){
+  int clover = iClover+1;
+  CryEN.clear();
+  SegEN.clear();
+  CryE.clear();
+  SegE.clear();
+  //Energy
+  if(Singles_CloverMap_CryEN.find(iClover+1) != Singles_CloverMap_CryEN.end()){
+    CryEN  = Singles_CloverMap_CryEN[iClover+1];
+    CryE   = Singles_CloverMap_CryE[iClover+1];
   }
-  //segments
-  // For geometry parameters (clover, crystal, segment) we use the segment vectors
-  for(unsigned int iSeg = 0 ; iSeg < sizes ; iSeg++){
-    clover = m_PreTreatedData->GetSegmentCloverNbrE(iSeg);
-    crystal = m_PreTreatedData->GetSegmentCloverNbrE(iSeg);
-    segment = m_PreTreatedData->GetSegmentSegmentNbrE(iSeg);
-    Energy = m_PreTreatedData->GetSegmentEnergy(iSeg);
-    //avoid impossible cases
-    if ( !(segment==0 && (crystal==2||crystal==3)) || (segment==2 && (crystal==1||crystal==4)) ){
-      Singles_Clover.push_back(clover);
-      Singles_Crystal.push_back(crystal);
-      Singles_Segment.push_back(segment);
-      Singles_SegE.push_back(Energy);
+  else
+    continue; // no need to go further if Cores energies are non existant
+
+  if(Singles_CloverMap_SegEN.find(iClover+1) != Singles_CloverMap_SegEN.end()){
+    SegEN   = Singles_CloverMap_SegEN[iClover+1];
+    SegE    = Singles_CloverMap_SegE[iClover+1];
+  }
+
+/*
+//calculate total energy for sum rules
+  double totCryE =0;
+  double totSegE =0;
+  for (unsigned i=0 ; i < CryE.size();i++)    totCryE+=CryE[i];
+  for (unsigned i=0 ; i < SegE.size();i++)   totSegE+=SegE[i];
+*/
+
+//sort the crystal energies;
+  int swapEN;
+  double swapE;
+  int size = (int) CryE.size();
+for (int i=0; i< (size -1); i++){    // element to be compared
+  for(int j = (i+1); j < size; j++){   // rest of the elements
+    if (CryE[i] < CryE[j]){          // descending order
+      swapE= CryE[i];    CryE[i] = CryE[j];   CryE[j] = swapE;
+      swapEN= CryEN[i];   CryEN[i] = CryEN[j];   CryEN[j] = swapEN;
     }
   }
+}
 
+  //create hit matrix
+  int pixel[4][3];
+  memset(pixel,0,sizeof(pixel)); // filling the matrix with zeros
+  //fill the cores
+  for (unsigned iCry = 0 ; iCry < CryEN.size() ; iCry++){
+    int Cry = CryEN[iCry] ;
+    for (unsigned iSeg = 1 ; iSeg <= 3 ; iSeg++){
+      // (Segment 3 (right) + Cry 1 or 4) or (Segment 1 (Left)  + Cry 2 or 3) are impossible
+      if( ((Cry==1 || Cry==4) && iSeg!=3) ||  ((Cry==2 || Cry==3) && iSeg!=1) ){
+        pixel[Cry-1][iSeg-1]++;
+      }
+    }
+  }
+  //fill the segments
+  for (unsigned iSeg = 0 ; iSeg < SegEN.size() ; iSeg++){
+    int Seg = SegEN[iSeg];
+    for (unsigned iCry = 1 ; iCry <= 4 ; iCry++){
+      if( pixel[iCry-1][Seg-1] != 0){ // Only if the pixel was filled by a Cry signal previously, access it
+        pixel[iCry-1][Seg-1]++;
+      }
+    }
+  }
+
+  // show
+  /*
+  int* apixel = *pixel;
+  cout <<endl<<"------- Clover "<< clover << endl;
+  for (unsigned i = 0 ; i < 4 ; i++){
+    for (unsigned j = 0 ; j < 3 ; j++)
+        cout << *(apixel+(i*3)+j) << " ";
+    cout << endl;
+  }
+   cout <<"----------------------- " <<endl;
+  */
+   //Calculate the singles
+  for (unsigned i = 0 ; i < CryEN.size() ; i++){
+    int segment = -1;
+    unsigned crystal = CryEN[i];
+    unsigned segmentA = 2;
+    unsigned segmentB = 3;
+    if (crystal==1 || crystal==4){ // if Core 1 or 4 change the segments to segment 1 and 2
+      segmentA = 1;
+      segmentB = 2;
+    }
+    //pick between segment A or B for each case
+    if (pixel[crystal-1][segmentA-1] == pixel[crystal-1][segmentB-1])
+      segment = 0; // system can't be resolved
+    else if (pixel[crystal-1][segmentA-1] > pixel[crystal-1][segmentB-1])
+      segment = segmentA;
+    else if (pixel[crystal-1][segmentA-1] < pixel[crystal-1][segmentB-1])
+       segment = segmentB;
+
+    //cout << i <<" picked: crystal " << crystal << "   segment " << segment << "  Energy " << CryE[i] << endl;
+    Singles_Clover.push_back(clover);
+    Singles_Crystal.push_back(CryEN[i]);
+    Singles_Segment.push_back(segment);
+    Singles_E.push_back(CryE[i]);
+    TVector3 Pos = GetSegmentPosition(clover,CryEN[i],segment);
+    Singles_X.push_back(Pos.X());
+    Singles_Y.push_back(Pos.Y());
+    Singles_Z.push_back(Pos.Z());
+    Singles_Theta.push_back(Pos.Theta());
+    //cout << " XYZ "<< Pos.X() << " "<< Pos.Y() << " "<< Pos.Z() << " Theta: " <<Pos.Theta()/deg<< endl ;
+    }
+  //cin.get();
+//Time
+
+  } // end of Clover loop on map
 
 }
 
@@ -296,21 +404,21 @@ TVector3 TGeTAMUPhysics::GetCloverPosition(int& CloverNbr){
 }
 /////////////////////////////////////////////////
 TVector3 TGeTAMUPhysics::GetCorePosition(int& CloverNbr,int& CoreNbr){
-  static double offset = 33.4; // mm
-  static double depth = 45; // mm
-  static TVector3 Pos;
-  TVector3 CloverPos = m_CloverPosition[CloverNbr];
+  double offset = 20; // mm
+  double depth = 40; // mm
+  TVector3 Pos;
+  TVector3 CloverPos = GetCloverPosition(CloverNbr);
 
   if(CoreNbr==1)
-    Pos.SetXYZ(-offset,offset,depth);
-  else if(CoreNbr==2)
     Pos.SetXYZ(offset,offset,depth);
+  else if(CoreNbr==2)
+    Pos.SetXYZ(-offset,offset,depth);
   else if(CoreNbr==3)
-    Pos.SetXYZ(offset,-offset,depth);
-  else if(CoreNbr==4)
     Pos.SetXYZ(-offset,-offset,depth);
+  else if(CoreNbr==4)
+    Pos.SetXYZ(offset,-offset,depth);
   else
-    cout << "Warning: GeTAMU crystal number " << CoreNbr << " is out of range (1 to 4)" << endl;
+    cout << "Warning in GetCorePosition: GeTAMU crystal number " << CoreNbr << " is out of range (1 to 4)" << endl;
 
   // Define reference axis as the clover direction
   Pos.RotateUz(CloverPos.Unit());
@@ -320,28 +428,22 @@ TVector3 TGeTAMUPhysics::GetCorePosition(int& CloverNbr,int& CoreNbr){
 /////////////////////////////////////////////////
 TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& SegmentNbr){
 
-  static double offsetX = 33.4; // mm assumed the same as TIGRESS, CHECK
-  static double offsetY = 33.4; //mm in case of left and right segments, CHECK
-  static double depth = 45; // mm in crystal, assumed the same as TIGRESS, CHECK
+  if (SegmentNbr==0) return GetCorePosition(CloverNbr,CoreNbr);
+
+   double offsetX = 20; // 20mm in GeTAMU according to measurments, 33.4 mm in TIGRESS
+   double offsetY = 20;
+   double depth = 40;  // 40mm in GeTAMU according to measurments, 45 mm in TIGRESS
 
   // Changes signs with segment/core combinations
-  if (CoreNbr==1||CoreNbr==4)
+  if (CoreNbr==3||CoreNbr==2)
     offsetX = -offsetX;
   if (CoreNbr==3||CoreNbr==4)
     offsetY = -offsetY;
 
-  //impossible cases
-  if ( (SegmentNbr==0 && (CoreNbr==2||CoreNbr==3)) || (SegmentNbr==2 && (CoreNbr==1||CoreNbr==4)) )
-    cout << "Warning: GeTAMU segment number " << SegmentNbr
-    << " and core number " << CoreNbr << " are not compatible " << endl;
-
-  TVector3 CorePos = GetCorePosition(CloverNbr,CoreNbr);
   TVector3 CloverPos = GetCloverPosition(CloverNbr);
-  static TVector3 Pos;
+  TVector3 Pos;
 
-  if(SegmentNbr == 0)
-    return CorePos;
-  else if(SegmentNbr==1 || SegmentNbr==3){ // LEFT OR RIGHT
+  if(SegmentNbr==1 || SegmentNbr==3){ // LEFT OR RIGHT
     offsetX = 1.5*offsetX ;
     Pos.SetXYZ(offsetX,offsetY,depth);
     }
@@ -350,8 +452,9 @@ TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& Se
     Pos.SetXYZ(offsetX,offsetY,depth);
     }
   else
+
     cout << "Warning: GeTAMU segment number " << SegmentNbr
-    << " is out of range\n accepted range: 0 (core) or 1-3 (L,M,R) " << endl;
+    << " 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
@@ -415,8 +518,15 @@ void TGeTAMUPhysics::InitializeRootOutput()    {
 ///////////////////////////////////////////////////////////////////////////
 void TGeTAMUPhysics::Clear() {
 
-  Singles_CoreE.clear();
-  Singles_SegE.clear();
+  Singles_CloverMap_CryEN.clear(); // cry number energy
+  Singles_CloverMap_SegEN.clear(); // seg number
+  Singles_CloverMap_CryE.clear(); // cry energy
+  Singles_CloverMap_SegE.clear(); // seg energy
+  Singles_CloverMap_CryTN.clear(); // cry number time
+  Singles_CloverMap_SegTN.clear(); // seg number
+  Singles_CloverMap_CryT.clear(); // cry energy
+  Singles_CloverMap_SegT.clear(); // seg energy
+  Singles_E.clear();
   Singles_DC.clear();
   Singles_Theta.clear();
   Singles_X.clear();
@@ -443,14 +553,19 @@ void TGeTAMUPhysics::Clear() {
 ///////////////////////////////////////////////////////////////////////////
 void TGeTAMUPhysics::AddParameterToCalibrationManager(){
   CalibrationManager* Cal = CalibrationManager::getInstance();
-  for(int i = 0 ; i < 16; ++i){
-    for(int cry = 0 ; cry < 4 ; cry++){
-      // core are 0 and 9 , segment 1 to 8
-      for( int j = 0 ; j < 10 ; ++j){
-        Cal->AddParameter("GETAMU", "D"+ NPL::itoa(i+1)+"_CRY"+NPL::itoa(cry+1)+"_SEG"+ NPL::itoa(j)+"_E","GETAMU_D"+ NPL::itoa(i+1)+"_CRY"+NPL::itoa(cry+1)+"_SEG"+NPL::itoa(j)+"_E");
-      }
+  for(int det = 0 ; det < 4; det++){ //4 detectors
+
+    for(int cry = 0 ; cry < 4 ; cry++){ // 4 crystals
+      Cal->AddParameter("GETAMU", "D"+ NPL::itoa(det+1)+"_CRY"+NPL::itoa(cry+1)+"_E","GETAMU_D"+ NPL::itoa(det+1)+"_CRY"+NPL::itoa(cry+1)+"_E");
+      Cal->AddParameter("GETAMU", "D"+ NPL::itoa(det+1)+"_CRY"+NPL::itoa(cry+1)+"_E","GETAMU_D"+ NPL::itoa(det+1)+"_CRY"+NPL::itoa(cry+1)+"_T");
+    }
+    for( int seg = 0 ; seg < 3 ; seg++){ // 3 segments
+      Cal->AddParameter("GETAMU", "D"+ NPL::itoa(det+1)+"_SEG"+ NPL::itoa(seg+1)+"_E","GETAMU_D"+ NPL::itoa(det+1)+"_SEG"+NPL::itoa(seg+1)+"_E");
+      Cal->AddParameter("GETAMU", "D"+ NPL::itoa(det+1)+"_SEG"+ NPL::itoa(seg+1)+"_T","GETAMU_D"+ NPL::itoa(det+1)+"_SEG"+NPL::itoa(seg+1)+"_T");
     }
+
   }
+
   return;
 }
 
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
index 31e86b0494118a47978cccb6af6edbec060c8193..05ea52d9197150469d1db7f157cfa91ec357ed8a 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
@@ -87,10 +87,19 @@ class TGeTAMUPhysics :  public TObject, public NPL::VDetector{
     TGeTAMUPhysics* m_EventPhysics;//!
 
   public: // Data Member
-    //singles
-    vector<double> Singles_CoreE; 
-    vector<double> Singles_SegE;   
-    vector<double> Singles_DC;   
+    //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 
+    //sorting parameters
+    vector<double> Singles_E;    
+    vector<double> Singles_T;    
+    vector<double> Singles_DC;   // Doppler Corrected Energy
     vector<double> Singles_Theta;
     vector<double> Singles_X;
     vector<double> Singles_Y;
@@ -98,11 +107,10 @@ class TGeTAMUPhysics :  public TObject, public NPL::VDetector{
     vector<int> Singles_Clover;
     vector<int> Singles_Crystal;
     vector<int> Singles_Segment;
-
-    // add back by clover
+    // add back parameters
     vector<double> AddBack_E;   
     vector<double> AddBack_T;   
-    vector<double> AddBack_DC;   
+    vector<double> AddBack_DC;    // Doppler Corrected Energy
     vector<double> AddBack_Theta;
     vector<double> AddBack_X;
     vector<double> AddBack_Y;
@@ -128,7 +136,7 @@ class TGeTAMUPhysics :  public TObject, public NPL::VDetector{
 
   private:
     map<unsigned int,TVector3> m_CloverPosition;//!
-
+    void FillSingles(void);
   public: // Static constructor to be passed to the Detector Factory
     static NPL::VDetector* Construct();
     ClassDef(TGeTAMUPhysics,1)  // GeTAMUPhysics structure
diff --git a/Projects/T40/Analysis.cxx b/Projects/T40/Analysis.cxx
index fcca74d16991b730644fc52624930d12ce7a8f93..372228a4b26f72ad1bf946f17f8e90a92113d464 100644
--- a/Projects/T40/Analysis.cxx
+++ b/Projects/T40/Analysis.cxx
@@ -97,11 +97,11 @@ void Analysis::Init(){
   // energy losses
   string light=NPL::ChangeNameToG4Standard(myReaction->GetNucleus3().GetName());
   string beam=NPL::ChangeNameToG4Standard(myReaction->GetNucleus1().GetName());
-  LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".SRIM","SRIM",10 );
-  LightAl = NPL::EnergyLoss(light+"_Al.SRIM","SRIM",10);
-  //LightSi = NPL::EnergyLoss(light+"_Si.SRIM","SRIM",10);
+  LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",10 );
+  LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",10);
+  //LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",10);
   LightSi = NPL::EnergyLoss("He4_Si.SRIM","SRIM",10);
-  BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".SRIM","SRIM",10);
+  BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",10);
   FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness*0.5, 0);
   myReaction->SetBeamEnergy(FinalBeamEnergy);
   cout << "Final Beam energy (middle of target): " << FinalBeamEnergy << endl;
@@ -126,6 +126,7 @@ void Analysis::Init(){
   InitInputBranch();
   
   //Ge
+  GammaSinglesE=0;
   
   //FPD
   Delta_E = 0; // Energy ionisation chamber
@@ -157,6 +158,7 @@ void Analysis::Init(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::TreatEvent(){
+
   // Reinitiate calculated variable
   ReInitValue();
   ////////////////////////////////////////// LOOP on TiaraHyball + SSSD Hit //////////////////////////////////////////
@@ -266,49 +268,20 @@ void Analysis::TreatEvent(){
   } // end loop TiaraBarrel
 
   /////////////////////////// LOOP on Ge TAMU /////////////////////////////
-  /*
-for(unsigned int countGe = 0 ; countGe < TG->something.size() ; countGe++) // multiplicity treated for now is zero 
+  
+for(unsigned int countGe = 0 ; countGe < TG->Singles_E.size() ; countGe++) // multiplicity treated for now is zero 
   { 
-
-  unsigned int c_size_e = TG->m_PreTreatedData->GetMultiplicityCoreE();
-  unsigned int s_size_e = m_PreTreatedData->GetMultiplicitySegmentE();
-  // map for add back
-  map<int,double> clv_energy;   
-  map<int,int> clv_segment;
-  map<int,int> clv_crystal;
-  map<int,double> max_core;
-  map<int,double> max_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_core[clv]){
-      max_core[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);
-      }
-    }
-  }
-
   // Singles spectra 
   
   // Addback spectra 
+  GammaSinglesE+= TG->Singles_E[countGe];
 
   // calculate angle
 
   // calculate doppler corrected spectra 
   	
   }
-*/
+
 
  ////////////////////////////////////////// LOOP on FPD  //////////////////////////////////////////
 	// Micromega energy
@@ -472,6 +445,8 @@ void Analysis::InitOutputBranch() {
 
   //GeTamu
   // stuff goes here 
+  RootOutput::getInstance()->GetTree()->Branch("GammaSinglesE",&GammaSinglesE,"GammaSinglesE/D");
+
 
   //FPD
   RootOutput::getInstance()->GetTree()->Branch("Delta_E",&Delta_E,"Delta_E/D");
diff --git a/Projects/T40/Analysis.h b/Projects/T40/Analysis.h
index 47de3961488fddc64e35949615011161cb7e2457..ff0ca6689f2d7d2023717143510c76913934bffc 100644
--- a/Projects/T40/Analysis.h
+++ b/Projects/T40/Analysis.h
@@ -94,6 +94,9 @@ class Analysis: public NPL::VAnalysis{
 	double YTarget ;
 	TVector3 BeamDirection ;
   
+  //Gamma
+  double GammaSinglesE;
+  
   //FPD
 	static const Int_t kNumAw = 4; // number of wires
   double Delta_E      ;