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

Compare commits

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

No commits in common. "19a567f8fc1e31597cba4143bf4bc156d914130a" and "5bb5f5895653d540dd830552f8ab07cbd8eed4ee" have entirely different histories.
19a567f8fc ... 5bb5f58956

10 changed files with 245 additions and 808 deletions
Show Stats Expand all files Collapse all files

1
.gitignore vendored
View File

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

103
Armory/ClassCorrParas.h
View File

@ -3,71 +3,20 @@
//************************************** Correction parameters;
class CorrParas {
private:
int defaultSize;
public:
CorrParas(int defaultSize = 100){
is_xn_OK = false;
is_xfxne_OK = false;
is_e_OK = false;
is_xScale_OK = false;
is_rdt_OK = false;
this->defaultSize = defaultSize;
CorrParas(){
xnCorr.clear();
xScale.clear();
xfxneCorr.clear();
eCorr.clear();
rdtCorr.clear();
};
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> 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>> 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
//~========================================= xf = xn correction
@ -79,17 +28,15 @@ public:
if( file.is_open() ){
float a;
while( file >> a ) xnCorr.push_back(a);
printf(".......... done. size:%zu\n", xnCorr.size());
is_xn_OK = true;
printf(".......... done.\n");
}else{
for( int i = 0; i < defaultSize; i++ ) xnCorr.push_back(1.0);
printf(".......... fail. xnCorr[1..99] = 1.0\n");
is_xn_OK = false;
printf(".......... fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xnCorr.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
}
//~========================================= X-Scale correction
void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_scaleX.dat"){
printf(" loading x-Scale correction.");
@ -97,14 +44,11 @@ public:
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a;
float a, b;
while( file >> a ) xScale.push_back(a);
printf("........ done. size:%zu\n", xScale.size());
is_xScale_OK = true;
printf("........ done.\n");
}else{
for( int i = 0; i < defaultSize; i++ ) xScale.push_back(1.0);
printf("........ fail. xScale[1..99] = 1.0\n");
is_xScale_OK = false;
printf("........ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
@ -119,12 +63,9 @@ public:
if( file.is_open() ){
float a, b;
while( file >> a >> b) xfxneCorr.push_back({a, b});
printf("........ done. size:%zu\n", xfxneCorr.size());
is_xfxne_OK = true;
printf("........ done.\n");
}else{
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;
printf("........ fail.\n");
}
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]);
@ -139,12 +80,9 @@ public:
if( file.is_open() ){
float a, b;
while( file >> a >> b) eCorr.push_back( {a, b} ); // 1/a1, a0 , e' = e * a1 + a0
printf(".............. done. size:%zu\n", eCorr.size());
is_e_OK = true;
printf(".............. done.\n");
}else{
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;
printf(".............. fail.\n");
}
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]);
@ -159,12 +97,9 @@ public:
if( file.is_open() ){
float a, b;
while( file >> a >> b) rdtCorr.push_back({a, b});
printf("............ done. size:%zu\n", rdtCorr.size());
is_rdt_OK = true;
printf("............ done.\n");
}else{
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;
printf("............ fail.\n");
}
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]);

7
Armory/ClassDetGeo.h
View File

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

4
Armory/GeneralSort.h
View File

@ -182,7 +182,7 @@ void GeneralSort::SetUpTree(){
newSaveTree->SetDirectory(saveFile);
newSaveTree->AutoSave();
newSaveTree->Branch( "evID", &evID, "evID/l"); // simply copy
newSaveTree->Branch( "evID", &evID, "EventID/l"); // simply copy
eE = new Float_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]+"_t").c_str(), eT[i], Form("%s_t[%d]/l", 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]));
}

8
Armory/Process_Download
View File

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

248
working/ClassMonPlotter.h
View File

