From cdbdbee478eea7954d94f05e9850ca4252823b59 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Thu, 9 Jul 2020 16:13:41 +0200
Subject: [PATCH] * Adding randomization option before applying calibration    
     - avoid systematic error du to ADC floor behaviour

---
 NPLib/Core/NPCalibrationManager.cxx       | 26 +++++++++++----
 NPLib/Core/NPCalibrationManager.h         |  4 +--
 NPLib/Detectors/Mugast/TMugastPhysics.cxx | 12 +++----
 Projects/MUGAST/Analysis.cxx              | 39 +++++++++++++----------
 Projects/MUGAST/e793.detector             |  6 ++--
 5 files changed, 53 insertions(+), 34 deletions(-)

diff --git a/NPLib/Core/NPCalibrationManager.cxx b/NPLib/Core/NPCalibrationManager.cxx
index b998df204..45864d6e0 100644
--- a/NPLib/Core/NPCalibrationManager.cxx
+++ b/NPLib/Core/NPCalibrationManager.cxx
@@ -24,6 +24,7 @@
 #include "TAsciiFile.h"
 #include "RootOutput.h"
 #include "NPCore.h"
+#include "TRandom.h"
 //   STL
 #include <cstdlib>
 #include <limits>
@@ -199,7 +200,7 @@ void CalibrationManager::LoadParameterFromFile(){
 }
 
 //////////////////////////////////////////////////////////////////
-double CalibrationManager::ApplyCalibration(const std::string& ParameterPath , const double& RawValue) const {
+double CalibrationManager::ApplyCalibration(const std::string& ParameterPath , const double& RawValue, double random) const {
   std::map< std::string , std::vector<double> >::const_iterator it ;
   static std::map< std::string , std::vector<double> >::const_iterator ite = fCalibrationCoeff.end();
 
@@ -208,24 +209,29 @@ double CalibrationManager::ApplyCalibration(const std::string& ParameterPath , c
   it = fCalibrationCoeff.find(ParameterPath)  ;
   // If the find methods return the end iterator it's mean the parameter was not found
   if(it == ite ){
-//by Shuya 170222
-//std::cout << ParameterPath << "!" << std::endl;
     return RawValue ;
   }
 
+  double val ;
+  if(random){
+    val=RawValue + gRandom->Uniform(random);
+    }
+  else
+    val=RawValue;
+
   // The std::vector size give the degree of calibration
   // We just apply the coeff it->second and returned the calibrated value
   double CalibratedValue = 0 ;
   unsigned int mysize = it->second.size(); 
   for(unsigned int i = 0 ; i < mysize ; i++){
-    CalibratedValue += it->second[i]*pow(RawValue, (double)i);
+    CalibratedValue += it->second[i]*pow(val, (double)i);
   }
 
   return CalibratedValue ;
 
 }
 //////////////////////////////////////////////////////////////////
-double CalibrationManager::ApplyCalibrationDebug(const std::string& ParameterPath , const double& RawValue) const{
+double CalibrationManager::ApplyCalibrationDebug(const std::string& ParameterPath , const double& RawValue, double random) const{
   std::map< std::string , std::vector<double> >::const_iterator it ;
   static std::map< std::string , std::vector<double> >::const_iterator ite = fCalibrationCoeff.end();
 
@@ -239,9 +245,17 @@ double CalibrationManager::ApplyCalibrationDebug(const std::string& ParameterPat
     return RawValue ;
   }
 
+  double val ;
+  if(random){
+    val=RawValue + gRandom->Uniform(random);
+    }
+  else
+    val=RawValue;
+
+
   // Else we take the second part of the element (first is index, ie: parameter path)
   // Second is the std::vector of Coeff
-  std::cout << it->first << " :  raw = " << RawValue << " coeff = "  ;
+  std::cout << it->first << " :  raw = " << RawValue << "  randomize = " << val <<  "  coeff = "  ;
   std::vector<double> Coeff = it->second  ;
 
   // The std::vector size give the degree of calibration
diff --git a/NPLib/Core/NPCalibrationManager.h b/NPLib/Core/NPCalibrationManager.h
index 4ea2fd5ce..b1981418f 100644
--- a/NPLib/Core/NPCalibrationManager.h
+++ b/NPLib/Core/NPCalibrationManager.h
@@ -56,10 +56,10 @@ class CalibrationManager{
 
     // call like : myCalibrationManager->ApplyCalibration( "MUST2/Telescope5_Si_X38_E" , RawEnergy )
     // return the Calibrated value
-    double ApplyCalibration (const std::string& ParameterPath , const double& RawValue) const ;
+    double ApplyCalibration (const std::string& ParameterPath , const double& RawValue, double random=0) const ;
     double ApplyResistivePositionCalibration (const std::string& ParameterPath , const double& RawValue) const ;
     // Same but with debug information outputs
-    double ApplyCalibrationDebug (const std::string& ParameterPath , const double& RawValue) const ;
+    double ApplyCalibrationDebug (const std::string& ParameterPath , const double& RawValue, double random=0) const ;
     double ApplyResistivePositionCalibrationDebug (const std::string& ParameterPath , const double& RawValue) const ;
 
     bool ApplyThreshold (const std::string& ParameterPath, const double& RawValue) const ;
diff --git a/NPLib/Detectors/Mugast/TMugastPhysics.cxx b/NPLib/Detectors/Mugast/TMugastPhysics.cxx
index fb0690af6..137669a5a 100644
--- a/NPLib/Detectors/Mugast/TMugastPhysics.cxx
+++ b/NPLib/Detectors/Mugast/TMugastPhysics.cxx
@@ -1040,7 +1040,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetDSSDXEStripNbr(i));
     name += "_E";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetDSSDXEEnergy(i));
+        name, m_EventData->GetDSSDXEEnergy(i),1);
   }
 
   double fDSSD_X_T(const TMugastData* m_EventData, const int& i) {
@@ -1051,7 +1051,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetDSSDXTStripNbr(i));
     name += "_T";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetDSSDXTTime(i));
