diff --git a/protoduneana/verticaldrift/singlehit/SingleHit_module.cc b/protoduneana/verticaldrift/singlehit/SingleHit_module.cc
index 5ad15d5..ec858ab 100644
--- a/protoduneana/verticaldrift/singlehit/SingleHit_module.cc
+++ b/protoduneana/verticaldrift/singlehit/SingleHit_module.cc
@@ -53,13 +53,14 @@
 #include "TPaveStats.h"
 #include "TGraph.h"
 #include "TEllipse.h"
+#include "TGraph2D.h"
 
 using std::vector;
 using std::string;
 
 typedef struct 
 { 
-  float y, z; 
+  float y, z, t; 
   int group, index;  
 } point_t, *point;
 
@@ -70,7 +71,7 @@ typedef struct
   std::list<int> lChannelCol , lChannelInd1 , lChannelInd2;
   std::vector<int> vMCPDG , vMCMOMpdg , vNOF;
   std::vector<float> vMCWEI;
-  std::vector<float> vMCX , vMCY , vMCZ;
+  std::vector<float> vMCX , vMCY , vMCZ, vMCE, vMCElec;
   std::vector<std::string> vMCGenTag;
 } TempCluster ;                            
                               
@@ -80,7 +81,7 @@ typedef struct
   int Npoint, NCol , NInd1 , NInd2 , NOF;
   std::vector<int> vMCPDG , vMCMOMpdg;
   std::vector<float> vMCWEI;
-  std::vector<float> vMCX , vMCY , vMCZ;
+  std::vector<float> vMCX , vMCY , vMCZ, vMCE, vMCElec;
   std::vector<std::string> vMCGenTag;
 } Cluster ; 
 
