adding task for subtracting background after jet finding, used for all clustering...
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtCheckDecays.cxx
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>
30 using namespace std;
31
32
33
34 void 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
301 bool 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 }