JCalibrateNB.cc

Go to the documentation of this file.

#include <iostream>

#include "JDAQ/JDAQEventIO.hh"
#include "Jeep/JParser.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JROOT/JManager.hh"
#include "JPhysics/JDispersion.hh"
#include "JTrigger/JPMTSelector.hh"
#include "JTools/JCombinatorics.hh"
#include "JMath/JMathToolkit.hh"

#include "JNanobeacon.hh"

#include "TH1.h"
#include "TH2.h"
#include "TFile.h"

using namespace std;
using namespace JPP;
using namespace JDETECTOR;
using namespace KM3NETDAQ;
using namespace JNANOBEACON;

int main(int argc , char** argv){

  JMultipleFileScanner<JDAQEvent>  inputFiles;  
  JLimit_t&                        numberOfEvents = inputFiles.getLimit();
  string                           detectorFile;
  string                           outFile;
  JDetector                        detector;
  double                           depth;
  double                           wavelength;
  JPMTSelector                     selector;
  JPMTSelector                     default_pmts;

  default_pmts.push_back(JPMTIdentifier_t(-1,0));
  default_pmts.push_back(JPMTIdentifier_t(-1,3));
  default_pmts.push_back(JPMTIdentifier_t(-1,4));
    
  try {

    JParser<> zap;
    
    zap['f'] = make_field(inputFiles                           );
    zap['a'] = make_field(detectorFile                         );
    zap['n'] = make_field(numberOfEvents                       ) = JLimit::max();
    zap['o'] = make_field(outFile                              ) = "out.root";
    zap['h'] = make_field(depth         , "detector depth [m]" ) = 2500.0;
    zap['l'] = make_field(wavelength    , "LED wavelength [nm]") = 400.0;
    zap['T'] = make_field(selector                             ) = default_pmts; // for example '-1 0 -1 3 -1 4', to select PMTs 0, 3 and 4 on every DOM

    zap(argc,argv);
  }
  catch(const exception &error) {
    ERROR(error.what() << endl);
  }

  typedef JHitL0 JHit_t;
  typedef vector<JHit_t> JHitBuffer;

  JBuildL0<JHit_t> builder;
  
  JHitBuffer triggerBuffer;
  JHitBuffer snapshotBuffer;
  
  load(detectorFile , detector);

  JModuleRouter moduleRouter(detector);

  JDispersion jd (0.1 * depth);

  double cw = getSpeedOfLight() / jd.getIndexOfRefractionGroup(wavelength);

  int num_of_floors = getNumberOfFloors(detector);
  JCombinatorics c(num_of_floors);
  int npairs = c.getNumberOfPairs();
  c.sort(JNANOBEACON::comparepair);

  JManager < int, TH2D >* timeDifferences;
  timeDifferences = new JManager < int, TH2D > ( new TH2D("%", "", npairs, 0.5, npairs+0.5, 501, -250.5, 250.5) );
  
  while(inputFiles.hasNext()){
    
    JDAQEvent* event = inputFiles.next();
  
    triggerBuffer .clear();
    snapshotBuffer.clear();
    
    builder(*event, moduleRouter, true  , back_inserter(snapshotBuffer));
    builder(*event, moduleRouter, false , back_inserter( triggerBuffer)); 

    JHitBuffer::iterator __end = partition(triggerBuffer.begin(), triggerBuffer.end(), selector);
    if(triggerBuffer.begin() != __end){

      sort(triggerBuffer.begin(), __end, less<JHit>());
      JHitBuffer::iterator triggeredHit = triggerBuffer.begin();
                
      if( moduleRouter.hasModule( triggeredHit -> getModuleID() ) ){
        const JModule& triggeredModule = moduleRouter.getModule( triggeredHit -> getModuleID() );

        for(JHitBuffer::const_iterator snapshotHit = snapshotBuffer.begin(); snapshotHit != snapshotBuffer.end(); ++snapshotHit){
                         
          if( moduleRouter.hasModule( snapshotHit -> getModuleID() ) ){
            const JModule& snapshotModule = moduleRouter.getModule( snapshotHit -> getModuleID() );

            if( triggeredModule.getString() == snapshotModule.getString() ){

              //double ToF = triggeredHit->getDistance(*snapshotHit) / cw;
              double Dt  = snapshotHit->getT() - triggeredHit->getT();
              double ToF = getDistance( triggeredHit->getPosition(), snapshotHit->getPosition() ) / cw;
              //indexes are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153
              int xbin = c.getIndex( triggeredModule.getFloor()-1 , snapshotModule.getFloor()-1 ) + 1;
              if( triggeredModule.getFloor() < snapshotModule.getFloor() ){
                (*timeDifferences)[triggeredModule.getString()]->Fill(xbin,Dt-ToF);
              }
            }
          } else {
            FATAL("JModuleRouter trying to access non existing identifier: " << snapshotHit->getModuleID() );
          }
        }

      } else {
        FATAL("JModuleRouter trying to access non existing identifier: " << triggeredHit->getModuleID() );
      }     
    }
  }
  TFile output(outFile.c_str() , "recreate") ;
  timeDifferences -> Write(output);
}