TEvtGen/EvtGenBase/EvtMultiChannelParser.cxx

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtMultiChannelParser.cc,v 1.17 2009/02/19 03:22:30 ryd 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 #include <assert.h>
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <string>
  23 #include <vector>
  24 #include "EvtGenBase/EvtParser.hh"
  25 #include "EvtGenBase/EvtMultiChannelParser.hh"
  26 #include "EvtGenBase/EvtDecayMode.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28
  29 using std::string;
  30 using std::vector;
  31
  32 EvtDecayMode EvtMultiChannelParser::getDecayMode(const char* file)
  33 {
  34   // Open file, read tokens
  35
  36   EvtParser parser;
  37   parser.read(file);
  38
  39   // Seek Decay
  40
  41   int i = 0;
  42   int N = parser.getNToken();
  43   while(i<N) {
  44
  45     std::string tok = parser.getToken(i++);
  46     if(tok == std::string("Decay")) break;
  47   }
  48
  49   // Get mother
  50
  51   string mother = string(parser.getToken(i++).c_str());
  52   std::string bf = parser.getToken(i++);
  53
  54   vector<string> dauV;
  55   // Get daughters
  56
  57   while(1) {
  58
  59     std::string d = parser.getToken(i++);
  60
  61     if(EvtPDL::getStdHep(EvtPDL::getId(d.c_str())) == 0) break;
  62
  63     dauV.push_back(string(d.c_str()));
  64   }
  65
  66   EvtDecayMode mode(mother,dauV);
  67   printf("Decay File defines mode %s\n",mode.mode().c_str());
  68
  69   return mode;
  70 }
  71
  72
  73
  74 void EvtMultiChannelParser::parse(const char* file, const char* model)
  75 {
  76   // Open file, read tokens
  77
  78   EvtParser parser;
  79   parser.read(file);
  80
  81
  82   // Get parameters (tokens between the model name and ;)
  83
  84   int i = 0;
  85   int N = parser.getNToken();
  86
  87   // Seek the model name
  88
  89   while(i<N) {
  90
  91     std::string tok = parser.getToken(i++);
  92     if(tok == std::string(model)) break;
  93   }
  94   if(i == N) {
  95
  96     printf("No model %s found in decay file %s",model,file);
  97     exit(0);
  98   }
  99
 100
 101   // Add all tokens up to a semicolon to vector
 102
 103   std::vector<std::string> v;
 104   while(i<N) {
 105
 106     std::string tok = parser.getToken(i++);
 107     if(tok == std::string(";")) break;
 108     else v.push_back(tok);
 109   }
 110
 111   if(i == N) {
 112
 113     printf("No terminating ; found in decay file %s",file);
 114     assert(0);
 115   }
 116
 117   parse(v);
 118 }
 119
 120
 121 void EvtMultiChannelParser::parse(const std::vector<std::string>& v)
 122 {
 123   // place holder for strtod
 124   char** tc = 0;
 125
 126
 127   // Get PDF maximum or number of points to
 128   // use in the scan.
 129
 130   if(v[0] == std::string("MAXPDF")) {
 131
 132     _pdfMax = strtod(v[1].c_str(),tc);
 133     if(_pdfMax <= 0) { printf("Bad pdfMax=%f\n",_pdfMax); assert(0); }
 134   }
 135   else
 136     if(v[0] == std::string("SCANPDF")) {
 137
 138       _nScan = atoi(v[1].c_str());
 139     }
 140     else {
 141
 142       printf("Error parsing decay file\n");
 143       assert(0);
 144     }
 145
 146
 147   // Now parse the rest of file for amplitude specifications.
 148
 149   bool conjugate = false;
 150   size_t i = 2;
 151   assert(isKeyword(v[2]));
 152
 153   while(i  < v.size()) {
 154
 155     size_t i0 = i;
 156
 157     // Switch to conjugate amplitudes after keyword
 158     if(v[i] == std::string("CONJUGATE")) {
 159
 160       assert(conjugate == false);
 161       conjugate = true;
 162       assert(!isKeyword(v[++i]));
 163       _dm =  strtod(v[i++].c_str(),tc);
 164       _mixAmpli = strtod(v[i++].c_str(),tc);
 165       _mixPhase = strtod(v[i++].c_str(),tc);
 166     }
 167
 168     if (i >= v.size()) break;
 169     std::vector<std::string> params;
 170     EvtComplex c;
 171     int format;
 172
 173     if(!conjugate && v[i] == std::string("AMPLITUDE")) {
 174
 175       while(!isKeyword(v[++i])) params.push_back(v[i]);
 176       _amp.push_back(params);
 177
 178       parseComplexCoef(i,v,c,format);
 179       _ampCoef.push_back(c);
 180       _coefFormat.push_back(format);
 181       continue;
 182     }
 183     else
 184       if(conjugate && v[i] == std::string("AMPLITUDE")) {
 185
 186         while(!isKeyword(v[++i])) params.push_back(v[i]);
 187         _ampConj.push_back(params);
 188         parseComplexCoef(i,v,c,format);
 189         _ampConjCoef.push_back(c);
 190         _coefConjFormat.push_back(format);
 191         continue;
 192       }
 193       else {
 194
 195         printf("Expect keyword, found parameter %s\n",v[i].c_str());
 196         assert(0);
 197       }
 198
 199
 200     assert(i > i0);
 201   }
 202
 203   printf("PARSING SUCCESSFUL\n");
 204   printf("%d amplitude terms\n",(int)_amp.size());
 205   printf("%d conj amplitude terms\n",(int)_ampConj.size());
 206 }
 207
 208
 209
 210 void EvtMultiChannelParser::parseComplexCoef(size_t& i, const std::vector<std::string>& v,
 211                                              EvtComplex& c, int& format)
 212 {
 213   // place holder for strtod
 214   char** tc = 0;
 215
 216   assert(v[i++] == std::string("COEFFICIENT"));
 217
 218   if(v[i] == std::string("POLAR_DEG")) {
 219
 220       double mag = strtod(v[i+1].c_str(),tc);
 221       double phaseRad = strtod(v[i+2].c_str(),tc)*EvtConst::pi/180.0;
 222       i += 3;
 223       c = EvtComplex(mag*cos(phaseRad),mag*sin(phaseRad));
 224       format = POLAR_DEG;
 225   }
 226   else if(v[i] == std::string("POLAR_RAD")) {
 227
 228     double mag = strtod(v[i+1].c_str(),tc);
 229     double phaseRad = strtod(v[i+2].c_str(),tc);
 230     i += 3;
 231     c = EvtComplex(mag*cos(phaseRad),mag*sin(phaseRad));
 232     format = POLAR_RAD;
 233   }
 234   else if(v[i] == std::string("CARTESIAN")) {
 235
 236     double re = strtod(v[i+1].c_str(),tc);
 237     double im = strtod(v[i+2].c_str(),tc);
 238     i += 3;
 239     c = EvtComplex(re,im);
 240     format = CARTESIAN;
 241   }
 242   else {
 243
 244     printf("Invalid format %s for complex coefficient\n",v[i].c_str());
 245     exit(0);
 246   }
 247 }
 248
 249
 250 double EvtMultiChannelParser::parseRealCoef(int& i, const std::vector<std::string>& v)
 251 {
 252   // place holder for strtod
 253   char** tc = 0;
 254   double value = 0;
 255
 256   if(v[i] == std::string("COEFFICIENT")) {
 257
 258     value = strtod(v[i+1].c_str(),tc);
 259   }
 260   else assert(0);
 261
 262   i += 2;
 263
 264   assert(value > 0.);
 265   return value;
 266 }
 267
 268
 269 bool EvtMultiChannelParser::isKeyword(const std::string& s)
 270 {
 271   if(s == std::string("AMPLITUDE")) return true;
 272   if(s == std::string("CONJUGATE")) return true;
 273   if(s == std::string("COEFFICIENT")) return true;
 274   return false;
 275 }
 276
 277
 278