LArPhysWaveHECTool.cxx File Reference

#include "LArCalibUtils/LArPhysWaveHECTool.h"
#include "LArIdentifier/LArOnline_SuperCellID.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <TGraphSmooth.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>

Include dependency graph for LArPhysWaveHECTool.cxx:

Go to the source code of this file.

Classes
struct	gsl_f_params
struct	usertim

Functions
static Double_t	Tp (Double_t *t, Double_t *par)
static Double_t	Tp_gsl (Double_t *t, Double_t *par)
static Double_t	DTp (Double_t *t, Double_t *par)
static Double_t	Tp4 (Double_t t)
static Double_t	Tp5_gsl (Double_t t)
static Double_t	Tp4_gsl (Double_t t)
static Double_t	Tp5 (Double_t t)
static Double_t	Tc (Double_t t)
Double_t	normalize (TF1 *func, Double_t *rampl=NULL, Double_t from=0., Double_t to=0., Double_t step=1.)
Double_t	normalize_prof (TProfile *histo, Double_t *rampl=NULL)
static Double_t	f4_g (Double_t *x, Double_t *par)
static Double_t	f5_g (Double_t *x, Double_t *par)
static double	f4_gsl (double x, void *par)
static double	f5_gsl (double x, void *par)
void	Smooth_new (Double_t *x, Double_t *y, Double_t *yref, Int_t ymax)
void	SmootKern (Double_t *x, Double_t *y, Int_t ymax, Double_t band=10.)
static Double_t	P (Double_t *tt, Double_t *par)
static Double_t	Rd (Double_t t)
static Double_t	Ro (Double_t t)
static double	f4_gsl (Double_t x, void *par)
static double	f5_gsl (Double_t x, void *par)

Variables
	ATLAS_NO_CHECK_FILE_THREAD_SAFETY
const int	NMAX =1000
static Double_t	tau0
static Double_t	taup
static Double_t	tauz
static Double_t	taus
static Double_t	rc
static Double_t	rz
static Double_t	rs
static Double_t	sp
static Double_t	s0
static Double_t	sz
static Double_t	ss
static Double_t	sc
static Double_t	fs
static Double_t	tdr
static Double_t	tc
static Double_t	a
static Double_t	t0
static Double_t	al
static Double_t	x [NMAX]
static Double_t	y [NMAX]
static Double_t	STEP
static Double_t *	pref = NULL
static TSpline3 *	splin = NULL
int	LArWaveFlag =0
int	Tp4_gsl_err_counter =0
int	Tp5_gsl_err_counter =0

Function Documentation

DTp()

Double_t DTp ( Double_t * t, Double_t * par )

static

Definition at line 325 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res,t;
 
 static TF1 *ff=NULL;
 static TGraph *gg=NULL;
 static Double_t xx[21],yy[21];
 static Int_t j;
 
 if(*tt < t0) return 0.;
 t = *tt;
 
 if(ff==NULL) ff=new TF1("ff","pol2",0.,800.);
 if(gg==NULL) gg=new TGraph(11,xx,yy);
 
 
 ff->SetRange(*tt-5.,*tt+5.);
 for(j=-10; j<=10; ++j) {
           res = t + j/2.;
           gg->SetPoint(j+10,t+j/2.,Tp(&res,par));
          }
 if(gg->GetFunction("ff")) gg->GetFunction("ff")->Delete();
 gg->Fit("ff","QRI0");
 res = 2*t*ff->GetParameter(2) + ff->GetParameter(1);
 
 
 return res; 
}

f4_g()

Double_t f4_g ( Double_t * x, Double_t * par )

static

Definition at line 625 of file LArPhysWaveHECTool.cxx.

{
  static Double_t t;
 
  t = par[0];
 
  return  Tc(*x) * Rd(t - *x);
}

f4_gsl() [1/2]

double f4_gsl ( double x, void * par )

static

f4_gsl() [2/2]

double f4_gsl ( Double_t x, void * par )

