41 MERGE_THRESH_DIST = 40.0;
42 MERGE_THRESH_TIME = 2.5;
43 MERGE_THRESH_ZDIST = 30.0;
44 MERGE_THRESH_XYDIST = 40.0;
47 BREAK_THRESH_TRMS= 5.0;
49 gPARMS->SetDefaultParameter(
"BCALRECON:CLUST_THRESH", CLUST_THRESH );
50 gPARMS->SetDefaultParameter(
"BCALRECON:MERGE_THRESH_DIST", MERGE_THRESH_DIST );
51 gPARMS->SetDefaultParameter(
"BCALRECON:MERGE_THRESH_TIME", MERGE_THRESH_TIME );
52 gPARMS->SetDefaultParameter(
"BCALRECON:MERGE_THRESH_ZDIST", MERGE_THRESH_ZDIST );
53 gPARMS->SetDefaultParameter(
"BCALRECON:MERGE_THRESH_XYDIST", MERGE_THRESH_XYDIST );
54 gPARMS->SetDefaultParameter(
"BCALRECON:BREAK_THRESH_TRMS", BREAK_THRESH_TRMS );
78 m_scaleZ_p0 = 0.99798;
79 m_scaleZ_p1 = 0.000361096;
80 m_scaleZ_p2 = -2.17338e-06;
81 m_scaleZ_p3 = 1.32201e-09;
83 m_nonlinZ_p0 = -0.0201272;
84 m_nonlinZ_p1 = 0.000103649;
102 vector<const DBCALGeometry*> bcalGeomVect;
103 loop->Get( bcalGeomVect );
104 bcalGeom = bcalGeomVect[0];
110 ATTEN_LENGTH = bcalGeom->GetBCAL_attenutation_length();
111 C_EFFECTIVE = bcalGeom->GetBCAL_c_effective();
113 fiberLength = bcalGeom->GetBCAL_length();
114 zOffset = bcalGeom->GetBCAL_center();
118 int modmax = bcalGeom->GetBCAL_Nmodules();
120 int rowmax1= bcalGeom->GetBCAL_NInnerLayers();
121 int rowmin2= rowmax1;
122 int rowmax2= bcalGeom->GetBCAL_NOuterLayers()+rowmin2;
124 int colmax1=bcalGeom->GetBCAL_NInnerSectors();
126 int colmax2=bcalGeom->GetBCAL_NOuterSectors();
128 float r_inner= bcalGeom->GetBCAL_inner_rad();
130 for (
int i = (rowmin1+1); i < (rowmax1+1); i++){
132 int cellId = bcalGeom->cellId(1,i,1);
134 rt[i]=bcalGeom->r(cellId)-r_inner;
137 for (
int i = (rowmin2+1); i < (rowmax2+1); i++){
138 int cellId = bcalGeom->cellId(1,i,1);
139 rt[i]=bcalGeom->r(cellId)-r_inner;
147 for (
int k = modmin; k < modmax; k++){
148 for (
int i = rowmin1; i < rowmax1; i++){
149 for (
int j = colmin1; j < colmax1; j++){
152 int cellId = bcalGeom->cellId(k+1,i+1,j+1);
153 r[k][i][j]=bcalGeom->r(cellId);
154 phi[k][i][j]=bcalGeom->phi(cellId);
156 xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]);
157 yy[k][i][j]=r[k][i][j]*
sin(phi[k][i][j]);
161 for (
int i = rowmin2; i < rowmax2; i++){
162 for (
int j = colmin2; j < colmax2; j++){
163 int cellId = bcalGeom->cellId(k+1,i+1,j+1);
164 r[k][i][j]=bcalGeom->r(cellId);
165 phi[k][i][j]=bcalGeom->phi(cellId);
167 xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]);
168 yy[k][i][j]=r[k][i][j]*
sin(phi[k][i][j]);
195 vector<DBCALShower*> clusters;
197 for (
int i = 1; i < (clstot+1); i++){
201 if( e_cls[j] < CLUST_THRESH )
continue;
234 vector<const DBCALPoint*> pointsInShower;
235 FindPointsInShower(j, loop, pointsInShower);
244 double E=0,
x=0,
y=0,z=0,t=0;
246 double sig_x=0,sig_y=0,sig_z=0,sig_t=0;
248 for(
unsigned int j=0; j<pointsInShower.size(); j++){
249 double cell_E = pointsInShower[j]->E();
250 double cell_r = pointsInShower[j]->r();
251 double cell_phi = pointsInShower[j]->phi();
253 x += cell_E*cell_r*cos(cell_phi);
254 sig_x += cell_E*cell_r*cos(cell_phi)*cell_r*cos(cell_phi);
255 y += cell_E*cell_r*
sin(cell_phi);
256 sig_y += cell_E*cell_r*
sin(cell_phi)*cell_r*
sin(cell_phi);
259 double z_wt = 1/(pointsInShower[j]->sigZ()*pointsInShower[j]->sigZ());
260 double cell_z = pointsInShower[j]->z();
261 double cell_t = pointsInShower[j]->t();
265 sig_z += z_wt*cell_z*cell_z;
268 sig_t += cell_E*cell_t*cell_t;
271 shower->AddAssociatedObject(pointsInShower[j]);
282 sig_z =
sqrt(sig_z - z*z)/
sqrt(N_cell);
285 sig_t =
sqrt(sig_t - t*t)/
sqrt(N_cell);
291 shower->
z = z + m_z_target_center;
307 float r =
sqrt( shower->
x * shower->
x + shower->
y * shower->
y );
309 float zEntry = ( shower->
z - m_z_target_center ) * ( bcalGeom->GetBCAL_inner_rad() / r );
311 float scale = m_scaleZ_p0 + m_scaleZ_p1*zEntry +
312 m_scaleZ_p2*(zEntry*zEntry) + m_scaleZ_p3*(zEntry*zEntry*zEntry);
313 float nonlin = m_nonlinZ_p0 + m_nonlinZ_p1*zEntry +
314 m_nonlinZ_p2*(zEntry*zEntry) + m_nonlinZ_p3*(zEntry*zEntry*zEntry);
316 shower->
E = pow( (shower->
E_raw ) / scale, 1 / ( 1 + nonlin ) );
324 _data.push_back(shower);
352 vector<const DBCALPoint*> points;
355 int start_indx = indx;
357 int module = narr[1][indx];
358 int layer = narr[2][indx];
359 int sector = narr[3][indx];
362 for(
unsigned int i=0; i<points.size(); i++){
363 if(points[i]->module() !=module)
continue;
365 if(points[i]->sector() !=sector)
continue;
366 pointsInShower.push_back(points[i]);
371 }
while(indx != start_indx);
416 vector<const DBCALPoint*> points;
418 if(points.size() <=0)
return;
420 for (vector<const DBCALPoint*>::const_iterator point_iter = points.begin();
421 point_iter != points.end();
424 int module = point.
module();
426 int sector = point.
sector();
427 double r = point.
r();
428 double phi = point.
phi();
429 double x = r*cos(phi);
430 double y = r*
sin(phi);
432 xcel[module-1][layer-1][sector-1] =
x;
433 ycel[module-1][layer-1][sector-1] =
y;
436 zcel[module-1][layer-1][sector-1] = point.
z()+m_z_target_center-zOffset;
437 tcel[module-1][layer-1][sector-1] = point.
t();
438 ecel[module-1][layer-1][sector-1] = point.
E();
442 double EUp=0,EDown=0,tUp=0,tDown=0;
443 vector<const DBCALUnifiedHit*> assoc_hits;
444 point.Get(assoc_hits);
445 for (
unsigned int i=0; i<assoc_hits.size(); i++) {
447 EUp = assoc_hits[i]->E;
448 tUp = assoc_hits[i]->t;
451 EDown = assoc_hits[i]->E;
452 tDown = assoc_hits[i]->t;
457 ecel_a[module-1][layer-1][sector-1] = EUp;
459 tcel_a[module-1][layer-1][sector-1] = tUp;
460 tcell_anor[module-1][layer-1][sector-1] = tUp;
462 ecel_b[module-1][layer-1][sector-1] = EDown;
463 tcel_b[module-1][layer-1][sector-1] = tDown;
464 tcell_bnor[module-1][layer-1][sector-1] = tDown;
516 float ea = ecel_a[k][i][j];
517 float eb = ecel_b[k][i][j];
518 float ta = tcel_a[k][i][j];
519 float tb = tcel_b[k][i][j];
521 if( (
min(ea,eb)>ethr_cell) & (fabs(ta-tb)<35.) & (ta!=0.) & (tb!=0.)) {
538 celdata[1][celtot]=ea/0.145;
539 celdata[2][celtot]=eb/0.145;
541 nclus[celtot] = celtot;
542 next[celtot] = celtot;
544 e_cel[celtot] = ecel[k][i][j];
545 x_cel[celtot] = xcel[k][i][j];
546 y_cel[celtot] = ycel[k][i][j];
547 z_cel[celtot] = zcel[k][i][j];
548 t_cel[celtot] = tcel[k][i][j];
550 ta_cel[celtot]=tcell_anor[k][i][j];
551 tb_cel[celtot]=tcell_bnor[k][i][j];
582 int modmax = bcalGeom->GetBCAL_Nmodules();
586 int rowmax1= bcalGeom->GetBCAL_NInnerLayers();
587 int rowmin2= rowmax1+1;
589 int colmax1=bcalGeom->GetBCAL_NInnerSectors();
590 int colmax2=bcalGeom->GetBCAL_NOuterSectors();
592 float r_middle= bcalGeom->GetBCAL_middle_rad();
595 float thick_inner=bcalGeom->rSize(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers(),1));
597 float thick_outer=bcalGeom->rSize(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers()+1,1));
600 float dis_in_out=bcalGeom->r(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers()+1,1))-bcalGeom->r(bcalGeom->cellId(1,bcalGeom->GetBCAL_NInnerLayers(),1));
602 float degree_permodule=360.0/(modmax-modmin);
603 float half_degree_permodule=degree_permodule/2.0;
606 float width_1=2.0*(r_middle-thick_inner/2.0)*
607 sin(half_degree_permodule*3.141593/180)/colmax1;
609 float width_2=2.0*(r_middle+thick_outer/2.0)*
610 sin(half_degree_permodule*3.141593/180)/colmax2;
614 float disthres=width_2*1.5-width_1*0.5+0.0001;
616 for (
int i = 1; i < (celtot+1); i++){
622 for (
int j = 1; j < (celtot+1); j++){
623 if ( (j!=i) & (nclus[j]!=nclus[i]) & (e_cel[j]>emin)) {
632 int amodif = abs(modiff);
635 if ( (abs(i1-i2)<=k) & ((amodif<=1) || (amodif==47)) ) {
642 if ( (i1<=rowmax1) & (i2<=rowmax1) & (abs(j2-j1)<=k) ) {
649 if ( (i1>=rowmin2) & (i2>=rowmin2) & (abs(j2-j1)<=k) ) {
656 if( (modiff==1) || (modiff==-47) ) {
657 if ( (i1<=rowmax1) & (i2<=rowmax1) ){
658 if(abs((j1+colmax1)-j2)<=k){
664 if ( (i1>=rowmin2) & (i2>=rowmin2) ) {
665 if(abs((j1+colmax2)-j2)<=k){
672 if ( (modiff==-1) || (modiff==47) ) {
674 if ( (i1<=rowmax1) & (i2<=rowmax1) ){
675 if(abs((j2+colmax1)-j1)<=k){
681 if ( (i1>=rowmin2) & (i2>=rowmin2) ){
682 if(abs((j2+colmax2)-j1)<=k){
693 if( ( (i1 == rowmax1) & (i2 == rowmin2) ) ||
694 ( (i1 == rowmin2) & (i2 == rowmax1) ) ) {
696 float delta_xx=xx[k1-1][i1-1][j1-1]-xx[k2-1][i2-1][j2-1];
697 float delta_yy=yy[k1-1][i1-1][j1-1]-yy[k2-1][i2-1][j2-1];
700 float dis =
sqrt( delta_xx * delta_xx + delta_yy * delta_yy );
702 dis =
sqrt( dis*dis - dis_in_out * dis_in_out );
704 if( dis < disthres ){
760 if(nclus[n]!=nclus[m]){
779 memset( e_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
780 memset( x_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
781 memset( y_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
782 memset( z_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
783 memset( t_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
784 memset( ea_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
785 memset( eb_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
786 memset( ta_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
787 memset( tb_cls, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
788 memset( tsqr_a, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
789 memset( tsqr_b, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
790 memset( trms_a, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
791 memset( trms_b, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
792 memset( e2_a, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
793 memset( e2_b, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
794 memset( clspoi, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
795 memset( ncltot, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
796 memset( ntopol, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
809 for (
int ix = 1; ix < (celtot+1); ix++){
818 for (
int i = 1; i < (clstot+1); i++){
819 if(n==clspoi[i]) j=i;
830 if(e_cel[ix]<0.000000001)
continue;
832 x_cls[n]=(e_cls[n]*x_cls[n]+e_cel[ix]*x_cel[ix])
833 /(e_cls[n]+e_cel[ix]);
835 y_cls[n]=(e_cls[n]*y_cls[n]+e_cel[ix]*y_cel[ix])
836 /(e_cls[n]+e_cel[ix]);
838 z_cls[n]=(e_cls[n]*z_cls[n]+e_cel[ix]*z_cel[ix])
839 /(e_cls[n]+e_cel[ix]);
841 t_cls[n]=(e_cls[n]*t_cls[n]+e_cel[ix]*t_cel[ix])
842 /(e_cls[n]+e_cel[ix]);
844 e_cls[n]=e_cls[n]+e_cel[ix];
851 ta_cls[n]=(ea_cls[n]*ta_cls[n]+celdata[1][ix]*ta_cel[ix])
852 /(ea_cls[n]+celdata[1][ix]);
853 tsqr_a[n]=(ea_cls[n]*tsqr_a[n]+celdata[1][ix]*ta_cel[ix]*
854 ta_cel[ix])/(ea_cls[n]+celdata[1][ix]);
855 ea_cls[n]=ea_cls[n]+celdata[1][ix];
856 e2_a[n]=e2_a[n]+celdata[1][ix]*celdata[1][ix];
857 tb_cls[n]=(eb_cls[n]*tb_cls[n]+celdata[2][ix]*tb_cel[ix])/
858 (eb_cls[n]+celdata[2][ix]);
859 tsqr_b[n]=(eb_cls[n]*tsqr_b[n]+celdata[2][ix]*tb_cel[ix]*
860 tb_cel[ix])/(eb_cls[n]+celdata[2][ix]);
861 eb_cls[n]=eb_cls[n]+celdata[2][ix];
862 e2_b[n]=e2_b[n]+celdata[2][ix]*celdata[2][ix];
870 if( narr[1][n] != narr[1][ix] || narr[2][n] != narr[2][ix] )
877 for (
int n = 1; n < (clstot+1); n++){
880 if( ncltot[ix] > 1) {
882 float effnum = ea_cls[ix] * ea_cls[ix] / e2_a[ix];
883 trms_a[ix] = effnum / ( effnum - 1 ) *
884 ( tsqr_a[ix] - ta_cls[ix] * ta_cls[ix] );
886 effnum = eb_cls[ix] * eb_cls[ix] / e2_b[ix];
887 trms_b[ix] = effnum / ( effnum - 1 ) *
888 ( tsqr_b[ix] - tb_cls[ix] * tb_cls[ix] );
890 if( trms_a[ix] <= 0.0 ) trms_a[ix] = 0.;
891 if( trms_b[ix] <= 0.0 ) trms_b[ix] = 0.;
892 trms_a[ix] =
sqrt( trms_a[ix] );
893 trms_b[ix] =
sqrt( trms_b[ix] );
909 bool newClust =
false;
915 for (
int i = 0; i < 2; i++){
916 for (
int j = 1; j < (clstot+1); j++){
919 float dist=
sqrt(trms_a[ix]*trms_a[ix]+trms_b[ix]*trms_b[ix]);
920 if(dist>BREAK_THRESH_TRMS) {
939 for (
int i = 1; i < clstot; i++){
942 for (
int j = (i+1); j < (clstot+1); j++){
948 if ( (e_cls[ix]>0.0) & (e_cls[iy]>0.0) ) {
950 float delta_x=x_cls[ix]-x_cls[iy];
951 float delta_y=y_cls[ix]-y_cls[iy];
952 float delta_z=z_cls[ix]-z_cls[iy];
953 float dist=
sqrt(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z);
955 float tdif=fabs(t_cls[ix]-t_cls[iy]);
960 if ( (dist<MERGE_THRESH_DIST) & (tdif<MERGE_THRESH_TIME) ){
961 float zdif=fabs(z_cls[ix]-z_cls[iy]);
962 float distran=
sqrt(delta_x*delta_x+delta_y*delta_y);
964 if ( (zdif<MERGE_THRESH_ZDIST) & (distran<MERGE_THRESH_XYDIST) ){
965 if(e_cls[ix]>=e_cls[iy]) {
978 if(
min(icls[1],icls[2])>0){
980 Connect(icls[1],icls[2]);
998 float tdif,tdif_a,tdif_b,tseed[5];
1001 for (
int i =0; i < 5; i++){
1009 tdif_a=ta_cel[n]-ta_cls[nclust];
1010 tdif_b=tb_cel[n]-tb_cls[nclust];
1035 float tdif=
sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
1036 if(tdif>tseed[selnum]){
1045 while(next[n]!=nclust) {
1047 tdif_a=ta_cel[n]-ta_cls[nclust];
1048 tdif_b=tb_cel[n]-tb_cls[nclust];
1078 tdif=
sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
1080 if(tdif>tseed[selnum]){
1097 for (
int i =1; i < 5; i++){
1101 nclus[nseed[i]]=nseed[i];
1102 next[nseed[i]]=nseed[i];
1104 for (
int j =1; j < (celtot+1); j++){
1105 if ( (nclus[j]==nclust) & (j!=nseed[i]) ){
1109 Connect(nseed[i],j);
1133 memset( apx, 0, (
clsmax_bcal + 1 ) * 4 *
sizeof(
float ) );
1134 memset( eapx, 0, (
clsmax_bcal + 1 ) * 4 *
sizeof(
float ) );
1135 memset( ctrk, 0, (
clsmax_bcal + 1 ) * 4 *
sizeof(
float ) );
1136 memset( ectrk, 0, (
clsmax_bcal + 1 ) * 4 *
sizeof(
float ) );
1138 for (
int ix = 1; ix < (celtot+1); ix++){
1146 int lyr = narr[2][ix];
1149 clslyr[xlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[xlyr][lyr][n]
1150 +e_cel[ix]*x_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1152 clslyr[ylyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[ylyr][lyr][n]
1153 +e_cel[ix]*y_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1156 clslyr[zlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[zlyr][lyr][n]
1157 +e_cel[ix]*z_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1160 clslyr[tlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[tlyr][lyr][n]
1161 +e_cel[ix]*t_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
1164 clslyr[elyr][lyr][n]=clslyr[elyr][lyr][n]+e_cel[ix];
1169 memset( nlrtot, 0, (
clsmax_bcal + 1 ) *
sizeof(
float ) );
1171 for (
int n = 1; n < ( clstot + 1 ); n++){
1177 if( clslyr[elyr][i][ix] > 0.0 ) nlrtot[ix]++;
1195 if(clslyr[elyr][il][ix]>emin) {
1198 x[nltot]= clslyr[xlyr][il][ix];
1199 y[nltot]= clslyr[ylyr][il][ix];
1200 z[nltot]= clslyr[zlyr][il][ix];
1201 e[nltot]= clslyr[elyr][il][ix];
1203 sigy[nltot] = 1.0/
e[nltot];
1204 sigx[nltot] = 1.0/
e[nltot];
1205 sigz[nltot] = 1.0/
sqrt(
e[nltot]);
1251 for (
int i = 1; i < 4; i++){
1253 ctrk[i][ix]=ctrk_ix[i];
1254 ectrk[i][ix]=ectrk_ix[i];
1255 apx[i][ix]=apx_ix[i];
1256 eapx[i][ix]=eapx_ix[i];
1266 eapx[1][ix] = sigx[1];
1267 eapx[2][ix] = sigy[1];
1268 eapx[3][ix] = sigz[1];
1272 ctrk[1][ix] = 999.0;
1273 ctrk[2][ix] = 999.0;
1274 ctrk[3][ix] = 999.0;
1285 float d,
e,
f,chi2,q,norm;
1286 float siga,sigb,sigc,sigd,sige,sigf;
1287 float sigb2,sigd2,sigf2;
1290 Linefit(1,1,a,b,siga,sigb,chi2,q);
1292 Linefit(2,1,c,d,sigc,sigd,chi2,q);
1294 Linefit(3,1,e,f,sige,sigf,chi2,q);
1311 norm=
sqrt(b*b+d*d+f*f);
1317 float norm3=norm*norm*norm;
1319 ectrk_ix[1]=
sqrt((d*d+f*f)*(d*d+f*f)*sigb2+b*b*d*d*sigd2+b*b*f*f*sigf2)/norm3;
1320 ectrk_ix[2]=
sqrt((b*b+f*f)*(b*b+f*f)*sigd2+d*d*b*b*sigb2+d*d*f*f*sigf2)/norm3;
1321 ectrk_ix[3]=
sqrt((b*b+d*d)*(b*b+d*d)*sigf2+f*f*b*b*sigb2+f*f*d*d*sigd2)/norm3;
1331 float &b,
float &siga,
float &sigb,
float &chi2,
float &q)
1353 float sigdat,ss,st2,sx,sxoss,sy,t,wt;
1369 if(etemp>0.0001)ndata=ndata+1;
1379 if(etemp>0.000001)ndata=ndata+1;
1388 if(etemp>0.000001)ndata=ndata+1;
1394 for (
int i = 1; i < (ndata+1); i++){
1395 wt=1.0/(sig[i]*sig[i]);
1403 for (
int i = 1; i < (ndata+1); i++){
1413 for (
int i = 1; i < (ndata+1); i++){
1414 t=(xtemp[i]-sxoss)/sig[i];
1416 b=b+t*ytemp[i]/sig[i];
1423 for (
int i = 1; i < (ndata+1); i++){
1432 siga=
sqrt((1.0+sx*sx/(ss*st2))/ss);
1437 for (
int i = 1; i < (ndata+1); i++){
1438 chi2=chi2+(ytemp[i]-a-b*xtemp[i])*(ytemp[i]-a-b*xtemp[i]);
1441 sigdat=
sqrt(chi2/(ndata-2));
1446 for (
int i = 1; i < (ndata+1); i++){
1447 chi2=chi2+((ytemp[i]-a-b*xtemp[i])/
1448 sig[i])*((ytemp[i]-a-b*xtemp[i])/sig[i]);
1450 q=Gammq(0.5*(ndata-2),0.5*chi2);
1472 float gammcf=0,gamser;
1481 if(x_gammq<0. || a_gammq<= 0.0) {
1486 if(x_gammq<(a_gammq+1.)) {
1487 Gser(gamser,a_gammq,x_gammq);
1491 Gcf(gammcf,a_gammq,x_gammq);
1517 if(x_gser<0.0) cout<<
"x_gser<0 in gser"<<
"\n";
1527 for (
int n = 1; n < (itmax+1); n++){
1532 if(fabs(del)<fabs(sum)*eps) {
1533 gamser=sum*exp(-x_gser+a_gser*log(x_gser)-gln);
1553 float fpmin=1.0e-30;
1566 float an,b,
c,d,del,
h;
1575 for (
int i = 1; i < (itmax+1); i++){
1579 if(fabs(d)<fpmin)d=fpmin;
1581 if(fabs(c)<fpmin)c=fpmin;
1585 if(fabs(del-1.0)<eps) {
1586 gammcf=exp(-x_gcf+a_gcf*log(x_gcf)-gln)*
h;
1601 float ser,stp,tmp,x_gln,y_gln;
1612 stp=2.5066282746310005;
1613 cof[1]=76.18009172947146;
1614 cof[2]=-86.50532032941677;
1615 cof[3]=24.01409824083091;
1616 cof[4]=-1.231739572450155;
1617 cof[5]=.1208650973866179e-2;
1618 cof[6]=-.5395239384953e-5;
1623 tmp=(x_gln+0.5)*log(tmp)-tmp;
1625 ser=1.000000000190015;
1627 for (
int j = 1; j < 7; j++){
1629 ser=ser+cof[j]/y_gln;
1633 gammln=tmp+log(stp*ser/x_gln);
void FindPointsInShower(int indx, JEventLoop *loop, vector< const DBCALPoint * > &pointsInShower)
void Gser(float &gamser, float a, float x)
float Gammq(float a, float x)
void Gcf(float &gammcf, float a, float x)
void PreCluster(JEventLoop *loop)
void Clus_Break(int nclust)
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via JEventProcessor virtual method.
DGeometry * GetDGeometry(unsigned int run_number)
void Linefit(int ixyz, int mwt, float &a, float &b, float &siga, float &sigb, float &chi2, float &q)
float Gammln(float xx_gln)
void CellRecon(JEventLoop *loop)
bool GetTargetZ(double &z_target) const
z-location of center of target
jerror_t brun(JEventLoop *loop, int32_t runnumber)