Rivet analyses referenceMC_HFBRANCHINGMonte Carlo analysis to compute semi-leptonic branching ratios of heavy-flavour hadronsExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Plots to study semi-leptonic decays of heavy-flavour hadrons (branching ratios, hadron and lepton $p_T$. Source code: MC_HFBRANCHING.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/HeavyHadrons.hh"
5#include "Rivet/Math/LorentzTrans.hh"
6
7namespace Rivet {
8
9
10 /// @brief MC analysis to compute semi-leptonic branching ratios of heavy-flavour hadrons
11 class MC_HFBRANCHING : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(MC_HFBRANCHING);
16
17 /// @name Analysis methods
18 /// @{
19
20 const string hadron_id(const int pid) const {
21 switch (pid) {
22 case PID::B0: return "B0";
23 case PID::BPLUS: return "BPLUS";
24 case PID::B0S: return "B0S";
25 case PID::LAMBDAB: return "LAMBDAB";
26 case PID::D0: return "D0";
27 case PID::DPLUS: return "DPLUS";
28 case PID::DSPLUS: return "DSPLUS";
29 case PID::LAMBDACPLUS: return "LAMBDACPLUS";
30 default: return "";
31 }
32 }
33
34 //Semi-leptonic decays with one hadron
35 const vector<int> decay_modes_3body(const int pid) const {
36 switch (pid) {
37 case PID::B0: return {413,411,10413,10411,20413,415};
38 case PID::BPLUS: return {423,421,10423,10421,20423,425};
39 case PID::B0S: return {433,431,10433,10431,20433,435};
40 case PID::LAMBDAB: return {4122,14122,4124};
41 case PID::D0: return {323,321,10323,325,211,213};
42 case PID::DPLUS: return {313,311,10313,315,111,113,221,331,231};
43 case PID::DSPLUS: return {333,221,331,311,313};
44 case PID::LAMBDACPLUS: return {3122,3212,3214,2112,2114};
45 default: return {};
46 }
47 }
48 //Semi-leptonic decays with two hadrons
49 const vector<int> decay_modes_4body(const int pid) const {
50 switch (pid) {
51 case PID::D0: return {321,111,211,311};
52 case PID::DPLUS: return {311,111,321,211};
53 case PID::DSPLUS: return {};
54 case PID::LAMBDACPLUS: return {211,211,111,2112};
55 default: return {};
56 }
57 }
58
59 void fill_Histos(const string &hadron_type, const Particle &p) {
60
61 //********Find decay products of hadron*******//
62
63 vector<int> decay_par; //Vector with direct descendants in hadron decay
64 vector<double> child_pt_LAB, child_pt_COM; //Vector with pT of children from hadron decay
65 decay_par.clear(), decay_par.resize(0); //vector with hadron's direct descendants
66 child_pt_LAB.clear(), child_pt_LAB.resize(0); //vector with pT of hadron's direct descendants
67 child_pt_COM.clear(), child_pt_COM.resize(0); //vector with pT of hadron's direct descendants in hadron's COM
68
69 for( const Particle& child : p.children()){
70 decay_par.push_back(child.abspid()); //Get a list of the direct descendants from the current particle
71 child_pt_LAB.push_back(child.pT()/GeV);
72 // Reset the boost
73 _boost = combine(_boost, _boost.inverse());
74 _boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
75 Particle temp = p;
76 Particle temp_child = child;
77 // temp.setMomentum(_boost.transform(temp.momentum())); // to test that the boost works and gives hadron's pt=0
78 temp_child.setMomentum(_boost.transform(temp_child.momentum())); //transform child's pt in hadron's COM
79
80 child_pt_COM.push_back(temp_child.pT()/GeV);
81 }
82 //remove photons from vector to consider QED radiation effects
83 decay_par.erase(std::remove(decay_par.begin(), decay_par.end(), 22), decay_par.end());
84
85
86 //********Compute branching fractions and fill histograms*******//
87
88 vector<int> Modes_3body = decay_modes_3body(p.abspid()); //vector with decay modes of 3-body decays
89 vector<int> Modes_4body = decay_modes_4body(p.abspid()); //some 4-body decays considered for c-hadrons
90 //Bins to be filled, different for b- and c-hadrons
91 int bin_position = -1;
92 int last_bin_position = -1;
93 bool found_decay_mode = false;
94
95 if (decay_par.size() == 3){ //Semileptonic decays with exactly three decay products
96 int lepton_position = -1; //Lepton position in decay_par vector
97 for (unsigned int i=0; i < Modes_3body.size(); i++){
98 if(contains(decay_par, Modes_3body[i])){
99
100 _h["pt_"+hadron_id(p.abspid())]->fill(p.pT()/GeV); //hadron's pT
101
102 if(hadron_type=="b"){
103 bin_position = 3*i; //For b-hadrons there are e, mu and tau decays
104 last_bin_position = 3*Modes_3body.size()+1;
105 }
106 else if(hadron_type=="c"){
107 bin_position = 2*i; //For c-hadrons there are e and mu decays only
108 last_bin_position = 2*Modes_3body.size()+Modes_4body.size()+1;
109 }
110 else cout<<"I compute decays of heavy-flavour hadrons, pass b or c"<<endl;
111
112 //electron decays
113 if((contains(decay_par, 11) && contains(decay_par, 12))){
114 _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(bin_position+1); //branching fraction
115 //Place the e, mu and tau decays for each mode successively
116 found_decay_mode = true;
117
118 lepton_position = std::find(decay_par.begin(), decay_par.end(), 11)-decay_par.begin();
119 _h[hadron_id(p.abspid())+"_e_pT"] -> fill(child_pt_COM[lepton_position]); //lepton's pT
120 _h[hadron_id(p.abspid())+"_lepton_pT_COM"] -> fill(child_pt_COM[lepton_position]);
121 _h[hadron_id(p.abspid())+"_lepton_pT_LAB"] -> fill(child_pt_LAB[lepton_position]);
122 }
123 //muon decays
124 else if((contains(decay_par, 13) && contains(decay_par, 14))){
125 _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(bin_position+2);
126 found_decay_mode = true;
127 lepton_position = std::find(decay_par.begin(), decay_par.end(), 13)-decay_par.begin();
128 _h[hadron_id(p.abspid())+"_mu_pT"] -> fill(child_pt_COM[lepton_position]);
129 _h[hadron_id(p.abspid())+"_lepton_pT_COM"] -> fill(child_pt_COM[lepton_position]);
130 _h[hadron_id(p.abspid())+"_lepton_pT_LAB"] -> fill(child_pt_LAB[lepton_position]);
131 }
132 //tau decays
133 else if((contains(decay_par, 15) && contains(decay_par, 16))){
134 _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(bin_position+3);
135 found_decay_mode = true;
136 lepton_position = std::find(decay_par.begin(), decay_par.end(), 15)-decay_par.begin();
137 _h[hadron_id(p.abspid())+"_tau_pT"] -> fill(child_pt_COM[lepton_position]);
138 _h[hadron_id(p.abspid())+"_lepton_pT_COM"] -> fill(child_pt_COM[lepton_position]);
139 _h[hadron_id(p.abspid())+"_lepton_pT_LAB"] -> fill(child_pt_LAB[lepton_position]);
140 }
141 }
142 }
143 }
144 else if (decay_par.size() == 4 && hadron_type=="c"){ //Semileptonic decays with exactly four decay products
145 int lepton_position = -1;
146 for (unsigned int i=0; i < Modes_4body.size(); i++){
147 if(contains(decay_par, Modes_4body[i]) && contains(decay_par, Modes_4body[i+1])){
148
149 _h["pt_"+hadron_id(p.abspid())]->fill(p.pT()/GeV); //hadron's pT
150
151 bin_position = 2*Modes_3body.size()+i; //place the 4-body decay modes after the 3-body ones
152 last_bin_position = 2*Modes_3body.size()+Modes_4body.size()+1;
153
154 //electron decays
155 if((contains(decay_par, 11) && contains(decay_par, 12))){
156 _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(bin_position+1);
157 found_decay_mode = true;
158 lepton_position = std::find(decay_par.begin(), decay_par.end(), 11)-decay_par.begin();
159 _h[hadron_id(p.abspid())+"_e_pT"] -> fill(child_pt_COM[lepton_position]); //lepton's pT
160 _h[hadron_id(p.abspid())+"_lepton_pT_COM"] -> fill(child_pt_COM[lepton_position]);
161 _h[hadron_id(p.abspid())+"_lepton_pT_LAB"] -> fill(child_pt_LAB[lepton_position]);
162 }
163 //muon decays
164 else if((contains(decay_par, 13) && contains(decay_par, 14))){
165 _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(bin_position+2);
166 found_decay_mode = true;
167 lepton_position = std::find(decay_par.begin(), decay_par.end(), 13)-decay_par.begin();
168 _h[hadron_id(p.abspid())+"_mu_pT"] -> fill(child_pt_COM[lepton_position]);
169 _h[hadron_id(p.abspid())+"_lepton_pT_COM"] -> fill(child_pt_COM[lepton_position]);
170 _h[hadron_id(p.abspid())+"_lepton_pT_LAB"] -> fill(child_pt_LAB[lepton_position]);
171 }
172 }
173 }
174 }
175 //Fill last bin of branching fractions with decays that don't fall into any category
176 if(!found_decay_mode) _h[hadron_id(p.abspid())+"_frac_clnu"]->fill(last_bin_position);
177 }
178
179 /// Book histograms and initialise projections before the run
180 void init() {
181
182 declare(HeavyHadrons(Cuts::pT > 5.*GeV && Cuts::abseta < 2.5),"HA");
183
184 // histograms
185 //Branching ratios
186 //Include semi-leptonic decay modes only
187 book(_h["B0_frac_clnu"], "BR_B0_clnu", 19, 0.5,19.5);
188 book(_h["B0S_frac_clnu"], "BR_B0S_clnu", 19, 0.5,19.5);
189 book(_h["BPLUS_frac_clnu"], "BR_BPLUS_clnu", 19, 0.5,19.5);
190 book(_h["LAMBDAB_frac_clnu"], "BR_LAMBDAB_clnu", 10, 0.5,10.5);
191 book(_h["D0_frac_clnu"], "BR_D0_clnu", 17, 0.5,17.5);
192 book(_h["DPLUS_frac_clnu"], "BR_DPLUS_clnu", 23, 0.5,23.5);
193 book(_h["DSPLUS_frac_clnu"], "BR_DSPLUS_clnu", 11, 0.5,11.5);
194 book(_h["LAMBDACPLUS_frac_clnu"], "BR_LAMBDACPLUS_clnu", 15, 0.5,15.5);
195
196 //Hadron momentum
197 book(_h["pt_B0" ], "B0_pT", 40, 0., 200.);
198 book(_h["pt_BPLUS" ], "BPLUS_pT", 40, 0., 200.);
199 book(_h["pt_B0S" ], "B0S_pT", 40, 0., 200.);
200 book(_h["pt_D0" ], "D0_pT", 40, 0., 200.);
201 book(_h["pt_DPLUS" ], "DPLUS_pT", 40, 0., 200.);
202 book(_h["pt_DSPLUS" ], "DSPLUS_pT", 40, 0., 200.);
203 book(_h["pt_LAMBDAB" ], "LAMBDAB_pT", 40, 0., 200.);
204 book(_h["pt_LAMBDACPLUS"], "LAMBDACPLUS_pT", 40, 0., 200.);
205
206 //Chared lepton momentum in LAB
207 book(_h["B0_lepton_pT_LAB" ], "B0_lepton_pT_LAB", 20, 0., 100.);
208 book(_h["BPLUS_lepton_pT_LAB" ], "BPLUS_lepton_pT_LAB", 20, 0., 100.);
209 book(_h["B0S_lepton_pT_LAB" ], "B0S_lepton_pT_LAB", 20, 0., 100.);
210 book(_h["D0_lepton_pT_LAB" ], "D0_lepton_pT_LAB", 20, 0., 100.);
211 book(_h["DPLUS_lepton_pT_LAB" ], "DPLUS_lepton_pT_LAB", 20, 0., 100.);
212 book(_h["DSPLUS_lepton_pT_LAB" ], "DSPLUS_lepton_pT_LAB", 20, 0., 100.);
213 book(_h["LAMBDAB_lepton_pT_LAB" ], "LAMBDAB_lepton_pT_LAB", 20, 0., 100.);
214 book(_h["LAMBDACPLUS_lepton_pT_LAB"], "LAMBDACPLUS_lepton_pT_LAB", 20, 0., 100.);
215
216
217 //Chared lepton momentum in COM
218 book(_h["B0_e_pT"], "B0_e_pT", 25, 0., 2.5);
219 book(_h["B0_mu_pT"], "B0_mu_pT", 25, 0., 2.5);
220 book(_h["B0_tau_pT"], "B0_tau_pT", 25, 0., 2.5);
221 book(_h["B0_lepton_pT_COM"], "B0_lepton_pT_COM", 25, 0., 2.5);
222
223 book(_h["B0S_e_pT"], "B0S_e_pT", 25, 0., 2.5);
224 book(_h["B0S_mu_pT"], "B0S_mu_pT", 25, 0., 2.5);
225 book(_h["B0S_tau_pT"], "B0S_tau_pT", 25, 0., 2.5);
226 book(_h["B0S_lepton_pT_COM"], "B0S_lepton_pT_COM", 25, 0., 2.5);
227
228 book(_h["BPLUS_e_pT"], "BPLUS_e_pT", 25, 0., 2.5);
229 book(_h["BPLUS_mu_pT"], "BPLUS_mu_pT", 25, 0., 2.5);
230 book(_h["BPLUS_tau_pT"], "BPLUS_tau_pT", 25, 0., 2.5);
231 book(_h["BPLUS_lepton_pT_COM"], "BPLUS_lepton_pT_COM", 25, 0., 2.5);
232
233 book(_h["LAMBDAB_e_pT"], "LAMBDAB_e_pT", 25, 0., 2.5);
234 book(_h["LAMBDAB_mu_pT"], "LAMBDAB_mu_pT", 25, 0., 2.5);
235 book(_h["LAMBDAB_tau_pT"], "LAMBDAB_tau_pT", 25, 0., 2.5);
236 book(_h["LAMBDAB_lepton_pT_COM"], "LAMBDAB_lepton_pT_COM", 25, 0., 2.5);
237
238 book(_h["D0_e_pT"], "D0_e_pT", 15, 0., 1.5);
239 book(_h["D0_mu_pT"], "D0_mu_pT", 15, 0., 1.5);
240 book(_h["D0_lepton_pT_COM"], "D0_lepton_pT_COM", 15, 0., 1.5);
241
242 book(_h["DSPLUS_e_pT"], "DSPLUS_e_pT", 15, 0., 1.5);
243 book(_h["DSPLUS_mu_pT"], "DSPLUS_mu_pT", 15, 0., 1.5);
244 book(_h["DSPLUS_lepton_pT_COM"], "DSPLUS_lepton_pT_COM", 15, 0., 1.5);
245
246 book(_h["DPLUS_e_pT"], "DPLUS_e_pT", 15, 0., 1.5);
247 book(_h["DPLUS_mu_pT"], "DPLUS_mu_pT", 15, 0., 1.5);
248 book(_h["DPLUS_lepton_pT_COM"], "DPLUS_lepton_pT_COM", 15, 0., 1.5);
249
250 book(_h["LAMBDACPLUS_e_pT"], "LAMBDACPLUS_e_pT", 25, 0., 2.5);
251 book(_h["LAMBDACPLUS_mu_pT"], "LAMBDACPLUS_mu_pT", 25, 0., 2.5);
252 book(_h["LAMBDACPLUS_lepton_pT_COM"], "LAMBDACPLUS_lepton_pT_COM", 25, 0., 2.5);
253 }
254
255
256 /// Perform the per-event analysis
257 void analyze(const Event& event) {
258
259 const HeavyHadrons &ha = apply<HeavyHadrons>(event,"HA");
260
261 if (ha.bHadrons().empty() && ha.cHadrons().empty()) vetoEvent;
262
263 //b hadrons branching ratios
264 for (const Particle &hadron : ha.bHadrons()) {
265 //Compute branching ratios for listed b-hadrons
266 if (hadron_id(hadron.abspid())!="") fill_Histos("b", hadron);
267 }
268 //c hadrons branching ratios
269 for (const Particle &hadron : ha.cHadrons()) {
270 if( !hadron.fromBottom()){ //take into account only c-hadrons that don't come from a b-hadron decay
271 //Compute branching ratios for listed c-hadrons
272 if (hadron_id(hadron.abspid())!="") fill_Histos("c", hadron);
273 }
274 }
275 } // Close Event
276
277 /// Normalise histograms etc., after the run
278 void finalize() {
279 normalize(_h);
280 }
281
282 /// @}
283
284
285 private:
286 /// @name Histograms
287 /// @{
288 map<string, Histo1DPtr> _h;
289 /// @}
290 LorentzTransform _boost;
291
292 };
293
294
295 RIVET_DECLARE_PLUGIN(MC_HFBRANCHING);
296
297}
|