diff --git a/Cleopatra/InFileCreator.h b/Cleopatra/InFileCreator.h index ce64422..b6cf622 100644 --- a/Cleopatra/InFileCreator.h +++ b/Cleopatra/InFileCreator.h @@ -47,6 +47,7 @@ int GetLValue(string spdf){ return -1; } + int InFileCreator(string readFile, string infile, double angMin, double angMax, double angStep) { //================= read infile. extract the reactions, write pptolemy infile for each reaction @@ -89,28 +90,37 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, str2[0] = "(" + str2[0]; - //TODO use mass table for d, p, t, 3He int lenStr20 = str2[0].length(); size_t posTemp1 = str2[0].find(","); - string ma = str2[0].substr(1, posTemp1-1); + string mass_a = str2[0].substr(1, posTemp1-1); size_t posTemp2 = str2[0].find(")"); - string mb = str2[0].substr(posTemp1+1, posTemp2-posTemp1-1); - ///printf(" ma : |%s| , mb : |%s| \n", ma.c_str(), mb.c_str()); - Isotope isotopea(ma); - Isotope isotopeb(mb); - + string mass_b = str2[0].substr(posTemp1+1, posTemp2-posTemp1-1); + ///printf(" mass_a : |%s| , mass_b : |%s| \n", mass_a.c_str(), mass_b.c_str()); + Isotope iso_a(mass_a); + Isotope iso_b(mass_b); + + // Check is the Reaction supported bool isReactionSupported = false; - if( isotopea.A <= 4 && isotopea.Z <= 2 && isotopeb.A <=4 && isotopeb.Z <=2 ) isReactionSupported = true; - - ///if( ma == "p" || ma == "d" || ma == "t" || ma == "3He" || mb == "n"){ - /// if( mb == "p" || mb == "d" || mb == "t" || mb == "3He" || mb == "n" ) isReactionSupported = true; - ///} - + bool isTransferReaction = true; + + if( iso_a.A <= 4 && iso_a.Z <= 2 && iso_b.A <=4 && iso_b.Z <=2 ) isReactionSupported = true; + + ///======= elastics-ish scattering + if( iso_a.Mass == iso_b.Mass ) isTransferReaction = false; + + ///======= p/n-exchange is not supported + if( iso_a.A == iso_b.A && iso_a.Z != iso_b.Z ) isReactionSupported = false; + + ///======= 3-nucleons transfer is not supported. e.g. (n,a), (p,a), (a,n), (a,p) + int numNucleonsTransfer = iso_a.A - iso_b.A; + if( abs(numNucleonsTransfer) >= 3 ) isReactionSupported = false; + if( isReactionSupported == false ){ printf(" ===> Ignored. Reaction type not supported. \n"); continue; } + // Continues to decode the input string string gsSpinA = str0[1]; string orbital = str0[2]; @@ -127,28 +137,25 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, string isoB = str2[1]; string reactionType = str2[0]; - Isotope isotopeA(str1[0]); - Isotope isotopeB(str2[1]); + Isotope iso_A(str1[0]); + Isotope iso_B(str2[1]); - ///check reaction valid by balancing the A and Z number; - ///printf("A: %d + %d = %d + %d \n", isotopeA.A, isotopea.A, isotopeB.A, isotopeb.A); - ///printf("Z: %d + %d = %d + %d \n", isotopeA.Z, isotopea.Z, isotopeB.Z, isotopeb.Z); - - if( isotopeA.A + isotopea.A != isotopeB.A + isotopeb.A || isotopeA.Z + isotopea.Z != isotopeB.Z + isotopeb.Z ) { - printf(" ====> ERROR! A-number or Z-number not balanced. \n"); - - Isotope isotopeK(isotopeA.A + isotopea.A - isotopeb.A, isotopeA.Z + isotopea.Z - isotopeb.Z); - - printf(" try : %s ??\n", isotopeK.Name.c_str()); - + /// check is iso_A or iso_B exist in the mass table + if( iso_A.Mass == -404 || iso_B.Mass == -404 ){ + printf(" ===> Error! mass does not found. \n"); continue; } - bool isTransferReaction = true; - if( ma == mb ) isTransferReaction = false; + /// check reaction valid by balancing the A and Z number; + if( iso_A.A + iso_a.A != iso_B.A + iso_b.A || iso_A.Z + iso_a.Z != iso_B.Z + iso_b.Z ) { + printf("====> ERROR! A-number or Z-number not balanced. \n"); + Isotope isotopeK(iso_A.A + iso_a.A - iso_b.A, iso_A.Z + iso_a.Z - iso_b.Z); + printf(" try : %s(%s,%s)%s ??\n", iso_A.Name.c_str(), iso_a.Name.c_str(), iso_b.Name.c_str(), isotopeK.Name.c_str()); + continue; + } if( isTransferReaction && potential.length() != 2 ){ - printf(" ===> ERROR! Potential input should be 2 charaters! skipped. \n"); + printf("====> ERROR! Potential input should be 2 charaters! skipped. \n"); continue; } @@ -158,13 +165,11 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, int spdf = -1; if( isTransferReaction ) { - - ///printf("------------ %d nucleon(s) transfer \n", abs(isotopea.A - isotopeb.A)); - + ///printf("------------ %d nucleon(s) transfer \n", abs(iso_a.A - iso_b.A)); node = orbital.substr(0,1); // single nucleon transfer - if( abs(isotopea.A - isotopeb.A) == 1 ){ + if( abs(iso_a.A - iso_b.A) == 1 ){ lValue = orbital.substr(1,1); jValue = orbital.substr(2); ///printf(" l : %s, j : %s \n", lValue.c_str(), jValue.c_str()); @@ -172,13 +177,13 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, } // two-nucleons transfer - if( abs(isotopea.A - isotopeb.A) == 2 ){ + if( abs(iso_a.A - iso_b.A) == 2 ){ size_t posEq = orbital.find('='); lValue = orbital.substr(posEq+1,1); spdf=atoi(lValue.c_str()); } - if( abs(isotopea.A - isotopeb.A) == 0 ){ + if( abs(iso_a.A - iso_b.A) == 0 ){ printf(" ===? skipped. p-n exchange reaction does not support. \n"); } @@ -186,7 +191,6 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, printf(" ===> skipped. Not reconginzed orbital-label. (user input : l=%s | %s) \n", lValue.c_str(), orbital.c_str()); continue; } - } //get Beam energy, distingusih MeV or MeV/u @@ -198,9 +202,8 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, } } string unit = reactionEnergy.substr(pos+1); - //string beamParticle = reactionType.substr(1,1); int factor = 1; - if( unit == "MeV/u") factor = isotopea.A; + if( unit == "MeV/u") factor = iso_a.A; double totalBeamEnergy = atof(reactionEnergy.substr(0, pos+1).c_str()) * factor; ///printf("unit : |%s| , %f\n", unit.c_str(), totalBeamEnergy); @@ -212,59 +215,122 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, ///printf(" orbital : %s \n", orbital.c_str()); ///printf(" Ex [MeV] : %s \n", Ex.c_str()); - double Qvalue = isotopea.Mass + isotopeA.Mass - isotopeb.Mass - isotopeB.Mass; + double Qvalue = iso_a.Mass + iso_A.Mass - iso_b.Mass - iso_B.Mass; ///printf("Q-Value = %f MeV\n", Qvalue); - ///printf("A: %d, Z: %d, mass: %f MeV/c2 \n", isotopeA.A, isotopeA.Z, isotopeA.Mass); - ///printf("A: %d, Z: %d, mass: %f MeV/c2 \n", isotopeB.A, isotopeB.Z, isotopeB.Mass); - if( isotopeA.Mass == -404 || isotopeB.Mass == -404 ){ - printf(" ===> Error! mass does not found. \n"); - continue; - } - - //write ptolmey infile + //########################################################## + //============ write ptolmey infile numOfReaction ++ ; + //================ elastic-ish transfer + if( isTransferReaction == false ){ + if ( atof(Ex.c_str()) == 0.0 ) { + fprintf(file_out, "$============================================ ELab=%5.2f(%s+%s)%s\n", + totalBeamEnergy, mass_a.c_str(), isoA.c_str(), potential.c_str()); + fprintf(file_out, "reset\n"); + fprintf(file_out, "CHANNEL %s + %s\n", mass_a.c_str(), isoA.c_str()); + fprintf(file_out, "r0target\n"); + fprintf(file_out, "ELAB = %f\n", totalBeamEnergy); + fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); + string pot1Name = potential.substr(0,1); + string pot1Ref = potentialRef(pot1Name); + fprintf(file_out, "$%s\n", pot1Ref.c_str()); + CallPotential(pot1Name, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.Z); + fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); + fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); + fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); + fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); + fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); + fprintf(file_out, "ELASTIC SCATTERING\n"); + fprintf(file_out, ";\n"); + }else{ + fprintf(file_out, "$============================================ Ex=%s(%s+%s|%s%s)%s,ELab=%5.2f\n", + Ex.c_str(), mass_a.c_str(), isoA.c_str(), spin.c_str(), parity.c_str(), potential.c_str(),totalBeamEnergy); + fprintf(file_out, "reset\n"); + fprintf(file_out, "REACTION: %s%s%s(%s%s %s) ELAB=%7.3f\n", + isoA.c_str(), reactionType.c_str(), isoB.c_str(), spin.c_str(), parity.c_str(), Ex.c_str(), totalBeamEnergy); + fprintf(file_out, "PARAMETERSET ineloca2 r0target\n"); + fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); + if( str0.size() >= 8 ){ + fprintf(file_out, "BETA=%s\n", str0[7].c_str()); //deformation length + } + string pot1Name = potential.substr(0,1); + string pot1Ref = potentialRef(pot1Name); + fprintf(file_out, "$%s\n", pot1Ref.c_str()); + CallPotential(pot1Name, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.Z); + fprintf(file_out, "INCOMING\n"); + fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); + fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); + fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0); + ///fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); + ///fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); + fprintf(file_out, ";\n"); + + fprintf(file_out, "OUTGOING\n"); + fprintf(file_out, "$%s\n", pot1Ref.c_str()); + CallPotential(pot1Name, iso_A.A, iso_A.Z, totalBeamEnergy - atof(Ex.c_str()), iso_a.Z); + fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); + fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); + fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0); + ///fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); + ///fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); + ///fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); + fprintf(file_out, ";\n"); + } + } + + //================ Transfer reaction if( isTransferReaction ){ fprintf(file_out, "$============================================ Ex=%s(%s)%s\n", Ex.c_str(), orbital.c_str(), potential.c_str()); fprintf(file_out, "reset\n"); fprintf(file_out, "REACTION: %s%s%s(%s%s %s) ELAB=%7.3f\n", isoA.c_str(), reactionType.c_str(), isoB.c_str(), spin.c_str(), parity.c_str(), Ex.c_str(), totalBeamEnergy); - if( isotopea.A <= 2 && isotopea.Z <= 1){ // incoming d or p - fprintf(file_out, "PARAMETERSET dpsb r0target \n"); - fprintf(file_out, "$lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n"); - fprintf(file_out, "\n"); - fprintf(file_out, "PROJECTILE \n"); - fprintf(file_out, "wavefunction av18 \n"); - fprintf(file_out, "r0=1 a=0.5 l=0 rc0=1.2\n"); - } - if( 3 <= isotopea.A && isotopea.A <= 4 && abs(isotopea.A-isotopeb.A) == 1){ - fprintf(file_out, "PARAMETERSET alpha3 r0target \n"); - fprintf(file_out, "$lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n"); - fprintf(file_out, "\n"); - fprintf(file_out, "PROJECTILE \n"); - fprintf(file_out, "wavefunction phiffer \n"); - if( isotopeb.A <= 3 ){ - fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.31 v=179.94 r=0.54 a=0.68 param1=0.64 param2=1.13 rc=2.0\n"); + + //-------- Projectile (the light particle) + if( abs(numNucleonsTransfer) == 1 ){ + if( iso_a.A <= 2 && iso_a.Z <= 1 && iso_b.A <=2 && iso_b.Z <= 1){ // incoming d or p + fprintf(file_out, "PARAMETERSET dpsb r0target \n"); + fprintf(file_out, "lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n"); + fprintf(file_out, "\n"); + fprintf(file_out, "PROJECTILE \n"); + fprintf(file_out, "wavefunction av18 \n"); + fprintf(file_out, "r0=1 a=0.5 l=0 rc0=1.2\n"); } - if( isotopeb.A == 4 ){ - fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.61 v=202.21 r=.93 a=.66 param1=.81 param2=.87 rc=2.0 $ rc=2 is a quirk\n"); - } - } - if( abs(isotopea.A-isotopeb.A) == 2 ){ // 2 neutron transfer + + if( (3 <= iso_a.A && iso_a.A <= 4) || (3 <= iso_b.A && iso_b.A <= 4) ){ + fprintf(file_out, "PARAMETERSET alpha3 r0target \n"); + fprintf(file_out, "lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n"); + fprintf(file_out, "\n"); + fprintf(file_out, "PROJECTILE \n"); + fprintf(file_out, "wavefunction phiffer \n"); + + if( iso_a.Z + iso_b.Z == 2){ // (t,d) or (d,t) + fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.30 v=172.88 r=0.56 a=0.69 param1=0.64 param2=1.15 rc=2.0\n"); + } + + if( iso_a.Z + iso_b.Z == 3){ // (3He,d) or (d, 3He) + fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.31 v=179.94 r=0.54 a=0.68 param1=0.64 param2=1.13 rc=2.0\n"); + } + + if( iso_b.A == 4 ){ + fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.61 v=202.21 r=.93 a=.66 param1=.81 param2=.87 rc=2.0 $ rc=2 is a quirk\n"); + } + } + }else if( abs(numNucleonsTransfer) == 2 ){ // 2 nucleons transfer fprintf(file_out, "PARAMETERSET alpha3 r0target\n"); fprintf(file_out, "lstep=1 lmin=0 lmax=30 maxlextrap=0 ASYMPTOPIA=40\n"); fprintf(file_out, "\n"); fprintf(file_out, "PROJECTILE\n"); + fprintf(file_out, "wavefunction phiffer\n"); fprintf(file_out, "L = 0 NODES=0 R0 = 1.25 A = .65 RC0 = 1.25\n"); - } + } fprintf(file_out, ";\n"); //===== TARGET fprintf(file_out, "TARGET\n"); ///check Ex is above particle threshold - double nThreshold = isotopeB.CalSp(0,1); - double pThreshold = isotopeB.CalSp(1,0); + double nThreshold = iso_B.CalSp(0,1); + double pThreshold = iso_B.CalSp(1,0); bool isAboveThreshold = false; double ExEnergy = atof(Ex.c_str()); if( ExEnergy > nThreshold || ExEnergy > pThreshold ) { @@ -272,7 +338,7 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, printf(" Ex = %.3f MeV is above thresholds; Sp = %.3f MeV, Sn = %.3f MeV\n", ExEnergy, pThreshold, nThreshold); } - if( abs(isotopea.A-isotopeb.A) == 1 ){ + if( abs(iso_a.A-iso_b.A) == 1 ){ fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); if( isAboveThreshold ) { fprintf(file_out, "nodes=%s l=%d jp=%s E=-.2 $node is n-1, set binding 200 keV \n", node.c_str(), spdf, jValue.c_str()); @@ -284,7 +350,7 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, fprintf(file_out, "rc0=1.3 \n"); } - if( abs(isotopea.A-isotopeb.A) == 2 ){ + if( abs(iso_a.A-iso_b.A) == 2 ){ fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); if( isAboveThreshold ){ fprintf(file_out, "nodes=%s L=%d E=-.2 $node is n-1, binding by 200 keV \n", node.c_str(), spdf); @@ -299,7 +365,7 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, string pot1Name = potential.substr(0,1); string pot1Ref = potentialRef(pot1Name); fprintf(file_out, "INCOMING $%s\n", pot1Ref.c_str()); - CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z); + CallPotential(pot1Name, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.Z); fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); @@ -312,7 +378,7 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, fprintf(file_out, "OUTGOING $%s\n", pot2Ref.c_str()); //printf(" total Beam Energy : %f | Qvalue : %f | Ex : %f \n", totalBeamEnergy, Qvalue, atof(Ex.c_str())); double eBeam = totalBeamEnergy + Qvalue - atof(Ex.c_str()); - CallPotential(pot2Name, isotopeB.A, isotopeB.Z, eBeam, isotopeb.Z); + CallPotential(pot2Name, iso_B.A, iso_B.Z, eBeam, iso_b.Z); fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); @@ -321,60 +387,7 @@ int InFileCreator(string readFile, string infile, double angMin, double angMax, fprintf(file_out, ";\n"); } - if( isTransferReaction == false ){ - if ( atof(Ex.c_str()) == 0.0 ) { - fprintf(file_out, "$============================================ ELab=%5.2f(%s+%s)%s\n", - totalBeamEnergy, ma.c_str(), isoA.c_str(), potential.c_str()); - fprintf(file_out, "reset\n"); - fprintf(file_out, "CHANNEL %s + %s\n", ma.c_str(), isoA.c_str()); - fprintf(file_out, "ELAB = %f\n", totalBeamEnergy); - fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); - string pot1Name = potential.substr(0,1); - string pot1Ref = potentialRef(pot1Name); - fprintf(file_out, "$%s\n", pot1Ref.c_str()); - CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z); - fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); - fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); - fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); - fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); - fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); - fprintf(file_out, "ELASTIC SCATTERING\n"); - fprintf(file_out, ";\n"); - }else{ - fprintf(file_out, "$============================================ Ex=%s(%s+%s|%s%s)%s,ELab=%5.2f\n", - Ex.c_str(), ma.c_str(), isoA.c_str(), spin.c_str(), parity.c_str(), potential.c_str(),totalBeamEnergy); - fprintf(file_out, "reset\n"); - fprintf(file_out, "REACTION: %s%s%s(%s%s %s) ELAB=%7.3f\n", - isoA.c_str(), reactionType.c_str(), isoB.c_str(), spin.c_str(), parity.c_str(), Ex.c_str(), totalBeamEnergy); - fprintf(file_out, "PARAMETERSET ineloca2 r0target\n"); - fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str()); - if( str0.size() >= 8 ){ - fprintf(file_out, "BETA=%s\n", str0[7].c_str()); //deformation length - } - string pot1Name = potential.substr(0,1); - string pot1Ref = potentialRef(pot1Name); - fprintf(file_out, "$%s\n", pot1Ref.c_str()); - CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z); - fprintf(file_out, "INCOMING\n"); - fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); - fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); - fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0); - //fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); - //fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); - fprintf(file_out, ";\n"); - - fprintf(file_out, "OUTGOING\n"); - fprintf(file_out, "$%s\n", pot1Ref.c_str()); - CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy - atof(Ex.c_str()), isotopea.Z); - fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a); - fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai); - fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0); - //fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi); - //fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso); - //fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0); - fprintf(file_out, ";\n"); - } - } + fprintf(file_out, "anglemin=%f anglemax=%f anglestep=%f\n", angMin, angMax, angStep); diff --git a/Cleopatra/Isotope.h b/Cleopatra/Isotope.h index c00353c..afabaa9 100644 --- a/Cleopatra/Isotope.h +++ b/Cleopatra/Isotope.h @@ -1,534 +1,535 @@ -/*********************************************************************** - * - * This is Isotope.h, To extract the isotope mass from massXX.txt - * - *------------------------------------------------------- - * created by Ryan (Tsz Leung) Tang, Nov-18, 2018 - * email: goluckyryan@gmail.com - * ********************************************************************/ - - -#ifndef ISOTOPE_H -#define ISOTOPE_H - -#include -#include -#include -#include -#include -#include -#include "constant.h" // amu -#include //atoi -#include -using namespace std; - -string massData="/Cleopatra/mass20.txt"; - -// about the mass**.txt -// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu -// mass excess uncertaintly -// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons -// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons - -class Isotope { -public: - int A, Z; - double Mass, MassError, BEA; - string Name, Symbol; - string dataSource; - - Isotope(){findHeliosPath();}; - Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; - Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; - - void SetIso(int a, int z); - void SetIsoByName(string name); - - double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal - double CalSp2(int a, int z); // this is for (a,z) nucleus removal - - double CalBeta(double T){ - double Etot = Mass + T; - double gamma = 1 + T/Mass; - double beta = sqrt(1 - 1 / gamma / gamma ) ; - return beta; - } - - void Print(); - void ListShell(); - -private: - void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; - void FindMassByName(string name); // give Z, mass, massError, BEA; - - int TwoJ(int nShell); - string Orbital(int nShell); - int magic(int i){ - switch (i){ - case 0: return 2; break; - case 1: return 8; break; - case 2: return 20; break; - case 3: return 28; break; - case 4: return 40; break; - case 5: return 50; break; - case 6: return 82; break; - case 7: return 128; break; - } - return 0; - } - - int magicShellID(int i){ - switch (i){ - case 0: return 0; break; - case 1: return 2; break; - case 2: return 5; break; - case 3: return 6; break; - case 4: return 9; break; - case 5: return 10; break; - case 6: return 15; break; - case 7: return 21; break; - } - return 0; - } - - int fileStartLine; - int fileEndLine; - int lineMass050_099; - int lineMass100_149; - int lineMass150_199; - int lineMass200; - - - void setFileLines(){ - fileStartLine = 37; - fileEndLine = 3594; - - lineMass050_099 = 466; - lineMass100_149 = 1160; - lineMass150_199 = 1994; - lineMass200 = 2774; - } - - char * heliosPath; - bool isFindOnce; - - void findHeliosPath(){ - dataSource = massData; - heliosPath = getenv("PtolemyPath"); - if( heliosPath ){ - dataSource = heliosPath; - dataSource += "/.." + massData; - } - } - - - - -}; - -void Isotope::SetIso(int a, int z){ - this->A = a; - this->Z = z; - FindMassByAZ(a,z); -} - -void Isotope::SetIsoByName(string name){ - FindMassByName(name); -} - -void Isotope::FindMassByAZ(int A, int Z){ - string line; - int lineNum=0; - int list_A, list_Z; - - ifstream myfile; - int flag=0; - - setFileLines(); - - int numLineStart = fileStartLine; - int numLineEnd = fileEndLine; - - if ( A >= 50 && A < 100) numLineStart = lineMass050_099; - if ( A >=100 && A < 150) numLineStart = lineMass100_149; - if ( A >=150 && A < 200) numLineStart = lineMass150_199; - if ( A >=200 ) numLineStart = lineMass200; - - myfile.open(dataSource.c_str()); - - if (myfile.is_open()) { - while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ - list_Z = atoi((line.substr(10,5)).c_str()); - list_A = atoi((line.substr(15,5)).c_str()); - - if ( A == list_A && Z == list_Z) { - this->BEA = atof((line.substr(54,11)).c_str()); - this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A; - this->MassError = atof((line.substr(65,7)).c_str()); - string str = line.substr(20,2); - str.erase(remove(str.begin(), str.end(), ' '), str.end()); - this->Symbol = str; - - ostringstream ss; - ss << A << this->Symbol; - this->Name = ss.str(); - flag = 1; - }else if ( list_A > A) { - this->BEA = -404; - this->Mass = -404; - this->MassError = -404; - this->Symbol = "non"; - this->Name = "non"; - break; - } - - } - } - - if( this->Name == "1H" ) this->Name = "p"; - if( this->Name == "2H" ) this->Name = "d"; - if( this->Name == "3H" ) this->Name = "t"; - if( this->Name == "4He" ) this->Name = "a"; - - myfile.close(); - }else { - printf("Unable to open %s\n", dataSource.c_str()); - } -} - -void Isotope::FindMassByName(string name){ - - // done seperate the Mass number and the name - if( name == "n" ) { - this->Name = "1n"; - this->BEA = 0; - this->Mass = mn; - this->MassError = 0; - this->Name = "n"; - this->A = 1; - this->Z = 0; - return; - } - if( name == "p" ) name = "1H"; - if( name == "d" ) name = "2H"; - if( name == "t" ) name = "3H"; - if( name == "a" ) name = "4He"; - - string temp = name; - int lastDigit = 0; - - for(int i=0; temp[i]; i++){ - if(temp[i] == '0') lastDigit = i; - if(temp[i] == '1') lastDigit = i; - if(temp[i] == '2') lastDigit = i; - if(temp[i] == '3') lastDigit = i; - if(temp[i] == '4') lastDigit = i; - if(temp[i] == '5') lastDigit = i; - if(temp[i] == '6') lastDigit = i; - if(temp[i] == '7') lastDigit = i; - if(temp[i] == '8') lastDigit = i; - if(temp[i] == '9') lastDigit = i; - } - - this->Symbol = temp.erase(0, lastDigit +1); - //check is Symbol is 2 charaters, if not, add " " at the end - if( this->Symbol.length() == 1 ){ - this->Symbol = this->Symbol + " "; - } - - - temp = name; - int len = temp.length(); - temp = temp.erase(lastDigit+1, len); - - this->A = atoi(temp.c_str()); - //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); - - // find the nucleus in the data - string line; - int lineNum=0; - int list_A; - string list_symbol; - - ifstream myfile; - int flag=0; - - setFileLines(); - - int numLineStart = fileStartLine; - int numLineEnd = fileEndLine; - - if ( A >= 50 && A < 100) numLineStart = lineMass050_099; - if ( A >=100 && A < 150) numLineStart = lineMass100_149; - if ( A >=150 && A < 200) numLineStart = lineMass150_199; - if ( A >=200 ) numLineStart = lineMass200; - - myfile.open(dataSource.c_str()); - - if (myfile.is_open()) { - while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ - list_symbol = line.substr(20,2); - list_A = atoi((line.substr(15,5)).c_str()); - - //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str()); - - if ( this->A == list_A && this->Symbol == list_symbol) { - this->Z = atoi((line.substr(10,5)).c_str()); - this->BEA = atof((line.substr(54,11)).c_str()); - this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A; - this->MassError = atof((line.substr(65,7)).c_str()); - - string str = line.substr(20,2); - str.erase(remove(str.begin(), str.end(), ' '), str.end()); - this->Symbol = str; - - ostringstream ss; - ss << this->A << this->Symbol; - this->Name = ss.str(); - flag = 1; - }else if ( list_A > this->A) { - this->BEA = -404; - this->Mass = -404; - this->MassError = -404; - this->Symbol = "non"; - this->Name = "non"; - break; - } - - } - } - myfile.close(); - }else { - printf("Unable to open %s\n", dataSource.c_str()); - } -} - -double Isotope::CalSp(int Np, int Nn){ - Isotope nucleusD(A - Np - Nn, Z - Np); - - if( nucleusD.Mass != -404){ - return nucleusD.Mass + Nn*mn + Np*mp - this->Mass; - }else{ - return -404; - } -} - -double Isotope::CalSp2(int a, int z){ - Isotope nucleusD(A - a , Z - z); - Isotope nucleusS(a,z); - - if( nucleusD.Mass != -404 && nucleusS.Mass != -404){ - return nucleusD.Mass + nucleusS.Mass - this->Mass; - }else{ - return -404; - } -} - -int Isotope::TwoJ(int nShell){ - - switch(nShell){ - case 0: return 1; break; // 0s1/2 - case 1: return 3; break; // 0p3/2 - case 2: return 1; break; // 0p1/2 -- 8 - case 3: return 5; break; // 0d5/2 - case 4: return 1; break; // 1s1/2 - case 5: return 3; break; // 0d3/2 -- 20 - case 6: return 7; break; // 0f7/2 -- 28 - case 7: return 3; break; // 1p3/2 - case 8: return 1; break; // 1p1/2 - case 9: return 5; break; // 0f5/2 -- 40 - case 10: return 9; break; // 0g9/2 -- 50 - case 11: return 7; break; // 0g7/2 - case 12: return 5; break; // 1d5/2 - case 13: return 11; break; // 0h11/2 - case 14: return 3; break; // 1d3/2 - case 15: return 1; break; // 2s1/2 -- 82 - case 16: return 9; break; // 0h9/2 - case 17: return 7; break; // 1f7/2 - case 18: return 13; break; // 0i13/2 - case 19: return 3; break; // 2p3/2 - case 20: return 5; break; // 1f5/2 - case 21: return 1; break; // 1p1/2 -- 126 - case 22: return 9; break; // 1g9/2 - case 23: return 11; break; // 0i11/2 - case 24: return 15; break; // 0j15/2 - case 25: return 5; break; // 2d5/2 - case 26: return 1; break; // 3s1/2 - case 27: return 3; break; // 2d3/2 - case 28: return 7; break; // 1g7/2 - } - - return 0; -} - -string Isotope::Orbital(int nShell){ - - switch(nShell){ - case 0: return "0s1 "; break; // - case 1: return "0p3 "; break; // - case 2: return "0p1 "; break; //-- 8 - case 3: return "0d5 "; break; // - case 4: return "1s1 "; break; // - case 5: return "0d3 "; break; //-- 20 - case 6: return "0f7 "; break; //-- 28 - case 7: return "1p3 "; break; // - case 8: return "1p1 "; break; // - case 9: return "0f5 "; break; //-- 40 - case 10: return "0g9 "; break; //-- 50 - case 11: return "0g7 "; break; // - case 12: return "1d5 "; break; // - case 13: return "0h11"; break; // - case 14: return "1d3 "; break; // - case 15: return "2s1 "; break; //-- 82 - case 16: return "0h9 "; break; // - case 17: return "1f7 "; break; // - case 18: return "0i13"; break; // - case 19: return "2p3 "; break; // - case 20: return "1f5 "; break; // - case 21: return "1p1 "; break; //-- 126 - case 22: return "1g9 "; break; // - case 23: return "0i11"; break; // - case 24: return "0j15"; break; // - case 25: return "2d5 "; break; // - case 26: return "3s1 "; break; // - case 27: return "2d3 "; break; // - case 28: return "1g7 "; break; // - } - - return "nan"; -} - -void Isotope::ListShell(){ - - if( Mass < 0 ) return; - - int n = A-Z; - int p = Z; - - int k = min(n,p); - int nMagic = 0; - for( int i = 0; i < 7; i++){ - if( magic(i) < k && k <= magic(i+1) ){ - nMagic = i; - break; - } - } - - int coreShell = magicShellID(nMagic-1); - int coreZ1 = magic(nMagic-1); - int coreZ2 = magic(nMagic); - - Isotope core1( 2*coreZ1, coreZ1); - Isotope core2( 2*coreZ2, coreZ2); - - printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str()); - printf(" || "); - int t = max(n,p); - int nShell = 0; - do{ - int occ = TwoJ(nShell)+1; - if( nShell > coreShell) { - printf("%4s", Orbital(nShell).c_str()); - if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ - printf("|"); - }else{ - printf(","); - } - } - t = t - occ; - nShell++; - }while( t > 0 && nShell < 29); - for( int i = 1; i <= 6; i++){ - if (nShell < 28) { - printf("%4s,", Orbital(nShell).c_str()); - }else if( nShell == 28 ) { - printf("%4s", Orbital(nShell).c_str()); - } - nShell ++; - } - if (nShell < 29) printf("%4s", Orbital(nShell).c_str()); - printf("\n"); - - - printf(" Z = %3d || ", p); - nShell = 0; - do{ - int occ = TwoJ(nShell)+1; - if( nShell > coreShell ){ - if( p > occ ) { - printf("%-4d", occ); - if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ - printf("|"); - }else{ - printf(","); - } - }else{ - printf("%-4d", p); - } - } - p = p - occ; - nShell++; - }while( p > 0 && nShell < 29); - printf("\n"); - - printf(" N = %3d || ", n); - nShell = 0; - do{ - int occ = TwoJ(nShell)+1; - if ( nShell > coreShell ){ - if( n > occ ) { - printf("%-4d", occ); - if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ - printf("|"); - }else{ - printf(","); - } - }else{ - printf("%-4d", n); - } - } - n = n - occ; - nShell++; - }while( n > 0 && nShell < 29); - printf("\n"); - - printf("------------------ \n"); -} - - - -void Isotope::Print(){ - - if (Mass > 0){ - - findHeliosPath(); - - printf(" using mass data : %s \n", dataSource.c_str()); - printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.); - printf(" total BE : %12.5f MeV\n",BEA*A/1000.); - printf(" mass in amu : %12.5f u\n",Mass/amu); - printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu); - printf("-------------- Seperation energy \n"); - printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1)); - printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1)); - printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n", CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2)); - printf(" S4p: %8.4f| S4n: %8.4f| \n", CalSp(4, 0), CalSp(0, 4)); - - }else{ - printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A); - } - - -} - -#endif +/*********************************************************************** + * + * This is Isotope.h, To extract the isotope mass from massXX.txt + * + *------------------------------------------------------- + * created by Ryan (Tsz Leung) Tang, Nov-18, 2018 + * email: goluckyryan@gmail.com + * ********************************************************************/ + + +#ifndef ISOTOPE_H +#define ISOTOPE_H + +#include +#include +#include +#include +#include +#include +#include "constant.h" // amu +#include //atoi +#include +using namespace std; + +string massData="/Cleopatra/mass20.txt"; + +// about the mass**.txt +// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu +// mass excess uncertaintly +// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons +// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons + +class Isotope { +public: + int A, Z; + double Mass, MassError, BEA; + string Name, Symbol; + string dataSource; + + Isotope(){findHeliosPath();}; + Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; + Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; + + void SetIso(int a, int z); + void SetIsoByName(string name); + + double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal + double CalSp2(int a, int z); // this is for (a,z) nucleus removal + + double CalBeta(double T){ + double Etot = Mass + T; + double gamma = 1 + T/Mass; + double beta = sqrt(1 - 1 / gamma / gamma ) ; + return beta; + } + + void Print(); + void ListShell(); + +private: + void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; + void FindMassByName(string name); // give Z, mass, massError, BEA; + + int TwoJ(int nShell); + string Orbital(int nShell); + int magic(int i){ + switch (i){ + case 0: return 2; break; + case 1: return 8; break; + case 2: return 20; break; + case 3: return 28; break; + case 4: return 40; break; + case 5: return 50; break; + case 6: return 82; break; + case 7: return 128; break; + } + return 0; + } + + int magicShellID(int i){ + switch (i){ + case 0: return 0; break; + case 1: return 2; break; + case 2: return 5; break; + case 3: return 6; break; + case 4: return 9; break; + case 5: return 10; break; + case 6: return 15; break; + case 7: return 21; break; + } + return 0; + } + + int fileStartLine; + int fileEndLine; + int lineMass050_099; + int lineMass100_149; + int lineMass150_199; + int lineMass200; + + + void setFileLines(){ + fileStartLine = 37; + fileEndLine = 3594; + + lineMass050_099 = 466; + lineMass100_149 = 1160; + lineMass150_199 = 1994; + lineMass200 = 2774; + } + + char * heliosPath; + bool isFindOnce; + + void findHeliosPath(){ + heliosPath = getenv("HELIOSSYS"); + if( heliosPath ){ + dataSource = heliosPath; + dataSource += "/analysis" + massData; + }else{ + dataSource = ".." + massData; + } + } + + + + +}; + +void Isotope::SetIso(int a, int z){ + this->A = a; + this->Z = z; + FindMassByAZ(a,z); +} + +void Isotope::SetIsoByName(string name){ + FindMassByName(name); +} + +void Isotope::FindMassByAZ(int A, int Z){ + string line; + int lineNum=0; + int list_A, list_Z; + + ifstream myfile; + int flag=0; + + setFileLines(); + + int numLineStart = fileStartLine; + int numLineEnd = fileEndLine; + + if ( A >= 50 && A < 100) numLineStart = lineMass050_099; + if ( A >=100 && A < 150) numLineStart = lineMass100_149; + if ( A >=150 && A < 200) numLineStart = lineMass150_199; + if ( A >=200 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_Z = atoi((line.substr(10,5)).c_str()); + list_A = atoi((line.substr(15,5)).c_str()); + + if ( A == list_A && Z == list_Z) { + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + break; + } + + } + } + + if( this->Name == "1H" ) this->Name = "p"; + if( this->Name == "2H" ) this->Name = "d"; + if( this->Name == "3H" ) this->Name = "t"; + if( this->Name == "4He" ) this->Name = "a"; + + myfile.close(); + }else { + printf("Unable to open %s\n", dataSource.c_str()); + } +} + +void Isotope::FindMassByName(string name){ + + // done seperate the Mass number and the name + if( name == "n" ) { + this->Name = "1n"; + this->BEA = 0; + this->Mass = mn; + this->MassError = 0; + this->Name = "n"; + this->A = 1; + this->Z = 0; + return; + } + if( name == "p" ) name = "1H"; + if( name == "d" ) name = "2H"; + if( name == "t" ) name = "3H"; + if( name == "a" ) name = "4He"; + + string temp = name; + int lastDigit = 0; + + for(int i=0; temp[i]; i++){ + if(temp[i] == '0') lastDigit = i; + if(temp[i] == '1') lastDigit = i; + if(temp[i] == '2') lastDigit = i; + if(temp[i] == '3') lastDigit = i; + if(temp[i] == '4') lastDigit = i; + if(temp[i] == '5') lastDigit = i; + if(temp[i] == '6') lastDigit = i; + if(temp[i] == '7') lastDigit = i; + if(temp[i] == '8') lastDigit = i; + if(temp[i] == '9') lastDigit = i; + } + + this->Symbol = temp.erase(0, lastDigit +1); + //check is Symbol is 2 charaters, if not, add " " at the end + if( this->Symbol.length() == 1 ){ + this->Symbol = this->Symbol + " "; + } + + + temp = name; + int len = temp.length(); + temp = temp.erase(lastDigit+1, len); + + this->A = atoi(temp.c_str()); + //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); + + // find the nucleus in the data + string line; + int lineNum=0; + int list_A; + string list_symbol; + + ifstream myfile; + int flag=0; + + setFileLines(); + + int numLineStart = fileStartLine; + int numLineEnd = fileEndLine; + + if ( A >= 50 && A < 100) numLineStart = lineMass050_099; + if ( A >=100 && A < 150) numLineStart = lineMass100_149; + if ( A >=150 && A < 200) numLineStart = lineMass150_199; + if ( A >=200 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_symbol = line.substr(20,2); + list_A = atoi((line.substr(15,5)).c_str()); + + //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str()); + + if ( this->A == list_A && this->Symbol == list_symbol) { + this->Z = atoi((line.substr(10,5)).c_str()); + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << this->A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > this->A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + break; + } + + } + } + myfile.close(); + }else { + printf("Unable to open %s\n", dataSource.c_str()); + } +} + +double Isotope::CalSp(int Np, int Nn){ + Isotope nucleusD(A - Np - Nn, Z - Np); + + if( nucleusD.Mass != -404){ + return nucleusD.Mass + Nn*mn + Np*mp - this->Mass; + }else{ + return -404; + } +} + +double Isotope::CalSp2(int a, int z){ + Isotope nucleusD(A - a , Z - z); + Isotope nucleusS(a,z); + + if( nucleusD.Mass != -404 && nucleusS.Mass != -404){ + return nucleusD.Mass + nucleusS.Mass - this->Mass; + }else{ + return -404; + } +} + +int Isotope::TwoJ(int nShell){ + + switch(nShell){ + case 0: return 1; break; // 0s1/2 + case 1: return 3; break; // 0p3/2 + case 2: return 1; break; // 0p1/2 -- 8 + case 3: return 5; break; // 0d5/2 + case 4: return 1; break; // 1s1/2 + case 5: return 3; break; // 0d3/2 -- 20 + case 6: return 7; break; // 0f7/2 -- 28 + case 7: return 3; break; // 1p3/2 + case 8: return 1; break; // 1p1/2 + case 9: return 5; break; // 0f5/2 -- 40 + case 10: return 9; break; // 0g9/2 -- 50 + case 11: return 7; break; // 0g7/2 + case 12: return 5; break; // 1d5/2 + case 13: return 11; break; // 0h11/2 + case 14: return 3; break; // 1d3/2 + case 15: return 1; break; // 2s1/2 -- 82 + case 16: return 9; break; // 0h9/2 + case 17: return 7; break; // 1f7/2 + case 18: return 13; break; // 0i13/2 + case 19: return 3; break; // 2p3/2 + case 20: return 5; break; // 1f5/2 + case 21: return 1; break; // 1p1/2 -- 126 + case 22: return 9; break; // 1g9/2 + case 23: return 11; break; // 0i11/2 + case 24: return 15; break; // 0j15/2 + case 25: return 5; break; // 2d5/2 + case 26: return 1; break; // 3s1/2 + case 27: return 3; break; // 2d3/2 + case 28: return 7; break; // 1g7/2 + } + + return 0; +} + +string Isotope::Orbital(int nShell){ + + switch(nShell){ + case 0: return "0s1 "; break; // + case 1: return "0p3 "; break; // + case 2: return "0p1 "; break; //-- 8 + case 3: return "0d5 "; break; // + case 4: return "1s1 "; break; // + case 5: return "0d3 "; break; //-- 20 + case 6: return "0f7 "; break; //-- 28 + case 7: return "1p3 "; break; // + case 8: return "1p1 "; break; // + case 9: return "0f5 "; break; //-- 40 + case 10: return "0g9 "; break; //-- 50 + case 11: return "0g7 "; break; // + case 12: return "1d5 "; break; // + case 13: return "0h11"; break; // + case 14: return "1d3 "; break; // + case 15: return "2s1 "; break; //-- 82 + case 16: return "0h9 "; break; // + case 17: return "1f7 "; break; // + case 18: return "0i13"; break; // + case 19: return "2p3 "; break; // + case 20: return "1f5 "; break; // + case 21: return "1p1 "; break; //-- 126 + case 22: return "1g9 "; break; // + case 23: return "0i11"; break; // + case 24: return "0j15"; break; // + case 25: return "2d5 "; break; // + case 26: return "3s1 "; break; // + case 27: return "2d3 "; break; // + case 28: return "1g7 "; break; // + } + + return "nan"; +} + +void Isotope::ListShell(){ + + if( Mass < 0 ) return; + + int n = A-Z; + int p = Z; + + int k = min(n,p); + int nMagic = 0; + for( int i = 0; i < 7; i++){ + if( magic(i) < k && k <= magic(i+1) ){ + nMagic = i; + break; + } + } + + int coreShell = magicShellID(nMagic-1); + int coreZ1 = magic(nMagic-1); + int coreZ2 = magic(nMagic); + + Isotope core1( 2*coreZ1, coreZ1); + Isotope core2( 2*coreZ2, coreZ2); + + printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str()); + printf(" || "); + int t = max(n,p); + int nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if( nShell > coreShell) { + printf("%4s", Orbital(nShell).c_str()); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + } + t = t - occ; + nShell++; + }while( t > 0 && nShell < 29); + for( int i = 1; i <= 6; i++){ + if (nShell < 28) { + printf("%4s,", Orbital(nShell).c_str()); + }else if( nShell == 28 ) { + printf("%4s", Orbital(nShell).c_str()); + } + nShell ++; + } + if (nShell < 29) printf("%4s", Orbital(nShell).c_str()); + printf("\n"); + + + printf(" Z = %3d || ", p); + nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if( nShell > coreShell ){ + if( p > occ ) { + printf("%-4d", occ); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + }else{ + printf("%-4d", p); + } + } + p = p - occ; + nShell++; + }while( p > 0 && nShell < 29); + printf("\n"); + + printf(" N = %3d || ", n); + nShell = 0; + do{ + int occ = TwoJ(nShell)+1; + if ( nShell > coreShell ){ + if( n > occ ) { + printf("%-4d", occ); + if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){ + printf("|"); + }else{ + printf(","); + } + }else{ + printf("%-4d", n); + } + } + n = n - occ; + nShell++; + }while( n > 0 && nShell < 29); + printf("\n"); + + printf("------------------ \n"); +} + + + +void Isotope::Print(){ + + if (Mass > 0){ + + findHeliosPath(); + + printf(" using mass data : %s \n", dataSource.c_str()); + printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.); + printf(" total BE : %12.5f MeV\n",BEA*A/1000.); + printf(" mass in amu : %12.5f u\n",Mass/amu); + printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu); + printf("-------------- Seperation energy \n"); + printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1)); + printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1)); + printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n", CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2)); + printf(" S4p: %8.4f| S4n: %8.4f| \n", CalSp(4, 0), CalSp(0, 4)); + + }else{ + printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A); + } + + +} + +#endif diff --git a/Cleopatra/PlotTGraphTObjArray.C b/Cleopatra/PlotTGraphTObjArray.C index 74c3810..c65fbdc 100644 --- a/Cleopatra/PlotTGraphTObjArray.C +++ b/Cleopatra/PlotTGraphTObjArray.C @@ -29,8 +29,9 @@ int main (int argc, char *argv[]) { printf("==================== Plot Results in ROOT =======================\n"); printf("=================================================================\n"); - if(argc != 2) { - printf("Usage: ./PlotTGraphTObjArray root_file\n"); + if(argc < 2) { + printf("Usage: ./PlotTGraphTObjArray root_file [savePNG]\n"); + printf(" savePNG : 1 or 0, default is 0\n"); exit(0); }else{ printf("From file : %s \n", argv[1]); @@ -38,13 +39,20 @@ int main (int argc, char *argv[]) { } string readFile = argv[1]; - - TApplication app ("app", &argc, argv); + bool isSavePNG = false; + if( argc >= 3) isSavePNG = atoi(argv[2]); - PlotTGraphTObjArray(readFile); - - app.Run(); //anything after this line is not running + if( isSavePNG ){ + + PlotTGraphTObjArray(readFile, true); + + }else{ + TApplication app ("app", &argc, argv); + + PlotTGraphTObjArray(readFile); + app.Run(); //anything after this line is not running + } return 0; } diff --git a/Cleopatra/PlotTGraphTObjArray.h b/Cleopatra/PlotTGraphTObjArray.h index 95b2a21..8c3858d 100644 --- a/Cleopatra/PlotTGraphTObjArray.h +++ b/Cleopatra/PlotTGraphTObjArray.h @@ -21,12 +21,13 @@ #include #include #include +#include #include #include #include #include -void PlotTGraphTObjArray(TString rootFileName){ +void PlotTGraphTObjArray(TString rootFileName, bool isSavePNG = false){ TFile * file = new TFile(rootFileName, "READ"); @@ -40,19 +41,22 @@ void PlotTGraphTObjArray(TString rootFileName){ TCanvas * cPlots = new TCanvas("cPlots", "Ptolemy Results", 0, 0, 800, 600); cPlots->SetLogy(); - TLegend * legend = new TLegend( 0.6, 0.2, 0.9, 0.4); + TLegend * legend = new TLegend( 0.6, 0.6, 0.9, 0.9); //x1, y1, x2, y2 const int n = gList->GetLast() + 1 ; TGraph * gr[n]; + const std::vector color = {kBlack, kRed, kOrange+7, kYellow+1, kGreen+2, kGreen+3, kBlue, kBlue+3, kMagenta, kMagenta+2, kGray+2, kRed+2, kYellow+3 }; + short nColor = color.size(); + //Get minimum, maximum Y double maxY = 0; double minY = 10000000; for ( int i = 0; i < n ; i++){ gr[i] = (TGraph *) gList->At(i); - gr[i]->SetLineColor(i+1); + gr[i]->SetLineColor(color[ i % nColor]); gr[i]->GetXaxis()->SetTitle("#theta_{CM} [deg]"); gr[i]->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/sr]"); @@ -68,8 +72,6 @@ void PlotTGraphTObjArray(TString rootFileName){ if( maxy > maxY ) maxY = maxy; } - - for ( int i = 0; i < n ; i++){ gr[i]->Draw("same"); @@ -82,7 +84,19 @@ void PlotTGraphTObjArray(TString rootFileName){ } legend->Draw(); + cPlots->SetGrid(1,1); + cPlots->Update(); cPlots->Draw(); + + if( isSavePNG ){ + TDatime dateTime; + TString outPNGName = Form("Xsec_%d%02d%02d_%06d.png", dateTime.GetYear(), dateTime.GetMonth(), dateTime.GetDay(), dateTime.GetTime()); + + cPlots->SaveAs(outPNGName); + printf("%s\n", outPNGName.Data()); + + gROOT->ProcessLine(".q"); + } }