Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CWB_Plugin_SimMDC_SimData.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 "gnetwork.hh"
16 
17 #define WAVEFORM_NAME "Waveforms/SG554Q8d9.txt"
18 
19 
20 void
21 CWB_Plugin(TFile* jfile, CWB::config* cfg, network* net, WSeries<double>* x, TString ifo, int type) {
22 
23  cout << endl;
24  cout << "-----> plugins/CWB_Plugin_SimMDC_SimData.C" << endl;
25  cout << "ifo " << ifo.Data() << endl;
26  cout << "type " << type << endl;
27  cout << endl;
28 
29 
30  if(type==CWB_PLUGIN_CONFIG) {
31  cfg->dataPlugin=true; // disable read data from frames
32  cfg->mdcPlugin=true; // disable read mdc from frames
33  }
34 
35  if(type==CWB_PLUGIN_DATA) {
36 
37  CWB::Toolbox TB;
38 
39  int seed;
40  if(ifo.CompareTo("L1")==0) seed=1000;
41  if(ifo.CompareTo("H1")==0) seed=2000;
42  if(ifo.CompareTo("V1")==0) seed=3000;
43  if(ifo.CompareTo("J1")==0) seed=4000;
44  if(ifo.CompareTo("A2")==0) seed=5000;
45 
46  TString fName;
47  if(ifo.CompareTo("L1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
48  if(ifo.CompareTo("H1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
49  if(ifo.CompareTo("V1")==0) fName="plugins/strains/advVIRGO_sensitivity_12May09_8khz_one_side.txt";
50  if(ifo.CompareTo("J1")==0) fName="plugins/strains/LCGT_sensitivity_8khz_one_side.txt";
51  if(ifo.CompareTo("A2")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
52 
53  int size=x->size();
54  double start=x->start();
55  TB.getSimNoise(*x, fName, seed, net->nRun);
56  x->resize(size);
57  x->start(start);
58 /*
59  // dump spectrum
60  char file[1024];
61  sprintf(file,"%s/sensitivity_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
62  cout << endl << "Dump Sensitivity : " << file << endl << endl;
63  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
64  if(TString(ifo).CompareTo("V1")==0) gSystem->Exit(0);
65 */
66  }
67 
68  if(type==CWB_PLUGIN_MDC) {
69 
70  double dt = 1./x->rate();
71 
72  // log burstMDC parameters
73  double burstMDC_theta = 7.799320e-01;
74  double burstMDC_phi = 7.720941e-01;
75  double burstMDC_psi = 3.397006;
76 
77  /* ------------------------------
78  earthCenterTime = _EarthCtrGPS;
79  phi = _External_phi;
80  theta = _External_x;
81  psi = _External_psi;
82  --------------------------------- */
83  double theta,phi,psi;
84  double Pi = TMath::Pi();
85  if(true) {
86  theta = acos(burstMDC_theta);
87  theta*= 180/Pi;
88  phi = burstMDC_phi > 0 ? burstMDC_phi : 2*Pi+burstMDC_phi;
89  phi*= 180/Pi;
90 // phi = sm.phi2RA(phi,_EarthCtrGPS);
91  psi = burstMDC_psi*180/Pi;
92  }
93  cout << "theta : " << theta << " phi : " << phi << " psi " << psi << endl;
94 
95  int nIFO=net->ifoListSize();
96  TString ifos[NIFO_MAX];
97  for(int n=0;n<nIFO;n++)ifos[n]= net->ifoName[n];
98  gnetwork gNET(nIFO,ifos);
99  for(int i=0;i<nIFO;i++) {
100  for(int j=i+1;j<nIFO;j++) {
101  cout << ifos[i].Data() << " " << ifos[j].Data() << " -> "
102  << gNET.GetDelay(ifos[i].Data(),ifos[j].Data(),phi,theta) << " sec " << endl;
103  }
104  }
105 
106 
107  char logString[1024]="";
108 
109  char GravEn_SimID[1024]=WAVEFORM_NAME;
110  double SimHrss = 1.213000e-21;
111  double SimEgwR2 = 1.011301e-47;
112  double GravEn_Ampl = 1.213000e-21;
113  double Internal_x = 1.0;
114  double Internal_phi = 0.0;
115  double External_x = burstMDC_theta;
116  double External_phi = burstMDC_phi;
117  double External_psi = burstMDC_psi;
118 
119  double FrameGPS = x->start();
120  //double EarthCtrGPS = FrameGPS+dt*x->size()/2.; // injected in the center of buffer
121  double EarthCtrGPS = 968654066.616913;
122  char SimName[64] = "SG554Q8d9";
123  double SimHpHp = 1.471369e-42;
124  double SimHcHc = 0;
125  double SimHpHc = 0;
126 
127  sprintf(logString,"%s",GravEn_SimID);
128  sprintf(logString,"%s %e",logString,SimHrss);
129  sprintf(logString,"%s %e",logString,SimEgwR2);
130  sprintf(logString,"%s %e",logString,GravEn_Ampl);
131  sprintf(logString,"%s %e",logString,Internal_x);
132  sprintf(logString,"%s %e",logString,Internal_phi);
133  sprintf(logString,"%s %e",logString,External_x);
134  sprintf(logString,"%s %e",logString,External_phi);
135  sprintf(logString,"%s %e",logString,External_psi);
136 
137  sprintf(logString,"%s %10.6f",logString,FrameGPS);
138  sprintf(logString,"%s %10.6f",logString,EarthCtrGPS);
139  sprintf(logString,"%s %s",logString,SimName);
140  sprintf(logString,"%s %e",logString,SimHpHp);
141  sprintf(logString,"%s %e",logString,SimHcHc);
142  sprintf(logString,"%s %e",logString,SimHpHc);
143 
144  double tShift=0;
145  double fPlus=0;
146  double fCross=0;
147  for(int i=0;i<nIFO;i++) {
148  double IFOctrGPS = EarthCtrGPS;
149  if(i>0) IFOctrGPS += gNET.GetDelay(ifos[i].Data(),ifos[0].Data(),phi,theta);
150  double IFOfPlus = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, true);
151  double IFOfCross = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, false);
152  if(ifos[i].CompareTo(ifo)==0) {
153  tShift=IFOctrGPS-EarthCtrGPS;
154  fPlus=IFOfPlus;
155  fCross=IFOfCross;
156  }
157  sprintf(logString,"%s %s %10.6f %e %e",logString,ifos[i].Data(),IFOctrGPS,IFOfPlus,IFOfCross);
158  }
159 
160  if(ifo.CompareTo(net->ifoName[0])==0) {
161  cout << logString << endl;
162  net->mdcList.clear();
163  net->mdcList.push_back(logString);
164  net->mdcType.clear();
165  net->mdcType.push_back(SimName);
166  net->mdcTime.clear();
167  net->mdcTime.push_back(EarthCtrGPS);
168  }
169 /*
170  cout << endl << endl;
171  for (int nmdc=0; nmdc<(int)net->mdcListSize(); nmdc++) {
172  TString mdcstring(net->getmdcList(nmdc));
173 // if(nmdc==72) {
174  cout << endl << endl;
175  cout << "--------------------------------> " << nmdc << endl;
176  cout << endl << endl;
177  cout << mdcstring.Data() << endl;
178 // }
179  }
180  //gSystem->Exit(0);
181 */
182  // read waveform
183  ifstream in;
184  in.open(WAVEFORM_NAME, ios::in);
185  if (!in.good()) {cout << "Error Opening File : " << WAVEFORM_NAME << endl;gSystem->Exit(1);}
186 
187  // fill mdc vector
188  (*x)=0;
189  int offset = (EarthCtrGPS-FrameGPS)*x->rate();
190  int size=0;
191  double hrss=0;
192  while (1) {
193  in >> x->data[offset+size];
194  hrss+=pow(x->data[offset+size],2);
195  if (!in.good()) break;
196  size++;
197  if ((offset+size)>=(int)x->size()) break;
198  }
199  hrss=sqrt(hrss*dt);
200  for(int i=0;i<(int)x->size();i++) x->data[i] = x->data[i]*fPlus*SimHrss/hrss;
201 
202  in.close();
203 
204  if(tShift==0) return;
205 
206  // apply time shift to mdc vector
207  x->FFTW(1);
208  TComplex C;
209  double df = x->rate()/x->size();
210  //double tShift = 10./x->rate();
211  cout << "tShift : " << tShift << endl;
212  for (int ii=0;ii<(int)x->size()/2;ii++) {
213  TComplex X(x->data[2*ii],x->data[2*ii+1]);
214  X=X*C.Exp(TComplex(0.,-2*Pi*ii*df*tShift)); // Time Shift
215  x->data[2*ii]=X.Re();
216  x->data[2*ii+1]=X.Im();
217  }
218  x->FFTW(-1);
219 
220  }
221 
222  if(type==CWB_PLUGIN_WHITE) {
223 /*
224  CWB::Toolbox TB;
225  int level=x->getLevel();
226  x->Inverse(-1);
227  // dump spectrum
228  char file[1024];
229  sprintf(file,"%s/sensitivity_white_%s_%d_%s_job%lu.txt",cfg->dump_dir,ifo.Data(),int(x->start()),cfg->data_label,net->nRun);
230  cout << endl << "Dump Sensitivity : " << file << endl << endl;
231  TB.makeSpectrum(file, *x, 8, cfg->segEdge);
232  if(TString(ifo).CompareTo("V1")==0) gSystem->Exit(0);
233  x->Forward(level);
234 */
235  }
236 
237  return;
238 }
network::ifoName
std::vector< char * > ifoName
Definition: network.hh:591
cfg
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
WSeries::resize
virtual void resize(unsigned int)
Definition: wseries.cc:883
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
C
static const double C
Definition: GNGen.cc:10
gNET
gnetwork * gNET
Definition: cwb_draw_antpat.C:47
CWB::config::mdcPlugin
bool mdcPlugin
Definition: config.hh:347
offset
int offset
Definition: TestSTFT_2.C:19
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:123
CWB::config::dataPlugin
bool dataPlugin
Definition: config.hh:346
n
int n
Definition: cwb_net.C:10
TString
TString("c")
Definition: cwb_report_skymap.C:111
network::nRun
size_t nRun
Definition: network.hh:554
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<< "\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
Definition: BaseFunctionWDM.C:51
theta
float theta
Definition: cwb_condor_create_ced.C:51
network::mdcType
std::vector< std::string > mdcType
Definition: network.hh:595
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:21
size
Long_t size
Definition: cwb_condor_check.C:48
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:119
j
int j
Definition: cwb_net.C:10
i
i drho i
Definition: cwb_epparameters.C:85
network::mdcTime
std::vector< double > mdcTime
Definition: network.hh:596
TB
CWB::Toolbox TB
Definition: ComputeSNR.C:5
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:413
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
phi
float phi
Definition: cwb_condor_create_ced.C:51
hrss
double hrss
Definition: TestMDC.C:70
psi
float psi
Definition: cwb_condor_create_ced.C:51
jfile
jfile
Definition: cwb_job_obj.C:25
CWB::Toolbox::getSimNoise
static void getSimNoise(wavearray< double > &u, TString fName, int seed, int run)
Definition: Toolbox.cc:4809
seed
int seed
Definition: CWB_Plugin_HEN_BKG_Config.C:13
int
i() int(T_cor *100))
Pi
double Pi
Definition: DrawPhaseShift.C:12
network::mdcList
std::vector< std::string > mdcList
Definition: network.hh:594
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:4
WAVEFORM_NAME
#define WAVEFORM_NAME
Definition: CWB_Plugin_SimMDC_SimData.C:17
dt
double dt
Definition: DrawPhaseShift.C:13
wavearray::FFTW
virtual void FFTW(int=1)
Definition: wavearray.cc:878
gnetwork::GetDelay
double GetDelay(TString ifo1, TString ifo2, double phi, double theta)
Definition: gnetwork.cc:907
in
ifstream in
Definition: CWB_Plugin_HEN_BKG_Config.C:20
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:545
wavearray::data
DataType_t * data
Definition: wavearray.hh:301
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:52