ArrivalTime.cxx

Go to the documentation of this file.

// $Id: ArrivalTime.cxx,v 1.5 2006/06/15 18:40:51 rhatcher Exp $
//
// ArrivalTime.cxx
//
// Author:  R. Lee
//

#include <cmath>

using namespace std;

#include "RecoBase/ArrivalTime.h"
#include "Conventions/Munits.h"
#include "math.h"

#include "MessageService/MsgService.h"

ClassImp(ArrivalTime)

CVSID("$Id: ArrivalTime.cxx,v 1.5 2006/06/15 18:40:51 rhatcher Exp $");

ArrivalTime::ArrivalTime()
{
  minprob = .01;
  timewindow[0] = -100.;
  timewindow[1] =  200.;

  pe[ 0] = 1;
  pe[ 1] = 2;
  pe[ 2] = 3;
  pe[ 3] = 4;
  pe[ 4] = 5;
  pe[ 5] = 6;
  pe[ 6] = 7;
  pe[ 7] = 8;
  pe[ 8] = 9;
  pe[ 9] = 10;
  pe[10] = 15;
  pe[11] = 20;
  pe[12] = 25;
  pe[13] = 50;
  pe[14] = 100;
  pe[15] = 250;

  tab[ 0] = 2.0;
  tab[ 1] = 1.8;
  tab[ 2] = 1.9;
  tab[ 3] = 2.0;
  tab[ 4] = 2.1;
  tab[ 5] = 2.2;
  tab[ 6] = 2.3;
  tab[ 7] = 2.4;
  tab[ 8] = 2.4;
  tab[ 9] = 2.4;
  tab[10] = 2.0;
  tab[11] = 1.8;
  tab[12] = 1.8;
  tab[13] = 1.7;
  tab[14] = 1.05;
  tab[15] = 0.55;

/*

  a0[ 0] = 38.18;
  a0[ 1] = -0.52543;
  a0[ 2] = -10.35;
  a0[ 3] = -11.864;
  a0[ 4] = -6.9668;
  a0[ 5] = -19.561;
  a0[ 6] = -8.7371;
  a0[ 7] = -18.341;
  a0[ 8] = -10.916;
  a0[ 9] = -16.503;
  a0[10] = -5.0758;
  a0[11] = -3.5335;
  a0[12] = -8.4337;
  a0[13] = -10.632;
  a0[14] = -15.17;
  a0[15] = -12.953;

  a1[ 0] = 486.1;
  a1[ 1] = 646.86;
  a1[ 2] = 723.69;
  a1[ 3] = 796.31;
  a1[ 4] = 773.12;
  a1[ 5] = 838.15;
  a1[ 6] = 861.29;
  a1[ 7] = 886.74;
  a1[ 8] = 860.61;
  a1[ 9] = 822.71;
  a1[10] = 950.38;
  a1[11] = 1018.3;
  a1[12] = 1061.3;
  a1[13] = 1660.9;
  a1[14] = 2444.1;
  a1[15] = 3592.2;

  a2[ 0] = 0.;
  a2[ 1] = 0.;
  a2[ 2] = -94.136;
  a2[ 3] = -176.51;
  a2[ 4] = -172.8;
  a2[ 5] = -236.41;
  a2[ 6] = -262.85;
  a2[ 7] = -289.62;
  a2[ 8] = -295.83;
  a2[ 9] = -293.62;
  a2[10] = -397.63;
  a2[11] = -496.42;
  a2[12] = -766.46;
  a2[13] = -1889.4;
  a2[14] = -3903.9;
  a2[15] = -8204.4;

  a3[ 0] = 0.;
  a3[ 1] = 0.;
  a3[ 2] = 0.;
  a3[ 3] = 0.;
  a3[ 4] = 0.;
  a3[ 5] = 0.;
  a3[ 6] = 0.;
  a3[ 7] = 0.;
  a3[ 8] = 0.;
  a3[ 9] = 0.;
  a3[10] = 0.;
  a3[11] = 0.;
  a3[12] = 125.01;
  a3[13] = 547.77;
  a3[14] = 1579.8;
  a3[15] = 4121.3;

  b0[ 0] = 7.288;
  b0[ 1] = 7.5345;
  b0[ 2] = 7.6679;
  b0[ 3] = 7.87;
  b0[ 4] = 8.0165;
  b0[ 5] = 8.2466;
  b0[ 6] = 8.284;
  b0[ 7] = 8.5896;
  b0[ 8] = 8.5376;
  b0[ 9] = 8.7022;
  b0[10] = 9.6703;
  b0[11] = 10.092;
  b0[12] = 10.573;
  b0[13] = 15.05;
  b0[14] = 11.199;
  b0[15] = 10.245;

  b1[ 0] = -0.12743;
  b1[ 1] = -0.24693;
  b1[ 2] = -0.36993;
  b1[ 3] = -0.50634;
  b1[ 4] = -0.64232;
  b1[ 5] = -0.7944;
  b1[ 6] = -0.88006;
  b1[ 7] = -1.0483;
  b1[ 8] = -1.1270;
  b1[ 9] = -1.2819;
  b1[10] = -1.946;
  b1[11] = -2.4642;
  b1[12] = -3.0376;
  b1[13] = -6.6758;
  b1[14] = -6.4324;
  b1[15] = -8.943;

*/


// the following have already been renormalized

  a0[0] = 0.00380099;
  a0[1] = -8.80699e-05;
  a0[2] = -0.00260464;
  a0[3] = -0.00394667;
  a0[4] = -0.00289652;
  a0[5] = -0.0097678;
  a0[6] = -0.00484191;
  a0[7] = -0.0114235;
  a0[8] = -0.0077634;
  a0[9] = -0.0136922;
  a0[10] = -0.00504345;
  a0[11] = -0.00440376;
  a0[12] = -0.0139975;
  a0[13] = -0.0240564;
  a0[14] = -0.0471868;
  a0[15] = -0.0646101;
  
  a1[0] = 0.0483934;
  a1[1] = 0.108423;
  a1[2] = 0.182121;
  a1[3] = 0.2649;
  a1[4] = 0.321433;
  a1[5] = 0.418531;
  a1[6] = 0.477308;
  a1[7] = 0.552298;
  a1[8] = 0.612061;
  a1[9] = 0.682585;
  a1[10] = 0.944322;
  a1[11] = 1.26909;
  a1[12] = 1.76145;
  a1[13] = 3.75802;
  a1[14] = 7.60245;
  a1[15] = 17.9181;
  
  a2[0] = 0;
  a2[1] = 0;
  a2[2] = -0.0236899;
  a2[3] = -0.0587177;
  a2[4] = -0.0718435;
  a2[5] = -0.118051;
  a2[6] = -0.145666;
  a2[7] = -0.180387;
  a2[8] = -0.210393;
  a2[9] = -0.24361;
  a2[10] = -0.395095;
  a2[11] = -0.618682;
  a2[12] = -1.2721;
  a2[13] = -4.27503;
  a2[14] = -12.1432;
  a2[15] = -40.9239;
  
  a3[0] = 0;
  a3[1] = 0;
  a3[2] = 0;
  a3[3] = 0;
  a3[4] = 0;
  a3[5] = 0;
  a3[6] = 0;
  a3[7] = 0;
  a3[8] = 0;
  a3[9] = 0;
  a3[10] = 0;
  a3[11] = 0;
  a3[12] = 0.207481;
  a3[13] = 1.23941;
  a3[14] = 4.91402;
  a3[15] = 20.5572;
  
  b0[0] = -1.92681;
  b0[1] = -1.15934;
  b0[2] = -0.619546;
  b0[3] = -0.138392;
  b0[4] = 0.2311;
  b0[5] = 0.644398;
  b0[6] = 0.785976;
  b0[7] = 1.20838;
  b0[8] = 1.28904;
  b0[9] = 1.60773;
  b0[10] = 2.75615;
  b0[11] = 3.40441;
  b0[12] = 4.17189;
  b0[13] = 8.95878;
  b0[14] = 5.42604;
  b0[15] = 4.94429;

  b1[ 0] = -0.12743;
  b1[ 1] = -0.24693;
  b1[ 2] = -0.36993;
  b1[ 3] = -0.50634;
  b1[ 4] = -0.64232;
  b1[ 5] = -0.7944;
  b1[ 6] = -0.88006;
  b1[ 7] = -1.0483;
  b1[ 8] = -1.1270;
  b1[ 9] = -1.2819;
  b1[10] = -1.946;
  b1[11] = -2.4642;
  b1[12] = -3.0376;
  b1[13] = -6.6758;
  b1[14] = -6.4324;
  b1[15] = -8.943;


//  Renormalize();

}

