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