Alpha ana

.gitlab-ci.yml

@@ -26,6 +26,9 @@ stages:
   # Stage to validate restG4 package examples
   - restG4
 
+  # Stage to validate framework examples
+  - examples
+
   # Stage to validate data chain processing constructions (TREX-DM data, PandaX-III topology, etc).
   - restManager_process
 
@@ -201,6 +204,32 @@ listRESTMacros:
       - ${CI_PROJECT_DIR}/pipeline/pandaxiii_MC/Xe136bb0n_n2E06.root
     expire_in: 1 day
 
+08_alphas:
+  type: restG4
+  script:
+    - . ${CI_PROJECT_DIR}/install_latest/thisREST.sh
+    - cd ${CI_PROJECT_DIR}/install_latest/examples/restG4/08.Alphas/
+    - mkdir data
+    - export REST_ENERGY=5
+    - export REST_FOIL=1
+    - restG4 alphas.rml
+  artifacts:
+    paths:
+      - ${CI_PROJECT_DIR}/install_latest/examples/restG4/08.Alphas/data/Run_5MeV_1um.root
+    expire_in: 1 day
+
+01_alphaTrack:
+    type: examples
+    script:
+    - . ${CI_PROJECT_DIR}/install_latest/thisREST.sh
+    - cd ${CI_PROJECT_DIR}/install_latest/examples/01.alphaTrack/
+    - mkdir data
+    - export RAW_PROCESSING=ON
+    - export HIT_PROCESSING=ON
+    - export TRACK_PROCESSING=ON
+    - restManager --c processing.rml --f ${CI_PROJECT_DIR}/install_latest/examples/restG4/08.Alphas/data/Run_5MeV_1um.root
+    - restRoot -b -q Validate.C'("data/Response_5MeV_1um.root")'
+
 LoadGas:
   type: metadata
   script:

.gitmodules

@@ -64,3 +64,6 @@
 [submodule "doc/userguide"]
 	path = doc/userguide
 	url = https://github.com/rest-for-physics/rest-docs.github.io.git
+[submodule "projects/alphacamm"]
+	path = projects/alphacamm
+	url = git@lfna.unizar.es:gifna/alphacamm.git

cmake/FindGarfieldOld.cmake

@@ -15,14 +15,17 @@
 set(Garfield_DIR $ENV{GARFIELD_HOME})
 
 find_path(Garfield_INCLUDE_DIRS Sensor.hh
-        HINTS ${Garfield_DIR}/include/ ${Garfield_DIR}/include/Garfield ${Garfield_INC_DIR}
-        ${Garfield_DIR}/Include)
+        HINTS ${Garfield_DIR}/Include/ ${Garfield_DIR}/Include/Garfield
+        ${Garfield_DIR}/include/ ${Garfield_DIR}/include/Garfield)
+
+find_path(Garfield_INCLUDE_Heed_DIRS HeedChamber.hh
+        HINTS ${Garfield_DIR}/Heed/)
 
 #message(STATUS Garfield_INCLUDE_DIRS ${Garfield_INCLUDE_DIRS})
 
 find_library(Garfield_LIBRARIES NAMES libGarfield.so Garfield
         HINTS ${Garfield_DIR}/lib ${Garfield_LIB_DIR}
-        $ENV{GARFIELD_HOME}/lib)
+        $ENV{GARFIELD_DIR}/lib)
 #message(STATUS Garfield_LIBRARIES ${Garfield_LIBRARIES})
 
 #if (${Garfield_LIBRARY_DIR})

examples/01.alphaTrack/README.md