static

Definition at line 504 of file LArPhysWaveHECTool.cxx.

{
  struct gsl_f_params * params = (struct gsl_f_params *)par;
  static Double_t t;
  
  t = params -> a;
 
 
  return  Tc(x) * Rd(t - x);
}

f5_g()

Double_t f5_g ( Double_t * x, Double_t * par )

static

Definition at line 651 of file LArPhysWaveHECTool.cxx.

{
  static Double_t t;
 
  t = par[0];
 
  return  Tc(*x) * Ro(t - *x);
}

f5_gsl() [1/2]

double f5_gsl ( double x, void * par )

static

f5_gsl() [2/2]

double f5_gsl ( Double_t x, void * par )

static

Definition at line 594 of file LArPhysWaveHECTool.cxx.

{
  struct gsl_f_params * params = (struct gsl_f_params *)par;
  static Double_t t;
  
  t = params -> a;
 
  return  Tc(x) * Ro(t - x);
}

normalize()

Double_t normalize	(	TF1 *	func,
		Double_t *	rampl = NULL,
		Double_t	from = 0.,
		Double_t	to = 0.,
		Double_t	step = 1. )

Definition at line 825 of file LArPhysWaveHECTool.cxx.

{
 static Double_t amax, tmax,t, ampl;
 static Double_t t1,t2,tt1,tt2;
 static Int_t i;
 
 if((from == to) || (from > to)) 
     func->GetRange(t1,t2);
 else {
    func->GetRange(tt1,tt2);
    if(from>tt1) t1 = from; else t1 = tt1;
    if(to<tt2) t2 = to; else t2 = tt2;
  }
 
 amax = 0; tmax = t1;
 for(t=t1; t<=t2; t+=step)
  {
   if(func->Eval(t) > amax) {
     amax = func->Eval(t); tmax = t; 
    }
   }
 Double_t x[20], y[20];
 i=0;
 for(t=tmax-5; t<tmax+5; ++t) {
    x[i] = t;
    y[i] = func->Eval(t);
    ++i;
   }
 TGraph *g = new TGraph(10,x,y);
 g->Fit("pol2","QIE0");
 tmax = - g->GetFunction("pol2")->GetParameter(1) / (2 *g->GetFunction("pol2")->GetParameter(2));
 amax = g->GetFunction("pol2")->Eval(tmax);
 
 ampl = func->GetParameter(1);
 if(amax != 0) func->SetParameter(1,ampl/amax);
   else tmax = 0;
 cout << "Normalize: tmax="<<tmax << " ampl=" << ampl<<" amax="<<amax << endl;
 if(rampl != NULL) *rampl = amax;
 g->Delete();
 return tmax;
}

normalize_prof()

Double_t normalize_prof	(	TProfile *	histo,
		Double_t *	rampl = NULL )

Definition at line 868 of file LArPhysWaveHECTool.cxx.

{
 Double_t amax, tmax, x, y, err;
 Int_t nbins,i, j, nent;
 amax = histo->GetBinContent(0);
 tmax = histo->GetBinCenter(0);
 nbins = histo->GetNbinsX();
 for(i=1; i<nbins; i++)
  {
   if(amax < histo->GetBinContent(i))
     {
      amax = histo->GetBinContent(i);
      tmax =  histo->GetBinCenter(i);
     }
  }
 
 histo->Fit("pol2","0QI","",tmax-10., tmax+10.);
 tmax = - histo->GetFunction("pol2")->GetParameter(1) / (2 * histo->GetFunction("pol2")->GetParameter(2));
 amax = histo->GetFunction("pol2")->Eval(tmax);
 TProfile *htmp = new TProfile(*histo);
 histo->Reset();
 for(i=0; i<=nbins; i++)
  {
   if(htmp->GetBinEntries(i) == 0) continue;
   x = htmp->GetBinCenter(i);
   y = htmp->GetBinContent(i);
   nent=int(htmp->GetBinEntries(i));
   err = htmp->GetBinError(i);
// Hack for stupid root v.3, which do not like the one entries in profile histo
   for(j=0; j<7; ++j) {
      if(nent != 0) {
           histo->Fill(x,y/amax+sqrt(2)*err/nent);
           histo->Fill(x,y/amax-sqrt(2)*err/nent);
      } else {
           histo->Fill(x,y/amax+sqrt(2)*err);
           histo->Fill(x,y/amax-sqrt(2)*err);
      }
    }
  }
 
 if(rampl != NULL) *rampl = amax;
 
 return tmax;
}

