Logo coherent WaveBurst  
Library Reference Guide
Logo
ProbabilityIntegratedOverCircles.C
Go to the documentation of this file.
1 // Integrate likelihood skymap over the big circles
2 // Author : Gabriele Vedovato
3 //
4 
5 {
6  //#define SKYMAP_FILE "data/ced_968654036_60_sg554q8d9_obj_20_job1/L1H1V1_968654042.750_968654042.750_968654042.750/probability.root"
7  #define SKYMAP_FILE "/home/vedovato/waveburst/coherence/offline/SIMULATION/ADV_SIM_SGQ9_L1H1V1_run5/data/ced_931158378_44_ADV_SIM_SGQ9_L1H1V1_run5_30_job1/L1H1V1_931158395.102_931158395.102_931158395.102/probability.root"
8 
9  #define SpeedOfLightInVacuo 2.99792458000000000e+08
10 
11  #define ALPHA_FACTOR 4
12  #define BETA_FACTOR 4
13 
14  #define H1L1_AXIS
15 
16  using namespace ROOT::Math;
17 
18  gSystem->Load("libMathCore");
19 
21 
22 
23  int nIFO=3;
24  TString ifo[3]={"L1","H1","V1"};
25  gnetwork* gNET = new gnetwork(nIFO,ifo);
26  gskymap* gSM = gNET->GetGskymap();
27  gSM->SetOptions("","Geographic");
28  *gSM = *iSM;
29 
30  skymap* osm[3];
31  for(int n=0;n<nIFO;n++) {
32  osm[n] = new skymap(*gSM);
33  int L = osm[n]->size();
34  for(int l=0;l<L;l++) osm[n]->set(l,0);
35  }
36 
37  detector *pD[3];
38  for(int n=0; n<nIFO; n++) pD[n]=gNET->getifo(n);
39 
40  XYZVector L1(pD[0]->Rv[0],pD[0]->Rv[1],pD[0]->Rv[2]);
41  XYZVector H1(pD[1]->Rv[0],pD[1]->Rv[1],pD[1]->Rv[2]);
42  XYZVector V1(pD[2]->Rv[0],pD[2]->Rv[1],pD[2]->Rv[2]);
43 
47 
48  XYZVector D1 = H1-L1;
49  XYZVector D3 = D1.Cross(V1-L1);
50  XYZVector D2 = D1.Cross(D3);
51 
52  XYZVector eD1 = D1.Unit();
53  XYZVector eD2 = D2.Unit();
54  XYZVector eD3 = D3.Unit();
55 
56  Rotation3D R(eD1,eD2,eD3);
57 
58  XYZVector H1L1 = H1-L1;
59  H1L1.SetRho(1);
60  XYZVector V1L1 = V1-L1;
61  V1L1.SetRho(1);
62  XYZVector V1H1 = V1-H1;
63  V1H1.SetRho(1);
64 
65  TVector3 vH1L1(H1L1.X(),H1L1.Y(),H1L1.Z());
66  TVector3 vV1L1(V1L1.X(),V1L1.Y(),V1L1.Z());
67  TVector3 vV1H1(V1H1.X(),V1H1.Y(),V1H1.Z());
68  TVector3 veD3(eD3.X(),eD3.Y(),eD3.Z());
69 
70  double r2d = 180./TMath::Pi();
71 
72  double _PI=TMath::Pi();
73 
74  int N=360*BETA_FACTOR;
75  wavearray<double> th(N);
77  double dalpha = _PI/180./(double)ALPHA_FACTOR;
78  for(int n=0;n<nIFO;n++) {
79  for(int k=0;k<ALPHA_FACTOR*180;k++) {
80 
81  double alpha = k*dalpha;
82  TRotation oR;
83  oR.Rotate(alpha,veD3.Unit());
84  TVector3 roSD;
85  if(n==0) roSD = oR*vH1L1;
86  if(n==1) roSD = oR*vV1L1;
87  if(n==2) roSD = oR*vV1H1;
88 
89  double like=0;
90  for (int i=0;i<N;i++) {
91 
92  double angle = TMath::Pi()*i/180./BETA_FACTOR;
93  TRotation vR;
94  if(n==0) vR.Rotate(angle,vH1L1.Unit());
95  if(n==1) vR.Rotate(angle,vV1L1.Unit());
96  if(n==2) vR.Rotate(angle,vV1H1.Unit());
97 
98  TVector3 rvSD = vR*roSD;
99  double theta = rvSD.Theta()*r2d;
100  double phi = rvSD.Phi()*r2d;
101  if (phi<0) phi+=360; if (phi>360) phi-=360;
102  if (theta<0) theta+=180; if (theta>180) theta-=180;
103 
104  th[i]=theta;
105  ph[i]=phi;
106  }
107 
108  for (int i=0;i<N;i++) {
109  int l=gSM->getSkyIndex(th[i],ph[i]);
110  like+=gSM->get(l);
111  }
112  like/=N;
113 //like*=alpha;
114 //if(alpha>_PI/2) like*=(_PI-alpha);
115  for (int i=0;i<N;i++) {
116  int l=osm[n]->getSkyIndex(th[i],ph[i]);
117  osm[n]->set(l,like);
118  }
119  }
120  }
121 
122  int L = gSM->size();
123  for(int l=0;l<L;l++) {
124  double like=0;
125  for (int n=0;n<nIFO;n++) like+=osm[n]->get(l);
126  //like=osm[0]->get(l);
127  gSM->set(l,like);
128  }
129 
130  gSM->Draw(0);
131 }
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
#define SKYMAP_FILE
int n
Definition: cwb_net.C:28
double angle
TString("c")
void set(size_t i, double a)
Definition: gskymap.hh:128
float theta
TH2F * ph
std::vector< double > vR
void Draw(int dpaletteId=1, Option_t *option="colfz")
Definition: gskymap.cc:460
i drho i
#define N
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:720
skymap * osm[3]
float phi
detector V1((char *)"V1")
Definition: DrawGNET.C:11
detector D1
double Pi
int l
gNET add & L1
Definition: DrawGNET.C:9
int k
Definition: skymap.hh:63
detector * pD[3]
gskymap * GetGskymap()
Definition: gnetwork.hh:44
#define BETA_FACTOR
#define SpeedOfLightInVacuo
void set(size_t i, double a)
param: sky index param: value to set
Definition: skymap.hh:122
detector H1((char *)"H1")
Definition: DrawGNET.C:10
double get(size_t i)
param: sky index
Definition: skymap.cc:699
#define ALPHA_FACTOR
size_t size()
Definition: skymap.hh:136
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:84