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

Compare commits

base: rtang:5bb5f5895653d540dd830552f8ab07cbd8eed4ee
Branches Tags
rtang:Master
..
compare: rtang:19a567f8fc1e31597cba4143bf4bc156d914130a
Branches Tags
rtang:Master

4 Commits
5bb5f58956 ... 19a567f8fc

Author SHA1 Message Date
Ryan@Home 19a567f8fc adding RDT Cuts... not tested 2024-07-09 19:13:04 -04:00
Ryan@Home 9c2f52009a MonAnalyzer: finished array. need to do RDT 2024-07-09 17:13:28 -04:00
Ryan@Home a80e5e2b64 making progress... need to convert on of helios run to have better test 2024-07-08 22:24:44 -04:00
Ryan@Home 56285cef62 still working on Monitor. Shift to TTreeReader 2024-07-08 18:04:25 -04:00
10 changed files with 808 additions and 245 deletions
Show Stats Expand all files Collapse all files

1
.gitignore vendored
View File

@ -15,6 +15,7 @@ root_data
*.out *.out
*.txt *.txt
*.csv *.csv
*.dat
Cleopatra/ExtractXSec Cleopatra/ExtractXSec
Cleopatra/ExtractXSecFromText Cleopatra/ExtractXSecFromText

103
Armory/ClassCorrParas.h
View File

@ -3,20 +3,71 @@
//************************************** Correction parameters; //************************************** Correction parameters;
class CorrParas { class CorrParas {
private:
int defaultSize;
public: public:
CorrParas(){ CorrParas(int defaultSize = 100){
xnCorr.clear(); is_xn_OK = false;
xScale.clear(); is_xfxne_OK = false;
xfxneCorr.clear(); is_e_OK = false;
eCorr.clear(); is_xScale_OK = false;
rdtCorr.clear(); is_rdt_OK = false;
this->defaultSize = defaultSize;
}; };
void LoadAllCorrections(){
LoadXNCorr();
LoadXFXN2ECorr();
LoadECorr();
LoadXScaleCorr();
LoadRDTCorr();
}
void CheckCorrParasSize(size_t arraySize, size_t rdtSize){
printf("------------ Check Correction parameter sizes\n");
if( is_xn_OK && xnCorr.size() < arraySize ) {
printf(" xnCorr [%zu] < array size %zu. Set xnCorr[1..99] = 1.0 \n", xnCorr.size(), arraySize);
for( int i = 0; i < defaultSize; i++ ) xnCorr.push_back(1.0);
is_xn_OK = false;
}
if( is_xfxne_OK && xfxneCorr.size() < arraySize ) {
printf(" xfxneCorr [%zu] < array size %zu. Set xfxneCorr[1..99] = (0.0, 1.0) \n", xfxneCorr.size(), arraySize);
for( int i = 0; i < defaultSize; i++ ) xfxneCorr.push_back({0.0, 1.0});
is_xScale_OK = false;
}
if( is_e_OK && eCorr.size() < arraySize ) {
printf(" eCorr [%zu] < array size %zu. Set eCorr[1..99] = (1.0, 0.0) \n", xnCorr.size(), arraySize);
for( int i = 0; i < defaultSize; i++ ) eCorr.push_back({1.0, 0.0});
is_e_OK = false;
}
if( is_xScale_OK && xScale.size() < arraySize ) {
printf(" xScale [%zu] < array size %zu. Set xScale[1..99] = 1.0 \n", xScale.size(), arraySize);
for( int i = 0; i < defaultSize; i++ ) xScale.push_back(1.0);
is_xScale_OK = false;
}
if( is_rdt_OK && rdtCorr.size() < rdtSize ) {
printf(" rdtCorr [%zu] < array size %zu. Set rdtScale[1..99] = (0.0, 1.0) \n", rdtCorr.size(), arraySize);
for( int i = 0; i < defaultSize; i++ ) rdtCorr.push_back({0.0, 1.0});
is_rdt_OK = false;
}
printf("------------ Done Check Corr. Para. Size.\n");
}
bool is_xn_OK;
bool is_xfxne_OK;
bool is_e_OK;
bool is_xScale_OK;
bool is_rdt_OK;
std::vector<float> xnCorr; //correction of xn to match xf std::vector<float> xnCorr; //correction of xn to match xf
std::vector<float> xScale; // correction of x to be (0,1)
std::vector<std::vector<float>> xfxneCorr; //correction of xn and xf to match e std::vector<std::vector<float>> xfxneCorr; //correction of xn and xf to match e
std::vector<std::vector<float>> eCorr; // correction to e, ch -> MeV std::vector<std::vector<float>> eCorr; // correction to e, ch -> MeV
std::vector<float> xScale; // correction of x to be (0,1)
std::vector<std::vector<float>> rdtCorr; // correction of rdt, ch -> MeV std::vector<std::vector<float>> rdtCorr; // correction of rdt, ch -> MeV
//~========================================= xf = xn correction //~========================================= xf = xn correction
@ -28,15 +79,17 @@ public:
if( file.is_open() ){ if( file.is_open() ){
float a; float a;
while( file >> a ) xnCorr.push_back(a); while( file >> a ) xnCorr.push_back(a);
printf(".......... done.\n"); printf(".......... done. size:%zu\n", xnCorr.size());
is_xn_OK = true;
}else{ }else{
printf(".......... fail.\n"); for( int i = 0; i < defaultSize; i++ ) xnCorr.push_back(1.0);
printf(".......... fail. xnCorr[1..99] = 1.0\n");
is_xn_OK = false;
} }
file.close(); file.close();
if( verbose ) for(int i = 0; i < (int) xnCorr.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]); if( verbose ) for(int i = 0; i < (int) xnCorr.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
} }
//~========================================= X-Scale correction //~========================================= X-Scale correction
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.");
@ -44,11 +97,14 @@ public:
std::ifstream file; std::ifstream file;
file.open(fileName); file.open(fileName);
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a;
while( file >> a ) xScale.push_back(a); while( file >> a ) xScale.push_back(a);
printf("........ done.\n"); printf("........ done. size:%zu\n", xScale.size());
is_xScale_OK = true;
}else{ }else{
printf("........ fail.\n"); for( int i = 0; i < defaultSize; i++ ) xScale.push_back(1.0);
printf("........ fail. xScale[1..99] = 1.0\n");
is_xScale_OK = false;
} }
file.close(); file.close();
if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]); if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
@ -63,9 +119,12 @@ public:
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
while( file >> a >> b) xfxneCorr.push_back({a, b}); while( file >> a >> b) xfxneCorr.push_back({a, b});
printf("........ done.\n"); printf("........ done. size:%zu\n", xfxneCorr.size());
is_xfxne_OK = true;
}else{ }else{
printf("........ fail.\n"); for( int i = 0; i < defaultSize; i++ ) xfxneCorr.push_back({0.0, 1.0});
printf("........ fail. xfxneCorr[1..99] = (0.0, 1.0)\n");
is_xfxne_OK = false;
} }
file.close(); file.close();
if( verbose ) for(int i = 0; i < (int) xfxneCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, xfxneCorr[i][0], xfxneCorr[i][1]); if( verbose ) for(int i = 0; i < (int) xfxneCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, xfxneCorr[i][0], xfxneCorr[i][1]);
@ -80,9 +139,12 @@ public:
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
while( file >> a >> b) eCorr.push_back( {a, b} ); // 1/a1, a0 , e' = e * a1 + a0 while( file >> a >> b) eCorr.push_back( {a, b} ); // 1/a1, a0 , e' = e * a1 + a0
printf(".............. done.\n"); printf(".............. done. size:%zu\n", eCorr.size());
is_e_OK = true;
}else{ }else{
printf(".............. fail.\n"); for( int i = 0; i < defaultSize; i++ ) eCorr.push_back( {1.0, 0.0} );
printf(".............. fail. eCorr[1..99] = (1.0, 0.0) \n");
is_e_OK = false;
} }
file.close(); file.close();
if( verbose ) for(int i = 0; i < (int) eCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, eCorr[i][0], eCorr[i][1]); if( verbose ) for(int i = 0; i < (int) eCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, eCorr[i][0], eCorr[i][1]);
@ -97,9 +159,12 @@ public:
if( file.is_open() ){ if( file.is_open() ){
float a, b; float a, b;
while( file >> a >> b) rdtCorr.push_back({a, b}); while( file >> a >> b) rdtCorr.push_back({a, b});
printf("............ done.\n"); printf("............ done. size:%zu\n", rdtCorr.size());
is_rdt_OK = true;
}else{ }else{
printf("............ fail.\n"); for( int i = 0; i < defaultSize; i++ ) rdtCorr.push_back( {0.0, 1.0} );
printf("............ fail. rdtCorr[1..99] = (0.0, 1.0)\n");
is_rdt_OK = false;
} }
file.close(); file.close();
if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]); if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]);

