diff --git a/.github/workflows/clang.yml b/.github/workflows/clang.yml index 7d38565..183cc70 100644 --- a/.github/workflows/clang.yml +++ b/.github/workflows/clang.yml @@ -11,5 +11,5 @@ jobs: - uses: actions/checkout@v2 - uses: RafikFarhad/clang-format-github-action@v2.0.0 with: - sources: "src/*.cpp src/*.h pid_new/*/*.cpp pid_new/*/*.h macro/*.C tasks/*.cpp input/*.C at_interface/*pp" + sources: "src_tof/*.cpp src_tof/*.h src_trd/*.cpp src_trd/*.h pid_new/*/*.cpp pid_new/*/*.h macro/*.C tasks/*.cpp input/*.C at_interface/*pp interface/*pp" style: file diff --git a/.gitignore b/.gitignore index 53f726d..44ccbdd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ Doxygen/* *build*/* *install*/* +.DS_Store +*/.DS_Store diff --git a/CMakeLists.txt b/CMakeLists.txt index 588e82b..6c9fff1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,13 +47,15 @@ set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -g -DDEBUG -D__DEBUG -Wall -Wextra") 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}}") -add_subdirectory(src) +add_subdirectory(src_tof) +add_subdirectory(src_trd) +add_subdirectory(interface) if(${PID_BUNDLED_AT}) add_subdirectory(at_interface) endif() add_subdirectory(tasks) - # Packaging routine +# Packaging routine # For complete explanation # @see https://cmake.org/cmake/help/git-master/manual/cmake-packages.7.html#creating-packages include(GenerateExportHeader) diff --git a/README.md b/README.md index bbd7866..2c4a5bf 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,15 @@ # Particle Identification Framework ## General information -The procedure of particle identification is described in Sec. 4.4 of this [thesis](https://publikationen.ub.uni-frankfurt.de/opus4/frontdoor/deliver/index/docId/51779/file/main.pdf). +The Pid framework identifies hadrons based on the Tof-*m*2 and/or the Trd-*dE/dx*. It assigns a probability to be a certain particle species to +every track depending on its *m*2 and its momentum and/or on its Trd-dEdx and its momentum. + +The procedure of particle identification with the Tof-*m*2 is described in Sec. 4.4 of this [thesis](https://publikationen.ub.uni-frankfurt.de/opus4/frontdoor/deliver/index/docId/51779/file/main.pdf). It consists of two steps - determination of particle type hypothesis (fitting) and application of hypothesis to a set of particles (filling). -Running of these steps is described below, in the section "Examples". + +The procedure of particle identification with the the Trd-*dE/dx* is described in [presentation](https://next.hessenbox.de/index.php/s/DJrpKEaH7zZzNDx). It constists of three steps - creating of MC-input, calculating MC-probabilities and assigning probabilities and pid hypothesis to a set of particles (filling). In addition, the separtion of electrons and pions with the RICH detector is applied. + +Running of these steps is described below, in section "First Run". ## Installation ROOT6 is needed for installation. @@ -37,18 +43,89 @@ Now perform following commands for the installation: *path-to-root-installation* and *path-to-analysistree-installation* must be replaced with your actual location of Root and AnalysisTree install directories respectively. -## Examples -Before running the Pid framework you need to set up the environment variables: +## First run + +Before running the Pid framework the environment variables need to be set: source path-to-pid-installation/bin/PidConfig.sh -### Fitting -Let us consider the fitting procedure with an example of 5M Au+Au collisions generated with DCM-QGSM-SMM at beam momentum 12*A* GeV/*c*. -For this you need to prepare an input file with *m*2-*p* histograms, like input/m2_pq_vtx.apr20.dcmqgsm.12agev.root. -In the folder "reco_info" there is a histogram with all reconstructed tracks matched to TOF hits. -In folders "reco_vs_sim_info_from_tof", "reco_vs_sim_info_via_vtx" there are histograms for separate particles species (MC-true) - determined from "TOF hit - simulated particles" and "TOF hit - reconstructed track - simulated particles" matchings respectively.\ +The following steps need to be performed (for more details see below): + +- Preparation: creating mc-histograms + + -- for Tof pid: preparation outside of the Pid-framework + + -- for Trd pid run: + + ./create_mcinput_trd filelist.sh outfilename (for one analysistree.root) + + and merge output-files for runs into one file: + + hadd mcfilename.root *.mcfilename.root + +- Step 1: producing probabilities that are required to identify the tracks. + + -- for Tof pid: + + ./run_pid_dcm12 + + -- for Trd pid: + + ./calculate_mcprob_trd mcfilename + +- Step 2: assigning probabilities for particle species and a pid-hypothesis to reconstructed tracks from an analysistree and writing them into a new analysistree. + + -- for Tof and Trd pid simoultanously: + + ./fill_pid 0 filelist.txt outputfilename pid_file_tof pid_file_trd truncation_mode probability_mode min_hits (purity) + + -- for Tof pid only: + + ./fill_pid 1 filelist.txt outputfilename pid_file_tof + + -- for Trd pid only: + + ./fill_pid 2 filelist.txt outputfilename pid_file_trd truncation_mode probability_mode min_hits (purity) + +The preparation and the first step only need to be performed once. + +In the following the above steps are described in more detail. + +### Preparation: Creating MC-Input + +#### TOF +The MC-input file for a specific system and energy needs to be created by the user before running the Pid framework. The input file needs to contain *m*2-*p* histograms in the format like input/m2_pq_vtx.apr20.dcmqgsm.12agev.root. In the folder "reco_info" there is a histogram with all reconstructed tracks matched to TOF hits. In folders "reco_vs_sim_info_from_tof", "reco_vs_sim_info_via_vtx" there are histograms for separate particles species (MC-true) - determined from "TOF hit - simulated particles" and "TOF hit - reconstructed track - simulated particles" matchings respectively.\ + Another file which you need contains graphic cuts of certain particles species - it allows to ignore during fitting the particles entries in a non-specific parts of the *m*2-*p* histogram, which mainly originate from mismatching between track and TOF hit. This file is input/cuts.root.\ + +#### TRD +The MC-input file for a specific system and energy needs to be created by running tasks/create_mcinput_trd.cpp. The MC-input file contains *dE/dx*-*p* histograms for all reconstructed tracks and hits as well as for reconstructed tracks and hits separated by particle species according to the MC-matching information. Separate histograms are createdf for: + +- tracks and hits +- number of Trd-hits (1-4) +- truncation (1-4) (including average) +- all tracks/hits and separated by particles species (MC-matching information) + +An example file for 5 million minbias Au+Au events at *p*beam = 12 AGeV reconstructed with cbmroot jul25 can be found in input/dEdx_p.jul25.phqmd.auau.12agev.root. + +To run the executable type: + + ./create_mcinput_trd filelist.sh outfilename + +The input for the creation of the MC-histrograms must be an analysistree including mc-information in the format jobId.analysistree.root. + +The path and name of the analysistree is read in via filelist.txt. Only one analysistree is read in one run. A separate filelist.sh is required for every run. + +The ouput filename will be in the format jobId.mcfilename.root + +For reliable MC-histograms a set of run needs to be performed and the resulting files for every run need to get merged into one file with: + + hadd mcfilename.root *.mcfilename.root + +### Calculating probabilities + +#### TOF: Fitting The task which runs the fitting process is tasks/run_pid_dcm12.cpp (look for comments in this file for understanding the configuration example), and it produces a C++ executable install/bin/run_pid_dcm12. To run this executable just type: @@ -66,32 +143,130 @@ Graphs pionpos_0, pionpos_1 and pionpos_2, kaonpos_0, ..., bgpos_1, bgpos_2 repr Finally, the chi2 graph shows the chi2 dependence on momentum when *m*2 fits were performed.\ Running the macro/drawFits.C macro on the allpos.root (allneg.root) file allows to build all these distributions on the same canvas. -Another file produced during the fitting procedure is pid_getter.root. -There is a Pid::Getter which contains results of fitting. +Another file produced during the fitting procedure is pidtof_getter.root. +There is a Pid::GetterTof which contains results of fitting. In order to determine bayesian probabilities (purities) of particle species in a certain *m*2-*p* point (*p*=3 GeV/*c*, *m*2=1 (GeV/*c*2)2), one runs - pid_getter->GetBayesianProbability(3, 1) + pid_getter_tof->GetBayesianProbability(3, 1) and obtains result in the following form: (std::map) { -321 => 0.0000000, -211 => 0.0000000, -1 => 0.0000000, 1 => 0.0011928893, 211 => 8.0296922e-81, 321 => 7.1999527e-46, 2212 => 0.99880711 } // Here 1 and -1 stand for background in the right and left side of the histogram respectively. -Running the macro/RunGetter.C on the pid_getter.root file produces a set of plots showing purity distribution of different particles species and background. +Running the macro/RunGetterTof.C on the pid_getter_tof.root file produces a set of plots showing purity distribution of different particles species and background. + +#### TRD + +The task which runs the calculation of mc-probabilities is tasks/calculate_prob_trd.cpp. To run this executable just type: + + ./calculate_mcprob_trd mcfilename + +Mc-probabilities for all particle species are calculated for every *dE/dx*-*p* for all histograms with two methods: Total probability and likelyhood method (see PID_with_Trd_dEdx.pdf for more details). The probabilities are saved in pid_getter_trd.root. The getter is required for the filling procedure. The probabilites for all histograms in mcfilename.root are also saved in *mcfilename*_probabilities.root. + +**** Additional information on how to work the pid_getter_trd.root independently of the filling procedure: **** + +The Pid::GetterTrd contains the results of the probability calculations. For a reconstructed track with certain *dE/dx* for each Trd-hit and *p* value, probabilities can get obtained for the total proability method with: + + pid_getter_trd->GetTrdProbabilities(trdtrack); + +and for the likelihood method with: + + pid_getter_trd->GetTrdProbabilitiesMulti(trdtrack); + +The trdtrack is an object of the folllowing form: + + TrdContainer trdtrack(mom, pT, charge, nhits_trd, dEdx_hits); + +Options for the application of probabilities and pid hypothesis (see below for a more detailed explanation of options) can be set with: + + pid_getter_trd->SetProbabilityMode(mode); // default is total probability method + pid_getter_trd->SetTruncationMode(mode); // default is average + pid_getter_trd->SetMinHits(minhits); // default is 1 hit + pid_getter_trd->SetPurity(minpurity); // default is 0.0 + +For example, for a postive track with an average *dE/dx* = 25 keV/cm and *p* = 5 GeV/*c*, the obtained result with the total probability method is in the form: + + (std::map) {(2212, 0.5803), (211, 0.2162), (321, 0.0366), (1000010020, 0.0281), (1000010030, 0.0011), (1000020030, 0.0087), (1000020040 = 0.0009, (-11; 0.1277, (-13, 0.0004), (1, 0.0001) } + // Here 1 stands for background + +Alternatively, probabilites can get applied for the average mode and for the total probability method direcly with: + + getter->GetTrdProbabilitiesTotalAverage(mom, dEdx, charge, nhits_trd); + // with dEdx = average dE/dx of a track + +and for the likelihood method with: + + getter->GetTrdProbabilitiesLikelihoodAverage(mom, dEdx_hits, charge); + // with dEdx_hits = dE/dx vector for hits of a track + +The macro/RunGetterTrd.C shows some examples on how to access the pid_getter_trd.root to assign probabilities to a track. I also produces a set of plots showing probability distributions of different particles species and background. ### Filling -Once fitting is preformed and the pid_getter.root file is produced, filling the root-file containing reconstructed tracks can be done - each track will be assigned with probabilities of belonging to different particle species and (optionally) its particle type hypothesis. -This is done in the at_interface/PidFiller.*pp, which are managed by the task tasks/fill_pid.cpp (look for comments in this file for understanding the configuration example). -It produces an executable fill_pid which can be run: +Once fitting / probability calculatoin is preformed and the pid_getter_tof.root and/or pid_getter_trd.root files are produced, filling the root-file containing reconstructed tracks can be done - each track will be assigned with probabilities of belonging to different particle species and (optionally) its particle type hypothesis according to the Tof and/or Trd measurement. + +This is done in the at_interface/PidFiller.*pp, which are managed by the task tasks/fill_pid.cpp. It produces an executable fill_pid which can be run for either Tof and Trd pid simoultanously or for both separately. + +The first argument manages the selection of the detector: =0: Tof and Trd, =1: Tof, =2: Trd. + +The following arguments depend on the selection of the detector. + +To run Tof and Trd pid simoultanously type: + + ./fill_pid 0 filelist.txt outputfilename pid_file_tof pid_file_trd truncation_mode probability_mode min_hits (purity) + +To run Tof pid only type: - ./fill_pid filelist.list pid_getter.root + ./fill_pid 1 filelist.txt outputfilename pid_file_tof + +To run Trd pid only type + + ./fill_pid 2 filelist.txt outputfilename pid_file_trd truncation_mode probability_mode min_hits (purity) + +- filelist.list is a text file containing names of the AnalysisTree root-files to be worked on (an example of AnalysisTree file can be downloaded [here](https://sf.gsi.de/f/3ba5a9e3ff5248edba2c/?dl=1)). + +- pid_file_tof / pid_file_tof are the getter-files that were produced in the step before + +- for the Trd pid the following options need to get selected: + + -- truncation_mode: modes for calculation of energy loss for up to 4 +layers: + + =0: average over all hits (default) + + =1-4: Select hits with lowest dEdx: =1: 1 hit, =2: 2 hits, =3: 3 hits, =4: 4 hits + + -- probability_mode: + + =0: total probability - probability based on particle multiplicites (default) + + =1: likelihood - probability based on dEdx-distribution of particle +species + + -- min_hits: minimum number of required hits per track + + optional: + + -- purity: minimum purity for the pid-hypothesis (most probable particle species) (default purity is 0) -where filelist.list is a text file containing names of the AnalysisTree root-files to be worked on (en example of AnalysisTree file can be downloaded [here](https://sf.gsi.de/f/3ba5a9e3ff5248edba2c/?dl=1)). After running this exacutable a pid.analysistree.root file is produced. -It is based on the input file with following changes: -branches RichRings and TrdTracks are missing, and branch VtxTracks is replaced with RecParticles branch. -The RecParticles differs from VtxTracks in, firstly, the type of branch (Particles instead of Track), and secondly - in few additional fields: prob_K, prob_d, prob_p, prob_pi, prob_bg - a probability of the particle to be a kaon, deuteron, proton, pion or background (undefined type). -Also there is a field 'pid' which shows the most probable particle type hypothesis. +It is based on the input file. The branch VtxTracks is replaced with RecParticles branch. +The RecParticles differs from VtxTracks in, firstly, the type of branch (Particles instead of Track), and secondly - in few additional fields: +- for Tof: + + prob_K, prob_d, prob_p, prob_pi, prob_bg - probability of the particle to be a kaon, deuteron, proton, pion or background (undefined type). + + pid - Tof pid hypothesis which is the most probable particle type + +- for Trd: + + prob_trd_p, prob_trd_pi, prob_trd_K, prob_trd_d, prob_trd_t, prob_trd_He3, prob_trd_He4, prob_trd_e, prob_trd_mu, prob_trd_bg - probability of the particle to be a proton, pion, kaon, deuteron, triton, He3, He4, electron, myon or background (undefined type). + + pid_trd - Trd pid hypothesis which is the most probable particle type + + In addition, when running with Trd, the RICH electron hyposthesis will be added: + + electron_rich - RICH electron hypothesis: true = electron, false = no electron ### Doxygen documentation doxygen docs/Doxyfile diff --git a/at_interface/CMakeLists.txt b/at_interface/CMakeLists.txt index 032ecdc..8918c9c 100644 --- a/at_interface/CMakeLists.txt +++ b/at_interface/CMakeLists.txt @@ -2,16 +2,17 @@ include(AnalysisTree) set(SOURCES PidFiller.cpp + PidTrdMcInput.cpp ) string(REPLACE ".cpp" ".hpp" HEADERS "${SOURCES}") -include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src) +include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src_tof ${CMAKE_SOURCE_DIR}/src_trd) add_library(PidAT SHARED ${SOURCES} G__PidAT.cxx) message("${PROJECT_INCLUDE_DIRECTORIES}") ROOT_GENERATE_DICTIONARY(G__PidAT ${HEADERS} LINKDEF PidATLinkDef.h) -target_link_libraries(PidAT ${ROOT_LIBRARIES} AnalysisTreeBase AnalysisTreeInfra Pid) +target_link_libraries(PidAT ${ROOT_LIBRARIES} AnalysisTreeBase AnalysisTreeInfra Pid PidTrd) install(TARGETS PidAT EXPORT AnalysisTreeFillerTargets LIBRARY DESTINATION lib diff --git a/at_interface/PidFiller.cpp b/at_interface/PidFiller.cpp index 9233f71..6a583ed 100644 --- a/at_interface/PidFiller.cpp +++ b/at_interface/PidFiller.cpp @@ -1,17 +1,44 @@ #include "PidFiller.hpp" +using std::string; using namespace AnalysisTree; -PidFiller::PidFiller(const std::string& file, const std::string& getter) { +PidFiller::PidFiller(const std::string& pid_file_name_tof, const std::string& pid_file_name_trd, const std::string& getter_name_tof, const std::string& getter_name_trd, int pid_mode) { std::cout << "PidFiller::PidFiller" << std::endl; - auto pid_file = std::unique_ptr(TFile::Open(file.c_str(), "read")); - if ((!pid_file) || (pid_file->IsZombie())) { - throw std::runtime_error("No file or file is zombie: " + file); + std::unique_ptr pid_file_tof; + std::unique_ptr pid_file_trd; + + if (pid_mode == 0) { + is_run_pid_tof_ = true; + is_run_pid_trd_ = true; } - getter_ = pid_file->Get(getter.c_str()); - if (getter_ == nullptr) { - throw std::runtime_error("PID Getter is nullptr: " + getter); + + if (pid_mode == 1) + is_run_pid_tof_ = true; + + if (pid_mode == 2) + is_run_pid_trd_ = true; + if (is_run_pid_tof_ == true) { + pid_file_tof = std::unique_ptr(TFile::Open(pid_file_name_tof.c_str(), "read")); + if ((!pid_file_tof) || (pid_file_tof->IsZombie())) { + throw std::runtime_error("No file or file is zombie: " + pid_file_name_tof); + } + getter_tof_ = pid_file_tof->Get(getter_name_tof.c_str()); + if (getter_tof_ == nullptr) { + throw std::runtime_error("PID Tof Getter is nullptr: " + getter_name_tof); + } + } + + if (is_run_pid_trd_ == true) { + pid_file_trd = std::unique_ptr(TFile::Open(pid_file_name_trd.c_str(), "read")); + if ((!pid_file_trd) || (pid_file_trd->IsZombie())) { + throw std::runtime_error("No file or file is zombie: " + pid_file_name_trd); + } + getter_trd_ = pid_file_trd->Get(getter_name_trd.c_str()); + if (getter_trd_ == nullptr) { + throw std::runtime_error("PID Trd Getter is nullptr: " + getter_name_trd); + } } } @@ -24,23 +51,128 @@ int PidFiller::signum(int x) const { return -1; } +float PidFiller::GetMomentumTrd(const AnalysisTree::BranchChannel& trd_track) { + float momentum = 0.0; + if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("p")) > 0)) + momentum = trd_track.Value(trd_tracks_.GetField("p")); + else if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("p_out")) > 0)) + momentum = trd_track.Value(trd_tracks_.GetField("p_out")); + else + Warning("Exec", "Could not assign any momentum to the track, use p=0."); + + return momentum; +} + +float PidFiller::GetPtTrd(const AnalysisTree::BranchChannel& trd_track) { + float pT = 0; + if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("pT")) > 0)) + pT = trd_track.Value(trd_tracks_.GetField("pT")); + else if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("pT_out")) > 0)) + pT = trd_track.Value(trd_tracks_.GetField("pT_out")); + else + Warning("Exec", "Could not assign any pT to the track, use pT=0."); + return pT; +} + +int PidFiller::GetCharge(const AnalysisTree::BranchChannel& rec_track) { + return rec_track.Value(rec_tracks_.GetField("q")); +} + +void PidFiller::GetEnergyLossHitsTrd(const AnalysisTree::BranchChannel& trd_track, int& nhits_trd, std::array& dEdx) { + + nhits_trd = 0; + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) { + if (trd_track.Value(trd_dEdx_field_.at(ihit)) > 0) { + dEdx.at(ihit) = trd_track.Value(trd_dEdx_field_.at(ihit)); + nhits_trd++; + } else + dEdx.at(ihit) = 0; + } +} + +bool PidFiller::IsRichElectron(const AnalysisTree::BranchChannel& rec_track) { + bool isElectron = false; + int i_rich = rec_to_rich_->GetMatch(rec_track.GetId()); + if (i_rich > 0) { + const auto& rich_ring = rich_rings_[i_rich]; + float Aaxis = rich_ring.Value(rich_rings_.GetField("axis_a")); + float Baxis = rich_ring.Value(rich_rings_.GetField("axis_b")); + Double_t dist = 0;// richRing->GetDistance(); + + Float_t mom = GetMomentumTrd(rec_track); + Double_t MeanA = 4.95; + Double_t MeanB = 4.54; + Double_t RmsA = 0.30; + Double_t RmsB = 0.22; + Double_t RmsCoeff = 3.5; + Double_t DistCut = 1.; + + if (mom < 5.) { + if (fabs(Aaxis - MeanA) < RmsCoeff * RmsA && fabs(Baxis - MeanB) < RmsCoeff * RmsB && dist < DistCut) + isElectron = true; + } else { + ///2 sigma + Double_t polAaxis = 5.64791 - 4.24077 / (mom - 3.65494); + Double_t polBaxis = 5.41106 - 4.49902 / (mom - 3.52450); + if (Aaxis < (MeanA + RmsCoeff * RmsA) && Aaxis > polAaxis && Baxis < (MeanB + RmsCoeff * RmsB) && Baxis > polBaxis && dist < DistCut) + isElectron = true; + } + } + return isElectron; +} + void PidFiller::Init() { + // Settings for Trd Pid + if (is_run_pid_trd_ == true) { + getter_trd_->SetMinHits(nhits_min_); + getter_trd_->SetTruncationMode(trunc_mode_); + getter_trd_->SetProbabiltyMode(prob_mode_); + } + auto man = TaskManager::GetInstance(); auto chain = man->GetChain(); - rec_tracks_ = chain->GetBranch(rec_tracks_name_); - tof_hits_ = chain->GetBranch(tof_hits_name_); - pid_match_ = chain->GetMatchPointers().find({rec_tracks_name_ + "2" + tof_hits_name_})->second; + rec_tracks_ = chain->GetBranch(rec_tracks_name_); in_branches_.emplace(rec_tracks_name_); - in_branches_.emplace(tof_hits_name_); - auto conf = rec_tracks_.GetConfig().Clone(out_branch_name_, AnalysisTree::DetType::kParticle); - std::vector names{}; - for (const auto& pid : pid_codes_) { - names.push_back("prob_" + pid.second); + // Set and configure branches for Tof Pid + if (is_run_pid_tof_ == true) { + tof_hits_ = chain->GetBranch(tof_hits_name_); + rec_to_tof_ = chain->GetMatchPointers().find({rec_tracks_name_ + "2" + tof_hits_name_})->second; + + in_branches_.emplace(tof_hits_name_); + + std::vector names{}; + for (const auto& pid : pid_codes_tof_) { + names.push_back("prob_" + pid.second); + } + conf.AddFields(names, "probability to be proton, pion, kaon etc"); + } + + // Set and configure branches for Trd Pid + if (is_run_pid_trd_ == true) { + + trd_tracks_ = chain->GetBranchObject(trd_tracks_name_); + rich_rings_ = chain->GetBranchObject(rich_rings_name_); + rec_to_trd_ = chain->GetMatching(rec_tracks_name_, trd_tracks_name_); + rec_to_rich_ = chain->GetMatching(rec_tracks_name_, rich_rings_name_); + + in_branches_.emplace(trd_tracks_name_); + in_branches_.emplace(rich_rings_name_); + + conf.AddField("electron_rich", "rich electron hypothesis"); + conf.AddField("pid_trd", "trd pid hypothesis"); + conf.AddField("nhits_trd_eloss", "number of trd hits dEdx > 0"); + + std::vector names{}; + for (const auto& pid : pid_codes_trd_) { + string name = Form("prob_trd_%s", pid.second.Data()); + names.push_back(Form("prob_trd_%s", pid.second.Data())); + } + conf.AddFields(names, "trd probability to be proton, pion, kaon etc"); } - conf.AddFields(names, "probability to be proton, pion, kaon etc"); ana_tracks_ = Branch(conf); ana_tracks_.SetMutable(); @@ -50,23 +182,49 @@ void PidFiller::Init() { man->AddBranch(&ana_tracks_); - int i{0}; - auto match_br = {"SimParticles", "TofHits"}; + // Add matching for ouput + int imatch{0}; + std::vector match_br; + match_br.push_back("SimParticles"); + + if (is_run_pid_tof_ == true) + match_br.push_back("TofHits"); + + if (is_run_pid_trd_ == true) { + match_br.push_back("TrdTracks"); + match_br.push_back("RichRings"); + } + out_matches_.assign(match_br.size(), nullptr); for (const auto& br : match_br) { in_matches_.emplace_back(chain->GetMatchPointers().find({rec_tracks_name_ + "2" + br})->second); - man->AddMatching(out_branch_name_, br, out_matches_.at(i)); - i++; + man->AddMatching(out_branch_name_, br, out_matches_.at(imatch)); + imatch++; + } + + // Set input variables for Tof Pid + if (is_run_pid_tof_ == true) { + qp_tof_field_ = tof_hits_.GetField("qp_tof"); + q_field_ = rec_tracks_.GetField("q"); + mass2_field_ = tof_hits_.GetField("mass2"); + pid_tof_field_ = ana_tracks_.GetField("pid"); + for (const auto& pid : pid_codes_tof_) { + prob_tof_field_.push_back(ana_tracks_.GetField("prob_" + pid.second)); + } } - qp_tof_field_ = tof_hits_.GetField("qp_tof"); - q_field_ = rec_tracks_.GetField("q"); - mass2_field_ = tof_hits_.GetField("mass2"); - pid_field_ = ana_tracks_.GetField("pid"); + // Set input variables for Trd Pid + if (is_run_pid_trd_ == true) { + nhits_trd_eloss_field_ = ana_tracks_.GetField("nhits_trd_eloss"); + pid_trd_field_ = ana_tracks_.GetField("pid_trd"); + el_rich_field_ = ana_tracks_.GetField("electron_rich"); - for (const auto& pid : pid_codes_) { - prob_field_.push_back(ana_tracks_.GetField("prob_" + pid.second)); + for (const auto& pid : pid_codes_trd_) { + prob_trd_field_.push_back(ana_tracks_.GetField(Form("prob_trd_%s", pid.second.Data()))); + } + for (int i = 0; i < NumberOfTrdLayers; i++) + trd_dEdx_field_.push_back(trd_tracks_.GetField(("energy_loss_" + std::to_string(i)).c_str())); } } @@ -77,32 +235,87 @@ void PidFiller::Exec() { const auto& track = rec_tracks_[i]; auto particle = ana_tracks_.NewChannel(); particle.CopyContent(track); - auto q = track[q_field_]; - auto hit_id = pid_match_->GetMatch(i); - if (hit_id >= 0) { - const auto& tof_hit = tof_hits_[hit_id]; - auto pq = tof_hit[qp_tof_field_]; - auto m2 = tof_hit[mass2_field_]; + // Fill in Tof pid + if (is_run_pid_tof_ == true) { + auto q = track[q_field_]; - auto pid = getter_->GetPid(pq, m2, 0.5); - if (pid == 1 && q < 0) { - pid = -1; - } + auto hit_id = rec_to_tof_->GetMatch(i); + if (hit_id >= 0) { + const auto& tof_hit = tof_hits_[hit_id]; + auto pq = tof_hit[qp_tof_field_]; + auto m2 = tof_hit[mass2_field_]; + + auto pid = getter_tof_->GetPid(pq, m2, 0.5); + if (pid == 1 && q < 0) { + pid = -1; + } - particle.SetValue(pid_field_, pid); - auto prob = getter_->GetBayesianProbability(pq, m2); + particle.SetValue(pid_tof_field_, pid); + auto prob = getter_tof_->GetBayesianProbability(pq, m2); - int specie{0}; - for (const auto& pdg : pid_codes_) { - particle.SetValue(prob_field_.at(specie++), prob[pdg.first * signum(q)]);// Think what to do in case of electrons and muons + int specie{0}; + for (const auto& pdg : pid_codes_tof_) { + particle.SetValue(prob_tof_field_.at(specie++), prob[pdg.first * signum(q)]);// Think what to do in case of electrons and muons + } + } else { + int specie{0}; + for (const auto& pdg : pid_codes_tof_) { + particle.SetValue(prob_tof_field_.at(specie++), -1.f); + } + particle.SetValue(pid_tof_field_, 2 * signum(q)); } - } else { - int specie{0}; - for (const auto& pdg : pid_codes_) { - particle.SetValue(prob_field_.at(specie++), -1.f); + } + + // Fill in Trd pid + if (is_run_pid_trd_ == true) { + + int itrd = rec_to_trd_->GetMatch(i); + int nhits_trd = 0; + + if (itrd > -1) { + const auto& trd_track = trd_tracks_[itrd]; + float mom = GetMomentumTrd(trd_track); + float pT = GetPtTrd(trd_track); + int charge = GetCharge(track); + + std::array dEdx_hits = {0.0, 0.0, 0.0, 0.0}; + GetEnergyLossHitsTrd(trd_track, nhits_trd, dEdx_hits); + + TrdContainer trdtrack(mom, pT, charge, nhits_trd, dEdx_hits); + trdtrack.ScaleEnergyLossLength(); + + int pidtrd_hypo; + std::map prob_trd; + if (nhits_trd >= nhits_min_) { + if (prob_mode_ == 0) + prob_trd = getter_trd_->GetTrdProbabilities(trdtrack); + if (prob_mode_ == 1) + prob_trd = getter_trd_->GetTrdProbabilitiesMulti(trdtrack); + + int specie{0}; + for (const auto& pdg : pid_codes_trd_) + particle.SetValue(prob_trd_field_.at(specie++), prob_trd[pdg.first]); + + pidtrd_hypo = getter_trd_->GetTrdPid(prob_trd, purity_, charge); + particle.SetValue(pid_trd_field_, pidtrd_hypo); + } else { + particle.SetValue(pid_trd_field_, -2); + int specie{0}; + for (const auto& pdg : pid_codes_trd_) + particle.SetValue(prob_trd_field_.at(specie++), -1.f); + } + } else { + particle.SetValue(pid_trd_field_, -2); + int specie{0}; + for (const auto& pdg : pid_codes_trd_) + particle.SetValue(prob_trd_field_.at(specie++), -1.f); } - particle.SetValue(pid_field_, 2 * signum(q)); + + particle.SetValue(nhits_trd_eloss_field_, nhits_trd); + + bool isElectron = IsRichElectron(track); + particle.SetValue(el_rich_field_, isElectron); } } int i{0}; diff --git a/at_interface/PidFiller.hpp b/at_interface/PidFiller.hpp index d766a36..7960829 100644 --- a/at_interface/PidFiller.hpp +++ b/at_interface/PidFiller.hpp @@ -1,17 +1,25 @@ -#ifndef PID_AT_INTERFACE_PIDFILLER_HPP_ -#define PID_AT_INTERFACE_PIDFILLER_HPP_ +#ifndef PID_INTERFACE_PIDFILLER_HPP_ +#define PID_INTERFACE_PIDFILLER_HPP_ -#include "Constants.h" -#include "Getter.h" +#include "ConstantsTof.h" +#include "ConstantsTrd.h" +#include "ContainerTrd.h" +#include "GetterTof.h" +#include "GetterTrd.h" #include "AnalysisTree/Task.hpp" #include "AnalysisTree/TaskManager.hpp" +#include "TH2F.h" + class PidFiller : public AnalysisTree::Task { public: - PidFiller(const std::string& file, const std::string& getter); - ~PidFiller() override { delete getter_; } + PidFiller(const std::string& pid_file_name_tof, const std::string& pid_file_name_trd, const std::string& getter_name_tof, const std::string& getter_name_trd, int pid_mode); + ~PidFiller() override { + delete getter_tof_; + delete getter_trd_; + } void Init() override; void Exec() override; @@ -20,34 +28,99 @@ class PidFiller : public AnalysisTree::Task { // man->RemoveBranch(rec_tracks_name_); } + void SetRecTracksName(const std::string& name) { rec_tracks_name_ = name; } + void SetTofHitsName(const std::string& name) { tof_hits_name_ = name; } + void SetTrdTracksName(const std::string& name) { trd_tracks_name_ = name; }; + void SetRichRingsName(const std::string& name) { rich_rings_name_ = name; }; + + // Settings for Trd Pid + void SetMinHits(int nhits_min) {// Min. number of hits per track + if (nhits_min < 1 || nhits_min > 4) + throw std::runtime_error("Wrong number of minimum hits: Number of minimum hits need to be between 1 and 4!"); + else + nhits_min_ = nhits_min; + } + void SetTruncationMode(int trunc_mode) {// Calculation of energy loss for up to 4 layers: + // =0: average over all hits + if (trunc_mode < 0 || trunc_mode > 4) // =1-4: Select hits with lowest dEdx: 1 hit, =2: 2 hits, =3: 3 hits, =4: 4 hits + throw std::runtime_error("Wrong index for truncation mode: Truncation mode needs to be between 0 and 4!"); + else + trunc_mode_ = trunc_mode; + } + + void SetProbabilityMode(int prob_mode) {// Probability for particle species i: + // =0: total probability - probability based on particle multiplicites i + if (prob_mode < 0 || prob_mode > 1) // =1: likelihood - probability based on dEdx-distribution of particle + throw std::runtime_error("Wrong index for probability mode: Truncation mode needs to be between 0 and 4!"); + else + prob_mode_ = prob_mode; + } + + void SetPurity(const float purity) {// Minium purity for pid-hypothesis + if (purity < 0.0 || purity > 1.0) + throw std::runtime_error("Wrong value " + std::to_string(purity) + " for purity: Purity needs to be between 0.0 and 1.0!"); + else + purity_ = purity; + } + protected: int signum(int x) const; + float GetMomentumTrd(const AnalysisTree::BranchChannel& trd_particle); + float GetPtTrd(const AnalysisTree::BranchChannel& trd_particle); + int GetCharge(const AnalysisTree::BranchChannel& trd_particle); + void GetEnergyLossHitsTrd(const AnalysisTree::BranchChannel& trd_track, int& nhits_trd, std::array& dEdx); + bool IsRichElectron(const AnalysisTree::BranchChannel& rec_particle); + AnalysisTree::Branch rec_tracks_; AnalysisTree::Branch tof_hits_; + AnalysisTree::Branch trd_tracks_; + AnalysisTree::Branch rich_rings_; AnalysisTree::Branch ana_tracks_; - AnalysisTree::Matching* pid_match_{nullptr}; + AnalysisTree::Matching* rec_to_tof_{nullptr}; + AnalysisTree::Matching* rec_to_trd_{nullptr}; + AnalysisTree::Matching* rec_to_rich_{nullptr}; + // Tof fields AnalysisTree::Field qp_tof_field_; AnalysisTree::Field q_field_; AnalysisTree::Field mass2_field_; - AnalysisTree::Field pid_field_; - std::vector prob_field_{}; + AnalysisTree::Field pid_tof_field_; + std::vector prob_tof_field_{}; + + // Trd fields + std::vector trd_dEdx_field_{}; + AnalysisTree::Field el_rich_field_; + AnalysisTree::Field pid_trd_field_; + AnalysisTree::Field nhits_trd_eloss_field_; + std::vector prob_trd_field_{}; std::vector in_matches_{}; std::vector out_matches_{}; - std::string rec_tracks_name_{"VtxTracks"}; // Branch with input tracks - std::string tof_hits_name_{"TofHits"}; // Branch with TOF info (m2) - std::string out_branch_name_{"RecParticles"}; // Output branch (based on VtxTracks) with pid info: probabilities and particle type hypothesis - - Pid::Getter* getter_{nullptr}; - std::vector> pid_codes_{ - {PidParticles::kProton, "p"}, - {PidParticles::kPionPos, "pi"}, - {PidParticles::kKaonPos, "K"}, - {PidParticles::kDeutron, "d"}, - {PidParticles::kBgPos, "bg"}}; + std::string rec_tracks_name_{"VtxTracks"}; // Branch with input tracks + std::string tof_hits_name_{"TofHits"}; // Branch with TOF info (m2) + std::string trd_tracks_name_{"TrdTracks"}; // Branch with Trd info (dEdx) + std::string rich_rings_name_{"RichRings"}; // Branch with Rich info + std::string out_branch_name_{"RecParticles"};// Output branch (based on VtxTracks) with pid info: probabilities and particle type hypothesis + + bool is_run_pid_tof_{false}; + bool is_run_pid_trd_{false}; + + int trunc_mode_{0}; + int prob_mode_{0}; + float purity_{0.0}; + int nhits_min_{1}; + + Pid::GetterTof* getter_tof_{nullptr}; + std::vector> pid_codes_tof_{ + {PidTofParticles::kProton, "p"}, + {PidTofParticles::kPionPos, "pi"}, + {PidTofParticles::kKaonPos, "K"}, + {PidTofParticles::kDeutron, "d"}, + {PidTofParticles::kBgPos, "bg"}}; + + PidTrd::GetterTrd* getter_trd_{nullptr}; }; -#endif// PID_AT_INTERFACE_PIDFILLER_HPP_ +#endif// PID_INTERFACE_PIDFILLER_HPP_ diff --git a/at_interface/PidTrdMcInput.cpp b/at_interface/PidTrdMcInput.cpp new file mode 100644 index 0000000..64e96c1 --- /dev/null +++ b/at_interface/PidTrdMcInput.cpp @@ -0,0 +1,397 @@ +#include "PidTrdMcInput.hpp" +#include "ConstantsTrd.h" + +using namespace AnalysisTree; + +using std::to_string; + +float PidTrdMcInput::GetMomentum(const AnalysisTree::BranchChannel& trd_track) { + float momentum = 0; + if (TMath::Abs(trd_track.Value(rec_tracks_.GetField("p")) > 0)) + momentum = trd_track.Value(rec_tracks_.GetField("p")); + else if (TMath::Abs(trd_track.Value(rec_tracks_.GetField("p_out")) > 0)) + momentum = trd_track.Value(rec_tracks_.GetField("p_out")); + else + Warning("Exec", "Could not assign any momentum to the track, use p=0."); + return momentum; +} + +float PidTrdMcInput::GetPt(const AnalysisTree::BranchChannel& trd_track) { + float pT = 0; + if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("pT")) > 0)) + pT = trd_track.Value(trd_tracks_.GetField("pT")); + else if (TMath::Abs(trd_track.Value(trd_tracks_.GetField("pT_out")) > 0)) + pT = trd_track.Value(trd_tracks_.GetField("pT_out")); + else + Warning("Exec", "Could not assign any pT to the track, use pT=0."); + return pT; +} + +int PidTrdMcInput::GetMcPdg(const AnalysisTree::BranchChannel& rec_track) { + const int sim_id = rec_to_sim_->GetMatch(rec_track.GetId()); + int pdg; + if (sim_id < 0) pdg = -1; + else + pdg = sim_tracks_[sim_id].Value(sim_tracks_.GetField("pid")); + return pdg; +} + +int PidTrdMcInput::GetCharge(const AnalysisTree::BranchChannel& rec_track) { + int charge = rec_track.Value(rec_tracks_.GetField("q")); + return charge; +} + +void PidTrdMcInput::GetEnergyLossHits(const AnalysisTree::BranchChannel& trd_track, int& nhits_trd, array& dEdx) { + + nhits_trd = 0; + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) { + if (trd_track.Value(trd_dEdx_field_.at(ihit)) > 0) { + dEdx.at(ihit) = trd_track.Value(trd_dEdx_field_.at(ihit)); + nhits_trd++; + } else + dEdx.at(ihit) = 0; + } +} + +void PidTrdMcInput::FillHistogram(TrdContainer track) { + + float mom = track.GetP(); + int nhits_trd = track.GetNhitsTrd(); + int mc_pdg = track.GetMcPdg(); + int charge = track.GetCharge(); + + // Histograms for dEdx of track + for (int imode = 0; imode < nhits_trd; imode++) { + + float dEdx = track.GetdEdxTrack(imode); + if (dEdx <= 0) continue; + + outFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits_trd - 1) + "/reco_info"); + if (charge > 0) + h2dEdx_p_pos_[imode]->Fill(mom, dEdx); + if (charge < 0) + h2dEdx_p_neg_[imode]->Fill(mom, dEdx); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + if (TMath::Abs(mc_pdg) == pid_codes_trd_.at(ipdg).first) { + outFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits_trd - 1) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + if (charge > 0) + h2dEdx_p_pdg_[imode].at(ipdg)->Fill(mom, dEdx); + if (charge < 0) + h2dEdx_p_pdg_[imode].at(NumberOfPidsTrd - 1 + ipdg)->Fill(mom, dEdx); + } + } + } + + // Histograms for dEdx of hits + std::array dEdx_hits_sorted = track.GetdEdxHitsSorted(); + + // All hits in one plot for every particle + outFile_->cd(dirname_hits_ + "/reco_info"); + for (int ihit = 0; ihit < nhits_trd; ihit++) { + if (charge > 0) + h2dEdx_hits_all_p_pos_->Fill(mom, dEdx_hits_sorted.at(ihit)); + if (charge < 0) + h2dEdx_hits_all_p_neg_->Fill(mom, dEdx_hits_sorted.at(ihit)); + } + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + if (TMath::Abs(mc_pdg) == pid_codes_trd_.at(ipdg).first) { + outFile_->cd(dirname_hits_ + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + for (int ihit = 0; ihit < nhits_trd; ihit++) { + if (charge > 0) + h2dEdx_hits_all_p_pdg_.at(ipdg)->Fill(mom, dEdx_hits_sorted.at(ihit)); + if (charge < 0) + h2dEdx_hits_all_p_pdg_.at(NumberOfPidsTrd - 1 + ipdg)->Fill(mom, dEdx_hits_sorted.at(ihit)); + } + } + } + for (int imode = 0; imode < nhits_trd; imode++) { + outFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits_trd - 1) + "/reco_info"); + for (int ihit = 0; ihit < imode + 1; ihit++) { + if (charge > 0) + h2dEdx_hits_p_pos_[imode]->Fill(mom, dEdx_hits_sorted.at(ihit)); + if (charge < 0) + h2dEdx_hits_p_neg_[imode]->Fill(mom, dEdx_hits_sorted.at(ihit)); + } + + // Hits seperated by number of hits per track & truncation mode + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + if (TMath::Abs(mc_pdg) == pid_codes_trd_.at(ipdg).first) { + outFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits_trd - 1) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + for (int ihit = 0; ihit < imode + 1; ihit++) { + if (charge > 0) + h2dEdx_hits_p_pdg_[imode].at(ipdg)->Fill(mom, dEdx_hits_sorted.at(ihit)); + if (charge < 0) + h2dEdx_hits_p_pdg_[imode].at(NumberOfPidsTrd - 1 + ipdg)->Fill(mom, dEdx_hits_sorted.at(ihit)); + } + } + } + } +} + +void PidTrdMcInput::CreateMcHistograms(int nhits) { + + TString histname; + + // Histgrams for dEdx of tracks + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + outFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits - 1) + "/reco_info"); + + h2dEdx_p_pos_[imode] = (TH2F*) gDirectory->Get(histnames_all_pos_.at(imode)); + h2dEdx_p_neg_[imode] = (TH2F*) gDirectory->Get(histnames_all_neg_.at(imode)); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + outFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits - 1) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + + histname = histnames_pos_.at(imode).at(ipdg); + h2dEdx_p_pdg_[imode].at(ipdg) = (TH2F*) gDirectory->Get(histname); + + histname = histnames_neg_.at(imode).at(ipdg); + h2dEdx_p_pdg_[imode].at(NumberOfPidsTrd - 1 + ipdg) = (TH2F*) gDirectory->Get(histname); + } + } + + // Histgrams for dEdx of hits + outFile_->cd(dirname_hits_ + "/reco_info"); + h2dEdx_hits_all_p_pos_ = (TH2F*) gDirectory->Get(histnames_hits_all_pos_); + h2dEdx_hits_all_p_neg_ = (TH2F*) gDirectory->Get(histnames_hits_all_neg_); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + outFile_->cd(dirname_hits_ + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + + histname = histnames_hits_pos_.at(ipdg); + h2dEdx_hits_all_p_pdg_.at(ipdg) = (TH2F*) gDirectory->Get(histname); + + histname = histnames_hits_neg_.at(ipdg); + h2dEdx_hits_all_p_pdg_.at(NumberOfPidsTrd - 1 + ipdg) = (TH2F*) gDirectory->Get(histname); + } + + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + outFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits - 1) + "/reco_info"); + + h2dEdx_hits_p_pos_[imode] = (TH2F*) gDirectory->Get(histnames_all_pos_.at(imode)); + h2dEdx_hits_p_neg_[imode] = (TH2F*) gDirectory->Get(histnames_all_neg_.at(imode)); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + outFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits - 1) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + + histname = histnames_pos_.at(imode).at(ipdg); + h2dEdx_hits_p_pdg_[imode].at(ipdg) = (TH2F*) gDirectory->Get(histname); + + histname = histnames_neg_.at(imode).at(ipdg); + h2dEdx_hits_p_pdg_[imode].at(NumberOfPidsTrd - 1 + ipdg) = (TH2F*) gDirectory->Get(histname); + } + } + + auto it = trd_container_.find(nhits); + if (it != trd_container_.end()) { + for (auto track : it->second) { + FillHistogram(track); + } + } + outFile_->cd(); + outFile_->Write("", TObject::kOverwrite); +} + +void PidTrdMcInput::InitMcHistograms() { + + outFile_ = new TFile(outfilename_, "RECREATE"); + + // Create directories + TDirectory *directory, *directory1, *directory2, *directory3; + for (int idir = 0; idir < 2; idir++) { + if (idir == 0) directory = outFile_->mkdir(dirname_tracks_); + if (idir == 1) { + directory = outFile_->mkdir(dirname_hits_); + directory1 = directory->mkdir("reco_info"); + directory1 = directory->mkdir("reco_vs_sim_info"); + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + directory2 = directory1->mkdir(pid_codes_trd_.at(ipdg).second); + } + } + for (int inhits = 0; inhits < NumberOfTrdLayers; inhits++) { + directory1 = directory->mkdir(dirname_nhits_.at(inhits)); + directory2 = directory1->mkdir("reco_info"); + directory2 = directory1->mkdir("reco_vs_sim_info"); + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + directory3 = directory2->mkdir(pid_codes_trd_.at(ipdg).second); + } + } + } + + // Initialize histograms + TString dirname, histname, histtitle; + TH2F* histogram; + + for (int idir = 0; idir < 2; idir++) { + + if (idir == 0) dirname = dirname_tracks_; + + if (idir == 1) { + + dirname = dirname_hits_; + + outFile_->cd(dirname + "/reco_info"); + histname = histnames_hits_all_pos_; + histtitle = "dEdx hits vs. p_{rec} all+"; + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + histname = histnames_hits_all_neg_; + histtitle = "dEdx hits vs. p_{rec} all- "; + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + + outFile_->cd(dirname + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + + histname = histnames_hits_pos_.at(ipdg); + histtitle = Form("dEdx hits : p_{rec} (%s+)", particlename.Data()); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + histname = histnames_hits_neg_.at(ipdg); + histtitle = Form("dEdx hits : p_{rec} (%s-)", particlename.Data()); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + } + } + + for (int inhits = 0; inhits < NumberOfTrdLayers; inhits++) { + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + + if (imode > inhits) continue; + + outFile_->cd(dirname + "/" + dirname_nhits_.at(inhits) + "/reco_info"); + + histname = histnames_all_pos_.at(imode); + histtitle = "dEdx vs. p_{rec} all+ " + to_string(inhits + 1) + " Trd Hits - " + histtitle_mode_.at(imode); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + histname = histnames_all_neg_.at(imode); + histtitle = "dEdx vs. p_{rec} all- " + to_string(inhits + 1) + " Trd Hits - " + histtitle_mode_.at(imode); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + + outFile_->cd(dirname + "/" + dirname_nhits_.at(inhits) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + + histname = histnames_pos_.at(imode).at(ipdg); + histtitle = Form("dEdx : p_{rec} (%s+) %d Trd Hits - %s", particlename.Data(), inhits + 1, histtitle_mode_.at(imode).Data()); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + + histname = histnames_neg_.at(imode).at(ipdg); + histtitle = Form("dEdx : p_{rec} (%s-) %d Trd Hits - %s", particlename.Data(), inhits + 1, histtitle_mode_.at(imode).Data()); + histogram = new TH2F(histname, histtitle, nbins_mom_, bins_mom_, nbins_dEdx_, bins_dEdx_); + histogram->Write(); + } + } + } + } +} + +void PidTrdMcInput::OpenMcHistograms() { + outFile_ = new TFile(outfilename_, "UPDATE"); +} + +void PidTrdMcInput::Init() { + auto man = TaskManager::GetInstance(); + auto chain = man->GetChain(); + + chain->InitPointersToBranches({}); + rec_tracks_ = chain->GetBranchObject(rec_tracks_name_); + trd_tracks_ = chain->GetBranchObject(trd_tracks_name_); + sim_tracks_ = chain->GetBranchObject(sim_tracks_name_); + rec_to_trd_ = chain->GetMatching(rec_tracks_name_, trd_tracks_name_); + rec_to_sim_ = chain->GetMatching(rec_tracks_name_, sim_tracks_name_); + + for (int i = 0; i < NumberOfTrdLayers; i++) + trd_dEdx_field_.push_back(trd_tracks_.GetField(("energy_loss_" + to_string(i)).c_str())); + + Float_t bw_mom, bw_dEdx; + + if (update_mchisto_ == kFALSE) { + bw_mom = BwMom; + nbins_mom_ = NbinsMom; + bins_mom_[0] = 0.; + for (int i = 0; i < nbins_mom_; i++) + bins_mom_[i + 1] = TMath::Power(10, -1.0 + i * bw_mom); + + bw_dEdx = BwdEdx; + nbins_dEdx_ = NbinsdEdx; + bins_dEdx_[0] = 0.; + for (int i = 0; i < nbins_dEdx_; i++) + bins_dEdx_[i + 1] = TMath::Power(10, -1.0 + i * bw_dEdx); + } + + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + histnames_all_pos_.at(imode) = "h2dEdx_p_pos_" + to_string(imode); + histnames_all_neg_.at(imode) = "h2dEdx_p_neg_" + to_string(imode); + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + histnames_pos_.at(imode).at(ipdg) = "h2dEdx_p_" + particlename + "_pos_" + to_string(imode); + histnames_neg_.at(imode).at(ipdg) = "h2dEdx_p_" + particlename + "_neg_" + to_string(imode); + } + + histnames_hits_all_pos_ = "h2dEdx_p_pos"; + histnames_hits_all_neg_ = "h2dEdx_p_neg"; + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + histnames_hits_pos_.at(ipdg) = "h2dEdx_p_" + particlename + "_pos"; + histnames_hits_neg_.at(ipdg) = "h2dEdx_p_" + particlename + "_neg"; + } + } + + if (update_mchisto_ == kTRUE) + OpenMcHistograms(); + else + InitMcHistograms(); +} + +void PidTrdMcInput::Exec() { + trd_container_.clear(); + + for (int itrack = 0; itrack < rec_tracks_.size(); ++itrack) { + const auto& rec_track = rec_tracks_[itrack]; + int itrd = rec_to_trd_->GetMatch(rec_track.GetId()); + if (itrd < 0) continue; + const auto& trd_track = trd_tracks_[itrd]; + + float mom = GetMomentum(trd_track); + float pT = GetPt(trd_track); + int mc_pdg = GetMcPdg(rec_track); + int charge = GetCharge(rec_track); + array dEdx_hits = {0.0, 0.0, 0.0, 0.0}; + int nhits_trd = 0; + + GetEnergyLossHits(trd_track, nhits_trd, dEdx_hits); + + if (nhits_trd < 1) continue; + + TrdContainer trdtrack(mom, pT, charge, nhits_trd, dEdx_hits, mc_pdg); + trdtrack.ScaleEnergyLossLength(); + trdtrack.CalculateEnergyLossTrackAllModes(); + + auto it = trd_container_.find(nhits_trd); + if (it != trd_container_.end()) { + it->second.emplace_back(trdtrack); + } else { + trd_container_[nhits_trd] = {static_cast(trdtrack)}; + } + } + + for (int inhits = 0; inhits < NumberOfTrdLayers; inhits++) { + CreateMcHistograms(inhits + 1); + } +} + +void PidTrdMcInput::Finish() { + outFile_->Close(); +} diff --git a/at_interface/PidTrdMcInput.hpp b/at_interface/PidTrdMcInput.hpp new file mode 100644 index 0000000..4ca28de --- /dev/null +++ b/at_interface/PidTrdMcInput.hpp @@ -0,0 +1,89 @@ +#ifndef PIDTRD_MCINPUT_HPP_ +#define PIDTRD_MCINPUT_HPP_ + +#include "TH2F.h" + +#include "AnalysisTree/Task.hpp" +#include "AnalysisTree/TaskManager.hpp" + +#include "ContainerTrd.h" + +using std::array; + +class PidTrdMcInput : public AnalysisTree::Task { + + public: + PidTrdMcInput(const std::string& outfilename) { + outfilename_ = outfilename; + }; + ~PidTrdMcInput() override = default; + + void Init() override; + void InitMcHistograms(); + void OpenMcHistograms(); + void Exec() override; + void Finish() override; + + void SetRecTracksName(const std::string& name) { rec_tracks_name_ = name; }; + void SetSimTracksName(const std::string& name) { sim_tracks_name_ = name; }; + void SetTrdTracksName(const std::string& name) { trd_tracks_name_ = name; }; + void SetUpdateMcHistos(const Bool_t update_mchisto) { update_mchisto_ = update_mchisto; }; + + protected: + float GetMomentum(const AnalysisTree::BranchChannel& trd_track); + float GetPt(const AnalysisTree::BranchChannel& trd_particle); + int GetMcPdg(const AnalysisTree::BranchChannel& rec_track); + int GetCharge(const AnalysisTree::BranchChannel& rec_track); + void GetEnergyLossHits(const AnalysisTree::BranchChannel& trd_track, int& nhits_trd, array& dEdx); + void CalculateEnergyLoss(TrdContainer track, float& dEdx_track, int trunc_mode); + void CreateMcHistograms(int trunc_mode); + void FillHistogram(TrdContainer track, float dEdx); + void FillHistogram(TrdContainer track); + bool IsRichElectron(const AnalysisTree::BranchChannel& rec_particle); + + AnalysisTree::Branch rec_tracks_; + AnalysisTree::Branch sim_tracks_; + AnalysisTree::Branch trd_tracks_; + + AnalysisTree::Matching* rec_to_sim_{nullptr}; + AnalysisTree::Matching* rec_to_trd_{nullptr}; + + std::string rec_tracks_name_{"RecTracks"}; // Branch with input tracks + std::string trd_tracks_name_{"TrdTracks"}; // Branch with Trd info (dEdx) + std::string sim_tracks_name_{"SimParticles"};// Branch with Mc info + + std::vector trd_dEdx_field_{}; + + TString outfilename_; + TFile* outFile_; + + // MC histograms + + TString histnames_hits_all_pos_, histnames_hits_all_neg_; + array histnames_hits_pos_, histnames_hits_neg_; + array histnames_all_pos_, histnames_all_neg_; + array, NumberOfTrdLayers> histnames_pos_, histnames_neg_; + + TString outfile_info_{""}; + + TH2F* h2dEdx_p_pos_[NumberOfTruncMode - 1]; + TH2F* h2dEdx_p_neg_[NumberOfTruncMode - 1]; + array h2dEdx_p_pdg_[NumberOfTrdLayers]; + TH2F* h2dEdx_hits_all_p_pos_; + TH2F* h2dEdx_hits_all_p_neg_; + array h2dEdx_hits_all_p_pdg_; + TH2F* h2dEdx_hits_p_pos_[NumberOfTruncMode - 1]; + TH2F* h2dEdx_hits_p_neg_[NumberOfTruncMode - 1]; + array h2dEdx_hits_p_pdg_[NumberOfTrdLayers]; + + Int_t nbins_mom_; + Double_t bins_mom_[NbinsMax]; + Int_t nbins_dEdx_; + Double_t bins_dEdx_[NbinsMax]; + + Bool_t update_mchisto_{kFALSE}; + + std::map> trd_container_;//container with trd tracks +}; + +#endif//PIDTRD_AT_INTERFACE_PIDFILLER_HPP_ diff --git a/input/dEdx_p.jul25.phqmd.auau.12agev.root b/input/dEdx_p.jul25.phqmd.auau.12agev.root new file mode 100644 index 0000000..97b5761 Binary files /dev/null and b/input/dEdx_p.jul25.phqmd.auau.12agev.root differ diff --git a/input/dEdx_p.jul25.phqmd.auau.12agev_probabilities.root b/input/dEdx_p.jul25.phqmd.auau.12agev_probabilities.root new file mode 100644 index 0000000..a3af9c3 Binary files /dev/null and b/input/dEdx_p.jul25.phqmd.auau.12agev_probabilities.root differ diff --git a/interface/CMakeLists.txt b/interface/CMakeLists.txt new file mode 100644 index 0000000..e6f4aea --- /dev/null +++ b/interface/CMakeLists.txt @@ -0,0 +1,32 @@ +set(SOURCES + PidTrdRunGetter.cpp + ) + +string(REPLACE ".cpp" ".hpp" HEADERS "${SOURCES}") +include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/src_trd) +add_library(PidInt SHARED ${SOURCES} G__PidInt.cxx) + +message("${PROJECT_INCLUDE_DIRECTORIES}") + +ROOT_GENERATE_DICTIONARY(G__PidInt ${HEADERS} LINKDEF PidIntLinkDef.h) +target_link_libraries(PidInt ${ROOT_LIBRARIES} PidTrd) + +install( + FILES + ${HEADERS} + DESTINATION + include + COMPONENT + Devel +) + +set(PCM_FILE_NAME libPidInt) + +install( + FILES + "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}_rdict.pcm" + "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}.rootmap" + DESTINATION + lib + OPTIONAL +) diff --git a/interface/PidIntLinkDef.h b/interface/PidIntLinkDef.h new file mode 100644 index 0000000..90de817 --- /dev/null +++ b/interface/PidIntLinkDef.h @@ -0,0 +1,8 @@ +#ifndef PID_INTERFACE_PIDLINKDEF_H_ +#define PID_INTERFACE_PIDLINKDEF_H_ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#endif //PID_INTERFACE_PIDLINKDEF_H_ diff --git a/interface/PidTrdRunGetter.cpp b/interface/PidTrdRunGetter.cpp new file mode 100644 index 0000000..ac60d33 --- /dev/null +++ b/interface/PidTrdRunGetter.cpp @@ -0,0 +1,163 @@ +#include "PidTrdRunGetter.hpp" +#include "ConstantsTrd.h" +using std::to_string; + +void PidTrdRunGetter::CalculateProbabilities(int nhits) { + + TH2F* h2pos; + TH2F* h2neg; + TH2F* h2tmp; + TH2F* h2prob; + TString histname; + PidTrd::ParticleProb particleprob; + + for (int truncmode = 0; truncmode < nhits + 1; truncmode++) { + inFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits) + "/reco_info"); + h2pos = (TH2F*) gDirectory->Get(histnames_all_pos_.at(truncmode)); + h2neg = (TH2F*) gDirectory->Get(histnames_all_neg_.at(truncmode)); + + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second.Data(); + int charge = 1;// charge = 1: positive particles + + int probmode = 0;// total probability - probability based on particle multiplicites + inFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + h2tmp = (TH2F*) gDirectory->Get(histnames_pos_.at(truncmode).at(ipdg)); + h2prob = (TH2F*) h2tmp->Clone(histnames_prob_pos_.at(truncmode).at(ipdg)); + CalculateProbabilitiesTot(h2pos, h2tmp, h2prob); + h2prob->SetTitle(Form("dEdx : p_{rec} probability (T) (%s+) %d Trd Hits - %s", pid_codes_trd_.at(ipdg).second.Data(), nhits + 1, histtitle_mode_.at(truncmode).Data())); + particleprob.Update(ipdg, charge, nhits, truncmode, probmode, h2prob); + getter_trd_.AddParticleProb(particleprob); + if (write_mchistograms_out_ == kTRUE) + outFile_->cd(dirname_nhits_.at(nhits) + "/" + particlename + "/" + probname_.at(0)) && h2prob->Write("", TObject::kOverwrite); + + probmode = 1;// likelihood - probability based on dEdx-distribution of particle + inFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + h2tmp = (TH2F*) gDirectory->Get(histnames_pos_.at(truncmode).at(ipdg)); + h2prob = (TH2F*) h2tmp->Clone(histnames_prob_pos_.at(truncmode).at(ipdg)); + CalculateProbabilitiesLike(h2tmp, h2prob); + h2prob->SetTitle(Form("dEdx : p_{rec} probability (L) (%s+) %d Trd Hits - %s", pid_codes_trd_.at(ipdg).second.Data(), nhits + 1, histtitle_mode_.at(truncmode).Data())); + particleprob.Update(ipdg, charge, nhits, truncmode, probmode, h2prob); + getter_trd_.AddParticleProb(particleprob); + if (write_mchistograms_out_ == kTRUE) + outFile_->cd(dirname_nhits_.at(nhits) + "/" + particlename + "/" + probname_.at(1)) && h2prob->Write("", TObject::kOverwrite); + + charge = -1;// charge = -1: negative particles + probmode = 0; + inFile_->cd(dirname_tracks_ + "/" + dirname_nhits_.at(nhits) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + h2tmp = (TH2F*) gDirectory->Get(histnames_neg_.at(truncmode).at(ipdg)); + h2prob = (TH2F*) h2tmp->Clone(histnames_prob_neg_.at(truncmode).at(ipdg)); + CalculateProbabilitiesTot(h2neg, h2tmp, h2prob); + h2prob->SetTitle(Form("dEdx : p_{rec} probability (T) (%s-) %d Trd Hits - %s", pid_codes_trd_.at(ipdg).second.Data(), nhits + 1, histtitle_mode_.at(truncmode).Data())); + particleprob.Update(ipdg, charge, nhits, truncmode, probmode, h2prob); + getter_trd_.AddParticleProb(particleprob); + if (write_mchistograms_out_ == kTRUE) + outFile_->cd(dirname_nhits_.at(nhits) + "/" + particlename + "/" + probname_.at(0)) && h2prob->Write("", TObject::kOverwrite); + + probmode = 1; + inFile_->cd(dirname_hits_ + "/" + dirname_nhits_.at(nhits) + "/reco_vs_sim_info/" + pid_codes_trd_.at(ipdg).second); + h2tmp = (TH2F*) gDirectory->Get(histnames_neg_.at(truncmode).at(ipdg)); + h2prob = (TH2F*) h2tmp->Clone(histnames_prob_neg_.at(truncmode).at(ipdg)); + CalculateProbabilitiesLike(h2tmp, h2prob); + h2prob->SetTitle(Form("dEdx : p_{rec} probability (L) (%s-) %d Trd Hits - %s", pid_codes_trd_.at(ipdg).second.Data(), nhits + 1, histtitle_mode_.at(truncmode).Data())); + particleprob.Update(ipdg, charge, nhits, truncmode, probmode, h2prob); + getter_trd_.AddParticleProb(particleprob); + if (write_mchistograms_out_ == kTRUE) + outFile_->cd(dirname_nhits_.at(nhits) + "/" + particlename + "/" + probname_.at(1)) && h2prob->Write("", TObject::kOverwrite); + } + } + h2tmp->Delete(); +} + +void PidTrdRunGetter::CalculateProbabilitiesTot(TH2F* h2all, TH2F* h2tmp, TH2F*& h2prob) { + h2prob->Divide(h2all); +} + +void PidTrdRunGetter::CalculateProbabilitiesLike(TH2F* h2tmp, TH2F*& h2prob) { + for (Int_t ibinx = 1; ibinx < h2tmp->GetNbinsX() + 1; ibinx++) { + Double_t sumy = 0.; + for (Int_t ibiny = 1; ibiny < h2tmp->GetNbinsY() + 1; ibiny++) + sumy += h2tmp->GetBinContent(ibinx, ibiny); + for (Int_t ibiny = 1; ibiny < h2tmp->GetNbinsY() + 1; ibiny++) + if (sumy > 0) h2prob->SetBinContent(ibinx, ibiny, h2tmp->GetBinContent(ibinx, ibiny) / sumy); + } +} + +void PidTrdRunGetter::InitMcHistogramsOut() { + + TDirectory *directory, *directory1, *directory2; + + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + histnames_prob_pos_.at(imode).at(ipdg) = "h2dEdx_p_prob_" + particlename + "_pos_" + to_string(imode); + histnames_prob_neg_.at(imode).at(ipdg) = "h2dEdx_p_prob_" + particlename + "_neg_" + to_string(imode); + } + } + + for (int inhits = 0; inhits < NumberOfTrdLayers; inhits++) { + directory = outFile_->mkdir(dirname_nhits_.at(inhits)); + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + outFile_->cd(); + directory1 = directory->mkdir(pid_codes_trd_.at(ipdg).second); + directory2 = directory1->mkdir(probname_.at(0)); + directory2 = directory1->mkdir(probname_.at(1)); + } + } +} + +void PidTrdRunGetter::OpenMcHistograms() { + + for (int imode = 0; imode < NumberOfTruncMode - 1; imode++) { + histnames_all_pos_.at(imode) = "h2dEdx_p_pos_" + to_string(imode); + histnames_all_neg_.at(imode) = "h2dEdx_p_neg_" + to_string(imode); + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + histnames_pos_.at(imode).at(ipdg) = "h2dEdx_p_" + particlename + "_pos_" + to_string(imode); + histnames_neg_.at(imode).at(ipdg) = "h2dEdx_p_" + particlename + "_neg_" + to_string(imode); + } + } + + histnames_hits_all_pos_ = "h2dEdx_p_pos"; + histnames_hits_all_neg_ = "h2dEdx_p_neg"; + for (int ipdg = 0; ipdg < NumberOfPidsTrd - 1; ipdg++) { + TString particlename = pid_codes_trd_.at(ipdg).second; + histnames_hits_pos_.at(ipdg) = "h2dEdx_p_" + particlename + "_pos"; + histnames_hits_neg_.at(ipdg) = "h2dEdx_p_" + particlename + "_neg"; + } +} + +void PidTrdRunGetter::Init() { + + inFile_ = new TFile(mcfile_name_, "READ"); + if (!inFile_ || !inFile_->IsOpen()) { + throw std::runtime_error("Could not open input file: " + mcfile_name_); + } + + OpenMcHistograms(); + + if (write_mchistograms_out_ == kTRUE) { + TString name = mcfile_name_; + name.Replace(name.Index(".root"), 5, ""); + TString outfilename = Form("%s_probabilities.root", name.Data()); + outFile_ = new TFile(outfilename, "RECREATE"); + InitMcHistogramsOut(); + } +} + +void PidTrdRunGetter::Exec() { + + for (int inhits = 0; inhits < NumberOfTrdLayers; inhits++) + CalculateProbabilities(inhits); + + std::unique_ptr outFile_getter{TFile::Open(getter_file_, "recreate")}; + outFile_getter->WriteObject(&getter_trd_, getter_name_); + outFile_getter->Close(); + + if (write_mchistograms_out_ == kTRUE) outFile_->Close(); +} + +void PidTrdRunGetter::Finish() { + inFile_->Close(); + if (write_mchistograms_out_ == kTRUE) outFile_->Close(); +} diff --git a/interface/PidTrdRunGetter.hpp b/interface/PidTrdRunGetter.hpp new file mode 100644 index 0000000..5da4812 --- /dev/null +++ b/interface/PidTrdRunGetter.hpp @@ -0,0 +1,60 @@ +#ifndef PIDTRD_RUNGETTER_HPP_ +#define PIDTRD_RUNGETTER_HPP_ + +#include "TH2F.h" + +#include "GetterTrd.h" +#include "ParticleProb.h" + +using std::array; + +class PidTrdRunGetter { + + public: + PidTrdRunGetter(const std::string& mcfile_name, const std::string& getter_file, TString getter_name) { + mcfile_name_ = mcfile_name; + getter_file_ = getter_file; + getter_name_ = getter_name; + }; + virtual ~PidTrdRunGetter() = default; + + void Init(); + void OpenMcHistograms(); + void InitMcHistogramsOut(); + void Exec(); + void Finish(); + void SetWriteMcHistogramsOut(const Bool_t write_mchistograms_out) { write_mchistograms_out_ = write_mchistograms_out; };// Probabilities are written and saved in histograms in addition to Getter + + protected: + void CalculateProbabilities(int trunc_mode); + void CalculateProbabilitiesTot(TH2F* h2all, TH2F* h2tmp, TH2F*& h2prob);// total probability - probability based on particle multiplicites + void CalculateProbabilitiesLike(TH2F* h2tmp, TH2F*& h2prob); // likelihood - probability based on dEdx-distribution of particle + + TString mcfile_name_{""}; + TFile* inFile_; + TFile* outFile_; + TString getter_file_; + TString getter_name_; + + // MC histograms + + array probname_ = {"probT", "probL"}; + + array histnames_all_pos_, histnames_all_neg_; + array, NumberOfTrdLayers> histnames_pos_, histnames_neg_; + array, NumberOfTrdLayers> histnames_prob_pos_, histnames_prob_neg_; + + TString histnames_hits_all_pos_, histnames_hits_all_neg_; + array histnames_hits_pos_, histnames_hits_neg_; + + Int_t nbins_mom_; + Double_t bins_mom_[NbinsMax]; + Int_t nbins_dEdx_; + Double_t bins_dEdx_[NbinsMax]; + + Bool_t write_mchistograms_out_{kFALSE}; + + PidTrd::GetterTrd getter_trd_{}; +}; + +#endif//PIDTRD_RUNGETTER_HPP_ diff --git a/macro/RunGetter.C b/macro/RunGetter.C deleted file mode 100644 index e697f63..0000000 --- a/macro/RunGetter.C +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef __CLING__ - -#include "Getter.h" -#include -#include -#include -#include -#include -#include - -void RunGetter(TString InputFile); - -int main(int argc, char** argv) { - if (argc == 2) RunGetter(argv[1]); - else - RunGetter("pid_getter.root"); - return 1; -} - -#endif// __CLING__ - -void RunGetter(TString InputFile = "pid_getter.root") { - - TFile* f2{TFile::Open(InputFile)}; - std::unique_ptr getter{(Pid::Getter*) f2->Get("pid_getter")}; - std::unique_ptr r{new TRandom}; - -// std::vector fittedPids = {PidParticles::kBgPos, PidParticles::kPionPos, PidParticles::kKaonPos, PidParticles::kProton, PidParticles::kHe3, PidParticles::kDeutron, PidParticles::kBgNeg, PidParticles::kPionNeg, PidParticles::kKaonNeg}; - std::vector fittedPids = {PidParticles::kBgPos, PidParticles::kPionPos, PidParticles::kKaonPos, PidParticles::kProton, PidParticles::kBgNeg, PidParticles::kPionNeg, PidParticles::kKaonNeg}; - - std::map h2Purity; - for (auto pid : fittedPids) { - h2Purity.emplace(pid, new TH2F(Form("h2Purity_%d", pid), Form("Purity for %d;q*p (GeV/#it{c});m^{2} (GeV^{2}/#it{c}^{4})", pid), 1000, -12, 20, 1000, -6., 6)); - } - TAxis* xAxis = h2Purity.at(PidParticles::kBgPos)->GetXaxis(); - TAxis* yAxis = h2Purity.at(PidParticles::kBgPos)->GetYaxis(); - uint nBinsX = xAxis->GetNbins(); - uint nBinsY = yAxis->GetNbins(); - for (uint binX = 1; binX <= nBinsX; binX++) { - float qp = xAxis->GetBinCenter(binX); - for (uint binY = 1; binY <= nBinsY; binY++) { - float m2 = yAxis->GetBinCenter(binY); - auto prob = getter->GetBayesianProbability(qp, m2); - for (auto pid : fittedPids) - h2Purity.at(pid)->Fill(qp, m2, prob.at(pid)); - } - } - TCanvas* c = new TCanvas(); - c->Divide(3, 3); - int i = 0; - for (auto pid : fittedPids) { - c->cd(++i); - gPad->SetLogz(); - h2Purity.at(pid)->Draw("colz"); - h2Purity.at(pid)->GetYaxis()->SetTitleOffset(1.4); - h2Purity.at(pid)->SetMinimum(0.5); - } - /* - for (int i=0; i<10; ++i) - { - const float m2 = r->Uniform(-0.1, 1); - const float qp = r->Uniform(-5, 5); - - auto prob = getter->GetBayesianProbability(qp, m2); - - std::cout << "m2 = " << m2 << " qp = " << qp << " proton probability = " << prob[PidParticles::kProton] << std::endl; - std::cout << "m2 = " << m2 << " qp = " << qp << " kaon probability = " << prob[PidParticles::kKaonPos] << std::endl; - std::cout << "m2 = " << m2 << " qp = " << qp << " pion probability = " << prob[PidParticles::kPionPos] << std::endl; - std::cout << "m2 = " << m2 << " qp = " << qp << " bg probability = " << prob[PidParticles::kBg] << std::endl; - } - - f2->Close(); -*/ -} diff --git a/macro/RunGetterTof.C b/macro/RunGetterTof.C new file mode 100644 index 0000000..23c20ba --- /dev/null +++ b/macro/RunGetterTof.C @@ -0,0 +1,77 @@ +#ifndef __CLING__ + +#include "GetterTof.h" +#include +#include +#include +#include +#include +#include + +void RunGetterTof(TString InputFile); + +int main(int argc, char** argv) { + if (argc == 2) RunGetter(argv[1]); + else + RunGetter("pid_getter_tof.root"); + return 1; +} + +#endif// __CLING__ + +void RunGetterTof(TString InputFile = "pid_getter_tof.root") { + + TFile* f2{TFile::Open(InputFile)}; + + std::unique_ptr getter{(Pid::GetterTof*) f2->Get("pid_getter_tof")}; + std::unique_ptr r{new TRandom}; + + // std::vector fittedPids = {PidTofParticles::kBgPos, PidTofParticles::kPionPos, PidTofParticles::kKaonPos, PidTofParticles::kProton, PidTofParticles::kHe3, PidTofParticles::kDeutron, PidTofParticles::kBgNeg, PidTofParticles::kPionNeg, PidTofParticles::kKaonNeg}; + std::vector fittedPids = {PidTofParticles::kBgPos, PidTofParticles::kPionPos, PidTofParticles::kKaonPos, PidTofParticles::kProton, PidTofParticles::kBgNeg, PidTofParticles::kPionNeg, PidTofParticles::kKaonNeg}; + + std::map h2Purity; + for (auto pid : fittedPids) { + h2Purity.emplace(pid, new TH2F(Form("h2Purity_%d", pid), Form("Purity for %d;q*p (GeV/#it{c});m^{2} (GeV^{2}/#it{c}^{4})", pid), 1000, -12, 20, 1000, -6., 6)); + } + TAxis* xAxis = h2Purity.at(PidTofParticles::kBgPos)->GetXaxis(); + TAxis* yAxis = h2Purity.at(PidTofParticles::kBgPos)->GetYaxis(); + uint nBinsX = xAxis->GetNbins(); + uint nBinsY = yAxis->GetNbins(); + for (uint binX = 1; binX <= nBinsX; binX++) { + float qp = xAxis->GetBinCenter(binX); + for (uint binY = 1; binY <= nBinsY; binY++) { + float m2 = yAxis->GetBinCenter(binY); + auto prob = getter->GetBayesianProbability(qp, m2); + for (auto pid : fittedPids) + h2Purity.at(pid)->Fill(qp, m2, prob.at(pid)); + } + } + TCanvas* c = new TCanvas(); + c->Divide(3, 3); + int i = 0; + for (auto pid : fittedPids) { + c->cd(++i); + gPad->SetLogz(); + h2Purity.at(pid)->Draw("colz"); + h2Purity.at(pid)->GetYaxis()->SetTitleOffset(1.4); + h2Purity.at(pid)->SetMinimum(0.5); + } + + /* + for (int i=0; i<10; ++i) + { + const float m2 = r->Uniform(-0.1, 1); + //const float qp = r->Uniform(-5, 5); + const float qp = r->Uniform(0, 5); + + auto prob = getter->GetBayesianProbability(qp, m2); + + std::cout << "m2 = " << m2 << " qp = " << qp << " proton probability = " << prob[PidTofParticles::kProton] << std::endl; + std::cout << "m2 = " << m2 << " qp = " << qp << " kaon probability = " << prob[PidTofParticles::kKaonPos] << std::endl; + std::cout << "m2 = " << m2 << " qp = " << qp << " pion probability = " << prob[PidTofParticles::kPionPos] << std::endl; + std::cout << "m2 = " << m2 << " qp = " << qp << " bg probability = " << prob[PidTofParticles::kBgPos] << std::endl; + } + + f2->Close(); + */ +} diff --git a/macro/RunGetterTrd.C b/macro/RunGetterTrd.C new file mode 100644 index 0000000..7814503 --- /dev/null +++ b/macro/RunGetterTrd.C @@ -0,0 +1,159 @@ +#ifndef __CLING__ + +#include "GetterTrd.h" +#include +#include +#include +#include +#include +#include + +void RunGetterTrd(TString InputFile); + +int main(int argc, char** argv) { + if (argc == 2) RunGetter(argv[1]); + else + RunGetter("pid_getter_trd.root"); + return 1; +} + +#endif// __CLING__ + +void RunGetterTrd(TString InputFile = "pid_getter_trd.root") { + + TFile* f2{TFile::Open(InputFile)}; + std::unique_ptr getter{(PidTrd::GetterTrd*) f2->Get("pid_getter_trd")}; + + std::unique_ptr r{new TRandom}; + + std::vector Pids = {PidTrdParticles::kBgPos, PidTrdParticles::kPionPos, PidTrdParticles::kKaonPos, PidTrdParticles::kProton, PidTrdParticles::kDeutron, PidTrdParticles::kTriton, PidTrdParticles::kHe3}; + + std::map h2Purity; + for (auto pid : Pids) { + h2Purity.emplace(pid, new TH2F(Form("h2Purity_%d", pid), Form("Purity for %d;p (GeV/#it{c});dE/dx (keV/cm)", pid), 1000, 0., 100, 1000, 0, 200)); + } + TAxis* xAxis = h2Purity.at(PidTrdParticles::kBgPos)->GetXaxis(); + TAxis* yAxis = h2Purity.at(PidTrdParticles::kBgPos)->GetYaxis(); + uint nBinsX = xAxis->GetNbins(); + uint nBinsY = yAxis->GetNbins(); + int nhits = 2; + int charge = 1; + for (uint binX = 1; binX <= nBinsX; binX++) { + float mom = xAxis->GetBinCenter(binX); + for (uint binY = 1; binY <= nBinsY; binY++) { + float dEdx = yAxis->GetBinCenter(binY); + auto prob = getter->GetTrdProbabilitiesTotalAverage(mom, dEdx, charge, nhits); + for (auto pid : Pids) + h2Purity.at(pid)->Fill(mom, dEdx, prob.at(pid)); + } + } + TCanvas* c = new TCanvas(); + c->Divide(3, 3); + int i = 0; + for (auto pid : Pids) { + c->cd(++i); + gPad->SetLogz(); + h2Purity.at(pid)->Draw("colz"); + h2Purity.at(pid)->GetYaxis()->SetTitleOffset(1.4); + //h2Purity.at(pid)->SetMinimum(0.5); + } + + // Calcluation of probabilities for test tracks + + for (int i = 0; i < 10; ++i)// Create test tracks + { + const float pT = 0.0; + const float mom = r->Uniform(0.0, 30.0); + const int charge = 1; + + std::array dEdx_hits = {0.0, 0.0, 0.0, 0.0}; + float dEdx = 0.0; + int nhits_trd = 0; + for (int i = 0; i < NumberOfTrdLayers; i++) { + dEdx_hits.at(i) = r->Uniform(1.5, 40); + if (dEdx_hits.at(i) > 0) { + dEdx += dEdx_hits.at(i); + nhits_trd++; + } + } + + dEdx /= nhits_trd; + + std::cout << std::endl; + std::cout << "-- Track " << i << " -- " << std::endl; + std::cout << "momentum = " << mom << std::endl; + std::cout << "dEdx (hits) = " << dEdx_hits.at(0) << ", " << dEdx_hits.at(1) << ", " << dEdx_hits.at(2) << ", " << dEdx_hits.at(3) << std::endl; + std::cout << "dEdx (average) = " << dEdx << std::endl; + std::cout << std::endl; + + // Simple functions for "average" mode (no truncation) + + float min_purity = 0.2;// minimum required purity for pid hypothesis + + std::cout << "- Total probability method -" << std::endl; + + auto probT = getter->GetTrdProbabilitiesTotalAverage(mom, dEdx, charge, nhits_trd); + + for (const auto& pdg : pid_codes_trd_) + printf("mom = %.4f dEdx (average) = %.4f probability %-4s %-10d = %.4f\n", mom, dEdx, pdg.second.Data(), pdg.first, probT[pdg.first]); + std::cout << std::endl; + + int pid_hypothesis = getter->GetTrdPid(probT, min_purity, charge); + + std::cout << "pid hypothesis : " << pid_hypothesis << std::endl; + std::cout << std::endl; + + std::cout << "- Likelihood method -" << std::endl; + + auto probL = getter->GetTrdProbabilitiesLikelihoodAverage(mom, dEdx_hits, charge); + + for (const auto& pdg : pid_codes_trd_) + printf("mom = %.4f dEdx (hits) = %.4f, %.4f, %.4f, %.4f probability %-4s %-10d = %.4f\n", mom, dEdx_hits.at(0), dEdx_hits.at(1), dEdx_hits.at(2), dEdx_hits.at(3), pdg.second.Data(), pdg.first, probL[pdg.first]); + std::cout << std::endl; + + pid_hypothesis = getter->GetTrdPid(probL, min_purity, charge); + + std::cout << "pid hypothesis : " << pid_hypothesis << std::endl; + std::cout << std::endl; + + // Alternative, more flexible way to run getter: + // Functions with options for truncation mode + /* + float min_purity = 0.2; // minimum required purity for pid hypothesis + + //Options: + //getter->SetTruncationMode(0); // Default truncation mode is "average", optional: set truncation mode (see README or GetterTrd.h) + //getter->SetMinHits(1); // Default minimum number of required hits per track is 1 + + TrdContainer trdtrack(mom, pT, charge, nhits_trd, dEdx_hits); + + std::cout << "- Total probability method -" << std::endl; + + getter->SetProbabiltyMode(0); + auto probT = getter->GetTrdProbabilities(trdtrack); + + for (const auto& pdg : pid_codes_trd_) + printf("mom = %.4f dEdx (average) = %.4f probability %-4s %-10d = %.4f\n", mom, dEdx, pdg.second.Data(), pdg.first, probT[pdg.first]); + std::cout << std::endl; + + int pid_hypothesis = getter->GetTrdPid(probT, min_purity, charge); + + std::cout << "pid hypothesis : " << pid_hypothesis << std::endl; + std::cout << std::endl; + + std::cout << "- Likelihood method -" << std::endl; + + getter->SetProbabiltyMode(1); + auto probL = getter->GetTrdProbabilitiesMulti(trdtrack); + + for (const auto& pdg : pid_codes_trd_) + printf("mom = %.4f dEdx (hits) = %.4f, %.4f, %.4f, %.4f probability %-4s %-10d = %.4f\n", mom, dEdx_hits.at(0), dEdx_hits.at(1), dEdx_hits.at(2), dEdx_hits.at(3), pdg.second.Data(), pdg.first, probL[pdg.first]); + std::cout << std::endl; + + pid_hypothesis = getter->GetTrdPid(probL, min_purity, charge); + + std::cout << "pid hypothesis : " << pid_hypothesis << std::endl; + std::cout << std::endl; + */ + } +} diff --git a/macro/fit_botvina12.C b/macro/fit_botvina12.C index f38f259..a3502e5 100644 --- a/macro/fit_botvina12.C +++ b/macro/fit_botvina12.C @@ -47,7 +47,7 @@ void fit_botvina12() { std::unique_ptr fCuts{TFile::Open(cutsFileName, "read")}; Pid::Fitter tof; - Pid::Getter getter; + Pid::GetterTof getter; float xmin, xmax, ymin, ymax; cout << "\n\npionpos\n"; @@ -65,13 +65,13 @@ void fit_botvina12() { TF1 pionpos_1("pionpos_1", "pol9", xmin, xmax); TF1 pionpos_2("pionpos_2", "pol9", xmin, xmax); - Pid::ParticleFit pionpos(PidParticles::kPionPos); + Pid::ParticleFit pionpos(PidTofParticles::kPionPos); pionpos.SetParametrization({pionpos_0, pionpos_1, pionpos_2}); pionpos.SetFitFunction(fit_pionpos); pionpos.SetRange(xmin, xmax); pionpos.SetIsFitted(); - tof.AddParticle(pionpos, PidParticles::kPionPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hpionpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -95,13 +95,13 @@ void fit_botvina12() { TF1 kaonpos_1("kaonpos_1", "pol9", xmin, xmax); TF1 kaonpos_2("kaonpos_2", "pol8", xmin, xmax); - Pid::ParticleFit kaonpos(PidParticles::kKaonPos); + Pid::ParticleFit kaonpos(PidTofParticles::kKaonPos); kaonpos.SetParametrization({kaonpos_0, kaonpos_1, kaonpos_2}); kaonpos.SetFitFunction(fit_kaonpos); kaonpos.SetRange(xmin, xmax); kaonpos.SetIsFitted(); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); tof.SetHisto2D(std::move(hkaonpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -125,13 +125,13 @@ void fit_botvina12() { TF1 proton_1("proton_1", "pol9", 1.5, 19.); TF1 proton_2("proton_2", "pol12", 1., 20.); - Pid::ParticleFit proton(PidParticles::kProton); + Pid::ParticleFit proton(PidTofParticles::kProton); proton.SetParametrization({proton_0, proton_1, proton_2}); proton.SetFitFunction(fit_proton); proton.SetRange(xmin, xmax); proton.SetIsFitted(); - tof.AddParticle(proton, PidParticles::kPionPos); + tof.AddParticle(proton, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hproton_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -146,7 +146,7 @@ void fit_botvina12() { hhe3_cut->Rebin2D(10, 5); xmin = 1.8, xmax = 8.7, ymin = 0.5, ymax = 3.5; tof.SetChi2Max(10); - Pid::ParticleFit he3(PidParticles::kHe3); + Pid::ParticleFit he3(PidTofParticles::kHe3); TF1 he3_fit("fit_he3", "gaus", ymin, ymax); he3_fit.SetParNames("p9", "p10", "p11"); he3_fit.SetParLimits(0, 0., 5.e2); @@ -160,7 +160,7 @@ void fit_botvina12() { he3.SetFitFunction(he3_fit); he3.SetRange(xmin, xmax); he3.SetIsFitted(); - tof.AddParticle(he3, PidParticles::kHe3); + tof.AddParticle(he3, PidTofParticles::kHe3); tof.SetHisto2D(std::move(hhe3_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -184,13 +184,13 @@ void fit_botvina12() { TF1 deutron_1("deutron_1", "pol6", xmin, xmax); TF1 deutron_2("deutron_2", "pol10", xmin, xmax); - Pid::ParticleFit deutron(PidParticles::kProton); + Pid::ParticleFit deutron(PidTofParticles::kProton); deutron.SetParametrization({deutron_0, deutron_1, deutron_2}); deutron.SetFitFunction(fit_deutron); deutron.SetRange(xmin, xmax); deutron.SetIsFitted(); - tof.AddParticle(deutron, PidParticles::kDeutron); + tof.AddParticle(deutron, PidTofParticles::kDeutron); tof.SetHisto2D(std::move(hdeutron_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -201,7 +201,7 @@ void fit_botvina12() { cout << "\n\nbgpos\n"; xmin = 0.3, xmax = 20., ymin = -3., ymax = 3.; - Pid::ParticleFit bgpos(PidParticles::kBgPos); + Pid::ParticleFit bgpos(PidTofParticles::kBgPos); TF1 bgpos_fit("fit_bgpos", "pol2", ymin, ymax); bgpos_fit.SetParNames("p15", "p16", "p17"); TF1 bgpos_0("bgpos_0", "pol3", xmin, xmax); @@ -223,12 +223,12 @@ void fit_botvina12() { proton.SetIsFixed({false, true, true}); deutron.SetIsFixed({false, true, true}); he3.SetIsFixed({false, true, true}); - tof.AddParticle(pionpos, PidParticles::kPionPos); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); - tof.AddParticle(proton, PidParticles::kProton); - tof.AddParticle(deutron, PidParticles::kDeutron); - tof.AddParticle(he3, PidParticles::kHe3); - tof.AddParticle(bgpos, PidParticles::kBgPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); + tof.AddParticle(proton, PidTofParticles::kProton); + tof.AddParticle(deutron, PidTofParticles::kDeutron); + tof.AddParticle(he3, PidTofParticles::kHe3); + tof.AddParticle(bgpos, PidTofParticles::kBgPos); tof.SetHisto2D(std::move(hpos)); tof.SetRangeX(xmin, xmax); @@ -236,12 +236,12 @@ void fit_botvina12() { tof.SetOutputFileName("allpos.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionPos), PidParticles::kPionPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonPos), PidParticles::kKaonPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kProton), PidParticles::kProton); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kHe3), PidParticles::kHe3); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kDeutron), PidParticles::kDeutron); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgPos), PidParticles::kBgPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionPos), PidTofParticles::kPionPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonPos), PidTofParticles::kKaonPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kProton), PidTofParticles::kProton); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kHe3), PidTofParticles::kHe3); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kDeutron), PidTofParticles::kDeutron); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgPos), PidTofParticles::kBgPos); tof.Clear(); @@ -260,13 +260,13 @@ void fit_botvina12() { TF1 pionneg_1("pionneg_1", "pol14", -11.7, xmax); TF1 pionneg_2("pionneg_2", "pol12", -11.7, xmax); - Pid::ParticleFit pionneg(PidParticles::kPionNeg); + Pid::ParticleFit pionneg(PidTofParticles::kPionNeg); pionneg.SetParametrization({pionneg_0, pionneg_1, pionneg_2}); pionneg.SetFitFunction(fit_pionneg); pionneg.SetRange(xmin, xmax); pionneg.SetIsFitted(); - tof.AddParticle(pionneg, PidParticles::kPionNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); tof.SetHisto2D(std::move(hpionneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -290,13 +290,13 @@ void fit_botvina12() { TF1 kaonneg_1("kaonneg_1", "pol11", xmin, xmax); TF1 kaonneg_2("kaonneg_2", "pol8", xmin, xmax); - Pid::ParticleFit kaonneg(PidParticles::kKaonNeg); + Pid::ParticleFit kaonneg(PidTofParticles::kKaonNeg); kaonneg.SetParametrization({kaonneg_0, kaonneg_1, kaonneg_2}); kaonneg.SetFitFunction(fit_kaonneg); kaonneg.SetRange(xmin, xmax); kaonneg.SetIsFitted(); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); tof.SetHisto2D(std::move(hkaonneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -307,7 +307,7 @@ void fit_botvina12() { cout << "\n\nbgneg\n"; xmin = -10., xmax = -0.25, ymin = -1., ymax = 2.; - Pid::ParticleFit bgneg(PidParticles::kBgNeg); + Pid::ParticleFit bgneg(PidTofParticles::kBgNeg); TF1 bgneg_0("bgneg_0", "pol3", xmin, xmax);//bgneg_0.SetParameters(100, 0, 0); TF1 bgneg_1("bgneg_1", "pol5", xmin, xmax);//bgneg_1.SetParameters(0, 0, 0); TF1 bgneg_2("bgneg_2", "pol5", xmin, xmax);//bgneg_2.SetParameters(0.0, 0.0, 0); @@ -330,9 +330,9 @@ void fit_botvina12() { tof.SetChi2Max(1e6); pionneg.SetIsFixed({false, true, true}); kaonneg.SetIsFixed({false, true, true}); - tof.AddParticle(pionneg, PidParticles::kPionNeg); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); - tof.AddParticle(bgneg, PidParticles::kBgNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); + tof.AddParticle(bgneg, PidTofParticles::kBgNeg); tof.SetHisto2D(std::move(hneg)); tof.SetRangeX(xmin, xmax); @@ -340,9 +340,9 @@ void fit_botvina12() { tof.SetOutputFileName("allneg.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionNeg), PidParticles::kPionNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonNeg), PidParticles::kKaonNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgNeg), PidParticles::kBgNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionNeg), PidTofParticles::kPionNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonNeg), PidTofParticles::kKaonNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgNeg), PidTofParticles::kBgNeg); tof.Clear(); diff --git a/macro/fit_botvina3.C b/macro/fit_botvina3.C index 22d71ea..d2245f1 100644 --- a/macro/fit_botvina3.C +++ b/macro/fit_botvina3.C @@ -47,7 +47,7 @@ void fit_botvina3() { std::unique_ptr fCuts{TFile::Open(cutsFileName, "read")}; Pid::Fitter tof; - Pid::Getter getter; + Pid::GetterTof getter; float xmin, xmax, ymin, ymax; cout << "\n\npionpos\n"; @@ -65,13 +65,13 @@ void fit_botvina3() { TF1 pionpos_1("pionpos_1", "pol2", xmin, xmax); TF1 pionpos_2("pionpos_2", "pol5", xmin, xmax); - Pid::ParticleFit pionpos(PidParticles::kPionPos); + Pid::ParticleFit pionpos(PidTofParticles::kPionPos); pionpos.SetParametrization({pionpos_0, pionpos_1, pionpos_2}); pionpos.SetFitFunction(fit_pionpos); pionpos.SetRange(xmin, xmax); pionpos.SetIsFitted(); - tof.AddParticle(pionpos, PidParticles::kPionPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hpionpos_cut)); tof.SetRangeX(xmin, 8.); tof.SetRangeY(ymin, ymax); @@ -95,13 +95,13 @@ void fit_botvina3() { TF1 kaonpos_1("kaonpos_1", "pol2", xmin, xmax); TF1 kaonpos_2("kaonpos_2", "pol6", xmin, xmax); - Pid::ParticleFit kaonpos(PidParticles::kKaonPos); + Pid::ParticleFit kaonpos(PidTofParticles::kKaonPos); kaonpos.SetParametrization({kaonpos_0, kaonpos_1, kaonpos_2}); kaonpos.SetFitFunction(fit_kaonpos); kaonpos.SetRange(xmin, xmax); kaonpos.SetIsFitted(); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); tof.SetHisto2D(std::move(hkaonpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -125,13 +125,13 @@ void fit_botvina3() { TF1 proton_1("proton_1", "pol11", .3, 12.); TF1 proton_2("proton_2", "pol8", .3, 12.); - Pid::ParticleFit proton(PidParticles::kProton); + Pid::ParticleFit proton(PidTofParticles::kProton); proton.SetParametrization({proton_0, proton_1, proton_2}); proton.SetFitFunction(fit_proton); proton.SetRange(xmin, xmax); proton.SetIsFitted(); - tof.AddParticle(proton, PidParticles::kPionPos); + tof.AddParticle(proton, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hproton_cut)); tof.SetRangeX(xmin, 8.); tof.SetRangeY(ymin, ymax); @@ -146,7 +146,7 @@ void fit_botvina3() { hhe3_cut->Rebin2D(10, 10); xmin = 1.5, xmax = 7., ymin = 1., ymax = 3.; tof.SetChi2Max(10); - Pid::ParticleFit he3(PidParticles::kHe3); + Pid::ParticleFit he3(PidTofParticles::kHe3); TF1 he3_fit("fit_he3", "gaus", ymin, ymax); he3_fit.SetParNames("p9", "p10", "p11"); he3_fit.SetParLimits(0, 0., 1.e6); @@ -160,7 +160,7 @@ void fit_botvina3() { he3.SetFitFunction(he3_fit); he3.SetRange(xmin, xmax); he3.SetIsFitted(); - tof.AddParticle(he3, PidParticles::kHe3); + tof.AddParticle(he3, PidTofParticles::kHe3); tof.SetHisto2D(std::move(hhe3_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -184,13 +184,13 @@ void fit_botvina3() { TF1 deutron_1("deutron_1", "pol10", xmin, xmax); TF1 deutron_2("deutron_2", "pol3", xmin, xmax); - Pid::ParticleFit deutron(PidParticles::kProton); + Pid::ParticleFit deutron(PidTofParticles::kProton); deutron.SetParametrization({deutron_0, deutron_1, deutron_2}); deutron.SetFitFunction(fit_deutron); deutron.SetRange(xmin, xmax); deutron.SetIsFitted(); - tof.AddParticle(deutron, PidParticles::kDeutron); + tof.AddParticle(deutron, PidTofParticles::kDeutron); tof.SetHisto2D(std::move(hdeutron_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -201,7 +201,7 @@ void fit_botvina3() { cout << "\n\nbgpos\n"; xmin = 0.3, xmax = 20., ymin = -3., ymax = 3.; - Pid::ParticleFit bgpos(PidParticles::kBgPos); + Pid::ParticleFit bgpos(PidTofParticles::kBgPos); TF1 bgpos_fit("fit_bgpos", "pol2", ymin, ymax); bgpos_fit.SetParNames("p15", "p16", "p17"); TF1 bgpos_0("bgpos_0", "pol3", xmin, xmax); @@ -223,12 +223,12 @@ void fit_botvina3() { proton.SetIsFixed({false, true, true}); deutron.SetIsFixed({false, true, true}); he3.SetIsFixed({false, true, true}); - tof.AddParticle(pionpos, PidParticles::kPionPos); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); - tof.AddParticle(proton, PidParticles::kProton); - tof.AddParticle(deutron, PidParticles::kDeutron); - tof.AddParticle(he3, PidParticles::kHe3); - tof.AddParticle(bgpos, PidParticles::kBgPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); + tof.AddParticle(proton, PidTofParticles::kProton); + tof.AddParticle(deutron, PidTofParticles::kDeutron); + tof.AddParticle(he3, PidTofParticles::kHe3); + tof.AddParticle(bgpos, PidTofParticles::kBgPos); tof.SetHisto2D(std::move(hpos)); tof.SetRangeX(xmin, xmax); @@ -236,12 +236,12 @@ void fit_botvina3() { tof.SetOutputFileName("allpos.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionPos), PidParticles::kPionPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonPos), PidParticles::kKaonPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kProton), PidParticles::kProton); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kHe3), PidParticles::kHe3); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kDeutron), PidParticles::kDeutron); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgPos), PidParticles::kBgPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionPos), PidTofParticles::kPionPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonPos), PidTofParticles::kKaonPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kProton), PidTofParticles::kProton); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kHe3), PidTofParticles::kHe3); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kDeutron), PidTofParticles::kDeutron); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgPos), PidTofParticles::kBgPos); tof.Clear(); @@ -260,13 +260,13 @@ void fit_botvina3() { TF1 pionneg_1("pionneg_1", "pol2", xmin, xmax); TF1 pionneg_2("pionneg_2", "pol5", xmin, xmax); - Pid::ParticleFit pionneg(PidParticles::kPionNeg); + Pid::ParticleFit pionneg(PidTofParticles::kPionNeg); pionneg.SetParametrization({pionneg_0, pionneg_1, pionneg_2}); pionneg.SetFitFunction(fit_pionneg); pionneg.SetRange(xmin, xmax); pionneg.SetIsFitted(); - tof.AddParticle(pionneg, PidParticles::kPionNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); tof.SetHisto2D(std::move(hpionneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -290,13 +290,13 @@ void fit_botvina3() { TF1 kaonneg_1("kaonneg_1", "pol2", xmin, xmax); TF1 kaonneg_2("kaonneg_2", "pol2", xmin, xmax); - Pid::ParticleFit kaonneg(PidParticles::kKaonNeg); + Pid::ParticleFit kaonneg(PidTofParticles::kKaonNeg); kaonneg.SetParametrization({kaonneg_0, kaonneg_1, kaonneg_2}); kaonneg.SetFitFunction(fit_kaonneg); kaonneg.SetRange(xmin, xmax); kaonneg.SetIsFitted(); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); tof.SetHisto2D(std::move(hkaonneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -307,7 +307,7 @@ void fit_botvina3() { cout << "\n\nbgneg\n"; xmin = -10., xmax = -0.25, ymin = -1., ymax = 2.; - Pid::ParticleFit bgneg(PidParticles::kBgNeg); + Pid::ParticleFit bgneg(PidTofParticles::kBgNeg); TF1 bgneg_0("bgneg_0", "pol3", xmin, xmax);//bgneg_0.SetParameters(100, 0, 0); TF1 bgneg_1("bgneg_1", "pol5", xmin, xmax);//bgneg_1.SetParameters(0, 0, 0); TF1 bgneg_2("bgneg_2", "pol5", xmin, xmax);//bgneg_2.SetParameters(0.0, 0.0, 0); @@ -330,9 +330,9 @@ void fit_botvina3() { tof.SetChi2Max(1e6); pionneg.SetIsFixed({false, true, true}); kaonneg.SetIsFixed({false, true, true}); - tof.AddParticle(pionneg, PidParticles::kPionNeg); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); - tof.AddParticle(bgneg, PidParticles::kBgNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); + tof.AddParticle(bgneg, PidTofParticles::kBgNeg); tof.SetHisto2D(std::move(hneg)); tof.SetRangeX(xmin, xmax); @@ -340,9 +340,9 @@ void fit_botvina3() { tof.SetOutputFileName("allneg.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionNeg), PidParticles::kPionNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonNeg), PidParticles::kKaonNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgNeg), PidParticles::kBgNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionNeg), PidTofParticles::kPionNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonNeg), PidTofParticles::kKaonNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgNeg), PidTofParticles::kBgNeg); tof.Clear(); diff --git a/pid_new/examples/shine/FitBinnedTracks.cpp b/pid_new/examples/shine/FitBinnedTracks.cpp index 9b9bab6..032b216 100644 --- a/pid_new/examples/shine/FitBinnedTracks.cpp +++ b/pid_new/examples/shine/FitBinnedTracks.cpp @@ -4,219 +4,218 @@ #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include +#include +#include +#include #include #include -#include -#include #include -#include +#include #include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include +#include #include #include -#include using namespace std; -int main(int argc, char **argv) { +int main(int argc, char** argv) { - string input_file{argv[1]}; - string charge_str = "negative"; - string bin_name{argv[2]}; + string input_file{argv[1]}; + string charge_str = "negative"; + string bin_name{argv[2]}; - std::vector ncl_binning{30., 40., 50., 60., 70., 90., 110., 130., 150., 170., 200., 250.}; - TAxis ncl_axis(ncl_binning.size() - 1, ncl_binning.data()); + std::vector ncl_binning{30., 40., 50., 60., 70., 90., 110., 130., 150., 170., 200., 250.}; + TAxis ncl_axis(ncl_binning.size() - 1, ncl_binning.data()); - size_t min_entries_ncl = 1000ul; - std::map nn_ncl_bin; + size_t min_entries_ncl = 1000ul; + std::map nn_ncl_bin; - ROOT::RDataFrame rdf(charge_str + "/" + bin_name + "/track_tree", input_file); + ROOT::RDataFrame rdf(charge_str + "/" + bin_name + "/track_tree", input_file); - auto rdf_ncl_binned = rdf - .Define("ncl_bin", [ncl_axis](int ncl) -> int { return ncl_axis.FindBin(ncl); }, {"ncl"}) - .Define("sigma_factor", "TMath::Power(ncl, -0.5)"); - rdf_ncl_binned.Display(".*").GetValue().Print(); - std::cout << "NE = " << rdf_ncl_binned.Count().GetValue() << std::endl; + auto rdf_ncl_binned = rdf + .Define("ncl_bin", [ncl_axis](int ncl) -> int { return ncl_axis.FindBin(ncl); }, {"ncl"}) + .Define("sigma_factor", "TMath::Power(ncl, -0.5)"); + rdf_ncl_binned.Display(".*").GetValue().Print(); + std::cout << "NE = " << rdf_ncl_binned.Count().GetValue() << std::endl; - rdf_ncl_binned.Foreach([&nn_ncl_bin](int ncl_bin) { - auto insert_result = nn_ncl_bin.emplace(ncl_bin, 0); - if (!insert_result.second) { /* map entry already exists */ - insert_result.first->second++; - } - }, {"ncl_bin"}); + rdf_ncl_binned.Foreach([&nn_ncl_bin](int ncl_bin) { + auto insert_result = nn_ncl_bin.emplace(ncl_bin, 0); + if (!insert_result.second) { /* map entry already exists */ + insert_result.first->second++; + } + }, + {"ncl_bin"}); + + auto rdf_filtered = rdf_ncl_binned.Filter([=](int ncl_bin) { + return ncl_bin != -1 && nn_ncl_bin.at(ncl_bin) > min_entries_ncl; + }, + {"ncl_bin"}, "nn_cut"); + + rdf_filtered.Display(".*").GetValue().Print(); + std::cout << "NE = " << rdf_filtered.Count().GetValue() << std::endl; + + /* make snapshot */ + std::string tmp_name = std::tmpnam(nullptr); + tmp_name.append(".root"); + std::cout << tmp_name << std::endl; + rdf_filtered.Snapshot("track_tree", tmp_name); + + /* prepare model */ + auto w = RooWorkspace("w"); + + /* observables */ + RooRealVar v_dedx("dedx", "dE/dx", 0, 2.5, "MIP"); + RooCategory v_ncl_bin("ncl_bin", ""); + for (auto p : nn_ncl_bin) { + if (p.second >= min_entries_ncl) v_ncl_bin.defineType(Form("ncl_bin_%d", p.first), p.first); + } + + /* input data */ + RooDataSet ds("ds", "", RooArgSet(v_dedx, v_ncl_bin), RooFit::ImportFromFile(tmp_name.c_str(), "track_tree")); + v_ncl_bin.Print(); + ds.Print(); + + /* model parameters */ + RooRealVar p_delta("delta", "#delta", 0.05); + RooRealVar p_sigma0_pion_neg("sigma0_pion_neg", "#sigma_0 (#pi^{-})", 0., 0.8, "MIP"); + RooRealVar p_yield_pion_neg("yield_pion_neg", "Y (#pi^{-})", 0., 1.); + RooRealVar p_yield_electron("yield_electron", "Y (e)", 0., 1.); + + /* Bethe-Blochs */ + auto p_mean = rdf_filtered.Mean("p").GetValue(); + RooRealVar bb_pion_neg("bb_pion_neg", "BB (#pi^{-})", 0.85 * BetheBlochAntoniMod(p_mean / 0.139), 1.15 * BetheBlochAntoniMod(p_mean / 0.139)); + RooRealVar bb_electron("bb_electron", "BB (e)", 0.85 * BetheBlochAntoniMod(p_mean / 0.000511), 1.15 * BetheBlochAntoniMod(p_mean / 0.000511)); + RooConstVar bb_kaon_neg("bb_kaon_neg", "BB (K^{-})", BetheBlochAntoniMod(p_mean / 0.494)); + + /* Derivatives */ + RooFormulaVar p_sigma0_electron("sigma0_electron", "#sigma_0 (e^{-})", + "@0 * TMath::Power(@1 / @2, 0.65)", + RooArgList(p_sigma0_pion_neg, bb_electron, bb_pion_neg)); + + RooFormulaVar p_sigma0_kaon_neg("sigma0_kaon_neg", "#sigma_0 (K^{-})", + "@0 * TMath::Power(@1 / @2, 0.65)", + RooArgList(p_sigma0_pion_neg, bb_kaon_neg, bb_pion_neg)); + + std::map sim_fit_components; + for (auto p : nn_ncl_bin) { + if (p.second < min_entries_ncl) continue; + + auto ncl_bin = p.first; + auto mean_sigma_factor = rdf_filtered + .Filter([ncl_bin](int _ncl_bin) { return ncl_bin == _ncl_bin; }, {"ncl_bin"}) + .Mean("sigma_factor") + .GetValue(); + + auto sigma_factor = new RooConstVar("sigma_factor", "", mean_sigma_factor); + auto p_sigma_pion_neg = new RooFormulaVar("sigma_pion_neg", "#sigma (#pi^{-})", + "@0 * @1", + RooArgList(p_sigma0_pion_neg, *sigma_factor)); + + auto p_sigma_electron = new RooFormulaVar("sigma_electron", "#sigma (e^{-})", + "@0 * @1", + RooArgList(p_sigma0_electron, *sigma_factor)); + + auto p_sigma_kaon_neg = new RooFormulaVar("sigma_kaon_neg", "#sigma (K^{-})", + "@0 * @1", + RooArgList(p_sigma0_kaon_neg, *sigma_factor)); + + auto pdf_pion_neg = new AsymmetricGaussianPDF(Form("pdf_pion_neg_%d", ncl_bin), "", v_dedx, bb_pion_neg, + *p_sigma_pion_neg, p_delta); + auto pdf_electron = new AsymmetricGaussianPDF(Form("pdf_electron_%d", ncl_bin), "", v_dedx, bb_electron, + *p_sigma_electron, p_delta); + auto pdf_kaon_neg = new AsymmetricGaussianPDF(Form("pdf_kaon_neg_%d", ncl_bin), "", v_dedx, bb_kaon_neg, + *p_sigma_kaon_neg, p_delta); + + auto composite_pdf = new RooAddPdf(Form("pdf_composite_%d", ncl_bin), "", + RooArgList(*pdf_pion_neg, *pdf_electron), + RooArgList(p_yield_pion_neg)); + composite_pdf->Print(); + + sim_fit_components.emplace(std::string(Form("ncl_bin_%d", ncl_bin)), composite_pdf); + } + + if (!sim_fit_components.empty()) { + + auto sim_pdf = new RooSimultaneous("sim_pdf", "", sim_fit_components, v_ncl_bin); + sim_pdf->Print(); + + /** + * FITTING + */ + auto fit_result = sim_pdf->fitTo(ds, RooFit::Minos(), RooFit::Save()); + TFile qa_file(input_file.c_str(), "update"); + string qa_dir_name = charge_str + "/" + bin_name + "/qa"; + qa_file.mkdir(qa_dir_name.c_str(), "dEdx fit qa", true); - auto rdf_filtered = rdf_ncl_binned.Filter([=](int ncl_bin) { - return ncl_bin != -1 && nn_ncl_bin.at(ncl_bin) > min_entries_ncl; - }, {"ncl_bin"}, "nn_cut"); + /** + * Extracting fit parameters + */ + map par_map_results; + auto fit_params = sim_pdf->getParameters(ds); + for (auto par : *fit_params) { + string par_name{par->GetName()}; + + double value = dynamic_cast(par)->getVal(); + double value_err = dynamic_cast(par)->getError(); + par_map_results.emplace(par_name, value); + par_map_results.emplace(par_name + "_err", value_err); + } + auto c = new TCanvas; + c->Print("fit_qa.pdf(", "pdf"); - rdf_filtered.Display(".*").GetValue().Print(); - std::cout << "NE = " << rdf_filtered.Count().GetValue() << std::endl; + qa_file.cd(qa_dir_name.c_str()); + fit_result->Write("fit_result"); - /* make snapshot */ - std::string tmp_name = std::tmpnam(nullptr); - tmp_name.append(".root"); - std::cout << tmp_name << std::endl; - rdf_filtered.Snapshot("track_tree", tmp_name); + gDirectory->WriteObject(&par_map_results, "fit_parameters"); - /* prepare model */ - auto w = RooWorkspace("w"); + for (const auto& st : v_ncl_bin) { + std::cout << st->GetName() << std::endl; + auto frame_dedx = v_dedx.frame(); + ds.plotOn(frame_dedx, RooFit::Cut(std::string("ncl_bin==ncl_bin::").append(st->GetName()).c_str())); + sim_pdf->plotOn(frame_dedx, RooFit::Slice(v_ncl_bin, st->GetName()), RooFit::ProjWData(v_ncl_bin, ds)); - /* observables */ - RooRealVar v_dedx("dedx", "dE/dx", 0, 2.5, "MIP"); - RooCategory v_ncl_bin("ncl_bin", ""); - for (auto p: nn_ncl_bin) { - if (p.second >= min_entries_ncl) v_ncl_bin.defineType(Form("ncl_bin_%d", p.first), p.first); - } + qa_file.cd(qa_dir_name.c_str()); + frame_dedx->Write(Form("fit_%s", st->GetName())); - /* input data */ - RooDataSet ds("ds", "", RooArgSet(v_dedx, v_ncl_bin), RooFit::ImportFromFile(tmp_name.c_str(), "track_tree")); - v_ncl_bin.Print(); - ds.Print(); - - /* model parameters */ - RooRealVar p_delta("delta", "#delta", 0.05); - RooRealVar p_sigma0_pion_neg("sigma0_pion_neg", "#sigma_0 (#pi^{-})", 0., 0.8, "MIP"); - RooRealVar p_yield_pion_neg("yield_pion_neg", "Y (#pi^{-})", 0., 1.); - RooRealVar p_yield_electron("yield_electron", "Y (e)", 0., 1.); - - /* Bethe-Blochs */ - auto p_mean = rdf_filtered.Mean("p").GetValue(); - RooRealVar bb_pion_neg("bb_pion_neg", "BB (#pi^{-})", 0.85*BetheBlochAntoniMod(p_mean / 0.139), 1.15*BetheBlochAntoniMod(p_mean / 0.139)); - RooRealVar bb_electron("bb_electron", "BB (e)", 0.85*BetheBlochAntoniMod(p_mean / 0.000511), 1.15*BetheBlochAntoniMod(p_mean / 0.000511)); - RooConstVar bb_kaon_neg("bb_kaon_neg", "BB (K^{-})", BetheBlochAntoniMod(p_mean / 0.494)); - - /* Derivatives */ - RooFormulaVar p_sigma0_electron("sigma0_electron", "#sigma_0 (e^{-})", - "@0 * TMath::Power(@1 / @2, 0.65)", - RooArgList(p_sigma0_pion_neg, bb_electron, bb_pion_neg)); - - RooFormulaVar p_sigma0_kaon_neg("sigma0_kaon_neg", "#sigma_0 (K^{-})", - "@0 * TMath::Power(@1 / @2, 0.65)", - RooArgList(p_sigma0_pion_neg, bb_kaon_neg, bb_pion_neg)); - - std::map sim_fit_components; - for (auto p : nn_ncl_bin) { - if (p.second < min_entries_ncl) continue; - - auto ncl_bin = p.first; - auto mean_sigma_factor = rdf_filtered - .Filter([ncl_bin](int _ncl_bin) { return ncl_bin == _ncl_bin; }, {"ncl_bin"}) - .Mean("sigma_factor").GetValue(); - - auto sigma_factor = new RooConstVar("sigma_factor", "", mean_sigma_factor); - auto p_sigma_pion_neg = new RooFormulaVar("sigma_pion_neg", "#sigma (#pi^{-})", - "@0 * @1", - RooArgList(p_sigma0_pion_neg, *sigma_factor)); - - auto p_sigma_electron = new RooFormulaVar("sigma_electron", "#sigma (e^{-})", - "@0 * @1", - RooArgList(p_sigma0_electron, *sigma_factor)); - - auto p_sigma_kaon_neg = new RooFormulaVar("sigma_kaon_neg", "#sigma (K^{-})", - "@0 * @1", - RooArgList(p_sigma0_kaon_neg, *sigma_factor)); - - auto pdf_pion_neg = new AsymmetricGaussianPDF(Form("pdf_pion_neg_%d", ncl_bin), "", v_dedx, bb_pion_neg, - *p_sigma_pion_neg, p_delta); - auto pdf_electron = new AsymmetricGaussianPDF(Form("pdf_electron_%d", ncl_bin), "", v_dedx, bb_electron, - *p_sigma_electron, p_delta); - auto pdf_kaon_neg = new AsymmetricGaussianPDF(Form("pdf_kaon_neg_%d", ncl_bin), "", v_dedx, bb_kaon_neg, - *p_sigma_kaon_neg, p_delta); - - auto composite_pdf = new RooAddPdf(Form("pdf_composite_%d", ncl_bin), "", - RooArgList(*pdf_pion_neg, *pdf_electron), - RooArgList(p_yield_pion_neg)); - composite_pdf->Print(); - - sim_fit_components.emplace(std::string(Form("ncl_bin_%d", ncl_bin)), composite_pdf); + c->Clear(); + c->Divide(2, 1); + c->cd(1); + frame_dedx->Draw(); + c->cd(2); + gPad->SetLogy(); + frame_dedx->Draw(); + c->Print("fit_qa.pdf", "pdf"); } - if (!sim_fit_components.empty()) { + auto frame_dedx = v_dedx.frame(); + ds.plotOn(frame_dedx); + sim_pdf->plotOn(frame_dedx, RooFit::ProjWData(v_ncl_bin, ds)); + c->Clear(); + frame_dedx->Draw(); + c->Print("fit_qa.pdf", "pdf"); - auto sim_pdf = new RooSimultaneous("sim_pdf", "", sim_fit_components, v_ncl_bin); - sim_pdf->Print(); - - /** - * FITTING - */ - auto fit_result = sim_pdf->fitTo(ds, RooFit::Minos(), RooFit::Save()); - - TFile qa_file(input_file.c_str(), "update"); - string qa_dir_name = charge_str + "/" + bin_name + "/qa"; - qa_file.mkdir(qa_dir_name.c_str(), "dEdx fit qa", true); - - /** - * Extracting fit parameters - */ - map par_map_results; - auto fit_params = sim_pdf->getParameters(ds); - for (auto par : *fit_params) { - string par_name {par->GetName()}; - - double value = dynamic_cast(par)->getVal(); - double value_err = dynamic_cast(par)->getError(); - par_map_results.emplace(par_name, value); - par_map_results.emplace(par_name + "_err", value_err); - } - - auto c = new TCanvas; - c->Print("fit_qa.pdf(", "pdf"); - - qa_file.cd(qa_dir_name.c_str()); - fit_result->Write("fit_result"); - - gDirectory->WriteObject(&par_map_results, "fit_parameters"); - - - for (const auto &st : v_ncl_bin) { - std::cout << st->GetName() << std::endl; - auto frame_dedx = v_dedx.frame(); - ds.plotOn(frame_dedx, RooFit::Cut(std::string("ncl_bin==ncl_bin::").append(st->GetName()).c_str())); - sim_pdf->plotOn(frame_dedx, RooFit::Slice(v_ncl_bin, st->GetName()), RooFit::ProjWData(v_ncl_bin, ds)); - - qa_file.cd(qa_dir_name.c_str()); - frame_dedx->Write(Form("fit_%s", st->GetName())); - - c->Clear(); - c->Divide(2, 1); - c->cd(1); - frame_dedx->Draw(); - c->cd(2); - gPad->SetLogy(); - frame_dedx->Draw(); - c->Print("fit_qa.pdf", "pdf"); - } - - auto frame_dedx = v_dedx.frame(); - ds.plotOn(frame_dedx); - sim_pdf->plotOn(frame_dedx, RooFit::ProjWData(v_ncl_bin, ds)); - c->Clear(); - frame_dedx->Draw(); - c->Print("fit_qa.pdf", "pdf"); - - - c->Clear(); - c->Print("fit_qa.pdf)", "pdf"); - } + c->Clear(); + c->Print("fit_qa.pdf)", "pdf"); + } - cout << "Removing " << tmp_name << "..." << endl; - std::remove(tmp_name.c_str()); - return 0; + cout << "Removing " << tmp_name << "..." << endl; + std::remove(tmp_name.c_str()); + return 0; } diff --git a/pid_new/examples/shine/FitResultsExplorer.cpp b/pid_new/examples/shine/FitResultsExplorer.cpp index 94ec477..b88dad3 100644 --- a/pid_new/examples/shine/FitResultsExplorer.cpp +++ b/pid_new/examples/shine/FitResultsExplorer.cpp @@ -3,161 +3,148 @@ // #include +#include +#include #include +#include #include #include -#include -#include -#include #include using namespace std; - struct BinKinematics_t { - bool is_initialized {false}; + bool is_initialized{false}; - double p_lo{0.}; - double p_hi{0.}; + double p_lo{0.}; + double p_hi{0.}; - double pt_lo{0.}; - double pt_hi{0.}; + double pt_lo{0.}; + double pt_hi{0.}; - double phi_lo{0.}; - double phi_hi{0.}; + double phi_lo{0.}; + double phi_hi{0.}; }; struct Bin_t { - string name{""}; - TDirectory *bin_dir{nullptr}; + string name{""}; + TDirectory* bin_dir{nullptr}; - TTree *track_tree{nullptr}; - TDirectory *qa_dir{nullptr}; + TTree* track_tree{nullptr}; + TDirectory* qa_dir{nullptr}; - std::map fit_parameters; + std::map fit_parameters; - BinKinematics_t bin_kin_info; + BinKinematics_t bin_kin_info; }; -void readBinKinInfo(Bin_t &bin) { - assert(bin.track_tree); - - ROOT::RDataFrame rdf(*bin.track_tree); - - auto p_lo = rdf.Min("p"); - auto p_hi = rdf.Max("p"); - auto pt_lo = rdf.Min("pt"); - auto pt_hi = rdf.Max("pt"); - auto phi_lo = rdf.Min("phi"); - auto phi_hi = rdf.Max("phi"); - - bin.bin_kin_info.p_lo = *p_lo; - bin.bin_kin_info.p_hi = *p_hi; - bin.bin_kin_info.pt_lo = *pt_lo; - bin.bin_kin_info.pt_hi = *pt_hi; - bin.bin_kin_info.phi_lo = *phi_lo; - bin.bin_kin_info.phi_hi = *phi_hi; - bin.bin_kin_info.is_initialized = true; +void readBinKinInfo(Bin_t& bin) { + assert(bin.track_tree); + + ROOT::RDataFrame rdf(*bin.track_tree); + + auto p_lo = rdf.Min("p"); + auto p_hi = rdf.Max("p"); + auto pt_lo = rdf.Min("pt"); + auto pt_hi = rdf.Max("pt"); + auto phi_lo = rdf.Min("phi"); + auto phi_hi = rdf.Max("phi"); + + bin.bin_kin_info.p_lo = *p_lo; + bin.bin_kin_info.p_hi = *p_hi; + bin.bin_kin_info.pt_lo = *pt_lo; + bin.bin_kin_info.pt_hi = *pt_hi; + bin.bin_kin_info.phi_lo = *phi_lo; + bin.bin_kin_info.phi_hi = *phi_hi; + bin.bin_kin_info.is_initialized = true; } +vector lookupBins(TDirectory* dir) { + auto keys_list = dir->GetListOfKeys(); + vector lookup_result; -vector lookupBins(TDirectory *dir) { - - auto keys_list = dir->GetListOfKeys(); - vector lookup_result; - - for (auto key_tobj : *keys_list) { - auto key = static_cast(key_tobj); - string key_name {key_tobj->GetName()}; - - if ( - key_name.substr(0, strlen("bin_")) == "bin_" && - key->ReadObject() != nullptr - ) { - std::cout << "Looking up " << key_name << std::endl; - Bin_t bin; - bin.name = key_name; - bin.bin_dir = key->ReadObject(); + for (auto key_tobj : *keys_list) { + auto key = static_cast(key_tobj); + string key_name{key_tobj->GetName()}; - bin.track_tree = bin.bin_dir->Get("track_tree"); - bin.qa_dir = bin.bin_dir->Get("qa"); + if ( + key_name.substr(0, strlen("bin_")) == "bin_" && key->ReadObject() != nullptr) { + std::cout << "Looking up " << key_name << std::endl; + Bin_t bin; + bin.name = key_name; + bin.bin_dir = key->ReadObject(); - if (bin.qa_dir) { - auto fit_parameters_ptr = bin.qa_dir->Get>("fit_parameters"); + bin.track_tree = bin.bin_dir->Get("track_tree"); + bin.qa_dir = bin.bin_dir->Get("qa"); - if (fit_parameters_ptr) { - bin.fit_parameters = *fit_parameters_ptr; - } + if (bin.qa_dir) { + auto fit_parameters_ptr = bin.qa_dir->Get>("fit_parameters"); - } - - readBinKinInfo(bin); - lookup_result.push_back(bin); + if (fit_parameters_ptr) { + bin.fit_parameters = *fit_parameters_ptr; } + } - + readBinKinInfo(bin); + lookup_result.push_back(bin); } + } - return lookup_result; + return lookup_result; } -vector getMapKeys(const map &map) { - vector result; - transform(map.begin(), map.end(), back_inserter(result), [] (const pair &e) {return e.first; }); - return result; +vector getMapKeys(const map& map) { + vector result; + transform(map.begin(), map.end(), back_inserter(result), [](const pair& e) { return e.first; }); + return result; } int main() { - string input_file_name{"/home/eugene/STORAGE/data/pid/vtx_tack_dumps/pbpb_13agev_16_026/tracks_binned.root"}; - string charge_str {"negative"}; - - TFile input_file(input_file_name.c_str(), "read"); - auto tgt_dir = input_file.Get(charge_str.c_str()); - - auto bins_list = lookupBins(tgt_dir); - + string input_file_name{"/home/eugene/STORAGE/data/pid/vtx_tack_dumps/pbpb_13agev_16_026/tracks_binned.root"}; + string charge_str{"negative"}; + TFile input_file(input_file_name.c_str(), "read"); + auto tgt_dir = input_file.Get(charge_str.c_str()); + auto bins_list = lookupBins(tgt_dir); + BinKinematics_t bin_kin_info; - BinKinematics_t bin_kin_info; + TFile output_file("fit_results_explorer.root", "recreate"); + /* init output tree */ + map> fit_param_trees; + for (auto& bin : bins_list) { + bin_kin_info = bin.bin_kin_info; - TFile output_file("fit_results_explorer.root", "recreate"); - /* init output tree */ - map> fit_param_trees; - for (auto &bin : bins_list) { - bin_kin_info = bin.bin_kin_info; + for (auto& p : bin.fit_parameters) { + auto par_name = p.first; - for (auto &p : bin.fit_parameters) { - auto par_name = p.first; + auto par_tree_it = fit_param_trees.find(par_name); + if (fit_param_trees.end() == par_tree_it) { + auto par_tree = new TTree(par_name.c_str(), ""); + par_tree->Branch("p_lo", &bin_kin_info.p_lo); + par_tree->Branch("p_hi", &bin_kin_info.p_hi); + par_tree->Branch("pt_lo", &bin_kin_info.pt_lo); + par_tree->Branch("pt_hi", &bin_kin_info.pt_hi); + par_tree->Branch("phi_lo", &bin_kin_info.phi_lo); + par_tree->Branch("phi_hi", &bin_kin_info.phi_hi); - auto par_tree_it = fit_param_trees.find(par_name); - if (fit_param_trees.end() == par_tree_it) { - auto par_tree = new TTree(par_name.c_str(), ""); - par_tree->Branch("p_lo", &bin_kin_info.p_lo); - par_tree->Branch("p_hi", &bin_kin_info.p_hi); - par_tree->Branch("pt_lo", &bin_kin_info.pt_lo); - par_tree->Branch("pt_hi", &bin_kin_info.pt_hi); - par_tree->Branch("phi_lo", &bin_kin_info.phi_lo); - par_tree->Branch("phi_hi", &bin_kin_info.phi_hi); + par_tree_it = fit_param_trees.emplace(par_name, make_pair(par_tree, 0.)).first; + par_tree->Branch("par_value", &(par_tree_it->second.second)); + } - par_tree_it = fit_param_trees.emplace(par_name, make_pair(par_tree, 0.)).first; - par_tree->Branch("par_value", &(par_tree_it->second.second)); - } - - auto par_value = p.second; - par_tree_it->second.second = par_value; - par_tree_it->second.first->Fill(); - } + auto par_value = p.second; + par_tree_it->second.second = par_value; + par_tree_it->second.first->Fill(); } + } + for (auto& t : fit_param_trees) { + t.second.first->Write(); + } - for (auto &t : fit_param_trees) { - t.second.first->Write(); - } - - return 0; + return 0; } diff --git a/pid_new/examples/shine/RunFitShine.cpp b/pid_new/examples/shine/RunFitShine.cpp index 5019cbb..c157259 100644 --- a/pid_new/examples/shine/RunFitShine.cpp +++ b/pid_new/examples/shine/RunFitShine.cpp @@ -5,397 +5,398 @@ #include #include -#include -#include #include -#include #include +#include +#include #include #include +#include #include #include #include -#include #include +#include -double poly(const std::vector &pv, double x) { - double result = 0.; - double powX = 1.; +double poly(const std::vector& pv, double x) { + double result = 0.; + double powX = 1.; - for (auto p : pv) { - result += p * powX; - powX *= x; - } + for (auto p : pv) { + result += p * powX; + powX *= x; + } - return result; + return result; } FitParameter::ConstraintFct_t polyF(std::initializer_list params) { - std::vector pv(params.begin(), params.end()); - return [pv](double x) { return poly(pv, x); }; + std::vector pv(params.begin(), params.end()); + return [pv](double x) { return poly(pv, x); }; } - FitParameter::ConstraintFct_t wrapToX(FitParameter::ConstraintFct_t f, int charge) { - return [=](double x) { return f(TMath::Exp(charge * x) / 20); }; + return [=](double x) { return f(TMath::Exp(charge * x) / 20); }; } +int main(int argc, char** argv) { + + /* Building RooFit model */ + + /* observable */ + RooRealVar v_dEdx("x", "dEdx", 0, 2.5, "MIP"); + RooDataHist dh_dEdx("dh_dEdx", "", v_dEdx); + + /* common for all particle species asymmetry parameter */ + RooRealVar v_d("d", "#delta", -0.5, 0.5, ""); + + /* negative pions model */ + RooRealVar v_pion_neg_bb("pion_neg_bb", "BB(#pi^{-})", 0, "MIP"); + RooRealVar v_pion_neg_sigma("pion_neg_sigma", "#sigma(#pi^{-})", 0, "MIP"); + RooRealVar v_pion_neg_n("pion_neg_n", "N(#pi^{-})", 0, ""); + AsymmetricGaussianPDF pdf_pion_neg("pion_neg", "", v_dEdx, v_pion_neg_bb, v_pion_neg_sigma, v_d); + RooExtendPdf epdf_pion_neg("ext_pion_neg", "", pdf_pion_neg, v_pion_neg_n); + + /* negative kaons model */ + RooRealVar v_kaon_neg_bb("kaon_neg_bb", "BB(K^{-})", 0, "MIP"); + RooFormulaVar v_kaon_neg_sigma("kaon_neg_sigma", "#sigma(K^{-})", "@0*TMath::Power(@1/@2, 0.65)", + RooArgSet(v_pion_neg_sigma, v_kaon_neg_bb, v_pion_neg_bb)); + RooRealVar v_kaon_neg_n("kaon_neg_n", "N(K^{-})", 0, ""); + AsymmetricGaussianPDF pdf_kaon_neg("kaon_neg", "", v_dEdx, v_kaon_neg_bb, v_kaon_neg_sigma, v_d); + RooExtendPdf epdf_kaon_pdf("ext_kaon_neg", "", pdf_kaon_neg, v_kaon_neg_n); + + /* electrons model */ + RooRealVar v_electron_bb("electron_bb", "BB(e^{-})", 0, "MIP"); + RooFormulaVar v_electron_sigma("electron_sigma", "#sigma(e^{-})", "@0*TMath::Power(@1/@2, 0.65)", + RooArgSet(v_pion_neg_sigma, v_electron_bb, v_pion_neg_bb)); + RooRealVar v_electron_n("electron_n", "N(e^{-})", 0, ""); + AsymmetricGaussianPDF pdf_electron("electron", "", v_dEdx, v_electron_bb, v_electron_sigma, v_d); + RooExtendPdf epdf_electron("ext_electron", "", pdf_electron, v_electron_n); + + /* positive pions */ + RooRealVar v_pion_pos_bb("pion_pos_bb", "BB(#pi^{+})", 0, "MIP"); + RooRealVar v_pion_pos_sigma("pion_pos_sigma", "#sigma(#pi^{+})", 0, "MIP"); + RooRealVar v_pion_pos_n("pion_pos_n", "N(#pi^{+})", 0, ""); + AsymmetricGaussianPDF pdf_pion_pos("pion_pos", "", v_dEdx, v_pion_pos_bb, v_pion_pos_sigma, v_d); + RooExtendPdf epdf_pion_pos("ext_pion_pos", "", pdf_pion_pos, v_pion_pos_n); + + /* protons model */ + RooRealVar v_proton_bb("proton_bb", "BB(p)", 0, "MIP"); + RooFormulaVar v_proton_sigma("proton_sigma", "#sigma(p)", "@0*TMath::Power(@1/@2, 0.65)", + RooArgSet(v_pion_pos_sigma, v_proton_bb, v_pion_pos_bb)); + RooRealVar v_proton_n("proton_n", "N(p)", 0, ""); + AsymmetricGaussianPDF pdf_proton("proton", "", v_dEdx, v_proton_bb, v_proton_sigma, v_d); + RooExtendPdf epdf_proton("ext_proton", "", pdf_proton, v_proton_n); + + /* positrons model */ + RooRealVar v_positron_bb("positron_bb", "BB(e^{-})", 0, "MIP"); + RooFormulaVar v_positron_sigma("positron_sigma", "#sigma(e^{-})", "@0*TMath::Power(@1/@2, 0.65)", + RooArgSet(v_pion_pos_sigma, v_positron_bb, v_pion_pos_bb)); + RooRealVar v_positron_n("positron_n", "N(e^{-})", 0, ""); + AsymmetricGaussianPDF pdf_positron("positron", "", v_dEdx, v_positron_bb, v_positron_sigma, v_d); + RooExtendPdf epdf_positron("ext_positron", "", pdf_positron, v_positron_n); + + /* kaons model */ + RooRealVar v_kaon_pos_bb("kaon_pos_bb", "BB(K^{+})", 0, "MIP"); + RooFormulaVar v_kaon_pos_sigma("kaon_pos_sigma", "#sigma(K^{+})", "@0*TMath::Power(@1/@2, 0.65)", + RooArgSet(v_pion_pos_sigma, v_kaon_pos_bb, v_pion_pos_bb)); + RooRealVar v_kaon_pos_n("kaon_pos_n", "N(K^{+})", 0, ""); + AsymmetricGaussianPDF pdf_kaon_pos("kaon_pos", "", v_dEdx, v_kaon_pos_bb, v_kaon_pos_sigma, v_d); + RooExtendPdf epdf_kaon_pos_pdf("ext_kaon_pos", "", pdf_kaon_pos, v_kaon_pos_n); + + std::list fit_models; + + /* negative */ + fit_models.emplace_back("pion_neg", epdf_pion_neg, dh_dEdx, "pion_neg_"); + fit_models.back().setRange(-5.5, -0.7); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-211), -1)); + fit_models.back().parameter("sigma").fix(polyF({1.17732e-01, 8.11273e-03}), -5.5, -3.0); + fit_models.back().parameter("sigma").fix(polyF({1.28975e-01, 1.14362e-02}), -3.0, -2.6); + fit_models.back().parameter("sigma").fix(polyF({4.81931e-01, 7.07309e-01, 5.30289e-01, 1.81315e-01, 2.31737e-02}), -2.6, -0.4); + fit_models.back().parameter("n").fix(polyF({5.98777e+05, 2.17637e+05, 2.35418e+04, 1.51676e+03, 1.53068e+02}), -5.5, -4.5); + fit_models.back().parameter("n").fix(polyF({-5.56917e+04, -1.64847e+05, -1.78915e+05, -8.46502e+04, -1.34296e+04}), -2.2, -0.4); + + fit_models.emplace_back("kaon_neg", epdf_kaon_pdf, dh_dEdx, "kaon_neg_"); + fit_models.back().setRange(-5.3, -2.5); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-321), -1)); + fit_models.back().parameter("n").range(0, 2000); + + fit_models.emplace_back("electron", epdf_electron, dh_dEdx, "electron_"); + fit_models.back().setRange(-5., -0.4); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-11), -1)); + + /* positive */ + fit_models.emplace_back("pion_pos", epdf_pion_pos, dh_dEdx, "pion_pos_"); + fit_models.back().setRange(0.4, 5.5); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(211), 1)); + fit_models.back().parameter("sigma").fix(polyF({4.81931e-01, -7.07309e-01, 5.30289e-01, -1.81315e-01, 2.31737e-02}), 0.4, 2.6); + fit_models.back().parameter("sigma").fix(polyF({1.28975e-01, -1.14362e-02}), 2.6, 3.0); + fit_models.back().parameter("sigma").fix(polyF({1.17732e-01, -8.11273e-03}), 3.0, 5.5); + + fit_models.emplace_back("proton", epdf_proton, dh_dEdx, "proton_"); + fit_models.back().setRange(2.2, 5.5); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(2212), 1)); + + fit_models.emplace_back("positron", epdf_positron, dh_dEdx, "positron_"); + fit_models.back().setRange(0.4, 5.0); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(11), 1)); + + fit_models.emplace_back("kaon_pos", epdf_kaon_pos_pdf, dh_dEdx, "kaon_pos_"); + fit_models.back().setRange(2.5, 5.5); + fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(321), 1)); + fit_models.back().parameter("n").range(0, 2000); + + FitParameterRegistry::instance().par_find("d")->fix(0.05); + + /* IO */ + auto inputFile = TFile::Open(argv[1]); + TFile outputFile("output.root", "RECREATE"); + + TH2D* inputHistogram = nullptr; + inputFile->GetObject("reco_info/hTrackdEdx_allTPC", inputHistogram); + inputHistogram->SetDirectory(nullptr); + inputHistogram->Print(); + + inputHistogram->RebinX(1); + inputHistogram->RebinY(1); + + auto c = new TCanvas; + c->Print("output.pdf(", "pdf"); + + inputHistogram->Draw("colz"); + gPad->SetLogz(); + inputHistogram->GetYaxis()->SetRangeUser(0, 4.); + + { + TF1 bbPionNeg( + "bbPionNeg", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(-211), -1).operator()(x[0]); + }, + -6., -0.5, 0); + bbPionNeg.SetLineColor(kRed); + bbPionNeg.Draw("same"); + + TF1 bbPionPos( + "bbPionPos", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(211), 1).operator()(x[0]); + }, + 0.5, 6, 0); + bbPionPos.SetLineColor(kMagenta + 1); + bbPionPos.Draw("same"); + + TF1 bbProton( + "bbProton", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(2212), 1).operator()(x[0]); + }, + 2.5, 6, 0); + bbProton.SetLineColor(kBlack); + bbProton.Draw("same"); + + TF1 bbPositron( + "bbPositron", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(11), 1).operator()(x[0]); + }, + 0., 5., 0); + bbPositron.SetLineColor(kBlue + 1); + bbPositron.Draw("same"); + + TF1 bbElectron( + "bbElectron", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(-11), -1).operator()(x[0]); + }, + -5., 0., 0); + bbElectron.SetLineColor(kBlue + 1); + bbElectron.Draw("same"); + + TF1 bbDeuteron( + "bbDeuteron", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(1000010020), 1).operator()(x[0]); + }, + 3., 6., 0); + bbDeuteron.SetLineColor(kOrange + 1); + bbDeuteron.Draw("same"); + + TF1 bbKaonPos( + "bbKaonPos", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(321), 1).operator()(x[0]); + }, + 2., 5., 0); + bbKaonPos.SetLineColor(kCyan + 1); + bbKaonPos.Draw("same"); + + TF1 bbKaonNeg( + "bbKaonNeg", [](const Double_t* x, const Double_t* p) { + return wrapToX(BetheBlochHelper::makeBBForPdg(-321), -1).operator()(x[0]); + }, + -5., -2., 0); + bbKaonNeg.SetLineColor(kCyan + 1); + bbKaonNeg.Draw("same"); -int main(int argc, char **argv) { - - /* Building RooFit model */ - - /* observable */ - RooRealVar v_dEdx("x", "dEdx", 0, 2.5, "MIP"); - RooDataHist dh_dEdx("dh_dEdx", "", v_dEdx); - - /* common for all particle species asymmetry parameter */ - RooRealVar v_d("d", "#delta", -0.5, 0.5, ""); - - /* negative pions model */ - RooRealVar v_pion_neg_bb("pion_neg_bb", "BB(#pi^{-})", 0, "MIP"); - RooRealVar v_pion_neg_sigma("pion_neg_sigma", "#sigma(#pi^{-})", 0, "MIP"); - RooRealVar v_pion_neg_n("pion_neg_n", "N(#pi^{-})", 0, ""); - AsymmetricGaussianPDF pdf_pion_neg("pion_neg", "", v_dEdx, v_pion_neg_bb, v_pion_neg_sigma, v_d); - RooExtendPdf epdf_pion_neg("ext_pion_neg", "", pdf_pion_neg, v_pion_neg_n); - - /* negative kaons model */ - RooRealVar v_kaon_neg_bb("kaon_neg_bb", "BB(K^{-})", 0, "MIP"); - RooFormulaVar v_kaon_neg_sigma("kaon_neg_sigma", "#sigma(K^{-})", "@0*TMath::Power(@1/@2, 0.65)", - RooArgSet(v_pion_neg_sigma, v_kaon_neg_bb, v_pion_neg_bb)); - RooRealVar v_kaon_neg_n("kaon_neg_n", "N(K^{-})", 0, ""); - AsymmetricGaussianPDF pdf_kaon_neg("kaon_neg", "", v_dEdx, v_kaon_neg_bb, v_kaon_neg_sigma, v_d); - RooExtendPdf epdf_kaon_pdf("ext_kaon_neg", "", pdf_kaon_neg, v_kaon_neg_n); - - /* electrons model */ - RooRealVar v_electron_bb("electron_bb", "BB(e^{-})", 0, "MIP"); - RooFormulaVar v_electron_sigma("electron_sigma", "#sigma(e^{-})", "@0*TMath::Power(@1/@2, 0.65)", - RooArgSet(v_pion_neg_sigma, v_electron_bb, v_pion_neg_bb)); - RooRealVar v_electron_n("electron_n", "N(e^{-})", 0, ""); - AsymmetricGaussianPDF pdf_electron("electron", "", v_dEdx, v_electron_bb, v_electron_sigma, v_d); - RooExtendPdf epdf_electron("ext_electron", "", pdf_electron, v_electron_n); - - /* positive pions */ - RooRealVar v_pion_pos_bb("pion_pos_bb", "BB(#pi^{+})", 0, "MIP"); - RooRealVar v_pion_pos_sigma("pion_pos_sigma", "#sigma(#pi^{+})", 0, "MIP"); - RooRealVar v_pion_pos_n("pion_pos_n", "N(#pi^{+})", 0, ""); - AsymmetricGaussianPDF pdf_pion_pos("pion_pos", "", v_dEdx, v_pion_pos_bb, v_pion_pos_sigma, v_d); - RooExtendPdf epdf_pion_pos("ext_pion_pos", "", pdf_pion_pos, v_pion_pos_n); - - /* protons model */ - RooRealVar v_proton_bb("proton_bb", "BB(p)", 0, "MIP"); - RooFormulaVar v_proton_sigma("proton_sigma", "#sigma(p)", "@0*TMath::Power(@1/@2, 0.65)", - RooArgSet(v_pion_pos_sigma, v_proton_bb, v_pion_pos_bb)); - RooRealVar v_proton_n("proton_n", "N(p)", 0, ""); - AsymmetricGaussianPDF pdf_proton("proton", "", v_dEdx, v_proton_bb, v_proton_sigma, v_d); - RooExtendPdf epdf_proton("ext_proton", "", pdf_proton, v_proton_n); - - /* positrons model */ - RooRealVar v_positron_bb("positron_bb", "BB(e^{-})", 0, "MIP"); - RooFormulaVar v_positron_sigma("positron_sigma", "#sigma(e^{-})", "@0*TMath::Power(@1/@2, 0.65)", - RooArgSet(v_pion_pos_sigma, v_positron_bb, v_pion_pos_bb)); - RooRealVar v_positron_n("positron_n", "N(e^{-})", 0, ""); - AsymmetricGaussianPDF pdf_positron("positron", "", v_dEdx, v_positron_bb, v_positron_sigma, v_d); - RooExtendPdf epdf_positron("ext_positron", "", pdf_positron, v_positron_n); - - /* kaons model */ - RooRealVar v_kaon_pos_bb("kaon_pos_bb", "BB(K^{+})", 0, "MIP"); - RooFormulaVar v_kaon_pos_sigma("kaon_pos_sigma", "#sigma(K^{+})", "@0*TMath::Power(@1/@2, 0.65)", - RooArgSet(v_pion_pos_sigma, v_kaon_pos_bb, v_pion_pos_bb)); - RooRealVar v_kaon_pos_n("kaon_pos_n", "N(K^{+})", 0, ""); - AsymmetricGaussianPDF pdf_kaon_pos("kaon_pos", "", v_dEdx, v_kaon_pos_bb, v_kaon_pos_sigma, v_d); - RooExtendPdf epdf_kaon_pos_pdf("ext_kaon_pos", "", pdf_kaon_pos, v_kaon_pos_n); - - std::list fit_models; - - /* negative */ - fit_models.emplace_back("pion_neg", epdf_pion_neg, dh_dEdx, "pion_neg_"); - fit_models.back().setRange(-5.5, -0.7); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-211), -1)); - fit_models.back().parameter("sigma").fix(polyF({1.17732e-01, 8.11273e-03}), -5.5, -3.0); - fit_models.back().parameter("sigma").fix(polyF({1.28975e-01, 1.14362e-02}), -3.0, -2.6); - fit_models.back().parameter("sigma").fix(polyF({4.81931e-01,7.07309e-01,5.30289e-01,1.81315e-01,2.31737e-02}), -2.6, -0.4); - fit_models.back().parameter("n").fix(polyF({5.98777e+05, 2.17637e+05, 2.35418e+04, 1.51676e+03, 1.53068e+02}), -5.5, -4.5); - fit_models.back().parameter("n").fix(polyF({-5.56917e+04, -1.64847e+05, -1.78915e+05, -8.46502e+04, -1.34296e+04}), -2.2, -0.4); - - fit_models.emplace_back("kaon_neg", epdf_kaon_pdf, dh_dEdx, "kaon_neg_"); - fit_models.back().setRange(-5.3, -2.5); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-321), -1)); - fit_models.back().parameter("n").range(0,2000); - - fit_models.emplace_back("electron", epdf_electron, dh_dEdx, "electron_"); - fit_models.back().setRange(-5., -0.4); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(-11), -1)); - - /* positive */ - fit_models.emplace_back("pion_pos", epdf_pion_pos, dh_dEdx, "pion_pos_"); - fit_models.back().setRange(0.4, 5.5); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(211), 1)); - fit_models.back().parameter("sigma").fix(polyF({4.81931e-01,-7.07309e-01,5.30289e-01,-1.81315e-01,2.31737e-02}), 0.4, 2.6); - fit_models.back().parameter("sigma").fix(polyF({1.28975e-01, -1.14362e-02}), 2.6, 3.0); - fit_models.back().parameter("sigma").fix(polyF({1.17732e-01, -8.11273e-03}), 3.0, 5.5); - - fit_models.emplace_back("proton", epdf_proton, dh_dEdx, "proton_"); - fit_models.back().setRange(2.2, 5.5); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(2212), 1)); - - fit_models.emplace_back("positron", epdf_positron, dh_dEdx, "positron_"); - fit_models.back().setRange(0.4, 5.0); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(11), 1)); - - fit_models.emplace_back("kaon_pos", epdf_kaon_pos_pdf, dh_dEdx, "kaon_pos_"); - fit_models.back().setRange(2.5, 5.5); - fit_models.back().parameter("bb").fix(wrapToX(BetheBlochHelper::makeBBForPdg(321), 1)); - fit_models.back().parameter("n").range(0,2000); - - FitParameterRegistry::instance().par_find("d")->fix(0.05); - - - /* IO */ - auto inputFile = TFile::Open(argv[1]); - TFile outputFile("output.root", "RECREATE"); - - TH2D *inputHistogram = nullptr; - inputFile->GetObject("reco_info/hTrackdEdx_allTPC", inputHistogram); - inputHistogram->SetDirectory(nullptr); - inputHistogram->Print(); - - - inputHistogram->RebinX(1); - inputHistogram->RebinY(1); - - auto c = new TCanvas; c->Print("output.pdf(", "pdf"); - - inputHistogram->Draw("colz"); - gPad->SetLogz(); - inputHistogram->GetYaxis()->SetRangeUser(0, 4.); - - { - TF1 bbPionNeg("bbPionNeg", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(-211), -1).operator()(x[0]); - }, -6., -0.5, 0); - bbPionNeg.SetLineColor(kRed); - bbPionNeg.Draw("same"); - - TF1 bbPionPos("bbPionPos", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(211), 1).operator()(x[0]); - }, 0.5, 6, 0); - bbPionPos.SetLineColor(kMagenta + 1); - bbPionPos.Draw("same"); - - TF1 bbProton("bbProton", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(2212), 1).operator()(x[0]); - }, 2.5, 6, 0); - bbProton.SetLineColor(kBlack); - bbProton.Draw("same"); - - TF1 bbPositron("bbPositron", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(11), 1).operator()(x[0]); - }, 0., 5., 0); - bbPositron.SetLineColor(kBlue + 1); - bbPositron.Draw("same"); - - - TF1 bbElectron("bbElectron", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(-11), -1).operator()(x[0]); - }, -5., 0., 0); - bbElectron.SetLineColor(kBlue + 1); - bbElectron.Draw("same"); - - TF1 bbDeuteron("bbDeuteron", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(1000010020), 1).operator()(x[0]); - }, 3., 6., 0); - bbDeuteron.SetLineColor(kOrange + 1); - bbDeuteron.Draw("same"); - - TF1 bbKaonPos("bbKaonPos", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(321), 1).operator()(x[0]); - }, 2., 5., 0); - bbKaonPos.SetLineColor(kCyan + 1); - bbKaonPos.Draw("same"); - - TF1 bbKaonNeg("bbKaonNeg", [](const Double_t *x, const Double_t *p) { - return wrapToX(BetheBlochHelper::makeBBForPdg(-321), -1).operator()(x[0]); - }, -5., -2., 0); - bbKaonNeg.SetLineColor(kCyan + 1); - bbKaonNeg.Draw("same"); - - c->Print("output.pdf(", "pdf"); - c->Clear(); + c->Clear(); + } + + for (int i = 1; i < inputHistogram->GetXaxis()->GetNbins(); ++i) { + double x = inputHistogram->GetXaxis()->GetBinCenter(i); + double xlo = inputHistogram->GetXaxis()->GetBinLowEdge(i); + double xhi = inputHistogram->GetXaxis()->GetBinUpEdge(i); + + std::cout << "Bin: " << i << "; x = " << x << std::endl; + + /* building composite model */ + RooArgList composite_model_components; + std::for_each(fit_models.begin(), fit_models.end(), [x, &composite_model_components](ParticleFitModel& m) { + if (m.isDefinedAt(x)) { + composite_model_components.add(m.pdf()); + } + }); + + if (composite_model_components.empty()) { + std::cout << "Nothing to fit in this momentum bin" << std::endl; + continue; } + RooAddPdf pdf_composite(Form("composite_%d", i), "", composite_model_components); + ParticleFitModel m_composite("composite", pdf_composite, dh_dEdx); + + /* preparing input */ + auto py = inputHistogram->ProjectionY("tmp", i, i); + py->SetDirectory(nullptr); + + dh_dEdx.reset(); + RooDataHist slice_data("slice_dh", "", v_dEdx, py); + dh_dEdx.add(slice_data); + + auto qa_dir = outputFile.mkdir(Form("bin_%d", i)); + + auto slice_integral = dh_dEdx.sumEntries(); + + /* set range for the component integrals */ + auto integral_params = m_composite.parameters_matching(".*_n$"); + std::for_each(integral_params.begin(), integral_params.end(), [=](FitParameter* p) { + auto& constr_lo = p->findConstraint(xlo); + auto& constr_hi = p->findConstraint(xhi); + + if ( + constr_lo.type_ == FitParameter::EConstraintType::kNone && constr_hi.type_ == FitParameter::EConstraintType::kNone) { + p->range(0., slice_integral, xlo, xhi); + } + }); + /* constraints */ + auto model_params = m_composite.parameters(); + std::for_each(model_params.begin(), model_params.end(), [=](FitParameter* p) { p->applyConstraint(x); }); + + m_composite.pdf().fitTo(dh_dEdx, RooFit::Extended()); + std::for_each(model_params.begin(), model_params.end(), [=](FitParameter* p) { p->pickFitResultAt(x); }); + + auto frame_data_fit = v_dEdx.frame(); + dh_dEdx.plotOn(frame_data_fit, RooFit::Name("data")); + m_composite.pdf().plotOn(frame_data_fit, RooFit::Name("fit")); + + for (auto& m : fit_models) { + if (m.isDefinedAt(x)) { + m.pdf().plotOn(frame_data_fit, + RooFit::Normalization(1., RooAbsReal::RelativeExpected), + RooFit::LineColor(m.getColor())); + } + } - for (int i = 1; i < inputHistogram->GetXaxis()->GetNbins(); ++i) { - double x = inputHistogram->GetXaxis()->GetBinCenter(i); - double xlo = inputHistogram->GetXaxis()->GetBinLowEdge(i); - double xhi = inputHistogram->GetXaxis()->GetBinUpEdge(i); - - std::cout << "Bin: " << i << "; x = " << x << std::endl; - - - - /* building composite model */ - RooArgList composite_model_components; - std::for_each(fit_models.begin(), fit_models.end(), [x, &composite_model_components](ParticleFitModel &m) { - if (m.isDefinedAt(x)) { - composite_model_components.add(m.pdf()); - } - }); - - if (composite_model_components.empty()) { - std::cout << "Nothing to fit in this momentum bin" << std::endl; - continue; - } - - RooAddPdf pdf_composite(Form("composite_%d", i), "", composite_model_components); - ParticleFitModel m_composite("composite", pdf_composite, dh_dEdx); - - /* preparing input */ - auto py = inputHistogram->ProjectionY("tmp", i, i); - py->SetDirectory(nullptr); - - dh_dEdx.reset(); - RooDataHist slice_data("slice_dh", "", v_dEdx, py); - dh_dEdx.add(slice_data); - - auto qa_dir = outputFile.mkdir(Form("bin_%d", i)); - - auto slice_integral = dh_dEdx.sumEntries(); - - /* set range for the component integrals */ - auto integral_params = m_composite.parameters_matching(".*_n$"); - std::for_each(integral_params.begin(), integral_params.end(), [=](FitParameter *p) { - auto &constr_lo = p->findConstraint(xlo); - auto &constr_hi = p->findConstraint(xhi); - - if ( - constr_lo.type_ == FitParameter::EConstraintType::kNone && - constr_hi.type_ == FitParameter::EConstraintType::kNone - ) { - p->range(0., slice_integral, xlo, xhi); - } - }); - /* constraints */ - auto model_params = m_composite.parameters(); - std::for_each(model_params.begin(), model_params.end(), [=](FitParameter *p) { p->applyConstraint(x); }); - - m_composite.pdf().fitTo(dh_dEdx, RooFit::Extended()); - std::for_each(model_params.begin(), model_params.end(), [=](FitParameter *p) { p->pickFitResultAt(x); }); - - auto frame_data_fit = v_dEdx.frame(); - dh_dEdx.plotOn(frame_data_fit, RooFit::Name("data")); - m_composite.pdf().plotOn(frame_data_fit, RooFit::Name("fit")); - - for (auto &m: fit_models) { - if (m.isDefinedAt(x)) { - m.pdf().plotOn(frame_data_fit, - RooFit::Normalization(1., RooAbsReal::RelativeExpected), - RooFit::LineColor(m.getColor()) - ); - } - } - - - qa_dir->WriteObject(frame_data_fit, "data_fit"); - - TPaveText pave_text(0.1, 0.9, 0.4, 0.98); - pave_text.AddText(Form("x = %f", x)); - - c->cd(); - c->Divide(2, 2); - - c->cd(1); - frame_data_fit->Draw(); - pave_text.DrawClone(); - - c->cd(2); - frame_data_fit->Draw(); - gPad->SetLogy(1); - - c->cd(3); - auto hist_data = dh_dEdx.createHistogram("tmp_data", v_dEdx, RooFit::IntrinsicBinning()); - auto hist_fit = m_composite.pdf().createHistogram("tmp_fit", v_dEdx, RooFit::IntrinsicBinning(), - RooFit::Extended(kTRUE), RooFit::Scaling(kTRUE)); - hist_fit->Scale(v_dEdx.getBinWidth(1)); - - TRatioPlot rp(hist_fit, hist_data, "divsym"); - rp.Draw("nohide,grid"); - rp.GetLowerRefGraph()->SetMinimum(0.5); - rp.GetLowerRefGraph()->SetMaximum(1.5); - - - c->cd(4); - - gStyle->SetOptStat(0); - gStyle->SetOptTitle(0); + qa_dir->WriteObject(frame_data_fit, "data_fit"); + + TPaveText pave_text(0.1, 0.9, 0.4, 0.98); + pave_text.AddText(Form("x = %f", x)); + + c->cd(); + c->Divide(2, 2); + + c->cd(1); + frame_data_fit->Draw(); + pave_text.DrawClone(); + + c->cd(2); + frame_data_fit->Draw(); + gPad->SetLogy(1); + + c->cd(3); + auto hist_data = dh_dEdx.createHistogram("tmp_data", v_dEdx, RooFit::IntrinsicBinning()); + auto hist_fit = m_composite.pdf().createHistogram("tmp_fit", v_dEdx, RooFit::IntrinsicBinning(), + RooFit::Extended(kTRUE), RooFit::Scaling(kTRUE)); + hist_fit->Scale(v_dEdx.getBinWidth(1)); + + TRatioPlot rp(hist_fit, hist_data, "divsym"); + rp.Draw("nohide,grid"); + rp.GetLowerRefGraph()->SetMinimum(0.5); + rp.GetLowerRefGraph()->SetMaximum(1.5); + + c->cd(4); + + gStyle->SetOptStat(0); + gStyle->SetOptTitle(0); + + TLegend species_legend(0.3, 0.85, 0.7, 0.95); + species_legend.SetNColumns(2); + + bool same_flag = false; + bool legend_flag = false; + /* calculate purity */ + for (auto& m : fit_models) { + if (m.isDefinedAt(x)) { + auto model_color = m.getColor(); + + auto particle_pdf_hist = (TH1F*) m.pdf().createHistogram("tmp_component", v_dEdx, + RooFit::IntrinsicBinning()); + auto composite_pdf_hist = (TH1F*) m_composite.pdf().createHistogram("tmp_composite", v_dEdx, + RooFit::IntrinsicBinning(), + RooFit::Extended()); + auto data_hist = (TH1F*) dh_dEdx.createHistogram("tmp_data", v_dEdx, RooFit::IntrinsicBinning()); + qa_dir->WriteObject(particle_pdf_hist, Form("pdf_%s", m.getName().c_str())); + + auto purity_hist = (*particle_pdf_hist) / (*composite_pdf_hist); + purity_hist.SetName(Form("purity_%s", m.getName().c_str())); + purity_hist.SetLineColor(model_color); + purity_hist.SetMarkerColor(model_color); + purity_hist.SetLineWidth(2); + qa_dir->WriteObject(&purity_hist, purity_hist.GetName()); + + particle_pdf_hist->Scale(v_dEdx.getBinWidth(1)); + auto purity_hist_data = (*particle_pdf_hist) / (*data_hist); + purity_hist_data.SetName(Form("purity_to_data_%s", m.getName().c_str())); + purity_hist_data.SetLineColor(model_color); + purity_hist_data.SetMarkerColor(model_color); + qa_dir->WriteObject(&purity_hist_data, purity_hist_data.GetName()); + + purity_hist_data.GetYaxis()->SetRangeUser(0., 2.); + purity_hist_data.DrawClone(same_flag ? "same,LE1" : "LE1"); + same_flag = true; + purity_hist.DrawClone("same"); + + species_legend.AddEntry(Form("purity_%s", m.getName().c_str()), m.getName().c_str(), "l"); + + delete particle_pdf_hist; + delete composite_pdf_hist; + delete data_hist; + } + } - TLegend species_legend(0.3, 0.85, 0.7, 0.95); - species_legend.SetNColumns(2); - - bool same_flag = false; - bool legend_flag = false; - /* calculate purity */ - for (auto &m : fit_models) { - if (m.isDefinedAt(x)) { - auto model_color = m.getColor(); - - auto particle_pdf_hist = (TH1F *) m.pdf().createHistogram("tmp_component", v_dEdx, - RooFit::IntrinsicBinning()); - auto composite_pdf_hist = (TH1F *) m_composite.pdf().createHistogram("tmp_composite", v_dEdx, - RooFit::IntrinsicBinning(), - RooFit::Extended()); - auto data_hist = (TH1F *) dh_dEdx.createHistogram("tmp_data", v_dEdx, RooFit::IntrinsicBinning()); - qa_dir->WriteObject(particle_pdf_hist, Form("pdf_%s", m.getName().c_str())); - - auto purity_hist = (*particle_pdf_hist) / (*composite_pdf_hist); - purity_hist.SetName(Form("purity_%s", m.getName().c_str())); - purity_hist.SetLineColor(model_color); - purity_hist.SetMarkerColor(model_color); - purity_hist.SetLineWidth(2); - qa_dir->WriteObject(&purity_hist, purity_hist.GetName()); - - particle_pdf_hist->Scale(v_dEdx.getBinWidth(1)); - auto purity_hist_data = (*particle_pdf_hist) / (*data_hist); - purity_hist_data.SetName(Form("purity_to_data_%s", m.getName().c_str())); - purity_hist_data.SetLineColor(model_color); - purity_hist_data.SetMarkerColor(model_color); - qa_dir->WriteObject(&purity_hist_data, purity_hist_data.GetName()); - - purity_hist_data.GetYaxis()->SetRangeUser(0., 2.); - purity_hist_data.DrawClone(same_flag ? "same,LE1" : "LE1"); - same_flag = true; - purity_hist.DrawClone("same"); - - species_legend.AddEntry(Form("purity_%s", m.getName().c_str()), m.getName().c_str(), "l"); - - - delete particle_pdf_hist; - delete composite_pdf_hist; - delete data_hist; - } - } - - species_legend.Draw(); + species_legend.Draw(); - c->Print("output.pdf", "pdf"); - c->Clear(); - - delete hist_fit; - delete hist_data; - delete py; - } + c->Print("output.pdf", "pdf"); + c->Clear(); + delete hist_fit; + delete hist_data; + delete py; + } - std::for_each(FitParameterRegistry::instance().par_begin(), FitParameterRegistry::instance().par_end(), - [&outputFile, c](FitParameter &p) { - p.dumpResult(outputFile, c); - c->Print("output.pdf", "pdf"); - c->Clear(); - }); + std::for_each(FitParameterRegistry::instance().par_begin(), FitParameterRegistry::instance().par_end(), + [&outputFile, c](FitParameter& p) { + p.dumpResult(outputFile, c); + c->Print("output.pdf", "pdf"); + c->Clear(); + }); - c->Print("output.pdf)", "pdf"); - return 0; + c->Print("output.pdf)", "pdf"); + return 0; } diff --git a/src/CMakeLists.txt b/src_tof/CMakeLists.txt similarity index 90% rename from src/CMakeLists.txt rename to src_tof/CMakeLists.txt index 57b75cf..5544359 100644 --- a/src/CMakeLists.txt +++ b/src_tof/CMakeLists.txt @@ -2,11 +2,11 @@ set(SOURCES Fitter.cpp ParticleFit.cpp Parameters.cpp - Getter.cpp + GetterTof.cpp ) string(REPLACE ".cpp" ".h" HEADERS "${SOURCES}") -list(APPEND HEADERS Constants.h) +list(APPEND HEADERS ConstantsTof.h) set(DICT_FILE_NAME G__${PROJECT_NAME}) set(PCM_FILE_NAME lib${PROJECT_NAME}) @@ -27,7 +27,7 @@ endforeach () target_link_libraries(Pid ${ROOT_LIBRARIES}) target_include_directories(Pid - PRIVATE ${CMAKE_SOURCE_DIR}/src + PRIVATE ${CMAKE_SOURCE_DIR}/src_tof PUBLIC $ $ @@ -42,4 +42,4 @@ install( OPTIONAL ) -install(FILES ${HEADERS} DESTINATION include) \ No newline at end of file +install(FILES ${HEADERS} DESTINATION include) diff --git a/src/Constants.h b/src_tof/ConstantsTof.h similarity index 85% rename from src/Constants.h rename to src_tof/ConstantsTof.h index 311c615..55f3d23 100644 --- a/src/Constants.h +++ b/src_tof/ConstantsTof.h @@ -1,15 +1,15 @@ -/** @file Constants.h +/** @file ConstantsTof.h @author Viktor Klochkov (klochkov44@gmail.com) @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com) @brief Some constants and enumerators */ -#ifndef PidConstants_H -#define PidConstants_H 1 +#ifndef PidConstantsTof_H +#define PidConstantsTof_H 1 #include -namespace PidParticles { +namespace PidTofParticles { enum eParticle { kBgPos = 1, kBgNeg = -1, @@ -38,7 +38,7 @@ const std::unordered_map masses = {kHe3, 2.793}, {kHe4, 3.724}, }; -}// namespace PidParticles +}// namespace PidTofParticles namespace PidFunction { enum eNames { diff --git a/src/Container.cpp b/src_tof/ContainerTof.cpp similarity index 67% rename from src/Container.cpp rename to src_tof/ContainerTof.cpp index 37a6465..9714cb2 100644 --- a/src/Container.cpp +++ b/src_tof/ContainerTof.cpp @@ -1,4 +1,4 @@ -#include "Container.h" +#include "ContainerTof.h" // ClassImp(Centrality::Container) namespace Pid { diff --git a/src/Container.h b/src_tof/ContainerTof.h similarity index 83% rename from src/Container.h rename to src_tof/ContainerTof.h index 1ee5909..7344792 100644 --- a/src/Container.h +++ b/src_tof/ContainerTof.h @@ -1,8 +1,8 @@ -/** @file Container.h - @class Pid::Container +/** @file ContainerTof.h + @class PidTof::Container @author Viktor Klochkov (klochkov44@gmail.com) @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com) - @brief Class to store PID information + @brief Class to store Tof PID information */ #ifndef CENTRALITY_CONTAINER_H @@ -12,9 +12,9 @@ typedef unsigned int uint; -namespace Pid { +namespace PidTof { -class Container { +class ContainerTof { public: Container() = default; @@ -38,6 +38,6 @@ class Container { // ClassDef(Container, 1); }; -}// namespace Pid +}// namespace PidTof #endif//CENTRALITY_CONTAINER_H diff --git a/src/FitHelper.h b/src_tof/FitHelper.h similarity index 98% rename from src/FitHelper.h rename to src_tof/FitHelper.h index eb62e01..640fff0 100644 --- a/src/FitHelper.h +++ b/src_tof/FitHelper.h @@ -88,7 +88,7 @@ Pid::ParticleFit FitPionsPos(Pid::Fitter& tof, const std::shared_ptr& hist pions.xmax = 10.5; pions.ymin = -2.; pions.ymax = 2.; - pions.pid = PidParticles::kPionPos; + pions.pid = PidTofParticles::kPionPos; pions.name = "piplus"; return FitParticle(tof, pions, hist); diff --git a/src/Fitter.cpp b/src_tof/Fitter.cpp similarity index 100% rename from src/Fitter.cpp rename to src_tof/Fitter.cpp diff --git a/src/Fitter.h b/src_tof/Fitter.h similarity index 100% rename from src/Fitter.h rename to src_tof/Fitter.h diff --git a/src/Getter.cpp b/src_tof/GetterTof.cpp similarity index 86% rename from src/Getter.cpp rename to src_tof/GetterTof.cpp index be9018a..31782ff 100644 --- a/src/Getter.cpp +++ b/src_tof/GetterTof.cpp @@ -1,6 +1,6 @@ -#include "Getter.h" +#include "GetterTof.h" -ClassImp(Pid::Getter) +ClassImp(Pid::GetterTof) ClassImp(Pid::CutGGetter) namespace Pid { @@ -11,7 +11,7 @@ ClassImp(Pid::Getter) * @param m2 square mass for TOF (or y-axis value, in general) * @return map with probabilities for all particle species */ - std::map Getter::GetBayesianProbability(double p, double m2) { + std::map GetterTof::GetBayesianProbability(double p, double m2) { std::map prob{}; // if (p > maxx_ || p < minx_) //NOTE should be done in cuts maybe? diff --git a/src/Getter.h b/src_tof/GetterTof.h similarity index 93% rename from src/Getter.h rename to src_tof/GetterTof.h index 57a0465..891bb88 100644 --- a/src/Getter.h +++ b/src_tof/GetterTof.h @@ -1,12 +1,12 @@ -/** @file Getter.h - @class Pid::Getter +/** @file GetterTof.h + @class Pid::GetterTof @author Viktor Klochkov (klochkov44@gmail.com) @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com) - @brief Class to calculate PID probabilities + @brief Class to calculate Tof-PID probabilities */ -#ifndef PidGetter_H -#define PidGetter_H 1 +#ifndef PidGetterTof_H +#define PidGetterTof_H 1 #include "ParticleFit.h" #include "TObject.h" @@ -37,7 +37,7 @@ class BaseGetter { /** * @brief Bayesian Pid getter */ -class Getter : public TObject, public BaseGetter { +class GetterTof : public TObject, public BaseGetter { public: void AddParticle(const ParticleFit& particle, uint id) { species_[id] = particle; } void AddParticles(std::map&& species) { species_ = species; } @@ -85,7 +85,7 @@ class Getter : public TObject, public BaseGetter { double minx_{-100000.}; double maxx_{100000.}; - ClassDefOverride(Getter, 1); + ClassDefOverride(GetterTof, 1); }; /** @@ -168,4 +168,4 @@ class CutGGetter : public TObject, public BaseGetter { }; }// namespace Pid -#endif// PidGetter_H +#endif// PidGetterTof_H diff --git a/src/Parameters.cpp b/src_tof/Parameters.cpp similarity index 100% rename from src/Parameters.cpp rename to src_tof/Parameters.cpp diff --git a/src/Parameters.h b/src_tof/Parameters.h similarity index 100% rename from src/Parameters.h rename to src_tof/Parameters.h diff --git a/src/ParticleFit.cpp b/src_tof/ParticleFit.cpp similarity index 100% rename from src/ParticleFit.cpp rename to src_tof/ParticleFit.cpp diff --git a/src/ParticleFit.h b/src_tof/ParticleFit.h similarity index 98% rename from src/ParticleFit.h rename to src_tof/ParticleFit.h index a9df58a..5b1239b 100644 --- a/src/ParticleFit.h +++ b/src_tof/ParticleFit.h @@ -11,7 +11,7 @@ #include #include -#include "Constants.h" +#include "ConstantsTof.h" #include "TF1.h" namespace Pid { diff --git a/src/PidLinkDef.h b/src_tof/PidLinkDef.h similarity index 87% rename from src/PidLinkDef.h rename to src_tof/PidLinkDef.h index 3a6eb0a..a4a6c08 100644 --- a/src/PidLinkDef.h +++ b/src_tof/PidLinkDef.h @@ -7,7 +7,7 @@ #pragma link C++ class Pid::Fitter + ; #pragma link C++ class Pid::ParticleFit + ; #pragma link C++ class Pid::BaseGetter--; -#pragma link C++ class Pid::Getter + ; +#pragma link C++ class Pid::GetterTof + ; #pragma link C++ class Pid::CutGGetter + ; #endif diff --git a/src_trd/CMakeLists.txt b/src_trd/CMakeLists.txt new file mode 100644 index 0000000..5712a85 --- /dev/null +++ b/src_trd/CMakeLists.txt @@ -0,0 +1,35 @@ +set(SOURCES + ParticleProb.cpp + ContainerTrd.cpp + GetterTrd.cpp + ) + +string(REPLACE ".cpp" ".h" HEADERS "${SOURCES}") +list(APPEND HEADERS ConstantsTrd.h) +include_directories(${PROJECT_INCLUDE_DIRECTORIES} ${CMAKE_CURRENT_SOURCE_DIR}) +add_library(PidTrd SHARED ${SOURCES} G__PidTrd.cxx) + +message("${PROJECT_INCLUDE_DIRECTORIES}") + +ROOT_GENERATE_DICTIONARY(G__PidTrd ${HEADERS} LINKDEF PidTrdLinkDef.h) +target_link_libraries(PidTrd ${ROOT_LIBRARIES}) + +install( + FILES + ${HEADERS} + DESTINATION + include + COMPONENT + Devel +) + +set(PCM_FILE_NAME libPidTrd) + +install( + FILES + "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}_rdict.pcm" + "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}.rootmap" + DESTINATION + lib + OPTIONAL +) diff --git a/src_trd/ConstantsTrd.h b/src_trd/ConstantsTrd.h new file mode 100644 index 0000000..087e86d --- /dev/null +++ b/src_trd/ConstantsTrd.h @@ -0,0 +1,75 @@ +/** @file ConstantsTrd.h + @author Susanne Glaessel (glaessel@ikf.uni-frankfurt.de) + @brief Some constants and enumerators +*/ + +#ifndef ConstantsTrd_H +#define ConstantsTrd_H + +#include "TString.h" +#include +#include + +using std::make_pair; + +namespace PidTrdParticles { +enum eParticle { + kProton = 2212, + kPionPos = 211, + kKaonPos = 321, + kDeutron = 1000010020, + kTriton = 1000010030, + kHe3 = 1000020030, + kHe4 = 1000020040, + kElectronPos = -11, + kMuonPos = -13, + kBgPos = 1, + kAntiProton = -2212, + kPionNeg = -211, + kKaonNeg = -321, + kDeutronNeg = -1000010020, + kTritonNeg = -1000010030, + kHe3Neg = -1000020030, + kHe4Neg = -1000020040, + kElectronNeg = 11, + kMuonNeg = 13, + kBgNeg = -1 +}; +}; +constexpr int NbinsMax = 100; +constexpr int NumberOfTrdLayers = 4; +constexpr int NumberOfPidsTrd = 10; +constexpr int NumberOfTruncMode = 5; +constexpr int NumberOfProbMode = 2; + +constexpr float BwMom = 0.05; +constexpr int NbinsMom = 60; +constexpr float BwdEdx = 0.1; +constexpr int NbinsdEdx = 35; + +inline std::array, NumberOfPidsTrd> pid_codes_trd_{ + make_pair(2212, "p"), + make_pair(211, "pi"), + make_pair(321, "K"), + make_pair(1000010020, "d"), + make_pair(1000010030, "t"), + make_pair(1000020030, "he3"), + make_pair(1000020040, "he4"), + make_pair(-11, "e"), + make_pair(-13, "mu"), + make_pair(1, "bg")}; + +inline TString dirname_tracks_ = "dEdx_tracks"; +inline TString dirname_hits_ = "dEdx_hits"; +inline std::array dirname_nhits_ = { + "1hit", + "2hits", + "3hits", + "4hits"}; +inline std::array histtitle_mode_ = { + "truncation mode: 1 hit", + "truncation mode: 2 hits", + "truncation mode: 3 hits", + "truncation mode: 4 hits"}; + +#endif// Constants_H diff --git a/src_trd/ContainerTrd.cpp b/src_trd/ContainerTrd.cpp new file mode 100644 index 0000000..d559281 --- /dev/null +++ b/src_trd/ContainerTrd.cpp @@ -0,0 +1,139 @@ +#include "ContainerTrd.h" + +/** +* Scales energy loss to the length of its trace through the Trd +*/ + +void TrdContainer::ScaleEnergyLossLength() { + + if (dEdx_is_scaled_ == true) { + Warning("Exec", "Energy loss is already scaled."); + return; + } + + float pz = 0.0; + if (TMath::Abs(mom_) > TMath::Abs(pT_)) + pz = TMath::Sqrt(mom_ * mom_ - pT_ * pT_); + else + Warning("Exec", "Could not assign pz to the track, use pz=0. Energy loss will not be scaled."); + + if (pz > 0 && mom_ > 0) { + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + dEdx_hits_.at(ihit) *= pz / mom_; + + dEdx_is_scaled_ = true; + } +} + +/** +* Sorts hits of a track depending on its dEdx-value from lowest to highest +*/ + +std::array TrdContainer::GetdEdxHitsSorted() { + + std::vector dEdx; + dEdx.clear(); + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + if (dEdx_hits_.at(ihit) > 0) + dEdx.push_back(dEdx_hits_.at(ihit)); + + std::sort(dEdx.begin(), dEdx.end()); + + for (int ihit = 0; ihit < dEdx.size(); ihit++) + hits_sorted_.at(ihit) = dEdx.at(ihit); + + for (int ihit = dEdx.size(); ihit < NumberOfTrdLayers; ihit++) + hits_sorted_.at(ihit) = 0.0; + + return hits_sorted_; +} + +/** +* Gets number of hits for a track according to the truncation mode +* @param truncation mode +* @return number of selected hits +*/ + +int TrdContainer::GetNHitsSel(int trunc_mode) { + + if (trunc_mode == 0) + nhits_sel_ = nhits_trd_; + else + nhits_sel_ = trunc_mode; + if (nhits_sel_ > nhits_trd_) + nhits_sel_ = nhits_trd_; + + return nhits_sel_; +} + +/** +* Selects hits of track according to the truncation mode +* @param truncation mode +*/ + +void TrdContainer::SelectHitIndices(int trunc_mode) { + + if (nhits_sel_ == 0) nhits_sel_ = GetNHitsSel(trunc_mode); + + std::array dEdx_hits_tmp = dEdx_hits_; + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + if (dEdx_hits_tmp.at(ihit) <= 0.0) dEdx_hits_tmp.at(ihit) = std::numeric_limits::max(); + + for (int ihit = 0; ihit < nhits_sel_; ihit++) { + auto dEdx_min = std::min_element(std::begin(dEdx_hits_tmp), std::end(dEdx_hits_tmp)); + auto dEdx_min_index = std::distance(dEdx_hits_tmp.begin(), dEdx_min); + hits_sel_index_.at(dEdx_min_index) = true; + dEdx_hits_tmp.at(dEdx_min_index) = std::numeric_limits::max(); + } + + indices_are_sel_ = true; +} + +/** +* Calculates energy loss for a track for a selected truncation mode +* @param truncation mode +*/ + +void TrdContainer::CalculateEnergyLossTrack(int trunc_mode) { + + if (nhits_sel_ == 0) nhits_sel_ = GetNHitsSel(trunc_mode); + + std::vector dEdx; + dEdx.clear(); + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + if (dEdx_hits_.at(ihit) > 0) + dEdx.push_back(dEdx_hits_.at(ihit)); + + std::sort(dEdx.begin(), dEdx.end()); + + float dEdx_sum = 0; + for (int ihit = 0; ihit < nhits_sel_; ihit++) { + dEdx_sum += dEdx.at(ihit); + } + dEdx_track_.at(trunc_mode) = dEdx_sum / nhits_sel_; +} + +/** +* Calculates energy loss for a track for a all truncation modes +*/ + +void TrdContainer::CalculateEnergyLossTrackAllModes() { + + std::vector dEdx; + dEdx.clear(); + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + if (dEdx_hits_.at(ihit) > 0) + dEdx.push_back(dEdx_hits_.at(ihit)); + + std::sort(dEdx.begin(), dEdx.end()); + + float dEdx_sum = 0; + for (int ihit = 0; ihit < dEdx.size(); ihit++) { + dEdx_sum += dEdx.at(ihit); + dEdx_track_.at(ihit) = dEdx_sum / (ihit + 1); + } +} diff --git a/src_trd/ContainerTrd.h b/src_trd/ContainerTrd.h new file mode 100644 index 0000000..a75cba9 --- /dev/null +++ b/src_trd/ContainerTrd.h @@ -0,0 +1,83 @@ +/** @file ContainerTrd.h + @author Susanne Glaessel (glaessel@ikf.uni-frankfurt.de) + @brief Class to store Trd PID information +*/ + +#ifndef ContainerTrd_H +#define ContainerTrd_H + +#include "ConstantsTrd.h" + +#include "TMath.h" +#include + +class TrdContainer { + public: + TrdContainer() = default; + + TrdContainer(float mom, float pT, int charge, int nhits_trd, std::array dEdx_hits, std::array dEdx_track, int nhits_sel, std::array hits_sel_index, int mc_pdg, bool dEdxIsScaled = false) : mom_(mom), pT_(pT), charge_(charge), nhits_trd_(nhits_trd), dEdx_hits_(dEdx_hits), dEdx_track_(dEdx_track), nhits_sel_(nhits_sel), hits_sel_index_(hits_sel_index), mc_pdg_(mc_pdg), dEdx_is_scaled_(dEdxIsScaled) {} + + TrdContainer(const float mom, float pT, int charge, int nhits_trd, std::array dEdx_hits, int mc_pdg, bool dEdxIsScaled = false) : mom_(mom), pT_(pT), charge_(charge), nhits_trd_(nhits_trd), dEdx_hits_(dEdx_hits), mc_pdg_(mc_pdg), dEdx_is_scaled_(dEdxIsScaled) { + nhits_sel_ = 0; + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + hits_sel_index_.at(ihit) = false; + } + + TrdContainer(float mom, float pT, int charge, int nhits_trd, std::array dEdx_hits, bool dEdxIsScaled = false) : mom_(mom), pT_(pT), charge_(charge), nhits_trd_(nhits_trd), dEdx_hits_(dEdx_hits), dEdx_is_scaled_(dEdxIsScaled) { + nhits_sel_ = 0; + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) + hits_sel_index_.at(ihit) = false; + mc_pdg_ = -2; + } + + virtual ~TrdContainer() = default; + + void ScaleEnergyLossLength(); + void CalculateEnergyLossTrack(int trunc_mode);// Calculates energy loss of a track for a selected truncation mode + void CalculateEnergyLossTrackAllModes(); // Calculates energy loss of a track for truncation modes 1-4 + void SelectHitIndices(int trunc_mode); // Returns hit indices with dEdx > 0 selected in truncation mode + int GetNHitsSel(int trunc_mode); // Returns number of hits with dEdx > 0 selected in truncation mode + + float GetP() const { return mom_; } + int GetCharge() const { return charge_; } + int GetNhitsTrd() const { return nhits_trd_; } + int GetMcPdg() const { + if (mc_pdg_ == -2) throw std::runtime_error("MC PDG not set."); + else + return mc_pdg_; + } + + std::array GetdEdxHits() const { return dEdx_hits_; } + float GetdEdxTrack(int trunc_mode) { + if (dEdx_track_.at(trunc_mode) == 0) + CalculateEnergyLossTrack(trunc_mode); + return dEdx_track_.at(trunc_mode); + } + std::array GetdEdxHitsSorted();// Orders hits from lowest to highest dEdx (hits with dEdx = 0 in last posisition) + std::array GetHitsSelIndex() const { + if (indices_are_sel_ == false) + throw std::runtime_error("Hits are not selected. Call function SelectHitIndices(trunc_mode) first.\n"); + else + return hits_sel_index_; + } + + protected: + float mom_{0.0}; + float pT_{0.0}; + ; + int charge_{0}; + int nhits_trd_{0}; + int mc_pdg_{-1}; + + std::array dEdx_hits_ = {0.0, 0.0, 0.0, 0.0}; // energy loss of 1-4 hits of a track + std::array dEdx_track_ = {0.0, 0.0, 0.0, 0.0, 0.0};// energy loss of a track for truncation modes 1-4 + + std::array hits_sorted_ = {0.0, 0.0, 0.0, 0.0};// dEdx values of hits sorted from lowest to higest dEdx + int nhits_sel_{0}; + std::array hits_sel_index_ = {false, false, false, false};// hits (in their original) order are selected according to truncation mode + + bool dEdx_is_scaled_{false}; + bool indices_are_sel_{false}; +}; + +#endif//ContainerTrd_HPP diff --git a/src_trd/GetterTrd.cpp b/src_trd/GetterTrd.cpp new file mode 100644 index 0000000..8fc0e8c --- /dev/null +++ b/src_trd/GetterTrd.cpp @@ -0,0 +1,198 @@ +#include "GetterTrd.h" +#include +using std::cout; +using std::endl; +using std::to_string; + +ClassImp(PidTrd::GetterTrd) + + namespace PidTrd { + + /** +* Gets probabilites for a hit or track for a selected truncation and probability mode. +* @param trdtrack +* @param ihit - for likelihood method (= -1 for total probability method) +* @return map with probabilities for a track or hit for all particle species +*/ + std::map GetterTrd::GetTrdProbabilities(TrdContainer trdtrack, int ihit) { + std::map prob; + + float prob_bg = 1.0; + float sum = 0.0; + + int nhits = trdtrack.GetNhitsTrd(); + float mom = trdtrack.GetP(); + float dEdx; + + if (prob_mode_ == 0) { + if (ihit != -1) + throw std::runtime_error("For probability mode 0 (total probability) function GetTrdProbabilities need to be used and ihit should not be set."); + dEdx = trdtrack.GetdEdxTrack(trunc_mode_); + } + if (prob_mode_ == 1) { + if (ihit == -1) + throw std::runtime_error("For probability mode 1 (likelihood) function GetTrdProbabilitiesMulti need to be used."); + dEdx = trdtrack.GetdEdxHits().at(ihit); + } + + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + Int_t ipid_pm = trdtrack.GetCharge() > 0 ? ipid : ipid + NumberOfPidsTrd - 1; + Int_t trunc_mode_getter; + if (trunc_mode_ == 0) trunc_mode_getter = trdtrack.GetNhitsTrd() - 1; + else if (trunc_mode_ > trdtrack.GetNhitsTrd()) + trunc_mode_getter = trdtrack.GetNhitsTrd() - 1;// if (truncation > nhits) truncation = truncation of nhits + else + trunc_mode_getter = trunc_mode_ - 1; + + Int_t id = ipid_pm + (trdtrack.GetNhitsTrd() - 1) * 100 + trunc_mode_getter * 1000 + prob_mode_ * 10000; + prob[pid_codes_trd_.at(ipid).first] = particles_prob_.find(id)->second.Eval(mom, dEdx); + prob_bg -= prob[pid_codes_trd_.at(ipid).first]; + } + + if (prob_bg < 0) prob_bg = 0.0; + prob[PidTrdParticles::kBgPos] = prob_bg; + + return prob; + } + + /** +* Multiplies probabilites of selected hits of a track for a selected truncation and probability mode. +* @param trdtrack +* @return map with probabilities for all particle species for a track +*/ + + std::map GetterTrd::GetTrdProbabilitiesMulti(TrdContainer trdtrack) { + std::map prob; + for (const auto& pdg : pid_codes_trd_) + prob[pdg.first] = 1.0; + + std::map prob_tmp; + + trdtrack.SelectHitIndices(trunc_mode_); + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) { + if (trdtrack.GetHitsSelIndex().at(ihit) == false) continue; + prob_tmp = GetTrdProbabilities(trdtrack, ihit); + + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + prob[pid_codes_trd_.at(ipid).first] *= prob_tmp[pid_codes_trd_.at(ipid).first]; + } + } + + float prob_tot = 0.0; + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) + if (prob[pid_codes_trd_.at(ipid).first] >= 0. && prob[pid_codes_trd_.at(ipid).first] <= 1.) prob_tot += prob[pid_codes_trd_.at(ipid).first]; + + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + if (prob_tot > 0) prob[pid_codes_trd_.at(ipid).first] /= prob_tot; + else + prob[pid_codes_trd_.at(ipid).first] = 0.0; + } + + if (prob_tot > 0) + prob[PidTrdParticles::kBgPos] = 0.0; + else + prob[PidTrdParticles::kBgPos] = 1.0; + + return prob; + } + + /** +* Gets probabilites for a track for the average dEdx with the total probability method. +* @param momentum +* @param average dEdx +* @param number of hits +* @return map with probabilities for a track for all particle species +*/ + + std::map GetterTrd::GetTrdProbabilitiesTotalAverage(float mom, float dEdx, int charge, int nhits, int prob_mode) { + + std::map prob; + float prob_bg = 1.0; + + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + Int_t ipid_pm = charge > 0 ? ipid : ipid + NumberOfPidsTrd - 1; + Int_t trunc_mode_getter = nhits - 1; + + Int_t id = ipid_pm + (nhits - 1) * 100 + trunc_mode_getter * 1000 + prob_mode * 10000; + prob[pid_codes_trd_.at(ipid).first] = particles_prob_.find(id)->second.Eval(mom, dEdx); + prob_bg -= prob[pid_codes_trd_.at(ipid).first]; + } + if (prob_bg < 0) prob_bg = 0.0; + prob[PidTrdParticles::kBgPos] = prob_bg; + return prob; + } + + /** +* Gets probabilites for a track for the average dEdx with the likelihood method. +* @param momentum +* @param vector of dEdx +* @param number of hits +* @return map with probabilities for a track for all particle species + +* Multiplies probabilites of selected hits of a track for a selected truncation and probability mode. +* @param trdtrack +* @return map with probabilities for all particle species for a track +*/ + + std::map GetterTrd::GetTrdProbabilitiesLikelihoodAverage(float mom, std::array dEdx, int charge) { + + std::map prob; + for (const auto& pdg : pid_codes_trd_) + prob[pdg.first] = 1.0; + + std::map prob_tmp; + int nhits = 0; + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) { + if (dEdx.at(ihit) <= 0.0) continue; + nhits++; + } + + for (int ihit = 0; ihit < NumberOfTrdLayers; ihit++) { + if (dEdx.at(ihit) <= 0.0) continue; + prob_tmp = GetTrdProbabilitiesTotalAverage(mom, dEdx.at(ihit), charge, nhits, 1); + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + prob[pid_codes_trd_.at(ipid).first] *= prob_tmp[pid_codes_trd_.at(ipid).first]; + } + } + + float prob_tot = 0.0; + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) + if (prob[pid_codes_trd_.at(ipid).first] >= 0. && prob[pid_codes_trd_.at(ipid).first] <= 1.) prob_tot += prob[pid_codes_trd_.at(ipid).first]; + for (size_t ipid = 0; ipid < NumberOfPidsTrd - 1; ipid++) { + if (prob_tot > 0) prob[pid_codes_trd_.at(ipid).first] /= prob_tot; + else + prob[pid_codes_trd_.at(ipid).first] = 0.0; + } + if (prob_tot > 0) + prob[PidTrdParticles::kBgPos] = 0.0; + else + prob[PidTrdParticles::kBgPos] = 1.0; + return prob; + } + + /** +* Gets pid hypothesis by selecting the particle specie with the highest probability +* @param map with probabilities for all particle species +* @param minium selected purity +* @param charge +* @return pid +*/ + + int GetterTrd::GetTrdPid(std::map prob, float purity, int charge) { + int pid; + + std::array prob_vec; + for (int i = 0; i < NumberOfPidsTrd; i++) + prob_vec[i] = prob[pid_codes_trd_.at(i).first]; + + auto prob_max = std::max_element(std::begin(prob_vec), std::end(prob_vec)); + auto prob_max_index = std::distance(prob_vec.begin(), prob_max); + if (*prob_max >= purity) + pid = TMath::Sign(pid_codes_trd_.at(prob_max_index).first, charge); + else + pid = TMath::Sign(1, charge); + return pid; + } +} diff --git a/src_trd/GetterTrd.h b/src_trd/GetterTrd.h new file mode 100644 index 0000000..c16a90d --- /dev/null +++ b/src_trd/GetterTrd.h @@ -0,0 +1,69 @@ +/** @file GetterTrd.h + @author Susanne Glaessel (glaessel@ikf.uni-frankfurt.de) + @brief Class to calculate Trd-PID probabilities +*/ + +#ifndef PidGetterTrd_H +#define PidGetterTrd_H 1 + +#include "ContainerTrd.h" +#include "ParticleProb.h" +#include "TObject.h" +#include "TString.h" + +namespace PidTrd { + +/** + * @brief Pid Trd getter + */ +class GetterTrd : public TObject { + + public: + GetterTrd() = default; + virtual ~GetterTrd() = default; + + void SetMinHits(int nhits_min) { nhits_min_ = nhits_min; } // Min. number of required hits per track + void SetTruncationMode(int trunc_mode) { trunc_mode_ = trunc_mode; }// Modes for calculating the energy loss for up to 4 layers: + // =0: average over all hits + // =1-4: Select hits with lowest dEdx: + // =1: 1 hit, =2: 2 hits, =3: 3 hits, =4: 4 hits + void SetProbabiltyMode(int prob_mode) { prob_mode_ = prob_mode; } // =0: total probability - probability based on particle multiplicites + // =1: likelihood - probability based on dEdx-distribution of particle species + + int GetMinHits() { return nhits_min_; } + + void AddParticlesProb(std::map particlesprob) { particles_prob_ = particlesprob; } + void AddParticleProb(ParticleProb particleprob) { + int type = particleprob.GetCharge() > 0 ? particleprob.GetType() : particleprob.GetType() + NumberOfPidsTrd - 1; + int nhits = particleprob.GetNhits(); + int truncmode = particleprob.GetTruncMode(); + int probmode = particleprob.GetProbMode(); + particles_prob_[type + nhits * 100 + truncmode * 1000 + probmode * 10000] = particleprob; + } + + std::map GetParticlesProb() { return particles_prob_; } + ParticleProb GetParticleProb(int id) { return particles_prob_[id]; } + ParticleProb GetParticleProb(int type, int charge, int nhits, int truncmode, int probmode) { + if (charge < 0) type += NumberOfPidsTrd - 1; + return particles_prob_[type + nhits * 100 + truncmode * 1000 + probmode * 10000]; + } + + std::map GetTrdProbabilities(TrdContainer trdtrack, int ihit = -1); + std::map GetTrdProbabilitiesMulti(TrdContainer trdtrack); + std::map GetTrdProbabilitiesTotalAverage(float mom, float dEdx, int charge, int nhits, int probmode = 0); + std::map GetTrdProbabilitiesLikelihoodAverage(float mom, std::array dEdx, int charge); + int GetTrdPid(std::map prob, float purity, int charge); + + private: + std::map particles_prob_{}; + + int nhits_min_{1}; + int trunc_mode_{0}; + int prob_mode_{0}; + + array probnames_ = {"probT", "probL"}; + + ClassDef(GetterTrd, 1); +}; +}// namespace PidTrd +#endif//PidTrd_Getter_H diff --git a/src_trd/ParticleProb.cpp b/src_trd/ParticleProb.cpp new file mode 100644 index 0000000..092cb23 --- /dev/null +++ b/src_trd/ParticleProb.cpp @@ -0,0 +1,19 @@ +#include "ParticleProb.h" + +ClassImp(PidTrd::ParticleProb) + + namespace PidTrd { + + /** +* Gets probability depending on momentum and dEdx +* @param momentum +* @param dEdx energy loss in the Trd +* @return probability for a particle specie +*/ + + float ParticleProb::Eval(float mom, float dEdx) { + Int_t binx = hprobabilities_->GetXaxis()->FindBin(mom); + Int_t biny = hprobabilities_->GetYaxis()->FindBin(dEdx); + return hprobabilities_->GetBinContent(binx, biny); + } +} diff --git a/src_trd/ParticleProb.h b/src_trd/ParticleProb.h new file mode 100644 index 0000000..1a1497d --- /dev/null +++ b/src_trd/ParticleProb.h @@ -0,0 +1,61 @@ +/** @file ParticleProb.h + @author Susanne Glaessel (glaessel@ikf.uni-frankfurt.de) + @brief Class to store probabilities for particle species, number of hits, truncation mode and probability mode +*/ + +#ifndef ParticleProb_H +#define ParticleProb_H + +#include "ConstantsTrd.h" + +#include "TFile.h" +#include "TH2F.h" +#include "TString.h" + +using std::array; + +namespace PidTrd { + +class ParticleProb : public TObject { + + public: + ParticleProb() = default; + virtual ~ParticleProb() = default; + + ParticleProb(int type, int charge, int nhits, int truncmode, int probmode, TH2F* hprobabilities) { + particle_type_ = type; + charge_ = charge; + nhits_ = nhits, truncmode_ = truncmode, probmode_ = probmode, hprobabilities_ = hprobabilities; + }; + + void Update(int type, int charge, int nhits, int truncmode, int probmode, TH2F* hprobabilities) { + particle_type_ = type; + charge_ = charge; + nhits_ = nhits, truncmode_ = truncmode, probmode_ = probmode, hprobabilities_ = hprobabilities; + }; + + float Eval(float mom, float dEdx); + + int GetId() { + int type_pm = charge_ > 0 ? particle_type_ : particle_type_ + NumberOfPidsTrd - 1; + return type_pm + nhits_ * 100 + truncmode_ * 1000 + probmode_ * 10000; + }; + int GetType() { return particle_type_; }; + int GetCharge() { return charge_; }; + int GetNhits() { return nhits_; }; + int GetTruncMode() { return truncmode_; }; + int GetProbMode() { return probmode_; }; + TH2F* GetProbabilities() { return hprobabilities_; }; + + private: + int particle_type_{-1}; + int charge_{0}; + int nhits_{0}; + int truncmode_{0}; + int probmode_{0}; + TH2F* hprobabilities_; + + ClassDef(ParticleProb, 1); +}; +};// namespace PidTrd +#endif diff --git a/src_trd/PidTrdLinkDef.h b/src_trd/PidTrdLinkDef.h new file mode 100644 index 0000000..3724f6e --- /dev/null +++ b/src_trd/PidTrdLinkDef.h @@ -0,0 +1,11 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class PidTrd::ParticleProb + ; +#pragma link C++ class PidTrd::GetterTrd + ; + + +#endif diff --git a/tasks/CMakeLists.txt b/tasks/CMakeLists.txt index 38885aa..b08dbb7 100644 --- a/tasks/CMakeLists.txt +++ b/tasks/CMakeLists.txt @@ -1,12 +1,13 @@ # Create a main program using the library -set(EXECUTABLES run_pid_dcm12 run_pid_dcm3) +set(EXECUTABLES run_pid_dcm12 run_pid_dcm3 calculate_mcprob_trd) foreach(EXE ${EXECUTABLES}) add_executable(${EXE} ${EXE}.cpp) - target_link_libraries(${EXE} Pid ${ROOT_LIBRARIES}) - target_include_directories(${EXE} PUBLIC ${CMAKE_SOURCE_DIR}/src) + target_link_libraries(${EXE} Pid PidTrd PidInt ${ROOT_LIBRARIES}) + target_include_directories(${EXE} PUBLIC + ${CMAKE_SOURCE_DIR}/src_tof ${CMAKE_SOURCE_DIR}/src_trd ${CMAKE_SOURCE_DIR}/interface) install(TARGETS ${EXE} RUNTIME DESTINATION bin) endforeach() @@ -16,7 +17,13 @@ if (${PID_BUNDLED_AT}) add_executable(fill_pid fill_pid.cpp) target_link_libraries(fill_pid Pid PidAT AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES}) - target_include_directories(fill_pid PUBLIC ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR}/at_interface ${PROJECT_INCLUDE_DIRECTORIES}) + target_include_directories(fill_pid PUBLIC + ${CMAKE_SOURCE_DIR}/src_tof ${CMAKE_SOURCE_DIR}/src_trd + ${CMAKE_SOURCE_DIR}/at_interface ${PROJECT_INCLUDE_DIRECTORIES}) + + add_executable(create_mcinput_trd create_mcinput_trd.cpp) + target_link_libraries(create_mcinput_trd PidAT AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES}) + target_include_directories(create_mcinput_trd PUBLIC ${CMAKE_SOURCE_DIR}/src_trd ${CMAKE_SOURCE_DIR}/at_interface ${PROJECT_INCLUDE_DIRECTORIES}) install(TARGETS fill_pid RUNTIME DESTINATION bin) endif () diff --git a/tasks/calculate_mcprob_trd.cpp b/tasks/calculate_mcprob_trd.cpp new file mode 100644 index 0000000..c467d1a --- /dev/null +++ b/tasks/calculate_mcprob_trd.cpp @@ -0,0 +1,67 @@ +#include "PidTrdRunGetter.hpp" +#include "TStopwatch.h" +#include +#include + +int calculate_mcprob_trd(const std::string& mcfile_name) { + + Bool_t write_mchistos_out = kTRUE; + + TString macroname = "create_getter"; + + const std::string getter_file = "pid_getter_trd.root"; + const std::string getter_name = "pid_getter_trd"; + + TString outFile; + if (write_mchistos_out == kTRUE) { + TString name = mcfile_name; + name.Replace(name.Index(".root"), 5, ""); + outFile = Form("%s_probabilities.root", name.Data()); + } + + std::cout << std::endl; + std::cout << "-I- " << macroname << ": Using input file " << mcfile_name << std::endl; + + TStopwatch timer; + timer.Start(); + + auto start = std::chrono::system_clock::now(); + ROOT::EnableImplicitMT(4); + + gROOT->SetBatch(true); + gROOT->SetEscape(true); + + auto* pid_trd_rungetter = new PidTrdRunGetter(mcfile_name, getter_file, getter_name); + pid_trd_rungetter->SetWriteMcHistogramsOut(write_mchistos_out); + pid_trd_rungetter->Init(); + pid_trd_rungetter->Exec(); + pid_trd_rungetter->Finish(); + + timer.Stop(); + Double_t rtime = timer.RealTime(); + Double_t ctime = timer.CpuTime(); + std::cout << std::endl + << std::endl; + std::cout << "Macro finished successfully." << std::endl; + std::cout << "Getter is " << getter_file << std::endl; + if (write_mchistos_out == kTRUE) std::cout << "Output file is " << outFile << std::endl; + std::cout << "Real time " << rtime << " s, CPU time " << ctime << "s" + << std::endl + << std::endl; + + return EXIT_SUCCESS; +} + +int main(int argc, char** argv) { + + if (argc < 2) { + std::cout << "Wrong number of arguments! Please use: " << std::endl; + std::cout << " ./create_getter filename_mc\n"; + return EXIT_FAILURE; + } + + const std::string& mcfile_name = argv[1]; + calculate_mcprob_trd(mcfile_name); + + return 0; +} diff --git a/tasks/create_mcinput_trd.cpp b/tasks/create_mcinput_trd.cpp new file mode 100644 index 0000000..07ef5cc --- /dev/null +++ b/tasks/create_mcinput_trd.cpp @@ -0,0 +1,49 @@ +#include + +#include "AnalysisTree/TaskManager.hpp" +#include "PidTrdMcInput.hpp" + +using namespace AnalysisTree; + +int create_mcinput_trd(const std::string& filelist, const std::string& outfilename) { + + auto start = std::chrono::system_clock::now(); + ROOT::EnableImplicitMT(4); + + gROOT->SetBatch(true); + gROOT->SetEscape(true); + + Bool_t update_mchistos = kFALSE; + + auto* pid_trd_mcinput = new PidTrdMcInput(outfilename); + //pid_trd_mcinput->SetRecTracksName("RecTracks"); //analysistree with tof-pid + pid_trd_mcinput->SetRecTracksName("VtxTracks");// analysistree without tof-pid + pid_trd_mcinput->SetSimTracksName("SimParticles"); + pid_trd_mcinput->SetTrdTracksName("TrdTracks"); + pid_trd_mcinput->SetUpdateMcHistos(update_mchistos); + + auto* man = TaskManager::GetInstance(); + man->AddTask(pid_trd_mcinput); + //man->Init({filelist}, {"pTree"}); //analysistree with tof-pid + man->Init({filelist}, {"rTree"});// analysistree without tof-pid + man->Run(-1); // -1 = all events + man->Finish(); + man->ClearTasks(); + + return EXIT_SUCCESS; +} + +int main(int argc, char** argv) { + + if (argc < 3) { + std::cout << "Wrong number of arguments! Please use: " << std::endl; + std::cout << " ./create_mcinput filelist.txt outfilename\n"; + return EXIT_FAILURE; + } + + const std::string& filelist = argv[1]; + const std::string& outfilename = argv[2]; + create_mcinput_trd(filelist, outfilename); + + return 0; +} diff --git a/tasks/fill_pid.cpp b/tasks/fill_pid.cpp index 65848c6..f42c730 100644 --- a/tasks/fill_pid.cpp +++ b/tasks/fill_pid.cpp @@ -4,7 +4,37 @@ using namespace AnalysisTree; -void fill_pid(const std::string& filelist, const std::string& pid_file, const std::string& output) { +void fill_pid(const std::string& filelist, int pid_mode, const std::string& pid_file_tof, const std::string& pid_file_trd, const std::string& output, int truncation_mode, int probability_mode, int trdhits_min, float purity) { + + if (pid_mode == 0) { + std::cout << "\n -I- fill_pid : Run with Tof and Trd pid.\n" + << std::endl; + std::cout << " -I- fill_pid : Using input files " << pid_file_tof << " " << pid_file_trd << std::endl; + std::cout << std::endl; + } + if (pid_mode == 1) { + std::cout << "\n -I- fill_pid : Run with Tof pid only.\n" + << std::endl; + std::cout << " -I- fill_pid : Using input file " << pid_file_tof << std::endl; + std::cout << std::endl; + } + if (pid_mode == 2) { + std::cout << "\n -I- fill_pid : Run with Trd pid only.\n" + << std::endl; + std::cout << " -I- fill_pid : Using input file " << pid_file_trd << std::endl; + std::cout << std::endl; + } + + if (pid_mode == 0 || pid_mode == 2) { + std::cout << "\n************************" << std::endl; + std::cout << "Settings for Trd pid: " << std::endl; + std::cout << " truncation mode : " << truncation_mode << std::endl; + std::cout << " probability mode : " << probability_mode << std::endl; + std::cout << " min hits trd : " << trdhits_min << std::endl; + if (purity != 0.0) std::cout << " min purtiy : " << purity << std::endl; + std::cout << "************************\n" + << std::endl; + } auto* man = TaskManager::GetInstance(); man->SetOutputName(output, "aTree"); @@ -12,12 +42,23 @@ void fill_pid(const std::string& filelist, const std::string& pid_file, const st // The output file will be based on the input one with additional branches (and possibly with removed other branches) man->SetWriteMode(eBranchWriteMode::kCopyTree); - // Branches TrdTracks and RichRings will be removed from the output file, since they are irrelevant for the following analysis. // The branch VtxTracks, which is an input branch, will be also removed, because a new one, based on it, will be created. - man->SetBranchesExclude({"TrdTracks", "RichRings", "VtxTracks"}); + man->SetBranchesExclude({"VtxTracks"}); // Initialize the Pid::Getter prepared at the fitting step: 1-st argument - file, 2-nd argument - name of the getter. - auto* pid_task = new PidFiller(pid_file, "pid_getter"); + auto* pid_task = new PidFiller(pid_file_tof, pid_file_trd, "pid_getter_tof", "pid_getter_trd", pid_mode); + + pid_task->SetRecTracksName("VtxTracks"); + pid_task->SetTofHitsName("TofHits"); + pid_task->SetTrdTracksName("TrdTracks"); + pid_task->SetRichRingsName("RichRings"); + + if (pid_mode == 0 || pid_mode == 2) { + pid_task->SetMinHits(trdhits_min); + pid_task->SetTruncationMode(truncation_mode); + pid_task->SetProbabilityMode(probability_mode); + pid_task->SetPurity(purity); + } man->AddTask(pid_task); @@ -27,16 +68,84 @@ void fill_pid(const std::string& filelist, const std::string& pid_file, const st } int main(int argc, char** argv) { - if (argc <= 2) { - std::cout << "Not enough arguments! Please use:" << std::endl; - std::cout << " ./fill_pid filelist pid_file" << std::endl; + + if (argc < 5) { + std::cout << std::endl; + std::cout << "Wrong number of arguments! Please use: " << std::endl; + std::cout << std::endl; + + std::cout << "For Tof and Trd pid: " << std::endl; + std::cout << " ./fill_pid 0 filelist.txt outputfile pid_file_tof pid_file_trd truncation_mode probability_mode min_hits\n"; + std::cout << "or: " << std::endl; + std::cout << " ./fill_pid 0 filelist.txt outputfile pid_file_tof pid_file_trd truncation_mode probability_mode min_hits purity\n"; + std::cout << std::endl; + + std::cout << "For Tof pid only: " << std::endl; + std::cout << " ./fill_pid 1 filelist.txt outputfile pid_file_tof\n"; + std::cout << std::endl; + + std::cout << "For Trd pid only: " << std::endl; + std::cout << " ./fill_pid 2 filelist.txt outputfile pid_file_trd truncation_mode probability_mode min_hits\n"; + std::cout << "or: " << std::endl; + std::cout << " ./fill_pid 2 filelist.txt outputfile pid_file_trd truncation_mode probability_mode min_hits purity\n"; + std::cout << std::endl; + return EXIT_FAILURE; } - const std::string filelist = argv[1]; - const std::string pid_file = argv[2]; - const std::string output_file = "pid.analysistree.root"; + int pid_mode = atoi(argv[1]); + + if (pid_mode == 0 && argc < 9) { + std::cout << std::endl; + std::cout << "Wrong number of arguments for running Tof and Trd pid simultanously! Please use: " << std::endl; + std::cout << " ./fill_pid 0 filelist.txt outputfile pid_file_tof pid_file_trd truncation_mode probability_mode min_hits\n"; + std::cout << "or: " << std::endl; + std::cout << " ./fill_pid 0 filelist.txt outputfile pid_file_tof pid_file_trd truncation_mode probability_mode min_hits purity\n"; + std::cout << std::endl; + return EXIT_FAILURE; + } + + if (pid_mode == 2 && argc < 8) { + std::cout << std::endl; + std::cout << "Wrong number of arguments for running Trd pid only! Please use: " << std::endl; + std::cout << " ./fill_pid 2 filelist.txt outputfile pid_file_trd truncation_mode probability_mode min_hits\n"; + std::cout << "or: " << std::endl; + std::cout << " ./fill_pid 2 filelist.txt outputfile pid_file_trd truncation_mode probability_mode min_hits purity\n"; + std::cout << std::endl; + return EXIT_FAILURE; + } + + const std::string filelist = argv[2]; + const std::string output_file = argv[3]; + + std::string pid_file_tof = ""; + std::string pid_file_trd = ""; + if (pid_mode == 0) { + pid_file_tof = argv[4]; + pid_file_trd = argv[5]; + } + if (pid_mode == 1) { + pid_file_tof = argv[4]; + } + if (pid_mode == 2) { + pid_file_trd = argv[4]; + } + + int truncation_mode = -1; + int probability_mode = -1; + int trdhits_min = -1; + double purity = 0.0; + if (pid_mode == 0 || pid_mode == 2) { + int i = 0; + if (pid_mode == 0) i = 1; + truncation_mode = atoi(argv[5 + i]); + probability_mode = atoi(argv[6 + i]); + trdhits_min = atoi(argv[7 + i]); + if (argv[8 + i]) purity = atof(argv[8 + i]); + else + purity = 0.0; + } + fill_pid(filelist, pid_mode, pid_file_tof, pid_file_trd, output_file, truncation_mode, probability_mode, trdhits_min, purity); - fill_pid(filelist, pid_file, output_file); return EXIT_SUCCESS; } diff --git a/tasks/run_pid_dcm12.cpp b/tasks/run_pid_dcm12.cpp index 169f84c..f1c93d0 100644 --- a/tasks/run_pid_dcm12.cpp +++ b/tasks/run_pid_dcm12.cpp @@ -1,7 +1,7 @@ -#include "Constants.h" +#include "ConstantsTof.h" #include "FitHelper.h" #include "Fitter.h" -#include "Getter.h" +#include "GetterTof.h" #include "ParticleFit.h" #include "TFile.h" @@ -31,7 +31,7 @@ void run_pid_dcm12() { std::unique_ptr fIn{TFile::Open(inputFileName, "read")}; std::unique_ptr fCuts{TFile::Open(cutsFileName, "read")}; - Pid::Getter getter; + Pid::GetterTof getter; Pid::Fitter tof; float xmin, xmax, ymin, ymax; @@ -43,30 +43,30 @@ void run_pid_dcm12() { // Fit positively charged pions std::cout << "\n\npionpos\n"; - xmin = 0.25, xmax = 10.5, ymin = -2., ymax = 2.; // set ranges for fitting: x - p; y - m2 - tof.SetChi2Max(10000); // set maximal chi2 of fitting to be visualized + xmin = 0.25, xmax = 10.5, ymin = -2., ymax = 2.;// set ranges for fitting: x - p; y - m2 + tof.SetChi2Max(10000); // set maximal chi2 of fitting to be visualized // Set input histogram. // In this task we use the reco_vs_sim_info_via_vtx folder - with "TOF hit - reconstructed track - simulated particles" matching, see ReadMe. std::shared_ptr hpionpos{(TH2D*) fIn->Get("reco_vs_sim_info_via_vtx/h2TofM2_piplus")}; - std::shared_ptr hpionpos_cut{(TH2D*) cutTH2(hpionpos, (TCutG*) fCuts->Get("piplus"))}; // Set graphic cuts for pions. + std::shared_ptr hpionpos_cut{(TH2D*) cutTH2(hpionpos, (TCutG*) fCuts->Get("piplus"))};// Set graphic cuts for pions. hpionpos_cut->Rebin2D(10, 1); - TF1 fit_pionpos("fit_pionpos", "gaus", ymin, ymax); // Fit pions yield as a function of m2 in a certain momentum p range with a gaussian function. - fit_pionpos.SetParNames("p0", "p1", "p2"); // p0 - height in a maximum, p1 - mean value, p2 - sigma. + TF1 fit_pionpos("fit_pionpos", "gaus", ymin, ymax);// Fit pions yield as a function of m2 in a certain momentum p range with a gaussian function. + fit_pionpos.SetParNames("p0", "p1", "p2"); // p0 - height in a maximum, p1 - mean value, p2 - sigma. fit_pionpos.SetParLimits(0, 0., 4.e6); fit_pionpos.SetParLimits(1, -.25, 0.03); fit_pionpos.SetParLimits(2, 0., 0.7); - TF1 pionpos_0("pionpos_0", "0", xmin, xmax); // "0" means do not fit the p0 as a function of momentum p, just remeber all the values. - TF1 pionpos_1("pionpos_1", "pol9", xmin, xmax); // Fit p1 as a function of momentum p with 9-th order polynomial. - TF1 pionpos_2("pionpos_2", "pol9", xmin, xmax); // Fit p2 as a function of momentum p with 9-th order polynomial. + TF1 pionpos_0("pionpos_0", "0", xmin, xmax); // "0" means do not fit the p0 as a function of momentum p, just remeber all the values. + TF1 pionpos_1("pionpos_1", "pol9", xmin, xmax);// Fit p1 as a function of momentum p with 9-th order polynomial. + TF1 pionpos_2("pionpos_2", "pol9", xmin, xmax);// Fit p2 as a function of momentum p with 9-th order polynomial. - Pid::ParticleFit pionpos(PidParticles::kPionPos); + Pid::ParticleFit pionpos(PidTofParticles::kPionPos); pionpos.SetParametrization({pionpos_0, pionpos_1, pionpos_2}); pionpos.SetFitFunction(fit_pionpos); pionpos.SetRange(xmin, xmax); pionpos.SetIsFitted(); - tof.AddParticle(pionpos, PidParticles::kPionPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hpionpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -92,13 +92,13 @@ void run_pid_dcm12() { TF1 kaonpos_1("kaonpos_1", "pol9", xmin, xmax); TF1 kaonpos_2("kaonpos_2", "pol8", xmin, xmax); - Pid::ParticleFit kaonpos(PidParticles::kKaonPos); + Pid::ParticleFit kaonpos(PidTofParticles::kKaonPos); kaonpos.SetParametrization({kaonpos_0, kaonpos_1, kaonpos_2}); kaonpos.SetFitFunction(fit_kaonpos); kaonpos.SetRange(xmin, xmax); kaonpos.SetIsFitted(); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); tof.SetHisto2D(std::move(hkaonpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -124,13 +124,13 @@ void run_pid_dcm12() { TF1 proton_1("proton_1", "pol9", 1.5, 19.); TF1 proton_2("proton_2", "pol12", 1., 20.); - Pid::ParticleFit proton(PidParticles::kProton); + Pid::ParticleFit proton(PidTofParticles::kProton); proton.SetParametrization({proton_0, proton_1, proton_2}); proton.SetFitFunction(fit_proton); proton.SetRange(xmin, xmax); proton.SetIsFitted(); - tof.AddParticle(proton, PidParticles::kProton); + tof.AddParticle(proton, PidTofParticles::kProton); tof.SetHisto2D(std::move(hproton_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -147,7 +147,7 @@ void run_pid_dcm12() { // hhe3_cut->Rebin2D(10,5); // xmin = 1.8, xmax = 8.7, ymin = 0.5, ymax = 3.5; // tof.SetChi2Max(10); - // Pid::ParticleFit he3( PidParticles::kHe3); + // Pid::ParticleFit he3( PidTofParticles::kHe3); // TF1 he3_fit ("fit_he3", "gaus", ymin, ymax); // he3_fit.SetParNames("p9", "p10", "p11"); // he3_fit.SetParLimits(0, 0., 5.e2); @@ -161,7 +161,7 @@ void run_pid_dcm12() { // he3.SetFitFunction( he3_fit ); // he3.SetRange( xmin, xmax ); // he3.SetIsFitted(); - // tof.AddParticle(he3, PidParticles::kHe3); + // tof.AddParticle(he3, PidTofParticles::kHe3); // tof.SetHisto2D( std::move(hhe3_cut) ); // tof.SetRangeX( xmin, xmax ); // tof.SetRangeY( ymin, ymax ); @@ -172,35 +172,36 @@ void run_pid_dcm12() { // End of fit of He-3 (not executed in this example, but can be uncommented) // Fit of deuterons (not executed in this example, but can be uncommented) - // std::cout << "\n\ndeutron\n"; - // xmin = 1.3, xmax = 20., ymin = -3., ymax = 6.; - // tof.SetChi2Max(1000); - // std::unique_ptr hdeutron {(TH2D*) fIn->Get("reco_vs_sim_info_via_vtx/h2TofM2_d")}; - // std::unique_ptr hdeutron_cut {(TH2D*) cutTH2 (hdeutron, (TCutG*) fCuts->Get("deutron"))}; - // hdeutron_cut->Rebin2D(10,25); - // TF1 fit_deutron ("fit_deutron", "gaus", ymin, ymax); - // fit_deutron.SetParNames ("p12", "p13", "p14"); - // fit_deutron.SetParLimits (0, 0., 6.e4); - // fit_deutron.SetParLimits (1, 2.5, 3.6); - // fit_deutron.SetParLimits (2, 0., 2.5); - // TF1 deutron_0 ("deutron_0", "0", xmin, xmax); - // TF1 deutron_1 ("deutron_1", "pol6", xmin, xmax); - // TF1 deutron_2 ("deutron_2", "pol10", xmin, xmax); - // - // Pid::ParticleFit deutron( PidParticles::kProton ); - // deutron.SetParametrization({ deutron_0, deutron_1, deutron_2 }); - // deutron.SetFitFunction(fit_deutron); - // deutron.SetRange( xmin, xmax ); - // deutron.SetIsFitted(); - // - // tof.AddParticle(deutron, PidParticles::kDeutron); - // tof.SetHisto2D( std::move(hdeutron_cut) ); - // tof.SetRangeX( xmin, xmax ); - // tof.SetRangeY( ymin, ymax ); - // tof.SetOutputFileName("deutron.root"); - // tof.Fit(); - // deutron = tof.GetParticle(0); - // tof.Clear(); + std::cout << "\n\ndeutron\n"; + xmin = 1.3, xmax = 20., ymin = -3., ymax = 6.; + tof.SetChi2Max(1000); + std::shared_ptr hdeutron{(TH2D*) fIn->Get("reco_vs_sim_info_via_vtx/h2TofM2_d")}; + std::shared_ptr hdeutron_cut{(TH2D*) cutTH2(hdeutron, (TCutG*) fCuts->Get("deutron"))}; + + hdeutron_cut->Rebin2D(10, 25); + TF1 fit_deutron("fit_deutron", "gaus", ymin, ymax); + fit_deutron.SetParNames("p12", "p13", "p14"); + fit_deutron.SetParLimits(0, 0., 6.e4); + fit_deutron.SetParLimits(1, 2.5, 3.6); + fit_deutron.SetParLimits(2, 0., 2.5); + TF1 deutron_0("deutron_0", "0", xmin, xmax); + TF1 deutron_1("deutron_1", "pol6", xmin, xmax); + TF1 deutron_2("deutron_2", "pol10", xmin, xmax); + + Pid::ParticleFit deutron(PidTofParticles::kProton); + deutron.SetParametrization({deutron_0, deutron_1, deutron_2}); + deutron.SetFitFunction(fit_deutron); + deutron.SetRange(xmin, xmax); + deutron.SetIsFitted(); + + tof.AddParticle(deutron, PidTofParticles::kDeutron); + tof.SetHisto2D(std::move(hdeutron_cut)); + tof.SetRangeX(xmin, xmax); + tof.SetRangeY(ymin, ymax); + tof.SetOutputFileName("deutron.root"); + tof.Fit(); + deutron = tof.GetParticle(0); + tof.Clear(); // End of fit of deuterons (not executed in this example, but can be uncommented) // Here we finished fitting of MC-true distributions fitting and can switch to reconstructed data fitting, using @@ -209,10 +210,10 @@ void run_pid_dcm12() { // Configure fit of background in the positively charged particles side std::cout << "\n\nbgpos\n"; xmin = 0.3, xmax = 20., ymin = -3., ymax = 3.; - Pid::ParticleFit bgpos(PidParticles::kBgPos); - TF1 bgpos_fit("fit_bgpos", "pol2", ymin, ymax); // Fit BG as a function of m2 in a certain momentum p range with a second order polynomial - bgpos_fit.SetParNames("p15", "p16", "p17"); // p15 is a free term, p16 - a linear term, p17 - a square term - TF1 bgpos_0("bgpos_0", "pol3", xmin, xmax); // Fit p15, p16, p17 as functions of momentum p with a third order polynomial + Pid::ParticleFit bgpos(PidTofParticles::kBgPos); + TF1 bgpos_fit("fit_bgpos", "pol2", ymin, ymax);// Fit BG as a function of m2 in a certain momentum p range with a second order polynomial + bgpos_fit.SetParNames("p15", "p16", "p17"); // p15 is a free term, p16 - a linear term, p17 - a square term + TF1 bgpos_0("bgpos_0", "pol3", xmin, xmax); // Fit p15, p16, p17 as functions of momentum p with a third order polynomial TF1 bgpos_1("bgpos_1", "pol3", xmin, xmax); TF1 bgpos_2("bgpos_2", "pol3", xmin, xmax); @@ -224,7 +225,7 @@ void run_pid_dcm12() { std::cout << "\n\nallpos\n"; xmin = 0.3, xmax = 20., ymin = -6., ymax = 6.; - std::unique_ptr hpos{(TH2D*) fIn->Get("reco_info/h2TofM2")}; // Set the input histogram + std::unique_ptr hpos{(TH2D*) fIn->Get("reco_info/h2TofM2")};// Set the input histogram hpos->Rebin2D(10, 1); tof.SetChi2Max(10000); @@ -234,14 +235,14 @@ void run_pid_dcm12() { pionpos.SetIsFixed({false, true, true}); kaonpos.SetIsFixed({false, true, true}); proton.SetIsFixed({false, true, true}); - // deutron.SetIsFixed( {false, true, true} ); + deutron.SetIsFixed({false, true, true}); // he3.SetIsFixed( {false, true, true} ); - tof.AddParticle(pionpos, PidParticles::kPionPos); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); - tof.AddParticle(proton, PidParticles::kProton); - // tof.AddParticle(deutron, PidParticles::kDeutron); - // tof.AddParticle(he3, PidParticles::kHe3); - tof.AddParticle(bgpos, PidParticles::kBgPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); + tof.AddParticle(proton, PidTofParticles::kProton); + tof.AddParticle(deutron, PidTofParticles::kDeutron); + // tof.AddParticle(he3, PidTofParticles::kHe3); + tof.AddParticle(bgpos, PidTofParticles::kBgPos); tof.SetHisto2D(std::move(hpos)); tof.SetRangeX(xmin, xmax); @@ -249,12 +250,12 @@ void run_pid_dcm12() { tof.SetOutputFileName("allpos.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionPos), PidParticles::kPionPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonPos), PidParticles::kKaonPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kProton), PidParticles::kProton); - // getter.AddParticle(tof.GetParticleSpecie(PidParticles::kHe3 ), PidParticles::kHe3 ); - // getter.AddParticle(tof.GetParticleSpecie(PidParticles::kDeutron), PidParticles::kDeutron); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgPos), PidParticles::kBgPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionPos), PidTofParticles::kPionPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonPos), PidTofParticles::kKaonPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kProton), PidTofParticles::kProton); + // getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kHe3 ), PidTofParticles::kHe3 ); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kDeutron), PidTofParticles::kDeutron); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgPos), PidTofParticles::kBgPos); tof.Clear(); // Fit of positively charged particles is finished. @@ -276,13 +277,13 @@ void run_pid_dcm12() { TF1 pionneg_1("pionneg_1", "pol14", -11.7, xmax); TF1 pionneg_2("pionneg_2", "pol12", -11.7, xmax); - Pid::ParticleFit pionneg(PidParticles::kPionNeg); + Pid::ParticleFit pionneg(PidTofParticles::kPionNeg); pionneg.SetParametrization({pionneg_0, pionneg_1, pionneg_2}); pionneg.SetFitFunction(fit_pionneg); pionneg.SetRange(xmin, xmax); pionneg.SetIsFitted(); - tof.AddParticle(pionneg, PidParticles::kPionNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); tof.SetHisto2D(std::move(hpionneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -306,13 +307,13 @@ void run_pid_dcm12() { TF1 kaonneg_1("kaonneg_1", "pol11", xmin, xmax); TF1 kaonneg_2("kaonneg_2", "pol8", xmin, xmax); - Pid::ParticleFit kaonneg(PidParticles::kKaonNeg); + Pid::ParticleFit kaonneg(PidTofParticles::kKaonNeg); kaonneg.SetParametrization({kaonneg_0, kaonneg_1, kaonneg_2}); kaonneg.SetFitFunction(fit_kaonneg); kaonneg.SetRange(xmin, xmax); kaonneg.SetIsFitted(); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); tof.SetHisto2D(std::move(hkaonneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -323,7 +324,7 @@ void run_pid_dcm12() { std::cout << "\n\nbgneg\n"; xmin = -10., xmax = -0.25, ymin = -1., ymax = 2.; - Pid::ParticleFit bgneg(PidParticles::kBgNeg); + Pid::ParticleFit bgneg(PidTofParticles::kBgNeg); TF1 bgneg_0("bgneg_0", "pol3", xmin, xmax);//bgneg_0.SetParameters(100, 0, 0); TF1 bgneg_1("bgneg_1", "pol5", xmin, xmax);//bgneg_1.SetParameters(0, 0, 0); TF1 bgneg_2("bgneg_2", "pol5", xmin, xmax);//bgneg_2.SetParameters(0.0, 0.0, 0); @@ -346,9 +347,9 @@ void run_pid_dcm12() { tof.SetChi2Max(1e6); pionneg.SetIsFixed({false, true, true}); kaonneg.SetIsFixed({false, true, true}); - tof.AddParticle(pionneg, PidParticles::kPionNeg); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); - tof.AddParticle(bgneg, PidParticles::kBgNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); + tof.AddParticle(bgneg, PidTofParticles::kBgNeg); tof.SetHisto2D(std::move(hneg)); tof.SetRangeX(xmin, xmax); @@ -356,16 +357,16 @@ void run_pid_dcm12() { tof.SetOutputFileName("allneg.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionNeg), PidParticles::kPionNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonNeg), PidParticles::kKaonNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgNeg), PidParticles::kBgNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionNeg), PidTofParticles::kPionNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonNeg), PidTofParticles::kKaonNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgNeg), PidTofParticles::kBgNeg); tof.Clear(); // Fit of negatively charged particles finished. - std::unique_ptr outfile{TFile::Open("pid_getter.root", "recreate")}; - getter.Write("pid_getter"); + std::unique_ptr outfile{TFile::Open("pid_getter_tof.root", "recreate")}; + getter.Write("pid_getter_tof"); outfile->Close(); } diff --git a/tasks/run_pid_dcm3.cpp b/tasks/run_pid_dcm3.cpp index 78bad96..76f7323 100644 --- a/tasks/run_pid_dcm3.cpp +++ b/tasks/run_pid_dcm3.cpp @@ -11,10 +11,10 @@ #include #include -#include "Constants.h" +#include "ConstantsTof.h" #include "FitHelper.h" #include "Fitter.h" -#include "Getter.h" +#include "GetterTof.h" #include "ParticleFit.h" void run_pid_dcm3() { @@ -29,7 +29,7 @@ void run_pid_dcm3() { std::shared_ptr fIn{TFile::Open(inputFileName, "read")}; std::shared_ptr fCuts{TFile::Open(cutsFileName, "read")}; - Pid::Getter getter; + Pid::GetterTof getter; Pid::Fitter tof; float xmin, xmax, ymin, ymax; @@ -54,13 +54,13 @@ void run_pid_dcm3() { TF1 pionpos_1("pionpos_1", "pol2", xmin, xmax); TF1 pionpos_2("pionpos_2", "pol5", xmin, xmax); - Pid::ParticleFit pionpos(PidParticles::kPionPos); + Pid::ParticleFit pionpos(PidTofParticles::kPionPos); pionpos.SetParametrization({pionpos_0, pionpos_1, pionpos_2}); pionpos.SetFitFunction(fit_pionpos); pionpos.SetRange(xmin, xmax); pionpos.SetIsFitted(); - tof.AddParticle(pionpos, PidParticles::kPionPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); tof.SetHisto2D(std::move(hpionpos_cut)); tof.SetRangeX(xmin, 8.); tof.SetRangeY(ymin, ymax); @@ -84,13 +84,13 @@ void run_pid_dcm3() { TF1 kaonpos_1("kaonpos_1", "pol2", xmin, xmax); TF1 kaonpos_2("kaonpos_2", "pol6", xmin, xmax); - Pid::ParticleFit kaonpos(PidParticles::kKaonPos); + Pid::ParticleFit kaonpos(PidTofParticles::kKaonPos); kaonpos.SetParametrization({kaonpos_0, kaonpos_1, kaonpos_2}); kaonpos.SetFitFunction(fit_kaonpos); kaonpos.SetRange(xmin, xmax); kaonpos.SetIsFitted(); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); tof.SetHisto2D(std::move(hkaonpos_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -114,13 +114,13 @@ void run_pid_dcm3() { TF1 proton_1("proton_1", "pol11", .3, 12.); TF1 proton_2("proton_2", "pol8", .3, 12.); - Pid::ParticleFit proton(PidParticles::kProton); + Pid::ParticleFit proton(PidTofParticles::kProton); proton.SetParametrization({proton_0, proton_1, proton_2}); proton.SetFitFunction(fit_proton); proton.SetRange(xmin, xmax); proton.SetIsFitted(); - tof.AddParticle(proton, PidParticles::kProton); + tof.AddParticle(proton, PidTofParticles::kProton); tof.SetHisto2D(std::move(hproton_cut)); tof.SetRangeX(xmin, 8.); tof.SetRangeY(ymin, ymax); @@ -135,7 +135,7 @@ void run_pid_dcm3() { // hhe3_cut->Rebin2D(10,5); // xmin = 1.8, xmax = 8.7, ymin = 0.5, ymax = 3.5; // tof.SetChi2Max(10); - // Pid::ParticleFit he3( PidParticles::kHe3); + // Pid::ParticleFit he3( PidTofParticles::kHe3); // TF1 he3_fit ("fit_he3", "gaus", ymin, ymax); // he3_fit.SetParNames("p9", "p10", "p11"); // he3_fit.SetParLimits(0, 0., 5.e2); @@ -149,7 +149,7 @@ void run_pid_dcm3() { // he3.SetFitFunction( he3_fit ); // he3.SetRange( xmin, xmax ); // he3.SetIsFitted(); - // tof.AddParticle(he3, PidParticles::kHe3); + // tof.AddParticle(he3, PidTofParticles::kHe3); // tof.SetHisto2D( std::move(hhe3_cut) ); // tof.SetRangeX( xmin, xmax ); // tof.SetRangeY( ymin, ymax ); @@ -159,39 +159,39 @@ void run_pid_dcm3() { // tof.Clear(); // // - // std::cout << "\n\ndeutron\n"; - // xmin = 1.3, xmax = 20., ymin = -3., ymax = 6.; - // tof.SetChi2Max(1000); - // std::shared_ptr hdeutron {(TH2D*) fIn->Get("reco_vs_sim_info_via_vtx/h2TofM2_d")}; - // std::shared_ptr hdeutron_cut {(TH2D*) cutTH2 (hdeutron, (TCutG*) fCuts->Get("deutron"))}; - // hdeutron_cut->Rebin2D(10,25); - // TF1 fit_deutron ("fit_deutron", "gaus", ymin, ymax); - // fit_deutron.SetParNames ("p12", "p13", "p14"); - // fit_deutron.SetParLimits (0, 0., 6.e4); - // fit_deutron.SetParLimits (1, 2.5, 3.6); - // fit_deutron.SetParLimits (2, 0., 2.5); - // TF1 deutron_0 ("deutron_0", "0", xmin, xmax); - // TF1 deutron_1 ("deutron_1", "pol6", xmin, xmax); - // TF1 deutron_2 ("deutron_2", "pol10", xmin, xmax); - // - // Pid::ParticleFit deutron( PidParticles::kProton ); - // deutron.SetParametrization({ deutron_0, deutron_1, deutron_2 }); - // deutron.SetFitFunction(fit_deutron); - // deutron.SetRange( xmin, xmax ); - // deutron.SetIsFitted(); - // - // tof.AddParticle(deutron, PidParticles::kDeutron); - // tof.SetHisto2D( std::move(hdeutron_cut) ); - // tof.SetRangeX( xmin, xmax ); - // tof.SetRangeY( ymin, ymax ); - // tof.SetOutputFileName("deutron.root"); - // tof.Fit(); - // deutron = tof.GetParticle(0); - // tof.Clear(); + std::cout << "\n\ndeutron\n"; + xmin = 1.3, xmax = 20., ymin = -3., ymax = 6.; + tof.SetChi2Max(1000); + std::shared_ptr hdeutron{(TH2D*) fIn->Get("reco_vs_sim_info_via_vtx/h2TofM2_d")}; + std::shared_ptr hdeutron_cut{(TH2D*) cutTH2(hdeutron, (TCutG*) fCuts->Get("deutron"))}; + hdeutron_cut->Rebin2D(10, 25); + TF1 fit_deutron("fit_deutron", "gaus", ymin, ymax); + fit_deutron.SetParNames("p12", "p13", "p14"); + fit_deutron.SetParLimits(0, 0., 6.e4); + fit_deutron.SetParLimits(1, 2.5, 3.6); + fit_deutron.SetParLimits(2, 0., 2.5); + TF1 deutron_0("deutron_0", "0", xmin, xmax); + TF1 deutron_1("deutron_1", "pol6", xmin, xmax); + TF1 deutron_2("deutron_2", "pol10", xmin, xmax); + + Pid::ParticleFit deutron(PidTofParticles::kProton); + deutron.SetParametrization({deutron_0, deutron_1, deutron_2}); + deutron.SetFitFunction(fit_deutron); + deutron.SetRange(xmin, xmax); + deutron.SetIsFitted(); + + tof.AddParticle(deutron, PidTofParticles::kDeutron); + tof.SetHisto2D(std::move(hdeutron_cut)); + tof.SetRangeX(xmin, xmax); + tof.SetRangeY(ymin, ymax); + tof.SetOutputFileName("deutron.root"); + tof.Fit(); + deutron = tof.GetParticle(0); + tof.Clear(); std::cout << "\n\nbgpos\n"; xmin = 0.3, xmax = 20., ymin = -3., ymax = 3.; - Pid::ParticleFit bgpos(PidParticles::kBgPos); + Pid::ParticleFit bgpos(PidTofParticles::kBgPos); TF1 bgpos_fit("fit_bgpos", "pol2", ymin, ymax); bgpos_fit.SetParNames("p15", "p16", "p17"); TF1 bgpos_0("bgpos_0", "pol3", xmin, xmax); @@ -211,14 +211,14 @@ void run_pid_dcm3() { pionpos.SetIsFixed({false, true, true}); kaonpos.SetIsFixed({false, true, true}); proton.SetIsFixed({false, true, true}); - // deutron.SetIsFixed( {false, true, true} ); + deutron.SetIsFixed({false, true, true}); // he3.SetIsFixed( {false, true, true} ); - tof.AddParticle(pionpos, PidParticles::kPionPos); - tof.AddParticle(kaonpos, PidParticles::kKaonPos); - tof.AddParticle(proton, PidParticles::kProton); - // tof.AddParticle(deutron, PidParticles::kDeutron); - // tof.AddParticle(he3, PidParticles::kHe3); - tof.AddParticle(bgpos, PidParticles::kBgPos); + tof.AddParticle(pionpos, PidTofParticles::kPionPos); + tof.AddParticle(kaonpos, PidTofParticles::kKaonPos); + tof.AddParticle(proton, PidTofParticles::kProton); + tof.AddParticle(deutron, PidTofParticles::kDeutron); + // tof.AddParticle(he3, PidTofParticles::kHe3); + tof.AddParticle(bgpos, PidTofParticles::kBgPos); tof.SetHisto2D(std::move(hpos)); tof.SetRangeX(xmin, xmax); @@ -226,12 +226,12 @@ void run_pid_dcm3() { tof.SetOutputFileName("allpos.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionPos), PidParticles::kPionPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonPos), PidParticles::kKaonPos); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kProton), PidParticles::kProton); - // getter.AddParticle(tof.GetParticleSpecie(PidParticles::kHe3 ), PidParticles::kHe3 ); - // getter.AddParticle(tof.GetParticleSpecie(PidParticles::kDeutron), PidParticles::kDeutron); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgPos), PidParticles::kBgPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionPos), PidTofParticles::kPionPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonPos), PidTofParticles::kKaonPos); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kProton), PidTofParticles::kProton); + // getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kHe3 ), PidTofParticles::kHe3 ); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kDeutron), PidTofParticles::kDeutron); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgPos), PidTofParticles::kBgPos); tof.Clear(); @@ -250,13 +250,13 @@ void run_pid_dcm3() { TF1 pionneg_1("pionneg_1", "pol2", xmin, xmax); TF1 pionneg_2("pionneg_2", "pol5", xmin, xmax); - Pid::ParticleFit pionneg(PidParticles::kPionNeg); + Pid::ParticleFit pionneg(PidTofParticles::kPionNeg); pionneg.SetParametrization({pionneg_0, pionneg_1, pionneg_2}); pionneg.SetFitFunction(fit_pionneg); pionneg.SetRange(xmin, xmax); pionneg.SetIsFitted(); - tof.AddParticle(pionneg, PidParticles::kPionNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); tof.SetHisto2D(std::move(hpionneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -280,13 +280,13 @@ void run_pid_dcm3() { TF1 kaonneg_1("kaonneg_1", "pol2", xmin, xmax); TF1 kaonneg_2("kaonneg_2", "pol2", xmin, xmax); - Pid::ParticleFit kaonneg(PidParticles::kKaonNeg); + Pid::ParticleFit kaonneg(PidTofParticles::kKaonNeg); kaonneg.SetParametrization({kaonneg_0, kaonneg_1, kaonneg_2}); kaonneg.SetFitFunction(fit_kaonneg); kaonneg.SetRange(xmin, xmax); kaonneg.SetIsFitted(); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); tof.SetHisto2D(std::move(hkaonneg_cut)); tof.SetRangeX(xmin, xmax); tof.SetRangeY(ymin, ymax); @@ -297,7 +297,7 @@ void run_pid_dcm3() { std::cout << "\n\nbgneg\n"; xmin = -10., xmax = -0.25, ymin = -1., ymax = 2.; - Pid::ParticleFit bgneg(PidParticles::kBgNeg); + Pid::ParticleFit bgneg(PidTofParticles::kBgNeg); TF1 bgneg_0("bgneg_0", "pol3", xmin, xmax);//bgneg_0.SetParameters(100, 0, 0); TF1 bgneg_1("bgneg_1", "pol5", xmin, xmax);//bgneg_1.SetParameters(0, 0, 0); TF1 bgneg_2("bgneg_2", "pol5", xmin, xmax);//bgneg_2.SetParameters(0.0, 0.0, 0); @@ -320,9 +320,9 @@ void run_pid_dcm3() { tof.SetChi2Max(1e6); pionneg.SetIsFixed({false, true, true}); kaonneg.SetIsFixed({false, true, true}); - tof.AddParticle(pionneg, PidParticles::kPionNeg); - tof.AddParticle(kaonneg, PidParticles::kKaonNeg); - tof.AddParticle(bgneg, PidParticles::kBgNeg); + tof.AddParticle(pionneg, PidTofParticles::kPionNeg); + tof.AddParticle(kaonneg, PidTofParticles::kKaonNeg); + tof.AddParticle(bgneg, PidTofParticles::kBgNeg); tof.SetHisto2D(std::move(hneg)); tof.SetRangeX(xmin, xmax); @@ -330,14 +330,14 @@ void run_pid_dcm3() { tof.SetOutputFileName("allneg.root"); tof.Fit(); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kPionNeg), PidParticles::kPionNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kKaonNeg), PidParticles::kKaonNeg); - getter.AddParticle(tof.GetParticleSpecie(PidParticles::kBgNeg), PidParticles::kBgNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kPionNeg), PidTofParticles::kPionNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kKaonNeg), PidTofParticles::kKaonNeg); + getter.AddParticle(tof.GetParticleSpecie(PidTofParticles::kBgNeg), PidTofParticles::kBgNeg); tof.Clear(); - std::shared_ptr outfile{TFile::Open("pid_getter.root", "recreate")}; - getter.Write("pid_getter"); + std::shared_ptr outfile{TFile::Open("pid_getter_tof.root", "recreate")}; + getter.Write("pid_getter_tof"); outfile->Close(); }