Public Member Functions
	ParticleFlow (const std::string &name, framework::Process &process)

virtual void	configure (framework::config::Parameters &ps)
	Callback for the EventProcessor to configure itself from the given set of parameters.

virtual void	produce (framework::Event &event)
	Process the event and put new data products into it.

virtual void	onProcessEnd ()
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.

void	fillCandTrack (ldmx::PFCandidate &cand, const ldmx::SimTrackerHit &tk)

void	fillCandEMCalo (ldmx::PFCandidate &cand, const ldmx::CaloCluster &em)

void	fillCandHadCalo (ldmx::PFCandidate &cand, const ldmx::CaloCluster &had)

Public Member Functions inherited from framework::Producer
	Producer (const std::string &name, Process &process)
	Class constructor.

virtual void	beforeNewRun (ldmx::RunHeader &header)
	Handle allowing producers to modify run headers before the run begins.

Public Member Functions inherited from framework::EventProcessor
	EventProcessor (const std::string &name, Process &process)
	Class constructor.

virtual	~EventProcessor ()
	Class destructor.

virtual void	onNewRun (const ldmx::RunHeader &runHeader)
	Callback for the EventProcessor to take any necessary action when the run being processed changes.

virtual void	onFileOpen (EventFile &eventFile)
	Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.

virtual void	onFileClose (EventFile &eventFile)
	Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.

virtual void	onProcessStart ()
	Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.

template<class T >
const T &	getCondition (const std::string &condition_name)
	Access a conditions object for the current event.

TDirectory *	getHistoDirectory ()
	Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.

void	setStorageHint (framework::StorageControl::Hint hint)
	Mark the current event as having the given storage control hint from this module.

void	setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
	Mark the current event as having the given storage control hint from this module and the given purpose string.

int	getLogFrequency () const
	Get the current logging frequency from the process.

int	getRunNumber () const
	Get the run number from the process.

std::string	getName () const
	Get the processor name.

void	createHistograms (const std::vector< framework::config::Parameters > &histos)
	Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.

Private Attributes
TGraph *	eCorr_ {0}

TGraph *	hCorr_ {0}

std::string	inputEcalCollName_

std::string	inputHcalCollName_

std::string	inputTrackCollName_

std::string	outputCollName_

bool	singleParticle_

Additional Inherited Members
Static Public Member Functions inherited from framework::EventProcessor
static void	declare (const std::string &classname, int classtype, EventProcessorMaker *)
	Internal function which is part of the PluginFactory machinery.

Static Public Attributes inherited from framework::Producer
static const int	CLASSTYPE {1}
	Constant used to track EventProcessor types by the PluginFactory.

Protected Member Functions inherited from framework::EventProcessor
void	abortEvent ()
	Abort the event immediately.

Protected Attributes inherited from framework::EventProcessor
HistogramHelper	histograms_
	Interface class for making and filling histograms.

NtupleManager &	ntuple_ {NtupleManager::getInstance()}
	Manager for any ntuples.

logging::logger	theLog_
	The logger for this EventProcessor.

Detailed Description

Definition at line 28 of file ParticleFlow.h.

Constructor & Destructor Documentation

◆ ParticleFlow()

recon::ParticleFlow::ParticleFlow	(	const std::string &	name,
		framework::Process &	process
	)

inline

Definition at line 30 of file ParticleFlow.h.

31 : framework::Producer(name, process) {}

framework::Producer

Base class for a module which produces a data product.

Definition EventProcessor.h:252

Member Function Documentation

◆ configure()

void recon::ParticleFlow::configure ( framework::config::Parameters & parameters )

virtual

Callback for the EventProcessor to configure itself from the given set of parameters.

The parameters a processor has access to are the member variables of the python class in the sequence that has className equal to the EventProcessor class name.

For an example, look at MyProcessor.

Parameters

parameters Parameters for configuration.

Reimplemented from framework::EventProcessor.

