Logo coherent WaveBurst  
Library Reference Guide
Logo
Read_nRMS.C
Go to the documentation of this file.
1 // This macro read the nRMS from the cstrain root files produced with
2 // the plugin CWB_Plugin_nRMS.C
3 // Author : Gabriele Vedovato
4 
5 //#define SAVE_PLOT // save plots
6 
8 void ReadFile(TString fName);
9 
10 void Read_nRMS() {
11 
13 
14  cout << "Starting reading output directory ..." << endl;
15  vector<TString> fileList = TB.getFileListFromDir("output",".root","cstrain_","",true);
16  int nfile = fileList.size();
17  for(int n=0;n<nfile;n++) {
18  ReadFile(fileList[n]);
19  }
20 
21  exit(0);
22 }
23 
25 
26  network NET;
27  CWB::config CFG;
28  detector* pD[NIFO_MAX]; //! pointers to detectors
29 
30  // open root file
31  TFile* ifile = new TFile(fName);
32  if(ifile==NULL) {cout << "Error opening root file : " << fName.Data() << endl;exit(1);}
33  //ifile->ls();
34 
35  // read network object from job file
36  network* pnet = (network*)ifile->Get("network");
37  if(pnet!=NULL) {
38  // copy network object into local NET object
39  NET = *pnet;
40  delete pnet;
41  } else {
42  cout << "cwb::InitNetwork - Error : net is not contained in root file " << fName.Data() << endl;
43  exit(1);
44  }
45 
46  // read config object
47  if(ifile->Get("config")!=NULL) {
48  // read config object
49  CWB::config* pcfg = (CWB::config*)ifile->Get("config");
50  CFG = *pcfg;
51  delete pcfg;
52  }
53 
54  int nIFO = NET.ifoListSize();
55 
56  // save detector pointers into local structures pD
57  for(int n=0; n<nIFO; n++) pD[n] = NET.getifo(n);
58 
59  char cdrms_name[32];sprintf(cdrms_name,"rms-f%d",0);
60 
61  // read strain rms
62  for(int n=0; n<nIFO; n++) {
63  //cout << n << " " << pD[n]->Name << endl;
64 
65  ifile->cd();
66  WSeries<double>* pws = (WSeries<double>*)ifile->Get(TString("rms/")+cdrms_name+"/"+pD[n]->Name);
67  if(pws==NULL)
68  {cout << "Error : strain rms not present, job terminated!!!" << endl;exit(1);}
69  pD[n]->nRMS=*pws;;
70  delete pws;
71  }
72 
73  ifile->Close();
74 
75  // dump to file
76  for(int n=0; n<nIFO; n++) {
77  DumpToFile(pD[n]->nRMS, -1., &CFG, &NET, pD[n]->Name, fName);
78  }
79 
80 }
81 
83 
84  // gps<0 : the rms is the the average over the all rms samples in the buffer
85  // gps=0 : the selected rms is in the middle of the buffer
86  // gps>0 : only the rms at that time is selected
87 
88  double GPS = gps;
89 
90  bool oneside = true;
91 
92  int levels = nRMS.getLevel()+1; // number of levels
93  int slices = nRMS.size()/levels; // number of nRMS samples
94  double length = slices*cfg->whiteStride; // nRMS len in sec
95 
96  // nRMS do not come from a TF transform
97  // nRMS params must be fixed before to pass to watplot
98  double rate = nRMS.rate();
99  double fNinquist = rate/2.;
100  nRMS.rate(nRMS.size()/length);
101  WDM<double>* wdm = (WDM<double>*) nRMS.pWavelet;
102  wdm->nSTS=nRMS.size();
103 
104  // Plot nRMS Scalogram
105  double start = nRMS.start();
106  double stop = nRMS.start()+nRMS.size()/nRMS.rate();
107  double flow = cfg->fLow;
108  double fhigh = cfg->fHigh;
109 
110  int M = nRMS.getLevel();
111  const double scale = 1./nRMS.wrate();
112 
113  int igps;
114  if(GPS<=0) {
115  igps = (nRMS.stop()-nRMS.start())/scale/2; // time index half buffer
116  gps = nRMS.start()+igps*scale;
117  } else {
118  igps = (gps-nRMS.start())/scale + 1; // time index at gps time
119  }
120 
121  int N = nRMS.size()/(M+1); // number of time indexes
122 /*
123  // dump nRMS samples at the same frequency
124  for(int j=0; j<N; j++) {
125  cout << j*scale << " " << nRMS.data[j*(M+1)+100] << endl;
126  }
127 */
128 
129  // dump nRMS to file
130  fName.ReplaceAll("cstrain_","nrms_");
131  fName.ReplaceAll(".root",".txt");
132  cout << endl << "Dump nRMS " << " : " << fName << endl;
133 
134  ofstream out;
135  out.open(fName.Data(),ios::out);
136  if(!out.good()) {cout << "Error Opening File : " << fName << endl;exit(1);}
137 
138  double inRate = cfg->fResample>0 ? cfg->fResample : cfg->inRate;
139  double rescale = oneside ? sqrt(2.)/sqrt(inRate) : 1./sqrt(inRate);
140 
141  double df=fNinquist/M;
142  for(int j=1; j<M; ++j) {
143  double f = j*df+df/2.;
144  if(f<flow || f>fhigh) continue;
145  double rms = 0;
146  if(GPS<0) {
147  for(int i=0; i<N; i++) rms += nRMS.data[i*(M+1)+j]*rescale;
148  rms/=N;
149  } else {
150  rms = nRMS.data[igps*(M+1)+j]*rescale;
151  }
152  out << j*df+df/2. << " " << rms << endl;
153  }
154 
155  out.close();
156 
157 #ifdef SAVE_PLOT
158  char cmd[1024];
159  // execute cwb_draw_sensitivity
160  sprintf(cmd,"export CWB_SENSITIVITY_FILE_NAME=\"%s\";",fName.Data());
161  sprintf(cmd,"%s export CWB_SENSITIVITY_SAVE_PLOT=%d;",cmd,1);
162  sprintf(cmd,"%s export CWB_SENSITIVITY_RANGE_FIX=%d;",cmd,1);
163  sprintf(cmd,"%s root -n -l -b ${CWB_MACROS}/cwb_draw_sensitivity.C",cmd);
164  cout << cmd << endl; gSystem->Exec(cmd);
165 #endif
166 
167 }
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
CWB::config * cfg
int slices
static vector< TString > getFileListFromDir(TString dir_name, TString endString="", TString beginString="", TString containString="", bool fast=false)
Definition: Toolbox.cc:5108
double M
Definition: DrawEBHH.C:13
virtual void rate(double r)
Definition: wavearray.hh:141
int n
Definition: cwb_net.C:28
TString("c")
ofstream out
Definition: cwb_merge.C:214
double fLow
Definition: config.hh:140
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
CWB::Toolbox TB
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
#define N
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
char ifo[NIFO_MAX][8]
network ** net
NOISE_MDC_SIMULATION.
size_t ifoListSize()
Definition: network.hh:431
#define nIFO
void wrate(double r)
Definition: wseries.hh:120
virtual size_t size() const
Definition: wavearray.hh:145
int getLevel()
Definition: wseries.hh:109
size_t fResample
Definition: config.hh:142
unsigned long nSTS
Definition: WaveDWT.hh:143
void ReadFile(TString fName)
Definition: Read_nRMS.C:24
WSeries< double > nRMS
Definition: detector.hh:358
network NET
Definition: cwb_dump_inj.C:30
const int NIFO_MAX
Definition: wat.hh:22
double fhigh
TFile * ifile
double flow
double gps
double whiteStride
Definition: config.hh:190
double fHigh
Definition: config.hh:141
cout<< "Starting reading output directory ..."<< endl;vector< TString > fileList
virtual void stop(double s)
Definition: wavearray.hh:139
#define GPS
char cmd[1024]
double df
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
DataType_t * data
Definition: wavearray.hh:319
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:456
void Read_nRMS()
Definition: Read_nRMS.C:10
char fName[256]
double length
Definition: TestBandPass.C:18
size_t inRate
Definition: config.hh:132
detector ** pD
void DumpToFile(WSeries< double > &nRMS, double gps, CWB::config *cfg, network *net, TString ifo, TString fName)
Definition: Read_nRMS.C:82
exit(0)