rtang/SOLARIS_Analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

simplfy Monitor.C

Browse Source
This commit is contained in:
Ryan Tang 2023-04-06 20:03:08 -04:00
parent fdaa35f28d
commit 8a60ae29de
5 changed files with 183 additions and 217 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
.vscode/settings.json vendored
View File

@ -97,7 +97,29 @@
"__split_buffer": "cpp", "__split_buffer": "cpp",
"__tree": "cpp", "__tree": "cpp",
"bitset": "cpp", "bitset": "cpp",
"stack": "cpp" "stack": "cpp",
"Monitor.C": "cpp",
"DataHoSei.C": "cpp",
"Monitors.C": "cpp",
"__bit_reference": "cpp",
"__bits": "cpp",
"__config": "cpp",
"__debug": "cpp",
"__errc": "cpp",
"__locale": "cpp",
"__mutex_base": "cpp",
"__node_handle": "cpp",
"__nullptr": "cpp",
"__string": "cpp",
"__threading_support": "cpp",
"__tuple": "cpp",
"complex": "cpp",
"forward_list": "cpp",
"ios": "cpp",
"locale": "cpp",
"__verbose_abort": "cpp",
"Monitors_Util.C": "cpp",
"Monitor_Util.C": "cpp"
}, },
"better-comments.multilineComments": true, "better-comments.multilineComments": true,

31
armory/AnalysisLib.h
View File

