Logo coherent WaveBurst  
Library Reference Guide 		Logo
gnetwork.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 
20 #include "gnetwork.hh"
21 #include "constants.hh"
22 
23 using namespace std;
24 using namespace ROOT::Math;
25 
26 ClassImp(gnetwork)
27 
28 
29 //______________________________________________________________________________
30 /* Begin_Html
31 <center><h2>gnetwork class</h2></center>
32 This class gnetwork is derived from the network class and uses the methods of the class gskymap
33 It is used to produce the skymap plots of the detector's network antenna pattern or the statistics produced in the cwb likelihood stage.
34 
35 <p>
36 <h3><a name="example">Example</a></h3>
37 <p>
38 The macro <a href="./tutorials/gwat/DrawAntennaPattern.C.html">DrawAntennaPattern.C</a> is an example which shown how to use the gnetwork class.<br>
39 The pictures below gives the macro output plot.<br>
40 <p>
41 
42 End_Html
43 
44 Begin_Macro
45 DrawAntennaPattern.C
46 End_Macro */
47 
48 
49 //______________________________________________________________________________
50 gnetwork::gnetwork(int nifo, TString* ifo, detectorParams* detParms) : network() {
51 //
52 // constructor
53 //
54 // Input: nifo - number of detectors
55 // ifo - array of built-in ifo names
56 // Ex : TString ifo[3] = {"L1","H1","V1"}
57 // detParms - array of usewr defined detector's parameters
58 // it is used when ifo entry is empty -> ""
59 //
60 
61  for(int n=0; n<nifo; n++) {
62  if((ifo!=NULL)&&(ifo[n]!="")) {
63  pD[n] = new detector(const_cast<char*>(ifo[n].Data())); // built in detector
64  } else {
65  if(detParms!=NULL) {
66  pD[n] = new detector(detParms[n]); // user define detector
67  } else {
68  cout << "gnetwork::gnetwork : Error - user define detector with no detParams" << endl;
69  exit(1);
70  }
71  }
72  }
73 
74  for(int n=0; n<nifo; n++) this->add(pD[n]);
75  Init();
76 }
77 
78 //______________________________________________________________________________
79 gnetwork::gnetwork(network* net) : network() {
80 //
81 // constructor
82 //
83 // Input: net - pointer to network
84 //
85 
86  int nifo = net->ifoListSize();
87  for(int n=0; n<nifo; n++) pD[n] = new detector(*net->getifo(n));
88 
89  for(int n=0; n<nifo; n++) this->add(pD[n]);
90  Init();
91 }
92 
93 //______________________________________________________________________________
94 gnetwork::~gnetwork() {
95 //
96 // destructor
97 //
98 
99  ClearSites();
100  ClearSitesArms();
101  ClearSitesLabel();
102  ClearCircles();
103  for(int n=0;n<NIFO_MAX;n++) if(pD[n]!=NULL) delete pD[n];
104  for(int i=0;i<2;i++) if(grP[i]!=NULL) delete grP[i];
105  if(canvasP!=NULL) delete canvasP;
106  for(size_t n=0; n<this->ifoList.size(); n++) if(this->ifoList[n]!=NULL) delete this->ifoList[n];
107 }
108 
109 //______________________________________________________________________________
110 void
111 gnetwork::Init() {
112 //
113 // initialize internal parameters
114 //
115 
116  siteLabelType=0;
117  for(int n=0;n<NIFO_MAX;n++) pD[n]=NULL;
118  for(int n=0;n<NIFO_MAX;n++) scale[n]=1;
119  for(int i=0;i<2;i++) grP[i]=NULL;
120  canvasP=NULL;
121 }
122 
123 //______________________________________________________________________________
124 void
125 gnetwork::Draw(TString smName, network* net) {
126 //
127 // Draw skymap network info
128 //
129 // Input: smName - skymap name
130 //
131 // nSensitivity : network sensitivity
132 // nAlignment : network alignment factor
133 // nCorrelation : network correlation coefficient
134 // nLikelihood : network likelihood
135 // nNullEnergy : network null energy
136 // nPenalty : signal * noise penalty factor
137 // nCorrEnergy : reduced correlated energy
138 // nNetIndex : network index
139 // nDisbalance : energy disbalance
140 // nSkyStat : sky optimization statistic
141 // nEllipticity : waveform ellipticity
142 // nPolarisation : polarisation angle
143 // nProbability : probability skymap
144 //
145 // index : theta, phi mask index array
146 // skyMask : index array for setting sky mask
147 // skyMaskCC : index array for setting sky mask Celestial Coordinates
148 // skyHole : static sky mask describing "holes"
149 // veto : veto array for pixel selection
150 // skyProb : sky probability
151 // skyENRG : energy skymap
152 //
153 // net - pointer to network object
154 // if net!=NULL then net is used to retrive the data to be plotted
155 //
156 
157  network* NET = net!=NULL ? net : this;
158 
159  if(smName=="nSensitivity") {Draw(NET->nSensitivity,smName);return;} // network sensitivity
160  if(smName=="nAlignment") {Draw(NET->nAlignment,smName);return;} // network alignment factor
161  if(smName=="nCorrelation") {Draw(NET->nCorrelation,smName);return;} // network correlation coefficient
162  if(smName=="nLikelihood") {Draw(NET->nLikelihood,smName);return;} // network likelihood
163  if(smName=="nNullEnergy") {Draw(NET->nNullEnergy,smName);return;} // network null energy
164  if(smName=="nPenalty") {Draw(NET->nPenalty,smName);return;} // signal * noise penalty factor
165  if(smName=="nCorrEnergy") {Draw(NET->nCorrEnergy,smName);return;} // reduced correlated energy
166  if(smName=="nNetIndex") {Draw(NET->nNetIndex,smName);return;} // network index
167  if(smName=="nDisbalance") {Draw(NET->nDisbalance,smName);return;} // energy disbalance
168  if(smName=="nSkyStat") {Draw(NET->nSkyStat,smName);return;} // sky optimization statistic
169  if(smName=="nEllipticity") {Draw(NET->nEllipticity,smName);return;} // waveform ellipticity
170  if(smName=="nPolarisation") {Draw(NET->nPolarisation,smName);return;} // polarisation angle
171  if(smName=="nProbability") {Draw(NET->nProbability,smName);return;} // probability skymap
172 
173  if(smName=="index") {Draw(NET->index,smName);return;} // theta, phi mask index array
174  if(smName=="skyMask") {Draw(NET->skyMask,smName);return;} // index array for setting sky mask
175  if(smName=="skyMaskCC") {Draw(NET->skyMaskCC,smName);return;} // index array for setting sky mask Celestial Coordinates
176  if(smName=="skyHole") {Draw(NET->skyHole,smName);return;} // static sky mask describing "holes"
177  if(smName=="veto") {Draw(NET->veto,smName);return;} // veto array for pixel selection
178  if(smName=="skyProb") {Draw(NET->skyProb,smName);return;} // sky probability
179  if(smName=="skyENRG") {Draw(NET->skyENRG,smName);return;} // energy skymap
180 
181  cout << "gnetwork::Draw : Warning - skymap '" << smName.Data() << "' not valid !" << endl;
182 }
183 
184 //______________________________________________________________________________
185 double
186 gnetwork::GetSite(TString ifo, TString type) {
187 //
188 // return detector infos
189 //
190 // Input: ifo - detector name
191 // type - detector info type
192 //
193 // THETA : theta full
194 // THETA_D : theta degrees
195 // THETA_M : theta minutes
196 // THETA_S : theta seconds
197 //
198 // PHI : phi full
199 // PHI_D : phi degrees
200 // PHI_M : phi minutes
201 // PHI_S : phi seconds
202 //
203 // "" : print out theta,phi infos
204 //
205 // if type invalid return 1111
206 //
207 
208  int nIFO = this->ifoListSize();
209  if(nIFO<1) {cout << "gnetwork::GetSite nIFO must be >= 1" << endl;return -1111;}
210 
211  detector *pD[NIFO];
212  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
213 
214  int ifoId=0;
215  int nifo=0;
216  for(int n=0; n<nIFO; n++) {
217  if(ifo.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId=n;}
218  }
219  if(nifo<1) {cout << "gnetwork::GetSite ifo not valid" << endl;return -1111;}
220 
221  detectorParams dP = pD[ifoId]->getDetectorParams();
222 
223  double phi=dP.longitude;
224  double theta=dP.latitude;
225 
226  char LAT;
227  double theta_t=theta;
228  if(theta_t>0) LAT='N'; else {LAT='S';theta_t=-theta_t;}
229  int theta_d = int(theta_t);
230  int theta_m = int((theta_t-theta_d)*60);
231  float theta_s = (theta_t-theta_d-theta_m/60.)*3600.;
232 
233  char LON;
234  double phi_t=phi;
235  if(phi_t>0) LON='E'; else {LON='W';phi_t=-phi_t;}
236  int phi_d = int(phi_t);
237  int phi_m = int((phi_t-phi_d)*60);
238  float phi_s = (phi_t-phi_d-phi_m/60.)*3600.;
239 
240  type.ToUpper();
241  if(type.CompareTo("THETA")==0) return theta;
242  if(type.CompareTo("THETA_D")==0) return theta_d;
243  if(type.CompareTo("THETA_M")==0) return theta_m;
244  if(type.CompareTo("THETA_S")==0) return theta_s;
245  if(type.CompareTo("PHI")==0) return phi;
246  if(type.CompareTo("PHI_D")==0) return phi_d;
247  if(type.CompareTo("PHI_M")==0) return phi_m;
248  if(type.CompareTo("PHI_S")==0) return phi_s;
249 
250  if(type.CompareTo("")==0) pD[ifoId]->print();
251 
252  return 1111;
253 }
254 
255 //______________________________________________________________________________
256 void
257 gnetwork::DrawSites(Color_t mcolor, Size_t msize, Style_t mstyle) {
258 //
259 // Draw a marker in the ifo site position (see ROOT TMarker input parameters)
260 //
261 // Input: mcolor - color of marker
262 // msize - size of marker
263 // mstyle - style of marker
264 //
265 
266  if(gSM.goff) return;
267  if(gSM.canvas==NULL) return;
268 
269  int nIFO = this->ifoListSize();
270  if(nIFO<1) return;
271 
272  ClearSites();
273 
274  detector *pD[NIFO];
275  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
276 
277  gSM.canvas->cd();
278 
279  double r2d = 180./TMath::Pi();
280  for(int n=0;n<nIFO;n++) {
281  XYZVector Rv(pD[n]->Rv[0],pD[n]->Rv[1],pD[n]->Rv[2]);
282  TVector3 vRv(Rv.X(),Rv.Y(),Rv.Z());
283 
284  double phi=vRv.Phi()*r2d;
285  double theta=vRv.Theta()*r2d;
286 
287  theta=90-theta;
288 
289  if (gSM.projection.CompareTo("hammer")==0) {
290  gSM.ProjectHammer(phi, theta, phi, theta);
291  } else if (gSM.projection.CompareTo("sinusoidal")==0) {
292  gSM.ProjectSinusoidal(phi, theta, phi, theta);
293  } else if (gSM.projection.CompareTo("parabolic")==0) {
294  gSM.ProjectParabolic(phi, theta, phi, theta);
295  }
296  if (gSM.coordinate.CompareTo("CWB")==0) GeographicToCwb(phi,theta,phi,theta);
297  TMarker* pM = new TMarker(phi,theta,mstyle);
298  pM->SetMarkerSize(msize);
299  pM->SetMarkerColor(mcolor);
300  pM->Draw();
301  siteM.push_back(pM);
302  pM = new TMarker(phi,theta,mstyle);
303  pM->SetMarkerSize(msize/2);
304  pM->SetMarkerColor(kWhite);
305  pM->Draw();
306  siteM.push_back(pM);
307  }
308 
309  return;
310 }
311 
312 //______________________________________________________________________________
313 void
314 gnetwork::DrawSitesArms(double mlength, Color_t lcolor, Size_t lwidth, Style_t lstyle) {
315 //
316 // Draw detector's arms in the ifo site position (see ROOT TPolyLine input parameters)
317 //
318 // Input: mlength - arm's length in KM
319 // mcolor - color of line
320 // msize - size of line
321 // mstyle - style of line
322 //
323 
324 
325  if(gSM.goff) return;
326  if(gSM.canvas==NULL) return;
327 
328  int nIFO = this->ifoListSize();
329  if(nIFO<1) return;
330 
331  ClearSitesArms();
332 
333  detector *pD[NIFO];
334  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
335 
336  gSM.canvas->cd();
337 
338  double mlength_deg = 360.*mlength/(TMath::TwoPi()*watconstants::EarthEquatorialRadius());
339 
340  double x[nIFO][2][2];
341  double y[nIFO][2][2];
342  double phi,theta;
343  double r2d = 180./TMath::Pi();
344  for(int n=0;n<nIFO;n++) {
345  XYZVector Rv(pD[n]->Rv[0],pD[n]->Rv[1],pD[n]->Rv[2]);
346  XYZVector Ex(pD[n]->Ex[0],pD[n]->Ex[1],pD[n]->Ex[2]);
347  XYZVector Ey(pD[n]->Ey[0],pD[n]->Ey[1],pD[n]->Ey[2]);
348  TVector3 vRv(Rv.X(),Rv.Y(),Rv.Z());
349  TVector3 vEx(Ex.X(),Ex.Y(),Ex.Z());
350  TVector3 vEy(Ey.X(),Ey.Y(),Ey.Z());
351  vEx=mlength*vEx+vRv;
352  vEy=mlength*vEy+vRv;
353 
354  phi=vRv.Phi()*r2d;
355  theta=vRv.Theta()*r2d;
356  theta=90-theta;
357  if (gSM.projection.CompareTo("hammer")==0) {
358  gSM.ProjectHammer(phi, theta, phi, theta);
359  } else if (gSM.projection.CompareTo("sinusoidal")==0) {
360  gSM.ProjectSinusoidal(phi, theta, phi, theta);
361  } else if (gSM.projection.CompareTo("parabolic")==0) {
362  gSM.ProjectParabolic(phi, theta, phi, theta);
363  }
364  if (gSM.coordinate.CompareTo("CWB")==0) GeographicToCwb(phi,theta,phi,theta);
365  x[n][0][0]=phi;y[n][0][0]=theta;
366  x[n][1][0]=phi;y[n][1][0]=theta;
367 
368  phi=vEx.Phi()*r2d;
369  theta=vEx.Theta()*r2d;
370  theta=90-theta;
371  if (gSM.projection.CompareTo("hammer")==0) {
372  gSM.ProjectHammer(phi, theta, phi, theta);
373  } else if (gSM.projection.CompareTo("sinusoidal")==0) {
374  gSM.ProjectSinusoidal(phi, theta, phi, theta);
375  } else if (gSM.projection.CompareTo("parabolic")==0) {
376  gSM.ProjectParabolic(phi, theta, phi, theta);
377  }
378  if (gSM.coordinate.CompareTo("CWB")==0) GeographicToCwb(phi,theta,phi,theta);
379  x[n][0][1]=phi;y[n][0][1]=theta;
380  // fix draw on the edge
381  if (gSM.coordinate.CompareTo("CWB")==0) {
382  if((x[n][0][0]>360-mlength_deg)&&(x[n][0][0]+x[n][0][1]-360<mlength_deg)) x[n][0][1]=360;
383  }
384  if (gSM.coordinate.CompareTo("GEOGRAPHIC")==0) {
385  if((x[n][0][0]>180-mlength_deg)&&(x[n][0][0]+x[n][0][1]<mlength_deg)) x[n][0][1]=180;
386  }
387 
388  phi=vEy.Phi()*r2d;
389  theta=vEy.Theta()*r2d;
390  theta=90-theta;
391  if (gSM.projection.CompareTo("hammer")==0) {
392  gSM.ProjectHammer(phi, theta, phi, theta);
393  } else if (gSM.projection.CompareTo("sinusoidal")==0) {
394  gSM.ProjectSinusoidal(phi, theta, phi, theta);
395  } else if (gSM.projection.CompareTo("parabolic")==0) {
396  gSM.ProjectParabolic(phi, theta, phi, theta);
397  }
398  if (gSM.coordinate.CompareTo("CWB")==0) GeographicToCwb(phi,theta,phi,theta);
399  x[n][1][1]=phi;y[n][1][1]=theta;
400  // fix draw on the edge
401  if (gSM.coordinate.CompareTo("CWB")==0) {
402  if((x[n][1][0]>360-mlength_deg)&&(x[n][1][0]+x[n][1][1]-360<mlength_deg)) x[n][1][1]=360;
403  }
404  if (gSM.coordinate.CompareTo("GEOGRAPHIC")==0) {
405  if((x[n][1][0]>180-mlength_deg)&&(x[n][1][0]+x[n][1][1]<mlength_deg)) x[n][1][1]=180;
406  }
407  }
408 
409  // Draw X,Y arms
410  for(int n=0;n<nIFO;n++) {
411  for(int m=0;m<2;m++) {
412  TPolyLine *pL = new TPolyLine(2,x[n][m],y[n][m]);
413  pL->SetLineColor(lcolor);
414  pL->SetLineStyle(lstyle);
415  pL->SetLineWidth(lwidth);
416  pL->Draw();
417  siteA.push_back(pL);
418  }
419  }
420 
421  return;
422 }
423 
424 //______________________________________________________________________________
425 void
426 gnetwork::DrawSitesShortLabel(Color_t tcolor, Size_t tsize, Font_t tfont) {
427 //
428 // Draw detector's short label in the ifo site position (see ROOT TLatex input parameters)
429 // Ex : L,H,V
430 //
431 // Input: tcolor - color of text
432 // tsize - size of text
433 // tfont - font of text
434 //
435 
436  siteLabelType=1;
437  DrawSitesLabel(tcolor, tsize, tfont);
438 }
439 
440 //______________________________________________________________________________
441 void
442 gnetwork::DrawSitesLabel(Color_t tcolor, Size_t tsize, Font_t tfont) {
443 //
444 // Draw detector's short label in the ifo site position (see ROOT TLatex input parameters)
445 // Ex : L1,H1,V1
446 //
447 // Input: tcolor - color of text
448 // tsize - size of text
449 // tfont - font of text
450 //
451 
452  if(gSM.goff) return;
453  if(gSM.canvas==NULL) return;
454 
455  int nIFO = this->ifoListSize();
456  if(nIFO<1) return;
457 
458  ClearSitesLabel();
459 
460  detector *pD[NIFO];
461  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
462 
463  gSM.canvas->cd();
464 
465  double r2d = 180./TMath::Pi();
466  for(int n=0;n<nIFO;n++) {
467  XYZVector Rv(pD[n]->Rv[0],pD[n]->Rv[1],pD[n]->Rv[2]);
468  TVector3 vRv(Rv.X(),Rv.Y(),Rv.Z());
469 
470  double phi=vRv.Phi()*r2d;
471  double theta=vRv.Theta()*r2d;
472 
473  theta=90-theta;
474 
475  if (gSM.projection.CompareTo("hammer")==0) {
476  gSM.ProjectHammer(phi, theta, phi, theta);
477  } else if (gSM.projection.CompareTo("sinusoidal")==0) {
478  gSM.ProjectSinusoidal(phi, theta, phi, theta);
479  } else if (gSM.projection.CompareTo("parabolic")==0) {
480  gSM.ProjectParabolic(phi, theta, phi, theta);
481  }
482  if (gSM.coordinate.CompareTo("CWB")==0) GeographicToCwb(phi,theta,phi,theta);
483  TString ifoName=this->getifo(n)->Name;
484  if(ifoName.CompareTo("H2")==0) {phi+=2;theta-=10;}
485  else if(ifoName.CompareTo("J1")==0) {phi+=6;theta-=1;}
486  else if(ifoName.CompareTo("H1")==0) {phi+=6;theta-=1;}
487  else if(ifoName.CompareTo("L1")==0) {phi+=3;theta+=5;}
488  else if(ifoName.CompareTo("V1")==0) {phi+=4;theta+=2;}
489  else {phi+=2;theta+=2;}
490  if(siteLabelType==1) {
491  if(ifoName.CompareTo("A1")==0) ifoName=TString("#tilde{A}");
492  else if(ifoName.CompareTo("H2")==0) ifoName=TString("#tilde{H}");
493  else if(ifoName.CompareTo("I2")==0) ifoName=TString("#tilde{I}");
494  else if(ifoName.CompareTo("R2")==0) ifoName=TString("#tilde{R}");
495  else ifoName.Resize(1);
496  }
497  TLatex *pL = new TLatex(phi+1.0,theta, ifoName);
498  pL->SetTextFont(tfont);
499  pL->SetTextSize(tsize);
500  pL->SetTextColor(tcolor);
501  pL->Draw();
502  siteL.push_back(pL);
503  }
504 
505  return;
506 }
507 
508 //______________________________________________________________________________
509 void
510 gnetwork::ClearCircles() {
511 //
512 // Clear all circles
513 //
514 
515  for(int i=0; i<(int)circleL.size(); i++) {if(circleL[i]) delete circleL[i];}
516  circleL.clear();
517  for(int i=0; i<(int)daxisM.size(); i++) {if(daxisM[i]) delete daxisM[i];}
518  daxisM.clear();
519  return;
520 }
521 
522 //______________________________________________________________________________
523 void
524 gnetwork::ClearSites() {
525 //
526 // Clear all site markers
527 //
528 
529  for(int i=0; i<(int)siteM.size(); i++) {if(siteM[i]) delete siteM[i];}
530  siteM.clear();
531  return;
532 }
533 
534 //______________________________________________________________________________
535 void
536 gnetwork::ClearSitesArms() {
537 //
538 // Clear all site arms
539 //
540 
541  for(int i=0; i<(int)siteA.size(); i++) {if(siteA[i]) delete siteA[i];}
542  siteA.clear();
543  return;
544 }
545 
546 //______________________________________________________________________________
547 void
548 gnetwork::ClearSitesLabel() {
549 //
550 // Clear all site labels
551 //
552 
553  for(int i=0; i<(int)siteT.size(); i++) {if(siteT[i]) delete siteT[i];}
554  siteT.clear();
555  for(int i=0; i<(int)siteL.size(); i++) {if(siteL[i]) delete siteL[i];}
556  siteL.clear();
557  for(int i=0; i<(int)siteP.size(); i++) {if(siteP[i]) delete siteP[i];}
558  siteP.clear();
559  return;
560 }
561 
562 //______________________________________________________________________________
563 double
564 gnetwork::GetAntennaPattern(double phi, double theta, double psi, int polarization) {
565 //
566 // Get Network Antenna Pattern value in the DPF : Dominant Polarization Frame
567 //
568 // Input: phi - phi (degrees)
569 // theta - theta (degrees)
570 // psi - polarization angle (degrees)
571 // polarization - DPF component
572 //
573 // 0 -> |Fx| - DPF
574 // 1 -> |F+| - DPF
575 // 2 -> |Fx|/|F+| - DPF
576 // 3 -> sqrt((|F+|^2+|Fx|^2)/nIFO) - DPF
577 // 4 -> |Fx|^2 - DPF
578 // 5 -> |F+|^2 - DPF
579 // 6 -> Fx - only with 1 detector
580 // 7 -> F+ - only with 1 detector
581 // 8 -> F1x/F2x - only with 2 detectors
582 // 9 -> F1+/F2+ - only with 2 detectors
583 // 10 -> sqrt(|F1+|^2+|F1x|^2)/sqrt(|F2+|^2+|F2x|^2) - only with 2 detectors
584 // 12 -> DPF Angle
585 //
586 
587  if((polarization<0||polarization>10)&&(polarization!=12))
588  {cout << "polarization must be [0-10][12]" << endl;exit(1);}
589 
590  int nIFO = this->ifoListSize();
591  if(nIFO==0)
592  {cout << "detectors are no declared" << endl;exit(1);}
593  if((nIFO!=1)&&(polarization==6||polarization==7))
594  {cout << "with polarization [6-7] works only with one detector" << endl;exit(1);}
595  if((nIFO!=2)&&(polarization>=8&&polarization<=11))
596  {cout << "with polarization [8-11] works only two one detectors" << endl;exit(1);}
597 
598  detector *pD[NIFO];
599  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
600 
601  double z=0.;
602  std::vector<wavecomplex> F;
603  F.resize(nIFO);
604 
605  for (int n=0;n<nIFO;n++) F[n] = pD[n]->antenna(theta,phi,psi);
606  double gp=0,gx=0,gI=0;
607  for (int n=0;n<nIFO;n++) {
608  gp+=scale[n]*(F[n].real()*F[n].real());
609  gx+=scale[n]*(F[n].imag()*F[n].imag());
610  gI+=scale[n]*(F[n].real()*F[n].imag());
611  }
612  double gR = (gp-gx)/2.;
613  double gr = (gp+gx)/2.;
614  double gc = sqrt(gR*gR+gI*gI);
615 
616  if (polarization==0) z = sqrt(fabs(gr-gc)<1e-8?0.:gr-gc);
617  if (polarization==1) z = sqrt(gr+gc);
618  if (polarization==2&&sqrt(gr+gc)>0) z=sqrt((gr-gc)/(gr+gc));
619  if (polarization==3) z = sqrt(2*gr/nIFO);
620  if (polarization==4) z = (fabs(gr-gc)<1e-8?0.:gr-gc);
621  if (polarization==5) z = (gr+gc);
622  if (polarization==6) z = F[0].imag();
623  if (polarization==7) z = F[0].real();
624  if (polarization==8) z = (F[1].imag()!=0) ? F[0].imag()/F[1].imag() : 0.;
625  if (polarization==9) z = (F[1].real()!=0) ? F[0].real()/F[1].real() : 0.;
626  if (polarization==10) {
627  double F0=sqrt(pow(F[0].real(),2)+pow(F[0].imag(),2));
628  double F1=sqrt(pow(F[1].real(),2)+pow(F[1].imag(),2));
629  z = (F1!=0) ? F0/F1 : 0.;
630  }
631  if (polarization==12) { // DPF Angle
632  double cos2G = gc!=0 ? gR/gc : 0; // (Fp2-Fc2)/Norm
633  double sin2G = gc!=0 ? gI/gc : 0; // 2*Fpc/Norm
634  double r2d = 180./TMath::Pi();
635  z = r2d*atan2(sin2G,cos2G)/2.;
636  }
637 
638  return z;
639 }
640 
641 //______________________________________________________________________________
642 double
643 gnetwork::GetAntennaPattern(TString ifo, double phi, double theta, double psi, bool plus) {
644 //
645 // Get Detector Antenna Pattern value
646 //
647 // Input: ifo - detector name (Ex: "L1","H1","V1")
648 // phi - phi (degrees)
649 // theta - theta (degrees)
650 // psi - polarization angle (degrees)
651 // plus - true/false -> f+/fx
652 //
653 
654  int nIFO = this->ifoListSize();
655  if(nIFO<1) {cout << "gnetwork::GetAntennaPattern nIFO must be >= 1" << endl;return -100;}
656 
657  detector *pD[NIFO];
658  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
659 
660  int ifoId=0;
661  int nifo=0;
662  for(int n=0; n<nIFO; n++) {
663  if(ifo.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId=n;}
664  }
665  if(nifo<1) {cout << "gnetwork::GetAntennaPattern ifo not valid" << endl;return -100;}
666 
667  wavecomplex F = pD[ifoId]->antenna(theta,phi,psi);
668  if(plus) return F.real(); else return F.imag();
669 
670 }
671 
672 //______________________________________________________________________________
673 void
674 gnetwork::DrawAntennaPattern(int polarization, int dpaletteId, bool btitle, int order) {
675 //
676 // Draw Network Antenna Pattern sky map in the DPF : Dominant Polarization Frame
677 //
678 // Input: polarization - DPF component
679 //
680 // 0 -> |Fx| - DPF
681 // 1 -> |F+| - DPF
682 // 2 -> |Fx|/|F+| - DPF
683 // 3 -> sqrt((|F+|^2+|Fx|^2)/nIFO) - DPF
684 // 4 -> |Fx|^2 - DPF
685 // 5 -> |F+|^2 - DPF
686 // 6 -> Fx - only with 1 detector
687 // 7 -> F+ - only with 1 detector
688 // 8 -> F1x/F2x - only with 2 detectors
689 // 9 -> F1+/F2+ - only with 2 detectors
690 // 10 -> sqrt(|F1+|^2+|F1x|^2)/sqrt(|F2+|^2+|F2x|^2) - only with 2 detectors
691 // 11 -> The same as (10) but averaged over psi - only with 2 detectors
692 // 12 -> DPF Angle
693 //
694 // dpaletteId - palette ID
695 // btitle - if true then a default title is provided
696 // order - healpix order used for Antenna Pattern
697 //
698 
699  if(gSM.goff) return;
700 
701  gSM=skymap(int(order));
702 
703  if(polarization<0||polarization>12) {
704  cout << "---------------------------------------------------------------------------------------" << endl;
705  cout << "polarization must be [0-12]" << endl;
706  cout << "---------------------------------------------------------------------------------------" << endl;
707  cout << "polarization=0 -> |Fx| DPF" << endl;
708  cout << "polarization=1 -> |F+| DPF" << endl;
709  cout << "polarization=2 -> |Fx|/|F+| DPF" << endl;
710  cout << "polarization=3 -> sqrt((|F+|^2+|Fx|^2)/nIFO) DPF" << endl;
711  cout << "polarization=4 -> |Fx|^2 DPF" << endl;
712  cout << "polarization=5 -> |F+|^2 DPF" << endl;
713  cout << "polarization=6 -> Fx only with 1 detector" << endl;
714  cout << "polarization=7 -> F+ only with 1 detector" << endl;
715  cout << "polarization=8 -> F1x/F2x only with 2 detectors" << endl;
716  cout << "polarization=9 -> F1+/F2+ only with 2 detectors" << endl;
717  cout << "polarization=10 -> sqrt(|F1+|^2+|F1x|^2)/sqrt(|F2+|^2+|F2x|^2) only with 2 detectors" << endl;
718  cout << "polarization=11 -> The same as (10) but averaged over psi only with 2 detectors" << endl;
719  cout << "polarization=12 -> DPF Angle" << endl;
720  cout << "---------------------------------------------------------------------------------------" << endl;
721  return;
722  }
723 
724  int nIFO = this->ifoListSize();
725  if(nIFO==0)
726  {cout << "detectors are no declared" << endl;return;}
727  if((nIFO!=1)&&(polarization==6||polarization==7))
728  {cout << "with polarization [6-7] works only with one detector" << endl;return;}
729  if((nIFO!=2)&&(polarization>=8&&polarization<=11))
730  {cout << "with polarization [8-11] works only two one detectors" << endl;return;}
731 
732  detector *pD[NIFO];
733  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
734 
735  TString pTitle="";
736  switch(polarization) {
737  case 0:
738  pTitle=TString("|F_{x}|");
739  break;
740  case 1:
741  pTitle=TString("|F_{+}|");
742  break;
743  case 2:
744  pTitle=TString("|F_{x}|/|F_{+}|");
745  break;
746  case 3:
747  pTitle=TString("#sqrt{(|F_{+}|^{2}+|F_{x}|^{2})/nIFO}");
748  break;
749  case 4:
750  pTitle=TString("|F_{x}|^{2}");
751  break;
752  case 5:
753  pTitle=TString("|F_{+}|^{2}");
754  break;
755  case 6:
756  pTitle=TString("F_{x}");
757  break;
758  case 7:
759  pTitle=TString("F_{+}");
760  break;
761  case 8:
762  pTitle=TString("F1_{x}/F2_{x}");
763  break;
764  case 9:
765  pTitle=TString("F1_{+}/F2_{+}");
766  break;
767  case 10:
768  pTitle=TString("#sqrt{|F1_{+}|^{2}+|F1_{x}|^{2})/sqrt(|F2_{+}|^{2}+|F2_{x}|^{2}}");
769  break;
770  case 12:
771  pTitle=TString("DPF Angle");
772  break;
773  default:
774  break;
775  }
776  TString netName="";
777  for(int n=0; n<nIFO; n++) {
778  TString ifoName=this->getifo(n)->Name;
779  if(ifoName.CompareTo("A1")==0) ifoName=TString("#tilde{A}");
780  else if(ifoName.CompareTo("H2")==0) ifoName=TString("#tilde{H}");
781  else if(ifoName.CompareTo("I2")==0) ifoName=TString("#tilde{I}");
782  else if(ifoName.CompareTo("R2")==0) ifoName=TString("#tilde{R}");
783  else ifoName.Resize(1);
784  netName+=ifoName;
785  }
786  if(btitle) gSM.h2->SetTitle("Network = "+netName+
787  " "+
788  "Antenna Pattern = "+pTitle);
789 
790  std::vector<wavecomplex> F(nIFO);
791  int size=int(180*2*gSM.resolution*360*2*gSM.resolution);
792  double* x = new double[size];
793  double* y = new double[size];
794  double* z = new double[size];
795  int cnt=0;
796  for (int i=0;i<360*2*gSM.resolution;i++) {
797  for (int j=0;j<180*2*gSM.resolution;j++) {
798  double phi = (double)(i+1)/(double)(2*gSM.resolution);
799  double theta = (double)(j+1)/(double)(2*gSM.resolution);
800  double psi=0;
801  x[cnt]=phi;
802  y[cnt]=theta;
803  if (polarization!=11) {
804 /*
805  psi=0;
806  for (int n=0;n<nIFO;n++) F[n] = pD[n]->antenna(theta,phi,psi);
807  double gp=0,gx=0,gI=0;
808  for (int n=0;n<nIFO;n++) {
809  gp+=F[n].real()*F[n].real();
810  gx+=F[n].imag()*F[n].imag();
811  gI+=F[n].real()*F[n].imag();
812  }
813  double gR = (gp-gx)/2.;
814  double gr = (gp+gx)/2.;
815  double gc = sqrt(gR*gR+gI*gI);
816 
817  if (polarization==0) z[cnt] = sqrt(fabs(gr-gc)<1e-8?0.:gr-gc);
818  if (polarization==1) z[cnt] = sqrt(gr+gc);
819  if (polarization==2&&sqrt(gr+gc)>0) z[cnt]=sqrt((gr-gc)/(gr+gc));
820  if (polarization==3) z[cnt] = sqrt(2*gr/nIFO);
821  //if (polarization==3) z[cnt] = sqrt(2*gr);
822  if (polarization==4) z[cnt] = (fabs(gr-gc)<1e-8?0.:gr-gc);
823  if (polarization==5) z[cnt] = (gr+gc);
824  if (polarization==6) z[cnt] = F[0].imag();
825  if (polarization==7) z[cnt] = F[0].real();
826  if (polarization==8) z[cnt] = (F[1].imag()!=0) ? F[0].imag()/F[1].imag() : 0.;
827  if (polarization==9) z[cnt] = (F[1].real()!=0) ? F[0].real()/F[1].real() : 0.;
828  //if (polarization==8) z[cnt] = F[0].imag()/F[1].imag()>0 ? 1. : -1.;
829  //if (polarization==9) z[cnt] = F[0].real()/F[1].real()>0 ? 1. : -1.;
830  if (polarization==10) {
831  double F0=sqrt(pow(F[0].real(),2)+pow(F[0].imag(),2));
832  double F1=sqrt(pow(F[1].real(),2)+pow(F[1].imag(),2));
833  z[cnt] = (F1!=0) ? F0/F1 : 0.;
834  }
835 */
836  z[cnt] = GetAntennaPattern(phi, theta, psi, polarization);
837  } else {
838  int counts=0;
839  z[cnt] = 0.;
840  //for (int k=0;k<360;k+=8) {
841  for (int k=0;k<1;k+=1) {
842  psi=(double)k;
843  for (int n=0;n<nIFO;n++) F[n] = pD[n]->antenna(theta,phi,psi);
844  //double F0=sqrt(pow(F[0].real(),2)+pow(F[0].imag(),2));
845  //double F1=sqrt(pow(F[1].real(),2)+pow(F[1].imag(),2));
846  double F0=F[0].real()+F[0].imag();
847  double F1=F[1].real()+F[1].imag();
848  z[cnt] += (F1!=0) ? fabs(F0/F1) : 0.;
849  }
850  z[cnt] /= (double)counts;
851  }
852  // fill skymap
853  int ind = gSM.getSkyIndex(theta,phi);
854  gSM.set(ind,z[cnt]);
855  cnt++;
856  }
857  }
858 
859  if (gSM.coordinate.CompareTo("GEOGRAPHIC")==0) {
860  for (int i=0;i<size;i++) CwbToGeographic(x[i],y[i],x[i],y[i]);
861  }
862  if (gSM.coordinate.CompareTo("CELESTIAL")==0) {
863  for (int i=0;i<size;i++) CwbToCelestial(x[i],y[i],x[i],y[i]);
864  }
865  gSM.FillData(size, x, y, z);
866  gSM.Draw(dpaletteId);
867 
868  delete [] x;
869  delete [] y;
870  delete [] z;
871 }
872 
873 //______________________________________________________________________________
874 double
875 gnetwork::GetDelay(TString ifo1, TString ifo2, TString mode) {
876 //
877 // Return the delay (sec) between two detector
878 //
879 // Input: ifo1 - detector name of the first detector (Ex: "L1")
880 // ifo2 - detector name of the second detector (Ex: "V1")
881 // mode - delay mde
882 // "" : (max-min)/2 - default
883 // min : minum delay value
884 // max : maximum delay value
885 // MAX : maximum absolute delay value
886 //
887 
888  int nIFO = this->ifoListSize();
889  if(nIFO<2) {cout << "nIFO must be >= 2" << endl;return -1.;}
890 
891  if(ifo1.CompareTo(ifo2)==0) return 0.0;
892 
893  detector *pD[NIFO];
894  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
895 
896  int ifoId1=0;
897  int ifoId2=0;
898  int nifo=0;
899  for(int n=0; n<nIFO; n++) {
900  if(ifo1.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId1=n;}
901  if(ifo2.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId2=n;}
902  }
903  if(nifo<2) {cout << "gnetwork::GetDelay ifo not valid" << endl;return -1.;}
904 
905  skymap sm0, sTau;
906  double maxTau = -1.;
907  double minTau = 1.;
908  double tmax, tmin;
909 
910  tmax = tmin = 0.;
911  for(int n=0; n<nifo; n++) {
912  sTau = pD[n]->tau;
913  tmax = sTau.max();
914  tmin = sTau.min();
915  if(tmax > maxTau) maxTau = tmax;
916  if(tmin < minTau) minTau = tmin;
917  }
918  if(strstr(mode,"min")) return minTau;
919  if(strstr(mode,"max")) return maxTau;
920  if(strstr(mode,"MAX")) return fabs(maxTau)>fabs(minTau) ? fabs(maxTau) : fabs(minTau);
921  return (maxTau-minTau)/2.;
922 }
923 
924 //______________________________________________________________________________
925 double
926 gnetwork::GetDelay(TString ifo1, TString ifo2, double phi, double theta) {
927 //
928 // Return the delay (sec) between two detector from the direction phi,theta
929 //
930 // Input: ifo1 - detector name of the first detector (Ex: "L1")
931 // ifo2 - detector name of the second detector (Ex: "V1")
932 // phi - phi direction (degrees)
933 // theta - theta direction (degrees)
934 //
935 
936  int nIFO = this->ifoListSize();
937  if(nIFO<2) {cout << "gnetwork::GetDelay nIFO must be >= 2" << endl;return -1.;}
938 
939  if(ifo1.CompareTo(ifo2)==0) return 0.0;
940 
941  detector *pD[NIFO];
942  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
943 
944  int ifoId1=0;
945  int ifoId2=0;
946  int nifo=0;
947  for(int n=0; n<nIFO; n++) {
948  if(ifo1.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId1=n;}
949  if(ifo2.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId2=n;}
950  }
951  if(nifo<2) {cout << "gnetwork::GetDelay ifo not valid" << endl;return -1.;}
952 
953  return pD[ifoId1]->getTau(theta,phi)-pD[ifoId2]->getTau(theta,phi);
954 }
955 
956 //______________________________________________________________________________
957 void
958 gnetwork::DrawDelay(TString ifo1, TString ifo2, double phi, double theta, double frequency) {
959 
960  if(gSM.goff) return;
961 
962  int nIFO = this->ifoListSize();
963  if(nIFO<2) {cout << "gnetwork::DrawDelay nIFO must be >= 2" << endl;return;}
964 
965  if(ifo1.CompareTo(ifo2)==0) {cout << "gnetwork::DrawDelay ifo1 must != ifo2" << endl;return;}
966 
967  double ifo1d2_pos_delay=0.;
968  // if ifo1d2_pos_delay!=0 draw delay respect to ifo1d2_pos_delay
969  if(phi!=1000) ifo1d2_pos_delay=GetDelay(ifo1, ifo2, phi, theta);
970 
971  //double ifo1d2_max_delay=0.;
972  //if(frequency>0) ifo1d2_max_delay=GetDelay(ifo1, ifo2, "max");
973 
974  detector *pD[NIFO];
975  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
976 
977  int ifoId1=0;
978  int ifoId2=0;
979  int nifo=0;
980  for(int n=0; n<nIFO; n++) {
981  if(ifo1.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId1=n;}
982  if(ifo2.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId2=n;}
983  }
984  if(nifo<2) {cout << "gnetwork::DrawDelay ifo not valid" << endl;return;}
985 
986  int size=int(180*2*gSM.resolution*360*2*gSM.resolution);
987  double* x = new double[size];
988  double* y = new double[size];
989  double* z = new double[size];
990  int cnt=0;
991  for (int i=0;i<360*2*gSM.resolution;i++) {
992  for (int j=0;j<180*2*gSM.resolution;j++) {
993  double phi = (double)(i+1)/(double)(2*gSM.resolution);
994  double theta = (double)(j+1)/(double)(2*gSM.resolution);
995 
996  double delay=pD[ifoId1]->getTau(theta,phi)-pD[ifoId2]->getTau(theta,phi);
997 
998  delay-=ifo1d2_pos_delay;
999  if(frequency>0) {
1000  double fdelay=fabs(fmod(delay,0.5/frequency));
1001  int idelay = int((fabs(delay)-fdelay)/(0.5/frequency));
1002  delay = idelay%2==0 ? fdelay : (0.5/frequency)-fdelay;
1003  }
1004 
1005  x[cnt]=phi;
1006  y[cnt]=theta;
1007  z[cnt]=delay;
1008  cnt++;
1009  }
1010  }
1011 
1012  if (gSM.coordinate.CompareTo("GEOGRAPHIC")==0) {
1013  for (int i=0;i<size;i++) CwbToGeographic(x[i],y[i],x[i],y[i]);
1014  }
1015  if (gSM.coordinate.CompareTo("CELESTIAL")==0) {
1016  for (int i=0;i<size;i++) CwbToCelestial(x[i],y[i],x[i],y[i]);
1017  }
1018  gSM.FillData(size, x, y, z);
1019  gSM.Draw(0);
1020 
1021  delete [] x;
1022  delete [] y;
1023  delete [] z;
1024 }
1025 
1026 //______________________________________________________________________________
1027 int
1028 gnetwork::Delay2Coordinates(double t1, double t2, double t3, double& ph1, double& th1, double& ph2,double& th2) {
1029 //
1030 // Return the direction & mirror-direction uding the input arrival times of the 3 detectors
1031 // NOTE: Works only with 3 detectors
1032 // With 3 detector there are 2 possible solutions
1033 //
1034 // Input: t1,t2,t3 - arrival times of the 3 detectors (sec)
1035 //
1036 // Output: ph1,th1 - phi,theta coordinates of the first direction
1037 // ph2,th2 - phi,theta coordinates of the second direction
1038 //
1039 // Return 0/1 - success/error
1040 //
1041 
1042  int nIFO = this->ifoListSize();
1043  if(nIFO!=3) {cout << "gnetwork::Delay2Source nIFO must be = 3" << endl;return 1;}
1044 
1045  TString ifoName[3];
1046  for(int n=0; n<nIFO; n++) ifoName[n]=TString(this->getifo(n)->Name);
1047  //for(int n=0; n<nIFO; n++) cout << n << " " << ifoName[n].Data() << " " << ifoId[n] << endl;
1048  for(int n=0; n<nIFO; n++) for(int m=n+1; m<nIFO; m++)
1049  if(ifoName[n].CompareTo(ifoName[m])==0)
1050  {cout << "gnetwork::Delay2Source ifos must differents" << endl;return 2;}
1051 
1052  detector *pD[NIFO];
1053  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
1054 
1055  if(!gROOT->GetClass("XYZVector")) gSystem->Load("libMathCore");
1056 
1057  XYZVector vD[3];
1058  for(int n=0; n<nIFO; n++) vD[n].SetXYZ(pD[n]->Rv[0],pD[n]->Rv[1],pD[n]->Rv[2]);
1059  for(int n=0; n<nIFO; n++) vD[n]/=speedlight;
1060 
1061  double a=sqrt((vD[1]-vD[0]).Mag2());
1062  double b=sqrt((vD[2]-vD[1]).Mag2());
1063  double c=sqrt((vD[2]-vD[0]).Mag2());
1064 
1065  double s=(a+b+c)/2;
1066  double b2=(2/a)*sqrt(s*(s-a)*(s-b)*(s-c));
1067  double b1=sqrt(c*c-b2*b2);
1068 
1069  double delta=pow(b1*(t1-t2)-a*(t1-t3),2)+pow(b2*(t1-t2),2);
1070 
1071  double stheta = sqrt(delta)/(a*b2);
1072  double sphi = -b2*(t1-t2)/sqrt(delta);
1073 
1074  if (stheta>1) stheta=1;
1075  if (stheta<-1) stheta=-1;
1076 
1077  double theta = asin(stheta);
1078  double phi = asin(sphi);
1079 
1080  double ctheta = sqrt(pow(a*b2,2)-delta)/(a*b2);
1081  double cphi = (a*(t1-t3)-b1*(t1-t2))/sqrt(delta);
1082 
1083  if (TMath::IsNaN(ctheta)) ctheta=-1;
1084  if (ctheta>1) ctheta=-1;
1085  if (ctheta<-1) ctheta=-1;
1086 
1087  double thetacn = acos(-ctheta);
1088 
1089  XYZVector D1 = vD[1]-vD[0];
1090  XYZVector D3 = D1.Cross(vD[2]-vD[0]);
1091  XYZVector D2 = D1.Cross(D3);
1092 
1093  XYZVector eD1 = D1.Unit();
1094  XYZVector eD2 = D2.Unit();
1095  XYZVector eD3 = D3.Unit();
1096 
1097  Rotation3D R(eD1,eD2,eD3);
1098 
1099  //XYZVector IS;
1100  //XYZVector S;
1101  //S.SetR(1);
1102  TVector3 IS;
1103  TVector3 S(1,0,0);
1104 
1105  if (cphi<0) {
1106 
1107  S.SetTheta(theta);
1108  S.SetPhi((TMath::PiOver2()+phi));
1109  IS = R*S;
1110  th1=IS.Theta()*180./TMath::Pi();
1111  ph1=IS.Phi()*180./TMath::Pi();
1112  if (ph1<0) ph1+=360;
1113  if (ph1>360) ph1-=360;
1114 
1115  S.SetTheta(thetacn);
1116  S.SetPhi((TMath::PiOver2()+phi));
1117  IS = R*S;
1118  th2=IS.Theta()*180./TMath::Pi();
1119  ph2=IS.Phi()*180./TMath::Pi();
1120  if (ph2<0) ph2+=360;
1121  if (ph2>360) ph2-=360;
1122 
1123  } else {
1124 
1125  S.SetTheta(theta);
1126  S.SetPhi(-(TMath::PiOver2()+phi));
1127  IS = R*S;
1128  th1=IS.Theta()*180./TMath::Pi();
1129  ph1=IS.Phi()*180./TMath::Pi();
1130  if (ph1<0) ph1+=360;
1131  if (ph1>360) ph1-=360;
1132 
1133  S.SetTheta(thetacn);
1134  S.SetPhi(-(TMath::PiOver2()+phi));
1135  IS = R*S;
1136  th2=IS.Theta()*180./TMath::Pi();
1137  ph2=IS.Phi()*180./TMath::Pi();
1138  if (ph2<0) ph2+=360;
1139  if (ph2>360) ph2-=360;
1140  }
1141 
1142  return 0;
1143 }
1144 
1145 //______________________________________________________________________________
1146 void
1147 gnetwork::MakeNetworkDetectorIndex(double gamma, int loop, double snr) {
1148 //
1149 // Make Network Detector Index (1G analysis) and store it to the internal gSM sky map
1150 //
1151 // Input: gamma - 1G gamma value
1152 // loop - increase the statistic
1153 // snr - add noise to signal -> snr
1154 //
1155 
1156  if(gamma>1) {cout << "gamma must be <=1" << endl;exit(1);}
1157  if(loop<1) {cout << "gamma must be >=1" << endl;exit(1);}
1158 
1159  int nIFO = this->ifoListSize();
1160  if(nIFO==0) {cout << "detectors are no declared" << endl;exit(1);}
1161 
1162  detector *pD[NIFO];
1163  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
1164 
1165  std::vector<wavecomplex> F(nIFO);
1166  double* X = new double[nIFO];
1167  double* u = new double[nIFO];
1168  double* v = new double[nIFO];
1169 
1170  int L = gSM.size();
1171 
1172  double h2max=0.;
1173  gRandom->SetSeed(0);
1174  for (int l=0;l<L;l++) {
1175  double ndi=0.;
1176  for (int i=0;i<loop;i++) {
1177 
1178  double phi=gSM.getPhi(l);
1179  double theta=gSM.getTheta(l);
1180  double psi=gRandom->Uniform(0,180);
1181 
1182  double hp=gRandom->Uniform(-1,1);
1183  double hx=gRandom->Uniform(-1,1);
1184 
1185  double E=0.;
1186  for (int n=0;n<nIFO;n++) {
1187  F[n] = pD[n]->antenna(theta,phi,psi);
1188  X[n] = hp*F[n].real()+hx*F[n].imag();
1189  if(snr>=0.0) E+=X[n]*X[n];
1190  }
1191  // add noise to signal -> snr
1192  E=sqrt(E);
1193  if(snr>=0.0) for (int n=0;n<nIFO;n++) X[n] = X[n]*snr/E + gRandom->Gaus(0.,1.);
1194 
1195  double gp=0,gx=0,gI=0;
1196  for (int n=0;n<nIFO;n++) {
1197  gp+=F[n].real()*F[n].real();
1198  gx+=F[n].imag()*F[n].imag();
1199  gI+=F[n].real()*F[n].imag();
1200  }
1201  gp+=1.e-12;
1202  gx+=1.e-12;
1203 
1204  double Xp=0;
1205  double Xx=0;
1206  for (int k=0;k<nIFO;k++) {
1207  Xp+=X[k]*F[k].real();
1208  Xx+=X[k]*F[k].imag();
1209  }
1210 
1211  double uc = (Xp*gx - Xx*gI); // u cos of rotation to PCF
1212  double us = (Xx*gp - Xp*gI); // u sin of rotation to PCF
1213  double vc = (gp*uc + gI*us); // (u*f+)/|u|^2 - 'cos' for v
1214  double vs = (gx*us + gI*uc); // (u*fx)/|u|^2 - 'sin' for v
1215  double um=0;
1216  double vm=0;
1217  for (int k=0;k<nIFO;k++) {
1218  u[k]=0;
1219  v[k]=0;
1220  u[k]=F[k].real()*uc+F[k].imag()*us;
1221  um+=u[k]*u[k];
1222  v[k]=-F[k].imag()*vc+F[k].real()*vs;
1223  vm+=v[k]*v[k];
1224  }
1225  vm+=1.e-24;
1226  if(gamma>0) {
1227  double GAMMA = 1.-gamma*gamma; // network regulator
1228  ndi+=this->netx(u,um,v,vm,GAMMA);
1229  } else ndi+=this->ndi(u,um,nIFO);
1230 /*
1231  } else {
1232  double xndi=this->ndi(u,um,nIFO);
1233  if(xndi>=-gamma) ndi+=1;
1234  }
1235 */
1236  }
1237  double xndi=ndi/loop;
1238  if(h2max<xndi) h2max=xndi;
1239  gSM.set(l,xndi);
1240  }
1241 
1242  //if(gamma<0) h2->GetZaxis()->SetRangeUser(1,nIFO);
1243  if(gamma<0) gSM.h2->GetZaxis()->SetRangeUser(1,h2max);
1244 
1245  delete [] X;
1246  delete [] u;
1247  delete [] v;
1248 }
1249 
1250 //______________________________________________________________________________
1251 void
1252 gnetwork::DrawNetworkDetectorIndex(double gamma, int loop, double snr, int dpaletteId, bool btitle) {
1253 //
1254 // Draw Network Detector Index sky map (1G analysis)
1255 //
1256 // Input: gamma - 1G gamma value
1257 // loop - increase the statistic
1258 // snr - add noise to signal -> snr
1259 //
1260 // dpaletteId - palette ID
1261 // btitle - if true then a default title is provided
1262 //
1263 
1264  if(gSM.size()==0) {
1265  cout << "gnetwork::DrawNetworkDetectorIndex : Error - network skymap not initialized" << endl;
1266  exit(1);
1267  }
1268 
1269  MakeNetworkDetectorIndex(gamma, loop, snr);
1270 
1271  int nIFO = this->ifoListSize();
1272 
1273  TString netName="";
1274  for(int n=0; n<nIFO; n++) {
1275  TString ifoName=this->getifo(n)->Name;
1276  if(ifoName.CompareTo("A1")==0) ifoName=TString("#tilde{A}");
1277  else if(ifoName.CompareTo("H2")==0) ifoName=TString("#tilde{H}");
1278  else if(ifoName.CompareTo("I2")==0) ifoName=TString("#tilde{I}");
1279  else if(ifoName.CompareTo("R2")==0) ifoName=TString("#tilde{R}");
1280  else ifoName.Resize(1);
1281  netName+=ifoName;
1282  }
1283  if(btitle) gSM.h2->SetTitle("Network = "+netName+
1284  " "+
1285  "Network Detector Index");
1286 
1287  gSM.FillData();
1288  gSM.Draw(dpaletteId);
1289 }
1290 
1291 //______________________________________________________________________________
1292 void
1293 gnetwork::DrawCircles(double phi, double theta, Color_t lcolor, Width_t lwidth, Style_t lstyle, bool labels) {
1294 //
1295 // Draw iso delay circles for each couple of detectors along the direction phi,theta
1296 //
1297 // Input: phi,theta - sky direction coordinates (degrees)
1298 // lcolor - line color
1299 // lwidth - line width
1300 // lstyle - line style
1301 //
1302 
1303  DrawCircles(phi, theta, 0., lcolor, lwidth, lstyle, labels);
1304 }
1305 
1306 //______________________________________________________________________________
1307 void
1308 gnetwork::DrawCircles(double phi, double theta, double gps, Color_t lcolor, Width_t lwidth, Style_t lstyle, bool labels) {
1309 //
1310 // Draw iso delay circles for each couple of detectors along the direction phi,theta at time gps
1311 //
1312 // Input: phi,theta - sky direction coordinates (degrees)
1313 // gps - gps time (sec)
1314 // lcolor - line color
1315 // lwidth - line width
1316 // lstyle - line style
1317 //
1318 
1319  if(gSM.goff) return;
1320  if(gSM.canvas==NULL) return;
1321 
1322  int n,m,l;
1323 
1324  int nIFO = this->ifoListSize();
1325  if(nIFO<2) return;
1326 
1327  detector *pD[NIFO];
1328  for(n=0; n<nIFO; n++) pD[n] = this->getifo(n);
1329 
1330  gSM.canvas->cd();
1331 
1332  double phi_off = gps>0 ? gSM.phi2RA(phi, gps) : phi;
1333  phi_off=fmod(phi_off-phi+360,360);
1334 
1335  double ph1,th1,ph2,th2;
1336  double r2d = 180./TMath::Pi();
1337  double d2r = TMath::Pi()/180.;
1338 //Color_t colors[6] = {kBlue,kRed,kYellow,kGreen,kBlack,kWhite};
1339  Style_t marker[3] = {20,21,22}; // circle,square,triangle
1340  int ncircle=0;
1341  for(n=0;n<nIFO;n++) {
1342  for(m=n+1;m<nIFO;m++) {
1343 
1344  XYZVector D1(pD[n]->Rv[0],pD[n]->Rv[1],pD[n]->Rv[2]);
1345  XYZVector D2(pD[m]->Rv[0],pD[m]->Rv[1],pD[m]->Rv[2]);
1346 
1347  D1=D1/speedlight;
1348  D2=D2/speedlight;
1349 
1350  XYZVector D12 = D1-D2;
1351  if(D12.Mag2()==0) continue; // the 2 detectors are located in the same site
1352  TVector3 vD12(D12.X(),D12.Y(),D12.Z());
1353 
1354  TVector3 vSD(1,0,0);
1355  // coordinates are transformed in the CWB frame
1356  if (gSM.coordinate.CompareTo("GEOGRAPHIC")==0) {
1357  GeographicToCwb(phi,theta,ph1,th1);
1358  vSD.SetTheta(th1*d2r);
1359  vSD.SetPhi(ph1*d2r);
1360  } else if (gSM.coordinate.CompareTo("CELESTIAL")==0) {
1361  CelestialToCwb(phi,theta,ph1,th1);
1362  ph1=fmod(ph1-phi_off+360,360);
1363  vSD.SetTheta(th1*d2r);
1364  vSD.SetPhi(ph1*d2r);
1365  } else {
1366  vSD.SetTheta(theta*d2r);
1367  vSD.SetPhi(phi*d2r);
1368  }
1369 
1370  TVector3 rvSD;
1371  double as;
1372  int L=4*360;
1373  wavearray<float> th(L);
1374  wavearray<float> ph(L);
1375  for (l=0;l<L;l++) {
1376  as = l*d2r/4.;
1377  TRotation vR;
1378  vR.Rotate(as,vD12);
1379  rvSD = vR*vSD;
1380  th[l] = rvSD.Theta()*r2d;
1381  ph[l] = rvSD.Phi()*r2d+phi_off;
1382  ph[l]=fmod(ph[l]+360,360);
1383  }
1384 
1385  wavearray<double> x(L);
1386  wavearray<double> y(L);
1387  int P=0;
1388  for (l=0;l<L;l++) {
1389  // the main circle is shown in sub-circles to avoid to draw fake lines from ph=0 to 360
1390  // if((fabs(ph[l]-ph[l==0?L-1:l-1])<180)&&(l<L-1)) {
1391  // TVector3 is in radians [phi=-Pi,+Pi] [theta=0,Pi]
1392  ph1=ph[l];
1393  ph2=ph[l==0?L-1:l-1];
1394  if (gSM.coordinate=="GEOGRAPHIC") {
1395  CwbToGeographic(ph1,th1,ph1,th1);
1396  CwbToGeographic(ph2,th2,ph2,th2);
1397  }
1398  if((fabs(ph1-ph2)<180)&&(l<L-1)) {
1399  x[P]=ph[l]; y[P]=th[l];
1400 
1401  double X=x[P];
1402  if (gSM.projection=="hammer") {
1403  if(X>180) {x[P]+=360; y[P]-=90;}
1404  else {y[P]-=90;}
1405  gSM.ProjectHammer(x[P], y[P], x[P], y[P]);
1406  y[P]=y[P]+90;
1407  } else if (gSM.projection=="sinusoidal") {
1408  if(X>180) {x[P]-=360; y[P]+=90;}
1409  else {y[P]-=90;}
1410  gSM.ProjectSinusoidal(x[P], y[P], x[P], y[P]);
1411  if(X>180) {x[P]=-x[P]; y[P]-=90;}
1412  else {y[P]+=90;}
1413  } else if (gSM.projection=="parabolic") {
1414  if(X>180) {x[P]-=360; y[P]-=90;}
1415  else {y[P]-=90;}
1416  gSM.ProjectParabolic(x[P], y[P], x[P], y[P]);
1417  y[P]+=90;
1418  } else {
1419 // x[P]=x[P]-180; y[P]=y[P]-90;
1420  }
1421  if (gSM.coordinate=="GEOGRAPHIC") CwbToGeographic(x[P],y[P],x[P],y[P]);
1422  if (gSM.coordinate=="CELESTIAL") CwbToCelestial(x[P],y[P],x[P],y[P]);
1423  P++;
1424  } else {
1425  TPolyLine *pL = new TPolyLine(P,x.data,y.data);
1426  //pL->SetLineColor(colors[ncircle%6]);
1427  pL->SetLineColor(lcolor);
1428  pL->SetLineStyle(lstyle);
1429  pL->SetLineWidth(lwidth);
1430  pL->Draw();
1431  circleL.push_back(pL);
1432  P=0;
1433  }
1434  }
1435  // draw detector axis
1436  double vph,vth;
1437  if(vSD.Dot(vD12)>0) {
1438  vph=D12.Phi()*r2d;
1439  vth=D12.Theta()*r2d;
1440  } else {
1441  vph=(-D12).Phi()*r2d;
1442  vth=(-D12).Theta()*r2d;
1443  }
1444  vph=fmod(vph+phi_off+360,360);
1445 
1446  double X=vph;
1447  if (gSM.projection=="hammer") {
1448  if(X>180) {vph+=360; vth-=90;}
1449  if(X<180) {vth-=90;}
1450  gSM.ProjectHammer(vph, vth, vph, vth);
1451  vth=vth+90;
1452  } else if (gSM.projection=="sinusoidal") {
1453  if(X>180) {vph-=360; vth+=90;}
1454  else {vth-=90;}
1455  gSM.ProjectSinusoidal(vph, vth, vph, vth);
1456  if(X>180) {vph=-vph; vth-=90;}
1457  else {vth+=90;}
1458  } else if (gSM.projection=="parabolic") {
1459  if(X>180) {vph-=360; vth-=90;}
1460  else {vth-=90;}
1461  gSM.ProjectParabolic(vph, vth, vph, vth);
1462  vth+=90;
1463  } else {
1464 // vph=vph-180; vth=vth-90;
1465  }
1466  if (gSM.coordinate=="GEOGRAPHIC") CwbToGeographic(vph,vth,vph,vth);
1467  if (gSM.coordinate=="CELESTIAL") CwbToCelestial(vph,vth,vph,vth);
1468  //if (coordinate.CompareTo("CWB")==0) {vph+=180;vth+=90;}
1469  TMarker* pM = new TMarker(vph,vth,marker[ncircle%3]);
1470  pM->SetMarkerSize(1.0);
1471  pM->SetMarkerColor(kBlack);
1472  if(labels) pM->Draw();
1473  daxisM.push_back(pM);
1474 
1475  // draw ifos pair labels on the detector axis
1476  char ifoPair[16];sprintf(ifoPair,"%c%c",pD[n]->Name[0],pD[m]->Name[0]);
1477  TLatex *pL = new TLatex(vph+4.0,vth,ifoPair);
1478  pL->SetTextFont(32);
1479  pL->SetTextSize(0.045);
1480  pL->SetTextColor(kBlack);
1481  if(labels) pL->Draw();
1482  siteP.push_back(pL);
1483 
1484  ncircle++;
1485  }
1486  }
1487 
1488  return;
1489 }
1490 
1491 //______________________________________________________________________________
1492 void gnetwork::DumpObject(char* file)
1493 {
1494 //
1495 // dump gnetwork object into root file
1496 //
1497 // Input: file - root file name
1498 //
1499 
1500  TFile* rfile = new TFile(file, "RECREATE");
1501  this->Write("gnetwork"); // write this object
1502  rfile->Close();
1503 }
1504 
1505 //______________________________________________________________________________
1506 void gnetwork::LoadObject(char* file)
1507 {
1508 //
1509 // load gnetwork object from root file
1510 //
1511 // Input: file - root file name
1512 //
1513 
1514  TFile* rfile = new TFile(file);
1515  gnetwork* gnet = (gnetwork*)rfile->Get("gnetwork;1");
1516  if(gnet==NULL) {
1517  cout << "gnetwork::LoadObject : Error - input file don't contains object gnetwork" << endl;
1518  exit(1);
1519  }
1520  *this=*gnet;
1521 
1522  siteM.clear();
1523  siteA.clear();
1524  siteT.clear();
1525  siteL.clear();
1526  siteP.clear();
1527  circleL.clear();
1528  daxisM.clear();
1529 
1530  gSM.FillData();
1531  rfile->Close();
1532 }
1533 
1534 //______________________________________________________________________________
1535 TCanvas* gnetwork::DrawPolargram(int ptype, network* net)
1536 {
1537 //
1538 // Plot polargrams of the event pixels in the network plane
1539 //
1540 // Input: ptype - plot type [0/1]
1541 // 0 : projection on network plane of 00/90 amplitudes
1542 // 1 : standard response of 00/90 amplitudes
1543 // net - pointer to network object
1544 // if net!=NULL then net is used to retrive the data pixels to be plotted
1545 //
1546 // Output: return canvas object
1547 //
1548 
1549  if(ptype!=0 && ptype!=1) {
1550  cout << "gnetwork::DrawPolargram : Error - wrong input parameter, must be [0/1]" << endl;
1551  exit(1);
1552  }
1553 
1554  network* NET = net!=NULL ? net : this;
1555 
1556  // get size
1557  int size = ptype==0 ? NET->p00_POL[0].size() : NET->r00_POL[0].size();
1558  if(size==0) {
1559  cout << "gnetwork::DrawPolargram : Warning - event pixels size=0" << endl;
1560  return NULL;
1561  }
1562 
1563  // find max radius
1564  double rmax=0;
1565  if(ptype==0) {
1566  for(int i=0;i<size;i++) if(NET->p00_POL[0][i]>rmax) rmax=NET->p00_POL[0][i];
1567  for(int i=0;i<size;i++) if(NET->p90_POL[0][i]>rmax) rmax=NET->p90_POL[0][i];
1568  }
1569  if(ptype==1) {
1570  for(int i=0;i<size;i++) if(NET->r00_POL[0][i]>rmax) rmax=NET->r00_POL[0][i];
1571  for(int i=0;i<size;i++) if(NET->r90_POL[0][i]>rmax) rmax=NET->r90_POL[0][i];
1572  }
1573  rmax=TMath::Nint(1.1*rmax+0.5);
1574 
1575  for(int n=0;n<2;n++) if(grP[n]!=NULL) delete grP[n];
1576  if(canvasP!=NULL) delete canvasP;
1577  canvasP = new TCanvas("canvasP","polargram",800,800);
1578  gStyle->SetLineColor(kBlack);
1579  double* rpol;
1580  double* apol;
1581  for(int n=0;n<2;n++) {
1582  // projection on network plane 00 amplitudes
1583  if(ptype==0 && n==0) {rpol = NET->p00_POL[0].data; apol = NET->p00_POL[1].data;}
1584  // projection on network plane 90 amplitudes
1585  if(ptype==0 && n==1) {rpol = NET->p90_POL[0].data; apol = NET->p90_POL[1].data;}
1586  // standard response 00 amplitudes
1587  if(ptype==1 && n==0) {rpol = NET->r00_POL[0].data; apol = NET->r00_POL[1].data;}
1588  // standard response 90 amplitudes
1589  if(ptype==1 && n==1) {rpol = NET->r90_POL[0].data; apol = NET->r90_POL[1].data;}
1590  grP[n] = new TGraphPolar(size,apol,rpol);
1591 
1592  grP[n]->SetMarkerStyle(20);
1593  grP[n]->SetMarkerSize(0.8);
1594  grP[n]->SetTitle(" ");
1595 
1596  if(n==0) {
1597  grP[n]->SetMarkerColor(kBlue);
1598  grP[n]->SetLineColor(kBlue);
1599  grP[n]->Draw("P");
1600  }
1601 
1602  if(n==1) {
1603  grP[n]->SetMarkerColor(kRed);
1604  grP[n]->SetLineColor(kRed);
1605  grP[n]->Draw("PSAME");
1606  }
1607 
1608  canvasP->Update();
1609  grP[n]->GetPolargram()->SetTextColor(8);
1610  grP[n]->GetPolargram()->SetRangePolar(-TMath::Pi(),TMath::Pi());
1611  grP[n]->SetMinRadial(0);
1612  grP[n]->SetMaxRadial(rmax);
1613  grP[n]->GetPolargram()->SetNdivPolar(612);
1614  grP[n]->GetPolargram()->SetNdivRadial(504);
1615  grP[n]->GetPolargram()->SetToDegree();
1616  grP[n]->GetPolargram()->SetPolarLabelSize(0.03);
1617  grP[n]->GetPolargram()->SetRadialLabelSize(0.03);
1618  }
1619 
1620  canvasP->Update();
1621 
1622  return canvasP;
1623 }
1624 
1625 //______________________________________________________________________________
1626 void
1627 gnetwork::SetScale(TString ifo, double scale) {
1628 //
1629 // set detector scale factor (PSD @ user given frequency)
1630 // return detector scale factor
1631 //
1632 // Input: ifo - detector name
1633 //
1634 // if ifo invalid return -1
1635 //
1636 
1637  int nIFO = this->ifoListSize();
1638  if(nIFO<1) {cout << "gnetwork::GetScale nIFO must be >= 1" << endl;exit(-1);}
1639 
1640  detector *pD[NIFO];
1641  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
1642 
1643  int ifoId=0;
1644  int nifo=0;
1645  for(int n=0; n<nIFO; n++) {
1646  if(ifo.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId=n;}
1647  }
1648 
1649  this->scale[ifoId]=scale;
1650 }
1651 
1652 //______________________________________________________________________________
1653 double
1654 gnetwork::GetScale(TString ifo) {
1655 //
1656 // return detector scale factor (PSD @ user given frequency)
1657 //
1658 // Input: ifo - detector name
1659 //
1660 // if ifo invalid return -1
1661 //
1662 
1663  int nIFO = this->ifoListSize();
1664  if(nIFO<1) {cout << "gnetwork::GetScale nIFO must be >= 1" << endl;return -1.;}
1665 
1666  detector *pD[NIFO];
1667  for(int n=0; n<nIFO; n++) pD[n] = this->getifo(n);
1668 
1669  int ifoId=0;
1670  int nifo=0;
1671  for(int n=0; n<nIFO; n++) {
1672  if(ifo.CompareTo(this->getifo(n)->Name)==0) {nifo++;ifoId=n;}
1673  }
1674 
1675  return scale[ifoId];
1676 }
1677 
network::ifoName
std::vector< char * > ifoName
Definition: network.hh:609
F1
#define F1
Definition: Regression_H1.C:13
gnetwork::DrawPolargram
TCanvas * DrawPolargram(int ptype, network *net=NULL)
Definition: gnetwork.cc:1535
network::veto
wavearray< short > veto
Definition: network.hh:644
network::getifo
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
gnetwork::circleL
std::vector< TPolyLine * > circleL
Definition: gnetwork.hh:122
detector::getDetectorParams
detectorParams getDetectorParams()
Definition: detector.cc:218
gnetwork::GetScale
double GetScale(TString ifo)
Definition: gnetwork.cc:1654
NIFO
#define NIFO
Definition: wat.hh:74
network::skyHole
wavearray< double > skyHole
Definition: network.hh:643
t1
TH1 * t1
Definition: cwb_report_cbc.C:1908
Init
void Init()
Definition: ChirpMass.C:284
network::add
size_t add(detector *)
param: detector structure return number of detectors in the network
Definition: network.cc:2559
gnetwork::DrawSitesArms
void DrawSitesArms(double mlength=600000., Color_t lcolor=kBlack, Size_t lwidth=1.0, Style_t lstyle=1)
Definition: gnetwork.cc:314
network::p90_POL
wavearray< double > p90_POL[2]
buffer for projection on network plane 00 ampl
Definition: network.hh:660
gnetwork::DrawAntennaPattern
void DrawAntennaPattern(int polarization=-1, int dpaletteId=0, bool btitle=true, int order=6)
Definition: gnetwork.cc:674
gnetwork::DrawDelay
void DrawDelay(TString ifo1, TString ifo2, double phi=1000., double theta=1000., double frequency=0.)
Definition: gnetwork.cc:958
gnetwork::canvasP
TCanvas * canvasP
Definition: gnetwork.hh:127
CwbToCelestial
void CwbToCelestial(double ilongitude, double ilatitude, double &olongitude, double &olatitude, double gps=0)
Definition: skycoord.hh:428
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:1308
ptype
TString ptype[nPLOT]
Definition: cwb_mkfad.C:112
a
wavearray< double > a(hp.size())
v
WSeries< float > v[nIFO]
Definition: cwb_net.C:80
n
int n
Definition: cwb_net.C:28
gnetwork::GetSite
double GetSite(TString ifo, TString type="")
Definition: gnetwork.cc:186
network::nSkyStat
skymap nSkyStat
Definition: network.hh:628
z
wavearray< double > z
Definition: Test10.C:32
wavecomplex::imag
double imag() const
Definition: wavecomplex.hh:70
gnetwork::pD
detector * pD[NIFO_MAX]
Definition: gnetwork.hh:125
TString
TString("c")
Definition: cwb_report_skymap.C:129
gnetwork::~gnetwork
virtual ~gnetwork()
Definition: gnetwork.cc:94
gskymap::set
void set(size_t i, double a)
Definition: gskymap.hh:128
gnetwork::Delay2Coordinates
int Delay2Coordinates(double t1, double t2, double t3, double &ph1, double &th1, double &ph2, double &th2)
Definition: gnetwork.cc:1028
network::nPenalty
skymap nPenalty
Definition: network.hh:624
x
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
Definition: BaseFunctionWDM.C:51
frequency
double frequency
Definition: DrawClusterWithSkyplot.C:20
theta
float theta
Definition: cwb_condor_create_ced.C:69
gskymap::coordinate
TString coordinate
Definition: gskymap.hh:229
ph
TH2F * ph
Definition: ReadHistFromJobFile.C:17
watconstants::EarthEquatorialRadius
double EarthEquatorialRadius()
Definition: constants.hh:186
skymap::min
double min()
Definition: skymap.cc:459
gnetwork::DrawSitesLabel
void DrawSitesLabel(Color_t tcolor=kBlack, Size_t tsize=0.045, Font_t tfont=32)
Definition: gnetwork.cc:442
network::p00_POL
wavearray< double > p00_POL[2]
buffers for cluster MRA energy
Definition: network.hh:659
std
STL namespace.
skymap::getTheta
double getTheta(size_t i)
Definition: skymap.hh:224
network::nPolarisation
skymap nPolarisation
Definition: network.hh:630
vR
std::vector< double > vR
Definition: cwb_report_cbc.C:1455
gnetwork::DumpObject
void DumpObject(char *file)
Definition: gnetwork.cc:1492
gskymap::Draw
void Draw(int dpaletteId=1, Option_t *option="colfz")
Definition: gskymap.cc:460
polarization
int polarization
Definition: cwb_draw_antpat.C:41
size
Long_t size
Definition: cwb_condor_check.C:66
hp
wavearray< double > hp
Definition: DrawInspiral.C:43
gskymap::canvas
TCanvas * canvas
Definition: gskymap.hh:212
gnetwork::siteP
std::vector< TLatex * > siteP
Definition: gnetwork.hh:121
m
int m
Definition: cwb_net.C:28
skymap::max
double max()
Definition: skymap.cc:438
j
int j
Definition: cwb_net.C:28
i
i drho i
Definition: cwb_epparameters.C:88
detectorParams::longitude
double longitude
Definition: detector.hh:52
detector::tau
skymap tau
Definition: detector.hh:346
network::ifoList
std::vector< detector * > ifoList
Definition: network.hh:608
gskymap::h2
TH2D * h2
`
Definition: gskymap.hh:213
gnetwork::siteM
std::vector< TMarker * > siteM
Definition: gnetwork.hh:117
network::skyENRG
wavearray< double > skyENRG
Definition: network.hh:646
add
r add(tfmap, const_cast< char *>("hchannel"))
network::nCorrEnergy
skymap nCorrEnergy
Definition: network.hh:625
ifo
char ifo[NIFO_MAX][8]
Definition: cwb1G_parameters.C:17
net
network ** net
NOISE_MDC_SIMULATION.
Definition: CWB_Plugin_HEN_BKG_Config.C:4
network::ifoListSize
size_t ifoListSize()
Definition: network.hh:431
network::netx
static int netx(double *, double, double *, double, double)
Definition: network.hh:830
network::nLikelihood
skymap nLikelihood
Definition: network.hh:622
mode
size_t mode
Definition: cwb1G_parameters.C:134
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
skymap::getSkyIndex
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:720
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:145
gnetwork::DrawNetworkDetectorIndex
void DrawNetworkDetectorIndex(double gamma, int loop, double snr=-1.0, int dpaletteId=DUMMY_PALETTE_ID, bool btitle=true)
Definition: gnetwork.cc:1252
gskymap::FillData
void FillData(int size, double *phi, double *theta, double *binc)
Definition: gskymap.cc:394
gnetwork::DrawSitesShortLabel
void DrawSitesShortLabel(Color_t tcolor=kBlack, Size_t tsize=0.052, Font_t tfont=32)
Definition: gnetwork.cc:426
Write
tlive_fix Write()
phi
float phi
Definition: cwb_condor_create_ced.C:69
detector::antenna
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
Definition: detector.cc:490
psi
float psi
Definition: cwb_condor_create_ced.C:69
gr
TGraph * gr
Definition: cwb_draw_sensitivity.C:22
hx
wavearray< double > hx
Definition: DrawInspiral.C:44
D1
detector D1
Definition: TestReadWriteDetectorObject.C:9
wavecomplex::real
double real() const
Definition: wavecomplex.hh:69
int
i() int(T_cor *100))
NET
network NET
Definition: cwb_dump_inj.C:30
Pi
double Pi
Definition: DrawPhaseShift.C:12
gx
gwavearray< double > * gx
Definition: CompareWDMvsTime.C:18
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:22
gnetwork::DrawSites
void DrawSites(Color_t mcolor=kBlack, Size_t msize=2.0, Style_t mstyle=20)
Definition: gnetwork.cc:257
l
int l
Definition: cwb_report_cbc.C:775
gnetwork::gnetwork
gnetwork()
Definition: gnetwork.hh:39
gnetwork::siteLabelType
unsigned int siteLabelType
default=1, used to scale the antenna pattern acccording the PSD contribution of each detector (at a f...
Definition: gnetwork.hh:115
speedlight
#define speedlight
Definition: watfun.hh:33
k
int k
Definition: cwb_report_prod_2.C:168
gnetwork::ndi
static double ndi(double *, double, int)
Definition: gnetwork.hh:133
gnetwork::LoadObject
void LoadObject(char *file)
Definition: gnetwork.cc:1506
network::index
wavearray< int > index
Definition: network.hh:640
F
double F
Definition: DrawPhaseShift.C:14
skymap
Definition: skymap.hh:63
detectorParams::latitude
double latitude
Definition: detector.hh:51
gnetwork::SetScale
void SetScale(TString ifo, double scale)
Definition: gnetwork.cc:1627
e
double e
Definition: CreateListEBBH.C:35
network::skyProb
wavearray< double > skyProb
Definition: network.hh:645
skymap::phi2RA
double phi2RA(double ph, double gps)
Definition: skymap.hh:212
gnetwork::scale
double scale[NIFO_MAX]
Definition: gnetwork.hh:113
ind
int * ind
Definition: cwb_fix_live_slag_missed.C:107
gnetwork::gSM
gskymap gSM
Definition: gnetwork.hh:111
gnetwork::GetDelay
double GetDelay(TString ifo1, TString ifo2, double phi, double theta)
Definition: gnetwork.cc:926
s
s s
Definition: cwb_net.C:155
gskymap::resolution
double resolution
Definition: gskymap.hh:224
nifo
int nifo
Definition: cwb_dump_events.C:50
skymap::getPhi
double getPhi(size_t i)
Definition: skymap.hh:164
gnetwork::siteT
std::vector< TText * > siteT
Definition: gnetwork.hh:119
CelestialToCwb
void CelestialToCwb(double ilongitude, double ilatitude, double &olongitude, double &olatitude, double gps=0)
Definition: skycoord.hh:437
gps
double gps
Definition: cwb_report_prod_2.C:212
network::nAlignment
skymap nAlignment
Definition: network.hh:620
network::print
void print()
Definition: network.cc:8183
gnetwork::MakeNetworkDetectorIndex
void MakeNetworkDetectorIndex(double gamma, int loop, double snr=-1.0)
Definition: gnetwork.cc:1147
gskymap::ProjectParabolic
void ProjectParabolic(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
Definition: gskymap.cc:893
network::nDisbalance
skymap nDisbalance
Definition: network.hh:627
detector::Name
char Name[16]
Definition: detector.hh:327
fabs
double fabs(const Complex &x)
Definition: numpy.cc:55
gskymap::ProjectHammer
void ProjectHammer(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
Definition: gskymap.cc:836
network::gamma
double gamma
Definition: network.hh:592
CwbToGeographic
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:414
cnt
int cnt
Definition: cwb_condor_cleanup.C:37
network::skyMaskCC
wavearray< double > skyMaskCC
Definition: network.hh:642
S
Meyer< double > S(1024, 2)
network::nNetIndex
skymap nNetIndex
Definition: network.hh:626
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
network::nSensitivity
skymap nSensitivity
list of wdm tranformations
Definition: network.hh:619
detector::getTau
double getTau(double, double)
param: source theta,phi angles in degrees
Definition: detector.cc:698
gnetwork::ClearSites
void ClearSites()
Definition: gnetwork.cc:524
gnetwork::siteA
std::vector< TPolyLine * > siteA
Definition: gnetwork.hh:118
gnetwork::ClearSitesLabel
void ClearSitesLabel()
Definition: gnetwork.cc:548
gnetwork::GetAntennaPattern
double GetAntennaPattern(double phi, double theta, double psi=0., int polarization=1)
Definition: gnetwork.cc:564
btitle
bool btitle
Definition: DrawGnetwork2.C:16
gskymap::ProjectSinusoidal
void ProjectSinusoidal(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
Definition: gskymap.cc:875
wavearray::data
DataType_t * data
Definition: wavearray.hh:319
gnetwork::Draw
void Draw(char *smName, network *net=NULL)
Definition: gnetwork.hh:77
gnetwork::ClearSitesArms
void ClearSitesArms()
Definition: gnetwork.cc:536
SetXYZ
xyz SetXYZ(xgc, ygc, zgc)
gnetwork::ClearCircles
void ClearCircles()
Definition: gnetwork.cc:510
gnetwork::grP
TGraphPolar * grP[2]
canvas used for the polargrams
Definition: gnetwork.hh:128
network::nNullEnergy
skymap nNullEnergy
Definition: network.hh:623
network::skyMask
wavearray< short > skyMask
Definition: network.hh:641
GeographicToCwb
void GeographicToCwb(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:421
network::delay
void delay(double theta, double phi)
Definition: network.cc:7931
network::nCorrelation
skymap nCorrelation
Definition: network.hh:621
gskymap::goff
bool goff
Definition: gskymap.hh:225
detParms
detectorParams detParms[4]
Definition: DrawInspiralUserDet.C:19
skymap::size
size_t size()
Definition: skymap.hh:136
gnetwork::siteL
std::vector< TLatex * > siteL
Definition: gnetwork.hh:120
y
wavearray< double > y
Definition: Test10.C:31
gnetwork::Init
void Init()
Definition: gnetwork.cc:111
network::delta
double delta
Definition: network.hh:591
gskymap::projection
TString projection
Definition: gskymap.hh:230
gnetwork::daxisM
std::vector< TMarker * > daxisM
Definition: gnetwork.hh:123
network::r90_POL
wavearray< double > r90_POL[2]
buffer for standard response 00 ampl
Definition: network.hh:662
network::nProbability
skymap nProbability
Definition: network.hh:631
pD
detector ** pD
Definition: cwb_dump_events.C:52
exit
exit(0)
GAMMA
#define GAMMA(TYPE)
Definition: xroot.hh:5
network::r00_POL
wavearray< double > r00_POL[2]
buffer for projection on network plane 90 ampl
Definition: network.hh:661
network::nEllipticity
skymap nEllipticity
Definition: network.hh:629