3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
22 #define nPEAK 10 // number of peaks to be extracted from the reconstructed waveform
25 float PB[][nPEAK],
float PA[][nPEAK],
float PE[][nPEAK]);
36 cout <<
"-----> CWB_Plugin_WavePeaks.C" << endl;
37 cout <<
"ifo " << ifo.Data() << endl;
38 cout <<
"type " << type << endl;
58 cout <<
"CWB_Plugin_WavePeaks.C -> "
59 <<
"CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
65 cout <<
"-----> CWB_Plugin_WavePeaks.C -> "
66 <<
" gIFACTOR : " << gIFACTOR << endl;
69 float* gSLAGSHIFT=NULL;
IMPORT(
float*,gSLAGSHIFT)
80 TList *
files = (TList*)gROOT->GetListOfFiles();
87 while ((file=(TSystemFile*)
next())) {
88 fname = file->GetName();
90 if(fname.Contains(
"wave_")) {
91 froot=(TFile*)file;froot->cd();
93 outDump.ReplaceAll(
".root.tmp",
".txt");
98 cout <<
"CWB_Plugin_WavePeaks.C : Error - output root file not found" << endl;
102 cout <<
"CWB_Plugin_WavePeaks.C : Error - output root file not found" << endl;
106 TTree*
net_tree = (TTree *) froot->Get(
"waveburst");
109 net_tree->SetBranchAddress(
"nP",&nP);
110 net_tree->SetBranchAddress(
"PF",PF);
111 net_tree->SetBranchAddress(
"PB",PB);
112 net_tree->SetBranchAddress(
"PA",PA);
113 net_tree->SetBranchAddress(
"PE",PE);
117 net_tree->Branch(
"nP",&nP,
"nP/I");
118 net_tree->Branch(
"PF",PF,TString::Format(
"PF[%i][%i]/F",cfg->
nIFO,
nPEAK));
119 net_tree->Branch(
"PB",PB,TString::Format(
"PB[%i][%i]/F",cfg->
nIFO,
nPEAK));
120 net_tree->Branch(
"PA",PA,TString::Format(
"PA[%i][%i]/F",cfg->
nIFO,
nPEAK));
121 net_tree->Branch(
"PE",PE,TString::Format(
"PE[%i][%i]/F",cfg->
nIFO,
nPEAK));
125 for(
int k=0;
k<
K;
k++) {
131 int ID = size_t(
id.data[
j]+0.5);
144 int idSize = pd->
RWFID.size();
147 for (
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
148 if(wfIndex==-1)
continue;
155 pwfREC[
n] = pd->
RWFP[wfIndex];
178 if(cfg->
dump) EVT->
dopen(outDump.Data(),
const_cast<char*
>(
"a"),
false);
179 EVT->
output2G(net_tree,NET,ID,
k,ofactor);
193 float PB[][nPEAK],
float PA[][nPEAK],
float PE[][nPEAK]) {
196 double Fs=((double)wf->
rate()/(double)(2*size));
207 for(
int j=0;
j<
size;
j++) ET+=ener[
j];
210 for(
int j=0;j<size;j++) if(ener[j]>MA) MA=ener[
j];
218 for(
int j=0;j<size;j++) if(ener[j]>ME) {ME=ener[
j];jME=
j;}
221 float fmean=ME*Fs*jME;
222 float frms=ME*pow(Fs*jME,2);
227 for(
int j=jME-1;
j>=0;
j--) {
230 if(E<ge) {EP+=E;finf=
F;fmean+=E*
F;frms+=E*F*
F;ener[
j]=0;}
else break;
238 if(E<ge) {EP+=E;fsup=
F;fmean+=E*
F;frms+=E*F*
F;ener[
j]=0;}
else break;
243 frms = frms/EP-fmean*fmean;
244 frms = frms>0 ? sqrt(frms) : 0;
247 PA[
ifoID][
k] = sqrt(ME)/MA;
249 PE[
ifoID][
k] = ET>0 ? EP/ET : 0;
260 double tmin = wfREC->
start();
267 PTS.
graph[0]->SetLineWidth(1);
271 if(fft)
sprintf(label,
"%s",
"fft");
272 else sprintf(label,
"%s",
"time");
273 if(strain)
sprintf(label,
"%s_%s",label,
"strain");
274 else sprintf(label,
"%s_%s",label,
"white");
279 cout <<
"write : " << fname << endl;
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
virtual size_t size() const
void setSLags(float *slag)
virtual void rate(double r)
void PlotWaveform(TString ifo, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
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...
std::vector< wavearray< double > * > RWFP
std::vector< TGraph * > graph
virtual void start(double s)
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 dopen(const char *fname, char *mode, bool header=true)
void GetGlitchParams(wavearray< double > *wf, int ifoID, float PF[][nPEAK], float PB[][nPEAK], float PA[][nPEAK], float PE[][nPEAK])
void output2G(TTree *, network *, size_t, int, double)
#define IMPORT(TYPE, VAR)
TIter next(twave->GetListOfBranches())
netevent EVT(itree, nifo)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double factors[FACTORS_MAX]
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.