@ -3,10 +3,13 @@
#include <cstdio> #include <cstdio>
#include <vector> #include <vector>
#include <fstream>
#include <string> #include <string>
#include <TMacro.h> #include <TMacro.h>
#include <TList.h> #include <TList.h>
#include <TObjArray.h>
#include <TCutG.h>
namespace AnalysisLib { namespace AnalysisLib {
@ -312,7 +315,7 @@ void PrintReactionConfig(ReactionConfig reaction){
DetGeo detGeo; DetGeo detGeo;
ReactionConfig reactionConfig; ReactionConfig reactionConfig;
void LoadDetGeoAndReactionConfigFile(string detGeoFileName = "detectorGeo.txt", string reactionConfigFileName = "reactionConfig.txt"){ void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt", std::string reactionConfigFileName = "reactionConfig.txt"){
printf("=======================\n"); printf("=======================\n");
printf(" loading detector geometery : %s.", detGeoFileName.c_str()); printf(" loading detector geometery : %s.", detGeoFileName.c_str());
TMacro * haha = new TMacro(); TMacro * haha = new TMacro();
@ -352,7 +355,7 @@ std::vector<std::vector<float>> rdtCorr; // correction of rdt, ch -> MeV
void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.dat"){ void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.dat"){
printf(" loading xf-xn correction."); printf(" loading xf-xn correction.");
xnCorr.clear(); xnCorr.clear();
ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a; float a;
@ -370,7 +373,7 @@ void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.
void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_scaleX.dat"){ void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_scaleX.dat"){
printf(" loading x-Scale correction."); printf(" loading x-Scale correction.");
xScale.clear(); xScale.clear();
ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
@ -387,7 +390,7 @@ void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_sc
void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xfxn_e.dat"){ void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xfxn_e.dat"){
printf(" loading xf/xn-e correction."); printf(" loading xf/xn-e correction.");
xfxneCorr.clear(); xfxneCorr.clear();
ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
@ -404,7 +407,7 @@ void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xf
void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat"){ void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat"){
printf(" loading e correction."); printf(" loading e correction.");
eCorr.clear(); eCorr.clear();
ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
@ -421,7 +424,7 @@ void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat")
void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.dat"){ void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.dat"){
printf(" loading rdt correction."); printf(" loading rdt correction.");
rdtCorr.clear(); rdtCorr.clear();
ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
@ -460,11 +463,11 @@ void LoadReactionParas(bool verbose = false){
// printf("########################## transfer.root updated\n"); // printf("########################## transfer.root updated\n");
//} //}
printf(" loading reaction parameters"); printf(" loading reaction parameters");
ifstream file; std::ifstream file;
file.open("reaction.dat"); file.open("reaction.dat");
hasReactionPara = false; hasReactionPara = false;
if( file.is_open() ){ if( file.is_open() ){
string x; std::string x;
int i = 0; int i = 0;
while( file >> x ){ while( file >> x ){
if( x.substr(0,2) == "//" ) continue; if( x.substr(0,2) == "//" ) continue;
@ -506,6 +509,9 @@ void LoadReactionParas(bool verbose = false){
std::vector<double> CalExTheta(double e, double z){ std::vector<double> CalExTheta(double e, double z){
if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()}; if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
double Ex = TMath::QuietNaN();
double thetaCM = TMath::QuietNaN();
double y = e + mass; // to give the KE + mass of proton; double y = e + mass; // to give the KE + mass of proton;
double Z = alpha * gamm * betRel * z; double Z = alpha * gamm * betRel * z;
double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ; double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ;
@ -537,16 +543,9 @@ std::vector<double> CalExTheta(double e, double z){
double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y; double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y;
double hahaha2 = gamm* betRel * momt; double hahaha2 = gamm* betRel * momt;
thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg(); thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg();
}else{
Ex = TMath::QuietNaN();
thetaCM = TMath::QuietNaN();
} }
}else{
Ex = TMath::QuietNaN();
thetaCM = TMath::QuietNaN();
} }
return {Ex, thetaCM}; return {Ex, thetaCM};
} }

13
working/Mapping.h
View File

@ -11,6 +11,9 @@
//^ Array-xf : 100 - 199 //^ Array-xf : 100 - 199
//^ Array-xn : 200 - 299 //^ Array-xn : 200 - 299
//^ Recoil : 300 - 399 //^ Recoil : 300 - 399
//^ Enum : 400 - 499
//^ EZERO : 500 - 599
//^ Apollo : 600 - 699
//^ //^
//^ line comment is '//^' //^ line comment is '//^'
//^ //^
@ -21,12 +24,12 @@
namespace mapping{ namespace mapping{
const std::vector<std::string> detTypeName = { "e", "xf", "xn", "rdt"}; //C= The comment "//C=" is an indicator DON't Remove const std::vector<std::string> detTypeName = { "e", "xf", "xn", "rdt", "eNum"}; //C= The comment "//C=" is an indicator DON't Remove
const std::vector<int> detGroupID = { 0, 0, 0, 1}; //C& The comment "//C&" is an indicator DON't Remove const std::vector<int> detGroupID = { 0, 0, 0, 1, 2}; //C& The comment "//C&" is an indicator DON't Remove
const std::vector<int> detMaxID = { 100, 200, 300, 400}; //C# The comment "//C#" is an indicator DON't Remove const std::vector<int> detMaxID = { 100, 200, 300, 400, 500}; //C# The comment "//C#" is an indicator DON't Remove
const std::vector<int> detParity = { 1, 1, 1, 1}; const std::vector<int> detParity = { 1, 1, 1, 1, 1};
const std::vector<std::string> groupName = { "Array", "Recoil"}; //C% The comment "//C%" is an indicator DON't Remove const std::vector<std::string> groupName = { "Array", "Recoil", "ELUM"}; //C% The comment "//C%" is an indicator DON't Remove
//!The mapping[i] must match as the IP setting in the DAQ //!The mapping[i] must match as the IP setting in the DAQ

233
working/Monitors.C → working/Monitor.C
View File

@ -1,4 +1,4 @@
#define Monitors_cxx #define Monitor_cxx
#include <TH2.h> #include <TH2.h>
#include <TH1.h> #include <TH1.h>
@ -86,70 +86,57 @@ TString ezCutFile = "";//"ezCut.root";
* * * *
*******************************************************************/ *******************************************************************/
//======== raw data //======== raw data
TH1F** he; TH1F ** he, ** hxf, ** hxn, * hMultiHit; //basic data
TH1F** hxf; TH2F ** hxfVxn, ** heVxs, ** heVx; // correlation
TH1F** hxn; TH2F * heVID, * hxfVID, * hxnVID; // vs ID
TH2F** hxfVxn;
TH2F** heVxs;
TH1F* hMultiHit;
TH2F** heVx; // e vs (xf-xn)/e
TH2F* heVID;
TH2F* hxfVID;
TH2F* hxnVID;
//====== cal data //====== cal data
TH1F** heCal; TH1F ** heCal;
TH2F** hxfCalVxnCal; TH2F ** hxfCalVxnCal;
TH2F** heVxsCal; // raw e vs xf TH2F ** heVxsCal; // raw e vs xf
TH2F** heCalVxCal; // eCal vs xCal TH2F ** heCalVxCal; // eCal vs xCal
TH2F** heCalVxCalG; // eCal vs xCal TH2F ** heCalVxCalG; // eCal vs xCal
TH2F* heCalID; // e vs detID TH2F * heCalID; // e vs detID
TH2F* heCalVz; TH2F * heCalVz;
TH2F* heCalVzGC; TH2F * heCalVzGC;
TH2F** hecalVzRow; TH2F ** hecalVzRow;
//====== Ex data //====== Ex data
TH1F* hEx; TH1F * hEx;
TH1F** hExi; TH1F ** hExi;
TH2F** hExVxCal; TH2F ** hExVxCal;
TH2F* hExThetaCM; TH2F * hExThetaCM;
TH1F* hExCut1; TH1F * hExCut1;
TH1F* hExCut2; TH1F * hExCut2;
//======= Recoil //======= Recoil
TH2F* hrdtID; TH2F * hrdtID;
TH1F** hrdt; // single recoil TH1F ** hrdt; // single recoil
TH1F** hrdtg; TH1F ** hrdtg;
TH2F** hrdt2D; TH2F ** hrdt2D;
TH2F** hrdt2Dg; TH2F ** hrdt2Dg;
TH2F** hrdt2Dsum; TH2F ** hrdt2Dsum;
TH1F* hrdtRate1; TH1F * hrdtRate1;
TH1F* hrdtRate2; TH1F * hrdtRate2;
//======= multi-Hit //======= multi-Hit
TH2I *hmult; TH2I * hmult;
TH1I *hmultEZ; TH1I * hmultEZ;
TH2I *hArrayRDTMatrix; TH2I * hArrayRDTMatrix;
TH2I *hArrayRDTMatrixG; TH2I * hArrayRDTMatrixG;
//======= ARRAY-RDT time diff //======= ARRAY-RDT time diff
TH1I *htdiff; TH1I * htdiff;
TH1I *htdiffg; TH1I * htdiffg;
/*************************** /***************************
***************************/ ***************************/
double zRange[2] = {-1000, 0}; // zMin, zMax double zRange[2] = {-1000, 0}; // zMin, zMax
double Ex, thetaCM;
int padID = 0;
TLatex text; TLatex text;
int numCol, numRow; int numCol, numRow;
@ -159,11 +146,8 @@ Float_t Frac = 0.1; ///Progress bar
TStopwatch StpWatch; TStopwatch StpWatch;
//======= Canvas //======= Canvas
TCanvas *cCanvas; TCanvas * cCanvas;
TString canvasTitle; TString canvasTitle;
int lastRunID;
bool contFlag;
double runTime=0;
//======= Recoil Cut //======= Recoil Cut
TCutG* cutG; //! //general temeprary pointer to cut TCutG* cutG; //! //general temeprary pointer to cut
@ -172,23 +156,19 @@ TObjArray * cutList1;
TObjArray * cutList2; TObjArray * cutList2;
//======= Other Cuts //======= Other Cuts
TCutG* EZCut; TCutG * EZCut;
Bool_t isEZCutFileOpen;
#include "Monitors.h" #include "Monitor.h"
//^########################################################### //^###########################################################
//^ * Begin //^ * Begin
//^########################################################### //^###########################################################
void Monitors::Begin(TTree *tree){ void Monitor::Begin(TTree *tree){
TString option = GetOption(); TString option = GetOption();
NumEntries = tree->GetEntries(); NumEntries = tree->GetEntries();
canvasTitle = GetCanvasTitle(); canvasTitle = GetCanvasTitle();
lastRunID = -1;
contFlag = false;
printf("###########################################################\n"); printf("###########################################################\n");
printf("########## SOLARIS Monitors.C #########\n"); printf("########## SOLARIS Monitors.C #########\n");
@ -229,36 +209,22 @@ void Monitors::Begin(TTree *tree){
gROOT->cd(); gROOT->cd();
he = new TH1F * [mapping::NARRAY]; CreateListOfHist1D(he, mapping::NARRAY, "he", "Raw e (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]);
hxf = new TH1F * [mapping::NARRAY]; CreateListOfHist1D(hxf, mapping::NARRAY, "hxf", "Raw xf (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]);
hxn = new TH1F * [mapping::NARRAY]; CreateListOfHist1D(hxn, mapping::NARRAY, "hxn", "Raw xn (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]);
hxfVxn = new TH2F * [mapping::NARRAY];
heVxs = new TH2F * [mapping::NARRAY];
heVx = new TH2F * [mapping::NARRAY];
heCal = new TH1F * [mapping::NARRAY]; CreateListOfHist2D(hxfVxn, mapping::NARRAY, "hxfVxn", "Raw xf vs. xn (ch=%d);xf (channel);xn (channel)", 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]);
hxfCalVxnCal = new TH2F * [mapping::NARRAY]; CreateListOfHist2D(heVxs, mapping::NARRAY, "heVxs", "Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
heVxsCal = new TH2F * [mapping::NARRAY]; CreateListOfHist2D(heVx, mapping::NARRAY, "heVx", "Raw PSD E vs. X (ch=%d);X (channel);E (channel)", 500, -0.1, 1.1, 500, rawEnergyRange[0], rawEnergyRange[1]);
heCalVxCal = new TH2F * [mapping::NARRAY];
heCalVxCalG = new TH2F * [mapping::NARRAY];
for (Int_t i = 0; i < mapping::NARRAY; i++) {//array loop CreateListOfHist1D(heCal, mapping::NARRAY, "heCal", "Corrected e (ch=%d); e (MeV); count", 2000, energyRange[0], energyRange[1]);
he[i] = new TH1F(Form("he%d", i), Form("Raw e (ch=%d); e (channel); count", i), 2000, rawEnergyRange[0], rawEnergyRange[1]);
hxf[i] = new TH1F(Form("hxf%d", i), Form("Raw xf (ch=%d); xf (channel); count", i), 200, rawEnergyRange[0], rawEnergyRange[1]);
hxn[i] = new TH1F(Form("hxn%d", i), Form("Raw xn (ch=%d); xn (channel); count", i), 200, rawEnergyRange[0], rawEnergyRange[1]);
hxfVxn[i] = new TH2F(Form("hxfVxn%d",i), Form("Raw xf vs. xn (ch=%d);xf (channel);xn (channel)",i), 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]);
heVxs[i] = new TH2F(Form("heVxs%d", i), Form("Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", i), 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
heVx[i] = new TH2F(Form("heVx%d",i), Form("Raw PSD E vs. X (ch=%d);X (channel);E (channel)",i), 500, -0.1, 1.1, 500, rawEnergyRange[0], rawEnergyRange[1]);
heCal[i] = new TH1F(Form("heCal%d", i), Form("Corrected e (ch=%d); e (MeV); count", i), 2000, energyRange[0], energyRange[1]);
hxfCalVxnCal[i] = new TH2F(Form("hxfCalVxnCal%d",i), Form("Corrected XF vs. XN (ch=%d);XF (channel);XN (channel)",i), 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]);
heVxsCal[i] = new TH2F(Form("heVxsCal%d", i), Form("Raw e vs Corrected xf+xn (ch=%d); corrected xf+xn (channel); Raw e (channel)", i), 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
heCalVxCal[i] = new TH2F(Form("heCalVxCal%d",i), Form("Cal PSD E vs. X (ch=%d);X (cm);E (MeV)",i), 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]);
heCalVxCalG[i] = new TH2F(Form("heCalVxCalG%d",i), Form("Cal PSD E vs. X (ch=%d);X (cm);E (MeV)",i), 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]);
}
heVID = new TH2F("heVID", "Raw e vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); CreateListOfHist2D(hxfCalVxnCal, mapping::NARRAY, "hxfCalVxnCal", "Corrected XF vs. XN (ch=%d);XF (channel);XN (channel)", 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]);
CreateListOfHist2D(heVxsCal , mapping::NARRAY, "heVxsCal", "Raw e vs Corrected xf+xn (ch=%d); corrected xf+xn (channel); Raw e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
CreateListOfHist2D(heCalVxCal , mapping::NARRAY, "heCalVxCal", "Cal PSD E vs. X (ch=%d);X (cm);E (MeV)", 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]);
CreateListOfHist2D(heCalVxCalG , mapping::NARRAY, "heCalVxCalG", "Cal PSD E vs. X (ch=%d);X (cm);E (MeV)", 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]);
heVID = new TH2F("heVID", "Raw e vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]);
hxfVID = new TH2F("hxfVID", "Raw xf vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); hxfVID = new TH2F("hxfVID", "Raw xf vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]);
hxnVID = new TH2F("hxnVID", "Raw xn vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); hxnVID = new TH2F("hxnVID", "Raw xn vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]);
@ -269,7 +235,7 @@ void Monitors::Begin(TTree *tree){
//====================== E-Z plot //====================== E-Z plot
heCalVz = new TH2F("heCalVz", "E vs. Z;Z (mm);E (MeV)" , 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]); heCalVz = new TH2F("heCalVz", "E vs. Z;Z (mm);E (MeV)" , 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]);
heCalVzGC = new TH2F("heCalVzGC","E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, 0, energyRange[1]); heCalVzGC = new TH2F("heCalVzGC","E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, 0, energyRange[1]);
hecalVzRow = new TH2F * [numRow]; hecalVzRow = new TH2F * [numRow];
for( int i = 0; i < numRow; i++){ for( int i = 0; i < numRow; i++){
hecalVzRow[i] = new TH2F(Form("heCalVzRow%d", i), Form("E vs. Z (ch=%d-%d); Z (cm); E (MeV)", numCol*i, numCol*(i+1)-1), 500, zRange[0], zRange[1], 500, energyRange[0], energyRange[1]); hecalVzRow[i] = new TH2F(Form("heCalVzRow%d", i), Form("E vs. Z (ch=%d-%d); Z (cm); E (MeV)", numCol*i, numCol*(i+1)-1), 500, zRange[0], zRange[1], 500, energyRange[0], energyRange[1]);
@ -282,13 +248,10 @@ void Monitors::Begin(TTree *tree){
hExCut2 = new TH1F("hExCut2",Form("excitation spectrum w/ goodFlag; Ex [MeV] ; Count / %4.0f keV", exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); hExCut2 = new TH1F("hExCut2",Form("excitation spectrum w/ goodFlag; Ex [MeV] ; Count / %4.0f keV", exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
hExCut1->SetLineColor(2); hExCut1->SetLineColor(2);
hExCut2->SetLineColor(4); hExCut2->SetLineColor(4);
hExi = new TH1F * [mapping::NARRAY]; TString haha = "Ex (det=%i) w/goodFlag; Ex [MeV]; Count / " +std::to_string(exRange[0]) + "keV";
hExVxCal = new TH2F * [mapping::NARRAY]; CreateListOfHist1D(hExi, mapping::NARRAY, "hExi", haha.Data(), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
for(int i = 0 ; i < mapping::NARRAY; i++){ CreateListOfHist2D(hExVxCal, mapping::NARRAY, "hExVxCal", "Ex vs X (ch=%d); X (cm); Ex (MeV)", 500, -0.1, 1.1, (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
hExi[i] = new TH1F(Form("hExi%02d", i), Form("Ex (det=%i) w/goodFlag; Ex [MeV]; Count / %4.0f keV",i, exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
hExVxCal[i] = new TH2F(Form("hExVxCal%d",i), Form("Ex vs X (ch=%d); X (cm); Ex (MeV)", i), 500, -0.1, 1.1, (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
}
hExThetaCM = new TH2F("hExThetaCM", "Ex vs ThetaCM; ThetaCM [deg]; Ex [MeV]", 200, thetaCMRange[0], thetaCMRange[1], (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); hExThetaCM = new TH2F("hExThetaCM", "Ex vs ThetaCM; ThetaCM [deg]; Ex [MeV]", 200, thetaCMRange[0], thetaCMRange[1], (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]);
@ -332,7 +295,6 @@ void Monitors::Begin(TTree *tree){
htdiff = new TH1I("htdiff" ,"Coincident time (recoil-dE - array); time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]); htdiff = new TH1I("htdiff" ,"Coincident time (recoil-dE - array); time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]);
htdiffg = new TH1I("htdiffg","Coincident time (recoil-dE - array) w/ recoil gated; time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]); htdiffg = new TH1I("htdiffg","Coincident time (recoil-dE - array) w/ recoil gated; time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]);
printf("======================================== End of histograms Declaration\n"); printf("======================================== End of histograms Declaration\n");
StpWatch.Start(); StpWatch.Start();
@ -341,7 +303,7 @@ void Monitors::Begin(TTree *tree){
//^########################################################### //^###########################################################
//^ * Process //^ * Process
//^########################################################### //^###########################################################
Bool_t Monitors::Process(Long64_t entry){ Bool_t Monitor::Process(Long64_t entry){
if( entry == 0 ) printf("========== %s \n", __func__); if( entry == 0 ) printf("========== %s \n", __func__);
@ -547,7 +509,7 @@ Bool_t Monitors::Process(Long64_t entry){
if( !isTimeGateOn ) coinFlag = true; if( !isTimeGateOn ) coinFlag = true;
//================ E-Z gate //================ E-Z gate
if( isEZCutFileOpen ) { if( EZCut ) {
if( EZCut->IsInside(z[detID], eCal[detID]) ) ezGate = true; if( EZCut->IsInside(z[detID], eCal[detID]) ) ezGate = true;
}else{ }else{
ezGate = true; ezGate = true;
@ -562,9 +524,9 @@ Bool_t Monitors::Process(Long64_t entry){
isGoodEventFlag = true; isGoodEventFlag = true;
} }
}//end of array loop
if( !isEZCutFileOpen ) ezGate = true; }//end of array loop
if( EZCut == nullptr ) ezGate = true;
//@*********** RECOILS ***********************************************/ //@*********** RECOILS ***********************************************/
for( int i = 0; i < mapping::NRDT ; i++){ for( int i = 0; i < mapping::NRDT ; i++){
@ -594,54 +556,14 @@ Bool_t Monitors::Process(Long64_t entry){
if( eCal[detID] < eCalCut[0] ) continue ; if( eCal[detID] < eCalCut[0] ) continue ;
if( eCal[detID] > eCalCut[1] ) continue ; if( eCal[detID] > eCalCut[1] ) continue ;
double Ex, thetaCM;
if( AnalysisLib::hasReactionPara ){ if( AnalysisLib::hasReactionPara ){
std::vector<double> ExThetaCM = AnalysisLib::CalExTheta(eCal[detID], x[detID]); std::vector<double> ExThetaCM = AnalysisLib::CalExTheta(eCal[detID], x[detID]);
Ex = ExThetaCM[0]; Ex = ExThetaCM[0];
thetaCM = ExThetaCM[1]; thetaCM = ExThetaCM[1];
// ///======== Ex calculation by Ryan
// double y = eCal[detID] + mass; // to give the KE + mass of proton;
// double Z = alpha * gamm * betRel * z[detID];
// double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ;
// if( TMath::Abs(Z) < H ) {
// ///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi)
// double tolerrence = 0.001;
// double phi = 0; ///initial phi = 0 -> ensure the solution has f'(phi) > 0
// double nPhi = 0; /// new phi
// int iter = 0;
// do{
// phi = nPhi;
// nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2));
// iter ++;
// if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break;
// }while( TMath::Abs(phi - nPhi ) > tolerrence);
// phi = nPhi;
// /// check f'(phi) > 0
// double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2);
// if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){
// double K = H * TMath::Sin(phi);
// double x = TMath::ACos( mass / ( y * gamm - K));
// double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame
// double EB = TMath::Sqrt(mass*mass + Et*Et - 2*Et*TMath::Sqrt(momt*momt + mass * mass));
// Ex = EB - massB;
// double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y;
// double hahaha2 = gamm* betRel * momt;
// thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg();
// }else{
// Ex = TMath::QuietNaN();
// thetaCM = TMath::QuietNaN();
// }
// }else{
// Ex = TMath::QuietNaN();
// thetaCM = TMath::QuietNaN();
// }
}else{ }else{
Ex = TMath::QuietNaN(); Ex = TMath::QuietNaN();
thetaCM = TMath::QuietNaN(); thetaCM = TMath::QuietNaN();
@ -674,19 +596,10 @@ Bool_t Monitors::Process(Long64_t entry){
//^########################################################### //^###########################################################
//^ * Terminate //^ * Terminate
//^########################################################### //^###########################################################
void Monitors::Terminate(){ void Monitor::Terminate(){
printf("============================== finishing.\n"); printf("============================== finishing.\n");
gROOT->cd(); gROOT->cd();
int strLen = canvasTitle.Sizeof();
canvasTitle.Remove(strLen-3);
TString runTimeStr = "";
if( runTime > 0. ) {
runTimeStr = Form("%.0f min", runTime);
canvasTitle += " | " + runTimeStr;
}
//############################################ User is free to edit this section //############################################ User is free to edit this section
//--- Canvas Size //--- Canvas Size
@ -867,19 +780,19 @@ void Monitors::Terminate(){
/************************************/ /************************************/
StpWatch.Start(kFALSE); StpWatch.Start(kFALSE);
gROOT->ProcessLine(".L ../Armory/Monitor_Util.C"); //gROOT->ProcessLine(".L ../Armory/Monitor_Util.C"); //TODO some pointer is empty
printf("=============== loaded Monitor_Utils.C\n"); //printf("=============== loaded Monitor_Utils.C\n");
gROOT->ProcessLine(".L ../Armory/AutoFit.C"); gROOT->ProcessLine(".L ../Armory/AutoFit.C");
printf("=============== loaded Armory/AutoFit.C\n"); printf("=============== loaded Armory/AutoFit.C\n");
gROOT->ProcessLine(".L ../Armory/RDTCutCreator.C"); // gROOT->ProcessLine(".L ../Armory/RDTCutCreator.C");
printf("=============== loaded Armory/RDTCutCreator.C\n"); // printf("=============== loaded Armory/RDTCutCreator.C\n");
gROOT->ProcessLine(".L ../Armory/Check_rdtGate.C"); // gROOT->ProcessLine(".L ../Armory/Check_rdtGate.C");
printf("=============== loaded Armory/Check_rdtGate.C\n"); // printf("=============== loaded Armory/Check_rdtGate.C\n");
gROOT->ProcessLine(".L ../Armory/readTrace.C"); // gROOT->ProcessLine(".L ../Armory/readTrace.C");
printf("=============== loaded Armory/readTrace.C\n"); // printf("=============== loaded Armory/readTrace.C\n");
//gROOT->ProcessLine(".L ../Armory/readRawTrace.C"); // gROOT->ProcessLine(".L ../Armory/readRawTrace.C");
//printf("=============== loaded Armory/readRawTrace.C\n"); // printf("=============== loaded Armory/readRawTrace.C\n");
gROOT->ProcessLine("listDraws()"); // gROOT->ProcessLine("listDraws()");
/************************* Save histograms to root file*/ /************************* Save histograms to root file*/

99
working/Monitors.h → working/Monitor.h
View File

@ -1,5 +1,5 @@
#ifndef Monitors_h #ifndef Monitor_h
#define Monitors_h #define Monitor_h
#include <TROOT.h> #include <TROOT.h>
#include <TChain.h> #include <TChain.h>
@ -13,7 +13,7 @@
#include "Mapping.h" #include "Mapping.h"
#include "../armory/AnalysisLib.h" #include "../armory/AnalysisLib.h"
class Monitors : public TSelector { class Monitor : public TSelector {
public : public :
TTree *fChain; //!pointer to the analyzed TTree or TChain TTree *fChain; //!pointer to the analyzed TTree or TChain
@ -38,19 +38,19 @@ public :
TBranch *b_RDTTimestamp; //! TBranch *b_RDTTimestamp; //!
// trace analysis data // trace analysis data
Float_t * we; //! // Float_t * we; //!
Float_t * weR; //! // Float_t * weR; //!
Float_t * weT; //! // Float_t * weT; //!
Float_t * wrdt; //! // Float_t * wrdt; //!
Float_t * wrdtT; //! // Float_t * wrdtT; //!
Float_t * wrdtR; //! // Float_t * wrdtR; //!
TBranch *b_Trace_Energy; //! // TBranch *b_Trace_Energy; //!
TBranch *b_Trace_Energy_RiseTime; //! // TBranch *b_Trace_Energy_RiseTime; //!
TBranch *b_Trace_Energy_Time; //! // TBranch *b_Trace_Energy_Time; //!
TBranch *b_Trace_RDT; //! // TBranch *b_Trace_RDT; //!
TBranch *b_Trace_RDT_Time; //! // TBranch *b_Trace_RDT_Time; //!
TBranch *b_Trace_RDT_RiseTime; //! // TBranch *b_Trace_RDT_RiseTime; //!
bool isArrayTraceExist; bool isArrayTraceExist;
bool isRDTTraceExist; bool isRDTTraceExist;
@ -70,7 +70,11 @@ public :
ULong64_t startTime ; ULong64_t startTime ;
ULong64_t endTime ; ULong64_t endTime ;
Monitors(TTree * /*tree*/ =0) : fChain(0) { int padID;
Monitor(TTree * /*tree*/ =0) : fChain(0) {
printf("------ %s\n", __func__);
e = new Float_t [mapping::NARRAY]; e = new Float_t [mapping::NARRAY];
xf = new Float_t [mapping::NARRAY]; xf = new Float_t [mapping::NARRAY];
@ -88,15 +92,18 @@ public :
xnCal = new float [mapping::NARRAY]; xnCal = new float [mapping::NARRAY];
eCal = new float [mapping::NARRAY]; eCal = new float [mapping::NARRAY];
isXNCorrOK = true; isXNCorrOK = true;
isXFXNCorrOK = true; isXFXNCorrOK = true;
isXScaleCorrOK = true; isXScaleCorrOK = true;
isECorrOK = true; isECorrOK = true;
isRDTCorrOK = true; isRDTCorrOK = true;
padID = 0;
} }
virtual ~Monitors() { virtual ~Monitor() {
printf("------ %s\n", __func__);
delete e ; delete e ;
delete xf ; delete xf ;
@ -114,6 +121,8 @@ public :
delete xnCal; delete xnCal;
delete eCal; delete eCal;
printf("------end of %s\n", __func__);
} }
virtual Int_t Version() const { return 2; } virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree); virtual void Begin(TTree *tree);
@ -144,16 +153,17 @@ public :
void PlotRDT(int id, bool isRaw); void PlotRDT(int id, bool isRaw);
//void PlotCRDTPolar(); //void PlotCRDTPolar();
//template<typename T> T ** CreateListOfHist(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax, int binY = 0, float yMin = 0, float yMax = 0); template<typename T> void CreateListOfHist1D(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax);
template<typename T> void CreateListOfHist2D(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax, int binY, float yMin, float yMax);
ClassDef(Monitors,0); ClassDef(Monitor,0);
}; };
#endif #endif
#ifdef Monitors_cxx #ifdef Monitor_cxx
void Monitors::Init(TTree *tree){ void Monitor::Init(TTree *tree){
printf("========== %s \n", __func__); printf("========== %s \n", __func__);
@ -220,7 +230,7 @@ void Monitors::Init(TTree *tree){
printf("=================================== End of Branch Pointer Inititization. \n"); printf("=================================== End of Branch Pointer Inititization. \n");
} }
Bool_t Monitors::Notify(){ Bool_t Monitor::Notify(){
return kTRUE; return kTRUE;
} }
@ -232,20 +242,39 @@ void DrawLine(TH1 * hist, double pos){
line->Draw(""); line->Draw("");
} }
void Monitor::SlaveBegin(TTree * /*tree*/){
/// not use, if use, place in Monitor.C
void Monitors::SlaveBegin(TTree * /*tree*/){
/// not use, if use, place in Monitors.C
TString option = GetOption(); TString option = GetOption();
} }
void Monitors::SlaveTerminate(){ void Monitor::SlaveTerminate(){
/// not use, if use, place in Monitors.C /// not use, if use, place in Monitor.C
} }
template<typename T> void Monitor::CreateListOfHist1D(T ** &histList,
int size,
const char * namePrefix,
const char * TitleForm,
int binX, float xMin, float xMax){
histList = new T * [size];
for(int i = 0; i < size; i++) histList[i] = new T(Form("%s%d", namePrefix, i), Form(TitleForm, i), binX, xMin, xMax);
}
template<typename T> void Monitor::CreateListOfHist2D(T ** &histList,
int size,
const char * namePrefix,
const char * TitleForm,
int binX, float xMin, float xMax,
int binY, float yMin, float yMax){
histList = new T * [size];
for(int i = 0; i < size; i++) histList[i] = new T(Form("%s%d", namePrefix, i), Form(TitleForm, i), binX, xMin, xMax, binY, yMin, yMax);
}
/*########################################################### /*###########################################################
* Plotting Function * Plotting Function
###########################################################*/ ###########################################################*/
@ -259,7 +288,7 @@ void DrawBox(TH1* hist, double x1, double x2, Color_t color, float alpha){
} }
void Monitors::Draw2DHist(TH2F * hist){ void Monitor::Draw2DHist(TH2F * hist){
if( hist->Integral() < 3000 ){ if( hist->Integral() < 3000 ){
hist->SetMarkerStyle(20); hist->SetMarkerStyle(20);
@ -271,7 +300,7 @@ void Monitors::Draw2DHist(TH2F * hist){
} }
void Monitors::PlotEZ(bool isRaw){ void Monitor::PlotEZ(bool isRaw){
padID++; cCanvas->cd(padID); padID++; cCanvas->cd(padID);
if( isRaw ) { if( isRaw ) {
@ -321,7 +350,7 @@ void Monitors::PlotEZ(bool isRaw){
} }
void Monitors::PlotTDiff(bool isGated, bool isLog){ void Monitor::PlotTDiff(bool isGated, bool isLog){
padID++; cCanvas->cd(padID); padID++; cCanvas->cd(padID);
if( isLog ) cCanvas->cd(padID)->SetLogy(1); if( isLog ) cCanvas->cd(padID)->SetLogy(1);
double yMax = 0; double yMax = 0;
@ -344,7 +373,7 @@ void Monitors::PlotTDiff(bool isGated, bool isLog){
DrawBox(htdiff, timeGate[0], timeGate[1], kGreen, 0.2); DrawBox(htdiff, timeGate[0], timeGate[1], kGreen, 0.2);
} }
void Monitors::PlotRDT(int id, bool isRaw){ void Monitor::PlotRDT(int id, bool isRaw){
padID++; cCanvas->cd(padID); padID++; cCanvas->cd(padID);
if( isRaw ){ if( isRaw ){
@ -359,10 +388,10 @@ void Monitors::PlotRDT(int id, bool isRaw){
} }
//void Monitors::PlotCRDTPolar(){ //void Monitor::PlotCRDTPolar(){
// padID++; cCanvas->cd(padID); // padID++; cCanvas->cd(padID);
// cCanvas->cd(padID)->DrawFrame(-50, -50, 50, 50); // cCanvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// hcrdtPolar->Draw("same colz pol"); // hcrdtPolar->Draw("same colz pol");
//} //}
#endif // #ifdef Monitors_cxx #endif // #ifdef Monitor_cxx