P()

Double_t P ( Double_t * tt, Double_t * par )

static

Definition at line 280 of file LArPhysWaveHECTool.cxx.

{
   return par[1]*splin->Eval(*tt+par[0]);
}

Rd()

Double_t Rd ( Double_t t )

static

Definition at line 285 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res;
 
 if(t<0 || t>tdr) return 0;
 
 res = rc*((1+sc*al*tdr)*exp(-t*sc*al)-1) + rz*((1+sz*tdr)*exp(-t*sz)-1);
 res += rs*((1+ss*tdr)*exp(-t*ss)-1) - 1;
 res *= fs/tdr;
 
 return res;
}

Ro()

Double_t Ro ( Double_t t )

static

Definition at line 298 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res;
 
 if(t<tdr) return 0;
 
 res = rc*((1+sc*al*tdr)*exp(-t*sc*al) - exp((tdr-t)*sc*al));
 res += rz*((1+sz*tdr)*exp(-t*sz) - exp((tdr-t)*sz));
 res += rs*((1+ss*tdr)*exp(-t*ss) - exp((tdr-t)*ss));
 res *= fs/tdr;
 
 return res;
}

Smooth_new()

void Smooth_new	(	Double_t *	x,
		Double_t *	y,
		Double_t *	yref,
		Int_t	ymax )

Definition at line 925 of file LArPhysWaveHECTool.cxx.

{
// Smoothing the difference between the yref and y, 
// then y = y + smoothed differnce
// repeating few times
  static Double_t *dif = NULL;
  static Int_t i,j;
  
  if(dif) delete[] dif;
  dif = new Double_t[ymax];
  for(j=0; j<4; ++j) {
    for(i=0; i<ymax; ++i) dif[i] = y[i] - yref[i];
    SmootKern(x,dif,ymax, 10.);
    for(i=0; i<ymax; ++i) yref[i] += dif[i];
  }
}

SmootKern()

void SmootKern	(	Double_t *	x,
		Double_t *	y,
		Int_t	ymax,
		Double_t	band = 10. )

Definition at line 912 of file LArPhysWaveHECTool.cxx.

{
   Double_t *xx, *yy;
   TGraph *g=new TGraph(ymax,x,y);
   TGraphSmooth *gs = new TGraphSmooth();
   TGraph *g1 = gs->SmoothKern(g,"normal",band);
   xx = g1->GetX(); yy = g1->GetY();
   for(Int_t i=0; i<g->GetN(); ++i) {
      x[i] = xx[i]; y[i] = yy[i];
   }
   delete g; delete gs;
}

Tc()

Double_t Tc ( Double_t t )

static

Definition at line 312 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res;
 
 if(t<t0) return 0;
 
 res = splin->Eval(t);
 
 
 return res;
}

Tp()

Double_t Tp ( Double_t * t, Double_t * par )

static

Definition at line 354 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res;
 static Double_t t,t1,t4,t5;
 
 if(*tt < t0 + par[0]) return 0;
// if(par[0] < 0) return 0;
 
 t = *tt - par[0];
 taus = par[2];
 
// Redefinition of poles
 
  sp = 1/taup; s0 = 1/tau0; sz = 1/tauz; ss = 1/taus; sc = 1/tc;
  //
