MonopoleAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    MonopoleAna.cxx
// \brief   Analyse monopole monte carlo events
// \version $Id: MonopoleAna.cxx,v 1.6 2012-08-23 01:52:17 bckhouse Exp $
// \author  $Author: bckhouse $
////////////////////////////////////////////////////////////////////////
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "Utilities/func/MathUtil.h"
#include "nusimdata/SimulationBase/MCTrajectory.h"
#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/RecoHit.h"

#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Core/ModuleMacros.h"

#include "TLorentzVector.h"
#include "TTree.h"
#include "TH1I.h"
#include "TH1F.h"

#include <iostream>
#include <string>
class TTree;

namespace mcchk {
  class MonopoleAna : public art::EDAnalyzer
  {
  public:
    explicit MonopoleAna(fhicl::ParameterSet const &pset);
    virtual ~MonopoleAna();
                        
    void analyze(const art::Event &evt);
    void beginJob();
    void reconfigure(fhicl::ParameterSet const &pset);
      
  private:

    TTree*    MM_MC;
    TTree*    MMH_cheater;

    std::string fClusterInput;
    std::string fFLSHitListLabel;
    std::string fRawDataLabel;

    double x0,y0,z0;
    double cosx,cosy,cosz,momentum;
    double pathLength;       //pathlength inside FD
    double mcTi, mcTf;       //MC entry and exit time of monopole to FD
    double rcTi, rcTf;       //RC entry and exit time of monopole to FD
    int nModule,nCell,nSCell;
    double dE, dPE;

    int plane_id, cell_id;
    double entry_t;
    double exit_t;
    double entry_x,exit_x,entry_y,exit_y,entry_z,exit_z;
    double dEnergy, chpl;
  };
}


namespace mcchk{

  MonopoleAna::MonopoleAna(fhicl::ParameterSet const& pset) 
    : EDAnalyzer(pset)
  {
    reconfigure(pset);
  }


  //......................................................................
  void MonopoleAna::reconfigure(fhicl::ParameterSet const &pset)
  {
    fClusterInput    = pset.get< std::string >("ClusterInput") ;
    fFLSHitListLabel = pset.get< std::string >("FLSHitListLabel") ;
    fRawDataLabel    = pset.get< std::string >("RawDataLabel");
  }

  //......................................................................

  MonopoleAna::~MonopoleAna(){}
  

  //......................................................................

  void MonopoleAna::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    std::cout<<"begin job!"<<std::endl;
    //The following two ntuples can give info: 
    // efficiency VS track length 
 
    MM_MC = tfs->make<TTree>("MM_MC","monopole_mc_info");
    MM_MC -> Branch("x0",&x0,"x0/D");
    MM_MC -> Branch("y0",&y0,"y0/D");
    MM_MC -> Branch("z0",&z0,"z0/D");
    MM_MC -> Branch("pathLength",&pathLength,"pathLength/D");
    MM_MC -> Branch("cosx",&cosx,"cosx/D");
    MM_MC -> Branch("cosy",&cosy,"cosy/D");
    MM_MC -> Branch("cosz",&cosz,"cosz/D");   
    MM_MC -> Branch("momentum",&momentum,"momentum/D");
    MM_MC -> Branch("mcTi",&mcTi,"mcTi/D");
    MM_MC -> Branch("mcTf",&mcTf,"mcTf/D");
    MM_MC -> Branch("rcTi",&rcTi,"rcTi/D");
    MM_MC -> Branch("dE",&dE,"dE/D");
    MM_MC -> Branch("dPE",&dPE,"dPE/D");
    MM_MC -> Branch("nModule",&nModule,"nModule/I"); 
    MM_MC -> Branch("nCell",&nCell,"nCell/I"); 
    MM_MC -> Branch("nSCell",&nSCell,"nSCell/I"); 


