3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
42 cout <<
"-----> CWB_Plugin_QLveto.C" << endl;
43 cout <<
"ifo " << ifo.Data() << endl;
44 cout <<
"type " << type << endl;
59 TList *
files = (TList*)gROOT->GetListOfFiles();
68 while ((file=(TSystemFile*)
next())) {
69 fname = file->GetName();
71 if(fname.Contains(
"wave_")) {
72 froot=(TFile*)file;froot->cd();
74 outDump.ReplaceAll(
".root.tmp",
".txt");
79 cout <<
"CWB_Plugin_QLveto.C : Error - output root file not found" << endl;
83 cout <<
"CWB_Plugin_QLveto.C : Error - output root file not found" << endl;
87 TTree*
net_tree = (TTree *) froot->Get(
"waveburst");
91 net_tree->Branch(
"Qveto",Qveto,TString::Format(
"Qveto[%i]/F",2));
92 net_tree->Branch(
"Lveto",Lveto,TString::Format(
"Lveto[%i]/F",3));
95 bool qveto_exists=
false;
96 bool lveto_exists=
false;
97 TIter
next(net_tree->GetListOfBranches());
98 while ((branch=(TBranch*)
next())) {
99 if(
TString(
"Qveto").CompareTo(branch->GetName())==0) qveto_exists=
true;
100 if(
TString(
"Lveto").CompareTo(branch->GetName())==0) lveto_exists=
true;
103 if(!qveto_exists) net_tree->Branch(
"Qveto",Qveto,TString::Format(
"Qveto[%i]/F",2));
104 if(!lveto_exists) net_tree->Branch(
"Lveto",Lveto,TString::Format(
"Lveto[%i]/F",3));
111 cout <<
"CWB_Plugin_QLveto.C -> "
112 <<
"CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
118 cout <<
"-----> CWB_Plugin_QLveto.C -> "
119 <<
" gIFACTOR : " << gIFACTOR << endl;
122 float* gSLAGSHIFT=NULL;
IMPORT(
float*,gSLAGSHIFT)
134 TList *
files = (TList*)gROOT->GetListOfFiles();
141 while ((file=(TSystemFile*)
next())) {
142 fname = file->GetName();
144 if(fname.Contains(
"wave_")) {
145 froot=(TFile*)file;froot->cd();
147 outDump.ReplaceAll(
".root.tmp",
".txt");
152 cout <<
"CWB_Plugin_QLveto.C : Error - output root file not found" << endl;
156 cout <<
"CWB_Plugin_QLveto.C : Error - output root file not found" << endl;
160 TTree*
net_tree = (TTree *) froot->Get(
"waveburst");
163 net_tree->SetBranchAddress(
"Qveto",Qveto);
164 net_tree->SetBranchAddress(
"Lveto",Lveto);
168 net_tree->Branch(
"Qveto",Qveto,TString::Format(
"Qveto[%i]/F",2));
169 net_tree->Branch(
"Lveto",Lveto,TString::Format(
"Lveto[%i]/F",3));
173 for(
int k=0;
k<
K;
k++) {
179 int ID = size_t(
id.data[
j]+0.5);
192 int idSize = pd->
RWFID.size();
195 for (
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
196 if(wfIndex==-1)
continue;
199 cout << endl <<
"----------------------------------------------------------------" << endl;
202 Qveto[0]=Qveto[1]=1.e20;
207 pwfREC[
n] = pd->
RWFP[wfIndex];
210 #ifdef PLOT_WHITENED_WAVEFORMS
220 if(qveto<Qveto[0]) Qveto[0]=
qveto;
221 if(qfactor<Qveto[1]) Qveto[1]=qfactor;
226 if(qveto<Qveto[0]) Qveto[0]=
qveto;
227 if(qfactor<Qveto[1]) Qveto[1]=qfactor;
229 cout <<
"Qveto : " << pd->
Name <<
" Qveto[0] = " << Qveto[0]
230 <<
" Qveto[1] = " << Qveto[1] << endl;
236 std::vector<netpixel> vPIX;
242 cout <<
"Lveto : " <<
"fmean : " << Lveto[0] <<
" frms : " << Lveto[1]
243 <<
" Energy Ratio : " << Lveto[2] << endl << endl;
244 cout <<
"----------------------------------------------------------------" << endl;
253 if(cfg->
dump) EVT->
dopen(outDump.Data(),
const_cast<char*
>(
"a"),
false);
254 EVT->
output2G(net_tree,NET,ID,
k,ofactor);
262 for(
int i=0; i<3; i++)
fprintf(EVT->
fP,
"%f ",Lveto[i]);
286 for(
int j=xsize;
j<x.
size();
j++) x[
j]=0;
309 if(
a[
i]>amax) {amax=
a[
i];imax=
i;}
326 if(abs(imax-
i)<=
NTHR) {
334 Qveto = ein>0 ? eout/ein : 0.;
338 float R = (
a[imax-1]+
a[imax+1])/
a[imax]/2.;
358 Lveto[0] = Lveto[1] = Lveto[2] = 0;
360 std::vector<int>* vint = &(pwc->
cList[cid-1]);
361 int V = vint->size();
372 for(
int n=0;
n<V;
n++) {
375 cout <<
"CWB_Plugin_QLveto.C - Error : is enabled only for WDM (2G)" << endl;
378 if(!pix->
core)
continue;
382 likePix += pow(pix->
getdata(
'S',
m),2);
383 likePix += pow(pix->
getdata(
'P',
m),2);
390 if(likePix>likeMax) {likeMax=likePix;freqMax=freq;dfreqMax=
df;}
402 for(
int n=0;
n<V;
n++) {
404 if(!pix->
core)
continue;
408 likePix += pow(pix->
getdata(
'S',
m),2);
409 likePix += pow(pix->
getdata(
'P',
m),2);
415 if(
fabs(freq-freqMax)<=dfreqMax) {
417 fmean += likePix*freq;
418 frms += likePix*freq*freq;
422 fmean = fmean/likeLin;
423 frms = frms/likeLin-fmean*fmean;
424 frms = frms>0 ? sqrt(frms) : 0.;
426 if(frms<dfreqMax/2.) frms=dfreqMax/2.;
434 Lveto[2] = likeTot>0. ? likeLin/likeTot : 0.;
440 #if defined PLOT_LIKELIHOOD
443 WTS.
canvas->Print(
"l_tfmap_scalogram.png");
456 double tmin = wfREC->
start();
463 PTS.
graph[0]->SetLineWidth(1);
467 if(fft)
sprintf(label,
"%s",
"fft");
468 else sprintf(label,
"%s",
"time");
469 if(strain)
sprintf(label,
"%s_%s",label,
"strain");
470 else sprintf(label,
"%s_%s",label,
"white");
476 cout <<
"write : " << fname << endl;
485 n = ifo->
IWFP.size();
486 for (
int i=0;
i<
n;
i++) {
493 n = ifo->
RWFP.size();
494 for (
int i=0;
i<
n;
i++) {
515 std::vector<netpixel> vPIX;
519 if(V>cfg->
BATCH)
return vPIX;
524 __m128* _xi = (__m128*) xi.
data;
525 __m128* _XI = (__m128*) XI.
data;
527 __m128* _aa = (__m128*) NET->
a_00.
data;
528 __m128* _AA = (__m128*) NET->
a_90.
data;
531 for(
int j=0;
j<V;
j++) {
536 if(ee>En) nPC++;
else ee=0.;
543 for(
int j=0;
j<V;
j++) {
545 vPIX.push_back(*pix);
563 if(V>cfg->
BATCH)
return;
564 for(
int j=0;
j<V;
j++) {
569 while(!vPIX->empty()) vPIX->pop_back();
570 vPIX->clear(); std::vector<netpixel>().swap(*vPIX);
monster wdmMRA
list of pixel pointers for MRA
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
static float _sse_abs_ps(__m128 *_a)
virtual size_t size() const
void GetLveto(netcluster *pwc, int cid, int nifo, float *Lveto)
void setSLags(float *slag)
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
virtual void rate(double r)
wavearray< double > a(hp.size())
void ResetPCA(network *NET, CWB::config *cfg, netcluster *pwc, std::vector< netpixel > *vPIX, int ID)
static void _sse_zero_ps(__m128 *_p)
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
wavearray< float > a_90
buffer for cluster sky 00 amplitude
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
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< wavearray< double > * > RWFP
int _sse_mra_ps(float *, float *, float, int)
std::vector< TGraph * > graph
WSeries< double > waveBand
std::vector< vector_int > cList
virtual void start(double s)
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
std::vector< detector * > ifoList
std::vector< netpixel > DoPCA(network *NET, CWB::config *cfg, int lag, int id)
void dopen(const char *fname, char *mode, bool header=true)
wavearray< float > pNRG
buffers for cluster residual energy
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
void output2G(TTree *, network *, size_t, int, double)
#define IMPORT(TYPE, VAR)
xtalk getXTalk(int nLay1, size_t indx1, int nLay2, size_t indx2)
param: numbers of layers identifying the resolution of the first pixel param: TF map index of the fir...
TCut qveto("qveto","Qveto[0]>0.3 && Qveto[1]>0.3")
void ClearWaveforms(detector *ifo)
wavearray< float > a_00
wdm multi-resolution analysis
TIter next(twave->GetListOfBranches())
std::vector< wavearray< double > * > IWFP
netpixel * getPixel(size_t n, size_t i)
FILE * fP
injected reconstructed xcor waveform
void PlotWaveform(TString ifo, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
netevent EVT(itree, nifo)
WSeries< double > waveForm
double fabs(const Complex &x)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void GetQveto(wavearray< double > *wf, float &Qveto, float &Qfactor)
wavearray< float > rNRG
buffer for cluster sky 90 amplitudes
double factors[FACTORS_MAX]
double getdata(char type='R', size_t n=0)
virtual void resize(unsigned int)
bool setdata(double a, char type='R', size_t n=0)