Logo coherent WaveBurst  
Library Reference Guide
Logo
cwb_mkeff.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 // used by the cwb_mkeff command
20 {
21 
22  #define NINJ_MAX 50
23  #define NTYPE_MAX 20
24  #define NMDC_MAX 64
25 
26  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_ROOTLOGON_FILE"));
27  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_PARAMETERS_FILE"));
28  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_UPARAMETERS_FILE"));
29  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_PPARAMETERS_FILE"));
30  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_UPPARAMETERS_FILE"));
31  CWB::Toolbox::checkFile(gSystem->Getenv("CWB_EPPARAMETERS_FILE"));
32 
33  gStyle->SetTitleOffset(1.0,"X");
34  gStyle->SetTitleOffset(1.2,"Y");
35 // gStyle->SetLabelOffset(0.014,"X");
36 // gStyle->SetLabelOffset(0.010,"Y");
37  gStyle->SetLabelFont(42,"X");
38  gStyle->SetLabelFont(42,"Y");
39  gStyle->SetTitleFont(42,"X");
40  gStyle->SetTitleFont(42,"Y");
41 // gStyle->SetLabelSize(0.03,"X");
42 // gStyle->SetLabelSize(0.03,"Y");
43 
44  gStyle->SetTitleH(0.050);
45  gStyle->SetTitleW(0.95);
46  gStyle->SetTitleY(0.98);
47  gStyle->SetTitleFont(12,"D");
48  gStyle->SetTitleColor(kBlue,"D");
49  gStyle->SetTextFont(12);
50  gStyle->SetTitleFillColor(kWhite);
51  gStyle->SetLineColor(kWhite);
52  gStyle->SetNumberContours(256);
53  gStyle->SetCanvasColor(kWhite);
54  gStyle->SetStatBorderSize(1);
55  gStyle->SetOptStat(kFALSE);
56 
57  // remove the red box around canvas
58  gStyle->SetFrameBorderMode(0);
59  gROOT->ForceStyle();
60 
61  // read injection file types
62  char imdc_set[NMDC_MAX][128]; // injection set
63  size_t imdc_type[NMDC_MAX]; // injection type
64  char imdc_name[NMDC_MAX][128]; // injection name
65  double imdc_fcentral[NMDC_MAX]; // injection central frequencies
66  double imdc_fbandwidth[NMDC_MAX]; // injection bandwidth frequencies
67  size_t imdc_index[NMDC_MAX]; // type reference array
68  size_t imdc_iset[NMDC_MAX]; // injection set index
69 
70  int ninj=ReadInjType(mdc_inj_file,NMDC_MAX,imdc_set,imdc_type,
71  imdc_name,imdc_fcentral,imdc_fbandwidth);
72  if(ninj==0) {cout << "Error - no injection - terminated" << endl;exit(1);}
73 
75  int nset=0;
76  for(int i=0;i<ninj;i++) {
77  bool bnew=true;
78  for(int j=0;j<nset;j++) if(imdc_set[i]==imdc_set_name[j]) bnew=false;
79  if(bnew) imdc_set_name[nset++]=imdc_set[i];
80  }
81  cout << "nset : " << nset << endl;
82 
83  for(int i=0;i<nset;i++) {
84  for(int j=0;j<ninj;j++) if(imdc_set[j]==imdc_set_name[i]) imdc_iset[j]=i;
85  }
86  for (int iset=0;iset<nset;iset++) cout << iset << " " << imdc_set_name[iset].Data() << endl;
87 
88  char etitle[256];
89  char ofile[256];
90  TCanvas* canvas[NTYPE_MAX];
91  int ecount[NINJ_MAX];
92  TString piumeno[NINJ_MAX];
93  float chi2[NINJ_MAX], err[NINJ_MAX], par1[NINJ_MAX], par2[NINJ_MAX], par3[NINJ_MAX];
94  double ehrss10[NINJ_MAX], ehrss50[NINJ_MAX], ehrss90[NINJ_MAX];
95  double hrss50_bis[NINJ_MAX];
97 
98  TF1* fFit[NINJ_MAX];
100 
101  double inf = simulation==2 ? log10(factors[0]) : -25;
102  double sup = simulation==2 ? log10(factors[nfactor-1]) : -18.5;
103 
104  if(simulation==1 && pp_factor2distance) {
105  inf = log10(pp_factor2distance/factors[nfactor-1]);
106  sup = log10(pp_factor2distance/factors[0]);
107  }
108 
109  int k=0;
110  for (int iset=0;iset<nset;iset++) {
111 
112  char file[256];
113  sprintf(file,"%s/fit_parameters_%s.txt",netdir, imdc_set_name[iset].Data());
114  cout << file << endl;
115 
116  ifstream in2;
117  in2.open(file,ios::in);
118  if (!in2.good()) {cout << "Error Opening File : " << file << endl;exit(1);}
119 
120  for (int j=0; j<NINJ_MAX; j++) {
121  hrss50_bis[j]=0;
122  hrss10[iset][j]=0;
123  hrss50[iset][j]=0;
124  hrss90[iset][j]=0;
125  ecount[j]=0;
126  ewaveform[j]="";
127  }
128 
129  for (int l=0; l<NINJ_MAX; l++) {
130 
131  in2 >> ecount[k] >> chi2[k] >> hrss50[iset][k] >> piumeno[k]
132  >> err[k] >> par1[k] >> par2[k] >> par3[k] >> ewaveform[k];
133  if (!in2.good()) break;
134  cout << ewaveform[k].Data() << endl;
135  double par0=TMath::Log10(hrss50[iset][k]);
136  fFit[k] = new TF1("logNfit",logNfit,pow(10.0,inf),pow(10.0,sup),5);
137  fFit[k]->SetNpx(100000);
138  fFit[k]->SetParameters(par0,par1[k],par2[k],par3[k],pp_factor2distance);
139  hrss10[iset][k]=fFit[k]->GetX(.1,pow(10.0,inf),pow(10.0,sup));
140  hrss50_bis[k] =fFit[k]->GetX(.5,pow(10.0,inf),pow(10.0,sup));
141  hrss90[iset][k]=fFit[k]->GetX(.9,pow(10.0,inf),pow(10.0,sup));
142  if(fFit[k]->Eval(hrss90[iset][k])<0.89) hrss90[iset][k]=1e-10;
143 
144  par0=TMath::Log10(hrss50[iset][k]+err[k]);
145  fFit[k]->SetParameters(par0,par1[k],par2[k],par3[k],pp_factor2distance);
146  ehrss10[k]=fabs(fFit[k]->GetX(.1,pow(10.0,inf),pow(10.0,sup))-hrss10[iset][k]);
147  ehrss50[k]=fabs(fFit[k]->GetX(.5,pow(10.0,inf),pow(10.0,sup))-hrss50[iset][k]);
148  ehrss90[k]=fabs(fFit[k]->GetX(.9,pow(10.0,inf),pow(10.0,sup))-hrss90[iset][k]);
149 
150  cout << hrss10[iset][k] << " " << hrss50[iset][k] << " " << hrss90[iset][k] << endl;
151 
152  k++;
153  }
154  }
155 
156  TGraphErrors* gr10[NTYPE_MAX];
157  TGraphErrors* gr50[NTYPE_MAX];
158  TGraphErrors* gr90[NTYPE_MAX];
159  TMultiGraph* mg[NTYPE_MAX];
160  TLegend* legend[NTYPE_MAX];
161 
162  for(int iset=0;iset<nset;iset++) {
163 
164  sprintf(etitle,"%s",imdc_set_name[iset].Data());
165 
166  mg[iset] = new TMultiGraph();
167  legend[iset] = new TLegend(0.732,0.125,0.956,0.354,NULL,"brNDC");
168 
169 
170  gr10[iset] = new TGraphErrors();
171  gr10[iset]->SetLineColor(1);
172  gr10[iset]->SetLineWidth(1);
173  gr10[iset]->SetMarkerColor(1);
174  gr10[iset]->SetMarkerStyle(20);
175  gr10[iset]->SetLineStyle(7);
176 
177  gr50[iset] = new TGraphErrors();
178  gr50[iset]->SetLineColor(2);
179  gr50[iset]->SetLineWidth(1);
180  gr50[iset]->SetMarkerColor(2);
181  gr50[iset]->SetMarkerStyle(20);
182  gr50[iset]->SetLineStyle(7);
183 
184  gr90[iset] = new TGraphErrors();
185  gr90[iset]->SetLineColor(kBlue);
186  gr90[iset]->SetLineWidth(1);
187  gr90[iset]->SetMarkerColor(kBlue);
188  gr90[iset]->SetMarkerStyle(20);
189  gr90[iset]->SetLineStyle(7);
190 
191  int npoint=0;
192  for(int i=0; i<ninj; i++) {
193  if(imdc_iset[i]!=iset) continue; // skip unwanted injection types
194  if(pp_show_eff_fit_curve) { // set hrss@10/50/90 from efficiency fits
195  gr10[iset]->SetPoint(npoint,imdc_fcentral[i],hrss10[iset][i]);
196  gr50[iset]->SetPoint(npoint,imdc_fcentral[i],hrss50[iset][i]);
197  gr90[iset]->SetPoint(npoint,imdc_fcentral[i],hrss90[iset][i]);
198  } else {
199  gr10[iset]->SetPoint(npoint,imdc_fcentral[i],0.);
200  gr50[iset]->SetPoint(npoint,imdc_fcentral[i],0.);
201  gr90[iset]->SetPoint(npoint,imdc_fcentral[i],0.);
202  }
203  npoint++;
204  }
205 
206  mg[iset]->Add(gr10[iset]);
207  mg[iset]->Add(gr50[iset]);
208  mg[iset]->Add(gr90[iset]);
209 
210  char namecanvas[8];
211  sprintf(namecanvas,"c%i",iset);
212  canvas[iset] = new TCanvas(namecanvas,namecanvas,125+iset*10,82,976,576);
213  canvas[iset]->Clear();
214  canvas[iset]->ToggleEventStatus();
215  canvas[iset]->SetLogy();
216 #ifdef CWB_MKEFF_LINX
217  canvas[iset]->SetLogx(false);
218 #else
219  canvas[iset]->SetLogx(true);
220 #endif
221  canvas[iset]->SetGridx();
222  canvas[iset]->SetGridy();
223  canvas[iset]->SetFillColor(kWhite);
224  canvas[iset]->cd();
225 
226  mg[iset]->SetTitle(etitle);
227  mg[iset]->Paint("alp");
228  mg[iset]->GetHistogram()->GetXaxis()->SetTitle("Frequency (Hz)");
229  mg[iset]->GetHistogram()->GetXaxis()->CenterTitle(true);
230  mg[iset]->GetHistogram()->GetXaxis()->SetLabelFont(42);
231  mg[iset]->GetHistogram()->GetXaxis()->SetTitleFont(42);
232  mg[iset]->GetHistogram()->GetYaxis()->SetLabelFont(42);
233  mg[iset]->GetHistogram()->GetYaxis()->SetTitleFont(42);
234  if(simulation==1 && pp_factor2distance) {
235  mg[iset]->GetHistogram()->GetYaxis()->SetTitle("distance (Kpc)");
236  mg[iset]->GetHistogram()->GetYaxis()->SetRangeUser(pp_factor2distance/factors[nfactor-1],pp_factor2distance/factors[0]);
237  } else if(simulation==2) {
238  mg[iset]->GetHistogram()->GetYaxis()->SetTitle("snr");
239  mg[iset]->GetHistogram()->GetYaxis()->SetRangeUser(factors[0],factors[nfactor-1]);
240  } else {
241  mg[iset]->GetHistogram()->GetYaxis()->SetTitle("hrss (1/Hz^{-1/2})");
242  mg[iset]->GetHistogram()->GetYaxis()->SetRangeUser(pp_hrss_min,pp_hrss_max);
243  }
244 
245  mg[iset]->Draw("alp");
246 
247  legend[iset]->SetBorderSize(1);
248  legend[iset]->SetTextAlign(22);
249  legend[iset]->SetTextFont(12);
250  legend[iset]->SetLineColor(1);
251  legend[iset]->SetLineStyle(1);
252  legend[iset]->SetLineWidth(1);
253  legend[iset]->SetFillColor(0);
254  legend[iset]->SetFillStyle(1001);
255  legend[iset]->SetTextSize(0.04);
256  legend[iset]->SetLineColor(kBlack);
257  legend[iset]->SetFillColor(kWhite);
258 
259  if(simulation==1 && pp_factor2distance) {
260  legend[iset]->AddEntry(gr90[iset],"dstance (Kpc) 90%","lp");
261  legend[iset]->AddEntry(gr50[iset],"dstance (Kpc) 50%","lp");
262  legend[iset]->AddEntry(gr10[iset],"dstance (Kpc) 10%","lp");
263  } else if(simulation==2) {
264  legend[iset]->AddEntry(gr90[iset],"snr 90%","lp");
265  legend[iset]->AddEntry(gr50[iset],"snr 50%","lp");
266  legend[iset]->AddEntry(gr10[iset],"snr 10%","lp");
267  } else {
268  legend[iset]->AddEntry(gr90[iset],"Hrss 90%","lp");
269  legend[iset]->AddEntry(gr50[iset],"Hrss 50%","lp");
270  legend[iset]->AddEntry(gr10[iset],"Hrss 10%","lp");
271  }
272  legend[iset]->Draw();
273 
274  sprintf(ofile,"%s/eff_freq_%s.gif",netdir,imdc_set_name[iset].Data());
275  cout << ofile << endl;
276  canvas[iset]->SaveAs(ofile);
277  }
278 
279  exit(0);
280 }
char ofile[1024]
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:132
int ecount
TString("c")
int nset
Definition: cwb_mkeff.C:75
return wmap canvas
double pp_hrss_min
int j
Definition: cwb_net.C:28
double pp_hrss_max
i drho i
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:4670
#define NMDC_MAX
size_t imdc_iset[NMDC_MAX]
Definition: cwb_mkeff.C:68
#define NINJ_MAX
size_t imdc_index[NMDC_MAX]
Definition: cwb_mkeff.C:67
char mdc_inj_file[1024]
Definition: cwb_dump_inj.C:98
exit(0)
char netdir[1024]
int ninj
Definition: cwb_mkeff.C:70
int l
hrss10[iset][k]
Definition: cwb_mkeff.C:139
int k
char imdc_name[NMDC_MAX][128]
Definition: cwb_mkeff.C:64
double e
int ReadInjType(TString ifName, int ntype_max, char set[][128], size_t type[], char name[][128], double fcentral[], double fbandwidth[])
Definition: Toolfun.hh:816
hrss90[iset][k]
Definition: cwb_mkeff.C:141
#define NTYPE_MAX
ifstream in
double fabs(const Complex &x)
Definition: numpy.cc:55
hrss50_bis[k]
Definition: cwb_mkeff.C:140
int nfactor
Definition: test_config1.C:83
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
size_t imdc_type[NMDC_MAX]
Definition: cwb_mkeff.C:63
simulation
Definition: cwb_eced.C:26
factors[0]
Definition: cwb_eced.C:27
TString * imdc_set_name
Definition: cwb_mkeff.C:74
double imdc_fcentral[NMDC_MAX]
Definition: cwb_mkeff.C:65
TMultiGraph * mg
Double_t logNfit(Double_t *x, Double_t *par)
Definition: Toolfun.hh:42
char imdc_set[NMDC_MAX][128]
Definition: cwb_mkeff.C:62
double imdc_fbandwidth[NMDC_MAX]
Definition: cwb_mkeff.C:66