71 #define SHOW_3SIGMA // plot 3 sigma probability line 72 #define PROB_3SIGMA (1./370.398) // percentage outside 3 sigma gaussian probability 86 TGraph*
GetGraph(vector<double>
x, vector<double>
y,
96 double refline=0,
TString reflabel=
"reference",
TString pfile=
"") {
99 if(pfile==
"batch") {gROOT->SetBatch(
true);pfile=
"";batch=
true;}
101 if(pfile!=
"" && !pfile.EndsWith(
".png")) {
103 cout <<
"cwb_mkcfad : Error in pfile name : " << pfile << endl;
104 cout <<
" file name must ends with .png" << endl << endl;
118 unsigned int FAP_TYPE[
MAX_FAP];
123 in.open(fapList.Data(),
ios::in);
124 if (!in.good()) {cout <<
"Error Opening File : " << fapList.Data() << endl;
exit(1);}
130 in.getline(str,1024);
131 if (!in.good())
break;
132 if(str[0] !=
'#') size++;
134 cout <<
"size " << size << endl;
135 in.clear(ios::goodbit);
136 in.seekg(0, ios::beg);
137 if (size==0) {cout <<
"Error : File " << fapList.Data() <<
" is empty" << endl;
exit(1);}
139 char iFARvsRHO[1024];
140 char iOBSTIMEvsFAR[1024];
141 char iFADvsRHO[1024];
142 char iPRODvsFAD[1024];
143 char iFAP_LABEL[1024];
150 in.getline(line,1024);
152 std::stringstream linestream(line);
153 if(
TString(line).BeginsWith(
"#"))
continue;
154 if(
TString(line)==
"")
continue;
155 if(
TString(line).BeginsWith(
"$SEARCH")) {
156 if(FAP_TYPE[nFAP]>1) nFAP++;
157 linestream >> iFAP_TAG >> iFAP_LABEL >> iFAP_COLOR >> iFAP_STYLE;
158 FAP_LABEL[nFAP]=iFAP_LABEL;
159 FAP_COLOR[nFAP]=iFAP_COLOR;
160 FAP_STYLE[nFAP]=iFAP_STYLE;
163 if(
TString(line).BeginsWith(
"$FAR")) {
164 linestream >> iFAP_TAG >> iFARvsRHO >> iOBSTIMEvsFAR;
165 FARvsRHO[nFAP]=iFARvsRHO;
166 OBSTIMEvsFAR[nFAP]=iOBSTIMEvsFAR;
169 if(
TString(line).BeginsWith(
"$FAD")) {
170 linestream >> iFAP_TAG >> iFADvsRHO >> iPRODvsFAD;
171 FADvsRHO[nFAP]=iFADvsRHO;
172 PRODvsFAD[nFAP]=iPRODvsFAD;
176 if(FAP_TYPE[nFAP]>1) nFAP++;
179 for(
int i=0;
i<nFAP;
i++) {
180 cout <<
"SEARCH : " << FAP_LABEL[
i] <<
" " << FAP_COLOR[
i] <<
" " << FAP_STYLE[
i] << endl;
181 if(FAP_TYPE[
i]&2) cout << FARvsRHO[
i] <<
" " << OBSTIMEvsFAR[
i] << endl;
182 if(FAP_TYPE[i]&4) cout << FADvsRHO[
i] <<
" " << PRODvsFAD[
i] << endl;
184 for(
int j=i;
j<nFAP;
j++)
if(FAP_TYPE[i]!=FAP_TYPE[
j]) {
185 cout <<
"cwb_mkcfap : Error - searches not consistent !!!" << endl;
186 cout <<
" Missing FAR or FAD declaration" << endl << endl;
189 if((!fapMode.Contains(
"RHO"))&&(!fapMode.Contains(
"FAR"))&&(!fapMode.Contains(
"FAD"))) {
190 cout <<
"cwb_mkcfap : Error - fapMode not contains valid declarations (FAD/FAR/RHO)" << endl;
193 if((!FAP_TYPE[i]&2)&&(fapMode.Contains(
"RHO"))) {
194 cout <<
"cwb_mkcfap : Error - fapMode=RHO needs FAR declarations" << endl;
197 if((!FAP_TYPE[i]&2)&&(fapMode.Contains(
"FAR"))) {
198 cout <<
"cwb_mkcfap : Error - fapMode=FAR needs FAR declarations" << endl;
201 if((!FAP_TYPE[i]&2)&&(fapMode.Contains(
"FAD"))) {
202 cout <<
"cwb_mkcfap : Error - fapMode=FAD needs FAD declarations" << endl;
207 if((fapMode.Contains(
"RHO"))&&(!fapMode.Contains(
"FAR"))&&(!fapMode.Contains(
"FAD"))) {
208 for(
int i=0;
i<nFAP;
i++) FAP_COLOR[
i]=2;
213 cout << endl <<
"---------------------------------------------------------------------" << endl;
214 for(
int i=0;
i<nFAP;
i++) {
216 if(FAP_TYPE[
i]&2) FAP =
GetFAP(refline, i, nFAP, FARvsRHO, OBSTIMEvsFAR);
217 if(FAP_TYPE[i]&4) FAP =
GetFAP(refline, i, nFAP, FADvsRHO, PRODvsFAD);
218 cout <<
"Search : " << FAP_LABEL[
i] <<
"\t -> FAP at RHO = " << refline <<
" is " << FAP << endl;
220 cout <<
"---------------------------------------------------------------------" << endl << endl;
226 gCANVAS =
new TCanvas(
"canvas",
"canvas",16,30,825,546);
227 gCANVAS->Range(-19.4801,-9.25,-17.4775,83.25);
228 gCANVAS->SetBorderSize(2);
229 gCANVAS->SetFrameFillColor(0);
236 gStyle->SetOptFit(kTRUE);
246 for(
int i=0;
i<nFAP;
i++) {
247 if(fapMode.Contains(
"FAD")) {
250 int xsize = xFAP[
i].size();
251 if(xsize) {
for(
int j=0;
j<xsize;
j++)
if(yFAD[i][
j]<yFAP[i][
j]) yFAP[
i][
j]=yFAD[
i][
j];}
252 else {xFAP[
i]=xFAD[
i];yFAP[
i]=yFAD[
i];}
254 if(fapMode.Contains(
"FAR")) {
256 GetFAPvsRHO(
i, xFAR[
i], yFAR[i], nFAP, FARvsRHO, OBSTIMEvsFAR);
258 int xsize = xFAP[
i].size();
259 if(xsize) {
for(
int j=0;
j<xsize;
j++)
if(yFAR[i][
j]<yFAP[i][
j]) yFAP[
i][
j]=yFAR[
i][
j];}
260 else {xFAP[
i]=xFAR[
i];yFAP[
i]=yFAR[
i];}
262 if(fapMode.Contains(
"RHO")) {
264 GetFAPvsRHO(xRHO[
i], yRHO[i], nFAP, FARvsRHO, OBSTIMEvsFAR);
266 int xsize = xFAP[
i].size();
267 if(xsize) {
for(
int j=0;
j<xsize;
j++)
if(yRHO[i][
j]<yFAP[i][
j]) yFAP[
i][
j]=yRHO[
i][
j];}
268 else {xFAP[
i]=xRHO[
i];yFAP[
i]=yRHO[
i];}
272 for(
int i=0;
i<nFAP;
i++) {
275 if(fapMode.Contains(
"FAD")) {subtitle+=
" FAD ";nMODE++;}
276 if(fapMode.Contains(
"FAR")) {subtitle+=
" FAR ";nMODE++;}
277 if(fapMode.Contains(
"RHO")) {subtitle+=
" RHO ";nMODE++;}
280 title = TString::Format(
"False Alarm Probability - ranked by %s",subtitle.Data());
282 title = TString::Format(
"False Alarm Probability - minimum FAP ranked by (%s)",subtitle.Data());
285 gr[
i] =
GetGraph(xFAP[
i],yFAP[i],title,
"IFAD",
"FAP");
287 gr[
i] =
GetGraph(xFAP[i],yFAP[i],title,
"rho",
"FAP");
289 gr[
i]->SetLineColor(FAP_COLOR[i]);
290 gr[
i]->SetLineStyle(FAP_STYLE[i]);
291 if(xFAP[i][xFAP[i].
size()-1] > rho_max) rho_max=xFAP[
i][xFAP[
i].size()-1];
295 TMultiGraph*
mg =
new TMultiGraph();
297 for(
int i=0;
i<nFAP;
i++) mg->Add(gr[
i]);
301 gr[nFAP] =
new TGraph;
302 gr[nFAP]->SetPoint(0,refline,0);
303 gr[nFAP]->SetPoint(1,refline,1);
304 gr[nFAP]->SetLineColor(kGreen);
305 gr[nFAP]->SetLineWidth(2);
311 gr[nFAP+1] =
new TGraph;
314 gr[nFAP+1]->SetLineColor(kBlue);
315 gr[nFAP+1]->SetLineWidth(2);
316 gr[nFAP+1]->SetLineStyle(10);
322 mg->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
323 mg->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
324 mg->GetHistogram()->GetXaxis()->SetTitleSize(0.04);
325 mg->GetHistogram()->GetYaxis()->SetTitleSize(0.04);
326 mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
328 mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
329 mg->GetHistogram()->GetYaxis()->SetLabelOffset(0.01);
330 mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.5);
333 mg->GetXaxis()->SetLabelFont(42);
334 mg->GetYaxis()->SetLabelFont(42);
335 mg->GetXaxis()->SetTitleFont(42);
336 mg->GetYaxis()->SetTitleFont(42);
337 mg->GetXaxis()->SetTitleOffset(1.20);
338 mg->GetYaxis()->SetTitleOffset(1.05);
339 mg->GetXaxis()->SetTitleSize(0.04);
340 mg->GetYaxis()->SetTitleSize(0.04);
342 mg->GetXaxis()->SetTitle(
"IFAD");
344 mg->GetXaxis()->SetTitle(
"rho");
346 mg->GetYaxis()->SetTitle(
"FAP");
348 mg->GetXaxis()->SetMoreLogLabels();
354 double hleg = 0.84-nFAP*0.03;
355 leg =
new TLegend(0.7369062,hleg,0.9914738,0.9265385,NULL,
"brNDC");
356 leg->SetBorderSize(1);
357 leg->SetTextAlign(22);
358 leg->SetTextFont(12);
359 leg->SetLineColor(1);
360 leg->SetLineStyle(1);
361 leg->SetLineWidth(1);
362 leg->SetFillColor(0);
363 leg->SetFillStyle(1001);
364 leg->SetTextSize(0.03);
365 leg->SetLineColor(kBlack);
366 leg->SetFillColor(kWhite);
368 for(
int i=0;i<nFAP;i++) leg->AddEntry(gr[i],FAP_LABEL[i].Data(),
"lp");
369 if(refline>0) leg->AddEntry(gr[nFAP],reflabel.Data(),
"lp");
371 leg->AddEntry(gr[nFAP+1],
"3 sigma",
"lp");
378 gfile.ReplaceAll(
".png",
".gif");
379 gCANVAS->Print(gfile);
381 sprintf(cmd,
"convert %s %s",gfile.Data(),pfile.Data());
384 sprintf(cmd,
"rm %s",gfile.Data());
396 double MU = FAX*PROD;
397 double FAP = 1.-exp(-MU);
412 int nrho = TMath::Nint((rho_max-rho_min)/drho)+2;
414 for(
int i=0;
i<nrho;
i++) {
415 double RHO = rho_min+
i*drho;
419 double MU = FAX*PROD;
420 double FAP = 1.-exp(-MU);
440 for(
int j=0;
j<nFAP;
j++) {
443 if(_rho_min<rho_min) rho_min=_rho_min;
444 if(_rho_max>rho_max) rho_max=_rho_max;
450 int nrho = TMath::Nint((rho_max-rho_min)/drho)+2;
454 for(
int i=0;
i<nrho;
i++) {
455 double RHO = rho_min+
i*drho;
457 for(
int j=0;
j<nFAP;
j++) {
460 if(FAR>0)
if(FAR<MIN_FAR[j]) MIN_FAR[
j]=
FAR;
461 MU += MIN_FAR[
j]*OBSTIME;
463 double FAP = 1.-exp(-MU);
478 TGraph*
gr =
new TGraph;
483 if(y[
i]==0)
continue;
484 double dx = (x[
i]-x[
i-1])/10000.;
485 gr->SetPoint(cnt++,x[
i]-dx,y[
i-1]);
486 gr->SetPoint(cnt++,x[
i]+dx,y[
i]);
489 gr->GetHistogram()->GetXaxis()->SetTitle(xtitle.Data());
490 gr->GetHistogram()->GetYaxis()->SetTitle(ytitle.Data());
491 gr->GetHistogram()->GetYaxis()->SetTitleOffset(1.4);
493 double max=0;
double min=1e80;
495 if(y[
i]==0)
continue;
496 if(y[
i]>max) max=y[
i];
if(y[
i]<min && y[
i]!=0) min=y[
i];
499 gr->GetHistogram()->GetYaxis()->SetMoreLogLabels();
501 gr->SetTitle(ptitle);
502 gr->SetLineColor(kRed);
TGraph * GetGraph(vector< double > x, vector< double > y, TString ptitle, TString xtitle, TString ytitle)
double min(double x, double y)
double GetFAP(double rho, int n, int nFAP, TString *FAXvsRHO, TString *PRODvsFAX)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
void cwb_mkcfad(TString fapList, TString fapMode, double rho_min=0, double rho_max=0, double refline=0, TString reflabel="reference", TString pfile="")
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void GetFAPvsRHO(int n, vector< double > &x, vector< double > &y, int nFAP, TString *FAXvsRHO, TString *PRODvsFAX)