+        name, m_EventData->GetDSSDXTTime(i),1);
   }
 
   //   Y
@@ -1063,7 +1063,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetDSSDYEStripNbr(i));
     name += "_E";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetDSSDYEEnergy(i));
+        name, m_EventData->GetDSSDYEEnergy(i),1);
   }
 
   double fDSSD_Y_T(const TMugastData* m_EventData, const int& i) {
@@ -1074,7 +1074,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetDSSDYTStripNbr(i));
     name += "_T";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetDSSDYTTime(i));
+        name, m_EventData->GetDSSDYTTime(i),1);
   }
 
   //   SecondLayer
@@ -1086,7 +1086,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetSecondLayerEStripNbr(i));
     name += "_E";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetSecondLayerEEnergy(i));
+        name, m_EventData->GetSecondLayerEEnergy(i),1);
   }
 
   double fSecondLayer_T(const TMugastData* m_EventData, const int& i) {
@@ -1097,7 +1097,7 @@ namespace MUGAST_LOCAL {
     name += NPL::itoa(m_EventData->GetSecondLayerTStripNbr(i));
     name += "_T";
     return CalibrationManager::getInstance()->ApplyCalibration(
-        name, m_EventData->GetSecondLayerTTime(i));
+        name, m_EventData->GetSecondLayerTTime(i),1);
   }
 }
 
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index fbe025ba0..acbb74d52 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -109,6 +109,7 @@ void Analysis::Init() {
   Z = 0;
   dE = 0;
   BeamDirection = TVector3(0,0,1);
+  nbTrack=0;
 
 }
 
@@ -280,6 +281,7 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D");
   RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
   RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
+  if(!simulation){
   // Vamos 
   RootOutput::getInstance()->GetTree()->Branch("LTS",&LTS,"LTS/l");
 
@@ -302,7 +304,8 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("coreTS",coreTS,"coreTS[nbCores]/l");
   RootOutput::getInstance()->GetTree()->Branch("coreE0",coreE0,"coreE0[nbCores]/F");
   //
-  if(simulation){
+  }
+  else{
     RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
     RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
   }
@@ -312,22 +315,24 @@ void Analysis::InitOutputBranch() {
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitInputBranch(){
   // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
-  // Vamos
-  RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",&LTS);
-  // Agata
-  RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("nbCores",&nbCores);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
-  if(simulation){
+  if(!simulation){
+    // Vamos
+    RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",&LTS);
+    // Agata
+    RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("nbCores",&nbCores);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
+    RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
+  }
+  else{
     RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
     RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
     RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
diff --git a/Projects/MUGAST/e793.detector b/Projects/MUGAST/e793.detector
index bb0bad911..178f0b0ba 100644
--- a/Projects/MUGAST/e793.detector
+++ b/Projects/MUGAST/e793.detector
@@ -1,14 +1,14 @@
 %%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 Target
- THICKNESS= 20 micrometer
+ THICKNESS= 0.1 micrometer
  ANGLE= 0 deg
  RADIUS= 10 mm
- MATERIAL= CD2
+ MATERIAL= Vacuum
  X= 0 mm
  Y= 0 mm
  Z= 0 mm
- NbSlices= 100
+ NbSlices= 10
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CATSDetector
-- 
GitLab