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