@ -9,8 +9,6 @@
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TStyle.h"
/******************************************************************
* This is Plotter for Monitor.C. It contains
@ -55,23 +53,10 @@ public:
~MonPlotter();
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 PlotRaw(bool isLog = false);
void PlotCal();
void PlotEZ();
void PlotEx();
TCanvas * canvas;
//====================== Histograms
@ -88,7 +73,6 @@ public:
TH2F ** hxfCal_xnCal;
TH2F ** he_xsCal; // raw e vs xf
TH2F ** he_x; // raw e vs x
TH2F * heCal_ID;
//===== eCal V z
TH2F * heCal_z;
@ -98,14 +82,8 @@ public:
TH2F * hrdt_ID;
TH1F ** hrdt; // single recoil
TH1I * hrdtMulti;
TH2F ** hrdt2D;
TH2F ** hrdt2Dg; // gated
//====== tDiff
TH1F * htDiff;
TH1F * htDiffg;
TH2F ** hrdt2Dg;
//====== Ex data
TH1F * hEx;
@ -118,9 +96,6 @@ public:
TH2F * hEx_ThetaCM;
//=======================
//======= Recoil Cut
TObjArray * cutList;
private:
unsigned short aID;
@ -130,7 +105,6 @@ private:
float recoilOutter;
double zRange[2] ; // zMin, zMax
TString canvasTitle;
TString suffix;
int numPad;
@ -155,7 +129,6 @@ MonPlotter::MonPlotter(unsigned short arrayID, DetGeo * detGeo, int numRDT){
zRange[1] = detGeo->array[aID].zMax + 50;
canvas = nullptr;
cutList = nullptr;
}
@ -170,19 +143,13 @@ MonPlotter::~MonPlotter(){
delete hxn_ID;
delete hArrayMulti;
delete heCal_ID;
delete heCal_zGC;
delete heCal_z;
delete hEx_ThetaCM;
delete hExCut1;
delete hExCut2;
delete heCal_zGC;
delete hrdt_ID;
delete hrdtMulti;
delete htDiff;
delete htDiffg;
for( int i = 0; i < numDet ; i++ ){
delete he[i];
@ -218,17 +185,13 @@ MonPlotter::~MonPlotter(){
delete [] hrdt2D;
delete [] hrdt2Dg;
delete cutG;
delete cutList;
}
void MonPlotter::SetUpCanvas(TString title, int padSize, int divX, int divY){
canvas = new TCanvas("canavs" + suffix, title, 500 * aID, 0, divX * padSize, divY * padSize);
canvas = new TCanvas("canavs" + suffix, title, 200 * aID, 200 * aID, divX * padSize, divY * padSize);
canvas->Divide(divX, divY);
numPad = divX * divY;
canvasTitle = title;
}
@ -259,13 +222,7 @@ template<typename T> void MonPlotter::CreateListOfHist2D(T ** &histList,
}
}
void MonPlotter::SetUpHistograms(int * rawEnergyRange,
int * energyRange,
double * exRange,
int * thetaCMRange,
int * rdtDERange,
int * rdtERange,
int * coinTimeRange){
void MonPlotter::SetUpHistograms(int * rawEnergyRange, int * energyRange, double * exRange, int * thetaCMRange, int * rdtDERange, int * rdtERange){
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]);
@ -280,14 +237,14 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange,
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_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]);
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]);
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);
//====================== 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_zGC = new TH2F("heCal_zGC" + suffix ,"E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]);
@ -298,8 +255,6 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange,
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];
hrdt2D = new TH2F * [numRDT/2];
hrdt2Dg = new TH2F * [numRDT/2];
@ -316,11 +271,6 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange,
}
}
//===================== 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
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]);
@ -347,187 +297,9 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange,
}
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;
}
}
}
void MonPlotter::LoadRDTGate(TString rdtCutFile){
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();
for( int i = 0; i < numPad; i++ ){
canvas->cd(i+1);
}
}

513
working/MonAnalyzer.C
View File

@ -1,443 +1,88 @@
#include "../Armory/AnalysisLib.h"
#include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h"
#include "../Armory/ClassCorrParas.h"
// #include "../Cleopatra/ClassHelios.h"
#include "../Cleopatra/ClassTransfer.h"
#include "ClassMonPlotter.h"
#include "Mapping.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"
#include "TH2.h"
#include "TStyle.h"
#include "TStopwatch.h"
#include "TMath.h"
#include "vector"
//^############################################ User setting
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};
double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV]
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;
}
//*==========================================
TransferReaction * transfer = new TransferReaction[numGeo];
int tempCount = 0;
for( int i = 0; i < (int) detGeo->array.size() ; i++){
if( !detGeo->array[i].enable ) continue;
transfer[tempCount].SetReactionFromFile("reactionConfig.txt", i);
tempCount ++;
}
//*==========================================
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");
#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.
// The following methods are defined in this file:
// Begin(): called every time a loop on the tree starts,
// 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 "MonAnalyzer.h"
#include <TH2.h>
#include <TStyle.h>
void MonAnalyzer::Begin(TTree * /*tree*/)
{
// 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();
}
void MonAnalyzer::SlaveBegin(TTree * /*tree*/)
{
// The SlaveBegin() function is called after the Begin() function.
// When running with PROOF SlaveBegin() is called on each slave server.
// The tree argument is deprecated (on PROOF 0 is passed).
TString option = GetOption();
}
//%=============================================
void raw(bool isLog = false, int arrayID = -1){
if( arrayID < 0 ){
for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotRaw(isLog);
}else{
if( arrayID < numGeo) plotter[arrayID]->PlotRaw(isLog);
}
bool MonAnalyzer::Process(Long64_t entry)
{
// The Process() function is called for each entry in the tree (or possibly
// keyed object in the case of PROOF) to be processed. The entry argument
// specifies which entry in the currently loaded tree is to be processed.
// When processing keyed objects with PROOF, the object is already loaded
// 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 cal(int arrayID = -1){
if( arrayID < 0 ){
for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotCal();
}else{
if( arrayID < numGeo) plotter[arrayID]->PlotCal();
}
void MonAnalyzer::SlaveTerminate()
{
// The SlaveTerminate() function is called after all entries or objects
// have been processed. When running with PROOF SlaveTerminate() is called
// on each slave server.
}
void ez(int arrayID = -1){
if( arrayID < 0 ){
for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotEZ();
}else{
if( arrayID < numGeo) plotter[arrayID]->PlotEZ();
}
}
void MonAnalyzer::Terminate()
{
// The Terminate() function is the last function to be called during
// a query. It always runs on the client, it can be used to present
// the results graphically or save the results to file.
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 Normal file
View File

@ -0,0 +1,115 @@
//////////////////////////////////////////////////////////
// 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,6 +317,12 @@ Bool_t Monitor::Process(Long64_t entry){
}
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
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];

56
working/test.C
View File

@ -1,18 +1,9 @@
#include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h"
#include "../Armory/ClassCorrParas.h"
#include "../Cleopatra/ClassHelios.h"
#include "../Cleopatra/ClassTransfer.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(){
@ -57,48 +48,25 @@ void test(){
// delete helios;
// delete transfer;root
// DetGeo dd("detectorGeo.txt");
// MonPlotter * pp = new MonPlotter(0, &dd, 8);
DetGeo dd("detectorGeo.txt");
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 energyRange[2] = { 0, 10}; /// in the E-Z plot
int rawEnergyRange[2] = { 100, 4000}; /// 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 rdtDERange[2] = { 0, 80};
int rdtERange[2] = { 0, 80};
// double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV]
// int thetaCMRange[2] = {0, 80};
double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV]
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();
}