29 #include <xmmintrin.h> 30 #include <immintrin.h> 47 nRun(0), nLag(1),
nSky(0), mIFO(0), rTDF(0), Step(0.), Edge(0.),
gNET(0.), aNET(0.), iNET(0),
48 eCOR(0.), norm(1.), e2or(0.),
acor(sqrt(2.)), pOUT(false),
EFEC(true), local(true),
optim(true),
52 this->ifoList.clear();
53 this->ifoName.clear();
55 this->mdcList.clear();
56 this->livTime.clear();
57 this->mdcTime.clear();
58 this->mdcType.clear();
59 this->mdc__ID.clear();
87 size_t nIFO = this->ifoList.size();
90 cout<<
"network::getNetworkPixels(): " 91 <<
"invalid number of detectors or\n";
94 if(getifo(0)->
getTFmap()->w_mode != 1) {
95 cout<<
"network::getNetworkPixels(): invalid whitening mode.\n";
103 int i,
j,
k,
m,
n,NN,jj,nM,jE,jb,je,J,
K;
107 double R = pTF->
wrate();
108 double r = hTS->
rate();
113 int jB =
int(this->Edge*R+0.001);
114 if(jB&1) {cout<<
"getNetworkPixels(1): WDM parity violation\n";
exit(1);}
117 cout<<
"network::getNetworkPixels(): insufficient data edge length.\n";
132 for(n=0; n<
nIFO; n++) {
133 pdata[
n] = getifo(n)->getTFmap()->data;
138 double a,b,E,Ct,Cb,Ht,Hb;
140 if(hist) {pixeLHood = *
pTF; pixeLHood=-1.;}
141 if(this->veto.size() !=
M) {
144 short* pveto = this->veto.data;
146 this->wc_List[LAG].clear();
147 this->livTime[LAG] = 0.;
148 this->wc_List[LAG].setlow(pTF->
getlow());
149 this->wc_List[LAG].sethigh(pTF->
gethigh());
152 for(n=0; n<
nIFO; n++) {
153 b = this->getifo(n)->lagShift.data[LAG];
154 if(a>b) { a = b; nM =
n; }
157 for(n=0; n<
nIFO; n++) {
158 b = this->getifo(n)->lagShift.data[LAG];
159 K =
int((b-a)*R+0.001);
160 if(K&1) {cout<<
"getNetworkPixels(2): WDM parity violation\n";
exit(1);}
161 in[
n] = IN[
n] = K+jB;
176 if(jE&1) {cout<<
"getNetworkPixels(3): WDM parity violation\n";
exit(1);}
181 for(jj=0; jj<NN; jj++) {
184 pmap = MAP.
data+(jj+jB)*I;
185 for(n=0; n<
nIFO; n++) {
186 if(in[n] >= jE) in[
n] -= NN;
187 jb =
int(in[n]*r/R+0.01);
188 je =
int((in[n]+1)*r/R+0.01);
189 while(jb<je)
if(!pveto[jb++]) VETO=0.;
190 PDATA = &(pdata[
n][in[
n]*I]);
191 for(i=0; i<I; i++) pmap[i]+=*PDATA++;
197 if(pmap[i]<Eo || i<ib) pmap[
i]=0.;
198 if(pmap[i]>Em) pmap[
i]=Em+0.1;
203 for(jj=0; jj<NN; jj++) {
205 pmap = MAP.
data+(jj+jB)*I;
206 for(n=0; n<
nIFO; n++) {
207 if(IN[n] >= jE) IN[
n] -= NN;
209 for(n=0; n<5; n++) pp[n]=pmap+(n-2)*I;
210 for(i=ib; i<ie; i++) {
211 if((E=pp[2][i])<Eo)
continue;
212 Ct = pp[2][i+1]+pp[3][
i ]+pp[3][i+1];
213 Cb = pp[2][i-1]+pp[1][
i ]+pp[1][i-1];
215 Ht+= i<II? pp[4][i+2]+pp[3][i+2]:0.;
217 Hb+= i>1 ? pp[0][i-2]+pp[1][i-2]:0.;
225 for(n=0; n<
nIFO; n++) {
228 pix.
data[
n].asnr = sqrt(pdata[n][j]);
232 if(hist) hist->Fill(E);
233 if(hist) pixeLHood.data[
j] = E;
238 wc_List[LAG].append(pix);
241 for(n=0; n<
nIFO; n++) IN[n]++;
245 this->wc_List[LAG].start = pTF->
start();
246 this->wc_List[LAG].stop = pTF->
stop();
247 this->wc_List[LAG].rate = pTF->
rate();
248 this->livTime[LAG] = count/R;
250 if(nPix) this->setRMS();
266 __m128* _pa = (__m128*) pa;
267 __m128* _pb = (__m128*) pb;
268 __m128* _px = (__m128*) px;
270 for(
int i=0;
i<n*4;
i++) {
273 x.
data[
i]=
i*float(m)/float(n*4);
276 for(
int i=0;
i<
n;
i++) {
277 *(_pa+
i) = _mm_sqrt_ps(*(_px+
i));
278 *(_pa+
i) = _mm_div_ps(_mm_set1_ps(1.),*(_pa+
i));
279 *(_pb+
i) = _mm_rsqrt_ps(*(_px+
i));
280 _mm_storeu_ps(qq,*(_px+
i));
281 _mm_storeu_ps(sq,*(_pa+
i));
282 _mm_storeu_ps(rq,*(_pb+
i));
283 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[0],sq[0],rq[0]);
284 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[1],sq[1],rq[1]);
285 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[2],sq[2],rq[2]);
286 printf(
"%d %9.7f %9.7f %9.7f \n",
i,qq[3],sq[3],rq[3]);
313 if(!this->wc_List[lag].
size())
return 0;
319 bool cirwave = mode==
'g' || mode==
'G' || mode==
'c' || mode==
'C';
320 bool linwave = mode==
'l' || mode==
'L' || mode==
's' || mode==
'S';
321 bool iotwave = mode==
'i' || mode==
'l' || mode==
'e' || mode==
'c' ||
322 mode==
'I' || mode==
'L' || mode==
'E' || mode==
'C';
323 bool psiwave = mode==
'l' || mode==
'e' || mode==
'p' ||
324 mode==
'L' || mode==
'E' || mode==
'P';
325 bool mureana = mode==
'i' || mode==
'e' || mode==
'c' ||
326 mode==
'r' || mode==
'p' || mode==
'b' ||
327 mode==
'l' || mode==
's' || mode==
'g';
328 bool rndwave = mode==
'r' || mode==
'R';
330 bool prior = this->gamma<0 ? true :
false;
331 bool m_chirp = this->
optim ? false : mureana;
333 if(!this->
optim) mureana =
true;
335 size_t nIFO = this->ifoList.size();
336 size_t ID = abs(iID);
339 cout<<
"network::likelihoodAVX(): invalid network.\n";
344 float gama = this->gamma*this->gamma*2./3.;
345 float deta =
fabs(this->
delta);
if(deta>1) deta=1;
346 float REG[2]; REG[0] = deta*sqrt(2);
349 static const __m128 _oo = _mm_set1_ps(1.
e-16);
350 static const __m128 _sm = _mm_set1_ps(-0.
f);
351 static const __m128 _En = _mm_set1_ps(En);
353 float aa,AA,Lo,Eo,Co,No,Ep,Lp,Np,Cp,Ec,Dc,To,Fo,Em,Lm,Rc,Mo,Mw,Eh;
354 float STAT,ee,EE,cc,ff,FF,Lw,Ew,Cw,Nw,Gn,
rho,norm,ch,CH,Cr,Mp,
N;
356 size_t i,
j,
k,
l,
m,Vm,lm,V,V4,V44,
id,
K,
M;
359 short* mm = this->skyMask.data;
360 short* MM = skyMMcc.
data;
361 bool skymaskcc = (skyMaskCC.size()==
L);
378 std::vector<float*> _vtd;
379 std::vector<float*> _vTD;
380 std::vector<float*> _eTD;
381 std::vector<float*> _AVX;
382 std::vector<float*> _APN;
383 std::vector<float*> _DAT;
384 std::vector<float*> _SIG;
385 std::vector<float*> _NUL;
386 std::vector<float*> _TMP;
388 for(i=0; i<
NIFO; i++) {
390 ml[
i] = getifo(i)->index.
data;
391 FP[
i] = getifo(i)->fp.data;
392 FX[
i] = getifo(i)->fx.data;
395 ml[
i] = getifo(0)->index.data;
396 FP[
i] = getifo(0)->fp.data;
397 FX[
i] = getifo(0)->fx.data;
410 std::vector<int>* vint;
411 std::vector<int>* vtof;
417 std::map<int,float> vLr;
421 int csize = precision%65536;
422 int healpix = this->nSkyStat.getOrder();
423 int order = (precision-csize)/65536;
426 if(healpix && csize && order && order<healpix) {
428 for(
int l=0;l<rsm.
size();l++) {
429 int m = this->nSkyStat.getSkyIndex(rsm.
getTheta(l),rsm.
getPhi(l));
432 for(
int l=0;l<
L;l++) BB[l] = BB[l] ? 0 : 1;
439 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
440 cTo = pwc->
get((
char*)
"time",0,
'L',0);
444 id = size_t(cid.
data[k]+0.1);
446 if(pwc->
sCuts[
id-1] != -2)
continue;
448 vint = &(pwc->
cList[
id-1]);
449 vtof = &(pwc->
nTofF[
id-1]);
453 pI = wdmMRA.getXTalk(pwc,
id);
457 bBB = (V>this->wdmMRA.nRes*csize) ?
true :
false;
460 this->nSensitivity = 0.;
461 this->nAlignment = 0.;
462 this->nNetIndex = 0.;
463 this->nDisbalance = 0.;
464 this->nLikelihood = 0.;
465 this->nNullEnergy = 0.;
466 this->nCorrEnergy = 0.;
467 this->nCorrelation = 0.;
469 this->nEllipticity = 0.;
470 this->nPolarisation = 0.;
471 this->nProbability = 0.;
473 this->nAntenaPrior = 0.;
476 tsize = pix->
tdAmp[0].size();
477 if(!tsize || tsize&1) {
478 cout<<
"network::likelihoodWP() error: wrong pixel TD data\n";
484 if(!(V=pI.size()))
continue;
485 V4 = V + (V%4 ? 4 - V%4 : 0);
488 for(j=0; j<V4; j++) pJ.push_back(0);
498 for(i=0; i<
NIFO; i++) {
499 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
500 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vtd.push_back(ptmp);
501 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
502 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vTD.push_back(ptmp);
503 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
504 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _eTD.push_back(ptmp);
505 ptmp = (
float*)_mm_malloc((V4*3+16)*
sizeof(float),32);
506 for(j=0; j<(V4*3+16); j++) ptmp[j]=0; _APN.push_back(ptmp);
507 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
508 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _DAT.push_back(ptmp);
509 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
510 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _SIG.push_back(ptmp);
511 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
512 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _NUL.push_back(ptmp);
517 this->p00_POL[
i].resize(V4); this->p00_POL[
i]=0.;
518 this->p90_POL[
i].resize(V4); this->p90_POL[
i]=0.;
519 this->r00_POL[
i].resize(V4); this->r00_POL[
i]=0.;
520 this->r90_POL[
i].resize(V4); this->r90_POL[
i]=0.;
523 for(i=0; i<
NIFO; i++) {
524 pa[
i] = _vtd[
i] + (tsize/2)*V4;
525 pA[
i] = _vTD[
i] + (tsize/2)*V4;
526 pe[
i] = _eTD[
i] + (tsize/2)*V4;
527 pd[
i] = _DAT[
i]; pD[
i] = _DAT[
i]+V4;
528 ps[
i] = _SIG[
i]; pS[
i] = _SIG[
i]+V4;
529 pn[
i] = _NUL[
i]; pN[
i] = _NUL[
i]+V4;
532 this->a_00.resize(NIFO*V44); this->a_00=0.;
533 this->a_90.resize(NIFO*V44); this->a_90=0.;
534 this->rNRG.resize(V4); this->rNRG=0.;
535 this->pNRG.resize(V4); this->pNRG=0.;
537 __m128* _aa = (__m128*) this->a_00.data;
538 __m128* _AA = (__m128*) this->a_90.data;
541 float* p_et = (
float*)_mm_malloc(V4*
sizeof(
float),32);
542 for(j=0; j<V4; j++) p_et[j]=0; _AVX.push_back(p_et);
543 float* pMSK = (
float*)_mm_malloc(V44*
sizeof(
float),32);
544 for(j=0; j<V44; j++) pMSK[j]=0; _AVX.push_back(pMSK); pMSK[V4]=
nIFO;
545 float* p_fp = (
float*)_mm_malloc(V44*
sizeof(
float),32);
546 for(j=0; j<V44; j++) p_fp[j]=0; _AVX.push_back(p_fp);
547 float* p_fx = (
float*)_mm_malloc(V44*
sizeof(
float),32);
548 for(j=0; j<V44; j++) p_fx[j]=0; _AVX.push_back(p_fx);
549 float* p_si = (
float*)_mm_malloc(V4*
sizeof(
float),32);
550 for(j=0; j<V4; j++) p_si[j]=0; _AVX.push_back(p_si);
551 float* p_co = (
float*)_mm_malloc(V4*
sizeof(
float),32);
552 for(j=0; j<V4; j++) p_co[j]=0; _AVX.push_back(p_co);
553 float* p_uu = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
554 for(j=0; j<V4+16; j++) p_uu[j]=0; _AVX.push_back(p_uu);
555 float* p_UU = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
556 for(j=0; j<V4+16; j++) p_UU[j]=0; _AVX.push_back(p_UU);
557 float* p_vv = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
558 for(j=0; j<V4+16; j++) p_vv[j]=0; _AVX.push_back(p_vv);
559 float* p_VV = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
560 for(j=0; j<V4+16; j++) p_VV[j]=0; _AVX.push_back(p_VV);
561 float* p_au = (
float*)_mm_malloc(V4*
sizeof(
float),32);
562 for(j=0; j<V4; j++) p_au[j]=0; _AVX.push_back(p_au);
563 float* p_AU = (
float*)_mm_malloc(V4*
sizeof(
float),32);
564 for(j=0; j<V4; j++) p_AU[j]=0; _AVX.push_back(p_AU);
565 float* p_av = (
float*)_mm_malloc(V4*
sizeof(
float),32);
566 for(j=0; j<V4; j++) p_av[j]=0; _AVX.push_back(p_av);
567 float* p_AV = (
float*)_mm_malloc(V4*
sizeof(
float),32);
568 for(j=0; j<V4; j++) p_AV[j]=0; _AVX.push_back(p_AV);
569 float* p_uv = (
float*)_mm_malloc(V4*4*
sizeof(
float),32);
570 for(j=0; j<V4*4; j++) p_uv[j]=0; _AVX.push_back(p_uv);
571 float* p_ee = (
float*)_mm_malloc(V4*
sizeof(
float),32);
572 for(j=0; j<V4; j++) p_ee[j]=0; _AVX.push_back(p_ee);
573 float* p_EE = (
float*)_mm_malloc(V4*
sizeof(
float),32);
574 for(j=0; j<V4; j++) p_EE[j]=0; _AVX.push_back(p_EE);
575 float* pTMP=(
float*)_mm_malloc(V4*4*NIFO*
sizeof(
float),32);
576 for(j=0; j<V4*4*
NIFO; j++) pTMP[j]=0; _AVX.push_back(pTMP);
577 float* p_ni = (
float*)_mm_malloc(V4*
sizeof(
float),32);
578 for(j=0; j<V4; j++) p_ni[j]=0; _AVX.push_back(p_ni);
579 float* p_ec = (
float*)_mm_malloc(V4*
sizeof(
float),32);
580 for(j=0; j<V4; j++) p_ec[j]=0; _AVX.push_back(p_ec);
581 float* p_gn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
582 for(j=0; j<V4; j++) p_gn[j]=0; _AVX.push_back(p_gn);
583 float* p_ed = (
float*)_mm_malloc(V4*
sizeof(
float),32);
584 for(j=0; j<V4; j++) p_ed[j]=0; _AVX.push_back(p_ed);
585 float* p_rn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
586 for(j=0; j<V4; j++) p_rn[j]=0; _AVX.push_back(p_rn);
592 this->pList.push_back(pix);
595 for(i=0; i<
nIFO; i++) {
596 xx[
i] = 1./pix->
data[
i].noiserms;
601 for(i=0; i<
nIFO; i++) {
602 _APN[
i][V4*2+
j] =(float)xx[i]/rms;
603 for(l=0; l<tsize; l++) {
605 AA = pix->
tdAmp[
i].data[l+tsize];
606 _vtd[
i][l*V4+
j] = aa;
607 _vTD[
i][l*V4+
j] = AA;
608 _eTD[
i][l*V4+
j] = aa*aa+AA*AA;
621 STAT=-1.e12; lm=0; Em=Lm=ff=FF=0;
624 for(l=lb; l<=le; l++) {
626 if(bBB && !BB[l])
continue;
633 if (!skyMaskCC.data[lc])
continue;
637 if(aa > gama) ff += 1;
639 REG[1] = (FF*FF/(ff*ff+1.e-9)-1)*En;
644 for(l=lb; l<=le; l++) {
645 skyProb.data[
l] = -1.e12;
647 pnt_(v00, pa, ml, (
int)l, (
int)V4);
648 pnt_(v90, pA, ml, (
int)l, (
int)V4);
655 if(lb==le) _avx_saveGW_ps(ps,pS,V);
659 _mm256_storeu_ps(vvv,_CC);
664 CH = No/(nIFO*Mo+sqrt(Mo));
666 Co = Ec/(Ec+No*cc-Mo*(nIFO-1));
668 if(Cr<
netCC)
continue;
670 aa = Eo>0. ? Eo-No : 0.;
672 skyProb.data[
l] = this->
delta<0 ? aa : AA;
676 if(pMSK[j]<=0)
continue;
678 ff += p_fp[
j]*p_et[
j];
679 FF += p_fx[
j]*p_et[
j];
681 ff = ee>0. ? ff/ee : 0.;
682 FF = ee>0. ? FF/ee : 0.;
683 this->nAntenaPrior.set(l, sqrt(ff+FF));
686 this->nSensitivity.set(l, sqrt(ff+FF));
687 this->nAlignment.set(l, ff>0 ? sqrt(FF/ff):0);
688 this->nLikelihood.set(l, Eo-No);
689 this->nNullEnergy.set(l, No);
690 this->nCorrEnergy.set(l, Ec);
691 this->nCorrelation.set(l,Co);
692 this->nSkyStat.set(l,AA);
693 this->nProbability.set(l, skyProb.data[l]);
694 this->nDisbalance.set(l,CH);
695 this->nNetIndex.set(l,cc);
696 this->nEllipticity.set(l,Cr);
697 this->nPolarisation.set(l,Mp);
700 if(AA>=STAT) {STAT=AA; lm=
l; Em=Eo-Eh;}
701 if(skyProb.data[l]>sky) sky=skyProb.data[
l];
707 D_snr = _avx_norm_ps(pd,pD,_AVX,V4);
708 S_snr = _avx_norm_ps(pS,pD,p_ec,V4);
713 _mm256_storeu_ps(vvv,_CC);
724 D_snr = _avx_norm_ps(pd,pD,_AVX,-V4);
725 N_snr = _avx_norm_ps(pn,pN,_AVX,-V4);
729 norm = Em>0 ? (Eo-Eh)/Em : 1.e9;
733 ch = (Np+Gn)/(N*nIFO);
735 rho = Ec>0 ? sqrt(Ec*Rc/2.) : 0.;
744 if(le-lb) {lb=le=lm;
goto optsky;}
746 if(Lm<=0.||(Eo-Eh)<=0.||Ec*Rc/cc<netEC||N<1) {
747 pwc->
sCuts[
id-1]=1; count=0;
748 pwc->
clean(
id);
continue;
755 vint = &(pwc->
cList[
id-1]);
756 for(j=0; j<vint->size(); j++) {
774 for(i=0; i<
nIFO; i++) {
775 pix->
setdata(
double(pd[i][j]),
'W',i);
776 pix->
setdata(
double(pD[i][j]),
'U',i);
777 pix->
setdata(
double(ps[i][j]),
'S',i);
778 pix->
setdata(
double(pS[i][j]),
'P',i);
784 if(!pix->
core)
continue;
785 if(p_gn[j]<=0)
continue;
790 if(!xpix->
core || p_gn[k]<=0 || xt.
CC[0]>2)
continue;
791 for(i=0; i<
nIFO; i++) {
799 if(p_ec[j]<=0)
continue;
805 if(!xpix->
core || p_ec[k]<=0 || xt.
CC[0]>2)
continue;
806 for(i=0; i<
nIFO; i++) {
818 for(j=1; j<=
nIFO; j++) {
819 if(S_snr[j]>Emax) Emax=S_snr[
j];
821 double Esub = S_snr.
data[0]-Emax;
822 Esub = Esub*(1+2*Rc*Esub/Emax);
823 Nmax = Gn+Np-N*(nIFO-1);
834 NETX (vtof->push_back(ml[0][lm]); ,
835 vtof->push_back(ml[1][lm]); ,
836 vtof->push_back(ml[2][lm]); ,
837 vtof->push_back(ml[3][lm]); ,
838 vtof->push_back(ml[4][lm]); ,
839 vtof->push_back(ml[5][lm]); ,
840 vtof->push_back(ml[6][lm]); ,
841 vtof->push_back(ml[7][lm]); )
844 if((wfsave)||(mdcListSize() && !lag)) {
845 if(this->getMRAwave(
id,lag,
'S',0,
true)) {
847 for(i=0; i<
nIFO; i++) {
848 pd = this->getifo(i);
849 pd->
RWFID.push_back(
id);
853 pd->
RWFP.push_back(wf);
856 if(this->getMRAwave(
id,lag,
's',0,
true)) {
858 for(i=0; i<
nIFO; i++) {
859 pd = this->getifo(i);
860 pd->
RWFID.push_back(-
id);
864 pd->
RWFP.push_back(wf);
869 Lw = Ew = To = Fo = Nw = ee = norm = 0.;
870 for(i=0; i<
nIFO; i++) {
875 this->getMRAwave(
id,lag,
'W',0);
876 this->getMRAwave(
id,lag,
'S',0);
877 for(i=0; i<
nIFO; i++) {
896 ch = (Nw+Gn)/(N*nIFO);
897 cc = ch>1 ? 1+(ch-1)*2*(1-Rc) : 1;
898 Cr = Ec*Rc/(Ec*Rc+(Dc+Nw+Gn)*cc-N*(nIFO-1));
900 Cp = Ec*Rc/(Ec*Rc+(Dc+Nw+Gn)-N*(nIFO-1));
903 pwc->
cData[
id-1].norm = norm*2;
904 pwc->
cData[
id-1].skyStat = 0;
905 pwc->
cData[
id-1].skySize = Mw;
906 pwc->
cData[
id-1].netcc = Cp;
907 pwc->
cData[
id-1].skycc = Cr;
908 pwc->
cData[
id-1].subnet = Esub/(Esub+Nmax);
909 pwc->
cData[
id-1].SUBNET = Co;
910 pwc->
cData[
id-1].likenet = Lw;
911 pwc->
cData[
id-1].netED = Nw+Gn+Dc-N*
nIFO;
912 pwc->
cData[
id-1].netnull = Nw+Gn;
913 pwc->
cData[
id-1].energy = Ew;
914 pwc->
cData[
id-1].likesky = Em;
915 pwc->
cData[
id-1].enrgsky = Eo;
916 pwc->
cData[
id-1].netecor = Ec;
917 pwc->
cData[
id-1].normcor = Ec*Rc;
918 pwc->
cData[
id-1].netRHO = rho/sqrt(cc);
920 pwc->
cData[
id-1].cTime = To;
921 pwc->
cData[
id-1].cFreq = Fo;
922 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
923 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
924 pwc->
cData[
id-1].gNET = sqrt(ff+FF);
925 pwc->
cData[
id-1].aNET = sqrt(FF/ff);
926 pwc->
cData[
id-1].iNET = 0;
927 pwc->
cData[
id-1].nDoF =
N;
928 pwc->
cData[
id-1].skyChi2 = CH;
929 pwc->
cData[
id-1].Gnoise = Gn;
930 pwc->
cData[
id-1].iota = 0;
931 pwc->
cData[
id-1].psi = 0;
932 pwc->
cData[
id-1].ellipticity = 0.;
934 cc = pwc->
cData[
id-1].netcc;
936 printf(
"rho=%4.2f|%4.2f cc: %5.3f|%5.3f|%5.3f subnet=%4.3f|%4.3f \n",
937 rho,rho*sqrt(Cp),Co,Cp,Cr,pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
938 printf(
" L: %5.1f|%5.1f|%5.1f E: %5.1f|%5.1f|%5.1f|%5.1f N: %4.1f|%4.1f|%4.1f|%4.1f|%4.1f \n",
939 Lw,Lp,Lo,Ew,Ep,Eo,Em,Nw,Np,Rc,Eh,No);
940 printf(
"id|lm %3d|%6d Vm|m=%3d|%3d|%3d|%3d T|F: %6.3f|%4.1f (t,p)=(%4.1f|%4.1f) \n",
941 int(
id),
int(lm),
int(V),
int(Mo),
int(Mw),
int(M),To,Fo,nLikelihood.getTheta(lm),nLikelihood.getPhi(lm));
942 cout<<
" L: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",S_snr[i]);}
943 cout<<
" E: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",D_snr[i]);}
944 cout<<
" N: |";
for(i=1; i<nIFO+1; i++) {
printf(
"%5.1f|",N_snr[i]);}
945 cout<<endl<<
" dof|G|G+R ";
printf(
"%5.1f|%5.1f|%5.1f r[1]=%4.1f",N,Gn,Nw+Gn,REG[1]);
946 printf(
" norm=%3.1f chi2 %3.2f|%3.2f Rc=%3.2f, Dc=%4.1f\n",norm,ch,CH,Rc,Dc);
954 pwc->
p_Ind[
id-1].push_back(Mo);
956 std::vector<float> sArea;
957 pwc->
sArea.push_back(sArea);
958 pwc->
p_Map.push_back(sArea);
960 double var = norm*Rc*sqrt(Mo)*(1+
fabs(1-CH));
962 if(iID<=0 || ID==
id) {
963 getSkyArea(
id,lag,T,var);
967 pwc->
cData[
id-1].mchirp = 0;
968 pwc->
cData[
id-1].mchirperr = 0;
969 pwc->
cData[
id-1].tmrgr = 0;
970 pwc->
cData[
id-1].tmrgrerr = 0;
971 pwc->
cData[
id-1].chi2chirp = 0;
975 cc = Ec/(
fabs(Ec)+ee);
976 printf(
"mchirp : %d %g %.2e %.3f %.3f %.3f %.3f \n\n",
977 int(
id),cc,pwc->
cData[
id-1].mchirp,
978 pwc->
cData[
id-1].mchirperr, pwc->
cData[
id-1].tmrgr,
979 pwc->
cData[
id-1].tmrgrerr, pwc->
cData[
id-1].chi2chirp);
982 if(ID==
id && !
EFEC) {
983 this->nSensitivity.gps =
T;
984 this->nAlignment.gps =
T;
985 this->nDisbalance.gps =
T;
986 this->nLikelihood.gps =
T;
987 this->nNullEnergy.gps =
T;
988 this->nCorrEnergy.gps =
T;
989 this->nCorrelation.gps =
T;
990 this->nSkyStat.gps =
T;
991 this->nEllipticity.gps =
T;
992 this->nPolarisation.gps=
T;
993 this->nNetIndex.gps =
T;
996 pwc->
sCuts[
id-1] = -1;
1023 if(!this->wc_List[lag].
size())
return 0;
1025 size_t nIFO = this->ifoList.size();
1028 cout<<
"network::subNetCut(): invalid network.\n";
1035 subnet =
fabs(subnet);
1036 subcut =
fabs(subcut);
1038 __m128 _En = _mm_set1_ps(En);
1039 __m128 _Es = _mm_set1_ps(Es);
1040 __m128 _oo = _mm_set1_ps(1.
e-12);
1041 __m128 _0 = _mm_set1_ps(0.);
1042 __m128 _05 = _mm_set1_ps(0.5);
1043 __m128 _1 = _mm_set1_ps(1.);
1048 float Lm,Em,Am,Lo,Eo,Co,Lr,Er,ee,em,To;
1049 float cc,aa,AA,rHo,stat,Ls,Ln,EE;
1053 short* mm = this->skyMask.data;
1066 for(i=0; i<
NIFO; i++) {
1068 ml[
i] = getifo(i)->index.data;
1069 FP[
i] = getifo(i)->fp.data;
1070 FX[
i] = getifo(i)->fx.data;
1073 ml[
i] = getifo(0)->index.data;
1074 FP[
i] = getifo(0)->fp.data;
1075 FX[
i] = getifo(0)->fx.data;
1080 std::vector<int> pI;
1084 std::vector<int>* vint;
1094 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
1095 cTo = pwc->
get((
char*)
"time",0,
'L',0);
1098 for(k=0; k<
K; k++) {
1099 id = size_t(cid.
data[k]+0.1);
1100 if(pwc->
sCuts[
id-1] != -2)
continue;
1101 vint = &(pwc->
cList[
id-1]);
1105 pI = wdmMRA.getXTalk(pwc,
id);
1111 tsize = pix->
tdAmp[0].size();
1112 if(!tsize || tsize&1) {
1113 cout<<
"network::subNetCut() error: wrong pixel TD data\n";
1117 V4 = V + (V%4 ? 4 - V%4 : 0);
1121 std::vector<wavearray<float> > vtd;
1122 std::vector<wavearray<float> > vTD;
1123 std::vector<wavearray<float> > eTD;
1142 __m128* _Fp = (__m128*) Fp.
data;
1143 __m128* _Fx = (__m128*) Fx.
data;
1144 __m128* _am = (__m128*) am.
data;
1145 __m128* _AM = (__m128*) AM.
data;
1146 __m128* _xi = (__m128*) xi.
data;
1147 __m128* _XI = (__m128*) XI.
data;
1148 __m128* _fp = (__m128*) fp.
data;
1149 __m128* _fx = (__m128*) fx.
data;
1150 __m128* _nr = (__m128*) nr.
data;
1151 __m128* _ww = (__m128*) ww.
data;
1152 __m128* _WW = (__m128*) WW.
data;
1153 __m128* _bb = (__m128*) bb.
data;
1154 __m128* _BB = (__m128*) BB.
data;
1156 for(i=0; i<NIFO; i++) {
1162 for(i=0; i<
NIFO; i++) {
1163 pa[
i] = vtd[
i].data + (tsize/2)*V4;
1164 pA[
i] = vTD[
i].data + (tsize/2)*V4;
1165 pe[
i] = eTD[
i].data + (tsize/2)*V4;
1168 this->a_00.resize(NIFO*V4); this->a_00=0.;
1169 this->a_90.resize(NIFO*V4); this->a_90=0.;
1170 this->rNRG.resize(V4); this->rNRG=0.;
1171 this->pNRG.resize(V4); this->pNRG=0.;
1173 __m128* _aa = (__m128*) this->a_00.data;
1174 __m128* _AA = (__m128*) this->a_90.data;
1176 this->pList.clear();
1177 for(j=0; j<V; j++) {
1179 pList.push_back(pix);
1182 for(i=0; i<
nIFO; i++) {
1187 for(i=0; i<
nIFO; i++) {
1188 nr.
data[j*NIFO+
i]=(float)
xx[i]/sqrt(rms);
1189 for(l=0; l<tsize; l++) {
1191 AA = pix->
tdAmp[
i].data[l+tsize];
1192 vtd[
i].data[l*V4+
j] = aa;
1193 vTD[
i].data[l*V4+
j] = AA;
1194 eTD[
i].data[l*V4+
j] = aa*aa+AA*AA;
1208 bool skymaskcc = (skyMaskCC.size()==Lsky);
1210 stat=Lm=Em=Am=EE=0.; lm=Vm= -1;
1214 for(l=lb; l<=le; l++) {
1215 if(!mm[l] || l<0)
continue;
1222 if (!skyMaskCC.data[lc])
continue;
1228 __m128 _E_o = _mm_setzero_ps();
1229 __m128 _E_n = _mm_setzero_ps();
1230 __m128 _E_s = _mm_setzero_ps();
1231 __m128 _M_m = _mm_setzero_ps();
1232 __m128* _rE = (__m128*) rNRG.data;
1233 __m128* _pE = (__m128*) pNRG.data;
1235 for(j=0; j<V4; j+=4) {
1237 _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_rE,_En));
1238 _M_m = _mm_add_ps(_M_m,_msk);
1239 *_pE = _mm_mul_ps(*_rE,_msk);
1240 _E_o = _mm_add_ps(_E_o,*_pE);
1242 _E_s = _mm_add_ps(_E_s,*_pE);
1243 _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_pE++,_Es));
1244 _E_n = _mm_add_ps(_E_n,_mm_mul_ps(*_rE++,_msk));
1247 _mm_storeu_ps(vvv,_E_n);
1248 Ln = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1249 _mm_storeu_ps(vvv,_E_o);
1250 Eo = vvv[0]+vvv[1]+vvv[2]+vvv[3]+0.01;
1251 _mm_storeu_ps(vvv,_E_s);
1252 Ls = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1253 _mm_storeu_ps(vvv,_M_m);
1254 m = 2*(vvv[0]+vvv[1]+vvv[2]+vvv[3])+0.01;
1257 if((aa-m)/(aa+
m)<subcut)
continue;
1259 pnt_(v00, pa, ml, (
int)l, (
int)V4);
1260 pnt_(v90, pA, ml, (
int)l, (
int)V4);
1261 float* pfp = fp.
data;
1262 float* pfx = fx.
data;
1263 float* p00 = this->a_00.data;
1264 float* p90 = this->a_90.data;
1267 for(j=0; j<V; j++) {
1269 cpp_(p00,v00); cpp_(p90,v90);
1270 cpf_(pfp,FP,l); cpf_(pfx,FX,l);
1275 if(rNRG.data[j]>En) m++;
1278 __m128* _pp = (__m128*) am.
data;
1279 __m128* _PP = (__m128*) AM.
data;
1283 _pp = (__m128*) xi.
data;
1284 _PP = (__m128*) XI.
data;
1288 for(j=0; j<V; j++) {
1303 Ls += ee-em; Eo += ee;
1304 if(ee-em>Es) Ln += ee;
1308 size_t m4 = m + (m%4 ? 4 - m%4 : 0);
1309 _E_n = _mm_setzero_ps();
1311 for(j=0; j<m4; j+=4) {
1315 _E_n = _mm_add_ps(_E_n,_E_s);
1317 _mm_storeu_ps(vvv,_E_n);
1319 Lo = vvv[0]+vvv[1]+vvv[2]+vvv[3];
1320 AA = aa/(
fabs(aa)+
fabs(Eo-Lo)+2*m*(Eo-Ln)/Eo);
1322 em =
fabs(Eo-Lo)+2*
m;
1326 if(AA>stat && !mra) {
1327 stat=AA; Lm=Lo; Em=Eo; Am=aa; lm=
l; Vm=
m; suball=ee; EE=em;
1331 if(!mra && lm>=0) {mra=
true; le=lb=lm;
goto skyloop;}
1333 pwc->
sCuts[
id-1] = -1;
1334 pwc->
cData[
id-1].likenet = Lm;
1335 pwc->
cData[
id-1].energy = Em;
1336 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
1337 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
1338 pwc->
cData[
id-1].skyIndex = lm;
1342 submra = Ls*Eo/(Eo-Ls);
1343 submra/=
fabs(submra)+
fabs(Eo-Lo)+2*(m+6);
1345 pwc->
p_Ind[
id-1].push_back(lm);
1346 for(j=0; j<vint->size(); j++) {
1348 pix->
theta = nLikelihood.getTheta(lm);
1349 pix->
phi = nLikelihood.getPhi(lm);
1357 rHo = sqrt(Lo*Lo/(Eo+2*m)/2);
1360 if(hist && rHo>this->
netRHO)
1361 for(j=0;j<vint->size();j++) hist->Fill(suball,submra);
1363 if(fmin(suball,submra)>subnet && rHo>this->
netRHO) {
1364 count += vint->size();
1366 printf(
"lag|id %3d|%3d rho=%5.2f To=%5.1f stat: %5.3f|%5.3f|%5.3f ",
1367 int(lag),
int(
id),rHo,To,suball,submra,stat);
1368 printf(
"E: %6.1f|%6.1f L: %6.1f|%6.1f|%6.1f pix: %4d|%4d|%3d|%2d \n",
1369 Em,Eo,Lm,Lo,Ls,
int(vint->size()),
int(V),Vm,
int(m));
1372 else pwc->
sCuts[
id-1]=1;
1377 for(j=0; j<V; j++) {
1412 if(!this->wc_List[lag].
size())
return 0;
1418 bool cirwave = mode==
'g' || mode==
'G' || mode==
'c' || mode==
'C';
1419 bool linwave = mode==
'l' || mode==
'L' || mode==
's' || mode==
'S';
1420 bool iotwave = mode==
'i' || mode==
'l' || mode==
'e' || mode==
'c' ||
1421 mode==
'I' || mode==
'L' || mode==
'E' || mode==
'C';
1422 bool psiwave = mode==
'l' || mode==
'e' || mode==
'p' ||
1423 mode==
'L' || mode==
'E' || mode==
'P';
1424 bool mureana = mode==
'i' || mode==
'e' || mode==
'c' ||
1425 mode==
'r' || mode==
'p' || mode==
'b' ||
1426 mode==
'l' || mode==
's' || mode==
'g';
1427 bool rndwave = mode==
'r' || mode==
'R';
1429 bool prior = this->gamma<0 ? true :
false;
1430 bool m_chirp = this->
optim ? false : mureana;
1432 if(!this->
optim) mureana =
true;
1434 size_t nIFO = this->ifoList.size();
1435 size_t ID = abs(iID);
1438 cout<<
"network::likelihood2G(): invalid network.\n";
1444 float gama =
fabs(this->gamma);
1446 if(gama<=0) gama = 1.e-24;
1447 if(gama>=1) gama = 0.999999;
1455 static const __m128 _D0 = _mm_set1_ps(deta<0.5?1-deta:0.5);
1456 static const __m128 _D9 = _mm_set1_ps(deta<0.5?1:1.5-deta);
1458 static const __m128 _oo = _mm_set1_ps(1.
e-16);
1459 static const __m128 _sm = _mm_set1_ps(-0.
f);
1460 static const __m128 _En = _mm_set1_ps(En);
1461 static const __m128 _rG = _mm_set1_ps(-1./
log(gama));
1462 static const __m128 _kG = _mm_set1_ps(gama);
1463 static const __m128 _PW = _mm_set1_ps(psiwave?0:1);
1464 static const __m128 _01 = _mm_set1_ps(0.1);
1465 static const __m128 _05 = _mm_set1_ps(0.5);
1466 static const __m128 _09 = _mm_set1_ps(0.9);
1467 static const __m128 _1 = _mm_set1_ps(1.0+1.
e-16);
1468 static const __m128 _2 = _mm_set1_ps(2.0);
1469 static const __m128 _4 = _mm_set1_ps(4.0);
1474 float NRG,Lm,Em,Lo,Eo,No,Nm,cc,Cm,Co,Do,To,Fo,Ln,Ns;
1475 float STAT,CHR,aa,AA,ee,em,EE,ff,FF,Lr,Cr,
ss,Ls,Nc,gg;
1476 double eLp, s2p, c2p;
1479 size_t i,
j,
k,
l,
m,Vm,lm,V,V4,V44,
id,
K,
M;
1481 short* mm = this->skyMask.data;
1482 bool skymaskcc = (skyMaskCC.size()==
L);
1496 for(i=0; i<
NIFO; i++) {
1498 ml[
i] = getifo(i)->index.data;
1499 FP[
i] = getifo(i)->fp.data;
1500 FX[
i] = getifo(i)->fx.data;
1503 ml[
i] = getifo(0)->index.data;
1504 FP[
i] = getifo(0)->fp.data;
1505 FX[
i] = getifo(0)->fx.data;
1510 std::vector<int> pI;
1511 std::vector<int> pJ;
1515 std::vector<int>* vint;
1516 std::vector<int>* vtof;
1517 std::vector<int> pRate;
1523 std::map<int,float> vLr;
1529 cid = pwc->
get((
char*)
"ID", 0,
'S',0);
1530 cTo = pwc->
get((
char*)
"time",0,
'L',0);
1533 for(k=0; k<
K; k++) {
1534 id = size_t(cid.
data[k]+0.1);
1536 if(pwc->
sCuts[
id-1] != -2)
continue;
1538 vint = &(pwc->
cList[
id-1]);
1539 vtof = &(pwc->
nTofF[
id-1]);
1543 pI = wdmMRA.getXTalk(pwc,
id);
1548 this->nSensitivity = 0.;
1549 this->nAlignment = 0.;
1550 this->nNetIndex = 0.;
1551 this->nDisbalance = 0.;
1552 this->nLikelihood = 0.;
1553 this->nNullEnergy = 0.;
1554 this->nCorrEnergy = 0.;
1555 this->nCorrelation = 0.;
1556 this->nSkyStat = 0.;
1557 this->nEllipticity = 0.;
1558 this->nPolarisation = 0.;
1559 this->nProbability = 0.;
1561 this->nAntenaPrior = 0.;
1564 tsize = pix->
tdAmp[0].size();
1565 if(!tsize || tsize&1) {
1566 cout<<
"network::likelihood2G() error: wrong pixel TD data\n";
1572 if(!(V=pI.size()))
continue;
1573 V4 = V + (V%4 ? 4 - V%4 : 0);
1576 for(j=0; j<V4; j++) pJ.push_back(0);
1580 std::vector<wavearray<float> > vtd;
1581 std::vector<wavearray<float> > vTD;
1582 std::vector<wavearray<float> > eTD;
1605 this->p00_POL[
i].
resize(V4); this->p00_POL[
i]=0.;
1606 this->p90_POL[
i].resize(V4); this->p90_POL[
i]=0.;
1607 this->r00_POL[
i].resize(V4); this->r00_POL[
i]=0.;
1608 this->r90_POL[
i].resize(V4); this->r90_POL[
i]=0.;
1611 __m128* _Fp = (__m128*) Fp.
data;
1612 __m128* _Fx = (__m128*) Fx.
data;
1613 __m128* _am = (__m128*) am.
data;
1614 __m128* _AM = (__m128*) AM.
data;
1615 __m128* _xi = (__m128*) xi.
data;
1616 __m128* _XI = (__m128*) XI.
data;
1617 __m128* _xp = (__m128*) xp.
data;
1618 __m128* _XP = (__m128*) XP.
data;
1619 __m128* _ww = (__m128*) ww.
data;
1620 __m128* _WW = (__m128*) WW.
data;
1621 __m128* _bb = (__m128*) bb.
data;
1622 __m128* _BB = (__m128*) BB.
data;
1623 __m128* _fp = (__m128*) fp.
data;
1624 __m128* _fx = (__m128*) fx.
data;
1625 __m128* _nr = (__m128*) nr.
data;
1626 __m128* _ep = (__m128*) ep.
data;
1627 __m128* _ex = (__m128*) ex.
data;
1629 __m128* _fp4 = _fp+V4*f_;
1630 __m128* _fx4 = _fx+V4*f_;
1631 __m128* _uu4 = _am+V4*f_;
1632 __m128* _vv4 = _AM+V4*f_;
1633 __m128* _bb4 = _bb+V4*f_;
1634 __m128* _BB4 = _BB+V4*f_;
1636 for(i=0; i<NIFO; i++) {
1642 for(i=0; i<
NIFO; i++) {
1643 pa[
i] = vtd[
i].data + (tsize/2)*V4;
1644 pA[
i] = vTD[
i].data + (tsize/2)*V4;
1645 pe[
i] = eTD[
i].data + (tsize/2)*V4;
1663 this->a_00.
resize(NIFO*V44); this->a_00=0.;
1664 this->a_90.resize(NIFO*V44); this->a_90=0.;
1665 this->rNRG.resize(V4); this->rNRG=0.;
1666 this->pNRG.resize(V4); this->pNRG=0.;
1668 __m128* _aa = (__m128*) this->a_00.data;
1669 __m128* _AA = (__m128*) this->a_90.data;
1672 this->pList.clear(); pRate.clear();
1673 for(j=0; j<V; j++) {
1675 this->pList.push_back(pix);
1676 pRate.push_back(
int(pix->
rate+0.5));
1678 if(vLr.find(pRate[j]) == vLr.end())
1682 for(i=0; i<
nIFO; i++) {
1688 for(i=0; i<
nIFO; i++) {
1689 nr.
data[j*NIFO+
i]=(float)
xx[i]/rms;
1690 for(l=0; l<tsize; l++) {
1692 AA = pix->
tdAmp[
i].data[l+tsize];
1693 vtd[
i].data[l*V4+
j] = aa;
1694 vTD[
i].data[l*V4+
j] = AA;
1695 eTD[
i].data[l*V4+
j] = aa*aa+AA*AA;
1704 STAT=0.; lm=0; Vm=0;
1705 double skystat = 0.;
1713 for(l=lb; l<=le; l++) {
1714 if(!mra) skyProb.data[
l] = 0.;
1715 if(!mm[l])
continue;
1722 if (!skyMaskCC.data[lc])
continue;
1725 pnt_(v00, pa, ml, (
int)l, (
int)V4);
1726 pnt_(v90, pA, ml, (
int)l, (
int)V4);
1727 float* pfp = fp.
data;
1728 float* pfx = fx.
data;
1729 float* p00 = this->a_00.data;
1730 float* p90 = this->a_90.data;
1734 for(j=0; j<V; j++) {
1735 cpp_(p00,v00); cpp_(p90,v90);
1736 cpf_(pfp,FP,l); cpf_(pfx,FX,l);
1737 if(!this->
optim || !mra)
continue;
1738 if(vLr[pRate[j]] <= mxLr)
continue;
1739 mxLr = vLr[pRate[
j]];
1744 for(j=0; j<V; j++) {
1749 if(optR && optR!=pRate[j]) {
1757 if(ee>En) m++;
else ee=0.;
1759 pNRG.data[
j] = rNRG.data[
j];
1762 __m128* _pp = (__m128*) am.
data;
1763 __m128* _PP = (__m128*) AM.
data;
1765 if(mra && mureana) {
1767 _pp = (__m128*) xi.
data;
1768 _PP = (__m128*) XI.
data;
1772 for(j=0; j<V; j++) {
1789 size_t m4 = m + (m%4 ? 4 - m%4 : 0);
1790 __m128 _ll,_LL,_ec,_EC,_ee,_EE,_NI,_s2,_c2,_AX,_NN,_FF,_QQ,_ie,_gg;
1791 __m128 _en,_EN,_ed,_ED,_cc,_ss,_ni,_si,_co,_ax,_nn,_ff,_mm,_IE,_GG;
1793 __m128* _siP = (__m128*) siPHS.
data;
1794 __m128* _coP = (__m128*) coPHS.
data;
1795 __m128* _siO = (__m128*) siORT.
data;
1796 __m128* _coO = (__m128*) coORT.
data;
1797 __m128* _siD = (__m128*) siDPF.
data;
1798 __m128* _coD = (__m128*) coDPF.
data;
1799 __m128* _nrg = (__m128*) q2Q2.
data;
1800 __m128* _chr = (__m128*) chir.
data;
1803 _pp = (__m128*) xi.
data;
1804 _PP = (__m128*) XI.
data;
1809 Lo = Ln = Co = Eo = 0.;
1810 for(j=0; j<m4; j+=4) {
1813 __m128* _pbb = _bb+jf;
1814 __m128* _pBB = _BB+jf;
1815 __m128* _pxi = _pp+jf;
1816 __m128* _pXI = _PP+jf;
1817 __m128* _pxp = _xp+jf;
1818 __m128* _pXP = _XP+jf;
1819 __m128* _pww = _ww+jf;
1820 __m128* _pWW = _WW+jf;
1821 __m128* _pfp = _fp+jf;
1822 __m128* _pfx = _fx+jf;
1823 __m128* _pFp = _Fp+jf;
1824 __m128* _pFx = _Fx+jf;
1837 _coO++; _siO++; _coD++; _siD++;
1840 if(le==lb && (optR==0)) {
1841 _sse_pol4_ps(_pfp, _pfx, _pxp, p00_POL[0].data+j, p00_POL[1].data+j);
1842 _sse_pol4_ps(_pfp, _pfx, _pXP, p90_POL[0].data+j, p90_POL[1].data+j);
1843 _sse_pol4_ps(_pfp, _pfx, _pxi, r00_POL[0].data+j, r00_POL[1].data+j);
1844 _sse_pol4_ps(_pfp, _pfx, _pXI, r90_POL[0].data+j, r90_POL[1].data+j);
1855 _cc = _mm_and_ps(_mm_cmpge_ps(_ec,_05),_1);
1856 _ss = _mm_and_ps(_mm_cmpgt_ps(_EC,_05),_cc);
1861 _nn = _mm_add_ps(_mm_add_ps(_ee,_EE),_2);
1862 _cc = _mm_add_ps(_ec,_mm_sub_ps(_nn,_ll));
1863 _cc = _mm_div_ps(_ec,_mm_add_ps(_cc,_mm));
1865 _mm_storeu_ps(vvv,_mm_add_ps(_ee,_EE));
1866 Eo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1867 _mm_storeu_ps(vvv,_ll);
1868 Lo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1869 _mm_storeu_ps(vvv,_mm_mul_ps(_ll,_cc));
1870 Ln += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1871 _mm_storeu_ps(vvv,_ec);
1872 Co += vvv[0]+vvv[1]+vvv[2]+vvv[3];
1873 _mm_storeu_ps(vvv,_mm_mul_ps(_ec,_ll));
1874 if(le==lb && !mra) {
1875 vLr[pRate[pJ[j+0]]] += vvv[0];
1876 vLr[pRate[pJ[j+1]]] += vvv[1];
1877 vLr[pRate[pJ[j+2]]] += vvv[2];
1878 vLr[pRate[pJ[j+3]]] += vvv[3];
1881 Ln = Eo>0 ? Ln/Eo : 0;
1882 if(Ln<this->
netCC)
continue;
1883 _IE = _mm_set1_ps(1-Co/Lo);
1888 __m128* _xx = (__m128*) xxx.
data;
1889 __m128* _yy = (__m128*) yyy.
data;
1890 __m128* _zz = (__m128*) zzz.
data;
1891 __m128* _rr = (__m128*) rrr.
data;
1893 for(j=0; j<m4; j+=4) {
1895 __m128* _pfp = _fp+jf;
1896 __m128* _pfx = _fx+jf;
1899 _ee = _mm_add_ps(_mm_sqrt_ps(*_xx++),_oo);
1900 _EE = _mm_add_ps(_mm_sqrt_ps(*_yy++),_oo);
1905 _xx = (__m128*) xxx.
data;
1906 _yy = (__m128*) yyy.
data;
1907 _zz = (__m128*) zzz.
data;
1908 _rr = (__m128*) rrr.
data;
1909 _siP = (__m128*) siPHS.
data;
1910 _coP = (__m128*) coPHS.
data;
1911 _siD = (__m128*) siDPF.
data;
1912 _coD = (__m128*) coDPF.
data;
1913 _nrg = (__m128*) q2Q2.
data;
1915 __m128 linw = linwave ? _oo : _1;
1916 __m128 cirw = cirwave ? _oo : _1;
1917 __m128 _ch = _mm_setzero_ps();
1918 __m128 _eqQ = _mm_setzero_ps();
1919 __m128 _qQ2 = _mm_setzero_ps();
1921 _c2 = _mm_setzero_ps();
1922 _s2 = _mm_setzero_ps();
1923 _cc = _mm_setzero_ps();
1924 _ss = _mm_setzero_ps();
1926 for(j=0; j<m4; j+=4) {
1928 __m128* _pfp = _fp+jf;
1929 __m128* _pfx = _fx+jf;
1930 __m128* _pep = _ep+jf;
1931 __m128* _pex = _ex+jf;
1932 __m128* _pxp = _xp+jf;
1933 __m128* _pXP = _XP+jf;
1934 __m128* _pxi = _pp+jf;
1935 __m128* _pXI = _PP+jf;
1937 _ee = _mm_div_ps(_1,_mm_add_ps(*_xx,_oo));
1943 _nn = _mm_div_ps(_1,_mm_add_ps(_NN,_oo));
1944 *_siP = _mm_mul_ps(_si,_nn);
1945 *_coP = _mm_mul_ps(_co,_nn);
1947 if(le==lb && (optR==0)) {
1948 __m128 _snn = _mm_sqrt_ps(_nn);
1951 _sse_pol4_ps(_ep+jf, _ex+jf, _uu4, r00_POL[0].data+j, r00_POL[1].data+j);
1952 _sse_pol4_ps(_ep+jf, _ex+jf, _vv4, r90_POL[0].data+j, r90_POL[1].data+j);
1957 _ch = _mm_add_ps(_ch,_EN);
1958 _ll = _mm_and_ps(_mm_cmpgt_ps(_EN,_oo),_1);
1959 *_chr = _mm_sub_ps(_mm_mul_ps(_ll,_2),_1);
1965 _ff = _mm_div_ps(_mm_add_ps(*_xx,*_yy),_ni);
1968 _LL = _mm_mul_ps(_LL,_nn);
1969 _FF = _mm_mul_ps(*_yy,_ee);
1970 _gg = _mm_mul_ps(_01,_ff);
1971 _NN = _mm_mul_ps(_NN,*_xx);
1972 _ll = _mm_mul_ps(_NN,_mm_add_ps(_FF,_gg));
1973 _co = _mm_sub_ps(_ll,_LL);
1976 *_zz = _mm_add_ps(_mm_sqrt_ps(*_zz),_oo);
1977 _ll = _mm_mul_ps(_NN,_FF);
1978 _QQ = _mm_add_ps(_mm_add_ps(_ll,_LL),*_zz);
1979 _ff = _mm_sqrt_ps(_mm_mul_ps(_ff,_05));
1980 _FF = _mm_sqrt_ps(_FF);
1982 _ax = _mm_add_ps(_1,_mm_div_ps(_co,*_zz));
1983 _ax = _mm_sqrt_ps(_mm_mul_ps(_ax,_05));
1985 _qQ2 = _mm_add_ps(_qQ2,_QQ);
1986 _eqQ = _mm_add_ps(_eqQ,_AX);
1987 _AX = _mm_div_ps(_mm_mul_ps(_2,_AX),_QQ);
1989 _gg = _mm_mul_ps(_mm_sub_ps(_05,_ni),
1990 _mm_sub_ps(_1,_FF));
1991 _GG = _mm_sub_ps(_mm_sub_ps(_09,_GG),_gg);
1992 _GG = _mm_mul_ps(_mm_mul_ps(_IE,_rG),_GG);
1993 _GG = _mm_mul_ps(_mm_mul_ps(_IE,_4),_GG);
1995 _gg = _mm_mul_ps(_mm_sub_ps(_IE,_ff),_kG);
1996 _gg = _mm_mul_ps(_ni,_gg);
1997 _nn = _mm_and_ps(_mm_cmpgt_ps(_ff,_gg),_1);
2001 _mm = _mm_and_ps(_mm_cmpgt_ps(_ff,_GG),_nn);
2003 _gg = _mm_andnot_ps(_sm,_mm_sub_ps(_NI,_ni));
2004 _gg = _mm_sub_ps(_IE,_mm_mul_ps(_gg,_FF));
2005 _nn = _mm_and_ps(_mm_cmplt_ps(_gg,_D0),_mm);
2006 _mm = _mm_and_ps(_mm_cmplt_ps(_gg,_D9),_mm);
2007 _nn = _mm_mul_ps(_nn,cirw);
2008 _mm = _mm_mul_ps(_mm,linw);
2009 _EE = _mm_mul_ps(_en,_nn);
2013 _nn = _mm_mul_ps(_nn,_PW);
2014 *_rr = _mm_add_ps(_mm_add_ps(_mm,_nn),_1);
2015 _coP++;_siP++;_chr++;_xx++;_yy++;_zz++,_rr++;
2017 _cc = _mm_add_ps(_cc,_ee);
2018 _ss = _mm_add_ps(_ss,_EE);
2024 _c2 = _mm_add_ps(_c2,_ec);
2025 _s2 = _mm_sub_ps(_s2,_EC);
2029 _mm_storeu_ps(vvv,_c2);
2030 c2p = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2031 _mm_storeu_ps(vvv,_s2);
2032 s2p = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2033 gg = sqrt(c2p*c2p+s2p*s2p+1.
e-16);
2034 _si = _mm_set1_ps(s2p/gg);
2035 _co = _mm_set1_ps(c2p/gg);
2038 _zz = (__m128*) zzz.
data;
2039 _siD = (__m128*) siDPF.
data;
2040 _coD = (__m128*) coDPF.
data;
2042 _mm_storeu_ps(vvv,_cc);
2043 cc = (vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2044 _mm_storeu_ps(vvv,_ss);
2045 ss = (vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2046 gg = sqrt(cc*cc+ss*ss+1.
e-16);
2047 _si = _mm_set1_ps(ss/gg);
2048 _co = _mm_set1_ps(cc/gg);
2050 for(j=0; j<m4; j+=4) {
2052 __m128* _pxi = _pp+jf;
2053 __m128* _pXI = _PP+jf;
2071 _siO = (__m128*) siORT.
data;
2072 _coO = (__m128*) coORT.
data;
2073 _siP = (__m128*) siPHS.
data;
2074 _coP = (__m128*) coPHS.
data;
2075 _chr = (__m128*) chir.
data;
2077 _mm_storeu_ps(vvv,_eqQ);
2078 _mm_storeu_ps(uuu,_qQ2);
2079 eLp = 2.*(vvv[0]+vvv[1]+vvv[2]+vvv[3]);
2080 eLp/= uuu[0]+uuu[1]+uuu[2]+uuu[3]+1.e-16;
2081 _mm_storeu_ps(vvv,_ch);
2082 CHR = vvv[0]+vvv[1]+vvv[2]+vvv[3];
2083 ff = CHR>0. ? 1. : -1.;
2084 _ch = _mm_set1_ps(ff);
2085 _gg = rndwave ? _oo : _mm_set1_ps(0.5);
2087 for(j=0; j<m4; j+=4) {
2089 __m128* _pbb = _bb+jf;
2090 __m128* _pBB = _BB+jf;
2091 __m128* _pxp = _xp+jf;
2092 __m128* _pXP = _XP+jf;
2093 __m128* _pxi = _pp+jf;
2094 __m128* _pXI = _PP+jf;
2096 _ee = _mm_sub_ps(_mm_mul_ps(_ch,*_chr),_1);
2108 _coO++; _siO++; _coP++; _siP++; _chr++;
2115 _rr = (__m128*) rrr.
data;
2117 Lo = Co = Eo = Lr = Cr = Do = 0.;
2118 for(j=0; j<m4; j+=4) {
2121 __m128* _pbb = _bb+jf;
2122 __m128* _pBB = _BB+jf;
2123 __m128* _pxi = _pp+jf;
2124 __m128* _pXI = _PP+jf;
2125 __m128* _pww = _ww+jf;
2126 __m128* _pWW = _WW+jf;
2138 _ll = _mm_add_ps(_ll,_LL);
2139 _ec = _mm_add_ps(_ec,_EC);
2140 _ed = _mm_add_ps(_ed,_ED);
2141 _ee = _mm_add_ps(_ee,_EE);
2143 _mm_storeu_ps(vvv,_ee);
2144 Eo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2145 _mm_storeu_ps(vvv,_ec);
2146 Co += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2147 _mm_storeu_ps(vvv,_ed);
2148 Do += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2149 _mm_storeu_ps(vvv,_ll);
2150 Lo += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2154 _en = _mm_andnot_ps(_sm,_mm_sub_ps(_ee,_ll));
2155 _en = _mm_add_ps(_en,*_rr); _rr++;
2156 _EC = _mm_andnot_ps(_sm,_ec);
2157 _cc = _mm_add_ps(_mm_add_ps(_EC,_ed),_en);
2158 _cc = _mm_div_ps(_mm_sub_ps(_ec,_ed),_cc);
2160 _mm_storeu_ps(vvv,_mm_mul_ps(_ll,_cc));
2161 Lr += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2162 _mm_storeu_ps(vvv,_mm_mul_ps(_ec,_cc));
2163 Cr += vvv[0]+vvv[1]+vvv[2]+vvv[3];
2173 aa = Eo>0. ? Lo/Eo : 0.;
2174 AA = Eo>0. ? Lr/Eo : 0.;
2175 if(!mra) skyProb.data[
l] = this->
delta<0 ? aa : AA;
2179 ff = FF = Et = Nm = 0.;
2180 for(j=0; j<
m; j++) {
2191 Nm = Et>0.&&Nm>0 ? Et/Nm : 0.;
2192 ff = Eo>0. ? 2*ff/Eo : 0.;
2193 FF = Eo>0. ? 2*FF/Eo : 0.;
2196 if(ID==
id && !mra) {
2197 Eo += 0.001; Cr += 0.001;
2198 this->nAntenaPrior.set(l, sqrt(ff+FF));
2199 this->nSensitivity.set(l, sqrt(ff+FF));
2200 this->nAlignment.set(l, sqrt(FF/ff));
2201 this->nLikelihood.set(l, Lo/Eo);
2202 this->nNullEnergy.set(l, (Eo-Lo)/Eo);
2203 this->nCorrEnergy.set(l, Cr/Eo);
2204 this->nCorrelation.set(l, Ln);
2205 this->nSkyStat.set(l, AA);
2206 this->nProbability.set(l, skyProb.data[l]);
2207 this->nDisbalance.set(l, 2*Do/Eo);
2208 this->nEllipticity.set(l, eLp);
2209 this->nPolarisation.set(l, -atan2(s2p,c2p)*180./
PI/4.);
2210 this->nNetIndex.set(l, Nm);
2213 if(prior && !mra && ID!=
id) {
2215 for(j=0; j<
m; j++) {
2222 ff = Eo>0. ? 2*ff/Eo : 0.;
2223 FF = Eo>0. ? 2*FF/Eo : 0.;
2224 this->nAntenaPrior.set(l, sqrt(ff+FF));
2227 if(AA>STAT && !mra) {STAT=AA; lm=
l; Vm=
m;}
2230 if(STAT==0. || (mra && AA<=0.)) {
2231 pwc->
sCuts[
id-1]=1; count=0;
2232 pwc->
clean(
id);
continue;
2235 if(le-lb) {lb=le=lm;
goto optsky;}
2237 double Em_all,Ln_all,Ls_all,Lr_all;
2238 double Eo_all,Lo_all,Co_all,Do_all,cc_all;
2242 Em_all=Ls_all=Ln_all = 0;
2248 pwc->
cData[
id-1].skySize =
m;
2249 pwc->
cData[
id-1].likesky = Lo;
2250 vint = &(pwc->
cList[
id-1]);
2251 for(j=0; j<vint->size(); j++) {
2257 for(j=0; j<Vm; j++) {
2265 if(ee-em>Es) Ln_all += ee;
2266 GNoise += rrr.
data[
j];
2269 pwc->
cData[
id-1].skyStat = Lr/Eo;
2272 Ns = Eo_all-Lo_all+Do+GNoise;
2273 gg = Ls_all*Ln_all/Em_all;
2274 pwc->
cData[
id-1].SUBNET = gg/(
fabs(gg)+Ns);
2276 mra=
true;
goto optsky;
2279 if(AA<this->
netCC || !m) {
2280 pwc->
sCuts[
id-1]=1; count=0;
2281 pwc->
clean(
id);
continue;
2290 double Em_mra,Ln_mra,Ls_mra;
2291 double Eo_mra,Lo_mra,Co_mra;
2293 M=
m; m=0; GNoise=0.;
2294 Em_mra=Ln_mra=Ls_mra = 0;
2295 for(j=0; j<
M; j++) {
2299 float* pco = coORT.
data+
j;
2300 __m128* _pxi = _xi+jf;
2301 __m128* _pXI = _XI+jf;
2306 if(ee-em>Es) Ln_mra += ee;
2307 GNoise += rrr.
data[
j];
2314 for(i=0; i<
nIFO; i++) {
2323 pwc->
sCuts[
id-1]=1; count=0;
2324 pwc->
clean(
id);
continue;
2327 Em=Eo; Lm=Lo; Do*=2;
2328 Eo_mra=Eo; Lo_mra=Lo; Co_mra=Co;
2330 pwc->
cData[
id-1].netcc = Lr/Eo;
2331 Nc = Eo-Lo+Do/2+GNoise;
2332 gg = Ls_mra*Ln_mra/Em_mra;
2333 pwc->
cData[
id-1].subnet = gg/(
fabs(gg)+Nc);
2334 pwc->
cData[
id-1].skycc = Co/(
fabs(Co)+Nc);
2345 NETX (vtof->push_back(ml[0][lm]); ,
2346 vtof->push_back(ml[1][lm]); ,
2347 vtof->push_back(ml[2][lm]); ,
2348 vtof->push_back(ml[3][lm]); ,
2349 vtof->push_back(ml[4][lm]); ,
2350 vtof->push_back(ml[5][lm]); ,
2351 vtof->push_back(ml[6][lm]); ,
2352 vtof->push_back(ml[7][lm]); )
2355 if((wfsave)||(mdcListSize() && !lag)) {
2356 int m0d = mureana ? 0 : 1;
2357 if(this->getMRAwave(
id,lag,
'S',m0d,
true)) {
2359 for(i=0; i<
nIFO; i++) {
2360 pd = this->getifo(i);
2361 pd->
RWFID.push_back(
id);
2365 pd->
RWFP.push_back(wf);
2368 if(this->getMRAwave(
id,lag,
's',m0d,
true)) {
2370 for(i=0; i<
nIFO; i++) {
2371 pd = this->getifo(i);
2372 pd->
RWFID.push_back(-
id);
2376 pd->
RWFP.push_back(wf);
2381 Lo = Eo = To = Fo = No = 0.;
2382 for(i=0; i<
nIFO; i++) {
2387 int two = mureana ? 1 : 2;
2388 int m0d = mureana ? 0 : -1;
2390 this->getMRAwave(
id,lag,
'W',m0d);
2391 this->getMRAwave(
id,lag,
'S',m0d);
2392 for(i=0; i<
nIFO; i++) {
2396 float sSNR = d->
get_SS()/two;
2397 float xSNR = d->
get_XS()/two;
2398 float null = d->
get_NN()/two;
2399 float enrg = d->
get_XX()/two;
2420 pwc->
cData[
id-1].likenet = Lo;
2421 pwc->
cData[
id-1].energy = Eo;
2422 pwc->
cData[
id-1].enrgsky = Eo_all;
2423 pwc->
cData[
id-1].netecor = Co;
2424 pwc->
cData[
id-1].netnull = No+GNoise/2;
2425 pwc->
cData[
id-1].netED = Do;
2426 pwc->
cData[
id-1].netRHO = sqrt(Co*cc_all/(nIFO-1.));
2427 pwc->
cData[
id-1].netrho = sqrt(Cr/(nIFO-1.));
2428 pwc->
cData[
id-1].cTime = To;
2429 pwc->
cData[
id-1].cFreq = Fo;
2430 pwc->
cData[
id-1].theta = nLikelihood.getTheta(lm);
2431 pwc->
cData[
id-1].phi = nLikelihood.getPhi(lm);
2432 pwc->
cData[
id-1].gNET = sqrt(ff+FF);
2433 pwc->
cData[
id-1].aNET = sqrt(FF/ff);
2434 pwc->
cData[
id-1].iNET = Nm;
2435 pwc->
cData[
id-1].iota = Ns;
2436 pwc->
cData[
id-1].psi = -atan2(s2p,c2p)*180./
PI/4.;
2437 pwc->
cData[
id-1].ellipticity = eLp;
2441 if(sqrt(Co/(nIFO-1.))<this->
netRHO || pwc->
cData[
id-1].skycc<this->netCC) {
2442 pwc->
sCuts[
id-1]=1; count=0;
2443 pwc->
clean(
id);
continue;
2446 cc = pwc->
cData[
id-1].skycc;
2448 printf(
"id|lm %3d|%6d rho=%4.2f cc: %5.3f|%5.3f|%5.3f|%5.3f \n",
2449 int(
id),
int(lm),sqrt(Co/(nIFO-1)),STAT,cc,pwc->
cData[
id-1].netcc,AA);
2450 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",
2451 nLikelihood.getTheta(l),nLikelihood.getPhi(l),To,Fo,Lo,Lo_mra,Lo_all,Eo,Em,Eo_all);
2452 printf(
" D|N: %4.1f|%4.1f|%4.1f Vm|m=%3d|%3d subnet=%4.3f|%4.3f \n",
2453 Do,No,Nc,
int(Vm),
int(M),pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
2454 hist->Fill(pwc->
cData[
id-1].subnet,pwc->
cData[
id-1].SUBNET);
2460 pwc->
p_Ind[
id-1].push_back(m);
2461 double T = To+pwc->
start;
2462 std::vector<float> sArea;
2463 pwc->
sArea.push_back(sArea);
2464 pwc->
p_Map.push_back(sArea);
2468 double rMs = this->
delta<0 ? 0 : 2;
2469 if(iID<=0 || ID==
id) getSkyArea(
id,lag,T,rMs);
2473 pwc->
cData[
id-1].mchirp = 0;
2474 pwc->
cData[
id-1].mchirperr = 0;
2475 pwc->
cData[
id-1].tmrgr = 0;
2476 pwc->
cData[
id-1].tmrgrerr = 0;
2477 pwc->
cData[
id-1].chi2chirp = 0;
2481 cc = Co_all/(
fabs(Co_all)+ee);
2482 printf(
"mchirp : %d %g %.2e %.3f %.3f %.3f %.3f \n\n",
2483 int(
id),cc,pwc->
cData[
id-1].mchirp,
2484 pwc->
cData[
id-1].mchirperr, pwc->
cData[
id-1].tmrgr,
2485 pwc->
cData[
id-1].tmrgrerr, pwc->
cData[
id-1].chi2chirp);
2488 if(ID==
id && !
EFEC) {
2489 this->nSensitivity.gps =
T;
2490 this->nAlignment.gps =
T;
2491 this->nDisbalance.gps =
T;
2492 this->nLikelihood.gps =
T;
2493 this->nNullEnergy.gps =
T;
2494 this->nCorrEnergy.gps =
T;
2495 this->nCorrelation.gps =
T;
2496 this->nSkyStat.gps =
T;
2497 this->nEllipticity.gps =
T;
2498 this->nPolarisation.gps=
T;
2499 this->nNetIndex.gps =
T;
2502 pwc->
sCuts[
id-1] = -1;
2513 this->wfsave = value.
wfsave;
2514 this->MRA = value.
MRA;
2515 this->nRun = value.
nRun;
2516 this->nLag = value.
nLag;
2518 this->mIFO = value.
mIFO;
2519 this->Step = value.
Step;
2520 this->Edge = value.
Edge;
2522 this->aNET = value.
aNET;
2523 this->iNET = value.
iNET;
2524 this->eCOR = value.
eCOR;
2525 this->e2or = value.
e2or;
2527 this->norm = value.
norm;
2529 this->local = value.
local;
2533 this->gamma = value.
gamma;
2534 this->penalty = value.
penalty;
2537 this->pSigma = value.
pSigma;
2538 this->ifoList = value.
ifoList;
2541 this->NDM.clear(); this->NDM=value.
NDM;
2542 this->ifoList.clear(); this->ifoList=value.
ifoList;
2543 this->ifoName.clear(); this->ifoName=value.
ifoName;
2544 this->wc_List.clear(); this->wc_List=value.
wc_List;
2546 this->mdcList.clear(); this->mdcList=value.
mdcList;
2547 this->livTime.clear(); this->livTime=value.
livTime;
2548 this->mdcTime.clear(); this->mdcTime=value.
mdcTime;
2549 this->mdcType.clear(); this->mdcType=value.
mdcType;
2550 this->mdc__ID.clear(); this->mdc__ID=value.
mdc__ID;
2561 if(ifoList.size()==
NIFO) {
2562 cout <<
"network::add - Error : max number of detectors is " <<
NIFO << endl;
2568 this->ifoList.push_back(d);
2569 this->ifoName.push_back(d->
Name);
2573 for(i=0; i<
n; i++) {
2575 if(i<n-1) this->NDM[
i].push_back(0);
2576 else this->NDM.push_back(v);
2580 return ifoList.size();
2589 int N = ifoListSize();
2592 size_t nL = size_t(Edge*pw->
wrate()*
M);
2593 size_t nR = pw->
size() - nL - 1;
2595 for(
int i=1;
i<
N;
i++) w += getifo(
i)->TFmap;
2596 double amp,
avr, bbb, alp;
2599 for(
int i=nL;
i<nR;
i++) {
2600 amp = (double)w.
data[
i];
2601 if(amp>N*100) amp = N*100.;
2602 if(amp>0.001) {avr+=amp; bbb+=
log(amp); nn++;}
2605 alp =
log(avr)-bbb/nn;
2606 alp = (3-alp+sqrt((alp-3)*(alp-3)+24*alp))/12./alp;
2609 return avr*
iGamma(alp,bbb)/alp/2;
2618 int N = ifoListSize();
2621 size_t nL = size_t(Edge*pw->
wrate()*
M);
2622 size_t nR = pw->
size() - nL;
2624 for(
int i=1;
i<
N;
i++) w += getifo(
i)->TFmap;
2628 double v10 = w.
waveSplit(nL,nR,nR-
int(p10*fff*(nR-nL)));
2630 double med = w.
waveSplit(nL,nR,nR-
int(0.2*fff*(nR-nL)));
2632 while(p00<0.2) {p00 = 1-
Gamma(N*m,med); m+=0.01;}
2634 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");
2635 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",
2636 m,M,p,med,
iGamma(N*m,0.2),v10,
iGamma(N*m,p10),val,
iGamma(N*m,p),N*
log(m),fff);
2637 return (
iGamma(N*m,p)+val)*0.3+N*
log(m);
2647 size_t K = ifoListSize();
2649 if(getifo(0) && t>0. && K) {
2650 I = getifo(0)->TFmap.maxLayer()+1;
2651 n = 2*t*getifo(0)->TFmap.rate()/I + 1;
2652 if(!getifo(0)->TFmap.size()) n = 1.;
2654 else if(t < 0.) n = -
t;
2655 return sqrt(
iGamma1G(K/2.,1.-(p/2.)/pow(n,K/2.))/K);
2660 int iTYPE = this->MRA ? 0 : 1;
2673 cout<<
"wrong size "<<cid.
size()<<
" "<<rho.
size()<<endl;
2675 for(
size_t i=0;
i<cid.
size();
i++){
2676 printf(
"%2d %5.0f vol=%4.0f size=%4.0f like=%5.1e rho=%5.1f ",
2677 int(n),cid[
i],vol[i],siz[i],lik[i],rho[i]);
2678 printf(
"sub=%3.2f rate=%4.0f time=%8.3f To=%8.3f freq=%5.0f\n",
2679 sub[i],rat[i],tim[i],T_o[i],frq[i]);
2691 skymap temp(sms,t1,t2,p1,p2);
2692 size_t m = temp.
size();
2693 size_t n = this->ifoList.size();
2695 nSensitivity = temp;
2697 nCorrelation = temp;
2705 nEllipticity = temp;
2706 nProbability = temp;
2707 nPolarisation= temp;
2708 nAntenaPrior = temp;
2710 for(i=0; i<
n; i++) {
2712 d->
setTau(sms,t1,t2,p1,p2);
2717 skyMask.resize(m); skyMask = 1;
2718 skyMaskCC.resize(0);
2719 skyHole.resize(m); skyHole = 1.;
2731 skymap temp(healpix_order);
2732 size_t m = temp.
size();
2733 size_t n = this->ifoList.size();
2735 nSensitivity = temp;
2737 nCorrelation = temp;
2745 nEllipticity = temp;
2746 nProbability = temp;
2747 nPolarisation= temp;
2748 nAntenaPrior = temp;
2750 for(i=0; i<
n; i++) {
2752 d->
setTau(healpix_order);
2757 skyMask.resize(m); skyMask = 1;
2758 skyMaskCC.resize(0);
2759 skyHole.resize(m); skyHole = 1.;
2770 size_t N = this->ifoList.size();
2777 if(strstr(frame,
"FL") || strstr(frame,
"FS")) {
2778 tm = strstr(frame,
"FS") ? 1. : 0.;
2779 gg = strstr(frame,
"FS") ? 1. : -1.;
2782 for(n=0; n<
N; n++) {
2783 for(m=n+1; m<
N; m++) {
2784 s = ifoList[
n]->tau;
2785 s -= ifoList[
m]->tau;
2786 t = gg*(s.
max()-s.
min());
2787 if(t < tm) { tm=
t; nn =
n; mm =
m; }
2791 s = ifoList[nn]->tau;
2792 s+= ifoList[mm]->tau;
2797 else if(strstr(frame,
"BC")) {
2798 for(n=1; n<
N; n++) s += ifoList[n]->
tau;
2804 for(n=0; n<
N; n++) {
2805 if(strstr(frame,getifo(n)->Name)) this->mIFO = n;
2807 s = ifoList[this->mIFO]->tau;
2811 for(n=0; n<
N; n++) ifoList[n]->
tau -= s;
2820 size_t n = this->ifoList.size();
2821 double maxTau = -1.;
2825 if(n < 2)
return 0.;
2828 for(i=0; i<
n; i++) {
2829 tmax = ifoList[
i]->tau.max();
2830 tmin = ifoList[
i]->tau.min();
2831 if(tmax > maxTau) maxTau = tmax;
2832 if(tmin < minTau) minTau = tmin;
2834 if(strstr(name,
"min"))
return minTau;
2835 if(strstr(name,
"max"))
return maxTau;
2836 if(strstr(name,
"MAX"))
return fabs(maxTau)>
fabs(minTau) ?
fabs(maxTau) :
fabs(minTau);
2837 return (maxTau-minTau)/2.;
2875 size_t M = this->ifoList.size();
2876 if(M >
NIFO)
return;
2880 for(
size_t m=0;
m<
M;
m++) {
2881 D[
m] = this->getifo(
m);
2882 if(D[
m]->mFp.size() != D[0]->
mFp.
size()) {
2883 cout<<
"network::setIndex(): invalid detector skymaps\n";
2886 this->setAntenna(D[
m]);
2901 size_t N = ifoList.size();
2908 cout<<
"network::setDelayIndex(): invalid network\n";
2913 for(n=0; n<
N; n++) dr[n] = ifoList[n];
2921 for(n=0; n<
N; n++) {
2932 this->nPenalty = dr[0]->
tau;
2933 this->nNetIndex = dr[0]->
tau;
2942 for(n=0; n<
N; n++) {
2943 for(m=0; m<
N; m++) {
2945 i = t>0 ?
int(t*rTDF+0.5) :
int(t*rTDF-0.5);
2951 for(n=0; n<
N; n++) {
2953 for(m=0; m<
N; m++) {
2954 for(k=0; k<
N; k++) {
2955 t =
fabs(mm[n][k]-mm[n][m]-tt[m][k]);
2956 if(TT[n] < t) TT[
n] =
t;
2962 for(m=0; m<
N; m++) {
2963 if(t>TT[m]) { t = TT[
m]; k =
m; }
2965 this->nPenalty.set(l,
double(t));
2968 if(mIFO<9) i = mm[
k][this->mIFO];
2969 else i = t>0 ?
int(t*rTDF+0.5) :
int(t*rTDF-0.5);
2975 for(m=0; m<
N; m++) {
3000 size_t M = this->ifoList.size();
3006 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3008 int type = this->
optim ? R : TYPE;
3013 vector<wavearray<double> > snr(M);
3014 vector<wavearray<double> > nul(M);
3016 double xcut,xnul,
xsnr,amin,emin,like,LIKE;
3019 for(j=0; j<nLag; j++) {
3021 if(!this->wc_List[j].
size())
continue;
3023 cid = this->wc_List[
j].
get((
char*)
"ID",0,
'S',type);
3027 siz = this->wc_List[
j].
get((
char*)
"size",0,
'S',type);
3028 SIZ = this->wc_List[
j].
get((
char*)
"SIZE",0,
'S',type);
3030 for(m=0; m<
M; m++) {
3031 snr[
m] = this->wc_List[
j].get((
char*)
"energy",m+1,
'S',type);
3032 nul[
m] = this->wc_List[
j].get((
char*)
"null",m+1,
'W',type);
3035 for(k=0; k<
K; k++) {
3037 if(siz.
data[k] < core_size)
continue;
3038 i = ID = size_t(cid.
data[k]+0.1);
3040 if(tYPe==
'i' || tYPe==
's' || tYPe==
'g' || tYPe==
'r' ||
3041 tYPe==
'I' || tYPe==
'S' || tYPe==
'G' || tYPe==
'R') {
3042 skip = !this->SETNDM(i,j,
true,type);
3044 else if(tYPe==
'b' || tYPe==
'B' || tYPe==
'E') {
3045 skip = !this->setndm(i,j,
true,type);
3052 cout<<
"network::netcut: - undefined NDM matrix"<<endl;
3053 this->wc_List[
j].ignore(ID);
3054 count += size_t(SIZ.
data[k]+0.1);
3058 xnul = xsnr = like = LIKE = 0;
3059 emin = amin = 1.e99;
3060 for(n=0; n<
M; n++) {
3061 xnul += this->getifo(n)->null;
3062 xsnr += snr[
n].data[
k];
3063 like += snr[
n].data[
k]-this->getifo(n)->null;
3064 xcut = snr[
n].data[
k]-this->getifo(n)->null;
3065 if(xcut < emin) emin = xcut;
3066 xcut = (xcut-this->getifo(n)->null)/snr[n].data[k];
3067 if(xcut < amin) amin = xcut;
3068 for(m=0; m<
M; m++) {
3069 LIKE += this->getNDM(n,m);
3073 if(cut ==
'A' || cut ==
'E') {
3074 emin = amin = 1.e99;
3075 for(n=0; n<
M; n++) {
3076 xcut = snr[
n].data[
k] - this->getifo(n)->null;
3077 if(like-xcut < emin) emin = like-xcut;
3078 xcut = 2*(like-xcut)/(xsnr-xcut) - 1.;
3079 if(xcut < amin) amin = xcut;
3084 if(cut ==
'c' || cut ==
'C') xcut = this->eCOR/(xnul +
fabs(this->eCOR));
3085 if(cut ==
'x' || cut ==
'X') xcut = this->eCOR/sqrt(
fabs(this->eCOR)*M);
3086 if(cut ==
'l' || cut ==
'L') xcut = like;
3087 if(cut ==
'a' || cut ==
'A') xcut = amin;
3088 if(cut ==
'e' || cut ==
'E') xcut = emin;
3089 if(cut ==
'r' || cut ==
'R') xcut = this->eCOR/sqrt((xnul +
fabs(this->eCOR))*M);
3092 this->wc_List[
j].ignore(ID);
3093 count += size_t(SIZ.
data[k]+0.1);
3097 if(cut==
'X') cout<<xcut<<cut<<
": size="<<SIZ.
data[
k]<<
" L="<<LIKE<<
" ID="<<ID<<endl;
3101 printf(
"%3d %4d %7.2e %7.2e %7.2e %7.2e %7.2e %7.2e %7.2e\n",
3102 (
int)n,(
int)i,getNDM(0,0),getNDM(1,1),getNDM(2,2),getNDM(0,1),getNDM(0,2),getNDM(1,2),
3103 getNDM(0,0)+getNDM(1,1)+getNDM(2,2)+2*(getNDM(0,1)+getNDM(0,2)+getNDM(1,2)));
3121 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3122 size_t R = size_t(this->ifoList[0]->
getTFmap()->
rate()/I+0.5);
3123 size_t N = this->ifoList[0]->getTFmap()->size();
3124 size_t M = this->ifoList.size();
3126 int window =
int(T*R+0.5)*I/2;
3130 int inD, frst, last, nB, nT;
3135 std::vector<detector*> pDet; pDet.clear();
3136 std::vector<double*> pDat; pDat.clear();
3138 for(m=0; m<
M; m++) {
3139 pDet.push_back(ifoList[m]);
3140 pDat.push_back(ifoList[m]->
getTFmap()->data);
3143 size_t wmode = pDet[0]->getTFmap()->w_mode;
3144 if(wmode != 1)
return 0;
3151 for(n=0; n<nLag; n++) {
3152 V = this->wc_List[
n].pList.
size();
3153 nB =
int(this->wc_List[n].getbpp()*2*window+0.5);
3157 for(j=0; j<V; j++) {
3159 pix = &(this->wc_List[
n].pList[
j]);
3160 if(R !=
size_t(pix->
rate+0.5))
continue;
3164 for(m=0; m<
M; m++) {
3166 inD = size_t(pix->
getdata(
'I',m)+0.1);
3171 if(inD-window < offset) {
3173 last = frst + 2*window;
3175 else if(inD+window >
int(N-offset)) {
3177 frst = last - 2*window;
3184 if( inD>
int(N-offset) || inD<offset ||
3185 frst>
int(N-offset) || frst<offset ||
3186 last>
int(N-offset) || last<offset) {
3187 cout<<
"network::setRank() error\n";
3193 rank = pDet[
m]->getTFmap()->getSampleRankE(inD,frst,last);
3196 if(rank > nT-nB) rank =
log(nB/
double(nT+1-rank));
3215 size_t L = this->skyHole.size();
3224 cout<<endl<<
"network::setSkyMask() - skymap size L=0"<<endl<<endl;
3227 cout<<endl<<
"network::setSkyMask() - NULL input skymask file"<<endl<<endl;
3229 }
else if(!strlen(file)) {
3230 cout<<endl<<
"network::setSkyMask() - input skymask file not defined"<<endl<<endl;
3232 }
else if( (in=fopen(file,
"r"))==NULL ) {
3233 cout << endl <<
"network::setSkyMask() - input skymask file '" 3234 << file <<
"' not exist" << endl << endl;;
3238 while(fgets(str,1024,in) != NULL){
3240 if(str[0] ==
'#')
continue;
3241 if((pc = strtok(str,
" \t")) == NULL)
continue;
3242 if(pc) i = atoi(pc);
3243 if((pc = strtok(NULL,
" \t")) == NULL)
continue;
3244 if(pc && i>=0 && i<
int(L)) {
3245 this->skyHole.data[
i] = atof(pc);
3246 this->nSkyStat.set(i, atof(pc));
3249 cout<<endl<<
"network::setSkyMask() - " 3250 <<
"skymask file contains index > max L="<<L<<endl<<endl;
3255 cout<<endl<<
"network::setSkyMask() - " 3256 <<
"the number of indexes in the skymask file != L="<<L<<endl<<endl;
3259 a = skyHole.mean()*skyHole.size();
3260 skyHole *= a>0. ? 1./
a : 0.;
3261 if(f==0.) { skyHole = 1.;
return n; }
3263 double* p = this->skyHole.data;
3264 double** pp = (
double **)malloc(L*
sizeof(
double*));
3265 for(l=0; l<
L; l++) pp[l] = p + l;
3267 skyProb.waveSort(pp,0,L-1);
3270 for(l=0; l<
L; l++) {
3272 *pp[
l] = a/L<f ? 0. : 1.;
3273 if(*pp[l] == 0.) this->nSkyStat.set(pp[l]-p,*pp[l]);
3283 size_t L = this->skyHole.size();
3290 if(skycoord!=
'e' && skycoord!=
'c') {
3291 cout <<
"network::setSkyMask() - wrong input sky coordinates " 3292 <<
" must be 'e'/'c' earth/celestial" << endl;;
3297 cout<<endl<<
"network::setSkyMaskCC() - skymap size L=0"<<endl<<endl;
3300 cout<<endl<<
"network::setSkyMaskCC() - NULL input skymask file"<<endl<<endl;
3302 }
else if(!strlen(file)) {
3303 cout<<endl<<
"network::setSkyMaskCC() - input skymask file not defined"<<endl<<endl;
3305 }
else if( (in=fopen(file,
"r"))==NULL ) {
3306 cout << endl <<
"network::setSkyMaskCC() - input skymask file '" 3307 << file <<
"' not exist" << endl << endl;;
3311 if(skycoord==
'e') {skyMask.resize(L); skyMask = 1;}
3312 if(skycoord==
'c') {skyMaskCC.resize(L); skyMaskCC = 1;}
3314 while(fgets(str,1024,in) != NULL){
3316 if(str[0] ==
'#')
continue;
3317 if((pc = strtok(str,
" \t")) == NULL)
continue;
3318 if(pc) i = atoi(pc);
3319 if((pc = strtok(NULL,
" \t")) == NULL)
continue;
3320 if(pc && i>=0 && i<
int(L)){
3322 if(skycoord==
'e') this->skyHole.data[
i]=data;
3323 if(skycoord==
'c') this->skyMaskCC.data[
i]=data;
3326 cout<<endl<<
"network::setSkyMask() - " 3327 <<
"skymask file contains index > max L="<<L<<endl<<endl;
3332 cout<<endl<<
"network::setSkyMask() - " 3333 <<
"the number of indexes in the skymask file != L="<<L<<endl<<endl;
3344 if(skycoord!=
'e' && skycoord!=
'c') {
3345 cout <<
"network::setSkyMask() - wrong input sky coordinates " 3346 <<
" must be 'e'/'c' earth/celestial" << endl;;
3350 size_t L = this->skyHole.size();
3351 if((
int)sm.
size()!=
L) {
3352 cout <<
"network::setSkyMask() - wrong input skymap size " 3353 << sm.
size() <<
" instead of " << L << endl;;
3359 for(
int i=0;
i<
L;
i++) this->skyHole.data[
i]=sm.
get(
i);
3362 skyMaskCC.resize(L);
3363 for(
int i=0;
i<
L;
i++) this->skyMaskCC.data[
i]=sm.
get(
i);
3381 std::map <string, int> mdcMap;
3383 if( (in=fopen(file,
"r"))==NULL ) {
3384 cout<<
"network::readMDClog() - no file is found \n";
3388 while(fgets(str,1024,in) != NULL){
3390 if(str[0] ==
'#')
continue;
3395 if((p = strtok(STR,
" \t")) == NULL)
continue;
3397 for(i=1; i<nTime; i++) {
3398 p = strtok(NULL,
" \t");
3404 if(gps==0. ||
fabs(GPS-gps)<7200.) {
3405 this->mdcList.push_back(str);
3406 this->mdcTime.push_back(GPS);
3412 if((p = strtok(str,
" \t")) == NULL)
continue;
3414 for(i=1; i<nName; i++) {
3415 p = strtok(NULL,
" \t");
3419 if(p)
if(mdcMap.find(p)==mdcMap.end()) mdcMap[p]=imdcMap++;
3424 this->mdcType.resize(mdcMap.size());
3425 std::map<std::string, int>::iterator iter;
3426 for (iter=mdcMap.begin(); iter!=mdcMap.end(); iter++) {
3427 this->mdcType[iter->second]=iter->first;
3430 for(
int j=0;j<this->mdcType.size();j++) {
3433 else if(j<10000) step=100;
3436 printf(
"type %3d\t",(
int)j);
3437 cout<<
" has been assigned to waveform "<<mdcType[
j]<<endl;
3441 return this->mdcList.size();
3454 if( (in=fopen(file,
"r"))==NULL ) {
3455 cout<<
"network::readSEGlist(): specified segment file "<<file<<
" does not exist\n";
3459 while(fgets(str,1024,in) != NULL){
3461 if(str[0] ==
'#')
continue;
3465 if((p = strtok(str,
" \t")) == NULL)
continue;
3467 for(i=1; i<
n; i++) {
3468 p = strtok(NULL,
" \t");
3474 SEG.
start = atof(p);
3475 p = strtok(NULL,
" \t");
3497 this->mdc__ID.clear();
3499 size_t I = this->ifoList.
size();
3510 size_t L = this->mdcList.size();
3511 size_t n = size_t(this->Edge*R+0.5);
3512 int W =
int(Tw*R/2.+0.5);
3515 if(!I || !N)
return 0.;
3517 if(this->veto.size() != size_t(N)) {
3518 this->veto.resize(N);
3523 for(k=0; k<
K; k++) {
3534 for(j=jb; j<je; j++) this->veto.data[j] = 1;
3537 if(!K) this->veto = 1;
3539 for(k=0; k<
L; k++) {
3541 if(gps == 0.)
continue;
3545 for(i=0; i<I; i++) {
3546 d = this->ifoList[
i];
3554 jm =
int((gps-S)*R);
3555 jb = jm-W; je = jm+W;
3557 if(jb >=N)
continue;
3559 if(je <=0)
continue;
3560 if(je-jb <
int(R))
continue;
3561 if(jm<jb || jm>je)
continue;
3563 for(j=jb; j<je; j++) w.
data[j] = 1;
3565 if(veto.data[jm]) this->mdc__ID.push_back(k);
3568 if(L) this->veto *=
w;
3570 for(k=n; k<N-
n; k++) live+=this->veto.
data[k];
3581 size_t M = this->ifoList.size();
3595 for(j=0; j<
id.size(); j++) {
3596 if(
id.data[j] == ID) { flag=
true;
break; }
3599 if(!flag)
return false;
3602 v = wc->
nTofF[ID-1];
3604 l = size_t(
log(k*1.)/
log(2.)+0.1);
3608 for(i=0; i<
M; i++) {
3609 pd = this->getifo(i);
3610 if(pd->
getwave(ID,*wc,atype,i) == 0.) { flag=
false;
break; }
3659 if(m > (
int)w.
size()) m = w.
size();
3674 size_t nIFO = this->ifoList.size();
3679 bool signal = (abs(atype)==
'W' || abs(atype)==
'w') ?
false :
true;
3682 for(j=0; j<
id.size(); j++) {
3683 if(
size_t(
id.data[j]+0.1) ==
ID) flag=
true;
3685 if(!flag)
return false;
3690 v = pwc->
nTofF[ID-1];
3694 for(i=0; i<
nIFO; i++) {
3697 if(x.
size() == 0.) {cout<<
"zero length\n";
return false;}
3702 double R = this->rTDF;
3703 double tShift = -v[
i]/R;
3708 for (
int ii=0;ii<(
int)x.
size()/2;ii++) {
3709 TComplex X(x.
data[2*ii],x.
data[2*ii+1]);
3710 X=X*C.Exp(TComplex(0.,-2*
PI*ii*df*tShift));
3711 x.
data[2*ii]=X.Re();
3712 x.
data[2*ii+1]=X.Im();
3717 if(signal) this->getifo(i)->waveForm =
x;
3718 else this->getifo(i)->waveBand =
x;
3734 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3735 size_t R = size_t(this->ifoList[0]->
getTFmap()->
rate()/I+0.5);
3736 size_t N = this->ifoList[0]->getTFmap()->size();
3737 size_t M = this->ifoList.size();
3738 size_t jB = size_t(this->Edge*R)*I;
3742 std::vector<detector*> pDet; pDet.clear();
3743 std::vector<double*> pDat; pDat.clear();
3744 std::vector<int> pLag; pLag.clear();
3748 pix.
rate = float(R);
3752 for(m=0; m<
M; m++) {
3753 pDet.push_back(ifoList[m]);
3754 pDat.push_back(ifoList[m]->
getTFmap()->data);
3755 pLag.push_back(
int(ifoList[m]->sHIFt*R*I+0.5));
3758 size_t il = size_t(2.*pDet[0]->TFmap.getlow()/R);
3759 size_t ih = size_t(2.*pDet[0]->TFmap.gethigh()/R);
3760 if(ih==0 || ih>=I) ih = I;
3762 size_t J = size_t(sTARt*R+0.1);
3763 size_t K = size_t(duration*R+0.1);
3766 this->wc_List[0].clear();
3770 for(i=il; i<ih; i++){
3772 S = pDet[0]->TFmap.getSlice(i);
3778 cout<<
"network::initwc() error - index out of limit \n";
3784 for(m=0; m<
M; m++) {
3785 k = pix.
time+pLag[
m];
3787 if(!this->veto.data[k]) save =
false;
3792 pix.
setdata(pDet[m]->getNoise(i,k),
'N',m);
3795 if(save) this->wc_List[0].append(pix);
3800 wc_List[0].start = pDet[0]->TFmap.start();
3801 wc_List[0].stop = N/R/I;
3802 wc_List[0].rate = pDet[0]->TFmap.rate();
3805 this->wc_List[0].pList[npix-1].neighbors[0]=0;
3806 this->wc_List[0].cluster();
3817 size_t nIFO = this->ifoList.size();
3819 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
3820 cout<<
"network::coherence(): \n" 3821 <<
"invalid number of detectors or\n" 3822 <<
"delay filter is not set\n";
3825 if(getifo(0)->getTFmap()->w_mode != 1) {
3826 cout<<
"network::coherence(): invalid whitening mode.\n";
3830 if(factor > 1.) factor = float(nIFO-1)/
nIFO;
3835 size_t jS,jj,nM,jE,LL;
3837 double R = this->ifoList[0]->getTFmap()->rate();
3838 size_t N = this->ifoList[0]->getTFmap()->size();
3839 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
3840 size_t M = this->
filter.size()/I;
3841 size_t K = this->
filter[0].index.size();
3843 size_t jB = size_t(this->Edge*R/I)*I;
3844 double band = this->ifoList[0]->TFmap.rate()/I/2.;
3845 slice S = getifo(0)->getTFmap()->getSlice(0);
3857 double* pdata[
NIFO];
3859 for(n=0; n<
NIFO; n++) {
3861 if(n >= nIFO)
continue;
3862 pdata[
n] = getifo(n)->getTFmap()->
data;
3877 for(l=0; l<
L; l++)
if(skyMask.data[l]) LL++;
3878 for(n=0; n<
NIFO; n++) {
3880 if(n>=nIFO) inDEx[
n] = 0;
3883 ina = this->getifo(n)->index.
data;
3884 for(l=0; l<
L; l++) {
3885 if(skyMask.data[l]) inDEx[
n].
data[k++] = ina[
l];
3891 short* m0 = inDEx[0].data; ,
3892 short*
m1 = inDEx[1].
data; ,
3893 short*
m2 = inDEx[2].
data; ,
3894 short* m3 = inDEx[3].
data; ,
3895 short* m4 = inDEx[4].
data; ,
3896 short* m5 = inDEx[5].
data; ,
3897 short* m6 = inDEx[6].
data; ,
3898 short* m7 = inDEx[7].
data; )
3901 double* pF = (
double*)malloc(K*M*
sizeof(
double));
3902 int* pJ = (
int*)malloc(K*M*
sizeof(
int));
3911 size_t KK = this->
filter.size()/I;
3912 for(n=0; n<
nIFO; n++) {
3913 t1 = getifo(n)->tau.min()*R*getifo(0)->nDFS/I;
3914 t2 = getifo(n)->tau.max()*R*getifo(0)->nDFS/I;
3915 M1[
n] = short(t1+KK/2+0.5)-1;
3916 M2[
n] = short(t2+KK/2+0.5)+1;
3926 this->pixeLHood = getifo(0)->TFmap;
3927 this->pixeLHood = 0.;
3930 if(this->veto.size() !=
N) { veto.
resize(N); veto = 1; }
3934 for(k=0; k<nLag; k++) {
3935 this->wc_List[
k].clear();
3936 this->livTime[
k] = 0.;
3937 this->wc_List[
k].setlow(getifo(0)->TFmap.getlow());
3938 this->wc_List[
k].sethigh(getifo(0)->TFmap.gethigh());
3941 for(i=1; i<I; i++) {
3945 if(a >= getifo(0)->TFmap.gethigh())
continue;
3947 if(a <= getifo(0)->TFmap.getlow())
continue;
3955 pF[k+m*
K] = double(pv->
value[k]);
3960 S = getifo(0)->getTFmap()->getSlice(i);
3963 for(n=0; n<
NIFO; n++) {
3965 emax[
n] = 0; in[
n] = 0;
3975 for(n=0; n<
nIFO; n++) {
3979 for(j=jS; j<
N; j+=I) {
3981 emax[
n].
data[jj] = 0.;
3989 for(m=M1[n]; m<M2[
n]; m++){
3995 emax[
n].
data[jj++] = b;
4003 for(k=0; k<nLag; k++) {
4008 for(n=0; n<
nIFO; n++) {
4009 b = this->getifo(n)->lagShift.data[
k];
4010 if(a>b) { a = b; nM =
n; }
4013 for(n=0; n<
nIFO; n++) {
4014 b = this->getifo(n)->lagShift.data[
k];
4015 in[
n] = (
int((b-a)*R)+jB)/I - 1;
4018 for(jj=jB/I; jj<NN; jj++) {
4021 for(n=0; n<
nIFO; n++) {
4022 if((++in[n]) >=
int(NN)) in[n] -= jE;
4023 m += this->veto.data[inTF.
data[in[
n]]];
4024 E += emax[
n].
data[in[
n]];
4026 this->livTime[
k] += float(m/nIFO);
4027 if(E<Eo || m<nIFO)
continue;
4030 for(n=0; n<
nIFO; n++) {
4031 b = E - emax[
n].
data[in[
n]];
4032 if(b<Eo*factor) skip =
true;
4038 for(n=0; n<
nIFO; n++) {
4041 pq = pdata[
n]+inTF.
data[in[
n]];
4043 for(m=M1[n]; m<M2[
n]; m++){
4051 for(l=0; l<LL; l++) {
4054 pet+=pe[0][m0[l]]; ,
4055 pet+=pe[1][m1[
l]]; ,
4056 pet+=pe[2][m2[
l]]; ,
4057 pet+=pe[3][m3[
l]]; ,
4058 pet+=pe[4][m4[
l]]; ,
4059 pet+=pe[5][m5[
l]]; ,
4060 pet+=pe[6][m6[
l]]; ,
4061 pet+=pe[7][m7[
l]]; )
4063 { skip =
false;
break; }
4069 j = inTF.
data[in[nM]];
4070 pix.
rate = float(R/I);
4072 pix.
core = E>Es ? true :
false;
4076 for(n=0; n<
nIFO; n++) {
4081 wc_List[
k].append(pix);
4082 if(!k) this->pixeLHood.data[
j] = sqrt(E/nIFO);
4089 for(k=0; k<nLag; k++) {
4090 a = getifo(0)->getTFmap()->start();
4091 b = getifo(0)->getTFmap()->size()/R;
4092 this->wc_List[
k].start =
a;
4093 this->wc_List[
k].
stop = a+b;
4094 this->wc_List[
k].rate = R;
4095 this->livTime[
k] *= double(I)/(Io*R);
4098 if(nZero) this->setRMS();
4116 size_t N = this->wc_List[lag].csize();
4117 size_t M = this->mdc__IDSize();
4118 size_t L = this->skyProb.size();
4119 size_t Lm = L-
int(0.9999*L);
4120 size_t nIFO = this->ifoList.size();
4121 bool prior = this->gamma<0?
true:
false;
4136 std::vector<float>* vf = &(this->wc_List[lag].sArea[
id-1]);
4139 double* p = this->skyProb.
data;
4140 double** pp = (
double **)malloc(L*
sizeof(
double*));
4141 for(l=0; l<
L; l++) pp[l] = p + l;
4143 skyProb.waveSort(pp,0,L-1);
4145 double Po = *pp[L-1];
4146 double rms =
fabs(rMs);
4148 double smax=nAntenaPrior.max();
4150 for(l=0; l<
L; l++) {
4151 if(*pp[l] <= 0.) {*pp[
l]=0.;
continue;}
4152 *pp[
l] = exp(-(Po - *pp[l])/2./rms);
4153 if(prior) *pp[
l] *= pow(nAntenaPrior.get(
int(pp[l]-p))/smax,4);
4156 if(prior) skyProb.waveSort(pp,0,L-1);
4158 for(l=0; l<
L; l++) {
4160 nProbability.set(l,p[l]);
4163 if(pOUT) cout<<rMs<<
" "<<*pp[L-1]<<
" "<<*pp[L-2]<<
" "<<*pp[L-3]<<
"\n";
4166 for(m=0; m<11; m++) { v[
m] = 0; vf->push_back(0.); }
4169 for(l=L-1; l>Lm; l--){
4171 for(m=
size_t(vol*10.)+1; m<10; m++) v[m] += 1;
4172 if(vol >= 0.9)
break;
4175 for(m=1; m<10; m++) {
4176 (*vf)[
m] = sqrt(v[m]*s);
4177 if(pOUT && !M) cout<<m<<
" error region: "<<(*vf)[
m]<<endl;
4183 std::vector<float>* vP = &(this->wc_List[lag].p_Map[
id-1]);
4184 std::vector<int>* vI = &(this->wc_List[lag].p_Ind[
id-1]);
4192 if(
nSky<0) {
char spthr[1024];
sprintf(spthr,
"0.%d",
int(abs(
nSky)));pthr=atof(spthr);}
4193 for(l=L-1; l>Lm; l--){
4195 if(
nSky==0 && (K==1000 || sum > 0.99) && K>0)
break;
4196 else if(nSky<0 && sum > pthr && K>0)
break;
4197 else if(
nSky>0 && K==
nSky && K>0)
break;
4199 vI->push_back(
int(pp[l]-p));
4200 vP->push_back(
float(*pp[l]));
4205 if(!M) { free(pp);
return; }
4208 double injTime = 1.e12;
4213 for(m=0; m<
M; m++) {
4214 mdcID = this->getmdc__ID(m);
4215 dT =
fabs(To - this->getmdcTime(mdcID));
4216 if(dT<injTime && INJ.fill_in(
this,mdcID)) {
4219 if(pOUT)
printf(
"getSkyArea: %4d %12.4f %7.3f %f \n",
int(m),To,dT,s);
4223 if(INJ.fill_in(
this,injID)) {
4227 i = this->getIndex(th,ph);
4229 vI->push_back(
int(i));
4230 vP->push_back(
float(p[i]));
4233 for(l=L-1; l>Lm; l--){
4236 if(pp[l]-p ==
int(i))
break;
4238 (*vf)[0] = sqrt(vol);
4243 printf(
"getSkyArea: %5d %12.4f %6.1f %6.1f %6.1f %6.1f %6.2f %6.2f %6.2f %7.5f, %e %d \n",
4244 int(
id),INJ.time[0]-this->getifo(0)->TFmap.start(),INJ.theta[0],INJ.phi[0],
4258 size_t N = this->wc_List[lag].csize();
4259 size_t L = this->skyProb.size();
4260 size_t Lm = L-
int(0.9999*L);
4267 double a,th,
ph,st,sp;
4268 double TH,PH,ST,SP,Eo;
4276 std::vector<float>* vf = &(this->wc_List[lag].sArea[
id-1]);
4279 double* p = this->skyProb.
data;
4280 double** pp = (
double **)malloc(L*
sizeof(
double*));
4281 for(l=0; l<
L; l++) pp[l] = p + l;
4283 skyProb.waveSplit(pp,0,L-1,Lm);
4284 skyProb.waveSort(pp,Lm,L-1);
4287 for(l=L-1; l>=Lm; l--) *pp[l] -= Eo;
4288 for(l=0; l<Lm; l++) *pp[l] = -1.e10;
4293 for(l=Lm; l<
L; l++){
4299 for(in=-1; in<2; in++) {
4301 for(im=-1; im<2; im++) {
4302 i = this->getIndex(th+in*st,ph+im*sp);
4303 if(p[i] >= p[j] || p[i] < -1.e9)
continue;
4309 for(IN=-1; IN<2; IN++) {
4311 for(IM=-1; IM<2; IM++) {
4312 k = this->getIndex(TH+IN*ST,PH+IM*SP);
4313 if(p[i] >=p[k])
continue;
4321 skyProb.waveSplit(pp,0,L-1,Lm);
4322 skyProb.waveSort(pp,Lm,L-1);
4328 for(l=L-2; l>=Lm; l--) {
4330 skyENRG.data[L-l-1] =
a;
4331 if(a < pSigma || vol<1.) {
4338 if(pOUT) cout<<sum/vol<<
" "<<*pp[L-2]<<
" "<<*pp[L-3]<<
" "<<*pp[0]<<
"\n";
4342 for(l=L-1; l>0; l--) {
4343 a = l>Lm ? (L-l-1)*Eo : 30.;
4344 if(a > 30.) a = 30.;
4348 if(pOUT) cout<<
"norm: "<<(skyProb.mean()*
L)<<endl;
4349 skyProb *= 1./(skyProb.mean()*
L);
4350 for(l=0; l<
L; l++) nProbability.set(l,log10(p[l]));
4353 for(m=0; m<11; m++) { v[
m] = 0; vf->push_back(0.); }
4356 for(l=L-1; l>Lm; l--){
4357 for(m=
size_t(sum*10.)+1; m<10; m++) v[m] += 1;
4359 if(sum >= 0.9)
break;
4362 for(m=1; m<10; m++) (*vf)[
m] = sqrt(v[m]*s);
4366 std::vector<float>* vP = &(this->wc_List[lag].p_Map[
id-1]);
4367 std::vector<int>* vI = &(this->wc_List[lag].p_Ind[
id-1]);
4375 if(
nSky<0) {
char spthr[1024];
sprintf(spthr,
"0.%d",
int(abs(
nSky)));pthr=atof(spthr);}
4376 for(l=L-1; l>Lm; l--){
4378 if(
nSky==0 && (K==1000 || sum > 0.99) && K>0)
break;
4379 else if(nSky<0 && sum > pthr && K>0)
break;
4380 else if(
nSky>0 && K==
nSky && K>0)
break;
4382 vI->push_back(
int(pp[l]-p));
4383 vP->push_back(
float(*pp[l]));
4388 size_t M = this->mdc__IDSize();
4390 if(!M) { free(pp);
return; }
4393 double injTime = 1.e12;
4398 for(m=0; m<
M; m++) {
4399 mdcID = this->getmdc__ID(m);
4400 dT =
fabs(To - this->getmdcTime(mdcID));
4401 if(dT<injTime && INJ.fill_in(
this,mdcID)) {
4404 if(pOUT)
printf(
"getSkyArea: %4d %12.4f %7.3f %f \n",
int(m),To,dT,s);
4408 if(INJ.fill_in(
this,injID)) {
4412 i = this->getIndex(th,ph);
4414 vI->push_back(
int(i));
4415 vP->push_back(
float(p[i]));
4418 for(l=L-1; l>Lm; l--){
4421 if(pp[l]-p ==
int(i))
break;
4423 (*vf)[0] = sqrt(vol);
4428 printf(
"getSkyArea: %5d %12.4f %6.1f %6.1f %6.1f %6.1f %6.2f %6.2f %6.2f %7.5f, %e %d \n",
4429 int(
id),INJ.time[0]-this->getifo(0)->TFmap.start(),INJ.theta[0],INJ.phi[0],
4448 size_t nIFO = this->ifoList.size();
4451 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
4452 cout<<
"network::likelihood(): invalid number of detectors or delay filter is not set.\n";
4455 if(type==
'b' || type==
'B' || type==
'E')
4456 return likelihoodB(type, Ao, ID, lag, ind, core);
4458 return likelihoodI(type, Ao, ID, lag, ind, core);
4474 size_t nIFO = this->ifoList.size();
4475 size_t ID = abs(iID);
4476 int N_1 = nIFO>2 ?
int(nIFO)-1 : 2;
4477 int N_2 = nIFO>2 ?
int(nIFO)-2 : 1;
4478 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
4479 cout<<
"network::likelihood(): invalid number of detectors or delay filter is not set.\n";
4487 double Eo = nIFO*Ao*Ao;
4489 double GAMMA = 1.-gamma*gamma;
4491 int LOCAL = local ? 1 : 0;
4494 double gr,gp,
gx,gR,gI,gc,gg,gP,gX,
T,rm,Le,E,LPm;
4495 double STAT,P,EE,
Et,Lo,Lm,
hh,Xp,Xx,XX,Lp,Em,Ep;
4496 double S,co,si,cc,em,um,vm,uc,us,vc,vs,Co,Cp;
4497 NETX(
double c0;,
double c1;,
double c2;,
double c3;,
double c4;,
double c5;,
double c6;,
double c7;)
4500 if(type==
'E' && Eo<nIFO) Eo = double(nIFO);
4501 if(!ID && ind>=0) ind = -1;
4507 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
4508 size_t M = this->
filter.size()/I;
4509 size_t L = ind<0 ? this->
index.
size() : ind+1;
4513 double* pdata[
NIFO];
4514 double* qdata[
NIFO];
4516 for(n=0; n<
nIFO; n++) {
4517 pdata[
n] = getifo(n)->getTFmap()->data;
4518 qdata[
n] = getifo(n)->getTFmap()->data;
4522 std::vector<float>*
F;
4523 std::vector<short>* J;
4533 for(n=0; n<
NIFO; n++) {
4534 fp[
n] = n<nIFO ? getifo(n)->fp.
data : f00.
data;
4535 fx[
n] = n<nIFO ? getifo(n)->fx.data : f00.
data;
4536 ffp[
n] = n<nIFO ? getifo(n)->ffp.data : f00.
data;
4537 ffm[
n] = n<nIFO ? getifo(n)->ffm.data : f00.
data;
4538 fpx[
n] = n<nIFO ? getifo(n)->fpx.data : f00.
data;
4541 short* mm = this->skyMask.data;
4553 std::vector<size_t> pI;
4585 for(n=0; n<
NIFO; n++) {
4594 std::vector<int>* vint;
4595 std::vector<int>* vtof;
4603 this->pixeLHood = getifo(0)->TFmap;
4604 this->pixeLHood = 0.;
4605 this->pixeLNull = getifo(0)->TFmap;
4606 this->pixeLNull = 0.;
4607 this->nSensitivity = 0.;
4608 this->nAlignment = 0.;
4609 this->nNetIndex = 0.;
4610 this->nDisbalance = 0.;
4611 this->nLikelihood = 0.;
4612 this->nNullEnergy = 0.;
4613 this->nCorrEnergy = 0.;
4614 this->nCorrelation = 0.;
4615 this->nSkyStat = 0.;
4616 this->nPenalty = 0.;
4617 this->nProbability = 0.;
4618 this->nAntenaPrior = 0.;
4625 for(n=lag; n<nLag; n++) {
4627 if(!this->wc_List[n].
size())
continue;
4628 logbpp = -
log(this->wc_List[n].getbpp());
4630 cid = this->wc_List[
n].
get((
char*)
"ID",0,
'S',
optim ? R : -R);
4631 cTo = this->wc_List[
n].
get((
char*)
"time",0,
'L',
optim ? R : -R);
4635 for(k=0; k<
K; k++) {
4637 id = size_t(cid.
data[k]+0.1);
4639 if(ID &&
id!=ID)
continue;
4641 vint = &(this->wc_List[
n].cList[
id-1]);
4642 vtof = &(this->wc_List[
n].nTofF[
id-1]);
4649 for(j=0; j<V; j++) {
4651 pix = this->wc_List[
n].getPixel(
id,j);
4654 cout<<
"network::likelihood() error: NULL pointer"<<endl;
4658 if(R !=
int(pix->
rate+0.5))
continue;
4659 if(!pix->
core && core)
continue;
4668 if(NV[0].
size() < V) {
4670 for(i=0; i<
NIFO; i++) {
4680 for(j=0; j<V; j++) {
4682 pix = this->wc_List[
n].getPixel(
id,pI[j]);
4685 for(i=0; i<
nIFO; i++) {
4686 ee[
i] = 1./pix->
data[
i].noiserms;
4690 for(i=0; i<
nIFO; i++) {
4691 nv[
i][
j] = ee[
i]*ee[
i];
4692 nr[
i][
j] = ee[
i]/sqrt(cc);
4693 qdata[
i] = pdata[
i] + pix->
data[
i].index;
4698 for(i=0; i<
nIFO; i++) {
4714 STAT = -1.e64; m=IIm=0; Lm=Em=LPm= 0.;
4721 l = ind<0 ? 0 :
ind;
4723 if(!mm[l])
continue;
4726 pe[0] = eS[0].data + m0[l]*V; ,
4727 pe[1] = eS[1].
data +
m1[
l]*V; ,
4728 pe[2] = eS[2].
data +
m2[
l]*V; ,
4729 pe[3] = eS[3].
data + m3[
l]*V; ,
4730 pe[4] = eS[4].
data + m4[
l]*V; ,
4731 pe[5] = eS[5].
data + m5[
l]*V; ,
4732 pe[6] = eS[6].
data + m6[
l]*V; ,
4733 pe[7] = eS[7].
data + m7[
l]*V; )
4735 NETX(c0=0.;,c1=0.;,c2=0.;,c3=0.;,c4=0.;,c5=0.;,c6=0.;,c7=0.;)
4737 for(j=0; j<V; j++) {
4762 E = 0.; NETX(E+=c0;,E+=
c1;,E+=
c2;,E+=
c3;,E+=c4;,E+=c5;,E+=c6;,E+=c7;) E-=i*nIFO;
4763 if(E>STAT) { m =
l; STAT = E; Lm = E; }
4765 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;
4767 if(skip) { this->wc_List[
n].ignore(
id);
continue; }
4774 l = ind<0 ? 0 :
ind;
4777 skyProb.data[
l] = 0.;
4779 if(skyMaskCC.size()==
L) {
4784 double gpsT = cTo.
data[
k]+getifo(0)->TFmap.start();
4786 int cc_l = this->getIndex(th,ra);
4787 if (skyMaskCC.data[cc_l]<=0)
continue;
4790 if(!mm[l] && ind<0)
continue;
4793 pa[0] = aS[0].data + m0[l]*V; ,
4794 pa[1] = aS[1].
data +
m1[
l]*V; ,
4795 pa[2] = aS[2].
data +
m2[
l]*V; ,
4796 pa[3] = aS[3].
data + m3[
l]*V; ,
4797 pa[4] = aS[4].
data + m4[
l]*V; ,
4798 pa[5] = aS[5].
data + m5[
l]*V; ,
4799 pa[6] = aS[6].
data + m6[
l]*V; ,
4800 pa[7] = aS[7].
data + m7[
l]*V; )
4809 EE=Lo=Co=Ep=Lp=Cp=inet = 0.; II=0;
4811 for(j=0; j<V; j++) {
4813 Et = dotx(pa,j,pa,j);
4820 gp = dotx(Fp,Fp)+1.e-12;
4821 gx = dotx(Fx,Fx)+1.e-12;
4827 um = rotx(Fp,uc,Fx,us,u);
4828 hh = dotx(u,pa,j,e);
4830 cc = Le-dotx(e,e)/um;
4833 Lp+= Le*cc*ii/(Et-Le+1+cc);
4838 mulx(u,hh*ii/um,xi.
data+j*NIFO);
4842 cc = Ep>0 ? Co/(Ep-Lo+II+
fabs(Co)) : 0.;
4843 if(type==
'B') Lp = Lo*cc;
4845 skyProb.data[
l] = Lp/EE;
4847 if(Lp>LPm) { LPm=Lp; Em=EE; IIm=II; }
4849 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
4850 if(cc<this->
netCC || cc*Co<rho)
continue;
4858 for(j=0; j<V; j++) {
4870 px = xi.
data+j*NIFO;
4871 if(dotx(px,px)<=0)
continue;
4877 gp = dotx(Fp,Fp)+1.e-12;
4878 gx = dotx(Fx,Fx)+1.e-12;
4882 gc = sqrt(gR*gR+gI*gI);
4890 um = rotx(Fp,uc,Fx,us,u);
4891 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
4895 ii = netx(u,um,v,vm,GAMMA);
4897 if(ii<N_2 && gamma<0)
continue;
4902 uc = Xp*(gx+gg) - Xx*gI;
4903 us = Xx*(gp+gg) - Xp*gI;
4905 if(ii<N_1 && gamma!=0) {
4906 uc = Xp*(gc+gR)+Xx*gI;
4907 us = Xx*(gc-gR)+Xp*gI;
4912 um = rotx(Fp,uc,Fx,us,u);
4913 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
4918 rr[0] = LOCAL*am[0]/(am[0]*am[0]+2.)+1-LOCAL; ,
4919 rr[1] = LOCAL*am[1]/(am[1]*am[1]+2.)+1-LOCAL; ,
4920 rr[2] = LOCAL*am[2]/(am[2]*am[2]+2.)+1-LOCAL; ,
4921 rr[3] = LOCAL*am[3]/(am[3]*am[3]+2.)+1-LOCAL; ,
4922 rr[4] = LOCAL*am[4]/(am[4]*am[4]+2.)+1-LOCAL; ,
4923 rr[5] = LOCAL*am[5]/(am[5]*am[5]+2.)+1-LOCAL; ,
4924 rr[6] = LOCAL*am[6]/(am[6]*am[6]+2.)+1-LOCAL; ,
4925 rr[7] = LOCAL*am[7]/(am[7]*am[7]+2.)+1-LOCAL; )
4929 pp[0] = rr[0]*(am[0]-hh*u[0])*u[0]; ,
4930 pp[1] = rr[1]*(am[1]-hh*u[1])*u[1]; ,
4931 pp[2] = rr[2]*(am[2]-hh*u[2])*u[2]; ,
4932 pp[3] = rr[3]*(am[3]-hh*u[3])*u[3]; ,
4933 pp[4] = rr[4]*(am[4]-hh*u[4])*u[4]; ,
4934 pp[5] = rr[5]*(am[5]-hh*u[5])*u[5]; ,
4935 pp[6] = rr[6]*(am[6]-hh*u[6])*u[6]; ,
4936 pp[7] = rr[7]*(am[7]-hh*u[7])*u[7]; )
4941 qq[0] = rr[0]*((hh*u[0]-am[0])*v[0]+u[0]*u[0]*gg); ,
4942 qq[1] = rr[1]*((hh*u[1]-am[1])*v[1]+u[1]*u[1]*gg); ,
4943 qq[2] = rr[2]*((hh*u[2]-am[2])*v[2]+u[2]*u[2]*gg); ,
4944 qq[3] = rr[3]*((hh*u[3]-am[3])*v[3]+u[3]*u[3]*gg); ,
4945 qq[4] = rr[4]*((hh*u[4]-am[4])*v[4]+u[4]*u[4]*gg); ,
4946 qq[5] = rr[5]*((hh*u[5]-am[5])*v[5]+u[5]*u[5]*gg); ,
4947 qq[6] = rr[6]*((hh*u[6]-am[6])*v[6]+u[6]*u[6]*gg); ,
4948 qq[7] = rr[7]*((hh*u[7]-am[7])*v[7]+u[7]*u[7]*gg); )
4950 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
4953 em = rotx(u,co,v,si/vm,e);
4957 rm = rotx(v,co,u,-si/um,r)+1.e-24;
4961 pp[0] = rr[0]*(am[0]-hh*e[0])*e[0]; ,
4962 pp[1] = rr[1]*(am[1]-hh*e[1])*e[1]; ,
4963 pp[2] = rr[2]*(am[2]-hh*e[2])*e[2]; ,
4964 pp[3] = rr[3]*(am[3]-hh*e[3])*e[3]; ,
4965 pp[4] = rr[4]*(am[4]-hh*e[4])*e[4]; ,
4966 pp[5] = rr[5]*(am[5]-hh*e[5])*e[5]; ,
4967 pp[6] = rr[6]*(am[6]-hh*e[6])*e[6]; ,
4968 pp[7] = rr[7]*(am[7]-hh*e[7])*e[7]; )
4973 qq[0] = rr[0]*((hh*e[0]-am[0])*r[0]+e[0]*e[0]*gg); ,
4974 qq[1] = rr[1]*((hh*e[1]-am[1])*r[1]+e[1]*e[1]*gg); ,
4975 qq[2] = rr[2]*((hh*e[2]-am[2])*r[2]+e[2]*e[2]*gg); ,
4976 qq[3] = rr[3]*((hh*e[3]-am[3])*r[3]+e[3]*e[3]*gg); ,
4977 qq[4] = rr[4]*((hh*e[4]-am[4])*r[4]+e[4]*e[4]*gg); ,
4978 qq[5] = rr[5]*((hh*e[5]-am[5])*r[5]+e[5]*e[5]*gg); ,
4979 qq[6] = rr[6]*((hh*e[6]-am[6])*r[6]+e[6]*e[6]*gg); ,
4980 qq[7] = rr[7]*((hh*e[7]-am[7])*r[7]+e[7]*e[7]*gg); )
4982 co = dotx(qq,qq)/rm+dotx(pp,pp)/em+1.e-24;
4985 em = rotx(e,co,r,si/rm,e);
4988 Co+= (hh*hh-dotx(ee,ee))/em;
4994 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.;)
4995 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;)
4996 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.;)
4998 for(j=0; j<V; j++) {
5007 S+=
fabs(pp[0]-qq[0]); ,
5008 S+=
fabs(pp[1]-qq[1]); ,
5009 S+=
fabs(pp[2]-qq[2]); ,
5010 S+=
fabs(pp[3]-qq[3]); ,
5011 S+=
fabs(pp[4]-qq[4]); ,
5012 S+=
fabs(pp[5]-qq[5]); ,
5013 S+=
fabs(pp[6]-qq[6]); ,
5014 S+=
fabs(pp[7]-qq[7]); )
5017 pp[0] /= sqrt(qq[0]); ,
5018 pp[1] /= sqrt(qq[1]); ,
5019 pp[2] /= sqrt(qq[2]); ,
5020 pp[3] /= sqrt(qq[3]); ,
5021 pp[4] /= sqrt(qq[4]); ,
5022 pp[5] /= sqrt(qq[5]); ,
5023 pp[6] /= sqrt(qq[6]); ,
5024 pp[7] /= sqrt(qq[7]); )
5048 cc = Co/(EE-Lo+
fabs(Co));
5051 skyProb.data[
l] *= P;
5053 if(XX>=STAT) { m=
l; STAT=XX; Lm=Lo; }
5056 this->nAntenaPrior.set(l, gP/EE);
5057 this->nSensitivity.set(l, gP/EE);
5058 this->nAlignment.set(l, gX/EE);
5059 this->nNetIndex.set(l, inet/EE);
5060 this->nDisbalance.set(l, S);
5061 this->nLikelihood.set(l, Lo);
5062 this->nNullEnergy.set(l, (EE-Lo));
5063 this->nCorrEnergy.set(l, Co);
5064 this->nCorrelation.set(l, cc);
5065 this->nSkyStat.set(l, XX);
5066 this->nPenalty.set(l, P);
5067 this->nProbability.set(l,skyProb.data[l]);
5071 if(STAT<=0.001 && ind<0) {
5072 this->wc_List[
n].ignore(
id);
5080 l = ind<0 ?
m :
ind;
5083 pa[0] = aS[0].data + m0[l]*V; ,
5084 pa[1] = aS[1].
data +
m1[
l]*V; ,
5085 pa[2] = aS[2].
data +
m2[
l]*V; ,
5086 pa[3] = aS[3].
data + m3[
l]*V; ,
5087 pa[4] = aS[4].
data + m4[
l]*V; ,
5088 pa[5] = aS[5].
data + m5[
l]*V; ,
5089 pa[6] = aS[6].
data + m6[
l]*V; ,
5090 pa[7] = aS[7].
data + m7[
l]*V; )
5094 for(j=0; j<V; j++) {
5097 for(i=0; i<
nIFO; i++) {
5099 Fp[
i] = fp[
i][
l]*nr[
i][
j];
5100 Fx[
i] = fx[
i][
l]*nr[
i][
j];
5105 gp = dotx(Fp,Fp)+1.e-12;
5106 gx = dotx(Fx,Fx)+1.e-12;
5111 gc = sqrt(gR*gR+gI*gI);
5116 um = rotx(Fp,uc,Fx,us,u);
5117 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5120 if((hh*hh-dotx(e,e))/um<=0.) U=0;
5125 for(i=0; i<
nIFO; i++) {
5126 if(u[i]*u[i]/um > 1-GAMMA) ii++;
5127 if(u[i]*u[i]/um+v[i]*v[i]/vm > GAMMA) ii--;
5130 if(ii<N_2 && gamma<0.) U = 0;
5133 uc = Xp*(gx+gg) - Xx*gI;
5134 us = Xx*(gp+gg) - Xp*gI;
5136 if(ii<N_1 && gamma!=0) {
5137 uc = Xp*(gc+gR)+Xx*gI;
5138 us = Xx*(gc-gR)+Xp*gI;
5141 vc = (gp*uc + gI*us);
5142 vs = (gx*us + gI*uc);
5143 um = rotx(Fp,uc,Fx,us,u);
5144 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5152 for(i=0; i<
nIFO; i++) {
5162 for(i=0; i<
nIFO; i++) {
5163 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
5164 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
5165 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
5166 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
5169 cc = sqrt(si*si+co*co)+1.e-24;
5177 for(i=0; i<
nIFO; i++) {
5178 e[
i] = u[
i]*co + v[
i]*si;
5179 r[
i] = v[
i]*co - u[
i]*si;
5187 for(i=0; i<
nIFO; i++) {
5188 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
5189 pp[
i] = cc*(am[
i]-hh*e[
i])*e[i];
5190 qq[
i] = cc*((2.*hh*e[
i]-am[
i])*r[i] + e[i]*e[i]*gg);
5191 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
5194 cc = sqrt(si*si+co*co)+1.e-24;
5199 if(type==
'E') U = 0;
5201 for(i=0; i<
nIFO; i++) {
5202 e[
i] = e[
i]*co+r[
i]*si;
5211 Lo += (type==
'E') ? Et-nIFO : Lp;
5215 pix = this->wc_List[
n].getPixel(
id,pI[j]);
5216 pix->
likelihood = (type==
'E') ? Et-nIFO : Lp;
5217 pix->
theta = nLikelihood.getTheta(l);
5218 pix->
phi = nLikelihood.getPhi(l);
5219 if(!core) pix->
core = Et<Eo ? false :
true;
5221 for(i=0; i<
nIFO; i++) {
5223 pix->
setdata(e[i]*hh/sqrt(nv[i][j]),
'W',i);
5227 this->pixeLHood.data[pix->
time] = Lp;
5228 this->pixeLNull.data[pix->
time] = Et-Lp+1.;
5230 cout<<j<<
" SNR="<<Et<<
" ";
5231 for(i=0; i<
nIFO; i++) {
5232 cout<<e[
i]*e[
i]<<
":"<<am[
i]*am[
i]/Et<<
" ";
5244 vtof->push_back(
int(M/2)-
int(m0[l])); ,
5245 vtof->push_back(
int(M/2)-
int(
m1[l])); ,
5246 vtof->push_back(
int(M/2)-
int(
m2[l])); ,
5247 vtof->push_back(
int(M/2)-
int(m3[l])); ,
5248 vtof->push_back(
int(M/2)-
int(m4[l])); ,
5249 vtof->push_back(
int(M/2)-
int(m5[l])); ,
5250 vtof->push_back(
int(M/2)-
int(m6[l])); ,
5251 vtof->push_back(
int(M/2)-
int(m7[l])); )
5254 T = cTo.
data[
k]+getifo(0)->TFmap.start();
5255 if(iID<=0 && type!=
'E') getSkyArea(
id,n,T);
5257 if(
fabs((Lm-Lo)/Lo)>1.e-4)
5258 cout<<
"likelihood: incorrect likelihood : "<<Lm<<
" "<<Lo<<
" "<<Em<<endl;
5261 cout<<
"max value: "<<STAT<<
" at (theta,phi) = ("<<nLikelihood.getTheta(l)
5262 <<
","<<nLikelihood.getPhi(l)<<
") Likelihood: loop: "<<Lm
5263 <<
", final: "<<Lo<<
", eff: "<<Em<<endl;
5268 this->nSensitivity.gps =
T;
5269 this->nAlignment.gps =
T;
5270 this->nNetIndex.gps =
T;
5271 this->nDisbalance.gps =
T;
5272 this->nLikelihood.gps =
T;
5273 this->nNullEnergy.gps =
T;
5274 this->nCorrEnergy.gps =
T;
5275 this->nCorrelation.gps =
T;
5276 this->nSkyStat.gps =
T;
5277 this->nPenalty.gps =
T;
5278 this->nProbability.gps =
T;
5301 size_t nIFO = this->ifoList.size();
5302 size_t ID = abs(iID);
5303 int N_1 = nIFO>2 ?
int(nIFO)-1 : 2;
5304 int N_2 = nIFO>2 ?
int(nIFO)-2 : 1;
5305 if(nIFO>
NIFO || !this->
filter.size() || this->
filter[0].index.size()!=32) {
5306 cout<<
"network::likelihoodE(): invalid number of detectors or delay filter is not set.\n";
5314 double Eo = nIFO*Ao*Ao;
5316 double GAMMA = 1.-gamma*gamma;
5317 int LOCAL = local ? 1 : 0;
5320 double one = (type==
'g' || type==
'G') ? 0. : 1.;
5321 double En = (type!=
'i' && type!=
'I') ? 1.e9 : 1.;
5325 if(type==
'I' || type==
'S' || type==
'G') ISG =
true;
5326 if(type==
'i' || type==
's' || type==
'g') ISG =
true;
5327 if(type==
'i' || type==
's' || type==
'g' || type==
'r') isgr =
true;
5329 if(!ID && ind>=0) ind = -1;
5334 double gr,gg,gc,gp,
gx,cc,gP,gX,gR,gI,EE,
Et,
T,Em,LPm,HH;
5335 double em,
STAT,Lo,Lm,Lp,Co,E00,E90,L00,C00,AA,CO,SI;
5336 double co,si,stat,XP,XX,L90,C90,
sI,cO,aa,bb,as,ac,ap,La;
5337 double xp,
xx,
hh,um,vm,uc,us,vc,vs,
hp,
hx,HP,HX,wp,wx,WP,WX;
5338 size_t i,
j,
k,
l,
m,
n,V,U,
K,
id,j5;
5340 size_t I = this->ifoList[0]->TFmap.maxLayer()+1;
5341 size_t M = this->
filter.size()/I;
5342 size_t L = ind<0 ? this->
index.
size() : ind+1;
5346 double* pdata[
NIFO];
5347 double* qdata[
NIFO];
5349 for(n=0; n<
nIFO; n++) {
5350 pdata[
n] = getifo(n)->getTFmap()->data;
5351 qdata[
n] = getifo(n)->getTFmap()->data;
5355 std::vector<float>* F00;
5356 std::vector<short>* J00;
5357 std::vector<float>* F90;
5358 std::vector<short>* J90;
5368 for(n=0; n<
NIFO; n++) {
5369 fp[
n] = n<nIFO ? getifo(n)->fp.
data : f00.
data;
5370 fx[
n] = n<nIFO ? getifo(n)->fx.data : f00.
data;
5371 ffp[
n] = n<nIFO ? getifo(n)->ffp.data : f00.
data;
5372 ffm[
n] = n<nIFO ? getifo(n)->ffm.data : f00.
data;
5373 fpx[
n] = n<nIFO ? getifo(n)->fpx.data : f00.
data;
5375 short* mm = skyMask.data;
5387 std::vector<size_t> pI;
5428 for(n=0; n<
NIFO; n++) {
5439 std::vector<int>* vint;
5440 std::vector<int>* vtof;
5447 this->pixeLHood = getifo(0)->TFmap;
5448 this->pixeLHood = 0.;
5449 this->pixeLNull = getifo(0)->TFmap;
5450 this->pixeLNull = 0.;
5451 this->nSensitivity = 0.;
5452 this->nAlignment = 0.;
5453 this->nNetIndex = 0.;
5454 this->nDisbalance = 0.;
5455 this->nLikelihood = 0.;
5456 this->nNullEnergy = 0.;
5457 this->nCorrEnergy = 0.;
5458 this->nCorrelation = 0.;
5459 this->nSkyStat = 0.;
5460 this->nEllipticity = 0.;
5461 this->nPolarisation = 0.;
5462 this->nProbability = 0.;
5463 this->nAntenaPrior = 0.;
5470 for(n=lag; n<nLag; n++) {
5472 if(!this->wc_List[n].
size())
continue;
5473 logbpp = -
log(this->wc_List[n].getbpp());
5475 cid = this->wc_List[
n].
get((
char*)
"ID",0,
'S',
optim ? R : -R);
5476 cTo = this->wc_List[
n].
get((
char*)
"time",0,
'L',
optim ? R : -R);
5480 for(k=0; k<
K; k++) {
5482 id = size_t(cid.
data[k]+0.1);
5484 if(ID &&
id!=ID)
continue;
5486 vint = &(this->wc_List[
n].cList[
id-1]);
5487 vtof = &(this->wc_List[
n].nTofF[
id-1]);
5494 for(j=0; j<V; j++) {
5496 pix = this->wc_List[
n].getPixel(
id,j);
5499 cout<<
"network::likelihood() error: NULL pointer"<<endl;
5503 if(R !=
int(pix->
rate+0.5))
continue;
5504 if(!pix->
core && core)
continue;
5513 if(NV[0].
size() < V) {
5515 for(i=0; i<
NIFO; i++) {
5529 for(j=0; j<V; j++) {
5531 pix = this->wc_List[
n].getPixel(
id,pI[j]);
5534 for(i=0; i<
nIFO; i++) {
5535 ee[
i] = 1./pix->
data[
i].noiserms;
5539 GG.
data[
j] = sqrt(cc);
5540 for(i=0; i<
nIFO; i++) {
5541 nv[
i][
j] = ee[
i]*ee[
i];
5544 qdata[
i] = pdata[
i] + pix->
data[
i].index;
5549 for(i=0; i<
nIFO; i++) {
5559 gg = dot32(F00,pq,J00);
5563 gg = dot32(F90,pq,J90);
5570 STAT = -1.e64; m=U=IIm=0; Em=Lm=LPm=AA=SI=CO=0.;
5576 l = ind<0 ? 0 :
ind;
5579 skyProb.data[
l] = 0.;
5581 if(skyMaskCC.size()==
L) {
5586 double gpsT = cTo.
data[
k]+getifo(0)->TFmap.start();
5588 int cc_l = this->getIndex(th,ra);
5589 if (!skyMaskCC.data[cc_l])
continue;
5592 if(!mm[l] && ind<0)
continue;
5595 pa00[0] = a00[0].data + m0[l]*V; ,
5596 pa00[1] = a00[1].
data +
m1[
l]*V; ,
5597 pa00[2] = a00[2].
data +
m2[
l]*V; ,
5598 pa00[3] = a00[3].
data + m3[
l]*V; ,
5599 pa00[4] = a00[4].
data + m4[
l]*V; ,
5600 pa00[5] = a00[5].
data + m5[
l]*V; ,
5601 pa00[6] = a00[6].
data + m6[
l]*V; ,
5602 pa00[7] = a00[7].
data + m7[
l]*V; )
5604 EE=E00=L00=C00=Lp=0.; II=0;
5606 for(j=0; j<V; j++) {
5620 gp = dotx(Fp,Fp)+1.e-12;
5621 gx = dotx(Fx,Fx)+1.e-12;
5627 um = rotx(Fp,xp*gx-xx*gI,
5631 Co = La - dotx(e,e)/um;
5634 Lp+= La*Co*ii/(Et-La+1+Co);
5642 cc = E00>0 ? C00/(E00-L00+
fabs(C00)) : 0.;
5643 if(cc<this->
netCC || cc*C00<rho)
continue;
5645 E90=L90=C90=inet=0.;
5648 pa90[0] = a90[0].data + m0[l]*V; ,
5649 pa90[1] = a90[1].
data +
m1[
l]*V; ,
5650 pa90[2] = a90[2].
data +
m2[
l]*V; ,
5651 pa90[3] = a90[3].
data + m3[
l]*V; ,
5652 pa90[4] = a90[4].
data + m4[
l]*V; ,
5653 pa90[5] = a90[5].
data + m5[
l]*V; ,
5654 pa90[6] = a90[6].
data + m6[
l]*V; ,
5655 pa90[7] = a90[7].
data + m7[
l]*V; )
5657 for(j=0; j<V; j++) {
5668 gp = dotx(Fp,Fp)+1.e-12;
5669 gx = dotx(Fx,Fx)+1.e-12;
5675 vm = rotx(Fp,XP*gx-XX*gI,
5679 Co = La - dotx(e,e)/vm;
5682 Lp+= La*Co*ii/(Et-La+1+Co);
5684 LL.
data[
j] += La*ii;
5690 cc = E90>0 ? C90/(E90-L90+
fabs(C90)) : 0.;
5691 if(cc<this->
netCC || cc*C90<rho)
continue;
5706 for(j=0; j<V; j++) {
5711 gp = dotx(Fp,Fp)+1.e-12;
5712 gx = dotx(Fx,Fx)+1.e-12;
5716 wp = dotx(Fp,xi00.
data+j5);
5717 wx = dotx(Fx,xi00.
data+j5);
5718 WP = dotx(Fp,xi90.
data+j5);
5719 WX = dotx(Fx,xi90.
data+j5);
5721 xp = (wp*wp + WP*WP);
5722 xx = (wx*wx + WX*WX);
5725 sI += 2*(wp*wx+WP*WX-La*gI)/gg;
5726 cO += (xp-xx - La*(gp-
gx))/gg;
5727 bb += La - (xp+
xx)/gg;
5728 aa += 2*(wp*WX - wx*WP)/gg;
5731 gg = sqrt(cO*cO+sI*sI);
5732 cO = cO/gg; sI = sI/gg;
5734 if(aa> one) aa = 1.;
5735 if(aa<-one) aa =-1.;
5737 for(j=0; j<V; j++) {
5741 um = rotx(Fp,1+cO,Fx, sI,u);
5742 vm = rotx(Fx,1+cO,Fp,-sI,v);
5743 gp = dotx(u,u)+1.e-12;
5744 gx = dotx(v,v)+1.e-12;
5750 um = dotx(xi,xi)+1.e-6;
5751 hh = dotx(xi,am,pp);
5754 bb = (xp*xp+xx*
xx)/(um*gg*um);
5756 cc = bb*dotx(pp,pp);
5760 vm = dotx(XI,XI)+1.e-6;
5761 HH = dotx(XI,AM,qq);
5764 bb = (XP*XP+XX*XX)/(vm*gg*vm);
5766 cc+= bb*dotx(qq,qq);
5768 bb = 2*(xp*XX-xx*XP)/(um*gg*vm);
5770 cc+= bb*dotx(pp,qq);
5774 Et = dotx(am,am) + dotx(AM,AM);
5775 Lp+= La*cc/(Et-La+2+
fabs(cc));
5779 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
5781 if(!isgr) Lp = Lo*cc;
5782 if(Lp>LPm) {LPm=Lp; Em=EE; IIm=II;}
5784 skyProb.data[
l] = Lp/EE;
5786 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
5787 if(cc<this->
netCC || cc*Co<rho*2)
continue;
5795 for(j=0; j<V; j++) {
5799 if(dotx(xi,xi)<=0.)
continue;
5805 gp = dotx(Fp,Fp)+1.e-12;
5806 gx = dotx(Fx,Fx)+1.e-12;
5814 um = rotx(Fp,uc,Fx,us,u);
5815 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5819 ii = netx(u,um,v,vm,GAMMA);
5821 if(ii<N_2 && gamma<0.) {inix(xi,0);
continue;}
5827 if(ii<N_1 && gamma!=0) {
5829 gc = sqrt(gR*gR+gI*gI);
5830 uc = xp*(gc+gR)+xx*gI;
5831 us = xx*(gc-gR)+xp*gI;
5836 um = rotx(Fp,uc,Fx,us,u);
5837 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5842 rr[0] = LOCAL*am[0]/(am[0]*am[0]+2.)+1-LOCAL; ,
5843 rr[1] = LOCAL*am[1]/(am[1]*am[1]+2.)+1-LOCAL; ,
5844 rr[2] = LOCAL*am[2]/(am[2]*am[2]+2.)+1-LOCAL; ,
5845 rr[3] = LOCAL*am[3]/(am[3]*am[3]+2.)+1-LOCAL; ,
5846 rr[4] = LOCAL*am[4]/(am[4]*am[4]+2.)+1-LOCAL; ,
5847 rr[5] = LOCAL*am[5]/(am[5]*am[5]+2.)+1-LOCAL; ,
5848 rr[6] = LOCAL*am[6]/(am[6]*am[6]+2.)+1-LOCAL; ,
5849 rr[7] = LOCAL*am[7]/(am[7]*am[7]+2.)+1-LOCAL; )
5853 pp[0] = rr[0]*(am[0]-hh*u[0])*u[0]; ,
5854 pp[1] = rr[1]*(am[1]-hh*u[1])*u[1]; ,
5855 pp[2] = rr[2]*(am[2]-hh*u[2])*u[2]; ,
5856 pp[3] = rr[3]*(am[3]-hh*u[3])*u[3]; ,
5857 pp[4] = rr[4]*(am[4]-hh*u[4])*u[4]; ,
5858 pp[5] = rr[5]*(am[5]-hh*u[5])*u[5]; ,
5859 pp[6] = rr[6]*(am[6]-hh*u[6])*u[6]; ,
5860 pp[7] = rr[7]*(am[7]-hh*u[7])*u[7]; )
5865 qq[0] = rr[0]*((hh*u[0]-am[0])*v[0]+u[0]*u[0]*gg); ,
5866 qq[1] = rr[1]*((hh*u[1]-am[1])*v[1]+u[1]*u[1]*gg); ,
5867 qq[2] = rr[2]*((hh*u[2]-am[2])*v[2]+u[2]*u[2]*gg); ,
5868 qq[3] = rr[3]*((hh*u[3]-am[3])*v[3]+u[3]*u[3]*gg); ,
5869 qq[4] = rr[4]*((hh*u[4]-am[4])*v[4]+u[4]*u[4]*gg); ,
5870 qq[5] = rr[5]*((hh*u[5]-am[5])*v[5]+u[5]*u[5]*gg); ,
5871 qq[6] = rr[6]*((hh*u[6]-am[6])*v[6]+u[6]*u[6]*gg); ,
5872 qq[7] = rr[7]*((hh*u[7]-am[7])*v[7]+u[7]*u[7]*gg); )
5874 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
5877 em = rotx(u,co,v,si/vm,e);
5881 vm = rotx(v,co,u,-si/um,v)+1.e-24;
5885 pp[0] = rr[0]*(am[0]-hh*e[0])*e[0]; ,
5886 pp[1] = rr[1]*(am[1]-hh*e[1])*e[1]; ,
5887 pp[2] = rr[2]*(am[2]-hh*e[2])*e[2]; ,
5888 pp[3] = rr[3]*(am[3]-hh*e[3])*e[3]; ,
5889 pp[4] = rr[4]*(am[4]-hh*e[4])*e[4]; ,
5890 pp[5] = rr[5]*(am[5]-hh*e[5])*e[5]; ,
5891 pp[6] = rr[6]*(am[6]-hh*e[6])*e[6]; ,
5892 pp[7] = rr[7]*(am[7]-hh*e[7])*e[7]; )
5897 qq[0] = rr[0]*((hh*e[0]-am[0])*v[0]+e[0]*e[0]*gg); ,
5898 qq[1] = rr[1]*((hh*e[1]-am[1])*v[1]+e[1]*e[1]*gg); ,
5899 qq[2] = rr[2]*((hh*e[2]-am[2])*v[2]+e[2]*e[2]*gg); ,
5900 qq[3] = rr[3]*((hh*e[3]-am[3])*v[3]+e[3]*e[3]*gg); ,
5901 qq[4] = rr[4]*((hh*e[4]-am[4])*v[4]+e[4]*e[4]*gg); ,
5902 qq[5] = rr[5]*((hh*e[5]-am[5])*v[5]+e[5]*e[5]*gg); ,
5903 qq[6] = rr[6]*((hh*e[6]-am[6])*v[6]+e[6]*e[6]*gg); ,
5904 qq[7] = rr[7]*((hh*e[7]-am[7])*v[7]+e[7]*e[7]*gg); )
5906 co = dotx(qq,qq)/vm+dotx(pp,pp)/em+1.e-24;
5909 em = rotx(e,co,v,si/vm,e);
5912 cc = La-dotx(ee,ee)/em;
5915 Lp+= La*cc*
int(cc>0)/(Et-La+1+cc);
5923 for(j=0; j<V; j++) {
5927 if(dotx(XI,XI)<=0)
continue;
5933 gp = dotx(Fp,Fp)+1.e-12;
5934 gx = dotx(Fx,Fx)+1.e-12;
5942 um = rotx(Fp,uc,Fx,us,u);
5943 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5947 ii = netx(u,um,v,vm,GAMMA);
5949 if(ii<N_2 && gamma<0.) {inix(XI,0.);
continue;}
5955 if(ii<N_1 && gamma!=0) {
5957 gc = sqrt(gR*gR+gI*gI);
5958 uc = XP*(gc+gR)+XX*gI;
5959 us = XX*(gc-gR)+XP*gI;
5964 um = rotx(Fp,uc,Fx,us,u);
5965 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
5970 rr[0] = LOCAL*AM[0]/(AM[0]*AM[0]+2.)+1-LOCAL; ,
5971 rr[1] = LOCAL*AM[1]/(AM[1]*AM[1]+2.)+1-LOCAL; ,
5972 rr[2] = LOCAL*AM[2]/(AM[2]*AM[2]+2.)+1-LOCAL; ,
5973 rr[3] = LOCAL*AM[3]/(AM[3]*AM[3]+2.)+1-LOCAL; ,
5974 rr[4] = LOCAL*AM[4]/(AM[4]*AM[4]+2.)+1-LOCAL; ,
5975 rr[5] = LOCAL*AM[5]/(AM[5]*AM[5]+2.)+1-LOCAL; ,
5976 rr[6] = LOCAL*AM[6]/(AM[6]*AM[6]+2.)+1-LOCAL; ,
5977 rr[7] = LOCAL*AM[7]/(AM[7]*AM[7]+2.)+1-LOCAL; )
5981 pp[0] = rr[0]*(AM[0]-hh*u[0])*u[0]; ,
5982 pp[1] = rr[1]*(AM[1]-hh*u[1])*u[1]; ,
5983 pp[2] = rr[2]*(AM[2]-hh*u[2])*u[2]; ,
5984 pp[3] = rr[3]*(AM[3]-hh*u[3])*u[3]; ,
5985 pp[4] = rr[4]*(AM[4]-hh*u[4])*u[4]; ,
5986 pp[5] = rr[5]*(AM[5]-hh*u[5])*u[5]; ,
5987 pp[6] = rr[6]*(AM[6]-hh*u[6])*u[6]; ,
5988 pp[7] = rr[7]*(AM[7]-hh*u[7])*u[7]; )
5993 qq[0] = rr[0]*((hh*u[0]-AM[0])*v[0]+u[0]*u[0]*gg); ,
5994 qq[1] = rr[1]*((hh*u[1]-AM[1])*v[1]+u[1]*u[1]*gg); ,
5995 qq[2] = rr[2]*((hh*u[2]-AM[2])*v[2]+u[2]*u[2]*gg); ,
5996 qq[3] = rr[3]*((hh*u[3]-AM[3])*v[3]+u[3]*u[3]*gg); ,
5997 qq[4] = rr[4]*((hh*u[4]-AM[4])*v[4]+u[4]*u[4]*gg); ,
5998 qq[5] = rr[5]*((hh*u[5]-AM[5])*v[5]+u[5]*u[5]*gg); ,
5999 qq[6] = rr[6]*((hh*u[6]-AM[6])*v[6]+u[6]*u[6]*gg); ,
6000 qq[7] = rr[7]*((hh*u[7]-AM[7])*v[7]+u[7]*u[7]*gg); )
6002 co = dotx(qq,qq)/vm+dotx(pp,pp)/um+1.e-24;
6005 em = rotx(u,co,v,si/vm,e);
6009 vm = rotx(v,co,u,-si/um,v)+1.e-24;
6013 pp[0] = rr[0]*(AM[0]-hh*e[0])*e[0]; ,
6014 pp[1] = rr[1]*(AM[1]-hh*e[1])*e[1]; ,
6015 pp[2] = rr[2]*(AM[2]-hh*e[2])*e[2]; ,
6016 pp[3] = rr[3]*(AM[3]-hh*e[3])*e[3]; ,
6017 pp[4] = rr[4]*(AM[4]-hh*e[4])*e[4]; ,
6018 pp[5] = rr[5]*(AM[5]-hh*e[5])*e[5]; ,
6019 pp[6] = rr[6]*(AM[6]-hh*e[6])*e[6]; ,
6020 pp[7] = rr[7]*(AM[7]-hh*e[7])*e[7]; )
6025 qq[0] = rr[0]*((hh*e[0]-AM[0])*v[0]+e[0]*e[0]*gg); ,
6026 qq[1] = rr[1]*((hh*e[1]-AM[1])*v[1]+e[1]*e[1]*gg); ,
6027 qq[2] = rr[2]*((hh*e[2]-AM[2])*v[2]+e[2]*e[2]*gg); ,
6028 qq[3] = rr[3]*((hh*e[3]-AM[3])*v[3]+e[3]*e[3]*gg); ,
6029 qq[4] = rr[4]*((hh*e[4]-AM[4])*v[4]+e[4]*e[4]*gg); ,
6030 qq[5] = rr[5]*((hh*e[5]-AM[5])*v[5]+e[5]*e[5]*gg); ,
6031 qq[6] = rr[6]*((hh*e[6]-AM[6])*v[6]+e[6]*e[6]*gg); ,
6032 qq[7] = rr[7]*((hh*e[7]-AM[7])*v[7]+e[7]*e[7]*gg); )
6034 co = dotx(qq,qq)/vm+dotx(pp,pp)/em+1.e-24;
6037 em = rotx(e,co,v,si/vm,e);
6040 cc = La-dotx(ee,ee)/em;
6043 Lp+= La*cc*
int(cc>0)/(Et-La+1+cc);
6055 for(j=0; j<V; j++) {
6059 um = rotx(Fp,1+cO,Fx, sI,u);
6060 vm = rotx(Fx,1+cO,Fp,-sI,v);
6061 gp = dotx(u,u)+1.e-12;
6062 gx = dotx(v,v)+1.e-12;
6067 um = dotx(xi,xi)+1.e-6;
6068 hh = dotx(xi,am,pp);
6071 bb = (xp*xp+xx*
xx)/(um*gg*um);
6073 cc = bb*dotx(pp,pp);
6077 vm = dotx(XI,XI)+1.e-6;
6078 HH = dotx(XI,AM,qq);
6081 bb = (XP*XP+XX*XX)/(vm*gg*vm);
6083 cc+= bb*dotx(qq,qq);
6085 bb = 2*(xp*XX-xx*XP)/(um*gg*vm);
6087 cc+= bb*dotx(pp,qq);
6091 Et = dotx(am,am)+dotx(AM,AM);
6092 Lp+= La*cc*
int(cc>0)/(Et-La+2+cc);
6098 cc = EE>0 ? Co/(EE-Lo+II+
fabs(Co)) : 0.;
6099 stat = isgr ? Lp/EE : Lo*cc/EE;
6103 if(stat>=STAT) { m=
l; STAT=stat; AA=aa; CO=cO; SI=
sI; Lm=Lo/2.; }
6106 for(j=0; j<V; j++) {
6109 Et = dotx(pa00,j,pa00,j);
6110 gp = dotx(Fp,Fp)+1.e-12;
6111 gx = dotx(Fx,Fx)+1.e-12;
6115 gc = sqrt(gR*gR+gI*gI);
6120 this->nAntenaPrior.set(l, gP/V);
6121 this->nSensitivity.set(l, gP/V);
6122 this->nAlignment.set(l, gX/V);
6123 this->nLikelihood.set(l, Lo/2);
6124 this->nNullEnergy.set(l, (EE-Lo)/2.);
6125 this->nCorrEnergy.set(l, Co/2);
6126 this->nCorrelation.set(l, cc);
6127 this->nSkyStat.set(l, stat);
6128 this->nProbability.set(l, skyProb.data[l]);
6129 this->nDisbalance.set(l,
fabs(L00-L90)/
fabs(L00+L90));
6130 this->nEllipticity.set(l, aa);
6131 this->nPolarisation.set(l, atan2(sI,cO));
6132 this->nNetIndex.set(l, inet);
6136 if(STAT<0.001 && ind<0) {
6137 this->wc_List[
n].ignore(
id);
6145 l = ind<0 ?
m :
ind;
6148 pa00[0] = a00[0].data + m0[l]*V; ,
6149 pa00[1] = a00[1].
data +
m1[
l]*V; ,
6150 pa00[2] = a00[2].
data +
m2[
l]*V; ,
6151 pa00[3] = a00[3].
data + m3[
l]*V; ,
6152 pa00[4] = a00[4].
data + m4[
l]*V; ,
6153 pa00[5] = a00[5].
data + m5[
l]*V; ,
6154 pa00[6] = a00[6].
data + m6[
l]*V; ,
6155 pa00[7] = a00[7].
data + m7[
l]*V; )
6157 pa90[0] = a90[0].data + m0[l]*V; ,
6158 pa90[1] = a90[1].
data +
m1[
l]*V; ,
6159 pa90[2] = a90[2].
data +
m2[
l]*V; ,
6160 pa90[3] = a90[3].
data + m3[
l]*V; ,
6161 pa90[4] = a90[4].
data + m4[
l]*V; ,
6162 pa90[5] = a90[5].
data + m5[
l]*V; ,
6163 pa90[6] = a90[6].
data + m6[
l]*V; ,
6164 pa90[7] = a90[7].
data + m7[
l]*V; )
6168 for(j=0; j<V; j++) {
6170 pix = this->wc_List[
n].getPixel(
id,pI[j]);
6171 pix->
theta = nLikelihood.getTheta(l);
6172 pix->
phi = nLikelihood.getPhi(l);
6176 E00=E90=gp=gx=gI=xp=xx=XP=XX = 0.;
6184 for(i=0; i<
nIFO; i++) {
6185 Fp[
i] = fp[
i][
l]*nr[
i][
j];
6186 Fx[
i] = fx[
i][
l]*nr[
i][
j];
6192 gp += Fp[
i]*Fp[
i]+1.e-12;
6193 gx += Fx[
i]*Fx[
i]+1.e-12;
6209 gc = sqrt(gR*gR+gI*gI);
6216 uc = (xp*gx - xx*gI);
6217 us = (xx*gp - xp*gI);
6220 um = rotx(Fp,uc,Fx,us,u);
6221 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6225 if((hh*hh-cc)/um<=0. || E00<Eo) U=0;
6229 ii = netx(u,um,v,vm,GAMMA);
6230 if(ii<N_2 && gamma<0.) U=0;
6235 if(ii<N_1 && gamma!=0) {
6236 uc = xp*(gc+gR)+xx*gI;
6237 us = xx*(gc-gR)+xp*gI;
6242 um = rotx(Fp,uc,Fx,us,u);
6243 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6251 for(i=0; i<
nIFO; i++) {
6261 for(i=0; i<
nIFO; i++) {
6262 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
6263 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
6264 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
6265 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6268 cc = sqrt(si*si+co*co)+1.e-24;
6276 for(i=0; i<
nIFO; i++) {
6277 e[
i] = u[
i]*co + v[
i]*si;
6278 v[
i] = v[
i]*co - u[
i]*si;
6287 for(i=0; i<
nIFO; i++) {
6288 cc = local ? am[
i]/(am[
i]*am[
i]+2.) : 1.;
6289 pp[
i] = cc*(am[
i]-hh*u[
i])*u[i];
6290 qq[
i] = cc*((2.*hh*u[
i]-am[
i])*v[i] + u[i]*u[i]*gg);
6291 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6294 cc = sqrt(si*si+co*co)+1.e-24;
6300 for(i=0; i<
nIFO; i++) {
6301 u[
i] = u[
i]*co + v[
i]*si;
6308 for(i=0; i<
nIFO; i++) {
6309 pix->
setdata(u[i]*hh/sqrt(nv[i][j]),
'W',i);
6310 wp += Fp[
i]*u[
i]*
hh;
6311 wx += Fx[
i]*u[
i]*
hh;
6324 um = rotx(Fp,uc,Fx,us,u);
6325 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6329 if((hh*hh-cc)/um<=0. || E90<Eo) U=0;
6333 ii = netx(u,um,v,vm,GAMMA);
6334 if(ii<N_2 && gamma<0.) U=0;
6339 if(ii<N_1 && gamma!=0) {
6340 uc = XP*(gc+gR)+XX*gI;
6341 us = XX*(gc-gR)+XP*gI;
6346 um = rotx(Fp,uc,Fx,us,u);
6347 vm = rotx(Fx,-vc,Fp,vs,v)+1.e-24;
6355 for(i=0; i<
nIFO; i++) {
6365 for(i=0; i<
nIFO; i++) {
6366 cc = local ? AM[
i]/(AM[
i]*AM[
i]+2.) : 1.;
6367 pp[
i] = cc*(AM[
i]-hh*u[
i])*u[i];
6368 qq[
i] = cc*((2.*hh*u[
i]-AM[
i])*v[i] + u[i]*u[i]*gg);
6369 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6372 cc = sqrt(si*si+co*co)+1.e-24;
6380 for(i=0; i<
nIFO; i++) {
6381 e[
i] = u[
i]*co + v[
i]*si;
6382 v[
i] = v[
i]*co - u[
i]*si;
6391 for(i=0; i<
nIFO; i++) {
6392 cc = local ? AM[
i]/(AM[
i]*AM[
i]+2.) : 1.;
6393 pp[
i] = cc*(AM[
i]-hh*u[
i])*u[i];
6394 qq[
i] = cc*((2.*hh*u[
i]-AM[
i])*v[i] + u[i]*u[i]*gg);
6395 co += pp[
i]*pp[
i] + qq[
i]*qq[
i];
6398 cc = sqrt(si*si+co*co)+1.e-24;
6405 for(i=0; i<
nIFO; i++) {
6406 u[
i] = u[
i]*co + v[
i]*si;
6413 for(i=0; i<
nIFO; i++) {
6414 pix->
setdata(u[i]*hh/sqrt(nv[i][j]),
'U',i);
6415 WP += Fp[
i]*u[
i]*
hh;
6416 WX += Fx[
i]*u[
i]*
hh;
6424 gg = (gp+
gx)*ap + (gp-gx)*ac + 2*gI*as;
6425 hp = (wp*(ap+ac) + 2*AA*WX + wx*as)/gg;
6426 hx = (wx*(ap-ac) - 2*AA*WP + wp*as)/gg;
6427 HP = (WP*(ap+ac) - 2*AA*wx + WX*as)/gg;
6428 HX = (WX*(ap-ac) + 2*AA*wp + WP*as)/gg;
6429 La = rotx(Fp,hp,Fx,hx,e)+rotx(Fp,HP,Fx,HX,e);
6431 for(i=0; i<
nIFO; i++) {
6432 pix->
setdata((hp*fp[i][l]+hx*fx[i][l])/GG.
data[j],
'W',i);
6433 pix->
setdata((HP*fp[i][l]+HX*fx[i][l])/GG.
data[j],
'U',i);
6440 if(!core) pix->
core = (E00<Eo && E90<Eo) ?
false :
true;
6442 this->pixeLHood.data[pix->
time] = La/2.;
6443 this->pixeLNull.data[pix->
time] = (E00+E90-La)/2+1.;
6453 vtof->push_back(
int(M/2)-
int(m0[l])); ,
6454 vtof->push_back(
int(M/2)-
int(
m1[l])); ,
6455 vtof->push_back(
int(M/2)-
int(
m2[l])); ,
6456 vtof->push_back(
int(M/2)-
int(m3[l])); ,
6457 vtof->push_back(
int(M/2)-
int(m4[l])); ,
6458 vtof->push_back(
int(M/2)-
int(m5[l])); ,
6459 vtof->push_back(
int(M/2)-
int(m6[l])); ,
6460 vtof->push_back(
int(M/2)-
int(m7[l])); )
6465 T = cTo.
data[
k]+getifo(0)->TFmap.start();
6466 if(iID<=0 && type!=
'E') getSkyArea(
id,n,T);
6468 if(
fabs(Lm-Lo)/Lo>1.e-4)
6469 cout<<
"likelihood warning: "<<Lm<<
" "<<Lo<<endl;
6472 cout<<
"max value: "<<STAT<<
" at (theta,phi) = ("<<nLikelihood.getTheta(l)
6473 <<
","<<nLikelihood.getPhi(l)<<
") Likelihood: loop: " 6474 <<Lm<<
", final: "<<Lo<<
", sky: "<<LPm/2<<
", energy: "<<Em/2<<endl;
6478 if (mdcListSize() && n==0) {
6479 if (this->getwave(
id, 0,
'W')) {
6482 size_t M = this->ifoList.size();
6483 for(
int i=0; i<(
int)M; i++) {
6484 pd = this->getifo(i);
6485 pd->
RWFID.push_back(
id);
6489 double gps = this->getifo(0)->getTFmap()->start();
6497 pd->
RWFP.push_back(wf);
6503 this->nSensitivity.gps =
T;
6504 this->nAlignment.gps =
T;
6505 this->nDisbalance.gps =
T;
6506 this->nLikelihood.gps =
T;
6507 this->nNullEnergy.gps =
T;
6508 this->nCorrEnergy.gps =
T;
6509 this->nCorrelation.gps =
T;
6510 this->nSkyStat.gps =
T;
6511 this->nEllipticity.gps =
T;
6512 this->nPolarisation.gps =
T;
6513 this->nNetIndex.gps =
T;
6529 size_t N = this->ifoList.size();
6530 int N_1 = N>2 ?
int(N)-1 : 2;
6531 int N_2 = N>2 ?
int(N)-2 : 1;
6532 if(!N)
return false;
6537 vector<wavearray<double> > snr(N);
6538 vector<wavearray<double> > nul(N);
6543 std::vector<int>* vi;
6544 std::vector<wavecomplex>
A;
6558 double um,vm,cc,
hh,
gr,gc,gI,
Et,E,Xp,Xx,gp,
gx;
6559 double gg,co,si,uc,us,vc,vs,gR,
a,nr;
6567 double response = 0.;
6574 double GAMMA = 1.-gamma*gamma;
6575 bool status =
false;
6577 this->
gNET = this->aNET = this->eCOR = E = 0.;
6580 for(n=0; n<
N; n++) {
6581 nul[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'W',type);
6582 snr[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'S',type);
6583 this->getifo(n)->sSNR = 0.;
6584 this->getifo(n)->xSNR = 0.;
6585 this->getifo(n)->ekXk = 0.;
6586 this->getifo(n)->null = 0.;
6587 for(m=0; m<
NIFO; m++) { this->getifo(n)->ED[
m] = 0.; }
6588 for(m=0; m<
N; m++) { NDM[
n][
m] = 0.; }
6589 esnr[
n] = xsnr[
n] = 0.;
6592 if(!this->wc_List[lag].
size())
return status;
6594 cid = this->wc_List[lag].
get((
char*)
"ID",0,
'S',type);
6595 rat = this->wc_List[lag].
get((
char*)
"rate",0,
'S',type);
6596 lik = this->wc_List[lag].
get((
char*)
"like",0,
'S',type);
6599 for(k=0; k<
K; k++) {
6601 id = size_t(cid[k]+0.1);
6602 if(
id != ID)
continue;
6604 vi = &(this->wc_List[lag].cList[ID-1]);
6612 for(j=0; j<V; j++) {
6614 pix = this->wc_List[lag].getPixel(ID,j);
6616 cout<<
"network::setndm() error: NULL pointer"<<endl;
6619 if(!pix->
core && core)
continue;
6620 if(pix->
rate != rat.
data[k])
continue;
6623 gr=Et=gg=Xp=Xx = 0.;
6626 for(n=0; n<
N; n++) {
6632 for(n=0; n<
N; n++) {
6636 pd = this->getifo(n);
6640 A[
n].set(Fp[n],Fx[n]);
6642 gr += A[
n].abs()/2.;
6654 this->
gNET += (gr+gc)*Et;
6655 this->aNET += (gr-gc)*Et;
6664 um = vm = hh = cc = 0.;
6665 for(n=0; n<
N; n++) {
6666 u[
n] = Fp[
n]*uc + Fx[
n]*us;
6667 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
6671 cc += u[
n]*am[
n]*u[
n]*am[
n];
6675 if((hh*hh-cc)/um <= 0.)
continue;
6680 for(n=0; n<
N; n++) {
6681 if(u[n]*u[n]/um > 1-GAMMA) ii++;
6682 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
6684 this->iNET += ii*
Et;
6686 if(ii<N_2 && gamma<0.)
continue;
6689 uc = Xp*(gx+gg) - Xx*gI;
6690 us = Xx*(gp+gg) - Xp*gI;
6692 if(ii<N_1 && gamma!=0) {
6693 uc = Xp*(gc+gR)+Xx*gI;
6694 us = Xx*(gc-gR)+Xp*gI;
6700 for(n=0; n<
N; n++) {
6701 u[
n] = Fp[
n]*uc + Fx[
n]*us;
6702 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
6711 for(n=0; n<
N; n++) {
6721 for(n=0; n<
N; n++) {
6722 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
6723 qq[
n] = cc*((2*hh*u[
n]-am[
n])*v[n]+u[n]*u[n]*gg);
6724 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
6725 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
6728 cc = atan2(si,co+1.e-24);
6734 for(n=0; n<
N; n++) {
6735 e[
n] = u[
n]*co+v[
n]*si;
6736 r[
n] = v[
n]*co-u[
n]*si;
6744 for(n=0; n<
N; n++) {
6745 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
6746 pp[
n] = cc*(am[
n]-hh*e[
n])*e[n];
6747 qq[
n] = cc*((2*hh*e[
n]-am[
n])*r[n]+e[n]*e[n]*gg);
6748 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
6751 cc = atan2(si,co+1.e-24);
6756 for(n=0; n<
N; n++) {
6757 e[
n] = e[
n]*co+r[
n]*si;
6762 for(n=0; n<
N; n++) {
6764 this->getifo(n)->null += e[
n]*e[
n];
6766 for(m=0; m<
N; m++) {
6767 gg = e[
n]*am[
n]*e[
m]*am[
m];
6769 response += e[
n]*e[
m]*am[
m];
6770 if(n!=m) this->eCOR += gg;
6774 esnr[
n] += response*response;
6775 xsnr[
n] += am[
n]*response;
6776 e_SNR += response*response;
6777 x_SNR += am[
n]*response;
6778 this->getifo(n)->ED[3] +=
fabs(response*(am[n]-response));
6779 this->getifo(n)->ED[4] +=
fabs(response*(am[n]-response));
6785 this->norm = x_SNR/e_SNR;
6786 if(
fabs(this->norm-1) > 1.e-4)
6787 cout<<
"network::setndm(): incorrect likelihood normalization: "<<this->norm<<endl;
6789 for(n=0; n<
N; n++) {
6790 bIAS += this->getifo(n)->null;
6791 this->getifo(n)->null += nul[
n].data[
k];
6792 this->getifo(n)->sSNR = esnr[
n];
6793 this->getifo(n)->xSNR = xsnr[
n];
6794 S_NUL += this->getifo(n)->null;
6795 S_NIL += nul[
n].data[
k];
6796 S_SNR += snr[
n].data[
k];
6798 this->getifo(n)->ED[0] = (esnr[
n]-xsnr[
n]);
6799 this->getifo(n)->ED[1] =
fabs(esnr[n]-xsnr[n]);
6800 this->getifo(n)->ED[2] =
fabs(esnr[n]-xsnr[n]);
6802 for(m=0; m<
N; m++) S_NDM += NDM[n][m];
6806 if(count) { this->
gNET /= E; this->aNET /= E; this->iNET /= E;}
6808 a = S_NDM - lik.
data[
k];
6809 if(
fabs(a)/S_SNR>1.e-6)
6810 cout<<
"ndm-likelihood mismatch: "<<a<<
" "<<S_NDM<<
" "<<lik.
data[
k]<<
" "<<norm<<endl;
6811 a =
fabs(1 - (S_NDM+S_NIL)/S_SNR)/count;
6813 cout<<
"biased energy disbalance: "<<a<<
" "<<S_SNR-S_NDM<<
" "<<S_NIL<<
" size="<<count<<endl;
6828 size_t N = this->ifoList.size();
6829 int N_1 = N>2 ?
int(N)-1 : 2;
6830 int N_2 = N>2 ?
int(N)-2 : 1;
6831 if(!N)
return false;
6836 vector<wavearray<double> > snr(N);
6837 vector<wavearray<double> > nul(N);
6838 vector<wavearray<double> >
SNR(N);
6839 vector<wavearray<double> > NUL(N);
6844 std::vector<int>* vint;
6845 std::vector<wavecomplex>
A;
6865 double a,b,aa,
psi,gg,
gr,gc,gI,gR,E90,E00,E,fp,fx,vc,vs;
6866 double xx,xp,XX,XP,uc,us,co,
hh,gp,
gx,xi00,xi90,o00,o90;
6867 double hgw,HGW,wp,WP,wx,WX,HH,um,vm,si,cc;
6878 double GAMMA = 1.-gamma*gamma;
6881 double nC = this->MRA ? 1. : 2.;
6883 bool status =
false;
6886 if(tYPe==
'I' || tYPe==
'S' || tYPe==
'G') ISG =
true;
6887 if(tYPe==
'i' || tYPe==
's' || tYPe==
'g') ISG =
true;
6897 for(n=0; n<
N; n++) {
6898 nul[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'W',type);
6899 NUL[
n] = this->wc_List[lag].get((
char*)
"null",n+1,
'U',type);
6900 snr[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'S',type);
6901 SNR[
n] = this->wc_List[lag].get((
char*)
"energy",n+1,
'P',type);
6902 this->getifo(n)->sSNR = 0.;
6903 this->getifo(n)->xSNR = 0.;
6904 this->getifo(n)->ekXk = 0.;
6905 this->getifo(n)->null = 0.;
6906 for(m=0; m<5; m++) { this->getifo(n)->ED[
m] = 0.; }
6907 for(m=0; m<
N; m++) { NDM[
n][
m] = 0.; }
6908 esnr[
n]=xsnr[
n]=ssnr[
n]=SSNR[
n]=0.;
6911 if(!this->wc_List[lag].
size())
return false;
6913 cid = this->wc_List[lag].
get((
char*)
"ID",0,
'S',type);
6914 rat = this->wc_List[lag].
get((
char*)
"rate",0,
'S',type);
6915 lik = this->wc_List[lag].
get((
char*)
"like",0,
'S',type);
6918 for(k=0; k<
K; k++) {
6920 if(
size_t(cid[k]+0.1) != ID)
continue;
6922 vint = &(this->wc_List[lag].cList[ID-1]);
6930 for(j=0; j<V; j++) {
6932 pix = this->wc_List[lag].getPixel(ID,j);
6934 cout<<
"network::SETNDM() error: NULL pointer"<<endl;
6937 if(!pix->
core && core)
continue;
6938 if((pix->
rate != rat.
data[k]) && type)
continue;
6942 Z.
set(cos(psi),-sin(psi));
6945 gr=gg=xp=xx=XP=XX=E00=E90 = 0.;
6949 for(n=0; n<
N; n++) {
6955 for(n=0; n<
N; n++) {
6961 pd = this->getifo(n);
6964 A[
n].set(fp/b,fx/b);
6966 Fp[
n] = A[
n].real();
6967 Fx[
n] = A[
n].imag();
6968 gr += A[
n].abs()/2.;
6985 this->
gNET += (gr+gc)*(E00+E90);
6986 this->aNET += (gr-gc)*(E00+E90);
6999 um = vm = hh = cc = 0.;
7000 for(n=0; n<
N; n++) {
7001 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7002 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7006 cc += u[
n]*am[
n]*u[
n]*am[
n];
7010 if((hh*hh-cc)/um<=0. || E00<Eo) o00=0.;
7015 for(n=0; n<
N; n++) {
7016 if(u[n]*u[n]/um > 1-GAMMA) ii++;
7017 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
7019 this->iNET += ii*E00;
7020 if(ii<N_2 && gamma<0.) o00=0.;
7023 uc = xp*(gx+gg) - xx*gI;
7024 us = xx*(gp+gg) - xp*gI;
7026 if(ii<N_1 && gamma!=0) {
7027 uc = xp*(gc+gR)+xx*gI;
7028 us = xx*(gc-gR)+xp*gI;
7034 for(n=0; n<
N; n++) {
7035 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7036 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7045 for(n=0; n<
N; n++) {
7055 for(n=0; n<
N; n++) {
7056 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
7057 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
7058 qq[
n] = cc*((2.*hh*u[
n]-am[
n])*v[n] + u[n]*u[n]*gg);
7059 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7062 cc = sqrt(si*si+co*co)+1.e-24;
7070 for(n=0; n<
N; n++) {
7071 e[
n] = u[
n]*co + v[
n]*si;
7072 v[
n] = v[
n]*co - u[
n]*si;
7081 for(n=0; n<
N; n++) {
7082 cc = local ? am[
n]/(am[
n]*am[
n]+2.) : 1.;
7083 pp[
n] = cc*(am[
n]-hh*u[
n])*u[n];
7084 qq[
n] = cc*((2.*hh*u[
n]-am[
n])*v[n] + u[n]*u[n]*gg);
7085 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7088 cc = sqrt(si*si+co*co)+1.e-24;
7094 for(n=0; n<
N; n++) {
7095 u[
n] = u[
n]*co+v[
n]*si;
7098 wx += Fx[
n]*u[
n]*aa;
7100 for(n=0; n<
N; n++) {
7101 h00[
n] = u[
n]*am[
n]*o00;
7103 am[
n] = hh*u[
n]*o00;
7117 um = vm = hh = cc = 0.;
7118 for(n=0; n<
N; n++) {
7119 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7120 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7124 cc += u[
n]*AM[
n]*u[
n]*AM[
n];
7128 if((hh*hh-cc)/um<=0. || E90<Eo) o90=0.;
7133 for(n=0; n<
N; n++) {
7134 if(u[n]*u[n]/um > 1-GAMMA) ii++;
7135 if(u[n]*u[n]/um+v[n]*v[n]/vm > GAMMA) ii--;
7137 this->iNET += ii*E90;
7138 if(ii<N_2 && gamma<0.) o90=0.;
7141 uc = XP*(gx+gg) - XX*gI;
7142 us = XX*(gp+gg) - XP*gI;
7144 if(ii<N_1 && gamma!=0) {
7145 uc = XP*(gc+gR)+XX*gI;
7146 us = XX*(gc-gR)+XP*gI;
7152 for(n=0; n<
N; n++) {
7153 u[
n] = Fp[
n]*uc + Fx[
n]*us;
7154 v[
n] = Fp[
n]*vs - Fx[
n]*vc;
7163 for(n=0; n<
N; n++) {
7173 for(n=0; n<
N; n++) {
7174 cc = local ? AM[
n]/(AM[
n]*AM[
n]+2.) : 1.;
7175 pp[
n] = cc*(AM[
n]-hh*u[
n])*u[n];
7176 qq[
n] = cc*((2.*hh*u[
n]-AM[
n])*v[n] + u[n]*u[n]*gg);
7177 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7180 cc = sqrt(si*si+co*co)+1.e-24;
7188 for(n=0; n<
N; n++) {
7189 e[
n] = u[
n]*co + v[
n]*si;
7190 v[
n] = v[
n]*co - u[
n]*si;
7199 for(n=0; n<
N; n++) {
7200 cc = local ? AM[
n]/(AM[
n]*AM[
n]+2.) : 1.;
7201 pp[
n] = cc*(AM[
n]-hh*u[
n])*u[n];
7202 qq[
n] = cc*((2.*hh*u[
n]-AM[
n])*v[n] + u[n]*u[n]*gg);
7203 co += pp[
n]*pp[
n] + qq[
n]*qq[
n];
7206 cc = sqrt(si*si+co*co)+1.e-24;
7212 for(n=0; n<
N; n++) {
7213 u[
n] = u[
n]*co+v[
n]*si;
7216 WX += Fx[
n]*u[
n]*aa;
7218 for(n=0; n<
N; n++) {
7219 h90[
n] = u[
n]*AM[
n]*o90;
7221 AM[
n] = HH*u[
n]*o90;
7225 cc = ISG ? (wp*WX-wx*WP)/gg : 0.0;
7226 hh = ISG ? (wp*wp+wx*wx)/gg/nC : 1./nC;
7227 HH = ISG ? (WP*WP+WX*WX)/gg/nC : 1./nC;
7235 for(n=0; n<
N; n++) {
7236 hgw += (am[
n]*Fp[
n] + aa*AM[
n]*Fx[
n])/gg;
7237 HGW += (AM[
n]*Fp[
n] - aa*am[
n]*Fx[
n])/gg;
7240 for(n=0; n<
N; n++) {
7244 for(m=0; m<
N; m++) {
7246 a = h00[
n]*h00[
m]*hh
7250 xi00 += u00[
n]*h00[
m];
7251 xi90 += u90[
n]*h90[
m];
7254 if(n!=m) this->eCOR +=
a;
7256 a = h00[
n]*h00[
m] - h90[
n]*h90[
m];
7257 this->getifo(n)->ED[3] +=
a;
7261 a = u00[
n]*u00[
n]*hh + u90[
n]*u90[
n]*HH;
7262 this->getifo(n)->null +=
a;
7265 if(ISG) xi00 = hgw*Fp[
n]-aa*HGW*Fx[
n];
7266 esnr[
n] += xi00*xi00;
7269 this->getifo(n)->ED[1] += xi00*(xx-xi00);
7270 this->getifo(n)->ED[4] +=
fabs(xi00*(xx-xi00));
7273 if(ISG) xi90 = HGW*Fp[
n]+aa*hgw*Fx[
n];
7274 esnr[
n] += xi90*xi90;
7277 this->getifo(n)->ED[2] += xi90*(XX-xi90);
7278 this->getifo(n)->ED[4] +=
fabs(xi90*(XX-xi90));
7285 for(n=0; n<
N; n++) {
7287 b = (SSNR[
n]+ssnr[
n] - esnr[
n])/nC;
7288 this->getifo(n)->null += b;
7289 this->getifo(n)->sSNR = esnr[
n]/nC;
7290 this->getifo(n)->xSNR = xsnr[
n]/nC;
7295 for(m=0; m<
N; m++) S_NDM += NDM[n][m];
7297 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)
7298 cout<<
"SETNDM()-likelihoodI() SNR mismatch: "<<n<<
" "<<ID<<
" " 7299 <<ssnr[
n]<<
":"<<snr[
n].data[
k]<<
" "<<SSNR[
n]<<
":"<<SNR[
n].data[
k]<<endl;
7302 if(count) { this->
gNET /= E; this->aNET /= E; this->iNET /= E; }
7304 a = S_NDM - lik.
data[
k];
7305 if(
fabs(a)/S_NDM>1.e-6)
7306 cout<<
"NDM-likelihood mismatch: "<<a<<
" "<<S_NDM<<
" "<<lik.
data[
k]<<endl;
7308 a =
fabs(1 - nC*(S_NDM+S_NUL)/(s_snr+S_SNR))/count;
7310 cout<<
"biased energy disbalance: "<<a<<
" "<<S_NDM+S_NUL<<
" "<<(s_snr+S_SNR)/nC<<endl;
7323 const char*
fname,
const char* fmode,
size_t*
lagSite) {
7325 size_t nIFO = this->ifoList.size();
7326 this->wc_List.clear(); this->livTime.clear();
7329 cout <<
"network::setTimeShifts : lagStep must be positive" << endl;
7333 if(lagSize<1) lagSize=1;
7335 if(strcmp(fmode,
"r") && strcmp(fmode,
"w") && strcmp(fmode,
"s")) {
7336 cout <<
"network::setTimeShifts : bad fmode : must be r/w/s" << endl;
7340 if(fname) {
if(strlen(fname)<1) fname = NULL; }
7353 int maxIter = 10000000;
7356 for(n=0;n<
NIFO;n++) {
7368 lagIDS +=
int(getifo(0)->sHIFt/lagStep);
7369 maxList = lagSize+lagIDS;
7371 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7372 for(m=0; m<maxList; m++) {
7373 for(n=0; n<
nIFO; n++) {
7374 pd = this->getifo(n);
7375 lagList[
n][
m] = n==0 ?
m :
int(pd->
sHIFt/lagStep);
7384 if(fname && (!strcmp(fmode,
"r") || !strcmp(fmode,
"s"))) {
7385 if(!strcmp(fmode,
"r")) {
7390 cout <<
"network::setTimeShifts : Error Opening File : " << fname << endl;
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);
7412 for(n=0; n<nIFO; n++) in >> lagList[
n][
m];
7413 if (!in.good())
break;
7415 in.seekg(fpos+1, ios::beg);
7421 if(!strcmp(fmode,
"s")) {
7430 in.getline(str,1024);
7431 if (!in.good())
break;
7432 if(str[0] !=
'#') maxList++;
7435 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7436 in.clear(ios::goodbit);
7437 in.seekg(0, ios::beg);
7440 in.getline(str,1024);
7441 if(str[0] ==
'#')
continue;
7442 in.seekg(fpos, ios::beg);
7444 for(n=0; n<nIFO; n++) in >> lagList[
n][
m];
7445 if (!in.good())
break;
7447 in.seekg(fpos+1, ios::beg);
7455 for(m=0; m<maxList; m++){
7457 for (n=0; n<
nIFO; n++)
id[n]=lagList[n][m];
7461 for (n=0; n<nIFO; n++) if(id[n]<0||id[n]>
int(lagMax)) check=
false;
7465 for (
int i=nIFO-1;
i>=0;
i--) {
7466 for (
int j=
i-1;
j>=0;
j--) {
7467 if (!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7468 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7475 cout <<
"network::setTimeShifts : no lags in the list" << endl;
7476 cout <<
"lagP : " << lagP <<
" " << lagSize << endl;
7479 if(lagP!=
int(maxList)) {
7480 cout <<
"network::setTimeShifts : lags out of lagMax" << endl;
7481 cout <<
"lagP : " << lagP <<
" " << lagSize << endl;
7490 if(lagSite!=NULL)
for(n=0; n<
nIFO; n++) {
7491 if(lagSite[n] >= nIFO) {
7492 cout <<
"network::setTimeShifts : Error lagSite - value out of range " << endl;
7497 for(n=1; n<
nIFO; n++) N[n]=lagMax;
7501 for(n=0; n<
nIFO; n++) lagList[n] =
new int[maxList];
7502 for(n=0; n<
nIFO; n++) lagList[n][nList]=0;
7506 for (
int k=0;k<maxIter;k++) {
7507 for(n=0; n<
nIFO; n++) ID[n] = TMath::Nint(rnd.Uniform(-(N[n]+0.5),N[n]+0.5));
7508 for(n=0; n<
nIFO; n++)
id[n] = (lagSite==NULL) ? ID[
n] : ID[lagSite[
n]];
7510 for(
int i=nIFO-1;
i>=0;
i--) {
7511 for(
int j=
i-1;
j>=0;
j--) {
7512 if(!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7513 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7515 if(lagSite[
i]!=lagSite[j] &&
id[
i]==
id[j]) check=
false;
7517 if(
id[
i]==
id[j]) check=
false;
7523 for(m=0;m<nList;m++) {
7525 for(n=0; n<
nIFO; n++)
if(lagList[n][m]!=
id[n]) pass=
false;
7526 if(pass) check=
false;
7530 if(NETX(
id[0]||,
id[1]||,
id[2]||,
id[3]||,
id[4]||,
id[5]||,
id[6]||,
id[7]||)
false) {
7531 for(n=0; n<
nIFO; n++) lagList[n][nList]=
id[n];
7535 if (nList>=maxList)
break;
7543 for(m=0; m<nList; m++) {
7544 int lagMin = kMaxInt;
7545 for(n=0; n<
nIFO; n++)
if (lagList[n][m]<lagMin) lagMin=lagList[
n][
m];
7546 for(n=0; n<
nIFO; n++) lagList[n][m]-=lagMin;
7549 if(lagIDS+lagSize>nList) {
7550 cout <<
"network::setTimeShifts : lagOff+lagSize > nList of lags : " << nList << endl;
7554 for(n=0; n<
nIFO; n++){
7555 pd = this->getifo(n);
7564 double R = this->getifo(0)->getTFmap()->
rate();
7565 double segLen = this->getifo(0)->getTFmap()->size();
7566 double edge = this->Edge;
7573 segLen = (segLen/R-2*edge)/lagStep;
7574 lagMaxSeg =
int(segLen)-1;
7576 for(n=0; n<
nIFO; n++) {
7578 lagH[
n] = lagMaxSeg;
7583 for (n=0; n<
nIFO; n++)
id[n]=lagList[n][m+lagIDS];
7586 for(n=0; n<nIFO; n++) if(id[n]<0||id[n]>
int(lagMaxSeg)) check=
false;
7589 for(
int i=nIFO-1;
i>=0;
i--) {
7590 for(
int j=
i-1;
j>=0;
j--) {
7591 if (!(((
id[
i]-
id[
j])>=(lagL[
i]-lagH[j]))&&
7592 ((
id[
i]-
id[j])<=(lagH[
i]-lagL[j])))) check=
false;
7600 for(n=0; n<
nIFO; n++) {
7601 k = lagList[
n][m+lagIDS];
7602 if(k) check =
false;
7603 this->getifo(n)->lagShift.data[selSize] = k*
lagStep;
7607 k = lagList[0][m+lagIDS];
7608 this->getifo(0)->lagShift.data[selSize] = k*
lagStep;
7610 for(n=1; n<
nIFO; n++) {
7611 pd = this->getifo(n);
7615 if(zero>0.1) check =
false;
7618 wc_List.push_back(wc);
7619 livTime.push_back(0.);
7625 cout <<
"network::setTimeShifts error: no lag was selected" << endl;
7629 for(n=0; n<
nIFO; n++) {
7630 m = this->getifo(n)->lagShift.size();
7631 if(m!=selSize) this->getifo(n)->lagShift.resize(selSize);
7636 if(fname && !strcmp(fmode,
"w") && lagMax) {
7639 if((fP = fopen(fname,
"w")) == NULL) {
7640 cout <<
"network::setTimeShifts error: cannot open file " << fname << endl;
7646 fprintf(fP,
"#total %10d lags \n",
int(nList));
7648 fprintf(fP,
"#%13s%14s%14s\n",
" nIFO",
"lagStep",
" lagMax");
7649 fprintf(fP,
"#%13d%14.3f%14d\n",
int(nIFO),lagStep,
int(lagMax));
7652 for(n=0; n<
nIFO; n++)
fprintf(fP,
"%12s-%1d",
"lagShift",
int(n));
7657 for(m=0; m<nList; m++){
7659 for (n=0; n<
nIFO; n++)
fprintf(fP,
"%14d",lagList[n][m]);
7668 for(n=0; n<
nIFO; n++)
delete [] lagList[n];
7673 for(n=0; n<
nIFO; n++)
printf(
"%12.12s%2s",
"ifo",getifo(n)->Name);
7675 for(m=0; m<selSize; m++){
7677 for(n=0; n<
nIFO; n++)
printf(
"%14.5f",this->getifo(n)->lagShift.data[m]);
7693 std::vector<delayFilter>().swap(this->
filter);
7695 if(!d) d = this->getifo(0);
7698 if(ifoList.size()<2 || !d || rate==0.) {
7699 cout<<
"network::setFilter() error: incomplete network initialization"<<endl;
7703 double T = this->getDelay((
char*)
"MAX")+0.002;
7710 int J =
int((
fabs(T)*rate*N)/M+0.5);
7712 if(N<M) { cout<<
"network::setFilter() error"<<endl;
return 0; }
7715 this->getifo(0)->nDFS =
N;
7716 this->
filter.reserve((2*J-1)*M);
7718 for(i=0; i<
M; i++) {
7719 for(j=-(J-1); j<J; j++) {
7720 m = j>0 ? (j+N/2-1)/N : (j-N/2)/
N;
7725 for(k=0; k<
K; k++) v.
index[k] -= m*M;
7726 this->filter.push_back(v);
7741 size_t N = this->ifoList.size();
7743 if(d) this->getifo(0)->setFilter(*d);
7745 this->getifo(0)->clearFilter();
7756 size_t N = this->ifoList.size();
7760 this->getifo(0)->readFilter(gname);
7762 this->filter90.clear();
7763 std::vector<delayFilter>().swap(this->filter90);
7764 this->filter90 = this->
filter;
7766 this->getifo(0)->readFilter(fname);
7768 this->getifo(0)->clearFilter();
7778 size_t N = this->ifoList.size();
7781 this->readFilter(gname);
7782 this->filter90 = this->
filter;
7784 this->readFilter(fname);
7796 if ( (fp=fopen(fname,
"wb")) == NULL ) {
7797 cout <<
" network::writeFilter() error : cannot open file " << fname <<
". \n";
7801 size_t M = size_t(getifo(0)->TFmap.maxLayer()+1);
7802 size_t N = size_t(
filter.size()/
M);
7804 size_t n = K *
sizeof(float);
7805 size_t m = K *
sizeof(short);
7810 fwrite(&K,
sizeof(
size_t), 1, fp);
7811 fwrite(&M,
sizeof(
size_t), 1, fp);
7812 fwrite(&N,
sizeof(
size_t), 1, fp);
7814 for(i=0; i<
M; i++) {
7815 for(j=0; j<
N; j++) {
7816 for(k=0; k<
K; k++) {
7820 fwrite(value.data, n, 1, fp);
7821 fwrite(index.
data, m, 1, fp);
7835 if ( (fp=fopen(fname,
"rb")) == NULL ) {
7836 cout <<
" network::readFilter() error : cannot open file " << fname <<
". \n";
7844 fread(&K,
sizeof(
size_t), 1, fp);
7845 fread(&M,
sizeof(
size_t), 1, fp);
7846 fread(&N,
sizeof(
size_t), 1, fp);
7848 size_t n = K *
sizeof(float);
7849 size_t m = K *
sizeof(short);
7858 for(k=0; k<
K; k++) {
7859 v.
value.push_back(0.);
7860 v.
index.push_back(0);
7863 for(i=0; i<
M; i++) {
7864 for(j=0; j<
N; j++) {
7865 fread(value.data, n, 1, fp);
7866 fread(index.
data, m, 1, fp);
7867 for(k=0; k<
K; k++) {
7888 int nres = this->wdmList.size();
7889 int nIFO = this->ifoList.size();
7890 int V =
int(this->pList.size());
7898 for(
int i=0;
i<nres;
i++) {
7900 pd = this->getifo(
k);
7901 pd->
vSS[
i].Expand(
false);
7904 pd->
vSS[
i].Inverse();
7906 pd->
vSS[
i].getLayer(x,0);
7910 pd->
vSS[
i].Forward(x);
7911 layers = pd->
vSS[
i].maxLayer()+1;
7913 for(
int j=0;
j<V;
j++) {
7914 pix = this->pList[
j];
7915 if(pix->
layers != layers)
continue;
7917 v00[
k][
j] = pd->
vSS[
i].GetMap00(ind);
7918 v90[
k][
j] = pd->
vSS[
i].GetMap90(ind);
7921 pd->
vSS[
i].Shrink();
7934 size_t N = this->ifoList.size();
7935 size_t k = this->getIndex(theta,phi);
7938 for(
size_t n=1;
n<
N;
n++){
7939 d = this->getifo(
n);
7957 size_t M = this->
filter.size()/I;
7958 size_t K = this->
filter[0].index.size();
7959 size_t jB = size_t(this->Edge*R/I)*I;
7966 double*
F = (
double*)malloc(K*
sizeof(
double));
7967 int* J = (
int*)malloc(K*
sizeof(
int));
7977 double* b0 = temp.
data;
7979 for(i=0; i<I; i++) {
7984 F[
k] = double(pv->
value[k]);
7991 for(j=jS; j<
N; j+=I) {
8011 size_t N = ifoList.size();
8018 cout<<
"network::setDelayIndex(): invalid network\n";
8023 for(n=0; n<
N; n++) dr[n] = ifoList[n];
8025 size_t I = dr[0]->
nDFL;
8026 size_t K = this->
filter.size()/I;
8030 rate *= dr[0]->
nDFS/I;
8032 if(pOUT) cout<<
"filter size="<<this->
filter.size()
8033 <<
" layers="<<I<<
" delays="<<K<<
" samples="<<dr[0]->
nDFS<<endl;
8035 if(!(K&1) || rate == 0.) {
8036 cout<<
"network::setDelayIndex(): invalid network\n";
8040 for(n=0; n<
N; n++) {
8051 this->nPenalty = dr[0]->
tau;
8060 for(n=0; n<
N; n++) {
8061 for(m=0; m<
N; m++) {
8063 i =
int(t*rate+2*K+0.5) - 2*
K;
8069 for(n=0; n<
N; n++) {
8071 for(m=0; m<
N; m++) {
8072 for(k=0; k<
N; k++) {
8073 t =
fabs(mm[n][k]-mm[n][m]-tt[m][k]);
8074 if(TT[n] < t) TT[
n] =
t;
8080 for(m=0; m<
N; m++) {
8081 if(t>TT[m]) { t = TT[
m]; k =
m; }
8083 this->nPenalty.set(l,
double(t));
8086 i = mIFO<9 ? mm[
k][this->mIFO] :
int(t*rate+2*K+0.5)-2*
K;
8092 for(m=0; m<
N; m++) {
8093 ii = (K/2+mm[
k][
m])-i;
8095 if(ii < 0) cout<<
"network::setDelayIndex error: sky index<0: "<<k<<endl;
8109 if(ifoList.size()<2 || !dr->
tau.
size()) {
8110 cout<<
"network::setIndex() - invalid network"<<endl;
8116 size_t N = ifoList.size();
8117 size_t I = mode ? dr->
nDFL : 1;
8120 size_t K = this->
filter.size()/I;
8122 size_t M = mIFO<9 ? mIFO : 0;
8126 if(this->skyMask.size()!=
L) this->skyMask.resize(L);
8127 if(this->skyHole.size()!=
L) { this->skyHole.resize(L); this->skyHole = 1.; }
8128 for(j=0; j<
L; j++) {
8130 skyMask.data[
j] = size_t(skyHole.data[j]+0.1);
8134 if(mode==2 || mode==4) {
8140 long long **pp = (
long long**)malloc(L*
sizeof(
long long*));
8143 for(n=0; n<
N; n++) {
8144 if(!this->getifo(n)->index.size()) {
8145 cout<<
"network::setIndex() - invalid network"<<endl;
8154 for(n=0; n<
N; n++) {
8155 if(n == M)
continue;
8156 ll = this->getifo(n)->index.data[
i];
8157 if(this->mIFO==99) ll += K/2 - this->getifo(0)->index.data[
i];
8158 delay.
data[
i] += ll<<(m*12);
8165 for(i=1; i<
L; i++) {
8166 j = pp[
i] - delay.
data;
8167 if(ll == *(pp[i])) {
8168 skyMask.data[
j] = 0;
8174 if(pOUT) cout<<
"\n ll="<<ll<<
" "<<j<<
"|"<<sm->
getTheta(j)<<
"|"<<sm->
getPhi(j);
8186 int nIFO = ifoListSize();
8187 for(
int n=0;
n<
nIFO;
n++) getifo(
n)->print();
8191 cout <<
"----------------------------------------------" << endl;
8192 cout <<
" INJECTIONS : " << this->mdcListSize() << endl;
8193 cout <<
"----------------------------------------------" << endl;
8194 for(
size_t k=0;
k<this->mdcListSize();
k++) {
8195 string str(this->getmdcList(
k));
8196 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)
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...
size_t readMDClog(char *, double=0., int=11, int=12)
param: MDC log file param: approximate gps time
virtual void rate(double r)
std::vector< wavearray< float > > tdAmp
static void _avx_pol_ps(float **p, float **q, wavearray< double > *pol00, wavearray< double > *pol90, std::vector< float *> &pAPN, std::vector< float *> &pAVX, int II)
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)
static float _avx_dpf_ps(double **Fp, double **Fx, int l, std::vector< float *> &pAPN, std::vector< float *> &pAVX, int I)
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
std::vector< 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
static __m256 _avx_stat_ps(float **x, float **X, float **s, float **S, std::vector< float *> &pAVX, int I)
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
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
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
static float _avx_GW_ps(float **p, float **q, std::vector< float *> &pAPN, float *rr, std::vector< float *> &pAVX, int II)
std::vector< vector_float > p_Map
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
void writeFilter(const char *fname)
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
watplot p2(const_cast< char *>("TFMap2"))
virtual void waveSplit(DataType_t **pp, size_t l, size_t r, size_t m) const
virtual size_t size() const
size_t wavecount(double x, int n=0)
long likelihoodWP(char mode, int lag, int ID, TH2F *hist=NULL)
double setVeto(double=5.)
param: time window around injections
static float _avx_packet_ps(float **p, float **q, std::vector< float *> &pAVX, int I)
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 __m128 _sse_ecoh4_ps(__m128 *_p, __m128 *_q, __m128 _L)
wavearray< double > * getHoT()
param: no parameters
gwavearray< double > * gx
watplot p1(const_cast< char *>("TFMap1"))
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)
static void _avx_free_ps(std::vector< float *> &v)
WSeries< double > pTF[nRES]
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
printf("total live time: non-zero lags = %10.1f \, liveTot)
std::vector< vector_int > p_Ind
static void _sse_pol4_ps(__m128 *_fp, __m128 *_fx, __m128 *_v, double *r, double *a)
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")
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)
static float _avx_ort_ps(float **p, float **q, std::vector< float *> &pAVX, 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)
static float _avx_setAMP_ps(float **p, float **q, std::vector< float *> &pAVX, int I)
WSeries< double > waveForm
WSeries< double > * getTFmap()
param: no parameters
virtual void delay(double T)
virtual void stop(double s)
double fabs(const Complex &x)
virtual void waveSort(DataType_t **pp, size_t l=0, size_t r=0) const
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)
static float _avx_loadata_ps(float **p, float **q, float **u, float **v, float En, std::vector< float *> &pAVX, int I)
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\", 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)
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
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)
static __m256 _avx_noise_ps(float **p, float **q, std::vector< float *> &pAVX, int I)
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...