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

reclaibrate energy, add back in Analyzer

Browse Source
This commit is contained in:
Ryan Tang 2022-01-22 20:30:39 -05:00
parent e281afa456
commit 54c9de513f
5 changed files with 198 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
Analyzer.C
View File

@ -4,6 +4,7 @@
#include <TH2.h> #include <TH2.h>
#include <TStyle.h> #include <TStyle.h>
#include <TH1.h> #include <TH1.h>
#include <TCutG.h>
#include <TCanvas.h> #include <TCanvas.h>
#include <TMath.h> #include <TMath.h>
#include <vector> #include <vector>
@ -12,26 +13,26 @@
//############################################ User setting //############################################ User setting
int rawEnergyRange[2] = {0, 6000}; // in ch 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 double BGO_threshold = 0; // in ch
TString e_corr = "correction_e.dat"; TString e_corr = "correction_e.dat";
bool save_ev2 = true; TString cutFileName = "proton209Cut.root";
bool save_ev2 = false;
//############################################ end of user setting //############################################ end of user setting
//############################################ histogram declaration //############################################ histogram declaration
TH2F * heVID; TH2F * heVID;
TH1F * he[NCRYSTAL]; TH1F * he[NCRYSTAL];
TH2F * hgg[NCRYSTAL][NCRYSTAL]; TH2F * hgg[NCRYSTAL][NCRYSTAL];
TH2F * hcoin; TH2F * hcoin;
TH2F * hcrystalBGO; TH2F * hcrystalBGO;
TH1F * hTDiff; TH1F * hTDiff;
@ -44,8 +45,17 @@ TH2F * hPID219;
TH2F * heCalVID; TH2F * heCalVID;
TH1F * heCal[NCRYSTAL]; TH1F * heCal[NCRYSTAL];
TH2F * hcoinBGO; TH2F * hcoinBGO;
TH2F * hcrystalBGO_G; TH2F * hcrystalBGO_G;
TH1F * hg[NCLOVER];
///----- after particle gate
TH1F * heCal_g[NCRYSTAL];
TH2F * heCalVID_g;
TH1F * hg_g[NCLOVER];
///============= cut
TCutG * cut;
//############################################ BEGIN //############################################ BEGIN
void Analyzer::Begin(TTree * tree){ void Analyzer::Begin(TTree * tree){
@ -63,6 +73,7 @@ void Analyzer::Begin(TTree * tree){
heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); 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 = 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); hTDiff = new TH1F("hTDiff", "data time different within an event; tick [10 ns]", 110, 0, 110);
@ -72,11 +83,17 @@ void Analyzer::Begin(TTree * tree){
for( int i = 0; i < NCRYSTAL; i ++){ 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]); 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[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); for( int i = 0; i < NCLOVER; i++){
hPID209 = new TH2F("hPID209", "PID detID = 209; tail; peak ", 400, -20, 300, 400, -50, 800); hg[i] = new TH1F(Form("hg%02d", i), Form("Clover-%02d (added-back)", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
hPID219 = new TH2F("hPID219", "PID detID = 219; tail; peak ", 400, -20, 300, 400, -50, 800); 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++){ for( int i = 0; i < NCRYSTAL; i++){
@ -97,6 +114,15 @@ void Analyzer::Begin(TTree * tree){
eCorr = LoadCorrectionParameters(e_corr); 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; saveEV2 = save_ev2;
} }
@ -126,12 +152,47 @@ Bool_t Analyzer::Process(Long64_t entry){
if( multi == 0 ) return kTRUE; if( multi == 0 ) return kTRUE;
for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); 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); ///printf("---------------------------- %d \n", multi);
double bgC[2]={0}, peakC[2]={0}, tailC[2] = {0}; double bgC[2]={0}, peakC[2]={0}, tailC[2] = {0};
int count = 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 ///=========== Looping data for the event
for( int i = 0; i < multi ; i ++){ 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]); 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( dropflag ) return kTRUE;
if( 0<= id && id < NCRYSTAL ) { if( 0<= id && id < NCRYSTAL ) {
if( e_corr == "" ){ if( e_corr == "" ){
eCal[id] = e[i]; eCal[id] = e[i];
@ -215,12 +257,31 @@ Bool_t Analyzer::Process(Long64_t entry){
if( 100 <= detID[j] && detID[j] < 200 ) hcrystalBGO_G->Fill(id, detID[j]-100); /// crystal - BGO coincident 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(); if ( saveEV2) Save2ev2();
@ -246,25 +307,44 @@ void Analyzer::Terminate(){
cCanvas->cd(1); cCanvas->cd(1);
cCanvas->cd(1)->SetLogz(1); cCanvas->cd(1)->SetLogz(1);
hPID209->Draw("colz"); heCalVID->Draw("colz");
cCanvas->cd(2); cCanvas->cd(2);
cCanvas->cd(2)->SetLogz(1); cCanvas->cd(2)->SetLogz(1);
hPID219->Draw("colz"); hPID->Draw("colz");
cut->Draw("same");
cCanvas->cd(3); cCanvas->cd(3);
cCanvas->cd(3)->SetLogz(1); cCanvas->cd(3)->SetLogz(1);
hPID->Draw("colz"); heCalVID_g->Draw("colz");
cCanvas->cd(4); cCanvas->cd(4);
cCanvas->cd(4)->SetLogz(1); cCanvas->cd(4)->SetLogy(1);
//gROOT->ProcessLine(".x script.C"); //gROOT->ProcessLine(".x script.C");
//hcrystalBGO_G->Draw("colz"); //hcrystalBGO_G->Draw("colz");
//he[0]->SetLineColor(2); hg[0]->SetLineColor(2);
//he[0]->Draw(); hg[0]->Draw();
//heCal[0]->Draw("same"); hg_g[0]->Draw("same");
//hcoinBGO->Draw("colz"); //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"); printf("=============== loaded AutoFit.C, try showFitMethos()\n");
gROOT->ProcessLine(".L armory/AutoFit.C"); gROOT->ProcessLine(".L armory/AutoFit.C");
printf("=============== Analyzer Utility\n"); printf("=============== Analyzer Utility\n");

9
Analyzer.h
View File

@ -30,21 +30,25 @@ public :
// Declaration of leaf types // Declaration of leaf types
ULong64_t evID; ULong64_t evID;
Int_t runID;
Int_t detID[MAX_MULTI]; Int_t detID[MAX_MULTI];
Double_t e[MAX_MULTI]; Double_t e[MAX_MULTI];
ULong64_t e_t[MAX_MULTI]; ULong64_t e_t[MAX_MULTI];
Int_t qdc[MAX_MULTI][8]; Int_t qdc[MAX_MULTI][8];
Int_t multi; Int_t multi;
Int_t multiCry; Int_t multiCry;
Int_t multiGagg;
// List of branches // List of branches
TBranch *b_event_ID; //! TBranch *b_event_ID; //!
TBranch *b_runID; //!
TBranch *b_detID; //! TBranch *b_detID; //!
TBranch *b_energy; //! TBranch *b_energy; //!
TBranch *b_time; //! TBranch *b_time; //!
TBranch *b_qdc; //! TBranch *b_qdc; //!
TBranch *b_multi; //! TBranch *b_multi; //!
TBranch *b_multiCry; //! TBranch *b_multiCry; //!
TBranch *b_multiGagg; //!
Analyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0; Analyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0;
Frac = 0.1; Frac = 0.1;
@ -81,7 +85,10 @@ public :
bool saveEV2; bool saveEV2;
FILE * outEV2; FILE * outEV2;
TString outEV2Name; TString outEV2Name;
double eCal[NCRYSTAL]; double eCal[NCRYSTAL];
double gamma[NCLOVER]; // added-back energy
}; };
#endif #endif
@ -103,12 +110,14 @@ void Analyzer::Init(TTree *tree)
fChain->SetMakeClass(1); fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_event_ID); fChain->SetBranchAddress("evID", &evID, &b_event_ID);
fChain->SetBranchAddress("runID", &runID, &b_runID);
fChain->SetBranchAddress("detID", detID, &b_detID); fChain->SetBranchAddress("detID", detID, &b_detID);
fChain->SetBranchAddress("e", e, &b_energy); fChain->SetBranchAddress("e", e, &b_energy);
fChain->SetBranchAddress("e_t", e_t, &b_time); fChain->SetBranchAddress("e_t", e_t, &b_time);
fChain->SetBranchAddress("qdc", qdc, &b_qdc); fChain->SetBranchAddress("qdc", qdc, &b_qdc);
fChain->SetBranchAddress("multi", &multi, &b_multi); fChain->SetBranchAddress("multi", &multi, &b_multi);
fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry); fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry);
fChain->SetBranchAddress("multiGagg", &multiGagg, &b_multiGagg);
TString option = GetOption(); TString option = GetOption();
if ( option != "" ) outEV2Name = option; if ( option != "" ) outEV2Name = option;

44
armory/Analyzer_Utili.c
View File

@ -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++){ ///for( Int_t j = 0; j < nCrystalPerClover; j++){
int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1; /// int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1;
cRawE->cd(canvasID); /// cRawE->cd(canvasID);
cRawE->cd(canvasID)->SetGrid(); /// cRawE->cd(canvasID)->SetGrid();
cRawE->cd(canvasID)->SetTickx(2); /// cRawE->cd(canvasID)->SetTickx(2);
cRawE->cd(canvasID)->SetTicky(2); /// cRawE->cd(canvasID)->SetTicky(2);
cRawE->cd(canvasID)->SetBottomMargin(0.06); /// cRawE->cd(canvasID)->SetBottomMargin(0.06);
if ( i == nClover -1 ) cRawE->cd(canvasID)->SetRightMargin(0.002); /// if ( i == nClover -1 ) cRawE->cd(canvasID)->SetRightMargin(0.002);
if( isLogy ) cRawE->cd(canvasID)->SetLogy(); /// if( isLogy ) cRawE->cd(canvasID)->SetLogy();
int hID = CloverID < 0 ? nCrystalPerClover*i+ j : nCrystalPerClover * CloverID + j ; /// int hID = CloverID < 0 ? nCrystalPerClover*i+ j : nCrystalPerClover * CloverID + j ;
if( cali ) { /// if( cali ) {
if ( xMin != 0 || xMax != 0 ) heCal[hID]->GetXaxis()->SetRangeUser(xMin, xMax); /// if ( xMin != 0 || xMax != 0 ) heCal[hID]->GetXaxis()->SetRangeUser(xMin, xMax);
heCal[hID]->SetMaximum(maxY); /// heCal[hID]->SetMaximum(maxY);
heCal[hID]->Draw(""); /// heCal[hID]->Draw("");
}else{ /// }else{
if ( xMin != 0 || xMax != 0 ) he[hID]->GetXaxis()->SetRangeUser(xMin, xMax); /// if ( xMin != 0 || xMax != 0 ) he[hID]->GetXaxis()->SetRangeUser(xMin, xMax);
he[hID]->SetMaximum(maxY); /// he[hID]->SetMaximum(maxY);
he[hID]->Draw(""); /// he[hID]->Draw("");
} /// }
}*/ ///}
} }
cRawE->SetCrosshair(1); cRawE->SetCrosshair(1);
@ -144,7 +144,7 @@ void energyCalibration(int detID = -1, int BG = 10, double threshold = 0.1, doub
867.390, 867.390,
964.055, 964.055,
1085.842, 1085.842,
1089.700, ///1089.700,
1112.087, 1112.087,
1408.022}; 1408.022};

