Logo Coherent WaveBurst  
Reference Guide
Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CWB_Plugin_PE.C
Go to the documentation of this file.
1 #define XIFO 4
2 
3 #pragma GCC system_header
4 
5 #include "cwb.hh"
6 #include "config.hh"
7 #include "network.hh"
8 #include "wavearray.hh"
9 #include "TString.h"
10 #include "TObjArray.h"
11 #include "TObjString.h"
12 #include "TRandom.h"
13 #include "TComplex.h"
14 #include "TGraphAsymmErrors.h"
15 #include "TMath.h"
16 #include "mdc.hh"
17 #include "cwb2G.hh"
18 #include "watplot.hh"
19 #include "gwavearray.hh"
20 #include "network.hh"
21 #include <fstream>
22 #include <vector>
23 
24 // ---------------------------------------------------------------------------------
25 // TO BE DONE
26 // ---------------------------------------------------------------------------------
27 
28 // MEDIAN SKYMAP, reformat of xCED SkyMap TAB
29 // FIX WAVEFORM INJECTION in RedoAnalysis, must be excluded CAT2 , we can use TH1D for random selection
30 // ADD SVN Version to DumpUserOptions ---> (DONE)
31 // FIX output format of DumpUserOptions ---> (DONE)
32 // FIX xCED Spectrograms when pe_noise=0, only detected waveform is displayed ---> (DONE)
33 // ADD GPS time offset to the time plots ---> (DONE)
34 // ADD WDM TF maps ?
35 // ADD descriptions for each TAB in the xCED
36 // SAVE injected whitened signal ---> (DONE)
37 // SAVE all PE/waveform in the output ROOT file ---> (DONE)
38 // For simulation=0 set window +/- iwindow/2 around the event ---> (DONE)
39 // SAVE pe_config to root output file
40 // DUMP to CED dir all PE parameters ---> (DONE)
41 // ADD Checks of the input PE config parameters ---> (DONE)
42 // FIX pe_noise=2 when cedDump=false -> the TF is cleaned and whitened noise is deleted.
43 // USE pe_ced_options to select tabs to be included into xCED ---> (DONE)
44 // IMPLEMENT command : cwb_condor mtpe jobID ---> (DONE)
45 // FIX iSNR,oSNR for non simulation PE : the output SNRnet plot is zero
46 // add check on ISTAGE, PE can only be executed with FULL stage ---> (DONE)
47 // add TF PCA likelihood/null ---> (DONE)
48 
49 // ---------------------------------------------------------------------------------
50 // HOW TO CONFIGURE PLUGIN
51 // the following is an example : must be included in the config/user_parameters.C file
52 // see the DumpUserOptions function for the full description of parameters
53 // ---------------------------------------------------------------------------------
54 /*
55  plugin = TMacro(gSystem->ExpandPathName("$HOME_CWB/plugins/CWB_Plugin_PE.C")); // Macro source
56 
57  TString optpe = ""; // NOTE : add space at the end of each line
58  optpe += "pe_retry=100 "; // number of trials
59  optpe += "pe_gps=1167559936 "; // only gps +/- iwindow is analyzed
60  optpe += "pe_noise=0 "; // signal used for trials : 0/1/2 reconstructed/injected/(original=reconstructed+noise)
61  optpe += "pe_signal=0 "; // waveform is injected in the whitened HoT and apply Coherent+SuperCluster+Likelihood stage
62  optpe += "pe_amp_cal_err=0.1 "; // max percentage of amplitude miscalibration (uniform in -0.9,+1.1) -> det 1
63  optpe += "pe_phs_cal_err=10 "; // max phase (degrees) of phase miscalibration (uniform in -10,+10) -> det 1
64  optpe += "pe_amp_cal_err=0.1 "; // ... -> det 2
65  optpe += "pe_phs_cal_err=10 "; // ... -> det 2
66  optpe += "pe_ced_dump=-1 "; // dump CED at gLRETRY=PE_CED_DUMP
67  optpe += "pe_skymask=0 "; // disable -> search over all sky
68 
69  //optpe += "pe_multitask=true "; // enable/disable multitask (only for 1 job and simulation=0/4, dag is generated with cwb_condor mtpe jobID)
70  // for simulation=0 user must define : optpe += "pe_retry=100 ";
71  // for simulation=4 user must define : nfactor=xx, factor[0]=yy : where xx>0 and yy>0
72 
73  // CED options (only if cedDump=true)
74  optpe += "pe_ced_enable=tfmap ";
75  optpe += "pe_ced_enable=rdr ";
76  optpe += "pe_ced_enable=skymap ";
77  optpe += "pe_ced_enable=rec ";
78  optpe += "pe_ced_enable=inj ";
79  optpe += "pe_ced_disable=rinj ";
80  optpe += "pe_ced_enable=cm ";
81  optpe += "pe_ced_enable=distr ";
82  optpe += "pe_ced_enable=null ";
83  optpe += "pe_ced_disable=pca ";
84 
85  // OUTPUT options (only if cedDump=false)
86  optpe += "pe_output_enable=inj "; // save injection to the output root file
87  optpe += "pe_output_enable=rec "; // save reconstructed waveform to the output root file
88  optpe += "pe_output_enable=med "; // save median to the output root file
89  optpe += "pe_output_enable=p50 "; // save percentile 50 to the output root file
90  optpe += "pe_output_enable=p90 "; // save percentile 90 to the output root file
91  optpe += "pe_output_enable=avr "; // save averaged waveform to the output root file
92  optpe += "pe_output_enable=rms "; // save RMS to the output root file
93 
94  strcpy(parPlugin,optpe.Data()); // set PE plugin parameters
95  strcpy(comment,"pe configuration example");
96 
97 */
98 
99 // ---------------------------------------------------------------------------------
100 // DEFINES
101 // ---------------------------------------------------------------------------------
102 
103 #define nTRIALS 100 // number of trials for noise statistic (GetNoiseStatistic)
104 
105 #define PE_INDEX "$HOME_CWB/www/pe_index.html" // html template used for pe index report page
106 
107 //#define PLOT_TYPE "root"
108 #define PLOT_TYPE "png"
109 //#define SAVE_WAVE2SPARSE_PLOT
110 //#define SAVE_REDUCED_INJ
111 
112 #define PLOT_MEDIAN // plot median,lower50 boud,upper50 bound,lower90 boud,upper90 bound
113  // if commented we plot average,rms,2*rms
114 
115 //#define SAVE_SKYPROB // save gSKYPROB
116 #define SKYMASK_RADIUS 0.1 // used to define skymask
117 
118 #define MAX_TRIALS 1000
119 
120 #define EFRACTION_CUT 0.98 // energy threshold for time and frequency cuts
121 //#define SET_WAVEFORM_CUTS // enable time and frequency cuts
122 
123 // default pe user options
124 #define PE_TRIALS 100 // number of trials
125 #define PE_SIGNAL 0 // signal used for trials : 0/1/2 reconstructed/injected/(original=reconstructed+noise)
126  // note : option 1 can be used only in simulation mode
127 #define PE_NOISE 1 // noise used for trials : 0/1/2
128  // 0 - waveform is injected in the whitened HoT and apply Coherent+SuperCluster+Likelihood stages
129  // 1 - add gaussian noise to likelihood sparse map and apply Likelihood stage
130  // 2 - add whitened HoT to likelihood sparse map and apply Likelihood stage
131 #define PE_AMP_CAL_ERR 0 // max percentage of amplitude miscalibration if>0 (uniform in -max,+max) else (gaus(max))
132 #define PE_PHS_CAL_ERR 0 // max phase (degrees) of phase miscalibration if>0 (uniform in -max,+max) else (gaus(max))
133 #define PE_SKYMASK 0 // skymask used for trials : 0(def)/1/2
134  // 0 : disable -> search over all sky
135  // 1 : skymask with SKYMASK_RADIUS and source selected according the reconstructed skymap probability
136  // 2 : skymask select only the sky position where the waveform is reconstructed
137 #define PE_SEED 0 // 0 : seed used by PE for ramdom generation
138 #define PE_GPS 0 // default 0 - if >0 only gps +/- iwindow is analyzed
139 
140 #define PE_MULTITASK false // if enabled, PE trials are executed by different jobs in multitask : only simulation=0 and single event (gps>0)
141 
142 #define NOISE_SIGMA 1.0 // noise sigma used in AddNoiseAndCalErrToSparse
143 
144 #define PE_CED_DUMP -1 // dump CED at gLRETRY=PE_CED_DUMP
145 
146 #define PE_CED_TFMAP true // Shows the Time-Frequency Maps -> Spectrograms, Scalograms, TF Like,Null
147 #define PE_CED_RDR true // Shows Reconstructed Detector Response tab
148 #define PE_CED_PSD true // Shows Power Spectral Density tab
149 #define PE_CED_SKYMAP true // Shows SkyMaps
150 #define PE_CED_REC true // Shows Reconstructed Waveforms/Instantaneous-Frequency with errors
151 #define PE_CED_INJ true // Shows Injected vs Reconstructed Waveforms/Instantaneous-Frequency with errors
152 #define PE_CED_rINJ false // Shows Reduced Injected vs Reconstructed Waveforms with errors
153 #define PE_CED_CM true // Shows Chirp Mass Value/Error Distributions
154 #define PE_CED_DISTR false // Shows Residuals Distributions
155 #define PE_CED_NULL true // Shows Null Pixels Distributions
156 #define PE_CED_PCA false // Shows PCA likelihood TF plot
157 
158 #define PE_OUTPUT_INJ false // save injected into the output root file
159 #define PE_OUTPUT_REC false // save reconstructed into the output root file
160 #define PE_OUTPUT_WHT false // save whitened data into the output root file
161 #define PE_OUTPUT_DAT false // save rec+nois data into the output root file
162 #define PE_OUTPUT_MED false // save median into the output root file
163 #define PE_OUTPUT_P50 false // save percentile 50 into the output root file
164 #define PE_OUTPUT_P90 false // save percentile 90 into the output root file
165 #define PE_OUTPUT_AVR false // save averaged into the output root file
166 #define PE_OUTPUT_RMS false // save rms into the output root file
167 
168 // ---------------------------------------------------------------------------------
169 // FUNCTIONS
170 // ---------------------------------------------------------------------------------
171 
172 void ClearVectors();
174 
175 double GetSparseMapData(SSeries<double>* SS, bool phase, int index);
176 
177 void SetOutputFile(network* NET, netevent* &EVT, CWB::config* cfg, bool dump_pe);
178 void DumpOutputFile(network* NET, netevent* &EVT, CWB::config* cfg, int ID, int k, int factor);
180 void CreateIndexHTML(TString dirCED, int nIFO, TString* ifo, bool sim=false);
181 void WriteBodyHTML(TString html_template, TString html_tag_beg, TString html_tag_end, ofstream* out, int nIFO=1, TString* ifo=NULL);
182 
183 void Wave2Sparse(network* NET, CWB::config* cfg, char type);
185 int RedoAnalysis(TFile* jfile, CWB::config* cfg, network* NET);
186 void ReplaceSuperclusterData(TFile*& jfile, CWB::config* cfg, network* NET, double gps=0);
187 
188 void GetChirpMassStatistic(std::vector<double>* vCHIRP);
189 void GetFactorsStatistic(int nIFO);
190 void GetSNRnetStatistic(int nIFO);
191 void GetNullPixels(std::vector<double>* aNUL, std::vector<double>* ANUL,
192  network* NET, CWB::config* cfg, int lag, int id);
193 void GetNoisePixels(std::vector<double>* aNSE, std::vector<double>* ANSE,
194  network* NET, CWB::config* cfg, int lag, int id);
195 std::vector<int> ComputeStatisticalErrors(network* NET, CWB::config* cfg, int ID);
196 void SaveSkyProb(network* NET, CWB::config* cfg, int id);
197 void SetEventWindow(CWB::config* cfg, double gps);
198 skymap GetSkyProb(network* NET, int id);
200 void DumpSkyProb(skymap* skyprob, network* NET, netevent* &EVT, TString odir);
201 
202 void SetWaveformCuts(wavearray<double>* x, double bT, double eT, double bF, double eF);
203 wavearray<double> GetCutWaveform(wavearray<double> x, double bT, double eT, double bF, double eF);
204 
206 double GetTimeBoundaries(wavearray<double> x, double P, double& bT, double& eT);
207 double GetFrequencyBoundaries(wavearray<double> x, double P, double& bF, double& eF);
209 
210 std::vector<netpixel> DoPCA(network* NET, CWB::config* cfg, int lag, int id);
211 void ResetPCA(network* NET, CWB::config* cfg, netcluster* pwc, std::vector<netpixel>* vPIX, int ID);
212 
213 std::vector<wavearray<double> > GetPCAWaveform(network* NET, CWB::config* cfg, int lag, int id);
214 std::vector<wavearray<double> > GetInjWaveform(network* NET, netevent* EVT, int id, double factor);
215 std::vector<wavearray<double> > GetRecWaveform(network* NET, netevent* EVT, int id);
216 std::vector<wavearray<double> > GetWhitenedData(network* NET, CWB::config* cfg);
217 std::vector<wavearray<double> > GetWaveform(network* NET, int lag, int id, char type, bool shift=true);
218 std::vector<wavearray<double> > GetSigWaveform(network* NET, CWB::config* cfg, int lag, int id);
219 
220 void SetSkyMask(network* NET, CWB::config* cfg, double theta, double phi, double radius);
221 
226 
228 
230 
233  bool fft=false, TString pdir="", double P=0.99, EColor col1=kRed, EColor col2=kBlue, EColor col3=(EColor)418);
234 void PlotWaveformErrors(TString ofname, TString title, CWB::config* cfg, wavearray<double>* wrec,
235  wavearray<double>* wavr, wavearray<double>* werr, wavearray<double>* wref, TString pdir="", double P=0.99);
238  wavearray<double>* wl90, wavearray<double>* wu90, wavearray<double>* wref, TString pdir, double P, bool freq=false);
239 void PlotFrequencyErrors(TString ofname, TString title, CWB::config* cfg, wavearray<double>* frec,
240  wavearray<double>* favr, wavearray<double>* ferr, wavearray<double>* wref, TString pdir, double P);
241 void PlotWaveforms(network* NET, CWB::config* cfg, int ID, TString pdir="");
242 void PlotResiduals(int ID, TString pdir="", int sim=true, char type='w');
243 void PlotSNRnet(int nIFO, TString pdir, int sim=true);
244 void PlotChirpMass(int gtype, TString pdir, int sim=true);
245 void PlotFactors(int gtype, int nIFO, TString pdir);
246 void GetNullStatistic(std::vector<double>* vNUL, std::vector<double>* vNSE, int ifoId, TString ifoName, TString gtype, TString pdir="");
247 void PlotTimeFrequencyPCA(network* NET, netevent* EVT, CWB::config* cfg, int id, int lag, TString pdir);
248 void PlotSpectrogram(TString type, network* NET, netevent* &EVT, CWB::config* cfg, TString pdir);
249 
250 void DumpRecWavePE(network* NET, TString pdir="");
251 void DumpInjWavePE(network* NET, TString pdir="");
252 
253 void LoadFromMultiTaskJobOutput(int runID, CWB::config* cfg);
254 
255 // ---------------------------------------------------------------------------------
256 // USER CONFIG OPTIONS
257 // ---------------------------------------------------------------------------------
258 
259 struct uoptions {
260  int trials;
261  int signal;
262  int noise;
263  float amp_cal_err[NIFO_MAX];
264  float phs_cal_err[NIFO_MAX];
265  int id;
266  int skymask;
267  int seed; // default 0
268  double gps; // default 0 - if >0 only gps +/- iwindow is analyzed
269 
270  bool multitask; // if enabled, PE trials are executed by different jobs in multitask : only simulation=0 and single event (gps>0)
271 
272  int ced_dump; // dump CED at gLRETRY=ced_dump (-1 = disable)
273 
274  bool ced_tfmap; // Shows the Time-Frequency Maps -> Spectrograms, Scalograms, TF Like,Null
275  bool ced_rdr; // Shows Reconstructed Detector Response
276  bool ced_psd; // Show Power Spectral Density
277  bool ced_skymap; // Shows SkyMaps
278  bool ced_rec; // Shows Reconstructed Waveforms/Instantaneous-Frequency with errors
279  bool ced_inj; // Shows Injected vs Reconstructed Waveforms/Instantaneous-Frequency with errors
280  bool ced_rinj; // Shows Reduced Injected vs Reconstructed Waveforms with errors
281  bool ced_cm; // Shows Chirp Mass Value/Error Distributions
282  bool ced_distr; // Shows Residuals Distributions
283  bool ced_null; // Shows Null Pixels Distributions
284  bool ced_pca; // Shows PCA likelihood TF plot
285 
286  bool output_inj; // save injected into the output root file
287  bool output_rec; // save reconstructed into the output root file
288  bool output_wht; // save whitened data into the output root file
289  bool output_dat; // save rec+nois data into the output root file
290  bool output_med; // save median into the output root file
291  bool output_p50; // save percentile 50 into the output root file
292  bool output_p90; // save percentile 90 into the output root file
293  bool output_avr; // save averaged into the output root file
294  bool output_rms; // save rms into the output root file
295 };
296 
297 void ResetUserOptions();
299 void PrintUserOptions(CWB::config* cfg);
301 
302 // ----------------------------------------------------
303 // ROOT Output PE Parameters
304 // ----------------------------------------------------
305 
306 struct PE { // structure for output for estimated parameters
307 
308  int trials; // number of effective trials
309  float erR[11]; // probability distribution of residuals
310  float erF[11]; // probability distribution of frequency residuals
311  float nstat[2*7*NIFO_MAX]; // null pixel statistic, for each detector
312  // 0->7 : Null 00
313  // 7->13 : Null 90
314  // gPE.stat[0] = Null Pixels Mean
315  // gPE.stat[1] = Null Pixels RMS
316  // gPE.stat[2] = Noise Pixels Mean
317  // gPE.stat[3] = Noise Pixels RMS
318  // gPE.stat[4] = Noise Pixels Chi2
319  // gPE.stat[5] = KolmogorovTest
320  // gPE.stat[6] = AndersonDarlingTest
321  float snet[2]; // SNRnet statistic, 0 -> avr, 1 -> rms
322  float ff[2]; // Fitting Factor statistic, 0 -> avr, 1 -> rms
323  float of[2]; // Overlap Factor statistic, 0 -> avr, 1 -> rms
324  float mch[2]; // chirp mass statistic, 0 -> avr, 1 -> rms
325 
331  wavearray<double>* wL50[NIFO_MAX];
337 };
338 
339 // ---------------------------------------------------------------------------------
340 // Global Variables
341 // ---------------------------------------------------------------------------------
342 
343 uoptions gOPT; // global User Options
344 PE gPE; // global output for estimated parameters
345 TTree* gTREE; // output tree file name
346 TString gOUTPUT; // output root file name
348 bool gCEDDUMP; // save user_parameters cedDump
349 
350 double gITHETA; // injected theta sky position (EVT->theta[1])
351 double gIPHI; // injected phy sky position (EVT->phi[1])
352 double gOTHETA; // reconstructed theta sky position (EVT->theta[3]) : the maximum of detection statistic
353 double gOPHI; // reconstructed phy sky position (EVT->phi[3]) : the maximum of detection statistic
354 
355 skymap gSKYPROB; // saved reconstructed probability skymap
356 TH1D gHSKYPROB; // used to extract random skyloc for skymask
357 
358 double gBT, gET; // time ranges used for cutted injections (EFRACTION_CUT)
359 double gBF, gEF; // freq ranges used for cutted injections (EFRACTION_CUT)
360 
361 wavearray<double> gHOT[NIFO_MAX]; // HoT time series used to save whitened data
362 
363 double gINRATE; // input data sampling rate
364 int gRATEANA; // analysis data rate
365 
366 int gMTRIAL; // trial ID, used to process a single trial when PE is executed in multitask mode (only with simulation=0)
367 int gID; // ID of the detected event, used in multitask mode to exclude from the processing the ID!=gID
368 
369 double gSEGGPS; // segment start time
370 
371 // ---------------------------------------------------------------------------------
372 // waveforms (vectors index is ifos)
373 // ---------------------------------------------------------------------------------
374 
375 std::vector<wavearray<double> > vINJ; // injected
376 std::vector<wavearray<double> > vREC; // signal reconstructed by wave packet
377 std::vector<wavearray<double> > vWHT; // whitened data (in the vREC time range)
378 std::vector<wavearray<double> > vPCA; // reconstructed by PCA
379 std::vector<wavearray<double> > vDAT; // noise+signal reconstructed by wave packet
380 std::vector<wavearray<double> > vNUL; // vDAT-vREC
381 std::vector<wavearray<double> > cINJ; // injected (cutted) on TF sub-space of the reconstructed waveform vREC
382 
383 std::vector<wavearray<double> > vAVR; // average reconstructed amplitude
384 std::vector<wavearray<double> > vRMS; // rms
385 
386 std::vector<wavearray<double> > vMED; // median reconstructed amplitude
387 std::vector<wavearray<double> > vL50; // = vMED - 50% Lower Bound
388 std::vector<wavearray<double> > vU50; // = vMED - 50% Upper Bound
389 std::vector<wavearray<double> > vL90; // = 90% Lower Bound - vMED
390 std::vector<wavearray<double> > vU90; // = 90% Upper Bound - vMED
391 
392 std::vector<wavearray<double> > fREC; // instantaneous reconstructed frequency
393 std::vector<wavearray<double> > fINJ; // instantaneous cutted injected frequency
394 std::vector<wavearray<double> > fAVR; // instantaneous averaged reconstructed frequency
395 std::vector<wavearray<double> > fRMS; // rms of instantaneous averaged frequency
396 
397 std::vector<wavearray<double> > fMED; // median instantaneous reconstructed frequency
398 std::vector<wavearray<double> > fL50; // 50% Lower Bound
399 std::vector<wavearray<double> > fU50; // 50% Upper Bound
400 std::vector<wavearray<double> > fL90; // 90% Lower Bound
401 std::vector<wavearray<double> > fU90; // 90% Upper Bound
402 
403 std::vector<double > vRES; // normalized residual energy
404 std::vector<double > fRES; // normalized frequency weighted residuals
405 
406 std::vector<wavearray<double> > wREC[MAX_TRIALS]; // reconstructed signal in the trials
407 
408 std::vector<double> iSNR[NIFO_MAX]; // input SNR^2
409 std::vector<double> oSNR[NIFO_MAX]; // reconstructed SNR^2
410 std::vector<double> ioSNR[NIFO_MAX]; // iSNR*oSNR
411 std::vector<double> vLIKELIHOOD; // likelihood
412 
413 std::vector<double> aNUL[NIFO_MAX]; // null 00 phase amplitudes
414 std::vector<double> ANUL[NIFO_MAX]; // null 90 phase amplitudes
415 
416 std::vector<double> aNSE[NIFO_MAX]; // noise 00 phase amplitudes
417 std::vector<double> ANSE[NIFO_MAX]; // noise 90 phase amplitudes
418 
419 std::vector<double> vCHIRP[6]; // 0 -> injected chirp mass
420  // 1 -> reconstructed chirp mass
421  // 2 -> reconstructed chirp mass error
422  // 3 -> ellipticity parameter
423  // 4 -> pixel fraction
424  // 5 -> energy fraction
425 
426 skymap wSKYPROB[MAX_TRIALS]; // sky probability trials
427 skymap mSKYPROB; // sky probability median
428 
429 // ---------------------------------------------------------------------------------
430 // sparse maps
431 // ---------------------------------------------------------------------------------
432 
433 std::vector<SSeries<double> > vSS[NIFO_MAX]; // original sparse maps
434 std::vector<SSeries<double> > jSS[NIFO_MAX]; // injected
435 std::vector<SSeries<double> > rSS[NIFO_MAX]; // reconstructed
436 std::vector<SSeries<double> > dSS[NIFO_MAX]; // noise+signal reconstructed by wave packet
437 
438 void
440 //!MISCELLANEA
441 // Plugin used for Parameters Estimation
442 
443  if(type==CWB_PLUGIN_CONFIG) {
444  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
445 
446  if(gCWB2G->istage!=CWB_STAGE_FULL) {cout<< "CWB_Plugin_PE - Error : PE can be executed only in CWB_STAGE_FULL mode" << endl;exit(1);}
447 
448  if(cfg->pattern<=0) {
449  cout << "CWB_Plugin_PE Error : PE enable only with pattern>0 !!!" << endl;
450  exit(1);
451  }
452  cfg->outPlugin=true; // disable built-in output root file
453  gINRATE=cfg->inRate; // input data sampling rate
454  gRATEANA=gCWB2G->rateANA; // analysis data rate
455  gID=0; // used in multitask mode
456 
457  ResetUserOptions(); // set default config options
458  ReadUserOptions(cfg->parPlugin); // user config options : read from parPlugin
459  TString cwb_inet_options=TString(gSystem->Getenv("CWB_INET_OPTIONS")); // overwrite parPlugin
460  ReadUserOptions(cwb_inet_options); // user config options : read from command line
461 
462  if(gOPT.ced_inj && gOPT.ced_rinj) gOPT.ced_rinj=false;
463 
464  if(gOPT.multitask) {
465  if(cfg->simulation==4 && cfg->nfactor!=1) {
466  cout<< "CWB_Plugin_PE - Error : when simulation=4, PE multitask can be executed only with nfactors=1" << endl;exit(1);
467  }
468  }
469 
470  if(cfg->cedDump) { // if at least one of the output options is enabled
471  // then we enable the output root file even when cedDump=true
472  if(gOPT.output_inj || gOPT.output_rec ||
473  gOPT.output_med || gOPT.output_p50 ||
474  gOPT.output_p90 || gOPT.output_avr ||
475  gOPT.output_rms || gOPT.output_wht || gOPT.output_dat) cfg->online=true;
476  }
477 
478  // get gMTRIAL
479  // NOTE : trial ID, used to process a single trial when PE is executed in multitask mode
480  // can be used only for simulation=0 and gOPT.gps>0 (single event) !!!
481  gMTRIAL=gOPT.trials;
482  if(gOPT.multitask) {
483  TString gtrial=TString(gSystem->Getenv("CWB_MDC_FACTOR"));
484  if(gtrial.CompareTo("")!=0) {
485  if(!gtrial.IsDigit()) {cout<< "CWB_Plugin_PE - Error : when simulation=0, CWB_MDC_FACTOR must be an interger > 0" << endl;exit(1);}
486  gMTRIAL=gtrial.Atoi();
487  if(gMTRIAL==0) {cout<< "CWB_Plugin_PE - Error : when simulation=0, CWB_MDC_FACTOR must be an interger > 0" << endl;exit(1);}
488  if(gMTRIAL!=gOPT.trials) { // multitask jobs
489  // add trial to data_label -> different output files for each trial
490  sprintf(cfg->data_label,"%s_trial%d",cfg->data_label,gMTRIAL);
491  cout << "CWB_Plugin_PE - MultiTask Mode : new data_label = " << cfg->data_label << endl;
492  // disable ooptions not used in multitask mode
493  cfg->cedDump=false;
494  gOPT.output_inj=false;
495  gOPT.output_rec=false;
496  gOPT.output_wht=false;
497  gOPT.output_dat=false;
498  gOPT.output_med=false;
499  gOPT.output_p50=false;
500  gOPT.output_p90=false;
501  gOPT.output_avr=false;
502  gOPT.output_rms=false;
503  } else { // last multitask job : collect all output from multitask jobs
504  LoadFromMultiTaskJobOutput(gCWB2G->runID, cfg); // read wREC,gSKYPROB from multitalk output root files
505  }
506  bool last_trial = (gMTRIAL==gOPT.trials) ? true : false;
507  if(!last_trial) {gOPT.trials=1;gMTRIAL*=-1;}
508  }
509  }
510 
511  gCEDDUMP=cfg->cedDump; // save user_parameter cedDump
512 
513  gRandom->SetSeed(gOPT.seed); // set seed for PE random generation
514 
515  SetEventWindow(cfg, gOPT.gps); // if gps>0 only gps +/- iwindow is analyzed
516 
517  // check input PE config parameters
518 
519  if(gOPT.signal!=0 && gOPT.signal!=1 && gOPT.signal!=2) {
520  cout << endl;
521  cout << "CWB_Plugin_PE Error : option pe_signal not allowed : must be (0/1/2) !!!" << endl;
522  cout << endl;
523  exit(1);
524  }
525 
526  if(cfg->simulation==0 && gOPT.signal==1) {
527  cout << endl;
528  cout << "CWB_Plugin_PE Error : option pe_signal=1 is allowed only in simulation mode !!!" << endl;
529  cout << endl;
530  exit(1);
531  }
532 
533  if(cfg->simulation==0 && gOPT.skymask==3) {
534  cout << endl;
535  cout << "CWB_Plugin_PE Error : option pe_skymask=3 is allowed only in simulation mode !!!" << endl;
536  cout << endl;
537  exit(1);
538  }
539  }
540 
541  if(type==CWB_PLUGIN_INIT_JOB) {
542  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
543 
544  if(gOPT.gps<gCWB2G->GetSegBegin() || gOPT.gps>gCWB2G->GetSegEnd()) {
545  cout.precision(14);
546  cout << "CWB_Plugin_PE - Error : gps time must be within the segment interval: "
547  << gCWB2G->GetSegBegin() << "\t" << gCWB2G->GetSegEnd() << endl;exit(1);
548  cout.precision(6);
549  }
550  }
551 
552  if(type==CWB_PLUGIN_NETWORK) {
553  PrintUserOptions(cfg); // print config options
554  }
555 
557  if(NET->mdcTime.size()>1) {
558  cout << endl;
559  cout << "CWB_Plugin_PE Error : detected " << NET->mdcTime.size() << " injections."
560  << " Must be only one injection per segment !!!" << endl;
561  cout << endl;
562  exit(1);
563  }
564  }
565 
567  // save whitened HoT
568  int ifoID =0; for(int n=0;n<cfg->nIFO;n++) if(ifo==NET->getifo(n)->Name) {ifoID=n;break;}
569  gHOT[ifoID] = *x;
570  }
571 
572  if(type==CWB_PLUGIN_ILIKELIHOOD) { // INIT CWB LIKELIHOOD STAGE (once for each lag)
573  NET->wfsave=true; // enable waveform save
574 
575  // export gLRETRY
576  if(gOPT.multitask) EXPORT(int,gLRETRY,"gLRETRY = 1;")
577  else EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",gOPT.trials).Data())
578 
579  ReplaceSuperclusterData(jfile, cfg, NET, gOPT.gps); // save in jfile max size supercluster
580  }
581 
582  if(type==CWB_PLUGIN_XLIKELIHOOD) { // BEFORE EVENT RECONSTRUCTION (for each event repeated gOPT.trials)
583 
584  wavearray<double> cid = NET->getwc(0)->get((char*)"ID", 0,'S',0); // get cluster ID
585  if(cid.size()==0) return;
586  int ID = size_t(cid.data[cid.size()-1]+0.1);
587 
588  if(gOPT.trials==0) return;
589 
590  // import gLRETRY
591  int gLRETRY=-1; IMPORT(int,gLRETRY)
592  cout << endl;
593  cout << "-------------------------------------------------------------------" << endl;
594  cout << "-----> CWB_Plugin_PE.C -> ID : " << ID
595  << " -> gLRETRY : " << gLRETRY << "/" << gOPT.trials << endl;
596  cout << "-------------------------------------------------------------------" << endl;
597 
598  // -------------------------------------------------
599  // FIRST TRIAL
600  // reconstruct signal with original noise (full sky)
601  // -------------------------------------------------
602  if(gLRETRY==gOPT.trials) {
603  ClearVectors();
604  cfg->cedDump=false;
605  }
606 
607  // -----------------------------------------------------------------------------------
608  // AFTER FIRST TRIAL
609  // reconstruct signal with reconstructed signal + gaussian noise & (calib. errors)
610  // -----------------------------------------------------------------------------------
611  if(gLRETRY!=gOPT.trials) {
612 
613  // SKYMASK
614  if(gOPT.skymask==1) {
615  // select random sky position from sky probability gHSKYPROB
616  int isky = (int)gHSKYPROB.GetRandom();
617  double th = 90-gSKYPROB.getTheta(isky);
618  double ph = gSKYPROB.getPhi(isky);
619 
620  SetSkyMask(NET, cfg, th, ph, SKYMASK_RADIUS);
621  }
622  if(gOPT.skymask==2) {
623  // select reconstructed sky position (the maximum of detection statistic : waveform is reconstructed in such position)
624  SetSkyMask(NET, cfg, gOTHETA, gOPHI, SKYMASK_RADIUS);
625  }
626  if(gOPT.skymask==3) {
627  // select injected sky position : only in simulation mode
628  SetSkyMask(NET, cfg, gITHETA, gIPHI, SKYMASK_RADIUS);
629  }
630 
631  // ADD NOISE (& ADD CALIBRATION ERRORS)
632  // add noise to the sparse map of the reconstructed event
633  // Note : calibration errors are applied only if gOPT.noise is enabled !!!
634  char jtype='r';
635  if(gOPT.signal==0) jtype='r'; // use reconstructed waveform for trials
636  if(gOPT.signal==1) jtype='j'; // use injected waveform for trials
637  if(gOPT.signal==2) jtype='v'; // use data for trials
638  if(gOPT.noise) AddNoiseAndCalErrToSparse(NET, cfg, jtype);
639 
640  // DUMP CED
641  cfg->cedDump=false;
642  if((gOPT.ced_dump>=0) && (gLRETRY==gOPT.ced_dump)) cfg->cedDump=true;
643  }
644 
645  // -------------------------------------------------
646  // LAST TRIAL
647  // reconstruct signal with original noise (full sky)
648  // -------------------------------------------------
649  if((gLRETRY==0)&&(gMTRIAL==gOPT.trials)) { // use original injection setup
650 
651  // SKYMASK
652  // select detected sky position (waveform is reconstructed in such poosition)
653  //SetSkyMask(NET, cfg, gOTHETA, gOPHI, SKYMASK_RADIUS);
654  // restore full sky search for the last trial
655  SetSkyMask(NET, cfg, 0, 0, 360); // full sky search
656 
657  if(gOPT.noise!=0) {
658  // restore original sparse maps
659  int nIFO = NET->ifoListSize(); // number of detectors
660  for(int n=0;n<nIFO;n++) {
661  detector* pD = NET->getifo(n);
662  pD->vSS=vSS[n];
663  }
664  }
665 
666  // DUMP CED
667  cfg->cedDump=gCEDDUMP; // restore user_parameter cedDump (executed in the last trial)
668  }
669  }
670 
671  if(type==CWB_PLUGIN_OLIKELIHOOD) { // AFTER EVENT RECONSTRUCTION (for each event repeated gOPT.trials)
672 
673  // import trials
674  int gLRETRY=-1; IMPORT(int,gLRETRY)
675 
676  // import ifactor
677  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
678  double factor = cfg->simulation==3||cfg->simulation==4 ? -gIFACTOR : cfg->factors[gIFACTOR];
679  double ofactor=0;
680  if(cfg->simulation==4) ofactor=-gIFACTOR;
681  else if(cfg->simulation==3) ofactor=-gIFACTOR;
682  else ofactor=factor;
683 
684  int nIFO = NET->ifoListSize(); // number of detectors
685  int K = NET->nLag; // number of time lag
686  int rate = 0; // select all resolutions
687  netevent* EVT;
689 
690  SetOutputFile(NET, EVT, cfg, false); // set output root file
691 
692  for(int k=0; k<K; k++) { // loop over the lags
693 
694  id = NET->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
695 
696  if(id.size()==0) cfg->cedDump=gCEDDUMP; // rejected -> restore user_parameter cedDump
697 
698  for(int j=0; j<(int)id.size(); j++) { // loop over cluster index
699 
700  int ID = size_t(id.data[j]+0.5);
701 
702  if(NET->getwc(k)->sCuts[ID-1]!=-1) continue; // skip rejected/processed clusters
703 
704  EVT->output2G(NULL,NET,ID,k,ofactor); // get reconstructed parameters
705 
706  for(int n=0; n<nIFO; n++) { // save SNR
707  iSNR[n].push_back(EVT->iSNR[n]);
708  oSNR[n].push_back(EVT->oSNR[n]);
709  ioSNR[n].push_back(EVT->ioSNR[n]);
710  }
711  for(int n=0; n<6; n++) { // save chirp mass parameters
712  vCHIRP[n].push_back(EVT->chirp[n]);
713  }
714  vLIKELIHOOD.push_back(EVT->likelihood); // save likelihood
715 
716  // print event parameters
717  cout << endl;
718  cout << "event parameters : ID -> " << ID << endl;
719  for(int n=0;n<nIFO;n++) printf("rec_time %s : %.4f\n",NET->ifoName[n], EVT->time[n]);
720  cout << "rec_theta : " << EVT->theta[0] << " rec_phi : " << EVT->phi[0] << endl;
721  cout << "SNRnet : " << sqrt(EVT->likelihood) << " netcc[0] : " << EVT->netcc[0]
722  << " rho[0] : " << EVT->rho[0] << " size : " << EVT->size[0] << endl;
723 
724  if(gOPT.trials==0) {
725  DumpOutputFile(NET, EVT, cfg, ID, k, ofactor); // dump event to output root file
726  continue;
727  }
728 
729  // -------------------------------------------------
730  // FIRST TRIAL
731  // reconstruct signal with original noise (full sky)
732  // -------------------------------------------------
733  if((gLRETRY==gOPT.trials) || (gOPT.multitask && gLRETRY==1)) {
734  // get waveforms
735  if(cfg->simulation) {
736  vINJ = GetInjWaveform(NET, EVT, ID, factor); // get injected waveform
737  if(vINJ.size()!=nIFO) {
738  cout << "CWB_Plugin_PE Error : Injection Waveform Not Found !!!" << endl;
739  exit(1);
740  }
741  cINJ = vINJ; // cINJ is overwritten by GetSigWaveform
742  }
743  vREC = GetRecWaveform(NET, EVT, ID); // get reconstructed waveform
744  vWHT = GetWhitenedData(NET, cfg); // get whitened data (in the vREC time range)
745  vDAT = GetWaveform(NET, k, ID, 'W'); // get reconstructed+noise waveform
746  vNUL = GetWaveform(NET, k, ID, 'N'); // get noise-reconstructed waveform
747  for(int n=0;n<nIFO;n++) { // compute the istantaneous frequency of reconstructed waveform
748  fREC.push_back(CWB::Toolbox::getHilbertIFrequency(vREC[n]));
749  }
750 
751  Wave2Sparse(NET,cfg,'v'); // save original sparse map to vSS
752  Wave2Sparse(NET,cfg,'r'); // rec to rSS
753  Wave2Sparse(NET,cfg,'d'); // dat to dSS
754  if(cfg->simulation) Wave2Sparse(NET,cfg,'j'); // inj to jSS
755  // get injected waveform on TF sub-space of the reconstructed waveform vREC
756  if(cfg->simulation) {
757  cINJ = GetSigWaveform(NET, cfg, k, ID);
758  for(int n=0;n<nIFO;n++) { // compute the istantaneous frequency of cutted injected waveform
759  fINJ.push_back(CWB::Toolbox::getHilbertIFrequency(cINJ[n]));
760  }
761  }
762 
763  if(gOPT.skymask==1) SaveSkyProb(NET,cfg,ID); // save sky probability, used for trials
764 
765  GetNullPixels(aNUL, ANUL, NET, cfg, k, ID); // get null pixels in TF domain
766  GetNoisePixels(aNSE, ANSE, NET, cfg, k, ID); // get noise pixels in TF domain
767  for(int n=0; n<nIFO; n++) { // fill gPE.nstat output
768  GetNullStatistic(&aNUL[n], &aNSE[n], n, NET->ifoName[n], "null_00");
769  GetNullStatistic(&ANUL[n], &ANSE[n], n, NET->ifoName[n], "null_90");
770  }
771 
772  // save event parameters
773  gSEGGPS = EVT->gps[0]; // save segment start time (sec)
774  gITHETA = EVT->theta[1]; gIPHI = EVT->phi[1]; // save injected sky position
775  gOTHETA = EVT->theta[3]; gOPHI = EVT->phi[3]; // save reconstructed sky poosition (the maximum of detection statistic)
776 
777  //vPCA = GetPCAWaveform(NET, cfg, k, ID); // get pca waveform
778  }
779 
780  // -----------------------------------------------------------------------------------
781  // AFTER FIRST TRIAL
782  // reconstruct signal with reconstructed signal + gaussian noise & (calib. errors)
783  // -----------------------------------------------------------------------------------
784  if(gLRETRY!=gOPT.trials) {
785  wREC[gLRETRY] = GetRecWaveform(NET, EVT, ID); // get reconstructed waveform
786  wSKYPROB[gLRETRY] = GetSkyProb(NET, ID); // get sky probability
787  gID=ID; // used in multitask mode
788  }
789 
790  // -------------------------------------------------
791  // LAST TRIAL
792  // reconstruct signal with original noise (full sky)
793  // -------------------------------------------------
794  if(gLRETRY==0) {
795  std::vector<int> nrec = ComputeStatisticalErrors(NET, cfg, ID);
796  if(cfg->simulation) GetFactorsStatistic(nIFO); // fill gPE.ff, gPE.of
797  GetSNRnetStatistic(nIFO); // fill gPE.snet
798  GetChirpMassStatistic(vCHIRP); // fill gPE.mch
799  mSKYPROB = GetMedianSkyProb(NET); // get median sky probability
800 
801  SetOutputFile(NET, EVT, cfg, true); // set output root file
802  DumpOutputFile(NET, EVT, cfg, ID, k, ofactor); // dump event to output root file
803  if(cfg->cedDump) {
804  TString dirCED = DumpCED(NET, EVT, cfg, ofactor); // dump CED
805  cout << "dirCED : " << dirCED << endl;
806  DumpSkyProb(&mSKYPROB, NET, EVT, dirCED); // dump median mSKYPROB to fits file
807  TString ifo[NIFO_MAX];
808  for(int n=0; n<nIFO; n++) ifo[n]=NET->ifoName[n];
809  CreateIndexHTML(dirCED, nIFO,ifo,cfg->simulation); // create html file
810  if(nrec.size()) {
811  PlotWaveforms(NET,cfg,ID,dirCED); // plot waveforms
812  if(gOPT.ced_distr) {
813  PlotResiduals(ID,dirCED,cfg->simulation,'w'); // plot residual enegy
814  PlotResiduals(ID,dirCED,cfg->simulation,'f'); // plot frequency residual enegy
815  gStyle->SetOptFit(1111);
816  PlotSNRnet(nIFO,dirCED,cfg->simulation); // plot SNRnet
817  }
818  if(gOPT.ced_cm) {
819  for(int n=1;n<6;n++) PlotChirpMass(n,dirCED,cfg->simulation); // plot mchirp
820  }
821  PlotFactors(0,nIFO,dirCED); // plot fitting factor
822  PlotFactors(1,nIFO,dirCED); // plot overlap factor;
823  if(gOPT.ced_null) {
824  for(int n=0; n<nIFO; n++) { // plot null statistic
825  GetNullStatistic(&aNUL[n], &aNSE[n], n, ifo[n], "null_00", dirCED);
826  GetNullStatistic(&ANUL[n], &ANSE[n], n, ifo[n], "null_90", dirCED);
827  }
828  }
829  gStyle->SetOptFit(0000);
830  DumpRecWavePE(NET,dirCED); // dumps rec waveform/time/freq/errors array in ASCII format.
831  if(cfg->simulation) DumpInjWavePE(NET,dirCED); // dumps inj waveform/time in ASCII format.
832  }
833  DumpUserOptions(dirCED, cfg); // dump PE user options
834 
835  if(gOPT.ced_pca) PlotTimeFrequencyPCA(NET, EVT, cfg, ID, k, dirCED); // plot tf likelihood pca
836  if(gOPT.ced_null) PlotSpectrogram("null", NET, EVT, cfg, dirCED); // plot null spectrograms
837  if(gOPT.ced_null) PlotSpectrogram("signal", NET, EVT, cfg, dirCED); // plot reconstructed signal spectrograms
838  if(cfg->simulation) PlotSpectrogram("injection", NET, EVT, cfg, dirCED); // plot injected signal spectrograms
839  }
840 
841  for(int n=0;n<nIFO;n++) {
842  detector* pD = NET->getifo(n);
843  pD->vSS=vSS[n]; // restore original sparse maps
844  if(!cfg->simulation) ClearWaveforms(pD); // release waveform memory
845  }
846  }
847  }
848  }
849 
850  jfile->cd();
851  if(EVT) delete EVT;
852  }
853 
854  // if gOPT.noise=0 the reconstructed is injected in the whitened data and Coherence+SuperCluster is done
855  if(type==CWB_PLUGIN_OLIKELIHOOD) {
856  if(gOPT.noise==0) {
857  int gLRETRY=-1; IMPORT(int,gLRETRY)
858  int csize=0;
859  do {
860  cout << endl;
861  cout << "-------------------------------------------------------------------" << endl;
862  cout << "-----> CWB_Plugin_PE.C -> RedoAnalysis : "
863  << " -> gLRETRY : " << gLRETRY << "/" << gOPT.trials << endl;
864  cout << "-------------------------------------------------------------------" << endl;
865  csize = RedoAnalysis(jfile, cfg, NET);
866  EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",--gLRETRY).Data())
867  } while(csize==0 && gLRETRY>0);
868  EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",++gLRETRY).Data())
869  }
870  }
871 
872  return;
873 }
874 
875 std::vector<netpixel> DoPCA(network* NET, CWB::config* cfg, int lag, int id) {
876 
877  double ee;
878 
879  size_t nIFO = NET->ifoList.size();
880 
881  float En = 2*NET->acor*NET->acor*nIFO; // network energy threshold in the sky loop
882 
883  int size = NET->a_00.size();
884  int f_ = NIFO/4;
885 
886  netpixel* pix;
887  netcluster* pwc = NET->getwc(lag);
888  std::vector<netpixel> vPIX;
889 
890  std::vector<int> pI = NET->wdmMRA.getXTalk(pwc, id, false); // buffer for pixel IDs
891  int V = pI.size();
892  if(V>cfg->BATCH) return vPIX; // attach TD amp to pixels < V
893 
894  wavearray<float> xi(size); xi=0; // PC 00 array
895  wavearray<float> XI(size); XI=0; // PC 90 array
896 
897  __m128* _xi = (__m128*) xi.data; // set pointer to PC 00 array
898  __m128* _XI = (__m128*) XI.data; // set pointer to PC 90 array
899 
900  __m128* _aa = (__m128*) NET->a_00.data; // set pointer to 00 array
901  __m128* _AA = (__m128*) NET->a_90.data; // set pointer to 90 array
902 
903  int nPC = 0;
904  for(int j=0; j<V; j++) {
905  int jf = j*f_; // source sse pointer increment
906  _sse_zero_ps(_xi+jf); // zero MRA amplitudes
907  _sse_zero_ps(_XI+jf); // zero MRA amplitudes
908  ee = _sse_abs_ps(_aa+jf,_AA+jf); // total pixel energy / quadrature
909  if(ee>En) nPC++; else ee=0.; // count core pixels
910  NET->rNRG.data[j] = ee; // init residual energy array
911  NET->pNRG.data[j] = NET->rNRG.data[j]; // init residual energy array
912  }
913 
914  nPC = NET->_sse_mra_ps(xi.data,XI.data,En,nPC); // get principal components
915 
916  for(int j=0; j<V; j++) { // loop over principal components
917  pix = pwc->getPixel(id,pI[j]);
918  vPIX.push_back(*pix); // save original pixels
919  pix->core = false;
920  ee = NET->pNRG.data[j]; // total pixel energy
921  if(ee<En) continue;
922  pix->core = true;
923  for(int i=0; i<nIFO; i++) {
924  pix->setdata(double(xi.data[j*NIFO+i]),'S',i); // store 00 whitened response PC
925  pix->setdata(double(XI.data[j*NIFO+i]),'P',i); // store 90 whitened response PC
926  }
927  }
928 
929  return vPIX;
930 }
931 
932 void ResetPCA(network* NET, CWB::config* cfg, netcluster* pwc, std::vector<netpixel>* vPIX, int ID) {
933 
934  std::vector<int> pI = NET->wdmMRA.getXTalk(pwc, ID, false); // buffer for pixel IDs
935  int V = pI.size();
936  if(V>cfg->BATCH) return; // attach TD amp to pixels < V
937  for(int j=0; j<V; j++) {
938  netpixel* pix = pwc->getPixel(ID,pI[j]);
939  *pix = (*vPIX)[j];
940  }
941 
942  while(!vPIX->empty()) vPIX->pop_back();
943  vPIX->clear(); std::vector<netpixel>().swap(*vPIX);
944 }
945 
946 std::vector<wavearray<double> > GetPCAWaveform(network* NET, CWB::config* cfg, int lag, int id) {
947 
948  std::vector<wavearray<double> > xPCA; // reconstructed
949 
950  std::vector<netpixel> vPIX;
951  vPIX = DoPCA(NET, cfg, lag, id); // do PCA analysis
952 
953  int nIFO = NET->ifoListSize(); // number of detectors
954  netcluster* pwc = NET->getwc(lag);
955 
956  if(NET->getMRAwave(id,lag,'S',0,true)) { // reconstruct whitened shifted pd->waveForm
957  detector* pd;
958  for(int i=0; i<nIFO; i++) { // loop over detectors
959  pd = NET->getifo(i);
960  wavearray<double>* wf = &pd->waveForm;
961  wf->start(pwc->start+pd->waveForm.start());
962  xPCA.push_back(*wf);
963  }
964  }
965 
966  ResetPCA(NET, cfg, pwc, &vPIX, id); // restore WP pwc
967 
968  return xPCA;
969 }
970 
971 void PlotTimeFrequencyPCA(network* NET, netevent* EVT, CWB::config* cfg, int id, int lag, TString pdir) {
972 
973  std::vector<netpixel> vPIX;
974  vPIX = DoPCA(NET, cfg, lag, id); // do PCA analysis
975 
976  int nIFO = NET->ifoListSize(); // number of detectors
977  netcluster* pwc = NET->getwc(lag);
978 
979  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",EVT->gps[0]);
980 
981  watplot WTS(const_cast<char*>("wts"));
982  WTS.canvas->cd();
983  char fname[1024];
984  sprintf(fname, "%s/pca_l_tfmap_scalogram.png",pdir.Data());
985  cout << "write " << fname << endl;
986  WTS.plot(pwc, 1, nIFO, 'L', 0, const_cast<char*>("COLZ"));
987  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[0],EVT->high[0]);
988  WTS.hist2D->GetXaxis()->SetTitle(xtitle);
989  WTS.canvas->Print(fname);
990  WTS.clear();
991 /*
992  sprintf(fname, "%s/pca_n_tfmap_scalogram.png",pdir.Data());
993  cout << "write " << fname << endl;
994  WTS.plot(pwc, 1, nIFO, 'N', 0, const_cast<char*>("COLZ"));
995  WTS.canvas->Print(fname);
996  WTS.clear();
997 */
998 
999  ResetPCA(NET, cfg, pwc, &vPIX, id); // restore WP pwc
1000 
1001 }
1002 
1003 std::vector<wavearray<double> > GetInjWaveform(network* NET, netevent* EVT, int ID, double factor) {
1004 
1005  std::vector<wavearray<double> > xINJ; // injected
1006 
1007  int nIFO = NET->ifoListSize(); // number of detectors
1008 
1009  double recTime = EVT->time[0]; // rec event time det=0
1010 
1011  injection INJ(nIFO);
1012  // get inj ID
1013  int M = NET->mdc__IDSize();
1014  double injTime = 1.e12;
1015  int injID = -1;
1016  for(int m=0; m<M; m++) {
1017  int mdcID = NET->getmdc__ID(m);
1018  double T = fabs(recTime - NET->getmdcTime(mdcID));
1019  if(T<injTime && INJ.fill_in(NET,mdcID)) {
1020  injTime = T;
1021  injID = mdcID;
1022  }
1023  }
1024 
1025  if(INJ.fill_in(NET,injID)) { // get inj parameters
1026 
1027  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
1028 
1029  // extract whitened injected & reconstructed waveforms
1030  for(int n=0; n<nIFO; n++) {
1031 
1032  detector* pd = NET->getifo(n);
1033 
1034  pwfINJ[n] = INJ.pwf[n];
1035  if (pwfINJ[n]==NULL) {
1036  cout << "CWB_Plugin_Boostrap.C : Error : Injected waveform not saved !!! : detector "
1037  << NET->ifoName[n] << endl;
1038  continue;
1039  }
1040 
1041  double rFactor = 1.;
1042  rFactor *= factor>0 ? factor : 1;
1043  wavearray<double>* wf = pwfINJ[n];
1044  *wf*=rFactor; // injected wf is multiplied by the custom factor
1045  xINJ.push_back(*wf);
1046  *wf*=1./rFactor; // restore injected amplitude
1047  }
1048  delete [] pwfINJ;
1049  }
1050 
1051  return xINJ;
1052 }
1053 
1054 std::vector<wavearray<double> > GetRecWaveform(network* NET, netevent* EVT, int ID) {
1055 
1056  std::vector<wavearray<double> > xREC; // reconstructed
1057 
1058  int nIFO = NET->ifoListSize(); // number of detectors
1059 
1060  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
1061  for(int n=0; n<nIFO; n++) {
1062 
1063  detector* pd = NET->getifo(n);
1064  int idSize = pd->RWFID.size();
1065  int wfIndex=-1;
1066  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
1067 
1068  if(wfIndex<0) {
1069  cout << "CWB_Plugin_Boostrap.C : Error : Reconstructed waveform not saved !!! : ID -> "
1070  << ID << " : detector " << NET->ifoName[n] << endl;
1071  continue;
1072  }
1073  if(wfIndex>=0) pwfREC[n] = pd->RWFP[wfIndex];
1074 
1075  wavearray<double>* wf = pwfREC[n];
1076  xREC.push_back(*wf);
1077  }
1078  delete [] pwfREC;
1079 
1080  return xREC;
1081 }
1082 
1083 std::vector<wavearray<double> > GetWaveform(network* NET, int lag, int id, char type, bool shift) {
1084 
1085  if(type!='S' && type!='s' && type!='W' && type!='w' && type!='N' && type!='n') {
1086  cout << "GetWaveform : not valid type : Abort" << endl; exit(1);
1087  }
1088 
1089  std::vector<wavearray<double> > vWAV; // reconstructed
1090 
1091  int nIFO = NET->ifoListSize(); // number of detectors
1092  netcluster* pwc = NET->getwc(lag);
1093 
1094  if(type=='S' || type=='s') {
1095  NET->getMRAwave(id,lag,'S',0,shift); // reconstruct whitened shifted pd->waveForm
1096  detector* pd;
1097  for(int i=0; i<nIFO; i++) { // loop over detectors
1098  pd = NET->getifo(i);
1099  wavearray<double>* wf = &pd->waveForm;
1100  wf->start(pwc->start+pd->waveForm.start());
1101  vWAV.push_back(*wf);
1102  }
1103  }
1104 
1105  if(type=='W' || type=='w') {
1106  NET->getMRAwave(id,lag,type,0,shift); // reconstruct+noise whitened shifted pd->waveBand
1107  detector* pd;
1108  for(int i=0; i<nIFO; i++) { // loop over detectors
1109  pd = NET->getifo(i);
1110  wavearray<double>* wf = &pd->waveBand;
1111  wf->start(pwc->start+pd->waveBand.start());
1112  vWAV.push_back(*wf);
1113  }
1114  }
1115 
1116  if(type=='N' || type=='n') {
1117  NET->getMRAwave(id,lag,'W',0,shift); // reconstruct+noise whitened shifted pd->waveBand
1118  NET->getMRAwave(id,lag,'S',0,shift); // reconstruct whitened shifted pd->waveForm
1119  detector* pd;
1120  for(int i=0; i<nIFO; i++) { // loop over detectors
1121  pd = NET->getifo(i);
1122  pd->waveNull = pd->waveBand;
1123  pd->waveNull-= pd->waveForm;
1124  wavearray<double>* wf = &pd->waveNull;
1125  wf->start(pwc->start+pd->waveNull.start());
1126  vWAV.push_back(*wf);
1127  }
1128  }
1129 
1130  return vWAV;
1131 }
1132 
1133 std::vector<wavearray<double> > GetSigWaveform(network* NET, CWB::config* cfg, int lag, int id) {
1134 // get injected waveform on TF sub-space of the reconstructed waveform vREC
1135 
1136  std::vector<wavearray<double> > vSIG; // reconstructed
1137 
1138  int nIFO = NET->ifoListSize(); // number of detectors
1139  netcluster* pwc = NET->getwc(lag);
1140 
1141  std::vector<int> pI = NET->wdmMRA.getXTalk(pwc, id, false); // buffer for pixel IDs
1142 
1143  netpixel* pix = pwc->getPixel(id,pI[0]);
1144 
1145  int tsize=1;
1146 
1147  int V = pI.size();
1148  if(!V) return vSIG;
1149  int V4 = V + (V%4 ? 4 - V%4 : 0);
1150  int V44 = V4 + 4;
1151 
1152  float En = 2*NET->acor*NET->acor*nIFO; // network energy threshold in the sky loop
1153 
1154  float *pd[NIFO], *pD[NIFO];
1155  float *v00[NIFO],*v90[NIFO];
1156  float *pa[NIFO], *pA[NIFO];
1157 
1158  float* ptmp;
1159  int i,j;
1160 
1161  std::vector<float*> _DAT; // vectors for packet amplitudes
1162  std::vector<float*> _vtd; // vectors of TD amplitudes
1163  std::vector<float*> _vTD; // vectors of TD amplitudes
1164 
1165  if(_DAT.size()) _avx_free_ps(_DAT); // container for data packet amplitudes
1166  if(_vtd.size()) _avx_free_ps(_vtd); // array for 00 amplitudes
1167  if(_vTD.size()) _avx_free_ps(_vTD); // array for 90 amplitudes
1168  for(i=0; i<NIFO; i++) {
1169  ptmp = (float*)_mm_malloc((V4*3+8)*sizeof(float),32);
1170  for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _DAT.push_back(ptmp); // concatenated arrays {amp}{AMP}{norm}{n,N,c,s}
1171  ptmp = (float*)_mm_malloc(tsize*V4*sizeof(float),32);
1172  for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vtd.push_back(ptmp); // array of aligned vectors
1173  ptmp = (float*)_mm_malloc(tsize*V4*sizeof(float),32);
1174  for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vTD.push_back(ptmp); // array of aligned vectors
1175  }
1176 
1177  std::vector<float*> _AVX; // vectors for network pixel statistics
1178  if(_AVX.size()) _avx_free_ps(_AVX);
1179  float* p_et = (float*)_mm_malloc(V4*sizeof(float),32); // 0
1180  for(j=0; j<V4; j++) p_et[j]=0; _AVX.push_back(p_et);
1181  float* pMSK = (float*)_mm_malloc(V44*sizeof(float),32); // 1 - pixel mask
1182  for(j=0; j<V44; j++) pMSK[j]=0; _AVX.push_back(pMSK); pMSK[V4]=nIFO;
1183  float* p_fp = (float*)_mm_malloc(V44*sizeof(float),32); // 2- |f+|^2 (0:V4), +norm (V4:V4+4)
1184  for(j=0; j<V44; j++) p_fp[j]=0; _AVX.push_back(p_fp);
1185  float* p_fx = (float*)_mm_malloc(V44*sizeof(float),32); // 3- |fx|^2 (0:V4), xnorm (V4:V4+4)
1186  for(j=0; j<V44; j++) p_fx[j]=0; _AVX.push_back(p_fx);
1187  float* p_si = (float*)_mm_malloc(V4*sizeof(float),32); // 4
1188  for(j=0; j<V4; j++) p_si[j]=0; _AVX.push_back(p_si);
1189  float* p_co = (float*)_mm_malloc(V4*sizeof(float),32); // 5
1190  for(j=0; j<V4; j++) p_co[j]=0; _AVX.push_back(p_co);
1191  float* p_uu = (float*)_mm_malloc((V4+16)*sizeof(float),32); // 6 - 00+ unit vector(0:V4), norm(V4), cos(V4+4)
1192  for(j=0; j<V4+16; j++) p_uu[j]=0; _AVX.push_back(p_uu);
1193  float* p_UU = (float*)_mm_malloc((V4+16)*sizeof(float),32); // 7 - 90+ unit vector(0:V4), norm(V4), sin(V4+4)
1194  for(j=0; j<V4+16; j++) p_UU[j]=0; _AVX.push_back(p_UU);
1195  float* p_vv = (float*)_mm_malloc((V4+16)*sizeof(float),32); // 8- 00x unit vector(0:V4), norm(V4), cos(V4+4)
1196  for(j=0; j<V4+16; j++) p_vv[j]=0; _AVX.push_back(p_vv);
1197  float* p_VV = (float*)_mm_malloc((V4+16)*sizeof(float),32); // 9- 90x unit vector(0:V4), norm(V4), sin(V4+4)
1198  for(j=0; j<V4+16; j++) p_VV[j]=0; _AVX.push_back(p_VV);
1199  float* p_au = (float*)_mm_malloc(V4*sizeof(float),32); // 10
1200  for(j=0; j<V4; j++) p_au[j]=0; _AVX.push_back(p_au);
1201  float* p_AU = (float*)_mm_malloc(V4*sizeof(float),32); // 11
1202  for(j=0; j<V4; j++) p_AU[j]=0; _AVX.push_back(p_AU);
1203  float* p_av = (float*)_mm_malloc(V4*sizeof(float),32); // 12
1204  for(j=0; j<V4; j++) p_av[j]=0; _AVX.push_back(p_av);
1205  float* p_AV = (float*)_mm_malloc(V4*sizeof(float),32); // 13
1206  for(j=0; j<V4; j++) p_AV[j]=0; _AVX.push_back(p_AV);
1207  float* p_uv = (float*)_mm_malloc(V4*4*sizeof(float),32); // 14 special array for GW norm calculation
1208  for(j=0; j<V4*4; j++) p_uv[j]=0; _AVX.push_back(p_uv);
1209  float* p_ee = (float*)_mm_malloc(V4*sizeof(float),32); // 15 + energy array
1210  for(j=0; j<V4; j++) p_ee[j]=0; _AVX.push_back(p_ee);
1211  float* p_EE = (float*)_mm_malloc(V4*sizeof(float),32); // 16 x energy array
1212  for(j=0; j<V4; j++) p_EE[j]=0; _AVX.push_back(p_EE);
1213  float* pTMP=(float*)_mm_malloc(V4*4*NIFO*sizeof(float),32); // 17 temporary array for _avx_norm_ps()
1214  for(j=0; j<V4*4*NIFO; j++) pTMP[j]=0; _AVX.push_back(pTMP);
1215  float* p_ni = (float*)_mm_malloc(V4*sizeof(float),32); // 18 + network index
1216  for(j=0; j<V4; j++) p_ni[j]=0; _AVX.push_back(p_ni);
1217  float* p_ec = (float*)_mm_malloc(V4*sizeof(float),32); // 19 + coherent energy
1218  for(j=0; j<V4; j++) p_ec[j]=0; _AVX.push_back(p_ec);
1219  float* p_gn = (float*)_mm_malloc(V4*sizeof(float),32); // 20 + Gaussian noise correction
1220  for(j=0; j<V4; j++) p_gn[j]=0; _AVX.push_back(p_gn);
1221  float* p_ed = (float*)_mm_malloc(V4*sizeof(float),32); // 21 + energy disbalance
1222  for(j=0; j<V4; j++) p_ed[j]=0; _AVX.push_back(p_ed);
1223  float* p_rn = (float*)_mm_malloc(V4*sizeof(float),32); // 22 + sattelite noise in TF domain
1224  for(j=0; j<V4; j++) p_rn[j]=0; _AVX.push_back(p_rn);
1225 
1226  double R = NET->getifo(0)->TFmap.rate();
1227  for(int j=0; j<V; j++) {
1228  netpixel* pix = pwc->getPixel(id,pI[j]);
1229  if(!pix->core) continue; // select only core pixels
1230  for(int n=0; n<nIFO; n++) {
1231  int index = pix->data[n].index;
1232  int ires = int(TMath::Log2(R/pix->rate))-cfg->l_low;
1233  double aa = GetSparseMapData(&jSS[n][ires], true, index);
1234  double AA = GetSparseMapData(&jSS[n][ires], false , index);
1235  //aa = GetSparseMapData(&rSS[n][ires], true, index);
1236  //AA = GetSparseMapData(&rSS[n][ires], false , index);
1237  _vtd[n][j] = aa; // copy 00 data
1238  _vTD[n][j] = AA; // copy 90 data
1239  }
1240  }
1241 
1242  for(int i=0; i<NIFO; i++) { // set up zero delay and packet pointers
1243  pd[i] = _DAT[i];
1244  pD[i] = _DAT[i]+V4;
1245  pa[i] = _vtd[i];
1246  pA[i] = _vTD[i];
1247  v00[i] = pa[i];
1248  v90[i] = pA[i];
1249  }
1250 
1251  En=0.0001;
1252  int M;
1253  double Eo=0;
1254  wavearray<float> D_snr;
1255  Eo = _avx_loadata_ps(v00,v90,pd,pD,En,_AVX,V4); // calculate data stats and store in _AVX
1256  Eo = _avx_packet_ps(pd,pD,_AVX,V4); // get data packet
1257  D_snr = NET->_avx_norm_ps(pd,pD,_AVX,V4); // data packet energy snr
1258  M = _avx_setAMP_ps(pd,pD,_AVX,V4); // set data packet amplitudes for data
1259 
1260  for(int j=0; j<V; j++) { // loop over principle components
1261  pix = pwc->getPixel(id,pI[j]);
1262  pix->likelihood = (pMSK[j]>0) ? (p_ee[j]+p_EE[j])/2 : 0; // total pixel energy
1263  if(pMSK[j]>0) pix->core = true;
1264  //pix->core = true;
1265  for(int i=0; i<nIFO; i++) {
1266  pix->setdata(double(pd[i][j]),'S',i); // 00 whitened
1267  pix->setdata(double(pD[i][j]),'P',i); // 90 whitened
1268  }
1269  }
1270  vSIG = GetWaveform(NET, lag, id, 'S', false); // get reconstructed+noise waveform
1271 
1272 #ifdef SAVE_REDUCED_INJ
1273  for(int i=0; i<nIFO; i++) {
1274  char title[256];
1275  char ofname[256];
1276  sprintf(title,"%s (TIME) : vREC(red) - vSIG(blue)",NET->ifoName[i]);
1277  sprintf(ofname,"%s_vREC_vSIG_time_id_%d.%s",NET->ifoName[i],id,"root");
1278  PlotWaveform(ofname, title, cfg,&vREC[i], &vSIG[i], NULL, NULL, false); // time
1279  sprintf(title,"%s (TIME) : vINJ(red) - vSIG(blue)",NET->ifoName[i]);
1280  sprintf(ofname,"%s_vINJ_vSIG_time_id_%d.%s",NET->ifoName[i],id,"root");
1281  PlotWaveform(ofname, title, cfg,&vINJ[i], &vSIG[i], NULL, NULL, false); // time
1282  }
1283 #endif
1284 
1285  if(_AVX.size()) _avx_free_ps(_AVX);
1286  if(_DAT.size()) _avx_free_ps(_DAT); // container for data packet amplitudes
1287  if(_vtd.size()) _avx_free_ps(_vtd); // array for 00 amplitudes
1288  if(_vTD.size()) _avx_free_ps(_vTD); // array for 90 amplitudes
1289 
1290  return vSIG;
1291 }
1292 
1294 
1295  gPE.ff[0]=0; gPE.ff[1]=0;
1296  gPE.of[0]=0; gPE.of[1]=0;
1297 
1298  int size = iSNR[0].size();
1299  if(size==0) return;
1300 
1301  for(int i=0;i<size;i++) {
1302  double isnr=0; for(int n=0;n<nIFO;n++) isnr+= iSNR[n][i];
1303  double osnr=0; for(int n=0;n<nIFO;n++) osnr+= oSNR[n][i];
1304  double iosnr=0; for(int n=0;n<nIFO;n++) iosnr+=ioSNR[n][i];
1305  double ff = iosnr/sqrt(isnr*isnr); // fitting factor
1306  double of = iosnr/sqrt(isnr*osnr); // overlap factor
1307 
1308  gPE.ff[0] += ff;
1309  gPE.ff[1] += pow(ff,2);
1310 
1311  gPE.of[0] += of;
1312  gPE.of[1] += pow(of,2);
1313  }
1314 
1315  gPE.ff[0] /= size;
1316  gPE.ff[1] = sqrt(gPE.ff[1]/size-gPE.ff[0]*gPE.ff[0]);
1317 
1318  gPE.of[0] /= size;
1319  gPE.of[1] = sqrt(gPE.of[1]/size-gPE.of[0]*gPE.of[0]);
1320 }
1321 
1322 void GetChirpMassStatistic(std::vector<double>* vCHIRP) {
1323 
1324  gPE.mch[0]=0;
1325  gPE.mch[1]=0;
1326 
1327  int size = vCHIRP[1].size();
1328  if(size==0) return;
1329 
1330  for(int i=0;i<size;i++) {
1331  gPE.mch[0] += vCHIRP[1][i];
1332  gPE.mch[1] += pow(vCHIRP[1][i],2);
1333  }
1334 
1335  gPE.mch[0] /= size;
1336  gPE.mch[1] = sqrt(gPE.mch[1]/size-gPE.mch[0]*gPE.mch[0]);
1337 }
1338 
1340 
1341  gPE.snet[0]=0;
1342  gPE.snet[1]=0;
1343 
1344  int size = iSNR[0].size();
1345  if(size==0) return;
1346 
1347  for(int i=0;i<size;i++) {
1348  double osnr=0;
1349  for(int n=0;n<nIFO;n++) osnr+=oSNR[n][i];
1350  gPE.snet[0] += sqrt(osnr);
1351  gPE.snet[1] += osnr;
1352  }
1353 
1354  gPE.snet[0] /= size;
1355  gPE.snet[1] = sqrt(gPE.snet[1]/size-gPE.snet[0]*gPE.snet[0]);
1356 }
1357 
1358 void GetNullPixels(std::vector<double>* aNUL, std::vector<double>* ANUL,
1359  network* NET, CWB::config* cfg, int lag, int id) {
1360 // get null pixels in TF domain
1361 // we use delayed vSS amplitudes and the reconstructed TF pixels save in the NET->a_00,NET->a_90 arrays
1362 
1363  size_t nIFO = NET->ifoList.size();
1364 
1365  netcluster* pwc = NET->getwc(lag);
1366 
1367  // the TD amplitudes are cleaned by the likelihoodWP function, must be computed again
1368  int sCuts = pwc->sCuts[id-1]; // save cluster status
1369  pwc->sCuts[id-1] = 0; // cluster must be enabled (mark as done by likelihoodWP)
1370  pwc->setcore(false,id);
1371  pwc->loadTDampSSE(*NET, 'a', cfg->BATCH, cfg->BATCH); // attach TD amp to pixels
1372 
1373  std::vector<int> pI = NET->wdmMRA.getXTalk(pwc, id, true); // buffer for pixel IDs
1374  int V = pI.size();
1375 
1376  netpixel* pix = pwc->getPixel(id,pI[0]);
1377  int tsize = pix->tdAmp[0].size();
1378  if(!tsize || tsize&1) { // tsize%1 = 1/0 = power/amplitude
1379  cout<<"GetNullPixels error: wrong pixel TD data\n";
1380  exit(1);
1381  }
1382  tsize /= 2;
1383  cout << "tsize :" << tsize << endl;
1384  int V4 = V + (V%4 ? 4 - V%4 : 0);
1385  std::vector<int>* vtof = &(pwc->nTofF[id-1]);
1386 
1387  for(int n=0; n<2*7*nIFO; n++) gPE.nstat[n]=0;
1388 
1389  int cnt[NIFO_MAX];
1390  for(int m=0; m<nIFO; m++) cnt[m]=0;
1391 
1392  float* a_00 = NET->a_00.data; // set pointer to 00 array of reconstructed signal
1393  float* a_90 = NET->a_90.data; // set pointer to 90 array of reconstructed signal
1394 
1395  double R = NET->getifo(0)->TFmap.rate();
1396  for(int j=0; j<V; j++) {
1397  netpixel* pix = pwc->getPixel(id,pI[j]);
1398  //if(!pix->core) continue; // select only core pixels
1399  float ee=0;
1400  for(int m=0; m<nIFO; m++) {
1401  double ss = a_00[j*NIFO+m];
1402  double SS = a_90[j*NIFO+m];
1403  ee += pow(ss,2)+pow(SS,2);
1404  }
1405  if(ee==0) continue; // only core pixels are selected
1406  for(int m=0; m<nIFO; m++) {
1407  int index = pix->data[m].index;
1408  int ires = int(TMath::Log2(R/pix->rate))-cfg->l_low;
1409  // get original sparse map amplitudes
1410  //double dd = GetSparseMapData(&vSS[m][ires], true, index);
1411  //double DD = GetSparseMapData(&vSS[m][ires], false , index);
1412  //double dd = GetSparseMapData(&dSS[m][ires], true, index);
1413  //double DD = GetSparseMapData(&dSS[m][ires], false , index);
1414 
1415  // get shifted data (signal+noise) amplitudes
1416  int l = (*vtof)[m]+tsize/2;
1417  double dd = pix->tdAmp[m].data[l]; // copy TD 00 data
1418  double DD = pix->tdAmp[m].data[l+tsize]; // copy TD 90 data
1419 
1420  //double ss = GetSparseMapData(&rSS[m][ires], true, index);
1421  //double SS = GetSparseMapData(&rSS[m][ires], false , index);
1422  // get reconstructed amplitudes (shifted)
1423  double ss = a_00[j*NIFO+m];
1424  double SS = a_90[j*NIFO+m];
1425 
1426  // 00 phase
1427  aNUL[m].push_back(dd-ss);
1428 
1429  // 90 phase
1430  ANUL[m].push_back(DD-SS);
1431 
1432  cnt[m]++;
1433  }
1434  }
1435 
1436  pwc->sCuts[id-1] = sCuts; // restore cluster status
1437  pwc->clean(id); // cean TD ammplitudes
1438 
1439  return;
1440 }
1441 
1443 
1444  int layer = SS->GetLayer(index);
1445  int start = SS->sparseLookup[layer]; // sparse table layer offset
1446  int end = SS->sparseLookup[layer+1]-1; // sparse table layer+1 offset
1447 
1448  int i = SS->binarySearch(SS->sparseIndex.data,start,end,index);
1449 
1450  if(i<0) return 0;
1451 
1452  return phase ? SS->sparseMap00[i] : SS->sparseMap90[i];
1453 }
1454 
1457  bool fft, TString pdir, double P, EColor col1, EColor col2, EColor col3) {
1458 
1459  watplot PTS(const_cast<char*>("ptspe"),200,20,800,500);
1460 
1461  //cout << "Print " << ofname << endl;
1462  double tmin=1.e20;
1463  if(wref==NULL) return;
1464  if(wf1==NULL) return;
1465  else tmin=wf1->start();
1466  if(wf2!=NULL) if(wf2->start()<tmin) tmin=wf2->start();
1467  if(wf3!=NULL) if(wf3->start()<tmin) tmin=wf3->start();
1468 
1469  tmin=gSEGGPS;
1470 
1471  wf1->start(wf1->start()-tmin);
1472  if(wf2!=NULL) if(wf2!=wf1) wf2->start(wf2->start()-tmin);
1473  if(wf3!=NULL) if(wf3!=wf1 && wf3!=wf2) wf3->start(wf3->start()-tmin);
1474  if(wref!=wf1 && wref!=wf2 && wref!=wf3) wref->start(wref->start()-tmin);
1475 
1476  double bT, eT;
1477  GetTimeBoundaries(*wref, P, bT, eT);
1478 
1479  if(fft) {
1480  wavearray<double> wf=*wf1; int k=0;
1481  for(int i=0;i<wf1->size();i++) {
1482  double time=i/wf1->rate()+wf1->start();
1483  if(time>bT && time<eT) wf[k++]=wf1->data[i];
1484  }
1485  wf.resize(k);
1486  PTS.plot(&wf, const_cast<char*>("AL"), col1, 0, 0, true, cfg->fLow, cfg->fHigh);
1487  } else {
1488  PTS.plot(wf1, const_cast<char*>("AL"), col1, bT, eT);
1489  PTS.graph[0]->GetXaxis()->SetRangeUser(bT, eT);
1490  }
1491  PTS.graph[0]->SetLineWidth(1);
1492  PTS.graph[0]->SetTitle(title);
1493 
1494  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",tmin);
1495  PTS.graph[0]->GetXaxis()->SetTitle(xtitle);
1496 
1497  if(wf2!=NULL) {
1498  if(fft) {
1499  if(wf2!=NULL) {
1500  wavearray<double> wf=*wf2; int k=0;
1501  for(int i=0;i<wf2->size();i++) {
1502  double time=i/wf2->rate()+wf2->start();
1503  if(time>bT && time<eT) wf[k++]=wf2->data[i];
1504  }
1505  wf.resize(k);
1506  PTS.plot(&wf, const_cast<char*>("SAME"), col2, 0, 0, true, cfg->fLow, cfg->fHigh);
1507  }
1508  } else {
1509  if(wf2!=NULL) PTS.plot(wf2, const_cast<char*>("SAME"), col2, 0, 0);
1510  }
1511  if(wf2!=NULL) PTS.graph[1]->SetLineWidth(2);
1512  }
1513 
1514  if(wf3!=NULL) {
1515  if(fft) {
1516  if(wf3!=NULL) {
1517  wavearray<double> wf=*wf3; int k=0;
1518  for(int i=0;i<wf3->size();i++) {
1519  double time=i/wf3->rate()+wf3->start();
1520  if(time>bT && time<eT) wf[k++]=wf3->data[i];
1521  }
1522  wf.resize(k);
1523  PTS.plot(&wf, const_cast<char*>("SAME"), col3, 0, 0, true, cfg->fLow, cfg->fHigh);
1524  }
1525  } else {
1526  if(wf3!=NULL) PTS.plot(wf3, const_cast<char*>("SAME"), col3, 0, 0);
1527  }
1528  if(wf3!=NULL) PTS.graph[2]->SetLineWidth(1);
1529  }
1530 
1531  wf1->start(wf1->start()+tmin);
1532  if(wf2!=NULL) if(wf2!=wf1) wf2->start(wf2->start()+tmin);
1533  if(wf3!=NULL) if(wf3!=wf1 && wf3!=wf2) wf3->start(wf3->start()+tmin);
1534  if(wref!=wf1 && wref!=wf2 && wref!=wf3) wref->start(wref->start()+tmin);
1535 
1536  if(fft) {
1537  PTS.canvas->SetLogx();
1538  PTS.canvas->SetLogy();
1539  }
1540 
1541  if(ofname!="") {
1542  ofname = TString(pdir)+TString("/")+ofname;
1543  PTS.canvas->Print(ofname);
1544  cout << "write : " << ofname << endl;
1545  }
1546 }
1547 
1549 
1550  double R=wf1->rate();
1551 
1552  double b_wf1 = wf1->start();
1553  double e_wf1 = wf1->start()+wf1->size()/R;
1554  double b_wf2 = wf2->start();
1555  double e_wf2 = wf2->start()+wf2->size()/R;
1556 
1557  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
1558  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
1559 
1560  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
1561  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
1562  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1563 
1564  double wf11=0.;
1565  double wf22=0.;
1566  double wf12=0.;
1567  for(int i=0;i<sizeXCOR;i++) wf12 += wf1->data[i+o_wf1] * wf2->data[i+o_wf2];
1568  for(int i=0;i<wf1->size();i++) wf11 += pow(wf1->data[i],2);
1569  for(int i=0;i<wf2->size();i++) wf22 += pow(wf2->data[i],2);
1570 
1571  double FF = wf12/sqrt(wf11*wf22);
1572 
1573  return FF;
1574 }
1575 
1577 
1578  double R=wf1->rate();
1579 
1580  double b_wf1 = wf1->start();
1581  double e_wf1 = wf1->start()+wf1->size()/R;
1582  double b_wf2 = wf2->start();
1583  double e_wf2 = wf2->start()+wf2->size()/R;
1584 
1585  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
1586  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
1587 
1588  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
1589  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
1590  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1591 
1592  wavearray<double> wfdif(sizeXCOR);
1593  wfdif=0.;
1594  wfdif.rate(R);
1595  wfdif.start(b_wf1+double(o_wf1)/R);
1596 
1597  for(int i=0;i<sizeXCOR;i++) wfdif[i] = wf1->data[i+o_wf1] - wf2->data[i+o_wf2];
1598 
1599  return wfdif;
1600 }
1601 
1603 
1604  detector* pD = NET->getifo(ifoID);
1605 
1606  int nRES = cfg->l_high-cfg->l_low+1; // number of frequency resolution levels
1607  for(int i=0;i<nRES;i++) {
1608  // rebuild wseries from sparse table only with core pixels
1609  bool core = true;
1610  SSeries<double> ss = pD->vSS[i];
1611  ss.Expand(core);
1612  ss.Inverse();
1613 
1614  char title[256];
1615  char ofname[256];
1616  sprintf(title,"%s (TIME) : INJ(red) - SM(blue)",NET->ifoName[ifoID]);
1617  sprintf(ofname,"%s_INJ_SM_time_id_%d_res_%d.%s",NET->ifoName[ifoID],ID,i,PLOT_TYPE);
1618  PlotWaveform(ofname, title, cfg, wf, &ss, NULL, NULL, false); // time
1619  }
1620 }
1621 
1623 
1624  int n_amp_cal_err=0;
1625  int n_phs_cal_err=0;
1626  if(options.CompareTo("")!=0) {
1627  cout << options << endl;
1628  if(!options.Contains("--")) { // parameters are used only by cwb_inet
1629 
1630  TObjArray* token = TString(options).Tokenize(TString(' '));
1631  for(int j=0;j<token->GetEntries();j++){
1632 
1633  TObjString* tok = (TObjString*)token->At(j);
1634  TString stok = tok->GetString();
1635 
1636  if(stok.Contains("pe_trials=") || stok.Contains("pe_retry=")) {
1637  TString pe_trials=stok;
1638  pe_trials.Remove(0,pe_trials.Last('=')+1);
1639  if(pe_trials.IsDigit()) gOPT.trials=pe_trials.Atoi();
1640  if(gOPT.trials>MAX_TRIALS) {
1641  cout << "ReadUserOptions Error : trials must be <= MAX_TRIALS : " << MAX_TRIALS << endl;exit(1);
1642  }
1643  }
1644 
1645  if(stok.Contains("pe_ced_dump=")) {
1646  TString pe_ced_dump=stok;
1647  pe_ced_dump.Remove(0,pe_ced_dump.Last('=')+1);
1648  if(pe_ced_dump.IsDigit()) gOPT.ced_dump=pe_ced_dump.Atoi();
1649  }
1650 
1651  if(stok.Contains("pe_ced_enable=")) {
1652  TString pe_ced_enable=stok;
1653  pe_ced_enable.Remove(0,pe_ced_enable.Last('=')+1);
1654  if(pe_ced_enable=="tfmap") gOPT.ced_tfmap = true;
1655  if(pe_ced_enable=="rdr") gOPT.ced_rdr = true;
1656  if(pe_ced_enable=="psd") gOPT.ced_psd = true;
1657  if(pe_ced_enable=="skymap") gOPT.ced_skymap = true;
1658  if(pe_ced_enable=="rec") gOPT.ced_rec = true;
1659  if(pe_ced_enable=="inj") gOPT.ced_inj = true;
1660  if(pe_ced_enable=="rinj") gOPT.ced_rinj = true;
1661  if(pe_ced_enable=="cm") gOPT.ced_cm = true;
1662  if(pe_ced_enable=="distr") gOPT.ced_distr = true;
1663  if(pe_ced_enable=="null") gOPT.ced_null = true;
1664  if(pe_ced_enable=="pca") gOPT.ced_pca = true;
1665  }
1666 
1667  if(stok.Contains("pe_ced_disable=")) {
1668  TString pe_ced_disable=stok;
1669  pe_ced_disable.Remove(0,pe_ced_disable.Last('=')+1);
1670  if(pe_ced_disable=="tfmap") gOPT.ced_tfmap = false;
1671  if(pe_ced_disable=="rdr") gOPT.ced_rdr = false;
1672  if(pe_ced_disable=="psd") gOPT.ced_psd = false;
1673  if(pe_ced_disable=="skymap") gOPT.ced_skymap = false;
1674  if(pe_ced_disable=="rec") gOPT.ced_rec = false;
1675  if(pe_ced_disable=="inj") gOPT.ced_inj = false;
1676  if(pe_ced_disable=="rinj") gOPT.ced_rinj = false;
1677  if(pe_ced_disable=="cm") gOPT.ced_cm = false;
1678  if(pe_ced_disable=="distr") gOPT.ced_distr = false;
1679  if(pe_ced_disable=="null") gOPT.ced_null = false;
1680  if(pe_ced_disable=="pca") gOPT.ced_pca = false;
1681  }
1682 
1683  if(stok.Contains("pe_seed=")) {
1684  TString pe_seed=stok;
1685  pe_seed.Remove(0,pe_seed.Last('=')+1);
1686  if(pe_seed.IsDigit()) gOPT.seed=pe_seed.Atoi();
1687  }
1688 
1689  if(stok.Contains("pe_gps=")) {
1690  TString pe_gps=stok;
1691  pe_gps.Remove(0,pe_gps.Last('=')+1);
1692  if(pe_gps.IsFloat()) gOPT.gps=pe_gps.Atof();
1693  }
1694 
1695  if(stok.Contains("pe_amp_cal_err=")) {
1696  TString pe_amp_cal_err=stok;
1697  pe_amp_cal_err.Remove(0,pe_amp_cal_err.Last('=')+1);
1698  if(pe_amp_cal_err.IsFloat()) gOPT.amp_cal_err[n_amp_cal_err]=pe_amp_cal_err.Atof();
1699  if(n_amp_cal_err<(NIFO_MAX-1)) n_amp_cal_err++;
1700  }
1701 
1702  if(stok.Contains("pe_phs_cal_err=")) {
1703  TString pe_phs_cal_err=stok;
1704  pe_phs_cal_err.Remove(0,pe_phs_cal_err.Last('=')+1);
1705  if(pe_phs_cal_err.IsFloat()) gOPT.phs_cal_err[n_phs_cal_err]=pe_phs_cal_err.Atof();
1706  if(n_phs_cal_err<(NIFO_MAX-1)) n_phs_cal_err++;
1707  }
1708 
1709  if(stok.Contains("pe_skymask=")) {
1710  TString pe_skymask=stok;
1711  pe_skymask.Remove(0,pe_skymask.Last('=')+1);
1712  if(pe_skymask=="false") gOPT.skymask=0;
1713  if(pe_skymask=="true") gOPT.skymask=1;
1714  if(pe_skymask.IsDigit()) gOPT.skymask=pe_skymask.Atoi();
1715  }
1716 
1717  if(stok.Contains("pe_signal=")) {
1718  TString pe_signal=stok;
1719  pe_signal.Remove(0,pe_signal.Last('=')+1);
1720  if(pe_signal.IsDigit()) gOPT.signal=pe_signal.Atoi();
1721  }
1722 
1723  if(stok.Contains("pe_noise=")) {
1724  TString pe_noise=stok;
1725  pe_noise.Remove(0,pe_noise.Last('=')+1);
1726  if(pe_noise=="false") gOPT.noise=0;
1727  if(pe_noise=="true") gOPT.noise=1;
1728  if(pe_noise.IsDigit()) gOPT.noise=pe_noise.Atoi();
1729  }
1730 
1731  if(stok.Contains("pe_multitask=")) {
1732  TString pe_multitask=stok;
1733  pe_multitask.Remove(0,pe_multitask.Last('=')+1);
1734  if(pe_multitask=="true") gOPT.multitask=true;
1735  if(pe_multitask=="false") gOPT.multitask=false;
1736  }
1737 
1738  if(stok.Contains("pe_output_enable=")) {
1739  TString pe_output_enable=stok;
1740  pe_output_enable.Remove(0,pe_output_enable.Last('=')+1);
1741  if(pe_output_enable=="inj") gOPT.output_inj=true;
1742  if(pe_output_enable=="rec") gOPT.output_rec=true;
1743  if(pe_output_enable=="wht") gOPT.output_wht=true;
1744  if(pe_output_enable=="dat") gOPT.output_dat=true;
1745  if(pe_output_enable=="med") gOPT.output_med=true;
1746  if(pe_output_enable=="p50") gOPT.output_p50=true;
1747  if(pe_output_enable=="p90") gOPT.output_p90=true;
1748  if(pe_output_enable=="avr") gOPT.output_avr=true;
1749  if(pe_output_enable=="rms") gOPT.output_rms=true;
1750  }
1751 
1752  if(stok.Contains("pe_output_disable=")) {
1753  TString pe_output_disable=stok;
1754  pe_output_disable.Remove(0,pe_output_disable.Last('=')+1);
1755  if(pe_output_disable=="inj") gOPT.output_inj=false;
1756  if(pe_output_disable=="rec") gOPT.output_rec=false;
1757  if(pe_output_disable=="wht") gOPT.output_wht=false;
1758  if(pe_output_disable=="dat") gOPT.output_dat=false;
1759  if(pe_output_disable=="med") gOPT.output_med=false;
1760  if(pe_output_disable=="p50") gOPT.output_p50=false;
1761  if(pe_output_disable=="p90") gOPT.output_p90=false;
1762  if(pe_output_disable=="avr") gOPT.output_avr=false;
1763  if(pe_output_disable=="rms") gOPT.output_rms=false;
1764  }
1765 
1766  }
1767  }
1768  }
1769 }
1770 
1772 
1773  cout << "-----------------------------------------" << endl;
1774  cout << "PE config options " << endl;
1775  cout << "-----------------------------------------" << endl << endl;
1776  cout << "PE_TRIALS " << gOPT.trials << endl;
1777  cout << "PE_SIGNAL " << gOPT.signal << endl;
1778  cout << "PE_NOISE " << gOPT.noise << endl;
1779  for(int n=0;n<cfg->nIFO;n++) {
1780  cout << "PE_AMP_CAL_ERR " << cfg->ifo[n] << " " << gOPT.amp_cal_err[n] << endl;
1781  cout << "PE_PHS_CAL_ERR " << cfg->ifo[n] << " " << gOPT.phs_cal_err[n] << endl;
1782  }
1783  cout << "PE_SKYMASK " << gOPT.skymask << endl;
1784  cout << "PE_SEED " << gOPT.seed << endl;
1785  cout.precision(14);
1786  cout << "PE_GPS " << gOPT.gps << endl;
1787  cout.precision(6);
1788 
1789  cout << "PE_MULTITASK " << gOPT.multitask << endl;
1790 
1791  cout << "PE_CED_DUMP " << gOPT.ced_dump << endl;
1792 
1793  cout << "PE_CED_TFMAP " << gOPT.ced_tfmap << endl;
1794  cout << "PE_CED_RDR " << gOPT.ced_rdr << endl;
1795  cout << "PE_CED_PSD " << gOPT.ced_psd << endl;
1796  cout << "PE_CED_SKYMAP " << gOPT.ced_skymap << endl;
1797  cout << "PE_CED_REC " << gOPT.ced_rec << endl;
1798  cout << "PE_CED_INJ " << gOPT.ced_inj << endl;
1799  cout << "PE_CED_rINJ " << gOPT.ced_rinj << endl;
1800  cout << "PE_CED_CM " << gOPT.ced_cm << endl;
1801  cout << "PE_CED_DISTR " << gOPT.ced_distr << endl;
1802  cout << "PE_CED_NULL " << gOPT.ced_null << endl;
1803  cout << "PE_CED_PCA " << gOPT.ced_pca << endl;
1804 
1805  cout << "PE_OUTPUT_INJ " << gOPT.output_inj << endl;
1806  cout << "PE_OUTPUT_REC " << gOPT.output_rec << endl;
1807  cout << "PE_OUTPUT_WHT " << gOPT.output_wht << endl;
1808  cout << "PE_OUTPUT_DAT " << gOPT.output_dat << endl;
1809  cout << "PE_OUTPUT_MED " << gOPT.output_med << endl;
1810  cout << "PE_OUTPUT_P50 " << gOPT.output_p50 << endl;
1811  cout << "PE_OUTPUT_P90 " << gOPT.output_p90 << endl;
1812  cout << "PE_OUTPUT_AVR " << gOPT.output_avr << endl;
1813  cout << "PE_OUTPUT_RMS " << gOPT.output_rms << endl;
1814 
1815  cout << endl;
1816 }
1817 
1819 
1820  TString ofName = odir+"/pe_config.txt";
1821 
1822  ofstream out;
1823  out.open(ofName,ios::out);
1824  if(!out.good()) {cout << "DumpUserOptions : Error Opening File : " << ofName << endl;exit(1);}
1825 
1826  TString info="";
1827  TString tabs="\t\t\t\t";
1828 
1829  char version[128];
1830  sprintf(version,"WAT Version : %s - SVN Revision : %s - Tag/Branch : %s",watversion('f'),watversion('r'),watversion('b'));
1831 
1832  out << endl;
1833  out << "--------------------------------" << endl;
1834  out << "PE version " << endl;
1835  out << "--------------------------------" << endl;
1836  out << endl;
1837  out << version << endl;
1838  out << endl;
1839 
1840  out << "--------------------------------" << endl;
1841  out << "PE config options " << endl;
1842  out << "--------------------------------" << endl;
1843  out << endl;
1844 
1845  out << "pe_trials " << gOPT.trials << endl;
1846  info = "// number of trials (def=100)";
1847  out << tabs << info << endl;
1848 
1849  out << "pe_signal " << gOPT.signal << endl;
1850  info = "// signal used for trials : 0(def)/1/2";
1851  out << tabs << info << endl;
1852  info = "// 0 - reconstructed waveform";
1853  out << tabs << info << endl;
1854  info = "// 1 - injected waveform (available only in simulation mode)";
1855  out << tabs << info << endl;
1856  info = "// 2 - original waveform = (reconstructed+null)";
1857  out << tabs << info << endl;
1858 
1859  out << "pe_noise " << gOPT.noise << endl;
1860  info = "// noise used for trials : 0/1(def)/2 ";
1861  out << tabs << info << endl;
1862  info = "// 0 - waveform is injected in the whitened HoT and apply Coherent+SuperCluster+Likelihood stages";
1863  out << tabs << info << endl;
1864  info = "// 1 - add gaussian noise to likelihood sparse maps and apply Likelihood stage";
1865  out << tabs << info << endl;
1866  info = "// 2 - add whitened HoT to likelihood sparse maps and apply Likelihood stage";
1867  out << tabs << info << endl;
1868 
1869  for(int n=0;n<cfg->nIFO;n++) {
1870  out << "pe_amp_cal_err " << cfg->ifo[n] << " " << gOPT.amp_cal_err[n] << endl;
1871  }
1872  info = "// max percentage of amplitude miscalibration : def(0) -> disabled";
1873  out << tabs << info << endl;
1874  info = "// if>0 -> uniform in [1-pe_amp_cal_err,1+pe_amp_cal_err] else gaus(1,pe_amp_cal_err)";
1875  out << tabs << info << endl;
1876 
1877  for(int n=0;n<cfg->nIFO;n++) {
1878  out << "pe_phs_cal_err " << cfg->ifo[n] << " " << gOPT.phs_cal_err[n] << endl;
1879  }
1880  info = "// max phase (degrees) miscalibration : def(0) -> disabled";
1881  out << tabs << info << endl;
1882  info = "// if>0 -> uniform in [-pe_phs_cal_err,+pe_phs_cal_err] else gaus(pe_phs_cal_err)";
1883  out << tabs << info << endl;
1884 
1885  out << "pe_multitask " << (gOPT.multitask ? "enable" : "disable") << endl;
1886  info = "// if enabled, PE trials are executed by different jobs in multitask : only simulation=0 and single event (gps>0)";
1887  out << tabs << info << endl;
1888 
1889  out << "pe_ced_dump " << gOPT.ced_dump << endl;
1890  info = "// dump CED at gLRETRY=PE_CED_DUMP (def(-1) -> disabled)";
1891  out << tabs << info << endl;
1892 
1893  out << "pe_skymask " << gOPT.skymask << endl;
1894  info = "// skymask used for trials : 0(def)/1/2/3 ";
1895  out << tabs << info << endl;
1896  info = "// 0 - disable -> search over all sky";
1897  out << tabs << info << endl;
1898  out << tabs << "// 1 - skymask with radius " << SKYMASK_RADIUS << "(deg) and source selected according the reconstructed skymap probability " << endl;
1899  info = "// 2 - skymask select only the sky position where the waveform is reconstructed (the maximum of detection statistic)";
1900  out << tabs << info << endl;
1901  info = "// 3 - skymask select only the sky position where the waveform has been injected";
1902  out << tabs << info << endl;
1903 
1904  out << "pe_seed " << gOPT.seed << endl;
1905  info = "// seed used by PE for random generation - 0(def) -> random seed";
1906  out << tabs << info << endl;
1907 
1908  out.precision(14);
1909  out << "pe_gps " << gOPT.gps << endl;
1910  info = "// if >0 only gps +/- iwindow is analyzed - 0(def) -> disabled";
1911  out << tabs << info << endl;
1912  out.precision(6);
1913 
1914  out << endl;
1915  out << "pe_ced tfmap " << (gOPT.ced_tfmap ? "enable" : "disable") << endl;
1916  info = "// tfmap - pe_ced_enable(def)/disable : Shows/Hides the Time-Frequency Maps Tab -> Spectrograms, Scalograms, TF Likelihood,Null";
1917  out << tabs << info << endl;
1918  out << "pe_ced rdr " << (gOPT.ced_rdr ? "enable" : "disable") << endl;
1919  info = "// rdr - pe_ced_enable(def)/disable : Shows/Hides Reconstructed Detector Response Tab";
1920  out << tabs << info << endl;
1921  out << "pe_ced psd " << (gOPT.ced_psd ? "enable" : "disable") << endl;
1922  info = "// psd - pe_ced_enable(def)/disable : Shows/Hides Power Spectral Density Tab";
1923  out << tabs << info << endl;
1924  out << "pe_ced skymap " << (gOPT.ced_skymap ? "enable" : "disable") << endl;
1925  info = "// skymap - pe_ced_enable(def)/disable : Shows/Hides SkyMaps Tab";
1926  out << tabs << info << endl;
1927  out << "pe_ced rec " << (gOPT.ced_rec ? "enable" : "disable") << endl;
1928  info = "// rec - pe_ced_enable(def)/disable : Shows/Hides Reconstructed Waveforms/Instantaneous-Frequency with errors Tab";
1929  out << tabs << info << endl;
1930  out << "pe_ced inj " << (gOPT.ced_inj ? "enable" : "disable") << endl;
1931  info = "// inj - pe_ced_enable(def)/disable : Showsi/Hides Injection' tab reports the comparison with the injected whitened waveform";
1932  out << tabs << info << endl;
1933  out << "pe_ced rinj " << (gOPT.ced_rinj ? "enable" : "disable") << endl;
1934  info = "// rinj - pe_ced_enable/disable(def) : Shows/Hides Injection' tab reports the comparison with the injected whitened waveform";
1935  out << tabs << info << endl;
1936  info = " in time-frequency subspace of the PointEstimated waveform";
1937  out << tabs << info << endl;
1938  out << "pe_ced cm " << (gOPT.ced_cm ? "enable" : "disable") << endl;
1939  info = "// cm - pe_ced_enable(def)/disable : Shows/Hides Chirp Mass Value/Error Distributions Tab";
1940  out << tabs << info << endl;
1941  out << "pe_ced distr " << (gOPT.ced_distr ? "enable" : "disable") << endl;
1942  info = "// distr - pe_ced_enable/disable(def) : Shows/Hides Residuals Distributions Tab";
1943  out << tabs << info << endl;
1944  out << "pe_ced null " << (gOPT.ced_null ? "enable" : "disable") << endl;
1945  info = "// null - pe_ced_enable/disable(def) : Shows/Hides Null Pixels Distributions Tab";
1946  out << tabs << info << endl;
1947  out << "pe_ced pca " << (gOPT.ced_pca ? "enable" : "disable") << endl;
1948  info = "// pca - pe_ced_enable/disable(def) : Shows/Hides PCA likelihood TF plot in the Time-Frequency Maps Tab";
1949  out << tabs << info << endl;
1950 
1951  out.close();
1952 }
1953 
1955 
1956  gOPT.trials = PE_TRIALS;
1957  gOPT.signal = PE_SIGNAL;
1958  gOPT.noise = PE_NOISE;
1959  for(int n=0;n<NIFO_MAX;n++) gOPT.amp_cal_err[n] = PE_AMP_CAL_ERR;
1960  for(int n=0;n<NIFO_MAX;n++) gOPT.phs_cal_err[n] = PE_PHS_CAL_ERR;
1961  gOPT.skymask = PE_SKYMASK;
1962  gOPT.seed = PE_SEED;
1963  gOPT.gps = PE_GPS;
1964 
1965  gOPT.multitask = PE_MULTITASK;
1966 
1967  gOPT.ced_dump = PE_CED_DUMP;
1968 
1969  gOPT.ced_tfmap = PE_CED_TFMAP;
1970  gOPT.ced_rdr = PE_CED_RDR;
1971  gOPT.ced_psd = PE_CED_PSD;
1972  gOPT.ced_skymap = PE_CED_SKYMAP;
1973  gOPT.ced_rec = PE_CED_REC;
1974  gOPT.ced_inj = PE_CED_INJ;
1975  gOPT.ced_rinj = PE_CED_rINJ;
1976  gOPT.ced_cm = PE_CED_CM;
1977  gOPT.ced_distr = PE_CED_DISTR;
1978  gOPT.ced_null = PE_CED_NULL;
1979  gOPT.ced_pca = PE_CED_PCA;
1980 
1981  gOPT.output_inj = PE_OUTPUT_INJ;
1982  gOPT.output_rec = PE_OUTPUT_REC;
1983  gOPT.output_wht = PE_OUTPUT_WHT;
1984  gOPT.output_dat = PE_OUTPUT_DAT;
1985  gOPT.output_med = PE_OUTPUT_MED;
1986  gOPT.output_p50 = PE_OUTPUT_P50;
1987  gOPT.output_p90 = PE_OUTPUT_P90;
1988  gOPT.output_avr = PE_OUTPUT_AVR;
1989  gOPT.output_rms = PE_OUTPUT_RMS;
1990 }
1991 
1992 void SetOutputFile(network* NET, netevent* &EVT, CWB::config* cfg, bool dump_pe) {
1993 
1994  // import slagShift
1995  float* gSLAGSHIFT=NULL; IMPORT(float*,gSLAGSHIFT)
1996 
1997  int nIFO = NET->ifoListSize(); // number of detectors
1998 
1999  // search output root file in the system list
2000  TFile* froot = NULL;
2001  TList *files = (TList*)gROOT->GetListOfFiles();
2002  gOUTPUT="";
2003  if (files) {
2004  TIter next(files);
2005  TSystemFile *file;
2006  TString fname;
2007  bool check=false;
2008  while ((file=(TSystemFile*)next())) {
2009  fname = file->GetName();
2010  // set output root file as the current file
2011  if(fname.Contains("wave_")) {
2012  froot=(TFile*)file;froot->cd();
2013  gOUTPUT=fname;
2014  gOUTPUT.ReplaceAll(".root.tmp",".txt");
2015  //cout << "output file name : " << fname << endl;
2016  }
2017  }
2018  if(!froot) {
2019  cout << "CWB_Plugin_Boostrap.C : Error - output root file not found" << endl;
2020  gSystem->Exit(1);
2021  }
2022  } else {
2023  cout << "CWB_Plugin_Boostrap.C : Error - output root file not found" << endl;
2024  gSystem->Exit(1);
2025  }
2026 
2027  if(dump_pe) {
2028  if(gOPT.output_inj) if(cfg->simulation) for(int n=0;n<nIFO;n++) gPE.wINJ[n] = &vINJ[n];
2029  if(gOPT.output_rec) for(int n=0;n<nIFO;n++) gPE.wREC[n] = &vREC[n];
2030  if(gOPT.output_wht) for(int n=0;n<nIFO;n++) gPE.wWHT[n] = &vWHT[n];
2031  if(gOPT.output_dat) for(int n=0;n<nIFO;n++) gPE.wDAT[n] = &vDAT[n];
2032  if(gOPT.output_med) for(int n=0;n<nIFO;n++) gPE.wMED[n] = &vMED[n];
2033  if(gOPT.output_p50) for(int n=0;n<nIFO;n++) gPE.wL50[n] = &vL50[n];
2034  if(gOPT.output_p50) for(int n=0;n<nIFO;n++) gPE.wU50[n] = &vU50[n];
2035  if(gOPT.output_p90) for(int n=0;n<nIFO;n++) gPE.wL90[n] = &vL90[n];
2036  if(gOPT.output_p90) for(int n=0;n<nIFO;n++) gPE.wU90[n] = &vU90[n];
2037  if(gOPT.output_avr) for(int n=0;n<nIFO;n++) gPE.wAVR[n] = &vAVR[n];
2038  if(gOPT.output_rms) for(int n=0;n<nIFO;n++) gPE.wRMS[n] = &vRMS[n];
2039  }
2040 
2041  gTREE = (TTree *) froot->Get("waveburst");
2042  if(gTREE!=NULL) {
2043  EVT = new netevent(gTREE,nIFO);
2044  if(dump_pe) {
2045  gTREE->Branch("pe_trials",&gPE.trials,"pe_trials/I");
2046  gTREE->Branch("pe_erR",gPE.erR,TString::Format("pe_erR[%i]/F",11));
2047  gTREE->Branch("pe_erF",gPE.erF,TString::Format("pe_erF[%i]/F",11));
2048  gTREE->Branch("pe_nstat",gPE.nstat,TString::Format("pe_nstat[%i]/F",2*7*cfg->nIFO));
2049  gTREE->Branch("pe_snet",gPE.snet,TString::Format("pe_snet[%i]/F",2));
2050  gTREE->Branch("pe_ff",gPE.ff,TString::Format("pe_ff[%i]/F",2));
2051  gTREE->Branch("pe_of",gPE.of,TString::Format("pe_of[%i]/F",2));
2052  gTREE->Branch("pe_mch",gPE.mch,TString::Format("pe_mch[%i]/F",2));
2053 
2054  for(int n=0;n<nIFO;n++) {
2055  if(gOPT.output_inj) if(cfg->simulation) gTREE->Branch(TString::Format("pe_wINJ_%d",n).Data(),"wavearray<double>",&gPE.wINJ[n],32000,0);
2056  if(gOPT.output_rec) gTREE->Branch(TString::Format("pe_wREC_%d",n).Data(),"wavearray<double>",&gPE.wREC[n],32000,0);
2057  if(gOPT.output_wht) gTREE->Branch(TString::Format("pe_wWHT_%d",n).Data(),"wavearray<double>",&gPE.wWHT[n],32000,0);
2058  if(gOPT.output_dat) gTREE->Branch(TString::Format("pe_wDAT_%d",n).Data(),"wavearray<double>",&gPE.wDAT[n],32000,0);
2059  if(gOPT.output_med) gTREE->Branch(TString::Format("pe_wMED_%d",n).Data(),"wavearray<double>",&gPE.wMED[n],32000,0);
2060  if(gOPT.output_p50) gTREE->Branch(TString::Format("pe_wL50_%d",n).Data(),"wavearray<double>",&gPE.wL50[n],32000,0);
2061  if(gOPT.output_p50) gTREE->Branch(TString::Format("pe_wU50_%d",n).Data(),"wavearray<double>",&gPE.wU50[n],32000,0);
2062  if(gOPT.output_p90) gTREE->Branch(TString::Format("pe_wL90_%d",n).Data(),"wavearray<double>",&gPE.wL90[n],32000,0);
2063  if(gOPT.output_p90) gTREE->Branch(TString::Format("pe_wU90_%d",n).Data(),"wavearray<double>",&gPE.wU90[n],32000,0);
2064  if(gOPT.output_avr) gTREE->Branch(TString::Format("pe_wAVR_%d",n).Data(),"wavearray<double>",&gPE.wAVR[n],32000,0);
2065  if(gOPT.output_rms) gTREE->Branch(TString::Format("pe_wRMS_%d",n).Data(),"wavearray<double>",&gPE.wRMS[n],32000,0);
2066  }
2067 
2068  if(gOPT.multitask&&(gMTRIAL!=gOPT.trials)) { // save reconstricted waveforms and skyprob when multitask mode
2069  for(int n=0;n<nIFO;n++) gTREE->Branch(TString::Format("mtpe_wREC_%d",n).Data(),"wavearray<double>",&wREC[0][n],32000,0);
2070  gTREE->Branch("mtpe_skyprob","skymap",&wSKYPROB[0], 32000,0);
2071  }
2072  }
2073  } else {
2074  EVT = new netevent(nIFO);
2075  gTREE = EVT->setTree();
2076  }
2077  EVT->setSLags(gSLAGSHIFT); // set slags into netevent
2078  EVT->Psave=cfg->Psave;
2079 }
2080 
2081 void DumpOutputFile(network* NET, netevent* &EVT, CWB::config* cfg, int ID, int k, int factor) {
2082 
2083  if(cfg->dump) EVT->dopen(gOUTPUT.Data(),const_cast<char*>("a"),false);
2084  EVT->output2G(gTREE,NET,ID,k,factor); // get reconstructed parameters
2085  if(cfg->dump) EVT->dclose();
2086  if(!cfg->cedDump) NET->getwc(k)->sCuts[ID-1]=1; // mark as processed
2087 }
2088 
2090 
2091  while(!vINJ.empty()) vINJ.pop_back();
2092  vINJ.clear(); std::vector<wavearray<double> >().swap(vINJ);
2093  while(!vREC.empty()) vREC.pop_back();
2094  vREC.clear(); std::vector<wavearray<double> >().swap(vREC);
2095  while(!vWHT.empty()) vWHT.pop_back();
2096  vWHT.clear(); std::vector<wavearray<double> >().swap(vWHT);
2097  while(!vPCA.empty()) vPCA.pop_back();
2098  vPCA.clear(); std::vector<wavearray<double> >().swap(vPCA);
2099  while(!vDAT.empty()) vDAT.pop_back();
2100  vDAT.clear(); std::vector<wavearray<double> >().swap(vDAT);
2101  while(!vNUL.empty()) vNUL.pop_back();
2102  vNUL.clear(); std::vector<wavearray<double> >().swap(vNUL);
2103  while(!cINJ.empty()) cINJ.pop_back();
2104  cINJ.clear(); std::vector<wavearray<double> >().swap(cINJ);
2105 
2106  while(!vAVR.empty()) vAVR.pop_back();
2107  vAVR.clear(); std::vector<wavearray<double> >().swap(vAVR);
2108  while(!vRMS.empty()) vRMS.pop_back();
2109  vRMS.clear(); std::vector<wavearray<double> >().swap(vRMS);
2110 
2111  while(!vMED.empty()) vMED.pop_back();
2112  vMED.clear(); std::vector<wavearray<double> >().swap(vMED);
2113  while(!vL50.empty()) vL50.pop_back();
2114  vL50.clear(); std::vector<wavearray<double> >().swap(vL50);
2115  while(!vU50.empty()) vU50.pop_back();
2116  vU50.clear(); std::vector<wavearray<double> >().swap(vU50);
2117  while(!vL90.empty()) vL90.pop_back();
2118  vL90.clear(); std::vector<wavearray<double> >().swap(vL90);
2119  while(!vU90.empty()) vU90.pop_back();
2120  vU90.clear(); std::vector<wavearray<double> >().swap(vU90);
2121 
2122  while(!fREC.empty()) fREC.pop_back();
2123  fREC.clear(); std::vector<wavearray<double> >().swap(fREC);
2124  while(!fINJ.empty()) fINJ.pop_back();
2125  fINJ.clear(); std::vector<wavearray<double> >().swap(fINJ);
2126  while(!fAVR.empty()) fAVR.pop_back();
2127  fAVR.clear(); std::vector<wavearray<double> >().swap(fAVR);
2128  while(!fRMS.empty()) fRMS.pop_back();
2129  fRMS.clear(); std::vector<wavearray<double> >().swap(fRMS);
2130 
2131  while(!fMED.empty()) fMED.pop_back();
2132  fMED.clear(); std::vector<wavearray<double> >().swap(fMED);
2133  while(!fL50.empty()) fL50.pop_back();
2134  fL50.clear(); std::vector<wavearray<double> >().swap(fL50);
2135  while(!fU50.empty()) fU50.pop_back();
2136  fU50.clear(); std::vector<wavearray<double> >().swap(fU50);
2137  while(!fL90.empty()) fL90.pop_back();
2138  fL90.clear(); std::vector<wavearray<double> >().swap(fL90);
2139  while(!fU90.empty()) fU90.pop_back();
2140  fU90.clear(); std::vector<wavearray<double> >().swap(fU90);
2141 
2142  while(!vRES.empty()) vRES.pop_back();
2143  vRES.clear(); std::vector<double >().swap(vRES);
2144  while(!fRES.empty()) fRES.pop_back();
2145  fRES.clear(); std::vector<double >().swap(fRES);
2146 
2147  for(int n=0;n<MAX_TRIALS;n++) {
2148  while(!wREC[n].empty()) wREC[n].pop_back();
2149  wREC[n].clear(); std::vector<wavearray<double> >().swap(wREC[n]);
2150  }
2151 
2152  for(int n=0;n<NIFO_MAX;n++) {
2153  while(!iSNR[n].empty()) iSNR[n].pop_back();
2154  iSNR[n].clear(); std::vector<double>().swap(iSNR[n]);
2155  while(!oSNR[n].empty()) oSNR[n].pop_back();
2156  oSNR[n].clear(); std::vector<double>().swap(oSNR[n]);
2157  while(!ioSNR[n].empty()) ioSNR[n].pop_back();
2158  ioSNR[n].clear(); std::vector<double>().swap(ioSNR[n]);
2159  }
2160 
2161  for(int n=0;n<6;n++) {
2162  while(!vCHIRP[n].empty()) vCHIRP[n].pop_back();
2163  vCHIRP[n].clear(); std::vector<double>().swap(vCHIRP[n]);
2164  }
2165 
2166  while(!vLIKELIHOOD.empty()) vLIKELIHOOD.pop_back();
2167  vLIKELIHOOD.clear(); std::vector<double >().swap(vLIKELIHOOD);
2168 
2169  for(int n=0;n<NIFO_MAX;n++) {
2170  while(!aNUL[n].empty()) aNUL[n].pop_back();
2171  aNUL[n].clear(); std::vector<double>().swap(aNUL[n]);
2172  while(!ANUL[n].empty()) ANUL[n].pop_back();
2173  ANUL[n].clear(); std::vector<double>().swap(ANUL[n]);
2174  }
2175 
2176  for(int n=0;n<NIFO_MAX;n++) {
2177  while(!aNSE[n].empty()) aNSE[n].pop_back();
2178  aNSE[n].clear(); std::vector<double>().swap(aNSE[n]);
2179  while(!ANSE[n].empty()) ANSE[n].pop_back();
2180  ANSE[n].clear(); std::vector<double>().swap(ANSE[n]);
2181  }
2182 
2183  for(int n=0;n<NIFO_MAX;n++) {
2184  while(!vSS[n].empty()) vSS[n].pop_back();
2185  vSS[n].clear(); std::vector<SSeries<double> >().swap(vSS[n]);
2186  while(!rSS[n].empty()) rSS[n].pop_back();
2187  rSS[n].clear(); std::vector<SSeries<double> >().swap(rSS[n]);
2188  while(!jSS[n].empty()) jSS[n].pop_back();
2189  jSS[n].clear(); std::vector<SSeries<double> >().swap(jSS[n]);
2190  while(!dSS[n].empty()) dSS[n].pop_back();
2191  dSS[n].clear(); std::vector<SSeries<double> >().swap(dSS[n]);
2192  }
2193 }
2194 
2195 void PlotWaveforms(network* NET, CWB::config* cfg, int ID, TString pdir) {
2196 
2197  int nIFO = NET->ifoListSize(); // number of detectors
2198 
2199  wavearray<double> wDIF,wALI,wENV;
2200 
2201  char title[256]; char ofname[256];
2202  for(int n=0; n<nIFO; n++) {
2203 
2204  // COMPUTE SN, TString pdirR
2205 /*
2206  double SNR=0;
2207  for(int j=0;j<vREC[n].size();j++) SNR+=pow(vREC[n].data[j],2);
2208  cout << "SNR " << NET->ifoName[n] << " " << sqrt(SNR) << endl;
2209  double NUL=0;
2210  for(int j=0;j<vNUL[n].size();j++) NUL+=pow(vNUL[n].data[j],2);
2211  cout << "NUL " << NET->ifoName[n] << " " << sqrt(NUL) << endl;
2212 */
2213 /*
2214  // PLOT -> SPARSE MAP
2215  PlotSparse(n, NET, cfg, ID, &vINJ[n]);
2216 */
2217 /*
2218  // PLOT -> INJ : PCA : REC
2219  sprintf(title,"%s (TIME) : vINJ(red) - vPCA(blue) - vREC(green)",NET->ifoName[n]);
2220  sprintf(ofname,"%s_vINJ_vPCA_vREC_time_id_%d.%s",NET->ifoName[n],ID,PLOT_TYPE);
2221  PlotWaveform(ofname, title, cfg, &vINJ[n], &vPCA[n], &vREC[n],NULL, false,pdir); // time
2222  sprintf(title,"%s (FFT) : vINJ(red) - vPCA(blue) - vREC(green)",NET->ifoName[n]);
2223  //sprintf(ofname,"%s_vINJ_vPCA_vREC_fft_id_%d.%s",NET->ifoName[n],ID,PLOT_TYPE);
2224  //PlotWaveform(ofname, title, cfg, &vINJ[n], &vPCA[n], &vREC[n],NULL, true,pdir); // fft
2225 */
2226 
2227 #ifdef SET_WAVEFORM_CUTS
2228  // PLOT -> cINJ : cINJ-vARV : vRMS
2229  wDIF = GetDifWaveform(&cINJ[n], &vAVR[n]);
2230  sprintf(title,"%s (TIME) : cINJ(red) - cINJ-vAVR(blue) - vRMS(green)",NET->ifoName[n]);
2231  sprintf(ofname,"%s_cINJ_cINJvAVR_vRMS_time.%s",NET->ifoName[n],PLOT_TYPE);
2232  PlotWaveform(ofname, title, cfg, &cINJ[n], &wDIF, &vRMS[n], NULL, false, pdir);
2233 #endif
2234 /*
2235  // PLOT -> cINJ : wREC
2236  sprintf(title,"%s (TIME) : cINJ(red) - wREC(blue)",NET->ifoName[n]);
2237  sprintf(ofname,"%s_cINJ_wREC_time.%s",NET->ifoName[n],PLOT_TYPE);
2238  PlotWaveform(ofname, title, cfg, &cINJ[n], &wREC[0][n], NULL, NULL, false, pdir);
2239 */
2240 
2241 #ifdef SET_WAVEFORM_CUTS
2242  // PLOT -> cINJ : vREC : vAVR
2243  sprintf(title,"%s (TIME) : cINJ(red) - vREC(blue) - vAVR(green)",NET->ifoName[n]);
2244  sprintf(ofname,"%s_cINJ_vREC_vAVR_time.%s",NET->ifoName[n],PLOT_TYPE);
2245  PlotWaveform(ofname, title, cfg, &cINJ[n], &vREC[n], &vAVR[n], NULL, false, pdir);
2246 #endif
2247 
2248 #ifdef SET_WAVEFORM_CUTS
2249  // PLOT -> cINJ : cINJ
2250  sprintf(title,"%s (TIME) : cINJ(red) - cINJ(blue)",NET->ifoName[n]);
2251  sprintf(ofname,"%s_cINJ_cINJ_time.%s",NET->ifoName[n],PLOT_TYPE);
2252  PlotWaveform(ofname, title, cfg, &cINJ[n], &cINJ[n], NULL, NULL, false, pdir);
2253 #endif
2254 
2255  if(gOPT.ced_rec) {
2256  // PLOT -> vWHT vs vREC
2257  // NOTE : this plot overwrite the standard IFO_wf_signal.png/IFO_wf_signal_fft.png/ plots
2258  // the band signal is the now the whitend data
2259  double scale = sqrt(1<<cfg->levelR);
2260  vWHT[n]*=1/scale;
2261  vREC[n]*=1/scale;
2262  sprintf(ofname,"%s_wf_signal.%s",NET->ifoName[n],PLOT_TYPE);
2263  PlotWaveform(ofname, "", cfg, &vWHT[n], &vREC[n], NULL, &vREC[n], false, pdir,0.999,(EColor)16,kRed);
2264  sprintf(ofname,"%s_wf_signal_fft.%s",NET->ifoName[n],PLOT_TYPE);
2265  PlotWaveform(ofname, "", cfg, &vWHT[n], &vREC[n], NULL, &vREC[n], true, pdir,0.999,(EColor)16,kRed);
2266  vWHT[n]*=scale;
2267  vREC[n]*=scale;
2268  }
2269 
2270  for(int j=0;j<2;j++) {
2271 
2272  TString tag = j==0 ? "time" : "time_zoom"; // file tag
2273  double P = j==0 ? 0.995 : 0.9; // time zoom
2274 
2275  // RECONSTRUCTED
2276 
2277  if(gOPT.ced_rec) {
2278 
2279  // PLOT -> vAVR(vMED) : vRMS(PERC)
2280  //wENV = GetWaveformEnvelope(&vREC[n]);
2281  sprintf(title,"%s (TIME) : vREC (red) - vRMS:2vRMS (grey:light_gray)",NET->ifoName[n]);
2282  sprintf(ofname,"%s_vREC_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2283 #ifdef PLOT_MEDIAN
2284  PlotWaveformAsymmErrors(ofname, "", cfg, &vREC[n], &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2285 #else
2286  PlotWaveformErrors(ofname, "", cfg, &vREC[n], &vAVR[n], &vRMS[n], &vREC[n], pdir, P);
2287 #endif
2288 
2289  // PLOT -> vREC-vAVR(vMED) : vRMS(PERC)
2290 #ifdef PLOT_MEDIAN
2291  wDIF = GetDifWaveform(&vREC[n], &vMED[n]);
2292 #else
2293  wDIF = GetDifWaveform(&vREC[n], &vAVR[n]);
2294 #endif
2295  sprintf(title,"%s (TIME) : vREC-vAVR(red) - 2*vRMS(gray)",NET->ifoName[n]);
2296  sprintf(ofname,"%s_vRECvAVR_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2297 #ifdef PLOT_MEDIAN
2298  PlotWaveformAsymmErrors(ofname, "", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2299 #else
2300  PlotWaveformErrors(ofname, "", cfg, &wDIF, &wDIF, &vRMS[n], &vREC[n], pdir, P);
2301 #endif
2302 
2303  // PLOT -> fAVR : fRMS(fMED)
2304  sprintf(title,"%s (FREQ) : fREC (red) - fRMS:2fRMS (grey:light_gray)",NET->ifoName[n]);
2305  sprintf(ofname,"%s_fREC_fRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2306 #ifdef PLOT_MEDIAN
2307  PlotWaveformAsymmErrors(ofname, "", cfg, &fREC[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P, true);
2308 #else
2309  PlotFrequencyErrors(ofname, "", cfg, &fREC[n], &fAVR[n], &fRMS[n], &vREC[n], pdir, P);
2310 #endif
2311 
2312  // PLOT -> vREC : vDAT-vREC : vRMS
2313  wDIF = GetDifWaveform(&vDAT[n], &vREC[n]);
2314  sprintf(title,"%s (TIME) : vREC(red) - vDAT-vREC(blue) - vRMS(green)",NET->ifoName[n]);
2315  sprintf(ofname,"%s_vREC_vDATvREC_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2316  PlotWaveform(ofname, "", cfg, &vREC[n], &wDIF, &vRMS[n], &vREC[n], false, pdir, P);
2317  }
2318 
2319  if(!cfg->simulation) continue;
2320 
2321  // RECONSTRUCTED vs INJECTED
2322 
2323  if(gOPT.ced_inj) {
2324  // PLOT -> vINJ : vAVR(vMED) : vRMS(PERC)
2325  wALI = GetAlignedWaveform(&vINJ[n], &vREC[n]);
2326  sprintf(title,"%s (TIME) : vINJ (red) - vRMS:2vRMS (grey:light_gray)",NET->ifoName[n]);
2327  sprintf(ofname,"%s_vINJ_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2328 #ifdef PLOT_MEDIAN
2329  PlotWaveformAsymmErrors(ofname, "", cfg, &wALI, &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2330 #else
2331  PlotWaveformErrors(ofname, "", cfg, &wALI, &vAVR[n], &vRMS[n], &vREC[n], pdir, P);
2332 #endif
2333 
2334  // PLOT -> fINJ : fAVE(fMED) : fRMS(PERC)
2335  sprintf(title,"%s (FREQ) : fINJ (red) - fRMS:2fRMS (grey:light_gray)",NET->ifoName[n]);
2336  sprintf(ofname,"%s_fINJ_fRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2337 #ifdef PLOT_MEDIAN
2338  PlotWaveformAsymmErrors(ofname, "", cfg, &fINJ[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P, true);
2339 #else
2340  PlotFrequencyErrors(ofname, "", cfg, &fINJ[n], &fAVR[n], &fRMS[n], &vREC[n], pdir, P);
2341 #endif
2342 
2343  // PLOT -> vAVR(vMED)-vINJ : vRMS(PERC)
2344 #ifdef PLOT_MEDIAN
2345  wDIF = GetDifWaveform(&vMED[n], &vINJ[n]);
2346 #else
2347  wDIF = GetDifWaveform(&vAVR[n], &vINJ[n]);
2348 #endif
2349  wDIF = GetAlignedWaveform(&wDIF, &vREC[n]);
2350  sprintf(title,"%s (TIME) : vAVR-vINJ(red) - 2*vRMS(gray)",NET->ifoName[n]);
2351  sprintf(ofname,"%s_vAVRvINJ_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2352 #ifdef PLOT_MEDIAN
2353  PlotWaveformAsymmErrors(ofname, "", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2354 #else
2355  PlotWaveformErrors(ofname, "", cfg, &wDIF, &wDIF, &vRMS[n], &vREC[n], pdir, P);
2356 #endif
2357  }
2358 
2359  // RECONSTRUCTED vs reduced INJECTED
2360 
2361  if(gOPT.ced_rinj) {
2362 
2363  // PLOT -> vAVR(vMED)-cINJ : vRMS(PERC)
2364 #ifdef PLOT_MEDIAN
2365  wDIF = GetDifWaveform(&vMED[n], &cINJ[n]);
2366 #else
2367  wDIF = GetDifWaveform(&vAVR[n], &cINJ[n]);
2368 #endif
2369  wDIF = GetAlignedWaveform(&wDIF, &vREC[n]);
2370  sprintf(title,"%s (TIME) : vAVR-cINJ(red) - 2*vRMS(gray)",NET->ifoName[n]);
2371  sprintf(ofname,"%s_vAVRcINJ_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2372 #ifdef PLOT_MEDIAN
2373  PlotWaveformAsymmErrors(ofname, "", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2374 #else
2375  PlotWaveformErrors(ofname, "", cfg, &wDIF, &wDIF, &vRMS[n], &vREC[n], pdir, P);
2376 #endif
2377 
2378  // PLOT -> cINJ : vAVR(vMED) : vRMS(PERC)
2379  sprintf(title,"%s (TIME) : cINJ (red) - vRMS:2vRMS (grey:light_gray)",NET->ifoName[n]);
2380  sprintf(ofname,"%s_cINJ_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2381 #ifdef PLOT_MEDIAN
2382  PlotWaveformAsymmErrors(ofname, "", cfg, &cINJ[n], &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2383 #else
2384  PlotWaveformErrors(ofname, "", cfg, &cINJ[n], &vAVR[n], &vRMS[n], &vREC[n], pdir, P);
2385 #endif
2386 
2387  // PLOT -> fINJ : fAVR(PERC) : fRMS
2388  sprintf(title,"%s (FREQ) : fINJ (red) - fRMS:2fRMS (grey:light_gray)",NET->ifoName[n]);
2389  sprintf(ofname,"%s_fINJ_fRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2390 #ifdef PLOT_MEDIAN
2391  PlotWaveformAsymmErrors(ofname, "", cfg, &fINJ[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P, true);
2392 #else
2393  PlotFrequencyErrors(ofname, "", cfg, &fINJ[n], &fAVR[n], &fRMS[n], &vREC[n], pdir, P);
2394 #endif
2395 
2396  // PLOT -> cINJ : vREC : vRMS
2397  sprintf(title,"%s (TIME) : cINJ(red) - vREC(blue) - vRMS(green)",NET->ifoName[n]);
2398  sprintf(ofname,"%s_cINJ_vREC_vRMS_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2399  PlotWaveform(ofname, "", cfg, &cINJ[n], &vREC[n], &vRMS[n], &vREC[n], false, pdir, P);
2400 
2401  // PLOT -> vINJ : vINJ-cINJ : vINJcINJ
2402  wDIF = GetDifWaveform(&vINJ[n], &cINJ[n]);
2403  wDIF = GetAlignedWaveform(&wDIF, &vREC[n]);
2404  wALI = GetAlignedWaveform(&vINJ[n], &vREC[n]);
2405  sprintf(title,"%s (TIME) : vINJ(red) - cINJ(blue) - vINJ-cINJ(green)",NET->ifoName[n]);
2406  sprintf(ofname,"%s_vINJ_cINJ_vINJcINJ_%s.%s",NET->ifoName[n],tag.Data(),PLOT_TYPE);
2407  PlotWaveform(ofname, "", cfg, &wALI, &cINJ[n], &wDIF, &vREC[n], false, pdir, P);
2408  }
2409  }
2410 /*
2411  // FITTING FACTOR
2412  double FF1 = FittingFactor(&cINJ[n], &vREC[n]);
2413  double FF2 = FittingFactor(&cINJ[n], &vPCA[n]);
2414  cout << NET->ifoName[n] << " -> Fitting Factor - inj/rec = " << FF1 << " inj/pca = " << FF2 << endl;
2415 */
2416  }
2417 }
2418 
2419 void SetSkyMask(network* NET, CWB::config* cfg, double theta, double phi, double radius) {
2420 
2421  // --------------------------------------------------------
2422  // define SetSkyMask
2423  // --------------------------------------------------------
2424  sprintf(cfg->skyMaskFile,"--theta %f --phi %f --radius %f",90-theta,phi,radius);
2425  cout << endl << "SetSkyMask : " << cfg->skyMaskFile << endl << endl;
2426  gCWB.SetSkyMask(NET,cfg,cfg->skyMaskFile,'e');
2427 
2428 }
2429 
2431 
2432  if(type!='j' && type!='r' && type!='v') {
2433  cout << "AddNoiseAndCalErrToSparse - Error : wrong wave type" << endl; exit(1);
2434  }
2435 
2436  double d2r = TMath::Pi()/180.;
2437 
2438  int m;
2439  int nIFO = NET->ifoListSize(); // number of detectors
2440  int nRES = NET->wdmListSize(); // number of frequency resolution levels
2441 
2442  // update vSS maps
2445  SSeries<double>* p;
2446  WDM<double>* pwdm = NULL;
2447 
2448  for(int n=0; n<nIFO; n++) { // create random time series
2449  p = NET->getifo(n)->getSTFmap(0);
2450  WF[n].resize(p->wdm_nSTS);
2451  WF[n].rate(p->wdm_rate);
2452  WF[n].start(p->wdm_start);
2453 
2454  int size = WF[n].size();
2455  double rate = WF[n].rate();
2456  int jS,jE,jW,k;
2457  double dt,df,ts;
2458  TComplex C;
2459  switch(gOPT.noise) {
2460  case 1: // gaussian noise
2461  for(int i=0; i<WF[n].size(); i++) {
2462  WF[n].data[i] = gRandom->Gaus(0,NOISE_SIGMA);
2463  }
2464  cout << "AddNoiseAndCalErrToSparse - Info : add gaussian noise to waveform type : " << type << endl;
2465  break;
2466  case 2: // data noise
2467  w = *(NET->getifo(n)->getTFmap()); // get whitened data
2468  w.Inverse(); // get level 0
2469  // apply time shift to input whitened data (integer number of sammples)
2470  jS = cfg->segEdge*rate;
2471  jE = size-jS;
2472  jW = rate*cfg->iwindow/2;
2473  k = gRandom->Uniform(jS+jW,jE-jW); // select random offset (excludes +/- iwindow/2 around the event)
2474  WF[n] = w;
2475  for(int i=jS; i<jE; i++) {
2476  WF[n].data[k++] = w.data[i];
2477  if(k==jE) k=jS;
2478  }
2479  // apply time shift to input whitened data (fraction of dt)
2480  WF[n].FFTW(1);
2481  dt = 1./rate;
2482  df = rate/size;
2483  ts = gRandom->Uniform(0.,dt); // select a random time shift within the sample time range
2484  for (int j=0;j<size/2;j++) {
2485  TComplex X(WF[n].data[2*j],WF[n].data[2*j+1]);
2486  X=X*C.Exp(TComplex(0.,-2*PI*j*df*ts));
2487  WF[n].data[2*j]=X.Re();
2488  WF[n].data[2*j+1]=X.Im();
2489  }
2490  WF[n].FFTW(-1);
2491  printf("Info : %s add data noise with time shift : %.5f sec to waveform type : %c\n",NET->getifo(n)->Name,(double)k/rate+ts,type);
2492  break;
2493  default:
2494  WF[n]=0;
2495  break;
2496  }
2497  }
2498 
2499  // fill waveform sparse maps
2500  for(int n=0; n<nIFO; n++) {
2501  double A = 1;
2502  if(gOPT.amp_cal_err[n]) {
2503  double amp_cal_err = 1;
2504  if(gOPT.amp_cal_err[n]>0) amp_cal_err = gRandom->Uniform(1-gOPT.amp_cal_err[n],1+gOPT.amp_cal_err[n]);
2505  else amp_cal_err = gRandom->Gaus(1,fabs(gOPT.amp_cal_err[n]));
2506  cout << NET->ifoName[n] << " -> amp_cal_err : " << amp_cal_err << endl;
2507  A = amp_cal_err;
2508  }
2509  double C=1,S=0;
2510  if(gOPT.phs_cal_err[n]) {
2511  double phs_cal_err = 0;
2512  if(gOPT.phs_cal_err[n]>0) phs_cal_err = gRandom->Uniform(-gOPT.phs_cal_err[n],gOPT.phs_cal_err[n]);
2513  else phs_cal_err = gRandom->Gaus(0,fabs(gOPT.phs_cal_err[n]));
2514  cout << NET->ifoName[n] << " -> phs_cal_err : " << phs_cal_err << endl;
2515  C = cos(-phs_cal_err*d2r);
2516  S = sin(-phs_cal_err*d2r);
2517  }
2518  for(int j=0; j<nRES; j++) {
2519  w.Forward(WF[n],*(NET->wdmList[j]));
2520  int k = NET->getifo(n)->getSTFind(w.wrate()); // pointer to sparse TF array
2521  p = NET->getifo(n)->getSTFmap(k); // pointer to sparse TF array
2522  for(int l=0; l<p->sparseMap00.size(); l++) {
2523  m = p->sparseIndex[l];
2524  double aa,AA;
2525  if(type=='r') { // reconstructed waveform
2526  aa = rSS[n][k].sparseMap00[l];
2527  AA = rSS[n][k].sparseMap90[l];
2528  }
2529  if(type=='j') { // injected waveform
2530  aa = jSS[n][k].sparseMap00[l];
2531  AA = jSS[n][k].sparseMap90[l];
2532  }
2533  if(type=='v') { // original sparse map (reconstructed+noise)
2534  aa = vSS[n][k].sparseMap00[l];
2535  AA = vSS[n][k].sparseMap90[l];
2536  }
2537  // add calibration error
2538  aa*=A; AA*=A;
2539  // add phase error
2540  double bb = C*aa + S*AA;
2541  double BB = -S*aa + C*AA ;
2542  // add gaussian noise
2543  p->sparseMap00[l]=bb+w.data[m];
2544  p->sparseMap90[l]=BB+w.data[m+w.maxIndex()+1];
2545  }
2546  }
2547  }
2548  return;
2549 }
2550 
2551 void Wave2Sparse(network* NET, CWB::config* cfg, char type) {
2552 
2553  if(type!='j' && type!='r' && type!='v' && type!='d') {
2554  cout << "Wave2Sparse - Error : wrong wave type" << endl; exit(1);
2555  }
2556 
2557  // init waveform sparse maps
2558  int nIFO = NET->ifoListSize(); // number of detectors
2559  for(int n=0;n<nIFO;n++) {
2560  detector* pD = NET->getifo(n);
2561  if(type=='v') vSS[n]=pD->vSS;
2562  if(type=='r') rSS[n]=pD->vSS;
2563  if(type=='j') jSS[n]=pD->vSS;
2564  if(type=='d') dSS[n]=pD->vSS;
2565  }
2566  if(type=='v') return;
2567 
2568  int nRES = cfg->l_high-cfg->l_low+1; // number of frequency resolution levels
2569 
2570  // build waveform array
2572  for(int n=0; n<nIFO; n++) {
2573  WF[n].resize(vSS[n][0].wdm_nSTS);
2574  WF[n].rate(vSS[n][0].wdm_rate);
2575  WF[n].start(vSS[n][0].wdm_start);
2576  WF[n]=0;
2577  if(type=='j') {
2578  int wos = double(vINJ[n].start()-WF[n].start())*WF[n].rate();
2579  for (int i=0;i<vINJ[n].size();i++) WF[n][wos+i] = vINJ[n][i];
2580  }
2581  if(type=='r') {
2582  int wos = double(vREC[n].start()-WF[n].start())*WF[n].rate();
2583  for (int i=0;i<vREC[n].size();i++) WF[n][wos+i] = vREC[n][i];
2584  }
2585  if(type=='d') {
2586  int wos = double(vDAT[n].start()-WF[n].start())*WF[n].rate();
2587  for (int i=0;i<vDAT[n].size();i++) WF[n][wos+i] = vDAT[n][i];
2588  }
2589 
2590 #ifdef SAVE_WAVE2SPARSE_PLOT
2591  gwavearray<double> gw(&WF[n]);
2592  gw.Draw(GWAT_TIME);
2593  watplot* plot = gw.GetWATPLOT();
2594  TString gfile;
2595  if(type=='r') gfile=TString::Format("%s/WAVE2SPARSE_REC_%s.root",".",NET->ifoName[n]);
2596  if(type=='j') gfile=TString::Format("%s/WAVE2SPARSE_INJ_%s.root",".",NET->ifoName[n]);
2597  if(type=='d') gfile=TString::Format("%s/WAVE2SPARSE_DAT_%s.root",".",NET->ifoName[n]);
2598  (*plot) >> gfile;
2599 #endif
2600  }
2601 
2602  // fill waveform sparse maps
2603  WSeries<double> w;
2604  WDM<double>* pwdm = NULL;
2605  double r00,r90;
2606  for(int n=0; n<nIFO; n++) {
2607  for(int i=0; i<nRES; i++) {
2608  w.Forward(WF[n],*(NET->wdmList[i]));
2609  int k = NET->getifo(n)->getSTFind(w.wrate()); // pointer to sparse TF array
2610 
2611  int size;
2612  if(type=='r') size = rSS[n][k].sparseMap00.size();
2613  if(type=='j') size = jSS[n][k].sparseMap00.size();
2614  if(type=='d') size = dSS[n][k].sparseMap00.size();
2615 
2616  for(int j=0; j<size; j++) {
2617  int index;
2618  if(type=='r') index = rSS[n][k].sparseIndex[j];
2619  if(type=='j') index = jSS[n][k].sparseIndex[j];
2620  if(type=='d') index = dSS[n][k].sparseIndex[j];
2621  double v00 = w.pWavelet->pWWS[index];
2622  double v90 = w.pWavelet->pWWS[index+w.maxIndex()+1];
2623  if(type=='r') {
2624  rSS[n][k].sparseMap00[j]=v00;
2625  rSS[n][k].sparseMap90[j]=v90;
2626  }
2627  if(type=='j') {
2628  jSS[n][k].sparseMap00[j]=v00;
2629  jSS[n][k].sparseMap90[j]=v90;
2630  }
2631  if(type=='d') {
2632  dSS[n][k].sparseMap00[j]=v00;
2633  dSS[n][k].sparseMap90[j]=v90;
2634  }
2635  }
2636  }
2637  }
2638 }
2639 
2640 std::vector<int> ComputeStatisticalErrors(network* NET, CWB::config* cfg, int ID) {
2641 
2642  int nIFO = NET->ifoListSize(); // number of detectors
2643  std::vector<int> nrec(nIFO); // number of averaged
2644 
2645 #ifdef SET_WAVEFORM_CUTS
2646  // get time and frequency ranges
2647  gBT=1e20; gET=0;
2648  gBF=1e20; gEF=0;
2649  for(int n=0;n<nIFO;n++) {
2650  double bT,eT,bF,eF;
2651  GetTimeBoundaries(vREC[n], EFRACTION_CUT, bT, eT);
2652  GetFrequencyBoundaries(vREC[n], EFRACTION_CUT, bF, eF);
2653  if(bT<gBT) gBT=bT; if(eT>gET) gET=eT;
2654  if(bF<gBF) gBF=bF; if(eF>gEF) gEF=eF;
2655  }
2656  cout << endl;
2657  cout << "CUT TIME RANGE : " << gBT-gBT << " " << gET-gBT << endl;
2658  cout << "CUT FREQ RANGE : " << gBF << " " << gEF << endl;
2659  cout << endl;
2660 #endif
2661 
2662  // copy cINJ (only the vREC time range) into winj
2663  wavearray<double> winj[NIFO_MAX];
2665  if(cfg->simulation) {
2666  for(int n=0;n<nIFO;n++) {
2667  winj[n] = GetAlignedWaveform(&cINJ[n], &vREC[n]);
2668  finj[n] = CWB::Toolbox::getHilbertIFrequency(winj[n]);
2669 #ifdef SET_WAVEFORM_CUTS
2670  SetWaveformCuts(&winj[n], gBT, gET, gBF, gEF);
2671  cINJ.push_back(GetAlignedWaveform(&winj[n], &vREC[n]));
2672 #endif
2673  }
2674  }
2675 
2676  // compute average of wREC (only the vREC time range) into wavr
2679  for(int n=0;n<nIFO;n++) {
2680  nrec[n]=0;
2681  for(int i=0;i<gOPT.trials;i++) {
2682  if(wREC[i].size()) { // event detected
2683  wrec[i][n] = GetAlignedWaveform(&wREC[i][n], &vREC[n]);
2684  frec[i][n] = CWB::Toolbox::getHilbertIFrequency(wrec[i][n]);
2685 #ifdef SET_WAVEFORM_CUTS
2686  SetWaveformCuts(&wrec[i][n], gBT, gET, gBF, gEF);
2687 #endif
2688  nrec[n]++;
2689  } else { // event rejected (set to 0)
2690  wrec[i][n] = vREC[n]; wrec[i][n]=0;
2691  frec[i][n] = wrec[i][n];
2692  }
2693  }
2694  if(nrec[n]==0) {nrec.clear();return nrec;} // return if no events are detects
2695  }
2696 
2697  gPE.trials = nrec[0]; // save number of effective trials
2698  cout << endl << "ComputeStatisticalErrors - detected/trials : " << nrec[0] << "/" << gOPT.trials << endl << endl;
2699 
2700  // compute vAVR & vRMS
2701  wavearray<double> wavr[NIFO_MAX];
2702  wavearray<double> wrms[NIFO_MAX];
2703  for(int n=0;n<nIFO;n++) {
2704  wrms[n] = wrec[0][n]; wrms[n] = 0;
2705  wavr[n] = wrec[0][n]; wavr[n] = 0;
2706  for(int i=0;i<gOPT.trials;i++) {
2707  wavr[n] += wrec[i][n];
2708  for(int j=0;j< wrec[i][n].size();j++) wrms[n][j] += wrec[i][n][j]*wrec[i][n][j];
2709  }
2710  wavr[n] *= 1./nrec[n];
2711 
2712  wrms[n] *= 1./nrec[n];
2713  for(int j=0;j< wrms[n].size();j++) {
2714  wrms[n][j] -= wavr[n][j]*wavr[n][j];
2715  wrms[n][j] = sqrt(wrms[n][j]);
2716  }
2717 
2718  vAVR.push_back(wavr[n]);
2719  vRMS.push_back(wrms[n]);
2720  }
2721 
2722  // compute vMED, vL50, vU50, vL90, vU90
2723  wavearray<double> wmed[NIFO_MAX];
2724  wavearray<double> wl50[NIFO_MAX];
2725  wavearray<double> wu50[NIFO_MAX];
2726  wavearray<double> wl90[NIFO_MAX];
2727  wavearray<double> wu90[NIFO_MAX];
2728  for(int n=0;n<nIFO;n++) {
2729  wmed[n] = wrec[0][n]; wmed[n] = 0;
2730  wl50[n] = wrec[0][n]; wl50[n] = 0;
2731  wu50[n] = wrec[0][n]; wu50[n] = 0;
2732  wl90[n] = wrec[0][n]; wl90[n] = 0;
2733  wu90[n] = wrec[0][n]; wu90[n] = 0;
2734 
2735  int ntry = gPE.trials; // number of detected events in the trials
2736  int *index = new int[ntry];
2737  float *value = new float[ntry];
2738  for(int j=0;j<wrec[0][n].size();j++) {
2739 
2740  int k=0; for(int i=0;i<gOPT.trials;i++) if(wREC[i].size()) value[k++] = wrec[i][n][j]; // select detected events
2741  TMath::Sort(ntry,value,index,false);
2742 
2743  int imed = (ntry*50.)/100.; if(imed>=ntry) imed=ntry-1;
2744  wmed[n][j] = value[index[imed]];
2745 
2746  int il50 = (ntry*25.)/100.; if(il50>=ntry) il50=ntry-1;
2747  int iu50 = (ntry*75.)/100.; if(iu50>=ntry) iu50=ntry-1;
2748  int il90 = (ntry*5.)/100.; if(il90>=ntry) il90=ntry-1;
2749  int iu90 = (ntry*95.)/100.; if(iu90>=ntry) iu90=ntry-1;
2750 
2751  wl50[n][j] = wmed[n][j]>0 ? wmed[n][j]-value[index[il50]] : value[index[iu50]]-wmed[n][j];
2752  wu50[n][j] = wmed[n][j]>0 ? value[index[iu50]]-wmed[n][j] : wmed[n][j]-value[index[il50]];
2753  wl90[n][j] = wmed[n][j]>0 ? wmed[n][j]-value[index[il90]] : value[index[iu90]]-wmed[n][j];
2754  wu90[n][j] = wmed[n][j]>0 ? value[index[iu90]]-wmed[n][j] : wmed[n][j]-value[index[il90]];
2755  }
2756  delete [] index;
2757  delete [] value;
2758 
2759  vMED.push_back(wmed[n]);
2760  vL50.push_back(wl50[n]);
2761  vU50.push_back(wu50[n]);
2762  vL90.push_back(wl90[n]);
2763  vU90.push_back(wu90[n]);
2764  }
2765 
2766  // compute fAVR & fRMS
2767  wavearray<double> favr[NIFO_MAX];
2768  wavearray<double> frms[NIFO_MAX];
2769  for(int n=0;n<nIFO;n++) {
2770  frms[n] = frec[0][n]; frms[n] = 0;
2771  favr[n] = frec[0][n]; favr[n] = 0;
2772  for(int i=0;i<gOPT.trials;i++) {
2773  favr[n] += frec[i][n];
2774  for(int j=0;j< frec[i][n].size();j++) frms[n][j] += frec[i][n][j]*frec[i][n][j];
2775  }
2776  favr[n] *= 1./nrec[n];
2777 
2778  frms[n] *= 1./nrec[n];
2779  for(int j=0;j< frms[n].size();j++) {
2780  frms[n][j] -= favr[n][j]*favr[n][j];
2781  frms[n][j] = sqrt(frms[n][j]);
2782  }
2783 
2784  fAVR.push_back(favr[n]);
2785  fRMS.push_back(frms[n]);
2786  }
2787 
2788  // compute fMED, fL50, fU50, fL90, fU90
2789  wavearray<double> fmed[NIFO_MAX];
2790  wavearray<double> fl50[NIFO_MAX];
2791  wavearray<double> fu50[NIFO_MAX];
2792  wavearray<double> fl90[NIFO_MAX];
2793  wavearray<double> fu90[NIFO_MAX];
2794  for(int n=0;n<nIFO;n++) {
2795  fmed[n] = frec[0][n]; fmed[n] = 0;
2796  fl50[n] = frec[0][n]; fl50[n] = 0;
2797  fu50[n] = frec[0][n]; fu50[n] = 0;
2798  fl90[n] = frec[0][n]; fl90[n] = 0;
2799  fu90[n] = frec[0][n]; fu90[n] = 0;
2800 
2801  int ntry = gPE.trials; // number of detected events in the trials
2802  int *index = new int[ntry];
2803  float *value = new float[ntry];
2804  for(int j=0;j<frec[0][n].size();j++) {
2805 
2806  int k=0; for(int i=0;i<gOPT.trials;i++) if(wREC[i].size()) value[k++] = frec[i][n][j]; // select detected events
2807  TMath::Sort(ntry,value,index,false);
2808 
2809  int imed = (ntry*50.)/100.; if(imed>=ntry) imed=ntry-1;
2810  fmed[n][j] = value[index[imed]];
2811 
2812  int il50 = (ntry*25.)/100.; if(il50>=ntry) il50=ntry-1;
2813  int iu50 = (ntry*75.)/100.; if(iu50>=ntry) iu50=ntry-1;
2814  int il90 = (ntry*5.)/100.; if(il90>=ntry) il90=ntry-1;
2815  int iu90 = (ntry*95.)/100.; if(iu90>=ntry) iu90=ntry-1;
2816 
2817  fl50[n][j] = fmed[n][j]>0 ? fmed[n][j]-value[index[il50]] : value[index[iu50]]-fmed[n][j];
2818  fu50[n][j] = fmed[n][j]>0 ? value[index[iu50]]-fmed[n][j] : fmed[n][j]-value[index[il50]];
2819  fl90[n][j] = fmed[n][j]>0 ? fmed[n][j]-value[index[il90]] : value[index[iu90]]-fmed[n][j];
2820  fu90[n][j] = fmed[n][j]>0 ? value[index[iu90]]-fmed[n][j] : fmed[n][j]-value[index[il90]];
2821  }
2822  delete [] index;
2823  delete [] value;
2824 
2825  fMED.push_back(fmed[n]);
2826  fL50.push_back(fl50[n]);
2827  fU50.push_back(fu50[n]);
2828  fL90.push_back(fl90[n]);
2829  fU90.push_back(fu90[n]);
2830  }
2831 
2832  // ---------------------------------------------------
2833  // compute residuals energy
2834  // ---------------------------------------------------
2835  double eavr=0; // energy avr^2
2836  for(int n=0;n<nIFO;n++) {
2837  for(int j=0;j< wavr[n].size();j++) eavr += pow(wavr[n][j],2);
2838  }
2839  for(int i=0;i<gOPT.trials;i++) {
2840  if(wREC[i].size()) { // event detected
2841  double eres=0; // residual energy (rec-avr)^2
2842  for(int n=0;n<nIFO;n++) {
2843  for(int j=0;j< wavr[n].size();j++) {
2844  eres += pow(wrec[i][n][j]-wavr[n][j],2);
2845  }
2846  }
2847  // add normalized residual energy of wrec to vRES
2848  vRES.push_back(eres/eavr);
2849  }
2850  }
2851  double jres=0; // residual energy (inj-avr)^2
2852  if(cfg->simulation) {
2853  for(int n=0;n<nIFO;n++) {
2854  for(int j=0;j< wavr[n].size();j++) {
2855  jres += pow(winj[n][j]-wavr[n][j],2);
2856  }
2857  }
2858  jres/=eavr;
2859  }
2860  // add normalized residual energy of winj to vRES
2861  vRES.push_back(jres);
2862 
2863  // compute probability distribution of residuals
2864  int size = vRES.size()-1;
2865  wavearray<int> index(size);
2866  wavearray<double> eres(size);
2867  for(int i=0;i<size;i++) eres[i] = vRES[i];
2868  TMath::Sort(size,const_cast<double*>(eres.data),index.data,false);
2869  gPE.erR[0] = jres;
2870  for(int i=1;i<10;i++) gPE.erR[i]=eres[index[int(i*size/10.)]];
2871  gPE.erR[10]=0;
2872  for(int i=0;i<size;i++) {
2873  if(eres[i]<=jres) gPE.erR[10]+=1.;
2874  }
2875  gPE.erR[10]/=size;
2876  // print probability distribution of residuals
2877  cout.precision(3);
2878  cout << endl << "gPE.erR : ";
2879  for(int i=0;i<10;i++) cout << gPE.erR[i] << "/";
2880  cout << gPE.erR[10] << endl << endl;
2881 
2882  // ---------------------------------------------------
2883  // compute frequency residuals
2884  // ---------------------------------------------------
2885  eavr=0; // energy wavr^4*favr^2
2886  for(int n=0;n<nIFO;n++) {
2887  for(int j=0;j< wavr[n].size();j++) eavr += pow(wavr[n][j],4)*pow(favr[n][j],2);
2888  }
2889  for(int i=0;i<gOPT.trials;i++) {
2890  if(wREC[i].size()) { // event detected
2891  double eres=0; // weighted frequency residual energy wavr^4*(frec-favr)^2
2892  for(int n=0;n<nIFO;n++) {
2893  for(int j=0;j<wavr[n].size();j++) {
2894  eres += pow(wavr[n][j],4)*pow(frec[i][n][j]-favr[n][j],2);
2895  }
2896  }
2897  // add normalized frequency residual energy of frec to fRES
2898  fRES.push_back(eres/eavr);
2899  }
2900  }
2901  jres=0; // weighted frequency residual energy wavr^4*(finj-favr)^2
2902  if(cfg->simulation) {
2903  for(int n=0;n<nIFO;n++) {
2904  for(int j=0;j< wavr[n].size();j++) {
2905  jres += pow(wavr[n][j],4)*pow(finj[n][j]-favr[n][j],2);
2906  }
2907  }
2908  jres/=eavr;
2909  }
2910  // add normalized frequency residual energy of finj to fRES
2911  fRES.push_back(jres);
2912 
2913  // compute probability distribution of residuals
2914  for(int i=0;i<size;i++) eres[i] = fRES[i];
2915  TMath::Sort(size,const_cast<double*>(eres.data),index.data,false);
2916  gPE.erF[0] = jres;
2917  for(int i=1;i<10;i++) gPE.erF[i]=eres[index[int(i*size/10.)]];
2918  gPE.erF[10]=0;
2919  for(int i=0;i<size;i++) {
2920  if(eres[i]<=jres) gPE.erF[10]+=1.;
2921  }
2922  gPE.erF[10]/=size;
2923  // print probability distribution of frequency residuals
2924  cout.precision(3);
2925  cout << endl << "gPE.erF : ";
2926  for(int i=0;i<10;i++) cout << gPE.erF[i] << "/";
2927  cout << gPE.erF[10] << endl << endl;
2928 
2929  return nrec;
2930 }
2931 
2933 
2934  wavearray<double> wf = *wf2;
2935  wf=0;
2936 
2937  if(wf1==NULL) return wf;
2938  if(wf1->size()==0) return wf;
2939 
2940  double R=wf1->rate();
2941 
2942  double b_wf1 = wf1->start();
2943  double e_wf1 = wf1->start()+wf1->size()/R;
2944  double b_wf2 = wf2->start();
2945  double e_wf2 = wf2->start()+wf2->size()/R;
2946 
2947  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
2948  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
2949 
2950  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
2951  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
2952  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
2953 
2954  for(int i=0;i<sizeXCOR;i++) wf[i+o_wf2] = wf1->data[i+o_wf1];
2955 
2956  return wf;
2957 }
2958 
2960 
2961  wavearray<double> wfq;
2962 
2963  if(wf==NULL) return wfq;
2964  if(wf->size()==0) return wfq;
2965 
2966  wfq = *wf;
2967 
2968  CWB::mdc::PhaseShift(wfq,90);
2969  for(int i=0;i<wf->size();i++) wfq[i] = sqrt(pow(wf->data[i],2)+pow(wfq[i],2));
2970 
2971  return wfq;
2972 }
2973 
2975 
2976  std::vector<float>* vP;
2977  std::vector<int>* vI;
2978 
2979  skymap skyprob = NET->getifo(0)->tau;
2980  skyprob = 0;
2981 
2982  // save the probability skymap
2983  if(NET->wc_List[0].p_Map.size()) {
2984 
2985  vP = &(NET->wc_List[0].p_Map[id-1]);
2986  vI = &(NET->wc_List[0].p_Ind[id-1]);
2987 
2988  skyprob = NET->getifo(0)->tau;
2989  skyprob = 0.;
2990 
2991  for(int j=0; j<int(vP->size()); j++) {
2992  int i = (*vI)[j];
2993  double th = skyprob.getTheta(i);
2994  double ph = skyprob.getPhi(i);
2995  int k=skyprob.getSkyIndex(th, ph);
2996  skyprob.set(k,(*vP)[j]);
2997  }
2998  }
2999 
3000  return skyprob;
3001 }
3002 
3003 void SaveSkyProb(network* NET, CWB::config* cfg, int id) {
3004 
3005  std::vector<float>* vP;
3006  std::vector<int>* vI;
3007 
3008  // save the probability skymap
3009  if(NET->wc_List[0].p_Map.size()) {
3010 
3011  vP = &(NET->wc_List[0].p_Map[id-1]);
3012  vI = &(NET->wc_List[0].p_Ind[id-1]);
3013 
3014  gSKYPROB = NET->getifo(0)->tau;
3015  gSKYPROB = 0.;
3016 
3017  for(int j=0; j<int(vP->size()); j++) {
3018  int i = (*vI)[j];
3019  double th = gSKYPROB.getTheta(i);
3020  double ph = gSKYPROB.getPhi(i);
3021  int k=gSKYPROB.getSkyIndex(th, ph);
3022  gSKYPROB.set(k,(*vP)[j]);
3023  }
3024  }
3025 
3026  gHSKYPROB.SetBins(gSKYPROB.size(),0,gSKYPROB.size()-1);
3027 
3028  double PROB=0;
3029  for(int l=0;l<gSKYPROB.size();l++) PROB+=gSKYPROB.get(l);
3030  cout << "PROB : " << PROB << endl;
3031  for(int l=0;l<gSKYPROB.size();l++) gHSKYPROB.SetBinContent(l,gSKYPROB.get(l));
3032 
3033 /*
3034  wavearray<int> index(gSKYPROB.size());
3035  wavearray<float> skyprob(gSKYPROB.size());
3036  for(int l=0;l<gSKYPROB.size();l++) skyprob[l]=gSKYPROB.get(l);
3037  TMath::Sort(int(gSKYPROB.size()),const_cast<float*>(skyprob.data),index.data,true);
3038  double THR = PROB>0.998 ? 0.998 : PROB;
3039  double prob=0;
3040  int L=0;
3041  skymap SkyMask = gSKYPROB;
3042  for(int l=0;l<gSKYPROB.size();l++) {
3043  prob+=skyprob[index[l]];
3044  if(prob>THR && L>1000) SkyMask.set(index[l],0); else {SkyMask.set(index[l],1);L++;}
3045  }
3046  cout << "PROB(0.99) : " << prob << " " << L << endl;
3047  //NET->setSkyMask(SkyMask,'e');
3048  //NET->setIndexMode(0);
3049 */
3050 
3051 #ifdef SAVE_SKYPROB
3052  // plot gSKYPROB
3054  //gSM.SetOptions("hammer","Celestial",2);
3055  gSM.SetOptions("","Geographic");
3056  gSM.SetZaxisTitle("probability");
3057  gSM.Draw(1);
3058  //TH2D* hsm = gSM.GetHistogram();
3059  //hsm->GetZaxis()->SetTitleOffset(0.85);
3060  //hsm->GetZaxis()->SetTitleSize(0.03);
3061  gSM.Print(TString::Format("%s/gSkyprob_%d.png",".",id));
3062 #endif
3063 }
3064 
3065 void PlotResiduals(int ID, TString pdir, int sim, char type) {
3066 
3067  std::vector<double > xRES = type=='f' ? fRES : vRES;
3068 
3069  int size = xRES.size()-1; // last index of xRES contains the residual of injection
3070  double jres = xRES[size]; // residual of injection
3071  if(size==0) return;
3072 
3073  double hmin=1e20;
3074  double hmax=0;
3075  for(int i=0;i<size;i++) {if(xRES[i]>hmax) hmax=xRES[i]; if(xRES[i]<hmin) hmin=xRES[i];}
3076  hmin-=0.4*(hmax-hmin);
3077  hmax+=0.4*(hmax-hmin);
3078  hmax*=1.2;if(hmax>1) hmax=1;
3079 
3080  // plot residuals
3081  TCanvas* cres = new TCanvas("residuals", "residuals", 200, 20, 600, 600);
3082  TH1F* hres = new TH1F("residuals","residuals",100,hmin,hmax);
3083  for(int i=0;i<size;i++) hres->Fill(xRES[i]);
3084 
3085  // make it cumulative
3086  double integral = hres->ComputeIntegral();
3087  if (integral==0) {cout << "PlotResiduals : Empty histogram !!!" << endl;return;}
3088  double* cumulative = hres->GetIntegral();
3089  for (int i=0;i<hres->GetNbinsX();i++) hres->SetBinContent(i,cumulative[i]/integral);
3090  //cout << "integral " << integral << endl;
3091  hres->SetLineWidth(2);
3092 
3093  hres->Draw();
3094 
3095  hres->GetXaxis()->SetTitle("residuals");
3096  hres->GetXaxis()->SetTitleOffset(1.4);
3097  hres->GetXaxis()->SetLabelOffset(0.02);
3098  hres->GetXaxis()->SetNoExponent();
3099  hres->GetXaxis()->SetMoreLogLabels();
3100  if(type=='f') {
3101  hres->GetXaxis()->SetRangeUser(0.001,1.0);
3102  } else {
3103  hres->GetXaxis()->SetRangeUser(0.01,1.0);
3104  }
3105 
3106  hres->GetYaxis()->SetTitle("probability");
3107  hres->GetYaxis()->SetTitleOffset(1.6);
3108  hres->GetYaxis()->SetLabelOffset(0.02);
3109  hres->GetYaxis()->SetNoExponent();
3110  hres->GetYaxis()->SetRangeUser(0.01,1.0);
3111 
3112  cres->SetLogx();
3113  cres->SetLogy();
3114  cres->SetGridx();
3115  cres->SetGridy();
3116 
3117  // draw vertical line at the jres value
3118  float ymax = hres->GetMaximum();
3119  TLine *line = new TLine(jres,0,jres,ymax);
3120  line->SetLineWidth(2);
3121  line->SetLineColor(kRed);
3122  if(sim) line->Draw();
3123 
3124  // print probability distribution of residuals
3125  cout.precision(3);
3126  if(type=='f') {
3127  cout << endl << "---------> gPE.erF : ";
3128  for(int i=0;i<10;i++) cout << gPE.erF[i] << "/";
3129  cout << gPE.erF[10] << endl << endl;
3130  } else {
3131  cout << endl << "---------> gPE.erR : ";
3132  for(int i=0;i<10;i++) cout << gPE.erR[i] << "/";
3133  cout << gPE.erR[10] << endl << endl;
3134  }
3135 
3136  gStyle->SetLineColor(kBlack);
3137  char title[256];
3138  if(type=='f') {
3139  sprintf(title,"Cumulative distribution of freq residuals errors (entries=%d - prob inj = %2.2g)",size,gPE.erF[10]);
3140  } else {
3141  sprintf(title,"Cumulative distribution of residuals errors (entries=%d - prob inj = %2.2g)",size,gPE.erR[10]);
3142  }
3143  hres->SetTitle(title);
3144  hres->SetStats(kFALSE);
3145  if(type=='f') {
3146  cres->Print(TString::Format("%s/fresiduals_errors.%s",pdir.Data(),PLOT_TYPE));
3147  } else {
3148  cres->Print(TString::Format("%s/residuals_errors.%s",pdir.Data(),PLOT_TYPE));
3149  }
3150 
3151  delete line;
3152  delete hres;
3153  delete cres;
3154 }
3155 
3156 void PlotSNRnet(int nIFO, TString pdir, int sim) {
3157 
3158  int size = sim ? iSNR[0].size() : vLIKELIHOOD.size();
3159  if(size==0) return;
3160 
3161  double hmin=1e20;
3162  double hmax=0;
3163  for(int i=0;i<size;i++) {
3164  double osnr=0;
3165  if(sim) for(int n=0;n<nIFO;n++) osnr+=oSNR[n][i];
3166  else osnr+=vLIKELIHOOD[i];
3167  osnr=sqrt(osnr);
3168  if(osnr>hmax) hmax=osnr;
3169  if(osnr<hmin) hmin=osnr;
3170  }
3171  hmin-=0.4*(hmax-hmin);
3172  hmax+=0.4*(hmax-hmin);
3173 
3174  double isnr=0;
3175  if(sim) {
3176  for(int n=0;n<nIFO;n++) isnr+=iSNR[n][0];
3177  isnr=sqrt(isnr);
3178  if(isnr<hmin) hmin=0.8*isnr;
3179  if(isnr>hmax) hmax=1.2*isnr;
3180  }
3181 
3182  // plot SNRnet
3183  TCanvas* csnr = new TCanvas("SNRnet", "SNRnet", 200, 20, 600, 600);
3184  TH1F* hsnr = new TH1F("SNRnet","SNRnet",10,hmin,hmax);
3185  for(int i=0;i<size;i++) {
3186  double osnr=0;
3187  if(sim) for(int n=0;n<nIFO;n++) osnr+=oSNR[n][i];
3188  else osnr+=vLIKELIHOOD[i];
3189  hsnr->Fill(sqrt(osnr));
3190  }
3191  hsnr->SetLineWidth(2);
3192 
3193  hsnr->Draw();
3194 
3195  hsnr->GetXaxis()->SetTitle("SNRnet");
3196  hsnr->GetXaxis()->SetTitleOffset(1.4);
3197  hsnr->GetXaxis()->SetLabelOffset(0.02);
3198  hsnr->GetXaxis()->SetNoExponent();
3199 
3200  hsnr->GetYaxis()->SetTitle("counts");
3201  hsnr->GetYaxis()->SetTitleOffset(1.6);
3202  hsnr->GetYaxis()->SetLabelOffset(0.02);
3203  hsnr->GetYaxis()->SetNoExponent();
3204 
3205  csnr->SetLogy();
3206  csnr->SetGridx();
3207  csnr->SetGridy();
3208 
3209  // draw vertical line at the jres value
3210  float ymax = hsnr->GetMaximum();
3211  TLine *line = new TLine(isnr,0,isnr,ymax);
3212  line->SetLineWidth(2);
3213  line->SetLineColor(kRed);
3214  if(sim) line->Draw();
3215 
3216  hsnr->Fit("gaus","q0");
3217  TF1* gauss = (TF1*)hsnr->GetFunction("gaus");
3218  if(gauss) {
3219  gauss->Draw("SAME");
3220  gauss->SetLineColor(kGreen);
3221  gauss->SetLineWidth(2);
3222  }
3223 
3224  char title[256];
3225  if(sim) sprintf(title,"SNRnet Distribution (entries=%d - SNRnet inj = %2.2g)",size,isnr);
3226  else sprintf(title,"SNRnet Distribution (entries=%d)",size);
3227  hsnr->SetTitle(title);
3228  hsnr->SetStats(kTRUE);
3229  csnr->Print(TString::Format("%s/SNRnet.%s",pdir.Data(),PLOT_TYPE));
3230  hsnr->SetStats(kFALSE);
3231 
3232  delete line;
3233  delete hsnr;
3234  delete csnr;
3235 }
3236 
3237 void PlotChirpMass(int gtype, TString pdir, int sim) {
3238 
3239  // gtype=1 -> reconstructed chirp mass
3240  // gtype=2 -> reconstructed chirp mass error
3241  // gtype=3 -> ellipticity parameter
3242  // gtype=4 -> pixel fraction
3243  // gtype=5 -> energy fraction
3244 
3245  if(gtype<1 || gtype>5) return;
3246 
3247  TString gname = "";
3248  if(gtype==1) gname="chirp_mass";
3249  if(gtype==2) gname="chirp_mass_err";
3250  if(gtype==3) gname="chirp_ell";
3251  if(gtype==4) gname="chirp_pfrac";
3252  if(gtype==5) gname="chirp_efrac";
3253 
3254  int size = vCHIRP[gtype].size();
3255  if(size==0) return;
3256 
3257  double hmin=1e20;
3258  double hmax=0;
3259  for(int i=0;i<size;i++) {
3260  double omchirp = vCHIRP[gtype][i]; // reconstructed mchirp
3261  if(omchirp>hmax) hmax=omchirp;
3262  if(omchirp<hmin) hmin=omchirp;
3263  }
3264  hmin-=0.4*(hmax-hmin);
3265  hmax+=0.4*(hmax-hmin);
3266 
3267  double imchirp=vCHIRP[0][0]; // injected mchirp
3268  if(sim && gtype==1) {
3269  if(imchirp<hmin) hmin=0.8*imchirp;
3270  if(imchirp>hmax) hmax=1.2*imchirp;
3271  }
3272 
3273  // plot mchirp
3274  TCanvas* cmchirp = new TCanvas("mchirp", "mchirp", 200, 20, 600, 600);
3275  TH1F* hmchirp = new TH1F("mchirp","mchirp",10,hmin,hmax);
3276  for(int i=0;i<size;i++) hmchirp->Fill(vCHIRP[gtype][i]);
3277  hmchirp->SetLineWidth(2);
3278 
3279  hmchirp->Draw();
3280 
3281  hmchirp->GetXaxis()->SetTitle(gname);
3282  hmchirp->GetXaxis()->SetTitleOffset(1.4);
3283  hmchirp->GetXaxis()->SetLabelOffset(0.02);
3284  hmchirp->GetXaxis()->SetNoExponent();
3285 
3286  hmchirp->GetYaxis()->SetTitle("counts");
3287  hmchirp->GetYaxis()->SetTitleOffset(1.6);
3288  hmchirp->GetYaxis()->SetLabelOffset(0.02);
3289  hmchirp->GetYaxis()->SetNoExponent();
3290 
3291  cmchirp->SetLogy();
3292  cmchirp->SetGridx();
3293  cmchirp->SetGridy();
3294 
3295  // draw vertical line at the jres value
3296  float ymax = hmchirp->GetMaximum();
3297  TLine *line = new TLine(imchirp,0,imchirp,ymax);
3298  line->SetLineWidth(2);
3299  line->SetLineColor(kRed);
3300  if(sim && gtype==1) line->Draw();
3301 /*
3302  hmchirp->Fit("gaus","q0");
3303  TF1* gaus = (TF1*)hmchirp->GetFunction("gaus");
3304  if(gaus) {
3305  gaus->Draw("SAME");
3306  gaus->SetLineColor(kGreen);
3307  gaus->SetLineWidth(2);
3308  }
3309 */
3310  gStyle->SetLineColor(kBlack);
3311  char title[256];
3312  if(sim && gtype==1) sprintf(title,"%s distribution (entries=%d - inj chirp mass = %2.2g)",gname.Data(),size,imchirp);
3313  else sprintf(title,"%s distribution (entries=%d)",gname.Data(),size);
3314  hmchirp->SetTitle(title);
3315  hmchirp->SetStats(kTRUE);
3316  cmchirp->Print(TString::Format("%s/%s.%s",pdir.Data(),gname.Data(),PLOT_TYPE));
3317  hmchirp->SetStats(kFALSE);
3318 
3319  delete line;
3320  delete hmchirp;
3321  delete cmchirp;
3322 }
3323 
3324 void PlotFactors(int gtype, int nIFO, TString pdir) {
3325 
3326  // if gtype=0 -> Fitting Factor
3327  // if gtype=1 -> Overlap Factor
3328 
3329  if(gtype!=0 && gtype!=1) return;
3330 
3331  int size = iSNR[0].size();
3332  if(size==0) return;
3333 
3334  double hmin=1e20;
3335  double hmax=0;
3336  for(int i=0;i<size;i++) {
3337  double isnr=0; for(int n=0;n<nIFO;n++) isnr+= iSNR[n][i];
3338  double osnr=0; for(int n=0;n<nIFO;n++) osnr+= oSNR[n][i];
3339  double iosnr=0; for(int n=0;n<nIFO;n++) iosnr+=ioSNR[n][i];
3340  double ff = iosnr/sqrt(isnr*isnr); // fitting factor
3341  double of = iosnr/sqrt(isnr*osnr); // overlap factor
3342  double factor = gtype==0 ? ff : of;
3343  if(factor>hmax) hmax=factor; if(factor<hmin) hmin=factor;
3344  }
3345  hmin-=0.4*(hmax-hmin);
3346  hmax+=0.4*(hmax-hmin);
3347 
3348  // plot Factors
3349  TCanvas* cf = new TCanvas("Factors", "Factors", 200, 20, 600, 600);
3350  TH1F* hf = new TH1F("Factors","Factors",10,hmin,hmax);
3351  for(int i=0;i<size;i++) {
3352  double isnr=0; for(int n=0;n<nIFO;n++) isnr+= iSNR[n][i];
3353  double osnr=0; for(int n=0;n<nIFO;n++) osnr+= oSNR[n][i];
3354  double iosnr=0; for(int n=0;n<nIFO;n++) iosnr+=ioSNR[n][i];
3355  double ff = iosnr/sqrt(isnr*isnr); // fitting factor
3356  double of = iosnr/sqrt(isnr*osnr); // overlap factor
3357  double factor = gtype==0 ? ff : of;
3358  hf->Fill(factor);
3359  }
3360  hf->SetLineWidth(2);
3361 
3362  hf->Draw();
3363 
3364  if(gtype==0) hf->GetXaxis()->SetTitle("Fitting Factor");
3365  if(gtype==1) hf->GetXaxis()->SetTitle("Overlap Factor");
3366  hf->GetXaxis()->SetTitleOffset(1.4);
3367  hf->GetXaxis()->SetLabelOffset(0.02);
3368  hf->GetXaxis()->SetNoExponent();
3369 
3370  hf->GetYaxis()->SetTitle("counts");
3371  hf->GetYaxis()->SetTitleOffset(1.6);
3372  hf->GetYaxis()->SetLabelOffset(0.02);
3373  hf->GetYaxis()->SetNoExponent();
3374 
3375  cf->SetLogy();
3376  cf->SetGridx();
3377  cf->SetGridy();
3378 
3379  char title[256];
3380  if(gtype==0) sprintf(title,"Fitting Factor Distribution (entries=%d)",size);
3381  if(gtype==1) sprintf(title,"Overlap Factor Distribution (entries=%d)",size);
3382  hf->SetTitle(title);
3383  gStyle->SetLineColor(kBlack);
3384  hf->SetStats(kTRUE);
3385  if(gtype==0) cf->Print(TString::Format("%s/FittingFactor.%s",pdir.Data(),PLOT_TYPE));
3386  if(gtype==1) cf->Print(TString::Format("%s/OverlapFactor.%s",pdir.Data(),PLOT_TYPE));
3387  hf->SetStats(kFALSE);
3388 
3389  delete hf;
3390  delete cf;
3391 }
3392 
3393 void GetNullStatistic(std::vector<double>* vNUL, std::vector<double>* vNSE, int ifoId, TString ifoName, TString gtype, TString pdir) {
3394 
3395  if(vNUL->size()==0 || vNSE->size()==0) return;
3396 
3397  double xmin=1e20;
3398  double xmax=0;
3399  for(int i=0;i<vNUL->size();i++) {
3400  if((*vNUL)[i]>xmax) xmax=(*vNUL)[i]; if((*vNUL)[i]<xmin) xmin=(*vNUL)[i];
3401  }
3402  xmin *= xmin>0 ? 0.5 : 1.5;
3403  xmax *= xmax>0 ? 1.5 : 0.5;
3404  if(fabs(xmin)>fabs(xmax)) xmax=fabs(xmin)*xmax/fabs(xmax); else xmin=fabs(xmax)*xmin/fabs(xmin);
3405 
3406  TH1F* hnul = new TH1F("null","null",10,xmin,xmax);
3407  for(int i=0;i<vNUL->size();i++) hnul->Fill((*vNUL)[i]);
3408  hnul->SetLineWidth(2);
3409  hnul->SetLineColor(kRed);
3410 
3411  TH1F* hnse = new TH1F("noise","noise",10,xmin,xmax);
3412  for(int i=0;i<vNSE->size();i++) hnse->Fill((*vNSE)[i]);
3413  double scale = double(hnul->GetEntries())/hnse->GetEntries();
3414  hnse->Scale(scale);
3415  hnse->SetLineWidth(2);
3416  hnse->SetLineColor(kBlack);
3417 
3418  double ymin=0.9;
3419  double ymax=0;
3420  for(int i=0;i<=hnul->GetNbinsX();i++) {
3421  double binc = hnul->GetBinContent(i);
3422  if(binc>ymax) ymax=binc;
3423  }
3424  for(int i=0;i<=hnse->GetNbinsX();i++) {
3425  double binc = hnse->GetBinContent(i);
3426  if(binc>ymax) ymax=binc;
3427  }
3428  ymax *= 1.1;
3429  hnse->GetYaxis()->SetRangeUser(ymin,ymax);
3430 
3431  double KS = hnul->Integral() ? hnul->KolmogorovTest(hnse,"N") : 0.;
3432  double AD = hnul->Integral() ? hnul->AndersonDarlingTest(hnse) : 0;
3433 
3434  hnul->Fit("gaus","q0");
3435  TF1* gaus = (TF1*)hnul->GetFunction("gaus");
3436  double chi2 = gaus ? gaus->GetChisquare()/gaus->GetNDF() : 0;
3437 
3438  if(pdir!="") {
3439  // plot null
3440  TCanvas* cnul = new TCanvas("cnull", "cnull", 200, 20, 600, 600);
3441 
3442  hnse->Draw();
3443  hnul->Draw("same");
3444 
3445  hnse->GetXaxis()->SetTitle(gtype);
3446  hnse->GetXaxis()->SetTitleOffset(1.4);
3447  hnse->GetXaxis()->SetLabelOffset(0.02);
3448  hnse->GetXaxis()->SetNoExponent();
3449 
3450  hnse->GetYaxis()->SetTitle("counts");
3451  hnse->GetYaxis()->SetTitleOffset(1.6);
3452  hnse->GetYaxis()->SetLabelOffset(0.02);
3453  hnse->GetYaxis()->SetNoExponent();
3454 
3455  cnul->SetLogy();
3456  cnul->SetGridx();
3457  cnul->SetGridy();
3458 
3459  hnse->Fit("gaus","q0");
3460  gaus = (TF1*)hnse->GetFunction("gaus");
3461  if(gaus) {
3462  gaus->Draw("SAME");
3463  gaus->SetLineColor(kGreen);
3464  gaus->SetLineWidth(2);
3465  }
3466 
3467  char title[256];
3468  sprintf(title,"%s : %s (entries=%d) : KS=%0.2f : AD=%0.2f : CHI2=%0.2f",ifoName.Data(),gtype.Data(),vNUL->size(),KS,AD,chi2);
3469  hnse->SetTitle(title);
3470  hnse->SetStats(kTRUE);
3471  cnul->Print(TString::Format("%s/%s_%s_dist.%s",pdir.Data(),ifoName.Data(),gtype.Data(),PLOT_TYPE));
3472  hnse->SetStats(kFALSE);
3473 
3474  delete hnul;
3475  delete hnse;
3476  delete cnul;
3477 
3478  } else { // fill gPE.nstat
3479 
3480  int I = 2*7*ifoId; // offset
3481  if(gtype=="null_90") I+=7;
3482 
3483  gPE.nstat[I+0] = hnul->GetMean();
3484  gPE.nstat[I+1] = hnul->GetRMS();
3485  gPE.nstat[I+2] = hnse->GetMean();
3486  gPE.nstat[I+3] = hnse->GetRMS();
3487  gPE.nstat[I+4] = chi2;
3488  gPE.nstat[I+5] = KS;
3489  gPE.nstat[I+6] = AD;
3490 
3491  cout << endl;
3492  cout << " Pixels Statistic : " << ifoName << " " << gtype << endl;
3493  cout << " Null Pixels Mean : " << gPE.nstat[I+0] << endl;
3494  cout << " Null Pixels RMS : " << gPE.nstat[I+1] << endl;
3495  cout << " Noise Pixels Mean : " << gPE.nstat[I+2] << endl;
3496  cout << " Noise Pixels RMS : " << gPE.nstat[I+3] << endl;
3497  cout << " Noise Pixels Chi2 : " << gPE.nstat[I+4] << endl;
3498  cout << " KolmogorovTest : " << gPE.nstat[I+5] << endl;
3499  cout << " AndersonDarlingTest : " << gPE.nstat[I+6] << endl;
3500  cout << endl;
3501 
3502  delete hnul;
3503  delete hnse;
3504  }
3505 }
3506 
3508 
3509  double a;
3510  double E=0.,T=0.;
3511  int size=(int)x.size();
3512  double rate=x.rate();
3513  for(int j=0;j<size;j++) {
3514  a = x[j];
3515  T += a*a*j/rate; // central time
3516  E += a*a; // energy
3517  }
3518  T = E>0 ? T/E : 0.5*size/rate;
3519 
3520  return T;
3521 }
3522 
3523 double GetTimeBoundaries(wavearray<double> x, double P, double& bT, double& eT) {
3524 
3525  if(P<0) P=0;
3526  if(P>1) P=1;
3527 
3528  int N = x.size();
3529 
3530  double E = 0; // signal energy
3531  double avr = 0; // average
3532  for(int i=0;i<N;i++) {avr+=i*x[i]*x[i]; E+=x[i]*x[i];}
3533  int M=int(avr/E); // central index
3534 
3535  // search range which contains percentage P of the total energy E
3536  int jB=0;
3537  int jE=N-1;
3538  double a,b;
3539  double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
3540  for(int j=1; j<N; j++) {
3541  a = ((M-j>=0)&&(M-j<N)) ? x[M-j] : 0.;
3542  b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
3543  if(a) jB=M-j;
3544  if(b) jE=M+j;
3545  sum += a*a+b*b;
3546  if(sum/E > P) break;
3547  }
3548 
3549  bT = x.start()+jB/x.rate();
3550  eT = x.start()+jE/x.rate();
3551 
3552  return eT-bT;
3553 }
3554 
3556 
3557  double a;
3558  double E=0.,F=0.;
3559  int size=(int)x.size();
3560  double rate=x.rate();
3561  x.FFTW(1);
3562  double dF=(rate/(double)size)/2.;
3563  for(int j=0;j<size;j+=2) {
3564  a = x[j]*x[j]+x[j+1]*x[j+1];
3565  F += a*j*dF; // central frequency
3566  E += a; // energy
3567  }
3568  F = E>0 ? F/E : 0.5*rate;
3569 
3570  return F;
3571 }
3572 
3573 double GetFrequencyBoundaries(wavearray<double> x, double P, double& bF, double& eF) {
3574 
3575  if(P<0) P=0;
3576  if(P>1) P=1;
3577 
3578  int N = x.size();
3579 
3580  x.FFTW(1);
3581 
3582  double E = 0; // signal energy
3583  double avr = 0; // average
3584  for(int i=0;i<N;i++) {avr+=i*x[i]*x[i]; E+=x[i]*x[i];}
3585  int M=int(avr/E); // central index
3586 
3587  // search range which contains percentage P of the total energy E
3588  int jB=0;
3589  int jE=N-1;
3590  double a,b;
3591  double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
3592  for(int j=1; j<N; j++) {
3593  a = ((M-j>=0)&&(M-j<N)) ? x[M-j] : 0.;
3594  b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
3595  if(a) jB=M-j;
3596  if(b) jE=M+j;
3597  sum += a*a+b*b;
3598  if(sum/E > P) break;
3599  }
3600 
3601  double dF=(x.rate()/(double)x.size())/2.;
3602 
3603  bF = jB*dF;
3604  eF = jE*dF;
3605 
3606  return eF-bF;
3607 }
3608 
3609 wavearray<double> GetCutWaveform(wavearray<double> x, double bT, double eT, double bF, double eF) {
3610 
3611  bT-=x.start();
3612  eT-=x.start();
3613 
3614  int size=(int)x.size();
3615 
3616  // cut time range bT,eT
3617  double T=0.;
3618  double dT = 1./x.rate();
3619  for(int j=0;j<size;j++) {
3620  T = j*dT;
3621  if(T<bT || T>eT) x[j]=0;
3622  }
3623 
3624  // cut frequency range bF,eF
3625  double F=0.;
3626  double dF=(x.rate()/(double)x.size())/2.;
3627  x.FFTW(1);
3628  for(int j=0;j<size;j+=2) {
3629  F = j*dF;
3630  if(F<bF || F>eF) {x[j]=0;x[j+1]=0;}
3631  }
3632  x.FFTW(-1);
3633 
3634  return x;
3635 }
3636 
3637 void SetWaveformCuts(wavearray<double>* x, double bT, double eT, double bF, double eF) {
3638 
3639  bT-=x->start();
3640  eT-=x->start();
3641 
3642  int size=(int)x->size();
3643 
3644  // cut time range bT,eT
3645  double T=0.;
3646  double dT = 1./x->rate();
3647  for(int j=0;j<size;j++) {
3648  T = j*dT;
3649  if(T<bT || T>eT) x->data[j]=0;
3650  }
3651 
3652  // cut frequency range bF,eF
3653  double F=0.;
3654  double dF=(x->rate()/(double)x->size())/2.;
3655  x->FFTW(1);
3656  for(int j=0;j<size;j+=2) {
3657  F = j*dF;
3658  if(F<bF || F>eF) {x->data[j]=0;x->data[j+1]=0;}
3659  }
3660  x->FFTW(-1);
3661 }
3662 
3663 TString DumpCED(network* NET, netevent* &EVT, CWB::config* cfg, double ofactor) {
3664 
3665  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
3666  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
3667 
3668  int nIFO = NET->ifoListSize(); // number of detectors
3669  double factor = cfg->factors[gIFACTOR];
3670 
3671  char sfactor[32];
3672  if(cfg->simulation==3) {
3673  if(factor<0) sprintf(sfactor,"n%g",fabs(factor));
3674  if(factor==0) sprintf(sfactor,"z%g",factor);
3675  if(factor>0) sprintf(sfactor,"p%g",factor);
3676  } else sprintf(sfactor,"%g",factor);
3677 
3678  char out_CED[1024];
3679  if(cfg->simulation) {
3680  char sim_label[512];
3681  sprintf(sim_label,"%d_%d_%s_%s_job%d",int(gCWB2G->Tb),int(gCWB2G->dT),cfg->data_label,sfactor,gCWB2G->runID);
3682  sprintf(out_CED,"%s/ced_%s_%d",cfg->nodedir,sim_label,gSystem->GetPid());
3683  } else {
3684  char prod_label[512];
3685  sprintf(prod_label,"%d_%d_%s_slag%d_lag%lu_%lu_job%d",
3686  int(gCWB2G->Tb),int(gCWB2G->dT),cfg->data_label,gCWB2G->slagID,cfg->lagOff,cfg->lagSize,gCWB2G->runID);
3687  sprintf(out_CED,"%s/ced_%s_%d",cfg->nodedir,prod_label,gSystem->GetPid());
3688  }
3689  cout << out_CED << endl;
3690 
3691  // restore whitened data into the detector TF map
3692  if(gOPT.noise==0) {
3693  for(int n=0; n<nIFO; n++) {
3694  WSeries<double>* pTF = NET->getifo(n)->getTFmap();
3695  *pTF = gHOT[n];
3696  }
3697  }
3698 
3699  cout<<"DumpCED : dump ced to disk ..." <<endl;
3700  CWB::ced ced(NET, EVT, out_CED);
3701  cfg->cedRHO=0;
3702  switch(gCWB2G->istage) {
3703  case CWB_STAGE_FULL:
3704  case CWB_STAGE_INIT:
3705  // use TF map & inj signals
3706  ced.SetOptions(cfg->simulation,cfg->cedRHO,cfg->inRate);
3707  for(int n=0; n<nIFO; n++) {gCWB2G->pTF[n]->setlow(cfg->fLow); gCWB2G->pTF[n]->sethigh(cfg->fHigh);}
3708  break;
3709  default:
3710  // use sparse map & inj signals
3711  ced.SetOptions(cfg->simulation,cfg->cedRHO,cfg->inRate,true);
3712  }
3713  if(gCWB2G->singleDetector) ced.SetChannelName(cfg->channelNamesRaw[0]);
3714  bool fullCED = gCWB2G->singleDetector ? false : true;
3715  ced.Write(ofactor,fullCED);
3716 
3717  char ifostr[20]="";
3718  if(gCWB2G->singleDetector) {
3719  sprintf(ifostr,"%s%s",ifostr,NET->ifoName[0]);
3720  } else {
3721  for(int n=0; n<nIFO; n++) sprintf(ifostr,"%s%s",ifostr,NET->ifoName[n]);
3722  }
3723 
3724  char dirCED[1024];
3725  sprintf(dirCED, "%s/%s", out_CED, ifostr);
3726  if(gCWB2G->singleDetector) {
3727  sprintf(dirCED, "%s_%.3f",dirCED,EVT->start[0]);
3728  } else {
3729  for(int n=0; n<nIFO; n++) sprintf(dirCED, "%s_%.3f",dirCED,EVT->start[n]);
3730  }
3731 
3732  return dirCED;
3733 }
3734 
3735 void CreateIndexHTML(TString dirCED, int nIFO, TString* ifo, bool sim) {
3736 
3737  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
3738  if(gCWB2G->singleDetector) nIFO=1;
3739 
3740  // ------------------------------------------------------------------------------------
3741  // expand html file paths
3742  // ------------------------------------------------------------------------------------
3743 
3744  TString pe_index = gSystem->ExpandPathName(PE_INDEX);
3745 
3746  // ------------------------------------------------------------------------------------
3747  // check if environmental variables are defined
3748  // ------------------------------------------------------------------------------------
3749  if(gSystem->Getenv("HOME_CED_WWW")==NULL) {
3750  cout << "Error : environment HOME_CED_WWW is not defined!!! " << endl;exit(1);}
3751 
3752  char ofileName[256]="";
3753  sprintf(ofileName,"%s/pe_index.html",dirCED.Data());
3754  cout << "ofileName : " << ofileName << endl;
3755 
3756  // -------------------------------------------------------------------------------
3757  // open output file
3758  // -------------------------------------------------------------------------------
3759  ofstream out;
3760  out.open(ofileName,ios::out);
3761  if(!out.good()) {cout << "Error Opening File : " << ofileName << endl;exit(1);}
3762  // HEADER
3763  WriteBodyHTML(pe_index, "PE_HEADER_BEG", "PE_HEADER_END", &out);
3764  // BODY
3765  WriteBodyHTML(pe_index, "PE_BODY_BEG", "PE_BODY_END", &out);
3766  // TABBERS
3767  int sbody_height = 0;
3768  out << "<div class=\"tabber\">" << endl;
3769  // -------------------------------------------------------------------------------
3770  // Time-Frequency Maps
3771  // -------------------------------------------------------------------------------
3772  if(gOPT.ced_tfmap) {
3773  out << "<div class=\"tabbertab\">" << endl;
3774  out << " <h2>Time-Frequency Maps</h2>" << endl;
3775  sbody_height = 1000+nIFO*400;
3776  if(gOPT.ced_pca) sbody_height += 600;
3777  out << " <iframe src=\"tfmap_body.html\" width=\"100%\" "
3778  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3779  out << "</div>" << endl;
3780  }
3781  // -------------------------------------------------------------------------------
3782  // Reconstructed Detector Responses
3783  // -------------------------------------------------------------------------------
3784  if(gOPT.ced_rdr) {
3785  out << "<div class=\"tabbertab\">" << endl;
3786  out << " <h2>Reconstructed Detector Responses</h2>" << endl;
3787  if(!sim) sbody_height = 350+nIFO*400;
3788  else sbody_height = 350+nIFO*1340;
3789  out << " <iframe src=\"rec_signal_body.html\" width=\"100%\" "
3790  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3791  out << "</div>" << endl;
3792  }
3793  // -------------------------------------------------------------------------------
3794  // Power Spectral Density
3795  // -------------------------------------------------------------------------------
3796  if(gOPT.ced_psd) {
3797  out << "<div class=\"tabbertab\">" << endl;
3798  out << " <h2>PSD</h2>" << endl;
3799  sbody_height = 350+(nIFO+1)*750;
3800  out << " <iframe src=\"psd_body.html\" width=\"100%\" "
3801  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3802  out << "</div>" << endl;
3803  }
3804  // -------------------------------------------------------------------------------
3805  // SkyMaps
3806  // -------------------------------------------------------------------------------
3807  sbody_height = 2050;
3808  if(gOPT.ced_skymap && !gCWB2G->singleDetector) {
3809  out << "<div class=\"tabbertab\">" << endl;
3810  out << " <h2>SkyMaps</h2>" << endl;
3811  out << " <iframe src=\"skymap_body.html\" width=\"100%\" "
3812  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3813  out << "</div>" << endl;
3814  }
3815  // -------------------------------------------------------------------------------
3816  // Reconstructed vs Data/Average
3817  // -------------------------------------------------------------------------------
3818  if(gOPT.ced_rec) {
3819  out << "<div class=\"tabbertab\">" << endl;
3820  out << " <h2>Waveform Errors</h2>" << endl;
3821  sbody_height = 370+nIFO*880;
3822  out << " <iframe src=\"rec_avr_signal_body.html\" width=\"100%\" "
3823  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3824  out << "</div>" << endl;
3825  }
3826  // -------------------------------------------------------------------------------
3827  // Reconstructed vs Injection
3828  // -------------------------------------------------------------------------------
3829  if(sim && gOPT.ced_inj) {
3830  out << "<div class=\"tabbertab\">" << endl;
3831  out << " <h2>Injection</h2>" << endl;
3832  sbody_height = 390+nIFO*880;
3833  out << " <iframe src=\"rec_inj_signal_body.html\" width=\"100%\" "
3834  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3835  out << "</div>" << endl;
3836  }
3837  // -------------------------------------------------------------------------------
3838  // -------------------------------------------------------------------------------
3839  // Reconstructed vs Reduced Injection
3840  // -------------------------------------------------------------------------------
3841  if(sim && gOPT.ced_rinj) {
3842  out << "<div class=\"tabbertab\">" << endl;
3843  out << " <h2>Injection</h2>" << endl;
3844  sbody_height = 390+nIFO*320;
3845  out << " <iframe src=\"rec_inj_signal_body.html\" width=\"100%\" "
3846  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3847  out << "</div>" << endl;
3848  }
3849  // -------------------------------------------------------------------------------
3850  // Chirp Mass Distributions
3851  // -------------------------------------------------------------------------------
3852  if(gOPT.ced_cm) {
3853  sbody_height = 2100;
3854  out << "<div class=\"tabbertab\">" << endl;
3855  out << " <h2>Chirp Mass</h2>" << endl;
3856  out << " <iframe src=\"mchirp_body.html\" width=\"100%\" "
3857  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3858  out << "</div>" << endl;
3859  }
3860  // -------------------------------------------------------------------------------
3861  // Distributions
3862  // -------------------------------------------------------------------------------
3863  if(gOPT.ced_distr) {
3864  sbody_height = 2000;
3865  out << "<div class=\"tabbertab\">" << endl;
3866  out << " <h2>Distributions</h2>" << endl;
3867  out << " <iframe src=\"distribution_body.html\" width=\"100%\" "
3868  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3869  out << "</div>" << endl;
3870  }
3871  // -------------------------------------------------------------------------------
3872  // Null Pixels Distribution
3873  // -------------------------------------------------------------------------------
3874  if(gOPT.ced_null) {
3875  sbody_height = 100+nIFO*1300;
3876  out << "<div class=\"tabbertab\">" << endl;
3877  out << " <h2>Null Pixels</h2>" << endl;
3878  out << " <iframe src=\"nstat_body.html\" width=\"100%\" "
3879  << " height=\"" << sbody_height << "px\" frameborder=\"0\"></iframe>" << endl;
3880  out << "</div>" << endl;
3881  }
3882  // -------------------------------------------------------------------------------
3883  // close output file
3884  // -------------------------------------------------------------------------------
3885  out << "</div>" << endl;
3886  out.close();
3887 
3888 
3889  if(gOPT.ced_tfmap) {
3890  out.open(TString::Format("%s/tfmap_body.html",dirCED.Data()).Data(),ios::out);
3891  // TFMAPS
3892  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3893  WriteBodyHTML(pe_index, "PE_TFMAP_HEADER_BEG", "PE_TFMAP_HEADER_END", &out);
3894  WriteBodyHTML(pe_index, "PE_TFMAP_BEG", "PE_TFMAP_END", &out, nIFO, ifo);
3895  // LIKELIHOOD
3896  out << "<p><br /></p> " << endl;
3897  out << "<p><br /></p> " << endl;
3898  WriteBodyHTML(pe_index, "PE_LIKELIHOOD_BEG", "PE_LIKELIHOOD_END", &out);
3899  // POLARGRAMS : NOT IMPLEMENTED IN WP !!!
3900  //out << "<p><br /></p> " << endl;
3901  //out << "<p><br /></p> " << endl;
3902  //WriteBodyHTML(pe_index, "PE_POLARGRAM_BEG", "PE_POLARGRAM_END", &out);
3903  if(gOPT.ced_pca) WriteBodyHTML(pe_index, "PE_PCA_BEG", "PE_PCA_END", &out);
3904  out.close();
3905  }
3906 
3907  // RECONSTRUCTED SIGNAL
3908  if(gOPT.ced_rdr) {
3909  out.open(TString::Format("%s/rec_signal_body.html",dirCED.Data()).Data(),ios::out);
3910  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3911  WriteBodyHTML(pe_index, "PE_REC_SIGNAL_HEADER_BEG", "PE_REC_SIGNAL_HEADER_END", &out);
3912  for(int n=0;n<nIFO;n++) {
3913  WriteBodyHTML(pe_index, "PE_IFO_SIGNAL_BEG", "PE_IFO_SIGNAL_END", &out, 1, &ifo[n]);
3914  if(sim) {
3915  WriteBodyHTML(pe_index, "PE_INJ_SIGNAL_BEG", "PE_INJ_SIGNAL_END", &out, 1, &ifo[n]);
3916  }
3917  WriteBodyHTML(pe_index, "PE_REC_SIGNAL_BEG", "PE_REC_SIGNAL_END", &out, 1, &ifo[n]);
3918  out << "<p><br /></p> " << endl;
3919  out << "<p><br /></p> " << endl;
3920  }
3921  out.close();
3922  }
3923 
3924  // PSD
3925  if(gOPT.ced_psd) {
3926  out.open(TString::Format("%s/psd_body.html",dirCED.Data()).Data(),ios::out);
3927  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3928  WriteBodyHTML(pe_index, "PE_PSD_HEADER_BEG", "PE_PSD_HEADER_END", &out);
3929  for(int n=0;n<nIFO;n++) {
3930  WriteBodyHTML(pe_index, "PE_IFO_PSD_BEG", "PE_IFO_PSD_END", &out, 1, &ifo[n]);
3931  WriteBodyHTML(pe_index, "PE_PSD_BEG", "PE_PSD_END", &out, 1, &ifo[n]);
3932  out << "<p><br /></p> " << endl;
3933  out << "<p><br /></p> " << endl;
3934  }
3935  out.close();
3936  }
3937 
3938  // SKYMAP
3939  if(gOPT.ced_skymap && !gCWB2G->singleDetector) {
3940  out.open(TString::Format("%s/skymap_body.html",dirCED.Data()).Data(),ios::out);
3941  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3942  WriteBodyHTML(pe_index, "PE_SKYMAP_BEG", "PE_SKYMAP_END", &out);
3943  out.close();
3944  }
3945 
3946  // RECONSTRUCTED vs DATA/AVERAGE
3947  if(gOPT.ced_rec) {
3948  out.open(TString::Format("%s/rec_avr_signal_body.html",dirCED.Data()).Data(),ios::out);
3949  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3950  WriteBodyHTML(pe_index, "PE_REC_AVR_SIGNAL_HEADER_BEG", "PE_REC_AVR_SIGNAL_HEADER_END", &out);
3951  for(int n=0;n<nIFO;n++) {
3952  WriteBodyHTML(pe_index, "PE_REC_AVR_SIGNAL_BEG", "PE_REC_AVR_SIGNAL_END", &out, 1, &ifo[n]);
3953  out << "<p><br /></p> " << endl;
3954  out << "<p><br /></p> " << endl;
3955  }
3956  out.close();
3957  }
3958 
3959  // RECONSTRUCTED vs INJECTED
3960  if(sim && (gOPT.ced_inj || gOPT.ced_rinj)) {
3961  out.open(TString::Format("%s/rec_inj_signal_body.html",dirCED.Data()).Data(),ios::out);
3962  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3963  if(gOPT.ced_inj) WriteBodyHTML(pe_index, "PE_REC_INJ_SIGNAL_0_HEADER_BEG", "PE_REC_INJ_SIGNAL_0_HEADER_END", &out);
3964  if(gOPT.ced_rinj) WriteBodyHTML(pe_index, "PE_REC_INJ_SIGNAL_HEADER_BEG", "PE_REC_INJ_SIGNAL_HEADER_END", &out);
3965  if(sim) {
3966  for(int n=0;n<nIFO;n++) {
3967  if(gOPT.ced_inj) WriteBodyHTML(pe_index, "PE_REC_INJ_SIGNAL_0_BEG", "PE_REC_INJ_SIGNAL_0_END", &out, 1, &ifo[n]);
3968  if(gOPT.ced_rinj) WriteBodyHTML(pe_index, "PE_REC_INJ_SIGNAL_BEG", "PE_REC_INJ_SIGNAL_END", &out, 1, &ifo[n]);
3969  out << "<p><br /></p> " << endl;
3970  out << "<p><br /></p> " << endl;
3971  }
3972  }
3973  out.close();
3974  }
3975 
3976 
3977  // DISTRIBUTION
3978  if(gOPT.ced_distr) {
3979  out.open(TString::Format("%s/distribution_body.html",dirCED.Data()).Data(),ios::out);
3980  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3981  WriteBodyHTML(pe_index, "PE_DISTRIBUTION_BEG", "PE_DISTRIBUTION_END", &out);
3982  if(sim) WriteBodyHTML(pe_index, "PE_FACTORS_DISTRIBUTION_BEG", "PE_FACTORS_DISTRIBUTION_END", &out);
3983  out.close();
3984  }
3985 
3986  // NULL PIXELS DISTRIBUTION
3987  if(gOPT.ced_null) {
3988  out.open(TString::Format("%s/nstat_body.html",dirCED.Data()).Data(),ios::out);
3989  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
3990  WriteBodyHTML(pe_index, "PE_NULPIX_DISTRIBUTION_HEADER_BEG", "PE_NULPIX_DISTRIBUTION_HEADER_END", &out);
3991  for(int n=0;n<nIFO;n++) {
3992  WriteBodyHTML(pe_index, "PE_NULPIX_DISTRIBUTION_BEG", "PE_NULPIX_DISTRIBUTION_END", &out, 1, &ifo[n]);
3993  out << "<p><br /></p> " << endl;
3994  out << "<p><br /></p> " << endl;
3995  }
3996  out.close();
3997  }
3998 
3999  // MCHIRP DISTRIBUTION
4000  if(gOPT.ced_cm) {
4001  out.open(TString::Format("%s/mchirp_body.html",dirCED.Data()).Data(),ios::out);
4002  WriteBodyHTML(pe_index, "PE_JSCRIPT_BEG", "PE_JSCRIPT_END", &out);
4003  WriteBodyHTML(pe_index, "PE_MCHIRP_DISTRIBUTION_HEADER_BEG", "PE_MCHIRP_DISTRIBUTION_HEADER_END", &out);
4004  WriteBodyHTML(pe_index, "PE_MCHIRP_DISTRIBUTION_BEG", "PE_MCHIRP_DISTRIBUTION_END", &out);
4005  out.close();
4006  }
4007 
4008  // copy tabber javascript to CED directory
4009  char cmd[1024];
4010  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.css %s",dirCED.Data());
4011  gSystem->Exec(cmd);
4012  sprintf(cmd,"cp ${HOME_WAT}//html/etc/html/tabber.js %s",dirCED.Data());
4013  gSystem->Exec(cmd);
4014  // overwrite CED index.html
4015  sprintf(cmd,"mv %s/pe_index.html %s/index.html",dirCED.Data(),dirCED.Data());
4016  gSystem->Exec(cmd);
4017 
4018  return;
4019 }
4020 
4021 void WriteBodyHTML(TString html_template, TString html_tag_beg, TString html_tag_end, ofstream* out, int nIFO, TString* ifo) {
4022 
4023  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
4024  if(gCWB2G->singleDetector) nIFO=1;
4025 
4026  TString home_ced_www=TString(gSystem->Getenv("HOME_CED_WWW"));
4027 
4028  // if CWB_DOC_URL is define then man infos are added to web pages
4029  TString cwb_doc_url="";
4030  if(gSystem->Getenv("CWB_DOC_URL")!=NULL) {
4031  cwb_doc_url=TString(gSystem->Getenv("CWB_DOC_URL"));
4032  }
4033 
4034  // get HOME_WWW
4035  TString home_www="~waveburst/LSC/waveburst"; // default value
4036  if(gSystem->Getenv("HOME_WWW")!=NULL) {
4037  home_www=TString(gSystem->Getenv("HOME_WWW"));
4038  }
4039 
4040  char line[1024];
4041  bool output=false;
4042  for(int n=0;n<nIFO;n++) {
4043  ifstream in;
4044  in.open(html_template.Data(),ios::in);
4045  if(!in.good()) {cout << "Error Opening File : " << html_template.Data() << endl;exit(1);}
4046  while(1) {
4047  in.getline(line,1024);
4048  if(!in.good()) break;
4049  TString sline = line;
4050  if(sline.Contains(html_tag_beg)&&!output) output=true;
4051  if(sline.Contains(html_tag_end)&& output) {output=false;break;}
4052  if(ifo!=NULL) sline.ReplaceAll("IFO",ifo[n]);
4053  sline.ReplaceAll("xNTRIALS",TString::Format("%d",(int)vCHIRP[0].size()-1).Data());
4054  sline.ReplaceAll("xINRATE",TString::Format("%d",(int)gINRATE));
4055  sline.ReplaceAll("xRATEANA",TString::Format("%d",(int)gRATEANA));
4056 #ifdef PLOT_MEDIAN
4057  sline.ReplaceAll("xMEAN","Median");
4058  sline.ReplaceAll("xSIGMA","50% Percentile");
4059  sline.ReplaceAll("x2SIGMA","90% Percentile");
4060 #else
4061  sline.ReplaceAll("xMEAN","Mean");
4062  sline.ReplaceAll("xSIGMA","RMS");
4063  sline.ReplaceAll("x2SIGMA","2RMS");
4064 #endif
4065  sline.ReplaceAll("HOME_CED_WWW",home_ced_www.Data());
4066  sline.ReplaceAll("CWB_DOC_URL",cwb_doc_url.Data());
4067  sline.ReplaceAll("HOME_WWW",home_www.Data());
4068  if(cwb_doc_url!="") {
4069  sline.ReplaceAll("<!--CWB_DOC_URL","");
4070  sline.ReplaceAll("CWB_DOC_URL-->","");
4071  sline.ReplaceAll("XCWB_DOC_URL",cwb_doc_url.Data());
4072  }
4073  if(output) (*out) << sline.Data() << endl;
4074  }
4075  in.close();
4076  }
4077 }
4078 
4079 void
4081 
4082  int n;
4083 
4084  n = ifo->IWFP.size();
4085  for (int i=0;i<n;i++) {
4087  delete wf;
4088  }
4089  ifo->IWFP.clear();
4090  ifo->IWFID.clear();
4091 
4092  n = ifo->RWFP.size();
4093  for (int i=0;i<n;i++) {
4095  delete wf;
4096  }
4097  ifo->RWFP.clear();
4098  ifo->RWFID.clear();
4099 }
4100 
4102  wavearray<double>* favr, wavearray<double>* ferr, wavearray<double>* wref, TString pdir, double P) {
4103 
4104  int size = frec->size();
4105 
4106  wavearray<double> time(size);
4107  wavearray<double> etime(size); etime=0;
4108  for (int i=0; i<size; i++) time.data[i] = i/frec->rate()+(wref->start()-gSEGGPS);
4109 
4110  wavearray<double> f2err=*ferr; f2err*=2.;
4111 
4112  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4113 
4114  double bT, eT;
4115  GetTimeBoundaries(*wref, P, bT, eT);
4116  bT-=gSEGGPS;
4117  eT-=gSEGGPS;
4118  for(int i=0;i<frec->size();i++) {
4119  double xtime = i/wref->rate();
4120  if(xtime>bT && xtime<eT) continue;
4121  frec->data[i]=0;
4122  favr->data[i]=0;
4123  ferr->data[i]=0;
4124  f2err.data[i]=0;
4125  }
4126 
4127  TGraphErrors* e2gr = new TGraphErrors(size,time.data,favr->data,etime.data,f2err.data);
4128  e2gr->SetLineColor(17);
4129  e2gr->SetFillStyle(1001);
4130  e2gr->SetFillColor(17);
4131  e2gr->GetXaxis()->SetTitle(xtitle);
4132  e2gr->GetYaxis()->SetTitle("frequency (hz)");
4133  e2gr->SetTitle(title);
4134  e2gr->GetXaxis()->SetTitleFont(42);
4135  e2gr->GetXaxis()->SetLabelFont(42);
4136  e2gr->GetXaxis()->SetLabelOffset(0.012);
4137  e2gr->GetXaxis()->SetTitleOffset(1.5);
4138  e2gr->GetYaxis()->SetTitleFont(42);
4139  e2gr->GetYaxis()->SetLabelFont(42);
4140  e2gr->GetYaxis()->SetLabelOffset(0.01);
4141  e2gr->GetYaxis()->SetTitleOffset(1.4);
4142 
4143 
4144  TGraphErrors* egr = new TGraphErrors(size,time.data,favr->data,etime.data,ferr->data);
4145  egr->SetLineColor(15);
4146  egr->SetFillStyle(1001);
4147  egr->SetFillColor(15);
4148 
4149  TGraph* agr = new TGraph(size,time.data,favr->data);
4150  agr->SetLineWidth(1);
4151  agr->SetLineColor(kWhite);
4152  agr->SetLineStyle(1);
4153 
4154  TGraph* gr = new TGraph(size,time.data,frec->data);
4155  gr->SetLineWidth(1);
4156  gr->SetLineColor(2);
4157 
4158  TCanvas* canvas = new TCanvas("frec", "frec",200,20,800,500);
4159  canvas->cd();
4160  canvas->SetGridx();
4161  canvas->SetGridy();
4162 
4163  e2gr->GetXaxis()->SetRangeUser(bT, eT);
4164  e2gr->Draw("AP4");
4165  egr->Draw("P4same");
4166  agr->Draw("Lsame");
4167  gr->Draw("Lsame");
4168 
4169  if(ofname!="") {
4170  ofname = TString(pdir)+TString("/")+ofname;
4171  canvas->Print(ofname);
4172  cout << "write : " << ofname << endl;
4173  }
4174 
4175  delete canvas;
4176  delete e2gr;
4177  delete egr;
4178  delete agr;
4179  delete gr;
4180 }
4181 
4183  wavearray<double>* wavr, wavearray<double>* werr, wavearray<double>* wref, TString pdir, double P) {
4184 
4185  int size = wrec->size();
4186 
4187  wavearray<double> time(size);
4188  wavearray<double> etime(size); etime=0;
4189  for (int i=0; i<size; i++) time.data[i] = i/wrec->rate()+(wref->start()-gSEGGPS);
4190 
4191  wavearray<double> w2err=*werr; w2err*=2.;
4192 
4193  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4194 
4195  TGraphErrors* e2gr = new TGraphErrors(size,time.data,wavr->data,etime.data,w2err.data);
4196  e2gr->SetLineColor(17);
4197  e2gr->SetFillStyle(1001);
4198  e2gr->SetFillColor(17);
4199  e2gr->GetXaxis()->SetTitle(xtitle);
4200  e2gr->GetYaxis()->SetTitle("magnitude");
4201  e2gr->SetTitle(title);
4202  e2gr->GetXaxis()->SetTitleFont(42);
4203  e2gr->GetXaxis()->SetLabelFont(42);
4204  e2gr->GetXaxis()->SetLabelOffset(0.012);
4205  e2gr->GetXaxis()->SetTitleOffset(1.5);
4206  e2gr->GetYaxis()->SetTitleFont(42);
4207  e2gr->GetYaxis()->SetLabelFont(42);
4208  e2gr->GetYaxis()->SetLabelOffset(0.01);
4209  e2gr->GetYaxis()->SetTitleOffset(1.4);
4210 
4211  TGraphErrors* egr = new TGraphErrors(size,time.data,wavr->data,etime.data,werr->data);
4212  egr->SetLineColor(15);
4213  egr->SetFillStyle(1001);
4214  egr->SetFillColor(15);
4215 
4216  TGraph* agr = new TGraph(size,time.data,wavr->data);
4217  agr->SetLineWidth(1);
4218  agr->SetLineColor(kWhite);
4219  agr->SetLineStyle(1);
4220 
4221  TGraph* gr = new TGraph(size,time.data,wrec->data);
4222  gr->SetLineWidth(1);
4223  gr->SetLineColor(2);
4224 
4225  TCanvas* canvas = new TCanvas("wrec", "wrec",200,20,800,500);
4226  canvas->cd();
4227  canvas->SetGridx();
4228  canvas->SetGridy();
4229 
4230  double bT, eT;
4231  GetTimeBoundaries(*wref, P, bT, eT);
4232  bT-=gSEGGPS;
4233  eT-=gSEGGPS;
4234  e2gr->GetXaxis()->SetRangeUser(bT, eT);
4235  e2gr->Draw("A3");
4236  egr->Draw("3same");
4237  agr->Draw("Lsame");
4238  gr->Draw("Lsame");
4239 
4240  if(ofname!="") {
4241  ofname = TString(pdir)+TString("/")+ofname;
4242  canvas->Print(ofname);
4243  cout << "write : " << ofname << endl;
4244  }
4245 
4246  delete canvas;
4247  delete e2gr;
4248  delete egr;
4249  delete agr;
4250  delete gr;
4251 }
4252 
4255  wavearray<double>* wl90, wavearray<double>* wu90, wavearray<double>* wref, TString pdir, double P, bool freq) {
4256 
4257  int size = wrec->size();
4258 
4259  wavearray<double> time(size);
4260  wavearray<double> etime(size); etime=0;
4261  for (int i=0; i<size; i++) time[i] = i/wrec->rate()+(wref->start()-gSEGGPS);
4262 
4263  double bT, eT;
4264  GetTimeBoundaries(*wref, P, bT, eT);
4265  bT-=gSEGGPS;
4266  eT-=gSEGGPS;
4267 
4268  // set to 0 the frequency values outside the time range -> fix the y scale autoscale
4269  // info : this procedure modify the frequency input data but it is not relevant
4270  if(freq) {
4271  for(int i=0;i<wrec->size();i++) {
4272  if(time[i]>bT && time[i]<eT) continue;
4273  wrec->data[i]=0; wmed->data[i]=0; wl50->data[i]=0; wu50->data[i]=0; wl90->data[i]=0; wu90->data[i]=0;
4274  }
4275  }
4276 
4277  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4278 
4279  TGraphAsymmErrors* egr90 = new TGraphAsymmErrors(size,time.data,wmed->data,etime.data,etime.data,wl90->data,wu90->data);
4280  egr90->SetLineColor(17);
4281  egr90->SetFillStyle(1001);
4282  egr90->SetFillColor(17);
4283  egr90->GetXaxis()->SetTitle(xtitle);
4284  if(freq) egr90->GetYaxis()->SetTitle("frequency (hz)");
4285  else egr90->GetYaxis()->SetTitle("magnitude");
4286  egr90->SetTitle(title);
4287  egr90->GetXaxis()->SetTitleFont(42);
4288  egr90->GetXaxis()->SetLabelFont(42);
4289  egr90->GetXaxis()->SetLabelOffset(0.012);
4290  egr90->GetXaxis()->SetTitleOffset(1.5);
4291  egr90->GetYaxis()->SetTitleFont(42);
4292  egr90->GetYaxis()->SetLabelFont(42);
4293  egr90->GetYaxis()->SetLabelOffset(0.01);
4294  egr90->GetYaxis()->SetTitleOffset(1.4);
4295 
4296  TGraphAsymmErrors* egr50 = new TGraphAsymmErrors(size,time.data,wmed->data,etime.data,etime.data,wl50->data,wu50->data);
4297  egr50->SetLineColor(15);
4298  egr50->SetFillStyle(1001);
4299  egr50->SetFillColor(15);
4300 
4301  TGraph* agr = new TGraph(size,time.data,wmed->data);
4302  agr->SetLineWidth(1);
4303  agr->SetLineColor(kWhite);
4304  agr->SetLineStyle(1);
4305 
4306  TGraph* gr = new TGraph(size,time.data,wrec->data);
4307  gr->SetLineWidth(1);
4308  gr->SetLineColor(2);
4309 
4310  TCanvas* canvas = new TCanvas("wrec", "wrec",200,20,800,500);
4311  canvas->cd();
4312  canvas->SetGridx();
4313  canvas->SetGridy();
4314 
4315  egr90->GetXaxis()->SetRangeUser(bT, eT);
4316  egr90->Draw("A3");
4317  egr50->Draw("3same");
4318  agr->Draw("Lsame");
4319  gr->Draw("Lsame");
4320 
4321  if(ofname!="") {
4322  ofname = TString(pdir)+TString("/")+ofname;
4323  canvas->Print(ofname);
4324  cout << "write : " << ofname << endl;
4325  }
4326 
4327  delete canvas;
4328  delete egr50;
4329  delete egr90;
4330  delete agr;
4331  delete gr;
4332 }
4333 
4334 // Dumps reconstructed waveform/time/freq/errors array in ASCII format.
4336 
4337  int nIFO = NET->ifoListSize(); // number of detectors
4338 
4339  char ofname[256];
4340  for(int n=0; n<nIFO; n++) {
4341 
4342  sprintf(ofname,"%s/%s_pe_wave.dat",pdir.Data(),NET->ifoName[n]);
4343 
4344  ofstream out;
4345  out.open(ofname,ios::out);
4346  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
4347  cout << "Create Output File : " << ofname << endl;
4348  out.precision(19);
4349 
4350  // write header
4351  out << "#whitened data : time," <<
4352  " amp_point, amp_mean, amp_rms," <<
4353  " amp_median, amp_lower_50_perc, amp_lower_90_perc, amp_upper_50_perc, amp_upper_90_perc," <<
4354  " frq_point, frq_mean, frq_rms," <<
4355  " frq_median, frq_lower_50_perc, frq_lower_90_perc, frq_upper_50_perc, frq_upper_90_perc" << endl;
4356 
4357  // write data
4358  int size = vREC[n].size();
4359  double dt=1./vREC[n].rate();
4360  //netcluster* pwc = NET->getwc(0);
4361  //detector* pD = NET->getifo(n);
4362  //double toffset = pwc->start+pD->waveForm.start();
4363  for (int i=0; i<size; i++) {
4364  double time = i*dt+vREC[n].start();
4365 
4366  double vl50 = vMED[n][i]-fabs(vL50[n][i]);
4367  double vu50 = vMED[n][i]+fabs(vU50[n][i]);
4368  double vl90 = vMED[n][i]-fabs(vL90[n][i]);
4369  double vu90 = vMED[n][i]+fabs(vU90[n][i]);
4370 
4371  double fl50 = fMED[n][i]-fabs(fL50[n][i]);
4372  double fu50 = fMED[n][i]+fabs(fU50[n][i]);
4373  double fl90 = fMED[n][i]-fabs(fL90[n][i]);
4374  double fu90 = fMED[n][i]+fabs(fU90[n][i]);
4375 
4376  out << time
4377  << " " << vREC[n][i] << " " << vAVR[n][i] << " " << vRMS[n][i]
4378  << " " << vMED[n][i] << " " << vl50 << " " << vl90 << " " << vu50 << " " << vu90
4379  << " " << fREC[n][i] << " " << fAVR[n][i] << " " << fRMS[n][i]
4380  << " " << fMED[n][i] << " " << fl50 << " " << fl90 << " " << fu50 << " " << fu90
4381  << endl;
4382  }
4383 
4384  out.close();
4385  }
4386 }
4387 
4388 // Dumps injected waveform/time array in ASCII format.
4390 
4391  int nIFO = NET->ifoListSize(); // number of detectors
4392 
4393  char ofname[256];
4394  for(int n=0; n<nIFO; n++) {
4395 
4396  sprintf(ofname,"%s/%s_inj_pe.dat",pdir.Data(),NET->ifoName[n]);
4397 
4398  ofstream out;
4399  out.open(ofname,ios::out);
4400  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
4401  cout << "Create Output File : " << ofname << endl;
4402  out.precision(19);
4403 
4404  // write header
4405  out << "# time white_amp" << endl;
4406 
4407  // write data
4408  int size = vINJ[n].size();
4409  double dt=1./vINJ[n].rate();
4410  for (int i=0; i<size; i++) {
4411  double time = i*dt+vINJ[n].start();
4412  out << time << " " << vINJ[n][i] << endl;
4413  }
4414 
4415  out.close();
4416  }
4417 }
4418 
4420 
4421  // import global variables
4422  int gLRETRY=-1; IMPORT(int,gLRETRY)
4423  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
4424  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
4425 
4426  if(gLRETRY==0) return 1;
4427 
4428  bool use_original_data = (gLRETRY==1) ? true : false;
4429  // when multitask only in the last trial the original data is used
4430  if(gOPT.multitask) use_original_data = (gMTRIAL==gOPT.trials) ? true : false;
4431 
4432  if(use_original_data) cout << endl << "Last Trial : Analysis of the original data" << endl << endl;
4433 
4434  int nRES = cfg->l_high-cfg->l_low+1; // number of frequency resolution levels
4435  int nIFO = NET->ifoListSize(); // number of detectors
4436 
4437  vector<TString> delObjList;
4438  // supercluster clusters and parse maps are removed
4439  delObjList.push_back("supercluster");
4440  delObjList.push_back("sparse");
4441  TString jname = jfile->GetPath();
4442  jname.ReplaceAll(":/","");
4443  jfile->Close();
4444  gCWB2G->FileGarbageCollector(jname,"",delObjList);
4445 
4446  // circular shift in the range [jS,jE] randomly whitened ifos HoT & add rec/inj waveforms into @ the same rec time
4447  int size,rate;
4448  int jS,jE,jW,k;
4450  for(int n=0; n<nIFO; n++) { // create random time series
4451 
4452  // select waveform (reconstructed/injected) to be added to the whitened HoT
4454  if(gOPT.signal==0) wf=vREC[n]; // reconstructed
4455  if(gOPT.signal==1) wf=GetAlignedWaveform(&vINJ[n], &vREC[n]); // injected
4456  if(gOPT.signal==2) wf=vDAT[n]; // reconstructed+null
4457 
4458  // apply amplitude mis-calibration amp_cal_err
4459  if((!use_original_data) && gOPT.amp_cal_err[n]) { // in the last try (gLRETRY==1) we use the original data
4460  double amp_cal_err = 1;
4461  if(gOPT.amp_cal_err[n]>0) amp_cal_err = gRandom->Uniform(1-gOPT.amp_cal_err[n],1+gOPT.amp_cal_err[n]);
4462  else amp_cal_err = gRandom->Gaus(1,fabs(gOPT.amp_cal_err[n]));
4463  cout << "Apply Amplitude Mis-Calibration (" << gOPT.amp_cal_err[n] <<"%) "
4464  << NET->ifoName[n] << " -> amp_cal_err : " << amp_cal_err << endl;
4465  wf*=amp_cal_err;
4466  }
4467  // apply phase mis-calibration phs_cal_err
4468  if((!use_original_data) && gOPT.phs_cal_err[n]) { // in the last try (gLRETRY==1) we use the original data
4469  double phs_cal_err = 0;
4470  if(gOPT.phs_cal_err[n]>0) phs_cal_err = gRandom->Uniform(-gOPT.phs_cal_err[n],gOPT.phs_cal_err[n]);
4471  else phs_cal_err = gRandom->Gaus(0,fabs(gOPT.phs_cal_err[n]));
4472  cout << "Apply Phase Mis-Calibration (" << gOPT.phs_cal_err[n] <<" deg) "
4473  << NET->ifoName[n] << " -> phs_cal_err : " << phs_cal_err << " deg" << endl;
4474  CWB::mdc::PhaseShift(wf,phs_cal_err);
4475  }
4476 
4477  pTF = NET->getifo(n)->getTFmap();
4478  size = gHOT[n].size();
4479  rate = gHOT[n].rate();
4480  // apply time shift to input whitened data (integer number of sammples)
4481  jS = cfg->segEdge*rate;
4482  jE = size-jS;
4483  jW = rate*cfg->iwindow/2;
4484  k = gRandom->Uniform(jS+jW,jE-jW); // select random offset (excludes +/- iwindow/2 around the event)
4485  if(use_original_data) { // in the last try we restore the original data
4486  *gCWB2G->hot[n] = gHOT[n];
4487  } else {
4488  printf("Info : %s data time shift : %.5f sec\n",NET->getifo(n)->Name,(double)(k-jS)/rate);
4489  for(int i=jS; i<jE; i++) {
4490  gCWB2G->hot[n]->data[k++] = gHOT[n].data[i];
4491  if(k==jE) k=jS;
4492  }
4493  *(gCWB2G->hot[n]) = AddWaveforms(gCWB2G->hot[n],&wf); // add waveform (reconstructed/injected) to whitened HoT
4494  }
4495  pTF->Forward(*(gCWB2G->hot[n]),*(NET->wdmList[nRES-1]));
4496  }
4497 
4498  // perform coherence and supercluster stages
4499  gCWB2G->Coherence(gIFACTOR);
4500  gCWB2G->SuperCluster(gIFACTOR);
4501 
4502  NET->setDelayIndex(gCWB2G->TDRate);
4503  if(gCWB2G->singleDetector) NET->setIndexMode(cfg->mode); // when nIFO=1 exclude duplicate delay configurations
4504 
4505  // set low-rate TD filters
4506  for(int k=0; k<nRES; k++) gCWB2G->pwdm[k]->setTDFilter(cfg->TDSize, cfg->upTDF);
4507  NET->setDelayIndex(gCWB2G->TDRate);
4508 
4509  jfile = new TFile(jname);
4510  ReplaceSuperclusterData(jfile, cfg, NET, gOPT.gps); // save in jfile only max size supercluster
4511 
4512  // read sparse map from job file
4513  cout << "Loading sparse TF map ... " << endl;
4514  for(int n=0; n<nIFO; n++) {
4515  detector* pD = NET->getifo(n);
4516  pD->sclear(); // clear vector with sparse maps
4517  TString ifo=NET->ifoName[n];
4518  for(int i=0; i<nRES; i++) {
4519  char swname[32];
4520  if(cfg->simulation) sprintf(swname,"sparse/%s-level:%d:%d",ifo.Data(),gIFACTOR,i+cfg->l_low);
4521  else sprintf(swname,"sparse/%s-level:%d",ifo.Data(),i+cfg->l_low);
4522  SSeries<double>* psw = (SSeries<double>*)jfile->Get(swname);
4523  if(psw==NULL) {
4524  cout << "CWB_Plugin_PE - Likelihood : sparse map " << swname
4525  << " not exist in job file" << endl;exit(1);
4526  }
4527  SSeries<double> SS = *psw;
4528  pD->vSS.push_back(SS);
4529  delete psw;
4530  }
4531  cout<<endl;
4532  }
4533 
4534  netcluster* pwc = NET->getwc(0);
4535  pwc->cData.clear();
4536  // read cluster list & metadata netcluster object
4537  int cycle = cfg->simulation ? gIFACTOR : Long_t(NET->wc_List[0].shift);
4538  vector<int> clist = pwc->read(jfile,"supercluster","clusters",0,cycle);
4539  return clist.size();
4540 }
4541 
4543 
4544  wavearray<double> wf = *wf1;
4545 
4546  if(wf1==NULL) return wf;
4547  if(wf1->size()==0) return wf;
4548 
4549  double R=wf1->rate();
4550 
4551  double b_wf1 = wf1->start();
4552  double e_wf1 = wf1->start()+wf1->size()/R;
4553  double b_wf2 = wf2->start();
4554  double e_wf2 = wf2->start()+wf2->size()/R;
4555 
4556  int o_wf1 = b_wf1>b_wf2 ? 0 : int((b_wf2-b_wf1)*R+0.5);
4557  int o_wf2 = b_wf1<b_wf2 ? 0 : int((b_wf1-b_wf2)*R+0.5);
4558 
4559  double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
4560  double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
4561  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
4562 
4563  for(int i=0;i<sizeXCOR;i++) wf[i+o_wf1] += wf2->data[i+o_wf2];
4564 
4565  return wf;
4566 }
4567 
4568 void SetEventWindow(CWB::config* cfg, double gps) {
4569 
4570  if(gps<0) return;
4571 
4572  // dq file list
4573  // {ifo, dqcat_file, dqcat[0/1/2], shift[sec], inverse[false/true], 4columns[true/false]}
4574 
4575  for(int n=0; n<cfg->nIFO; n++) {
4576 
4577  strcpy(cfg->DQF[cfg->nDQF].ifo, cfg->ifo[n]);
4578  sprintf(cfg->DQF[cfg->nDQF].file, "%s/%s_%s.gps_%d",cfg->tmp_dir,cfg->ifo[n],cfg->data_label,int(gps));
4579  cfg->DQF[cfg->nDQF].cat = CWB_CAT2;
4580  cfg->DQF[cfg->nDQF].shift = 0.;
4581  cfg->DQF[cfg->nDQF].invert = false;
4582  cfg->DQF[cfg->nDQF].c4 = true;
4583  cfg->nDQF++;
4584 
4585  cout << cfg->DQF[cfg->nDQF-1].file << endl;
4586 
4587  ofstream out;
4588  out.open(cfg->DQF[cfg->nDQF-1].file,ios::out);
4589  cout << "Write file : " << cfg->DQF[cfg->nDQF-1].file << endl;
4590  if (!out.good()) {cout << "Error Opening File : " << cfg->DQF[cfg->nDQF-1].file << endl;exit(1);}
4591  out.precision(14);
4592  int istart = int(gps)-cfg->iwindow;
4593  int istop = int(gps)+cfg->iwindow;
4594  out << "1 " << istart << " " << istop << " " << 2*cfg->iwindow << endl;
4595  out.close();
4596  }
4597 }
4598 
4600 // compute median sky probability
4601 
4602  skymap skyprob = NET->getifo(0)->tau;
4603  skyprob=0;
4604 
4605  int ntry=0; for(int i=0;i<gOPT.trials;i++) if(wREC[i].size()) ntry++; // number of detected events
4606 
4607  int *index = new int[ntry];
4608  float *value = new float[ntry];
4609 
4610  int L = skyprob.size();
4611  double sum=0;
4612  for(int l=0;l<L;l++) { // for each sky location compute median sky probability
4613 
4614  int k=0; for(int i=0;i<gOPT.trials;i++) if(wREC[i].size()) value[k++] = wSKYPROB[i].get(l);
4615  TMath::Sort(ntry,value,index,false);
4616 
4617  int imed = (ntry*50.)/100.; if(imed>=ntry) imed=ntry-1;
4618 
4619  skyprob.set(l,value[index[imed]]); // save median sky probability
4620 
4621  sum+=value[index[imed]];
4622  }
4623 
4624  // normalize sky probability
4625  if(sum>0) for(int l=0;l<L;l++) skyprob.set(l,skyprob.get(l)/sum);
4626 
4627  delete [] index;
4628  delete [] value;
4629 
4630  return skyprob;
4631 }
4632 
4634 
4635  skymap skyprobcc = NET->getifo(0)->tau;
4636 
4637  // convert skyprob in celestial coordinates
4638  int L = skyprobcc.size();
4639  for(int l=0; l<L; l++) {
4640  double th = skyprob->getTheta(l);
4641  double ph = skyprob->getPhi(l);
4642 
4643  double ra = skyprob->getRA(l);
4644  int k=skyprob->getSkyIndex(th, ra);
4645  skyprobcc.set(k,skyprob->get(l));
4646  }
4647 
4648  // dump skyprobcc to fits file
4649  char fname[1024];
4650  if(skyprobcc.getType()) { // check if it is healpix
4651  sprintf(fname, "%s/mskyprobcc.%s", odir.Data(), "fits");
4652 
4653  // build fits configur info
4654  TString analysis = "2G";
4655  if(NET->MRA) analysis+=":MRA";
4656  if(NET->pattern==0) analysis+=":Packet(0)";
4657  if((NET->pattern!=0 && NET->pattern<0)) analysis+=TString::Format(":Packet(%d)",NET->pattern);
4658  if((NET->pattern!=0 && NET->pattern>0)) analysis+=TString::Format(":Packet(+%d)",NET->pattern);
4659 
4660  char configur[64]="";
4661  char search = NET->tYPe;
4662  if (search=='r') sprintf(configur,"%s un-modeled",analysis.Data());
4663  if (search=='i') sprintf(configur,"%s iota-wave",analysis.Data());
4664  if (search=='p') sprintf(configur,"%s psi-wave",analysis.Data());
4665  if((search=='l')||(search=='s')) sprintf(configur,"%s linear",analysis.Data());
4666  if((search=='c')||(search=='g')) sprintf(configur,"%s circular",analysis.Data());
4667  if((search=='e')||(search=='b')) sprintf(configur,"%s elliptical",analysis.Data());
4668  skyprobcc.Dump2fits(fname,EVT->time[0],configur,const_cast<char*>("PROB"),const_cast<char*>("pix-1"),'C');
4669  }
4670 
4671 }
4672 
4674 
4675  int nIFO = cfg->nIFO;
4676 
4677  char beginString[1024];
4678  sprintf(beginString,"wave_");
4679  char endString[1024];
4680  sprintf(endString,"_job%d.root",runID);
4681  char containString[1024];
4682  sprintf(containString,"%s_trial",cfg->data_label);
4683 
4684  vector<TString> fileList = CWB::Toolbox::getFileListFromDir(cfg->output_dir,endString,beginString,containString,containString);
4685 
4686  wavearray<double>* mtpe_wREC[NIFO_MAX];
4687  for(int i=0;i<nIFO;i++) mtpe_wREC[i] = new wavearray<double>;
4688  skymap* mtpe_skyprob = new skymap(int(0));
4689  float* chirp = new float[6];
4690  float* isnr = new float[nIFO];
4691  float* osnr = new float[nIFO];
4692  float* iosnr = new float[nIFO];
4693  float likelihood;
4694 
4695  char command[1024];
4696  int nfile = fileList.size();
4697  for(int n=0;n<nfile;n++) {
4698  cout << n << " " << fileList[n].Data() << endl;
4699 
4700  TFile* froot = new TFile(fileList[n].Data(),"READ");
4701  if(froot==NULL) {
4702  cout << "CWB_Plugin_PE Error : Failed to open file : " << fileList[n].Data() << endl;
4703  gSystem->Exit(1);
4704  }
4705  TTree* itree = (TTree*)froot->Get("waveburst");
4706  if(itree==NULL) {
4707  cout << "CWB_Plugin_PE Error : Failed to open tree waveburst from file : " << fileList[n].Data() << endl;
4708  gSystem->Exit(1);
4709  }
4710 
4711  for(int i=0;i<nIFO;i++) itree->SetBranchAddress(TString::Format("mtpe_wREC_%d",i).Data(),&mtpe_wREC[i]);
4712  itree->SetBranchAddress("mtpe_skyprob",&mtpe_skyprob);
4713  itree->SetBranchAddress("chirp",chirp);
4714  itree->SetBranchAddress("likelihood",&likelihood);
4715 
4716  itree->GetEntry(0); // read wREC,skyprob objects
4717 
4718  // check if objects are not null
4719  for(int i=0;i<nIFO;i++) {
4720  if(mtpe_wREC[i]==NULL) {
4721  cout << "CWB_Plugin_PE Error : Object wavearray not exist !!! " << endl;
4722  gSystem->Exit(1);
4723  }
4724  }
4725  if(mtpe_skyprob==NULL) {
4726  cout << "CWB_Plugin_PE Error : Object skymap not exist !!! " << endl;
4727  gSystem->Exit(1);
4728  }
4729 
4730  // fill object vectors
4731  wSKYPROB[gOPT.trials-n-1] = *mtpe_skyprob; // restore sky probability
4732 
4733  for(int i=0;i<nIFO;i++)
4734  wREC[gOPT.trials-n-1].push_back(*mtpe_wREC[i]); // restore reconstructed waveform
4735 
4736  for(int i=0;i<nIFO;i++) { // restore SNR
4737  iSNR[i].push_back(isnr[i]);
4738  oSNR[i].push_back(osnr[i]);
4739  ioSNR[i].push_back(iosnr[i]);
4740  }
4741 
4742  for(int j=0; j<6; j++) vCHIRP[j].push_back(chirp[j]); // restore chirp mass parameters
4743 
4744  vLIKELIHOOD.push_back(likelihood); // restore likelihood
4745 
4746  froot->Close();
4747 
4748  // remove temporary root file
4749  sprintf(command,"/bin/rm %s",fileList[n].Data());
4750  cout << command << endl;
4751  gSystem->Exec(command);
4752  }
4753 
4754  for(int i=0;i<nIFO;i++) delete mtpe_wREC[i];
4755  delete mtpe_skyprob;
4756  delete [] chirp;
4757  delete [] isnr;
4758  delete [] osnr;
4759  delete [] iosnr;
4760 }
4761 
4762 void GetNoisePixels(std::vector<double>* aNSE, std::vector<double>* ANSE,
4763  network* NET, CWB::config* cfg, int lag, int id) {
4764 // get noise statistic in TF domain using random event TF patterns
4765 // this is used for Kolmogorov test of the residuals
4766 
4767  size_t nIFO = NET->ifoList.size();
4768 
4769  netcluster* pwc = NET->getwc(lag);
4770  std::vector<int> pI = NET->wdmMRA.getXTalk(pwc, id, false); // buffer for pixel IDs
4771  netpixel* pix = pwc->getPixel(id,pI[0]);
4772  int V = pI.size();
4773  if(!V) return;
4774 
4775  int nRES = NET->wdmListSize(); // number of frequency resolution levels
4776 
4777  double R = NET->getifo(0)->TFmap.rate();
4778 
4779  double rate = gHOT[0].rate();
4780  int size = gHOT[0].size();
4781 
4782  int layers_high = 1<<cfg->l_high;
4783 
4784  // compute the minimum index in the super cluster
4785  int index_min=2*gHOT[0].size();
4786  for(int k=0; k<V; k++) {
4787  netpixel* pix = pwc->getPixel(id,pI[k]);
4788  for(int n=0; n<nIFO; n++) {
4789  int index = pix->data[n].index;
4790  if(index<index_min) index_min=index;
4791  }
4792  }
4793  index_min = index_min-(index_min%layers_high); // set index_min a multiple of layers_high
4794 
4795  // find range where to select noise data
4796  int jS = cfg->segEdge*rate;
4797  int jE = size-jS;
4798  int jW = rate*cfg->iwindow;
4799 
4801  for(int n=0; n<nIFO; n++) {
4802  for(int j=0;j<nRES;j++) { // j=0 -> l_high : j=nRES-1 -> l_low
4803  w.Forward(gHOT[n],*(NET->wdmList[nRES-1-j]));
4804  for(int m=0; m<nTRIALS; m++) {
4805  int index_off = gRandom->Uniform(jS,jE-jW); // select random offset (excludes last iwindow range)
4806  index_off = index_off-(index_off%layers_high); // set index_off a multiple of layers_high
4807  for(int k=0; k<V; k++) {
4808  netpixel* pix = pwc->getPixel(id,pI[k]);
4809  int index = pix->data[n].index;
4810  int ires = int(TMath::Log2(R/pix->rate))-cfg->l_low;
4811  int ilayers = 1<<(ires+cfg->l_low);
4812  index = index_off+(index-index_min); // shift pixels
4813  if(ires==j) {
4814  // get original sparse map amplitudes
4815  //double dd = GetSparseMapData(&vSS[n][ires], true, index);
4816  //double DD = GetSparseMapData(&vSS[n][ires], false , index);
4817  // get the shifted amplitudes
4818  double dd = w.data[index];
4819  double DD = w.data[index+w.maxIndex()+1];
4820  aNSE[n].push_back(dd);
4821  ANSE[n].push_back(DD);
4822  }
4823  }
4824  }
4825  }
4826  }
4827  return;
4828 }
4829 
4831 
4832  if(type!="data" && type!="null" && type!="signal" && type!="injection") {
4833  cout << "CWB_Plugin_PE Error : wrong PlotSpectrogram type, must be data/null/signal" << endl;
4834  exit(1);
4835  }
4836 
4837  int nIFO = NET->ifoListSize(); // number of detectors
4838  char fname[1024];
4839 
4840  for(int n=0; n<nIFO; n++) {
4841  WSeries<double>* pTF = NET->getifo(n)->getTFmap();
4842  if(type=="null") {
4843  wavearray<double> wf = vREC[n];
4844  wf*=-1;
4845  *pTF = AddWaveforms(&gHOT[n],&wf); // sub waveform (reconstructed/injected) to whitened HoT
4846  }
4847  if(type=="signal") {
4848  wavearray<double> wf = vREC[n];
4849  *pTF=0;
4850  *pTF = AddWaveforms(pTF,&wf); // reconstructed signal
4851  }
4852  if(type=="injection") {
4853  wavearray<double> wf = vINJ[n];
4854  *pTF=0;
4855  *pTF = AddWaveforms(pTF,&wf); // injected signal
4856  }
4857 
4858  *pTF*=1./sqrt(1<<cfg->levelR);
4859 
4860  TString xtitle = TString::Format("Time (sec) : GPS OFFSET = %.3f",EVT->gps[n]);
4861 
4862  if(pTF->size()==0) continue;
4863  if(pTF->getLevel()>0) pTF->Inverse();
4864 
4865  int nfact=4;
4866  int nfft=nfact*512;
4867  int noverlap=nfft-10;
4868  double fparm=nfact*6;
4869  int ystart = int((EVT->start[n]-EVT->gps[n]-1)*pTF->rate());
4870  int ystop = int((EVT->stop[n]-EVT->gps[n]+1)*pTF->rate());
4871  ystart-=nfft;
4872  ystop+=nfft;
4873  int ysize=ystop-ystart;
4874  wavearray<double> y;y.resize(ysize);y.rate(pTF->rate());y.start(ystart/pTF->rate());
4875 
4876  // stft use dt=y.rate() to normalize data but whitened data are already normalized by dt
4877  // so before stft data must be divided by 1./sqrt(dt)
4878  for(int i=0;i<(int)y.size();i++) y.data[i]=pTF->data[i+ystart]/sqrt(y.rate());
4879 
4880  CWB::STFT stft(y,nfft,noverlap,"energy","gauss",fparm);
4881 
4882  TCanvas* canvas;
4883  double tstart = nfft/pTF->rate()+ystart/pTF->rate();
4884  double tstop = (ysize-nfft)/pTF->rate()+ystart/pTF->rate();
4885 
4886  tstart+=0.9;tstop-=0.9;
4887  stft.Draw(tstart,tstop,pTF->getlow(),pTF->gethigh(),0,0,1);
4888  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle);
4889  sprintf(fname, "%s/%s_%s_spectrogram_0.png", pdir.Data(), NET->ifoName[n], type.Data());
4890  canvas = stft.GetCanvas();
4891  stft.Print(fname);
4892  canvas->SetLogy(true);
4893  stft.GetHistogram()->GetXaxis()->SetTitle(xtitle);
4894  sprintf(fname, "%s/%s_%s_spectrogram_logy_0.png", pdir.Data(), NET->ifoName[n], type.Data());
4895  stft.Print(fname);
4896 
4897  y.resize(0);
4898  }
4899 }
4900 
4901 std::vector<wavearray<double> > GetWhitenedData(network* NET, CWB::config* cfg) {
4902  // get whitened data (in the vREC time range)
4903 
4904  std::vector<wavearray<double> > xWHT; // temporary stuff
4905 
4906  int nIFO = NET->ifoListSize(); // number of detectors
4907 
4908  for(int n=0; n<nIFO; n++) {
4909  xWHT.push_back(GetAlignedWaveform(&gHOT[n], &vREC[n]));
4910  }
4911  return xWHT;
4912 }
4913 
4914 // the following function replace the supercluster clusters with the max size supercluster
4915 void ReplaceSuperclusterData(TFile*& jfile, CWB::config* cfg, network* NET, double gps) {
4916  // gps=0 -> select cluster with max size
4917  // gps>0 -> select cluster nearest to gps time
4918 
4919  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
4920  cwb2G* gCWB2G; IMPORT(cwb2G*,gCWB2G)
4921 
4922  int cycle = cfg->simulation ? gIFACTOR : Long_t(NET->wc_List[0].shift);
4923 
4924  netcluster* pwc = NET->getwc(0);
4925  pwc->clear();
4926 
4927  // read cluster list & metadata netcluster object
4928  vector<int> clist = pwc->read(jfile,"supercluster","clusters",0,cycle);
4929  //pwc->print();
4930 
4931  if(clist.size()==0) return;
4932 
4933  // find max size supercluster index
4934  int maxk=0;
4935  int npixels=0; // total loaded pixels per lag
4936  int msize=0;
4937  for(int k=0;k<(int)clist.size();k++) { // loop over the cluster list
4938  // read pixels & tdAmp into netcluster pwc
4939  pwc->read(jfile,"supercluster","clusters",-2,cycle,0,clist[k]);
4940  int psize = pwc->size()-npixels;
4941  if(psize>msize) {maxk=k;msize=psize;}
4942  npixels+=psize;
4943  }
4944 
4945  // find min abs(gps-time) supercluster index
4946  int mink=0;
4947  double mtime=1.e20;
4948  wavearray<double> ctime = pwc->get((char*)"time",0,'L',0); // get cluster time
4949  for(int k=0; k<ctime.size(); k++) { // loop over clusters
4950  if(fabs(ctime[k]+pwc->start-gps)<mtime) {mink=k;mtime=fabs(ctime[k]+pwc->start-gps);}
4951  }
4952 
4953  int kindex = gps>0 ? mink : maxk;
4954 
4955  pwc->clear();
4956  // read max supercluster (pixels & tdAmp) into netcluster pwc
4957  pwc->read(jfile,"supercluster","clusters",-1,cycle,0,clist[kindex]);
4958 
4959  // supercluster are removed from jfile
4960  vector<TString> delObjList;
4961  delObjList.push_back("supercluster");
4962  TString jname = jfile->GetPath();
4963  jname.ReplaceAll(":/","");
4964  jfile->Close();
4965  gCWB2G->FileGarbageCollector(jname,"",delObjList);
4966 
4967  // max supercluster is stored in jfile
4968  jfile = new TFile(jname,"UPDATE");
4969  gCWB2G->jfile=jfile;
4970  if(jfile==NULL||!jfile->IsOpen())
4971  {cout << "CWB_Plugin_PE.C - Error opening : " << jname << endl;exit(1);}
4972  pwc->write(jfile,"supercluster","clusters",0,cycle);
4973  pwc->write(jfile,"supercluster","clusters",-1,cycle);
4974 
4975  jfile->Write();
4976 }
monster wdmMRA
list of pixel pointers for MRA
Definition: network.hh:633
std::vector< char * > ifoName
Definition: network.hh:591
std::vector< wavearray< double > > fL50
void sethigh(double f)
Definition: wseries.hh:114
detector * getifo(size_t n)
param: detector index
Definition: network.hh:418
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
size_t rateANA
Definition: cwb.hh:267
gskymap * gSM
static float _sse_abs_ps(__m128 *_a)
Definition: watsse.hh:119
uoptions gOPT
double iwindow
Definition: config.hh:182
skymap wSKYPROB[MAX_TRIALS]
double FittingFactor(wavearray< double > *wf1, wavearray< double > *wf2)
virtual size_t size() const
Definition: wavearray.hh:127
static const double C
Definition: GNGen.cc:10
std::vector< SSeries< double > > dSS[NIFO_MAX]
#define NIFO
Definition: wat.hh:56
size_t write(const char *file, int app=0)
Definition: netcluster.cc:2989
wavearray< double > GetWaveformEnvelope(wavearray< double > *wf)
TCut chirp
#define PE_CED_TFMAP
int noverlap
Definition: TestDelta.C:20
size_t nLag
Definition: network.hh:555
void GetNullPixels(std::vector< double > *aNUL, std::vector< double > *ANUL, network *NET, CWB::config *cfg, int lag, int id)
Float_t * rho
biased null statistics
Definition: netevent.hh:93
#define PLOT_TYPE
std::vector< wavearray< double > > fINJ
std::vector< wavearray< double > > vREC
char xtitle[1024]
TString ofName
Float_t * high
min frequency
Definition: netevent.hh:85
static vector< TString > getFileListFromDir(TString dir_name, TString endString="", TString beginString="", TString containString="", bool fast=false)
Definition: Toolbox.cc:4333
#define PE_GPS
size_t TDSize
Definition: config.hh:199
std::vector< wavearray< double > > GetSigWaveform(network *NET, CWB::config *cfg, int lag, int id)
void FileGarbageCollector(TString ifName, TString ofName="", vector< TString > delObjList=vector< TString >())
Definition: cwb.cc:2319
std::vector< wavearray< double > > wREC[MAX_TRIALS]
std::vector< double > vCHIRP[6]
std::vector< netcluster > wc_List
Definition: network.hh:592
par[0] value
void setSLags(float *slag)
Definition: netevent.cc:404
void PlotWaveforms(network *NET, CWB::config *cfg, int ID, TString pdir="")
int wdm_start
Definition: sseries.hh:178
void SetSkyMask(network *NET, CWB::config *cfg, double theta, double phi, double radius)
TString gOUTPUT
printf("total live time: non-zero lags = %10.1f \n", liveTot)
void GetFactorsStatistic(int nIFO)
bool singleDetector
used for the stage stuff
Definition: cwb.hh:265
void DumpRecWavePE(network *NET, TString pdir="")
Double_t * start
cluster duration = stopW-startW
Definition: netevent.hh:80
char cmd[1024]
WSeries< double > * pTF[NIFO_MAX]
pointers to detectors
Definition: cwb.hh:252
std::vector< double > * getmdcTime()
Definition: network.hh:404
WDM< double > * pwdm[NRES_MAX]
Definition: cwb2G.hh:60
virtual void rate(double r)
Definition: wavearray.hh:123
void PlotChirpMass(int gtype, TString pdir, int sim=true)
std::vector< wavearray< float > > tdAmp
Definition: netpixel.hh:105
Float_t * low
average center_of_snr frequency
Definition: netevent.hh:84
double cedRHO
Definition: config.hh:280
void Print(TString pname)
Definition: STFT.cc:297
char skyMaskFile[1024]
Definition: config.hh:290
size_t upTDF
Definition: config.hh:198
wavearray< double > a(hp.size())
std::vector< wavearray< double > > fU50
float likelihood
Definition: netpixel.hh:96
std::vector< wavearray< double > > fREC
std::vector< wavearray< double > > vDAT
wavearray< double > GetCutWaveform(wavearray< double > x, double bT, double eT, double bF, double eF)
wavearray< double > AddWaveforms(wavearray< double > *wf1, wavearray< double > *wf2)
#define EFRACTION_CUT
int n
Definition: cwb_net.C:10
void ResetUserOptions()
wavearray< double > gHOT[NIFO_MAX]
static void _sse_zero_ps(__m128 *_p)
Definition: watsse.hh:26
void PlotFactors(int gtype, int nIFO, TString pdir)
double dT
Definition: cwb.hh:271
bool invert
Definition: Toolbox.hh:70
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
wavearray< float > a_90
buffer for cluster sky 00 amplitude
Definition: network.hh:635
TTree * setTree()
Definition: netevent.cc:412
char * watversion(char c='s')
Definition: watversion.hh:21
TString("c")
#define PE_OUTPUT_MED
size_t setIndexMode(size_t=0)
Definition: network.cc:8072
void SetEventWindow(CWB::config *cfg, double gps)
char ifo[32]
Definition: Toolbox.hh:66
#define PE_CED_DISTR
TFile * jfile
output root file
Definition: cwb.hh:247
double Tb
Definition: cwb.hh:269
std::vector< wavearray< double > > vRMS
int ID
Definition: TestMDC.C:70
Int_t * size
cluster volume
Definition: netevent.hh:61
double gSEGGPS
void ResetPCA(network *NET, CWB::config *cfg, netcluster *pwc, std::vector< netpixel > *vPIX, int ID)
size_t maxIndex()
Definition: wseries.hh:131
#define PE_CED_PCA
#define PE_CED_NULL
std::vector< double > aNUL[NIFO_MAX]
static void PhaseShift(wavearray< double > &x, double pShift=0.)
Definition: mdc.cc:2926
void DumpSkyProb(skymap *skyprob, network *NET, netevent *&EVT, TString odir)
#define PE_OUTPUT_P90
#define PE_OUTPUT_DAT
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<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t 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
std::vector< pixdata > data
Definition: netpixel.hh:104
double phase
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2188
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:364
std::vector< wavearray< double > > GetInjWaveform(network *NET, netevent *EVT, int id, double factor)
bool c4
Definition: Toolbox.hh:71
int _sse_mra_ps(float *, float *, float, int)
Definition: network.hh:1258
double fLow
Definition: config.hh:122
void PlotWaveform(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wf1, wavearray< double > *wf2, wavearray< double > *wf3, wavearray< double > *wref, bool fft=false, TString pdir="", double P=0.99, EColor col1=kRed, EColor col2=kBlue, EColor col3=(EColor) 418)
char odir[1024]
float theta
int nfact
Definition: TestDelta.C:18
int wdm_nSTS
Definition: sseries.hh:179
double gIPHI
int pattern
Definition: config.hh:223
netpixel pix(nifo)
delete[] radius
netcluster * pwc
Definition: cwb_job_obj.C:20
int nfft
Definition: TestDelta.C:19
#define PE_OUTPUT_P50
TH2F * ph
ofstream finj
int GetLayer(int index)
Definition: sseries.hh:104
std::vector< wavearray< double > > vMED
std::vector< TGraph * > graph
Definition: watplot.hh:176
return wmap canvas
CWB_STAGE istage
Definition: cwb.hh:235
double gOPHI
#define PE_CED_DUMP
bool cedDump
Definition: config.hh:279
void setlow(double f)
Definition: wseries.hh:107
TRandom3 P
Definition: compare_bkg.C:296
int slagID
Definition: cwb.hh:279
std::vector< SSeries< double > > vSS[NIFO_MAX]
std::vector< wavearray< double > > cINJ
char ifostr[64]
dqfile DQF[DQF_MAX]
Definition: config.hh:331
std::vector< wavearray< double > > vINJ
void Coherence(int ifactor)
Definition: cwb2G.cc:755
double getTheta(size_t i)
Definition: skymap.hh:206
Float_t * ioSNR
reconstructed snr waveform
Definition: netevent.hh:120
std::vector< wavearray< double > > vNUL
static wavearray< double > getHilbertIFrequency(wavearray< double > x)
Definition: Toolbox.cc:6151
#define PE_CED_SKYMAP
waveform wf
TString DumpCED(network *NET, netevent *&EVT, CWB::config *cfg, double factor)
std::vector< double > fRES
size_t setcore(bool core, int id=0)
Definition: netcluster.cc:224
Long_t size
WSeries< double > waveBand
Definition: detector.hh:338
void SetOptions(int simulation, double rho, double inRate, bool useSparse=false, char *gtype=const_cast< char * >("png"), int paletteId=0)
Definition: ced.hh:68
double gITHETA
#define SKYMASK_RADIUS
#define M
Definition: UniqSLagsList.C:3
std::vector< wavearray< double > > fL90
int m
Definition: cwb_net.C:10
#define PE_INDEX
std::vector< double > oSNR[NIFO_MAX]
DataType_t * pWWS
Definition: WaveDWT.hh:123
int wdm_rate
Definition: sseries.hh:177
virtual void start(double s)
Definition: wavearray.hh:119
int j
Definition: cwb_net.C:10
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
Definition: network.cc:3635
i drho i
std::vector< double > mdcTime
Definition: network.hh:596
void setTDFilter(int nCoeffs, int L=1)
Definition: WDM.cc:621
char nodedir[1024]
Definition: config.hh:334
skymap tau
Definition: detector.hh:328
search
std::vector< double > ANUL[NIFO_MAX]
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:132
void Draw(int dpaletteId=0, Option_t *option="colfz")
Definition: gskymap.cc:442
wavearray< int > sparseIndex
Definition: sseries.hh:162
std::vector< SSeries< double > > jSS[NIFO_MAX]
std::vector< detector * > ifoList
Definition: network.hh:590
TH1D gHSKYPROB
std::vector< double > ioSNR[NIFO_MAX]
#define N
double gBF
bool outPlugin
Definition: config.hh:351
#define PE_CED_CM
tuple ff
Definition: cwb_online.py:394
size_t wdmListSize()
Definition: network.hh:428
bool core
Definition: netpixel.hh:102
#define PE_OUTPUT_RMS
char ifo[NIFO_MAX][8]
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:291
Definition: cwb2G.hh:15
TTree * gTREE
TH2F * hist2D
Definition: watplot.hh:175
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
Definition: STFT.cc:76
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
size_t ifoListSize()
Definition: network.hh:413
static float _avx_loadata_ps(float **p, float **q, float **u, float **v, float En, std::vector< float * > &pAVX, int I)
Definition: watavx.hh:412
bool online
Definition: config.hh:100
std::vector< wavearray< double > > vL90
iseg push_back(seg)
#define PI
Definition: watfun.hh:14
void clean(int cID=0)
Definition: netcluster.hh:433
int Psave
Definition: config.hh:175
void clear()
Definition: watplot.hh:40
TChain sim("waveburst")
wavearray< double > w
Definition: Test1.C:27
double segEdge
Definition: config.hh:146
#define nIFO
void wrate(double r)
Definition: wseries.hh:102
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:702
std::vector< wavearray< double > > GetRecWaveform(network *NET, netevent *EVT, int id)
ofstream out
Definition: cwb_merge.C:196
CWB_CAT cat
Definition: Toolbox.hh:68
wavearray< float > pNRG
buffers for cluster residual energy
Definition: network.hh:637
skymap mSKYPROB
size_t mode
Definition: config.hh:257
double tstart
TCanvas * canvas
Definition: watplot.hh:174
int getLevel()
Definition: wseries.hh:91
float phi
#define PE_OUTPUT_AVR
double GetSegEnd()
Definition: cwb.hh:171
gWSeries< double > gw(w)
void PlotWaveformAsymmErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wrec, wavearray< double > *wmed, wavearray< double > *wl50, wavearray< double > *wu50, wavearray< double > *wl90, wavearray< double > *wu90, wavearray< double > *wref, TString pdir, double P, bool freq=false)
double factor
double ra
Definition: ConvertGWGC.C:46
std::vector< wavearray< double > > vU50
char file[1024]
Definition: Toolbox.hh:67
x plot
double time[6]
Definition: cbc_plots.C:435
jfile
Definition: cwb_job_obj.C:25
double getlow() const
Definition: wseries.hh:111
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:687
size_t lagOff
Definition: config.hh:152
Int_t Psave
number of detectors
Definition: netevent.hh:48
TGraph * gr
int simulation
Definition: config.hh:181
int l_high
Definition: config.hh:138
int nDQF
Definition: config.hh:330
#define PE_SKYMASK
std::vector< double > wRMS[NIFO_MAX]
void PlotSparse(int ifoID, network *NET, CWB::config *cfg, int ID, wavearray< double > *wf)
#define PE_CED_rINJ
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:51
#define PE_CED_PSD
watplot * WTS
Definition: ChirpMass.C:115
xtalk getXTalk(int nLay1, size_t indx1, int nLay2, size_t indx2)
param: numbers of layers identifying the resolution of the first pixel param: TF map index of the fir...
Definition: monster.cc:351
Double_t * gps
average center_of_gravity time
Definition: netevent.hh:76
double getRA(size_t i)
Definition: skymap.hh:197
void SuperCluster(int ifactor)
Definition: cwb2G.cc:939
wavearray< int > sparseLookup
Definition: sseries.hh:160
i() int(T_cor *100))
#define nTRIALS
double GetTimeBoundaries(wavearray< double > x, double P, double &bT, double &eT)
char output_dir[1024]
Definition: config.hh:300
double TDRate
Definition: cwb2G.hh:58
network NET
Definition: cwb_dump_inj.C:12
double Pi
void setDelayIndex(double rate)
param: MDC log file
Definition: network.cc:2865
const int NIFO_MAX
Definition: wat.hh:4
double gEF
void PlotFrequencyErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *frec, wavearray< double > *favr, wavearray< double > *ferr, wavearray< double > *wref, TString pdir, double P)
int BATCH
Definition: config.hh:227
TH2D * GetHistogram()
Definition: STFT.hh:53
std::vector< netpixel > DoPCA(network *NET, CWB::config *cfg, int lag, int id)
char tmp_dir[1024]
Definition: config.hh:307
Double_t * stop
GPS start time of the cluster.
Definition: netevent.hh:81
double GetSegBegin()
Definition: cwb.hh:170
char parPlugin[1024]
Definition: config.hh:345
wavearray< float > a_00
wdm multi-resolution analysis
Definition: network.hh:634
std::vector< WDM< double > * > wdmList
Definition: network.hh:599
WSeries< double > pTF[nRES]
Definition: revMonster.cc:8
int gIFACTOR
double start
Definition: netcluster.hh:361
static void _avx_free_ps(std::vector< float * > &v)
Definition: watavx.hh:21
TIter next(twave->GetListOfBranches())
char fname[1024]
std::vector< wavearray< double > > fRMS
double shift
Definition: Toolbox.hh:69
#define PE_OUTPUT_INJ
void GetNoisePixels(std::vector< double > *aNSE, std::vector< double > *ANSE, network *NET, CWB::config *cfg, int lag, int id)
std::vector< int > RWFID
Definition: detector.hh:363
int Write(double factor, bool fullCED=true)
Definition: ced.cc:577
std::vector< wavearray< double > > fAVR
void SetOutputFile(network *NET, netevent *&EVT, CWB::config *cfg, bool dump_pe)
static float _avx_setAMP_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
Definition: watavx.hh:734
size_t mdc__IDSize()
Definition: network.hh:396
char channelNamesRaw[NIFO_MAX][50]
Definition: config.hh:292
void SetChannelName(char *chName)
Definition: ced.hh:79
std::vector< double > iSNR[NIFO_MAX]
void ClearVectors()
bool gCEDDUMP
char output[256]
int k
int pattern
Definition: network.hh:583
double tstop
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:94
void ClearWaveforms(detector *ifo)
int nfactor
Definition: config.hh:183
std::vector< wavearray< double > * > IWFP
Definition: detector.hh:361
void PlotSNRnet(int nIFO, TString pdir, int sim=true)
std::vector< double > wAVR[NIFO_MAX]
void DumpOutputFile(network *NET, netevent *&EVT, CWB::config *cfg, int ID, int k, int factor)
#define PE_CED_INJ
float chi2
Definition: cbc_plots.C:603
double gINRATE
void PlotResiduals(int ID, TString pdir="", int sim=true, char type='w')
static double A
Definition: geodesics.cc:8
TObjArray * token
double acor
Definition: network.hh:567
#define PE_SIGNAL
size_t lagSize
Definition: config.hh:150
double F
#define PE_CED_RDR
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
Definition: netcluster.cc:3298
void DumpInjWavePE(network *NET, TString pdir="")
Definition: skymap.hh:45
#define PE_MULTITASK
string command
Definition: cwb_online.py:382
std::vector< int > IWFID
Definition: detector.hh:360
#define PE_TRIALS
#define PE_AMP_CAL_ERR
PE gPE
double GetFrequencyBoundaries(wavearray< double > x, double P, double &bF, double &eF)
size_t size()
Definition: netcluster.hh:129
TString ofileName
Definition: MergeTrees.C:37
skymap GetMedianSkyProb(network *NET)
char tag[256]
Definition: cwb_merge.C:74
TFile * froot
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:395
void WriteBodyHTML(TString html_template, TString html_tag_beg, TString html_tag_end, ofstream *out, int nIFO=1, TString *ifo=NULL)
cwb gCWB
double gethigh() const
Definition: wseries.hh:118
double dt
virtual void FFTW(int=1)
Definition: wavearray.cc:878
std::vector< wavearray< double > > vWHT
std::vector< wavearray< double > > vAVR
void PrintUserOptions(CWB::config *cfg)
int l_low
Definition: config.hh:137
char options[256]
#define EXPORT(TYPE, VAR, CMD)
Definition: config.cc:74
double getPhi(size_t i)
Definition: skymap.hh:146
WSeries< double > TFmap
Definition: detector.hh:336
double GetCentralFrequency(wavearray< double > x)
std::vector< double > vLIKELIHOOD
#define PE_SEED
skymap gSKYPROB
char title[256]
Definition: SSeriesExample.C:1
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:68
bool wfsave
Definition: network.hh:582
double gps
#define PE_CED_REC
void sclear()
Definition: detector.hh:173
netevent EVT(itree, nifo)
double GetCentralTime(wavearray< double > x)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:69
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:75
double T
Definition: testWDM_4.C:11
double gET
void Expand(bool bcore=true)
Definition: sseries.cc:397
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
TLine * line
Definition: compare_bkg.C:482
void PlotWaveformErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wrec, wavearray< double > *wavr, wavearray< double > *werr, wavearray< double > *wref, TString pdir="", double P=0.99)
ifstream in
std::vector< int > sCuts
Definition: netcluster.hh:374
WSeries< double > waveForm
Definition: detector.hh:337
void Wave2Sparse(network *NET, CWB::config *cfg, char type)
void dclose()
Definition: netevent.hh:300
double xtime
Definition: cwb_dump_sjob.C:10
wavearray< double > ts(N)
double fHigh
Definition: config.hh:123
void GetChirpMassStatistic(std::vector< double > *vCHIRP)
std::vector< wavearray< double > > GetWhitenedData(network *NET, CWB::config *cfg)
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:161
cout<< "Starting reading output directory ..."<< endl;vector< TString > fileList
int getType()
Definition: skymap.hh:278
wavearray< int > index
wavearray< float > _avx_norm_ps(float **, float **, std::vector< float * > &, int)
Definition: network.hh:1347
char Name[16]
Definition: detector.hh:309
char ifo[NIFO_MAX][8]
Definition: config.hh:106
double fabs(const Complex &x)
Definition: numpy.cc:37
#define MAX_TRIALS
double GetSparseMapData(SSeries< double > *SS, bool phase, int index)
void GetNullStatistic(std::vector< double > *vNUL, std::vector< double > *vNSE, int ifoId, TString ifoName, TString gtype, TString pdir="")
string file
Definition: cwb_online.py:385
wavearray< float > sparseMap00
Definition: sseries.hh:163
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
Definition: cwb.cc:2465
wavearray< float > sparseMap90
Definition: sseries.hh:164
int gRATEANA
std::vector< wavearray< double > > vL50
TCanvas * GetCanvas()
Definition: STFT.hh:52
std::vector< int > ComputeStatisticalErrors(network *NET, CWB::config *cfg, int ID)
TString cwb_doc_url
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:228
strcpy(RunLabel, RUN_LABEL)
int cnt
Float_t likelihood
Definition: netevent.hh:105
Meyer< double > S(1024, 2)
double df
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:421
void DumpUserOptions(TString odir, CWB::config *cfg)
int l
Definition: cbc_plots.C:434
std::vector< SSeries< double > > rSS[NIFO_MAX]
void PlotSpectrogram(TString type, network *NET, netevent *&EVT, CWB::config *cfg, TString pdir)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void CreateIndexHTML(TString dirCED, int nIFO, TString *ifo, bool sim=false)
float erR[11]
std::vector< clusterdata > cData
Definition: netcluster.hh:373
wavearray< double > GetDifWaveform(wavearray< double > *wf1, wavearray< double > *wf2)
TTree * itree
void SaveSkyProb(network *NET, CWB::config *cfg, int id)
TString gfile
void set(size_t i, double a)
param: sky index param: value to set
Definition: skymap.hh:104
std::vector< double > aNSE[NIFO_MAX]
std::vector< wavearray< double > > vPCA
int nIFO
Definition: config.hh:102
std::vector< wavearray< double > > GetPCAWaveform(network *NET, CWB::config *cfg, int lag, int id)
watplot * GetWATPLOT()
Definition: gwavearray.hh:70
int RedoAnalysis(TFile *jfile, CWB::config *cfg, network *NET)
Definition: cwb.hh:119
DataType_t * data
Definition: wavearray.hh:301
wavearray< float > rNRG
buffer for cluster sky 90 amplitudes
Definition: network.hh:636
const int nRES
Definition: revMonster.cc:7
TString jname
Long_t id
Definition: ced.hh:26
double factors[FACTORS_MAX]
Definition: config.hh:184
double get(size_t i)
param: sky index
Definition: skymap.cc:681
#define PE_NOISE
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:438
skymap GetSkyProb(network *NET, int id)
double dT
Definition: testWDM_5.C:12
#define NOISE_SIGMA
wavearray< double > * hot[NIFO_MAX]
wavelet pointers: pwdm[0] - l_high, wdm[nRES-1] l_low
Definition: cwb2G.hh:61
wavearray< double > wINJ[NIFO_MAX]
void ReplaceSuperclusterData(TFile *&jfile, CWB::config *cfg, network *NET, double gps=0)
SSeries< double > * getSTFmap(size_t n)
Definition: detector.hh:169
double ctime
wavearray< double > GetAlignedWaveform(wavearray< double > *wf1, wavearray< double > *wf2)
string version
Definition: cWB_conf.py:136
char data_label[1024]
Definition: config.hh:314
Bool_t fill_in(network *, int, bool=true)
Definition: injection.cc:325
list files
Definition: cwb_online.py:529
float rate
Definition: netpixel.hh:95
std::vector< SSeries< double > > vSS
Definition: detector.hh:334
int runID
Definition: cwb.hh:234
std::vector< wavearray< double > > vU90
int gMTRIAL
size_t read(const char *)
Definition: netcluster.cc:3096
size_t size()
Definition: skymap.hh:118
WSeries< double > waveNull
Definition: detector.hh:339
#define PE_PHS_CAL_ERR
size_t getmdc__ID(size_t n)
Definition: network.hh:408
std::vector< wavearray< double > > GetWaveform(network *NET, int lag, int id, char type, bool shift=true)
double shift[NIFO_MAX]
void AddNoiseAndCalErrToSparse(network *NET, CWB::config *cfg, char type)
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
Definition: injection.hh:82
virtual void resize(unsigned int)
Definition: wavearray.cc:445
double gOTHETA
void Print(TString pname)
Definition: gskymap.cc:1104
bool MRA
Definition: network.hh:581
int check
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:273
wavearray< double > y
Definition: Test10.C:31
int gID
int getSTFind(double r)
Definition: detector.hh:164
static float _avx_packet_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
Definition: watavx.hh:532
char tYPe
Definition: network.hh:570
std::vector< wavearray< double > > fU90
size_t inRate
Definition: config.hh:114
void PlotTimeFrequencyPCA(network *NET, netevent *EVT, CWB::config *cfg, int id, int lag, TString pdir)
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:66
double fparm
Definition: TestDelta.C:22
std::vector< double > ANSE[NIFO_MAX]
#define PE_OUTPUT_REC
Float_t * chirp
range to source: [0/1]-rec/inj
Definition: netevent.hh:109
int levelR
Definition: config.hh:134
void GetSNRnetStatistic(int nIFO)
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
Definition: gwavearray.cc:71
detector ** pD
void ReadUserOptions(TString options)
std::vector< double > vRES
#define PE_OUTPUT_WHT
double gBT
void LoadFromMultiTaskJobOutput(int runID, CWB::config *cfg)
CWB::STFT * stft
Definition: ChirpMass.C:117
std::vector< wavearray< double > > fMED
void SetZaxisTitle(TString zAxisTitle)
Definition: gskymap.hh:132
int binarySearch(int array[], int start, int end, int key)
Definition: sseries.cc:443
std::vector< vector_int > nTofF
Definition: netcluster.hh:386
exit(0)
bool setdata(double a, char type='R', size_t n=0)
Definition: netpixel.hh:40
Float_t * iSNR
energy of reconstructed responses Sk*Sk
Definition: netevent.hh:118
void SetWaveformCuts(wavearray< double > *x, double bT, double eT, double bF, double eF)
double avr
void clear()
Definition: netcluster.hh:409
Float_t * oSNR
injected snr waveform
Definition: netevent.hh:119
bool dump
Definition: config.hh:277