Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends 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=8
6 // -------------------------------------------------------------------------------------
7 
8 #define IFAR 8
9 
10 //#define MODE "Inclusive"
11 //#define IFILE "ifit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST1_LF_rMRA_run0a.txt"
12 //#define MDC_INJ_FILE "O1_BRST1_LF_L1H1.inj"
13 //#define IFILE "ifit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST2_LF_rMRA_run0a.txt"
14 //#define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
15 //#define IFILE "ifit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST3_LF_rMRA_run0a.txt"
16 //#define MDC_INJ_FILE "O1_BRST3_LF_L1H1.inj"
17 
18 #define MODE "Exclusive"
19 //#define IFILE "efit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST1_LF_rMRA_run0a.txt"
20 //#define MDC_INJ_FILE "O1_BRST1_LF_L1H1.inj"
21 //#define IFILE "efit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST2_LF_rMRA_run0a.txt"
22 //#define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
23 //#define IFILE "efit_parameters_ifar8_ER8b_12Sep20Oct_C01_SIM1_BRST3_LF_rMRA_run0a.txt"
24 //#define MDC_INJ_FILE "O1_BRST3_LF_L1H1.inj"
25 #define IFILE "xfit_parameters_ifar8_ER8b_12Sep20Oct_C0101_SIM1_LALBRST2_LF_rMRA_run0a.txt"
26 #define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
27 
28 
29 // -------------------------------------------------------------------------------------
30 // IFAR=100
31 // -------------------------------------------------------------------------------------
32 /*
33 #define IFAR 100
34 
35 //#define MODE "Inclusive"
36 //#define IFILE "ifit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST1_LF_rMRA_run0a.txt"
37 //#define MDC_INJ_FILE "O1_BRST1_LF_L1H1.inj"
38 //#define IFILE "ifit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST2_LF_rMRA_run0a.txt"
39 //#define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
40 //#define IFILE "ifit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST3_LF_rMRA_run0a.txt"
41 //#define MDC_INJ_FILE "O1_BRST3_LF_L1H1.inj"
42 
43 #define MODE "Exclusive"
44 //#define IFILE "efit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST1_LF_rMRA_run0a.txt"
45 //#define MDC_INJ_FILE "O1_BRST1_LF_L1H1.inj"
46 //#define IFILE "efit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST2_LF_rMRA_run0a.txt"
47 //#define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
48 //#define IFILE "efit_parameters_ifar100_ER8b_12Sep20Oct_C01_SIM1_BRST3_LF_rMRA_run0a.txt"
49 //#define MDC_INJ_FILE "O1_BRST3_LF_L1H1.inj"
50 #define IFILE "xfit_parameters_ifar100_ER8b_12Sep20Oct_C0101_SIM1_LALBRST2_LF_rMRA_run0a.txt"
51 #define MDC_INJ_FILE "O1_BRST2_LF_L1H1.inj"
52 */
53 
54 void PlotEgw(int ninj, double* freq, double* egw, int nset, size_t* iset, TString* set);
55 
56 void MakePlotEgw() {
57 
58  double C = watconstants::SpeedOfLightInVacuo(); // m s^-1
59  double G = watconstants::GravitationalConstant(); // N m^2 kg^-2 (N=Kg m s^-2)
60  double Mo = watconstants::SolarMass(); // Kg
61  double Pc = watconstants::Parsec(); // m
62  double Pi = TMath::Pi();
63 
64 
65  // read injection file types
66  char imdc_set[nINJ_MAX][128]; // injection set
67  size_t imdc_type[nINJ_MAX]; // injection type
68  char imdc_name[nINJ_MAX][128]; // injection name
69  double imdc_fcentral[nINJ_MAX]; // injection central frequencies
70  double imdc_fbandwidth[nINJ_MAX]; // injection bandwidth frequencies
71  size_t imdc_index[nINJ_MAX]; // type reference array
72  size_t imdc_iset[nINJ_MAX]; // injection set index
73  double imdc_egw[nINJ_MAX]; // injection Egw
74 
75  int ninj=ReadInjType(MDC_INJ_FILE,nINJ_MAX,imdc_set,imdc_type,
76  imdc_name,imdc_fcentral,imdc_fbandwidth);
77 
78  TString imdc_sset[nINJ_MAX];
79  for(int i=0;i<ninj;i++) imdc_sset[i] = imdc_set[i];
80 
81  TString* imdc_set_name = new TString[ninj];
82  int nset=0;
83  for(int i=0;i<ninj;i++) {
84  bool bnew=true;
85  for(int j=0;j<nset;j++) if(imdc_set[i]==imdc_set_name[j]) bnew=false;
86  if(bnew) imdc_set_name[nset++]=imdc_set[i];
87  }
88  cout << "nset : " << nset << endl;
89  for(int i=0;i<nset;i++) {
90  for(int j=0;j<ninj;j++) if(imdc_set[j]==imdc_set_name[i]) imdc_iset[j]=i;
91  }
92 
93  int ecount[nINJ_MAX];
94  TString piumeno[nINJ_MAX];
95  float chi2[nINJ_MAX], err[nINJ_MAX], par1[nINJ_MAX], par2[nINJ_MAX], par3[nINJ_MAX];
96  double hrss50[nINJ_MAX];
97  TString ewaveform[nINJ_MAX];
98 
99  ifstream in;
100  in.open(IFILE,ios::in);
101  if (!in.good()) {cout << "Error Opening File : " << IFILE << endl;exit(1);}
102 
103  for(int k=0; k<ninj; k++) {
104 
105  in >> ecount[k] >> chi2[k] >> hrss50[k] >> piumeno[k]
106  >> err[k] >> par1[k] >> par2[k] >> par3[k] >> ewaveform[k];
107  if (!in.good()) break;
108 
109 // cout << k << "\t" << ewaveform[k] << "\t" << hrss50[k] << "\t" << imdc_fcentral[k] << "\t" << imdc_name[k] << endl;
110 
111  double Ro = 10000*Pc; // m
112  double Fo = imdc_fcentral[k]; // s^-1
113  double Ho = hrss50[k] ; // s^(1/2)
114  double Egw = pow(Pi,2) * pow(C,3)/G * pow(Ro*Fo*Ho,2);
115  cout << imdc_name[k] << "\tEgw(Mo) : " << Egw/(Mo*C*C)
116  << "\thrss50 : " << hrss50[k] << "\t" << imdc_fcentral[k] << endl;
117 
118  imdc_egw[k] = Egw/(Mo*C*C);
119  }
120 
121  PlotEgw(ninj, imdc_fcentral, imdc_egw, nset, imdc_iset, imdc_sset);
122 
123 }
124 
125 void PlotEgw(int ninj, double* freq, double* egw, int nset, size_t* iset, TString* set) {
126 
127  gStyle->SetTitleOffset(1.0,"X");
128  gStyle->SetTitleOffset(1.2,"Y");
129  gStyle->SetLabelFont(42,"X");
130  gStyle->SetLabelFont(42,"Y");
131  gStyle->SetTitleFont(42,"X");
132  gStyle->SetTitleFont(42,"Y");
133 
134  gStyle->SetTitleH(0.050);
135  gStyle->SetTitleW(0.95);
136  gStyle->SetTitleY(0.98);
137  gStyle->SetTitleFont(12,"D");
138  gStyle->SetTitleColor(kBlue,"D");
139  gStyle->SetTextFont(12);
140  gStyle->SetTitleFillColor(kWhite);
141  gStyle->SetLineColor(kWhite);
142  gStyle->SetNumberContours(256);
143  gStyle->SetCanvasColor(kWhite);
144 
145  // remove the red box around canvas
146  gStyle->SetFrameBorderMode(0);
147  gROOT->ForceStyle();
148 
149  TCanvas* canvas;
150 
151  TGraphErrors* gr50[nINJ_MAX];
152  TMultiGraph* mg;
153  TLegend* legend;
154 
155  TString set_name[nINJ_MAX];
156 
157  mg = new TMultiGraph();
158  //legend = new TLegend(0.6676955,0.1163636,0.8919753,0.4436364,NULL,"brNDC");
159  legend = new TLegend(0.1111111,0.5563636,0.3353909,0.8836364,NULL,"brNDC");
160 
161  for(int k=0;k<nset;k++) {
162 
163  gr50[k] = new TGraphErrors();
164  gr50[k]->SetLineColor(2);
165  gr50[k]->SetLineWidth(1);
166  gr50[k]->SetLineStyle(7);
167 
168  if(k==0) gr50[k]->SetMarkerColor(kBlack);
169  if(k==1) gr50[k]->SetMarkerColor(kRed);
170  if(k==2) gr50[k]->SetMarkerColor(kGreen);
171  if(k==3) gr50[k]->SetMarkerColor(kBlue);
172  if(k==4) gr50[k]->SetMarkerColor(kMagenta);
173 
174  if(k==0) gr50[k]->SetMarkerStyle(20);
175  if(k==1) gr50[k]->SetMarkerStyle(21);
176  if(k==2) gr50[k]->SetMarkerStyle(22);
177  if(k==3) gr50[k]->SetMarkerStyle(23);
178  if(k==4) gr50[k]->SetMarkerStyle(28);
179 
180  if(k==0) gr50[k]->SetMarkerSize(1.4);
181  if(k==1) gr50[k]->SetMarkerSize(1.2);
182  if(k==2) gr50[k]->SetMarkerSize(1.5);
183  if(k==3) gr50[k]->SetMarkerSize(1.5);
184  if(k==4) gr50[k]->SetMarkerSize(1.5);
185 
186  for(int i=0; i<ninj; i++) {
187  if(iset[i]==k) gr50[k]->SetPoint(i,freq[i],egw[i]);
188  if(iset[i]==k) set_name[k]=set[i];
189  }
190  mg->Add(gr50[k]);
191  }
192 
193  canvas = new TCanvas("Egw","Egw",125,82,976,576);
194  canvas->Clear();
195  canvas->ToggleEventStatus();
196  canvas->SetLogy();
197  canvas->SetLogx(true);
198  canvas->SetGridx();
199  canvas->SetGridy();
200  canvas->SetFillColor(kWhite);
201  canvas->cd();
202 
203  mg->SetTitle(TString::Format("Egw vs Frequency : %s - IFAR = %d",MODE,IFAR));
204  mg->Paint("alp");
205  mg->GetHistogram()->GetXaxis()->SetTitle("Frequency (Hz)");
206  mg->GetHistogram()->GetXaxis()->CenterTitle(true);
207  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
208  mg->GetHistogram()->GetXaxis()->SetTitleFont(42);
209  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
210  mg->GetHistogram()->GetYaxis()->SetTitleFont(42);
211  mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.20);
212  mg->GetYaxis()->SetNdivisions(10);
213  mg->GetHistogram()->GetYaxis()->SetTitle("Egw (Mo)");
214  mg->GetHistogram()->GetXaxis()->SetRangeUser(32,1024+256);
215 #if (IFAR == 8)
216  mg->GetHistogram()->GetYaxis()->SetRangeUser(2e-10,2e-6);
217 #endif
218 #if (IFAR == 100)
219  mg->GetHistogram()->GetYaxis()->SetRangeUser(5e-10,5e-6);
220 #endif
221 
222  mg->Draw("ap");
223 
224  legend->SetBorderSize(1);
225  legend->SetTextAlign(22);
226  legend->SetTextFont(12);
227  legend->SetLineColor(1);
228  legend->SetLineStyle(1);
229  legend->SetLineWidth(1);
230  legend->SetFillColor(0);
231  legend->SetFillStyle(1001);
232  legend->SetTextSize(0.04);
233  legend->SetLineColor(kBlack);
234  legend->SetFillColor(kWhite);
235 
236  for(int k=0;k<nset;k++) {
237  legend->AddEntry(gr50[k],set_name[k].Data(),"lp");
238  }
239 
240  legend->Draw();
241 
242 
243  TString ofname = IFILE;
244  ofname.ReplaceAll(".txt",".gif");
245  cout << ofname << endl;
246  canvas->Print(ofname);
247  char cmd[256];
248  TString pfname(ofname);
249  pfname.ReplaceAll(".gif",".png");
250  sprintf(cmd,"convert %s %s",ofname.Data(),pfname.Data());
251  cout << cmd << endl;
252  gSystem->Exec(cmd);
253  sprintf(cmd,"rm %s",ofname.Data());
254  cout << cmd << endl;
255  gSystem->Exec(cmd);
256 }
C
static const double C
Definition: GNGen.cc:10
ewaveform
cout<< "nset : "<< nset<< endl;for(int i=0;i< nset;i++){for(int j=0;j< ninj;j++) if(imdc_set[j]==imdc_set_name[i]) imdc_iset[j]=i;}for(int iset=0;iset< nset;iset++) cout<< iset<< " "<< imdc_set_name[iset].Data()<< endl;char etitle[256];char ofile[256];TCanvas *canvas[NTYPE_MAX];int ecount[NINJ_MAX];TString piumeno[NINJ_MAX];float chi2[NINJ_MAX], err[NINJ_MAX], par1[NINJ_MAX], par2[NINJ_MAX], par3[NINJ_MAX];double ehrss10[NINJ_MAX], ehrss50[NINJ_MAX], ehrss90[NINJ_MAX];double hrss50_bis[NINJ_MAX];TString ewaveform[NINJ_MAX];TF1 *fFit[NINJ_MAX];double hrss50[NTYPE_MAX][NINJ_MAX], hrss90[NTYPE_MAX][NINJ_MAX], hrss10[NTYPE_MAX][NINJ_MAX];double inf=simulation==2?log10(factors[0]):-25;double sup=simulation==2?log10(factors[nfactor-1]):-18.5;if(simulation==1 &&pp_factor2distance){inf=log10(pp_factor2distance/factors[nfactor-1]);sup=log10(pp_factor2distance/factors[0]);}int k=0;for(int iset=0;iset< nset;iset++){char file[256];sprintf(file,"%s/fit_parameters_%s.txt", netdir, imdc_set_name[iset].Data());cout<< file<< endl;ifstream in2;in2.open(file, ios::in);if(!in2.good()){cout<< "Error Opening File : "<< file<< endl;exit(1);}for(int j=0;j< NINJ_MAX;j++){hrss50_bis[j]=0;hrss10[iset][j]=0;hrss50[iset][j]=0;hrss90[iset][j]=0;ecount[j]=0;ewaveform[j]="";}for(int l=0;l< NINJ_MAX;l++){in2 > ecount[k] chi2[k] hrss50[iset][k] piumeno[k] err[k] par1[k] par2[k] par3[k] ewaveform[k]
Definition: cwb_mkeff.C:114
cmd
char cmd[1024]
Definition: cwb_condor_cleanup.C:20
ecount
int ecount
Definition: cwb_mkhtml_prod.C:210
TString
TString("c")
Definition: cwb_report_skymap.C:111
nset
int nset
Definition: cwb_mkeff.C:57
watconstants::SolarMass
double SolarMass()
Definition: constants.hh:184
canvas
return wmap canvas
Definition: DrawClusterWithSkyplot.C:80
MakePlotEgw
void MakePlotEgw()
Definition: MakePlotEgw.C:56
j
int j
Definition: cwb_net.C:10
i
i drho i
Definition: cwb_epparameters.C:85
imdc_iset
size_t imdc_iset[NMDC_MAX]
Definition: cwb_mkeff.C:50
imdc_index
size_t imdc_index[NMDC_MAX]
Definition: cwb_mkeff.C:49
G
double G
Definition: DrawEBHH.C:12
IFILE
#define IFILE
Definition: MakePlotEgw.C:25
MODE
#define MODE
Definition: MakePlotEgw.C:18
Pi
double Pi
Definition: DrawPhaseShift.C:12
ninj
int ninj
Definition: cwb_mkeff.C:52
k
int k
Definition: cwb_report_prod_2.C:150
watconstants::Parsec
double Parsec()
Definition: constants.hh:188
chi2
float chi2
Definition: cbc_plots.C:603
imdc_name
char imdc_name[NMDC_MAX][128]
Definition: cwb_mkeff.C:46
e
double e
Definition: CreateListEBBH.C:35
ReadInjType
int ReadInjType(TString ifName, int ntype_max, char set[][128], size_t type[], char name[][128], double fcentral[], double fbandwidth[])
Definition: Toolfun.hh:740
watconstants::GravitationalConstant
double GravitationalConstant()
Definition: constants.hh:113
in
ifstream in
Definition: CWB_Plugin_HEN_BKG_Config.C:20
IFAR
#define IFAR
Definition: MakePlotEgw.C:8
nINJ_MAX
#define nINJ_MAX
Definition: MakePlotEgw.C:1
PlotEgw
void PlotEgw(int ninj, double *freq, double *egw, int nset, size_t *iset, TString *set)
Definition: MakePlotEgw.C:125
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
set
sm set(30, 10)
MDC_INJ_FILE
#define MDC_INJ_FILE
Definition: MakePlotEgw.C:26
imdc_type
size_t imdc_type[NMDC_MAX]
Definition: cwb_mkeff.C:45
imdc_set_name
TString * imdc_set_name
Definition: cwb_mkeff.C:56
imdc_fcentral
double imdc_fcentral[NMDC_MAX]
Definition: cwb_mkeff.C:47
mg
TMultiGraph * mg
Definition: cwb_draw_sensitivity.C:5
imdc_set
char imdc_set[NMDC_MAX][128]
Definition: cwb_mkeff.C:44
imdc_fbandwidth
double imdc_fbandwidth[NMDC_MAX]
Definition: cwb_mkeff.C:48
exit
exit(0)
watconstants::SpeedOfLightInVacuo
double SpeedOfLightInVacuo()
Definition: constants.hh:96