3
armory/to2root.cpp
View File

@ -139,14 +139,13 @@ int main(int argc, char **argv) {
//Check for end of event, rewind, and break out of while loop //Check for end of event, rewind, and break out of while loop
if (tdif > timeWindow) { if (tdif > timeWindow) {
//Gate Gate
if( multiCry > 0 && multiGagg > 0 ) { if( multiCry > 0 && multiGagg > 0 ) {
outRootFile->cd(); outRootFile->cd();
tree->Fill(); tree->Fill();
countGP++; countGP++;
} }
evID ++; evID ++;

73
correction_e.dat
View File

@ -1,36 +1,37 @@
-0.11851105 0.73602844 0.2198244789 0.3679598044
-0.05329699 0.30617141 0.1820861672 0.1519049184
0.21919925 0.30256406 0.2826720086 0.1495011052
0.19549630 0.29988911 0.2276197020 0.1487333470
0.23576185 0.31232252 0.2701041460 0.1548212887
0.15862498 0.31719203 0.2373135736 0.1572502143
0.17456490 0.30495670 0.1611211260 0.1511348420
0.07937138 0.31393499 0.0763625382 0.1554911646
-0.17752085 0.30734768 -0.0311746466 0.1515947085
0.47543250 0.30864178 0.7135668017 0.1529901094
0.08375230 0.30843852 0.2757900177 0.1520869728
0.28037171 0.31263324 0.5408119147 0.1549171095
-0.04410183 0.74143159 0.3290766322 0.3706024656
0.07905962 0.73641543 0.3681932261 0.3689368402
-0.05892825 0.71786084 0.3604623220 0.3603678504
0.07476386 0.36785529 0.3857242086 0.1838693331
0.07951184 0.26260823 0.2175140528 0.1314356363
0.02161385 0.25884364 0.0874480331 0.1290573470
0.25371149 0.29681333 0.3928400202 0.1479762423
0.23290589 0.34255969 0.2514077943 0.1707848626
0.20677949 0.30504413 0.1375653253 0.1509230900
0.16341696 0.30761409 0.1447086122 0.1531835401
0.04406586 0.30595347 0.1784763570 0.1519575184
0.07292338 0.30758425 0.1400936263 0.1518167665
0.00136881 0.17925931 0.0915977450 0.0904661792
0.23758200 0.31520725 0.2911458476 0.1577982901
-0.47281914 0.17676788 98.4907657327 0.0816230149
0.04230014 0.17917457 0.2636618573 0.0901480589
0.20654489 0.30340233 -0.3874033379 0.1507092013
0.20762807 0.30960594 0.2204322169 0.1533291004
0.19673688 0.30110502 0.4647959240 0.1508393248
0.07362825 0.29715807 0.1381054743 0.1485480928
0.17023147 0.30259114 0.2238173432 0.1502901467
0.23642061 0.29846387 0.2951614679 0.1481891679
0.15627111 0.31411607 0.2229254902 0.1551385739
-0.01255732 0.15930220 0.0663532207 0.0800650041