FmmTriggerAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       FmmTriggerAna
// Module Type: analyzer
// File:        FmmTriggerAna_module.cc
//
// Generated at Wed Dec 18 15:41:32 2013 by Zukai Wang using artmod
// from cetpkgsupport v1_04_02.
////////////////////////////////////////////////////////////////////////
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"

#include "MCCheater/BackTracker.h"
#include "Calibrator/Calibrator.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "Utilities/func/MathUtil.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"

#include "TVector3.h"
#include "TTree.h"
#include "TProfile.h"

#include <algorithm> 
#include <utility>
#include <cmath>
#include <iostream>

namespace zcl {
  class FmmTriggerAna;
}

class zcl::FmmTriggerAna : public art::EDAnalyzer {
  public:
    explicit FmmTriggerAna(fhicl::ParameterSet const & pset);
    virtual ~FmmTriggerAna();

    void analyze(const art::Event & evt);
    void beginJob();
    void endJob(){};

  private:
    std::string fClusterInput;    ///< Input folder from cluster reco

    //Each entry is corresponding to 1 reconstructed slice
    TTree    * RC_Tree;
    //TProfile * hitTdiff;
    //TProfile * hitTdiff_log;
    //TProfile * DCS_beta;
    //TProfile * DCS_logbeta;

    int    evtID, subRunID, RunID;

    //Trigger decision branches
    int   _sliced;
    bool   _penetrated;
    int    _trueHits;
    int    _nSatHits;
    double _meanADC;
    int    _deltaP;
    double _deltaTNS;
    int    _nxhits_slice, _nyhits_slice;

    double _tmin, _tmax;  //The time window of the slice
    int _pmin, _pmax;

    double beta_MC;
    int    nxhits_MC, nyhits_MC; 

};

zcl::FmmTriggerAna::FmmTriggerAna(fhicl::ParameterSet const & p) :
  EDAnalyzer(p),
  fClusterInput      (p.get< std::string > ("ClusterInput")      )
{
}

zcl::FmmTriggerAna::~FmmTriggerAna()
{
  // Clean up dynamic memory and other resources here.
}

void zcl::FmmTriggerAna::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;

  //Combining cheater info and the RC info to make the data structure better
  //One RC track info per entry, if no tracks reconstructed in this event, fill in 1 entry anyway tagging the 
  //reconstruction failed
  RC_Tree = tfs->make<TTree>("RC_Tree","Information from RC track info");
  RC_Tree -> Branch("event",  &evtID,      "event/I");
  RC_Tree -> Branch("subRun", &subRunID,   "subRun/I");     
  RC_Tree -> Branch("Run",    &RunID,      "Run/I");

  //Trigger threshold parameters
  RC_Tree -> Branch("sliced",      &_sliced,      "sliced/I");
  RC_Tree -> Branch("penetrated",  &_penetrated,  "penetrated/O");
  RC_Tree -> Branch("nSatHits",    &_nSatHits,    "nSatHits/I");
  RC_Tree -> Branch("meanADC",     &_meanADC,     "meanADC/D");
  RC_Tree -> Branch("deltaP",      &_deltaP,      "deltaP/I");
  RC_Tree -> Branch("deltaTNS",    &_deltaTNS,    "deltaTNS/D");
  RC_Tree -> Branch("nxhS",        &_nxhits_slice,    "nxhS/I");
  RC_Tree -> Branch("nyhS",        &_nyhits_slice,    "nyhS/I");
  RC_Tree -> Branch("trueHits",    &_trueHits,    "trueHits/I");

  //Cheater info ntuple part
  RC_Tree -> Branch("betaMC",  &beta_MC,  "betaMC/D");
}