Double_t ArrivalTime::Weight(Double_t npe) const
{
  Double_t weight = 1./Uncertainty(npe);
  return weight*weight;
}

Double_t ArrivalTime::Uncertainty(Double_t npe) const
{
  Double_t sigma = 0.; // sigma is uncertainty in arrival of first p.e.

// a simple monte carlo was done to calculate the uncertainty in the
// arrival time of the first p.e. as a function of # of p.e.
//
// the only assumption made was about the shape of the pulse, a 2 ns risetime
// and an 8 ns decay time of the fluor
//
// rlee feb 2002
//

  if (npe<10.) {
    if (npe>0.) sigma = 0.90+8.0/npe;
    else {
      sigma = 1.0e+5;
      MSG("ArrivalTime", Msg::kWarning)
        << "Uncertainty was passed npe = " << npe
        << ", use sigma = " << sigma << endl;
    }
  }
  else {
    sigma = 5.0*exp(-0.5*log(npe));
  }

  return sigma;
}

void ArrivalTime::SetTimeWindow(Double_t t0, Double_t t1)
{
  if (t0<t1) {
    timewindow[0] = t0;
    timewindow[1] = t1;
  }
  else {
    timewindow[0] = t1;
    timewindow[1] = t0;
  }
  Renormalize();
}