@@ -213,6 +214,9 @@ class pdvdana::SingleHit : public art::EDAnalyzer {
   std::vector<int>   vNInd1Cluster; 
   std::vector<int>   vNInd2Cluster;
 
+  float sumMCEnergy = 0 ;
+  float sumMCElec   = 0;
+
   std::vector<std::string> vMCGenTagCluster;
   std::vector<int>         vMCMOMpdgCluster;
   std::vector<int>         vMCPDGCluster;
@@ -220,6 +224,8 @@ class pdvdana::SingleHit : public art::EDAnalyzer {
   std::vector<float>       vMCXCluster;
   std::vector<float>       vMCYCluster;
   std::vector<float>       vMCZCluster;
+  std::vector<float>       vMCECluster;
+  std::vector<float>       vMCElecCluster;
 
   //Input variables
   std::string fSpacePointLabel;
@@ -298,6 +304,8 @@ class pdvdana::SingleHit : public art::EDAnalyzer {
   std::vector<float>       vMCXByEvent;
   std::vector<float>       vMCYByEvent;
   std::vector<float>       vMCZByEvent;
+  std::vector<float>       vMCEnergyByEvent;
+  std::vector<float>       vMCElecByEvent;
 
   std::vector<std::string> vGeneratorTagByEvent;
 
@@ -361,6 +369,7 @@ class pdvdana::SingleHit : public art::EDAnalyzer {
   // CLUSTER FUNCTIONS
   point gen_yz(int size , std::vector<int> vIndex , std::vector<float> vY , std::vector<float> vZ , std::vector<int> vNOF);
 
+  point gen_yzt(int size , std::vector<int> vIndex , std::vector<float> vY , std::vector<float> vZ , std::vector<float> vT , std::vector<int> vNOF , float fElectronVelocity , float fTickToMus);
   float dist2(point a, point b);
 
   void Plot(TCanvas * canvas , int event , int nbin , float ymin , float ymax , float zmin , float zmax , std::vector<std::vector<float> > &data,std::vector<std::vector<float> > &cluster,double RMS );
@@ -392,7 +401,8 @@ class pdvdana::SingleHit : public art::EDAnalyzer {
                                    std::vector<int> vMCPDGByEvent , std::vector<int> vMCMOMpdgByEvent ,std::vector<float> vMCWeightByEvent , 
                                    std::vector<std::string> vGeneratorTagByEvent ,
                                    std::vector<float> vMCXByEvent , std::vector<float> vMCYByEvent , std::vector<float> vMCZByEvent , 
-                                   std::vector<int> vNoFByEvent);
+                                   std::vector<int> vNoFByEvent,
+                                   std::vector<float> vMCEnergyByEvent, std::vector<float> vMCElecByEvent);
 };
 
 
@@ -620,6 +630,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
     vMCXCluster.push_back(-9999);
     vMCYCluster.push_back(-9999);
     vMCZCluster.push_back(-9999);
+    vMCECluster.push_back( -999 );
+    vMCElecCluster.push_back( -999 );
 
     vVetoTrackStartX.push_back(-9999.0);
     vVetoTrackStartY.push_back(-9999.0);
@@ -663,6 +675,9 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   if ( !vMCYByEvent.empty() )          vMCYByEvent.clear();
   if ( !vMCZByEvent.empty() )          vMCZByEvent.clear();
   if ( !vGeneratorTagByEvent.empty() ) vGeneratorTagByEvent.clear();
+  if ( !vMCEnergyByEvent.empty() )     vMCEnergyByEvent.clear();
+  if ( !vMCElecByEvent.empty() )       vMCElecByEvent.clear();
+
 
   for(int index =0 ; index<fNHits; index++)
   {
@@ -683,7 +698,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
       vMCPart_Starty.clear();
       vMCPart_Startz.clear();
 
-      float sumMCEnergy = 0 ;
+      sumMCEnergy = 0 ;
+      sumMCElec   = 0;
 
       using weightedMCPair = std::pair<const simb::MCParticle*, float>;
       std::vector<weightedMCPair> tempMCPair;
@@ -693,6 +709,7 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
         const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
         tempMCPair.push_back(std::make_pair(curr_part,ide.energy));
         sumMCEnergy += ide.energy;
+        sumMCElec   += ide.numElectrons;
       }
 
       
@@ -758,7 +775,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
 	 
     fNearOrFarToTheBeam = NearOrFar(bIsPDVD , bIsPDHD , hit);
 
-    fEnergy         = hit.ROISummedADC();///fADCtoEl;
+    //fEnergy         = hit.ROISummedADC();///fADCtoEl;
+    fEnergy         = hit.Integral();
     fPeakTime       = hit.PeakTime();//*ftick_in_mus;
     fSigmaPeakTime  = hit.SigmaPeakTime();//*ftick_in_mus;
     fRMS            = hit.RMS();
@@ -938,6 +956,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
       vMCZByEvent.push_back( vMCPart_Endz[0] );
 
       vGeneratorTagByEvent.push_back( vGenerator_tag[0] );
+      vMCEnergyByEvent.push_back( sumMCEnergy ) ;
+      vMCElecByEvent.push_back( sumMCElec );
 
       z++;
       ch1++;
@@ -1037,7 +1057,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
     vMCXCluster.push_back(-9999);
     vMCYCluster.push_back(-9999);
     vMCZCluster.push_back(-9999);
-
+    vMCECluster.push_back( -999 );
+    vMCElecCluster.push_back( -999 );
     vVetoTrackStartX.push_back(-9999.0);
     vVetoTrackStartZ.push_back(-9999.0);
     vVetoTrackStartY.push_back(-9999.0);
@@ -1056,12 +1077,28 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   }
 
   point v;
-  v = gen_yz( PTSIsolated , vIso , vYPointByEvent , vZPointByEvent , vNoFByEvent );
-
+  //v = gen_yz( PTSIsolated , vIso , vYPointByEvent , vZPointByEvent , vNoFByEvent );
+  v = gen_yzt( PTSIsolated , vIso , vYPointByEvent , vZPointByEvent , vPeakTimeColByEvent, vNoFByEvent , fElectronVelocity , fTickTimeInMus);
 
   std::vector<std::vector<float> > dataPos = GetData(PTSIsolated,v);
   std::vector<std::vector<float> > clustersPos;
 
+  std::vector<float> dataPosZ = dataPos[0];
+  std::vector<float> dataPosY = dataPos[1];
+  std::vector<float> dataPosT = dataPos[2];
+
+  if ( LogLevel > 20)
+  {
+    art::ServiceHandle<art::TFileService> tfs;
+    TString canvasName = Form("isopoint_%i", fEventID);
+    TCanvas* canvas1 = tfs->make<TCanvas>(canvasName, canvasName);
+    int N = dataPosZ.size();
+    TGraph2D *graph = new TGraph2D(N, &dataPosZ[0] , &dataPosY[0] , &dataPosT[0]);
+    canvas1->cd();
+    graph->Draw("AP");
+    canvas1->Write(canvasName);
+  }
+
   std::vector<int> vchecks(2,0);
 
   int K = fNumberInitClusters;
@@ -1093,7 +1130,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
       }
       else
       {
-        K = vchecks[1] + 1;
+        threshold = fMaxSizeCluster;
+	K = vchecks[1] + 5;
         if ( LogLevel > 2) printf("Threshold Max reached \n");
       }
     }