// Coefficients
 
 rc = pow(al,3)*sc*sc + al*al*(-sc*sc-sc*s0-sc*sp);
 rc += al*(sc*s0+sc*sp+sp*s0) - sp*s0;
 rc /= -al*(al*al*sc*sc + al*(-sc*ss-sc*sz) + sz*ss);
 
 rz = sz*sz*(sc+s0+sp) - pow(sz,3) - sz*sc*(s0+sp) + sp*s0*(sc-sz);
 rz /= sz*(al*(-sc*sz+sc*ss) + sz*sz - sz*ss);
 
 rs = -ss*ss*(sc+s0+sp-ss) + sc*ss*(s0+sp) + ss*sp*s0 - sc*sp*s0;
 rs /= ss*(al*(-sc*sz+sc*ss) + sz*ss - ss*ss);
 
 fs = sz*ss / (s0*sp);
 
 t1 = fs*Tc(t);
 
 t4 = Tp4(t); t5 = Tp5(t);
 res = t1 + t4 + t5;
 
 return par[1]*res;
}

Tp4()

Double_t Tp4 ( Double_t t )

static

Definition at line 607 of file LArPhysWaveHECTool.cxx.

{
  static Double_t a,b;//, par[1];
  static const Double_t epsrel = 0.001;
  static TF1 *fun4 = NULL;
 
  if(fun4 == NULL) fun4 = new TF1("fun4",&f4_g,0.,NMAX,1);
 
  a = t - tdr;
  b = t;
 
 
  fun4->SetParameter(0,t);
  return fun4->Integral(a,b,epsrel);
 
}

Tp4_gsl()

Double_t Tp4_gsl ( Double_t t )

static

Definition at line 429 of file LArPhysWaveHECTool.cxx.

{
  static Double_t a,b;
  
 
 
  //BEG: variables for QAWO
  static double L;
  static const double omega=0.0;
  static const size_t n=100;
  static const size_t limit=1001; // "n" must be < limit
  static const double epsabs=0.0;
  static double epsrel=1e-10;
  double abserr;
  static gsl_integration_workspace  *w=0;
  static gsl_integration_qawo_table *wf=0;
  //END: variables for QAWO
 
  struct gsl_f_params params = {t};
  gsl_function gsl_func;
  
  a = t - tdr;
  if(a < t0) a = t0;
  b = t;
  
  //for gsl_integration_qawo:
  L=b-a;
  if(!w) w = gsl_integration_workspace_alloc(limit);
  if(!wf) wf=gsl_integration_qawo_table_alloc(omega,L,GSL_INTEG_COSINE,n); 
 
  gsl_integration_qawo_table_set_length(wf, L);
  //end of gsl_integration_qawo
 
  gsl_func.function = &f4_gsl;
  gsl_func.params   = &params;
  double result;
 
  // set off automatic GSL error handling
  gsl_set_error_handler_off();
  // set integration precision 
  if(b>t0+130) epsrel=1e-8;
  else         epsrel=1e-10;
 
  unsigned short int INTEGRAL_ITER=0;
  const unsigned short int INTEGRAL_ITER_MAX=10;
  int status=-1;
 
    // compute integral with the starting high precision
  while((status!=0 || INTEGRAL_ITER==0) && INTEGRAL_ITER<INTEGRAL_ITER_MAX){
    status = gsl_integration_qawo (&gsl_func,a, epsabs, epsrel,limit,w,wf,&result,&abserr);
 
    // if failes, descrease the precision
    if(status == GSL_EROUND)      epsrel*=10;
    // if different error from GSL_EROUND, annonce integration ERROR
    else if(status!=0){
      if(Tp4_gsl_err_counter==0)
    std::cout<< "WARNING Tp4_gsl:QAWO ERROR (not roundoff): status="<<gsl_strerror(status)<<" ERRcode="<<status << " LArWaveFlag set to -1 !!"<< std::endl ;
      status=0;
      LArWaveFlag=-1;
      Tp4_gsl_err_counter++;
    }
    INTEGRAL_ITER++;
  }
  //if iteration of integration reaches the limit (=> bad precision) set LArWaveFlag to -1 
  if(INTEGRAL_ITER==INTEGRAL_ITER_MAX && status!=0){
    std::cout<<"WARNING  Integration of Tp4_gsl FAILED !!!! #iterations reached the MAXIMUM!!! status="<<gsl_strerror(status)<<" ERRcode="<<status << " LArWaveFlag set to -1 !!"<<std::endl ;
    LArWaveFlag=-1;
  }
 
 
  return result;
 
}

