JClusterBuilder.cc

Go to the documentation of this file.

/**
 * \file
   Example program to demonstrate the usage of JClusterBuilder

   If the option -v ("verbose") is used, the number of coincidence clusters
   is printed for each DOM for each time slice.
**/

// c++ standard library
#include <string>

// JPP
#include "Jeep/JParser.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JMonitor/JCluster.hh"
#include "JMonitor/JClusterBuilder.hh"

// root
#include "TH2D.h"

// namespaces
using namespace std ; 
using namespace KM3NETDAQ ; 
using namespace JPP;   
using namespace JDETECTOR ;
using namespace JTRIGGER ;
using namespace JMONITOR ;

int main(int argc, char **argv) {
  const unsigned int max_multiplicity = 31 ;

  //------------------------------------------------------------
  //  READ COMMAND-LINE OPTIONS
  //------------------------------------------------------------
  JMultipleFileScanner<> input ;
  string detectorFile ;
  double window ;
  int maxnslices ;
  bool verbose ;
  string ofname ;

  try { 
    JParser<> zap;
    
    zap['f'] = make_field(input) ; 
    zap['m'] = make_field(maxnslices) = 100 ;
    zap['a'] = make_field(detectorFile) ;
    zap['w'] = make_field(window) ;
    zap['v'] = make_field(verbose) ;
    zap['o'] = make_field(ofname) = "out.root" ;
    zap(argc, argv);
  }
  catch(const exception &error) {
    cerr << error.what() << endl ;
    cout << "Use option -h! to display command line options." << endl ;
    exit(1) ;
  }

  //------------------------------------------------------------
  //  OPEN DETECTOR FILE
  //------------------------------------------------------------
  JDetector detector;
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl ;
    exit(1) ;
  }
  JModuleRouter moduleRouter(detector) ;

  //------------------------------------------------------------
  //  ALLOCATE HISTOGRAMS
  //------------------------------------------------------------
  const int ncolors = 5 ;
  const int nice_colors[ncolors] = { kRed, kBlue, kBlack, kViolet, kCyan } ;

  // demonstration plots: the ToT distribution for hits in clusters 
  // of different sizes
  vector<TH1D*> hToT(max_multiplicity+1,NULL) ;
  for( unsigned int m=0; m<max_multiplicity+1; ++m ) {
    char hname[200] ;
    sprintf( hname, "hToT_m%u", m ) ;
    char htitle[300] ;
    sprintf( htitle, "Exclusive %u-fold coincidence clusters;ToT [ns];a.u. (normalized)", m ) ;
    hToT[m] = new TH1D( hname, htitle, 256, -0.5, 255.5 ) ;
    hToT[m]->SetLineColor( nice_colors[m%ncolors] ) ;
    hToT[m]->SetLineWidth(2)   ;
  }

  // demonstration plots: the time distribution of the hits in the clusters
  vector<TH1D*> ht(max_multiplicity+1,NULL) ;
  for( unsigned int m=0; m<max_multiplicity+1; ++m ) {
    char hname[200] ;
    sprintf( hname, "ht_m%u", m ) ;
    char htitle[300] ;
    sprintf( htitle, "Exclusive %u-fold coincidence clusters;Time after first hits [ns];a.u. (normalized)", m ) ;
    const int margin = 5 ; // ns
    double xmin = -margin ;
    double xmax = ceil(window) + margin ;
    int nbins = (int) round(xmax-xmin) ;
    ht[m] = new TH1D( hname, htitle, nbins, xmin, xmax ) ;
    ht[m]->SetLineColor( nice_colors[m%ncolors] ) ;
    ht[m]->SetLineWidth(2)   ;
  }


  // another demonstration plot showing the relation between cluster size 
  // and cluster multiplicity (they are practically always the same)
  TH2D hSizeVsMultiplicity("hSizeVsMultiplicity",";cluster size;cluster multiplicity",
                           32,-0.5,31.5,
                           32,-0.5,31.5 ) ;
  hSizeVsMultiplicity.SetOption("colz") ;

  //------------------------------------------------------------
  //  READ FILES
  //------------------------------------------------------------
  cout << endl ;
  cout << "---------- Reading file(s) ----------" << endl ;

  JMultipleFileScanner<JDAQTimeslice> scan(input) ;

  JMONITOR::JClusterBuilder cluster_builder(window,true) ;

  unsigned int nTS = 0 ; // time slice counter
  while( scan.hasNext() ) {
    JDAQTimeslice* ts = scan.next() ;
    if( (int)nTS == maxnslices ) break ;
    ++nTS ;

    if( verbose ) {
      cout << "------ Frame index = " << ts->getFrameIndex() << endl ;
    }

    // loop over the JDAQSuperFrames
    for(JDAQTimeslice::const_iterator sf = ts->begin() ; sf!=ts->end() ; ++sf ) {
      // ignore frames without hits
      if( sf->size() == 0 ) continue ;

      int moduleID   = sf->getModuleID() ;
      int localID    = moduleRouter.getAddress(moduleID).first ;
      JModule module = detector[localID] ;

      // build clusters
      cluster_builder.reset(*sf,module) ;

      // print number of clusters found in this frame
      if( verbose ) {
        cout << "--- " << "S" << module.getString() << "F" << module.getFloor() << endl ;
        cout << setw(10) << "multiplicity"
             << setw(20) << "excl. nclusters"
             << setw(20) << "incl. nclusters"
             << endl ;
        for( unsigned int m=2; m<=max_multiplicity; ++m) {
          if( cluster_builder.getNclusters(m) != 0 ) {
            cout << setw(10) << m 
                 << setw(20) << cluster_builder.getNclusters(m) 
                 << setw(20) << cluster_builder.getInclusiveNclusters(m)
                 << endl ;
          }
        }
      }

      /**
         loop over all clusters
         the begin/end_inclusive_m(M) methods return the begin and end iterator
         for clusters of multiplicity M or higher
      **/
      for( vector<JCluster>::const_iterator it=cluster_builder.begin_inclusive_m(0); it!=cluster_builder.end_inclusive_m(0); ++it ) {
        hSizeVsMultiplicity.Fill( it->size(), it->getMultiplicity() ) ;
      }

      /**
         loop over clusters of a given multiplicity
         the begin/end_m(M) methods return the begin and end iterator for
         clusters of multiplicity exactly M
      **/
      for( unsigned int m=0; m<=max_multiplicity; ++m) {
        for( vector<JCluster>::const_iterator it=cluster_builder.begin_m(m); it!=cluster_builder.end_m(m); ++it ) {
          // loop over the hits in the cluster to fill ToT histogram
          for( JHitL1::const_iterator hit=it->begin(); hit!=it->end(); ++hit ) {
            hToT[m]->Fill( hit->getToT() ) ;
          }
          // loop over the hits after the first hit to fill dt histogram
          if( it->size() > 1 ) {
            JHitL1::const_iterator first_hit = it->begin() ;
            JHitL1::const_iterator hit(first_hit) ;
            for( ++hit; hit!=it->end(); ++hit ) {
              ht[m]->Fill( hit->getT() - first_hit->getT() ) ;
            }
          }
        }
      }

    } // end of loop over JDAQSuperFrames

    if( verbose ) cout << endl ;

  } // end of loop over JDAQTimeSlices

  cout << "Read " << nTS << " time slices." << endl ;

  TFile* f = new TFile( ofname.c_str(), "recreate" ) ;

  hSizeVsMultiplicity.Write() ;

  for( unsigned int m=0; m<=max_multiplicity; ++m) {
    if( hToT[m]->GetEntries()>0 ) {
      // normalize
      hToT[m]->Scale( 1.0/hToT[m]->Integral() ) ;
      // write
      hToT[m]->Write() ;
    }
    if( ht[m]->GetEntries()>0 ) {
      // normalize
      ht[m]->Scale( 1.0/ht[m]->Integral() ) ;
      // write
      ht[m]->Write() ;
    }
  }

  f->Close() ;
  delete f ;
  cout << "Output in '" << ofname << "'." << endl ;

  cout << endl ;
}