8 #include <gsl/gsl_integration.h>
9 #include <gsl/gsl_errno.h>
11 #include <TGraphSmooth.h>
24 static Double_t Tp(Double_t *
t, Double_t *
par);
25 static Double_t Tp_gsl(Double_t *
t, Double_t *
par);
26 static Double_t DTp(Double_t *
t, Double_t *
par);
27 static Double_t Tp4(Double_t
t);
28 static Double_t Tp5_gsl(Double_t
t);
29 static Double_t Tp4_gsl(Double_t
t);
30 static Double_t Tp5(Double_t
t);
31 static Double_t Tc(Double_t
t);
32 Double_t
normalize(TF1 *func, Double_t *rampl=NULL,
33 Double_t from=0., Double_t
to = 0., Double_t
step = 1.);
37 static Double_t tau0, taup, tauz, taus, rc, rz, rs, sp, s0, sz, ss, sc, fs,
41 static Double_t f4_g(Double_t *
x, Double_t *
par);
42 static Double_t f5_g(Double_t *
x, Double_t *
par);
44 static double f4_gsl(
double x,
void *
par);
45 static double f5_gsl(
double x,
void *
par);
51 static TSpline3 *splin = NULL;
54 void SmootKern(Double_t *
x, Double_t *
y, Int_t
ymax, Double_t band = 10.);
64 declareInterface<LArPhysWaveHECTool>(
this);
93 return StatusCode::SUCCESS;
100 int & LArPhysWaveFlag) {
119 <<
LArWaveFlag <<
" LArPhysWaveFlag="<< LArPhysWaveFlag );
123 return StatusCode::SUCCESS;
147 ATH_MSG_VERBOSE (
"*** Normalisation \t|-> YES (CaliWave peak = " << peak <<
")" );
151 gCaliMB = gCaliMB * (1./peak) ;
171 unsigned int CaliWaveBins = gCaliMB.
getSize();
172 double const BinWidth=25.0/24;
173 double const HalfBinWidth=BinWidth/2.0;
174 double par[3]={0.0,1.0,0.10};
175 bool uset0=
false,
norm=
false;
177 double const QUAL_REQ=peak_tmp*0.04;
178 double const QUAL_REQ_AMPL=peak_tmp*0.3;
180 double DIFF=0,DIFF_AMPL=0;
181 unsigned short int repro_count=0;
182 int idx_bad_time=-1,idx_bad_time_ampl=-1;
184 const unsigned short int ITER_MAX = 10;
185 double Ampl_problem=-1, Time_problem=-1, Ampl_problem_ampl=-1, Time_problem_ampl=-1;
186 double CALIWAVE_SHIFT=0;
187 TProfile pcal(
"pcal",
"Amplitude vs Time",CaliWaveBins,0-HalfBinWidth,CaliWaveBins*BinWidth-HalfBinWidth);
188 while(( (
DIFF>=QUAL_REQ || DIFF_AMPL>=QUAL_REQ_AMPL) && repro_count<ITER_MAX)
193 CALIWAVE_SHIFT=0.0005*repro_count;
195 CALIWAVE_SHIFT=1.0*repro_count;
200 ATH_MSG_INFO (repro_count<<
". Iteration of INTEGRATION: CALIWAVE IS MOVED UP by "<<CALIWAVE_SHIFT<<
" ADC units");
201 if(DIFF_AMPL>=QUAL_REQ_AMPL)
202 ATH_MSG_INFO (
"Problematic DIFF_AMPL bin="<<idx_bad_time_ampl<<
" AmplPhysGSL="<<Ampl_problem_ampl<<
" Time="<< Time_problem_ampl <<
" Deviation="<<DIFF_AMPL<<
" ADC units"<<
" Peak="<<peak_tmp);
204 ATH_MSG_INFO (
"Problematic DIFF bin="<<idx_bad_time<<
" AmplPhysGSL="<<Ampl_problem<<
" Time="<<Time_problem<<
" Deviation="<<
DIFF<<
" ADC units"<<
" Peak="<<peak_tmp);
208 for(
unsigned tbin=1;tbin<=gCaliMB.
getSize();++tbin){
209 double AmplCalib = gCaliMB.
getSample(tbin-1)+CALIWAVE_SHIFT;
211 pcal.Fill(
Time,AmplCalib,1);
215 double parCL[2]={wfParam.
tcal(),wfParam.
fstep()};
218 if (deriv)
delete deriv;
223 PhysWaveFunc= (TF1*)pwf->Clone();
224 DIFF=0;DIFF_AMPL=0;Ampl_problem=-1; Time_problem=-1; Ampl_problem_ampl=-1; Time_problem_ampl=-1;
227 for(
int i=1;
i<pcal.GetNbinsX(); ++
i){
228 double TimePhys1 = pcal.GetBinCenter(
i+1);
229 double TimePhys0 = TimePhys1-BinWidth;
230 double DIFF_TMP = PhysWaveFunc->Eval(TimePhys1)-PhysWaveFunc->Eval(TimePhys0);
231 if(TimePhys1>200 &&
DIFF<fabs(DIFF_TMP)){
234 Ampl_problem=PhysWaveFunc->Eval(TimePhys1);
235 Time_problem=TimePhys1;
237 else if(TimePhys1<=200 && DIFF_AMPL<fabs(DIFF_TMP)){
238 DIFF_AMPL=fabs(DIFF_TMP);
240 Ampl_problem_ampl=PhysWaveFunc->Eval(TimePhys1);
241 Time_problem_ampl=TimePhys1;
246 if(repro_count>=ITER_MAX && (
DIFF>=QUAL_REQ || DIFF_AMPL>=QUAL_REQ_AMPL)){
247 ATH_MSG_WARNING (
"FT="<<
FT<<
" Slot="<<Slot<<
" Ch="<<
Channel<<
" Gain="<<
gain<<
" #iterations for CALIWAVE increasing reached the limit! LArWaveFlag set to -1" );
254 for(
unsigned int i=0;
i<gCaliMB.
getSize(); ++
i){
260 predLArPhysWave.
setSample(
i,PhysWaveFunc->Eval(
Time)-CALIWAVE_SHIFT);
269 MphysMcali = predLArPhysWave.
getSample( wHelper.
getMax(predLArPhysWave) ) ;
270 std::cout<<
"Normalized: MphysMcali="<<MphysMcali<<std::endl;
272 MphysMcali = predLArPhysWave.
getSample( wHelper.
getMax(predLArPhysWave) ) /
275 ATH_MSG_VERBOSE (
"*** Physics waveform\t|-> m_MphysMcali = " << MphysMcali <<
" adc=" <<
adc);
280 static Double_t P(Double_t *
tt, Double_t *
par)
282 return par[1]*splin->Eval(*
tt+
par[0]);
285 static Double_t Rd(Double_t
t)
289 if(t<0 || t>tdr)
return 0;
291 res = rc*((1+
sc*al*tdr)*
exp(-
t*sc*al)-1) + rz*((1+sz*tdr)*
exp(-
t*sz)-1);
298 static Double_t Ro(Double_t
t)
312 Double_t Tc(Double_t
t)
318 res = splin->Eval(
t);
325 Double_t DTp(Double_t *
tt, Double_t *
par)
327 static Double_t
res,
t;
330 static TGraph *gg=NULL;
331 static Double_t xx[21],
yy[21];
334 if(*
tt < t0)
return 0.;
337 if(
ff==NULL)
ff=
new TF1(
"ff",
"pol2",0.,800.);
338 if(gg==NULL) gg=
new TGraph(11,xx,
yy);
341 ff->SetRange(*
tt-5.,*
tt+5.);
342 for(j=-10; j<=10; ++j) {
344 gg->SetPoint(j+10,
t+j/2.,Tp(&
res,
par));
346 if(gg->GetFunction(
"ff")) gg->GetFunction(
"ff")->Delete();
347 gg->Fit(
"ff",
"QRI0");
348 res = 2*
t*
ff->GetParameter(2) +
ff->GetParameter(1);
354 Double_t Tp(Double_t *
tt, Double_t *
par)
359 if(*
tt < t0 +
par[0])
return 0;
367 sp = 1/taup; s0 = 1/tau0;
sz = 1/tauz;
ss = 1/taus;
sc = 1/tc;
372 rc += al*(
sc*s0+
sc*sp+sp*s0) - sp*s0;
375 rz =
sz*
sz*(
sc+s0+sp) -
pow(sz,3) -
sz*
sc*(s0+sp) + sp*s0*(sc-sz);
378 rs = -
ss*
ss*(
sc+s0+sp-
ss) + sc*ss*(s0+sp) +
ss*sp*s0 -
sc*sp*s0;
391 Double_t Tp_gsl(Double_t *
tt, Double_t *
par)
396 if(*
tt < t0 +
par[0])
return 0;
403 sp = 1/taup; s0 = 1/tau0;
sz = 1/tauz;
ss = 1/taus;
sc = 1/tc;
408 rc += al*(
sc*s0+
sc*sp+sp*s0) - sp*s0;
411 rz =
sz*
sz*(
sc+s0+sp) -
pow(sz,3) -
sz*
sc*(s0+sp) + sp*s0*(sc-sz);
414 rs = -
ss*
ss*(
sc+s0+sp-
ss) + sc*ss*(s0+sp) +
ss*sp*s0 -
sc*sp*s0;
421 t4 = Tp4_gsl(
t);
t5 = Tp5_gsl(
t);
429 Double_t Tp4_gsl(Double_t
t)
437 static const double omega=0.0;
438 static const size_t n=100;
439 static const size_t limit=1001;
440 static const double epsabs=0.0;
441 static double epsrel=1
e-10;
443 static gsl_integration_workspace *
w=0;
444 static gsl_integration_qawo_table *wf=0;
448 gsl_function gsl_func;
456 if(!
w)
w = gsl_integration_workspace_alloc(
limit);
457 if(!wf) wf=gsl_integration_qawo_table_alloc(omega,L,GSL_INTEG_COSINE,
n);
459 gsl_integration_qawo_table_set_length(wf, L);
462 gsl_func.function = &f4_gsl;
463 gsl_func.params = &
params;
467 gsl_set_error_handler_off();
469 if(
b>t0+130) epsrel=1
e-8;
472 unsigned short int INTEGRAL_ITER=0;
473 const unsigned short int INTEGRAL_ITER_MAX=10;
477 while((
status!=0 || INTEGRAL_ITER==0) && INTEGRAL_ITER<INTEGRAL_ITER_MAX){
481 if(
status == GSL_EROUND) epsrel*=10;
485 std::cout<<
"WARNING Tp4_gsl:QAWO ERROR (not roundoff): status="<<gsl_strerror(
status)<<
" ERRcode="<<
status <<
" LArWaveFlag set to -1 !!"<< std::endl ;
493 if(INTEGRAL_ITER==INTEGRAL_ITER_MAX &&
status!=0){
494 std::cout<<
"WARNING Integration of Tp4_gsl FAILED !!!! #iterations reached the MAXIMUM!!! status="<<gsl_strerror(
status)<<
" ERRcode="<<
status <<
" LArWaveFlag set to -1 !!"<<std::endl ;
504 static double f4_gsl(Double_t
x,
void *
par)
512 return Tc(
x) * Rd(
t -
x);
517 Double_t Tp5_gsl(Double_t
t)
519 if(
t<t0 + tdr)
return 0.;
526 static const double omega=0.0;
527 static const size_t n=100;
528 static const size_t limit=1001;
529 static const double epsabs=0.0;
530 static double epsrel = 1
e-10;
532 static gsl_integration_workspace *
w=0;
533 static gsl_integration_qawo_table *wf=0;
537 gsl_function gsl_func;
543 if(!
w)
w = gsl_integration_workspace_alloc(
limit);
544 if(!wf) wf=gsl_integration_qawo_table_alloc(omega,L,GSL_INTEG_COSINE,
n);
546 gsl_integration_qawo_table_set_length(wf, L);
549 gsl_func.function = &f5_gsl;
550 gsl_func.params = &
params;
555 gsl_set_error_handler_off();
557 if(
b>t0+60) epsrel=1
e-8;
561 unsigned short int INTEGRAL_ITER=0;
562 const unsigned short int INTEGRAL_ITER_MAX=10;
565 while((
status!=0 || INTEGRAL_ITER==0) && INTEGRAL_ITER<INTEGRAL_ITER_MAX){
568 if(
status == GSL_EROUND) epsrel*=10;
572 std::cout <<
"WARNING: Tp5_gsl:QAWO ERROR (not roundoff): status="<<gsl_strerror(
status)<<
" ERRcode="<<
status <<
" LArWaveFlag set to -1 !!"<< std::endl ;
580 if(INTEGRAL_ITER==INTEGRAL_ITER_MAX &&
status!=0){
581 std::cout<<
"WARNING: Integration of Tp5_gsl FAILED: #iterations reached the MAXIMUM!!!! status="<<gsl_strerror(
status)<<
" ERRcode="<<
status <<
" LArWaveFlag set to -1 !!"<<std::endl ;
594 static double f5_gsl(Double_t
x,
void *
par)
601 return Tc(
x) * Ro(
t -
x);
607 Double_t Tp4(Double_t
t)
610 static const Double_t epsrel = 0.001;
611 static TF1 *fun4 = NULL;
613 if(fun4 == NULL) fun4 =
new TF1(
"fun4",&f4_g,0.,
NMAX,1);
619 fun4->SetParameter(0,
t);
620 return fun4->Integral(
a,
b,epsrel);
625 static Double_t f4_g(Double_t *
x, Double_t *
par)
631 return Tc(*
x) * Rd(
t - *
x);
635 Double_t Tp5(Double_t
t)
638 static const Double_t
a = 0., epsrel = 0.001;
639 static TF1 *fun5 = NULL;
641 if(fun5 == NULL) fun5 =
new TF1(
"fun5",&f5_g,0.,
NMAX,1);
644 fun5->SetParameter(0,
t);
645 return fun5->Integral(
a,
b,epsrel);
651 static Double_t f5_g(Double_t *
x, Double_t *
par)
657 return Tc(*
x) * Ro(
t - *
x);
663 Bool_t uset0, Bool_t
norm, Int_t , Double_t *,
bool gsl_flag){
665 TF1 *fitfun{}, *fu{};
666 Int_t
i{},nbin{},nbino{},i5bin{},
count{} ;
667 Double_t pampl{},
t5{}, mcut{}, rshift{}, rmult{};
689 nbino = nbin = pcal->GetNbinsX();
697 mcut = 0.05*pcal->GetMaximum();
700 cout<<
"mcut: "<<mcut<<endl;
702 for(
i=1;
i<=nbin; ++
i) {
703 x[
i-1] = pcal->GetBinCenter(
i);
704 y[
i-1] = pcal->GetBinContent(
i);
705 if(i5bin == 0 &&
y[
i-1] > mcut) i5bin =
i;
709 printf(
"\n number of bins too small !\n");
712 t5 = (mcut-
y[i5bin-2])*(
x[i5bin-1]-
x[i5bin-2])/(
y[i5bin-1]-
y[i5bin-2]);
720 if(fabs((
x[1] -
x[0]) - (
x[3] -
x[2])) > 0.000001 ) {
721 cout<<
x[0] <<
" " <<
x[1] <<
" " <<
x[2] <<
" " <<
x[3] << endl;
722 printf(
"\n Nonuniform step !\n");
729 if(pref)
delete [] pref;
730 pref =
new Double_t[nbin];
733 for(
i=0;
i<nbin; ++
i) {
734 if(
x[
i] < t0) {pref[
i] = 0.;
continue;}
740 if(
i<nbin) nbin =
i+1;
743 if(splin)
delete splin;
744 splin =
new TSpline3(
"sp",
x,pref,nbin,
"b1",0.);
745 fu =
new TF1(
"fu",P,
x[15],
x[nbin-15],2);
746 fu->SetParName(0,
"t0 shift");
747 fu->SetParLimits(0,-15.,15.);
748 fu->SetParName(1,
"mult");
749 fu->SetParLimits(1,0.,10000.);
750 fu->SetParameters(5.,5.);
752 if( pcal->Fit(
"fu",
"RQ0") == 0) {
753 rshift = pcal->GetFunction(
"fu")->GetParameter(0);
754 rmult = pcal->GetFunction(
"fu")->GetParameter(1);
762 for(
i=0;
i<nbin; ++
i) {
763 if(
x[
i] < t0-rshift) {pref[
i] = 0.;
continue;}
769 if(
i<nbin) nbin =
i+1;
772 for(
i=1;
i<=nbin; ++
i)
x[
i-1] -= t0;
776 for(
i=1;
i<=nbin; ++
i)
x[
i-1] +=
par[0];
781 if(splin)
delete splin;
782 splin =
new TSpline3(
"sp",
x,pref,nbin,
"b1",0.);
787 fitfun =
new TF1(
"fitfun",&Tp,
x[0],
x[nbin-1],3);
789 fitfun =
new TF1(
"fitfun",&Tp_gsl,
x[0],
x[nbin-1],3);
791 fitfun->SetParName(0,
"t0 shift");
792 fitfun->SetParLimits(0,0.,100.);
793 fitfun->SetParName(1,
"mult");
794 fitfun->SetParLimits(1,0.,100.);
795 fitfun->SetParName(2,
"Taus");
796 fitfun->SetParLimits(2,0.1,10.);
797 fitfun->SetNpx(nbin+1);
798 fitfun->SetParameters(
par[0],
par[1],
par[2]);
807 deriv =
new TF1(
"derivfun",&DTp,
x[0],
x[nbin-1],3);
808 deriv->SetParName(0,
"t0 shift");
809 deriv->SetParLimits(0,0.,100.);
810 deriv->SetParName(1,
"mult");
811 deriv->SetParLimits(1,0.,100.);
812 deriv->SetParName(2,
"Taus");
813 deriv->SetParLimits(2,0.1,10.);
814 deriv->SetNpx(nbin+1);
815 deriv->SetParameters(
par[0],
par[1],
par[2]);
826 Double_t from, Double_t
to, Double_t
step)
828 static Double_t amax, tmax,
t, ampl;
829 static Double_t
t1,
t2,tt1,tt2;
832 if((from ==
to) || (from >
to))
833 func->GetRange(
t1,
t2);
835 func->GetRange(tt1,tt2);
836 if(from>tt1)
t1 = from;
else t1 = tt1;
843 if(func->Eval(
t) > amax) {
844 amax = func->Eval(
t); tmax =
t;
847 Double_t
x[20],
y[20];
849 for(
t=tmax-5;
t<tmax+5; ++
t) {
851 y[
i] = func->Eval(
t);
854 TGraph *
g =
new TGraph(10,
x,
y);
855 g->Fit(
"pol2",
"QIE0");
856 tmax = -
g->GetFunction(
"pol2")->GetParameter(1) / (2 *
g->GetFunction(
"pol2")->GetParameter(2));
857 amax =
g->GetFunction(
"pol2")->Eval(tmax);
859 ampl = func->GetParameter(1);
860 if(amax != 0) func->SetParameter(1,ampl/amax);
862 cout <<
"Normalize: tmax="<<tmax <<
" ampl=" << ampl<<
" amax="<<amax << endl;
863 if(rampl != NULL) *rampl = amax;
870 Double_t amax, tmax,
x,
y,
err;
872 amax =
histo->GetBinContent(0);
873 tmax =
histo->GetBinCenter(0);
877 if(amax < histo->GetBinContent(
i))
879 amax =
histo->GetBinContent(
i);
880 tmax =
histo->GetBinCenter(
i);
884 histo->Fit(
"pol2",
"0QI",
"",tmax-10., tmax+10.);
885 tmax = -
histo->GetFunction(
"pol2")->GetParameter(1) / (2 *
histo->GetFunction(
"pol2")->GetParameter(2));
886 amax =
histo->GetFunction(
"pol2")->Eval(tmax);
891 if(htmp->GetBinEntries(
i) == 0)
continue;
892 x = htmp->GetBinCenter(
i);
893 y = htmp->GetBinContent(
i);
894 nent=
int(htmp->GetBinEntries(
i));
895 err = htmp->GetBinError(
i);
908 if(rampl != NULL) *rampl = amax;
915 TGraph *
g=
new TGraph(
ymax,
x,
y);
916 TGraphSmooth *gs =
new TGraphSmooth();
917 TGraph *
g1 = gs->SmoothKern(
g,
"normal",band);
918 xx =
g1->GetX();
yy =
g1->GetY();
919 for(Int_t
i=0;
i<
g->GetN(); ++
i) {
930 static Double_t *dif = NULL;
933 if(dif)
delete[] dif;
934 dif =
new Double_t[
ymax];