diff --git a/Inputs/DetectorConfiguration/Fatima4GREAT_32c.detector b/Inputs/DetectorConfiguration/Fatima.detector
similarity index 54%
rename from Inputs/DetectorConfiguration/Fatima4GREAT_32c.detector
rename to Inputs/DetectorConfiguration/Fatima.detector
index 592e15c38a1b6591e36e2fe4551d83afc917dbaa..fa85842a1e40ae6c345be075574a03eec2aba255 100644
--- a/Inputs/DetectorConfiguration/Fatima4GREAT_32c.detector
+++ b/Inputs/DetectorConfiguration/Fatima.detector
@@ -28,13 +28,6 @@ Target
% Y= 0
% Z= 0
-
-%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
-%%Position and R given in mm
-%%Angle given in degree
-%%Option: 0,1 for Si SiLi and CsI
-%%Option: all or sensible for VISualisation
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GeneralChamber
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -63,269 +56,230 @@ Target
Fatima
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1
FatimaDetector
- THETA= 19.02390167
- PHI= 319.5114947
+ Theta= 19.02390167
+ Phi= 319.5114947
R= 182.5428674
- BETA= -3.272242219 -3.272242219 0
- VIS= all
+ Beta= -3.272242219 -3.272242219 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2
FatimaDetector
- THETA= 37.70329404
- PHI= 340.5573012
+ Theta= 37.70329404
+ Phi= 340.5573012
R= 189.784275
- BETA= -2.75768229 -2.75768229 0
- VIS= all
+ Beta= -2.75768229 -2.75768229 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3
FatimaDetector
- THETA= 57.52298638
- PHI= 346.4393719
+ Theta= 57.52298638
+ Phi= 346.4393719
R= 195.3169113
- BETA= -1.93008888 -1.93008888 0
- VIS= all
+ Beta= -1.93008888 -1.93008888 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4
FatimaDetector
- THETA= 77.29968198
- PHI= 348.489585
+ Theta= 77.29968198
+ Phi= 348.489585
R= 198.4676409
- BETA= -0.989888226 -0.989888226 0
- VIS= all
+ Beta= -0.989888226 -0.989888226 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5
FatimaDetector
- THETA= 25.83989928
- PHI= 11.03538389
+ Theta= 25.83989928
+ Phi= 11.03538389
R= 186.4716798
- BETA= -3.052133684 3.052133684 0
- VIS= all
+ Beta= -3.052133684 3.052133684 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6
FatimaDetector
- THETA= 46.55103678
- PHI= 6.372873059
+ Theta= 46.55103678
+ Phi= 6.372873059
R= 193.0582918
- BETA= -2.336257006 2.336257006 0
- VIS= all
+ Beta= -2.336257006 2.336257006 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7
FatimaDetector
- THETA= 67.02589098
- PHI= 4.906399705
+ Theta= 67.02589098
+ Phi= 4.906399705
R= 197.5675029
- BETA= -1.342490056 1.342490056 0
- VIS= all
+ Beta= -1.342490056 1.342490056 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8
FatimaDetector
- THETA= 107.1442293
- PHI= 4.707364909
+ Theta= 107.1442293
+ Phi= 4.707364909
R= 198.3842185
- BETA= -1.028278389 1.028278389 0
- VIS= all
+ Beta= -1.028278389 1.028278389 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9
FatimaDetector
- THETA= 26.85689097
- PHI= 57.83150522
+ Theta= 26.85689097
+ Phi= 57.83150522
R= 182.3833519
- BETA= -3.278549609 3.278549609 0
- VIS= all
+ Beta= -3.278549609 3.278549609 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10
FatimaDetector
- THETA= 41.56673389
- PHI= 33.74691678
+ Theta= 41.56673389
+ Phi= 33.74691678
R= 189.2241657
- BETA= -2.815552839 2.815552839 0
- VIS= all
+ Beta= -2.815552839 2.815552839 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11
FatimaDetector
- THETA= 59.43554149
- PHI= 24.61616818
+ Theta= 59.43554149
+ Phi= 24.61616818
R= 194.4599021
- BETA= -2.099135856 2.099135856 0
- VIS= all
+ Beta= -2.099135856 2.099135856 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 12
FatimaDetector
- THETA= 77.97363101
- PHI= 21.17191049
+ Theta= 77.97363101
+ Phi= 21.17191049
R= 197.4448302
- BETA= -1.382760452 1.382760452 0
- VIS= all
+ Beta= -1.382760452 1.382760452 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13
FatimaDetector
- THETA= 115.2395483
- PHI= 23.21965799
+ Theta= 115.2395483
+ Phi= 23.21965799
R= 195.5744725
- BETA= -1.874825744 1.874825744 0
- VIS= all
+ Beta= -1.874825744 1.874825744 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 14
FatimaDetector
- THETA= 46.94600569
- PHI= 61.13966047
+ Theta= 46.94600569
+ Phi= 61.13966047
R= 184.8217039
- BETA= -3.160159266 3.160159266 0
- VIS= all
+ Beta= -3.160159266 3.160159266 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15
FatimaDetector
- THETA= 72.54645084
- PHI= 39.89896247
+ Theta= 72.54645084
+ Phi= 39.89896247
R= 193.2944032
- BETA= -2.299216098 2.299216098 0
- VIS= all
+ Beta= -2.299216098 2.299216098 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 16
FatimaDetector
- THETA= 103.1039385
- PHI= 38.77251614
+ Theta= 103.1039385
+ Phi= 38.77251614
R= 193.922275
- BETA= -2.195219842 2.195219842 0
- VIS= all
+ Beta= -2.195219842 2.195219842 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 17
%FatimaDetector
-% THETA= 90
-% PHI= 180
+% Theta= 90
+% Phi= 180
% R= 21
-% BETA= #DIV/0! #DIV/0! 0
-% VIS= all
+% Beta= #DIV/0! #DIV/0! 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 18
%FatimaDetector
-% THETA= 90
-% PHI= 180
+% Theta= 90
+% Phi= 180
% R= 21
-% BETA= #DIV/0! #DIV/0! 0
-% VIS= all
+% Beta= #DIV/0! #DIV/0! 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 19
%FatimaDetector
-% THETA= 90
-% PHI= 180
+% Theta= 90
+% Phi= 180
% R= 21
-% BETA= #DIV/0! #DIV/0! 0
-% VIS= all
+% Beta= #DIV/0! #DIV/0! 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 20
FatimaDetector
- THETA= 19.02390167
- PHI= 220.4885053
+ Theta= 19.02390167
+ Phi= 220.4885053
R= 182.5428674
- BETA= -3.272242219 -3.272242219 0
- VIS= all
+ Beta= -3.272242219 -3.272242219 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 21
FatimaDetector
- THETA= 37.70329404
- PHI= 199.4426988
+ Theta= 37.70329404
+ Phi= 199.4426988
R= 189.784275
- BETA= -2.75768229 -2.75768229 0
- VIS= all
+ Beta= -2.75768229 -2.75768229 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 22
FatimaDetector
- THETA= 57.52298638
- PHI= 193.5606281
+ Theta= 57.52298638
+ Phi= 193.5606281
R= 195.3169113
- BETA= -1.93008888 -1.93008888 0
- VIS= all
+ Beta= -1.93008888 -1.93008888 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 23
FatimaDetector
- THETA= 77.29968198
- PHI= 191.510415
+ Theta= 77.29968198
+ Phi= 191.510415
R= 198.4676409
- BETA= -0.989888226 -0.989888226 0
- VIS= all
+ Beta= -0.989888226 -0.989888226 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 24
FatimaDetector
- THETA= 25.83989928
- PHI= 168.9646161
+ Theta= 25.83989928
+ Phi= 168.9646161
R= 186.4716798
- BETA= -3.052133684 -3.052133684 0
- VIS= all
+ Beta= -3.052133684 -3.052133684 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 25
FatimaDetector
- THETA= 46.55103678
- PHI= 173.6271269
+ Theta= 46.55103678
+ Phi= 173.6271269
R= 193.0582918
- BETA= -2.336257006 -2.336257006 0
- VIS= all
+ Beta= -2.336257006 -2.336257006 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 26
FatimaDetector
- THETA= 67.02589098
- PHI= 175.0936003
+ Theta= 67.02589098
+ Phi= 175.0936003
R= 197.5675029
- BETA= -1.342490056 -1.342490056 0
- VIS= all
+ Beta= -1.342490056 -1.342490056 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 27
FatimaDetector
- THETA= 107.1442293
- PHI= 175.2926351
+ Theta= 107.1442293
+ Phi= 175.2926351
R= 198.3842185
- BETA= -1.028278389 -1.028278389 0
- VIS= all
+ Beta= -1.028278389 -1.028278389 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 28
FatimaDetector
- THETA= 26.85689097
- PHI= 122.1684948
+ Theta= 26.85689097
+ Phi= 122.1684948
R= 182.3833519
- BETA= -3.278549609 -3.278549609 0
- VIS= all
+ Beta= -3.278549609 -3.278549609 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 29
FatimaDetector
- THETA= 41.56673389
- PHI= 146.2530832
+ Theta= 41.56673389
+ Phi= 146.2530832
R= 189.2241657
- BETA= -2.815552839 -2.815552839 0
- VIS= all
+ Beta= -2.815552839 -2.815552839 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 30
FatimaDetector
- THETA= 59.43554149
- PHI= 155.3838318
+ Theta= 59.43554149
+ Phi= 155.3838318
R= 194.4599021
- BETA= -2.099135856 -2.099135856 0
- VIS= all
+ Beta= -2.099135856 -2.099135856 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 31
FatimaDetector
- THETA= 77.97363101
- PHI= 158.8280895
+ Theta= 77.97363101
+ Phi= 158.8280895
R= 197.4448302
- BETA= -1.382760452 -1.382760452 0
- VIS= all
+ Beta= -1.382760452 -1.382760452 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 32
FatimaDetector
- THETA= 115.2395483
- PHI= 156.780342
+ Theta= 115.2395483
+ Phi= 156.780342
R= 195.5744725
- BETA= -1.874825744 -1.874825744 0
- VIS= all
+ Beta= -1.874825744 -1.874825744 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 33
FatimaDetector
- THETA= 46.94600569
- PHI= 118.8603395
+ Theta= 46.94600569
+ Phi= 118.8603395
R= 184.8217039
- BETA= -3.160159266 -3.160159266 0
- VIS= all
+ Beta= -3.160159266 -3.160159266 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 34
FatimaDetector
- THETA= 72.54645084
- PHI= 140.1010375
+ Theta= 72.54645084
+ Phi= 140.1010375
R= 193.2944032
- BETA= -2.299216098 -2.299216098 0
- VIS= all
+ Beta= -2.299216098 -2.299216098 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 35
FatimaDetector
- THETA= 103.1039385
- PHI= 141.2274839
+ Theta= 103.1039385
+ Phi= 141.2274839
R= 193.922275
- BETA= -2.195219842 -2.195219842 0
- VIS= all
+ Beta= -2.195219842 -2.195219842 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36
%FatimaDetector
-% THETA= 6.709836808
-% PHI= 180
+% Theta= 6.709836808
+% Phi= 180
% R= 179.7310491
-% BETA= -3.354918404 -3.354918404 0
-% VIS= all
+% Beta= -3.354918404 -3.354918404 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37
%FatimaDetector
-% THETA= 6.709836808
-% PHI= 180
+% Theta= 6.709836808
+% Phi= 180
% R= 179.7310491
-% BETA= -3.354918404 -3.354918404 0
-% VIS= all
+% Beta= -3.354918404 -3.354918404 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 38
%FatimaDetector
-% THETA= 6.709836808
-% PHI= 180
+% Theta= 6.709836808
+% Phi= 180
% R= 179.7310491
-% BETA= -3.354918404 -3.354918404 0
-% FIRSTSTAGE= 1
-% VIS= 0
+% Beta= -3.354918404 -3.354918404 0
diff --git a/NPLib/Fatima/Fatima.cxx b/NPLib/Fatima/Fatima.cxx
deleted file mode 100644
index 9ca03bc56eaa68b724895ab1e761353561638818..0000000000000000000000000000000000000000
--- a/NPLib/Fatima/Fatima.cxx
+++ /dev/null
@@ -1,822 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/2014 *
- *---------------------------------------------------------------------------*
- * Decription: This class is mainly an interface to the *
- * TFatimaPhysics class *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#include "Fatima.h"
-
-// C++ headers
-#include <iostream>
-#include <fstream>
-#include <string>
-#include <cmath>
-#include <stdlib.h>
-
-// NPL headers
-#include "RootInput.h"
-#include "RootOutput.h"
-
-// ROOT headers
-#include "TChain.h"
-
-using namespace std ;
-
-// Default Constructor
-
-Fatima::Fatima()
-{
-
- m_NumberOfDetector = 0; //
- m_EventData = new TFatimaData();
- m_EventPhysics = new TFatimaPhysics();
-}
-
-
-
-Fatima::~Fatima()
-{
-
- m_NumberOfDetector = 0;
- delete m_EventData;
- delete m_EventPhysics;
-}
-
-
-
-// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
-void Fatima::ReadConfiguration(string Path)
-{
- ifstream ConfigFile ;
- ConfigFile.open(Path.c_str()) ;
- string LineBuffer ;
- string DataBuffer ;
-
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
- double Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Az, Bz, Cz, Dz;
- TVector3 A, B, C, D;
- double Theta = 0, Phi = 0, R = 0, beta_u = 0 , beta_v = 0 , beta_w = 0;
-
- bool check_A = false;
- bool check_C = false;
- bool check_B = false;
- bool check_D = false;
-
- bool check_Theta = false;
- bool check_Phi = false;
- bool check_R = false;
- bool check_beta = false;
-
- bool ReadingStatus = false;
-
- bool isDetector = false;
-
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
-
- // If line is a Fatima bloc, reading toggle to true
- // and toggle to true flags indicating which shape is treated.
- if (LineBuffer.compare(0, 14, "FatimaDetector") == 0 ) {
- cout << "///////////////////////" << endl;
- cout << "Module found:" << endl;
-
- if (LineBuffer.compare(0, 14, "FatimaDetector") == 0) isDetector = true;
- ReadingStatus = true;
- }
- // Else don't toggle to Reading Block Status
- else ReadingStatus = false;
-
- // Reading Block
- while (ReadingStatus) {
- if (isDetector) { // FatimaDetector shape
- ConfigFile >> DataBuffer ;
- // Comment Line
- if (DataBuffer.compare(0, 1, "%") == 0) {
- ConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
- }
- // Finding another telescope (safety), toggle out
- else if (DataBuffer.compare(0, 14, "FatimaDetector") == 0) {
- cout << "WARNING: Another Module is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
- ReadingStatus = false;
- }
-
- // Position method
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- Ax = atof(DataBuffer.c_str()) ;
- Ax = Ax ;
- ConfigFile >> DataBuffer ;
- Ay = atof(DataBuffer.c_str()) ;
- Ay = Ay ;
- ConfigFile >> DataBuffer ;
- Az = atof(DataBuffer.c_str()) ;
- Az = Az ;
-
- A = TVector3(Ax, Ay, Az);
- cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
- }
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- Bx = atof(DataBuffer.c_str()) ;
- Bx = Bx ;
- ConfigFile >> DataBuffer ;
- By = atof(DataBuffer.c_str()) ;
- By = By ;
- ConfigFile >> DataBuffer ;
- Bz = atof(DataBuffer.c_str()) ;
- Bz = Bz ;
-
- B = TVector3(Bx, By, Bz);
- cout << "X128 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
- }
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- Cx = atof(DataBuffer.c_str()) ;
- Cx = Cx ;
- ConfigFile >> DataBuffer ;
- Cy = atof(DataBuffer.c_str()) ;
- Cy = Cy ;
- ConfigFile >> DataBuffer ;
- Cz = atof(DataBuffer.c_str()) ;
- Cz = Cz ;
-
- C = TVector3(Cx, Cy, Cz);
- cout << "X1 Y128 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
- }
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- Dx = atof(DataBuffer.c_str()) ;
- Dx = Dx ;
- ConfigFile >> DataBuffer ;
- Dy = atof(DataBuffer.c_str()) ;
- Dy = Dy ;
- ConfigFile >> DataBuffer ;
- Dz = atof(DataBuffer.c_str()) ;
- Dz = Dz ;
-
- D = TVector3(Dx, Dy, Dz);
- cout << "X128 Y128 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
- } // End Position Method
-
- // Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- Theta = Theta ;
- cout << "Theta: " << Theta << endl;
- }
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- Phi = Phi ;
- cout << "Phi: " << Phi << endl;
- }
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
- R = R ;
- cout << "R: " << R << endl;
- }
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
- check_beta = true;
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- beta_u = beta_u ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- beta_v = beta_v ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- beta_w = beta_w ;
- cout << "Beta: " << beta_u << " " << beta_v << " " << beta_w << endl ;
- }
-
- /////////////////////////////////////////////////
- // If All necessary information there, toggle out
- if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) ) {
- ReadingStatus = false;
-
- // Add The previously define telescope
- // With position method
- if ( check_A && check_B && check_C && check_D ) {
- AddModuleDummyShape(A ,
- B ,
- C ,
- D ) ;
- }
-
- // with angle method
- else if ( check_Theta && check_Phi && check_R && check_beta ) {
- AddModuleDummyShape(Theta,
- Phi,
- R,
- beta_u,
- beta_v,
- beta_w);
- }
-
- // reset boolean flag for point positioning
- check_A = false;
- check_B = false;
- check_C = false;
- check_D = false;
-
- // reset boolean flag for angle positioning
- check_Theta = false;
- check_Phi = false;
- check_R = false;
- check_beta = false;
-
- // reset boolean flag for shape determination
- isDetector = false;
-
- } // end test for adding a module
- } // end test for FatimaDetector shape
-
-
- } // end while for reading block
- } // end while for reading file
-
- cout << endl << "/////////////////////////////" << endl<<endl;
-}
-
-
-// Read stream at Path and pick-up calibration parameter using Token
-// If argument is "Simulation" no change calibration is loaded
-void Fatima::ReadCalibrationFile(string Path)
-{
- // Order of Polynom function used for calibration
- int Calibration_Si_E_Order;
- int Calibration_Si_T_Order;
- int Calibration_SiLi_E_Order;
- int Calibration_CsI_E_Order;
-
- // Calibration_Si_X_E[DetectorNumber][StripNumber][Order of Coeff]
- vector< vector< vector< double > > > Calibration_Si_X_E ;
- vector< vector< vector< double > > > Calibration_Si_X_T ;
- vector< vector< vector< double > > > Calibration_Si_Y_E ;
- vector< vector< vector< double > > > Calibration_Si_Y_T ;
-
- // Calibration_SiLi_E[DetectorNumber][PadNumber][Order of Coeff]
- vector< vector< vector< double > > > Calibration_SiLi_E ;
-
- // Calibration_SiLi_E[DetectorNumber][CrystalNumber][Order of Coeff]
- vector< vector< vector< double > > > Calibration_CsI_E ;
-
- if (Path == "Simulation") { // Simulation case: data already calibrated
- Calibration_Si_E_Order = 1;
- Calibration_Si_T_Order = 1;
- Calibration_SiLi_E_Order = 1;
- Calibration_CsI_E_Order = 1;
-
- vector<double> Coef;
- // Order 0 Order 1
- Coef.push_back(0) ; Coef.push_back(1) ;
-
- vector< vector<double> > StripLine ;
- StripLine.resize( 128 , Coef) ;
-
- Calibration_Si_X_E.resize( m_NumberOfDetector , StripLine) ;
- Calibration_Si_X_T.resize( m_NumberOfDetector , StripLine) ;
- Calibration_Si_Y_E.resize( m_NumberOfDetector , StripLine) ;
- Calibration_Si_Y_T.resize( m_NumberOfDetector , StripLine) ;
-
- Calibration_SiLi_E.resize( m_NumberOfDetector , StripLine) ;
- Calibration_CsI_E .resize( m_NumberOfDetector , StripLine) ;
- }
- else {
- }
-}
-
-
-
-// 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
-void Fatima::InitializeRootInput()
-{
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus("FATIMA", true);
- inputChain->SetBranchStatus("fFATIMA*", true);
- inputChain->SetBranchAddress("FATIMA", &m_EventData);
-}
-
-
-
-// Create associated branches and associated private member DetectorPhysics address
-void Fatima::InitializeRootOutput()
-{
- TTree* outputTree = RootOutput::getInstance()->GetTree();
- outputTree->Branch("FATIMA", "TFatimaPhysics", &m_EventPhysics);
-}
-
-
-
-// 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 Fatima::BuildPhysicalEvent()
-{
- m_EventPhysics -> BuildPhysicalEvent(m_EventData);
-}
-
-
-
-// Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
-// This method aimed to be used for analysis performed during experiment, when speed is requiered.
-// NB: This method can eventually be the same as BuildPhysicalEvent.
-void Fatima::BuildSimplePhysicalEvent()
-{
- m_EventPhysics -> BuildSimplePhysicalEvent(m_EventData);
-}
-
-
-//
-// Used for 3x3 clusters (FatimaCluster):
-//
-void Fatima::AddModuleSquare(TVector3 C_X1_Y1,
- TVector3 C_X128_Y1,
- TVector3 C_X1_Y128,
- TVector3 C_X128_Y128)
-{
- m_NumberOfDetector++;
-
- cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
-
- // remove warning using C_X128_Y128
- C_X128_Y128.Unit();
-
- // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U = C_X128_Y1 - C_X1_Y1;
- U = U.Unit();
-
- // Vector V on Module Face (parallele to X Strip)
- TVector3 V = C_X1_Y128 - C_X1_Y1;
- V = V.Unit();
-
- // Position Vector of Strip Center
- TVector3 StripCenter = TVector3(0,0,0);
- // Position Vector of X=1 Y=1 Strip
- TVector3 Strip_1_1;
-
-
- // !!!!!! for cluster 3x3 only !!!!!
- //
- // Geometry Parameter in case of 3x3 cluster
- double Face = 169; // mm cf: FatimaCluster.hh
- double NumberOfStrip = 3; // number of crystals per raw or per column
- double StripPitch = Face/NumberOfStrip; // mm
- //
- //double Face = 98; // mm
- //double NumberOfStrip = 128;
- //double StripPitch = Face/NumberOfStrip; // mm
- //
- //
-
- // Buffer object to fill Position Array
- vector<double> lineX;
- vector<double> lineY;
- vector<double> lineZ;
-
- vector< vector< double > > OneModuleStripPositionX;
- vector< vector< double > > OneModuleStripPositionY;
- vector< vector< double > > OneModuleStripPositionZ;
-
- // Moving StripCenter to 1.1 corner:
- Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.);
-
- for (int i = 0; i < NumberOfStrip; i++) {
- lineX.clear();
- lineY.clear();
- lineZ.clear();
-
- for (int j = 0; j < NumberOfStrip; j++) {
- StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
-// StripCenter += -TargetPosition;
-
- lineX.push_back( StripCenter.X() );
- lineY.push_back( StripCenter.Y() );
- lineZ.push_back( StripCenter.Z() );
- }
-
- OneModuleStripPositionX.push_back(lineX);
- OneModuleStripPositionY.push_back(lineY);
- OneModuleStripPositionZ.push_back(lineZ);
- }
-
- m_DetPositionX.push_back( OneModuleStripPositionX );
- m_DetPositionY.push_back( OneModuleStripPositionY );
- m_DetPositionZ.push_back( OneModuleStripPositionZ );
-}
-
-
-
-void Fatima::AddModuleSquare(double theta,
- double phi,
- double distance,
- double beta_u,
- double beta_v,
- double beta_w)
-{
- m_NumberOfDetector++;
-
-
- // convert from degree to radian:
- double Pi = 3.141592654;
- theta = theta * Pi/180. ;
- phi = phi * Pi/180. ;
-
- // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U ;
- // Vector V on Module Face (parallele to X Strip)
- TVector3 V ;
- // Vector W normal to Module Face (pointing CsI)
- TVector3 W ;
- // Vector position of Module Face center
- TVector3 C ;
-
- C = TVector3(distance * sin(theta) * cos(phi),
- distance * sin(theta) * sin(phi),
- distance * cos(theta));
-
- TVector3 YperpC = TVector3( cos(theta) * cos(phi),
- cos(theta) * sin(phi),
- -sin(theta));
-
- W = C.Unit();
- U = W.Cross(YperpC);
- V = W.Cross(U);
-
- U = U.Unit();
- V = V.Unit();
-
- U.Rotate( beta_u * Pi/180. , U ) ;
- V.Rotate( beta_u * Pi/180. , U ) ;
-
- U.Rotate( beta_v * Pi/180. , V ) ;
- V.Rotate( beta_v * Pi/180. , V ) ;
-
- U.Rotate( beta_w * Pi/180. , W ) ;
- V.Rotate( beta_w * Pi/180. , W ) ;
-
-
- //
- // !!!!!! for cluster 3x3 only !!!!!
- //
- // Geometry Parameter in case of 3x3 cluster
- double Face = 169; // mm cf: FatimaCluster.hh
- double NumberOfStrip = 3; // number of crystals per raw or per column
- double StripPitch = Face/NumberOfStrip; // mm
- //
- //double Face = 98; // mm
- //double NumberOfStrip = 128;
- //double StripPitch = Face/NumberOfStrip; // mm
- //
- //
-
- vector<double> lineX;
- vector<double> lineY;
- vector<double> lineZ;
-
- vector< vector< double > > OneModuleStripPositionX;
- vector< vector< double > > OneModuleStripPositionY;
- vector< vector< double > > OneModuleStripPositionZ;
-
- double X, Y, Z;
-
- // Moving C to the 1.1 corner:
- C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
- C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
- C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
-
- for (int i = 0; i < NumberOfStrip; i++) {
- lineX.clear();
- lineY.clear();
- lineZ.clear();
-
- for (int j = 0; j < NumberOfStrip; j++) {
- X = C.X() + StripPitch * ( U.X()*i + V.X()*j );
- Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j );
- Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j );
-
- lineX.push_back(X);
- lineY.push_back(Y);
- lineZ.push_back(Z);
- }
-
- OneModuleStripPositionX.push_back(lineX);
- OneModuleStripPositionY.push_back(lineY);
- OneModuleStripPositionZ.push_back(lineZ);
- }
-
- m_DetPositionX.push_back( OneModuleStripPositionX );
- m_DetPositionY.push_back( OneModuleStripPositionY );
- m_DetPositionZ.push_back( OneModuleStripPositionZ );
-
-
-}
-
-
-//
-// Used for single detector (FatimaDetector):
-//
-void Fatima::AddModuleDummyShape(TVector3 C_X1_Y1,
- TVector3 C_X128_Y1,
- TVector3 C_X1_Y128,
- TVector3 C_X128_Y128)
-{
- m_NumberOfDetector++;
-
- cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
-
- // remove warning using C_X128_Y128
- C_X128_Y128.Unit();
-
- // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U = C_X128_Y1 - C_X1_Y1;
- U = U.Unit();
-
- // Vector V on Module Face (parallele to X Strip)
- TVector3 V = C_X1_Y128 - C_X1_Y1;
- V = V.Unit();
-
- // Position Vector of Strip Center
- TVector3 StripCenter = TVector3(0,0,0);
- // Position Vector of X=1 Y=1 Strip
- TVector3 Strip_1_1;
-
- //
- // !!!!!! for for single phoswich !!!!
- // Geometry Parameter in case of 3x3 cluster
- double Face = 57; // mm cf: FatimaCluster.hh
- double NumberOfStrip = 1; // number of crystals per raw or per column
- double StripPitch = Face/NumberOfStrip; // mm
- //
- //
-
- // Buffer object to fill Position Array
- vector<double> lineX;
- vector<double> lineY;
- vector<double> lineZ;
-
- vector< vector< double > > OneModuleStripPositionX;
- vector< vector< double > > OneModuleStripPositionY;
- vector< vector< double > > OneModuleStripPositionZ;
-
- // Moving StripCenter to 1.1 corner:
- Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.);
-
- for (int i = 0; i < NumberOfStrip; i++) {
- lineX.clear();
- lineY.clear();
- lineZ.clear();
-
- for (int j = 0; j < NumberOfStrip; j++) {
- StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
-// StripCenter += -TargetPosition;
-
- lineX.push_back( StripCenter.X() );
- lineY.push_back( StripCenter.Y() );
- lineZ.push_back( StripCenter.Z() );
- }
-
- OneModuleStripPositionX.push_back(lineX);
- OneModuleStripPositionY.push_back(lineY);
- OneModuleStripPositionZ.push_back(lineZ);
- }
-
- m_DetPositionX.push_back( OneModuleStripPositionX );
- m_DetPositionY.push_back( OneModuleStripPositionY );
- m_DetPositionZ.push_back( OneModuleStripPositionZ );
-}
-
-
-
-
-void Fatima::AddModuleDummyShape(double theta,
- double phi,
- double distance,
- double beta_u,
- double beta_v,
- double beta_w)
-{
- m_NumberOfDetector++;
-
- cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
-
- // convert from degree to radian:
- double Pi = 3.141592654;
- theta = theta * Pi/180. ;
- phi = phi * Pi/180. ;
-
- // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U ;
- // Vector V on Module Face (parallele to X Strip)
- TVector3 V ;
- // Vector W normal to Module Face (pointing CsI)
- TVector3 W ;
- // Vector position of Module Face center
- TVector3 C ;
-
- C = TVector3(distance * sin(theta) * cos(phi),
- distance * sin(theta) * sin(phi),
- distance * cos(theta));
-
- TVector3 YperpW = TVector3( cos(theta) * cos(phi),
- cos(theta) * sin(phi),
- -sin(theta));
-
- W = C.Unit();
- U = W.Cross(YperpW);
- V = W.Cross(U);
-
- U = U.Unit();
- V = V.Unit();
-
- U.Rotate( beta_u * Pi/180. , U ) ;
- V.Rotate( beta_u * Pi/180. , U ) ;
-
- U.Rotate( beta_v * Pi/180. , V ) ;
- V.Rotate( beta_v * Pi/180. , V ) ;
-
- U.Rotate( beta_w * Pi/180. , W ) ;
- V.Rotate( beta_w * Pi/180. , W ) ;
-
- double Face = 50; // mm
- double NumberOfStrip = 25;
- double StripPitch = Face/NumberOfStrip; // mm
-
- vector<double> lineX;
- vector<double> lineY;
- vector<double> lineZ;
-
- vector< vector< double > > OneModuleStripPositionX;
- vector< vector< double > > OneModuleStripPositionY;
- vector< vector< double > > OneModuleStripPositionZ;
-
- double X, Y, Z;
-
- // Moving C to the 1.1 corner:
- C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
- C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
- C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
-
- for (int i = 0; i < NumberOfStrip; i++) {
- lineX.clear();
- lineY.clear();
- lineZ.clear();
-
- for (int j = 0; j < NumberOfStrip; j++) {
- X = C.X() + StripPitch * ( U.X()*i + V.X()*j );
- Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j );
- Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j );
-
- lineX.push_back(X);
- lineY.push_back(Y);
- lineZ.push_back(Z);
- }
-
- OneModuleStripPositionX.push_back(lineX);
- OneModuleStripPositionY.push_back(lineY);
- OneModuleStripPositionZ.push_back(lineZ);
- }
-
- m_DetPositionX.push_back( OneModuleStripPositionX );
- m_DetPositionY.push_back( OneModuleStripPositionY );
- m_DetPositionZ.push_back( OneModuleStripPositionZ );
-}
-
-
-
-double Fatima::GetEnergyDeposit()
-{
- if (m_EventPhysics->FatimaTotalEnergy.size() > 0)
- return m_EventPhysics->FatimaTotalEnergy[0];
- else
- return -1000;
-}
-
-double Fatima::GetEnergyInDeposit() // inner Layer
-{
- if (m_EventPhysics->FatimaInTotalEnergy.size() > 0)
- return m_EventPhysics->FatimaInTotalEnergy[0];
- else
- return -1000;
-}
-
-double Fatima::GetEnergyOutDeposit() // Outer Layer
-{
- if (m_EventPhysics->FatimaOutTotalEnergy.size() > 0)
- return m_EventPhysics->FatimaOutTotalEnergy[0];
- else
- return -1000;
-}
-
-
-TVector3 Fatima::GetPositionOfInteraction()
-{
- TVector3 Position = TVector3(-1000,-1000,-1000);
-
-
- if(m_EventPhysics->DetectorNumber[0]>100 ){ // ie: cluster
-
- cout << "Cluster#: " << m_EventPhysics->DetectorNumber[0] << endl;
- cout << "Crystal column (x) #: " << m_EventPhysics->FirstStage_X[0] << endl;
- cout << "Crystal raw (y) #: " << m_EventPhysics->FirstStage_Y[0] << endl;
-
- if(m_EventPhysics->FirstStage_X[0]>-1 && m_EventPhysics->FirstStage_Y[0]>-1){ // ie: LaBr hit first
- Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
- GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
- GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]));
- }/*
- else
- {
- Position = TVector3(GetCrystPositionX(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
- GetCrystPositionY(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
- GetCrystPositionZ(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]));
-
- }
- */
- }
-
- if(m_EventPhysics->DetectorNumber[0]>0 ){ // ie: phoswich
-
- cout << "Detector #: " << m_EventPhysics->DetectorNumber[0] << endl;
- cout << "Crystal column (x) #: " << m_EventPhysics->FirstStage_X[0] << endl;
- cout << "Crystal raw (y) #: " << m_EventPhysics->FirstStage_Y[0] << endl;
-
-
- if(m_EventPhysics->FirstStage_X[0]>-1 && m_EventPhysics->FirstStage_Y[0]>-1){ // ie: LaBr hit first
- Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
- GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
- GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]));
- }/*else
- {
- Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
- GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
- GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]));
-
- }*/
-
- }
-
-
- return(Position);
-}
-
-
-
-void Fatima::Print()
-{
- cout << "Number of Detectors: " << m_NumberOfDetector << endl;
-
- for (int f = 0; f < m_NumberOfDetector; f++) {
- cout << "Detector " << f+1 << endl;
-
- for (int i = 0; i < 3; i++) {
- for (int j = 0; j < 3; j++) {
- cout << i+1 << " "<< j+1 << " "
- << m_DetPositionX[f][i][j] << " "
- << m_DetPositionY[f][i][j] << " "
- << m_DetPositionZ[f][i][j] << " "
- << endl ;
- }
- }
-
- /*
- for (int i = 0; i < 128; i++) {
- for (int j = 0; j < 128; j++) {
- cout << i+1 << " "<< j+1 << " "
- << m_StripPositionX[f][i][j] << " "
- << m_StripPositionY[f][i][j] << " "
- << m_StripPositionZ[f][i][j] << " "
- << endl ;
- }
- }
- */
-
- }
-}
diff --git a/NPLib/Fatima/Fatima.h b/NPLib/Fatima/Fatima.h
deleted file mode 100644
index 903740901643e5587772c81f0f6a15096ad96e0f..0000000000000000000000000000000000000000
--- a/NPLib/Fatima/Fatima.h
+++ /dev/null
@@ -1,166 +0,0 @@
-
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/2014 *
- *---------------------------------------------------------------------------*
- * Decription: This class is mainly an interface to the *
- * TFatimaPhysics class *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#ifndef Fatima_H
-#define Fatima_H
-
-// NPL
-#include "../include/VDetector.h"
-#include "TFatimaData.h"
-#include "TFatimaPhysics.h"
-
-// Root
-#include "TVector3.h"
-
-class Fatima : public NPA::VDetector
-{
-public:
- Fatima();
- virtual ~Fatima();
-
-public:
- /////////////////////////////////////
- // Innherited from VDetector Class //
- /////////////////////////////////////
- // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
- void ReadConfiguration(string);
-
- // Read stream at CalibFile and pick-up calibration parameter using Token
- // If argument is "Simulation" no change calibration is loaded
- void ReadCalibrationFile(string);
-
- // 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
- void InitializeRootInput();
-
- // Create associated branches and associated private member DetectorPhysics address
- void InitializeRootOutput();
-
- // This method is called at each event read from the Input Tree.
- // The aim 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 ...).
- // This method aimed to be used for analysis performed during experiment, when speed is requiered.
- // NB: This method can eventually be the same as BuildPhysicalEvent.
- void BuildSimplePhysicalEvent();
-
- // Those two method all to clear the Event Physics or Data
- void ClearEventPhysics() {m_EventPhysics->Clear();}
- void ClearEventData() {m_EventData->Clear();}
-
-
-public:
- ////////////////////////////////
- // Specific to Fatima //
- ////////////////////////////////
- // Case of a Square module
- // Add a Module using Corner Coordinate information
- void AddModuleSquare(TVector3 C_X1_Y1,
- TVector3 C_X128_Y1,
- TVector3 C_X1_Y128,
- TVector3 C_X128_Y128);
-
- // Add a Module using R Theta Phi of Si center information
- void AddModuleSquare(double theta,
- double phi,
- double distance,
- double beta_u,
- double beta_v,
- double beta_w);
-
- // Case of a DummyShape module
- // Add a Module using Corner Coordinate information
- void AddModuleDummyShape(TVector3 C_X1_Y1,
- TVector3 C_X128_Y1,
- TVector3 C_X1_Y128,
- TVector3 C_X128_Y128);
-
- // Add a Module using R Theta Phi of Si center information
- void AddModuleDummyShape(double theta,
- double phi,
- double distance,
- double beta_u,
- double beta_v,
- double beta_w);
-
- // Getters to retrieve the (X,Y,Z) coordinates of a pixel (crystal) defined by (X,Y) for Detector
-
- double GetDetPositionX(int N ,int X ,int Y) { return m_DetPositionX[N-1][X][Y]; } // remember index=0 for Fatima Detector (cf FatimaModule.cxx)
- double GetDetPositionY(int N ,int X ,int Y) { return m_DetPositionY[N-1][X][Y]; }
- double GetDetPositionZ(int N ,int X ,int Y) { return m_DetPositionZ[N-1][X][Y]; }
-
- double GetNumberOfDetector() { return m_NumberOfDetector; }
-
- /*
- // Getters to retrieve the (X,Y,Z) coordinates of a pixel (crystal) defined by (X,Y) for clusters
-
- double GetCrystPositionX(int N ,int X ,int Y) { return m_CrystPositionX[N-101][X][Y]; } // remember index=100 for Fatima Cluster (cf FatimaModule.cxx)
- double GetCrystPositionY(int N ,int X ,int Y) { return m_CrystPositionY[N-101][X][Y]; }
- double GetCrystPositionZ(int N ,int X ,int Y) { return m_CrystPositionZ[N-101][X][Y]; }
- */
-
-
- //double GetStripPositionX(int N ,int X ,int Y) { return m_StripPositionX[N-1][X-1][Y-1]; }
- //double GetStripPositionY(int N ,int X ,int Y) { return m_StripPositionY[N-1][X-1][Y-1]; }
- //double GetStripPositionZ(int N ,int X ,int Y) { return m_StripPositionZ[N-1][X-1][Y-1]; }
-
- //double GetNumberOfModule() { return m_NumberOfModule; }
-
-
- // Get Root input and output objects
- TFatimaData* GetEventData() {return m_EventData;}
- TFatimaPhysics* GetEventPhysics() {return m_EventPhysics;}
-
- // To be called after a build Physical Event
- double GetEnergyDeposit();
- double GetEnergyInDeposit();
- double GetEnergyOutDeposit();
- TVector3 GetPositionOfInteraction();
-
- void Print();
-
-
-private:
- ////////////////////////////////////////
- // Root Input and Output tree classes //
- ////////////////////////////////////////
- TFatimaData* m_EventData;
- TFatimaPhysics* m_EventPhysics;
-
-
-private:
- // Spatial Position of Crystal Calculated on basis of detector(cluster) position
- /*
- int m_NumberOfModule; // cluster
- vector< vector < vector < double > > > m_CrystPositionX;
- vector< vector < vector < double > > > m_CrystPositionY;
- vector< vector < vector < double > > > m_CrystPositionZ;
- */
- int m_NumberOfDetector; // single detector
- vector< vector < vector < double > > > m_DetPositionX;
- vector< vector < vector < double > > > m_DetPositionY;
- vector< vector < vector < double > > > m_DetPositionZ;
-
-};
-
-#endif
diff --git a/NPLib/Fatima/Makefile b/NPLib/Fatima/Makefile
index e25e9da81bcd88dfec2c030dca76eaaf5da7cac3..96f038876be8c9c1798d076d04ebc282decc128d 100644
--- a/NPLib/Fatima/Makefile
+++ b/NPLib/Fatima/Makefile
@@ -1,5 +1,5 @@
include ../Makefile.arch
-
+
#------------------------------------------------------------------------------
SHARELIB = libFatima.so
all: $(SHARELIB)
@@ -7,8 +7,8 @@ all: $(SHARELIB)
############### Detector ##############
## Fatima ##
-libFatima.so: TFatimaData.o TFatimaDataDict.o Fatima.o TFatimaPhysics.o TFatimaPhysicsDict.o
- $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+libFatima.so: TFatimaData.o TFatimaDataDict.o TFatimaPhysics.o TFatimaPhysicsDict.o
+ $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
TFatimaDataDict.cxx: TFatimaData.h
rootcint -f $@ -c $^
@@ -16,12 +16,10 @@ TFatimaDataDict.cxx: TFatimaData.h
TFatimaPhysicsDict.cxx: TFatimaPhysics.h
rootcint -f $@ -c $^
+
# dependances
-Fatima.o: Fatima.cxx Fatima.h
TFatimaData.o: TFatimaData.cxx TFatimaData.h
-TFatimaDataDict.o: TFatimaData.cxx TFatimaData.h
TFatimaPhysics.o: TFatimaPhysics.cxx TFatimaPhysics.h
-TFatimaPhysicsDict.o: TFatimaPhysics.cxx TFatimaPhysics.h
#######################################
diff --git a/NPLib/Fatima/TFatimaData.cxx b/NPLib/Fatima/TFatimaData.cxx
index cd68697d943cb79c037000730ff22c3f21a47a77..16b69f1667bd951abba72bd6e7872b3aeefa9c70 100644
--- a/NPLib/Fatima/TFatimaData.cxx
+++ b/NPLib/Fatima/TFatimaData.cxx
@@ -1,7 +1,5 @@
-#include <iostream>
-#include "TFatimaData.h"
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2013 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 *
@@ -10,8 +8,8 @@
/*****************************************************************************
* Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/2013 *
- * Last update : 02/07/2014 *
+ * Creation Date : 04/12/2009 *
+ * Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class described the raw data of the Fatima detector *
@@ -21,87 +19,34 @@
* *
* *
*****************************************************************************/
-
+#include <iostream>
+#include "TFatimaData.h"
ClassImp(TFatimaData)
-TFatimaData::TFatimaData()
-{
- // Default constructor
- Clear();
+////////////////////////////////////////////////////////////////////////////////
+TFatimaData::TFatimaData(){
}
-
-
-TFatimaData::~TFatimaData()
-{
+TFatimaData::~TFatimaData(){
}
-
-
-void TFatimaData::Clear()
-{
-
- /*
- fFatima_Energy.clear();
- fFatima_Number.clear();
- fFatima_Time.clear();
- */
-
-
- fFATIMA_LaBr3Stage_E_DetectorNbr.clear();
- fFATIMA_LaBr3Stage_E_CrystalNbr.clear();
- fFATIMA_LaBr3Stage_E_Energy.clear();
- // fFATIMA_LaBr3Stage_Eff_phpeak.clear();
- // Time
- fFATIMA_LaBr3Stage_T_DetectorNbr.clear();
- fFATIMA_LaBr3Stage_T_CrystalNbr.clear();
- fFATIMA_LaBr3Stage_T_Time.clear();
-
-
- /*
- // Second Stage CsI
- // CsI
- // Energy
- fFATIMA_CsIStage_E_DetectorNbr.clear();
- fFATIMA_CsIStage_E_CrystalNbr.clear();
- fFATIMA_CsIStage_E_Energy.clear();
- // fFATIMA_CsIStage_Eff_phpeak.clear();
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaData::Clear(){
+ fFATIMA_LaBr3_E_DetectorNbr.clear();
+ fFATIMA_LaBr3_E_Energy.clear();
// Time
- fFATIMA_CsIStage_T_DetectorNbr.clear();
- fFATIMA_CsIStage_T_CrystalNbr.clear();
- fFATIMA_CsIStage_T_Time.clear();
- */
+ fFATIMA_LaBr3_T_DetectorNbr.clear();
+ fFATIMA_LaBr3_T_Time.clear();
}
-
-
-void TFatimaData::Dump() const
-{
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaData::Dump() const{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
- /*
- for (unsigned short i = 0; i<fFatima_Energy.size(); i ++) {
- cout << "Fatima Number " << fFatima_Number[i] << " Energy: " << fFatima_Energy[i] << " Time: "<< fFatima_Time[i] << endl;
- }
- */
-
- for (UShort_t i = 0; i < fFATIMA_LaBr3Stage_E_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fFATIMA_LaBr3Stage_E_DetectorNbr[i] << " Crystal: " << fFATIMA_LaBr3Stage_E_CrystalNbr[i] << " Energy: " << fFATIMA_LaBr3Stage_E_Energy[i] << endl;
- // (X,T)
- cout << "FATIMA_LaBr3Stage_T_Mult = " << fFATIMA_LaBr3Stage_T_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fFATIMA_LaBr3Stage_T_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fFATIMA_LaBr3Stage_T_DetectorNbr[i] << " Crystal: " << fFATIMA_LaBr3Stage_T_CrystalNbr[i] << " Time: " << fFATIMA_LaBr3Stage_T_Time[i] << endl;
-
- /*
- // Second Stage
- // Energy
- cout << "FATIMA_CsIStage_E_Mult = " << fFATIMA_CsIStage_E_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fFATIMA_CsIStage_E_DetectorNbr.size(); i++)
- cout << "Det: " << fFATIMA_CsIStage_E_DetectorNbr[i] << " Crystal: " << fFATIMA_CsIStage_E_CrystalNbr[i] << " Energy: " << fFATIMA_CsIStage_E_Energy[i] << endl;
- // Time
- cout << "FATIMA_CsIStage_T_Mult = " << fFATIMA_CsIStage_T_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fFATIMA_CsIStage_T_DetectorNbr.size(); i++)
- cout << "Det: " << fFATIMA_CsIStage_T_DetectorNbr[i] << " Crystal: " << fFATIMA_CsIStage_T_CrystalNbr[i] << " Time: " << fFATIMA_CsIStage_T_Time[i] << endl;
-
- */
-
+ // E
+ for (UShort_t i = 0; i < fFATIMA_LaBr3_E_DetectorNbr.size(); i++)
+ cout << "DetNbr: " << fFATIMA_LaBr3_E_DetectorNbr[i] << " Energy: " << fFATIMA_LaBr3_E_Energy[i] << endl;
+
+ // T
+ for (UShort_t i = 0; i < fFATIMA_LaBr3_T_DetectorNbr.size(); i++)
+ cout << "DetNbr: " << fFATIMA_LaBr3_T_DetectorNbr[i] << " Time: " << fFATIMA_LaBr3_T_Time[i] << endl;
}
diff --git a/NPLib/Fatima/TFatimaData.h b/NPLib/Fatima/TFatimaData.h
index 1019cc0f2f6bcd8c0c590b771e27f841d89d207c..d6cf7adb05412f45097aca7f663e702bb7404507 100644
--- a/NPLib/Fatima/TFatimaData.h
+++ b/NPLib/Fatima/TFatimaData.h
@@ -1,7 +1,7 @@
#ifndef __FatimaDATA__
#define __FatimaDATA__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -10,8 +10,8 @@
/*****************************************************************************
* Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/2013 *
- * Last update : 02/07/2014 *
+ * Creation Date : 04/12/2009 *
+ * Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class described the raw data of the Fatima detector *
@@ -29,215 +29,77 @@ using namespace std ;
class TFatimaData : public TObject {
- protected:
- // First Stage LaBr
- // LaBr3
- // Energy
- vector<UShort_t> fFATIMA_LaBr3Stage_E_DetectorNbr;
- vector<UShort_t> fFATIMA_LaBr3Stage_E_CrystalNbr;
- vector<Double_t> fFATIMA_LaBr3Stage_E_Energy;
- // vector<Double_t> fFATIMA_LaBr3Stage_Eff_phpeak;
- // Time
- vector<UShort_t> fFATIMA_LaBr3Stage_T_DetectorNbr;
- vector<UShort_t> fFATIMA_LaBr3Stage_T_CrystalNbr;
- vector<Double_t> fFATIMA_LaBr3Stage_T_Time;
-
-
- /*
- // Second Stage CsI
- // CsI
- // Energy
- vector<UShort_t> fFATIMA_CsIStage_E_DetectorNbr;
- vector<UShort_t> fFATIMA_CsIStage_E_CrystalNbr;
- vector<Double_t> fFATIMA_CsIStage_E_Energy;
- // vector<Double_t> fFATIMA_CsIStage_Eff_phpeak;
- // Time
- vector<UShort_t> fFATIMA_CsIStage_T_DetectorNbr;
- vector<UShort_t> fFATIMA_CsIStage_T_CrystalNbr;
- vector<Double_t> fFATIMA_CsIStage_T_Time;
- */
-
- /*
- private:
- vector<double> fFatima_Energy;
- vector<double> fFatima_Time;
- vector<short> fFatima_Number;
- */
-
- public:
- TFatimaData();
- virtual ~TFatimaData();
-
- void Clear();
- void Clear(const Option_t*) {};
- void Dump() const;
-
- ///////////////////// GETTERS ////////////////////////
-
- //
- // First stage
- //
- // (E)
- UShort_t GetFATIMALaBr3StageEMult() {
- return fFATIMA_LaBr3Stage_E_DetectorNbr.size(); // TODO: Maybe change to CrystalNbr for multiplicity
- }
- UShort_t GetFATIMALaBr3StageEDetectorNbr(Int_t i) {
- return fFATIMA_LaBr3Stage_E_DetectorNbr.at(i);
- }
- UShort_t GetFATIMALaBr3StageECrystalNbr(Int_t i) {
- return fFATIMA_LaBr3Stage_E_CrystalNbr.at(i);
- }
- Double_t GetFATIMALaBr3StageEEnergy(Int_t i) {
- return fFATIMA_LaBr3Stage_E_Energy.at(i);
- }
- /*
- Double_t GetFATIMALaBr3StageEffphpeak(Int_t i) {
- return fFATIMA_LaBr3Stage_Eff_phpeak.at(i);
- }
- */
-
- // (T)
- UShort_t GetFATIMALaBr3StageTMult() {
- return fFATIMA_LaBr3Stage_E_DetectorNbr.size();
- }
- UShort_t GetFATIMALaBr3StageTDetectorNbr(Int_t i) {
- return fFATIMA_LaBr3Stage_T_DetectorNbr.at(i);
- }
- UShort_t GetFATIMALaBr3StageTCrystalNbr(Int_t i) {
- return fFATIMA_LaBr3Stage_T_CrystalNbr.at(i);
- }
- Double_t GetFATIMALaBr3StageTTime(Int_t i) {
- return fFATIMA_LaBr3Stage_T_Time.at(i);
- }
-
- /*
- //
- // Second stage (CsI
- //
- // (E)
- UShort_t GetFATIMACsIStageEMult() {
- return fFATIMA_CsIStage_E_DetectorNbr.size();
- }
- UShort_t GetFATIMACsIStageEDetectorNbr(Int_t i) {
- return fFATIMA_CsIStage_E_DetectorNbr.at(i);
- }
- UShort_t GetFATIMACsIStageECrystalNbr(Int_t i) {
- return fFATIMA_CsIStage_E_CrystalNbr.at(i);
- }
- Double_t GetFATIMACsIStageEEnergy(Int_t i) {
- return fFATIMA_CsIStage_E_Energy.at(i);
- }
- */
- /*
- Double_t GetFATIMACsIStageEffphpeak(Int_t i) {
- return fFATIMA_CsIStage_Eff_phpeak.at(i);
- }
- */
- /*
- // (T)
- UShort_t GetFATIMACsIStageTMult() {
- return fFATIMA_CsIStage_E_DetectorNbr.size();
- }
- UShort_t GetFATIMACsIStageTDetectorNbr(Int_t i) {
- return fFATIMA_CsIStage_T_DetectorNbr.at(i);
- }
- UShort_t GetFATIMACsIStageTCrystalNbr(Int_t i) {
- return fFATIMA_CsIStage_T_CrystalNbr.at(i);
- }
- Double_t GetFATIMACsIStageTTime(Int_t i) {
- return fFATIMA_CsIStage_T_Time.at(i);
- }
- */
-
- /*
- // (E)
- //double GetEnergy(int i) {return fFatima_Energy[i];}
- // (T)
- //double GetTime(int i) {return fFatima_Time[i];}
- // (N)
- int GetFatimaNumber(int i) {return fFatima_Number[i];}
- double GetEnergySize() {return fFatima_Energy.size();}
- // (T)
- double GetTimeSize() {return fFatima_Time.size();}
- // (N)
- int GetFatimaNumberSize() {return fFatima_Number.size();}
- */
-
- ///////////////////// SETTERS ////////////////////////
-
- //
- // First stage
- //
- // (E)
-
- void SetFATIMALaBr3StageEDetectorNbr(UShort_t DetNbr) {
- fFATIMA_LaBr3Stage_E_DetectorNbr.push_back(DetNbr);
- }
- void SetFATIMALaBr3StageECrystalNbr(UShort_t CrystalNbr) {
- fFATIMA_LaBr3Stage_E_CrystalNbr.push_back(CrystalNbr);
- }
- void SetFATIMALaBr3StageEEnergy(Double_t Energy) {
- fFATIMA_LaBr3Stage_E_Energy.push_back(Energy);
- }
- /*
- void SetFATIMALaBr3StageEffphpeak(Double_t Energy) {
- fFATIMA_LaBr3Stage_Eff_phpeak.push_back(Energy);
- }
- */
-
- // (T)
- void SetFATIMALaBr3StageTDetectorNbr(UShort_t DetNbr) {
- fFATIMA_LaBr3Stage_T_DetectorNbr.push_back(DetNbr);
- }
- void SetFATIMALaBr3StageTCrystalNbr(UShort_t CrystalNbr) {
- fFATIMA_LaBr3Stage_T_CrystalNbr.push_back(CrystalNbr);
- }
- void SetFATIMALaBr3StageTTime(Double_t Time) {
- fFATIMA_LaBr3Stage_T_Time.push_back(Time);
- }
-
- /*
- //
- // Second stage (CsI
- //
- // (E)
- void SetFATIMACsIStageEDetectorNbr(UShort_t DetNbr) {
- fFATIMA_CsIStage_E_DetectorNbr.push_back(DetNbr);
- }
- void SetFATIMACsIStageECrystalNbr(UShort_t CrystalNbr) {
- fFATIMA_CsIStage_E_CrystalNbr.push_back(CrystalNbr);
- }
- void SetFATIMACsIStageEEnergy(Double_t Energy) {
- fFATIMA_CsIStage_E_Energy.push_back(Energy);
- }
- */
-
- /*
- void SetFATIMACsIStageEffphpeak(Double_t Energy) {
- fFATIMA_CsIStage_Eff_phpeak.push_back(Energy);
- }
- */
-
- /*
- // (T)
- void SetFATIMACsIStageTDetectorNbr(UShort_t DetNbr) {
- fFATIMA_CsIStage_T_DetectorNbr.push_back(DetNbr);
- }
- void SetFATIMACsIStageTCrystalNbr(UShort_t CrystalNbr) {
- fFATIMA_CsIStage_T_CrystalNbr.push_back(CrystalNbr);
- }
- void SetFATIMACsIStageTTime(Double_t Time) {
- fFATIMA_CsIStage_T_Time.push_back(Time);
- }
- */
-
- /*
- // (E)
- void SetEnergy(double E) {fFatima_Energy.push_back(E);}
- void SetTime(double T) {fFatima_Time.push_back(T);}
- void SetFatimaNumber(int N) {fFatima_Number.push_back(N);}
- */
- ClassDef(TFatimaData,1) // FatimaData structure
+ protected:
+ // LaBr3
+ // Energy
+ vector<UShort_t> fFATIMA_LaBr3_E_DetectorNbr;
+ vector<Double_t> fFATIMA_LaBr3_E_Energy;
+ // Time
+ vector<UShort_t> fFATIMA_LaBr3_T_DetectorNbr;
+ vector<Double_t> fFATIMA_LaBr3_T_Time;
+
+ public:
+ TFatimaData();
+ virtual ~TFatimaData();
+
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+ ///////////////////// GETTERS ////////////////////////
+ // (E)
+ UShort_t GetFatimaLaBr3EMult() {
+ return fFATIMA_LaBr3_E_DetectorNbr.size();
+ }
+
+ UShort_t GetFatimaLaBr3EDetectorNbr(Int_t i) {
+ return fFATIMA_LaBr3_E_DetectorNbr.at(i);
+ }
+ Double_t GetFatimaLaBr3EEnergy(Int_t i) {
+ return fFATIMA_LaBr3_E_Energy.at(i);
+ }
+
+
+ // (T)
+ UShort_t GetFatimaLaBr3TMult() {
+ return fFATIMA_LaBr3_E_DetectorNbr.size();
+ }
+
+ UShort_t GetFatimaLaBr3TDetectorNbr(Int_t i) {
+ return fFATIMA_LaBr3_T_DetectorNbr.at(i);
+ }
+ Double_t GetFatimaLaBr3TTime(Int_t i) {
+ return fFATIMA_LaBr3_T_Time.at(i);
+ }
+
+ ///////////////////// SETTERS ////////////////////////
+ // (E)
+ void SetFatimaLaBr3E(UShort_t DetectorNbr, Double_t E){
+ fFATIMA_LaBr3_E_DetectorNbr.push_back(DetectorNbr);
+ fFATIMA_LaBr3_E_Energy.push_back(E);
+ }
+
+ void SetFatimaLaBr3EDetectorNbr(UShort_t DetectorNbr) {
+ fFATIMA_LaBr3_E_DetectorNbr.push_back(DetectorNbr);
+ }
+ void SetFatimaLaBr3EEnergy(Double_t Energy) {
+ fFATIMA_LaBr3_E_Energy.push_back(Energy);
+ }
+
+ // (T)
+ void SetFatimaLaBr3T(UShort_t DetectorNbr, Double_t T){
+ fFATIMA_LaBr3_T_DetectorNbr.push_back(DetectorNbr);
+ fFATIMA_LaBr3_T_Time.push_back(T);
+ }
+
+ void SetFatimaLaBr3TDetectorNbr(UShort_t DetectorNbr) {
+ fFATIMA_LaBr3_T_DetectorNbr.push_back(DetectorNbr);
+ }
+ void SetFatimaLaBr3TTime(Double_t Time) {
+ fFATIMA_LaBr3_T_Time.push_back(Time);
+ }
+
+ ClassDef(TFatimaData,1) // FatimaData structure
};
#endif
diff --git a/NPLib/Fatima/TFatimaPhysics.cxx b/NPLib/Fatima/TFatimaPhysics.cxx
index 199d022012349af24ddacde804561a3f23698248..286e6e386e0f905cc04563df829ded31e53e011f 100644
--- a/NPLib/Fatima/TFatimaPhysics.cxx
+++ b/NPLib/Fatima/TFatimaPhysics.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2013 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 *
@@ -8,8 +8,8 @@
/*****************************************************************************
* Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/2014 *
+ * Creation Date : 04/12/09 *
+ * Last update : 8/12/14 by A. MATTA *
*---------------------------------------------------------------------------*
* Decription: This class stores the physical results after NPAnalysis is run*
* for the FATIMA detector. *
@@ -17,340 +17,552 @@
* stored in the output TTree of NPAnalysis. *
*---------------------------------------------------------------------------*
* Comment: *
- * *
+ * Last update: Merging of Fatima and TFatimaPhyscis class *
* *
*****************************************************************************/
+// NPL
#include "TFatimaPhysics.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+
+// STL
#include <vector>
#include <iostream>
+#include <fstream>
#include <cstdlib>
+#include <limits>
+#include <cmath>
+// Root
+#include "TChain.h"
+#include "TTree.h"
+
+using namespace std;
ClassImp(TFatimaPhysics)
+ ////////////////////////////////////////////////////////////////////////////////
+ TFatimaPhysics::TFatimaPhysics(){
+ m_NumberOfModule = 0;
+ m_EventData = new TFatimaData();
+ m_EventPhysics = new TFatimaPhysics();
+ }
-TFatimaPhysics::TFatimaPhysics()
-{
- //EventMultiplicity = 0;
+////////////////////////////////////////////////////////////////////////////////
+TFatimaPhysics::~TFatimaPhysics(){
+ Clear();
}
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::BuildSimplePhysicalEvent(){
+ BuildPhysicalEvent();
+}
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::BuildPhysicalEvent(){
-TFatimaPhysics::~TFatimaPhysics()
-{
- Clear();
+ int multLaBrE = m_EventData->GetFatimaLaBr3EMult();
+
+ FatimaEventMult=multLaBrE;
+ // get energy from strips and store it
+ if(FatimaEventMult>=1){
+ double EnergyTot=0.;
+ double EnergyStripFront;
+ double EnergyStrip;
+
+ for(int j=0;j<multLaBrE;j++){
+ EnergyStripFront= m_EventData->GetFatimaLaBr3EEnergy(j);
+
+ EnergyStrip = EnergyStripFront;
+ FatimaLaBr_E.push_back(EnergyStrip);
+
+ EnergyTot += EnergyStrip;
+ //cout << "Energytot LaBr=" << EnergyTot << endl;
+ }
+
+ // Fill total energy in inner shell
+ FatimaInTotalEnergy.push_back(EnergyTot);
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::Clear(){
+ //EventMultiplicity= 0;
+ FatimaEventMult= 0;
+ //ModuleNumber.clear();
+ //EventType.clear();
+ FatimaInTotalEnergy.clear(); // inner shell
+ FatimaOutTotalEnergy.clear(); // outter shell
+ FatimaTotalEnergy.clear();
+
+ // LaBr
+ FatimaLaBr_E.clear();
+ //First_T.clear();
+ //First_X.clear();
+ //First_Y.clear();
+
+ // NaI
+ FatimaNaI_E.clear();
+ //Second_T.clear();
+ //Second_N.clear();
+
+ /*
+ // NaI
+ Third_E.clear();
+ Third_T.clear();
+ Third_N.clear();
+ */
}
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::ReadConfiguration(string Path) {
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+ string LineBuffer ;
+ string DataBuffer ;
+
+ double Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Az, Bz, Cz, Dz;
+ TVector3 A, B, C, D;
+ double Theta = 0, Phi = 0, R = 0, beta_u = 0 , beta_v = 0 , beta_w = 0;
+
+ bool check_A = false;
+ bool check_C = false;
+ bool check_B = false;
+ bool check_D = false;
+
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+ bool check_beta = false;
+
+ bool ReadingStatus = false;
+
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+
+ // If line is a Fatima bloc, reading toggle to true
+ // and toggle to true flags indicating which shape is treated.
+ if (LineBuffer.compare(0, 12, "FatimaCluster") == 0 ) {
+ cout << "///////////////////////" << endl;
+ cout << "Detector found:" << endl;
+ ReadingStatus = true;
+ }
+ // Reading Block
+ while (ReadingStatus) {
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {
+ ConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
+ }
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer.compare(0, 13, "FatimaDetector") == 0) {
+ cout << "WARNING: Another Module is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
+ ReadingStatus = false;
+ }
+ // Position method
+ else if (DataBuffer=="A=") {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az ;
+
+ A = TVector3(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+ }
+ else if (DataBuffer=="B=") {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz ;
+
+ B = TVector3(Bx, By, Bz);
+ cout << "X128 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+ }
+ else if (DataBuffer == "C=") {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz ;
+
+ C = TVector3(Cx, Cy, Cz);
+ cout << "X1 Y128 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
+ }
+ else if (DataBuffer=="D=") {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz ;
+
+ D = TVector3(Dx, Dy, Dz);
+ cout << "X128 Y128 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+ } // End Position Method
+
+ // Angle method
+ else if (DataBuffer.compare(0, 6, "Theta=") == 0) {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta ;
+ cout << "Theta: " << Theta << endl;
+ }
+ else if (DataBuffer.compare(0, 4, "Phi=") == 0) {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi ;
+ cout << "Phi: " << Phi << endl;
+ }
+ else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R ;
+ cout << "R: " << R << endl;
+ }
+ else if (DataBuffer.compare(0, 5, "Beta=") == 0) {
+ check_beta = true;
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w ;
+ cout << "Beta: " << beta_u << " " << beta_v << " " << beta_w << endl ;
+ }
-void TFatimaPhysics::BuildSimplePhysicalEvent(TFatimaData* Data)
-{
- BuildPhysicalEvent(Data);
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) ) {
+ ReadingStatus = false;
+
+ // Add The previously define telescope
+ // With position method
+ if ( check_A && check_B && check_C && check_D ) {
+ AddDetector(A ,
+ B ,
+ C ,
+ D ) ;
+ }
+
+ // with angle method
+ else if ( check_Theta && check_Phi && check_R && check_beta ) {
+ AddDetector(Theta,
+ Phi,
+ R,
+ beta_u,
+ beta_v,
+ beta_w);
+ }
+
+ // reset boolean flag for point positioning
+ check_A = false;
+ check_B = false;
+ check_C = false;
+ check_D = false;
+
+ // reset boolean flag for angle positioning
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
+ check_beta = false;
+
+ } // end test for adding a module
+
+
+ } // end while for reading block
+ } // end while for reading file
+
+ cout << endl << "/////////////////////////////" << endl<<endl;
}
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::ReadCalibrationFile(string Path){
+ // Order of Polynom function used for calibration
+ /* int Calibration_Si_E_Order;
+ int Calibration_Si_T_Order;
+ int Calibration_SiLi_E_Order;
+ int Calibration_NaI_E_Order;
+ */
+ // Calibration_Si_X_E[DetectorNumber][StripNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_Si_X_E ;
+ vector< vector< vector< double > > > Calibration_Si_X_T ;
+ vector< vector< vector< double > > > Calibration_Si_Y_E ;
+ vector< vector< vector< double > > > Calibration_Si_Y_T ;
+
+ // Calibration_SiLi_E[DetectorNumber][PadNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_SiLi_E ;
+
+ // Calibration_SiLi_E[DetectorNumber][CrystalNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_NaI_E ;
+
+ if (Path == "Simulation") { // Simulation case: data already calibrated
+ /* Calibration_Si_E_Order = 1;
+ Calibration_Si_T_Order = 1;
+ Calibration_SiLi_E_Order = 1;
+ Calibration_NaI_E_Order = 1;
+ */
+ vector<double> Coef;
+ // Order 0 Order 1
+ Coef.push_back(0) ; Coef.push_back(1) ;
+
+ vector< vector<double> > StripLine ;
+ StripLine.resize( 128 , Coef) ;
+
+ Calibration_Si_X_E.resize( m_NumberOfModule , StripLine) ;
+ Calibration_Si_X_T.resize( m_NumberOfModule , StripLine) ;
+ Calibration_Si_Y_E.resize( m_NumberOfModule , StripLine) ;
+ Calibration_Si_Y_T.resize( m_NumberOfModule , StripLine) ;
+
+ Calibration_SiLi_E.resize( m_NumberOfModule , StripLine) ;
+ Calibration_NaI_E .resize( m_NumberOfModule , StripLine) ;
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::InitializeRootInputRaw(){
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("FATIMA", true);
+ inputChain->SetBranchStatus("fFATIMA*", true);
+ inputChain->SetBranchAddress("FATIMA", &m_EventData);
+}
-void TFatimaPhysics::BuildPhysicalEvent(TFatimaData* Data)
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::InitializeRootOutput(){
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("FATIMA", "TFatimaPhysics", &m_EventPhysics);
+}
+void TFatimaPhysics::AddDetector(TVector3 A,
+ TVector3 B,
+ TVector3 C,
+ TVector3 D)
{
+ m_NumberOfModule++;
- int multLaBrE = Data->GetFATIMALaBr3StageEMult();
- // int multCsIE = Data->GetFATIMACsIStageEMult();
+ // remove warning using D
+ D.Unit();
- cout << "%%%%%%%%%%%%%%%%%%%%%%%" << endl;
- cout << "multLaBr= " << multLaBrE << endl;
- //cout << "multCsI= " << multCsIE << endl;
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U = B - A;
+ U = U.Unit();
- FatimaEventMult=multLaBrE;
- //FatimaEventMult=multLaBrE+multCsIE;
-
- // Get energy from strips and store it
-
- if(FatimaEventMult>=1)
- {
-
- double EnergyTot=0.;
-
- double HighestEnergy;
- int DetID_HighestE;
- double CrystID_HighestE;
- // Initialisation
- HighestEnergy=0;
- DetID_HighestE=-1;
- CrystID_HighestE=-1;
-
- /*
- double HighestEnergy, Outer_HighestEnergy;
- int DetID_HighestE, OuterDetID_HighestE;
- double CrystID_HighestE, OuterCrystID_HighestE;
-
- // Initialisation
- HighestEnergy=Outer_HighestEnergy=0;
- DetID_HighestE=OuterDetID_HighestE=-1;
- CrystID_HighestE=OuterCrystID_HighestE=-1;
- */
-
- if(multLaBrE>=1){
- //cout << "cava1b" <<endl;
- //cout << Data->GetFATIMALaBr3StageEEnergy(0) <<endl;
- //cout << "cava1b" <<endl;
-
- double EnergyTotFirst;
- double EnergyStripFront;
- double EnergyStrip;
- //double HighestEnergyStrip;
-
- for(int j=0;j<multLaBrE;j++)
- {
- EnergyStripFront= Data->GetFATIMALaBr3StageEEnergy(j);
-
- EnergyStrip = EnergyStripFront;
- cout << "Energy LaBr=" << EnergyStrip << endl;
- FatimaLaBr_E.push_back(EnergyStrip);
-
- EnergyTotFirst += EnergyStrip;
- EnergyTot += EnergyStrip;
-
- //Let's find the first detector hit:
- if(HighestEnergy<EnergyStripFront)
- {
- HighestEnergy=EnergyStripFront;
- DetID_HighestE=Data->GetFATIMALaBr3StageEDetectorNbr(j);
- CrystID_HighestE=Data->GetFATIMALaBr3StageECrystalNbr(j);
- }
-
-
- }
-
- // Fill total energy in inner shell and the first detector/crystal hit
-
- FatimaInTotalEnergy.push_back(EnergyTotFirst);
-
- cout << " Id of Fatima detector hit= " << DetID_HighestE << endl;
- cout << " EnergyTotFirst= " << EnergyTotFirst << endl;
- cout << " EnergyTot= " << EnergyTot << endl;
-
- if(DetID_HighestE>100)// Cluster case
- {
- //ModuleNumber.push_back(DetID_HighestE);
- DetectorNumber.push_back(-1);
- }else // Detector case
- {
- //ModuleNumber.push_back(-1);
- DetectorNumber.push_back(DetID_HighestE);
- }
-
- if(DetID_HighestE<100) // Detector case
- {
- FirstStage_X.push_back(0);
- FirstStage_Y.push_back(0);
- }
-
-
- /*
- if(DetID_HighestE>100) // Cluster case
- {
-
- if(CrystID_HighestE==0)
- {
- FirstStage_X.push_back(0);
- FirstStage_Y.push_back(0);
- }
- if(CrystID_HighestE==1)
- {
- FirstStage_X.push_back(1);
- FirstStage_Y.push_back(0);
- }
- if(CrystID_HighestE==2)
- {
- FirstStage_X.push_back(2);
- FirstStage_Y.push_back(0);
- }
- if(CrystID_HighestE==3)
- {
- FirstStage_X.push_back(0);
- FirstStage_Y.push_back(1);
- }
- if(CrystID_HighestE==4)
- {
- FirstStage_X.push_back(1);
- FirstStage_Y.push_back(1);
- }
- if(CrystID_HighestE==5)
- {
- FirstStage_X.push_back(2);
- FirstStage_Y.push_back(1);
- }
- if(CrystID_HighestE==6)
- {
- FirstStage_X.push_back(0);
- FirstStage_Y.push_back(2);
- }
- if(CrystID_HighestE==7)
- {
- FirstStage_X.push_back(1);
- FirstStage_Y.push_back(2);
- }
- if(CrystID_HighestE==8)
- {
- FirstStage_X.push_back(2);
- FirstStage_Y.push_back(2);
- }
- } // end of cluster case
- */
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V = C - A;
+ V = V.Unit();
- }
+ // Position Vector of Strip Center
+ TVector3 StripCenter = TVector3(0,0,0);
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1;
- /*
- if(multCsIE>=1){
- double EnergySecond;
- double EnergyTotSecond;
- for(int j=0;j<multCsIE;j++)
- {
- EnergySecond = Data->GetFATIMACsIStageEEnergy(j);
- FatimaCsI_E.push_back(EnergySecond);
- EnergyTotSecond +=EnergySecond;
-
- EnergyTot += EnergySecond;
- cout << "Energy CsI=" << EnergySecond << endl;
- //cout << "Energytot CsI=" << EnergyTot << endl;
-
- //Let's find the first outer detector hit:
- if(Outer_HighestEnergy<EnergySecond)
- {
- Outer_HighestEnergy=EnergySecond;
- OuterDetID_HighestE=Data->GetFATIMACsIStageEDetectorNbr(j);
- OuterCrystID_HighestE=Data->GetFATIMACsIStageECrystalNbr(j);
- }
- }
-
- // Fill total energy in outer shell
- FatimaOutTotalEnergy.push_back(EnergyTotSecond);
-
- cout << " Id of Fatima outer detector hit= " << OuterDetID_HighestE << endl;
- cout << " EnergyTotSecond= " << EnergyTotSecond << endl;
- cout << " EnergyTot= " << EnergyTot << endl;
-
- if(multLaBrE==0 || (EnergyTotSecond>EnergyTot) ) { // = only if Outer layer is hit first
-
- FirstStage_X.push_back(-1);
- FirstStage_Y.push_back(-1);
-
- cout << " Id of Fatima (CsI) detector hit= " << OuterDetID_HighestE << endl;
-
- if(OuterDetID_HighestE>100)
- {
- ModuleNumber.push_back(OuterDetID_HighestE);
- DetectorNumber.push_back(-1);
- }else
- {
- ModuleNumber.push_back(-1);
- DetectorNumber.push_back(OuterDetID_HighestE);
- }
-
-
- if(OuterDetID_HighestE<100) // Detector case
- {
- SecondStage_X.push_back(0);
- SecondStage_Y.push_back(0);
- }
-
-
- if(OuterDetID_HighestE>100) // Cluster case
- {
- if(OuterCrystID_HighestE==0)
- {
- SecondStage_X.push_back(0);
- SecondStage_Y.push_back(0);
- }
- if(OuterCrystID_HighestE==1)
- {
- SecondStage_X.push_back(1);
- SecondStage_Y.push_back(0);
- }
- if(OuterCrystID_HighestE==2)
- {
- SecondStage_X.push_back(2);
- SecondStage_Y.push_back(0);
- }
- if(OuterCrystID_HighestE==3)
- {
- SecondStage_X.push_back(0);
- SecondStage_Y.push_back(1);
- }
- if(OuterCrystID_HighestE==4)
- {
- SecondStage_X.push_back(1);
- SecondStage_Y.push_back(1);
- }
- if(OuterCrystID_HighestE==5)
- {
- SecondStage_X.push_back(2);
- SecondStage_Y.push_back(1);
- }
- if(OuterCrystID_HighestE==6)
- {
- SecondStage_X.push_back(0);
- SecondStage_Y.push_back(2);
- }
- if(OuterCrystID_HighestE==7)
- {
- SecondStage_X.push_back(1);
- SecondStage_Y.push_back(2);
- }
- if(OuterCrystID_HighestE==8)
- {
- SecondStage_X.push_back(2);
- SecondStage_Y.push_back(2);
- }
-
- }
- }
+ // Geometry Parameter
+ double Face = 98; // mm
+ double NumberOfStrip = 128;
+ double StripPitch = Face/NumberOfStrip; // mm
- }
- */
+ // Buffer object to fill Position Array
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModulePositionX;
+ vector< vector< double > > OneModulePositionY;
+ vector< vector< double > > OneModulePositionZ;
+
+ // Moving StripCenter to 1.1 corner:
+ Strip_1_1 = A + (U+V) * (StripPitch/2.);
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
+ // StripCenter += -TargetPosition;
- // Fill total energy
- FatimaTotalEnergy.push_back(EnergyTot);
+ lineX.push_back( StripCenter.X() );
+ lineY.push_back( StripCenter.Y() );
+ lineZ.push_back( StripCenter.Z() );
}
-}
+ OneModulePositionX.push_back(lineX);
+ OneModulePositionY.push_back(lineY);
+ OneModulePositionZ.push_back(lineZ);
+ }
+ m_PositionX.push_back( OneModulePositionX );
+ m_PositionY.push_back( OneModulePositionY );
+ m_PositionZ.push_back( OneModulePositionZ );
+}
-void TFatimaPhysics::Clear()
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::AddDetector(double theta,
+ double phi,
+ double R,
+ double beta_u,
+ double beta_v,
+ double beta_w)
{
- //EventMultiplicity= 0;
- FatimaEventMult= 0;
- DetectorNumber.clear(); // phoswich #
-
- //EventType.clear();
- FatimaInTotalEnergy.clear(); // inner shell
- //FatimaOutTotalEnergy.clear(); // outter shell
- FatimaTotalEnergy.clear();
-
-
- // LaBr
- FatimaLaBr_E.clear();
- //FatimaIn_1stDetId.clear();
- //FatimaIn_1stCrystId.clear();
- //FirstStage_T.clear();
- FirstStage_X.clear();
- FirstStage_Y.clear();
-
- /* // CsI
- FatimaCsI_E.clear();
- //FatimaOut_1stDetId.clear();
- //FatimaOut_1stCrystId.clear();
- //SecondStage_T.clear();
- //SecondStage_N.clear();
- SecondStage_X.clear();
- SecondStage_Y.clear();
- */
-
- /*
- // CsI
- ThirdStage_E.clear();
- ThirdStage_T.clear();
- ThirdStage_N.clear();
- */
+ m_NumberOfModule++;
+
+ // convert from degree to radian:
+ double Pi = 3.141592654;
+ theta = theta * Pi/180. ;
+ phi = phi * Pi/180. ;
+
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U ;
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V ;
+ // Vector W normal to Module Face (pointing NaI)
+ TVector3 W ;
+ // Vector position of Module Face center
+ TVector3 C ;
+
+ C = TVector3(R * sin(theta) * cos(phi),
+ R * sin(theta) * sin(phi),
+ R * cos(theta));
+
+ TVector3 YperpC = TVector3( cos(theta) * cos(phi),
+ cos(theta) * sin(phi),
+ -sin(theta));
+
+ W = C.Unit();
+ U = W.Cross(YperpC);
+ V = W.Cross(U);
+
+ U = U.Unit();
+ V = V.Unit();
+
+ U.Rotate( beta_u * Pi/180. , U ) ;
+ V.Rotate( beta_u * Pi/180. , U ) ;
+
+ U.Rotate( beta_v * Pi/180. , V ) ;
+ V.Rotate( beta_v * Pi/180. , V ) ;
+
+ U.Rotate( beta_w * Pi/180. , W ) ;
+ V.Rotate( beta_w * Pi/180. , W ) ;
+
+ double Face = 98; // mm
+ double NumberOfStrip = 128;
+ double StripPitch = Face/NumberOfStrip; // mm
+
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModulePositionX;
+ vector< vector< double > > OneModulePositionY;
+ vector< vector< double > > OneModulePositionZ;
+
+ double X, Y, Z;
+
+ // Moving C to the 1.1 corner:
+ C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
+ C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
+ C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j );
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j );
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j );
+
+ lineX.push_back(X);
+ lineY.push_back(Y);
+ lineZ.push_back(Z);
+ }
+
+ OneModulePositionX.push_back(lineX);
+ OneModulePositionY.push_back(lineY);
+ OneModulePositionZ.push_back(lineZ);
+ }
+
+ m_PositionX.push_back( OneModulePositionX );
+ m_PositionY.push_back( OneModulePositionY );
+ m_PositionZ.push_back( OneModulePositionZ );
}
+
+////////////////////////////////////////////////////////////////////////////////
+double TFatimaPhysics::GetEnergyDeposit(){
+ if (m_EventPhysics->FatimaTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaTotalEnergy[0];
+ else
+ return -1000;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double TFatimaPhysics::GetEnergyInDeposit(){// inner Layer
+ if (m_EventPhysics->FatimaInTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaInTotalEnergy[0];
+ else
+ return -1000;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double TFatimaPhysics::GetEnergyOutDeposit(){// Outer Layer
+ if (m_EventPhysics->FatimaOutTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaOutTotalEnergy[0];
+ else
+ return -1000;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TVector3 TFatimaPhysics::GetPositionOfInteraction(){
+ TVector3 Position = TVector3(-1000,-1000,-1000);
+ return(Position);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TFatimaPhysics::Print(){
+ cout << "Number of Modules: " << m_NumberOfModule << endl;
+
+ for (int f = 0; f < m_NumberOfModule; f++) {
+ cout << "Module " << f+1 << endl;
+
+ for (int i = 0; i < 128; i++) {
+ for (int j = 0; j < 128; j++) {
+ cout << i+1 << " "<< j+1 << " "
+ << m_PositionX[f][i][j] << " "
+ << m_PositionY[f][i][j] << " "
+ << m_PositionZ[f][i][j] << " "
+ << endl ;
+ }
+ }
+ }
+}
+
diff --git a/NPLib/Fatima/TFatimaPhysics.h b/NPLib/Fatima/TFatimaPhysics.h
index 4a4592765cdf20c9955f34dc15e392f846eab862..713d26b483c5954565c45156986d30b42e8f5623 100644
--- a/NPLib/Fatima/TFatimaPhysics.h
+++ b/NPLib/Fatima/TFatimaPhysics.h
@@ -1,5 +1,7 @@
+#ifndef TFATIMAPHYSICS_H
+#define TFATIMAPHYSICS_H
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -8,8 +10,8 @@
/*****************************************************************************
* Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/2014 *
+ * Creation Date : 04/12/09 *
+ * Last update : *
*---------------------------------------------------------------------------*
* Decription: This class stores the physical results after NPAnalysis is run*
* for the FATIMA detector. *
@@ -21,71 +23,127 @@
* *
*****************************************************************************/
-#ifndef TFATIMAPHYSICS_H
-#define TFATIMAPHYSICS_H
-
+// STD
#include <vector>
-#include "TObject.h"
-#include "TFatimaData.h"
#include <cstdlib>
-
using namespace std ;
-class TFatimaPhysics : public TObject
-{
-public:
- TFatimaPhysics();
- ~TFatimaPhysics();
-
-public:
- void Clear();
- void Clear(const Option_t*) {};
- void BuildPhysicalEvent(TFatimaData* Data);
- void BuildSimplePhysicalEvent(TFatimaData* Data);
-
-public:
- // Provide Physical Multiplicity
- Int_t FatimaEventMult;
-
- // Provide a Classification of Event
- //vector<int> EventType;
-
-
- // Telescope
- //vector<int> ModuleNumber; // Id of cluster
-
- vector<int> DetectorNumber; // Id of Detector
-
- // FirstStage
- vector<double> FatimaLaBr_E;
- //vector<double> FirstStage_T;
- vector<int> FirstStage_X; // column # of LaBr crystal in 3x3 cluster
- vector<int> FirstStage_Y; // raw # of LaBr crystal in 3x3 cluster
-
- /*
- // SecondStage
- vector<double> FatimaCsI_E;
- //vector<double> SecondStage_T;
- //vector<int> SecondStage_N;
- vector<int> SecondStage_X; // column # of CsI crystal in 3x3 cluster
- vector<int> SecondStage_Y; // raw # of CsI crystal in 3x3 cluster
- */
-
-
- /*
- // ThirdStage
- vector<double> ThirdStage_E;
- vector<double> ThirdStage_T;
- vector<int> ThirdStage_N;
- */
-
- // Physical Value
- vector<double> FatimaTotalEnergy;
- vector<double> FatimaInTotalEnergy;
- vector<double> FatimaOutTotalEnergy;
-
+// NPL
+#include "../include/VDetector.h"
+#include "TFatimaData.h"
- ClassDef(TFatimaPhysics,1) // GaspardTrackerPHysics structure
+// Root
+#include "TObject.h"
+#include "TVector3.h"
+
+
+class TFatimaPhysics : public TObject, public NPA::VDetector{
+ public:
+ TFatimaPhysics();
+ ~TFatimaPhysics();
+
+ public:
+ void Clear();
+ void Clear(const Option_t*) {};
+
+ public:
+ /////////////////////////////////////
+ // Innherited from VDetector Class //
+ /////////////////////////////////////
+ // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+ void ReadConfiguration(string);
+
+ // Read stream at CalibFile and pick-up calibration parameter using Token
+ // If argument is "Simulation" no change calibration is loaded
+ void ReadCalibrationFile(string);
+
+ // 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
+ void InitializeRootInputRaw();
+
+ // Create associated branches and associated private member DetectorPhysics address
+ void InitializeRootOutput();
+
+ // This method is called at each event read from the Input Tree.
+ // The aim 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 ...).
+ // This method aimed to be used for analysis performed during experiment, when speed is requiered.
+ // NB: This method can eventually be the same as BuildPhysicalEvent.
+ void BuildSimplePhysicalEvent();
+
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {m_EventPhysics->Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+ public:
+ ////////////////////////////////
+ // Specific to Fatima //
+ ////////////////////////////////
+ // Case of a Square module
+ // Add a Module using Corner Coordinate information
+ void AddDetector(TVector3 A,TVector3 B,TVector3 C,TVector3 D);//!
+
+ // Add a Module using R Theta Phi of Si center information
+ void AddDetector(double theta,
+ double phi,
+ double r,
+ double beta_u,
+ double beta_v,
+ double beta_w);//!
+
+ // Getters to retrieve the (X,Y,Z) coordinates of a pixel defined by strips (X,Y)
+ double GetPositionX(int N ,int X ,int Y) { return m_PositionX[N-1][X-1][Y-1]; }
+ double GetPositionY(int N ,int X ,int Y) { return m_PositionY[N-1][X-1][Y-1]; }
+ double GetPositionZ(int N ,int X ,int Y) { return m_PositionZ[N-1][X-1][Y-1]; }
+ double GetNumberOfModule() { return m_NumberOfModule; }
+
+ // Get Root input and output objects
+ TFatimaData* GetEventData() {return m_EventData;}//!
+ TFatimaPhysics* GetEventPhysics() {return m_EventPhysics;}//!
+
+ // To be called after a build Physical Event
+ double GetEnergyDeposit();
+ double GetEnergyInDeposit();
+ double GetEnergyOutDeposit();
+ TVector3 GetPositionOfInteraction();
+
+ void Print();
+
+
+ private:
+ ////////////////////////////////////////
+ // Root Input and Output tree classes //
+ ////////////////////////////////////////
+ TFatimaData* m_EventData;//!
+ TFatimaPhysics* m_EventPhysics;//!
+
+
+ private:
+ // Spatial Position of Strip Calculated on basis of detector position
+ int m_NumberOfModule;
+ vector< vector < vector < double > > > m_PositionX;//!
+ vector< vector < vector < double > > > m_PositionY;//!
+ vector< vector < vector < double > > > m_PositionZ;//!
+
+
+ public: // Data Member
+ // Provide Physical Multiplicity
+ Int_t FatimaEventMult;
+
+ // FirstStage
+ vector<double> FatimaLaBr_E;
+
+ // SecondStage
+ vector<double> FatimaNaI_E;
+
+ // Physical Value
+ vector<double> FatimaTotalEnergy;
+ vector<double> FatimaInTotalEnergy;
+ vector<double> FatimaOutTotalEnergy;
+
+ ClassDef(TFatimaPhysics,1)
};
#endif
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 3147c9d51a41c4b561a69dd94663644514331ce7..48390f0ce23eb58135cac78d2af18355e3466d1f 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -31,8 +31,6 @@
#include "../DetectorList.inc"
#include "GaspardTracker.h"
#include "Hyde2Tracker.h"
-#include "TParisPhysics.h"
-#include "Fatima.h"
#include "TAnnularS1Physics.h"
#include "TCATSPhysics.h"
#include "TCharissaPhysics.h"
@@ -41,8 +39,10 @@
#include "TChio_digPhysics.h"
#include "TExlPhysics.h"
#include "TExogamPhysics.h"
+#include "TFatimaPhysics.h"
#include "TLaBr3Physics.h"
#include "TMust2Physics.h"
+#include "TParisPhysics.h"
#include "TPlasticPhysics.h"
#include "TS2Physics.h"
#include "TSSSDPhysics.h"
@@ -262,7 +262,7 @@ void DetectorManager::ReadConfigurationFile(string Path) {
cout << "//////// FATIMA ////////" << endl << endl;
// Instantiate the new array as a VDetector Object
- VDetector* myDetector = new Fatima();
+ VDetector* myDetector = new TFatimaPhysics();
// Read Position of Telescope
ConfigFile.close();
myDetector->ReadConfiguration(Path);
diff --git a/NPSimulation/Fatima/Fatima.cc b/NPSimulation/Fatima/Fatima.cc
old mode 100755
new mode 100644
index d7951824ea28d522dabc383e5c1c19b9aeb2a9fa..f2112cefea3ba4069d41640fd1077e54dc068ffc
--- a/NPSimulation/Fatima/Fatima.cc
+++ b/NPSimulation/Fatima/Fatima.cc
@@ -1,133 +1,488 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/13 *
- * Last update : *
+ * Creation Date : December 2009 *
+ * Last update : December 2014 *
*---------------------------------------------------------------------------*
- * Decription: This class manages different shapes of module for the Fatima *
- * detector. It allows to have Fatima geometries with an *
- * heterogeneous set of modules *
+ * Decription: *
+ * This class describe the Fatima scintillator array *
+ * *
*---------------------------------------------------------------------------*
* Comment: *
* *
- * *
*****************************************************************************/
// C++ headers
-#include <fstream>
#include <sstream>
-#include <string>
#include <cmath>
+#include <limits>
+using namespace std;
-// NPTool headers
+//Geant4
+#include "G4VSolid.hh"
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4UnionSolid.hh"
+#include "G4SubtractionSolid.hh"
+#include "G4SDManager.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4Colour.hh"
+
+// NPS
#include "Fatima.hh"
-#include "FatimaDetector.hh"
+using namespace FATIMA;
-using namespace std;
+#include "CalorimeterScorers.hh"
+#include "MaterialManager.hh"
+// NPL
+#include "NPOptionManager.h"
+#include "RootOutput.h"
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+using namespace CLHEP;
-Fatima::Fatima()
-{
-}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Fatima Specific Method
+Fatima::Fatima(){
+ m_Event = new TFatimaData();
+
+ // Blue
+ m_LaBr3VisAtt = new G4VisAttributes(G4Colour(0, 0.5, 1));
+ // Dark Grey
+ m_PMTVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1));
+ // Grey wireframe
+ m_DetectorCasingVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5,0.2));
+
+ m_LogicalDetector = 0;
+ m_LaBr3Scorer = 0 ;
+}
-Fatima::~Fatima()
-{
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+Fatima::~Fatima(){
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Fatima::AddDetector(G4ThreeVector Pos1, G4ThreeVector Pos2, G4ThreeVector Pos3, G4ThreeVector Pos4){
+ G4ThreeVector Pos=(Pos1+Pos2+Pos3+Pos4)/4.;
+ G4ThreeVector u = Pos1-Pos2;
+ G4ThreeVector v = Pos1-Pos4;
+ u = u.unit(); v = v.unit();
+ G4ThreeVector w = Pos.unit();
+ Pos = Pos + w*Length*0.5;
+
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(new G4RotationMatrix(u,v,w));
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Fatima::AddDetector(G4ThreeVector Pos, double beta_u, double beta_v, double beta_w){
+ double Theta = Pos.theta();
+ double Phi = Pos.phi();
+
+ // vector parallel to one axis of silicon plane
+ G4double ii = cos(Theta / rad) * cos(Phi / rad);
+ G4double jj = cos(Theta / rad) * sin(Phi / rad);
+ G4double kk = -sin(Theta / rad);
+ G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
+
+ G4ThreeVector w = Pos.unit();
+ G4ThreeVector u = w.cross(Y);
+ G4ThreeVector v = w.cross(u);
+ v = v.unit();
+ u = u.unit();
+
+ G4RotationMatrix* r = new G4RotationMatrix(u,v,w);
+ r->rotate(beta_u,u);
+ r->rotate(beta_v,v);
+ r->rotate(beta_w,w);
+
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(r);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of VDetector class
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void Fatima::ReadConfiguration(string Path)
-{
- // open configuration file
- ifstream ConfigFile;
- ConfigFile.open(Path.c_str());
+void Fatima::ReadConfiguration(string Path){
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer, DataBuffer;
+
+ // A,B,C,D are the four corner of the detector
+
+ G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
+ G4ThreeVector A , B , C , D ;
+ G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
+
+ bool ReadingStatus = false;
+
+ bool check_A = false;
+ bool check_C = false;
+ bool check_B = false;
+ bool check_D = false;
+
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+ if (LineBuffer.compare(0, 14, "FatimaDetector") == 0) {
+ G4cout << "///" << G4endl ;
+ G4cout << "Detector found: " << G4endl ;
+ ReadingStatus = true ;
+ }
+
+
+ while (ReadingStatus) {
+ ConfigFile >> DataBuffer;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
+
+ // Position method
+ else if (DataBuffer == "A=") {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax * mm ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay * mm ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az * mm ;
+
+ A = G4ThreeVector(Ax, Ay, Az);
+ G4cout << "Corner A position : " << A << G4endl;
+ }
+ else if (DataBuffer == "B=") {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx * mm ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By * mm ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz * mm ;
+
+ B = G4ThreeVector(Bx, By, Bz);
+ G4cout << "Corner B position : " << B << G4endl;
+ }
- bool bFatimaDetector = false;
+ else if (DataBuffer == "C=") {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx * mm ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy * mm ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz * mm ;
+
+ C = G4ThreeVector(Cx, Cy, Cz);
+ G4cout << "Corner C position : " << C << G4endl;
+ }
+ else if (DataBuffer == "D=") {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx * mm ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy * mm ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz * mm ;
+
+ D = G4ThreeVector(Dx, Dy, Dz);
+ G4cout << "Corner D position : " << D << G4endl;
+ }
- string LineBuffer;
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
- if (LineBuffer.compare(0, 14, "FatimaDetector") == 0 && bFatimaDetector == false) {
- bFatimaDetector = true;
+ // Angle method
+ else if (DataBuffer=="Theta=") {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta * deg;
+ G4cout << "Theta: " << Theta / deg << G4endl;
+ }
+ else if (DataBuffer=="Phi=") {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi * deg;
+ G4cout << "Phi: " << Phi / deg << G4endl;
+ }
+ else if (DataBuffer=="R=") {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R * mm;
+ G4cout << "R: " << R / mm << G4endl;
+ }
+ else if (DataBuffer=="Beta=") {
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w * deg ;
+ G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
+ }
- // instantiate a new "detector" corresponding to the Square elements
- FatimaModule* myDetector = new FatimaDetector();
+ else G4cout << "WARNING: Wrong Token, Fatima: Dector not added" << G4endl;
+
+ // Add The previously define telescope
+ // With position method
+ if ((check_A && check_B && check_C && check_D) &&
+ !(check_Theta && check_Phi && check_R)) {
+ ReadingStatus = false;
+ check_A = false;
+ check_C = false;
+ check_B = false;
+ check_D = false;
+
+ AddDetector(A, B, C, D);
+ }
- // read part of the configuration file corresponding to square elements
- ConfigFile.close();
- myDetector->ReadConfiguration(Path);
- ConfigFile.open(Path.c_str());
+ // With angle method
+ if ((check_Theta && check_Phi && check_R ) &&
+ !(check_A && check_B && check_C && check_D)) {
+ ReadingStatus = false;
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
+
+ R = R + 0.5*Length;
+ G4ThreeVector Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
+ AddDetector(Pos, beta_u, beta_v, beta_w);
+ }
+ }
+ }
+}
- // ms_InterCoord comes from VDetector
- myDetector->SetInterCoordPointer(ms_InterCoord);
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Construct detector and inialise sensitive part.
+// Called After DetecorConstruction::AddDetector Method
+void Fatima::ConstructDetector(G4LogicalVolume* world){
+ unsigned int mysize = m_Pos.size();
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ new G4PVPlacement(G4Transform3D(*m_Rot[i], m_Pos[i]), ConstructDetector(), "FatimaDetector", world, false, i+1);
+ }
+}
- // store FatimaDetector "detector"
- m_Modules.push_back(myDetector);
- }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Fatima::ConstructDetector(){
+ if(!m_LogicalDetector){
+
+ G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
+ G4Material* Alu = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4Material* Lead = MaterialManager::getInstance()->GetMaterialFromLibrary("Pb");
+ G4Material* LaBr3 = MaterialManager::getInstance()->GetMaterialFromLibrary("LaBr3");
+
+ // Mother Volume
+ G4Tubs* solidFatimaDetector =
+ new G4Tubs("Fatima",0, 0.5*FaceFront, 0.5*Length, 0.*deg, 360.*deg);
+ m_LogicalDetector =
+ new G4LogicalVolume(solidFatimaDetector, Vacuum, "Fatima", 0, 0, 0);
+
+ m_LogicalDetector->SetVisAttributes(G4VisAttributes::Invisible);
+
+ // Detector construction
+ // LaBr3
+ G4ThreeVector positionLaBr3 = G4ThreeVector(0, 0, LaBr3_PosZ);
+
+ G4Tubs* solidLaBr3 = new G4Tubs("solidLaBr3", 0., 0.5*LaBr3Face, 0.5*LaBr3Thickness, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3 = new G4LogicalVolume(solidLaBr3, LaBr3, "logicLaBr3", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionLaBr3,
+ logicLaBr3,
+ "Fatima_LaBr3",
+ m_LogicalDetector,
+ false,
+ 0);
+
+ // Set LaBr3 sensible
+ logicLaBr3->SetSensitiveDetector(m_LaBr3Scorer);
+
+ // Visualisation of LaBr3 Strip
+ logicLaBr3->SetVisAttributes(m_LaBr3VisAtt);
+
+ // Aluminium can around LaBr3
+ // LaBr3 Can
+ G4ThreeVector positionLaBr3Can = G4ThreeVector(0, 0, LaBr3Can_PosZ);
+
+ G4Tubs* solidLaBr3Can = new G4Tubs("solidLaBr3Can", 0.5*CanInnerDiameter, 0.5*CanOuterDiameter, 0.5*CanLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3Can = new G4LogicalVolume(solidLaBr3Can, Alu, "logicLaBr3Can", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionLaBr3Can,
+ logicLaBr3Can,
+ "Fatima_LaBr3Can",
+ m_LogicalDetector,
+ false,
+ 0);
+
+ // Visualisation of LaBr3Can
+ logicLaBr3Can->SetVisAttributes(m_DetectorCasingVisAtt);
+
+ // Aluminium window in front of LaBr3
+ // LaBr3 Window
+ G4ThreeVector positionLaBr3Win = G4ThreeVector(0, 0, LaBr3Win_PosZ);
+
+ G4Tubs* solidLaBr3Win = new G4Tubs("solidLaBr3Win", 0.5*WinInnerDiameter, 0.5*WinOuterDiameter, 0.5*WinLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3Win = new G4LogicalVolume(solidLaBr3Win, Alu, "logicLaBr3Win", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionLaBr3Win,
+ logicLaBr3Win,
+ "Fatima_LaBr3Win",
+ m_LogicalDetector,
+ false,
+ 0);
+
+ // Visualisation of LaBr3Win
+ logicLaBr3Win->SetVisAttributes(m_DetectorCasingVisAtt);
+
+ // PMT
+ G4ThreeVector positionPMT = G4ThreeVector(0, 0, PMT_PosZ);
+
+ G4Tubs* solidPMout = new G4Tubs("solidPMOut", 0.5*LaBr3Face, 0.5*PMTFace, 0.5*PMTThickness, 0.*deg, 360.*deg);
+ G4Tubs* solidPMin = new G4Tubs("solidPMIn", 0.5*LaBr3Face-0.1*cm, 0.5*PMTFace-0.5*cm, 0.5*(PMTThickness-2.*cm)-0.1*cm, 0.*deg, 360.*deg);
+ G4RotationMatrix* RotMat=NULL;
+ const G4ThreeVector &Trans= G4ThreeVector(0.,0.,1.*cm);
+ G4SubtractionSolid* solidPMT = new G4SubtractionSolid("solidPMT", solidPMout,solidPMin, RotMat, Trans);
+
+ G4LogicalVolume* logicPMT = new G4LogicalVolume(solidPMT, Alu, "logicPMT", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionPMT,
+ logicPMT,
+ "Fatima_PMT",
+ m_LogicalDetector,
+ false,
+ 0);
+
+ // Visualisation of PMT Strip
+ logicPMT->SetVisAttributes(m_PMTVisAtt);
+
+ // Lead shielding
+ // A
+ G4ThreeVector positionLeadAShield = G4ThreeVector(0, 0, LeadAShield_PosZ);
+ G4Tubs* solidLeadA = new G4Tubs("solidLead", 0.5*LeadAMinR, 0.5*LeadAMaxR, 0.5*LeadALength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLeadAShield = new G4LogicalVolume(solidLeadA, Lead, "logicLeadAShield", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionLeadAShield,
+ logicLeadAShield,
+ "Fatima_LeadAShield",
+ m_LogicalDetector,
+ false,
+ 0);
+ // B
+ G4ThreeVector positionLeadBShield = G4ThreeVector(0, 0, LeadBShield_PosZ);
+ G4Tubs* solidLeadB = new G4Tubs("solidLead", 0.5*LeadBMinR, 0.5*LeadBMaxR, 0.5*LeadBLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLeadBShield = new G4LogicalVolume(solidLeadB, Lead, "logicLeadBShield", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionLeadBShield,
+ logicLeadBShield,
+ "Fatima_LeadBShield",
+ m_LogicalDetector,
+ false,
+ 0);
+
+ // Visualisation of PMT Strip
+ G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(1., 1., 0.));
+ logicLeadAShield->SetVisAttributes(LeadVisAtt);
+ logicLeadBShield->SetVisAttributes(LeadVisAtt);
+ }
+
+ return m_LogicalDetector;
+}
- }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void Fatima::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ pTree->Branch("Fatima", "TFatimaData", &m_Event) ;
+ pTree->SetBranchAddress("Fatima", &m_Event) ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void Fatima::ReadSensitive(const G4Event* event){
+ m_Event->Clear();
-// Construct detector and initialize sensitive part.
-// Called After DetectorConstruction::AddDetector Method
-void Fatima::ConstructDetector(G4LogicalVolume* world)
-{
- // loop on sub-detectors belonging to Fatima
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ConstructDetector(world);
-}
+ ///////////
+ // LaBr3
+ G4THitsMap<G4double*>* LaBr3HitMap;
+ std::map<G4int, G4double**>::iterator LaBr3_itr;
+ G4int LaBr3CollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Fatima_LaBr3Scorer/FatimaLaBr3");
+ LaBr3HitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(LaBr3CollectionID));
+ // Loop on the LaBr3 map
+ for (LaBr3_itr = LaBr3HitMap->GetMap()->begin() ; LaBr3_itr != LaBr3HitMap->GetMap()->end() ; LaBr3_itr++){
-// Connect the FatimaData class to the output TTree
-// of the simulation
-void Fatima::InitializeRootOutput()
-{
- // loop on sub-detectors belonging to Fatima
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeRootOutput();
-}
+ G4double* Info = *(LaBr3_itr->second);
+ double Energy = RandGauss::shoot(Info[0], EnergyResolution);
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int DetectorNbr = (int) Info[2];
-// Initialize all scorers necessary for each detector
-void Fatima::InitializeScorers()
-{
- // loop on sub-detectors belonging to Fatima
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeScorers();
+ m_Event->SetFatimaLaBr3E(DetectorNbr,Energy);
+ m_Event->SetFatimaLaBr3T(DetectorNbr,Time);
+ }
+ }
+ // clear map for next event
+ LaBr3HitMap->clear();
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Fatima::InitializeScorers(){
+ vector<G4int> NestingLevel;
+ NestingLevel.push_back(1);
+ // LaBr3 Associate Scorer
+ bool already_exist = false;
+ m_LaBr3Scorer = CheckScorer("Fatima_LaBr3Scorer",already_exist);
-// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAction
-void Fatima::ReadSensitive(const G4Event* event)
-{
- // Before looping on each sub-detector, clear the static variable
- // ms_InterCoord
- // This is done on the first element of the m_Modules vector.
- // This should be done here since this variable (of type TIneractionCoordinates)
- // deals with multiplicity of events > 1.
- m_Modules[0]->GetInterCoordPointer()->Clear();
-
- // We do the same for the static variable ms_Event
- m_Modules[0]->GetEventPointer()->Clear();
-
- // loop on sub-detectors belonging to Fatima
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ReadSensitive(event);
+ // if the scorer were created previously nothing else need to be made
+ if(already_exist) return;
+
+ G4VPrimitiveScorer* LaBr3Scorer =
+ new CALORIMETERSCORERS::PS_Calorimeter("FatimaLaBr3",NestingLevel);
+ //and register it to the multifunctionnal detector
+ m_LaBr3Scorer->RegisterPrimitive(LaBr3Scorer);
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3Scorer) ;
}
+
diff --git a/NPSimulation/Fatima/Fatima.hh b/NPSimulation/Fatima/Fatima.hh
old mode 100755
new mode 100644
index 545369180849e4a7ce7c4efa8470888c86d59bd1..254fc3ae16d59a64b4346fec5d3ce5492f85c4ca
--- a/NPSimulation/Fatima/Fatima.hh
+++ b/NPSimulation/Fatima/Fatima.hh
@@ -1,5 +1,7 @@
+#ifndef Fatima_h
+#define Fatima_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -8,66 +10,158 @@
/*****************************************************************************
* Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
* *
- * Creation Date : 04/01/13 *
- * Last update : *
+ * Creation Date : December 2009 *
+ * Last update : December 2014 *
*---------------------------------------------------------------------------*
- * Decription: This class manages different shapes of module for the Fatima *
- * detector. It allows to have Fatima geometries with an *
- * heterogeneous set of modules *
+ * Decription: *
+ * This class describe the Fatima scintillator array *
+ * *
*---------------------------------------------------------------------------*
* Comment: *
* *
- * *
*****************************************************************************/
+// C++ header
+#include <string>
+#include <vector>
-#ifndef Fatima_h
-#define Fatima_h 1
+// G4 header defining G4 types
+#include "globals.hh"
-// C++ headers
-#include <vector>
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4VisAttributes.hh"
+#include "G4MultiFunctionalDetector.hh"
-// NPTool header
+// NPSimulation header
#include "VDetector.hh"
-#include "FatimaModule.hh"
+// NPLib
+#include "TFatimaData.h"
using namespace std;
+using namespace CLHEP;
+class Fatima : public VDetector{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+public:
+ Fatima() ;
+ ~Fatima() ;
+
+ ////////////////////////////////////////////////////
+ //////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // To add a Detector
+ // By corner position
+ void AddDetector(G4ThreeVector Pos1,G4ThreeVector Pos2,G4ThreeVector Pos3,G4ThreeVector Pos4);
+ // By Center position
+ void AddDetector(G4ThreeVector Pos, double beta_u=0, double beta_v=0, double beta_w=0);
+
+ // Return a logical volume of the detector
+ G4LogicalVolume* ConstructDetector();
+private: // Guarranty that each volume is created only once
+ G4LogicalVolume* m_LogicalDetector;
-class Fatima : public VDetector
-{
- ////////////////////////////////////////////////////
- /////// Default Constructor and Destructor /////////
- ////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////
+ ///////// Inherite from VDetector class ///////////
+ ////////////////////////////////////////////////////
public:
- Fatima();
- virtual ~Fatima();
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path) ;
+
+ // 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////////////////
+ ////////////////////////////////////////////////////
+private:
+ TFatimaData* m_Event ;
+
+ ////////////////////////////////////////////////////
+ ///////////////// Scorer Related ///////////////////
+ ////////////////////////////////////////////////////
+
+private:
+ // Initialize all Scorer
+ void InitializeScorers() ;
+
+ // Scorer Associate to the Silicon
+ G4MultiFunctionalDetector* m_LaBr3Scorer ;
- ////////////////////////////////////////////////////
- ///////// Inherite from VDetector class ///////////
- ////////////////////////////////////////////////////
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
+private:
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private:
+ // Detector position
+ vector<G4ThreeVector> m_Pos;
+ vector<G4RotationMatrix*> m_Rot;
+
+private:/// Visualisation Attribute:
+ G4VisAttributes* m_LaBr3VisAtt;
+ G4VisAttributes* m_DetectorCasingVisAtt ;
+ G4VisAttributes* m_PMTVisAtt;
+};
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world);
+namespace FATIMA{
+ // Resolution
+ const G4double EnergyResolution = 0.0099; // = 3.5% at .662MeV of Resolution // Unit is MeV/2.35
+ const G4double EnergyThreshold = 100*keV;
+ // Geometry for the mother volume
+ const G4double FaceFront = 7.5*cm;
+ const G4double Length = 26.33*cm;
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput();
+ // LaBr3
+ const G4double LaBr3Face = 3.81*cm;
+ const G4double LaBr3Thickness = 5.08*cm;
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
+ // Al Can
+ const G4double CanOuterDiameter = 4.3*cm;
+ const G4double CanInnerDiameter = 4.15*cm;
+ const G4double CanLength = 4.33*cm;
+ // Al front Window
+ const G4double WinOuterDiameter = 4.15*cm;
+ const G4double WinInnerDiameter = 0*cm;
+ const G4double WinLength = 0.1*cm; //
-public:
- // Initialize all scorers necessary for each detector
- void InitializeScorers();
+ // PMT
+ const G4double PMTFace = FaceFront;
+ const G4double PMTThickness = 22.*cm; // for detector
+
+ // Lead tube
+ const G4double LeadAMinR = CanOuterDiameter;
+ const G4double LeadAMaxR = CanOuterDiameter+1*cm;
+ const G4double LeadALength = 4.33*cm;
+
+ const G4double LeadBMinR = CanOuterDiameter;
+ const G4double LeadBMaxR = FaceFront;
+ const G4double LeadBLength = 1.*cm;
+
+ // Position
+ const G4double LaBr3_PosZ = -Length*0.5 + 0.5*LaBr3Thickness + 0.1*cm;
+ const G4double LaBr3Can_PosZ = -Length*0.5 + 0.5*CanLength;
+ const G4double LaBr3Win_PosZ = -Length*0.5 + 0.5*WinLength;
+
+ const G4double PMT_PosZ = -Length*0.5 + (Length-PMTThickness) + 0.5*PMTThickness;
+
+ const G4double LeadAShield_PosZ = -Length*0.5 + 0.5*LeadALength;
+ const G4double LeadBShield_PosZ = -Length*0.5 + LeadALength - 0.5*LeadBLength ;
+
+}
-private:
- vector<FatimaModule*> m_Modules;
-};
#endif
diff --git a/NPSimulation/Fatima/FatimaDetector.cc b/NPSimulation/Fatima/FatimaDetector.cc
deleted file mode 100755
index 98bc7d6a8d83f15e89dc56aa951e4da06193bfe2..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaDetector.cc
+++ /dev/null
@@ -1,1237 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/14 *
- *---------------------------------------------------------------------------*
- * Decription: Define a detector module for the Fatima detector. *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-// C++ headers
-#include <sstream>
-#include <string>
-#include <cmath>
-
-// G4 Geometry headers
-#include "G4Trd.hh"
-#include "G4Box.hh"
-#include "G4Tubs.hh"
-#include "G4Trap.hh"
-
-#include "G4SubtractionSolid.hh"
-
-// G4 various headers
-#include "G4MaterialTable.hh"
-#include "G4Element.hh"
-#include "G4ElementTable.hh"
-#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
-#include "G4RotationMatrix.hh"
-#include "G4Transform3D.hh"
-#include "G4PVPlacement.hh"
-#include "G4PVDivision.hh"
-
-// G4 sensitive
-#include "G4SDManager.hh"
-#include "G4MultiFunctionalDetector.hh"
-
-// NPTool headers
-#include "FatimaDetector.hh"
-#include "ObsoleteGeneralScorers.hh"
-#include "FatimaScorers.hh"
-#include "RootOutput.h"
-
-// CLHEP
-#include "CLHEP/Random/RandGauss.h"
-
-using namespace std;
-using namespace CLHEP;
-using namespace FATIMADETECTOR;
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-FatimaDetector::FatimaDetector()
-{
- ms_InterCoord = 0;
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-FatimaDetector::~FatimaDetector()
-{
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void FatimaDetector::AddModule(G4ThreeVector X1_Y1,
- G4ThreeVector X128_Y1,
- G4ThreeVector X1_Y128,
- G4ThreeVector X128_Y128)
-{
- m_DefinitionType.push_back(true);
-
- m_X1_Y1.push_back(X1_Y1);
- m_X128_Y1.push_back(X128_Y1);
- m_X1_Y128.push_back(X1_Y128);
- m_X128_Y128.push_back(X128_Y128);
-
- m_R.push_back(0);
- m_Theta.push_back(0);
- m_Phi.push_back(0);
- m_beta_u.push_back(0);
- m_beta_v.push_back(0);
- m_beta_w.push_back(0);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void FatimaDetector::AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w)
-{
- G4ThreeVector empty = G4ThreeVector(0, 0, 0);
-
- m_DefinitionType.push_back(false);
-
- m_X1_Y1.push_back(empty);
- m_X128_Y1.push_back(empty);
- m_X1_Y128.push_back(empty);
- m_X128_Y128.push_back(empty);
-
- m_R.push_back(R);
- m_Theta.push_back(Theta);
- m_Phi.push_back(Phi);
- m_beta_u.push_back(beta_u);
- m_beta_v.push_back(beta_v);
- m_beta_w.push_back(beta_w);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void FatimaDetector::VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world)
-{
- G4double NbrTelescopes = DetecNumber;
- G4String DetectorNumber;
- ostringstream Number;
- Number << NbrTelescopes;
- DetectorNumber = Number.str();
-
- /////////////////////////////////////////////////////////////////
- /////////////////Element Definition ///////////////////////////
- ////////////////////////////////////////////////////////////////
- G4String symbol, name;
- G4double density = 0. , a = 0, z = 0;
- G4int ncomponents = 0;
- G4int nel = 0, natoms = 0;
-
- ////////////////////////////////////////////////////////////////
- /////////////////Material Definition ///////////////////////////
- ////////////////////////////////////////////////////////////////
- a=137.327*g/mole;
- G4Element* Ba = new G4Element(name="Barium",symbol="Ba",z=56.,a);
- a=18.9984032*g/mole;
- G4Element* F = new G4Element(name="Fluorine",symbol="F",z=9.,a);
- a=22.99*g/mole;
- G4Element* Na = new G4Element(name="Sodium",symbol="Na",z=11.,a);
- a=79.904*g/mole;
- G4Element* Br = new G4Element(name="Bromine",symbol="Br",z=35.,a);
- a=126.90477*g/mole;
- G4Element* I = new G4Element(name="Iodine",symbol="I",z=53.,a);
- a=132.90545*g/mole;
- G4Element* Cs = new G4Element(name="Cesium",symbol="Cs",z=55.,a);
- a=138.9055*g/mole;
- G4Element* La = new G4Element(name="Lanthanum",symbol="La",z=57.,a);
-
- // Vacuum
- G4Element* N = new G4Element("Nitrogen" , symbol = "N" , z = 7 , a = 14.01 * g / mole);
- G4Element* O = new G4Element("Oxigen" , symbol = "O" , z = 8 , a = 16.00 * g / mole);
-
- density = 0.000000001 * mg / cm3;
- G4Material* Vacuum = new G4Material("Vacuum", density, ncomponents = 2);
- Vacuum->AddElement(N, .7);
- Vacuum->AddElement(O, .3);
-
- // Germanium
- density = 5.32*g/cm3;
- G4Material* Ge;
- Ge = new G4Material(name="Germanium",32,72.59*g/mole, density);
-
- // Aluminum
- density = 2.7*g/cm3;
- G4Material* Alu;
- Alu = new G4Material(name="Aluminum",13,26.98154*g/mole, density);
-
- // lead
- density =11.35*g/cm3;
- G4Material* Lead;
- Lead = new G4Material(name="Lead",82,207.2*g/mole, density);
-
- // NaI
- density = 3.67*g/cm3, nel = 2;
- G4Material* NaI = new G4Material(name="NaI",density,nel);
- NaI->AddElement(Na, natoms = 1);
- NaI->AddElement(I, natoms = 1);
-
- // CsI
- density = 4.51*g/cm3, nel = 2;
- G4Material* CsI = new G4Material(name="CsI", density, nel);
- CsI->AddElement(Cs, natoms = 1);
- CsI->AddElement(I, natoms = 1);
-
- // LaBr3
- density = 5.29*g/cm3, nel = 2;
- G4Material* LaBr3 = new G4Material(name="LaBr3",density,nel);
- LaBr3->AddElement(La, natoms = 1);
- LaBr3->AddElement(Br, natoms = 3);
-
- // BaF2
- density = 4.89*g/cm3, nel = 2;
- G4Material* BaF2 = new G4Material(name="BaF2", density, nel);
- BaF2->AddElement(Ba, natoms = 1);
- BaF2->AddElement(F, natoms = 2);
-
- ////////////////////////////////////////////////////////////////
- ////////////// Starting Volume Definition //////////////////////
- ////////////////////////////////////////////////////////////////
- // Little trick to avoid warning in compilation: Use a PVPlacement "buffer".
- // If don't you will have a Warning unused variable 'myPVP'
- G4PVPlacement* PVPBuffer;
- G4String Name = "FatimaDetector" + DetectorNumber ;
-
- // Mother Volume
- //G4Box* solidFatimaDetector = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
- G4Tubs* solidFatimaDetector = new G4Tubs(Name,0, 0.5*FaceFront, 0.5*Length, 0.*deg, 360.*deg);
- G4LogicalVolume* logicFatimaDetector = new G4LogicalVolume(solidFatimaDetector, Vacuum, Name, 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos) ,
- logicFatimaDetector ,
- Name ,
- world ,
- false ,
- 0);
-
- logicFatimaDetector->SetVisAttributes(G4VisAttributes::Invisible);
- //logicFatimaDetector->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
- //if (m_non_sensitive_part_visiualisation) logicFatimaDetector->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
-
- // Detector construction
- // LaBr3
- G4ThreeVector positionLaBr3Stage = G4ThreeVector(0, 0, LaBr3Stage_PosZ);
-
- //G4Box* solidLaBr3Stage = new G4Box("solidLaBr3Stage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*LaBr3Thickness);
- G4Tubs* solidLaBr3Stage = new G4Tubs("solidLaBr3Stage", 0., 0.5*LaBr3Face, 0.5*LaBr3Thickness, 0.*deg, 360.*deg);
- G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, LaBr3, "logicLaBr3Stage", 0, 0, 0);
- //G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, Ge, "logicLaBr3Stage", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionLaBr3Stage,
- logicLaBr3Stage,
- Name + "_LaBr3Stage",
- logicFatimaDetector,
- false,
- 0);
-
- // Set LaBr3 sensible
- logicLaBr3Stage->SetSensitiveDetector(m_LaBr3StageScorer);
-
- // Visualisation of LaBr3Stage Strip
- G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(1., 0., 0.));
- logicLaBr3Stage->SetVisAttributes(LaBr3VisAtt);
-
-
- // Aluminium can around LaBr3
- // LaBr3 Can
- G4ThreeVector positionLaBr3Can = G4ThreeVector(0, 0, LaBr3Can_PosZ);
-
- G4Tubs* solidLaBr3Can = new G4Tubs("solidLaBr3Can", 0.5*CanInnerDiameter, 0.5*CanOuterDiameter, 0.5*CanLength, 0.*deg, 360.*deg);
- G4LogicalVolume* logicLaBr3Can = new G4LogicalVolume(solidLaBr3Can, Alu, "logicLaBr3Can", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionLaBr3Can,
- logicLaBr3Can,
- Name + "_LaBr3Can",
- logicFatimaDetector,
- false,
- 0);
-
-
- // Visualisation of LaBr3Can
- G4VisAttributes* CanVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
- logicLaBr3Can->SetVisAttributes(CanVisAtt);
-
- // Aluminium window in front of LaBr3
- // LaBr3 Window
- G4ThreeVector positionLaBr3Win = G4ThreeVector(0, 0, LaBr3Win_PosZ);
-
- G4Tubs* solidLaBr3Win = new G4Tubs("solidLaBr3Win", 0.5*WinInnerDiameter, 0.5*WinOuterDiameter, 0.5*WinLength, 0.*deg, 360.*deg);
- G4LogicalVolume* logicLaBr3Win = new G4LogicalVolume(solidLaBr3Win, Alu, "logicLaBr3Win", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionLaBr3Win,
- logicLaBr3Win,
- Name + "_LaBr3Win",
- logicFatimaDetector,
- false,
- 0);
- // Visualisation of LaBr3Win
- logicLaBr3Win->SetVisAttributes(CanVisAtt);
-
-
- // CsI stage replaced by Aluminum housing (TODO: replace CsI by Alu in variable names to avoid confusion)
- G4ThreeVector positionCsIStage = G4ThreeVector(0, 0, CsIStage_PosZ);
-
- //G4Box* solidCsIStage = new G4Box("solidCsIStage", 0.5*CsIFace, 0.5*CsIFace, 0.5*CsIThickness);
-
- G4Tubs* solidPMout = new G4Tubs("solidPMOut", 0.5*LaBr3Face, 0.5*CsIFace, 0.5*CsIThickness, 0.*deg, 360.*deg);
- G4Tubs* solidPMin = new G4Tubs("solidPMIn", 0.5*LaBr3Face-0.1*cm, 0.5*CsIFace-0.5*cm, 0.5*(CsIThickness-2.*cm)-0.1*cm, 0.*deg, 360.*deg);
- G4RotationMatrix* RotMat=NULL;
- const G4ThreeVector &Trans= G4ThreeVector(0.,0.,1.*cm);
- G4SubtractionSolid* solidCsIStage = new G4SubtractionSolid("solidCsIStage", solidPMout,solidPMin, RotMat, Trans);
-
- //G4Tubs* solidCsIStage = new G4Tubs("solidCsIStage", 0.5*LaBr3Face, 0.5*CsIFace, 0.5*CsIThickness, 0.*deg, 360.*deg);
- //G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, CsI, "logicCsIStage", 0, 0, 0);
- G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, Alu, "logicCsIStage", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionCsIStage,
- logicCsIStage,
- Name + "_CsIStage",
- logicFatimaDetector,
- false,
- 0);
-
-
- // Visualisation of CsIStage Strip
- G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
- logicCsIStage->SetVisAttributes(CsIVisAtt);
-
-
-
- // Lead shielding
- // A
- G4ThreeVector positionLeadAShield = G4ThreeVector(0, 0, LeadAShield_PosZ);
- G4Tubs* solidLeadA = new G4Tubs("solidLead", 0.5*LeadAMinR, 0.5*LeadAMaxR, 0.5*LeadALength, 0.*deg, 360.*deg);
- G4LogicalVolume* logicLeadAShield = new G4LogicalVolume(solidLeadA, Lead, "logicLeadAShield", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionLeadAShield,
- logicLeadAShield,
- Name + "_LeadAShield",
- logicFatimaDetector,
- false,
- 0);
- // B
- G4ThreeVector positionLeadBShield = G4ThreeVector(0, 0, LeadBShield_PosZ);
- G4Tubs* solidLeadB = new G4Tubs("solidLead", 0.5*LeadBMinR, 0.5*LeadBMaxR, 0.5*LeadBLength, 0.*deg, 360.*deg);
- G4LogicalVolume* logicLeadBShield = new G4LogicalVolume(solidLeadB, Lead, "logicLeadBShield", 0, 0, 0);
-
- PVPBuffer = new G4PVPlacement(0,
- positionLeadBShield,
- logicLeadBShield,
- Name + "_LeadBShield",
- logicFatimaDetector,
- false,
- 0);
-
-
- // Visualisation of CsIStage Strip
- G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(1., 1., 0.));
- logicLeadAShield->SetVisAttributes(LeadVisAtt);
- logicLeadBShield->SetVisAttributes(LeadVisAtt);
-
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Virtual Method of VDetector class
-
-// Read stream at Configfile to pick-up parameters of detector (Position,...)
-// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void FatimaDetector::ReadConfiguration(string Path)
-{
- ifstream ConfigFile;
- ConfigFile.open(Path.c_str());
- string LineBuffer, DataBuffer;
-
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
- G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
- G4ThreeVector A , B , C , D ;
- G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
-
- bool ReadingStatus = false;
-
- bool check_A = false;
- bool check_C = false;
- bool check_B = false;
- bool check_D = false;
-
- bool check_Theta = false;
- bool check_Phi = false;
- bool check_R = false;
- bool check_beta = false;
-
- bool checkVis = false;
-
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
- if (LineBuffer.compare(0, 14, "FatimaDetector") == 0) {
- G4cout << "///" << G4endl ;
- G4cout << "Detector element found: " << G4endl ;
- ReadingStatus = true ;
- }
-
- while (ReadingStatus) {
- ConfigFile >> DataBuffer;
- // Comment Line
- if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
-
- // Position method
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- Ax = atof(DataBuffer.c_str()) ;
- Ax = Ax * mm ;
- ConfigFile >> DataBuffer ;
- Ay = atof(DataBuffer.c_str()) ;
- Ay = Ay * mm ;
- ConfigFile >> DataBuffer ;
- Az = atof(DataBuffer.c_str()) ;
- Az = Az * mm ;
-
- A = G4ThreeVector(Ax, Ay, Az);
- cout << "X1 Y1 corner position : " << A << endl;
- }
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- Bx = atof(DataBuffer.c_str()) ;
- Bx = Bx * mm ;
- ConfigFile >> DataBuffer ;
- By = atof(DataBuffer.c_str()) ;
- By = By * mm ;
- ConfigFile >> DataBuffer ;
- Bz = atof(DataBuffer.c_str()) ;
- Bz = Bz * mm ;
-
- B = G4ThreeVector(Bx, By, Bz);
- cout << "X128 Y1 corner position : " << B << endl;
- }
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- Cx = atof(DataBuffer.c_str()) ;
- Cx = Cx * mm ;
- ConfigFile >> DataBuffer ;
- Cy = atof(DataBuffer.c_str()) ;
- Cy = Cy * mm ;
- ConfigFile >> DataBuffer ;
- Cz = atof(DataBuffer.c_str()) ;
- Cz = Cz * mm ;
-
- C = G4ThreeVector(Cx, Cy, Cz);
- cout << "X1 Y128 corner position : " << C << endl;
- }
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- Dx = atof(DataBuffer.c_str()) ;
- Dx = Dx * mm ;
- ConfigFile >> DataBuffer ;
- Dy = atof(DataBuffer.c_str()) ;
- Dy = Dy * mm ;
- ConfigFile >> DataBuffer ;
- Dz = atof(DataBuffer.c_str()) ;
- Dz = Dz * mm ;
-
- D = G4ThreeVector(Dx, Dy, Dz);
- cout << "X128 Y128 corner position : " << D << endl;
- }
-
- // Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- Theta = Theta * deg;
- cout << "Theta: " << Theta / deg << endl;
- }
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- Phi = Phi * deg;
- cout << "Phi: " << Phi / deg << endl;
- }
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
- R = R * mm;
- cout << "R: " << R / mm << endl;
- }
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
- check_beta = true;
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- beta_u = beta_u * deg ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- beta_v = beta_v * deg ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- beta_w = beta_w * deg ;
- G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
- }
-
- else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
- checkVis = true ;
- ConfigFile >> DataBuffer;
- if (DataBuffer.compare(0, 3, "all") == 0) m_non_sensitive_part_visiualisation = true;
- }
-
- else G4cout << "WARNING: Wrong Token, FatimaDetector: Detector Element not added" << G4endl;
-
- // Add The previously define telescope
- // With position method
- if ((check_A && check_B && check_C && check_D && checkVis) &&
- !(check_Theta && check_Phi && check_R)) {
- ReadingStatus = false;
- check_A = false;
- check_C = false;
- check_B = false;
- check_D = false;
- checkVis = false;
-
- AddModule(A, B, C, D);
- }
-
- // With angle method
- if ((check_Theta && check_Phi && check_R && checkVis) &&
- !(check_A && check_B && check_C && check_D)) {
- ReadingStatus = false;
- check_Theta = false;
- check_Phi = false;
- check_R = false;
- check_beta = false;
- checkVis = false;
-
- AddModule(R, Theta, Phi, beta_u, beta_v, beta_w);
- }
- }
- }
-}
-
-
-
-// Construct detector and inialise sensitive part.
-// Called After DetectorConstruction::AddDetector Method
-void FatimaDetector::ConstructDetector(G4LogicalVolume* world)
-{
- G4RotationMatrix* MMrot = NULL ;
- G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
-
- G4int NumberOfTelescope = m_DefinitionType.size() ;
-
- for (G4int i = 0; i < NumberOfTelescope; i++) {
- // By Point
- if (m_DefinitionType[i]) {
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- MMu = m_X128_Y1[i] - m_X1_Y1[i];
- MMu = MMu.unit();
-
- MMv = m_X1_Y128[i] - m_X1_Y1[i];
- MMv = MMv.unit();
-
- // G4double MMscal = MMu.dot(MMv);
-
- MMw = MMu.cross(MMv);
-// if (MMw.z() > 0) MMw = MMv.cross(MMu) ;
- MMw = MMw.unit();
-
- MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4;
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + MMCenter;
- }
-
- // By Angle
- else {
- G4double Theta = m_Theta[i] ;
- G4double Phi = m_Phi[i] ;
-
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- // Phi is angle between X axis and projection in (X,Y) plan
- // Theta is angle between position vector and z axis
- G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
- G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
- G4double wZ = m_R[i] * cos(Theta / rad);
- MMw = G4ThreeVector(wX, wY, wZ);
-
- // vector corresponding to the center of the module
- G4ThreeVector CT = MMw;
-
- // vector parallel to one axis of silicon plane
- G4double ii = cos(Theta / rad) * cos(Phi / rad);
- G4double jj = cos(Theta / rad) * sin(Phi / rad);
- G4double kk = -sin(Theta / rad);
- G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
-
- MMw = MMw.unit();
- MMu = MMw.cross(Y);
- MMv = MMw.cross(MMu);
- MMv = MMv.unit();
- MMu = MMu.unit();
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // Telescope is rotate of Beta angle around MMv axis.
- MMrot->rotate(m_beta_u[i], MMu);
- MMrot->rotate(m_beta_v[i], MMv);
- MMrot->rotate(m_beta_w[i], MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + CT ;
-
- G4cout << "%%%%%%%%%%%%%%" << G4endl;
- G4cout << MMpos << " " << CT << G4endl;
- }
-
- VolumeMaker(i + 1, MMpos, MMrot, world);
- }
-
- delete MMrot ;
-}
-
-
-
-// Connect the FatimaData class to the output TTree
-// of the simulation
-void FatimaDetector::InitializeRootOutput()
-{
-}
-
-
-
-// Set the TinteractionCoordinates object from VDetector to the present class
-void FatimaDetector::SetInterCoordPointer(TInteractionCoordinates* interCoord)
-{
- ms_InterCoord = interCoord;
-}
-
-
-
-// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAvtion
-void FatimaDetector::ReadSensitive(const G4Event* event)
-{
- //////////////////////////////////////////////////////////////////////////////////////
- //////////////////////// Used to Read Event Map of detector //////////////////////////
- //////////////////////////////////////////////////////////////////////////////////////
-
- momentum = event->GetPrimaryVertex()->GetPrimary()->GetMomentum();
- G4double EGamma;
- EGamma = momentum.getR(); // for photon E=p
- //G4double EGammaMin = EGamma-4*ResoFirstStage;
- //G4double EGammaMax = EGamma+4*ResoFirstStage;
-
-
- // First Stage (LaBr3)
- std::map<G4int, G4int*>::iterator DetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CrystalNumber_itr;
- std::map<G4int, G4double*>::iterator Energy_itr;
- std::map<G4int, G4double*>::iterator Time_itr;
- //std::map<G4int, G4double*>::iterator X_itr;
- //std::map<G4int, G4double*>::iterator Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_X_itr;
- //std::map<G4int, G4double*>::iterator Pos_Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_Z_itr;
- //std::map<G4int, G4double*>::iterator Ang_Theta_itr;
- //std::map<G4int, G4double*>::iterator Ang_Phi_itr;
-
- G4THitsMap<G4int>* DetectorNumberHitMap;
- G4THitsMap<G4int>* CrystalNumberHitMap;
- G4THitsMap<G4double>* EnergyHitMap;
- G4THitsMap<G4double>* TimeHitMap;
- // G4THitsMap<G4double>* XHitMap;
- // G4THitsMap<G4double>* YHitMap;
- //G4THitsMap<G4double>* PosXHitMap;
- //G4THitsMap<G4double>* PosYHitMap;
- //G4THitsMap<G4double>* PosZHitMap;
- //G4THitsMap<G4double>* AngThetaHitMap;
- //G4THitsMap<G4double>* AngPhiHitMap;
-
- // NULL pointer are given to avoid warning at compilation
- // Second Stage (CsI)
- /*
- std::map<G4int, G4int*>::iterator CsIDetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CsICrystalNumber_itr; // added by Marc
- std::map<G4int, G4double*>::iterator CsIStageEnergy_itr ;
- G4THitsMap<G4int>* CsIDetectorNumberHitMap = NULL ;
- G4THitsMap<G4int>* CsICrystalNumberHitMap = NULL ; // added by Marc
- G4THitsMap<G4double>* CsIStageEnergyHitMap = NULL ;
-
- */
-
- // Read the Scorer associate to the LaBr
- //Detector Number
- G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/DetectorNumber") ;
- DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ;
- DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ;
-
- //Crystal Number
- G4int CrystDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/CrystalNumber") ;
- CrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CrystDetCollectionID)) ;
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin() ;
-
- //Energy
- G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripEnergy") ;
- EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
- Energy_itr = EnergyHitMap->GetMap()->begin() ;
-
- //Time of Flight
- G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripTime") ;
- TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ;
- Time_itr = TimeHitMap->GetMap()->begin() ;
-
- /*
- //Strip Number X
- G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripIDFront") ;
- XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ;
- X_itr = XHitMap->GetMap()->begin() ;
-
- //Strip Number Y
- G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripIDBack");
- YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ;
- Y_itr = YHitMap->GetMap()->begin() ;
-
-
- //Interaction Coordinate X
- G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordX") ;
- PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ;
- Pos_X_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Y
- G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordY") ;
- PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ;
- Pos_Y_itr = PosYHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Z
- G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordZ") ;
- PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ;
- Pos_Z_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Theta
- G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordAngTheta") ;
- AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ;
- Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Phi
- G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordAngPhi") ;
- AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ;
- Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ;
- */
-
- /*
- // Read the Scorer associate to the SecondStage
- //CsI Detector Number
- G4int CsIDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsIDetectorNumber") ;
- CsIDetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsIDetCollectionID)) ;
- CsIDetectorNumber_itr = CsIDetectorNumberHitMap->GetMap()->begin() ;
-
- //CsI Crystal Number // added by Marc
- G4int CsICryCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsICrystalNumber") ; // added by Marc
- CsICrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICryCollectionID)) ; // added by Anna
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin() ; // added by Anna
-
- //CsI Energy
- G4int CsIStageEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsIStageEnergy") ;
- CsIStageEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIStageEnergyCollectionID)) ;
- CsIStageEnergy_itr = CsIStageEnergyHitMap->GetMap()->begin() ;
- */
-
- // Check the size of different map
- G4int sizeN = DetectorNumberHitMap->entries(); // number of objects hit by trackID=1 (can be the same object hit several time)
- //G4int sizeC = CrystalNumberHitMap->entries();
- G4int sizeE = EnergyHitMap->entries(); // = number of steps with edep non null
- G4int sizeT = TimeHitMap->entries();
- //G4int sizeX = PosXHitMap->entries();
- //G4int sizeY = PosYHitMap->entries();
- //G4int sizeX = XHitMap->entries();
- //G4int sizeY = YHitMap->entries();
-
- /*
- G4int sizeNCsI= CsIDetectorNumberHitMap->entries();
- //G4int sizeCCsI= CsICrystalNumberHitMap->entries();
- //G4int sizeECsI= CsIStageEnergyHitMap->entries();
- */
-
- //sizeC *= 1; // remove warning at compilation added by Marc
- //sizeECsI *= 1; // remove warning at compilation added by Marc
- //G4cout <<"SizeN=" << sizeN << endl;
- //G4cout <<"SizeC=" << sizeC << endl;
- //G4cout <<"SizeN CsI =" << sizeNCsI << endl;
- //G4cout <<"SizeE CsI =" << sizeECsI << endl;
-
- //DetectorNumberHitMap->PrintAllHits();
-
-
- if (sizeE != sizeT) {
- G4cout << "No match size FATIMA Event Map: sE:"
- << sizeE << " sT:" << sizeT << endl ;
-
- // if (sizeE != sizeX) {
- //G4cout << "No match size FATIMA Event Map: sE:"
- //<< sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << endl ;
- return;
- }
-
-
- //G4cout <<"SizeN=" << sizeN << G4endl;
-
-
- //if at least 1 detector is hit:
- if(sizeN>0)
- {
-
- // Deal with trackID=1:
- G4int N_first= *(DetectorNumber_itr->second); // ID of first det hit
- G4int NTrackID = DetectorNumber_itr->first - N_first; // first trackID dealt with (not always =1)
- G4double E = *(Energy_itr->second);
- G4double T = *(Time_itr->second);
- G4int NCryst= *(CrystalNumber_itr->second);
-
- NCryst *= 1; //remove warning at compilation // added by Marc
- /*
- G4cout <<"NTrackID=" << NTrackID << G4endl;
- G4cout <<"N_first=" << N_first << G4endl;
- G4cout <<"CrystalNumber_first=" << NCryst << G4endl;
- G4cout <<"Energy first=" << E << G4endl;
- G4cout <<"Time first =" << T << G4endl;
- cout<<"*******"<<endl;
-
- //added by Nicolas on the model of gaspard scorers 8/7/11
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin();
-
- for (G4int l = 0 ; l < sizeC ; l++) {
- G4double NCryst = *(CrystalNumber_itr->second);
- G4int CTrackID = CrystalNumber_itr->first - N_first - NCryst;
- G4cout << N_first << " " << NTrackID << " " << CTrackID << " " << NCryst << G4endl;
- if (CTrackID == NTrackID) {
- G4cout << "****************** NCryst " << NCryst << G4endl;
- //ms_Event->SetGPDTrkFirstStageFrontEEnergy(RandGauss::shoot(E, ResoFirstStage));
- //ms_Event->SetGPDTrkFirstStageBackEEnergy(RandGauss::shoot(E, ResoFirstStage));
- }
- CrystalNumber_itr++;
- }
- //rewind:
- for (G4int l = 0 ; l < sizeC ; l++) {
- CrystalNumber_itr--;
- }
- */
-
- //if there is more than 1 interaction in crystals (1 interaction = 1 hit)
- if(sizeN>1)
- {
-
- //At clusters'level
- for (G4int l = 1; l < sizeN ; l++) { // loop on all the other tracks
-
- Energy_itr++;
- Time_itr++;
- DetectorNumber_itr++;
- CrystalNumber_itr++;
-
- G4int N= *(DetectorNumber_itr->second); // ID of Detector hit
- G4int Nc= *(CrystalNumber_itr->second); // ID of hit crystal (always =0 for Detector)
-
- NTrackID = DetectorNumber_itr->first - N; // ID of the track
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"N=" << N << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << DetectorNumber_itr->first << G4endl;
- //G4cout <<"NTrackID=" << NTrackID << G4endl;
-
- if(N==N_first)
- {
- // At crystals'level:
- if(Nc==NCryst) // if we are still in the same crystal than the first hit crystal
- {
- E += *(Energy_itr->second);
-
- }else // we fill the tree for the first detector hit and move to the next detector hit
- {
- if(E!=0)
- {
- //G4cout <<" -1- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Marc
- // Fill Energy
- // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
- NCryst=Nc; // continue with the next hit crystal
- E=*(Energy_itr->second);
- }
-
-
- }else // if it's a new detector, one fills the tree with first hit and inspect this new cluster
- {
-
- if(E!=0)
- {
- //G4cout <<" -2- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first );
- ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first );
- ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Anna
- // Fill Energy
- //ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree with energy deposited in first detector
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
- N_first=N;
- NCryst=Nc; // and continue with the new hit crystal in this new cluster
- E=*(Energy_itr->second);
-
- }
-
-
- //G4cout <<"Energy=" << E << G4endl;
- //G4cout <<"Time =" << T << G4endl;
-
- // Always fill the tree at the end of the loop:
- if(l==(sizeN-1) && E!=0)
- {
- //G4cout <<"-3- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Anna
- //cout<<NTrackID<<" filled at the end "<<NCryst<<endl;
- // Fill Energy
- // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
-
- }
-
- }else
- {
- // Fill the tree if sizeN=1:
- if(E!=0)
- {
- //G4cout <<" -4- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
- ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst);
- //cout<<NTrackID<<" filled with sizeN=1 "<<NCryst<<endl;
- // Fill Energy
- // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
- }
-
-
-
- }
-
-
- ///////// For CsI stage:
- //EGammaMin = EGamma-4*ResoSecondStage;
- //EGammaMax = EGamma+4*ResoSecondStage;
- /* if(sizeNCsI>0)
- {
- // Deal with trackID=1:
- G4int NCsI_first= *(CsIDetectorNumber_itr->second); // ID of first det hit
- G4int NCsITrackID = CsIDetectorNumber_itr->first - NCsI_first; // first trackID dealt with (not always =1)
- G4double E_CsI=*(CsIStageEnergy_itr->second); // Energy second stage
- G4int NCsICryst= *(CsICrystalNumber_itr->second); // added by Marc
-
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
- //G4cout <<"NCsI_first=" << NCsI_first << G4endl;
- //G4cout <<"CsI Energy first=" << E_CsI << G4endl;
-
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin(); // added by Marc
-
- //for (G4int l = 0 ; l < sizeCCsI ; l++) {
- // G4double NCsICryst = *(CsICrystalNumber_itr->second);
- // G4int CCsITrackID = CsICrystalNumber_itr->first - NCsI_first - NCsICryst;
- // G4cout << NCsI_first << " " << NCsITrackID << " " << CCsITrackID << " " << NCsICryst << G4endl;
- // if (CCsITrackID == NCsITrackID) {
- // G4cout << "****************** NCsICryst " << NCsICryst << G4endl;
- // }
- // CsICrystalNumber_itr++;
- // }
- ////rewind:
- //for (G4int l = 0 ; l < sizeCCsI ; l++) {
- // CsICrystalNumber_itr--;
- //}
-
-
- if(sizeNCsI>1)
- {
- for (G4int l = 1; l < sizeNCsI ; l++) { // loop on all the other tracks
-
- CsIStageEnergy_itr++;
- CsIDetectorNumber_itr++;
- CsICrystalNumber_itr++;
-
-
- G4int NCsI= *(CsIDetectorNumber_itr->second); // ID of Detector hit
- G4int NcCsI= *(CsICrystalNumber_itr->second); // ID of crystal hit (always =0 for detector)
- NCsITrackID = CsIDetectorNumber_itr->first - NCsI; // ID of the track
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"NCsI=" << NCsI << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << CsIDetectorNumber_itr->first << G4endl;
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
-
- // At detectors'level:
- if(NCsI==NCsI_first)
- {
- // At crystals'level:
- if(NcCsI==NCsICryst) // if we are still in the same crystal than the first hit CsI crystal
- {
-
- E_CsI += *(CsIStageEnergy_itr->second);
-
- }else // we fill the tree for the first detector hit and move to the next detector hit
- {
- if(E_CsI!=0)
- {
- // Fill detector number
- ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- //ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
- }
-
- NCsICryst=NcCsI;
- E_CsI=*(CsIStageEnergy_itr->second);
- }
-
- }else // if it's a new cluster, one fills the tree with first hit and inspect the new cluster
- {
- if(E_CsI!=0)
- {
- //G4cout <<" -2- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first );
- ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first );
- ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst); // added by Marc
- // Fill Energy
- //ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoFirstStage);
- ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree with energy deposited in first detector
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- }
-
- NCsI_first=NCsI; // continue with the new hit cluster
- NCsICryst=NcCsI; // and continue with the new hit CsI crystal in this new cluster
- E_CsI=*(CsIStageEnergy_itr->second);
-
- }
-
- // Always fill the tree at the end of the loop:
-
- if(l==(sizeNCsI-1) && E_CsI!=0)
- {
- //G4cout <<" -3- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- // ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- }
-
- }
-
- }else
- {
- // Fill the tree if sizeN=1:
- if(E_CsI!=0)
- {
- //////G4cout <<" -4- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
- ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- //ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
- }
- }
-
- }
-
- */
-
-
- // clear map for next event
- DetectorNumberHitMap -> clear();
- CrystalNumberHitMap -> clear(); // added by Marc
- EnergyHitMap -> clear();
- TimeHitMap -> clear();
- //XHitMap -> clear();
- //YHitMap -> clear();
- //PosXHitMap -> clear();
- //PosYHitMap -> clear();
- //PosZHitMap -> clear();
- //AngThetaHitMap -> clear();
- //AngPhiHitMap -> clear();
-
- /*
- CsIDetectorNumberHitMap -> clear();
- CsICrystalNumberHitMap -> clear(); // added by Marc
- CsIStageEnergyHitMap -> clear();
- */
-
-
-}
-
-
-
-void FatimaDetector::InitializeScorers()
-{
- // LaBr3 Associate Scorer
- m_LaBr3StageScorer = new G4MultiFunctionalDetector("LaBr3StageScorerFatimaDetector");
-
- // G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber", "FatimaDetector", 0);
- G4VPrimitiveScorer* DetNbr = new FATIMAScorerLaBr3StageDetectorNumber("DetectorNumber", "FatimaDetector", 0);
- // G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime","FatimaDetector", 0);
- G4VPrimitiveScorer* TOF = new FATIMAScorerLaBr3StageTOF("StripTime","FatimaDetector", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","FatimaDetector", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","FatimaDetector", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","FatimaDetector", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","FatimaDetector", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","FatimaDetector", 0);
-
- G4VPrimitiveScorer* Energy = new FATIMAScorerLaBr3StageEnergy("StripEnergy", "FatimaDetector", 0);
- G4VPrimitiveScorer* CrystNbr = new FATIMAScorerLaBr3StageCrystal("CrystalNumber", "FatimaDetector", 0);
-
- // G4VPrimitiveScorer* StripPositionX = new FATIMAcorerLaBr3StageFrontStripDummyShape("StripIDFront", 0, NumberOfStrips);
- // G4VPrimitiveScorer* StripPositionY = new FATIMAScorerLaBr3StageBackStripDummyShape("StripIDBack", 0, NumberOfStrips);
-
- //and register it to the multifunctionnal detector
- m_LaBr3StageScorer->RegisterPrimitive(DetNbr);
- m_LaBr3StageScorer->RegisterPrimitive(CrystNbr);
- m_LaBr3StageScorer->RegisterPrimitive(Energy);
- m_LaBr3StageScorer->RegisterPrimitive(TOF);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionX);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesX);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesZ);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
-
-
-
- /*
-
- // Second stage Associate Scorer
- m_CsIStageScorer = new G4MultiFunctionalDetector("CsIStageScorerFatimaDetector");
-
- G4VPrimitiveScorer* CsIDetNbr = new FATIMAScorerCsIStageDetectorNumber("CsIDetectorNumber", "FatimaDetector", 0);
- G4VPrimitiveScorer* CsICryNbr = new FATIMAScorerCsIStageCrystalNumber("CsICrystalNumber", "FatimaDetector", 0); // added by Marc
- G4VPrimitiveScorer* CsIStageEnergy = new FATIMAScorerCsIStageEnergy("CsIStageEnergy", "FatimaDetector", 0);
-
- m_CsIStageScorer->RegisterPrimitive(CsIDetNbr);
- m_CsIStageScorer->RegisterPrimitive(CsICryNbr); // Added by Marc
- m_CsIStageScorer->RegisterPrimitive(CsIStageEnergy);
-
- */
-
- // Add All Scorer to the Global Scorer Manager
- G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3StageScorer);
- //G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIStageScorer);
-
-}
diff --git a/NPSimulation/Fatima/FatimaDetector.hh b/NPSimulation/Fatima/FatimaDetector.hh
deleted file mode 100755
index d4e8d22a35923f8821bba2fa258293ee1948be89..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaDetector.hh
+++ /dev/null
@@ -1,191 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: Define a phoswich module for the Fatima detector. *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#ifndef FatimaDetector_h
-#define FatimaDetector_h 1
-
-// C++ headers
-#include <vector>
-
-// NPTool header
-#include "FatimaModule.hh"
-#include "TInteractionCoordinates.h"
-
-using namespace CLHEP;
-
-
-
-class FatimaDetector : public FatimaModule
-{
- ////////////////////////////////////////////////////
- /////// Default Constructor and Destructor /////////
- ////////////////////////////////////////////////////
-public:
- FatimaDetector();
- virtual ~FatimaDetector();
-
- ////////////////////////////////////////////////////
- //////// Specific Function of this Class ///////////
- ////////////////////////////////////////////////////
-public:
- // By Position Method
- void AddModule(G4ThreeVector TL,
- G4ThreeVector BL,
- G4ThreeVector BR,
- G4ThreeVector CT);
-
- // By Angle Method
- void AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w);
-
- // Effectively construct Volume
- // Avoid to have two time same code for Angle and Point definition
- void VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world);
-
-
- ///////////////////////////////////////////
- //// Inherite from FatimaModule class /////
- ///////////////////////////////////////////
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
-
- // 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();
-
- // Initialize all scorers necessary for the detector
- void InitializeScorers();
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
-
- // Give the static TInteractionCoordinates from VDetector to the classes
- // deriving from FatimaModule
- // This is mandatory since the Fatima*** does not derive from VDetector
- void SetInterCoordPointer(TInteractionCoordinates* interCoord);
- TInteractionCoordinates* GetInterCoordPointer() {return ms_InterCoord;};
-
-
- ////////////////////////////////////////////////////
- ///////////////Private intern Data//////////////////
- ////////////////////////////////////////////////////
-private:
- // Interaction Coordinates coming from VDetector through the
- // SetInteractionCoordinatesPointer method
- TInteractionCoordinates* ms_InterCoord;
-
-
-
- // True if Define by Position, False is Define by angle
- vector<bool> m_DefinitionType ;
-
- // Used for "By Point Definition"
- vector<G4ThreeVector> m_X1_Y1 ; // Top Left Corner Position Vector
- vector<G4ThreeVector> m_X1_Y128 ; // Bottom Left Corner Position Vector
- vector<G4ThreeVector> m_X128_Y1 ; // Bottom Right Corner Position Vector
- vector<G4ThreeVector> m_X128_Y128 ; // Center Corner Position Vector
-
- // Used for "By Angle Definition"
- vector<G4double> m_R ; // |
- vector<G4double> m_Theta ; // > Spherical coordinate of Strips Silicium Plate
- vector<G4double> m_Phi ; // |
-
- vector<G4double> m_beta_u ; // |
- vector<G4double> m_beta_v ; // > Tilt angle of the Telescope
- vector<G4double> m_beta_w ; // |
-
- G4ThreeVector momentum;
-
-};
-
-
-
-namespace FATIMADETECTOR
-{
-
-
-
- // Resolution
-// const G4double ResoFirstStage = 0; // = 50 keV of Resolution // Unit is MeV/2.35
- const G4double ResoFirstStage = 0.0099; // = 3.5% at .662MeV of Resolution // Unit is MeV/2.35
- // const G4double ResoSecondStage = 0.0; // = 13% at .662 MeV of resolution // Unit is MeV/2.35
- // const G4double ResoThirdStage = 0.0; // = 50 keV of resolution // Unit is MeV/2.35
- const G4double ResoTimeGpd = 0.212765957;// = 500ps // Unit is ns/2.35
-
- // Geometry for the mother volume containing the different layers of your dummy shape module
- const G4double FaceFront = 7.5*cm; //;
- const G4double Length = 26.33*cm; //
-
- // Geometry for the phoswitch
- // LaBr3
- const G4double LaBr3Face = 3.81*cm;
- const G4double LaBr3Thickness = 5.08*cm; // for detector 2x2x2
-
- // Al Can
- const G4double CanOuterDiameter = 4.3*cm;
- const G4double CanInnerDiameter = 4.15*cm;
- const G4double CanLength = 4.33*cm; // = (5.08-(22-0.85)+0.1 with 0.85=(26.33-0.1-5.08)-22
- // Al front Window
- const G4double WinOuterDiameter = 4.15*cm;
- const G4double WinInnerDiameter = 0*cm;
- const G4double WinLength = 0.1*cm; //
-
- // PM (TODO Change CsI to Alu)
- const G4double CsIFace = FaceFront;
- const G4double CsIThickness = 22.*cm; // for detector
- //const G4double CsIThickness = 0.001*cm; // For LaBr long
-
- // Lead tube (TODO Change CsI to Alu)
- const G4double LeadAMinR = CanOuterDiameter;
- const G4double LeadAMaxR = CanOuterDiameter+1*cm; // for detector
- const G4double LeadALength = 4.33*cm; // For LaBr long
-
- const G4double LeadBMinR = CanOuterDiameter;
- const G4double LeadBMaxR = FaceFront; // for detector
- const G4double LeadBLength = 1.*cm; // For LaBr long
-
- // Position
- const G4double LaBr3Stage_PosZ = -Length*0.5 + 0.5*LaBr3Thickness + 0.1*cm;
- const G4double LaBr3Can_PosZ = -Length*0.5 + 0.5*CanLength;
- const G4double LaBr3Win_PosZ = -Length*0.5 + 0.5*WinLength;
-
- const G4double CsIStage_PosZ = -Length*0.5 + (Length-CsIThickness) + 0.5*CsIThickness;
-
- const G4double LeadAShield_PosZ = -Length*0.5 + 0.5*LeadALength;
- const G4double LeadBShield_PosZ = -Length*0.5 + LeadALength - 0.5*LeadBLength ;
-
-}
-
-#endif
diff --git a/NPSimulation/Fatima/FatimaModule.cc b/NPSimulation/Fatima/FatimaModule.cc
deleted file mode 100755
index b53bfa4268646016f7b6819d9948dc8265d8332e..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaModule.cc
+++ /dev/null
@@ -1,60 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/14 *
- *---------------------------------------------------------------------------*
- * Decription: This class is an Abstract Base Class (ABC) from which should *
- * derive all different modules from the Fatima detector. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#include "FatimaModule.hh"
-#include "RootOutput.h"
-
-
-TFatimaData *FatimaModule::ms_Event = 0;
-
-
-
-FatimaModule::FatimaModule()
-{
- if (ms_Event == 0) ms_Event = new TFatimaData();
-
- InitializeRootOutput();
- InitializeIndex();
-}
-
-
-
-FatimaModule::~FatimaModule()
-{
-}
-
-
-
-void FatimaModule::InitializeRootOutput()
-{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- // if the branch does not exist yet, create it
- if (!pTree->GetBranch("FATIMA"))
- pTree->Branch("FATIMA", "TFatimaData", &ms_Event);
-}
-
-
-
-void FatimaModule::InitializeIndex()
-{
- m_index["LaBr"] = 0; // 24 LaBr max
-}
diff --git a/NPSimulation/Fatima/FatimaModule.hh b/NPSimulation/Fatima/FatimaModule.hh
deleted file mode 100755
index 45bd3a5b4d0da6cdac5b2406cbba32abf3b6ab70..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaModule.hh
+++ /dev/null
@@ -1,94 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class is an Abstract Base Class (ABC) from which should *
- * derive all different modules from the Fatima detector. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#ifndef FatimaModule_h
-#define FatimaModule_h 1
-
-// C++ headers
-#include <string>
-#include <vector>
-
-// G4 headers
-#include "G4LogicalVolume.hh"
-#include "G4Event.hh"
-#include "G4MultiFunctionalDetector.hh"
-
-// NPTool - ROOT headers
-#include "TInteractionCoordinates.h"
-#include "TFatimaData.h"
-
-using namespace std;
-
-
-
-class FatimaModule
-{
-public:
- FatimaModule();
- virtual ~FatimaModule();
-
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- virtual void ReadConfiguration(string Path) = 0;
-
- // Construct detector and inialise sensitive part.
- virtual void ConstructDetector(G4LogicalVolume* world) = 0;
-
- // Read sensitive part of a the FatimaModule detector and fill the Root tree.
- virtual void ReadSensitive(const G4Event* event) = 0;
-
- // Add Detector branch to the ROOT output tree
- virtual void InitializeRootOutput();
-
- // Initialize all scorers necessary for each detector
- virtual void InitializeScorers() = 0;
-
- // Give the static TInteractionCoordinates from VDetector to the classes
- // deriving from FatimaModule
- // This is mandatory since the Fatima*** does not derive from VDetector
- virtual void SetInterCoordPointer(TInteractionCoordinates* interCoord) = 0;
- virtual TInteractionCoordinates* GetInterCoordPointer() = 0;
-
- // Initialize the Index map for the different modules of Fatima
- void InitializeIndex();
-
-public:
- TFatimaData* GetEventPointer() {return ms_Event;};
-
-protected:
- // Class to store the result of simulation
- static TFatimaData* ms_Event;
-
- // Set to true if you want to see Telescope Frame in your visualisation
- bool m_non_sensitive_part_visiualisation;
-
-protected:
- // First stage Associate Scorer
- G4MultiFunctionalDetector* m_LaBr3StageScorer;
-
- // Second stage Associate Scorer
- G4MultiFunctionalDetector* m_CsIStageScorer;
-
-protected:
- map<string, int> m_index;
-};
-
-#endif
diff --git a/NPSimulation/Fatima/FatimaScorers.cc b/NPSimulation/Fatima/FatimaScorers.cc
deleted file mode 100755
index a869412d6add588e2ba835da970eaa95ac02e3ef..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaScorers.cc
+++ /dev/null
@@ -1,514 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : 02/07/13 *
- *---------------------------------------------------------------------------*
- * Decription: This class holds all the scorers needed by the *
- * Fatima*** objects. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-// G4 headers
-#include "G4UnitsTable.hh"
-
-// NPTool headers
-#include "ObsoleteGeneralScorers.hh"
-#include "FatimaScorers.hh"
-
-#include "FatimaDetector.hh"
-
-using namespace FATIMADETECTOR;
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Added by Adrien MATTA:
-// Those Scorer use TrackID as map index. This way ones can rebuild energy deposit,
-// time of flight or position,... particle by particle for each event. Because standard
-// scorer provide by G4 don't work this way but using a global ID for each event you should
-// not use those scorer with some G4 provided ones or being very carefull doing so.
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage Detector Number Scorer (deal with multiple particle hit)
-FATIMAScorerLaBr3StageDetectorNumber::FATIMAScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerLaBr3StageDetectorNumber::~FATIMAScorerLaBr3StageDetectorNumber()
-{
-}
-
-G4bool FATIMAScorerLaBr3StageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
- return TRUE;
-}
-void FATIMAScorerLaBr3StageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerLaBr3StageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerLaBr3StageDetectorNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerLaBr3StageDetectorNumber::DrawAll()
-{
-}
-
-void FATIMAScorerLaBr3StageDetectorNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
- << G4endl;
- }
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage Energy Scorer (deal with multiple particle hit)
-FATIMAScorerLaBr3StageEnergy::FATIMAScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerLaBr3StageEnergy::~FATIMAScorerLaBr3StageEnergy()
-{
-}
-
-G4bool FATIMAScorerLaBr3StageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
- POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, edep);
- return TRUE;
-}
-
-void FATIMAScorerLaBr3StageEnergy::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerLaBr3StageEnergy::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerLaBr3StageEnergy::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerLaBr3StageEnergy::DrawAll()
-{
-}
-
-void FATIMAScorerLaBr3StageEnergy::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
-FATIMAScorerLaBr3StageCrystal::FATIMAScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerLaBr3StageCrystal::~FATIMAScorerLaBr3StageCrystal()
-{
-}
-
-G4bool FATIMAScorerLaBr3StageCrystal::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
-
- //G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1);
- //Adde by Anna:
- G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); // daughter crystal volume
- //depth:0 is the world, 1 is the mother...
- //AC instead of GetReplicaNumber(1)
- //cout<<"****replica****"<<endl;
- //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
- //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- //EvtMap->add(DetNbr + index, CrystalNbr); // Marc
- //cout<<"index "<<index<<endl;
- //cout<<"DetNbr + index "<<DetNbr + index<<endl;
- EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);// modified by Anna
- return TRUE;
-
-}
-
-void FATIMAScorerLaBr3StageCrystal::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerLaBr3StageCrystal::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerLaBr3StageCrystal::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerLaBr3StageCrystal::DrawAll()
-{
-}
-
-void FATIMAScorerLaBr3StageCrystal::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
- << G4endl;
- }
-}
-
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage ToF Scorer (deal with multiple particle hit)
-FATIMAScorerLaBr3StageTOF::FATIMAScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerLaBr3StageTOF::~FATIMAScorerLaBr3StageTOF()
-{
-}
-
-G4bool FATIMAScorerLaBr3StageTOF::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get TOF
- G4double TOF = aStep->GetPreStepPoint()->GetGlobalTime();
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, TOF);
- return TRUE;
-}
-
-void FATIMAScorerLaBr3StageTOF::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerLaBr3StageTOF::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerLaBr3StageTOF::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerLaBr3StageTOF::DrawAll()
-{
-}
-
-void FATIMAScorerLaBr3StageTOF::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
-}
-
-
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// SecondStage (CsI) Detector Number Scorer (deal with multiple particle hit)
-FATIMAScorerCsIStageDetectorNumber::FATIMAScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerCsIStageDetectorNumber::~FATIMAScorerCsIStageDetectorNumber()
-{
-}
-
-G4bool FATIMAScorerCsIStageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
- return TRUE;
-}
-void FATIMAScorerCsIStageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerCsIStageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerCsIStageDetectorNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerCsIStageDetectorNumber::DrawAll()
-{
-}
-
-void FATIMAScorerCsIStageDetectorNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// CsIStage Energy Scorer (deal with multiple particle hit)
-FATIMAScorerCsIStageEnergy::FATIMAScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerCsIStageEnergy::~FATIMAScorerCsIStageEnergy()
-{
-}
-
-G4bool FATIMAScorerCsIStageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
- POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- //if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, edep);
- return TRUE;
-}
-
-void FATIMAScorerCsIStageEnergy::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerCsIStageEnergy::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerCsIStageEnergy::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerCsIStageEnergy::DrawAll()
-{
-}
-
-void FATIMAScorerCsIStageEnergy::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// all added by Anna:
-
-// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
-FATIMAScorerCsIStageCrystalNumber::FATIMAScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-FATIMAScorerCsIStageCrystalNumber::~FATIMAScorerCsIStageCrystalNumber()
-{
-}
-
-G4bool FATIMAScorerCsIStageCrystalNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
-
- // get energy
- G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
- //depth:0 is the world, 1 is the mother...
- //AC instead of GetReplicaNumber(1)
- //cout<<"****replica****"<<endl;
- //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
- //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- //cout<<"index "<<index<<endl;
- //cout<<"DetNbr + index "<<DetNbr + index<<endl;
- EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);
- return TRUE;
-
-}
-
-void FATIMAScorerCsIStageCrystalNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void FATIMAScorerCsIStageCrystalNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void FATIMAScorerCsIStageCrystalNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void FATIMAScorerCsIStageCrystalNumber::DrawAll()
-{
-}
-
-void FATIMAScorerCsIStageCrystalNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
- << G4endl;
- }
-}
diff --git a/NPSimulation/Fatima/FatimaScorers.hh b/NPSimulation/Fatima/FatimaScorers.hh
deleted file mode 100755
index 2db2cb6cae09c5844c9159146e83c570fe743525..0000000000000000000000000000000000000000
--- a/NPSimulation/Fatima/FatimaScorers.hh
+++ /dev/null
@@ -1,198 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 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: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/01/13 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class holds all the scorers needed by the *
- * GaspardTracker*** objects. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#ifndef FATIMAScorer_h
-#define FATIMAScorer_h 1
-
-#include "G4VPrimitiveScorer.hh"
-#include "G4THitsMap.hh"
-
-namespace FATIMASCORERS
-{
-
-
-class FATIMAScorerLaBr3StageDetectorNumber : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerLaBr3StageDetectorNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-
-class FATIMAScorerLaBr3StageEnergy : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerLaBr3StageEnergy();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-
-
-class FATIMAScorerLaBr3StageCrystal : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerLaBr3StageCrystal();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-class FATIMAScorerLaBr3StageTOF : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerLaBr3StageTOF();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-
-
-class FATIMAScorerCsIStageDetectorNumber : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerCsIStageDetectorNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-// Added by Anna
-class FATIMAScorerCsIStageCrystalNumber : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerCsIStageCrystalNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-class FATIMAScorerCsIStageEnergy : public G4VPrimitiveScorer
-{
-public: // with description
- FATIMAScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~FATIMAScorerCsIStageEnergy();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-}
-
-using namespace FATIMASCORERS;
-#endif
diff --git a/NPSimulation/Sharc/Sharc.hh b/NPSimulation/Sharc/Sharc.hh
index 944f2eff534d365a5e4e68b6df1a422781b2b788..8371d3e8e529e1795747fab247cc2fca24710d56 100644
--- a/NPSimulation/Sharc/Sharc.hh
+++ b/NPSimulation/Sharc/Sharc.hh
@@ -40,7 +40,7 @@
// NPLib
#include "TSharcData.h"
using namespace std;
-
using namespace CLHEP;
+using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace SHARC{