From ce5eba3408f55a9843901b0ef6a968a392eceaa2 Mon Sep 17 00:00:00 2001 From: Lubynets Date: Fri, 24 Jul 2020 21:37:45 +0200 Subject: [PATCH 1/4] try to match with mc --- AnalysisTreeInterface/ConverterOut.cpp | 99 +++++++-- AnalysisTreeInterface/ConverterOut.h | 36 +++- AnalysisTreeInterface/KFManLinkDef.h | 1 + AnalysisTreeInterface/PFTaskManager.cpp | 3 +- AnalysisTreeInterface/RecoSimMatch.cpp | 259 ++++++++++++++++++++++++ AnalysisTreeInterface/RecoSimMatch.h | 124 ++++++++++++ AnalysisTreeInterface/main.cxx | 4 +- 7 files changed, 507 insertions(+), 19 deletions(-) create mode 100644 AnalysisTreeInterface/RecoSimMatch.cpp create mode 100644 AnalysisTreeInterface/RecoSimMatch.h diff --git a/AnalysisTreeInterface/ConverterOut.cpp b/AnalysisTreeInterface/ConverterOut.cpp index 859cc81..de60d74 100644 --- a/AnalysisTreeInterface/ConverterOut.cpp +++ b/AnalysisTreeInterface/ConverterOut.cpp @@ -16,18 +16,13 @@ void ConverterOut::Exec() const KFParticle& particle = candidate.GetParticle(); float mass, masserr; particle.GetMass(mass, masserr); -// lambdarec->SetField(mass, mass_field_id_); lambdarec->SetField(particle.DaughterIds()[0], daughter1_id_field_id_); lambdarec->SetField(particle.DaughterIds()[1], daughter2_id_field_id_); -// std::cout << "FC\t" << particle.DaughterIds()[0] << "\t" << particle.DaughterIds()[1] << "\n"; lambdarec->SetField(particle.X(), x_field_id_); lambdarec->SetField(particle.Y(), y_field_id_); lambdarec->SetField(particle.Z(), z_field_id_); lambdarec->SetMomentum(particle.GetPx(), particle.GetPy(), particle.GetPz()); -// lambdarec->SetField( particle.GetRapidity(), rap_lab_field_id_); - lambdarec->SetField( particle.GetRapidity() - data_header_->GetBeamRapidity(), rap_cm_field_id_); -// lambdarec->SetField( particle.GetPDG(), pdg_field_id_w_ ); lambdarec->SetPid( particle.GetPDG()); lambdarec->SetMass(mass); @@ -43,6 +38,44 @@ void ConverterOut::Exec() lambdarec->SetField( candidate.GetIsFromPV(), isfrompv_field_id_); lambdarec->SetField( candidate.GetCosineTopo(), cosinetopo_field_id_); lambdarec->SetField( candidate.GetChi2Topo(), chi2topo_field_id_); + + if(is_match_with_mc_){ + const int simtrackid1 = kf2sim_tracks_->GetMatch(particle.DaughterIds()[0]); + const int simtrackid2 = kf2sim_tracks_->GetMatch(particle.DaughterIds()[1]); + int is_signal = 0; + int mother_id = -999; + int mother_pdg = -1; + if(!(simtrackid1 < 0 || simtrackid2 < 0)) + { + const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); + const AnalysisTree::Particle& simtrack2 = sim_tracks_->GetChannel(simtrackid2); + if(simtrack1.GetField(mother_id_field_id_) == simtrack2.GetField(mother_id_field_id_)) + { + mother_id = simtrack1.GetField(mother_id_field_id_); + if(mother_id < 0) continue; + const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); + mother_pdg = simtrackmother.GetPid(); + if(mother_pdg == pdg_lambda) + is_signal = 1; + } + } + lambdarec->SetField( is_signal, is_signal_field_id_w_); + + if(is_signal!=1) continue; + + auto* lambdasim = lambda_sim_->AddChannel(); + lambdasim->Init(out_config_->GetBranchConfig(lambda_sim_->GetId())); + const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); + const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); + + lambdasim -> SetMass(mass_lambda); + lambdasim -> SetPid(mother_pdg); + lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); + lambdasim -> SetField(mother_id, sim_id_w_); + + lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); + } + } out_tree_->Fill(); } @@ -70,12 +103,11 @@ void ConverterOut::Init(std::map& branches) LambdaRecoBranch.AddField("x"); LambdaRecoBranch.AddField("y"); LambdaRecoBranch.AddField("z"); -// LambdaRecoBranch.AddField("invmass"); -// LambdaRecoBranch.AddField("rapidity_lab"); LambdaRecoBranch.AddField("rapidity_cm"); -// LambdaRecoBranch.AddField("pdg"); LambdaRecoBranch.AddField("daughter1id"); LambdaRecoBranch.AddField("daughter2id"); + if(is_match_with_mc_) + LambdaRecoBranch.AddField("issignal"); out_config_->AddBranchConfig( LambdaRecoBranch ); lambda_reco_ = new AnalysisTree::Particles(out_config_->GetLastId()); @@ -84,12 +116,10 @@ void ConverterOut::Init(std::map& branches) x_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("x"); y_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("y"); z_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("z"); -// mass_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("invmass"); -// rap_lab_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_lab"); - rap_cm_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_cm"); -// pdg_field_id_w_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("pdg"); daughter1_id_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter1id"); daughter2_id_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter2id"); + if(is_match_with_mc_) + is_signal_field_id_w_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("issignal"); chi2primpos_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primpos"); chi2primneg_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primneg"); @@ -102,4 +132,49 @@ void ConverterOut::Init(std::map& branches) isfrompv_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("isfrompv"); cosinetopo_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinetopo"); chi2topo_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2topo"); + + if(is_match_with_mc_){ + AnalysisTree::BranchConfig LambdaSimBranch("LambdaSimulated", AnalysisTree::DetType::kParticle); + LambdaSimBranch.AddField("x"); + LambdaSimBranch.AddField("y"); + LambdaSimBranch.AddField("z"); + LambdaSimBranch.AddField("simid"); + LambdaSimBranch.AddField("daughter1id"); + LambdaSimBranch.AddField("daughter2id"); + + out_config_->AddBranchConfig( LambdaSimBranch ); + lambda_sim_ = new AnalysisTree::Particles(out_config_->GetLastId()); + lambda_reco2sim_ = new AnalysisTree::Matching(out_config_->GetBranchConfig("LambdaCandidates").GetId(), out_config_->GetBranchConfig("LambdaSimulated").GetId()); + out_config_->AddMatch("LambdaCandidates", "LambdaSimulated", "LambdaCand2LambdaSim"); + out_tree_ -> Branch("LambdaSimulated", "AnalysisTree::Particles", &lambda_sim_); + out_tree_ -> Branch("LambdaCand2LambdaSim", "AnalysisTree::Matching", &lambda_reco2sim_); + + sim_id_w_ = out_config_->GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("simid"); + } +} + +void ConverterOut::InitSimAT(const std::string sim_file_name, const std::string sim_tree_name) +// void ConverterOut::InitSimAT(std::map& branches) +{ + in_file_sim_ = TFile::Open(sim_file_name.c_str(), "read"); + config_sim_ = (AnalysisTree::Configuration*) in_file_sim_->Get("Configuration"); + in_chain_sim_ = new TChain(sim_tree_name.c_str()); + in_chain_sim_ -> Add(sim_file_name.c_str()); +// in_chain_sim_ -> SetBranchAddress("KfpfTracks", &kf_tracks_); +// in_chain_sim_ -> SetBranchAddress("SimTracks", &sim_tracks_); +// in_chain_sim_ -> SetBranchAddress(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str(), &kf2sim_tracks_); + kf_tracks_ = (AnalysisTree::TrackDetector*)in_chain_sim_->GetBranch("KfpfTracks"); + sim_tracks_ = (AnalysisTree::Particles*)in_chain_sim_->GetBranch("SimTracks"); + kf2sim_tracks_ = (AnalysisTree::Matching*)in_chain_sim_->GetBranch(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str()); + +// kf_tracks_ = (AnalysisTree::TrackDetector*) branches.find("KfpfTracks")->second; +// sim_tracks_ = (AnalysisTree::Particles*) branches.find("SimTracks")->second; +// kf2sim_tracks_ = (AnalysisTree::Matching*) branches.find( +// config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str())->second; + + auto branch_conf_kftr = config_sim_->GetBranchConfig( "KfpfTracks" ); + par_field_id_ = branch_conf_kftr.GetFieldId("x"); + + auto branch_conf_simtr = config_sim_->GetBranchConfig( "SimTracks" ); + mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); } diff --git a/AnalysisTreeInterface/ConverterOut.h b/AnalysisTreeInterface/ConverterOut.h index f5355d3..e45a766 100644 --- a/AnalysisTreeInterface/ConverterOut.h +++ b/AnalysisTreeInterface/ConverterOut.h @@ -1,6 +1,10 @@ #ifndef KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ #define KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" + #include #include "AnalysisTree/FillTask.hpp" @@ -12,14 +16,31 @@ class ConverterOut : public AnalysisTree::FillTask { ~ConverterOut() override = default; void Init(std::map& branches) override; + void InitSimAT(const std::string sim_file_name, const std::string sim_tree_name); void Exec() override; void Finish() override {} void SetCandidates(const std::vector& canditates) { canditates_ = canditates; } + void SetIsMatchWithMc(bool ismatch = true){ is_match_with_mc_ = ismatch; }; protected: - + + // input info for matching with mc + TFile* in_file_sim_{nullptr}; + TChain* in_chain_sim_{nullptr}; + AnalysisTree::Configuration* config_sim_{nullptr}; + AnalysisTree::TrackDetector* kf_tracks_{nullptr}; + AnalysisTree::Particles* sim_tracks_{nullptr}; + AnalysisTree::Matching* kf2sim_tracks_{nullptr}; + + // field ids for input parameters of simulated lambda + int par_field_id_ {-1}; + int mother_id_field_id_{-1}; + + // output info AnalysisTree::Particles* lambda_reco_{nullptr}; + AnalysisTree::Particles* lambda_sim_{nullptr}; + AnalysisTree::Matching* lambda_reco2sim_{nullptr}; std::vector canditates_; // field ids for selected lambda candidates kinematic parameters @@ -27,9 +48,7 @@ class ConverterOut : public AnalysisTree::FillTask { int y_field_id_{-1}; int z_field_id_{-1}; int mass_field_id_{-1}; -// int rap_lab_field_id_{-1}; - int rap_cm_field_id_{-1}; -// int pdg_field_id_w_{-1}; + int is_signal_field_id_w_{-1}; int daughter1_id_field_id_{-1}; int daughter2_id_field_id_{-1}; @@ -45,6 +64,15 @@ class ConverterOut : public AnalysisTree::FillTask { int isfrompv_field_id_{-1}; int cosinetopo_field_id_{-1}; int chi2topo_field_id_{-1}; + + // field ids for simulated lamdas kinematic parameters + int sim_id_w_{-1}; + + bool is_match_with_mc_{true}; + + + const int pdg_lambda = 3122; + const float mass_lambda = 1.115683; }; #endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ diff --git a/AnalysisTreeInterface/KFManLinkDef.h b/AnalysisTreeInterface/KFManLinkDef.h index 306dd62..a83c0bc 100644 --- a/AnalysisTreeInterface/KFManLinkDef.h +++ b/AnalysisTreeInterface/KFManLinkDef.h @@ -7,5 +7,6 @@ #pragma link C++ class PFTaskManager+; #pragma link C++ class ConverterIn+; #pragma link C++ class ConverterOut+; +#pragma link C++ class RecoSimMatch+; #endif diff --git a/AnalysisTreeInterface/PFTaskManager.cpp b/AnalysisTreeInterface/PFTaskManager.cpp index aa0dce3..54b4699 100644 --- a/AnalysisTreeInterface/PFTaskManager.cpp +++ b/AnalysisTreeInterface/PFTaskManager.cpp @@ -19,9 +19,10 @@ void PFTaskManager::Run(long long int nEvents) { FCFinder.FindParticles(); auto* converter_out = ((ConverterOut*)tasks_.at(kOutConverter)); + converter_out->SetIsMatchWithMc(true); + converter_out->InitSimAT("1.analysistree.root", "aTree"); // TODO avoid hardcoding converter_out->SetCandidates(FCFinder.GetLambdas()); converter_out->Exec(); -// out_tree_->Fill(); } // Event loop } diff --git a/AnalysisTreeInterface/RecoSimMatch.cpp b/AnalysisTreeInterface/RecoSimMatch.cpp new file mode 100644 index 0000000..fd1fcc6 --- /dev/null +++ b/AnalysisTreeInterface/RecoSimMatch.cpp @@ -0,0 +1,259 @@ +#include "RecoSimMatch.h" + +void RecoSimMatch::InitSimAT(const std::string sim_file_name, const std::string sim_tree_name) +{ + in_file_sim_ = TFile::Open(sim_file_name.c_str(), "read"); + config_sim_ = (AnalysisTree::Configuration*) in_file_sim_->Get("Configuration"); + in_chain_sim_ = new TChain(sim_tree_name.c_str()); + in_chain_sim_ -> Add(sim_file_name.c_str()); + in_chain_sim_ -> SetBranchAddress("KfpfTracks", &kf_tracks_); + in_chain_sim_ -> SetBranchAddress("SimTracks", &sim_tracks_); + in_chain_sim_ -> SetBranchAddress("RecEventHeader", &rec_event_header_); + in_chain_sim_ -> SetBranchAddress("SimEventHeader", &sim_event_header_); + in_chain_sim_ -> SetBranchAddress(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str(), &kf2sim_tracks_); + + auto branch_conf_kftr = config_sim_->GetBranchConfig( "KfpfTracks" ); + par_field_id_ = branch_conf_kftr.GetFieldId("x"); + + auto branch_conf_simtr = config_sim_->GetBranchConfig( "SimTracks" ); + mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); +} + +void RecoSimMatch::InitRecoAT(const std::string reco_file_name, const std::string reco_tree_name) +{ + in_file_reco_ = TFile::Open(reco_file_name.c_str(), "read"); + in_chain_reco_ = new TChain(reco_tree_name.c_str()); + in_chain_reco_ -> Add(reco_file_name.c_str()); + in_chain_reco_ -> SetBranchAddress("LambdaCandidates", &reco_tracks_); + + config_reco_ = (AnalysisTree::Configuration*) in_file_reco_->Get("Configuration"); + + auto branch_conf_recotr = config_reco_->GetBranchConfig("LambdaCandidates"); + + mass_field_id_ = branch_conf_recotr.GetFieldId("invmass"); + rap_lab_field_id_ = branch_conf_recotr.GetFieldId("rapidity_lab"); + rap_cm_field_id_ = branch_conf_recotr.GetFieldId("rapidity_cm"); + pdg_field_id_ = branch_conf_recotr.GetFieldId("pdg"); + x_field_id_ = branch_conf_recotr.GetFieldId("x"); + y_field_id_ = branch_conf_recotr.GetFieldId("y"); + z_field_id_ = branch_conf_recotr.GetFieldId("z"); + + daughter1_id_field_id_ = branch_conf_recotr.GetFieldId("daughter1id"); + daughter2_id_field_id_ = branch_conf_recotr.GetFieldId("daughter2id"); + + chi2primpos_field_id_ = branch_conf_recotr.GetFieldId("chi2primpos"); + chi2primneg_field_id_ = branch_conf_recotr.GetFieldId("chi2primneg"); + distance_field_id_ = branch_conf_recotr.GetFieldId("distance"); + cosinepos_field_id_ = branch_conf_recotr.GetFieldId("cosinepos"); + cosineneg_field_id_ = branch_conf_recotr.GetFieldId("cosineneg"); + chi2geo_field_id_ = branch_conf_recotr.GetFieldId("chi2geo"); + l_field_id_ = branch_conf_recotr.GetFieldId("l"); + ldl_field_id_ = branch_conf_recotr.GetFieldId("ldl"); + isfrompv_field_id_ = branch_conf_recotr.GetFieldId("isfrompv"); + cosinetopo_field_id_ = branch_conf_recotr.GetFieldId("cosinetopo"); + chi2topo_field_id_ = branch_conf_recotr.GetFieldId("chi2topo"); +} + +void RecoSimMatch::InitOutput(const std::string out_file_name) +{ + out_file_ = TFile::Open(out_file_name.c_str(), "recreate"); + +// ***** Lambda Candidates ******* + + AnalysisTree::BranchConfig LambdaRecoBranch("LambdaCandidates", AnalysisTree::DetType::kParticle); + + LambdaRecoBranch.AddField("chi2primpos"); + LambdaRecoBranch.AddField("chi2primneg"); + LambdaRecoBranch.AddField("distance"); + LambdaRecoBranch.AddField("cosinepos"); + LambdaRecoBranch.AddField("cosineneg"); + LambdaRecoBranch.AddField("chi2geo"); + LambdaRecoBranch.AddField("l"); + LambdaRecoBranch.AddField("ldl"); + LambdaRecoBranch.AddField ("isfrompv"); + LambdaRecoBranch.AddField("cosinetopo"); + LambdaRecoBranch.AddField("chi2topo"); + LambdaRecoBranch.AddField ("iszok"); + + LambdaRecoBranch.AddField("x"); + LambdaRecoBranch.AddField("y"); + LambdaRecoBranch.AddField("z"); + LambdaRecoBranch.AddField("invmass"); + LambdaRecoBranch.AddField("rapidity_lab"); + LambdaRecoBranch.AddField("rapidity_cm"); + LambdaRecoBranch.AddField("pdg"); + LambdaRecoBranch.AddField("daughter1id"); + LambdaRecoBranch.AddField("daughter2id"); + LambdaRecoBranch.AddField("issignal"); + + config_out_.AddBranchConfig( LambdaRecoBranch ); + lambda_reco_ = new AnalysisTree::Particles(config_out_.GetLastId()); + +// ***** Lambda simulated ******* + + AnalysisTree::BranchConfig LambdaSimBranch("LambdaSimulated", AnalysisTree::DetType::kParticle); + + LambdaSimBranch.AddField("x"); + LambdaSimBranch.AddField("y"); + LambdaSimBranch.AddField("z"); + LambdaSimBranch.AddField("invmass"); + LambdaSimBranch.AddField("rapidity_lab"); + LambdaSimBranch.AddField("rapidity_cm"); + LambdaSimBranch.AddField("pdg"); + LambdaSimBranch.AddField("simid"); + LambdaSimBranch.AddField("daughter1id"); + LambdaSimBranch.AddField("daughter2id"); + + config_out_.AddBranchConfig( LambdaSimBranch ); + lambda_sim_ = new AnalysisTree::Particles(config_out_.GetLastId()); + +// *********** Matching ********** + lambda_reco2sim_ = new AnalysisTree::Matching(config_out_.GetBranchConfig("LambdaCandidates").GetId(), config_out_.GetBranchConfig("LambdaSimulated").GetId()); + config_out_.AddMatch("LambdaCandidates", "LambdaSimulated", "LambdaCand2LambdaSim"); + + + out_tree_ = new TTree("aTree", "AnalysisTree Lambda"); + out_tree_ -> Branch("LambdaCandidates", "AnalysisTree::Particles", &lambda_reco_); + out_tree_ -> Branch("LambdaSimulated", "AnalysisTree::Particles", &lambda_sim_); + out_tree_ -> Branch("LambdaCand2LambdaSim", "AnalysisTree::Matching", &lambda_reco2sim_); + + config_out_.Write("Configuration"); + + x_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("x"); + y_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("y"); + z_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("z"); + mass_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("invmass"); + rap_lab_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_lab"); + rap_cm_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_cm"); + pdg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("pdg"); + daughter1_id_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter1id"); + daughter2_id_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter2id"); + is_signal_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("issignal"); + + chi2primpos_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primpos"); + chi2primneg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primneg"); + distance_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("distance"); + cosinepos_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinepos"); + cosineneg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosineneg"); + chi2geo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2geo"); + l_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("l"); + ldl_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("ldl"); + isfrompv_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("isfrompv"); + cosinetopo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinetopo"); + chi2topo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2topo"); + iszok_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("iszok"); + + // ------------------------------- + mass_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("invmass"); + rap_lab_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("rapidity_lab"); + rap_cm_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("rapidity_cm"); + pdg_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("pdg"); + sim_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("simid"); +} + +void RecoSimMatch::Run(int nEvents) +{ + if(in_chain_reco_->GetEntries() != in_chain_sim_->GetEntries()) + std::cout << "Error! Files are not matchable!\n"; + if (nEvents == -1 || nEvents > in_chain_reco_->GetEntries()) + nEvents = in_chain_reco_->GetEntries(); + + for(int iEvent=0; iEventGetEntry(iEvent); + in_chain_sim_->GetEntry(iEvent); + + lambda_reco_->ClearChannels(); + lambda_sim_->ClearChannels(); + lambda_reco2sim_->Clear(); + + for(int iCh=0; iChGetNumberOfChannels(); iCh++) + { + const AnalysisTree::Particle& recotrack = reco_tracks_->GetChannel(iCh); + + auto* lambdarec = lambda_reco_->AddChannel(); + lambdarec->Init(config_out_.GetBranchConfig(lambda_reco_->GetId())); + lambdarec->SetField(recotrack.GetField(mass_field_id_), mass_field_id_w_); + lambdarec->SetField(recotrack.GetField(daughter1_id_field_id_), daughter1_id_field_id_w_); + lambdarec->SetField(recotrack.GetField(daughter2_id_field_id_), daughter2_id_field_id_w_); + + lambdarec->SetMomentum(recotrack.GetPx(), recotrack.GetPy(), recotrack.GetPz()); + lambdarec->SetField(recotrack.GetField(rap_lab_field_id_), rap_lab_field_id_w_); + lambdarec->SetField(recotrack.GetField(rap_cm_field_id_), rap_cm_field_id_w_); + lambdarec->SetField(recotrack.GetField(pdg_field_id_), pdg_field_id_w_ ); + + lambdarec->SetPid(recotrack.GetField(pdg_field_id_)); + lambdarec->SetMass(recotrack.GetMass()); + lambdarec->SetField(recotrack.GetField(x_field_id_), x_field_id_w_); + lambdarec->SetField(recotrack.GetField(y_field_id_), y_field_id_w_); + lambdarec->SetField(recotrack.GetField(z_field_id_), z_field_id_w_); + + lambdarec->SetField(recotrack.GetField(chi2primpos_field_id_), chi2primpos_field_id_w_); + lambdarec->SetField(recotrack.GetField(chi2primneg_field_id_), chi2primneg_field_id_w_); + lambdarec->SetField(recotrack.GetField(distance_field_id_), distance_field_id_w_); + lambdarec->SetField(recotrack.GetField(cosinepos_field_id_), cosinepos_field_id_w_); + lambdarec->SetField(recotrack.GetField(cosineneg_field_id_), cosineneg_field_id_w_); + lambdarec->SetField(recotrack.GetField(chi2geo_field_id_), chi2geo_field_id_w_); + lambdarec->SetField(recotrack.GetField(l_field_id_), l_field_id_w_); + lambdarec->SetField(recotrack.GetField(ldl_field_id_), ldl_field_id_w_); + lambdarec->SetField(recotrack.GetField(isfrompv_field_id_), isfrompv_field_id_w_); + lambdarec->SetField(recotrack.GetField(cosinetopo_field_id_), cosinetopo_field_id_w_); + lambdarec->SetField(recotrack.GetField(chi2topo_field_id_), chi2topo_field_id_w_); + + const AnalysisTree::Track& reco_daughter1 = kf_tracks_->GetChannel(recotrack.GetField(daughter1_id_field_id_)); + const AnalysisTree::Track& reco_daughter2 = kf_tracks_->GetChannel(recotrack.GetField(daughter2_id_field_id_)); + + float z_daughter_min = std::min(reco_daughter1.GetField(par_field_id_+2), reco_daughter2.GetField(par_field_id_+2)); + if(recotrack.GetField(z_field_id_) < z_daughter_min) + lambdarec->SetField(1, iszok_field_id_w_); + else + lambdarec->SetField(0, iszok_field_id_w_); + + const int simtrackid1 = kf2sim_tracks_->GetMatch(recotrack.GetField(daughter1_id_field_id_)); + const int simtrackid2 = kf2sim_tracks_->GetMatch(recotrack.GetField(daughter2_id_field_id_)); + int is_signal = 0; + int mother_id = -999; + int mother_pdg = -1; + if(!(simtrackid1 < 0 || simtrackid2 < 0)) + { + const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); + const AnalysisTree::Particle& simtrack2 = sim_tracks_->GetChannel(simtrackid2); + if(simtrack1.GetField(mother_id_field_id_) == simtrack2.GetField(mother_id_field_id_)) + { + mother_id = simtrack1.GetField(mother_id_field_id_); + if(mother_id < 0) continue; + const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); + mother_pdg = simtrackmother.GetPid(); + if(mother_pdg == pdg_lambda) + is_signal = 1; + } + } + lambdarec->SetField( is_signal, is_signal_field_id_w_); + + if(is_signal!=1) continue; + + auto* lambdasim = lambda_sim_->AddChannel(); + lambdasim->Init(config_out_.GetBranchConfig(lambda_sim_->GetId())); + const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); + const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); + + lambdasim -> SetField(mass_lambda, mass_sim_field_id_w_); + lambdasim -> SetField(mother_pdg, pdg_sim_field_id_w_); + lambdasim -> SetPid(mother_pdg); + lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); + lambdasim -> SetField(simtrackmother.GetRapidity(), rap_lab_sim_field_id_w_); + lambdasim -> SetField(simtrackmother.GetRapidity() - RapidityShift(12.), rap_cm_sim_field_id_w_); + lambdasim -> SetField(mother_id, sim_id_w_); + + lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); + } + out_tree_->Fill(); + } + out_tree_->Write(); +} + +float RecoSimMatch::RapidityShift(float pbeam) +{ + const float nucleonmass = 0.939; + return TMath::Log(((sqrt(nucleonmass*nucleonmass + pbeam*pbeam) + pbeam) / (sqrt(nucleonmass*nucleonmass + pbeam*pbeam) - pbeam)))/4.; +} diff --git a/AnalysisTreeInterface/RecoSimMatch.h b/AnalysisTreeInterface/RecoSimMatch.h new file mode 100644 index 0000000..5d80c6f --- /dev/null +++ b/AnalysisTreeInterface/RecoSimMatch.h @@ -0,0 +1,124 @@ +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" + +#include "AnalysisTree/Track.h" +#include "AnalysisTree/Hit.h" +#include "AnalysisTree/Detector.h" +#include "AnalysisTree/Configuration.h" +#include "AnalysisTree/EventHeader.h" +#include "AnalysisTree/Matching.h" + +class RecoSimMatch +{ +public: + + // Constructors/Destructors --------- + RecoSimMatch() = default; + virtual ~RecoSimMatch() = default; + + void InitSimAT(const std::string sim_file_name, const std::string sim_tree_name); + void InitRecoAT(const std::string reco_file_name, const std::string reco_tree_name); + void InitOutput(const std::string out_file_name); + void Run(int nEvents=-1); + +private: + + float RapidityShift(float pbeam); + + TFile* in_file_sim_{nullptr}; + TFile* in_file_reco_{nullptr}; + TChain* in_chain_sim_{nullptr}; + TChain* in_chain_reco_{nullptr}; + TFile* out_file_{nullptr}; + TFile* out_file_tree_{nullptr}; + + AnalysisTree::Configuration* config_sim_{nullptr}; + AnalysisTree::Configuration* config_reco_{nullptr}; + AnalysisTree::Configuration config_out_; + + AnalysisTree::Particles* lambda_reco_{nullptr}; + AnalysisTree::Particles* lambda_sim_{nullptr}; + AnalysisTree::Matching* lambda_reco2sim_{nullptr}; + + TTree* out_tree_{nullptr}; + + //sims + AnalysisTree::TrackDetector* kf_tracks_{nullptr}; + AnalysisTree::Particles* sim_tracks_{nullptr}; + AnalysisTree::Matching* kf2sim_tracks_{nullptr}; + AnalysisTree::EventHeader* rec_event_header_{nullptr}; + AnalysisTree::EventHeader* sim_event_header_{nullptr}; + + //recos + AnalysisTree::Particles* reco_tracks_{nullptr}; + + //-----------INPUT field ids---------------------------------- + // field ids for input parameters of simulated lambda + int par_field_id_ {-1}; + int mother_id_field_id_{-1}; + + // field ids for input kinematic parameters of selected lambda candidates + int mass_field_id_{-1}; + int rap_lab_field_id_{-1}; + int rap_cm_field_id_{-1}; + int pdg_field_id_{-1}; + int x_field_id_{-1}; + int y_field_id_{-1}; + int z_field_id_{-1}; + + int daughter1_id_field_id_{-1}; + int daughter2_id_field_id_{-1}; + + // field ids for input selection parameters of lambda candidates (cutting) + int chi2primpos_field_id_{-1}; + int chi2primneg_field_id_{-1}; + int distance_field_id_{-1}; + int cosinepos_field_id_{-1}; + int cosineneg_field_id_{-1}; + int chi2geo_field_id_{-1}; + int l_field_id_{-1}; + int ldl_field_id_{-1}; + int isfrompv_field_id_{-1}; + int cosinetopo_field_id_{-1}; + int chi2topo_field_id_{-1}; + //------------------------------------------------------------ + + //--------OUTPUT field ids------------------------------------ + // field ids for selected lambda candidates kinematic parameters + int mass_field_id_w_{-1}; + int rap_lab_field_id_w_{-1}; + int rap_cm_field_id_w_{-1}; + int pdg_field_id_w_{-1}; + int daughter1_id_field_id_w_{-1}; + int daughter2_id_field_id_w_{-1}; + int is_signal_field_id_w_{-1}; + int x_field_id_w_{-1}; + int y_field_id_w_{-1}; + int z_field_id_w_{-1}; + + // field ids for simulated lamdas kinematic parameters + int mass_sim_field_id_w_{-1}; + int rap_lab_sim_field_id_w_{-1}; + int rap_cm_sim_field_id_w_{-1}; + int pdg_sim_field_id_w_{-1}; + int sim_id_w_{-1}; + + // field ids for lambda candidate cutting parameters + int chi2primpos_field_id_w_{-1}; + int chi2primneg_field_id_w_{-1}; + int distance_field_id_w_{-1}; + int cosinepos_field_id_w_{-1}; + int cosineneg_field_id_w_{-1}; + int chi2geo_field_id_w_{-1}; + int l_field_id_w_{-1}; + int ldl_field_id_w_{-1}; + int isfrompv_field_id_w_{-1}; + int cosinetopo_field_id_w_{-1}; + int chi2topo_field_id_w_{-1}; + int iszok_field_id_w_{-1}; + + const int pdg_lambda = 3122; + const float mass_lambda = 1.115683; + +}; \ No newline at end of file diff --git a/AnalysisTreeInterface/main.cxx b/AnalysisTreeInterface/main.cxx index ff99cb4..8d690de 100644 --- a/AnalysisTreeInterface/main.cxx +++ b/AnalysisTreeInterface/main.cxx @@ -9,7 +9,7 @@ int main(int argc, char** argv) return EXIT_FAILURE; } - const bool make_plain_tree{true}; + const bool make_plain_tree{false}; CutsContainer cuts; cuts.CancelCuts(); @@ -31,7 +31,7 @@ int main(int argc, char** argv) man.AddTasks(in_converter, new ConverterOut); man.Init(); - man.Run(100); // -1 = all events + man.Run(-1); // -1 = all events man.Finish(); if(make_plain_tree){ -- GitLab From 3c39550f6d70a3558e87e4a98383de1d31ce5ff4 Mon Sep 17 00:00:00 2001 From: Lubynets Date: Sun, 26 Jul 2020 20:44:27 +0200 Subject: [PATCH 2/4] force replace by master --- AnalysisTreeInterface/ConverterOut.cpp | 149 +++++--------- AnalysisTreeInterface/ConverterOut.h | 43 +--- AnalysisTreeInterface/KFManLinkDef.h | 1 - AnalysisTreeInterface/PFTaskManager.cpp | 3 +- AnalysisTreeInterface/RecoSimMatch.cpp | 259 ------------------------ AnalysisTreeInterface/RecoSimMatch.h | 124 ------------ AnalysisTreeInterface/main.cxx | 11 +- 7 files changed, 66 insertions(+), 524 deletions(-) delete mode 100644 AnalysisTreeInterface/RecoSimMatch.cpp delete mode 100644 AnalysisTreeInterface/RecoSimMatch.h diff --git a/AnalysisTreeInterface/ConverterOut.cpp b/AnalysisTreeInterface/ConverterOut.cpp index de60d74..34c0806 100644 --- a/AnalysisTreeInterface/ConverterOut.cpp +++ b/AnalysisTreeInterface/ConverterOut.cpp @@ -16,13 +16,18 @@ void ConverterOut::Exec() const KFParticle& particle = candidate.GetParticle(); float mass, masserr; particle.GetMass(mass, masserr); +// lambdarec->SetField(mass, mass_field_id_); lambdarec->SetField(particle.DaughterIds()[0], daughter1_id_field_id_); lambdarec->SetField(particle.DaughterIds()[1], daughter2_id_field_id_); +// std::cout << "FC\t" << particle.DaughterIds()[0] << "\t" << particle.DaughterIds()[1] << "\n"; lambdarec->SetField(particle.X(), x_field_id_); lambdarec->SetField(particle.Y(), y_field_id_); lambdarec->SetField(particle.Z(), z_field_id_); lambdarec->SetMomentum(particle.GetPx(), particle.GetPy(), particle.GetPz()); +// lambdarec->SetField( particle.GetRapidity(), rap_lab_field_id_); + lambdarec->SetField( particle.GetRapidity() - data_header_->GetBeamRapidity(), rap_cm_field_id_); +// lambdarec->SetField( particle.GetPDG(), pdg_field_id_w_ ); lambdarec->SetPid( particle.GetPDG()); lambdarec->SetMass(mass); @@ -38,45 +43,8 @@ void ConverterOut::Exec() lambdarec->SetField( candidate.GetIsFromPV(), isfrompv_field_id_); lambdarec->SetField( candidate.GetCosineTopo(), cosinetopo_field_id_); lambdarec->SetField( candidate.GetChi2Topo(), chi2topo_field_id_); - - if(is_match_with_mc_){ - const int simtrackid1 = kf2sim_tracks_->GetMatch(particle.DaughterIds()[0]); - const int simtrackid2 = kf2sim_tracks_->GetMatch(particle.DaughterIds()[1]); - int is_signal = 0; - int mother_id = -999; - int mother_pdg = -1; - if(!(simtrackid1 < 0 || simtrackid2 < 0)) - { - const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); - const AnalysisTree::Particle& simtrack2 = sim_tracks_->GetChannel(simtrackid2); - if(simtrack1.GetField(mother_id_field_id_) == simtrack2.GetField(mother_id_field_id_)) - { - mother_id = simtrack1.GetField(mother_id_field_id_); - if(mother_id < 0) continue; - const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); - mother_pdg = simtrackmother.GetPid(); - if(mother_pdg == pdg_lambda) - is_signal = 1; - } - } - lambdarec->SetField( is_signal, is_signal_field_id_w_); - - if(is_signal!=1) continue; - - auto* lambdasim = lambda_sim_->AddChannel(); - lambdasim->Init(out_config_->GetBranchConfig(lambda_sim_->GetId())); - const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); - const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); - - lambdasim -> SetMass(mass_lambda); - lambdasim -> SetPid(mother_pdg); - lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); - lambdasim -> SetField(mother_id, sim_id_w_); - - lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); - } - } + MatchWithMc(); out_tree_->Fill(); } @@ -86,6 +54,11 @@ void ConverterOut::Init(std::map& branches) assert(out_config_ && out_tree_); out_branch_ = "LambdaCandidates"; + if(!in_branches_.empty()){ + mc_particles_ = (AnalysisTree::Particles*) branches.find(in_branches_[0])->second; + rec_tracks_ = (AnalysisTree::TrackDetector*) branches.find(in_branches_[1])->second; + } + AnalysisTree::BranchConfig LambdaRecoBranch(out_branch_, AnalysisTree::DetType::kParticle); LambdaRecoBranch.AddField("chi2primpos"); @@ -106,75 +79,47 @@ void ConverterOut::Init(std::map& branches) LambdaRecoBranch.AddField("rapidity_cm"); LambdaRecoBranch.AddField("daughter1id"); LambdaRecoBranch.AddField("daughter2id"); - if(is_match_with_mc_) - LambdaRecoBranch.AddField("issignal"); + if(mc_particles_) { + LambdaRecoBranch.AddField("is_signal"); + } out_config_->AddBranchConfig( LambdaRecoBranch ); lambda_reco_ = new AnalysisTree::Particles(out_config_->GetLastId()); out_tree_ -> Branch(out_branch_.c_str(), "AnalysisTree::Particles", &lambda_reco_); - x_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("x"); - y_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("y"); - z_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("z"); - daughter1_id_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter1id"); - daughter2_id_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter2id"); - if(is_match_with_mc_) - is_signal_field_id_w_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("issignal"); - - chi2primpos_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primpos"); - chi2primneg_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primneg"); - distance_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("distance"); - cosinepos_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinepos"); - cosineneg_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosineneg"); - chi2geo_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2geo"); - l_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("l"); - ldl_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("ldl"); - isfrompv_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("isfrompv"); - cosinetopo_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinetopo"); - chi2topo_field_id_ = out_config_->GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2topo"); - - if(is_match_with_mc_){ - AnalysisTree::BranchConfig LambdaSimBranch("LambdaSimulated", AnalysisTree::DetType::kParticle); - LambdaSimBranch.AddField("x"); - LambdaSimBranch.AddField("y"); - LambdaSimBranch.AddField("z"); - LambdaSimBranch.AddField("simid"); - LambdaSimBranch.AddField("daughter1id"); - LambdaSimBranch.AddField("daughter2id"); - - out_config_->AddBranchConfig( LambdaSimBranch ); - lambda_sim_ = new AnalysisTree::Particles(out_config_->GetLastId()); - lambda_reco2sim_ = new AnalysisTree::Matching(out_config_->GetBranchConfig("LambdaCandidates").GetId(), out_config_->GetBranchConfig("LambdaSimulated").GetId()); - out_config_->AddMatch("LambdaCandidates", "LambdaSimulated", "LambdaCand2LambdaSim"); - out_tree_ -> Branch("LambdaSimulated", "AnalysisTree::Particles", &lambda_sim_); - out_tree_ -> Branch("LambdaCand2LambdaSim", "AnalysisTree::Matching", &lambda_reco2sim_); - - sim_id_w_ = out_config_->GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("simid"); - } + InitIndexes(); } -void ConverterOut::InitSimAT(const std::string sim_file_name, const std::string sim_tree_name) -// void ConverterOut::InitSimAT(std::map& branches) -{ - in_file_sim_ = TFile::Open(sim_file_name.c_str(), "read"); - config_sim_ = (AnalysisTree::Configuration*) in_file_sim_->Get("Configuration"); - in_chain_sim_ = new TChain(sim_tree_name.c_str()); - in_chain_sim_ -> Add(sim_file_name.c_str()); -// in_chain_sim_ -> SetBranchAddress("KfpfTracks", &kf_tracks_); -// in_chain_sim_ -> SetBranchAddress("SimTracks", &sim_tracks_); -// in_chain_sim_ -> SetBranchAddress(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str(), &kf2sim_tracks_); - kf_tracks_ = (AnalysisTree::TrackDetector*)in_chain_sim_->GetBranch("KfpfTracks"); - sim_tracks_ = (AnalysisTree::Particles*)in_chain_sim_->GetBranch("SimTracks"); - kf2sim_tracks_ = (AnalysisTree::Matching*)in_chain_sim_->GetBranch(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str()); - -// kf_tracks_ = (AnalysisTree::TrackDetector*) branches.find("KfpfTracks")->second; -// sim_tracks_ = (AnalysisTree::Particles*) branches.find("SimTracks")->second; -// kf2sim_tracks_ = (AnalysisTree::Matching*) branches.find( -// config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str())->second; - - auto branch_conf_kftr = config_sim_->GetBranchConfig( "KfpfTracks" ); - par_field_id_ = branch_conf_kftr.GetFieldId("x"); - - auto branch_conf_simtr = config_sim_->GetBranchConfig( "SimTracks" ); - mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); +void ConverterOut::MatchWithMc(){ + + std::cout << mc_particles_->GetNumberOfChannels() << " " << rec_tracks_->GetNumberOfChannels() << std::endl; + } + + +void ConverterOut::InitIndexes(){ + + const auto& out_branch = out_config_->GetBranchConfig( lambda_reco_->GetId()); + + x_field_id_ = out_branch.GetFieldId("x"); + y_field_id_ = out_branch.GetFieldId("y"); + z_field_id_ = out_branch.GetFieldId("z"); +// mass_field_id_ = out_branch.GetFieldId("invmass"); +// rap_lab_field_id_ = out_branch.GetFieldId("rapidity_lab"); + rap_cm_field_id_ = out_branch.GetFieldId("rapidity_cm"); +// pdg_field_id_w_ = out_branch.GetFieldId("pdg"); + daughter1_id_field_id_ = out_branch.GetFieldId("daughter1id"); + daughter2_id_field_id_ = out_branch.GetFieldId("daughter2id"); + + chi2primpos_field_id_ = out_branch.GetFieldId("chi2primpos"); + chi2primneg_field_id_ = out_branch.GetFieldId("chi2primneg"); + distance_field_id_ = out_branch.GetFieldId("distance"); + cosinepos_field_id_ = out_branch.GetFieldId("cosinepos"); + cosineneg_field_id_ = out_branch.GetFieldId("cosineneg"); + chi2geo_field_id_ = out_branch.GetFieldId("chi2geo"); + l_field_id_ = out_branch.GetFieldId("l"); + ldl_field_id_ = out_branch.GetFieldId("ldl"); + isfrompv_field_id_ = out_branch.GetFieldId("isfrompv"); + cosinetopo_field_id_ = out_branch.GetFieldId("cosinetopo"); + chi2topo_field_id_ = out_branch.GetFieldId("chi2topo"); +} \ No newline at end of file diff --git a/AnalysisTreeInterface/ConverterOut.h b/AnalysisTreeInterface/ConverterOut.h index e45a766..f24b418 100644 --- a/AnalysisTreeInterface/ConverterOut.h +++ b/AnalysisTreeInterface/ConverterOut.h @@ -1,10 +1,6 @@ #ifndef KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ #define KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ -#include "TFile.h" -#include "TTree.h" -#include "TChain.h" - #include #include "AnalysisTree/FillTask.hpp" @@ -16,31 +12,21 @@ class ConverterOut : public AnalysisTree::FillTask { ~ConverterOut() override = default; void Init(std::map& branches) override; - void InitSimAT(const std::string sim_file_name, const std::string sim_tree_name); void Exec() override; void Finish() override {} void SetCandidates(const std::vector& canditates) { canditates_ = canditates; } - void SetIsMatchWithMc(bool ismatch = true){ is_match_with_mc_ = ismatch; }; protected: - - // input info for matching with mc - TFile* in_file_sim_{nullptr}; - TChain* in_chain_sim_{nullptr}; - AnalysisTree::Configuration* config_sim_{nullptr}; - AnalysisTree::TrackDetector* kf_tracks_{nullptr}; - AnalysisTree::Particles* sim_tracks_{nullptr}; - AnalysisTree::Matching* kf2sim_tracks_{nullptr}; - - // field ids for input parameters of simulated lambda - int par_field_id_ {-1}; - int mother_id_field_id_{-1}; - - // output info + + void InitIndexes(); + void MatchWithMc(); + AnalysisTree::Particles* lambda_reco_{nullptr}; - AnalysisTree::Particles* lambda_sim_{nullptr}; - AnalysisTree::Matching* lambda_reco2sim_{nullptr}; + + AnalysisTree::Particles* mc_particles_{nullptr}; + AnalysisTree::TrackDetector* rec_tracks_{nullptr}; + std::vector canditates_; // field ids for selected lambda candidates kinematic parameters @@ -48,7 +34,9 @@ class ConverterOut : public AnalysisTree::FillTask { int y_field_id_{-1}; int z_field_id_{-1}; int mass_field_id_{-1}; - int is_signal_field_id_w_{-1}; +// int rap_lab_field_id_{-1}; + int rap_cm_field_id_{-1}; +// int pdg_field_id_w_{-1}; int daughter1_id_field_id_{-1}; int daughter2_id_field_id_{-1}; @@ -64,15 +52,6 @@ class ConverterOut : public AnalysisTree::FillTask { int isfrompv_field_id_{-1}; int cosinetopo_field_id_{-1}; int chi2topo_field_id_{-1}; - - // field ids for simulated lamdas kinematic parameters - int sim_id_w_{-1}; - - bool is_match_with_mc_{true}; - - - const int pdg_lambda = 3122; - const float mass_lambda = 1.115683; }; #endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ diff --git a/AnalysisTreeInterface/KFManLinkDef.h b/AnalysisTreeInterface/KFManLinkDef.h index a83c0bc..306dd62 100644 --- a/AnalysisTreeInterface/KFManLinkDef.h +++ b/AnalysisTreeInterface/KFManLinkDef.h @@ -7,6 +7,5 @@ #pragma link C++ class PFTaskManager+; #pragma link C++ class ConverterIn+; #pragma link C++ class ConverterOut+; -#pragma link C++ class RecoSimMatch+; #endif diff --git a/AnalysisTreeInterface/PFTaskManager.cpp b/AnalysisTreeInterface/PFTaskManager.cpp index 54b4699..aa0dce3 100644 --- a/AnalysisTreeInterface/PFTaskManager.cpp +++ b/AnalysisTreeInterface/PFTaskManager.cpp @@ -19,10 +19,9 @@ void PFTaskManager::Run(long long int nEvents) { FCFinder.FindParticles(); auto* converter_out = ((ConverterOut*)tasks_.at(kOutConverter)); - converter_out->SetIsMatchWithMc(true); - converter_out->InitSimAT("1.analysistree.root", "aTree"); // TODO avoid hardcoding converter_out->SetCandidates(FCFinder.GetLambdas()); converter_out->Exec(); +// out_tree_->Fill(); } // Event loop } diff --git a/AnalysisTreeInterface/RecoSimMatch.cpp b/AnalysisTreeInterface/RecoSimMatch.cpp deleted file mode 100644 index fd1fcc6..0000000 --- a/AnalysisTreeInterface/RecoSimMatch.cpp +++ /dev/null @@ -1,259 +0,0 @@ -#include "RecoSimMatch.h" - -void RecoSimMatch::InitSimAT(const std::string sim_file_name, const std::string sim_tree_name) -{ - in_file_sim_ = TFile::Open(sim_file_name.c_str(), "read"); - config_sim_ = (AnalysisTree::Configuration*) in_file_sim_->Get("Configuration"); - in_chain_sim_ = new TChain(sim_tree_name.c_str()); - in_chain_sim_ -> Add(sim_file_name.c_str()); - in_chain_sim_ -> SetBranchAddress("KfpfTracks", &kf_tracks_); - in_chain_sim_ -> SetBranchAddress("SimTracks", &sim_tracks_); - in_chain_sim_ -> SetBranchAddress("RecEventHeader", &rec_event_header_); - in_chain_sim_ -> SetBranchAddress("SimEventHeader", &sim_event_header_); - in_chain_sim_ -> SetBranchAddress(config_sim_->GetMatchName("KfpfTracks", "SimTracks").c_str(), &kf2sim_tracks_); - - auto branch_conf_kftr = config_sim_->GetBranchConfig( "KfpfTracks" ); - par_field_id_ = branch_conf_kftr.GetFieldId("x"); - - auto branch_conf_simtr = config_sim_->GetBranchConfig( "SimTracks" ); - mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); -} - -void RecoSimMatch::InitRecoAT(const std::string reco_file_name, const std::string reco_tree_name) -{ - in_file_reco_ = TFile::Open(reco_file_name.c_str(), "read"); - in_chain_reco_ = new TChain(reco_tree_name.c_str()); - in_chain_reco_ -> Add(reco_file_name.c_str()); - in_chain_reco_ -> SetBranchAddress("LambdaCandidates", &reco_tracks_); - - config_reco_ = (AnalysisTree::Configuration*) in_file_reco_->Get("Configuration"); - - auto branch_conf_recotr = config_reco_->GetBranchConfig("LambdaCandidates"); - - mass_field_id_ = branch_conf_recotr.GetFieldId("invmass"); - rap_lab_field_id_ = branch_conf_recotr.GetFieldId("rapidity_lab"); - rap_cm_field_id_ = branch_conf_recotr.GetFieldId("rapidity_cm"); - pdg_field_id_ = branch_conf_recotr.GetFieldId("pdg"); - x_field_id_ = branch_conf_recotr.GetFieldId("x"); - y_field_id_ = branch_conf_recotr.GetFieldId("y"); - z_field_id_ = branch_conf_recotr.GetFieldId("z"); - - daughter1_id_field_id_ = branch_conf_recotr.GetFieldId("daughter1id"); - daughter2_id_field_id_ = branch_conf_recotr.GetFieldId("daughter2id"); - - chi2primpos_field_id_ = branch_conf_recotr.GetFieldId("chi2primpos"); - chi2primneg_field_id_ = branch_conf_recotr.GetFieldId("chi2primneg"); - distance_field_id_ = branch_conf_recotr.GetFieldId("distance"); - cosinepos_field_id_ = branch_conf_recotr.GetFieldId("cosinepos"); - cosineneg_field_id_ = branch_conf_recotr.GetFieldId("cosineneg"); - chi2geo_field_id_ = branch_conf_recotr.GetFieldId("chi2geo"); - l_field_id_ = branch_conf_recotr.GetFieldId("l"); - ldl_field_id_ = branch_conf_recotr.GetFieldId("ldl"); - isfrompv_field_id_ = branch_conf_recotr.GetFieldId("isfrompv"); - cosinetopo_field_id_ = branch_conf_recotr.GetFieldId("cosinetopo"); - chi2topo_field_id_ = branch_conf_recotr.GetFieldId("chi2topo"); -} - -void RecoSimMatch::InitOutput(const std::string out_file_name) -{ - out_file_ = TFile::Open(out_file_name.c_str(), "recreate"); - -// ***** Lambda Candidates ******* - - AnalysisTree::BranchConfig LambdaRecoBranch("LambdaCandidates", AnalysisTree::DetType::kParticle); - - LambdaRecoBranch.AddField("chi2primpos"); - LambdaRecoBranch.AddField("chi2primneg"); - LambdaRecoBranch.AddField("distance"); - LambdaRecoBranch.AddField("cosinepos"); - LambdaRecoBranch.AddField("cosineneg"); - LambdaRecoBranch.AddField("chi2geo"); - LambdaRecoBranch.AddField("l"); - LambdaRecoBranch.AddField("ldl"); - LambdaRecoBranch.AddField ("isfrompv"); - LambdaRecoBranch.AddField("cosinetopo"); - LambdaRecoBranch.AddField("chi2topo"); - LambdaRecoBranch.AddField ("iszok"); - - LambdaRecoBranch.AddField("x"); - LambdaRecoBranch.AddField("y"); - LambdaRecoBranch.AddField("z"); - LambdaRecoBranch.AddField("invmass"); - LambdaRecoBranch.AddField("rapidity_lab"); - LambdaRecoBranch.AddField("rapidity_cm"); - LambdaRecoBranch.AddField("pdg"); - LambdaRecoBranch.AddField("daughter1id"); - LambdaRecoBranch.AddField("daughter2id"); - LambdaRecoBranch.AddField("issignal"); - - config_out_.AddBranchConfig( LambdaRecoBranch ); - lambda_reco_ = new AnalysisTree::Particles(config_out_.GetLastId()); - -// ***** Lambda simulated ******* - - AnalysisTree::BranchConfig LambdaSimBranch("LambdaSimulated", AnalysisTree::DetType::kParticle); - - LambdaSimBranch.AddField("x"); - LambdaSimBranch.AddField("y"); - LambdaSimBranch.AddField("z"); - LambdaSimBranch.AddField("invmass"); - LambdaSimBranch.AddField("rapidity_lab"); - LambdaSimBranch.AddField("rapidity_cm"); - LambdaSimBranch.AddField("pdg"); - LambdaSimBranch.AddField("simid"); - LambdaSimBranch.AddField("daughter1id"); - LambdaSimBranch.AddField("daughter2id"); - - config_out_.AddBranchConfig( LambdaSimBranch ); - lambda_sim_ = new AnalysisTree::Particles(config_out_.GetLastId()); - -// *********** Matching ********** - lambda_reco2sim_ = new AnalysisTree::Matching(config_out_.GetBranchConfig("LambdaCandidates").GetId(), config_out_.GetBranchConfig("LambdaSimulated").GetId()); - config_out_.AddMatch("LambdaCandidates", "LambdaSimulated", "LambdaCand2LambdaSim"); - - - out_tree_ = new TTree("aTree", "AnalysisTree Lambda"); - out_tree_ -> Branch("LambdaCandidates", "AnalysisTree::Particles", &lambda_reco_); - out_tree_ -> Branch("LambdaSimulated", "AnalysisTree::Particles", &lambda_sim_); - out_tree_ -> Branch("LambdaCand2LambdaSim", "AnalysisTree::Matching", &lambda_reco2sim_); - - config_out_.Write("Configuration"); - - x_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("x"); - y_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("y"); - z_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("z"); - mass_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("invmass"); - rap_lab_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_lab"); - rap_cm_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("rapidity_cm"); - pdg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("pdg"); - daughter1_id_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter1id"); - daughter2_id_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("daughter2id"); - is_signal_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("issignal"); - - chi2primpos_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primpos"); - chi2primneg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2primneg"); - distance_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("distance"); - cosinepos_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinepos"); - cosineneg_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosineneg"); - chi2geo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2geo"); - l_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("l"); - ldl_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("ldl"); - isfrompv_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("isfrompv"); - cosinetopo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("cosinetopo"); - chi2topo_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("chi2topo"); - iszok_field_id_w_ = config_out_.GetBranchConfig( lambda_reco_->GetId() ).GetFieldId("iszok"); - - // ------------------------------- - mass_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("invmass"); - rap_lab_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("rapidity_lab"); - rap_cm_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("rapidity_cm"); - pdg_sim_field_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("pdg"); - sim_id_w_ = config_out_.GetBranchConfig( lambda_sim_->GetId() ).GetFieldId("simid"); -} - -void RecoSimMatch::Run(int nEvents) -{ - if(in_chain_reco_->GetEntries() != in_chain_sim_->GetEntries()) - std::cout << "Error! Files are not matchable!\n"; - if (nEvents == -1 || nEvents > in_chain_reco_->GetEntries()) - nEvents = in_chain_reco_->GetEntries(); - - for(int iEvent=0; iEventGetEntry(iEvent); - in_chain_sim_->GetEntry(iEvent); - - lambda_reco_->ClearChannels(); - lambda_sim_->ClearChannels(); - lambda_reco2sim_->Clear(); - - for(int iCh=0; iChGetNumberOfChannels(); iCh++) - { - const AnalysisTree::Particle& recotrack = reco_tracks_->GetChannel(iCh); - - auto* lambdarec = lambda_reco_->AddChannel(); - lambdarec->Init(config_out_.GetBranchConfig(lambda_reco_->GetId())); - lambdarec->SetField(recotrack.GetField(mass_field_id_), mass_field_id_w_); - lambdarec->SetField(recotrack.GetField(daughter1_id_field_id_), daughter1_id_field_id_w_); - lambdarec->SetField(recotrack.GetField(daughter2_id_field_id_), daughter2_id_field_id_w_); - - lambdarec->SetMomentum(recotrack.GetPx(), recotrack.GetPy(), recotrack.GetPz()); - lambdarec->SetField(recotrack.GetField(rap_lab_field_id_), rap_lab_field_id_w_); - lambdarec->SetField(recotrack.GetField(rap_cm_field_id_), rap_cm_field_id_w_); - lambdarec->SetField(recotrack.GetField(pdg_field_id_), pdg_field_id_w_ ); - - lambdarec->SetPid(recotrack.GetField(pdg_field_id_)); - lambdarec->SetMass(recotrack.GetMass()); - lambdarec->SetField(recotrack.GetField(x_field_id_), x_field_id_w_); - lambdarec->SetField(recotrack.GetField(y_field_id_), y_field_id_w_); - lambdarec->SetField(recotrack.GetField(z_field_id_), z_field_id_w_); - - lambdarec->SetField(recotrack.GetField(chi2primpos_field_id_), chi2primpos_field_id_w_); - lambdarec->SetField(recotrack.GetField(chi2primneg_field_id_), chi2primneg_field_id_w_); - lambdarec->SetField(recotrack.GetField(distance_field_id_), distance_field_id_w_); - lambdarec->SetField(recotrack.GetField(cosinepos_field_id_), cosinepos_field_id_w_); - lambdarec->SetField(recotrack.GetField(cosineneg_field_id_), cosineneg_field_id_w_); - lambdarec->SetField(recotrack.GetField(chi2geo_field_id_), chi2geo_field_id_w_); - lambdarec->SetField(recotrack.GetField(l_field_id_), l_field_id_w_); - lambdarec->SetField(recotrack.GetField(ldl_field_id_), ldl_field_id_w_); - lambdarec->SetField(recotrack.GetField(isfrompv_field_id_), isfrompv_field_id_w_); - lambdarec->SetField(recotrack.GetField(cosinetopo_field_id_), cosinetopo_field_id_w_); - lambdarec->SetField(recotrack.GetField(chi2topo_field_id_), chi2topo_field_id_w_); - - const AnalysisTree::Track& reco_daughter1 = kf_tracks_->GetChannel(recotrack.GetField(daughter1_id_field_id_)); - const AnalysisTree::Track& reco_daughter2 = kf_tracks_->GetChannel(recotrack.GetField(daughter2_id_field_id_)); - - float z_daughter_min = std::min(reco_daughter1.GetField(par_field_id_+2), reco_daughter2.GetField(par_field_id_+2)); - if(recotrack.GetField(z_field_id_) < z_daughter_min) - lambdarec->SetField(1, iszok_field_id_w_); - else - lambdarec->SetField(0, iszok_field_id_w_); - - const int simtrackid1 = kf2sim_tracks_->GetMatch(recotrack.GetField(daughter1_id_field_id_)); - const int simtrackid2 = kf2sim_tracks_->GetMatch(recotrack.GetField(daughter2_id_field_id_)); - int is_signal = 0; - int mother_id = -999; - int mother_pdg = -1; - if(!(simtrackid1 < 0 || simtrackid2 < 0)) - { - const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); - const AnalysisTree::Particle& simtrack2 = sim_tracks_->GetChannel(simtrackid2); - if(simtrack1.GetField(mother_id_field_id_) == simtrack2.GetField(mother_id_field_id_)) - { - mother_id = simtrack1.GetField(mother_id_field_id_); - if(mother_id < 0) continue; - const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); - mother_pdg = simtrackmother.GetPid(); - if(mother_pdg == pdg_lambda) - is_signal = 1; - } - } - lambdarec->SetField( is_signal, is_signal_field_id_w_); - - if(is_signal!=1) continue; - - auto* lambdasim = lambda_sim_->AddChannel(); - lambdasim->Init(config_out_.GetBranchConfig(lambda_sim_->GetId())); - const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); - const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); - - lambdasim -> SetField(mass_lambda, mass_sim_field_id_w_); - lambdasim -> SetField(mother_pdg, pdg_sim_field_id_w_); - lambdasim -> SetPid(mother_pdg); - lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); - lambdasim -> SetField(simtrackmother.GetRapidity(), rap_lab_sim_field_id_w_); - lambdasim -> SetField(simtrackmother.GetRapidity() - RapidityShift(12.), rap_cm_sim_field_id_w_); - lambdasim -> SetField(mother_id, sim_id_w_); - - lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); - } - out_tree_->Fill(); - } - out_tree_->Write(); -} - -float RecoSimMatch::RapidityShift(float pbeam) -{ - const float nucleonmass = 0.939; - return TMath::Log(((sqrt(nucleonmass*nucleonmass + pbeam*pbeam) + pbeam) / (sqrt(nucleonmass*nucleonmass + pbeam*pbeam) - pbeam)))/4.; -} diff --git a/AnalysisTreeInterface/RecoSimMatch.h b/AnalysisTreeInterface/RecoSimMatch.h deleted file mode 100644 index 5d80c6f..0000000 --- a/AnalysisTreeInterface/RecoSimMatch.h +++ /dev/null @@ -1,124 +0,0 @@ -#include "TFile.h" -#include "TTree.h" -#include "TChain.h" - -#include "AnalysisTree/Track.h" -#include "AnalysisTree/Hit.h" -#include "AnalysisTree/Detector.h" -#include "AnalysisTree/Configuration.h" -#include "AnalysisTree/EventHeader.h" -#include "AnalysisTree/Matching.h" - -class RecoSimMatch -{ -public: - - // Constructors/Destructors --------- - RecoSimMatch() = default; - virtual ~RecoSimMatch() = default; - - void InitSimAT(const std::string sim_file_name, const std::string sim_tree_name); - void InitRecoAT(const std::string reco_file_name, const std::string reco_tree_name); - void InitOutput(const std::string out_file_name); - void Run(int nEvents=-1); - -private: - - float RapidityShift(float pbeam); - - TFile* in_file_sim_{nullptr}; - TFile* in_file_reco_{nullptr}; - TChain* in_chain_sim_{nullptr}; - TChain* in_chain_reco_{nullptr}; - TFile* out_file_{nullptr}; - TFile* out_file_tree_{nullptr}; - - AnalysisTree::Configuration* config_sim_{nullptr}; - AnalysisTree::Configuration* config_reco_{nullptr}; - AnalysisTree::Configuration config_out_; - - AnalysisTree::Particles* lambda_reco_{nullptr}; - AnalysisTree::Particles* lambda_sim_{nullptr}; - AnalysisTree::Matching* lambda_reco2sim_{nullptr}; - - TTree* out_tree_{nullptr}; - - //sims - AnalysisTree::TrackDetector* kf_tracks_{nullptr}; - AnalysisTree::Particles* sim_tracks_{nullptr}; - AnalysisTree::Matching* kf2sim_tracks_{nullptr}; - AnalysisTree::EventHeader* rec_event_header_{nullptr}; - AnalysisTree::EventHeader* sim_event_header_{nullptr}; - - //recos - AnalysisTree::Particles* reco_tracks_{nullptr}; - - //-----------INPUT field ids---------------------------------- - // field ids for input parameters of simulated lambda - int par_field_id_ {-1}; - int mother_id_field_id_{-1}; - - // field ids for input kinematic parameters of selected lambda candidates - int mass_field_id_{-1}; - int rap_lab_field_id_{-1}; - int rap_cm_field_id_{-1}; - int pdg_field_id_{-1}; - int x_field_id_{-1}; - int y_field_id_{-1}; - int z_field_id_{-1}; - - int daughter1_id_field_id_{-1}; - int daughter2_id_field_id_{-1}; - - // field ids for input selection parameters of lambda candidates (cutting) - int chi2primpos_field_id_{-1}; - int chi2primneg_field_id_{-1}; - int distance_field_id_{-1}; - int cosinepos_field_id_{-1}; - int cosineneg_field_id_{-1}; - int chi2geo_field_id_{-1}; - int l_field_id_{-1}; - int ldl_field_id_{-1}; - int isfrompv_field_id_{-1}; - int cosinetopo_field_id_{-1}; - int chi2topo_field_id_{-1}; - //------------------------------------------------------------ - - //--------OUTPUT field ids------------------------------------ - // field ids for selected lambda candidates kinematic parameters - int mass_field_id_w_{-1}; - int rap_lab_field_id_w_{-1}; - int rap_cm_field_id_w_{-1}; - int pdg_field_id_w_{-1}; - int daughter1_id_field_id_w_{-1}; - int daughter2_id_field_id_w_{-1}; - int is_signal_field_id_w_{-1}; - int x_field_id_w_{-1}; - int y_field_id_w_{-1}; - int z_field_id_w_{-1}; - - // field ids for simulated lamdas kinematic parameters - int mass_sim_field_id_w_{-1}; - int rap_lab_sim_field_id_w_{-1}; - int rap_cm_sim_field_id_w_{-1}; - int pdg_sim_field_id_w_{-1}; - int sim_id_w_{-1}; - - // field ids for lambda candidate cutting parameters - int chi2primpos_field_id_w_{-1}; - int chi2primneg_field_id_w_{-1}; - int distance_field_id_w_{-1}; - int cosinepos_field_id_w_{-1}; - int cosineneg_field_id_w_{-1}; - int chi2geo_field_id_w_{-1}; - int l_field_id_w_{-1}; - int ldl_field_id_w_{-1}; - int isfrompv_field_id_w_{-1}; - int cosinetopo_field_id_w_{-1}; - int chi2topo_field_id_w_{-1}; - int iszok_field_id_w_{-1}; - - const int pdg_lambda = 3122; - const float mass_lambda = 1.115683; - -}; \ No newline at end of file diff --git a/AnalysisTreeInterface/main.cxx b/AnalysisTreeInterface/main.cxx index 8d690de..4da6167 100644 --- a/AnalysisTreeInterface/main.cxx +++ b/AnalysisTreeInterface/main.cxx @@ -5,11 +5,11 @@ int main(int argc, char** argv) { if(argc < 2){ - std::cout << "Wrong number of arguments! Please use:\n ./main filename\n"; + std::cout << "Wrong number of arguments! Please use:\n ./main filelist.txt\n"; return EXIT_FAILURE; } - const bool make_plain_tree{false}; + const bool make_plain_tree{true}; CutsContainer cuts; cuts.CancelCuts(); @@ -29,9 +29,12 @@ int main(int argc, char** argv) in_converter->SetTrackCuts(new AnalysisTree::Cuts("Cut to reproduce KFPF", {{{"KfpfTracks", "pass_cuts"}, 1}})); in_converter->SetIsShine(false); //TODO maybe change name - man.AddTasks(in_converter, new ConverterOut); + auto* out_converter = new ConverterOut; + out_converter->SetInputBranchNames({"SimTracks", "KfpfTracks"}); + + man.AddTasks(in_converter, out_converter); man.Init(); - man.Run(-1); // -1 = all events + man.Run(10); // -1 = all events man.Finish(); if(make_plain_tree){ -- GitLab From a6b9a29f541d81b7a1fc157b86a737a2321ac9fd Mon Sep 17 00:00:00 2001 From: Lubynets Date: Sun, 26 Jul 2020 20:47:40 +0200 Subject: [PATCH 3/4] try no2 of mc-matching --- AnalysisTreeInterface/ConverterOut.cpp | 58 +++++++++++++++++++++----- AnalysisTreeInterface/ConverterOut.h | 11 +++-- AnalysisTreeInterface/main.cxx | 2 +- 3 files changed, 55 insertions(+), 16 deletions(-) diff --git a/AnalysisTreeInterface/ConverterOut.cpp b/AnalysisTreeInterface/ConverterOut.cpp index 34c0806..747af73 100644 --- a/AnalysisTreeInterface/ConverterOut.cpp +++ b/AnalysisTreeInterface/ConverterOut.cpp @@ -16,18 +16,13 @@ void ConverterOut::Exec() const KFParticle& particle = candidate.GetParticle(); float mass, masserr; particle.GetMass(mass, masserr); -// lambdarec->SetField(mass, mass_field_id_); lambdarec->SetField(particle.DaughterIds()[0], daughter1_id_field_id_); lambdarec->SetField(particle.DaughterIds()[1], daughter2_id_field_id_); -// std::cout << "FC\t" << particle.DaughterIds()[0] << "\t" << particle.DaughterIds()[1] << "\n"; lambdarec->SetField(particle.X(), x_field_id_); lambdarec->SetField(particle.Y(), y_field_id_); lambdarec->SetField(particle.Z(), z_field_id_); lambdarec->SetMomentum(particle.GetPx(), particle.GetPy(), particle.GetPz()); -// lambdarec->SetField( particle.GetRapidity(), rap_lab_field_id_); - lambdarec->SetField( particle.GetRapidity() - data_header_->GetBeamRapidity(), rap_cm_field_id_); -// lambdarec->SetField( particle.GetPDG(), pdg_field_id_w_ ); lambdarec->SetPid( particle.GetPDG()); lambdarec->SetMass(mass); @@ -57,6 +52,7 @@ void ConverterOut::Init(std::map& branches) if(!in_branches_.empty()){ mc_particles_ = (AnalysisTree::Particles*) branches.find(in_branches_[0])->second; rec_tracks_ = (AnalysisTree::TrackDetector*) branches.find(in_branches_[1])->second; + kf2sim_tracks_ = (AnalysisTree::Matching*) branches.find(in_branches_[2])->second; } AnalysisTree::BranchConfig LambdaRecoBranch(out_branch_, AnalysisTree::DetType::kParticle); @@ -76,7 +72,6 @@ void ConverterOut::Init(std::map& branches) LambdaRecoBranch.AddField("x"); LambdaRecoBranch.AddField("y"); LambdaRecoBranch.AddField("z"); - LambdaRecoBranch.AddField("rapidity_cm"); LambdaRecoBranch.AddField("daughter1id"); LambdaRecoBranch.AddField("daughter2id"); if(mc_particles_) { @@ -93,7 +88,47 @@ void ConverterOut::Init(std::map& branches) void ConverterOut::MatchWithMc(){ std::cout << mc_particles_->GetNumberOfChannels() << " " << rec_tracks_->GetNumberOfChannels() << std::endl; - + + for(int i_ch=0; i_chGetNumberOfChannels(); i_ch++){ + AnalysisTree::Particle* lambdarec = &(lambda_reco_->GetChannel(i_ch)); + + const int simtrackid1 = kf2sim_tracks_->GetMatch(lambdarec->GetField(daughter1_id_field_id_)); + const int simtrackid2 = kf2sim_tracks_->GetMatch(lambdarec->GetField(daughter2_id_field_id_)); + bool is_signal = false; + int mother_id = -999; + int mother_pdg = -1; + if(!(simtrackid1 < 0 || simtrackid2 < 0)) + { + const AnalysisTree::Particle& simtrack1 = mc_particles_->GetChannel(simtrackid1); + const AnalysisTree::Particle& simtrack2 = mc_particles_->GetChannel(simtrackid2); + if(simtrack1.GetField(mother_id_field_id_) == simtrack2.GetField(mother_id_field_id_)) + { + mother_id = simtrack1.GetField(mother_id_field_id_); + if(mother_id < 0) continue; + const AnalysisTree::Particle& simtrackmother = mc_particles_->GetChannel(mother_id); + mother_pdg = simtrackmother.GetPid(); + if(mother_pdg == pdg_lambda) + is_signal = true; + } + } + lambdarec->SetField( is_signal, is_signal_field_id_); +// +// if(is_signal!=1) continue; +// +// auto* lambdasim = lambda_sim_->AddChannel(); +// lambdasim->Init(config_out_.GetBranchConfig(lambda_sim_->GetId())); +// const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); +// const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); +// +// lambdasim -> SetField(mass_lambda, mass_sim_field_id_w_); +// lambdasim -> SetField(mother_pdg, pdg_sim_field_id_w_); +// lambdasim -> SetPid(mother_pdg); +// lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); +// lambdasim -> SetField(simtrackmother.GetRapidity(), rap_lab_sim_field_id_w_); +// lambdasim -> SetField(mother_id, sim_id_w_); +// +// lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); + } } @@ -104,12 +139,13 @@ void ConverterOut::InitIndexes(){ x_field_id_ = out_branch.GetFieldId("x"); y_field_id_ = out_branch.GetFieldId("y"); z_field_id_ = out_branch.GetFieldId("z"); -// mass_field_id_ = out_branch.GetFieldId("invmass"); -// rap_lab_field_id_ = out_branch.GetFieldId("rapidity_lab"); - rap_cm_field_id_ = out_branch.GetFieldId("rapidity_cm"); -// pdg_field_id_w_ = out_branch.GetFieldId("pdg"); daughter1_id_field_id_ = out_branch.GetFieldId("daughter1id"); daughter2_id_field_id_ = out_branch.GetFieldId("daughter2id"); + if(mc_particles_){ + is_signal_field_id_ = out_branch.GetFieldId("is_signal"); + auto branch_conf_kftr = config_->GetBranchConfig(in_branches_[0]); + mother_id_field_id_ = branch_conf_kftr.GetFieldId("mother_id"); + } chi2primpos_field_id_ = out_branch.GetFieldId("chi2primpos"); chi2primneg_field_id_ = out_branch.GetFieldId("chi2primneg"); diff --git a/AnalysisTreeInterface/ConverterOut.h b/AnalysisTreeInterface/ConverterOut.h index f24b418..0555205 100644 --- a/AnalysisTreeInterface/ConverterOut.h +++ b/AnalysisTreeInterface/ConverterOut.h @@ -26,19 +26,20 @@ class ConverterOut : public AnalysisTree::FillTask { AnalysisTree::Particles* mc_particles_{nullptr}; AnalysisTree::TrackDetector* rec_tracks_{nullptr}; + AnalysisTree::Matching* kf2sim_tracks_{nullptr}; std::vector canditates_; + // field ids of input simulated lambdas + int mother_id_field_id_{-1}; + // field ids for selected lambda candidates kinematic parameters int x_field_id_{-1}; int y_field_id_{-1}; int z_field_id_{-1}; - int mass_field_id_{-1}; -// int rap_lab_field_id_{-1}; - int rap_cm_field_id_{-1}; -// int pdg_field_id_w_{-1}; int daughter1_id_field_id_{-1}; int daughter2_id_field_id_{-1}; + int is_signal_field_id_{-1}; // field ids for lambda candidate cutting parameters int chi2primpos_field_id_{-1}; @@ -52,6 +53,8 @@ class ConverterOut : public AnalysisTree::FillTask { int isfrompv_field_id_{-1}; int cosinetopo_field_id_{-1}; int chi2topo_field_id_{-1}; + + const int pdg_lambda = 3122; }; #endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ diff --git a/AnalysisTreeInterface/main.cxx b/AnalysisTreeInterface/main.cxx index 4da6167..f83ce4b 100644 --- a/AnalysisTreeInterface/main.cxx +++ b/AnalysisTreeInterface/main.cxx @@ -30,7 +30,7 @@ int main(int argc, char** argv) in_converter->SetIsShine(false); //TODO maybe change name auto* out_converter = new ConverterOut; - out_converter->SetInputBranchNames({"SimTracks", "KfpfTracks"}); + out_converter->SetInputBranchNames({"SimTracks", "KfpfTracks", "KfpfTracks2SimTracks"}); man.AddTasks(in_converter, out_converter); man.Init(); -- GitLab From d3a6cb7ffabca71e3c2c976c8df12ac11a885975 Mon Sep 17 00:00:00 2001 From: Lubynets Date: Mon, 27 Jul 2020 18:41:30 +0200 Subject: [PATCH 4/4] matching LamCand with SimTracks added --- AnalysisTreeInterface/ConverterOut.cpp | 27 +++++++++----------------- AnalysisTreeInterface/ConverterOut.h | 4 ++++ AnalysisTreeInterface/main.cxx | 4 ++-- 3 files changed, 15 insertions(+), 20 deletions(-) diff --git a/AnalysisTreeInterface/ConverterOut.cpp b/AnalysisTreeInterface/ConverterOut.cpp index e5e8f5e..5ba89ca 100644 --- a/AnalysisTreeInterface/ConverterOut.cpp +++ b/AnalysisTreeInterface/ConverterOut.cpp @@ -7,6 +7,7 @@ void ConverterOut::Exec() { lambda_reco_->ClearChannels(); + lambda_reco2sim_->Clear(); for(const auto& candidate : canditates_) { @@ -81,14 +82,17 @@ void ConverterOut::Init(std::map& branches) out_config_->AddBranchConfig( LambdaRecoBranch ); lambda_reco_ = new AnalysisTree::Particles(out_config_->GetLastId()); + + lambda_reco2sim_ = new AnalysisTree::Matching(out_config_->GetBranchConfig(out_branch_.c_str()).GetId(), config_->GetBranchConfig(in_branches_[0]).GetId()); + out_config_->AddMatch(out_branch_.c_str(), in_branches_[0], "LambdaCand2LambdaSim"); + out_tree_ -> Branch(out_branch_.c_str(), "AnalysisTree::Particles", &lambda_reco_); + out_tree_ -> Branch("LambdaCand2LambdaSim", "AnalysisTree::Matching", &lambda_reco2sim_); InitIndexes(); } void ConverterOut::MatchWithMc(){ - -// std::cout << mc_particles_->GetNumberOfChannels() << " " << rec_tracks_->GetNumberOfChannels() << std::endl; for(int i_ch=0; i_chGetNumberOfChannels(); i_ch++){ AnalysisTree::Particle* lambdarec = &(lambda_reco_->GetChannel(i_ch)); @@ -113,22 +117,9 @@ void ConverterOut::MatchWithMc(){ } } lambdarec->SetField( is_signal, is_signal_field_id_); -// -// if(is_signal!=1) continue; -// -// auto* lambdasim = lambda_sim_->AddChannel(); -// lambdasim->Init(config_out_.GetBranchConfig(lambda_sim_->GetId())); -// const AnalysisTree::Particle& simtrackmother = sim_tracks_->GetChannel(mother_id); -// const AnalysisTree::Particle& simtrack1 = sim_tracks_->GetChannel(simtrackid1); -// -// lambdasim -> SetField(mass_lambda, mass_sim_field_id_w_); -// lambdasim -> SetField(mother_pdg, pdg_sim_field_id_w_); -// lambdasim -> SetPid(mother_pdg); -// lambdasim -> SetMomentum(simtrackmother.GetPx(), simtrackmother.GetPy(), simtrackmother.GetPz()); -// lambdasim -> SetField(simtrackmother.GetRapidity(), rap_lab_sim_field_id_w_); -// lambdasim -> SetField(mother_id, sim_id_w_); -// -// lambda_reco2sim_->AddMatch(lambdarec->GetId(), lambdasim->GetId()); + + if(is_signal==true) + lambda_reco2sim_->AddMatch(lambdarec->GetId(), mother_id); } } diff --git a/AnalysisTreeInterface/ConverterOut.h b/AnalysisTreeInterface/ConverterOut.h index 18c2fc4..67c90f6 100644 --- a/AnalysisTreeInterface/ConverterOut.h +++ b/AnalysisTreeInterface/ConverterOut.h @@ -22,12 +22,16 @@ class ConverterOut : public AnalysisTree::FillTask { void InitIndexes(); void MatchWithMc(); + // output branches AnalysisTree::Particles* lambda_reco_{nullptr}; + AnalysisTree::Matching* lambda_reco2sim_{nullptr}; + // input branches AnalysisTree::Particles* mc_particles_{nullptr}; AnalysisTree::TrackDetector* rec_tracks_{nullptr}; AnalysisTree::Matching* rec_to_mc_{nullptr}; + std::vector canditates_; // field ids of input simulated lambdas diff --git a/AnalysisTreeInterface/main.cxx b/AnalysisTreeInterface/main.cxx index f83ce4b..be78fd4 100644 --- a/AnalysisTreeInterface/main.cxx +++ b/AnalysisTreeInterface/main.cxx @@ -30,11 +30,11 @@ int main(int argc, char** argv) in_converter->SetIsShine(false); //TODO maybe change name auto* out_converter = new ConverterOut; - out_converter->SetInputBranchNames({"SimTracks", "KfpfTracks", "KfpfTracks2SimTracks"}); + out_converter->SetInputBranchNames({"SimTracks", "KfpfTracks"}); man.AddTasks(in_converter, out_converter); man.Init(); - man.Run(10); // -1 = all events + man.Run(-1); // -1 = all events man.Finish(); if(make_plain_tree){ -- GitLab