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 }