Inheritance diagram for dune::CalibrationdEdXPDSP:

Public Member Functions
	CalibrationdEdXPDSP (fhicl::ParameterSet const &p)

	CalibrationdEdXPDSP (CalibrationdEdXPDSP const &)=delete

	CalibrationdEdXPDSP (CalibrationdEdXPDSP &&)=delete

CalibrationdEdXPDSP &	operator= (CalibrationdEdXPDSP const &)=delete

CalibrationdEdXPDSP &	operator= (CalibrationdEdXPDSP &&)=delete

void	produce (art::Event &e) override

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
std::string	fTrackModuleLabel

std::string	fCalorimetryModuleLabel

std::string	fHitModuleLabel

calo::CalorimetryAlg	caloAlg

double	fModBoxA

double	fModBoxB

bool	fSCE

bool	fApplyNormCorrection

bool	fApplyXCorrection

bool	fApplyYZCorrection

bool	fApplyLifetimeCorrection

bool	fUseLifetimeFromDatabase

std::vector< double >	fReferencedQdx

double	fLifetime

double	vDrift

double	xAnode

bool	first

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 45 of file CalibrationdEdXPDSP_module.cc.

Constructor & Destructor Documentation

dune::CalibrationdEdXPDSP::CalibrationdEdXPDSP ( fhicl::ParameterSet const & p )

explicit

Definition at line 88 of file CalibrationdEdXPDSP_module.cc.

   : EDProducer(p)
   , fTrackModuleLabel      (p.get< std::string >("TrackModuleLabel"))
   , fCalorimetryModuleLabel(p.get< std::string >("CalorimetryModuleLabel"))
   , fHitModuleLabel        (p.get< std::string >("HitModuleLabel"))
   , caloAlg                (p.get< fhicl::ParameterSet >("CaloAlg"))
   , fModBoxA               (p.get< double >("ModBoxA"))
   , fModBoxB               (p.get< double >("ModBoxB"))
   , fSCE                   (p.get< bool >("CorrectSCE"))
   , fApplyNormCorrection   (p.get< bool >("ApplyNormCorrection"))
   , fApplyXCorrection      (p.get< bool >("ApplyXCorrection"))
   , fApplyYZCorrection     (p.get< bool >("ApplyYZCorrection"))
   , fApplyLifetimeCorrection(p.get< bool >("ApplyLifetimeCorrection"))
   , fUseLifetimeFromDatabase(p.get< bool >("UseLifetimeFromDatabase"))
   , fReferencedQdx         (p.get<std::vector<double>>("ReferencedQdx"))
 {
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataForJob();
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataForJob(clockData);
   vDrift = detProp.DriftVelocity(); //cm/us
   xAnode = std::abs(detProp.ConvertTicksToX(trigger_offset(clockData),0,0,0));
   //create calorimetry product and its association with track
   produces< std::vector<anab::Calorimetry>              >();
   produces< art::Assns<recob::Track, anab::Calorimetry> >();
   first = true;
 }

dune::CalibrationdEdXPDSP::CalibrationdEdXPDSP ( CalibrationdEdXPDSP const & )

delete

dune::CalibrationdEdXPDSP::CalibrationdEdXPDSP ( CalibrationdEdXPDSP && )

delete

Member Function Documentation

CalibrationdEdXPDSP& dune::CalibrationdEdXPDSP::operator= ( CalibrationdEdXPDSP const & )

delete

CalibrationdEdXPDSP& dune::CalibrationdEdXPDSP::operator= ( CalibrationdEdXPDSP && )

delete

void dune::CalibrationdEdXPDSP::produce ( art::Event & e )

overridevirtual

Implements art::EDProducer.

