From 08577871ee5c89ae9cc2fa80a180b82f5ec0d83d Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Wed, 10 Jul 2024 16:18:46 -0400 Subject: [PATCH] MonAnalyser.C basically done. Need to link to ChainMonitors.C, need RDTCutCreator, need Calibrations stuffs. --- working/ClassMonPlotter.h | 87 +++++++++++++------ working/MonAnalyzer.C | 93 ++++++++++++--------- working/{ => obsolete}/DataHoSei.C | 0 working/{ => obsolete}/DataHoSei.h | 0 working/{ => obsolete}/Monitor.C | 0 working/{ => obsolete}/Monitor.h | 0 {Armory => working/obsolete}/Monitor_Util.C | 0 working/{ => obsolete}/script_multi.C | 0 working/{ => obsolete}/script_single.C | 0 9 files changed, 114 insertions(+), 66 deletions(-) rename working/{ => obsolete}/DataHoSei.C (100%) rename working/{ => obsolete}/DataHoSei.h (100%) rename working/{ => obsolete}/Monitor.C (100%) rename working/{ => obsolete}/Monitor.h (100%) rename {Armory => working/obsolete}/Monitor_Util.C (100%) rename working/{ => obsolete}/script_multi.C (100%) rename working/{ => obsolete}/script_single.C (100%) diff --git a/working/ClassMonPlotter.h b/working/ClassMonPlotter.h index 2717407..c002088 100644 --- a/working/ClassMonPlotter.h +++ b/working/ClassMonPlotter.h @@ -72,6 +72,8 @@ public: void PlotEZ(); void PlotEx(); + void PlotRDT(bool isLog = false); + TCanvas * canvas; //====================== Histograms @@ -148,7 +150,7 @@ MonPlotter::MonPlotter(unsigned short arrayID, DetGeo * detGeo, int numRDT){ suffix = Form("_%d", arrayID); - this->numRDT = numDet; + this->numRDT = numRDT; recoilOutter = detGeo->aux[aID].outerRadius; zRange[0] = detGeo->array[aID].zMin - 50; @@ -218,7 +220,6 @@ MonPlotter::~MonPlotter(){ delete [] hrdt2D; delete [] hrdt2Dg; - delete cutG; delete cutList; } @@ -271,7 +272,7 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, 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", numDet, 0, numDet); + 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]); @@ -280,7 +281,7 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, CreateListOfHist2D(hxf_xn, numDet, "hxf_xn", "Raw xf vs. xn (ch=%d);xf (channel);xn (channel)" , 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); CreateListOfHist2D(he_xs, numDet, "he_xs", "Raw e vs xf+xn (ch=%d); xf+xn (channel); e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, rawEnergyRange[0], rawEnergyRange[1]); - CreateListOfHist2D(he_x , numDet, "he_x", "Raw e vs x (ch=%d); x (mm); Raw e (channel)", 500, rawEnergyRange[0], rawEnergyRange[1], 500, -1, detLength +1); + 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]); @@ -347,12 +348,19 @@ void MonPlotter::SetUpHistograms(int * rawEnergyRange, } 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("");break; + case 3: { + htDiff->Draw(""); + htDiffg->Draw("same"); + }break; case 4: hEx->Draw("colz");break; default:break; } @@ -366,7 +374,7 @@ void MonPlotter::LoadRDTGate(TString rdtCutFile){ TFile * fCut = new TFile(rdtCutFile); bool isCutFileOpen = fCut->IsOpen(); if(!isCutFileOpen) { - printf( "Failed to open rdt-cutfile 1 : %s\n" , fileName.Data()); + printf( "Failed to open rdt-cutfile 1 : %s\n" , rdtCutFile.Data()); }else{ cutList = (TObjArray *) fCut->FindObjectAny("cutList"); @@ -461,12 +469,13 @@ void MonPlotter::PlotCal(){ TCanvas *cxfxneC = new TCanvas("cxfxneC" + suffix,Form("Raw E - Corrected XF+XN | %s", canvasTitle.Data()), 200 + 500 * aID, 200, canvasSize[0], canvasSize[1]); cxfxneC->Clear(); cxfxneC->Divide(colDet,rowDet); - TLine line(0,0, 4000, 4000); line.SetLineColor(2); + 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"); + line->Draw("same"); } TCanvas *cEC = new TCanvas("cEC" + suffix,Form("E corrected | %s", canvasTitle.Data()), 300 + 500 * aID, 300, canvasSize[0], canvasSize[1]); @@ -502,25 +511,7 @@ void MonPlotter::PlotEZ(){ void MonPlotter::PlotEx(){ - TCanvas *cex = new TCanvas("cex" + suffix,Form("EX : %s", canvasTitle.Data()),0, 0, 1000,650); - cex->Clear(); - gStyle->SetOptStat("neiou"); - hEx->Draw(""); - - TCanvas *cexI = new TCanvas("cexI" + suffix,Form("EX : %s", canvasTitle.Data()),500, 0, 1600,1000); - cexI->Clear();cexI->Divide(colDet,rowDet); - gStyle->SetOptStat("neiou"); - for( int i = 0; i < numDet; i++){ - cexI->cd(i+1); - hExi[i]->Draw(""); - } - - TCanvas *cExThetaCM = new TCanvas("cExThetaCM" + suffix,Form("EX - ThetaCM | %s", canvasTitle.Data()), 500, 500, 650,650); - cExThetaCM->Clear(); - gStyle->SetOptStat("neiou"); - hEx_ThetaCM->Draw("colz"); - - TCanvas *cExVxCal = new TCanvas("cExVxCal" + suffix,Form("EX | %s", canvasTitle.Data()),200, 200, 1600,1000); + 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); @@ -529,6 +520,48 @@ void MonPlotter::PlotEx(){ 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(""); + } } diff --git a/working/MonAnalyzer.C b/working/MonAnalyzer.C index 418ffd0..1996482 100644 --- a/working/MonAnalyzer.C +++ b/working/MonAnalyzer.C @@ -28,7 +28,7 @@ int rawEnergyRange[2] = { 0, 3000}; /// share with e, xf, xn int energyRange[2] = { 0, 10}; /// in the E-Z plot int rdtDERange[2] = { 0, 80}; int rdtERange[2] = { 0, 80}; -int thetaCMRange[2] = {0, 80}; +int thetaCMRange[2] = { 0, 50}; /// deg double exRange[3] = { 100, -2, 10}; /// bin [keV], low[MeV], high[MeV] @@ -51,6 +51,7 @@ std::vector rdtCutFile1 = {"", ""}; /// {reaction-0, reaction-1}, can a MonPlotter ** plotter = nullptr; int numGeo = 1; +TChain *gen_tree = nullptr; void MonAnalyzer(){ @@ -58,20 +59,20 @@ void MonAnalyzer(){ printf("####################### MonAnalyzer.C #######################\n"); printf("#####################################################################\n"); - TChain *chain = new TChain("gen_tree"); - //chain->Add("../root_data/gen_run043.root"); - chain->Add("../root_data/trace_run029.root"); + gen_tree = new TChain("gen_tree"); + //gen_tree->Add("../root_data/gen_run043.root"); + gen_tree->Add("../root_data/trace_run033.root"); - TObjArray * fileList = chain->GetListOfFiles(); + 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", chain->GetEntries()); + printf(" Total Number of entries : %llu \n", gen_tree->GetEntries()); printf("///////////////////////////////////////////////////////////////////\n"); - TTreeReader reader(chain); + TTreeReader reader(gen_tree); TTreeReaderValue evID = {reader, "evID"}; TTreeReaderArray e = {reader, "e"}; @@ -85,7 +86,7 @@ void MonAnalyzer(){ //TODO // TTreeReaderArray array = {reader, "trace"}; - ULong64_t NumEntries = chain->GetEntries(); + ULong64_t NumEntries = gen_tree->GetEntries(); //*========================================== DetGeo * detGeo = new DetGeo("detectorGeo.txt"); @@ -134,9 +135,13 @@ void MonAnalyzer(){ 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; @@ -148,6 +153,9 @@ void MonAnalyzer(){ 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); @@ -180,11 +188,12 @@ void MonAnalyzer(){ //@==================== 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 - if( corr->xnCorr.size() && corr->xfxneCorr.size() ) xnCal[id] = xn[id] * corr->xnCorr[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; - if( corr->xfxneCorr.size() ) xfCal[id] = xf[id] * corr->xfxneCorr[id][1] + corr->xfxneCorr[id][0]; - if( corr->eCorr.size() ) eCal[id] = e[id] / corr->eCorr[id][0] + corr->eCorr[id][1]; + 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; @@ -193,7 +202,7 @@ void MonAnalyzer(){ plotter[aID]->hxfCal_xnCal[id]->Fill(xfCal[id], xnCal[id]); plotter[aID]->he_xsCal[id]->Fill(xnCal[id] + xfCal[id], e[id]); - //@===================== calculate X + //@===================== 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; @@ -205,12 +214,12 @@ void MonAnalyzer(){ if ( TMath::IsNaN(xf[id]) && !TMath::IsNaN(xn[id]) ) xCal[id] = 1.0 - xnCal[id]/ e[id]; //@=================== Fill in histogram - plotter[aID]->he_x[id]->Fill(x[id],e[id]); + plotter[aID]->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) - if( corr->xScale.size() ) xCal[id] = (xCal[id]-0.5)/corr->xScale[id] + 0.5; /// if include this scale, need to also inclused in Cali_littleTree + 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; @@ -231,10 +240,10 @@ void MonAnalyzer(){ //@=================== Recoil Gate if( plotter[aID]->cutList ){ - for(int i = 0 ; i < cutList1->GetEntries() ; i++ ){ - TCutG * cutG = (TCutG *)cutList1->At(i) ; + 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; + rdtgate1 = true; break; /// only one is enough } } @@ -263,26 +272,27 @@ void MonAnalyzer(){ if( j%2 == 1) { plotter[aID]->htDiff->Fill(tdiff); // if((rdtgate1 || rdtgate2) && (eCalCut[1] > eCal[id] && eCal[id]>eCalCut[0])) { - // plotter[aID]->htdiffg->Fill(tdiff); - // } + 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 ) 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); + 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); - // hrdtg[j]->Fill(rdt[j]); - // coinFlag = true; + // plotter[aID]->hrdtg[j]->Fill(rdt[j]); + coinFlag = true; - // } + } } } - // if( !isTimeGateOn ) coinFlag = true; + if( !isTimeGateOn ) coinFlag = true; //@================ E-Z gate // if( EZCut ) { @@ -290,15 +300,14 @@ void MonAnalyzer(){ // }else{ // ezGate = true; // } - + + //@================ is Good Event ? // if( coinFlag && (rdtgate1 || rdtgate2) && ezGate){ - // plotter[arrayID]->heCalVzGC->Fill( z[id] , eCal[id] ); + if( coinFlag && rdtgate1){ + plotter[aID]->heCal_zGC->Fill( z[id] , eCal[id] ); - // heCalVxCalG[id]->Fill(xCal[id]*detGeo->array[arrayID].detLength, eCal[id]); - - // multiEZ ++; - // isGoodEventFlag = true; - // } + isGoodEventFlag = true; + } }//*====== end of array loop @@ -321,8 +330,9 @@ void MonAnalyzer(){ 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 ; @@ -363,10 +373,11 @@ void MonAnalyzer(){ 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; + // if( processedEntries > 1000 ) break; }//^############################################## End of Process gStyle->SetOptStat("neiou"); @@ -386,13 +397,13 @@ void MonAnalyzer(){ // printf("-----------------------------------------------------\n"); printf(" raw() - Raw data\n"); printf(" cal() - Calibrated data\n"); + printf(" rdt() - Raw RDT data\n"); printf("-----------------------------------------------------\n"); printf(" ez() - Energy vs. Z\n"); + printf("-----------------------------------------------------\n"); printf(" excite() - Excitation Energy\n"); - // printf(" recoils() - Raw DE vs. E Recoil spectra\n"); //printf(" elum() - Luminosity Energy Spectra\n"); //printf(" ic() - Ionization Chamber Spectra\n"); - // printf("-----------------------------------------------------\n"); // printf(" eCalVzRow() - Energy vs. Z for each row\n"); // printf(" ExThetaCM() - Ex vs ThetaCM\n"); // printf(" ExVxCal() - Ex vs X for all %d detectors\n", numDet); @@ -426,6 +437,10 @@ void cal(int arrayID = -1){ } } +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(); diff --git a/working/DataHoSei.C b/working/obsolete/DataHoSei.C similarity index 100% rename from working/DataHoSei.C rename to working/obsolete/DataHoSei.C diff --git a/working/DataHoSei.h b/working/obsolete/DataHoSei.h similarity index 100% rename from working/DataHoSei.h rename to working/obsolete/DataHoSei.h diff --git a/working/Monitor.C b/working/obsolete/Monitor.C similarity index 100% rename from working/Monitor.C rename to working/obsolete/Monitor.C diff --git a/working/Monitor.h b/working/obsolete/Monitor.h similarity index 100% rename from working/Monitor.h rename to working/obsolete/Monitor.h diff --git a/Armory/Monitor_Util.C b/working/obsolete/Monitor_Util.C similarity index 100% rename from Armory/Monitor_Util.C rename to working/obsolete/Monitor_Util.C diff --git a/working/script_multi.C b/working/obsolete/script_multi.C similarity index 100% rename from working/script_multi.C rename to working/obsolete/script_multi.C diff --git a/working/script_single.C b/working/obsolete/script_single.C similarity index 100% rename from working/script_single.C rename to working/obsolete/script_single.C