Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:31 

0001 #include "RecoParticleFlow/PFProducer/interface/PFCandConnector.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 using namespace reco;
0008 using namespace std;
0009 
0010 const double PFCandConnector::pion_mass2 = 0.0194;
0011 
0012 const reco::PFCandidate::Flags PFCandConnector::fT_TO_DISP_ = PFCandidate::T_TO_DISP;
0013 const reco::PFCandidate::Flags PFCandConnector::fT_FROM_DISP_ = PFCandidate::T_FROM_DISP;
0014 
0015 void PFCandConnector::setParameters(bool bCorrect,
0016                                     bool bCalibPrimary,
0017                                     double dptRel_PrimaryTrack,
0018                                     double dptRel_MergedTrack,
0019                                     double ptErrorSecondary,
0020                                     const std::vector<double>& nuclCalibFactors) {
0021   bCorrect_ = bCorrect;
0022   bCalibPrimary_ = bCalibPrimary;
0023   dptRel_PrimaryTrack_ = dptRel_PrimaryTrack;
0024   dptRel_MergedTrack_ = dptRel_MergedTrack;
0025   ptErrorSecondary_ = ptErrorSecondary;
0026 
0027   if (nuclCalibFactors.size() == 5) {
0028     fConst_[0] = nuclCalibFactors[0];
0029     fConst_[1] = nuclCalibFactors[1];
0030 
0031     fNorm_[0] = nuclCalibFactors[2];
0032     fNorm_[1] = nuclCalibFactors[3];
0033 
0034     fExp_[0] = nuclCalibFactors[4];
0035   } else {
0036     edm::LogWarning("PFCandConnector")
0037         << "Wrong calibration factors for nuclear interactions. The calibration procedure would not be applyed."
0038         << std::endl;
0039     bCalibPrimary_ = false;
0040   }
0041 
0042   std::string sCorrect = bCorrect_ ? "On" : "Off";
0043   edm::LogInfo("PFCandConnector") << " ====================== The PFCandConnector is switched " << sCorrect.c_str()
0044                                   << " ==================== " << std::endl;
0045   std::string sCalibPrimary = bCalibPrimary_ ? "used for calibration" : "not used for calibration";
0046   if (bCorrect_)
0047     edm::LogInfo("PFCandConnector") << "Primary Tracks are " << sCalibPrimary.c_str() << std::endl;
0048   if (bCorrect_ && bCalibPrimary_)
0049     edm::LogInfo("PFCandConnector") << "Under the condition that the precision on the Primary track is better than "
0050                                     << dptRel_PrimaryTrack_ << " % " << std::endl;
0051   if (bCorrect_ && bCalibPrimary_)
0052     edm::LogInfo("PFCandConnector") << "      and on merged tracks better than " << dptRel_MergedTrack_ << " % "
0053                                     << std::endl;
0054   if (bCorrect_ && bCalibPrimary_)
0055     edm::LogInfo("PFCandConnector") << "      and secondary tracks in some cases more precise than "
0056                                     << ptErrorSecondary_ << " GeV" << std::endl;
0057   if (bCorrect_ && bCalibPrimary_)
0058     edm::LogInfo("PFCandConnector") << "factor = (" << fConst_[0] << " + " << fConst_[1] << "*cFrac) - (" << fNorm_[0]
0059                                     << " - " << fNorm_[1] << "cFrac)*exp( " << -1 * fExp_[0] << "*pT )" << std::endl;
0060   edm::LogInfo("PFCandConnector") << " =========================================================== " << std::endl;
0061 }
0062 
0063 reco::PFCandidateCollection PFCandConnector::connect(PFCandidateCollection& pfCand) const {
0064   /// Collection of primary PFCandidates to be transmitted to the Event
0065   PFCandidateCollection pfC{};
0066   /// A mask to define the candidates which shall not be transmitted
0067   std::vector<bool> bMask;
0068   bMask.resize(pfCand.size(), false);
0069 
0070   // loop on primary
0071   if (bCorrect_) {
0072     LogTrace("PFCandConnector|connect") << "pfCand.size()=" << pfCand.size() << "bCalibPrimary_=" << bCalibPrimary_;
0073 
0074     for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1) {
0075       if (isPrimaryNucl(pfCand.at(ce1))) {
0076         LogTrace("PFCandConnector|connect")
0077             << "" << endl
0078             << "Nuclear Interaction w Primary Candidate " << ce1 << " " << pfCand.at(ce1) << endl
0079             << " based on the Track " << pfCand.at(ce1).trackRef().key()
0080             << " w pT = " << pfCand.at(ce1).trackRef()->pt() << " #pm "
0081             << pfCand.at(ce1).trackRef()->ptError() / pfCand.at(ce1).trackRef()->pt() * 100 << " %"
0082             << " ECAL = " << pfCand.at(ce1).ecalEnergy() << " HCAL = " << pfCand.at(ce1).hcalEnergy() << endl;
0083 
0084 #ifdef EDM_ML_DEBUG
0085         (pfCand.at(ce1)).displacedVertexRef(fT_TO_DISP_)->Dump();
0086 #endif
0087 
0088         analyseNuclearWPrim(pfCand, bMask, ce1);
0089 
0090 #ifdef EDM_ML_DEBUG
0091         LogTrace("PFCandConnector|connect")
0092             << "After Connection the candidate " << ce1 << " is " << pfCand.at(ce1) << endl
0093             << endl;
0094 
0095         PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce1).elementsInBlocks();
0096         for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
0097           if (blockElem == 0)
0098             LogTrace("PFCandConnector|connect") << *(elementsInBlocks[blockElem].first) << endl;
0099           LogTrace("PFCandConnector|connect") << " position " << elementsInBlocks[blockElem].second;
0100         }
0101 #endif
0102       }
0103     }
0104 
0105     for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1) {
0106       if (!bMask[ce1] && isSecondaryNucl(pfCand.at(ce1))) {
0107         LogTrace("PFCandConnector|connect")
0108             << "" << endl
0109             << "Nuclear Interaction w no Primary Candidate " << ce1 << " " << pfCand.at(ce1) << endl
0110             << " based on the Track " << pfCand.at(ce1).trackRef().key()
0111             << " w pT = " << pfCand.at(ce1).trackRef()->pt() << " #pm " << pfCand.at(ce1).trackRef()->ptError() << " %"
0112             << " ECAL = " << pfCand.at(ce1).ecalEnergy() << " HCAL = " << pfCand.at(ce1).hcalEnergy()
0113             << " dE(Trk-CALO) = "
0114             << pfCand.at(ce1).trackRef()->p() - pfCand.at(ce1).ecalEnergy() - pfCand.at(ce1).hcalEnergy()
0115             << " Nmissing hits = "
0116             << pfCand.at(ce1).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
0117 
0118 #ifdef EDM_ML_DEBUG
0119         (pfCand.at(ce1)).displacedVertexRef(fT_FROM_DISP_)->Dump();
0120 #endif
0121 
0122         analyseNuclearWSec(pfCand, bMask, ce1);
0123 
0124 #ifdef EDM_ML_DEBUG
0125         LogTrace("PFCandConnector|connect") << "After Connection the candidate " << ce1 << " is " << pfCand.at(ce1)
0126                                             << " and elements connected to it are: " << endl;
0127 
0128         PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce1).elementsInBlocks();
0129         for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
0130           if (blockElem == 0)
0131             LogTrace("PFCandConnector|connect") << *(elementsInBlocks[blockElem].first) << endl;
0132           LogTrace("PFCandConnector|connect") << " position " << elementsInBlocks[blockElem].second;
0133         }
0134 #endif
0135       }
0136     }
0137   }
0138 
0139   for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1)
0140     if (!bMask[ce1])
0141       pfC.push_back(pfCand.at(ce1));
0142 
0143   LogTrace("PFCandConnector|connect") << "end of function";
0144 
0145   return pfC;
0146 }
0147 
0148 void PFCandConnector::analyseNuclearWPrim(PFCandidateCollection& pfCand,
0149                                           std::vector<bool>& bMask,
0150                                           unsigned int ce1) const {
0151   PFDisplacedVertexRef ref1, ref2, ref1_bis;
0152 
0153   PFCandidate primaryCand = pfCand.at(ce1);
0154 
0155   // ------- look for the little friends -------- //
0156 
0157   const math::XYZTLorentzVectorD& momentumPrim = primaryCand.p4();
0158 
0159   math::XYZTLorentzVectorD momentumSec(0., 0., 0., 0.);
0160 
0161   if (momentumPrim.E() > 0)
0162     momentumSec = momentumPrim / momentumPrim.E() * (primaryCand.ecalEnergy() + primaryCand.hcalEnergy());
0163 
0164   map<double, math::XYZTLorentzVectorD> candidatesWithTrackExcess;
0165   map<double, math::XYZTLorentzVectorD> candidatesWithoutCalo;
0166 
0167   ref1 = primaryCand.displacedVertexRef(fT_TO_DISP_);
0168 
0169   for (unsigned int ce2 = 0; ce2 < pfCand.size(); ++ce2) {
0170     if (ce2 != ce1 && isSecondaryNucl(pfCand.at(ce2))) {
0171       ref2 = (pfCand.at(ce2)).displacedVertexRef(fT_FROM_DISP_);
0172 
0173       if (ref1 == ref2) {
0174         LogTrace("PFCandConnector|analyseNuclearWPrim")
0175             << "\t here is a Secondary Candidate " << ce2 << " " << pfCand.at(ce2) << endl
0176             << "\t based on the Track " << pfCand.at(ce2).trackRef().key()
0177             << " w p = " << pfCand.at(ce2).trackRef()->p() << " w pT = " << pfCand.at(ce2).trackRef()->pt() << " #pm "
0178             << pfCand.at(ce2).trackRef()->ptError() << " %"
0179             << " ECAL = " << pfCand.at(ce2).ecalEnergy() << " HCAL = " << pfCand.at(ce2).hcalEnergy()
0180             << " dE(Trk-CALO) = "
0181             << pfCand.at(ce2).trackRef()->p() - pfCand.at(ce2).ecalEnergy() - pfCand.at(ce2).hcalEnergy()
0182             << " Nmissing hits = "
0183             << pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
0184 
0185         if (isPrimaryNucl(pfCand.at(ce2))) {
0186           LogTrace("PFCandConnector|analyseNuclearWPrim") << "\t\t but it is also a Primary Candidate " << ce2 << endl;
0187 
0188           ref1_bis = (pfCand.at(ce2)).displacedVertexRef(fT_TO_DISP_);
0189           if (ref1_bis.isNonnull())
0190             analyseNuclearWPrim(pfCand, bMask, ce2);
0191         }
0192 
0193         // Take now the parameters of the secondary track that are relevant and use them to construct the NI candidate
0194 
0195         PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce2).elementsInBlocks();
0196         PFCandidate::ElementsInBlocks elementsAlreadyInBlocks = pfCand.at(ce1).elementsInBlocks();
0197         for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
0198           bool isAlreadyHere = false;
0199           for (unsigned alreadyBlock = 0; alreadyBlock < elementsAlreadyInBlocks.size(); alreadyBlock++) {
0200             if (elementsAlreadyInBlocks[alreadyBlock].second == elementsInBlocks[blockElem].second)
0201               isAlreadyHere = true;
0202           }
0203           if (!isAlreadyHere)
0204             pfCand.at(ce1).addElementInBlock(elementsInBlocks[blockElem].first, elementsInBlocks[blockElem].second);
0205         }
0206 
0207         double candE = pfCand.at(ce2).p4().E();
0208         double caloEn = pfCand.at(ce2).ecalEnergy() + pfCand.at(ce2).hcalEnergy();
0209         double deltaEn = candE - caloEn;
0210         int nMissOuterHits =
0211             pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0212 
0213         // Check if the difference Track Calo is not too large and if we can trust the track, ie it doesn't miss too much hits.
0214         if (deltaEn > 1 && nMissOuterHits > 1 && candE > 0) {
0215           math::XYZTLorentzVectorD momentumToAdd = pfCand.at(ce2).p4() * caloEn / candE;
0216           momentumSec += momentumToAdd;
0217           LogTrace("PFCandConnector|analyseNuclearWPrim")
0218               << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may "
0219                  "have happened. Let's trust the calo energy"
0220               << endl
0221               << "add " << momentumToAdd << endl;
0222 
0223         } else {
0224           // Check if the difference Track Calo is not too large and if we can trust the track, ie it doesn't miss too much hits.
0225           if (caloEn > 0.01 && deltaEn > 1 && nMissOuterHits > 0 && candE > 0) {
0226             math::XYZTLorentzVectorD momentumExcess = pfCand.at(ce2).p4() * deltaEn / candE;
0227             candidatesWithTrackExcess[pfCand.at(ce2).trackRef()->pt() / pfCand.at(ce2).trackRef()->ptError()] =
0228                 momentumExcess;
0229           } else if (caloEn < 0.01)
0230             candidatesWithoutCalo[pfCand.at(ce2).trackRef()->pt() / pfCand.at(ce2).trackRef()->ptError()] =
0231                 pfCand.at(ce2).p4();
0232           momentumSec += (pfCand.at(ce2)).p4();
0233         }
0234 
0235         bMask[ce2] = true;
0236       }
0237     }
0238   }
0239 
0240   // We have more primary energy than secondary: reject all secondary tracks which have no calo energy attached.
0241 
0242   if (momentumPrim.E() < momentumSec.E()) {
0243     LogTrace("PFCandConnector|analyseNuclearWPrim")
0244         << "Size of 0 calo Energy secondary candidates" << candidatesWithoutCalo.size() << endl;
0245     for (map<double, math::XYZTLorentzVectorD>::iterator iter = candidatesWithoutCalo.begin();
0246          iter != candidatesWithoutCalo.end() && momentumPrim.E() < momentumSec.E();
0247          iter++)
0248       if (momentumSec.E() > iter->second.E() + 0.1) {
0249         momentumSec -= iter->second;
0250 
0251         LogTrace("PFCandConnector|analyseNuclearWPrim")
0252             << "\t Remove a SecondaryCandidate with 0 calo energy " << iter->second << endl;
0253         LogTrace("PFCandConnector|analyseNuclearWPrim")
0254             << "momentumPrim.E() = " << momentumPrim.E() << " and momentumSec.E() = " << momentumSec.E() << endl;
0255       }
0256   }
0257 
0258   if (momentumPrim.E() < momentumSec.E()) {
0259     LogTrace("PFCandConnector|analyseNuclearWPrim")
0260         << "0 Calo Energy rejected but still not sufficient. Size of not enough calo Energy secondary candidates"
0261         << candidatesWithTrackExcess.size() << endl;
0262     for (map<double, math::XYZTLorentzVectorD>::iterator iter = candidatesWithTrackExcess.begin();
0263          iter != candidatesWithTrackExcess.end() && momentumPrim.E() < momentumSec.E();
0264          iter++)
0265       if (momentumSec.E() > iter->second.E() + 0.1)
0266         momentumSec -= iter->second;
0267   }
0268 
0269   double dpt = pfCand.at(ce1).trackRef()->ptError() / pfCand.at(ce1).trackRef()->pt() * 100;
0270 
0271   if (momentumSec.E() < 0.1) {
0272     bMask[ce1] = true;
0273     return;
0274   }
0275 
0276   // Rescale the secondary candidates to account for the loss of energy, but only if we can trust the primary track:
0277   // if it has more energy than secondaries and is precise enough and secondary exist and was not eaten or rejected during the PFAlgo step.
0278 
0279   if (((ref1->isTherePrimaryTracks() && dpt < dptRel_PrimaryTrack_) ||
0280        (ref1->isThereMergedTracks() && dpt < dptRel_MergedTrack_)) &&
0281       momentumPrim.E() > momentumSec.E() && momentumSec.E() > 0.1) {
0282     if (bCalibPrimary_) {
0283       double factor = momentumPrim.E() > 0 ? rescaleFactor(momentumPrim.Pt(), momentumSec.E() / momentumPrim.E()) : 1.;
0284       LogTrace("PFCandConnector|analyseNuclearWPrim") << "factor = " << factor << endl;
0285       if (factor * momentumPrim.Pt() < momentumSec.Pt())
0286         momentumSec = momentumPrim;
0287       else
0288         momentumSec += (1 - factor) * momentumPrim;
0289     }
0290 
0291     double px = 0, py = 0, pz = 0.;
0292     if (momentumPrim.P() > 0) {
0293       px = momentumPrim.Px() * momentumSec.P() / momentumPrim.P();
0294       py = momentumPrim.Py() * momentumSec.P() / momentumPrim.P();
0295       pz = momentumPrim.Pz() * momentumSec.P() / momentumPrim.P();
0296     }
0297     double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
0298     math::XYZTLorentzVectorD momentum(px, py, pz, E);
0299     pfCand.at(ce1).setP4(momentum);
0300 
0301     return;
0302 
0303   } else {
0304     math::XYZVector primDir = ref1->primaryDirection();
0305 
0306     if (primDir.Mag2() < 0.1) {
0307       // It might be 0 but this situation should never happend. Throw a warning if it happens.
0308       edm::LogWarning("PFCandConnector") << "A Nuclear Interaction do not have primary direction" << std::endl;
0309       pfCand.at(ce1).setP4(momentumSec);
0310       return;
0311     } else {
0312       // rescale the primary direction to the optimal momentum. But take care of the factthat it shall not be completly 0 to avoid a warning if Jet Area.
0313       double momentumS = momentumSec.P();
0314       if (momentumS < 1e-4)
0315         momentumS = 1e-4;
0316       double px = momentumS * primDir.x();
0317       double py = momentumS * primDir.y();
0318       double pz = momentumS * primDir.z();
0319       double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
0320 
0321       math::XYZTLorentzVectorD momentum(px, py, pz, E);
0322       pfCand.at(ce1).setP4(momentum);
0323       return;
0324     }
0325   }
0326 }
0327 
0328 void PFCandConnector::analyseNuclearWSec(PFCandidateCollection& pfCand,
0329                                          std::vector<bool>& bMask,
0330                                          unsigned int ce1) const {
0331   PFDisplacedVertexRef ref1, ref2;
0332 
0333   // Check if the track excess was not too large and track may miss some outer hits. This may point to a secondary NI.
0334 
0335   double caloEn = pfCand.at(ce1).ecalEnergy() + pfCand.at(ce1).hcalEnergy();
0336   double deltaEn = pfCand.at(ce1).p4().E() - caloEn;
0337   int nMissOuterHits = pfCand.at(ce1).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0338 
0339   ref1 = pfCand.at(ce1).displacedVertexRef(fT_FROM_DISP_);
0340 
0341   // ------- check if an electron or a muon vas spotted as incoming track -------- //
0342   // ------- this mean probably that the NI was fake thus we do not correct it -------- /
0343 
0344   if (ref1->isTherePrimaryTracks() || ref1->isThereMergedTracks()) {
0345     std::vector<reco::Track> refittedTracks = ref1->refittedTracks();
0346     for (unsigned it = 0; it < refittedTracks.size(); it++) {
0347       reco::TrackBaseRef primaryBaseRef = ref1->originalTrack(refittedTracks[it]);
0348       if (ref1->isIncomingTrack(primaryBaseRef))
0349         LogTrace("PFCandConnector|analyseNuclearWSec")
0350             << "There is a Primary track ref with pt = " << primaryBaseRef->pt() << endl;
0351 
0352       for (unsigned int ce = 0; ce < pfCand.size(); ++ce) {
0353         //    cout << "PFCand Id = " << (pfCand.at(ce)).particleId() << endl;
0354         if ((pfCand.at(ce)).particleId() == reco::PFCandidate::e ||
0355             (pfCand.at(ce)).particleId() == reco::PFCandidate::mu) {
0356           LogTrace("PFCandConnector|analyseNuclearWSec")
0357               << " It is an electron and it has a ref to a track " << (pfCand.at(ce)).trackRef().isNonnull() << endl;
0358 
0359           if ((pfCand.at(ce)).trackRef().isNonnull()) {
0360             reco::TrackRef tRef = (pfCand.at(ce)).trackRef();
0361             reco::TrackBaseRef bRef(tRef);
0362             LogTrace("PFCandConnector|analyseNuclearWSec")
0363                 << "With Track Ref pt = " << (pfCand.at(ce)).trackRef()->pt() << endl;
0364 
0365             if (bRef == primaryBaseRef) {
0366               if ((pfCand.at(ce)).particleId() == reco::PFCandidate::e)
0367                 LogTrace("PFCandConnector|analyseNuclearWSec")
0368                     << "It is a NI from electron. NI Discarded. Just release the candidate." << endl;
0369               if ((pfCand.at(ce)).particleId() == reco::PFCandidate::mu)
0370                 LogTrace("PFCandConnector|analyseNuclearWSec")
0371                     << "It is a NI from muon. NI Discarded. Just release the candidate" << endl;
0372 
0373               // release the track but take care of not overcounting bad tracks. In fact those tracks was protected against destruction in
0374               // PFAlgo. Now we treat them as if they was treated in PFAlgo
0375 
0376               if (caloEn < 0.1 && pfCand.at(ce1).trackRef()->ptError() > ptErrorSecondary_) {
0377                 edm::LogInfo("PFCandConnector|analyseNuclearWSec")
0378                     << "discarded track since no calo energy and ill measured" << endl;
0379                 bMask[ce1] = true;
0380               }
0381               if (caloEn > 0.1 && deltaEn > ptErrorSecondary_ &&
0382                   pfCand.at(ce1).trackRef()->ptError() > ptErrorSecondary_) {
0383                 edm::LogInfo("PFCandConnector|analyseNuclearWSec")
0384                     << "rescaled momentum of the track since no calo energy and ill measured" << endl;
0385 
0386                 double factor = caloEn / pfCand.at(ce1).p4().E();
0387                 pfCand.at(ce1).rescaleMomentum(factor);
0388               }
0389 
0390               return;
0391             }
0392           }
0393         }
0394       }
0395     }
0396   }
0397 
0398   PFCandidate secondaryCand = pfCand.at(ce1);
0399 
0400   math::XYZTLorentzVectorD momentumSec = secondaryCand.p4();
0401 
0402   if (deltaEn > ptErrorSecondary_ && nMissOuterHits > 1) {
0403     math::XYZTLorentzVectorD momentumToAdd(0., 0., 0., 0.);
0404     float candE = pfCand.at(ce1).p4().E();
0405     if (candE > 0)
0406       momentumToAdd = pfCand.at(ce1).p4() * caloEn / candE;
0407     momentumSec = momentumToAdd;
0408     LogTrace("PFCandConnector|analyseNuclearWSec")
0409         << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may have "
0410            "happened. Let's trust the calo energy"
0411         << endl
0412         << "add " << momentumToAdd << endl;
0413   }
0414 
0415   // ------- look for the little friends -------- //
0416   for (unsigned int ce2 = ce1 + 1; ce2 < pfCand.size(); ++ce2) {
0417     if (isSecondaryNucl(pfCand.at(ce2))) {
0418       ref2 = (pfCand.at(ce2)).displacedVertexRef(fT_FROM_DISP_);
0419 
0420       if (ref1 == ref2) {
0421         LogTrace("PFCandConnector|analyseNuclearWSec")
0422             << "\t here is a Secondary Candidate " << ce2 << " " << pfCand.at(ce2) << endl
0423             << "\t based on the Track " << pfCand.at(ce2).trackRef().key()
0424             << " w pT = " << pfCand.at(ce2).trackRef()->pt() << " #pm " << pfCand.at(ce2).trackRef()->ptError() << " %"
0425             << " ECAL = " << pfCand.at(ce2).ecalEnergy() << " HCAL = " << pfCand.at(ce2).hcalEnergy()
0426             << " dE(Trk-CALO) = "
0427             << pfCand.at(ce2).trackRef()->p() - pfCand.at(ce2).ecalEnergy() - pfCand.at(ce2).hcalEnergy()
0428             << " Nmissing hits = "
0429             << pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
0430 
0431         // Take now the parameters of the secondary track that are relevant and use them to construct the NI candidate
0432         PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce2).elementsInBlocks();
0433         PFCandidate::ElementsInBlocks elementsAlreadyInBlocks = pfCand.at(ce1).elementsInBlocks();
0434         for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
0435           bool isAlreadyHere = false;
0436           for (unsigned alreadyBlock = 0; alreadyBlock < elementsAlreadyInBlocks.size(); alreadyBlock++) {
0437             if (elementsAlreadyInBlocks[alreadyBlock].second == elementsInBlocks[blockElem].second)
0438               isAlreadyHere = true;
0439           }
0440           if (!isAlreadyHere)
0441             pfCand.at(ce1).addElementInBlock(elementsInBlocks[blockElem].first, elementsInBlocks[blockElem].second);
0442         }
0443 
0444         double candE = pfCand.at(ce2).p4().E();
0445         double caloEn = pfCand.at(ce2).ecalEnergy() + pfCand.at(ce2).hcalEnergy();
0446         double deltaEn = candE - caloEn;
0447         int nMissOuterHits =
0448             pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0449         if (deltaEn > ptErrorSecondary_ && nMissOuterHits > 1) {
0450           math::XYZTLorentzVectorD momentumToAdd = pfCand.at(ce2).p4() * caloEn / candE;
0451           momentumSec += momentumToAdd;
0452           LogTrace("PFCandConnector|analyseNuclearWSec")
0453               << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may "
0454                  "have happened. Let's trust the calo energy"
0455               << endl
0456               << "add " << momentumToAdd << endl;
0457         } else {
0458           momentumSec += (pfCand.at(ce2)).p4();
0459         }
0460 
0461         bMask[ce2] = true;
0462       }
0463     }
0464   }
0465 
0466   math::XYZVector primDir = ref1->primaryDirection();
0467 
0468   if (primDir.Mag2() < 0.1) {
0469     // It might be 0 but this situation should never happend. Throw a warning if it happens.
0470     pfCand.at(ce1).setP4(momentumSec);
0471     edm::LogWarning("PFCandConnector") << "A Nuclear Interaction do not have primary direction" << std::endl;
0472     return;
0473   } else {
0474     // rescale the primary direction to the optimal momentum. But take care of the factthat it shall not be completly 0 to avoid a warning if Jet Area.
0475     double momentumS = momentumSec.P();
0476     if (momentumS < 1e-4)
0477       momentumS = 1e-4;
0478     double px = momentumS * primDir.x();
0479     double py = momentumS * primDir.y();
0480     double pz = momentumS * primDir.z();
0481     double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
0482 
0483     math::XYZTLorentzVectorD momentum(px, py, pz, E);
0484 
0485     pfCand.at(ce1).setP4(momentum);
0486     return;
0487   }
0488 }
0489 
0490 bool PFCandConnector::isSecondaryNucl(const PFCandidate& pf) const {
0491   PFDisplacedVertexRef ref1;
0492   // nuclear
0493   if (pf.flag(fT_FROM_DISP_)) {
0494     ref1 = pf.displacedVertexRef(fT_FROM_DISP_);
0495     //    ref1->Dump();
0496     if (!ref1.isNonnull())
0497       return false;
0498     else if (ref1->isNucl() || ref1->isNucl_Loose() || ref1->isNucl_Kink())
0499       return true;
0500   }
0501 
0502   return false;
0503 }
0504 
0505 bool PFCandConnector::isPrimaryNucl(const PFCandidate& pf) const {
0506   PFDisplacedVertexRef ref1;
0507 
0508   // nuclear
0509   if (pf.flag(fT_TO_DISP_)) {
0510     ref1 = pf.displacedVertexRef(fT_TO_DISP_);
0511     //ref1->Dump();
0512 
0513     if (!ref1.isNonnull())
0514       return false;
0515     else if (ref1->isNucl() || ref1->isNucl_Loose() || ref1->isNucl_Kink())
0516       return true;
0517   }
0518 
0519   return false;
0520 }
0521 
0522 double PFCandConnector::rescaleFactor(const double pt, const double cFrac) const {
0523   /*
0524     LOG NORMAL FIT
0525  FCN=35.8181 FROM MIGRAD    STATUS=CONVERGED     257 CALLS         258 TOTAL
0526  EDM=8.85763e-09    STRATEGY= 1      ERROR MATRIX ACCURATE
0527   EXT PARAMETER                                   STEP         FIRST
0528   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
0529    1  p0           7.99434e-01   2.77264e-02   6.59108e-06   9.80247e-03
0530    2  p1           1.51303e-01   2.89981e-02   1.16775e-05   6.99035e-03
0531    3  p2          -5.03829e-01   2.87929e-02   1.90070e-05   1.37015e-03
0532    4  p3           4.54043e-01   5.00908e-02   3.17625e-05   3.86622e-03
0533    5  p4          -4.61736e-02   8.07940e-03   3.25775e-06  -1.37247e-02
0534   */
0535 
0536   /*
0537     FCN=34.4051 FROM MIGRAD    STATUS=CONVERGED     221 CALLS         222 TOTAL
0538     EDM=1.02201e-09    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.3 per cent
0539 
0540    fConst
0541    1  p0           7.99518e-01   2.23519e-02   1.41523e-06   4.05975e-04
0542    2  p1           1.44619e-01   2.39398e-02  -7.68117e-07  -2.55775e-03
0543 
0544    fNorm
0545    3  p2          -5.16571e-01   3.12362e-02   5.74932e-07   3.42292e-03
0546    4  p3           4.69055e-01   5.09665e-02   1.94353e-07   1.69031e-03
0547 
0548    fExp
0549    5  p4          -5.18044e-02   8.13458e-03   4.29815e-07  -1.07624e-02
0550   */
0551 
0552   double fConst, fNorm, fExp;
0553 
0554   fConst = fConst_[0] + fConst_[1] * cFrac;
0555   fNorm = fNorm_[0] - fNorm_[1] * cFrac;
0556   fExp = fExp_[0];
0557 
0558   double factor = fConst - fNorm * exp(-fExp * pt);
0559 
0560   return factor;
0561 }
0562 
0563 void PFCandConnector::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0564   iDesc.add<bool>("bCorrect", true);
0565   iDesc.add<bool>("bCalibPrimary", true);
0566   iDesc.add<double>("dptRel_PrimaryTrack", 10.0);
0567   iDesc.add<double>("dptRel_MergedTrack", 5.0);
0568   iDesc.add<double>("ptErrorSecondary", 1.0);
0569   iDesc.add<std::vector<double>>("nuclCalibFactors", {0.8, 0.15, 0.5, 0.5, 0.05});
0570 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!