Definition at line 114 of file CalibrationdEdXPDSP_module.cc.

 {
 
   art::ServiceHandle<calib::XYZCalibService> xyzcalibHandler;
   calib::XYZCalibService & xyzcalibService = *xyzcalibHandler;
   calib::XYZCalib *xyzcalib = xyzcalibService.provider();
 
   // Electron lifetime from database calibration service provider
   art::ServiceHandle<calib::LifetimeCalibService> lifetimecalibHandler;
   calib::LifetimeCalibService & lifetimecalibService = *lifetimecalibHandler;
   calib::LifetimeCalib *lifetimecalib = lifetimecalibService.provider();
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
 
   if (fApplyLifetimeCorrection) {
     if (fUseLifetimeFromDatabase) {
       fLifetime = lifetimecalib->GetLifetime()*1000.0; // [ms]*1000.0 -> [us]
       //std::cout << "use lifetime from database   " << fLifetime << std::endl;
     }
     else {
       fLifetime = detProp.ElectronLifetime(); // [us]
     }
   }
 
   //int run = evt.run();
   //int subrun = evt.subRun();
   //int event = evt.id().event();
   // evttime
   art::Timestamp ts = evt.time();
   TTimeStamp tts(ts.timeHigh(), ts.timeLow());
   uint64_t evttime = tts.AsDouble();
 
   std::cout << "run: " << evt.run() << " ; subrun: " << evt.subRun() << " ; event: " << evt.id().event() << std::endl;
   std::cout << "evttime: " << evttime << std::endl;
   if (fApplyLifetimeCorrection) std::cout << "fLifetime: " << fLifetime << " [us]" << std::endl;
   if (first){
     std::cout<<"fApplyNormCorrection: "<<fApplyNormCorrection<<std::endl;
     std::cout<<"fApplyXCorrection: "<<fApplyXCorrection<<std::endl;
     std::cout<<"fApplyYZCorrection: "<<fApplyYZCorrection<<std::endl;
     std::cout<<"fApplyLifetimeCorrection: "<<fApplyLifetimeCorrection<<std::endl;
     std::cout<<"fUseLifetimeFromDatabase: "<<fUseLifetimeFromDatabase<<std::endl;
     std::cout<<"Plane\tReference dQ/dx\tGlobal median dQ/dx\tCalorimetry constant"<<std::endl;
     for (int i = 0; i<3; ++i){
       printf("Plane[%d]\t%f\t%f\t%f\n",i,fReferencedQdx[i], xyzcalib->GetNormCorr(i), 1./caloAlg.ElectronsFromADCArea(1, i));
     }
     first = false;
   }
   //Spacecharge services provider
   auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
 
   //create anab::Calorimetry objects and make association with recob::Track
   std::unique_ptr< std::vector<anab::Calorimetry> > calorimetrycol(new std::vector<anab::Calorimetry>);
   std::unique_ptr< art::Assns<recob::Track, anab::Calorimetry> > assn(new art::Assns<recob::Track, anab::Calorimetry>);
 
   //get existing track/calorimetry objects
   auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
 
   std::vector<art::Ptr<recob::Track> > tracklist;
   art::fill_ptr_vector(tracklist, trackListHandle);
 
   art::FindManyP<anab::Calorimetry> fmcal(trackListHandle, evt, fCalorimetryModuleLabel);
 
   //get hits
   auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitModuleLabel);
 
   std::vector<art::Ptr<recob::Hit> > hitlist;
   art::fill_ptr_vector(hitlist, hitListHandle);
 
   if (!fmcal.isValid()){
     throw art::Exception(art::errors::ProductNotFound)
       <<"Could not get assocated Calorimetry objects";
   }
 
   for (size_t trkIter = 0; trkIter < tracklist.size(); ++trkIter){
     for (size_t i = 0; i<fmcal.at(trkIter).size(); ++i){
       auto & calo = fmcal.at(trkIter)[i];
 
       if (!(calo->dEdx()).size()){
         //empty calorimetry product, just copy it
         calorimetrycol->push_back(*calo);
         util::CreateAssn(*this, evt, *calorimetrycol, tracklist[trkIter], *assn);
       }
       else{
         //start calibrating dQdx
 
         //get original calorimetry information
         //double                Kin_En     = calo->KineticEnergy();
         std::vector<float>   vdEdx      = calo->dEdx();
         std::vector<float>   vdQdx      = calo->dQdx();
         std::vector<float>   vresRange  = calo->ResidualRange();
         std::vector<float>   deadwire   = calo->DeadWireResRC();
         float                Trk_Length = calo->Range();
         std::vector<float>   fpitch     = calo->TrkPitchVec();
         const auto&          fHitIndex  = calo->TpIndices();
         const auto&          vXYZ       = calo->XYZ();
         geo::PlaneID         planeID    = calo->PlaneID();
 
         //make sure the vectors are of the same size
         if (vdEdx.size()!=vXYZ.size()||
             vdQdx.size()!=vXYZ.size()||
             vresRange.size()!=vXYZ.size()||
             fpitch.size()!=vXYZ.size()){
           throw art::Exception(art::errors::Configuration)
       <<"Vector sizes mismatch for vdEdx, vdQdx, vresRange, fpitch, vXYZ";
         }
 
         //make sure the planeID is reasonable
         if (!planeID.isValid){
           throw art::Exception(art::errors::Configuration)
       <<"planeID is invalid";
         }
         if (planeID.Plane>2){
           throw art::Exception(art::errors::Configuration)
             <<"plane is invalid "<<planeID.Plane;
         }
         // update the kinetic energy
         double EkinNew = 0.;
 
         for (size_t j = 0; j<vdQdx.size(); ++j){
           auto & hit = hitlist[fHitIndex[j]];
           if (hit->WireID().Plane != planeID.Plane){
             throw art::Exception(art::errors::Configuration)
               <<"Hit plane = "<<hit->WireID().Plane<<" calo plane = "<<planeID.Plane;
           }
           double normcorrection = 1;
           if (fApplyNormCorrection){
             normcorrection = xyzcalib->GetNormCorr(planeID.Plane);
             if (normcorrection) normcorrection = fReferencedQdx[planeID.Plane]/normcorrection;
             if (!normcorrection) normcorrection = 1.;
           }
           double xcorrection = 1;
           if (fApplyXCorrection){
             xcorrection = xyzcalib->GetXCorr(planeID.Plane, vXYZ[j].X());
             if (!xcorrection) xcorrection = 1.;
           }
           double yzcorrection = 1;
           if (fApplyYZCorrection){
             yzcorrection = xyzcalib->GetYZCorr(planeID.Plane, vXYZ[j].X()>0, vXYZ[j].Y(), vXYZ[j].Z());
             if (!yzcorrection) yzcorrection = 1.;
           }
           if (fApplyLifetimeCorrection){
             xcorrection *= exp((xAnode-std::abs(vXYZ[j].X()))/(fLifetime*vDrift));
           }
           //if (planeID.Plane == 2 && tracklist[trkIter]->ID() == 0) std::cout<<"plane = "<<planeID.Plane<<" x = "<<vXYZ[j].X()<<" y = "<<vXYZ[j].Y()<<" z = "<<vXYZ[j].Z()<<" normcorrection = "<<normcorrection<<" xcorrection = "<<xcorrection<<" yzcorrection = "<<yzcorrection<<" "<<vdQdx[j]<<std::endl;
 
           vdQdx[j] = normcorrection*xcorrection*yzcorrection*vdQdx[j];
 
 
           //set time to be trgger time so we don't do lifetime correction
           //we will turn off lifetime correction in caloAlg, this is just to be double sure
           //vdEdx[j] = caloAlg.dEdx_AREA(vdQdx[j], detProp.TriggerOffset(), planeID.Plane, 0);
 
 
           //Calculate dE/dx uisng the new recombination constants
           double dQdx_e = caloAlg.ElectronsFromADCArea(vdQdx[j], planeID.Plane);
           double rho = detProp.Density();                       // LAr density in g/cm^3
           double Wion = 1000./util::kGeVToElectrons;    // 23.6 eV = 1e, Wion in MeV/e
           double E_field_nominal = detProp.Efield();   // Electric Field in the drift region in KV/cm
 
           //correct Efield for SCE
           geo::Vector_t E_field_offsets = {0., 0., 0.};
 
           if(sce->EnableCalEfieldSCE()&&fSCE) E_field_offsets = sce->GetCalEfieldOffsets(geo::Point_t{vXYZ[j].X(), vXYZ[j].Y(), vXYZ[j].Z()},hit->WireID().TPC);
 
           TVector3 E_field_vector = {E_field_nominal*(1 + E_field_offsets.X()), E_field_nominal*E_field_offsets.Y(), E_field_nominal*E_field_offsets.Z()};
           double E_field = E_field_vector.Mag();
 
           //calculate recombination factors
           double Beta = fModBoxB / (rho * E_field);
           double Alpha = fModBoxA;
           //double old_vdEdx = vdEdx[j];
           vdEdx[j] = (exp(Beta * Wion * dQdx_e) - Alpha) / Beta;
           //if (planeID.Plane == 2 && tracklist[trkIter]->ID() == 1) std::cout<<calo->dQdx()[j]<<" "<<vdEdx[j]<<" "<<vdQdx[j]<<" "<<vresRange[j]<<" "<<E_field<<" "<<rho<<" "<<Beta<<" "<<Alpha<<" "<<E_field_nominal<<" "<<E_field<<" "<<E_field_offsets.X()<<" "<<E_field_offsets.Y()<<" "<<E_field_offsets.Z()<<" "<<vXYZ[j].X()<<" "<<vXYZ[j].Y()<<" "<<vXYZ[j].Z()<<" "<<hit->WireID().TPC<<" "<<fHitIndex[j]<<" "<<normcorrection<<" "<<xcorrection<<" "<<yzcorrection<<std::endl;
 
          /*if (planeID.Plane==2){
          std::cout << sce->EnableCalEfieldSCE() << " " << fSCE << std::endl;
          std::cout << E_field << " " << E_field_nominal << std::endl;
          std::cout << vdQdx[j] << " " << dQdx_e << std::endl;
          std::cout << old_vdEdx << " " << vdEdx[j] << std::endl;
          std::cout << rho << " " << Wion << " " << Beta << "\n" << std::endl; }
          */
           //update kinetic energy calculation
           if (j>=1) {
             if ( (vresRange[j] < 0) || (vresRange[j-1] < 0) ) continue;
             EkinNew += fabs(vresRange[j]-vresRange[j-1]) * vdEdx[j];
           }
           if (j==0){
             if ( (vresRange[j] < 0) || (vresRange[j+1] < 0) ) continue;
             EkinNew += fabs(vresRange[j]-vresRange[j+1]) * vdEdx[j];
           }
 
         }
         //save new calorimetry information
         calorimetrycol->push_back(anab::Calorimetry(EkinNew,// Kin_En, // change by David C. to update kinetic energy calculation
                                                     vdEdx,
                                                     vdQdx,
                                                     vresRange,
                                                     deadwire,
                                                     Trk_Length,
                                                     fpitch,
                                                     vXYZ,
                                                     fHitIndex,
                                                     planeID));
         util::CreateAssn(*this, evt, *calorimetrycol, tracklist[trkIter], *assn);
       }//calorimetry object not empty
     }//loop over calorimetry objects
   }//loop over tracks
 
   evt.put(std::move(calorimetrycol));
   evt.put(std::move(assn));
 
   return;
 }