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