12 #include "TRotation.h"
13 #include "Math/Rotation3D.h"
14 #include "Math/Vector3Dfwd.h"
23 extern const double _rH1[3] = {-2.161414928e6, -3.834695183e6, 4.600350224e6};
24 extern const double _eH1x[3] = {-0.223891216, 0.799830697, 0.556905359};
25 extern const double _eH1y[3] = {-0.913978490, 0.026095321, -0.404922650};
28 extern const double _rH2[3] = {-2.161414928e6, -3.834695183e6, 4.600350224e6};
29 extern const double _eH2x[3] = {-0.223891216, 0.799830697, 0.556905359};
30 extern const double _eH2y[3] = {-0.913978490, 0.026095321, -0.404922650};
33 extern const double _rL1[3] = {-7.427604192e4, -5.496283721e6, 3.224257016e6};
34 extern const double _eL1x[3] = {-0.954574615, -0.141579994, -0.262187738};
35 extern const double _eL1y[3] = { 0.297740169, -0.487910627, -0.820544948};
43 extern const double _rG1[3] = {3.8563112e6, 6.665978e5, 5.0196406e6};
44 extern const double _eG1x[3] = {-0.445184239, 0.866534205, 0.225675575};
45 extern const double _eG1y[3] = {-0.626000687, -0.552167273, 0.550667271};
48 extern const double _rV1[3] = {4.5463741e6, 8.429897e5, 4.378577e6};
49 extern const double _eV1x[3] = {-0.700458309, 0.208487795, 0.682562083};
50 extern const double _eV1y[3] = {-0.053791331, -0.969082169, 0.240803326};
53 extern const double _rT1[3] = {-3.946409e6, 3.366259e6, 3.6991507e6};
54 extern const double _eT1x[3] = {0.648969405, 0.760814505, 0};
55 extern const double _eT1y[3] = {-0.443713769, 0.378484715, -0.812322234};
58 extern const double _rJ1[3] = {-3.7685777e6, 3.49218552e6, 3.76722999e6};
59 extern const double _eJ1x[3] = {-0.5014606643, -0.824380550, 0.262552681483};
60 extern const double _eJ1y[3] = {0.6313793503,-0.1412130545258,0.762508353526};
63 extern const double _rA1[3] = {-2.3567784e6, 4.8970238e6, -3.3173147e6};
64 extern const double _eA1x[3] = {-9.01077021322091554e-01, -4.33659084587544319e-01, 0};
65 extern const double _eA1y[3] = {-2.25940560005277846e-01, 4.69469807139026196e-01, 8.53550797275327455e-01};
68 extern const double _rO1[3] = {4.392467e6, 0.9295086e6, 4.515029e6};
69 extern const double _eO1x[3] = {-0.644504130, 0.573655377, 0.50550364};
70 extern const double _eO1y[3] = {0., 0., 0.};
73 extern const double _rN1[3] = {4.64410999868e6, 1.04425342477e6, 4.23104713307e6};
74 extern const double _eN1x[3] = {-0.62792641437, 0.56480832712, 0.53544371484};
75 extern const double _eN1y[3] = {0., 0., 0.};
78 extern const double _rE1[3] = {4.37645395452e6, 4.75435044067e5, 4.59985274450e6};
79 extern const double _eE1x[3] = {-0.62792641437, 0.56480832712, 0.53544371484};
80 extern const double _eE1y[3] = {0., 0., 0.};
92 for(
int i=0;
i<3;
i++){
105 this->TFmap.rate(4096);
128 if(strstr(name,
"Virgo")) {pRv=
_rV1; pEx=
_eV1x; pEy=
_eV1y;xifo=
true;}
129 if(strstr(name,
"VIRGO")) {pRv=
_rV1; pEx=
_eV1x; pEy=
_eV1y;xifo=
true;}
130 if(strstr(name,
"GEO")) {pRv=
_rG1; pEx=
_eG1x; pEy=
_eG1y;xifo=
true;}
131 if(strstr(name,
"TAMA")) {pRv=
_rT1; pEx=
_eT1x; pEy=
_eT1y;xifo=
true;}
132 if(strstr(name,
"LCGT")) {pRv=
_rJ1; pEx=
_eJ1x; pEy=
_eJ1y;xifo=
true;}
139 cout <<
"detector::detector - Error : detector " << name
140 <<
" is not in present in the builtin list" << endl;
146 for(
int i=0;
i<3;
i++){
147 this->Rv[
i] = pRv[
i];
148 this->Ex[
i] = pEx[
i];
149 this->Ey[
i] = pEy[
i];
151 if(strstr(name,
"A2")) this->rotate(45);
160 this->TFmap.rate(4096);
182 for(
int i=0;
i<3;
i++){
183 this->Rv[
i] = uRv[
i];
184 this->Ex[
i] = uEx[
i];
185 this->Ey[
i] = uEy[
i];
195 this->TFmap.rate(4096);
206 if(strlen(this->dP.name)>0)
return this->dP;
214 XYZVector iRv(this->Rv[0],this->Rv[1],this->Rv[2]);
215 TVector3 vRv(iRv.X(),iRv.Y(),iRv.Z());
218 XYZVector iEx(this->Ex[0],this->Ex[1],this->Ex[2]);
219 XYZVector iEy(this->Ey[0],this->Ey[1],this->Ey[2]);
220 XYZVector iEZ(0.,0.,1.);
222 TVector3 vEx(iEx.X(),iEx.Y(),iEx.Z());
223 TVector3 vEy(iEy.X(),iEy.Y(),iEy.Z());
224 TVector3 vEZ(iEZ.X(),iEZ.Y(),iEZ.Z());
231 TVector3 vEe=vEZ.Cross(vRv);
234 double cosExAngle=vEe.Dot(vEx);
235 if(cosExAngle>1.) cosExAngle=1.;
236 if(cosExAngle<-1.) cosExAngle=-1.;
237 double cosEyAngle=vEe.Dot(vEy);
238 if(cosEyAngle>1.) cosEyAngle=1.;
239 if(cosEyAngle<-1.) cosEyAngle=-1.;
241 double ExAngle=acos(cosExAngle)*
rad2deg;
242 double EyAngle=acos(cosEyAngle)*
rad2deg;
247 if(vEn.Dot(vRv)<0) ExAngle*=-1;
249 if(vEn.Dot(vRv)<0) EyAngle*=-1;
252 ExAngle=fmod(90-ExAngle,360.);
254 EyAngle=fmod(90-EyAngle,360.);
256 double latitude,longitude,elevation;
261 sprintf(idP.name,
"%s",this->Name);
262 idP.latitude = latitude*
rad2deg;
263 idP.longitude = longitude*
rad2deg;
264 idP.elevation = elevation;
278 double si = sin(a*
PI/180.);
279 double co = cos(a*
PI/180.);
281 for(
int i=0;
i<3;
i++) {
295 for(
int i=0;
i<3;
i++) axy+=ax[
i]*ay[
i];
299 for(
int i=0;
i<3;
i++) {aw[
i]=-ax[
i]*axy+ay[
i]; aww+=aw[
i]*aw[
i];}
300 for(
int i=0;
i<3;
i++) aw[
i]/=sqrt(aww);
301 for(
int i=0;
i<3;
i++) this->Ex[
i]=ax[
i]*co+aw[
i]*si;
305 for(
int i=0;
i<3;
i++) {aw[
i]=-ax[
i]+ay[
i]*axy; aww+=aw[
i]*aw[
i];}
306 for(
int i=0;
i<3;
i++) aw[
i]/=sqrt(aww);
307 for(
int i=0;
i<3;
i++) this->Ey[
i]=ay[
i]*co+aw[
i]*si;
312 if(strlen(this->dP.name)>0) {
313 this->dP.AzX = fmod(this->dP.AzX-a,360.);
314 this->dP.AzY = fmod(this->dP.AzY-a,360.);
331 for (
int i=0;
i<
n;
i++) {
339 for (
int i=0;
i<
n;
i++) {
355 for(
int i=0;
i<3;
i++){
356 this->Rv[
i] = value.
Rv[
i];
357 this->Ex[
i] = value.
Ex[
i];
358 this->Ey[
i] = value.
Ey[
i];
385 double tsRate =TFmap.wavearray<double>::rate();
387 waveForm.resize(
size_t(tsRate));
388 waveForm.rate(tsRate);
389 waveForm.setWavelet(*(TFmap.pWavelet));
406 for (
int i=0;
i<IWFP.size();
i++) {
411 for(
int i=0;
i<value.
IWFP.size();
i++) {
413 *wf = *value.
IWFP[
i];
433 for (
int i=0;
i<value.
IWFP.size();
i++) {
438 for(
int i=0;
i<IWFP.size();
i++) {
441 value.
IWFP.push_back(wf);
453 DT[0] = Ex[0]*Ex[0]-Ey[0]*Ey[0];
454 DT[1] = Ex[0]*Ex[1]-Ey[0]*Ey[1];
455 DT[2] = Ex[0]*Ex[2]-Ey[0]*Ey[2];
458 DT[4] = Ex[1]*Ex[1]-Ey[1]*Ey[1];
459 DT[5] = Ex[1]*Ex[2]-Ey[1]*Ey[2];
463 DT[8] = Ex[2]*Ex[2]-Ey[2]*Ey[2];
465 if (strcmp(Name,
"O1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
466 if (strcmp(Name,
"N1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
467 if (strcmp(Name,
"E1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
477 theta *=
PI/180.; phi *=
PI/180.; psi *=
PI/180.;
479 double cT = cos(theta);
480 double sT = sin(theta);
481 double cP = cos(phi);
482 double sP = sin(phi);
502 fp = (cT*cP*d11 + cT*sP*d21 - sT*d31)*cT*cP
503 + (cT*cP*d12 + cT*sP*d22 - sT*d32)*cT*sP
504 - (cT*cP*d13 + cT*sP*d23 - sT*d33)*sT
506 - (cP*d22-sP*d12)*cP;
508 fx = -(cT*cP*d11 + cT*sP*d21 - sT*d31)*sP
509 + (cT*cP*d12 + cT*sP*d22 - sT*d32)*cP
510 + (cP*d21-sP*d11)*cT*cP
511 + (cP*d22-sP*d12)*cT*sP
512 - (cP*d23-sP*d13)*sT;
517 a = fp*cos(2*psi)+fx*sin(2*psi);
518 b = -fp*sin(2*psi)+fx*cos(2*psi);
525 fp = -(cT*cP*d11 + cT*sP*d21 - sT*d31)*cT*cP
526 - (cT*cP*d12 + cT*sP*d22 - sT*d32)*cT*sP
527 + (cT*cP*d13 + cT*sP*d23 - sT*d33)*sT
529 - (cP*d22-sP*d12)*cP;
547 double a,b,rms,fl,fh;
548 double R = this->TFmap.
rate();
549 int L =
int(this->TFmap.maxLayer()+1);
551 waveForm.setWavelet(*(TFmap.pWavelet));
553 rms = wc.
getwave(ID,waveForm,atype,index);
555 if(rms==0.)
return rms;
559 waveBand = waveForm; waveBand = 0.;
560 fh = waveBand.gethigh();
561 fl = waveBand.getlow();
562 l =
int(this->TFmap.getLevel()) -
int(waveBand.getLevel());
566 if(l<0) { waveBand.Inverse(-l); }
567 else { waveBand.Forward(l); }
569 a = waveBand.start()*R;
570 i =
int(a + ((a>0) ? 0.1 : -0.1));
571 k = waveBand.size(); j = 0;
572 if(i<0) { k +=
i; j = -
i; i = 0; }
573 if((i/L)*L != i) cout<<
"detector::getwave() time mismatch: "<<L<<
" "<<i<<
"\n";
575 waveBand.cpf(this->TFmap,k,i,j);
578 n =
int(2*L*fl/R+0.1)-1;
579 m =
int(2*L*fh/R+0.1)+1;
582 if(m<=n) { n=0; m=
L; }
590 for(k=n; k<
m; k++) { w.
getLayer(x,k); waveBand.putLayer(x,k); }
595 waveForm *= atype==
'w' ? 1./sqrt(this->
rate/R) : 1.;
600 double sTARt= waveForm.start();
601 size_t I = waveForm.size();
603 double sum = waveForm.data[
M]*waveForm.data[
M];
608 for(i=1; i<
int(M); i++) {
609 a = waveForm.data[M-
i];
610 b = waveForm.data[M+
i];
612 if(sum/hrss > 0.999 && i/R>0.05)
break;
615 if(n <
int(M-1)) i = size_t(n);
620 waveForm.cpf(w,2*i,M-i);
621 waveForm.resize(2*i);
622 waveForm.start(sTARt+(M-i)/R);
625 waveBand.
cpf(w,2*i,M-i);
626 waveBand.resize(2*i);
627 waveBand.start(sTARt+(M-i)/R);
641 size_t n = SM.
size();
647 z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
664 size_t n = SM.
size();
670 z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
683 double x = theta*
PI/180.;
684 double y = phi*
PI/180.;
685 double z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
695 skymap Sp(sms,t1,t2,p1,p2);
696 skymap Sx(sms,t1,t2,p1,p2);
697 size_t n = Sp.
size();
722 size_t n = Sp.
size();
745 cout<<
"wseries::setFilter(): not applicable to WDM TFmaps\n";
748 size_t i,
j,
k,
n,ii,jj;
749 size_t M = TFmap.maxLayer()+1;
750 size_t L = TFmap.getLevel();
751 size_t m = M*TFmap.pWavelet->m_H;
752 size_t N = M*(1<<upL);
755 std::vector<delayFilter>
F;
757 for(i=0; i<
K; i++) v.
value[i] = 0.;
758 for(i=0; i<M; i++) F.push_back(v);
781 cout<<
"wsize: "<<w.
size()<<endl;
786 double* pb = (
double * )malloc(m*
sizeof(
double));
787 double** pp = (
double **)malloc(m*
sizeof(
double*));
788 for(i=0; i<
m; i++) pp[i] = w.
data + i + (
int(j) -
int(m/2));
819 for(k=2; k<W.
size(); k+=2) {
849 for(k=m-1; k>=0; k--) {
851 if(abs(offst) >32767)
continue;
852 inDex = short(offst);
853 vaLue = float(pb[pp[k]-p0]);
854 if(
fabs(vaLue)<1.
e-4)
break;
856 offst -= (offst/
M)*M;
857 if(offst<0) offst +=
M;
862 offst -= (offst/
M)*M;
863 if(offst<0) offst +=
M;
864 if(offst !=
int(s)) cout<<
"setFilter error 2: "<<offst<<
" "<<s<<endl;
868 for(
size_t kk=0; kk<
K; kk++)
871 if(jj>=K) {cout<<jj<<endl;
continue;}
873 F[ii].value[jj] = vaLue;
874 F[ii].index[jj] = -inDex;
896 sum += F[
i].value[
k]*F[
i].value[
k];
897 if(n && F[i].
value[k] == 0.)
printf(
"%4d %4d %d4 \n",
int(n),
int(i),
int(k));
898 if(v.
value[k] != F[i].value[k] ||
899 v.
index[k] != F[i].index[k]) cout<<
"setFilter error 4\n";
902 if(sum<0.97)
printf(
"%4d %4d %8.5f \n",
int(n),
int(i),sum);
919 std::vector<delayFilter>().swap(
filter);
922 for(
size_t k=0;
k<
K;
k++) {
935 if ( (fp=fopen(fname,
"wb")) == NULL ) {
936 cout <<
" DumpBinary() error : cannot open file " << fname <<
". \n";
940 size_t M = size_t(TFmap.maxLayer()+1);
942 size_t N = this->nDFS;
943 size_t n = K *
sizeof(float);
944 size_t m = K *
sizeof(short);
949 fwrite(&K,
sizeof(
size_t), 1, fp);
950 fwrite(&M,
sizeof(
size_t), 1, fp);
951 fwrite(&N,
sizeof(
size_t), 1, fp);
959 fwrite(value.data, n, 1, fp);
960 fwrite(index.
data, m, 1, fp);
974 if ( (fp=fopen(fname,
"rb")) == NULL ) {
975 cout <<
" DumpBinary() error : cannot open file " << fname <<
". \n";
983 fread(&K,
sizeof(
size_t), 1, fp);
984 fread(&M,
sizeof(
size_t), 1, fp);
985 fread(&N,
sizeof(
size_t), 1, fp);
987 size_t n = K *
sizeof(float);
988 size_t m = K *
sizeof(short);
995 this->clearFilter();
filter.reserve(N*M);
1000 v.
value.push_back(0.);
1001 v.
index.push_back(0);
1004 for(i=0; i<
N; i++) {
1005 for(j=0; j<
M; j++) {
1006 fread(value.data, n, 1, fp);
1007 fread(index.
data, m, 1, fp);
1008 for(k=0; k<
K; k++) {
1039 cout<<
"wseries::delay(): not applicable to WDM TFmaps\n";
1049 int n =
int(t*TFmap.wavearray<double>::rate());
1050 int m = n>0 ? (n+N/2-1)/N : (n-N/2)/
N;
1051 int l = TFmap.pWavelet->m_H/4+2;
1057 cout<<
"delay="<<t<<
" m="<<m<<
" n="<<n<<
" M="<<M<<
" K="<<K<<
" N="<<N<<endl;
1060 double*
A = (
double*)malloc(K*
sizeof(
double));
1061 int* I = (
int*)malloc(K*
sizeof(
int));
1063 for(i=0; i<
M; i++) {
1068 je = m<0 ? w.
size() +(m-
l)*M : w.
size() -l*
M;
1073 for(j=jb; j<je; j+=
M) {
1075 TFmap.data[
j] = A[0]*p[I[0]] + A[1]*p[I[1]] + A[2]*p[I[2]] + A[3]*p[I[3]]
1076 + A[4]*p[I[4]] + A[5]*p[I[5]] + A[6]*p[I[6]] + A[7]*p[I[7]]
1077 + A[8]*p[I[8]] + A[9]*p[I[9]] + A[10]*p[I[10]] + A[11]*p[I[11]]
1078 + A[12]*p[I[12]] + A[13]*p[I[13]] + A[14]*p[I[14]] + A[15]*p[I[15]];
1083 for(j=jb; j<je; j+=
M) {
1085 TFmap.data[
j] = A[0]*p[I[0]] + A[1]*p[I[1]] + A[2]*p[I[2]] + A[3]*p[I[3]]
1086 + A[4]*p[I[4]] + A[5]*p[I[5]] + A[6]*p[I[6]] + A[7]*p[I[7]]
1087 + A[8]*p[I[8]] + A[9]*p[I[9]] + A[10]*p[I[10]] + A[11]*p[I[11]]
1088 + A[12]*p[I[12]] + A[13]*p[I[13]] + A[14]*p[I[14]] + A[15]*p[I[15]]
1089 + A[16]*p[I[16]] + A[17]*p[I[17]] + A[18]*p[I[18]] + A[19]*p[I[19]]
1090 + A[20]*p[I[20]] + A[21]*p[I[21]] + A[22]*p[I[22]] + A[23]*p[I[23]]
1091 + A[24]*p[I[24]] + A[25]*p[I[25]] + A[26]*p[I[26]] + A[27]*p[I[27]]
1092 + A[28]*p[I[28]] + A[29]*p[I[29]] + A[30]*p[I[30]] + A[31]*p[I[31]];
1097 for(j=jb; j<je; j+=
M) {
1101 while(k-- > 0) { *q += *(A++) * p[*(I++)]; }
1114 if(!nRMS.size())
return 0.;
1115 if(
int(I)>TFmap.maxLayer())
return 0.;
1119 size_t N = nRMS.maxLayer()+1;
1120 size_t M = TFmap.maxLayer()+1;
1122 size_t m = (I+1)*(N/M);
1123 double t = TFmap.
start();
1124 double T = nRMS.start();
1125 slice S = TFmap.getSlice(I);
1126 double rATe = TFmap.wrate();
1128 double rESn = nRMS.frequency(1)-nRMS.frequency(0);
1133 cout<<
"detector::getNoise(): invalid noise rms array nRMS\n";
1137 size_t L = size_t(TFmap.getlow()/rESn+0.1);
1138 size_t H = size_t(TFmap.gethigh()/rESn+0.1);
1139 if(n>=H || m<=L)
return 0.;
1144 cout<<
"detector::getNoise():: invalid noise rms array nRMS\n";
1148 S = nRMS.getSlice(n);
1149 size_t K = S.
size();
1157 for(l=n; l<
m; l++) {
1158 S = nRMS.getSlice(l);
1159 g = (n<L || m>H) ? 1.e10 : 1.;
1160 for(j=0; j<rms.
size(); j++) {
1166 for(j=0; j<rms.
size(); j++) {
1167 rms.
data[
j] = sqrt((
double(m)-
double(n))/rms.
data[j]);
1179 int inRMS =
int((t-nRMS.start())*nRMS.rate());
1180 int inVAR = nVAR.size() ?
int((t-nVAR.start())*nVAR.rate()) : 0;
1182 if(inRMS < 0) inRMS = 0;
1183 else if(
size_t(inRMS)>=K &&
K) inRMS = K-1;
1185 if(inVAR <= 0) inVAR = 0;
1186 else if(
size_t(inVAR)>=nVAR.size()) inVAR = nVAR.size()-1;
1188 for(l=n; l<
m; l++) {
1189 S = nRMS.getSlice(l);
1190 g = (n<L || m>H) ? 1.e10 : 1.;
1195 RMS /= double(m)-double(n);
1198 if(!nVAR.size() || rESn*n<nVAR.getlow() || rESn*m>nVAR.gethigh())
return RMS;
1200 return RMS*double(nVAR.data[inVAR]);
1211 size_t M = wc->
size();
1213 if(!M)
return false;
1215 if(!nRMS.size())
return false;
1217 size_t max_layer = nRMS.maxLayer();
1222 int K = nRMS.size()/(max_layer+1);
1223 double To = nRMS.start();
1224 double Ro = nRMS.wrate();
1225 double Fo = nRMS.frequency(0);
1226 double dF = nRMS.frequency(1)-Fo;
1227 double fl = wc->
getlow()-0.1;
1228 double fh = wc->
gethigh()+0.1;
1231 Fo = Fo==0. ? 0.5 : 0.;
1237 int(p->
rate/Ro+0.01) < 1 ||
1239 cout<<
"detector::setrms() - illegal pixel from zero level\n";
1245 n = size_t(f/dF+0.6);
1246 F = (x+1)*p->
rate/2.;
1247 m =
size_t(F/dF+0.6);
1248 if(m>max_layer) m=max_layer+1;
1253 if(k>=K) k -= k ? 1 : 0;
1254 if(k<0 || n>=m || k>=K) {
1255 cout<<
"detector::setrms() - invalid input: ";
1256 cout<<k<<
" "<<n<<
" "<<m<<
" "<<f<<
" "<<F<<
" "<<t<<endl;
1257 cout<<p->
frequency<<
" "<<p->
rate/2.<<
" "<<dF<<
" "<<fl<<endl;
1264 for(j=n; j<
m; j++) {
1265 S = nRMS.getSlice(j);
1266 g = (f<fl || F>fh) ? 1.e10 : 1.;
1271 if(nVAR.size() && f<nVAR.gethigh() && F>nVAR.getlow()) r *= nVAR.get(t);
1284 double dF = TFmap.frequency(1)-TFmap.frequency(0);
1285 double fl =
fabs(f1)>0. ?
fabs(f1) : this->TFmap.getlow();
1286 double fh =
fabs(f2)>0. ?
fabs(f2) : this->TFmap.gethigh();
1287 size_t n = TFmap.pWavelet->m_WaveType==
WDMT ? size_t((fl+dF/2.)/dF+0.1) : size_t(fl/dF+0.1);
1288 size_t m = TFmap.pWavelet->m_WaveType==
WDMT ? size_t((fh+dF/2.)/dF+0.1)-1 : size_t(fh/dF+0.1)-1;
1289 size_t M = this->TFmap.maxLayer()+1;
1294 for(i=0; i<
int(M); i++) {
1296 if((f1>=0 && i>=n) && (f2>=0 && i<=
m))
continue;
1297 if((f1<0 && i<n) || (f2<0 && i>
m))
continue;
1298 if((f1<0 && f2>=0 && i<n))
continue;
1299 if((f1>=0 && f2<0 && i>=m))
continue;
1301 this->TFmap.getLayer(w,i+0.01); w=0.; this->TFmap.putLayer(w,i+0.01);
1302 this->TFmap.getLayer(w,-i-0.01); w=0.; this->TFmap.putLayer(w,-i-0.01);
1350 int j,nstop,nstrt,
n,
m,J;
1352 double a,b,
T,E,
D,H,
f;
1353 double R = this->
rate;
1354 size_t K = pT->size();
1357 bool pOUT = dT>0. ?
false :
true;
1361 if(wi.
size() != TFmap.size()) {
1362 cout<<
"setsim(): mismatch between MDC size "
1363 <<wi.
size()<<
" and data size "<<TFmap.size()<<
"\n";
1368 if(!this->nRMS.size() || !this->TFmap.size())
return 0;
1370 if(this->HRSS.size() !=
K) {
1371 this->HRSS.resize(K);
1372 this->ISNR.resize(K);
1373 this->FREQ.resize(K);
1374 this->BAND.resize(K);
1375 this->TIME.resize(K);
1376 this->TDUR.resize(K);
1394 if(!W.
white(this->nRMS,1)) {
1395 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1397 if(!W.
white(this->nRMS,-1)) {
1398 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1407 if(!W.
white(this->nRMS)) {
1408 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1417 size_t N = w.
size();
1418 double rate = w.wavearray<double>::rate();
1419 double bgps = w.
start()+offset+1.;
1420 double sgps = w.
start()-offset+N/rate-1.;
1422 for(k=0; k<
K; k++) {
1424 T = (*pT)[
k] - w.
start();
1426 nstrt =
int((T - dT)*rate);
1427 nstop =
int((T + dT)*rate);
1428 if(nstrt<=0) nstrt = 0;
1429 if(nstop>=
int(N)) nstop =
N;
1430 if(nstop<=0)
continue;
1431 if(nstrt>=
int(N))
continue;
1434 for(j=nstrt; j<nstop; j++) {
1443 nstrt =
int((T - dT)*rate);
1444 nstop =
int((T + dT)*rate);
1445 if(nstrt<=0) nstrt = 0;
1446 if(nstop>=
int(N)) nstop =
N;
1447 if(nstop<=0)
continue;
1448 if(nstrt>=
int(N))
continue;
1451 for(j=nstrt; j<nstop; j++) {
1462 if(n>=
int(S.size())) n = S.
size()-1;
1464 for(j=nstrt; j<nstop; j++) {
1465 a = w.
data[
j]*(j/rate-
T);
1474 if(T<bgps || T>sgps)
continue;
1476 this->ISNR.data[
k] = E;
1477 this->TIME.data[
k] =
T;
1478 this->TDUR.data[
k] =
D;
1479 this->HRSS.data[
k] = H;
1482 for(i=0; i<I; i++) {
1487 this->FREQ.data[
k] += f*a*
a;
1488 this->BAND.
data[
k] += f*f*a*
a;
1492 this->FREQ.
data[
k] /= E;
1493 a = this->FREQ.data[
k];
1494 b = this->BAND.data[
k]/E - a*
a;
1495 this->BAND.
data[
k] = b>0. ? sqrt(b) : 0.;
1499 double time = this->TIME.data[
k];
1500 double tdur = this->TDUR.data[
k];
1501 double tDur = 6*tdur;
1502 if (tDur > dT) tDur =
dT;
1503 int nStrt =
int((time-tDur-w.
start())*rate);
1504 int nStop =
int((time+tDur-w.
start())*rate);
1505 if(nStrt<0) nStrt = 0;
1506 if(nStop>
int(N)) nStop =
N;
1510 size_t I =
int(2.*dT*rate);
1514 for (j=0;j<I;j++) ms += ((OS+j>=0)&&(OS+j<N)) ? w.
data[OS+j]*w.
data[OS+j] : 0.;
1518 double sum = ((OS>=0)&&(OS<
N)) ? w.
data[
OS]*w.
data[
OS] : 0.;
1519 for(j=1; j<
int(I/2); j++) {
1520 a = ((OS-j>=0)&&(OS-j<N)) ? w.
data[OS-j] : 0.;
1521 b = ((OS+j>=0)&&(OS+j<
N)) ? w.
data[OS+
j] : 0.;
1523 if(sum/ms > 0.999)
break;
1528 if(nStrt<0) nStrt = 0;
1529 if(nStop>
int(N)) nStop = N;
1535 for(j=nStrt; j<nStop; j++) {
1541 double f_low = this->TFmap.getlow();
1542 double f_high = this->TFmap.gethigh();
1545 double Fs=((double)wf->
rate()/(double)wf->
size())/2.;
1546 for (j=0;j<wf->
size()/2;j+=2) {
1548 if((f<f_low)||(f>f_high)) {wf->
data[
j]=0.;wf->
data[j+1]=0.;}
1555 for(j=0;j<wf->
size();j++) {
1562 for(j=0;j<wf->
size();j++) {
1563 a = wf->
data[
j]*(j/rate-
T);
1570 this->ISNR.data[
k] = E;
1571 this->TIME.data[
k] =
T;
1572 this->TDUR.data[
k] =
D;
1587 for(j=nstrt; j<nstop; j++) {
1595 I =
int(2.*dT*rate);
1596 OS =
int((T - dT)*rate);
1599 for (j=0;j<I;j++) ms += ((OS+j>=0)&&(OS+j<N)) ? hot.
data[OS+j]*hot.
data[OS+j] : 0.;
1603 for(j=1; j<
int(I/2); j++) {
1604 a = ((OS-j>=0)&&(OS-j<N)) ? hot.
data[OS-j] : 0.;
1605 b = ((OS+j>=0)&&(OS+j<
N)) ? hot.
data[OS+
j] : 0.;
1607 if(sum/ms > 0.999)
break;
1610 nStrt =
int(T*rate)-
j;
1611 nStop =
int(T*rate)+
j;
1612 if(nStrt<0) nStrt = 0;
1613 if(nStop>
int(N)) nStop =
N;
1620 for(j=nStrt; j<nStop; j++) WF->
data[j-nStrt] = hot.
data[j];
1623 IWFID.push_back(-(k+1));
1628 printf(
"%3d T+-dT: %8.3f +-%5.3f, F+-dF: %4.0f +-%4.0f, SNR: %6.1e, hrss: %6.1e\n",
1629 int(M), T-bgps, D, FREQ.data[k], BAND.data[k], E, H);
1646 size_t K = pT->size();
1647 size_t N = wi.
size();
1656 for(k=0; k<
K; k++) {
1659 T = (*pT)[
k] - w.
start();
1661 nstrt =
int((T - dT)*rate);
1662 nstop =
int((T + dT)*rate);
1663 if(nstrt<=0) nstrt = 0;
1664 if(nstop>=
int(N)) nstop =
N;
1665 if(nstop<=0)
continue;
1666 if(nstrt>=
int(N))
continue;
1668 for(j=nstrt; j<nstop; j++) {
1684 if(!TFmap.size())
return;
1685 double R = this->TFmap.wavearray<double>::rate();
1686 double T = this->getTau(theta,phi);
1687 size_t n = size_t(
fabs(T)*R);
1688 size_t m = this->TFmap.size();
1693 if(T<0) TFmap.
cpf(w,m-n,0,n);
1694 else TFmap.cpf(w,m-n,n,0);
1703 if(!x.
size())
return;
1704 double R = x.
rate();
1705 double T = this->getTau(theta,phi);
1706 size_t n = size_t(
fabs(T)*R);
1707 size_t m = x.
size();
1712 if(T<0) x.
cpf(w,m-n,0,n);
1713 else x.
cpf(w,m-n,n,0);
1722 if(!x.
size())
return;
1723 double R = x.
rate();
1724 size_t n = size_t(
fabs(T)*R+0.5);
1725 size_t m = x.
size();
1730 if(T<0) x.
cpf(w,m-n,0,n);
1731 else x.
cpf(w,m-n,n,0);
1742 for(
size_t i=0;
i<wf->
size();
i++) {
1747 return E>0. ? wf->
start()+T/E/wf->
rate() : 0.;
1758 for(
size_t i=0;
i<wf->
size();
i+=2) {
1764 return E>0. ? 0.5*F*wf->
rate()/E/wf->
size() : 0.;
1774 if(theta_t>0) LAT=
'N';
else {LAT=
'S';theta_t=-theta_t;}
1775 int theta_d =
int(theta_t);
1776 int theta_m =
int((theta_t-theta_d)*60);
1777 float theta_s = (theta_t-theta_d-theta_m/60.)*3600.;
1781 if(phi_t>0) LON=
'E';
else {LON=
'W';phi_t=-phi_t;}
1782 int phi_d =
int(phi_t);
1783 int phi_m =
int((phi_t-phi_d)*60);
1784 float phi_s = (phi_t-phi_d-phi_m/60.)*3600.;
1787 cout <<
"----------------------------------------------" << endl;
1788 cout <<
"IFO : " << xdP.
name <<
" (Geographic Coordinates) " << endl;
1789 cout <<
"----------------------------------------------" << endl;
1791 cout <<
"latitude : " << xdP.
latitude <<
" longitude : " << xdP.
longitude << endl;
1793 cout <<
"LAT : " << LAT <<
" " << theta_d <<
", " << theta_m <<
", " << theta_s << endl;
1794 cout <<
"LON : " << LON <<
" " << phi_d <<
", " << phi_m <<
", " << phi_s << endl;
1798 cout <<
"Rv : " << Rv[0] <<
" " << Rv[1] <<
" " << Rv[2] <<
" " << endl;
1800 cout <<
"Ex : " << Ex[0] <<
" " << Ex[1] <<
" " << Ex[2] <<
" " << endl;
1802 cout <<
"Ey : " << Ey[0] <<
" " << Ey[1] <<
" " << Ey[2] <<
" " << endl;
1804 cout.precision(precision);
1805 cout <<
"Ex-North Angle Clockwise (deg) : " << xdP.
AzX << endl;
1806 cout <<
"Ey-North Angle Clockwise (deg) : " << xdP.
AzY << endl;
1816 void detector::Streamer(TBuffer &R__b)
1821 if (R__b.IsReading()) {
1822 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
1823 TNamed::Streamer(R__b);
1824 R__b.ReadStaticArray((
char*)Name);
1826 R__b.ReadStaticArray((
double*)Rv);
1827 R__b.ReadStaticArray((
double*)Ex);
1828 R__b.ReadStaticArray((
double*)Ey);
1829 R__b.ReadStaticArray((
double*)DT);
1830 R__b.ReadStaticArray((
double*)ED);
1831 if(R__v > 2) R__b >>
ifoID;
1835 R__b >> *
reinterpret_cast<Int_t*
>(ptr_polarization);
1841 HRSS.Streamer(R__b);
1842 ISNR.Streamer(R__b);
1843 FREQ.Streamer(R__b);
1844 BAND.Streamer(R__b);
1845 TIME.Streamer(R__b);
1846 TDUR.Streamer(R__b);
1847 if(wfSAVE==1||wfSAVE==3) {
1849 vector<int> &R__stl = IWFID;
1853 R__stl.reserve(R__n);
1854 for (R__i = 0; R__i < R__n; R__i++) {
1857 R__stl.push_back(R__t);
1861 vector<wavearray<double>*> &R__stl = IWFP;
1865 R__stl.reserve(R__n);
1866 for (R__i = 0; R__i < R__n; R__i++) {
1868 R__t->Streamer(R__b);
1869 R__stl.push_back(R__t);
1873 if(wfSAVE==2||wfSAVE==3) {
1875 vector<int> &R__stl = RWFID;
1879 R__stl.reserve(R__n);
1880 for (R__i = 0; R__i < R__n; R__i++) {
1883 R__stl.push_back(R__t);
1887 vector<wavearray<double>*> &R__stl = RWFP;
1891 R__stl.reserve(R__n);
1892 for (R__i = 0; R__i < R__n; R__i++) {
1894 R__t->Streamer(R__b);
1895 R__stl.push_back(R__t);
1902 R__b.CheckByteCount(R__s, R__c, detector::IsA());
1904 R__c = R__b.WriteVersion(detector::IsA(), kTRUE);
1905 TNamed::Streamer(R__b);
1906 R__b.WriteArray(Name, 16);
1908 R__b.WriteArray(Rv, 3);
1909 R__b.WriteArray(Ex, 3);
1910 R__b.WriteArray(Ey, 3);
1911 R__b.WriteArray(DT, 9);
1912 R__b.WriteArray(ED, 5);
1919 HRSS.Streamer(R__b);
1920 ISNR.Streamer(R__b);
1921 FREQ.Streamer(R__b);
1922 BAND.Streamer(R__b);
1923 TIME.Streamer(R__b);
1924 TDUR.Streamer(R__b);
1925 if(wfSAVE==1||wfSAVE==3) {
1927 vector<int> &R__stl = IWFID;
1928 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1931 vector<int>::iterator R__k;
1932 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1938 vector<wavearray<double> > IWF(IWFP.size());
1939 for(
int i=0;
i<IWFP.size();
i++) IWF[
i] = *IWFP[
i];
1940 vector<wavearray<double> > &R__stl = IWF;
1941 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1944 vector<wavearray<double> >::iterator R__k;
1945 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1951 if(wfSAVE==2||wfSAVE==3) {
1953 vector<int> &R__stl = RWFID;
1954 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1957 vector<int>::iterator R__k;
1958 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1964 vector<wavearray<double> > RWF(RWFP.size());
1965 for(
int i=0;
i<RWFP.size();
i++) RWF[
i] = *RWFP[
i];
1966 vector<wavearray<double> > &R__stl = RWF;
1967 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1970 vector<wavearray<double> >::iterator R__k;
1971 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1979 R__b.SetByteCount(R__c, kTRUE);
wavearray< double > t(hp.size())
void writeFilter(const char *)
param: file name
detectorParams getDetectorParams()
virtual void resize(unsigned int)
virtual size_t size() const
double getWFfreq(char atype='S')
static double g(double e)
void setFpFx(double, double=0., double=180., double=0., double=360.)
param - step on phi and theta param - theta begin param - theta end param - phi begin param - phi end...
printf("total live time: non-zero lags = %10.1f \n", liveTot)
size_t setsnr(wavearray< double > &, std::vector< double > *, std::vector< double > *, double=5., double=8.)
virtual void rate(double r)
plot hist2D SetName("WSeries-1")
wavearray< double > a(hp.size())
double getWFtime(char atype='S')
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
std::vector< delayFilter > filter
double getTheta(size_t i)
std::slice getSlice(double n)
size_t setFilter(size_t, double=0., size_t=0)
param: filter length param: filter phase delay in degrees param: number of up-sample layers return fi...
bool setrms(netcluster *, size_t=0)
watplot p1(const_cast< char * >("TFMap1"))
virtual void start(double s)
size_t minDFindex(delayFilter &F)
void delay(double, double)
param: theta param: phi
virtual double mean() const
void init()
param: Rx,Ry,Rz in ECEF frame
double getNoise(size_t, int=-1)
param: wavelet layer index (frequency) param: wavelet time index if param 2 is specified - return noi...
void readFilter(const char *)
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
void set(double x, double y)
detector & operator=(const detector &)
void bandPass1G(double f1=0., double f2=0.)
wavearray< double > hot[2]
double getwave(int, WSeries< double > &, char='W', size_t=0)
param: cluster ID param: WSeries where to put the waveform return: noise rms
void GeocentricToGeodetic(double X, double Y, double Z, double &latitude, double &longitude, double &elevation)
virtual int convertF2L(int, int)
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
POLARIZATION polarization
detector & operator>>(detector &)
std::vector< wavearray< double > * > IWFP
virtual int convertO2F(int, int)
netpixel * getPixel(size_t n, size_t i)
void GeodeticToGeocentric(double latitude, double longitude, double elevation, double &X, double &Y, double &Z)
virtual Wavelet * Clone() const
return: Wavelet* - duplicate of *this, allocated on heap
virtual void waveSort(DataType_t **pp, size_t l=0, size_t r=0) const
WSeries< double > waveForm
watplot p2(const_cast< char * >("TFMap2"))
double fabs(const Complex &x)
std::vector< short > index
virtual int getOffset(int, int)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Meyer< double > S(1024, 2)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double getTau(double, double)
param: source theta,phi angles in degrees
void set(size_t i, double a)
param: sky index param: value to set
double getwave(int, netcluster &, char, size_t)
param: no parameters
WaveDWT< DataType_t > * pWavelet
void GetCartesianComponents(double u[3], double Alt, double Az, double Lat, double Lon)
double getdata(char type='R', size_t n=0)
std::vector< float > value
virtual void resize(unsigned int)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
size_t setsim(WSeries< double > &, std::vector< double > *, double=5., double=8., bool saveWF=false)
detector & operator<<(detector &)
void setTau(double, double=0., double=180., double=0., double=360.)
param - step on phi and theta param - theta begin param - theta end param - phi begin param - phi end...
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
virtual WSeries< double > white(double, int, double=0., double=0.)
what it does: each wavelet layer is devided into k intervals.
bool setdata(double a, char type='R', size_t n=0)