Definition at line 7 of file ParticleFlow.cxx.

                                                          {
  // I/O
  inputEcalCollName_ = ps.getParameter<std::string>("inputEcalCollName");
  inputHcalCollName_ = ps.getParameter<std::string>("inputHcalCollName");
  inputTrackCollName_ = ps.getParameter<std::string>("inputTrackCollName");
  outputCollName_ = ps.getParameter<std::string>("outputCollName");
  // Algorithm configuration
  singleParticle_ = ps.getParameter<bool>("singleParticle");
 
  // Calibration factors, from jason, temperary
  std::vector<float> em1{250.0,  750.0,  1250.0, 1750.0, 2250.0, 2750.0,
                         3250.0, 3750.0, 4250.0, 4750.0, 5250.0, 5750.};
  std::vector<float> em2{1.175, 1.02, 0.99, 0.985, 0.975, 0.975,
                         0.96,  0.94, 0.87, 0.8,   0.73,  0.665};
  std::vector<float> h1{25.0,  75.0,  125.0, 175.0, 225.0,
                        275.0, 325.0, 375.0, 425.0};
  std::vector<float> h2{8.44,  7.38,  7.76, 8.535, 9.47,
                        10.45, 10.47, 9.71, 8.87};
  eCorr_ = new TGraph(em1.size(), em1.data(), em2.data());
  hCorr_ = new TGraph(h1.size(), h1.data(), h2.data());
}

References framework::config::Parameters::getParameter().

◆ fillCandEMCalo()

void recon::ParticleFlow::fillCandEMCalo	(	ldmx::PFCandidate &	cand,
		const ldmx::CaloCluster &	em
	)

Definition at line 53 of file ParticleFlow.cxx.

                                                             {
  float corr = 1.;
  float e = em.getEnergy();
  if (e < eCorr_->GetX()[0]) {
    corr = eCorr_->GetX()[0];
  } else if (e > eCorr_->GetX()[eCorr_->GetN() - 1]) {
    corr = eCorr_->GetX()[eCorr_->GetN() - 1];
  } else {
    corr = eCorr_->Eval(e);
  }
  cand.setEcalEnergy(e * corr);
  cand.setEcalRawEnergy(e);
  cand.setEcalClusterXYZ(em.getCentroidX(), em.getCentroidY(),
                         em.getCentroidZ());
  cand.setEcalClusterEXYZ(em.getRMSX(), em.getRMSY(), em.getRMSZ());
  cand.setEcalClusterDXDZ(em.getDXDZ());
  cand.setEcalClusterDYDZ(em.getDYDZ());
  cand.setEcalClusterEDXDZ(em.getEDXDZ());
  cand.setEcalClusterEDYDZ(em.getEDYDZ());
  cand.setPID(cand.getPID() | 2);  // OR with 010
}

◆ fillCandHadCalo()

void recon::ParticleFlow::fillCandHadCalo	(	ldmx::PFCandidate &	cand,
		const ldmx::CaloCluster &	had
	)

Definition at line 76 of file ParticleFlow.cxx.

                                                               {
  float corr = 1.;
  float e = had.getEnergy();
  if (e < hCorr_->GetX()[0]) {
    corr = hCorr_->GetX()[0];
  } else if (e > hCorr_->GetX()[hCorr_->GetN() - 1]) {
    corr = hCorr_->GetX()[hCorr_->GetN() - 1];
  } else {
    corr = hCorr_->Eval(e);
  }
  cand.setHcalEnergy(e * corr);
  cand.setHcalRawEnergy(e);
  cand.setHcalClusterXYZ(had.getCentroidX(), had.getCentroidY(),
                         had.getCentroidZ());
  cand.setHcalClusterEXYZ(had.getRMSX(), had.getRMSY(), had.getRMSZ());
  cand.setHcalClusterDXDZ(had.getDXDZ());
  cand.setHcalClusterDYDZ(had.getDYDZ());
  cand.setHcalClusterEDXDZ(had.getEDXDZ());
  cand.setHcalClusterEDYDZ(had.getEDYDZ());
  cand.setPID(cand.getPID() | 4);  // OR with 100
}

◆ fillCandTrack()