Tp5()

Double_t Tp5 ( Double_t t )

static

Definition at line 635 of file LArPhysWaveHECTool.cxx.

{
  static Double_t b;//, par[1];
  static const Double_t a = 0., epsrel = 0.001;
  static TF1 *fun5 = NULL;
 
  if(fun5 == NULL) fun5 = new TF1("fun5",&f5_g,0.,NMAX,1);
 
  b = t - tdr;
  fun5->SetParameter(0,t);
  return fun5->Integral(a,b,epsrel);
}

Tp5_gsl()

Double_t Tp5_gsl ( Double_t t )

static

Definition at line 517 of file LArPhysWaveHECTool.cxx.

{
   if(t<t0 + tdr) return 0.;
  Double_t b;
  
 
 
  //BEG: variables for QAWO
  static double L;
  static const double omega=0.0;
  static const size_t n=100;
  static const size_t limit=1001; // "n" must be < limit
  static const double epsabs=0.0;
  static double epsrel = 1e-10;
  double abserr;
  static gsl_integration_workspace  *w=0;
  static gsl_integration_qawo_table *wf=0;
  //END: variables for QAWO
 
  struct gsl_f_params params = {t};
  gsl_function gsl_func;
  double a = 0.;
  b = t - tdr;
  
  //for gsl_integration_qawo:
  L=b-a;
  if(!w) w = gsl_integration_workspace_alloc(limit);
  if(!wf) wf=gsl_integration_qawo_table_alloc(omega,L,GSL_INTEG_COSINE,n); 
 
  gsl_integration_qawo_table_set_length(wf, L);
  //end of gsl_integration_qawo
 
  gsl_func.function = &f5_gsl;
  gsl_func.params   = &params;
 
  
 
  // set off automatic GSL error handling
  gsl_set_error_handler_off();
  // set integration precision 
  if(b>t0+60) epsrel=1e-8;    // behinde the peak = tale region
  else         epsrel=1e-10;  // peak/amplitude region
 
  double result;
  unsigned short int INTEGRAL_ITER=0;
  const unsigned short int INTEGRAL_ITER_MAX=10;
  int status=-1;
 
  while((status!=0 || INTEGRAL_ITER==0) && INTEGRAL_ITER<INTEGRAL_ITER_MAX){
    status = gsl_integration_qawo (&gsl_func,a, epsabs, epsrel,limit,w,wf,&result,&abserr);
    // if failes, descrease the precision
    if(status == GSL_EROUND)      epsrel*=10;
    // if different error from GSL_EROUND, annonce integration ERROR
    else if(status!=0){
      if(Tp5_gsl_err_counter==0)
    std::cout <<"WARNING: Tp5_gsl:QAWO ERROR (not roundoff): status="<<gsl_strerror(status)<<" ERRcode="<<status << " LArWaveFlag set to -1 !!"<< std::endl ;
      status=0;
      LArWaveFlag=-1;
      Tp5_gsl_err_counter++;
    }
    INTEGRAL_ITER++;
  }
  //if iteration of integration reaches the limit (=> bad precision) set LArWaveFlag to -1 
  if(INTEGRAL_ITER==INTEGRAL_ITER_MAX && status!=0){
    std::cout<<"WARNING: Integration of Tp5_gsl FAILED: #iterations reached the MAXIMUM!!!! status="<<gsl_strerror(status)<<" ERRcode="<<status << " LArWaveFlag set to -1 !!"<<std::endl ;
    LArWaveFlag=-1;
  }
 
 
 
  return result;
 
}

