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