3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
25 #define SAVE_TF_RES_ENERGY //save TF residual energy (reconstructed -injected)
28 #define SAVE_EVT_CLUSTER // save cluster to the event output tree
30 void ComputeResidualEnergyOptTF(
wavearray<double>* wfREC,
network*
NET,
netevent*
EVT,
int lag,
int ID,
int ifoID,
double& entf_REC,
double& TFtime_REC,
double& sig_TFtime_REC,
double& TFfreq_REC,
double& sig_TFfreq_REC,
double& TFduration,
double& TFband);
45 cout <<
"-----> CWB_Plugin_WRC.C" << endl;
46 cout <<
"ifo " << ifo.Data() << endl;
47 cout <<
"type " << type << endl;
67 sprintf(cmd,
"network* net = (network*)%p;",NET);
68 gROOT->ProcessLine(cmd);
95 if(ifo.CompareTo(NET->
ifoName[0])==0) {
105 for(
int k=0;
k<(
int)NET->
mdcList.size();
k++) cout <<
k <<
" mdcList " << MDC\
107 for(
int k=0;
k<(
int)NET->
mdcTime.size();
k++) cout <<
k <<
" mdcTime " << MDC\
109 for(
int k=0;
k<(
int)NET->
mdcType.size();
k++) cout <<
k <<
" mdcType " << MDC\
119 cout <<
"CWB_Plugin_WRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
130 cout <<
"-----> CWB_Plugin_WRC.C -> " <<
" gIFACTOR : " << gIFACTOR << endl;
146 TList *
files = (TList*)gROOT->GetListOfFiles();
153 while ((file=(TSystemFile*)
next())) {
154 fname = file->GetName();
156 if(fname.Contains(
"wave_")) {
157 froot=(TFile*)file;froot->cd();
159 outDump.ReplaceAll(
".root.tmp",
".txt");
160 cout <<
"output file name : " << fname << endl;
164 cout <<
"CWB_Plugin_WRC.C : Error - output root file not found" << endl;
168 cout <<
"CWB_Plugin_WRC.C : Error - output root file not found" << endl;
173 char cEnTf_REC[24];
sprintf(cEnTf_REC,
"_EnTf_REC[%1d]/D",nIFO);
174 double* EnTf_REC =
new double[
nIFO];
175 char cEnTfWh_REC[24];
sprintf(cEnTfWh_REC,
"_EnTfWh_REC[%1d]/D",nIFO);
176 double* EnTfWh_REC =
new double[
nIFO];
180 char cTFtauREC_tot[24];
sprintf(cTFtauREC_tot,
"_TFtauREC_tot[%1d]/D",nIFO);
181 char cTFfreqREC_tot[24];
sprintf(cTFfreqREC_tot,
"_TFfreqREC_tot[%1d]/D",nIFO);
182 char cTFtauWhREC_tot[24];
sprintf(cTFtauWhREC_tot,
"_TFtauWhREC_tot[%1d]/D",nIFO);
183 char cTFfreqWhREC_tot[24];
sprintf(cTFfreqWhREC_tot,
"_TFfreqWhREC_tot[%1d]/D",nIFO);
185 char csig_TFtauREC_tot[48];
sprintf(csig_TFtauREC_tot,
"_sig_TFtauREC_tot[%1d]/D",nIFO);
186 char csig_TFfreqREC_tot[48];
sprintf(csig_TFfreqREC_tot,
"_sig_TFfreqREC_tot[%1d]/D",nIFO);
187 char csig_TFtauWhREC_tot[48];
sprintf(csig_TFtauWhREC_tot,
"_sig_TFtauWhREC_tot[%1d]/D",nIFO);
188 char csig_TFfreqWhREC_tot[48];
sprintf(csig_TFfreqWhREC_tot,
"_sig_TFfreqWhREC_tot[%1d]/D",nIFO);
192 char cTFbandREC_tot[24];
sprintf(cTFbandREC_tot,
"_TFbandREC_tot[%1d]/D",nIFO);
193 char cTFdurationREC_tot[24];
sprintf(cTFdurationREC_tot,
"_TFdurationREC_tot[%1d]/D",nIFO);
194 char cTFbandWhREC_tot[24];
sprintf(cTFbandWhREC_tot,
"_TFbandWhREC_tot[%1d]/D",nIFO);
195 char cTFdurationWhREC_tot[24];
sprintf(cTFdurationWhREC_tot,
"_TFdurationWhREC_tot[%1d]/D",nIFO);
198 double* TFtauREC_tot =
new double[
nIFO];
199 double* TFfreqREC_tot =
new double[
nIFO];
200 double* TFtauWhREC_tot =
new double[
nIFO];
201 double* TFfreqWhREC_tot =
new double[
nIFO];
203 double* sig_TFtauREC_tot =
new double[
nIFO];
204 double* sig_TFfreqREC_tot =
new double[
nIFO];
205 double* sig_TFtauWhREC_tot =
new double[
nIFO];
206 double* sig_TFfreqWhREC_tot =
new double[
nIFO];
210 double* TFbandREC_tot =
new double[
nIFO];
211 double* TFdurationREC_tot =
new double[
nIFO];
212 double* TFbandWhREC_tot =
new double[
nIFO];
213 double* TFdurationWhREC_tot =
new double[
nIFO];
216 TTree*
net_tree = (TTree *) froot->Get(
"waveburst");
219 #ifdef SAVE_EVT_CLUSTER
220 net_tree->SetBranchAddress(
"netcluster",&pwc);
223 #ifdef SAVE_TF_RES_ENERGY
224 net_tree->SetBranchAddress(
"_EnTf_REC", EnTf_REC);
225 net_tree->SetBranchAddress(
"_EnTfWh_REC", EnTfWh_REC);
227 net_tree->SetBranchAddress(
"_TFtauREC_tot",TFtauREC_tot);
228 net_tree->SetBranchAddress(
"_TFfreqREC_tot",TFfreqREC_tot);
229 net_tree->SetBranchAddress(
"_TFtauWhREC_tot",TFtauWhREC_tot);
230 net_tree->SetBranchAddress(
"_TFfreqWhREC_tot",TFfreqWhREC_tot);
232 net_tree->SetBranchAddress(
"_sig_TFtauREC_tot",sig_TFtauREC_tot);
233 net_tree->SetBranchAddress(
"_sig_TFfreqREC_tot",sig_TFfreqREC_tot);
234 net_tree->SetBranchAddress(
"_sig_TFtauWhREC_tot",sig_TFtauWhREC_tot);
235 net_tree->SetBranchAddress(
"_sig_TFfreqWhREC_tot",sig_TFfreqWhREC_tot);
237 net_tree->SetBranchAddress(
"_TFbandREC_tot",TFbandREC_tot);
238 net_tree->SetBranchAddress(
"_TFdurationREC_tot",TFdurationREC_tot);
239 net_tree->SetBranchAddress(
"_TFbandWhREC_tot",TFbandWhREC_tot);
240 net_tree->SetBranchAddress(
"_TFdurationWhREC_tot",TFdurationWhREC_tot);
246 #ifdef SAVE_EVT_CLUSTER
248 net_tree->Branch(
"netcluster",
"netcluster",&pwc,32000,0);
251 #ifdef SAVE_TF_RES_ENERGY
252 net_tree->Branch(
"_EnTf_REC", EnTf_REC,cEnTf_REC);
253 net_tree->Branch(
"_EnTfWh_REC", EnTfWh_REC,cEnTfWh_REC);
255 net_tree->Branch(
"_TFtauREC_tot",TFtauREC_tot, cTFtauREC_tot);
256 net_tree->Branch(
"_TFfreqREC_tot",TFfreqREC_tot,cTFfreqREC_tot);
257 net_tree->Branch(
"_TFtauWhREC_tot",TFtauWhREC_tot,cTFtauWhREC_tot);
258 net_tree->Branch(
"_TFfreqWhREC_tot",TFfreqWhREC_tot,cTFfreqWhREC_tot);
260 net_tree->Branch(
"_sig_TFtauREC_tot",sig_TFtauREC_tot, csig_TFtauREC_tot);
261 net_tree->Branch(
"_sig_TFfreqREC_tot",sig_TFfreqREC_tot,csig_TFfreqREC_tot);
262 net_tree->Branch(
"_sig_TFtauWhREC_tot",sig_TFtauWhREC_tot,csig_TFtauWhREC_tot);
263 net_tree->Branch(
"_sig_TFfreqWhREC_tot",sig_TFfreqWhREC_tot,csig_TFfreqWhREC_tot);
265 net_tree->Branch(
"_TFbandREC_tot",TFbandREC_tot, cTFbandREC_tot);
266 net_tree->Branch(
"_TFdurationREC_tot",TFdurationREC_tot,cTFdurationREC_tot);
267 net_tree->Branch(
"_TFbandWhREC_tot",TFbandWhREC_tot,cTFbandWhREC_tot);
268 net_tree->Branch(
"_TFdurationWhREC_tot",TFdurationWhREC_tot,cTFdurationWhREC_tot);
277 for(
int k=0;
k<
K;
k++) {
281 int ID = size_t(
id.data[
j]+0.5);
290 double recTime = EVT->
time[0];
291 double recTh = EVT->
theta[0];
292 double recPh = EVT->
phi[0];
298 int idSize = pd->
RWFID.size();
300 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
304 cout <<
"CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> "
305 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
310 if (wfIndex>=0) pwfREC[
n] = pd->
RWFP[wfIndex];
314 double entf_REC, tftime_REC,sig_tftime_REC,tffreq_REC,sig_tffreq_REC, tfband_rec, tfduration_rec;
317 ComputeResidualEnergyOptTF( wfREC, NET,EVT,
k, ID,
n, entf_REC, tftime_REC, sig_tftime_REC, tffreq_REC, sig_tffreq_REC,tfduration_rec, tfband_rec);
318 EnTfWh_REC[
n] = entf_REC;
319 TFtauWhREC_tot[
n] = wfREC->
start()+tftime_REC;
320 TFfreqWhREC_tot[
n] = tffreq_REC;
321 sig_TFtauWhREC_tot[
n] = sig_tftime_REC;
322 sig_TFfreqWhREC_tot[
n] = sig_tffreq_REC;
323 TFbandWhREC_tot[
n] = tfband_rec;
324 TFdurationWhREC_tot[
n] = tfduration_rec;
325 cout<<
"WHITE "<<tfduration_rec<<
" "<<tfband_rec<<endl;
334 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==-ID) wfIndex=mm;
339 cout <<
"CWB_Plugin_WRC.C : Error : Reconstructed waveform not \
341 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
345 if (wfIndex>=0) pwfREC_s[
n] = pd->
RWFP[wfIndex];
349 double entf_REC, tftime_REC,sig_tftime_REC,tffreq_REC,sig_tffreq_REC, tfband_rec, tfduration_rec;
350 cout<<
"STRAIN "<<endl;
353 ComputeResidualEnergyOptTF( wfREC_s, NET,EVT,
k, ID,
n, entf_REC, tftime_REC, sig_tftime_REC, tffreq_REC, sig_tffreq_REC,tfduration_rec, tfband_rec);
356 EnTf_REC[
n] = entf_REC;
357 TFtauREC_tot[
n] = wfREC_s->
start()+tftime_REC;
358 TFfreqREC_tot[
n] = tffreq_REC;
359 sig_TFtauREC_tot[
n] = sig_tftime_REC;
360 sig_TFfreqREC_tot[
n] = sig_tffreq_REC;
361 TFbandREC_tot[
n] = tfband_rec;
362 TFdurationREC_tot[
n] = tfduration_rec;
381 if(cfg->
dump) EVT->
dopen(outDump.Data(),
const_cast<char*
>(
"a"),
false);
382 EVT->
output2G(net_tree,NET,ID,
k,ofactor);
395 if( TFtauREC_tot)
delete [] TFtauREC_tot;
396 if( TFfreqREC_tot)
delete [] TFfreqREC_tot;
397 if( TFtauWhREC_tot)
delete [] TFtauWhREC_tot;
398 if( TFfreqWhREC_tot)
delete [] TFfreqWhREC_tot;
412 void ComputeResidualEnergyOptTF(
wavearray<double>* wfREC,
network*
NET,
netevent*
EVT,
int lag,
int ID,
int ifoID,
double& entf_REC,
double& TFtime_REC,
double& sig_TFtime_REC,
double& TFfreq_REC,
double& sig_TFfreq_REC,
double& TFduration,
double& TFband)
418 double optRate = (pwc->
cRate[ID-1])[0];
419 double optLayer = pwc->
rate/optRate;
420 double optLevel =
int(log10(optLayer)/log10(2));
422 double bREC = wfREC->
start();
423 double eREC = wfREC->
stop();
425 double R=wfREC->
rate();
434 int sizeXCOR =
int((eREC-bREC)*R+0.5);
437 int xcor_length = sizeXCOR/wfREC->
rate();
438 int wts_size =xcor_length<8 ? 16 : 16*
int(1+xcor_length/8.);
439 wts_size*=wfREC->
rate();
441 cout<<
"SIXX "<<sizeXCOR<<
" "<<wts_size<<endl;
445 double wts_start,wts_stop;
459 for (
int m=0;
m<sizeXCOR;
m++)
465 cout<<
"LLLLL "<<
wREC.start()+(double)(wts_size/2)/
wREC.rate() <<
" "<<wts_start+sizeXCOR/
wREC.rate()<<endl;
468 cout <<
"TF : " <<
"enINJ : "<<
" enREC : " << entf_REC
470 double tm, fq,stm,sfq, bd, dr;
480 cout<<
"TEST freqtime"<<TFfreq_REC<<
" "<< sig_TFfreq_REC<<
" TFtime " <<TFtime_REC<<
" band "<<TFband<< endl;
494 float dt = 1./(2*
df);
502 double pixel_frequency;
507 if(
j==1) pixel_frequency = df/2;
508 if(
j==layers) pixel_frequency = (layers-1)*df;
509 if((
j!=1)&&(
j!=
layers)) pixel_frequency = df/2 + (
j-1)*df;
511 if(
j==1) EE += ER/2.;
512 if(
j==layers) EE += ER/2.;
529 float dt = 1./(2*
df);
536 double timeEn =0;
double freqEn=0;
537 double EE=0;
double EE4=0;
538 double pixel_frequency;
547 if(
j==1) {pixel_frequency = df/4; }
548 if(
j==layers) {pixel_frequency = (layers-1)*df -df/4;}
549 if((
j!=1)&&(
j!=
layers)) pixel_frequency = df/2 + (
j-1)*df +df/2;
552 f = en*pixel_frequency;
554 if(
j==1) { EE += en/2.; EE4 += (pow(en,2))/2.; timeEn += t/2; freqEn +=t/2;}
555 if(
j==layers){ EE += en/2.;EE4 += (pow(en,2))/2.; timeEn +=t/2; freqEn +=t/2 ;}
556 if((
j!=1)&&(
j!=
layers)){ EE += en; EE4 += (pow(en,2));timeEn +=
t; freqEn +=
f;}
563 sig_TFtime = dt*sqrt(1./12.)*sqrt(EE4)/EE;
565 sig_TFfreq= df*sqrt(1./12.)*sqrt(EE4)/EE;
567 double bandEn =0;
double durationEn=0;
575 if(
j==1) {pixel_frequency = df/4; }
576 if(
j==layers) {pixel_frequency = (layers-1)*df -df/4;}
577 if((
j!=1)&&(
j!=
layers)) pixel_frequency = df/2 + (
j-1)*df +df/2;
580 band = en*(pixel_frequency-TFfreq)*(pixel_frequency-TFfreq);
581 duration = en*(pixel_time-TFtime)*(pixel_time-TFtime);
582 if(
j==1) { EE += en/2.;durationEn +=duration/2; bandEn += band/2; }
583 if(
j==layers){ EE += en/2.; durationEn +=duration/2; bandEn +=band/2 ;}
584 if((
j!=1)&&(
j!=layers)){ EE += en; durationEn +=
duration; bandEn +=band;}
589 TFband = sqrt(bandEn/EE);
590 TFduration= sqrt(durationEn/EE);
591 cout<<
" BBB DDF "<<TFband <<
" "<<TFduration<<endl;
wavearray< double > t(hp.size())
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
virtual size_t size() const
std::vector< vector_int > cRate
std::vector< wavearray< double > > wREC[MAX_TRIALS]
TString Get(wavearray< double > &x, TString ifo)
std::vector< std::string > mdcList
virtual void rate(double r)
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
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<< "\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< std::string > mdcType
virtual void start(double s)
std::vector< double > mdcTime
void dopen(const char *fname, char *mode, bool header=true)
void ComputeTFtimefreq(WSeries< double > *WS, double &TFtime, double &sig_TFtime, double &TFfreq, double &sig_TFfreq, double &TFduration, double &TFband)
WDM< double > * getwdm(size_t M)
param: number of wdm layers
void output2G(TTree *, network *, size_t, int, double)
#define IMPORT(TYPE, VAR)
Double_t * gps
average center_of_gravity time
std::vector< std::string > mdcList
TIter next(twave->GetListOfBranches())
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
double ComputeEnergy(WSeries< double > *WS)
size_t cpf(const netcluster &, bool=false, int=0)
void ComputeResidualEnergyOptTF(wavearray< double > *wfREC, network *NET, netevent *EVT, int lag, int ID, int ifoID, double &entf_REC, double &TFtime_REC, double &sig_TFtime_REC, double &TFfreq_REC, double &sig_TFfreq_REC, double &TFduration, double &TFband)
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)
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double factors[FACTORS_MAX]
DataType_t getSample(int n, double m)
std::vector< double > mdcTime
virtual void resize(unsigned int)