void recon::ParticleFlow::fillCandTrack	(	ldmx::PFCandidate &	cand,
		const ldmx::SimTrackerHit &	tk
	)

Definition at line 30 of file ParticleFlow.cxx.

                                                              {
  // TODO: smear
  std::vector<float> xyz = tk.getPosition();
  std::vector<double> pxyz = tk.getMomentum();
  float ecalZ = 248;
  float ecalX =
      xyz[0] + pxyz[0] / pxyz[2] * (ecalZ - xyz[2]);  // project onto ecal face
  float ecalY = xyz[1] + pxyz[1] / pxyz[2] * (ecalZ - xyz[2]);
  cand.setEcalPositionXYZ(ecalX, ecalY, ecalZ);
  cand.setTrackPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
  // also use this object to set truth info
  cand.setTruthEcalXYZ(ecalX, ecalY, ecalZ);
  cand.setTruthPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
  float m2 = pow(tk.getEnergy(), 2) - pow(pxyz[0], 2) - pow(pxyz[1], 2) -
             pow(pxyz[2], 2);
  if (m2 < 0) m2 = 0;
  cand.setTruthMass(sqrt(m2));
  cand.setTruthEnergy(tk.getEnergy());
  cand.setTruthPdgId(tk.getPdgID());
  cand.setPID(cand.getPID() | 1);  // OR with 001
}

◆ onProcessEnd()

void recon::ParticleFlow::onProcessEnd ( )

virtual

Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.

Reimplemented from framework::EventProcessor.

Definition at line 398 of file ParticleFlow.cxx.

                                {
  ldmx_log(debug) << "Process ends!";
  delete eCorr_;
  delete hCorr_;
 
  return;
}

◆ produce()

void recon::ParticleFlow::produce ( framework::Event & event )

virtual

Process the event and put new data products into it.

Parameters

event The Event to process.

Implements framework::Producer.

