Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ced.cc
Go to the documentation of this file.
1 #include "ced.hh"
2 #include "watplot.hh"
3 #include "wseries.hh"
4 
5 #include "THtml.h"
6 #include "TDocOutput.h"
7 #include "TH2.h"
8 #include "TStyle.h"
9 #include "TCanvas.h"
10 #include "TSystem.h"
11 #include "TMath.h"
12 #include "TMarker.h"
13 #include "TVector3.h"
14 #include "TRotation.h"
15 #include "TPolyLine.h"
16 #include "TMacro.h"
17 #include "Math/Rotation3D.h"
18 #include "Math/Vector3Dfwd.h"
19 #include "TLegend.h"
20 
21 #include "Meyer.hh"
22 #include "gskymap.hh"
23 #include "gnetwork.hh"
24 #include "STFT.hh"
25 #include <string.h>
26 
27 #define REPLACE(STRING,DIR,EXT) TString(STRING).ReplaceAll("/","").ReplaceAll(DIR,"").ReplaceAll("."+TString(EXT),"")
28 
29 using namespace ROOT::Math;
30 
31 ClassImp(CWB::ced) // used by THtml doc
32 
33 void
34 CWB::ced::Write(double factor, size_t iID, int LAG, char* dirCED) {
35 
36  if(!NET || !NET->nLag) exit(1);
37 
38  int i,j,k,ind;
39  int m,K,M,injID;
40  detector* pd;
41  netcluster* p;
42  int nIFO = (int)NET->ifoListSize(); // number of detectors
43  int bLag = LAG<0 ? 0 : LAG;
44  int eLag = LAG<0 ? NET->nLag : LAG+1;
45  wavecomplex Aa, gC;
46  factor = factor<=0 ? 1 : fabs(factor); // used to rescale injected events
47 
48  if(!nIFO) return;
49 
50  // check if it is one detector network
51  if(nIFO==2) if(strcmp(NET->ifoName[0],NET->ifoName[1])==0) nIFO=1;
52 
53  watplot WTS(const_cast<char*>("wtswrc"));
54  watplot PTS(const_cast<char*>("ptswrc"),200,20,800,500);
55  gskymap gSM;
56 
57 // injections
58 
59  injection INJ(nIFO);
60  size_t mdcID, ID;
61  double T, injTime;
62  skymap* psm;
63  std::vector<float>* vP;
64  std::vector<int>* vI;
65 
66  bool ellips = NET->tYPe=='i' || NET->tYPe=='I' ||
67  NET->tYPe=='g' || NET->tYPe=='G' ||
68  NET->tYPe=='s' || NET->tYPe=='S' ||
69  NET->tYPe=='l' || NET->tYPe=='L' ||
70  NET->tYPe=='c' || NET->tYPe=='C' ||
71  NET->tYPe=='e' || NET->tYPe=='E' ||
72  NET->tYPe=='p' || NET->tYPe=='P' ||
73  NET->tYPe=='r' || NET->tYPe=='R';
74 
75  bool burst = NET->tYPe=='b' || NET->tYPe=='B';
76 
77  wavearray<double> skSNR(nIFO);
78  wavearray<double> xkSNR(nIFO);
79 
80 // arrays for cluster parameters
81 
82  wavearray<double> clusterID_net;
83 
84  EVT->run = NET->nRun;
85  EVT->nevent = 0;
86 
87  for(int lag=bLag; lag<eLag; lag++){ // loop on time lags
88 
89  p = NET->getwc(lag); // pointer to netcluster
90  if(!p->size()) continue;
91 
92  if(NET->MRA) clusterID_net = p->get((char*)"ID",0,'S',0);
93  else clusterID_net = p->get((char*)"ID",0,'S',1);
94  K = clusterID_net.size();
95 
96  if(!K) continue;
97 
98 // read cluster parameters
99 
100  for(k=0; k<K; k++) // loop on events
101  {
102  ID = size_t(clusterID_net.data[k]+0.1);
103  if((iID)&&(ID!=iID)) continue;
104  if(ellips) {
105  if(NET->like()!='2') NET->SETNDM(ID,lag,true,NET->MRA ? 0 : 1);
106  }
107  else if(burst)
108  NET->setndm(ID,lag,true,1);
109  else
110  {cout<<"netevent::output(): incorrect search option"<<endl; exit(1);}
111 
112  psm = &(NET->getifo(0)->tau);
113  vI = &(NET->wc_List[lag].p_Ind[ID-1]);
114  ind = (*vI)[0]; // reconstructed sky index
115 
116 // set injections if there are any
117 
118  M = NET->mdc__IDSize();
119  if(!lag) { // only for zero lag
120  injTime = 1.e12;
121  injID = -1;
122  for(m=0; m<M; m++) {
123  mdcID = NET->getmdc__ID(m);
124  T = fabs(EVT->time[0] - NET->getmdcTime(mdcID));
125  if(T<injTime && INJ.fill_in(NET,mdcID)) {
126  injTime = T;
127  injID = mdcID;
128  }
129 // printf("%d %12.3f %12.3f\n",mdcID,NET->getmdcTime(mdcID),T);
130  }
131 
132  if(INJ.fill_in(NET,injID)) { // set injections
133 
134 // printf("injection type: %d %12.3f\n",INJ.type,INJ.time[0]);
135 
136  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
137  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
138  pd = NET->getifo(0);
139  int idSize = pd->RWFID.size();
140  int wfIndex=-1;
141  for (int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==(int)ID) wfIndex=mm;
142 
143  for(j=0; j<int(nIFO); j++) {
144  pd = NET->getifo(j);
145  Aa = pd->antenna(EVT->theta[1],EVT->phi[1]); // inj antenna pattern
146  EVT->hrss[j+nIFO] = INJ.hrss[j]*factor;
147  EVT->bp[j+nIFO] = Aa.real();
148  EVT->bx[j+nIFO] = Aa.imag();
149  EVT->time[j+nIFO] = INJ.time[j];
150  EVT->Deff[j] = INJ.Deff[j]/factor;
151  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",EVT->gps[j]);
152 
153  pwfINJ[j] = INJ.pwf[j];
154  if (pwfINJ[j]==NULL) {
155  cout << "Error : Injected waveform not saved !!! : detector "
156  << NET->ifoName[j] << endl;
157  continue;
158  }
159  if (wfIndex<0) {
160  cout << "Error : Reconstructed waveform not saved !!! : ID -> "
161  << ID << " : detector " << NET->ifoName[j] << endl;
162  continue;
163  }
164  if (wfIndex>=0) pwfREC[j] = pd->RWFP[wfIndex];
165  double R = pd->TFmap.rate();
166  double rFactor = 1.;
167  rFactor *= factor;
168  wavearray<double>* wfINJ = pwfINJ[j];
169  *wfINJ*=rFactor;
170  wavearray<double>* wfREC = pwfREC[j];
171 
172  // rescale data when data are resampled (resample with wavelet change the amplitude)
173  double rescale = 1./pow(sqrt(2.),TMath::Log2(inRate/wfINJ->rate()));
174  *wfINJ*=rescale; *wfREC*=rescale; // rescale data
175 
176  double bINJ = wfINJ->start();
177  double eINJ = wfINJ->start()+wfINJ->size()/R;
178  double bREC = wfREC->start();
179  double eREC = wfREC->start()+wfREC->size()/R;
180  //cout.precision(14);
181  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
182  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
183 
184  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
185  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
186  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
187 
188  double startXCOR = bINJ>bREC ? bINJ : bREC;
189  double endXCOR = eINJ<eREC ? eINJ : eREC;
190  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
191  //cout << "startXCOR : " << startXCOR << " endXCOR : " << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
192 
193  if (sizeXCOR<=0) continue;
194 
195  // the enINJ, enREC, xcorINJ_REC are computed in the INJ range
196 
197  double enINJ=0;
198  for (int i=0;i<(int)wfINJ->size();i++) enINJ+=wfINJ->data[i]*wfINJ->data[i];
199  //for (int i=0;i<sizeXCOR;i++) enINJ+=wfINJ->data[i+oINJ]*wfINJ->data[i+oINJ];
200 
201  double enREC=0;
202  //for (int i=0;i<wfREC->size();i++) enREC+=wfREC->data[i]*wfREC->data[i];
203  for (int i=0;i<sizeXCOR;i++) enREC+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
204 
205  double xcorINJ_REC=0;
206  for (int i=0;i<sizeXCOR;i++) xcorINJ_REC+=wfINJ->data[i+oINJ]*wfREC->data[i+oREC];
207 
208  WSeries<double> wfREC_SUB_INJ;
209  wfREC_SUB_INJ.resize(sizeXCOR);
210  for (int i=0;i<sizeXCOR;i++) wfREC_SUB_INJ.data[i]=wfREC->data[i+oREC]-wfINJ->data[i+oINJ];
211  wfREC_SUB_INJ.start(startXCOR);
212  wfREC_SUB_INJ.rate(wfREC->rate());
213 
214 
215  EVT->iSNR[j] = enINJ;
216  EVT->oSNR[j] = enREC;
217  EVT->ioSNR[j] = xcorINJ_REC;
218 
219  //double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
220  //cout << "enINJ : " << enINJ << " enREC : " << enREC << " xcorINJ_REC : " << xcorINJ_REC << endl;
221  //cout << "erINJ_REC/enINJ : " << erINJ_REC/enINJ << endl;
222 
223  bool batch = gROOT->IsBatch();
224  gROOT->SetBatch(true);
225 
226  // find TF maps to optimal resolution level
227  double optRate = (p->cRate[ID-1])[0];
228  double optLayer = p->rate/optRate;
229  double optLevel = int(log10(optLayer)/log10(2));
230 
231  // DUMP WRC !!!
232  if (wfIndex>=0) {
233  double gmin,gmax;
234  char fname[1024];
235 
236  PTS.canvas->cd();
237 
238  sprintf(fname, "%s/%s_wf_white_inj.dat", dirCED, NET->ifoName[j]);
239  if(sbasedirCED!=NULL) wfINJ->wavearray<double>::Dump(const_cast<char*>(fname),2); // time, amp format
240 // if(sbasedirCED!=NULL) wfINJ->Dump(fname);
241  //cout << "Dump " << fname << endl;
242 
243  sprintf(fname, "%s/%s_wf_white_rec.dat", dirCED, NET->ifoName[j]);
244  if(sbasedirCED!=NULL) wfREC->wavearray<double>::Dump(const_cast<char*>(fname),2); // time, amp format
245 // if(sbasedirCED!=NULL) wfREC->Dump(fname);
246  //cout << "Dump " << fname << endl;
247 
248  sprintf(fname, "%s/%s_wf_white_inj_rec.png", dirCED, NET->ifoName[j]);
249  //cout << "Print " << fname << endl;
250  wfINJ->start(wfINJ->start()-EVT->gps[0]);
251  wfREC->start(wfREC->start()-EVT->gps[0]);
252  double bT,eT;
253  GetBoundaries((wavearray<double>&)*wfINJ, 0.995, bT, eT);
254  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1, bT, eT);
255  // startXCOR-EVT->gps[0], endXCOR-EVT->gps[0]);
256  PTS.graph[0]->SetLineWidth(1);
257  PTS.graph[0]->GetXaxis()->SetTitle(xtitle);
258  if(bT<(startXCOR-EVT->gps[0])) bT=startXCOR-EVT->gps[0];
259  if(eT>(endXCOR-EVT->gps[0])) eT=endXCOR-EVT->gps[0];
260  PTS.plot(wfREC, const_cast<char*>("SAME"), 2, bT, eT);
261  // startXCOR-EVT->gps[0], endXCOR-EVT->gps[0]);
262  PTS.graph[1]->SetLineWidth(2);
263  gmin = TMath::Min(PTS.graph[0]->GetHistogram()->GetMinimum(),PTS.graph[1]->GetHistogram()->GetMinimum());
264  gmax = TMath::Max(PTS.graph[0]->GetHistogram()->GetMaximum(),PTS.graph[1]->GetHistogram()->GetMaximum());
265  PTS.graph[0]->GetYaxis()->SetRangeUser(0.9*gmin,1.1*gmax);
266  wfINJ->start(wfINJ->start()+EVT->gps[0]);
267  wfREC->start(wfREC->start()+EVT->gps[0]);
268  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
269  //sprintf(fname, "%s/%s_wf_white_inj_rec.root", dirCED, NET->ifoName[j]);
270  //PTS.canvas->Print(fname);
271  PTS.clear();
272 
273  if((pd->TFmap.gethigh()-pd->TFmap.getlow())>200) PTS.canvas->SetLogx(true);
274  PTS.canvas->SetLogy(true);
275  sprintf(fname, "%s/%s_wf_white_inj_rec_fft.png", dirCED, NET->ifoName[j]);
276  //cout << "Print " << fname << endl;
277  wfINJ->start(wfINJ->start()-EVT->gps[0]);
278  wfREC->start(wfREC->start()-EVT->gps[0]);
279  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1,
280  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0],
281  true, pd->TFmap.getlow(), pd->TFmap.gethigh());
282  PTS.graph[0]->SetLineWidth(1);
283  PTS.plot(wfREC, const_cast<char*>("SAME"), 2,
284  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0],
285  true, pd->TFmap.getlow(), pd->TFmap.gethigh());
286  PTS.graph[1]->SetLineWidth(2);
287  gmin = TMath::Min(PTS.graph[0]->GetHistogram()->GetMinimum(),PTS.graph[1]->GetHistogram()->GetMinimum());
288  gmax = TMath::Max(PTS.graph[0]->GetHistogram()->GetMaximum(),PTS.graph[1]->GetHistogram()->GetMaximum());
289  PTS.graph[0]->GetYaxis()->SetRangeUser(0.9*gmin,1.1*gmax);
290  wfINJ->start(wfINJ->start()+EVT->gps[0]);
291  wfREC->start(wfREC->start()+EVT->gps[0]);
292  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
293  PTS.clear();
294  PTS.canvas->SetLogx(false);
295  PTS.canvas->SetLogy(false);
296  sprintf(fname, "%s/%s_wf_white_rec_sub_inj.png", dirCED, NET->ifoName[j]);
297  //cout << "Print " << fname << endl;
298  wfREC_SUB_INJ.start(wfREC_SUB_INJ.start()-EVT->gps[0]);
299  wfREC->start(wfREC->start()-EVT->gps[0]);
300  PTS.plot(wfREC, const_cast<char*>("ALP"), 2,
301  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0]);
302  PTS.graph[0]->SetLineWidth(1);
303  PTS.graph[0]->GetXaxis()->SetTitle(xtitle);
304  PTS.plot(wfREC_SUB_INJ, const_cast<char*>("SAME"), 1,
305  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0]);
306  PTS.graph[1]->SetLineWidth(2);
307  gmin = TMath::Min(PTS.graph[0]->GetHistogram()->GetMinimum(),PTS.graph[1]->GetHistogram()->GetMinimum());
308  gmax = TMath::Max(PTS.graph[0]->GetHistogram()->GetMaximum(),PTS.graph[1]->GetHistogram()->GetMaximum());
309  PTS.graph[0]->GetYaxis()->SetRangeUser(0.9*gmin,1.1*gmax);
310  wfREC_SUB_INJ.start(wfREC_SUB_INJ.start()+EVT->gps[0]);
311  wfREC->start(wfREC->start()+EVT->gps[0]);
312  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
313  //sprintf(fname, "%s/%s_wf_white_rec_sub_inj.root", dirCED, NET->ifoName[j]);
314  //PTS.canvas->Print(fname);
315  PTS.clear();
316 
317  if((pd->TFmap.gethigh()-pd->TFmap.getlow())>200) PTS.canvas->SetLogx(true);
318  PTS.canvas->SetLogy(true);
319  sprintf(fname, "%s/%s_wf_white_rec_sub_inj_fft.png", dirCED, NET->ifoName[j]);
320  //cout << "Print " << fname << endl;
321  wfREC_SUB_INJ.start(wfREC_SUB_INJ.start()-EVT->gps[0]);
322  wfREC->start(wfREC->start()-EVT->gps[0]);
323  PTS.plot(wfREC, const_cast<char*>("ALP"), 2,
324  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0],
325  true, pd->TFmap.getlow(), pd->TFmap.gethigh());
326  PTS.graph[0]->SetLineWidth(1);
327  PTS.plot(wfREC_SUB_INJ, const_cast<char*>("SAME"), 1,
328  startXCOR-EVT->gps[0], endXCOR-EVT->gps[0],
329  true, pd->TFmap.getlow(), pd->TFmap.gethigh());
330  PTS.graph[1]->SetLineWidth(2);
331  gmin = TMath::Min(PTS.graph[0]->GetHistogram()->GetMinimum(),PTS.graph[1]->GetHistogram()->GetMinimum());
332  gmax = TMath::Max(PTS.graph[0]->GetHistogram()->GetMaximum(),PTS.graph[1]->GetHistogram()->GetMaximum());
333  PTS.graph[0]->GetYaxis()->SetRangeUser(0.9*gmin,1.1*gmax);
334  wfREC_SUB_INJ.start(wfREC_SUB_INJ.start()+EVT->gps[0]);
335  wfREC->start(wfREC->start()+EVT->gps[0]);
336  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
337  PTS.clear();
338  PTS.canvas->SetLogx(false);
339  PTS.canvas->SetLogy(false);
340 
341  WTS.canvas->cd();
342  //if((pd->TFmap.gethigh()-pd->TFmap.getlow())>200) WTS.canvas->SetLogy(true);
343  // length of temporary buffer for tf plots
344  int xcor_length = sizeXCOR/wfINJ->rate();
345  int wts_size = xcor_length<8 ? 16 : 16*int(1+xcor_length/8.);
346  wts_size*=wfINJ->rate();
347 /*
348  // find optimal level
349  int oRate = int(wfINJ->rate()/EVT->rate[0]);
350  int oLevel=0;
351  while(oRate>1) {oRate/=2;oLevel++;}
352  //cout << "EVT->rate[0] : " << EVT->rate[0] << " oLevel : " << oLevel << endl;
353 */
354  Meyer<double> S(1024,2); // set wavelet for production
355 
356  wavearray<double> xINJ(wts_size);
357  xINJ.start(wfINJ->start()-EVT->gps[0]+double(oINJ+wts_size/2)/wfINJ->rate());
358  xINJ.rate(wfINJ->rate());
359  xINJ=0.;
360  for (int m=0;m<sizeXCOR;m++)
361  if(m<(int)xINJ.size()/2) xINJ.data[m+wts_size/2]=wfINJ->data[m+oINJ];
362  WSeries<double> wINJ(xINJ,S);
363  xINJ.resize(0);
364 
365  double wts_start,wts_stop;
366 
367  //scalogram maps
368  sprintf(fname, "%s/%s_wf_white_inj_tf.png", dirCED, NET->ifoName[j]);
369  //cout << "Print " << fname << endl;
370 // wINJ.Forward(oLevel);
371  if(NET->wdm()) {if(NET->getwdm(optLayer+1)) wINJ.Forward(wINJ,*NET->getwdm(optLayer+1));}
372  else wINJ.Forward(optLevel);
373  wts_start = wINJ.start()+(double)(wts_size/2)/wINJ.rate();
374  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wINJ.rate() : wts_start+(wts_size/2)/wINJ.rate();
375  WTS.plot(&wINJ, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
376  WTS.hist2D->GetYaxis()->SetRangeUser(pd->TFmap.getlow(), pd->TFmap.gethigh());
377  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
378  WTS.clear();
379 
380  wavearray<double> xREC(wts_size);
381  xREC.start(wfREC->start()-EVT->gps[0]+double(oREC+wts_size/2)/wfREC->rate());
382  xREC.rate(wfREC->rate());
383  xREC=0.;
384  for (int m=0;m<sizeXCOR;m++)
385  if(m<(int)xREC.size()/2) xREC.data[m+wts_size/2]=wfREC->data[m+oREC];
386  WSeries<double> wREC(xREC,S);
387  xREC.resize(0);
388 
389  //scalogram maps
390  sprintf(fname, "%s/%s_wf_white_rec_tf.png", dirCED, NET->ifoName[j]);
391  //cout << "Print " << fname << endl;
392 // wREC.Forward(oLevel);
393  if(NET->wdm()) {if(NET->getwdm(optLayer+1)) wREC.Forward(wREC,*NET->getwdm(optLayer+1));}
394  else wREC.Forward(optLevel);
395  wts_start = wREC.start()+(double)(wts_size/2)/wREC.rate();
396  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wREC.rate() : wts_start+(wts_size/2)/wREC.rate();
397  WTS.plot(&wREC, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
398  WTS.hist2D->GetYaxis()->SetRangeUser(pd->TFmap.getlow(), pd->TFmap.gethigh());
399  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
400  WTS.clear();
401 
402  wavearray<double> xDIF(wts_size);
403  xDIF.start(wfREC->start()-EVT->gps[0]+double(oREC+wts_size/2)/wfREC->rate());
404  xDIF.rate(wfREC->rate());
405  xDIF=0.;
406  for (int m=0;m<sizeXCOR;m++)
407  if(m<(int)xDIF.size()/2) xDIF.data[m+wts_size/2]=wfREC->data[m+oREC]-wfINJ->data[m+oINJ];
408  WSeries<double> wDIF(xDIF,S);
409  xDIF.resize(0);
410 
411  //scalogram maps
412  sprintf(fname, "%s/%s_wf_white_dif_tf.png", dirCED, NET->ifoName[j]);
413  //cout << "Print " << fname << endl;
414 // wDIF.Forward(oLevel);
415  if(NET->wdm()) {if(NET->getwdm(optLayer+1)) wDIF.Forward(wDIF,*NET->getwdm(optLayer+1));}
416  else wDIF.Forward(optLevel);
417  wts_start = wDIF.start()+(double)(wts_size/2)/wDIF.rate();
418  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wDIF.rate() : wts_start+(wts_size/2)/wDIF.rate();
419  WTS.plot(&wDIF, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
420  WTS.hist2D->GetYaxis()->SetRangeUser(pd->TFmap.getlow(), pd->TFmap.gethigh());
421  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
422  WTS.clear();
423 
424  double t1 = wDIF.start();
425  double t2 = wDIF.start()+wDIF.size()/wDIF.rate();
426 
427  int ni = 1<<wDIF.pWavelet->m_Level;
428  int nb = int((t1-wDIF.start())*wDIF.rate()/ni);
429  int nj = int((t2-t1)*wDIF.rate())/ni;
430  int ne = nb+nj;
431 
432  int nsize=0;
433  double eeINJ=0;
434  double eeREC=0;
435  double eeDIF=0;
436  wavearray<double> wlINJ;
437  wavearray<double> wlREC;
438  wavearray<double> wlDIF;
439  for(int m=0;m<ni;m++) {
440  wINJ.getLayer(wlINJ,m);
441  wREC.getLayer(wlREC,m);
442  wDIF.getLayer(wlDIF,m);
443  for(int k=nb;k<ne;k++) {
444  if(fabs(wlREC.data[k])>0.1) {
445  eeINJ+=wlINJ.data[k]*wlINJ.data[k];
446  eeREC+=wlREC.data[k]*wlREC.data[k];
447  eeDIF+=wlDIF.data[k]*wlDIF.data[k];
448  nsize++;
449  }
450  }
451  }
452 
453  //cout << j << " tfNorm : " << eeINJ << " " << " " << eeREC << " " << eeDIF/eeINJ << " " << nsize << endl;
454 
455  wINJ.resize(0);
456  wREC.resize(0);
457  wDIF.resize(0);
458 
459  WTS.canvas->SetLogy(false);
460  gROOT->SetBatch(batch);
461  }
462  *wfINJ*=1./rFactor;
463  *wfINJ*=1./rescale; *wfREC*=1./rescale; // restore original rescaled data
464  }
465  }
466  }
467 /*
468  if(EVT->fP!=NULL) {
469  fprintf(EVT->fP,"\n# trigger %d in lag %d for \n",int(ID),int(n));
470  EVT->Dump();
471  vP = &(NET->wc_List[n].p_Map[ID-1]);
472  vI = &(NET->wc_List[n].p_Ind[ID-1]);
473  x = cos(psm->theta_1*PI/180.)-cos(psm->theta_2*PI/180.);
474  x*= (psm->phi_2-psm->phi_1)*180/PI/psm->size();
475  fprintf(EVT->fP,"sky_res: %f\n",sqrt(fabs(x)));
476  fprintf(EVT->fP,"map_lenght: %d\n",int(vP->size()));
477  fprintf(EVT->fP,"#skyID theta DEC step phi R.A step probability cumulative\n");
478  x = 0.;
479  for(j=0; j<int(vP->size()); j++) {
480  i = (*vI)[j];
481  x+= (*vP)[j];
482  fprintf(EVT->fP,"%6d %5.1f %5.1f %6.2f %5.1f %5.1f %6.2f %e %e\n",
483  int(i),psm->getTheta(i),psm->getDEC(i),psm->getThetaStep(i),
484  psm->getPhi(i),psm->getRA(i),psm->getPhiStep(i),(*vP)[j],x);
485  }
486  }
487 */
488  // dump to fits the probability skymap (only for healpix skymap)
489  if(nIFO>1) {
490  bool batch = gROOT->IsBatch();
491  gROOT->SetBatch(true);
492 
493  // Dump2fits probability skymap (healpix)
494  skymap skyprobcc = *psm;
495  skyprobcc=0.;
496  skymap skyprob = *psm;
497  skyprob=1.e-12;
498 
499  vP = &(NET->wc_List[lag].p_Map[ID-1]);
500  vI = &(NET->wc_List[lag].p_Ind[ID-1]);
501  double th,ph,ra;
502  int k;
503  for(j=0; j<int(vP->size()); j++) {
504  i = (*vI)[j];
505  th = skyprob.getTheta(i);
506  ph = skyprob.getPhi(i);
507 
508  k=skyprob.getSkyIndex(th, ph);
509  skyprob.set(k,(*vP)[j]);
510 
511  ra = skyprob.getRA(i);
512  k=skyprob.getSkyIndex(th, ra);
513  skyprobcc.set(k,(*vP)[j]);
514  }
515 
516  char fname[1024];
517 #ifdef _USE_HEALPIX
518  if(skyprobcc.getType()&&(sbasedirCED!=NULL)) { // healpix in celestial coordinates
519  sprintf(fname, "%s/skyprobcc.%s", dirCED, "fits");
520 
521  // build fits configur info
522  TString analysis = "1G";
523  if(NET->like()=='2') analysis="2G";
524  if(NET->MRA) analysis+=":MRA";
525  if(NET->pattern==0) analysis+=":Packet(0)";
526  if((NET->pattern!=0 && NET->pattern<0)) analysis+=TString::Format(":Packet(%d)",NET->pattern);
527  if((NET->pattern!=0 && NET->pattern>0)) analysis+=TString::Format(":Packet(+%d)",NET->pattern);
528 
529  char configur[64]="";
530  char search = NET->tYPe;
531  if (search=='r') sprintf(configur,"%s un-modeled",analysis.Data());
532  if(analysis=="1G") {
533  if((search=='i')||(search=='I')) sprintf(configur,"%s elliptical",analysis.Data());
534  if((search=='s')||(search=='S')) sprintf(configur,"%s linear",analysis.Data());
535  if((search=='g')||(search=='G')) sprintf(configur,"%s circular",analysis.Data());
536  } else {
537  if (search=='i') sprintf(configur,"%s iota-wave",analysis.Data());
538  if (search=='p') sprintf(configur,"%s psi-wave",analysis.Data());
539  if((search=='l')||(search=='s')) sprintf(configur,"%s linear",analysis.Data());
540  if((search=='c')||(search=='g')) sprintf(configur,"%s circular",analysis.Data());
541  if((search=='e')||(search=='b')) sprintf(configur,"%s elliptical",analysis.Data());
542  }
543  skyprobcc.Dump2fits(fname,EVT->time[0],configur,const_cast<char*>("PROB"),const_cast<char*>("pix-1"),'C');
544  }
545 #endif
546 
547  // dump skyprob plot in cc coordinates
548  gSM=skyprobcc;
549  gSM.SetOptions("hammer","Celestial",2);
550 
551  int L = gSM.size(); // number of pixels in the sphere
552  double pi = TMath::Pi();
553  double S = 4*pi*pow(180/pi,2); // solid angle of a sphere
554  double dS = S/L; // solid angle of a pixel
555 
556  for(int l=0;l<L;l++) { // loop over the sky grid
557  double prob = gSM.get(l); // get probability per pixel
558  prob/=dS; // normalize probability to prob per deg^2
559  if(prob==0) gSM.set(l,1e-40); // set !=0 to force dark blue background
560  }
561 
562  gSM.SetZaxisTitle("prob. per deg^{2}");
563  gSM.Draw(57);
564  TH2D* hsm = gSM.GetHistogram();
565  hsm->GetZaxis()->SetTitleOffset(0.85);
566  hsm->GetZaxis()->SetTitleSize(0.03);
567  sprintf(fname, "%s/skyprobcc.%s", dirCED, "png");
568  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
569 
570  gROOT->SetBatch(batch);
571  }
572  }
573  }
574 }
575 
576 int
577 CWB::ced::Write(double factor, bool fullCED) {
578 
579  int i,j,k,n,m,l;
580  int nIFO = NET->ifoListSize();
581  int K = NET->nLag;
582  int M = NET->mdc__ID.size();
583  if(M==0&&simulation) {
584  cout << "CWB::ced::Write : No injected events found -> force simulation=0" << endl;
585  simulation=0;
586  }
587  int ID;
588  char home_ced_path[256]="~waveburst/WWW/LSC/waveburst/ced";
589  FILE* hP = NULL; // html file
590  char search = NET->tYPe;
591  TString analysis = "1G";
592  if(NET->like()=='2') analysis="2G";
593  if(NET->optim) analysis+=":SRA";
594  else analysis+=":MRA";
595  if(NET->pattern==0) analysis+=":Packet(0)";
596  if((NET->pattern!=0 && NET->pattern<0)) analysis+=TString::Format(":Packet(%d)",NET->pattern);
597  if((NET->pattern!=0 && NET->pattern>0)) analysis+=TString::Format(":Packet(+%d)",NET->pattern);
598  int wmod = NET->optim ? 1 : 0; // gwtMRAwave mode
599 
600  wavearray<double> id;
601  wavearray<double> tm;
602 
603  bool batch = gROOT->IsBatch();
604  gROOT->SetBatch(true);
605 
606  watplot PCH(const_cast<char*>("chirp"));
607  watplot WTS(const_cast<char*>("wts"));
608  watplot PTS(const_cast<char*>("pts"),200,20,800,500);
609  gskymap gSM; gSM.SetOptions("","Geographic",2);
610  gnetwork gNET(NET); // used to plot circles
611 
612  int optRate,optLayer,optLevel;
613 
614  TMarker inj; inj.SetMarkerStyle(29); // injected pos (white star)
615  TMarker INJ; INJ.SetMarkerStyle(30); // injected pos (empty black star)
616  TMarker rec; rec.SetMarkerStyle(29); // recostructed pos (black star)
617  TMarker REC; REC.SetMarkerStyle(30); // recostructed pos (emply white star)
618  TMarker det; det.SetMarkerStyle(20); // detected pos (black dot )
619  TMarker DET; DET.SetMarkerStyle(24); // detected pos (empty white dot )
620 
621  double r2d = 180./TMath::Pi();
622  double d2r = TMath::Pi()/180.;
623 
624  char dirCED[1024]="";
625  char command[1024];
626  char ifostr[20]="";
627  char fname[1024];
628 
629  int minTimeDet=nIFO;
630  double minTime=1e40;
631  double eventTime[NIFO];
632  double lagTime[NIFO];
633  double startSegTime,stopSegTime;
634  int ifoid[NIFO],sortid[NIFO];
635  TString xtitle[NIFO];
636  detector *pD[NIFO];
637 
638  int rate = 1;
639  if(NET->like()=='2') rate = NET->MRA ? 0 : 1; // likelihood2G
640  else rate = int(2*NET->getifo(0)->TFmap.resolution(0)+0.5); // likelihoodI
641 
642  //Fill in all skymaps
643  double old_cc = NET->netCC;
644  double old_rho = NET->netRHO;
645  NET->netCC = -1;
646  NET->netRHO = 0;
647 
648  // check if it is one detector network
649  if(nIFO==2) if(strcmp(NET->ifoName[0],NET->ifoName[1])==0) nIFO=1;
650 
651  for(n=0; n<nIFO; n++) {
652  sprintf(ifostr,"%s%s",ifostr,NET->ifoName[n]);
653  pD[n] = NET->getifo(n);
654  ifoid[n]=n;
655  }
656  TMath::Sort(nIFO,ifoid,sortid,false);
657 
658  // create sbasedirCED if selected
659  if(sbasedirCED!=NULL) {
660  sprintf(command,"mkdir -p %s",sbasedirCED);
661  if(gSystem->Exec(command)) return 0;
662  } else {
663  rbasedirCED->cd();
664  }
665 
666  for(k=0; k<K; k++) { // loop over the lags
667 
668  id = NET->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
669  tm = NET->getwc(k)->get(const_cast<char*>("TIME"), 0, 'L', 0);
670 
671 
672  for(j=0; j<(int)id.size(); j++) { // loop over cluster index
673 
674  ID = size_t(id.data[j]+0.5);
675 
676  if(NET->wdm() && NET->getwc(k)->sCuts[ID-1]!=-1) continue; // skip rejected/processed clusters
677 
678  if(NET->like()=='2') EVT->output2G(NULL,NET,ID,k,factor);
679  else EVT->output(NULL,NET,factor,ID,k);
680 
681  if(EVT->rho[analysis=="1G"?1:0] < rho) continue;
682 
683  int masterDet=0;
684  int lagMin=2147483647;
685  for(n=0; n<nIFO; n++) if(EVT->lag[n]<lagMin) {lagMin=int(EVT->lag[n]);masterDet=n;}
686 
687  // create ced directory & index.html link
688  strcpy(ifostr,"");
689  for(n=0; n<nIFO; n++) sprintf(ifostr,"%s%s",ifostr,NET->ifoName[n]);
690  if(sbasedirCED!=NULL) {
691  //if(nIFO==1&&strlen(chName)>1) sprintf(dirCED, "%s/%s_%s", sbasedirCED, ifostr,chName);
692  //else sprintf(dirCED, "%s/%s", sbasedirCED, ifostr);
693  sprintf(dirCED, "%s/%s", sbasedirCED, ifostr);
694  for(n=0; n<nIFO; n++) sprintf(dirCED, "%s_%.3f",dirCED,EVT->start[n]);
695  sprintf(command,"mkdir -p %s",dirCED);
696  if(gSystem->Exec(command)) return 0;
697 
698  if (getenv("HOME_CED_PATH")!=NULL) {
699  sprintf(home_ced_path,"%s",getenv("HOME_CED_PATH"));
700  }
701  strcpy(ifostr,"");
702  for(n=0; n<nIFO; n++) sprintf(ifostr,"%s%s",ifostr,NET->ifoName[sortid[n]]);
703  if(!simulation) sprintf(command,"ln -s %s/index/cedindex_%s.html %s/index.html",home_ced_path,ifostr,dirCED);
704  else sprintf(command,"ln -s %s/index/cedindex_sim_%s.html %s/index.html",home_ced_path,ifostr,dirCED);
705  if(gSystem->Exec(command)) return 0;
706  } else {
707  //if(nIFO==1&&strlen(chName)>1) sprintf(dirCED, "%s_%s", ifostr,chName);
708  //else sprintf(dirCED, "%s", ifostr);
709  sprintf(dirCED, "%s", ifostr);
710  for(n=0; n<nIFO; n++) sprintf(dirCED, "%s_%.3f",dirCED,EVT->start[n]);
711  rbasedirCED->cd(); rbasedirCED->mkdir(dirCED)->cd();
712  }
713 
714  if(NET->like()=='2') EVT->output2G(NULL,NET,ID,k,factor);
715  else EVT->output(NULL,NET,factor,ID,k);
716 
717  //Write(factor,ID,k);
718  if(EVT->rho[analysis=="1G"?1:0] < rho) continue;
719 
720  double bPP =0.01;
721  double TH = 0.2;
722  TH2F* hist = NULL;
723  if(NET->like()=='2') ; // likelihood2G
724  else NET->likelihood(search, NET->acor, ID, k);
725 
726  // dump event file & get event parameters
727  if(sbasedirCED!=NULL) {
728  sprintf(fname, "%s/eventDump.txt", dirCED);
729  } else {
730  gRandom->SetSeed(0);
731  int rnID = int(gRandom->Rndm(13)*1.e9);
732  UserGroup_t* uinfo = gSystem->GetUserInfo();
733  TString uname = uinfo->fUser;
734  gSystem->Exec(TString("mkdir -p /dev/shm/")+uname);
735  sprintf(fname,"/dev/shm/%s/eventDump-%d.txt",uname.Data(),rnID);
736  }
737  EVT->dopen(fname,const_cast<char*>("w"));
738 
739  if(NET->like()=='2') EVT->output2G(NULL,NET,ID,k,factor);
740  else EVT->output(NULL,NET,factor,ID,k);
741 
742  Write(factor,ID,k,dirCED);
743  EVT->dclose();
744  if(sbasedirCED==NULL) {
745  TMacro macro(fname);
746  macro.Write("eventDump.txt");
747  gSystem->Exec(TString("rm ")+fname);
748  }
749 
750  // compute event time & lags time
751  for(n=0; n<nIFO; n++) eventTime[n]=(EVT->start[n]+EVT->stop[n])/2.;
752  for(n=0; n<nIFO; n++) minTime = ((eventTime[n]-EVT->gps[n])<minTime) ? eventTime[n]-EVT->gps[n] : minTime;
753  for(n=0; n<nIFO; n++) minTimeDet = ((eventTime[n]-EVT->gps[n])<=minTime) ? n : minTimeDet;
754  for(n=0; n<nIFO; n++) lagTime[n]=eventTime[n]-EVT->gps[minTimeDet]-minTime;
755  startSegTime= EVT->gps[minTimeDet];
756  // NOTE : we use EVT->duration[1] because it caintains the event stop-start
757  // while EVT->duration[0] contains the duration estimated with rms
758  stopSegTime = EVT->gps[minTimeDet]+EVT->left[minTimeDet]+EVT->right[minTimeDet]+EVT->duration[1];
759  for(n=0; n<nIFO; n++) xtitle[n] = TString::Format("Time (sec) : GPS OFFSET = %.3f",EVT->gps[n]);
760 
761  // create jobSummary.html & eventSummary.html
762  if(sbasedirCED!=NULL) sprintf(fname, "%s/jobSummary.html", dirCED);
763  else sprintf(fname, "/dev/null");
764  if((hP = fopen(fname, "w")) != NULL) {
765  fprintf(hP,"<table border=1 cellspacing=0 width=100%% height=100%%>\n");
766  fprintf(hP,"<tr align=\"center\">\n");
767  if(nIFO==1) {
768  fprintf(hP,"<th>DETECTOR</th>\n");
769  } else {
770  fprintf(hP,"<th>NETWORK</th>\n");
771  }
772  fprintf(hP,"<td>%s</td>\n",ifostr);
773  fprintf(hP,"</tr>\n");
774  fprintf(hP,"<tr align=\"center\">\n");
775  if(nIFO==1) {
776  fprintf(hP,"<th>CHANNEL</th>\n");
777  fprintf(hP,"<td>%s</td>\n",chName);
778  fprintf(hP,"</tr>\n");
779  fprintf(hP,"<tr align=\"center\">\n");
780  fprintf(hP,"<th>SEARCH</th>\n");
781  fprintf(hP,"<td>%s</td>\n",analysis.Data());
782  } else {
783  fprintf(hP,"<th>SEARCH</th>\n");
784  if(analysis=="1G") {
785  if(search=='E' || search=='E') fprintf(hP,"<td>%s un-modeled (%c)</td>\n",analysis.Data(),search);
786  if(search=='b' || search=='B') fprintf(hP,"<td>%s un-modeled (%c)</td>\n",analysis.Data(),search);
787  if(search=='r' || search=='R') fprintf(hP,"<td>%s un-modeled (%c)</td>\n",analysis.Data(),search);
788  if(search=='i' || search=='I') fprintf(hP,"<td>%s elliptical (%c)</td>\n",analysis.Data(),search);
789  if(search=='g' || search=='G') fprintf(hP,"<td>%s circular (%c)</td>\n",analysis.Data(),search);
790  if(search=='s' || search=='S') fprintf(hP,"<td>%s linear (%c)</td>\n",analysis.Data(),search);
791  } else {
792  char _search = std::tolower(search);
793  if(_search=='r') fprintf(hP,"<td>%s un-modeled(%c)</td>\n",analysis.Data(),search);
794  if(_search=='i') fprintf(hP,"<td>%s iota-wave(%c)</td>\n",analysis.Data(),search);
795  if(_search=='p') fprintf(hP,"<td>%s psi-wave(%c)</td>\n",analysis.Data(),search);
796  if(_search=='e'||_search=='b') fprintf(hP,"<td>%s elliptical(%c)</td>\n",analysis.Data(),search);
797  if(_search=='c'||_search=='g') fprintf(hP,"<td>%s circular(%c)</td>\n",analysis.Data(),search);
798  if(_search=='l'||_search=='s') fprintf(hP,"<td>%s linear(%c)</td>\n",analysis.Data(),search);
799  }
800  }
801  fprintf(hP,"</tr>\n");
802  fprintf(hP,"<tr align=\"center\">\n");
803  fprintf(hP,"<th>START SEGMENT</th>");
804  fprintf(hP,"<td>%9.3f</td>\n",startSegTime);
805  fprintf(hP,"</tr>\n");
806  fprintf(hP,"<tr align=\"center\">\n");
807  fprintf(hP,"<th>STOP SEGMENT</th>");
808  fprintf(hP,"<td>%9.3f</td>\n",stopSegTime);
809  fprintf(hP,"</tr>\n");
810  if(simulation) {
811  fprintf(hP,"<tr align=\"center\">\n");
812  fprintf(hP,"<th>MDC</th>\n");
813  fprintf(hP,"<td>%s</td>\n",NET->getmdcType(EVT->type[1]-1).c_str());
814  fprintf(hP,"</tr>\n");
815  }
816  if(!simulation&&nIFO>1) {
817  fprintf(hP,"<tr align=\"center\">\n");
818  fprintf(hP,"<th>LAG</th>\n");
819  fprintf(hP,"<td>");
820  for(n=0; n<nIFO-1; n++) fprintf(hP,"%3.3f / ",lagTime[n]);
821  fprintf(hP,"%3.3f",lagTime[nIFO-1]);
822  fprintf(hP,"</td>\n");
823  fprintf(hP,"</tr>\n");
824  }
825  fprintf(hP,"</table>\n");
826  fclose(hP);hP = NULL;
827  } else {
828  cout << "netevent::ced() error: cannot open file " << fname <<". \n";
829  }
830 
831  // convert user macro into html and copy to ced dir
832  if(sbasedirCED!=NULL) {
833  THtml html;
834  html.SetEtcDir(gSystem->ExpandPathName("$HOME_WAT/html/etc/html"));
835  html.SetProductName("CED");
836  TString html_input_dir=dirCED;
837  html.SetInputDir(html_input_dir.Data());
838 
839  char cmd[1024];
840  sprintf(cmd,"cp %s/html/etc/html/ROOT.css %s/",gSystem->ExpandPathName("$HOME_WAT"),dirCED);
841  gSystem->Exec(cmd);
842  sprintf(cmd,"cp %s/html/etc/html/ROOT.js %s/",gSystem->ExpandPathName("$HOME_WAT"),dirCED);
843  gSystem->Exec(cmd);
844 
845  TString cwb_parameters_file;
846  if(gSystem->Getenv("CWB_PARAMETERS_FILE")==NULL) {
847  cout << "Error : environment CWB_PARAMETERS_FILE is not defined!!!" << endl;exit(1);
848  } else {
849  cwb_parameters_file=TString(gSystem->Getenv("CWB_PARAMETERS_FILE"));
850  }
851  TString cwb_uparameters_file;
852  if(gSystem->Getenv("CWB_UPARAMETERS_FILE")==NULL) {
853  cout << "Error : environment CWB_UPARAMETERS_FILE is not defined!!!" << endl;exit(1);
854  } else {
855  cwb_uparameters_file=TString(gSystem->Getenv("CWB_UPARAMETERS_FILE"));
856  }
857 
858  // redirect stderr to /dev/null to getrid of messages produced by html.Convert
859  fpos_t poserr; fflush(stderr); fgetpos(stderr, &poserr);
860  int fderr = dup(fileno(stderr)); freopen("/dev/null", "w", stderr);
861  // redirect stdout to /dev/null to getrid of messages produced by html.Convert
862  fpos_t posout; fflush(stdout); fgetpos(stdout, &posout);
863  int fdout = dup(fileno(stdout)); freopen("/dev/null", "w", stdout);
864 
865  html.Convert(cwb_parameters_file.Data(),cwb_parameters_file.Data(),dirCED,"");
866  if(analysis=="1G") sprintf(cmd,"mv %s/cwb1G_parameters.C.html %s/cwb_parameters.C.html",dirCED,dirCED);
867  else sprintf(cmd,"mv %s/cwb2G_parameters.C.html %s/cwb_parameters.C.html",dirCED,dirCED);
868  gSystem->Exec(cmd);
869 
870  char work_dir[1024]; sprintf(work_dir,"%s",gSystem->WorkingDirectory()); // working dir
871  TString cwb_uparameters_path = TString(work_dir)+"/"+cwb_uparameters_file;
872  if(cwb_uparameters_file!="") {
873  html.Convert(cwb_uparameters_path.Data(),cwb_uparameters_file.Data(),dirCED,"");
874  // rename file to the standard user_parameters.C name
875  TString cwb_iparameters_path = TString(dirCED)+"/"+gSystem->BaseName(cwb_uparameters_path.Data())+".html";
876  TString cwb_oparameters_path = TString(dirCED)+"/user_parameters.C.html";
877  gSystem->Rename(cwb_iparameters_path.Data(),cwb_oparameters_path.Data());
878  }
879 
880  // restore the stderr output
881  fflush(stderr); dup2(fderr, fileno(stderr)); close(fderr);
882  clearerr(stderr); fsetpos(stderr, &poserr);
883  // restore the stdout output
884  fflush(stdout); dup2(fdout, fileno(stdout)); close(fdout);
885  clearerr(stdout); fsetpos(stdout, &posout);
886  }
887 
888  // create jobSummary.html & eventSummary.html
889  if(sbasedirCED!=NULL) sprintf(fname, "%s/eventSummary.html", dirCED);
890  else sprintf(fname, "/dev/null");
891  if((hP = fopen(fname, "w")) != NULL) {
892  fprintf(hP,"<table border=1 cellspacing=0 width=100%% height=100%%>\n");
893  fprintf(hP,"<tr align=\"center\">\n");
894  fprintf(hP,"<th>GPS TIME</th>\n");
895  fprintf(hP,"<th>SNR</th>\n");
896  if(nIFO==1) {
897  fprintf(hP,"<th>DURATION</th>\n");
898  fprintf(hP,"<th>FREQUENCY</th>\n");
899  fprintf(hP,"<th>BANDWIDTH</th>\n");
900  } else {
901  fprintf(hP,"<th>RHO[0/1]</th>\n");
902  fprintf(hP,"<th>CC[0/1/2/3]</th>\n");
903  fprintf(hP,"<th>ED</th>\n");
904  fprintf(hP,"<th>PHI</th>\n");
905  fprintf(hP,"<th>THETA</th>\n");
906  }
907  fprintf(hP,"</tr>\n");
908  fprintf(hP,"<tr align=\"center\">\n");
909  fprintf(hP,"<td>%9.3f</td>\n",EVT->time[minTimeDet]);
910  if(nIFO==1) {
911  fprintf(hP,"<td>%4.1f</td>\n",sqrt(EVT->likelihood/2.));
912  fprintf(hP,"<td>%3.3f</td>\n",EVT->duration[0]);
913  fprintf(hP,"<td>%3.1f</td>\n",EVT->frequency[0]);
914  fprintf(hP,"<td>%3.1f</td>\n",EVT->bandwidth[0]);
915  } else {
916  fprintf(hP,"<td>%4.1f</td>\n",sqrt(EVT->likelihood));
917  fprintf(hP,"<td>%3.1f / %3.1f</td>\n",EVT->rho[0],EVT->rho[1]);
918  fprintf(hP,"<td>%1.2f / %1.2f / %1.2f / %1.2f</td>\n",
919  EVT->netcc[0],EVT->netcc[1],EVT->netcc[2],EVT->netcc[3]);
920  fprintf(hP,"<td>%1.2f</td>\n",EVT->neted[0]/EVT->ecor);
921  fprintf(hP,"<td>%3.1f</td>\n",EVT->phi[0]);
922  fprintf(hP,"<td>%3.1f</td>\n",90.-EVT->theta[0]);
923  }
924  fprintf(hP,"</tr>\n");
925  fprintf(hP,"</table>\n");
926  fclose(hP);hP = NULL;
927  } else {
928  cout << "netevent::ced() error: cannot open file " << fname <<". \n";
929  }
930 
931  // plots nRMS
932  int nIFO_RMS=0;
933  std::vector<TGraph*> mgraph(nIFO);
934  Color_t psd_color[8] = {kBlue, kRed, kGreen, kBlack, 6, 3, 8, 43};
935  for(n=0; n<nIFO; n++) {
936 
937  if(pD[n]->nRMS.size()==0) continue; // it is zero when CED is produced starting from SUPERCLUSTER stage
938 
939  nIFO_RMS++;
940 
941  WSeries<double> nRMS = pD[n]->nRMS;
942 
943  double fLow = pD[n]->getTFmap()->getlow();
944  double fHigh = pD[n]->getTFmap()->gethigh();
945  double rate = pD[n]->getTFmap()->rate();
946  if(fLow<16) fLow=16;
947  if(fHigh>rate/2.) fHigh=rate/2.;
948 
949  WDM<double>* wdm = (WDM<double>*) pD[n]->nRMS.pWavelet;
950  int M = pD[n]->nRMS.getLevel();
951  double* map00 = wdm->pWWS;
952  double dF = rate/M/2./2.; // psd frequency resolution
953 
954  wavearray<double> psd(2*M); // PSD @ event time
955  psd.start(0);
956  psd.rate(1./dF);
957 
958  wavearray<double> PSD=psd; // average PSD
959  PSD=0;
960 
961  int nPSD = pD[n]->nRMS.size()/(M+1); // number of psd in the nRMS object
962  double etime = eventTime[n]-EVT->gps[n]; // event time from the beginning of the buffer
963 
964  double segLen = stopSegTime-startSegTime;
965  int I = (int)etime/(segLen/nPSD); // psd event index
966 
967  for(int i=0; i<nPSD; ++i){
968  if(i==I) psd[0] = map00[0]; // first half-band
969  PSD[0]+= pow(map00[0],2); // first half-band
970  for(int j=1; j<M; j++){
971  if(i==I) psd[2*j-1] = map00[j];
972  if(i==I) psd[2*j] = map00[j];
973  PSD[2*j-1]+= pow(map00[j],2);
974  PSD[2*j] += pow(map00[j],2);
975  }
976  if(i==I) psd[2*M-1] = map00[M]; // last half-band
977  PSD[2*M-1]+= pow(map00[M],2); // last half-band
978  map00+=M+1;
979  }
980  // NOTE : check if the following normalization is correct when cfg->fResample>0
981  psd*=sqrt(2.)/sqrt(inRate); // oneside psd normalization
982  for(int j=0; j<PSD.size(); j++) PSD[j] = sqrt(PSD[j]/nPSD)*sqrt(2.)/sqrt(inRate); // oneside PSD normalization
983 
984  PTS.canvas->cd();
985  gStyle->SetTitleFont(12,"D");
986  sprintf(fname, "%s/%s_psd.%s", dirCED, NET->ifoName[n], gtype);
987  PTS.plot((wavearray<double>&)psd, const_cast<char*>("ALP"), 1, fLow+2*dF, fHigh-2*dF);
988  //PTS.plot((wavearray<double>&)PSD, const_cast<char*>("ALP"), 1, fLow+2*dF, fHigh-2*dF);
989  PTS.graph[0]->SetTitle(TString::Format("Power Spectral Density after lines' removal @ Time (sec) = %.0f",EVT->time[minTimeDet]));
990  PTS.graph[0]->GetXaxis()->SetTitle("Frequency (Hz) ");
991  PTS.graph[0]->GetYaxis()->SetTitle("[strain / #sqrt{Hz}] ");
992  PTS.graph[0]->SetMarkerStyle(1);
993  PTS.graph[0]->SetMinimum(1.e-24);
994  PTS.graph[0]->SetMaximum(1.e-20);
995  PTS.graph[0]->SetLineWidth(2);
996  PTS.graph[0]->SetLineColor(psd_color[n%8]);
997  if(TString(NET->ifoName[n])=="L1") PTS.graph[0]->SetLineColor(TColor::GetColor("#66ccff"));
998  if(TString(NET->ifoName[n])=="L1") PTS.graph[0]->SetMarkerColor(TColor::GetColor("#66ccff"));
999  if(TString(NET->ifoName[n])=="H1") PTS.graph[0]->SetLineColor(kRed);
1000  if(TString(NET->ifoName[n])=="H1") PTS.graph[0]->SetMarkerColor(kRed);
1001  if(TString(NET->ifoName[n])=="V1") PTS.graph[0]->SetLineColor(TColor::GetColor("#9900cc"));
1002  if(TString(NET->ifoName[n])=="V1") PTS.graph[0]->SetMarkerColor(TColor::GetColor("#9900cc"));
1003  PTS.canvas->SetLogx(true);
1004  PTS.canvas->SetLogy(true);
1005  // draw the legend
1006  TLegend leg(0.753012,0.8-0.01,0.8885542,0.8793706,NULL,"brNDC");
1007  leg.SetTextAlign(22);
1008  leg.SetLineColor(kBlack);
1009  leg.AddEntry(PTS.graph[0],TString::Format("%s",NET->ifoName[n]),"lp");
1010  leg.Draw();
1011  //PTS.plot((wavearray<double>&)psd, const_cast<char*>("SAME"), 2, fLow+2*dF, fHigh-2*dF);
1012  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1013  mgraph[n] = PTS.graph[0];
1014  PTS.graph.clear();
1015  //PTS.clear();
1016 
1017  sprintf(fname, "%s/%s_psd.dat", dirCED, NET->ifoName[n]);
1018  if(sbasedirCED!=NULL) psd.wavearray<double>::Dump(const_cast<char*>(fname),2); // save freq, psd format into ascii file
1019  }
1020  if(nIFO_RMS==nIFO) { // nRMS is present for all detectors -> we can print the combined PSD
1021  PTS.graph=mgraph;;
1022  PTS.graph[0]->Draw("ALP"); for(n=1; n<nIFO; n++) PTS.graph[n]->Draw("SAME");
1023  // draw the legend
1024  TLegend leg(0.753012,0.8-nIFO*0.01,0.8885542,0.8793706,NULL,"brNDC");
1025  leg.SetTextAlign(22);
1026  leg.SetLineColor(kBlack);
1027  for(int n=0;n<nIFO;n++) leg.AddEntry(PTS.graph[n],TString::Format("%s",NET->ifoName[n]),"lp");
1028  leg.Draw();
1029  sprintf(fname, "%s/NET_psd.%s", dirCED, gtype);
1030  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1031  }
1032  PTS.canvas->SetLogx(false);
1033  PTS.canvas->SetLogy(false);
1034  PTS.clear();
1035  mgraph.clear();
1036 
1037  // get TF maps @ optimal resolution level
1038  optRate = (NET->getwc(k)->cRate[ID-1])[0];
1039  optLayer = NET->getwc(k)->rate/optRate;
1040  optLevel = int(log10(optLayer)/log10(2));
1041  WSeries<double>* pTF[NIFO_MAX];
1042  for(int n=0; n<nIFO; n++) {
1043  pTF[n] = pD[n]->getTFmap();
1044  if(useSparse&&(analysis!="1G")) { // data not present in TF map (2G phase2), the sparse TF maps are used
1045  int optSparse = pD[n]->getSTFind(optRate); // get optimal sparse map index
1046  pD[n]->vSS[optSparse].Expand(true); // rebuild wseries from sparse table
1047  pTF[n] = new WSeries<double>(*(WSeries<double>*)&(pD[n]->vSS[optSparse]));
1048  pTF[n]->setlow(pD[n]->getTFmap()->getlow());
1049  pTF[n]->sethigh(pD[n]->getTFmap()->gethigh());
1050  }
1051  // rescale data when data are resampled (resample with wavelet change the amplitude)
1052  double rescale = 1./pow(sqrt(2.),TMath::Log2(inRate/pTF[n]->rate()));
1053  *pTF[n]*=rescale; // rescale TF map
1054  }
1055 
1056  // plot spectrograms
1057  for(int n=0; n<nIFO; n++) {
1058  if(pTF[n]->size()==0) continue;
1059  if(pTF[n]->getLevel()>0) pTF[n]->Inverse();
1060 
1061  int nfact=4;
1062  int nfft=nfact*512;
1063  int noverlap=nfft-10;
1064  double fparm=nfact*6;
1065  //int ystart = (EVT->start[n]-pTF[n]->start()-1)*pTF[n]->rate();
1066  //int ystop = (EVT->stop[n]-pTF[n]->start()+1)*pTF[n]->rate();
1067  int ystart = int((EVT->start[n]-EVT->gps[n]-1)*pTF[n]->rate());
1068  int ystop = int((EVT->stop[n]-EVT->gps[n]+1)*pTF[n]->rate());
1069  ystart-=nfft;
1070  ystop+=nfft;
1071  int ysize=ystop-ystart;
1072  wavearray<double> y;y.resize(ysize);y.rate(pTF[n]->rate());y.start(ystart/pTF[n]->rate());
1073 
1074  // stft use dt=y.rate() to normalize data but whitened data are already normalized by dt
1075  // so before stft data must be divided by 1./sqrt(dt)
1076  for(int i=0;i<(int)y.size();i++) y.data[i]=pTF[n]->data[i+ystart]/sqrt(y.rate());
1077 
1078  CWB::STFT stft(y,nfft,noverlap,"energy","gauss",fparm);
1079 
1080  TCanvas* canvas;
1081  double tstart = nfft/pTF[n]->rate()+ystart/pTF[n]->rate();
1082  double tstop = (ysize-nfft)/pTF[n]->rate()+ystart/pTF[n]->rate();
1083  stft.Draw(tstart,tstop,pTF[n]->getlow(),pTF[n]->gethigh(),0,0,1);
1084  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle[n]);
1085  sprintf(fname, "%s/%s_spectrogram_1.%s", dirCED, NET->ifoName[n], gtype);
1086  canvas = stft.GetCanvas();
1087  if(sbasedirCED!=NULL) stft.Print(fname); else canvas->Write(REPLACE(fname,dirCED,gtype));
1088  canvas->SetLogy(true);
1089  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle[n]);
1090  sprintf(fname, "%s/%s_spectrogram_logy_1.%s", dirCED, NET->ifoName[n], gtype);
1091  if(sbasedirCED!=NULL) stft.Print(fname); else canvas->Write(REPLACE(fname,dirCED,gtype));
1092 
1093  tstart+=0.9;tstop-=0.9;
1094  stft.Draw(tstart,tstop,pTF[n]->getlow(),pTF[n]->gethigh(),0,0,1);
1095  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle[n]);
1096  sprintf(fname, "%s/%s_spectrogram_0.%s", dirCED, NET->ifoName[n], gtype);
1097  canvas = stft.GetCanvas();
1098  if(sbasedirCED!=NULL) stft.Print(fname); else canvas->Write(REPLACE(fname,dirCED,gtype));
1099  canvas->SetLogy(true);
1100  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle[n]);
1101  sprintf(fname, "%s/%s_spectrogram_logy_0.%s", dirCED, NET->ifoName[n], gtype);
1102  if(sbasedirCED!=NULL) stft.Print(fname); else canvas->Write(REPLACE(fname,dirCED,gtype));
1103 
1104  y.resize(0);
1105  }
1106 
1107  // set TF maps to optimal resolution level
1108  for(int n=0; n<nIFO; n++) {
1109  float fLow = pTF[n]->getlow();
1110  float fHigh = pTF[n]->gethigh();
1111  if(useSparse&&(analysis!="1G")) { // data not present in TF map (2G phase2), use sparse TP map
1112  int optSparse = pD[n]->getSTFind(optRate); // get optimal sparse map index
1113  delete pTF[n];
1114  pTF[n] = new WSeries<double>(*(WSeries<double>*)&(pD[n]->vSS[optSparse]));
1115  } else {
1116  if(pTF[n]->size()==0) continue;
1117  if(NET->wdm()) {if(NET->getwdm(optLayer+1)) pTF[n]->Forward(*pTF[n],*NET->getwdm(optLayer+1));}
1118  else pTF[n]->Forward(optLevel);
1119  }
1120  pTF[n]->setlow(fLow);
1121  pTF[n]->sethigh(fHigh);
1122  }
1123 
1124  WTS.canvas->cd();
1125 
1126  // plot scalogram maps at optimal resolution level
1127  for(int n=0; n<nIFO; n++) {
1128  if(pTF[n]->size()==0) continue;
1129  sprintf(fname, "%s/%s_scalogram_1.%s", dirCED, NET->ifoName[n], gtype);
1130  WTS.plot(pTF[n], 2, EVT->start[n]-EVT->slag[n]-1,
1131  EVT->stop[n]-EVT->slag[n]+1,const_cast<char*>("COLZ"));
1132  WTS.hist2D->GetYaxis()->SetRangeUser(pTF[n]->getlow(),pTF[n]->gethigh());
1133  WTS.hist2D->GetXaxis()->SetTitle(xtitle[n]);
1134  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1135  WTS.clear();
1136  }
1137 
1138  // plot zoomed scalogram maps at optimal resolution level
1139  for(int n=0; n<nIFO; n++) {
1140  if(pTF[n]->size()==0) continue;
1141  WTS.plot(pTF[n], 2, EVT->start[n]-EVT->slag[n]-0.1,
1142  EVT->stop[n]-EVT->slag[n]+0.1, const_cast<char*>("COLZ"));
1143  WTS.hist2D->GetYaxis()->SetRangeUser(pTF[n]->getlow(),pTF[n]->gethigh());
1144  WTS.hist2D->GetXaxis()->SetTitle(xtitle[n]);
1145  sprintf(fname, "%s/%s_scalogram_0.%s", dirCED, NET->ifoName[n], gtype);
1146  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1147  //sprintf(fname, "%s/%s_scalogram_0.%s", dirCED, NET->ifoName[n], "root");
1148  //WTS.canvas->Print(fname);
1149  WTS.clear();
1150  }
1151 
1152  // resize to 0 sseries
1153  for(int n=0; n<nIFO; n++) {
1154  if(useSparse&&(analysis!="1G")) { // data not present in TF map (2G phase2), the sparse TF maps are used
1155  int optSparse = pD[n]->getSTFind(optRate); // get optimal sparse map index
1156  pD[n]->vSS[optSparse].Shrink(); // resize 0 wseries : leave only sparse table
1157  delete pTF[n];
1158  } else {
1159  // rescale data when data are resampled (resample with wavelet change the amplitude)
1160  double rescale = 1./pow(sqrt(2.),TMath::Log2(inRate/pTF[n]->rate()));
1161  *pTF[n]*=1./rescale; // restore original rescaled TF map
1162  }
1163  }
1164 
1165  // get reconstructed wave forms, signal
1166 
1167  int type = 1;
1168 
1169  PTS.canvas->cd();
1170 
1171  if(NET->like()=='2') {NET->getMRAwave(ID,k,'S',wmod,true);NET->getMRAwave(ID,k,'W',wmod,true);}
1172  else NET->getwave(ID, k, 'W');
1173 
1174  // rescale data when data are resampled (resample with wavelet change the amplitude)
1175  double rescale = 1./pow(sqrt(2.),TMath::Log2(inRate/pD[0]->waveForm.rate()));
1176  for(int n=0; n<nIFO; n++) {
1177  pD[n]->waveForm*=rescale; pD[n]->waveBand*=rescale; // rescale data
1178  double bT,eT;
1179  GetBoundaries((wavearray<double>&)pD[n]->waveForm, 0.999, bT, eT);
1180  sprintf(fname, "%s/%s_wf_signal.%s", dirCED, NET->ifoName[n], gtype);
1181  if(NET->wdm()) {
1182  PTS.plot((wavearray<double>&)pD[n]->waveBand, const_cast<char*>("ALP"), 1, bT, eT);
1183  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("SAME"), 2, bT, eT);
1184  } else {
1185  PTS.plot((wavearray<double>&)pD[n]->waveBand, const_cast<char*>("ALP"), 1, bT, eT);
1186  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("SAME"), 2, bT, eT);
1187  }
1188  PTS.graph[0]->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1189  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1190  PTS.clear();
1191 
1192  sprintf(fname, "%s/%s_wf_signal_fft.%s", dirCED, NET->ifoName[n], gtype);
1193  double flow = EVT->low[0];
1194  double fhigh = EVT->high[0];
1195  if(NET->wdm()) {
1196  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("ALP"), 2, 0., 0., true, flow, fhigh);
1197  } else {
1198  PTS.plot((wavearray<double>&)pD[n]->waveBand, const_cast<char*>("ALP"), 1, 0., 0., true, flow, fhigh);
1199  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("SAME"), 2, 0., 0., true, flow, fhigh);
1200  }
1201  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1202  PTS.clear();
1203  }
1204 
1205  //save whitenet waveforms (SNR^2 = the sum of the squares of the samples)
1206 
1207  for(int n=0; n<nIFO; n++) {
1208  pD[n]->waveForm*=1./rescale; pD[n]->waveBand*=1./rescale; // restore original rescaled data
1209  //pD[n]->waveForm*=sqrt(pD[n]->waveForm.rate()); // normalize data to get snr=sqrt(sum(amp*amp*dt))
1210  sprintf(fname, "%s/%s_wf_signal.dat", dirCED, NET->ifoName[n]);
1211  double wstart = pD[n]->waveForm.wavearray<double>::start(); // save relative start time
1212  pD[n]->waveForm.wavearray<double>::start(EVT->gps[n]+pD[n]->waveForm.wavearray<double>::start()); // set absolute time
1213  if(sbasedirCED!=NULL) pD[n]->waveForm.wavearray<double>::Dump(const_cast<char*>(fname),2); // time, strain format
1214  //if(sbasedirCED!=NULL) pD[n]->waveForm.Dump(const_cast<char*>(fname));
1215  //pD[n]->waveForm*=1./sqrt(pD[n]->waveForm.rate()); // restore original rescaled data
1216  }
1217 
1218  // get reconstructed waveforms, signal + noise
1219 
1220  for(int n=0; n<nIFO; n++) pD[n]->waveBand.sethigh(0);
1221 
1222  if(NET->like()=='2') {NET->getMRAwave(ID,k,'s',wmod,true);NET->getMRAwave(ID,k,'w',wmod,true);}
1223  else NET->getwave(ID, k, 'S');
1224  for(int n=0; n<nIFO; n++) {
1225  pD[n]->waveForm*=rescale; pD[n]->waveBand*=rescale; // rescale data
1226  sprintf(fname, "%s/%s_wf_noise.%s", dirCED, NET->ifoName[n], gtype);
1227  if(NET->wdm()) {
1228  PTS.plot((wavearray<double>&)pD[n]->waveBand, const_cast<char*>("ALP"), 1);
1229  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("SAME"), 2);
1230  } else {
1231  PTS.plot((wavearray<double>&)pD[n]->waveBand, const_cast<char*>("ALP"), 1);
1232  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("SAME"), 2);
1233  }
1234  PTS.graph[0]->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1235  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1236  PTS.clear();
1237  }
1238 
1239  // get reconstructed waveforms, strain
1240  double bT,eT;
1241  for(int n=0; n<nIFO; n++) pD[n]->waveBand.sethigh(0);
1242  if(NET->like()=='2') NET->getMRAwave(ID,k,'s',wmod,true);
1243  else NET->getwave(ID, k, 'w');
1244  for(int n=0; n<nIFO; n++) {
1245  pD[n]->waveForm*=rescale; // rescale data
1246  GetBoundaries((wavearray<double>&)pD[n]->waveForm, 0.999, bT, eT);
1247  sprintf(fname, "%s/%s_wf_strain.%s", dirCED, NET->ifoName[n], gtype);
1248  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("ALP"), 2, bT, eT);
1249  PTS.graph[0]->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1250  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1251  PTS.clear();
1252 
1253  GetBoundaries((wavearray<double>&)pD[n]->waveForm, 0.95, bT, eT);
1254  sprintf(fname, "%s/%s_wf_strain_zoom.%s", dirCED, NET->ifoName[n], gtype);
1255  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("ALP"), 2, bT, eT);
1256  PTS.graph[0]->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1257  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1258  PTS.clear();
1259 
1260  double flow = EVT->low[0];
1261  double fhigh = EVT->high[0];
1262  sprintf(fname, "%s/%s_wf_strain_fft.%s", dirCED, NET->ifoName[n], gtype);
1263  PTS.plot((wavearray<double>&)pD[n]->waveForm, const_cast<char*>("ALP"), 2, 0., 0., true, flow, fhigh);
1264  if(sbasedirCED!=NULL) PTS.canvas->Print(fname); else PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1265  PTS.clear();
1266  }
1267 
1268  //save waveforms
1269 
1270  for(int n=0; n<nIFO; n++) {
1271  sprintf(fname, "%s/%s_wf_strain.dat", dirCED, NET->ifoName[n]);
1272  double wstart = pD[n]->waveForm.wavearray<double>::start(); // save relative start time
1273  pD[n]->waveForm.wavearray<double>::start(EVT->gps[n]+pD[n]->waveForm.wavearray<double>::start()); // set absolute time
1274  if(sbasedirCED!=NULL) pD[n]->waveForm.wavearray<double>::Dump(const_cast<char*>(fname),2); // time, strain format
1275  pD[n]->waveForm.wavearray<double>::start(wstart); // restore relative start time
1276  //if(sbasedirCED!=NULL) pD[n]->waveForm.Dump(const_cast<char*>(fname));
1277  }
1278 
1279  if(nIFO==1 || !fullCED) {
1280  if(NET->wdm()) NET->getwc(k)->sCuts[ID-1]=1; // mark clusters as processed (rejected)
1281  goto skip_skymap; // skip skymaps
1282  }
1283 
1284  // plot polargrams
1285 
1286  if(analysis.Contains("2G")) {
1287  for(m=0;m<2;m++) {
1288  gStyle->SetLineColor(kBlack);
1289  TCanvas* CPol = gNET.DrawPolargram(m,NET);
1290  if(CPol!=NULL) {
1291  sprintf(fname, "%s/polargram_%d.%s", dirCED, m+1, gtype);
1292  if(sbasedirCED!=NULL) CPol->Print(fname); else CPol->Write(REPLACE(fname,dirCED,gtype));
1293  }
1294  gStyle->SetLineColor(kWhite);
1295  }
1296  }
1297 
1298  //likelihood skymaps
1299 
1300  inj.SetX(EVT->phi[1]<180?EVT->phi[1]:EVT->phi[1]-360); inj.SetY(90.-EVT->theta[1]); // injected pos (white star)
1301  INJ.SetX(EVT->phi[1]<180?EVT->phi[1]:EVT->phi[1]-360); INJ.SetY(90.-EVT->theta[1]); // injected pos (ewhite star)
1302  rec.SetX(EVT->phi[0]<180?EVT->phi[0]:EVT->phi[0]-360); rec.SetY(90.-EVT->theta[0]); // recostr. pos (black star)
1303  REC.SetX(EVT->phi[0]<180?EVT->phi[0]:EVT->phi[0]-360); REC.SetY(90.-EVT->theta[0]); // recostr. pos (ewhite star)
1304  det.SetX(EVT->phi[3]<180?EVT->phi[3]:EVT->phi[3]-360); det.SetY(90.-EVT->theta[3]); // detected pos (black dot )
1305  DET.SetX(EVT->phi[3]<180?EVT->phi[3]:EVT->phi[3]-360); DET.SetY(90.-EVT->theta[3]); // detected pos (ewhite dot )
1306 
1307  inj.SetMarkerSize(2.5); inj.SetMarkerColor(kWhite);
1308  INJ.SetMarkerSize(2.5); INJ.SetMarkerColor(kBlack);
1309  rec.SetMarkerSize(2.5); rec.SetMarkerColor(kBlack);
1310  REC.SetMarkerSize(2.5); REC.SetMarkerColor(kWhite);
1311  det.SetMarkerSize(1.5); det.SetMarkerColor(kBlack);
1312  DET.SetMarkerSize(1.5); DET.SetMarkerColor(kWhite);
1313 
1314  gSM=NET->nSensitivity;
1315  sprintf(fname, "%s/sensitivity_plus.%s", dirCED, gtype);
1316  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1317  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1318  gSM=NET->nAlignment;
1319  sprintf(fname, "%s/sensitivity_cross.%s", dirCED, gtype);
1320  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1321  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1322  gSM=NET->nSkyStat;;
1323  sprintf(fname, "%s/skystat.%s", dirCED, gtype);
1324  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1325  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1326  gSM=NET->nLikelihood;;
1327  sprintf(fname, "%s/likelihood.%s", dirCED, gtype);
1328  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1329  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1330  gSM=NET->nNullEnergy;;
1331  sprintf(fname, "%s/null_energy.%s", dirCED, gtype);
1332  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1333  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1334  gSM=NET->nCorrEnergy;;
1335  sprintf(fname, "%s/corr_energy.%s", dirCED, gtype);
1336  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1337  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1338  gSM=NET->nPenalty;;
1339  sprintf(fname, "%s/penalty.%s", dirCED, gtype);
1340  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1341  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1342  gSM=NET->nDisbalance;;
1343  sprintf(fname, "%s/disbalance.%s", dirCED, gtype);
1344  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1345  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1346  gSM=NET->nCorrelation;;
1347  sprintf(fname, "%s/correlation.%s", dirCED, gtype);
1348  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1349  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1350  gSM=NET->nNetIndex;;
1351  gSM.GetHistogram()->GetZaxis()->SetRangeUser(1,nIFO);
1352  sprintf(fname, "%s/netindex.%s", dirCED, gtype);
1353  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1354  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1355  gSM.GetHistogram()->GetZaxis()->UnZoom();
1356  gSM=NET->nEllipticity;;
1357  sprintf(fname, "%s/ellipticity.%s", dirCED, gtype);
1358  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1359  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1360  gSM=NET->nPolarisation;;
1361  sprintf(fname, "%s/polarisation.%s", dirCED, gtype);
1362  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1363  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1364  gSM=NET->nProbability;;
1365  gSM.GetHistogram()->GetZaxis()->SetRangeUser(0,gSM.max());
1366  sprintf(fname, "%s/probability.%s", dirCED, gtype);
1367  gSM.Draw(); rec.Draw();REC.Draw(); det.Draw();DET.Draw(); if(simulation) {inj.Draw();INJ.Draw();}
1368  if(sbasedirCED!=NULL) gSM.Print(fname); else gSM.Write(REPLACE(fname,dirCED,gtype));
1369 
1370  // Dump probability skymap
1371  //int L = NET->nProbability.size();
1372  //wavearray<float> xprob(L);
1373  //for(int ll=0;ll<L;ll++) xprob.data[ll]=NET->nProbability.get(ll);
1374  sprintf(fname, "%s/probability.%s", dirCED, "root");
1375  if(sbasedirCED!=NULL) {
1376  gskymap gSM(NET->nProbability);
1377  gSM.SetOptions("","Geographic");
1378  gSM.DumpObject(fname);
1379  }
1380  //xprob.resize(0);
1381 
1382 #ifdef _USE_HEALPIX
1383  // Dump2fits probability skymap (healpix)
1384  if(sbasedirCED!=NULL) {
1385  sprintf(fname, "%s/probability.%s", dirCED, "fits");
1386  if(NET->nProbability.getType())
1387  NET->nProbability.Dump2fits(fname,EVT->time[0],const_cast<char*>(""),
1388  const_cast<char*>("PROB"),const_cast<char*>("pix-1"),'C');
1389  }
1390 #endif
1391 
1392  // add great circles (rec pos) to probability skymap
1393  sprintf(fname, "%s/probability_circles.%s", dirCED, gtype);
1394  gNET.SetGskymap(gSM);
1395  gNET.GetGskymap()->Draw();
1396  double phi,theta;
1397  CwbToGeographic(EVT->phi[0],EVT->theta[0],phi,theta);
1398  gNET.DrawCircles(phi,theta,(Color_t)kBlack,1,1,true);
1399  rec.Draw(); det.Draw();
1400  if(sbasedirCED!=NULL) gNET.GetGskymap()->Print(fname);
1401  else gNET.GetGskymap()->Write(REPLACE(fname,dirCED,gtype));
1402 
1403  skip_skymap:
1404 
1405  // get network time
1406  double gps_start = EVT->time[masterDet]-EVT->duration[1];
1407  double gps_stop = EVT->time[masterDet]+EVT->duration[1];
1408 
1409  // draw chirp : f^(-8/3) vs time
1410  if(NET->like()=='2') {
1411  PCH.canvas->cd();
1412  clusterdata* pcd = &(NET->getwc(k)->cData[ID-1]);
1413  sprintf(fname, "%s/mchirp.%s", dirCED, gtype);
1414  PCH.plot(pcd, simulation ? EVT->chirp[0] : 0);
1415  if(sbasedirCED!=NULL) PCH.canvas->Print(fname);
1416  else PCH.canvas->Write(REPLACE(fname,dirCED,gtype));
1417  PCH.clear();
1418  }
1419 
1420  // in likelihood2G pixeLHood,pixeLNull are not defined
1421  // use monster event display (multi resolution analysis)
1422  if(NET->like()=='2') {
1423  bool isPCs = !(NET->optim&&std::isupper(search)); // are Principal Components ?
1424  WTS.canvas->cd();
1425  netcluster* pwc = NET->getwc(k);
1426  sprintf(fname, "%s/l_tfmap_scalogram.%s", dirCED, gtype);
1427  WTS.plot(pwc, ID, nIFO, isPCs?'L':'l', 0, const_cast<char*>("COLZ"),256,NET->pattern>0);
1428  WTS.hist2D->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1429  if(sbasedirCED!=NULL) WTS.canvas->Print(fname);
1430  else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1431  WTS.clear();
1432  sprintf(fname, "%s/n_tfmap_scalogram.%s", dirCED, gtype);
1433  WTS.plot(pwc, ID, nIFO, isPCs?'N':'n', 0, const_cast<char*>("COLZ"),256,NET->pattern>0);
1434  WTS.hist2D->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1435  if(sbasedirCED!=NULL) WTS.canvas->Print(fname);
1436  else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1437  WTS.clear();
1438  } else {
1439  // plot likelihood map
1440  WTS.canvas->cd();
1441  sprintf(fname, "%s/l_tfmap_scalogram.%s", dirCED, gtype);
1442  WTS.plot(NET->pixeLHood, 0, gps_start-EVT->slag[masterDet],
1443  gps_stop-EVT->slag[masterDet], const_cast<char*>("COLZ"));
1444  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[0],EVT->high[0]);
1445  WTS.hist2D->SetTitle("Scalogram");
1446  WTS.hist2D->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1447  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1448  WTS.clear();
1449 
1450  // plot null map
1451  sprintf(fname, "%s/n_tfmap_scalogram.%s", dirCED, gtype);
1452  WTS.plot(NET->pixeLNull, 0, gps_start-EVT->slag[masterDet],
1453  gps_stop-EVT->slag[masterDet], const_cast<char*>("COLZ"));
1454  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[0],EVT->high[0]);
1455  WTS.hist2D->SetTitle("Scalogram");
1456  WTS.hist2D->GetXaxis()->SetTitle(xtitle[minTimeDet]);
1457  if(sbasedirCED!=NULL) WTS.canvas->Print(fname); else WTS.canvas->Write(REPLACE(fname,dirCED,gtype));
1458  WTS.clear();
1459  }
1460 
1461  // mark clusters as processed (rejected)
1462  if(NET->wdm()) NET->getwc(k)->sCuts[ID-1]=1;
1463 
1464  } // End loop on found events
1465  } // End loop on lags
1466 
1467  // restore NET thresholds
1468  NET->netCC = old_cc;
1469  NET->netRHO = old_rho;
1470 
1471  gROOT->SetBatch(batch); // restore batch status
1472 
1473  return 1;
1474 }
1475 
1476 double
1477 CWB::ced::GetBoundaries(wavearray<double> x, double P, double& bT, double& eT) {
1478 
1479  if(P<0) P=0;
1480  if(P>1) P=1;
1481 
1482  int N = x.size();
1483 
1484  double E = 0; // signal energy
1485  double avr = 0; // average
1486  for(int i=0;i<N;i++) {avr+=i*x[i]*x[i]; E+=x[i]*x[i];}
1487  int M=int(avr/E); // central index
1488 
1489  // search range which contains percentage P of the total energy E
1490  int jB=0;
1491  int jE=N-1;
1492  double a,b;
1493  double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
1494  for(int j=1; j<N; j++) {
1495  a = ((M-j>=0)&&(M-j<N)) ? x[M-j] : 0.;
1496  b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
1497  if(a) jB=M-j;
1498  if(b) jE=M+j;
1499  sum += a*a+b*b;
1500  if(sum/E > P) break;
1501  }
1502 
1503  bT = x.start()+jB/x.rate();
1504  eT = x.start()+jE/x.rate();
1505 
1506  return eT-bT;
1507 }
network::ifoName
std::vector< char * > ifoName
Definition: network.hh:591
gnetwork::DrawPolargram
TCanvas * DrawPolargram(int ptype, network *net=NULL)
Definition: gnetwork.cc:1516
WSeries::sethigh
void sethigh(double f)
Definition: wseries.hh:114
network::getifo
detector * getifo(size_t n)
param: detector index
Definition: network.hh:418
gskymap::GetHistogram
TH2D * GetHistogram()
Definition: gskymap.hh:120
gSM
gskymap * gSM
Definition: cwb_draw_antpat.C:49
Forward
tfmap Forward(ts, wdm)
rho
double rho
Definition: cwb_report_prod_2.C:149
network::netCC
double netCC
Definition: network.hh:578
WSeries::resize
virtual void resize(unsigned int)
Definition: wseries.cc:883
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
injection::time
Double_t * time
beam pattern coefficients for hx
Definition: injection.hh:70
NIFO
#define NIFO
Definition: wat.hh:56
wavecomplex::imag
double imag() const
Definition: wavecomplex.hh:52
noverlap
int noverlap
Definition: TestDelta.C:20
network::nLag
size_t nLag
Definition: network.hh:555
netcluster::cRate
std::vector< vector_int > cRate
Definition: netcluster.hh:380
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:93
xtitle
char xtitle[1024]
Definition: cwb_epparameters.C:182
netevent::high
Float_t * high
min frequency
Definition: netevent.hh:85
network::pixeLHood
WSeries< double > pixeLHood
Definition: network.hh:616
gNET
gnetwork * gNET
Definition: cwb_draw_antpat.C:47
wREC
std::vector< wavearray< double > > wREC[MAX_TRIALS]
Definition: CWB_Plugin_PE.C:406
network::wc_List
std::vector< netcluster > wc_List
Definition: network.hh:592
netevent::ecor
Float_t ecor
Definition: netevent.hh:101
netevent::start
Double_t * start
cluster duration = stopW-startW
Definition: netevent.hh:80
cmd
char cmd[1024]
Definition: cwb_condor_cleanup.C:20
network::getmdcTime
std::vector< double > * getmdcTime()
Definition: network.hh:404
fprintf
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
CWB
Definition: ced.hh:24
fHigh
double fHigh
Definition: cwb2G_parameters.C:20
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:123
netevent::duration
Float_t * duration
max cluster time relative to segment start
Definition: netevent.hh:79
netevent::low
Float_t * low
average center_of_snr frequency
Definition: netevent.hh:84
CWB::STFT::Print
void Print(TString pname)
Definition: STFT.cc:297
network::getwave
bool getwave(size_t, size_t, char='W')
param: cluster ID param: delay index param: time series type return: true if time series are extracte...
Definition: network.cc:3545
gnetwork::DrawCircles
void DrawCircles(double phi, double theta, double gps, Color_t lcolor=kBlack, Width_t lwidth=1, Style_t lstyle=1, bool labels=false)
Definition: gnetwork.cc:1289
a
wavearray< double > a(hp.size())
gnetwork::SetGskymap
void SetGskymap(gskymap &gSM)
Definition: gnetwork.hh:27
n
int n
Definition: cwb_net.C:10
network::nSkyStat
skymap nSkyStat
Definition: network.hh:610
const_cast< char * >
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
TString
TString("c")
Definition: cwb_report_skymap.C:111
gskymap::set
void set(size_t i, double a)
Definition: gskymap.hh:110
network::nRun
size_t nRun
Definition: network.hh:554
CWB::ced::GetBoundaries
double GetBoundaries(wavearray< double > x, double P, double &bT, double &eT)
Definition: ced.cc:1477
ID
int ID
Definition: TestMDC.C:70
netevent::right
Float_t * right
segment start GPS time
Definition: netevent.hh:77
network::nPenalty
skymap nPenalty
Definition: network.hh:606
leg
TLegend * leg
Definition: compare_bkg.C:246
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
netcluster::get
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2188
detector::RWFP
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:364
theta
float theta
Definition: cwb_condor_create_ced.C:51
nfact
int nfact
Definition: TestDelta.C:18
wdm
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
psd
wavearray< double > psd(33)
pwc
netcluster * pwc
Definition: cwb_job_obj.C:20
nfft
int nfft
Definition: TestDelta.C:19
ph
TH2F * ph
Definition: ReadHistFromJobFile.C:17
network::getmdcType
string getmdcType(size_t n)
Definition: network.hh:402
watplot::graph
std::vector< TGraph * > graph
Definition: watplot.hh:176
canvas
return wmap canvas
Definition: DrawClusterWithSkyplot.C:80
netevent::left
Float_t * left
min cluster time relative to segment start
Definition: netevent.hh:78
WSeries::setlow
void setlow(double f)
Definition: wseries.hh:107
P
TRandom3 P
Definition: compare_bkg.C:296
vSS
std::vector< SSeries< double > > vSS[NIFO_MAX]
Definition: CWB_Plugin_PE.C:433
html
THtml html
Definition: create_htmlpage.C:117
wavecomplex::real
double real() const
Definition: wavecomplex.hh:51
ifostr
char ifostr[64]
Definition: cwb_report_skymap.C:124
skymap::getTheta
double getTheta(size_t i)
Definition: skymap.hh:206
netevent::ioSNR
Float_t * ioSNR
reconstructed snr waveform
Definition: netevent.hh:120
network::nPolarisation
skymap nPolarisation
Definition: network.hh:612
size
Long_t size
Definition: cwb_condor_check.C:48
detector::waveBand
WSeries< double > waveBand
Definition: detector.hh:338
netevent::type
Int_t * type
event ID
Definition: netevent.hh:56
M
#define M
Definition: UniqSLagsList.C:3
m
int m
Definition: cwb_net.C:10
WaveDWT::pWWS
DataType_t * pWWS
Definition: WaveDWT.hh:123
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:119
skymap::max
double max()
Definition: skymap.cc:420
network::mdc__ID
std::vector< size_t > mdc__ID
Definition: network.hh:597
j
int j
Definition: cwb_net.C:10
detector::rate
double rate
Definition: detector.hh:323
network::getMRAwave
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
Definition: network.cc:3635
i
i drho i
Definition: cwb_epparameters.C:85
netcluster::rate
double rate
Definition: netcluster.hh:360
detector::tau
skymap tau
Definition: detector.hh:328
search
search
Definition: user_parameters.C:3
watplot::plot
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:132
gskymap::Draw
void Draw(int dpaletteId=0, Option_t *option="colfz")
Definition: gskymap.cc:442
N
#define N
GetBoundaries
double GetBoundaries(wavearray< double > x, double P, double &bT, double &eT)
Definition: CWB_Plugin_waveform.C:229
pi
double pi
Definition: TestChirp.C:18
network::nCorrEnergy
skymap nCorrEnergy
Definition: network.hh:607
netevent::dopen
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:291
watplot::hist2D
TH2F * hist2D
Definition: watplot.hh:175
CWB::STFT::Draw
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
Definition: STFT.cc:76
network::ifoListSize
size_t ifoListSize()
Definition: network.hh:413
netevent::run
Int_t run
max size used by allocate() for the probability maps
Definition: netevent.hh:53
network::getwdm
WDM< double > * getwdm(size_t M)
param: number of wdm layers
Definition: network.hh:430
network::nLikelihood
skymap nLikelihood
Definition: network.hh:604
watplot::clear
void clear()
Definition: watplot.hh:40
network::optim
bool optim
Definition: network.hh:572
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
skymap::getSkyIndex
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:702
tstart
double tstart
Definition: ReadCStrainFromJobFile.C:69
watplot::canvas
TCanvas * canvas
Definition: watplot.hh:174
Write
tlive_fix Write()
WSeries::getLevel
int getLevel()
Definition: wseries.hh:91
phi
float phi
Definition: cwb_condor_create_ced.C:51
netevent::frequency
Float_t * frequency
GPS stop time of the cluster.
Definition: netevent.hh:83
chName
TString chName[NIFO_MAX]
Definition: CWB_Plugin_Recolor.C:24
factor
double factor
Definition: ReadCStrainFromJobFile.C:18
ra
double ra
Definition: ConvertGWGC.C:46
detector::antenna
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
Definition: detector.cc:473
uname
TString uname
Definition: create_htmlpage.C:166
WSeries::getlow
double getlow() const
Definition: wseries.hh:111
netevent::output2G
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:687
netevent::nevent
Int_t nevent
run ID
Definition: netevent.hh:54
injection::hrss
Double_t * hrss
estimated bandwidth
Definition: injection.hh:76
detector::nRMS
WSeries< double > nRMS
Definition: detector.hh:340
WTS
watplot * WTS
Definition: ChirpMass.C:115
netevent::Deff
Float_t * Deff
eBBH array
Definition: netevent.hh:111
netevent::gps
Double_t * gps
average center_of_gravity time
Definition: netevent.hh:76
skymap::getRA
double getRA(size_t i)
Definition: skymap.hh:197
PCH
watplot * PCH
Definition: ChirpMass.C:116
int
i() int(T_cor *100))
NET
network NET
Definition: cwb_dump_inj.C:12
Pi
double Pi
Definition: DrawPhaseShift.C:12
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:4
fhigh
double fhigh
Definition: ReadCStrainFromJobFile.C:68
CWB::STFT::GetHistogram
TH2D * GetHistogram()
Definition: STFT.hh:53
netevent::stop
Double_t * stop
GPS start time of the cluster.
Definition: netevent.hh:81
pTF
WSeries< double > pTF[nRES]
Definition: revMonster.cc:8
netevent::lag
Float_t * lag
time between consecutive events
Definition: netevent.hh:65
fname
char fname[1024]
Definition: cwb_report_prod_1.C:110
segLen
segLen
Definition: cwb_eced.C:7
nsize
int nsize
Definition: cwb_report_prod_2.C:174
detector::RWFID
std::vector< int > RWFID
Definition: detector.hh:363
CWB::ced::Write
int Write(double factor, bool fullCED=true)
Definition: ced.cc:577
network::mdc__IDSize
size_t mdc__IDSize()
Definition: network.hh:396
k
int k
Definition: cwb_report_prod_2.C:150
network::pattern
int pattern
Definition: network.hh:583
tstop
double tstop
Definition: ReadCStrainFromJobFile.C:70
uinfo
UserGroup_t * uinfo
Definition: cwb_frdisplay.C:73
netevent::netcc
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:94
network::acor
double acor
Definition: network.hh:567
netevent::slag
Float_t * slag
time lag [sec]
Definition: netevent.hh:66
netevent::output
void output(TTree *=NULL, network *=NULL, double=0., size_t=0, int=-1)
Definition: netevent.cc:1166
skymap
Definition: skymap.hh:45
cwb_online.command
string command
Definition: cwb_online.py:382
e
double e
Definition: CreateListEBBH.C:35
network::likelihood
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
Definition: network.cc:4415
netcluster::size
size_t size()
Definition: netcluster.hh:129
netevent::hrss
Double_t * hrss
high-low
Definition: netevent.hh:87
ind
int * ind
Definition: cwb_fix_live_slag_missed.C:89
WSeries::gethigh
double gethigh() const
Definition: wseries.hh:118
network::like
char like()
Definition: network.hh:494
flow
double flow
Definition: ReadCStrainFromJobFile.C:67
gnetwork::GetGskymap
gskymap * GetGskymap()
Definition: gnetwork.hh:26
skymap::getPhi
double getPhi(size_t i)
Definition: skymap.hh:146
detector::TFmap
WSeries< double > TFmap
Definition: detector.hh:336
cwb_parameters_file
TString cwb_parameters_file
Definition: cwb_mkhtml_header.C:122
netevent::phi
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:68
network::pixeLNull
WSeries< double > pixeLNull
Definition: network.hh:617
network::nAlignment
skymap nAlignment
Definition: network.hh:602
EVT
netevent EVT(itree, nifo)
netevent::theta
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:69
netevent::time
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:75
T
double T
Definition: testWDM_4.C:11
netcluster::sCuts
std::vector< int > sCuts
Definition: netcluster.hh:374
detector::waveForm
WSeries< double > waveForm
Definition: detector.hh:337
netevent::dclose
void dclose()
Definition: netevent.hh:300
fLow
double fLow
Definition: cwb2G_parameters.C:19
detector::getTFmap
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:161
skymap::getType
int getType()
Definition: skymap.hh:278
network::nDisbalance
skymap nDisbalance
Definition: network.hh:609
Meyer
Definition: Meyer.hh:18
fabs
double fabs(const Complex &x)
Definition: numpy.cc:37
WSeries::resolution
double resolution(int=0)
Definition: wseries.hh:137
CWB::STFT::GetCanvas
TCanvas * GetCanvas()
Definition: STFT.hh:52
WSeries::Forward
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:228
CwbToGeographic
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:396
strcpy
strcpy(RunLabel, RUN_LABEL)
gskymap::DumpObject
void DumpObject(const char *file, const char *name="gskymap")
Definition: gskymap.cc:1213
netevent::likelihood
Float_t likelihood
Definition: netevent.hh:105
S
Meyer< double > S(1024, 2)
network::getwc
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:421
netevent::bx
Float_t * bx
beam pattern coefficients for hp
Definition: netevent.hh:73
l
int l
Definition: cbc_plots.C:434
network::nNetIndex
skymap nNetIndex
Definition: network.hh:608
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
network::nSensitivity
skymap nSensitivity
list of wdm tranformations
Definition: network.hh:601
netcluster::cData
std::vector< clusterdata > cData
Definition: netcluster.hh:373
skymap::set
void set(size_t i, double a)
param: sky index param: value to set
Definition: skymap.hh:104
cwb_uparameters_file
TString cwb_uparameters_file
Definition: cwb_mkhtml_header.C:128
wavearray::data
DataType_t * data
Definition: wavearray.hh:301
K
int K
Definition: CWB_Plugin_Sim4_BRST_Config.C:87
id
Long_t id
Definition: cwb_condor_check.C:48
CWB::ced
Definition: ced.hh:26
t1
TH1 * t1
Definition: cbc_plots_sim4_ifar.C:1896
network::nNullEnergy
skymap nNullEnergy
Definition: network.hh:605
work_dir
char work_dir[512]
Definition: test_config1.C:143
skymap::get
double get(size_t i)
param: sky index
Definition: skymap.cc:681
WSeries::pWavelet
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:438
wINJ
wavearray< double > wINJ[NIFO_MAX]
Definition: CWB_Plugin_xWRC.C:73
close
in close()
rnID
int rnID
Definition: create_htmlpage.C:164
network::nCorrelation
skymap nCorrelation
Definition: network.hh:603
REPLACE
#define REPLACE(STRING, DIR, EXT)
Definition: ced.cc:27
getTFmap
cout<< "live time after cat 2 : "<< detSegs_ctime<< endl;if(detSegs_ctime< segTHR){cout<< "job segment live time after cat2 < "<< segTHR<< " sec, job terminated !!!"<< endl;exit(1);}double Tb=detSegs[0].start;double Te=detSegs[0].stop;double dT=Te-Tb;char file[512], tdf00[512], tdf90[512], buFFer[1024];int rnID=int(gRandom->Rndm(13)*1.e9);if(simulation){i=NET.readMDClog(injectionList, double(long(Tb))-mdcShift);printf("GPS: %16.6f saved, injections: %d\n", double(long(Tb)), i);frTB[nIFO].shiftBurstMDCLog(NET.mdcList, ifos, mdcShift);for(int i=0;i< NET.mdcTime.size();i++) NET.mdcTime[i]+=mdcShift;vector< waveSegment > mdcSegs(NET.mdcTime.size());for(int k=0;k< NET.mdcTime.size();k++){mdcSegs[k].start=NET.mdcTime[k]-gap;mdcSegs[k].stop=NET.mdcTime[k]+gap;}vector< waveSegment > mdcSegs_dq2=slagTB.mergeSegLists(detSegs_dq2, mdcSegs);double mdcSegs_ctime=slagTB.getTimeSegList(mdcSegs_dq2);cout<< "live time in zero lag after cat2+inj : "<< mdcSegs_ctime<< endl;if(mdcSegs_ctime==0){cout<< "job segment with zero cat2+inj live time in zero lag, job terminated !!!"<< endl;exit(1);}}if(dump_infos_and_exit) exit(0);if(mask >0.) NET.setSkyMask(mask, skyMaskFile);for(i=0;i< nIFO;i++){frTB[i].readFrames(FRF[i], channelNamesRaw[i], x);x.start(x.start()+dataShift[i]);x.start(x.start()-segLen *(segID[i]-segID[0]));if(singleDetector) TB.resampleToPowerOfTwo(x);sprintf(file,"%s/%s_%d_%s_%d_%d.dat", nodedir, ifo[i], int(Tb), data_label, runID, rnID);if(dump_sensitivity_and_exit){sprintf(file,"%s/sensitivity_%s_%d_%s_job%d.txt", dump_dir, ifo[i], int(Tb), data_label, runID);cout<< endl<< "Dump Sensitivity : "<< file<< endl<< endl;TB.makeSpectrum(file, x);continue;}if(dcCal[i]>0.) x *=dcCal[i];if(fResample >0){x.FFT(1);x.resize(fResample/x.rate()*x.size());x.FFT(-1);x.rate(fResample);}pTF[i]=pD[i]-> getTFmap()
injection::fill_in
Bool_t fill_in(network *, int, bool=true)
Definition: injection.cc:325
netevent::neted
Float_t * neted
network correlation coefficients: 0-net,1-pc,2-cc,3-net2
Definition: netevent.hh:95
simulation
simulation
Definition: cwb_eced.C:9
fclose
fclose(ftrig)
detector::vSS
std::vector< SSeries< double > > vSS
Definition: detector.hh:334
skymap::size
size_t size()
Definition: skymap.hh:118
network::getmdc__ID
size_t getmdc__ID(size_t n)
Definition: network.hh:408
for
for(int i=0;i< 101;++i) Cos2[2][i]=0
Definition: CWB_Plugin_MDC_OTF_Config_BRST.C:36
injection::pwf
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
Definition: injection.hh:82
wavearray::resize
virtual void resize(unsigned int)
Definition: wavearray.cc:445
gskymap::Print
void Print(TString pname)
Definition: gskymap.cc:1104
network::MRA
bool MRA
Definition: network.hh:581
WSeries::Inverse
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:273
y
wavearray< double > y
Definition: Test10.C:31
network::wdm
bool wdm()
Definition: network.hh:492
detector::getSTFind
int getSTFind(double r)
Definition: detector.hh:164
network::tYPe
char tYPe
Definition: network.hh:570
network::netRHO
double netRHO
Definition: network.hh:579
network::setndm
bool setndm(size_t, size_t, bool=true, int=1)
param: cluster ID param: lag index param: statistic identificator param: resolution idenificator retu...
Definition: network.cc:6494
gskymap::SetOptions
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:66
network::nProbability
skymap nProbability
Definition: network.hh:613
network::SETNDM
bool SETNDM(size_t, size_t, bool=true, int=1)
Definition: network.cc:6793
fparm
double fparm
Definition: TestDelta.C:22
netevent::chirp
Float_t * chirp
range to source: [0/1]-rec/inj
Definition: netevent.hh:109
pD
detector ** pD
Definition: cwb_dump_events.C:34
injection::Deff
Float_t * Deff
injected snr in the detectors
Definition: injection.hh:78
stft
CWB::STFT * stft
Definition: ChirpMass.C:117
gskymap::SetZaxisTitle
void SetZaxisTitle(TString zAxisTitle)
Definition: gskymap.hh:132
exit
exit(0)
netevent::bp
Float_t * bp
[0]-reconstructed iota angle, [1]-injected iota angle
Definition: netevent.hh:72
netevent::iSNR
Float_t * iSNR
energy of reconstructed responses Sk*Sk
Definition: netevent.hh:118
avr
double avr
Definition: TestNormalization.C:46
netevent::bandwidth
Float_t * bandwidth
max frequency
Definition: netevent.hh:86
network::nEllipticity
skymap nEllipticity
Definition: network.hh:611
netevent::oSNR
Float_t * oSNR
injected snr waveform
Definition: netevent.hh:119