11 #include "TObjString.h"
12 #include "TObjArray.h"
26 #include "TMultiGraph.h"
31 #include "Math/VectorUtil.h"
61 #define LAL_USE_OLD_COMPLEX_STRUCTS
62 #define restrict // this is to getrid of restrict which is defined only in c99
64 #include <lal/LIGOLwXMLInspiralRead.h>
65 #include <lal/LALConfig.h>
66 #include <lal/LALStdio.h>
67 #include <lal/LALStdlib.h>
68 #include <lal/LALError.h>
69 #include <lal/LALDatatypes.h>
70 #include <lal/LIGOMetadataUtils.h>
71 #include <lal/LIGOMetadataTables.h>
72 #include <lal/AVFactories.h>
73 #include <lal/NRWaveIO.h>
74 #include <lal/NRWaveInject.h>
75 #include <lal/Inject.h>
76 #include <lal/FileIO.h>
77 #include <lal/Units.h>
78 #include <lal/FrequencySeries.h>
79 #include <lal/TimeSeries.h>
80 #include <lal/TimeFreqFFT.h>
81 #include <lal/VectorOps.h>
82 #include <lal/LALDetectors.h>
83 #include <lal/FindChirp.h>
84 #include <lal/Random.h>
85 #include <lal/LALNoiseModels.h>
87 #include <lal/LALStatusMacros.h>
88 #include <lal/LALSimulation.h>
89 #include <lal/LALInspiral.h>
90 #include <lal/LALSimInspiralWaveformFlags.h>
91 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \
92 (LAL_VERSION_MINOR > 14 || (LAL_VERSION_MINOR == 14 && \
93 LAL_VERSION_MICRO >= 0 ))) // LAL_VERSION >= 6.14.0
94 #include <lal/SnglBurstUtils.h>
96 #include <lal/LIGOMetadataBurstUtils.h>
98 #include <lal/LIGOLwXMLBurstRead.h>
99 #include <lal/GenerateBurst.h>
100 #include <lal/LIGOLwXML.h>
109 #define MDC_INJ_RATE 0.01 // sec^-1
110 #define MDC_INJ_JITTER 10. // sec
111 #define MDC_INJ_LENGTH 1. // sec
112 #define MDC_INJ_HRSS 2.5e-21
114 #define MDC_SAMPLE_RATE 16384. // sample/sec
116 #define EPZOOM 0.9999999
167 default:
return "MDC_RANDOM";
216 class mdc :
public TNamed {
234 double srate, vector<mdcpar> par=vector<mdcpar>());
236 double srate, vector<mdcpar> par=vector<mdcpar>());
268 double& iota,
double&
hrss,
int&
ID,
int&
id);
270 double& iota,
double&
hrss,
int&
ID,
int&
id);
333 void Print(
int level=0);
349 size_t length=1000,
bool log=
false, vector<TString>
chName=vector<TString>());
379 TString GenInspiralXML(
int gps_start_time,
int gps_end_time,
bool rmFile=
false);
386 void CreateBurstXML(
TString fName, vector<mdcpar> xml_parms);
388 void CloseBurstXML();
391 static double e2cosi(
double e);
392 static double cosi2e(
double cosi);
406 void exit(
int err) {gSystem->Exit(err);}
455 ProcessTable *xml_process_table_head;
456 ProcessTable *xml_process;
457 ProcessParamsTable *xml_process_params_table_head;
458 SearchSummaryTable *xml_search_summary_table_head;
459 SearchSummaryTable *xml_search_summary;
460 TimeSlide *xml_time_slide_table_head;
461 SimBurst *xml_sim_burst_table_head;
462 SimBurst **xml_sim_burst;
void Dump(TString fname, int ID, int id, TString polarization)
waveform GetSourceWaveform(int &ID, int &id)
vector< mdcpar > GetSkyParms()
watplot * Draw(TString name, int id=0, TString polarization="hp", MDC_DRAW type=MDC_TIME, TString options="ALP", Color_t color=kBlack)
vector< mdcpar > sky_parms
static double GetTimeRange(wavearray< double > x, double &tMin, double &tMax, double efraction=EPZOOM)
std::vector< float > psiList
TString Get(wavearray< double > &x, TString ifo)
std::vector< std::string > mdcList
void SetInjRate(double inj_rate=MDC_INJ_RATE)
wavearray< double > a(hp.size())
static double GetCentralTime(wavearray< double > x)
void SetSourceListSeed(unsigned int srcList_seed)
std::vector< std::string > mdcType
wavearray< double > GetSGQ(double frequency, double Q)
static void PhaseShift(wavearray< double > &x, double pShift=0.)
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
unsigned int srcList_seed
void GetSourceCoordinates(double &theta, double &phi, double &psi, double &rho, double &iota, double &hrss, int &ID, int &id)
std::vector< double > gpsList
void SetInjLength(double inj_length=MDC_INJ_LENGTH)
static void AddExp(wavearray< double > &td, double v, int M)
wavearray< double > GetGA(double duration)
std::vector< float > phList
std::vector< int > IDList
vector< waveform > GetWaveformList()
double GetDelay(TString ifo1, TString ifo2, double phi, double theta)
watplot * DrawTime(wavearray< double > &x, TString options="ALP", Color_t color=kBlack)
mdc & operator=(const mdc &)
mdcid AddWaveform(MDC_TYPE mdc_type, vector< mdcpar > par, TString uname="")
virtual void Browse(TBrowser *b)
std::vector< std::string > nameList
void SetSkyDistribution(MDC_DISTRIBUTION sky_distribution, vector< mdcpar > par, int seed=0, bool add=false)
TString GetTemporaryFileName(TString tag="mdc", TString ext="txt", TString dir="/tmp", bool mkdir=false)
double GetPar(TString name, vector< mdcpar > par, bool &error)
int GetWaveformID(TString name)
std::vector< std::string > xmlType
void DrawSkyDistribution(TString name="skymap", TString projection="", TString coordinate="Geographic", double resolution=2, bool background=true)
wavearray< double > GetCGQ(double frequency, double Q)
static void AddGauss(wavearray< double > &td, double v, double u=0.)
void SetInjOffset(double inj_offset=0)
TString WriteFrameFile(TString frDir, TString frLabel, size_t gps, size_t length=1000, bool log=false, vector< TString > chName=vector< TString >())
void DumpLogHeader(TString fName, TString label="", int size=0)
cout<< "Injected signals: "<< mdc.GetEntries()<< endl;cout<< "Injected signals in histogram factor_events_inj: "<< NEVTS<< endl;float myifar, ecor, m1, m2, netcc[3], neted, penalty;float rho[2];float chirp[6];float range[2];float frequency[2];float iSNR[3], sSNR[3];sim.SetBranchAddress("mass", mass);sim.SetBranchAddress("factor",&factor);sim.SetBranchAddress("range", range);sim.SetBranchAddress("chirp", chirp);sim.SetBranchAddress("rho", rho);sim.SetBranchAddress("netcc", netcc);sim.SetBranchAddress("neted",&neted);sim.SetBranchAddress("ifar",&myifar);sim.SetBranchAddress("ecor",&ecor);sim.SetBranchAddress("penalty",&penalty);sim.SetBranchAddress("time", mytime);sim.SetBranchAddress("iSNR", iSNR);sim.SetBranchAddress("sSNR", sSNR);sim.SetBranchAddress("spin", spin);sim.SetBranchAddress("frequency", frequency);float **volume=new float *[NBINS_mass1];float **volume_first_shell=new float *[NBINS_mass1];float **radius=new float *[NBINS_mass1];float **error_volume=new float *[NBINS_mass1];float **error_volume_first_shell=new float *[NBINS_mass1];float **error_radius=new float *[NBINS_mass1];for(int i=0;i< NBINS_mass1;i++){volume[i]=new float[NBINS_mass2];volume_first_shell[i]=new float[NBINS_mass2];radius[i]=new float[NBINS_mass2];error_volume[i]=new float[NBINS_mass2];error_volume_first_shell[i]=new float[NBINS_mass2];error_radius[i]=new float[NBINS_mass2];for(int j=0;j< NBINS_mass2;j++){volume[i][j]=0.;volume_first_shell[i][j]=0.;radius[i][j]=0.;error_volume[i][j]=0.;error_volume_first_shell[i][j]=0.;error_radius[i][j]=0.;}}float **spin_mtot_volume=new float *[NBINS_MTOT+1];float **spin_mtot_radius=new float *[NBINS_MTOT+1];float **error_spin_mtot_volume=new float *[NBINS_MTOT+1];float **error_spin_mtot_radius=new float *[NBINS_MTOT+1];for(int i=0;i< NBINS_MTOT+1;i++){spin_mtot_volume[i]=new float[NBINS_SPIN+1];spin_mtot_radius[i]=new float[NBINS_SPIN+1];error_spin_mtot_volume[i]=new float[NBINS_SPIN+1];error_spin_mtot_radius[i]=new float[NBINS_SPIN+1];for(int j=0;j< NBINS_SPIN+1;j++){spin_mtot_volume[i][j]=0.;error_spin_mtot_volume[i][j]=0.;spin_mtot_radius[i][j]=0.;error_spin_mtot_radius[i][j]=0.;}}char fname[1024];sprintf(fname,"%s/recovered_signals.txt", netdir);ofstream fev;fev.open(fname, std::ofstream::out);sprintf(line,"#GPS@L1 FAR[Hz] eFAR[Hz] Pval ""ePval factor rho frequency iSNR sSNR \n");fev<< line<< endl;ofstream *fev_single=new ofstream[nfactor];for(int l=1;l< nfactor+1;l++){sprintf(fname,"%s/recovered_signals_%d.txt", netdir, l);fev_single[l-1].open(fname, std::ofstream::out);fev_single[l-1]<< line<< endl;}double Vrho[RHO_NBINS], eVrho[RHO_NBINS], Rrho[RHO_NBINS], eRrho[RHO_NBINS], Trho[RHO_NBINS];for(int i=0;i< RHO_NBINS;i++){Vrho[i]=0.;eVrho[i]=0.;Rrho[i]=0.;eRrho[i]=0.;Trho[i]=RHO_MIN+i *RHO_BIN;}double dV, dV1, dV_spin_mtot, nevts, internal_volume;int nT;int countv=0;int cnt=0;int cnt2=0;int cntfreq=0;bool bcut=false;double liveTot=sim.GetMaximum("ifar");double BKG_LIVETIME_yr=liveTot/CYS;double BKG_LIVETIME_Myr=BKG_LIVETIME_yr/(1.e6);cout.precision(14);cout<< "Total live time ---> background
const char * DistributionToString(MDC_DISTRIBUTION n)
TString GetBurst(wavearray< double > &x, TString ifo)
wavearray< double > GetWNB(double frequency, double bandwidth, double duration, int seed=0, bool mode=0)
MDC_COORDINATES GetCoordinatesSystem()
void SetInjJitter(double inj_jitter=MDC_INJ_JITTER)
std::vector< int > idList
void DrawTF(wavearray< double > &x, TString options="")
watplot * DrawFFT(wavearray< double > &x, TString options="ALP", Color_t color=kBlack)
vector< waveform > wfList
static double cosi2e(double cosi)
static void AddCGBurst(wavearray< double > &td, double a, double f, double s, double d=0.)
TString GetBurstLog(source src, double FrameGPS, double SimHpHp, double SimHcHc, double SimHpHc)
static void TimeShift(wavearray< double > &x, double tShift=0.)
double GetSourceListSeed()
MDC_COORDINATES mdc_coordinates
static void AddWGNoise(wavearray< double > &td, double a, double s)
std::vector< std::string > mdcName
MDC_DISTRIBUTION sky_distribution
std::vector< float > rhoList
vector< source > GetSourceList(double start, double stop)
MDC_DISTRIBUTION GetSkyDistribution()
waveform GetWaveform(int ID, int id=0)
void ReadWaveform(wavearray< double > &x, TString fName, double srate)
static double e2cosi(double e)
void SetZoom(double epzoom=EPZOOM)
TString GetParString(TString name, vector< mdcpar > par, bool &error)
std::vector< double > hrssList
std::vector< float > iotaList
void DumpLog(TString fName, TString label="", bool append=false)
void SetInjHrss(double inj_hrss=MDC_INJ_HRSS)
std::vector< source > srcList
static double GetCentralFrequency(wavearray< double > x)
double GetAntennaPattern(TString ifo, double phi, double theta, double psi=0., TString polarization="hp")
std::vector< double > mdcTime
wavearray< double > GetRD(double frequency, double tau, double iota, bool polarization=0)
void SetCoordinatesSystem(MDC_COORDINATES mdc_coordinates=MDC_EARTH)
TString frLabel[NIFO_MAX]
static void AddSGBurst(wavearray< double > &td, double a, double f, double s, double d=0.)
std::vector< float > thList
MDC SetInspiral("inspNameTEST", inspOptions)