Tp_gsl()

Double_t Tp_gsl ( Double_t * t, Double_t * par )

static

Definition at line 391 of file LArPhysWaveHECTool.cxx.

{
 static Double_t res;
 static Double_t t,t1,t4,t5;
 
 if(*tt < t0 + par[0]) return 0;
 
 t = *tt - par[0];
 taus = par[2];
 
// Redefinition of poles
 
  sp = 1/taup; s0 = 1/tau0; sz = 1/tauz; ss = 1/taus; sc = 1/tc;
  //
// Coefficients
 
 rc = pow(al,3)*sc*sc + al*al*(-sc*sc-sc*s0-sc*sp);
 rc += al*(sc*s0+sc*sp+sp*s0) - sp*s0;
 rc /= -al*(al*al*sc*sc + al*(-sc*ss-sc*sz) + sz*ss);
 
 rz = sz*sz*(sc+s0+sp) - pow(sz,3) - sz*sc*(s0+sp) + sp*s0*(sc-sz);
 rz /= sz*(al*(-sc*sz+sc*ss) + sz*sz - sz*ss);
 
 rs = -ss*ss*(sc+s0+sp-ss) + sc*ss*(s0+sp) + ss*sp*s0 - sc*sp*s0;
 rs /= ss*(al*(-sc*sz+sc*ss) + sz*ss - ss*ss);
 
 fs = sz*ss / (s0*sp);
 
 t1 = fs*Tc(t);
 
 t4 = Tp4_gsl(t); t5 = Tp5_gsl(t);
 res = t1 + t4 + t5;
 
 return par[1]*res;
}

Variable Documentation

a

Double_t a

static

Definition at line 38 of file LArPhysWaveHECTool.cxx.

al

Double_t al

static

Definition at line 38 of file LArPhysWaveHECTool.cxx.

ATLAS_NO_CHECK_FILE_THREAD_SAFETY

ATLAS_NO_CHECK_FILE_THREAD_SAFETY

Definition at line 17 of file LArPhysWaveHECTool.cxx.

fs

Double_t fs

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

LArWaveFlag

int LArWaveFlag =0

Definition at line 57 of file LArPhysWaveHECTool.cxx.

NMAX

const int NMAX =1000

Definition at line 22 of file LArPhysWaveHECTool.cxx.

pref

Double_t * pref = NULL

static

Definition at line 39 of file LArPhysWaveHECTool.cxx.

rc

Double_t rc

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

rs

Double_t rs

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

rz

Double_t rz

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

s0

Double_t s0

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

sc

Double_t sc

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

sp

Double_t sp

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

splin

TSpline3* splin = NULL

static

Definition at line 51 of file LArPhysWaveHECTool.cxx.

ss

Double_t ss

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

STEP

Double_t STEP

static

Definition at line 39 of file LArPhysWaveHECTool.cxx.

sz

Double_t sz

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

t0

Double_t t0

static

Definition at line 38 of file LArPhysWaveHECTool.cxx.

tau0

Double_t tau0

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

taup

Double_t taup

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

taus

Double_t taus

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

tauz

Double_t tauz

static

Definition at line 37 of file LArPhysWaveHECTool.cxx.

tc

Double_t tc

static

Definition at line 38 of file LArPhysWaveHECTool.cxx.

tdr

Double_t tdr

static

Definition at line 38 of file LArPhysWaveHECTool.cxx.

Tp4_gsl_err_counter

int Tp4_gsl_err_counter =0

Definition at line 57 of file LArPhysWaveHECTool.cxx.

Tp5_gsl_err_counter

int Tp5_gsl_err_counter =0

Definition at line 57 of file LArPhysWaveHECTool.cxx.

x

Double_t x[NMAX]

static

Definition at line 39 of file LArPhysWaveHECTool.cxx.

y

Double_t y[NMAX]

static

Definition at line 39 of file LArPhysWaveHECTool.cxx.