11 #include <xmmintrin.h>
12 #include <immintrin.h>
29 nRun(0), nLag(1),
nSky(0), mIFO(0), rTDF(0), Step(0.), Edge(0.),
gNET(0.), aNET(0.), iNET(0),
30 eCOR(0.), norm(1.), e2or(0.),
acor(sqrt(2.)), pOUT(false),
EFEC(true), local(true),
optim(true),
34 this->ifoList.clear();
35 this->ifoName.clear();
37 this->mdcList.clear();
38 this->livTime.clear();
39 this->mdcTime.clear();
40 this->mdcType.clear();
41 this->mdc__ID.clear();
69 size_t nIFO = this->ifoList.size();
72 cout<<
"network::getNetworkPixels(): "
73 <<
"invalid number of detectors or\n";
76 if(getifo(0)->
getTFmap()->w_mode != 1) {
77 cout<<
"network::getNetworkPixels(): invalid whitening mode.\n";
85 int i,
j,
k,
m,
n,NN,jj,nM,jE,jb,je,J,
K;
89 double R = pTF->
wrate();
90 double r = hTS->
rate();
95 int jB =
int(this->Edge*R+0.001);
96 if(jB&1) {cout<<
"getNetworkPixels(1): WDM parity violation\n";
exit(1);}
99 cout<<
"network::getNetworkPixels(): insufficient data edge length.\n";
114 for(n=0; n<
nIFO; n++) {
115 pdata[
n] = getifo(n)->getTFmap()->data;
120 double a,b,E,Ct,Cb,Ht,Hb;
122 if(hist) {pixeLHood = *
pTF; pixeLHood=-1.;}
123 if(this->veto.size() !=
M) {
126 short* pveto = this->veto.data;
128 this->wc_List[LAG].clear();
129 this->livTime[LAG] = 0.;
130 this->wc_List[LAG].setlow(pTF->
getlow());
131 this->wc_List[LAG].sethigh(pTF->
gethigh());
134 for(n=0; n<
nIFO; n++) {
135 b = this->getifo(n)->lagShift.data[LAG];
136 if(a>b) { a = b; nM =
n; }
139 for(n=0; n<
nIFO; n++) {
140 b = this->getifo(n)->lagShift.data[LAG];
141 K =
int((b-a)*R+0.001);
142 if(K&1) {cout<<
"getNetworkPixels(2): WDM parity violation\n";
exit(1);}
143 in[
n] = IN[
n] = K+jB;
158 if(jE&1) {cout<<
"getNetworkPixels(3): WDM parity violation\n";
exit(1);}
163 for(jj=0; jj<NN; jj++) {
166 pmap = MAP.
data+(jj+jB)*I;
167 for(n=0; n<
nIFO; n++) {
168 if(in[n] >= jE) in[
n] -= NN;
169 jb =
int(in[n]*r/R+0.01);
170 je =
int((in[n]+1)*r/R+0.01);
171 while(jb<je)
if(!pveto[jb++]) VETO=0.;
172 PDATA = &(pdata[
n][in[
n]*I]);
173 for(i=0; i<I; i++) pmap[i]+=*PDATA++;
179 if(pmap[i]<Eo || i<ib) pmap[
i]=0.;
180 if(pmap[i]>Em) pmap[
i]=Em+0.1;
185 for(jj=0; jj<NN; jj++) {
187 pmap = MAP.
data+(jj+jB)*I;
188 for(n=0; n<
nIFO; n++) {
189 if(IN[n] >= jE) IN[
n] -= NN;
191 for(n=0; n<5; n++) pp[n]=pmap+(n-2)*I;
192 for(i=ib; i<ie; i++) {
193 if((E=pp[2][i])<Eo)
continue;
194 Ct = pp[2][i+1]+pp[3][
i ]+pp[3][i+1];
195 Cb = pp[2][i-1]+pp[1][
i ]+pp[1][i-1];
197 Ht+= i<II? pp[4][i+2]+pp[3][i+2]:0.;
199 Hb+= i>1 ? pp[0][i-2]+pp[1][i-2]:0.;
207 for(n=0; n<
nIFO; n++) {
210 pix.
data[
n].asnr = sqrt(pdata[n][j]);
214 if(hist) hist->Fill(E);
215 if(hist) pixeLHood.data[
j] = E;
220 wc_List[LAG].append(pix);
223 for(n=0; n<
nIFO; n++) IN[n]++;
227 this->wc_List[LAG].start = pTF->
start();
228 this->wc_List[LAG].stop = pTF->
stop();
229 this->wc_List[LAG].rate = pTF->
rate();
230 this->livTime[LAG] = count/R;
232 if(nPix) this->setRMS();
248 __m128* _pa = (__m128*) pa;
249 __m128* _pb = (__m128*) pb;
250 __m128* _px = (__m128*) px;
252 for(
int i=0;
i<n*4;
i++) {
255 x.
data[
i]=
i*float(m)/float(n*4);
258 for(
int i=0;
i<
n;
i++) {
259 *(_pa+
i) = _mm_sqrt_ps(*(_px+
i));
260 *(_pa+
i) = _mm_div_ps(_mm_set1_ps(1.),*(_pa+
i));
261 *(_pb+
i) = _mm_rsqrt_ps(*(_px+
i));
262 _mm_storeu_ps(qq,*(_px+
i));
263 _mm_storeu_ps(sq,*(_pa+
i));
264 _mm_storeu_ps(rq,*(_pb+
i));
265 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[0],sq[0],rq[0]);
266 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[1],sq[1],rq[1]);
267 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[2],sq[2],rq[2]);
268 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[3],sq[3],rq[3]);
295 if(!this->wc_List[lag].
size())
return 0;
301 bool cirwave = mode==
'g' || mode==
'G' || mode==
'c' || mode==
'C';
302 bool linwave = mode==
'l' || mode==
'L' || mode==
's' || mode==
'S';
303 bool iotwave = mode==
'i' || mode==
'l' || mode==
'e' || mode==
'c' ||
304 mode==
'I' || mode==
'L' || mode==
'E' || mode==
'C';
305 bool psiwave = mode==
'l' || mode==
'e' || mode==
'p' ||
306 mode==
'L' || mode==
'E' || mode==
'P';
307 bool mureana = mode==
'i' || mode==
'e' || mode==
'c' ||
308 mode==
'r' || mode==
'p' || mode==
'b' ||
309 mode==
'l' || mode==
's' || mode==
'g';
310 bool rndwave = mode==
'r' || mode==
'R';
312 bool prior = this->gamma<0 ?
true :
false;
313 bool m_chirp = this->
optim ?
false : mureana;
315 if(!this->
optim) mureana =
true;
317 size_t nIFO = this->ifoList.size();
318 size_t ID = abs(iID);
321 cout<<
"network::likelihoodAVX(): invalid network.\n";
326 float gama = this->gamma*this->gamma*2./3.;
327 float deta =
fabs(this->
delta);
if(deta>1) deta=1;
328 float REG[2]; REG[0] = deta*sqrt(2);
331 static const __m128 _oo = _mm_set1_ps(1.
e-16);
332 static const __m128 _sm = _mm_set1_ps(-0.
f);
333 static const __m128 _En = _mm_set1_ps(En);
335 float aa,AA,Lo,Eo,Co,No,Ep,Lp,Np,Cp,Ec,Dc,To,Fo,Em,Lm,Rc,Mo,Mw,Eh;
336 float STAT,ee,EE,cc,
ff,FF,Lw,Ew,Cw,Nw,Gn,
rho,norm,ch,CH,Cr,Mp,
N;
338 size_t i,
j,
k,
l,
m,Vm,lm,V,V4,V44,
id,
K,
M;
341 short* mm = this->skyMask.data;
342 short* MM = skyMMcc.
data;
343 bool skymaskcc = (skyMaskCC.size()==
L);
360 std::vector<float*> _vtd;
361 std::vector<float*> _vTD;
362 std::vector<float*> _eTD;
363 std::vector<float*> _AVX;
364 std::vector<float*> _APN;
365 std::vector<float*> _DAT;
366 std::vector<float*> _SIG;
367 std::vector<float*> _NUL;
368 std::vector<float*> _TMP;
370 for(i=0; i<
NIFO; i++) {
372 ml[
i] = getifo(i)->index.
data;
373 FP[
i] = getifo(i)->fp.data;
374 FX[
i] = getifo(i)->fx.data;
377 ml[
i] = getifo(0)->index.data;
378 FP[
i] = getifo(0)->fp.data;
379 FX[
i] = getifo(0)->fx.data;
392 std::vector<int>* vint;
393 std::vector<int>* vtof;
399 std::map<int,float> vLr;
403 int csize = precision%65536;
404 int healpix = this->nSkyStat.getOrder();
405 int order = (precision-csize)/65536;
408 if(healpix && csize && order && order<healpix) {
410 for(
int l=0;l<rsm.
size();l++) {
411 int m = this->nSkyStat.getSkyIndex(rsm.
getTheta(l),rsm.
getPhi(l));
414 for(
int l=0;l<
L;l++) BB[l] = BB[l] ? 0 : 1;
421 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
422 cTo = pwc->
get((
char*)
"time",0,
'L',0);
426 id = size_t(cid.
data[k]+0.1);
428 if(pwc->
sCuts[
id-1] != -2)
continue;
430 vint = &(pwc->
cList[
id-1]);
431 vtof = &(pwc->
nTofF[
id-1]);
435 pI = wdmMRA.getXTalk(pwc,
id);
439 bBB = (V>this->wdmMRA.nRes*csize) ?
true :
false;
442 this->nSensitivity = 0.;
443 this->nAlignment = 0.;
444 this->nNetIndex = 0.;
445 this->nDisbalance = 0.;
446 this->nLikelihood = 0.;
447 this->nNullEnergy = 0.;
448 this->nCorrEnergy = 0.;
449 this->nCorrelation = 0.;
451 this->nEllipticity = 0.;
452 this->nPolarisation = 0.;
453 this->nProbability = 0.;
455 this->nAntenaPrior = 0.;
458 tsize = pix->
tdAmp[0].size();
459 if(!tsize || tsize&1) {
460 cout<<
"network::likelihoodWP() error: wrong pixel TD data\n";
466 if(!(V=pI.size()))
continue;
467 V4 = V + (V%4 ? 4 - V%4 : 0);
470 for(j=0; j<V4; j++) pJ.push_back(0);
480 for(i=0; i<
NIFO; i++) {
481 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
482 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vtd.push_back(ptmp);
483 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
484 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vTD.push_back(ptmp);
485 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
486 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _eTD.push_back(ptmp);
487 ptmp = (
float*)_mm_malloc((V4*3+16)*
sizeof(float),32);
488 for(j=0; j<(V4*3+16); j++) ptmp[j]=0; _APN.push_back(ptmp);
489 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
490 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _DAT.push_back(ptmp);
491 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
492 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _SIG.push_back(ptmp);
493 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
494 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _NUL.push_back(ptmp);
497 for(i=0; i<
NIFO; i++) {
498 pa[
i] = _vtd[
i] + (tsize/2)*V4;
499 pA[
i] = _vTD[
i] + (tsize/2)*V4;
500 pe[
i] = _eTD[
i] + (tsize/2)*V4;
501 pd[
i] = _DAT[
i]; pD[
i] = _DAT[
i]+V4;
502 ps[
i] = _SIG[
i]; pS[
i] = _SIG[
i]+V4;
503 pn[
i] = _NUL[
i]; pN[
i] = _NUL[
i]+V4;
506 this->a_00.resize(NIFO*V44); this->a_00=0.;
507 this->a_90.resize(NIFO*V44); this->a_90=0.;
508 this->rNRG.resize(V4); this->rNRG=0.;
509 this->pNRG.resize(V4); this->pNRG=0.;
511 __m128* _aa = (__m128*) this->a_00.data;
512 __m128* _AA = (__m128*) this->a_90.data;
515 float* p_et = (
float*)_mm_malloc(V4*
sizeof(
float),32);
516 for(j=0; j<V4; j++) p_et[j]=0; _AVX.push_back(p_et);
517 float* pMSK = (
float*)_mm_malloc(V44*
sizeof(
float),32);
518 for(j=0; j<V44; j++) pMSK[j]=0; _AVX.push_back(pMSK); pMSK[V4]=
nIFO;
519 float* p_fp = (
float*)_mm_malloc(V44*
sizeof(
float),32);
520 for(j=0; j<V44; j++) p_fp[j]=0; _AVX.push_back(p_fp);
521 float* p_fx = (
float*)_mm_malloc(V44*
sizeof(
float),32);
522 for(j=0; j<V44; j++) p_fx[j]=0; _AVX.push_back(p_fx);
523 float* p_si = (
float*)_mm_malloc(V4*
sizeof(
float),32);
524 for(j=0; j<V4; j++) p_si[j]=0; _AVX.push_back(p_si);
525 float* p_co = (
float*)_mm_malloc(V4*
sizeof(
float),32);
526 for(j=0; j<V4; j++) p_co[j]=0; _AVX.push_back(p_co);
527 float* p_uu = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
528 for(j=0; j<V4+16; j++) p_uu[j]=0; _AVX.push_back(p_uu);
529 float* p_UU = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
530 for(j=0; j<V4+16; j++) p_UU[j]=0; _AVX.push_back(p_UU);
531 float* p_vv = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
532 for(j=0; j<V4+16; j++) p_vv[j]=0; _AVX.push_back(p_vv);
533 float* p_VV = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
534 for(j=0; j<V4+16; j++) p_VV[j]=0; _AVX.push_back(p_VV);
535 float* p_au = (
float*)_mm_malloc(V4*
sizeof(
float),32);
536 for(j=0; j<V4; j++) p_au[j]=0; _AVX.push_back(p_au);
537 float* p_AU = (
float*)_mm_malloc(V4*
sizeof(
float),32);
538 for(j=0; j<V4; j++) p_AU[j]=0; _AVX.push_back(p_AU);
539 float* p_av = (
float*)_mm_malloc(V4*
sizeof(
float),32);
540 for(j=0; j<V4; j++) p_av[j]=0; _AVX.push_back(p_av);
541 float* p_AV = (
float*)_mm_malloc(V4*
sizeof(
float),32);
542 for(j=0; j<V4; j++) p_AV[j]=0; _AVX.push_back(p_AV);
543 float* p_uv = (
float*)_mm_malloc(V4*4*
sizeof(
float),32);
544 for(j=0; j<V4*4; j++) p_uv[j]=0; _AVX.push_back(p_uv);
545 float* p_ee = (
float*)_mm_malloc(V4*
sizeof(
float),32);
546 for(j=0; j<V4; j++) p_ee[j]=0; _AVX.push_back(p_ee);
547 float* p_EE = (
float*)_mm_malloc(V4*
sizeof(
float),32);
548 for(j=0; j<V4; j++) p_EE[j]=0; _AVX.push_back(p_EE);
549 float* pTMP=(
float*)_mm_malloc(V4*4*NIFO*
sizeof(
float),32);
550 for(j=0; j<V4*4*
NIFO; j++) pTMP[j]=0; _AVX.push_back(pTMP);
551 float* p_ni = (
float*)_mm_malloc(V4*
sizeof(
float),32);
552 for(j=0; j<V4; j++) p_ni[j]=0; _AVX.push_back(p_ni);
553 float* p_ec = (
float*)_mm_malloc(V4*
sizeof(
float),32);
554 for(j=0; j<V4; j++) p_ec[j]=0; _AVX.push_back(p_ec);
555 float* p_gn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
556 for(j=0; j<V4; j++) p_gn[j]=0; _AVX.push_back(p_gn);
557 float* p_ed = (
float*)_mm_malloc(V4*
sizeof(
float),32);
558 for(j=0; j<V4; j++) p_ed[j]=0; _AVX.push_back(p_ed);
559 float* p_rn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
560 for(j=0; j<V4; j++) p_rn[j]=0; _AVX.push_back(p_rn);
566 this->pList.push_back(pix);
569 for(i=0; i<
nIFO; i++) {
570 xx[
i] = 1./pix->
data[
i].noiserms;
575 for(i=0; i<
nIFO; i++) {
576 _APN[
i][V4*2+
j] =(float)xx[i]/rms;
577 for(l=0; l<tsize; l++) {
579 AA = pix->
tdAmp[
i].data[l+tsize];
580 _vtd[
i][l*V4+
j] = aa;
581 _vTD[
i][l*V4+
j] = AA;
582 _eTD[
i][l*V4+
j] = aa*aa+AA*AA;
595 STAT=-1.e12; lm=0; Em=Lm=ff=FF=0;
598 for(l=lb; l<=le; l++) {
600 if(bBB && !BB[l])
continue;
607 if (!skyMaskCC.data[lc])
continue;
611 if(aa > gama) ff += 1;
613 REG[1] = (FF*FF/(ff*ff+1.e-9)-1)*En;
618 for(l=lb; l<=le; l++) {
619 skyProb.data[
l] = -1.e12;
621 pnt_(v00, pa, ml, (
int)l, (
int)V4);
622 pnt_(v90, pA, ml, (
int)l, (
int)V4);
629 if(lb==le) _avx_saveGW_ps(ps,pS,V);
633 _mm256_storeu_ps(vvv,_CC);
638 CH = No/(nIFO*Mo+sqrt(Mo));
640 Co = Ec/(Ec+No*cc-Mo*(nIFO-1));
642 if(Cr<
netCC)
continue;
644 aa = Eo>0. ? Eo-No : 0.;
646 skyProb.data[
l] = this->
delta<0 ? aa : AA;
650 if(pMSK[j]==0)
continue;
652 ff += p_fp[
j]*p_et[
j];
653 FF += p_fx[
j]*p_et[
j];
655 ff = ee>0. ? ff/ee : 0.;
656 FF = ee>0. ? FF/ee : 0.;
657 this->nAntenaPrior.set(l, sqrt(ff+FF));
660 this->nSensitivity.set(l, sqrt(ff+FF));
661 this->nAlignment.set(l, ff>0 ? sqrt(FF/ff):0);
662 this->nLikelihood.set(l, Eo-No);
663 this->nNullEnergy.set(l, No);
664 this->nCorrEnergy.set(l, Ec);
665 this->nCorrelation.set(l,Co);
666 this->nSkyStat.set(l,AA);
667 this->nProbability.set(l, skyProb.data[l]);
668 this->nDisbalance.set(l,CH);
669 this->nNetIndex.set(l,cc);
670 this->nEllipticity.set(l,Cr);
671 this->nPolarisation.set(l,Mp);
674 if(AA>=STAT) {STAT=AA; lm=
l; Em=Eo-Eh;}
675 if(skyProb.data[l]>sky) sky=skyProb.data[
l];
681 D_snr = _avx_norm_ps(pd,pD,_AVX,V4);
682 S_snr = _avx_norm_ps(pS,pD,p_ec,V4);
687 _mm256_storeu_ps(vvv,_CC);
694 norm = Em>0 ? Ep/Em : 0;
695 norm*= norm<1 ? 1. : 1./norm;
698 rho= Ec>0 ? sqrt(Ec*Rc/2.) : 0.;
704 D_snr = _avx_norm_ps(pd,pD,_AVX,-V4);
705 N_snr = _avx_norm_ps(pn,pN,_AVX,-V4);
709 ch = (Np+Gn)/(N*nIFO);
713 if(le-lb) {lb=le=lm;
goto optsky;}
715 if(Lm<=0. || Em<=0. || Ec*Rc/cc<netEC || N<1) {
716 pwc->
sCuts[
id-1]=1; count=0;
717 pwc->
clean(
id);
continue;
724 vint = &(pwc->
cList[
id-1]);
725 for(j=0; j<vint->size(); j++) {
743 for(i=0; i<
nIFO; i++) {
744 pix->
setdata(
double(pd[i][j]),
'W',i);
745 pix->
setdata(
double(pD[i][j]),
'U',i);
746 pix->
setdata(
double(ps[i][j]),
'S',i);
747 pix->
setdata(
double(pS[i][j]),
'P',i);
753 if(!pix->
core)
continue;
754 if(p_gn[j]<=0)
continue;
759 if(!xpix->
core || p_gn[k]<=0 || xt.
CC[0]>2)
continue;
760 for(i=0; i<
nIFO; i++) {
768 if(p_ec[j]<=0)
continue;
774 if(!xpix->
core || p_ec[k]<=0 || xt.
CC[0]>2)
continue;
775 for(i=0; i<
nIFO; i++) {
787 for(j=1; j<=
nIFO; j++) {
788 if(S_snr[j]>Emax) Emax=S_snr[
j];
790 double Esub = S_snr.
data[0]-Emax;
791 Esub = Esub*(1+2*Rc*Esub/Emax);
792 Nmax = Gn+Np-N*(nIFO-1);
803 NETX (vtof->push_back(ml[0][lm]); ,
804 vtof->push_back(ml[1][lm]); ,
805 vtof->push_back(ml[2][lm]); ,
806 vtof->push_back(ml[3][lm]); ,
807 vtof->push_back(ml[4][lm]); ,
808 vtof->push_back(ml[5][lm]); ,
809 vtof->push_back(ml[6][lm]); ,
810 vtof->push_back(ml[7][lm]); )
813 if((wfsave)||(mdcListSize() && !lag)) {
814 if(this->getMRAwave(
id,lag,
'S',0,
true)) {
816 for(i=0; i<
nIFO; i++) {
817 pd = this->getifo(i);
818 pd->
RWFID.push_back(
id);
822 pd->
RWFP.push_back(wf);
825 if(this->getMRAwave(
id,lag,
's',0,
true)) {
827 for(i=0; i<
nIFO; i++) {
828 pd = this->getifo(i);
829 pd->
RWFID.push_back(-
id);
833 pd->
RWFP.push_back(wf);
838 Lw = Ew = To = Fo = Nw = ee = norm = 0.;
839 for(i=0; i<
nIFO; i++) {
844 this->getMRAwave(
id,lag,
'W',0);
845 this->getMRAwave(
id,lag,
'S',0);
846 for(i=0; i<
nIFO; i++) {
865 ch = (Nw+Gn)/(N*nIFO);
866 cc = ch>1 ? 1+(ch-1)*2*(1-Rc) : 1;
867 Cr = Ec*Rc/(Ec*Rc+(Dc+Nw+Gn)*cc-N*(nIFO-1));
869 Cp = Ec*Rc/(Ec*Rc+(Dc+Nw+Gn)-N*(nIFO-1));
872 pwc->
cData[
id-1].norm = norm*2;
873 pwc->
cData[
id-1].skyStat = 0;
874 pwc->
cData[
id-1].skySize = Mw;
875 pwc->
cData[
id-1].netcc = Cp;
876 pwc->
cData[
id-1].skycc = Cr;
877 pwc->
cData[
id-1].subnet = Esub/(Esub+Nmax);
878 pwc->
cData[
id-1].SUBNET = Co;
879 pwc->
cData[
id-1].likenet = Lw;
880 pwc->
cData[
id-1].netED = Nw+Gn+Dc-N*
nIFO;
881 pwc->
cData[
id-1].netnull = Nw+Gn;
882 pwc->
cData[
id-1].energy = Ew;
883 pwc->
cData[
id-1].likesky = Eo-No;
884 pwc->
cData[
id-1].enrgsky = Eo;
885 pwc->
cData[
id-1].netecor = Ec;
886 pwc->
cData[
id-1].normcor = Ec*Rc;
887 pwc->
cData[
id-1].netRHO = rho/sqrt(cc);
889 pwc->
cData[
id-1].cTime = To;
890 pwc->
cData[
id-1].cFreq = Fo;
891 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
892 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
893 pwc->
cData[
id-1].gNET = sqrt(ff+FF);
894 pwc->
cData[
id-1].aNET = sqrt(FF/ff);
895 pwc->
cData[
id-1].iNET = 0;
896 pwc->
cData[
id-1].nDoF =
N;
897 pwc->
cData[
id-1].skyChi2 = CH;
898 pwc->
cData[
id-1].Gnoise = Gn;
899 pwc->
cData[
id-1].iota = 0;
900 pwc->
cData[
id-1].psi = 0;
901 pwc->
cData[
id-1].ellipticity = 0.;
903 cc = pwc->
cData[
id-1].netcc;
905 printf(
"rho=%4.2f|%4.2f cc: %5.3f|%5.3f|%5.3f subnet=%4.3f|%4.3f \n",
906 rho,rho*sqrt(Cp),Co,Cp,Cr,pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
907 printf(
" L: %5.1f|%5.1f|%5.1f E: %5.1f|%5.1f|%5.1f N: %4.1f|%4.1f|%4.1f|%4.1f|%4.1f \n",
908 Lw,Lp,Lo,Ew,Ep,Eo,Nw,Np,Rc,Eh,No);
909 printf(
"id|lm %3d|%6d Vm|m=%3d|%3d|%3d|%3d T|F: %6.3f|%4.1f (t,p)=(%4.1f|%4.1f) \n",
910 int(
id),
int(lm),
int(V),
int(Mo),
int(Mw),
int(M),To,Fo,nLikelihood.getTheta(lm),nLikelihood.getPhi(lm));
911 cout<<
" L: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",S_snr[i]);}
912 cout<<
" E: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",D_snr[i]);}
913 cout<<
" N: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",N_snr[i]);}
914 cout<<endl<<
" dof|G|G+R ";
printf(
"%5.1f|%5.1f|%5.1f r[1]=%4.1f",N,Gn,Nw+Gn,REG[1]);
915 printf(
" norm=%3.1f chi2 %3.2f|%3.2f Rc=%3.2f, Dc=%4.1f\n",norm,ch,CH,Rc,Dc);
923 pwc->
p_Ind[
id-1].push_back(Mo);
925 std::vector<float> sArea;
926 pwc->
sArea.push_back(sArea);
927 pwc->
p_Map.push_back(sArea);
929 double var = norm*Rc*sqrt(Mo)*(1+
fabs(1-CH));
931 if(iID<=0 || ID==
id) {
932 getSkyArea(
id,lag,T,var);
936 pwc->
cData[
id-1].mchirp = 0;
937 pwc->
cData[
id-1].mchirperr = 0;
938 pwc->
cData[
id-1].tmrgr = 0;
939 pwc->
cData[
id-1].tmrgrerr = 0;
940 pwc->
cData[
id-1].chi2chirp = 0;
944 cc = Ec/(
fabs(Ec)+ee);
945 printf(
"mchirp : %d %g %.2e %.3f %.3f %.3f %.3f \n\n",
946 int(
id),cc,pwc->
cData[
id-1].mchirp,
947 pwc->
cData[
id-1].mchirperr, pwc->
cData[
id-1].tmrgr,
948 pwc->
cData[
id-1].tmrgrerr, pwc->
cData[
id-1].chi2chirp);
951 if(ID==
id && !
EFEC) {
952 this->nSensitivity.gps =
T;
953 this->nAlignment.gps =
T;
954 this->nDisbalance.gps =
T;
955 this->nLikelihood.gps =
T;
956 this->nNullEnergy.gps =
T;
957 this->nCorrEnergy.gps =
T;
958 this->nCorrelation.gps =
T;
959 this->nSkyStat.gps =
T;
960 this->nEllipticity.gps =
T;
961 this->nPolarisation.gps=
T;
962 this->nNetIndex.gps =
T;
965 pwc->
sCuts[
id-1] = -1;
992 if(!this->wc_List[lag].
size())
return 0;
994 size_t nIFO = this->ifoList.size();
997 cout<<
"network::subNetCut(): invalid network.\n";
1004 subnet =
fabs(subnet);
1005 subcut =
fabs(subcut);
1007 __m128 _En = _mm_set1_ps(En);
1008 __m128 _Es = _mm_set1_ps(Es);
1009 __m128 _oo = _mm_set1_ps(1.
e-12);
1010 __m128 _0 = _mm_set1_ps(0.);
1011 __m128 _05 = _mm_set1_ps(0.5);
1012 __m128 _1 = _mm_set1_ps(1.);
1017 float Lm,Em,Am,Lo,Eo,Co,Lr,Er,ee,em,To;
1018 float cc,aa,AA,rHo,stat,Ls,Ln,EE;
1022 short* mm = this->skyMask.data;
1035 for(i=0; i<
NIFO; i++) {
1037 ml[
i] = getifo(i)->index.data;
1038 FP[
i] = getifo(i)->fp.data;
1039 FX[
i] = getifo(i)->fx.data;
1042 ml[
i] = getifo(0)->index.data;
1043 FP[
i] = getifo(0)->fp.data;
1044 FX[
i] = getifo(0)->fx.data;
1049 std::vector<int> pI;
1053 std::vector<int>* vint;
1063 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
1064 cTo = pwc->
get((
char*)
"time",0,
'L',0);
1067 for(k=0; k<
K; k++) {
1068 id = size_t(cid.
data[k]+0.1);
1069 if(pwc->
sCuts[
id-1] != -2)
continue;
1070 vint = &(pwc->
cList[
id-1]);
1074 pI = wdmMRA.getXTalk(pwc,
id);
1080 tsize = pix->
tdAmp[0].size();
1081 if(!tsize || tsize&1) {
1082 cout<<
"network::subNetCut() error: wrong pixel TD data\n";
1086 V4 = V + (V%4 ? 4 - V%4 : 0);
1090 std::vector<wavearray<float> > vtd;
1091 std::vector<wavearray<float> > vTD;
1092 std::vector<wavearray<float> > eTD;
1111 __m128* _Fp = (__m128*) Fp.
data;
1112 __m128* _Fx = (__m128*) Fx.
data;
1113 __m128* _am = (__m128*) am.
data;
1114 __m128* _AM = (__m128*) AM.
data;
1115 __m128* _xi = (__m128*) xi.
data;
1116 __m128* _XI = (__m128*) XI.
data;
1117 __m128* _fp = (__m128*) fp.
data;
1118 __m128* _fx = (__m128*) fx.
data;
1119 __m128* _nr = (__m128*) nr.
data;
1120 __m128* _ww = (__m128*) ww.
data;
1121 __m128* _WW = (__m128*) WW.
data;
1122 __m128* _bb = (__m128*) bb.
data;
1123 __m128* _BB = (__m128*) BB.
data;
1125 for(i=0; i<NIFO; i++) {
1131 for(i=0; i<
NIFO; i++) {
1132 pa[
i] = vtd[
i].data + (tsize/2)*V4;
1133 pA[
i] = vTD[
i].data + (tsize/2)*V4;
1134 pe[
i] = eTD[
i].data + (tsize/2)*V4;
1137 this->a_00.resize(NIFO*V4); this->a_00=0.;
1138 this->a_90.resize(NIFO*V4); this->a_90=0.;
1139 this->rNRG.resize(V4); this->rNRG=0.;
1140 this->pNRG.resize(V4); this->pNRG=0.;
1142 __m128* _aa = (__m128*) this->a_00.data;
1143 __m128* _AA = (__m128*) this->a_90.data;
1145 this->pList.clear();
1146 for(j=0; j<V; j++) {
1148 pList.push_back(pix);
1151 for(i=0; i<
nIFO; i++) {
1156 for(i=0; i<
nIFO; i++) {
1157 nr.
data[j*NIFO+
i]=(float)
xx[i]/sqrt(rms);
1158 for(l=0; l<tsize; l++) {
1160 AA = pix->
tdAmp[
i].data[l+tsize];
1161 vtd[
i].data[l*V4+
j] = aa;
1162 vTD[
i].data[l*V4+
j] = AA;
1163 eTD[
i].data[l*V4+
j] = aa*aa+AA*AA;
1177 bool skymaskcc = (skyMaskCC.size()==Lsky);
1179 stat=Lm=Em=Am=EE=0.; lm=Vm= -1;
1183 for(l=lb; l<=le; l++) {
1184 if(!mm[l] || l<0)
continue;
1191 if (!skyMaskCC.data[lc])
continue;
1197 __m128 _E_o = _mm_setzero_ps();
1198 __m128 _E_n = _mm_setzero_ps();
1199 __m128 _E_s = _mm_setzero_ps();
1200 __m128 _M_m = _mm_setzero_ps();
1201 __m128* _rE = (__m128*) rNRG.data;
1202 __m128* _pE = (__m128*) pNRG.data;
1204 for(j=0; j<V4; j+=4) {
1206 _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_rE,_En));
1207 _M_m = _mm_add_ps(_M_m,_msk);
1208 *_pE = _mm_mul_ps(*_rE,_msk);
1209 _E_o = _mm_add_ps(_E_o,*_pE);
1211 _E_s = _mm_add_ps(_E_s,*_pE);
1212 _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_pE++,_Es));
1213 _E_n = _mm_add_ps(_E_n,_mm_mul_ps(*_rE++,_msk));
1216 _mm_storeu_ps(vvv,_E_n);
1217 Ln = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1218 _mm_storeu_ps(vvv,_E_o);
1219 Eo = vvv[0]+vvv[1]+vvv[2]+vvv[3]+0.01;
1220 _mm_storeu_ps(vvv,_E_s);
1221 Ls = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1222 _mm_storeu_ps(vvv,_M_m);
1223 m = 2*(vvv[0]+vvv[1]+vvv[2]+vvv[3])+0.01;
1226 if((aa-m)/(aa+
m)<subcut)
continue;
1228 pnt_(v00, pa, ml, (
int)l, (
int)V4);
1229 pnt_(v90, pA, ml, (
int)l, (
int)V4);
1230 float* pfp = fp.
data;
1231 float* pfx = fx.
data;
1232 float* p00 = this->a_00.data;
1233 float* p90 = this->a_90.data;
1236 for(j=0; j<V; j++) {
1238 cpp_(p00,v00); cpp_(p90,v90);
1239 cpf_(pfp,FP,l); cpf_(pfx,FX,l);
1244 if(rNRG.data[j]>En) m++;
1247 __m128* _pp = (__m128*) am.
data;
1248 __m128* _PP = (__m128*) AM.
data;
1252 _pp = (__m128*) xi.
data;
1253 _PP = (__m128*) XI.
data;
1257 for(j=0; j<V; j++) {
1272 Ls += ee-em; Eo += ee;
1273 if(ee-em>Es) Ln += ee;
1277 size_t m4 = m + (m%4 ? 4 - m%4 : 0);
1278 _E_n = _mm_setzero_ps();
1280 for(j=0; j<m4; j+=4) {
1284 _E_n = _mm_add_ps(_E_n,_E_s);
1286 _mm_storeu_ps(vvv,_E_n);
1288 Lo = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1289 AA = aa/(
fabs(aa)+
fabs(Eo-Lo)+2*m*(Eo-Ln)/Eo);
1291 em =
fabs(Eo-Lo)+2*
m;
1295 if(AA>stat && !mra) {
1296 stat=AA; Lm=Lo; Em=Eo; Am=aa; lm=
l; Vm=
m; suball=ee; EE=em;
1300 if(!mra && lm>=0) {mra=
true; le=lb=lm;
goto skyloop;}
1302 pwc->
sCuts[
id-1] = -1;
1303 pwc->
cData[
id-1].likenet = Lm;
1304 pwc->
cData[
id-1].energy = Em;
1305 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
1306 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
1307 pwc->
cData[
id-1].skyIndex = lm;
1311 submra = Ls*Eo/(Eo-Ls);
1312 submra/=
fabs(submra)+
fabs(Eo-Lo)+2*(m+6);
1314 pwc->
p_Ind[
id-1].push_back(lm);
1315 for(j=0; j<vint->size(); j++) {
1317 pix->
theta = nLikelihood.getTheta(lm);
1318 pix->
phi = nLikelihood.getPhi(lm);
1326 rHo = sqrt(Lo*Lo/(Eo+2*m)/2);
1329 if(hist && rHo>this->
netRHO)
1330 for(j=0;j<vint->size();j++) hist->Fill(suball,submra);
1332 if(fmin(suball,submra)>subnet && rHo>this->
netRHO) {
1333 count += vint->size();
1335 printf(
"lag|id %3d|%3d rho=%5.2f To=%5.1f stat: %5.3f|%5.3f|%5.3f ",
1336 int(lag),
int(
id),rHo,To,suball,submra,stat);
1337 printf(
"E: %6.1f|%6.1f L: %6.1f|%6.1f|%6.1f pix: %4d|%4d|%3d|%2d \n",
1338 Em,Eo,Lm,Lo,Ls,
int(vint->size()),
int(V),Vm,
int(m));
1341 else pwc->
sCuts[
id-1]=1;
1346 for(j=0; j<V; j++) {
1381 if(!this->wc_List[lag].
size())
return 0;
1387 bool cirwave = mode==
'g' || mode==
'G' || mode==
'c' || mode==
'C';
1388 bool linwave = mode==
'l' || mode==
'L' || mode==
's' || mode==
'S';
1389 bool iotwave = mode==
'i' || mode==
'l' || mode==
'e' || mode==
'c' ||
1390 mode==
'I' || mode==
'L' || mode==
'E' || mode==
'C';
1391 bool psiwave = mode==
'l' || mode==
'e' || mode==
'p' ||
1392 mode==
'L' || mode==
'E' || mode==
'P';
1393 bool mureana = mode==
'i' || mode==
'e' || mode==
'c' ||
1394 mode==
'r' || mode==
'p' || mode==
'b' ||
1395 mode==
'l' || mode==
's' || mode==
'g';
1396 bool rndwave = mode==
'r' || mode==
'R';
1398 bool prior = this->gamma<0 ?
true :
false;
1399 bool m_chirp = this->
optim ?
false : mureana;
1401 if(!this->
optim) mureana =
true;
1403 size_t nIFO = this->ifoList.size();
1404 size_t ID = abs(iID);
1407 cout<<
"network::likelihood2G(): invalid network.\n";
1413 float gama =
fabs(this->gamma);
1415 if(gama<=0) gama = 1.e-24;
1416 if(gama>=1) gama = 0.999999;
1424 static const __m128 _D0 = _mm_set1_ps(deta<0.5?1-deta:0.5);
1425 static const __m128 _D9 = _mm_set1_ps(deta<0.5?1:1.5-deta);
1427 static const __m128 _oo = _mm_set1_ps(1.
e-16);
1428 static const __m128 _sm = _mm_set1_ps(-0.
f);
1429 static const __m128 _En = _mm_set1_ps(En);
1430 static const __m128 _rG = _mm_set1_ps(-1./
log(gama));
1431 static const __m128 _kG = _mm_set1_ps(gama);
1432 static const __m128 _PW = _mm_set1_ps(psiwave?0:1);
1433 static const __m128 _01 = _mm_set1_ps(0.1);
1434 static const __m128 _05 = _mm_set1_ps(0.5);
1435 static const __m128 _09 = _mm_set1_ps(0.9);
1436 static const __m128 _1 = _mm_set1_ps(1.0+1.
e-16);
1437 static const __m128 _2 = _mm_set1_ps(2.0);
1438 static const __m128 _4 = _mm_set1_ps(4.0);
1443 float NRG,Lm,Em,Lo,Eo,No,Nm,cc,Cm,Co,Do,To,Fo,Ln,Ns;
1444 float STAT,CHR,aa,AA,ee,em,EE,
ff,FF,Lr,Cr,
ss,Ls,Nc,gg;
1445 double eLp, s2p, c2p;
1448 size_t i,
j,
k,
l,
m,Vm,lm,V,V4,V44,
id,
K,
M;
1450 short* mm = this->skyMask.data;
1451 bool skymaskcc = (skyMaskCC.size()==
L);
1465 for(i=0; i<
NIFO; i++) {
1467 ml[
i] = getifo(i)->index.data;
1468 FP[
i] = getifo(i)->fp.data;
1469 FX[
i] = getifo(i)->fx.data;
1472 ml[
i] = getifo(0)->index.data;
1473 FP[
i] = getifo(0)->fp.data;
1474 FX[
i] = getifo(0)->fx.data;
1479 std::vector<int> pI;
1480 std::vector<int> pJ;
1484 std::vector<int>* vint;
1485 std::vector<int>* vtof;
1486 std::vector<int> pRate;
1492 std::map<int,float> vLr;
1498 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
1499 cTo = pwc->
get((
char*)
"time",0,
'L',0);
1502 for(k=0; k<
K; k++) {
1503 id = size_t(cid.
data[k]+0.1);
1505 if(pwc->
sCuts[
id-1] != -2)
continue;
1507 vint = &(pwc->
cList[
id-1]);
1508 vtof = &(pwc->
nTofF[
id-1]);
1512 pI = wdmMRA.getXTalk(pwc,
id);
1517 this->nSensitivity = 0.;
1518 this->nAlignment = 0.;
1519 this->nNetIndex = 0.;
1520 this->nDisbalance = 0.;
1521 this->nLikelihood = 0.;
1522 this->nNullEnergy = 0.;
1523 this->nCorrEnergy = 0.;
1524 this->nCorrelation = 0.;
1525 this->nSkyStat = 0.;
1526 this->nEllipticity = 0.;
1527 this->nPolarisation = 0.;
1528 this->nProbability = 0.;
1530 this->nAntenaPrior = 0.;
1533 tsize = pix->
tdAmp[0].size();
1534 if(!tsize || tsize&1) {
1535 cout<<
"network::likelihood2G() error: wrong pixel TD data\n";
1541 if(!(V=pI.size()))
continue;
1542 V4 = V + (V%4 ? 4 - V%4 : 0);
1545 for(j=0; j<V4; j++) pJ.push_back(0);
1549 std::vector<wavearray<float> > vtd;
1550 std::vector<wavearray<float> > vTD;
1551 std::vector<wavearray<float> > eTD;
1574 this->p00_POL[
i].
resize(V4); this->p00_POL[
i]=0.;
1575 this->p90_POL[
i].resize(V4); this->p90_POL[
i]=0.;
1576 this->r00_POL[
i].resize(V4); this->r00_POL[
i]=0.;
1577 this->r90_POL[
i].resize(V4); this->r90_POL[
i]=0.;
1580 __m128* _Fp = (__m128*) Fp.
data;
1581 __m128* _Fx = (__m128*) Fx.
data;
1582 __m128* _am = (__m128*) am.
data;
1583 __m128* _AM = (__m128*) AM.
data;
1584 __m128* _xi = (__m128*) xi.
data;
1585 __m128* _XI = (__m128*) XI.
data;
1586 __m128* _xp = (__m128*) xp.
data;
1587 __m128* _XP = (__m128*) XP.
data;
1588 __m128* _ww = (__m128*) ww.
data;
1589 __m128* _WW = (__m128*) WW.
data;
1590 __m128* _bb = (__m128*) bb.
data;
1591 __m128* _BB = (__m128*) BB.
data;
1592 __m128* _fp = (__m128*) fp.
data;
1593 __m128* _fx = (__m128*) fx.
data;
1594 __m128* _nr = (__m128*) nr.
data;
1595 __m128* _ep = (__m128*) ep.
data;
1596 __m128* _ex = (__m128*) ex.
data;
1598 __m128* _fp4 = _fp+V4*f_;
1599 __m128* _fx4 = _fx+V4*f_;
1600 __m128* _uu4 = _am+V4*f_;
1601 __m128* _vv4 = _AM+V4*f_;
1602 __m128* _bb4 = _bb+V4*f_;
1603 __m128* _BB4 = _BB+V4*f_;
1605 for(i=0; i<NIFO; i++) {
1611 for(i=0; i<
NIFO; i++) {
1612 pa[
i] = vtd[
i].data + (tsize/2)*V4;
1613 pA[
i] = vTD[
i].data + (tsize/2)*V4;
1614 pe[
i] = eTD[
i].data + (tsize/2)*V4;
1632 this->a_00.
resize(NIFO*V44); this->a_00=0.;
1633 this->a_90.resize(NIFO*V44); this->a_90=0.;
1634 this->rNRG.resize(V4); this->rNRG=0.;
1635 this->pNRG.resize(V4); this->pNRG=0.;
1637 __m128* _aa = (__m128*) this->a_00.data;
1638 __m128* _AA = (__m128*) this->a_90.data;
1641 this->pList.clear(); pRate.clear();
1642 for(j=0; j<V; j++) {
1644 this->pList.push_back(pix);
1645 pRate.push_back(
int(pix->
rate+0.5));
1647 if(vLr.find(pRate[j]) == vLr.end())
1651 for(i=0; i<
nIFO; i++) {
1657 for(i=0; i<
nIFO; i++) {
1658 nr.
data[j*NIFO+
i]=(float)
xx[i]/rms;
1659 for(l=0; l<tsize; l++) {
1661 AA = pix->
tdAmp[
i].data[l+tsize];
1662 vtd[
i].data[l*V4+
j] = aa;
1663 vTD[
i].data[l*V4+
j] = AA;
1664 eTD[
i].data[l*V4+
j] = aa*aa+AA*AA;
1673 STAT=0.; lm=0; Vm=0;
1674 double skystat = 0.;
1682 for(l=lb; l<=le; l++) {
1683 if(!mra) skyProb.data[
l] = 0.;
1684 if(!mm[l])
continue;
1691 if (!skyMaskCC.data[lc])
continue;
1694 pnt_(v00, pa, ml, (
int)l, (
int)V4);
1695 pnt_(v90, pA, ml, (
int)l, (
int)V4);
1696 float* pfp = fp.
data;
1697 float* pfx = fx.
data;
1698 float* p00 = this->a_00.data;
1699 float* p90 = this->a_90.data;
1703 for(j=0; j<V; j++) {
1704 cpp_(p00,v00); cpp_(p90,v90);
1705 cpf_(pfp,FP,l); cpf_(pfx,FX,l);
1706 if(!this->
optim || !mra)
continue;
1707 if(vLr[pRate[j]] <= mxLr)
continue;
1708 mxLr = vLr[pRate[
j]];
1713 for(j=0; j<V; j++) {
1718 if(optR && optR!=pRate[j]) {
1726 if(ee>En) m++;
else ee=0.;
1728 pNRG.data[
j] = rNRG.data[
j];
1731 __m128* _pp = (__m128*) am.
data;
1732 __m128* _PP = (__m128*) AM.
data;
1734 if(mra && mureana) {
1736 _pp = (__m128*) xi.
data;
1737 _PP = (__m128*) XI.
data;
1741 for(j=0; j<V; j++) {
1758 size_t m4 = m + (m%4 ? 4 - m%4 : 0);
1759 __m128 _ll,_LL,_ec,_EC,_ee,_EE,_NI,_s2,_c2,_AX,_NN,_FF,_QQ,_ie,_gg;
1760 __m128 _en,_EN,_ed,_ED,_cc,_ss,_ni,_si,_co,_ax,_nn,_ff,_mm,_IE,_GG;
1762 __m128* _siP = (__m128*) siPHS.
data;
1763 __m128* _coP = (__m128*) coPHS.
data;
1764 __m128* _siO = (__m128*) siORT.
data;
1765 __m128* _coO = (__m128*) coORT.
data;
1766 __m128* _siD = (__m128*) siDPF.
data;
1767 __m128* _coD = (__m128*) coDPF.
data;
1768 __m128* _nrg = (__m128*) q2Q2.
data;
1769 __m128* _chr = (__m128*) chir.
data;
1772 _pp = (__m128*) xi.
data;
1773 _PP = (__m128*) XI.
data;
1778 Lo = Ln = Co = Eo = 0.;
1779 for(j=0; j<m4; j+=4) {
1782 __m128* _pbb = _bb+jf;
1783 __m128* _pBB = _BB+jf;
1784 __m128* _pxi = _pp+jf;
1785 __m128* _pXI = _PP+jf;
1786 __m128* _pxp = _xp+jf;
1787 __m128* _pXP = _XP+jf;
1788 __m128* _pww = _ww+jf;
1789 __m128* _pWW = _WW+jf;
1790 __m128* _pfp = _fp+jf;
1791 __m128* _pfx = _fx+jf;
1792 __m128* _pFp = _Fp+jf;
1793 __m128* _pFx = _Fx+jf;
1806 _coO++; _siO++; _coD++; _siD++;
1809 if(le==lb && (optR==0)) {
1810 _sse_pol4_ps(_pfp, _pfx, _pxp, p00_POL[0].data+j, p00_POL[1].data+j);
1811 _sse_pol4_ps(_pfp, _pfx, _pXP, p90_POL[0].data+j, p90_POL[1].data+j);
1812 _sse_pol4_ps(_pfp, _pfx, _pxi, r00_POL[0].data+j, r00_POL[1].data+j);
1813 _sse_pol4_ps(_pfp, _pfx, _pXI, r90_POL[0].data+j, r90_POL[1].data+j);
1824 _cc = _mm_and_ps(_mm_cmpge_ps(_ec,_05),_1);
1825 _ss = _mm_and_ps(_mm_cmpgt_ps(_EC,_05),_cc);
1830 _nn = _mm_add_ps(_mm_add_ps(_ee,_EE),_2);
1831 _cc = _mm_add_ps(_ec,_mm_sub_ps(_nn,_ll));
1832 _cc = _mm_div_ps(_ec,_mm_add_ps(_cc,_mm));
1834 _mm_storeu_ps(vvv,_mm_add_ps(_ee,_EE));
1835 Eo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1836 _mm_storeu_ps(vvv,_ll);
1837 Lo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1838 _mm_storeu_ps(vvv,_mm_mul_ps(_ll,_cc));
1839 Ln += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1840 _mm_storeu_ps(vvv,_ec);
1841 Co += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1842 _mm_storeu_ps(vvv,_mm_mul_ps(_ec,_ll));
1843 if(le==lb && !mra) {
1844 vLr[pRate[pJ[j+0]]] += vvv[0];
1845 vLr[pRate[pJ[j+1]]] += vvv[1];
1846 vLr[pRate[pJ[j+2]]] += vvv[2];
1847 vLr[pRate[pJ[j+3]]] += vvv[3];
1850 Ln = Eo>0 ? Ln/Eo : 0;
1851 if(Ln<this->
netCC)
continue;
1852 _IE = _mm_set1_ps(1-Co/Lo);
1857 __m128* _xx = (__m128*) xxx.
data;
1858 __m128* _yy = (__m128*) yyy.
data;
1859 __m128* _zz = (__m128*) zzz.
data;
1860 __m128* _rr = (__m128*) rrr.
data;
1862 for(j=0; j<m4; j+=4) {
1864 __m128* _pfp = _fp+jf;
1865 __m128* _pfx = _fx+jf;
1868 _ee = _mm_add_ps(_mm_sqrt_ps(*_xx++),_oo);
1869 _EE = _mm_add_ps(_mm_sqrt_ps(*_yy++),_oo);
1874 _xx = (__m128*) xxx.
data;
1875 _yy = (__m128*) yyy.
data;
1876 _zz = (__m128*) zzz.
data;
1877 _rr = (__m128*) rrr.
data;
1878 _siP = (__m128*) siPHS.
data;
1879 _coP = (__m128*) coPHS.
data;
1880 _siD = (__m128*) siDPF.
data;
1881 _coD = (__m128*) coDPF.
data;
1882 _nrg = (__m128*) q2Q2.
data;
1884 __m128 linw = linwave ? _oo : _1;
1885 __m128 cirw = cirwave ? _oo : _1;
1886 __m128 _ch = _mm_setzero_ps();
1887 __m128 _eqQ = _mm_setzero_ps();
1888 __m128 _qQ2 = _mm_setzero_ps();
1890 _c2 = _mm_setzero_ps();
1891 _s2 = _mm_setzero_ps();
1892 _cc = _mm_setzero_ps();
1893 _ss = _mm_setzero_ps();
1895 for(j=0; j<m4; j+=4) {
1897 __m128* _pfp = _fp+jf;
1898 __m128* _pfx = _fx+jf;
1899 __m128* _pep = _ep+jf;
1900 __m128* _pex = _ex+jf;
1901 __m128* _pxp = _xp+jf;
1902 __m128* _pXP = _XP+jf;
1903 __m128* _pxi = _pp+jf;
1904 __m128* _pXI = _PP+jf;
1906 _ee = _mm_div_ps(_1,_mm_add_ps(*_xx,_oo));
1912 _nn = _mm_div_ps(_1,_mm_add_ps(_NN,_oo));
1913 *_siP = _mm_mul_ps(_si,_nn);
1914 *_coP = _mm_mul_ps(_co,_nn);
1916 if(le==lb && (optR==0)) {
1917 __m128 _snn = _mm_sqrt_ps(_nn);
1920 _sse_pol4_ps(_ep+jf, _ex+jf, _uu4, r00_POL[0].data+j, r00_POL[1].data+j);
1921 _sse_pol4_ps(_ep+jf, _ex+jf, _vv4, r90_POL[0].data+j, r90_POL[1].data+j);
1926 _ch = _mm_add_ps(_ch,_EN);
1927 _ll = _mm_and_ps(_mm_cmpgt_ps(_EN,_oo),_1);
1928 *_chr = _mm_sub_ps(_mm_mul_ps(_ll,_2),_1);
1934 _ff = _mm_div_ps(_mm_add_ps(*_xx,*_yy),_ni);
1937 _LL = _mm_mul_ps(_LL,_nn);
1938 _FF = _mm_mul_ps(*_yy,_ee);
1939 _gg = _mm_mul_ps(_01,_ff);
1940 _NN = _mm_mul_ps(_NN,*_xx);
1941 _ll = _mm_mul_ps(_NN,_mm_add_ps(_FF,_gg));
1942 _co = _mm_sub_ps(_ll,_LL);
1945 *_zz = _mm_add_ps(_mm_sqrt_ps(*_zz),_oo);
1946 _ll = _mm_mul_ps(_NN,_FF);
1947 _QQ = _mm_add_ps(_mm_add_ps(_ll,_LL),*_zz);
1948 _ff = _mm_sqrt_ps(_mm_mul_ps(_ff,_05));
1949 _FF = _mm_sqrt_ps(_FF);
1951 _ax = _mm_add_ps(_1,_mm_div_ps(_co,*_zz));
1952 _ax = _mm_sqrt_ps(_mm_mul_ps(_ax,_05));
1954 _qQ2 = _mm_add_ps(_qQ2,_QQ);
1955 _eqQ = _mm_add_ps(_eqQ,_AX);
1956 _AX = _mm_div_ps(_mm_mul_ps(_2,_AX),_QQ);
1958 _gg = _mm_mul_ps(_mm_sub_ps(_05,_ni),
1959 _mm_sub_ps(_1,_FF));
1960 _GG = _mm_sub_ps(_mm_sub_ps(_09,_GG),_gg);
1961 _GG = _mm_mul_ps(_mm_mul_ps(_IE,_rG),_GG);
1962 _GG = _mm_mul_ps(_mm_mul_ps(_IE,_4),_GG);
1964 _gg = _mm_mul_ps(_mm_sub_ps(_IE,_ff),_kG);
1965 _gg = _mm_mul_ps(_ni,_gg);
1966 _nn = _mm_and_ps(_mm_cmpgt_ps(_ff,_gg),_1);
1970 _mm = _mm_and_ps(_mm_cmpgt_ps(_ff,_GG),_nn);
1972 _gg = _mm_andnot_ps(_sm,_mm_sub_ps(_NI,_ni));
1973 _gg = _mm_sub_ps(_IE,_mm_mul_ps(_gg,_FF));
1974 _nn = _mm_and_ps(_mm_cmplt_ps(_gg,_D0),_mm);
1975 _mm = _mm_and_ps(_mm_cmplt_ps(_gg,_D9),_mm);
1976 _nn = _mm_mul_ps(_nn,cirw);
1977 _mm = _mm_mul_ps(_mm,linw);
1978 _EE = _mm_mul_ps(_en,_nn);
1982 _nn = _mm_mul_ps(_nn,_PW);
1983 *_rr = _mm_add_ps(_mm_add_ps(_mm,_nn),_1);
1984 _coP++;_siP++;_chr++;_xx++;_yy++;_zz++,_rr++;
1986 _cc = _mm_add_ps(_cc,_ee);
1987 _ss = _mm_add_ps(_ss,_EE);
1993 _c2 = _mm_add_ps(_c2,_ec);
1994 _s2 = _mm_sub_ps(_s2,_EC);
1998 _mm_storeu_ps(vvv,_c2);
1999 c2p = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2000 _mm_storeu_ps(vvv,_s2);
2001 s2p = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2002 gg = sqrt(c2p*c2p+s2p*s2p+1.
e-16);
2003 _si = _mm_set1_ps(s2p/gg);
2004 _co = _mm_set1_ps(c2p/gg);
2007 _zz = (__m128*) zzz.
data;
2008 _siD = (__m128*) siDPF.
data;
2009 _coD = (__m128*) coDPF.
data;
2011 _mm_storeu_ps(vvv,_cc);
2012 cc = (vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2013 _mm_storeu_ps(vvv,_ss);
2014 ss = (vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2015 gg = sqrt(cc*cc+ss*ss+1.
e-16);
2016 _si = _mm_set1_ps(ss/gg);
2017 _co = _mm_set1_ps(cc/gg);
2019 for(j=0; j<m4; j+=4) {
2021 __m128* _pxi = _pp+jf;
2022 __m128* _pXI = _PP+jf;
2040 _siO = (__m128*) siORT.
data;
2041 _coO = (__m128*) coORT.
data;
2042 _siP = (__m128*) siPHS.
data;
2043 _coP = (__m128*) coPHS.
data;
2044 _chr = (__m128*) chir.
data;
2046 _mm_storeu_ps(vvv,_eqQ);
2047 _mm_storeu_ps(uuu,_qQ2);
2048 eLp = 2.*(vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2049 eLp/= uuu[0]+uuu[1]+uuu[2]+uuu[3]+1.e-16;
2050 _mm_storeu_ps(vvv,_ch);
2051 CHR = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2052 ff = CHR>0. ? 1. : -1.;
2053 _ch = _mm_set1_ps(ff);
2054 _gg = rndwave ? _oo : _mm_set1_ps(0.5);
2056 for(j=0; j<m4; j+=4) {
2058 __m128* _pbb = _bb+jf;
2059 __m128* _pBB = _BB+jf;
2060 __m128* _pxp = _xp+jf;
2061 __m128* _pXP = _XP+jf;
2062 __m128* _pxi = _pp+jf;
2063 __m128* _pXI = _PP+jf;
2065 _ee = _mm_sub_ps(_mm_mul_ps(_ch,*_chr),_1);
2077 _coO++; _siO++; _coP++; _siP++; _chr++;
2084 _rr = (__m128*) rrr.
data;
2086 Lo = Co = Eo = Lr = Cr = Do = 0.;
2087 for(j=0; j<m4; j+=4) {
2090 __m128* _pbb = _bb+jf;
2091 __m128* _pBB = _BB+jf;
2092 __m128* _pxi = _pp+jf;
2093 __m128* _pXI = _PP+jf;
2094 __m128* _pww = _ww+jf;
2095 __m128* _pWW = _WW+jf;
2107 _ll = _mm_add_ps(_ll,_LL);
2108 _ec = _mm_add_ps(_ec,_EC);
2109 _ed = _mm_add_ps(_ed,_ED);
2110 _ee = _mm_add_ps(_ee,_EE);
2112 _mm_storeu_ps(vvv,_ee);
2113 Eo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2114 _mm_storeu_ps(vvv,_ec);
2115 Co += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2116 _mm_storeu_ps(vvv,_ed);
2117 Do += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2118 _mm_storeu_ps(vvv,_ll);
2119 Lo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2123 _en = _mm_andnot_ps(_sm,_mm_sub_ps(_ee,_ll));
2124 _en = _mm_add_ps(_en,*_rr); _rr++;
2125 _EC = _mm_andnot_ps(_sm,_ec);
2126 _cc = _mm_add_ps(_mm_add_ps(_EC,_ed),_en);
2127 _cc = _mm_div_ps(_mm_sub_ps(_ec,_ed),_cc);
2129 _mm_storeu_ps(vvv,_mm_mul_ps(_ll,_cc));
2130 Lr += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2131 _mm_storeu_ps(vvv,_mm_mul_ps(_ec,_cc));
2132 Cr += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2142 aa = Eo>0. ? Lo/Eo : 0.;
2143 AA = Eo>0. ? Lr/Eo : 0.;
2144 if(!mra) skyProb.data[
l] = this->
delta<0 ? aa : AA;
2148 ff = FF = Et = Nm = 0.;
2149 for(j=0; j<
m; j++) {
2160 Nm = Et>0.&&Nm>0 ? Et/Nm : 0.;
2161 ff = Eo>0. ? 2*ff/Eo : 0.;
2162 FF = Eo>0. ? 2*FF/Eo : 0.;
2165 if(ID==
id && !mra) {
2166 Eo += 0.001; Cr += 0.001;
2167 this->nAntenaPrior.set(l, sqrt(ff+FF));
2168 this->nSensitivity.set(l, sqrt(ff+FF));
2169 this->nAlignment.set(l, sqrt(FF/ff));
2170 this->nLikelihood.set(l, Lo/Eo);
2171 this->nNullEnergy.set(l, (Eo-Lo)/Eo);
2172 this->nCorrEnergy.set(l, Cr/Eo);
2173 this->nCorrelation.set(l, Ln);
2174 this->nSkyStat.set(l, AA);
2175 this->nProbability.set(l, skyProb.data[l]);
2176 this->nDisbalance.set(l, 2*Do/Eo);
2177 this->nEllipticity.set(l, eLp);
2178 this->nPolarisation.set(l, -atan2(s2p,c2p)*180./
PI/4.);
2179 this->nNetIndex.set(l, Nm);
2182 if(prior && !mra && ID!=
id) {
2184 for(j=0; j<
m; j++) {
2191 ff = Eo>0. ? 2*ff/Eo : 0.;
2192 FF = Eo>0. ? 2*FF/Eo : 0.;
2193 this->nAntenaPrior.set(l, sqrt(ff+FF));
2196 if(AA>STAT && !mra) {STAT=AA; lm=
l; Vm=
m;}
2199 if(STAT==0. || (mra && AA<=0.)) {
2200 pwc->
sCuts[
id-1]=1; count=0;
2201 pwc->
clean(
id);
continue;
2204 if(le-lb) {lb=le=lm;
goto optsky;}
2206 double Em_all,Ln_all,Ls_all,Lr_all;
2207 double Eo_all,Lo_all,Co_all,Do_all,cc_all;
2211 Em_all=Ls_all=Ln_all = 0;
2217 pwc->
cData[
id-1].skySize =
m;
2218 pwc->
cData[
id-1].likesky = Lo;
2219 vint = &(pwc->
cList[
id-1]);
2220 for(j=0; j<vint->size(); j++) {
2226 for(j=0; j<Vm; j++) {
2234 if(ee-em>Es) Ln_all += ee;
2235 GNoise += rrr.
data[
j];
2238 pwc->
cData[
id-1].skyStat = Lr/Eo;
2241 Ns = Eo_all-Lo_all+Do+GNoise;
2242 gg = Ls_all*Ln_all/Em_all;
2243 pwc->
cData[
id-1].SUBNET = gg/(
fabs(gg)+Ns);
2245 mra=
true;
goto optsky;
2248 if(AA<this->
netCC || !m) {
2249 pwc->
sCuts[
id-1]=1; count=0;
2250 pwc->
clean(
id);
continue;
2259 double Em_mra,Ln_mra,Ls_mra;
2260 double Eo_mra,Lo_mra,Co_mra;
2262 M=
m; m=0; GNoise=0.;
2263 Em_mra=Ln_mra=Ls_mra = 0;
2264 for(j=0; j<
M; j++) {
2268 float* pco = coORT.
data+
j;
2269 __m128* _pxi = _xi+jf;
2270 __m128* _pXI = _XI+jf;
2275 if(ee-em>Es) Ln_mra += ee;
2276 GNoise += rrr.
data[
j];
2283 for(i=0; i<
nIFO; i++) {
2292 pwc->
sCuts[
id-1]=1; count=0;
2293 pwc->
clean(
id);
continue;
2296 Em=Eo; Lm=Lo; Do*=2;
2297 Eo_mra=Eo; Lo_mra=Lo; Co_mra=Co;
2299 pwc->
cData[
id-1].netcc = Lr/Eo;
2300 Nc = Eo-Lo+Do/2+GNoise;
2301 gg = Ls_mra*Ln_mra/Em_mra;
2302 pwc->
cData[
id-1].subnet = gg/(
fabs(gg)+Nc);
2303 pwc->
cData[
id-1].skycc = Co/(
fabs(Co)+Nc);
2314 NETX (vtof->push_back(ml[0][lm]); ,
2315 vtof->push_back(ml[1][lm]); ,
2316 vtof->push_back(ml[2][lm]); ,
2317 vtof->push_back(ml[3][lm]); ,
2318 vtof->push_back(ml[4][lm]); ,
2319 vtof->push_back(ml[5][lm]); ,
2320 vtof->push_back(ml[6][lm]); ,
2321 vtof->push_back(ml[7][lm]); )
2324 if((wfsave)||(mdcListSize() && !lag)) {
2325 int m0d = mureana ? 0 : 1;
2326 if(this->getMRAwave(
id,lag,
'S',m0d,
true)) {
2328 for(i=0; i<
nIFO; i++) {
2329 pd = this->getifo(i);
2330 pd->
RWFID.push_back(
id);
2334 pd->
RWFP.push_back(wf);
2337 if(this->getMRAwave(
id,lag,
's',m0d,
true)) {
2339 for(i=0; i<
nIFO; i++) {
2340 pd = this->getifo(i);
2341 pd->
RWFID.push_back(-
id);
2345 pd->
RWFP.push_back(wf);
2350 Lo = Eo = To = Fo = No = 0.;
2351 for(i=0; i<
nIFO; i++) {
2356 int two = mureana ? 1 : 2;
2357 int m0d = mureana ? 0 : -1;
2359 this->getMRAwave(
id,lag,
'W',m0d);
2360 this->getMRAwave(
id,lag,
'S',m0d);
2361 for(i=0; i<
nIFO; i++) {
2365 float sSNR = d->
get_SS()/two;
2366 float xSNR = d->
get_XS()/two;
2367 float null = d->
get_NN()/two;
2368 float enrg = d->
get_XX()/two;
2389 pwc->
cData[
id-1].likenet = Lo;
2390 pwc->
cData[
id-1].energy = Eo;
2391 pwc->
cData[
id-1].enrgsky = Eo_all;
2392 pwc->
cData[
id-1].netecor = Co;
2393 pwc->
cData[
id-1].netnull = No+GNoise/2;
2394 pwc->
cData[
id-1].netED = Do;
2395 pwc->
cData[
id-1].netRHO = sqrt(Co*cc_all/(nIFO-1.));
2396 pwc->
cData[
id-1].netrho = sqrt(Cr/(nIFO-1.));
2397 pwc->
cData[
id-1].cTime = To;
2398 pwc->
cData[
id-1].cFreq = Fo;
2399 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
2400 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
2401 pwc->
cData[
id-1].gNET = sqrt(ff+FF);
2402 pwc->
cData[
id-1].aNET = sqrt(FF/ff);
2403 pwc->
cData[
id-1].iNET = Nm;
2404 pwc->
cData[
id-1].iota = Ns;
2405 pwc->
cData[
id-1].psi = -atan2(s2p,c2p)*180./
PI/4.;
2406 pwc->
cData[
id-1].ellipticity = eLp;
2410 if(sqrt(Co/(nIFO-1.))<this->
netRHO || pwc->
cData[
id-1].skycc<this->netCC) {
2411 pwc->
sCuts[
id-1]=1; count=0;
2412 pwc->
clean(
id);
continue;
2415 cc = pwc->
cData[
id-1].skycc;
2417 printf(
"id|lm %3d|%6d rho=%4.2f cc: %5.3f|%5.3f|%5.3f|%5.3f \n",
2418 int(
id),
int(lm),sqrt(Co/(nIFO-1)),STAT,cc,pwc->
cData[
id-1].netcc,AA);
2419 printf(
" (t,p)=(%4.1f|%4.1f) T|F: %6.3f|%4.1f L: %5.1f|%5.1f|%5.1f E: %5.1f|%5.1f|%5.1f \n",
2420 nLikelihood.getTheta(l),nLikelihood.getPhi(l),To,Fo,Lo,Lo_mra,Lo_all,Eo,Em,Eo_all);
2421 printf(
" D|N: %4.1f|%4.1f|%4.1f Vm|m=%3d|%3d subnet=%4.3f|%4.3f \n",
2422 Do,No,Nc,
int(Vm),
int(M),pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
2423 hist->Fill(pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
2429 pwc->
p_Ind[
id-1].push_back(m);
2430 double T = To+pwc->
start;
2431 std::vector<float> sArea;
2432 pwc->
sArea.push_back(sArea);
2433 pwc->
p_Map.push_back(sArea);
2437 double rMs = this->
delta<0 ? 0 : 2;
2438 if(iID<=0 || ID==
id) getSkyArea(
id,lag,T,rMs);
2442 pwc->
cData[
id-1].mchirp = 0;
2443 pwc->
cData[
id-1].mchirperr = 0;
2444 pwc->
cData[
id-1].tmrgr = 0;
2445 pwc->
cData[
id-1].tmrgrerr = 0;
2446 pwc->
cData[
id-1].chi2chirp = 0;
2450 cc = Co_all/(
fabs(Co_all)+ee);
2451 printf(
"mchirp : %d %g %.2e %.3f %.3f %.3f %.3f \n\n",
2452 int(
id),cc,pwc->
cData[
id-1].mchirp,
2453 pwc->
cData[
id-1].mchirperr, pwc->
cData[
id-1].tmrgr,
2454 pwc->
cData[
id-1].tmrgrerr, pwc->
cData[
id-1].chi2chirp);
2457 if(ID==
id && !
EFEC) {
2458 this->nSensitivity.gps =
T;
2459 this->nAlignment.gps =
T;
2460 this->nDisbalance.gps =
T;
2461 this->nLikelihood.gps =
T;
2462 this->nNullEnergy.gps =
T;
2463 this->nCorrEnergy.gps =
T;
2464 this->nCorrelation.gps =
T;
2465 this->nSkyStat.gps =
T;
2466 this->nEllipticity.gps =
T;
2467 this->nPolarisation.gps=
T;
2468 this->nNetIndex.gps =
T;
2471 pwc->
sCuts[
id-1] = -1;
2482 this->wfsave = value.
wfsave;
2483 this->MRA = value.
MRA;
2484 this->nRun = value.
nRun;
2485 this->nLag = value.
nLag;
2487 this->mIFO = value.
mIFO;
2488 this->Step = value.
Step;
2489 this->Edge = value.
Edge;
2491 this->aNET = value.
aNET;
2492 this->iNET = value.
iNET;
2493 this->eCOR = value.
eCOR;
2494 this->e2or = value.
e2or;
2496 this->norm = value.
norm;
2498 this->local = value.
local;
2502 this->gamma = value.
gamma;
2503 this->penalty = value.
penalty;
2506 this->pSigma = value.
pSigma;
2507 this->ifoList = value.
ifoList;
2510 this->NDM.clear(); this->NDM=value.
NDM;
2511 this->ifoList.clear(); this->ifoList=value.
ifoList;
2512 this->ifoName.clear(); this->ifoName=value.
ifoName;
2513 this->wc_List.clear(); this->wc_List=value.
wc_List;
2515 this->mdcList.clear(); this->mdcList=value.
mdcList;
2516 this->livTime.clear(); this->livTime=value.
livTime;
2517 this->mdcTime.clear(); this->mdcTime=value.
mdcTime;
2518 this->mdcType.clear(); this->mdcType=value.
mdcType;
2519 this->mdc__ID.clear(); this->mdc__ID=value.
mdc__ID;
2530 if(ifoList.size()==
NIFO) {
2531 cout <<
"network::add - Error : max number of detectors is " <<
NIFO << endl;
2537 this->ifoList.push_back(d);
2538 this->ifoName.push_back(d->
Name);
2542 for(i=0; i<
n; i++) {
2544 if(i<n-1) this->NDM[
i].push_back(0);
2545 else this->NDM.push_back(v);
2549 return ifoList.size();
2558 int N = ifoListSize();
2561 size_t nL = size_t(Edge*pw->
wrate()*
M);
2562 size_t nR = pw->
size() - nL - 1;
2564 for(
int i=1;
i<
N;
i++) w += getifo(
i)->TFmap;
2565 double amp,
avr, bbb, alp;
2568 for(
int i=nL;
i<nR;
i++) {
2569 amp = (double)w.
data[
i];
2570 if(amp>N*100) amp = N*100.;
2571 if(amp>0.001) {avr+=amp; bbb+=
log(amp); nn++;}
2574 alp =
log(avr)-bbb/nn;
2575 alp = (3-alp+sqrt((alp-3)*(alp-3)+24*alp))/12./alp;
2578 return avr*
iGamma(alp,bbb)/alp/2;
2587 int N = ifoListSize();
2590 size_t nL = size_t(Edge*pw->
wrate()*
M);
2591 size_t nR = pw->
size() - nL;
2593 for(
int i=1;
i<
N;
i++) w += getifo(
i)->TFmap;
2597 double v10 = w.
waveSplit(nL,nR,nR-
int(p10*fff*(nR-nL)));
2599 double med = w.
waveSplit(nL,nR,nR-
int(0.2*fff*(nR-nL)));
2601 while(p00<0.2) {p00 = 1-
Gamma(N*m,med); m+=0.01;}
2603 printf(
"\nm\tM\tbpp\t0.2(D)\t0.2(G)\t0.01(D)\t0.01(G)\tbpp(D)\tbpp(G)\tN*log(m)\tfff\n");
2604 printf(
"%g\t%d\t%g\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t\t%.3f\n\n",
2605 m,M,p,med,
iGamma(N*m,0.2),v10,
iGamma(N*m,p10),val,
iGamma(N*m,p),N*
log(m),fff);
2606 return (
iGamma(N*m,p)+val)*0.3+N*
log(m);
2616 size_t K = ifoListSize();
2618 if(getifo(0) && t>0. && K) {
2619 I = getifo(0)->TFmap.maxLayer()+1;
2620 n = 2*t*getifo(0)->TFmap.rate()/I + 1;
2621 if(!getifo(0)->TFmap.size()) n = 1.;
2623 else if(t < 0.) n = -
t;
2624 return sqrt(
iGamma1G(K/2.,1.-(p/2.)/pow(n,K/2.))/K);
2629 int iTYPE = this->MRA ? 0 : 1;
2642 cout<<
"wrong size "<<cid.
size()<<
" "<<rho.
size()<<endl;
2644 for(
size_t i=0;
i<cid.
size();
i++){
2645 printf(
"%2d %5.0f vol=%4.0f size=%4.0f like=%5.1e rho=%5.1f ",
2646 int(n),cid[
i],vol[i],siz[i],lik[i],rho[i]);
2647 printf(
"sub=%3.2f rate=%4.0f time=%8.3f To=%8.3f freq=%5.0f\n",
2648 sub[i],rat[i],tim[i],T_o[i],frq[i]);
2660 skymap temp(sms,t1,t2,p1,p2);
2661 size_t m = temp.
size();
2662 size_t n = this->ifoList.size();
2664 nSensitivity = temp;
2666 nCorrelation = temp;
2674 nEllipticity = temp;
2675 nProbability = temp;
2676 nPolarisation= temp;
2677 nAntenaPrior = temp;
2679 for(i=0; i<
n; i++) {
2681 d->
setTau(sms,t1,t2,p1,p2);
2686 skyMask.resize(m); skyMask = 1;
2687 skyMaskCC.resize(0);
2688 skyHole.resize(m); skyHole = 1.;
2700 skymap temp(healpix_order);
2701 size_t m = temp.
size();
2702 size_t n = this->ifoList.size();
2704 nSensitivity = temp;
2706 nCorrelation = temp;
2714 nEllipticity = temp;
2715 nProbability = temp;
2716 nPolarisation= temp;
2717 nAntenaPrior = temp;
2719 for(i=0; i<
n; i++) {
2721 d->
setTau(healpix_order);
2726 skyMask.resize(m); skyMask = 1;
2727 skyMaskCC.resize(0);
2728 skyHole.resize(m); skyHole = 1.;
2739 size_t N = this->ifoList.size();
2746 if(strstr(frame,
"FL") || strstr(frame,
"FS")) {
2747 tm = strstr(frame,
"FS") ? 1. : 0.;
2748 gg = strstr(frame,
"FS") ? 1. : -1.;
2751 for(n=0; n<
N; n++) {
2752 for(m=n+1; m<
N; m++) {
2753 s = ifoList[
n]->tau;
2754 s -= ifoList[
m]->tau;
2755 t = gg*(s.
max()-s.
min());
2756 if(t < tm) { tm=
t; nn =
n; mm =
m; }
2760 s = ifoList[nn]->tau;
2761 s+= ifoList[mm]->tau;
2766 else if(strstr(frame,
"BC")) {
2767 for(n=1; n<
N; n++) s += ifoList[n]->
tau;
2773 for(n=0; n<
N; n++) {
2774 if(strstr(frame,getifo(n)->Name)) this->mIFO = n;
2776 s = ifoList[this->mIFO]->tau;
2780 for(n=0; n<
N; n++) ifoList[n]->
tau -= s;
2789 size_t n = this->ifoList.size();
2790 double maxTau = -1.;
2794 if(n < 2)
return 0.;
2797 for(i=0; i<
n; i++) {
2798 tmax = ifoList[
i]->tau.max();
2799 tmin = ifoList[
i]->tau.min();
2800 if(tmax > maxTau) maxTau = tmax;
2801 if(tmin < minTau) minTau = tmin;
2803 if(strstr(name,
"min"))
return minTau;
2804 if(strstr(name,
"max"))
return maxTau;
2805 if(strstr(name,
"MAX"))
return fabs(maxTau)>
fabs(minTau) ?
fabs(maxTau) :
fabs(minTau);
2806 return (maxTau-minTau)/2.;
2844 size_t M = this->ifoList.size();
2845 if(M >
NIFO)
return;
2849 for(
size_t m=0;
m<
M;
m++) {
2850 D[
m] = this->getifo(
m);
2851 if(D[
m]->mFp.size() != D[0]->
mFp.
size()) {
2852 cout<<
"network::setIndex(): invalid detector skymaps\n";
2855 this->setAntenna(D[
m]);
2870 size_t N = ifoList.size();
2877 cout<<
"network::setDelayIndex(): invalid network\n";
2882 for(n=0; n<
N; n++) dr[n] = ifoList[n];
2890 for(n=0; n<
N; n++) {
2901 this->nPenalty = dr[0]->
tau;
2902 this->nNetIndex = dr[0]->
tau;
2911 for(n=0; n<
N; n++) {
2912 for(m=0; m<
N; m++) {
2914 i = t>0 ?
int(t*rTDF+0.5) :
int(t*rTDF-0.5);
2920 for(n=0; n<
N; n++) {
2922 for(m=0; m<
N; m++) {
2923 for(k=0; k<
N; k++) {
2924 t =
fabs(mm[n][k]-mm[n][m]-tt[m][k]);
2925 if(TT[n] < t) TT[
n] =
t;
2931 for(m=0; m<
N; m++) {
2932 if(t>TT[m]) { t = TT[
m]; k =
m; }
2934 this->nPenalty.set(l,
double(t));
2937 if(mIFO<9) i = mm[
k][this->mIFO];
2938 else i = t>0 ?
int(t*rTDF+0.5) :
int(t*rTDF-0.5);
2944 for(m=0; m<
N; m++) {
2969 size_t M = this->ifoList.size();
2975 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
2977 int type = this->
optim ? R : TYPE;
2982 vector<wavearray<double> >
snr(M);
2983 vector<wavearray<double> > nul(M);
2985 double xcut,xnul,xsnr,amin,emin,like,LIKE;
2988 for(j=0; j<nLag; j++) {
2990 if(!this->wc_List[j].
size())
continue;
2992 cid = this->wc_List[
j].
get((
char*)
"ID",0,
'S',type);
2996 siz = this->wc_List[
j].
get((
char*)
"size",0,
'S',type);
2997 SIZ = this->wc_List[
j].
get((
char*)
"SIZE",0,
'S',type);
2999 for(m=0; m<
M; m++) {
3000 snr[
m] = this->wc_List[
j].get((
char*)
"energy",m+1,
'S',type);
3001 nul[
m] = this->wc_List[
j].get((
char*)
"null",m+1,
'W',type);
3004 for(k=0; k<
K; k++) {
3006 if(siz.
data[k] < core_size)
continue;
3007 i = ID = size_t(cid.
data[k]+0.1);
3009 if(tYPe==
'i' || tYPe==
's' || tYPe==
'g' || tYPe==
'r' ||
3010 tYPe==
'I' || tYPe==
'S' || tYPe==
'G' || tYPe==
'R') {
3011 skip = !this->SETNDM(i,j,
true,type);
3013 else if(tYPe==
'b' || tYPe==
'B' || tYPe==
'E') {
3014 skip = !this->setndm(i,j,
true,type);
3021 cout<<
"network::netcut: - undefined NDM matrix"<<endl;
3022 this->wc_List[
j].ignore(ID);
3023 count += size_t(SIZ.
data[k]+0.1);
3027 xnul = xsnr = like = LIKE = 0;
3028 emin = amin = 1.e99;
3029 for(n=0; n<
M; n++) {
3030 xnul += this->getifo(n)->null;
3031 xsnr += snr[
n].data[
k];
3032 like += snr[
n].data[
k]-this->getifo(n)->null;
3033 xcut = snr[
n].data[
k]-this->getifo(n)->null;
3034 if(xcut < emin) emin = xcut;
3035 xcut = (xcut-this->getifo(n)->null)/snr[n].data[k];
3036 if(xcut < amin) amin = xcut;
3037 for(m=0; m<
M; m++) {
3038 LIKE += this->getNDM(n,m);
3042 if(cut ==
'A' || cut ==
'E') {
3043 emin = amin = 1.e99;
3044 for(n=0; n<
M; n++) {
3045 xcut = snr[
n].data[
k] - this->getifo(n)->null;
3046 if(like-xcut < emin) emin = like-xcut;
3047 xcut = 2*(like-xcut)/(xsnr-xcut) - 1.;
3048 if(xcut < amin) amin = xcut;
3053 if(cut ==
'c' || cut ==
'C') xcut = this->eCOR/(xnul +
fabs(this->eCOR));
3054 if(cut ==
'x' || cut ==
'X') xcut = this->eCOR/sqrt(
fabs(this->eCOR)*M);
3055 if(cut ==
'l' || cut ==
'L') xcut = like;
3056 if(cut ==
'a' || cut ==
'A') xcut = amin;
3057 if(cut ==
'e' || cut ==
'E') xcut = emin;
3058 if(cut ==
'r' || cut ==
'R') xcut = this->eCOR/sqrt((xnul +
fabs(this->eCOR))*M);
3061 this->wc_List[
j].ignore(ID);
3062 count += size_t(SIZ.
data[k]+0.1);
3066 if(cut==
'X') cout<<xcut<<cut<<
": size="<<SIZ.
data[
k]<<
" L="<<LIKE<<
" ID="<<ID<<endl;
3070 printf(
"%3d %4d %7.2e %7.2e %7.2e %7.2e %7.2e %7.2e %7.2e\n",
3071 (
int)n,(
int)i,getNDM(0,0),getNDM(1,1),getNDM(2,2),getNDM(0,1),getNDM(0,2),getNDM(1,2),
3072 getNDM(0,0)+getNDM(1,1)+getNDM(2,2)+2*(getNDM(0,1)+getNDM(0,2)+getNDM(1,2)));
3090 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3091 size_t R = size_t(this->ifoList[0]->
getTFmap()->
rate()/I+0.5);
3092 size_t N = this->ifoList[0]->getTFmap()->size();
3093 size_t M = this->ifoList.size();
3095 int window =
int(T*R+0.5)*I/2;
3099 int inD, frst, last, nB,
nT;
3104 std::vector<detector*> pDet; pDet.clear();
3105 std::vector<double*> pDat; pDat.clear();
3107 for(m=0; m<
M; m++) {
3108 pDet.push_back(ifoList[m]);
3109 pDat.push_back(ifoList[m]->
getTFmap()->data);
3112 size_t wmode = pDet[0]->getTFmap()->w_mode;
3113 if(wmode != 1)
return 0;
3120 for(n=0; n<nLag; n++) {
3121 V = this->wc_List[
n].pList.
size();
3122 nB =
int(this->wc_List[n].getbpp()*2*window+0.5);
3126 for(j=0; j<V; j++) {
3128 pix = &(this->wc_List[
n].pList[
j]);
3129 if(R !=
size_t(pix->
rate+0.5))
continue;
3133 for(m=0; m<
M; m++) {
3135 inD = size_t(pix->
getdata(
'I',m)+0.1);
3140 if(inD-window < offset) {
3142 last = frst + 2*window;
3144 else if(inD+window >
int(N-offset)) {
3146 frst = last - 2*window;
3153 if( inD>
int(N-offset) || inD<offset ||
3154 frst>
int(N-offset) || frst<offset ||
3155 last>
int(N-offset) || last<offset) {
3156 cout<<
"network::setRank() error\n";
3162 rank = pDet[
m]->getTFmap()->getSampleRankE(inD,frst,last);
3165 if(rank > nT-nB) rank =
log(nB/
double(nT+1-rank));
3184 size_t L = this->skyHole.size();
3193 cout<<endl<<
"network::setSkyMask() - skymap size L=0"<<endl<<endl;
3196 cout<<endl<<
"network::setSkyMask() - NULL input skymask file"<<endl<<endl;
3198 }
else if(!strlen(file)) {
3199 cout<<endl<<
"network::setSkyMask() - input skymask file not defined"<<endl<<endl;
3201 }
else if( (in=fopen(file,
"r"))==NULL ) {
3202 cout << endl <<
"network::setSkyMask() - input skymask file '"
3203 << file <<
"' not exist" << endl << endl;;
3207 while(fgets(str,1024,in) != NULL){
3209 if(str[0] ==
'#')
continue;
3210 if((pc = strtok(str,
" \t")) == NULL)
continue;
3211 if(pc) i = atoi(pc);
3212 if((pc = strtok(NULL,
" \t")) == NULL)
continue;
3213 if(pc && i>=0 && i<
int(L)) {
3214 this->skyHole.data[
i] = atof(pc);
3215 this->nSkyStat.set(i, atof(pc));
3218 cout<<endl<<
"network::setSkyMask() - "
3219 <<
"skymask file contains index > max L="<<L<<endl<<endl;
3224 cout<<endl<<
"network::setSkyMask() - "
3225 <<
"the number of indexes in the skymask file != L="<<L<<endl<<endl;
3228 a = skyHole.mean()*skyHole.size();
3229 skyHole *= a>0. ? 1./a : 0.;
3230 if(f==0.) { skyHole = 1.;
return n; }
3232 double* p = this->skyHole.data;
3233 double** pp = (
double **)malloc(L*
sizeof(
double*));
3234 for(l=0; l<
L; l++) pp[l] = p + l;
3236 skyProb.waveSort(pp,0,L-1);
3239 for(l=0; l<
L; l++) {
3241 *pp[
l] = a/L<f ? 0. : 1.;
3242 if(*pp[l] == 0.) this->nSkyStat.set(pp[l]-p,*pp[l]);
3252 size_t L = this->skyHole.size();
3259 if(skycoord!=
'e' && skycoord!=
'c') {
3260 cout <<
"network::setSkyMask() - wrong input sky coordinates "
3261 <<
" must be 'e'/'c' earth/celestial" << endl;;
3266 cout<<endl<<
"network::setSkyMaskCC() - skymap size L=0"<<endl<<endl;
3269 cout<<endl<<
"network::setSkyMaskCC() - NULL input skymask file"<<endl<<endl;
3271 }
else if(!strlen(file)) {
3272 cout<<endl<<
"network::setSkyMaskCC() - input skymask file not defined"<<endl<<endl;
3274 }
else if( (in=fopen(file,
"r"))==NULL ) {
3275 cout << endl <<
"network::setSkyMaskCC() - input skymask file '"
3276 << file <<
"' not exist" << endl << endl;;
3280 if(skycoord==
'e') {skyMask.resize(L); skyMask = 1;}
3281 if(skycoord==
'c') {skyMaskCC.resize(L); skyMaskCC = 1;}
3283 while(fgets(str,1024,in) != NULL){
3285 if(str[0] ==
'#')
continue;
3286 if((pc = strtok(str,
" \t")) == NULL)
continue;
3287 if(pc) i = atoi(pc);
3288 if((pc = strtok(NULL,
" \t")) == NULL)
continue;
3289 if(pc && i>=0 && i<
int(L)){
3291 if(skycoord==
'e') this->skyHole.data[
i]=data;
3292 if(skycoord==
'c') this->skyMaskCC.data[
i]=data;
3295 cout<<endl<<
"network::setSkyMask() - "
3296 <<
"skymask file contains index > max L="<<L<<endl<<endl;
3301 cout<<endl<<
"network::setSkyMask() - "
3302 <<
"the number of indexes in the skymask file != L="<<L<<endl<<endl;
3313 if(skycoord!=
'e' && skycoord!=
'c') {
3314 cout <<
"network::setSkyMask() - wrong input sky coordinates "
3315 <<
" must be 'e'/'c' earth/celestial" << endl;;
3319 size_t L = this->skyHole.size();
3320 if((
int)sm.
size()!=
L) {
3321 cout <<
"network::setSkyMask() - wrong input skymap size "
3322 << sm.
size() <<
" instead of " << L << endl;;
3328 for(
int i=0;
i<
L;
i++) this->skyHole.data[
i]=sm.
get(
i);
3331 skyMaskCC.resize(L);
3332 for(
int i=0;
i<
L;
i++) this->skyMaskCC.data[
i]=sm.
get(
i);
3350 std::map <string, int> mdcMap;
3352 if( (in=fopen(file,
"r"))==NULL ) {
3353 cout<<
"network::readMDClog() - no file is found \n";
3357 while(fgets(str,1024,in) != NULL){
3359 if(str[0] ==
'#')
continue;
3364 if((p = strtok(STR,
" \t")) == NULL)
continue;
3366 for(i=1; i<nTime; i++) {
3367 p = strtok(NULL,
" \t");
3373 if(gps==0. ||
fabs(GPS-gps)<7200.) {
3374 this->mdcList.push_back(str);
3375 this->mdcTime.push_back(GPS);
3381 if((p = strtok(str,
" \t")) == NULL)
continue;
3383 for(i=1; i<nName; i++) {
3384 p = strtok(NULL,
" \t");
3388 if(p)
if(mdcMap.find(p)==mdcMap.end()) mdcMap[p]=imdcMap++;
3393 this->mdcType.resize(mdcMap.size());
3394 std::map<std::string, int>::iterator iter;
3395 for (iter=mdcMap.begin(); iter!=mdcMap.end(); iter++) {
3396 this->mdcType[iter->second]=iter->first;
3399 for(
int j=0;j<this->mdcType.size();j++) {
3402 else if(j<10000) step=100;
3405 printf(
"type %3d\t",(
int)j);
3406 cout<<
" has been assigned to waveform "<<mdcType[
j]<<endl;
3410 return this->mdcList.size();
3423 if( (in=fopen(file,
"r"))==NULL ) {
3424 cout<<
"network::readSEGlist(): specified segment file "<<file<<
" does not exist\n";
3428 while(fgets(str,1024,in) != NULL){
3430 if(str[0] ==
'#')
continue;
3434 if((p = strtok(str,
" \t")) == NULL)
continue;
3436 for(i=1; i<
n; i++) {
3437 p = strtok(NULL,
" \t");
3443 SEG.
start = atof(p);
3444 p = strtok(NULL,
" \t");
3466 this->mdc__ID.clear();
3468 size_t I = this->ifoList.
size();
3479 size_t L = this->mdcList.size();
3480 size_t n = size_t(this->Edge*R+0.5);
3481 int W =
int(Tw*R/2.+0.5);
3484 if(!I || !N)
return 0.;
3486 if(this->veto.size() != size_t(N)) {
3487 this->veto.resize(N);
3492 for(k=0; k<
K; k++) {
3503 for(j=jb; j<je; j++) this->veto.data[j] = 1;
3506 if(!K) this->veto = 1;
3508 for(k=0; k<
L; k++) {
3510 if(gps == 0.)
continue;
3514 for(i=0; i<I; i++) {
3515 d = this->ifoList[
i];
3523 jm =
int((gps-S)*R);
3524 jb = jm-W; je = jm+W;
3526 if(jb >=N)
continue;
3528 if(je <=0)
continue;
3529 if(je-jb <
int(R))
continue;
3530 if(jm<jb || jm>je)
continue;
3532 for(j=jb; j<je; j++) w.
data[j] = 1;
3534 if(veto.data[jm]) this->mdc__ID.push_back(k);
3537 if(L) this->veto *=
w;
3539 for(k=n; k<N-
n; k++) live+=this->veto.
data[k];
3550 size_t M = this->ifoList.size();
3564 for(j=0; j<
id.size(); j++) {
3565 if(
id.data[j] == ID) { flag=
true;
break; }
3568 if(!flag)
return false;
3571 v = wc->
nTofF[ID-1];
3573 l = size_t(
log(k*1.)/
log(2.)+0.1);
3577 for(i=0; i<
M; i++) {
3578 pd = this->getifo(i);
3579 if(pd->
getwave(ID,*wc,atype,i) == 0.) { flag=
false;
break; }
3628 if(m > (
int)w.
size()) m = w.
size();
3643 size_t nIFO = this->ifoList.size();
3648 bool signal = (abs(atype)==
'W' || abs(atype)==
'w') ?
false :
true;
3651 for(j=0; j<
id.size(); j++) {
3652 if(
size_t(
id.data[j]+0.1) ==
ID) flag=
true;
3654 if(!flag)
return false;
3659 v = pwc->
nTofF[ID-1];
3663 for(i=0; i<
nIFO; i++) {
3666 if(x.
size() == 0.) {cout<<
"zero length\n";
return false;}
3671 double R = this->rTDF;
3672 double tShift = -v[
i]/R;
3677 for (
int ii=0;ii<(
int)x.
size()/2;ii++) {
3678 TComplex X(x.
data[2*ii],x.
data[2*ii+1]);
3679 X=X*C.Exp(TComplex(0.,-2*
PI*ii*df*tShift));
3680 x.
data[2*ii]=X.Re();
3681 x.
data[2*ii+1]=X.Im();
3686 if(signal) this->getifo(i)->waveForm =
x;
3687 else this->getifo(i)->waveBand =
x;
3703 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3704 size_t R = size_t(this->ifoList[0]->
getTFmap()->
rate()/I+0.5);
3705 size_t N = this->ifoList[0]->getTFmap()->size();
3706 size_t M = this->ifoList.size();
3707 size_t jB = size_t(this->Edge*R)*I;
3711 std::vector<detector*> pDet; pDet.clear();
3712 std::vector<double*> pDat; pDat.clear();
3713 std::vector<int> pLag; pLag.clear();
3717 pix.
rate = float(R);
3721 for(m=0; m<
M; m++) {
3722 pDet.push_back(ifoList[m]);
3723 pDat.push_back(ifoList[m]->
getTFmap()->data);
3724 pLag.push_back(
int(ifoList[m]->sHIFt*R*I+0.5));
3727 size_t il = size_t(2.*pDet[0]->TFmap.getlow()/R);
3728 size_t ih = size_t(2.*pDet[0]->TFmap.gethigh()/R);
3729 if(ih==0 || ih>=I) ih = I;
3731 size_t J = size_t(sTARt*R+0.1);
3732 size_t K = size_t(duration*R+0.1);
3735 this->wc_List[0].clear();
3739 for(i=il; i<ih; i++){
3741 S = pDet[0]->TFmap.getSlice(i);
3747 cout<<
"network::initwc() error - index out of limit \n";
3753 for(m=0; m<
M; m++) {
3754 k = pix.
time+pLag[
m];
3756 if(!this->veto.data[k]) save =
false;
3761 pix.
setdata(pDet[m]->getNoise(i,k),
'N',m);
3764 if(save) this->wc_List[0].append(pix);
3769 wc_List[0].start = pDet[0]->TFmap.start();
3770 wc_List[0].stop = N/R/I;
3771 wc_List[0].rate = pDet[0]->TFmap.rate();
3774 this->wc_List[0].pList[npix-1].neighbors[0]=0;
3775 this->wc_List[0].cluster();
3786 size_t nIFO = this->ifoList.size();
3788 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
3789 cout<<
"network::coherence(): \n"
3790 <<
"invalid number of detectors or\n"
3791 <<
"delay filter is not set\n";
3794 if(getifo(0)->getTFmap()->w_mode != 1) {
3795 cout<<
"network::coherence(): invalid whitening mode.\n";
3799 if(factor > 1.) factor = float(nIFO-1)/
nIFO;
3804 size_t jS,jj,nM,jE,LL;
3806 double R = this->ifoList[0]->getTFmap()->rate();
3807 size_t N = this->ifoList[0]->getTFmap()->size();
3808 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3809 size_t M = this->
filter.size()/I;
3810 size_t K = this->
filter[0].index.size();
3812 size_t jB = size_t(this->Edge*R/I)*I;
3813 double band = this->ifoList[0]->TFmap.rate()/I/2.;
3814 slice S = getifo(0)->getTFmap()->getSlice(0);
3826 double* pdata[
NIFO];
3828 for(n=0; n<
NIFO; n++) {
3830 if(n >= nIFO)
continue;
3831 pdata[
n] = getifo(n)->getTFmap()->
data;
3846 for(l=0; l<
L; l++)
if(skyMask.data[l]) LL++;
3847 for(n=0; n<
NIFO; n++) {
3849 if(n>=nIFO) inDEx[
n] = 0;
3852 ina = this->getifo(n)->index.
data;
3853 for(l=0; l<
L; l++) {
3854 if(skyMask.data[l]) inDEx[
n].
data[k++] = ina[
l];
3860 short* m0 = inDEx[0].data; ,
3861 short*
m1 = inDEx[1].
data; ,
3862 short*
m2 = inDEx[2].
data; ,
3863 short* m3 = inDEx[3].
data; ,
3864 short* m4 = inDEx[4].
data; ,
3865 short* m5 = inDEx[5].
data; ,
3866 short* m6 = inDEx[6].
data; ,
3867 short* m7 = inDEx[7].
data; )
3870 double* pF = (
double*)malloc(K*M*
sizeof(
double));
3871 int* pJ = (
int*)malloc(K*M*
sizeof(
int));
3880 size_t KK = this->
filter.size()/I;
3881 for(n=0; n<
nIFO; n++) {
3882 t1 = getifo(n)->tau.min()*R*getifo(0)->nDFS/I;
3883 t2 = getifo(n)->tau.max()*R*getifo(0)->nDFS/I;
3884 M1[
n] = short(t1+KK/2+0.5)-1;
3885 M2[
n] = short(t2+KK/2+0.5)+1;
3895 this->pixeLHood = getifo(0)->TFmap;
3896 this->pixeLHood = 0.;
3899 if(this->veto.size() !=
N) { veto.
resize(N); veto = 1; }
3903 for(k=0; k<nLag; k++) {
3904 this->wc_List[
k].clear();
3905 this->livTime[
k] = 0.;
3906 this->wc_List[
k].setlow(getifo(0)->TFmap.getlow());
3907 this->wc_List[
k].sethigh(getifo(0)->TFmap.gethigh());
3910 for(i=1; i<I; i++) {
3914 if(a >= getifo(0)->TFmap.gethigh())
continue;
3916 if(a <= getifo(0)->TFmap.getlow())
continue;
3924 pF[k+m*
K] = double(pv->
value[k]);
3929 S = getifo(0)->getTFmap()->getSlice(i);
3932 for(n=0; n<
NIFO; n++) {
3934 emax[
n] = 0; in[
n] = 0;
3944 for(n=0; n<
nIFO; n++) {
3948 for(j=jS; j<
N; j+=I) {
3950 emax[
n].
data[jj] = 0.;
3958 for(m=M1[n]; m<M2[
n]; m++){
3964 emax[
n].
data[jj++] = b;
3972 for(k=0; k<nLag; k++) {
3977 for(n=0; n<
nIFO; n++) {
3978 b = this->getifo(n)->lagShift.data[
k];
3979 if(a>b) { a = b; nM =
n; }
3982 for(n=0; n<
nIFO; n++) {
3983 b = this->getifo(n)->lagShift.data[
k];
3984 in[
n] = (
int((b-a)*R)+jB)/I - 1;
3987 for(jj=jB/I; jj<NN; jj++) {
3990 for(n=0; n<
nIFO; n++) {
3991 if((++in[n]) >=
int(NN)) in[n] -= jE;
3992 m += this->veto.data[inTF.
data[in[
n]]];
3993 E += emax[
n].
data[in[
n]];
3995 this->livTime[
k] += float(m/nIFO);
3996 if(E<Eo || m<nIFO)
continue;
3999 for(n=0; n<
nIFO; n++) {
4000 b = E - emax[
n].
data[in[
n]];
4001 if(b<Eo*factor) skip =
true;
4007 for(n=0; n<
nIFO; n++) {
4010 pq = pdata[
n]+inTF.
data[in[
n]];
4012 for(m=M1[n]; m<M2[
n]; m++){
4020 for(l=0; l<LL; l++) {
4023 pet+=pe[0][m0[l]]; ,
4024 pet+=pe[1][m1[
l]]; ,
4025 pet+=pe[2][m2[
l]]; ,
4026 pet+=pe[3][m3[
l]]; ,
4027 pet+=pe[4][m4[
l]]; ,
4028 pet+=pe[5][m5[
l]]; ,
4029 pet+=pe[6][m6[
l]]; ,
4030 pet+=pe[7][m7[
l]]; )
4032 { skip =
false;
break; }
4038 j = inTF.
data[in[nM]];
4039 pix.
rate = float(R/I);
4041 pix.
core = E>Es ?
true :
false;
4045 for(n=0; n<
nIFO; n++) {
4050 wc_List[
k].append(pix);
4051 if(!k) this->pixeLHood.data[
j] = sqrt(E/nIFO);
4058 for(k=0; k<nLag; k++) {
4059 a = getifo(0)->getTFmap()->start();
4060 b = getifo(0)->getTFmap()->size()/R;
4061 this->wc_List[
k].start =
a;
4062 this->wc_List[
k].
stop = a+b;
4063 this->wc_List[
k].rate = R;
4064 this->livTime[
k] *= double(I)/(Io*R);
4067 if(nZero) this->setRMS();
4085 size_t N = this->wc_List[lag].csize();
4086 size_t M = this->mdc__IDSize();
4087 size_t L = this->skyProb.size();
4088 size_t Lm = L-
int(0.9999*L);
4089 size_t nIFO = this->ifoList.size();
4090 bool prior = this->gamma<0?
true:
false;
4105 std::vector<float>* vf = &(this->wc_List[lag].sArea[
id-1]);
4108 double* p = this->skyProb.
data;
4109 double** pp = (
double **)malloc(L*
sizeof(
double*));
4110 for(l=0; l<
L; l++) pp[l] = p + l;
4112 skyProb.waveSort(pp,0,L-1);
4114 double Po = *pp[L-1];
4115 double rms =
fabs(rMs);
4117 double smax=nAntenaPrior.max();
4119 for(l=0; l<
L; l++) {
4120 if(*pp[l] <= 0.) {*pp[
l]=0.;
continue;}
4121 *pp[
l] = exp(-(Po - *pp[l])/2./rms);
4122 if(prior) *pp[
l] *= pow(nAntenaPrior.get(
int(pp[l]-p))/smax,4);
4125 if(prior) skyProb.waveSort(pp,0,L-1);
4127 for(l=0; l<
L; l++) {
4129 nProbability.set(l,p[l]);
4132 if(pOUT) cout<<rMs<<
" "<<*pp[L-1]<<
" "<<*pp[L-2]<<
" "<<*pp[L-3]<<
"\n";
4135 for(m=0; m<11; m++) { v[
m] = 0; vf->push_back(0.); }
4138 for(l=L-1; l>Lm; l--){
4140 for(m=
size_t(vol*10.)+1; m<10; m++) v[m] += 1;
4141 if(vol >= 0.9)
break;
4144 for(m=1; m<10; m++) {
4145 (*vf)[
m] = sqrt(v[m]*s);
4146 if(pOUT && !M) cout<<m<<
" error region: "<<(*vf)[
m]<<endl;
4152 std::vector<float>* vP = &(this->wc_List[lag].p_Map[
id-1]);
4153 std::vector<int>* vI = &(this->wc_List[lag].p_Ind[
id-1]);
4161 if(
nSky<0) {
char spthr[1024];
sprintf(spthr,
"0.%d",
int(abs(
nSky)));pthr=atof(spthr);}
4162 for(l=L-1; l>Lm; l--){
4164 if(
nSky==0 && (K==1000 || sum > 0.99) && K>0)
break;
4165 else if(nSky<0 && sum > pthr && K>0)
break;
4166 else if(
nSky>0 && K==
nSky && K>0)
break;
4168 vI->push_back(
int(pp[l]-p));
4169 vP->push_back(
float(*pp[l]));
4174 if(!M) { free(pp);
return; }
4177 double injTime = 1.e12;
4182 for(m=0; m<
M; m++) {
4183 mdcID = this->getmdc__ID(m);
4184 dT =
fabs(To - this->getmdcTime(mdcID));
4185 if(dT<injTime && INJ.fill_in(
this,mdcID)) {
4188 if(pOUT)
printf(
"getSkyArea: %4d %12.4f %7.3f %f \n",
int(m),To,dT,s);
4192 if(INJ.fill_in(
this,injID)) {
4196 i = this->getIndex(th,ph);
4198 vI->push_back(
int(i));
4199 vP->push_back(
float(p[i]));
4202 for(l=L-1; l>Lm; l--){
4205 if(pp[l]-p ==
int(i))
break;
4207 (*vf)[0] = sqrt(vol);
4212 printf(
"getSkyArea: %5d %12.4f %6.1f %6.1f %6.1f %6.1f %6.2f %6.2f %6.2f %7.5f, %e %d \n",
4213 int(
id),INJ.time[0]-this->getifo(0)->TFmap.start(),INJ.theta[0],INJ.phi[0],
4227 size_t N = this->wc_List[lag].csize();
4228 size_t L = this->skyProb.size();
4229 size_t Lm = L-
int(0.9999*L);
4236 double a,th,
ph,st,sp;
4237 double TH,PH,ST,SP,Eo;
4245 std::vector<float>* vf = &(this->wc_List[lag].sArea[
id-1]);
4248 double* p = this->skyProb.
data;
4249 double** pp = (
double **)malloc(L*
sizeof(
double*));
4250 for(l=0; l<
L; l++) pp[l] = p + l;
4252 skyProb.waveSplit(pp,0,L-1,Lm);
4253 skyProb.waveSort(pp,Lm,L-1);
4256 for(l=L-1; l>=Lm; l--) *pp[l] -= Eo;
4257 for(l=0; l<Lm; l++) *pp[l] = -1.e10;
4262 for(l=Lm; l<
L; l++){
4268 for(in=-1; in<2; in++) {
4270 for(im=-1; im<2; im++) {
4271 i = this->getIndex(th+in*st,ph+im*sp);
4272 if(p[i] >= p[j] || p[i] < -1.e9)
continue;
4278 for(IN=-1; IN<2; IN++) {
4280 for(IM=-1; IM<2; IM++) {
4281 k = this->getIndex(TH+IN*ST,PH+IM*SP);
4282 if(p[i] >=p[k])
continue;
4290 skyProb.waveSplit(pp,0,L-1,Lm);
4291 skyProb.waveSort(pp,Lm,L-1);
4297 for(l=L-2; l>=Lm; l--) {
4299 skyENRG.data[L-l-1] =
a;
4300 if(a < pSigma || vol<1.) {
4307 if(pOUT) cout<<sum/vol<<
" "<<*pp[L-2]<<
" "<<*pp[L-3]<<
" "<<*pp[0]<<
"\n";
4311 for(l=L-1; l>0; l--) {
4312 a = l>Lm ? (L-l-1)*Eo : 30.;
4313 if(a > 30.) a = 30.;
4317 if(pOUT) cout<<
"norm: "<<(skyProb.mean()*
L)<<endl;
4318 skyProb *= 1./(skyProb.mean()*
L);
4319 for(l=0; l<
L; l++) nProbability.set(l,log10(p[l]));
4322 for(m=0; m<11; m++) { v[
m] = 0; vf->push_back(0.); }
4325 for(l=L-1; l>Lm; l--){
4326 for(m=
size_t(sum*10.)+1; m<10; m++) v[m] += 1;
4328 if(sum >= 0.9)
break;
4331 for(m=1; m<10; m++) (*vf)[
m] = sqrt(v[m]*s);
4335 std::vector<float>* vP = &(this->wc_List[lag].p_Map[
id-1]);
4336 std::vector<int>* vI = &(this->wc_List[lag].p_Ind[
id-1]);
4344 if(
nSky<0) {
char spthr[1024];
sprintf(spthr,
"0.%d",
int(abs(
nSky)));pthr=atof(spthr);}
4345 for(l=L-1; l>Lm; l--){
4347 if(
nSky==0 && (K==1000 || sum > 0.99) && K>0)
break;
4348 else if(nSky<0 && sum > pthr && K>0)
break;
4349 else if(
nSky>0 && K==
nSky && K>0)
break;
4351 vI->push_back(
int(pp[l]-p));
4352 vP->push_back(
float(*pp[l]));
4357 size_t M = this->mdc__IDSize();
4359 if(!M) { free(pp);
return; }
4362 double injTime = 1.e12;
4367 for(m=0; m<
M; m++) {
4368 mdcID = this->getmdc__ID(m);
4369 dT =
fabs(To - this->getmdcTime(mdcID));
4370 if(dT<injTime && INJ.fill_in(
this,mdcID)) {
4373 if(pOUT)
printf(
"getSkyArea: %4d %12.4f %7.3f %f \n",
int(m),To,dT,s);
4377 if(INJ.fill_in(
this,injID)) {
4381 i = this->getIndex(th,ph);
4383 vI->push_back(
int(i));
4384 vP->push_back(
float(p[i]));
4387 for(l=L-1; l>Lm; l--){
4390 if(pp[l]-p ==
int(i))
break;
4392 (*vf)[0] = sqrt(vol);
4397 printf(
"getSkyArea: %5d %12.4f %6.1f %6.1f %6.1f %6.1f %6.2f %6.2f %6.2f %7.5f, %e %d \n",
4398 int(
id),INJ.time[0]-this->getifo(0)->TFmap.start(),INJ.theta[0],INJ.phi[0],
4417 size_t nIFO = this->ifoList.size();
4420 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
4421 cout<<
"network::likelihood(): invalid number of detectors or delay filter is not set.\n";
4424 if(type==
'b' || type==
'B' || type==
'E')
4425 return likelihoodB(type, Ao, ID, lag, ind, core);
4427 return likelihoodI(type, Ao, ID, lag, ind, core);
4443 size_t nIFO = this->ifoList.size();
4444 size_t ID = abs(iID);
4445 int N_1 = nIFO>2 ?
int(nIFO)-1 : 2;
4446 int N_2 = nIFO>2 ?
int(nIFO)-2 : 1;
4447 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
4448 cout<<
"network::likelihood(): invalid number of detectors or delay filter is not set.\n";
4456 double Eo = nIFO*Ao*Ao;
4458 double GAMMA = 1.-gamma*gamma;
4460 int LOCAL = local ? 1 : 0;
4463 double gr,gp,
gx,gR,gI,gc,gg,gP,gX,
T,rm,Le,E,LPm;
4464 double STAT,
P,EE,
Et,Lo,Lm,hh,Xp,Xx,XX,Lp,Em,Ep;
4465 double S,co,si,cc,em,um,vm,uc,us,vc,vs,Co,Cp;
4466 NETX(
double c0;,
double c1;,
double c2;,
double c3;,
double c4;,
double c5;,
double c6;,
double c7;)
4469 if(type==
'E' && Eo<nIFO) Eo = double(nIFO);
4470 if(!ID && ind>=0) ind = -1;
4476 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
4477 size_t M = this->
filter.size()/I;
4478 size_t L = ind<0 ? this->
index.
size() : ind+1;
4482 double* pdata[
NIFO];
4483 double* qdata[
NIFO];
4485 for(n=0; n<
nIFO; n++) {
4486 pdata[
n] = getifo(n)->getTFmap()->data;
4487 qdata[
n] = getifo(n)->getTFmap()->data;
4491 std::vector<float>*
F;
4492 std::vector<short>* J;
4502 for(n=0; n<
NIFO; n++) {
4503 fp[
n] = n<nIFO ? getifo(n)->fp.
data : f00.
data;
4504 fx[
n] = n<nIFO ? getifo(n)->fx.data : f00.
data;
4505 ffp[
n] = n<nIFO ? getifo(n)->ffp.data : f00.
data;
4506 ffm[
n] = n<nIFO ? getifo(n)->ffm.data : f00.
data;
4507 fpx[
n] = n<nIFO ? getifo(n)->fpx.data : f00.
data;
4510 short* mm = this->skyMask.data;
4522 std::vector<size_t> pI;
4554 for(n=0; n<
NIFO; n++) {
4563 std::vector<int>* vint;
4564 std::vector<int>* vtof;
4572 this->pixeLHood = getifo(0)->TFmap;
4573 this->pixeLHood = 0.;
4574 this->pixeLNull = getifo(0)->TFmap;
4575 this->pixeLNull = 0.;
4576 this->nSensitivity = 0.;
4577 this->nAlignment = 0.;
4578 this->nNetIndex = 0.;
4579 this->nDisbalance = 0.;
4580 this->nLikelihood = 0.;
4581 this->nNullEnergy = 0.;
4582 this->nCorrEnergy = 0.;
4583 this->nCorrelation = 0.;
4584 this->nSkyStat = 0.;
4585 this->nPenalty = 0.;
4586 this->nProbability = 0.;
4587 this->nAntenaPrior = 0.;
4594 for(n=lag; n<nLag; n++) {
4596 if(!this->wc_List[n].
size())
continue;
4597 logbpp = -
log(this->wc_List[n].getbpp());
4599 cid = this->wc_List[
n].
get((
char*)
"ID",0,
'S',
optim ? R : -R);
4600 cTo = this->wc_List[
n].
get((
char*)
"time",0,
'L',
optim ? R : -R);
4604 for(k=0; k<
K; k++) {
4606 id = size_t(cid.
data[k]+0.1);
4608 if(ID &&
id!=ID)
continue;
4610 vint = &(this->wc_List[
n].cList[
id-1]);
4611 vtof = &(this->wc_List[
n].nTofF[
id-1]);
4618 for(j=0; j<V; j++) {
4620 pix = this->wc_List[
n].getPixel(
id,j);
4623 cout<<
"network::likelihood() error: NULL pointer"<<endl;
4627 if(R !=
int(pix->
rate+0.5))
continue;
4628 if(!pix->
core && core)
continue;
4637 if(NV[0].
size() < V) {
4639 for(i=0; i<
NIFO; i++) {
4649 for(j=0; j<V; j++) {
4651 pix = this->wc_List[
n].getPixel(
id,pI[j]);
4654 for(i=0; i<
nIFO; i++) {
4655 ee[
i] = 1./pix->
data[
i].noiserms;
4659 for(i=0; i<
nIFO; i++) {
4660 nv[
i][
j] = ee[
i]*ee[
i];
4661 nr[
i][
j] = ee[
i]/sqrt(cc);
4662 qdata[
i] = pdata[
i] + pix->
data[
i].index;
4667 for(i=0; i<
nIFO; i++) {
4683 STAT = -1.e64; m=IIm=0; Lm=Em=LPm= 0.;
4690 l = ind<0 ? 0 :
ind;
4692 if(!mm[l])
continue;
4695 pe[0] = eS[0].data + m0[l]*V; ,
4696 pe[1] = eS[1].
data +
m1[
l]*V; ,
4697 pe[2] = eS[2].
data +
m2[
l]*V; ,
4698 pe[3] = eS[3].
data + m3[
l]*V; ,
4699 pe[4] = eS[4].
data + m4[
l]*V; ,
4700 pe[5] = eS[5].
data + m5[
l]*V; ,
4701 pe[6] = eS[6].
data + m6[
l]*V; ,
4702 pe[7] = eS[7].
data + m7[
l]*V; )
4704 NETX(c0=0.;,c1=0.;,c2=0.;,c3=0.;,c4=0.;,c5=0.;,c6=0.;,c7=0.;)
4706 for(j=0; j<V; j++) {
4731 E = 0.; NETX(E+=c0;,E+=
c1;,E+=
c2;,E+=
c3;,E+=c4;,E+=c5;,E+=c6;,E+=c7;) E-=i*nIFO;
4732 if(E>STAT) { m =
l; STAT = E; Lm = E; }
4734 if(NETX(E-c0>e2or &&,E-c1>e2or &&,E-c2>e2or &&,E-c3>e2or &&,E-c4>e2or &&,E-c5>e2or &&,E-c6>e2or &&,E-c7>e2or &&)
true) skip =
false;
4736 if(skip) { this->wc_List[
n].ignore(
id);
continue; }
4743 l = ind<0 ? 0 :
ind;
4746 skyProb.data[
l] = 0.;
4748 if(skyMaskCC.size()==
L) {
4753 double gpsT = cTo.
data[
k]+getifo(0)->TFmap.start();
4755 int cc_l = this->getIndex(th,ra);
4756 if (skyMaskCC.data[cc_l]<=0)
continue;
4759 if(!mm[l] && ind<0)
continue;
4762 pa[0] = aS[0].data + m0[l]*V; ,
4763 pa[1] = aS[1].
data +
m1[
l]*V; ,
4764 pa[2] = aS[2].
data +
m2[
l]*V; ,
4765 pa[3] = aS[3].
data + m3[
l]*V; ,
4766 pa[4] = aS[4].
data + m4[
l]*V; ,
4767 pa[5] = aS[5].
data + m5[
l]*V; ,
4768 pa[6] = aS[6].
data + m6[
l]*V; ,
4769 pa[7] = aS[7].
data + m7[
l]*V; )
4778 EE=Lo=Co=Ep=Lp=Cp=inet = 0.; II=0;
4780 for(j=0; j<V; j++) {
4782 Et = dotx(pa,j,pa,j);
4789 gp = dotx(Fp,Fp)+1.e-12;
4790 gx = dotx(Fx,Fx)+1.e-12;
4796 um = rotx(Fp,uc,Fx,us,u);
4797 hh = dotx(u,pa,j,e);
4799 cc = Le-dotx(e,e)/um;
4802 Lp+= Le*cc*ii/(Et-Le+1+cc);
4807 mulx(u,hh*ii/um,xi.
data+j*NIFO);
4811 cc = Ep>0 ? Co/(Ep-Lo+II+
fabs(Co)) : 0.;
4812 if(type==
'B') Lp = Lo*cc;
4814 skyProb.data[
l] = Lp/EE;
4816 if(Lp>LPm) { LPm=Lp; Em=EE; IIm=II; }
4818 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
4819 if(cc<this->
netCC || cc*Co<rho)
continue;
4827 for(j=0; j<V; j++) {
4839 px = xi.
data+j*NIFO;
4840 if(dotx(px,px)<=0)
continue;
4846 gp = dotx(Fp,Fp)+1.e-12;
4847 gx = dotx(Fx,Fx)+1.e-12;
4851 gc = sqrt(gR*gR+gI*gI);
4859 um = rotx(Fp,uc,Fx,us,u);
4860 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
4864 ii = netx(u,um,v,vm,GAMMA);
4866 if(ii<N_2 && gamma<0)
continue;
4871 uc = Xp*(gx+gg) - Xx*gI;
4872 us = Xx*(gp+gg) - Xp*gI;
4874 if(ii<N_1 && gamma!=0) {
4875 uc = Xp*(gc+gR)+Xx*gI;
4876 us = Xx*(gc-gR)+Xp*gI;
4881 um = rotx(Fp,uc,Fx,us,u);
4882 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
4887 rr[0] = LOCAL*am[0]/(am[0]*am[0]+2.)+1-LOCAL; ,
4888 rr[1] = LOCAL*am[1]/(am[1]*am[1]+2.)+1-LOCAL; ,
4889 rr[2] = LOCAL*am[2]/(am[2]*am[2]+2.)+1-LOCAL; ,
4890 rr[3] = LOCAL*am[3]/(am[3]*am[3]+2.)+1-LOCAL; ,
4891 rr[4] = LOCAL*am[4]/(am[4]*am[4]+2.)+1-LOCAL; ,
4892 rr[5] = LOCAL*am[5]/(am[5]*am[5]+2.)+1-LOCAL; ,
4893 rr[6] = LOCAL*am[6]/(am[6]*am[6]+2.)+1-LOCAL; ,
4894 rr[7] = LOCAL*am[7]/(am[7]*am[7]+2.)+1-LOCAL; )
4898 pp[0] = rr[0]*(am[0]-hh*u[0])*u[0]; ,
4899 pp[1] = rr[1]*(am[1]-hh*u[1])*u[1]; ,
4900 pp[2] = rr[2]*(am[2]-hh*u[2])*u[2]; ,
4901 pp[3] = rr[3]*(am[3]-hh*u[3])*u[3]; ,
4902 pp[4] = rr[4]*(am[4]-hh*u[4])*u[4]; ,
4903 pp[5] = rr[5]*(am[5]-hh*u[5])*u[5]; ,
4904 pp[6] = rr[6]*(am[6]-hh*u[6])*u[6]; ,
4905 pp[7] = rr[7]*(am[7]-hh*u[7])*u[7]; )
4910 qq[0] = rr[0]*((hh*u[0]-am[0])*v[0]+u[0]*u[0]*gg); ,
4911 qq[1] = rr[1]*((hh*u[1]-am[1])*v[1]+u[1]*u[1]*gg); ,
4912 qq[2] = rr[2]*((hh*u[2]-am[2])*v[2]+u[2]*u[2]*gg); ,
4913 qq[3] = rr[3]*((hh*u[3]-am[3])*v[3]+u[3]*u[3]*gg); ,
4914 qq[4] = rr[4]*((hh*u[4]-am[4])*v[4]+u[4]*u[4]*gg); ,
4915 qq[5] = rr[5]*((hh*u[5]-am[5])*v[5]+u[5]*u[5]*gg); ,
4916 qq[6] = rr[6]*((hh*u[6]-am[6])*v[6]+u[6]*u[6]*gg); ,
4917 qq[7] = rr[7]*((hh*u[7]-am[7])*v[7]+u[7]*u[7]*gg); )
4919 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
4922 em = rotx(u,co,v,si/vm,e);
4926 rm = rotx(v,co,u,-si/um,r)+1.e-24;
4930 pp[0] = rr[0]*(am[0]-hh*e[0])*e[0]; ,
4931 pp[1] = rr[1]*(am[1]-hh*e[1])*e[1]; ,
4932 pp[2] = rr[2]*(am[2]-hh*e[2])*e[2]; ,
4933 pp[3] = rr[3]*(am[3]-hh*e[3])*e[3]; ,
4934 pp[4] = rr[4]*(am[4]-hh*e[4])*e[4]; ,
4935 pp[5] = rr[5]*(am[5]-hh*e[5])*e[5]; ,
4936 pp[6] = rr[6]*(am[6]-hh*e[6])*e[6]; ,
4937 pp[7] = rr[7]*(am[7]-hh*e[7])*e[7]; )
4942 qq[0] = rr[0]*((hh*e[0]-am[0])*r[0]+e[0]*e[0]*gg); ,
4943 qq[1] = rr[1]*((hh*e[1]-am[1])*r[1]+e[1]*e[1]*gg); ,
4944 qq[2] = rr[2]*((hh*e[2]-am[2])*r[2]+e[2]*e[2]*gg); ,
4945 qq[3] = rr[3]*((hh*e[3]-am[3])*r[3]+e[3]*e[3]*gg); ,
4946 qq[4] = rr[4]*((hh*e[4]-am[4])*r[4]+e[4]*e[4]*gg); ,
4947 qq[5] = rr[5]*((hh*e[5]-am[5])*r[5]+e[5]*e[5]*gg); ,
4948 qq[6] = rr[6]*((hh*e[6]-am[6])*r[6]+e[6]*e[6]*gg); ,
4949 qq[7] = rr[7]*((hh*e[7]-am[7])*r[7]+e[7]*e[7]*gg); )
4951 co = dotx(qq,qq)/rm+dotx(pp,pp)/em+1.e-24;
4954 em = rotx(e,co,r,si/rm,e);
4957 Co+= (hh*hh-dotx(ee,ee))/em;
4963 NETX(pp[0]=0.;,pp[1]=0.;,pp[2]=0.;,pp[3]=0.;,pp[4]=0.;,pp[5]=0.;,pp[6]=0.;,pp[7]=0.;)
4964 NETX(qq[0]=0.001;,qq[1]=0.001;,qq[2]=0.001;,qq[3]=0.001;,qq[4]=0.001;,qq[5]=0.001;,qq[6]=0.001;,qq[7]=0.001;)
4965 NETX(ee[0]=0.;,ee[1]=0.;,ee[2]=0.;,ee[3]=0.;,ee[4]=0.;,ee[5]=0.;,ee[6]=0.;,ee[7]=0.;)
4967 for(j=0; j<V; j++) {
4976 S+=
fabs(pp[0]-qq[0]); ,
4977 S+=
fabs(pp[1]-qq[1]); ,
4978 S+=
fabs(pp[2]-qq[2]); ,
4979 S+=
fabs(pp[3]-qq[3]); ,
4980 S+=
fabs(pp[4]-qq[4]); ,
4981 S+=
fabs(pp[5]-qq[5]); ,
4982 S+=
fabs(pp[6]-qq[6]); ,
4983 S+=
fabs(pp[7]-qq[7]); )
4986 pp[0] /= sqrt(qq[0]); ,
4987 pp[1] /= sqrt(qq[1]); ,
4988 pp[2] /= sqrt(qq[2]); ,
4989 pp[3] /= sqrt(qq[3]); ,
4990 pp[4] /= sqrt(qq[4]); ,
4991 pp[5] /= sqrt(qq[5]); ,
4992 pp[6] /= sqrt(qq[6]); ,
4993 pp[7] /= sqrt(qq[7]); )
5017 cc = Co/(EE-Lo+
fabs(Co));
5020 skyProb.data[
l] *=
P;
5022 if(XX>=STAT) { m=
l; STAT=XX; Lm=Lo; }
5025 this->nAntenaPrior.set(l, gP/EE);
5026 this->nSensitivity.set(l, gP/EE);
5027 this->nAlignment.set(l, gX/EE);
5028 this->nNetIndex.set(l, inet/EE);
5029 this->nDisbalance.set(l, S);
5030 this->nLikelihood.set(l, Lo);
5031 this->nNullEnergy.set(l, (EE-Lo));
5032 this->nCorrEnergy.set(l, Co);
5033 this->nCorrelation.set(l, cc);
5034 this->nSkyStat.set(l, XX);
5035 this->nPenalty.set(l, P);
5036 this->nProbability.set(l,skyProb.data[l]);
5040 if(STAT<=0.001 && ind<0) {
5041 this->wc_List[
n].ignore(
id);
5049 l = ind<0 ? m :
ind;
5052 pa[0] = aS[0].data + m0[l]*V; ,
5053 pa[1] = aS[1].
data +
m1[
l]*V; ,
5054 pa[2] = aS[2].
data +
m2[
l]*V; ,
5055 pa[3] = aS[3].
data + m3[
l]*V; ,
5056 pa[4] = aS[4].
data + m4[
l]*V; ,
5057 pa[5] = aS[5].
data + m5[
l]*V; ,
5058 pa[6] = aS[6].
data + m6[
l]*V; ,
5059 pa[7] = aS[7].
data + m7[
l]*V; )
5063 for(j=0; j<V; j++) {
5066 for(i=0; i<
nIFO; i++) {
5068 Fp[
i] = fp[
i][
l]*nr[
i][
j];
5069 Fx[
i] = fx[
i][
l]*nr[
i][
j];
5074 gp = dotx(Fp,Fp)+1.e-12;
5075 gx = dotx(Fx,Fx)+1.e-12;
5080 gc = sqrt(gR*gR+gI*gI);
5085 um = rotx(Fp,uc,Fx,us,u);
5086 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5089 if((hh*hh-dotx(e,e))/um<=0.) U=0;
5094 for(i=0; i<
nIFO; i++) {
5095 if(u[i]*u[i]/um > 1-GAMMA) ii++;
5096 if(u[i]*u[i]/um+v[i]*v[i]/vm > GAMMA) ii--;
5099 if(ii<N_2 && gamma<0.) U = 0;
5102 uc = Xp*(gx+gg) - Xx*gI;
5103 us = Xx*(gp+gg) - Xp*gI;
5105 if(ii<N_1 && gamma!=0) {
5106 uc = Xp*(gc+gR)+Xx*gI;
5107 us = Xx*(gc-gR)+Xp*gI;
5110 vc = (gp*uc + gI*us);
5111 vs = (gx*us + gI*uc);
5112 um = rotx(Fp,uc,Fx,us,u);
5113 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5121 for(i=0; i<
nIFO; i++) {
5131 for(i=0; i<
nIFO; i++) {
5132 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
5133 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
5134 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
5135 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
5138 cc = sqrt(si*si+co*co)+1.e-24;
5146 for(i=0; i<
nIFO; i++) {
5147 e[
i] = u[
i]*co + v[
i]*si;
5148 r[
i] = v[
i]*co - u[
i]*si;
5156 for(i=0; i<
nIFO; i++) {
5157 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
5158 pp[
i] = cc*(am[
i]-hh*e[
i])*e[i];
5159 qq[
i] = cc*((2.*hh*e[
i]-am[
i])*r[i] + e[i]*e[i]*gg);
5160 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
5163 cc = sqrt(si*si+co*co)+1.e-24;
5168 if(type==
'E') U = 0;
5170 for(i=0; i<
nIFO; i++) {
5171 e[
i] = e[
i]*co+r[
i]*si;
5180 Lo += (type==
'E') ? Et-nIFO : Lp;
5184 pix = this->wc_List[
n].getPixel(
id,pI[j]);
5185 pix->
likelihood = (type==
'E') ? Et-nIFO : Lp;
5186 pix->
theta = nLikelihood.getTheta(l);
5187 pix->
phi = nLikelihood.getPhi(l);
5188 if(!core) pix->
core = Et<Eo ?
false :
true;
5190 for(i=0; i<
nIFO; i++) {
5192 pix->
setdata(e[i]*hh/sqrt(nv[i][j]),
'W',i);
5196 this->pixeLHood.data[pix->
time] = Lp;
5197 this->pixeLNull.data[pix->
time] = Et-Lp+1.;
5199 cout<<j<<
" SNR="<<Et<<
" ";
5200 for(i=0; i<
nIFO; i++) {
5201 cout<<e[
i]*e[
i]<<
":"<<am[
i]*am[
i]/Et<<
" ";
5213 vtof->push_back(
int(M/2)-
int(m0[l])); ,
5214 vtof->push_back(
int(M/2)-
int(
m1[l])); ,
5215 vtof->push_back(
int(M/2)-
int(
m2[l])); ,
5216 vtof->push_back(
int(M/2)-
int(m3[l])); ,
5217 vtof->push_back(
int(M/2)-
int(m4[l])); ,
5218 vtof->push_back(
int(M/2)-
int(m5[l])); ,
5219 vtof->push_back(
int(M/2)-
int(m6[l])); ,
5220 vtof->push_back(
int(M/2)-
int(m7[l])); )
5223 T = cTo.
data[
k]+getifo(0)->TFmap.start();
5224 if(iID<=0 && type!=
'E') getSkyArea(
id,n,T);
5226 if(
fabs((Lm-Lo)/Lo)>1.e-4)
5227 cout<<
"likelihood: incorrect likelihood : "<<Lm<<
" "<<Lo<<
" "<<Em<<endl;
5230 cout<<
"max value: "<<STAT<<
" at (theta,phi) = ("<<nLikelihood.getTheta(l)
5231 <<
","<<nLikelihood.getPhi(l)<<
") Likelihood: loop: "<<Lm
5232 <<
", final: "<<Lo<<
", eff: "<<Em<<endl;
5237 this->nSensitivity.gps =
T;
5238 this->nAlignment.gps =
T;
5239 this->nNetIndex.gps =
T;
5240 this->nDisbalance.gps =
T;
5241 this->nLikelihood.gps =
T;
5242 this->nNullEnergy.gps =
T;
5243 this->nCorrEnergy.gps =
T;
5244 this->nCorrelation.gps =
T;
5245 this->nSkyStat.gps =
T;
5246 this->nPenalty.gps =
T;
5247 this->nProbability.gps =
T;
5270 size_t nIFO = this->ifoList.size();
5271 size_t ID = abs(iID);
5272 int N_1 = nIFO>2 ?
int(nIFO)-1 : 2;
5273 int N_2 = nIFO>2 ?
int(nIFO)-2 : 1;
5274 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
5275 cout<<
"network::likelihoodE(): invalid number of detectors or delay filter is not set.\n";
5283 double Eo = nIFO*Ao*Ao;
5285 double GAMMA = 1.-gamma*gamma;
5286 int LOCAL = local ? 1 : 0;
5289 double one = (type==
'g' || type==
'G') ? 0. : 1.;
5290 double En = (type!=
'i' && type!=
'I') ? 1.e9 : 1.;
5294 if(type==
'I' || type==
'S' || type==
'G') ISG =
true;
5295 if(type==
'i' || type==
's' || type==
'g') ISG =
true;
5296 if(type==
'i' || type==
's' || type==
'g' || type==
'r') isgr =
true;
5298 if(!ID && ind>=0) ind = -1;
5303 double gr,gg,gc,gp,
gx,cc,gP,gX,gR,gI,EE,
Et,
T,Em,LPm,HH;
5304 double em,
STAT,Lo,Lm,Lp,Co,E00,E90,L00,C00,AA,CO,SI;
5305 double co,si,stat,XP,XX,L90,C90,
sI,cO,aa,bb,as,ac,ap,La;
5306 double xp,
xx,hh,um,vm,uc,us,vc,vs,
hp,
hx,HP,HX,wp,wx,WP,WX;
5307 size_t i,
j,
k,
l,
m,
n,V,U,
K,
id,j5;
5309 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
5310 size_t M = this->
filter.size()/I;
5311 size_t L = ind<0 ? this->
index.
size() : ind+1;
5315 double* pdata[
NIFO];
5316 double* qdata[
NIFO];
5318 for(n=0; n<
nIFO; n++) {
5319 pdata[
n] = getifo(n)->getTFmap()->data;
5320 qdata[
n] = getifo(n)->getTFmap()->data;
5324 std::vector<float>* F00;
5325 std::vector<short>* J00;
5326 std::vector<float>* F90;
5327 std::vector<short>* J90;
5337 for(n=0; n<
NIFO; n++) {
5338 fp[
n] = n<nIFO ? getifo(n)->fp.
data : f00.
data;
5339 fx[
n] = n<nIFO ? getifo(n)->fx.data : f00.
data;
5340 ffp[
n] = n<nIFO ? getifo(n)->ffp.data : f00.
data;
5341 ffm[
n] = n<nIFO ? getifo(n)->ffm.data : f00.
data;
5342 fpx[
n] = n<nIFO ? getifo(n)->fpx.data : f00.
data;
5344 short* mm = skyMask.data;
5356 std::vector<size_t> pI;
5397 for(n=0; n<
NIFO; n++) {
5408 std::vector<int>* vint;
5409 std::vector<int>* vtof;
5416 this->pixeLHood = getifo(0)->TFmap;
5417 this->pixeLHood = 0.;
5418 this->pixeLNull = getifo(0)->TFmap;
5419 this->pixeLNull = 0.;
5420 this->nSensitivity = 0.;
5421 this->nAlignment = 0.;
5422 this->nNetIndex = 0.;
5423 this->nDisbalance = 0.;
5424 this->nLikelihood = 0.;
5425 this->nNullEnergy = 0.;
5426 this->nCorrEnergy = 0.;
5427 this->nCorrelation = 0.;
5428 this->nSkyStat = 0.;
5429 this->nEllipticity = 0.;
5430 this->nPolarisation = 0.;
5431 this->nProbability = 0.;
5432 this->nAntenaPrior = 0.;
5439 for(n=lag; n<nLag; n++) {
5441 if(!this->wc_List[n].
size())
continue;
5442 logbpp = -
log(this->wc_List[n].getbpp());
5444 cid = this->wc_List[
n].
get((
char*)
"ID",0,
'S',
optim ? R : -R);
5445 cTo = this->wc_List[
n].
get((
char*)
"time",0,
'L',
optim ? R : -R);
5449 for(k=0; k<
K; k++) {
5451 id = size_t(cid.
data[k]+0.1);
5453 if(ID &&
id!=ID)
continue;
5455 vint = &(this->wc_List[
n].cList[
id-1]);
5456 vtof = &(this->wc_List[
n].nTofF[
id-1]);
5463 for(j=0; j<V; j++) {
5465 pix = this->wc_List[
n].getPixel(
id,j);
5468 cout<<
"network::likelihood() error: NULL pointer"<<endl;
5472 if(R !=
int(pix->
rate+0.5))
continue;
5473 if(!pix->
core && core)
continue;
5482 if(NV[0].
size() < V) {
5484 for(i=0; i<
NIFO; i++) {
5498 for(j=0; j<V; j++) {
5500 pix = this->wc_List[
n].getPixel(
id,pI[j]);
5503 for(i=0; i<
nIFO; i++) {
5504 ee[
i] = 1./pix->
data[
i].noiserms;
5508 GG.
data[
j] = sqrt(cc);
5509 for(i=0; i<
nIFO; i++) {
5510 nv[
i][
j] = ee[
i]*ee[
i];
5513 qdata[
i] = pdata[
i] + pix->
data[
i].index;
5518 for(i=0; i<
nIFO; i++) {
5528 gg = dot32(F00,pq,J00);
5532 gg = dot32(F90,pq,J90);
5539 STAT = -1.e64; m=U=IIm=0; Em=Lm=LPm=AA=SI=CO=0.;
5545 l = ind<0 ? 0 :
ind;
5548 skyProb.data[
l] = 0.;
5550 if(skyMaskCC.size()==
L) {
5555 double gpsT = cTo.
data[
k]+getifo(0)->TFmap.start();
5557 int cc_l = this->getIndex(th,ra);
5558 if (!skyMaskCC.data[cc_l])
continue;
5561 if(!mm[l] && ind<0)
continue;
5564 pa00[0] = a00[0].data + m0[l]*V; ,
5565 pa00[1] = a00[1].
data +
m1[
l]*V; ,
5566 pa00[2] = a00[2].
data +
m2[
l]*V; ,
5567 pa00[3] = a00[3].
data + m3[
l]*V; ,
5568 pa00[4] = a00[4].
data + m4[
l]*V; ,
5569 pa00[5] = a00[5].
data + m5[
l]*V; ,
5570 pa00[6] = a00[6].
data + m6[
l]*V; ,
5571 pa00[7] = a00[7].
data + m7[
l]*V; )
5573 EE=E00=L00=C00=Lp=0.; II=0;
5575 for(j=0; j<V; j++) {
5589 gp = dotx(Fp,Fp)+1.e-12;
5590 gx = dotx(Fx,Fx)+1.e-12;
5596 um = rotx(Fp,xp*gx-xx*gI,
5600 Co = La - dotx(e,e)/um;
5603 Lp+= La*Co*ii/(Et-La+1+Co);
5611 cc = E00>0 ? C00/(E00-L00+
fabs(C00)) : 0.;
5612 if(cc<this->
netCC || cc*C00<rho)
continue;
5614 E90=L90=C90=inet=0.;
5617 pa90[0] = a90[0].data + m0[l]*V; ,
5618 pa90[1] = a90[1].
data +
m1[
l]*V; ,
5619 pa90[2] = a90[2].
data +
m2[
l]*V; ,
5620 pa90[3] = a90[3].
data + m3[
l]*V; ,
5621 pa90[4] = a90[4].
data + m4[
l]*V; ,
5622 pa90[5] = a90[5].
data + m5[
l]*V; ,
5623 pa90[6] = a90[6].
data + m6[
l]*V; ,
5624 pa90[7] = a90[7].
data + m7[
l]*V; )
5626 for(j=0; j<V; j++) {
5637 gp = dotx(Fp,Fp)+1.e-12;
5638 gx = dotx(Fx,Fx)+1.e-12;
5644 vm = rotx(Fp,XP*gx-XX*gI,
5648 Co = La - dotx(e,e)/vm;
5651 Lp+= La*Co*ii/(Et-La+1+Co);
5653 LL.
data[
j] += La*ii;
5659 cc = E90>0 ? C90/(E90-L90+
fabs(C90)) : 0.;
5660 if(cc<this->
netCC || cc*C90<rho)
continue;
5675 for(j=0; j<V; j++) {
5680 gp = dotx(Fp,Fp)+1.e-12;
5681 gx = dotx(Fx,Fx)+1.e-12;
5685 wp = dotx(Fp,xi00.
data+j5);
5686 wx = dotx(Fx,xi00.
data+j5);
5687 WP = dotx(Fp,xi90.
data+j5);
5688 WX = dotx(Fx,xi90.
data+j5);
5690 xp = (wp*wp + WP*WP);
5691 xx = (wx*wx + WX*WX);
5694 sI += 2*(wp*wx+WP*WX-La*gI)/gg;
5695 cO += (xp-xx - La*(gp-
gx))/gg;
5696 bb += La - (xp+
xx)/gg;
5697 aa += 2*(wp*WX - wx*WP)/gg;
5700 gg = sqrt(cO*cO+sI*sI);
5701 cO = cO/gg; sI = sI/gg;
5703 if(aa> one) aa = 1.;
5704 if(aa<-one) aa =-1.;
5706 for(j=0; j<V; j++) {
5710 um = rotx(Fp,1+cO,Fx, sI,u);
5711 vm = rotx(Fx,1+cO,Fp,-sI,v);
5712 gp = dotx(u,u)+1.e-12;
5713 gx = dotx(v,v)+1.e-12;
5719 um = dotx(xi,xi)+1.e-6;
5720 hh = dotx(xi,am,pp);
5723 bb = (xp*xp+xx*
xx)/(um*gg*um);
5725 cc = bb*dotx(pp,pp);
5729 vm = dotx(XI,XI)+1.e-6;
5730 HH = dotx(XI,AM,qq);
5733 bb = (XP*XP+XX*XX)/(vm*gg*vm);
5735 cc+= bb*dotx(qq,qq);
5737 bb = 2*(xp*XX-xx*XP)/(um*gg*vm);
5739 cc+= bb*dotx(pp,qq);
5743 Et = dotx(am,am) + dotx(AM,AM);
5744 Lp+= La*cc/(Et-La+2+
fabs(cc));
5748 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
5750 if(!isgr) Lp = Lo*cc;
5751 if(Lp>LPm) {LPm=Lp; Em=EE; IIm=II;}
5753 skyProb.data[
l] = Lp/EE;
5755 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
5756 if(cc<this->
netCC || cc*Co<rho*2)
continue;
5764 for(j=0; j<V; j++) {
5768 if(dotx(xi,xi)<=0.)
continue;
5774 gp = dotx(Fp,Fp)+1.e-12;
5775 gx = dotx(Fx,Fx)+1.e-12;
5783 um = rotx(Fp,uc,Fx,us,u);
5784 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5788 ii = netx(u,um,v,vm,GAMMA);
5790 if(ii<N_2 && gamma<0.) {inix(xi,0);
continue;}
5796 if(ii<N_1 && gamma!=0) {
5798 gc = sqrt(gR*gR+gI*gI);
5799 uc = xp*(gc+gR)+xx*gI;
5800 us = xx*(gc-gR)+xp*gI;
5805 um = rotx(Fp,uc,Fx,us,u);
5806 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5811 rr[0] = LOCAL*am[0]/(am[0]*am[0]+2.)+1-LOCAL; ,
5812 rr[1] = LOCAL*am[1]/(am[1]*am[1]+2.)+1-LOCAL; ,
5813 rr[2] = LOCAL*am[2]/(am[2]*am[2]+2.)+1-LOCAL; ,
5814 rr[3] = LOCAL*am[3]/(am[3]*am[3]+2.)+1-LOCAL; ,
5815 rr[4] = LOCAL*am[4]/(am[4]*am[4]+2.)+1-LOCAL; ,
5816 rr[5] = LOCAL*am[5]/(am[5]*am[5]+2.)+1-LOCAL; ,
5817 rr[6] = LOCAL*am[6]/(am[6]*am[6]+2.)+1-LOCAL; ,
5818 rr[7] = LOCAL*am[7]/(am[7]*am[7]+2.)+1-LOCAL; )
5822 pp[0] = rr[0]*(am[0]-hh*u[0])*u[0]; ,
5823 pp[1] = rr[1]*(am[1]-hh*u[1])*u[1]; ,
5824 pp[2] = rr[2]*(am[2]-hh*u[2])*u[2]; ,
5825 pp[3] = rr[3]*(am[3]-hh*u[3])*u[3]; ,
5826 pp[4] = rr[4]*(am[4]-hh*u[4])*u[4]; ,
5827 pp[5] = rr[5]*(am[5]-hh*u[5])*u[5]; ,
5828 pp[6] = rr[6]*(am[6]-hh*u[6])*u[6]; ,
5829 pp[7] = rr[7]*(am[7]-hh*u[7])*u[7]; )
5834 qq[0] = rr[0]*((hh*u[0]-am[0])*v[0]+u[0]*u[0]*gg); ,
5835 qq[1] = rr[1]*((hh*u[1]-am[1])*v[1]+u[1]*u[1]*gg); ,
5836 qq[2] = rr[2]*((hh*u[2]-am[2])*v[2]+u[2]*u[2]*gg); ,
5837 qq[3] = rr[3]*((hh*u[3]-am[3])*v[3]+u[3]*u[3]*gg); ,
5838 qq[4] = rr[4]*((hh*u[4]-am[4])*v[4]+u[4]*u[4]*gg); ,
5839 qq[5] = rr[5]*((hh*u[5]-am[5])*v[5]+u[5]*u[5]*gg); ,
5840 qq[6] = rr[6]*((hh*u[6]-am[6])*v[6]+u[6]*u[6]*gg); ,
5841 qq[7] = rr[7]*((hh*u[7]-am[7])*v[7]+u[7]*u[7]*gg); )
5843 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
5846 em = rotx(u,co,v,si/vm,e);
5850 vm = rotx(v,co,u,-si/um,v)+1.e-24;
5854 pp[0] = rr[0]*(am[0]-hh*e[0])*e[0]; ,
5855 pp[1] = rr[1]*(am[1]-hh*e[1])*e[1]; ,
5856 pp[2] = rr[2]*(am[2]-hh*e[2])*e[2]; ,
5857 pp[3] = rr[3]*(am[3]-hh*e[3])*e[3]; ,
5858 pp[4] = rr[4]*(am[4]-hh*e[4])*e[4]; ,
5859 pp[5] = rr[5]*(am[5]-hh*e[5])*e[5]; ,
5860 pp[6] = rr[6]*(am[6]-hh*e[6])*e[6]; ,
5861 pp[7] = rr[7]*(am[7]-hh*e[7])*e[7]; )
5866 qq[0] = rr[0]*((hh*e[0]-am[0])*v[0]+e[0]*e[0]*gg); ,
5867 qq[1] = rr[1]*((hh*e[1]-am[1])*v[1]+e[1]*e[1]*gg); ,
5868 qq[2] = rr[2]*((hh*e[2]-am[2])*v[2]+e[2]*e[2]*gg); ,
5869 qq[3] = rr[3]*((hh*e[3]-am[3])*v[3]+e[3]*e[3]*gg); ,
5870 qq[4] = rr[4]*((hh*e[4]-am[4])*v[4]+e[4]*e[4]*gg); ,
5871 qq[5] = rr[5]*((hh*e[5]-am[5])*v[5]+e[5]*e[5]*gg); ,
5872 qq[6] = rr[6]*((hh*e[6]-am[6])*v[6]+e[6]*e[6]*gg); ,
5873 qq[7] = rr[7]*((hh*e[7]-am[7])*v[7]+e[7]*e[7]*gg); )
5875 co = dotx(qq,qq)/vm+dotx(pp,pp)/em+1.e-24;
5878 em = rotx(e,co,v,si/vm,e);
5881 cc = La-dotx(ee,ee)/em;
5884 Lp+= La*cc*
int(cc>0)/(Et-La+1+cc);
5892 for(j=0; j<V; j++) {
5896 if(dotx(XI,XI)<=0)
continue;
5902 gp = dotx(Fp,Fp)+1.e-12;
5903 gx = dotx(Fx,Fx)+1.e-12;
5911 um = rotx(Fp,uc,Fx,us,u);
5912 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5916 ii = netx(u,um,v,vm,GAMMA);
5918 if(ii<N_2 && gamma<0.) {inix(XI,0.);
continue;}
5924 if(ii<N_1 && gamma!=0) {
5926 gc = sqrt(gR*gR+gI*gI);
5927 uc = XP*(gc+gR)+XX*gI;
5928 us = XX*(gc-gR)+XP*gI;
5933 um = rotx(Fp,uc,Fx,us,u);
5934 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5939 rr[0] = LOCAL*AM[0]/(AM[0]*AM[0]+2.)+1-LOCAL; ,
5940 rr[1] = LOCAL*AM[1]/(AM[1]*AM[1]+2.)+1-LOCAL; ,
5941 rr[2] = LOCAL*AM[2]/(AM[2]*AM[2]+2.)+1-LOCAL; ,
5942 rr[3] = LOCAL*AM[3]/(AM[3]*AM[3]+2.)+1-LOCAL; ,
5943 rr[4] = LOCAL*AM[4]/(AM[4]*AM[4]+2.)+1-LOCAL; ,
5944 rr[5] = LOCAL*AM[5]/(AM[5]*AM[5]+2.)+1-LOCAL; ,
5945 rr[6] = LOCAL*AM[6]/(AM[6]*AM[6]+2.)+1-LOCAL; ,
5946 rr[7] = LOCAL*AM[7]/(AM[7]*AM[7]+2.)+1-LOCAL; )
5950 pp[0] = rr[0]*(AM[0]-hh*u[0])*u[0]; ,
5951 pp[1] = rr[1]*(AM[1]-hh*u[1])*u[1]; ,
5952 pp[2] = rr[2]*(AM[2]-hh*u[2])*u[2]; ,
5953 pp[3] = rr[3]*(AM[3]-hh*u[3])*u[3]; ,
5954 pp[4] = rr[4]*(AM[4]-hh*u[4])*u[4]; ,
5955 pp[5] = rr[5]*(AM[5]-hh*u[5])*u[5]; ,
5956 pp[6] = rr[6]*(AM[6]-hh*u[6])*u[6]; ,
5957 pp[7] = rr[7]*(AM[7]-hh*u[7])*u[7]; )
5962 qq[0] = rr[0]*((hh*u[0]-AM[0])*v[0]+u[0]*u[0]*gg); ,
5963 qq[1] = rr[1]*((hh*u[1]-AM[1])*v[1]+u[1]*u[1]*gg); ,
5964 qq[2] = rr[2]*((hh*u[2]-AM[2])*v[2]+u[2]*u[2]*gg); ,
5965 qq[3] = rr[3]*((hh*u[3]-AM[3])*v[3]+u[3]*u[3]*gg); ,
5966 qq[4] = rr[4]*((hh*u[4]-AM[4])*v[4]+u[4]*u[4]*gg); ,
5967 qq[5] = rr[5]*((hh*u[5]-AM[5])*v[5]+u[5]*u[5]*gg); ,
5968 qq[6] = rr[6]*((hh*u[6]-AM[6])*v[6]+u[6]*u[6]*gg); ,
5969 qq[7] = rr[7]*((hh*u[7]-AM[7])*v[7]+u[7]*u[7]*gg); )
5971 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
5974 em = rotx(u,co,v,si/vm,e);
5978 vm = rotx(v,co,u,-si/um,v)+1.e-24;
5982 pp[0] = rr[0]*(AM[0]-hh*e[0])*e[0]; ,
5983 pp[1] = rr[1]*(AM[1]-hh*e[1])*e[1]; ,
5984 pp[2] = rr[2]*(AM[2]-hh*e[2])*e[2]; ,
5985 pp[3] = rr[3]*(AM[3]-hh*e[3])*e[3]; ,
5986 pp[4] = rr[4]*(AM[4]-hh*e[4])*e[4]; ,
5987 pp[5] = rr[5]*(AM[5]-hh*e[5])*e[5]; ,
5988 pp[6] = rr[6]*(AM[6]-hh*e[6])*e[6]; ,
5989 pp[7] = rr[7]*(AM[7]-hh*e[7])*e[7]; )
5994 qq[0] = rr[0]*((hh*e[0]-AM[0])*v[0]+e[0]*e[0]*gg); ,
5995 qq[1] = rr[1]*((hh*e[1]-AM[1])*v[1]+e[1]*e[1]*gg); ,
5996 qq[2] = rr[2]*((hh*e[2]-AM[2])*v[2]+e[2]*e[2]*gg); ,
5997 qq[3] = rr[3]*((hh*e[3]-AM[3])*v[3]+e[3]*e[3]*gg); ,
5998 qq[4] = rr[4]*((hh*e[4]-AM[4])*v[4]+e[4]*e[4]*gg); ,
5999 qq[5] = rr[5]*((hh*e[5]-AM[5])*v[5]+e[5]*e[5]*gg); ,
6000 qq[6] = rr[6]*((hh*e[6]-AM[6])*v[6]+e[6]*e[6]*gg); ,
6001 qq[7] = rr[7]*((hh*e[7]-AM[7])*v[7]+e[7]*e[7]*gg); )
6003 co = dotx(qq,qq)/vm+dotx(pp,pp)/em+1.e-24;
6006 em = rotx(e,co,v,si/vm,e);
6009 cc = La-dotx(ee,ee)/em;
6012 Lp+= La*cc*
int(cc>0)/(Et-La+1+cc);
6024 for(j=0; j<V; j++) {
6028 um = rotx(Fp,1+cO,Fx, sI,u);
6029 vm = rotx(Fx,1+cO,Fp,-sI,v);
6030 gp = dotx(u,u)+1.e-12;
6031 gx = dotx(v,v)+1.e-12;
6036 um = dotx(xi,xi)+1.e-6;
6037 hh = dotx(xi,am,pp);
6040 bb = (xp*xp+xx*
xx)/(um*gg*um);
6042 cc = bb*dotx(pp,pp);
6046 vm = dotx(XI,XI)+1.e-6;
6047 HH = dotx(XI,AM,qq);
6050 bb = (XP*XP+XX*XX)/(vm*gg*vm);
6052 cc+= bb*dotx(qq,qq);
6054 bb = 2*(xp*XX-xx*XP)/(um*gg*vm);
6056 cc+= bb*dotx(pp,qq);
6060 Et = dotx(am,am)+dotx(AM,AM);
6061 Lp+= La*cc*
int(cc>0)/(Et-La+2+cc);
6067 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
6068 stat = isgr ? Lp/EE : Lo*cc/EE;
6072 if(stat>=STAT) { m=
l; STAT=stat; AA=aa; CO=cO; SI=
sI; Lm=Lo/2.; }
6075 for(j=0; j<V; j++) {
6078 Et = dotx(pa00,j,pa00,j);
6079 gp = dotx(Fp,Fp)+1.e-12;
6080 gx = dotx(Fx,Fx)+1.e-12;
6084 gc = sqrt(gR*gR+gI*gI);
6089 this->nAntenaPrior.set(l, gP/V);
6090 this->nSensitivity.set(l, gP/V);
6091 this->nAlignment.set(l, gX/V);
6092 this->nLikelihood.set(l, Lo/2);
6093 this->nNullEnergy.set(l, (EE-Lo)/2.);
6094 this->nCorrEnergy.set(l, Co/2);
6095 this->nCorrelation.set(l, cc);
6096 this->nSkyStat.set(l, stat);
6097 this->nProbability.set(l, skyProb.data[l]);
6098 this->nDisbalance.set(l,
fabs(L00-L90)/
fabs(L00+L90));
6099 this->nEllipticity.set(l, aa);
6100 this->nPolarisation.set(l, atan2(sI,cO));
6101 this->nNetIndex.set(l, inet);
6105 if(STAT<0.001 && ind<0) {
6106 this->wc_List[
n].ignore(
id);
6114 l = ind<0 ? m :
ind;
6117 pa00[0] = a00[0].data + m0[l]*V; ,
6118 pa00[1] = a00[1].
data +
m1[
l]*V; ,
6119 pa00[2] = a00[2].
data +
m2[
l]*V; ,
6120 pa00[3] = a00[3].
data + m3[
l]*V; ,
6121 pa00[4] = a00[4].
data + m4[
l]*V; ,
6122 pa00[5] = a00[5].
data + m5[
l]*V; ,
6123 pa00[6] = a00[6].
data + m6[
l]*V; ,
6124 pa00[7] = a00[7].
data + m7[
l]*V; )
6126 pa90[0] = a90[0].data + m0[l]*V; ,
6127 pa90[1] = a90[1].
data +
m1[
l]*V; ,
6128 pa90[2] = a90[2].
data +
m2[
l]*V; ,
6129 pa90[3] = a90[3].
data + m3[
l]*V; ,
6130 pa90[4] = a90[4].
data + m4[
l]*V; ,
6131 pa90[5] = a90[5].
data + m5[
l]*V; ,
6132 pa90[6] = a90[6].
data + m6[
l]*V; ,
6133 pa90[7] = a90[7].
data + m7[
l]*V; )
6137 for(j=0; j<V; j++) {
6139 pix = this->wc_List[
n].getPixel(
id,pI[j]);
6140 pix->
theta = nLikelihood.getTheta(l);
6141 pix->
phi = nLikelihood.getPhi(l);
6145 E00=E90=gp=gx=gI=xp=xx=XP=XX = 0.;
6153 for(i=0; i<
nIFO; i++) {
6154 Fp[
i] = fp[
i][
l]*nr[
i][
j];
6155 Fx[
i] = fx[
i][
l]*nr[
i][
j];
6161 gp += Fp[
i]*Fp[
i]+1.e-12;
6162 gx += Fx[
i]*Fx[
i]+1.e-12;
6178 gc = sqrt(gR*gR+gI*gI);
6185 uc = (xp*gx - xx*gI);
6186 us = (xx*gp - xp*gI);
6189 um = rotx(Fp,uc,Fx,us,u);
6190 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6194 if((hh*hh-cc)/um<=0. || E00<Eo) U=0;
6198 ii = netx(u,um,v,vm,GAMMA);
6199 if(ii<N_2 && gamma<0.) U=0;
6204 if(ii<N_1 && gamma!=0) {
6205 uc = xp*(gc+gR)+xx*gI;
6206 us = xx*(gc-gR)+xp*gI;
6211 um = rotx(Fp,uc,Fx,us,u);
6212 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6220 for(i=0; i<
nIFO; i++) {
6230 for(i=0; i<
nIFO; i++) {
6231 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
6232 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
6233 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
6234 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6237 cc = sqrt(si*si+co*co)+1.e-24;
6245 for(i=0; i<
nIFO; i++) {
6246 e[
i] = u[
i]*co + v[
i]*si;
6247 v[
i] = v[
i]*co - u[
i]*si;
6256 for(i=0; i<
nIFO; i++) {
6257 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
6258 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
6259 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
6260 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6263 cc = sqrt(si*si+co*co)+1.e-24;
6269 for(i=0; i<
nIFO; i++) {
6270 u[
i] = u[
i]*co + v[
i]*si;
6277 for(i=0; i<
nIFO; i++) {
6278 pix->
setdata(u[i]*hh/sqrt(nv[i][j]),
'W',i);
6279 wp += Fp[
i]*u[
i]*hh;
6280 wx += Fx[
i]*u[
i]*hh;
6293 um = rotx(Fp,uc,Fx,us,u);
6294 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6298 if((hh*hh-cc)/um<=0. || E90<Eo) U=0;
6302 ii = netx(u,um,v,vm,GAMMA);
6303 if(ii<N_2 && gamma<0.) U=0;
6308 if(ii<N_1 && gamma!=0) {
6309 uc = XP*(gc+gR)+XX*gI;
6310 us = XX*(gc-gR)+XP*gI;
6315 um = rotx(Fp,uc,Fx,us,u);
6316 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6324 for(i=0; i<
nIFO; i++) {
6334 for(i=0; i<
nIFO; i++) {
6335 cc = local ? AM[
i]/(AM[
i]*AM[
i]+2.) : 1.;
6336 pp[
i] = cc*(AM[
i]-hh*u[
i])*u[i];
6337 qq[
i] = cc*((2.*hh*u[
i]-AM[
i])*v[i] + u[i]*u[i]*gg);
6338 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6341 cc = sqrt(si*si+co*co)+1.e-24;
6349 for(i=0; i<
nIFO; i++) {
6350 e[
i] = u[
i]*co + v[
i]*si;
6351 v[
i] = v[
i]*co - u[
i]*si;
6360 for(i=0; i<
nIFO; i++) {
6361 cc = local ? AM[
i]/(AM[
i]*AM[
i]+2.) : 1.;
6362 pp[
i] = cc*(AM[
i]-hh*u[
i])*u[i];
6363 qq[
i] = cc*((2.*hh*u[
i]-AM[
i])*v[i] + u[i]*u[i]*gg);
6364 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6367 cc = sqrt(si*si+co*co)+1.e-24;
6374 for(i=0; i<
nIFO; i++) {
6375 u[
i] = u[
i]*co + v[
i]*si;
6382 for(i=0; i<
nIFO; i++) {
6383 pix->
setdata(u[i]*hh/sqrt(nv[i][j]),
'U',i);
6384 WP += Fp[
i]*u[
i]*hh;
6385 WX += Fx[
i]*u[
i]*hh;
6393 gg = (gp+
gx)*ap + (gp-gx)*ac + 2*gI*as;
6394 hp = (wp*(ap+ac) + 2*AA*WX + wx*as)/gg;
6395 hx = (wx*(ap-ac) - 2*AA*WP + wp*as)/gg;
6396 HP = (WP*(ap+ac) - 2*AA*wx + WX*as)/gg;
6397 HX = (WX*(ap-ac) + 2*AA*wp + WP*as)/gg;
6398 La = rotx(Fp,hp,Fx,hx,e)+rotx(Fp,HP,Fx,HX,e);
6400 for(i=0; i<
nIFO; i++) {
6401 pix->
setdata((hp*fp[i][l]+hx*fx[i][l])/GG.
data[j],
'W',i);
6402 pix->
setdata((HP*fp[i][l]+HX*fx[i][l])/GG.
data[j],
'U',i);
6409 if(!core) pix->
core = (E00<Eo && E90<Eo) ?
false :
true;
6411 this->pixeLHood.data[pix->
time] = La/2.;
6412 this->pixeLNull.data[pix->
time] = (E00+E90-La)/2+1.;
6422 vtof->push_back(
int(M/2)-
int(m0[l])); ,
6423 vtof->push_back(
int(M/2)-
int(
m1[l])); ,
6424 vtof->push_back(
int(M/2)-
int(
m2[l])); ,
6425 vtof->push_back(
int(M/2)-
int(m3[l])); ,
6426 vtof->push_back(
int(M/2)-
int(m4[l])); ,
6427 vtof->push_back(
int(M/2)-
int(m5[l])); ,
6428 vtof->push_back(
int(M/2)-
int(m6[l])); ,
6429 vtof->push_back(
int(M/2)-
int(m7[l])); )
6434 T = cTo.
data[
k]+getifo(0)->TFmap.start();
6435 if(iID<=0 && type!=
'E') getSkyArea(
id,n,T);
6437 if(
fabs(Lm-Lo)/Lo>1.e-4)
6438 cout<<
"likelihood warning: "<<Lm<<
" "<<Lo<<endl;
6441 cout<<
"max value: "<<STAT<<
" at (theta,phi) = ("<<nLikelihood.getTheta(l)
6442 <<
","<<nLikelihood.getPhi(l)<<
") Likelihood: loop: "
6443 <<Lm<<
", final: "<<Lo<<
", sky: "<<LPm/2<<
", energy: "<<Em/2<<endl;
6447 if (mdcListSize() && n==0) {
6448 if (this->getwave(
id, 0,
'W')) {
6451 size_t M = this->ifoList.size();
6452 for(
int i=0; i<(
int)M; i++) {
6453 pd = this->getifo(i);
6454 pd->
RWFID.push_back(
id);
6458 double gps = this->getifo(0)->getTFmap()->start();
6466 pd->
RWFP.push_back(wf);
6472 this->nSensitivity.gps =
T;
6473 this->nAlignment.gps =
T;
6474 this->nDisbalance.gps =
T;
6475 this->nLikelihood.gps =
T;
6476 this->nNullEnergy.gps =
T;
6477 this->nCorrEnergy.gps =
T;
6478 this->nCorrelation.gps =
T;
6479 this->nSkyStat.gps =
T;
6480 this->nEllipticity.gps =
T;
6481 this->nPolarisation.gps =
T;
6482 this->nNetIndex.gps =
T;
6498 size_t N = this->ifoList.size();
6499 int N_1 = N>2 ?
int(N)-1 : 2;
6500 int N_2 = N>2 ?
int(N)-2 : 1;
6501 if(!N)
return false;
6506 vector<wavearray<double> >
snr(N);
6507 vector<wavearray<double> > nul(N);
6512 std::vector<int>* vi;
6513 std::vector<wavecomplex>
A;
6527 double um,vm,cc,hh,
gr,gc,gI,
Et,E,Xp,Xx,gp,
gx;
6528 double gg,co,si,uc,us,vc,vs,gR,
a,nr;
6536 double response = 0.;
6543 double GAMMA = 1.-gamma*gamma;
6544 bool status =
false;
6546 this->
gNET = this->aNET = this->eCOR = E = 0.;
6549 for(n=0; n<
N; n++) {
6550 nul[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'W',type);
6551 snr[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'S',type);
6552 this->getifo(n)->sSNR = 0.;
6553 this->getifo(n)->xSNR = 0.;
6554 this->getifo(n)->ekXk = 0.;
6555 this->getifo(n)->null = 0.;
6556 for(m=0; m<
NIFO; m++) { this->getifo(n)->ED[
m] = 0.; }
6557 for(m=0; m<
N; m++) { NDM[
n][
m] = 0.; }
6558 esnr[
n] = xsnr[
n] = 0.;
6561 if(!this->wc_List[lag].
size())
return status;
6563 cid = this->wc_List[lag].
get((
char*)
"ID",0,
'S',type);
6564 rat = this->wc_List[lag].
get((
char*)
"rate",0,
'S',type);
6565 lik = this->wc_List[lag].
get((
char*)
"like",0,
'S',type);
6568 for(k=0; k<
K; k++) {
6570 id = size_t(cid[k]+0.1);
6571 if(
id != ID)
continue;
6573 vi = &(this->wc_List[lag].cList[ID-1]);
6581 for(j=0; j<V; j++) {
6583 pix = this->wc_List[lag].getPixel(ID,j);
6585 cout<<
"network::setndm() error: NULL pointer"<<endl;
6588 if(!pix->
core && core)
continue;
6589 if(pix->
rate != rat.
data[k])
continue;
6592 gr=Et=gg=Xp=Xx = 0.;
6595 for(n=0; n<
N; n++) {
6601 for(n=0; n<
N; n++) {
6605 pd = this->getifo(n);
6609 A[
n].set(Fp[n],Fx[n]);
6611 gr += A[
n].abs()/2.;
6623 this->
gNET += (gr+gc)*Et;
6624 this->aNET += (gr-gc)*Et;
6633 um = vm = hh = cc = 0.;
6634 for(n=0; n<
N; n++) {
6635 u[
n] = Fp[
n]*uc + Fx[
n]*us;
6636 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
6640 cc += u[
n]*am[
n]*u[
n]*am[
n];
6644 if((hh*hh-cc)/um <= 0.)
continue;
6649 for(n=0; n<
N; n++) {
6650 if(u[n]*u[n]/um > 1-GAMMA) ii++;
6651 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
6653 this->iNET += ii*
Et;
6655 if(ii<N_2 && gamma<0.)
continue;
6658 uc = Xp*(gx+gg) - Xx*gI;
6659 us = Xx*(gp+gg) - Xp*gI;
6661 if(ii<N_1 && gamma!=0) {
6662 uc = Xp*(gc+gR)+Xx*gI;
6663 us = Xx*(gc-gR)+Xp*gI;
6669 for(n=0; n<
N; n++) {
6670 u[
n] = Fp[
n]*uc + Fx[
n]*us;
6671 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
6680 for(n=0; n<
N; n++) {
6690 for(n=0; n<
N; n++) {
6691 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
6692 qq[
n] = cc*((2*hh*u[
n]-am[
n])*v[n]+u[n]*u[n]*gg);
6693 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
6694 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
6697 cc = atan2(si,co+1.e-24);
6703 for(n=0; n<
N; n++) {
6704 e[
n] = u[
n]*co+v[
n]*si;
6705 r[
n] = v[
n]*co-u[
n]*si;
6713 for(n=0; n<
N; n++) {
6714 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
6715 pp[
n] = cc*(am[
n]-hh*e[
n])*e[n];
6716 qq[
n] = cc*((2*hh*e[
n]-am[
n])*r[n]+e[n]*e[n]*gg);
6717 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
6720 cc = atan2(si,co+1.e-24);
6725 for(n=0; n<
N; n++) {
6726 e[
n] = e[
n]*co+r[
n]*si;
6731 for(n=0; n<
N; n++) {
6733 this->getifo(n)->null += e[
n]*e[
n];
6735 for(m=0; m<
N; m++) {
6736 gg = e[
n]*am[
n]*e[
m]*am[
m];
6738 response += e[
n]*e[
m]*am[
m];
6739 if(n!=m) this->eCOR += gg;
6743 esnr[
n] += response*response;
6744 xsnr[
n] += am[
n]*response;
6745 e_SNR += response*response;
6746 x_SNR += am[
n]*response;
6747 this->getifo(n)->ED[3] +=
fabs(response*(am[n]-response));
6748 this->getifo(n)->ED[4] +=
fabs(response*(am[n]-response));
6754 this->norm = x_SNR/e_SNR;
6755 if(
fabs(this->norm-1) > 1.e-4)
6756 cout<<
"network::setndm(): incorrect likelihood normalization: "<<this->norm<<endl;
6758 for(n=0; n<
N; n++) {
6759 bIAS += this->getifo(n)->null;
6760 this->getifo(n)->null += nul[
n].data[
k];
6761 this->getifo(n)->sSNR = esnr[
n];
6762 this->getifo(n)->xSNR = xsnr[
n];
6763 S_NUL += this->getifo(n)->null;
6764 S_NIL += nul[
n].data[
k];
6765 S_SNR += snr[
n].data[
k];
6767 this->getifo(n)->ED[0] = (esnr[
n]-xsnr[
n]);
6768 this->getifo(n)->ED[1] =
fabs(esnr[n]-xsnr[n]);
6769 this->getifo(n)->ED[2] =
fabs(esnr[n]-xsnr[n]);
6771 for(m=0; m<
N; m++) S_NDM += NDM[n][m];
6775 if(count) { this->
gNET /= E; this->aNET /= E; this->iNET /= E;}
6777 a = S_NDM - lik.
data[
k];
6778 if(
fabs(a)/S_SNR>1.e-6)
6779 cout<<
"ndm-likelihood mismatch: "<<a<<
" "<<S_NDM<<
" "<<lik.
data[
k]<<
" "<<norm<<endl;
6780 a =
fabs(1 - (S_NDM+S_NIL)/S_SNR)/
count;
6782 cout<<
"biased energy disbalance: "<<a<<
" "<<S_SNR-S_NDM<<
" "<<S_NIL<<
" size="<<count<<endl;
6797 size_t N = this->ifoList.size();
6798 int N_1 = N>2 ?
int(N)-1 : 2;
6799 int N_2 = N>2 ?
int(N)-2 : 1;
6800 if(!N)
return false;
6805 vector<wavearray<double> >
snr(N);
6806 vector<wavearray<double> > nul(N);
6807 vector<wavearray<double> >
SNR(N);
6808 vector<wavearray<double> > NUL(N);
6813 std::vector<int>* vint;
6814 std::vector<wavecomplex>
A;
6834 double a,b,aa,
psi,gg,
gr,gc,gI,gR,E90,E00,E,fp,fx,vc,vs;
6835 double xx,xp,XX,XP,uc,us,co,hh,gp,
gx,xi00,xi90,o00,o90;
6836 double hgw,HGW,wp,WP,wx,WX,HH,um,vm,si,cc;
6847 double GAMMA = 1.-gamma*gamma;
6850 double nC = this->MRA ? 1. : 2.;
6852 bool status =
false;
6855 if(tYPe==
'I' || tYPe==
'S' || tYPe==
'G') ISG =
true;
6856 if(tYPe==
'i' || tYPe==
's' || tYPe==
'g') ISG =
true;
6866 for(n=0; n<
N; n++) {
6867 nul[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'W',type);
6868 NUL[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'U',type);
6869 snr[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'S',type);
6870 SNR[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'P',type);
6871 this->getifo(n)->sSNR = 0.;
6872 this->getifo(n)->xSNR = 0.;
6873 this->getifo(n)->ekXk = 0.;
6874 this->getifo(n)->null = 0.;
6875 for(m=0; m<5; m++) { this->getifo(n)->ED[
m] = 0.; }
6876 for(m=0; m<
N; m++) { NDM[
n][
m] = 0.; }
6877 esnr[
n]=xsnr[
n]=ssnr[
n]=SSNR[
n]=0.;
6880 if(!this->wc_List[lag].
size())
return false;
6882 cid = this->wc_List[lag].
get((
char*)
"ID",0,
'S',type);
6883 rat = this->wc_List[lag].
get((
char*)
"rate",0,
'S',type);
6884 lik = this->wc_List[lag].
get((
char*)
"like",0,
'S',type);
6887 for(k=0; k<
K; k++) {
6889 if(
size_t(cid[k]+0.1) != ID)
continue;
6891 vint = &(this->wc_List[lag].cList[ID-1]);
6899 for(j=0; j<V; j++) {
6901 pix = this->wc_List[lag].getPixel(ID,j);
6903 cout<<
"network::SETNDM() error: NULL pointer"<<endl;
6906 if(!pix->
core && core)
continue;
6907 if((pix->
rate != rat.
data[k]) && type)
continue;
6911 Z.
set(cos(psi),-sin(psi));
6914 gr=gg=xp=xx=XP=XX=E00=E90 = 0.;
6918 for(n=0; n<
N; n++) {
6924 for(n=0; n<
N; n++) {
6930 pd = this->getifo(n);
6933 A[
n].set(fp/b,fx/b);
6935 Fp[
n] = A[
n].real();
6936 Fx[
n] = A[
n].imag();
6937 gr += A[
n].abs()/2.;
6954 this->
gNET += (gr+gc)*(E00+E90);
6955 this->aNET += (gr-gc)*(E00+E90);
6968 um = vm = hh = cc = 0.;
6969 for(n=0; n<
N; n++) {
6970 u[
n] = Fp[
n]*uc + Fx[
n]*us;
6971 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
6975 cc += u[
n]*am[
n]*u[
n]*am[
n];
6979 if((hh*hh-cc)/um<=0. || E00<Eo) o00=0.;
6984 for(n=0; n<
N; n++) {
6985 if(u[n]*u[n]/um > 1-GAMMA) ii++;
6986 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
6988 this->iNET += ii*E00;
6989 if(ii<N_2 && gamma<0.) o00=0.;
6992 uc = xp*(gx+gg) - xx*gI;
6993 us = xx*(gp+gg) - xp*gI;
6995 if(ii<N_1 && gamma!=0) {
6996 uc = xp*(gc+gR)+xx*gI;
6997 us = xx*(gc-gR)+xp*gI;
7003 for(n=0; n<
N; n++) {
7004 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7005 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7014 for(n=0; n<
N; n++) {
7024 for(n=0; n<
N; n++) {
7025 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
7026 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
7027 qq[
n] = cc*((2.*hh*u[
n]-am[
n])*v[n] + u[n]*u[n]*gg);
7028 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7031 cc = sqrt(si*si+co*co)+1.e-24;
7039 for(n=0; n<
N; n++) {
7040 e[
n] = u[
n]*co + v[
n]*si;
7041 v[
n] = v[
n]*co - u[
n]*si;
7050 for(n=0; n<
N; n++) {
7051 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
7052 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
7053 qq[
n] = cc*((2.*hh*u[
n]-am[
n])*v[n] + u[n]*u[n]*gg);
7054 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7057 cc = sqrt(si*si+co*co)+1.e-24;
7063 for(n=0; n<
N; n++) {
7064 u[
n] = u[
n]*co+v[
n]*si;
7067 wx += Fx[
n]*u[
n]*aa;
7069 for(n=0; n<
N; n++) {
7070 h00[
n] = u[
n]*am[
n]*o00;
7072 am[
n] = hh*u[
n]*o00;
7086 um = vm = hh = cc = 0.;
7087 for(n=0; n<
N; n++) {
7088 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7089 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7093 cc += u[
n]*AM[
n]*u[
n]*AM[
n];
7097 if((hh*hh-cc)/um<=0. || E90<Eo) o90=0.;
7102 for(n=0; n<
N; n++) {
7103 if(u[n]*u[n]/um > 1-GAMMA) ii++;
7104 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
7106 this->iNET += ii*E90;
7107 if(ii<N_2 && gamma<0.) o90=0.;
7110 uc = XP*(gx+gg) - XX*gI;
7111 us = XX*(gp+gg) - XP*gI;
7113 if(ii<N_1 && gamma!=0) {
7114 uc = XP*(gc+gR)+XX*gI;
7115 us = XX*(gc-gR)+XP*gI;
7121 for(n=0; n<
N; n++) {
7122 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7123 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7132 for(n=0; n<
N; n++) {
7142 for(n=0; n<
N; n++) {
7143 cc = local ? AM[
n]/(AM[
n]*AM[
n]+2.) : 1.;
7144 pp[
n] = cc*(AM[
n]-hh*u[
n])*u[n];
7145 qq[
n] = cc*((2.*hh*u[
n]-AM[
n])*v[n] + u[n]*u[n]*gg);
7146 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7149 cc = sqrt(si*si+co*co)+1.e-24;
7157 for(n=0; n<
N; n++) {
7158 e[
n] = u[
n]*co + v[
n]*si;
7159 v[
n] = v[
n]*co - u[
n]*si;
7168 for(n=0; n<
N; n++) {
7169 cc = local ? AM[
n]/(AM[
n]*AM[
n]+2.) : 1.;
7170 pp[
n] = cc*(AM[
n]-hh*u[
n])*u[n];
7171 qq[
n] = cc*((2.*hh*u[
n]-AM[
n])*v[n] + u[n]*u[n]*gg);
7172 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7175 cc = sqrt(si*si+co*co)+1.e-24;
7181 for(n=0; n<
N; n++) {
7182 u[
n] = u[
n]*co+v[
n]*si;
7185 WX += Fx[
n]*u[
n]*aa;
7187 for(n=0; n<
N; n++) {
7188 h90[
n] = u[
n]*AM[
n]*o90;
7190 AM[
n] = HH*u[
n]*o90;
7194 cc = ISG ? (wp*WX-wx*WP)/gg : 0.0;
7195 hh = ISG ? (wp*wp+wx*wx)/gg/nC : 1./nC;
7196 HH = ISG ? (WP*WP+WX*WX)/gg/nC : 1./nC;
7204 for(n=0; n<
N; n++) {
7205 hgw += (am[
n]*Fp[
n] + aa*AM[
n]*Fx[
n])/gg;
7206 HGW += (AM[
n]*Fp[
n] - aa*am[
n]*Fx[
n])/gg;
7209 for(n=0; n<
N; n++) {
7213 for(m=0; m<
N; m++) {
7215 a = h00[
n]*h00[
m]*hh
7219 xi00 += u00[
n]*h00[
m];
7220 xi90 += u90[
n]*h90[
m];
7223 if(n!=m) this->eCOR +=
a;
7225 a = h00[
n]*h00[
m] - h90[
n]*h90[
m];
7226 this->getifo(n)->ED[3] +=
a;
7230 a = u00[
n]*u00[
n]*hh + u90[
n]*u90[
n]*HH;
7231 this->getifo(n)->null +=
a;
7234 if(ISG) xi00 = hgw*Fp[
n]-aa*HGW*Fx[
n];
7235 esnr[
n] += xi00*xi00;
7238 this->getifo(n)->ED[1] += xi00*(xx-xi00);
7239 this->getifo(n)->ED[4] +=
fabs(xi00*(xx-xi00));
7242 if(ISG) xi90 = HGW*Fp[
n]+aa*hgw*Fx[
n];
7243 esnr[
n] += xi90*xi90;
7246 this->getifo(n)->ED[2] += xi90*(XX-xi90);
7247 this->getifo(n)->ED[4] +=
fabs(xi90*(XX-xi90));
7254 for(n=0; n<
N; n++) {
7256 b = (SSNR[
n]+ssnr[
n] - esnr[
n])/nC;
7257 this->getifo(n)->null += b;
7258 this->getifo(n)->sSNR = esnr[
n]/nC;
7259 this->getifo(n)->xSNR = xsnr[
n]/nC;
7264 for(m=0; m<
N; m++) S_NDM += NDM[n][m];
7266 if(
fabs(ssnr[n]-snr[n].data[k])/ssnr[n] > 1.e-6 ||
fabs(SSNR[n]-SNR[n].data[k])/SSNR[n]>1.e-6)
7267 cout<<
"SETNDM()-likelihoodI() SNR mismatch: "<<n<<
" "<<ID<<
" "
7268 <<ssnr[
n]<<
":"<<snr[
n].data[
k]<<
" "<<SSNR[
n]<<
":"<<SNR[
n].data[
k]<<endl;
7271 if(count) { this->
gNET /= E; this->aNET /= E; this->iNET /= E; }
7273 a = S_NDM - lik.
data[
k];
7274 if(
fabs(a)/S_NDM>1.e-6)
7275 cout<<
"NDM-likelihood mismatch: "<<a<<
" "<<S_NDM<<
" "<<lik.
data[
k]<<endl;
7277 a =
fabs(1 - nC*(S_NDM+S_NUL)/(s_snr+S_SNR))/
count;
7279 cout<<
"biased energy disbalance: "<<a<<
" "<<S_NDM+S_NUL<<
" "<<(s_snr+S_SNR)/nC<<endl;
7292 const char*
fname,
const char* fmode,
size_t*
lagSite) {
7294 size_t nIFO = this->ifoList.size();
7295 this->wc_List.clear(); this->livTime.clear();
7298 cout <<
"network::setTimeShifts : lagStep must be positive" << endl;
7302 if(lagSize<1) lagSize=1;
7304 if(strcmp(fmode,
"r") && strcmp(fmode,
"w") && strcmp(fmode,
"s")) {
7305 cout <<
"network::setTimeShifts : bad fmode : must be r/w/s" << endl;
7309 if(fname) {
if(strlen(fname)<1) fname = NULL; }
7322 int maxIter = 10000000;
7325 for(n=0;n<
NIFO;n++) {
7337 lagIDS +=
int(getifo(0)->sHIFt/lagStep);
7338 maxList = lagSize+lagIDS;
7340 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7341 for(m=0; m<maxList; m++) {
7342 for(n=0; n<
nIFO; n++) {
7343 pd = this->getifo(n);
7344 lagList[
n][
m] = n==0 ? m :
int(pd->
sHIFt/lagStep);
7353 if(fname && (!strcmp(fmode,
"r") || !strcmp(fmode,
"s"))) {
7354 if(!strcmp(fmode,
"r")) {
7359 cout <<
"network::setTimeShifts : Error Opening File : " << fname << endl;
7367 in.getline(str,1024);
7368 if (!in.good())
break;
7369 if(str[0] !=
'#') maxList++;
7372 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7373 in.clear(ios::goodbit);
7374 in.seekg(0, ios::beg);
7377 in.getline(str,1024);
7378 if(str[0] ==
'#')
continue;
7379 in.seekg(fpos, ios::beg);
7382 for(n=0; n<nIFO; n++) in >> lagList[
n][
m];
7383 if (!in.good())
break;
7389 if(!strcmp(fmode,
"s")) {
7398 in.getline(str,1024);
7399 if (!in.good())
break;
7400 if(str[0] !=
'#') maxList++;
7403 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7404 in.clear(ios::goodbit);
7405 in.seekg(0, ios::beg);
7408 in.getline(str,1024);
7409 if(str[0] ==
'#')
continue;
7410 in.seekg(fpos, ios::beg);
7413 for(n=0; n<nIFO; n++) in >> lagList[
n][
m];
7414 if (!in.good())
break;
7422 for(m=0; m<maxList; m++){
7424 for (n=0; n<
nIFO; n++)
id[n]=lagList[n][m];
7428 for (n=0; n<nIFO; n++) if(id[n]<0||id[n]>
int(lagMax)) check=
false;
7432 for (
int i=nIFO-1;
i>=0;
i--) {
7433 for (
int j=
i-1;
j>=0;
j--) {
7434 if (!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7435 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7442 cout <<
"network::setTimeShifts : no lags in the list" << endl;
7443 cout <<
"lagP : " << lagP <<
" " << lagSize << endl;
7446 if(lagP!=
int(maxList)) {
7447 cout <<
"network::setTimeShifts : lags out of lagMax" << endl;
7448 cout <<
"lagP : " << lagP <<
" " << lagSize << endl;
7457 if(lagSite!=NULL)
for(n=0; n<
nIFO; n++) {
7458 if(lagSite[n] >= nIFO) {
7459 cout <<
"network::setTimeShifts : Error lagSite - value out of range " << endl;
7464 for(n=1; n<
nIFO; n++) N[n]=lagMax;
7468 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7469 for(n=0; n<
nIFO; n++) lagList[n][nList]=0;
7473 for (
int k=0;k<maxIter;k++) {
7474 for(n=0; n<
nIFO; n++) ID[n] = TMath::Nint(rnd.Uniform(-(N[n]+0.5),N[n]+0.5));
7475 for(n=0; n<
nIFO; n++)
id[n] = (lagSite==NULL) ? ID[
n] : ID[lagSite[
n]];
7477 for(
int i=nIFO-1;
i>=0;
i--) {
7478 for(
int j=
i-1;
j>=0;
j--) {
7479 if(!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7480 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7482 if(lagSite[
i]!=lagSite[j] &&
id[
i]==
id[j]) check=
false;
7484 if(
id[
i]==
id[j]) check=
false;
7490 for(m=0;m<nList;m++) {
7492 for(n=0; n<
nIFO; n++)
if(lagList[n][m]!=
id[n]) pass=
false;
7493 if(pass) check=
false;
7497 if(NETX(
id[0]||,
id[1]||,
id[2]||,
id[3]||,
id[4]||,
id[5]||,
id[6]||,
id[7]||)
false) {
7498 for(n=0; n<
nIFO; n++) lagList[n][nList]=
id[n];
7502 if (nList>=maxList)
break;
7510 for(m=0; m<nList; m++) {
7511 int lagMin = kMaxInt;
7512 for(n=0; n<
nIFO; n++)
if (lagList[n][m]<lagMin) lagMin=lagList[
n][
m];
7513 for(n=0; n<
nIFO; n++) lagList[n][m]-=lagMin;
7516 if(lagIDS+lagSize>nList) {
7517 cout <<
"network::setTimeShifts : lagOff+lagSize > nList of lags : " << nList << endl;
7521 for(n=0; n<
nIFO; n++){
7522 pd = this->getifo(n);
7531 double R = this->getifo(0)->getTFmap()->
rate();
7532 double segLen = this->getifo(0)->getTFmap()->size();
7533 double edge = this->Edge;
7540 segLen = (segLen/R-2*edge)/lagStep;
7541 lagMaxSeg =
int(segLen)-1;
7543 for(n=0; n<
nIFO; n++) {
7545 lagH[
n] = lagMaxSeg;
7550 for (n=0; n<
nIFO; n++)
id[n]=lagList[n][m+lagIDS];
7553 for(n=0; n<nIFO; n++) if(id[n]<0||id[n]>
int(lagMaxSeg)) check=
false;
7556 for(
int i=nIFO-1;
i>=0;
i--) {
7557 for(
int j=
i-1;
j>=0;
j--) {
7558 if (!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7559 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7567 for(n=0; n<
nIFO; n++) {
7568 k = lagList[
n][m+lagIDS];
7569 if(k) check =
false;
7570 this->getifo(n)->lagShift.data[selSize] = k*
lagStep;
7574 k = lagList[0][m+lagIDS];
7575 this->getifo(0)->lagShift.data[selSize] = k*
lagStep;
7577 for(n=1; n<
nIFO; n++) {
7578 pd = this->getifo(n);
7582 if(zero>0.1) check =
false;
7585 wc_List.push_back(wc);
7586 livTime.push_back(0.);
7592 cout <<
"network::setTimeShifts error: no lag was selected" << endl;
7596 for(n=0; n<
nIFO; n++) {
7597 m = this->getifo(n)->lagShift.size();
7598 if(m!=selSize) this->getifo(n)->lagShift.resize(selSize);
7603 if(fname && !strcmp(fmode,
"w") && lagMax) {
7606 if((fP = fopen(fname,
"w")) == NULL) {
7607 cout <<
"network::setTimeShifts error: cannot open file " << fname << endl;
7613 fprintf(fP,
"#total %10d lags \n",
int(nList));
7615 fprintf(fP,
"#%13s%14s%14s\n",
" nIFO",
"lagStep",
" lagMax");
7616 fprintf(fP,
"#%13d%14.3f%14d\n",
int(nIFO),lagStep,
int(lagMax));
7619 for(n=0; n<
nIFO; n++)
fprintf(fP,
"%12s-%1d",
"lagShift",
int(n));
7624 for(m=0; m<nList; m++){
7626 for (n=0; n<
nIFO; n++)
fprintf(fP,
"%14d",lagList[n][m]);
7635 for(n=0; n<
nIFO; n++)
delete [] lagList[n];
7640 for(n=0; n<
nIFO; n++)
printf(
"%12.12s%2s",
"ifo",getifo(n)->Name);
7642 for(m=0; m<selSize; m++){
7644 for(n=0; n<
nIFO; n++)
printf(
"%14.5f",this->getifo(n)->lagShift.data[m]);
7660 std::vector<delayFilter>().swap(this->
filter);
7662 if(!d) d = this->getifo(0);
7665 if(ifoList.size()<2 || !d || rate==0.) {
7666 cout<<
"network::setFilter() error: incomplete network initialization"<<endl;
7670 double T = this->getDelay((
char*)
"MAX")+0.002;
7677 int J =
int((
fabs(T)*rate*N)/M+0.5);
7679 if(N<M) { cout<<
"network::setFilter() error"<<endl;
return 0; }
7682 this->getifo(0)->nDFS =
N;
7683 this->
filter.reserve((2*J-1)*M);
7685 for(i=0; i<
M; i++) {
7686 for(j=-(J-1); j<J; j++) {
7687 m = j>0 ? (j+N/2-1)/N : (j-N/2)/
N;
7692 for(k=0; k<
K; k++) v.
index[k] -= m*M;
7693 this->filter.push_back(v);
7708 size_t N = this->ifoList.size();
7710 if(d) this->getifo(0)->setFilter(*d);
7712 this->getifo(0)->clearFilter();
7723 size_t N = this->ifoList.size();
7727 this->getifo(0)->readFilter(gname);
7729 this->filter90.clear();
7730 std::vector<delayFilter>().swap(this->filter90);
7731 this->filter90 = this->
filter;
7733 this->getifo(0)->readFilter(fname);
7735 this->getifo(0)->clearFilter();
7745 size_t N = this->ifoList.size();
7748 this->readFilter(gname);
7749 this->filter90 = this->
filter;
7751 this->readFilter(fname);
7763 if ( (fp=fopen(fname,
"wb")) == NULL ) {
7764 cout <<
" network::writeFilter() error : cannot open file " << fname <<
". \n";
7768 size_t M = size_t(getifo(0)->TFmap.maxLayer()+1);
7769 size_t N = size_t(
filter.size()/
M);
7771 size_t n = K *
sizeof(float);
7772 size_t m = K *
sizeof(short);
7777 fwrite(&K,
sizeof(
size_t), 1, fp);
7778 fwrite(&M,
sizeof(
size_t), 1, fp);
7779 fwrite(&N,
sizeof(
size_t), 1, fp);
7781 for(i=0; i<
M; i++) {
7782 for(j=0; j<
N; j++) {
7783 for(k=0; k<
K; k++) {
7787 fwrite(value.data, n, 1, fp);
7788 fwrite(index.
data, m, 1, fp);
7802 if ( (fp=fopen(fname,
"rb")) == NULL ) {
7803 cout <<
" network::readFilter() error : cannot open file " << fname <<
". \n";
7811 fread(&K,
sizeof(
size_t), 1, fp);
7812 fread(&M,
sizeof(
size_t), 1, fp);
7813 fread(&N,
sizeof(
size_t), 1, fp);
7815 size_t n = K *
sizeof(float);
7816 size_t m = K *
sizeof(short);
7825 for(k=0; k<
K; k++) {
7826 v.
value.push_back(0.);
7827 v.
index.push_back(0);
7830 for(i=0; i<
M; i++) {
7831 for(j=0; j<
N; j++) {
7832 fread(value.data, n, 1, fp);
7833 fread(index.
data, m, 1, fp);
7834 for(k=0; k<
K; k++) {
7855 int nres = this->wdmList.size();
7856 int nIFO = this->ifoList.size();
7857 int V =
int(this->pList.size());
7865 for(
int i=0;
i<nres;
i++) {
7867 pd = this->getifo(
k);
7868 pd->
vSS[
i].Expand(
false);
7871 pd->
vSS[
i].Inverse();
7873 pd->
vSS[
i].getLayer(x,0);
7877 pd->
vSS[
i].Forward(x);
7878 layers = pd->
vSS[
i].maxLayer()+1;
7880 for(
int j=0;
j<V;
j++) {
7881 pix = this->pList[
j];
7882 if(pix->
layers != layers)
continue;
7884 v00[
k][
j] = pd->
vSS[
i].GetMap00(ind);
7885 v90[
k][
j] = pd->
vSS[
i].GetMap90(ind);
7888 pd->
vSS[
i].Shrink();
7901 size_t N = this->ifoList.size();
7902 size_t k = this->getIndex(theta,phi);
7905 for(
size_t n=1;
n<
N;
n++){
7906 d = this->getifo(
n);
7924 size_t M = this->
filter.size()/I;
7925 size_t K = this->
filter[0].index.size();
7926 size_t jB = size_t(this->Edge*R/I)*I;
7933 double*
F = (
double*)malloc(K*
sizeof(
double));
7934 int* J = (
int*)malloc(K*
sizeof(
int));
7944 double* b0 = temp.
data;
7946 for(i=0; i<I; i++) {
7951 F[
k] = double(pv->
value[k]);
7958 for(j=jS; j<
N; j+=I) {
7978 size_t N = ifoList.size();
7985 cout<<
"network::setDelayIndex(): invalid network\n";
7990 for(n=0; n<
N; n++) dr[n] = ifoList[n];
7992 size_t I = dr[0]->
nDFL;
7993 size_t K = this->
filter.size()/I;
7997 rate *= dr[0]->
nDFS/I;
7999 if(pOUT) cout<<
"filter size="<<this->
filter.size()
8000 <<
" layers="<<I<<
" delays="<<K<<
" samples="<<dr[0]->
nDFS<<endl;
8002 if(!(K&1) || rate == 0.) {
8003 cout<<
"network::setDelayIndex(): invalid network\n";
8007 for(n=0; n<
N; n++) {
8018 this->nPenalty = dr[0]->
tau;
8027 for(n=0; n<
N; n++) {
8028 for(m=0; m<
N; m++) {
8030 i =
int(t*rate+2*K+0.5) - 2*
K;
8036 for(n=0; n<
N; n++) {
8038 for(m=0; m<
N; m++) {
8039 for(k=0; k<
N; k++) {
8040 t =
fabs(mm[n][k]-mm[n][m]-tt[m][k]);
8041 if(TT[n] < t) TT[
n] =
t;
8047 for(m=0; m<
N; m++) {
8048 if(t>TT[m]) { t = TT[
m]; k =
m; }
8050 this->nPenalty.set(l,
double(t));
8053 i = mIFO<9 ? mm[
k][this->mIFO] :
int(t*rate+2*K+0.5)-2*
K;
8059 for(m=0; m<
N; m++) {
8060 ii = (K/2+mm[
k][
m])-i;
8062 if(ii < 0) cout<<
"network::setDelayIndex error: sky index<0: "<<k<<endl;
8076 if(ifoList.size()<2 || !dr->
tau.
size()) {
8077 cout<<
"network::setIndex() - invalid network"<<endl;
8083 size_t N = ifoList.size();
8084 size_t I = mode ? dr->
nDFL : 1;
8087 size_t K = this->
filter.size()/I;
8089 size_t M = mIFO<9 ? mIFO : 0;
8093 if(this->skyMask.size()!=
L) this->skyMask.resize(L);
8094 if(this->skyHole.size()!=
L) { this->skyHole.resize(L); this->skyHole = 1.; }
8095 for(j=0; j<
L; j++) {
8097 skyMask.data[
j] = size_t(skyHole.data[j]+0.1);
8101 if(mode==2 || mode==4) {
8107 long long **pp = (
long long**)malloc(L*
sizeof(
long long*));
8110 for(n=0; n<
N; n++) {
8111 if(!this->getifo(n)->index.size()) {
8112 cout<<
"network::setIndex() - invalid network"<<endl;
8121 for(n=0; n<
N; n++) {
8122 if(n == M)
continue;
8123 ll = this->getifo(n)->index.data[
i];
8124 if(this->mIFO==99) ll += K/2 - this->getifo(0)->index.data[
i];
8125 delay.
data[
i] += ll<<(m*12);
8132 for(i=1; i<
L; i++) {
8133 j = pp[
i] - delay.
data;
8134 if(ll == *(pp[i])) {
8135 skyMask.data[
j] = 0;
8141 if(pOUT) cout<<
"\n ll="<<ll<<
" "<<j<<
"|"<<sm->
getTheta(j)<<
"|"<<sm->
getPhi(j);
8153 int nIFO = ifoListSize();
8154 for(
int n=0;
n<
nIFO;
n++) getifo(
n)->print();
8158 cout <<
"----------------------------------------------" << endl;
8159 cout <<
" INJECTIONS : " << this->mdcListSize() << endl;
8160 cout <<
"----------------------------------------------" << endl;
8161 for(
size_t k=0;
k<this->mdcListSize();
k++) {
8162 string str(this->getmdcList(
k));
8163 cout << endl << str.c_str() << endl;
wavearray< double > t(hp.size())
std::vector< char * > ifoName
static float _sse_abs_ps(__m128 *_a)
virtual void resize(unsigned int)
virtual size_t size() const
double getWFfreq(char atype='S')
std::vector< vector_int > cRate
std::vector< netcluster > wc_List
wavearray< double > getMRAwave(network *net, int ID, size_t n, char atype='S', int mode=0)
size_t add(detector *)
param: detector structure return number of detectors in the network
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 readMDClog(char *, double=0., int=11, int=12)
param: MDC log file param: approximate gps time
static float _avx_dpf_ps(double **Fp, double **Fx, int l, std::vector< float * > &pAPN, std::vector< float * > &pAVX, int I)
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
virtual void rate(double r)
std::vector< wavearray< float > > tdAmp
bool getwave(size_t, size_t, char='W')
param: cluster ID param: delay index param: time series type return: true if time series are extracte...
wavearray< double > a(hp.size())
double getWFtime(char atype='S')
long getNetworkPixels(int LAG, double Eo, double DD=1., TH1F *hist=NULL)
static void _sse_add4_ps(__m128 *_a, __m128 *_b, __m128 _c)
static __m128 _sse_abs4_ps(__m128 *_p)
static void _sse_zero_ps(__m128 *_p)
static __m128 _sse_ed4_ps(__m128 *_p, __m128 *_q, __m128 _L)
std::vector< vectorD > NDM
void getSkyArea(size_t id, size_t lag, double T)
param: cluster id param: lag param: cluster time
std::vector< int > neighbors
double iGamma(double r, double p)
size_t setIndexMode(size_t=0)
static __m128 _sse_dot4_ps(__m128 *_p, __m128 *_q)
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< pixdata > data
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
std::vector< delayFilter > filter
std::vector< wavearray< double > * > RWFP
static void _sse_cpf4_ps(__m128 *_aa, __m128 *_pp)
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
void add(const wavearray< DataType_t > &, int=0, int=0, int=0)
std::vector< vector_float > sArea
virtual size_t initwc(double, double)
param: cluster start time relative to segment start param: cluster duration return cluster list size ...
size_t setSkyMask(double f, char *fname)
std::vector< std::string > mdcType
double getTheta(size_t i)
std::slice getSlice(double n)
double getThetaStep(size_t i)
static void _sse_cpf_ps(float *a, __m128 *_p)
WSeries< double > waveBand
watplot p1(const_cast< char * >("TFMap1"))
size_t setFilter(detector *=NULL)
param: detector
double getPhiStep(size_t i)
std::vector< vector_int > cList
virtual void start(double s)
std::vector< size_t > mdc__ID
void downsample(wavearray< short > &, size_t=4)
static void _sse_rotm_ps(__m128 *u, float *c, __m128 *v, float *s, __m128 *a)
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
void updateTDamp(int, float **, float **)
std::vector< double > mdcTime
long subNetCut(int lag, float subnet=0.6, float subcut=0.33, TH2F *hist=NULL)
long coherence(double, double=0., double=0.)
param: threshold on lognormal pixel energy (in units of noise rms) param: threshold on total pixel en...
std::vector< double > vectorD
long likelihood2G(char mode, int lag, int ID, TH2F *hist=NULL)
std::vector< detector * > ifoList
std::vector< vector_float > p_Map
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
void writeFilter(const char *fname)
static float _avx_loadata_ps(float **p, float **q, float **u, float **v, float En, std::vector< float * > &pAVX, int I)
size_t wavecount(double x, int n=0)
long likelihoodWP(char mode, int lag, int ID, TH2F *hist=NULL)
static __m256 _avx_noise_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
double setVeto(double=5.)
param: time window around injections
int _sse_MRA_ps(network *net, float *amp, float *AMP, float Eo, int K)
void set(double x, double y)
std::vector< double > livTime
static __m256 _avx_stat_ps(float **x, float **X, float **s, float **S, std::vector< float * > &pAVX, int I)
static __m128 _sse_ecoh4_ps(__m128 *_p, __m128 *_q, __m128 _L)
wavearray< double > * getHoT()
param: no parameters
gwavearray< double > * gx
void setDelayIndex(double rate)
param: MDC log file
std::vector< std::string > mdcList
double getDelay(const char *c="")
static __m128 _sse_ind4_ps(__m128 *_p, __m128 _L)
WSeries< double > pTF[nRES]
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
static void _avx_free_ps(std::vector< float * > &v)
std::vector< vector_int > p_Ind
static void _sse_pol4_ps(__m128 *_fp, __m128 *_fx, __m128 *_v, double *r, double *a)
static float _avx_setAMP_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
static void _avx_loadNULL_ps(float **n, float **N, float **d, float **D, float **h, float **H, int I)
std::vector< waveSegment > segList
int _sse_mra_ps(network *NET, float *amp, float *AMP, float Eo, int K)
static void _sse_rot4m_ps(__m128 *_u, __m128 *_c, __m128 *_v, __m128 *_s, __m128 *_a)
double THRESHOLD(double bpp)
param: selected fraction of LTF pixels assuming Gaussian noise
double mchirp(int ID, double=2.5, double=1.e20, double=0)
void readFilter(const char *)
void setDelay(const char *="L1")
static float _avx_ort_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
static float _avx_GW_ps(float **p, float **q, std::vector< float * > &pAPN, float *rr, std::vector< float * > &pAVX, int II)
wavearray< double > lagShift
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
static void _avx_cpf_ps(float **p, float **q, float **u, float **v, int I)
long likelihoodI(char='P', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
double phi2RA(double ph, double gps)
netpixel * getPixel(size_t n, size_t i)
virtual void waveSort(DataType_t **pp, size_t l=0, size_t r=0) const
WSeries< double > waveForm
watplot p2(const_cast< char * >("TFMap2"))
WSeries< double > * getTFmap()
param: no parameters
virtual void delay(double T)
virtual void stop(double s)
double fabs(const Complex &x)
std::vector< short > index
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Meyer< double > S(1024, 2)
static __m128 _sse_sum_ps(__m128 **_p)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
std::vector< clusterdata > cData
static void _sse_point_ps(__m128 **_p, float **p, short **m, int l, int n)
void setDelayFilters(detector *=NULL)
static void _sse_mul_ps(__m128 *_a, float b)
double getwave(int, netcluster &, char, size_t)
param: no parameters
double iGamma1G(double r, double p)
double get(size_t i)
param: sky index
WaveDWT< DataType_t > * pWavelet
size_t netcut(double, char='L', size_t=0, int=1)
param: threshold param: minimum cluster size processed by the corrcut param: cluster type return: num...
void delay(double theta, double phi)
long likelihoodB(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
param: maximized statistic: param: threshold to define core pixels (in units of noise rms) ...
static float _sse_maxE_ps(__m128 *_a, __m128 *_A)
cout<< "live time after cat 2 : "<< detSegs_ctime<< endl;if(detSegs_ctime< segTHR){cout<< "job segment live time after cat2 < "<< segTHR<< " sec, job terminated !!!"<< endl;exit(1);}double Tb=detSegs[0].start;double Te=detSegs[0].stop;double dT=Te-Tb;char file[512], tdf00[512], tdf90[512], buFFer[1024];int rnID=int(gRandom->Rndm(13)*1.e9);if(simulation){i=NET.readMDClog(injectionList, double(long(Tb))-mdcShift);printf("GPS: %16.6f saved, injections: %d\n", double(long(Tb)), i);frTB[nIFO].shiftBurstMDCLog(NET.mdcList, ifos, mdcShift);for(int i=0;i< NET.mdcTime.size();i++) NET.mdcTime[i]+=mdcShift;vector< waveSegment > mdcSegs(NET.mdcTime.size());for(int k=0;k< NET.mdcTime.size();k++){mdcSegs[k].start=NET.mdcTime[k]-gap;mdcSegs[k].stop=NET.mdcTime[k]+gap;}vector< waveSegment > mdcSegs_dq2=slagTB.mergeSegLists(detSegs_dq2, mdcSegs);double mdcSegs_ctime=slagTB.getTimeSegList(mdcSegs_dq2);cout<< "live time in zero lag after cat2+inj : "<< mdcSegs_ctime<< endl;if(mdcSegs_ctime==0){cout<< "job segment with zero cat2+inj live time in zero lag, job terminated !!!"<< endl;exit(1);}}if(dump_infos_and_exit) exit(0);if(mask >0.) NET.setSkyMask(mask, skyMaskFile);for(i=0;i< nIFO;i++){frTB[i].readFrames(FRF[i], channelNamesRaw[i], x);x.start(x.start()+dataShift[i]);x.start(x.start()-segLen *(segID[i]-segID[0]));if(singleDetector) TB.resampleToPowerOfTwo(x);sprintf(file,"%s/%s_%d_%s_%d_%d.dat", nodedir, ifo[i], int(Tb), data_label, runID, rnID);if(dump_sensitivity_and_exit){sprintf(file,"%s/sensitivity_%s_%d_%s_job%d.txt", dump_dir, ifo[i], int(Tb), data_label, runID);cout<< endl<< "Dump Sensitivity : "<< file<< endl<< endl;TB.makeSpectrum(file, x);continue;}if(dcCal[i]>0.) x *=dcCal[i];if(fResample >0){x.FFT(1);x.resize(fResample/x.rate()*x.size());x.FFT(-1);x.rate(fResample);}pTF[i]=pD[i]-> getTFmap()
std::vector< SSeries< double > > vSS
static void _sse_minSNE_ps(__m128 *_pE, __m128 **_pe, __m128 *_es)
network & operator=(const network &)
WSeries< double > waveNull
double getdata(char type='R', size_t n=0)
size_t setRank(double, double=0.)
static __m128 _sse_like4_ps(__m128 *_f, __m128 *_a, __m128 *_A)
static void _sse_rot4p_ps(__m128 *_u, __m128 *_c, __m128 *_v, __m128 *_s, __m128 *_a)
std::vector< float > value
static void _sse_rotp_ps(__m128 *u, float *c, __m128 *v, float *s, __m128 *a)
virtual void resize(unsigned int)
size_t readSEGlist(char *, int=1)
param: segment list file param: start time collumn number
int setTimeShifts(size_t=1, double=1., size_t=0, size_t=0, const char *=NULL, const char *="w", size_t *=NULL)
param number of time lags param time shift step in seconds param first lag ID param maximum lag ID pa...
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
virtual void waveSplit(DataType_t **pp, size_t l, size_t r, size_t m) const
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
static float _avx_packet_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
bool setndm(size_t, size_t, bool=true, int=1)
param: cluster ID param: lag index param: statistic identificator param: resolution idenificator retu...
static void _sse_dpf4_ps(__m128 *_Fp, __m128 *_Fx, __m128 *_fp, __m128 *_fx)
bool SETNDM(size_t, size_t, bool=true, int=1)
static void _sse_pnp4_ps(__m128 *_fp, __m128 *_fx, __m128 *_am, __m128 *_AM, __m128 *_u, __m128 *_v)
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...
static void _sse_ort4_ps(__m128 *_u, __m128 *_v, __m128 *_s, __m128 *_c)
std::vector< vector_int > nTofF
void setSkyMaps(double, double=0., double=180., double=0., double=360.)
param: sky map granularity step, degrees param: theta begin, degrees param: theta end...
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
bool setdata(double a, char type='R', size_t n=0)
double threshold(double, double)
param: selected fraction of LTF pixels assuming Gaussian noise param: maximum time delay between dete...