#include "TimeRecoTool.h"#include "Event/CalibReadout.h"#include "Event/CalibReadoutPmtCrate.h"#include "Event/CalibReadoutPmtChannel.h"#include "Event/RecTrigger.h"#include "Event/UserDataHeader.h"#include "Conventions/Reconstruction.h"#include "Conventions/Detectors.h"#include "DataSvc/ICableSvc.h"#include "DataSvc/ICalibDataSvc.h"#include "DybKernel/IPmtCalibSvc.h"#include "DetHelpers/IPmtGeomInfoSvc.h"#include "CLHEP/Vector/LorentzVector.h"#include "DybKernel/IChannelQualitySvc.h"#include <iostream>#include <fstream>#include <math.h>#include "TROOT.h"#include "TStyle.h"#include "TFile.h"#include "TChain.h"#include "TTree.h"#include "TF1.h"#include "TH2.h"#include "TMath.h"#include "TCanvas.h"#include "TMinuit.h"
Include dependency graph for TimeRecoTool.cc:
Go to the source code of this file.
Defines | |
| #define | NOCQ_ERRORS 10 |
Functions | |
| Double_t | chi2_sum_t (Double_t x, Double_t y, Double_t z, Double_t t, Int_t iSelect) |
| Double_t | chi2_time_t (Double_t x, Double_t y, Double_t z, Double_t t, Int_t iRing, Int_t iColumn) |
| void | IterativeFit (Int_t n, Int_t maxstep, Double_t x, Double_t y, Double_t z, Int_t iSelect) |
| void | SimpleStart () |
| Double_t | SimpleR (Int_t ring, Int_t column) |
| void | CalculateDeltime (Double_t x, Double_t y, Double_t z) |
| void | SelectPMT () |
Variables | |
| const Int_t | maxstep_t = 20 |
| const Int_t | maxselect_t = 5 |
| Int_t | nPMT_t [maxselect_t] |
| Double_t | z_phys_min_t = -2400 |
| Double_t | z_phys_max_t = 2400 |
| Double_t | x_phys_max_t = 2300 |
| Double_t | x_phys_min_t = -2300 |
| Double_t | y_phys_max_t = 2300 |
| Double_t | y_phys_min_t = -2300 |
| Double_t | r_phys_max_t = 2300 |
| Double_t | total_pe_t |
| Double_t | early_deltime_t |
| Double_t | median_deltime_t |
| Double_t | deltime_t [8][24] |
| Double_t | pmt_charge_t [8][24] |
| Double_t | pmt_time_t [8][24] |
| Double_t | pmt_x_t [8][24] |
| Double_t | pmt_y_t [8][24] |
| Double_t | pmt_z_t [8][24] |
| Double_t | pmt_error_t [8][24] |
| Int_t | pmt_flag_t [8][24] |
| Double_t | initial_pos_t [4] |
| Double_t | initial_energy_t |
| Int_t | totalsteps |
| Double_t | final_pos_t [4] |
| Int_t | entry |
Define Documentation
| #define NOCQ_ERRORS 10 |
Definition at line 34 of file TimeRecoTool.cc.
Function Documentation
| Double_t chi2_sum_t | ( | Double_t | x, |
| Double_t | y, | ||
| Double_t | z, | ||
| Double_t | t, | ||
| Int_t | iSelect | ||
| ) |
Definition at line 652 of file TimeRecoTool.cc.
{
Double_t chi2_tmp = 0;
Int_t nHits = 0;
for (Int_t iring = 0; iring < 8; iring++){
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
//This is the time portion
if (pmt_flag_t[iring][icolumn] != -1 && pmt_flag_t[iring][icolumn] <= iSelect){
chi2_tmp += chi2_time_t(x, y, z, t, iring, icolumn);
nHits ++;
}
}
}
//std::cout << "Sum Chi2: " << chi2_tmp << " nHits: " << nHits << std::endl;
if (nHits == 0) chi2_tmp = 1e10;
else chi2_tmp = chi2_tmp / nHits;
return chi2_tmp;
}
| Double_t chi2_time_t | ( | Double_t | x, |
| Double_t | y, | ||
| Double_t | z, | ||
| Double_t | t, | ||
| Int_t | iRing, | ||
| Int_t | iColumn | ||
| ) |
Definition at line 642 of file TimeRecoTool.cc.
{
Double_t t_exp0 = t +(1.465/300)*sqrt(pow((x-pmt_x_t[iRing][iColumn]),2)+pow((y-pmt_y_t[iRing][iColumn]),2)+pow((z-pmt_z_t[iRing][iColumn]),2));
Double_t t_obs = pmt_time_t[iRing][iColumn];
Double_t chi2 = pow((t_exp0 - t_obs)/pmt_error_t[iRing][iColumn],2);
//std::cout << "iring: " << iRing << " iColumn: " << iColumn << " Theory: " << t_exp0 << " Observed: " << t_obs << " Error: " << pmt_error_t[iRing][iColumn] << " Charge: " << pmt_charge_t[iRing][iColumn] << " chi2: " << chi2 << std::endl;
return chi2;
}
| void IterativeFit | ( | Int_t | n, |
| Int_t | maxstep, | ||
| Double_t | x, | ||
| Double_t | y, | ||
| Double_t | z, | ||
| Int_t | iSelect | ||
| ) |
Definition at line 796 of file TimeRecoTool.cc.
{
Double_t x0 = x;
Double_t y0 = y;
Double_t z0 = z;
Double_t t0 = 0.0;
Double_t x1 = 1e10;
Double_t y1 = 1e10;
Double_t z1 = 1e10;
Double_t t1 = 1e10;
Double_t dis[n];
Double_t tempFunc;
Double_t Func_t, Func_x, Func_y, Func_z;
Double_t tempB, tempC;
Double_t A, Bx, By, Bz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz;
Double_t I00, I01, I02, I03, I11, I22, I33, I12, I13, I23;
Double_t DetH;
//Initializing
Double_t pmtinfo[n][5]; //{pmt x, pmt y, pmt z, pmt time, 1/pmt error^2}
Int_t i = 0;
Double_t T1=0;
Double_t T2=0;
Double_t T3=0;
for (Int_t iring = 0; iring < 8; iring++){
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
if (pmt_flag_t[iring][icolumn] != -1 && pmt_flag_t[iring][icolumn] <= iSelect){
pmtinfo[i][0] = pmt_x_t[iring][icolumn];
pmtinfo[i][1] = pmt_y_t[iring][icolumn];
pmtinfo[i][2] = pmt_z_t[iring][icolumn];
pmtinfo[i][3] = pmt_time_t[iring][icolumn] - median_deltime_t;
pmtinfo[i][4] = 1/(pmt_error_t[iring][icolumn]*pmt_error_t[iring][icolumn]);
dis[i] = sqrt((x0-pmtinfo[i][0])*(x0-pmtinfo[i][0])+(y0-pmtinfo[i][1])*(y0-pmtinfo[i][1])+(z0-pmtinfo[i][2])*(z0-pmtinfo[i][2]));
T1 += pmtinfo[i][4];
T2 += pmtinfo[i][3]*pmtinfo[i][4];
T3 += dis[i]*pmtinfo[i][4];
i++;
}
if(i==n) break; //Get out of loop when done filling info
}
if(i==n) break; //Get out of loop when done filling info
}
Int_t step = 0;
//If T1 is 0, meaningless to carry on
if (T1==0){
t1 = 1e10;
x1 = 1e10;
y1 = 1e10;
z1 = 1e10;
}
else{
t0 = (T2 - (1.465/300)*T3) / T1;
A=2*T1;
while (step < maxstep){
Func_t=0;
Func_x=0;
Func_y=0;
Func_z=0;
Bx=0;
By=0;
Bz=0;
Cxx=0;
Cyy=0;
Czz=0;
Cxy=0;
Cxz=0;
Cyz=0;
//Since this was already done for step 0
if (step>0){
for (Int_t m = 0; m<n; m++){
dis[m] = sqrt((x0-pmtinfo[m][0])*(x0-pmtinfo[m][0])+(y0-pmtinfo[m][1])*(y0-pmtinfo[m][1])+(z0-pmtinfo[m][2])*(z0-pmtinfo[m][2]));
}
}
for (Int_t m = 0; m<n; m++){
Func_t += 2 * (t0 + (1.465/300)*dis[m] - pmtinfo[m][3]) * pmtinfo[m][4];
tempFunc = 2 * (1.465/300) * (t0 + (1.465/300)*dis[m] - pmtinfo[m][3])/dis[m] * pmtinfo[m][4];
Func_x += tempFunc * (x0 - pmtinfo[m][0]);
Func_y += tempFunc * (y0 - pmtinfo[m][1]);
Func_z += tempFunc * (z0 - pmtinfo[m][2]);
tempB = 2 * (1.465/300) * 1/dis[m] * pmtinfo[m][4];
Bx += tempB * (x0 - pmtinfo[m][0]);
By += tempB * (y0 - pmtinfo[m][1]);
Bz += tempB * (z0 - pmtinfo[m][2]);
tempC = -2 * (1.465/300) * (t0 - pmtinfo[m][3])/(dis[m]*dis[m]*dis[m]) * pmtinfo[m][4];
Cxx += tempFunc + (x0 - pmtinfo[m][0])*(x0 - pmtinfo[m][0])*tempC;
Cyy += tempFunc + (y0 - pmtinfo[m][1])*(y0 - pmtinfo[m][1])*tempC;
Czz += tempFunc + (z0 - pmtinfo[m][2])*(z0 - pmtinfo[m][2])*tempC;
Cxy += (x0 - pmtinfo[m][0])*(y0 - pmtinfo[m][1])*tempC;
Cxz += (x0 - pmtinfo[m][0])*(z0 - pmtinfo[m][2])*tempC;
Cyz += (y0 - pmtinfo[m][1])*(z0 - pmtinfo[m][2])*tempC;
}
I00 = Cxx*Cyy*Czz - Cxx*Cyz*Cyz - Cyy*Cxz*Cxz - Czz*Cxy*Cxy + 2*Cxy*Cyz*Cxz;
I01 = Bx*(Cyz*Cyz - Cyy*Czz) + By*(Cxy*Czz - Cxz*Cyz) + Bz*(Cyy*Cxz - Cxy*Cyz);
I02 = By*(Cxz*Cxz - Cxx*Czz) + Bz*(Cyz*Cxx - Cxy*Cxz) + Bx*(Czz*Cxy - Cyz*Cxz);
I03 = Bz*(Cxy*Cxy - Cxx*Cyy) + Bx*(Cxz*Cyy - Cyz*Cxy) + By*(Cxx*Cyz - Cxz*Cxy);
I11 = A*(Cyy*Czz - Cyz*Cyz) - By*By*Czz - Bz*Bz*Cyy + 2*By*Bz*Cyz;
I22 = A*(Czz*Cxx - Cxz*Cxz) - Bz*Bz*Cxx - Bx*Bx*Czz + 2*Bz*Bx*Cxz;
I33 = A*(Cxx*Cyy - Cxy*Cxy) - Bx*Bx*Cyy - By*By*Cxx + 2*Bx*By*Cxy;
I13 = A*(Cxy*Cyz - Cyy*Cxz) + By*By*Cxz - Bx*By*Cyz + Bx*Bz*Cyy - By*Bz*Cxy;
I12 = A*(Cyz*Cxz - Czz*Cxy) + Bz*Bz*Cxy - By*Bz*Cxz + By*Bx*Czz - Bz*Bx*Cyz;
I23 = A*(Cxz*Cxy - Cxx*Cyz) + Bx*Bx*Cyz - Bz*Bx*Cxy + Bz*By*Cxx - Bx*By*Cxz;
DetH = A*I00 + Bx*I01 + By*I02 + Bz*I03;
t1 = t0 - (I00*Func_t + I01*Func_x + I02*Func_y + I03*Func_z)/DetH;
x1 = x0 - (I01*Func_t + I11*Func_x + I12*Func_y + I13*Func_z)/DetH;
y1 = y0 - (I02*Func_t + I12*Func_x + I22*Func_y + I23*Func_z)/DetH;
z1 = z0 - (I03*Func_t + I13*Func_x + I23*Func_y + I33*Func_z)/DetH;
//Convergence reached!
if (fabs(x1-x0)<0.1 && fabs(y1-y0)<0.1 && fabs(z1-z0)<0.1 && fabs(t1-t0)<0.01) break;
//Setting conditions for next iteration
t0 = t1;
x0 = x1;
y0 = y1;
z0 = z1;
step ++;
}
//If for loop did not break, assume it did not converge
if (step == maxstep){
t1 = 1e10;
x1 = 1e10;
y1 = 1e10;
z1 = 1e10;
}
}
final_pos_t[0] = x1;
final_pos_t[1] = y1;
final_pos_t[2] = z1;
final_pos_t[3] = t1 + median_deltime_t;
totalsteps = step;
}//End of IterativeFit
| void SimpleStart | ( | ) |
Definition at line 675 of file TimeRecoTool.cc.
{
//If AdSimple did not provide a good vertex, have to produce an initial guess
Double_t Qrow[8] = {0}; //Total charge in a row
Double_t Qcol[24] = {0}; //Total charge in a column
Double_t Qcol3[24] = {0};//Total charge in 3 adjacent columns
for (Int_t iring = 0; iring < 8; iring++){
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
Qrow[iring] += pmt_charge_t[iring][icolumn];
Qcol[icolumn] += pmt_charge_t[iring][icolumn];
Qcol3[icolumn] += pmt_charge_t[iring][icolumn];
Qcol3[(icolumn+1)%24] += pmt_charge_t[iring][icolumn];
Qcol3[(icolumn+23)%24] += pmt_charge_t[iring][icolumn];
}
}
Int_t maxring = 0, maxcolumn = 0;
Double_t maxQrow = 0.0, maxQcol3 = 0.0;
for (Int_t iring = 0; iring < 8; iring++){
if (Qrow[iring] > maxQrow){
maxQrow = Qrow[iring];
maxring = iring;
}
}
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
if (Qcol3[icolumn] > maxQcol3){
maxQcol3 = Qcol3[icolumn];
maxcolumn = icolumn;
}
}
Double_t phiMinus1 = atan2(pmt_y_t[maxring][(maxcolumn+23)%24],pmt_x_t[maxring][(maxcolumn+23)%24]);
Double_t phi0 = atan2(pmt_y_t[maxring][maxcolumn],pmt_x_t[maxring][maxcolumn]);
Double_t phiPlus1 = atan2(pmt_y_t[maxring][(maxcolumn+1)%24],pmt_x_t[maxring][(maxcolumn+1)%24]);
Double_t phi = ( Qcol[(maxcolumn+23)%24]*phiMinus1 + Qcol[maxcolumn]*phi0 + Qcol[(maxcolumn+1)%24]*phiPlus1 );
phi = phi / ( Qcol[(maxcolumn+23)%24] + Qcol[maxcolumn] + Qcol[(maxcolumn+1)%24]);
Double_t r = SimpleR(maxring,maxcolumn);
initial_pos_t[0] = r*cos(phi);
initial_pos_t[1] = r*sin(phi);
initial_pos_t[2] = pmt_z_t[maxring][maxcolumn];
initial_pos_t[3] = 1e10;
}
| Double_t SimpleR | ( | Int_t | ring, |
| Int_t | column | ||
| ) |
Definition at line 725 of file TimeRecoTool.cc.
{
Double_t r = 1100.0;
std::vector<double> cluster_t1;
Double_t median_t1, early_t1;
//Search for usable PMT corresponding to highest charge row & column
for (Int_t i = 0; i < 3; i++){
Int_t di = int(pow(-1,i) * (i+1)/2);
if (pmt_flag_t[ring][(column+24+di)%24] == 1){
cluster_t1.push_back(pmt_time_t[ring][(column+24+di)%24]);
}
}
if (cluster_t1.size()!=0){
std::sort(cluster_t1.begin(),cluster_t1.end());
median_t1 = cluster_t1[int(cluster_t1.size()/2)];
//Find early_t1 such that it is consistent with the median_t1 to within 2ns
early_t1 = median_t1;
for (Int_t i=0;i<cluster_t1.size();i++){
if (median_t1 - cluster_t1[i] < 2){
early_t1 = cluster_t1[i];
break;
}
}
//Scan 3x3 cluster of PMT on the opposite side
std::vector<double> cluster_t2;
Double_t median_t2, early_t2;
for (Int_t dring = -1; dring <= 1; dring++){
if (ring+dring < 0 || ring+dring > 7) continue;
for (Int_t dcolumn = -1; dcolumn <= 1; dcolumn++){
if (pmt_flag_t[ring+dring][(column+12+dcolumn)%24] == 1){
if (pmt_charge_t[ring+dring][(column+12+dcolumn)%24]>0.5 && fabs(pmt_time_t[ring+dring][(column+12+dcolumn)%24]-early_t1) < 25){
cluster_t2.push_back(pmt_time_t[ring+dring][(column+12+dcolumn)%24]);
}
}
}
}
if (cluster_t2.size()!=0){
//Determine the median time of ths cluster
std::sort(cluster_t2.begin(),cluster_t2.end());
if (cluster_t2.size()%2 == 0){ //even
median_t2 = (cluster_t2[int(cluster_t2.size()/2)-1] + cluster_t2[int(cluster_t2.size()/2)]) / 2;
}
else { //odd
median_t2 = cluster_t2[int(cluster_t2.size()/2)];
}
//Find early_t2 such that it is consistent with the median_t2
early_t2 = median_t2;
for (Int_t i=0;i<cluster_t2.size();i++){
if (median_t2 - cluster_t2[i] < 3){
early_t2 = cluster_t2[i];
break;
}
}
r = 0.5*(early_t2 - early_t1) * (300/1.465);
}
}
return r;
}
| void CalculateDeltime | ( | Double_t | x, |
| Double_t | y, | ||
| Double_t | z | ||
| ) |
Definition at line 556 of file TimeRecoTool.cc.
{
//Calculate deltime = PMT Observed time - Predicted Expected time and store into a vector for chronological sort
std::vector<double> sorted_deltime;
for (Int_t iring = 0; iring < 8; iring++){
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
if (pmt_flag_t[iring][icolumn] != -1){
deltime_t[iring][icolumn] = pmt_time_t[iring][icolumn] - (1.465/300)*sqrt(pow((x-pmt_x_t[iring][icolumn]),2)+pow((y-pmt_y_t[iring][icolumn]),2)+pow((z-pmt_z_t[iring][icolumn]),2));
sorted_deltime.push_back(deltime_t[iring][icolumn]);
}
}
}
std::sort(sorted_deltime.begin(),sorted_deltime.end());
Double_t median_t = sorted_deltime[int(sorted_deltime.size()/2)];
Double_t early_t = 1e10;
//Find early_deltime by finding a 3ns window containing (nPMT/8) number of PMT
for (Int_t i = 0; i<int(sorted_deltime.size()/2); i++){
if (sorted_deltime[i+int(sorted_deltime.size()/8)] - sorted_deltime[i] < 3.0){
early_t = sorted_deltime[i];
break;
}
}
//If first attempt to find early_deltime_t failed ...
if (early_t > 9e9){
//Find early_deltime by finding a 4ns window containing 5 PMT
for (Int_t i = 0; i<int(sorted_deltime.size()/2); i++){
if (sorted_deltime[i+4] - sorted_deltime[i] < 4.0){
early_t = sorted_deltime[i];
break;
}
}
}
//If attempts yield an absurd early_t
if (median_t - early_t > 10 || median_t - early_t < 0){
early_t = median_t - 10.0;
}
early_deltime_t = early_t;
median_deltime_t = median_t;
}
| void SelectPMT | ( | ) |
Definition at line 605 of file TimeRecoTool.cc.
{
for (Int_t i=0; i<maxselect_t; i++){
nPMT_t[i] = 0;
}
Double_t min_window_t = median_deltime_t - early_deltime_t;
Double_t max_window_t = 2.25 - 1.45e-3*sqrt(initial_pos_t[0]*initial_pos_t[0] + initial_pos_t[1]*initial_pos_t[1])*sqrt(initial_energy_t/2.2);
//Ensure that at least half the PMTs are utilized in fit
if (max_window_t <0.0) max_window_t = 0.0;
for (Int_t iring = 0; iring < 8; iring++){
for (Int_t icolumn = 0; icolumn < 24; icolumn++){
if (pmt_flag_t[iring][icolumn] != -1){
pmt_flag_t[iring][icolumn] = -1;
//Too Early
if (deltime_t[iring][icolumn] < median_deltime_t - min_window_t) continue;
for (Int_t i=0; i<maxselect_t; i++){
if (deltime_t[iring][icolumn] < median_deltime_t + max_window_t + i*0.5){
pmt_flag_t[iring][icolumn] = i;
nPMT_t[i]++;
break;
}
}
} //pmt flag
} //icolumn
} //iring
}
Variable Documentation
| const Int_t maxstep_t = 20 |
Definition at line 37 of file TimeRecoTool.cc.
| const Int_t maxselect_t = 5 |
Definition at line 40 of file TimeRecoTool.cc.
| Int_t nPMT_t[maxselect_t] |
Definition at line 41 of file TimeRecoTool.cc.
| Double_t z_phys_min_t = -2400 |
Definition at line 45 of file TimeRecoTool.cc.
| Double_t z_phys_max_t = 2400 |
Definition at line 46 of file TimeRecoTool.cc.
| Double_t x_phys_max_t = 2300 |
Definition at line 47 of file TimeRecoTool.cc.
| Double_t x_phys_min_t = -2300 |
Definition at line 48 of file TimeRecoTool.cc.
| Double_t y_phys_max_t = 2300 |
Definition at line 49 of file TimeRecoTool.cc.
| Double_t y_phys_min_t = -2300 |
Definition at line 50 of file TimeRecoTool.cc.
| Double_t r_phys_max_t = 2300 |
Definition at line 51 of file TimeRecoTool.cc.
| Double_t total_pe_t |
Definition at line 53 of file TimeRecoTool.cc.
| Double_t early_deltime_t |
Definition at line 54 of file TimeRecoTool.cc.
| Double_t median_deltime_t |
Definition at line 55 of file TimeRecoTool.cc.
| Double_t deltime_t[8][24] |
Definition at line 56 of file TimeRecoTool.cc.
| Double_t pmt_charge_t[8][24] |
Definition at line 58 of file TimeRecoTool.cc.
| Double_t pmt_time_t[8][24] |
Definition at line 59 of file TimeRecoTool.cc.
| Double_t pmt_x_t[8][24] |
Definition at line 60 of file TimeRecoTool.cc.
| Double_t pmt_y_t[8][24] |
Definition at line 61 of file TimeRecoTool.cc.
| Double_t pmt_z_t[8][24] |
Definition at line 62 of file TimeRecoTool.cc.
| Double_t pmt_error_t[8][24] |
Definition at line 63 of file TimeRecoTool.cc.
| Int_t pmt_flag_t[8][24] |
Definition at line 64 of file TimeRecoTool.cc.
| Double_t initial_pos_t[4] |
Definition at line 66 of file TimeRecoTool.cc.
| Double_t initial_energy_t |
Definition at line 67 of file TimeRecoTool.cc.
| Int_t totalsteps |
Definition at line 69 of file TimeRecoTool.cc.
| Double_t final_pos_t[4] |
Definition at line 70 of file TimeRecoTool.cc.
| Int_t entry |
Definition at line 80 of file TimeRecoTool.cc.
