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

finished the PID alignment script.

Browse Source
This commit is contained in:
Ryan Tang 2022-06-22 19:33:20 -04:00
parent 6db36958ec
commit 7614580770
7 changed files with 737 additions and 344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
armory/AnalysisLibrary.h
View File

@ -91,6 +91,28 @@ double* sumMeanVar(std::vector<double> data){
return output; return output;
} }
double Rsquared(std::vector<double> x, std::vector<double> 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<double> dataX, std::vector<double> dataY){ double* fitSlopeIntercept(std::vector<double> dataX, std::vector<double> dataY){
double * smvY = sumMeanVar(dataY); double * smvY = sumMeanVar(dataY);

143
compareTH2D.C Normal file
View File

@ -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");
}

20
correction_PID.dat Normal file
View File

@ -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

790
peachCake.C
View File

@ -9,6 +9,13 @@
#include <TCutG.h> #include <TCutG.h>
#include "./armory/AnalysisLibrary.h" #include "./armory/AnalysisLibrary.h"
//============== User settings
int tofRange[3] = {500, -280, -210};
int dERange[3] = {500, 0, 5500};
//============== histograms //============== histograms
TH2F * hCloverID; TH2F * hCloverID;
@ -16,6 +23,10 @@ TH2F * hPID;
TH1F * hZ; TH1F * hZ;
TH2F * hPID0; // raw TH2F * hPID0; // raw
TH2F * hPID2; // z vs A/Q TH2F * hPID2; // z vs A/Q
TH1F * hTOF;
TH1F * hdE;
TH2F * hImplantIons; TH2F * hImplantIons;
TH2F * hImplantBeta; //high gain TH2F * hImplantBeta; //high gain
@ -31,6 +42,9 @@ TH1F * hDecay_veto;
TH1F * hFlag; TH1F * hFlag;
TH1F * hvetoFlag; TH1F * hvetoFlag;
TH2F * heVrun;
TH2F * htVrun;
//============ parameters //============ parameters
//TString cutFileName="PIDCuts.root"; //TString cutFileName="PIDCuts.root";
//TFile * cutFile; //TFile * cutFile;
@ -40,385 +54,464 @@ TH1F * hvetoFlag;
vector<vector<double>> eCorr; vector<vector<double>> 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 // for run 250
+0.0000710839 * energy*energy double par[6] = {10.9179,
-6.28402e-8 * TMath::Power(energy,3) 0.00416034,
+2.90563e-11 * TMath::Power(energy,4) -233.306,
-6.70137e-15 * TMath::Power(energy,5) 3607.16,
+6.08418e-19 * TMath::Power(energy,6) ); 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) );
//return 34./TMath::Sqrt(energy - 93) -235.5;
} }
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*/){ 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); hPID = createTH2F("hPID", "PID; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]);
hPID0 = new TH2F("hPID0", "PID raw; ns; msx100", 1000, -265, -220, 1000, 0, 3500); hPID0 = createTH2F("hPID0", "PID raw; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]);
hPID2 = new TH2F("hPID2", "PID2; A/Q; Z", 500, 2, 4.2, 500, 0, 16); 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); printf("============ start \n");
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");
} }
Bool_t peachCake::Process(Long64_t entry){ 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(); nProcessed++;
Long64_t startTimeL = 0;
Long64_t startTimeR = 0;
Long64_t stopTime = 0;
double cfdL = TMath::QuietNaN();
double cfdR = TMath::QuietNaN();
double cfd2 = TMath::QuietNaN();
Long64_t startTime40 = 0; b_multi->GetEntry(entry);
double cfd40 = TMath::QuietNaN(); 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 energy = TMath::QuietNaN();
double a_H[5] = {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(); Long64_t startTime40 = 0;
veto_r = TMath::QuietNaN(); double cfd40 = 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 double a_L[5] = {TMath::QuietNaN()}; ///0 = dy, 1-4 = anode
if( 0 <= detID[i] && detID[i] < 100 ){ double a_H[5] = {TMath::QuietNaN()};
if( e[i] > 0 ) {
int id = detID[i]; veto_f = TMath::QuietNaN();
double eCal = 0; veto_r = TMath::QuietNaN();
for( int k = 0; k < int(eCorr[id].size()); k++){
eCal += eCorr[id][k] * TMath::Power(e[i], k); veto_f_time = 0;
} veto_r_time = 0;
hCloverID->Fill(detID[i], eCal);
} crossEnergy = 0; /// for analog signal, cross T2
} crossTime = 0;
//----- dy Beta
if( detID[i] == 200 ) { for( int i = 0; i < 4 ; i++) {
if ( multiDyBeta < 4 ){ dyIonsTime[i] = 0;
dyBetaTime[multiDyBeta] = e_t[i]; dyBetaTime[i] = 0;
dyBetaEnergy[multiDyBeta] = e[i]; dyIonsEnergy[i] = 0;
multiDyBeta ++; 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);
} }
} hCloverID->Fill(detID[i], eCal);
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;
} }
} }
//----- 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 //================================== PID
if ( flag == 0 ) return kTRUE; 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(); //================================== Implant
newTree->Fill(); double sumA_L = 0;
} double sumA_H = 0;
for( int i = 1; i <= 4 ; i++){
return kTRUE; 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(){ void peachCake::Terminate(){
printf("\n===============\n"); printf("\n===============\n");
if( saveNewTree ){
saveFile->cd(); saveFile->cd();
newTree->Write(); newTree->Write();
saveFile->Close(); saveFile->Close();
}
gStyle->SetOptStat("neiou"); if( plotHists ){
TCanvas * c1 = new TCanvas("c1", "c1", 1400, 1400);
c1->Divide(2,2);
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 ++; int padID=1;
c1->cd(padID); //c1->cd(padID)->SetLogz();
hPID->Draw("colz");
padID ++; c1->cd(padID); //c1->cd(padID)->SetLogz();
c1->cd(padID); //c1->cd(padID)->SetLogz(); hPID0->Draw("colz");
hPID2->Draw("colz");
padID ++; //padID ++;
c1->cd(padID); //c1->cd(padID)->SetLogz(); //c1->cd(padID); //c1->cd(padID)->SetLogz();
hZ->Draw("colz"); //hTOF->Draw();
//
//padID ++; //padID ++;
//c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz(); //c1->cd(padID); //c1->cd(padID)->SetLogz();
//hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13); //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 ++; padID ++;
@ -447,4 +540,21 @@ void peachCake::Terminate(){
hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same"); 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());
}
} }

16
peachCake.h
View File

@ -49,6 +49,9 @@ public :
peachCake(TTree * /*tree*/ =0) : fChain(0) { peachCake(TTree * /*tree*/ =0) : fChain(0) {
saveNewTree = false; saveNewTree = false;
pidCorrFileName = "";
fHistRootName = "";
plotHists = true;
} }
virtual ~peachCake() { } virtual ~peachCake() { }
virtual Int_t Version() const { return 2; } virtual Int_t Version() const { return 2; }
@ -66,9 +69,17 @@ public :
virtual void Terminate(); virtual void Terminate();
void SaveNewTree(bool YesOrNo){ saveNewTree = YesOrNo;} 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); ClassDef(peachCake,0);
TString pidCorrFileName;
vector<vector<double>> pidCorr;
//=========== varibales in new tree //=========== varibales in new tree
bool saveNewTree; bool saveNewTree;
@ -107,6 +118,10 @@ public :
Long64_t nProcessed; Long64_t nProcessed;
TString fHistRootName;
bool plotHists;
}; };
#endif #endif
@ -193,6 +208,7 @@ void peachCake::Init(TTree *tree){
nProcessed = 0; nProcessed = 0;
printf("==========================================\n");
} }
Bool_t peachCake::Notify(){ Bool_t peachCake::Notify(){

42
plotPID.C Normal file
View File

@ -0,0 +1,42 @@
#include "peachCake.C+"
void plotPID() {
vector<int> 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;
}
}

48
script.C
View File

@ -5,30 +5,70 @@
#include "peachCake.C+" #include "peachCake.C+"
void script() { void script() {
int runNum = 250;
TChain * chain = new TChain("tree"); 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-0237-*.root");
//chain->Add("root_data/run-0238-*.root"); //chain->Add("root_data/run-0238-*.root");
//chain->Add("root_data/run-0241-*.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-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; bool isSaveNewTree = false;
TString histRootFileName = "";// Form("PID_%03d.root", runNum);
TString pidCorrFileName = "correction_PID.dat";
chain->GetListOfFiles()->Print(); chain->GetListOfFiles()->Print();
printf("================================\n");
peachCake * selector = new peachCake(); peachCake * selector = new peachCake();
selector->SaveNewTree(isSaveNewTree); selector->SaveNewTree(isSaveNewTree);
selector->SetPIDCorrectionFile(pidCorrFileName);
selector->SetHistRootFileName(histRootFileName);
chain->Process(selector, ""); chain->Process(selector, "");
///gROOT->ProcessLine("armory/nsclEvtReader.h"); ///gROOT->ProcessLine("armory/nsclEvtReader.h");
///evt->ReadBlock(2); ///evt->ReadBlock(2);
} }