Update master to aliroot
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Parser.cxx
CommitLineData
2e967919 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>
f8777da0 33#include <cstring>
2e967919 34
35extern ReadPar *sRPInstance;
36
37const double factorials[7] = {1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0 };
38
39Parser::Parser(ParticleDB *aDB)
40{
41 mDB = aDB;
42 ReadParameters();
43}
44
45Parser::~Parser()
46{
47}
48
49int 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
61int Parser::GetParticleCount()
62{
63 return mParticleCount;
64}
65
66char * Parser::GetParticleName(int i)
67{
68 return mNameBuffer[i];
69}
70
71double Parser::GetParticleMass(int i)
72{
73 return mMassBuffer[i];
74}
75
76double Parser::GetParticleGamma(int i)
77{
78 return mGammaBuffer[i];
79}
80
81double Parser::GetParticleSpin(int i)
82{
83 return mSpinBuffer[i];
84}
85
86int Parser::GetParticleBarionN(int i)
87{
88 return mBarionBuffer[i];
89}
90
91int Parser::GetParticleStrangeN(int i)
92{
93 return mStrangeBuffer[i];
94}
95
96double Parser::GetParticleI3(int i)
97{
98 return mI3Buffer[i];
99}
100
101int Parser::GetParticleDecayChannelCount(int i,int j)
102{
103 return mDecayBuffer[i][j];
104}
105
106int Parser::GetParticleNumber(int i)
107{
108 return mTypeCountBuffer[i];
109}
110
111void Parser::ReadInput()
112{
c61a7285 113 int j,tPartIter=0,l,tIter2, tIter; //variables
09b7c66c 114 char str[200];
115 char str1[200];
2e967919 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
f7134a95 302 char *tLBrackert=NULL, *tFirstComma=NULL, *tSecondComma=NULL, *tThirdComma=NULL, *tRBracket=NULL;
2e967919 303
304 tLBrackert = strchr(str,'[');
305 tFirstComma = strchr(str,',');
aec4187c 306 if (!tFirstComma) exit(0);
2e967919 307 tSecondComma = strchr(tFirstComma+1,',');
aec4187c 308 if (!tSecondComma) exit(0);
2e967919 309 tThirdComma = strchr(tSecondComma+1,',');
aec4187c 310 if (!tThirdComma) exit(0);
2e967919 311 tRBracket = strchr(tThirdComma,']');
312
f5e2bbd0 313 if (!((tLBrackert != NULL) && (tFirstComma != NULL) && (tSecondComma != NULL) && ( tThirdComma != NULL) && (tRBracket!= NULL))) {
2e967919 314 PRINT_DEBUG_1("Malformed line!: " << str);
f5e2bbd0 315 exit(0);
316 }
317
f8777da0 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);
2e967919 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 {
f8777da0 335 tRatComponent = new char[tMiddle-tBRatio+1];
336 strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
2e967919 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,' ');
f8777da0 344 delete [] tRatComponent;
2e967919 345 }
346 if (strchr(tBRatio,'/'))
347 tRatio *= atof(tBRatio)/atof(tBRatio+2);
348 else
349 tRatio *= atof(tBRatio);
09b7c66c 350
8e85ac49 351 delete [] tFather;
352 delete [] tDaughter1;
353 delete [] tDaughter2;
354 delete [] tBRatio;
2e967919 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;
2e967919 413
414 tLBrackert = strchr(str,'[');
37273086 415 if (!tLBrackert) exit(0);
2e967919 416 tFirstComma = strchr(str,',');
37273086 417 if (!tFirstComma) exit(0);
2e967919 418 tSecondComma = strchr(tFirstComma+1,',');
37273086 419 if (!tSecondComma) exit(0);
2e967919 420 tThirdComma = strchr(tSecondComma+1,',');
37273086 421 if (!tThirdComma) exit(0);
2e967919 422 tRBracket = strchr(tThirdComma,']');
423 if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tRBracket))
424 PRINT_DEBUG_1("Malformed line!: " << str);
425
f8777da0 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);
2e967919 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 {
f8777da0 443 tRatComponent = new char[tMiddle-tBRatio];
444 strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
2e967919 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,' ');
f8777da0 452 delete [] tRatComponent;
2e967919 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);
09b7c66c 466
0cf15e65 467 delete [] tFather;
468 delete [] tDaughter1;
469 delete [] tDaughter2;
0cf15e65 470 delete [] tBRatio;
471
09b7c66c 472 delete newChannel;
2e967919 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;
2e967919 479
480 tLBrackert = strchr(str,'[');
481 tFirstComma = strchr(str,',');
adc5ea82 482 if (!tFirstComma) exit(0);
2e967919 483 tSecondComma = strchr(tFirstComma+1,',');
adc5ea82 484 if (!tSecondComma) exit(0);
2e967919 485 tThirdComma = strchr(tSecondComma+1,',');
adc5ea82 486 if (!tThirdComma) exit(0);
2e967919 487 tFourthComma = strchr(tThirdComma+1,',');
488 tRBracket = strchr(tThirdComma,']');
489
490 if (!(tLBrackert && tFirstComma && tSecondComma && tThirdComma && tFourthComma && tRBracket))
491 PRINT_DEBUG_1("Malformed line!: " << str);
492
f8777da0 493 char *tFather = new char[tFirstComma-tLBrackert];
494 strncpy(tFather, tLBrackert+1, tFirstComma-tLBrackert-1);
495 char *tDaughter1 = new char[tSecondComma-tFirstComma];
496 strncpy(tDaughter1, tFirstComma+1, tSecondComma-tFirstComma-1);
497 char *tDaughter2 = new char[tThirdComma-tSecondComma];
498 strncpy(tDaughter2, tSecondComma+1, tThirdComma-tSecondComma-1);
499 char *tDaughter3 = new char[tFourthComma-tThirdComma];
500 strncpy(tDaughter3, tThirdComma+1, tFourthComma-tThirdComma-1);
501 char *tBRatio = new char[tRBracket-tFourthComma];
502 strncpy(tBRatio, tFourthComma+1, tRBracket-tFourthComma-1);
2e967919 503
504 // Getting the ratio
505 char *tMiddle, *tRatComponent;
506 double tRatio = 1.0;
507
508 tMiddle = strchr(tBRatio,'_');
509 if (!tMiddle) tMiddle = strchr(tBRatio,' ');
510 while(tMiddle)
511 {
f8777da0 512 tRatComponent = new char[tMiddle-tBRatio+1];
513 strncpy(tRatComponent, tBRatio, tMiddle-tBRatio);
2e967919 514 if (strchr(tRatComponent,'/'))
515 tRatio *= atof(tRatComponent)/atof(tRatComponent+2);
516 else
517 tRatio *= atof(tRatComponent);
518 tBRatio = tMiddle+1;
519 tMiddle = strchr(tBRatio,'_');
520 if (!tMiddle) tMiddle = strchr(tBRatio,' ');
f8777da0 521 delete [] tRatComponent;
2e967919 522 }
523 if (strchr(tBRatio,'/'))
524 tRatio *= atof(tBRatio)/atof(tBRatio+2);
525 else
526 tRatio *= atof(tBRatio);
527
528 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
529 // if (mDB->GetParticleType(tDaughter1)->GetMass() +
530 // mDB->GetParticleType(tDaughter2)->GetMass() +
531 // mDB->GetParticleType(tDaughter3)->GetMass()
532 // < mDB->GetParticleType(tFather)->GetMass())
533 (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
d58024dc 534
0cf15e65 535 delete [] tFather;
536 delete [] tDaughter1;
537 delete [] tDaughter2;
538 delete [] tDaughter3;
539 delete [] tBRatio;
6750cf94 540
541 delete newChannel;
2e967919 542 }
543 }
544
545 ifstream in3("pdgcodes.m");
546 while (in3.getline(str,200))
547 {
548 string tName;
549 int tCode;
550
551 std::stringstream tIS(str);
552 tIS >> tName >> tCode;
553 mDB->GetParticleType(tName)->SetPDGCode(tCode);
554 }
555
556}
557
558
559void Parser::ReadShare()
560{
561 char str[50];
562 char str1[200];
563
564 ParticleType *tPartBuf;
565 int tNum;
566
567 PRINT_DEBUG_1("Reading from |"<<(mInputDir+"/"+"particles.data").Data()<<"|");
568 ifstream in((mInputDir+"/"+"particles.data").Data());
569
570 int number=0;
571 if ((in) && (in.is_open()))
572 {
573 //START OF HEAD-LINE
574 in.ignore(200,'\n');
575 in.ignore(200,'\n');
576 in.ignore(200,'\n');
577 //END OF HEAD-LINE
578
579 while (in>>str)
580 {
581 if (/*(*str == '#')||*/(*str<65)||(*str>122))
582 {
583 in.getline(str1,200);
584 continue;
585 }
586 double mass, gamma, spin, I3, I, q, s, aq, as, c, ac, mc;
587
588 in>>mass>>gamma>>spin>>I>>I3>>q>>s>>aq>>as>>c>>ac>>mc;
589 number++;
590 PRINT_DEBUG_2(number<<" "<<str<<" "<<mass<<" "<<gamma<<" "<<spin<<" "<<I<<" "<<I3<<" "<<q<<" "<<aq<<" "<<s<<" "<<as<<" "<<c<<" "<<ac<<" "<<mc);
591
592 tPartBuf = new ParticleType();
593 tPartBuf->SetName(str);
594 tPartBuf->SetMass(mass);
595 tPartBuf->SetGamma(gamma);
596 tPartBuf->SetSpin(spin);
597 tPartBuf->SetBarionN((int) ((q+s+c)/3. - (aq+as+ac)/3.) );
598 tPartBuf->SetCharmN((int) (c - ac));
599 tPartBuf->SetStrangeness((int) (as-s));
600 tPartBuf->SetI(I);
601 tPartBuf->SetI3(I3);
602 tPartBuf->SetPDGCode((int) mc);
603 tPartBuf->SetNumber(number);
604 tNum = mDB->AddParticleType(tPartBuf);
d779ed42 605
606 delete tPartBuf;
2e967919 607 }
608 in.close();
609 }
610 else
611 {
612 PRINT_MESSAGE("File "<<(mInputDir+"/"+"particles.data").Data()<<" containing particle data not found!");
613 PRINT_MESSAGE("Please set the correct path to this file in the input parameter file");
614 PRINT_MESSAGE("Aborting!");
615 exit(0);
616 }
617
618 ifstream in2((mInputDir+"/decays.data").Data());
619 if ((in2) && (in2.is_open()))
620 {
621 //START OF HEAD-LINE
622 in2.ignore(200,'\n');
623 in2.ignore(200,'\n');
624 in2.ignore(200,'\n');
625 //END OF HEAD-LINE
626
f12dabf0 627 char tFather[50], tDaughter1[50], tDaughter2[50], tDaughter3[50];
2e967919 628 double tBRatio, tRatio;
629 int CGcoeff; // complete branching ratio by Clebsch-Gordan coefficient: 0-no 1-yes
630
631 while (in2>>str)
632 {
633 if (*str == '#')
634 {
635 in2.getline(str1,200);
636 continue;
637 }
638 in2>>tDaughter1>>tDaughter2>>tDaughter3;
639 if (!mDB->ExistsParticleType(tDaughter1)) {
640 PRINT_MESSAGE("Did not find the daughter 1 particle: " << tDaughter1);
641 PRINT_MESSAGE("Not adding channel");
642 in2.getline(str1,200);
643 continue;
644 }
645 if (!mDB->ExistsParticleType(tDaughter2)) {
646 PRINT_MESSAGE("Did not find the daughter 2 particle: " << tDaughter2);
647 PRINT_MESSAGE("Not adding channel");
648 in2.getline(str1,200);
649 continue;
650 }
651 if ((*tDaughter3>65)&&(*tDaughter3<122)&&(!mDB->ExistsParticleType(tDaughter3))) {
652 PRINT_MESSAGE("Did not find the daughter 3 particle: " << tDaughter3);
653 PRINT_MESSAGE("Not adding channel");
654 in2.getline(str1,200);
655 continue;
656 }
657
658 strcpy(tFather,str);
659 PRINT_DEBUG_2(tFather<<"\t"<<tDaughter1<<"\t"<<tDaughter2<<"\t");
660 if ((*tDaughter3>65)&&(*tDaughter3<122)) // check if first char is a letter - if yes then 3-body decay
661 {
662 in2>>tBRatio>>CGcoeff;
663 PRINT_DEBUG_2(tDaughter3<<" (3-body decay)\t");
664 if (mDB->ExistsParticleType(tFather)) {
665 mDB->GetParticleType(tFather)->SetDecayChannelCount3(mDB->GetParticleType(tFather)->GetDecayChannelCount3()+1);
666
667 tRatio=tBRatio;
668 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), mDB->GetParticleTypeIndex(tDaughter3));
669 if (mDB->GetParticleType(tDaughter1)->GetMass() +
670 mDB->GetParticleType(tDaughter2)->GetMass() +
671 mDB->GetParticleType(tDaughter3)->GetMass()
672 < mDB->GetParticleType(tFather)->GetMass())
673 (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
09b7c66c 674
675 delete newChannel;
2e967919 676
677 tRatio=tBRatio;
678 PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
679 }
680 else {
681 PRINT_MESSAGE("Did not find the father particle: " << tFather);
682 PRINT_MESSAGE("Not adding channel");
683 }
684 }
685 else // 2-body decay
686 {
687 tBRatio=atof(tDaughter3);
688
689 in2>>CGcoeff;
690 PRINT_DEBUG_2(" (2-body decay)\t");
691 if (mDB->ExistsParticleType(tFather)) {
692 mDB->GetParticleType(tFather)->SetDecayChannelCount2(mDB->GetParticleType(tFather)->GetDecayChannelCount2()+1);
693
694 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
695 if (CGcoeff) // complete branching ratio by Clebsch-Gordan coefficient
696 {
697 double j1, m1, j2, m2, J, M, CB;
698 J=mDB->GetParticleType(tFather)->GetI();
699 M=mDB->GetParticleType(tFather)->GetI3();
700 j1=mDB->GetParticleType(tDaughter1)->GetI();
701 m1=mDB->GetParticleType(tDaughter1)->GetI3();
702 j2=mDB->GetParticleType(tDaughter2)->GetI();
703 m2=mDB->GetParticleType(tDaughter2)->GetI3();
704 PRINT_DEBUG_2(" "<<J<<" "<<M<<" "<<j1<<" "<<m1<<" "<<j2<<" "<<m2);
705
706 CB = ClebschGordan(J, M, j1, m1, j2, m2);
707 tRatio = CB*CB * tBRatio;
708
709 // Multiply the Clebsh by two?
710 // The same spin, mass, strangeness
711 // and different I3?
712 if ((fabs(mDB->GetParticleType(tDaughter1)->GetSpin() - mDB->GetParticleType(tDaughter2)->GetSpin()) < 0.01) &&
713 (fabs(mDB->GetParticleType(tDaughter1)->GetMass() - mDB->GetParticleType(tDaughter2)->GetMass()) < 0.01) &&
714 (mDB->GetParticleType(tDaughter1)->GetStrangeness() == mDB->GetParticleType(tDaughter2)->GetStrangeness()) &&
715 (fabs(mDB->GetParticleType(tDaughter1)->GetI3() - mDB->GetParticleType(tDaughter2)->GetI3()) > 0.01))
716 {
717 PRINT_DEBUG_2("Multuplying Clebsch by two for " << tFather << "->" << tDaughter1 << "+" << tDaughter2);
718 tRatio *= 2.0;
719 }
720
721 PRINT_DEBUG_2(CB << '\t' << tBRatio << '\t' << tRatio<<"\t"<<CGcoeff);
722 }
723
724 else
725 {
726 tRatio=tBRatio;
727 PRINT_DEBUG_2(tBRatio << '\t' << tRatio);
728 }
729 DecayChannel *newChannel = new DecayChannel(tRatio, mDB->GetParticleTypeIndex(tDaughter1), mDB->GetParticleTypeIndex(tDaughter2), -1);
730 if (mDB->GetParticleType(tDaughter1)->GetMass() +
731 mDB->GetParticleType(tDaughter2)->GetMass()
732 < mDB->GetParticleType(tFather)->GetMass())
733 {
734 (mDB->GetParticleType(tFather))->AddDecayChannel(*newChannel);
735 PRINT_DEBUG_2("Added channel " << newChannel << " " << mDB->GetParticleTypeIndex(tFather) << " " << mDB->GetParticleTypeIndex(tDaughter1) << " " << mDB->GetParticleTypeIndex(tDaughter2));
6750cf94 736 delete newChannel;
2e967919 737 }
738 else
739 {
740
741 PRINT_DEBUG_2("Masses do not match! Not adding channel " << newChannel);
d58024dc 742
743 delete newChannel;
2e967919 744 }
745 }
746 else {
747 PRINT_MESSAGE("Did not find the father particle: " << tFather);
748 PRINT_MESSAGE("Not adding channel");
749 }
750 }
751 }
752 in2.close();
753 }
754 else {
755 PRINT_MESSAGE("File "<<(mInputDir+"/decays.data").Data()<<" with particle decay channels not found!");
756 PRINT_MESSAGE("No particle decays will be simulated");
757 }
758}
759
760double
761Parser::ClebschGordan(double aJot, double aEm, double aJot1, double aEm1, double aJot2, double aEm2)
762{
763 int mint, maxt;
764 double cgc = 0.0;
765 int titer;
766 double coef;
767
768 maxt = lrint(aJot1 + aJot2 - aJot);
769 mint = 0;
770 if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
771 if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
772 if (lrint(-(aJot-aJot2+aEm1)) > mint) mint = lrint(-(aJot-aJot2+aEm1));
773 if (lrint(-(aJot-aJot1-aEm2)) > mint) mint = lrint(-(aJot-aJot1-aEm2));
774
775 PRINT_DEBUG_3("mint " << mint << " " << aJot1 << " " << aEm1);
776 PRINT_DEBUG_3("maxt " << maxt << " " << aJot2 << " " << aEm2);
777
778 for (titer = mint; titer<=maxt; titer ++)
779 {
780 coef = TMath::Power(-1, titer);
781 PRINT_DEBUG_3("coef1 " << coef);
782 coef *= TMath::Sqrt((2*aJot+1)*
783 factorials[lrint(aJot1+aEm1)] *
784 factorials[lrint(aJot1-aEm1)] *
785 factorials[lrint(aJot2+aEm2)] *
786 factorials[lrint(aJot2-aEm2)] *
787 factorials[lrint(aJot+aEm)] *
788 factorials[lrint(aJot-aEm)]);
789 PRINT_DEBUG_3("coef2 " << coef);
790 coef /= (factorials[titer] *
791 factorials[lrint(aJot1+aJot2-aJot-titer)] *
792 factorials[lrint(aJot1-aEm1-titer)] *
793 factorials[lrint(aJot2+aEm2-titer)] *
794 factorials[lrint(aJot-aJot2+aEm1+titer)] *
795 factorials[lrint(aJot-aJot1-aEm2+titer)]);
796 PRINT_DEBUG_3("coef3 " << coef);
797
798 cgc += coef;
799 }
800
801 cgc *= DeltaJ(aJot1, aJot2, aJot);
802
803 return cgc;
804}
805
806double
807Parser::DeltaJ(double aJot1, double aJot2, double aJot)
808{
809 double res = TMath::Sqrt(1.0 *
810 factorials[lrint(aJot1+aJot2-aJot)] *
811 factorials[lrint(aJot1-aJot2+aJot)] *
812 factorials[lrint(-aJot1+aJot2+aJot)] /
813 factorials[lrint(aJot1+aJot2+aJot+1)]);
814
815 return res;
816}
817
818void
819Parser::ReadParameters()
820{
821 // Read the input directory
822 try {
823 mInputDir = sRPInstance->getPar("InputDirSHARE");
824 }
825 catch (STR tError) {
826 PRINT_DEBUG_1("Parser::ReadParameters - Caught exception " << tError);
827 PRINT_MESSAGE("Did not find SHARE input file location.");
828 PRINT_MESSAGE("Using default: '../share'");
829 mInputDir = "../share";
830 }
831}