TEvtGen/EvtGen/EvtGenBase/EvtPto3PAmpFactory.cpp

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtPto3PAmpFactory.cpp,v 1.3 2009-03-16 15:44:04 robbep Exp $
   4 //
   5 // Environment:
   6 //      This software is part of the EvtGen package developed jointly
   7 //      for the BaBar and CLEO collaborations. If you use all or part
   8 //      of it, please give an appropriate acknowledgement.
   9 //
  10 // Copyright Information:
  11 //      Copyright (C) 1998 Caltech, UCSB
  12 //
  13 // Module creator:
  14 //      Alexei Dvoretskii, Caltech, 2001-2002.
  15 //-----------------------------------------------------------------------
  16 #include "EvtGenBase/EvtPatches.hh"
  17
  18 // AmpFactory for building a P -> 3P decay
  19 // (pseudoscalar to three pseudoscalars)
  20
  21 #include <assert.h>
  22 #include <math.h>
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25
  26 #include "EvtGenBase/EvtId.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtConst.hh"
  29 #include "EvtGenBase/EvtComplex.hh"
  30 #include "EvtGenBase/EvtCyclic3.hh"
  31 #include "EvtGenBase/EvtSpinType.hh"
  32 #include "EvtGenBase/EvtPto3PAmp.hh"
  33 #include "EvtGenBase/EvtNonresonantAmp.hh"
  34 #include "EvtGenBase/EvtFlatAmp.hh"
  35 #include "EvtGenBase/EvtLASSAmp.hh"
  36 #include "EvtGenBase/EvtPto3PAmpFactory.hh"
  37 #include "EvtGenBase/EvtPropBreitWigner.hh"
  38 #include "EvtGenBase/EvtPropFlatte.hh"
  39 #include "EvtGenBase/EvtPropBreitWignerRel.hh"
  40 #include "EvtGenBase/EvtDalitzResPdf.hh"
  41 #include "EvtGenBase/EvtDalitzFlatPdf.hh"
  42
  43 using namespace EvtCyclic3;
  44 #include <iostream>
  45
  46 void EvtPto3PAmpFactory::processAmp(EvtComplex c, std::vector<std::string> vv, bool conj)
  47 {
  48   if(_verbose) {
  49
  50     printf("Make %samplitude\n",conj ? "CP conjugate" : "");
  51     unsigned i;
  52     for(i=0;i<vv.size();i++) printf("%s\n",vv[i].c_str());
  53     printf("\n");
  54   }
  55
  56   EvtAmplitude<EvtDalitzPoint>* amp = 0;
  57   EvtPdf<EvtDalitzPoint>* pdf = 0;
  58   std::string name;
  59   Pair pairRes=AB;
  60
  61   size_t i;
  62   /*
  63          Experimental amplitudes
  64   */
  65   if(vv[0] == "PHASESPACE") {
  66
  67     pdf = new EvtDalitzFlatPdf(_dp);
  68     amp = new EvtFlatAmp<EvtDalitzPoint>();
  69     name = "NR";
  70   }
  71   else if (!vv[0].find("NONRES")) {
  72     double alpha=0;
  73     EvtPto3PAmp::NumType typeNRes=EvtPto3PAmp::NONRES;
  74     if      (vv[0]=="NONRES_LIN") {
  75       typeNRes=EvtPto3PAmp::NONRES_LIN;
  76       pairRes=strToPair(vv[1].c_str());
  77     }
  78     else if (vv[0]=="NONRES_EXP") {
  79       typeNRes=EvtPto3PAmp::NONRES_EXP;
  80       pairRes = strToPair(vv[1].c_str());
  81       alpha   = strtod(vv[2].c_str(),0);
  82     }
  83     else assert(0);
  84     pdf = new EvtDalitzFlatPdf(_dp);
  85     amp = new EvtNonresonantAmp( &_dp, typeNRes, pairRes, alpha);
  86   }
  87   else if (vv[0]=="LASS" || vv[0]=="LASS_ELASTIC" || vv[0]=="LASS_RESONANT") {
  88     pairRes = strToPair(vv[1].c_str());
  89     double m0 = strtod(vv[2].c_str(),0);
  90     double g0 = strtod(vv[3].c_str(),0);
  91     double a  = strtod(vv[4].c_str(),0);
  92     double r  = strtod(vv[5].c_str(),0);
  93     double cutoff  = strtod(vv[6].c_str(),0);
  94     pdf = new EvtDalitzResPdf(_dp,m0,g0,pairRes);
  95     amp = new EvtLASSAmp( &_dp, pairRes, m0, g0, a, r, cutoff, vv[0]);
  96   }
  97
  98   /*
  99       Resonant amplitudes
 100   */
 101   else if(vv[0] == "RESONANCE") {
 102     EvtPto3PAmp* partAmp = 0;
 103
 104     // RESONANCE stanza
 105
 106     pairRes = strToPair(vv[1].c_str());
 107     EvtSpinType::spintype spinR = EvtSpinType::SCALAR;
 108     double mR, gR;
 109     name = vv[2];
 110     EvtId resId = EvtPDL::getId(vv[2]);
 111     if(_verbose) printf("Particles %s form %sresonance %s\n",
 112                         vv[1].c_str(),vv[2].c_str(), conj ? "(conj) " : "");
 113
 114     // If no valid particle name is given, assume that
 115     // it is the spin, the mass and the width of the particle.
 116
 117     if(resId.getId() == -1) {
 118
 119       switch(atoi(vv[2].c_str())) {
 120
 121       case 0: { spinR = EvtSpinType::SCALAR; break; }
 122       case 1: { spinR = EvtSpinType::VECTOR; break; }
 123       case 2: { spinR = EvtSpinType::TENSOR; break; }
 124       case 3: { spinR = EvtSpinType::SPIN3; break; }
 125       case 4: { spinR = EvtSpinType::SPIN4; break; }
 126       default: { assert(0); break; }
 127       }
 128
 129       mR = strtod(vv[3].c_str(),0);
 130       gR = strtod(vv[4].c_str(),0);
 131       i = 4;
 132     }
 133     else {
 134
 135       // For a valid particle get spin, mass and width
 136
 137       spinR = EvtPDL::getSpinType(resId);
 138       mR = EvtPDL::getMeanMass(resId);
 139       gR = EvtPDL::getWidth(resId);
 140       i = 2;
 141
 142       // It's possible to specify mass and width of a particle
 143       // explicitly
 144
 145       if(vv[3] != "ANGULAR") {
 146
 147         if(_verbose)
 148           printf("Setting m(%s)=%s g(%s)=%s\n",
 149                  vv[2].c_str(),vv[3].c_str(),vv[2].c_str(),vv[4].c_str());
 150
 151         mR = strtod(vv[3].c_str(),0);
 152         gR = strtod(vv[4].c_str(),0);
 153         i = 4;
 154       }
 155     }
 156
 157     // ANGULAR stanza
 158
 159     if(vv[++i] != "ANGULAR") {
 160
 161       printf("%s instead of ANGULAR\n",vv[i].c_str());
 162       exit(0);
 163     }
 164     Pair pairAng = strToPair(vv[++i].c_str());
 165     if(_verbose) printf("Angle is measured between particles %s\n",vv[i].c_str());
 166
 167     // TYPE stanza
 168
 169     std::string typeName = vv[++i];
 170     assert(typeName == "TYPE");
 171     std::string type = vv[++i];
 172     if(_verbose) printf("Propagator type %s\n",vv[i].c_str());
 173
 174     if(type == "NBW") {
 175
 176       EvtPropBreitWigner prop(mR,gR);
 177       partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::NBW);
 178     }
 179     else if(type == "RBW_ZEMACH") {
 180
 181       EvtPropBreitWignerRel prop(mR,gR);
 182       partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_ZEMACH);
 183     }
 184     else if(type == "RBW_KUEHN") {
 185
 186       EvtPropBreitWignerRel prop(mR,gR);
 187       partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_KUEHN);
 188     }
 189     else if(type == "RBW_CLEO") {
 190
 191       EvtPropBreitWignerRel prop(mR,gR);
 192       partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_CLEO);
 193     }
 194     else if(type == "FLATTE") {
 195
 196       double m1a = _dp.m( first(pairRes) );
 197       double m1b = _dp.m( second(pairRes) );
 198       // 2nd channel
 199       double g2  = strtod(vv[++i].c_str(),0);
 200       double m2a = strtod(vv[++i].c_str(),0);
 201       double m2b = strtod(vv[++i].c_str(),0);
 202       EvtPropFlatte  prop( mR, gR, m1a, m1b, g2, m2a, m2b );
 203       partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::FLATTE);
 204     }
 205     else assert(0);
 206
 207     // Optional DVFF, BVFF stanzas
 208
 209     if(i < vv.size() - 1) {
 210       if(vv[i+1] == "DVFF") {
 211         i++;
 212         if(vv[++i] == "BLATTWEISSKOPF") {
 213
 214           double R = strtod(vv[++i].c_str(),0);
 215           partAmp->set_fd(R);
 216         }
 217         else assert(0);
 218       }
 219     }
 220
 221     if(i < vv.size() - 1) {
 222       if(vv[i+1] == "BVFF") {
 223         i++;
 224         if(vv[++i] == "BLATTWEISSKOPF") {
 225
 226           if(_verbose) printf("BVFF=%s\n",vv[i].c_str());
 227           double R = strtod(vv[++i].c_str(),0);
 228           partAmp->set_fb(R);
 229         }
 230         else assert(0);
 231       }
 232     }
 233
 234     const int minwidths=5;
 235     //Optional resonance minimum and maximum
 236     if(i < vv.size() - 1) {
 237       if(vv[i+1] == "CUTOFF") {
 238         i++;
 239         if(vv[i+1] == "MIN") {
 240           i++;
 241           double min = strtod(vv[++i].c_str(),0);
 242           if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<" "<<minwidths<<std::endl;
 243           //ensure against cutting off too close to the resonance
 244           assert( min<(mR-minwidths*gR) );
 245           partAmp->setmin(min);
 246         }
 247         else if (vv[i+1] == "MAX") {
 248           i++;
 249           double max = strtod(vv[++i].c_str(),0);
 250           if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<" "<<minwidths<<std::endl;
 251           //ensure against cutting off too close to the resonance
 252           assert( max>(mR+minwidths*gR) );
 253           partAmp->setmax(max);
 254         }
 255         else assert(0);
 256       }
 257     }
 258
 259     //2nd iteration in case min and max are both specified
 260     if(i < vv.size() - 1) {
 261       if(vv[i+1] == "CUTOFF") {
 262         i++;
 263         if(vv[i+1] == "MIN") {
 264           i++;
 265           double min = strtod(vv[++i].c_str(),0);
 266           if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<std::endl;
 267           //ensure against cutting off too close to the resonance
 268           assert( min<(mR-minwidths*gR) );
 269           partAmp->setmin(min);
 270         }
 271         else if (vv[i+1] == "MAX") {
 272           i++;
 273           double max = strtod(vv[++i].c_str(),0);
 274           if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<std::endl;
 275           //ensure against cutting off too close to the resonance
 276           assert( max>(mR+minwidths*gR) );
 277           partAmp->setmax(max);
 278         }
 279         else assert(0);
 280       }
 281     }
 282
 283
 284     i++;
 285
 286     pdf = new EvtDalitzResPdf(_dp,mR,gR,pairRes);
 287     amp = partAmp;
 288   }
 289
 290   assert(amp);
 291   assert(pdf);
 292
 293   if(!conj) {
 294     _amp->addOwnedTerm(c,amp);
 295   }
 296   else {
 297     _ampConj->addOwnedTerm(c,amp);
 298   }
 299
 300   double scale = matchIsobarCoef(_amp, pdf, pairRes);
 301   _pc->addOwnedTerm(abs2(c)*scale,pdf);
 302
 303   _names.push_back(name);
 304 }
 305
 306 double EvtPto3PAmpFactory::matchIsobarCoef(EvtAmplitude<EvtDalitzPoint>* amp,
 307                                            EvtPdf<EvtDalitzPoint>* pdf,
 308                                            EvtCyclic3::Pair ipair) {
 309
 310   // account for differences in the definition of amplitudes by matching
 311   //        Integral( c'*pdf ) = Integral( c*|A|^2 )
 312   // to improve generation efficiency ...
 313
 314   double Ipdf  = pdf->compute_integral(10000).value();
 315   double Iamp2 = 0;
 316
 317
 318   EvtCyclic3::Pair jpair = EvtCyclic3::next(ipair);
 319   EvtCyclic3::Pair kpair = EvtCyclic3::next(jpair);
 320
 321   // Trapezoidal integral
 322   int N=10000;
 323
 324   double di = (_dp.qAbsMax(ipair) - _dp.qAbsMin(ipair))/((double) N);
 325
 326   double siMin = _dp.qAbsMin(ipair);
 327
 328   double s[3]; // playing with fire
 329   for(int i=1; i<N; i++) {
 330
 331     s[ipair] = siMin + di*i;
 332     s[jpair] = _dp.q(jpair, 0.9999, ipair, s[ipair]);
 333     s[kpair] = _dp.bigM()*_dp.bigM() - s[ipair] - s[jpair]
 334       + _dp.mA()*_dp.mA() + _dp.mB()*_dp.mB() + _dp.mC()*_dp.mC();
 335
 336     EvtDalitzPoint point( _dp.mA(), _dp.mB(), _dp.mC(),
 337                           s[EvtCyclic3::AB], s[EvtCyclic3::BC], s[EvtCyclic3::CA]);
 338
 339     if (!point.isValid()) continue;
 340
 341     double p = point.p(other(ipair), ipair);
 342     double q = point.p(first(ipair), ipair);
 343
 344     double itg = abs2( amp->evaluate(point) )*di*4*q*p;
 345     Iamp2 += itg;
 346
 347   }
 348   if (_verbose) std::cout << "integral = " << Iamp2 << "  pdf="<<Ipdf << std::endl;
 349
 350   assert(Ipdf>0 && Iamp2>0);
 351
 352   return Iamp2/Ipdf;
 353 }