diff --git a/Inputs/DetectorConfiguration/hira.detector b/Inputs/DetectorConfiguration/hira.detector
index 5ab9fcd6e020ca3ac6e30e80f8919f54164b03c3..d331d4440a181956a53dd040614df43df0608376 100755
--- a/Inputs/DetectorConfiguration/hira.detector
+++ b/Inputs/DetectorConfiguration/hira.detector
@@ -7,7 +7,7 @@ GeneralTarget
% Radius in mm
% Temperature in K, Pressure in bar
Target
- THICKNESS= 25
+ THICKNESS= 37
ANGLE= 0
RADIUS= 10
MATERIAL= CH2
diff --git a/Inputs/DetectorConfiguration/lassaFullArray.detector b/Inputs/DetectorConfiguration/lassaFullArray.detector
new file mode 100644
index 0000000000000000000000000000000000000000..3b775e59640deff180f39819d9feafa7c7a4503f
--- /dev/null
+++ b/Inputs/DetectorConfiguration/lassaFullArray.detector
@@ -0,0 +1,84 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GeneralTarget
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Thickness in micrometer
+% Radius in mm
+% Temperature in K, Pressure in bar
+Target
+ THICKNESS= 0
+ ANGLE= 0
+ RADIUS= 10
+ MATERIAL= CD2
+ X= 0
+ Y= 0
+ Z= 0
+
+%%%%%%%%%%%%%%%%%%%%%
+LASSAArray
+%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%% Telescope 0 %%%%%%%
+LassaTelescope
+A= 7.14047 29.6923 196.049
+B= 4.37927 73.977 182.139
+C= 53.0511 30.2433 188.69
+D= 50.2899 74.528 174.78
+
+
+%%%%%%% Telescope 1 %%%%%%%
+LassaTelescope
+A= 69.3517 30.7385 183.164
+B= 62.4168 75.0505 170.892
+C= 110.505 31.1849 161.52
+D= 103.57 75.4969 149.248
+
+
+%%%%%%% Telescope 2 %%%%%%%
+LassaTelescope
+A= 123.949 31.5221 151.611
+B= 113.682 75.8459 142.008
+C= 155.962 31.6308 117.885
+D= 145.694 75.9546 108.282
+
+
+%%%%%%% Telescope 3 %%%%%%%
+LassaTelescope
+A= 39.8538 -28.4586 194.329
+B= 39.2873 18.0371 194.054
+C= 83.9296 -28.0093 179.517
+D= 83.3631 18.4864 179.242
+
+
+%%%%%%% Telescope 4 %%%%%%%
+LassaTelescope
+A= 99.6612 -27.5979 171.242
+B= 99.2749 18.8996 170.947
+C= 136.8 -27.4665 143.261
+D= 136.414 19.0309 142.966
+
+
+%%%%%%% Telescope 5 %%%%%%%
+LassaTelescope
+A= 5.43374 -85.5137 183.786
+B= 7.42935 -40.9223 196.819
+C= 51.2967 -85.3249 176.117
+D= 53.2923 -40.7334 189.151
+
+
+%%%%%%% Telescope 6 %%%%%%%
+LassaTelescope
+A= 63.6721 -85.4182 171.982
+B= 69.4239 -40.7797 183.668
+C= 104.796 -85.0359 150.281
+D= 110.548 -40.3974 161.967
+
+
+%%%%%%% Telescope 7 %%%%%%%
+LassaTelescope
+A= 115.146 -84.6098 142.665
+B= 124.503 -39.9765 151.751
+C= 147.179 -84.4639 108.959
+D= 156.536 -39.8307 118.045
+
+
diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index 66b22a966b951845f8a627fb5d72aa51efb728fd..2c4bf8b610a8ea35d9f21ecca33a4d4f7e939488 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 0
- EnergyHigh= 26
+ EnergyLow= 120
+ EnergyHigh= 200
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 45
+ HalfOpenAngleMax= 12
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 6f387b1eaaaa59f1b1c3a6614aa61827d5d222a6..02b1c4dfeae0658076e3c8cff6632c4995e62693 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -7,7 +7,7 @@ Isotropic
EnergyLow= 5
EnergyHigh= 5
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 180
+ HalfOpenAngleMax= 5
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index cf3a617525e74a14a0cd143ba0045201d79f8599..fe74556d683666950e9bdeb3785c9f507634b5af 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 17
+ EnergyHigh= 150
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 45
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/NPLib/Detectors/Lassa/TLassaData.h b/NPLib/Detectors/Lassa/TLassaData.h
index a5165b478b15f4f55fd4e3f0ace3da46932bc656..85ecb3b62a6625584dccd5e58df82df825c1c97a 100644
--- a/NPLib/Detectors/Lassa/TLassaData.h
+++ b/NPLib/Detectors/Lassa/TLassaData.h
@@ -162,7 +162,7 @@ class TLassaData : public TObject {
UShort_t GetLassaCsITCristalNbr(const Int_t i) const {return fLassa_CsIT_CristalNbr[i];}
Double_t GetLassaCsITTime(const Int_t i) const {return fLassa_CsIT_Time[i];}
- ClassDef(TLassaData,2) // HiraData structure
+ ClassDef(TLassaData,2) // LassaData structure
};
#endif
diff --git a/NPLib/Detectors/Lassa/TLassaPhysics.cxx b/NPLib/Detectors/Lassa/TLassaPhysics.cxx
index 861c928e53ed40a9b580532a7a896e576388f124..55d053c425f4b78fc0c503bb3f7d1ff345481f14 100644
--- a/NPLib/Detectors/Lassa/TLassaPhysics.cxx
+++ b/NPLib/Detectors/Lassa/TLassaPhysics.cxx
@@ -277,7 +277,8 @@ void TLassaPhysics::ReadConfiguration(string Path){
}
}
- ReadAnalysisConfig();
+ InitializeStandardParameter();
+ ReadAnalysisConfig();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -398,10 +399,6 @@ void TLassaPhysics::BuildSimplePhysicalEvent(){
m_ThickSi_EYMult = m_PreTreatedData->GetLassaStripYEMult();
m_CsI_EMult = m_PreTreatedData->GetLassaCsIEMult();
-// cout << "Si_X mult: " << m_ThickSi_EXMult << endl;
-// cout << "Si_Y mult: " << m_ThickSi_EYMult << endl;
-// cout << "CsI mult: " << m_CsI_EMult << endl;
-
/*for(unsigned int i = 0 ; i < m_ThinSi_EMult ; ++i){
TelescopeNumber.push_back(m_PreTreatedData->GetLassaThinSiStripEDetectorNbr(i));
@@ -499,22 +496,21 @@ void TLassaPhysics::PreTreat(){
////////////////////////////////////////////////////////////////////////////
bool TLassaPhysics :: IsValidChannel(const string DetectorType, const int telescope , const int channel){
- if(DetectorType == "CsI")
- *(m_CsIChannelStatus[telescope].begin()+channel); //return *(m_CsIChannelStatus[telescope-1].begin()+channel-1);
+ if(DetectorType == "CsI")
+ return *(m_CsIChannelStatus[telescope].begin()+channel);
- else if(DetectorType == "X")
- return *(m_XChannelStatus[telescope].begin()+channel); //return *(m_XChannelStatus[telescope-1].begin()+channel-1);
+ else if(DetectorType == "X")
+ return *(m_XChannelStatus[telescope].begin()+channel);
- else if(DetectorType == "Y")
- return *(m_YChannelStatus[telescope].begin()+channel); //return *(m_YChannelStatus[telescope-1].begin()+channel-1);
+ else if(DetectorType == "Y")
+ return *(m_YChannelStatus[telescope].begin()+channel);
- return false;
+ else return false;
}
////////////////////////////////////////////////////////////////////////////
int TLassaPhysics :: CheckEvent(const double XEMult, const double YEMult){
- //Check that the multiplicity in the X-strip is the same as the y-strip
if( XEMult == YEMult ){
return 1; //Regular Event
@@ -526,122 +522,107 @@ int TLassaPhysics :: CheckEvent(const double XEMult, const double YEMult){
///////////////////////////////////////////////////////////////////////////
void TLassaPhysics::ReadAnalysisConfig(){
- bool ReadingStatus = false;
-
- for(int i = 0; i < m_NumberOfTelescope; i++){ //for(int i = 0; i < 8; i++){
- for(int j = 0; j < 16; j++){
- m_XChannelStatus[i][j] = true;
- m_YChannelStatus[i][j] = true;
- }
- for(int k = 0; k < 4; k++){
- m_CsIChannelStatus[i][k] = true;
- }
- }
+ bool ReadingStatus = false;
- // path to file
- string FileName = "./configs/ConfigLassaTest.dat"; //"./configs/ConfigLassa.dat";
+ // path to file
+ string FileName = "./configs/ConfigLassa.dat";
- // open analysis config file
- ifstream AnalysisConfigFile;
- AnalysisConfigFile.open(FileName.c_str());
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
- if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigLassa.dat found: Default parameters loaded for Analysis " << FileName << endl;
- return;
- }
- cout << " Loading user parameters for Analysis from ConfigLassa.dat " << endl;
-
- // Save it in a TAsciiFile
- TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
- asciiConfig->AppendLine("%%% ConfigLassa.dat %%%");
- asciiConfig->Append(FileName.c_str());
- asciiConfig->AppendLine("");
- // read analysis config file
- string LineBuffer,DataBuffer,whatToDo;
- while (!AnalysisConfigFile.eof()) {
- // Pick-up next line
- getline(AnalysisConfigFile, LineBuffer);
-
- // search for "header"
- if (LineBuffer.compare(0, 11, "ConfigLassa") == 0) ReadingStatus = true;
-
- // loop on tokens and data
- while (ReadingStatus ) {
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigLassa.dat found: Default parameters loaded for Analysis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameters for Analysis from ConfigLassa.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigLassa.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer,DataBuffer,whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ if (LineBuffer.compare(0, 11, "ConfigLassa") == 0) ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus ) {
- whatToDo="";
- AnalysisConfigFile >> whatToDo;
-
- // Search for comment symbol (%)
- if (whatToDo.compare(0, 1, "%") == 0) {
- AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
- }
+ whatToDo="";
+ AnalysisConfigFile >> whatToDo;
- else if (whatToDo=="MAX_STRIP_MULTIPLICITY") {
- AnalysisConfigFile >> DataBuffer;
- m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() );
- cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
- }
-
- else if (whatToDo=="STRIP_ENERGY_MATCHING_SIGMA") {
- AnalysisConfigFile >> DataBuffer;
- m_StripEnergyMatchingSigma = atof(DataBuffer.c_str() );
- cout << "STRIP ENERGY MATCHING SIGMA " << m_StripEnergyMatchingSigma <<endl;
- }
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
- else if (whatToDo=="STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA") {
- AnalysisConfigFile >> DataBuffer;
- m_StripEnergyMatchingNumberOfSigma = atoi(DataBuffer.c_str() );
- cout << "STRIP ENERGY MATCHING NUMBER OF SIGMA " << m_StripEnergyMatchingNumberOfSigma << endl;
- }
+ else if (whatToDo=="MAX_STRIP_MULTIPLICITY") {
+ AnalysisConfigFile >> DataBuffer;
+ m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() );
+ cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
+ }
- else if (whatToDo== "DISABLE_ALL") {
- AnalysisConfigFile >> DataBuffer;
- cout << whatToDo << " " << DataBuffer << endl;
- int telescope = atoi(DataBuffer.substr(2,1).c_str());
- vector< bool > ChannelStatus;
- ChannelStatus.resize(128,false);
- m_XChannelStatus[telescope] = ChannelStatus; //m_XChannelStatus[telescope-1] = ChannelStatus;
- m_YChannelStatus[telescope] = ChannelStatus; //m_YChannelStatus[telescope-1] = ChannelStatus;
- ChannelStatus.resize(16,false);
- m_CsIChannelStatus[telescope] = ChannelStatus; //m_CsIChannelStatus[telescope-1] = ChannelStatus;
- }
+ else if (whatToDo=="STRIP_ENERGY_MATCHING_SIGMA") {
+ AnalysisConfigFile >> DataBuffer;
+ m_StripEnergyMatchingSigma = atof(DataBuffer.c_str() );
+ cout << "STRIP ENERGY MATCHING SIGMA " << m_StripEnergyMatchingSigma <<endl;
+ }
- else if (whatToDo == "DISABLE_CHANNEL") {
- AnalysisConfigFile >> DataBuffer;
- cout << whatToDo << " " << DataBuffer << endl;
- int telescope = atoi(DataBuffer.substr(2,1).c_str());
- int channel = -1;
- if (DataBuffer.compare(3,4,"STRX") == 0) {
- channel = atoi(DataBuffer.substr(7).c_str());
- *(m_XChannelStatus[telescope].begin()+channel) = false; //*(m_XChannelStatus[telescope-1].begin()+channel-1) = false;
- }
+ else if (whatToDo=="STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA") {
+ AnalysisConfigFile >> DataBuffer;
+ m_StripEnergyMatchingNumberOfSigma = atoi(DataBuffer.c_str() );
+ cout << "STRIP ENERGY MATCHING NUMBER OF SIGMA " << m_StripEnergyMatchingNumberOfSigma << endl;
+ }
- else if (DataBuffer.compare(3,4,"STRY") == 0) {
- channel = atoi(DataBuffer.substr(7).c_str());
- *(m_YChannelStatus[telescope].begin()+channel) = false; //*(m_YChannelStatus[telescope-1].begin()+channel-1) = false;
- }
+ else if (whatToDo== "DISABLE_ALL") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ int telescope = atoi(DataBuffer.substr(2,1).c_str());
+ vector< bool > ChannelStatus;
+ ChannelStatus.resize(16,false);
+ m_XChannelStatus[telescope] = ChannelStatus;
+ m_YChannelStatus[telescope] = ChannelStatus;
+ ChannelStatus.resize(4,false);
+ m_CsIChannelStatus[telescope] = ChannelStatus;
+ }
- else if (DataBuffer.compare(3,3,"CSI") == 0) {
- channel = atoi(DataBuffer.substr(6).c_str());
- *(m_CsIChannelStatus[telescope].begin()+channel) = false; //*(m_CsIChannelStatus[telescope-1].begin()+channel-1) = false;
- }
+ else if (whatToDo == "DISABLE_CHANNEL") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ int telescope = atoi(DataBuffer.substr(2,1).c_str());
+ int channel = -1;
+ if (DataBuffer.compare(3,4,"STRX") == 0) {
+ channel = atoi(DataBuffer.substr(7).c_str());
+ *(m_XChannelStatus[telescope].begin()+channel) = false;
+ }
- else cout << "Warning: detector type for Lassa unknown!" << endl;
+ else if (DataBuffer.compare(3,4,"STRY") == 0) {
+ channel = atoi(DataBuffer.substr(7).c_str());
+ *(m_YChannelStatus[telescope].begin()+channel) = false;
+ }
- }
+ else if (DataBuffer.compare(3,3,"CSI") == 0) {
+ channel = atoi(DataBuffer.substr(6).c_str());
+ *(m_CsIChannelStatus[telescope].begin()+channel) = false;
+ }
- //I want to print out the values of my map to see if they are filled correctly
+ else cout << "Warning: detector type for Lassa unknown!" << endl;
- else {
- ReadingStatus = false;
- }
- }
- }
+ }
+
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
-
- //I want to print out the values of my map to see if they are filled correctly
-
- for(int i = 0; i < m_NumberOfTelescope; i++){
+ /*for(int i = 0; i < m_NumberOfTelescope; i++){
for(int j = 0; j < 16; j++){
// cout << "Hello!" << endl;
cout << "Telescope: " << i << " " << "XStrip: " << j << " " << m_XChannelStatus[i][j] << " " << "YStrip: " << j << " " << m_YChannelStatus[i][j] << endl;
@@ -649,11 +630,33 @@ void TLassaPhysics::ReadAnalysisConfig(){
for(int k = 0; k < 4; k++){
cout << "Telescope: " << i << " " << "Crystal: " << k << " " << m_CsIChannelStatus[i][k] << endl;
}
- }
+ }*/
}
+///////////////////////////////////////////////////////////////////////////
+void TLassaPhysics::InitializeStandardParameter(){
+ // Enable all channel
+ vector< bool > ChannelStatus;
+ m_XChannelStatus.clear() ;
+ m_YChannelStatus.clear() ;
+ m_CsIChannelStatus.clear() ;
+
+ ChannelStatus.resize(16,true);
+ for(int i = 0 ; i < m_NumberOfTelescope ; ++i){
+ m_XChannelStatus[i] = ChannelStatus;
+ m_YChannelStatus[i] = ChannelStatus;
+ }
+
+ ChannelStatus.resize(4,true);
+ for(int i = 0 ; i < m_NumberOfTelescope ; ++i){
+ m_CsIChannelStatus[i] = ChannelStatus;
+ }
+
+ return;
+}
+
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Lassa/TLassaPhysics.h b/NPLib/Detectors/Lassa/TLassaPhysics.h
index f69a7abfce453482d21f5137ac9c2914f7cb9b87..2873a5e58e61bddcecbdadf4ca7f5b71eb75a669 100644
--- a/NPLib/Detectors/Lassa/TLassaPhysics.h
+++ b/NPLib/Detectors/Lassa/TLassaPhysics.h
@@ -116,6 +116,8 @@ class TLassaPhysics : public TObject, public NPL::VDetector
// Read the user configuration file; if no file found, load standard one
void ReadAnalysisConfig();
+
+ void InitializeStandardParameter();
// Use to access the strip position
// double GetStripPositionX( const int N , const int X , const int Y ) const{ return m_StripPositionX[N-1][X-1][Y-1] ; } ;
diff --git a/NPLib/Physics/nubtab03.asc b/NPLib/Physics/nubtab03.asc
index eda8276023eb0e957b4a2c2e9d163b31afc46178..edc1ab002d26d544c762b48b953df004ce42b2ec 100644
--- a/NPLib/Physics/nubtab03.asc
+++ b/NPLib/Physics/nubtab03.asc
@@ -10,7 +10,7 @@
005 0010 5H 32890 100 >910 ys (1/2+) 02 03Go11t 2n=100
005 0020 5He 11390 50 700 ys 30 3/2- 02 n=100
005 0030 5Li 11680 50 370 ys 30 3/2- 02 p=100
-005 0040 5Be 38000# 4000# 1/2+# 02 p ?
+005 0040 5Be 37140# 4000# 1/2+# 02 p ?
006 0010 6H 41860 260 290 ys 70 2-# 02 n ?;3n ?
006 0020 6He 17595.1 0.8 806.7 ms 1.5 0+ 02 90Ri01d B-=100;B-d=0.00028 5
006 0030 6Li 14086.793 0.015 stbl 1+ 02 IS=7.59 4
diff --git a/NPLib/scripts/Style_nptool.C b/NPLib/scripts/Style_nptool.C
index 9cd992297c2ab326d20a7c72eae61ffb40e1bd61..d77f62f5c19a08643fc1fbd12bca45a65d5a2f03 100644
--- a/NPLib/scripts/Style_nptool.C
+++ b/NPLib/scripts/Style_nptool.C
@@ -53,7 +53,7 @@ void Style_nptool(){
// y axis
//style_nptool->GetYaxis()->CenterTitle();
style_nptool->SetTitleYSize(0.07);
- style_nptool->SetTitleYOffset(1.1);
+ style_nptool->SetTitleYOffset(1.02);
style_nptool->SetLabelOffset(0.1,"Y");
style_nptool->SetLabelSize(0.04,"Y");
style_nptool->SetLabelOffset(0.006,"Y");
diff --git a/NPSimulation/Detectors/Hira/Hira.hh b/NPSimulation/Detectors/Hira/Hira.hh
index 95bd19375477ec1daa3bbe0ac348ddd830323bad..56edbb193f420a2e31eeed05eed376615b6a72c6 100644
--- a/NPSimulation/Detectors/Hira/Hira.hh
+++ b/NPSimulation/Detectors/Hira/Hira.hh
@@ -69,9 +69,9 @@ namespace HIRA
const G4double CsIThickness = 10.*cm;// + 2*MylarCsIThickness ;
const G4double CsIXFront = 34.93*mm;
- const G4double CsIXBack = 44.62*mm;
+ const G4double CsIXBack = 44.62*mm;//39.*mm
const G4double CsIYFront = 34.93*mm;
- const G4double CsIYBack = 44.62*mm;
+ const G4double CsIYBack = 44.62*mm;//39.*mm
const G4double DistInterCsI = 0.6*mm;
const G4double ClusterFaceFront = 2*CsIXFront+DistInterCsI;
diff --git a/Projects/Hira2/Analysis.cxx b/Projects/Hira2/Analysis.cxx
index b2a54d0a2fe9463bb1707ecbe2137b1c5ea8b51a..60d0f38879e42fea579f06eed13510ef4407db9b 100644
--- a/Projects/Hira2/Analysis.cxx
+++ b/Projects/Hira2/Analysis.cxx
@@ -41,7 +41,7 @@ void Analysis::Init(){
InitOutputBranch();
InitInputBranch();
Rand = TRandom3();
- TransferReaction= new NPL::Reaction("8B(p,t)6B@400");
+ TransferReaction= new NPL::Reaction("8B(p,4He)5Be@400");
//TransferReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
E_ThinSi = 0;
@@ -70,6 +70,7 @@ void Analysis::Init(){
Triton_CH2 = EnergyLoss("triton_CD2.G4table","G4Table",100 );
Deuteron_CH2 = EnergyLoss("deuteron_CH2.G4table","G4Table",100 );
He3_CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
+ He4_CH2 = EnergyLoss("alpha_CH2.G4table","G4Table",100 );
f_proton = new TF1("f_proton","1 - TMath::Exp(-(x-0.0538203)/28.3125)",0,10);
@@ -149,8 +150,9 @@ void Analysis::TreatEvent(){
else ELab = E_ThinSi + E_ThickSi;
if(ELab>0){
- //ELab = Triton_CD2.EvaluateInitialEnergy(ELab,(double)TargetThickness/2*micrometer,0);
- ELab = Triton_CH2.EvaluateInitialEnergy(ELab,(double)TargetThickness/2*micrometer,0);
+ //ELab = Deuteron_CH2.EvaluateInitialEnergy(ELab,(double)TargetThickness/2*micrometer,0);
+ //ELab = Triton_CH2.EvaluateInitialEnergy(ELab,(double)TargetThickness/2*micrometer,0);
+ ELab = He4_CH2.EvaluateInitialEnergy(ELab,(double)TargetThickness/2*micrometer,0);
//cout << TargetThickness << endl;
// ********************** Angle in the CM frame *****************************
TransferReaction -> SetNuclei3(ELab, ThetaLab*deg);
diff --git a/Projects/Hira2/Analysis.h b/Projects/Hira2/Analysis.h
index 4825c5f167208933f2f4673ba596afb12957a4ca..3263899d6cd1e9737e3fc5031d931141b5828859 100644
--- a/Projects/Hira2/Analysis.h
+++ b/Projects/Hira2/Analysis.h
@@ -78,6 +78,7 @@ private:
NPL::EnergyLoss Triton_CH2 ;
NPL::EnergyLoss Deuteron_CH2 ;
NPL::EnergyLoss He3_CD2;
+ NPL::EnergyLoss He4_CH2;
THiraPhysics* Hira;
TInitialConditions* InitialConditions;
diff --git a/Projects/Hira2/RunToTreat.txt b/Projects/Hira2/RunToTreat.txt
index de6e4f63bdd4bf060e8ad4f9803bbe05e0e1cd77..f91784f3c00f836b5e922f7e5282c9b18daae412 100755
--- a/Projects/Hira2/RunToTreat.txt
+++ b/Projects/Hira2/RunToTreat.txt
@@ -1,6 +1,6 @@
TTreeName
SimulatedTree
RootFileName
- ../../Outputs/Simulation/test.root
+ ../../Outputs/Simulation/dumb.root
%../../Outputs/Simulation/12Be_pt_iso.root
%../../Outputs/Simulation/12Be_pt_1st.root
diff --git a/Projects/Lassa/Analysis.cxx b/Projects/Lassa/Analysis.cxx
index eda5b6648d904e2f205018cea48c19fb31d65458..d8bb65d16c2bc3b7bd113b1caab3c4da9c3532a7 100644
--- a/Projects/Lassa/Analysis.cxx
+++ b/Projects/Lassa/Analysis.cxx
@@ -6,8 +6,7 @@
*****************************************************************************/
/*****************************************************************************
- ConfigLassa
- MAX_STRIP_MULTIPLICITY 1
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
* *
* Creation Date : march 2025 *
* Last update : *
@@ -20,25 +19,14 @@
* *
* *
*****************************************************************************/
-
-#include "Analysis.h"
-#include "NPAnalysisFactory.h"
-#include "NPDetectorManager.h"
-#include "NPOptionManager.h"
-#include "RootOutput.h"
-#include "RootInput.h"
-#include "TROOT.h"
-#include "TSystem.h"
-#include "TFile.h"
-#include "TString.h"
-#include "TTree.h"
-#include "TBranch.h"
-#include "TMath.h"
-#include "TInteractionCoordinates.h"
-#include <algorithm>
-
+#include<iostream>
using namespace std;
-
+#include"Analysis.h"
+#include"NPAnalysisFactory.h"
+#include"NPDetectorManager.h"
+#include"NPOptionManager.h"
+#include"RootOutput.h"
+#include"RootInput.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -48,292 +36,79 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
- Lassa= (TLassaPhysics*) m_DetectorManager->GetDetector("LASSAArray");
-
- Initial = new TInitialConditions();
- InitOutputBranch();
- InitInputBranch();
-
- totalEvents = 0;
- detectedEvents = 0;
- peakEvents = 0;
-
- maxEnergy = -1.;
- minEnergy = 10000.;
-
-
- for(int i = 0; i < 16; i++){
- for(int j = 0; j < 16; j++){
- detectedEvs[i][j] = 0;
- peakEvs[i][j] = 0;
- peakEffs[i][j] = -1.;
-
- if(i < 8 && j < 8){
- peakEffsAvg[i][j] = -1.;
- }
- }
- }
-
- for(int i = 0; i < 5; i++){
-
- fitFunction = "";
- Beta.str("");
- Mu.str("");
-
- fitFunction = "1-exp(-1*";
- Beta << par[i][0];
- fitFunction += Beta.str();
-
- fitFunction += "*(x-";
- Mu << par[i][1];
- fitFunction += Mu.str();
- fitFunction += "))";
-
- ReactionLoss[i] = new TF1("Reaction Loss Fit Function",fitFunction.c_str(),0,100);
-
- }
-
-
-
- for(int i = 0; i < 5; i++){
- for(int j = 0; j < 2; j++){
-
- EnergyRange[i][j].clear();
-
- }
- }
-
-
- // Declare the histograms for the geometric
- // and peak efficiency
-
- // Peak Efficiencies
-
- h_peakEff = new TH2F("Peak Efficiencies", "Peak Efficiencies",16,0,16,16,0,16);
- h_peakEff->GetXaxis()->SetTitle("X-Strip");
- h_peakEff->GetYaxis()->SetTitle("Y-Strip");
-
- h_peakEffAvg = new TH2F("Average Peak Efficiencies","Average Peak Efficiencies",8,0,8,8,0,8);
- h_peakEffAvg->GetXaxis()->SetTitle("X-Strip");
- h_peakEffAvg->GetYaxis()->SetTitle("Y-Strip");
-
- // Geometric Efficiencies
-
- hDetecTheta = new TH1F("hDetecTheta","DetecTheta",180,0,180);
- hDetecTheta->GetXaxis()->SetTitle("#Theta (deg)");
- hDetecTheta->GetYaxis()->SetTitle("Counts");
-
- hEmittTheta = new TH1F("hEmittTheta","EmittTheta",180,0,180);
- hEmittTheta->GetXaxis()->SetTitle("#Theta (deg)");
- hEmittTheta->GetYaxis()->SetTitle("Counts");
-
- hDetecThetaVsPhi = new TH2F("hDetecThetaVsPhi", "hDetecThetaVsPhi", 180,0,180,360,-180,180);
- hDetecThetaVsPhi->GetXaxis()->SetTitle("#Theta (deg)");
- hDetecThetaVsPhi->GetYaxis()->SetTitle("#Phi (deg)");
-
- hGeometricEfficiency = new TH1F("hEfficiency","Efficiency",180,0,180);
- hGeometricEfficiency->GetXaxis()->SetTitle("#Theta (deg)");
- hGeometricEfficiency->GetYaxis()->SetTitle("#epsilon (%)");
-
- // Initialize the Energy loss table for each species use.
-
- for(int i = 0; i < 5; i++){
-
- EnergyLoss_Table = "_CsI.G4table";
- EnergyLoss_Table = speciesNameArray[i] + EnergyLoss_Table;
-
- EnergyLoss_CsI[i] = EnergyLoss(EnergyLoss_Table.c_str(),"G4Table",100);
- }
-
- peakEffSum_meat = 0.;
- peakEffSquaredSum_meat = 0.;
- peakEffAvg_meat = 0.;
- peakEffSquaredAvg_meat = 0.;
- N_meat = 0.;
- sigma_meat = 0.;
-
- peakEffSum_edge = 0.;
- peakEffSquaredSum_edge = 0.;
- peakEffAvg_edge = 0.;
- peakEffSquaredAvg_edge = 0.;
- N_edge = 0.;
- sigma_edge = 0.;
+ Lassa = (TLassaPhysics*) m_DetectorManager->GetDetector("LASSAArray");
+ InitialConditions = new TInitialConditions();
+ InitOutputBranch();
+ InitInputBranch();
+ EDelta = 0.;
+ totalEvents = 0;
+ detectedEvents = 0;
+ peakEvents = 0;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
- // Reinitiate calculated variable
- ReInitValue();
-
- totalEvents++;
-
- totalEnergy = Initial->GetKineticEnergy(0);
- speciesName = Initial->GetParticleName(0);
-
- for(int i = 0; i < 5; i++){
- if(speciesName == speciesNameArray[i]){
- species = i;
- }
- }
-
- if(maxEnergy < totalEnergy){
- maxEnergy = totalEnergy;
- }
-
- if(minEnergy > totalEnergy){
- minEnergy = totalEnergy;
- }
-
- if(minEnergy < .0001){
- minEnergy = 0.;
- }
+ // Reinitiate calculated variable
+ ReInitValue();
+
+ totalEvents++;
- EDelta = .05*totalEnergy;
+ double BeamEnergy = 120.0;
+ //double BeamEnergy = InitialConditions->GetIncidentInitialKineticEnergy();
- ThetaLab_Init = Initial->GetThetaLab_WorldFrame(0);
+ //cout << BeamEnergy << endl;
+
+ EDelta = 2.0;
- hEmittTheta->Fill(ThetaLab_Init);
+ //cout << Lassa->ThickSi_E.size() << endl;
+ ///////////////////////////LOOP on Lassa Hit//////////////////////////////////
+ if(Lassa->ThickSi_E.size() == 1){
- ///////////////////////////LOOP on Lassa Hit//////////////////////////////////
-// if(Lassa -> ThickSi_E.size() == 1){
- if(Lassa -> ThickSi_E.size() >= 1){
+ detectedEvents++;
- TelescopeNumber = Lassa->TelescopeNumber[0];
- XStrip = Lassa->ThickSi_StripNumberX[0];
- YStrip = Lassa->ThickSi_StripNumberY[0];
+ //for(unsigned int countLassa = 0 ; countLassa < Lassa->ThickSi_E.size(); countLassa++){
+ TelescopeNumber = Lassa->TelescopeNumber[0];
- X = Lassa->GetPositionOfInteraction(0).X();
- Y = Lassa->GetPositionOfInteraction(0).Y();
- Z = Lassa->GetPositionOfInteraction(0).Z();
+ X = Lassa->GetPositionOfInteraction(0).X();
+ Y = Lassa->GetPositionOfInteraction(0).Y();
+ Z = Lassa->GetPositionOfInteraction(0).Z();
- TVector3 PositionOnLassa = TVector3(X,Y,Z);
- TVector3 ZUnit = TVector3(0,0,1);
+ TVector3 PositionOnLassa = TVector3(X,Y,Z);
+ TVector3 ZUnit = TVector3(0,0,1);
- double X_target = Initial->GetIncidentPositionX();
- double Y_target = Initial->GetIncidentPositionY();
- double Z_target = Initial->GetIncidentPositionZ();
+ double X_target = InitialConditions->GetIncidentPositionX();
+ double Y_target = InitialConditions->GetIncidentPositionY();
+ double Z_target = InitialConditions->GetIncidentPositionZ();
- TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
- TVector3 HitDirection = PositionOnLassa-PositionOnTarget;
- TVector3 HitDirectionUnit = HitDirection.Unit();
+ TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
+ TVector3 HitDirection = PositionOnLassa-PositionOnTarget;
+ TVector3 HitDirectionUnit = HitDirection.Unit();
- TVector3 BeamDirection = Initial->GetBeamDirection();
- double XBeam = BeamDirection.X();
- double YBeam = BeamDirection.Y();
- double ZBeam = BeamDirection.Z();
+ //TVector3 BeamDirection = InitialConditions->GetBeamDirection();
+ //double XBeam = BeamDirection.X();
+ //double YBeam = BeamDirection.Y();
+ //double ZBeam = BeamDirection.Z();
- ThetaLab = BeamDirection.Angle(HitDirection);
- ThetaLab = ThetaLab/deg;
+ //ThetaLab = BeamDirection.Angle(HitDirection);
+ ThetaLab = 0;//ThetaLab/deg;
- PhiLab = HitDirection.Phi()/deg;
+ PhiLab = HitDirection.Phi()/deg;
- E_ThickSi = Lassa->ThickSi_E[0];
-
- hDetecTheta->Fill(ThetaLab);
-
- hDetecThetaVsPhi->Fill(ThetaLab,PhiLab);
-
- if(Lassa->CsI_E.size()>=1){
-
- detectedEvents++;
- detectedEvs[XStrip][YStrip]++;
-
- E_CsI_Temp = Lassa->CsI_E[0];
-
- if(E_CsI_Temp > E_nucThres[species]){
-
- zeroth_thres = .01;
- first_thres = .05;
- x_0 = -1;
- x_1 = -1;
- x_neg1 = -1;
- zeroth_upperbound = -1.;
- first_lowerbound = -1.;
- slope = 10000;
+ E_ThickSi = Lassa->ThickSi_E[0];
- // We want to check to see if the energy is within 2% of any of the stored energy values
-
- for(int i = 0; i < EnergyRange[species][1].size(); i++){
-
- zeroth_upperbound = (1+zeroth_thres)*EnergyRange[species][1][i];
- zeroth_lowerbound = (1-zeroth_thres)*EnergyRange[species][1][i];
-
- first_upperbound = (1+first_thres)*EnergyRange[species][1][i];
- first_lowerbound = (1-first_thres)*EnergyRange[species][1][i];
-
- if(zeroth_upperbound > E_CsI_Temp && E_CsI_Temp > zeroth_lowerbound) {x_0 = i;}
-
- if(first_upperbound > E_CsI_Temp && E_CsI_Temp > zeroth_upperbound){x_1 = i;}
-
- if(first_lowerbound < E_CsI_Temp && E_CsI_Temp < zeroth_lowerbound){x_neg1 = i;}
-
- if(x_0 != -1 && x_1 != -1 && x_neg1 != -1){break;}
-
- }
-
- // If there is no matching point, then a new (energy,range) point
- // will be added to the vector.
+ if(Lassa->CsI_E.size()>=1){
- if(x_0 == -1){
-
- Range = EnergyLoss_CsI[species].EvaluateMaterialThickness(.0000001*MeV,Lassa->CsI_E[0]*MeV,15*cm,10*micrometer);
- Range_0 = Range;
- Range_1 = Range;
- // Add a new point at the end
- EnergyRange[species][1].push_back(E_CsI_Temp);
-
- EnergyRange[species][2].push_back(Range);
- }
-
- else{
-
- if(x_1 == -1 || x_neg1 == -1){Range_1 = EnergyRange[species][2][x_0];}
-
- else{
-
- if(E_CsI_Temp > EnergyRange[species][1][x_0]){
- slope = (EnergyRange[species][2][x_1]-EnergyRange[species][2][x_0])/(EnergyRange[species][1][x_1]-EnergyRange[species][1][x_0]);
- }
-
- else {
- slope = (EnergyRange[species][2][x_0]-EnergyRange[species][2][x_neg1])/(EnergyRange[species][1][x_0]-EnergyRange[species][1][x_neg1]);
- }
-
- Range_0 = EnergyRange[species][2][x_0];
-
- Range_1 = slope*(E_CsI_Temp-EnergyRange[species][1][x_0]) + Range_0;
- }
+ E_CsI = Lassa->CsI_E[0];
- }
-
-// Range = EnergyLoss_CsI[species].EvaluateMaterialThickness(.0000001*MeV,Lassa->CsI_E[0]*MeV,15*cm,10*micrometer);
- Range = Range_1;
-
- eval = ReactionLoss[species]->Eval(Range/10);
- Random_value = Rand.Uniform(0,1);
-
- if(Random_value>eval)E_CsI = Lassa->CsI_E[0];
- else E_CsI = Rand.Uniform(0,Lassa->CsI_E[0]);
-
- }
-
- else {E_CsI = E_CsI_Temp;}
-
- //E_CsI = Lassa->CsI_E[0];
-
- ELab = E_ThickSi + E_CsI;
+ ELab = E_ThickSi + E_CsI;
- //Check if the energy of the event is within the peak of the beam energy
-
- if(ELab > (totalEnergy-EDelta) && ELab < (totalEnergy+EDelta) && ELab > 0.){
- peakEvents++;
- peakEvs[XStrip][YStrip]++;
- }
- }
- }
+ if(ELab > (BeamEnergy-EDelta) && ELab < (BeamEnergy+EDelta) && ELab > 0.){
+ peakEvents++;
+ }
+ }
+ }
}
-
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
@@ -341,99 +116,6 @@ void Analysis::End(){
peakEff = 100*((double)peakEvents)/((double)detectedEvents);
- hist_Title = "";
- canvas_Title = "";
-
- maxEnergy_String.str("");
-
- maxEnergy_String << maxEnergy;
-
- minEnergy_String.str("");
-
- minEnergy_String << minEnergy;
-
- cout << endl;
- cout << "Max Energy: " << maxEnergy << endl;
- cout << "Min Energy: " << minEnergy << endl;
-
- if(minEnergy == maxEnergy){
-
- hist_Title += speciesNameArray[species];
- hist_Title += " ";
- hist_Title += maxEnergy_String.str();
- hist_Title += " MeV";
-
- canvas_Title += "_";
- canvas_Title += speciesNameArray[species];
- canvas_Title += maxEnergy_String.str();
- canvas_Title += "MeV";
- }
-
- else{
-
- cout << endl;
- cout << "Else!" << endl;
-
- hist_Title += speciesNameArray[species];
- hist_Title += " ";
- hist_Title += minEnergy_String.str();
- hist_Title += "->";
- hist_Title += maxEnergy_String.str();
- hist_Title += " MeV";
-
- canvas_Title += "_";
- canvas_Title += speciesNameArray[species];
- canvas_Title += minEnergy_String.str();
- canvas_Title += "to";
- canvas_Title += maxEnergy_String.str();
- canvas_Title += "MeV";
- }
-
- // Files
-
- simData_Title = "simData";
- simData_Title += canvas_Title;
- simData_Title += ".csv";
- simData.open(simData_Title.c_str());
-
- peakEffMatrix_Title = "peakEffMatrix";
- peakEffMatrix_Title += canvas_Title;
- peakEffMatrix_Title += ".csv";
- peakEffMatrix.open(peakEffMatrix_Title.c_str());
-
- peakEffAvgMatrix_Title = "peakEffAvgMatrix";
- peakEffAvgMatrix_Title += canvas_Title;
- peakEffAvgMatrix_Title += ".csv";
- peakEffAvgMatrix.open(peakEffAvgMatrix_Title.c_str());
-
- detectEvMatrix_Title = "detectEvMatrix";
- detectEvMatrix_Title += canvas_Title;
- detectEvMatrix_Title += ".csv";
- detectEvMatrix.open(detectEvMatrix_Title.c_str());
-
- // Fill the histograms and files
-
- for(int i = 0; i < 16; i++){
- for(int j = 0; j < 16; j++){
-
- if(detectedEvs[i][j] != 0){
- peakEffs[i][j] = 100*((double)peakEvs[i][j])/((double)detectedEvs[i][j]);
- h_peakEff->Fill(i,j,peakEffs[i][j]);
- }
-
- else{
- peakEffs[i][j] = 0;
- }
- }
- }
-
- for(int i = 0; i < 8; i++){
- for(int j = 0; j < 8; j++){
- peakEffsAvg[i][j] = .25*(peakEffs[i][j] + peakEffs[i+8][j] + peakEffs[i][j+8] + peakEffs[i+8][j+8]);
- h_peakEffAvg->Fill(i,j,peakEffsAvg[i][j]);
- }
- }
-
cout << endl;
cout << "Total Events: " << totalEvents << endl;
cout << "Detected Events: " << detectedEvents << endl;
@@ -442,247 +124,38 @@ void Analysis::End(){
cout << "Geometric Efficiency: " << geomEff << endl;
cout << "Peak Efficiency: " << peakEff << endl;
- cout << "Range List Length: " << EnergyRange[species][2].size() << endl;
-
- simData << "Species Name: " << speciesNameArray[species] << endl;
- simData << "Maximum Energy: " << maxEnergy << endl;
- simData << endl << endl;
-
- simData << "Total Events: " << totalEvents << endl;
- simData << "Detected Events: " << detectedEvents << endl;
- simData << "PeakEvents: " << peakEvents << endl;
-
- simData << "Geometric Efficiency: " << geomEff << endl;
- simData << "Peak Efficiency: " << peakEff << endl;
-
- for(int i = 0; i < 16; i++){
- for(int j = 0; j < 16; j++){
-
- if(j == 0) peakEffMatrix << peakEffs[i][j] << ",";
- else{peakEffMatrix<< peakEffs[i][j] << ",";}
- if(j == 15) peakEffMatrix << endl;
-
- if(j == 0) detectEvMatrix << detectedEvs[i][j] << ",";
- else{detectEvMatrix << detectedEvs[i][j] << ",";}
- if(j == 15) detectEvMatrix << endl;
-
- if(i < 8 && j < 8){
-
- if(j == 0) peakEffAvgMatrix << peakEffsAvg[i][j] << ",";
- else{peakEffAvgMatrix<< peakEffsAvg[i][j] << ",";}
- if(j == 15) peakEffAvgMatrix << endl;
-
- }
- }
- }
-
-
- // I want to calculate the mean and the standard deviation of
- // the peak efficiency in two regions (the meat of the crystal
- // and the edges).
-
- // We will be averaging over the four quadrants of one detector.
-
- for(int i = 0; i < 16; i++){
- for(int j = 0; j < 16; j++){
-
- if(0 < i && i < 7 && 0 < j && j < 7 || 0 < i && i < 7 && 8 < j && j < 15 || 8 < i && i < 15 && 0 < j && j < 7 || 8 < i && i < 15 && 8 < j && j < 15){
- peakEffSum_meat += peakEffs[i][j];
- peakEffSquaredSum_meat += peakEffs[i][j]*peakEffs[i][j];
- N_meat++;
-
- cout << endl;
- cout << "x: " << i << " " << "y: " << j << endl;
- cout << "Meat!" << endl;
- cout << "Peak Efficiency: " << peakEffs[i][j] << endl;
-
- }
-
- else{
- peakEffSum_edge += peakEffs[i][j];
- peakEffSquaredSum_edge += peakEffs[i][j]*peakEffs[i][j];
- N_edge++;
-
- cout << endl;
- cout << "x: " << i << " " << "y: " << j << endl;
- cout << "Edge!" << endl;
- cout << "Peak Efficiency: " << peakEffs[i][j] << endl;
- }
- }
- }
-
-
- peakEffAvg_meat = peakEffSum_meat/N_meat;
- peakEffAvg_edge = peakEffSum_edge/N_edge;
-
- peakEffSquaredAvg_meat = peakEffSquaredSum_meat/N_meat;
- peakEffSquaredAvg_edge = peakEffSquaredSum_edge/N_edge;
-
- sigma_meat = sqrt(peakEffSquaredAvg_meat-peakEffAvg_meat*peakEffAvg_meat);
- sigma_edge = sqrt(peakEffSquaredAvg_edge-peakEffAvg_edge*peakEffAvg_edge);
-
- cout << endl;
- cout << "N_meat: " << N_meat << endl;
- cout << "N_edge: " << N_edge << endl;
-
- cout << "Peak Eff Avg Meat: " << peakEffAvg_meat << endl;
- cout << "Peak Eff Avg Edge: " << peakEffAvg_edge << endl;
-
- cout << "sigma_meat: " << sigma_meat << endl;
- cout << "sigma_edge: " << sigma_edge << endl;
-
- // Close the files
-
- detectEvMatrix.close();
- peakEffMatrix.close();
- peakEffAvgMatrix.close();
- simData.close();
-
- // Draw the histograms on canvases
-
- // Canvases
-
- // Peak Efficiencies
-
- canvas_peakEffMatrix = new TCanvas("Peak Eff Matrix","Peak Eff Matrix",1000,1000);
- canvas_peakEffMatrix->cd();
- h_peakEff->Draw("colz");
- h_peakEff->SetStats(0);
- h_peakEff->SetMaximum(100);
- h_peakEff->SetMinimum(30);
- peakEff_hTitle = "Peak Eff Matrix ";
- peakEff_hTitle += hist_Title;
- h_peakEff->SetTitle(peakEff_hTitle.c_str());
- peakEff_cTitle = "canvas_peakEffMatrix";
- peakEff_cTitle += canvas_Title;
- peakEff_cTitle += ".jpg";
- canvas_peakEffMatrix->SaveAs(peakEff_cTitle.c_str());
-
- canvas_peakEffAvgMatrix = new TCanvas("Peak Eff Average Matrix","Peak Eff Average",1000,1000);
- canvas_peakEffAvgMatrix->cd();
- h_peakEffAvg->Draw("colz");
- h_peakEffAvg->SetStats(0);
- h_peakEffAvg->SetMaximum(100);
- h_peakEffAvg->SetMinimum(30);
- peakEffAvg_hTitle = "Peak Eff Average Matrix";
- peakEffAvg_hTitle += hist_Title;
- h_peakEff->SetTitle(peakEffAvg_hTitle.c_str());
- peakEffAvg_cTitle = "canvas_peakEffAvgMatrix";
- peakEffAvg_cTitle += canvas_Title;
- peakEffAvg_cTitle += ".jpg";
- canvas_peakEffAvgMatrix->SaveAs(peakEffAvg_cTitle.c_str());
-
- // Geometric Efficiency
-
- canvas_Emitt = new TCanvas("canvas_Emitt", "Distrib",900,900);
- canvas_Emitt->cd();
- hEmittTheta->Draw("");
- hDetecTheta->SetLineColor(kGreen);
- hDetecTheta->Draw("same");
-
- canvas_Detec = new TCanvas("canvas_Detec", "Distrib",900,900);
- canvas_Detec->cd();
- hDetecTheta->Draw("");
-
- canvas_DetecThetaVsPhi = new TCanvas("canvas_DetecThetaVsPhi", "Distrib",900,900);
- canvas_DetecThetaVsPhi->cd();
- hDetecThetaVsPhi->Draw("colz");
-
- canvas_GeometricEfficiency = new TCanvas("canvas_GeometricEfficiency", "Distrib",900,900);
- canvas_GeometricEfficiency->cd();
- hGeometricEfficiency->Divide(hDetecTheta, hEmittTheta, 100, 1);
- hGeometricEfficiency->Draw("");
-
- // Save canvases to jpg's
-
- emitt_hTitle = "Emitted Angle";
- emitt_hTitle += hist_Title;
- hEmittTheta->SetTitle(emitt_hTitle.c_str());
-
- detec_hTitle = "Detected Angle";
- detec_hTitle += hist_Title;
- hDetecTheta->SetTitle(detec_hTitle.c_str());
-
- detec_hTitle = "Detected Angle (Theta Vs Phi)";
- detec_hTitle += hist_Title;
- hDetecThetaVsPhi->SetTitle(detec_hTitle.c_str());
-
- geomEff_hTitle = "Detected Angle (Theta Vs Phi)";
- geomEff_hTitle += hist_Title;
- hGeometricEfficiency->SetTitle(geomEff_hTitle.c_str());
-
- emitt_cTitle = "canvas_Emitt";
- emitt_cTitle += canvas_Title;
- emitt_cTitle += ".jpg";
- canvas_Emitt->SaveAs(emitt_cTitle.c_str());
-
- detec_cTitle = "canvas_Detec";
- detec_cTitle += canvas_Title;
- detec_cTitle += ".jpg";
- canvas_Detec->SaveAs(detec_cTitle.c_str());
-
- detecThetaVsPhi_cTitle = "canvas_DetecThetaVsPhi";
- detecThetaVsPhi_cTitle += canvas_Title;
- detecThetaVsPhi_cTitle += ".jpg";
- canvas_DetecThetaVsPhi->SaveAs(detecThetaVsPhi_cTitle.c_str());
-
- geomEff_cTitle = "canvas_GeometricEfficiency";
- geomEff_cTitle += canvas_Title;
- geomEff_cTitle += ".jpg";
- canvas_GeometricEfficiency->SaveAs(geomEff_cTitle.c_str());
-
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
- //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
- RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
- RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
- RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
-// RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
-// RootOutput::getInstance()->GetTree()->Branch("totalEvents",&totalEvents,"totalEvents/I");
-// RootOutput::getInstance()->GetTree()->Branch("detectedEvents",&detectedEvents,"detectedEvents/I");
-// RootOutput::getInstance()->GetTree()->Branch("peakEvents",&peakEvents,"peakEvents/I");
+ //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
+ // RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+ // RootOutput::getInstance()->GetTree()->Branch("totalEvents",&totalEvents,"totalEvents/I");
+ // RootOutput::getInstance()->GetTree()->Branch("detectedEvents",&detectedEvents,"detectedEvents/I");
+ // RootOutput::getInstance()->GetTree()->Branch("peakEvents",&peakEvents,"peakEvents/I");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
- RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
-
- E_ThickSi = -1.;
- E_CsI = -1.;
- ELab = -1.;
- E_CsI_Temp = -1.;
- ThetaLab = -1000;
- PhiLab = -1000;
- X = -1000;
- Y = -1000;
- Z = -1000;
- TelescopeNumber = -1;
- EDelta = -1.;
-
-
- ThetaLab_Init = -1000;
-
- fitFunction = "";
- Mu.str("");
- Beta.str("");
- species = -1;
-
- Range = -1.;
- Range_0 = -1.;
- Range_1 = -1.;
- devRange_0 = -1000.;
- devRange_1 = -1000.;
- eval = -1.;
- Random_value = -1.;
-
+ E_ThickSi = 0.;
+ E_CsI = 0.;
+ ELab = 0.;
+ ThetaLab = -1000;
+ PhiLab = -1000;
+ X = -1000;
+ Y = -1000;
+ Z = -1000;
+ TelescopeNumber = -1;
}
////////////////////////////////////////////////////////////////////////////////
@@ -696,13 +169,12 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy{
- public:
- proxy(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
-};
-
-proxy p;
+ class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
+ proxy p;
}
diff --git a/Projects/Lassa/Analysis.h b/Projects/Lassa/Analysis.h
index 33c2590c6015a10ce2dbf502b112c41e0b445f23..a4de214c71bde1f5f8a55d542b94251cb5127317 100644
--- a/Projects/Lassa/Analysis.h
+++ b/Projects/Lassa/Analysis.h
@@ -21,24 +21,12 @@
* *
* *
*****************************************************************************/
-#include"NPVAnalysis.h"
-#include"TLassaPhysics.h"
+#include "NPVAnalysis.h"
+#include "TLassaPhysics.h"
#include "TInitialConditions.h"
#include "NPEnergyLoss.h"
#include "NPReaction.h"
#include "TRandom3.h"
-#include "TMath.h"
-#include <sstream>
-#include <string>
-#include "TCanvas.h"
-#include "TH1F.h"
-#include "TH2F.h"
-#include "TF1.h"
-#include <iostream>
-#include <fstream>
-#include <vector>
-
-
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -60,153 +48,21 @@ class Analysis: public NPL::VAnalysis{
double EDelta;
double PhiLab;
double ThetaLab;
- double ThetaLab_Init;
double X,Y,Z;
double TelescopeNumber;
double thresholdEnergy;
- double totalEnergy;
- double maxEnergy;
- double minEnergy;
- ostringstream maxEnergy_String;
- ostringstream minEnergy_String;
- double E_CsI_Temp;
- double devRange_0;
- double devRange_1;
- double devRange_Sum_0;
- double devRange_Sum_1;
- double devRangeSquared_Sum_0;
- double devRangeSquared_Sum_1;
- double devRange_Avg_0;
- double devRangeSquared_Avg_0;
- double devRange_Avg_1;
- double devRangeSquared_Avg_1;
- double sigma_Range_0;
- double sigma_Range_1;
- double devRangeMax_0;
- double devRangeMax_1;
-
- string hist_Title;
- string canvas_Title;
-
- int XStrip;
- int YStrip;
-
- TLassaPhysics* Lassa;
- TInitialConditions* Initial;
- string speciesName;
- // intermediate variable
- TRandom3 Rand;
-
- // Geometric Efficiency
- double geomEff;
- TH1F* hEmittTheta;
- TH1F* hDetecTheta;
- TH2F* hDetecThetaVsPhi;
- TH1F* hGeometricEfficiency;
- TCanvas *canvas_Emitt;
- TCanvas *canvas_Detec;
- TCanvas *canvas_DetecThetaVsPhi;
- TCanvas *canvas_GeometricEfficiency;
- string emitt_Title;
- string Detec_Title;
- string DetecThetaVsPhi_Title;
- string GeometricEfficiency_Title;
-
- // Peak Efficiencies
int totalEvents;
int detectedEvents;
int peakEvents;
- double peakEff;
- int detectedEvs[16][16];
- int peakEvs[16][16];
- double peakEffs[16][16];
- double peakEffsAvg[8][8];
- TH2F* h_peakEff;
- TH2F* h_peakEffAvg;
- TCanvas* canvas_peakEffMatrix;
- TCanvas* canvas_peakEffAvgMatrix;
- string peakEffs_Title;
- string peakEffsAvg_Title;
-
- // Files
- ofstream simData;
- ofstream peakEffMatrix;
- ofstream peakEffAvgMatrix;
- ofstream detectEvMatrix;
- string simData_Title;
- string peakEffMatrix_Title;
- string peakEffAvgMatrix_Title;
- string detectEvMatrix_Title;
+ double geomEff;
+ double peakEff;
- // Range calculation
- double zeroth_upperbound;
- double zeroth_lowerbound;
- double first_upperbound;
- double first_lowerbound;
- double slope;
- double zeroth_thres;
- double first_thres;
- double Energy_Temp;
- double Range_0;
- double Range_1;
- int x_0;
- int x_1;
- int x_neg1;
-
- //Histogram and Canvas Titles
- string peakEff_hTitle;
- string peakEff_cTitle;
- string peakEffAvg_hTitle;
- string peakEffAvg_cTitle;
- string emitt_hTitle;
- string emitt_cTitle;
- string detec_hTitle;
- string detec_cTitle;
- string detecThetaVsPhi_hTitle;
- string detecThetaVsPhi_cTitle;
- string geomEff_hTitle;
- string geomEff_cTitle;
-
- // CsI Nuclear Reaction Correction Fit Function
- double Range;
- double maxRange;
- double eval;
- double Random_value;
- TF1* ReactionLoss[5];
- NPL::EnergyLoss EnergyLoss_CsI[5];
- string EnergyLoss_Name;
- string EnergyLoss_Table;
- string fitFunction;
- ostringstream Mu;
- ostringstream Beta;
- int species;
- string speciesNameArray[5] = {"proton","deuteron","triton","He3","alpha"};
- double par[5][2] = {{0.0360046, 0.0921675},
- {0.0473272, 0.0433612},
- {0.0558374, 0.0327295},
- {0.0478917, 0.0475857},
- {0.0518713, 0.036529}};
+ // intermediate variable
+ TRandom3 Rand;
- vector <double> EnergyRange[5][2];
- double E_nucThres[5] = {13.617,7.022,5.106,11.064,8.936};
-// double E_nucThres[5] = {0.,0.,0.,0.,0.};
-
- int N_meat;
- int N_edge;
- double peakEffSum_meat;
- double peakEffSum_edge;
- double peakEffAvg_meat;
- double peakEffAvg_edge;
- double peakEffSquaredSum_meat;
- double peakEffSquaredSum_edge;
- double peakEffSquaredAvg_meat;
- double peakEffSquaredAvg_edge;
- double sigma_meat;
- double sigma_edge;
-
-
-
+ TLassaPhysics* Lassa;
+ TInitialConditions* InitialConditions;
};
#endif
diff --git a/Projects/Lassa/RunToTreat.txt b/Projects/Lassa/RunToTreat.txt
index 82908222e590c8760e20a5e47347834912b6c19c..a50a6ea6708811cf5ff58f7e738e4e9b6aa4d9e9 100755
--- a/Projects/Lassa/RunToTreat.txt
+++ b/Projects/Lassa/RunToTreat.txt
@@ -1,5 +1,5 @@
TTreeName
SimulatedTree
RootFileName
- ../../Outputs/Simulation/Example1.root
+ ../../Outputs/Simulation/dumb.root