    MMH_cheater = tfs->make<TTree>("MMH_cheater","monopole_hits_mc_info");
    MMH_cheater -> Branch("plane",&plane_id,"plane/I");
    MMH_cheater -> Branch("cell",&cell_id,"cell/I");
    MMH_cheater -> Branch("x0",&entry_x,"x0/D");
    MMH_cheater -> Branch("y0",&entry_y,"y0/D");
    MMH_cheater -> Branch("z0",&entry_z,"z0/D");
    MMH_cheater -> Branch("t0",&entry_t,"t0/D");
    MMH_cheater -> Branch("x1",&exit_x,"x1/D");
    MMH_cheater -> Branch("y1",&exit_y,"y1/D");
    MMH_cheater -> Branch("z1",&exit_z,"z1/D");
    MMH_cheater -> Branch("t1",&exit_t,"t1/D");
    MMH_cheater -> Branch("energydump",&dEnergy,"dEnergy/D");
    MMH_cheater -> Branch("pathlength",&chpl,"pathlength/D");

  }

  void MonopoleAna::analyze(const art::Event& evt)
  {
    //Initialization of MM_MC (each entry corresponds to a single MM):
    mcTi = 9999999;
    mcTf = -9999999;
    pathLength = 0;
    nCell = 0;
    nSCell = 0;
    nModule = 0;
    dE =0;
    dPE =0;
    x0 = 0;
    y0 = 0;
    z0 =0;
    momentum = 0;
    cosx = 0;
    cosy = 0;
    cosz = 0;


    //make a backtracker
    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

    for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
      const sim::Particle *p = (*i).second;
      if (p->PdgCode()!=42 ) continue;

      momentum = p->P();
      x0 = p->Vx();
      y0 = p->Vy();
      z0 = p->Vz();
      cosx = (p->Px())/momentum;
      cosy = (p->Py())/momentum;
      cosz = (p->Pz())/momentum; 

      const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
      for (unsigned h =0; h!= flshits.size(); ++h) {
        plane_id = flshits[h].GetPlaneID();
        cell_id  = flshits[h].GetCellID();
        entry_x  = flshits[h].GetEntryX();
        exit_x   = flshits[h].GetExitX();
        entry_y  = flshits[h].GetEntryY();
        exit_y   = flshits[h].GetExitY();
        entry_z  = flshits[h].GetEntryZ();
        exit_z   = flshits[h].GetExitZ();
        entry_t  = flshits[h].GetEntryT();
        exit_t   = flshits[h].GetExitT();
        dEnergy  = flshits[h].GetEdep();
        chpl     = flshits[h].GetTotalPathLength();
        MMH_cheater->Fill();

        pathLength += chpl;
        if (mcTi>entry_t) mcTi = entry_t;
        if (mcTf<exit_t) mcTf = exit_t;
      }
    }

    // define vector holding the clusters
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fClusterInput, slices);
    art::ServiceHandle<calib::Calibrator> calib;

    const int sliceMax = slices->size();
    for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){

      const rb::Cluster& slice = (*slices)[sliceIdx];

      if (slice.IsNoise()) continue;

      nCell = slice.NCell();
      nModule =0;
      unsigned short plane,cell;
      int pStore = -99;
      int mStore = -99;

      for(int hitIdx = 0; hitIdx < nCell; ++hitIdx){
        const art::Ptr<rb::CellHit> chit = slice.Cell(hitIdx);
        plane = chit->Plane();
        cell  = chit->Cell();
        int module = cell/32;
        if (plane!=pStore or module!=mStore) nModule++;

        try {
          if(!bt->IsNoise(chit) || (chit->ADC()>3100 && chit->TNS() > mcTi-1000 && chit->TNS() < mcTf+1000) ) {

            std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(chit);
            TVector3                     position = bt->HitToXYZ(chit);

            int mcPE=0;

            for (unsigned int i=0; i!=(bt->HitToPhotonSignal(chit)).size();++i)
              mcPE += (bt->HitToPhotonSignal(chit))[i].NPhoton();
            
            double w = position[0];
            if (chit->View() == geo::kX) w= position[1];
            
            rb::RecoHit recohit(calib->MakeRecoHit(*chit,w));

            const double rbTime  = recohit.T();

            if (!hitIdx) rcTi = rbTime;
            rcTf = rbTime;

            dE  += (trackides[0].energy) / (trackides[0].energyFrac);
            dPE += mcPE;

            if (chit->ADC(3)==4095) nSCell++;

          }

          pStore = plane;
          mStore = module;


        }

        catch (...) {continue;}
      
      }
    }

    MM_MC->Fill();

    return;
  }





  DEFINE_ART_MODULE(MonopoleAna)

}