@@ -1114,7 +1152,7 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   if (LogLevel > 3) std::cout<< " reallocation done " << std::endl;
 
 
-  if ( LogLevel > 2)
+  if ( LogLevel > 20)
   {
     art::ServiceHandle<art::TFileService> tfs;
     TString canvasName = Form("clusterisation_%i", fEventID);
@@ -1123,7 +1161,7 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
     canvas->Write(canvasName);
   } 
 
-  std::vector<Cluster> vCluster = GetCluster( PTSIsolated , clustersPos[0].size() , v , vEInd1PointByEvent , vEInd2PointByEvent , vChInd1PointByEvent , vChInd2PointByEvent , vEnergyColByEvent , vPeakTimeColByEvent , vChannelColByEvent , vMCPDGByEvent , vMCMOMpdgByEvent , vMCWeightByEvent , vGeneratorTagByEvent , vMCXByEvent , vMCYByEvent , vMCZByEvent , vNoFByEvent);
+  std::vector<Cluster> vCluster = GetCluster( PTSIsolated , clustersPos[0].size() , v , vEInd1PointByEvent , vEInd2PointByEvent , vChInd1PointByEvent , vChInd2PointByEvent , vEnergyColByEvent , vPeakTimeColByEvent , vChannelColByEvent , vMCPDGByEvent , vMCMOMpdgByEvent , vMCWeightByEvent , vGeneratorTagByEvent , vMCXByEvent , vMCYByEvent , vMCZByEvent , vNoFByEvent, vMCEnergyByEvent, vMCElecByEvent);
   NCluster = vCluster.size();
 
   if ( !vYCluster.empty() )        vYCluster.clear();
@@ -1144,7 +1182,9 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   if ( !vMCYCluster.empty() )      vMCYCluster.clear();
   if ( !vMCZCluster.empty() )      vMCZCluster.clear();
   if ( !vMCGenTagCluster.empty() ) vMCGenTagCluster.clear();
-  
+  if ( !vMCECluster.empty() )      vMCECluster.clear();
+  if ( !vMCElecCluster.empty() )   vMCElecCluster.clear();  
+
   int NEmpty_cluster = 0;
   for( int j = 0 ; j < NCluster ; j++ )
   {
@@ -1177,6 +1217,9 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
     vMCYCluster.insert( vMCYCluster.end() , (vCluster[j].vMCY).begin() , (vCluster[j].vMCY).end() );
     vMCZCluster.insert( vMCZCluster.end() , (vCluster[j].vMCZ).begin() , (vCluster[j].vMCZ).end() );
 
+    vMCECluster.insert( vMCECluster.end() , (vCluster[j].vMCE).begin() , (vCluster[j].vMCE).end() );
+    vMCElecCluster.insert( vMCElecCluster.end() , (vCluster[j].vMCElec).begin() , (vCluster[j].vMCElec).end() );
+
     vMCGenTagCluster.insert( vMCGenTagCluster.end()    , (vCluster[j].vMCGenTag).begin() , (vCluster[j].vMCGenTag).end() );
   }
 
@@ -1288,6 +1331,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   vMCXCluster.clear();
   vMCYCluster.clear();
   vMCZCluster.clear();
