#include "hits.h" // return absolute value, if that value is >1, // otherwise return 1/absolute value inline double hits::absrat(double el) { double ael=fabs(el); if (ael>1) return(ael); if (ael<1e-10) return(1e10); return(1/ael); } // ********************************************** // find ``smallest'' matrix element ``above'' // (pass,pass) // ********************************************** inline double hits::findmin(double *mat,int *rowind,int *colind,int pass) { int row,col,ind; int saverow=rowind[pass],savecol=colind[pass]; double test,min=absrat(mat[5*saverow+savecol]); for(row=pass; row<3; row++) { ind=5*rowind[row]; for(col=pass; col<3; col++) if ((row!=pass) || (col!=pass)) if ((test=absrat(mat[ind+colind[col]])) =1e10) && (col<5); col++) { ind=colind[col]; for(row=pass; row<3; row++) if ((test=absrat(mat[5*rowind[row]+ind])) 0; nfirst--) list[firstp++]=*(first++); } // ********************************************** // * create space for all hits * // ********************************************** inline void hits::create_arrays(int n_raw,int n_hit) { places=new short int[n_raw]; cables=new int[n_hit]; times=new float[n_hit]; charges=new float[n_hit]; positions=new float[3*n_hit]; } // ********************************************** // ********************************************** // public functions // ********************************************** // constructor: // Stores hit times, charges, cable numbers, and // positions // ********************************************** hits::hits(int ns,float *set,float *pmt_loc,int *bad_ch, int n_raw,int *cable,float *tim_raw,float *chg_raw) { int i,j,number,pos,*nh2; // count #hits per charge bin and create offset 'pointers' nchargeset=ns; nhit=new int[nchargeset]; avcharges=new float[nchargeset]; for(j=0; j=set[j]); j++); if (--j>=0) nhit[j]++; } } for(j=1; j=set[j]); j++); if (--j>=0) { if (j>0) pos=nhit[j-1]+nh2[j]; else pos=*nh2; places[i]=pos; cables[pos]=number+1; times[pos]=tim_raw[number]; charges[pos]=chg_raw[number]; avcharges[j]+=chg_raw[number]; positions[3*pos]=pmt_loc[3*number]; positions[3*pos+1]=pmt_loc[3*number+1]; positions[3*pos+2]=pmt_loc[3*number+2]; nh2[j]++; } else places[i]=-1; } } for(j=0; j0) avcharges[j]/=nh2[j]; delete(nh2); } hits::hits(int ns,float *set,pmt_geometry *geom,int nh,int *cab,float *t,float *q) { int i,j,number,pos,*nh2; // count #hits per charge bin and create offset 'pointers' nchargeset=ns; nhit=new int[nchargeset]; avcharges=new float[nchargeset]; for(j=0; j=set[j]); j++); if (--j>=0) nhit[j]++; } for(j=1; j=set[j]); j++); if (--j>=0) { if (j>0) pos=nhit[j-1]+nh2[j]; else pos=*nh2; places[i]=pos; number=(int)(cab[i]); cables[pos]=number--; times[pos]=t[i]; charges[pos]=q[i]; avcharges[j]+=charges[pos]; geom->position(positions+3*pos,number); nh2[j]++; } else places[i]=-1; } for(j=0; j0) avcharges[j]/=nh2[j]; delete(nh2); } hits::hits(int ns,float *set,pmt_geometry *geom,comtype2 *itevent) { int i,j,number,pos,*nh2; // count #hits per charge bin and create offset 'pointers' nchargeset=ns; nhit=new int[nchargeset]; avcharges=new float[nchargeset]; for(j=0; jit_index; for(i=0; ihits[i][1]>=set[j]); j++); if (--j>=0) nhit[j]++; } for(j=1; jhits[i][1]>=set[j]); j++); if (--j>=0) { if (j>0) pos=nhit[j-1]+nh2[j]; else pos=*nh2; places[i]=pos; number=(int)fabs(itevent->hits[i][0]); cables[pos]=number--; times[pos]=itevent->hits[i][2]; charges[pos]=itevent->hits[i][1]; avcharges[j]+=charges[pos]; geom->position(positions+3*pos,number); nh2[j]++; } else places[i]=-1; } for(j=0; j0) avcharges[j]/=nh2[j]; delete(nh2); } // frees hit storage arrays hits::~hits(void) { delete(nhit); delete(places); delete(positions); delete(times); delete(charges); delete(avcharges); delete(cables); } // ********************************************** // average hit times between tmin and tmax // ********************************************** float hits::time_av(float tmin,float tmax) { int i,n; float t_av; if (nhit[nchargeset-1]<=0) return(-1); for(t_av=0,n=0,i=0; itmin) && (times[i]0) { if ((j==3) && (sprod+root<0)) return(0); vert[j]=-(root+sprod)/sq1; vert[j+4]=(root-sprod)/sq1; } else { if ((j==3) && (sprod>root)) return(0); vert[j]=(root-sprod)/sq1; vert[j+4]=-(root+sprod)/sq1; } if ((j==3) && (vert[7]>0)) for(i=0; i<3; i++) vert[i]=-(vert[3]*vert[i]+vert[i+4])+positions[3*imin+i]; else for(i=0; i<3; i++) { ii=colind[i]; dist=vert[ii]; vert[ii]=-(vert[j]*vert[ii]+vert[ii+4]); vert[ii+4]=-(vert[j+4]*dist+vert[ii+4]); if (ii<3) { vert[ii]+=positions[3*imin+ii]; vert[ii+4]+=positions[3*imin+ii]; } else { vert[3]=CM_TO_NS*vert[3]+times[imin]; vert[7]=CM_TO_NS*vert[7]+times[imin]; } } if (j==3) { vert[3]=CM_TO_NS*vert[3]+times[imin]; vert[7]=CM_TO_NS*vert[7]+times[imin]; } else { vert[j]+=positions[3*imin+j]; vert[j+4]+=positions[3*imin+j]; } if (vert[7]>times[imin]) return(1); else return(2); } // ********************************************** // find vertex from four-hit combination // ********************************************** int hits::vertex4(int *fourcombo,int *hitlist,int listsize, int &hbegin,int &hend,float tsig, float cyl_radius,float cyl_height,float dwallmin, double *vert,float *gdn) { float rad,height,tmin,tmax; int h,thishit,s,ns=vertex4(fourcombo,vert); for(s=0; s0) && (times[hitlist[hbegin]]>=tmin)) hbegin--; while((hbegin0) && (times[hitlist[hend]]>tmax)) hend--; //printf("begin %6.1f end %6.1f ",times[hitlist[hbegin]],times[hitlist[hend]]); for(h=hbegin; h<=hend; h++) { thishit=hitlist[h]; if ((thishit!=fourcombo[0]) && (thishit!=fourcombo[1]) && (thishit!=fourcombo[2]) && (thishit!=fourcombo[3])) { t=times[thishit]-tof(vertex,thishit)-vert[4*s+3]; t/=tsig; t*=-0.5*t; if (t>=-4.5) gdn[s]+=exp(t); } } //printf("%8.4f\n",gdn[s]); } return(ns); }