Logo coherent WaveBurst  
Library Reference Guide
Logo
CWB_Plugin_supercluster.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #define XIFO 4
20 
21 #pragma GCC system_header
22 
23 #include "cwb.hh"
24 #include "cwb2G.hh"
25 #include "config.hh"
26 #include "network.hh"
27 #include "wavearray.hh"
28 #include "TString.h"
29 #include "TObjArray.h"
30 #include "TObjString.h"
31 #include "TRandom.h"
32 #include "TComplex.h"
33 
34 //!SUPERCLUSTER
35 
36 void
38 
39 // This plugin implements the standard supercluster stage (only 2G)
40 
41  cout << endl;
42  cout << "-----> CWB_Plugin_supercluster.C" << endl;
43  cout << "ifo " << ifo.Data() << endl;
44  cout << "type " << type << endl;
45  cout << endl;
46 
47  if(type==CWB_PLUGIN_CONFIG) {
48  cfg->scPlugin=true; // disable built-in supercluster function
49  }
50 
51  if(type==CWB_PLUGIN_ISUPERCLUSTER) {
52 
53  cout << "type==CWB_PLUGIN_ISUPERCLUSTER" << endl;
54 
55  // import ifile
56  void* gIFILE=NULL; IMPORT(void*,gIFILE)
57  cout << "-----> CWB_Plugin_wavegraph.C -> " << " gIFILE : " << gIFILE << endl;
58  TFile* ifile = (TFile*)gIFILE;
59 
60  // import ifactor
61  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
62  cout << "-----> CWB_Plugin_wavegraph.C -> " << " gIFACTOR : " << gIFACTOR << endl;
63  int ifactor = gIFACTOR;
64 
65  int nIFO = net->ifoListSize(); // number of detectors
66  int rateANA=cfg->inRate>>cfg->levelR;
67  int nRES = net->wdmListSize(); // number of resolution levels
68  int lags = net->getifo(0)->lagShift.size();
69 
70  wavearray<double>* hot[NIFO_MAX]; // temporary time series
71  for(int i=0; i<nIFO; i++) hot[i] = net->getifo(i)->getHoT();
72 
73  int nevt = 0;
74  int nnn = 0;
75  int mmm = 0;
76  size_t count = 0;
77  netcluster wc;
78  netcluster* pwc;
79 
80  for(int j=0; j<(int)lags; j++) {
81 
82  int cycle = cfg->simulation ? ifactor : Long_t(net->wc_List[j].shift);
83 
84  // read cluster metadata
85  if(ifile!=NULL) wc.read(ifile,"coherence","clusters",0,cycle);
86  else wc.read(jfile,"coherence","clusters",0,cycle);
87  // read clusters from temporary job file, loop over TF resolutions
88  if(ifile!=NULL) {
89  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
90  wc.read(ifile,"coherence","clusters",-1,cycle,rateANA>>(i+cfg->l_low));
91  } else {
92  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
93  wc.read(jfile,"coherence","clusters",-1,cycle,rateANA>>(i+cfg->l_low));
94  }
95  if(!cfg->simulation) cout<<"process lag " <<cycle<<" ..."<<endl;
96  cout<<"loaded clusters|pixels: "<<wc.csize()<<"|"<<wc.size()<<endl;
97 
98  // supercluster analysis
99  wc.supercluster('L',net->e2or,cfg->TFgap,false); //likehood2G
100  cout<<"super clusters|pixels: "<<wc.esize(0)<<"|"<<wc.psize(0)<<endl;
101 
102  // release all pixels
103  pwc = net->getwc(j);
104  pwc->cpf(wc, false);
105 
106  net->setDelayIndex(hot[0]->rate());
107  pwc->setcore(false);
108 
109  // apply cuts
110  int psel = 0;
111  while(1) {
112  count = pwc->loadTDampSSE(*net, 'a', cfg->BATCH, cfg->LOUD);
113  psel += net->subNetCut((int)j,cfg->subnet,NULL);
114  int ptot = pwc->psize(1)+pwc->psize(-1);
115  double pfrac = ptot>0 ? double(psel)/double(ptot) : 0.;
116  cout<<"selected pixels: "<<psel<<", fraction: "<<pfrac<< endl;
117  if(count<10000) break;
118  }
119 
120  pwc->defragment(cfg->Tgap,cfg->Fgap); // SK added defragmentation
121 
122  nevt = net->events();
123  nnn += pwc->psize(-1);
124  mmm += pwc->psize(1)+pwc->psize(-1);
125 
126  if(mmm) cout<<"events in the buffer: "<<net->events()<<"|"<<nnn<<"|"<<nnn/double(mmm)<<"\n";
127  else cout<<"events in the buffer: "<<net->events()<<"\n";
128 
129  // store cluster into temporary job file [NEWSS]
130  pwc->write(jfile,"supercluster","clusters",0,cycle);
131  pwc->write(jfile,"supercluster","clusters",-1,cycle);
132  cout<<cycle<<"|"<<pwc->csize()<<"|"<<pwc->size()<<" ";cout.flush();
133 
134  pwc->clear();
135  cout<<endl;
136  }
137  }
138 
139  return;
140 }
141 
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
CWB::config * cfg
virtual size_t defragment(double T, double F, TH2F *=NULL)
T - maximum time gap in seconds F - maximum frequency gap in Hz.
Definition: netcluster.cc:1294
size_t write(const char *file, int app=0)
Definition: netcluster.cc:3008
std::vector< netcluster > wc_List
Definition: network.hh:610
TString("c")
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
netcluster * pwc
Definition: cwb_job_obj.C:38
size_t setcore(bool core, int id=0)
Definition: netcluster.cc:243
int nevt
virtual size_t supercluster(char atype, double S, bool core)
param: statistic: E - excess power, L - likelihood param: selection threshold T for likelihood cluste...
Definition: netcluster.cc:808
int j
Definition: cwb_net.C:28
i drho i
long subNetCut(int lag, float subnet=0.6, float subcut=0.33, TH2F *hist=NULL)
Definition: network.cc:1014
size_t wdmListSize()
Definition: network.hh:446
char ifo[NIFO_MAX][8]
network ** net
NOISE_MDC_SIMULATION.
size_t ifoListSize()
Definition: network.hh:431
#define nIFO
virtual size_t size() const
Definition: wavearray.hh:145
size_t esize(int k=2)
Definition: netcluster.hh:153
jfile
Definition: cwb_job_obj.C:43
int simulation
Definition: config.hh:199
double Fgap
Definition: config.hh:136
wavearray< double > hot[2]
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:69
i() int(T_cor *100))
wavearray< double > * getHoT()
param: no parameters
Definition: detector.hh:175
void setDelayIndex(double rate)
param: MDC log file
Definition: network.cc:2896
int ifactor
const int NIFO_MAX
Definition: wat.hh:22
int BATCH
Definition: config.hh:245
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:117
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
Definition: netcluster.cc:3317
wavearray< double > lagShift
Definition: detector.hh:369
double e2or
Definition: network.hh:584
TFile * ifile
size_t rateANA
Definition: test_config1.C:21
size_t size()
Definition: netcluster.hh:147
double subnet
Definition: config.hh:247
size_t events()
Definition: network.hh:329
int l_low
Definition: config.hh:155
int LOUD
Definition: config.hh:246
size_t csize()
Definition: netcluster.hh:151
int gIFACTOR
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
SUPERCLUSTER.
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
double Tgap
Definition: config.hh:134
netcluster wc
const int nRES
Definition: revMonster.cc:7
size_t read(const char *)
Definition: netcluster.cc:3115
double TFgap
Definition: config.hh:138
size_t psize(int k=2)
Definition: netcluster.hh:163
size_t inRate
Definition: config.hh:132
int levelR
Definition: config.hh:152
bool scPlugin
Definition: config.hh:368
void clear()
Definition: netcluster.hh:427