adding task for subtracting background after jet finding, used for all clustering...
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtCheckDecays.cxx
CommitLineData
da0e9ce3 1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed
5// for the BaBar collaboration. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Module: EvtCheckDecays
9//
10// Description: Holds code to conduct various checks on the
11// EvtDecayTable::decaytable()
12//
13// Modification history:
14// Abi Soffer Nov 29, 2007, created
15//
16//------------------------------------------------------------------------
17//
18#include "EvtGenBase/EvtPatches.hh"
19
20#include "EvtGen/EvtCheckDecays.hh"
21
22#include "EvtGenBase/EvtDecayTable.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtDecayBase.hh"
25#include "EvtGenBase/EvtReport.hh"
26#include "EvtGenBase/EvtId.hh"
27#include <vector>
28#include <string.h>
29#include <iostream>
30using namespace std;
31
32
33
34void EvtCheckDecays::checkConj(bool compareArguments) {
35 std::ostream & str = report(INFO, "EvtGen: checkConj")
36 << "Checking CP-conjugation of all decays:" << endl;
37
38 const std::vector<EvtParticleDecayList> & decaytable
39 = EvtDecayTable::decaytable();
40
41 vector<unsigned int> skipEntries;
42
43 // Loop over the entris (=particles) in the decay table:
44 for(size_t e = 0; e < decaytable.size(); ++e){
45 bool diffFound = false;
46
47 // Skip entries that were already found to be conjugates of other
48 // particles:
49 bool skipThis = false;
50 for (size_t s = 0; s < skipEntries.size(); ++s) {
51 if (skipEntries[s] == e) {
52 skipThis = true;
53 break;
54 }
55 }
56 if (skipThis) {
57 continue;
58 }
59
60 // Start working with this particle's decays list:
61 const EvtParticleDecayList & decList = decaytable[e];
62
63 if (decList.getNMode() == 0) { // an undecaying entry in the table. ignore:
64 continue;
65 }
66
67 // Get the decaying particle's ID (somewhat non-intuitively, this
68 // is available from the individual decays, but not the decay list):
69 const EvtId parentId = decList.getDecayModel(0)->getParentId();
70
71 if (selfConj(parentId)) { // ignore self-conjugate particles:
72 continue;
73 }
74
75 // Find the charge conjugate particle by looping over decaytable again:
76 const EvtParticleDecayList * decList2 = 0;
77 bool conjFound = false;
78
79 // Need to create parentId2 out of the loop in which it gets filled witn
80 // meaningful information. Since EVtId has no argumentless constructor,
81 // initialize it to parentId. This won't cause logic problems, since
82 // the variable conjFound is used to determine if parentId2 is meaningful:
83 EvtId parentId2 = parentId;
84
85 // the loop starts at e + 1, so as not to double count:
86 for(size_t e2 = e + 1; e2 < EvtPDL::entries(); ++e2){
87 decList2 = &(decaytable[e2]);
88
89 if (decList2->getNMode() == 0) { // an undecaying entry. ignore:
90 continue;
91 }
92
93 parentId2 = decList2->getDecayModel(0)->getParentId();
94 if (parentId2.isConjugate(parentId)) {
95 conjFound = true;
96 skipEntries.push_back(e2);
97 break; // found the conjugate. decList2 is the conj of decList
98 }
99 }
100
101 // Check if conjugate was found:
102 if (false == conjFound) {
103 report(INFO, "EvtGen: checkConj")
104 << "No conjugate particle decays found for " << EvtPDL::name(parentId)
105 << endl;
106 }
107 else { // conjugate found.
108 const std::string name = EvtPDL::name(parentId);
109 const std::string name2 = EvtPDL::name(parentId2);
110
111 // Analyze the two decays and compare them:
112 // First, compare the number of decay modes:
113 if (decList.getNMode() != decList2->getNMode()) {
114 report(INFO, "EvtGen: checkConj")
115 << "*** Unequal numbers of decays: "
116 << decList.getNMode() << " for " << name
117 << " and " << decList2->getNMode() << " for " << name2
118 << "." << endl;
119
120 // Find the particle with the more decays:
121 const EvtParticleDecayList * decListLonger = &decList;
122 string nameLonger = name;
123 string nameShorter = name2;
124 if (decList.getNMode() < decList2->getNMode()) {
125 decListLonger = decList2;
126 nameLonger = name2;
127 nameShorter = name;
128 }
129
130 const int ndiff = abs(decList.getNMode() - decList2->getNMode());
131 str << " The last " << ndiff << " " << nameLonger
132 << " decays (the ones missing from " << nameShorter
133 << ", assuming" << endl
134 << " there are no other problems) are " << endl;
135
136 for (int m = decListLonger->getNMode() - ndiff;
137 m < decListLonger->getNMode();
138 ++m) {
139 const EvtDecayBase * dec = decListLonger->getDecayModel(m);
140 str << " ";
141 dec->printInfo();
142 }
143
144 diffFound = true;
145 }
146
147 // Compare each decay mode. Expect them to be entered in the same order:
148 for (int m=0; m < decList.getNMode() && m < decList2->getNMode(); ++m) {
149 const EvtDecayBase * dec = decList.getDecayModel(m);
150 const EvtDecayBase * dec2 = decList2->getDecayModel(m);
151
152 // Compare BRs:
153 if (dec->getBranchingFraction() != dec2->getBranchingFraction()) {
154 report(INFO, "EvtGen: checkConj")
155 << "*** Unequal BRs: " << endl << " "
156 << dec->getBranchingFraction() << " for decay ";
157 dec->printInfo();
158 str << "O and " << endl << " "
159 << dec2->getBranchingFraction() << " for decay ";
160 dec2->printInfo();
161 str << endl;
162 diffFound = true;
163 }
164
165 // Compare number of daughters:
166 if (dec->getNDaug() != dec2->getNDaug()) {
167 report(INFO, "EvtGen: checkConj")
168 << "*** Unequal #'s of daughters: " << endl << " "
169 << dec->getNDaug() << " for decay ";
170 dec->printInfo();
171 str << " and " << endl << " "
172 << dec2->getNDaug() << " for decay ";
173 dec2->printInfo();
174 str << endl;
175 diffFound = true;
176 }
177
178
179 // Compare daughter ID's. Here we allow for the daughter order
180 // not to match, since it's a common practice to write DECAY.DEC
181 // without much regard to carefully conjugating the daughter order:
182
183 vector<int> usedDtrIndices; // to avoid double-considering a dtr2
184 for (int d = 0; d < dec->getNDaug(); ++d) { // loop on dtrs of dec
185 const EvtId dtr = dec->getDaug(d);
186 const bool sConj = selfConj(dtr);
187
188 bool dtrMatchFound = false;
189
190 for (int d2 = 0; d2 < dec2->getNDaug(); ++d2) { // on dtrs of dec2
191 // Skip if this index has been already found for this decay:
192 bool skipDtr = false;
193 for (size_t s = 0; s < usedDtrIndices.size(); ++s) {
194 if (usedDtrIndices[s] == d2) {
195 skipDtr = true;
196 break;
197 }
198 if (skipDtr) {
199 continue; // to next value of d2
200 }
201 } // end loop on s
202
203 const EvtId dtr2 = dec2->getDaug(d2);
204 // See if there is a matching daughter: Either dtr is
205 // self-conjugate and we find an identical particle among
206 // dec2's daughters, or it is no self-conjugate and we
207 // find its conjugate:
208 if (dtr.isConjugate(dtr2) && !sConj ||
209 dec->getDaug(d) == dec2->getDaug(d2) && sConj) {
210 dtrMatchFound = true;
211 usedDtrIndices.push_back(d2); // don't consider this d2 again
212 break; // out of d2 loop
213 }
214 } // end d2 loop
215
216 if (false == dtrMatchFound) {
217 // Couldn't find a match for dtr among the daughters of dec2, so:
218 report(INFO, "EvtGen: checkConj")
219 << "*** Daughter #" << d << " in decay"
220 << endl << " ";
221 dec->printInfo();
222 str << " has no conjugate in decay " << endl << " ";
223 dec2->printInfo();
224 str << endl;
225 diffFound = true;
226 }
227 }
228
229 // Compare models:
230 if (dec->getModelName() != dec2->getModelName()) {
231 report(INFO, "EvtGen: checkConj")
232 << "*** Unequal model names in decays: " << endl << " ";
233 dec->printInfo();
234 str << " and " << endl << " ";
235 dec2->printInfo();
236 str << endl;
237 diffFound = true;
238 }
239
240 // Compare numbers of arguments:
241 if (dec->getNArg() != dec2->getNArg()) {
242 report(INFO, "EvtGen: checkConj")
243 << "*** Unequal numbers of arguments: " << endl << " "
244 << dec->getNArg() << " in decay ";
245 dec->printInfo();
246 str << " and " << endl << " "
247 << dec2->getNArg() << " in decay ";
248 dec2->printInfo();
249 str << endl;
250 diffFound = true;
251 }
252
253 // Argument value comparison may not always be desired, since
254 // the argument values aren't always stored by EvtDecayBase, causing
255 // many printouts:
256 if (compareArguments) {
257 // Compare argument values. First, check if we can tell
258 // the numbers of arguments correctly:
259 if ( dec->getNArg() != dec->getNStoredArg() ||
260 dec2->getNArg() != dec2->getNStoredArg()) {
261 report(INFO, "EvtGen: checkConj")
262 << "=== Not sure about number of arguments in decay"
263 << endl << " ";
264 dec->printInfo();
265 str << " and " << endl << " ";
266 dec2->printInfo();
267 str << " Argument value comparison may not be complete." << endl;
268 }
269
270 // Now do the actual argument values comparison:
271 for (int a = 0; a < dec->getNArg() && a < dec2->getNArg()
272 && a < dec->getNStoredArg() && a < dec2->getNStoredArg();
273 ++a) {
274 if (dec->getStoredArg(a) != dec2->getStoredArg(a)) {
275 report(INFO, "EvtGen: checkConj")
276 << "*** Unequal arguments #: " << a << ":"
277 << endl << " "
278 << dec->getStoredArg(a) << " for decay ";
279 dec->printInfo();
280 str << " and " << endl << " "
281 << dec2->getStoredArg(a) << " for decay ";
282 dec2->printInfo();
283 str << endl;
284 diffFound = true;
285 } // end if arguments unequal
286 } // end loop on arguments
287 } // end if to compare arguments
288 } // end loop on decays
289 } // end if conjugate found
290
291 if (diffFound) { // mark diff between particle decays
292 report(INFO, "EvtGen: checkConj")
293 << "-------------------------------------------------------" << endl;
294 }
295
296 } // end loop on particle entries
297}
298
299
300
301bool EvtCheckDecays::selfConj(const EvtId & id) {
302 string name = EvtPDL::name(id);
303 if (name == "vpho" ||
304 name == "gamma" ||
305 name == "g" ||
306 name == "specflav" ||
307 name == "phasespa" ||
308 name == "pi0" ||
309 name == "pi(2S)0" ||
310 name == "eta" ||
311 name == "eta(2S)" ||
312 name == "eta(1405)" ||
313 name == "eta(1475)" ||
314 name == "phi(1680)" ||
315 name == "eta'" ||
316 name == "rho0" ||
317 name == "rho(2S)0" ||
318 name == "rho(3S)0" ||
319 name == "omega" ||
320 name == "phi" ||
321 name == "a_00" ||
322 name == "f_0" ||
323 name == "f'_0" ||
324 name == "b_10" ||
325 name == "h_1" ||
326 name == "h'_1" ||
327 name == "a_10" ||
328 name == "f_1" ||
329 name == "f'_1" ||
330 name == "a_10" ||
331 name == "a_20" ||
332 name == "f_2" ||
333 name == "f_0(1500)" ||
334 name == "f'_2" ||
335 name == "K_S0" ||
336 name == "K_L0" ||
337 name == "eta_c" ||
338 name == "eta_c(2S)" ||
339 name == "J/psi" ||
340 name == "psi(2S)" ||
341 name == "psi(3770)" ||
342 name == "psi(4040)" ||
343 name == "psi(4160)" ||
344 name == "psi(4415)" ||
345 name == "h_c" ||
346 name == "chi_c0" ||
347 name == "chi_c1" ||
348 name == "chi_c2" ||
349 name == "eta_b" ||
350 name == "eta_b(2S)" ||
351 name == "eta_b(3S)" ||
352 name == "Upsilon" ||
353 name == "Upsilon(2S)" ||
354 name == "Upsilon(3S)" ||
355 name == "Upsilon(4S)" ||
356 name == "Upsilon(5S)" ||
357 name == "h_b" ||
358 name == "h_b(2P)" ||
359 name == "h_b(3P)" ||
360 name == "chi_b0" ||
361 name == "chi_b1" ||
362 name == "chi_b2" ||
363 name == "chi_b0(2P)" ||
364 name == "chi_b1(2P)" ||
365 name == "chi_b2(2P)" ||
366 name == "chi_b0(3P)" ||
367 name == "chi_b1(3P)" ||
368 name == "chi_b2(3P)" ||
369 name == "eta_b2(1D)" ||
370 name == "eta_b2(2D)" ||
371 name == "Upsilon_1(1D)" ||
372 name == "Upsilon_2(1D)" ||
373 name == "Upsilon_3(1D)" ||
374 name == "Upsilon_1(2D)" ||
375 name == "Upsilon_2(2D)" ||
376 name == "Upsilon_3(2D)" ||
377 name == "Xu0" ||
378 name == "sigma_0") {
379 return true;
380 }
381
382 return false;
383}