Logo coherent WaveBurst  
Library Reference Guide
Logo
netevent.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 // netevent class to process and store cWB triggers in root file
20 // S.Klimenko, University of Florida, Gainesville, FL
21 // G.Vedovato, INFN, Sezione di Padova, Italy
22 
23 #include "netevent.hh"
24 #include "wseries.hh"
25 
26 #include "TH2.h"
27 #include "TStyle.h"
28 #include "TCanvas.h"
29 #include "TSystem.h"
30 #include "TMath.h"
31 #include "TMarker.h"
32 #include "TVector3.h"
33 #include "TRotation.h"
34 #include "TPolyLine.h"
35 #include "Math/Rotation3D.h"
36 #include "Math/Vector3Dfwd.h"
37 
38 #include "Meyer.hh"
39 #include <string.h>
40 
41 #define WAVE_TREE_NAME "waveburst"
42 
43 ClassImp(netevent) // used by THtml doc
44 
45 using namespace ROOT::Math;
46 
48 {
49  iFile = TFile::Open(fName);
50  if((iFile==NULL) || (iFile!=NULL && !iFile->IsOpen())) {
51  cout << "netevent::Init : Error opening root file " << fName.Data() << endl;
52  exit(1);
53  }
54 
55  TTree* tree = (TTree *) iFile->Get(WAVE_TREE_NAME);
56  if(tree) {
57  ndim = tree->GetUserInfo()->GetSize(); // get number of detectors
58  if(ndim>0) {
59  if(n>0 && ndim!=n) {
60  cout << "netevent::Init : number of detectors declared in the constructor (" << n
61  << ") are not equals to the one ("<<ndim<<") declared in the root file : "
62  << fName.Data() << endl;
63  exit(1);
64  }
65  } else ndim=n;
66  } else {
67  cout << "netevent::Init : object tree " << WAVE_TREE_NAME
68  << " not present in the root file " << fName.Data() << endl;
69  exit(1);
70  }
71 
72 
73  if(ndim==0) {
74  cout << "netevent::Init : detector number is not declared in the constructor or"
75  << " not present in the root file " << fName.Data() << endl;
76  exit(1);
77  }
78 
79 
80  return tree;
81 }
82 
83 // Set branch addresses
84 void netevent::Init(TTree *tree)
85 {
86  if (tree == 0) return;
87  fChain = tree;
88  fCurrent = -1;
89  fChain->SetMakeClass(1);
90 
91  fChain->SetBranchAddress("ndim",&ndim);
92  fChain->SetBranchAddress("run",&run);
93  fChain->SetBranchAddress("nevent",&nevent);
94  fChain->SetBranchAddress("eventID",eventID);
95  fChain->SetBranchAddress("type",type);
96  fChain->SetBranchAddress("name",&name);
97  fChain->SetBranchAddress("log",&log);
98  fChain->SetBranchAddress("rate",rate);
99 
100  fChain->SetBranchAddress("volume",volume);
101  fChain->SetBranchAddress("size",size);
102  fChain->SetBranchAddress("usize",&usize);
103 
104  fChain->SetBranchAddress("gap",gap);
105  fChain->SetBranchAddress("lag",lag);
106  fChain->SetBranchAddress("slag",slag);
107  fChain->SetBranchAddress("strain",strain);
108  fChain->SetBranchAddress("phi",phi);
109  fChain->SetBranchAddress("theta",theta);
110  fChain->SetBranchAddress("psi",psi);
111  fChain->SetBranchAddress("iota",iota);
112  fChain->SetBranchAddress("bp",bp);
113  fChain->SetBranchAddress("bx",bx);
114 
115 
116  fChain->SetBranchAddress("time",time);
117  fChain->SetBranchAddress("gps",gps);
118  fChain->SetBranchAddress("right",right);
119  fChain->SetBranchAddress("left",left);
120  fChain->SetBranchAddress("duration",duration);
121  fChain->SetBranchAddress("start",start);
122  fChain->SetBranchAddress("stop",stop);
123 
124  fChain->SetBranchAddress("frequency",frequency);
125  fChain->SetBranchAddress("low",low);
126  fChain->SetBranchAddress("high",high);
127  fChain->SetBranchAddress("bandwidth",bandwidth);
128 
129  fChain->SetBranchAddress("hrss",hrss);
130  fChain->SetBranchAddress("noise",noise);
131  fChain->SetBranchAddress("erA",erA);
132  if(fChain->GetBranch("Psm")!=NULL) fChain->SetBranchAddress("Psm",&Psm);
133  fChain->SetBranchAddress("null",null);
134  fChain->SetBranchAddress("nill",nill);
135  fChain->SetBranchAddress("netcc",netcc);
136  fChain->SetBranchAddress("neted",neted);
137  fChain->SetBranchAddress("rho",rho);
138 
139  fChain->SetBranchAddress("gnet",&gnet);
140  fChain->SetBranchAddress("anet",&anet);
141  fChain->SetBranchAddress("inet",&inet);
142  fChain->SetBranchAddress("ecor",&ecor);
143  fChain->SetBranchAddress("norm",&norm);
144  fChain->SetBranchAddress("ECOR",&ECOR);
145  fChain->SetBranchAddress("penalty",&penalty);
146  fChain->SetBranchAddress("likelihood",&likelihood);
147 
148  fChain->SetBranchAddress("factor",&factor);
149  fChain->SetBranchAddress("range",range);
150  fChain->SetBranchAddress("chirp",chirp);
151  fChain->SetBranchAddress("eBBH",eBBH);
152  fChain->SetBranchAddress("Deff",Deff);
153  fChain->SetBranchAddress("mass",mass);
154  fChain->SetBranchAddress("spin",spin);
155 
156  fChain->SetBranchAddress("snr",snr);
157  fChain->SetBranchAddress("xSNR",xSNR);
158  fChain->SetBranchAddress("sSNR",sSNR);
159  fChain->SetBranchAddress("iSNR",iSNR);
160  fChain->SetBranchAddress("oSNR",oSNR);
161  fChain->SetBranchAddress("ioSNR",ioSNR);
162 
163  Notify();
164 }
165 
166 // allocate memory
168 {
169  if(ndim==0) {
170  cout << "netevent::allocate : Error - number of detectors must be > 0" << endl;
171  exit(1);
172  }
173 
174  if (!eventID) eventID= (Int_t*)malloc(2*sizeof(Int_t));
175  else eventID= (Int_t*)realloc(eventID,2*sizeof(Int_t));
176  if (!type) type= (Int_t*)malloc(2*sizeof(Int_t));
177  else type= (Int_t*)realloc(type,2*sizeof(Int_t));
178  if (!name) name= new string();
179  else {delete name;name= new string();}
180  if (!log) log= new string();
181  else {delete log; log= new string();}
182  if (!rate) rate= (Int_t*)malloc(ndim*sizeof(Int_t));
183  else rate= (Int_t*)realloc(rate,ndim*sizeof(Int_t));
184 
185  if (!volume) volume= (Int_t*)malloc(ndim*sizeof(Int_t));
186  else volume= (Int_t*)realloc(volume,ndim*sizeof(Int_t));
187  if (!size) size= (Int_t*)malloc(ndim*sizeof(Int_t));
188  else size= (Int_t*)realloc(size,ndim*sizeof(Int_t));
189 
190  if (!gap) gap= (Float_t*)malloc(ndim*sizeof(Float_t));
191  else gap= (Float_t*)realloc(gap,ndim*sizeof(Float_t));
192  if (!lag) lag= (Float_t*)malloc((ndim+1)*sizeof(Float_t));
193  else lag= (Float_t*)realloc(lag,(ndim+1)*sizeof(Float_t));
194  if (!slag) slag= (Float_t*)malloc((ndim+1)*sizeof(Float_t));
195  else slag= (Float_t*)realloc(slag,(ndim+1)*sizeof(Float_t));
196  if (!gps) gps= (Double_t*)malloc((ndim)*sizeof(Double_t));
197  else gps= (Double_t*)realloc(gps,(ndim)*sizeof(Double_t));
198  if (!strain) strain= (Double_t*)malloc(2*sizeof(Double_t));
199  else strain= (Double_t*)realloc(strain,2*sizeof(Double_t));
200  if (!phi) phi= (Float_t*)malloc(4*sizeof(Float_t));
201  else phi= (Float_t*)realloc(phi,4*sizeof(Float_t));
202  if (!theta) theta= (Float_t*)malloc(4*sizeof(Float_t));
203  else theta= (Float_t*)realloc(theta,4*sizeof(Float_t));
204  if (!psi) psi= (Float_t*)malloc(2*sizeof(Float_t));
205  else psi= (Float_t*)realloc(psi,2*sizeof(Float_t));
206  if (!iota) iota= (Float_t*)malloc(2*sizeof(Float_t));
207  else iota= (Float_t*)realloc(iota,2*sizeof(Float_t));
208  if (!bp) bp= (Float_t*)malloc(ndim*2*sizeof(Float_t));
209  else bp= (Float_t*)realloc(bp,ndim*2*sizeof(Float_t));
210  if (!bx) bx= (Float_t*)malloc(ndim*2*sizeof(Float_t));
211  else bx= (Float_t*)realloc(bx,ndim*2*sizeof(Float_t));
212 
213  if (!time) time= (Double_t*)malloc(ndim*2*sizeof(Double_t));
214  else time= (Double_t*)realloc(time,ndim*2*sizeof(Double_t));
215  if (!right) right= (Float_t*)malloc(ndim*sizeof(Float_t));
216  else right= (Float_t*)realloc(right,ndim*sizeof(Float_t));
217  if (!left) left= (Float_t*)malloc(ndim*sizeof(Float_t));
218  else left= (Float_t*)realloc(left,ndim*sizeof(Float_t));
219  if (!duration) duration= (Float_t*)malloc(ndim*sizeof(Float_t));
220  else duration= (Float_t*)realloc(duration,ndim*sizeof(Float_t));
221  if (!start) start= (Double_t*)malloc(ndim*sizeof(Double_t));
222  else start= (Double_t*)realloc(start,ndim*sizeof(Double_t));
223  if (!stop) stop= (Double_t*)malloc(ndim*sizeof(Double_t));
224  else stop= (Double_t*)realloc(stop,ndim*sizeof(Double_t));
225 
226  if (!frequency) frequency=(Float_t*)malloc(ndim*sizeof(Float_t));
227  else frequency=(Float_t*)realloc(frequency,ndim*sizeof(Float_t));
228  if (!low) low= (Float_t*)malloc(ndim*sizeof(Float_t));
229  else low= (Float_t*)realloc(low,ndim*sizeof(Float_t));
230  if (!high) high= (Float_t*)malloc(ndim*sizeof(Float_t));
231  else high= (Float_t*)realloc(high,ndim*sizeof(Float_t));
232  if (!bandwidth) bandwidth=(Float_t*)malloc(ndim*sizeof(Float_t));
233  else bandwidth=(Float_t*)realloc(bandwidth,ndim*sizeof(Float_t));
234  if (!hrss) hrss= (Double_t*)malloc(ndim*2*sizeof(Double_t));
235  else hrss= (Double_t*)realloc(hrss,ndim*2*sizeof(Double_t));
236  if (!noise) noise= (Double_t*)malloc(ndim*sizeof(Double_t));
237  else noise= (Double_t*)realloc(noise,ndim*sizeof(Double_t));
238  if (!null) null= (Float_t*)malloc(ndim*sizeof(Float_t));
239  else null= (Float_t*)realloc(null,ndim*sizeof(Float_t));
240  if (!nill) nill= (Float_t*)malloc(ndim*sizeof(Float_t));
241  else nill= (Float_t*)realloc(nill,ndim*sizeof(Float_t));
242  if (!netcc) netcc= (Float_t*)malloc(4*sizeof(Float_t));
243  else netcc= (Float_t*)realloc(netcc,4*sizeof(Float_t));
244  if (!neted) neted= (Float_t*)malloc(5*sizeof(Float_t));
245  else neted= (Float_t*)realloc(neted,5*sizeof(Float_t));
246  if (!rho) rho= (Float_t*)malloc(2*sizeof(Float_t));
247  else rho= (Float_t*)realloc(rho,2*sizeof(Float_t));
248  if (!erA) erA= (Float_t*)malloc(11*sizeof(Float_t));
249  else erA= (Float_t*)realloc(erA,11*sizeof(Float_t));
250  if (!range) range= (Float_t*)malloc(2*sizeof(Float_t));
251  else range= (Float_t*)realloc(range,2*sizeof(Float_t));
252  if (!chirp) chirp= (Float_t*)malloc(6*sizeof(Float_t));
253  else chirp= (Float_t*)realloc(chirp,6*sizeof(Float_t));
254  if (!eBBH) eBBH= (Float_t*)malloc(4*sizeof(Float_t));
255  else eBBH= (Float_t*)realloc(eBBH,4*sizeof(Float_t));
256  if (!Deff) Deff= (Float_t*)malloc(ndim*sizeof(Float_t));
257  else Deff= (Float_t*)realloc(Deff,ndim*sizeof(Float_t));
258  if (!mass) mass= (Float_t*)malloc(2*sizeof(Float_t));
259  else mass= (Float_t*)realloc(mass,ndim*sizeof(Float_t));
260  if (!spin) spin= (Float_t*)malloc(6*sizeof(Float_t));
261  else spin= (Float_t*)realloc(spin,6*sizeof(Float_t));
262 
263  if (!snr) snr= (Float_t*)malloc(ndim*sizeof(Float_t));
264  else snr= (Float_t*)realloc(snr,ndim*sizeof(Float_t));
265  if (!xSNR) xSNR= (Float_t*)malloc(ndim*sizeof(Float_t));
266  else xSNR= (Float_t*)realloc(xSNR,ndim*sizeof(Float_t));
267  if (!sSNR) sSNR= (Float_t*)malloc(ndim*sizeof(Float_t));
268  else sSNR= (Float_t*)realloc(sSNR,ndim*sizeof(Float_t));
269  if (!iSNR) iSNR= (Float_t*)malloc(ndim*sizeof(Float_t));
270  else iSNR= (Float_t*)realloc(iSNR,ndim*sizeof(Float_t));
271  if (!oSNR) oSNR= (Float_t*)malloc(ndim*sizeof(Float_t));
272  else oSNR= (Float_t*)realloc(oSNR,ndim*sizeof(Float_t));
273  if (!ioSNR) ioSNR= (Float_t*)malloc(ndim*sizeof(Float_t));
274  else ioSNR= (Float_t*)realloc(ioSNR,ndim*sizeof(Float_t));
275 
276  if (!Psm) Psm= WAT::USE_HEALPIX() ? new skymap(int(0)) : new skymap(0.);
277  else {delete Psm; Psm= WAT::USE_HEALPIX() ? new skymap(int(0)) : new skymap(0.);}
278 
279  return;
280 }
281 
282 // init array
284 {
285  for(int i=0;i<2;i++) eventID[i]=0;
286  for(int i=0;i<2;i++) type[i]=0;
287  for(int i=0;i<2;i++) strain[i]=0;
288  for(int i=0;i<4;i++) phi[i]=0;
289  for(int i=0;i<2;i++) psi[i]=0;
290  for(int i=0;i<4;i++) theta[i]=0;
291  for(int i=0;i<2;i++) psi[i]=0;
292  for(int i=0;i<2;i++) iota[i]=0;
293  for(int i=0;i<4;i++) netcc[i]=0;
294  for(int i=0;i<5;i++) neted[i]=0;
295  for(int i=0;i<2;i++) rho[i]=0;
296  for(int i=0;i<11;i++) erA[i]=0;
297  for(int i=0;i<2;i++) mass[i]=0;
298  for(int i=0;i<2;i++) range[i]=0;
299  for(int i=0;i<6;i++) chirp[i]=0;
300  for(int i=0;i<4;i++) eBBH[i]=0;
301  for(int i=0;i<6;i++) spin[i]=0;
302 
303  if(Psm) *Psm=0;
304 
305  for(int i=0;i<ndim+1;i++) lag[i]=0;
306  for(int i=0;i<ndim+1;i++) slag[i]=0;
307 
308  for(int i=0;i<2*ndim;i++) bp[i]=0;
309  for(int i=0;i<2*ndim;i++) bx[i]=0;
310  for(int i=0;i<2*ndim;i++) hrss[i]=0;
311  for(int i=0;i<2*ndim;i++) time[i]=0;
312 
313  for(int i=0;i<ndim;i++) {
314  rate[i]=0; volume[i]=0; size[i]=0; gap[i]=0;
315  gps[i]=0; snr[i]=0; right[i]=0; left[i]=0;
316  duration[i]=0; start[i]=0; stop[i]=0; frequency[i]=0; low[i]=0;
317  high[i]=0; bandwidth[i]=0; noise[i]=0; null[i]=0; nill[i]=0;
318  sSNR[i]=0; Deff[i]=0; iSNR[i]=0; oSNR[i]=0; ioSNR[i]=0; xSNR[i]=0;
319  }
320 
321  return;
322 }
323 
324 
326 {
327  // Called when loading a new file.
328  // Get branch pointers.
329  b_ndim = fChain->GetBranch("ndim");
330  b_run = fChain->GetBranch("run");
331  b_nevent = fChain->GetBranch("nevent");
332  b_eventID = fChain->GetBranch("eventID");
333  b_type = fChain->GetBranch("type");
334  b_name = fChain->GetBranch("name");
335  b_log = fChain->GetBranch("log");
336  b_rate = fChain->GetBranch("rate");
337 
338  b_volume = fChain->GetBranch("volume");
339  b_size = fChain->GetBranch("size");
340  b_usize = fChain->GetBranch("usize");
341 
342  b_gap = fChain->GetBranch("gap");
343  b_lag = fChain->GetBranch("lag");
344  b_slag = fChain->GetBranch("slag");
345  b_strain = fChain->GetBranch("strain");
346  b_phi = fChain->GetBranch("phi");
347  b_theta = fChain->GetBranch("theta");
348  b_psi = fChain->GetBranch("psi");
349  b_iota = fChain->GetBranch("iota");
350  b_bp= fChain->GetBranch("bx");
351  b_bx = fChain->GetBranch("bp");
352 
353  b_time = fChain->GetBranch("time");
354  b_gps = fChain->GetBranch("gps");
355  b_right = fChain->GetBranch("right");
356  b_left = fChain->GetBranch("left");
357  b_start = fChain->GetBranch("start");
358  b_stop = fChain->GetBranch("stop");
359  b_duration = fChain->GetBranch("duration");
360 
361  b_frequency = fChain->GetBranch("frequency");
362  b_low = fChain->GetBranch("low");
363  b_high = fChain->GetBranch("high");
364  b_bandwidth = fChain->GetBranch("bandwidth");
365 
366  b_hrss = fChain->GetBranch("hrss");
367  b_noise = fChain->GetBranch("noise");
368  b_erA = fChain->GetBranch("erA");
369  b_Psm = fChain->GetBranch("Psm");
370  b_null = fChain->GetBranch("null");
371  b_nill = fChain->GetBranch("nill");
372  b_netcc = fChain->GetBranch("netcc");
373  b_neted = fChain->GetBranch("neted");
374  b_rho = fChain->GetBranch("rho");
375 
376  b_gnet = fChain->GetBranch("gnet");
377  b_anet = fChain->GetBranch("anet");
378  b_inet = fChain->GetBranch("inet");
379  b_ecor = fChain->GetBranch("ecor");
380  b_norm = fChain->GetBranch("norm");
381  b_ECOR = fChain->GetBranch("ECOR");
382  b_penalty = fChain->GetBranch("penalty");
383  b_likelihood = fChain->GetBranch("likelihood");
384 
385  b_factor = fChain->GetBranch("factor");
386  b_range = fChain->GetBranch("range");
387  b_chirp = fChain->GetBranch("chirp");
388  b_eBBH = fChain->GetBranch("eBBH");
389  b_Deff = fChain->GetBranch("Deff");
390  b_mass = fChain->GetBranch("mass");
391  b_spin = fChain->GetBranch("spin");
392 
393  b_snr = fChain->GetBranch("snr");
394  b_xSNR = fChain->GetBranch("xSNR");
395  b_sSNR = fChain->GetBranch("sSNR");
396  b_iSNR = fChain->GetBranch("iSNR");
397  b_oSNR = fChain->GetBranch("oSNR");
398  b_ioSNR = fChain->GetBranch("ioSNR");
399 
400  return kTRUE;
401 }
402 
404 {
405  if (!fChain) return 0;
406  return fChain->GetEntries();
407 };
408 
410 {
411  if (!fChain) return 0;
412  return fChain->GetEntry(entry);
413 };
414 
416 {
417 // Print contents of entry.
418 // If entry is not specified, print current entry
419  if (!fChain) return;
420  fChain->Show(entry);
421 }
422 
423 //++++++++++++++++++++++++++++++++++++++++++++++
424 // set super lags
425 //++++++++++++++++++++++++++++++++++++++++++++++
427 {
428  for(int n=0;n<=ndim;n++) this->slag[n]=slag[n];
429 }
430 
431 //++++++++++++++++++++++++++++++++++++++++++++++
432 // set single event tree
433 //++++++++++++++++++++++++++++++++++++++++++++++
435 {
436  TTree* waveTree = new TTree(WAVE_TREE_NAME,WAVE_TREE_NAME);
437 
438  int i;
439 
440  char crate[16];
441 
442  char cvolume[16];
443  char csize[16];
444 
445  char cgap[16];
446  char clag[16];
447  char cslag[16];
448  char cgps[16];
449  char cbp[16];
450  char cbx[16];
451 
452  char ctime[16];
453  char cright[16];
454  char cleft[16];
455  char cduration[16];
456  char cstart[16];
457  char cstop[16];
458 
459  char cfrequency[16];
460  char clow[16];
461  char chigh[16];
462  char cbandwidth[16];
463  char chrss[16];
464  char cnoise[16];
465  char cnull[16];
466  char cnill[16];
467  char crange[16];
468  char cchirp[16];
469  char ceBBH[16];
470  char cDeff[16];
471  char cmass[16];
472  char cspin[16];
473  char csnr[16];
474  char csSNR[16];
475  char cxSNR[16];
476  char ciSNR[16];
477  char coSNR[16];
478  char cioSNR[16];
479 
480 
481  sprintf(crate, "rate[%1d]/I",ndim);
482 
483  sprintf(cvolume, "volume[%1d]/I",ndim);
484  sprintf(csize, "size[%1d]/I",ndim);
485 
486  sprintf(cgap, "gap[%1d]/F",ndim);
487  sprintf(clag, "lag[%1d]/F",ndim+1);
488  sprintf(cslag, "slag[%1d]/F",ndim+1);
489  sprintf(cgps, "gps[%1d]/D",ndim);
490 
491  if(ndim<5) sprintf(cbp, "bp[%1d]/F",ndim*2);
492  else sprintf(cbp, "bp[%2d]/F",ndim*2);
493  if(ndim<5) sprintf(cbx, "bx[%1d]/F",ndim*2);
494  else sprintf(cbx, "bx[%2d]/F",ndim*2);
495 
496  if(ndim<5) sprintf(ctime, "time[%1d]/D",ndim*2);
497  else sprintf(ctime, "time[%2d]/D",ndim*2);
498 
499  sprintf(cright, "right[%1d]/F",ndim);
500  sprintf(cleft, "left[%1d]/F",ndim);
501  sprintf(cduration, "duration[%1d]/F",ndim);
502  sprintf(cstart, "start[%1d]/D",ndim);
503  sprintf(cstop, "stop[%1d]/D",ndim);
504 
505  sprintf(cfrequency, "frequency[%1d]/F",ndim);
506  sprintf(clow, "low[%1d]/F",ndim);
507  sprintf(chigh, "high[%1d]/F",ndim);
508  sprintf(cbandwidth, "bandwidth[%1d]/F",ndim);
509 
510  if(ndim<5) sprintf(chrss, "hrss[%1d]/D",ndim*2);
511  else sprintf(chrss, "hrss[%2d]/D",ndim*2);
512 
513  sprintf(cnoise, "noise[%1d]/D",ndim);
514  sprintf(cnull, "null[%1d]/F",ndim);
515  sprintf(cnill, "nill[%1d]/F",ndim);
516  sprintf(cDeff, "Deff[%1d]/F",ndim);
517  sprintf(crange, "range[2]/F");
518  sprintf(ceBBH, "eBBH[4]/F");
519  sprintf(cchirp, "chirp[6]/F");
520  sprintf(cmass, "mass[2]/F");
521  sprintf(cspin, "spin[6]/F");
522  sprintf(csnr, "snr[%1d]/F",ndim);
523  sprintf(csSNR, "sSNR[%1d]/F",ndim);
524  sprintf(cxSNR, "xSNR[%1d]/F",ndim);
525  sprintf(ciSNR, "iSNR[%1d]/F",ndim);
526  sprintf(coSNR, "oSNR[%1d]/F",ndim);
527  sprintf(cioSNR, "ioSNR[%1d]/F",ndim);
528 
529  //==================================
530  // Define trigger tree
531  //==================================
532 
533  waveTree->Branch("ndim", &ndim, "ndim/I");
534  waveTree->Branch("run", &run, "run/I");
535  waveTree->Branch("nevent", &nevent, "nevent/I");
536  waveTree->Branch("eventID", eventID, "eventID[2]/I");
537  waveTree->Branch("type", type, "type[2]/I");
538  waveTree->Branch("name", name);
539  waveTree->Branch("log", log);
540  waveTree->Branch("rate", rate, crate);
541 
542  waveTree->Branch("usize", &usize, "usize/I");
543  waveTree->Branch("volume", volume, cvolume);
544  waveTree->Branch("size", size, csize);
545 
546  waveTree->Branch("gap", gap, cgap);
547  waveTree->Branch("lag", lag, clag);
548  waveTree->Branch("slag", slag, cslag);
549  waveTree->Branch("strain", strain, "strain[2]/D");
550  waveTree->Branch("phi", phi, "phi[4]/F");
551  waveTree->Branch("theta", theta, "theta[4]/F");
552  waveTree->Branch("psi", psi, "psi[2]/F");
553  waveTree->Branch("iota", iota, "iota[2]/F");
554  waveTree->Branch("bp", bp, cbp);
555  waveTree->Branch("bx", bx, cbx);
556 
557  waveTree->Branch("time", time, ctime);
558  waveTree->Branch("gps", gps, cgps);
559  waveTree->Branch("left", left, cleft);
560  waveTree->Branch("right", right, cright);
561  waveTree->Branch("start", start, cstart);
562  waveTree->Branch("stop", stop, cstop);
563  waveTree->Branch("duration", duration, cduration);
564 
565  waveTree->Branch("frequency", frequency, cfrequency);
566  waveTree->Branch("low", low, clow);
567  waveTree->Branch("high", high, chigh);
568  waveTree->Branch("bandwidth", bandwidth, cbandwidth);
569 
570  waveTree->Branch("hrss", hrss, chrss);
571  waveTree->Branch("noise", noise, cnoise);
572  //waveTree->Branch("ndm", ndm, cndm);
573  waveTree->Branch("erA", erA, "erA[11]/F");
574 
575  waveTree->Branch("Psm","skymap",&Psm, 32000,0);
576 
577  waveTree->Branch("null", null, cnull);
578  waveTree->Branch("nill", nill, cnill);
579  waveTree->Branch("netcc", netcc, "netcc[4]/F");
580  waveTree->Branch("neted", neted, "neted[5]/F");
581  waveTree->Branch("rho", rho, "rho[2]/F");
582 
583  waveTree->Branch("gnet", &gnet, "gnet/F");
584  waveTree->Branch("anet", &anet, "anet/F");
585  waveTree->Branch("inet", &inet, "inet/F");
586  waveTree->Branch("ecor", &ecor, "ecor/F");
587  waveTree->Branch("norm", &norm, "norm/F");
588  waveTree->Branch("ECOR", &ECOR, "ECOR/F");
589  waveTree->Branch("penalty", &penalty, "penalty/F");
590  waveTree->Branch("likelihood", &likelihood, "likelihood/F");
591 
592  waveTree->Branch("factor", &factor, "factor/F");
593  waveTree->Branch("chirp", chirp, cchirp);
594  waveTree->Branch("range", range, crange);
595  waveTree->Branch("eBBH", eBBH, ceBBH);
596  waveTree->Branch("Deff", Deff, cDeff);
597  waveTree->Branch("mass", mass, cmass);
598  waveTree->Branch("spin", spin, cspin);
599 
600  waveTree->Branch("snr", snr, csnr);
601  waveTree->Branch("xSNR", xSNR, cxSNR);
602  waveTree->Branch("sSNR", sSNR, csSNR);
603  waveTree->Branch("iSNR", iSNR, ciSNR);
604  waveTree->Branch("oSNR", oSNR, coSNR);
605  waveTree->Branch("ioSNR", ioSNR, cioSNR);
606 
607  return waveTree;
608 }
609 
611 {
612  int i,j;
613 
614  if(ndim != a.ndim) { ndim=a.ndim; this->allocate(); }
615 
616  ifoList= a.ifoList;
617 
618  run= a.run;
619  nevent= a.nevent;
620  eventID[0]= a.eventID[0];
621  eventID[1]= a.eventID[1];
622  type[0]= a.type[0];
623  type[1]= a.type[1];
624  *name= *a.name;
625  *log= *a.log;
626  strain[0]= a.strain[0];
627  strain[1]= a.strain[1];
628  usize= a.usize;
629  psi[0]= a.psi[0];
630  psi[1]= a.psi[1];
631  iota[0]= a.iota[0];
632  iota[1]= a.iota[1];
633  mass[0]= a.mass[0];
634  mass[1]= a.mass[1];
635  gnet= a.gnet;
636  anet= a.anet;
637  inet= a.inet;
638  ecor= a.ecor;
639  norm= a.norm;
640  ECOR= a.ECOR;
641  penalty= a.penalty;
642  likelihood= a.likelihood;
643  factor= a.factor;
644  Psave= a.Psave;
645  *Psm = *(a.Psm);
646 
647  for(i=0; i<6; i++) spin[i] = a.spin[i];
648  for(i=0; i<4; i++) netcc[i] = a.netcc[i];
649  for(i=0; i<5; i++) neted[i] = a.neted[i];
650  for(i=0; i<2; i++) rho[i] = a.rho[i];
651  for(i=0; i<4; i++) theta[i] = a.theta[i];
652  for(i=0; i<4; i++) phi[i] = a.phi[i];
653  for(i=0; i<11; i++) erA[i] = a.erA[i];
654  for(i=0; i<2; i++) range[i] = a.range[i];
655  for(i=0; i<4; i++) eBBH[i] = a.eBBH[i];
656  for(i=0; i<6; i++) chirp[i] = a.chirp[i];
657 
658  lag[ndim]= a.lag[ndim];
659  slag[ndim]= a.slag[ndim];
660  for(i=0; i<ndim; i++){
661 
662  gps[i]= a.gps[i];
663 
664  rate[i]= a.rate[i];
665 
666  gap[i]= a.gap[i];
667  lag[i]= a.lag[i];
668  slag[i]= a.slag[i];
669  volume[i]= a.volume[i];
670  size[i]= a.size[i];
671  bp[i]= a.bp[i];
672  bp[i+ndim]= a.bp[i+ndim];
673  bx[i]= a.bx[i];
674  bx[i+ndim]= a.bx[i+ndim];
675 
676  time[i]= a.time[i];
677  time[i+ndim]= a.time[i+ndim];
678  right[i]= a.right[i];
679  left[i]= a.left[i];
680  duration[i]= a.duration[i];
681  start[i]= a.start[i];
682  stop[i]= a.stop[i];
683 
684  frequency[i]= a.frequency[i];
685  low[i]= a.low[i];
686  high[i]= a.high[i];
687  bandwidth[i]= a.bandwidth[i];
688 
689  hrss[i]= a.hrss[i];
690  hrss[i+ndim]= a.hrss[i+ndim];
691  noise[i]= a.noise[i];
692  null[i]= a.null[i];
693  nill[i]= a.nill[i];
694  Deff[i]= a.Deff[i];
695 
696  snr[i]= a.snr[i];
697  sSNR[i]= a.sSNR[i];
698  xSNR[i]= a.xSNR[i];
699  iSNR[i]= a.iSNR[i];
700  oSNR[i]= a.oSNR[i];
701  ioSNR[i]= a.ioSNR[i];
702 
703  }
704  return *this;
705 }
706 
707 
708 //++++++++++++++++++++++++++++++++++++++++++++++++++++
709 // output network events for 2G analysis
710 //++++++++++++++++++++++++++++++++++++++++++++++++++++
711 void netevent::output2G(TTree* waveTree, network* net, size_t ID, int LAG, double factor)
712 {
713 // output event ID in lag LAG to root file
714  if(!net || !net->nLag) exit(1);
715 
716  int i,j,k,ind,kid;
717  int n,m,K,M,injID;
718  detector* pd;
719  netcluster* p;
720  int nIFO = (int)net->ifoListSize(); // number of detectors
721  wavecomplex Aa, gC;
722  double Pi = 3.14159265358979312;
723  this->factor = fabs(factor); // used to tag events in root file
724  factor = factor<=0 ? 1 : fabs(factor); // used to rescale injected events
725  double inRate = net->getifo(0)->rate; // get original rate
726  bool pat0 = net->pattern==0 ? true : false; // pattern flag
727 
728  if(!nIFO) return;
729 
730 // injections
731 
732  injection INJ(nIFO);
733  size_t mdcID, id;
734  double T, injTime, TAU, pcc;
735  double x, ndmLike, Etot;
736  skymap* psm;
737  netcluster* pwc;
738  std::vector<float>* vP;
739  std::vector<int>* vI;
740 
741  bool ellips = net->tYPe=='i' || net->tYPe=='I' ||
742  net->tYPe=='g' || net->tYPe=='G' ||
743  net->tYPe=='s' || net->tYPe=='S' ||
744  net->tYPe=='r' || net->tYPe=='R';
745 
746  bool burst = net->tYPe=='b' || net->tYPe=='B';
747 
748  // arrays for cluster parameters
749 
750  wavearray<double> clusterID_net;
751  wavearray<double> vol0_net;
752  wavearray<double> vol1_net;
753  wavearray<double> size_net;
754  wavearray<double> start_net;
755  wavearray<double> stop_net;
756  wavearray<double> low_net;
757  wavearray<double> high_net;
758  wavearray<double> noise_net;
759  wavearray<double> NOISE_net;
760  wavearray<double> cFreq_net;
761  wavearray<double> rate_net;
762  wavearray<double> duration_net;
763  wavearray<double> bandwidth_net;
764 
765  // why do we need all this?
766 
767  this->ifoList = net->ifoList; // add detectors to tree user info
768  if(waveTree!=NULL) {
769  int dsize = waveTree->GetUserInfo()->GetSize();
770  if(dsize!=0 && dsize!=nIFO) {
771  cout<<"netevent::output2G(): wrong user detector list in header tree"<<endl; exit(1);
772  }
773  if(dsize==0) {
774  for(int n=0;n<nIFO;n++) {
775  // must be done a copy because detector object
776  // is destroyed when waveTree is destroyed
777  detector* pD = new detector(*net->getifo(n));
778  waveTree->GetUserInfo()->Add(pD);
779  }
780  }
781  }
782 
783  this->run = net->nRun;
784  this->nevent = 0;
785 
786  pwc = net->getwc(LAG); // pointer to netcluster
787  if(!pwc->size()) return;
788 
789  clusterID_net = pwc->get((char*)"ID",0,'S',0);
790  K = clusterID_net.size();
791  kid = -1;
792  for(k=0; k<K; k++) { // find cluster
793  id = size_t(clusterID_net.data[k]+0.1);
794  if(ID&&(ID==id)) {kid=k; break;}
795  }
796  if(kid<0) return;
797 
798 // read cluster parameters
799 
800  start_net.resize(0);
801  stop_net.resize(0);
802  noise_net.resize(0);
803  NOISE_net.resize(0);
804 
805  rate_net = pwc->get((char*)"rate",0,'R',0);
806  vol0_net = pwc->get((char*)"volume",0,'R',0); // stored in volume[0]
807  vol1_net = pwc->get((char*)"VOLUME",0,'R',0); // stored in volume[1]
808  size_net = pwc->get((char*)"size",0,'R',0); // stored in size[0]
809  low_net = pwc->get((char*)"low",0,'R',0);
810  high_net = pwc->get((char*)"high",0,'R',0);
811  cFreq_net = pwc->get((char*)"freq",0,'L',0,false);
812  duration_net = pwc->get((char*)"duration",0,'L',0,false);
813  bandwidth_net = pwc->get((char*)"bandwidth",0,'L',0,false);
814 
815  for(i=1; i<=int(nIFO); i++) { // loop on detectors
816  start_net.append(pwc->get((char*)"start",i,'L',0));
817  stop_net.append(pwc->get((char*)"stop",i,'L',0));
818  noise_net.append(pwc->get((char*)"noise",i,'S',0));
819  NOISE_net.append(pwc->get((char*)"NOISE",i,'S',0));
820  }
821 
822  psm = &(net->getifo(0)->tau);
823  vI = &(net->wc_List[LAG].p_Ind[ID-1]);
824  ind = (*vI)[0]; // reconstructed sky index
825 
826  for(i=0; i<int(nIFO); i++) // store gps time of data intervals
827  this->gps[i] = pwc->start+(slag[i]-slag[0]);
828 
829  clusterdata* pcd = &(pwc->cData[ID-1]);
830 
831  this->ndim = nIFO;
832  psm->gps = pcd->cTime + this->gps[0];
833  this->ecor = pcd->netecor;
834  this->nevent += 1;
835  this->eventID[0] = ID;
836  this->eventID[1] = 0;
837  this->iota[0] = pcd->iota; // ellipticity : 2*cos(iota)/(1+cos(iota)^2)
838  //this->iota[0] = pcd->skyChi2; // sky chi2 for all resolutions (temporary)
839  this->psi[0] = pcd->psi;
840  //this->psi[0] = pcd->Gnoise; // estimated Gaussian noise (temporary)
841  this->phi[0] = psm->getPhi(ind); // reconstructed phi
842  this->phi[2] = psm->getRA(ind);
843  this->phi[3] = pcd->phi; // detection phi
844  this->theta[0] = psm->getTheta(ind); // reconstructed theta
845  this->theta[2] = psm->getDEC(ind);
846  this->theta[3] = pcd->theta; // detection theta
847  this->gnet = pcd->gNET;
848  this->anet = pcd->aNET;
849  this->inet = pcd->iNET; // network index
850  this->norm = pcd->norm; // packet norm
851  this->likelihood = pcd->likenet; // total likelihood in sky loop
852  this->volume[0] = int(vol0_net.data[kid]+0.5); // event volume
853  this->volume[1] = int(vol1_net.data[kid]+0.5); // selected pixels volume
854  this->size[0] = int(size_net.data[kid]+0.5); // event size
855  this->size[1] = pcd->skySize; // signal size
856  this->chirp[1] = pcd->mchirp; // reconstructed chirp mass
857  this->chirp[2] = pcd->mchirperr; // reconstructed chirp mass error
858  this->chirp[3] = pcd->chirpEllip; // ellipticity parameter
859  this->chirp[4] = pcd->chirpPfrac; // pixel fraction
860  this->chirp[5] = pcd->chirpEfrac; // energy fraction
861  //this->chirp[6] = pcd->chi2chirp; // chirp chi2/DoF
862  //this->chirp[7] = pcd->tmrgr; // t merger parameter
863  //this->chirp[8] = pcd->tmrgrerr; // t merger error
864 
865  this->range[0] = 0.; // reconstructed distance to source
866 
867  TAU = psm->get(this->theta[0],this->phi[0]);
868 
869  M = 0; gC = 0.;
870  this->strain[0] = 0.;
871  this->penalty = 0.;
872 
873  for(i=0; i<5; i++) { this->neted[i] = 0.; } // zero neted array
874 
875  this->lag[nIFO] = pwc->shift;
876 
877  net->getMRAwave(ID,LAG,'s',net->optim?1:0); // get signal strain
878 
879  for(i=0; i<int(nIFO); i++) { // loop on detectors
880  pd = net->getifo(i);
881  Aa = pd->antenna(this->theta[0],this->phi[0],this->psi[0]); // antenna pattern
882  m = i*K+kid;
883 
884  this->type[0] = 1;
885  this->rate[i] = net->optim ? int(rate_net.data[kid]+0.1) : 0;
886 
887  this->gap[i] = 0.;
888  this->lag[i] = pd->lagShift.data[LAG];
889 
890  this->snr[i] = pd->enrg;
891  this->nill[i] = pd->xSNR-pd->sSNR;
892  this->null[i] = pd->null; // null per detector
893  this->xSNR[i] = pd->xSNR;
894  this->sSNR[i] = pd->sSNR;
895  this->time[i] = pcd->cTime + this->gps[i];
896 
897  if(i) { // take delays into account
898  psm = &(net->getifo(i)->tau);
899  this->time[i] += psm->get(this->theta[0],this->phi[0])-TAU;
900  }
901 
902  this->left[i] = start_net.data[m];
903  this->right[i] = pwc->stop-pwc->start-stop_net.data[m];
904  this->duration[i] = stop_net.data[m] - start_net.data[m];
905  this->start[i] = start_net.data[m] + this->gps[i];
906  this->stop[i] = stop_net.data[m] + this->gps[i];
907 
908  // take lag shift into account
909  double xstart = this->gps[i]+net->Edge; // start data
910  double xstop = this->gps[i]+pwc->stop-pwc->start-net->Edge; // end data
911  this->time[i] += lag[i];
912  if(this->time[i]>xstop) this->time[i] = xstart+(this->time[i]-xstop); // circular buffer
913 
914  this->frequency[i] = cFreq_net.data[k];
915  this->low[i] = low_net.data[k];
916  this->high[i] = high_net.data[k];
917  this->bandwidth[i] = high_net.data[k] - low_net.data[k];
918 
919  this->hrss[i] = sqrt(pd->get_SS()/inRate);
920  this->noise[i] = pow(10.,noise_net.data[m])/sqrt(inRate);
921  this->bp[i]= Aa.real();
922  this->bx[i]= Aa.imag();
923  this->strain[0] += this->hrss[i]*this->hrss[i];
924 
925  Aa /= pow(10.,NOISE_net.data[m]);
926  gC += Aa*Aa;
927 
928  psm->gps = pcd->cTime+this->gps[0];
929  }
930 
931  // Fix start,stop,duration,left,right when simulation!=0.
932  // Due to circular buffer the events detected on the edges
933  // of the segment could have a wrong start,stop values
934  // The start,stop of det with lag=0 are used to fix the wrong values of the other detectors
935  int mdet=0; for(i=0; i<int(nIFO); i++) if(this->lag[i]==0) mdet=i;
936  for(i=0; i<int(nIFO); i++) {
937  if(this->duration[i]!=this->duration[mdet]) {
938  double xstart = this->gps[i]+net->Edge; // start data
939  double xstop = this->gps[i]+pwc->stop-pwc->start-net->Edge; // end data
940 
941  this->start[i] = this->start[mdet]+this->lag[i]+(slag[i]-slag[mdet]);
942  if(this->start[i]>xstop) this->start[i] = xstart+(this->start[i]-xstop); // circular buffer
943 
944  this->stop[i] = this->stop[mdet]+this->lag[i]+(slag[i]-slag[mdet]);
945  if(this->stop[i]>xstop) this->stop[i] = xstart+(this->stop[i]-xstop); // circular buffer
946 
947  this->duration[i] = this->duration[mdet];
948  this->left[i] = this->start[i]-pwc->start;
949  this->right[i] = pwc->stop-this->stop[i];
950  }
951  }
952  this->duration[0] = duration_net.data[kid];
953  this->bandwidth[0] = bandwidth_net.data[kid];
954  this->frequency[0] = pcd->cFreq;
955 
956  ind = pwc->sArea[ID-1].size();
957  for(i=0; i<11; i++)
958  this->erA[i] = i<ind ? pwc->sArea[ID-1][i] : 0.;
959 
960  this->ECOR = pcd->normcor; // normalized coherent energy
961  this->netcc[0] = pcd->netcc; // MRA or SRA cc statistic
962  this->netcc[1] = pcd->skycc; // all-resolution cc statistic
963  this->netcc[2] = pcd->subnet; // MRA or SRA sub-network statistic
964  this->netcc[3] = pcd->SUBNET; // all-resolution sub-network statistic
965 
966  this->neted[0] = pcd->netED; // network ED
967  this->neted[1] = pcd->netnull; // total null energy with Gaussian bias correction
968  this->neted[2] = pcd->energy; // total event energy
969  this->neted[3] = pcd->likesky; // total likelihood at all resolutions
970  this->neted[4] = pcd->enrgsky; // total energy at all resolutions
971  //this->neted[5] = pcd->Gnoise; // estimated contribution of Gaussian noise
972  //this->neted[6] = pcd->skyChi2; // sky chi2 for all resolutions
973 
974  double chrho = this->chirp[3]*sqrt(this->chirp[5]); // reduction factor for chirp events
975  this->rho[0] = pcd->netRHO; // reduced coherent SNR per detector
976  this->rho[1] = pat0 ?pcd->netrho:pcd->netRHO*chrho; // reduced coherent SNR per detector for chirp events
977 
978  this->strain[0] = sqrt(this->strain[0]);
979  if(!ellips) this->psi[0] = gC.arg()*180/Pi;
980  this->penalty = pcd->netnull/nIFO;
981  this->penalty /= pat0 ? this->size[0] : pcd->nDoF; // cluster chi2/nDoF
982 
983 
984 // set injections if there are any
985 
986  M = net->mdc__IDSize();
987  if(!LAG) { // only for zero lag
988  injTime = 1.e12;
989  injID = -1;
990  for(m=0; m<M; m++) {
991  mdcID = net->getmdc__ID(m);
992  T = fabs(this->time[0] - net->getmdcTime(mdcID));
993  if(T<injTime && INJ.fill_in(net,mdcID)) {
994  injTime = T;
995  injID = mdcID;
996  }
997  // printf("%d %12.3f %12.3f\n",mdcID,net->getmdcTime(mdcID),T);
998  }
999 
1000  if(INJ.fill_in(net,injID)) { // set injections
1001  this->range[1] = INJ.distance/factor;
1002  this->chirp[0] = INJ.mchirp;
1003  this->eBBH[1] = INJ.e0;
1004  this->eBBH[2] = INJ.rp0;
1005  this->eBBH[3] = INJ.redshift;
1006  this->strain[1] = INJ.strain*factor;
1007  this->type[1] = INJ.type;
1008  *this->name = net->getmdcType(this->type[1]-1);
1009  *this->log = net->getmdcList(injID);
1010  this->log->erase(std::remove(this->log->begin(), this->log->end(), '\n'), this->log->end()); // remove new line
1011  this->theta[1] = INJ.theta[0];
1012  this->phi[1] = INJ.phi[0];
1013  this->psi[1] = INJ.psi[0];
1014  // this->iota[1] = 2.*INJ.iota[0]/(1+INJ.iota[0]*INJ.iota[0]); // ellipticity
1015  this->iota[1] = INJ.iota[1]; // cos(iota)
1016  this->mass[0] = INJ.mass[0];
1017  this->mass[1] = INJ.mass[1];
1018 
1019  for(i=0; i<6; i++) this->spin[i] = INJ.spin[i];
1020 
1021 // printf("injection type: %d %12.3f\n",INJ.type,INJ.time[0]);
1022 
1023  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
1024  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
1025  pd = net->getifo(0);
1026  int idSize = pd->RWFID.size();
1027  int wfIndex=-1;
1028  for (int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
1029 
1030  for(j=0; j<nIFO; j++) {
1031  pd = net->getifo(j);
1032  Aa = pd->antenna(this->theta[1],this->phi[1],this->psi[1]); // inj antenna pattern
1033  this->hrss[j+nIFO] = INJ.hrss[j]*factor;
1034  this->bp[j+nIFO] = Aa.real();
1035  this->bx[j+nIFO] = Aa.imag();
1036  this->time[j+nIFO] = INJ.time[j];
1037  this->Deff[j] = INJ.Deff[j]/factor;
1038 
1039  pwfINJ[j] = INJ.pwf[j];
1040  if (pwfINJ[j]==NULL) {
1041  cout << "Error : Injected waveform not saved !!! : detector "
1042  << net->ifoName[j] << endl;
1043  continue;
1044  }
1045  if (wfIndex<0) {
1046  cout << "Error : Reconstructed waveform not saved !!! : ID -> "
1047  << ID << " : detector " << net->ifoName[j] << endl;
1048  continue;
1049  }
1050 
1051  if (wfIndex>=0) pwfREC[j] = pd->RWFP[wfIndex];
1052  double R = pd->TFmap.rate();
1053  //double rFactor = log(2.)/log(pd->rate/R); // rescale waveform
1054  double rFactor = 1.;
1055  rFactor *= factor;
1056  wavearray<double>* wfINJ = pwfINJ[j];
1057  *wfINJ*=rFactor;
1058  wavearray<double>* wfREC = pwfREC[j];
1059 
1060  double bINJ = wfINJ->start();
1061  double eINJ = wfINJ->start()+wfINJ->size()/R;
1062  double bREC = wfREC->start();
1063  double eREC = wfREC->start()+wfREC->size()/R;
1064  //cout.precision(14);
1065  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
1066  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
1067 
1068  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
1069  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
1070  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
1071 
1072  double startXCOR = bINJ>bREC ? bINJ : bREC;
1073  double endXCOR = eINJ<eREC ? eINJ : eREC;
1074  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1075  //cout << "startXCOR : " << startXCOR << " endXCOR : " << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
1076 
1077  if (sizeXCOR<=0) {*wfINJ*=1./rFactor; continue;}
1078 
1079  // the enINJ, enREC, xcorINJ_REC are computed in the INJ range
1080 
1081  double enINJ=0;
1082  for (int i=0;i<wfINJ->size();i++) enINJ+=wfINJ->data[i]*wfINJ->data[i];
1083  //for (int i=0;i<sizeXCOR;i++) enINJ+=wfINJ->data[i+oINJ]*wfINJ->data[i+oINJ];
1084 
1085  double enREC=0;
1086  for (int i=0;i<wfREC->size();i++) enREC+=wfREC->data[i]*wfREC->data[i];
1087  //for (int i=0;i<sizeXCOR;i++) enREC+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
1088 
1089  double xcorINJ_REC=0;
1090  for (int i=0;i<sizeXCOR;i++) xcorINJ_REC+=wfINJ->data[i+oINJ]*wfREC->data[i+oREC];
1091 
1092  WSeries<double> wfREC_SUB_INJ;
1093  wfREC_SUB_INJ.resize(sizeXCOR);
1094  for (int i=0;i<sizeXCOR;i++) wfREC_SUB_INJ.data[i]=wfREC->data[i+oREC]-wfINJ->data[i+oINJ];
1095  wfREC_SUB_INJ.start(startXCOR);
1096  wfREC_SUB_INJ.rate(wfREC->rate());
1097 
1098  this->iSNR[j] = enINJ;
1099  this->oSNR[j] = enREC;
1100  this->ioSNR[j] = xcorINJ_REC;
1101 
1102  //double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
1103  //cout << "enINJ : " << enINJ << " enREC : " << enREC << " xcorINJ_REC : " << xcorINJ_REC << endl;
1104  //cout << "erINJ_REC/enINJ : " << erINJ_REC/enINJ << endl;
1105 
1106  *wfINJ*=1./rFactor;
1107  }
1108  delete [] pwfINJ;
1109  delete [] pwfREC;
1110  }
1111  else { // no injections
1112  this->range[1] = 0.;
1113  this->chirp[0] = 0.;
1114  this->eBBH[1] = 0.;
1115  this->eBBH[2] = 0.;
1116  this->eBBH[3] = 0.;
1117  this->strain[1] = 0.;
1118  this->type[1] = 0;
1119  this->theta[1] = 0.;
1120  this->phi[1] = 0.;
1121  this->psi[1] = 0.;
1122  this->iota[1] = 0.;
1123  this->mass[0] = 0.;
1124  this->mass[1] = 0.;
1125 
1126  for(i=0; i<6; i++) this->spin[i] = 0.;
1127 
1128  for(j=0; j<nIFO; j++) {
1129  this->hrss[j+nIFO] = 0.;
1130  this->bp[j+nIFO] = 0.;
1131  this->bx[j+nIFO] = 0.;
1132  this->time[j+nIFO] = 0.;
1133  this->Deff[j] = 0.;
1134  this->iSNR[j] = 0.;
1135  this->oSNR[j] = 0.;
1136  this->ioSNR[j] = 0.;
1137  }
1138  }
1139  }
1140 
1141  if(this->fP!=NULL) {
1142  fprintf(fP,"\n# trigger %d in lag %d for \n",int(ID),int(LAG));
1143  this->Dump("2G");
1144  vP = &(net->wc_List[LAG].p_Map[ID-1]);
1145  vI = &(net->wc_List[LAG].p_Ind[ID-1]);
1146  x = cos(psm->theta_1*PI/180.)-cos(psm->theta_2*PI/180.);
1147  x*= (psm->phi_2-psm->phi_1)*180/PI/psm->size();
1148  fprintf(fP,"sky_res: %f\n",sqrt(fabs(x)));
1149  fprintf(fP,"map_lenght: %d\n",int(vP->size()));
1150  fprintf(fP,"#skyID theta DEC step phi R.A step probability cumulative\n");
1151  x = 0;
1152  for(j=0; j<int(vP->size()); j++) {
1153  i = (*vI)[j];
1154  if(net->mdc__IDSize()) { // simulation mode
1155  x = (j==int(vP->size())-1) ? 0 : x+(*vP)[j]; // last value is the inj sky index (x=0)
1156  } else x+=(*vP)[j];
1157  fprintf(fP,"%6d %5.1f %5.1f %6.2f %5.1f %5.1f %6.2f %e %e\n",
1158  int(i),psm->getTheta(i),psm->getDEC(i),psm->getThetaStep(i),
1159  psm->getPhi(i),psm->getRA(i),psm->getPhiStep(i),(*vP)[j],x);
1160  }
1161  }
1162 
1163  // save the probability skymap to tree
1164  if(waveTree!=NULL && net->wc_List[LAG].p_Map.size()) {
1165 
1166  vP = &(net->wc_List[LAG].p_Map[ID-1]);
1167  vI = &(net->wc_List[LAG].p_Ind[ID-1]);
1168  if(this->Psave) {
1169 
1170  *Psm = *psm;
1171  *Psm = 0.;
1172 
1173  int k;
1174  double th,ph;
1175  for(j=0; j<int(vP->size()); j++) {
1176  i = (*vI)[j];
1177  th = Psm->getTheta(i);
1178  ph = Psm->getPhi(i);
1179  k=Psm->getSkyIndex(th, ph);
1180  Psm->set(k,(*vP)[j]);
1181  }
1182  }
1183  }
1184 
1185  if(waveTree!=NULL) waveTree->Fill();
1186 
1187 }
1188 
1189 
1190 //++++++++++++++++++++++++++++++++++++++++++++++++++++
1191 // output network events
1192 //++++++++++++++++++++++++++++++++++++++++++++++++++++
1193 void netevent::output(TTree* waveTree, network* net, double factor, size_t iID, int LAG)
1194 {
1195  if(!net || !net->nLag) exit(1);
1196 
1197  int i,j,k,ind;
1198  int n,m,K,M,injID;
1199  detector* pd;
1200  netcluster* p;
1201  int N = (int)net->ifoListSize(); // number of detectors
1202  int bLag = LAG<0 ? 0 : LAG;
1203  int eLag = LAG<0 ? net->nLag : LAG+1;
1204  wavecomplex Aa, gC;
1205  double Pi = 3.14159265358979312;
1206  this->factor = fabs(factor); // used to tag events in root file
1207  factor = factor<=0 ? 1 : fabs(factor); // used to rescale injected events
1208  double inRate = net->getifo(0)->rate; // get original rate
1209 
1210  if(!N) return;
1211  double* ndm = (double*)malloc(N*N*sizeof(double));
1212  int iTYPE = net->MRA ? 0 : 1;
1213 
1214 // injections
1215 
1216  injection INJ(N);
1217  size_t mdcID, ID;
1218  double T, injTime, TAU, pcc;
1219  double x, tot_null, ndmLike, Etot;
1220  skymap* psm;
1221  std::vector<float>* vP;
1222  std::vector<int>* vI;
1223 
1224  bool ellips = net->tYPe=='i' || net->tYPe=='I' ||
1225  net->tYPe=='g' || net->tYPe=='G' ||
1226  net->tYPe=='s' || net->tYPe=='S' ||
1227  net->tYPe=='r' || net->tYPe=='R';
1228 
1229  bool burst = net->tYPe=='b' || net->tYPe=='B';
1230 
1231  wavearray<double> skSNR(N);
1232  wavearray<double> xkSNR(N);
1233 
1234 // arrays for cluster parameters
1235 
1236  wavearray<double> clusterID_net;
1237  wavearray<double> volume_net;
1238  wavearray<double> size_net;
1239  wavearray<double> start_net;
1240  wavearray<double> stop_net;
1241  wavearray<double> time_net;
1242  wavearray<double> frequency_net;
1243  wavearray<double> low_net;
1244  wavearray<double> high_net;
1245  wavearray<double> LH_net;
1246  wavearray<double> null_net;
1247  wavearray<double> rSNR_net;
1248  wavearray<double> gSNR_net;
1249  wavearray<double> gSNR_NET;
1250  wavearray<double> rate_net;
1251  wavearray<double> hrss_net;
1252  wavearray<double> hrss_NET;
1253  wavearray<double> noise_net;
1254  wavearray<double> NOISE_net;
1255  wavearray<double> psi_net;
1256  wavearray<double> ell_net;
1257  wavearray<double> phi_net;
1258  wavearray<double> theta_net;
1259  wavearray<double> energy_net;
1260  wavearray<double> energy_NET;
1261 
1262  this->ifoList = net->ifoList;
1263  // add detectors to tree user info
1264  if(waveTree!=NULL) {
1265  for(int n=0;n<N;n++) {
1266  // must be done a copy because detector object
1267  // is destroyed when waveTree is destroyed
1268  detector* pD = new detector(*net->getifo(n));
1269  waveTree->GetUserInfo()->Add(pD);
1270  }
1271  }
1272 
1273  this->run = net->nRun;
1274  this->nevent = 0;
1275 
1276  for(n=bLag; n<eLag; n++){ // loop on time lags
1277 
1278  p = net->getwc(n); // pointer to netcluster
1279  if(!p->size()) continue;
1280 
1281  clusterID_net = p->get((char*)"ID",0,'S',iTYPE);
1282  K = clusterID_net.size();
1283 
1284  if(!K) continue;
1285 
1286 // read cluster parameters
1287 
1288  time_net.resize(0);
1289  start_net.resize(0);
1290  stop_net.resize(0);
1291  frequency_net.resize(0);
1292  energy_net.resize(0);
1293  energy_NET.resize(0);
1294  rSNR_net.resize(0);
1295  gSNR_net.resize(0);
1296  gSNR_NET.resize(0);
1297  hrss_net.resize(0);
1298  hrss_NET.resize(0);
1299  noise_net.resize(0);
1300  NOISE_net.resize(0);
1301  null_net.resize(0);
1302 
1303  LH_net = p->get((char*)"likelihood",0,'R',iTYPE);
1304  rate_net = p->get((char*)"rate",0,'R',iTYPE);
1305  ell_net = p->get((char*)"ellipticity",0,'R',iTYPE);
1306  psi_net = p->get((char*)"psi",0,'R',iTYPE);
1307  phi_net = p->get((char*)"phi",0,'R',iTYPE);
1308  theta_net = p->get((char*)"theta",0,'R',iTYPE);
1309  size_net = p->get((char*)"size",0,'R',iTYPE);
1310  volume_net = p->get((char*)"volume",0,'R',iTYPE);
1311  low_net = p->get((char*)"low",0,'R',iTYPE);
1312  high_net = p->get((char*)"high",0,'R',iTYPE);
1313 
1314  for(i=1; i<int(N+1); i++) // loop on detectors
1315  { // read cluster parameters
1316  time_net.append(p->get((char*)"time",i,'L',0));
1317  start_net.append(p->get((char*)"start",i,'L',iTYPE));
1318  stop_net.append(p->get((char*)"stop",i,'L',iTYPE));
1319  frequency_net.append(p->get((char*)"FREQUENCY",i,'L',0));
1320  rSNR_net.append(p->get((char*)"SNR",i,'R',iTYPE));
1321  gSNR_net.append(p->get((char*)"SNR",i,'S',iTYPE));
1322  gSNR_NET.append(p->get((char*)"SNR",i,'P',iTYPE));
1323  hrss_net.append(p->get((char*)"hrss",i,'W',iTYPE));
1324  hrss_NET.append(p->get((char*)"hrss",i,'U',iTYPE));
1325  noise_net.append(p->get((char*)"noise",i,'S',iTYPE));
1326  NOISE_net.append(p->get((char*)"NOISE",i,'S',iTYPE));
1327  energy_net.append(p->get((char*)"energy",i,'S',iTYPE));
1328  energy_NET.append(p->get((char*)"energy",i,'P',iTYPE));
1329  null_net.append(p->get((char*)"null",i,'W',iTYPE));
1330  }
1331 
1332  if(ellips) {
1333  energy_net += energy_NET;
1334  energy_net *= net->MRA ? 1.0 : 0.5;
1335  }
1336 
1337  for(k=0; k<K; k++) // loop on events
1338  {
1339  ID = size_t(clusterID_net.data[k]+0.1);
1340  if((iID)&&(ID!=iID)) continue;
1341  if(ellips) {
1342  net->SETNDM(ID,n,true,net->MRA?0:1);
1343  }
1344  else if(burst)
1345  net->setndm(ID,n,true,1);
1346 // else
1347 // cout<<"netevent::output(): incorrect search option"<<endl; //exit(1);}
1348 
1349  psm = &(net->getifo(0)->tau);
1350  vI = &(net->wc_List[n].p_Ind[ID-1]);
1351  ind = (*vI)[0]; // reconstructed sky index
1352 
1353  this->ndim = N;
1354  for(i=0; i<N; i++) // loop on detectors
1355  this->gps[i] = net->getifo(i)->getTFmap()->start()+(slag[i]-slag[0]);
1356  psm->gps = time_net.data[k]+this->gps[0];
1357  this->ecor = net->eCOR;
1358  this->nevent += 1;
1359  this->eventID[0] = ID;
1360  this->eventID[1] = n;
1361  this->psi[0] = psi_net.data[k];
1362  this->phi[0] = psm->getPhi(ind); // reconstructed phi
1363  this->phi[2] = psm->getRA(ind);
1364  this->phi[3] = phi_net.data[k]; // detection phi
1365  this->theta[0] = psm->getTheta(ind); // reconstructed theta
1366  this->theta[2] = psm->getDEC(ind);
1367  this->theta[3] = theta_net.data[k]; // detection theta
1368  this->gnet = net->gNET;
1369  this->anet = net->aNET;
1370  this->inet = net->iNET;
1371  this->iota[0] = net->norm;
1372  this->norm = net->norm;
1373  this->likelihood = 0.;
1374 
1375  TAU = psm->get(this->theta[0],this->phi[0]);
1376 
1377  M = 0; gC = 0.;
1378  this->strain[0] = 0.;
1379  this->penalty = 0.;
1380  tot_null = 0.;
1381 
1382  Etot = 0.;
1383  for(i=0; i<N; i++) Etot += energy_net.data[i*K+k];
1384 
1385  for(i=0; i<5; i++) { this->neted[i] = 0.; }
1386 
1387  this->lag[N] = p->shift;
1388 
1389  for(i=0; i<N; i++) // loop on detectors
1390  {
1391  pd = net->getifo(i);
1392  Aa = pd->antenna(this->theta[0],this->phi[0]); // antenna pattern
1393  m = i*K+k;
1394 
1395  if(i) { // take delays into account
1396  psm = &(net->getifo(i)->tau);
1397  time_net.data[m] += psm->get(this->theta[0],this->phi[0])-TAU;
1398  }
1399 
1400  this->type[0] = 1;
1401  this->rate[i] = int(rate_net.data[k]+0.1);
1402 
1403  this->gap[i] = 0.;
1404  this->lag[i] = pd->lagShift.data[n];
1405 
1406  this->volume[i] = int(volume_net.data[k]+0.5);
1407  this->size[i] = int(size_net.data[k]+0.5);
1408 
1409  this->snr[i] = energy_net.data[m];
1410  this->nill[i] = pd->xSNR-pd->sSNR;
1411  this->null[i] = pd->null;
1412  this->xSNR[i] = pd->xSNR;
1413  this->sSNR[i] = pd->sSNR;
1414  this->neted[0] += fabs(pd->xSNR-pd->sSNR);
1415  this->neted[1] += fabs(pd->ED[1]);
1416  this->neted[2] += fabs(pd->ED[2]);
1417  this->neted[3] += fabs(pd->ED[3]);
1418  this->neted[4] += fabs(pd->ED[4]);
1419  this->likelihood += snr[i] - pd->null;
1420 
1421  this->time[i] = time_net.data[m] + this->gps[i];
1422  this->left[i] = start_net.data[m];
1423  this->right[i] = p->stop-p->start-stop_net.data[m];
1424  this->duration[i] = stop_net.data[m] - start_net.data[m];
1425  this->start[i] = start_net.data[m] + this->gps[i];
1426  this->stop[i] = stop_net.data[m] + this->gps[i];
1427 
1428  this->frequency[i] = frequency_net.data[m];
1429  this->low[i] = low_net.data[k];
1430  this->high[i] = high_net.data[k];
1431  this->bandwidth[i] = high_net.data[k] - low_net.data[k];
1432 
1433  if(net->MRA) {
1434  this->hrss[i] = pow(pow(10.,hrss_net.data[m]),2)+pow(pow(10.,hrss_NET.data[m]),2);
1435  this->hrss[i] = sqrt(this->hrss[i]/inRate);
1436  } else {
1437  this->hrss[i] = pow(10.,hrss_net.data[m])/sqrt(inRate);
1438  }
1439  this->noise[i] = pow(10.,noise_net.data[m])/sqrt(inRate);
1440  this->bp[i]= Aa.real();
1441  this->bx[i]= Aa.imag();
1442  this->strain[0]+= this->hrss[i]*this->hrss[i];
1443 
1444  tot_null += this->null[i];
1445  Aa /= pow(10.,NOISE_net.data[m]);
1446  gC += Aa*Aa;
1447 
1448  skSNR.data[i] = pd->sSNR;
1449  xkSNR.data[i] = pd->xSNR;
1450 
1451  x = 1. - pd->sSNR/(this->snr[i]+net->precision*Etot);
1452  if(x<this->penalty) this->penalty = x;
1453 
1454  psm->gps = time_net.data[m]+this->gps[0];
1455 
1456  }
1457 
1458  ndmLike = 0.;
1459  this->ECOR = 0.;
1460  this->netcc[2] = 0.;
1461  this->netcc[3] = 0.;
1462 
1463  for(i=0; i<N; i++) { // loop on detectors
1464  for(j=i; j<int(N); j++) {
1465  ndm[M] = net->getNDM(i,j);
1466 
1467  if(i!=j) {
1468  ndm[M] *= 2.;
1469  pcc = 2*sqrt(net->getNDM(i,i)*net->getNDM(j,j));
1470  pcc = pcc>0. ? ndm[M]/pcc : 0.;
1471  this->ECOR += ndm[M]*fabs(pcc); // effective correlated energy
1472  x = skSNR.data[i] + skSNR.data[j];
1473  this->netcc[2] += x*pcc;
1474  x = xkSNR.data[i] + xkSNR.data[j];
1475  this->netcc[3] += x*pcc;
1476  }
1477 
1478  ndmLike += ndm[M++];
1479  }
1480  }
1481 
1482  if(ndmLike>0.) { // temporary 2G condition
1483  ind = p->sArea[ID-1].size();
1484  for(i=0; i<11; i++)
1485  this->erA[i] = i<ind ? p->sArea[ID-1][i] : 0.;
1486 
1487  x = this->ecor;
1488  this->netcc[0] = x/(tot_null+x); // network (ecor) correlation
1489 
1490  x = this->ECOR;
1491  this->netcc[1] = x/(tot_null+x); // network (ECOR) correlation
1492 
1493  this->netcc[2] /= (N-1)*ndmLike; // network Pearson's correlation
1494  this->netcc[3] /= (N-1)*ndmLike; // network Pearson's correlation
1495  }
1496 
1497  this->rho[0] = sqrt(ecor*this->netcc[0]/N);
1498  this->rho[1] = sqrt(ECOR*this->netcc[0]/N);
1499 
1500  this->strain[0] = sqrt(this->strain[0]);
1501  this->penalty = sqrt(1./(1-this->penalty));
1502  if(!ellips) this->psi[0] = gC.arg()*180/Pi;
1503 
1504 // set injections if there are any
1505 
1506  M = net->mdc__IDSize();
1507  if(!n) { // only for zero lag
1508  injTime = 1.e12;
1509  injID = -1;
1510  for(m=0; m<M; m++) {
1511  mdcID = net->getmdc__ID(m);
1512  T = fabs(this->time[0] - net->getmdcTime(mdcID));
1513  if(T<injTime && INJ.fill_in(net,mdcID)) {
1514  injTime = T;
1515  injID = mdcID;
1516  }
1517 // printf("%d %12.3f %12.3f\n",mdcID,net->getmdcTime(mdcID),T);
1518  }
1519 
1520  if(INJ.fill_in(net,injID)) { // set injections
1521  this->range[1] = INJ.distance/factor;
1522  this->chirp[0] = INJ.mchirp;
1523  this->eBBH[1] = INJ.e0;
1524  this->eBBH[2] = INJ.rp0;
1525  this->eBBH[3] = INJ.redshift;
1526  this->strain[1] = INJ.strain*factor;
1527  this->type[1] = INJ.type;
1528  *this->name = net->getmdcType(this->type[1]-1);
1529  this->theta[1] = INJ.theta[0];
1530  this->phi[1] = INJ.phi[0];
1531  this->psi[1] = INJ.psi[0];
1532  this->iota[1] = INJ.iota[0];
1533  this->mass[0] = INJ.mass[0];
1534  this->mass[1] = INJ.mass[1];
1535 
1536  for(i=0; i<6; i++) this->spin[i] = INJ.spin[i];
1537 
1538 // printf("injection type: %d %12.3f\n",INJ.type,INJ.time[0]);
1539 
1540  wavearray<double>** pwfINJ = new wavearray<double>*[N];
1541  wavearray<double>** pwfREC = new wavearray<double>*[N];
1542  pd = net->getifo(0);
1543  int idSize = pd->RWFID.size();
1544  int wfIndex=-1;
1545  for (int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
1546 
1547  for(j=0; j<int(N); j++) {
1548  pd = net->getifo(j);
1549  Aa = pd->antenna(this->theta[1],this->phi[1]); // inj antenna pattern
1550  this->hrss[j+N] = INJ.hrss[j]*factor;
1551  this->bp[j+N] = Aa.real();
1552  this->bx[j+N] = Aa.imag();
1553  this->time[j+N] = INJ.time[j];
1554  this->Deff[j] = INJ.Deff[j]/factor;
1555 
1556  pwfINJ[j] = INJ.pwf[j];
1557  if (pwfINJ[j]==NULL) {
1558  cout << "Error : Injected waveform not saved !!! : detector "
1559  << net->ifoName[j] << endl;
1560  continue;
1561  }
1562  if (wfIndex<0) {
1563  cout << "Error : Reconstructed waveform not saved !!! : ID -> "
1564  << ID << " : detector " << net->ifoName[j] << endl;
1565  continue;
1566  }
1567 
1568  if (wfIndex>=0) pwfREC[j] = pd->RWFP[wfIndex];
1569  double R = pd->TFmap.rate();
1570  //double rFactor = log(2.)/log(pd->rate/R); // rescale waveform
1571  double rFactor = 1.;
1572  rFactor *= factor;
1573  wavearray<double>* wfINJ = pwfINJ[j];
1574  *wfINJ*=rFactor;
1575  wavearray<double>* wfREC = pwfREC[j];
1576 
1577  double bINJ = wfINJ->start();
1578  double eINJ = wfINJ->start()+wfINJ->size()/R;
1579  double bREC = wfREC->start();
1580  double eREC = wfREC->start()+wfREC->size()/R;
1581  //cout.precision(14);
1582  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
1583  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
1584 
1585  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
1586  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
1587  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
1588 
1589  double startXCOR = bINJ>bREC ? bINJ : bREC;
1590  double endXCOR = eINJ<eREC ? eINJ : eREC;
1591  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1592  //cout << "startXCOR : " << startXCOR << " endXCOR : " << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
1593 
1594  if (sizeXCOR<=0) continue;
1595 
1596  // the enINJ, enREC, xcorINJ_REC are computed in the INJ range
1597 
1598  double enINJ=0;
1599  for (int i=0;i<wfINJ->size();i++) enINJ+=wfINJ->data[i]*wfINJ->data[i];
1600  //for (int i=0;i<sizeXCOR;i++) enINJ+=wfINJ->data[i+oINJ]*wfINJ->data[i+oINJ];
1601 
1602  double enREC=0;
1603  //for (int i=0;i<wfREC->size();i++) enREC+=wfREC->data[i]*wfREC->data[i];
1604  for (int i=0;i<sizeXCOR;i++) enREC+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
1605 
1606  double xcorINJ_REC=0;
1607  for (int i=0;i<sizeXCOR;i++) xcorINJ_REC+=wfINJ->data[i+oINJ]*wfREC->data[i+oREC];
1608 
1609  WSeries<double> wfREC_SUB_INJ;
1610  wfREC_SUB_INJ.resize(sizeXCOR);
1611  for (int i=0;i<sizeXCOR;i++) wfREC_SUB_INJ.data[i]=wfREC->data[i+oREC]-wfINJ->data[i+oINJ];
1612  wfREC_SUB_INJ.start(startXCOR);
1613  wfREC_SUB_INJ.rate(wfREC->rate());
1614  this->iSNR[j] = enINJ;
1615  this->oSNR[j] = enREC;
1616  this->ioSNR[j] = xcorINJ_REC;
1617 
1618  //double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
1619  //cout << "enINJ : " << enINJ << " enREC : " << enREC << " xcorINJ_REC : " << xcorINJ_REC << endl;
1620  //cout << "erINJ_REC/enINJ : " << erINJ_REC/enINJ << endl;
1621 
1622  *wfINJ*=1./rFactor;
1623  }
1624  delete [] pwfINJ;
1625  delete [] pwfREC;
1626  }
1627  else { // no injections
1628  this->range[1] = 0.;
1629  this->chirp[0] = 0.;
1630  this->eBBH[1] = 0.;
1631  this->eBBH[2] = 0.;
1632  this->eBBH[3] = 0.;
1633  this->strain[1] = 0.;
1634  this->type[1] = 0;
1635  this->theta[1] = 0.;
1636  this->phi[1] = 0.;
1637  this->psi[1] = 0.;
1638  this->iota[1] = 0.;
1639  this->mass[0] = 0.;
1640  this->mass[1] = 0.;
1641 
1642  for(i=0; i<6; i++) this->spin[i] = 0.;
1643 
1644  for(j=0; j<int(N); j++) {
1645  this->hrss[j+N] = 0.;
1646  this->bp[j+N] = 0.;
1647  this->bx[j+N] = 0.;
1648  this->time[j+N] = 0.;
1649  this->Deff[j] = 0.;
1650  this->iSNR[j] = 0.;
1651  this->oSNR[j] = 0.;
1652  this->ioSNR[j] = 0.;
1653  }
1654  }
1655  }
1656 
1657  if(this->fP!=NULL) {
1658  fprintf(fP,"\n# trigger %d in lag %d for \n",int(ID),int(n));
1659  this->Dump("1G");
1660  vP = &(net->wc_List[n].p_Map[ID-1]);
1661  vI = &(net->wc_List[n].p_Ind[ID-1]);
1662  x = cos(psm->theta_1*PI/180.)-cos(psm->theta_2*PI/180.);
1663  x*= (psm->phi_2-psm->phi_1)*180/PI/psm->size();
1664  fprintf(fP,"sky_res: %f\n",sqrt(fabs(x)));
1665  fprintf(fP,"map_lenght: %d\n",int(vP->size()));
1666  fprintf(fP,"#skyID theta DEC step phi R.A step probability cumulative\n");
1667  x = 0.;
1668  for(j=0; j<int(vP->size()); j++) {
1669  i = (*vI)[j];
1670  if(net->mdc__IDSize()) { // simulation mode
1671  x = (j==int(vP->size())-1) ? 0 : x+(*vP)[j]; // last value is the inj sky index (x=0)
1672  } else x+=(*vP)[j];
1673  fprintf(fP,"%6d %5.1f %5.1f %6.2f %5.1f %5.1f %6.2f %e %e\n",
1674  int(i),psm->getTheta(i),psm->getDEC(i),psm->getThetaStep(i),
1675  psm->getPhi(i),psm->getRA(i),psm->getPhiStep(i),(*vP)[j],x);
1676  }
1677  }
1678 
1679  // save the probability skymap to tree
1680  if(waveTree!=NULL && net->wc_List[n].p_Map.size()) {
1681 
1682  vP = &(net->wc_List[n].p_Map[ID-1]);
1683  vI = &(net->wc_List[n].p_Ind[ID-1]);
1684  if(this->Psave) {
1685 
1686  *Psm = *psm;
1687  *Psm = 0.;
1688 
1689  int k;
1690  double th,ra;
1691  for(j=0; j<int(vP->size()); j++) {
1692  i = (*vI)[j];
1693  th = Psm->getTheta(i);
1694  ra = Psm->getRA(i);
1695  k=Psm->getSkyIndex(th, ra);
1696  Psm->set(k,(*vP)[j]);
1697  }
1698  }
1699  }
1700 
1701  if(waveTree!=NULL) waveTree->Fill();
1702 
1703  }
1704  }
1705  if(ndm) free(ndm);
1706 }
1707 
1708 
TTree * tree
Definition: TimeSortTree.C:20
std::vector< char * > ifoName
Definition: network.hh:609
Float_t anet
Definition: netevent.hh:118
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
Float_t * mass
effective range for each detector
Definition: netevent.hh:131
size_t append(const wavearray< DataType_t > &)
Definition: wavearray.cc:793
float SUBNET
Definition: netcluster.hh:66
double rho
double stop
Definition: netcluster.hh:380
virtual void resize(unsigned int)
Definition: wseries.cc:901
Float_t * eBBH
chirp array: 0-injmass,1-recmass,2-merr,3-tmrgr,4-terr,5-chi2
Definition: netevent.hh:129
Double_t * time
beam pattern coefficients for hx
Definition: injection.hh:89
Int_t ndim
current Tree number in a TChain
Definition: netevent.hh:65
Float_t rp0
Definition: injection.hh:76
size_t nLag
Definition: network.hh:573
Float_t * rho
biased null statistics
Definition: netevent.hh:112
Float_t * high
min frequency
Definition: netevent.hh:104
Float_t e0
Definition: injection.hh:77
double M
Definition: DrawEBHH.C:13
Float_t * range
Definition: netevent.hh:127
double duration
std::vector< netcluster > wc_List
Definition: network.hh:610
Float_t ecor
Definition: netevent.hh:120
void setSLags(float *slag)
Definition: netevent.cc:426
Bool_t Notify()
Definition: netevent.cc:325
Double_t * start
cluster duration = stopW-startW
Definition: netevent.hh:99
std::vector< double > * getmdcTime()
Definition: network.hh:422
double arg() const
Definition: wavecomplex.hh:71
float likenet
Definition: netcluster.hh:55
float gNET
Definition: netcluster.hh:77
virtual void rate(double r)
Definition: wavearray.hh:141
Float_t * duration
max cluster time relative to segment start
Definition: netevent.hh:98
Float_t * low
average center_of_snr frequency
Definition: netevent.hh:103
float factor
CWB run(runID)
double sSNR
Definition: detector.hh:338
wavearray< double > a(hp.size())
par [0] name
float energy
Definition: netcluster.hh:53
int n
Definition: cwb_net.C:28
TH1F * ecor
double imag() const
Definition: wavecomplex.hh:70
TTree * setTree()
Definition: netevent.cc:434
TString("c")
string * log
event name: "" - prod, mdc_name - sim
Definition: netevent.hh:76
float netED
Definition: netcluster.hh:59
size_t nRun
Definition: network.hh:572
float chirpEfrac
Definition: netcluster.hh:87
Float_t inet
Definition: netevent.hh:119
int ID
Definition: TestMDC.C:70
Int_t * size
cluster volume
Definition: netevent.hh:80
double gps
Definition: skymap.hh:325
Float_t * right
segment start GPS time
Definition: netevent.hh:96
float mchirp
Definition: netcluster.hh:84
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
double frequency
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:2207
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:382
Int_t GetEntry(Int_t)
Definition: netevent.cc:409
Float_t ECOR
Definition: netevent.hh:122
float chirpPfrac
Definition: netcluster.hh:88
float theta
Float_t * gap
cluster union size
Definition: netevent.hh:83
float normcor
Definition: netcluster.hh:57
std::vector< vector_float > sArea
Definition: netcluster.hh:401
netcluster * pwc
Definition: cwb_job_obj.C:38
TH2F * ph
string getmdcType(size_t n)
Definition: network.hh:420
double get_SS()
Definition: detector.hh:309
Float_t * left
min cluster time relative to segment start
Definition: netevent.hh:97
double getTheta(size_t i)
Definition: skymap.hh:224
Float_t * ioSNR
reconstructed snr waveform
Definition: netevent.hh:139
Float_t * iota
source psi angle
Definition: injection.hh:83
double & gap
string getmdcList(size_t n)
Definition: network.hh:418
double getThetaStep(size_t i)
Definition: skymap.hh:242
double theta_1
Definition: skymap.hh:321
Long_t size
TTree * Init(TString fName, int n)
Definition: netevent.cc:47
Int_t * type
event ID
Definition: netevent.hh:74
int m
Definition: cwb_net.C:28
double getPhiStep(size_t i)
Definition: skymap.hh:182
std::vector< double > oSNR[NIFO_MAX]
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
double rate
Definition: detector.hh:341
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
Definition: network.cc:3666
Double_t strain
Definition: injection.hh:81
i drho i
skymap tau
Definition: detector.hh:346
std::vector< detector * > ifoList
Definition: network.hh:608
std::vector< double > ioSNR[NIFO_MAX]
#define N
Float_t * theta
source phi angle
Definition: injection.hh:85
double Edge
Definition: network.hh:578
network ** net
NOISE_MDC_SIMULATION.
double theta_2
Definition: skymap.hh:322
size_t ifoListSize()
Definition: network.hh:431
Int_t run
max size used by allocate() for the probability maps
Definition: netevent.hh:71
int Psave
Definition: test_config1.C:75
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
#define PI
Definition: watfun.hh:32
float likesky
Definition: netcluster.hh:61
double aNET
Definition: network.hh:580
bool optim
Definition: network.hh:590
#define nIFO
virtual size_t size() const
Definition: wavearray.hh:145
float phi
Float_t * frequency
GPS stop time of the cluster.
Definition: netevent.hh:102
double precision
Definition: network.hh:593
double ra
Definition: ConvertGWGC.C:46
Float_t redshift
Definition: injection.hh:78
Float_t factor
Definition: netevent.hh:126
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
Definition: detector.cc:490
double hrss
Definition: TestMDC.C:70
float psi
void init()
Definition: netevent.cc:283
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:711
Float_t * psi
Definition: injection.hh:82
Int_t nevent
run ID
Definition: netevent.hh:72
double shift
Definition: netcluster.hh:382
Int_t Psave
number of detectors
Definition: netevent.hh:66
Double_t * hrss
estimated bandwidth
Definition: injection.hh:95
#define WAVE_TREE_NAME
Definition: netevent.cc:41
Int_t * volume
1/rate - wavelet time resolution
Definition: netevent.hh:79
double real() const
Definition: wavecomplex.hh:69
double ED[5]
Definition: detector.hh:334
Float_t * Deff
eBBH array
Definition: netevent.hh:130
Double_t * gps
average center_of_gravity time
Definition: netevent.hh:95
double getRA(size_t i)
Definition: skymap.hh:215
i() int(T_cor *100))
double Pi
Float_t mchirp
Definition: injection.hh:74
bool log
Definition: WaveMDC.C:41
Float_t penalty
Definition: netevent.hh:123
Double_t * stop
GPS start time of the cluster.
Definition: netevent.hh:100
float enrgsky
Definition: netcluster.hh:54
double eCOR
Definition: network.hh:582
Float_t gnet
network energy disbalance: 0 - total, 1 - 00-phase, 2 - 90-phase
Definition: netevent.hh:117
double start
Definition: netcluster.hh:379
double gNET
Definition: network.hh:579
double getDEC(size_t i)
Definition: skymap.hh:251
Float_t * lag
time between consecutive events
Definition: netevent.hh:84
std::vector< int > RWFID
Definition: detector.hh:381
float netcc
Definition: netcluster.hh:63
Float_t * mass
detector specific effective distance
Definition: injection.hh:98
Float_t * xSNR
energy/noise_variance
Definition: netevent.hh:135
size_t mdc__IDSize()
Definition: network.hh:414
std::vector< double > iSNR[NIFO_MAX]
std::vector< detector * > ifoList
dump file
Definition: netevent.hh:143
void Show(Int_t entry=-1)
Definition: netevent.cc:415
int k
int pattern
Definition: network.hh:601
float spin[6]
Double_t * strain
time slag [sec]
Definition: netevent.hh:86
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:113
void Dump(WSeries< double > &w, double t1, double t2, const char *fname)
Float_t * slag
time lag [sec]
Definition: netevent.hh:85
void output(TTree *=NULL, network *=NULL, double=0., size_t=0, int=-1)
Definition: netevent.cc:1193
Definition: skymap.hh:63
wavearray< double > lagShift
Definition: detector.hh:369
double * entry
Definition: cwb_setcuts.C:224
Float_t * phi
source iota angle
Definition: injection.hh:84
Float_t * iota
[0]-reconstructed psi or phase of gc, [1]-injected psi angle
Definition: netevent.hh:90
size_t size()
Definition: netcluster.hh:147
bool USE_HEALPIX()
Definition: wat.hh:31
Float_t * psi
[0]-reconstructed, [1]-injected theta angle, [2]-DEC
Definition: netevent.hh:89
Double_t * hrss
high-low
Definition: netevent.hh:106
Int_t type
Definition: injection.hh:68
float iota
Definition: netcluster.hh:72
int * xstart
TH1F * ECOR
float norm
Definition: netcluster.hh:80
FILE * fP
float aNET
Definition: netcluster.hh:78
float mchirperr
Definition: netcluster.hh:85
double getPhi(size_t i)
Definition: skymap.hh:164
WSeries< double > TFmap
Definition: detector.hh:354
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:87
Float_t norm
Definition: netevent.hh:121
double gps
float iNET
Definition: netcluster.hh:79
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:88
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:94
double T
Definition: testWDM_4.C:11
Int_t * rate
event log: "" - prod, mdc_log - sim
Definition: netevent.hh:77
float netRHO
Definition: netcluster.hh:68
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:179
Float_t distance
Definition: injection.hh:73
double xSNR
Definition: detector.hh:339
double getNDM(size_t i, size_t j)
param: first detector param: second detector
Definition: network.hh:140
double fabs(const Complex &x)
Definition: numpy.cc:55
Float_t * nill
un-biased null statistics
Definition: netevent.hh:111
Float_t * spin
mass[2], binary mass parameters
Definition: netevent.hh:132
Float_t * sSNR
data-signal correlation Xk*Sk
Definition: netevent.hh:136
double norm
Definition: network.hh:583
void allocate()
Definition: netevent.cc:167
Float_t * erA
noise rms
Definition: netevent.hh:108
Float_t * null
probability cc skymap
Definition: netevent.hh:110
Float_t likelihood
Definition: netevent.hh:124
Definition: Toolbox.hh:99
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
Float_t * bx
beam pattern coefficients for hp
Definition: netevent.hh:92
float subnet
Definition: netcluster.hh:65
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
float cFreq
Definition: netcluster.hh:76
float skycc
Definition: netcluster.hh:62
std::vector< clusterdata > cData
Definition: netcluster.hh:391
virtual netevent & operator=(const netevent &)
Definition: netevent.cc:610
float cTime
Definition: netcluster.hh:75
float netnull
Definition: netcluster.hh:58
Double_t * noise
hrss
Definition: netevent.hh:107
double phi_2
Definition: skymap.hh:324
string * name
event type: [0] - prod, [1]-sim
Definition: netevent.hh:75
DataType_t * data
Definition: wavearray.hh:319
Float_t * snr
spin[6], binary spin parameters
Definition: netevent.hh:134
double null
Definition: detector.hh:336
Long_t id
float nDoF
Definition: netcluster.hh:81
double phi_1
Definition: skymap.hh:323
double get(size_t i)
param: sky index
Definition: skymap.cc:699
float mass[2]
double ctime
Int_t GetEntries()
Definition: netevent.cc:403
Bool_t fill_in(network *, int, bool=true)
Definition: injection.cc:349
Float_t * neted
network correlation coefficients: 0-net,1-pc,2-cc,3-net2
Definition: netevent.hh:114
size_t size()
Definition: skymap.hh:136
size_t getmdc__ID(size_t n)
Definition: network.hh:426
char fName[256]
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
Definition: injection.hh:101
virtual void resize(unsigned int)
Definition: wavearray.cc:463
bool MRA
Definition: network.hh:599
char tYPe
Definition: network.hh:588
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:6525
Float_t * spin
[m1,m2], binary mass parameters
Definition: injection.hh:99
float chirpEllip
Definition: netcluster.hh:89
bool SETNDM(size_t, size_t, bool=true, int=1)
Definition: network.cc:6824
float netecor
Definition: netcluster.hh:56
Float_t * chirp
range to source: [0/1]-rec/inj
Definition: netevent.hh:128
float theta
Definition: netcluster.hh:70
double TAU[nWF]
detector ** pD
double iNET
Definition: network.hh:581
skymap * Psm
error angle
Definition: netevent.hh:109
Float_t * Deff
injected snr in the detectors
Definition: injection.hh:97
float netrho
Definition: netcluster.hh:69
exit(0)
double enrg
Definition: detector.hh:337
Float_t * bp
[0]-reconstructed iota angle, [1]-injected iota angle
Definition: netevent.hh:91
Float_t * iSNR
energy of reconstructed responses Sk*Sk
Definition: netevent.hh:137
Int_t * eventID
event count
Definition: netevent.hh:73
Float_t * bandwidth
max frequency
Definition: netevent.hh:105
Int_t usize
cluster size (black pixels only)
Definition: netevent.hh:81
Float_t * oSNR
injected snr waveform
Definition: netevent.hh:138