30 #include "TRotation.h" 31 #include "Math/Rotation3D.h" 32 #include "Math/Vector3Dfwd.h" 41 extern const double _rH1[3] = {-2.161414928e6, -3.834695183e6, 4.600350224e6};
42 extern const double _eH1x[3] = {-0.223891216, 0.799830697, 0.556905359};
43 extern const double _eH1y[3] = {-0.913978490, 0.026095321, -0.404922650};
46 extern const double _rH2[3] = {-2.161414928e6, -3.834695183e6, 4.600350224e6};
47 extern const double _eH2x[3] = {-0.223891216, 0.799830697, 0.556905359};
48 extern const double _eH2y[3] = {-0.913978490, 0.026095321, -0.404922650};
51 extern const double _rL1[3] = {-7.427604192e4, -5.496283721e6, 3.224257016e6};
52 extern const double _eL1x[3] = {-0.954574615, -0.141579994, -0.262187738};
53 extern const double _eL1y[3] = { 0.297740169, -0.487910627, -0.820544948};
61 extern const double _rG1[3] = {3.8563112e6, 6.665978e5, 5.0196406e6};
62 extern const double _eG1x[3] = {-0.445184239, 0.866534205, 0.225675575};
63 extern const double _eG1y[3] = {-0.626000687, -0.552167273, 0.550667271};
66 extern const double _rV1[3] = {4.5463741e6, 8.429897e5, 4.378577e6};
67 extern const double _eV1x[3] = {-0.700458309, 0.208487795, 0.682562083};
68 extern const double _eV1y[3] = {-0.053791331, -0.969082169, 0.240803326};
71 extern const double _rT1[3] = {-3.946409e6, 3.366259e6, 3.6991507e6};
72 extern const double _eT1x[3] = {0.648969405, 0.760814505, 0};
73 extern const double _eT1y[3] = {-0.443713769, 0.378484715, -0.812322234};
76 extern const double _rK1[3] = {-3.7770908e+06, 3.4846722e+06, 3.7650676e+06};
77 extern const double _eK1x[3] = {-0.3740476967, -0.8378767376, 0.3975561509};
78 extern const double _eK1y[3] = {0.7142972255, 0.01312009275, 0.6997194701};
81 extern const double _rA1[3] = {-2.3567784e6, 4.8970238e6, -3.3173147e6};
82 extern const double _eA1x[3] = {-9.01077021322091554e-01, -4.33659084587544319e-01, 0};
83 extern const double _eA1y[3] = {-2.25940560005277846e-01, 4.69469807139026196e-01, 8.53550797275327455e-01};
86 extern const double _rO1[3] = {4.392467e6, 0.9295086e6, 4.515029e6};
87 extern const double _eO1x[3] = {-0.644504130, 0.573655377, 0.50550364};
88 extern const double _eO1y[3] = {0., 0., 0.};
91 extern const double _rN1[3] = {4.64410999868e6, 1.04425342477e6, 4.23104713307e6};
92 extern const double _eN1x[3] = {-0.62792641437, 0.56480832712, 0.53544371484};
93 extern const double _eN1y[3] = {0., 0., 0.};
96 extern const double _rE1[3] = {4.37645395452e6, 4.75435044067e5, 4.59985274450e6};
97 extern const double _eE1x[3] = {-0.62792641437, 0.56480832712, 0.53544371484};
98 extern const double _eE1y[3] = {0., 0., 0.};
110 for(
int i=0;
i<3;
i++){
123 this->TFmap.rate(4096);
146 if(strstr(name,
"Virgo")) {pRv=
_rV1; pEx=
_eV1x; pEy=
_eV1y;xifo=
true;}
147 if(strstr(name,
"VIRGO")) {pRv=
_rV1; pEx=
_eV1x; pEy=
_eV1y;xifo=
true;}
148 if(strstr(name,
"GEO")) {pRv=
_rG1; pEx=
_eG1x; pEy=
_eG1y;xifo=
true;}
149 if(strstr(name,
"TAMA")) {pRv=
_rT1; pEx=
_eT1x; pEy=
_eT1y;xifo=
true;}
156 cout <<
"detector::detector - Error : detector " << name
157 <<
" is not in present in the builtin list" << endl;
163 for(
int i=0;
i<3;
i++){
164 this->Rv[
i] = pRv[
i];
165 this->Ex[
i] = pEx[
i];
166 this->Ey[
i] = pEy[
i];
168 if(strstr(name,
"A2")) this->rotate(45);
177 this->TFmap.rate(4096);
199 for(
int i=0;
i<3;
i++){
200 this->Rv[
i] = uRv[
i];
201 this->Ex[
i] = uEx[
i];
202 this->Ey[
i] = uEy[
i];
212 this->TFmap.rate(4096);
223 if(strlen(this->dP.name)>0)
return this->dP;
231 XYZVector iRv(this->Rv[0],this->Rv[1],this->Rv[2]);
232 TVector3 vRv(iRv.X(),iRv.Y(),iRv.Z());
235 XYZVector iEx(this->Ex[0],this->Ex[1],this->Ex[2]);
236 XYZVector iEy(this->Ey[0],this->Ey[1],this->Ey[2]);
237 XYZVector iEZ(0.,0.,1.);
239 TVector3 vEx(iEx.X(),iEx.Y(),iEx.Z());
240 TVector3 vEy(iEy.X(),iEy.Y(),iEy.Z());
241 TVector3 vEZ(iEZ.X(),iEZ.Y(),iEZ.Z());
248 TVector3 vEe=vEZ.Cross(vRv);
251 double cosExAngle=vEe.Dot(vEx);
252 if(cosExAngle>1.) cosExAngle=1.;
253 if(cosExAngle<-1.) cosExAngle=-1.;
254 double cosEyAngle=vEe.Dot(vEy);
255 if(cosEyAngle>1.) cosEyAngle=1.;
256 if(cosEyAngle<-1.) cosEyAngle=-1.;
258 double ExAngle=acos(cosExAngle)*
rad2deg;
259 double EyAngle=acos(cosEyAngle)*
rad2deg;
264 if(vEn.Dot(vRv)<0) ExAngle*=-1;
266 if(vEn.Dot(vRv)<0) EyAngle*=-1;
269 ExAngle=fmod(90-ExAngle,360.);
271 EyAngle=fmod(90-EyAngle,360.);
273 double latitude,longitude,elevation;
278 sprintf(idP.name,
"%s",this->Name);
279 idP.latitude = latitude*
rad2deg;
280 idP.longitude = longitude*
rad2deg;
281 idP.elevation = elevation;
295 double si = sin(a*
PI/180.);
296 double co = cos(a*
PI/180.);
298 for(
int i=0;
i<3;
i++) {
312 for(
int i=0;
i<3;
i++) axy+=ax[
i]*ay[
i];
316 for(
int i=0;
i<3;
i++) {aw[
i]=-ax[
i]*axy+ay[
i]; aww+=aw[
i]*aw[
i];}
317 for(
int i=0;
i<3;
i++) aw[
i]/=sqrt(aww);
318 for(
int i=0;
i<3;
i++) this->Ex[
i]=ax[
i]*co+aw[
i]*si;
322 for(
int i=0;
i<3;
i++) {aw[
i]=-ax[
i]+ay[
i]*axy; aww+=aw[
i]*aw[
i];}
323 for(
int i=0;
i<3;
i++) aw[
i]/=sqrt(aww);
324 for(
int i=0;
i<3;
i++) this->Ey[
i]=ay[
i]*co+aw[
i]*si;
329 if(strlen(this->dP.name)>0) {
330 this->dP.AzX = fmod(this->dP.AzX-a,360.);
331 this->dP.AzY = fmod(this->dP.AzY-a,360.);
348 for (
int i=0;
i<
n;
i++) {
356 for (
int i=0;
i<
n;
i++) {
372 for(
int i=0;
i<3;
i++){
373 this->Rv[
i] = value.
Rv[
i];
374 this->Ex[
i] = value.
Ex[
i];
375 this->Ey[
i] = value.
Ey[
i];
402 double tsRate =TFmap.wavearray<double>::rate();
404 waveForm.resize(
size_t(tsRate));
405 waveForm.rate(tsRate);
406 waveForm.setWavelet(*(TFmap.pWavelet));
423 for (
int i=0;
i<IWFP.size();
i++) {
428 for(
int i=0;
i<value.
IWFP.size();
i++) {
430 *wf = *value.
IWFP[
i];
450 for (
int i=0;
i<value.
IWFP.size();
i++) {
455 for(
int i=0;
i<IWFP.size();
i++) {
458 value.
IWFP.push_back(wf);
470 DT[0] = Ex[0]*Ex[0]-Ey[0]*Ey[0];
471 DT[1] = Ex[0]*Ex[1]-Ey[0]*Ey[1];
472 DT[2] = Ex[0]*Ex[2]-Ey[0]*Ey[2];
475 DT[4] = Ex[1]*Ex[1]-Ey[1]*Ey[1];
476 DT[5] = Ex[1]*Ex[2]-Ey[1]*Ey[2];
480 DT[8] = Ex[2]*Ex[2]-Ey[2]*Ey[2];
482 if (strcmp(Name,
"O1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
483 if (strcmp(Name,
"N1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
484 if (strcmp(Name,
"E1")==0)
for (
int i=0;
i<9;
i++) DT[
i]*=2;
494 theta *=
PI/180.; phi *=
PI/180.; psi *=
PI/180.;
496 double cT = cos(theta);
497 double sT = sin(theta);
498 double cP = cos(phi);
499 double sP = sin(phi);
519 fp = (cT*cP*d11 + cT*sP*d21 - sT*d31)*cT*cP
520 + (cT*cP*d12 + cT*sP*d22 - sT*d32)*cT*sP
521 - (cT*cP*d13 + cT*sP*d23 - sT*d33)*sT
523 - (cP*d22-sP*d12)*cP;
525 fx = -(cT*cP*d11 + cT*sP*d21 - sT*d31)*sP
526 + (cT*cP*d12 + cT*sP*d22 - sT*d32)*cP
527 + (cP*d21-sP*d11)*cT*cP
528 + (cP*d22-sP*d12)*cT*sP
529 - (cP*d23-sP*d13)*sT;
534 a = fp*cos(2*psi)+fx*sin(2*psi);
535 b = -fp*sin(2*psi)+fx*cos(2*psi);
542 fp = -(cT*cP*d11 + cT*sP*d21 - sT*d31)*cT*cP
543 - (cT*cP*d12 + cT*sP*d22 - sT*d32)*cT*sP
544 + (cT*cP*d13 + cT*sP*d23 - sT*d33)*sT
546 - (cP*d22-sP*d12)*cP;
564 double a,b,rms,fl,fh;
565 double R = this->TFmap.
rate();
566 int L =
int(this->TFmap.maxLayer()+1);
568 waveForm.setWavelet(*(TFmap.pWavelet));
570 rms = wc.
getwave(ID,waveForm,atype,index);
572 if(rms==0.)
return rms;
576 waveBand = waveForm; waveBand = 0.;
577 fh = waveBand.gethigh();
578 fl = waveBand.getlow();
579 l =
int(this->TFmap.getLevel()) -
int(waveBand.getLevel());
583 if(l<0) { waveBand.Inverse(-l); }
584 else { waveBand.Forward(l); }
586 a = waveBand.start()*R;
587 i =
int(a + ((a>0) ? 0.1 : -0.1));
588 k = waveBand.size(); j = 0;
589 if(i<0) { k +=
i; j = -
i; i = 0; }
590 if((i/L)*L != i) cout<<
"detector::getwave() time mismatch: "<<L<<
" "<<i<<
"\n";
592 waveBand.cpf(this->TFmap,k,i,j);
595 n =
int(2*L*fl/R+0.1)-1;
596 m =
int(2*L*fh/R+0.1)+1;
599 if(m<=n) { n=0; m=
L; }
607 for(k=n; k<
m; k++) { w.
getLayer(x,k); waveBand.putLayer(x,k); }
612 waveForm *= atype==
'w' ? 1./sqrt(this->
rate/R) : 1.;
617 double sTARt= waveForm.start();
618 size_t I = waveForm.size();
620 double sum = waveForm.data[
M]*waveForm.data[
M];
625 for(i=1; i<
int(M); i++) {
626 a = waveForm.data[M-
i];
627 b = waveForm.data[M+
i];
629 if(sum/hrss > 0.999 && i/R>0.05)
break;
632 if(n <
int(M-1)) i = size_t(n);
637 waveForm.cpf(w,2*i,M-i);
638 waveForm.resize(2*i);
639 waveForm.start(sTARt+(M-i)/R);
642 waveBand.
cpf(w,2*i,M-i);
643 waveBand.resize(2*i);
644 waveBand.start(sTARt+(M-i)/R);
658 size_t n = SM.
size();
664 z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
681 size_t n = SM.
size();
687 z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
700 double x = theta*
PI/180.;
701 double y = phi*
PI/180.;
702 double z = Rv[0]*sin(x)*cos(y) + Rv[1]*sin(x)*sin(y) + Rv[2]*cos(x);
712 skymap Sp(sms,t1,t2,p1,p2);
713 skymap Sx(sms,t1,t2,p1,p2);
714 size_t n = Sp.
size();
739 size_t n = Sp.
size();
762 cout<<
"wseries::setFilter(): not applicable to WDM TFmaps\n";
765 size_t i,
j,
k,
n,ii,jj;
766 size_t M = TFmap.maxLayer()+1;
767 size_t L = TFmap.getLevel();
768 size_t m = M*TFmap.pWavelet->m_H;
769 size_t N = M*(1<<upL);
772 std::vector<delayFilter>
F;
774 for(i=0; i<
K; i++) v.
value[i] = 0.;
775 for(i=0; i<M; i++) F.push_back(v);
798 cout<<
"wsize: "<<w.
size()<<endl;
803 double* pb = (
double * )malloc(m*
sizeof(
double));
804 double** pp = (
double **)malloc(m*
sizeof(
double*));
805 for(i=0; i<
m; i++) pp[i] = w.
data + i + (
int(j) -
int(m/2));
836 for(k=2; k<W.
size(); k+=2) {
866 for(k=m-1; k>=0; k--) {
868 if(abs(offst) >32767)
continue;
869 inDex = short(offst);
870 vaLue = float(pb[pp[k]-p0]);
871 if(
fabs(vaLue)<1.
e-4)
break;
873 offst -= (offst/
M)*M;
874 if(offst<0) offst +=
M;
879 offst -= (offst/
M)*M;
880 if(offst<0) offst +=
M;
881 if(offst !=
int(s)) cout<<
"setFilter error 2: "<<offst<<
" "<<s<<endl;
885 for(
size_t kk=0; kk<
K; kk++)
888 if(jj>=K) {cout<<jj<<endl;
continue;}
890 F[ii].value[jj] = vaLue;
891 F[ii].index[jj] = -inDex;
913 sum += F[
i].value[
k]*F[
i].value[
k];
914 if(n && F[i].
value[k] == 0.)
printf(
"%4d %4d %d4 \n",
int(n),
int(i),
int(k));
915 if(v.
value[k] != F[i].value[k] ||
916 v.
index[k] != F[i].index[k]) cout<<
"setFilter error 4\n";
919 if(sum<0.97)
printf(
"%4d %4d %8.5f \n",
int(n),
int(i),sum);
936 std::vector<delayFilter>().swap(
filter);
939 for(
size_t k=0;
k<
K;
k++) {
952 if ( (fp=fopen(fname,
"wb")) == NULL ) {
953 cout <<
" DumpBinary() error : cannot open file " << fname <<
". \n";
957 size_t M = size_t(TFmap.maxLayer()+1);
959 size_t N = this->nDFS;
960 size_t n = K *
sizeof(float);
961 size_t m = K *
sizeof(short);
966 fwrite(&K,
sizeof(
size_t), 1, fp);
967 fwrite(&M,
sizeof(
size_t), 1, fp);
968 fwrite(&N,
sizeof(
size_t), 1, fp);
976 fwrite(value.data, n, 1, fp);
977 fwrite(index.
data, m, 1, fp);
991 if ( (fp=fopen(fname,
"rb")) == NULL ) {
992 cout <<
" DumpBinary() error : cannot open file " << fname <<
". \n";
1000 fread(&K,
sizeof(
size_t), 1, fp);
1001 fread(&M,
sizeof(
size_t), 1, fp);
1002 fread(&N,
sizeof(
size_t), 1, fp);
1004 size_t n = K *
sizeof(float);
1005 size_t m = K *
sizeof(short);
1012 this->clearFilter();
filter.reserve(N*M);
1016 for(k=0; k<
K; k++) {
1017 v.
value.push_back(0.);
1018 v.
index.push_back(0);
1021 for(i=0; i<
N; i++) {
1022 for(j=0; j<
M; j++) {
1023 fread(value.data, n, 1, fp);
1024 fread(index.
data, m, 1, fp);
1025 for(k=0; k<
K; k++) {
1056 cout<<
"wseries::delay(): not applicable to WDM TFmaps\n";
1066 int n =
int(t*TFmap.wavearray<double>::rate());
1067 int m = n>0 ? (n+N/2-1)/N : (n-N/2)/
N;
1068 int l = TFmap.pWavelet->m_H/4+2;
1074 cout<<
"delay="<<t<<
" m="<<m<<
" n="<<n<<
" M="<<M<<
" K="<<K<<
" N="<<N<<endl;
1077 double*
A = (
double*)malloc(K*
sizeof(
double));
1078 int* I = (
int*)malloc(K*
sizeof(
int));
1080 for(i=0; i<
M; i++) {
1085 je = m<0 ? w.
size() +(m-
l)*M : w.
size() -l*
M;
1090 for(j=jb; j<je; j+=
M) {
1092 TFmap.data[
j] = A[0]*p[I[0]] + A[1]*p[I[1]] + A[2]*p[I[2]] + A[3]*p[I[3]]
1093 + A[4]*p[I[4]] + A[5]*p[I[5]] + A[6]*p[I[6]] + A[7]*p[I[7]]
1094 + A[8]*p[I[8]] + A[9]*p[I[9]] + A[10]*p[I[10]] + A[11]*p[I[11]]
1095 + A[12]*p[I[12]] + A[13]*p[I[13]] + A[14]*p[I[14]] + A[15]*p[I[15]];
1100 for(j=jb; j<je; j+=
M) {
1102 TFmap.data[
j] = A[0]*p[I[0]] + A[1]*p[I[1]] + A[2]*p[I[2]] + A[3]*p[I[3]]
1103 + A[4]*p[I[4]] + A[5]*p[I[5]] + A[6]*p[I[6]] + A[7]*p[I[7]]
1104 + A[8]*p[I[8]] + A[9]*p[I[9]] + A[10]*p[I[10]] + A[11]*p[I[11]]
1105 + A[12]*p[I[12]] + A[13]*p[I[13]] + A[14]*p[I[14]] + A[15]*p[I[15]]
1106 + A[16]*p[I[16]] + A[17]*p[I[17]] + A[18]*p[I[18]] + A[19]*p[I[19]]
1107 + A[20]*p[I[20]] + A[21]*p[I[21]] + A[22]*p[I[22]] + A[23]*p[I[23]]
1108 + A[24]*p[I[24]] + A[25]*p[I[25]] + A[26]*p[I[26]] + A[27]*p[I[27]]
1109 + A[28]*p[I[28]] + A[29]*p[I[29]] + A[30]*p[I[30]] + A[31]*p[I[31]];
1114 for(j=jb; j<je; j+=
M) {
1118 while(k-- > 0) { *q += *(A++) * p[*(I++)]; }
1131 if(!nRMS.size())
return 0.;
1132 if(
int(I)>TFmap.maxLayer())
return 0.;
1136 size_t N = nRMS.maxLayer()+1;
1137 size_t M = TFmap.maxLayer()+1;
1139 size_t m = (I+1)*(N/M);
1140 double t = TFmap.
start();
1141 double T = nRMS.start();
1142 slice S = TFmap.getSlice(I);
1143 double rATe = TFmap.wrate();
1145 double rESn = nRMS.frequency(1)-nRMS.frequency(0);
1150 cout<<
"detector::getNoise(): invalid noise rms array nRMS\n";
1154 size_t L = size_t(TFmap.getlow()/rESn+0.1);
1155 size_t H = size_t(TFmap.gethigh()/rESn+0.1);
1156 if(n>=H || m<=L)
return 0.;
1161 cout<<
"detector::getNoise():: invalid noise rms array nRMS\n";
1165 S = nRMS.getSlice(n);
1166 size_t K = S.
size();
1174 for(l=n; l<
m; l++) {
1175 S = nRMS.getSlice(l);
1176 g = (n<L || m>H) ? 1.e10 : 1.;
1177 for(j=0; j<rms.
size(); j++) {
1183 for(j=0; j<rms.
size(); j++) {
1184 rms.
data[
j] = sqrt((
double(m)-
double(n))/rms.
data[j]);
1196 int inRMS =
int((t-nRMS.start())*nRMS.rate());
1197 int inVAR = nVAR.size() ?
int((t-nVAR.start())*nVAR.rate()) : 0;
1199 if(inRMS < 0) inRMS = 0;
1200 else if(
size_t(inRMS)>=K &&
K) inRMS = K-1;
1202 if(inVAR <= 0) inVAR = 0;
1203 else if(
size_t(inVAR)>=nVAR.size()) inVAR = nVAR.size()-1;
1205 for(l=n; l<
m; l++) {
1206 S = nRMS.getSlice(l);
1207 g = (n<L || m>H) ? 1.e10 : 1.;
1212 RMS /= double(m)-double(n);
1215 if(!nVAR.size() || rESn*n<nVAR.getlow() || rESn*m>nVAR.gethigh())
return RMS;
1217 return RMS*double(nVAR.data[inVAR]);
1228 size_t M = wc->
size();
1230 if(!M)
return false;
1232 if(!nRMS.size())
return false;
1234 size_t max_layer = nRMS.maxLayer();
1239 int K = nRMS.size()/(max_layer+1);
1240 double To = nRMS.start();
1241 double Ro = nRMS.wrate();
1242 double Fo = nRMS.frequency(0);
1243 double dF = nRMS.frequency(1)-Fo;
1244 double fl = wc->
getlow()-0.1;
1245 double fh = wc->
gethigh()+0.1;
1248 Fo = Fo==0. ? 0.5 : 0.;
1254 int(p->
rate/Ro+0.01) < 1 ||
1256 cout<<
"detector::setrms() - illegal pixel from zero level\n";
1262 n = size_t(f/dF+0.6);
1263 F = (x+1)*p->
rate/2.;
1264 m =
size_t(F/dF+0.6);
1265 if(m>max_layer) m=max_layer+1;
1270 if(k>=K) k -= k ? 1 : 0;
1271 if(k<0 || n>=m || k>=K) {
1272 cout<<
"detector::setrms() - invalid input: ";
1273 cout<<k<<
" "<<n<<
" "<<m<<
" "<<f<<
" "<<F<<
" "<<t<<endl;
1274 cout<<p->
frequency<<
" "<<p->
rate/2.<<
" "<<dF<<
" "<<fl<<endl;
1281 for(j=n; j<
m; j++) {
1282 S = nRMS.getSlice(j);
1283 g = (f<fl || F>fh) ? 1.e10 : 1.;
1290 ff = f<nVAR.getlow() ? nVAR.getlow() :
f;
1291 if(ff>=nVAR.gethigh()) ff=nVAR.gethigh();
1292 FF = F>nVAR.gethigh() ? nVAR.gethigh() :
F;
1293 if(FF<=nVAR.getlow()) FF=nVAR.getlow();
1294 ff = 2*(FF-ff)/p->
rate;
1295 FF = nVAR.get(t,0.5/p->
rate);
1300 if(r <= 0) cout<<
"detector:setrms() error!\n";
1313 double dF = TFmap.frequency(1)-TFmap.frequency(0);
1314 double fl =
fabs(f1)>0. ?
fabs(f1) : this->TFmap.getlow();
1315 double fh =
fabs(f2)>0. ?
fabs(f2) : this->TFmap.gethigh();
1316 size_t n = TFmap.pWavelet->m_WaveType==
WDMT ? size_t((fl+dF/2.)/dF+0.1) : size_t(fl/dF+0.1);
1317 size_t m = TFmap.pWavelet->m_WaveType==
WDMT ? size_t((fh+dF/2.)/dF+0.1)-1 : size_t(fh/dF+0.1)-1;
1318 size_t M = this->TFmap.maxLayer()+1;
1323 for(i=0; i<
int(M); i++) {
1325 if((f1>=0 && i>=n) && (f2>=0 && i<=
m))
continue;
1326 if((f1<0 && i<n) || (f2<0 && i>
m))
continue;
1327 if((f1<0 && f2>=0 && i<n))
continue;
1328 if((f1>=0 && f2<0 && i>=m))
continue;
1330 this->TFmap.getLayer(w,i+0.01); w=0.; this->TFmap.putLayer(w,i+0.01);
1331 this->TFmap.getLayer(w,-i-0.01); w=0.; this->TFmap.putLayer(w,-i-0.01);
1379 int j,nstop,nstrt,
n,
m,J;
1381 double a,b,
T,E,
D,H,
f;
1382 double R = this->
rate;
1383 size_t K = pT->size();
1386 bool pOUT = dT>0. ? false :
true;
1390 if(wi.
size() != TFmap.size()) {
1391 cout<<
"setsim(): mismatch between MDC size " 1392 <<wi.
size()<<
" and data size "<<TFmap.size()<<
"\n";
1397 if(!this->nRMS.size() || !this->TFmap.size())
return 0;
1399 if(this->HRSS.size() !=
K) {
1400 this->HRSS.resize(K);
1401 this->ISNR.resize(K);
1402 this->FREQ.resize(K);
1403 this->BAND.resize(K);
1404 this->TIME.resize(K);
1405 this->TDUR.resize(K);
1423 if(!W.
white(this->nRMS,1)) {
1424 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1426 if(!W.
white(this->nRMS,-1)) {
1427 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1436 if(!W.
white(this->nRMS)) {
1437 cout<<
"detector::setsim error: invalid noise array\n";
exit(1);
1446 size_t N = w.
size();
1447 double rate = w.wavearray<double>::rate();
1448 double bgps = w.
start()+offset+1.;
1449 double sgps = w.
start()-offset+N/rate-1.;
1451 for(k=0; k<
K; k++) {
1453 T = (*pT)[
k] - w.
start();
1455 nstrt =
int((T - dT)*rate);
1456 nstop =
int((T + dT)*rate);
1457 if(nstrt<=0) nstrt = 0;
1458 if(nstop>=
int(N)) nstop =
N;
1459 if(nstop<=0)
continue;
1460 if(nstrt>=
int(N))
continue;
1463 for(j=nstrt; j<nstop; j++) {
1472 nstrt =
int((T - dT)*rate);
1473 nstop =
int((T + dT)*rate);
1474 if(nstrt<=0) nstrt = 0;
1475 if(nstop>=
int(N)) nstop =
N;
1476 if(nstop<=0)
continue;
1477 if(nstrt>=
int(N))
continue;
1480 for(j=nstrt; j<nstop; j++) {
1491 if(n>=
int(S.size())) n = S.
size()-1;
1493 for(j=nstrt; j<nstop; j++) {
1494 a = w.
data[
j]*(j/rate-
T);
1503 if(T<bgps || T>sgps)
continue;
1505 this->ISNR.data[
k] = E;
1506 this->TIME.data[
k] =
T;
1507 this->TDUR.data[
k] =
D;
1508 this->HRSS.data[
k] = H;
1511 for(i=0; i<I; i++) {
1516 this->FREQ.data[
k] += f*a*
a;
1517 this->BAND.
data[
k] += f*f*a*
a;
1521 this->FREQ.
data[
k] /= E;
1522 a = this->FREQ.data[
k];
1523 b = this->BAND.data[
k]/E - a*
a;
1524 this->BAND.
data[
k] = b>0. ? sqrt(b) : 0.;
1528 double time = this->TIME.data[
k];
1529 double tdur = this->TDUR.data[
k];
1530 double tDur = 6*tdur;
1531 if (tDur > dT) tDur =
dT;
1532 int nStrt =
int((time-tDur-w.
start())*rate);
1533 int nStop =
int((time+tDur-w.
start())*rate);
1534 if(nStrt<0) nStrt = 0;
1535 if(nStop>
int(N)) nStop =
N;
1539 size_t I =
int(2.*dT*rate);
1543 for (j=0;j<I;j++) ms += ((OS+j>=0)&&(OS+j<N)) ? w.
data[OS+j]*w.
data[OS+j] : 0.;
1547 double sum = ((OS>=0)&&(OS<
N)) ? w.
data[
OS]*w.
data[
OS] : 0.;
1548 for(j=1; j<
int(I/2); j++) {
1549 a = ((OS-j>=0)&&(OS-j<N)) ? w.
data[OS-j] : 0.;
1550 b = ((OS+j>=0)&&(OS+j<
N)) ? w.
data[OS+
j] : 0.;
1552 if(sum/ms > 0.999)
break;
1557 if(nStrt<0) nStrt = 0;
1558 if(nStop>
int(N)) nStop = N;
1564 for(j=nStrt; j<nStop; j++) {
1570 double f_low = this->TFmap.getlow();
1571 double f_high = this->TFmap.gethigh();
1574 double Fs=((double)wf->
rate()/(double)wf->
size())/2.;
1575 for (j=0;j<wf->
size()/2;j+=2) {
1577 if((f<f_low)||(f>f_high)) {wf->
data[
j]=0.;wf->
data[j+1]=0.;}
1584 for(j=0;j<wf->
size();j++) {
1591 for(j=0;j<wf->
size();j++) {
1592 a = wf->
data[
j]*(j/rate-
T);
1599 this->ISNR.data[
k] = E;
1600 this->TIME.data[
k] =
T;
1601 this->TDUR.data[
k] =
D;
1616 for(j=nstrt; j<nstop; j++) {
1624 I =
int(2.*dT*rate);
1625 OS =
int((T - dT)*rate);
1628 for (j=0;j<I;j++) ms += ((OS+j>=0)&&(OS+j<N)) ? hot.
data[OS+j]*hot.
data[OS+j] : 0.;
1632 for(j=1; j<
int(I/2); j++) {
1633 a = ((OS-j>=0)&&(OS-j<N)) ? hot.
data[OS-j] : 0.;
1634 b = ((OS+j>=0)&&(OS+j<
N)) ? hot.
data[OS+
j] : 0.;
1636 if(sum/ms > 0.999)
break;
1639 nStrt =
int(T*rate)-
j;
1640 nStop =
int(T*rate)+
j;
1641 if(nStrt<0) nStrt = 0;
1642 if(nStop>
int(N)) nStop =
N;
1649 for(j=nStrt; j<nStop; j++) WF->
data[j-nStrt] = hot.
data[j];
1652 IWFID.push_back(-(k+1));
1657 printf(
"%3d T+-dT: %8.3f +-%5.3f, F+-dF: %4.0f +-%4.0f, SNR: %6.1e, hrss: %6.1e\n",
1658 int(M), T-bgps, D, FREQ.data[k], BAND.data[k], E, H);
1675 size_t K = pT->size();
1676 size_t N = wi.
size();
1685 for(k=0; k<
K; k++) {
1688 T = (*pT)[
k] - w.
start();
1690 nstrt =
int((T - dT)*rate);
1691 nstop =
int((T + dT)*rate);
1692 if(nstrt<=0) nstrt = 0;
1693 if(nstop>=
int(N)) nstop =
N;
1694 if(nstop<=0)
continue;
1695 if(nstrt>=
int(N))
continue;
1697 for(j=nstrt; j<nstop; j++) {
1713 if(!TFmap.size())
return;
1714 double R = this->TFmap.wavearray<double>::rate();
1715 double T = this->getTau(theta,phi);
1716 size_t n = size_t(
fabs(T)*R);
1717 size_t m = this->TFmap.size();
1722 if(T<0) TFmap.
cpf(w,m-n,0,n);
1723 else TFmap.cpf(w,m-n,n,0);
1732 if(!x.
size())
return;
1733 double R = x.
rate();
1734 double T = this->getTau(theta,phi);
1735 size_t n = size_t(
fabs(T)*R);
1736 size_t m = x.
size();
1741 if(T<0) x.
cpf(w,m-n,0,n);
1742 else x.
cpf(w,m-n,n,0);
1751 if(!x.
size())
return;
1752 double R = x.
rate();
1753 size_t n = size_t(
fabs(T)*R+0.5);
1754 size_t m = x.
size();
1759 if(T<0) x.
cpf(w,m-n,0,n);
1760 else x.
cpf(w,m-n,n,0);
1771 for(
size_t i=0;
i<wf->
size();
i++) {
1776 return E>0. ? wf->
start()+T/E/wf->
rate() : 0.;
1787 for(
size_t i=0;
i<wf->
size();
i+=2) {
1793 return E>0. ? 0.5*F*wf->
rate()/E/wf->
size() : 0.;
1803 if(theta_t>0) LAT=
'N';
else {LAT=
'S';theta_t=-theta_t;}
1804 int theta_d =
int(theta_t);
1805 int theta_m =
int((theta_t-theta_d)*60);
1806 float theta_s = (theta_t-theta_d-theta_m/60.)*3600.;
1810 if(phi_t>0) LON=
'E';
else {LON=
'W';phi_t=-phi_t;}
1811 int phi_d =
int(phi_t);
1812 int phi_m =
int((phi_t-phi_d)*60);
1813 float phi_s = (phi_t-phi_d-phi_m/60.)*3600.;
1816 cout <<
"----------------------------------------------" << endl;
1817 cout <<
"IFO : " << xdP.
name <<
" (Geographic Coordinates) " << endl;
1818 cout <<
"----------------------------------------------" << endl;
1820 cout <<
"latitude : " << xdP.
latitude <<
" longitude : " << xdP.
longitude << endl;
1822 cout <<
"LAT : " << LAT <<
" " << theta_d <<
", " << theta_m <<
", " << theta_s << endl;
1823 cout <<
"LON : " << LON <<
" " << phi_d <<
", " << phi_m <<
", " << phi_s << endl;
1827 cout <<
"Rv : " << Rv[0] <<
" " << Rv[1] <<
" " << Rv[2] <<
" " << endl;
1829 cout <<
"Ex : " << Ex[0] <<
" " << Ex[1] <<
" " << Ex[2] <<
" " << endl;
1831 cout <<
"Ey : " << Ey[0] <<
" " << Ey[1] <<
" " << Ey[2] <<
" " << endl;
1833 cout.precision(precision);
1834 cout <<
"Ex-North Angle Clockwise (deg) : " << xdP.
AzX << endl;
1835 cout <<
"Ey-North Angle Clockwise (deg) : " << xdP.
AzY << endl;
1845 void detector::Streamer(TBuffer &R__b)
1850 if (R__b.IsReading()) {
1851 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
1852 TNamed::Streamer(R__b);
1853 R__b.ReadStaticArray((
char*)Name);
1855 R__b.ReadStaticArray((
double*)Rv);
1856 R__b.ReadStaticArray((
double*)Ex);
1857 R__b.ReadStaticArray((
double*)Ey);
1858 R__b.ReadStaticArray((
double*)DT);
1859 R__b.ReadStaticArray((
double*)ED);
1860 if(R__v > 2) R__b >>
ifoID;
1864 R__b >> *
reinterpret_cast<Int_t*
>(ptr_polarization);
1870 HRSS.Streamer(R__b);
1871 ISNR.Streamer(R__b);
1872 FREQ.Streamer(R__b);
1873 BAND.Streamer(R__b);
1874 TIME.Streamer(R__b);
1875 TDUR.Streamer(R__b);
1876 if(wfSAVE==1||wfSAVE==3) {
1878 vector<int> &R__stl = IWFID;
1882 R__stl.reserve(R__n);
1883 for (R__i = 0; R__i < R__n; R__i++) {
1886 R__stl.push_back(R__t);
1890 vector<wavearray<double>*> &R__stl = IWFP;
1894 R__stl.reserve(R__n);
1895 for (R__i = 0; R__i < R__n; R__i++) {
1897 R__t->Streamer(R__b);
1898 R__stl.push_back(R__t);
1902 if(wfSAVE==2||wfSAVE==3) {
1904 vector<int> &R__stl = RWFID;
1908 R__stl.reserve(R__n);
1909 for (R__i = 0; R__i < R__n; R__i++) {
1912 R__stl.push_back(R__t);
1916 vector<wavearray<double>*> &R__stl = RWFP;
1920 R__stl.reserve(R__n);
1921 for (R__i = 0; R__i < R__n; R__i++) {
1923 R__t->Streamer(R__b);
1924 R__stl.push_back(R__t);
1931 R__b.CheckByteCount(R__s, R__c, detector::IsA());
1933 R__c = R__b.WriteVersion(detector::IsA(), kTRUE);
1934 TNamed::Streamer(R__b);
1935 R__b.WriteArray(Name, 16);
1937 R__b.WriteArray(Rv, 3);
1938 R__b.WriteArray(Ex, 3);
1939 R__b.WriteArray(Ey, 3);
1940 R__b.WriteArray(DT, 9);
1941 R__b.WriteArray(ED, 5);
1948 HRSS.Streamer(R__b);
1949 ISNR.Streamer(R__b);
1950 FREQ.Streamer(R__b);
1951 BAND.Streamer(R__b);
1952 TIME.Streamer(R__b);
1953 TDUR.Streamer(R__b);
1954 if(wfSAVE==1||wfSAVE==3) {
1956 vector<int> &R__stl = IWFID;
1957 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1960 vector<int>::iterator R__k;
1961 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1967 vector<wavearray<double> > IWF(IWFP.size());
1968 for(
int i=0;
i<IWFP.size();
i++) IWF[
i] = *IWFP[
i];
1969 vector<wavearray<double> > &R__stl = IWF;
1970 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1973 vector<wavearray<double> >::iterator R__k;
1974 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1980 if(wfSAVE==2||wfSAVE==3) {
1982 vector<int> &R__stl = RWFID;
1983 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1986 vector<int>::iterator R__k;
1987 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1993 vector<wavearray<double> > RWF(RWFP.size());
1994 for(
int i=0;
i<RWFP.size();
i++) RWF[
i] = *RWFP[
i];
1995 vector<wavearray<double> > &R__stl = RWF;
1996 int R__n=(&R__stl) ?
int(R__stl.size()) : 0;
1999 vector<wavearray<double> >::iterator R__k;
2000 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
2008 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)
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...
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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...
virtual Wavelet * Clone() const
return: Wavelet* - duplicate of *this, allocated on heap
bool setrms(netcluster *, size_t=0)
virtual void start(double s)
size_t minDFindex(delayFilter &F)
void delay(double, double)
param: theta param: phi
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
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...
watplot p2(const_cast< char *>("TFMap2"))
virtual size_t size() const
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)
watplot p1(const_cast< char *>("TFMap1"))
virtual int convertF2L(int, int)
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
printf("total live time: non-zero lags = %10.1f \, liveTot)
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)
WSeries< double > waveForm
virtual double mean() const
double fabs(const Complex &x)
virtual void waveSort(DataType_t **pp, size_t l=0, size_t r=0) const
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)