Logo coherent WaveBurst  
Library Reference Guide
Logo
cWB_Plugin_PE_on.C
Go to the documentation of this file.
1 #define XIFO 4
2 
3 #pragma GCC system_header
4 
5 #include "cwb.hh"
6 #include "config.hh"
7 #include "network.hh"
8 #include "wavearray.hh"
9 #include "TString.h"
10 #include "TObjArray.h"
11 #include "TObjString.h"
12 #include "TRandom.h"
13 #include "TComplex.h"
14 #include "TMath.h"
15 #include "mdc.hh"
16 #include "watplot.hh"
17 #include "gwavearray.hh"
18 #include <vector>
19 
20 //#define SAVE_MRA_PLOT // enable event MRA plots
21 //#define SAVE_WHT_PLOT // enable event WHITE plots
22 //#define SAVE_STR_PLOT // enable event STRAIN plots
23 //#define SAVE_TF_PLOT // enable event WHITE TF plots
24 
25 #define SAVE_TF_RES_ENERGY //save TF residual energy (reconstructed -injected)
26 
27 
28 #define SAVE_EVT_CLUSTER // save cluster to the event output tree
29 
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);
31 
32 double ComputeEnergy(WSeries<double> *WS);
33 
34 void ComputeTFtimefreq(WSeries<double> *WS, double& TFtime, double& sig_TFtime, double& TFfreq, double& sig_TFfreq,double& TFduration, double& TFband);
35 
36 
37 void
39 //!NOISE_MDC_SIMULATION
40 // Extract whitened/strain injected & reconstructed waveforms, compute residual energy
41 // Save residual energy to the output wave root file
42 // This plugin can be use for Waveform Reconstruction Challenge (WRC) studies
43 
44  cout << endl;
45  cout << "-----> CWB_Plugin_WRC.C" << endl;
46  cout << "ifo " << ifo.Data() << endl;
47  cout << "type " << type << endl;
48  cout << endl;
49 
50  if(type==CWB_PLUGIN_CONFIG) {
51  cfg->outPlugin=true; // disable built-in output root file
52 
53  // gCUTS.resize(cfg->lagSize); // resize gCUTS
54  }
55 
56  if(type==CWB_PLUGIN_ILIKELIHOOD) {
57 
58  NET->wfsave=true; // enable waveform save
59 
60  }
61 
62 
63  cfg->mdcPlugin=true;
64  if(type==CWB_PLUGIN_MDC) {
65 
66  char cmd[128];
67  sprintf(cmd,"network* net = (network*)%p;",NET);
68  gROOT->ProcessLine(cmd);
69 
70  CWB::mdc MDC(NET);
71 
72  // ---------------------------------
73  // read plugin config
74  // ---------------------------------
75 
76  cfg->configPlugin.Exec();
77 
78  // ---------------------------------
79  // set list of mdc waveforms
80  // ---------------------------------
81 
82  IMPORT(CWB::mdc,MDC)
83  MDC.Print();
84 
85 // ---------------------------------
86 // get mdc data
87 // ---------------------------------
88 
89  MDC.Get(*x,ifo);
90 
91  // ---------------------------------
92  // set mdc list in the network class
93  // ---------------------------------
94 
95  if(ifo.CompareTo(NET->ifoName[0])==0) {
96  NET->mdcList.clear();
97  NET->mdcType.clear();
98  NET->mdcTime.clear();
99  NET->mdcList=MDC.mdcList;
100  NET->mdcType=MDC.mdcType;
101  NET->mdcTime=MDC.mdcTime;
102  }
103 
104  cout.precision(14);
105  for(int k=0;k<(int)NET->mdcList.size();k++) cout << k << " mdcList " << MDC\
106  .mdcList[k] << endl;
107  for(int k=0;k<(int)NET->mdcTime.size();k++) cout << k << " mdcTime " << MDC\
108  .mdcTime[k] << endl;
109  for(int k=0;k<(int)NET->mdcType.size();k++) cout << k << " mdcType " << MDC\
110  .mdcType[k] << endl;
111 }
112 
113 
114 
115 
116  if(type==CWB_PLUGIN_OLIKELIHOOD) {
117 
118  if(TString(cfg->analysis)!="2G") {
119  cout << "CWB_Plugin_WRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
120  gSystem->Exit(1);
121  }
122 
123 
124  // bool sim = false;
125  //if(cfg->simulation>0) sim=true;
126 
127 
128  // import ifactor
129  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
130  cout << "-----> CWB_Plugin_WRC.C -> " << " gIFACTOR : " << gIFACTOR << endl;
131 
132  // if(gIFACTOR!=sIFACTOR) { // reset gCUTS when gIFACTOR change
133  // for(int i=0;i<gCUTS.size();i++) gCUTS[i].clear();
134  // sIFACTOR=gIFACTOR;
135  // }
136 
137  int nIFO = NET->ifoListSize(); // number of detectors
138  int K = NET->nLag; // number of time lag
139  netevent* EVT;
141  double factor = cfg->simulation==3||cfg->simulation==4 ? -gIFACTOR : cfg->factors[gIFACTOR];
142  int rate = 0; // select all resolutions
143 
144  // search output root file in the system list
145  TFile* froot = NULL;
146  TList *files = (TList*)gROOT->GetListOfFiles();
147  TString outDump="";
148  if (files) {
149  TIter next(files);
150  TSystemFile *file;
151  TString fname;
152  bool check=false;
153  while ((file=(TSystemFile*)next())) {
154  fname = file->GetName();
155  // set output root file as the current file
156  if(fname.Contains("wave_")) {
157  froot=(TFile*)file;froot->cd();
158  outDump=fname;
159  outDump.ReplaceAll(".root.tmp",".txt");
160  cout << "output file name : " << fname << endl;
161  }
162  }
163  if(!froot) {
164  cout << "CWB_Plugin_WRC.C : Error - output root file not found" << endl;
165  gSystem->Exit(1);
166  }
167  } else {
168  cout << "CWB_Plugin_WRC.C : Error - output root file not found" << endl;
169  gSystem->Exit(1);
170  }
171 
172  netcluster* pwc = new netcluster();
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];
177 
178 
179  // leaf
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);
184 
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);
189 
190 
191 
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);
196 
197 
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];
202 
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];
207 
208 
209 
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];
214 
215 
216  TTree* net_tree = (TTree *) froot->Get("waveburst");
217  if(net_tree!=NULL) {
218  EVT = new netevent(net_tree,nIFO);
219 #ifdef SAVE_EVT_CLUSTER
220  net_tree->SetBranchAddress("netcluster",&pwc);
221 #endif
222 
223 #ifdef SAVE_TF_RES_ENERGY
224  net_tree->SetBranchAddress("_EnTf_REC", EnTf_REC);
225  net_tree->SetBranchAddress("_EnTfWh_REC", EnTfWh_REC);
226 
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);
231 
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);
236 
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);
241  #endif
242 
243  } else {
244  EVT = new netevent(nIFO);
245  net_tree = EVT->setTree();
246  #ifdef SAVE_EVT_CLUSTER
247  // add new netcluster leaf
248  net_tree->Branch("netcluster","netcluster",&pwc,32000,0);
249  #endif
250 
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);
254 
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);
259 
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);
264 
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);
269 
270 
271 #endif
272 
273 
274  }
275 
276 
277  for(int k=0; k<K; k++) { // loop over the lags
278 
279  id = NET->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
280  for(int j=0; j<(int)id.size(); j++) { // loop over cluster index
281  int ID = size_t(id.data[j]+0.5);
282  if(NET->getwc(k)->sCuts[ID-1]!=-1) continue; // skip rejected/processed clusters
283  double ofactor=0;
284  if(cfg->simulation==4) ofactor=-factor;
285  else if(cfg->simulation==3) ofactor=-gIFACTOR;
286  else ofactor=factor;
287 
288  if(k==0) { // only for zero lag
289  EVT->output2G(NULL,NET,ID,k,ofactor); // get reconstructed parameters
290  double recTime = EVT->time[0]; // rec event time det=0
291  double recTh = EVT->theta[0]; // rec theta
292  double recPh = EVT->phi[0]; // rec phi
293 
294 
295 
296  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
297  detector* pd = NET->getifo(0);
298  int idSize = pd->RWFID.size();
299  int wfIndex=-1;
300  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
301  // extract whitened injected & reconstructed waveforms
302  for(int n=0; n<nIFO; n++) {
303  if (wfIndex<0) {
304  cout << "CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> "
305  << ID << " : detector " << NET->ifoName[n] << endl;
306  continue;
307  }
308 
309  pd = NET->getifo(n);
310  if (wfIndex>=0) pwfREC[n] = pd->RWFP[wfIndex];
311  double R = pd->TFmap.rate();
312  wavearray<double>* wfREC = pwfREC[n]; // array of reconstructed waveforms
313 
314  double entf_REC, tftime_REC,sig_tftime_REC,tffreq_REC,sig_tffreq_REC, tfband_rec, tfduration_rec;
315 
316  //compute residual energy in time domain
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;
326  }
327  delete [] pwfREC;
328 
329 
330  wavearray<double>** pwfREC_s = new wavearray<double>*[nIFO];
331  // detector* pd = NET->getifo(0);
332  // int idSize = pd->RWFID.size();
333  //int wfIndex=-1;
334  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==-ID) wfIndex=mm;
335 
336  for(int n=0; n<nIFO; n++) {
337  pd = NET->getifo(n);
338  if (wfIndex<0) {
339  cout << "CWB_Plugin_WRC.C : Error : Reconstructed waveform not \
340 saved !!! : ID -> "
341  << ID << " : detector " << NET->ifoName[n] << endl;
342  continue;
343  }
344 
345  if (wfIndex>=0) pwfREC_s[n] = pd->RWFP[wfIndex];
346  double R = pd->TFmap.rate();
347  wavearray<double>* wfREC_s = pwfREC_s[n]; // array of reconstructed waveforms
348 
349  double entf_REC, tftime_REC,sig_tftime_REC,tffreq_REC,sig_tffreq_REC, tfband_rec, tfduration_rec;
350  cout<<"STRAIN "<<endl;
351 
352  //compute residual energy in time domain
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);
354 
355 
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;
363  }
364 
365 
366  delete [] pwfREC_s;
367 
368  }
369 
370 
371 
372  std::vector<int> sCuts = NET->getwc(k)->sCuts; // save sCuts
373  // set sCuts=1 for the events which must be not copied with cps to pwc
374  for(int i=0; i<(int)sCuts.size(); i++) if(i!=ID-1) NET->getwc(k)->sCuts[i]=1;
375 
376  // after cpf, pwc contains only one event
377  // ID can not be used to get the event, to get event use ID=1 (ex: for watplot)
378  pwc->cpf(*(NET->getwc(k))); // note: likelihood2G delete tdAmp
379  NET->getwc(k)->sCuts = sCuts; // restore sCuts
380 
381  if(cfg->dump) EVT->dopen(outDump.Data(),const_cast<char*>("a"),false);
382  EVT->output2G(net_tree,NET,ID,k,ofactor); // get reconstructed parameters
383  if(cfg->dump) EVT->dclose();
384 
385  if(!cfg->cedDump) NET->getwc(k)->sCuts[ID-1]=1; // mark as processed
386  }
387 
388  }
389  jfile->cd();
390 
391  if(pwc) delete pwc;
392  if(EVT) delete EVT;
393 
394 
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;
399 
400  }
401 
402 
403  return;
404 
405 
406 
407 
408 }
409 
410 
411 
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)
413  {
414 
415  entf_REC=0;
416  // find TF optimal resolution level
417  netcluster* pwc = NET->getwc(lag);
418  double optRate = (pwc->cRate[ID-1])[0];
419  double optLayer = pwc->rate/optRate;
420  double optLevel = int(log10(optLayer)/log10(2));
421 
422  double bREC = wfREC->start();
423  double eREC = wfREC->stop();
424 
425  double R=wfREC->rate();
426  detector* pd = NET->getifo(0);
427 
428  // int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
429  //int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
430  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
431 
432  // double startXCOR = bINJ>bREC ? bINJ : bREC;
433  // double endXCOR = eINJ<eREC ? eINJ : eREC;
434  int sizeXCOR = int((eREC-bREC)*R+0.5);
435 
436  // length of temporary buffer for tf plots
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();
440  //wts_size=wfREC->size();
441  cout<<"SIXX "<<sizeXCOR<<" "<<wts_size<<endl;
442 
443  char fname[1024];
444  Meyer<double> S(1024,2); // set wavelet for production
445  double wts_start,wts_stop;
446 
447 
448  // cout<<wfREC->start()-EVT->gps[0]<<" bbb "<<bREC<< " "<<double(bREC+wts_size/2)<<" "<<wfREC->rate()<<endl;
449 
450  // ------------------------------------------------------------------------
451  // RECONSTRUCTED WAVEFORM
452  // ------------------------------------------------------------------------
453 
454  wavearray<double> xREC(wts_size);
455  xREC.start(wfREC->start()-EVT->gps[0]);
456  // xREC.start(wfREC->start()-EVT->gps[0]+double(bREC+wts_size/2)/wfREC->rate());
457  xREC.rate(wfREC->rate());
458  xREC=0.;
459  for (int m=0;m<sizeXCOR;m++)
460  if(m<(int)xREC.size()/2) xREC.data[m+wts_size/2]=wfREC->data[m];
461  WSeries<double> wREC(xREC,S);
462  xREC.resize(0);
463 
464  if(NET->getwdm(optLayer+1)) wREC.Forward(wREC,*NET->getwdm(optLayer+1));
465  cout<<"LLLLL "<< wREC.start()+(double)(wts_size/2)/wREC.rate() <<" "<<wts_start+sizeXCOR/wREC.rate()<<endl;
466 
467  entf_REC = ComputeEnergy(&wREC);
468  cout << "TF : " << "enINJ : "<< " enREC : " << entf_REC
469  << endl;
470  double tm, fq,stm,sfq, bd, dr;
471  ComputeTFtimefreq(&wREC,tm,stm,fq,sfq,dr,bd);
472  TFtime_REC= tm;
473  sig_TFtime_REC =stm;
474  TFfreq_REC =fq;
475  sig_TFfreq_REC=sfq;
476 
477  TFduration=dr;
478  TFband =bd;
479 
480  cout<<"TEST freqtime"<<TFfreq_REC<<" "<< sig_TFfreq_REC<<" TFtime " <<TFtime_REC<<" band "<<TFband<< endl;
481 
482  return;
483  }
484 
485 
486 
487 
489 
490  int layers = WS->maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
491  int slices = WS->sizeZero(); // number of time bins
492 
493  float df = WS->resolution(); // frequency bin resolution (hz)
494  float dt = 1./(2*df); // time bin resolution (sec)
495 
496  int rate = int(1./dt);
497 
498 // cout << "layers : " << layers << "\t slices : " << slices << "\t rate : " << rate
499 // << "\t dt : " << dt << "\t df : " << df << endl;
500 
501  double EE=0;
502  double pixel_frequency;
503  double pixel_time;
504  for(int i=1;i<=slices;i++) {
505  pixel_time = i*dt;
506  for(int j=1;j<=layers;j++) {
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;
510  double ER = pow(WS->getSample(i-1,j-1),2);
511  if(j==1) EE += ER/2.;
512  if(j==layers) EE += ER/2.;
513  if((j!=1)&&(j!=layers)) EE += ER;
514  //cout << pixel_time << " : " << pixel_frequency << endl;
515  }
516  }
517 
518  return EE;
519 }
520 
521 
522 
523 void ComputeTFtimefreq(WSeries<double> *WS, double& TFtime, double& sig_TFtime, double& TFfreq, double& sig_TFfreq, double& TFduration, double& TFband) {
524 
525  int layers = WS->maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
526  int slices = WS->sizeZero(); // number of time bins
527 
528  float df = WS->resolution(); // frequency bin resolution (hz)
529  float dt = 1./(2*df); // time bin resolution (sec)
530 
531  int rate = int(1./dt);
532 
533  // cout << "layers : " << layers << "\t slices : " << slices << "\t rate : " << r\
534 ate
535  // << "\t dt : " << dt << "\t df : " << df << endl;
536  double timeEn =0; double freqEn=0;
537  double EE=0; double EE4=0;
538  double pixel_frequency;
539  double pixel_time;
540 
541  for(int i=1;i<=slices;i++) {
542  pixel_time = i*dt;
543  for(int j=1;j<=layers;j++) {
544  double f = 0;
545  double t = 0;
546  double en=0;
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;
550 
551  en = pow(WS->getSample(i-1,j-1),2);
552  f = en*pixel_frequency;
553  t = en*pixel_time;
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;}
557  // if(en !=0) cout<<"AAIUTO "<< pixel_time << " : " << pixel_frequency << " "<<f<<" "<<t<<" "<<en<<" Sig "<<EE4<<" "<<endl;
558  }
559  }
560 
561  // cout<<"??? "<<timeEn<<" "<<EE<<" "<<freqEn<<" "<<EE4<<endl;
562  TFtime = timeEn/EE;
563  sig_TFtime = dt*sqrt(1./12.)*sqrt(EE4)/EE;
564  TFfreq = freqEn/EE;
565  sig_TFfreq= df*sqrt(1./12.)*sqrt(EE4)/EE;
566 
567 double bandEn =0; double durationEn=0;
568 
569  for(int i=1;i<=slices;i++) {
570  pixel_time = i*dt;
571  for(int j=1;j<=layers;j++) {
572  double band = 0;
573  double duration = 0;
574  double en=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;
578 
579  en = pow(WS->getSample(i-1,j-1),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;}
585  // if(en !=0) cout<<"AAIUTO "<< pixel_time << " : " << pixel_frequency << " "<<f<<" "<<t<<" "<<en<<" Sig "<<EE4<<" "<<endl;
586  }
587  }
588 
589  TFband = sqrt(bandEn/EE);
590  TFduration= sqrt(durationEn/EE);
591  cout<<" BBB DDF "<<TFband <<" "<<TFduration<<endl;
592 
593 
594 
595 
596  return;
597  }
598 
wavearray< double > t(hp.size())
std::vector< char * > ifoName
Definition: network.hh:609
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
CWB::config * cfg
char analysis[8]
Definition: config.hh:117
TString outDump
size_t nLag
Definition: network.hh:573
std::vector< vector_int > cRate
Definition: netcluster.hh:398
int slices
std::vector< wavearray< double > > wREC[MAX_TRIALS]
double duration
TMacro configPlugin
Definition: config.hh:362
TString Get(wavearray< double > &x, TString ifo)
Definition: mdc.cc:1529
bool mdcPlugin
Definition: config.hh:365
std::vector< std::string > mdcList
Definition: mdc.hh:389
virtual void rate(double r)
Definition: wavearray.hh:141
float factor
int n
Definition: cwb_net.C:28
TTree * setTree()
Definition: netevent.cc:434
TString("c")
std::vector< std::string > mdcType
Definition: mdc.hh:390
int ID
Definition: TestMDC.C:70
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...
Definition: netcluster.cc:2207
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:382
netcluster * pwc
Definition: cwb_job_obj.C:38
bool cedDump
Definition: config.hh:297
int layers
std::vector< std::string > mdcType
Definition: network.hh:613
CWB::mdc * MDC
Long_t size
int m
Definition: cwb_net.C:28
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
double rate
Definition: netcluster.hh:378
std::vector< double > mdcTime
Definition: network.hh:614
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
bool outPlugin
Definition: config.hh:369
void Print(int level=0)
Definition: mdc.cc:2736
char ifo[NIFO_MAX][8]
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:312
size_t ifoListSize()
Definition: network.hh:431
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
Definition: network.hh:448
#define nIFO
virtual size_t size() const
Definition: wavearray.hh:145
jfile
Definition: cwb_job_obj.C:43
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:711
TTree * net_tree
Definition: mdc.hh:248
int simulation
Definition: config.hh:199
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:69
Double_t * gps
average center_of_gravity time
Definition: netevent.hh:95
i() int(T_cor *100))
network NET
Definition: cwb_dump_inj.C:30
std::vector< std::string > mdcList
Definition: network.hh:612
TIter next(twave->GetListOfBranches())
char fname[1024]
std::vector< int > RWFID
Definition: detector.hh:381
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
int k
double ComputeEnergy(WSeries< double > *WS)
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:117
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)
TFile * froot
double dt
WSeries< double > TFmap
Definition: detector.hh:354
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:87
bool wfsave
Definition: network.hh:600
netevent EVT(itree, nifo)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:88
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:94
std::vector< int > sCuts
Definition: netcluster.hh:392
void dclose()
Definition: netevent.hh:321
Definition: Meyer.hh:36
virtual void stop(double s)
Definition: wavearray.hh:139
double resolution(int=0)
Definition: wseries.hh:155
int gIFACTOR
char cmd[1024]
Meyer< double > S(1024, 2)
double df
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
DataType_t * data
Definition: wavearray.hh:319
Long_t id
double factors[FACTORS_MAX]
Definition: config.hh:202
DataType_t getSample(int n, double m)
Definition: wseries.hh:185
std::vector< double > mdcTime
Definition: mdc.hh:392
virtual void resize(unsigned int)
Definition: wavearray.cc:463
int check
size_t sizeZero()
Definition: wseries.hh:144
int maxLayer()
Definition: wseries.hh:139
bool dump
Definition: config.hh:295