7
Armory/ClassDetGeo.h
View File

@ -124,7 +124,7 @@ public:
std::vector<Array> array; std::vector<Array> array;
std::vector<Auxillary> aux; std::vector<Auxillary> aux;
void Print( bool printArray = false) ; void Print( int printArray = 0) ; // 0 = no print, -1 = print all, 1 = print only enabled
short GetArrayID(int id){ short GetArrayID(int id){
int detCount = 0; int detCount = 0;
@ -244,7 +244,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
} }
inline void DetGeo::Print(bool printArray){ inline void DetGeo::Print(int printArray){
printf("#####################################################\n"); printf("#####################################################\n");
printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan"); printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan");
@ -253,8 +253,9 @@ inline void DetGeo::Print(bool printArray){
printf(" z-Min : %8.2f mm\n", zMin); printf(" z-Min : %8.2f mm\n", zMin);
printf(" z-Max : %8.2f mm\n", zMax); printf(" z-Max : %8.2f mm\n", zMax);
if( printArray ) { if( printArray != 0 ) {
for( size_t i = 0; i < array.size() ; i++){ for( size_t i = 0; i < array.size() ; i++){
if( printArray > 0 && !array[i].enable ) continue;
printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled"); printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled");
array[i].Print(); array[i].Print();
aux[i].Print(); aux[i].Print();

4
Armory/GeneralSort.h
View File

@ -182,7 +182,7 @@ void GeneralSort::SetUpTree(){
newSaveTree->SetDirectory(saveFile); newSaveTree->SetDirectory(saveFile);
newSaveTree->AutoSave(); newSaveTree->AutoSave();
newSaveTree->Branch( "evID", &evID, "EventID/l"); // simply copy newSaveTree->Branch( "evID", &evID, "evID/l"); // simply copy
eE = new Float_t * [mapping::nDetType]; eE = new Float_t * [mapping::nDetType];
eT = new ULong64_t * [mapping::nDetType]; eT = new ULong64_t * [mapping::nDetType];
@ -198,7 +198,7 @@ void GeneralSort::SetUpTree(){
} }
newSaveTree->Branch( mapping::detTypeName[i].c_str(), eE[i], Form("%s[%d]/F", mapping::detTypeName[i].c_str(), mapping::detNum[i])); newSaveTree->Branch( mapping::detTypeName[i].c_str(), eE[i], Form("%s[%d]/F", mapping::detTypeName[i].c_str(), mapping::detNum[i]));
newSaveTree->Branch( (mapping::detTypeName[i]+"_t").c_str(), eT[i], Form("%s_Timestamp[%d]/l", mapping::detTypeName[i].c_str(), mapping::detNum[i])); newSaveTree->Branch( (mapping::detTypeName[i]+"_t").c_str(), eT[i], Form("%s_t[%d]/l", mapping::detTypeName[i].c_str(), mapping::detNum[i]));
} }

2
Armory/Process_Download
View File

@ -88,9 +88,11 @@ else
echo -e "$CYAN############### Only in SOLARIS MAC can donwload data. skip.$NC" echo -e "$CYAN############### Only in SOLARIS MAC can donwload data. skip.$NC"
fi fi
if [ -z ${rawDataPath}/RunTimeStamp.dat ]; then
echo -e "$YELLOW=============================================$NC" echo -e "$YELLOW=============================================$NC"
tail -10 ${rawDataPath}/RunTimeStamp.dat tail -10 ${rawDataPath}/RunTimeStamp.dat
echo -e "$YELLOW=============================================$NC" echo -e "$YELLOW=============================================$NC"
fi
count=`ls -1 ${rawDataPath}/${expName}_${RUN}_*.sol 2>/dev/null | wc -l` count=`ls -1 ${rawDataPath}/${expName}_${RUN}_*.sol 2>/dev/null | wc -l`
echo -e "========================= Number of Files : ${count}${YELLOW}" echo -e "========================= Number of Files : ${count}${YELLOW}"

248
working/ClassMonPlotter.h
View File

@ -9,6 +9,8 @@
#include "TH1.h" #include "TH1.h"
#include "TH2.h" #include "TH2.h"
#include "TCanvas.h" #include "TCanvas.h"
#include "TLine.h"
#include "TStyle.h"
/****************************************************************** /******************************************************************
* This is Plotter for Monitor.C. It contains * This is Plotter for Monitor.C. It contains
@ -53,10 +55,23 @@ public:
~MonPlotter(); ~MonPlotter();
void SetUpCanvas(TString title, int padSize, int divX, int divY); void SetUpCanvas(TString title, int padSize, int divX, int divY);
void SetUpHistograms(int * rawEnergyRange,
int * energyRange,
double * exRange,
int * thetaCMRange,
int * rdtDERange,
int * rdtERange,
int * coinTimeRange);
void LoadRDTGate(TString rdtCutFile);
void SetUpHistograms(int * rawEnergyRange, int * energyRange, double * exRange, int * thetaCMRange, int * rdtDERange, int * rdtERange);
void Plot(); void Plot();
void PlotRaw(bool isLog = false);
void PlotCal();
void PlotEZ();
void PlotEx();
TCanvas * canvas; TCanvas * canvas;
//====================== Histograms //====================== Histograms
@ -73,6 +88,7 @@ public:
TH2F ** hxfCal_xnCal; TH2F ** hxfCal_xnCal;
TH2F ** he_xsCal; // raw e vs xf TH2F ** he_xsCal; // raw e vs xf
TH2F ** he_x; // raw e vs x TH2F ** he_x; // raw e vs x
TH2F * heCal_ID;
//===== eCal V z //===== eCal V z
TH2F * heCal_z; TH2F * heCal_z;
@ -82,8 +98,14 @@ public:
TH2F * hrdt_ID; TH2F * hrdt_ID;
TH1F ** hrdt; // single recoil TH1F ** hrdt; // single recoil
TH1I * hrdtMulti;
TH2F ** hrdt2D; TH2F ** hrdt2D;
TH2F ** hrdt2Dg; TH2F ** hrdt2Dg; // gated
//====== tDiff
TH1F * htDiff;
TH1F * htDiffg;
//====== Ex data //====== Ex data
TH1F * hEx; TH1F * hEx;
@ -96,6 +118,9 @@ public:
TH2F * hEx_ThetaCM; TH2F * hEx_ThetaCM;
//======================= //=======================
//======= Recoil Cut
TObjArray * cutList;
private: private:
unsigned short aID; unsigned short aID;
@ -105,6 +130,7 @@ private:
float recoilOutter; float recoilOutter;
double zRange[2] ; // zMin, zMax double zRange[2] ; // zMin, zMax
TString canvasTitle;
TString suffix; TString suffix;
int numPad; int numPad;
@ -129,6 +155,7 @@ MonPlotter::MonPlotter(unsigned short arrayID, DetGeo * detGeo, int numRDT){
zRange[1] = detGeo->array[aID].zMax + 50; zRange[1] = detGeo->array[aID].zMax + 50;
canvas = nullptr; canvas = nullptr;
cutList = nullptr;
} }
@ -143,13 +170,19 @@ MonPlotter::~MonPlotter(){
delete hxn_ID; delete hxn_ID;
delete hArrayMulti; delete hArrayMulti;
delete heCal_ID;
delete heCal_zGC;
delete heCal_z; delete heCal_z;
delete hEx_ThetaCM; delete hEx_ThetaCM;
delete hExCut1; delete hExCut1;
delete hExCut2; delete hExCut2;
delete heCal_zGC;
delete hrdt_ID; delete hrdt_ID;
delete hrdtMulti;
delete htDiff;
delete htDiffg;
for( int i = 0; i < numDet ; i++ ){ for( int i = 0; i < numDet ; i++ ){
delete he[i]; delete he[i];
@ -185,13 +218,17 @@ MonPlotter::~MonPlotter(){
delete [] hrdt2D; delete [] hrdt2D;
delete [] hrdt2Dg; delete [] hrdt2Dg;
delete cutG;
delete cutList;
} }
void MonPlotter::SetUpCanvas(TString title, int padSize, int divX, int divY){ void MonPlotter::SetUpCanvas(TString title, int padSize, int divX, int divY){
canvas = new TCanvas("canavs" + suffix, title, 200 * aID, 200 * aID, divX * padSize, divY * padSize); canvas = new TCanvas("canavs" + suffix, title, 500 * aID, 0, divX * padSize, divY * padSize);
canvas->Divide(divX, divY); canvas->Divide(divX, divY);
numPad = divX * divY; numPad = divX * divY;
canvasTitle = title;
} }
@ -222,7 +259,13 @@ template<typename T> void MonPlotter::CreateListOfHist2D(T ** &histList,
} }
} }
void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double * exRange, int * thetaCMRange, int * rdtDERange, int * rdtERange){ void MonPlotter::SetUpHistograms(int * rawEnergyRange,
int * energyRange,
double * exRange,
int * thetaCMRange,
int * rdtDERange,
int * rdtERange,
int * coinTimeRange){
he_ID = new TH2F("he_ID" + suffix, "Raw e vs array ID; Array ID; Raw e", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]); he_ID = new TH2F("he_ID" + suffix, "Raw e vs array ID; Array ID; Raw e", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]);
hxf_ID = new TH2F("hxf_ID" + suffix, "Raw xf vs array ID; Array ID; Raw xf", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]); hxf_ID = new TH2F("hxf_ID" + suffix, "Raw xf vs array ID; Array ID; Raw xf", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]);
@ -237,14 +280,14 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double
CreateListOfHist2D(hxf_xn, numDet, "hxf_xn", "Raw xf vs. xn (ch=%d);xf (channel);xn (channel)" , 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); CreateListOfHist2D(hxf_xn, numDet, "hxf_xn", "Raw xf vs. xn (ch=%d);xf (channel);xn (channel)" , 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
CreateListOfHist2D(he_xs, numDet, "he_xs", "Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); CreateListOfHist2D(he_xs, numDet, "he_xs", "Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]);
CreateListOfHist1D(heCal, numDet, "heCal", "Corrected e (ch=%d); e (MeV); count", 2000, energyRange[0], energyRange[1]); CreateListOfHist2D(he_x , numDet, "he_x", "Raw e vs x (ch=%d); x (mm); Raw e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, -1, detLength +1);
CreateListOfHist2D(hxfCal_xnCal, numDet, "hxfCal_xnCal", "Corrected XF vs. XN (ch=%d);XF (channel);XN (channel)", 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]); CreateListOfHist2D(hxfCal_xnCal, numDet, "hxfCal_xnCal", "Corrected XF vs. XN (ch=%d);XF (channel);XN (channel)", 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]);
CreateListOfHist2D(he_xsCal , numDet, "he_xsCal", "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(he_xsCal , numDet, "he_xsCal", "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(he_x , numDet, "he_x", "Raw e vs x (ch=%d); x (mm); Raw e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, -1, detLength +1); CreateListOfHist1D(heCal, numDet, "heCal", "Corrected e (ch=%d); e (MeV); count", 2000, energyRange[0], energyRange[1]);
//====================== E-Z plot //====================== E-Z plot
heCal_ID = new TH2F("heCal_ID" + suffix , "E vs. ID; ID;E (MeV)" , numDet, 0, numDet, 400, energyRange[0], energyRange[1]);
heCal_z = new TH2F("heCal_z" + suffix , "E vs. Z;Z (mm);E (MeV)" , 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]); heCal_z = new TH2F("heCal_z" + suffix , "E vs. Z;Z (mm);E (MeV)" , 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]);
heCal_zGC = new TH2F("heCal_zGC" + suffix ,"E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]); heCal_zGC = new TH2F("heCal_zGC" + suffix ,"E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]);
@ -255,6 +298,8 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double
hrdt_ID = new TH2F("hrdt_ID" + suffix, "Raw RDT vs ID; ID; Raw RDT", numRDT, 0, numRDT, 400, rdtRange[0], rdtRange[1]); hrdt_ID = new TH2F("hrdt_ID" + suffix, "Raw RDT vs ID; ID; Raw RDT", numRDT, 0, numRDT, 400, rdtRange[0], rdtRange[1]);
hrdtMulti = new TH1I("hrdtMulti" + suffix, "RDT Multiplicity", numRDT, 0, numRDT);
hrdt = new TH1F * [numRDT]; hrdt = new TH1F * [numRDT];
hrdt2D = new TH2F * [numRDT/2]; hrdt2D = new TH2F * [numRDT/2];
hrdt2Dg = new TH2F * [numRDT/2]; hrdt2Dg = new TH2F * [numRDT/2];
@ -271,6 +316,11 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double
} }
} }
//===================== tDiff = array_t - rdt_t
htDiff = new TH1F("htDiff" + suffix, "tDiff = e_t - rdt_t", (coinTimeRange[1]-coinTimeRange[0]), coinTimeRange[0], coinTimeRange[1]);
htDiffg = new TH1F("htDiffg" + suffix, "tDiff = e_t - rdt_t (gated)", (coinTimeRange[1]-coinTimeRange[0]), coinTimeRange[0], coinTimeRange[1]);
htDiffg->SetLineColor(2);
//===================== energy spectrum //===================== energy spectrum
hEx = new TH1F("hEx" + suffix, 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]); hEx = new TH1F("hEx" + suffix, 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]);
@ -297,9 +347,187 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double
} }
void MonPlotter::Plot(){ void MonPlotter::Plot(){
for( int i = 1; i <= numPad; i++ ){
canvas->cd(i);
switch (i){
case 1: heCal_z->Draw("colz");break;
case 2: heCal_zGC->Draw("colz");break;
case 3: htDiff->Draw("");break;
case 4: hEx->Draw("colz");break;
default:break;
}
}
}
for( int i = 0; i < numPad; i++ ){ void MonPlotter::LoadRDTGate(TString rdtCutFile){
canvas->cd(i+1);
if( rdtCutFile == "" ) return ;
TFile * fCut = new TFile(rdtCutFile);
bool isCutFileOpen = fCut->IsOpen();
if(!isCutFileOpen) {
printf( "Failed to open rdt-cutfile 1 : %s\n" , fileName.Data());
}else{
cutList = (TObjArray *) fCut->FindObjectAny("cutList");
if( cutList ){
int numCut = cutList->GetEntries();
printf("=========== found %d cutG in %s \n", numCut, fCut->GetName());
for(int i = 0; i < numCut ; i++){
printf("cut name : %s , VarX: %s, VarY: %s, numPoints: %d \n",
cutList->At(i)->GetName(),
((TCutG*)cutList->At(i))->GetVarX(),
((TCutG*)cutList->At(i))->GetVarY(),
((TCutG*)cutList->At(i))->GetN()
);
}
}
}
}
//^#######################################################
void MonPlotter::PlotRaw(bool isLog){
TCanvas * cRawID = new TCanvas("cRawID", Form("Raw e, Ring, xf, xn vs ID | %s", canvasTitle.Data()), 100 + 500 * aID, 100, 1200, 800);
cRawID->Clear(); cRawID->Divide(2,2);
cRawID->cd(1); he_ID->Draw("colz");
cRawID->cd(2); hArrayMulti->Draw();
cRawID->cd(3); hxf_ID->Draw("colz");
cRawID->cd(4); hxn_ID->Draw("colz");
int padSize = 200;
int canvasSize[2] = {padSize * colDet, padSize * rowDet};
TCanvas * cRawE = new TCanvas("cRawE" + suffix,Form("E raw | %s", canvasTitle.Data()), 200 + 500 * aID, 200, canvasSize[0], canvasSize[1]);
cRawE->Clear(); cRawE->Divide(colDet,rowDet);
for (Int_t i=0; i < numDet; i++) {
cRawE->cd(i+1);
cRawE->cd(i+1)->SetGrid();
if( isLog ) cRawE->cd(i+1)->SetLogy();
he[i]->Draw("");
}
TCanvas *cRawXf = new TCanvas("cRawXf" + suffix,Form("Xf raw | %s", canvasTitle.Data()), 300 + 500 * aID, 300, canvasSize[0], canvasSize[1]);
cRawXf->Clear(); cRawXf->Divide(colDet,rowDet);
for (Int_t i=0; i<numDet; i++) {
cRawXf->cd(i+1);
cRawXf->cd(i+1)->SetGrid();
if( isLog ) cRawXf->cd(i+1)->SetLogy();
hxf[i]->Draw("");
}
TCanvas *cRawXn = new TCanvas("cRawXn" + suffix,Form("Xn raw | %s", canvasTitle.Data()), 400 + 500 * aID, 400, canvasSize[0], canvasSize[1]);
cRawXn->Clear();cRawXn->Divide(colDet,rowDet);
for (Int_t i=0; i<numDet; i++) {
cRawXn->cd(i+1);
cRawXn->cd(i+1)->SetGrid();
if( isLog ) cRawXn->cd(i+1)->SetLogy();
hxn[i]->Draw("");
}
TCanvas *cxfxn = new TCanvas("cxfxn" + suffix,Form("XF vs. XN | %s", canvasTitle.Data()), 500 + 500 * aID, 500, canvasSize[0], canvasSize[1]);
cxfxn->Clear(); cxfxn->Divide(colDet,rowDet);
for (Int_t i=0;i<numDet;i++) {
cxfxn->cd(i+1);
cxfxn->cd(i+1)->SetGrid();
hxf_xn[i]->Draw("col");
}
TCanvas *cxfxne = new TCanvas("cxfxne" + suffix,Form("E - XF+XN | %s", canvasTitle.Data()), 600 + 500 * aID, 600, canvasSize[0], canvasSize[1]);
cxfxne->Clear(); cxfxne->Divide(colDet,rowDet);
TLine line(0,0, 4000, 4000); line.SetLineColor(2);
for (Int_t i=0;i<numDet;i++) {
cxfxne->cd(i+1);
cxfxne->cd(i+1)->SetGrid();
he_xs[i]->Draw("col");
line.Draw("same");
}
}
void MonPlotter::PlotCal(){
int padSize = 200;
int canvasSize[2] = {padSize * colDet, padSize * rowDet};
TCanvas *ceVx = new TCanvas("ceVx" + suffix, Form("E vs. X = (xf-xn)/e | %s", canvasTitle.Data()), 100 + 500 * aID, 100, canvasSize[0], canvasSize[1]);
ceVx->Clear(); ceVx->Divide(colDet,rowDet);
for (Int_t i=0;i<numDet;i++) {
ceVx->cd(i+1); he_x[i]->Draw("col");
}
TCanvas *cxfxneC = new TCanvas("cxfxneC" + suffix,Form("Raw E - Corrected XF+XN | %s", canvasTitle.Data()), 200 + 500 * aID, 200, canvasSize[0], canvasSize[1]);
cxfxneC->Clear(); cxfxneC->Divide(colDet,rowDet);
TLine line(0,0, 4000, 4000); line.SetLineColor(2);
for (Int_t i=0;i<numDet;i++) {
cxfxneC->cd(i+1);
cxfxneC->cd(i+1)->SetGrid();
he_xsCal[i]->Draw("col");
line.Draw("same");
}
TCanvas *cEC = new TCanvas("cEC" + suffix,Form("E corrected | %s", canvasTitle.Data()), 300 + 500 * aID, 300, canvasSize[0], canvasSize[1]);
cEC->Clear();cEC->Divide(colDet,rowDet);
for (Int_t i=0; i<numDet; i++) {
cEC->cd(i+1);
cEC->cd(i+1)->SetGrid();
heCal[i]->Draw("");
}
TCanvas *cEC2 = new TCanvas("cEC2" + suffix,Form("E corrected | %s", canvasTitle.Data()), 400 + 500 * aID, 400, canvasSize[0], canvasSize[1]);
cEC2->Clear();
heCal_ID->Draw("colz");
TCanvas *cxfxnC = new TCanvas("cxfxnC" + suffix,Form("XF vs XN corrected | %s", canvasTitle.Data()), 500 + 500 * aID, 500, canvasSize[0], canvasSize[1]);
cxfxnC->Clear(); cxfxnC->Divide(colDet,rowDet);
for (Int_t i=0;i<numDet;i++) {
cxfxnC->cd(i+1);
cxfxnC->cd(i+1)->SetGrid();
hxfCal_xnCal[i]->Draw("col");
}
}
void MonPlotter::PlotEZ(){
TCanvas *cecalVz = new TCanvas("cevalVz" + suffix,Form("ECALVZ : %s", canvasTitle.Data()),1000, 650);
cecalVz->Clear(); cecalVz->Divide(2,1);
gStyle->SetOptStat("neiou");
cecalVz->cd(1);heCal_z->Draw("col");
cecalVz->cd(2);heCal_zGC->Draw("col");
}
void MonPlotter::PlotEx(){
TCanvas *cex = new TCanvas("cex" + suffix,Form("EX : %s", canvasTitle.Data()),0, 0, 1000,650);
cex->Clear();
gStyle->SetOptStat("neiou");
hEx->Draw("");
TCanvas *cexI = new TCanvas("cexI" + suffix,Form("EX : %s", canvasTitle.Data()),500, 0, 1600,1000);
cexI->Clear();cexI->Divide(colDet,rowDet);
gStyle->SetOptStat("neiou");
for( int i = 0; i < numDet; i++){
cexI->cd(i+1);
hExi[i]->Draw("");
}
TCanvas *cExThetaCM = new TCanvas("cExThetaCM" + suffix,Form("EX - ThetaCM | %s", canvasTitle.Data()), 500, 500, 650,650);
cExThetaCM->Clear();
gStyle->SetOptStat("neiou");
hEx_ThetaCM->Draw("colz");
TCanvas *cExVxCal = new TCanvas("cExVxCal" + suffix,Form("EX | %s", canvasTitle.Data()),200, 200, 1600,1000);
cExVxCal->Clear();
gStyle->SetOptStat("neiou");
cExVxCal->Divide(colDet,rowDet);
for( int i = 0; i < numDet; i++){
cExVxCal->cd(i+1);
hEx_xCal[i]->SetMarkerStyle(7);
hEx_xCal[i]->Draw();
} }
} }

497
working/MonAnalyzer.C
View File

@ -1,88 +1,443 @@
#define MonAnalyzer_cxx
// The class definition in MonAnalyzer.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
#include "../Armory/AnalysisLib.h"
#include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h"
#include "../Armory/ClassCorrParas.h"
// #include "../Cleopatra/ClassHelios.h"
#include "../Cleopatra/ClassTransfer.h"
// The following methods are defined in this file: #include "ClassMonPlotter.h"
// Begin(): called every time a loop on the tree starts, #include "Mapping.h"
// a convenient place to create your histograms.
// SlaveBegin(): called after Begin(), when on PROOF called only on the
// slave servers.
// Process(): called for each event, in this function you decide what
// to read and fill your histograms.
// SlaveTerminate: called at the end of the loop on the tree, when on PROOF
// called only on the slave servers.
// Terminate(): called at the end of the loop on the tree,
// a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// root> T->Process("MonAnalyzer.C")
// root> T->Process("MonAnalyzer.C","some options")
// root> T->Process("MonAnalyzer.C+")
//
#include "TFile.h"
#include "TChain.h"
#include "TH1F.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TTreeReaderArray.h"
#include "TClonesArray.h"
#include "TGraph.h"
#include "TH2.h"
#include "TStyle.h"
#include "TStopwatch.h"
#include "TMath.h"
#include "MonAnalyzer.h" #include "vector"
#include <TH2.h> //^############################################ User setting
#include <TStyle.h> int rawEnergyRange[2] = { 0, 3000}; /// share with e, xf, xn
int energyRange[2] = { 0, 10}; /// in the E-Z plot
int rdtDERange[2] = { 0, 80};
int rdtERange[2] = { 0, 80};
int thetaCMRange[2] = {0, 80};
void MonAnalyzer::Begin(TTree * /*tree*/) double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV]
{
// The Begin() function is called at the start of the query.
// When running with PROOF Begin() is only called on the client.
// The tree argument is deprecated (on PROOF 0 is passed).
TString option = GetOption(); int coinTimeRange[2] = { -200, 200};
//---Gate
bool isTimeGateOn = true;
int timeGate[2] = {-20, 12}; /// min, max, 1 ch = 10 ns
double eCalCut[2] = {0.5, 20}; /// lower & higher limit for eCal
double xGate = 0.9; ///cut out the edge
double thetaCMGate = 10; /// deg
std::vector<int> skipDetID = {11} ;
std::vector<TString> rdtCutFile1 = {"", ""}; /// {reaction-0, reaction-1}, can add more for more reactions
// TString rdtCutFile2 = "";
// TString ezCutFile = "";//"ezCut.root";
//^############################################ end of user setting
MonPlotter ** plotter = nullptr;
int numGeo = 1;
void MonAnalyzer(){
printf("#####################################################################\n");
printf("####################### MonAnalyzer.C #######################\n");
printf("#####################################################################\n");
TChain *chain = new TChain("gen_tree");
//chain->Add("../root_data/gen_run043.root");
chain->Add("../root_data/trace_run029.root");
TObjArray * fileList = chain->GetListOfFiles();
printf("\033[0;31m========================================== Number of Files : %2d\n",fileList->GetEntries());
fileList->Print();
printf("========================================== Number of Files : %2d\033[0m\n",fileList->GetEntries());
printf("///////////////////////////////////////////////////////////////////\n");
printf(" Total Number of entries : %llu \n", chain->GetEntries());
printf("///////////////////////////////////////////////////////////////////\n");
TTreeReader reader(chain);
TTreeReaderValue<ULong64_t> evID = {reader, "evID"};
TTreeReaderArray<Float_t> e = {reader, "e"};
TTreeReaderArray<ULong64_t> e_t = {reader, "e_t"};
TTreeReaderArray<Float_t> xf = {reader, "xf"};
TTreeReaderArray<Float_t> xn = {reader, "xn"};
TTreeReaderArray<Float_t> rdt = {reader, "rdt"};
TTreeReaderArray<ULong64_t> rdt_t = {reader, "rdt_t"};
//TODO
// TTreeReaderArray<TGraph> array = {reader, "trace"};
ULong64_t NumEntries = chain->GetEntries();
//*==========================================
DetGeo * detGeo = new DetGeo("detectorGeo.txt");
numGeo = detGeo->numGeo;
printf("================== num. of Arrays : %d\n", numGeo);
int numTotArray = 0;
detGeo->Print(1);
for( size_t i = 0; i < detGeo->array.size(); i++ ){
if( detGeo->array[i].enable ) numTotArray += detGeo->array[i].numDet;
} }
void MonAnalyzer::SlaveBegin(TTree * /*tree*/) //*==========================================
{ TransferReaction * transfer = new TransferReaction[numGeo];
// The SlaveBegin() function is called after the Begin() function. int tempCount = 0;
// When running with PROOF SlaveBegin() is called on each slave server. for( int i = 0; i < (int) detGeo->array.size() ; i++){
// The tree argument is deprecated (on PROOF 0 is passed). if( !detGeo->array[i].enable ) continue;
transfer[tempCount].SetReactionFromFile("reactionConfig.txt", i);
tempCount ++;
}
TString option = GetOption(); //*==========================================
CorrParas * corr = new CorrParas;
corr->LoadAllCorrections();
corr->CheckCorrParasSize(numTotArray, mapping::NRDT);
plotter = new MonPlotter *[numGeo];
for( int i = 0; i < numGeo; i++ ) {
plotter[i] = new MonPlotter(i, detGeo, mapping::NRDT);
plotter[i]->SetUpCanvas("haha", 500, 3, 2); //TODO canvaseTitle
plotter[i]->SetUpHistograms(rawEnergyRange, energyRange, exRange, thetaCMRange, rdtDERange, rdtERange, coinTimeRange);
}
//*========================================== Load RDT Cuts
tempCount = 0;
for( int i = 0; i < (int) detGeo->array.size() ; i++){
if( !detGeo->array[i].enable ) continue;
plotter[tempCount]->LoadRDTGate(rdtCutFile1[i]);
tempCount ++;
}
//TODO make the data class.
std::vector<double>eCal (numTotArray);
std::vector<double>xfCal (numTotArray);
std::vector<double>xnCal (numTotArray);
std::vector<double>x (numTotArray);
std::vector<double>xCal (numTotArray);
std::vector<double>z (numTotArray);
//^###########################################################
//^ * Process
//^###########################################################
ULong64_t processedEntries = 0;
float Frac = 0.1;
TStopwatch StpWatch;
StpWatch.Start();
while (reader.Next()) {
//*============================================= Array;
int arrayMulti[numGeo] ; //array multiplicity, when any is calculated.
for( int i = 0; i < numGeo; i++ ) arrayMulti[i] = 0;
bool rdtgate1 = false;
for( int id = 0; id < (int) e.GetSize() ; id++ ){
short aID = detGeo->GetArrayID(id);
if( aID < 0 ) continue;
//@================== Filling raw data
plotter[aID]->he_ID->Fill(id, e[id]);
plotter[aID]->hxf_ID->Fill(id, xf[id]);
plotter[aID]->hxn_ID->Fill(id, xn[id]);
plotter[aID]->he[id]->Fill(e[id]);
plotter[aID]->hxf[id]->Fill(xf[id]);
plotter[aID]->hxn[id]->Fill(xn[id]);
plotter[aID]->hxf_xn[id]->Fill(xf[id],xn[id]);
plotter[aID]->he_xs[id]->Fill(xf[id]+xn[id], e[id]);
//@==================== Basic gate
if( TMath::IsNaN(e[id]) ) continue ;
if( TMath::IsNaN(xn[id]) && TMath::IsNaN(xf[id]) ) continue ;
//@==================== Skip detector
bool skipFlag = false;
for( unsigned int kk = 0; kk < skipDetID.size() ; kk++){
if( id == skipDetID[kk] ) {
skipFlag = true;
break;
}
}
if (skipFlag ) continue;
//@==================== When e, xn, or xf is valid.
arrayMulti[aID] ++;
//@==================== Calibrations go here
if( corr->xnCorr.size() && corr->xfxneCorr.size() ) xnCal[id] = xn[id] * corr->xnCorr[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0];
if( corr->xfxneCorr.size() ) xfCal[id] = xf[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0];
if( corr->eCorr.size() ) eCal[id] = e[id] / corr->eCorr[id][0] + corr->eCorr[id][1];
if( eCal[id] < eCalCut[0] || eCalCut[1] < eCal[id] ) continue;
//@===================== fill Calibrated data
plotter[aID]->heCal[id]->Fill(eCal[id]);
plotter[aID]->hxfCal_xnCal[id]->Fill(xfCal[id], xnCal[id]);
plotter[aID]->he_xsCal[id]->Fill(xnCal[id] + xfCal[id], e[id]);
//@===================== calculate X
if( (xf[id] > 0 || !TMath::IsNaN(xf[id])) && ( xn[id] > 0 || !TMath::IsNaN(xn[id])) ) {
///x[id] = 0.5*((xf[id]-xn[id]) / (xf[id]+xn[id]))+0.5;
x[id] = 0.5*((xf[id]-xn[id]) / e[id])+0.5;
}
/// range of x is (0, 1)
if ( !TMath::IsNaN(xf[id]) && !TMath::IsNaN(xn[id]) ) xCal[id] = 0.5 + 0.5 * (xfCal[id] - xnCal[id] ) / e[id];
if ( !TMath::IsNaN(xf[id]) && TMath::IsNaN(xn[id]) ) xCal[id] = xfCal[id]/ e[id];
if ( TMath::IsNaN(xf[id]) && !TMath::IsNaN(xn[id]) ) xCal[id] = 1.0 - xnCal[id]/ e[id];
//@=================== Fill in histogram
plotter[aID]->he_x[id]->Fill(x[id],e[id]);
plotter[aID]->hxfCal_xnCal[id]->Fill(xfCal[id],xnCal[id]);
plotter[aID]->he_xsCal[id]->Fill(e[id],xnCal[id] + xfCal[id]);
//@======= Scale xcal from (0,1)
if( corr->xScale.size() ) xCal[id] = (xCal[id]-0.5)/corr->xScale[id] + 0.5; /// if include this scale, need to also inclused in Cali_littleTree
if( abs(xCal[id] - 0.5) > xGate/2. ) continue;
//@==================== calculate Z
if( aID >= 0 ){
int colIndex = id % detGeo->array[aID].colDet;
if( detGeo->array[aID].firstPos > 0 ) {
z[id] = detGeo->array[aID].detLength*(1.0-xCal[id]) + detGeo->array[aID].detPos[colIndex];
}else{
z[id] = detGeo->array[aID].detLength*(xCal[id]-1.0) + detGeo->array[aID].detPos[colIndex];
}
}
//@=================== Fill histogram
plotter[aID]->heCal[id]->Fill(eCal[id]);
plotter[aID]->heCal_ID->Fill(id, eCal[id]);
plotter[aID]->heCal_z->Fill(z[id],eCal[id]);
//@=================== Recoil Gate
if( plotter[aID]->cutList ){
for(int i = 0 ; i < cutList1->GetEntries() ; i++ ){
TCutG * cutG = (TCutG *)cutList1->At(i) ;
if(cutG->IsInside(rdt[2*i],rdt[2*i+1])) {
rdtgate1= true;
break; /// only one is enough
}
}
// for(int i = 0 ; i < cutList2->GetEntries() ; i++ ){
// cutG = (TCutG *)cutList2->At(i) ;
// if(cutG->IsInside(rdt[2*i],rdt[2*i+1])) {
// //if(cutG->IsInside(rdt[2*i]+ rdt[2*i+1],rdt[2*i+1])) {
// rdtgate2= true;
// break; /// only one is enough
// }
// }
}else{
rdtgate1 = true;
// rdtgate2 = true;
}
//@================ coincident with Recoil when z is calculated.
if( !TMath::IsNaN(z[id]) ) {
for( int j = 0; j < mapping::NRDT ; j++){
if( TMath::IsNaN(rdt[j]) ) continue;
int tdiff = rdt_t[j] - e_t[id];
if( j%2 == 1) {
plotter[aID]->htDiff->Fill(tdiff);
// if((rdtgate1 || rdtgate2) && (eCalCut[1] > eCal[id] && eCal[id]>eCalCut[0])) {
// plotter[aID]->htdiffg->Fill(tdiff);
// }
}
// hArrayRDTMatrix->Fill(id, j);
// if( isTimeGateOn && timeGate[0] < tdiff && tdiff < timeGate[1] ) {
// if(j % 2 == 0 ) hrdt2Dg[j/2]->Fill(rdt[j],rdt[j+1]); /// x=E, y=dE
// ///if(j % 2 == 0 ) hrdt2Dg[j/2]->Fill(rdt[j+1],rdt[j]); /// x=dE, y=E
// hArrayRDTMatrixG->Fill(id, j);
// ///if( rdtgate1) hArrayRDTMatrixG->Fill(id, j);
// hrdtg[j]->Fill(rdt[j]);
// coinFlag = true;
// }
}
}
// if( !isTimeGateOn ) coinFlag = true;
//@================ E-Z gate
// if( EZCut ) {
// if( EZCut->IsInside(z[id], eCal[id]) ) ezGate = true;
// }else{
// ezGate = true;
// }
// if( coinFlag && (rdtgate1 || rdtgate2) && ezGate){
// plotter[arrayID]->heCalVzGC->Fill( z[id] , eCal[id] );
// heCalVxCalG[id]->Fill(xCal[id]*detGeo->array[arrayID].detLength, eCal[id]);
// multiEZ ++;
// isGoodEventFlag = true;
// }
}//*====== end of array loop
for( int i = 0 ; i < numGeo ; i++ ) plotter[i]->hArrayMulti->Fill(arrayMulti[i]);
//*********** RECOILS ***********************************************/
//Fill both plotter
int recoilMulti = 0;
for( int j = 0; j < numGeo; j++ ){
for( int i = 0; i < mapping::NRDT ; i++){
plotter[j]->hrdt_ID->Fill(i, rdt[i]);
plotter[j]->hrdt[i]->Fill(rdt[i]);
recoilMulti++;
if( i % 2 == 0 ){
plotter[j]->hrdt2D[i/2]->Fill(rdt[i],rdt[i+1]); //E-dE
}
}
plotter[j]->hrdtMulti->Fill(recoilMulti);
}
//@*********** Ex and thetaCM ****************************************/
for(Int_t id = 0; id < numTotArray ; id++){
if( TMath::IsNaN(e[id]) ) continue ;
if( TMath::IsNaN(z[id]) ) continue ;
if( eCal[id] < eCalCut[0] ) continue ;
if( eCal[id] > eCalCut[1] ) continue ;
short aID = detGeo->GetArrayID(id);
if( aID < 0 ) continue;
std::pair<double, double> ExThetaCM = transfer[aID].CalExThetaCM(eCal[id], z[id], detGeo->Bfield, detGeo->array[aID].detPerpDist);
double Ex = ExThetaCM.first;
double thetaCM = ExThetaCM.second;
if( thetaCM > thetaCMGate ) {
plotter[aID]->hEx->Fill(Ex);
plotter[aID]->hExi[id]->Fill(Ex);
plotter[aID]->hEx_xCal[id]->Fill(xCal[id], Ex);
plotter[aID]->hEx_ThetaCM->Fill(thetaCM, Ex);
// if( rdtgate1 ) {
// plotter[arrayID]->hExCut1->Fill(Ex);
// plotter[arrayID]->hExThetaCM->Fill(thetaCM, Ex);
// }
// if( rdtgate2 ) {
// plotter[arrayID]->hExCut2->Fill(Ex);
// plotter[arrayID]->hExThetaCM->Fill(thetaCM, Ex);
// }
}
}
//*============================================ Progress Bar
processedEntries ++;
if (processedEntries >= NumEntries*Frac - 1 ) {
TString msg; msg.Form("%llu", NumEntries/1000);
int len = msg.Sizeof();
printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
Frac*100, len, processedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
StpWatch.Start(kFALSE);
Frac += 0.1;
}
if( processedEntries > 1000 ) break;
}//^############################################## End of Process
gStyle->SetOptStat("neiou");
gStyle->GetAttDate()->SetTextSize(0.02);
gStyle->SetOptDate(1);
gStyle->SetDateX(0);
gStyle->SetDateY(0);
//TODO, an easy method to config plotter::Plot()
for( int i = 0; i < detGeo->numGeo ; i++){
plotter[i]->Plot();
}
//^###############################################
printf("------------------- List of Plots -------------------\n");
// printf(" newCanvas() - Create a new Canvas\n");
// printf("-----------------------------------------------------\n");
printf(" raw() - Raw data\n");
printf(" cal() - Calibrated data\n");
printf("-----------------------------------------------------\n");
printf(" ez() - Energy vs. Z\n");
printf(" excite() - Excitation Energy\n");
// printf(" recoils() - Raw DE vs. E Recoil spectra\n");
//printf(" elum() - Luminosity Energy Spectra\n");
//printf(" ic() - Ionization Chamber Spectra\n");
// printf("-----------------------------------------------------\n");
// printf(" eCalVzRow() - Energy vs. Z for each row\n");
// printf(" ExThetaCM() - Ex vs ThetaCM\n");
// printf(" ExVxCal() - Ex vs X for all %d detectors\n", numDet);
// printf("-----------------------------------------------------\n");
// printf(" ShowFitMethod() - Shows various fitting methods \n");
// printf(" RDTCutCreator() - Create RDT Cuts [May need to edit]\n");
// printf(" Check_rdtGate() - Check RDT Cuts. \n");
// printf(" readTrace() - read trace from gen_runXXX.root \n");
// printf(" readRawTrace() - read trace from runXXX.root \n");
// printf(" Check1D() - Count Integral within a range\n");
// printf("-----------------------------------------------------\n");
// printf(" %s\n", canvasTitle.Data());
printf("-----------------------------------------------------\n");
} }
bool MonAnalyzer::Process(Long64_t entry) //%=============================================
{ void raw(bool isLog = false, int arrayID = -1){
// The Process() function is called for each entry in the tree (or possibly if( arrayID < 0 ){
// keyed object in the case of PROOF) to be processed. The entry argument for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotRaw(isLog);
// specifies which entry in the currently loaded tree is to be processed. }else{
// When processing keyed objects with PROOF, the object is already loaded if( arrayID < numGeo) plotter[arrayID]->PlotRaw(isLog);
// and is available via the fObject pointer. }
//
// This function should contain the \"body\" of the analysis. It can contain
// simple or elaborate selection criteria, run algorithms on the data
// of the event and typically fill histograms.
//
// The processing can be stopped by calling Abort().
//
// Use fStatus to set the return value of TTree::Process().
//
// The return value is currently not used.
fReader.SetLocalEntry(entry);
return true;
} }
void MonAnalyzer::SlaveTerminate() void cal(int arrayID = -1){
{ if( arrayID < 0 ){
// The SlaveTerminate() function is called after all entries or objects for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotCal();
// have been processed. When running with PROOF SlaveTerminate() is called }else{
// on each slave server. if( arrayID < numGeo) plotter[arrayID]->PlotCal();
}
} }
void MonAnalyzer::Terminate() void ez(int arrayID = -1){
{ if( arrayID < 0 ){
// The Terminate() function is the last function to be called during for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotEZ();
// a query. It always runs on the client, it can be used to present }else{
// the results graphically or save the results to file. if( arrayID < numGeo) plotter[arrayID]->PlotEZ();
}
}
void excited(int arrayID = -1){
if( arrayID < 0 ){
for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotEx();
}else{
if( arrayID < numGeo) plotter[arrayID]->PlotEx();
}
} }

