diff --git a/AnalysisTreeInterface/CMakeLists.txt b/AnalysisTreeInterface/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a73aaa0ab6c6ca2cce3de5a4505c2336e598060d --- /dev/null +++ b/AnalysisTreeInterface/CMakeLists.txt @@ -0,0 +1,53 @@ +set(SOURCES + ConverterIn.cpp + ConverterOut.cpp + PFTaskManager.cpp + ) + +string(REPLACE ".cpp" ".h" HEADERS "${SOURCES}") + +include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/Interface ${CMAKE_SOURCE_DIR}/KFSimple ${AnalysisTree_INCLUDE_DIR}) + +if(PROJECT_LINK_DIRECTORIES) + include_directories(${PROJECT_INCLUDE_DIRECTORIES}) + link_directories(${PROJECT_LINK_DIRECTORIES}) +endif() + +add_library(KFMan SHARED ${SOURCES} G__KFMan.cxx) + +ROOT_GENERATE_DICTIONARY(G__KFMan ${HEADERS} LINKDEF KFManLinkDef.h OPTIONS "-DDO_TPCCATRACKER_EFF_PERFORMANCE" "-DNonhomogeneousField" "-DCBM" "-DUSE_TIMERS") + +target_link_libraries(KFMan ${ROOT_LIBRARIES} AnalysisTreeBase AnalysisTreeInfra KFParticleInterface KFParticleSimple KFParticle) +add_dependencies(KFMan KFParticleSimple) +add_target_property(KFMan COMPILE_FLAGS "-DDO_TPCCATRACKER_EFF_PERFORMANCE -DNonhomogeneousField -DCBM -DUSE_TIMERS") + +add_executable(main main.cxx) +add_dependencies(main KFMan) +target_link_libraries(main KFMan KFParticleSimple KFParticleInterface KFParticle Vc) + +install(TARGETS KFMan EXPORT KFManTargets + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib + RUNTIME DESTINATION bin + INCLUDES DESTINATION include + ) + +install( + FILES + ${HEADERS} + DESTINATION + include + COMPONENT + Devel +) + +set(PCM_FILE_NAME libKFMan) + +install( + FILES + "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}_rdict.pcm" + DESTINATION + lib + OPTIONAL +) +install (TARGETS main RUNTIME DESTINATION bin) \ No newline at end of file diff --git a/AnalysisTreeInterface/ConverterIn.cpp b/AnalysisTreeInterface/ConverterIn.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a04d33ac25301ae8c1dbd74cb814c78a22bc464c --- /dev/null +++ b/AnalysisTreeInterface/ConverterIn.cpp @@ -0,0 +1,164 @@ +#include +#include +#include "ConverterIn.h" + +void ConverterIn::FillParticle(const AnalysisTree::Track& rec_particle, InputContainer& input_info) const +{ + std::vector mf(kNumberOfFieldPars, 0.f); + for(int iF=0; iF(mf_field_id_+iF); + } + auto cov_matrix = is_shine_ ? GetCovMatrixShine(rec_particle) : GetCovMatrixCbm(rec_particle); + std::vector par(kNumberOfTrackPars,0.f); + + par[kX] = rec_particle.GetField(par_field_id_); + par[kY] = rec_particle.GetField(par_field_id_ + 1); + par[kZ] = rec_particle.GetField(par_field_id_ + 2); + par[kPx] = rec_particle.GetPx(); + par[kPy] = rec_particle.GetPy(); + par[kPz] = rec_particle.GetPz(); + + const int pdg = rec_particle.GetField(pdg_field_id_); //TODO +// const int pdg = rec_particle.GetPid(); + + input_info.AddTrack(par, cov_matrix, mf, rec_particle.GetField(q_field_id_), pdg, rec_particle.GetId()); +} + + +void ConverterIn::Init(std::map& branches) { + rec_event_header_ = (AnalysisTree::EventHeader*) branches.find(in_branches_[kRecEventHeader])->second; + sim_event_header_ = (AnalysisTree::EventHeader*) branches.find(in_branches_[kSimEventHeader])->second; + kf_tracks_ = (AnalysisTree::TrackDetector*) branches.find(in_branches_[kKfpfTracks])->second; + sim_tracks_ = (AnalysisTree::Particles*) branches.find(in_branches_[kSimTracks])->second; + kf2sim_tracks_ = (AnalysisTree::Matching*) branches.find( + config_->GetMatchName(in_branches_[kKfpfTracks], in_branches_[kSimTracks]))->second; + + auto branch_conf_kftr = config_->GetBranchConfig(in_branches_[kKfpfTracks]); + + q_field_id_ = branch_conf_kftr.GetFieldId("q"); + par_field_id_ = branch_conf_kftr.GetFieldId("x"); // par0 + mf_field_id_ = branch_conf_kftr.GetFieldId("cx0"); // magnetic field par0 + cov_field_id_ = branch_conf_kftr.GetFieldId("cov1"); // cov matrix 0 + passcuts_field_id_ = branch_conf_kftr.GetFieldId("pass_cuts"); + pdg_field_id_ = branch_conf_kftr.GetFieldId("mc_pdg"); + + auto branch_conf_simtr = config_->GetBranchConfig(in_branches_[kSimTracks]); + mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); + sim_pdg_field_id_ = branch_conf_simtr.GetFieldId("pdg"); + + if (track_cuts_){ + track_cuts_->Init(*config_); + } +} + +InputContainer ConverterIn::CreateInputContainer() const { + + InputContainer input_container; + const int n_tracks = kf_tracks_->GetNumberOfChannels(); + std::cout << " Ntracks = " << n_tracks << std::endl; + input_container.SetCuts(cuts_); + input_container.SetPV(rec_event_header_->GetVertexX(), rec_event_header_->GetVertexY(), rec_event_header_->GetVertexZ()); + + int n_good_tracks{0}; + input_container.Reserve(n_tracks); + for(int i_track=0; i_trackGetChannel(i_track); + if(!IsGoodTrack(rec_track)) continue; + FillParticle(rec_track, input_container); + n_good_tracks++; + } + std::cout << "Good tracks = " << n_good_tracks << "\n" << std::endl; + return input_container; +} + +std::vector ConverterIn::GetCovMatrixCbm(const AnalysisTree::Track& particle) const { + const auto tx = particle.GetField(par_field_id_ + 3); + const auto ty = particle.GetField(par_field_id_ + 4); + const auto qp = particle.GetField(par_field_id_ + 5); + const auto q = particle.GetField(q_field_id_); + + //calculate covariance matrix + const auto t = sqrt(1.f + tx * tx + ty * ty); + const auto t3 = t * t * t; + const auto dpxdtx = q / qp * (1.f + ty * ty) / t3; + const auto dpxdty = -q / qp * tx * ty / t3; + const auto dpxdqp = -q / (qp * qp) * tx / t; + const auto dpydtx = -q / qp * tx * ty / t3; + const auto dpydty = q / qp * (1.f + tx * tx) / t3; + const auto dpydqp = -q / (qp * qp) * ty / t; + const auto dpzdtx = -q / qp * tx / t3; + const auto dpzdty = -q / qp * ty / t3; + const auto dpzdqp = -q / (qp * qp * t3); + + const float F[kNumberOfTrackPars][5] = {{1.f, 0.f, 0.f, 0.f, 0.f}, + {0.f, 1.f, 0.f, 0.f, 0.f}, + {0.f, 0.f, 0.f, 0.f, 0.f}, + {0.f, 0.f, dpxdtx, dpxdty, dpxdqp}, + {0.f, 0.f, dpydtx, dpydty, dpydqp}, + {0.f, 0.f, dpzdtx, dpzdty, dpzdqp}}; + + float VFT[5][kNumberOfTrackPars]; + for(int i = 0; i < 5; i++){ + for(int j = 0; j < kNumberOfTrackPars; j++) { + VFT[i][j] = 0; + for(int k = 0; k < 5; k++) { + VFT[i][j] += particle.GetField(cov_field_id_ + std::min(i, k) + std::max(i, k) * (std::max(i, k) + 1) / 2) * F[j][k]; //parameters->GetCovariance(i,k) * F[j][k]; +// if(k <= i) +// VFT[i][j] += particle.GetField(cov_field_id_ + k + i * (i + 1) / 2) * F[j][k]; //parameters->GetCovariance(i,k) * F[j][k]; +// else +// VFT[i][j] += particle.GetField(cov_field_id_ + i + k * (k + 1) / 2) * F[j][k]; //parameters->GetCovariance(i,k) * F[j][k]; + } + } + } + + std::vector cov(21, 0); + for(int i=0, l=0; i ConverterIn::GetCovMatrixShine(const AnalysisTree::Track& particle) const +{ + std::vector cov(21, 0.); + + for (int iCov=0; iCov<21; ++iCov) + cov[iCov] = particle.GetField(cov_field_id_+iCov); + + return cov; +} + +bool ConverterIn::IsGoodTrack(const AnalysisTree::Track& rec_track) const +{ + if(!track_cuts_){ + return true; + } + return track_cuts_->Apply(rec_track); + +// return true; + bool is_good{false}; + + const int sim_id = kf2sim_tracks_->GetMatch(rec_track.GetId()); +// std::cout<< "sim " << sim_id << std::endl; + if (sim_id >= 0 && sim_id < sim_tracks_->GetNumberOfChannels()) { + const AnalysisTree::Track& sim_track = sim_tracks_->GetChannel(sim_id); + const int mother_id = sim_track.GetField(mother_id_field_id_); +// std::cout << "mother " << mother_id << std::endl; + + if (mother_id >= 0 && mother_id < sim_tracks_->GetNumberOfChannels() ) + { + const AnalysisTree::Track& mother_track = sim_tracks_->GetChannel(mother_id); + const int mother_pdg = mother_track.GetField(sim_pdg_field_id_); + std::cout << "mother pdg " << mother_pdg << std::endl; + + if(mother_pdg == 3122) + is_good = true; + } + } + return is_good; +} diff --git a/AnalysisTreeInterface/ConverterIn.h b/AnalysisTreeInterface/ConverterIn.h new file mode 100644 index 0000000000000000000000000000000000000000..aef4ebc8a963c4f72c7a5384485cd326dccde960 --- /dev/null +++ b/AnalysisTreeInterface/ConverterIn.h @@ -0,0 +1,79 @@ +#ifndef KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTERIN_H_ +#define KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTERIN_H_ + +#include +#include +#include +#include +#include "AnalysisTree/FillTask.hpp" + +#include "Interface/InputContainer.h" + +class ConverterIn : public AnalysisTree::FillTask { + + enum eInputBranches { + kRecEventHeader = 0, + kSimEventHeader, + kKfpfTracks, + kSimTracks, + kNumberOfInputBranches + }; + + public: + explicit ConverterIn(const CutsContainer& cuts) : cuts_(cuts) + { // Use default names for the input branches + in_branches_.resize(kNumberOfInputBranches); + in_branches_[kRecEventHeader] = "RecEventHeader"; + in_branches_[kSimEventHeader] = "SimEventHeader"; + in_branches_[kKfpfTracks] = "KfpfTracks"; + in_branches_[kSimTracks] = "SimTracks"; + } + ~ConverterIn() override = default; + + void Init(std::map& branches) override; + + void Exec() override {}; + + void Finish() override {}; + + const CutsContainer& GetCuts() const { return cuts_; } +// const InputContainer& GetInputContainer() const { return input_container_; } + + void SetIsShine(bool is=true) { is_shine_ = is; } + + void SetCuts(const CutsContainer& cuts) { cuts_ = cuts; }; + void SetTrackCuts(AnalysisTree::Cuts* const cuts) { track_cuts_ = cuts; }; + + InputContainer CreateInputContainer() const; + protected: + + std::vector GetCovMatrixCbm(const AnalysisTree::Track&) const; + std::vector GetCovMatrixShine(const AnalysisTree::Track&) const; + void FillParticle(const AnalysisTree::Track&, InputContainer&) const; + bool IsGoodTrack(const AnalysisTree::Track& rec_track) const; + + 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}; + AnalysisTree::Cuts* track_cuts_{nullptr}; + + CutsContainer cuts_; +// InputContainer input_container_; + + // field ids for input parameters + int q_field_id_{AnalysisTree::UndefValueInt}; + int pdg_field_id_{AnalysisTree::UndefValueInt}; + int par_field_id_ {AnalysisTree::UndefValueInt}; + int mf_field_id_ {AnalysisTree::UndefValueInt}; + int cov_field_id_{AnalysisTree::UndefValueInt}; + int mother_id_field_id_{AnalysisTree::UndefValueInt}; + int sim_pdg_field_id_{AnalysisTree::UndefValueInt}; + int passcuts_field_id_{AnalysisTree::UndefValueInt}; + + bool is_shine_{false}; + +}; + +#endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTERIN_H_ diff --git a/AnalysisTreeInterface/ConverterOut.cpp b/AnalysisTreeInterface/ConverterOut.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d1ff09127d4463c7be6ccc071f231a410a9f85f2 --- /dev/null +++ b/AnalysisTreeInterface/ConverterOut.cpp @@ -0,0 +1,104 @@ +#include "ConverterOut.h" + +#include "TTree.h" + +#include "AnalysisTree/DataHeader.hpp" + +void ConverterOut::Exec() +{ + lambda_reco_->ClearChannels(); + + for(const auto& candidate : canditates_) + { + auto* lambdarec = lambda_reco_->AddChannel(); + lambdarec->Init(out_config_->GetBranchConfig(lambda_reco_->GetId())); + + 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); + + lambdarec->SetField( candidate.GetChi2PrimPos(), chi2primpos_field_id_); + lambdarec->SetField( candidate.GetChi2PrimNeg(), chi2primneg_field_id_); + lambdarec->SetField( candidate.GetDistance(), distance_field_id_); + lambdarec->SetField( candidate.GetCosineDaughterPos(), cosinepos_field_id_); + lambdarec->SetField( candidate.GetCosineDaughterNeg(), cosineneg_field_id_); + lambdarec->SetField( candidate.GetChi2Geo(), chi2geo_field_id_); + lambdarec->SetField( candidate.GetL(), l_field_id_); + lambdarec->SetField( candidate.GetLdL(), ldl_field_id_); + lambdarec->SetField( candidate.GetIsFromPV(), isfrompv_field_id_); + lambdarec->SetField( candidate.GetCosineTopo(), cosinetopo_field_id_); + lambdarec->SetField( candidate.GetChi2Topo(), chi2topo_field_id_); + } + out_tree_->Fill(); +} + + +void ConverterOut::Init(std::map& branches) +{ + assert(out_config_ && out_tree_); + + 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("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"); + + out_config_->AddBranchConfig( LambdaRecoBranch ); + lambda_reco_ = new AnalysisTree::Particles(out_config_->GetLastId()); + out_tree_ -> Branch("LambdaCandidates", "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"); + 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"); + + 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"); +} diff --git a/AnalysisTreeInterface/ConverterOut.h b/AnalysisTreeInterface/ConverterOut.h new file mode 100644 index 0000000000000000000000000000000000000000..1e5eda3c046ab5985d548b2faeb37ece33c06869 --- /dev/null +++ b/AnalysisTreeInterface/ConverterOut.h @@ -0,0 +1,50 @@ +#ifndef KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ +#define KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ + +#include + +#include "AnalysisTree/FillTask.hpp" +#include "AnalysisTree/Detector.hpp" + +class ConverterOut : public AnalysisTree::FillTask { + public: + ConverterOut() = default; + ~ConverterOut() override = default; + + void Init(std::map& branches) override; + void Exec() override; + void Finish() override {} + + void SetCandidates(const std::vector& canditates) { canditates_ = canditates; } + + protected: + + AnalysisTree::Particles* lambda_reco_{nullptr}; + std::vector canditates_; + + // 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}; + + // field ids for lambda candidate cutting parameters + 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}; +}; + +#endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_CONVERTEROUT_H_ diff --git a/AnalysisTreeInterface/KFManLinkDef.h b/AnalysisTreeInterface/KFManLinkDef.h new file mode 100644 index 0000000000000000000000000000000000000000..306dd625dd3432642c7ee88d1ff1396a12d99f47 --- /dev/null +++ b/AnalysisTreeInterface/KFManLinkDef.h @@ -0,0 +1,11 @@ +#ifdef __MAKECINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class PFTaskManager+; +#pragma link C++ class ConverterIn+; +#pragma link C++ class ConverterOut+; + +#endif diff --git a/AnalysisTreeInterface/PFTaskManager.cpp b/AnalysisTreeInterface/PFTaskManager.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6dd7f14a0e361406ce51acdd2181aa0c709b8f7b --- /dev/null +++ b/AnalysisTreeInterface/PFTaskManager.cpp @@ -0,0 +1,27 @@ +#include +#include "PFTaskManager.h" + +void PFTaskManager::Run(long long int nEvents) { + std::cout << "PFTaskManager::Run" << std::endl; + nEvents = nEvents < 0 ? in_tree_->GetEntries() : nEvents; + + for (long long iEvent = 0; iEvent < nEvents; ++iEvent) { + in_tree_->GetEntry(iEvent); + std::cout << "Event # " << iEvent << " out of " << nEvents << std::endl; + + auto* converter_in = ((ConverterIn*)tasks_.at(kInConverter)); + SetCuts(converter_in->GetCuts()); + + SimpleFinder FCFinder; + FCFinder.Init(converter_in->CreateInputContainer()); + FCFinder.SetCuts(converter_in->GetCuts()); + FCFinder.SortTracks(); + FCFinder.FindParticles(); + + auto* converter_out = ((ConverterOut*)tasks_.at(kOutConverter)); + converter_out->SetCandidates(FCFinder.GetLambdas()); + converter_out->Exec(); + +// out_tree_->Fill(); + } // Event loop +} diff --git a/AnalysisTreeInterface/PFTaskManager.h b/AnalysisTreeInterface/PFTaskManager.h new file mode 100644 index 0000000000000000000000000000000000000000..efcc197ec8266d81888efb59e4ebb2893c159d14 --- /dev/null +++ b/AnalysisTreeInterface/PFTaskManager.h @@ -0,0 +1,39 @@ +#ifndef KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_PFTASKMANAGER_H_ +#define KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_PFTASKMANAGER_H_ + +#include "AnalysisTree/TaskManager.hpp" + +#include "ConverterOut.h" +#include "ConverterIn.h" + +class PFTaskManager : public AnalysisTree::TaskManager { + + enum eTasks{ + kInConverter = 0, + kOutConverter + }; + + public: + PFTaskManager(const std::string& filelist, const std::string& in_tree) : + AnalysisTree::TaskManager({filelist}, {in_tree}) {} + + void Run(long long nEvents) final; + + void AddTasks(ConverterIn* in_task, ConverterOut* out_task){ + assert(tasks_.empty()); + tasks_.emplace_back(in_task); + tasks_.emplace_back(out_task); + } + + void SetCuts(const CutsContainer& cuts) { cuts_ = cuts; }; + + + void AddTask(AnalysisTree::FillTask *task) = delete; //TODO make it virtual in AT + + protected: + CutsContainer cuts_; + + +}; + +#endif //KFPARTICLESIMPLE_ANALYSISTREEINTERFACE_PFTASKMANAGER_H_ diff --git a/AnalysisTreeInterface/main.cxx b/AnalysisTreeInterface/main.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6daa673827692160ff3a4f084554a59c49573c22 --- /dev/null +++ b/AnalysisTreeInterface/main.cxx @@ -0,0 +1,34 @@ +#include "PFTaskManager.h" + +int main(int argc, char** argv) +{ + if(argc < 2){ + std::cout << "Wrong number of arguments! Please use:\n ./main filename\n"; + return EXIT_FAILURE; + } + + CutsContainer cuts; + cuts.CancelCuts(); + cuts.SetCutChi2PrimPos(3.); + cuts.SetCutChi2PrimNeg(3.); + cuts.SetCutDistance(1.); + cuts.SetCutChi2Geo(3.); + cuts.SetCutLdL(5.); + + const std::string& filename = argv[1]; + + PFTaskManager man(filename, "aTree"); + man.SetOutFileName("PFSimpleOutput.root"); + man.SetOutTreeName("sTree"); + + auto* in_converter = new ConverterIn(cuts); + 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); + man.Init(); + man.Run(-1); // -1 = all events + man.Finish(); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/AnalysisTreeInterface/old/InConverter.cxx b/AnalysisTreeInterface/old/InConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..646905ee1bca8d6c8d355e779dd1afe68d156692 --- /dev/null +++ b/AnalysisTreeInterface/old/InConverter.cxx @@ -0,0 +1,231 @@ +#include "InConverter.h" + +//KF Particle headers +#include "KFParticle.h" +#include "KFPTrackVector.h" + +//KF Simple headers +#include "SimpleFinder.h" +#include "Constants.h" + +//c++ and std headers +#include +#include +#include + +void InConverter::InitAnalysisTree(const std::string& file_name, const std::string& tree_name) +{ + in_file_ = TFile::Open(file_name.c_str(), "read"); + config_ = (AnalysisTree::Configuration*) in_file_->Get("Configuration"); + + in_chain_ = new TChain(tree_name.c_str()); + in_chain_ -> Add(file_name.c_str()); + in_chain_ -> SetBranchAddress("KfpfTracks", &kf_tracks_); + in_chain_ -> SetBranchAddress("SimTracks", &sim_tracks_); + in_chain_ -> SetBranchAddress("RecEventHeader", &rec_event_header_); + in_chain_ -> SetBranchAddress("SimEventHeader", &sim_event_header_); + in_chain_ -> SetBranchAddress(config_->GetMatchName("KfpfTracks", "SimTracks").c_str(), &kf2sim_tracks_); + + auto branch_conf_kftr = config_->GetBranchConfig( "KfpfTracks" ); + q_field_id_ = branch_conf_kftr.GetFieldId("q"); + + par_field_id_ = branch_conf_kftr.GetFieldId("x"); // par0 + mf_field_id_ = branch_conf_kftr.GetFieldId("cx0"); // magnetic field par0 + cov_field_id_ = branch_conf_kftr.GetFieldId("cov1"); // cov matrix 0 + + passcuts_field_id_ = branch_conf_kftr.GetFieldId("pass_cuts"); + pdg_field_id_ = branch_conf_kftr.GetFieldId("mc_pdg"); + + auto branch_conf_simtr = config_->GetBranchConfig( "SimTracks" ); + mother_id_field_id_ = branch_conf_simtr.GetFieldId("mother_id"); + sim_pdg_field_id_ = branch_conf_simtr.GetFieldId("pdg"); + + if (track_cuts_ != nullptr) + track_cuts_->Init(*config_); +} + +int InConverter::GetTrackPid(const AnalysisTree::Track& rec_track) +{ + int pdg = -1; + if(is_shine_) { + pdg = rec_track.GetField(q_field_id_) > 0 ? 2212 : -211; + } + else { + pdg = rec_track.GetField(pdg_field_id_); + } + return pdg; +} + +void InConverter::FillTrack(const AnalysisTree::Track& rec_track, InputContainer& input_info) +{ + std::vector mf(kNumberOfFieldPars, 0.f); + for(int iF=0; iF(mf_field_id_+iF); + + auto cov_matrix = is_shine_ ? GetCovMatrixShine(rec_track) : GetCovMatrixCbm(rec_track); + std::vector par(kNumberOfTrackPars,0.f); + + par[kX] = rec_track.GetField(par_field_id_); + par[kY] = rec_track.GetField(par_field_id_ + 1); + par[kZ] = rec_track.GetField(par_field_id_ + 2); + par[kPx] = rec_track.GetPx(); + par[kPy] = rec_track.GetPy(); + par[kPz] = rec_track.GetPz(); + + const int pdg = GetTrackPid(rec_track); + + input_info.AddTrack(par, cov_matrix, mf, rec_track.GetField(q_field_id_), pdg, rec_track.GetId(), rec_track.GetField(passcuts_field_id_)); +} + + +InputContainer InConverter::CreateInputContainer(int iEvent) +{ + in_chain_->GetEntry(iEvent); + const int n_tracks = kf_tracks_->GetNumberOfChannels(); + std::cout << "Event # " << iEvent << " with Ntracks = " << n_tracks << std::endl; + + InputContainer input_info; + input_info.SetCuts(cuts_); + input_info.SetPV(rec_event_header_->GetVertexX(), rec_event_header_->GetVertexY(), rec_event_header_->GetVertexZ()); + + int n_good_tracks{0}; + for(int i_track=0; i_trackGetChannel(i_track); + if(!IsGoodTrack(rec_track)) continue; + FillTrack(rec_track, input_info); + n_good_tracks++; +// ======= +// const AnalysisTree::Track& rec_track = kf_tracks_->GetChannel(i_track); +// +// if (track_cuts_ != nullptr) +// if(!track_cuts_->Apply(rec_track)) +// continue; +// +// std::vector par; +// for(int iP=0; iP<3; iP++) +// par.push_back(rec_track.GetField(par_field_id_ + iP)); +// par.push_back(rec_track.GetPx()); +// par.push_back(rec_track.GetPy()); +// par.push_back(rec_track.GetPz()); +// +// std::vector mf; +// for(int iF=0; iF<10; iF++) +// mf.push_back(rec_track.GetField(mf_field_id_+iF)); +// +// int pdg = -1; +// if(is_shine_) { +// pdg = rec_track.GetField(q_field_id_) > 0 ? 2212 : -211; +// } +// else { +// pdg = rec_track.GetField(pdg_field_id_); +// } +// auto cov_matrix = is_shine_ ? GetCovMatrixShine(rec_track) : GetCovMatrixCbm(rec_track); +// +// // if(ststrack.GetField(m_pdg_field_id_)!=3122) continue; +// inputInfo.AddTrack( par, +// cov_matrix, +// mf, +// rec_track.GetField(q_field_id_), +// pdg, +// rec_track.GetId(), +// rec_track.GetField(passcuts_field_id_)); +// // ststrack.GetField(nhits_field_id_), +// // ststrack.GetField(passcuts_field_id_)); +// >>>>>>> lubynets-order + } + + std::cout << "Good tracks = " << n_good_tracks << std::endl; + return input_info; +} + +bool InConverter::IsGoodTrack(const AnalysisTree::Track& rec_track) +{ + return true; + bool is_good{false}; + + const int sim_id = kf2sim_tracks_->GetMatch(rec_track.GetId()); +// std::cout<< "sim " << sim_id << std::endl; + if (sim_id >= 0 && sim_id < sim_tracks_->GetNumberOfChannels()) { + const AnalysisTree::Track& sim_track = sim_tracks_->GetChannel(sim_id); + const int mother_id = sim_track.GetField(mother_id_field_id_); +// std::cout << "mother " << mother_id << std::endl; + + if (mother_id >= 0 && mother_id < sim_tracks_->GetNumberOfChannels() ) + { + const AnalysisTree::Track& mother_track = sim_tracks_->GetChannel(mother_id); + const int mother_pdg = mother_track.GetField(sim_pdg_field_id_); + std::cout << "mother pdg " << mother_pdg << std::endl; + + if(mother_pdg == 3122) + is_good = true; + } + } + return is_good; +} + +std::vector InConverter::GetCovMatrixCbm(const AnalysisTree::Track& track) const +{ + const double tx = track.GetField(par_field_id_+3); + const double ty = track.GetField(par_field_id_+4); + const double qp = track.GetField(par_field_id_+5); + const Int_t q = track.GetField(q_field_id_); + + //calculate covariance matrix + const double t = sqrt(1.f + tx*tx + ty*ty); + const double t3 = t*t*t; + const double dpxdtx = q/qp*(1.f+ty*ty)/t3; + const double dpxdty = -q/qp*tx*ty/t3; + const double dpxdqp = -q/(qp*qp)*tx/t; + const double dpydtx = -q/qp*tx*ty/t3; + const double dpydty = q/qp*(1.f+tx*tx)/t3; + const double dpydqp = -q/(qp*qp)*ty/t; + const double dpzdtx = -q/qp*tx/t3; + const double dpzdty = -q/qp*ty/t3; + const double dpzdqp = -q/(qp*qp*t3); + + const double F[6][5] = { {1.f, 0.f, 0.f, 0.f, 0.f}, + {0.f, 1.f, 0.f, 0.f, 0.f}, + {0.f, 0.f, 0.f, 0.f, 0.f}, + {0.f, 0.f, dpxdtx, dpxdty, dpxdqp}, + {0.f, 0.f, dpydtx, dpydty, dpydqp}, + {0.f, 0.f, dpzdtx, dpzdty, dpzdqp} }; + + double VFT[5][6]; + for(int i=0; i<5; i++) + for(int j=0; j<6; j++) + { + VFT[i][j] = 0; + for(int k=0; k<5; k++) + { + if (k<=i) + VFT[i][j] += track.GetField(cov_field_id_ + k + i*(i+1)/2) * F[j][k]; //parameters->GetCovariance(i,k) * F[j][k]; + else + VFT[i][j] += track.GetField(cov_field_id_ + i + k*(k+1)/2) * F[j][k]; //parameters->GetCovariance(i,k) * F[j][k]; + } + } + + + std::vector cov(21, 0); + for(int i=0, l=0; i<6; i++) + for(int j=0; j<=i; j++, l++) + { + cov[l] = 0; + for(int k=0; k<5; k++) + { + cov[l] += F[i][k] * VFT[k][j]; + } + } + return cov; +} + +std::vector InConverter::GetCovMatrixShine(const AnalysisTree::Track& track) const +{ + std::vector cov(21, 0.); + + for (int iCov=0; iCov<21; ++iCov) + cov[iCov] = track.GetField(cov_field_id_+iCov); + + return cov; +} diff --git a/AnalysisTreeInterface/old/InConverter.h b/AnalysisTreeInterface/old/InConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..6ceb729344564d92b2bc0f679ca9008f12a301c7 --- /dev/null +++ b/AnalysisTreeInterface/old/InConverter.h @@ -0,0 +1,76 @@ +#ifndef AT_KFC_in_HH +#define AT_KFC_in_HH + +#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" +#include "AnalysisTree/Cuts.h" + +#include "InputContainer.h" + +class InConverter +{ +public: + + // Constructors/Destructors --------- + InConverter() = default; + virtual ~InConverter() = default; + + void InitAnalysisTree(const std::string& file_name, const std::string& tree_name); + InputContainer CreateInputContainer(int iEvent); + + void SetIsShine(bool is=true) { is_shine_ = is; } + int GetNEvents(){return in_chain_ -> GetEntries();}; + + void SetCuts(const CutsContainer& cuts) { cuts_ = cuts; }; + void SetTrackCuts(AnalysisTree::Cuts* const cuts) { track_cuts_ = cuts; }; + +protected: + std::vector GetCovMatrixCbm(const AnalysisTree::Track& track) const; + std::vector GetCovMatrixShine(const AnalysisTree::Track& track) const; + bool IsGoodTrack(const AnalysisTree::Track& rec_track); + void FillTrack(const AnalysisTree::Track& rec_track, InputContainer& input_info); + int GetTrackPid(const AnalysisTree::Track& rec_track); + + TFile* in_file_{nullptr}; + TChain* in_chain_{nullptr}; + + AnalysisTree::TrackDetector* kf_tracks_{nullptr}; + AnalysisTree::TrackDetector* sim_tracks_{nullptr}; + AnalysisTree::Matching* kf2sim_tracks_{nullptr}; + AnalysisTree::EventHeader* rec_event_header_{nullptr}; + AnalysisTree::EventHeader* sim_event_header_{nullptr}; + AnalysisTree::Configuration* config_{nullptr}; +// <<<<<<< HEAD + +// ======= +// + AnalysisTree::Cuts* track_cuts_{nullptr}; +// +// const int nParametersKF{6}; +// >>>>>>> lubynets-order + CutsContainer cuts_; + + // field ids for input parameters + int q_field_id_{AnalysisTree::UndefValueInt}; + int pdg_field_id_{AnalysisTree::UndefValueInt}; + int par_field_id_ {AnalysisTree::UndefValueInt}; + int mf_field_id_ {AnalysisTree::UndefValueInt}; + int cov_field_id_{AnalysisTree::UndefValueInt}; + int mother_id_field_id_{AnalysisTree::UndefValueInt}; + int sim_pdg_field_id_{AnalysisTree::UndefValueInt}; + int passcuts_field_id_{AnalysisTree::UndefValueInt}; + + bool is_shine_{false}; + + +}; + +#endif diff --git a/AnalysisTreeInterface/old/Manager.cxx b/AnalysisTreeInterface/old/Manager.cxx new file mode 100644 index 0000000000000000000000000000000000000000..3a8ea0b3769097fcd98e3bd05c1b9761dd9628f3 --- /dev/null +++ b/AnalysisTreeInterface/old/Manager.cxx @@ -0,0 +1,33 @@ +#include "Manager.h" + +void Manager::Run(int n_events) +{ + InConverter converter_in; + converter_in.SetTrackCuts(track_cuts_); + converter_in.InitAnalysisTree(in_file_name_, in_tree_name_); + converter_in.SetIsShine(is_shine_); + + OutConverter converter_out; + converter_out.SetOutFileName("KFPS_output.root"); + converter_out.InitAT(); + + std::cout << converter_in.GetNEvents() << std::endl; + + if(n_events<0 || n_events>converter_in.GetNEvents()) + n_events = converter_in.GetNEvents(); + + for(int iEvent=0; iEvent +#include + +#include "AnalysisTreeInterface/old/InConverter.h" +#include "AnalysisTreeInterface/old/OutConverter.h" +#include "SimpleFinder.h" + +class Manager +{ +public: + + void Init() {}; + void Run(int n_events=-1); + void Finish(){}; + + void SetInFileName(const std::string& name){in_file_name_ = name;}; + void SetInTreeName(const std::string& name){in_tree_name_ = name;}; + void SetIsShine(bool is=true) { is_shine_ = is; } + + void SetCuts(const CutsContainer& cuts) { cuts_ = cuts; }; + void SetTrackCuts(AnalysisTree::Cuts* const cuts) { track_cuts_ = cuts; }; + + protected: + std::string in_file_name_; + std::string in_tree_name_; + CutsContainer cuts_; + AnalysisTree::Cuts* track_cuts_{nullptr}; + + bool is_shine_{false}; + +}; +#endif \ No newline at end of file diff --git a/AnalysisTreeInterface/old/OutConverter.cxx b/AnalysisTreeInterface/old/OutConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..da1fbfd1a7da790ee9368e21ecc1ee884285c7e1 --- /dev/null +++ b/AnalysisTreeInterface/old/OutConverter.cxx @@ -0,0 +1,116 @@ +#include "OutConverter.h" + +void OutConverter::InitAT() +{ + 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("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"); + + out_config_.AddBranchConfig( LambdaRecoBranch ); + lambda_reco_ = new AnalysisTree::Particles(out_config_.GetLastId()); + + out_tree_ = new TTree("aTree", "AnalysisTree Lambda"); + out_tree_ -> Branch("LambdaCandidates", "AnalysisTree::Particles", &lambda_reco_); + out_config_.Write("Configuration"); + + 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"); + + 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"); +} + +void OutConverter::WriteCandidates(const std::vector& canditates) +{ + lambda_reco_->ClearChannels(); + + for(const auto& candidate : canditates) + { + auto* lambdarec = lambda_reco_->AddChannel(); + lambdarec->Init(out_config_.GetBranchConfig(lambda_reco_->GetId())); + + 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() - RapidityShift(12.), rap_cm_field_id_); + lambdarec->SetField( particle.GetPDG(), pdg_field_id_w_ ); + + lambdarec->SetPid( particle.GetPDG()); + lambdarec->SetMass(mass); + + lambdarec->SetField( candidate.GetChi2PrimPos(), chi2primpos_field_id_); + lambdarec->SetField( candidate.GetChi2PrimNeg(), chi2primneg_field_id_); + lambdarec->SetField( candidate.GetDistance(), distance_field_id_); + lambdarec->SetField( candidate.GetCosineDaughterPos(), cosinepos_field_id_); + lambdarec->SetField( candidate.GetCosineDaughterNeg(), cosineneg_field_id_); + lambdarec->SetField( candidate.GetChi2Geo(), chi2geo_field_id_); + lambdarec->SetField( candidate.GetL(), l_field_id_); + lambdarec->SetField( candidate.GetLdL(), ldl_field_id_); + lambdarec->SetField( candidate.GetIsFromPV(), isfrompv_field_id_); + lambdarec->SetField( candidate.GetCosineTopo(), cosinetopo_field_id_); + lambdarec->SetField( candidate.GetChi2Topo(), chi2topo_field_id_); + } + out_tree_->Fill(); +} + +void OutConverter::Finish() +{ + out_tree_->Write(); + out_file_->Close(); +} + +float OutConverter::RapidityShift(float pbeam) +{ + const float nucleonmass = 0.939; + return TMath::Log(((sqrt(nucleonmass*nucleonmass + pbeam*pbeam) + pbeam) / (sqrt(nucleonmass*nucleonmass + pbeam*pbeam) - pbeam)))/4.; +} \ No newline at end of file diff --git a/AnalysisTreeInterface/old/OutConverter.h b/AnalysisTreeInterface/old/OutConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..a3257c7f4e44aed68005859c8fa79bad4c64b72c --- /dev/null +++ b/AnalysisTreeInterface/old/OutConverter.h @@ -0,0 +1,70 @@ +#ifndef KFC_out_AT_HH +#define KFC_out_AT_HH + +#include +#include + +#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" + +#include "OutputContainer.h" + +class OutConverter +{ +public: + + // Constructors/Destructors --------- + OutConverter() = default; + virtual ~OutConverter() = default; + + void InitAT(); + void SetOutFileName(std::string name){out_file_name_ = std::move(name);}; + void WriteCandidates(const std::vector& canditates); + + void Finish(); + +protected: + + float RapidityShift(float pbeam); + + std::string out_file_name_{"KFPS_output.root"}; + TFile* out_file_{nullptr}; + TTree* out_tree_{nullptr}; + + AnalysisTree::Configuration out_config_; + AnalysisTree::Particles* lambda_reco_{nullptr}; + + // 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}; + + // field ids for lambda candidate cutting parameters + 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}; + +}; +#endif \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 92ccf4c95f8ae7f2fde0a584dc4156e1ae903cd1..5c1c75c0a3b3b7bb42a969e345ed1150702ec1f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,7 @@ message(STATUS "Using C++${CMAKE_CXX_STANDARD}") # by default build optimized code with debug symbols if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - set(CMAKE_BUILD_TYPE RELWITHDEBINFO) + set(CMAKE_BUILD_TYPE RELEASE) endif () # in DEBUG mode make verbose Makefile @@ -30,7 +30,8 @@ endif () # set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -g -DDEBUG -D__DEBUG -Wall") set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -DDEBUG -D__DEBUG -Wall") set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELEASE} -ggdb") -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ftree-vectorize -ffast-math -DNODEBUG") +# set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ftree-vectorize -ffast-math -DNODEBUG") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -ftree-vectorize -ffast-math -DNODEBUG") message(STATUS "Using CXX flags for ${CMAKE_BUILD_TYPE}: ${CMAKE_CXX_FLAGS_${CMAKE_BUILD_TYPE}}") list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS}) @@ -38,32 +39,33 @@ set(FIXTARGET TRUE CACHE BOOL "Compile for fix target geometry.") find_package(ROOT REQUIRED RIO) find_package(KFParticle QUIET) -find_package(Vc QUIET) +find_package(AnalysisTree QUIET) set(EXTERNAL_DIR ${CMAKE_BINARY_DIR}/external) set(EXTERNAL_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/external) -if(Vc_FOUND) - Message("Vc found") -Else() - Message("Vc is not found. It will be installed as external package") - include(cmake_modules/Vc.cmake) -EndIf() - if(KFParticle_FOUND) Message("KFParticle found") + find_package(Vc REQUIRED) + include_directories(${KFParticle_INCLUDE_DIR} ${Vc_INCLUDE_DIR}) Else() Message("KFParticle is not found. It will be installed as external package") include(cmake_modules/KFParticle.cmake) + include_directories(${CMAKE_SOURCE_DIR} ${ROOT_INCLUDE_DIRS} ${PROJECT_INCLUDE_DIRECTORIES}) EndIf() -if (ROOT_FOUND) - message(STATUS "Using ROOT: ${ROOT_VERSION} <${ROOT_CONFIG}>") - include_directories(${CMAKE_SOURCE_DIR} ${ROOT_INCLUDE_DIR} ${ROOT_INCLUDE_DIRS}) - include(${ROOT_USE_FILE}) -endif (ROOT_FOUND) +#message(STATUS "PROJECT_INCLUDE_DIRECTORIES ${PROJECT_INCLUDE_DIRECTORIES}") -include_directories(${CMAKE_SOURCE_DIR} ${KFParticle_INCLUDE_DIR} ${Vc_INCLUDE_DIR}) +message(STATUS "Using ROOT: ${ROOT_VERSION} <${ROOT_CONFIG}>") +include_directories(${CMAKE_SOURCE_DIR} ${ROOT_INCLUDE_DIRS}) +include(${ROOT_USE_FILE}) + +if(AnalysisTree_FOUND) + Message("AnalysisTree found. Corresponding interface will be compiled") + add_subdirectory(AnalysisTreeInterface) +Else() + Message("AnalysisTree is not found. Corresponding interface will not be compiled") +EndIf() add_subdirectory(KFSimple) add_subdirectory(Interface) diff --git a/Interface/CMakeLists.txt b/Interface/CMakeLists.txt index d00b47ef3c12272da2894dcdc607f6aab305454b..21831abb8efceb13fda889c8fb7eaa0514707b25 100644 --- a/Interface/CMakeLists.txt +++ b/Interface/CMakeLists.txt @@ -6,14 +6,18 @@ set(SOURCES string(REPLACE ".cxx" ".h" HEADERS "${SOURCES}") include_directories(${CMAKE_SOURCE_DIR}/Interface ) + +include_directories(${PROJECT_INCLUDE_DIRECTORIES}) link_directories(${PROJECT_LINK_DIRECTORIES}) add_library(KFParticleInterface SHARED ${SOURCES} G__KFParticleInterface.cxx) -if(ROOT_FOUND) - ROOT_GENERATE_DICTIONARY(G__KFParticleInterface ${HEADERS} LINKDEF KFParticleInterfaceLinkDef.h OPTIONS "-DDO_TPCCATRACKER_EFF_PERFORMANCE" "-DNonhomogeneousField" "-DCBM" "-DUSE_TIMERS") - target_link_libraries(KFParticleInterface ${ROOT_LIBRARIES} KFParticle ${Vc_LIBRARIES}) -endif(ROOT_FOUND) +if(PROJECT_DEPENDENCIES) + add_dependencies(KFParticleInterface ${PROJECT_DEPENDENCIES}) +endif() + +ROOT_GENERATE_DICTIONARY(G__KFParticleInterface ${HEADERS} LINKDEF KFParticleInterfaceLinkDef.h OPTIONS "-DDO_TPCCATRACKER_EFF_PERFORMANCE" "-DNonhomogeneousField" "-DCBM" "-DUSE_TIMERS") +target_link_libraries(KFParticleInterface ${ROOT_LIBRARIES} KFParticle ${Vc_LIBRARIES}) add_target_property(KFParticleInterface COMPILE_FLAGS "-DDO_TPCCATRACKER_EFF_PERFORMANCE -DNonhomogeneousField -DCBM -DUSE_TIMERS") diff --git a/Interface/InputContainer.cxx b/Interface/InputContainer.cxx index aff0079c3d9537d292f2a56b849a1c215eeff312..7d77f6ee3ec183c25dc2a1df6478d537bd75ba01 100644 --- a/Interface/InputContainer.cxx +++ b/Interface/InputContainer.cxx @@ -14,7 +14,7 @@ void InputContainer::SetPV(float x, float y, float z) primVtx_tmp.SetNContributors( 0 ); primVtx_tmp.SetChi2( -100 ); - vtx_ = KFVertex(primVtx_tmp); + vtx_ = KFVertex(primVtx_tmp); } void InputContainer::AddTrack(const std::vector& par, @@ -22,15 +22,13 @@ void InputContainer::AddTrack(const std::vector& par, const std::vector& field, int charge, int pdg, - int id, - int passcuts) + int id) { - if (passcuts==0 || pdg==0 || pdg==-2) + if (pdg==0 || pdg==-2) return; - + if( par.size() != kNumberOfTrackPars || cov.size() != NumberOfCovElements || field.size() != kNumberOfFieldPars){ - std::cout << "InputContainer::AddTrack - Wrong size of input vector!!" << std::endl; - exit(kError); + throw std::runtime_error("Wrong size of input vector!"); } KFParticle particle; @@ -50,8 +48,39 @@ void InputContainer::AddTrack(const std::vector& par, particle.Q() = char(charge); //NOTE: is not safe particle.SetPDG(pdg); particle.SetId(id); + + tracks_.push_back(particle); +} + +double InputContainer::InversedChi2Prob(double p, int ndf) +{ + const double epsilon = 1.e-14; + double chi2Left = 0.f; + double chi2Right = 10000.f; + + double probLeft = p - TMath::Prob(chi2Left, ndf); - tracks_.emplace_back(particle); + double chi2Centr = (chi2Left+chi2Right)/2.f; + double probCentr = p - TMath::Prob( chi2Centr, ndf); + + while( TMath::Abs(chi2Right-chi2Centr)/chi2Centr > epsilon ) { + if(probCentr * probLeft > 0.f) { + chi2Left = chi2Centr; + probLeft = probCentr; + } + else { + chi2Right = chi2Centr; + } + + chi2Centr = (chi2Left+chi2Right)/2.f; + probCentr = p - TMath::Prob( chi2Centr, ndf); + } + + return chi2Centr; +} + +void InputContainer::Reserve(size_t n) { + tracks_.reserve(n); } // KFParticleTopoReconstructor* InputContainer::CreateTopoReconstructor() @@ -59,17 +88,17 @@ void InputContainer::AddTrack(const std::vector& par, // /* // * Creates the pointer on the KFParticleTopoReconstructor object // * with all necessary input information in order to perform particle selection using -// * non-simplified "standard" KFParticle algorithm. +// * non-simplified "standard" KFParticle algorithm. // */ // auto* TR = new KFParticleTopoReconstructor; -// +// // // cuts setting // TR -> GetKFParticleFinder() -> SetChiPrimaryCut2D(cuts_.GetCutChi2PrimPos()); // TR -> GetKFParticleFinder() -> SetMaxDistanceBetweenParticlesCut(cuts_.GetCutDistance()); // TR -> GetKFParticleFinder() -> SetChi2Cut2D(cuts_.GetCutChi2Geo()); // TR -> GetKFParticleFinder() -> SetLCut(cuts_.GetCutLDown()); // TR -> GetKFParticleFinder() -> SetLdLCut2D(cuts_.GetCutLdL()); -// +// // KFPTrackVector track_tmp, track_empty; // track_tmp.Resize(tracks_.size()); // for(int iTr=0; iTr& par, // track_tmp.SetFieldCoefficient(tracks_[iTr].GetFieldCoeff()[iF], iF, iTr); // track_tmp.SetPDG(tracks_[iTr].GetPDG(), iTr); // track_tmp.SetQ(tracks_[iTr].GetQ(), iTr); -// track_tmp.SetPVIndex(-1, iTr); +// track_tmp.SetPVIndex(-1, iTr); // track_tmp.SetId(tracks_[iTr].Id(), iTr); // } // TR->Init(track_tmp, track_empty); -// +// // TR->AddPV(vtx_); -// +// // return TR; // } @@ -98,12 +127,12 @@ void InputContainer::AddTrack(const std::vector& par, // * Creates the SimpleFinder object with all necessary input information in order to // * perform particle selection using KFSimple algorithm. // */ -// +// // SimpleFinder FCF; -// +// // KFPTrackVector track_tmp; // track_tmp.Resize(tracks_.size()); -// +// // for(int iTr=0; iTr& par, // track_tmp.SetFieldCoefficient(tracks_[iTr].GetFieldCoeff()[iF], iF, iTr); // track_tmp.SetPDG(tracks_[iTr].GetPDG(), iTr); // track_tmp.SetQ(tracks_[iTr].GetQ(), iTr); -// track_tmp.SetPVIndex(-1, iTr); +// track_tmp.SetPVIndex(-1, iTr); // track_tmp.SetId(tracks_[iTr].Id(), iTr); // } // FCF.Init(track_tmp, vtx_); // FCF.SetCuts(cuts_); -// -// return FCF; +// +// return FCF; // } - -double InputContainer::InversedChi2Prob(double p, int ndf) const -{ - const double epsilon = 1.e-14; - double chi2Left = 0.f; - double chi2Right = 10000.f; - - double probLeft = p - TMath::Prob(chi2Left, ndf); - - double chi2Centr = (chi2Left+chi2Right)/2.f; - double probCentr = p - TMath::Prob( chi2Centr, ndf); - - while( TMath::Abs(chi2Right-chi2Centr)/chi2Centr > epsilon ) - { - if(probCentr * probLeft > 0.f) - { - chi2Left = chi2Centr; - probLeft = probCentr; - } - else - { - chi2Right = chi2Centr; - } - - chi2Centr = (chi2Left+chi2Right)/2.f; - probCentr = p - TMath::Prob( chi2Centr, ndf); - } - - return chi2Centr; -} diff --git a/Interface/InputContainer.h b/Interface/InputContainer.h index 5e053ad4f85d1522dd881542c42246dc59d3f0f5..97fbd387a430e44c8d0b5eef1505db8e4b4591ec 100644 --- a/Interface/InputContainer.h +++ b/Interface/InputContainer.h @@ -42,23 +42,28 @@ class InputContainer{ void SetPV(float x, float y, float z); void SetPV(KFVertex vertex); void SetPV(KFPVertex vertex); - void AddTrack(const std::vector& par, const std::vector& cov, const std::vector& field, int charge, int pdg, int id, int passcuts=1); + void AddTrack(const std::vector& par, const std::vector& cov, const std::vector& field, int charge, int pdg, int id); // KFParticleTopoReconstructor* CreateTopoReconstructor(); //^ not good - void SetCuts(const CutsContainer& cuts) { cuts_ = cuts; }; const KFVertex& GetVertex() const {return vtx_;}; const std::vector& GetTracks() const {return tracks_;}; const CutsContainer& GetCuts() const {return cuts_;}; - + void Clear() { + if(!tracks_.empty()){ + tracks_.clear(); + } + } + + void Reserve(size_t n); protected: - double InversedChi2Prob(double p, int ndf) const; + static double InversedChi2Prob(double p, int ndf) ; KFVertex vtx_; - std::vector tracks_; + std::vector tracks_{}; CutsContainer cuts_; }; diff --git a/KFSimple/CMakeLists.txt b/KFSimple/CMakeLists.txt index f89a2450fa2542b286ce6c93b77a5ee1b368b43f..1b944f583aeebc9388fd0d9910cdbe8933ff2f9c 100644 --- a/KFSimple/CMakeLists.txt +++ b/KFSimple/CMakeLists.txt @@ -10,11 +10,9 @@ link_directories(${PROJECT_LINK_DIRECTORIES}) add_library(KFParticleSimple SHARED ${SOURCES} G__KFParticleSimple.cxx) -if(ROOT_FOUND) - ROOT_GENERATE_DICTIONARY(G__KFParticleSimple ${HEADERS} LINKDEF KFSimpleLinkDef.h OPTIONS "-DDO_TPCCATRACKER_EFF_PERFORMANCE" "-DNonhomogeneousField" "-DCBM" "-DUSE_TIMERS") - add_dependencies(KFParticleSimple KFParticleInterface) - target_link_libraries(KFParticleSimple ${ROOT_LIBRARIES} KFParticle ${Vc_LIBRARIES} KFParticleInterface) -endif(ROOT_FOUND) +ROOT_GENERATE_DICTIONARY(G__KFParticleSimple ${HEADERS} LINKDEF KFSimpleLinkDef.h OPTIONS "-DDO_TPCCATRACKER_EFF_PERFORMANCE" "-DNonhomogeneousField" "-DCBM" "-DUSE_TIMERS") +add_dependencies(KFParticleSimple KFParticleInterface) +target_link_libraries(KFParticleSimple ${ROOT_LIBRARIES} KFParticle KFParticleInterface) add_target_property(KFParticleSimple COMPILE_FLAGS "-DDO_TPCCATRACKER_EFF_PERFORMANCE -DNonhomogeneousField -DCBM -DUSE_TIMERS") diff --git a/KFSimple/SimpleFinder.cxx b/KFSimple/SimpleFinder.cxx index 5083800e0b8957405dd68f631a49862b2b0c730f..901871a579dd0cd0831e135cb2e772240f2af997 100644 --- a/KFSimple/SimpleFinder.cxx +++ b/KFSimple/SimpleFinder.cxx @@ -14,7 +14,7 @@ void SimpleFinder::Init(const InputContainer& input) KFPTrackVector track_tmp; const std::vector& tracks = input.GetTracks(); track_tmp.Resize(tracks.size()); - + for(int iTr=0; iTr const float dx = pars1.at(0) - pars2.at(0); const float dy = pars1.at(1) - pars2.at(1); const float dz = pars1.at(2) - pars2.at(2); - const float dr = sqrt(dx*dx+dy*dy+dz*dz); + const float dr = std::sqrt(dx*dx+dy*dy+dz*dz); return dr; } @@ -132,7 +132,7 @@ float SimpleFinder::CalculateCosMomentumSum(const std::array &pars1, c const std::array PSum = {P1.at(0)+P2.at(0), P1.at(1)+P2.at(1), P1.at(2)+P2.at(2)}; return (P1.at(0)*PSum.at(0) + P1.at(1)*PSum.at(1) + P1.at(2)*PSum.at(2)) / - (sqrt(P1.at(0)*P1.at(0) + P1.at(1)*P1.at(1) + P1.at(2)*P1.at(2)) * sqrt(PSum.at(0)*PSum.at(0) + PSum.at(1)*PSum.at(1) + PSum.at(2)*PSum.at(2))); + (sqrt(P1.at(0)*P1.at(0) + P1.at(1)*P1.at(1) + P1.at(2)*P1.at(2)) * std::sqrt(PSum.at(0)*PSum.at(0) + PSum.at(1)*PSum.at(1) + PSum.at(2)*PSum.at(2))); } KFParticleSIMD SimpleFinder::ConstructMother(const KFPTrack &track1, const int pid1, const KFPTrack &track2, const int pid2) const @@ -202,7 +202,7 @@ float SimpleFinder::CalculateCosTopo(const KFParticleSIMD& mother) const const float delta_z = z_mother - prim_vx_.GetZ(); const float sp = delta_x*px_mother + delta_y*py_mother + delta_z*pz_mother; - const float norm = sqrt(delta_x*delta_x + delta_y*delta_y + delta_z*delta_z) * sqrt(px_mother*px_mother + py_mother*py_mother + pz_mother*pz_mother); + const float norm = std::sqrt(delta_x*delta_x + delta_y*delta_y + delta_z*delta_z) * std::sqrt(px_mother*px_mother + py_mother*py_mother + pz_mother*pz_mother); return sp/norm; } @@ -291,7 +291,7 @@ void SimpleFinder::FindParticles() if((cuts_.GetIsApplyCutLDown())&&(lambda.GetL() <= cuts_.GetCutLDown())) continue; lambda.SetChi2Topo(CalculateChi2Topo(mother)); - if((cuts_.GetIsApplyCutChi2Topo())&&(lambda.GetChi2Topo() > cuts_.GetCutChi2Topo()) || lambda.GetChi2Topo()!=lambda.GetChi2Topo()) continue; + if((cuts_.GetIsApplyCutChi2Topo() && lambda.GetChi2Topo() > cuts_.GetCutChi2Topo()) || lambda.GetChi2Topo()!=lambda.GetChi2Topo()) continue; KFParticle particle; mother.GetKFParticle(particle, 0); diff --git a/cmake_modules/KFParticle.cmake b/cmake_modules/KFParticle.cmake index ae4dc04ace824d423fca5938c2c04500d1c012f9..b3039c30898e384a98689c4224748ea9f392a5ad 100644 --- a/cmake_modules/KFParticle.cmake +++ b/cmake_modules/KFParticle.cmake @@ -10,6 +10,7 @@ ExternalProject_Add(KFParticle_Ext BINARY_DIR "${EXTERNAL_DIR}/KFParticle_build" INSTALL_DIR "${KFParticle_INSTALL_DIR}" CMAKE_ARGS + "-DEXTERNAL_INSTALL_DIR=${KFParticle_INSTALL_DIR}" "-DCMAKE_INSTALL_PREFIX=${KFParticle_INSTALL_DIR}" "-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}" "-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}" diff --git a/cmake_modules/Vc.cmake b/cmake_modules/Vc.cmake deleted file mode 100644 index 0e74b7002b0e72717e72842d03604e062f3cd12f..0000000000000000000000000000000000000000 --- a/cmake_modules/Vc.cmake +++ /dev/null @@ -1,21 +0,0 @@ -set(Vc_INSTALL_DIR ${EXTERNAL_INSTALL_DIR}) -set(Vc_INCLUDE_DIR ${Vc_INSTALL_DIR}/include) -set(Vc_LIBRARY_DIR ${Vc_INSTALL_DIR}/lib) - -ExternalProject_Add(Vc_Ext - GIT_REPOSITORY "https://github.com/VcDevel/Vc" - GIT_TAG "1.4.1" - UPDATE_DISCONNECTED ${UPDATE_DISCONNECTED} - SOURCE_DIR "${EXTERNAL_DIR}/Vc_src" - BINARY_DIR "${EXTERNAL_DIR}/Vc_build" - INSTALL_DIR "${Vc_INSTALL_DIR}" - CMAKE_ARGS - "-DCMAKE_INSTALL_PREFIX=${Vc_INSTALL_DIR}" - "-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}" - "-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}" - "-DBUILD_TESTING=OFF" - ) - -list(APPEND PROJECT_DEPENDENCIES Vc_Ext) -list(APPEND PROJECT_LINK_DIRECTORIES ${Vc_LIBRARY_DIR}) -list(APPEND PROJECT_INCLUDE_DIRECTORIES ${Vc_INCLUDE_DIR}) \ No newline at end of file