JTurbotIter.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <map>
#include <vector>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TMath.h"

#include "JDAQ/JDAQ.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDAQHitRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL1.hh"
#include "JLang/JSinglePointer.hh"
#include "JGizmo/JManager.hh"
#include "JTools/JWeight.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JTreeScannerInterface.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to search for correlations between triggered events and timeslice data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JSingleFileScanner<> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  double             TMaxLocal_ns;
  string             pmtFile;
  int                numberOfTimeslices;
  double             binWidth_ns;
  multimap<int, int> mask;
  double             Pmin;
  JROOTClassSelector selector;
  int                debug;
  string             peaks;
  int                multiPeaks;

  try { 

    JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(TMaxLocal_ns)        = 20.0;
    zap['P'] = make_field(pmtFile)             = "";
    zap['N'] = make_field(numberOfTimeslices)  = 100;
    zap['W'] = make_field(binWidth_ns)         = 10e3;
    zap['M'] = make_field(mask)                = JPARSER::initialised();
    zap['p'] = make_field(Pmin)                = 1.0e-7;
    zap['C'] = make_field(selector)            = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['d'] = make_field(debug)               =  2;
    zap['O'] = make_field(peaks)               = ""; //Make a small summary file, with the list of the peaks
    zap['m'] = make_field(multiPeaks)          =  1; //Enables final multi-peaks analysis

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  cout.tie(&cerr);

  
  map<int, int> MASK;                // mask PMTs with permament high-rate veto.

  const int WR     =  0x80000000;    // White Rabbit

  //MASK[808969848]  =  0x00000020;    // floor 03, pmt 05
  //MASK[809544061]  =  0x00000080;    // floor 05, pmt 07
  //MASK[808432835]  =  0x00004000;    // floor 18, pmt 14

  for (multimap<int, int>::const_iterator i = mask.begin(); i != mask.end(); ++i) {
    JBit(i->second).set(MASK[i->first]);
  }
  
  for (map<int, int>::const_iterator i = MASK.begin(); i != MASK.end(); ++i) {
    DEBUG("Mask " << setw(8) << i->first << " 0x" << hex << i->second << dec << endl);
  }
  
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  
  const JDAQHitRouter router(detector);

  const double TMax_ns = max(binWidth_ns, getMaximalTime(detector));

  typedef  double    hit_type;
  
  JBuildL0<hit_type> buildL0;
  JBuildL1<hit_type> buildL1(TMaxLocal_ns, true);
  vector  <hit_type> data;  
  
  typedef JManager<int, TH1D>  JManager_t;

  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
  const Int_t    nx   = (Int_t) ((xmax - xmin) / binWidth_ns);
  
  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));

  JTreeScanner<>::debug = debug;
  
  JSinglePointer< JTreeScannerInterface<JDAQTimeslice> > ps;

  JAutoTreeScanner<JDAQTimeslice> zmap = JType<JDAQTimesliceTypes_t>();
  
  if (selector == "") {

    if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice,   JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
    } else {
      FATAL("No timeslice data." << endl);
    }

    NOTICE("Selected class " << ps->getClassName() << endl);

  } else {

    ps = zmap[selector];
    
    ps->configure(inputFile);
  }

  //Vector with OOS DOMs IDs
  vector<int> oosDomsId;

  /*   Updated method  :
   *   OOS DOMs can have an impact on the histogram of other DOMs   
   *   especialy OOS DOMs on adjacent In Sync or OOS DOMs (presence of ghost peaks)
   *   To remove these ghost peaks, the idea is to iterate the method, removing 
   *   at each loop the most OOS DOM, and recomputing the histograms without 
   *   them (removing events where there are involved)
   *   Once there is no more OOS DOMs in the remaining DOMs, end of precedure.   
   *   Steps of procedure (from J. Hofestaedt, code updated by R. Le Breton)    
   *   1) Check if any DOM shows a peak away from 0 (OOS)   
   *      If no OOS -> finish procedure                                         
   *   2) If at least one OOS, then identify the DOM that is most OOS,          
   *      ie. highest peak (other than delta time = 0)                            
   *   3) Remove the most OOS DOM and calculate the histograms again but do not
   *      consider the OOS DOM (in later iterations more than one DOM is OOS and not   
   *      considered for histogram filling).
   *   4) start with 1) again
   */

  typedef map<int, vector<double> > map_type;    // frame index -> event times

  bool OOS = true; //Boolean: to continue or not the procedure
  JPMTParametersMap parameters;
  map_type buffer;
    
  vector<vector<double> > peaksPerDoms;

  if (pmtFile != "") {
    try {
      parameters.load(pmtFile.c_str());
    }
    catch(const JException& error) {}
  }
  
  do{
    for (JTreeScanner<JDAQEvent> in(inputFile.getFilename(), inputFile.getLimit()); in.hasNext(); ) {
      
      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
      
      JDAQEvent* event = in.next();
      
      // Average time of triggered hits
      
      double t0 = 0.0;
      bool test_OOS = false; //If an OOS DOM is involved in an event, remove this event
      
      for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
        t0 += getTime(*hit, router.getPMT(*hit));

        if(std::find(oosDomsId.begin(), oosDomsId.end(), hit->getModuleIdentifier()) != oosDomsId.end()) {
          test_OOS = true;
          continue;
        }
      }

      if(test_OOS) continue; //If true, event not taken into account

      t0 /= event->size<JDAQTriggeredHit>();
      t0 += event->getFrameIndex() * getFrameTime();
      
      buffer[event->getFrameIndex()].push_back(t0);
    }    
    STATUS(endl);
    
    //Events : start from the beginning
    ps->rewind();
    
    while (ps->hasNext()) {
      
      STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
      
      JDAQTimeslice* timeslice = ps->next();
      
      map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
      map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
      
      int number_of_events = 0;
      
      for (map_type::const_iterator i = p; i != q; ++i) {
        number_of_events += i->second.size();
      }
      
      if (number_of_events == 0) {
        continue;
      }
      
      for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
        
        // Test that none of the PMTs has high-rate veto.
        
        if ((frame->getStatus() & ~MASK[frame->getModuleID()] & ~WR) == 0) {
          
          data.clear();
          
          buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
          
          TH1D* h1 = manager[frame->getModuleID()];
          
          for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
            
            const double t1 = *hit + frame->getFrameIndex() * getFrameTime();
            
            for (map_type::const_iterator i = p; i != q; ++i) {
              for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
                
                const double t0 = *j;
                
                h1->Fill(t1 - t0);
              }
            }
          }
        }
      }
    }

    buffer.clear(); //Clear the buffer for the next loop, if any

    STATUS(endl);

    if (pmtFile != "") {
      
      Double_t most_OOS_peak = 0;
      Int_t most_OOS_ID = 0;
      
      TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
      
      string option = "L";
      
      if (debug < debug_t) {
        option += "Q";
      }
      
      for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
        
        bool status = false;
        
        JManager_t::iterator p = manager.find(module->getID());
        
        if (p != manager.end()) {
          
          TH1D* h1       = p->second;
          
          if (h1->GetEntries() != 0) {
            
            // start values
            
            Double_t ymax  =   0.0;
            Double_t x0    =   0.0;    // [ns]
            Double_t sigma = 250.0;    // [ns]
            Double_t ymin  =   0.0;
            
            for (int i = 1; i != h1->GetNbinsX(); ++i) {
              
              const Double_t x = h1->GetBinCenter (i);
              const Double_t y = h1->GetBinContent(i);
              
              if (y > ymax) {
                ymax = y;
                x0   = x;
              }
              
              ymin += y;
            }
            
            ymin /= h1->GetNbinsX();
            
            f1.SetParameter(0, ymax);
            f1.SetParameter(1, x0);
            if (binWidth_ns < sigma) 
              f1.SetParameter(2, sigma);
            else
              f1.FixParameter(2, binWidth_ns/sqrt(12.0));
            f1.SetParameter(3, ymin);
            
            // fit
            
            h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
            
            status = (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() &&
                      f1.GetParameter(0)       >= f1.GetParameter(3));
          }

          if (h1->GetEntries() != 0) {
            
            // look for peaks at offsets equal to a multiple of frame time
            JWeight           bg;    // background 
            map<int, JWeight> ps;    // signal(s)

            for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {

              const Double_t x = h1->GetBinCenter (i);
              const Double_t y = h1->GetBinContent(i);

              while (x > (ns + 1) * getFrameTime() - TMax_ns) { ++ns; }

              if (fabs(x - ns * getFrameTime()) < TMax_ns)
                ps[ns].put(y);
              else
                bg    .put(y);
            }

            for (map<int, JWeight>::iterator i = ps.begin(); i != ps.end(); ) {

              const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
              const Double_t P = TMath::PoissonI(i->second.getTotal(), noise);

              if (P < Pmin && i->second.getTotal() > noise) //Avoid edge effect (second bolean)
                ++i;               // keep 
              else
                ps.erase(i++);     // remove
            }

            if (ps.size() != 1) {
              ERROR("Module " << setw(8) << module->getID() << " number of peaks " << ps.size() << " != 1" << endl);
            }
            
            for (map<int, JWeight>::const_iterator i = ps.begin(); i != ps.end(); ++i) {
              
              vector<double> subPeaksPerDoms;
              const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
              
              if (ps.size() != 1) {
                WARNING("Peak at " << setw(4) << i->first << " [frame time] S/N = "
                        << FIXED(7,1) << i->second.getTotal()                                 << " / "
                        << FIXED(7,1) << noise << endl);
              }
              
              subPeaksPerDoms.push_back(module->getID());
              subPeaksPerDoms.push_back(ps.size());
              subPeaksPerDoms.push_back(i->first);
              subPeaksPerDoms.push_back(i->second.getTotal());
              subPeaksPerDoms.push_back(noise);
              peaksPerDoms.push_back(subPeaksPerDoms);
            }
          }
          
          //if DOM is OOS 
          if (!status) {
            //if OOS peak greater than the last and not in OOS DOMs vector
            if ((f1.GetParameter(0) > most_OOS_peak) && std::find(oosDomsId.begin(), oosDomsId.end(), module->getID()) == oosDomsId.end()){
              
              most_OOS_peak = f1.GetParameter(0);
              most_OOS_ID = module->getID();
            }
          }
        }
      }
      
      if (most_OOS_ID != 0) {
        oosDomsId.push_back(most_OOS_ID);//exclude the DOM to compute the histograms
        
        NOTICE("Module " << most_OOS_ID << " set QEs of all PMTs to zero." << endl);
        
        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          parameters[JPMTIdentifier(most_OOS_ID, pmt)].QE = 0.0;
        }
        manager.clear();
        peaksPerDoms.clear();
      }
      else{
        OOS = false;
      }
    }
  }while(OOS);
  
  //If multipeaks analysis is activated (multiPeaks = 1)
  //Make ratio : sum_OOS_peaks / sum_all_peaks
  //If ratio > 0.05, DOM is declared OOS

  if (multiPeaks == 1) {

    NOTICE("Multi-peaks analysis activated" << endl);
    
    vector<int> workingDOMs;
    
    //If DOM is working and not in the list of OOS DOMs
    for(int i = 0; i != int(peaksPerDoms.size()); i++) {
      if((std::find(workingDOMs.begin(), workingDOMs.end(), peaksPerDoms[i][0]) == workingDOMs.end()) &&
         (std::find(oosDomsId.begin(), oosDomsId.end(), peaksPerDoms[i][0]) == oosDomsId.end())) {
        workingDOMs.push_back(peaksPerDoms[i][0]);
      }
    }

    for(int i = 0; i != int(workingDOMs.size()); i++) {

      double ratio = 0.;
      double sumAllPeaks = 0.;

      for(int j = 0; j != int(peaksPerDoms.size()); j++) {
        if (peaksPerDoms[j][0] == workingDOMs[i]){
          if (int(peaksPerDoms[j][2]) != 0) {
            ratio += peaksPerDoms[j][3] - peaksPerDoms[j][4];
          }
          sumAllPeaks += peaksPerDoms[j][3] - peaksPerDoms[j][4];
        }
      }

      ratio /= sumAllPeaks;

      if (ratio > 0.05){
        NOTICE("Module " << workingDOMs[i] << " set QEs of all PMTs to zero." << endl);
        NOTICE("----> Ratio = " << 100*ratio << "%" << endl);
        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          parameters[JPMTIdentifier(workingDOMs[i], pmt)].QE = 0.0;
        }
      }

    }
  }
    
  //To manage non working DOMs 
  NOTICE("Managing non-working DOMs." << endl);
  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
    JManager_t::iterator p = manager.find(module->getID());
    if (p == manager.end()){
      NOTICE("Module " << module->getID() << " set QEs of all PMTs to zero." << endl);
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
      }
    }
  }
  
  ofstream out(pmtFile.c_str());
  
  parameters.comment.add(JMeta(argc, argv));
  
  out << parameters << endl;

  out.close();
  
  if (outputFile != "") {
    manager.Write(outputFile.c_str()); 
  }

  //print small summary file of the peaks
  if (peaks != ""){
    ofstream stream(peaks);

    stream << "#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
    
    for(int i = 0; i != int(peaksPerDoms.size()); i++) {
      for(int j = 0; j != int(peaksPerDoms[i].size()); j++){ 
        if (j == 0 || j == 1 || j == 2 || j == 3) stream << int(peaksPerDoms[i][j]) << " " ;
        else stream << peaksPerDoms[i][j] << " " ; 
      }
        stream << "\n";
      }
    stream.close();
  }

}