115
working/MonAnalyzer.h
View File

@ -1,115 +0,0 @@
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Mon Jul 8 13:26:58 2024 by ROOT version 6.32.02
// from TTree gen_tree/Tree After GeneralSort
// found on file: ../root_data/gen_run043.root
//////////////////////////////////////////////////////////
#ifndef MonAnalyzer_h
#define MonAnalyzer_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>
// Headers needed by this particular selector
#include "TClonesArray.h"
class MonAnalyzer : public TSelector {
public :
TTreeReader fReader; //!the tree reader
TTree *fChain = 0; //!pointer to the analyzed TTree or TChain
// Readers to access the data (delete the ones you do not need).
TTreeReaderValue<ULong64_t> EventID = {fReader, "evID"};
TTreeReaderArray<Float_t> e = {fReader, "e"};
TTreeReaderArray<ULong64_t> e_Timestamp = {fReader, "e_t"};
TTreeReaderArray<Float_t> xf = {fReader, "xf"};
TTreeReaderArray<ULong64_t> xf_Timestamp = {fReader, "xf_t"};
TTreeReaderArray<Float_t> xn = {fReader, "xn"};
TTreeReaderArray<ULong64_t> xn_Timestamp = {fReader, "xn_t"};
TTreeReaderArray<Float_t> rdt = {fReader, "rdt"};
TTreeReaderArray<ULong64_t> rdt_Timestamp = {fReader, "rdt_t"};
TTreeReaderArray<unsigned int> trace_fUniqueID = {fReader, "trace.fUniqueID"};
TTreeReaderArray<unsigned int> trace_fBits = {fReader, "trace.fBits"};
TTreeReaderArray<TString> trace_fName = {fReader, "trace.fName"};
TTreeReaderArray<TString> trace_fTitle = {fReader, "trace.fTitle"};
TTreeReaderArray<Short_t> trace_fLineColor = {fReader, "trace.fLineColor"};
TTreeReaderArray<Short_t> trace_fLineStyle = {fReader, "trace.fLineStyle"};
TTreeReaderArray<Short_t> trace_fLineWidth = {fReader, "trace.fLineWidth"};
TTreeReaderArray<Short_t> trace_fFillColor = {fReader, "trace.fFillColor"};
TTreeReaderArray<Short_t> trace_fFillStyle = {fReader, "trace.fFillStyle"};
TTreeReaderArray<Short_t> trace_fMarkerColor = {fReader, "trace.fMarkerColor"};
TTreeReaderArray<Short_t> trace_fMarkerStyle = {fReader, "trace.fMarkerStyle"};
TTreeReaderArray<Float_t> trace_fMarkerSize = {fReader, "trace.fMarkerSize"};
TTreeReaderArray<Double_t> trace_fMinimum = {fReader, "trace.fMinimum"};
TTreeReaderArray<Double_t> trace_fMaximum = {fReader, "trace.fMaximum"};
TTreeReaderArray<TString> trace_fOption = {fReader, "trace.fOption"};
TTreeReaderArray<Float_t> trace_e = {fReader, "we"};
TTreeReaderArray<Float_t> trace_e_time = {fReader, "weT"};
TTreeReaderArray<Float_t> trace_e_rise = {fReader, "weR"};
TTreeReaderArray<Float_t> trace_xf = {fReader, "wxf"};
TTreeReaderArray<Float_t> trace_xf_time = {fReader, "wxfT"};
TTreeReaderArray<Float_t> trace_xf_rise = {fReader, "wxfR"};
TTreeReaderArray<Float_t> trace_xn = {fReader, "wxn"};
TTreeReaderArray<Float_t> trace_xn_time = {fReader, "wxnT"};
TTreeReaderArray<Float_t> trace_xn_rise = {fReader, "wxnR"};
TTreeReaderArray<Float_t> trace_rdt = {fReader, "wrdt"};
TTreeReaderArray<Float_t> trace_rdt_time = {fReader, "wrdtT"};
TTreeReaderArray<Float_t> trace_rdt_rise = {fReader, "wrdtR"};
MonAnalyzer(TTree * /*tree*/ =0) { }
~MonAnalyzer() override { }
Int_t Version() const override { return 2; }
void Begin(TTree *tree) override;
void SlaveBegin(TTree *tree) override;
void Init(TTree *tree) override;
bool Notify() override;
bool Process(Long64_t entry) override;
Int_t GetEntry(Long64_t entry, Int_t getall = 0) override { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
void SetOption(const char *option) override { fOption = option; }
void SetObject(TObject *obj) override { fObject = obj; }
void SetInputList(TList *input) override { fInput = input; }
TList *GetOutputList() const override { return fOutput; }
void SlaveTerminate() override;
void Terminate() override;
ClassDefOverride(MonAnalyzer,0);
};
#endif
#ifdef MonAnalyzer_cxx
void MonAnalyzer::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the reader is initialized.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
fReader.SetTree(tree);
}
bool MonAnalyzer::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return true;
}
#endif // #ifdef MonAnalyzer_cxx

