Merge pull request #128 from ejangelico/Application

first version NNLS reco algorithm and improvement to LAPPDParseACC

Author: brichards64

DataModel/NnlsSolution.cpp

@@ -0,0 +1,31 @@
+#include "NnlsSolution.h"
+
+NnlsSolution::~NnlsSolution()
+{
+	component_times.clear();
+	component_scales.clear();
+	templateWaveform.ClearSamples(); //template used in the nnls algo, at the time binning used in the algorithm
+	templateTimes.clear(); //time series of length of templateWaveform
+	fullNnlsWaveform.ClearSamples(); //the full solution, sum of all components
+}
+
+//add an nnls component to the components vector. 
+//each of these components is a copy of the template waveform
+//at time "t" and scaled with factor "s". 
+void NnlsSolution::AddComponent(double t, double s)
+{
+	component_times.push_back(t);
+	component_scales.push_back(s);
+}
+
+//sum of all components, full waveform fit
+void NnlsSolution::SetFullSoln(Waveform<double> fullsoln)
+{
+	fullNnlsWaveform.SetSamples(fullsoln.Samples());
+}
+
+void NnlsSolution::SetTemplate(Waveform<double> tmpwfm, vector<double> temptimes)
+{
+	templateWaveform.SetSamples(tmpwfm.Samples());
+	templateTimes = temptimes;
+}

DataModel/NnlsSolution.h

@@ -0,0 +1,76 @@
+/*
+This class is a data structure to store the already solved 
+nnls algorithm outputs. It is meant to reduce the memory overhead
+used by the nnls matrices/vectors and just store info on 
+the full solution and components of the solution. Originally, 
+components of the solution were stored as full waveforms, but
+this would crash due to memory overusage. This class stores
+components as a tuple of time index and scale factor applied to template. 
+-Evan Angelico, 7/11/2019
+*/
+
+
+#ifndef NnlsSolution_H
+#define NnlsSolution_H
+
+#include <SerialisableObject.h>
+#include "Waveform.h"
+
+using namespace std;
+
+class NnlsSolution : public SerialisableObject
+{
+	friend class boost::serialization::access;
+
+public:
+	NnlsSolution() : component_times(vector<double>{}), component_scales(vector<double>{}), templateWaveform(Waveform<double>{}),templateTimes(vector<double>{}),fullNnlsWaveform(Waveform<double>{}) {serialise=true;}
+	~NnlsSolution();
+	void AddComponent(double t, double s);
+	void SetFullSoln(Waveform<double> fullsoln); //sum of all components, full waveform fit
+	void SetTemplate(Waveform<double> tmpwfm, vector<double> temptimes);
+	inline int GetNumberOfComponents() {return component_times.size();} //how many components are populated in components vector
+	inline double GetComponentTime(int index) {return component_times.at(index);}
+	inline double GetComponentScale(int index) {return component_scales.at(index);}
+
+	bool Print() {
+		int verbose=0;
+		bool filled = false;
+		if(component_times.size() != 0) filled = true;
+		cout << "Solution filled?: " << filled << endl;
+		if(filled) cout << "Number of components: " << component_times.size() << endl;
+		if(verbose){
+			cout<<"Times and scales:" << endl;
+			for(int i = 0; i < component_times.size(); i++){
+				cout << component_times.at(i) << ", " << component_scales.at(i) << endl;
+			}
+		}
+
+		return true;
+	}
+
+	
+
+protected:
+	//components of the solution are pairs of doubles
+	//first element is the time, second element is the scale
+	//at which a template waveform appears relative to the
+	//raw waveform in the time units of the raw waveform
+	vector<double> component_times;
+	vector<double> component_scales;
+	Waveform<double> templateWaveform; //template used in the nnls algo, at the time binning used in the algorithm
+	vector<double> templateTimes; //time series of length of templateWaveform
+	Waveform<double> fullNnlsWaveform; //the full solution, sum of all components
+
+	template<class Archive> void serialize(Archive & ar, const unsigned int version){
+		if(serialise){
+			ar & component_times;
+			ar & component_scales;
+			ar & templateWaveform;
+			ar & templateTimes;
+			ar & fullNnlsWaveform;
+		}
+	}
+};
+
+
+#endif

UserTools/Factory/Factory.cpp

