3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
20 #define OUTPUT_DIR "plots"
21 #define OUTPUT_TYPE "png"
31 cout <<
"-----> plugins/CWB_Plugin_pixeLHood.C" << endl;
32 cout <<
"ifo " << ifo.Data() << endl;
33 cout <<
"type " << type << endl;
48 bool batch = gROOT->IsBatch();
49 gROOT->SetBatch(
true);
55 char strtime[1024] =
"";
61 double eventTime[
NIFO];
74 for(j=0; j<(
int)
id.
size(); j++) {
76 ID = size_t(
id.data[j]+0.5);
78 EVT->
output(NULL,NET,factor,ID,k);
82 int lagMin=2147483647;
83 for(n=0; n<
nIFO; n++)
if(EVT->
lag[n]<lagMin) {lagMin=
int(EVT->
lag[n]);masterDet=
n;}
88 double gps_start = EVT->
time[masterDet]-EVT->
duration[masterDet];
89 double gps_stop = EVT->
time[masterDet]+EVT->
duration[masterDet];
95 sprintf(fname,
"%s/eventDump%s.%s", outdir, strtime,
"txt");
97 EVT->
output(NULL,NET,factor,ID,k);
101 for(n=0; n<
nIFO; n++) eventTime[n]=(EVT->
start[n]+EVT->
stop[n])/2.;
102 for(n=0; n<
nIFO; n++) minTime = (eventTime[n]<minTime) ? eventTime[
n] : minTime;
103 for(n=0; n<
nIFO; n++) lagTime[n]=eventTime[n]-minTime;
104 for(n=0; n<
nIFO; n++) minTimeDet = (eventTime[n]<=minTime) ? n : minTimeDet;
107 double optRate = (NET->
getwc(k)->
cRate[ID-1])[0];
108 double optLayer = NET->
getwc(k)->
rate/optRate;
109 double optLevel =
int(log10(optLayer)/log10(2));
110 for(
int n=0;n<2;n++) {
120 sprintf(fname,
"%s/l_tfmap_scalogram%s.%s", outdir, strtime, gtype);
124 WTS.
hist2D->SetTitle(
"Scalogram");
129 sprintf(fname,
"%s/n_tfmap_scalogram%s.%s", outdir, strtime, gtype);
133 WTS.
hist2D->SetTitle(
"Scalogram");
138 sprintf(fname,
"%s/l_tfmap_scalogram%s.%s", outdir, strtime,
"txt");
141 sprintf(fname,
"%s/n_tfmap_scalogram%s.%s", outdir, strtime,
"txt");
147 gROOT->SetBatch(batch);
159 int nj =
int((t2-t1)*w.
rate())/ni;
169 if ( (fp = fopen(fname,
"w")) == NULL ) {
170 cout <<
" Dump() error: cannot open file " << fname <<
". \n";
174 for(
int i=0;
i<ni;
i++) {
176 for(
int j=nb;
j<=ne;
j++) {
177 float t = (
j+0.5)*
dt;
178 float f = (
i+0.5)*
df;
180 if(x>0.1)
fprintf( fp,
"%-3.6f \t%-3.6f \t%-3.3f \t%-3.3f \t%-3.3f\n", t, dt, f, df, x);
wavearray< double > t(hp.size())
detector * getifo(size_t n)
param: detector index
std::vector< vector_int > cRate
Float_t * rho
biased null statistics
WSeries< double > pixeLHood
Double_t * start
cluster duration = stopW-startW
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
virtual void rate(double r)
Float_t * duration
max cluster time relative to segment start
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
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...
virtual void start(double s)
std::vector< size_t > mdc__ID
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
void dopen(const char *fname, char *mode, bool header=true)
Double_t * stop
GPS start time of the cluster.
WSeries< double > pTF[nRES]
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Float_t * lag
time between consecutive events
void Dump(WSeries< double > &w, double t1, double t2, const char *fname)
Float_t * slag
time lag [sec]
void output(TTree *=NULL, network *=NULL, double=0., size_t=0, int=-1)
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
WSeries< double > pixeLNull
netevent EVT(itree, nifo)
Double_t * time
beam pattern coefficients for hx
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
WaveDWT< DataType_t > * pWavelet
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)