Definition at line 100 of file ParticleFlow.cxx.

                                                {
  if (!event.exists(inputTrackCollName_)) return;
  if (!event.exists(inputEcalCollName_)) return;
  if (!event.exists(inputHcalCollName_)) return;
  // get the track and clustering info
  const auto ecalClusters =
      event.getCollection<ldmx::CaloCluster>(inputEcalCollName_);
  const auto hcalClusters =
      event.getCollection<ldmx::CaloCluster>(inputHcalCollName_);
  const auto tracks =
      event.getCollection<ldmx::SimTrackerHit>(inputTrackCollName_);
 
  std::vector<ldmx::PFCandidate> pfCands;
  // multi-particle case
  if (!singleParticle_) {
    /*
      1. Build links maps at the Tk/Ecal and Hcal/Hcal interfaces
      2. Categorize tracks as: Ecal-matched, (Side) Hcal-matched, unmatched
      3. Categorize Ecal clusters as: Hcal-matched, unmatched
      4a. (Upstream?) Categorize tracks with dedx?
      4b. (Upstream?) Categorize ecal clusters as: EM/Had-like
      4c. (Upstream?) Categorize hcal clusters as: EM/Had-like
      5. Build candidates by category, moving from Tk-Ecal-Hcal
    */
 
    //
    // track-calo linking
    //
    std::map<int, std::vector<int> > tkCaloMap;
    std::map<int, std::vector<int> > caloTkMap;
    std::map<std::pair<int, int>, float> tkEMDist;
    // std::vector<int> unmatchedTracks;
    for (int i = 0; i < tracks.size(); i++) {
      const auto& tk = tracks[i];
      const std::vector<float> xyz = tk.getPosition();
      const std::vector<double> pxyz = tk.getMomentum();
      const float p = sqrt(pow(pxyz[0], 2) + pow(pxyz[1], 2) + pow(pxyz[2], 2));
      // float bestMatchVal = 9e9;
      for (int j = 0; j < ecalClusters.size(); j++) {
        const auto& ecal = ecalClusters[j];
        // Matching logic
        const float ecalClusZ = ecal.getCentroidZ();
        const float tkXAtClus =
            xyz[0] + pxyz[0] / pxyz[2] * (ecalClusZ - xyz[2]);  // extrapolation
        const float tkYAtClus =
            xyz[1] + pxyz[1] / pxyz[2] * (ecalClusZ - xyz[2]);
        float dist = hypot(
            (tkXAtClus - ecal.getCentroidX()) / std::max(1.0, ecal.getRMSX()),
            (tkYAtClus - ecal.getCentroidY()) / std::max(1.0, ecal.getRMSY()));
        tkEMDist[{i, j}] = dist;
        bool isMatch =
            (dist < 2) && (ecal.getEnergy() > 0.3 * p &&
                           ecal.getEnergy() < 2 * p);  // matching criteria *
        // if (isMatch && dist < bestMatchVal) bestMatchVal = dist;
        if (isMatch) {
          if (tkCaloMap.count(i))
            tkCaloMap[i].push_back(j);
          else
            tkCaloMap[i] = {j};
          if (caloTkMap.count(j))
            caloTkMap[j].push_back(i);
          else
            caloTkMap[j] = {i};
        }
      }
    }
 
    // em-hadcalo linking
    std::map<int, std::vector<int> > emHadCaloMap;
    std::map<std::pair<int, int>, float> EMHadDist;
    for (int i = 0; i < ecalClusters.size(); i++) {
      const auto& ecal = ecalClusters[i];
      for (int j = 0; j < hcalClusters.size(); j++) {
        const auto& hcal = hcalClusters[j];
        // TODO: matching logic
        const float xAtHClus =
            ecal.getCentroidX() +
            ecal.getDXDZ() * (hcal.getCentroidZ() -
                              ecal.getCentroidZ());  // extrapolated position
        const float yAtHClus =
            ecal.getCentroidY() +
            ecal.getDYDZ() * (hcal.getCentroidZ() - ecal.getCentroidZ());
        float dist = sqrt(
            pow(xAtHClus - hcal.getCentroidX(), 2) /
                std::max(1.0, pow(hcal.getRMSX(), 2) + pow(ecal.getRMSX(), 2)) +
            pow(yAtHClus - hcal.getCentroidY(), 2) /
                std::max(1.0, pow(hcal.getRMSY(), 2) + pow(ecal.getRMSY(), 2)));
        EMHadDist[{i, j}] = dist;
        bool isMatch = (dist < 5);  // matching criteria, was 2
        if (isMatch) {
          if (emHadCaloMap.count(i))
            emHadCaloMap[i].push_back(j);
          else
            emHadCaloMap[i] = {j};
        }
      }
    }
 
    // NOT YET IMPLEMENTED...
    // tk-hadcalo linking (Side HCal)
    std::map<int, std::vector<int> > tkHadCaloMap;
    // for(int i=0; i<tracks.size(); i++){
    //   const auto& tk = tracks[i];
    //   for(int j=0; j<hcalClusters.size(); j++){
    //  const auto& hcal = hcalClusters[j];
    //  // TODO: add the matching logic here...
    //  bool isMatch = true;
    //  if(isMatch){
    //    if (tkHadCaloMap.count(i)) tkHadCaloMap[i].push_back(j);
    //    else tkHadCaloMap[i] = {j};
    //  }
    //   }
    // }
 
    //
    // track / ecal cluster arbitration
    //
    std::vector<bool> tkIsEMLinked(tracks.size(), false);
    std::vector<bool> EMIsTkLinked(ecalClusters.size(), false);
    std::map<int, int> tkEMPairs{};
    for (int i = 0; i < tracks.size(); i++) {
      if (tkCaloMap.count(i)) {
        // pick first (highest-energy) unused matching cluster
        for (int em_idx : tkCaloMap[i]) {
          if (!EMIsTkLinked[em_idx]) {
            EMIsTkLinked[em_idx] = true;
            tkIsEMLinked[i] = true;
            tkEMPairs[i] = em_idx;
            break;
          }
        }
      }
    }
 
    // track / hcal cluster arbitration
    std::vector<bool> EMIsHadLinked(ecalClusters.size(), false);
    std::vector<bool> HadIsEMLinked(hcalClusters.size(), false);
    std::map<int, int> EMHadPairs{};
    for (int i = 0; i < ecalClusters.size(); i++) {
      if (emHadCaloMap.count(i)) {
        // pick first (highest-energy) unused matching cluster
        for (int had_idx : emHadCaloMap[i]) {
          if (!HadIsEMLinked[had_idx]) {
            HadIsEMLinked[had_idx] = true;
            EMIsHadLinked[i] = true;
            EMHadPairs[i] = had_idx;
            break;
          }
        }
      }
    }
 
    // can consider combining satellite clusters here...
    // define some "primary cluster" ID criterion
    //   and can add fails to the primaries
 
    //
    // Begin building pf candidates from tracks
    //
 
    // std::vector<ldmx::PFCandidate> chargedMatch;
    // std::vector<ldmx::PFCandidate> chargedUnmatch;
    for (int i = 0; i < tracks.size(); i++) {
      ldmx::PFCandidate cand;
      fillCandTrack(cand, tracks[i]);  // append track info to candidate
 
      if (!tkIsEMLinked[i]) {
        // chargedUnmatch.push_back(cand);
      } else {  // if track is linked with ECal cluster
        fillCandEMCalo(cand, ecalClusters[tkEMPairs[i]]);
        if (EMIsHadLinked[tkEMPairs[i]]) {  // if ECal is linked with HCal
                                            // cluster
          fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]);
        }
        // chargedMatch.push_back(cand);
      }
      pfCands.push_back(cand);
    }
 
    // std::vector<ldmx::PFCandidate> emMatch;
    // std::vector<ldmx::PFCandidate> emUnmatch;
    for (int i = 0; i < ecalClusters.size(); i++) {
      // already linked with ECal in the previous step
      if (EMIsTkLinked[i]) continue;
 
      ldmx::PFCandidate cand;
      fillCandEMCalo(cand, ecalClusters[i]);
      if (EMIsHadLinked[tkEMPairs[i]]) {
        fillCandHadCalo(cand, hcalClusters[EMHadPairs[i]]);
        // emMatch.push_back(cand);
      } else {
        // emUnmatch.push_back(cand);
      }
      pfCands.push_back(cand);
    }
    std::vector<ldmx::PFCandidate> hadOnly;
    for (int i = 0; i < hcalClusters.size(); i++) {
      if (HadIsEMLinked[i]) continue;
      ldmx::PFCandidate cand;
      fillCandHadCalo(cand, hcalClusters[i]);
      // hadOnly.push_back(cand);
      pfCands.push_back(cand);
    }
 
    // // track / ecal cluster arbitration
    // std::vector<ldmx::PFCandidate> caloMatchedTks;
    // std::vector<ldmx::PFCandidate> unmatchedTks;
    // std::vector<bool> emUsed (ecalClusters.size(), false);
    // for(int i=0; i<tracks.size(); i++){
    //   ldmx::PFCandidate cand;
    //   fillCandTrack(cand, tracks[i]);
    //   if( tkCaloMap.count(i)==0 ){
    //  unmatchedTks.push_back(cand);
    //   } else {
    //  // look for the first (highest-energy) unused matching cluster
    //  bool linked=false;
    //  for(int em_idx : tkCaloMap[i]){
    //    if(!emUsed[em_idx]){
    //      fillCandEMCalo(cand, tkCaloMap[i][0]);
    //      caloMatchedTks.push_back(cand);
    //      emUsed[ em_idx ] = true;
    //      linked = true;
    //      break;
    //    }
    //  }
    //  if (!linked) unmatchedTks.push_back(cand);
    //   }
    // }
 
    // // ecal / hcal cluster arbitration
    // std::vector<bool> hadUsed (hcalClusters.size(), false);
    // for(int i=0; i<ecalClusters.size(); i++){
    //   if( emHadCaloMap.count(i)==0 ){
    //  unmatchedTks.push_back(cand);
    //   } else {
    //  // look for the first (highest-energy) unused matching cluster
    //  bool linked=false;
    //  for(int em_idx : tkCaloMap[i]){
    //    if(!emUsed[em_idx]){
    //      fillCandEMCalo(cand, tkCaloMap[i][0]);
    //      caloMatchedTks.push_back(cand);
    //      emUsed[ em_idx ] = true;
    //      linked = true;
    //      break;
    //    }
    //  }
    //  if (!linked) unmatchedTks.push_back(cand);
    //   }
    // }
 
    // std::vector<ldmx::PFCandidate> unmatchedEMs;
    // for(int i=0; i<ecalClusters.size(); i++){
    //   if(emUsed[i]) continue;
    //   ldmx::PFCandidate cand;
    //   fillCandEMCalo(cand, ecalClusters[i]);
    // }
 
    //  }
    //  if( tkCaloMap[i].size()==1 ){
    //  if(!emUsed[ tkCaloMap[i][0] ]){
    //    fillCandEMCalo(cand, tkCaloMap[i][0]);
    //    caloMatchedTks.push_back(cand);
    //    emUsed[ tkCaloMap[i][0] ] = true;
    //  }
    //   }
    //   } else if( tkCaloMap[i].size()==1 ){
    //  if(!emUsed[ tkCaloMap[i][0] ]){
    //    fillCandEMCalo(cand, tkCaloMap[i][0]);
    //    caloMatchedTks.push_back(cand);
    //    emUsed[ tkCaloMap[i][0] ] = true;
    //  }
    //   }
    // }
 
  } else {
    // Single-particle builder
    ldmx::PFCandidate pf;
    int pid = 0;  // initialize pid to add
    if (tracks.size()) {
      fillCandTrack(pf, tracks[0]);
      pid += 1;
    }
    if (ecalClusters.size()) {
      fillCandEMCalo(pf, ecalClusters[0]);
      pid += 2;
    }
    if (hcalClusters.size()) {
      fillCandHadCalo(pf, hcalClusters[0]);
      pid += 4;
    }
    pf.setPID(pid);
    pf.setEnergy(pf.getEcalEnergy() + pf.getHcalEnergy());
    pfCands.push_back(pf);
  }
 
  event.add(outputCollName_, pfCands);
}

