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