@@ -61,6 +61,7 @@ if (tool=="VtxExtendedVertexFinder") ret=new VtxExtendedVertexFinder;
 if (tool=="VtxPointDirectionFinder") ret=new VtxPointDirectionFinder;
 if (tool=="VtxPointVertexFinder") ret=new VtxPointVertexFinder;
 if (tool=="LoadCCData") ret=new LoadCCData;
+if (tool=="WaveformNNLS") ret=new WaveformNNLS;
 if (tool=="HitCleaner") ret=new HitCleaner;
 if (tool=="HitResiduals") ret=new HitResiduals;
 if (tool=="MonitorReceive") ret=new MonitorReceive;

UserTools/GenerateHits/GenerateHits.cpp

@@ -1,4 +1,6 @@
 #include "GenerateHits.h"
+#include <cstdlib>
+#include <time.h>
 
 GenerateHits::GenerateHits():Tool(){}
 
@@ -27,36 +29,36 @@ bool GenerateHits::Execute(){
   std::vector<double> relpos;
   double tc,Q;
 
-  // MC truth hit 1
-  Q=1.;
-  pos.push_back(10.);
-  pos.push_back(5.);
-  pos.push_back(0.);
-  relpos.push_back(10.);
-  relpos.push_back(5.);
-  tc = 1.050;
-  LAPPDHit ahit1(TubeID,tc,Q,pos,relpos);
-
-  // clear vectors
-  pos.clear();
-  relpos.clear();
-
-/*
-  // MC truth hit 2
-  tc = 1.140;
-  Q=1.;
-  pos.push_back(20.);
-  pos.push_back(10.);
-  pos.push_back(0.);
-  relpos.push_back(20.);
-  relpos.push_back(10.);
-  LAPPDHit ahit2(TubeID,tc,Q,pos,relpos);
-*/
+  int nhits;
+  m_variables.Get("nhits", nhits);
 
+  //initialize random seed
+  srand(time(NULL));
+
+  //all the hits
   vector<LAPPDHit> hits;
-  hits.push_back(ahit1);
-//  hits.push_back(ahit2);
 
+  //generate uniformly distributed hits
+  //accross the surface of the LAPPD with
+  //random values of position charge and time
+  for(int i = 0; i < nhits; i++)
+  {
+    Q=1; 
+    //abs pos is not used for my analysis at the moment
+    pos.push_back(0.);
+    pos.push_back(0.);
+    pos.push_back(0.);
+    relpos.push_back(fRand(-98, 98));
+    relpos.push_back(fRand(-98, 98));
+    tc = fRand(0, 24.9); //ns, converted to ps later in LAPPDSim processing
+    LAPPDHit temphit(TubeID,tc,Q,pos,relpos);
+    pos.clear();
+    relpos.clear();
+    hits.push_back(temphit);
+
+  }
+
+  cout << "inserting " << hits.size() << " synthetic hits " << endl;
   // stuff the two hits into MCLAPPHit
   MCLAPPDHit.insert(pair <int,vector<LAPPDHit>> (0,hits));
 
@@ -67,7 +69,14 @@ bool GenerateHits::Execute(){
 }
 
 
+double GenerateHits::fRand(double fMin, double fMax)
+{
+    double f = (double)rand() / RAND_MAX;
+    return fMin + f * (fMax - fMin);
+}
+
 bool GenerateHits::Finalise(){
+  cout << "Finalized generate alright " << endl;
 
   return true;
 }

UserTools/GenerateHits/GenerateHits.h

@@ -16,6 +16,7 @@ class GenerateHits: public Tool {
   bool Initialise(std::string configfile,DataModel &data);
   bool Execute();
   bool Finalise();
+  double fRand(double fMin, double fMax);
 
 
  private:

UserTools/LAPPDAnalysis/LAPPDAnalysis.cpp

@@ -23,8 +23,9 @@ bool LAPPDAnalysis::Execute(){
   //BaselineSubtract
     Waveform<double> bwav;
     // get raw lappd data
-    std::map<int,vector<Waveform<double>>> rawlappddata;
+    map<int,vector<Waveform<double>>> rawlappddata;
     m_data->Stores["ANNIEEvent"]->Get("RawLAPPDData",rawlappddata);
+    
 
     // the filtered Waveform
     std::map<int,vector<Waveform<double>>> blsublappddata;