6
working/Monitor.C
View File

@ -317,12 +317,6 @@ Bool_t Monitor::Process(Long64_t entry){
} }
if (skipFlag ) continue; if (skipFlag ) continue;
//@==================== Basic gate
if( TMath::IsNaN(e[id]) ) continue ;
///if( ring[id] < -100 || ring[id] > 100 ) continue;
///if( ring[id] > 300 ) continue;
if( TMath::IsNaN(xn[id]) && TMath::IsNaN(xf[id]) ) continue ;
//@==================== Calibrations go here //@==================== Calibrations go here
if( corr->xnCorr.size() >= id && corr->xfxneCorr.size() >= id ) xnCal[id] = xn[id] * corr->xnCorr[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; if( corr->xnCorr.size() >= id && corr->xfxneCorr.size() >= id ) xnCal[id] = xn[id] * corr->xnCorr[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0];
if( corr->xfxneCorr.size() >= id ) xfCal[id] = xf[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; if( corr->xfxneCorr.size() >= id ) xfCal[id] = xf[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0];

56
working/test.C
View File

@ -1,9 +1,18 @@
#include "../Armory/ClassDetGeo.h" #include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h" #include "../Armory/ClassReactionConfig.h"
#include "../Armory/ClassCorrParas.h"
#include "../Cleopatra/ClassHelios.h" #include "../Cleopatra/ClassHelios.h"
#include "../Cleopatra/ClassTransfer.h" #include "../Cleopatra/ClassTransfer.h"
#include "ClassMonPlotter.h" #include "ClassMonPlotter.h"
#include "TFile.h"
#include "TChain.h"
#include "TH1F.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TTreeReaderArray.h"
#include "TClonesArray.h"
#include "TGraph.h"
void test(){ void test(){
@ -48,25 +57,48 @@ void test(){
// delete helios; // delete helios;
// delete transfer;root // delete transfer;root
DetGeo dd("detectorGeo.txt"); // DetGeo dd("detectorGeo.txt");
MonPlotter * pp = new MonPlotter(0, &dd, 8); // MonPlotter * pp = new MonPlotter(0, &dd, 8);
pp->SetUpCanvas("haha", 500, 3, 2); // pp->SetUpCanvas("haha", 500, 3, 2);
int rawEnergyRange[2] = { 100, 4000}; /// share with e, xf, xn // int rawEnergyRange[2] = { 100, 4000}; /// share with e, xf, xn
int energyRange[2] = { 0, 10}; /// in the E-Z plot // int energyRange[2] = { 0, 10}; /// in the E-Z plot
int rdtDERange[2] = { 0, 80}; // int rdtDERange[2] = { 0, 80};
int rdtERange[2] = { 0, 80}; // int rdtERange[2] = { 0, 80};
double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV] // double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV]
int thetaCMRange[2] = {0, 80}; // int thetaCMRange[2] = {0, 80};
pp->SetUpHistograms(rawEnergyRange, energyRange, exRange, thetaCMRange, rdtDERange, rdtERange); // pp->SetUpHistograms(rawEnergyRange, energyRange, exRange, thetaCMRange, rdtDERange, rdtERange);
pp->Plot(); // pp->Plot();
delete pp; // delete pp;
// TChain *chain = new TChain("gen_tree");
// chain->Add("../root_data/gen_run043.root");
// // chain->Print();
// TTreeReader reader(chain);
// TTreeReaderArray<TGraph> trace = {reader, "trace"};
// ULong64_t processedEntries = 0;
// while (reader.Next()) {
// printf("%llu | %lu \n", processedEntries, trace.GetSize());
// for( int i = 0; i < trace.GetSize(); i++ ){
// printf( " %d| %d\n", i, trace.At(i).GetN());
// }
// processedEntries ++;
// if( processedEntries > 10 ) break;
// }
CorrParas * corr = new CorrParas();
} }