From 54c9de513febcbe809f82eeb10e8313fa686cacb Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Sat, 22 Jan 2022 20:30:39 -0500 Subject: [PATCH] reclaibrate energy, add back in Analyzer --- Analyzer.C | 162 ++++++++++++++++++++++++++++++---------- Analyzer.h | 23 ++++-- armory/Analyzer_Utili.c | 44 +++++------ armory/to2root.cpp | 5 +- correction_e.dat | 73 +++++++++--------- 5 files changed, 198 insertions(+), 109 deletions(-) diff --git a/Analyzer.C b/Analyzer.C index 1a19163..6842a0a 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -12,26 +13,26 @@ //############################################ User setting int rawEnergyRange[2] = {0, 6000}; // in ch -int energyRange[3] = {1, 50, 2000}; // keV {resol, min, max} +int energyRange[3] = {1, 50, 6000}; // keV {resol, min, max} double BGO_threshold = 0; // in ch TString e_corr = "correction_e.dat"; -bool save_ev2 = true; +TString cutFileName = "proton209Cut.root"; + +bool save_ev2 = false; //############################################ end of user setting //############################################ histogram declaration - TH2F * heVID; TH1F * he[NCRYSTAL]; TH2F * hgg[NCRYSTAL][NCRYSTAL]; TH2F * hcoin; - TH2F * hcrystalBGO; TH1F * hTDiff; @@ -44,8 +45,17 @@ TH2F * hPID219; TH2F * heCalVID; TH1F * heCal[NCRYSTAL]; TH2F * hcoinBGO; - TH2F * hcrystalBGO_G; +TH1F * hg[NCLOVER]; + +///----- after particle gate +TH1F * heCal_g[NCRYSTAL]; +TH2F * heCalVID_g; + +TH1F * hg_g[NCLOVER]; + +///============= cut +TCutG * cut; //############################################ BEGIN void Analyzer::Begin(TTree * tree){ @@ -61,8 +71,9 @@ void Analyzer::Begin(TTree * tree){ printf("======================== Histograms declaration\n"); - heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); - heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); + heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + heCalVID_g = new TH2F("heCalVID_g", Form("eCal vs ID (BGO veto > %.1f & particle); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); hTDiff = new TH1F("hTDiff", "data time different within an event; tick [10 ns]", 110, 0, 110); @@ -70,13 +81,19 @@ void Analyzer::Begin(TTree * tree){ heCalVID->SetNdivisions(-409, "X"); for( int i = 0; i < NCRYSTAL; i ++){ - he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); - heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal -%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); + heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal-%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + heCal_g[i] = new TH1F(Form("heCal_g%02d", i), Form("eCal-%02d (BGO veto > %.1f & particle); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); } - hPID = new TH2F("hPID", "PID; tail; peak ", 400, -20, 300, 400, -50, 800); - hPID209 = new TH2F("hPID209", "PID detID = 209; tail; peak ", 400, -20, 300, 400, -50, 800); - hPID219 = new TH2F("hPID219", "PID detID = 219; tail; peak ", 400, -20, 300, 400, -50, 800); + for( int i = 0; i < NCLOVER; i++){ + hg[i] = new TH1F(Form("hg%02d", i), Form("Clover-%02d (added-back)", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + hg_g[i] = new TH1F(Form("hg_g%02d", i), Form("Clover-%02d (added-back) particle", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + } + + hPID = new TH2F("hPID", "PID; tail; peak ", 400, -20, 300, 400, -50, 1200); + hPID209 = new TH2F("hPID209", "PID detID = 209; tail; peak ", 400, -20, 300, 400, -50, 1200); + hPID219 = new TH2F("hPID219", "PID detID = 219; tail; peak ", 400, -20, 300, 400, -50, 1200); for( int i = 0; i < NCRYSTAL; i++){ @@ -97,6 +114,15 @@ void Analyzer::Begin(TTree * tree){ eCorr = LoadCorrectionParameters(e_corr); + + if( cutFileName != ""){ + printf("======================== Load cuts.\n"); + + TFile * cutFile = new TFile(cutFileName); + cut = (TCutG *) cutFile->Get("protonCut"); + printf(" %s is loaded.\n", cut->GetName()); + } + saveEV2 = save_ev2; } @@ -126,12 +152,47 @@ Bool_t Analyzer::Process(Long64_t entry){ if( multi == 0 ) return kTRUE; for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); + for( int i = 0; i < NCLOVER; i++) gamma[i] = 0; ///printf("---------------------------- %d \n", multi); double bgC[2]={0}, peakC[2]={0}, tailC[2] = {0}; int count = 0 ; + ///=========== make the particle gate + for( int i = 0; i < multi ; i ++){ + int id = detID[i]; + + if( id == 209 ) { + double bg = (qdc[i][0] + qdc[i][1])/60.; + double peak = qdc[i][3]/20. - bg; + double tail = qdc[i][5]/55. - bg; + hPID209->Fill( tail, peak); + } + if( id == 219 ) { + double bg = (qdc[i][0] + qdc[i][1])/60.; + double peak = qdc[i][3]/20. - bg; + double tail = qdc[i][5]/55. - bg; + hPID219->Fill( tail, peak); + } + + if( count < 2 && (id == 209 || id == 219 )){ + bgC[count] = (qdc[i][0] + qdc[i][1])/60.; + peakC[count] = qdc[i][3]/20. - bgC[count]; + tailC[count] = qdc[i][5]/55. - bgC[count]; + count++; + } + + } + + ///######################################################### + ///================ coincident gate between proton and gamma + double tail = (tailC[0]+tailC[1])/2.; + double peak = (peakC[0]+peakC[1])/2.; + + hPID->Fill( tail, peak); + + ///=========== Looping data for the event for( int i = 0; i < multi ; i ++){ @@ -160,25 +221,7 @@ Bool_t Analyzer::Process(Long64_t entry){ } - if( id == 209 ) { - double bg = (qdc[i][0] + qdc[i][1])/60.; - double peak = qdc[i][3]/20. - bg; - double tail = qdc[i][5]/55. - bg; - hPID209->Fill( tail, peak); - } - if( id == 219 ) { - double bg = (qdc[i][0] + qdc[i][1])/60.; - double peak = qdc[i][3]/20. - bg; - double tail = qdc[i][5]/55. - bg; - hPID219->Fill( tail, peak); - } - if( count < 2 && (id == 209 || id == 219 )){ - bgC[count] = (qdc[i][0] + qdc[i][1])/60.; - peakC[count] = qdc[i][3]/20. - bgC[count]; - tailC[count] = qdc[i][5]/55. - bgC[count]; - count++; - } if ( i > 0 ) hTDiff->Fill( e_t[i] - e_t[0]); @@ -195,7 +238,6 @@ Bool_t Analyzer::Process(Long64_t entry){ } if( dropflag ) return kTRUE; - if( 0<= id && id < NCRYSTAL ) { if( e_corr == "" ){ eCal[id] = e[i]; @@ -213,14 +255,33 @@ Bool_t Analyzer::Process(Long64_t entry){ for ( int j = i + 1; j < multi; j++){ if( 100 <= detID[j] && detID[j] < 200 ) hcrystalBGO_G->Fill(id, detID[j]-100); /// crystal - BGO coincident + } + + + ///========== add back and remove cross talk + int cloverID = id /4; + + if( eCal[id] > 100 ) gamma[cloverID] += eCal[id]; + + ///====== particle coincidet + if( cut->IsInside(tail, peak) ){ + heCal_g[id]->Fill(eCal[id]); + heCalVID_g->Fill(id, eCal[id]); } } } - - hPID->Fill( (tailC[0]+tailC[1])/2., ( peakC[0] + peakC[1])/2.); + for( int i = 0 ; i < NCLOVER; i++){ + if( gamma[i] > 0 ) { + + hg[i]->Fill(gamma[i]); + + if( cut->IsInside(tail, peak) ) hg_g[i]->Fill(gamma[i]); + + } + } if ( saveEV2) Save2ev2(); @@ -246,25 +307,44 @@ void Analyzer::Terminate(){ cCanvas->cd(1); cCanvas->cd(1)->SetLogz(1); - hPID209->Draw("colz"); + heCalVID->Draw("colz"); cCanvas->cd(2); cCanvas->cd(2)->SetLogz(1); - hPID219->Draw("colz"); + hPID->Draw("colz"); + cut->Draw("same"); cCanvas->cd(3); cCanvas->cd(3)->SetLogz(1); - hPID->Draw("colz"); + heCalVID_g->Draw("colz"); cCanvas->cd(4); - cCanvas->cd(4)->SetLogz(1); + cCanvas->cd(4)->SetLogy(1); //gROOT->ProcessLine(".x script.C"); //hcrystalBGO_G->Draw("colz"); - //he[0]->SetLineColor(2); - //he[0]->Draw(); - //heCal[0]->Draw("same"); + hg[0]->SetLineColor(2); + hg[0]->Draw(); + hg_g[0]->Draw("same"); //hcoinBGO->Draw("colz"); - + + TCanvas *cAux = new TCanvas("cAux", "" ,1000, 0, 800,600); + cAux->cd(); + + TH1F * h0 = (TH1F *) heCal[0]->Clone("h0"); + h0->Add(heCal[1]); + h0->Add(heCal[2]); + h0->Add(heCal[3]); + + TH1F * h0_g = (TH1F *) heCal_g[0]->Clone("h0_g"); + h0_g->Add(heCal_g[1]); + h0_g->Add(heCal_g[2]); + h0_g->Add(heCal_g[3]); + + h0->SetLineColor(kGreen+3); + h0->Draw(); + h0_g->Draw("same"); + hg[0]->Draw("same"); + printf("=============== loaded AutoFit.C, try showFitMethos()\n"); gROOT->ProcessLine(".L armory/AutoFit.C"); printf("=============== Analyzer Utility\n"); diff --git a/Analyzer.h b/Analyzer.h index a3a03bd..f7300d3 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -30,21 +30,25 @@ public : // Declaration of leaf types ULong64_t evID; + Int_t runID; Int_t detID[MAX_MULTI]; Double_t e[MAX_MULTI]; ULong64_t e_t[MAX_MULTI]; Int_t qdc[MAX_MULTI][8]; Int_t multi; Int_t multiCry; + Int_t multiGagg; // List of branches TBranch *b_event_ID; //! + TBranch *b_runID; //! TBranch *b_detID; //! TBranch *b_energy; //! TBranch *b_time; //! TBranch *b_qdc; //! TBranch *b_multi; //! TBranch *b_multiCry; //! + TBranch *b_multiGagg; //! Analyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0; Frac = 0.1; @@ -81,7 +85,10 @@ public : bool saveEV2; FILE * outEV2; TString outEV2Name; + + double eCal[NCRYSTAL]; + double gamma[NCLOVER]; // added-back energy }; #endif @@ -102,13 +109,15 @@ void Analyzer::Init(TTree *tree) fChain = tree; fChain->SetMakeClass(1); - fChain->SetBranchAddress("evID", &evID, &b_event_ID); - fChain->SetBranchAddress("detID", detID, &b_detID); - fChain->SetBranchAddress("e", e, &b_energy); - fChain->SetBranchAddress("e_t", e_t, &b_time); - fChain->SetBranchAddress("qdc", qdc, &b_qdc); - fChain->SetBranchAddress("multi", &multi, &b_multi); - fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry); + fChain->SetBranchAddress("evID", &evID, &b_event_ID); + fChain->SetBranchAddress("runID", &runID, &b_runID); + fChain->SetBranchAddress("detID", detID, &b_detID); + fChain->SetBranchAddress("e", e, &b_energy); + fChain->SetBranchAddress("e_t", e_t, &b_time); + fChain->SetBranchAddress("qdc", qdc, &b_qdc); + fChain->SetBranchAddress("multi", &multi, &b_multi); + fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry); + fChain->SetBranchAddress("multiGagg", &multiGagg, &b_multiGagg); TString option = GetOption(); if ( option != "" ) outEV2Name = option; diff --git a/armory/Analyzer_Utili.c b/armory/Analyzer_Utili.c index cc3f21c..4549f22 100644 --- a/armory/Analyzer_Utili.c +++ b/armory/Analyzer_Utili.c @@ -77,27 +77,27 @@ void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMi } - /* - for( Int_t j = 0; j < nCrystalPerClover; j++){ - int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1; - cRawE->cd(canvasID); - cRawE->cd(canvasID)->SetGrid(); - cRawE->cd(canvasID)->SetTickx(2); - cRawE->cd(canvasID)->SetTicky(2); - cRawE->cd(canvasID)->SetBottomMargin(0.06); - if ( i == nClover -1 ) cRawE->cd(canvasID)->SetRightMargin(0.002); - if( isLogy ) cRawE->cd(canvasID)->SetLogy(); - int hID = CloverID < 0 ? nCrystalPerClover*i+ j : nCrystalPerClover * CloverID + j ; - if( cali ) { - if ( xMin != 0 || xMax != 0 ) heCal[hID]->GetXaxis()->SetRangeUser(xMin, xMax); - heCal[hID]->SetMaximum(maxY); - heCal[hID]->Draw(""); - }else{ - if ( xMin != 0 || xMax != 0 ) he[hID]->GetXaxis()->SetRangeUser(xMin, xMax); - he[hID]->SetMaximum(maxY); - he[hID]->Draw(""); - } - }*/ + + ///for( Int_t j = 0; j < nCrystalPerClover; j++){ + /// int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1; + /// cRawE->cd(canvasID); + /// cRawE->cd(canvasID)->SetGrid(); + /// cRawE->cd(canvasID)->SetTickx(2); + /// cRawE->cd(canvasID)->SetTicky(2); + /// cRawE->cd(canvasID)->SetBottomMargin(0.06); + /// if ( i == nClover -1 ) cRawE->cd(canvasID)->SetRightMargin(0.002); + /// if( isLogy ) cRawE->cd(canvasID)->SetLogy(); + /// int hID = CloverID < 0 ? nCrystalPerClover*i+ j : nCrystalPerClover * CloverID + j ; + /// if( cali ) { + /// if ( xMin != 0 || xMax != 0 ) heCal[hID]->GetXaxis()->SetRangeUser(xMin, xMax); + /// heCal[hID]->SetMaximum(maxY); + /// heCal[hID]->Draw(""); + /// }else{ + /// if ( xMin != 0 || xMax != 0 ) he[hID]->GetXaxis()->SetRangeUser(xMin, xMax); + /// he[hID]->SetMaximum(maxY); + /// he[hID]->Draw(""); + /// } + ///} } cRawE->SetCrosshair(1); @@ -144,7 +144,7 @@ void energyCalibration(int detID = -1, int BG = 10, double threshold = 0.1, doub 867.390, 964.055, 1085.842, - 1089.700, + ///1089.700, 1112.087, 1408.022}; diff --git a/armory/to2root.cpp b/armory/to2root.cpp index 7612ec0..30ef780 100644 --- a/armory/to2root.cpp +++ b/armory/to2root.cpp @@ -139,14 +139,13 @@ int main(int argc, char **argv) { //Check for end of event, rewind, and break out of while loop if (tdif > timeWindow) { - //Gate + Gate if( multiCry > 0 && multiGagg > 0 ) { outRootFile->cd(); tree->Fill(); - + countGP++; - } evID ++; diff --git a/correction_e.dat b/correction_e.dat index 3f491ab..2f34252 100644 --- a/correction_e.dat +++ b/correction_e.dat @@ -1,36 +1,37 @@ - -0.11851105 0.73602844 - -0.05329699 0.30617141 - 0.21919925 0.30256406 - 0.19549630 0.29988911 - 0.23576185 0.31232252 - 0.15862498 0.31719203 - 0.17456490 0.30495670 - 0.07937138 0.31393499 - -0.17752085 0.30734768 - 0.47543250 0.30864178 - 0.08375230 0.30843852 - 0.28037171 0.31263324 - -0.04410183 0.74143159 - 0.07905962 0.73641543 - -0.05892825 0.71786084 - 0.07476386 0.36785529 - 0.07951184 0.26260823 - 0.02161385 0.25884364 - 0.25371149 0.29681333 - 0.23290589 0.34255969 - 0.20677949 0.30504413 - 0.16341696 0.30761409 - 0.04406586 0.30595347 - 0.07292338 0.30758425 - 0.00136881 0.17925931 - 0.23758200 0.31520725 - -0.47281914 0.17676788 - 0.04230014 0.17917457 - 0.20654489 0.30340233 - 0.20762807 0.30960594 - 0.19673688 0.30110502 - 0.07362825 0.29715807 - 0.17023147 0.30259114 - 0.23642061 0.29846387 - 0.15627111 0.31411607 - -0.01255732 0.15930220 + 0.2198244789 0.3679598044 + 0.1820861672 0.1519049184 + 0.2826720086 0.1495011052 + 0.2276197020 0.1487333470 + 0.2701041460 0.1548212887 + 0.2373135736 0.1572502143 + 0.1611211260 0.1511348420 + 0.0763625382 0.1554911646 + -0.0311746466 0.1515947085 + 0.7135668017 0.1529901094 + 0.2757900177 0.1520869728 + 0.5408119147 0.1549171095 + 0.3290766322 0.3706024656 + 0.3681932261 0.3689368402 + 0.3604623220 0.3603678504 + 0.3857242086 0.1838693331 + 0.2175140528 0.1314356363 + 0.0874480331 0.1290573470 + 0.3928400202 0.1479762423 + 0.2514077943 0.1707848626 + 0.1375653253 0.1509230900 + 0.1447086122 0.1531835401 + 0.1784763570 0.1519575184 + 0.1400936263 0.1518167665 + 0.0915977450 0.0904661792 + 0.2911458476 0.1577982901 + 98.4907657327 0.0816230149 + 0.2636618573 0.0901480589 + -0.3874033379 0.1507092013 + 0.2204322169 0.1533291004 + 0.4647959240 0.1508393248 + 0.1381054743 0.1485480928 + 0.2238173432 0.1502901467 + 0.2951614679 0.1481891679 + 0.2229254902 0.1551385739 + 0.0663532207 0.0800650041 +