From 76145807708b41364d9ad16c9a09e6a6bac026ec Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Wed, 22 Jun 2022 19:33:20 -0400 Subject: [PATCH] finished the PID alignment script. --- armory/AnalysisLibrary.h | 22 ++ compareTH2D.C | 143 +++++++ correction_PID.dat | 20 + peachCake.C | 790 ++++++++++++++++++++++----------------- peachCake.h | 16 + plotPID.C | 42 +++ script.C | 48 ++- 7 files changed, 737 insertions(+), 344 deletions(-) create mode 100644 compareTH2D.C create mode 100644 correction_PID.dat create mode 100644 plotPID.C diff --git a/armory/AnalysisLibrary.h b/armory/AnalysisLibrary.h index b4f7eb8..416e0fb 100644 --- a/armory/AnalysisLibrary.h +++ b/armory/AnalysisLibrary.h @@ -91,6 +91,28 @@ double* sumMeanVar(std::vector data){ return output; } +double Rsquared(std::vector x, std::vector y){ + + if ( x.size() != y.size() ) return 0; + + int n = x.size(); + double * qX = sumMeanVar(x); + + double meanX = qX[1]; + double varX = qX[2]; + + double * qY = sumMeanVar(y); + double meanY = qY[1]; + double varY = qY[2]; + + double sumXY = 0; + for( int k = 0; k < n; k++) sumXY += x[k] * y[k]; + + double r2 = (sumXY - n * meanX * meanY)/ sqrt( varX * varY ); + + return r2; +} + double* fitSlopeIntercept(std::vector dataX, std::vector dataY){ double * smvY = sumMeanVar(dataY); diff --git a/compareTH2D.C b/compareTH2D.C new file mode 100644 index 0000000..5af0ab1 --- /dev/null +++ b/compareTH2D.C @@ -0,0 +1,143 @@ +#include "armory/AnalysisLibrary.h" + +double xMin1, xMax1; +TGraph * fitGraph1 = new TGraph(); +Double_t fitFunc1(Double_t *x, Double_t *par) { + double pos = par[1]*x[0] + par[0]; + return (xMax1 > pos && pos > xMin1 ? par[2] * fitGraph1->Eval(pos): 0); +} + +TF1 * FitHist1(TH1F * hist1, TH1F * hist2, double a0, double a1, double scale){ + + if( hist1->GetNbinsX() != hist2->GetNbinsX() ) return NULL; + + xMin1 = hist1->GetXaxis()->GetXmin(); + xMax1 = hist2->GetXaxis()->GetXmax(); + int nBin = hist1->GetNbinsX(); + + for( int i = 0; i < nBin; i++){ + fitGraph1->AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); + } + + double par[3] = {a0, a1, scale}; + + TF1 * fit = new TF1("fit", fitFunc1, xMin1, xMax1, 3); + fit->SetParameters(par); + fit->SetLineWidth(2); + fit->SetLineColor(4); + fit->SetNpx(1000); + + hist2->Fit("fit", "Rq"); + + return fit; +} + +double xMin2, xMax2; +TGraph * fitGraph2 = new TGraph(); +Double_t fitFunc2(Double_t *x, Double_t *par) { + double pos = par[1]*x[0] + par[0]; + + //return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos) * exp(pos/1000./par[3]): 0); + return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos) : 0); + + + //if( x[0] < 2200 ){ + // return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos): 0); + //}else{ + // return (xMax2 > pos && pos > xMin2 ? par[3] * fitGraph2->Eval(pos): 0); + //} + + +} + +TF1 * FitHist2(TH1F * hist1, TH1F * hist2, double a0, double a1, double scale, double scale2){ + + if( hist1->GetNbinsX() != hist2->GetNbinsX() ) return NULL; + + xMin2 = hist1->GetXaxis()->GetXmin(); + xMax2 = hist2->GetXaxis()->GetXmax(); + int nBin = hist1->GetNbinsX(); + + for( int i = 0; i < nBin; i++){ + fitGraph2->AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); + } + + double par[4] = {a0, a1, scale, scale2}; + + TF1 * fit = new TF1("fit", fitFunc2, xMin2, xMax2, 4); + fit->SetParameters(par); + fit->SetLineWidth(2); + fit->SetLineColor(4); + fit->SetNpx(1000); + + hist2->Fit("fit", "Rq"); + + return fit; +} + + + +void compareTH2D(int run1, int run2){ + + TFile * f1 = new TFile(Form("PID_%d.root", run1)); + TH1F * hTOF1 = (TH1F*) f1->Get("hTOF"); + TH1F * hdE1 = (TH1F*) f1->Get("hdE"); + TH2F * hPID1 = (TH2F*) f1->Get("hPID0"); + + TFile * f2 = new TFile(Form("PID_%d.root", run2)); + TH1F * hTOF2 = (TH1F*) f2->Get("hTOF"); + TH1F * hdE2 = (TH1F*) f2->Get("hdE"); + TH2F * hPID2 = (TH2F*) f2->Get("hPID0"); + + TCanvas * cc = new TCanvas("cc", Form("cc %d, %d", run1, run2), 3000, 1000); + cc->Divide(3,1); + + //================= TOF + cc->cd(1); + hTOF1->Draw(""); + hTOF2->SetLineColor(2); + hTOF2->Draw("same"); + + TF1 * fitTOF = FitHist1(hTOF1, hTOF2, 6, 1.0, 0.01); + fitTOF->Draw("same"); + + //TPad * chaha = (TPad*) (cc->cd(1))->Clone(); + //cc->cd(); + //chaha->Draw(); + + const Double_t* paraE_TOF = fitTOF->GetParErrors(); + const Double_t* paraA_TOF = fitTOF->GetParameters(); + + printf("=================== TOF \n"); + printf("a0 : %7.4f +- %7.4f \n", paraA_TOF[0], paraE_TOF[0]); + printf("a1 : %7.4f +- %7.4f \n", paraA_TOF[1], paraE_TOF[1]); + printf(" B : %7.4f +- %7.4f \n", paraA_TOF[2], paraE_TOF[2]); + + + //============== dE + cc->cd(2); cc->cd(2)->SetLogy(); + hdE1->Draw(""); + hdE2->SetLineColor(2); + hdE2->Draw("same"); + + TF1 * fitdE = FitHist2(hdE1, hdE2, -1, 1, 1, -2); + + fitdE->Draw("same"); + + const Double_t* paraE_dE = fitdE->GetParErrors(); + const Double_t* paraA_dE = fitdE->GetParameters(); + + printf("=================== dE \n"); + printf("a0 : %7.4f +- %7.4f \n", paraA_dE[0], paraE_dE[0]); + printf("a1 : %7.4f +- %7.4f \n", paraA_dE[1], paraE_dE[1]); + printf(" B : %7.4f +- %7.4f \n", paraA_dE[2], paraE_dE[2]); + printf("B2 : %7.4f +- %7.4f \n", paraA_dE[3], paraE_dE[3]); + + printf("################################\n"); + printf("%3d %9.6f %9.6f %9.6f %9.6f\n", run2, paraA_TOF[0], paraA_TOF[1], paraA_dE[0], paraA_dE[1]); + printf("\n\n\n"); + + cc->cd(3); + hPID1->Draw("colz"); + hPID2->Draw("same"); +} diff --git a/correction_PID.dat b/correction_PID.dat new file mode 100644 index 0000000..8906f79 --- /dev/null +++ b/correction_PID.dat @@ -0,0 +1,20 @@ +246 5.202374 0.972581 1.823645 0.982412 +247 5.474717 0.973615 2.228469 0.989046 +248 -2.016251 0.974793 3.005200 0.993577 +249 -2.670544 0.972111 2.883073 0.996861 +250 0.000000 1.000000 0.000000 1.000000 +251 -0.657032 0.997220 -0.209729 1.001230 +252 -0.440989 0.998095 -0.235229 1.002683 +253 -0.736023 0.996942 -0.291993 1.003505 +254 8.379539 1.005064 -14.304124 1.011162 +256 12.176462 1.001000 10.480354 0.988767 +257 6.429354 1.010119 10.926955 0.995218 +258 8.853158 1.019825 11.033525 0.998783 +259 15.584682 1.014310 10.904478 1.000129 +261 14.648562 1.010594 9.730146 1.004049 +262 6.638140 1.010922 9.764069 1.005226 +263 7.778790 0.987777 9.657852 1.005952 +264 6.860067 1.011819 9.671465 1.005891 +269 4.648145 1.006698 9.352800 0.921318 +270 12.961422 1.008092 9.634760 0.892950 +271 5.361827 1.010096 9.085634 0.884655 diff --git a/peachCake.C b/peachCake.C index 0ed317c..3f322be 100644 --- a/peachCake.C +++ b/peachCake.C @@ -9,6 +9,13 @@ #include #include "./armory/AnalysisLibrary.h" + +//============== User settings + +int tofRange[3] = {500, -280, -210}; +int dERange[3] = {500, 0, 5500}; + + //============== histograms TH2F * hCloverID; @@ -16,6 +23,10 @@ TH2F * hPID; TH1F * hZ; TH2F * hPID0; // raw TH2F * hPID2; // z vs A/Q + +TH1F * hTOF; +TH1F * hdE; + TH2F * hImplantIons; TH2F * hImplantBeta; //high gain @@ -31,6 +42,9 @@ TH1F * hDecay_veto; TH1F * hFlag; TH1F * hvetoFlag; +TH2F * heVrun; +TH2F * htVrun; + //============ parameters //TString cutFileName="PIDCuts.root"; //TFile * cutFile; @@ -40,385 +54,464 @@ TH1F * hvetoFlag; vector> eCorr; -double TOFCorrection(double energy){ +double TOFCorrection(double x){ + //for run 238 + //return (-225.8 - 0.040017 * energy + // +0.0000710839 * energy*energy + // -6.28402e-8 * TMath::Power(energy,3) + // +2.90563e-11 * TMath::Power(energy,4) + // -6.70137e-15 * TMath::Power(energy,5) + // +6.08418e-19 * TMath::Power(energy,6) ); - return (-225.8 - 0.040017 * energy - +0.0000710839 * energy*energy - -6.28402e-8 * TMath::Power(energy,3) - +2.90563e-11 * TMath::Power(energy,4) - -6.70137e-15 * TMath::Power(energy,5) - +6.08418e-19 * TMath::Power(energy,6) ); - - - //return 34./TMath::Sqrt(energy - 93) -235.5; + // for run 250 + double par[6] = {10.9179, + 0.00416034, + -233.306, + 3607.16, + 1538.37, + 0.000430609}; + return par[0] * exp( - par[1] * x ) + par[2] - par[5] * (par[3] - x) * exp( - pow( (par[3]-x)/par[4],2) ); + } +TH2F * createTH2F(const char* name, const char* title, Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy, Double_t ylow, Double_t yup){ + + TH2F * hist2 = (TH2F *) gROOT->FindObjectAny( name ); + if ( hist2 == NULL ) hist2 = new TH2F( name , title , nbinsx, xlow, xup, nbinsy, ylow, yup); + hist2->Reset(); + + return hist2; +} + +TH1F * createTH1F(const char* name, const char* title, Int_t nbinsx, Double_t xlow, Double_t xup){ + + TH1F * hist1 = (TH1F *) gROOT->FindObjectAny( name ); + if ( hist1 == NULL ) hist1 = new TH1F( name , title , nbinsx, xlow, xup); + hist1->Reset(); + + return hist1; +} void peachCake::Begin(TTree * /*tree*/){ - TString option = GetOption(); + TString option = GetOption(); - hCloverID = new TH2F("hCloverID", "Clover; ID; energy [keV]", 52, 0, 52, 400, 0, 2000); + hCloverID = createTH2F("hCloverID", "Clover; ID; energy [keV]", 52, 0, 52, 400, 0, 2000); - hPID = new TH2F("hPID", "PID; ns; msx100", 1000, -265, -220, 1000, 0, 3500); - hPID0 = new TH2F("hPID0", "PID raw; ns; msx100", 1000, -265, -220, 1000, 0, 3500); - hPID2 = new TH2F("hPID2", "PID2; A/Q; Z", 500, 2, 4.2, 500, 0, 16); + hPID = createTH2F("hPID", "PID; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]); + hPID0 = createTH2F("hPID0", "PID raw; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]); + hPID2 = createTH2F("hPID2", "PID2; A/Q; Z", tofRange[0], 2.4, 3.6, dERange[0], 0, 16.5); + + hTOF = createTH1F("hTOF", "TOF", tofRange[0], tofRange[1], tofRange[2]); + hdE = createTH1F("hdE", "dE", dERange[0], dERange[1], dERange[2]); + + hZ = createTH1F("hZ", "Z; Z; count", 500, 0, 16); + + hImplantIons = createTH2F("hImplantIons", "Implat low gain; x ; y", 400, 0, 1, 400, 0, 1); + hImplantBeta = createTH2F("hImplantBeta", "Implat high gain; x ; y", 400, 0, 1, 400, 0, 1); + + hImplantX = createTH2F("hImplantX", "Implat X; low ; high", 400, 0, 1, 400, 0, 1); + hImplantY = createTH2F("hImplantY", "Implat Y; low ; high", 400, 0, 1, 400, 0, 1); + + + hImplantDyIons = createTH2F("hImplantDyIonsow", "Implat low; sum_a; dy", 400, 0, 80000, 400, 0, 80000); + hImplantDyBeta = createTH2F("hImplantDyBetaigh", "Implat high; sum_a; dy", 400, 0, 8000, 400, 0, 8000); + + hDecay = createTH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000); + hDecay_veto = createTH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000); + + hFlag = createTH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8); + hvetoFlag = createTH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8); + + heVrun = createTH2F("heVrun", "energy vs run; run; energy", 200, 237, 257, 400, 0, 2000); + htVrun = createTH2F("hrVrun", "TOF vs run; run; energy", 200, 237, 257, 400, -265, -220); + + eCorr = LoadCorrectionParameters("correction_e.dat"); + + if( pidCorrFileName != "" ){ + pidCorr = LoadCorrectionParameters(pidCorrFileName, 1); + } + + /** + cutFile = new TFile(cutFileName, "READ"); + bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); + if( listExist ) { + cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); + numCut = cutList->GetLast()+1; + printf("----- found %d cuts \n", numCut); + for( int k = 0; k < numCut; k++){ + if( cutList->At(k) != NULL ){ + printf("found a cut at %2d \n", k); + } + } + } + * */ - hZ = new TH1F("hZ", "Z; dE; count", 500, 0, 16); - - hImplantIons = new TH2F("hImplantIons", "Implat low gain; x ; y", 400, 0, 1, 400, 0, 1); - hImplantBeta = new TH2F("hImplantBeta", "Implat high gain; x ; y", 400, 0, 1, 400, 0, 1); - - hImplantX = new TH2F("hImplantX", "Implat X; low ; high", 400, 0, 1, 400, 0, 1); - hImplantY = new TH2F("hImplantY", "Implat Y; low ; high", 400, 0, 1, 400, 0, 1); - - - hImplantDyIons = new TH2F("hImplantDyIonsow", "Implat low; sum_a; dy", 400, 0, 80000, 400, 0, 80000); - hImplantDyBeta = new TH2F("hImplantDyBetaigh", "Implat high; sum_a; dy", 400, 0, 8000, 400, 0, 8000); - - hDecay = new TH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000); - hDecay_veto = new TH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000); - - hFlag = new TH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8); - hvetoFlag = new TH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8); - - eCorr = LoadCorrectionParameters("correction_e.dat"); - - /** - cutFile = new TFile(cutFileName, "READ"); - bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); - if( listExist ) { - cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); - numCut = cutList->GetLast()+1; - printf("----- found %d cuts \n", numCut); - for( int k = 0; k < numCut; k++){ - if( cutList->At(k) != NULL ){ - printf("found a cut at %2d \n", k); - } - } - } - * */ - - printf("============ start \n"); + printf("============ start \n"); } Bool_t peachCake::Process(Long64_t entry){ - - nProcessed++; - - b_multi->GetEntry(entry); - b_detID->GetEntry(entry); - b_e->GetEntry(entry); - b_e_t->GetEntry(entry); - b_cfd->GetEntry(entry); - energy = TMath::QuietNaN(); - Long64_t startTimeL = 0; - Long64_t startTimeR = 0; - Long64_t stopTime = 0; - double cfdL = TMath::QuietNaN(); - double cfdR = TMath::QuietNaN(); - double cfd2 = TMath::QuietNaN(); + nProcessed++; - Long64_t startTime40 = 0; - double cfd40 = TMath::QuietNaN(); + b_multi->GetEntry(entry); + b_detID->GetEntry(entry); + b_e->GetEntry(entry); + b_e_t->GetEntry(entry); + b_cfd->GetEntry(entry); + b_runID->GetEntry(entry); - double a_L[5] = {TMath::QuietNaN()}; ///0 = dy, 1-4 = anode - double a_H[5] = {TMath::QuietNaN()}; + energy = TMath::QuietNaN(); + Long64_t startTimeL = 0; + Long64_t startTimeR = 0; + Long64_t stopTime = 0; + double cfdL = TMath::QuietNaN(); + double cfdR = TMath::QuietNaN(); + double cfd2 = TMath::QuietNaN(); - veto_f = TMath::QuietNaN(); - veto_r = TMath::QuietNaN(); - - veto_f_time = 0; - veto_r_time = 0; - - crossEnergy = 0; /// for analog signal, cross T2 - crossTime = 0; - - for( int i = 0; i < 4 ; i++) { - dyIonsTime[i] = 0; - dyBetaTime[i] = 0; - dyIonsEnergy[i] = 0; - dyBetaEnergy[i] = 0; - } - - int multiDyBeta = 0; - int multiDyIons = 0; - - flag = 0; - vetoFlag = 0; - - //======= Scanning the event - for( int i = 0; i < multi; i++){ + Long64_t startTime40 = 0; + double cfd40 = TMath::QuietNaN(); - //----- clover - if( 0 <= detID[i] && detID[i] < 100 ){ - if( e[i] > 0 ) { - int id = detID[i]; - double eCal = 0; - for( int k = 0; k < int(eCorr[id].size()); k++){ - eCal += eCorr[id][k] * TMath::Power(e[i], k); - } - hCloverID->Fill(detID[i], eCal); - } - } - //----- dy Beta - if( detID[i] == 200 ) { - if ( multiDyBeta < 4 ){ - dyBetaTime[multiDyBeta] = e_t[i]; - dyBetaEnergy[multiDyBeta] = e[i]; - multiDyBeta ++; + double a_L[5] = {TMath::QuietNaN()}; ///0 = dy, 1-4 = anode + double a_H[5] = {TMath::QuietNaN()}; + + veto_f = TMath::QuietNaN(); + veto_r = TMath::QuietNaN(); + + veto_f_time = 0; + veto_r_time = 0; + + crossEnergy = 0; /// for analog signal, cross T2 + crossTime = 0; + + for( int i = 0; i < 4 ; i++) { + dyIonsTime[i] = 0; + dyBetaTime[i] = 0; + dyIonsEnergy[i] = 0; + dyBetaEnergy[i] = 0; + } + + int multiDyBeta = 0; + int multiDyIons = 0; + + flag = 0; + vetoFlag = 0; + + //======= Scanning the event + for( int i = 0; i < multi; i++){ + //----- clover + if( 0 <= detID[i] && detID[i] < 100 ){ + if( e[i] > 0 ) { + int id = detID[i]; + double eCal = 0; + for( int k = 0; k < int(eCorr[id].size()); k++){ + eCal += eCorr[id][k] * TMath::Power(e[i], k); } - } - if( 201 <= detID[i] && detID[i] < 250){ - a_H[detID[i] - 200] = e[i]; - } - //----- dy Ions - if( detID[i] == 250 ) { - if ( multiDyIons < 4 ){ - dyIonsTime[multiDyIons] = e_t[i]; - dyIonsEnergy[multiDyIons] = e[i]; - multiDyIons ++; - } - } - if( 251 <= detID[i] && detID[i] < 260){ - a_L[detID[i] - 250] = e[i]; - } - //----- veto_front - if( detID[i] == 290 ) { - veto_f = e[i]; - veto_f_time = e_t[i]; - if( veto_f > 0 && (vetoFlag & 1) == 0) vetoFlag += 1; - } - //----- veto_rear - if( detID[i] == 291 ) { - veto_r = e[i]; - veto_r_time = e_t[i]; - if( veto_r > 0 && (vetoFlag & 2) == 0) vetoFlag += 2; - } - //----- image scint L - if( detID[i] == 300 ) { - if( e[i] > 1000 ) { - startTimeL = e_t[i]; - cfdL = cfd[i]/16384.; - } - } - //----- image scint R - if( detID[i] == 301 ) { - startTimeR = e_t[i]; - cfdR = cfd[i]/16384.; - } - //----- cross T2 analog - if( detID[i] == 306 ) { - if( (vetoFlag & 4) == 0) vetoFlag += 4; - crossEnergy = e[i]; - crossTime = e_t[i]; - } - //----- cross B2 logic - if( detID[i] == 307 ) { - stopTime = e_t[i]; - cfd2 = cfd[i]/16384.; - if( (vetoFlag & 4) == 0) vetoFlag += 4; - } - //----- MSX40 - ///if( detID[i] == 308 ) energy = e[i]; - //----- MSX100 - if( detID[i] == 309 ) energy = e[i]; - //----- MSX40 logic - if( detID[i] == 310 ) { - startTime40 = e_t[i]; - cfd40 = cfd[i]/16384.; - } - } - - //================================== PID - Long64_t startTime = startTimeL; - double startCFD = cfdL; - - TOF = TMath::QuietNaN(); - if( startTime > 0 && stopTime > 0 && energy > 0 ){ - if( stopTime > startTime ){ - TOF = double(stopTime - startTime) - startCFD + cfd2; - }else{ - TOF = -double(startTime - stopTime) - startCFD + cfd2 ; - } - TOF = TOF * 8; // to ns - - /// TOF correction, A/Q = 3 to have a fixed TOF. - hPID0->Fill(TOF, energy); - - double newTOF = TOF - TOFCorrection(energy) - 234.; - - hPID->Fill(newTOF, energy); - flag += 1; - - ///======== Z vs A/Q - double c = 299.792458; /// mm/ns - double FL = 56423 ; /// mm - double Brho = 9.04; /// T.mm - double ma = 931.5; - double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3 - double tofOffset = 504; /// ns - - double beta = FL/c/(newTOF+tofOffset); - double gamma = 1./TMath::Sqrt(1.-beta*beta); - - double Z = sqrt((energy + 19.63) / 15.14) * TMath::Power(beta / beta0,1.3) ; - double AoQ = c * Brho / gamma/ beta / ma; - - //printf("tof : %f, beta: %f (%f), Z : %f, A/Q : %f \n", TOF + tofOffset, beta, beta0, Z, AoQ); - - Z = -0.15938 + 1.01736 *Z + 0.000203316 * Z*Z; - - hPID2->Fill(AoQ, Z); - - hZ->Fill(Z); - - } - - //================================== Implant - double sumA_L = 0; - double sumA_H = 0; - for( int i = 1; i <= 4 ; i++){ - if( a_L[i] > 0 ) sumA_L += a_L[i]; - if( a_H[i] > 0 ) sumA_H += a_H[i]; - } - - xIons = (a_L[3]+a_L[2])/sumA_L; - yIons = (a_L[1]+a_L[2])/sumA_L; - xBeta = (a_H[3]+a_H[2])/sumA_H; - yBeta = (a_H[1]+a_H[2])/sumA_H; - - - if( (flag & 1) == 0 ) { - hImplantIons->Fill(xIons, yIons); - hImplantBeta->Fill(xBeta, yBeta); - hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); - hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); - - hImplantX->Fill(xIons, xBeta); - hImplantY->Fill(yIons, yBeta); - } - - if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){ - //hImplantIons->Fill(xIons, yIons); - //hImplantBeta->Fill(xBeta, yBeta); - //hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); - //hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); - // - //hImplantX->Fill(xIons, xBeta); - //hImplantY->Fill(yIons, yBeta); - } - - if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2; - if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4; - - hFlag->Fill(flag); - hvetoFlag->Fill(vetoFlag); - - - //============================== debug - if ( false ) { - - printf("################# %lld, multi:%d\n", entry, multi); - - printf("----- beam Line \n"); - printf(" MSX100 eneergy : %f \n", energy); - printf(" startTime : %llu, CFD : %f \n", startTime, startCFD); - printf(" stopTime : %llu, CFD : %f \n", stopTime, cfd2); - printf(" TOF : %f ns \n", TOF); - - printf("----- Ions, low gain \n"); - for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_L[i]); - printf("sum A : %f \n", sumA_L); - printf("x : %f , y : %f \n", xIons, yIons); - for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyIonsEnergy[i], dyIonsTime[i]); - - printf("----- Beta, high gain \n"); - for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_H[i]); - printf("sum A : %f \n", sumA_H); - printf("x : %f , y : %f \n", xBeta, yBeta); - for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyBetaEnergy[i], dyBetaTime[i]); - - printf("----- Veto \n"); - printf(" cross Time : %llu \n", crossTime); - printf(" cross energy : %d \n", crossEnergy); - printf("veto_f energy : %f \n", veto_f); - printf( "veto_f Time : %llu \n", veto_f_time); - printf("veto_r energy : %f \n", veto_r); - printf( "veto_r Time : %llu \n", veto_r_time); - - printf("============ flag : %d, vetoFlag : %d\n", flag, vetoFlag); - - } - - - //============================= Progress - clock.Stop("timer"); - Double_t time = clock.GetRealTime("timer"); - clock.Start("timer"); - - if ( !shown ) { - if (fmod(time, 10) < 1 ){ - printf( "%12llu[%2d%%]|%3.0f min %5.2f sec | expect:%5.2f min \n", - nProcessed, - TMath::Nint((nProcessed+1)*100./totnumEntry), - TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60., - totnumEntry*time/(nProcessed+1.)/60.); - shown = 1;} - }else{ - if (fmod(time, 10) > 9 ){ - shown = 0; + hCloverID->Fill(detID[i], eCal); } - } + } + //----- dy Beta + if( detID[i] == 200 ) { + if ( multiDyBeta < 4 ){ + dyBetaTime[multiDyBeta] = e_t[i]; + dyBetaEnergy[multiDyBeta] = e[i]; + multiDyBeta ++; + } + } + if( 201 <= detID[i] && detID[i] < 250){ + a_H[detID[i] - 200] = e[i]; + } + //----- dy Ions + if( detID[i] == 250 ) { + if ( multiDyIons < 4 ){ + dyIonsTime[multiDyIons] = e_t[i]; + dyIonsEnergy[multiDyIons] = e[i]; + multiDyIons ++; + } + } + if( 251 <= detID[i] && detID[i] < 260){ + a_L[detID[i] - 250] = e[i]; + } + //----- veto_front + if( detID[i] == 290 ) { + veto_f = e[i]; + veto_f_time = e_t[i]; + if( veto_f > 0 && (vetoFlag & 1) == 0) vetoFlag += 1; + } + //----- veto_rear + if( detID[i] == 291 ) { + veto_r = e[i]; + veto_r_time = e_t[i]; + if( veto_r > 0 && (vetoFlag & 2) == 0) vetoFlag += 2; + } + //----- image scint L + if( detID[i] == 300 ) { + if( e[i] > 1000 ) { + startTimeL = e_t[i]; + cfdL = cfd[i]/16384.; + } + } + //----- image scint R + if( detID[i] == 301 ) { + startTimeR = e_t[i]; + cfdR = cfd[i]/16384.; + } + //----- cross T2 analog + if( detID[i] == 306 ) { + if( (vetoFlag & 4) == 0) vetoFlag += 4; + crossEnergy = e[i]; + crossTime = e_t[i]; + } + //----- cross B2 logic + if( detID[i] == 307 ) { + stopTime = e_t[i]; + cfd2 = cfd[i]/16384.; + if( (vetoFlag & 4) == 0) vetoFlag += 4; + } + //----- MSX40 + ///if( detID[i] == 308 ) energy = e[i]; + //----- MSX100 + if( detID[i] == 309 ) energy = e[i]; + //----- MSX40 logic + if( detID[i] == 310 ) { + startTime40 = e_t[i]; + cfd40 = cfd[i]/16384.; + } + } - //============================= Rejection - if ( flag == 0 ) return kTRUE; + //================================== PID + Long64_t startTime = startTimeL; + double startCFD = cfdL; + + TOF = TMath::QuietNaN(); + if( startTime > 0 && stopTime > 0 && energy > 0 ){ + if( stopTime > startTime ){ + TOF = double(stopTime - startTime) - startCFD + cfd2; + }else{ + TOF = -double(startTime - stopTime) - startCFD + cfd2 ; + } + TOF = TOF * 8; // to ns + + ///============== PID correction + if( pidCorrFileName != "" ){ + for( int k = 0; k < (int) pidCorr.size() ; k++){ + if( int(pidCorr[k][0]) != runID/100 ) continue; + + TOF = pidCorr[k][1] + pidCorr[k][2] * TOF; + energy = pidCorr[k][3] + pidCorr[k][4] * energy; + + } + } + + + /// TOF correction, A/Q = 3 to have a fixed TOF. + + if( energy > 160./44.*(TOF + 260) +180 ) { + hPID0->Fill(TOF, energy); + hTOF->Fill(TOF); + hdE->Fill(energy); + } + + double newTOF = TOF - TOFCorrection(energy) - 240.; + + hPID->Fill(newTOF, energy); + + heVrun->Fill(runID/100., energy); + htVrun->Fill(runID/100., TOF); + + flag += 1; + + ///======== Z vs A/Q + double c = 299.792458; /// mm/ns + double FL = 37656 ; /// mm + double Brho = 4.90370; /// T.mm + double ma = 931.5; + double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3 + double tofOffset = 509.79; /// ns + + double beta = FL/c/(newTOF+tofOffset); + double gamma = 1./TMath::Sqrt(1.-beta*beta); + + double Z = sqrt((energy + 11.9473) / 23.097) * TMath::Power(beta / beta0,1.3) ; + double AoQ = c * Brho / gamma/ beta / ma; + + //printf("tof : %f, beta: %f (%f), Z : %f, A/Q : %f \n", TOF + tofOffset, beta, beta0, Z, AoQ); + + //Z = -0.15938 + 1.01736 *Z + 0.000203316 * Z*Z; + //Z = -0.254632 + 1.06285 * Z - 0.00539634 * Z * Z + 0.000169443 * Z * Z * Z; ///2022-06-09 + + hPID2->Fill(AoQ, Z); + + hZ->Fill(Z); - //============================= Save Tree - if( saveNewTree ){ - saveFile->cd(); - newTree->Fill(); - } - - return kTRUE; + } + + //================================== Implant + double sumA_L = 0; + double sumA_H = 0; + for( int i = 1; i <= 4 ; i++){ + if( a_L[i] > 0 ) sumA_L += a_L[i]; + if( a_H[i] > 0 ) sumA_H += a_H[i]; + } + + xIons = (a_L[3]+a_L[2])/sumA_L; + yIons = (a_L[1]+a_L[2])/sumA_L; + xBeta = (a_H[3]+a_H[2])/sumA_H; + yBeta = (a_H[1]+a_H[2])/sumA_H; + + + if( (flag & 1) == 0 ) { + hImplantIons->Fill(xIons, yIons); + hImplantBeta->Fill(xBeta, yBeta); + hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); + hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); + + hImplantX->Fill(xIons, xBeta); + hImplantY->Fill(yIons, yBeta); + } + + if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){ + //hImplantIons->Fill(xIons, yIons); + //hImplantBeta->Fill(xBeta, yBeta); + //hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); + //hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); + // + //hImplantX->Fill(xIons, xBeta); + //hImplantY->Fill(yIons, yBeta); + } + + if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2; + if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4; + + hFlag->Fill(flag); + hvetoFlag->Fill(vetoFlag); + + + //============================== debug + if ( false ) { + + printf("################# %lld, multi:%d\n", entry, multi); + + printf("----- beam Line \n"); + printf(" MSX100 eneergy : %f \n", energy); + printf(" startTime : %llu, CFD : %f \n", startTime, startCFD); + printf(" stopTime : %llu, CFD : %f \n", stopTime, cfd2); + printf(" TOF : %f ns \n", TOF); + + printf("----- Ions, low gain \n"); + for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_L[i]); + printf("sum A : %f \n", sumA_L); + printf("x : %f , y : %f \n", xIons, yIons); + for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyIonsEnergy[i], dyIonsTime[i]); + + printf("----- Beta, high gain \n"); + for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_H[i]); + printf("sum A : %f \n", sumA_H); + printf("x : %f , y : %f \n", xBeta, yBeta); + for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyBetaEnergy[i], dyBetaTime[i]); + + printf("----- Veto \n"); + printf(" cross Time : %llu \n", crossTime); + printf(" cross energy : %d \n", crossEnergy); + printf("veto_f energy : %f \n", veto_f); + printf( "veto_f Time : %llu \n", veto_f_time); + printf("veto_r energy : %f \n", veto_r); + printf( "veto_r Time : %llu \n", veto_r_time); + + printf("============ flag : %d, vetoFlag : %d\n", flag, vetoFlag); + + } + + + //============================= Progress + clock.Stop("timer"); + Double_t time = clock.GetRealTime("timer"); + clock.Start("timer"); + + if ( !shown ) { + if (fmod(time, 10) < 1 ){ + printf( "%12llu[%2d%%]|%3.0f min %5.2f sec | expect:%5.2f min \n", + nProcessed, + TMath::Nint((nProcessed+1)*100./totnumEntry), + TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60., + totnumEntry*time/(nProcessed+1.)/60.); + shown = 1;} + }else{ + if (fmod(time, 10) > 9 ){ + shown = 0; + } + } + + //============================= Rejection + if ( flag == 0 ) return kTRUE; + + //============================= Save Tree + if( saveNewTree ){ + saveFile->cd(); + newTree->Fill(); + } + + return kTRUE; } void peachCake::Terminate(){ printf("\n===============\n"); + + if( saveNewTree ){ saveFile->cd(); newTree->Write(); saveFile->Close(); + } - gStyle->SetOptStat("neiou"); - TCanvas * c1 = new TCanvas("c1", "c1", 1400, 1400); - c1->Divide(2,2); + if( plotHists ){ - int padID=1; + gStyle->SetOptStat(111111); + int div[2] = {3, 2} ; ///x, y + TCanvas * c1 = new TCanvas("c1", "c1", 700 * div[0], 700 * div[1]); + c1->Divide(div[0], div[1]); - c1->cd(padID); //c1->cd(padID)->SetLogz(); - hPID0->Draw("colz"); - padID ++; - c1->cd(padID); //c1->cd(padID)->SetLogz(); - hPID->Draw("colz"); + int padID=1; - padID ++; - c1->cd(padID); //c1->cd(padID)->SetLogz(); - hPID2->Draw("colz"); + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hPID0->Draw("colz"); - padID ++; - c1->cd(padID); //c1->cd(padID)->SetLogz(); - hZ->Draw("colz"); - - //padID ++; - //c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz(); - //hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13); + //padID ++; + //c1->cd(padID); //c1->cd(padID)->SetLogz(); + //hTOF->Draw(); + // + //padID ++; + //c1->cd(padID); //c1->cd(padID)->SetLogz(); + //hdE->Draw(); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hPID->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hPID2->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hZ->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + heVrun->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + htVrun->Draw("colz"); + + //padID ++; + //c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz(); + //hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13); /* padID ++; @@ -447,4 +540,21 @@ void peachCake::Terminate(){ hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same"); */ + + } + + //================ Save historgram + if( fHistRootName != "" ){ + TFile * fHist = new TFile(fHistRootName, "recreate"); + fHist->cd(); + + hTOF->Write(); + hdE->Write(); + hPID0->Write(); + + fHist->Close(); + + printf("---- Save PID histogram as %s \n", fHistRootName.Data()); + } + } diff --git a/peachCake.h b/peachCake.h index bab3200..03fc0ca 100644 --- a/peachCake.h +++ b/peachCake.h @@ -49,6 +49,9 @@ public : peachCake(TTree * /*tree*/ =0) : fChain(0) { saveNewTree = false; + pidCorrFileName = ""; + fHistRootName = ""; + plotHists = true; } virtual ~peachCake() { } virtual Int_t Version() const { return 2; } @@ -66,9 +69,17 @@ public : virtual void Terminate(); void SaveNewTree(bool YesOrNo){ saveNewTree = YesOrNo;} + void SetPIDCorrectionFile(TString corr_PID){ pidCorrFileName = corr_PID;} + + void SetHistRootFileName(TString fileName){ fHistRootName = fileName;} + + void SetPlotHist(bool onOff) { plotHists = onOff; } ClassDef(peachCake,0); + TString pidCorrFileName; + vector> pidCorr; + //=========== varibales in new tree bool saveNewTree; @@ -107,6 +118,10 @@ public : Long64_t nProcessed; + TString fHistRootName; + + bool plotHists; + }; #endif @@ -193,6 +208,7 @@ void peachCake::Init(TTree *tree){ nProcessed = 0; + printf("==========================================\n"); } Bool_t peachCake::Notify(){ diff --git a/plotPID.C b/plotPID.C new file mode 100644 index 0000000..89a196f --- /dev/null +++ b/plotPID.C @@ -0,0 +1,42 @@ + +#include "peachCake.C+" + +void plotPID() { + + vector skipRun = {0};//{239, 240, 242, 243, 244, 245, 255}; + + peachCake * selector = new peachCake(); + selector->SaveNewTree(false); + selector->SetPIDCorrectionFile(""); + selector->SetPlotHist(false); + + for( int run = 269 ; run < 272 ; run ++){ + + printf("\033[31m######################################### run : %d \033[0m\n", run); + + bool contFlag = false; + for( int k = 0; k < (int) skipRun.size(); k ++){ + if( run == skipRun[k] ) contFlag = true; + } + + if( contFlag ) continue; + + TChain * chain = new TChain("tree"); + chain->Add(Form("root_data/run-%04d-*.root", run)); + chain->GetListOfFiles()->Print(); + + int nFile = chain->GetListOfFiles()->GetEntries(); + + printf("================================ num of Files : %d \n", nFile); + + if( nFile == 0 ) continue; + + selector->SetHistRootFileName(Form("PID_%03d.root", run)); + + chain->Process(selector, ""); + + delete chain; + + } + +} diff --git a/script.C b/script.C index 7681756..47ccf92 100644 --- a/script.C +++ b/script.C @@ -5,30 +5,70 @@ #include "peachCake.C+" void script() { + + int runNum = 250; TChain * chain = new TChain("tree"); + //chain->Add(Form("root_data/run-%04d-*.root", runNum)); //chain->Add("root_data/run-0237-*.root"); //chain->Add("root_data/run-0238-*.root"); //chain->Add("root_data/run-0241-*.root"); - //chain->Add("root_data/run-0244-*.root"); + + // new beam //chain->Add("root_data/run-0246-*.root"); + //chain->Add("root_data/run-0247-*.root"); + //chain->Add("root_data/run-0248-*.root"); + //chain->Add("root_data/run-0249-*.root"); + //chain->Add("root_data/run-0250-*.root"); + //chain->Add("root_data/run-0251-*.root"); + //chain->Add("root_data/run-0252-*.root"); + //chain->Add("root_data/run-0253-*.root"); + //chain->Add("root_data/run-0254-*.root"); + //chain->Add("root_data/run-0257-*.root"); + //chain->Add("root_data/run-0258-*.root"); + //chain->Add("root_data/run-0259-*.root"); + //chain->Add("root_data/run-0261-*.root"); + //chain->Add("root_data/run-0262-*.root"); + //chain->Add("root_data/run-0263-*.root"); + //chain->Add("root_data/run-0264-*.root"); + //chain->Add("root_data/run-0269-*.root"); + //chain->Add("root_data/run-0270-*.root"); + chain->Add("root_data/run-0271-*.root"); + int nFile = chain->GetListOfFiles()->GetEntries(); + + printf("================================ num of Files : %d \n", nFile); - chain->Add("root_data/run-02[3-4][0-9]-00.root"); + if( nFile == 0 ) return; bool isSaveNewTree = false; - + + TString histRootFileName = "";// Form("PID_%03d.root", runNum); + + TString pidCorrFileName = "correction_PID.dat"; chain->GetListOfFiles()->Print(); + + printf("================================\n"); + peachCake * selector = new peachCake(); selector->SaveNewTree(isSaveNewTree); + selector->SetPIDCorrectionFile(pidCorrFileName); + selector->SetHistRootFileName(histRootFileName); + chain->Process(selector, ""); + ///gROOT->ProcessLine("armory/nsclEvtReader.h"); ///evt->ReadBlock(2); - + + + + + + }