Line | Branch | Exec | Source |
---|---|---|---|
1 | #include "Validator.h" | ||
2 | |||
3 | #include <TProfile.h> | ||
4 | #include <TStyle.h> | ||
5 | |||
6 | namespace iguana::clas12 { | ||
7 | |||
8 | REGISTER_IGUANA_VALIDATOR(MomentumCorrectionValidator); | ||
9 | |||
10 | 1 | void MomentumCorrectionValidator::Start(hipo::banklist& banks) | |
11 | { | ||
12 | // define the algorithm sequence | ||
13 | 2 | m_algo_seq = std::make_unique<AlgorithmSequence>(); | |
14 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | m_algo_seq->Add("clas12::EventBuilderFilter"); |
15 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | m_algo_seq->Add("clas12::SectorFinder"); |
16 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | m_algo_seq->Add("clas12::MomentumCorrection"); |
17 |
4/10✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
3 | m_algo_seq->SetOption<std::vector<int>>("clas12::EventBuilderFilter", "pids", u_pdg_list); |
18 | 1 | m_algo_seq->Start(banks); | |
19 | |||
20 | // get bank indices | ||
21 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | b_particle = GetBankIndex(banks, "REC::Particle"); |
22 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | b_sector = GetBankIndex(banks, "REC::Particle::Sector"); |
23 | |||
24 | // set an output file | ||
25 | 1 | auto output_dir = GetOutputDirectory(); | |
26 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | if(output_dir) { |
27 |
3/8✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2 | m_output_file_basename = output_dir.value() + "/momentum_corrections"; |
28 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | m_output_file = new TFile(m_output_file_basename + ".root", "RECREATE"); |
29 | } | ||
30 | |||
31 | // define plots | ||
32 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | gStyle->SetOptStat(0); |
33 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | for(auto const& pdg : u_pdg_list) { |
34 | std::vector<TH2D*> deltaPvsP; | ||
35 |
2/6✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
4 | TString particle_name = particle::name.at(particle::PDG(pdg)); |
36 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | TString particle_title = particle::title.at(particle::PDG(pdg)); |
37 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 4 times.
|
28 | for(int sec = 1; sec <= 6; sec++) { |
38 |
2/4✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
|
24 | TString sector_name = Form("sec%d", sec); |
39 |
2/4✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
|
24 | TString sector_title = Form("sector %d", sec); |
40 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | deltaPvsP.push_back(new TH2D( |
41 |
3/6✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
|
48 | "deltaPvsP_" + particle_name + "_" + sector_name, |
42 |
2/6✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
48 | particle_title + " momentum correction, " + sector_title + ";p [GeV];#Delta p [GeV]", |
43 | 24 | 30, 0, m_p_max, | |
44 |
2/4✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
|
48 | 100, -m_deltaP_max, m_deltaP_max)); |
45 | 24 | } | |
46 | 4 | u_deltaPvsP.insert({pdg, deltaPvsP}); | |
47 | 4 | } | |
48 | 1 | } | |
49 | |||
50 | |||
51 | 1000 | void MomentumCorrectionValidator::Run(hipo::banklist& banks) const | |
52 | { | ||
53 | // get the momenta before | ||
54 |
1/2✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
|
1000 | auto& particle_bank = GetBank(banks, b_particle, "REC::Particle"); |
55 |
2/4✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1000 times.
✗ Branch 3 not taken.
|
2000 | auto& sector_bank = GetBank(banks, b_sector, "REC::Particle::Sector"); |
56 | std::vector<double> p_measured; | ||
57 |
3/4✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 7143 times.
✓ Branch 3 taken 1000 times.
|
8143 | for(auto const& row : particle_bank.getRowList()) |
58 |
1/4✓ Branch 0 taken 7143 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
14286 | p_measured.push_back(std::hypot( |
59 | particle_bank.getFloat("px", row), | ||
60 | particle_bank.getFloat("py", row), | ||
61 | particle_bank.getFloat("pz", row))); | ||
62 | |||
63 | // run the momentum corrections | ||
64 |
1/2✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
|
1000 | m_algo_seq->Run(banks); |
65 | |||
66 | // lock the mutex, so we can mutate plots | ||
67 |
1/2✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
|
1000 | std::scoped_lock<std::mutex> lock(m_mutex); |
68 | |||
69 | // fill the plots | ||
70 |
3/4✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2348 times.
✓ Branch 3 taken 1000 times.
|
3348 | for(auto const& row : particle_bank.getRowList()) { |
71 | |||
72 | 2348 | auto pdg = particle_bank.getInt("pid", row); | |
73 | 2348 | auto sector = sector_bank.getInt("sector", row); | |
74 | |||
75 | // skip central particle, or unknown sector | ||
76 |
2/2✓ Branch 0 taken 878 times.
✓ Branch 1 taken 1470 times.
|
2348 | if(sector == 0) |
77 | 878 | continue; | |
78 | |||
79 | 1470 | double p_corrected = std::hypot( | |
80 | particle_bank.getFloat("px", row), | ||
81 | particle_bank.getFloat("py", row), | ||
82 | 1470 | particle_bank.getFloat("pz", row)); | |
83 |
2/4✓ Branch 0 taken 1470 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1470 times.
✗ Branch 3 not taken.
|
1470 | auto delta_p = p_corrected - p_measured.at(row); |
84 |
2/4✓ Branch 0 taken 1470 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1470 times.
✗ Branch 3 not taken.
|
1470 | u_deltaPvsP.at(pdg).at(sector - 1)->Fill(p_corrected, delta_p); |
85 | } | ||
86 | 1000 | } | |
87 | |||
88 | |||
89 | 1 | void MomentumCorrectionValidator::Stop() | |
90 | { | ||
91 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
2 | if(GetOutputDirectory()) { |
92 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | for(auto const& [pdg, plots] : u_deltaPvsP) { |
93 | int n_cols = 3; | ||
94 | int n_rows = 2; | ||
95 | 4 | TString canv_name = Form("canv%d", pdg); | |
96 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | auto canv = new TCanvas(canv_name, canv_name, n_cols * 800, n_rows * 600); |
97 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | canv->Divide(n_cols, n_rows); |
98 | int pad_num = 0; | ||
99 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 4 times.
|
28 | for(auto const& plot : plots) { |
100 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | auto pad = canv->GetPad(++pad_num); |
101 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->cd(); |
102 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->SetGrid(1, 1); |
103 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->SetLogz(); |
104 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->SetLeftMargin(0.12); |
105 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->SetRightMargin(0.12); |
106 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | pad->SetBottomMargin(0.12); |
107 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | plot->Draw("colz"); |
108 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | plot->GetYaxis()->SetRangeUser(-m_deltaP_zoom, m_deltaP_zoom); |
109 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | auto prof = plot->ProfileX("_pfx", 1, -1, "s"); |
110 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | prof->SetLineColor(kBlack); |
111 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | prof->SetLineWidth(5); |
112 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | prof->Draw("same"); |
113 | } | ||
114 |
5/10✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
|
8 | canv->SaveAs(m_output_file_basename + "_" + std::to_string(pdg) + ".png"); |
115 | 4 | } | |
116 | 1 | m_output_file->Write(); | |
117 | 1 | m_log->Info("Wrote output file {}", m_output_file->GetName()); | |
118 | 1 | m_output_file->Close(); | |
119 | } | ||
120 | 1 | } | |
121 | |||
122 | } | ||
123 |