Inheritance diagram for CTPPSProtonReconstructionEfficiencyEstimatorData:

Classes
struct	ArmData

Public Member Functions
	CTPPSProtonReconstructionEfficiencyEstimatorData (const edm::ParameterSet &)

Public Member Functions inherited from edm::one::EDAnalyzer<>
	EDAnalyzer ()=default

	EDAnalyzer (const EDAnalyzer &)=delete

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

const EDAnalyzer &	operator= (const EDAnalyzer &)=delete

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

bool	wantsInputProcessBlocks () const final

bool	wantsProcessBlocks () const final

Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDAnalyzerBase () override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

std::vector< ESProxyIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const

std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::array< std::vector< ModuleDescription const > , NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, ProductRegistry const &preg, std::map< std::string, ModuleDescription const *> const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Private Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &) override

void	endJob () override

Private Attributes
std::map< unsigned int, ArmData >	data_

std::unique_ptr< TF1 >	ff_

edm::ESGetToken< LHCInfo, LHCInfoRcd >	lhcInfoESToken_

double	localAngleXMax_

double	localAngleXMin_

double	localAngleYMax_

double	localAngleYMin_

unsigned int	n_exp_prot_max_

unsigned int	n_prep_events_

std::vector< double >	n_sigmas_

edm::ESGetToken< LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd >	opticsESToken_

edm::ESWatcher< CTPPSInterpolatedOpticsRcd >	opticsWatcher_

std::string	outputFile_

bool	pixelDiscardBXShiftedTracks_

edm::ESGetToken< PPSAssociationCuts, PPSAssociationCutsRcd >	ppsAssociationCutsToken_

edm::EDGetTokenT< reco::ForwardProtonCollection >	tokenRecoProtonsMultiRP_

edm::EDGetTokenT< CTPPSLocalTrackLiteCollection >	tokenTracks_

unsigned int	verbosity_

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

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

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()

template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)

template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

void	resetItemsToGetFrom (BranchType iType)

Detailed Description

Definition at line 43 of file CTPPSProtonReconstructionEfficiencyEstimatorData.cc.

Constructor & Destructor Documentation

◆ CTPPSProtonReconstructionEfficiencyEstimatorData()

CTPPSProtonReconstructionEfficiencyEstimatorData::CTPPSProtonReconstructionEfficiencyEstimatorData ( const edm::ParameterSet & iConfig )

explicit

Definition at line 239 of file CTPPSProtonReconstructionEfficiencyEstimatorData.cc.

