gwm17/Mask
1
0
Fork 0
mirror of https://github.com/gwm17/Mask.git synced 2024-11-22 18:28:51 -05:00
Code Issues Projects Releases Wiki Activity

Some updates to the SABRE detector classes. Minor refinement of RootPlot

Browse Source
This commit is contained in:
Gordon McCann 2022-02-08 11:44:32 -05:00
parent 736094c4e4
commit a59cd18d9a
6 changed files with 202 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
include/RootPlotter.h
View File

@ -5,6 +5,8 @@
#include <string> #include <string>
#include "Nucleus.h" #include "Nucleus.h"
#include "RxnType.h"
#include "MaskFile.h"
#include <THashTable.h> #include <THashTable.h>
#include <TH1.h> #include <TH1.h>
@ -29,6 +31,8 @@ public:
return table; return table;
}; };
void FillData(const Mask::Nucleus& nuc, double detKE = 0.0, const std::string& modifier = ""); void FillData(const Mask::Nucleus& nuc, double detKE = 0.0, const std::string& modifier = "");
void FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type);
void FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type);
private: private:
THashTable* table; THashTable* table;

2
include/SabreEfficiency.h
View File

@ -1,6 +1,7 @@
#ifndef SABREEFFICIENCY_H #ifndef SABREEFFICIENCY_H
#define SABREEFFICIENCY_H #define SABREEFFICIENCY_H
#include "MaskFile.h"
#include "DetectorEfficiency.h" #include "DetectorEfficiency.h"
#include "SabreDetector.h" #include "SabreDetector.h"
#include "Target.h" #include "Target.h"
@ -18,6 +19,7 @@ public:
private: private:
std::pair<bool,double> IsSabre(Mask::Nucleus& nucleus); std::pair<bool,double> IsSabre(Mask::Nucleus& nucleus);
void CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type);
std::vector<SabreDetector> detectors; std::vector<SabreDetector> detectors;

8
src/Detectors/DetectorEfficiency.cpp
View File

@ -15,15 +15,16 @@ int main(int argc, char** argv) {
std::string outputname = argv[2]; std::string outputname = argv[2];
std::string statsname = argv[3]; std::string statsname = argv[3];
/*
SabreEfficiency sabre; SabreEfficiency sabre;
std::string mapfile = "./etc/SabreDeadChannels.txt"; std::string mapfile = "./etc/SabreDeadChannels.txt";
sabre.SetDeadChannelMap(mapfile); sabre.SetDeadChannelMap(mapfile);
sabre.CalculateEfficiency(inputname, outputname, statsname); sabre.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl; //std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt"); //sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
*/
/*
try { try {
AnasenEfficiency anasen; AnasenEfficiency anasen;
std::string mapfile = "./etc/AnasenDeadChannels.txt"; std::string mapfile = "./etc/AnasenDeadChannels.txt";
@ -36,5 +37,6 @@ int main(int argc, char** argv) {
std::cerr<<"Terminating."<<std::endl; std::cerr<<"Terminating."<<std::endl;
return 1; return 1;
} }
*/
return 0; return 0;
} }

135
src/Detectors/SabreEfficiency.cpp
View File

