From ad2b482ad838c5385b1ca1256749604eab6b0cc5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?GIRARD=20ALCINDOR=20Val=C3=A9rian?=
<girard-alcindor@ijclab.in2p3.fr>
Date: Mon, 27 Jan 2025 11:14:03 +0100
Subject: [PATCH] nptuto_gdr
---
Projects/nptuto_gdr/.MugastShape | 5 +
Projects/nptuto_gdr/.RunToTreat.txt | 3 +
Projects/nptuto_gdr/.ls_return | 1 +
Projects/nptuto_gdr/Analysis.cxx | 242 ++++++++++++++++++
Projects/nptuto_gdr/Analysis.h | 114 +++++++++
Projects/nptuto_gdr/CMakeLists.txt | 5 +
Projects/nptuto_gdr/configs/ConfigMust2.dat | 5 +
.../nptuto_gdr/detector/MUGAST_LISE.detector | 56 ++++
Projects/nptuto_gdr/reactions/48Crpd.reaction | 32 +++
9 files changed, 463 insertions(+)
create mode 100644 Projects/nptuto_gdr/.MugastShape
create mode 100644 Projects/nptuto_gdr/.RunToTreat.txt
create mode 100644 Projects/nptuto_gdr/.ls_return
create mode 100644 Projects/nptuto_gdr/Analysis.cxx
create mode 100644 Projects/nptuto_gdr/Analysis.h
create mode 100644 Projects/nptuto_gdr/CMakeLists.txt
create mode 100644 Projects/nptuto_gdr/configs/ConfigMust2.dat
create mode 100644 Projects/nptuto_gdr/detector/MUGAST_LISE.detector
create mode 100644 Projects/nptuto_gdr/reactions/48Crpd.reaction
diff --git a/Projects/nptuto_gdr/.MugastShape b/Projects/nptuto_gdr/.MugastShape
new file mode 100644
index 000000000..cf4a3011b
--- /dev/null
+++ b/Projects/nptuto_gdr/.MugastShape
@@ -0,0 +1,5 @@
+1 1
+3 1
+4 1
+5 1
+7 1
diff --git a/Projects/nptuto_gdr/.RunToTreat.txt b/Projects/nptuto_gdr/.RunToTreat.txt
new file mode 100644
index 000000000..4f1adfe89
--- /dev/null
+++ b/Projects/nptuto_gdr/.RunToTreat.txt
@@ -0,0 +1,3 @@
+Tree SimulatedTree
+ /Users/valerian/Documents/Work/Conferences/GDR_Resanet/nptuto_v3/Software/nptool/Outputs/Simulation/SimuNptuto.root
+
diff --git a/Projects/nptuto_gdr/.ls_return b/Projects/nptuto_gdr/.ls_return
new file mode 100644
index 000000000..c5c1f557b
--- /dev/null
+++ b/Projects/nptuto_gdr/.ls_return
@@ -0,0 +1 @@
+/Users/valerian/Documents/Work/Conferences/GDR_Resanet/nptuto_v3/Software/nptool/Outputs/Simulation/SimuNptuto.root
diff --git a/Projects/nptuto_gdr/Analysis.cxx b/Projects/nptuto_gdr/Analysis.cxx
new file mode 100644
index 000000000..c71f7c79d
--- /dev/null
+++ b/Projects/nptuto_gdr/Analysis.cxx
@@ -0,0 +1,242 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include <iostream>
+using namespace std;
+
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPFunction.h"
+#include "NPOptionManager.h"
+// #include <unistd.h>
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis() {}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis() {}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init() {
+ if (NPOptionManager::getInstance()->HasDefinition("simulation")) {
+ cout << "Considering input data as simulation" << endl;
+ simulation = true;
+ }
+ else {
+ cout << "Considering input data as real" << endl;
+ simulation = false;
+ }
+ simulation = true;
+
+ // initialize input and output branches
+ if (simulation) {
+ Initial = new TInitialConditions();
+ ReactionConditions = new TReactionConditions();
+ }
+
+ InitOutputBranch();
+ InitInputBranch();
+ // get MUST2 objects
+ M2 = (TMust2Physics*)m_DetectorManager->GetDetector("M2Telescope");
+
+ // get reaction information
+ reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ OriginalBeamEnergy = reaction.GetBeamEnergy();
+ // target thickness
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+
+ // energy losses
+ string light = NPL::ChangeNameToG4Standard(reaction.GetNucleus3()->GetName());
+ string beam = NPL::ChangeNameToG4Standard(reaction.GetNucleus1()->GetName());
+ cout << light << " " << beam << " " << TargetMaterial << " " << TargetThickness << endl;
+ LightTarget = NPL::EnergyLoss(light + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+ LightAl = NPL::EnergyLoss(light + "_Al.G4table", "G4Table", 100);
+ LightSi = NPL::EnergyLoss(light + "_Si.G4table", "G4Table", 100);
+ BeamTarget = NPL::EnergyLoss(beam + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+
+ FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness * 0.5, 0);
+ // FinalBeamEnergy = OriginalBeamEnergy;
+ // cout << "Original beam energy: " << OriginalBeamEnergy << " MeV Mid-target beam energy: " << FinalBeamEnergy
+ // << "MeV " << endl;
+ reaction.SetBeamEnergy(FinalBeamEnergy);
+
+ // initialize various parameters
+ Rand = TRandom3();
+ DetectorNumber = 0;
+ ThetaNormalTarget = 0;
+ ThetaM2Surface = 0;
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ Energy = 0;
+ ThetaGDSurface = 0;
+ X = 0;
+ Y = 0;
+ Z = 0;
+ dE = 0;
+ BeamDirection = TVector3(0, 0, 1);
+ nbTrack = 0;
+
+ sleep(5);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent() {
+ // Reinitiate calculated variable
+ ReInitValue();
+ double XTarget, YTarget;
+ TVector3 BeamDirection;
+ XTarget = 0;
+ YTarget = 0;
+ BeamDirection = TVector3(0, 0, 1);
+ BeamImpact = TVector3(XTarget, YTarget, 0);
+
+ // determine beam energy for a randomized interaction point in target
+ // 1% FWHM randominastion (E/100)/2.35
+ // reaction.SetBeamEnergy(Rand.Gaus(BeamEnergy, BeamEnergy * 1. / 100. / 2.35));
+ // reaction.SetBeamEnergy(BeamEnergy);
+
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////// LOOP on MUST2 ////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ for (unsigned int countMust2 = 0; countMust2 < M2->Si_E.size(); countMust2++) {
+ /************************************************/
+ // Part 0 : Get the usefull Data
+ // MUST2
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ TVector3 HitDirection = M2->GetPositionOfInteraction(countMust2) - BeamImpact;
+ ThetaLab = HitDirection.Angle(BeamDirection);
+
+ X = M2->GetPositionOfInteraction(countMust2).X();
+ Y = M2->GetPositionOfInteraction(countMust2).Y();
+ Z = M2->GetPositionOfInteraction(countMust2).Z();
+
+ ThetaM2Surface = HitDirection.Angle(-M2->GetTelescopeNormal(countMust2));
+ ThetaNormalTarget = HitDirection.Angle(TVector3(0, 0, 1));
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2 = M2->CsI_E[countMust2];
+
+ // if CsI
+ if (CsI_E_M2 > 0) {
+ // The energy in CsI is calculate form dE/dx Table because
+ Energy = CsI_E_M2;
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
+ Energy += Si_E_M2;
+ }
+ else {
+ Energy = Si_E_M2;
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
+ }
+
+ // Evaluate energy using the thickness
+ ELab = Energy;
+
+ // Target Correction
+ ELab = LightTarget.EvaluateInitialEnergy(ELab, TargetThickness * 0.5, ThetaNormalTarget);
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ // Ex = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // reaction.SetBeamEnergy(Initial->GetIncidentFinalKineticEnergy());
+ Ex = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // ExNoBeam = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // reaction.SetBeamEnergy(FinalBeamEnergy);
+ // ExNoProton = reaction.ReconstructRelativistic(ReactionConditions->GetKineticEnergy(0),
+ // ReactionConditions->GetParticleDirection(0).Angle(TVector3(0, 0,
+ // 1)));
+ ThetaLab = ThetaLab / deg;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = reaction.EnergyLabToThetaCM(ELab, ThetaLab) / deg;
+ /************************************************/
+ } // end loop MUST2
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End() {}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex, "Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab, "ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab, "ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("Run", &Run, "Run/I");
+ RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D");
+ RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D");
+ RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+ RootOutput::getInstance()->GetTree()->Branch("dE", &dE, "dE/D");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch() {
+ // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
+ if (!simulation) {
+ }
+ else {
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fIC_*", true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions", &Initial);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*", true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions);
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue() {
+ Ex = -1000;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+ X = -1000;
+ Y = -1000;
+ Z = -1000;
+ dE = -1000;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the AnalysisFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct() { return (NPL::VAnalysis*)new Analysis(); }
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+class proxy_analysis {
+ public:
+ proxy_analysis() { NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); }
+};
+
+proxy_analysis p_analysis;
+}
diff --git a/Projects/nptuto_gdr/Analysis.h b/Projects/nptuto_gdr/Analysis.h
new file mode 100644
index 000000000..f314fb353
--- /dev/null
+++ b/Projects/nptuto_gdr/Analysis.h
@@ -0,0 +1,114 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "NPVAnalysis.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TCATSPhysics.h"
+#include "TInitialConditions.h"
+#include "TMust2Physics.h"
+#include "TReactionConditions.h"
+#include <TMath.h>
+#include <TRandom3.h>
+#include <TVector3.h>
+
+class Analysis : public NPL::VAnalysis {
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+
+ double TargetThickness;
+ NPL::Reaction reaction;
+ // Energy loss table: the G4Table are generated by the simulation
+ NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightAl;
+ NPL::EnergyLoss LightSi;
+ NPL::EnergyLoss BeamTarget;
+
+ // Beam Energy
+ double OriginalBeamEnergy; // AMEV
+ double FinalBeamEnergy;
+
+ // intermediate variable
+ TVector3 BeamDirection;
+ TVector3 BeamImpact;
+ TRandom3 Rand;
+ int Run;
+ int DetectorNumber;
+ double ThetaNormalTarget;
+ double ThetaM2Surface;
+ double Si_E_M2;
+ double CsI_E_M2;
+ double Energy;
+ double BeamEnergy;
+
+ double ThetaGDSurface;
+ double X;
+ double Y;
+ double Z;
+ // Vamos Branches
+ unsigned long long int LTS;
+
+ // Agata branches
+ double agata_zShift;
+ unsigned long long int TStrack;
+ int nbHits;
+ int nbTrack;
+ float* trackE = new float(100);
+ float* trackX1 = new float(100);
+ float* trackY1 = new float(100);
+ float* trackZ1 = new float(100);
+ float* trackT = new float(100);
+ int* trackCrystalID = new int(100);
+ int nbCores;
+ int* coreId = new int(100);
+ ULong64_t* coreTS = new ULong64_t(100);
+ float* coreE0 = new float(100);
+ //
+ double dE;
+ double dTheta;
+ // Branches and detectors
+ TMust2Physics* M2;
+ bool simulation;
+ TInitialConditions* Initial;
+ TReactionConditions* ReactionConditions;
+};
+#endif
diff --git a/Projects/nptuto_gdr/CMakeLists.txt b/Projects/nptuto_gdr/CMakeLists.txt
new file mode 100644
index 000000000..22c74affd
--- /dev/null
+++ b/Projects/nptuto_gdr/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required (VERSION 2.8)
+# Setting the policy to match Cmake version
+cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
+# include the default NPAnalysis cmake file
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/nptuto_gdr/configs/ConfigMust2.dat b/Projects/nptuto_gdr/configs/ConfigMust2.dat
new file mode 100644
index 000000000..216b16b66
--- /dev/null
+++ b/Projects/nptuto_gdr/configs/ConfigMust2.dat
@@ -0,0 +1,5 @@
+ConfigMust2
+ CSI_SIZE 40
+ CAL_PIXEL
+ PIXEL_SIZE 8
+ MAX_STRIP_MULTIPLICITY 10
diff --git a/Projects/nptuto_gdr/detector/MUGAST_LISE.detector b/Projects/nptuto_gdr/detector/MUGAST_LISE.detector
new file mode 100644
index 000000000..f03211b19
--- /dev/null
+++ b/Projects/nptuto_gdr/detector/MUGAST_LISE.detector
@@ -0,0 +1,56 @@
+%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Target
+ THICKNESS= 50 micrometer
+ ANGLE= 0 deg
+ RADIUS= 25 mm
+ MATERIAL= CH2
+ X= 0
+ Y= 0
+ Z= 0
+
+%%%%%%% Telescope 1 %%%%%%%
+M2Telescope
+ X128_Y128= 115.88 9.61 154.54 mm
+ X128_Y1= 104.8 101.89 125.09 mm
+ X1_Y1= 14.55 102.4 160.63 mm
+ X1_Y128= 25.63 10.12 190.08 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 2 %%%%%%%
+M2Telescope
+ X128_Y128= -11.23 102.42 160.87 mm
+ X128_Y1= -101.39 102.39 124.37 mm
+ X1_Y1= -113.17 10.36 153.56 mm
+ X1_Y128= -23.03 10.38 190.05 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 3 %%%%%%%
+M2Telescope
+ X128_Y128= -113.28 -12.52 153.32 mm
+ X128_Y1= -101.58 -104.77 124.82 mm
+ X1_Y1= -11.39 -104.58 161.48 mm
+ X1_Y128= -23.1 -12.34 189.98 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X128_Y128= 13.82 -104.92 160.72 mm
+ X128_Y1= 104.3 -104.95 125.08 mm
+ X1_Y1= 115.75 -12.73 153.76 mm
+ X1_Y128= 25.23 -12.65 189.43 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/nptuto_gdr/reactions/48Crpd.reaction b/Projects/nptuto_gdr/reactions/48Crpd.reaction
new file mode 100644
index 000000000..ef5c42de5
--- /dev/null
+++ b/Projects/nptuto_gdr/reactions/48Crpd.reaction
@@ -0,0 +1,32 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 48Cr
+ Energy= 1440 MeV
+ SigmaEnergy= 40 MeV
+ ExcitationEnergy= 0 MeV
+ SigmaThetaX= 0 mrad
+ SigmaPhiY= 0 mrad
+ SigmaX= 1 mm
+ SigmaY= 1 mm
+ MeanThetaX= 0 mrad
+ MeanPhiY= 0 mrad
+ MeanX= 0 mm
+ MeanY= 0 mm
+ %EnergyProfilePath= eLise.txt EL
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 48Cr
+ Target= 1H
+ Light= 2H
+ Heavy= 47Cr
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= flat.txt
+ ShootLight= 1
+ ShootHeavy= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
--
GitLab