rest-for-physics · juanangp · Apr 26, 2022 · Jan 31, 2022 · Jan 31, 2022 · Jan 31, 2022
@@ -28,6 +28,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
 
@@ -148,6 +151,32 @@ List REST Macros:
       - ${CI_PROJECT_DIR}/pipeline/pandaxiii_MC/Xe136bb0n_n2E06.root
     expire_in: 1 day
 
+08_alphas:
+  type: restG4
+  script:
+    - source ${CI_PROJECT_DIR}/install/thisREST.sh
+    - cd ${CI_PROJECT_DIR}/install/examples/restG4/08.Alphas/
+    - mkdir data
+    - export REST_ENERGY=5
+    - export REST_FOIL=1
+    - restG4 alphas.rml
+  artifacts:
+    paths:
+      - ${CI_PROJECT_DIR}/install/examples/restG4/08.Alphas/data/Run_5MeV_1um.root
+    expire_in: 1 day
+
+01_alphaTrack:
+    type: examples
+    script:
+    - source ${CI_PROJECT_DIR}/install/thisREST.sh
+    - cd ${CI_PROJECT_DIR}/install/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/examples/restG4/08.Alphas/data/Run_5MeV_1um.root
+    - restRoot -b -q Validate.C'("data/Response_5MeV_1um.root")'
+
 Load Gas:
   stage: metadata
   script:

@@ -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
@@ -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.
 

@@ -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;
+}
@@ -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" >
@@ -118,6 +122,12 @@
 		<addProcess type="TRestDetectorSignalToHitsProcess" name="signalToHits" value="${HIT_PROCESSING}"
 			observable="all" method="all" />
 
+                <addProcess type="TRestDetectorHitsToTrackFastProcess" name="hitsToTrackFast" title="" value="OFF">
+                       <parameter name="netOrigin" value="(-100,-100,-100)" units="mm" />
+                       <parameter name="netSize" value="200" units="mm" />
+                       <parameter name="cellResolution" value="2" units="mm" />
+                </addProcess>
+
 		<addProcess type="TRestDetectorHitsToTrackProcess" name="hitsToTrack" title="" value="${TRACK_PROCESSING}">
 			<parameter name="clusterDistance" value="1" units="cm" />
 		</addProcess>

@@ -38,7 +38,7 @@
       <parameter name="pointThreshold" value="3" />
       <parameter name="pointsOverThreshold" value="3" />
     </addProcess>
-    <addProcess type="TRestRawZeroSuppresionProcess" name="zS" value="ON" file="processes.rml"/>
+    <addProcess type="TRestRawToDetectorSignalProcess" name="zS" value="ON" file="processes.rml"/>
 
   </TRestProcessRunner>
 

@@ -439,18 +439,19 @@
         <observable name="Qratio_Y_R30" value="ON" />
         <parameter name="hitsToCheckFraction" value="0.25" />
     </TRestTrackBlobAnalysisProcess>
-    <TRestRawZeroSuppresionProcess name="zS" title="" verboseLevel="silent">
-        <!-- This parameter is used to define the baseline calculation -->
-        <parameter name="baseLineRange" value="(50,120)" />
+    <TRestRawToDetectorSignalProcess name="zS" title="" verboseLevel="silent">
         <!-- These parameters define the integral with threshold. 
          threshold : number of baseline noise sigmas to consider a point for integration.
          pointsOverThreshold : Number of consecutive points over threshold to be considered for integration.
          minPeakAmplitude : Minimum peak signal amplitude to be considered at all.
         -->
+        <!-- This parameter is used to define the baseline calculation -->
+        <parameter name="baseLineRange" value="(50,120)" />
+        <parameter name="zeroSuppression" value="true" />
         <parameter name="pointThreshold" value="1.9" />
-        <parameter name="pointsOverThreshold" value="10" />
+        <parameter name="nPointsOverThreshold" value="10" />
         <parameter name="signalThreshold" value="2.8" />
-    </TRestRawZeroSuppresionProcess>
+    </TRestRawToDetectorSignalProcess>
     <TRestRawSignalAnalysisProcess name="sAna" title="" verboseLevel="silent">
         <!-- This parameter defines every how many seconds a signal will be drawn.
          In case the value is zero. There will not be any drawing. -->
@@ -465,7 +466,7 @@
          minPeakAmplitude : Minimum peak signal amplitude to be considered at all.
         -->
         <parameter name="pointThreshold" value="1.9" />
-        <parameter name="pointsOverThreshold" value="10" />
+        <parameter name="nPointsOverThreshold" value="10" />
         <parameter name="signalThreshold" value="2.8" />
         <!-- ++++++++++++++++++++++
          Observable definitions

@@ -31,7 +31,7 @@
 
     <addProcess type="TRestRawMultiCoBoAsAdToSignalProcess" name="virtualDAQ" value="ON"/>
     <addProcess type="TRestRawSignalAnalysisProcess" name="sAna" value="ON">
-
+      <parameter name="pointsOverThreshold" value="3" />
       /// We define all observables except MinValue because it is not yet in validation chain
       <observable name="BaseLineMean" value="ON" />
       <observable name="BaseLineSigmaMean" value="ON" />
@@ -56,8 +56,12 @@
 
     </addProcess>
     <addProcess type="TRestEventRateAnalysisProcess" name="rateAna" value="ON" observable="all"/>
-    <addProcess type="TRestRawZeroSuppresionProcess" name="zS" value="ON" sampling="0.2us"/> 5MHz sampling
-
+    <addProcess type="TRestRawToDetectorSignalProcess" name="rawSignalToSignal" value="ON">
+       <parameter name="sampling" value="0.2us"/>
+       <parameter name="zeroSuppression" value="true"/>
+       <parameter name="nPointsOverThreshold" value="3" />
+    </addProcess>
+
     <addProcess type="TRestRealTimeDrawingProcess" name="rD" value="ON" drawInterval="1000" >
       <TRestAnalysisPlot name="TriggerRatePlot"  previewPlot="false">
         <canvas size="(1000,800)" save="TriggerRate.png"/>
@@ -78,7 +82,6 @@
     <parameter name="mainDataPath" value="" />
 
     <parameter name="pointThreshold" value="3" />
-    <parameter name="pointsOverThreshold" value="3" />
   </globals>
 
 </TRestManager>

@@ -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}"
+		       nPointsOverThreshold="${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" />

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

@@ -195,6 +195,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) const;
 
     class TRestHits_Iterator : public std::iterator<std::random_access_iterator_tag, TRestHits_Iterator> {

@@ -69,6 +69,8 @@ class TRestVolumeHits : public TRestHits {
         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();

@@ -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;

@@ -562,6 +562,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();

@@ -161,3 +161,54 @@ void TRestVolumeHits::PrintHits() const {
              << " 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);
+         }
+
+}
+
+