Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-14 22:21:04 

0001 #include <memory>
0002 #include "TTree.h"
0003 #include "TFile.h"
0004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0005 
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h"
0014 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0015 // #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0016 
0017 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0018 #include "SimDataFormats/Associations/interface/MtdSimLayerClusterToTPAssociatorBaseImpl.h"
0019 #include "SimDataFormats/Associations/interface/MtdRecoClusterToSimLayerClusterAssociationMap.h"
0020 
0021 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0022 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0023 #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h"
0024 
0025 using reco::TrackCollection;
0026 
0027 class MVATrainingNtuple : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0028   typedef math::XYZTLorentzVector LorentzVector;
0029 
0030   // auxiliary class holding simulated vertices (originally from Primary4DVertexValidation)
0031   struct simPrimaryVertex {
0032     simPrimaryVertex(double x1, double y1, double z1, double t1, int k1) : x(x1), y(y1), z(z1), t(t1), key(k1) {};
0033     double x, y, z, t;
0034     int key;
0035     int eventId;
0036     int bunchCrossing;
0037     TrackingVertexRef sim_vertex;
0038     int OriginalIndex = -1;
0039     bool is_LV;
0040   };
0041 
0042 public:
0043   explicit MVATrainingNtuple(const edm::ParameterSet&);
0044   ~MVATrainingNtuple() override;
0045 
0046   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047 
0048 private:
0049   void analyze(const edm::Event&, const edm::EventSetup&) override;
0050 
0051   const bool trkTPSelAll(const TrackingParticle&);
0052   const bool trkRecSel(const reco::TrackBase&);
0053 
0054   const edm::Ref<std::vector<TrackingParticle>>* getAnyMatchedTP(const reco::TrackBaseRef&);
0055   double timeFromTrueMass(double, double, double, double);
0056 
0057   std::vector<MVATrainingNtuple::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
0058 
0059   edm::Service<TFileService> fs_;
0060 
0061   // ----------member data ---------------------------
0062   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> theTTBToken;
0063   TrackFilterForPVFindingBase* theTrackFilter;
0064 
0065   static constexpr double simUnit_ = 1e9;     //sim time in s while reco time in ns
0066   static constexpr double c_ = 2.99792458e1;  //c in cm/ns
0067   static constexpr double etacutGEN_ = 4.;    // |eta| < 4;
0068   static constexpr double etacutREC_ = 3.;    // |eta| < 3;
0069   static constexpr double etacutBTL_ = 1.5;   // |eta| < 1.5;
0070   static constexpr double pTcutBTL_ = 0.7;    // PT > 0.7 GeV
0071   static constexpr double pTcutETL_ = 0.2;    // PT > 0.2 GeV
0072   static constexpr double rBTL_ = 110.0;
0073   static constexpr double zETL_ = 290.0;
0074   std::string fileName_;
0075 
0076   bool saveNtupleforBDT_;
0077   bool saveNtupleforGNN_;
0078 
0079   const reco::RecoToSimCollection* r2s_;
0080 
0081   // GNN input variables
0082   std::vector<double> gnn_pt, gnn_eta, gnn_phi, gnn_z_pca, gnn_dz, gnn_t_Pi, gnn_t_K, gnn_t_P, gnn_t0safe, gnn_t0pid,
0083       gnn_sigma_t0safe, gnn_mtdTime, gnn_sigma_tmtd, gnn_mva_qual, gnn_btlMatchChi2, gnn_btlMatchTimeChi2,
0084       gnn_etlMatchChi2, gnn_etlMatchTimeChi2, gnn_pathLength, gnn_probPi, gnn_probK, gnn_probP, gnn_trk_chi2,
0085       gnn_trk_ndof, gnn_sigma_tof_Pi, gnn_sigma_tof_K, gnn_sigma_tof_P, gnn_sim_vertex_z, gnn_sim_vertex_t, gnn_tp_tEst,
0086       gnn_outermostHitPosition;
0087   std::vector<int> gnn_npixBarrel, gnn_npixEndcap, gnn_sim_vertex_evID, gnn_sim_vertex_BX, gnn_sim_vertex_index,
0088       gnn_tp_pdgId, gnn_trk_validhits;
0089   std::vector<bool> gnn_is_matched_tp, gnn_sim_vertex_isLV;
0090 
0091   // BDT input variables
0092   std::vector<double> Ttrack_pt, Ttrack_eta, Ttrack_phi, Ttrack_dz, Ttrack_dxy, Ttrack_chi2, Ttrack_BTLchi2,
0093       Ttrack_BTLtime_chi2, Ttrack_ETLchi2, Ttrack_ETLtime_chi2, Ttrack_t0, Ttrack_sigmat0, Ttrack_Tmtd,
0094       Ttrack_sigmaTmtd, Ttrack_length, Ttrack_MtdMVA, Ttrack_lHitPos;
0095   std::vector<int> Ttrack_ndof, Ttrack_nValidHits, Ttrack_npixBarrelValidHits, Ttrack_npixEndcapValidHits;
0096   std::vector<bool> Ttrack_Signal, Ttrack_Associated, Ttrack_TPDirClu, Ttrack_TPOtherClu, Ttrack_isBTL, Ttrack_isETL,
0097       Ttrack_withMTD;
0098 
0099   edm::EDGetTokenT<edm::ValueMap<float>> btlMatchChi2Token_;
0100   edm::EDGetTokenT<edm::ValueMap<float>> btlMatchTimeChi2Token_;
0101   edm::EDGetTokenT<edm::ValueMap<float>> etlMatchChi2Token_;
0102   edm::EDGetTokenT<edm::ValueMap<float>> etlMatchTimeChi2Token_;
0103   edm::EDGetTokenT<edm::ValueMap<int>> npixBarrelToken_;
0104   edm::EDGetTokenT<edm::ValueMap<int>> npixEndcapToken_;
0105   edm::EDGetTokenT<edm::ValueMap<float>> outermostHitPositionToken_;
0106 
0107   edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
0108   edm::EDGetTokenT<reco::TrackCollection> RecMTDTrackToken_;
0109   edm::EDGetTokenT<std::vector<reco::Vertex>> RecVertexToken_;
0110   edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
0111   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
0112   edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
0113   edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
0114   edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
0115   edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
0116   edm::EDGetTokenT<reco::TPToSimCollectionMtd> tp2SimAssociationMapToken_;
0117   edm::EDGetTokenT<MtdRecoClusterToSimLayerClusterAssociationMap> r2sAssociationMapToken_;
0118   edm::EDGetTokenT<FTLClusterCollection> btlRecCluToken_;
0119   edm::EDGetTokenT<FTLClusterCollection> etlRecCluToken_;
0120 
0121   edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
0122   edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
0123   edm::EDGetTokenT<edm::ValueMap<float>> sigmatimeToken_;
0124   edm::EDGetTokenT<edm::ValueMap<float>> t0SrcToken_;
0125   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SrcToken_;
0126   edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
0127   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0PidToken_;
0128   edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
0129   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
0130   edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
0131   edm::EDGetTokenT<edm::ValueMap<float>> sigmatofpiToken_;
0132   edm::EDGetTokenT<edm::ValueMap<float>> sigmatofkToken_;
0133   edm::EDGetTokenT<edm::ValueMap<float>> sigmatofpToken_;
0134   edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
0135   edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
0136   edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
0137   edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
0138   edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
0139   edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
0140   edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
0141 };
0142 
0143 MVATrainingNtuple::MVATrainingNtuple(const edm::ParameterSet& iConfig)
0144     : theTTBToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0145       fileName_(iConfig.getUntrackedParameter<std::string>("fileName")) {
0146   RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTracks"));
0147   RecMTDTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTagT"));
0148   RecVertexToken_ = consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("inputTagV"));
0149   tp2SimAssociationMapToken_ =
0150       consumes<reco::TPToSimCollectionMtd>(iConfig.getParameter<edm::InputTag>("tp2SimAssociationMapTag"));
0151   r2sAssociationMapToken_ = consumes<MtdRecoClusterToSimLayerClusterAssociationMap>(
0152       iConfig.getParameter<edm::InputTag>("r2sAssociationMapTag"));
0153   trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
0154   RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
0155   trackingParticleCollectionToken_ =
0156       consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0157   recoToSimAssociationToken_ =
0158       consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0159   simToRecoAssociationToken_ =
0160       consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0161   trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0162   btlRecCluToken_ = consumes<FTLClusterCollection>(iConfig.getParameter<edm::InputTag>("recCluTagBTL"));
0163   etlRecCluToken_ = consumes<FTLClusterCollection>(iConfig.getParameter<edm::InputTag>("recCluTagETL"));
0164   pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
0165   momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
0166   sigmatimeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmaSrc"));
0167   t0SrcToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0Src"));
0168   sigmat0SrcToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0Src"));
0169   t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
0170   sigmat0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0PID"));
0171   t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
0172   sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
0173   trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
0174   tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
0175   tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
0176   tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
0177   tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
0178   probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
0179   probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
0180   probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
0181   sigmatofpiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofpiSrc"));
0182   sigmatofkToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofkSrc"));
0183   sigmatofpToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofpSrc"));
0184   btlMatchChi2Token_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("btlMatchChi2Src"));
0185   btlMatchTimeChi2Token_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("btlMatchTimeChi2Src"));
0186   etlMatchChi2Token_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("etlMatchChi2Src"));
0187   etlMatchTimeChi2Token_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("etlMatchTimeChi2Src"));
0188   npixBarrelToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("npixBarrelSrc"));
0189   npixEndcapToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("npixEndcapSrc"));
0190   outermostHitPositionToken_ =
0191       consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("outermostHitPositionSrc"));
0192   saveNtupleforBDT_ = iConfig.getParameter<bool>("ntupleforBDT");
0193   saveNtupleforGNN_ = iConfig.getParameter<bool>("ntupleforGNN");
0194   // select and configure the track selection
0195   std::string trackSelectionAlgorithm =
0196       iConfig.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
0197   if (trackSelectionAlgorithm == "filter") {
0198     theTrackFilter = new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0199   } else if (trackSelectionAlgorithm == "filterWithThreshold") {
0200     theTrackFilter = new HITrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0201   } else {
0202     edm::LogWarning("MVATrainingNtuple: unknown track selection algorithm: " + trackSelectionAlgorithm);
0203   }
0204 }
0205 
0206 MVATrainingNtuple::~MVATrainingNtuple() {
0207   if (theTrackFilter)
0208     delete theTrackFilter;
0209 }
0210 
0211 const edm::Ref<std::vector<TrackingParticle>>* MVATrainingNtuple::getAnyMatchedTP(const reco::TrackBaseRef& recoTrack) {
0212   auto found = r2s_->find(recoTrack);
0213 
0214   // reco track not matched to any TP
0215   if (found == r2s_->end())
0216     return nullptr;
0217 
0218   //matched TP equal to any TP
0219   for (const auto& tp : found->val) {
0220     return &tp.first;
0221   }
0222 
0223   // reco track not matched to any TP from vertex
0224   return nullptr;
0225 }
0226 
0227 double MVATrainingNtuple::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
0228   if (time > 0 && pathlength > 0 && mass > 0) {
0229     double gammasq = 1. + momentum * momentum / (mass * mass);
0230     double v = c_ * std::sqrt(1. - 1. / gammasq);  // cm / ns
0231     double t_est = time - (pathlength / v);
0232 
0233     return t_est;
0234   } else {
0235     return -1;
0236   }
0237 }
0238 
0239 std::vector<MVATrainingNtuple::simPrimaryVertex> MVATrainingNtuple::getSimPVs(
0240     const edm::Handle<TrackingVertexCollection>& tVC) {
0241   std::vector<MVATrainingNtuple::simPrimaryVertex> simpv;
0242 
0243   int current_event = -1;
0244   int s = -1;
0245   for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
0246     // LV is the first vertex in each event, keep only at BX=0
0247     int eventId = v->eventId().event();
0248     int bunchCrossing = v->eventId().bunchCrossing();
0249 
0250     if (bunchCrossing != 0)
0251       continue;
0252 
0253     bool is_LV = true;
0254     if (eventId != current_event) {
0255       current_event = eventId;
0256     } else {
0257       is_LV = false;
0258     }
0259     s++;
0260 
0261     // could be a new vertex, check  all primaries found so far to avoid multiple entries
0262     int key = std::distance(tVC->begin(), v);
0263     simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t(), key);
0264     sv.eventId = eventId;
0265     sv.bunchCrossing = bunchCrossing;
0266     sv.sim_vertex = TrackingVertexRef(tVC, key);
0267     sv.OriginalIndex = s;
0268     sv.is_LV = is_LV;
0269 
0270     simPrimaryVertex* vp = nullptr;  // will become non-NULL if a vertex is found and then point to it
0271     for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
0272       if ((sv.eventId == v0->eventId) && (sv.bunchCrossing == v0->bunchCrossing) && (std::abs(sv.x - v0->x) < 1e-5) &&
0273           (std::abs(sv.y - v0->y) < 1e-5) && (std::abs(sv.z - v0->z) < 1e-5)) {
0274         vp = &(*v0);
0275         break;
0276       }
0277     }
0278     if (!vp) {
0279       // this is a new vertex, add it to the list of sim-vertices
0280       simpv.push_back(sv);
0281     }
0282 
0283   }  // End of for loop on tracking vertices
0284 
0285   // In case of no simulated vertices, break here
0286   if (simpv.empty())
0287     return simpv;
0288 
0289   return simpv;
0290 }
0291 
0292 // ------------ method called for each event  ------------
0293 void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0294   using namespace edm;
0295   using namespace std;
0296   using edm::Handle;
0297   using std::vector;
0298   using namespace reco;
0299 
0300   edm::Handle<reco::TrackCollection> tracksH;
0301   iEvent.getByToken(RecTrackToken_, tracksH);
0302 
0303   const auto& theB = &iSetup.getData(theTTBToken);
0304   std::vector<reco::TransientTrack> t_tks;
0305 
0306   edm::Handle<TrackingParticleCollection> TPCollectionH;
0307   iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
0308   if (!TPCollectionH.isValid())
0309     edm::LogWarning("MVATrainingNtuple") << "TPCollectionH is not valid";
0310 
0311   edm::Handle<reco::RecoToSimCollection> recoToSimH;
0312   iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
0313   if (recoToSimH.isValid())
0314     r2s_ = recoToSimH.product();
0315   else
0316     edm::LogWarning("MVATrainingNtuple") << "recoToSimH is not valid";
0317 
0318   reco::BeamSpot beamSpot;
0319   edm::Handle<reco::BeamSpot> BeamSpotH;
0320   iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
0321   if (!BeamSpotH.isValid())
0322     edm::LogWarning("MVATrainingNtuple") << "BeamSpotH is not valid";
0323   beamSpot = *BeamSpotH;
0324 
0325   edm::Handle<TrackingVertexCollection> TVCollectionH;
0326   iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
0327   if (!TVCollectionH.isValid())
0328     edm::LogWarning("MVATrainingNtuple") << "TVCollectionH is not valid";
0329 
0330   std::vector<simPrimaryVertex> simpv;
0331   simpv = getSimPVs(TVCollectionH);
0332 
0333   const auto& trackAssoc = iEvent.get(trackAssocToken_);
0334 
0335   const auto& tp2SimAssociationMap = iEvent.get(tp2SimAssociationMapToken_);
0336   const auto& r2sAssociationMap = iEvent.get(r2sAssociationMapToken_);
0337 
0338   // auto RecTrackHandle = makeValid(iEvent.getHandle(RecTrackToken_));
0339   auto RecTrackHandle = iEvent.getHandle(RecTrackToken_);
0340 
0341   const auto& btlRecCluHandle = iEvent.getHandle(btlRecCluToken_);
0342   const auto& etlRecCluHandle = iEvent.getHandle(etlRecCluToken_);
0343 
0344   const auto& pathLength = iEvent.get(pathLengthToken_);
0345   const auto& momentum = iEvent.get(momentumToken_);
0346   const auto& sigmatimemtd = iEvent.get(sigmatimeToken_);
0347   const auto& t0Src = iEvent.get(t0SrcToken_);
0348   const auto& sigmat0Src = iEvent.get(sigmat0SrcToken_);
0349   const auto& t0Pid = iEvent.get(t0PidToken_);
0350   const auto& sigmat0Pid = iEvent.get(sigmat0PidToken_);
0351   const auto& t0Safe = iEvent.get(t0SafePidToken_);
0352   const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
0353   const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
0354   const auto& tMtd = iEvent.get(tmtdToken_);
0355   const auto& tofPi = iEvent.get(tofPiToken_);
0356   const auto& tofK = iEvent.get(tofKToken_);
0357   const auto& tofP = iEvent.get(tofPToken_);
0358   const auto& probPi = iEvent.get(probPiToken_);
0359   const auto& probK = iEvent.get(probKToken_);
0360   const auto& probP = iEvent.get(probPToken_);
0361   const auto& sigmatofpi = iEvent.get(sigmatofpiToken_);
0362   const auto& sigmatofk = iEvent.get(sigmatofkToken_);
0363   const auto& sigmatofp = iEvent.get(sigmatofpToken_);
0364   const auto& btlMatchChi2 = iEvent.get(btlMatchChi2Token_);
0365   const auto& btlMatchTimeChi2 = iEvent.get(btlMatchTimeChi2Token_);
0366   const auto& etlMatchChi2 = iEvent.get(etlMatchChi2Token_);
0367   const auto& etlMatchTimeChi2 = iEvent.get(etlMatchTimeChi2Token_);
0368   const auto& npixBarrel = iEvent.get(npixBarrelToken_);
0369   const auto& npixEndcap = iEvent.get(npixEndcapToken_);
0370   const auto& outermostHitPosition = iEvent.get(outermostHitPositionToken_);
0371 
0372   // Fill TTree with input variables for GNN
0373   if (saveNtupleforGNN_) {
0374     std::string GNNtreeName = "GNNtree_" + std::to_string(iEvent.id().event());
0375     TTree* GNNtree = fs_->make<TTree>(GNNtreeName.c_str(), "Tree for GNN tracks");
0376 
0377     GNNtree->Branch("gnn_pt", &gnn_pt);
0378     GNNtree->Branch("gnn_eta", &gnn_eta);
0379     GNNtree->Branch("gnn_phi", &gnn_phi);
0380     GNNtree->Branch("gnn_z_pca", &gnn_z_pca);
0381     GNNtree->Branch("gnn_dz", &gnn_dz);
0382     GNNtree->Branch("gnn_t_Pi", &gnn_t_Pi);
0383     GNNtree->Branch("gnn_t_K", &gnn_t_K);
0384     GNNtree->Branch("gnn_t_P", &gnn_t_P);
0385     GNNtree->Branch("gnn_sigma_t0safe", &gnn_sigma_t0safe);
0386     GNNtree->Branch("gnn_sigma_tmtd", &gnn_sigma_tmtd);
0387     GNNtree->Branch("gnn_t0safe", &gnn_t0safe);
0388     GNNtree->Branch("gnn_t0pid", &gnn_t0pid);
0389     GNNtree->Branch("gnn_mva_qual", &gnn_mva_qual);
0390     GNNtree->Branch("gnn_btlMatchChi2", &gnn_btlMatchChi2);
0391     GNNtree->Branch("gnn_btlMatchTimeChi2", &gnn_btlMatchTimeChi2);
0392     GNNtree->Branch("gnn_etlMatchChi2", &gnn_etlMatchChi2);
0393     GNNtree->Branch("gnn_etlMatchTimeChi2", &gnn_etlMatchTimeChi2);
0394     GNNtree->Branch("gnn_pathLength", &gnn_pathLength);
0395     GNNtree->Branch("gnn_npixBarrel", &gnn_npixBarrel);
0396     GNNtree->Branch("gnn_npixEndcap", &gnn_npixEndcap);
0397     GNNtree->Branch("gnn_outermostHitPosition", &gnn_outermostHitPosition);
0398     GNNtree->Branch("gnn_mtdTime", &gnn_mtdTime);
0399     GNNtree->Branch("gnn_is_matched_tp", &gnn_is_matched_tp);
0400     GNNtree->Branch("gnn_tp_tEst", &gnn_tp_tEst);
0401     GNNtree->Branch("gnn_tp_pdgId", &gnn_tp_pdgId);
0402     GNNtree->Branch("gnn_probPi", &gnn_probPi);
0403     GNNtree->Branch("gnn_probK", &gnn_probK);
0404     GNNtree->Branch("gnn_probP", &gnn_probP);
0405     GNNtree->Branch("gnn_sigma_tof_Pi", &gnn_sigma_tof_Pi);
0406     GNNtree->Branch("gnn_sigma_tof_K", &gnn_sigma_tof_K);
0407     GNNtree->Branch("gnn_sigma_tof_P", &gnn_sigma_tof_P);
0408     GNNtree->Branch("gnn_trk_chi2", &gnn_trk_chi2);
0409     GNNtree->Branch("gnn_trk_ndof", &gnn_trk_ndof);
0410     GNNtree->Branch("gnn_trk_validhits", &gnn_trk_validhits);
0411     GNNtree->Branch("gnn_sim_vertex_evID", &gnn_sim_vertex_evID);
0412     GNNtree->Branch("gnn_sim_vertex_BX", &gnn_sim_vertex_BX);
0413     GNNtree->Branch("gnn_sim_vertex_index", &gnn_sim_vertex_index);
0414     GNNtree->Branch("gnn_sim_vertex_z", &gnn_sim_vertex_z);
0415     GNNtree->Branch("gnn_sim_vertex_t", &gnn_sim_vertex_t);
0416     GNNtree->Branch("gnn_sim_vertex_isLV", &gnn_sim_vertex_isLV);
0417 
0418     gnn_pt.clear();
0419     gnn_eta.clear();
0420     gnn_phi.clear();
0421     gnn_z_pca.clear();
0422     gnn_dz.clear();
0423     gnn_t_Pi.clear();
0424     gnn_t_K.clear();
0425     gnn_t_P.clear();
0426     gnn_sigma_t0safe.clear();
0427     gnn_sigma_tmtd.clear();
0428     gnn_t0safe.clear();
0429     gnn_t0pid.clear();
0430     gnn_mva_qual.clear();
0431     gnn_btlMatchChi2.clear();
0432     gnn_btlMatchTimeChi2.clear();
0433     gnn_etlMatchChi2.clear();
0434     gnn_etlMatchTimeChi2.clear();
0435     gnn_pathLength.clear();
0436     gnn_npixBarrel.clear();
0437     gnn_npixEndcap.clear();
0438     gnn_outermostHitPosition.clear();
0439     gnn_mtdTime.clear();
0440     gnn_is_matched_tp.clear();
0441     gnn_tp_tEst.clear();
0442     gnn_tp_pdgId.clear();
0443     gnn_probPi.clear();
0444     gnn_probK.clear();
0445     gnn_probP.clear();
0446     gnn_sigma_tof_Pi.clear();
0447     gnn_sigma_tof_K.clear();
0448     gnn_sigma_tof_P.clear();
0449     gnn_trk_chi2.clear();
0450     gnn_trk_ndof.clear();
0451     gnn_trk_validhits.clear();
0452     gnn_sim_vertex_evID.clear();
0453     gnn_sim_vertex_BX.clear();
0454     gnn_sim_vertex_index.clear();
0455     gnn_sim_vertex_z.clear();
0456     gnn_sim_vertex_t.clear();
0457     gnn_sim_vertex_isLV.clear();
0458 
0459     // build TransientTracks
0460     t_tks = (*theB).build(tracksH, beamSpot, t0Safe, sigmat0Safe);
0461 
0462     // track filter
0463     std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
0464 
0465     for (std::vector<reco::TransientTrack>::const_iterator itk = seltks.begin(); itk != seltks.end(); itk++) {
0466       reco::TrackBaseRef trackref = (*itk).trackBaseRef();
0467 
0468       gnn_pt.push_back((*itk).track().pt());
0469       gnn_eta.push_back((*itk).track().eta());
0470       gnn_phi.push_back((*itk).track().phi());
0471       gnn_z_pca.push_back((*itk).track().vz());
0472       gnn_dz.push_back((*itk).track().dzError());
0473       gnn_t_Pi.push_back(tMtd[trackref] - tofPi[trackref]);
0474       gnn_t_K.push_back(tMtd[trackref] - tofK[trackref]);
0475       gnn_t_P.push_back(tMtd[trackref] - tofP[trackref]);
0476       gnn_sigma_t0safe.push_back(sigmat0Safe[trackref]);
0477       gnn_sigma_tmtd.push_back(sigmatimemtd[trackref]);
0478       gnn_t0safe.push_back(t0Safe[trackref]);
0479       gnn_t0pid.push_back(t0Pid[trackref]);
0480       gnn_mva_qual.push_back(mtdQualMVA[trackref]);
0481       gnn_btlMatchChi2.push_back(btlMatchChi2[trackref]);
0482       gnn_btlMatchTimeChi2.push_back(btlMatchTimeChi2[trackref]);
0483       gnn_etlMatchChi2.push_back(etlMatchChi2[trackref]);
0484       gnn_etlMatchTimeChi2.push_back(etlMatchTimeChi2[trackref]);
0485       gnn_pathLength.push_back(pathLength[trackref]);
0486       gnn_npixBarrel.push_back(npixBarrel[trackref]);
0487       gnn_npixEndcap.push_back(npixEndcap[trackref]);
0488       gnn_outermostHitPosition.push_back(outermostHitPosition[trackref]);
0489       gnn_mtdTime.push_back(tMtd[trackref]);
0490       gnn_probPi.push_back(probPi[trackref]);
0491       gnn_probK.push_back(probK[trackref]);
0492       gnn_probP.push_back(probP[trackref]);
0493       gnn_sigma_tof_Pi.push_back(sigmatofpi[trackref]);
0494       gnn_sigma_tof_K.push_back(sigmatofk[trackref]);
0495       gnn_sigma_tof_P.push_back(sigmatofp[trackref]);
0496       gnn_trk_chi2.push_back((*itk).track().chi2());
0497       gnn_trk_ndof.push_back((*itk).track().ndof());
0498       gnn_trk_validhits.push_back((*itk).track().numberOfValidHits());
0499 
0500       auto anytp_info = getAnyMatchedTP(trackref);
0501       if (anytp_info != nullptr) {
0502         gnn_is_matched_tp.push_back(true);
0503         double anytp_mass = (*anytp_info)->mass();
0504         gnn_tp_tEst.push_back(timeFromTrueMass(anytp_mass, pathLength[trackref], momentum[trackref], tMtd[trackref]));
0505         gnn_tp_pdgId.push_back(std::abs((*anytp_info)->pdgId()));
0506 
0507         TrackingVertexRef parentVertexRef = (*anytp_info)->parentVertex();
0508 
0509         // Loop on TV Collection to retrive info on sim vertices
0510         bool vertex_match = false;
0511         for (const auto& vsim : simpv) {
0512           if (vsim.sim_vertex == parentVertexRef) {
0513             vertex_match = true;
0514             // Found the matching simPrimaryVertex
0515             gnn_sim_vertex_z.push_back(vsim.z);
0516             gnn_sim_vertex_t.push_back(vsim.t * simUnit_);
0517             gnn_sim_vertex_evID.push_back(vsim.eventId);
0518             gnn_sim_vertex_BX.push_back(vsim.bunchCrossing);
0519             gnn_sim_vertex_index.push_back(vsim.key);
0520             gnn_sim_vertex_isLV.push_back(vsim.is_LV);
0521           }
0522         }
0523         if (vertex_match == false) {
0524           gnn_sim_vertex_z.push_back(-999.);
0525           gnn_sim_vertex_t.push_back(-999.);
0526           gnn_sim_vertex_evID.push_back(-999);
0527           gnn_sim_vertex_BX.push_back(-999);
0528           gnn_sim_vertex_index.push_back(-999);
0529           gnn_sim_vertex_isLV.push_back(false);
0530         }
0531 
0532       } else {
0533         gnn_is_matched_tp.push_back(false);
0534         gnn_tp_tEst.push_back(-999.);
0535         gnn_tp_pdgId.push_back(-999);
0536         gnn_sim_vertex_z.push_back(-999.);
0537         gnn_sim_vertex_t.push_back(-999.);
0538         gnn_sim_vertex_evID.push_back(-999);
0539         gnn_sim_vertex_BX.push_back(-999);
0540         gnn_sim_vertex_index.push_back(-999);
0541         gnn_sim_vertex_isLV.push_back(false);
0542       }
0543 
0544     }  // loop on sel tracks
0545 
0546     GNNtree->Fill();
0547 
0548   }  // ntuple for GNN
0549 
0550   // Fill TTree with input variables for BDT
0551   if (saveNtupleforBDT_) {
0552     std::string BDTtreeName = "BDTtree_" + std::to_string(iEvent.id().event());
0553     TTree* BDTtree = fs_->make<TTree>(BDTtreeName.c_str(), "Tree for BDT tracks");
0554 
0555     BDTtree->Branch("Track_pt", &Ttrack_pt);
0556     BDTtree->Branch("Track_eta", &Ttrack_eta);
0557     BDTtree->Branch("Track_phi", &Ttrack_phi);
0558     BDTtree->Branch("Track_dz", &Ttrack_dz);
0559     BDTtree->Branch("Track_dxy", &Ttrack_dxy);
0560     BDTtree->Branch("Track_chi2", &Ttrack_chi2);
0561     BDTtree->Branch("Track_ndof", &Ttrack_ndof);
0562     BDTtree->Branch("Track_nValidHits", &Ttrack_nValidHits);
0563     BDTtree->Branch("Track_npixBarrelValidHits", &Ttrack_npixBarrelValidHits);
0564     BDTtree->Branch("Track_npixEndcapValidHits", &Ttrack_npixEndcapValidHits);
0565     BDTtree->Branch("Track_Signal", &Ttrack_Signal);
0566     BDTtree->Branch("Track_Associated", &Ttrack_Associated);
0567     BDTtree->Branch("Track_BTLchi2", &Ttrack_BTLchi2);
0568     BDTtree->Branch("Track_BTLtime_chi2", &Ttrack_BTLtime_chi2);
0569     BDTtree->Branch("Track_ETLchi2", &Ttrack_ETLchi2);
0570     BDTtree->Branch("Track_ETLtime_chi2", &Ttrack_ETLtime_chi2);
0571     BDTtree->Branch("Track_t0", &Ttrack_t0);
0572     BDTtree->Branch("Track_sigmat0", &Ttrack_sigmat0);
0573     BDTtree->Branch("Track_Tmtd", &Ttrack_Tmtd);
0574     BDTtree->Branch("Track_MtdMVA", &Ttrack_MtdMVA);
0575     BDTtree->Branch("Track_lHitPos", &Ttrack_lHitPos);
0576     BDTtree->Branch("Track_sigmaTmtd", &Ttrack_sigmaTmtd);
0577     BDTtree->Branch("Track_length", &Ttrack_length);
0578     BDTtree->Branch("Track_DirClu", &Ttrack_TPDirClu);
0579     BDTtree->Branch("Track_OtherClu", &Ttrack_TPOtherClu);
0580     BDTtree->Branch("Track_isBTL", &Ttrack_isBTL);
0581     BDTtree->Branch("Track_isETL", &Ttrack_isETL);
0582     BDTtree->Branch("Track_withMTD", &Ttrack_withMTD);
0583 
0584     Ttrack_pt.clear();
0585     Ttrack_eta.clear();
0586     Ttrack_phi.clear();
0587     Ttrack_dz.clear();
0588     Ttrack_dxy.clear();
0589     Ttrack_chi2.clear();
0590     Ttrack_ndof.clear();
0591     Ttrack_nValidHits.clear();
0592     Ttrack_npixBarrelValidHits.clear();
0593     Ttrack_npixEndcapValidHits.clear();
0594     Ttrack_Signal.clear();
0595     Ttrack_Associated.clear();
0596     Ttrack_BTLchi2.clear();
0597     Ttrack_BTLtime_chi2.clear();
0598     Ttrack_ETLchi2.clear();
0599     Ttrack_ETLtime_chi2.clear();
0600     Ttrack_t0.clear();
0601     Ttrack_sigmat0.clear();
0602     Ttrack_Tmtd.clear();
0603     Ttrack_MtdMVA.clear();
0604     Ttrack_lHitPos.clear();
0605     Ttrack_sigmaTmtd.clear();
0606     Ttrack_length.clear();
0607     Ttrack_TPDirClu.clear();
0608     Ttrack_TPOtherClu.clear();
0609     Ttrack_isBTL.clear();
0610     Ttrack_isETL.clear();
0611     Ttrack_withMTD.clear();
0612 
0613     unsigned int index = 0;
0614 
0615     // --- Loop over all RECO tracks ---
0616     for (const auto& trackGen : *RecTrackHandle) {
0617       const reco::TrackRef trackref(iEvent.getHandle(RecTrackToken_), index);
0618       index++;
0619 
0620       if (trackAssoc[trackref] == -1) {
0621         LogInfo("mtdTracks") << "Extended track not associated";
0622         continue;
0623       }
0624 
0625       const reco::TrackRef mtdTrackref = reco::TrackRef(iEvent.getHandle(RecMTDTrackToken_), trackAssoc[trackref]);
0626       const reco::Track& track = *mtdTrackref;
0627 
0628       bool withMTD = false;
0629       bool isBTL = false;
0630       bool isETL = false;
0631       bool twoETLdiscs = false;
0632       bool isTPDirClu = false;
0633       bool isTPOtherClu = false;
0634 
0635       bool isTPmtdDirectBTL = false, isTPmtdOtherBTL = false, isTPmtdDirectETL = false, isTPmtdOtherETL = false,
0636            isTPmtdCorrect = false;
0637 
0638       if (trkRecSel(trackGen)) {
0639         if (std::round(sigmatimemtd[trackref] - sigmat0Pid[trackref]) != 0) {
0640           LogWarning("mtdTracks")
0641               << "TimeError associated to refitted track is different from TimeError stored in tofPID "
0642                  "sigmat0 ValueMap: this should not happen";
0643         }
0644 
0645         std::vector<edm::Ref<edmNew::DetSetVector<FTLCluster>, FTLCluster>> recoClustersRefs;
0646 
0647         if (std::abs(trackGen.eta()) < etacutBTL_) {
0648           // --- all BTL tracks (with and without hit in MTD) ---
0649 
0650           for (const auto hit : track.recHits()) {
0651             if (hit->isValid() == false)
0652               continue;
0653             MTDDetId Hit = hit->geographicalId();
0654             if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 1)) {
0655               isBTL = true;
0656               const auto* mtdhit = static_cast<const MTDTrackingRecHit*>(hit);
0657               const auto& hitCluster = mtdhit->mtdCluster();
0658               if (hitCluster.size() != 0) {
0659                 auto recoClusterRef = edmNew::makeRefTo(btlRecCluHandle, &hitCluster);
0660                 recoClustersRefs.push_back(recoClusterRef);
0661               }
0662             }
0663           }
0664         }  //loop over (geometrical) BTL tracks
0665         else {
0666           // --- all ETL tracks (with and without hit in MTD) ---
0667           for (const auto hit : track.recHits()) {
0668             if (hit->isValid() == false)
0669               continue;
0670             MTDDetId Hit = hit->geographicalId();
0671             if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 2)) {
0672               isETL = true;
0673 
0674               const auto* mtdhit = static_cast<const MTDTrackingRecHit*>(hit);
0675               const auto& hitCluster = mtdhit->mtdCluster();
0676               if (hitCluster.size() != 0) {
0677                 auto recoClusterRef = edmNew::makeRefTo(etlRecCluHandle, &hitCluster);
0678                 recoClustersRefs.push_back(recoClusterRef);
0679               }
0680             }
0681           }
0682         }
0683 
0684         LogDebug("MVATrainingNtuple") << "Track p/pt = " << trackGen.p() << " " << trackGen.pt() << " eta "
0685                                       << trackGen.eta() << " BTL " << isBTL << " ETL " << isETL << " 2disks "
0686                                       << twoETLdiscs;
0687 
0688         // == TrackingParticle based matching
0689         const reco::TrackBaseRef trkrefb(trackref);
0690         auto tp_info = getAnyMatchedTP(trkrefb);
0691         if (tp_info != nullptr && trkTPSelAll(**tp_info)) {
0692           // ==  MC truth matching
0693 
0694           auto simClustersRefsIt = tp2SimAssociationMap.find(*tp_info);
0695           withMTD = (simClustersRefsIt != tp2SimAssociationMap.end());
0696 
0697           // If there is a mtdSimLayerCluster from the tracking particle
0698           if (withMTD) {
0699             // -- Get the refs to MtdSimLayerClusters associated to the TP
0700             std::vector<edm::Ref<MtdSimLayerClusterCollection>> simClustersRefs;
0701             for (const auto& ref : simClustersRefsIt->val) {
0702               simClustersRefs.push_back(ref);
0703             }
0704             // === ETL
0705             // -- Sort ETL sim clusters by time
0706             if (std::abs(trackGen.eta()) >= etacutBTL_ && !simClustersRefs.empty()) {
0707               std::sort(simClustersRefs.begin(), simClustersRefs.end(), [](const auto& a, const auto& b) {
0708                 return a->simLCTime() < b->simLCTime();
0709               });
0710               // Check if TP has direct or other sim cluster for BTL
0711               for (const auto& simClusterRef : simClustersRefs) {
0712                 if (simClusterRef->trackIdOffset() == 0) {
0713                   isTPmtdDirectETL = true;
0714                 } else if (simClusterRef->trackIdOffset() != 0) {
0715                   isTPmtdOtherETL = true;
0716                 }
0717               }
0718             }
0719             // === BTL
0720             // -- Sort BTL sim clusters by time
0721             if (std::abs(trackGen.eta()) < etacutBTL_ && !simClustersRefs.empty()) {
0722               std::sort(simClustersRefs.begin(), simClustersRefs.end(), [](const auto& a, const auto& b) {
0723                 return a->simLCTime() < b->simLCTime();
0724               });
0725               // Check if TP has direct or other sim cluster for BTL
0726               for (const auto& simClusterRef : simClustersRefs) {
0727                 if (simClusterRef->trackIdOffset() == 0) {
0728                   isTPmtdDirectBTL = true;
0729                 } else if (simClusterRef->trackIdOffset() != 0) {
0730                   isTPmtdOtherBTL = true;
0731                 }
0732               }
0733             }
0734 
0735             // ==  Check if the track-cluster association is correct: Track->RecoClus->SimClus == Track->TP->SimClus
0736             for (const auto& recClusterRef : recoClustersRefs) {
0737               if (recClusterRef.isNonnull()) {
0738                 auto itp = r2sAssociationMap.equal_range(recClusterRef);
0739                 if (itp.first != itp.second) {
0740                   auto& simClustersRefs_RecoMatch = (*itp.first).second;
0741                   for (const auto& simClusterRef_RecoMatch : simClustersRefs_RecoMatch) {
0742                     // Check if simClusterRef_RecoMatch  exists in SimClusters
0743                     auto simClusterIt =
0744                         std::find(simClustersRefs.begin(), simClustersRefs.end(), simClusterRef_RecoMatch);
0745                     // SimCluster found in SimClusters
0746                     if (simClusterIt != simClustersRefs.end()) {
0747                       isTPmtdCorrect = true;
0748                     }
0749                   }
0750                 }
0751               }
0752             }  /// end loop over reco clusters associated to this track.
0753           }  // --- end "withMTD"
0754         }  // TP matching
0755         Ttrack_pt.push_back(trackGen.pt());
0756         Ttrack_phi.push_back(trackGen.phi());
0757         Ttrack_eta.push_back(trackGen.eta());
0758         Ttrack_dz.push_back(std::abs(trackGen.dz(beamSpot.position())));
0759         Ttrack_dxy.push_back(std::abs(trackGen.dxy(beamSpot.position())));
0760         Ttrack_chi2.push_back(trackGen.chi2());
0761         Ttrack_ndof.push_back(trackGen.ndof());
0762         Ttrack_nValidHits.push_back(trackGen.numberOfValidHits());
0763 
0764         Ttrack_npixBarrelValidHits.push_back(npixBarrel[trackref]);
0765         Ttrack_npixEndcapValidHits.push_back(npixEndcap[trackref]);
0766         Ttrack_BTLchi2.push_back(btlMatchChi2[trackref]);
0767         Ttrack_BTLtime_chi2.push_back(btlMatchTimeChi2[trackref]);
0768         Ttrack_ETLchi2.push_back(etlMatchChi2[trackref]);
0769         Ttrack_ETLtime_chi2.push_back(etlMatchTimeChi2[trackref]);
0770 
0771         Ttrack_t0.push_back(t0Src[trackref]);
0772         Ttrack_sigmat0.push_back(sigmat0Src[trackref]);
0773         Ttrack_Tmtd.push_back(tMtd[trackref]);
0774         Ttrack_sigmaTmtd.push_back(sigmatimemtd[trackref]);
0775         Ttrack_length.push_back(pathLength[trackref]);
0776         Ttrack_MtdMVA.push_back(mtdQualMVA[trackref]);
0777         Ttrack_lHitPos.push_back(outermostHitPosition[trackref]);
0778 
0779         Ttrack_isBTL.push_back(isBTL);
0780         Ttrack_isETL.push_back(isETL);
0781         Ttrack_withMTD.push_back(withMTD);
0782 
0783         if ((tp_info != nullptr) && (*tp_info)->eventId().bunchCrossing() == 0 &&
0784             (*tp_info)->eventId().event() == 0) {  // Signal vs PU seperation
0785           Ttrack_Signal.push_back(true);           // Signal track
0786         } else {
0787           Ttrack_Signal.push_back(false);  // PU track?
0788         }
0789 
0790         if (isTPmtdCorrect) {
0791           Ttrack_Associated.push_back(true);  // Associated track
0792         } else {
0793           Ttrack_Associated.push_back(false);  // Not associated track
0794         }
0795         if (isTPmtdDirectBTL || isTPmtdDirectETL) {
0796           isTPDirClu = true;  // Direct cluster
0797         }
0798         if (isTPmtdOtherBTL || isTPmtdOtherETL) {
0799           isTPOtherClu = true;  // Not direct cluster
0800         }
0801 
0802         Ttrack_TPDirClu.push_back(isTPDirClu);
0803         Ttrack_TPOtherClu.push_back(isTPOtherClu);
0804       }  // trkRecSel
0805 
0806     }  // RECO tracks loop
0807 
0808     BDTtree->Fill();
0809 
0810   }  // ntuple for BDT
0811 }
0812 
0813 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0814 void MVATrainingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0815   edm::ParameterSetDescription desc;
0816   desc.add<edm::InputTag>("inputTracks", edm::InputTag("generalTracks"));
0817   desc.add<edm::InputTag>("inputTagT", edm::InputTag("trackExtenderWithMTD"));
0818   desc.add<edm::InputTag>("inputTagV", edm::InputTag("offlinePrimaryVertices4D"));
0819   desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
0820   desc.add<edm::InputTag>("tp2SimAssociationMapTag", edm::InputTag("mtdSimLayerClusterToTPAssociation"));
0821   desc.add<edm::InputTag>("r2sAssociationMapTag", edm::InputTag("mtdRecoClusterToSimLayerClusterAssociation"));
0822   desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
0823   desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
0824   desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
0825       ->setComment("Association between General and MTD Extended tracks");
0826   desc.add<edm::InputTag>("recCluTagBTL", edm::InputTag("mtdClusters", "FTLBarrel"));
0827   desc.add<edm::InputTag>("recCluTagETL", edm::InputTag("mtdClusters", "FTLEndcap"));
0828   desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
0829   desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
0830   desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
0831   desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
0832   desc.add<edm::InputTag>("t0Src", edm::InputTag("trackExtenderWithMTD:generalTrackt0"));
0833   desc.add<edm::InputTag>("sigmat0Src", edm::InputTag("trackExtenderWithMTD:generalTracksigmat0"));
0834   desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
0835   desc.add<edm::InputTag>("sigmat0PID", edm::InputTag("tofPID:sigmat0"));
0836   desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
0837   desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
0838   desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
0839   desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
0840   desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
0841   desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
0842   desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
0843   desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
0844   desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
0845   desc.add<edm::InputTag>("sigmatofpiSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofPi"));
0846   desc.add<edm::InputTag>("sigmatofkSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofK"));
0847   desc.add<edm::InputTag>("sigmatofpSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofP"));
0848   desc.add<edm::InputTag>("btlMatchChi2Src", edm::InputTag("trackExtenderWithMTD", "btlMatchChi2"));
0849   desc.add<edm::InputTag>("btlMatchTimeChi2Src", edm::InputTag("trackExtenderWithMTD", "btlMatchTimeChi2"));
0850   desc.add<edm::InputTag>("etlMatchChi2Src", edm::InputTag("trackExtenderWithMTD", "etlMatchChi2"));
0851   desc.add<edm::InputTag>("etlMatchTimeChi2Src", edm::InputTag("trackExtenderWithMTD", "etlMatchTimeChi2"));
0852   desc.add<edm::InputTag>("npixBarrelSrc", edm::InputTag("trackExtenderWithMTD", "npixBarrel"));
0853   desc.add<edm::InputTag>("npixEndcapSrc", edm::InputTag("trackExtenderWithMTD", "npixEndcap"));
0854   desc.add<edm::InputTag>("outermostHitPositionSrc",
0855                           edm::InputTag("trackExtenderWithMTD", "generalTrackOutermostHitPosition"));
0856   desc.addUntracked<std::string>("fileName", "file.root");
0857   desc.add<bool>("ntupleforBDT", true);
0858   desc.add<bool>("ntupleforGNN", false);
0859   {
0860     edm::ParameterSetDescription psd0;
0861     HITrackFilterForPVFinding::fillPSetDescription(psd0);  // extension of TrackFilterForPVFinding
0862     desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
0863   }
0864 
0865   descriptions.add("mvaTrainingNtuple", desc);
0866 }
0867 
0868 const bool MVATrainingNtuple::trkTPSelAll(const TrackingParticle& tp) {
0869   bool match = false;
0870 
0871   auto x_pv = tp.parentVertex()->position().x();
0872   auto y_pv = tp.parentVertex()->position().y();
0873   auto z_pv = tp.parentVertex()->position().z();
0874 
0875   auto r_pv = std::sqrt(x_pv * x_pv + y_pv * y_pv);
0876 
0877   match = tp.charge() != 0 && std::abs(tp.eta()) < etacutGEN_ &&
0878           ((std::abs(tp.eta()) < etacutBTL_ && tp.pt() > pTcutBTL_) ||
0879            (std::abs(tp.eta()) > etacutBTL_ && tp.pt() > pTcutETL_)) &&
0880           r_pv < rBTL_ && std::abs(z_pv) < zETL_;
0881   return match;
0882 }
0883 
0884 const bool MVATrainingNtuple::trkRecSel(const reco::TrackBase& trk) {
0885   bool match = false;
0886   match = std::abs(trk.eta()) <= etacutREC_ && ((std::abs(trk.eta()) <= etacutBTL_ && trk.pt() > pTcutBTL_) ||
0887                                                 (std::abs(trk.eta()) > etacutBTL_ && trk.pt() > pTcutETL_));
0888   return match;
0889 }
0890 
0891 //define this as a plug-in
0892 DEFINE_FWK_MODULE(MVATrainingNtuple);
[ 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!