References framework::Event::exists(), ldmx::SimTrackerHit::getMomentum(), and ldmx::SimTrackerHit::getPosition().

Member Data Documentation

◆ eCorr_

TGraph* recon::ParticleFlow::eCorr_ {0}

private

Definition at line 44 of file ParticleFlow.h.

44{0};

◆ hCorr_

TGraph* recon::ParticleFlow::hCorr_ {0}

private

Definition at line 45 of file ParticleFlow.h.

45{0};

◆ inputEcalCollName_

std::string recon::ParticleFlow::inputEcalCollName_

private

Definition at line 48 of file ParticleFlow.h.

◆ inputHcalCollName_

std::string recon::ParticleFlow::inputHcalCollName_

private

Definition at line 49 of file ParticleFlow.h.

◆ inputTrackCollName_

std::string recon::ParticleFlow::inputTrackCollName_

private

Definition at line 50 of file ParticleFlow.h.

◆ outputCollName_

std::string recon::ParticleFlow::outputCollName_

private

Definition at line 52 of file ParticleFlow.h.

◆ singleParticle_

bool recon::ParticleFlow::singleParticle_

private

Definition at line 54 of file ParticleFlow.h.

The documentation for this class was generated from the following files:

Recon/include/Recon/ParticleFlow.h
Recon/src/Recon/ParticleFlow.cxx

Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ParticleFlow()

Member Function Documentation

◆ configure()

◆ fillCandEMCalo()

◆ fillCandHadCalo()

◆ fillCandTrack()

◆ onProcessEnd()

◆ produce()

Member Data Documentation

◆ eCorr_

◆ hCorr_

◆ inputEcalCollName_

◆ inputHcalCollName_

◆ inputTrackCollName_

◆ outputCollName_

◆ singleParticle_