diff --git a/Working/ChainMonitors.C b/Working/ChainMonitors.C new file mode 100644 index 0000000..83255a9 --- /dev/null +++ b/Working/ChainMonitors.C @@ -0,0 +1,37 @@ +#include "Monitor.C+" // the plus sign mean compilation +#include "TObjArray.h" +#include "TFile.h" +#include "TMacro.h" +#include "TChain.h" + +TChain *gen_tree = nullptr; + +void ChainMonitors(int RUNNUM = -1, int RUNNUM2 = -1) { + ///default saveCanvas = false, no save Cavas + /// = true, save Canvas + + gen_tree = new TChain("gen_tree"); + if( RUNNUM == -1){ + /// this list only for manual Chain sort + ///********** start Marker for AutoCalibration. + + gen_tree->Add("../root_data/trace_run033.root"); + + ///********** end Marker for AutoCalibration. + + }else{ + + TString fileName; + int endRUNNUM = RUNNUM2; + if( RUNNUM2 == -1) endRUNNUM = RUNNUM; + + for( int i = RUNNUM ; i <= endRUNNUM ; i++){ + fileName.Form("../root_data/gen_run%03d.root", i); + gen_tree->Add(fileName); + } + } + + //^============== should have other things, like calibrations. + Monitor(gen_tree); + +} diff --git a/Working/ClassMonPlotter.h b/Working/ClassMonPlotter.h new file mode 100644 index 0000000..9f40f4d --- /dev/null +++ b/Working/ClassMonPlotter.h @@ -0,0 +1,570 @@ +#ifndef ClassMonitorPlotter_H +#define ClassMonitorPlotter_H + +#include "../Armory/ClassDetGeo.h" +#include "../Armory/ClassReactionConfig.h" +#include "../Cleopatra/ClassTransfer.h" +#include "../Cleopatra/ClassIsotope.h" + +#include "TH1.h" +#include "TH2.h" +#include "TCanvas.h" +#include "TLine.h" +#include "TStyle.h" + +/****************************************************************** +* This is Plotter for Monitor.C. It contains +* 1) Tcanvas +* 2) various Histograms ( exclude raw data histogram ) +* +* The reason for having Plotter is suppert multiple arrays. +* contained the Canvas and Histogram in a class, have better memory management +* +*******************************************************************/ +/****************************************************************** +* variable and histogram naming rules * +* name are case sensitive, so as any C/C++ code * +* * +* ID is detector 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 "_" to seperate 2 variables, like e_x * +* * +* histogram with TCutG, add suffix "GC" for Graphical-Cut. * +* * +*******************************************************************/ + +class MonPlotter{ +public: + + MonPlotter(unsigned short arrayID, DetGeo * detGeo, int numRDT); + ~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 Plot(); + + void PlotRaw(bool isLog = false); + void PlotCal(); + void PlotEZ(); + void PlotEx(); + + void PlotRDT(bool isLog = false); + + TCanvas * canvas; + + //====================== Histograms + //======== raw data + TH2F * he_ID, * hxf_ID, * hxn_ID; // vs ID + + TH1I * hArrayMulti; + + TH1F ** he, ** hxf, ** hxn; //basic data + TH2F ** hxf_xn, ** he_xs; // correlation + + //====== cal data + TH1F ** heCal; + 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; + TH2F * heCal_zGC; + + //======= Recoil + TH2F * hrdt_ID; + TH1F ** hrdt; // single recoil + + TH1I * hrdtMulti; + + TH2F ** hrdt2D; + TH2F ** hrdt2Dg; // gated + + //====== tDiff + TH1F * htDiff; + TH1F * htDiffg; + + //====== Ex data + TH1F * hEx; + TH1F ** hExi; + TH2F ** hEx_xCal; + + TH1F * hExCut1; + TH1F * hExCut2; + + TH2F * hEx_ThetaCM; + //======================= + + //======= Recoil Cut + TObjArray * cutList; + +private: + + unsigned short aID; + int numDet, colDet, rowDet; //array + float detLength; + int numRDT; + float recoilOutter; + double zRange[2] ; // zMin, zMax + + TString canvasTitle; + TString suffix; + int numPad; + + template void CreateListOfHist1D(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax); + template void CreateListOfHist2D(T ** &histList, int size, const char * namePrefix, const char * TitleForm, int binX, float xMin, float xMax, int binY, float yMin, float yMax); + +}; +//^################################################################################# +MonPlotter::MonPlotter(unsigned short arrayID, DetGeo * detGeo, int numRDT){ + aID = arrayID; + numDet = detGeo->array[aID].numDet; + colDet = detGeo->array[aID].colDet; + rowDet = numDet/colDet; + detLength = detGeo->array[aID].detLength; + + suffix = Form("_%d", arrayID); + + this->numRDT = numRDT; + recoilOutter = detGeo->aux[aID].outerRadius; + + zRange[0] = detGeo->array[aID].zMin - 50; + zRange[1] = detGeo->array[aID].zMax + 50; + + canvas = nullptr; + cutList = nullptr; + +} + +MonPlotter::~MonPlotter(){ + + printf("=============== %s\n", __func__); + + delete canvas; + + delete he_ID; + delete hxf_ID; + delete hxn_ID; + delete hArrayMulti; + + delete heCal_ID; + delete heCal_zGC; + delete heCal_z; + + delete hEx_ThetaCM; + delete hExCut1; + delete hExCut2; + + delete hrdt_ID; + delete hrdtMulti; + + delete htDiff; + delete htDiffg; + + for( int i = 0; i < numDet ; i++ ){ + delete he[i]; + delete hxf[i]; + delete hxn[i]; + delete hxf_xn[i]; + delete he_xs[i]; + delete he_x[i]; + delete heCal[i]; + delete hExi[i]; + delete hEx_xCal[i]; + } + + for( int i = 0; i < numRDT; i++){ + delete hrdt[i]; + } + for( int i = 0; i < numRDT/2; i++){ + delete hrdt2D[i]; + delete hrdt2Dg[i]; + } + + delete [] he; + delete [] hxf; + delete [] hxn; + delete [] hxf_xn; + delete [] he_xs; + delete [] he_x; + delete [] heCal; + delete [] hExi; + delete [] hEx_xCal; + + delete [] hrdt; + delete [] hrdt2D; + delete [] hrdt2Dg; + + 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->Divide(divX, divY); + numPad = divX * divY; + canvasTitle = title; + +} + +template void MonPlotter::CreateListOfHist1D(T ** &histList, + int size, + const char * namePrefix, + const char * TitleForm, + int binX, float xMin, float xMax){ + + //printf(" Making %d of %s.\n", size, namePrefix); + histList = new T * [size]; + for(int i = 0; i < size; i++) { + histList[i] = new T(Form("%s%d", namePrefix, i) + suffix, Form(TitleForm, i), binX, xMin, xMax); + } +} + +template void MonPlotter::CreateListOfHist2D(T ** &histList, + int size, + const char * namePrefix, + const char * TitleForm, + int binX, float xMin, float xMax, + int binY, float yMin, float yMax){ + + //printf(" Making %d of %s.\n", size, namePrefix); + histList = new T * [size]; + for(int i = 0; i < size; i++) { + histList[i] = new T(Form("%s%d", namePrefix, i) + suffix, Form(TitleForm, i), binX, xMin, xMax, binY, yMin, yMax); + } +} + +void MonPlotter::SetUpHistograms(int * rawEnergyRange, + int * energyRange, + double * exRange, + int * thetaCMRange, + int * rdtDERange, + int * rdtERange, + int * coinTimeRange){ + + he_ID = new TH2F("he_ID" + suffix, "Raw e vs array ID; Array ID; Raw e", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]); + hxf_ID = new TH2F("hxf_ID" + suffix, "Raw xf vs array ID; Array ID; Raw xf", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]); + hxn_ID = new TH2F("hxn_ID" + suffix, "Raw xn vs array ID; Array ID; Raw xn", numDet, 0, numDet, 200, rawEnergyRange[0], rawEnergyRange[1]); + + hArrayMulti = new TH1I("hArrayMulti", "Array Multiplicity ( e and (xf or xn) )", numDet, 0, numDet); + + CreateListOfHist1D(he, numDet, "he", "Raw e (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]); + CreateListOfHist1D(hxf, numDet, "hxf", "Raw xf (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]); + CreateListOfHist1D(hxn, numDet, "hxn", "Raw xn (ch=%d); e (channel); count", 200, rawEnergyRange[0], rawEnergyRange[1]); + + CreateListOfHist2D(hxf_xn, numDet, "hxf_xn", "Raw xf vs. xn (ch=%d);xf (channel);xn (channel)" , 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); + CreateListOfHist2D(he_xs, numDet, "he_xs", "Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); + + CreateListOfHist2D(he_x , numDet, "he_x", "Raw e vs x (ch=%d); x (mm); Raw e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, -0.5, 1.5); + 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]); + + //====================== 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]); + + //===================== Recoil + int rdtRange[2]; + rdtRange[0] = rdtDERange[0] < rdtERange[0] ? rdtDERange[0] : rdtERange[0]; + rdtRange[1] = rdtDERange[1] > rdtERange[1] ? rdtDERange[1] : rdtERange[1]; + + hrdt_ID = new TH2F("hrdt_ID" + suffix, "Raw RDT vs ID; ID; Raw RDT", numRDT, 0, numRDT, 400, rdtRange[0], rdtRange[1]); + + hrdtMulti = new TH1I("hrdtMulti" + suffix, "RDT Multiplicity", numRDT, 0, numRDT); + + hrdt = new TH1F * [numRDT]; + hrdt2D = new TH2F * [numRDT/2]; + hrdt2Dg = new TH2F * [numRDT/2]; + + for (Int_t i = 0; i < numRDT ; 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 == 1 ) hrdt[i] = new TH1F(Form("hrdt%d",i), Form("Raw Recoil DE(ch=%d); 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]); + } + } + + //===================== 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]); + + TString haha = "Ex (det=%i) w/goodFlag; Ex [MeV]; Count / " +std::to_string(exRange[0]) + "keV"; + + hExi = new TH1F * [numDet]; + hEx_xCal = new TH2F * [numDet]; + for(int i = 0; i < numDet; i++ ){ + hExi[i] = new TH1F(Form("hExi%d", i) + suffix, haha, (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + + hEx_xCal[i] = new TH2F(Form("hEx_xCal%d", i) + suffix, + 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]); + } + + hExCut1 = new TH1F("hExCut1" + 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]); + hExCut2 = new TH1F("hExCut2" + 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]); + hExCut1->SetLineColor(2); + hExCut2->SetLineColor(4); + + hEx_ThetaCM = new TH2F("hExThetaCM" + suffix, "Ex vs ThetaCM; ThetaCM [deg]; Ex [MeV]", 200, thetaCMRange[0], thetaCMRange[1], (int) (exRange[2]-exRange[1])/exRange[0]*1000, exRange[1], exRange[2]); + +} + +//^####################################################### Plot +void MonPlotter::Plot(){ + + //TODO a more user-friendly way. + //TODO display text on the 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(""); + htDiffg->Draw("same"); + }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" , rdtCutFile.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; icd(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; icd(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;icd(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;icd(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;icd(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 = new TLine(0,0, 4000, 4000); + line->SetLineColor(2); + for (Int_t i=0;icd(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; icd(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;icd(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 *cExVxCal = new TCanvas("cExVxCal" + suffix,Form("EX | %s", canvasTitle.Data()), 200 + 1000 * aID, 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(); + } + + TCanvas *cexI = new TCanvas("cexI" + suffix,Form("EX : %s", canvasTitle.Data()),300 + 1000 * aID, 300, 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()), 400 + 1000 * aID, 400, 650,650); + cExThetaCM->Clear(); + gStyle->SetOptStat("neiou"); + hEx_ThetaCM->Draw("colz"); + + TCanvas *cex = new TCanvas("cex" + suffix,Form("EX : %s", canvasTitle.Data()), 500 + 1000 * aID, 500, 1000,650); + cex->Clear(); + gStyle->SetOptStat("neiou"); + hEx->Draw(""); + +} + +void MonPlotter::PlotRDT(bool isLog){ + + TCanvas *crdt = new TCanvas("crdt" + suffix,Form("raw RDT | %s", canvasTitle.Data()), 1000, 0, 1000,1000); + crdt->Clear();crdt->Divide(numRDT/4,2); + + for( int i = 0; i < numRDT/2; i++){ + if( isLog ) crdt->cd(i+1)->SetLogz(); crdt->cd(i+1); hrdt2D[i]->Draw("col"); + } + + TCanvas *crdtID = new TCanvas("crdtID" + suffix,Form("raw RDT ID | %s", canvasTitle.Data()),1100,1100, 500, 500); + crdtID->Clear(); + if( isLog ) crdtID->SetLogz(); + hrdt_ID->Draw("colz"); + + TCanvas *crdtS = new TCanvas("crdtS" + suffix,Form("raw RDT | %s", canvasTitle.Data()),1200, 1200, 1000, 1000); + crdtS->Clear(); crdtS->Divide(2,numRDT/2); + for( int i = 0; i < numRDT; i ++){ + crdtS->cd(i+1); + if( isLog ) crdtS->cd(i+1)->SetLogy(); + hrdt[i]->Draw(""); + } + +} + +#endif \ No newline at end of file diff --git a/Working/DWBA b/Working/DWBA new file mode 100644 index 0000000..d3986eb --- /dev/null +++ b/Working/DWBA @@ -0,0 +1,50 @@ +#========= Input for Cleopatra +#===== # for comment line, must be at the beginning of line +# This is a simplified input for creating Ptolmey in-file. +# Ptolemy has many control parameters for specific situation. Use with caution or ask Ben Kay. +#===== the potential contain two words +# one for incoming +# one for outgoing +#================================================= Potenital abberation +#========================= deuteron +# A = An, Cai, 2006 | E < 183 | 12 < A < 238 | http://dx.doi.org/10.1103/PhysRevC.73.054605 +# H = Han, Shi, Shen, 2006 | E < 200 | 12 < A < 209 | http://dx.doi.org/10.1103/PhysRevC.74.044615 +# B = Bojowald et al., 1988 | 50 < E < 80 | 27 < A < 208 | http://dx.doi.org/10.1103/PhysRevC.38.1153 +# D = Daehnick, Childs, Vrcelj, 1980 | 11.8 < E < 80 | 27 < A < 238 (REL) | http://dx.doi.org/10.1103/PhysRevC.21.2253 +# C = Daehnick, Childs, Vrcelj, 1980 | 11.8 < E < 80 | 27 < A < 238 (NON-REL) | http://dx.doi.org/10.1103/PhysRevC.21.2253 // not impletmented yet +# L = Lohr and Haeberli, 1974 | 9 < E < 13 | 40 < A | http://dx.doi.org/10.1016/0375-9474(74)90627-7 +# Q = Perey and Perey, 1963 | 12 < E < 25 | 40 < A | http://dx.doi.org/10.1016/0370-1573(91)90039-O +# Z = Zhang, Pang, Lou, 2016 | 5 < E < 170 | A < 18, spe 6-7Li | https://doi.org/10.1103/PhysRevC.94.014619 +#========================= proton +# K = Koning and Delaroche, 2009 | 0.001 < E < 200 | 24 < A < 209 | Iso. Dep. | http://dx.doi.org/10.1016/S0375-9474(02)01321-0 +# V = Varner et al., (CH89), 1991 | 16 < E < 65 | 4 < A < 209 | http://dx.doi.org/10.1016/0370-1573(91)90039-O +# M = Menet et al., 1971 | 30 < E < 60 | 40 < A | http://dx/doi.org/10.1016/0092-640X(76)90007-3 +# G = Becchetti and Greenlees, 1969 | E < 50 | 40 < A | http://dx.doi.org/10.1103/PhysRev.182.1190 +# P = Perey, 1963 | E < 20 | 30 < A < 100 | http://dx/doi.org/10.1016/0092-640X(76)90007-3 +#========================= A=3 +# x = Xu, Guo, Han, Shen, 2011 | E < 250 | 20 < A < 209 | 3He | http://dx.doi.org/10.1007/s11433-011-4488-5 +# l = Liang, Li, Cai, 2009 | E < 270 | All masses | http://dx.doi.org/10.1088/0954-3899/36/8/085104 +# p = Pang et al., 2009 | All E | All masses | Iso. Dep. | http://dx.doi.org/10.1103/PhysRevC.79.024615 +# c = Li, Liang, Cai, 2007 | E < 40 | 48 < A < 232 | Tritons | http://dx.doi.org/10.1016/j.nuclphysa.2007.03.004 +# t = Trost et al., 1987 | 10 < E < 220 | 10 < A < 208 | http://dx.doi.org/10.1016/0375-9474(87)90551-3 +# h = Hyakutake et al., 1980 | 90 < E < 120 | About 58 < A < 92 | http://dx.doi.org/10.1016/0375-9474(80)90013-5 +# b = Becchetti and Greenlees, 1971 | E < 40 | 40 < A | Iso. Dep. +#========================= alpha +# s = Su and Han, 2015 | E < 398 | 20 < A < 209 | http://dx.doi/org/10.1142/S0218301315500925 +# a = Avrigeanu et al., 2009 | E ??? | A ??? | http://dx.doi/org/10.1016/j.adt.2009.02.001 +# f = Bassani and Picard, 1969(FIXED)| 24 < E < 31 | A = 90 | https://doi.org/10.1016/0375-9474(69)90601-0 +#======================================================================= +#reaction gs-spin orbital spin-pi(Ex) Ex ELab Potentials +#206Hg(d,d)206Hg 0 none 9/2+ 0.000 7.39MeV/u AA #elastic +#206Hg(d,d)206Hg 0 none 9/2+ 1.000 7.39MeV/u AA 0.12 #inelastics_0.12=beta +#206Hg(d,p)207Hg 0 1g9/2 9/2+ 0.000 7.39MeV/u AK +#20F(d,t)19F 2 0d5/2 5/2+ 0.197 10MeV/u Vl +#16N(d,3He)15C 2 0p1/2 5/2+ 0.74 12MeV/u Ax +#10Be(t,p)12Be 0 1L=0 0+ 0.000 5MeV/u lA #two-nucleon_transfer +#32Si(t,p)34Si 0 0L=0 0+ 0.000 8MeV/u lA #two-nucleon_transfer +#36Ar(d,a)34Cl 0 4L=2 3+ 0.000 8MeV/u As # (d,a) reaction + + +30Si(d,p)31Si 0 1s1/2 1/2+ 0.000 10MeV/u AK +# 32Si(d,p)33Si 0 0d5/2 5/2+ 0.197 10MeV/u AK +# 32Si(d,3He)31Al 0 0d5/2 5/2+ 0.000 10MeV/u Ax diff --git a/Working/Mapping.h b/Working/Mapping.h new file mode 100644 index 0000000..6500020 --- /dev/null +++ b/Working/Mapping.h @@ -0,0 +1,133 @@ +#ifndef MAPPING_H +#define MAPPING_H + +//^=============================================================================== +//^ This is mapping file for SOLARIS +//^ This file is used to constructe the SOLARIS panel in the SOLARIS DAQ +//^ If this file is modified, please Close Digitizer and Open again +//^------------------------------------------------------------------------------- +//^ +//^ Array-e : 0 - 99 +//^ Array-xf : 100 - 199 +//^ Array-xn : 200 - 299 +//^ Recoil : 300 - 399 +//^ Enum : 400 - 499 +//^ EZERO : 500 - 599 +//^ Apollo : 600 - 699 +//^ +//^ line comment is line constains '//^' or '// //' or '////' +//^ +//^=============================================================================== + +#include +#include + +namespace mapping{ + +const std::vector detTypeName = { "e", "xf", "xn", "rdt", "eNum"}; //C= The comment "//C=" is an indicator DON't Remove +const std::vector detGroupID = { 0, 0, 0, 1, 2}; //C& The comment "//C&" is an indicator DON't Remove +const std::vector detMaxID = { 100, 200, 300, 400, 500}; //C# The comment "//C#" is an indicator DON't Remove +const std::vector detParity = { 1, 1, 1, 1, 1}; + +const std::vector groupName = { "Array", "Recoil", "ELUM"}; //C% The comment "//C%" is an indicator DON't Remove + +//!The mapping[i] must match as the IP setting in the DAQ + +const std::vector> map = { +{ +//C 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 // this line is an indicator DON'T Remove "//C " is an indcator + 0, 100, 200, 1, 101, 201, 2, 102, 202, 3, 103, 203, 4, 104, 204, -1, /// 0 - 15 + 5, 105, 205, 6, 106, 206, 7, 107, 207, 8, 108, 208, 9, 109, 209, -1, /// 16 - 31 + 10, 110, 210, 11, 111, 211, 12, 112, 212, 13, 113, 213, 14, 114, 214, -1, /// 32 - 47 + 15, 115, 215, 16, 116, 216, 17, 117, 217, 18, 118, 218, 19, 119, 219, -1 /// 48 - 63 +//C------------- end of a digitizer // this line is an indicator DON'T Remove "//C-" is an indcator +} +,{ +//C 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 // this line is an indicator DON'T Remove "//C " is an indcator + 20, 120, 220, 21, 121, 221, 22, 122, 222, 23, 123, 223, 24, 124, 224, -1, /// 0 - 15 + 25, 125, 225, 26, 126, 226, 27, 127, 227, 28, 128, 228, 29, 129, 229, -1, /// 16 - 31 + 30, 130, 230, 31, 131, 231, 32, 132, 232, 33, 133, 233, 34, 134, 234, -1, /// 32 - 47 + 35, 135, 235, 36, 136, 236, 37, 137, 237, 38, 138, 238, 39, 139, 239, -1 /// 48 - 63 +//C------------- end of a digitizer // this line is an indicator DON'T Remove "//C-" is an indcator +} +,{ +//C 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 // this line is an indicator DON'T Remove "//C " is an indcator + 40, 140, 240, 41, 141, 241, 42, 142, 242, 43, 143, 243, 44, 144, 244, -1, /// 0 - 15 + 45, 145, 245, 46, 146, 246, 47, 147, 247, 48, 148, 248, 49, 149, 249, -1, /// 16 - 31 + 50, 150, 250, 51, 151, 251, 52, 152, 252, 53, 153, 253, 54, 154, 254, -1, /// 32 - 47 + 55, 155, 255, 56, 156, 256, 57, 157, 257, 58, 158, 258, 59, 159, 259, -1 /// 48 - 63 +//C------------- end of a digitizer // this line is an indicator DON'T Remove "//C-" is an indcator +} +,{ +//C 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 // this line is an indicator DON'T Remove "//C " is an indcator + 300, 301, 302, 303, 304, 305, 306, 307, -1, -1, -1, -1, -1, -1, -1, -1, /// 0 - 15 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, /// 16 - 31 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, /// 32 - 47 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 /// 48 - 63 +//C------------- end of a digitizer // this line is an indicator DON'T Remove "//C-" is an indcator +} +}; + +//^=============================================================================== +int FindDetTypeIndex(int detID){ + for( int k = 0; k < (int) detMaxID.size(); k++){ + int low = (k == 0 ? 0 : detMaxID[k-1]); + int high = detMaxID[k]; + if( low <= detID && detID < high ) { + return k; + } + } + return -1; +} + +std::vector ExtractDetNum(){ + std::vector detTypeNum; + for( int i = 0; i < (int) detTypeName.size(); i ++) detTypeNum.push_back(0); + for( int i = 0; i < (int) map.size(); i ++){ + for( int j = 0; j < (int) map[i].size(); j++){ + if( map[i][j] < 0) continue; + for( int k = 0; k < (int) detTypeName.size() ; k ++ ){ + int low = (k == 0 ? 0 : detMaxID[k-1]); + int high = detMaxID[k]; + if( low <= map[i][j] && map[i][j] < high ) { + detTypeNum[k]++; + } + } + } + } + return detTypeNum; +} + +void PrintMapping(){ + + //------------ Red Green Yellow Cyan blue Magenta Gray + std::vector Color = {"\033[31m", "\033[32m", "\033[33m", "\033[36m", "\033[34m", "\033[35m", "\033[37m"}; + + printf("==================================== Mapping ===================================\n"); + std::vector detTypeNum = ExtractDetNum(); + for(int i = 0 ; i < (int) detTypeName.size(); i++) { + printf(" %2d | %7s | %3d | %3d - %3d\n", i, detTypeName[i].c_str(), detTypeNum[i], (i == 0 ? 0 : detMaxID[i-1]), detMaxID[i]); + } + for( int i = 0; i < (int) map.size(); i++){ + printf("Digi-%d ------------------------------------------------------------------------ \n", i); + for( int j = 0; j < (int) map[i].size(); j++){ + if( map[i][j] < 0 ){ + printf("%4d,", map[i][j]); + }else{ + int colorIndex = FindDetTypeIndex(map[i][j]); + printf("%s%4d\033[0m,", Color[colorIndex], map[i][j]); + } + if( j % 16 == 15 ) printf("\n"); + } + } + printf("================================================================================\n"); +} + +const std::vector detNum = ExtractDetNum(); +const int nDetType = detNum.size(); + +const int NARRAY = detNum[0]; //@ assumed +const int NRDT = detNum[3]; //@ assumed + +} // namespace solarismap +#endif diff --git a/Working/Monitor.C b/Working/Monitor.C new file mode 100644 index 0000000..a1a8a61 --- /dev/null +++ b/Working/Monitor.C @@ -0,0 +1,486 @@ + +#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, 50}; /// deg + +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 skipDetID = {11} ; + +std::vector 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 Monitor(TChain *gen_tree){ + + printf("#####################################################################\n"); + printf("####################### Monitor.C #######################\n"); + printf("#####################################################################\n"); + + TObjArray * fileList = gen_tree->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", gen_tree->GetEntries()); + printf("///////////////////////////////////////////////////////////////////\n"); + + if( gen_tree->GetEntries() == 0 ) { + printf("========= no events. Abort.\n"); + return; + } + + double totDuration = 0; + std::vector startTime; + std::vector stopTime; + std::vector runList; + + for( int i = 0; i < fileList->GetEntries(); i++){ + TString fileName = fileList->At(i)->GetTitle(); + TFile file(fileName); + TMacro * timeStamp = (TMacro*) file.FindObjectAny("timeStamp"); + //timeStamp->Print(); + + TString haha = timeStamp->GetListOfLines()->At(0)->GetName(); + ULong64_t t1 = haha.Atoll(); + + haha = timeStamp->GetListOfLines()->At(1)->GetName(); + ULong64_t t2 = haha.Atoll(); + + haha = timeStamp->GetListOfLines()->At(2)->GetName(); + int RunID = haha.Atoi(); + + totDuration += (t2-t1)*8./1e9; + startTime.push_back(t1); + stopTime.push_back(t2); + runList.push_back(RunID); + } + + //======== format CanvasTitle + std::sort(runList.begin(), runList.end()); + TString title = "Run:" + AnalysisLib::create_range_string(runList); + title += Form(" | %.0f min", totDuration/60.) ; + + //*=========================================================== + TTreeReader reader(gen_tree); + + TTreeReaderValue evID = {reader, "evID"}; + TTreeReaderArray e = {reader, "e"}; + TTreeReaderArray e_t = {reader, "e_t"}; + TTreeReaderArray xf = {reader, "xf"}; + TTreeReaderArray xn = {reader, "xn"}; + + TTreeReaderArray rdt = {reader, "rdt"}; + TTreeReaderArray rdt_t = {reader, "rdt_t"}; + + //TODO + // TTreeReaderArray array = {reader, "trace"}; + + ULong64_t NumEntries = gen_tree->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(title, 500, 3, 2); + 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::vectoreCal (numTotArray); + std::vectorxfCal (numTotArray); + std::vectorxnCal (numTotArray); + std::vectorx (numTotArray); + std::vectorxCal (numTotArray); + std::vectorz (numTotArray); + + + //^########################################################### + //^ * Process + //^########################################################### + printf("############################################### Processing...\n"); + fflush(stdout); // flush out any printf + + 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; + // bool rdtgate2 = false; + bool coinFlag = false; + bool isGoodEventFlag = 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] ++; + // printf("%8llu | %d, %f %f %f | \n", processedEntries, id, e[id], xn[id], xf[id] ); + + //@==================== Calibrations go here + xnCal[id] = xn[id] * corr->xnCorr[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; + xfCal[id] = xf[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; + 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 (0,1) + 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(e[id], x[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) + 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 < plotter[aID]->cutList->GetEntries() ; i++ ){ + TCutG * cutG = (TCutG *)plotter[aID]->cutList->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])) { + if((rdtgate1 ) && (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 ) plotter[aID]->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); + + // plotter[aID]->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; + // } + + //@================ is Good Event ? + // if( coinFlag && (rdtgate1 || rdtgate2) && ezGate){ + if( coinFlag && rdtgate1){ + plotter[aID]->heCal_zGC->Fill( z[id] , eCal[id] ); + + 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 ****************************************/ + if( !isGoodEventFlag ) continue; + + 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 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); + fflush(stdout); + 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(" rdt() - Raw RDT data\n"); + //printf(" elum() - Luminosity Energy Spectra\n"); + printf("-----------------------------------------------------\n"); + printf(" ez() - Energy vs. Z\n"); + printf("-----------------------------------------------------\n"); + printf(" excite() - Excitation Energy\n"); + // 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", title.Data()); + printf("-----------------------------------------------------\n"); + +} + +//%============================================= +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); + } +} + +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 rdt(bool isLog = false){ + plotter[0]->PlotRDT(isLog); +} + +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 excited(int arrayID = -1){ + if( arrayID < 0 ){ + for( int i = 0; i < numGeo; i++ ) plotter[i]->PlotEx(); + }else{ + if( arrayID < numGeo) plotter[arrayID]->PlotEx(); + } +} \ No newline at end of file diff --git a/Working/SimHelper.C b/Working/SimHelper.C new file mode 120000 index 0000000..a4868ba --- /dev/null +++ b/Working/SimHelper.C @@ -0,0 +1 @@ +../Cleopatra/SimHelper.C \ No newline at end of file diff --git a/Working/obsolete/DataHoSei.C b/Working/obsolete/DataHoSei.C new file mode 100644 index 0000000..6df3ce8 --- /dev/null +++ b/Working/obsolete/DataHoSei.C @@ -0,0 +1,57 @@ +#define DataHoSei_cxx + +#include "DataHoSei.h" +#include +#include +#include +#include + +//^***************************** User Settings + + + + +//******************************** end of User Settings + +TH1F * h1 = new TH1F("h1", "h1", 400, 0, 50000); + + +void DataHoSei::Begin(TTree * /*tree*/){ + printf("--------- %s \n", __func__); + TString option = GetOption(); + +} + + +Bool_t DataHoSei::Process(Long64_t entry){ + + if( entry == 0 ) printf("--------- %s \n", __func__); + + b_e[0]->GetEntry(entry); + //printf("%f \n", e[0][1]); + + h1->Fill(e[0][1]); + + // std::vector> dataList; // {detID, e, xf, xn} + // for( int i = 0; i < mapping::nDetType; i++){ + std::vector temp; + // for( int j = 0; j < mapping::detNum[i]; j++){ + // if( TMath::IsNaN( e[i][j] ) ) continue; + // } + // } + + return kTRUE; +} + + + + +void DataHoSei::Terminate(){ + printf("--------- %s \n", __func__); + + h1->Draw(); + + + + printf("----------- A \n"); +} diff --git a/Working/obsolete/DataHoSei.h b/Working/obsolete/DataHoSei.h new file mode 100644 index 0000000..50f4835 --- /dev/null +++ b/Working/obsolete/DataHoSei.h @@ -0,0 +1,122 @@ +////////////////////////////////////////////////////////// +// +// Made by Ryan Tang on 2023, April 6 +// +// +////////////////////////////////////////////////////////// + +#ifndef DataHoSei_h +#define DataHoSei_h + +#include +#include +#include +#include + +// Header file for the classes stored in the TTree if any. +#include "TClonesArray.h" +#include "TObject.h" +#include "TNamed.h" +#include "TAttLine.h" +#include "TAttFill.h" +#include "TAttMarker.h" +#include "TGraph.h" + +#include "Mapping.h" +#include "../armory/AnalysisLib.h" + +class DataHoSei : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + + // Declaration of leaf types + ULong64_t evID; + + Float_t ** e; //! all kind of energy + ULong64_t ** e_t; //! all kind of timestamp + + Float_t ** we; //! wave energy + Float_t ** weT; //! wave time + Float_t ** weR; //! wave rise time + + // List of branches + TBranch *b_evID; //! + TBranch **b_e; //! + TBranch **b_e_t; //! + + DataHoSei(TTree * /*tree*/ =0) : fChain(0) { + printf("--------- %s \n", __func__); + + e = new Float_t * [mapping::nDetType]; + e_t = new ULong64_t * [mapping::nDetType]; + b_e = new TBranch * [mapping::nDetType]; + b_e_t = new TBranch * [mapping::nDetType]; + for( int i = 0 ; i < mapping::nDetType; i++) { + e[i] = new Float_t[mapping::detNum[i]]; + e_t[i] = new ULong64_t[mapping::detNum[i]]; + } + + } + virtual ~DataHoSei() { + + printf("--------- %s \n", __func__); + printf("----------- B \n"); + } + 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(); + + ClassDef(DataHoSei,0); + +}; + +#endif + +#ifdef DataHoSei_cxx + +void DataHoSei::Init(TTree *tree){ + + printf("--------- %s \n", __func__); + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("evID", &evID, &b_evID); + + for( int i = 0; i < mapping::nDetType; i++){ + fChain->SetBranchAddress((mapping::detTypeName[i]).c_str(), e[i], &b_e[i]); + fChain->SetBranchAddress((mapping::detTypeName[i] + "_t").c_str(), e_t[i], &b_e_t[i]); + } + + // TObjArray * branchList = fChain->GetListOfBranches(); + + +} + +void DataHoSei::SlaveBegin(TTree * /*tree*/){ + printf("--------- %s \n", __func__); + TString option = GetOption(); +} + +void DataHoSei::SlaveTerminate(){ + printf("--------- %s \n", __func__); +} + +Bool_t DataHoSei::Notify(){ + printf("--------- %s \n", __func__); + return kTRUE; +} + +#endif // #ifdef DataHoSei_cxx diff --git a/Working/obsolete/Monitor.C b/Working/obsolete/Monitor.C new file mode 100644 index 0000000..86604e6 --- /dev/null +++ b/Working/obsolete/Monitor.C @@ -0,0 +1,710 @@ +#define Monitor_cxx + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "../Cleopatra/ClassIsotope.h" +#include "Mapping.h" + +#define tick2ns 8. // 1clock tick = 8 ns +#define tick2min tick2ns/1e9/60. + +using namespace std; + +//############################################ User setting +//---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 + + +// //======= Recoil +// TH2F * hrdtID; +// TH1F ** hrdt; // single recoil +// TH1F ** hrdtg; + +// TH2F ** hrdt2D; +// TH2F ** hrdt2Dg; + +// TH1F * hrdtRate1; +// TH1F * hrdtRate2; + +// //======= multi-Hit +// TH2I * hmult; +// TH1I * hmultEZ; +// TH2I * hArrayRDTMatrix; +// TH2I * hArrayRDTMatrixG; + +// //======= ARRAY-RDT time diff +// TH1I * htdiff; +// TH1I * htdiffg; + +/*************************** + ***************************/ +TLatex text; + +ULong64_t NumEntries = 0; +ULong64_t ProcessedEntries = 0; +Float_t Frac = 0.1; ///Progress bar +TStopwatch StpWatch; + +//======= Recoil Cut +TCutG* cutG; //! //general temeprary pointer to cut + +TObjArray * cutList1; +TObjArray * cutList2; + +//======= Other Cuts +TCutG * EZCut; + +#include "Monitor.h" + +//^########################################################### +//^ * Begin +//^########################################################### +void Monitor::Begin(TTree *tree){ + + TString option = GetOption(); + + NumEntries = tree->GetEntries(); + canvasTitle = GetCanvasTitle(); + + printf("###########################################################\n"); + printf("########## SOLARIS Monitors.C #########\n"); + printf("###########################################################\n"); + + for( int i = 0; i < detGeo->numGeo ; i++) plotter[i]->SetUpHistograms(energyRange, exRange, thetaCMRange); + + //===================================================== loading parameter + + // corr->LoadDetGeoAndReactionConfigFile(); + corr->LoadXNCorr(); + corr->LoadXFXN2ECorr(); + corr->LoadXScaleCorr(); + corr->LoadECorr(); + corr->LoadRDTCorr(); + + if( (int) corr->xnCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); } + if( (int) corr->xfxneCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); } + if( (int) corr->eCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of eCorr < NARRAY .\n"); } + if( (int) corr->xScale.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xScale < NARRAY .\n"); } + if( (int) corr->rdtCorr.size() < mapping::NRDT ) { printf(" !!!!!!!! size of rdtCorr < NRDT .\n"); } + + + // printf("=====================================================\n"); + // printf(" time Range : %5.0f - %5.0f min\n", timeRangeInMin[0], timeRangeInMin[1]); + // printf("=====================================================\n"); + + //================ 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(); + + + // int startIndex = 0; + // for( int i = 0; i < detGeo->numGeo; i++ ){ + // CreateListOfHist2D(heVx , startIndex, detGeo->array[i].numDet, "heVx", "Raw PSD E vs. X (ch=%d);X (channel);E (channel)", 500, -2.5, detGeo->array[i].detLength + 2.5, 500, rawEnergyRange[0], rawEnergyRange[1]); + // CreateListOfHist2D(heCalVxCal , startIndex, detGeo->array[i].numDet, "heCalVxCal", "Cal PSD E vs. X (ch=%d);X (cm);E (MeV)", 500, -2.5, detGeo->array[i].detLength + 2.5, 500, energyRange[0], energyRange[1]); + // CreateListOfHist2D(heCalVxCalG , startIndex, detGeo->array[i].numDet, "heCalVxCalG", "Cal PSD E vs. X (ch=%d);X (cm);E (MeV)", 500, -2.5, detGeo->array[i].detLength + 2.5, 500, energyRange[0], energyRange[1]); + // startIndex += detGeo->array[i].numDet; + // } + + // 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 id; id; E / 10 keV", mapping::NARRAY, 0, mapping::NARRAY, 2000, energyRange[0], energyRange[1]); + + // hMultiHit = new TH1F("hMultiHit", "Multi-Hit of Energy", 10, 0, 1); + + // //===================== 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]; + + // 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]); + // } + // } + + // hrdtRate1 = new TH1F("hrdtRate1", "recoil rate 1 / min; min; count / 1 min", timeRangeInMin[1] - timeRangeInMin[0], timeRangeInMin[0], timeRangeInMin[1]); + // hrdtRate2 = new TH1F("hrdtRate2", "recoil rate 2 / min; min; count / 1 min", timeRangeInMin[1] - timeRangeInMin[0], timeRangeInMin[0], timeRangeInMin[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 Monitor::Process(Long64_t entry){ + + if( entry == 0 ) { + treeID ++; + baseTimeStamp = (treeID == 0 ? 0 : endTime[treeID-1]); + printf("============================================ %s , treeID : %d\n", __func__, treeID); + } + + 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( corr->rdtCorr.size() >= mapping::NRDT ){ + for( int i = 0 ; i < mapping::NRDT; i++){ + rdt[i] = corr->rdtCorr[i][0] + rdt[i]*corr->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 id = 0; id < mapping::NARRAY; id++) { + + //@================== Filling raw data + he[id]->Fill(e[id]); + hxf[id]->Fill(xf[id]); + hxn[id]->Fill(xn[id]); + hxfVxn[id]->Fill(xf[id],xn[id]); + heVxs[id]->Fill(xf[id]+xn[id], e[id]); + + heVID->Fill(id, e[id]); + hxfVID->Fill(id, xf[id]); + hxnVID->Fill(id, xn[id]); + + //if( !TMath::IsNaN(e[id]) ) printf("%llu | %d | %f %f %f \n", entry, id, e[id], xf[id], xn[id]); + + //@==================== 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 ; + + //@==================== Skip detector + bool skipFlag = false; + for( unsigned int kk = 0; kk < skipDetID.size() ; kk++){ + if( id == skipDetID[kk] ) { + skipFlag = true; + break; + } + } + if (skipFlag ) 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]; + if( corr->eCorr.size() >= id) 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 + heCal[id]->Fill(eCal[id]); + heCalID->Fill(id, eCal[id]); + hxfCalVxnCal[id]->Fill(xfCal[id], xnCal[id]); + heVxsCal[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]; + + //@======= Scale xcal from (0,1) + if( corr->xScale.size() >= id ) 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 + short arrayID = detGeo->GetArrayID(id); + if( arrayID >= 0 ){ + int colIndex = id % detGeo->array[arrayID].colDet; + if( detGeo->array[arrayID].firstPos > 0 ) { + z[id] = detGeo->array[arrayID].detLength*(1.0-xCal[id]) + detGeo->array[arrayID].detPos[colIndex]; + }else{ + z[id] = detGeo->array[arrayID].detLength*(xCal[id]-1.0) + detGeo->array[arrayID].detPos[colIndex]; + } + } + + //@===================== multiplicity + arrayMulti++; /// multi-hit when both e, xf, xn are not NaN + + //@=================== Array fill + heVx[id]->Fill(x[id],e[id]); + + heCalVxCal[id]->Fill(xCal[id]*detGeo->array[arrayID].detLength, eCal[id]); + + plotter[arrayID]->heCalVz->Fill(z[id],eCal[id]); + + //@=================== 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[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) { + htdiff->Fill(tdiff); + if((rdtgate1 || rdtgate2) && (eCalCut[1] > eCal[id] && eCal[id]>eCalCut[0])) { + 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 + if( EZCut == nullptr ) 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 + } + } + + hrdtRate1->Fill( (e_t[1] + baseTimeStamp) * tick2min ); //incorrect + + //@******************* 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 id = 0; id < mapping::NARRAY ; 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 arrayID = detGeo->GetArrayID(id); + if( arrayID < 0 ) continue; + + std::pair ExThetaCM = transfer->CalExThetaCM(eCal[id], z[id], detGeo->Bfield, detGeo->array[arrayID].detPerpDist); + double Ex = ExThetaCM.first; + double thetaCM = ExThetaCM.second; + + if( thetaCM > thetaCMGate ) { + + plotter[arrayID]->hEx->Fill(Ex); + plotter[arrayID]->hExThetaCM->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); + } + + plotter[arrayID]->hExi[id]->Fill(Ex); + plotter[arrayID]->hExVxCal[id]->Fill(xCal[id], Ex); + + } + } + + return kTRUE; +} + +//^########################################################### +//^ * Terminate +//^########################################################### +void Monitor::Terminate(){ + printf("============================================ Drawing Canvas.\n"); + + gROOT->cd(); + + //############################################ User is free to edit this section + //--- Canvas Size + int canvasXY[2] = {1200 , 800} ;// x, y + int canvasDiv[2] = {3,2}; + + 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); + + for( int i = 0; i < detGeo->numGeo; i++ ){ + + plotter[i]->SetUpCanvas(canvasTitle + " | " + rdtCutFile1, canvasXY[0],canvasXY[1], canvasDiv[0], canvasDiv[1]); + // cCanvas[i] = new TCanvas("cCanvas",canvasTitle + " | " + rdtCutFile1,canvasXY[0],canvasXY[1]); + // cCanvas[i]->Modified(); cCanvas->Update(); + // cCanvas[i]->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]); + + plotter[i]->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", timeRangeInMin[1]-timeRangeInMin[0], timeRangeInMin[0], timeRangeInMin[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); + + + } + + /************************************/ + gStyle->GetAttDate()->SetTextSize(0.02); + gStyle->SetOptDate(1); + gStyle->SetDateX(0); + gStyle->SetDateY(0); + + /************************************/ + StpWatch.Start(kFALSE); + + // gROOT->ProcessLine(".L ../armory/Monitor_Util.C"); //TODO some pointer is empty + // 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/obsolete/Monitor.h b/Working/obsolete/Monitor.h new file mode 100644 index 0000000..c2bc231 --- /dev/null +++ b/Working/obsolete/Monitor.h @@ -0,0 +1,392 @@ +#ifndef Monitor_h +#define Monitor_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Mapping.h" +#include "../Armory/AnalysisLib.h" +#include "../Armory/ClassDetGeo.h" +#include "../Armory/ClassReactionConfig.h" +#include "../Armory/ClassCorrParas.h" +#include "../Cleopatra/ClassTransfer.h" + +#include "ClassMonPlotter.h" + +class Monitor : 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; + + CorrParas * corr; //! + DetGeo * detGeo; //! + TransferReaction * transfer; //! + + TString canvasTitle; + MonPlotter ** plotter; //! + + //==== global variable + float * x, * z; + float * xCal, * xfCal, * xnCal, * eCal; + + std::vector startTime ; + std::vector endTime ; + + double timeRangeInMin[2]; + ULong64_t baseTimeStamp; + int treeID; + + int padID; + + Monitor(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]; + + padID = 0; + + timeRangeInMin[0] = 0; + timeRangeInMin[1] = 100; + + startTime.clear(); + endTime.clear(); + + baseTimeStamp = 0; + treeID = -1; + + corr = new CorrParas(); + detGeo = new DetGeo(); + detGeo->LoadDetectorGeo("detectorGeo.txt"); + transfer = new TransferReaction(); + transfer->SetReactionFromFile("reactionConfig.txt"); + + plotter = new MonPlotter *[detGeo->numGeo]; + for( int i = 0; i < detGeo->numGeo; i++ ){ + plotter[i] = new MonPlotter(i, detGeo); + } + } + virtual ~Monitor() { + + 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; + + delete corr; + + for( int i = 0; i < detGeo->numGeo; i++ ) delete [] plotter[i]; + delete plotter; + delete detGeo; + delete transfer; + + } + 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;} + + // void Draw2DHist(TH2F * hist); + // void PlotEZ(bool isRaw); + // void PlotTDiff(bool isGated, bool isLog); + // void PlotRDT(int id, bool isRaw); + //void PlotCRDTPolar(); + + + ClassDef(Monitor,0); +}; + +#endif + + +#ifdef Monitor_cxx +void Monitor::Init(TTree *tree){ + + printf("============================================ Branch Pointer Inititization. \n"); + + // 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 Branch.\n"); + isRDTExist = false; + }else{ + printf(" ++++++++ Found Recoil Branch.\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); + } + */ + + printf("============================================ End of Branch Pointer Inititization. \n"); +} + +Bool_t Monitor::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 Monitor::SlaveBegin(TTree * /*tree*/){ + /// not use, if use, place in Monitor.C + TString option = GetOption(); +} + + +void Monitor::SlaveTerminate(){ + /// not use, if use, place in Monitor.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 Monitor::Draw2DHist(TH2F * hist){ + +// if( hist->Integral() < 3000 ){ +// hist->SetMarkerStyle(20); +// hist->SetMarkerSize(0.3); +// hist->Draw(""); +// }else{ +// hist->Draw("colz"); +// } +// } + +/* +void Monitor::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 Monitor::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 Monitor::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 Monitor::PlotCRDTPolar(){ +// padID++; cCanvas->cd(padID); +// cCanvas->cd(padID)->DrawFrame(-50, -50, 50, 50); +// hcrdtPolar->Draw("same colz pol"); +//} +*/ + +#endif // #ifdef Monitor_cxx diff --git a/Working/obsolete/Monitor_Util.C b/Working/obsolete/Monitor_Util.C new file mode 100644 index 0000000..fa817a7 --- /dev/null +++ b/Working/obsolete/Monitor_Util.C @@ -0,0 +1,359 @@ +#ifndef Utilities +#define Utilities + +#include +#include + +//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", numDet); + printf(" rawxf() - Raw \033[0;31mxf\033[0m for all %d detectors\n", numDet); + printf(" rawxn() - Raw \033[0;31mxn\033[0m for all %d detectors\n", numDet); + printf(" rawxfVxn() - Raw \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", numDet); + printf(" raweVxs() - Raw \033[0;31me\033[0m vs. Raw \033[0;31mxs = xf + xn\033[0m for all %d detectors\n", numDet); + printf(" raweVx() - Raw \033[0;31me\033[0m vs. RAW \033[0;31mx\033[0m for all %d detectors\n", numDet); + printf("-----------------------------------------------------\n"); + printf(" eVxsCal() - Raw \033[0;31me\033[0m vs. Corrected \033[0;31mxs\033[0m for all %d detectors\n", numDet); + printf(" ecal() - Calibrated \033[0;31me\033[0m for all %d detectors\n", numDet); + printf("xfCalVxnCal() - Calibrated \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", numDet); + printf(" eCalVxCal() - Cal \033[0;31me\033[0m vs. \033[0;31mx\033[0m for all %d detectors\n", numDet); + 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", 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"); +} + +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(); +} + +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 < numDet; 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 rawxfVxn(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 raweVxs(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 raweVx(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 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 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 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 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 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 *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(""); + } +} + +// 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 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 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 < numDet; 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 < numDet; i++){ + cExVxCal->cd(i+1); + if( drawOpt == "" )hExVxCal[i]->SetMarkerStyle(7); + hExVxCal[i]->Draw(drawOpt); + } +} + +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/obsolete/script_multi.C b/Working/obsolete/script_multi.C new file mode 100644 index 0000000..37aa464 --- /dev/null +++ b/Working/obsolete/script_multi.C @@ -0,0 +1,169 @@ +#include "../armory/SolReader.h" +#include "TH1.h" +#include "TMath.h" +#include "TH2.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TGraph.h" + + +void script_multi(std::string run = "033"){ + + + SolReader* reader0 = new SolReader("../data_raw/*_" + run + "_00_21245_000.sol"); + //SolReader* reader1 = new SolReader("/home/ryan/analysis/data_raw/test_" + run + "_01_21233_000.sol"); + + Hit * evt0 = reader0->hit; + //Event * evt1 = reader1->evt; + + printf("----------file size: %u Byte\n", reader0->GetFileSize()); + //printf("----------file size: %u Byte\n", reader1->GetFileSize()); + + reader0->ScanNumBlock(); + //reader1->ScanNumBlock(); + unsigned long startTime, endTime; + + if( reader0->GetTotalNumBlock() == 0 ) return; + reader0->ReadBlock(0); + startTime = evt0->timestamp; + reader0->ReadBlock(reader0->GetTotalNumBlock() - 1); + endTime = evt0->timestamp; + + double duration = double(endTime - startTime)*8./1e9; + printf("============== %lu ns = %.4f sec.\n", (endTime - startTime)*8, duration); + reader0->RewindFile(); + + // if( reader1->GetTotalNumBlock() == 0 ) return; + // reader1->ReadBlock(0); + // startTime = evt1->timestamp; + // reader1->ReadBlock(reader1->GetTotalNumBlock() - 1); + // endTime = evt1->timestamp; + + // duration = double(endTime - startTime)*8./1e9; + // printf("============== %lu ns = %.4f sec.\n", (endTime - startTime)*8, duration); + // reader1->RewindFile(); + + //int minBlock = std::min(reader0->GetTotalNumBlock(), reader1->GetTotalNumBlock()); + int minBlock = reader0->GetTotalNumBlock(); + + TH1F * hID = new TH1F("hID", "hID", 64, 0, 64); + TH1F * hTdiff = new TH1F("hTdiff", "tdiff", 300, -100, 200); + TH2F * hTdiff2 = new TH2F("hTdiff2", "tdiff vs evt", 400, 0, minBlock, 100, 0, 100); + + TH1F * hRate0 = new TH1F("hRate0", "Rate", 20, 0, 20); + TH1F * hRate1 = new TH1F("hRate1", "Rate", 20, 0, 20); hRate1->SetLineColor(2); + + TH1F * hE0 = new TH1F("hE0", "Energy", 400, 0, 30000); + TH1F * hE1 = new TH1F("hE1", "Energy", 400, 0, 30000); hE1->SetLineColor(2); + + TH1F *hMulti = new TH1F("hMulti", "Multiplicy", 10, 0, 10); + + std::vector> ts ; + + for( int i = 0; i < minBlock; i++){ + + + reader0->ReadNextBlock(); + //reader1->ReadNextBlock(); + + if( i < 10 ){ + printf("#################################################\n"); + evt0->PrintAll(); + } + + ts.push_back(std::pair(evt0->channel,evt0->timestamp)); + // printf("---------------------------------------------\n"); + // evt1->PrintAll(); + + //if( evt0->channel == 30 ) evt0->PrintAll(); + + hID->Fill(evt0->channel); + + if( evt0->channel == 0 ) { + hRate0->Fill( evt0->timestamp * 8 / 1e9); + hE0->Fill(evt0->energy); + } + if( evt0->channel == 30 ) { + hRate1->Fill( evt0->timestamp * 8 / 1e9); + hE1->Fill(evt0->energy); + } + //if( i < 10 ) printf("t0 : %10lu, t1 : %10lu, %10lu \n", evt0->timestamp, evt1->timestamp, + // evt0->timestamp > evt1->timestamp ? evt0->timestamp - evt1->timestamp : evt1->timestamp - evt0->timestamp); +// + //hTdiff->Fill(evt0->timestamp > evt1->timestamp ? evt0->timestamp - evt1->timestamp : evt1->timestamp - evt0->timestamp); + //hTdiff2->Fill(i, evt0->timestamp > evt1->timestamp ? evt0->timestamp - evt1->timestamp : evt1->timestamp - evt0->timestamp); + + //if( i > 10 ) break; + + } + + delete reader0; + //delete reader1; + + //build event + int coinWin = 200; + + std::vector>> events; + std::vector multi; + events.push_back({ts[0]}); + multi.push_back(1); + + // Iterate through the vector starting from the second pair + for (size_t i = 1; i < ts.size(); ++i) { + uint64_t currentTimestamp = ts[i].second; + uint64_t previousTimestamp = ts[i - 1].second; + + // Check if the timestamp difference is within coinWin + if (currentTimestamp - previousTimestamp <= coinWin) { + events.back().push_back(ts[i]); // Add to the current event + multi.back() += 1; + } else { + // Create a new event + events.push_back({ts[i]}); + multi.push_back(1); + } + } + + // Print the events + + double maxTD = -999, minTD = 999; + + for (size_t i = 0; i < events.size(); ++i) { + hMulti->Fill(multi[i]); + if( multi[i] > 1 ) { + printf("Event %zu, Multi : %d\n", i, multi[i]); + double haha[2]; + for (size_t j = 0; j < events[i].size(); ++j) { + printf("Channel: %2d, Timestamp: %lu\n", events[i][j].first, events[i][j].second); + if( events[i][j].first == 0 ) haha[0] = static_cast ( events[i][j].second ); + if( events[i][j].first == 30 ) haha[1] = static_cast ( events[i][j].second ); + } + + double TD = haha[1]-haha[0]; + + if(TD > maxTD && TD < 105) maxTD = TD; + if( TD < minTD && TD > -105) minTD = TD; + hTdiff->Fill( TD ); + printf("-------------------\n"); + } + } + + printf(" max TD : %f \n", maxTD); + printf(" min TD : %f \n", minTD); + printf(" spn TD : %f \n", maxTD - minTD); + + TCanvas * canvas = new TCanvas("canvas", "canvas", 1200, 800); + + gStyle->SetOptStat("neiou"); + canvas->Divide(3,2); + canvas->cd(1);hID->Draw(); + canvas->cd(2);hE0->Draw(); hE1->Draw("same"); + canvas->cd(3);hRate0->Draw(); hRate1->Draw("same"); + canvas->cd(4);hTdiff->Draw(); + canvas->cd(5); hMulti->Draw(); + canvas->cd(5)->SetLogy(true); + canvas->cd(6); hE1->Draw(); + //canvas->cd(2);hTdiff2->Draw(); + + +} \ No newline at end of file diff --git a/Working/obsolete/script_single.C b/Working/obsolete/script_single.C new file mode 100644 index 0000000..99c9c3a --- /dev/null +++ b/Working/obsolete/script_single.C @@ -0,0 +1,89 @@ +#include "../armory/SolReader.h" +#include "TH1.h" +#include "TH2.h" +#include "TCanvas.h" + +void script_single(std::string fileName){ + + SolReader * reader = new SolReader(fileName); + Hit * hit = reader->hit; + + reader->ScanNumBlock(); + + long numBlock = reader->GetTotalNumBlock(); + + for( int i = 1; i < 10; i ++ ){ + reader->ReadBlock(numBlock-i); + hit->PrintEnergyTimeStamp(); + } + + /* + SolReader * reader = new SolReader("../data_raw/Master_003_00_21245_000.sol"); + Event * hit = reader->hit; + + printf("========== file size: %u Byte\n", reader->GetFileSize()); + + reader->ScanNumBlock(); + + if( reader->GetTotalNumBlock() == 0 ) return; + + unsigned long startTime, endTime; + reader->ReadBlock(0); + startTime = hit->timestamp; + reader->ReadBlock(reader->GetTotalNumBlock() - 1); + endTime = hit->timestamp; + + double duration = double(endTime - startTime)*8./1e9; + printf("============== %lu ns = %.4f sec.\n", (endTime - startTime)*8, duration); + printf(" avarge rate (16ch): %f Hz\n", reader->GetTotalNumBlock()/duration/16); + reader->RewindFile(); + + + TCanvas * canvas = new TCanvas("c1", "c1", 600, 600); + + TH1F * h1 = new TH1F("h1", "h1", 1000, startTime, endTime); + TH2F * h2 = new TH2F("h2", "h2", 1000, startTime, endTime, 1000, 0, reader->GetTotalNumBlock()); + + uint64_t tOld = startTime; + + for( int i = 0; i < reader->GetTotalNumBlock() ; i++){ + //for( int i = 0; i < 8 ; i++){ + + reader->ReadNextBlock(); + + if( i < 8 ){ + printf("########################## nBlock : %u, %u/%u\n", reader->GetBlockID(), + reader->GetFilePos(), + reader->GetFileSize()); + hit->PrintAll(); + //hit->PrintAllTrace(); + } + + h1->Fill(hit->timestamp); + h2->Fill(hit->timestamp, i); + + if( i > 0 ){ + if( hit->timestamp < tOld) printf("-------- time not sorted."); + tOld = hit->timestamp; + } + + } + + + h2->Draw(); + */ + + //printf("reader traceLength : %lu \n", hit->traceLenght); + + /* + for( int i = 0; i < hit->traceLenght; i++){ + + printf("%4d| %d\n", i, hit->analog_probes[0][i]); + + } + */ + + hit = NULL; + delete reader; + +} \ No newline at end of file diff --git a/Working/test.C b/Working/test.C new file mode 100644 index 0000000..aa9d502 --- /dev/null +++ b/Working/test.C @@ -0,0 +1,104 @@ +#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(){ + + // DetGeo haha("detectorGeo.txt"); + // haha.Print(true); + + // ReactionConfig config("reactionConfig.txt"); + // config.Print(); + + // TransferReaction * transfer = new TransferReaction(); + // transfer->SetReactionSimple(32, 14, 2, 1, 1, 1, 8.8); + + // int ID = 0; + // transfer->SetReactionFromFile("reactionConfig.txt", ID); + + // transfer->PrintReaction(); + + + // transfer->Event(25 * TMath::DegToRad(), 0 * TMath::DegToRad()); + // transfer->PrintFourVectors(); + + // ReactionConfig config2 = transfer->GetRectionConfig(); + + // HELIOS * helios = new HELIOS(); + // helios->SetDetectorGeometry("detectorGeo.txt", 1); + // helios->PrintGeometry(); + + + // TLorentzVector Pb = transfer->GetPb(); + // printf("Charge : %d\n", Pb.GetUniqueID()); + // int hit = helios->CalArrayHit(Pb, false); + // helios->CheckDetAcceptance(); + + // // //helios->CalTrajectoryPara(Pb, config.recoilLightZ, true); + + // printf("\n hit = %d | %s | %s\n", hit, helios->GetHitMessage().Data(), helios->GetAcceptanceMessage().Data()); + + // trajectory orb = helios->GetTrajectory_b(); + + // orb.PrintTrajectory(); + + // delete helios; + // delete transfer;root + + // DetGeo dd("detectorGeo.txt"); + // MonPlotter * pp = new MonPlotter(0, &dd, 8); + + // 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 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}; + + // pp->SetUpHistograms(rawEnergyRange, energyRange, exRange, thetaCMRange, rdtDERange, rdtERange); + + // pp->Plot(); + + // delete pp; + + // TChain *chain = new TChain("gen_tree"); + // chain->Add("../root_data/gen_run043.root"); + + // // chain->Print(); + + // TTreeReader reader(chain); + + // TTreeReaderArray 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(); + +} \ No newline at end of file