@ -39,7 +39,7 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read); Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
Mask::MaskFile output(outputname, Mask::MaskFile::FileType::read); Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
std::ofstream stats(statsname); std::ofstream stats(statsname);
stats<<std::setprecision(5); stats<<std::setprecision(5);
@ -52,18 +52,23 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
Mask::Nucleus nucleus; Mask::Nucleus nucleus;
std::vector<int> counts; std::vector<int> counts;
std::vector<int> coinc_counts;
switch(header.rxn_type) { switch(header.rxn_type) {
case Mask::RxnType::PureDecay: case Mask::RxnType::PureDecay:
counts.resize(3, 0); counts.resize(3, 0);
coinc_counts.resize(1, 0);
break; break;
case Mask::RxnType::OneStepRxn: case Mask::RxnType::OneStepRxn:
counts.resize(4, 0); counts.resize(4, 0);
coinc_counts.resize(1, 0);
break; break;
case Mask::RxnType::TwoStepRxn: case Mask::RxnType::TwoStepRxn:
counts.resize(6, 0); counts.resize(6, 0);
coinc_counts.resize(4, 0);
break; break;
case Mask::RxnType::ThreeStepRxn: case Mask::RxnType::ThreeStepRxn:
counts.resize(8, 0); counts.resize(8, 0);
coinc_counts.resize(11, 0);
break; break;
default: default:
{ {
@ -96,12 +101,15 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
auto result = IsSabre(nucleus); auto result = IsSabre(nucleus);
if(result.first) { if(result.first) {
data.detect_flag[i] = true; data.detect_flag[i] = true;
data.KE[i] = result.second;
counts[i]++; counts[i]++;
} else if(data.detect_flag[i] == true) { } else if(data.detect_flag[i] == true) {
data.detect_flag[i] = false; data.detect_flag[i] = false;
} }
} }
CountCoincidences(data, coinc_counts, header.rxn_type);
output.WriteData(data); output.WriteData(data);
} }
@ -114,6 +122,54 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl; stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl; stats<<"---------------------"<<std::endl;
} }
stats<<"Coincidence Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
if(header.rxn_type == Mask::RxnType::PureDecay)
{
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::OneStepRxn)
{
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::TwoStepRxn)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::ThreeStepRxn)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats.close(); stats.close();
std::cout<<std::endl; std::cout<<std::endl;
@ -206,3 +262,80 @@ std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
return std::make_pair(false,0.0); return std::make_pair(false,0.0);
} }
void SabreEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type) {
if (rxn_type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
{
counts[0]++;
}
else if (rxn_type == Mask::RxnType::OneStepRxn && data.detect_flag[2] && data.detect_flag[3])
{
counts[0]++;
}
else if(rxn_type == Mask::RxnType::TwoStepRxn)
{
if(data.detect_flag[2] && data.detect_flag[4])
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[5])
{
counts[1]++;
}
if(data.detect_flag[4] && data.detect_flag[5])
{
counts[2]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5])
{
counts[3]++;
}
}
else if(rxn_type == Mask::RxnType::ThreeStepRxn)
{
if(data.detect_flag[2] && data.detect_flag[4])
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[6])
{
counts[1]++;
}
if(data.detect_flag[2] && data.detect_flag[7])
{
counts[2]++;
}
if(data.detect_flag[4] && data.detect_flag[6])
{
counts[3]++;
}
if(data.detect_flag[4] && data.detect_flag[7])
{
counts[4]++;
}
if(data.detect_flag[6] && data.detect_flag[7])
{
counts[5]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6])
{
counts[6]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7])
{
counts[7]++;
}
if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7])
{
counts[8]++;
}
if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
{
counts[9]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
{
counts[10]++;
}
}
}

2
src/MaskApp.cpp
View File

@ -76,7 +76,7 @@ namespace Mask {
case RxnType::ThreeStepRxn: case RxnType::ThreeStepRxn:
{ {
m_sys = new ThreeStepSystem(); m_sys = new ThreeStepSystem();
m_rxn_type = RxnType::TwoStepRxn; m_rxn_type = RxnType::ThreeStepRxn;
for(int i=0; i<5; i++) { for(int i=0; i<5; i++) {
input>>z>>a; input>>z>>a;
avec.push_back(a); avec.push_back(a);

57
src/Plotters/ROOT/RootPlotter.cpp
View File

@ -1,5 +1,4 @@
#include "RootPlotter.h" #include "RootPlotter.h"
#include "MaskFile.h"
#include <TFile.h> #include <TFile.h>
#include <iostream> #include <iostream>
@ -39,6 +38,60 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::st
} }
void RootPlotter::FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type)
{
std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor";
std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
if(type == Mask::RxnType::PureDecay)
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
}
else
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
}
if(type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn)
{
std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor";
std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
}
if(type == Mask::RxnType::ThreeStepRxn)
{
std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor";
std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
}
}
void RootPlotter::FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type)
{
std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor_detected";
std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
if(type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
}
else if(data.detect_flag[2] && data.detect_flag[3])
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
}
if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4] && data.detect_flag[5])
{
std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor_detected";
std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
}
if(type == Mask::RxnType::ThreeStepRxn && data.detect_flag[6] && data.detect_flag[7])
{
std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor_detected";
std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
}
}
void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val) { void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val) {
TH1F* h = (TH1F*) table->FindObject(name.c_str()); TH1F* h = (TH1F*) table->FindObject(name.c_str());
if(h) { if(h) {
@ -224,6 +277,8 @@ int main(int argc, char** argv) {
plotter.FillData(nucleus, data.KE[i], "detected"); plotter.FillData(nucleus, data.KE[i], "detected");
} }
} }
plotter.FillCorrelations(data, header.rxn_type);
plotter.FillCorrelationsDetected(data, header.rxn_type);
count++; count++;
} }
std::cout<<std::endl; std::cout<<std::endl;