From 242366792ca85c5d02c4dc244b466dd2307bae65 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS-MAC" Date: Thu, 6 Apr 2023 19:05:20 -0400 Subject: [PATCH] added AnalysisLib and AutoFit namespace, Greatly improve Monitor.C/h --- armory/AnalysisLib.h | 371 ++++++++++++++-- armory/AutoFit.C | 57 ++- armory/GeneralSortAgent.C | 1 + armory/Monitor_Util.C | 481 ++++++++++++++++++++ working/Monitors.C | 892 ++++++++++++++++++++++++++++++++++++++ working/Monitors.h | 368 ++++++++++++++++ 6 files changed, 2113 insertions(+), 57 deletions(-) create mode 100644 armory/Monitor_Util.C create mode 100644 working/Monitors.C create mode 100644 working/Monitors.h diff --git a/armory/AnalysisLib.h b/armory/AnalysisLib.h index ae16de5..99bebb2 100644 --- a/armory/AnalysisLib.h +++ b/armory/AnalysisLib.h @@ -8,6 +8,43 @@ #include #include +namespace AnalysisLib { + +std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector output; + + size_t pos; + do{ + pos = tempLine.find(splitter); /// fine splitter + if( pos == 0 ){ ///check if it is splitter again + tempLine = tempLine.substr(pos+1); + continue; + } + + std::string secStr; + if( pos == std::string::npos ){ + secStr = tempLine; + }else{ + secStr = tempLine.substr(0, pos+shift); + tempLine = tempLine.substr(pos+shift); + } + + ///check if secStr is begin with space + while( secStr.substr(0, 1) == " ") secStr = secStr.substr(1); + + ///check if secStr is end with space + while( secStr.back() == ' ') secStr = secStr.substr(0, secStr.size()-1); + + output.push_back(secStr); + ///printf(" |%s---\n", secStr.c_str()); + + }while(pos != std::string::npos ); + + return output; +} + + struct DetGeo{ double Bfield; /// T @@ -76,40 +113,6 @@ struct ReactionConfig{ }; -std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ - - std::vector output; - - size_t pos; - do{ - pos = tempLine.find(splitter); /// fine splitter - if( pos == 0 ){ ///check if it is splitter again - tempLine = tempLine.substr(pos+1); - continue; - } - - std::string secStr; - if( pos == std::string::npos ){ - secStr = tempLine; - }else{ - secStr = tempLine.substr(0, pos+shift); - tempLine = tempLine.substr(pos+shift); - } - - ///check if secStr is begin with space - while( secStr.substr(0, 1) == " ") secStr = secStr.substr(1); - - ///check if secStr is end with space - while( secStr.back() == ' ') secStr = secStr.substr(0, secStr.size()-1); - - output.push_back(secStr); - ///printf(" |%s---\n", secStr.c_str()); - - }while(pos != std::string::npos ); - - return output; -} - ///Using TMacro to load the detectorGeo frist, ///this indrect method is good for loading detectorGeo from TMacro in root file DetGeo LoadDetectorGeo(TMacro * macro){ @@ -306,7 +309,304 @@ void PrintReactionConfig(ReactionConfig reaction){ } +DetGeo detGeo; +ReactionConfig reactionConfig; +void LoadDetGeoAndReactionConfigFile(string detGeoFileName = "detectorGeo.txt", string reactionConfigFileName = "reactionConfig.txt"){ + printf("=======================\n"); + printf(" loading detector geometery : %s.", detGeoFileName.c_str()); + TMacro * haha = new TMacro(); + if( haha->ReadFile(detGeoFileName.c_str()) > 0 ) { + detGeo = AnalysisLib::LoadDetectorGeo(haha); + printf("... done.\n"); + AnalysisLib::PrintDetGeo(detGeo); + }else{ + printf("... fail\n"); + } + + printf("=======================\n"); + printf(" loading reaction config : %s.", reactionConfigFileName.c_str()); + TMacro * kaka = new TMacro(); + if( kaka->ReadFile(reactionConfigFileName.c_str()) > 0 ) { + reactionConfig = AnalysisLib::LoadReactionConfig(kaka); + printf("..... done.\n"); + AnalysisLib::PrintReactionConfig(reactionConfig); + }else{ + printf("..... fail\n"); + } + + delete haha; + delete kaka; + +} + +//************************************** Correction parameters; + +std::vector xnCorr; //correction of xn to match xf +std::vector xScale; // correction of x to be (0,1) +std::vector> xfxneCorr; //correction of xn and xf to match e +std::vector> eCorr; // correction to e, ch -> MeV +std::vector> rdtCorr; // correction of rdt, ch -> MeV + +//~========================================= xf = xn correction +void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.dat"){ + printf(" loading xf-xn correction."); + xnCorr.clear(); + ifstream file; + file.open(fileName); + if( file.is_open() ){ + float a; + while( file >> a ) xnCorr.push_back(a); + printf(".......... done.\n"); + }else{ + 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."); + xScale.clear(); + ifstream file; + file.open(fileName); + if( file.is_open() ){ + float a, b; + while( file >> a ) xScale.push_back(a); + printf("........ done.\n"); + }else{ + printf("........ fail.\n"); + } + file.close(); + if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]); +} + +//~========================================= e = xf + xn correction +void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xfxn_e.dat"){ + printf(" loading xf/xn-e correction."); + xfxneCorr.clear(); + ifstream file; + file.open(fileName); + if( file.is_open() ){ + float a, b; + while( file >> a >> b) xfxneCorr.push_back({a, b}); + printf("........ done.\n"); + }else{ + 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]); +} + +//~========================================= e correction +void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat"){ + printf(" loading e correction."); + eCorr.clear(); + ifstream file; + file.open(fileName); + 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.\n"); + }else{ + 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]); +} + +//~========================================= rdt correction +void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.dat"){ + printf(" loading rdt correction."); + rdtCorr.clear(); + ifstream file; + file.open(fileName); + if( file.is_open() ){ + float a, b; + while( file >> a >> b) rdtCorr.push_back({a, b}); + printf("............ done.\n"); + }else{ + 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]); +} + +double q, alpha, Et, betRel, gamm, G, massB, mass; //variables for Ex calculation +bool hasReactionPara = false; + +//~========================================= reaction parameters +void LoadReactionParas(bool verbose = false){ + + //check is the transfer.root is using the latest reactionConfig.txt + //sicne reaction.dat is generated as a by-product of transfer.root + //TFile * transfer = new TFile("transfer.root"); + //TString aaa1 = ""; + //TString aaa2 = ""; + //if( transfer->IsOpen() ){ + // TMacro * reactionConfig = (TMacro *) transfer->FindObjectAny("reactionConfig"); + // TMacro presentReactionConfig ("reactionConfig.txt"); + // aaa1 = ((TMD5*) reactionConfig->Checksum())->AsString(); + // aaa2 = ((TMD5*) presentReactionConfig.Checksum())->AsString(); + //} + //printf("%s\n", aaa1.Data()); + //printf("%s\n", aaa2.Data()); + + //if( aaa1 != aaa2 ) { + // printf("########################## recalculate transfer.root \n"); + // system("../Cleopatra/Transfer"); + // printf("########################## transfer.root updated\n"); + //} + printf(" loading reaction parameters"); + ifstream file; + file.open("reaction.dat"); + hasReactionPara = false; + if( file.is_open() ){ + string x; + int i = 0; + while( file >> x ){ + if( x.substr(0,2) == "//" ) continue; + if( i == 0 ) mass = atof(x.c_str()); + if( i == 1 ) q = atof(x.c_str()); + if( i == 2 ) betRel = atof(x.c_str()); + if( i == 3 ) Et = atof(x.c_str()); + if( i == 4 ) massB = atof(x.c_str()); + i = i + 1; + } + printf("........ done.\n"); + + hasReactionPara = true; + alpha = 299.792458 * abs(detGeo.Bfield) * q / TMath::TwoPi()/1000.; //MeV/mm + gamm = 1./TMath::Sqrt(1-betRel*betRel); + G = alpha * gamm * betRel * detGeo.detPerpDist ; + + if( verbose ){ + printf("\tmass-b : %f MeV/c2 \n", mass); + printf("\tcharge-b : %f \n", q); + printf("\tE-total : %f MeV \n", Et); + printf("\tmass-B : %f MeV/c2 \n", massB); + printf("\tbeta : %f \n", betRel); + printf("\tB-field : %f T \n", detGeo.Bfield); + printf("\tslope : %f MeV/mm \n", alpha * betRel); + printf("\tdet radius: %f mm \n", detGeo.detPerpDist); + printf("\tG-coeff : %f MeV \n", G); + printf("=================================\n"); + } + + }else{ + printf("........ fail.\n"); + hasReactionPara = false; + } + file.close(); + +} + +std::vector CalExTheta(double e, double z){ + if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()}; + + double y = e + mass; // to give the KE + mass of proton; + double Z = alpha * gamm * betRel * z; + double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ; + + if( TMath::Abs(Z) < H ) { + ///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi) + double tolerrence = 0.001; + double phi = 0; ///initial phi = 0 -> ensure the solution has f'(phi) > 0 + double nPhi = 0; /// new phi + + int iter = 0; + do{ + phi = nPhi; + nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2)); + iter ++; + if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break; + }while( TMath::Abs(phi - nPhi ) > tolerrence); + phi = nPhi; + + /// check f'(phi) > 0 + double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2); + if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){ + double K = H * TMath::Sin(phi); + double x = TMath::ACos( mass / ( y * gamm - K)); + double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame + double EB = TMath::Sqrt(mass*mass + Et*Et - 2*Et*TMath::Sqrt(momt*momt + mass * mass)); + Ex = EB - massB; + + double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y; + double hahaha2 = gamm* betRel * momt; + thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg(); + + }else{ + Ex = TMath::QuietNaN(); + thetaCM = TMath::QuietNaN(); + } + }else{ + Ex = TMath::QuietNaN(); + thetaCM = TMath::QuietNaN(); + } + + return {Ex, thetaCM}; + +} + +//************************************** TCutG + +TObjArray * LoadListOfTCut(TString fileName, TString cutName = "cutList"){ + + TObjArray * cutList = nullptr; + + TFile * fCut = new TFile(fileName); + bool isCutFileOpen = fCut->IsOpen(); + if(!isCutFileOpen) { + printf( "Failed to open rdt-cutfile 1 : %s\n" , fileName.Data()); + }else{ + cutList = (TObjArray *) fCut->FindObjectAny(cutName); + + 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() + ); + } + } + } + + return cutList; +} + +TCutG * LoadSingleTCut( TString fileName, TString cutName = "cutEZ"){ + + TCutG * cut = nullptr; + + TFile * fCut = new TFile(fileName); + bool isCutFileOpen = fCut->IsOpen(); + if( !isCutFileOpen) { + printf( "Failed to open E-Z cutfile : %s\n" , fileName.Data()); + }else{ + cut = (TCutG *) fCut->FindObjectAny(cutName); + if( cut != NULL ) { + printf("Found EZ cut| name : %s, VarX: %s, VarY: %s, numPoints: %d \n", + cut->GetName(), + cut->GetVarX(), + cut->GetVarY(), + cut->GetN() + ); + } + } + + return cut; +} + + +//************************************** Others std::vector> combination(std::vector arr, int r){ std::vector> output; @@ -477,7 +777,6 @@ std::vector> FindMatchingPair(std::vector enX, std:: } - printf("fitEnergy = ");for( int k = 0; k < (int) fitEnergy.size() ; k++){ printf("%7.2f, ", fitEnergy[k]); }; printf("\n"); printf("refEnergy = ");for( int k = 0; k < (int) refEnergy.size() ; k++){ printf("%7.2f, ", refEnergy[k]); }; printf("\n"); @@ -490,4 +789,6 @@ std::vector> FindMatchingPair(std::vector enX, std:: } +} + #endif \ No newline at end of file diff --git a/armory/AutoFit.C b/armory/AutoFit.C index 95e2fc5..4b94781 100644 --- a/armory/AutoFit.C +++ b/armory/AutoFit.C @@ -11,16 +11,28 @@ #define AutoFit_C #include +#include +#include #include #include #include #include #include #include +#include +#include +#include +#include + #include +#include +#include +#include +#include + +namespace AutoFit{ //Global fit paramaters - std::vector BestFitMean; std::vector BestFitCount; std::vector BestFitSigma; @@ -91,9 +103,9 @@ TColor RGBWheel(double ang){ ang = ang * TMath::DegToRad(); - double r = max(0., (1+2*cos(ang))/3.); - double g = max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.); - double b = max(0., (1 - cos(ang) - sqrt(3)* sin(ang))/3.); + double r = std::max(0., (1+2*cos(ang))/3.); + double g = std::max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.); + double b = std::max(0., (1 - cos(ang) - sqrt(3)* sin(ang))/3.); TColor col(r,g,b); @@ -273,8 +285,8 @@ void GoodnessofFit(TH1F * hist, TF1 * fit){ printf("################################################\n"); } -vector energy, height, sigma, lowE, highE ; -vector energyFlag, sigmaFlag; +std::vector energy, height, sigma, lowE, highE ; +std::vector energyFlag, sigmaFlag; bool loadFitParameters(TString fitParaFile){ @@ -289,7 +301,7 @@ bool loadFitParameters(TString fitParaFile){ printf("====================================================================== \n"); printf("----- loading fit parameters from : %s", fitParaFile.Data()); - ifstream file; + std::ifstream file; file.open(fitParaFile.Data()); if( !file){ @@ -298,7 +310,7 @@ bool loadFitParameters(TString fitParaFile){ } while( file.good()) { - string tempLine; + std::string tempLine; getline(file, tempLine); if( tempLine.substr(0, 1) == "#" ) continue; @@ -307,7 +319,7 @@ bool loadFitParameters(TString fitParaFile){ ///printf("%s\n", tempLine.c_str()); - vector temp = SplitStrAF(tempLine, " "); + std::vector temp = SplitStrAF(tempLine, " "); if( temp.size() < 7 ) continue; @@ -491,7 +503,7 @@ void fitGaussPol(TH1F * hist, double mean, double sigmaMax, int degPol, double x //######################################## //#### fit 2 gauss + pol-1 // not updated //######################################## -vector fit2GaussP1(TH1F * hist, double mean1, double sigma1, +std::vector fit2GaussP1(TH1F * hist, double mean1, double sigma1, double mean2, double sigma2, double xMin, double xMax, TString optStat = "", bool newCanvas = false){ @@ -503,7 +515,7 @@ vector fit2GaussP1(TH1F * hist, double mean1, double sigma1, recentFitMethod = "fit2GaussP1"; - vector output; + std::vector output; output.clear(); gStyle->SetOptStat(optStat); @@ -684,7 +696,7 @@ void fitGF3Pol(TH1F * hist, double mean, double sigmaMax, double ratio, double b //GoodnessofFit(hist, fit); /// 0 1 2 3 4 5 - string label[8] = {"Area", "mean", "sigma", "ratio", "beta", "step"}; + std::string label[8] = {"Area", "mean", "sigma", "ratio", "beta", "step"}; printf("---------- The detail\n"); for(int i = 0 ; i < 6 ; i++){ printf("%d | %8s | %f (%f) \n", i, label[i].c_str(), paraA[i], paraE[i]); @@ -749,7 +761,7 @@ void fitGF3Pol(TH1F * hist, double mean, double sigmaMax, double ratio, double b //############################################## //##### Auto Fit n-Gauss with estimated BG //############################################## -vector fitAuto(TH1F * hist, int bgEst = 10, +std::vector fitAuto(TH1F * hist, int bgEst = 10, double peakThreshold = 0.05, double sigmaMax = 0, int peakDensity = 4, @@ -817,7 +829,7 @@ vector fitAuto(TH1F * hist, int bgEst = 10, int * inX = new int[nPeaks]; TMath::Sort(nPeaks, xpos, inX, 0 ); - vector energy, height; + std::vector energy, height; for( int j = 0; j < nPeaks; j++){ energy.push_back(xpos[inX[j]]); height.push_back(ypos[inX[j]]); @@ -930,7 +942,7 @@ vector fitAuto(TH1F * hist, int bgEst = 10, double bw = specS->GetBinWidth(1); - vector exPos; + std::vector exPos; for(int i = 0; i < nPeaks ; i++){ exPos.push_back(paraA[numParPerPeak*i+1]); printf("%2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", @@ -1004,7 +1016,7 @@ vector fitAuto(TH1F * hist, int bgEst = 10, //######################################## //###### NOT tested //######################################## -vector fitNGF3(TH1 * hist, int bgEst = 10, +std::vector fitNGF3(TH1 * hist, int bgEst = 10, double peakThreshold = 0.1, double sigmaMax = 20, int peakDensity = 4, @@ -1067,7 +1079,7 @@ vector fitNGF3(TH1 * hist, int bgEst = 10, int * inX = new int[nPeaks]; TMath::Sort(nPeaks, xpos, inX, 0 ); - vector energy, height; + std::vector energy, height; for( int j = 0; j < nPeaks; j++){ energy.push_back(xpos[inX[j]]); height.push_back(ypos[inX[j]]); @@ -1137,7 +1149,7 @@ vector fitNGF3(TH1 * hist, int bgEst = 10, double bw = specS->GetBinWidth(1); - vector exPos; + std::vector exPos; for(int i = 0; i < nPeaks ; i++){ exPos.push_back(paraA[numParPerPeak*i+1]); double totCount = paraA[numParPerPeak*i] + paraA[numParPerPeak*i+3]; @@ -1987,10 +1999,10 @@ void fitNGaussPolSub(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.t int nClick = 0; bool peakFlag = 1; -vector xPeakList; -vector yPeakList; -vector xBGList; -vector yBGList; +std::vector xPeakList; +std::vector yPeakList; +std::vector xBGList; +std::vector yBGList; TH1F * tempHist; int markerStyle = 23; @@ -2844,5 +2856,6 @@ void fitSpecial(TH1F * hist, TString fitFile = "AutoFit_para.txt"){ } +} #endif diff --git a/armory/GeneralSortAgent.C b/armory/GeneralSortAgent.C index 4cf9489..16c8e10 100644 --- a/armory/GeneralSortAgent.C +++ b/armory/GeneralSortAgent.C @@ -49,5 +49,6 @@ void GeneralSortAgent(Int_t runNum, int nWorker = 1, int traceMethod = -1){ f1->Close(); f2->Close(); + delete chain; } \ No newline at end of file diff --git a/armory/Monitor_Util.C b/armory/Monitor_Util.C new file mode 100644 index 0000000..c1437e0 --- /dev/null +++ b/armory/Monitor_Util.C @@ -0,0 +1,481 @@ +#ifndef Utilities +#define Utilities + +#include +#include +#include "../working/Mapping.h" + +//This file runs after on Monitor.C +//This file is parasite on Monitor.C + +int canvasSize[2] = {2000, 1200}; + +void listDraws(void) { + printf("------------------- List of Plots -------------------\n"); + printf(" newCanvas() - Create a new Canvas\n"); + printf("-----------------------------------------------------\n"); + printf(" rawID() - Raw \033[0;31me\033[0m, \033[0;31mring\033[0m, \033[0;31mxf\033[0m, \033[0;31mxn\033[0m vs detID\n"); + printf(" rawe() - Raw \033[0;31me\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" rawxf() - Raw \033[0;31mxf\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" rawxn() - Raw \033[0;31mxn\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" xfVxn() - Raw \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" eVxs() - Raw \033[0;31me\033[0m vs. Raw \033[0;31mxs = xf + xn\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" eVx() - Raw \033[0;31me\033[0m vs. RAW \033[0;31mx\033[0m for all %d detectors\n", mapping::NARRAY); + printf("-----------------------------------------------------\n"); + printf(" eVxsCal() - Raw \033[0;31me\033[0m vs. Corrected \033[0;31mxs\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" ecal() - Calibrated \033[0;31me\033[0m for all %d detectors\n", mapping::NARRAY); + printf(" ecal2() - Calibrated \033[0;31me\033[0m for all %d detectors (same row or same col)\n", mapping::NARRAY); + printf("xfCalVxnCal() - Calibrated \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", mapping::NARRAY); + printf("-----------------------------------------------------\n"); + printf(" eCalVxCal() - Cal \033[0;31me\033[0m vs. \033[0;31mx\033[0m for all %d detectors\n", mapping::NARRAY); + printf("-----------------------------------------------------\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(" eCalVz() - Energy vs. Z\n"); + printf(" eCalVzRow() - Energy vs. Z for each row\n"); + printf(" excite() - Excitation Energy\n"); + printf(" ExThetaCM() - Ex vs ThetaCM\n"); + printf(" ExVxCal() - Ex vs X for all %d detectors\n", mapping::NARRAY); + //printf(" eSVeRaw() - e(Ex,z) vs eRaw for all %d detectors\n", mapping::NARRAY); + 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"); +} + + +int xD, yD; +void FindBesCanvasDivision(int nPad){ + for( int i = TMath::Sqrt(nPad); i >= 2 ; i--){ + if( nPad % i == 0 ) { + yD = i; + xD = nPad/i; + break; + } + } +} + +int nCanvas=0; +void newCanvas(int sizeX = 800, int sizeY = 600, int posX = 0, int posY = 0){ + TString name; name.Form("cNewCanvas%d | %s", nCanvas, canvasTitle.Data()); + TCanvas * cNewCanvas = new TCanvas(name, name, posX, posY, sizeX, sizeY); + nCanvas++; + cNewCanvas->cd(); +} + + +//TODO set histogram y-axis all the same heigh + + +void rawID(){ + TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); + if( cRawID == NULL ) cRawID = new TCanvas("cRawID", Form("Raw e, Ring, xf, xn vs ID | %s", canvasTitle.Data()), canvasSize[0], canvasSize[1]); + cRawID->Clear();cRawID->Divide(2,2); + cRawID->cd(1); cRawID->cd(1)->SetGrid(); heVID->Draw("colz"); + cRawID->cd(2); cRawID->cd(2)->SetGrid(); hMultiHit->Draw(); + cRawID->cd(3); cRawID->cd(3)->SetGrid(); hxfVID->Draw("colz"); + cRawID->cd(4); cRawID->cd(4)->SetGrid(); hxnVID->Draw("colz"); +} + +void rawe(Bool_t isLogy = false) { + TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE"); + if( cRawE == NULL ) cRawE = new TCanvas("cRawE",Form("E raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cRawE->Clear();cRawE->Divide(numCol,numRow); + for (Int_t i=0; i < mapping::NARRAY; i++) { + cRawE->cd(i+1); + cRawE->cd(i+1)->SetGrid(); + if( isLogy ) cRawE->cd(i+1)->SetLogy(); + he[i]->Draw(""); + } +} + +void rawxf(Bool_t isLogy = false) { + TCanvas *cRawXf = (TCanvas *) gROOT->FindObjectAny("cRawXf"); + if( cRawXf == NULL ) cRawXf = new TCanvas("cRawXf",Form("Xf raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cRawXf->Clear();cRawXf->Divide(numCol,numRow); + for (Int_t i=0; icd(i+1); + cRawXf->cd(i+1)->SetGrid(); + if( isLogy ) cRawXf->cd(i+1)->SetLogy(); + hxf[i]->Draw(""); + } +} + +void rawxn(Bool_t isLogy = false) { + TCanvas *cRawXn = (TCanvas *) gROOT->FindObjectAny("cRawXn"); + if( cRawXn == NULL ) cRawXn = new TCanvas("cRawXn",Form("Xn raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cRawXn->Clear();cRawXn->Divide(numCol,numRow); + for (Int_t i=0; icd(i+1); + cRawXn->cd(i+1)->SetGrid(); + if( isLogy ) cRawXn->cd(i+1)->SetLogy(); + hxn[i]->Draw(""); + } +} + +void xfVxn(void) { + TCanvas *cxfxn = (TCanvas *) gROOT->FindObjectAny("cxfxn"); + if( cxfxn == NULL ) cxfxn = new TCanvas("cxfxn",Form("XF vs. XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cxfxn->Clear(); cxfxn->Divide(numCol,numRow); + for (Int_t i=0;icd(i+1); + cxfxn->cd(i+1)->SetGrid(); + hxfVxn[i]->Draw("col"); + } +} + +void eVxs(void) { + TCanvas *cxfxne = (TCanvas *) gROOT->FindObjectAny("cxfxne"); + if( cxfxne == NULL ) cxfxne = new TCanvas("cxfxne",Form("E - XF+XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cxfxne->Clear(); cxfxne->Divide(numCol,numRow); + TLine line(0,0, 4000, 4000); line.SetLineColor(2); + for (Int_t i=0;icd(i+1); + cxfxne->cd(i+1)->SetGrid(); + heVxs[i]->Draw("col"); + line.Draw("same"); + } +} + +void ecal(void) { + TCanvas *cEC = (TCanvas *) gROOT->FindObjectAny("cEC"); + if(cEC == NULL) cEC = new TCanvas("cEC",Form("E corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cEC->Clear();cEC->Divide(numCol,numRow); + for (Int_t i=0; icd(i+1); + cEC->cd(i+1)->SetGrid(); + heCal[i]->Draw(""); + } + + TCanvas *cEC2 = (TCanvas *) gROOT->FindObjectAny("cEC2"); + if(cEC2 == NULL) cEC2 = new TCanvas("cEC2",Form("E corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cEC2->Clear(); + heCalID->Draw("colz"); +} + +void ecal2(void) { + TCanvas *cECall = (TCanvas *) gROOT->FindObjectAny("cECall"); + if(cECall == NULL) cECall = new TCanvas("cECall",Form("E corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + + int maxRC = TMath::Max(numRow, numCol); + cECall->Clear();cECall->Divide(maxRC,2); + //plot same position + for(int i = 0; i < numCol ; i++){ + cECall->cd(i+1); + cECall->cd(i+1)->SetGrid(); + heCal[i]->SetLineColor(1); heCal[i]->Draw(""); + for(int j = 1 ; j < numRow; j++){ + heCal[numCol*j + i]->SetLineColor(j+1); + heCal[numCol*j + i]->Draw("same"); + } + } + //plot same side + TH1F * heC2[mapping::NARRAY]; + for (Int_t i = 0; i< numRow; i++) { + cECall->cd(i+maxRC+1); + cECall->cd(i+maxRC+1)->SetGrid(); + heC2[numCol*i] = (TH1F* )heCal[numCol*i]->Clone(); + heC2[numCol*i]->SetLineColor(1); heC2[numCol*i]->Draw(""); + for( int j = 1; j < numCol; j++){ + heC2[numCol*i+j] = (TH1F* )heCal[numCol*i+j]->Clone(); + heC2[numCol*i+j]->SetLineColor(j+1); + heC2[numCol*i+j]->Draw("same"); + } + } +} + + +void xfCalVxnCal(void) { + TCanvas *cxfxnC = (TCanvas *) gROOT->FindObjectAny("cxfxnC"); + if(cxfxnC == NULL) cxfxnC = new TCanvas("cxfxnC",Form("XF vs XN corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cxfxnC->Clear(); cxfxnC->Divide(numCol,numRow); + for (Int_t i=0;icd(i+1); + cxfxnC->cd(i+1)->SetGrid(); + hxfCalVxnCal[i]->Draw("col"); + } +} + +void eVxsCal(void) { + TCanvas *cxfxneC = (TCanvas *) gROOT->FindObjectAny("cxfxneC"); + if(cxfxneC == NULL)cxfxneC = new TCanvas("cxfxneC",Form("Raw E - Corrected XF+XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cxfxneC->Clear(); cxfxneC->Divide(numCol,numRow); + TLine line(0,0, 4000, 4000); line.SetLineColor(2); + for (Int_t i=0;icd(i+1); + cxfxneC->cd(i+1)->SetGrid(); + heVxsCal[i]->Draw("col"); + line.Draw("same"); + } +} + +void eVx(void) { + TCanvas *ceVx = (TCanvas *) gROOT->FindObjectAny("ceVx"); + if(ceVx == NULL) ceVx = new TCanvas("ceVx",Form("E vs. X = (xf-xn)/e | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + ceVx->Clear(); ceVx->Divide(numCol,numRow); + for (Int_t i=0;icd(i+1); heVx[i]->Draw("col"); + } +} + +void eCalVxCal(void) { + TCanvas *cecalVxcal = (TCanvas *) gROOT->FindObjectAny("cecalVxcal"); + if( cecalVxcal == NULL ) cecalVxcal = new TCanvas("cecalVxcal",Form("ECALVXCAL | %s",canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cecalVxcal->Clear(); cecalVxcal->Divide(numCol,numRow); + for (Int_t i=0;icd(i+1); + heCalVxCal[i]->SetMarkerStyle(7); + heCalVxCal[i]->Draw(""); + } +} + +void eCalVxCalG(void) { + TCanvas *cecalVxcalG = (TCanvas *) gROOT->FindObjectAny("cecalVxcalG"); + if( cecalVxcalG == NULL ) cecalVxcalG = new TCanvas("cecalVxcalG",Form("ECALVXCAL | %s",canvasTitle.Data()),canvasSize[0], canvasSize[1]); + cecalVxcalG->Clear(); cecalVxcalG->Divide(numCol,numRow); + for (Int_t i=0;icd(i+1); + heCalVxCalG[i]->SetMarkerStyle(7); + heCalVxCalG[i]->Draw(""); + } +} + +// void elum(void) { +// TCanvas *celum = (TCanvas *) gROOT->FindObjectAny("celum"); +// if( celum == NULL ) celum = new TCanvas("celum",Form("ELUM | %s", canvasTitle.Data()),1000,1000); +// celum->Clear(); celum->Divide(4,4); +// for( int i = 0 ; i < 16 ; i++){ +// celum->cd(i+1); +// helum[i]->Draw(""); +// } + +// TCanvas *celumID = (TCanvas *) gROOT->FindObjectAny("celumID"); +// if( celumID == NULL ) celumID = new TCanvas("celumID",Form("ELUM-ID | %s", canvasTitle.Data()),1100, 0, 500,500); +// celumID->Clear(); +// helumID->Draw("colz"); + +// } + +// void apollo(void) { +// TCanvas *capollo = (TCanvas *) gROOT->FindObjectAny("capollo"); +// if( capollo == NULL ) capollo = new TCanvas("capollo",Form("APOLLO | %s", canvasTitle.Data()),1000,1000); +// capollo->Clear(); capollo->Divide(5,4); +// for( int i = 0 ; i < 20 ; i++){ +// capollo->cd(i+1); +// hApollo[i]->Draw(""); +// } +// } + + +void recoils(bool isLogz = false) { + TCanvas *crdt = (TCanvas *) gROOT->FindObjectAny("crdt"); + if( crdt == NULL ) crdt = new TCanvas("crdt",Form("raw RDT | %s", canvasTitle.Data()),1700, 0, 1000,1000); + crdt->Clear();crdt->Divide(2,2); + + if( isLogz ) crdt->cd(1)->SetLogz(); crdt->cd(1); hrdt2D[0]->Draw("col"); + if( isLogz ) crdt->cd(2)->SetLogz(); crdt->cd(2); hrdt2D[1]->Draw("col"); + if( isLogz ) crdt->cd(3)->SetLogz(); crdt->cd(3); hrdt2D[3]->Draw("col"); + if( isLogz ) crdt->cd(4)->SetLogz(); crdt->cd(4); hrdt2D[2]->Draw("col"); + + + //TCanvas *crdtSum = (TCanvas *) gROOT->FindObjectAny("crdtSum"); + //if( crdtSum == NULL ) crdtSum = new TCanvas("crdtSum",Form("raw RDT dE-Esum | %s", canvasTitle.Data()),100, 0, 1000,1000); + //crdtSum->Clear();crdtSum->Divide(2,2); + // + //if( isLogz ) crdtSum->cd(1)->SetLogz(); crdtSum->cd(1); hrdt2Dsum[0]->Draw("col"); + //if( isLogz ) crdtSum->cd(2)->SetLogz(); crdtSum->cd(2); hrdt2Dsum[1]->Draw("col"); + //if( isLogz ) crdtSum->cd(3)->SetLogz(); crdtSum->cd(3); hrdt2Dsum[3]->Draw("col"); + //if( isLogz ) crdtSum->cd(4)->SetLogz(); crdtSum->cd(4); hrdt2Dsum[2]->Draw("col"); + + TCanvas *crdtID = (TCanvas *) gROOT->FindObjectAny("crdtID"); + if( crdtID == NULL ) crdtID = new TCanvas("crdtID",Form("raw RDT ID | %s", canvasTitle.Data()),0,0, 500, 500); + crdtID->Clear(); + if( isLogz ) crdtID->SetLogz(); + hrdtID->Draw("colz"); + + TCanvas *crdtS = (TCanvas *) gROOT->FindObjectAny("crdtS"); + if( crdtS == NULL ) crdtS = new TCanvas("crdtS",Form("raw RDT | %s", canvasTitle.Data()),600, 0, 1000, 1000); + crdtS->Clear(); crdtS->Divide(2,4); + for( int i = 0; i < 8; i ++){ + crdtS->cd(i+1); + if( isLogz ) crdtS->cd(i+1)->SetLogy(); + hrdt[i]->Draw(""); + } + + //TCanvas *crdtTAC = (TCanvas *) gROOT->FindObjectAny("crdtTAC"); + //if( crdtTAC == NULL ) crdtTAC = new TCanvas("crdtTAC",Form("raw RDTtac | %s", canvasTitle.Data()),0,0, 1600, 1600); + //crdtTAC->Clear(); crdtTAC->Divide(2,4); + //for( int i = 0; i < 8; i ++){ + // crdtTAC->cd(i+1); + // htacRecoil[i]->Draw("colz"); + //} + //for( int i = 0; i < 4; i ++){ + // crdtTAC->cd(i+1+8); + // htacRecoilsum[i]->Draw("colz"); + //} + +} + +void eCalVz(void) { + TCanvas *cecalVz = (TCanvas *) gROOT->FindObjectAny("cecalVz"); + if( cecalVz == NULL ) cecalVz = new TCanvas("cevalVz",Form("ECALVZ : %s", canvasTitle.Data()),1000,650); + cecalVz->Clear(); cecalVz->Divide(2,1); + gStyle->SetOptStat("neiou"); + cecalVz->cd(1);heCalVz->Draw("col"); + cecalVz->cd(2);heCalVzGC->Draw("col"); +} + +void eCalVzRow() { + TCanvas *cecalVzRow = (TCanvas *) gROOT->FindObjectAny("cecalVzRow"); + if( cecalVzRow == NULL ) cecalVzRow = new TCanvas("cevalVzRow",Form("eCal - Z : %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); + FindBesCanvasDivision(numRow); + cecalVzRow->Clear(); cecalVzRow->Divide(xD,yD); + gStyle->SetOptStat("neiou"); + + for(int row = 0; row < numRow; row ++){ + cecalVzRow->cd(row+1); + cecalVzRow->cd(row+1)->SetGrid(); + hecalVzRow[row]->Draw("colz"); + } + +} + +void excite(void) { + TCanvas *cex = (TCanvas *) gROOT->FindObjectAny("cex"); + if( cex == NULL ) cex = new TCanvas("cex",Form("EX : %s", canvasTitle.Data()),0, 0, 1000,650); + cex->Clear(); + gStyle->SetOptStat("neiou"); + hEx->Draw(""); + + + TCanvas *cexI = (TCanvas *) gROOT->FindObjectAny("cexI"); + if( cexI == NULL ) cexI = new TCanvas("cexI",Form("EX : %s", canvasTitle.Data()),500, 0, 1600,1000); + cexI->Clear();cexI->Divide(numCol,numRow); + gStyle->SetOptStat("neiou"); + for( int i = 0; i < mapping::NARRAY; i++){ + cexI->cd(i+1); + hExi[i]->Draw(""); + } + +} + + +void ExThetaCM(void) { + TCanvas *cExThetaCM = (TCanvas *) gROOT->FindObjectAny("cExThetaCM"); + if( cExThetaCM == NULL ) cExThetaCM = new TCanvas("cExThetaCM",Form("EX - ThetaCM | %s", canvasTitle.Data()),650,650); + cExThetaCM->Clear(); + gStyle->SetOptStat("neiou"); + hExThetaCM->Draw("colz"); +} + +void ExVxCal(TString drawOpt = "") { + TCanvas *cExVxCal = (TCanvas *) gROOT->FindObjectAny("cExVxCal"); + if( cExVxCal == NULL ) cExVxCal = new TCanvas("cExVxCal",Form("EX | %s", canvasTitle.Data()),1600,1000); + cExVxCal->Clear(); + gStyle->SetOptStat("neiou"); + + cExVxCal->Divide(numCol,numRow); + for( int i = 0; i < mapping::NARRAY; i++){ + cExVxCal->cd(i+1); + if( drawOpt == "" )hExVxCal[i]->SetMarkerStyle(7); + hExVxCal[i]->Draw(drawOpt); + } + +} + + +//void eSVeRaw(void) { +// TCanvas *ceSVeRaw = (TCanvas *) gROOT->FindObjectAny("ceSVeRaw"); +// if( ceSVeRaw == NULL ) ceSVeRaw = new TCanvas("ceSVeRaw",Form("e(Ex,z) vs eRaw | %s", canvasTitle.Data()),1000,650); +// ceSVeRaw->Clear(); +// gStyle->SetOptStat("neiou"); +// +// ceSVeRaw->Divide(numCol,numRow); +// for( int i = 0; i < mapping::NARRAY; i++){ +// ceSVeRaw->cd(i+1); +// heSVeRaw[i]->SetMarkerStyle(7); +// heSVeRaw[i]->Draw(""); +// } +// +//} + +// void tac(void) { +// TCanvas *ctac = (TCanvas *) gROOT->FindObjectAny("ctac"); +// if( ctac == NULL ) ctac = new TCanvas("ctac",Form("ARRAY-RDT | %s", canvasTitle.Data() ),1000,650); +// ctac->Clear();ctac->SetGrid(0);ctac->Divide(numCol,numRow); +// gStyle->SetOptStat("neiou"); +// for (Int_t i=0;icd(i+1); htacArray[i]->Draw(""); +// // cutG = (TCutG *)cutList->At(i); +// // cutG->Draw("same"); +// } +// } + + +// void ic(){ + +// TCanvas *cic = (TCanvas *) gROOT->FindObjectAny("cic"); +// if( cic == NULL ) cic = new TCanvas("cic",Form("Ionization Chamber | %s", canvasTitle.Data() ),1200,800); + +// cic->Clear(); cic->SetGrid(0); cic->Divide(3,2); +// gStyle->SetOptStat("neiou"); + +// cic->cd(1); hic0->Draw(); +// cic->cd(2); hic1->Draw(); +// cic->cd(3); hic2->Draw(); +// cic->cd(4); hic01->Draw("colz"); +// cic->cd(5); hic02->Draw("colz"); +// cic->cd(6); hic12->Draw("colz"); + +// } + + + +// void Count1DH(TString name, TH1F * hist, TCanvas * canvas, int padID, double x1, double x2, Color_t color){ + +// int k1 = hist->FindBin(x1); +// int k2 = hist->FindBin(x2); + +// int hight = 0 ; +// for( int i = k1; i < k2 ; i ++){ +// int temp = hist->GetBinContent(i); +// if( temp > hight ) hight = temp; +// } +// hight = hight * 1.2; +// int max = hist->GetMaximum(); + +// canvas->cd(padID); + +// if( color != 0 ){ +// TBox box; +// box.SetFillColorAlpha(color, 0.1); +// box.DrawBox(x1, 0, x2, hight); +// } + +// int count = hist->Integral(k1, k2); + +// TLatex text; +// text.SetTextFont(82); +// text.SetTextSize(0.06); +// if( color != 0 ){ +// text.SetTextColor(color); +// text.DrawLatex(x1, hight, Form("%d", count)); +// }else{ +// text.DrawLatex((x1+x2)/2., max, Form("%d", count)); +// } + +// printf(" %s : %d \n", name.Data(), count); + +// } + + + +#endif diff --git a/working/Monitors.C b/working/Monitors.C new file mode 100644 index 0000000..94205bc --- /dev/null +++ b/working/Monitors.C @@ -0,0 +1,892 @@ +#define Monitors_cxx + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "../Cleopatra/Isotope.h" +#include "Mapping.h" + +using namespace std; + +//############################################ User setting +ULong64_t maxNumberEvent = 1000000000; + +//---histogram setting +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 thetaCMRange[2] = {0, 80}; + +double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV] + +int coinTimeRange[2] = { -200, 200}; +int timeRangeUser[2] = {0, 99999999}; /// min, use when cannot find time, this set the min and max + +bool isUseArrayTrace = false; +bool isUseRDTTrace = false; + +//---Gate +bool isTimeGateOn = true; +int timeGate[2] = {-20, 12}; /// min, max, 1 ch = 10 ns +double eCalCut[2] = {0.5, 50}; /// lower & higher limit for eCal +int dEgate[2] = { 500, 1500}; +int Eresgate[2] = { 1000, 4000}; +double thetaCMGate = 10; /// deg +double xGate = 0.9; ///cut out the edge +vector skipDetID = {11, 16, 23} ;//{2, 11, 17} + +TString rdtCutFile1 = ""; +TString rdtCutFile2 = ""; +TString ezCutFile = "";//"ezCut.root"; + +//############################################ end of user setting +/****************************************************************** +* variable and histogram naming rules * +* name are case sensitive, so as any C/C++ code * +* * +* ID is dettector ID * +* * +* raw data from gen_tree are e, xf, xn, ring. * +* the x from raw data is x * +* * +* xf + xn = xs, s for sum * +* * +* calibrated data are eCal, xfCal, xnCal, ringCal. * +* the x from cal data is xCal * +* * +* xfCal + xnCal = xsCal * +* * +* since the z is always from xCal, so it calls z. * +* * +* Excitation energy calls Ex * +* * +* * +* TH2D is always using "V" to seperate 2 variables, like eVx * +* * +* histogram with TCutG, add suffix "GC" for Graphical-Cut. * +* * +*******************************************************************/ +//======== raw data +TH1F** he; +TH1F** hxf; +TH1F** hxn; + +TH2F** hxfVxn; +TH2F** heVxs; + +TH1F* hMultiHit; +TH2F** heVx; // e vs (xf-xn)/e + +TH2F* heVID; +TH2F* hxfVID; +TH2F* hxnVID; + +//====== cal data +TH1F** heCal; +TH2F** hxfCalVxnCal; +TH2F** heVxsCal; // raw e vs xf +TH2F** heCalVxCal; // eCal vs xCal +TH2F** heCalVxCalG; // eCal vs xCal + +TH2F* heCalID; // e vs detID +TH2F* heCalVz; +TH2F* heCalVzGC; +TH2F** hecalVzRow; + +//====== Ex data +TH1F* hEx; +TH1F** hExi; +TH2F** hExVxCal; + +TH2F* hExThetaCM; + +TH1F* hExCut1; +TH1F* hExCut2; + +//======= Recoil +TH2F* hrdtID; +TH1F** hrdt; // single recoil +TH1F** hrdtg; + +TH2F** hrdt2D; +TH2F** hrdt2Dg; +TH2F** hrdt2Dsum; + +TH1F* hrdtRate1; +TH1F* hrdtRate2; + +//======= multi-Hit +TH2I *hmult; +TH1I *hmultEZ; +TH2I *hArrayRDTMatrix; +TH2I *hArrayRDTMatrixG; + +//======= ARRAY-RDT time diff +TH1I *htdiff; +TH1I *htdiffg; + +/*************************** + ***************************/ +double zRange[2] = {-1000, 0}; // zMin, zMax + +double Ex, thetaCM; +int padID = 0; +TLatex text; + +int numCol, numRow; +ULong64_t NumEntries = 0; +ULong64_t ProcessedEntries = 0; +Float_t Frac = 0.1; ///Progress bar +TStopwatch StpWatch; + +//======= Canvas +TCanvas *cCanvas; +TString canvasTitle; +int lastRunID; +bool contFlag; +double runTime=0; + +//======= Recoil Cut +TCutG* cutG; //! //general temeprary pointer to cut + +TObjArray * cutList1; +TObjArray * cutList2; + +//======= Other Cuts +TCutG* EZCut; +Bool_t isEZCutFileOpen; + +#include "Monitors.h" + +//^########################################################### +//^ * Begin +//^########################################################### +void Monitors::Begin(TTree *tree){ + + TString option = GetOption(); + + NumEntries = tree->GetEntries(); + + canvasTitle = GetCanvasTitle(); + lastRunID = -1; + contFlag = false; + + printf("###########################################################\n"); + printf("########## SOLARIS Monitors.C #########\n"); + printf("###########################################################\n"); + + //===================================================== loading parameter + printf("################## loading parameter files\n"); + + AnalysisLib::LoadDetGeoAndReactionConfigFile(); + AnalysisLib::LoadXNCorr(); + AnalysisLib::LoadXFXN2ECorr(); + AnalysisLib::LoadXScaleCorr(); + AnalysisLib::LoadECorr(); + AnalysisLib::LoadRDTCorr(); + AnalysisLib::LoadReactionParas(true); + + if( (int) AnalysisLib::xnCorr.size() < mapping::NARRAY ) { isXNCorrOK = false; printf("!!!!!!!! size of xnCorr < NARRAY .\n"); } + if( (int) AnalysisLib::xfxneCorr.size() < mapping::NARRAY ) { isXFXNCorrOK = false; printf("!!!!!!!! size of xfxneCorr < NARRAY .\n"); } + if( (int) AnalysisLib::eCorr.size() < mapping::NARRAY ) { isXScaleCorrOK = false; printf("!!!!!!!! size of eCorr < NARRAY .\n"); } + if( (int) AnalysisLib::xScale.size() < mapping::NARRAY ) { isECorrOK = false; printf("!!!!!!!! size of xScale < NARRAY .\n"); } + if( (int) AnalysisLib::rdtCorr.size() < mapping::NRDT ) { isRDTCorrOK = false; printf("!!!!!!!! size of rdtCorr < NRDT .\n"); } + + numRow = AnalysisLib::detGeo.nDet; + numCol = mapping::NARRAY/numRow; + + //================ Get Recoil cuts; + cutG = new TCutG(); + + cutList1 = AnalysisLib::LoadListOfTCut(rdtCutFile1, "cutList"); + cutList2 = AnalysisLib::LoadListOfTCut(rdtCutFile2, "cutList"); + + //================ Get EZ cuts; + EZCut = AnalysisLib::LoadSingleTCut(ezCutFile); + + + //========================= Generate all of the histograms needed for drawing later on + printf("======================================== Histograms declaration\n"); + + gROOT->cd(); + + he = new TH1F * [mapping::NARRAY]; + hxf = new TH1F * [mapping::NARRAY]; + hxn = new TH1F * [mapping::NARRAY]; + hxfVxn = new TH2F * [mapping::NARRAY]; + heVxs = new TH2F * [mapping::NARRAY]; + heVx = new TH2F * [mapping::NARRAY]; + + heCal = new TH1F * [mapping::NARRAY]; + hxfCalVxnCal = new TH2F * [mapping::NARRAY]; + heVxsCal = new TH2F * [mapping::NARRAY]; + heCalVxCal = new TH2F * [mapping::NARRAY]; + heCalVxCalG = new TH2F * [mapping::NARRAY]; + + for (Int_t i = 0; i < mapping::NARRAY; i++) {//array loop + + he[i] = new TH1F(Form("he%d", i), Form("Raw e (ch=%d); e (channel); count", i), 2000, rawEnergyRange[0], rawEnergyRange[1]); + hxf[i] = new TH1F(Form("hxf%d", i), Form("Raw xf (ch=%d); xf (channel); count", i), 200, rawEnergyRange[0], rawEnergyRange[1]); + hxn[i] = new TH1F(Form("hxn%d", i), Form("Raw xn (ch=%d); xn (channel); count", i), 200, rawEnergyRange[0], rawEnergyRange[1]); + hxfVxn[i] = new TH2F(Form("hxfVxn%d",i), Form("Raw xf vs. xn (ch=%d);xf (channel);xn (channel)",i), 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]); + heVxs[i] = new TH2F(Form("heVxs%d", i), Form("Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", i), 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); + heVx[i] = new TH2F(Form("heVx%d",i), Form("Raw PSD E vs. X (ch=%d);X (channel);E (channel)",i), 500, -0.1, 1.1, 500, rawEnergyRange[0], rawEnergyRange[1]); + + heCal[i] = new TH1F(Form("heCal%d", i), Form("Corrected e (ch=%d); e (MeV); count", i), 2000, energyRange[0], energyRange[1]); + hxfCalVxnCal[i] = new TH2F(Form("hxfCalVxnCal%d",i), Form("Corrected XF vs. XN (ch=%d);XF (channel);XN (channel)",i), 500, 0, rawEnergyRange[1], 500, 0, rawEnergyRange[1]); + heVxsCal[i] = new TH2F(Form("heVxsCal%d", i), Form("Raw e vs Corrected xf+xn (ch=%d); corrected xf+xn (channel); Raw e (channel)", i), 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); + heCalVxCal[i] = new TH2F(Form("heCalVxCal%d",i), Form("Cal PSD E vs. X (ch=%d);X (cm);E (MeV)",i), 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]); + heCalVxCalG[i] = new TH2F(Form("heCalVxCalG%d",i), Form("Cal PSD E vs. X (ch=%d);X (cm);E (MeV)",i), 500, -2.5, AnalysisLib::detGeo.detLength + 2.5, 500, energyRange[0], energyRange[1]); + } + + heVID = new TH2F("heVID", "Raw e vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); + hxfVID = new TH2F("hxfVID", "Raw xf vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); + hxnVID = new TH2F("hxnVID", "Raw xn vs channel", mapping::NARRAY, 0, mapping::NARRAY, 500, rawEnergyRange[0], rawEnergyRange[1]); + + heCalID = new TH2F("heCalID", "Corrected E vs detID; detID; E / 10 keV", mapping::NARRAY, 0, mapping::NARRAY, 2000, energyRange[0], energyRange[1]); + + hMultiHit = new TH1F("hMultiHit", "Multi-Hit of Energy", 10, 0, 1); + + //====================== E-Z plot + heCalVz = new TH2F("heCalVz", "E vs. Z;Z (mm);E (MeV)" , 400, zRange[0], zRange[1], 400, energyRange[0], energyRange[1]); + heCalVzGC = new TH2F("heCalVzGC","E vs. Z gated;Z (mm);E (MeV)", 400, zRange[0], zRange[1], 400, 0, energyRange[1]); + + hecalVzRow = new TH2F * [numRow]; + for( int i = 0; i < numRow; i++){ + hecalVzRow[i] = new TH2F(Form("heCalVzRow%d", i), Form("E vs. Z (ch=%d-%d); Z (cm); E (MeV)", numCol*i, numCol*(i+1)-1), 500, zRange[0], zRange[1], 500, energyRange[0], energyRange[1]); + } + + //===================== energy spectrum + hEx = new TH1F("hEx",Form("excitation spectrum w/ goodFlag; Ex [MeV] ; Count / %4.0f keV", exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + + hExCut1 = new TH1F("hExCut1",Form("excitation spectrum w/ goodFlag; Ex [MeV] ; Count / %4.0f keV", exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + hExCut2 = new TH1F("hExCut2",Form("excitation spectrum w/ goodFlag; Ex [MeV] ; Count / %4.0f keV", exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + hExCut1->SetLineColor(2); + hExCut2->SetLineColor(4); + + hExi = new TH1F * [mapping::NARRAY]; + hExVxCal = new TH2F * [mapping::NARRAY]; + for(int i = 0 ; i < mapping::NARRAY; i++){ + hExi[i] = new TH1F(Form("hExi%02d", i), Form("Ex (det=%i) w/goodFlag; Ex [MeV]; Count / %4.0f keV",i, exRange[0]), (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + hExVxCal[i] = new TH2F(Form("hExVxCal%d",i), Form("Ex vs X (ch=%d); X (cm); Ex (MeV)", i), 500, -0.1, 1.1, (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + } + + hExThetaCM = new TH2F("hExThetaCM", "Ex vs ThetaCM; ThetaCM [deg]; Ex [MeV]", 200, thetaCMRange[0], thetaCMRange[1], (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + + //===================== Recoils + hrdtID = new TH2F("hrdtID", "RDT vs ID; ID; energy [ch]", 8, 0, 8, 500, TMath::Min(rdtERange[0], rdtDERange[0]), TMath::Max(rdtERange[1], rdtDERange[1])); + + hrdt = new TH1F * [mapping::NRDT]; + hrdtg = new TH1F * [mapping::NRDT]; + + hrdt2D = new TH2F * [mapping::NRDT/2]; + hrdt2Dg = new TH2F * [mapping::NRDT/2]; + hrdt2Dsum = new TH2F * [mapping::NRDT/2]; + + for (Int_t i = 0; i < mapping::NRDT ; i++) { + if( i % 2 == 0 ) hrdt[i] = new TH1F(Form("hrdt%d",i), Form("Raw Recoil E(ch=%d); E (channel)",i), 500, rdtERange[0], rdtERange[1]); + if( i % 2 == 0 ) hrdtg[i] = new TH1F(Form("hrdt%dg",i),Form("Raw Recoil E(ch=%d) gated; E (channel)",i), 500, rdtERange[0], rdtERange[1]); + if( i % 2 == 1 ) hrdt[i] = new TH1F(Form("hrdt%d",i), Form("Raw Recoil DE(ch=%d); DE (channel)",i), 500, rdtDERange[0], rdtDERange[1]); + if( i % 2 == 1 ) hrdtg[i] = new TH1F(Form("hrdt%dg",i),Form("Raw Recoil DE(ch=%d) gated; DE (channel)",i), 500, rdtDERange[0], rdtDERange[1]); + + ///dE vs E + if( i % 2 == 0 ) { + int tempID = i / 2; + hrdt2D[tempID] = new TH2F(Form("hrdt2D%d",tempID), Form("Raw Recoil DE vs Eres (dE=%d, E=%d); Eres (channel); DE (channel)", i+1, i), 500, rdtERange[0], rdtERange[1],500,rdtDERange[0],rdtDERange[1]); + hrdt2Dg[tempID] = new TH2F(Form("hrdt2Dg%d",tempID), Form("Gated Raw Recoil DE vs Eres (dE=%d, E=%d); Eres (channel); DE (channel)",i+1, i), 500, rdtERange[0], rdtERange[1],500,rdtDERange[0], rdtDERange[1]); + hrdt2Dsum[tempID] = new TH2F(Form("hrdt2Dsum%d",tempID), Form("Raw Recoil DE vs Eres+DE (dE=%d, E=%d); Eres+DE (channel); DE (channel)", i+1, i), 500, rdtERange[0], rdtERange[1]+rdtDERange[1], 500, rdtDERange[0], rdtDERange[1]); + } + } + + hrdtRate1 = new TH1F("hrdtRate1", "recoil rate 1 / min; min; count / 1 min", timeRange[1] - timeRange[0], timeRange[0], timeRange[1]); + hrdtRate2 = new TH1F("hrdtRate2", "recoil rate 2 / min; min; count / 1 min", timeRange[1] - timeRange[0], timeRange[0], timeRange[1]); + hrdtRate1->SetLineColor(2); + hrdtRate2->SetLineColor(4); + + //===================== multiplicity + hmultEZ = new TH1I("hmultEZ", "Filled EZ with coinTime and recoil", 10, 0, 10); + hmult = new TH2I("hmult", "Array Multiplicity vs Recoil Multiplicity; Array ; Recoil",10, 0, 10, 10, 0, 10); + hArrayRDTMatrix = new TH2I("hArrayRDTMatrix", "Array ID vs Recoil ID; Array ID; Recoil ID", 30, 0, 30, 8, 0, 8); + hArrayRDTMatrixG = new TH2I("hArrayRDTMatrixG", "Array ID vs Recoil ID / g; Array ID; Recoil ID", 30, 0, 30, 8, 0, 8); + + //===================== coincident time + htdiff = new TH1I("htdiff" ,"Coincident time (recoil-dE - array); time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]); + htdiffg = new TH1I("htdiffg","Coincident time (recoil-dE - array) w/ recoil gated; time [ch = 10ns]; count", coinTimeRange[1] - coinTimeRange[0], coinTimeRange[0], coinTimeRange[1]); + + + printf("======================================== End of histograms Declaration\n"); + StpWatch.Start(); + +} + +//^########################################################### +//^ * Process +//^########################################################### +Bool_t Monitors::Process(Long64_t entry){ + + if( entry == 0 ) printf("========== %s \n", __func__); + + if( ProcessedEntries > maxNumberEvent ) return kTRUE; + ProcessedEntries++; + + //@*********** Progress Bar ******************************************/ + 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; + } + + //@********** Get Branch *********************************************/ + b_Energy->GetEntry(entry); + b_XF->GetEntry(entry); + b_XN->GetEntry(entry); + b_EnergyTimestamp->GetEntry(entry); + + if( isRDTExist ){ + b_RDT->GetEntry(entry); + b_RDTTimestamp->GetEntry(entry); + } + + // if( isArrayTraceExist ) { + // ///b_Trace_Energy->GetEntry(entry); + // b_Trace_Energy_RiseTime->GetEntry(entry); + // b_Trace_Energy_Time->GetEntry(entry); + // } + + // if( isRDTTraceExist ){ + // ///b_Trace_RDT->GetEntry(entry); + // b_Trace_RDT_Time->GetEntry(entry); + // b_Trace_RDT_RiseTime->GetEntry(entry); + // } + + //@*********** initization ******************************************/ + for( int i = 0 ; i < mapping::NARRAY; i++){ + z[i] = TMath::QuietNaN(); + x[i] = TMath::QuietNaN(); + xCal[i] = TMath::QuietNaN(); + eCal[i] = TMath::QuietNaN(); + } + + //@*********** Apply Recoil correction here *************************/ + if( isRDTCorrOK ){ + for( int i = 0 ; i < mapping::NRDT; i++){ + rdt[i] = rdt[i]*AnalysisLib::rdtCorr[i][0] + AnalysisLib::rdtCorr[i][1]; + } + } + + //@*********** Array ************************************************/ + //Do calculations and fill histograms + Int_t recoilMulti = 0; + Int_t arrayMulti = 0; + Int_t multiEZ = 0; + bool rdtgate1 = false; + bool rdtgate2 = false; + bool coinFlag = false; + bool ezGate = false; + bool isGoodEventFlag = false; + + for (Int_t detID = 0; detID < mapping::NARRAY; detID++) { + + //@================== Filling raw data + he[detID]->Fill(e[detID]); + hxf[detID]->Fill(xf[detID]); + hxn[detID]->Fill(xn[detID]); + hxfVxn[detID]->Fill(xf[detID],xn[detID]); + heVxs[detID]->Fill(xf[detID]+xn[detID], e[detID]); + + heVID->Fill(detID, e[detID]); + hxfVID->Fill(detID, xf[detID]); + hxnVID->Fill(detID, xn[detID]); + + //if( !TMath::IsNaN(e[detID]) ) printf("%llu | %d | %f %f %f \n", entry, detID, e[detID], xf[detID], xn[detID]); + + //@==================== Basic gate + if( TMath::IsNaN(e[detID]) ) continue ; + ///if( ring[detID] < -100 || ring[detID] > 100 ) continue; + ///if( ring[detID] > 300 ) continue; + if( TMath::IsNaN(xn[detID]) && TMath::IsNaN(xf[detID]) ) continue ; + + //@==================== Skip detector + bool skipFlag = false; + for( unsigned int kk = 0; kk < skipDetID.size() ; kk++){ + if( detID == skipDetID[kk] ) { + skipFlag = true; + break; + } + } + if (skipFlag ) continue; + + //@==================== Basic gate + if( TMath::IsNaN(e[detID]) ) continue ; + ///if( ring[detID] < -100 || ring[detID] > 100 ) continue; + ///if( ring[detID] > 300 ) continue; + if( TMath::IsNaN(xn[detID]) && TMath::IsNaN(xf[detID]) ) continue ; + + //@==================== Calibrations go here + if( isXNCorrOK && isXFXNCorrOK ) xnCal[detID] = xn[detID] * AnalysisLib::xnCorr[detID] * AnalysisLib::xfxneCorr[detID][1] + AnalysisLib::xfxneCorr[detID][0]; + if( isXNCorrOK && isXFXNCorrOK ) xfCal[detID] = xf[detID] * AnalysisLib::xfxneCorr[detID][1] + AnalysisLib::xfxneCorr[detID][0]; + if( isECorrOK ) eCal[detID] = e[detID] / AnalysisLib::eCorr[detID][0] + AnalysisLib::eCorr[detID][1]; + + if( eCal[detID] < eCalCut[0] ) continue; + if( eCal[detID] > eCalCut[1] ) continue; + + //@===================== fill Calibrated data + heCal[detID]->Fill(eCal[detID]); + heCalID->Fill(detID, eCal[detID]); + hxfCalVxnCal[detID]->Fill(xfCal[detID], xnCal[detID]); + heVxsCal[detID]->Fill(xnCal[detID] + xfCal[detID], e[detID]); + + //@===================== calculate X + if( (xf[detID] > 0 || !TMath::IsNaN(xf[detID])) && ( xn[detID]>0 || !TMath::IsNaN(xn[detID])) ) { + ///x[detID] = 0.5*((xf[detID]-xn[detID]) / (xf[detID]+xn[detID]))+0.5; + x[detID] = 0.5*((xf[detID]-xn[detID]) / e[detID])+0.5; + } + + /// range of x is (0, 1) + if ( !TMath::IsNaN(xf[detID]) && !TMath::IsNaN(xn[detID]) ) xCal[detID] = 0.5 + 0.5 * (xfCal[detID] - xnCal[detID] ) / e[detID]; + if ( !TMath::IsNaN(xf[detID]) && TMath::IsNaN(xn[detID]) ) xCal[detID] = xfCal[detID]/ e[detID]; + if ( TMath::IsNaN(xf[detID]) && !TMath::IsNaN(xn[detID]) ) xCal[detID] = 1.0 - xnCal[detID]/ e[detID]; + + //@======= Scale xcal from (0,1) + if( isXScaleCorrOK ) xCal[detID] = (xCal[detID]-0.5)/AnalysisLib::xScale[detID] + 0.5; /// if include this scale, need to also inclused in Cali_littleTree + + if( abs(xCal[detID] - 0.5) > xGate/2. ) continue; + + //@==================== calculate Z + if( AnalysisLib::detGeo.firstPos > 0 ) { + z[detID] = AnalysisLib::detGeo.detLength*(1.0-xCal[detID]) + AnalysisLib::detGeo.detPos[detID%numCol]; + }else{ + z[detID] = AnalysisLib::detGeo.detLength*(xCal[detID]-1.0) + AnalysisLib::detGeo.detPos[detID%numCol]; + } + + //@===================== multiplicity + arrayMulti++; /// multi-hit when both e, xf, xn are not NaN + + //@=================== Array fill + heVx[detID]->Fill(x[detID],e[detID]); + + heCalVxCal[detID]->Fill(xCal[detID]*AnalysisLib::detGeo.detLength,eCal[detID]); + heCalVz->Fill(z[detID],eCal[detID]); + + //@=================== Recoil Gate + if( isRDTExist && (cutList1 || cutList2)){ + for(int i = 0 ; i < cutList1->GetEntries() ; i++ ){ + cutG = (TCutG *)cutList1->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])) { + 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[detID]) ) { + for( int j = 0; j < mapping::NRDT ; j++){ + if( TMath::IsNaN(rdt[j]) ) continue; + + int tdiff = rdt_t[j] - e_t[detID]; + + if( j%2 == 1) { + htdiff->Fill(tdiff); + if((rdtgate1 || rdtgate2) && (eCalCut[1] > eCal[detID] && eCal[detID]>eCalCut[0])) { + htdiffg->Fill(tdiff); + } + } + + hArrayRDTMatrix->Fill(detID, 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(detID, j); + ///if( rdtgate1) hArrayRDTMatrixG->Fill(detID, j); + + hrdtg[j]->Fill(rdt[j]); + coinFlag = true; + + } + } + } + + if( !isTimeGateOn ) coinFlag = true; + + //================ E-Z gate + if( isEZCutFileOpen ) { + if( EZCut->IsInside(z[detID], eCal[detID]) ) ezGate = true; + }else{ + ezGate = true; + } + + if( coinFlag && (rdtgate1 || rdtgate2) && ezGate){ + heCalVzGC->Fill( z[detID] , eCal[detID] ); + + heCalVxCalG[detID]->Fill(xCal[detID]*AnalysisLib::detGeo.detLength,eCal[detID]); + + multiEZ ++; + isGoodEventFlag = true; + } + + }//end of array loop + + if( !isEZCutFileOpen ) ezGate = true; + + //@*********** RECOILS ***********************************************/ + for( int i = 0; i < mapping::NRDT ; i++){ + hrdtID->Fill(i, rdt[i]); + hrdt[i]->Fill(rdt[i]); + + if( i % 2 == 0 ){ + + recoilMulti++; // when both dE and E are hit + hrdt2D[i/2]->Fill(rdt[i],rdt[i+1]); //E-dE + } + } + + //@******************* Multi-hit *************************************/ + hmultEZ->Fill(multiEZ); + hmult->Fill(recoilMulti,arrayMulti); + hMultiHit->Fill(arrayMulti); + + //@*********** Good event Gate ***************************************/ + if( !isGoodEventFlag ) return kTRUE; + + //@*********** Ex and thetaCM ****************************************/ + for(Int_t detID = 0; detID < mapping::NARRAY ; detID++){ + + if( TMath::IsNaN(e[detID]) ) continue ; + if( TMath::IsNaN(z[detID]) ) continue ; + if( eCal[detID] < eCalCut[0] ) continue ; + if( eCal[detID] > eCalCut[1] ) continue ; + + if( AnalysisLib::hasReactionPara ){ + + std::vector ExThetaCM = AnalysisLib::CalExTheta(eCal[detID], x[detID]); + + Ex = ExThetaCM[0]; + thetaCM = ExThetaCM[1]; + + // ///======== Ex calculation by Ryan + // double y = eCal[detID] + mass; // to give the KE + mass of proton; + // double Z = alpha * gamm * betRel * z[detID]; + // double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ; + + // if( TMath::Abs(Z) < H ) { + // ///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi) + // double tolerrence = 0.001; + // double phi = 0; ///initial phi = 0 -> ensure the solution has f'(phi) > 0 + // double nPhi = 0; /// new phi + + // int iter = 0; + // do{ + // phi = nPhi; + // nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2)); + // iter ++; + // if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break; + // }while( TMath::Abs(phi - nPhi ) > tolerrence); + // phi = nPhi; + + // /// check f'(phi) > 0 + // double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2); + // if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){ + // double K = H * TMath::Sin(phi); + // double x = TMath::ACos( mass / ( y * gamm - K)); + // double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame + // double EB = TMath::Sqrt(mass*mass + Et*Et - 2*Et*TMath::Sqrt(momt*momt + mass * mass)); + // Ex = EB - massB; + + // double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y; + // double hahaha2 = gamm* betRel * momt; + // thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg(); + + // }else{ + // Ex = TMath::QuietNaN(); + // thetaCM = TMath::QuietNaN(); + // } + // }else{ + // Ex = TMath::QuietNaN(); + // thetaCM = TMath::QuietNaN(); + // } + }else{ + Ex = TMath::QuietNaN(); + thetaCM = TMath::QuietNaN(); + } + + if( thetaCM > thetaCMGate ) { + + hEx->Fill(Ex); + + hExThetaCM->Fill(thetaCM, Ex); + + if( rdtgate1 ) { + hExCut1->Fill(Ex); + hExThetaCM->Fill(thetaCM, Ex); + } + if( rdtgate2 ) { + hExCut2->Fill(Ex); + hExThetaCM->Fill(thetaCM, Ex); + } + + hExi[detID]->Fill(Ex); + hExVxCal[detID]->Fill(xCal[detID], Ex); + + } + } + + return kTRUE; +} + +//^########################################################### +//^ * Terminate +//^########################################################### +void Monitors::Terminate(){ + printf("============================== finishing.\n"); + + gROOT->cd(); + + int strLen = canvasTitle.Sizeof(); + canvasTitle.Remove(strLen-3); + + TString runTimeStr = ""; + if( runTime > 0. ) { + runTimeStr = Form("%.0f min", runTime); + canvasTitle += " | " + runTimeStr; + } + + //############################################ User is free to edit this section + //--- Canvas Size + int canvasXY[2] = {1200 , 800} ;// x, y + int canvasDiv[2] = {3,2}; + cCanvas = new TCanvas("cCanvas",canvasTitle + " | " + rdtCutFile1,canvasXY[0],canvasXY[1]); + cCanvas->Modified(); cCanvas->Update(); + cCanvas->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]); + + gStyle->SetOptStat("neiou"); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(2); + + double yMax = 0; + + Isotope hRecoil(AnalysisLib::reactionConfig.recoilHeavyA, AnalysisLib::reactionConfig.recoilHeavyZ); + double Sn = hRecoil.CalSp(0,1); + double Sp = hRecoil.CalSp(1,0); + double Sa = hRecoil.CalSp2(4,2); + + //TODO, Module each plot + ///----------------------------------- Canvas - 1 + PlotEZ(1); /// raw EZ + + ///----------------------------------- Canvas - 2 + PlotEZ(0); ///gated EZ + + ///----------------------------------- Canvas - 3 + PlotTDiff(1, 1); ///with Gated Tdiff, isLog + + ///----------------------------------- Canvas - 4 + padID++; cCanvas->cd(padID); + + //hEx->Draw(); + hExCut1->Draw(""); + hExCut2->Draw("same"); + DrawLine(hEx, Sn); + DrawLine(hEx, Sa); + + if(isTimeGateOn)text.DrawLatex(0.15, 0.8, Form("%d < coinTime < %d", timeGate[0], timeGate[1])); + if( xGate < 1 ) text.DrawLatex(0.15, 0.75, Form("with |x-0.5|<%.4f", xGate/2.)); + if( cutList1 ) text.DrawLatex(0.15, 0.7, "with recoil gated"); + + ///----------------------------------- Canvas - 5 + padID++; cCanvas->cd(padID); + + //Draw2DHist(hExThetaCM); + //heVIDG->Draw(); + //text.DrawLatex(0.15, 0.75, Form("#theta_{cm} > %.1f deg", thetaCMGate)); + + Draw2DHist(hrdt2D[0]); +// Draw2DHist(hrdt2Dsum[0]); + + if( cutList1 && cutList1->GetEntries() > 0 ) {cutG = (TCutG *)cutList1->At(0) ; cutG->Draw("same");} + if( cutList2 && cutList2->GetEntries() > 0 ) {cutG = (TCutG *)cutList2->At(0) ; cutG->Draw("same");} + + + //helum4D->Draw(); + //text.DrawLatex(0.25, 0.3, Form("gated from 800 to 1200 ch\n")); + + ///----------------------------------- Canvas - 6 + PlotRDT(0,0); + +// padID++; cCanvas->cd(padID); +// Draw2DHist(hrdtExGated); + + //padID++; cCanvas->cd(padID); + //Draw2DHist(htacEx); + + ///------------------------------------- Canvas - 7 + //PlotRDT(0, 0); + + ///----------------------------------- Canvas - 8 + //PlotRDT(1, 0); + + ///yMax = hic2->GetMaximum()*1.2; + ///hic2->GetYaxis()->SetRangeUser(0, yMax); + ///hic2->Draw(); + ///TBox * box14N = new TBox (-10, 0, -2, yMax); + ///box14N->SetFillColorAlpha(2, 0.1); + ///box14N->Draw(); + /// + ///TBox * box14C = new TBox (8, 0, 16, yMax); + ///box14C->SetFillColorAlpha(4, 0.1); + ///box14C->Draw(); + /// + ///text.SetTextColor(2); text.DrawLatex(0.38, 0.50, "14N"); + ///text.SetTextColor(4); text.DrawLatex(0.6, 0.45, "14C"); + ///text.SetTextColor(2); + ///----------------------------------- Canvas - 9 + //padID++; cCanvas->cd(padID); + + //Draw2DHist(hic01); + + ///----------------------------------- Canvas - 10 + //PlotRDT(3,0); + + //TH1F * helumDBIC = new TH1F("helumDBIC", "elum(d)/BIC; time [min]; count/min", timeRange[1]-timeRange[0], timeRange[0], timeRange[1]); + //helumDBIC = (TH1F*) helum4D->Clone(); + //helumDBIC->SetTitle("elum(d)/BIC; time [min]; count/min"); + //helumDBIC->SetName("helumDBIC"); + //helumDBIC->SetLineColor(2); + + //helumDBIC->Divide(hBIC); + + //yMax = helumDBIC->GetMaximum(); + //if( yMax < hBIC->GetMaximum() ) yMax = hBIC->GetMaximum(); + + //helumDBIC->SetMaximum(yMax * 1.2); + //hBIC->SetMaximum(yMax * 1.2); + + //hBIC->Draw(); + //helumDBIC->Draw("same"); + + //text.DrawLatex(0.15, 0.5, Form("Elum(D) / BIC \n")); + + ///----------------------------------- Canvas - 11 + //PlotRDT(2,0); + + ///----------------------------------- Canvas - 12 + //padID++; cCanvas->cd(padID); + //htac->Draw(); + + + /* + ///----------------------------------- Canvas - 13 + padID++; cCanvas->cd(padID); + + ///hicT14N->Draw(""); + ///hicT14C->Draw("same"); + /// + ///text.SetTextColor(2); text.DrawLatex(0.15, 0.60, "14N"); + ///text.SetTextColor(4); text.DrawLatex(0.15, 0.25, "14C"); + ///text.SetTextColor(2); + + ///----------------------------------- Canvas - 14 + padID++; cCanvas->cd(padID); + + ///hrdtRate1->Draw(""); + ///hrdtRate2->Draw("same"); + + ///----------------------------------- Canvas - 15 + padID++; cCanvas->cd(padID); + + ///----------------------------------- Canvas - 16 + padID++; cCanvas->cd(padID); + + + + ///----------------------------------- Canvas - 17 + padID++; cCanvas->cd(padID); + + + ///----------------------------------- Canvas - 18 + padID++; cCanvas->cd(padID); + + + ///----------------------------------- Canvas - 19 + padID++; cCanvas->cd(padID); + + + + ///----------------------------------- Canvas - 20 + padID++; cCanvas->cd(padID); + + htac->Draw(); + */ + + /************************************/ + gStyle->GetAttDate()->SetTextSize(0.02); + gStyle->SetOptDate(1); + gStyle->SetDateX(0); + gStyle->SetDateY(0); + + /************************************/ + StpWatch.Start(kFALSE); + + gROOT->ProcessLine(".L ../Armory/Monitor_Util.C"); + printf("=============== loaded Monitor_Utils.C\n"); + gROOT->ProcessLine(".L ../Armory/AutoFit.C"); + printf("=============== loaded Armory/AutoFit.C\n"); + gROOT->ProcessLine(".L ../Armory/RDTCutCreator.C"); + printf("=============== loaded Armory/RDTCutCreator.C\n"); + gROOT->ProcessLine(".L ../Armory/Check_rdtGate.C"); + printf("=============== loaded Armory/Check_rdtGate.C\n"); + gROOT->ProcessLine(".L ../Armory/readTrace.C"); + printf("=============== loaded Armory/readTrace.C\n"); + //gROOT->ProcessLine(".L ../Armory/readRawTrace.C"); + //printf("=============== loaded Armory/readRawTrace.C\n"); + gROOT->ProcessLine("listDraws()"); + + /************************* Save histograms to root file*/ + + gROOT->cd(); + + + /************************************/ + //gROOT->ProcessLine("recoils()"); + +} diff --git a/working/Monitors.h b/working/Monitors.h new file mode 100644 index 0000000..d46cff0 --- /dev/null +++ b/working/Monitors.h @@ -0,0 +1,368 @@ +#ifndef Monitors_h +#define Monitors_h + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Mapping.h" +#include "../armory/AnalysisLib.h" + +class Monitors : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + + // Declaration of leaf types + Float_t * e; //! + ULong64_t * e_t; //! + Float_t * xf; //! + ULong64_t * xf_t; //! + Float_t * xn; //! + ULong64_t * xn_t; //! + Float_t * rdt; //! + ULong64_t * rdt_t; //! + + // List of branches + TBranch *b_Energy; //! + TBranch *b_EnergyTimestamp; //! + TBranch *b_XF; //! + TBranch *b_XFTimestamp; //! + TBranch *b_XN; //! + TBranch *b_XNTimestamp; //! + TBranch *b_RDT; //! + TBranch *b_RDTTimestamp; //! + + // trace analysis data + Float_t * we; //! + Float_t * weR; //! + Float_t * weT; //! + Float_t * wrdt; //! + Float_t * wrdtT; //! + Float_t * wrdtR; //! + + TBranch *b_Trace_Energy; //! + TBranch *b_Trace_Energy_RiseTime; //! + TBranch *b_Trace_Energy_Time; //! + TBranch *b_Trace_RDT; //! + TBranch *b_Trace_RDT_Time; //! + TBranch *b_Trace_RDT_RiseTime; //! + + bool isArrayTraceExist; + bool isRDTTraceExist; + + bool isRDTExist; + + bool isXNCorrOK; + bool isXFXNCorrOK; + bool isXScaleCorrOK; + bool isECorrOK; + bool isRDTCorrOK; + + //==== global variable + float * x, * z; + float * xCal, * xfCal, * xnCal, * eCal; + + ULong64_t startTime ; + ULong64_t endTime ; + + Monitors(TTree * /*tree*/ =0) : fChain(0) { + + e = new Float_t [mapping::NARRAY]; + xf = new Float_t [mapping::NARRAY]; + xn = new Float_t [mapping::NARRAY]; + rdt = new Float_t [mapping::NRDT]; + e_t = new ULong64_t [mapping::NARRAY]; + xf_t = new ULong64_t [mapping::NARRAY]; + xn_t = new ULong64_t [mapping::NARRAY]; + rdt_t = new ULong64_t [mapping::NRDT]; + + x = new float [mapping::NARRAY]; + z = new float [mapping::NARRAY]; + xCal = new float [mapping::NARRAY]; + xfCal = new float [mapping::NARRAY]; + xnCal = new float [mapping::NARRAY]; + eCal = new float [mapping::NARRAY]; + + + isXNCorrOK = true; + isXFXNCorrOK = true; + isXScaleCorrOK = true; + isECorrOK = true; + isRDTCorrOK = true; + + } + virtual ~Monitors() { + + delete e ; + delete xf ; + delete xn ; + delete rdt ; + delete e_t ; + delete xf_t ; + delete xn_t ; + delete rdt_t; + + delete z ; + delete x ; + delete xCal; + delete xfCal; + delete xnCal; + delete eCal; + + } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + TString fCanvasTitle; + void SetCanvasTitle(TString title) {fCanvasTitle = title;} + TString GetCanvasTitle() const {return fCanvasTitle;} + + int timeRange[2]; + void SetTimeRange0(int minute){ timeRange[0] = minute;} + void SetTimeRange1(int minute){ timeRange[1] = minute;} + + void Draw2DHist(TH2F * hist); + + void PlotEZ(bool isRaw); + void PlotTDiff(bool isGated, bool isLog); + void PlotRDT(int id, bool isRaw); + //void PlotCRDTPolar(); + + //template T ** CreateListOfHist(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax, int binY = 0, float yMin = 0, float yMax = 0); + + ClassDef(Monitors,0); +}; + +#endif + + +#ifdef Monitors_cxx +void Monitors::Init(TTree *tree){ + + printf("========== %s \n", __func__); + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("e", e, &b_Energy); + fChain->SetBranchAddress("e_t", e_t, &b_EnergyTimestamp); + fChain->SetBranchAddress("xf", xf, &b_XF); + fChain->SetBranchAddress("xf_t", xf_t, &b_XFTimestamp); + fChain->SetBranchAddress("xn", xn, &b_XN); + fChain->SetBranchAddress("xn_t", xn_t, &b_XNTimestamp); + + TBranch * br = (TBranch *) fChain->GetListOfBranches()->FindObject("rdt"); + if( br == NULL ){ + printf(" ++++++++ no Recoil.\n"); + isRDTExist = false; + }else{ + printf(" ++++++++ FOund Recoil.\n"); + isRDTExist = true; + fChain->SetBranchAddress("rdt" , rdt, &b_RDT); + fChain->SetBranchAddress("rdt_t", rdt_t, &b_RDTTimestamp); + } + + /* + br = (TBranch *) fChain->GetListOfBranches()->FindObject("we"); + if( br == NULL ){ + printf(" ++++++++ no Array trace.\n"); + isArrayTraceExist = false; + }else{ + isArrayTraceExist = true; + if( isUseArrayTrace ){ + fChain->SetBranchAddress("te", e, &b_Energy);// replace e with te + }else{ + fChain->SetBranchAddress("te", te, &b_Trace_Energy); + } + fChain->SetBranchAddress("te_r", te_r, &b_Trace_Energy_RiseTime); + fChain->SetBranchAddress("te_t", te_t, &b_Trace_Energy_Time); + } + + br = (TBranch *) fChain->GetListOfBranches()->FindObject("wrdt"); + if( br == NULL ){ + printf(" ++++++++ no Recoil trace.\n"); + isRDTTraceExist = false; + }else{ + isRDTTraceExist = true; + if( isUseRDTTrace ) { + fChain->SetBranchAddress("trdt", rdt, &b_RDT); // replace rdt with trdt + printf("************ using Trace in recoil \n"); + }else{ + fChain->SetBranchAddress("trdt", trdt, &b_Trace_RDT); + } + fChain->SetBranchAddress("trdt_t", trdt_t, &b_Trace_RDT_Time); + fChain->SetBranchAddress("trdt_r", trdt_r, &b_Trace_RDT_RiseTime); + } + */ + + + startTime = 0; + endTime = 0; + + printf("=================================== End of Branch Pointer Inititization. \n"); +} + +Bool_t Monitors::Notify(){ + return kTRUE; +} + +void DrawLine(TH1 * hist, double pos){ + + double yMax = hist->GetMaximum(); + TLine * line = new TLine(pos, 0, pos, yMax); + line->SetLineColor(2); + line->Draw(""); + +} + + +void Monitors::SlaveBegin(TTree * /*tree*/){ + /// not use, if use, place in Monitors.C + TString option = GetOption(); +} + + +void Monitors::SlaveTerminate(){ + /// not use, if use, place in Monitors.C +} + + + +/*########################################################### + * Plotting Function +###########################################################*/ + +void DrawBox(TH1* hist, double x1, double x2, Color_t color, float alpha){ + + double yMax = hist->GetMaximum(); + TBox * box = new TBox (x1, 0, x2, yMax); + box->SetFillColorAlpha(color, alpha); + box->Draw(); + +} + +void Monitors::Draw2DHist(TH2F * hist){ + + if( hist->Integral() < 3000 ){ + hist->SetMarkerStyle(20); + hist->SetMarkerSize(0.3); + hist->Draw(""); + }else{ + hist->Draw("colz"); + } +} + + +void Monitors::PlotEZ(bool isRaw){ + padID++; cCanvas->cd(padID); + + if( isRaw ) { + Draw2DHist(heCalVz); + heCalVz->SetTitle("E vs Z | " + canvasTitle + " | " + rdtCutFile1); + if( skipDetID.size() > 0 ){ + text.DrawLatex(0.15, 0.3, "skipped Detector:"); + for( int i = 0; i < (int) skipDetID.size(); i++){ + text.DrawLatex(0.15 + 0.1*i, 0.25, Form("%d", skipDetID[i])); + } + } + + text.DrawLatex(0.15, 0.8, Form("%.1f < eCal < %.1f MeV", eCalCut[0], eCalCut[1])); + if( xGate < 1 ) text.DrawLatex(0.15, 0.75, Form("with |x-0.5|<%.4f", xGate/2.)); + + }else{ + Draw2DHist(heCalVzGC); + + if( cutList1 ) text.DrawLatex(0.15, 0.8, "with Recoil gate"); + if( xGate < 1 ) text.DrawLatex(0.15, 0.75, Form("with |x-0.5|<%.4f", xGate/2.)); + //if( isTACGate ) text.DrawLatex(0.15, 0.7, Form("%d < TAC < %d", tacGate[0], tacGate[1])); + if(isTimeGateOn)text.DrawLatex(0.15, 0.7, Form("%d < coinTime < %d", timeGate[0], timeGate[1])); + + } + + TFile * transfer = new TFile("transfer.root"); + TObjArray * gList = NULL ; + TObjArray * fxList = NULL ; + int nGList = 0; + int nFxList = 0; + if( transfer->IsOpen() ) { + gList = (TObjArray *) transfer->FindObjectAny("gList"); + nGList = gList->GetLast() + 1; + fxList = (TObjArray *) transfer->FindObjectAny("fxList"); + nFxList = fxList->GetLast() +1 ; + } + + ///the constant thetaCM line + if( transfer->IsOpen() ) gList->At(0)->Draw("same"); + ///the e-z line for excitation + if( transfer->IsOpen() ){ + for( int i = 0 ; i < nFxList ; i++){ + ((TF1*)fxList->At(i))->SetLineColor(6); + fxList->At(i)->Draw("same"); + } + } + +} + +void Monitors::PlotTDiff(bool isGated, bool isLog){ + padID++; cCanvas->cd(padID); + if( isLog ) cCanvas->cd(padID)->SetLogy(1); + double yMax = 0; + if( isGated ){ + yMax = htdiff->GetMaximum()*1.2; + if( isLog ){ + htdiff->GetYaxis()->SetRangeUser(1, yMax); + }else{ + htdiff->GetYaxis()->SetRangeUser(0, yMax); + } + } + htdiff->Draw(); + if( isGated ){ + htdiffg->SetLineColor(2); + htdiffg->Draw("same"); + } + + if( cutList1 ) text.DrawLatex(0.15, 0.8, "with Recoil gate"); + if(isTimeGateOn)text.DrawLatex(0.15, 0.7, Form("%d < coinTime < %d", timeGate[0], timeGate[1])); + DrawBox(htdiff, timeGate[0], timeGate[1], kGreen, 0.2); +} + +void Monitors::PlotRDT(int id, bool isRaw){ + padID++; cCanvas->cd(padID); + + if( isRaw ){ + Draw2DHist(hrdt2D[id]); + }else{ + Draw2DHist(hrdt2Dg[id]); + } + if(isTimeGateOn)text.DrawLatex(0.15, 0.8, Form("%d < coinTime < %d", timeGate[0], timeGate[1])); + //if( isTACGate ) text.DrawLatex(0.15, 0.7, Form("%d < TAC < %d", tacGate[0], tacGate[1])); + if( cutList1 && cutList1->GetEntries() > id ) {cutG = (TCutG *)cutList1->At(id) ; cutG->Draw("same");} + if( cutList2 && cutList2->GetEntries() > id ) {cutG = (TCutG *)cutList2->At(id) ; cutG->Draw("same");} + +} + +//void Monitors::PlotCRDTPolar(){ +// padID++; cCanvas->cd(padID); +// cCanvas->cd(padID)->DrawFrame(-50, -50, 50, 50); +// hcrdtPolar->Draw("same colz pol"); +//} + +#endif // #ifdef Monitors_cxx