MpdDecayerPyt8.cxx 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579
  1. // -------------------------------------------------------------------------
  2. // ----- MpdDecayerPyt8 header file -----
  3. // ----- Created 23/07/2020 -----
  4. // ----- A. Zinchenko -----
  5. // ----- External decayer for MPD using Pythia8 -----
  6. // ----- (adapted from TPythia8Decayer) -----
  7. // -------------------------------------------------------------------------
  8. #include "MpdDecayerPyt8.h"
  9. #include <Pythia8/Pythia.h>
  10. #include <TClonesArray.h>
  11. #include <TDatabasePDG.h>
  12. #include <TF1.h>
  13. #include <TLorentzVector.h>
  14. #include <TObjString.h>
  15. #include <TParticle.h>
  16. #include <TPythia8.h>
  17. #include <TString.h>
  18. #include <TSystem.h>
  19. #include <TVirtualMC.h>
  20. using namespace Pythia8;
  21. class MyRandomClass : public RndmEngine
  22. {
  23. double flat() { return gRandom->Rndm(); }
  24. };
  25. ClassImp(MpdDecayerPyt8);
  26. MpdDecayerPyt8* MpdDecayerPyt8::fgInstance = 0;
  27. ////////////////////////////////////////////////////////////////////////////////
  28. /// Get the singleton object.
  29. MpdDecayerPyt8* MpdDecayerPyt8::Instance()
  30. {
  31. if (!fgInstance) fgInstance = new MpdDecayerPyt8;
  32. return fgInstance;
  33. }
  34. ////////////////////////////////////////////////////////////////////////////////
  35. ///constructor
  36. MpdDecayerPyt8::MpdDecayerPyt8():
  37. fPythia8(new TPythia8()),
  38. fDebug(0),
  39. fGlobalPolar(0),
  40. fLambda(NULL),
  41. fRandom(NULL)
  42. {
  43. fPythia8->Pythia8()->readString("SoftQCD:elastic = on");
  44. // Disable cascade decays in one go
  45. fPythia8->Pythia8()->readString("ParticleDecays:limitRadius = on");
  46. fPythia8->Pythia8()->readString("ParticleDecays:rMax = 0.");
  47. fPythia8->Pythia8()->particleData.readString("59:m0 = 100.1396"); //??? - to avoid crash with Geant4
  48. fPythia8->Pythia8()->particleData.list(223);
  49. // Modify branching ratios
  50. ChangeBranchings();
  51. // Random number generator
  52. RndmEngine *rndm = new MyRandomClass();
  53. fPythia8->Pythia8()->setRndmEnginePtr(rndm);
  54. fPythia8->Pythia8()->init();
  55. Init();
  56. }
  57. ////////////////////////////////////////////////////////////////////////////////
  58. /// Initialize the decayer
  59. void MpdDecayerPyt8::Init()
  60. {
  61. static Bool_t init = kFALSE;
  62. if (init) return;
  63. init = kTRUE;
  64. //ForceDecay();
  65. fParticles = new TClonesArray("TParticle");
  66. // Lambda branching to p+\pi-
  67. fPythia8->Pythia8()->particleData.list(3122);
  68. //ParticleDataEntry* part = fPythia8->Pythia8()->particleData.particleDataEntryPtr(3122); // Lambda
  69. fLambda = fPythia8->Pythia8()->particleData.particleDataEntryPtr(3122); // Lambda
  70. DecayChannel& channel = fLambda->channel(0);
  71. fBranch = channel.bRatio();
  72. // Random number generator
  73. fRandom = new TRandom(0); // time-dependent seed
  74. //fRandom = gRandom;
  75. //fPythia8->Pythia8()->particleData.isResonance(113,kTRUE);
  76. fPythia8->Pythia8()->particleData.list(113);
  77. }
  78. ////////////////////////////////////////////////////////////////////////////////
  79. /// Decay a single particle
  80. void MpdDecayerPyt8::Decay(Int_t pdg, TLorentzVector* p)
  81. {
  82. //{ cout << " !!! Out !!! " << endl; exit(0); }
  83. // Reset internal particle container
  84. fParticles->Delete();
  85. fSourceFlag = kPythia;
  86. fMother.SetXYZT(0,0,0,0);
  87. if (!p) return;
  88. TRandom *saveRandom = gRandom;
  89. gRandom = fRandom;
  90. TParticle *part = gMC->GetStack()->GetCurrentTrack();
  91. if (fMothersPdg.find(pdg) == fMothersPdg.end()) {
  92. // Particle is not defined
  93. new ((*fParticles)[0]) TParticle(*part); // store mother particle
  94. fSourceFlag = kCustom;
  95. return;
  96. }
  97. TVector3 polar;
  98. part->GetPolarisation(polar);
  99. Bool_t polarFlag = kFALSE;
  100. if (TMath::Abs(polar.X()) > 0.0001 || TMath::Abs(polar.Z()) > 0.0001) {
  101. // Polarized lambda - simulate anysotropic decay only to p + \pi-
  102. polarFlag = kTRUE;
  103. if (gRandom->Rndm() < fBranch) {
  104. // Anysotropic decay
  105. //cout << " ----------------- " << fBranch << endl;
  106. new ((*fParticles)[0]) TParticle(*part); // store mother particle
  107. Gdecay (pdg, p);
  108. gRandom = saveRandom;
  109. return;
  110. } else {
  111. //cout << " ++++++++++ " << fBranch << endl;
  112. // Force the other channels
  113. DecayChannel& channel = fLambda->channel(0);
  114. channel.bRatio(0.0);
  115. }
  116. }
  117. ClearEvent();
  118. AppendParticle(pdg, p);
  119. Int_t idPart = fPythia8->Pythia8()->event[0].id();
  120. fPythia8->Pythia8()->particleData.mayDecay(idPart,kTRUE);
  121. fPythia8->Pythia8()->moreDecays();
  122. if (fDebug > 0) fPythia8->EventListing();
  123. if (polarFlag) {
  124. DecayChannel& channel = fLambda->channel(0);
  125. channel.bRatio(fBranch);
  126. }
  127. gRandom = saveRandom;
  128. }
  129. ////////////////////////////////////////////////////////////////////////////////
  130. ///import the decay products into particles array
  131. Int_t MpdDecayerPyt8::ImportParticles(TClonesArray *particles)
  132. {
  133. Int_t npart = 0;
  134. if (fSourceFlag == kPythia) {
  135. fParticles->Delete();
  136. npart = fPythia8->ImportParticles(fParticles, "All");
  137. TLorentzVector lvDaugh;
  138. for (Int_t j = 0; j < npart; ++j) {
  139. TParticle *part = (TParticle*) fParticles->UncheckedAt(j);
  140. part->ProductionVertex(lvDaugh);
  141. lvDaugh += fMother;
  142. part->SetProductionVertex(lvDaugh);
  143. new ((*particles)[j]) TParticle(*part);
  144. }
  145. } else {
  146. npart = fParticles->GetEntriesFast();
  147. for (Int_t j = 0; j < npart; ++j) new ((*particles)[j]) TParticle(*((TParticle*)fParticles->UncheckedAt(j)));
  148. }
  149. //return (fPythia8->ImportParticles(particles, "All"));
  150. return npart;
  151. }
  152. ////////////////////////////////////////////////////////////////////////////////
  153. /// Set forced decay mode
  154. void MpdDecayerPyt8::SetForceDecay(Int_t /*type*/)
  155. {
  156. printf("SetForceDecay not yet implemented !\n");
  157. }
  158. ////////////////////////////////////////////////////////////////////////////////
  159. /// ForceDecay not yet implemented
  160. void MpdDecayerPyt8::ForceDecay()
  161. {
  162. //printf("ForceDecay not yet implemented !\n");
  163. }
  164. ////////////////////////////////////////////////////////////////////////////////
  165. Float_t MpdDecayerPyt8::GetPartialBranchingRatio(Int_t /*ipart*/)
  166. {
  167. return 0.0;
  168. }
  169. ////////////////////////////////////////////////////////////////////////////////
  170. ///return lifetime in seconds of teh particle with PDG number pdg
  171. Float_t MpdDecayerPyt8::GetLifetime(Int_t pdg)
  172. {
  173. return (fPythia8->Pythia8()->particleData.tau0(pdg) * 3.3333e-12) ;
  174. }
  175. ////////////////////////////////////////////////////////////////////////////////
  176. ///to read a decay table (not yet implemented)
  177. void MpdDecayerPyt8::ReadDecayTable()
  178. {
  179. }
  180. ////////////////////////////////////////////////////////////////////////////////
  181. /// Append a particle to the stack
  182. void MpdDecayerPyt8::AppendParticle(Int_t pdg, TLorentzVector* p)
  183. {
  184. fPythia8->Pythia8()->event.append(pdg, 11, 0, 0, p->Px(), p->Py(), p->Pz(), p->E(), p->M());
  185. }
  186. ////////////////////////////////////////////////////////////////////////////////
  187. /// Clear the event stack
  188. void MpdDecayerPyt8::ClearEvent()
  189. {
  190. fPythia8->Pythia8()->event.clear();
  191. }
  192. ////////////////////////////////////////////////////////////////////////////////
  193. /// Simulate two-body decay process of Lambda-hyperon to p+\pi-
  194. void MpdDecayerPyt8::Gdecay(Int_t idpart, TLorentzVector* p)
  195. {
  196. Double_t pcm[2][4] = {{0},{0}};
  197. Double_t xm0 = p->M(), xm1 = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
  198. Double_t xm2 = TDatabasePDG::Instance()->GetParticle(-211)->Mass();
  199. Gdeca2 (xm0, xm1, xm2, pcm);
  200. // First, second decay products.
  201. TLorentzVector daughter1, daughter2;
  202. daughter1.SetPxPyPzE (pcm[0][0], pcm[0][1], pcm[0][2], pcm[0][3]);
  203. daughter2.SetPxPyPzE (pcm[1][0], pcm[1][1], pcm[1][2], pcm[1][3]);
  204. //Boost to lab frame
  205. //TVector3 boost = p->BoostVector(); //AZ - this is wrong !!!!
  206. TVector3 boost (0, 0, p->P()/p->E());
  207. daughter1.Boost(boost);
  208. daughter2.Boost(boost);
  209. daughter1.Transform(fRotation);
  210. daughter2.Transform(fRotation);
  211. TLorentzVector pos;
  212. gMC->TrackPosition(pos);
  213. //pos.Print();
  214. Int_t npart = fParticles->GetEntriesFast();
  215. new ((*fParticles)[npart++]) TParticle(2212, 1, 1, -1, 0, 0, daughter1, pos);
  216. new ((*fParticles)[npart]) TParticle(-211, 1, 1, -1, 0, 0, daughter2, pos);
  217. fSourceFlag = kCustom;
  218. // Modify mother particle parameters
  219. TParticle *mother = (TParticle*) fParticles->UncheckedAt(0);
  220. mother->SetFirstMother(1);
  221. mother->SetLastMother(-1);
  222. mother->SetFirstDaughter(2);
  223. mother->SetLastDaughter(3);
  224. mother->SetStatusCode(11);
  225. }
  226. ////////////////////////////////////////////////////////////////////////////////
  227. /// Simulate two-body decay process with anisotropic angular distribution in CMS.
  228. void MpdDecayerPyt8::Gdeca2(Double_t xm0, Double_t xm1, Double_t xm2, Double_t pcm[2][4])
  229. {
  230. Double_t random[2], pvert[3], costh = 0.0, sinth = 0.0, phi = 0.0;
  231. Double_t e1 = (xm0 * xm0 + xm1 * xm1 - xm2 * xm2) / (2. * xm0);
  232. Double_t p1 = TMath::Sqrt (TMath::Abs ((e1 - xm1) * (e1 + xm1)));
  233. static TRandom *myRandom = NULL;
  234. if (!myRandom) { myRandom = new TRandom(); myRandom->SetSeed(0); }
  235. TRandom *oldRandom = gRandom;
  236. gRandom = myRandom;
  237. gRandom->RndmArray (2, random);
  238. // Sanity check - should not happen (legacy code)
  239. TParticle *part = gMC->GetStack()->GetCurrentTrack();
  240. if (part->GetPdgCode() != 3122) { cout << " ??? Not Lambda - exit " << part->GetPdgCode() << endl; exit(0); }
  241. TLorentzVector lorV;
  242. part->Momentum(lorV);
  243. lorV.Vect().GetXYZ(pvert);
  244. TVector3 polar;
  245. part->GetPolarisation(polar);
  246. if (fGlobalPolar) polar.SetMag(part->GetWeight()); // particle polarization value is passed as its weight
  247. if (TMath::Abs(polar.X()) < 0.0001 && TMath::Abs(polar.Z()) < 0.0001 && TMath::Abs(polar.Y()) < 0.0001) {
  248. costh = 2. * random[0] - 1.0;
  249. phi = TMath::TwoPi() * random[1];
  250. } else {
  251. // Anisotropic decay angular distribution (0.5*(1+alpha*P*cos)).
  252. //Anisotropy (pvert, random, polar.X(), phi, costh);
  253. Anisotropy (pvert, random, polar, phi, costh);
  254. }
  255. if (TMath::Abs(costh) >= 1.0) costh = TMath::Sign (1.0, costh);
  256. else sinth = TMath::Sqrt ((1. - costh) * (1. + costh));
  257. // Polar co-ordinates to momentum components.
  258. pcm[0][0] = p1 * sinth * TMath::Cos(phi);
  259. pcm[0][1] = p1 * sinth * TMath::Sin(phi);
  260. pcm[0][2] = p1 * costh;
  261. pcm[0][3] = e1;
  262. // Generate second decay product.
  263. pcm[1][0] = -pcm[0][0];
  264. pcm[1][1] = -pcm[0][1];
  265. pcm[1][2] = -pcm[0][2];
  266. pcm[1][3] = TMath::Sqrt (pcm[1][0] * pcm[1][0] + pcm[1][1] * pcm[1][1] + pcm[1][2] * pcm[1][2] + xm2 * xm2);
  267. gRandom = oldRandom;
  268. }
  269. ////////////////////////////////////////////////////////////////////////////////
  270. /// Simulate anisotropic (due to polarization) decay of lambda
  271. //void MpdDecayerPyt8::Anisotropy (Double_t *pvert, Double_t *rndm, Double_t polar, Double_t &phi, Double_t &costh)
  272. void MpdDecayerPyt8::Anisotropy (Double_t *pvert, Double_t *rndm, TVector3& polar3, Double_t &phi, Double_t &costh)
  273. {
  274. // Simulate anisotropic (due to polarization) decay of lambda
  275. //exit(0);
  276. //std::ofstream outfile ("costh.txt");
  277. //freopen ("costh.txt", "w", stdout);
  278. Double_t polar = 0.0;
  279. if (fGlobalPolar) polar = polar3.Mag(); // global polarization - CHECK THIS
  280. else polar = polar3.X(); // transverse polar.
  281. //const Double_t alpha = 0.642;
  282. //static TF1 f("f","1+0.642*[0]*x",-1,1);
  283. static TF1 f("f","1+0.732*[0]*x",-1,1); // AZ 10.06.21 - new value
  284. f.SetParameter(0,polar);
  285. f.SetNpx(1000);
  286. Double_t costhe = f.GetRandom();
  287. Double_t sinth = 0.0;
  288. if (TMath::Abs(costhe) >= 1.0) costhe = TMath::Sign (1.0,costhe);
  289. else sinth = TMath::Sqrt (1.0 - costhe * costhe);
  290. phi = TMath::TwoPi() * rndm[1];
  291. //std::cout << " Cos: " << costhe << " " << phi << " " << pvert[0] << " " << pvert[1] << " " << pvert[2] << std::endl;
  292. // Compute normal vector
  293. TVector3 beam(0.0, 0.0, 1.0);
  294. TVector3 lambda(pvert[0], pvert[1], pvert[2]);
  295. TVector3 norm;
  296. if (fGlobalPolar) norm = polar3.Unit(); // global polarization
  297. else norm = beam.Cross(lambda).Unit(); // transverse polarization
  298. // Unit vector with theta and phi w.r.t. normal
  299. TVector3 unit(sinth*TMath::Cos(phi), sinth*TMath::Sin(phi), costhe);
  300. // Coordinate system transformation
  301. // (from lambda to lab)
  302. TVector3 lambU = lambda.Unit(), lambZ(0,0,1), lambX(1,0,0), lambY(0,1,0);
  303. lambZ.RotateUz(lambU);
  304. lambX.RotateUz(lambU);
  305. lambY.RotateUz(lambU);
  306. TRotation rotL;
  307. rotL.RotateAxes(lambX,lambY,lambZ); // transformation to lab. system
  308. fRotation = rotL;
  309. rotL.Invert(); // from lab. to lambda
  310. unit.RotateUz(norm); // to lab. system
  311. unit.Transform(rotL); // to lambda system
  312. /*
  313. TVector3 normZ(0,0,1), normX(1,0,0), normY(0,1,0);
  314. normZ.RotateUz(norm);
  315. normX.RotateUz(norm);
  316. normY.RotateUz(norm);
  317. TRotation rotN;
  318. rotN.RotateAxes(normX,normY,normZ); // transformation to lab. system
  319. unit.Transform(rotN); // to lab. system
  320. unit.Transform(rotL); // to lambda system
  321. */
  322. costh = TMath::Cos(unit.Theta());
  323. phi = unit.Phi();
  324. //std::cout << costh << " " << phi << std::endl;
  325. //outfile << i << " " << costh << endl;
  326. //outfile.close();
  327. }
  328. ///////////////////////////////////////////////////////////////////////////////
  329. /// Add PDG of particle to be decayed by this package
  330. void MpdDecayerPyt8::AddMotherPdg(Int_t pdg)
  331. {
  332. if (fMothersPdg.find(pdg) != fMothersPdg.end()) return;
  333. fMothersPdg.insert(pdg);
  334. gMC->SetUserDecay(pdg);
  335. }
  336. ////////////////////////////////////////////////////////////////////////////////
  337. /// Decay a TParticle
  338. void MpdDecayerPyt8::Decay(TParticle* part)
  339. {
  340. // Decay function
  341. //{ cout << " !!! Out !!! " << endl; exit(0); }
  342. // Reset internal particle container
  343. fParticles->Delete();
  344. fSourceFlag = kPythia;
  345. /*
  346. if (!p) return;
  347. p->ProductionVertex(fMother);
  348. TPythia6 *pythia = TPythia6::Instance();
  349. pythia->Py1ent(0, p->GetPdgCode(), p->Energy(), p->Theta(), p->Phi());
  350. pythia->GetPrimaries();
  351. */
  352. if (!part) return;
  353. //part->Print();
  354. TRandom *saveRandom = gRandom;
  355. gRandom = fRandom;
  356. part->ProductionVertex(fMother);
  357. Int_t pdg = part->GetPdgCode();
  358. TLorentzVector p;
  359. part->Momentum(p);
  360. if (part->GetFirstMother() < 0) {
  361. // Primary particle - select mass
  362. Double_t mass = fPythia8->Pythia8()->particleData.mSel(pdg);
  363. //cout << "mmm " << mass << " " << part->GetFirstMother() << endl;
  364. p.SetE (TMath::Sqrt (p.P()*p.P() + mass*mass));
  365. }
  366. //p.Print();
  367. ClearEvent();
  368. AppendParticle(pdg, &p);
  369. Int_t idPart = fPythia8->Pythia8()->event[0].id();
  370. fPythia8->Pythia8()->particleData.mayDecay(idPart,kTRUE);
  371. fPythia8->Pythia8()->moreDecays();
  372. if (fDebug > 0) fPythia8->EventListing();
  373. gRandom = saveRandom;
  374. }
  375. //__________________________________________________________________________
  376. void MpdDecayerPyt8::ChangeBranchings()
  377. {
  378. // Change branching ratios of some decay channels
  379. // (using description file)
  380. TString path = gSystem->ExpandPathName("$VMCWORKDIR");
  381. path += "/gconfig/MpdDecayConfig.txt";
  382. ifstream fin(path.Data());
  383. TString chline;
  384. TObjArray *tokens = NULL;
  385. while(1) {
  386. chline.ReadLine(fin);
  387. if (fin.eof()) break;
  388. cout << chline << endl;
  389. if (chline.Contains("#")) continue; // comment
  390. // Tokenize line
  391. tokens = chline.Tokenize(";"); // ";" - channel terminator / separator
  392. ChangeParticleBr(tokens);
  393. if (tokens) { tokens->Delete(); delete tokens; tokens = NULL; }
  394. }
  395. }
  396. //__________________________________________________________________________
  397. void MpdDecayerPyt8::ChangeParticleBr(TObjArray *channels)
  398. {
  399. // Change branching ratios of selected channels for single particle
  400. Int_t nchan = channels->GetEntriesFast(); // number of channels to be changed
  401. Int_t pdg = ((TObjString*)channels->UncheckedAt(0))->String().Atoi();
  402. //cout << pdg << endl;
  403. ParticleDataEntryPtr part = fPythia8->Pythia8()->particleData.particleDataEntryPtr(pdg); // mother particle
  404. Int_t ndecays = part->sizeChannels();
  405. TObjArray *tokens = NULL;
  406. Double_t changedBrs = 0.0, unchangedBrs = 1.0;
  407. vector<Int_t> decayInds;
  408. vector<Double_t> bRatios;
  409. for (Int_t ich = 0; ich < nchan; ++ich) {
  410. TString chline = ((TObjString*)channels->UncheckedAt(ich))->String();
  411. cout << chline << endl;
  412. // Tokenize channel
  413. tokens = chline.Tokenize(":");
  414. // Get daughters' PDGs
  415. Int_t nprongs = tokens->GetEntriesFast() - 2;
  416. vector<Int_t> pdgs;
  417. for (Int_t id = 0; id < nprongs; ++id) {
  418. pdgs.push_back(((TObjString*)tokens->UncheckedAt(1+id))->String().Atoi());
  419. }
  420. TString scale = ((TObjString*)tokens->Last())->String();
  421. scale.Remove(0,1);
  422. Double_t scaleBr = scale.Atof();
  423. // Pick channel requested
  424. for (Int_t idecay = 0; idecay < ndecays; ++idecay) {
  425. DecayChannel& channel = part->channel(idecay);
  426. if (channel.multiplicity() != nprongs) continue;
  427. // Check daughters
  428. Int_t ok = 0;
  429. for (Int_t id = 0; id < nprongs; ++id) {
  430. if (channel.product(id) != pdgs[id]) break;
  431. ++ok;
  432. }
  433. if (ok != nprongs) continue;
  434. changedBrs += channel.bRatio() * scaleBr;
  435. unchangedBrs -= channel.bRatio();
  436. decayInds.push_back(idecay);
  437. bRatios.push_back (channel.bRatio() * scaleBr);
  438. //cout << channel.bRatio() << endl;
  439. }
  440. //cout << changedBrs << endl;
  441. if (tokens) { tokens->Delete(); delete tokens; tokens = NULL; }
  442. } // for (Int_t ich = 0; ich < nchan;
  443. part->rescaleBR ((1.0 - changedBrs) / unchangedBrs);
  444. for (Int_t ich = 0; ich < nchan; ++ich) {
  445. DecayChannel& channel = part->channel(decayInds[ich]);
  446. channel.bRatio(bRatios[ich]);
  447. }
  448. // Check
  449. Double_t totalBr = 0;
  450. for (Int_t idecay = 0; idecay < ndecays; ++idecay) {
  451. DecayChannel& channel = part->channel(idecay);
  452. totalBr += channel.bRatio();
  453. }
  454. cout << " Total branching after rescaling: " << totalBr << endl;
  455. }