Logo coherent WaveBurst  
Config Reference Guide 		Logo
MakePlotEgw.C
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 #define nINJ_MAX 40
20 #define nSET_MAX 40
21 
22 // -------------------------------------------------------------------------------------
23 // IFAR=50
24 // cd report/dump
25 // ln -s ../postprod/M1.V_cat2LH_hvetoLH_cat3LH.S_bin1_cut.S_bin2_cut.R_rMRA_cat2_hveto_cat3_i0cc00_i0rho0_win01_ifar50_freq32_992/data/fit_parameters_ALL.txt efit_parameters_ifar50_O1_12Sep19Jan_C01_SIM_LALBRSTFRAME_LF_rMRA_run0a_watWP10_run18_bins_cut.txt
26 // -------------------------------------------------------------------------------------
27 
28 
29 #define MODE "Exclusive"
30 
31 #define IFAR 50
32 #define IFILE "report/dump/efit_parameters_ifar50_O1_12Sep19Jan_C01_SIM_LALBRSTFRAME_LF_rMRA_run0a_watWP10_run18_bins_cut.txt"
33 #define MDC_INJ_FILE "input/O1_LALBRSTFRAME_LF_L1H1_wnb_Egw.inj"
34 
35 
36 void PlotEgw(int ninj, double* freq, double* egw, int nset, size_t* iset, TString* set);
37 
38 void MakePlotEgw() {
39 
40  double C = watconstants::SpeedOfLightInVacuo(); // m s^-1
41  double G = watconstants::GravitationalConstant(); // N m^2 kg^-2 (N=Kg m s^-2)
42  double Mo = watconstants::SolarMass(); // Kg
43  double Pc = watconstants::Parsec(); // m
44  double Pi = TMath::Pi();
45 
46 
47  // read injection file types
48  char imdc_set[nINJ_MAX][128]; // injection set
49  size_t imdc_type[nINJ_MAX]; // injection type
50  char imdc_name[nINJ_MAX][128]; // injection name
51  double imdc_fcentral[nINJ_MAX]; // injection central frequencies
52  double imdc_fbandwidth[nINJ_MAX]; // injection bandwidth frequencies
53  size_t imdc_index[nINJ_MAX]; // type reference array
54  size_t imdc_iset[nINJ_MAX]; // injection set index
55  double imdc_egw[nINJ_MAX]; // injection Egw
56 
57  int ninj=ReadInjType(MDC_INJ_FILE,nINJ_MAX,imdc_set,imdc_type,
58  imdc_name,imdc_fcentral,imdc_fbandwidth);
59 
60  TString imdc_sset[nINJ_MAX];
61  for(int i=0;i<ninj;i++) imdc_sset[i] = imdc_set[i];
62 
63  TString* imdc_set_name = new TString[ninj];
64  int nset=0;
65  for(int i=0;i<ninj;i++) {
66  bool bnew=true;
67  for(int j=0;j<nset;j++) if(imdc_set[i]==imdc_set_name[j]) bnew=false;
68  if(bnew) imdc_set_name[nset++]=imdc_set[i];
69  }
70  cout << "nset : " << nset << endl;
71  for(int i=0;i<nset;i++) {
72  for(int j=0;j<ninj;j++) if(imdc_set[j]==imdc_set_name[i]) imdc_iset[j]=i;
73  }
74 
75  int ecount[nINJ_MAX];
76  TString piumeno[nINJ_MAX];
77  float chi2[nINJ_MAX], err[nINJ_MAX], par1[nINJ_MAX], par2[nINJ_MAX], par3[nINJ_MAX];
78  double hrss50[nINJ_MAX];
79  TString ewaveform[nINJ_MAX];
80 
81  ifstream in;
82  in.open(IFILE,ios::in);
83  if (!in.good()) {cout << "Error Opening File : " << IFILE << endl;exit(1);}
84 
85  for(int k=0; k<ninj; k++) {
86 
87  in >> ecount[k] >> chi2[k] >> hrss50[k] >> piumeno[k]
88  >> err[k] >> par1[k] >> par2[k] >> par3[k] >> ewaveform[k];
89  if (!in.good()) break;
90 
91 // cout << k << "\t" << ewaveform[k] << "\t" << hrss50[k] << "\t" << imdc_fcentral[k] << "\t" << imdc_name[k] << endl;
92 
93  double Ro = 10000*Pc; // m
94  double Fo = imdc_fcentral[k]; // s^-1
95  double Ho = hrss50[k] ; // s^(1/2)
96  double Egw = pow(Pi,2) * pow(C,3)/G * pow(Ro*Fo*Ho,2);
97  cout << imdc_name[k] << "\tEgw(Mo) : " << Egw/(Mo*C*C)
98  << "\thrss50 : " << hrss50[k] << "\t" << imdc_fcentral[k] << endl;
99 
100  imdc_egw[k] = Egw/(Mo*C*C);
101  }
102 
103  PlotEgw(ninj, imdc_fcentral, imdc_egw, nset, imdc_iset, imdc_sset);
104 
105  exit(0);
106 }
107 
108 void PlotEgw(int ninj, double* freq, double* egw, int nset, size_t* iset, TString* set) {
109 
110  gStyle->SetTitleOffset(1.0,"X");
111  gStyle->SetTitleOffset(1.2,"Y");
112  gStyle->SetLabelFont(42,"X");
113  gStyle->SetLabelFont(42,"Y");
114  gStyle->SetTitleFont(42,"X");
115  gStyle->SetTitleFont(42,"Y");
116 
117  gStyle->SetTitleH(0.050);
118  gStyle->SetTitleW(0.95);
119  gStyle->SetTitleY(0.98);
120  gStyle->SetTitleFont(12,"D");
121  gStyle->SetTitleColor(kBlue,"D");
122  gStyle->SetTextFont(12);
123  gStyle->SetTitleFillColor(kWhite);
124  gStyle->SetLineColor(kWhite);
125  gStyle->SetNumberContours(256);
126  gStyle->SetCanvasColor(kWhite);
127 
128  // remove the red box around canvas
129  gStyle->SetFrameBorderMode(0);
130  gROOT->ForceStyle();
131 
132  TCanvas* canvas;
133 
134  TGraphErrors* gr50[nINJ_MAX];
135  TMultiGraph* mg;
136  TLegend* legend;
137 
138  TString set_name[nINJ_MAX];
139 
140  mg = new TMultiGraph();
141  //legend = new TLegend(0.6676955,0.1163636,0.8919753,0.4436364,NULL,"brNDC");
142  legend = new TLegend(0.1111111,0.5563636,0.3353909,0.8836364,NULL,"brNDC");
143 
144  for(int k=0;k<nset;k++) {
145 
146  gr50[k] = new TGraphErrors();
147  gr50[k]->SetLineColor(2);
148  gr50[k]->SetLineWidth(1);
149  gr50[k]->SetLineStyle(7);
150 
151  if(k==0) gr50[k]->SetMarkerColor(kBlack);
152  if(k==1) gr50[k]->SetMarkerColor(kRed);
153  if(k==2) gr50[k]->SetMarkerColor(kGreen);
154  if(k==3) gr50[k]->SetMarkerColor(kBlue);
155  if(k==4) gr50[k]->SetMarkerColor(kMagenta);
156 
157  if(k==0) gr50[k]->SetMarkerStyle(20);
158  if(k==1) gr50[k]->SetMarkerStyle(21);
159  if(k==2) gr50[k]->SetMarkerStyle(22);
160  if(k==3) gr50[k]->SetMarkerStyle(23);
161  if(k==4) gr50[k]->SetMarkerStyle(28);
162 
163  if(k==0) gr50[k]->SetMarkerSize(1.4);
164  if(k==1) gr50[k]->SetMarkerSize(1.2);
165  if(k==2) gr50[k]->SetMarkerSize(1.5);
166  if(k==3) gr50[k]->SetMarkerSize(1.5);
167  if(k==4) gr50[k]->SetMarkerSize(1.5);
168 
169  for(int i=0; i<ninj; i++) {
170  if(iset[i]==k) gr50[k]->SetPoint(i,freq[i],egw[i]);
171  if(iset[i]==k) set_name[k]=set[i];
172  }
173  mg->Add(gr50[k]);
174  }
175 
176  canvas = new TCanvas("Egw","Egw",125,82,976,576);
177  canvas->Clear();
178  canvas->ToggleEventStatus();
179  canvas->SetLogy();
180  canvas->SetLogx(true);
181  canvas->SetGridx();
182  canvas->SetGridy();
183  canvas->SetFillColor(kWhite);
184  canvas->cd();
185 
186  mg->SetTitle(TString::Format("Egw vs Frequency : %s - IFAR = %d",MODE,IFAR));
187  mg->Paint("alp");
188  mg->GetHistogram()->GetXaxis()->SetTitle("Frequency (Hz)");
189  mg->GetHistogram()->GetXaxis()->CenterTitle(true);
190  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
191  mg->GetHistogram()->GetXaxis()->SetTitleFont(42);
192  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
193  mg->GetHistogram()->GetYaxis()->SetTitleFont(42);
194  mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.20);
195  mg->GetYaxis()->SetNdivisions(10);
196  mg->GetHistogram()->GetYaxis()->SetTitle("Egw (Mo)");
197  mg->GetHistogram()->GetXaxis()->SetRangeUser(32,1024+256);
198 #if (IFAR == 8)
199  mg->GetHistogram()->GetYaxis()->SetRangeUser(2e-10,2e-6);
200 #endif
201 #if (IFAR == 100)
202  mg->GetHistogram()->GetYaxis()->SetRangeUser(5e-10,5e-6);
203 #endif
204 
205  mg->Draw("ap");
206 
207  legend->SetBorderSize(1);
208  legend->SetTextAlign(22);
209  legend->SetTextFont(12);
210  legend->SetLineColor(1);
211  legend->SetLineStyle(1);
212  legend->SetLineWidth(1);
213  legend->SetFillColor(0);
214  legend->SetFillStyle(1001);
215  legend->SetTextSize(0.04);
216  legend->SetLineColor(kBlack);
217  legend->SetFillColor(kWhite);
218 
219  for(int k=0;k<nset;k++) {
220  legend->AddEntry(gr50[k],set_name[k].Data(),"lp");
221  }
222 
223  legend->Draw();
224 
225 
226  TString ofname = IFILE;
227  ofname.ReplaceAll(".txt",".gif");
228  cout << ofname << endl;
229  canvas->Print(ofname);
230  char cmd[1024];
231  TString pfname(ofname);
232  pfname.ReplaceAll(".gif",".png");
233  sprintf(cmd,"convert %s %s",ofname.Data(),pfname.Data());
234  cout << cmd << endl;
235  gSystem->Exec(cmd);
236  sprintf(cmd,"rm %s",ofname.Data());
237  cout << cmd << endl;
238  gSystem->Exec(cmd);
239 }
nINJ_MAX
#define nINJ_MAX
Definition: MakePlotEgw.C:19
MDC_INJ_FILE
#define MDC_INJ_FILE
Definition: MakePlotEgw.C:33
MODE
#define MODE
Definition: MakePlotEgw.C:29
IFAR
#define IFAR
Definition: MakePlotEgw.C:31
IFILE
#define IFILE
Definition: MakePlotEgw.C:32
PlotEgw
void PlotEgw(int ninj, double *freq, double *egw, int nset, size_t *iset, TString *set)
Definition: MakePlotEgw.C:108
MakePlotEgw
void MakePlotEgw()
Definition: MakePlotEgw.C:38
canvas
TCanvas * canvas
Definition: Make_PP_IFAR.C:166
cmd
char cmd[1024]
Definition: cwb_mklinks.C:27
exit
exit(0)
sprintf
sprintf(tag,"wave_%s", data_label)
C
static const double C
Definition: CreateListEBBH.C:7
G
static const double G
Definition: CreateListEBBH.C:5