Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CWB_Plugin_InjectMDC.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 void
20 CWB_Plugin(TFile* jfile, CWB::config* cfg, network* net, WSeries<double>* x, TString ifo, int type) {
21 
22 // Plugin to inject 'on the fly' the MDC
23 
24 /*
25 # GravEn_SimID SimHrss SimEgwR2 GravEn_Ampl Internal_x Internal_phi External_x Externa
26 l_phi External_psi FrameGPS EarthCtrGPS SimName SimHpHp SimHcHc SimHpHc
27 GEO GEOctrG
28 PS GEOfPlus GEOfCross H1 H1ctrGPS H1fPlus H1fCross H2 H2ctrGPS H2fPlus
29 H2fCross L1 L1ctrGPS L1fPlus L1fCross V1 V1ctrGPS V1fPlus V1fCross
30 
31 Waveforms/SG554Q8d9.txt 1.213000e-21 1.011301e-47 1.213000e-21 +1.000000e+00 +0.000000e+00 +7.799320e-01 +7.720941e-01 +3.397006e+00 968653615 968654066.616913 SG554Q8d9 1.471369e-42 0.000000e+00 0.000000e+00 GEO 968654066.597115 +8.635698e-01 -3.501556e-01 H1 968654066.613762 +1.334841e-01 +7.593384e-02 H2 968654066.613762 +1.334841e-01 +7.593384e-02 L1 968654066.616641 -1.870505e-02 -1.359386e-02 V1 968654066.597488 -4.561425e-01 -7.932736e-01
32 */
33 
34  if(type==CWB_PLUGIN_CONFIG) {
35  cfg->mdcPlugin=true; // disable read mdc from frames
36  }
37 
38  if(type==CWB_PLUGIN_MDC) { // mdc
39 
40  double dt = 1./x->rate();
41 
42  cout << endl;
43  cout << "-----> plugins/CWB_Plugin_InjectMDC.C" << endl;
44  cout << "ifo " << ifo.Data() << endl;
45  cout << "type " << type << endl;
46  cout << endl;
47 
48  // log burstMDC parameters
49  double burstMDC_theta = 7.799320e-01;
50  double burstMDC_phi = 7.720941e-01;
51  double burstMDC_psi = 3.397006;
52 
53  /* ------------------------------
54  earthCenterTime = _EarthCtrGPS;
55  phi = _External_phi;
56  theta = _External_x;
57  psi = _External_psi;
58  --------------------------------- */
59  double theta,phi,psi;
60  double Pi = TMath::Pi();
61  if(true) {
62  theta = acos(burstMDC_theta);
63  theta*= 180/Pi;
64  phi = burstMDC_phi > 0 ? burstMDC_phi : 2*Pi+burstMDC_phi;
65  phi*= 180/Pi;
66 // phi = sm.phi2RA(phi,_EarthCtrGPS);
67  psi = burstMDC_psi*180/Pi;
68  }
69  cout << "theta : " << theta << " phi : " << phi << " psi " << psi << endl;
70 
71  int nIFO=net->ifoListSize();
72  TString ifos[NIFO_MAX];
73  for(int n=0;n<nIFO;n++)ifos[n]= net->ifoName[n];
74  gnetwork gNET(nIFO,ifos);
75  for(int i=0;i<nIFO;i++) {
76  for(int j=i+1;j<nIFO;j++) {
77  cout << ifos[i].Data() << " " << ifos[j].Data() << " -> "
78  << gNET.GetDelay(ifos[i].Data(),ifos[j].Data(),phi,theta) << " sec " << endl;
79  }
80  }
81 
82 
83  char logString[1024]="";
84 
85  char GravEn_SimID[1024]="Waveforms/SG554Q8d9.txt";
86  double SimHrss = 1.213000e-21;
87  double SimEgwR2 = 1.011301e-47;
88  double GravEn_Ampl = 1.213000e-21;
89  double Internal_x = 1.0;
90  double Internal_phi = 0.0;
91  double External_x = burstMDC_theta;
92  double External_phi = burstMDC_phi;
93  double External_psi = burstMDC_psi;
94 
95  double FrameGPS = x->start();
96  //double EarthCtrGPS = FrameGPS+dt*x->size()/2.; // injected in the center of buffer
97  double EarthCtrGPS = 968654066.616913;
98  char SimName[64] = "SG554Q8d9";
99  double SimHpHp = 1.471369e-42;
100  double SimHcHc = 0;
101  double SimHpHc = 0;
102 
103  sprintf(logString,"%s",GravEn_SimID);
104  sprintf(logString,"%s %e",logString,SimHrss);
105  sprintf(logString,"%s %e",logString,SimEgwR2);
106  sprintf(logString,"%s %e",logString,GravEn_Ampl);
107  sprintf(logString,"%s %e",logString,Internal_x);
108  sprintf(logString,"%s %e",logString,Internal_phi);
109  sprintf(logString,"%s %e",logString,External_x);
110  sprintf(logString,"%s %e",logString,External_phi);
111  sprintf(logString,"%s %e",logString,External_psi);
112 
113  sprintf(logString,"%s %10.6f",logString,FrameGPS);
114  sprintf(logString,"%s %10.6f",logString,EarthCtrGPS);
115  sprintf(logString,"%s %s",logString,SimName);
116  sprintf(logString,"%s %e",logString,SimHpHp);
117  sprintf(logString,"%s %e",logString,SimHcHc);
118  sprintf(logString,"%s %e",logString,SimHpHc);
119 
120  double tShift=0;
121  double fPlus=0;
122  double fCross=0;
123  for(int i=0;i<nIFO;i++) {
124  double IFOctrGPS = EarthCtrGPS;
125  if(i>0) IFOctrGPS += gNET.GetDelay(ifos[i].Data(),ifos[0].Data(),phi,theta);
126  double IFOfPlus = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, true);
127  double IFOfCross = gNET.GetAntennaPattern(ifos[i], phi, theta, psi, false);
128  if(ifos[i].CompareTo(ifo)==0) {
129  tShift=IFOctrGPS-EarthCtrGPS;
130  fPlus=IFOfPlus;
131  fCross=IFOfCross;
132  }
133  sprintf(logString,"%s %s %10.6f %e %e",logString,ifos[i].Data(),IFOctrGPS,IFOfPlus,IFOfCross);
134  }
135  cout << logString << endl;
136 
137  net->mdcList.clear();
138  net->mdcList.push_back(logString);
139  net->mdcTime.clear();
140  net->mdcTime.push_back(EarthCtrGPS);
141 
142  cout << endl << endl;
143  for (int nmdc=0; nmdc<(int)net->mdcListSize(); nmdc++) {
144  TString mdcstring(net->getmdcList(nmdc));
145 // if(nmdc==72) {
146  cout << endl << endl;
147  cout << "--------------------------------> " << nmdc << endl;
148  cout << endl << endl;
149  cout << mdcstring.Data() << endl;
150 // }
151  }
152  //gSystem->Exit(0);
153 
154  // read waveform
155  ifstream in;
156  in.open(WAVEFORM_NAME, ios::in);
157  if (!in.good()) {cout << "Error Opening File : " << WAVEFORM_NAME << endl;gSystem->Exit(1);}
158 
159  // fill mdc vector
160  (*x)=0;
161  int offset = (EarthCtrGPS-FrameGPS)*x->rate();
162  int size=0;
163  double hrss=0;
164  while (1) {
165  in >> x->data[offset+size];
166  hrss+=pow(x->data[offset+size],2);
167  if (!in.good()) break;
168  size++;
169  if ((offset+size)>=(int)x->size()) break;
170  }
171  hrss=sqrt(hrss*dt);
172  for(int i=0;i<(int)x->size();i++) x->data[i] = x->data[i]*fPlus*SimHrss/hrss;
173 
174  in.close();
175 
176  if(tShift==0) return;
177 
178  // apply time shift to mdc vector
179  x->FFTW(1);
180  TComplex C;
181  double df = x->rate()/x->size();
182  //double tShift = 10./x->rate();
183  cout << "tShift : " << tShift << endl;
184  for (int ii=0;ii<(int)x->size()/2;ii++) {
185  TComplex X(x->data[2*ii],x->data[2*ii+1]);
186  X=X*C.Exp(TComplex(0.,-2*Pi*ii*df*tShift)); // Time Shift
187  x->data[2*ii]=X.Re();
188  x->data[2*ii+1]=X.Im();
189  }
190  x->FFTW(-1);
191  }
192 
193  return;
194 }
network::ifoName
std::vector< char * > ifoName
Definition: network.hh:591
cfg
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
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
n
int n
Definition: cwb_net.C:10
TString
TString("c")
Definition: cwb_report_skymap.C:111
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::getmdcList
string getmdcList(size_t n)
Definition: network.hh:400
size
Long_t size
Definition: cwb_condor_check.C:48
CWB_Plugin
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
Definition: CWB_Plugin_InjectMDC.C:20
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
network::mdcListSize
size_t mdcListSize()
Definition: network.hh:390
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
WAVEFORM_NAME
#define WAVEFORM_NAME
Definition: CWB_Plugin_InjectMDC.C:17
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
nmdc
int nmdc
Definition: cbc_plots.C:791
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
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
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
GravEn_SimID
char GravEn_SimID[1024]
Definition: CWB_Plugin_SimMDC_SimData.h:16
ifos
TString ifos[60]
Definition: cwb_draw_antpat.C:52