Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JClusterBuilder.cc
Go to the documentation of this file.
1 /**
2  * \file
3  Example program to demonstrate the usage of JClusterBuilder
4 
5  If the option -v ("verbose") is used, the number of coincidence clusters
6  is printed for each DOM for each time slice.
7 **/
8 
9 // c++ standard library
10 #include <string>
11 
12 // JPP
13 #include "Jeep/JParser.hh"
15 #include "JDAQ/JDAQTimesliceIO.hh"
18 #include "JSupport/JSupport.hh"
19 #include "JDetector/JDetector.hh"
22 #include "JMonitor/JCluster.hh"
24 
25 // root
26 #include "TH2D.h"
27 
28 // namespaces
29 using namespace std ;
30 using namespace KM3NETDAQ ;
31 using namespace JPP;
32 using namespace JDETECTOR ;
33 using namespace JTRIGGER ;
34 using namespace JMONITOR ;
35 
36 int main(int argc, char **argv) {
37  const unsigned int max_multiplicity = 31 ;
38 
39  //------------------------------------------------------------
40  // READ COMMAND-LINE OPTIONS
41  //------------------------------------------------------------
43  string detectorFile ;
44  double window ;
45  int maxnslices ;
46  bool verbose ;
47  string ofname ;
48 
49  try {
50  JParser<> zap;
51 
52  zap['f'] = make_field(input) ;
53  zap['m'] = make_field(maxnslices) = 100 ;
54  zap['a'] = make_field(detectorFile) ;
55  zap['w'] = make_field(window) ;
56  zap['v'] = make_field(verbose) ;
57  zap['o'] = make_field(ofname) = "out.root" ;
58  zap(argc, argv);
59  }
60  catch(const exception &error) {
61  cerr << error.what() << endl ;
62  cout << "Use option -h! to display command line options." << endl ;
63  exit(1) ;
64  }
65 
66  //------------------------------------------------------------
67  // OPEN DETECTOR FILE
68  //------------------------------------------------------------
70  try {
71  load(detectorFile, detector);
72  }
73  catch(const JException& error) {
74  cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl ;
75  exit(1) ;
76  }
77  JModuleRouter moduleRouter(detector) ;
78 
79  //------------------------------------------------------------
80  // ALLOCATE HISTOGRAMS
81  //------------------------------------------------------------
82  const int ncolors = 5 ;
83  const int nice_colors[ncolors] = { kRed, kBlue, kBlack, kViolet, kCyan } ;
84 
85  // demonstration plots: the ToT distribution for hits in clusters
86  // of different sizes
87  vector<TH1D*> hToT(max_multiplicity+1,NULL) ;
88  for( unsigned int m=0; m<max_multiplicity+1; ++m ) {
89  char hname[200] ;
90  sprintf( hname, "hToT_m%u", m ) ;
91  char htitle[300] ;
92  sprintf( htitle, "Exclusive %u-fold coincidence clusters;ToT [ns];a.u. (normalized)", m ) ;
93  hToT[m] = new TH1D( hname, htitle, 256, -0.5, 255.5 ) ;
94  hToT[m]->SetLineColor( nice_colors[m%ncolors] ) ;
95  hToT[m]->SetLineWidth(2) ;
96  }
97 
98  // demonstration plots: the time distribution of the hits in the clusters
99  vector<TH1D*> ht(max_multiplicity+1,NULL) ;
100  for( unsigned int m=0; m<max_multiplicity+1; ++m ) {
101  char hname[200] ;
102  sprintf( hname, "ht_m%u", m ) ;
103  char htitle[300] ;
104  sprintf( htitle, "Exclusive %u-fold coincidence clusters;Time after first hits [ns];a.u. (normalized)", m ) ;
105  const int margin = 5 ; // ns
106  double xmin = -margin ;
107  double xmax = ceil(window) + margin ;
108  int nbins = (int) round(xmax-xmin) ;
109  ht[m] = new TH1D( hname, htitle, nbins, xmin, xmax ) ;
110  ht[m]->SetLineColor( nice_colors[m%ncolors] ) ;
111  ht[m]->SetLineWidth(2) ;
112  }
113 
114 
115  // another demonstration plot showing the relation between cluster size
116  // and cluster multiplicity (they are practically always the same)
117  TH2D hSizeVsMultiplicity("hSizeVsMultiplicity",";cluster size;cluster multiplicity",
118  32,-0.5,31.5,
119  32,-0.5,31.5 ) ;
120  hSizeVsMultiplicity.SetOption("colz") ;
121 
122  //------------------------------------------------------------
123  // READ FILES
124  //------------------------------------------------------------
125  cout << endl ;
126  cout << "---------- Reading file(s) ----------" << endl ;
127 
129 
130  JMONITOR::JClusterBuilder cluster_builder(window,true) ;
131 
132  unsigned int nTS = 0 ; // time slice counter
133  while( scan.hasNext() ) {
134  JDAQTimeslice* ts = scan.next() ;
135  if( (int)nTS == maxnslices ) break ;
136  ++nTS ;
137 
138  if( verbose ) {
139  cout << "------ Frame index = " << ts->getFrameIndex() << endl ;
140  }
141 
142  // loop over the JDAQSuperFrames
143  for(JDAQTimeslice::const_iterator sf = ts->begin() ; sf!=ts->end() ; ++sf ) {
144  // ignore frames without hits
145  if( sf->size() == 0 ) continue ;
146 
147  int moduleID = sf->getModuleID() ;
148  int localID = moduleRouter.getAddress(moduleID).first ;
149  JModule module = detector[localID] ;
150 
151  // build clusters
152  cluster_builder.reset(*sf,module) ;
153 
154  // print number of clusters found in this frame
155  if( verbose ) {
156  cout << "--- " << "S" << module.getString() << "F" << module.getFloor() << endl ;
157  cout << setw(10) << "multiplicity"
158  << setw(20) << "excl. nclusters"
159  << setw(20) << "incl. nclusters"
160  << endl ;
161  for( unsigned int m=2; m<=max_multiplicity; ++m) {
162  if( cluster_builder.getNclusters(m) != 0 ) {
163  cout << setw(10) << m
164  << setw(20) << cluster_builder.getNclusters(m)
165  << setw(20) << cluster_builder.getInclusiveNclusters(m)
166  << endl ;
167  }
168  }
169  }
170 
171  /**
172  loop over all clusters
173  the begin/end_inclusive_m(M) methods return the begin and end iterator
174  for clusters of multiplicity M or higher
175  **/
176  for( vector<JCluster>::const_iterator it=cluster_builder.begin_inclusive_m(0); it!=cluster_builder.end_inclusive_m(0); ++it ) {
177  hSizeVsMultiplicity.Fill( it->size(), it->getMultiplicity() ) ;
178  }
179 
180  /**
181  loop over clusters of a given multiplicity
182  the begin/end_m(M) methods return the begin and end iterator for
183  clusters of multiplicity exactly M
184  **/
185  for( unsigned int m=0; m<=max_multiplicity; ++m) {
186  for( vector<JCluster>::const_iterator it=cluster_builder.begin_m(m); it!=cluster_builder.end_m(m); ++it ) {
187  // loop over the hits in the cluster to fill ToT histogram
188  for( JHitL1::const_iterator hit=it->begin(); hit!=it->end(); ++hit ) {
189  hToT[m]->Fill( hit->getToT() ) ;
190  }
191  // loop over the hits after the first hit to fill dt histogram
192  if( it->size() > 1 ) {
193  JHitL1::const_iterator first_hit = it->begin() ;
194  JHitL1::const_iterator hit(first_hit) ;
195  for( ++hit; hit!=it->end(); ++hit ) {
196  ht[m]->Fill( hit->getT() - first_hit->getT() ) ;
197  }
198  }
199  }
200  }
201 
202  } // end of loop over JDAQSuperFrames
203 
204  if( verbose ) cout << endl ;
205 
206  } // end of loop over JDAQTimeSlices
207 
208  cout << "Read " << nTS << " time slices." << endl ;
209 
210  TFile* f = new TFile( ofname.c_str(), "recreate" ) ;
211 
212  hSizeVsMultiplicity.Write() ;
213 
214  for( unsigned int m=0; m<=max_multiplicity; ++m) {
215  if( hToT[m]->GetEntries()>0 ) {
216  // normalize
217  hToT[m]->Scale( 1.0/hToT[m]->Integral() ) ;
218  // write
219  hToT[m]->Write() ;
220  }
221  if( ht[m]->GetEntries()>0 ) {
222  // normalize
223  ht[m]->Scale( 1.0/ht[m]->Integral() ) ;
224  // write
225  ht[m]->Write() ;
226  }
227  }
228 
229  f->Close() ;
230  delete f ;
231  cout << "Output in '" << ofname << "'." << endl ;
232 
233  cout << endl ;
234 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
ROOT TTree parameter settings of various packages.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:57
unsigned int getInclusiveNclusters(const unsigned int multiplicity) const
return the number of clusters with at least the given multiplicity
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
void reset(const JDAQSuperFrame &frame, const JModule &module)
This is a way to re-use the allocated memory.
exit
Definition: JPizza.sh:36
do set_array DAQHEADER JPrintDAQHeader f
Definition: JTuneHV.sh:79
Data structure for detector geometry and calibration.
vector< JCluster >::const_iterator end_m(unsigned int multiplicity) const
returns end iterator for clusters with exactly the given multiplicity
int getFrameIndex() const
Get frame index.
int first
index of module in detector data structure
virtual bool hasNext() override
Check availability of next element.
int verbose
Definition: elog.cc:70
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Data time slice.
vector< JCluster >::const_iterator end_inclusive_m(unsigned int multiplicity) const
returns end iterator for clusters with at least the given multiplicity
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
int getString() const
Get string number.
Definition: JLocation.hh:134
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Utility class to parse command line options.
unsigned int getNclusters(const unsigned int multiplicity) const
return the number of clusters with exactly the given multiplicity
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
virtual const pointer_type & next() override
Get next element.
vector< JCluster >::const_iterator begin_inclusive_m(unsigned int multiplicity) const
returns begin iterator for clusters with at least the given multiplicity
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
Local coincidence cluster builder.
vector< JCluster >::const_iterator begin_m(unsigned int multiplicity) const
returns begin iterator for clusters with exactly the given multiplicity
int main(int argc, char *argv[])
Definition: Main.cpp:15