References data_, edm::ParameterSet::getParameter(), n_exp_prot_max_, and n_sigmas_.

     : tokenTracks_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
 
       tokenRecoProtonsMultiRP_(
           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
 
       lhcInfoESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
       opticsESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
       ppsAssociationCutsToken_(
           esConsumes(ESInputTag("", iConfig.getParameter<std::string>("ppsAssociationCutsLabel")))),
 
       pixelDiscardBXShiftedTracks_(iConfig.getParameter<bool>("pixelDiscardBXShiftedTracks")),
 
       localAngleXMin_(iConfig.getParameter<double>("localAngleXMin")),
       localAngleXMax_(iConfig.getParameter<double>("localAngleXMax")),
       localAngleYMin_(iConfig.getParameter<double>("localAngleYMin")),
       localAngleYMax_(iConfig.getParameter<double>("localAngleYMax")),
 
       n_prep_events_(iConfig.getParameter<unsigned int>("n_prep_events")),
       n_exp_prot_max_(iConfig.getParameter<unsigned int>("n_exp_prot_max")),
       n_sigmas_(iConfig.getParameter<std::vector<double>>("n_sigmas")),
 
       outputFile_(iConfig.getParameter<string>("outputFile")),
 
       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity")),
 
       ff_(new TF1("ff", "[0] * exp(- (x-[1])*(x-[1]) / 2 / [2]/[2]) + [4]")) {
   data_[0].n_exp_prot_max = n_exp_prot_max_;
   data_[0].n_sigmas = n_sigmas_;
   data_[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
   data_[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
 
   data_[1].n_exp_prot_max = n_exp_prot_max_;
   data_[1].n_sigmas = n_sigmas_;
   data_[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
   data_[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
 }

Member Function Documentation

◆ analyze()

void CTPPSProtonReconstructionEfficiencyEstimatorData::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivatevirtual

Implements edm::one::EDAnalyzerBase.

Definition at line 319 of file CTPPSProtonReconstructionEfficiencyEstimatorData.cc.

References allShiftedPlanes, protons_cff::arm, CTPPSDetId::arm(), edm::ESWatcher< T >::check(), data_, alignCSCRings::de_x, alignCSCRings::de_y, protons_cff::decRPId, ff_, edm::EventSetup::getData(), mps_fire::i, heavyIonCSV_trainingSettings::idx, iEvent, dqmdumpme::indices, dqmiolumiharvest::j, lhcInfoESToken_, localAngleXMax_, localAngleXMin_, localAngleYMax_, localAngleYMin_, visualization-live-secondInstance_cfg::m, match(), mixedPlanes, n_prep_events_, OpticalFunctionsConfig_cfi::opticalFunctions, opticsESToken_, opticsWatcher_, AlCaHLTBitMon_ParallelJobs::p, pixelDiscardBXShiftedTracks_, HLTObjectsMonitor_cfi::plots, ppsAssociationCutsToken_, CTPPSDetId::sdTrackingPixel, CTPPSDetId::sdTrackingStrip, signature, tokenRecoProtonsMultiRP_, tokenTracks_, HLT_2022v12_cff::track, tier0::unique(), heppy_batch::val, verbosity_, x, and y.

                                                                                             {
   std::ostringstream os;
 
   // get conditions
   const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
   const auto &opticalFunctions = iSetup.getData(opticsESToken_);
   const auto &ppsAssociationCuts = iSetup.getData(ppsAssociationCutsToken_);
 
   // check optics change
   if (opticsWatcher_.check(iSetup)) {
     data_[0].UpdateOptics(opticalFunctions);
     data_[1].UpdateOptics(opticalFunctions);
   }
 
   // get input
   edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
   iEvent.getByToken(tokenTracks_, hTracks);
 
   Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
   iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
 
   // pre-selection of tracks
   std::vector<unsigned int> sel_track_idc;
 
   std::map<unsigned int, std::vector<unsigned int>> trackingSelection;
 
   for (unsigned int idx = 0; idx < hTracks->size(); ++idx) {
     const auto &tr = hTracks->at(idx);
 
     if (tr.tx() < localAngleXMin_ || tr.tx() > localAngleXMax_ || tr.ty() < localAngleYMin_ ||
         tr.ty() > localAngleYMax_)
       continue;
 
     if (pixelDiscardBXShiftedTracks_) {
       if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
           tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
         continue;
     }
 
     sel_track_idc.push_back(idx);
 
     const CTPPSDetId rpId(tr.rpId());
 
     const bool trackerRP =
         (rpId.subdetId() == CTPPSDetId::sdTrackingStrip || rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
 
     if (trackerRP)
       trackingSelection[rpId.arm()].push_back(idx);
   }
 
   // debug print
   if (verbosity_ > 1) {
     os << "* tracks (pre-selected):" << std::endl;
 
     for (const auto idx : sel_track_idc) {
       const auto &tr = hTracks->at(idx);
       CTPPSDetId rpId(tr.rpId());
       unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
       os << "    [" << idx << "] RP=" << decRPId << ", x=" << tr.x() << ", y=" << tr.y() << std::endl;
     }
 
     os << "* protons:" << std::endl;
     for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
       const auto &pr = (*hRecoProtonsMultiRP)[i];
       os << "    [" << i << "] track indices: ";
       for (const auto &trr : pr.contributingLocalTracks())
         os << trr.key() << ", ";
       os << std::endl;
     }
   }
 
   // make de_x and de_y plots
   map<unsigned int, bool> hasTracksInArm;
 
   for (const auto idx_i : sel_track_idc) {
     const auto &tr_i = hTracks->at(idx_i);
     CTPPSDetId rpId_i(tr_i.rpId());
     unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
 
     for (const auto idx_j : sel_track_idc) {
       const auto &tr_j = hTracks->at(idx_j);
       CTPPSDetId rpId_j(tr_j.rpId());
       unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
 
       // check whether desired RP combination
       unsigned int arm = 123;
       for (const auto &ap : data_) {
         if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
           arm = ap.first;
       }
 
       if (arm > 1)
         continue;
 
       // update flag
       hasTracksInArm[arm] = true;
 
       // update plots
       auto &ad = data_[arm];
       const double de_x = tr_j.x() - tr_i.x();
       const double de_y = tr_j.y() - tr_i.y();
 
       if (ad.n_events_with_tracks < n_prep_events_) {
         ad.h_de_x->Fill(de_x);
         ad.h_de_y->Fill(de_y);
       }
 
       ad.h2_de_y_vs_de_x->Fill(de_x, de_y);
     }
   }
 
   // update event counter
   for (auto &ap : data_) {
     if (hasTracksInArm[ap.first])
       ap.second.n_events_with_tracks++;
   }
 
   // if threshold reached do fits
   for (auto &ap : data_) {
     auto &ad = ap.second;
 
     if (ad.n_events_with_tracks == n_prep_events_) {
       if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.)
         continue;
 
       std::vector<std::pair<unsigned int, TH1D *>> m;
       m.emplace_back(0, ad.h_de_x.get());
       m.emplace_back(1, ad.h_de_y.get());
 
       for (const auto &p : m) {
         double max_pos = -1E100, max_val = -1.;
         for (int bi = 1; bi < p.second->GetNbinsX(); ++bi) {
           const double pos = p.second->GetBinCenter(bi);
           const double val = p.second->GetBinContent(bi);
 
           if (val > max_val) {
             max_val = val;
             max_pos = pos;
           }
         }
 
         const double sig = 0.2;
 
         ff_->SetParameters(max_val, max_pos, sig, 0.);
         p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
         p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
 
         if (p.first == 0) {
           ad.de_x_mean = ff_->GetParameter(1);
           ad.de_x_sigma = fabs(ff_->GetParameter(2));
         }
         if (p.first == 1) {
           ad.de_y_mean = ff_->GetParameter(1);
           ad.de_y_sigma = fabs(ff_->GetParameter(2));
         }
       }
 
       if (verbosity_) {
         os << "* fitting arm " << ap.first << std::endl
            << "    de_x: mean = " << ad.de_x_mean << ", sigma = " << ad.de_x_sigma << std::endl
            << "    de_y: mean = " << ad.de_y_mean << ", sigma = " << ad.de_y_sigma;
       }
     }
   }
 
   // logic below:
   //    1) evaluate "proton candiates" or "expected protons" = local tracks in the near RP which have a compatible
   //       partner in the far RP; the compatibility is evaluated at multiple n_sigma choices
   //    2) check how often each "proton candidate" can be matched to a reconstructed proton
   //    3) re-evaluate association cuts --> check details for each "proton candidate"
 
   // data structures for efficiency analysis
 
   struct ArmEventData {
     // n_sigma --> list of "proton candidates" (local track indices in the near RP)
     std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
 
     // list of valid reco-proton indices (in the chosen arm)
     std::set<unsigned int> reco_proton_idc;
 
     // n_sigma --> list of "proton candicates" which have/don't have matching reco proton
     std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
 
     // for re-evaluation of association cuts
     struct EvaluatedPair {
       unsigned idx_N, idx_F;
       bool x_cut, y_cut;
       bool match = false;
       bool unique = false;
     };
 
     vector<EvaluatedPair> evaluatedPairs;
   };
 
   std::map<unsigned int, ArmEventData> armEventData;
 
   // determine the number of proton candiates or expected protons
   for (const auto idx_i : sel_track_idc) {
     const auto &tr_i = hTracks->at(idx_i);
     CTPPSDetId rpId_i(tr_i.rpId());
     unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
 
     for (const auto idx_j : sel_track_idc) {
       const auto &tr_j = hTracks->at(idx_j);
       CTPPSDetId rpId_j(tr_j.rpId());
       unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
 
       // check whether desired RP combination
       unsigned int arm = 123;
       for (const auto &ap : data_) {
         if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
           arm = ap.first;
       }
 
       if (arm > 1)
         continue;
 
       // check near-far matching
       auto &ad = data_[arm];
       const double de_x = tr_j.x() - tr_i.x();
       const double de_y = tr_j.y() - tr_i.y();
 
       for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
         const double n_si = ad.n_sigmas[nsi];
         const bool match_x = fabs(de_x - ad.de_x_mean) < n_si * ad.de_x_sigma;
         const bool match_y = fabs(de_y - ad.de_y_mean) < n_si * ad.de_y_sigma;
         if (match_x && match_y)
           armEventData[arm].matched_track_idc[nsi].insert(idx_i);
       }
     }
   }
 
   // determine the number of reconstructed protons
   for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
     const auto &proton = (*hRecoProtonsMultiRP)[i];
 
     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
     unsigned int arm = rpId.arm();
 
     if (proton.validFit())
       armEventData[arm].reco_proton_idc.insert(i);
   }
 
   // compare matched tracks with reco protons
   if (verbosity_ > 1)
     os << "* cmp matched tracks vs. reco protons" << std::endl;
 
   for (auto &ap : armEventData) {
     auto &ad = data_[ap.first];
 
     if (verbosity_ > 1)
       os << "    arm " << ap.first << std::endl;
 
     for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
       if (verbosity_ > 1)
         os << "        nsi = " << nsi << std::endl;
 
       for (const auto &tri : ap.second.matched_track_idc[nsi]) {
         const auto &track = hTracks->at(tri);
 
         bool some_proton_matching = false;
 
         if (verbosity_ > 1)
           os << "            tri = " << tri << std::endl;
 
         for (const auto &pri : ap.second.reco_proton_idc) {
           const auto &proton = (*hRecoProtonsMultiRP)[pri];
 
           bool proton_matching = false;
 
           for (const auto &pr_tr : proton.contributingLocalTracks()) {
             CTPPSDetId rpId(pr_tr->rpId());
             unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
             if (decRPId != ad.rpId_N)
               continue;
 
             const double x = pr_tr->x();
             const double y = pr_tr->y();
             const double th = 1E-3;  // 1 um
 
             const bool match = (fabs(x - track.x()) < th) && (fabs(y - track.y()) < th);
 
             if (verbosity_ > 1)
               os << "                pri = " << pri << ": x_tr = " << track.x() << ", x_pr = " << x
                  << ", match = " << match << std::endl;
 
             if (match) {
               proton_matching = true;
               break;
             }
           }
 
           if (proton_matching) {
             some_proton_matching = true;
             break;
           }
         }
 
         if (verbosity_ > 1)
           os << "            --> some_proton_matching " << some_proton_matching << std::endl;
 
         if (some_proton_matching)
           ap.second.matched_track_with_prot_idc[nsi].insert(tri);
         else
           ap.second.matched_track_without_prot_idc[nsi].insert(tri);
       }
     }
   }
 
   // redo association cuts
   for (auto &ap : armEventData) {
     const auto &arm = ap.first;
 
     auto &ad = data_[arm];
 
     auto &aed = ap.second;
 
     auto &ass_cut = ppsAssociationCuts.getAssociationCuts(arm);
 
     const auto &indices = trackingSelection[arm];
 
     map<unsigned int, unsigned> matching_multiplicity;
 
     for (const auto &i : indices) {
       for (const auto &j : indices) {
         const auto &tr_i = hTracks->at(i);
         const auto &tr_j = hTracks->at(j);
 
         const CTPPSDetId id_i(tr_i.rpId());
         const CTPPSDetId id_j(tr_j.rpId());
 
         const unsigned rpDecId_i = id_i.arm() * 100 + id_i.station() * 10 + id_i.rp();
         const unsigned rpDecId_j = id_j.arm() * 100 + id_j.station() * 10 + id_j.rp();
 
         if (rpDecId_i != ad.rpId_N || rpDecId_j != ad.rpId_F)
           continue;
 
         ArmEventData::EvaluatedPair evp;
         evp.idx_N = i;
         evp.idx_F = j;
 
         evp.x_cut = ass_cut.isSatisfied(ass_cut.qX, tr_i.x(), tr_i.y(), lhcInfo.crossingAngle(), tr_i.x() - tr_j.x());
         evp.y_cut = ass_cut.isSatisfied(ass_cut.qY, tr_i.x(), tr_i.y(), lhcInfo.crossingAngle(), tr_i.y() - tr_j.y());
 
         evp.match = evp.x_cut && evp.y_cut;
 
         aed.evaluatedPairs.push_back(evp);
 
         if (evp.match) {
           matching_multiplicity[i]++;
           matching_multiplicity[j]++;
         }
       }
     }
 
     for (auto &evp : aed.evaluatedPairs)
       evp.unique = (matching_multiplicity[evp.idx_N] == 1) && (matching_multiplicity[evp.idx_F] == 1);
   }
 
   // debug print
   if (verbosity_ > 1) {
     for (auto &ap : armEventData) {
       auto &ad = data_[ap.first];
 
       if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
         continue;
 
       os << "* results for arm " << ap.first << std::endl;
 
       os << "    reco_proton_idc: ";
       for (const auto &idx : ap.second.reco_proton_idc)
         os << idx << ", ";
       os << std::endl;
 
       for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
         os << "    n_si = " << ad.n_sigmas[nsi] << std::endl;
 
         os << "        matched_track_idc: ";
         for (const auto &idx : ap.second.matched_track_idc[nsi])
           os << idx << ", ";
         os << std::endl;
 
         os << "        matched_track_with_prot_idc: ";
         for (const auto &idx : ap.second.matched_track_with_prot_idc[nsi])
           os << idx << ", ";
         os << std::endl;
 
         os << "        matched_track_without_prot_idc: ";
         for (const auto &idx : ap.second.matched_track_without_prot_idc[nsi])
           os << idx << ", ";
         os << std::endl;
       }
 
       os << "    evaluated pairs: ";
       for (const auto &p : ap.second.evaluatedPairs)
         os << "(" << p.idx_N << "-" << p.idx_F << ",M" << p.match << ",U" << p.unique << ")"
            << ", ";
       os << std::endl;
     }
   }
 
   // update efficiency plots
   for (auto &ap : armEventData) {
     const auto &arm = ap.first;
     auto &ad = data_[arm];
 
     // stop if sigmas not yet determined
     if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
       continue;
 
     // loop over n_sigma choices
     for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
       const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size();
       const unsigned int n_rec_prot = ap.second.reco_proton_idc.size();
 
       // stop if N(expected protons) out of range
       if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
         continue;
 
       // update method 1 plots
       const double eff = double(n_rec_prot) / n_exp_prot;
 
       for (unsigned int tri : ap.second.matched_track_idc[nsi]) {
         const double x_N = hTracks->at(tri).x();
         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
 
         ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
         ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
 
         ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
         ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
       }
 
       // update method 2 plots
       for (const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) {
         const double x_N = hTracks->at(tri).x();
         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
 
         ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
         ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
 
         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
       }
 
       for (const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) {
         const double x_N = hTracks->at(tri).x();
         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
 
         ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
         ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
 
         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
       }
 
       // update association cut plots
       for (const auto &cand : ap.second.matched_track_idc[nsi])  // loop over candidates
       {
         const double x_N = hTracks->at(cand).x();
 
         auto &plots = ad.effPlots[n_exp_prot][nsi];
 
         bool hasMatchingAndUniquePair = false;
         bool hasMatchingAndNonUniquePair = false;
         for (const auto &pair : ap.second.evaluatedPairs)  // loop over pairs
         {
           // stop if cand not compatible with pair
           if (cand != pair.idx_N)
             continue;
 
           if (pair.match && pair.unique)
             hasMatchingAndUniquePair = true;
 
           if (pair.match && !pair.unique)
             hasMatchingAndNonUniquePair = true;
 
           unsigned int signature = 999999;
           if (!pair.x_cut && !pair.y_cut)
             signature = 0;
           if (pair.x_cut && !pair.y_cut)
             signature = 1;
           if (!pair.x_cut && pair.y_cut)
             signature = 2;
           if (pair.x_cut && pair.y_cut)
             signature = 12;
 
           plots.p_fr_signature_0->Fill(x_N, (signature == 0) ? 1 : 0);
           plots.p_fr_signature_1->Fill(x_N, (signature == 1) ? 1 : 0);
           plots.p_fr_signature_2->Fill(x_N, (signature == 2) ? 1 : 0);
           plots.p_fr_signature_12->Fill(x_N, (signature == 12) ? 1 : 0);
         }
 
         plots.p_fr_match_unique->Fill(x_N, (hasMatchingAndUniquePair) ? 1 : 0);
         plots.p_fr_match_non_unique->Fill(x_N, (hasMatchingAndNonUniquePair) ? 1 : 0);
       }
     }
   }
 
   if (verbosity_)
     edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str();
 }

