TEvtGen/EvtGenModels/EvtDalitzTable.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtGen/EvtGenericDalitz.hh
  12 //
  13 // Description: Model to describe a generic dalitz decay
  14 //
  15 // Modification history:
  16 //
  17 //    DCC     16 December, 2011         Module created
  18 //
  19 //------------------------------------------------------------------------
  20
  21 #include "EvtGenModels/EvtDalitzTable.hh"
  22 #include "EvtGenBase/EvtReport.hh"
  23 #include "EvtGenBase/EvtParserXml.hh"
  24 #include "EvtGenBase/EvtPDL.hh"
  25 #include "EvtGenBase/EvtSpinType.hh"
  26 #include "EvtGenBase/EvtDalitzPlot.hh"
  27 #include "EvtGenBase/EvtCyclic3.hh"
  28
  29 #include <stdlib.h>
  30 #include <sstream>
  31
  32 using std::endl;
  33 using std::fstream;
  34 using std::ifstream;
  35
  36 EvtDalitzTable::EvtDalitzTable() {
  37   _dalitztable.clear();
  38   _readFiles.clear();
  39 }
  40
  41 EvtDalitzTable::~EvtDalitzTable() {
  42   _dalitztable.clear();
  43   _readFiles.clear();
  44 }
  45
  46 EvtDalitzTable* EvtDalitzTable::getInstance(const std::string dec_name, bool verbose) {
  47
  48   static EvtDalitzTable* theDalitzTable = 0;
  49
  50   if(theDalitzTable == 0) {
  51     theDalitzTable = new EvtDalitzTable();
  52   }
  53
  54   if(!theDalitzTable->fileHasBeenRead(dec_name)) {
  55     theDalitzTable->readXMLDecayFile(dec_name,verbose);
  56   }
  57
  58   return theDalitzTable;
  59
  60 }
  61
  62 bool EvtDalitzTable::fileHasBeenRead(const std::string dec_name) {
  63   std::vector<std::string>::iterator i = _readFiles.begin();
  64   for( ; i!=_readFiles.end(); i++) {
  65     if((*i).compare(dec_name) == 0) {
  66       return true;
  67     }
  68   }
  69   return false;
  70 }
  71
  72 void EvtDalitzTable::readXMLDecayFile(const std::string dec_name, bool verbose){
  73
  74   if (verbose) {
  75     report(INFO,"EvtGen")<<"EvtDalitzTable: Reading in xml parameter file "<<dec_name<<endl;
  76   }
  77
  78   _readFiles.push_back(dec_name);
  79
  80   EvtDalitzDecayInfo* dalitzDecay = 0;
  81   double probMax = 0;
  82   EvtId ipar;
  83   std::string decayParent = "";
  84   std::string daugStr = "";
  85   EvtId daughter[3];
  86
  87   EvtDalitzPlot dp;
  88   EvtComplex cAmp;
  89   std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> > angAndResPairs;
  90   std::string shape("");
  91   EvtSpinType::spintype spinType(EvtSpinType::SCALAR);
  92   double mass(0.), width(0.), FFp(0.), FFr(0.);
  93   std::vector<EvtFlatteParam> flatteParams;
  94   //Nonres parameters
  95   double alpha(0.);
  96   //LASS parameters
  97   double aLass(0.), rLass(0.), BLass(0.), phiBLass(0.), RLass(0.), phiRLass(0.), cutoffLass(-1.);
  98
  99   EvtParserXml parser;
 100   parser.open(dec_name);
 101
 102   bool endReached = false;
 103
 104   while(parser.readNextTag()) {
 105     //TAGS FOUND UNDER DATA
 106     if(parser.getParentTagTitle() == "data") {
 107       if(parser.getTagTitle() == "dalitzDecay") {
 108         int nDaughters = 0;
 109
 110         decayParent = parser.readAttribute("particle");
 111         daugStr = parser.readAttribute("daughters");
 112         probMax = parser.readAttributeDouble("probMax",-1);
 113
 114         checkParticle(decayParent);
 115         ipar=EvtPDL::getId(decayParent);
 116
 117         std::istringstream daugStream(daugStr);
 118
 119         std::string daugh;
 120         while(std::getline(daugStream, daugh, ' ')) {
 121           checkParticle(daugh);
 122           daughter[nDaughters++] = EvtPDL::getId(daugh);
 123         }
 124
 125         if(nDaughters!=3) {
 126           report(ERROR,"EvtGen") <<
 127                 "Expected to find three daughters for dalitzDecay of "<<decayParent<<" near line "
 128                 <<parser.getLineNumber()<<", "<<"found "<<nDaughters<<endl;
 129               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 130               ::abort();
 131         }
 132
 133         double m_d1 = EvtPDL::getMass(daughter[0]), m_d2 = EvtPDL::getMass(daughter[1]), m_d3 = EvtPDL::getMass(daughter[2]), M = EvtPDL::getMass(ipar);
 134
 135         dp = EvtDalitzPlot( m_d1, m_d2, m_d3, M );
 136
 137         dalitzDecay = new EvtDalitzDecayInfo(daughter[0],daughter[1],daughter[2]);
 138
 139       } else if(parser.getTagTitle() == "copyDalitz") {
 140         int nDaughters = 0;
 141         EvtId daughter[3];
 142         int nCopyDaughters = 0;
 143         EvtId copyDaughter[3];
 144
 145         decayParent = parser.readAttribute("particle");
 146         daugStr = parser.readAttribute("daughters");
 147
 148         std::string copyParent = parser.readAttribute("copy");
 149         std::string copyDaugStr = parser.readAttribute("copyDaughters");
 150
 151         checkParticle(decayParent);
 152         ipar=EvtPDL::getId(decayParent);
 153
 154         checkParticle(copyParent);
 155         EvtId copypar=EvtPDL::getId(copyParent);
 156
 157         std::istringstream daugStream(daugStr);
 158         std::istringstream copyDaugStream(copyDaugStr);
 159
 160         std::string daugh;
 161         while(std::getline(daugStream, daugh, ' ')) {
 162           checkParticle(daugh);
 163           daughter[nDaughters++] = EvtPDL::getId(daugh);
 164         }
 165         while(std::getline(copyDaugStream, daugh, ' ')) {
 166           checkParticle(daugh);
 167           copyDaughter[nCopyDaughters++] = EvtPDL::getId(daugh);
 168         }
 169
 170         if(nDaughters!=3 || nCopyDaughters!=3) {
 171           report(ERROR,"EvtGen") <<
 172                 "Expected to find three daughters for copyDecay of "<<decayParent<<
 173                 " from "<<copyParent<<" near line "<<parser.getLineNumber()<<
 174                 ", "<<"found "<<nDaughters<<" and "<<nCopyDaughters<<endl;
 175           report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 176               ::abort();
 177         }
 178
 179         copyDecay(ipar, daughter, copypar, copyDaughter);
 180
 181       } else if(parser.getTagTitle() == "/data") { //end of data
 182         endReached = true;
 183         parser.close();
 184         break;
 185       }
 186     //TAGS FOUND UNDER DALITZDECAY
 187     } else if(parser.getParentTagTitle() == "dalitzDecay") {
 188       if(parser.getTagTitle() == "resonance") {
 189
 190         flatteParams.clear();
 191
 192         //Amplitude
 193         EvtComplex ampFactor(parser.readAttributeDouble("ampFactorReal",1.),
 194                              parser.readAttributeDouble("ampFactorImag",0.));
 195         double mag = parser.readAttributeDouble("mag",-999.);
 196         double phase = parser.readAttributeDouble("phase",-999.);
 197         double real = parser.readAttributeDouble("real",-999.);
 198         double imag = parser.readAttributeDouble("imag",-999.);
 199
 200         if((real!=-999. || imag!=-999.) && mag==-999. && phase==-999.) {
 201           if(real==-999.) { real = 0; }
 202           if(imag==-999.) { imag = 0; }
 203           mag = sqrt(real*real + imag*imag);
 204           phase = atan2(imag,real) * EvtConst::radToDegrees;
 205         }
 206         if( mag==-999. ) {
 207           mag = 1.;
 208         }
 209         if( phase==-999. ) {
 210           phase = 0.;
 211         }
 212
 213         cAmp = ampFactor*EvtComplex(cos(phase*1.0/EvtConst::radToDegrees),sin(phase*1.0/EvtConst::radToDegrees))*mag;
 214
 215         //Resonance particle properties
 216         mass = 0.;
 217         width = 0.;
 218         spinType = EvtSpinType::SCALAR;
 219
 220         std::string particle = parser.readAttribute("particle");
 221         if(particle != "") {
 222           EvtId resId = EvtPDL::getId(particle);
 223           if(resId == EvtId(-1,-1)) {
 224             report(ERROR,"EvtGen") <<"Unknown particle name:"<<particle.c_str()<<endl;
 225             report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 226             ::abort();
 227           } else {
 228             mass = EvtPDL::getMeanMass(resId);
 229             width = EvtPDL::getWidth(resId);
 230             spinType = EvtPDL::getSpinType(resId);
 231           }
 232         }
 233
 234         width = parser.readAttributeDouble("width",width);
 235         mass = parser.readAttributeDouble("mass",mass);
 236         switch(parser.readAttributeInt("spin",-1)) {
 237         case -1://not set here
 238           break;
 239         case 0:
 240           spinType = EvtSpinType::SCALAR;
 241           break;
 242         case 1:
 243           spinType = EvtSpinType::VECTOR;
 244           break;
 245         case 2:
 246           spinType = EvtSpinType::TENSOR;
 247           break;
 248         default:
 249           report(ERROR,"EvtGen") << "Unsupported spin near line "<<parser.getLineNumber()<<" of XML file."<<endl;
 250           ::abort();
 251         }
 252
 253         //Shape and form factors
 254         shape = parser.readAttribute("shape");
 255         FFp = parser.readAttributeDouble("BlattWeisskopfFactorParent",0.0);
 256         FFr = parser.readAttributeDouble("BlattWeisskopfFactorResonance",1.5);
 257
 258         //Shape specific attributes
 259         if(shape=="NonRes_Exp") {
 260           alpha = parser.readAttributeDouble("alpha",0.0);
 261         }
 262         if(shape=="LASS") {
 263           aLass = parser.readAttributeDouble("a",2.07);
 264           rLass = parser.readAttributeDouble("r",3.32);
 265           BLass = parser.readAttributeDouble("B",1.0);
 266           phiBLass = parser.readAttributeDouble("phiB",0.0);
 267           RLass = parser.readAttributeDouble("R",1.0);
 268           phiRLass = parser.readAttributeDouble("phiR",0.0);
 269           cutoffLass = parser.readAttributeDouble("cutoff",-1.0);
 270         }
 271
 272         //Daughter pairs for resonance
 273         angAndResPairs.clear();
 274
 275         std::string resDaugStr = parser.readAttribute("resDaughters");
 276         if(resDaugStr != "") {
 277           std::istringstream resDaugStream(resDaugStr);
 278           std::string resDaug;
 279           int nResDaug(0);
 280           EvtId resDaughter[2];
 281           while(std::getline(resDaugStream, resDaug, ' ')) {
 282             checkParticle(resDaug);
 283             resDaughter[nResDaug++] = EvtPDL::getId(resDaug);
 284           }
 285           if(nResDaug != 2) {
 286             report(ERROR,"EvtGen") << "Resonance must have exactly 2 daughters near line "<<parser.getLineNumber()<<" of XML file."<<endl;
 287             ::abort();
 288           }
 289           int nRes = getDaughterPairs(resDaughter, daughter, angAndResPairs);
 290           if(nRes==0) {
 291             report(ERROR,"EvtGen") << "Resonance daughters must match decay daughters near line "<<parser.getLineNumber()<<" of XML file."<<endl;
 292             ::abort();
 293           }
 294           if(parser.readAttributeBool("normalise",true)) cAmp /= sqrt(nRes);
 295         }
 296         if(angAndResPairs.empty()) {
 297           switch(parser.readAttributeInt("daughterPair")) {
 298           case 1:
 299             angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::AB));
 300             break;
 301           case 2:
 302             angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::BC));
 303             break;
 304           case 3:
 305             angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::CA));
 306             break;
 307           default:
 308             if(shape == "NonRes") { //We don't expect a pair for non-resonant terms but add dummy values for convenience
 309               angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::AB));
 310               break;
 311             }
 312             report(ERROR,"EvtGen") << "Daughter pair must be 1, 2 or 3 near line "<<parser.getLineNumber()<<" of XML file."<<endl;
 313             ::abort();
 314           }
 315         }
 316
 317         if(parser.isTagInline()) {
 318           std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >::iterator it = angAndResPairs.begin();
 319           for( ; it != angAndResPairs.end(); it++) {
 320             std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> pairs = *it;
 321             EvtDalitzReso resonance = getResonance(shape, dp, pairs.first, pairs.second, spinType, mass, width, FFp, FFr, alpha, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass);
 322             dalitzDecay->addResonance(cAmp,resonance);
 323           }
 324         }
 325       } else if(parser.getTagTitle() == "/dalitzDecay") {
 326         if(probMax < 0) {
 327           report(INFO,"EvtGen") << "probMax is not defined for " << decayParent << " -> " << daugStr << endl;
 328           report(INFO,"EvtGen") << "Will now estimate probMax. This may take a while. Once probMax is calculated, update the XML file to skip this step in future." << endl;
 329           probMax = calcProbMax(dp,dalitzDecay);
 330         }
 331         dalitzDecay->setProbMax(probMax);
 332         addDecay(ipar, *dalitzDecay);
 333         delete dalitzDecay;
 334         dalitzDecay = 0;
 335       } else if(verbose) {
 336         report(INFO,"EvtGen") << "Unexpected tag "<<parser.getTagTitle()
 337                               <<" found in XML decay file near line "
 338                               <<parser.getLineNumber()<<". Tag will be ignored."<<endl;
 339       }
 340     //TAGS FOUND UNDER RESONANCE
 341     } else if(parser.getParentTagTitle() == "resonance"){
 342       if(parser.getTagTitle() == "flatteParam") {
 343         EvtFlatteParam param(parser.readAttributeDouble("mass1"),
 344                              parser.readAttributeDouble("mass2"),
 345                              parser.readAttributeDouble("g"));
 346         flatteParams.push_back(param);
 347       } else if(parser.getTagTitle() == "/resonance") {
 348         std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >::iterator it = angAndResPairs.begin();
 349         for( ; it != angAndResPairs.end(); it++) {
 350           std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> pairs = *it;
 351           EvtDalitzReso resonance = getResonance(shape, dp, pairs.first, pairs.second, spinType, mass, width, FFp, FFr, alpha, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass);
 352
 353           std::vector<EvtFlatteParam>::iterator flatteIt = flatteParams.begin();
 354           for( ; flatteIt != flatteParams.end(); flatteIt++) {
 355             resonance.addFlatteParam((*flatteIt));
 356           }
 357
 358           dalitzDecay->addResonance(cAmp,resonance);
 359         }
 360       }
 361     }
 362   }
 363
 364   if(!endReached) {
 365     report(ERROR,"EvtGen") << "Either the decay file ended prematurely or the file is badly formed.\n"
 366                           <<"Error occured near line"<<parser.getLineNumber()<<endl;
 367     ::abort();
 368   }
 369 }
 370
 371 void EvtDalitzTable::checkParticle(std::string particle) {
 372   if (EvtPDL::getId(particle)==EvtId(-1,-1)) {
 373     report(ERROR,"EvtGen") <<"Unknown particle name:"<<particle.c_str()<<endl;
 374     report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 375     ::abort();
 376   }
 377 }
 378
 379 void EvtDalitzTable::addDecay(EvtId parent, const EvtDalitzDecayInfo& dec) {
 380   if(_dalitztable.find(parent)!=_dalitztable.end()) {
 381     _dalitztable[parent].push_back(dec);
 382   } else {
 383     _dalitztable[parent].push_back(dec);
 384   }
 385 }
 386
 387 void EvtDalitzTable::copyDecay(EvtId parent, EvtId* daughters,
 388                                EvtId copy, EvtId* copyd) {
 389   EvtDalitzDecayInfo decay(daughters[0],daughters[1],daughters[2]);
 390   std::vector<EvtDalitzDecayInfo> copyTable = getDalitzTable(copy);
 391   std::vector<EvtDalitzDecayInfo>::iterator i = copyTable.begin();
 392   for( ; i != copyTable.end(); i++) {
 393     EvtId daughter1 = (*i).daughter1();
 394     EvtId daughter2 = (*i).daughter2();
 395     EvtId daughter3 = (*i).daughter3();
 396
 397     if((copyd[0] == daughter1 && copyd[1] == daughter2 && copyd[2] == daughter3) ||
 398        (copyd[0] == daughter1 && copyd[1] == daughter3 && copyd[2] == daughter2) ||
 399        (copyd[0] == daughter2 && copyd[1] == daughter1 && copyd[2] == daughter3) ||
 400        (copyd[0] == daughter2 && copyd[1] == daughter3 && copyd[2] == daughter1) ||
 401        (copyd[0] == daughter3 && copyd[1] == daughter1 && copyd[2] == daughter2) ||
 402        (copyd[0] == daughter3 && copyd[1] == daughter2 && copyd[2] == daughter1)) {
 403       decay.setProbMax((*i).getProbMax());
 404       std::vector<std::pair<EvtComplex, EvtDalitzReso> >::const_iterator j = (*i).getResonances().begin();
 405       for( ; j != (*i).getResonances().end(); j++) {
 406         decay.addResonance((*j));
 407       }
 408       addDecay(parent,decay);
 409       return;
 410     }
 411   }
 412   //if we get here then there was no match
 413   report(ERROR,"EvtGen") << "Did not find dalitz decays for particle:"
 414          <<copy<<"\n";
 415 }
 416
 417 std::vector<EvtDalitzDecayInfo> EvtDalitzTable::getDalitzTable(const EvtId& parent) {
 418   std::vector<EvtDalitzDecayInfo> table;
 419   if ( _dalitztable.find(parent)!=_dalitztable.end() ) {
 420     table=_dalitztable[parent];
 421   }
 422
 423   if (table.empty()){
 424     report(ERROR,"EvtGen") << "Did not find dalitz decays for particle:"
 425          <<parent<<"\n";
 426   }
 427
 428   return table;
 429 }
 430
 431
 432 EvtDalitzReso EvtDalitzTable::getResonance(std::string shape, EvtDalitzPlot dp, EvtCyclic3::Pair angPair, EvtCyclic3::Pair resPair,
 433                                            EvtSpinType::spintype spinType, double mass, double width, double FFp, double FFr, double alpha,
 434                                            double aLass, double rLass, double BLass, double phiBLass, double RLass, double phiRLass, double cutoffLass) {
 435   if( shape=="RBW" || shape=="RBW_CLEO") {
 436     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::RBW_CLEO, FFp, FFr );
 437   } else if( shape=="RBW_CLEO_ZEMACH" ) {
 438     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::RBW_CLEO_ZEMACH, FFp, FFr );
 439   } else if( shape=="GS" || shape=="GS_CLEO" ) {
 440     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GS_CLEO, FFp, FFr );
 441   } else if( shape=="GS_CLEO_ZEMACH" ) {
 442     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GS_CLEO_ZEMACH, FFp, FFr );
 443   } else if( shape=="GAUSS" || shape=="GAUSS_CLEO" ) {
 444     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GAUSS_CLEO, FFp, FFr );
 445   } else if( shape=="GAUSS_CLEO_ZEMACH" ) {
 446     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GAUSS_CLEO_ZEMACH, FFp, FFr );
 447   } else if( shape=="Flatte" ) {
 448     return EvtDalitzReso( dp, resPair, mass );
 449   } else if( shape=="LASS" ) {
 450     return EvtDalitzReso( dp, resPair, mass, width, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass, true );
 451   } else if( shape=="NonRes" ) {
 452     return EvtDalitzReso( );
 453   } else if( shape=="NonRes_Linear" ) {
 454     return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_LIN );
 455   } else if( shape=="NonRes_Exp" ) {
 456     return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_EXP, alpha );
 457   } else { //NBW
 458     if( shape!="NBW") report(WARNING,"EvtGen")<<"EvtDalitzTable: shape "<<shape<<" is unknown. Defaulting to NBW."<<endl;
 459     return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::NBW, FFp, FFr );
 460   }
 461 }
 462
 463 int EvtDalitzTable::getDaughterPairs(EvtId* resDaughter, EvtId* daughter, std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >& angAndResPairs) {
 464   int n(0);
 465   if(resDaughter[0]==daughter[0] && resDaughter[1]==daughter[1]) {
 466     angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::AB)); n++;
 467   } else if(resDaughter[0]==daughter[1] && resDaughter[1]==daughter[0]) {
 468     angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::AB)); n++;
 469   }
 470
 471   if(resDaughter[0]==daughter[1] && resDaughter[1]==daughter[2]) {
 472     angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::BC)); n++;
 473   } else if(resDaughter[0]==daughter[2] && resDaughter[1]==daughter[1]) {
 474     angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::BC)); n++;
 475   }
 476
 477   if(resDaughter[0]==daughter[2] && resDaughter[1]==daughter[0]) {
 478     angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::CA)); n++;
 479   } else if(resDaughter[0]==daughter[0] && resDaughter[1]==daughter[2]) {
 480     angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::CA)); n++;
 481   }
 482
 483   return n;
 484 }
 485
 486 double EvtDalitzTable::calcProbMax(EvtDalitzPlot dp, EvtDalitzDecayInfo* model) {
 487
 488   double factor = 1.2; //factor to increase our final answer by
 489   int nStep(1000);      //number of steps - total points will be 3*nStep*nStep
 490
 491   double maxProb(0);
 492   double min(0), max(0), step(0), min2(0), max2(0), step2(0);
 493
 494   //first do AB, BC
 495   min = dp.qAbsMin(EvtCyclic3::AB);
 496   max = dp.qAbsMax(EvtCyclic3::AB);
 497   step = (max-min)/nStep;
 498   for(int i=0; i<nStep; ++i) {
 499     double qAB = min + i*step;
 500     min2 = dp.qMin(EvtCyclic3::BC,EvtCyclic3::AB,qAB);
 501     max2 = dp.qMax(EvtCyclic3::BC,EvtCyclic3::AB,qAB);
 502     step2 = (max2-min2)/nStep;
 503     for(int j=0; j<nStep; ++j) {
 504       double qBC = min2+ j*step2;
 505       EvtDalitzCoord coord(EvtCyclic3::AB,qAB,EvtCyclic3::BC,qBC);
 506       EvtDalitzPoint point(dp,coord);
 507       double prob = calcProb(point,model);
 508       if(prob > maxProb) maxProb = prob;
 509     }
 510   }
 511
 512   //next do BC, CA
 513   min = dp.qAbsMin(EvtCyclic3::BC);
 514   max = dp.qAbsMax(EvtCyclic3::BC);
 515   step = (max-min)/nStep;
 516   for(int i=0; i<nStep; ++i) {
 517     double qBC = min + i*step;
 518     min2 = dp.qMin(EvtCyclic3::CA,EvtCyclic3::BC,qBC);
 519     max2 = dp.qMax(EvtCyclic3::CA,EvtCyclic3::BC,qBC);
 520     step2 = (max2-min2)/nStep;
 521     for(int j=0; j<nStep; ++j) {
 522       double qCA = min2+ j*step2;
 523       EvtDalitzCoord coord(EvtCyclic3::BC,qBC,EvtCyclic3::CA,qCA);
 524       EvtDalitzPoint point(dp,coord);
 525       double prob = calcProb(point,model);
 526       if(prob > maxProb) maxProb = prob;
 527     }
 528   }
 529
 530   //finally do CA, AB
 531   min = dp.qAbsMin(EvtCyclic3::CA);
 532   max = dp.qAbsMax(EvtCyclic3::CA);
 533   step = (max-min)/nStep;
 534   for(int i=0; i<nStep; ++i) {
 535     double qCA = min + i*step;
 536     min2 = dp.qMin(EvtCyclic3::AB,EvtCyclic3::CA,qCA);
 537     max2 = dp.qMax(EvtCyclic3::AB,EvtCyclic3::CA,qCA);
 538     step2 = (max2-min2)/nStep;
 539     for(int j=0; j<nStep; ++j) {
 540       double qAB = min2+ j*step2;
 541       EvtDalitzCoord coord(EvtCyclic3::CA,qCA,EvtCyclic3::AB,qAB);
 542       EvtDalitzPoint point(dp,coord);
 543       double prob = calcProb(point,model);
 544       if(prob > maxProb) maxProb = prob;
 545     }
 546   }
 547   report(INFO,"EvtGen") << "Largest probability found was " << maxProb << endl;
 548   report(INFO,"EvtGen") << "Setting probMax to " << factor*maxProb << endl;
 549   return factor*maxProb;
 550 }
 551
 552 double EvtDalitzTable::calcProb(EvtDalitzPoint point, EvtDalitzDecayInfo* model) {
 553
 554   std::vector<std::pair<EvtComplex,EvtDalitzReso> > resonances = model->getResonances();
 555
 556   EvtComplex amp(0,0);
 557   std::vector<std::pair<EvtComplex,EvtDalitzReso> >::iterator i = resonances.begin();
 558   for( ; i!= resonances.end(); i++) {
 559     std::pair<EvtComplex,EvtDalitzReso> res = (*i);
 560     amp += res.first * res.second.evaluate( point );
 561   }
 562   return abs2(amp);
 563 }