/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Defines | Functions | Variables
TimeRecoTool.cc File Reference
#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.

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.

| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:12:10 for TimeReco by doxygen 1.7.4