diff --git a/.gitignore b/.gitignore
index fdf84a2f46f5e76c147af7587628006b773148d9..52f234c603db5b906dcc961e3c0ff7cf240b2657 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,3 +40,8 @@ Documentation/user_guide.log
*.SRIM
*.sublime-*
PhysicsListOption.txt
+Projects/T40/Analysis.cxx
+Projects/T40/Analysis.h
+Projects/T40/T40.detector
+Projects/T40/19Fdp.reaction
+Projects/T40/configs/*.dat
diff --git a/Inputs/DetectorConfiguration/Tiara.detector b/Inputs/DetectorConfiguration/Tiara.detector
index af6a05d500fb6539a648f859900405cda2d4d998..3cb2707f9374b3df477332860801eb72ff449956 100644
--- a/Inputs/DetectorConfiguration/Tiara.detector
+++ b/Inputs/DetectorConfiguration/Tiara.detector
@@ -14,9 +14,9 @@ Target
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tiara
- TiaraInnerBarrel= 1
- TiaraOuterBarrel= 1
- TiaraChamber= 1
+ InnerBarrel= 1
+ OuterBarrel= 1
+ Chamber= 1
TiaraHyballWedge
Z= -147
R= 0
@@ -42,5 +42,3 @@ Tiara
R= 0
Phi= 300
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
diff --git a/NPLib/Detectors/AnnularCsI/TAnnularCsIPhysics.cxx b/NPLib/Detectors/AnnularCsI/TAnnularCsIPhysics.cxx
index 38d63673691ea4afab062499a8cbf4333422f50a..fffee0a21a7a851615ea1ed1f961469f9cfe8e23 100644
--- a/NPLib/Detectors/AnnularCsI/TAnnularCsIPhysics.cxx
+++ b/NPLib/Detectors/AnnularCsI/TAnnularCsIPhysics.cxx
@@ -36,6 +36,7 @@ using namespace std;
#include "RootOutput.h"
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"
+#include "NPSystemOfUnits.h"
// ROOT
#include "TChain.h"
@@ -194,14 +195,17 @@ int TAnnularCsIPhysics::MatchToSi(const TVector3& SiPosition) const {
// returns -1 if no detector found
// todo account for offset between detectors
- int i = 0;
- for (const auto& g : m_WedgeGeometry){
- if(SiPosition.Phi() >= g.Phi_min && SiPosition.Phi() <= g.Phi_max &&
- SiPosition.Pt() >= g.R_min && SiPosition.Pt() <= g.R_max) {
- // match!
- return i;
+ int hitno = 0;
+ for(const auto& detno : DetectorNumber){
+ const AnnularCsI_Utils::Geometry& g = m_WedgeGeometry.at(detno - 1);
+ if(SiPosition.Phi() >= g.Phi_min && SiPosition.Phi() < g.Phi_max &&
+ SiPosition.Pt() >= g.R_min && SiPosition.Pt() < g.R_max &&
+ fabs(SiPosition.Z()/NPUNITS::mm - g.Z/NPUNITS::mm) < 10) {
+ //
+ // Match!
+ return hitno;
}
- ++i;
+ ++hitno;
}
// no match
return -1;
diff --git a/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.cxx b/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c80adce16ebe4c77b6e3636849666d066dc9f694
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.cxx
@@ -0,0 +1,107 @@
+// STL
+#include <vector>
+#include <string>
+using namespace std;
+
+// NPL
+#include "NPSystemOfUnits.h"
+#include "NPOptionManager.h"
+#include "NPInputParser.h"
+using namespace NPUNITS;
+
+// ROOT
+#include <TMath.h>
+
+#include "AnnularTelescope_Utils.h"
+
+
+namespace { void do_print(AnnularTelescope_Utils::Geometry& g) {
+ std::cout
+ << "\tZ POSITION [mm]: " << g.Z / mm << "\n"
+ << "\tRADIUS [mm]: " << g.R_min / mm << " - " << g.R_max / mm << "\n"
+ << "\tCSI_THICKNESS [mm]: " << g.CsIThickness / mm << "\n"
+ << "\tCSI_WEDGES [deg]:\n";
+ for(size_t i=0; i< g.CsI_Wedge_Phi_Angle.size(); ++i){
+ cout <<
+ "\t\t" << i << " : " << (g.CsI_Wedge_Phi_Angle[i] - g.CsI_Wedge_Angle_Pitch/2.) / deg
+ << " - " << (g.CsI_Wedge_Phi_Angle[i] + g.CsI_Wedge_Angle_Pitch/2.) / deg << "\n";
+ }
+ std::cout
+ << "\tSI_THICKNESS [mm]: " << g.SiThickness / mm << "\n"
+ << "\tSI_PHI_STRIPS [deg]:\n";
+ for(size_t i=0; i< g.Si_Strip_Phi_Angle.size(); ++i){
+ cout <<
+ "\t\t" << i << " : " << (g.Si_Strip_Phi_Angle[i] - g.Si_Phi_Angle_Pitch/2.)/deg
+ << " - " << (g.Si_Strip_Phi_Angle[i] + g.Si_Phi_Angle_Pitch/2.)/deg << "\n";
+ }
+ std::cout << "\tSI_THETA_STRIPS [mm]:\n";
+ for(size_t i=0; i< g.Si_Strip_Theta_Radius.size(); ++i){
+ cout <<
+ "\t\t" << i << " : " << (g.Si_Strip_Theta_Radius[i] - g.Si_Theta_Radius_Pitch/2.)/mm
+ << " - " << (g.Si_Strip_Theta_Radius[i] + g.Si_Theta_Radius_Pitch/2.)/mm << "\n";
+ }
+ std::cout << "////////////////////////////////////////////////////////" << std::endl;
+} }
+
+std::vector<AnnularTelescope_Utils::Geometry>
+AnnularTelescope_Utils::ReadConfiguration(NPL::InputParser& parser){
+ vector<Geometry> out;
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("AnnularTelescope");
+ if(NPOptionManager::getInstance()->GetVerboseLevel()){
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+ }
+ vector<string> wedge = {
+ "Z", "R_MIN", "R_MAX",
+ "CSI_THICKNESS", "CSI_WEDGES",
+ "SI_THICKNESS", "SI_THETA_STRIPS", "SI_PHI_STRIPS"
+ };
+
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(wedge)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel()){
+ cout << endl << "//// AnnularTelescope " << i+1 << endl;
+ }
+ Geometry g;
+ g.R_min = blocks[i]->GetDouble("R_MIN", "mm");
+ g.R_max = blocks[i]->GetDouble("R_MAX", "mm");
+
+ int n_wedges = blocks[i]->GetInt("CSI_WEDGES");
+ g.CsI_Wedge_Angle_Pitch = (2*TMath::Pi()) / n_wedges;
+ for(int i=0; i< n_wedges; ++i){
+ g.CsI_Wedge_Phi_Angle.push_back(
+ -180*NPUNITS::deg + i*g.CsI_Wedge_Angle_Pitch + g.CsI_Wedge_Angle_Pitch/2. );
+ }
+
+ int n_theta = blocks[i]->GetInt("SI_THETA_STRIPS");
+ g.Si_Theta_Radius_Pitch = (g.R_max - g.R_min) / n_theta;
+ for(int i=0; i< n_theta; ++i){
+ g.Si_Strip_Theta_Radius.push_back(
+ g.R_min + i*g.Si_Theta_Radius_Pitch + g.Si_Theta_Radius_Pitch/2. );
+ }
+
+ int n_phi = blocks[i]->GetInt("SI_PHI_STRIPS");
+ g.Si_Phi_Angle_Pitch = (2*TMath::Pi()) / n_phi;
+ for(int i=0; i< n_phi; ++i){
+ g.Si_Strip_Phi_Angle.push_back(
+ -180*NPUNITS::deg + i*g.Si_Phi_Angle_Pitch + g.Si_Phi_Angle_Pitch/2. );
+ }
+
+ g.Z = blocks[i]->GetDouble("Z", "mm");
+ g.SiThickness = blocks[i]->GetDouble("SI_THICKNESS", "um");
+ g.CsIThickness = blocks[i]->GetDouble("CSI_THICKNESS", "mm");
+
+ out.push_back(g);
+ // if(NPOptionManager::getInstance()->GetVerboseLevel()){
+ // do_print(g);
+ // }
+ }
+ else{
+ cout << "ERROR (AnnularTelescope): "
+ << "check your input file formatting." << endl
+ << "Here is a dump of the problem block: " << endl;
+ blocks[i]->Dump();
+ exit(1);
+ }
+ }
+ return out;
+}
diff --git a/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.h b/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..9af0c1e650eef48c0643c4d56da1480a89f108cf
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/AnnularTelescope_Utils.h
@@ -0,0 +1,36 @@
+/// \file Define some utility functions and structs for the annular telescope
+/// detector
+#ifndef ANNULAR_TELESCOPE_UTILS_HEADER
+#define ANNULAR_TELESCOPE_UTILS_HEADER
+
+// STL
+#include <vector>
+// NPL
+namespace NPL { class InputParser; }
+
+namespace AnnularTelescope_Utils {
+
+/// Sescribes the geometry of an annular telescope
+struct Geometry{
+ double R_min; /// inner radius (mm)
+ double R_max; /// outer radius (mm)
+ double CsI_Wedge_Angle_Pitch; /// ANGULAR pitch of CsI wedges (rad)
+ double Si_Theta_Radius_Pitch; /// RADIAL pitch of Si theta strips (mm)
+ double Si_Phi_Angle_Pitch; /// ANGULAR pitch of Si Phi strips (rad)
+ std::vector<double> CsI_Wedge_Phi_Angle; /// CsI wedge phi central angles (rad)
+ std::vector<double> Si_Strip_Theta_Radius; /// Si Theta strips RADIUS (central) (mm)
+ std::vector<double> Si_Strip_Phi_Angle; /// Si Phi strips ANGLE (central) (rad)
+ double Z; /// z-position of detector (mm relative to target)
+ double SiThickness; /// Si detector thickness (mm)
+ double CsIThickness; /// CsI detector thickness (mm)
+};
+
+/// Read detector geometry from a config file
+std::vector<Geometry> ReadConfiguration(NPL::InputParser& parser);
+
+}
+
+#endif
+/* Local Variables: */
+/* mode: c++ */
+/* End: */
diff --git a/NPLib/Detectors/AnnularTelescope/CMakeLists.txt b/NPLib/Detectors/AnnularTelescope/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a84cdcfe3437b3ba1fec3a476cc3fb3585cd498e
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_custom_command(OUTPUT TAnnularTelescopePhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TAnnularTelescopePhysics.h TAnnularTelescopePhysicsDict.cxx TAnnularTelescopePhysics.rootmap libNPAnnularTelescope.dylib DEPENDS TAnnularTelescopePhysics.h)
+add_custom_command(OUTPUT TAnnularTelescopeDataDict.cxx COMMAND ../../scripts/build_dict.sh TAnnularTelescopeData.h TAnnularTelescopeDataDict.cxx TAnnularTelescopeData.rootmap libNPAnnularTelescope.dylib DEPENDS TAnnularTelescopeData.h)
+add_library(NPAnnularTelescope SHARED TAnnularTelescopeSpectra.cxx TAnnularTelescopeData.cxx TAnnularTelescopePhysics.cxx TAnnularTelescopeDataDict.cxx TAnnularTelescopePhysicsDict.cxx AnnularTelescope_Utils.cxx)
+target_link_libraries(NPAnnularTelescope ${ROOT_LIBRARIES} NPCore)
+install(FILES TAnnularTelescopeData.h TAnnularTelescopePhysics.h TAnnularTelescopeSpectra.h AnnularTelescope_Utils.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.cxx b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e4425d45f35fbd404144ab802ac3fe2366a8f8c8
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.cxx
@@ -0,0 +1,145 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include "TAnnularTelescopeData.h"
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+using namespace std;
+
+ClassImp(TAnnularTelescopeData)
+
+
+//////////////////////////////////////////////////////////////////////
+TAnnularTelescopeData::TAnnularTelescopeData() {
+ Clear();
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+TAnnularTelescopeData::~TAnnularTelescopeData() {
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeData::Clear() {
+ CsI_E_Detector.clear();
+ CsI_E_Wedge.clear();
+ CsI_E_Energy.clear();
+
+ CsI_T_Detector.clear();
+ CsI_T_Wedge.clear();
+ CsI_T_Time.clear();
+
+ Si_Theta_E_Detector.clear();
+ Si_Theta_E_Strip.clear();
+ Si_Theta_E_Energy.clear();
+
+ Si_Theta_T_Detector.clear();
+ Si_Theta_T_Strip.clear();
+ Si_Theta_T_Time.clear();
+
+ Si_Phi_E_Detector.clear();
+ Si_Phi_E_Strip.clear();
+ Si_Phi_E_Energy.clear();
+
+ Si_Phi_T_Detector.clear();
+ Si_Phi_T_Strip.clear();
+ Si_Phi_T_Time.clear();
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeData::Dump() const {
+ // This method is very useful for debuging and worth the dev.
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TAnnularTelescopeData::Dump()] XXXXXXXXXXXXXXXXX" << endl;
+
+ // CsI Energy
+ size_t mysize = CsI_E_Detector.size();
+ cout << "AnnularTelescope_E_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "DetNbr: " << CsI_E_Detector[i]
+ << " CsI WedgeNbr: " << CsI_E_Wedge[i]
+ << " CsI Energy: " << CsI_E_Energy[i];
+ }
+ cout << "\n";
+
+ // CsI Time
+ mysize = CsI_T_Detector.size();
+ cout << "AnnularTelescope_T_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "CsI DetNbr: " << CsI_T_Detector[i]
+ << " CsI WedgeNbr: " << CsI_T_Detector[i]
+ << " CsI Time: " << CsI_T_Time[i];
+ }
+ cout << "\n";
+
+ // Si Energy
+ mysize = Si_Theta_E_Detector.size();
+ cout << "Si:: AnnularTelescope_ThetaE_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "DetNbr: " << Si_Theta_E_Detector[i]
+ << " Si ThetaStripNbr: " << Si_Theta_E_Strip[i]
+ << " Si ThetaStripEnergy: " << Si_Theta_E_Energy[i];
+ }
+ cout << "\n";
+
+ // Si Time
+ mysize = Si_Theta_T_Detector.size();
+ cout << "Si:: AnnularTelescope_ThetaT_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "DetNbr: " << Si_Theta_T_Detector[i]
+ << " Si ThetaStripNbr: " << Si_Theta_T_Strip[i]
+ << " Si ThetaStripTime: " << Si_Theta_T_Time[i];
+ }
+ cout << "\n";
+
+ // Si Energy (phi)
+ mysize = Si_Phi_E_Detector.size();
+ cout << "Si:: AnnularTelescope_PhiE_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "DetNbr: " << Si_Phi_E_Detector[i]
+ << " Si PhiStripNbr: " << Si_Phi_E_Strip[i]
+ << " Si PhiStripEnergy: " << Si_Phi_E_Energy[i];
+ }
+ cout << "\n";
+
+ // Si Time (phi)
+ mysize = Si_Phi_T_Detector.size();
+ cout << "Si:: AnnularTelescope_PhiT_Mult: " << mysize << endl;
+
+ for (size_t i = 0 ; i < mysize ; i++){
+ cout << "DetNbr: " << Si_Phi_T_Detector[i]
+ << " Si PhiStripNbr: " << Si_Phi_T_Strip[i]
+ << " Si PhiStripTime: " << Si_Phi_T_Time[i];
+ }
+ cout << "\n";
+}
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.h b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.h
new file mode 100644
index 0000000000000000000000000000000000000000..24699a0bf83be1c6b0b007261b3695518ec13b72
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeData.h
@@ -0,0 +1,207 @@
+#ifndef AnnularTelescopeDATA_Include_Guard
+#define AnnularTelescopeDATA_Include_Guard
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <vector>
+#include <iostream>
+
+// ROOT
+#include "TObject.h"
+
+class TAnnularTelescopeData : public TObject {
+ //////////////////////////////////////////////////////////////
+ // data members are hold into vectors in order
+ // to allow multiplicity treatment
+private:
+// CsI Energy
+ std::vector<UShort_t> CsI_E_Detector; // detector number
+ std::vector<UShort_t> CsI_E_Wedge; // wedge number
+ std::vector<Double_t> CsI_E_Energy; // value (energy or time)
+// CsI Time
+ std::vector<UShort_t> CsI_T_Detector; // detector number
+ std::vector<UShort_t> CsI_T_Wedge; // wedge number
+ std::vector<Double_t> CsI_T_Time; // value (energy or time)
+// Si Energy (theta strips)
+ std::vector<UShort_t> Si_Theta_E_Detector; // detector number
+ std::vector<UShort_t> Si_Theta_E_Strip; // radial strip number
+ std::vector<Double_t> Si_Theta_E_Energy; // value (energy or time)
+// Si Time (theta strips)
+ std::vector<UShort_t> Si_Theta_T_Detector; // detector number
+ std::vector<UShort_t> Si_Theta_T_Strip; // radial strip number
+ std::vector<Double_t> Si_Theta_T_Time; // value (energy or time)
+// Si Energy (phi strips)
+ std::vector<UShort_t> Si_Phi_E_Detector; // detector number
+ std::vector<UShort_t> Si_Phi_E_Strip; // radial strip number
+ std::vector<Double_t> Si_Phi_E_Energy; // value (energy or time)
+// Si Time (phi strips)
+ std::vector<UShort_t> Si_Phi_T_Detector; // detector number
+ std::vector<UShort_t> Si_Phi_T_Strip; // radial strip number
+ std::vector<Double_t> Si_Phi_T_Time; // value (energy or time)
+
+ //////////////////////////////////////////////////////////////
+ // Constructor and destructor
+public:
+ TAnnularTelescopeData();
+ ~TAnnularTelescopeData();
+
+
+ //////////////////////////////////////////////////////////////
+ // Inherited from TObject and overriden to avoid warnings
+public:
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+
+ //////////////////////////////////////////////////////////////
+ // Getters and Setters
+ // Prefer inline declaration to avoid unnecessary called of
+ // frequently used methods
+ // add //! to avoid ROOT creating dictionnary for the methods
+public:
+ ////////////////////// SETTERS ////////////////////////
+ // CsI Energy
+ void SetCsIEnergy(UShort_t DetNbr, UShort_t WedgeNbr, Double_t Energy){
+ CsI_E_Detector.push_back(DetNbr);
+ CsI_E_Wedge.push_back(WedgeNbr);
+ CsI_E_Energy.push_back(Energy);
+ };//!
+
+ // CsI Time
+ void SetCsITime(UShort_t DetNbr, UShort_t WedgeNbr, Double_t Time) {
+ CsI_T_Detector.push_back(DetNbr);
+ CsI_T_Wedge.push_back(WedgeNbr);
+ CsI_T_Time.push_back(Time);
+ };//!
+
+ // Si Energy (theta strips)
+ void SetSiThetaEnergy(UShort_t DetNbr, UShort_t StripNbr, Double_t Energy) {
+ Si_Theta_E_Detector.push_back(DetNbr);
+ Si_Theta_E_Strip.push_back(StripNbr);
+ Si_Theta_E_Energy.push_back(Energy);
+ }
+
+ // Si Time (theta strips)
+ void SetSiThetaTime(UShort_t DetNbr, UShort_t StripNbr, Double_t Time) {
+ Si_Theta_T_Detector.push_back(DetNbr);
+ Si_Theta_T_Strip.push_back(StripNbr);
+ Si_Theta_T_Time.push_back(Time);
+ }
+
+ // Si Energy (phi strips)
+ void SetSiPhiEnergy(UShort_t DetNbr, UShort_t StripNbr, Double_t Energy) {
+ Si_Phi_E_Detector.push_back(DetNbr);
+ Si_Phi_E_Strip.push_back(StripNbr);
+ Si_Phi_E_Energy.push_back(Energy);
+ }
+
+ // Si Time (phi strips)
+ void SetSiPhiTime(UShort_t DetNbr, UShort_t StripNbr, Double_t Time) {
+ Si_Phi_T_Detector.push_back(DetNbr);
+ Si_Phi_T_Strip.push_back(StripNbr);
+ Si_Phi_T_Time.push_back(Time);
+ }
+
+ ////////////////////// GETTERS ////////////////////////
+ /////////////////////////////////////////////////////////////
+ // Helper function to catch errors better
+ template <class T>
+ T GetArrayVal(const std::vector<T>& v, size_t i, const std::string& name) const {
+ try { return v.at(i); }
+ catch(std::exception& e) {
+ std::cerr << "ERROR in AnnularTelescopeData, can't get value of "
+ << name << " at index " << i << std::endl;
+ exit(1);
+ }
+ } //!
+
+ // CsI Energy
+ UShort_t GetCsIMultEnergy() const
+ {return CsI_E_Detector.size();}
+ UShort_t GetCsIE_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(CsI_E_Detector,i,"CsI_E_Detector");}//!
+ UShort_t GetCsIE_WedgeNbr(const unsigned int &i) const
+ {return GetArrayVal(CsI_E_Wedge,i,"CsI_E_Wedge");}//!
+ Double_t GetCsIEnergy(const unsigned int &i) const
+ {return GetArrayVal(CsI_E_Energy,i,"CsI_E_Energy");}//!
+
+ // CsI Time
+ UShort_t GetCsIMultTime() const
+ {return CsI_T_Detector.size();}
+ UShort_t GetCsIT_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(CsI_T_Detector,i,"CsI_T_Detector");}//!
+ UShort_t GetCsIT_WedgeNbr(const unsigned int &i) const
+ {return GetArrayVal(CsI_T_Wedge,i,"CsI_T_Wedge");}//!
+ Double_t GetCsITime(const unsigned int &i) const
+ {return GetArrayVal(CsI_T_Time,i,"CsI_T_Time");}//!
+
+ // Si Energy (theta strips)
+ UShort_t GetSiThetaE_Mult() const
+ {return Si_Theta_E_Detector.size();}
+ UShort_t GetSiThetaE_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_E_Detector,i,"Si_Theta_E_Detector");}//!
+ UShort_t GetSiThetaE_StripNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_E_Strip,i,"Si_Theta_E_Strip");}//!
+ Double_t GetSiThetaE_Energy(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_E_Energy,i,"Si_Theta_E_Energy");}//!
+
+ // Si Time (theta strips)
+ UShort_t GetSiThetaT_Mult() const
+ {return Si_Theta_T_Detector.size();}
+ UShort_t GetSiThetaT_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_T_Detector,i,"Si_Theta_T_Detecor");}//!
+ UShort_t GetSiThetaT_StripNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_T_Strip,i,"Si_Theta_T_Strip");}//!
+ Double_t GetSiThetaT_Time(const unsigned int &i) const
+ {return GetArrayVal(Si_Theta_T_Time,i,"Si_Theta_T_Time");}//!
+
+ // Si Energy (phi strips)
+ UShort_t GetSiPhiE_Mult() const
+ {return Si_Phi_E_Detector.size();}
+ UShort_t GetSiPhiE_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_E_Detector,i,"Si_Phi_E_Detector");}//!
+ UShort_t GetSiPhiE_StripNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_E_Strip,i,"Si_Phi_E_Strip");}//!
+ Double_t GetSiPhiE_Energy(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_E_Energy,i,"Si_Phi_E_Energy");}//!
+
+ // Si Time (phi strips)
+ UShort_t GetSiPhiT_Mult() const
+ {return Si_Phi_T_Detector.size();}
+ UShort_t GetSiPhiT_DetectorNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_T_Detector,i,"Si_Phi_T_Detecor");}//!
+ UShort_t GetSiPhiT_StripNbr(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_T_Strip,i,"Si_Phi_T_Strip");}//!
+ Double_t GetSiPhiT_Time(const unsigned int &i) const
+ {return GetArrayVal(Si_Phi_T_Time,i,"Si_Phi_T_Time");}//!
+
+ //////////////////////////////////////////////////////////////
+ // Required for ROOT dictionnary
+ ClassDef(TAnnularTelescopeData,1) // AnnularTelescopeData structure
+};
+
+#endif
+/* Local Variables: */
+/* mode: c++ */
+/* End: */
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.cxx b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7e16d4633e4d771ef444b7c6dc4f8764d0c447f9
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.cxx
@@ -0,0 +1,625 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "TAnnularTelescopePhysics.h"
+
+// STL
+#include <cassert>
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+using namespace std;
+
+// NPL
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "NPDetectorFactory.h"
+#include "NPOptionManager.h"
+#include "NPSystemOfUnits.h"
+using namespace NPUNITS;
+
+// ROOT
+#include "TChain.h"
+#include "TRandom3.h"
+
+ClassImp(TAnnularTelescopePhysics)
+
+
+///////////////////////////////////////////////////////////////////////////
+TAnnularTelescopePhysics::TAnnularTelescopePhysics()
+ : m_EventData(new TAnnularTelescopeData),
+ m_PreTreatedData(new TAnnularTelescopeData),
+ m_EventPhysics(this),
+ m_Spectra(0)
+{
+ m_EnergyMatch.NSigma = 5;
+ m_EnergyMatch.Sigma = 0.0149;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::AddDetector(const AnnularTelescope_Utils::Geometry& g){
+ m_WedgeGeometry.push_back(g);
+}
+
+///////////////////////////////////////////////////////////////////////////
+TVector3 TAnnularTelescopePhysics::GetPositionOfInteraction(int i) const {
+ try {
+ const Int_t detector = m_SiHits.at(i).Detector;
+ const Int_t theta_n = m_SiHits.at(i).ThetaStrip;
+ const Int_t phi_n = m_SiHits.at(i).PhiStrip;
+ const AnnularTelescope_Utils::Geometry& geo =
+ m_WedgeGeometry.at(detector - 1);
+
+ TVector3 position(0,0,0);
+ if(theta_n > 0 && phi_n > 0) { // check we have both strips
+ // set position on the wafer
+ position.SetPtThetaPhi(geo.Si_Strip_Theta_Radius.at(theta_n - 1), 0,
+ geo.Si_Strip_Phi_Angle.at(phi_n - 1) );
+ // then translate to the Z-position (from target)
+ position.SetZ(geo.Z);
+ }
+ return position;
+ } catch(std::exception& e) {
+ std::cerr << "RANGE ERROR: GetPositionOfInteraction\n";
+ exit(1);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+TVector3 TAnnularTelescopePhysics::GetRandomisedPositionOfInteraction(int i) const {
+ try {
+ TVector3 pos = GetPositionOfInteraction(i);
+ double R = pos.Pt();
+ double Phi = pos.Phi();
+
+ double R_pitch = m_WedgeGeometry.at(i).Si_Theta_Radius_Pitch;
+ R += gRandom->Uniform(-R_pitch/2., +R_pitch/2.);
+
+ double Phi_pitch = m_WedgeGeometry.at(i).Si_Phi_Angle_Pitch;
+ double Arc_len = R*Phi_pitch;
+ double Arc_pos = R*Phi;
+ Arc_pos += gRandom->Uniform(-Arc_len/2., +Arc_len/2.);
+ Phi = Arc_pos / R;
+
+ TVector3 position(R*cos(Phi), R*sin(Phi), pos.Z());
+ if(pos.Z() == 0) { position.SetXYZ(0,0,0); } // no phi strip
+ return position;
+ } catch(std::exception& e) {
+ cerr << "RANGE ERROR: GetRandomisedPositionOfInteration\n";
+ exit(1);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::BuildSimplePhysicalEvent() {
+ BuildPhysicalEvent();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::BuildPhysicalEvent() {
+ // apply thresholds and calibration
+ PreTreat();
+ // build CsI event
+ BuildCsIEvent();
+ // build Si event
+ BuildSiEvent();
+ // // build total (Si + CsI event)
+ BuildTotalEvent();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::BuildCsIEvent() {
+ // Loop over raw hits
+ size_t mysizeE = m_PreTreatedData->GetCsIMultEnergy();
+ for (size_t i = 0; i < mysizeE ; ++i) {
+ // Read detector, wedge no, energy
+ const Int_t det = m_PreTreatedData->GetCsIE_DetectorNbr(i);
+ const Int_t wedge = m_PreTreatedData->GetCsIE_WedgeNbr(i);
+ const Double_t energy = m_PreTreatedData->GetCsIEnergy(i);
+ //
+ // Look for matching times
+ double time = -1000;
+ size_t mysizeT = m_PreTreatedData->GetCsIMultTime();
+ for(size_t j=0; j< mysizeT; ++j) {
+ const size_t detT = m_PreTreatedData->GetCsIT_DetectorNbr(j);
+ const size_t wedgeT = m_PreTreatedData->GetCsIT_WedgeNbr(j);
+ if(det == detT && wedge == wedgeT) {
+ // Match - set time
+ time = m_PreTreatedData->GetCsITime(j);
+ break;
+ }
+ }
+ // add event
+ SetCsIHit(det,wedge,energy,time);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::BuildSiEvent() {
+ // loop over theta strip energies
+ const size_t mysizeE = m_PreTreatedData->GetSiThetaE_Mult();
+ for (size_t iThetaE = 0; iThetaE < mysizeE; ++iThetaE) {
+ //
+ // Add theta strip energy info
+ const Int_t& det = m_PreTreatedData->GetSiThetaE_DetectorNbr(iThetaE);
+ const Int_t& thetaStrip = m_PreTreatedData->GetSiThetaE_StripNbr(iThetaE);
+ const Double_t& thetaEnergy = m_PreTreatedData->GetSiThetaE_Energy(iThetaE);
+ //
+ // Check for corresponding time
+ Double_t timeTheta = -1000;
+ const size_t mysizeT = m_PreTreatedData->GetSiThetaT_Mult();
+ for (size_t iThetaT = 0; iThetaT < mysizeT; ++iThetaT) {
+ size_t detT = m_PreTreatedData->GetSiThetaT_DetectorNbr(iThetaT);
+ size_t thT = m_PreTreatedData->GetSiThetaT_StripNbr(iThetaT);
+ if(det == detT && thetaStrip == thT) {
+ // Match - set time
+ timeTheta = m_PreTreatedData->GetSiThetaT_Time(iThetaT);
+ break;
+ }
+ }
+ //
+ // Now check for phi match
+ Int_t phiStrip = -1000;
+ Double_t phiEnergy = -1000;
+ Double_t timePhi = -1000;
+ // Check for match in energy window
+ const size_t mysizePhiE = m_PreTreatedData->GetSiPhiE_Mult();
+ for (size_t iPhiE = 0; iPhiE < mysizePhiE; ++iPhiE) {
+ size_t detPhi = m_PreTreatedData->GetSiPhiE_DetectorNbr(iPhiE);
+ Double_t e = m_PreTreatedData->GetSiPhiE_Energy(iPhiE);
+ if(det == detPhi &&
+ fabs(thetaEnergy - e) < m_EnergyMatch.NSigma*m_EnergyMatch.Sigma)
+ {
+ // Match!
+ phiStrip = m_PreTreatedData->GetSiPhiE_StripNbr(iPhiE);
+ phiEnergy = e;
+ //
+ // Now see if we have a time
+ const size_t mysizePhiT = m_PreTreatedData->GetSiPhiT_Mult();
+ for (size_t iPhiT = 0; iPhiT < mysizePhiT; ++iPhiT) {
+ const size_t detPhiT = m_PreTreatedData->GetSiPhiT_DetectorNbr(iPhiT);
+ const size_t phiStripT = m_PreTreatedData->GetSiPhiT_StripNbr(iPhiT);
+ if(det == detPhiT && phiStrip == phiStripT) {
+ // Match - set time
+ timePhi = m_PreTreatedData->GetSiPhiT_Time(iPhiT);
+ break;
+ }
+ }
+ break; // Done - phi energy match found
+ }
+ }
+ // Set hit info
+ SetSiHit(det,thetaStrip,phiStrip,thetaEnergy,phiEnergy,timeTheta,timePhi);
+ }
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+namespace { inline bool match_si_csi(
+ double siPitch, double siAngle, double csiPitch, double csiAngle) {
+ const double siLow = siAngle - siPitch/2.;
+ const double siHigh = siAngle + siPitch/2.;
+ const double csiLow = csiAngle - csiPitch/2.;
+ const double csiHigh = csiAngle + csiPitch/2.;
+
+ if ( (siLow >= csiLow && siLow < csiHigh) ||
+ (siHigh >= csiLow && siHigh < csiLow) ) {
+ return true;
+ }
+ else {
+ return false;
+ }
+}
+
+class SiCsIMatch {
+public:
+ SiCsIMatch(const SiHit_t& siHit, const AnnularTelescope_Utils::Geometry& geo, bool usePitch = true):
+ m_SiHit(siHit), m_Geo(geo), m_UsePitch(usePitch) { }
+ bool operator() (const CsIHit_t& csiHit) {
+ if(csiHit.Detector != m_SiHit.Detector) { return false; }
+ if(csiHit.Wedge < 0) { return false; }
+ if(m_SiHit.PhiStrip < 0) { return false; }
+ if(m_UsePitch) {
+ return match_si_csi(
+ m_Geo.Si_Phi_Angle_Pitch, m_Geo.Si_Strip_Phi_Angle.at(m_SiHit.PhiStrip-1),
+ m_Geo.CsI_Wedge_Angle_Pitch, m_Geo.CsI_Wedge_Phi_Angle.at(csiHit.Wedge-1) );
+ }
+ else { // no si pitch (except for f.p. errors)
+ return match_si_csi(
+ 1e-4, m_Geo.Si_Strip_Phi_Angle.at(m_SiHit.PhiStrip-1),
+ m_Geo.CsI_Wedge_Angle_Pitch, m_Geo.CsI_Wedge_Phi_Angle.at(csiHit.Wedge-1) );
+ }
+ }
+private:
+ const SiHit_t& m_SiHit;
+ const AnnularTelescope_Utils::Geometry& m_Geo;
+ bool m_UsePitch;
+};
+
+}
+void TAnnularTelescopePhysics::BuildTotalEvent(){
+ Int_t iSi = 0;
+ for(const auto& hitSi : m_SiHits) {
+ // Skip events w/o phi strip number
+ if(hitSi.PhiStrip < 0) {
+ continue;
+ }
+ //
+ // Add Si energy, time
+ const Double_t& siEnergy = hitSi.ThetaStripEnergy;
+ const Double_t& siTime = hitSi.ThetaStripTime;
+ TVector3 SiPos = GetPositionOfInteraction(iSi);
+ TVector3 SiPos_R = GetRandomisedPositionOfInteraction(iSi);
+ //
+ // Match CsI wedge and detector (by phi angle)
+ const AnnularTelescope_Utils::Geometry& geo =
+ m_WedgeGeometry.at(hitSi.Detector-1);
+ vector<const CsIHit_t*> csiMatch;
+ for(const auto& csiHit : m_CsIHits){
+ if(csiHit.Wedge < 0) {
+ continue;
+ }
+ if(SiCsIMatch(hitSi, geo)(csiHit)) {
+ csiMatch.push_back(&csiHit);
+ }
+ }
+ // Add CsI info
+ Double_t csiEnergy = -1000;
+ Double_t csiTime = -1000;
+ Int_t csiWedge = -1000;
+ if(!csiMatch.empty()) {
+ if(csiMatch.size() == 1) { // unambiguous
+ csiEnergy = csiMatch.front()->Energy;
+ csiTime = csiMatch.front()->Time;
+ csiWedge = csiMatch.front()->Wedge;
+ } else { // ambiguous, pick the closest one
+ for(const auto& csi: csiMatch) {
+ if(SiCsIMatch(hitSi, geo, false)(*csi)) {
+ csiEnergy = csi->Energy;
+ csiTime = csi->Time;
+ csiWedge = csi->Wedge;
+ }
+ }
+ }
+ }
+ // Set hit info
+ SetHit(siEnergy,siTime,csiEnergy,csiTime,
+ hitSi.Detector,hitSi.ThetaStrip,
+ hitSi.PhiStrip,csiWedge,
+ SiPos,SiPos_R);
+ ++iSi;
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::PreTreat() {
+/********************
+ TODO:: NEED TO CALIBRATE AND APPLY THRESHOLDS
+ (BUT NOT FOR SIMULATION).
+ Here is some example code for calibration
+
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal =
+ CalibrationManager::getInstance();
+
+ Double_t Energy =
+ Cal->ApplyCalibration(
+ "AnnularTelescope/CSI_ENERGY"+NPL::itoa(
+ m_EventData->GetCsIE_DetectorNbr(i) ),
+ m_EventData->GetCsIEnergy(i) );
+*********************/
+
+
+ // clear pre-treated object
+ ClearPreTreatedData();
+
+ // CsI energy
+ for(size_t i=0; i< m_EventData->GetCsIMultEnergy(); ++i){
+ m_PreTreatedData->SetCsIEnergy(
+ m_EventData->GetCsIE_DetectorNbr(i),
+ m_EventData->GetCsIE_WedgeNbr(i),
+ m_EventData->GetCsIEnergy(i) );
+ }
+
+ // CsI time
+ for(size_t i=0; i< m_EventData->GetCsIMultTime(); ++i){
+ m_PreTreatedData->SetCsITime(
+ m_EventData->GetCsIT_DetectorNbr(i),
+ m_EventData->GetCsIT_WedgeNbr(i),
+ m_EventData->GetCsITime(i) );
+ }
+
+ // Si energy (theta strips)
+ for(size_t i=0; i< m_EventData->GetSiThetaE_Mult(); ++i){
+ m_PreTreatedData->SetSiThetaEnergy(
+ m_EventData->GetSiThetaE_DetectorNbr(i),
+ m_EventData->GetSiThetaE_StripNbr(i),
+ m_EventData->GetSiThetaE_Energy(i) );
+ }
+
+ // Si time (theta strips)
+ for(size_t i=0; i< m_EventData->GetSiThetaT_Mult(); ++i){
+ m_PreTreatedData->SetSiThetaTime(
+ m_EventData->GetSiThetaT_DetectorNbr(i),
+ m_EventData->GetSiThetaT_StripNbr(i),
+ m_EventData->GetSiThetaT_Time(i) );
+ }
+
+ // Si energy (phi strips)
+ for(size_t i=0; i< m_EventData->GetSiPhiE_Mult(); ++i){
+ m_PreTreatedData->SetSiPhiEnergy(
+ m_EventData->GetSiPhiE_DetectorNbr(i),
+ m_EventData->GetSiPhiE_StripNbr(i),
+ m_EventData->GetSiPhiE_Energy(i) );
+ }
+
+ // Si time (phi strips)
+ for(size_t i=0; i< m_EventData->GetSiPhiT_Mult(); ++i){
+ m_PreTreatedData->SetSiPhiTime(
+ m_EventData->GetSiPhiT_DetectorNbr(i),
+ m_EventData->GetSiPhiT_StripNbr(i),
+ m_EventData->GetSiPhiT_Time(i) );
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::ReadAnalysisConfig() {
+ // TODO: Implement!
+/*********************
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigAnnularTelescope.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigAnnularTelescope.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigAnnularTelescope.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigAnnularTelescope.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer,DataBuffer,whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ string name = "ConfigAnnularTelescope";
+ if (LineBuffer.compare(0, name.length(), name) == 0)
+ ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus ) {
+ whatToDo="";
+ AnalysisConfigFile >> whatToDo;
+
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+
+ else if (whatToDo=="E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_RAW_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_RAW_Threshold << endl;
+ }
+
+ else if (whatToDo=="E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_Threshold << endl;
+ }
+
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
+**************************/
+}
+
+// ///////////////////////////////////////////////////////////////////////////
+// int TAnnularTelescopePhysics::MatchToSi(const TVector3& SiPosition) const {
+// // return detector number of detector that spans the Si position
+// // for when CsI backs an Si detector
+// // returns -1 if no detector found
+
+// // todo account for offset between detectors
+// int hitno = 0;
+// for(const auto& detno : DetectorNumber){
+// const AnnularTelescope_Utils::Geometry& g = m_WedgeGeometry.at(detno - 1);
+// if(SiPosition.Phi() >= g.Phi_min && SiPosition.Phi() < g.Phi_max &&
+// SiPosition.Pt() >= g.R_min && SiPosition.Pt() < g.R_max &&
+// fabs(SiPosition.Z()/NPUNITS::mm - g.Z/NPUNITS::mm) < 10) {
+// //
+// // Match!
+// return hitno;
+// }
+// ++hitno;
+// }
+// // no match
+// return -1;
+// }
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::Clear() {
+
+ m_CsIHits.clear();
+ m_SiHits.clear();
+ m_Hits.clear();
+
+ CsIHit_Detector.clear();
+ CsIHit_Wedge.clear();
+ CsIHit_Energy.clear();
+ CsIHit_Time.clear();
+
+ SiHit_Detector.clear();
+ SiHit_ThetaStrip.clear();
+ SiHit_PhiStrip.clear();
+ SiHit_ThetaStripEnergy.clear();
+ SiHit_PhiStripEnergy.clear();
+ SiHit_ThetaStripTime.clear();
+ SiHit_PhiStripTime.clear();
+
+ Hit_Si_Energy.clear();
+ Hit_Si_Time.clear();
+ Hit_CsI_Energy.clear();
+ Hit_CsI_Time.clear();
+ Hit_Detector.clear();
+ Hit_Si_ThetaStrip.clear();
+ Hit_Si_PhiStrip.clear();
+ Hit_CsI_Wedge.clear();
+ Hit_Position.clear();
+ Hit_Position_R.clear();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<AnnularTelescope_Utils::Geometry> geometries =
+ AnnularTelescope_Utils::ReadConfiguration(parser);
+ for(const auto& geo : geometries){ AddDetector(geo); }
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::InitSpectra() {
+ m_Spectra = new TAnnularTelescopeSpectra(GetNumberOfDetectors());
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::FillSpectra() {
+ m_Spectra -> FillRawSpectra(m_EventData);
+ m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
+ m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::CheckSpectra() {
+ m_Spectra->CheckSpectra();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::ClearSpectra() {
+ // To be done
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+map< string , TH1*> TAnnularTelescopePhysics::GetSpectra() {
+ if(m_Spectra)
+ return m_Spectra->GetMapHisto();
+ else{
+ map< string , TH1*> empty;
+ return empty;
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::WriteSpectra() {
+ m_Spectra->WriteSpectra();
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ for (int i = 0; i < GetNumberOfDetectors(); ++i) {
+ Cal->AddParameter("AnnularTelescope", "D"+ NPL::itoa(i+1)+"_ENERGY","AnnularTelescope_D"+ NPL::itoa(i+1)+"_ENERGY");
+ Cal->AddParameter("AnnularTelescope", "D"+ NPL::itoa(i+1)+"_TIME","AnnularTelescope_D"+ NPL::itoa(i+1)+"_TIME");
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("AnnularTelescope", true );
+ inputChain->SetBranchAddress("AnnularTelescope", &m_EventData );
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("AnnularTelescope", &m_EventPhysics);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopePhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("AnnularTelescope", "TAnnularTelescopePhysics", &m_EventPhysics);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TAnnularTelescopePhysics::Construct() {
+ return (NPL::VDetector*) new TAnnularTelescopePhysics();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+class proxy_AnnularTelescope{
+ public:
+ proxy_AnnularTelescope(){
+ NPL::DetectorFactory::getInstance()->AddToken("AnnularTelescope","AnnularTelescope");
+ NPL::DetectorFactory::getInstance()->AddDetector("AnnularTelescope",TAnnularTelescopePhysics::Construct);
+ }
+};
+
+proxy_AnnularTelescope p_AnnularTelescope;
+}
+
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.h b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..f312cdbec6a4375cbf8ea657bd19423766913279
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopePhysics.h
@@ -0,0 +1,332 @@
+#ifndef TAnnularTelescopePHYSICS_H
+#define TAnnularTelescopePHYSICS_H
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <vector>
+#include <map>
+#include <string>
+
+// ROOT headers
+#include "TObject.h"
+#include "TH1.h"
+#include "TVector3.h"
+// NPTool headers
+#include "AnnularTelescope_Utils.h"
+#include "TAnnularTelescopeData.h"
+#include "TAnnularTelescopeSpectra.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+// forward declaration
+class TAnnularTelescopeSpectra;
+
+struct CsIHit_t {
+ Int_t Detector;
+ Int_t Wedge;
+ Double_t Energy;
+ Double_t Time;
+};
+struct SiHit_t {
+ Int_t Detector;
+ Int_t ThetaStrip;
+ Int_t PhiStrip;
+ Double_t ThetaStripEnergy;
+ Double_t PhiStripEnergy;
+ Double_t ThetaStripTime;
+ Double_t SiHit_PhiStripTime;
+};
+struct Hit_t {
+ Double_t Si_Energy;
+ Double_t Si_Time;
+ Double_t CsI_Energy;
+ Double_t CsI_Time;
+ Int_t Detector;
+ Int_t Si_ThetaStrip;
+ Int_t Si_PhiStrip;
+ Int_t CsI_Wedge;
+ TVector3 Position;
+ TVector3 Position_R;
+};
+
+
+class TAnnularTelescopePhysics : public TObject, public NPL::VDetector {
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+public:
+ TAnnularTelescopePhysics();
+ ~TAnnularTelescopePhysics() {};
+
+
+ //////////////////////////////////////////////////////////////
+ // Inherited from TObject and overriden to avoid warnings
+public:
+ void Clear();
+ void Clear(const Option_t*) {};
+
+ void SetCsIHit(
+ Int_t detector, Int_t wedge, Double_t energy, Double_t time){
+ CsIHit_t hit = { detector, wedge, energy, time };
+ m_CsIHits.push_back(hit);
+ CsIHit_Detector.push_back(detector);
+ CsIHit_Wedge.push_back(wedge);
+ CsIHit_Energy.push_back(energy);
+ CsIHit_Time.push_back(time);
+ }
+ void SetSiHit(
+ Int_t detector, Int_t thetaStrip, Int_t phiStrip,
+ Double_t thetaStripEnergy, Double_t phiStripEnergy,
+ Double_t thetaStripTime, Double_t phiStripTime) {
+ SiHit_t hit = {
+ detector, thetaStrip, phiStrip,
+ thetaStripEnergy, phiStripEnergy,
+ thetaStripTime, phiStripTime
+ };
+ m_SiHits.push_back(hit);
+ SiHit_Detector. push_back(detector);
+ SiHit_ThetaStrip. push_back(thetaStrip);
+ SiHit_PhiStrip. push_back(phiStrip);
+ SiHit_ThetaStripEnergy.push_back(thetaStripEnergy);
+ SiHit_PhiStripEnergy. push_back(phiStripEnergy);
+ SiHit_ThetaStripTime. push_back(thetaStripTime);
+ SiHit_PhiStripTime. push_back(phiStripTime);
+ }
+ void SetHit(
+ Double_t siEnergy, Double_t siTime, Double_t csiEnergy, Double_t csiTime,
+ Int_t detector, Int_t siThetaStrip, Int_t siPhiStrip, Int_t csiWedge,
+ const TVector3& position, const TVector3& position_R){
+ Hit_t hit;
+ hit.Si_Energy = siEnergy; hit.Si_Time = siTime;
+ hit.CsI_Energy = csiEnergy; hit.CsI_Time = csiTime;
+ hit.Detector = detector; hit.Si_ThetaStrip = siThetaStrip;
+ hit.Si_PhiStrip = siPhiStrip; hit.CsI_Wedge = csiWedge;
+ hit.Position = position; hit.Position_R = position_R;
+
+ m_Hits.push_back(hit);
+ Hit_Si_Energy.push_back(siEnergy);
+ Hit_Si_Time.push_back(siTime);
+ Hit_CsI_Energy.push_back(csiEnergy);
+ Hit_CsI_Time.push_back(csiTime);
+ Hit_Position.push_back(position);
+ Hit_Position_R.push_back(position_R);
+
+ Hit_Detector.push_back(detector);
+ Hit_Si_ThetaStrip.push_back(siThetaStrip);
+ Hit_Si_PhiStrip.push_back(siPhiStrip);
+ Hit_CsI_Wedge.push_back(csiWedge);
+ }
+ const CsIHit_t& GetCsIHit(size_t i) const {
+ try { return m_CsIHits.at(i); }
+ catch(std::exception&) {
+ std::cerr
+ << "ERROR: Invalid index: " << i << " to TAnnularTelescopePhysics::GetCsIHit\n";
+ exit(1);
+ }
+ }
+ const SiHit_t& GetSiHit(size_t i) const {
+ try { return m_SiHits.at(i); }
+ catch(std::exception&) {
+ std::cerr
+ << "ERROR: Invalid index: " << i << " to TAnnularTelescopePhysics::GetSiHit\n";
+ exit(1);
+ }
+ }
+ const Hit_t& GetHit(size_t i) const {
+ try { return m_Hits.at(i); }
+ catch(std::exception&) {
+ std::cerr
+ << "ERROR: Invalid index: " << i << " to TAnnularTelescopePhysics::GetHit\n";
+ exit(1);
+ }
+ }
+ const std::vector<CsIHit_t>& GetCsIHits() const { return m_CsIHits; }
+ const std::vector<SiHit_t>& GetSiHits() const { return m_SiHits; }
+ const std::vector<Hit_t>& GetHits() const { return m_Hits; }
+
+private:
+ std::vector<CsIHit_t> m_CsIHits; //!
+ std::vector<SiHit_t> m_SiHits; //!
+ std::vector<Hit_t> m_Hits; //!
+ //////////////////////////////////////////////////////////////
+ // data obtained after BuildPhysicalEvent() and stored in
+ // output ROOT file
+ //
+ // CsI Hit
+ std::vector<Int_t> CsIHit_Detector;
+ std::vector<Int_t> CsIHit_Wedge;
+ std::vector<Double_t> CsIHit_Energy;
+ std::vector<Double_t> CsIHit_Time;
+ //
+ // Si Hit
+ std::vector<Int_t> SiHit_Detector;
+ std::vector<Int_t> SiHit_ThetaStrip;
+ std::vector<Int_t> SiHit_PhiStrip;
+ std::vector<Double_t> SiHit_ThetaStripEnergy;
+ std::vector<Double_t> SiHit_PhiStripEnergy;
+ std::vector<Double_t> SiHit_ThetaStripTime;
+ std::vector<Double_t> SiHit_PhiStripTime;
+ //
+ // Composed (Si+CsI) hit
+ std::vector<Double_t> Hit_Si_Energy;
+ std::vector<Double_t> Hit_Si_Time;
+ std::vector<Double_t> Hit_CsI_Energy;
+ std::vector<Double_t> Hit_CsI_Time;
+ std::vector<Int_t> Hit_Detector;
+ std::vector<Int_t> Hit_Si_ThetaStrip;
+ std::vector<Int_t> Hit_Si_PhiStrip;
+ std::vector<Int_t> Hit_CsI_Wedge;
+ std::vector<TVector3> Hit_Position;
+ std::vector<TVector3> Hit_Position_R;
+
+public:
+ // Add a detector
+ void AddDetector(const AnnularTelescope_Utils::Geometry& g);
+
+ //////////////////////////////////////////////////////////////
+ // methods inherited from the VDetector ABC class
+ //
+ // read stream from ConfigFile to pick-up detector parameters
+ void ReadConfiguration(NPL::InputParser);
+
+ // add parameters to the CalibrationManger
+ void AddParameterToCalibrationManager();
+
+ // build CsI part of event
+ void BuildCsIEvent();
+
+ // build Si part of event
+ void BuildSiEvent();
+
+ // combine Si + CsI
+ void BuildTotalEvent();
+
+ // method called event by event, aiming at extracting the
+ // physical information from detector
+ void BuildPhysicalEvent();
+
+ // same as BuildPhysicalEvent() method but with a simpler
+ // treatment
+ void BuildSimplePhysicalEvent();
+
+ // same as above but for online analysis
+ void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+
+ // activate raw data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputRaw();
+
+ // activate physics data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputPhysics();
+
+ // create branches of output ROOT file
+ void InitializeRootOutput();
+
+ // clear the raw and physical data objects event by event
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+ // methods related to the TAnnularTelescopeSpectra class
+ // instantiate the TAnnularTelescopeSpectra class and
+ // declare list of histograms
+ void InitSpectra();
+
+ // fill the spectra
+ void FillSpectra();
+
+ // used for Online mainly, sanity check for histograms and
+ // change their color if issues are found, for example
+ void CheckSpectra();
+
+ // used for Online only, clear all the spectra
+ void ClearSpectra();
+
+ // write spectra to ROOT output file
+ void WriteSpectra();
+
+
+ //////////////////////////////////////////////////////////////
+ // specific methods to AnnularTelescope array
+public:
+ // remove bad channels, calibrate the data and apply thresholds
+ void PreTreat();
+
+ // clear the pre-treated object
+ void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+
+ // read the user configuration file. If no file is found, load standard one
+ void ReadAnalysisConfig();
+
+ // give and external TAnnularTelescopeData object to TAnnularTelescopePhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TAnnularTelescopeData* rawDataPointer) {m_EventData = rawDataPointer;}
+
+ // objects are not written in the TTree
+private:
+ TAnnularTelescopeData* m_EventData; //!
+ TAnnularTelescopeData* m_PreTreatedData; //!
+ TAnnularTelescopePhysics* m_EventPhysics; //!
+
+ struct EnergyMatch_t {
+ Int_t NSigma;
+ Double_t Sigma;
+ } m_EnergyMatch; //!
+
+ // getters for raw and pre-treated data object
+public:
+ TAnnularTelescopeData* GetRawData() const {return m_EventData;}
+ TAnnularTelescopeData* GetPreTreatedData() const {return m_PreTreatedData;}
+
+ // getters for interaction position
+ TVector3 GetPositionOfInteraction(int i) const;
+ TVector3 GetRandomisedPositionOfInteraction(int i) const;
+
+ // Number of detectors
+ int GetNumberOfDetectors() const { return m_WedgeGeometry.size(); } //!
+
+ // parameters used in the analysis
+private:
+ // dimensions of wedge
+ std::vector<AnnularTelescope_Utils::Geometry> m_WedgeGeometry; //!
+
+ // spectra class
+private:
+ TAnnularTelescopeSpectra* m_Spectra; //!
+
+public:
+ // Getters
+ std::map<std::string, TH1*> GetSpectra();
+
+ // Static constructor to be passed to the Detector Factory
+public:
+ static NPL::VDetector* Construct();
+
+ ClassDef(TAnnularTelescopePhysics,1) // AnnularTelescopePhysics structure
+};
+#endif
+
+/* Local Variables: */
+/* mode: c++ */
+/* End: */
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.cxx b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b0f520ab98f3723419cd10fa6770a26e8d0ad540
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.cxx
@@ -0,0 +1,100 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Spectra *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// class header
+#include "TAnnularTelescopeSpectra.h"
+
+// STL
+#include <iostream>
+#include <string>
+using namespace std;
+
+// NPTool header
+#include "NPOptionManager.h"
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+TAnnularTelescopeSpectra::TAnnularTelescopeSpectra()
+ : fNumberOfDetectors(0) {
+ SetName("AnnularTelescope");
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+TAnnularTelescopeSpectra::TAnnularTelescopeSpectra(unsigned int NumberOfDetectors) {
+ // if(NPOptionManager::getInstance()->GetVerboseLevel()>0)
+ // cout << "************************************************" << endl
+ // << "TAnnularTelescopeSpectra : Initalizing control spectra for "
+ // << NumberOfDetectors << " Detectors" << endl
+ // << "************************************************" << endl ;
+ SetName("AnnularTelescope");
+ fNumberOfDetectors = NumberOfDetectors;
+
+ InitRawSpectra();
+ InitPreTreatedSpectra();
+ InitPhysicsSpectra();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+TAnnularTelescopeSpectra::~TAnnularTelescopeSpectra() {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::InitRawSpectra() {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::InitPreTreatedSpectra() {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::InitPhysicsSpectra() {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::FillRawSpectra(TAnnularTelescopeData* RawData) {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::FillPreTreatedSpectra(TAnnularTelescopeData* PreTreatedData) {
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TAnnularTelescopeSpectra::FillPhysicsSpectra(TAnnularTelescopePhysics* Physics) {
+}
+
diff --git a/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.h b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.h
new file mode 100644
index 0000000000000000000000000000000000000000..d2e75c0d2582f837a1684684fbd669bca4d8a6e8
--- /dev/null
+++ b/NPLib/Detectors/AnnularTelescope/TAnnularTelescopeSpectra.h
@@ -0,0 +1,62 @@
+#ifndef TAnnularTelescopeSPECTRA_H
+#define TAnnularTelescopeSPECTRA_H
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold AnnularTelescope Spectra *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// NPLib headers
+#include "NPVSpectra.h"
+#include "TAnnularTelescopeData.h"
+#include "TAnnularTelescopePhysics.h"
+
+// Forward Declaration
+class TAnnularTelescopePhysics;
+
+
+class TAnnularTelescopeSpectra : public VSpectra {
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+ public:
+ TAnnularTelescopeSpectra();
+ TAnnularTelescopeSpectra(unsigned int NumberOfDetectors);
+ ~TAnnularTelescopeSpectra();
+
+ //////////////////////////////////////////////////////////////
+ // Initialization methods
+ private:
+ void InitRawSpectra();
+ void InitPreTreatedSpectra();
+ void InitPhysicsSpectra();
+
+ //////////////////////////////////////////////////////////////
+ // Filling methods
+ public:
+ void FillRawSpectra(TAnnularTelescopeData*);
+ void FillPreTreatedSpectra(TAnnularTelescopeData*);
+ void FillPhysicsSpectra(TAnnularTelescopePhysics*);
+
+ //////////////////////////////////////////////////////////////
+ // Detector parameters
+ private:
+ unsigned int fNumberOfDetectors;
+};
+
+#endif
diff --git a/NPLib/Detectors/FPDTamu/TFPDTamuPhysics.cxx b/NPLib/Detectors/FPDTamu/TFPDTamuPhysics.cxx
index 73cee3f748fffddadbd44d873268ae397a73233d..e69b1d056de217bb9aff0053bf7e048ffda46632 100644
--- a/NPLib/Detectors/FPDTamu/TFPDTamuPhysics.cxx
+++ b/NPLib/Detectors/FPDTamu/TFPDTamuPhysics.cxx
@@ -199,7 +199,15 @@ void TFPDTamuPhysics::BuildPhysicalEvent() {
AWireRightCharge.push_back(EnergyR);
// calculate position in X and Z
double wire_length = 2*fabs(AWireLeftPos[det].X())/NPUNITS::cm;
- AWirePositionX.push_back(wire_length*(EnergyL-EnergyR)/(EnergyL+EnergyR));
+ double aWirePositionX_uncalib = wire_length*(EnergyL-EnergyR)/(EnergyL+EnergyR);
+ static string name;
+ name = "FPDTamu/AWire_R";
+ name+= NPL::itoa(det+1) ;
+ name+= "_POS" ;
+ double aWirePositionX_calib = CalibrationManager::getInstance()->ApplyCalibration(name, aWirePositionX_uncalib);
+ AWirePositionX.push_back(aWirePositionX_calib);
+ //cout << name << endl;
+ //cout << det << " " << wire_length << " " << AWirePositionX.at(r) << " " << EnergyL << " " << EnergyR << " " << relpos << " "<< aWirePositionX_uncalib << " " << aWirePositionX_calib << endl;
AWirePositionZ.push_back(AWireLeftPos[det].Z()/NPUNITS::cm);
}
}
@@ -884,6 +892,9 @@ void TFPDTamuPhysics::AddParameterToCalibrationManager() {
Cal->AddParameter("FPDTamu", "AWire_R"+ NPL::itoa(i+1)+"_C"+ NPL::itoa(iCol+1)+"_T",
"AWire_R"+ NPL::itoa(i+1)+"_C"+ NPL::itoa(iCol+1)+"_T");
}
+ Cal->AddParameter("FPDTamu", "AWire_R"+ NPL::itoa(i+1)+"_POS",
+ "AWire_R"+ NPL::itoa(i+1)+"_POS");
+
}
for (int i = 0; i < m_NumberOfPlast; ++i) { // Always 1
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
index 3e06689ca7d9d643340750ab7cec1d243c7678ce..e7c866605eca8a84fdc294fd9d7b529550b244b6 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
@@ -71,7 +71,7 @@ void TGeTAMUPhysics::InitializeStandardParameter(){
m_CryChannelStatus.clear() ;
m_SegChannelStatus.clear() ;
- //enable all channels
+ //enable all channels
vector< bool > CryChannelStatus;
vector< bool > SegChannelStatus;
CryChannelStatus.resize(4,true);
@@ -127,12 +127,12 @@ void TGeTAMUPhysics::ReadAnalysisConfig(){
}
else if (whatToDo== "LOW_GAIN_ENERGY_CRY") {
- m_LowGainCryIsSet = true ;
+ m_LowGainCryIsSet = true ;
cout << whatToDo << " " << m_LowGainCryIsSet << endl; // e.g. DataBuffer = CLOVER03
}
else if (whatToDo== "LOW_GAIN_ENERGY_SEG") {
- m_LowGainSegIsSet = true ;
+ m_LowGainSegIsSet = true ;
cout << whatToDo << " " << m_LowGainSegIsSet << endl; // e.g. DataBuffer = CLOVER03
}
@@ -205,7 +205,7 @@ void TGeTAMUPhysics::ReadAnalysisConfig(){
else if (whatToDo== "ADC_RANDOM_BIN") {
AnalysisConfigFile >> DataBuffer;
- m_ADCRandomBinIsSet = true ;
+ m_ADCRandomBinIsSet = true ;
cout << whatToDo << " " << m_ADCRandomBinIsSet << endl;
}
@@ -376,7 +376,7 @@ for(unsigned int i = 0 ; i < mysizeE ; i++){
}
if(Eraw>=m_Cry_E_Raw_Threshold && IsValidChannel(0, clover, crystal)){
name = "GETAMU/D"+ NPL::itoa(clover)+"_CRY"+ NPL::itoa(crystal);
- if(m_ADCRandomBinIsSet)
+ if(m_ADCRandomBinIsSet)
Eraw += Random->Rndm();
Energy = cal->ApplyCalibration(name+"_E", Eraw);
if(Energy>=m_Cry_E_Threshold){
@@ -395,7 +395,7 @@ for(unsigned int i = 0 ; i < mysizeE ; i++){
clover = m_EventData->GetCoreCloverNbrT(i);
crystal = m_EventData->GetCoreCrystalNbrT(i);
name = "GETAMU/D"+ NPL::itoa(clover)+"_CRY"+ NPL::itoa(crystal);
- if(m_ADCRandomBinIsSet)
+ if(m_ADCRandomBinIsSet)
Traw += Random->Rndm();
Time = cal->ApplyCalibration(name+"_T", Traw);
Singles_CloverMap_CryTN[clover].push_back(crystal);
@@ -423,7 +423,7 @@ for(unsigned int i = 0 ; i < mysizeE ; i++){
}
if(Eraw>=m_Seg_E_Raw_Threshold && IsValidChannel(1, clover, segment)){
name = "GETAMU/D"+ NPL::itoa(clover)+"_SEG"+ NPL::itoa(segment);
- if(m_ADCRandomBinIsSet)
+ if(m_ADCRandomBinIsSet)
Eraw += Random->Rndm();
Energy = cal->ApplyCalibration(name+"_E", Eraw);
if(Energy>=m_Seg_E_Threshold){
@@ -441,7 +441,7 @@ for(unsigned int i = 0 ; i < mysizeE ; i++){
clover = m_EventData->GetSegmentCloverNbrT(i);
segment = m_EventData->GetSegmentSegmentNbrT(i);
name = "GETAMU/D"+ NPL::itoa(clover)+"_SEG"+ NPL::itoa(segment);
- if(m_ADCRandomBinIsSet)
+ if(m_ADCRandomBinIsSet)
Traw += Random->Rndm();
Time = cal->ApplyCalibration(name+"_T", Traw);
Singles_CloverMap_CryTN[clover].push_back(segment);
@@ -480,10 +480,10 @@ void TGeTAMUPhysics::DCSingles( TVector3& BeamBeta){
int cry = Singles_Crystal[iPixel];
int seg = Singles_Segment[iPixel];
double energy = Singles_E[iPixel];
- TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
- // Fill The doppler corrected singles
+ TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
+ // Fill The doppler corrected singles
Singles_DC.push_back(GetDopplerCorrectedEnergy(energy, GammaLabDirection, BeamBeta)); // Doppler Corrected for highest energy
- Singles_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
+ Singles_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
}
}
@@ -495,9 +495,9 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
//////////////////////////////////////////////////////////////////////////////// PART 1
//Pick up the maximum energy from the core of each clover for Doppler Correction.
- // and sum up energies of the same clover
+ // and sum up energies of the same clover
vector<double> tot_E(4, 0); // size 4, all has zero
- vector<double> max_E(4, -1); // maximum stored energy
+ vector<double> max_E(4, -1); // maximum stored energy
vector<unsigned> max_E_pixel(4, 999); // pixel corresponding tot he maximum stored energy
for(unsigned int iPixel = 0 ; iPixel < Singles_E.size() ; iPixel++){
int clv = Singles_Clover[iPixel];
@@ -514,35 +514,35 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
/////////////////////////////////////////////////////////////////////////////// PART 2
// fill the doppler corrected data
for(unsigned int iClover = 0 ; iClover < tot_E.size() ; iClover++){
- //retrieve the total energy
+ //retrieve the total energy
double energy = tot_E[iClover];
if(energy>0){
- //retrieve the maximum pixel
+ //retrieve the maximum pixel
unsigned pixel = max_E_pixel[iClover];
int clv = Singles_Clover[pixel];
int cry = Singles_Crystal[pixel];
- int seg = Singles_Segment[pixel];
+ int seg = Singles_Segment[pixel];
// fill the addback data
AddBack_Clover.push_back(clv);
AddBack_Crystal.push_back(cry);
AddBack_Segment.push_back(seg);
TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
- AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
+ AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
AddBack_X.push_back(GammaLabDirection.X());
AddBack_Y.push_back(GammaLabDirection.Y());
AddBack_Z.push_back(GammaLabDirection.Z());
- AddBack_E.push_back(energy);
+ AddBack_E.push_back(energy);
//Addback energy
AddBack_DC.push_back(GetDopplerCorrectedEnergy(energy, GammaLabDirection, BeamBeta)); // DC(Etot)
}
}
} // end of scheme 1
-
+
else if (m_AddBackMode==2){ // facing clovers (1&3) or (2&4) add-back
- double max24 = -1; // maximum stored energy
+ double max24 = -1; // maximum stored energy
double max13 = -1;
- double totE24 = 0; // total energy
- double totE13 = 0;
+ double totE24 = 0; // total energy
+ double totE13 = 0;
unsigned max_E_pixel24 = -1;
unsigned max_E_pixel13 = -1;
@@ -579,7 +579,7 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
AddBack_Crystal.push_back(cry);
AddBack_Segment.push_back(seg);
TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
- AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
+ AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
AddBack_X.push_back(GammaLabDirection.X());
AddBack_Y.push_back(GammaLabDirection.Y());
AddBack_Z.push_back(GammaLabDirection.Z());
@@ -595,18 +595,18 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
AddBack_Crystal.push_back(cry);
AddBack_Segment.push_back(seg);
TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
- AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
+ AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
AddBack_X.push_back(GammaLabDirection.X());
AddBack_Y.push_back(GammaLabDirection.Y());
AddBack_Z.push_back(GammaLabDirection.Z());
AddBack_E.push_back(totE24);
//Doppler correction
AddBack_DC.push_back(GetDopplerCorrectedEnergy(totE24, GammaLabDirection, BeamBeta)); // Doppler Corrected for highest energy
- }
+ }
} // end of scheme 2
else if (m_AddBackMode==3){ // all clovers add-back
- double maxE = -1; // maximum stored energy
- double totE = 0; // total energy
+ double maxE = -1; // maximum stored energy
+ double totE = 0; // total energy
unsigned max_E_pixel = -1;
//////////////////////////////////////////////////////////////////////////////// PART 1
@@ -632,7 +632,7 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
AddBack_Crystal.push_back(cry);
AddBack_Segment.push_back(seg);
TVector3 GammaLabDirection = GetSegmentPosition(clv,cry,seg);
- AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
+ AddBack_Theta.push_back(GammaLabDirection.Angle(BeamBeta));
AddBack_X.push_back(GammaLabDirection.X());
AddBack_Y.push_back(GammaLabDirection.Y());
AddBack_Z.push_back(GammaLabDirection.Z());
@@ -640,17 +640,22 @@ void TGeTAMUPhysics::AddBack( TVector3& BeamBeta){
//Doppler correction
AddBack_DC.push_back(GetDopplerCorrectedEnergy(totE, GammaLabDirection, BeamBeta)); // Doppler Corrected for highest energy
}
- } // end of scheme 3
+ } // end of scheme 3
}
/////////////////////////////////////////////////
-void TGeTAMUPhysics::AddClover(unsigned int ID ,double R, double Theta, double Phi){
+void TGeTAMUPhysics::AddClover(unsigned int ID ,double R, double Theta, double Phi, double xShift, double yShift, double zShift){
TVector3 Pos(0,0,1);
Pos.SetTheta(Theta);
Pos.SetPhi(Phi);
Pos.SetMag(R);
+ cout << "Original Pos for Clover " << ID << " = (" << Pos.X() << ", " << Pos.Y() << ", " << Pos.Z() << ")" << endl;
+ Pos.SetX(Pos.X()+xShift);
+ Pos.SetY(Pos.Y()+yShift);
+ Pos.SetZ(Pos.Z()+zShift);
+ cout << "New Pos for Clover " << ID << " = (" << Pos.X() << ", " << Pos.Y() << ", " << Pos.Z() << ")" << endl;
m_CloverPosition[ID] = Pos;
return;
}
@@ -661,7 +666,7 @@ TVector3 TGeTAMUPhysics::GetCloverPosition(int& CloverNbr){
/////////////////////////////////////////////////
TVector3 TGeTAMUPhysics::GetCorePosition(int& CloverNbr,int& CoreNbr){
double offset = 20; // mm
- double depth = 40; // mm
+ double depth = 0; // was 40 mm
TVector3 Pos;
TVector3 CloverPos = GetCloverPosition(CloverNbr);
@@ -688,7 +693,7 @@ TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& Se
double offsetX = 20; // 20mm in GeTAMU according to measurments, 33.4 mm in TIGRESS
double offsetY = 20;
- double depth = 40; // 40mm in GeTAMU according to measurments, 45 mm in TIGRESS
+ double depth = 0; // 40mm in GeTAMU according to measurments, 45 mm in TIGRESS
// Changes signs with segment/core combinations
if (CoreNbr==3||CoreNbr==2)
@@ -712,7 +717,7 @@ TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& Se
cout << "Warning: GeTAMU segment number " << SegmentNbr
<< " is out of range\n accepted values: 0 (reserved for core) or 1-3 (Left, Middle, Right segment) " << endl;
- // Define reference axis as the Clover position,
+ // Define reference axis as the Clover position,
// This is a special case to GeTAMU where segmentation is with respect to clover
Pos.RotateUz(CloverPos.Unit());
Pos+=CloverPos;
@@ -728,17 +733,20 @@ void TGeTAMUPhysics::ReadConfiguration(NPL::InputParser parser) {
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " clovers found " << endl;
- vector<string> token = {"CloverID","R","Theta","Phi","Beta"};
+ vector<string> token = {"CloverID","R","Theta","Phi","xShift","yShift","zShift","Beta"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(token)){
double R = blocks[i]->GetDouble("R","mm");
double Theta = blocks[i]->GetDouble("Theta","deg");
double Phi = blocks[i]->GetDouble("Phi","deg");
+ double xShift = blocks[i]->GetDouble("xShift","mm");
+ double yShift = blocks[i]->GetDouble("yShift","mm");
+ double zShift = blocks[i]->GetDouble("zShift","mm");
int id = blocks[i]->GetInt("CloverID");
vector<double> Beta = blocks[i]->GetVectorDouble("Beta","deg");
cout << "WARNING: beta not used, need to be fixed!" << endl;
- AddClover(id,R,Theta,Phi);
+ AddClover(id,R,Theta,Phi,xShift,yShift,zShift);
}
else{
cout << "ERROR: check your input file formatting " << endl;
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
index 4608cbf999ddf0f2ac55f570c4baccb7420de52d..ad83b9708ad3146b135ee10c5add2b22aace09f8 100644
--- a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
@@ -50,7 +50,7 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
void ReadConfiguration(NPL::InputParser);
// Add Parameter to the CalibrationManger
- void AddParameterToCalibrationManager();
+ void AddParameterToCalibrationManager();
// Activated associated Branches and link it to the private member DetectorData address
// In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
@@ -74,7 +74,7 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
// Read the user configuration file; if no file found, load standard one
void ReadAnalysisConfig();
- // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
+ // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
void BuildPhysicalEvent() ;
// Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
@@ -83,9 +83,9 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
void BuildSimplePhysicalEvent(){BuildPhysicalEvent();} ;
// Clear the Event Physics
- void ClearEventPhysics() {Clear();}
- void ClearEventData() {m_EventData->Clear();}
-
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
public:
void PreTreat();
// clear the pre-treated object
@@ -102,16 +102,16 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
map<int, vector <int> > Singles_CloverMap_CryEN; //! cry number energy
map<int, vector <int> > Singles_CloverMap_SegEN; //1 seg number
map<int, vector <double> > Singles_CloverMap_CryE; //! cry energy
- map<int, vector <double> > Singles_CloverMap_SegE; //! seg energy
+ map<int, vector <double> > Singles_CloverMap_SegE; //! seg energy
map<int, vector <int> > Singles_CloverMap_CryTN; //! cry number time
map<int, vector <int> > Singles_CloverMap_SegTN; //! seg number
map<int, vector <double> > Singles_CloverMap_CryT; //! cry energy
- map<int, vector <double> > Singles_CloverMap_SegT; //! seg energy
-
+ map<int, vector <double> > Singles_CloverMap_SegT; //! seg energy
+
public: // Data Member
//sorting parameters
- vector<double> Singles_E;
- vector<double> Singles_T;
+ vector<double> Singles_E;
+ vector<double> Singles_T;
vector<double> Singles_DC; // Doppler Corrected Energy (filled externaly)
vector<double> Singles_Theta;
vector<double> Singles_X;
@@ -121,8 +121,8 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
vector<int> Singles_Crystal;
vector<int> Singles_Segment;
// add back parameters
- vector<double> AddBack_E;
- vector<double> AddBack_T;
+ vector<double> AddBack_E;
+ vector<double> AddBack_T;
vector<double> AddBack_DC; // Doppler Corrected Energy
vector<double> AddBack_Theta;
vector<double> AddBack_X;
@@ -140,11 +140,11 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
double m_Seg_E_Threshold;
int m_Cry_E_Raw_Threshold;
int m_Seg_E_Raw_Threshold;
- int m_AddBackMode;
- bool m_LowGainCryIsSet;
- bool m_LowGainSegIsSet;
+ int m_AddBackMode;
+ bool m_LowGainCryIsSet;
+ bool m_LowGainSegIsSet;
bool m_ADCRandomBinIsSet; //Randomise the raw energy in the Raw data within a bin
-
+
private: // use for anlysis
TLorentzVector m_GammaLV; //!
@@ -152,7 +152,7 @@ class TGeTAMUPhysics : public TObject, public NPL::VDetector{
TVector3 GetPositionOfInteraction(unsigned int& i);
double GetDopplerCorrectedEnergy(double& energy , TVector3 position, TVector3& beta);
// Add a detector and computes its coordinate
- void AddClover(unsigned int ID, double R, double Theta, double Phi);
+ void AddClover(unsigned int ID, double R, double Theta, double Phi, double xShift, double yShift, double zShift);
TVector3 GetCloverPosition(int& CloverNbr);
TVector3 GetCorePosition(int& CloverNbr, int& CoreNbr);
TVector3 GetSegmentPosition(int& CloverNbr, int& CoreNbr, int& SegmentNbr);
diff --git a/NPLib/Detectors/MDM/CMakeLists.txt b/NPLib/Detectors/MDM/CMakeLists.txt
index 5bfc2205bdfbcec0fea4f25dfe5af053816eb968..0d433a44bcb999665b7ab0e40efb96b7aed4199d 100644
--- a/NPLib/Detectors/MDM/CMakeLists.txt
+++ b/NPLib/Detectors/MDM/CMakeLists.txt
@@ -4,6 +4,16 @@ add_custom_command(OUTPUT TMDMPhysicsDict.cxx COMMAND ../../scripts/build_dict.s
add_custom_command(OUTPUT TMDMPhysicsMinimizerDict.cxx COMMAND ../../scripts/build_dict.sh TMDMPhysicsMinimizer.h TMDMPhysicsMinimizerDict.cxx TMDMPhysicsMinimizer.rootmap libNPMDM.dylib DEPENDS TMDMPhysicsMinimizer.h)
add_custom_command(OUTPUT TMDMDataDict.cxx COMMAND ../../scripts/build_dict.sh TMDMData.h TMDMDataDict.cxx TMDMData.rootmap libNPMDM.dylib DEPENDS TMDMData.h)
+##
+## Check if MINUIT2 is installed along with ROOT
+find_library(libMinuit2_FOUND NAMES Minuit2 HINTS "${ROOTSYS}/lib")
+if(libMinuit2_FOUND)
+ message(STATUS "Minuit2 support enabled for MDM.")
+ add_definitions(-DHAVE_MINUIT2)
+else()
+ message(STATUS "Minuit2 support disabled for MDM.")
+endif()
+
## Check for FORTRAN compiler, if so
## compile RAYTRACE code
if("${CMAKE_GENERATOR}" STREQUAL "Unix Makefiles")
@@ -12,7 +22,7 @@ if("${CMAKE_GENERATOR}" STREQUAL "Unix Makefiles")
message(STATUS "Compiling libMDM with RAYTRACE support included.")
add_definitions(-DUSE_RAYTRACE)
enable_language(Fortran)
- set (CMAKE_Fortran_FLAGS "-O3 -finit-local-zero -falign-commons -fno-automatic")
+ set (CMAKE_Fortran_FLAGS "-O3 -finit-local-zero -falign-commons -fno-automatic -w")
add_library(NPMDM SHARED TMDMSpectra.cxx TMDMData.cxx TMDMPhysics.cxx TMDMPhysicsMinimizer TMDMDataDict.cxx TMDMPhysicsDict.cxx TMDMPhysicsMinimizerDict.cxx MDMTrace.cpp RAYTKIN1.F)
else()
## No fortran support, compile "fake" c-version
@@ -24,6 +34,13 @@ else()
add_library(NPMDM SHARED TMDMSpectra.cxx TMDMData.cxx TMDMPhysics.cxx TMDMPhysicsMinimizer TMDMDataDict.cxx TMDMPhysicsDict.cxx TMDMPhysicsMinimizerDict.cxx MDMTrace.cpp)
endif()
-target_link_libraries(NPMDM ${ROOT_LIBRARIES} NPCore NPPhysics Minuit2)
+if(Minuit2_FOUND)
+## Link to Minuit2 library
+ target_link_libraries(NPMDM ${ROOT_LIBRARIES} NPCore NPPhysics Minuit2)
+else()
+## No Minuit2 library, skip linking
+ target_link_libraries(NPMDM ${ROOT_LIBRARIES} NPCore NPPhysics)
+endif()
+
install(FILES TMDMData.h TMDMPhysics.h TMDMPhysicsMinimizer.h TMDMSpectra.h MDMTrace.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/MDM/TMDMPhysics.cxx b/NPLib/Detectors/MDM/TMDMPhysics.cxx
index 1e21b623c457f3835e8dda9897e9d26e18cb92dd..b495ef65087aad6c9820a89b4be8206ec2835504 100644
--- a/NPLib/Detectors/MDM/TMDMPhysics.cxx
+++ b/NPLib/Detectors/MDM/TMDMPhysics.cxx
@@ -46,10 +46,11 @@ using namespace std;
#include "TGraph.h"
#include "TROOT.h"
#include "TVector3.h"
+#ifdef HAVE_MINUIT2
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
-
+#endif
ClassImp(TMDMPhysics)
@@ -285,7 +286,7 @@ void TMDMPhysics::ReadAnalysisConfig() {
DataBuffer == "false" ? false : atoi(DataBuffer.c_str());
cout << "\t" << whatToDo << " " << m_DoMinimization << endl;
}
-
+
else if (whatToDo=="MINIMIZER_NAME") {
AnalysisConfigFile >> DataBuffer;
m_MinimizerName = DataBuffer;
@@ -376,6 +377,9 @@ void TMDMPhysics::ReadConfiguration(NPL::InputParser parser) {
// Read analysis config file & initialize relavant variables
ReadAnalysisConfig();
+#ifndef HAVE_MINUIT2
+ m_DoMinimization = false;
+#endif
m_Particle = new NPL::Nucleus(m_ParticleZ, m_ParticleA);
@@ -503,6 +507,14 @@ void TMDMPhysics::SendRay(double thetaX,double thetaY,double Ekin) const {
}
+double TMDMPhysics::CalculateCentralEnergy(){
+ double brho = (m_Field/tesla)*1.6; // tesla*meter
+ m_Particle->SetBrho(brho);
+ m_Particle->BrhoToEnergy(m_ParticleQ); // charge state
+ return m_Particle->GetEnergy()/MeV;
+}
+
+
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
@@ -528,6 +540,16 @@ proxy_MDM p_MDM;
}
+
+/////////////////////////////////////////////////////
+//// The following routines depend on MINUIT2 //////
+
+
+#ifdef HAVE_MINUIT2
+#pragma message "Compiling TMDMPhysics with Minuit2 support"
+
+// Define real routines using minuit2
+//
void TMDMPhysics::MinimizeTarget(){
//
// check if we do the minimization
@@ -535,7 +557,7 @@ void TMDMPhysics::MinimizeTarget(){
return;
}
// Set up minimizer
- std::unique_ptr<ROOT::Math::Minimizer> min(
+ std::unique_ptr<Minimizer_t> min(
ROOT::Math::Factory::CreateMinimizer(
m_MinimizerName, m_MinimizerAlgorithm.c_str()
)
@@ -585,16 +607,24 @@ void TMDMPhysics::MinimizeTarget(){
Fit_Chi2 = m_MinimizerFunction->operator()(min->X());
}
-
-void TMDMPhysics::InitializeMinimizerWithDefaults(ROOT::Math::Minimizer& min){
+void TMDMPhysics::InitializeMinimizerWithDefaults(Minimizer_t& min){
min.SetMaxFunctionCalls(1000);
min.SetMaxIterations(1000);
min.SetTolerance(0.001);
}
-double TMDMPhysics::CalculateCentralEnergy(){
- double brho = (m_Field/tesla)*1.6; // tesla*meter
- m_Particle->SetBrho(brho);
- m_Particle->BrhoToEnergy(m_ParticleQ); // charge state
- return m_Particle->GetEnergy()/MeV;
+#else
+
+// Define fake routines without minuit2
+//
+#pragma message "Compiling TMDMPhysics without Minuit2 support"
+
+void TMDMPhysics::MinimizeTarget(){
+ // Do Nothing
}
+
+void TMDMPhysics::InitializeMinimizerWithDefaults(Minimizer_t& min){
+ // Do Nothing
+}
+
+#endif
diff --git a/NPLib/Detectors/MDM/TMDMPhysics.h b/NPLib/Detectors/MDM/TMDMPhysics.h
index 6c8d92545123a1a19a605695742c7cb5f5dd7e5a..6dd6b537bf50ade572b79cd66b3045a392930349 100644
--- a/NPLib/Detectors/MDM/TMDMPhysics.h
+++ b/NPLib/Detectors/MDM/TMDMPhysics.h
@@ -41,11 +41,24 @@
#include "NPReaction.h"
// MDM Trace
#include "MDMTrace.h"
-#include "TMDMPhysicsMinimizer.h"
// forward declaration
class TMDMSpectra;
+
+#ifdef HAVE_MINUIT2
+#include "TMDMPhysicsMinimizer.h"
+
namespace ROOT { namespace Math { class Minimizer; } }
+typedef ROOT::Math::Minimizer Minimizer_t;
+
+#else
+
+class TMDMPhysics;
+struct TMDMPhysicsMinimizer { };
+struct Minimizer_t { };
+
+#endif
+
class TMDMPhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
@@ -201,7 +214,7 @@ public:
// Minimization options & helpers
double CalculateCentralEnergy();
- void InitializeMinimizerWithDefaults(ROOT::Math::Minimizer& min);
+ void InitializeMinimizerWithDefaults(Minimizer_t& min);
private:
TMDMPhysicsMinimizer* m_MinimizerFunction; //!
diff --git a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
index 395185a13d648364ede296de53378296e41c7e3c..72c816b8dcc2580c9eae8ff1f7176f86b94c7ec4 100644
--- a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
+++ b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
@@ -77,18 +77,18 @@ void TTiaraBarrelPhysics::BuildPhysicalEvent(){
unsigned int sizeU = m_PreTreatedData->GetFrontUpstreamEMult();
unsigned int sizeO = m_PreTreatedData->GetOuterEMult();
- static string name; // token
+ static string name; // token
- for(unsigned int i = 0 ; i < sizeU ; i++){
- // detector and strip
+ for(unsigned int i = 0 ; i < sizeU ; i++){
+ // detector and strip
int det = m_PreTreatedData->GetFrontUpstreamEDetectorNbr(i);
int strip = m_PreTreatedData->GetFrontUpstreamEStripNbr(i) ;
- // calibrated energy
+ // calibrated energy
double EU = m_PreTreatedData->GetFrontUpstreamEEnergy(i) ;
- double ED = m_PreTreatedData->GetFrontDownstreamEEnergy(i);
+ double ED = m_PreTreatedData->GetFrontDownstreamEEnergy(i);
// matchsticked energy for position calibration
- double msU = m_PreTreatedMSData->GetFrontUpstreamEEnergy(i) ;
- double msD = m_PreTreatedMSData->GetFrontDownstreamEEnergy(i);
+ double msU = m_PreTreatedMSData->GetFrontUpstreamEEnergy(i) ;
+ double msD = m_PreTreatedMSData->GetFrontDownstreamEEnergy(i);
double RowPos = (msU-msD)/(msU+msD);
name = "TIARABARREL/B";
@@ -99,7 +99,7 @@ void TTiaraBarrelPhysics::BuildPhysicalEvent(){
double Pos = CalibrationManager::getInstance()->ApplyResistivePositionCalibration(name,RowPos); // returns ((RowPos-d)/k)
//Fix Balistic deficit
- // calibration is applied as: (U+D)*( 1 + BD*(pow(k,2)-pow(pos-d,2)) ),
+ // calibration is applied as: (U+D)*( 1 + BD*(pow(k,2)-pow(pos-d,2)) ),
// While BD > 0 and |k| >= |pos-d| for good events
// Get resistive shift and length, this will fix asymetries
double d = CalibrationManager::getInstance()->GetValue(name,0); // resistive strip length shift
@@ -109,7 +109,7 @@ void TTiaraBarrelPhysics::BuildPhysicalEvent(){
name+="_STRIP";
name+=NPL::itoa(strip);
name+="_BALLISTIC";
- double BD_x_k2 =CalibrationManager::getInstance()->ApplyCalibration(name, k );
+ double BD_x_k2 =CalibrationManager::getInstance()->ApplyCalibration(name, k );
double BD_x_Pos2 =CalibrationManager::getInstance()->ApplyCalibration(name, (RowPos-d) );
double BD = (BD_x_k2 - BD_x_Pos2);
@@ -118,20 +118,20 @@ void TTiaraBarrelPhysics::BuildPhysicalEvent(){
Strip_N.push_back(strip);
Strip_Pos.push_back(Pos); // position expressed in [-1;+1]
UpStream_E.push_back(EU);
- DownStream_E.push_back(ED);
+ DownStream_E.push_back(ED);
Strip_E.push_back( (EU+ED) * (1+BD) );
- //cout << det << " " << strip << " " << Pos << " " << EU << " " << ED << " " << endl ;
- }
+ //cout << det << " " << strip << " " << Pos << " " << EU << " " << ED << " " << endl ;
+ }
//Outer Barrel
- for(unsigned int i = 0 ; i < sizeO ; i++){
- // detector and strip
+ for(unsigned int i = 0 ; i < sizeO ; i++){
+ // detector and strip
int det = m_PreTreatedData->GetOuterEDetectorNbr(i);
int strip = m_PreTreatedData->GetOuterEStripNbr(i) ;
double EO = m_PreTreatedData->GetOuterEEnergy(i);
Outer_Detector_N.push_back(det);
Outer_Strip_N.push_back(strip);
- Outer_Strip_E.push_back(EO);
+ Outer_Strip_E.push_back(EO);
}
}
@@ -145,18 +145,18 @@ void TTiaraBarrelPhysics::PreTreat(){
unsigned int sizeB = m_EventData->GetBackEMult();
unsigned int sizeO = m_EventData->GetOuterEMult();
- for(unsigned int i = 0 ; i < sizeU ; i++){
+ for(unsigned int i = 0 ; i < sizeU ; i++){
double EU = Cal_Strip_Upstream_E(i) ;
int det = m_EventData->GetFrontUpstreamEDetectorNbr(i);
int strip = m_EventData->GetFrontUpstreamEStripNbr(i);
- int key = det*10+strip; // key of the map
+ int key = det*10+strip; // key of the map
if(EU > 0) { // threshold on energy is applied below
m_mapU[key].push_back(EU);
m_mapMSU[key].push_back(Match_Strip_Upstream_E(i));
}
}
- for(unsigned int i = 0 ; i < sizeD ; i++){
+ for(unsigned int i = 0 ; i < sizeD ; i++){
double ED = Cal_Strip_Downstream_E(i) ;
int det = m_EventData->GetFrontDownstreamEDetectorNbr(i);
int strip = m_EventData->GetFrontDownstreamEStripNbr(i);
@@ -167,24 +167,24 @@ void TTiaraBarrelPhysics::PreTreat(){
}
}
- for(unsigned int i = 0 ; i < sizeB ; i++){
+ for(unsigned int i = 0 ; i < sizeB ; i++){
double EB = Cal_Back_E(i) ;
int det = m_EventData->GetBackEDetectorNbr(i);
- int key = det; // key of the map
+ int key = det; // key of the map
if(EB > m_Back_E_Threshold) m_mapB[key].push_back(EB);
}
- for(unsigned int i = 0 ; i < sizeO ; i++){
+ for(unsigned int i = 0 ; i < sizeO ; i++){
double EO = m_EventData->GetOuterEEnergy(i);
int det = m_EventData->GetOuterEDetectorNbr(i);
int strip = m_EventData->GetOuterEStripNbr(i);
- int key = det*10+strip; // key of the map OuterStrip={1,2,3,4} => key
+ int key = det*10+strip; // key of the map OuterStrip={1,2,3,4} => key
if(EO > m_OuterBack_E_Threshold) m_mapO[key].push_back(EO);
}
-//Apply selection and matching
+//Apply selection and matching
//NOTE about Barrel Matching
- // Applying a strong strip-matching condition between the inner and outer barrel might not be adequate
+ // Applying a strong strip-matching condition between the inner and outer barrel might not be adequate
// in the case of a large beam spot, since some particles at specific angles can fire an Inner-strip (n)
// and an Outer-strip (n+/-1). The strip-matching can be addressed in the user Analysis.cxx
@@ -192,8 +192,8 @@ double EU, ED, EUms, EDms, EB, EO ;
map<int,vector <double> >::iterator it;
for (it= m_mapU.begin(); it!=m_mapU.end(); ++it){
- // Define the detector and strip
- int key = it->first ;
+ // Define the detector and strip
+ int key = it->first ;
int strip = (key)%10;
int det = (key)/10;
EU=ED=EUms=EDms=0;
@@ -206,7 +206,7 @@ for (it= m_mapU.begin(); it!=m_mapU.end(); ++it){
EDms = m_mapMSD[key][0];
m_PreTreatedData->SetFrontUpstreamE(det,strip,EU);
m_PreTreatedData->SetFrontDownstreamE(det,strip,ED);
- m_PreTreatedMSData->SetFrontUpstreamE(det,strip,EUms);
+ m_PreTreatedMSData->SetFrontUpstreamE(det,strip,EUms);
m_PreTreatedMSData->SetFrontDownstreamE(det,strip,EDms);
}
}
@@ -222,8 +222,8 @@ for (it= m_mapB.begin(); it!=m_mapB.end(); ++it){
}
for (it= m_mapO.begin(); it!=m_mapO.end(); ++it){
- EO=0;
- int key = it->first ;
+ EO=0;
+ int key = it->first ;
int strip = (key)%10;
int det = (key)/10;
if (m_mapO[key].size()==1){
@@ -295,7 +295,7 @@ void TTiaraBarrelPhysics::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("Tiara","Barrel");
if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
+ cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> token = {"InnerBarrel","OuterBarrel","Chamber","X","Y","Z"};
@@ -310,6 +310,9 @@ void TTiaraBarrelPhysics::ReadConfiguration(NPL::InputParser parser){
double Y = blocks[i]->GetDouble("Y","mm");
double Z = blocks[i]->GetDouble("Z","mm");
AddDetector(X,Y,Z);
+ BarrelXCoord = X;
+ BarrelYCoord = Y;
+ BarrelZCoord = Z;
}
else{
@@ -322,22 +325,22 @@ void TTiaraBarrelPhysics::ReadConfiguration(NPL::InputParser parser){
ReadAnalysisConfig();
}
///////////////////////////////////////////////////////////////////////////
-void TTiaraBarrelPhysics::InitSpectra(){
+void TTiaraBarrelPhysics::InitSpectra(){
m_Spectra = new TTiaraBarrelSpectra(m_NumberOfDetector);
}
///////////////////////////////////////////////////////////////////////////
-void TTiaraBarrelPhysics::FillSpectra(){
+void TTiaraBarrelPhysics::FillSpectra(){
m_Spectra -> FillRawSpectra(m_EventData);
m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
-void TTiaraBarrelPhysics::CheckSpectra(){
+void TTiaraBarrelPhysics::CheckSpectra(){
// To be done
}
///////////////////////////////////////////////////////////////////////////
-void TTiaraBarrelPhysics::ClearSpectra(){
+void TTiaraBarrelPhysics::ClearSpectra(){
// To be done
}
///////////////////////////////////////////////////////////////////////////
@@ -349,7 +352,7 @@ map< string,TH1* > TTiaraBarrelPhysics::GetSpectra() {
return empty;
}
-}
+}
///////////////////////////////////////////////////////////////////////////
void TTiaraBarrelPhysics::AddParameterToCalibrationManager(){
@@ -419,9 +422,9 @@ void TTiaraBarrelPhysics::AddDetector(double X,double Y,double Z){
m_NumberOfDetector+=8;
- /*
+ /*
double StripPitchSector = (Wedge_Phi_Max-Wedge_Phi_Min)/Wedge_Sector_NumberOfStrip ;
- double StripPitchRing = (Wedge_R_Max-Wedge_R_Min)/Wedge_Ring_NumberOfStrip ;
+ double StripPitchRing = (Wedge_R_Max-Wedge_R_Min)/Wedge_Ring_NumberOfStrip ;
TVector3 Strip_1_1;
@@ -468,30 +471,30 @@ TVector3 TTiaraBarrelPhysics::GetDetectorNormal( const int i) const{
}
///////////////////////////////////////////////////////////////////////////////
TVector3 TTiaraBarrelPhysics::GetPositionOfInteraction(const int i) const{
- // All in mm
- double INNERBARREL_PCB_Width = 27.76;
- double INNERBARREL_ActiveWafer_Length = 94.80;
- double INNERBARREL_ActiveWafer_Width = 24.0;
+ // All in mm
+ double INNERBARREL_PCB_Width = 27.10; // was 27.76
+ double INNERBARREL_ActiveWafer_Length = 94.80;
+ double INNERBARREL_ActiveWafer_Width = 22.60; // was 24.60
double StripPitch = INNERBARREL_ActiveWafer_Width/4.0;
//Calculate position locally as if it's detector 3 (at 12'oclock) that is hit
double X = (Strip_N[i]-0.5)*StripPitch - 0.5*INNERBARREL_ActiveWafer_Width;
double Y = INNERBARREL_PCB_Width*(0.5+sin(45*deg));
- double Z = Strip_Pos[i]*(0.5*INNERBARREL_ActiveWafer_Length);
+ double Z = Strip_Pos[i]*(0.5*INNERBARREL_ActiveWafer_Length);
TVector3 POS(X,Y,-Z); // since RowPos = (U-D)/(U+D) => Downstream hit (i.e. Z>0) has RowPos<0, thus the sign
- POS.RotateZ((3-Detector_N[i])*45*deg);// looking downstream, Detector 1 is at 3 o'clock
+ POS.RotateZ((3-Detector_N[i])*45*deg);// looking downstream, Detector 1 is at 3 o'clock
return( POS ) ;
}
///////////////////////////////////////////////////////////////////////////////
TVector3 TTiaraBarrelPhysics::GetRandomisedPositionOfInteraction(const int i) const{
TVector3 RandomPOS = GetPositionOfInteraction(i);
- TVector3 v1(-12.0, 27.76*(0.5+sin(45*deg)), 0.0);
- TVector3 v2(12.0, 27.76*(0.5+sin(45*deg)), 0.0);
+ TVector3 v1(-12.0, 27.10*(0.5+sin(45*deg)), 0.0);
+ TVector3 v2(12.0, 27.10*(0.5+sin(45*deg)), 0.0);
v1.RotateZ((3-Detector_N[i])*45*deg);
v2.RotateZ((3-Detector_N[i])*45*deg);
TVector3 u = (v2-v1).Unit();
double RandomNumber = Random->Rndm();
- TVector3 DeltaHolder((RandomNumber*6.0)-3.0,(RandomNumber*6.0)-3.0,0.0);
+ TVector3 DeltaHolder((RandomNumber*5.65)-(5.65/2),(RandomNumber*5.65)-(5.65/2),0.0);
TVector3 DeltaVector(u.X()*DeltaHolder.X(),u.Y()*DeltaHolder.Y(),u.Z()*DeltaHolder.Z());
RandomPOS.SetXYZ(RandomPOS.X()+DeltaVector.X(),RandomPOS.Y()+DeltaVector.Y(),RandomPOS.Z()+DeltaVector.Z());
return( RandomPOS );
@@ -556,7 +559,7 @@ double TTiaraBarrelPhysics::Match_Strip_Downstream_E(const int i){
name+= NPL::itoa( m_EventData->GetFrontDownstreamEDetectorNbr(i) ) ;
name+= "_DOWNSTREAM" ;
name+= NPL::itoa( m_EventData->GetFrontDownstreamEStripNbr(i) ) ;
- name+= "_MATCHSTICK";
+ name+= "_MATCHSTICK";
double RawEnergy = m_EventData->GetFrontDownstreamEEnergy(i);
double MSEnergy = CalibrationManager::getInstance()->ApplyCalibration(name,RawEnergy);
return MSEnergy;
@@ -596,4 +599,3 @@ class brlproxy{
brlproxy brlp;
}
-
diff --git a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.h b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.h
index bd3b7cd1069a66ca1f4c09ee9cf663b8f8370e09..ba66908d30b94480fdf28dc9a6e08f47a1595281 100644
--- a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.h
+++ b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.h
@@ -18,7 +18,7 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
- * *
+ * *
* *
* *
*****************************************************************************/
@@ -31,9 +31,9 @@
#include "NPCalibrationManager.h"
#include "NPVDetector.h"
#include "NPInputParser.h"
-// ROOT
-#include "TVector2.h"
-#include "TVector3.h"
+// ROOT
+#include "TVector2.h"
+#include "TVector3.h"
#include "TObject.h"
#include "TH1.h"
@@ -46,8 +46,8 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
TTiaraBarrelPhysics();
~TTiaraBarrelPhysics() {};
- public:
- void Clear();
+ public:
+ void Clear();
void Clear(const Option_t*) {};
public:
@@ -63,8 +63,8 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
vector<double> Strip_T;
vector<int> Strip_N;
vector<double> Strip_Pos;
-
- // Control stuff
+
+ // Control stuff
vector<double> DownStream_E;
vector<double> DownStream_T;
vector<double> UpStream_E;
@@ -79,12 +79,17 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
vector<double> Outer_Back_E;
vector<double> Outer_Back_T;
+ // Coordinates
+ double BarrelXCoord;
+ double BarrelYCoord;
+ double BarrelZCoord;
+
public: // Innherited from VDetector Class
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
void ReadConfiguration(NPL::InputParser) ;
// Add Parameter to the CalibrationManger
- void AddParameterToCalibrationManager() ;
+ void AddParameterToCalibrationManager() ;
// Activated associated Branches and link it to the private member DetectorData address
// In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
@@ -97,7 +102,7 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
// Create associated branches and associated private member DetectorPhysics address
void InitializeRootOutput() ;
- // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
+ // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
void BuildPhysicalEvent() ;
// Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
@@ -109,16 +114,16 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
// Those two method all to clear the Event Physics or Data
- void ClearEventPhysics() {Clear();}
- void ClearEventData() {m_EventData->Clear();}
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
- // Method related to the TSpectra classes, aimed at providing a framework
+ // Method related to the TSpectra classes, aimed at providing a framework
// for online applications
// Instantiate the Spectra class and the histogramm throught it
void InitSpectra();
// Fill the spectra hold by the spectra class
void FillSpectra();
- // Used for Online mainly, perform check on the histo and for example change
+ // Used for Online mainly, perform check on the histo and for example change
// their color if issues are found
void CheckSpectra();
// Used for Online only, clear all the spectra hold by the Spectra class
@@ -145,7 +150,7 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
// Add a Detector
void AddDetector( double X,double Y,double Z);
-
+
// Give and external TMustData object to TTiaraBarrelPhysics. Needed for online analysis for example.
void SetRawDataPointer(TTiaraBarrelData* rawDataPointer) {m_EventData = rawDataPointer;}
// Retrieve raw and pre-treated data
@@ -154,11 +159,15 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
double GetNumberOfDetector() const { return m_NumberOfDetector; };
- // To be called after a build Physical Event
+ // To be called after a build Physical Event
int GetEventMultiplicity() const { return EventMultiplicity; };
- TVector3 GetPositionOfInteraction(const int i) const;
- TVector3 GetRandomisedPositionOfInteraction(const int i) const;
+ double GetBarrelXCoord() const { return BarrelXCoord; };
+ double GetBarrelYCoord() const { return BarrelYCoord; };
+ double GetBarrelZCoord() const { return BarrelZCoord; };
+
+ TVector3 GetPositionOfInteraction(const int i) const;
+ TVector3 GetRandomisedPositionOfInteraction(const int i) const;
TVector3 GetDetectorNormal(const int i) const;
private: // Parameter used in the analysis
@@ -177,12 +186,12 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
TTiaraBarrelData* m_PreTreatedMSData;//! stores the intermediate Matchsticks calibrated Data
TTiaraBarrelPhysics* m_EventPhysics;//!
map<int, vector <double> > m_mapU;//! the maps sorts out the data before storing in m_PreTreatedData
- map<int, vector <double> > m_mapD;//!
- map<int, vector <double> > m_mapB;//!
- map<int, vector <double> > m_mapO;//!
- map<int, vector <double> > m_mapMSU;//!
- map<int, vector <double> > m_mapMSD;//!
-
+ map<int, vector <double> > m_mapD;//!
+ map<int, vector <double> > m_mapB;//!
+ map<int, vector <double> > m_mapO;//!
+ map<int, vector <double> > m_mapMSU;//!
+ map<int, vector <double> > m_mapMSD;//!
+
private: // Map of activated channel
map< int, vector<bool> > m_InnerBarrelStripUpstreamChannelStatus;//!
map< int, vector<bool> > m_InnerBarrelStripDownstreamChannelStatus;//!
@@ -200,7 +209,7 @@ class TTiaraBarrelPhysics : public TObject, public NPL::VDetector{
TTiaraBarrelSpectra* m_Spectra;//!
public:
- map< string,TH1* > GetSpectra();
+ map< string,TH1* > GetSpectra();
private: // Usefull method
// Calibrate data
diff --git a/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.cc b/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ca06a0f219e977fb54415a4fa6c867859cc34bb8
--- /dev/null
+++ b/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.cc
@@ -0,0 +1,454 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe AnnularTelescope simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <sstream>
+#include <cmath>
+#include <limits>
+//G4 Geometry object
+#include "G4Tubs.hh"
+#include "G4Box.hh"
+
+//G4 sensitive
+#include "G4SDManager.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+//G4 various object
+#include "G4Material.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+
+// NPTool header
+#include "AnnularTelescope.hh"
+#include "CalorimeterScorers.hh"
+#include "RootOutput.h"
+#include "MaterialManager.hh"
+#include "NPSDetectorFactory.hh"
+#include "NPOptionManager.h"
+#include "NPSHitsMap.hh"
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+
+using namespace std;
+using namespace CLHEP;
+using namespace AnnularTelescope_Utils;
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+namespace AnnularTelescope_NS{
+// Energy and time Resolution
+// const double EnergyThreshold = 0.*MeV;
+const double CsIResoTime = 10*ns ; // ??? (doesn't really matter)
+const double CsIResoEnergy = 0.01; // percent, from Wilton
+const double SiResoTime = 1*ns ; // ??? (doesn't really matter)
+const double SiResoEnergy = 0.0149; // absolute, from AnnularS1.hh
+// const double Radius = 50*mm ; // not used
+// const double Width = 100*mm ; // not used
+// const double Thickness = 10*mm ; // from Wilton
+// const string Material = "CsI";
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// AnnularTelescope Specific Method
+AnnularTelescope::AnnularTelescope(){
+ m_Event = new TAnnularTelescopeData() ;
+ m_AnnularTelescopeScorer_CsI= 0;
+ m_AnnularTelescopeScorer_Si = 0;
+
+ // RGB Color + Transparency
+ m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 0, 0.5));
+ m_VisCylinder = new G4VisAttributes(G4Colour(0, 0, 1, 0.5));
+}
+
+AnnularTelescope::~AnnularTelescope(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void AnnularTelescope::AddDetector(AnnularTelescope_Utils::Geometry& geo){
+ m_Geo.push_back(geo);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+std::pair<G4LogicalVolume*, G4LogicalVolume*>
+AnnularTelescope::BuildDetector(unsigned short i){
+
+ if(!(i<m_Geo.size())){
+ G4cerr << "ERROR (AnnularTelescope::BuildDetector):"
+ << "Invalid index: " << i << endl;
+ exit(1);
+ }
+
+ std::pair<G4LogicalVolume*, G4LogicalVolume*> logic; // output: <Si, CsI>
+ //
+ // Si Detector
+ {
+ std::stringstream name;
+ name << "AnnularTelescope_Si_" << i+1;
+ G4Tubs* tub = new G4Tubs(name.str().c_str(),
+ m_Geo[i].R_min, m_Geo[i].R_max,
+ m_Geo[i].SiThickness*0.5,
+ 0, 360);
+
+ G4Material* DetectorMaterial =
+ MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
+ logic.first = new G4LogicalVolume(
+ tub,DetectorMaterial,"logic_AnnularTelescope_Si_tub",0,0,0);
+ logic.first->SetVisAttributes(m_VisSquare);
+ logic.first->SetSensitiveDetector(m_AnnularTelescopeScorer_Si);
+ }
+ //
+ // CsI Detector
+ {
+ std::stringstream name;
+ name << "AnnularTelescope_CsI_" << i+1;
+ G4Tubs* tub = new G4Tubs(name.str().c_str(),
+ m_Geo[i].R_min, m_Geo[i].R_max,
+ m_Geo[i].CsIThickness*0.5,
+ 0, 360);
+
+ G4Material* DetectorMaterial =
+ MaterialManager::getInstance()->GetMaterialFromLibrary("CsI");
+ logic.second = new G4LogicalVolume(
+ tub,DetectorMaterial,"logic_AnnularTelescope_CsI_tub",0,0,0);
+ logic.second->SetVisAttributes(m_VisSquare);
+ logic.second->SetSensitiveDetector(m_AnnularTelescopeScorer_CsI);
+ }
+
+ return logic; // <Si, CsI>
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of NPS::VDetector class
+
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// Called in DetecorConstruction::ReadDetextorConfiguration Method
+void AnnularTelescope::ReadConfiguration(NPL::InputParser parser){
+ vector<Geometry> geometries =
+ AnnularTelescope_Utils::ReadConfiguration(parser);
+ for(auto& geo : geometries){ AddDetector(geo); }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Construct detector and inialise sensitive part.
+// Called After DetecorConstruction::AddDetector Method
+void AnnularTelescope::ConstructDetector(G4LogicalVolume* world){
+ for (unsigned short i = 0 ; i < m_Geo.size() ; i++) {
+
+ auto logic = BuildDetector(i); // <Si, CsI>
+ //
+ // Si
+ // Put FACE of detector at m_Geo[i].Z
+ double Z_Si = m_Geo[i].Z + m_Geo[i].SiThickness*0.5;
+ G4ThreeVector Si_pos (0, 0, Z_Si);
+ new G4PVPlacement(
+ 0, Si_pos, logic.first,"AnnularTelescope",world,false,i+1 );
+ //
+ // CsI
+ // Put FACE of detector at m_Geo[i].Z + Si thickness
+ double Z_CsI =
+ m_Geo[i].Z + m_Geo[i].SiThickness + m_Geo[i].CsIThickness*0.5;
+ G4ThreeVector CsI_pos (0, 0, Z_CsI);
+ new G4PVPlacement(
+ 0, CsI_pos, logic.second,"AnnularTelescope",world,false,i+1 );
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void AnnularTelescope::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ if(!pTree->FindBranch("AnnularTelescope")){
+ pTree->Branch("AnnularTelescope", "TAnnularTelescopeData", &m_Event) ;
+ }
+ pTree->SetBranchAddress("AnnularTelescope", &m_Event) ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read scorer
+// Helper for ReadSensitive()
+void AnnularTelescope::ReadScorer(
+ const G4Event* event, const char* scorerName,
+ std::vector<HitInfo_t>& hits) {
+
+ NPS::HitsMap<G4double*>* CaloHitMap;
+ std::map<G4int, G4double**>::iterator Calo_itr;
+
+ G4int CaloCollectionID =
+ G4SDManager::GetSDMpointer()->GetCollectionID(scorerName);
+ CaloHitMap = (NPS::HitsMap<G4double*>*)(
+ event->GetHCofThisEvent()->GetHC(CaloCollectionID) );
+
+ // Loop on the Calo map
+ for(const auto& hit : *(CaloHitMap->GetMap()) ) {
+ // Read Hit Information
+ // Infos[0] = aStep->GetTotalEnergyDeposit();
+ // Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+ // Infos[2] = m_Position.x();
+ // Infos[3] = m_Position.y();
+ // Infos[4] = m_Position.z();
+ // Infos[5] = m_Position.theta();
+ // Infos[6] = m_Position.phi();
+ // Infos[7] = Detector Number
+ hits.push_back(HitInfo_t());
+ G4double* Infos = *(hit.second);
+
+ hits.back().detector = int(Infos[7]);
+ hits.back().energy = Infos[0];
+ hits.back().time = Infos[1];
+ hits.back().x = Infos[2];
+ hits.back().y = Infos[3];
+ hits.back().z = Infos[4];
+ }
+ // clear map for next event
+ CaloHitMap->clear();
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void AnnularTelescope::ReadSensitive(const G4Event* event){
+ m_Event->Clear();
+
+ // CSI
+ std::vector<HitInfo_t> csiHits;
+ ReadScorer(event, "AnnularTelescopeScorer_CsI/Calorimeter_CsI", csiHits);
+ for(const auto& hit : csiHits) {
+ FillCsIData(hit.detector, hit.energy, hit.time, G4ThreeVector(hit.x,hit.y,hit.z));
+ }
+
+ // SI
+ std::vector<HitInfo_t> siHits;
+ ReadScorer(event, "AnnularTelescopeScorer_Si/Calorimeter_Si", siHits);
+ for(const auto& hit : siHits) {
+ FillSiData(hit.detector, hit.energy, hit.time, G4ThreeVector(hit.x,hit.y,hit.z));
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+void AnnularTelescope::FillCsIData(
+ int detector_number, double energy, double time, const G4ThreeVector& pos){
+ //
+ // Figure out wedge
+ int wedge_number = 0;
+ if(size_t(detector_number - 1) >= m_Geo.size()) {
+ std::cerr << "ERROR: AnnularTelescope::FillCsIData: "
+ << "invalid detector number: " << detector_number
+ << ", tried to access index: " << detector_number-1
+ << " in length " << m_Geo.size() << " array\n";
+ exit(1);
+ }
+ const Geometry& geo = m_Geo.at(detector_number - 1);
+ double pitch = geo.CsI_Wedge_Angle_Pitch;
+ for(const double& phi : geo.CsI_Wedge_Phi_Angle){
+ if(pos.phi() >= phi - pitch/2. && pos.phi() < phi + pitch/2.) {
+ // match!
+ break;
+ }
+ ++wedge_number;
+ }
+ // start counting at 1, not 0
+ ++wedge_number;
+ //
+ // Add resolutions
+ double eres = AnnularTelescope_NS::CsIResoEnergy*energy; // percent
+ double energy_res = -1000;
+ while(energy_res < 0) { // make sure it's a positive number
+ energy_res = RandGauss::shoot(energy, eres);
+ }
+ double tres = AnnularTelescope_NS::CsIResoTime; // absolute
+ double time_res = RandGauss::shoot(time, tres);
+ //
+ // Set CsI energy and time
+ m_Event->SetCsIEnergy(detector_number, wedge_number, energy_res);
+ m_Event->SetCsITime(detector_number, wedge_number, time_res);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+void AnnularTelescope::FillSiData(
+ int detector_number, double energy, double time, const G4ThreeVector& pos){
+ //
+ // Figure out strips
+ // Phi
+ if(size_t(detector_number - 1) >= m_Geo.size()) {
+ std::cerr << "ERROR: AnnularTelescope::FillSiData: "
+ << "invalid detector number: " << detector_number
+ << ", tried to access index: " << detector_number-1
+ << " in length " << m_Geo.size() << " array\n";
+ exit(1);
+ }
+ const Geometry& geo = m_Geo.at(detector_number - 1);
+ size_t phi_strip = 1;
+ {
+ double pitch = geo.Si_Phi_Angle_Pitch;
+ for(const double& phi : geo.Si_Strip_Phi_Angle){
+ if(pos.phi() >= phi - pitch/2. && pos.phi() < phi + pitch/2.) {
+ // match!
+ break;
+ }
+ ++phi_strip;
+ }
+ if(phi_strip == geo.Si_Strip_Phi_Angle.size() + 1) {
+ std::cerr << "\nWARNING: Phi strip number not found: "
+ << "Angle of hit [deg]: " << pos.phi()/deg << "\n";
+ std::cerr << "Strip angles:\n";
+ for(const double& R : geo.Si_Strip_Phi_Angle){
+ std::cerr << " " << (R-pitch/2.)/deg << " - " << (R+pitch/2.)/deg << "\n";
+ }
+ std::cerr << "---------\n";
+ }
+ }
+ // Theta
+ size_t theta_strip = 1;
+ {
+ double pitch = geo.Si_Theta_Radius_Pitch;
+ for(const double& R : geo.Si_Strip_Theta_Radius){
+ if(pos.perp() >= R - pitch/2. && pos.perp() < R + pitch/2.) {
+ // match!
+ break;
+ }
+ ++theta_strip;
+ }
+ if(theta_strip == geo.Si_Strip_Theta_Radius.size() + 1) {
+ // Strip not found
+ // Check for error in floating point
+ if(
+ pos.perp() < geo.Si_Strip_Theta_Radius.front() - pitch/2. &&
+ fabs(pos.perp() - (geo.Si_Strip_Theta_Radius.front() - pitch/2.)) < 1e-3) {
+ // okay!
+ theta_strip = 1; // inner strip
+ }
+ else if (
+ pos.perp() >= geo.Si_Strip_Theta_Radius.back() + pitch/2. &&
+ fabs(pos.perp() - (geo.Si_Strip_Theta_Radius.back() + pitch/2.)) < 1e-3) {
+ // okay!
+ --theta_strip; // outer strip
+ }
+ else {
+ //
+ // Genuine problem
+ std::cerr << "\nWARNING: Theta strip number not found: "
+ << "Radius of hit [mm]: " << pos.perp()/mm << "\n";
+ std::cerr << "Strip radii:\n";
+ for(const double& R : geo.Si_Strip_Theta_Radius){
+ std::cout << " " << (R-pitch/2.)/mm << " - " << (R+pitch/2.)/mm << "\n";
+ }
+ std::cerr << "---------\n";
+ }
+ }
+ }
+ //
+ // Add resolutions
+ double eres = AnnularTelescope_NS::SiResoEnergy; // absolute
+ double energy_theta = -1000;
+ while(energy_theta < 0) { // make sure it's a positive number
+ energy_theta = RandGauss::shoot(energy, eres);
+ }
+ double energy_phi = -1000;
+ while(energy_phi < 0) { // make sure it's a positive number
+ energy_phi = RandGauss::shoot(energy, eres);
+ }
+ double tres = AnnularTelescope_NS::SiResoTime; // absolute
+ double time_theta = RandGauss::shoot(time, tres);
+ double time_phi = RandGauss::shoot(time, tres);
+ //
+ // Set Si energy and time
+ m_Event->SetSiThetaEnergy(detector_number, theta_strip, energy_theta);
+ m_Event->SetSiPhiEnergy(detector_number, phi_strip, energy_phi);
+ m_Event->SetSiThetaTime(detector_number, theta_strip, time_theta);
+ m_Event->SetSiPhiTime(detector_number, phi_strip, time_phi);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+void AnnularTelescope::InitializeScorers() {
+ // CSI
+ //
+ // This check is necessary in case the geometry is reloaded
+ bool already_exist = false;
+ m_AnnularTelescopeScorer_CsI =
+ CheckScorer("AnnularTelescopeScorer_CsI",already_exist) ;
+
+ if(!already_exist) {
+ // Otherwise the scorer is initialised
+ vector<int> level; level.push_back(0);
+ G4VPrimitiveScorer* Calorimeter_CsI =
+ new CALORIMETERSCORERS::PS_CalorimeterWithInteraction(
+ "Calorimeter_CsI",level, 0) ;
+ //and register it to the multifunctionnal detector
+ m_AnnularTelescopeScorer_CsI->RegisterPrimitive(Calorimeter_CsI);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_AnnularTelescopeScorer_CsI) ;
+ }
+ // SI
+ //
+ // This check is necessary in case the geometry is reloaded
+ already_exist = false;
+ m_AnnularTelescopeScorer_Si =
+ CheckScorer("AnnularTelescopeScorer_Si",already_exist) ;
+
+ if(!already_exist) {
+ // Otherwise the scorer is initialised
+ vector<int> level; level.push_back(0);
+ G4VPrimitiveScorer* Calorimeter_Si =
+ new CALORIMETERSCORERS::PS_CalorimeterWithInteraction(
+ "Calorimeter_Si",level, 0) ;
+ //and register it to the multifunctionnal detector
+ m_AnnularTelescopeScorer_Si->RegisterPrimitive(Calorimeter_Si);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_AnnularTelescopeScorer_Si) ;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPS::VDetector* AnnularTelescope::Construct(){
+ return (NPS::VDetector*) new AnnularTelescope();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern"C" {
+class proxy_nps_AnnularTelescope{
+public:
+ proxy_nps_AnnularTelescope(){
+ NPS::DetectorFactory::getInstance()->AddToken("AnnularTelescope","AnnularTelescope");
+ NPS::DetectorFactory::getInstance()->AddDetector("AnnularTelescope",AnnularTelescope::Construct);
+ }
+};
+
+proxy_nps_AnnularTelescope p_nps_AnnularTelescope;
+}
diff --git a/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.hh b/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e240ae049cc5a993629a452c752cc44889d4d816
--- /dev/null
+++ b/NPSimulation/Detectors/AnnularTelescope/AnnularTelescope.hh
@@ -0,0 +1,126 @@
+#ifndef AnnularTELESCOPE_SIMULATION_HEADER
+#define AnnularTELESCOPE_SIMULATION_HEADER 1
+/*****************************************************************************
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * *
+ * Creation Date : March 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe AnnularTelescope simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ header
+#include <string>
+#include <vector>
+#include <utility>
+
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+// NPTool header
+#include "NPSVDetector.hh"
+#include "AnnularTelescope_Utils.h"
+#include "TAnnularTelescopeData.h"
+#include "NPInputParser.h"
+
+namespace AnnularTelescope_Utils { struct Geometry; }
+
+class AnnularTelescope : public NPS::VDetector{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+public:
+ AnnularTelescope() ;
+ virtual ~AnnularTelescope() ;
+
+ ////////////////////////////////////////////////////
+ /////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // Cartesian
+ void AddDetector(AnnularTelescope_Utils::Geometry& geo);
+
+ std::pair<G4LogicalVolume*, G4LogicalVolume*>
+ BuildDetector(unsigned short i);
+
+private:
+ ////////////////////////////////////////////////////
+ ////// Inherite from NPS::VDetector class /////////
+ ////////////////////////////////////////////////////
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(NPL::InputParser) ;
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world) ;
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event) ;
+
+ // Read scorer information
+ struct HitInfo_t {
+ G4int detector; G4double energy; G4double time;
+ G4double x; G4double y; G4double z;
+ };
+ void ReadScorer(const G4Event* event, const char* scorerName,
+ std::vector<HitInfo_t>& hits);
+
+ // Fill NPL CsI Data
+ // Called from ReadSensitive
+ void FillCsIData(int detector_number, double energy, double time, const G4ThreeVector& pos);
+
+ // Fill NPL Si Data
+ // Called from ReadSensitive
+ void FillSiData(int detector_number, double energy, double time, const G4ThreeVector& pos);
+
+public: // Scorer
+ // Initialize all Scorer used by the MUST2Array
+ void InitializeScorers() ;
+
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_AnnularTelescopeScorer_CsI ;
+ G4MultiFunctionalDetector* m_AnnularTelescopeScorer_Si ;
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+private:
+ TAnnularTelescopeData* m_Event ;
+
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private: // Geometry
+ // Detector Coordinates
+ std::vector<AnnularTelescope_Utils::Geometry> m_Geo;
+
+ // Visualisation Attribute
+ G4VisAttributes* m_VisSquare;
+ G4VisAttributes* m_VisCylinder;
+
+ // Needed for dynamic loading of the library
+public:
+ static NPS::VDetector* Construct();
+};
+#endif
diff --git a/NPSimulation/Detectors/AnnularTelescope/CMakeLists.txt b/NPSimulation/Detectors/AnnularTelescope/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..111a7c565e75d6122c0998c46e0f6c4e4ec19cb8
--- /dev/null
+++ b/NPSimulation/Detectors/AnnularTelescope/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_library(NPSAnnularTelescope SHARED AnnularTelescope.cc)
+target_link_libraries(NPSAnnularTelescope NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPAnnularTelescope)
diff --git a/NPSimulation/Detectors/MDM/MDM.cc b/NPSimulation/Detectors/MDM/MDM.cc
index 22d6b50b2018c809ff8be181f4c61318fc919e84..37bfcbc60c2c5844b8bb8db42c59d1e2b54ca82d 100644
--- a/NPSimulation/Detectors/MDM/MDM.cc
+++ b/NPSimulation/Detectors/MDM/MDM.cc
@@ -62,8 +62,8 @@ using namespace CLHEP;
namespace MDM_NS{
// Energy and time Resolution
const double Width = 250*mm ;
-const double Thickness = 10*mm;
-const double Zpos = 40*cm;
+const double Thickness = 1*m; // needs to be very thick to serve as calorimeter
+const double Zpos = 40*cm; // of FRONT face
const string Material = "BC400"; // fake!!
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -93,7 +93,7 @@ void MDM::AddDetector(double angle, double field, double xaccept, double yaccept
m_Rayin = new MDMTrace::Rayin(m_Rayin_file, false);
m_Trace = MDMTrace::Instance();
- m_Trace->SetMDMAngle(angle/mrad); // mrad
+ m_Trace->SetMDMAngle(angle/deg); // mrad
m_Trace->SetMDMDipoleField(field/gauss); // gauss
cout << "MDM::AddDetector :: Angle [mrad], Angle [deg], Field [G], Rayin File :: "
@@ -154,7 +154,7 @@ void MDM::ReadConfiguration(NPL::InputParser parser){
void MDM::ConstructDetector(G4LogicalVolume* world){
G4double wX = 0;
G4double wY = 0;
- G4double wZ = MDM_NS::Zpos;
+ G4double wZ = MDM_NS::Zpos + MDM_NS::Thickness;
G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ;
new G4PVPlacement(0, Det_pos, BuildSquareDetector(),
@@ -203,11 +203,11 @@ void MDM::ReadSensitive(const G4Event* event){
// check if within acceptance
// saves lots of time not tracking events outside of the
// acceptance
- if(fabs(thetaX) < m_Xaccept && fabs(thetaY) < m_Yaccept)
+ if((fabs(thetaX) < m_Xaccept && fabs(thetaY) < m_Yaccept))
{
// Calculate positions at TARGET
- double xTrgt = Pos.x()/mm - (MDM_NS::Zpos/mm - MDM_NS::Thickness*0.5/mm)*tan(thetaX);
- double yTrgt = Pos.y()/mm - (MDM_NS::Zpos/mm - MDM_NS::Thickness*0.5/mm)*tan(thetaY);
+ double xTrgt = Pos.x()/mm - (Pos.z()/mm)*tan(thetaX);
+ double yTrgt = Pos.y()/mm - (Pos.z()/mm)*tan(thetaY);
double zTrgt = 0.*mm;
// Send Through MDM
diff --git a/NPSimulation/Detectors/Tiara/Tiara.cc b/NPSimulation/Detectors/Tiara/Tiara.cc
index 97b0b4128cc57c424536058f40ec738c6ffd259e..70c84f65ffd48362ce6e93014ae94734e55eb9a5 100644
--- a/NPSimulation/Detectors/Tiara/Tiara.cc
+++ b/NPSimulation/Detectors/Tiara/Tiara.cc
@@ -27,7 +27,7 @@
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Cons.hh"
-#include "G4UnionSolid.hh"
+#include "G4UnionSolid.hh"
#include "G4ExtrudedSolid.hh"
#include "G4TwoVector.hh"
//G4 sensitive
@@ -74,14 +74,14 @@ Tiara::Tiara(){
// Light Grey
FrameVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
// Light Blue
- GuardRingVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ;
+ GuardRingVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ;
m_boolChamber = false;
m_boolInner = false;
m_boolOuter = false;
m_InnerBarrelScorer = 0 ;
- m_OuterBarrelScorer = 0 ;
- m_HyballScorer = 0 ;
+ m_OuterBarrelScorer = 0 ;
+ m_HyballScorer = 0 ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -96,7 +96,7 @@ void Tiara::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("Tiara","Barrel");
if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
+ cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> token = {"InnerBarrel","OuterBarrel","Chamber"};
@@ -116,10 +116,10 @@ void Tiara::ReadConfiguration(NPL::InputParser parser){
}
blocks.clear();
- blocks = parser.GetAllBlocksWithToken("HyballWedge");
+ blocks = parser.GetAllBlocksWithToken("TiaraHyballWedge");
if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
+ cout << "//// " << blocks.size() << " detectors found " << endl;
token.clear();
token = {"Z","R","Phi"};
@@ -172,29 +172,29 @@ void Tiara::ReadSensitive(const G4Event* event){
G4int InnerBarrelCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Tiara_InnerBarrelScorer/InnerBarrel");
InnerBarrelHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(InnerBarrelCollectionID));
- // Loop on the InnerBarrel map
+ // Loop on the InnerBarrel map
for (InnerBarrel_itr = InnerBarrelHitMap->GetMap()->begin() ; InnerBarrel_itr != InnerBarrelHitMap->GetMap()->end() ; InnerBarrel_itr++){
- G4double* Info = *(InnerBarrel_itr->second);
+ G4double* Info = *(InnerBarrel_itr->second);
// Upstream Energy
double EU = RandGauss::shoot(Info[0],ResoEnergyInnerBarrel);
if(EU>EnergyThreshold){
m_EventBarrel->SetFrontUpstreamE(Info[3],Info[4],EU);
- m_EventBarrel->SetFrontUpstreamT(Info[3],Info[4],Info[2]);
+ m_EventBarrel->SetFrontUpstreamT(Info[3],Info[4],Info[2]);
}
// Downstream Energy
- double ED = RandGauss::shoot(Info[1],ResoEnergyInnerBarrel);
+ double ED = RandGauss::shoot(Info[1],ResoEnergyInnerBarrel);
if(ED>EnergyThreshold){
m_EventBarrel->SetFrontDownstreamE(Info[3],Info[4],ED);
- m_EventBarrel->SetFrontDownstreamT(Info[3],Info[4],Info[2]);
+ m_EventBarrel->SetFrontDownstreamT(Info[3],Info[4],Info[2]);
}
// Back Energy
double EB = RandGauss::shoot(Info[1]+Info[0],ResoEnergyInnerBarrel);
if(EB>EnergyThreshold){
m_EventBarrel->SetBackE(Info[3],EB);
- m_EventBarrel->SetBackT(Info[3],Info[2]);
+ m_EventBarrel->SetBackT(Info[3],Info[2]);
}
// Interaction Coordinates
ms_InterCoord->SetDetectedPositionX(Info[5]) ;
@@ -212,14 +212,14 @@ void Tiara::ReadSensitive(const G4Event* event){
G4int OuterBarrelCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Tiara_OuterBarrelScorer/OuterBarrel");
OuterBarrelHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(OuterBarrelCollectionID));
- // Loop on the OuterBarrel map
+ // Loop on the OuterBarrel map
for (OuterBarrel_itr = OuterBarrelHitMap->GetMap()->begin() ; OuterBarrel_itr != OuterBarrelHitMap->GetMap()->end() ; OuterBarrel_itr++){
- G4double* Info = *(OuterBarrel_itr->second);
+ G4double* Info = *(OuterBarrel_itr->second);
double E = RandGauss::shoot(Info[0],ResoEnergyOuterBarrel);
if(E>EnergyThreshold){
m_EventBarrel->SetOuterE(Info[7],Info[9],E);
- m_EventBarrel->SetOuterT(Info[7],Info[9],Info[1]);
+ m_EventBarrel->SetOuterT(Info[7],Info[9],Info[1]);
}
// Interaction Coordinates
ms_InterCoord->SetDetectedPositionX(Info[5]) ;
@@ -238,20 +238,20 @@ void Tiara::ReadSensitive(const G4Event* event){
G4int HyballCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Tiara_HyballScorer/Hyball");
HyballHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(HyballCollectionID));
- // Loop on the Hyball map
+ // Loop on the Hyball map
for (Hyball_itr = HyballHitMap->GetMap()->begin() ; Hyball_itr != HyballHitMap->GetMap()->end() ; Hyball_itr++){
- G4double* Info = *(Hyball_itr->second);
+ G4double* Info = *(Hyball_itr->second);
// Front Energy
double EF = RandGauss::shoot(Info[0],ResoEnergyHyball);
if(EF>EnergyThreshold){
int RingNumber=Info[8];
//by Shuya 171009
- RingNumber=abs(RingNumber);
- //RingNumber=abs(RingNumber-17);
- Info[8]=RingNumber;
+ RingNumber=abs(RingNumber);
+ //RingNumber=abs(RingNumber-17);
+ Info[8]=RingNumber;
m_EventHyball->SetRingE(Info[7],Info[8],EF);
- m_EventHyball->SetRingT(Info[7],Info[8],Info[1]);
+ m_EventHyball->SetRingT(Info[7],Info[8],Info[1]);
}
// Back Energy
@@ -274,7 +274,7 @@ void Tiara::ReadSensitive(const G4Event* event){
//by Shuya 171009
//m_EventHyball->SetSectorE(Info[7],Info[9],EF);
m_EventHyball->SetSectorE(Info[7],Info[9],EB);
- m_EventHyball->SetSectorT(Info[7],Info[9],Info[1]);
+ m_EventHyball->SetSectorT(Info[7],Info[9],Info[1]);
}
// Interaction Coordinates
//by Shuya 171009
@@ -321,24 +321,24 @@ void Tiara::InitializeScorers(){
m_OuterBarrelScorer->RegisterPrimitive(OuterBarrel);
- G4VPrimitiveScorer* Hyball= new SILICONSCORERS::PS_Silicon_Annular("Hyball",1,
- HYBALL_ActiveWafer_InnerRadius,
- HYBALL_ActiveWafer_OuterRadius,
+ G4VPrimitiveScorer* Hyball= new SILICONSCORERS::PS_Silicon_Annular("Hyball",1,
+ HYBALL_ActiveWafer_InnerRadius,
+ HYBALL_ActiveWafer_OuterRadius,
-0.5*HYBALL_ActiveWafer_Angle,0.5*HYBALL_ActiveWafer_Angle,
HYBALL_NumberOfAnnularStrip,
HYBALL_NumberOfRadialStrip);
m_HyballScorer->RegisterPrimitive(Hyball);
- // Add All Scorer to the Global Scorer Manager
+ // Add All Scorer to the Global Scorer Manager
G4SDManager::GetSDMpointer()->AddNewDetector(m_InnerBarrelScorer) ;
- G4SDManager::GetSDMpointer()->AddNewDetector(m_OuterBarrelScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_OuterBarrelScorer) ;
G4SDManager::GetSDMpointer()->AddNewDetector(m_HyballScorer) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Tiara::InitializeRootOutput(){
- TTree *pTree = RootOutput::getInstance()->GetTree();
+ TTree *pTree = RootOutput::getInstance()->GetTree();
if(!pTree->FindBranch("TiaraBarrel")){
pTree->Branch("TiaraBarrel", "TTiaraBarrelData", &m_EventBarrel) ;
}
@@ -346,7 +346,7 @@ void Tiara::InitializeRootOutput(){
if(!pTree->FindBranch("TiaraHyball")){
pTree->Branch("TiaraHyball", "TTiaraHyballData", &m_EventHyball) ;
}
- // This insure that the object are correctly bind in case of
+ // This insure that the object are correctly bind in case of
// a redifinition of the geometry in the simulation
pTree->SetBranchAddress("TiaraBarrel", &m_EventBarrel) ;
pTree->SetBranchAddress("TiaraHyball", &m_EventHyball) ;
@@ -375,7 +375,7 @@ void Tiara::ConstructInnerBarrel(G4LogicalVolume* world){
PCBCrossSection.push_back(G4TwoVector(-INNERBARREL_PCB_Width/2.,0));
PCBCrossSection.push_back(G4TwoVector(-INNERBARREL_PCB_Width/2.+l2,-INNERBARREL_PCB_Thickness*0.5));
- G4ExtrudedSolid* PCBFull =
+ G4ExtrudedSolid* PCBFull =
new G4ExtrudedSolid("PCBFull",
PCBCrossSection,
INNERBARREL_PCB_Length/2.,
@@ -473,7 +473,7 @@ void Tiara::ConstructInnerBarrel(G4LogicalVolume* world){
G4ThreeVector WaferPosition(0,
0.5*(INNERBARREL_PCB_Thickness+INNERBARREL_InertWafer_Thickness)-INNERBARREL_PCB_WaferDepth
- ,0);
+ ,0);
G4ThreeVector DeadLayerPositionF = WaferPosition + G4ThreeVector(0,-INNERBARREL_ActiveWafer_Thickness*0.5-INNERBARREL_ActiveWafer_DeadLayerThickness*0.5,0);
@@ -504,21 +504,21 @@ void Tiara::ConstructInnerBarrel(G4LogicalVolume* world){
// The Distance from target is given by half the lenght of a detector
// plus the length of a detector inclined by 45 deg.
- G4double DistanceFromTarget = INNERBARREL_PCB_Width*(0.5+sin(45*deg)) ;
+ G4double DistanceFromTarget = INNERBARREL_PCB_Width*(0.5+sin(45*deg)) ;
for( unsigned int i = 0; i < 8; i ++){
// The following build the barrel, with detector one at the top
// and going clowise looking upstrea
- // Detector are rotate by 45deg with each other
- // Detector 3 [i=2] is perpendicular to positive y-axis, Detector 5 [i=4] is perpendicular to positive x-axis,
- G4RotationMatrix* DetectorRotation =
- new G4RotationMatrix(0*deg,0*deg,(2-i)*45*deg);
+ // Detector are rotate by 45deg with each other
+ // Detector 3 [i=2] is perpendicular to positive y-axis, Detector 5 [i=4] is perpendicular to positive x-axis,
+ G4RotationMatrix* DetectorRotation =
+ new G4RotationMatrix(0*deg,0*deg,(360+(2-i)*45)*deg);
// There center is also rotated by 45deg
G4ThreeVector DetectorPosition(0,DistanceFromTarget,0);
- DetectorPosition.rotate((2-i)*45*deg,G4ThreeVector(0,0,1));
+ DetectorPosition.rotate((360+(2-i)*-45)*deg,G4ThreeVector(0,0,1));
- // Place the Master volume with its two daugther volume at the final place
+ // Place the Master volume with its two daugther volume at the final place
new G4PVPlacement(G4Transform3D(*DetectorRotation,DetectorPosition),
logicBarrelDetector,"Tiara_InnerBarrel_Detector",
world,false,i+1);
@@ -544,7 +544,7 @@ void Tiara::ConstructOuterBarrel(G4LogicalVolume* world){
PCBCrossSection.push_back(G4TwoVector(-OUTERBARREL_PCB_Width/2.,0));
PCBCrossSection.push_back(G4TwoVector(-OUTERBARREL_PCB_Width/2.+l2,-OUTERBARREL_PCB_Thickness*0.5));
- G4ExtrudedSolid* PCBFull =
+ G4ExtrudedSolid* PCBFull =
new G4ExtrudedSolid("PCBFull",
PCBCrossSection,
OUTERBARREL_PCB_Length/2.,
@@ -625,7 +625,7 @@ void Tiara::ConstructOuterBarrel(G4LogicalVolume* world){
// The Distance from target is given by half the lenght of a detector
// plus the length of a detector inclined by 45 deg.
- G4double DistanceFromTarget = OUTERBARREL_PCB_Width*(0.5+sin(45*deg)) ;
+ G4double DistanceFromTarget = OUTERBARREL_PCB_Width*(0.5+sin(45*deg)) ;
// Place the sub volumes in the master volume
// Last argument is the detector number, used in the scorer to get the
@@ -638,7 +638,7 @@ void Tiara::ConstructOuterBarrel(G4LogicalVolume* world){
G4ThreeVector WaferPosition(0,
0.5*(OUTERBARREL_PCB_Thickness+OUTERBARREL_ActiveWafer_Thickness)-OUTERBARREL_PCB_WaferDepth
- ,0);
+ ,0);
new G4PVPlacement(new G4RotationMatrix(0,0,0),
WaferPosition,
@@ -654,15 +654,15 @@ void Tiara::ConstructOuterBarrel(G4LogicalVolume* world){
// The following build the barrel, with detector one at the top
// and going clowise looking upstrea
for( unsigned int i = 0; i < 8; i ++){
- // Detector are rotate by 45deg with each other
- G4RotationMatrix* DetectorRotation =
+ // Detector are rotate by 45deg with each other
+ G4RotationMatrix* DetectorRotation =
new G4RotationMatrix(0*deg,0*deg,i*45*deg);
// There center is also rotated by 45deg
G4ThreeVector DetectorPosition(0,DistanceFromTarget,0);
- DetectorPosition.rotate(i*45*deg,G4ThreeVector(0,0,-1));
+ DetectorPosition.rotate(i*45*deg,G4ThreeVector(0,0,-1));
- // Place the Master volume with its daugthers at the final place
+ // Place the Master volume with its daugthers at the final place
new G4PVPlacement(G4Transform3D(*DetectorRotation,DetectorPosition),
logicBarrelDetector,"Tiara_OuterBarrel_Detector",
world,false,i+1);
@@ -679,7 +679,7 @@ void Tiara::ConstructHyball(G4LogicalVolume* world){
PCBCrossSection.push_back(G4TwoVector(125.718*mm, 73.677*mm));
PCBCrossSection.push_back(G4TwoVector(28.108*mm , 16.473*mm));
- G4ExtrudedSolid* PCBFull =
+ G4ExtrudedSolid* PCBFull =
new G4ExtrudedSolid("PCBFull",
PCBCrossSection,
0.5*HYBALL_PCB_THICKNESS,
@@ -694,14 +694,14 @@ void Tiara::ConstructHyball(G4LogicalVolume* world){
WaferCrossSection.push_back(G4TwoVector(122.677*mm, 63.508*mm ));
WaferCrossSection.push_back(G4TwoVector(29.108*mm , 15.069*mm));
- G4ExtrudedSolid* WaferFull =
+ G4ExtrudedSolid* WaferFull =
new G4ExtrudedSolid("WaferFull",
WaferCrossSection,
0.5*HYBALL_ActiveWafer_Thickness,
G4TwoVector(0,0),1,
G4TwoVector(0,0),1);
- G4ExtrudedSolid* WaferShape =
+ G4ExtrudedSolid* WaferShape =
new G4ExtrudedSolid("WaferShape",
WaferCrossSection,
0.6*HYBALL_PCB_THICKNESS,
@@ -709,18 +709,18 @@ void Tiara::ConstructHyball(G4LogicalVolume* world){
G4TwoVector(0,0),1);
// Active Wafer
- G4Tubs* ActiveWafer =
+ G4Tubs* ActiveWafer =
new G4Tubs("HyballActiveWafer",HYBALL_ActiveWafer_InnerRadius,
HYBALL_ActiveWafer_OuterRadius,0.5*HYBALL_ActiveWafer_Thickness,
-0.5*HYBALL_ActiveWafer_Angle,HYBALL_ActiveWafer_Angle);
- G4Tubs* ActiveWaferShape =
+ G4Tubs* ActiveWaferShape =
new G4Tubs("HyballActiveWaferShape",HYBALL_ActiveWafer_InnerRadius,
HYBALL_ActiveWafer_OuterRadius,0.6*HYBALL_ActiveWafer_Thickness,
-0.5*HYBALL_ActiveWafer_Angle,HYBALL_ActiveWafer_Angle);
//by Shuya 180219. Hyball DeadLayer
- G4Tubs* Hyball_DeadLayer =
+ G4Tubs* Hyball_DeadLayer =
new G4Tubs("HyballDeadLayer",HYBALL_ActiveWafer_InnerRadius,
HYBALL_ActiveWafer_OuterRadius,0.5*HYBALL_DeadLayer_Thickness,
-0.5*HYBALL_ActiveWafer_Angle,HYBALL_ActiveWafer_Angle);
@@ -803,25 +803,25 @@ void Tiara::ConstructChamber(G4LogicalVolume* world){
// Hyball is hold on a back plate that close the Diabolo Shaped Chamber
// Making the Chamber //
- // We make the individual pieces, starting from the inside to the outside
+ // We make the individual pieces, starting from the inside to the outside
// Then we merge them together using the a G4AdditionSolid
- // The whole chamber is then placed
+ // The whole chamber is then placed
// Central Tube
- G4Tubs* solidCentralTube =
+ G4Tubs* solidCentralTube =
new G4Tubs("TiaraChamberCentralTube",CHAMBER_CentralTube_Inner_Radius,
CHAMBER_CentralTube_Outer_Radius,CHAMBER_CentralTube_Length/2.,
0*deg,360*deg);
// Forward-Backward Cones
- G4Cons* solidOuterCone =
+ G4Cons* solidOuterCone =
new G4Cons("TiaraChamberOuterCone",CHAMBER_CentralTube_Inner_Radius,
CHAMBER_CentralTube_Outer_Radius,CHAMBER_OuterCylinder_Inner_Radius,
CHAMBER_OuterCylinder_Outer_Radius,CHAMBER_OuterCone_Length/2.,
0*deg,360*deg);
// Outer Cylinder
- G4Tubs* solidOuterCylinder =
+ G4Tubs* solidOuterCylinder =
new G4Tubs("TiaraChamberOuterCylinder",CHAMBER_OuterCylinder_Inner_Radius,
CHAMBER_OuterCylinder_Outer_Radius,CHAMBER_OuterCylinder_Length/2.,
0*deg,360*deg);
@@ -830,22 +830,22 @@ void Tiara::ConstructChamber(G4LogicalVolume* world){
G4UnionSolid* solidTiaraChamberStep1 =
new G4UnionSolid("TiaraChamber", solidCentralTube, solidOuterCone,
new G4RotationMatrix,
- G4ThreeVector(0,0,CHAMBER_OuterCone_Z_Pos));
+ G4ThreeVector(0,0,CHAMBER_OuterCone_Z_Pos));
G4UnionSolid* solidTiaraChamberStep2 =
new G4UnionSolid("TiaraChamber", solidTiaraChamberStep1, solidOuterCone,
new G4RotationMatrix(0,180*deg,0),
- G4ThreeVector(0,0,-CHAMBER_OuterCone_Z_Pos));
+ G4ThreeVector(0,0,-CHAMBER_OuterCone_Z_Pos));
G4UnionSolid* solidTiaraChamberStep3 =
new G4UnionSolid("TiaraChamber", solidTiaraChamberStep2, solidOuterCylinder,
new G4RotationMatrix,
- G4ThreeVector(0,0,CHAMBER_OuterCylinder_Z_Pos));
+ G4ThreeVector(0,0,CHAMBER_OuterCylinder_Z_Pos));
G4UnionSolid* solidTiaraChamberStep4 =
new G4UnionSolid("TiaraChamber", solidTiaraChamberStep3, solidOuterCylinder,
new G4RotationMatrix,
- G4ThreeVector(0,0,-CHAMBER_OuterCylinder_Z_Pos));
+ G4ThreeVector(0,0,-CHAMBER_OuterCylinder_Z_Pos));
// Create Logic Volume
G4LogicalVolume* logicTiaraChamber =
@@ -865,7 +865,7 @@ void Tiara::ConstructChamber(G4LogicalVolume* world){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Tiara::InitializeMaterial(){
m_MaterialSilicon = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
- m_MaterialAl = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ m_MaterialAl = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
m_MaterialPCB = MaterialManager::getInstance()->GetMaterialFromLibrary("PCB");
m_MaterialVacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
}
@@ -891,8 +891,3 @@ class proxy_nps_tiara{
proxy_nps_tiara p_nps_tiara;
}
-
-
-
-
-
diff --git a/NPSimulation/Detectors/Tiara/Tiara.hh b/NPSimulation/Detectors/Tiara/Tiara.hh
index e31169cd500b7ee46516eeb765f22ccea738e4df..118de13014d80548f3dc74a76ba4c45b651c3525 100644
--- a/NPSimulation/Detectors/Tiara/Tiara.hh
+++ b/NPSimulation/Detectors/Tiara/Tiara.hh
@@ -48,9 +48,9 @@ using namespace CLHEP;
namespace TIARA{
// Energy and time Resolution
const G4double ResoTime = 0 ;
- //const G4double ResoEnergyInnerBarrel = 0.017*MeV ;// = 136keV FWHM
+ //const G4double ResoEnergyInnerBarrel = 0.017*MeV ;// = 136keV FWHM
//const G4double ResoEnergyOuterBarrel = 0.017*MeV ;// = 136keV FWHM
- const G4double ResoEnergyInnerBarrel = 0.050*MeV ;// = 136keV FWHM
+ const G4double ResoEnergyInnerBarrel = 0.050*MeV ;// = 136keV FWHM
const G4double ResoEnergyOuterBarrel = 0.050*MeV ;// = 136keV FWHM
const G4double ResoEnergyHyball = 0.017*MeV ;// = 70keV FWHM
@@ -62,7 +62,7 @@ namespace TIARA{
const G4double CHAMBER_CentralTube_Inner_Radius = 4.86*cm; //4.05->Original Value for the Single stage barrel
const G4double CHAMBER_CentralTube_Outer_Radius = 5.05*cm; //4.25->Original Value for the Single stage barrel
const G4double CHAMBER_CentralTube_Length = 8.24*cm;
-
+
// Outer Cone
const G4double CHAMBER_OuterCone_Length = 9.88*cm;
const G4double CHAMBER_OuterCone_Z_Pos = 9.06*cm;
@@ -75,26 +75,26 @@ namespace TIARA{
// Inner Barrel //
const G4double INNERBARREL_PCB_Length = 98.00*mm;
- const G4double INNERBARREL_PCB_Width = 27.76*mm;
+ const G4double INNERBARREL_PCB_Width = 27.1*mm;
const G4double INNERBARREL_PCB_Thickness = 1.60*mm;
const G4double INNERBARREL_PCB_HoleLength = 82*mm;
const G4double INNERBARREL_PCB_WaferDepth = 1.1*mm;
const G4double INNERBARREL_PCB_Bevel1_Theta = 50*deg ;
const G4double INNERBARREL_PCB_Bevel2_Theta = 67.5*deg;
// offset between the edge of the PCB and the Edge of the hole
- const G4double INNERBARREL_PCB_Offset = 15*mm;
+ const G4double INNERBARREL_PCB_Offset = 15*mm;
const G4double INNERBARREL_ActiveWafer_Length = 94.80*mm;
- const G4double INNERBARREL_ActiveWafer_Width = 24.0*mm;
+ const G4double INNERBARREL_ActiveWafer_Width = 22.6*mm;
const G4double INNERBARREL_ActiveWafer_Thickness =400*um;
const G4double INNERBARREL_InertWafer_Length = 97.00*mm;
const G4double INNERBARREL_InertWafer_Width = 24.80*mm;
//by Shuya 180219. Equivalent to thickness in Al (note material itself is Si).
//const G4double INNERBARREL_ActiveWafer_DeadLayerThickness = 1*um;
const G4double INNERBARREL_ActiveWafer_DeadLayerThickness = 0.3*um;
- const G4double INNERBARREL_InertWafer_Thickness =
+ const G4double INNERBARREL_InertWafer_Thickness =
INNERBARREL_ActiveWafer_Thickness+ 2*INNERBARREL_ActiveWafer_DeadLayerThickness;
const G4int INNERBARREL_NumberOfStrip = 4;
-
+
// Outer Barrel //
const G4double OUTERBARREL_PCB_Length = 98.00*mm;
const G4double OUTERBARREL_PCB_Width = 33.16*mm;
@@ -104,9 +104,9 @@ namespace TIARA{
const G4double OUTERBARREL_PCB_Bevel1_Theta = 50*deg ;
const G4double OUTERBARREL_PCB_Bevel2_Theta = 67.5*deg;
// offset between the edge of the PCB and the Edge of the hole
- const G4double OUTERBARREL_PCB_Offset = 15*mm;
+ const G4double OUTERBARREL_PCB_Offset = 15*mm;
// Different from Marc code, to be checked
- const G4double OUTERBARREL_ActiveWafer_Length = 94.80*mm;
+ const G4double OUTERBARREL_ActiveWafer_Length = 96.80*mm; //94.80 previously
const G4double OUTERBARREL_ActiveWafer_Width = 29.4*mm;
const G4double OUTERBARREL_ActiveWafer_Thickness =700*um;
const G4double OUTERBARREL_InertWafer_Length = 97.00*mm;
@@ -136,7 +136,7 @@ class Tiara : public NPS::VDetector
public:
Tiara() ;
~Tiara() ;
-
+
////////////////////////////////////////////////////
//////// Specific Function of this Class ///////////
////////////////////////////////////////////////////
@@ -155,20 +155,20 @@ public:
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
void ReadConfiguration(NPL::InputParser) ;
-
+
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
void ConstructDetector(G4LogicalVolume* world) ;
-
+
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
void InitializeRootOutput() ;
-
+
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void ReadSensitive(const G4Event* event) ;
-
-
+
+
////////////////////////////////////////////////////
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
@@ -183,26 +183,26 @@ private:
private:
// Initialize material used in detector definition
void InitializeMaterial();
-
+
// List of material
G4Material* m_MaterialSilicon ;
G4Material* m_MaterialAl ;
G4Material* m_MaterialVacuum ;
G4Material* m_MaterialPCB ;
-
+
////////////////////////////////////////////////////
///////////////// Scorer Related ///////////////////
////////////////////////////////////////////////////
-
+
private:
// Initialize all Scorer
void InitializeScorers() ;
-
+
// Scorer Associate with the Silicon
G4MultiFunctionalDetector* m_InnerBarrelScorer ;
G4MultiFunctionalDetector* m_OuterBarrelScorer ;
G4MultiFunctionalDetector* m_HyballScorer ;
-
+
////////////////////////////////////////////////////
///////////////Private intern Data//////////////////
////////////////////////////////////////////////////
@@ -225,7 +225,7 @@ private:/// Visualisation Attribute:
// Light Grey
G4VisAttributes* FrameVisAtt ;
// Light Blue
- G4VisAttributes* GuardRingVisAtt ;
+ G4VisAttributes* GuardRingVisAtt ;
public:
static NPS::VDetector* Construct();
diff --git a/NPSimulation/Scorers/MDMScorer.cc b/NPSimulation/Scorers/MDMScorer.cc
index 223f1f1b91148130660568a00a4e9695e7f8c896..b86ae0d91f4efefa2d5fcdb9a6b46f939672958a 100644
--- a/NPSimulation/Scorers/MDMScorer.cc
+++ b/NPSimulation/Scorers/MDMScorer.cc
@@ -67,7 +67,10 @@ G4bool MDMScorer::ProcessHits(G4Step* aStep, G4TouchableHistory*)
G4double edep = aStep->GetTotalEnergyDeposit();
G4double M = aStep->GetPreStepPoint()->GetMass();
- G4double Q = aStep->GetPreStepPoint()->GetCharge();
+ // This gives the correct charge for a fully stripped ion (proton number)
+ G4double Q = aStep->GetTrack()->GetDynamicParticle()->GetParticleDefinition()->GetPDGCharge();
+ // THIS gives something strange (non-integers that change every step - some sort of stripping?)
+ //aStep->GetPreStepPoint()->GetCharge();
G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
G4ThreeVector MOM = aStep->GetPreStepPoint()->GetMomentumDirection();
@@ -77,7 +80,7 @@ G4bool MDMScorer::ProcessHits(G4Step* aStep, G4TouchableHistory*)
info.Charge = Q;
info.Pos = POS;
info.Mom = MOM;
-
+ // NOTE: add() calls += operator
EvtMap->add(index+DetNumber, info);
return TRUE;
}
diff --git a/Projects/T40/Analysis.cxx b/Projects/T40/Analysis.cxx
index a414d062e9e675227625f5bd9cfa9fada90d161b..6e0ec09be5eaa0180827be9d96f1517557233482 100644
--- a/Projects/T40/Analysis.cxx
+++ b/Projects/T40/Analysis.cxx
@@ -671,6 +671,31 @@ void Analysis::TreatEvent(){
// do target parameter minimization
MDM->MinimizeTarget();
+
+ //GAC 180305
+ //Calculate excitation energy from lab energy of heavy particle
+ //See http://people.physics.tamu.edu/christian/files/transfer-momentum-conservation.pdf
+ //
+ double T4 = MDM->Target_Ekin; // kinetic energy, MeV
+ double M4 = myReaction->GetNucleus4().Mass(); // rest mass, MeV/c^2
+ double P4 = sqrt(pow(T4+M4,2) - M4*M4); // momentum, MeV/c
+
+ double Theta4 = // calculate from dummy vector w/ Pz == 1
+ TVector3(tan(MDM->Target_Xang), tan(MDM->Target_Yang), 1.).Theta();
+ ThetaLab_MDM = Theta4/deg;
+ ELab_MDM = T4;
+
+ double T1 = myReaction->GetBeamEnergy();
+ double M1 = myReaction->GetNucleus1().Mass();
+ double P1 = sqrt(pow(T1+M1,2) - M1*M1);
+ double M2 = myReaction->GetNucleus2().Mass();
+ double M3 = myReaction->GetNucleus3().Mass();
+
+ double P3_2 = P4*P4+P1*P1-2*P1*P4*cos(Theta4);
+ double E3 = sqrt(M3*M3 + P3_2);
+ double T3 = E3 - M3;
+
+ Ex_MDM = sqrt(pow(T1+M1+M2-E3, 2) - P4*P4) - M4;
}
@@ -705,11 +730,17 @@ void Analysis::ReInitValue(){
ThetaLab_Barrel = -1000;
ThetaLab = -1000;
+//GAC 180305
+ ThetaLab_MDM = -1000;
+ ELab_MDM = -1000;
+
ThetaCM = -1000;
LightParticleDetected = false ;
//by Shuya 170703
Ex_Hyball = -1000 ;
Ex_Barrel = -1000 ;
+//GAC 180305
+ Ex_MDM = -1000;
//by Shuya 171019
PhiLab = -1000;
//by Shuya 171208
@@ -791,6 +822,8 @@ void Analysis::InitOutputBranch() {
//by Shuya 170703
RootOutput::getInstance()->GetTree()->Branch("Ex_Hyball",&Ex_Hyball,"Ex_Hyball/D");
RootOutput::getInstance()->GetTree()->Branch("Ex_Barrel",&Ex_Barrel,"Ex_Barrel/D");
+//GAC 180305
+ RootOutput::getInstance()->GetTree()->Branch("Ex_MDM",&Ex_MDM,"Ex_MDM/D");
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
//by Shuya 171206
@@ -801,7 +834,10 @@ void Analysis::InitOutputBranch() {
//by Shuya 171206
RootOutput::getInstance()->GetTree()->Branch("ThetaLab_Hyball",&ThetaLab_Hyball,"ThetaLab_Hyball/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab_Barrel",&ThetaLab_Barrel,"ThetaLab_Barrel/D");
-
+//GAC 180305
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab_MDM",&ThetaLab_MDM,"ThetaLab_MDM/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab_MDM",&ELab_MDM,"ELab_MDM/D");
+
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
//by Shuya 171019
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
diff --git a/Projects/T40/Analysis.h b/Projects/T40/Analysis.h
index ff48385d6f89f5a7fdfa8e5e552a0f23e66773b8..eba1fc1dfc7b0a5007fc110a7d5759aeb9583747 100644
--- a/Projects/T40/Analysis.h
+++ b/Projects/T40/Analysis.h
@@ -59,7 +59,9 @@ class Analysis: public NPL::VAnalysis{
//by Shuya 170703
double Ex_Hyball;
double Ex_Barrel;
-
+//GAC 180305
+ double Ex_MDM;
+
double ELab;
//by Shuya 171206
double ELab_Hyball;
@@ -69,7 +71,10 @@ class Analysis: public NPL::VAnalysis{
//by Shuya 171206
double ThetaLab_Hyball;
double ThetaLab_Barrel;
-
+//GAC 180305
+ double ThetaLab_MDM;
+ double ELab_MDM;
+
double ThetaCM;
double TiaraIMX;
double TiaraIMY;