◆ endJob()

void CTPPSProtonReconstructionEfficiencyEstimatorData::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 827 of file CTPPSProtonReconstructionEfficiencyEstimatorData.cc.

References visDQMUpload::buf, data_, and outputFile_.

                                                               {
   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
 
   for (const auto &ait : data_) {
     char buf[100];
     sprintf(buf, "arm %u", ait.first);
     TDirectory *d_arm = f_out->mkdir(buf);
     gDirectory = d_arm;
 
     ait.second.write();
   }
 }

◆ fillDescriptions()

void CTPPSProtonReconstructionEfficiencyEstimatorData::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

Definition at line 280 of file CTPPSProtonReconstructionEfficiencyEstimatorData.cc.

References edm::ConfigurationDescriptions::add(), submitPVResolutionJobs::desc, HLT_2022v12_cff::InputTag, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                                                   {
   edm::ParameterSetDescription desc;
 
   desc.add<edm::InputTag>("tagTracks", edm::InputTag())->setComment("input tag for local lite tracks");
   desc.add<edm::InputTag>("tagRecoProtonsMultiRP", edm::InputTag())->setComment("input tag for multi-RP reco protons");
 
   desc.add<string>("lhcInfoLabel", "")->setComment("label of LHCInfo data");
   desc.add<string>("opticsLabel", "")->setComment("label of optics data");
   desc.add<string>("ppsAssociationCutsLabel", "")->setComment("label of PPSAssociationCuts data");
 
   desc.add<bool>("pixelDiscardBXShiftedTracks", false)
       ->setComment("whether to discard pixel tracks built from BX-shifted planes");
 
   desc.add<double>("localAngleXMin", -0.03)->setComment("minimal accepted value of local horizontal angle (rad)");
   desc.add<double>("localAngleXMax", +0.03)->setComment("maximal accepted value of local horizontal angle (rad)");
   desc.add<double>("localAngleYMin", -0.04)->setComment("minimal accepted value of local vertical angle (rad)");
   desc.add<double>("localAngleYMax", +0.04)->setComment("maximal accepted value of local vertical angle (rad)");
 
   desc.add<unsigned int>("n_prep_events", 1000)
       ->setComment("number of preparatory events (to determine de x and de y window)");
 
   desc.add<unsigned int>("n_exp_prot_max", 5)->setComment("maximum number of expected protons per event and per arm");
 
   desc.add<std::vector<double>>("n_sigmas", {3., 5., 7.})->setComment("list of n_sigma values");
 
   desc.add<unsigned int>("rpId_45_N", 0)->setComment("decimal RP id for 45 near");
   desc.add<unsigned int>("rpId_45_F", 0)->setComment("decimal RP id for 45 far");
   desc.add<unsigned int>("rpId_56_N", 0)->setComment("decimal RP id for 56 near");
   desc.add<unsigned int>("rpId_56_F", 0)->setComment("decimal RP id for 56 far");
 
   desc.add<std::string>("outputFile", "output.root")->setComment("output file name");
 
   desc.addUntracked<unsigned int>("verbosity", 0)->setComment("verbosity level");
 
   descriptions.add("ctppsProtonReconstructionEfficiencyEstimatorData", desc);
 }