+  vMCECluster.clear();
+  vMCElecCluster.clear();
 
   vYPointByEvent.clear();
   vZPointByEvent.clear();
@@ -1311,6 +1356,8 @@ void pdvdana::SingleHit::analyze(art::Event const& e)
   vMCXByEvent.clear();
   vMCYByEvent.clear();
   vMCZByEvent.clear();
+  vMCEnergyByEvent.clear();
+  vMCElecByEvent.clear();
 
   vVetoTrackStartX.clear();
   vVetoTrackStartY.clear();
@@ -1417,6 +1464,9 @@ void pdvdana::SingleHit::beginJob()
   tClusterTree->Branch("MCParticleY"        , &vMCYCluster      );
   tClusterTree->Branch("MCParticleZ"        , &vMCZCluster      );
   tClusterTree->Branch("HitGenerationTag"   , &vMCGenTagCluster );
+  tClusterTree->Branch("MCParticleEnergy"   , &vMCECluster      );
+  tClusterTree->Branch("MCParticleElec"     , &vMCElecCluster   );
+
   if (LogLevel > 2)
   {
     tClusterTree->Branch("Y_isolated_point" , &vY_point          );
@@ -1621,7 +1671,8 @@ void pdvdana::SingleHit::GetListOfTimeCoincidenceHit( bool IsPDVD , bool IsPDHD
 
       WireInd1.push_back(hit.WireID());
       ChannelInd1.push_back(hit.Channel());
-      EInd1.push_back(hit.ROISummedADC());
+      //EInd1.push_back(hit.ROISummedADC());
+      EInd1.push_back(hit.Integral());
       PTInd1.push_back(PeakTime);
       PAInd1.push_back(hit.PeakAmplitude());
       continue;
@@ -1632,7 +1683,8 @@ void pdvdana::SingleHit::GetListOfTimeCoincidenceHit( bool IsPDVD , bool IsPDHD
 
       WireInd2.push_back(hit.WireID());
       ChannelInd2.push_back(hit.Channel());
-      EInd2.push_back(hit.ROISummedADC());
+      //EInd2.push_back(hit.ROISummedADC());
+      EInd2.push_back(hit.Integral());
       PTInd2.push_back(PeakTime);
       PAInd2.push_back(hit.PeakAmplitude());
     }
@@ -1961,8 +2013,8 @@ std::vector<int> pdvdana::SingleHit::GetXYZIsolatedPoint( std::vector<float> vYP
 
 float pdvdana::SingleHit::dist2(point a, point b)
 {
-    float z = a->z - b->z, y = a->y - b->y;
-    return z*z + y*y;
+    float z = a->z - b->z, y = a->y - b->y , t = a->t - b->t;
+    return z*z + y*y + t*t;
 }
 
 float pdvdana::SingleHit::randf(float m)
@@ -1991,6 +2043,30 @@ point pdvdana::SingleHit::gen_yz(int size , std::vector<int> vIndex , std::vecto
   return pt;
 }
 
+point pdvdana::SingleHit::gen_yzt(int size , std::vector<int> vIndex , std::vector<float> vY , std::vector<float> vZ , std::vector<float> vT , std::vector<int> vNOF , float fElectronVelocity , float fTickToMus)
+{
+  int i = 0;
+  float electronDriftScale = fElectronVelocity * fTickToMus;
+  point p, pt = (point) malloc(sizeof(point_t) * size);
+
+  for (p = pt + size; p-- > pt;)
+  {
+    p->y = vY[vIndex[i]];
+    p->index = vIndex[i];
+    
+    float time_in_cm = vT[vIndex[i]]*electronDriftScale;
+    p->t = time_in_cm;
+    if (vNOF[vIndex[i]] == -1)
+    {
+      p->z = -1*vZ[vIndex[i]] - fgeoZmax;
+    }
+    else p->z = vZ[vIndex[i]];
+    i++;
+  }
+
+  return pt;
+}
+
 void pdvdana::SingleHit::Plot(TCanvas* canvas ,int event , int nbin , float ymin , float ymax , float zmin , float zmax , std::vector<std::vector<float> > &data,std::vector<std::vector<float> > &cluster,double RMS ){
 
     TGraph* grData = new TGraph();
@@ -2070,7 +2146,7 @@ int pdvdana::SingleHit::reallocate(point pt, std::vector<std::vector<float>> Clu
     for( int k = 0 ; k < (int) ClusterPosition[0].size() ; k++) 
     {
 
-      float dist = sqrt(GetDist2D(pt->z,pt->y,ClusterPosition[0][k],ClusterPosition[1][k]));
+      float dist = GetDist(pt->z,pt->y,pt->t,ClusterPosition[0][k],ClusterPosition[1][k],ClusterPosition[2][k]);
       if (min_d > dist ) 
       {
          min_d = dist; 
@@ -2082,13 +2158,13 @@ int pdvdana::SingleHit::reallocate(point pt, std::vector<std::vector<float>> Clu
     return min_i;
 }
 
-
+/*
 float pdvdana::SingleHit::GetDist2D(float y0,float z0,float y1,float z1){
     float z = z0-z1;
     float y = y0-y1;
     return z*z+y*y;
 }
-
+*/
 float pdvdana::SingleHit::mean(float y,float z){
     return (z+y)/2.;
 }
@@ -2134,13 +2210,20 @@ std::vector<std::vector<float>> pdvdana::SingleHit::lloyd(point pts, int len, in
 
     do {
         /* group element for centroids are used as counters */
-        for_n { c->group = 0; c->z = c->y = 0; }
+        for_n { c->group = 0; c->z = c->y = c->t = 0; }
         for_len {
             c = cent + p->group;
             c->group++;
-            c->z += p->z; c->y += p->y;
+            c->z += p->z;
+	    c->y += p->y;
+	    c->t += p->t;
         }
-        for_n { c->z /= c->group; c->y /= c->group; }
+        for_n 
+	{ 
+            c->z /= c->group; 
+	    c->y /= c->group; 
+	    c->t /= c->group;
+	}
 
         changed = 0;
         /* fInd closest centroid of each point */
@@ -2156,16 +2239,18 @@ std::vector<std::vector<float>> pdvdana::SingleHit::lloyd(point pts, int len, in
     for_n { c->group = i; }
 
     std::vector<std::vector<float> > clusterPos;
-    std::vector<float> clusterPosY,clusterPosZ;
+    std::vector<float> clusterPosY,clusterPosZ,clusterPosT;
 
     point result;
 
     for(i = 0, result = cent; i < n_cluster; i++, result++) {
         clusterPosZ.push_back(result->z);
         clusterPosY.push_back(result->y);
+	clusterPosT.push_back(result->t);
     }
     clusterPos.push_back(clusterPosZ);
     clusterPos.push_back(clusterPosY);
+    clusterPos.push_back(clusterPosT);
 
     return clusterPos;
 }
@@ -2174,15 +2259,17 @@ std::vector<std::vector<float>> pdvdana::SingleHit::lloyd(point pts, int len, in
 std::vector<std::vector<float>> pdvdana::SingleHit::GetData(int len,point data){
 
     std::vector<std::vector<float> > dataPos;
-    std::vector<float> dataPosZ,dataPosY;
+    std::vector<float> dataPosZ,dataPosY,dataPosT;
 
 
   for(int i = 0; i < len; i++, data++) {
         dataPosZ.push_back(data->z);
         dataPosY.push_back(data->y);
+	dataPosT.push_back(data->t);
     }
     dataPos.push_back(dataPosZ);
     dataPos.push_back(dataPosY);
+    dataPos.push_back(dataPosT);
 
     return dataPos;
 }
@@ -2203,7 +2290,7 @@ std::vector<int> pdvdana::SingleHit::CheckCompletude(std::vector<std::vector<flo
         {
             if(IDin[j] == 0)
             {
-                dist = sqrt(GetDist2D(data[0][j],data[1][j],cluster[0][i],cluster[1][i]));
+                dist = GetDist(data[0][j],data[1][j],data[2][j],cluster[0][i],cluster[1][i],cluster[2][i]);
                 // printf("Distance to cluster : %f %f \n",i,Nin,Nout);
                 if(dist <= RMS){ Nin++; IDin[j] = 1;}
                 else if(dist <= mult*RMS )
@@ -2246,7 +2333,7 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
     if ( LogLevel > 0) printf("Counting : %.03f %.03f %.03f sum : %.02f \n",float(Nin)/float(Npts),float(Nin2)/float(Npts),float(Nout)/float(Npts),float(Nin+Nin2+Nout)/float(Npts));
 
 
-    if((float(Nin+Nin2)/float(Npts) > tmp -0.04) && float(Nin)/float(Npts) < tmp)
+    if((float(Nin+Nin2)/float(Npts) > tmp) && float(Nin)/float(Npts) < tmp)
     {
       v[0] = 2;
       v[1] = cluster[0].size();
@@ -2271,7 +2358,7 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
         {
             if(j > i)
             {
-                dist = sqrt(GetDist2D(cluster[0][j],cluster[1][j],cluster[0][i],cluster[1][i]));
+                dist = GetDist(cluster[0][j],cluster[1][j],cluster[2][j],cluster[0][i],cluster[1][i],cluster[2][i]);
                 if(dist < 2.*RMS)
                 {
                     IDoverlap[i][j] = 1;
@@ -2282,19 +2369,22 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
     }
     int overlap_counter = 0;
 
-    vector<float> newclusterZ, newclusterY;
-    float meanZ, meanY;
+    vector<float> newclusterZ, newclusterY, newclusterT;
+    float meanZ, meanY, meanT;
     while( (overlap_counter<10)&&(overlap>0) )
     {
       newclusterZ.clear(); 
       newclusterY.clear();
+      newclusterT.clear();
       meanZ = 0;
       meanY = 0;
+      meanT = 0;
 
       for(int i = 0;i<Ncls;i++)
       {
         meanZ = 0.;
         meanY = 0.;
+	meanT = 0.;
         overlap = 0;
 
         for(int j = i;j<Ncls;j++)
@@ -2303,8 +2393,11 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
           {
             meanZ += mean(cluster[0][i],cluster[0][j]);
             meanY += mean(cluster[1][i],cluster[1][j]);
+	    meanT += mean(cluster[2][i],cluster[2][j]);
+
             cluster[0][j] = -999;
             cluster[1][j] = -999;
+	    cluster[2][j] = -999;
             overlap++;
           }
         }
@@ -2312,17 +2405,20 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
         {
           newclusterZ.push_back(cluster[0][i]);
           newclusterY.push_back(cluster[1][i]);
+	  newclusterT.push_back(cluster[2][i]);
         }
         else if(cluster[0][i] != -999)
         {
           newclusterZ.push_back(meanZ/float(overlap));
           newclusterY.push_back(meanY/float(overlap));
+	  newclusterT.push_back(meanT/float(overlap));
         }
       }
 
       cluster.clear();
       cluster.push_back(newclusterZ);
       cluster.push_back(newclusterY);
+      cluster.push_back(newclusterT);
 
       if ( LogLevel > 3) printf("%lu clusters has been removed at iteration %d \n",Ncls-cluster[0].size(),overlap_counter);
 
@@ -2338,8 +2434,8 @@ std::vector<int> pdvdana::SingleHit::CheckClusters(std::vector<std::vector<float
         {
           if(j > i)
           {
-            dist = sqrt(GetDist2D(cluster[0][j],cluster[1][j],cluster[0][i],cluster[1][i]));
-            if(dist < 2.*RMS)
+            dist = GetDist(cluster[0][j],cluster[1][j],cluster[2][j],cluster[0][i],cluster[1][i],cluster[2][i]);
+	    if(dist < 2.*RMS)
             {
               v[j] = 1;
               overlap++;
@@ -2385,7 +2481,8 @@ std::vector<Cluster> pdvdana::SingleHit::GetCluster( int n_point , int n_cluster
                                                      std::vector<int> vMCPDGByEvent , std::vector<int> vMCMOMpdgByEvent ,std::vector<float> vMCWeightByEvent ,  
                                                      std::vector<std::string> vGeneratorTagByEvent ,  
                                                      std::vector<float> vMCXByEvent ,  std::vector<float> vMCYByEvent , std::vector<float> vMCZByEvent , 
-                                                     std::vector<int> vNoFByEvent )
+                                                     std::vector<int> vNoFByEvent,
+                                                     std::vector<float> vMCEnergyByEvent, std::vector<float> vMCElecByEvent )
 {
 
   int k;
@@ -2405,6 +2502,8 @@ std::vector<Cluster> pdvdana::SingleHit::GetCluster( int n_point , int n_cluster
   NullCluster.vMCX.clear();
   NullCluster.vMCY.clear();
   NullCluster.vMCZ.clear();
+  NullCluster.vMCE.clear();
+  NullCluster.vMCElec.clear();
   NullCluster.lChannelCol.clear();
   NullCluster.lChannelInd1.clear();
   NullCluster.EInd1 = 0;
@@ -2438,7 +2537,8 @@ std::vector<Cluster> pdvdana::SingleHit::GetCluster( int n_point , int n_cluster
     float MCPart_x = vMCXByEvent[index];
     float MCPart_y = vMCYByEvent[index];
     float MCPart_z = vMCZByEvent[index];
-
+    float MCE      = vMCEnergyByEvent[index];
+    float MCElec   = vMCElecByEvent[index];
     std::string Generator_tag = vGeneratorTagByEvent[index];
 
     if (ClusterID >=  n_cluster)
@@ -2472,6 +2572,8 @@ std::vector<Cluster> pdvdana::SingleHit::GetCluster( int n_point , int n_cluster
       (vTempCluster[ClusterID].vMCZ).push_back( MCPart_z );
       (vTempCluster[ClusterID].lChannelCol).push_back( ChannelCol );
       (vTempCluster[ClusterID].vNOF).push_back( NoF );
+      (vTempCluster[ClusterID].vMCE).push_back( MCE );
+      (vTempCluster[ClusterID].vMCElec).push_back( MCElec );
     }
     if ( (ChannelInd1 != -1) && (!Inside( ChannelInd1 , vTempCluster[ClusterID].lChannelInd1 ) ) )
     {
@@ -2537,6 +2639,8 @@ std::vector<Cluster> pdvdana::SingleHit::GetCluster( int n_point , int n_cluster
     vCluster[j].vMCX = vTempCluster[j].vMCX;
     vCluster[j].vMCY = vTempCluster[j].vMCY;
     vCluster[j].vMCZ = vTempCluster[j].vMCZ;
+    vCluster[j].vMCE = vTempCluster[j].vMCE;
+    vCluster[j].vMCElec = vTempCluster[j].vMCElec;
   }
 
   if ( LogLevel > 2) std::cout << "WE HAVE " << (float) out/n_point << " POINTS OUTSIDE OF CLUSTERS" << std::endl;
diff --git a/protoduneana/verticaldrift/singlehit/runSingleHit_PDSP.fcl b/protoduneana/verticaldrift/singlehit/runSingleHit_PDSP.fcl
index 899b14a..316ba68 100644
--- a/protoduneana/verticaldrift/singlehit/runSingleHit_PDSP.fcl
+++ b/protoduneana/verticaldrift/singlehit/runSingleHit_PDSP.fcl
@@ -35,7 +35,7 @@ physics:
 {
   analyzers:
   {
-    ana:@Local::SingleHitParams
+    ana:@local::SingleHitParams
   }
 
   analysis: [ ana ]
diff --git a/protoduneana/verticaldrift/singlehit/runSingleHit_PDVD.fcl b/protoduneana/verticaldrift/singlehit/runSingleHit_PDVD.fcl
index 92faa6b..4eb534b 100644
--- a/protoduneana/verticaldrift/singlehit/runSingleHit_PDVD.fcl
+++ b/protoduneana/verticaldrift/singlehit/runSingleHit_PDVD.fcl
@@ -40,7 +40,7 @@ physics:
 {
   analyzers:
   {
-    ana:@Local::SingleHitParams
+    ana:@local::SingleHitParams
   }
 
   analysis: [ ana ]