JORCAShowerPositionPDF.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>

#include "evt/Evt.hh"

#include "JDAQ/JDAQEvent.hh"

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

#include "JFit/JPoint4D.hh"

#include "JLang/JTypeList.hh"

#include "JDAQ/JDAQEvent.hh"

#include "JSupport/JSupport.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JTreeScanner.hh"

#include "JTrigger/JTimeConverter.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL0.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JModuleRouter.hh"

#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"

#include <TFile.h>
#include <TH2D.h>

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

  using namespace JPP;

  typedef JTriggeredFileScanner<>                       JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type   multi_pointer_type;

  JTriggeredFileScanner_t   inputFile;
  string                    outputFile;
  string                    detectorFile;
  JRange<double>            t_res_ns;
  JRange<double>            D_m;
  JRange<double>            cosT;
  bool                      PMTweight;  // weight PDF to take into account PMT growing number 
                                        // with respect to vertex
  int                       debug;

  try{
    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)      = "JORCAShowerPositionFit_PDF.root";
    zap['a'] = make_field(detectorFile);
    zap['T'] = make_field(t_res_ns)        = JRange<double>(-100, 100);  // [ns]  
    zap['D'] = make_field(D_m)             = JRange<double>(0, 100);  // [m]  
    zap['C'] = make_field(cosT)            = JRange<double>(-1, 1);   
    zap['d'] = make_field(debug)           = 2;
    zap['P'] = make_field(PMTweight);

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

  cout.tie(&cerr);

  JDetector detector;

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

  const JPMTRouter     pmtRouter(detector);
  const JModuleRouter  moduleRouter(detector);

  typedef vector<JHitL0> JDataL0_t;
  const JBuildL0<JHitL0> buildL0;

  TFile *out = new TFile(outputFile.c_str(), "RECREATE");
  if ( !out->IsOpen() ) printf("Can't open outputfile\n");

  TH2D hPDF2Dist("hPDF2Dist", "PDF for ORCA Shower Position Fit; #delta t [ns]; D [m]", 
                 2000, t_res_ns.getLowerLimit(), t_res_ns.getUpperLimit(),
                 700, D_m.getLowerLimit(), D_m.getUpperLimit());

  while(inputFile.hasNext()){
    
    multi_pointer_type ps = inputFile.next();
    
    JDAQEvent* tev = ps;
    Evt*       event = ps; 
    JTimeConverter converter(*event, *tev);
    
    // Get MC truth neutrino and shower    
    Trk cascade;
    cascade.pos  = event->mc_trks[1].pos;
    cascade.pos -= Vec(0.0, 0.0, -117.164);

    for (unsigned int i = 1; i < event->mc_trks.size(); i++) {
      cascade.E   += event->mc_trks[i].E;
      cascade.dir += event->mc_trks[i].dir * event->mc_trks[i].E;
    }
    cascade.dir.normalize();

    JPosition3D vx = getPosition(cascade.pos);
    double t_vx = converter.putTime(event->mc_trks[1].t);
    JDirection3D dir = getDirection(cascade.dir);

    JDataL0_t dataL0;
    buildL0(JDAQTimeslice(*tev,  true),  moduleRouter, back_inserter(dataL0));
    
    //**************** PDF with real interaction vertex ***************
    const JRotation3D R(dir);
    vx.rotate(R);

    for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
      
      JHitL0 hit(*i);
      hit.transform(R, vx);
      
      const double D = hit.getLength(); 
      const double t_expected = t_vx + D * getInverseSpeedOfLight() * getIndexOfRefraction();
      const double dt = hit.getT() - t_expected; 
      const JDirection3D photonDir(hit.getPosition());
      const double ct = photonDir.getDot(hit.getDirection());
      
      if((D >= D_m.getLowerLimit() && D <= D_m.getUpperLimit()) &&
         (dt >= t_res_ns.getLowerLimit() && dt <= t_res_ns.getUpperLimit()) &&
         (ct >= cosT.getLowerLimit() && ct <= cosT.getUpperLimit())){

        double w = 1;

        if(PMTweight){ 
          w = 1 / (D*D);
        }
        
        hPDF2Dist.Fill(dt, D, w);
      }
    }
  } 
  STATUS(endl);
  
  out->cd();
  hPDF2Dist.Write();
  out->Close();
  
}