TTherminator/Therminator/Parser.cxx

   1 /******************************************************************************
   2  *                      T H E R M I N A T O R                                 *
   3  *                   THERMal heavy-IoN generATOR                              *
   4  *                           version 1.0                                      *
   5  *                                                                            *
   6  * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
   7  *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
   8  * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
   9  *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
  10  * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
  11  *                 Wojciech Florkowski                                        *
  12  *                                                                            *
  13  * For the detailed description of the program and furhter references         *
  14  * to the description of the model plesase refer to: nucl-th/0504047,         *
  15  * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
  16  *                                                                            *
  17  * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
  18  *                                                                            *
  19  * This code can be freely used and redistributed. However if you decide to   *
  20  * make modifications to the code, please contact the authors, especially     *
  21  * if you plan to publish the results obtained with such modified code.       *
  22  * Any publication of results obtained using this code must include the       *
  23  * reference to nucl-th/0504047 and the published version of it, when         *
  24  * available.                                                                 *
  25  *                                                                            *
  26  *****************************************************************************/
  27 #include <TMath.h>
  28 #include "THGlobal.h"
  29 #include "ReadPar.h"
  30 #include "Parser.h"
  31 #include "DecayChannel.h"
  32 #include <sstream>
  33 #include <cstring>
  34
  35 extern ReadPar *sRPInstance;
  36
  37 const double factorials[7] = {1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0 };
  38
  39 Parser::Parser(ParticleDB *aDB)
  40 {
  41   mDB = aDB;
  42   ReadParameters();
  43 }
  44
  45 Parser::~Parser()
  46 {
  47 }
  48
  49 int Parser::check(char *a,char *b)
  50 {
  51   int tIter3=0;
  52
  53   while(a[tIter3]!='\0')
  54     {
  55       if(a[tIter3]!=b[tIter3]) return 0;
  56       tIter3++;
  57     }
  58   return 1;
  59 }
  60
  61 int Parser::GetParticleCount()
  62 {
  63   return mParticleCount;
  64 }
  65
  66 char * Parser::GetParticleName(int i)
  67 {
  68   return mNameBuffer[i];
  69 }
  70
  71 double Parser::GetParticleMass(int i)
  72 {
  73   return mMassBuffer[i];
  74 }
  75
  76 double Parser::GetParticleGamma(int i)
  77 {
  78   return mGammaBuffer[i];
  79 }
  80
  81 double Parser::GetParticleSpin(int i)
  82 {
  83   return mSpinBuffer[i];
  84 }
  85
  86 int Parser::GetParticleBarionN(int i)
  87 {
  88   return mBarionBuffer[i];
  89 }
  90
  91 int Parser::GetParticleStrangeN(int i)
  92 {
  93   return mStrangeBuffer[i];
  94 }
  95
  96 double Parser::GetParticleI3(int i)
  97 {
  98   return mI3Buffer[i];
  99 }
 100
 101 int Parser::GetParticleDecayChannelCount(int i,int j)
 102 {
 103   return mDecayBuffer[i][j];
 104 }
 105
 106 int Parser::GetParticleNumber(int i)
 107 {
 108   return mTypeCountBuffer[i];
 109 }
 110
 111 void Parser::ReadInput()
 112 {
 113   int j,tPartIter=0,l,tIter2, tIter; //variables
 114   char str[50];
 115   char str1[20];
 116   double spin1,spin2,value;
 117
 118   //////  START OF "TABLES.M" /////////
 119   ifstream in("tables.m");
 120   if(in)
 121     {
 122       //START OF HEAD-LINE
 123       in.ignore(100,'\n');
 124       //END OF HEAD-LINE
 125
 126       for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
 127
 128       tPartIter=0;
 129
 130       //START OF DATA
 131       while(in>>str)
 132         {
 133           if(*str == 'x') break;
 134
 135           l=0;
 136           //PARTICLE NAME AND MASS
 137           if(str[0] == 'm' && str[1] == '[')
 138             {
 139               for(;;)
 140                 {
 141                   if(str[l+2] == ']') break;
 142                   mNameBuffer[tPartIter][l]=str[l+2];   //name
 143                   l++;
 144                 }
 145               //                mNameBuffer[tPartIter][l]='\0';
 146               in>>str;  // sign "="
 147               in>>value;        // mass
 148               mMassBuffer[tPartIter]=value;
 149               mTypeCountBuffer[tPartIter]=tPartIter;
 150               for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
 151             }
 152
 153           if(str[0] == 'G' && str[3] == '[')
 154             {
 155               in>>str;  // sign "="
 156               in>>value;        // gamma
 157               mGammaBuffer[tPartIter]=value;
 158               for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
 159             }
 160
 161           // spin
 162           if(str[0] == 'J' && str[1] == '[')
 163             {
 164               in>>str;  // sign "="
 165               in>>str;
 166               if(str[1] == '/')
 167                 {
 168                   *str1=str[0];spin1=atof(str1);
 169                   *str1=str[2];spin2=atof(str1);
 170                   mSpinBuffer[tPartIter]=spin1/spin2;
 171                 }
 172               if(str[0] == '-')
 173                 {
 174                   *str1=str[1];spin1=atof(str1);
 175                   mSpinBuffer[tPartIter]=-spin1;
 176                 }
 177               if(str[0]!='-' && str[1]!='/')
 178                 {
 179                   *str1=str[0];
 180                   spin1=atof(str1);
 181                   mSpinBuffer[tPartIter]=spin1;
 182                 }
 183               tPartIter++;      //next particle
 184               for(tIter2=0;tIter2<50;tIter2++) str[tIter2]='\0';
 185             }
 186         }
 187       //END OF DATA
 188     }
 189   mParticleCount=tPartIter;
 190   //////  END OF "TABLES.M" /////////
 191
 192
 193   //////        START OF "i200STAR.m"//////
 194
 195   double izospin1,izospin2;     //help to calculate izospin, if niecalkowity
 196   //    input=fopen("i200STAR.m","r+");
 197   ifstream in1("i200STAR.m");
 198
 199   for(tIter2=0;tIter2<369;tIter2++)
 200     for(int jj=0;jj<2;jj++) mDecayBuffer[tIter2][jj]=0;
 201
 202   for(;;)
 203     {
 204       j=3; //first letter of particle in line in file
 205       for(tIter=0;tIter<50;tIter++) str[tIter]='\0';
 206       for(tIter=0;tIter<20;tIter++) str1[tIter]='\0';
 207       tIter=0;
 208
 209       in1.getline(str,200);
 210
 211       if(str[0] == 'x') break;
 212
 213       if(str[0] == 'f' && str[1] == 'i')
 214         {
 215           //name
 216           for(;;)
 217             {
 218               if(str[j] == ',')
 219                 {
 220                   j++;  //to skip ","
 221                   break;
 222                 }
 223               str1[tIter++]=str[j++];
 224             }
 225           for(tIter=0;tIter<369;tIter++)
 226             {
 227               if(check(str1,mNameBuffer[tIter]))
 228                 {
 229                   //barion number
 230                   for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
 231                   if(str[j] == '-')
 232                     {
 233                       *str1=str[j+1];
 234                       mBarionBuffer[tIter]=atoi(str1);
 235                       mBarionBuffer[tIter]=-mBarionBuffer[tIter];
 236                     }
 237                   if(str[j] != '-')
 238                     {
 239                       *str1=str[j];
 240                       mBarionBuffer[tIter]=atoi(str1);
 241                     }
 242                   if(str[j] == '-') j+=3;
 243                   else j+=2;
 244
 245                   //strange number
 246                   for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
 247                   if(str[j] == '-')
 248                     {
 249                       *str1=str[j+1];
 250                       mStrangeBuffer[tIter]=atoi(str1);
 251                       mStrangeBuffer[tIter]=-mStrangeBuffer[tIter];
 252                     }
 253                   if(str[j] != '-')
 254                     {
 255                       *str1=str[j];
 256                       mStrangeBuffer[tIter]=atoi(str1);
 257                     }
 258                   if(str[j] == '-') j+=3;
 259                   else j+=2;
 260
 261                   //izospin3 number
 262                   for(tIter2=0;tIter2<20;tIter2++) str1[tIter2]='\0';
 263                   if(str[j+1] == '/' && str[j+3] == ']')
 264                     {
 265                       *str1=str[j];
 266                       izospin1=atoi(str1);
 267                       *str1=str[j+2];
 268                       izospin2=atoi(str1);
 269                       mI3Buffer[tIter]=izospin1/izospin2;
 270                     }
 271                   if(str[j] == '-' && str[j+2] == '/' && str[j+4] == ']')
 272                     {
 273                       *str1=str[j+1];
 274                       izospin1=atoi(str1);
 275                       *str1=str[j+3];
 276                       izospin2=atoi(str1);
 277                       mI3Buffer[tIter]=-izospin1/izospin2;
 278                     }
 279                   if(str[j] == '-' && str[j+2] == ']')
 280                     {
 281                       *str1=str[j+1];
 282                       izospin1=atoi(str1);
 283                       mI3Buffer[tIter]=-izospin1;
 284                     }
 285                   if(str[j+1] == ']')
 286                     {
 287                       *str1=str[j];
 288                       mI3Buffer[tIter]=atof(str1);
 289                     }
 290                   break;
 291                 }
 292             }
 293         }
 294
 295       //DECAY CHANNELS
 296
 297       tIter=0;
 298       //TWO-BODY DECAYS
 299       if(str[0] == 's' && str[1] == 'e' && str[2] == '[')
 300         {
 301           // Reading in the decay channels
 302           char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tRBracket;
 303
 304           tLBrackert = strchr(str,'[');
 305           tFirstComma = strchr(str,',');
 306           tSecondComma = strchr(tFirstComma+1,',');
 307           tThirdComma = strchr(tSecondComma+1,',');
 308           tRBracket = strchr(tThirdComma,']');
 309
 310           if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
 311             PRINT_DEBUG_1("Malformed line!: " << str);
 312
 313           char *tFather = new char[tFirstComma-tLBrackert];
 314           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 315           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 316           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 317           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 318           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 319           char *tBRatio = new char[tRBracket-tThirdComma];
 320           strncpy(tBRatio, tThirdComma+1,  tRBracket-tThirdComma-1);
 321
 322           // Getting the ratio
 323           char *tMiddle, *tRatComponent;
 324           double tRatio = 1.0;
 325
 326           tMiddle = strchr(tBRatio,'_');
 327           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 328           while(tMiddle)
 329             {
 330               tRatComponent = new char[tMiddle-tBRatio+1];
 331               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 332               if (strchr(tRatComponent,'/'))
 333                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 334               else
 335                 tRatio *= atof(tRatComponent);
 336               tBRatio = tMiddle+1;
 337               tMiddle = strchr(tBRatio,'_');
 338               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 339               delete [] tRatComponent;
 340             }
 341           if (strchr(tBRatio,'/'))
 342             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 343           else
 344             tRatio *= atof(tBRatio);
 345         }
 346
 347       //THREE-BODY DECAYS
 348       j++;      //because se3[]
 349       if(str[0] == 's' && str[1] == 'e' && str[2] == '3' && str[3] == '[')
 350         {
 351           for(;;)
 352             {
 353               if(str[j] == ',')
 354                 {
 355                   j++;  //to skip ","
 356                   break;
 357                 }
 358               str1[tIter++]=str[j++];   //name
 359             }
 360
 361           for(tIter=0;tIter<369;tIter++)
 362             {
 363               if(check(str1,mNameBuffer[tIter]))
 364                 {
 365                   mDecayBuffer[tIter][1]++;
 366                   break;
 367                 }
 368             }
 369         }
 370
 371
 372     }
 373   //////        END OF "i200STAR.m"//////
 374
 375   ParticleType *tPartBuf;
 376   int tNum, tPart;
 377
 378   for(tPart=0;tPart<mParticleCount;tPart++)
 379     {
 380       tPartBuf = new ParticleType();
 381       tPartBuf->SetName(mNameBuffer[tPart]);
 382       tPartBuf->SetMass(mMassBuffer[tPart]);
 383       tPartBuf->SetGamma(mGammaBuffer[tPart]);
 384       tPartBuf->SetSpin(mSpinBuffer[tPart]);
 385       tPartBuf->SetBarionN(mBarionBuffer[tPart]);
 386       tPartBuf->SetStrangeness(mStrangeBuffer[tPart]);
 387       tPartBuf->SetI3(mI3Buffer[tPart]);
 388       tPartBuf->SetDecayChannelCount2(mDecayBuffer[tPart][0]);
 389       tPartBuf->SetDecayChannelCount3(mDecayBuffer[tPart][1]);
 390       tPartBuf->SetNumber(tPart);
 391       tNum = mDB->AddParticleType(tPartBuf);
 392     }
 393
 394   ifstream in2("i200STAR.m");
 395   while (in2.getline(str,200))
 396     {
 397       tIter=0;
 398       //TWO-BODY DECAYS
 399       if(str[0] == 's' && str[1] == 'e' && str[2] == '[')
 400         {
 401           // Reading in the decay channels
 402           char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tRBracket;
 403
 404           tLBrackert = strchr(str,'[');
 405           tFirstComma = strchr(str,',');
 406           tSecondComma = strchr(tFirstComma+1,',');
 407           tThirdComma = strchr(tSecondComma+1,',');
 408           tRBracket = strchr(tThirdComma,']');
 409           if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
 410             PRINT_DEBUG_1("Malformed line!: " << str);
 411
 412           char *tFather = new char[tFirstComma-tLBrackert];
 413           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 414           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 415           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 416           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 417           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 418           char *tBRatio = new char[tRBracket-tThirdComma];
 419           strncpy(tBRatio, tThirdComma+1,  tRBracket-tThirdComma-1);
 420
 421           // Getting the ratio
 422           char *tMiddle, *tRatComponent;
 423           double tRatio = 1.0;
 424
 425           tMiddle = strchr(tBRatio,'_');
 426           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 427           while(tMiddle)
 428             {
 429               tRatComponent = new char[tMiddle-tBRatio];
 430               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 431               if (strchr(tRatComponent,'/'))
 432                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 433               else
 434                 tRatio *= atof(tRatComponent);
 435               tBRatio = tMiddle+1;
 436               tMiddle = strchr(tBRatio,'_');
 437               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 438               delete [] tRatComponent;
 439             }
 440           if (strchr(tBRatio,'/'))
 441             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 442           else
 443             tRatio *= atof(tBRatio);
 444
 445           DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
 446           //      if (mDB->GetParticleType(tDaughter1)->GetMass() +
 447           //          mDB->GetParticleType(tDaughter2)->GetMass()
 448           //          < mDB->GetParticleType(tFather)->GetMass()) {
 449           //        (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 450           //      }
 451           (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 452         }
 453
 454       if(str[0] == 's' && str[1] == 'e' && str[2] == '3')
 455         {
 456           // Reading in the decay channels
 457           char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tFourthComma, *tRBracket;
 458
 459           tLBrackert = strchr(str,'[');
 460           tFirstComma = strchr(str,',');
 461           tSecondComma = strchr(tFirstComma+1,',');
 462           tThirdComma = strchr(tSecondComma+1,',');
 463           tFourthComma = strchr(tThirdComma+1,',');
 464           tRBracket = strchr(tThirdComma,']');
 465
 466           if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tFourthComma && tRBracket))
 467             PRINT_DEBUG_1("Malformed line!: " << str);
 468
 469           char *tFather = new char[tFirstComma-tLBrackert];
 470           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 471           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 472           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 473           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 474           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 475           char *tDaughter3 = new char[tFourthComma-tThirdComma];
 476           strncpy(tDaughter3, tThirdComma+1,  tFourthComma-tThirdComma-1);
 477           char *tBRatio = new char[tRBracket-tFourthComma];
 478           strncpy(tBRatio, tFourthComma+1, tRBracket-tFourthComma-1);
 479
 480           // Getting the ratio
 481           char *tMiddle, *tRatComponent;
 482           double tRatio = 1.0;
 483
 484           tMiddle = strchr(tBRatio,'_');
 485           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 486           while(tMiddle)
 487             {
 488               tRatComponent = new char[tMiddle-tBRatio+1];
 489               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 490               if (strchr(tRatComponent,'/'))
 491                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 492               else
 493                 tRatio *= atof(tRatComponent);
 494               tBRatio = tMiddle+1;
 495               tMiddle = strchr(tBRatio,'_');
 496               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 497               delete [] tRatComponent;
 498             }
 499           if (strchr(tBRatio,'/'))
 500             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 501           else
 502             tRatio *= atof(tBRatio);
 503
 504           DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
 505           //      if (mDB->GetParticleType(tDaughter1)->GetMass() +
 506           //          mDB->GetParticleType(tDaughter2)->GetMass() +
 507           //          mDB->GetParticleType(tDaughter3)->GetMass()
 508           //          < mDB->GetParticleType(tFather)->GetMass())
 509           (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 510         }
 511     }
 512
 513   ifstream in3("pdgcodes.m");
 514   while (in3.getline(str,200))
 515     {
 516       string tName;
 517       int tCode;
 518
 519       std::stringstream tIS(str);
 520       tIS >> tName >> tCode;
 521       mDB->GetParticleType(tName)->SetPDGCode(tCode);
 522     }
 523
 524 }
 525
 526
 527 void Parser::ReadShare()
 528 {
 529   char str[50];
 530   char str1[200];
 531
 532   ParticleType *tPartBuf;
 533   int tNum;
 534
 535   PRINT_DEBUG_1("Reading from |"<<(mInputDir+"/"+"particles.data").Data()<<"|");
 536   ifstream in((mInputDir+"/"+"particles.data").Data());
 537
 538   int number=0;
 539   if ((in) && (in.is_open()))
 540     {
 541       //START OF HEAD-LINE
 542       in.ignore(200,'\n');
 543       in.ignore(200,'\n');
 544       in.ignore(200,'\n');
 545       //END OF HEAD-LINE
 546
 547       while (in>>str)
 548         {
 549           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 550             {
 551               in.getline(str1,200);
 552               continue;
 553             }
 554           double mass, gamma, spin, I3, I, q, s, aq, as, c, ac, mc;
 555
 556           in>>mass>>gamma>>spin>>I>>I3>>q>>s>>aq>>as>>c>>ac>>mc;
 557           number++;
 558           PRINT_DEBUG_2(number<<" "<<str<<" "<<mass<<" "<<gamma<<" "<<spin<<" "<<I<<" "<<I3<<" "<<q<<" "<<aq<<" "<<s<<" "<<as<<" "<<c<<" "<<ac<<" "<<mc);
 559
 560           tPartBuf = new ParticleType();
 561           tPartBuf->SetName(str);
 562           tPartBuf->SetMass(mass);
 563           tPartBuf->SetGamma(gamma);
 564           tPartBuf->SetSpin(spin);
 565           tPartBuf->SetBarionN((int) ((q+s+c)/3. - (aq+as+ac)/3.) );
 566           tPartBuf->SetCharmN((int) (c - ac));
 567           tPartBuf->SetStrangeness((int) (as-s));
 568           tPartBuf->SetI(I);
 569           tPartBuf->SetI3(I3);
 570           tPartBuf->SetPDGCode((int) mc);
 571           tPartBuf->SetNumber(number);
 572           tNum = mDB->AddParticleType(tPartBuf);
 573         }
 574       in.close();
 575     }
 576   else
 577     {
 578       PRINT_MESSAGE("File "<<(mInputDir+"/"+"particles.data").Data()<<" containing particle data not found!");
 579       PRINT_MESSAGE("Please set the correct path to this file in the input parameter file");
 580       PRINT_MESSAGE("Aborting!");
 581       exit(0);
 582     }
 583
 584   ifstream in2((mInputDir+"/decays.data").Data());
 585   if ((in2) && (in2.is_open()))
 586     {
 587       //START OF HEAD-LINE
 588       in2.ignore(200,'\n');
 589       in2.ignore(200,'\n');
 590       in2.ignore(200,'\n');
 591       //END OF HEAD-LINE
 592
 593       char tFather[20], tDaughter1[20], tDaughter2[20], tDaughter3[20];
 594       double tBRatio, tRatio;
 595       int CGcoeff; // complete branching ratio by Clebsch-Gordan coefficient: 0-no 1-yes
 596
 597       while (in2>>str)
 598         {
 599           if (*str == '#')
 600             {
 601               in2.getline(str1,200);
 602               continue;
 603             }
 604           in2>>tDaughter1>>tDaughter2>>tDaughter3;
 605           if (!mDB->ExistsParticleType(tDaughter1)) {
 606             PRINT_MESSAGE("Did not find the daughter 1 particle: " << tDaughter1);
 607             PRINT_MESSAGE("Not adding channel");
 608             in2.getline(str1,200);
 609             continue;
 610           }
 611           if (!mDB->ExistsParticleType(tDaughter2)) {
 612             PRINT_MESSAGE("Did not find the daughter 2 particle: " << tDaughter2);
 613             PRINT_MESSAGE("Not adding channel");
 614             in2.getline(str1,200);
 615             continue;
 616           }
 617           if ((*tDaughter3>65)&&(*tDaughter3<122)&&(!mDB->ExistsParticleType(tDaughter3))) {
 618             PRINT_MESSAGE("Did not find the daughter 3 particle: " << tDaughter3);
 619             PRINT_MESSAGE("Not adding channel");
 620             in2.getline(str1,200);
 621             continue;
 622           }
 623
 624           strcpy(tFather,str);
 625           PRINT_DEBUG_2(tFather<<"\t"<<tDaughter1<<"\t"<<tDaughter2<<"\t");
 626           if ((*tDaughter3>65)&&(*tDaughter3<122)) // check if first char is a letter - if yes then 3-body decay
 627             {
 628               in2>>tBRatio>>CGcoeff;
 629               PRINT_DEBUG_2(tDaughter3<<" (3-body decay)\t");
 630               if (mDB->ExistsParticleType(tFather)) {
 631                 mDB->GetParticleType(tFather)->SetDecayChannelCount3(mDB->GetParticleType(tFather)->GetDecayChannelCount3()+1);
 632
 633                 tRatio=tBRatio;
 634                 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
 635                 if (mDB->GetParticleType(tDaughter1)->GetMass() +
 636                     mDB->GetParticleType(tDaughter2)->GetMass() +
 637                     mDB->GetParticleType(tDaughter3)->GetMass()
 638                     < mDB->GetParticleType(tFather)->GetMass())
 639                   (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 640
 641                 tRatio=tBRatio;
 642                 PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
 643               }
 644               else {
 645                 PRINT_MESSAGE("Did not find the father particle: " << tFather);
 646                 PRINT_MESSAGE("Not adding channel");
 647               }
 648             }
 649           else // 2-body decay
 650             {
 651               tBRatio=atof(tDaughter3);
 652
 653               in2>>CGcoeff;
 654               PRINT_DEBUG_2(" (2-body decay)\t");
 655               if (mDB->ExistsParticleType(tFather)) {
 656                 mDB->GetParticleType(tFather)->SetDecayChannelCount2(mDB->GetParticleType(tFather)->GetDecayChannelCount2()+1);
 657
 658                 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 659                 if (CGcoeff) // complete branching ratio by Clebsch-Gordan coefficient
 660                   {
 661                     double j1, m1, j2, m2, J, M, CB;
 662                     J=mDB->GetParticleType(tFather)->GetI();
 663                     M=mDB->GetParticleType(tFather)->GetI3();
 664                     j1=mDB->GetParticleType(tDaughter1)->GetI();
 665                     m1=mDB->GetParticleType(tDaughter1)->GetI3();
 666                     j2=mDB->GetParticleType(tDaughter2)->GetI();
 667                     m2=mDB->GetParticleType(tDaughter2)->GetI3();
 668                     PRINT_DEBUG_2(" "<<J<<" "<<M<<" "<<j1<<" "<<m1<<" "<<j2<<" "<<m2);
 669
 670                     CB = ClebschGordan(J, M, j1, m1, j2, m2);
 671                     tRatio = CB*CB * tBRatio;
 672
 673                     // Multiply the Clebsh by two?
 674                     // The same spin, mass, strangeness
 675                     // and different I3?
 676                     if ((fabs(mDB->GetParticleType(tDaughter1)->GetSpin() - mDB->GetParticleType(tDaughter2)->GetSpin()) < 0.01) &&
 677                         (fabs(mDB->GetParticleType(tDaughter1)->GetMass() - mDB->GetParticleType(tDaughter2)->GetMass()) < 0.01) &&
 678                         (mDB->GetParticleType(tDaughter1)->GetStrangeness() == mDB->GetParticleType(tDaughter2)->GetStrangeness()) &&
 679                         (fabs(mDB->GetParticleType(tDaughter1)->GetI3() - mDB->GetParticleType(tDaughter2)->GetI3()) > 0.01))
 680                       {
 681                         PRINT_DEBUG_2("Multuplying Clebsch by two for " << tFather << "->" << tDaughter1 << "+" << tDaughter2);
 682                         tRatio *= 2.0;
 683                       }
 684
 685                     PRINT_DEBUG_2(CB << '\t' << tBRatio << '\t' << tRatio<<"\t"<<CGcoeff);
 686                   }
 687
 688                 else
 689                   {
 690                     tRatio=tBRatio;
 691                     PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
 692                   }
 693                 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
 694                 if (mDB->GetParticleType(tDaughter1)->GetMass() +
 695                     mDB->GetParticleType(tDaughter2)->GetMass()
 696                     < mDB->GetParticleType(tFather)->GetMass())
 697                   {
 698                     (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 699                     PRINT_DEBUG_2("Added channel " << newChannel << " " << mDB->GetParticleTypeIndex(tFather) << " " << mDB->GetParticleTypeIndex(tDaughter1) << " " << mDB->GetParticleTypeIndex(tDaughter2));
 700                   }
 701                 else
 702                   {
 703
 704                     PRINT_DEBUG_2("Masses do not match! Not adding channel " << newChannel);
 705                   }
 706               }
 707               else {
 708                 PRINT_MESSAGE("Did not find the father particle: " << tFather);
 709                 PRINT_MESSAGE("Not adding channel");
 710               }
 711             }
 712         }
 713       in2.close();
 714     }
 715   else {
 716     PRINT_MESSAGE("File "<<(mInputDir+"/decays.data").Data()<<" with particle decay channels not found!");
 717     PRINT_MESSAGE("No particle decays will be simulated");
 718   }
 719 }
 720
 721 double
 722 Parser::ClebschGordan(double aJot, double aEm, double aJot1, double aEm1, double aJot2, double aEm2)
 723 {
 724   int mint, maxt;
 725   double cgc = 0.0;
 726   int titer;
 727   double coef;
 728
 729   maxt = lrint(aJot1 + aJot2 - aJot);
 730   mint = 0;
 731   if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
 732   if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
 733   if (lrint(-(aJot-aJot2+aEm1)) > mint) mint = lrint(-(aJot-aJot2+aEm1));
 734   if (lrint(-(aJot-aJot1-aEm2)) > mint) mint = lrint(-(aJot-aJot1-aEm2));
 735
 736   PRINT_DEBUG_3("mint " << mint << " " <<  aJot1 << " " << aEm1);
 737   PRINT_DEBUG_3("maxt " << maxt << " " <<  aJot2 << " " << aEm2);
 738
 739   for (titer = mint; titer<=maxt; titer ++)
 740     {
 741       coef = TMath::Power(-1, titer);
 742       PRINT_DEBUG_3("coef1 " << coef);
 743       coef *= TMath::Sqrt((2*aJot+1)*
 744                           factorials[lrint(aJot1+aEm1)] *
 745                           factorials[lrint(aJot1-aEm1)] *
 746                           factorials[lrint(aJot2+aEm2)] *
 747                           factorials[lrint(aJot2-aEm2)] *
 748                           factorials[lrint(aJot+aEm)] *
 749                           factorials[lrint(aJot-aEm)]);
 750       PRINT_DEBUG_3("coef2 " << coef);
 751       coef /= (factorials[titer] *
 752                factorials[lrint(aJot1+aJot2-aJot-titer)] *
 753                factorials[lrint(aJot1-aEm1-titer)] *
 754                factorials[lrint(aJot2+aEm2-titer)] *
 755                factorials[lrint(aJot-aJot2+aEm1+titer)] *
 756                factorials[lrint(aJot-aJot1-aEm2+titer)]);
 757       PRINT_DEBUG_3("coef3 " << coef);
 758
 759       cgc += coef;
 760     }
 761
 762   cgc *= DeltaJ(aJot1, aJot2, aJot);
 763
 764   return cgc;
 765 }
 766
 767 double
 768 Parser::DeltaJ(double aJot1, double aJot2, double aJot)
 769 {
 770   double res = TMath::Sqrt(1.0 *
 771                            factorials[lrint(aJot1+aJot2-aJot)] *
 772                            factorials[lrint(aJot1-aJot2+aJot)] *
 773                            factorials[lrint(-aJot1+aJot2+aJot)] /
 774                            factorials[lrint(aJot1+aJot2+aJot+1)]);
 775
 776   return res;
 777 }
 778
 779 void
 780 Parser::ReadParameters()
 781 {
 782   // Read the input directory
 783   try {
 784     mInputDir = sRPInstance->getPar("InputDirSHARE");
 785   }
 786   catch (STR tError) {
 787     PRINT_DEBUG_1("Parser::ReadParameters - Caught exception " << tError);
 788     PRINT_MESSAGE("Did not find SHARE input file location.");
 789     PRINT_MESSAGE("Using default: '../share'");
 790     mInputDir = "../share";
 791   }
 792 }