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[200];
 115   char str1[200];
 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=NULL, *tFirstComma=NULL, *tSecondComma=NULL, *tThirdComma=NULL, *tRBracket=NULL;
 303
 304           tLBrackert = strchr(str,'[');
 305           tFirstComma = strchr(str,',');
 306           if (!tFirstComma) exit(0);
 307           tSecondComma = strchr(tFirstComma+1,',');
 308           if (!tSecondComma) exit(0);
 309           tThirdComma = strchr(tSecondComma+1,',');
 310           if (!tThirdComma) exit(0);
 311           tRBracket = strchr(tThirdComma,']');
 312
 313           if (!((tLBrackert != NULL) && (tFirstComma != NULL) && (tSecondComma != NULL) && ( tThirdComma != NULL) && (tRBracket!= NULL)))
 314             PRINT_DEBUG_1("Malformed line!: " << str);
 315
 316           char *tFather = new char[tFirstComma-tLBrackert];
 317           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 318           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 319           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 320           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 321           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 322           char *tBRatio = new char[tRBracket-tThirdComma];
 323           strncpy(tBRatio, tThirdComma+1,  tRBracket-tThirdComma-1);
 324
 325           // Getting the ratio
 326           char *tMiddle, *tRatComponent;
 327           double tRatio = 1.0;
 328
 329           tMiddle = strchr(tBRatio,'_');
 330           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 331           while(tMiddle)
 332             {
 333               tRatComponent = new char[tMiddle-tBRatio+1];
 334               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 335               if (strchr(tRatComponent,'/'))
 336                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 337               else
 338                 tRatio *= atof(tRatComponent);
 339               tBRatio = tMiddle+1;
 340               tMiddle = strchr(tBRatio,'_');
 341               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 342               delete [] tRatComponent;
 343             }
 344           if (strchr(tBRatio,'/'))
 345             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 346           else
 347             tRatio *= atof(tBRatio);
 348
 349           delete [] tFather;
 350           delete [] tDaughter1;
 351           delete [] tDaughter2;
 352           delete [] tBRatio;
 353         }
 354
 355       //THREE-BODY DECAYS
 356       j++;      //because se3[]
 357       if(str[0] == 's' && str[1] == 'e' && str[2] == '3' && str[3] == '[')
 358         {
 359           for(;;)
 360             {
 361               if(str[j] == ',')
 362                 {
 363                   j++;  //to skip ","
 364                   break;
 365                 }
 366               str1[tIter++]=str[j++];   //name
 367             }
 368
 369           for(tIter=0;tIter<369;tIter++)
 370             {
 371               if(check(str1,mNameBuffer[tIter]))
 372                 {
 373                   mDecayBuffer[tIter][1]++;
 374                   break;
 375                 }
 376             }
 377         }
 378
 379
 380     }
 381   //////        END OF "i200STAR.m"//////
 382
 383   ParticleType *tPartBuf;
 384   int tNum, tPart;
 385
 386   for(tPart=0;tPart<mParticleCount;tPart++)
 387     {
 388       tPartBuf = new ParticleType();
 389       tPartBuf->SetName(mNameBuffer[tPart]);
 390       tPartBuf->SetMass(mMassBuffer[tPart]);
 391       tPartBuf->SetGamma(mGammaBuffer[tPart]);
 392       tPartBuf->SetSpin(mSpinBuffer[tPart]);
 393       tPartBuf->SetBarionN(mBarionBuffer[tPart]);
 394       tPartBuf->SetStrangeness(mStrangeBuffer[tPart]);
 395       tPartBuf->SetI3(mI3Buffer[tPart]);
 396       tPartBuf->SetDecayChannelCount2(mDecayBuffer[tPart][0]);
 397       tPartBuf->SetDecayChannelCount3(mDecayBuffer[tPart][1]);
 398       tPartBuf->SetNumber(tPart);
 399       tNum = mDB->AddParticleType(tPartBuf);
 400     }
 401
 402   ifstream in2("i200STAR.m");
 403   while (in2.getline(str,200))
 404     {
 405       tIter=0;
 406       //TWO-BODY DECAYS
 407       if(str[0] == 's' && str[1] == 'e' && str[2] == '[')
 408         {
 409           // Reading in the decay channels
 410           char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tRBracket;
 411
 412           tLBrackert = strchr(str,'[');
 413           if (!tLBrackert) exit(0);
 414           tFirstComma = strchr(str,',');
 415           if (!tFirstComma) exit(0);
 416           tSecondComma = strchr(tFirstComma+1,',');
 417           if (!tSecondComma) exit(0);
 418           tThirdComma = strchr(tSecondComma+1,',');
 419           if (!tThirdComma) exit(0);
 420           tRBracket = strchr(tThirdComma,']');
 421           if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
 422             PRINT_DEBUG_1("Malformed line!: " << str);
 423
 424           char *tFather = new char[tFirstComma-tLBrackert];
 425           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 426           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 427           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 428           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 429           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 430           char *tBRatio = new char[tRBracket-tThirdComma];
 431           strncpy(tBRatio, tThirdComma+1,  tRBracket-tThirdComma-1);
 432
 433           // Getting the ratio
 434           char *tMiddle, *tRatComponent;
 435           double tRatio = 1.0;
 436
 437           tMiddle = strchr(tBRatio,'_');
 438           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 439           while(tMiddle)
 440             {
 441               tRatComponent = new char[tMiddle-tBRatio];
 442               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 443               if (strchr(tRatComponent,'/'))
 444                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 445               else
 446                 tRatio *= atof(tRatComponent);
 447               tBRatio = tMiddle+1;
 448               tMiddle = strchr(tBRatio,'_');
 449               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 450               delete [] tRatComponent;
 451             }
 452           if (strchr(tBRatio,'/'))
 453             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 454           else
 455             tRatio *= atof(tBRatio);
 456
 457           DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
 458           //      if (mDB->GetParticleType(tDaughter1)->GetMass() +
 459           //          mDB->GetParticleType(tDaughter2)->GetMass()
 460           //          < mDB->GetParticleType(tFather)->GetMass()) {
 461           //        (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 462           //      }
 463           (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 464
 465           delete [] tFather;
 466           delete [] tDaughter1;
 467           delete [] tDaughter2;
 468           delete [] tBRatio;
 469
 470           delete newChannel;
 471         }
 472
 473       if(str[0] == 's' && str[1] == 'e' && str[2] == '3')
 474         {
 475           // Reading in the decay channels
 476           char *tLBrackert, *tFirstComma, *tSecondComma, *tThirdComma, *tFourthComma, *tRBracket;
 477
 478           tLBrackert = strchr(str,'[');
 479           tFirstComma = strchr(str,',');
 480           tSecondComma = strchr(tFirstComma+1,',');
 481           tThirdComma = strchr(tSecondComma+1,',');
 482           tFourthComma = strchr(tThirdComma+1,',');
 483           tRBracket = strchr(tThirdComma,']');
 484
 485           if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tFourthComma && tRBracket))
 486             PRINT_DEBUG_1("Malformed line!: " << str);
 487
 488           char *tFather = new char[tFirstComma-tLBrackert];
 489           strncpy(tFather, tLBrackert+1,   tFirstComma-tLBrackert-1);
 490           char *tDaughter1 = new char[tSecondComma-tFirstComma];
 491           strncpy(tDaughter1, tFirstComma+1,  tSecondComma-tFirstComma-1);
 492           char *tDaughter2 = new char[tThirdComma-tSecondComma];
 493           strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
 494           char *tDaughter3 = new char[tFourthComma-tThirdComma];
 495           strncpy(tDaughter3, tThirdComma+1,  tFourthComma-tThirdComma-1);
 496           char *tBRatio = new char[tRBracket-tFourthComma];
 497           strncpy(tBRatio, tFourthComma+1, tRBracket-tFourthComma-1);
 498
 499           // Getting the ratio
 500           char *tMiddle, *tRatComponent;
 501           double tRatio = 1.0;
 502
 503           tMiddle = strchr(tBRatio,'_');
 504           if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 505           while(tMiddle)
 506             {
 507               tRatComponent = new char[tMiddle-tBRatio+1];
 508               strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
 509               if (strchr(tRatComponent,'/'))
 510                 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
 511               else
 512                 tRatio *= atof(tRatComponent);
 513               tBRatio = tMiddle+1;
 514               tMiddle = strchr(tBRatio,'_');
 515               if (!tMiddle) tMiddle = strchr(tBRatio,' ');
 516               delete [] tRatComponent;
 517             }
 518           if (strchr(tBRatio,'/'))
 519             tRatio *= atof(tBRatio)/atof(tBRatio+2);
 520           else
 521             tRatio *= atof(tBRatio);
 522
 523           DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
 524           //      if (mDB->GetParticleType(tDaughter1)->GetMass() +
 525           //          mDB->GetParticleType(tDaughter2)->GetMass() +
 526           //          mDB->GetParticleType(tDaughter3)->GetMass()
 527           //          < mDB->GetParticleType(tFather)->GetMass())
 528           (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 529
 530           delete [] tFather;
 531           delete [] tDaughter1;
 532           delete [] tDaughter2;
 533           delete [] tDaughter3;
 534           delete [] tBRatio;
 535
 536           delete newChannel;
 537         }
 538     }
 539
 540   ifstream in3("pdgcodes.m");
 541   while (in3.getline(str,200))
 542     {
 543       string tName;
 544       int tCode;
 545
 546       std::stringstream tIS(str);
 547       tIS >> tName >> tCode;
 548       mDB->GetParticleType(tName)->SetPDGCode(tCode);
 549     }
 550
 551 }
 552
 553
 554 void Parser::ReadShare()
 555 {
 556   char str[50];
 557   char str1[200];
 558
 559   ParticleType *tPartBuf;
 560   int tNum;
 561
 562   PRINT_DEBUG_1("Reading from |"<<(mInputDir+"/"+"particles.data").Data()<<"|");
 563   ifstream in((mInputDir+"/"+"particles.data").Data());
 564
 565   int number=0;
 566   if ((in) && (in.is_open()))
 567     {
 568       //START OF HEAD-LINE
 569       in.ignore(200,'\n');
 570       in.ignore(200,'\n');
 571       in.ignore(200,'\n');
 572       //END OF HEAD-LINE
 573
 574       while (in>>str)
 575         {
 576           if (/*(*str == '#')||*/(*str<65)||(*str>122))
 577             {
 578               in.getline(str1,200);
 579               continue;
 580             }
 581           double mass, gamma, spin, I3, I, q, s, aq, as, c, ac, mc;
 582
 583           in>>mass>>gamma>>spin>>I>>I3>>q>>s>>aq>>as>>c>>ac>>mc;
 584           number++;
 585           PRINT_DEBUG_2(number<<" "<<str<<" "<<mass<<" "<<gamma<<" "<<spin<<" "<<I<<" "<<I3<<" "<<q<<" "<<aq<<" "<<s<<" "<<as<<" "<<c<<" "<<ac<<" "<<mc);
 586
 587           tPartBuf = new ParticleType();
 588           tPartBuf->SetName(str);
 589           tPartBuf->SetMass(mass);
 590           tPartBuf->SetGamma(gamma);
 591           tPartBuf->SetSpin(spin);
 592           tPartBuf->SetBarionN((int) ((q+s+c)/3. - (aq+as+ac)/3.) );
 593           tPartBuf->SetCharmN((int) (c - ac));
 594           tPartBuf->SetStrangeness((int) (as-s));
 595           tPartBuf->SetI(I);
 596           tPartBuf->SetI3(I3);
 597           tPartBuf->SetPDGCode((int) mc);
 598           tPartBuf->SetNumber(number);
 599           tNum = mDB->AddParticleType(tPartBuf);
 600
 601           delete tPartBuf;
 602         }
 603       in.close();
 604     }
 605   else
 606     {
 607       PRINT_MESSAGE("File "<<(mInputDir+"/"+"particles.data").Data()<<" containing particle data not found!");
 608       PRINT_MESSAGE("Please set the correct path to this file in the input parameter file");
 609       PRINT_MESSAGE("Aborting!");
 610       exit(0);
 611     }
 612
 613   ifstream in2((mInputDir+"/decays.data").Data());
 614   if ((in2) && (in2.is_open()))
 615     {
 616       //START OF HEAD-LINE
 617       in2.ignore(200,'\n');
 618       in2.ignore(200,'\n');
 619       in2.ignore(200,'\n');
 620       //END OF HEAD-LINE
 621
 622       char tFather[50], tDaughter1[50], tDaughter2[50], tDaughter3[50];
 623       double tBRatio, tRatio;
 624       int CGcoeff; // complete branching ratio by Clebsch-Gordan coefficient: 0-no 1-yes
 625
 626       while (in2>>str)
 627         {
 628           if (*str == '#')
 629             {
 630               in2.getline(str1,200);
 631               continue;
 632             }
 633           in2>>tDaughter1>>tDaughter2>>tDaughter3;
 634           if (!mDB->ExistsParticleType(tDaughter1)) {
 635             PRINT_MESSAGE("Did not find the daughter 1 particle: " << tDaughter1);
 636             PRINT_MESSAGE("Not adding channel");
 637             in2.getline(str1,200);
 638             continue;
 639           }
 640           if (!mDB->ExistsParticleType(tDaughter2)) {
 641             PRINT_MESSAGE("Did not find the daughter 2 particle: " << tDaughter2);
 642             PRINT_MESSAGE("Not adding channel");
 643             in2.getline(str1,200);
 644             continue;
 645           }
 646           if ((*tDaughter3>65)&&(*tDaughter3<122)&&(!mDB->ExistsParticleType(tDaughter3))) {
 647             PRINT_MESSAGE("Did not find the daughter 3 particle: " << tDaughter3);
 648             PRINT_MESSAGE("Not adding channel");
 649             in2.getline(str1,200);
 650             continue;
 651           }
 652
 653           strcpy(tFather,str);
 654           PRINT_DEBUG_2(tFather<<"\t"<<tDaughter1<<"\t"<<tDaughter2<<"\t");
 655           if ((*tDaughter3>65)&&(*tDaughter3<122)) // check if first char is a letter - if yes then 3-body decay
 656             {
 657               in2>>tBRatio>>CGcoeff;
 658               PRINT_DEBUG_2(tDaughter3<<" (3-body decay)\t");
 659               if (mDB->ExistsParticleType(tFather)) {
 660                 mDB->GetParticleType(tFather)->SetDecayChannelCount3(mDB->GetParticleType(tFather)->GetDecayChannelCount3()+1);
 661
 662                 tRatio=tBRatio;
 663                 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
 664                 if (mDB->GetParticleType(tDaughter1)->GetMass() +
 665                     mDB->GetParticleType(tDaughter2)->GetMass() +
 666                     mDB->GetParticleType(tDaughter3)->GetMass()
 667                     < mDB->GetParticleType(tFather)->GetMass())
 668                   (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 669
 670                 delete newChannel;
 671
 672                 tRatio=tBRatio;
 673                 PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
 674               }
 675               else {
 676                 PRINT_MESSAGE("Did not find the father particle: " << tFather);
 677                 PRINT_MESSAGE("Not adding channel");
 678               }
 679             }
 680           else // 2-body decay
 681             {
 682               tBRatio=atof(tDaughter3);
 683
 684               in2>>CGcoeff;
 685               PRINT_DEBUG_2(" (2-body decay)\t");
 686               if (mDB->ExistsParticleType(tFather)) {
 687                 mDB->GetParticleType(tFather)->SetDecayChannelCount2(mDB->GetParticleType(tFather)->GetDecayChannelCount2()+1);
 688
 689                 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 690                 if (CGcoeff) // complete branching ratio by Clebsch-Gordan coefficient
 691                   {
 692                     double j1, m1, j2, m2, J, M, CB;
 693                     J=mDB->GetParticleType(tFather)->GetI();
 694                     M=mDB->GetParticleType(tFather)->GetI3();
 695                     j1=mDB->GetParticleType(tDaughter1)->GetI();
 696                     m1=mDB->GetParticleType(tDaughter1)->GetI3();
 697                     j2=mDB->GetParticleType(tDaughter2)->GetI();
 698                     m2=mDB->GetParticleType(tDaughter2)->GetI3();
 699                     PRINT_DEBUG_2(" "<<J<<" "<<M<<" "<<j1<<" "<<m1<<" "<<j2<<" "<<m2);
 700
 701                     CB = ClebschGordan(J, M, j1, m1, j2, m2);
 702                     tRatio = CB*CB * tBRatio;
 703
 704                     // Multiply the Clebsh by two?
 705                     // The same spin, mass, strangeness
 706                     // and different I3?
 707                     if ((fabs(mDB->GetParticleType(tDaughter1)->GetSpin() - mDB->GetParticleType(tDaughter2)->GetSpin()) < 0.01) &&
 708                         (fabs(mDB->GetParticleType(tDaughter1)->GetMass() - mDB->GetParticleType(tDaughter2)->GetMass()) < 0.01) &&
 709                         (mDB->GetParticleType(tDaughter1)->GetStrangeness() == mDB->GetParticleType(tDaughter2)->GetStrangeness()) &&
 710                         (fabs(mDB->GetParticleType(tDaughter1)->GetI3() - mDB->GetParticleType(tDaughter2)->GetI3()) > 0.01))
 711                       {
 712                         PRINT_DEBUG_2("Multuplying Clebsch by two for " << tFather << "->" << tDaughter1 << "+" << tDaughter2);
 713                         tRatio *= 2.0;
 714                       }
 715
 716                     PRINT_DEBUG_2(CB << '\t' << tBRatio << '\t' << tRatio<<"\t"<<CGcoeff);
 717                   }
 718
 719                 else
 720                   {
 721                     tRatio=tBRatio;
 722                     PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
 723                   }
 724                 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
 725                 if (mDB->GetParticleType(tDaughter1)->GetMass() +
 726                     mDB->GetParticleType(tDaughter2)->GetMass()
 727                     < mDB->GetParticleType(tFather)->GetMass())
 728                   {
 729                     (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
 730                     PRINT_DEBUG_2("Added channel " << newChannel << " " << mDB->GetParticleTypeIndex(tFather) << " " << mDB->GetParticleTypeIndex(tDaughter1) << " " << mDB->GetParticleTypeIndex(tDaughter2));
 731                     delete newChannel;
 732                   }
 733                 else
 734                   {
 735
 736                     PRINT_DEBUG_2("Masses do not match! Not adding channel " << newChannel);
 737
 738                     delete newChannel;
 739                   }
 740               }
 741               else {
 742                 PRINT_MESSAGE("Did not find the father particle: " << tFather);
 743                 PRINT_MESSAGE("Not adding channel");
 744               }
 745             }
 746         }
 747       in2.close();
 748     }
 749   else {
 750     PRINT_MESSAGE("File "<<(mInputDir+"/decays.data").Data()<<" with particle decay channels not found!");
 751     PRINT_MESSAGE("No particle decays will be simulated");
 752   }
 753 }
 754
 755 double
 756 Parser::ClebschGordan(double aJot, double aEm, double aJot1, double aEm1, double aJot2, double aEm2)
 757 {
 758   int mint, maxt;
 759   double cgc = 0.0;
 760   int titer;
 761   double coef;
 762
 763   maxt = lrint(aJot1 + aJot2 - aJot);
 764   mint = 0;
 765   if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
 766   if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
 767   if (lrint(-(aJot-aJot2+aEm1)) > mint) mint = lrint(-(aJot-aJot2+aEm1));
 768   if (lrint(-(aJot-aJot1-aEm2)) > mint) mint = lrint(-(aJot-aJot1-aEm2));
 769
 770   PRINT_DEBUG_3("mint " << mint << " " <<  aJot1 << " " << aEm1);
 771   PRINT_DEBUG_3("maxt " << maxt << " " <<  aJot2 << " " << aEm2);
 772
 773   for (titer = mint; titer<=maxt; titer ++)
 774     {
 775       coef = TMath::Power(-1, titer);
 776       PRINT_DEBUG_3("coef1 " << coef);
 777       coef *= TMath::Sqrt((2*aJot+1)*
 778                           factorials[lrint(aJot1+aEm1)] *
 779                           factorials[lrint(aJot1-aEm1)] *
 780                           factorials[lrint(aJot2+aEm2)] *
 781                           factorials[lrint(aJot2-aEm2)] *
 782                           factorials[lrint(aJot+aEm)] *
 783                           factorials[lrint(aJot-aEm)]);
 784       PRINT_DEBUG_3("coef2 " << coef);
 785       coef /= (factorials[titer] *
 786                factorials[lrint(aJot1+aJot2-aJot-titer)] *
 787                factorials[lrint(aJot1-aEm1-titer)] *
 788                factorials[lrint(aJot2+aEm2-titer)] *
 789                factorials[lrint(aJot-aJot2+aEm1+titer)] *
 790                factorials[lrint(aJot-aJot1-aEm2+titer)]);
 791       PRINT_DEBUG_3("coef3 " << coef);
 792
 793       cgc += coef;
 794     }
 795
 796   cgc *= DeltaJ(aJot1, aJot2, aJot);
 797
 798   return cgc;
 799 }
 800
 801 double
 802 Parser::DeltaJ(double aJot1, double aJot2, double aJot)
 803 {
 804   double res = TMath::Sqrt(1.0 *
 805                            factorials[lrint(aJot1+aJot2-aJot)] *
 806                            factorials[lrint(aJot1-aJot2+aJot)] *
 807                            factorials[lrint(-aJot1+aJot2+aJot)] /
 808                            factorials[lrint(aJot1+aJot2+aJot+1)]);
 809
 810   return res;
 811 }
 812
 813 void
 814 Parser::ReadParameters()
 815 {
 816   // Read the input directory
 817   try {
 818     mInputDir = sRPInstance->getPar("InputDirSHARE");
 819   }
 820   catch (STR tError) {
 821     PRINT_DEBUG_1("Parser::ReadParameters - Caught exception " << tError);
 822     PRINT_MESSAGE("Did not find SHARE input file location.");
 823     PRINT_MESSAGE("Using default: '../share'");
 824     mInputDir = "../share";
 825   }
 826 }