Logo coherent WaveBurst  
Library Reference Guide 		Logo
CWB_Plugin_SimMDC_SimData.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #define XIFO 4
20 
21 #pragma GCC system_header
22 
23 #include "cwb.hh"
24 #include "config.hh"
25 #include "network.hh"
26 #include "wavearray.hh"
27 #include "TString.h"
28 #include "TObjArray.h"
29 #include "TObjString.h"
30 #include "TRandom.h"
31 #include "TComplex.h"
32 #include "TMath.h"
33 #include "gnetwork.hh"
34 
35 #define WAVEFORM_NAME "Waveforms/SG554Q8d9.txt"
36 
37 
38 void
39 CWB_Plugin(TFile* jfile, CWB::config* cfg, network* net, WSeries<double>* x, TString ifo, int type) {
40 
41  cout << endl;
42  cout << "-----> plugins/CWB_Plugin_SimMDC_SimData.C" << endl;
43  cout << "ifo " << ifo.Data() << endl;
44  cout << "type " << type << endl;
45  cout << endl;
46 
47 
48  if(type==CWB_PLUGIN_CONFIG) {
49  cfg->dataPlugin=true; // disable read data from frames
50  cfg->mdcPlugin=true; // disable read mdc from frames
51  }
52 
53  if(type==CWB_PLUGIN_DATA) {
54 
55  CWB::Toolbox TB;
56 
57  int seed;
58  if(ifo.CompareTo("L1")==0) seed=1000;
59  if(ifo.CompareTo("H1")==0) seed=2000;
60  if(ifo.CompareTo("V1")==0) seed=3000;
61  if(ifo.CompareTo("J1")==0) seed=4000;
62  if(ifo.CompareTo("A2")==0) seed=5000;
63 
64  TString fName;
65  if(ifo.CompareTo("L1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
66  if(ifo.CompareTo("H1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
67  if(ifo.CompareTo("V1")==0) fName="plugins/strains/advVIRGO_sensitivity_12May09_8khz_one_side.txt";
68  if(ifo.CompareTo("J1")==0) fName="plugins/strains/LCGT_sensitivity_8khz_one_side.txt";
69  if(ifo.CompareTo("A2")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
70 
71  int size=x->size();
72  double start=x->start();
73  TB.getSimNoise(*x, fName, seed, net->nRun);
74  x->resize(size);
75  x->start(start);
76 /*
77  // dump spectrum
78  char file[1024];
79  sprintf(file,"%s/sensitivity_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
80  cout << endl << "Dump Sensitivity : " << file << endl << endl;
81  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
82  if(TString(ifo).CompareTo("V1")==0) gSystem->Exit(0);
83 */
84  }
85 
86  if(type==CWB_PLUGIN_MDC) {
87 
88  double dt = 1./x->rate();
89 
90  // log burstMDC parameters
91  double burstMDC_theta = 7.799320e-01;
92  double burstMDC_phi = 7.720941e-01;
93  double burstMDC_psi = 3.397006;
94 
95  /* ------------------------------
96  earthCenterTime = _EarthCtrGPS;
97  phi = _External_phi;
98  theta = _External_x;
99  psi = _External_psi;
100  --------------------------------- */
101  double theta,phi,psi;
102  double Pi = TMath::Pi();
103  if(true) {
104  theta = acos(burstMDC_theta);
105  theta*= 180/Pi;
106  phi = burstMDC_phi > 0 ? burstMDC_phi : 2*Pi+burstMDC_phi;
107  phi*= 180/Pi;
108 // phi = sm.phi2RA(phi,_EarthCtrGPS);
109  psi = burstMDC_psi*180/Pi;
110  }
111  cout << "theta : " << theta << " phi : " << phi << " psi " << psi << endl;
112 
113  int nIFO=net->ifoListSize();
114  TString ifos[NIFO_MAX];
115  for(int n=0;n<nIFO;n++)ifos[n]= net->ifoName[n];
116  gnetwork gNET(nIFO,ifos);
117  for(int i=0;i<nIFO;i++) {
118  for(int j=i+1;j<nIFO;j++) {
119  cout << ifos[i].Data() << " " << ifos[j].Data() << " -> "
120  << gNET.GetDelay(ifos[i].Data(),ifos[j].Data(),phi,theta) << " sec " << endl;
121  }
122  }
123 
124 
125  char logString[1024]="";
126 
127  char GravEn_SimID[1024]=WAVEFORM_NAME;
128  double SimHrss = 1.213000e-21;
129  double SimEgwR2 = 1.011301e-47;
130  double GravEn_Ampl = 1.213000e-21;
131  double Internal_x = 1.0;
132  double Internal_phi = 0.0;
133  double External_x = burstMDC_theta;
134  double External_phi = burstMDC_phi;
135  double External_psi = burstMDC_psi;
136 
137  double FrameGPS = x->start();
138  //double EarthCtrGPS = FrameGPS+dt*x->size()/2.; // injected in the center of buffer
139  double EarthCtrGPS = 968654066.616913;
140  char SimName[64] = "SG554Q8d9";
141  double SimHpHp = 1.471369e-42;
142  double SimHcHc = 0;
143  double SimHpHc = 0;
144 
145  sprintf(logString,"%s",GravEn_SimID);
146  sprintf(logString,"%s %e",logString,SimHrss);
147  sprintf(logString,"%s %e",logString,SimEgwR2);
148  sprintf(logString,"%s %e",logString,GravEn_Ampl);
149  sprintf(logString,"%s %e",logString,Internal_x);
150  sprintf(logString,"%s %e",logString,Internal_phi);
151  sprintf(logString,"%s %e",logString,External_x);
152  sprintf(logString,"%s %e",logString,External_phi);
153  sprintf(logString,"%s %e",logString,External_psi);
154 
155  sprintf(logString,"%s %10.6f",logString,FrameGPS);
156  sprintf(logString,"%s %10.6f",logString,EarthCtrGPS);
157  sprintf(logString,"%s %s",logString,SimName);
158  sprintf(logString,"%s %e",logString,SimHpHp);
159  sprintf(logString,"%s %e",logString,SimHcHc);
160  sprintf(logString,"%s %e",logString,SimHpHc);
161 
162  double tShift=0;
163  double fPlus=0;
164  double fCross=0;
165  for(int i=0;i<nIFO;i++) {
166  double IFOctrGPS = EarthCtrGPS;
167  if(i>0) IFOctrGPS += gNET.GetDelay(ifos[i].Data(),ifos[0].Data(),phi,theta);
168  double IFOfPlus = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, true);
169  double IFOfCross = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, false);
170  if(ifos[i].CompareTo(ifo)==0) {
171  tShift=IFOctrGPS-EarthCtrGPS;
172  fPlus=IFOfPlus;
173  fCross=IFOfCross;
174  }
175  sprintf(logString,"%s %s %10.6f %e %e",logString,ifos[i].Data(),IFOctrGPS,IFOfPlus,IFOfCross);
176  }
177 
178  if(ifo.CompareTo(net->ifoName[0])==0) {
179  cout << logString << endl;
180  net->mdcList.clear();
181  net->mdcList.push_back(logString);
182  net->mdcType.clear();
183  net->mdcType.push_back(SimName);
184  net->mdcTime.clear();
185  net->mdcTime.push_back(EarthCtrGPS);
186  }
187 /*
188  cout << endl << endl;
189  for (int nmdc=0; nmdc<(int)net->mdcListSize(); nmdc++) {
190  TString mdcstring(net->getmdcList(nmdc));
191 // if(nmdc==72) {
192  cout << endl << endl;
193  cout << "--------------------------------> " << nmdc << endl;
194  cout << endl << endl;
195  cout << mdcstring.Data() << endl;
196 // }
197  }
198  //gSystem->Exit(0);
199 */
200  // read waveform
201  ifstream in;
202  in.open(WAVEFORM_NAME, ios::in);
203  if (!in.good()) {cout << "Error Opening File : " << WAVEFORM_NAME << endl;gSystem->Exit(1);}
204 
205  // fill mdc vector
206  (*x)=0;
207  int offset = (EarthCtrGPS-FrameGPS)*x->rate();
208  int size=0;
209  double hrss=0;
210  while (1) {
211  in >> x->data[offset+size];
212  hrss+=pow(x->data[offset+size],2);
213  if (!in.good()) break;
214  size++;
215  if ((offset+size)>=(int)x->size()) break;
216  }
217  hrss=sqrt(hrss*dt);
218  for(int i=0;i<(int)x->size();i++) x->data[i] = x->data[i]*fPlus*SimHrss/hrss;
219 
220  in.close();
221 
222  if(tShift==0) return;
223 
224  // apply time shift to mdc vector
225  x->FFTW(1);
226  TComplex C;
227  double df = x->rate()/x->size();
228  //double tShift = 10./x->rate();
229  cout << "tShift : " << tShift << endl;
230  for (int ii=0;ii<(int)x->size()/2;ii++) {
231  TComplex X(x->data[2*ii],x->data[2*ii+1]);
232  X=X*C.Exp(TComplex(0.,-2*Pi*ii*df*tShift)); // Time Shift
233  x->data[2*ii]=X.Re();
234  x->data[2*ii+1]=X.Im();
235  }
236  x->FFTW(-1);
237 
238  }
239 
240  if(type==CWB_PLUGIN_WHITE) {
241 /*
242  CWB::Toolbox TB;
243  int level=x->getLevel();
244  x->Inverse(-1);
245  // dump spectrum
246  char file[1024];
247  sprintf(file,"%s/sensitivity_white_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
248  cout << endl << "Dump Sensitivity : " << file << endl << endl;
249  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
250  if(TString(ifo).CompareTo("V1")==0) gSystem->Exit(0);
251  x->Forward(level);
252 */
253  }
254 
255  return;
256 }
network::ifoName
std::vector< char * > ifoName
Definition: network.hh:609
cfg
CWB::config * cfg
Definition: TestCWB_Plugin.C:23
WSeries::resize
virtual void resize(unsigned int)
Definition: wseries.cc:901
C
static const double C
Definition: GNGen.cc:28
gNET
gnetwork * gNET
Definition: cwb_draw_antpat.C:65
CWB::config::mdcPlugin
bool mdcPlugin
Definition: config.hh:365
offset
int offset
Definition: TestSTFT_2.C:19
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:141
CWB::config::dataPlugin
bool dataPlugin
Definition: config.hh:364
n
int n
Definition: cwb_net.C:28
TString
TString("c")
Definition: cwb_report_skymap.C:129
network::nRun
size_t nRun
Definition: network.hh:572
x
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
Definition: BaseFunctionWDM.C:51
theta
float theta
Definition: cwb_condor_create_ced.C:69
TB
CWB::Toolbox TB
Definition: cwb_report_slags.C:29
network::mdcType
std::vector< std::string > mdcType
Definition: network.hh:613
CWB_Plugin
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
Definition: CWB_Plugin_SimMDC_SimData.C:39
size
Long_t size
Definition: cwb_condor_check.C:66
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:137
j
int j
Definition: cwb_net.C:28
i
i drho i
Definition: cwb_epparameters.C:88
network::mdcTime
std::vector< double > mdcTime
Definition: network.hh:614
ifo
char ifo[NIFO_MAX][8]
Definition: cwb1G_parameters.C:17
net
network ** net
NOISE_MDC_SIMULATION.
Definition: CWB_Plugin_HEN_BKG_Config.C:4
network::ifoListSize
size_t ifoListSize()
Definition: network.hh:431
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:145
phi
float phi
Definition: cwb_condor_create_ced.C:69
hrss
double hrss
Definition: TestMDC.C:70
psi
float psi
Definition: cwb_condor_create_ced.C:69
jfile
jfile
Definition: cwb_job_obj.C:43
CWB::Toolbox::getSimNoise
static void getSimNoise(wavearray< double > &u, TString fName, int seed, int run)
Definition: Toolbox.cc:5593
seed
int seed
Definition: CWB_Plugin_HEN_BKG_Config.C:19
int
i() int(T_cor *100))
Pi
double Pi
Definition: DrawPhaseShift.C:12
network::mdcList
std::vector< std::string > mdcList
Definition: network.hh:612
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:22
WAVEFORM_NAME
#define WAVEFORM_NAME
Definition: CWB_Plugin_SimMDC_SimData.C:35
dt
double dt
Definition: DrawPhaseShift.C:13
wavearray::FFTW
virtual void FFTW(int=1)
Definition: wavearray.cc:896
gnetwork::GetDelay
double GetDelay(TString ifo1, TString ifo2, double phi, double theta)
Definition: gnetwork.cc:926
in
ifstream in
Definition: CWB_Plugin_HEN_BKG_Config.C:26
df
double df
Definition: HowToAccessToWSeries.C:36
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
gnetwork::GetAntennaPattern
double GetAntennaPattern(double phi, double theta, double psi=0., int polarization=1)
Definition: gnetwork.cc:564
wavearray::data
DataType_t * data
Definition: wavearray.hh:319
fName
char fName[256]
Definition: CreateListEBBH.C:25
GravEn_SimID
char GravEn_SimID[1024]
Definition: CWB_Plugin_SimMDC_SimData.h:16
ifos
TString ifos[60]
Definition: cwb_draw_antpat.C:70