void zcl::FmmTriggerAna::analyze(const art::Event & evt)
{
  evtID    = evt.id().event();
  subRunID = evt.subRun();
  RunID    = evt.run();

  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<cheat::BackTracker> bt;
  art::ServiceHandle<calib::Calibrator> calib;


  //Get MC info first:
  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
    nxhits_MC = 0;
    nyhits_MC = 0;
    const sim::Particle *p = (*i).second;
    if (p->PdgCode()!=42 ) continue;
    const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
    unsigned ntruehits_Tot = flshits.size();

    if ( ntruehits_Tot==0 ) continue;
    double momentum2 = (p->P())*(p->P());
    double mass2     = (p->Mass())*(p->Mass());
    beta_MC = sqrt(momentum2 / (momentum2+mass2) );


    for (unsigned h = 0; h!= flshits.size(); ++h) {
      if (flshits[h].GetPlaneID() % 2) ++nxhits_MC;
      else ++nyhits_MC;
    }

    //1 monopole per event, so you can break out from the loop:
    break;
  }

  //=====================================================================================
  //Now, start pseudo trigger:
  const int _xMax       = 383; 
  const int _yMax       = 383; 
  const int _zMax       = 895; 
  const int _xdelt      = 35;
  const int _ydelt      = 35;
  const int _zdelt      = 15;

  int PX_max, PY_max, PX_min, PY_min;
  double TX_max, TY_max, TX_min, TY_min;

  art::Handle< std::vector<rb::Cluster> > slices;
  evt.getByLabel(fClusterInput, slices);

  _sliced = slices->size(); 
  //  std::cout<<"There are "<<slices->size()<<" slices in the event."<<std::endl;
  //If no slice reconstructed, just return the MC info
  if (_sliced == 0) {
    _deltaP      = 0;
    _deltaTNS    = 0;    
    _nxhits_slice = 0;
    _nyhits_slice = 0;
    _penetrated = false;
    _trueHits   = 0;
    _nSatHits   = 0;
    _meanADC    = 0;
    RC_Tree->Fill();
  }

  //Instead of splitting the hits into 2 views, just directly take the entire 3D slice:
  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
    const rb::Cluster& slice = (*slices)[sliceIdx];

    _nxhits_slice = 0;
    _nyhits_slice = 0;

    PX_max = 0; 
    PY_max = 0;
    PX_min = 999;
    PY_min = 999;
    TX_max = -9e7;      
    TY_max = -9e7;
    TX_min = 9e9;
    TY_min = 9e9;

    _penetrated = false;
    _trueHits   = 0;
    _nSatHits   = 0;
    _meanADC    = 0;

    bool Intrusion = false;
    int entry_id = -1;
    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);
      if (chit->ADC() < 501) continue;
      if (chit->ADC() > 3100) ++_nSatHits;
      _meanADC += chit->ADC();
      float tns = chit->TNS();

      //      int c = chit->Cell();
      //      int p = chit->Plane();

      if(chit->IsMC()==true) ++_trueHits;

      if (chit->View()==geo::kX) {       //X view
        ++_nxhits_slice;
        if (PX_max < chit->Plane() ) PX_max = chit->Plane();
        if (PX_min > chit->Plane() ) PX_min = chit->Plane();
        if (TX_max < tns   ) TX_max = tns;
        if (TX_min > tns   ) TX_min = tns;

        if (!Intrusion) {
          if (chit->Cell() < _xdelt){
            entry_id  = 1;
            Intrusion = true;
          }
          else if (chit->Cell() > _xMax - _xdelt){
            entry_id  = 2;
            Intrusion = true;
          }
        }
        else if (!_penetrated) {
          if      (chit->Cell() < _xdelt && entry_id!=1)       _penetrated = true;
          else if (chit->Cell() > _xMax-_xdelt && entry_id!=2) _penetrated = true;
        }
      }

      else {       //Y view
        ++_nyhits_slice;
        if (PY_max < chit->Plane() ) PY_max = chit->Plane();
        if (PY_min > chit->Plane() ) PY_min = chit->Plane();
        if (TY_max < chit->TNS()   ) TY_max = chit->TNS();
        if (TY_min > chit->TNS()   ) TY_min = chit->TNS();

        if (!Intrusion) {
          if (chit->Cell() < _ydelt) {
            entry_id  = 3;
            Intrusion = true;
          }
          else if (chit->Cell() > _yMax - _ydelt) {
            entry_id  = 4;
            Intrusion = true;
          }
        }

        else if (!_penetrated) {
          if      (chit->Cell() < _ydelt && entry_id!=3)       _penetrated = true;
          else if (chit->Cell() > _yMax-_ydelt && entry_id!=4) _penetrated = true;
        }
      }

      if (!Intrusion) {
        if (chit->Plane() < _zdelt) {
          entry_id  = 5;
          Intrusion = true;
        }
        else if (chit->Plane() > _zMax - _zdelt) {
          entry_id  = 6;
          Intrusion = true;
        }
      }
      else if (!_penetrated) {       //intrusion happened but not through yet
        if (chit->Plane() < _zdelt && entry_id!=5) _penetrated=true;
        else if (chit->Plane() > _zMax-_zdelt && entry_id!=6) _penetrated=true;
      }
    }
    //Record the info
    _meanADC     = _meanADC/(_nxhits_slice+_nyhits_slice+0.);
    _deltaP      = std::min(PX_max,PY_max) - std::max(PX_min,PY_min);
    _deltaTNS    = std::min(TX_max,TY_max) - std::max(TX_min,TY_min);    
    _tmin        = std::min(TX_min,TY_min);
    _tmax        = std::max(TX_max,TY_max);
    _pmin        = std::min(PX_min,PY_min);
    _pmax        = std::max(PX_max,PY_max);

    //   if (_nxhits_slice+_nyhits_slice < 10) _sliced = false;
    //   std::cout<<"xhits: "<<_nxhits_slice<<" yhits: "<<_nyhits_slice<<std::endl;

    // std::cout<<"There are "<<_trueHits<<" true hits in No."<<sliceIdx<<" slice in the event."<<std::endl;
    //if(_trueHits>2) RC_Tree->Fill();
    RC_Tree->Fill();
  }

} //end analyzing 


DEFINE_ART_MODULE(zcl::FmmTriggerAna)