Logo Coherent WaveBurst  
Reference Guide
Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CWB_Plugin_2G_DataConditioning.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 "regression.hh"
10 #include "TString.h"
11 #include "TObjArray.h"
12 #include "TObjString.h"
13 #include "TRandom.h"
14 #include "Toolbox.hh"
15 
16 void
18 //!DATA_CONDITIONING
19 // Plugin which implement the 2G data conditioning with regression
20 
21  cout << endl;
22  cout << "-----> CWB_Plugin_2G_DataConditioning.C" << endl;
23  cout << "ifo " << ifo.Data() << endl;
24  cout << "type " << type << endl;
25  cout << endl;
26 
27  if(type==CWB_PLUGIN_CONFIG) {
28  cfg->dcPlugin=true; // disable built-in dc
29  }
30 
31  if(type==CWB_PLUGIN_DATA) {
32 /*
33  CWB::Toolbox TB;
34  char file[1024];
35  sprintf(file,"%s/sensitivity_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
36  //sprintf(file,"%s/sensitivity_%s_%d_%s_job%lu.txt",".",ifo.Data(),int(x->start()),"prova",net->nRun);
37  cout << endl << "Dump Sensitivity : " << file << endl << endl;
38  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
39 
40  int nIFO=net->ifoListSize();
41  if(TString(ifo).CompareTo(net->ifoName[nIFO-1])==0) gSystem->Exit(0); // last ifo
42 */
43  }
44 
45  if(type==CWB_PLUGIN_REGRESSION) { // 2G like Data Conditioning
46 
47  if(TString(cfg->analysis)=="1G") {
48  cout << "CWB_Plugin_2G_DataConditioning.C : Error - this plugin works only with 2G" << endl;
49  gSystem->Exit(1);
50  }
51 
52  // get ifo id
53  int id=-1;
54  for(int n=0;n<cfg->nIFO;n++) if(ifo==net->ifoName[n]) {id=n;break;}
55  if(id<0) {cout << "Plugin : Error - bad ifo id" << endl; gSystem->Exit(1);}
56 
57  detector* pD = net->getifo(id);
58  WSeries<double>* pTF = pD->getTFmap();
59 
60  size_t rateANA=cfg->inRate>>cfg->levelR;
61 
62  WDM<double> WDMlpr(int(rateANA/4),int(rateANA/4),6,10);
63 
65  pTF = pD->getTFmap();
66  pTF->Forward(*hot,WDMlpr);
67  regression rr(*pTF,"target",1.,cfg->fHigh);
68  rr.add(*hot,"target");
69  rr.setFilter(8);
70  rr.setMatrix(net->Edge,0.95);
71  rr.solve(0.,4,'h');
72  rr.apply(0.35);
73  *hot = rr.getClean();
74  pTF->Forward(*hot,WDMlpr);
75  pTF->setlow(cfg->fLow);
76  pTF->sethigh(cfg->fHigh);
77  pD->white(cfg->whiteWindow,0,cfg->segEdge,cfg->whiteStride); // calculate noise rms
78  pTF->white(pD->nRMS,1); // whiten WSeries
79  pD->bandPass(16.); // high pass filtering at 16Hz
80  pTF->Inverse();
81  *hot = *pTF; // conditioned TS
82  }
83 
84  if(type==CWB_PLUGIN_WHITE) {
85 
87  int level=-1;
88  if(TString(cfg->analysis)=="1G") {
89  level=x->getLevel();
90  x->Inverse(-1);
91  }
92  // dump spectrum
93  char file[1024];
94  sprintf(file,"%s/sensitivity_white_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
95  cout << endl << "Dump Sensitivity : " << file << endl << endl;
96  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
97  int nIFO=net->ifoListSize();
98  if(TString(ifo).CompareTo(net->ifoName[nIFO-1])==0) gSystem->Exit(0); // last ifo
99  if(TString(cfg->analysis)=="1G") x->Forward(level);
100 
101  }
102 
103  return;
104 }
std::vector< char * > ifoName
Definition: network.hh:591
void sethigh(double f)
Definition: wseries.hh:114
detector * getifo(size_t n)
param: detector index
Definition: network.hh:418
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
char analysis[8]
Definition: config.hh:99
void white(double dT=0, int wtype=1, double offset=0., double stride=0.)
what it does: see algorithm description in wseries.hh
Definition: detector.hh:240
void setMatrix(double edge=0., double f=1.)
Definition: regression.cc:407
int n
Definition: cwb_net.C:10
TString("c")
size_t nRun
Definition: network.hh:554
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
double whiteWindow
Definition: config.hh:171
double fLow
Definition: config.hh:122
bool dcPlugin
Definition: config.hh:348
void setlow(double f)
Definition: wseries.hh:107
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
Definition: regression.cc:73
virtual void start(double s)
Definition: wavearray.hh:119
CWB::Toolbox TB
Definition: ComputeSNR.C:5
char ifo[NIFO_MAX][8]
double Edge
Definition: network.hh:560
network ** net
NOISE_MDC_SIMULATION.
size_t ifoListSize()
Definition: network.hh:413
double segEdge
Definition: config.hh:146
#define nIFO
int getLevel()
Definition: wseries.hh:91
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
jfile
Definition: cwb_job_obj.C:25
wavearray< double > hot[2]
WSeries< double > nRMS
Definition: detector.hh:340
i() int(T_cor *100))
void apply(double threshold=0., char c='a')
Definition: regression.cc:691
WSeries< double > pTF[nRES]
Definition: revMonster.cc:8
void solve(double th, int nE=0, char c='s')
Definition: regression.cc:592
size_t setFilter(size_t)
Definition: regression.cc:258
void bandPass(double f1, double f2, double a=0.)
Definition: detector.hh:265
wavearray< double > getClean()
Definition: regression.hh:117
size_t rateANA
Definition: test_config1.C:21
double whiteStride
Definition: config.hh:172
double fHigh
Definition: config.hh:123
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:161
string file
Definition: cwb_online.py:385
char dump_dir[1024]
Definition: config.hh:310
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:228
static void makeSpectrum(wavearray< double > &psd, wavearray< double > x, double chuncklen=8, double scratchlen=0, bool oneside=true)
Definition: Toolbox.cc:4705
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
int nIFO
Definition: config.hh:102
char data_label[1024]
Definition: config.hh:314
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:273
size_t inRate
Definition: config.hh:114
int levelR
Definition: config.hh:134
detector ** pD
virtual WSeries< double > white(double, int, double=0., double=0.)
what it does: each wavelet layer is devided into k intervals.
Definition: wseries.cc:1128