void ArrivalTime::SetMinimumProbability(Double_t prob)
{
  minprob = prob;
  Renormalize();
}

void ArrivalTime::Renormalize()
{
  for (int i=0; i<16; i++) {
    Double_t integral = a0[i]*tab[i]+a1[i]*tab[i]*tab[i]/2.+a2[i]*tab[i]*tab[i]*tab[i]/3.+a3[i]*tab[i]*tab[i]*tab[i]*tab[i]/4.-exp(b0[i]+b1[i]*tab[i])/b1[i];
    a0[i] *= (1.-minprob)/integral;
    a1[i] *= (1.-minprob)/integral;
    a2[i] *= (1.-minprob)/integral;
    a3[i] *= (1.-minprob)/integral;
    b0[i] += log((1.-minprob)/integral);
  }

}

Double_t ArrivalTime::PDF(Double_t npe, Double_t t) const
{
// input time is in seconds, convert into ns
  t /= Munits::ns;
  if (npe<=0. || t<timewindow[0] || t>timewindow[1]) {
    return 0.;
  }
  Double_t minpdf = minprob/(timewindow[1]-timewindow[0]);
  if (t<0.) {
    return minpdf;
  }
  Double_t p0[2] = { 0, 0 };
  Double_t pe0[2] = {-1.,1.};
  for (int i=0; i<16; i++) {
    if (npe>=pe[i]) {
      if (t<tab[i]) {
        p0[1] = a0[i]+a1[i]*t+a2[i]*t*t+a3[i]*t*t*t;
      }
      else {
        p0[1] = exp(b0[i]+b1[i]*t);
      }
      if (p0[1]<minpdf) p0[1] = minpdf;
      pe0[1] = (Double_t)(pe[i]);
    }
  }
  for (int i=15; i>=0; i--) {
    if (npe<=pe[i]) {
      if (t<tab[i]) {
        p0[0] = a0[i]+a1[i]*t+a2[i]*t*t+a3[i]*t*t*t;
      }
      else {
        p0[0] = exp(b0[i]+b1[i]*t);
      }
      if (p0[0]<minpdf) p0[0] = minpdf;
      pe0[0] = (Double_t)(pe[i]);
    }
  }
  if (pe0[0]==pe0[1]) {
    return (p0[0]+minpdf);
  }
  if (pe0[0]<0.) {
    return (p0[1]+minpdf);
  }
  if (pe0[1]<0.) {
    return (p0[0]+minpdf);
  }
  return (minpdf+ p0[0] + (p0[1]-p0[0])/(pe0[1]-pe0[0])*(npe-pe0[0]));
}

---