@@ -2,7 +2,7 @@ This example shows how to include the detector response and the reconstruction o
 
 The processes included have been adapted to facilitate the reconstruction of short energetic alpha tracks. Before using the example it is necessary to generate the Geant4 event MC data. You may use the following restG4 [example](https://github.com/rest-for-physics/restG4/examples/08.Alphas).
 
-A quick recipe to generate the necessary data.
+A quick recipe to generate the necessary data follows.
 
 Copy this example at your working directory. Then copy the restG4 generator example.
 

examples/01.alphaTrack/Validate.C

@@ -0,0 +1,12 @@
+Int_t Validate(string fname) {
+    TRestRun* run = new TRestRun(fname);
+    cout << "Run entries : " << run->GetEntries() << endl;
+    if (run->GetEntries() < 4350 || run->GetEntries() > 4500) {
+        cout << "The number of entries is not in the range (4350,4500)!" << endl;
+        return 1;
+    }
+    delete run;
+
+    cout << "[\033[92m OK \x1b[0m]" << endl;
+    return 0;
+}

examples/01.alphaTrack/processing.rml

@@ -101,12 +101,16 @@
 		<addProcess type="TRestRawSignalAnalysisProcess" name="rawAna" title="Signal noise" value="${RAW_PROCESSING}" observable="all">
 		</addProcess>
 
-		<addProcess type="TRestRawZeroSuppresionProcess" name="zS" value="${RAW_PROCESSING}"
-			baseLineRange="(5,55)"
-			integralRange="(10,500)"
-			pointThreshold="2"
-			pointsOverThreshold="5"
-			signalThreshold="3" />
+		<addProcess type="TRestRawToDetectorSignalProcess" name="zS" value="${RAW_PROCESSING}"
+                  triggerStarts="0us"
+                  gain="1"
+                  threshold="0.01"
+                  zeroSuppression="true"
+		  baseLineRange="(5,55)"
+		  integralRange="(10,500)"
+		  pointThreshold="2"
+		  nPointsOverThreshold="5"
+		  signalThreshold="3" />
 
 			<!-- Requires a bit of debugging 
 		<addProcess type="TRestDetectorSignalChannelActivityProcess" name="chActivity" value="${RAW_PROCESSING}" verboseLevel="info" observable="all" >
@@ -122,6 +126,9 @@
 			<parameter name="clusterDistance" value="1" units="cm" />
 		</addProcess>
 
+		<addProcess type="TRestDetectorHitsToTrackProcess" name="hitsToTrack" title="" value="${TRACK_PROCESSING}" 
+			clusterDistance="2.5mm" />
+
 		<addProcess type="TRestTrackReductionProcess" name="reduction" title="" value="${TRACK_PROCESSING}">
 			<parameter name="startingDistance" value="0.15mm" />
 			<parameter name="minimumDistance" value="0.5mm" />

pipeline/trex/02_signal.rml

@@ -13,14 +13,18 @@
 
 		<parameter name="eventsToProcess" value="0" />
 
-		<addProcess type="TRestRawZeroSuppresionProcess" name="zS" value="ON"
-			baseLineRange="(${BL_MIN},${BL_MAX})"
-			integralRange="(${INT_MIN},${INT_MAX})"
-			pointThreshold="${POINT_TH}"
-			pointsOverThreshold="${NPOINTS}"
-			sampling="${SAMPLING}"
-			signalThreshold="${SGNL_TH}"
-			verboseLevel="silent"	/>
+		<addProcess type="TRestRawToDetectorSignalProcess" name="detSignal" value="ON"
+		       threshold="0.01"
+		       triggerStarts="0us"
+		       gain="1"
+		       zeroSuppression="true"
+		       baseLineRange="(${BL_MIN},${BL_MAX})"
+		       integralRange="(${INT_MIN},${INT_MAX})"
+		       pointThreshold="${POINT_TH}"
+		       pointsOverThreshold="${NPOINTS}"
+		       sampling="${SAMPLING}"
+		       signalThreshold="${SGNL_TH}"
+		       verboseLevel="silent" />
 
 		<addProcess type="TRestDetectorSignalChannelActivityProcess" name="chActivity" value="ON" verboseLevel="info" observable="all" >
 		    <parameter name="daqChannels" value="72*4*4" />

source/framework/core/inc/TRestAnalysisPlot.h

@@ -70,6 +70,7 @@ class TRestAnalysisPlot : public TRestMetadata {
         string title;
 
         Double_t normalize;
+        string scale;
         Bool_t logX;
         Bool_t logY;
         Bool_t logZ;

source/framework/core/inc/TRestHits.h

@@ -194,6 +194,8 @@ class TRestHits : public TObject {
 
     void WriteHitsToTextFile(TString filename);
 
+    static void GetBoundaries(std::vector<double> &val, double &max, double &min, int &nBins, double offset=10);
+
     virtual void PrintHits(Int_t nHits = -1);
 
     class TRestHits_Iterator : public std::iterator<std::random_access_iterator_tag, TRestHits_Iterator> {

source/framework/core/inc/TRestVolumeHits.h

@@ -67,6 +67,8 @@ class TRestVolumeHits : public TRestHits {
     }
     Double_t GetXYSize(int n) { return TMath::Sqrt(fSigmaX[n] * fSigmaX[n] + fSigmaY[n] * fSigmaY[n]); }
 
+    static void kMeansClustering(TRestVolumeHits *hits, TRestVolumeHits &vHits, int maxIt =100);
+
     // Constructor & Destructor
     TRestVolumeHits();
     ~TRestVolumeHits();

source/framework/core/src/TRestAnalysisPlot.cxx

@@ -247,6 +247,7 @@ void TRestAnalysisPlot::InitFromConfigFile() {
             plot.gridY = StringToBool(GetParameter("gridY", plotele, "false"));
             plot.gridX = StringToBool(GetParameter("gridX", plotele, "false"));
             plot.normalize = StringToDouble(GetParameter("norm", plotele, ""));
+            plot.scale = GetParameter("scale", plotele, "");
             plot.labelX = GetParameter("xlabel", plotele, "");
             plot.labelY = GetParameter("ylabel", plotele, "");
             plot.ticksX = StringToInteger(GetParameter("xticks", plotele, "510"));
@@ -843,7 +844,7 @@ void TRestAnalysisPlot::PlotCombinedCanvas() {
             continue;
         }
 
-        // scale the histograms
+        // normalize the histograms
         if (plot.normalize > 0) {
             for (unsigned int i = 0; i < plot.histos.size(); i++) {
                 if (plot.histos[i].ptr == nullptr) continue;
@@ -855,6 +856,18 @@ void TRestAnalysisPlot::PlotCombinedCanvas() {
             }
         }
 
+       // scale the histograms
+       if (plot.scale != "") {
+            for (unsigned int i = 0; i < plot.histos.size(); i++) {
+                if (plot.histos[i].ptr == nullptr) continue;
+                Double_t scale = 1.;
+                if(plot.scale == "binSize")scale =  1./plot.histos[i].ptr->GetXaxis()->GetBinWidth(1);
+                else scale=StringToDouble(plot.scale);
+
+                plot.histos[i].ptr->Scale(scale);
+            }
+        }
+
         // draw to the pad
         targetPad = (TPad*)fCombinedCanvas->cd(n + 1 + fPanels.size());
         Double_t maxValue_Pad = 0;

source/framework/core/src/TRestHits.cxx

@@ -558,6 +558,26 @@ void TRestHits::WriteHitsToTextFile(TString filename) {
     fclose(fff);
 }
 
+void TRestHits::GetBoundaries(std::vector<double> &dist, double &max, double &min, int &nBins, double offset){
+
+  std::sort(dist.begin(),dist.end());
+  max=dist.back();
+  min=dist.front();
+
+  double minDiff = 1E6;
+  double prevVal=1E6;
+   for (const auto &h : dist){
+     double diff = std::abs( h - prevVal);
+     if(diff>0 && diff<minDiff )minDiff=diff;
+     prevVal=h;
+   }
+
+  max += offset*minDiff + minDiff/2.;
+  min -= offset*minDiff + minDiff/2.;
+  nBins = std::round((max-min)/minDiff);
+
+}
+
 Double_t TRestHits::GetGaussSigmaX() {
 	Double_t gausSigmaX = 0;
 	Int_t nHits = GetNumberOfHits();

source/framework/core/src/TRestVolumeHits.cxx

@@ -163,3 +163,54 @@ void TRestVolumeHits::PrintHits() {
              << " Energy: " << GetEnergy(n) << endl;
     }
 }
+
+
+//______________________________________________________________________________
+void TRestVolumeHits::kMeansClustering(TRestVolumeHits *hits, TRestVolumeHits &vHits, int maxIt){
+
+        const int nodes = vHits.GetNumberOfHits();
+        TRestVolumeHits volHits[nodes];
+        //std::cout<<"Nhits "<<hits->GetNumberOfHits()<<" Nodes "<<nodes<<std::endl;
+        TVector3 nullVector = TVector3(0,0,0);
+        std::vector<TVector3> centroid(nodes);
+        std::vector<TVector3> centroidOld(nodes,nullVector);
+
+          for(int h=0;h<nodes; h++)centroid[h] = vHits.GetPosition(h);
+
+           for(int it=0;it<maxIt;it++){
+             for(int n=0;n<nodes;n++)volHits[n].RemoveHits();
+               for (int i = 0; i < hits->GetNumberOfHits(); i++){
+                 double minDist=1E9;
+                 int clIndex=-1;
+                   for(int n=0;n<nodes;n++){
+                     TVector3 hitPos =hits->GetPosition(i);
+                     double dist = (centroid[n] - hitPos).Mag();
+                       if(dist < minDist){
+                         minDist=dist;
+                         clIndex=n;
+                       }
+                   }
+                 //cout<<minDist<<" "<<clIndex<<endl;
+                 volHits[clIndex].AddHit(*hits,i);
+               }
+              bool converge=true;
+                for(int n=0;n<nodes;n++){
+                  centroid[n] = volHits[n].GetMeanPosition();
+                  converge &= (centroid[n] == centroidOld[n]);
+                  centroidOld[n] = centroid[n];
+                }
+              if(converge){
+              break;
+              }
+           }
+
+       vHits.RemoveHits();
+       const TVector3 sigma(0., 0., 0.);
+         for(int n=0;n<nodes;n++){
+             if(volHits[n].GetNumberOfHits()>0)
+             vHits.AddHit(volHits[n].GetMeanPosition(),volHits[n].GetTotalEnergy(),0,volHits[n].GetType(0),sigma);
+         }
+
+}
+
+