21 #pragma GCC system_header 28 #include "TObjArray.h" 29 #include "TObjString.h" 43 #define MDC_OTF // enable MDC On The Fly 44 #define NOISE_OTF // enable NOISE On The Fly 52 double& enINJ,
double& enREC,
double& xcorINJ_REC);
69 cout <<
"-----> CWB_Plugin_WRC.C" << endl;
70 cout <<
"ifo " << ifo.Data() << endl;
71 cout <<
"type " << type << endl;
91 if(ifo.CompareTo(
"L1")==0) seed=1000;
92 if(ifo.CompareTo(
"H1")==0) seed=2000;
93 if(ifo.CompareTo(
"V1")==0) seed=3000;
94 if(ifo.CompareTo(
"J1")==0) seed=4000;
95 if(ifo.CompareTo(
"A2")==0) seed=5000;
96 if(ifo.CompareTo(
"Y2")==0) seed=6000;
97 if(ifo.CompareTo(
"Y3")==0) seed=7000;
100 if(ifo.CompareTo(
"L1")==0) fName=
"plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
101 if(ifo.CompareTo(
"H1")==0) fName=
"plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
102 if(ifo.CompareTo(
"V1")==0) fName=
"plugins/strains/advVIRGO_sensitivity_12May09_8khz_one_side.txt";
103 if(ifo.CompareTo(
"J1")==0) fName=
"plugins/strains/LCGT_sensitivity_8khz_one_side.txt";
104 if(ifo.CompareTo(
"A2")==0) fName=
"plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
105 if(ifo.CompareTo(
"Y2")==0) fName=
"plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
106 if(ifo.CompareTo(
"Y3")==0) fName=
"plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
120 sprintf(cmd,
"network* net = (network*)%p;",NET);
121 gROOT->ProcessLine(cmd);
148 if(ifo.CompareTo(NET->
ifoName[0])==0) {
159 for(
int k=0;
k<(
int)NET->
mdcTime.size();
k++) cout <<
k <<
" mdcTime " << MDC.
mdcTime[
k] << endl;
160 for(
int k=0;
k<(
int)NET->
mdcType.size();
k++) cout <<
k <<
" mdcType " << MDC.
mdcType[
k] << endl;
167 cout <<
"CWB_Plugin_WRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
173 cout <<
"-----> CWB_Plugin_WRC.C -> " <<
" gIFACTOR : " << gIFACTOR << endl;
176 float* gSLAGSHIFT=NULL;
IMPORT(
float*,gSLAGSHIFT)
187 TList *files = (TList*)gROOT->GetListOfFiles();
194 while ((file=(TSystemFile*)
next())) {
195 fname = file->GetName();
197 if(fname.Contains(
"wave_")) {
198 froot=(TFile*)file;froot->cd();
200 outDump.ReplaceAll(
".root.tmp",
".txt");
201 cout <<
"output file name : " << fname << endl;
205 cout <<
"CWB_Plugin_WRC.C : Error - output root file not found" << endl;
209 cout <<
"CWB_Plugin_WRC.C : Error - output root file not found" << endl;
214 char ciSNR[16];
sprintf(ciSNR,
"_iSNR[%1d]/F",nIFO);
215 char coSNR[16];
sprintf(coSNR,
"_oSNR[%1d]/F",nIFO);
216 char cioSNR[16];
sprintf(cioSNR,
"_ioSNR[%1d]/F",nIFO);
221 TTree*
net_tree = (TTree *) froot->Get(
"waveburst");
224 #ifdef SAVE_EVT_CLUSTER 225 net_tree->SetBranchAddress(
"netcluster",&pwc);
228 net_tree->SetBranchAddress(
"_iSNR",iSNR);
229 net_tree->SetBranchAddress(
"_oSNR",oSNR);
230 net_tree->SetBranchAddress(
"_ioSNR",ioSNR);
235 #ifdef SAVE_EVT_CLUSTER 237 net_tree->Branch(
"netcluster",
"netcluster",&pwc,32000,0);
240 net_tree->Branch(
"_iSNR",iSNR,ciSNR);
241 net_tree->Branch(
"_oSNR",oSNR,coSNR);
242 net_tree->Branch(
"_ioSNR",ioSNR,cioSNR);
247 for(
int k=0;
k<
K;
k++) {
249 id = NET->
getwc(
k)->
get(const_cast<char*>(
"ID"), 0,
'L', rate);
253 int ID = size_t(
id.data[
j]+0.5);
266 double recTime = EVT->
time[0];
267 double injTh = EVT->
theta[1];
268 double injPh = EVT->
phi[1];
269 double recTh = EVT->
theta[0];
270 double recPh = EVT->
phi[0];
275 double injTime = 1.e12;
277 for(
int m=0;
m<
M;
m++) {
280 if(T<injTime && INJ.
fill_in(NET,mdcID)) {
291 int idSize = pd->
RWFID.size();
293 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
300 pwfINJ[
n] = INJ.
pwf[
n];
301 if (pwfINJ[
n]==NULL) {
302 cout <<
"CWB_Plugin_WRC.C : Error : Injected waveform not saved !!! : detector " 307 cout <<
"CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> " 308 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
312 if (wfIndex>=0) pwfREC[
n] = pd->
RWFP[wfIndex];
320 double enINJ, enREC, xcorINJ_REC;
332 ioSNR[
n] = xcorINJ_REC;
341 if(INJ.
fill_in(NET,-(injID+1))) {
346 int idSize = pd->
RWFID.size();
348 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==-ID) wfIndex=mm;
354 pwfINJ[
n] = INJ.
pwf[
n];
355 if (pwfINJ[
n]==NULL) {
356 cout <<
"CWB_Plugin_WRC.C : Error : Injected waveform not saved !!! : detector " 361 cout <<
"CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> " 362 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
366 if (wfIndex>=0) pwfREC[
n] = pd->
RWFP[wfIndex];
380 double enINJ, enREC, xcorINJ_REC;
388 ioSNR[
n] = xcorINJ_REC;
406 if(cfg->
dump) EVT->
dopen(outDump.Data(),
const_cast<char*
>(
"a"),
false);
407 EVT->
output2G(net_tree,NET,ID,
k,ofactor);
411 #ifdef SAVE_MRA_PLOT // monster event display 420 if(iSNR)
delete []
iSNR;
421 if(oSNR)
delete []
oSNR;
422 if(ioSNR)
delete []
ioSNR;
430 double& enINJ,
double& enREC,
double& xcorINJ_REC) {
432 double bINJ = wfINJ->
start();
433 double eINJ = wfINJ->
stop();
434 double bREC = wfREC->
start();
435 double eREC = wfREC->
stop();
440 double R=wfINJ->
rate();
442 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
443 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
446 double startXCOR = bINJ>bREC ? bINJ : bREC;
447 double endXCOR = eINJ<eREC ? eINJ : eREC;
448 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
455 if (sizeXCOR<=0)
return;
459 for (
int i=0;
i<sizeXCOR;
i++) enREC+=wfREC->
data[
i+oREC]*wfREC->
data[
i+oREC];
460 for (
int i=0;
i<sizeXCOR;
i++) xcorINJ_REC+=wfINJ->
data[
i+oINJ]*wfREC->
data[
i+oREC];
462 double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
463 cout <<
"enINJ : " << enINJ <<
" enREC : " << enREC
464 <<
" xcorINJ_REC : " << xcorINJ_REC <<
" enINJREC : " << enINJ+enREC-2*xcorINJ_REC << endl;
472 double f_low = cfg->
fLow;
473 double f_high = cfg->
fHigh;
474 cout <<
"f_low : " << f_low <<
" f_high : " << f_high << endl;
476 double Fs=((double)wf->
rate()/(double)wf->
size())/2.;
477 for (
int j=0;
j<wf->
size();
j+=2) {
479 if((f<f_low)||(f>f_high)) {wf->
data[
j]=0.;wf->
data[
j+1]=0.;}
492 sprintf(fname,
"l_tfmap_scalogram_%d.png",ID);
493 cout <<
"write " << fname << endl;
494 WTS.
plot(pwc, ID, nIFO,
'L', 0, const_cast<char*>(
"COLZ"));
502 watplot PTS(const_cast<char*>(
"ptswrc"),200,20,800,500);
509 PTS.
plot(wfINJ, const_cast<char*>(
"ALP"), 1, 0, 0,
true, cfg->
fLow, cfg->
fHigh);
511 PTS.
plot(wfINJ, const_cast<char*>(
"ALP"), 1, 0, 0);
513 PTS.
graph[0]->SetLineWidth(1);
515 PTS.
plot(wfREC, const_cast<char*>(
"SAME"), 2, 0, 0,
true, cfg->
fLow, cfg->
fHigh);
517 PTS.
plot(wfREC, const_cast<char*>(
"SAME"), 2, 0, 0);
519 PTS.
graph[1]->SetLineWidth(2);
524 if(fft)
sprintf(label,
"%s",
"fft");
525 else sprintf(label,
"%s",
"time");
526 if(strain)
sprintf(label,
"%s_%s",label,
"strain");
527 else sprintf(label,
"%s_%s",label,
"white");
530 sprintf(fname,
"%s_wf_%s_inj_rec_gps_%d.root",ifo.Data(),
label,
int(tmin));
532 cout <<
"write : " << fname << endl;
545 double optRate = (pwc->
cRate[ID-1])[0];
546 double optLayer = pwc->
rate/optRate;
547 double optLevel =
int(log10(optLayer)/log10(2));
549 double bINJ = wfINJ->
start();
550 double eINJ = wfINJ->
stop();
551 double bREC = wfREC->
start();
552 double eREC = wfREC->
stop();
557 double R=wfINJ->
rate();
559 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
560 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
563 double startXCOR = bINJ>bREC ? bINJ : bREC;
564 double endXCOR = eINJ<eREC ? eINJ : eREC;
565 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
569 if (sizeXCOR<=0)
return;
577 int xcor_length = sizeXCOR/wfINJ->
rate();
578 int wts_size = xcor_length<8 ? 16 : 16*
int(1+xcor_length/8.);
579 wts_size*=wfINJ->
rate();
583 double wts_start,wts_stop;
593 for(
int m=0;
m<sizeXCOR;
m++)
603 cout <<
"Print " << fname << endl;
606 wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/
wINJ.
rate() : wts_start+(wts_size/2)/
wINJ.
rate();
607 WTS.
plot(&
wINJ, 2, wts_start, wts_stop,const_cast<char*>(
"COLZ"));
608 WTS.
hist2D->GetYaxis()->SetRangeUser(EVT->
low[ifoID], EVT->
high[ifoID]);
621 for (
int m=0;
m<sizeXCOR;
m++)
631 cout <<
"Print " << fname << endl;
632 wts_start =
wREC.start()+(double)(wts_size/2)/
wREC.rate();
633 wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/
wREC.rate() : wts_start+(wts_size/2)/
wREC.rate();
634 WTS.
plot(&
wREC, 2, wts_start, wts_stop,const_cast<char*>(
"COLZ"));
635 WTS.
hist2D->GetYaxis()->SetRangeUser(EVT->
low[ifoID], EVT->
high[ifoID]);
648 for (
int m=0;
m<sizeXCOR;
m++)
656 cout <<
"Print " << fname << endl;
657 if(NET->
getwdm(optLayer+1)) wDIF.Forward(wDIF,*NET->
getwdm(optLayer+1));
658 wts_start = wDIF.start()+(double)(wts_size/2)/wDIF.rate();
659 wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wDIF.rate() : wts_start+(wts_size/2)/wDIF.rate();
660 WTS.
plot(&wDIF, 2, wts_start, wts_stop,const_cast<char*>(
"COLZ"));
661 WTS.
hist2D->GetYaxis()->SetRangeUser(EVT->
low[ifoID], EVT->
high[ifoID]);
671 cout <<
"TF : " <<
"enINJ : " << enINJ <<
" enREC : " << enREC
672 <<
" enINJREC : " << enINJREC << endl;
683 float dt = 1./(2*
df);
691 double pixel_frequency;
696 if(
j==1) pixel_frequency = df/2;
697 if(
j==layers) pixel_frequency = (layers-1)*df;
698 if((
j!=1)&&(
j!=
layers)) pixel_frequency = df/2 + (
j-1)*df;
700 if(
j==1) EE += ER/2.;
701 if(
j==layers) EE += ER/2.;
716 cout <<
"ComputeResidualEnergy - Error : WS1 & WS2 are not compatible !!!" << endl;
721 float dt = 1./(2*
df);
729 double pixel_frequency;
734 if(
j==1) pixel_frequency = df/2;
735 if(
j==layers) pixel_frequency = (layers-1)*df;
736 if((
j!=1)&&(
j!=
layers)) pixel_frequency = df/2 + (
j-1)*df;
738 if(
j==1) EE += ER/2.;
739 if(
j==layers) EE += ER/2.;
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
virtual void resize(unsigned int)
std::vector< vector_int > cRate
Float_t * high
min frequency
std::vector< wavearray< double > > wREC[MAX_TRIALS]
void setSLags(float *slag)
TString Get(wavearray< double > &x, TString ifo)
std::vector< double > * getmdcTime()
std::vector< std::string > mdcList
virtual void rate(double r)
Float_t * low
average center_of_snr frequency
std::vector< std::string > mdcType
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
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
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
std::vector< std::string > mdcType
std::vector< double > oSNR[NIFO_MAX]
virtual void start(double s)
std::vector< double > mdcTime
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< double > ioSNR[NIFO_MAX]
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
void dopen(const char *fname, char *mode, bool header=true)
WDM< double > * getwdm(size_t M)
param: number of wdm layers
virtual size_t size() const
void output2G(TTree *, network *, size_t, int, double)
void ComputeResidualEnergyOptTF(wavearray< double > *wfINJ, wavearray< double > *wfREC, network *NET, netevent *EVT, int lag, int ID, int ifoID)
#define IMPORT(TYPE, VAR)
Double_t * gps
average center_of_gravity time
std::vector< std::string > mdcList
TIter next(twave->GetListOfBranches())
std::vector< double > iSNR[NIFO_MAX]
void PlotWaveform(TString ifo, wavearray< double > *wfINJ, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
size_t cpf(const netcluster &, bool=false, int=0)
void ComputeResidualEnergy(wavearray< double > *wfINJ, wavearray< double > *wfREC, double &enINJ, double &enREC, double &xcorINJ_REC)
Float_t * phi
sqrt(h+*h+ + hx*hx)
netevent EVT(itree, nifo)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Double_t * time
beam pattern coefficients for hx
virtual void stop(double s)
double fabs(const Complex &x)
void ApplyFrequencyCuts(wavearray< double > *wf, CWB::config *cfg)
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
void PlotMRA(netcluster *pwc, int nIFO, int ID)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double factors[FACTORS_MAX]
DataType_t getSample(int n, double m)
wavearray< double > wINJ[NIFO_MAX]
std::vector< double > mdcTime
Bool_t fill_in(network *, int, bool=true)
size_t getmdc__ID(size_t n)
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
virtual void resize(unsigned int)
double ComputeEnergy(WSeries< double > *WS)