diff --git a/.gitignore b/.gitignore index abdd51d..bb7ae6f 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ root_data *.in *.out *.txt +*.csv Cleopatra/ExtractXSec Cleopatra/ExtractXSecFromText @@ -25,4 +26,6 @@ Cleopatra/IsotopeShort Cleopatra/PlotSimulation Cleopatra/PlotTGraphTObjArray Cleopatra/SimAlpha -Cleopatra/SimTransfer \ No newline at end of file +Cleopatra/SimTransfer + +__pycache__ \ No newline at end of file diff --git a/Cleopatra/ClassIsotope.h b/Cleopatra/ClassIsotope.h index afabaa9..4416dba 100644 --- a/Cleopatra/ClassIsotope.h +++ b/Cleopatra/ClassIsotope.h @@ -20,9 +20,8 @@ #include "constant.h" // amu #include //atoi #include -using namespace std; -string massData="/Cleopatra/mass20.txt"; +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 @@ -34,15 +33,17 @@ class Isotope { public: int A, Z; double Mass, MassError, BEA; - string Name, Symbol; - string dataSource; + std::string Name, Symbol; + std::string dataSource; - Isotope(){findHeliosPath();}; - Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; - Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; + Isotope(){dataSource = massData;}; + Isotope(int a, int z){ dataSource = massData; SetIso(a,z); }; + Isotope(std::string name){ dataSource = massData; SetIsoByName(name); }; + + void SetMassTablePath(std::string path){ dataSource = path;} void SetIso(int a, int z); - void SetIsoByName(string name); + void SetIsoByName(std::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 @@ -59,10 +60,10 @@ public: private: void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; - void FindMassByName(string name); // give Z, mass, massError, BEA; + void FindMassByName(std::string name); // give Z, mass, massError, BEA; int TwoJ(int nShell); - string Orbital(int nShell); + std::string Orbital(int nShell); int magic(int i){ switch (i){ case 0: return 2; break; @@ -109,21 +110,7 @@ private: lineMass200 = 2774; } - char * heliosPath; bool isFindOnce; - - void findHeliosPath(){ - heliosPath = getenv("HELIOSSYS"); - if( heliosPath ){ - dataSource = heliosPath; - dataSource += "/analysis" + massData; - }else{ - dataSource = ".." + massData; - } - } - - - }; @@ -133,16 +120,16 @@ void Isotope::SetIso(int a, int z){ FindMassByAZ(a,z); } -void Isotope::SetIsoByName(string name){ +void Isotope::SetIsoByName(std::string name){ FindMassByName(name); } void Isotope::FindMassByAZ(int A, int Z){ - string line; + std::string line; int lineNum=0; int list_A, list_Z; - ifstream myfile; + std::ifstream myfile; int flag=0; setFileLines(); @@ -171,11 +158,11 @@ void Isotope::FindMassByAZ(int A, int 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); + std::string str = line.substr(20,2); str.erase(remove(str.begin(), str.end(), ' '), str.end()); this->Symbol = str; - ostringstream ss; + std::ostringstream ss; ss << A << this->Symbol; this->Name = ss.str(); flag = 1; @@ -202,7 +189,7 @@ void Isotope::FindMassByAZ(int A, int Z){ } } -void Isotope::FindMassByName(string name){ +void Isotope::FindMassByName(std::string name){ // done seperate the Mass number and the name if( name == "n" ) { @@ -220,7 +207,7 @@ void Isotope::FindMassByName(string name){ if( name == "t" ) name = "3H"; if( name == "a" ) name = "4He"; - string temp = name; + std::string temp = name; int lastDigit = 0; for(int i=0; temp[i]; i++){ @@ -251,12 +238,12 @@ void Isotope::FindMassByName(string name){ //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); // find the nucleus in the data - string line; + std::string line; int lineNum=0; int list_A; - string list_symbol; + std::string list_symbol; - ifstream myfile; + std::ifstream myfile; int flag=0; setFileLines(); @@ -289,11 +276,11 @@ void Isotope::FindMassByName(string name){ 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); + std::string str = line.substr(20,2); str.erase(remove(str.begin(), str.end(), ' '), str.end()); this->Symbol = str; - ostringstream ss; + std::ostringstream ss; ss << this->A << this->Symbol; this->Name = ss.str(); flag = 1; @@ -372,7 +359,7 @@ int Isotope::TwoJ(int nShell){ return 0; } -string Isotope::Orbital(int nShell){ +std::string Isotope::Orbital(int nShell){ switch(nShell){ case 0: return "0s1 "; break; // @@ -416,7 +403,7 @@ void Isotope::ListShell(){ int n = A-Z; int p = Z; - int k = min(n,p); + int k = std::min(n,p); int nMagic = 0; for( int i = 0; i < 7; i++){ if( magic(i) < k && k <= magic(i+1) ){ @@ -434,7 +421,7 @@ void Isotope::ListShell(){ printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str()); printf(" || "); - int t = max(n,p); + int t = std::max(n,p); int nShell = 0; do{ int occ = TwoJ(nShell)+1; @@ -512,8 +499,6 @@ 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.); diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index c3cf3a4..2534779 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -22,7 +22,7 @@ class TransferReaction { public: TransferReaction(){Inititization();}; - TransferReaction(string configFile, unsigned short ID = 0); + TransferReaction(std::string configFile, unsigned short ID = 0); TransferReaction(int beamA, int beamZ, int targetA, int targetZ, int recoilA, int recoilZ, float beamEnergy_AMeV); @@ -84,7 +84,7 @@ private: Recoil recoil; ReactionConfig config; - string nameA, namea, nameb, nameB; + std::string nameA, namea, nameb, nameB; double thetaIN, phiIN; double mA, ma, mb, mB; @@ -108,7 +108,7 @@ private: }; -TransferReaction::TransferReaction(string configFile, unsigned short ID){ +TransferReaction::TransferReaction(std::string configFile, unsigned short ID){ Inititization(); SetReactionFromFile(configFile, ID); } @@ -227,7 +227,7 @@ void TransferReaction::SetExB(double Ex){ isReady = false; } -void TransferReaction::SetReactionFromFile(string configFile, unsigned short ID){ +void TransferReaction::SetReactionFromFile(std::string configFile, unsigned short ID){ if( config.LoadReactionConfig(configFile) ){ diff --git a/Cleopatra/FindThetaCM.C b/Cleopatra/FindThetaCM.C index e2fc9e3..db33abc 100644 --- a/Cleopatra/FindThetaCM.C +++ b/Cleopatra/FindThetaCM.C @@ -34,8 +34,8 @@ int main(int argc, char *argv[]){ double Ex = 0; double xRatio = 0.95; int nDiv = 1; - string reactionTxt = "reactionConfig.txt"; - string detGeoTxt = "detectorGeo.txt"; + std::string reactionTxt = "reactionConfig.txt"; + std::string detGeoTxt = "detectorGeo.txt"; int ID = 0; if ( argc >= 2 ) Ex = atof(argv[1]); diff --git a/Cleopatra/FindThetaCM.h b/Cleopatra/FindThetaCM.h index 6956ab6..d4e4135 100644 --- a/Cleopatra/FindThetaCM.h +++ b/Cleopatra/FindThetaCM.h @@ -94,7 +94,7 @@ void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95, int iDet = array.nDet; double length = array.detLength; - vector midPos; + std::vector midPos; for(int i = 0; i < iDet; i++){ if( array.firstPos > 0 ){ diff --git a/Cleopatra/SimTransfer.C b/Cleopatra/SimTransfer.C index 0274076..bab3b06 100644 --- a/Cleopatra/SimTransfer.C +++ b/Cleopatra/SimTransfer.C @@ -135,9 +135,9 @@ void Transfer( dwbaExList = (TMacro *) distFile->FindObjectAny("ExList"); int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; for(int i = 1; i <= numEx ; i++){ - string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); + std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); if( temp[0] == '/' ) continue; - vector tempStr = AnalysisLib::SplitStr(temp, " "); + std::vector tempStr = AnalysisLib::SplitStr(temp, " "); exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); } @@ -330,7 +330,7 @@ void Transfer( //--- cal modified f TObjArray * fxList = new TObjArray(); TGraph ** fx = new TGraph*[numEx]; - vector px, py; + std::vector px, py; int countfx = 0; for( int j = 0 ; j < numEx; j++){ double a = helios.GetDetRadius(); @@ -692,8 +692,8 @@ int main (int argc, char *argv[]) { //run Armory/Check_Simulation if( isPlot ){ - ifstream file_in; - file_in.open("../Cleopatra/Check_Simulation.C", ios::in); + std::ifstream file_in; + file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in); if( file_in){ printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data()); TString cmd; diff --git a/Cleopatra/potentials.h b/Cleopatra/potentials.h index be623be..e39c935 100644 --- a/Cleopatra/potentials.h +++ b/Cleopatra/potentials.h @@ -26,7 +26,7 @@ void PrintPotential(){ /// A B C D E F G H I J K L M N O P Q R S T U V W X Y Z /// 1 1 1 1 0 0 1 1 0 0 1 1 1 0 0 1 1 0 0 0 0 1 0 1 1 1 -string potentialRef(string name){ +std::string potentialRef(std::string name){ //======== Deuteron if( name == "A" ){ @@ -1034,7 +1034,7 @@ bool BassaniPicardPotential(int A, int Z, double E){ return true; } -bool CallPotential(string potName, int A, int Z, double E, int Zproj){ +bool CallPotential(std::string potName, int A, int Z, double E, int Zproj){ bool okFlag = false; if( potName == "A") okFlag = AnCaiPotential(A, Z, E); diff --git a/WebSimHelper/EnergyLevelsPlot.js b/WebSimHelper/EnergyLevelsPlot.js new file mode 100644 index 0000000..31204c8 --- /dev/null +++ b/WebSimHelper/EnergyLevelsPlot.js @@ -0,0 +1,239 @@ +let energy = []; +let jpi = []; +let Name; +let A; +let Sym; + +function breakdownName(str) { + const match = str.match(/^(\d+)([a-zA-Z]+)$/); + if (match) { + const numberPart = parseInt(match[1]); + const stringPart = match[2]; + return { numberPart, stringPart }; + } else { + return null; // If the input string doesn't match the expected format + } +} + +let Sn = 999; +let Sp = 999; +let Sa = 999; + +function GetData(){ + + Name = document.getElementById('ASym').value; + let maxEx = parseFloat(document.getElementById('maxEx').value); + + let str = 'displayIsoData.py?ASym=' + Name + "&maxEx=" + maxEx; + + let client = new XMLHttpRequest(); + client.onreadystatechange = function() { + let haha = client.responseText.split('\n'); + + jpi = []; + energy = []; + haha.forEach(line =>{ + + // console.log(line) + if( line.includes("Sn:") ) { + let pos1 = line.indexOf("Sn:"); + let pos2 = line.indexOf("MeV"); + Sn = parseFloat(line.substring(pos1+3, pos2)); + } + if( line.includes("Sp:") ) { + let pos1 = line.indexOf("Sp:"); + let pos2 = line.indexOf("MeV"); + Sp = parseFloat(line.substring(pos1+3, pos2)); + } + if( line.includes("Sa:") ) { + let pos1 = line.indexOf("Sa:"); + let pos2 = line.indexOf("MeV"); + Sa = parseFloat(line.substring(pos1+3, pos2)); + } + + if( line.includes(" fontSizeMeV) { + count++; + } + } + if (count === l) { + noOverlap = true; + } + loop++; + } + + for (let i = 0; i < l; i++) { + fig.data.push({ + x: [0, 1], + y: [energy[i], energy[i]], + mode: 'lines', + line: { color: 'black', width: 1 } + }); + + fig.data.push({ + x: [1.03, 1.1, 1.19], + y: [energy[i], ypos[i], ypos[i]], + mode: 'lines', + line: { color: 'gray', width: 1 } + }); + + // console.log(energy[i]+ ", " + ypos[i]); + + fig.layout.annotations.push({ + x: 1.2, + y: ypos[i], + text: `${energy[i].toFixed(3)}, ${jpi[i]}`, + xanchor: 'left', + font: { size: fontSize }, + showarrow: false + }); + } + + console.log("Sn: " + Sn); + console.log("Sp: " + Sp); + console.log("Sa: " + Sa); + + let leftPos = -0.8; + + fig.data.push({ + x: [leftPos, 1], + y: [Sn, Sn], + mode: 'lines', + line: { color: 'blue', width: 1 } + }); + fig.layout.annotations.push({ + x: leftPos, + y: Sn + fontSizeMeV/2, + text: `${'Sn:'+Sn.toFixed(3)}`, + xanchor: 'left', + font: { size: fontSize, color: 'blue' }, + showarrow: false + }); + + fig.data.push({ + x: [leftPos, 1], + y: [Sp, Sp], + mode: 'lines', + line: { color: 'red', width: 1 } + }); + fig.layout.annotations.push({ + x: leftPos, + y: Sp + fontSizeMeV/2, + text: `${'Sp:'+Sp.toFixed(3)}`, + xanchor: 'left', + font: { size: fontSize, color: 'red' }, + showarrow: false + }); + + fig.data.push({ + x: [leftPos, 1], + y: [Sa, Sa], + mode: 'lines', + line: { color: 'purple', width: 1 } + }); + fig.layout.annotations.push({ + x: leftPos, + y: Sa + fontSizeMeV/2, + text: `${'Sa:'+Sa.toFixed(3)}`, + xanchor: 'left', + font: { size: fontSize, color: 'purple' }, + showarrow: false + }); + + // let NameYPos = (parseFloat(maxEx) + 2*fontSizeMeV); + // console.log(NameYPos); + + let name2 = breakdownName(Name); + + fig.layout.annotations.push({ + x: 0.5, + y: (maxEx + 1), + text: "" + name2.numberPart +"" + name2.stringPart, + font: { size: 2 * fontSize }, + showarrow: false + }); + + // Create the plot + Plotly.newPlot('Plot_Levels', fig.data, fig.layout); + +} \ No newline at end of file diff --git a/WebSimHelper/README.md b/WebSimHelper/README.md new file mode 100644 index 0000000..cbdb2e7 --- /dev/null +++ b/WebSimHelper/README.md @@ -0,0 +1,55 @@ +# Introduction + +This is a web inteface for the HELIOS/SOLARIS simulation. Its purpose is NOT to replace the Simulation_Helper.C in the origin digios repository. +It is simply provide a more easy accessible way to do simulation. + +# Installation in Apache2 + +Assume the parant SOLARIS_ANALYSIS is in the home folder +add a symbolic link + +```sh +$cd /var/www/html +$ln -s ~/SOLARIS_ANALYSIS SOLARIS +``` + +I want localhost/SOLARIS map to /var/www/html/SOLARIS/WebSimHelper, in the apache config + +```sh +$cd /etc/apache2/sit-available +$touch SOLARIS.conf +``` + +inside SOLARIS.conf + +```sh + + ServerAdmin rtang@anl.gov + DocumentRoot /var/www/html/ + ServerName localhost + + #map localhost/SOLARIS to /var/www/html/SOLARIS/WebSimHelper + Alias /SOLARIS /var/www/html/SOLARIS/WebSimHelper + + #set the directory properties + + Options Indexes FollowSymLinks + AllowOverride None + Require all granted + Options +ExecCGI + AddHandler cgi-script .cgi .py + + + ErrorLog ${APACHE_LOG_DIR}/error.log + CustomLog ${APACHE_LOG_DIR}/access.log combined + + +``` + +then enable the site + +```sh +$sudo a2ensite SOLARIS.conf +$sudo systemctl restart apach2.service +``` + diff --git a/WebSimHelper/displayIsoData.py b/WebSimHelper/displayIsoData.py new file mode 100755 index 0000000..e5b9bfc --- /dev/null +++ b/WebSimHelper/displayIsoData.py @@ -0,0 +1,29 @@ +#!/usr/bin/env /usr/bin/python3 + +import isotopeLib +import cgi + +form = cgi.FieldStorage() + +ASym = form.getvalue('ASym') +maxEx = form.getvalue('maxEx') + +print( "Content-type:text/html\r\n\r\n") +print("") +print("") +print("") + +isotopeLib.PrintIsoWeb(ASym) + +if maxEx == "can be omitted" or float(maxEx) <= 0: + maxEx = -1 +else: + isotopeLib.PrintIsoExWeb(ASym, float(maxEx)) + + +print("") +print("") + + + + diff --git a/WebSimHelper/heliosmatics.html b/WebSimHelper/heliosmatics.html new file mode 100644 index 0000000..57774d5 --- /dev/null +++ b/WebSimHelper/heliosmatics.html @@ -0,0 +1,470 @@ + + + +Heliosmatics + + + + + + + + + +

HELIOSmatics

+ + + +

24F(d,p)25F@10MeV/u

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
Beam (A): Beam Ex:MeV
Target (a):
Light (b): Q-value:2.057MeV
Heavy (B):25F
+ +

+ +

+ + + + + + + +
+ HELIOS + + SOLARIS + + ISS +
+

+ + + + + + + + + + + + + + + + + + + + + +
B-field (abs.):T
Beam Energy:MeV/u
+ +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Ex [MeV] θcm[deg]Eb[MeV]Zb0[mm]Zb[mm]2ρb[mm]θLabb[deg]Tb[ns]EB[MeV]θLabB[deg]ZB0/2[mm]2ρB[mm]
+ +

+ + + + + + + + + + + + + + + + + + + +
θCM:deg
Array Pos:mm
Recoil Pos:mm
+ + + + + + + + +
Recoil radius, inner [mm]: outter [mm]:
+ +

+ + +
+
+ + + + + + + + + + + + +
+
+
 
 zRange can be changed by Array position.
+ + + + + + + +
eRange:MeV
+
 
+
+
+ + + + + + + + + +
+
+
 
+ + + + + + + + + + + + + + + + + + + +
zRange(Min):mm
zRange(Max):mm
rRange:mm
+
 
+
+
+ +

+

+

+

+ +
+ +

θCM Calculator

+ +The calculation only give θCM after the bending. +

+ + + + + + + + + + + +
Ex [MeV] : θCM Gate [deg] : X Gate [%] :
+ + + + + + + + + + +
IDpos0(gated)pos1(gated)θ1[deg]θ2[deg]θavg[deg]Δθ[deg]sin(θavg)Δθ
+ + + + + + + + + +
+

+HELIOSmatics was first built by Ben P. Kay in MS Excel around 2010. It was modified by Ryan Tang later. And now it migrated to the web on Dec, 2022. +
+The calculation can be found in the source code (heliosmatics.js or press F12) + +

+ + + + + + + + + + + + diff --git a/WebSimHelper/heliosmatics.js b/WebSimHelper/heliosmatics.js new file mode 100644 index 0000000..2783ebc --- /dev/null +++ b/WebSimHelper/heliosmatics.js @@ -0,0 +1,1159 @@ + +function addRow() { + let table = document.getElementById("ExTable"); + let nRow = table.rows.length; + let row = table.insertRow(nRow-1); + + let energy = Math.random()*Math.min(heavy[4], heavy[5], heavy[6]); + let angle = Math.floor(Math.random() * 30) + 10; + + row.innerHTML = ' \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + '; + + CalculateEZ(); + CalculateRZ(); + + //let table2 = document.getElementById("DWBATable"); + //row = table2.insertRow(nRow-1); + //row.innerHTML = '' + (nRow-2) + ' \ + // \ + // '; +} + +function deleteRow(){ + let table = document.getElementById("ExTable"); + let nRow = table.rows.length; + if ( nRow > 3){ + table.deleteRow(nRow-2); + } + CalculateEZ(); + CalculateRZ(); + + //let table2 = document.getElementById("DWBATable"); + //if ( nRow > 3){ + // table2.deleteRow(nRow-2); + //} +} + +let beam = []; //A, Z, Mass, Name, Sn, Sp, Sa +let beamMass; +let yield; +let target= []; +let light=[]; +let heavy=[]; //A, Z, Mass, Name, Sn, Sp, Sa + +let beamEx; +let BField; +let KEA; +let KE; +let reactionName; + +let Qvalue; +let minKEA; + +let perpDistant = 11.5; //mm, detector prepdicular distance +let bore = 462.0; // mm +let arrayLen = 50 * 10 + 2 * 9 // SOLARIS +let detLen = 50; // SOLARIS +let nDet = 10; // SOLARIS +let detGap = 2; // SOLARIS +let pos = []; //1D array, store the 1st pos of detector +let arrayPos = []; // 2D array, store {tip, toe} of each detector + +const c = 299.792468; // mm/ns + +let beam_k_lab; +let E_tot_cm ; +let KE_cm ; +let max_Ex ; +let beta ; +let gamma ; + +let ez_slope ; // MeV/mm + +let alpha ; +let alpha_B ; + +let xList =[]; // for E-Z plot +let yList =[]; // for E-Z plot +let ExList = []; +let Zb0List = []; +let rhoList = []; + +let ZB0List = []; +let rhoBList = []; + +let xRange ; +let yRange = [0, 12]; + +let zList = []; // for R-Z plot +let rbList = []; // for R-Z plot +let rBList = []; // for R-Z plot + +let color = ['rgb(31,119,180)', // muted blue + 'rgb(255,127,14)', // safety orange + 'rgb(44,160,44)', // cooked asparagus green + 'rgb(214,39,40)', // brick red + 'rgb(148,103,189)', // muted purple + 'rgb(140,86,75)', // chestnut brown + 'rgb(277,119,194)', // raspberry yogurt pink + 'rgb(127,127,127)', // middle gray + 'rgb(118,189,34)', // curry yellow-green + 'rgb(23,190,207)']; // blue-teal + +function GetYield(A,Z){ +// let str = 'beamRate.py?A=' + A + '&Z=' + Z; +// let client = new XMLHttpRequest(); +// client.onreadystatechange = function() { +// let haha = client.responseText.split(","); +// yield = haha[0] +// } +// client.open('GET', str, false); +// client.send(); +// FRIB blocking the request....:( +} + +function GetMassFromSym(AZ, id){ + let str = 'massProxy.py?ASym=' + AZ; + + let client = new XMLHttpRequest(); + client.onreadystatechange = function() { + let mass = client.responseText.split(","); + if( id == 0 ){ + beam[0] = parseInt(mass[0]); + beam[1] = parseInt(mass[1]); + beam[2] = parseFloat(mass[2]); + beam[3] = AZ; + beam[4] = parseFloat(mass[4]); + beam[5] = parseFloat(mass[5]); + beam[6] = parseFloat(mass[6]); + GetYield(beam[0], beam[1]); + } + if( id == 1 ){ + target[0] = parseInt(mass[0]); + target[1] = parseInt(mass[1]); + target[2] = parseFloat(mass[2]); + target[3] = AZ; + } + if( id == 2 ){ + light[0] = parseInt(mass[0]); + light[1] = parseInt(mass[1]); + light[2] = parseFloat(mass[2]); + light[3] = AZ; + } + } + client.open('GET', str, false); + client.send(); + +} + +function GetMassFromAZ(A,Z){ + let str = 'massProxy.py?A=' + A + '&Z=' + Z; + + let client = new XMLHttpRequest(); + client.onreadystatechange = function() { + let mass = client.responseText.split(","); + heavy[2] = parseFloat(mass[2]); + heavy[3] = mass[3]?.trim(); + heavy[4] = parseFloat(mass[4]); + heavy[5] = parseFloat(mass[5]); + heavy[6] = parseFloat(mass[6]); + } + client.open('GET', str, false); + client.send(); + +} + +function GetMass(){ + + GetMassFromSym(document.getElementById('beam').value, 0); + GetMassFromSym(document.getElementById('target').value, 1); + GetMassFromSym(document.getElementById('light').value, 2); + + // console.log(beam); + // console.log(target); + // console.log(light); + + beamMass = beam[2]; + + heavy[0] = beam[0]+target[0]-light[0]; + heavy[1] = beam[1]+target[1]-light[1]; + GetMassFromAZ(heavy[0], heavy[1]); + + // console.log( heavy); + + document.getElementById('heavyName').innerHTML = heavy[3]; + document.getElementById('heavySp').innerHTML = 'Sn: ' + heavy[4] + ' MeV, Sp: ' + heavy[5] + ' MeV, Sa : ' + heavy[6] + ' MeV'; + + //document.getElementById('beamSp').innerHTML = "haah"; + //document.getElementById('beamYield').innerHTML = "FRIB ultimate yield : " + yield + " pps"; + //document.getElementById('n0').innerHTML = beam[0] + "," + beam[1] + "," + beam[2] + //document.getElementById('n1').innerHTML = target[0] + "," + target[1] + "," + target[2] + //document.getElementById('n2').innerHTML = light[0] + "," + light[1] + "," + light[2] + //document.getElementById('n3').innerHTML = heavy[0] + "," + heavy[1] + "," + heavy[2] + +} + +function CalConstants(){ + + beamEx = parseFloat(document.getElementById('beamEx').value); + + beam[2] = beamMass + beamEx; + + BField = parseFloat(document.getElementById('BField').value); + KEA = document.getElementById('KEA').value; + KE = KEA * beam[0]; + + reactionName = beam[3] +"(" + target[3] + "," + light[3] + ")" + heavy[3] + "@" + KEA + "MeV/u, " + BField.toFixed(2) + " T"; + + Qvalue = - heavy[2] - light[2] + beam[2] + target[2] ; + minKEA = (Math.pow(light[2] + heavy[2],2) - Math.pow(beam[2] + target[2],2))/2/target[2]/beam[0]; + + document.getElementById('reactionName').innerHTML = reactionName; + document.getElementById('minKEA').innerHTML = "min Beam Energy: " + minKEA.toFixed(3) + " MeV/u"; + document.getElementById('Q-value').innerHTML = Qvalue.toFixed(3); + + beam_k_lab = Math.sqrt(Math.pow(beam[2] + KE,2) - Math.pow(beam[2],2)); + E_tot_cm = Math.sqrt(Math.pow(target[2] + beam[2],2) + 2*target[2]*KE); + KE_cm = E_tot_cm - beam[2] - target[2]; + max_Ex = KE_cm + Qvalue; + beta = beam_k_lab/(beam[2] + target[2] + KE); + gamma = 1./Math.sqrt(1-beta*beta); + + ez_slope = Math.abs(BField) * c * light[1]*beta/2/Math.PI/1000; // MeV/mm + alpha = ez_slope/beta; + //alpha_B = alpha * heavy[1]/light[1]; + alpha_B = Math.abs(BField) * c * heavy[1]*beta/2/Math.PI/1000 / beta; + +} + +function CalArrayPos(){ + let helios = document.getElementById('HELIOS').checked; + + let haha = parseFloat(document.getElementById('posArray').value); + arrayPos = []; + for( let i = 0; i < nDet; i++){ + let kaka = []; + if( helios == true ){ + if( haha < 0 ){ + kaka.push(haha - pos[i] - detLen ); + kaka.push(haha - pos[i]); + }else{ + kaka.push(haha + pos[i] ); + kaka.push(haha + pos[i] + detLen ); + } + }else{ + if( haha < 0 ){ + kaka.push(haha - (i+1) * detLen - i * detGap); + kaka.push(haha - (i) * detLen - i * detGap); + }else{ + kaka.push(haha + (i) * detLen + i * detGap); + kaka.push(haha + (i+1) * detLen + i * detGap); + } + } + arrayPos.push(kaka); + } +} + +function SetSSType(){ + let helios = document.getElementById('HELIOS').checked; + let solaris = document.getElementById('SOLARIS').checked; + let iss = document.getElementById('ISS').checked; + + if ( helios == true ) { + perpDistant = 11.5; + detGap = 5; + detLen = 50; + nDet = 6; + pos = []; + pos.push(.0); + pos.push(58.6); + pos.push(117.9); + pos.push(176.8); + pos.push(235.8); + pos.push(294.0); + } + if ( solaris == true ) { + perpDistant = 11.5; + detGap = 5; + detLen = 50; + nDet = 10; + } + if ( iss == true ) { + perpDistant = 28.75; + detGap = 0.5; + detLen = 125; + nDet = 4; + } + arrayLen = detLen * nDet + detGap * (nDet-1); + + if( helios ) arrayLen = pos[5] + detLen; + + CalArrayPos(); + //document.getElementById('n1').innerHTML = perpDistant; +} + +function CalculateEZ(){ + + let tableEx = document.getElementById("ExTable"); + let nRow = tableEx.rows.length; + + xList = []; + yList = []; + ExList = []; + Zb0List = []; + rhoList = []; + ZB0List = []; + rhoBList = []; + + //alert("CalculateEZ called, nRow = " + nRow); + + for( let i = 1; i < nRow-1; i++){ + let Ex = parseFloat(document.getElementById("Ex" + i).value); + let theta = parseFloat(document.getElementById("theta" + i).value); + + ExList.push(Ex); + //alert( i, ", Ex : " + Ex); + + let haha1 = E_tot_cm*E_tot_cm - Math.pow(heavy[2] + Ex + light[2],2); + let haha2 = E_tot_cm*E_tot_cm - Math.pow(heavy[2] + Ex - light[2],2); + let k_cm = Math.sqrt(haha1*haha2)/2/E_tot_cm; + + let cs = Math.cos(theta*Math.PI/180.); + let ss = Math.sin(theta*Math.PI/180.); + let qb = Math.sqrt(light[2]*light[2]+k_cm*k_cm); + + let Eb = gamma * qb - gamma * beta * k_cm * cs - light[2]; + + let Zb0 = (gamma*beta* qb - gamma * k_cm * cs )/alpha; //mm + Zb0List.push(Zb0); + + let rho = k_cm * ss/c/light[1]/Math.abs(BField) * 1000; // mm + rhoList.push(rho); + let Zb = Zb0 * (1- Math.asin(perpDistant/rho/2)/Math.PI); + + let thetaLab = 180 - Math.atan2(k_cm * ss, gamma * (beta * qb - k_cm * cs)) * 180/Math.PI; + + let Tcyc = (light[2] + Eb)/alpha/c; + + tableEx.rows[i].cells[2].innerHTML = Eb.toFixed(3); + tableEx.rows[i].cells[3].innerHTML = Zb0.toFixed(1); + tableEx.rows[i].cells[4].innerHTML = Zb.toFixed(1); + tableEx.rows[i].cells[5].innerHTML = (2*rho).toFixed(1); + tableEx.rows[i].cells[6].innerHTML = thetaLab.toFixed(2); + tableEx.rows[i].cells[7].innerHTML = Tcyc.toFixed(2); + + let qB = Math.sqrt(heavy[2]*heavy[2]+k_cm*k_cm); + let EB = gamma * qB + gamma * beta * k_cm * cs - heavy[2]; + let thetaLab_B = Math.atan2(-k_cm * ss, gamma * (beta * qB + k_cm * cs)) * 180/Math.PI; + let ZB0 = (gamma*beta* qB + gamma * k_cm * cs )/alpha_B; //mm + let rhoB = k_cm * ss/c/heavy[1]/Math.abs(BField) * 1000; // mm + ZB0List.push(ZB0); + rhoBList.push(rhoB); + + tableEx.rows[i].cells[8].innerHTML = EB.toFixed(3); + tableEx.rows[i].cells[9].innerHTML = thetaLab_B.toFixed(2); + tableEx.rows[i].cells[10].innerHTML = (ZB0/2).toFixed(2); + tableEx.rows[i].cells[11].innerHTML = (2*rhoB).toFixed(2); + + let xxx = []; + let yyy = []; + + for( let j = 0; j < 100 ; j++){ + let thetaCM = Math.PI/Math.log10(100) * Math.log10(100/(100-j)) ; + let temp = Math.PI * ez_slope / beta / k_cm * perpDistant / Math.sin(thetaCM); // perpDistant / 2/ rho(thetaCM) + if( !isFinite(temp) ) continue; + let pxTemp = 1. /alpha * (gamma * beta * qb - gamma * k_cm * Math.cos(thetaCM)) * (1 - Math.asin(temp)/Math.PI) ; + let pyTemp = gamma * qb - light[2] - gamma * beta * k_cm * Math.cos(thetaCM); + if( !isNaN(pxTemp) && !isNaN(pyTemp) ){ + xxx.push(pxTemp); + yyy.push(pyTemp); + } + }; + + xList.push(xxx); + yList.push(yyy); + + }; + + PlotEZ(); + AdjustRangeEZ(); + +} + +function PlotEZ(){ + + SetSSType(); + Plotly.purge("Plot_EZ"); + + let nEx = xList.length; + let data = []; + + for( let i = 0; i < nEx; i++){ + let kaka = color[i%10]; + let newData = { + x : xList[i], + y : yList[i], + mode:"lines", + type:"scatter", + name:"Ex="+ExList[i], + marker : { color : kaka} + } + data.push(newData); + } + + let haha = parseFloat(document.getElementById('posArrayRange').value); + let xStart = (haha < 0 ? haha - arrayLen - 100 : haha - 100); + let xEnd = (haha < 0 ? haha + 100: haha + arrayLen + 100); + + xRange = [xStart, xEnd]; + //document.getElementById('n0').innerHTML = xRange; + + let haha2 = parseFloat(document.getElementById('eRangeSlider').value); + yRange = [0, haha2]; + + let layout = { + xaxis: {range: xRange, title: { text : "Z [mm]", standoff : 1}, mirror : "allticks", linewidth : "1"}, + yaxis: {range: yRange, title: "Energy [MeV]" , mirror : "allticks", linewidth : "1"}, + title: reactionName, + dragmode : "pan", + margin: { l: 40, r: 40, b : 60, t : 40}, + legend: {yanchor:"top", xanchor:"left", x:"0.01",y:"0.99" } + }; + + + Plotly.newPlot( "Plot_EZ", data, layout, {responsive: true}); + + PlotThetaCMLine(document.getElementById('thetaCM').value); + PlotBore(); + + PlotRZ(); +} + +function PlotThetaCMLine(thetaCM){ + let cs = Math.cos(thetaCM * Math.PI /180); + let ss = Math.sin(thetaCM * Math.PI /180); + + let zzList = []; + let eList = []; + + for( let z = -2000; z < 2000; z+=30){ + zzList.push(z); + + let kaka = Math.pow(gamma * ss,2); + let a1 = light[2]*light[2]*(1-kaka); + let haha = (cs*Math.sqrt(alpha*alpha*z*z + a1) - kaka * ez_slope * z)/(1-kaka)- light[2]; + eList.push(haha); + } + + Plotly.addTraces("Plot_EZ", {x : zzList, + y: eList, + name:"thetaCM=" + thetaCM , + marker : { color : 'rgb(100,100,100)'}, + line : {dash : 'solid', width : 1 } + } + , 0); +} + +function PlotBore(){ + let zzList = []; + let eList = []; + + for( let z = -2000; z < 2000; z+=30){ + zzList.push(z); + + let haha = Math.sqrt((z*z+Math.PI*Math.PI*bore*bore)*alpha*alpha + light[2]*light[2]) - light[2]; + eList.push(haha); + } + + Plotly.addTraces("Plot_EZ", {x : zzList, y: eList, name:"Bore", marker : { color : 'rgb(200,200, 200)'} }, 0); +} + +function CalculateRZ(){ + // this rquire ZB0List and rhoBList from CalculateEZ(); + + zList = []; + rbList = []; + rBList = []; + + for( let z = -2000; z < 4000; z += 2 ) zList.push(z); + + // cal the heavy recoil first + for( let i = 0; i < ZB0List.length; i++){ + let rrr = [] + for( let j = 0; j < zList.length; j++){ + if( zList[j] < 0 ){ + rrr.push(NaN); + }else{ + rrr.push(2*rhoBList[i] *Math.abs( Math.sin(zList[j] * Math.PI / ZB0List[i]))); + } + } + rBList.push(rrr); + } + + // cal the light recoil first + for( let i = 0; i < Zb0List.length; i++){ + let rrr = [] + for( let j = 0; j < zList.length; j++){ + if( target[0] > light[0] ){ + if( zList[j] > 0 ){ + rrr.push(NaN); + }else{ + rrr.push(2*rhoList[i] *Math.abs( Math.sin(zList[j] * Math.PI / Zb0List[i]))); + } + }else{ + if( zList[j] < 0 ){ + rrr.push(NaN); + }else{ + rrr.push(2*rhoList[i] *Math.abs( Math.sin(zList[j] * Math.PI / Zb0List[i]))); + } + } + } + rbList.push(rrr); + } + + PlotRZ(); + AdjustRecoilPos(); +} + +function PlotRZ(){ + Plotly.purge("Plot_RZ"); + + let data = []; + let nEx = ExList.length; + + for(let i = 0 ; i < nEx; i++ ){ + let kaka = color[i%10]; + let newData = { + x : zList, + y : rBList[i], + mode : "lines", + type : "scatter", + name : "Ex="+ExList[i] + ",theta=" + document.getElementById('theta'+(i+1)).value, + marker : { color : kaka} + } + data.push(newData); + } + for(let i = 0 ; i < nEx; i++ ){ + let kaka = color[i%10]; + let newData = { + x : zList, + y : rbList[i], + mode : "lines", + line : {dash : 'dashdot', width : 1 }, + type : "scatter", + name : "Ex="+ExList[i] + ",theta=" + document.getElementById('theta'+(i+1)).value, + marker : { color : kaka} + } + data.push(newData); + } + + let xxx = [parseInt(document.getElementById('zRange1').value), parseInt(document.getElementById('zRange2').value)]; + let yyy = [0, parseInt(document.getElementById('rRange').value)]; + + let layout = { + xaxis: {range: xxx, title: { text : "Z [mm]", standoff : 1}, mirror : "allticks", linewidth : "1"}, + yaxis: {range: yyy, title: "R [mm]" , mirror : "allticks", linewidth : "1"}, + title: reactionName, + showlegend : true, + dragmode : "pan", + margin: { l: 40, r: 40, b : 60, t : 40}, + legend: { yanchor:"top", xanchor:"left", x:"0.01",y:"0.99"} + }; + + Plotly.newPlot( "Plot_RZ", data, layout, {responsive: true}); + +} + +function CalZ(theta, k_cm, qb){ + let cs = Math.cos(theta*Math.PI/180.); + let ss = Math.sin(theta*Math.PI/180.); + + let Zb0 = (gamma*beta* qb - gamma * k_cm * cs )/alpha; //mm + let rho = k_cm * ss / c/ light[1]/ Math.abs(BField) * 1000; // mm + + return Zb0 * ( 1 - Math.sin( perpDistant/2/rho ) / Math.PI ); +} + +function SearchThetaCM(Z, thetaMin, k_cm, qb){ + + let step = 0.01; + + for( let yy = thetaMin; yy < 180; yy += step){ + if( yy > 12 ) step = 0.1; + if( yy > 20 ) step = 0.5; + if( yy > 40 ) step = 1.0; + let z1 = CalZ(yy, k_cm, qb); + let z2 = CalZ(yy+step, k_cm, qb); + //console.log(Z, ", ", yy.toFixed(3), ", ", z1.toFixed(3), ", ", z2.toFixed(3), ", ", (z2-z1).toFixed(3) ); + if( z1 > z2 ){ /// in the bending range + continue; + }else{ + if( Z < z1 ){ + // return yy; + break; + } + if( z1<= Z && Z <= z2 ){ // do a linear approximation + return (Z-z1)/(z2-z1)*step + yy; + } + } + } + return NaN; +} + +function CalThetaCM(){ + let ex = parseFloat(document.getElementById('Ex0').value); + let angGate = parseFloat(document.getElementById('thetaCMGate').value); + let xGate = parseFloat(document.getElementById('XGate').value); + + SetSSType(); + + let haha1 = E_tot_cm*E_tot_cm - Math.pow(heavy[2] + ex + light[2],2); + let haha2 = E_tot_cm*E_tot_cm - Math.pow(heavy[2] + ex - light[2],2); + let k_cm = Math.sqrt(haha1*haha2)/2/E_tot_cm; + let qb = Math.sqrt(light[2]*light[2]+k_cm*k_cm); + + let table = document.getElementById('thetaCMTable'); + let nRow = table.rows.length; + if( nRow > 1 ){ + for( let i = nRow-1; i > 0; i -- ){ + table.deleteRow(i); + } + } + + for( let i = 0; i < arrayPos.length; i++){ + let row = table.insertRow(i+1); + row.insertCell().innerHTML = i; + + let p1 = (arrayPos[i][0] + detLen * (100. - xGate)/200.); + let p2 = (arrayPos[i][1] - detLen * (100. - xGate)/200.); + + row.insertCell().innerHTML = arrayPos[i][0].toFixed(1) + "(" + p1.toFixed(1) + ")"; + row.insertCell().innerHTML = arrayPos[i][1].toFixed(1) + "(" + p2.toFixed(1) + ")"; + + ///Search thetaCM for Z + let a1 = SearchThetaCM(p1, angGate, k_cm, qb); + let a2 = SearchThetaCM(p2, angGate, k_cm, qb); + row.insertCell().innerHTML = a1.toFixed(2); + row.insertCell().innerHTML = a2.toFixed(2); + + let am = (a2+a1)/2; + let da = Math.abs(a2-a1); + + row.insertCell().innerHTML = am.toFixed(2); + row.insertCell().innerHTML = da.toFixed(2); + row.insertCell().innerHTML = (Math.sin(am*Math.PI/180) * da * Math.PI/180).toExponential(2); + + } + +} + +document.getElementById('beam').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + GetMass(); + CalConstants(); + } + }, false +); + +document.getElementById('beamEx').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + CalConstants(); + } + }, false +); + +document.getElementById('target').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + GetMass(); + CalConstants(); + } + }, false +); + +document.getElementById('light').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + GetMass(); + CalConstants(); + } + }, false +); + +document.getElementById('BField').addEventListener('keypress', + function(e){ + //document.getElementById('n0').innerHTML = e.keyCode; + if(e.keyCode == 13){ + document.getElementById('BRange').value = this.value; + CalConstants(); + CalculateEZ(); + CalculateRZ(); + CalThetaCM(); + } + }, false +); + +document.getElementById('BRange').oninput = function(){ + document.getElementById('BField').value = this.value; + CalConstants(); + CalculateEZ(); + CalculateRZ(); + AdjustRangeEZ(); + CalThetaCM(); +} + +document.getElementById('KEA').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('KEARange').value = this.value; + CalConstants(); + CalculateEZ(); + CalculateRZ(); + CalThetaCM(); + } + }, false +); + +document.getElementById('KEARange').oninput = function(){ + document.getElementById('KEA').value = this.value; + CalConstants(); + CalculateEZ(); + CalculateRZ(); + CalThetaCM(); +} + +document.getElementById('thetaCM').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('thetaCMRange').value = this.value; + CalConstants(); + CalculateEZ(); + AdjustRecoilPos(); + } + }, false +); + +document.getElementById('thetaCMRange').oninput = function(){ + document.getElementById('thetaCM').value = this.value; + CalConstants(); + CalculateEZ(); + AdjustRecoilPos(); +} + +function AdjustRangeEZ(){ + let haha = parseFloat(document.getElementById('posArray').value); + let xStart = (haha < 0 ? haha - arrayLen - 100 : haha - 100); + let xEnd = (haha < 0 ? haha + 100: haha + arrayLen + 100); + + SetSSType(); + + xRange = [xStart, xEnd]; + yRange = [0, parseFloat(document.getElementById('eRangeSlider').value)]; + + //document.getElementById('n0').innerHTML = pos; + //document.getElementById('n2').innerHTML = arrayLen; + //document.getElementById('n1').innerHTML = xRange; + + let shapeArray = []; + for( let i = 0; i < nDet; i++){ + let newBlock = { + type: 'rect', + xref: 'x', + yref: 'paper', + x0 : arrayPos[i][0], + x1 : arrayPos[i][1], + y0 : 0, + y1 : 1, + fillcolor : '#9999FF', + opacity : 0.1, + line : { width : 0} } + shapeArray.push(newBlock); + } + + let update = { + 'xaxis.range' : xRange, + 'yaxis.range' : yRange, + 'shapes' : shapeArray + } + Plotly.relayout("Plot_EZ", update); +} + +document.getElementById('posArray').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('posArrayRange').value = this.value; + AdjustRangeEZ(); + AdjustRecoilPos(); + CalThetaCM(); + } + }, false +); + +document.getElementById('posArrayRange').oninput = function(){ + document.getElementById('posArray').value = this.value; + AdjustRangeEZ(); + AdjustRecoilPos(); + CalThetaCM(); +} + +function AdjustRecoilPos(){ + let pos1 = parseInt(document.getElementById('posRecoil').value); + let pos2 = pos1+10; + + let radius1 = parseInt(document.getElementById('innerRecoil').value); + let radius2 = parseInt(document.getElementById('outterRecoil').value); + + //document.getElementById('n0').innerHTML = pos1; + + let shapeRecoil = { + type: 'rect', + xref: 'x', + yref: 'y', + x0 : pos1, + x1 : pos2, + y0 : radius1, + y1 : radius2, + fillcolor : '#FF0000', + opacity : 0.4, + line : { width : 0} + }; + + let totalShape = []; + for( let i = 0; i < nDet; i++){ + let newhaha = { + type: 'rect', + xref: 'x', + yref: 'y', + x0 : arrayPos[i][0], + x1 : arrayPos[i][1], + y0 : 0, + y1 : perpDistant, + fillcolor : '#9999FF', + opacity : 0.1, + line : { width : 0} }; + totalShape.push(newhaha); + } + totalShape.push(shapeRecoil); + + let update = { + 'shapes' : totalShape + } + Plotly.relayout("Plot_RZ", update); + +} + +document.getElementById('posRecoil').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('posRecoilRange').value = this.value; + AdjustRecoilPos(); + } + }, false +); + +document.getElementById('posRecoilRange').oninput = function(){ + document.getElementById('posRecoil').value = this.value; + AdjustRecoilPos(); +} + +document.getElementById('innerRecoil').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + AdjustRecoilPos(); + } + }, false +); + +document.getElementById('outterRecoil').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + AdjustRecoilPos(); + } + }, false +); + +function AdjustRangeRZ(){ + + let zmin = parseInt(document.getElementById('zRange1').value); + let zmax = parseInt(document.getElementById('zRange2').value); + let rmax = parseInt(document.getElementById('rRange').value); + + let update = { + 'xaxis.range' : [zmin, zmax], + 'yaxis.range' : [0, rmax] + } + Plotly.relayout("Plot_RZ", update); +} + +document.getElementById('zRange1').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('zRange1Slider').value = this.value; + AdjustRangeRZ(); + } + }, false +); + +document.getElementById('zRange1Slider').oninput = function(){ + document.getElementById('zRange1').value = this.value; + AdjustRangeRZ(); +} + +document.getElementById('zRange2').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + document.getElementById('zRange2Slider').value = this.value; + AdjustRangeRZ(); + } + }, false +); + +document.getElementById('zRange2Slider').oninput = function(){ + document.getElementById('zRange2').value = this.value; + AdjustRangeRZ(); +} + +document.getElementById('rRange').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + let rrrr = parseInt(this.value); + if ( rrrr < 1 ){ + document.getElementById('rRange').value = 1; + rrrr = 1; + } + document.getElementById('rRangeSlider').value = rrrr; + AdjustRangeRZ(); + + } + }, false +); + +document.getElementById('rRangeSlider').oninput = function(){ + document.getElementById('rRange').value = this.value; + AdjustRangeRZ(); +} + +document.getElementById('eRange').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + let rrrr = parseInt(this.value); + if ( rrrr < 1 ){ + document.getElementById('eRange').value = 1; + rrrr = 1; + } + document.getElementById('eRangeSlider').value = rrrr; + AdjustRangeEZ(); + } + }, false +); + +document.getElementById('eRangeSlider').oninput = function(){ + document.getElementById('eRange').value = this.value; + AdjustRangeEZ(); +} + +window.logMeThis = function(){ + SetSSType(); + CalculateEZ(); + CalculateRZ(); +} + +let FuncEx = window.logMeThis.bind(null, "Ex"); +window.addEventListener('keypress', FuncEx); + +let FuncThetaCM = window.logMeThis.bind(null, "thetaCM"); +window.addEventListener('keypress', FuncThetaCM); + +window.checkSSType = function(){ + SetSSType(); + CalculateEZ(); + CalThetaCM(); + CalculateRZ(); +} + +let FuncCheckSSType = window.checkSSType.bind(null, "SSType"); +document.getElementById('HELIOS').addEventListener('click', FuncCheckSSType); +document.getElementById('SOLARIS').addEventListener('click', FuncCheckSSType); +document.getElementById('ISS').addEventListener('click', FuncCheckSSType); + +document.getElementById('Ex0').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + CalThetaCM(); + } + }, false +); + +document.getElementById('thetaCMGate').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + CalThetaCM(); + } + }, false +); + +document.getElementById('XGate').addEventListener('keypress', + function(e){ + if(e.keyCode == 13){ + CalThetaCM(); + } + }, false +); + + +function HalfCylinder(r, y, d, up){ + let a = []; + let b = []; + let c = []; + + let nStep = 60; + let step = 2*r/nStep; + + for( let i = 0; i <= nStep ; i++){ + let x = -r + i * step; + a.push(x); + c.push(Math.sqrt(r*r - x*x) * up); + b.push(y); + } + + + for( let i = 0; i <= nStep; i ++){ + let x = r - i * step; + a.push(x); + c.push(Math.sqrt(r*r - x*x) * up); + b.push(y+d); + } + + + let haha = []; + haha.push(a); + haha.push(b); + haha.push(c); + + return haha; +} + +function detMesh(startPos, phi){ + let aaa = 110; // prepdist + let len = 500; // detLen + let www = 100; //det width + + let a =[]; + let b =[]; + let c =[]; + + let cs = Math.cos(phi * Math.PI / 180); + let ss = Math.sin(phi * Math.PI / 180); + + let A0 = aaa * cs - www / 2 * ss; + let A1 = aaa * ss + www / 2 * cs; + let B0 = aaa * cs + www / 2 * ss; + let B1 = aaa * ss - www / 2 * cs; + + a.push(A0); b.push(A1); c.push(startPos); + a.push(B0); b.push(B1); c.push(startPos); + a.push(B0); b.push(B1); c.push(startPos + len); + a.push(A0); b.push(A1); c.push(startPos + len); + + let haha = []; + haha.push(a); + haha.push(b); + haha.push(c); + + return haha; + +} + +function Helix(theta, phi, rho, sign, nCyc){ + // sign = B-field + // nCyc < 0 = updatream + + let a = []; + let b = []; + let c = []; + + let deg = Math.PI/180; + let nStep = 100; + let ts = Math.tan(theta * deg); + let zRange = nCyc * rho/ts * Math.PI * 2; + let zStep = zRange/nStep; + + + for( let i = 0; i < nStep; i++){ + let zzz; + zzz = i * zStep; + b.push(zzz); + a.push( rho * ( Math.sin( ts * zzz / rho - phi * deg ) + Math.sin(phi * deg) ) ); + c.push( rho * sign * (Math.cos( ts * zzz / rho - phi * deg ) - Math.cos(phi * deg)) ); + } + + let haha = []; + haha.push(a); + haha.push(b); + haha.push(c); + + return haha; + +} + + +function Plot3D(){ + + let u1 = HalfCylinder(bore, -1000, 2500., 1); + let d1 = HalfCylinder(bore, -1000, 2500., -1); + + let line1 = Helix(40, 0, 80, -1, -1); + let line2 = Helix(2, 180, 10, -1, 0.5); + + let data = [ + {type : 'mesh3d', x : u1[0], y : u1[1], z : u1[2], hoverinfo: 'none', opacity : 0.2, color : "#DDDDDD"} + ,{type : 'mesh3d', x : d1[0], y : d1[1], z : d1[2], hoverinfo: 'none', opacity : 0.2, color : "#DDDDDD"} + ,{type : 'scatter3d', mode : 'lines', x : line1[0], y : line1[1], z : line1[2], name : "haha"} + ,{type : 'scatter3d', mode : 'lines', x : line2[0], y : line2[1], z : line2[2], name : "kaka"} + ]; + + + /* + let data = []; + + for( let i = 0; i < 1; i++){ + let haha = detMesh(-200, 360/6*i + 90); + + console.log(haha[0], "| ", haha[1], "| ", haha[2]); + data.push( + {type : 'mesh3d', x : haha[0], y : haha[1], z : haha[2], hoverinfo: 'none', opacity : 1.0, color : "#FF0000"} + ); + }*/ + + let layout = { + aspectmode: 'manual', + aspectratio:{ x : 0.5, y : 0.5, z : 2.5 }, + margin: { l: 40, r: 40, b : 60, t : 40}, + }; + + Plotly.newPlot('Plot_3D', data, layout); + +} + +GetMass(); +CalConstants(); +SetSSType(); +CalculateEZ(); +CalculateRZ(); +AdjustRecoilPos(); +CalThetaCM(); + +//Plot3D(); diff --git a/WebSimHelper/index.html b/WebSimHelper/index.html new file mode 100644 index 0000000..f165f57 --- /dev/null +++ b/WebSimHelper/index.html @@ -0,0 +1,196 @@ + + + + SOLARIS Si-Array Mode Simulation + + + + + + + + + +
+
+
+ +
+
+ Simulation +
+
+
+ +
+ + + + +
+ +
+
+ + + + + + + + diff --git a/WebSimHelper/instruction.html b/WebSimHelper/instruction.html new file mode 100644 index 0000000..a5d3dd3 --- /dev/null +++ b/WebSimHelper/instruction.html @@ -0,0 +1,76 @@ + + + + +

Intructions:

+

+

    +
  • Use the HELIOSmatics to check the detector and field settings.
    • +
  • The easiest way to do DWBA and Monte Carlo Simulation is use the Simplied Interface +
  • The kenimatic calculation is documented in Here.
    • +
  • The DWBA calucation is using Ptolemy. Here for more detail.
    • +
  • Past calculations can be found Here. Clear every Monday.
    • +
+

+ +
+

Credits:

+ +

This page is created and hosted by Fox's Lab (FSU) in collabortion with SOLARIS (FRIB).

+

The SOLARIS project is based on HELIOS (ANL) and is leaded by ANL.

+

The ISS (ISOLDE Solenoidal Spectrometer) is located as CERN.

+ + + + + + + + + + +

+The simulation was started from Ben Kay (ANL) around 2010 using excel spreadsheet. +
+Ryan Tang (former ANL postdoc, now at FSU) developed a Monte Carlo simulation with DWBA using CERN ROOT framework. +
+And the whole simulation migrated to here at Dec, 2022. + +

+Contact: Ryan Tang (rtang at fsu.edu) +

any suggestions and feedback are very welcome.

+ +
+

Todo :

+
    +
  • File name of the images of past calculations include reaction
    • +
  • Change the legend position of the DWBA calcualtion
    • +
  • A 3-D simulation?
    • +
+ + + diff --git a/WebSimHelper/isotopeLib.py b/WebSimHelper/isotopeLib.py new file mode 100755 index 0000000..ed42a40 --- /dev/null +++ b/WebSimHelper/isotopeLib.py @@ -0,0 +1,185 @@ +#!/usr/bin/env /usr/bin/python3 + +import pandas as pd +import urllib.request +import re +import numpy + +me = 0.51099895000 # +- 15 +mp = 938.27208816 # +- 29 +mn = 939.56542052 # +- 54 +ma = 3727.37915 + +livechart = "https://nds.iaea.org/relnsd/v0/data?" + +def lc_read_csv(url): + req = urllib.request.Request(url) + req.add_header('User-Agent', 'Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:77.0) Gecko/20100101 Firefox/77.0') + return pd.read_csv(urllib.request.urlopen(req)) + +# Read the saved CSV file back into a DataFrame +try : + haha = pd.read_csv('IAEA_NuclearData.csv') +except FileNotFoundError: + # the service URL + url = livechart + "fields=ground_states&nuclides=all" + haha = lc_read_csv(url) + haha.insert(0, 'A', haha['z'] + haha['n']) + haha.to_csv('IAEA_NuclearData.csv', index=False) + haha = pd.read_csv('IAEA_NuclearData.csv') + +# print(haha.columns) +## ['A', 'z', 'n', 'symbol', 'radius', 'unc_r', 'abundance', 'unc_a', +## 'energy_shift', 'energy', 'unc_e', 'ripl_shift', 'jp', 'half_life', +## 'operator_hl', 'unc_hl', 'unit_hl', 'half_life_sec', 'unc_hls', +## 'decay_1', 'decay_1_%', 'unc_1', 'decay_2', 'decay_2_%', 'unc_2', +## 'decay_3', 'decay_3_%', 'unc_3', 'isospin', 'magnetic_dipole', 'unc_md', +## 'electric_quadrupole', 'unc_eq', 'qbm', 'unc_qb', 'qbm_n', 'unc_qbmn', +## 'qa', 'unc_qa', 'qec', 'unc_qec', 'sn', 'unc_sn', 'sp', 'unc_sp', +## 'binding', 'unc_ba', 'atomic_mass', 'unc_am', 'massexcess', 'unc_me', +## 'me_systematics', 'discovery', 'ENSDFpublicationcut-off', +## 'ENSDFauthors', 'Extraction_date'] + + +def GetExList(ASym : str, maxEx : float): + try: + exList = lc_read_csv(livechart + "fields=levels&nuclides=" + ASym) + exJpi = exList[['energy', 'jp']] + exJpi = exJpi[exJpi['energy'] < (maxEx * 1000)] + return exJpi + except: + return pd.DataFrame() + +def BreakDownName(ASym : str): + match = re.match(r'(\d+)(\D+)', ASym) + return [int(match.group(1)), match.group(2) ] + +def GetAZ(ASym : str): + [A, sym] = BreakDownName(ASym) + try: + dudu = haha[(haha['symbol']==sym) & (haha['A']==A)] + Z = int(dudu['z'].iloc[0]) + return [A, Z] + except: + return [A, numpy.nan] + +def GetBindingPerA(ASym : str) -> float: + [A, sym] = BreakDownName(ASym) + try: + dudu = haha[(haha['symbol']==sym) & (haha['A']==A)] + Z = int(dudu['z'].iloc[0]) + N = A - Z + return dudu['binding'].iloc[0]/1000 + except: + return numpy.nan + +def GetMassFromSym(ASym : str) -> float: + [A, sym] = BreakDownName(ASym) + try: + dudu = haha[(haha['symbol']==sym) & (haha['A']==A)] + Z = int(dudu['z'].iloc[0]) + N = A - Z + binding = dudu['binding'].iloc[0]/1000 + return Z*mp + N*mn - binding*A + except: + return numpy.nan + +def GetMassFromAZ(A : int, Z : int) -> float: + try: + dudu = haha[(haha['z']==Z) & (haha['A']==A)] + Z = int(dudu['z'].iloc[0]) + N = A - Z + binding = dudu['binding'].iloc[0]/1000 + return Z*mp + N*mn - binding*A + except: + return numpy.nan + +def GetSymbol(A : int, Z : int) -> str: + try: + dudu = haha[(haha['z']==Z) & (haha['A']==A)] + return "%d%s" % (A , dudu['symbol'].iloc[0]) + except: + return "0x" + +def GetJpi(ASym : str): + [A, sym] = BreakDownName(ASym) + try: + dudu = haha[(haha['symbol']==sym) & (haha['A']==A)] + return dudu['jp'].iloc[0] + except: + return "unknown" + +def GetHalfLife(ASym : str) -> float: + [A, sym] = BreakDownName(ASym) + try: + dudu = haha[(haha['symbol']==sym) & (haha['A']==A)] + return dudu['half_life_sec'].iloc[0] + except: + return "unknown" + +def GetSn(ASym : str) -> float: + [A, Z] = GetAZ(ASym) + if numpy.isnan(Z) : + return numpy.nan + else: + mass0 = GetMassFromAZ(A, Z) + mass1 = GetMassFromAZ(A-1, Z) + return mass1 + mn - mass0 + +def GetSp(ASym : str): + [A, Z] = GetAZ(ASym) + if numpy.isnan(Z) : + return numpy.nan + else: + mass0 = GetMassFromAZ(A, Z) + mass1 = GetMassFromAZ(A-1, Z-1) + return mass1 + mp - mass0 + +def GetSa(ASym : str): + [A, Z] = GetAZ(ASym) + if numpy.isnan(Z) : + return numpy.nan + else: + mass0 = GetMassFromAZ(A, Z) + mass1 = GetMassFromAZ(A-4, Z-2) + return mass1 + ma - mass0 + +def PrintIso(ASym: str): + [A, Z] = GetAZ(ASym) + print("========================= ", ASym) + print("A : %d, Z : %d, N : %d" % (A, Z, A-Z)) + print(" Jpi : ", GetJpi(ASym)) + print("half-live : %.2f sec" % (GetHalfLife(ASym))) + print(" Mass : %9.2f MeV" % (GetMassFromSym(ASym) )) + print(" Binding : %9.2f MeV/A" % (GetBindingPerA(ASym))) + print(" Binding : %9.2f MeV" % (GetBindingPerA(ASym) * A)) + print(" Sn: %9.2f MeV" % GetSn(ASym)) + print(" Sp: %9.2f MeV" % GetSp(ASym)) + print(" Sa: %9.2f MeV" % GetSa(ASym)) + print("=============================") + +def PrintIsoWeb(ASym : str): + [A, Z] = GetAZ(ASym) + print("
========================= ", ASym) + print("
A : %d, Z : %d, N : %d" % (A, Z, A-Z)) + print("
Jpi : ", GetJpi(ASym)) + print("
half-live : %.2f sec" % (GetHalfLife(ASym))) + print("
Mass : %9.2f MeV" % (GetMassFromSym(ASym) )) + print("
Binding : %9.2f MeV/A" % (GetBindingPerA(ASym))) + print("
Binding : %9.2f MeV" % (GetBindingPerA(ASym) * A)) + print("
Sn: %9.2f MeV" % GetSn(ASym)) + print("
Sp: %9.2f MeV" % GetSp(ASym)) + print("
Sa: %9.2f MeV" % GetSa(ASym)) + print("
=============================") + + +def PrintIsoExWeb(ASym : str, maxEx : float): + exList = GetExList(ASym, maxEx) + if exList.empty: + print("
cannot find Ex data") + else: + print("") + for i in range(0,len(exList)): + print("" % (exList['energy'].iloc[i]/1000, exList['jp'].iloc[i].replace(' ', ','))) + print("
%9.3f%s
") + diff --git a/WebSimHelper/logos/ANL_logo.gif b/WebSimHelper/logos/ANL_logo.gif new file mode 100644 index 0000000..e757086 Binary files /dev/null and b/WebSimHelper/logos/ANL_logo.gif differ diff --git a/WebSimHelper/logos/CERN_logo.svg b/WebSimHelper/logos/CERN_logo.svg new file mode 100644 index 0000000..02a64f0 --- /dev/null +++ b/WebSimHelper/logos/CERN_logo.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/WebSimHelper/logos/FRIB_logo.jpg b/WebSimHelper/logos/FRIB_logo.jpg new file mode 100644 index 0000000..d34be10 Binary files /dev/null and b/WebSimHelper/logos/FRIB_logo.jpg differ diff --git a/WebSimHelper/logos/FSU_logo_640.png b/WebSimHelper/logos/FSU_logo_640.png new file mode 100644 index 0000000..5778534 Binary files /dev/null and b/WebSimHelper/logos/FSU_logo_640.png differ diff --git a/WebSimHelper/logos/HELIOS_logo.jpg b/WebSimHelper/logos/HELIOS_logo.jpg new file mode 100644 index 0000000..6dcc11a Binary files /dev/null and b/WebSimHelper/logos/HELIOS_logo.jpg differ diff --git a/WebSimHelper/logos/ISS_logo.png b/WebSimHelper/logos/ISS_logo.png new file mode 100644 index 0000000..2dad6e8 Binary files /dev/null and b/WebSimHelper/logos/ISS_logo.png differ diff --git a/WebSimHelper/logos/SOLARIS_favicon.png b/WebSimHelper/logos/SOLARIS_favicon.png new file mode 100644 index 0000000..2f976e3 Binary files /dev/null and b/WebSimHelper/logos/SOLARIS_favicon.png differ diff --git a/WebSimHelper/logos/SOLARIS_logo.png b/WebSimHelper/logos/SOLARIS_logo.png new file mode 100644 index 0000000..bed86a6 Binary files /dev/null and b/WebSimHelper/logos/SOLARIS_logo.png differ diff --git a/WebSimHelper/massProxy.py b/WebSimHelper/massProxy.py new file mode 100755 index 0000000..c8ec507 --- /dev/null +++ b/WebSimHelper/massProxy.py @@ -0,0 +1,43 @@ +#!/usr/bin/python3 + +import isotopeLib +import cgi + +form = cgi.FieldStorage() + +ASym = form.getvalue('ASym') +A=form.getvalue('A') +Z=form.getvalue('Z') + + +print( "Content-type:text/html\r\n\r\n") + +if 'ASym' in form : + + if ASym == "p" : ASym = "1H" + if ASym == "d" : ASym = "2H" + if ASym == "t" : ASym = "3H" + if ASym == "a" : ASym = "4He" + + [A,Z] = isotopeLib.GetAZ(ASym) + +if 'A' in form: + ASym = isotopeLib.GetSymbol(float(A), float(Z)) + +#=================================== +# print(A) +# print(Z) +# print(isotopeLib.GetMassFromSym(ASym)) +# print(ASym) +# print(isotopeLib.GetSn(ASym)) +# print(isotopeLib.GetSp(ASym)) +# print(isotopeLib.GetSa(ASym)) + +print(A , ",", + Z , ",", + format(isotopeLib.GetMassFromSym(ASym), '.2f'), ",", + ASym, ",", + format(isotopeLib.GetSn(ASym), '.2f'), ",", + format(isotopeLib.GetSp(ASym), '.2f'), ",", + format(isotopeLib.GetSa(ASym), '.2f') + ) diff --git a/WebSimHelper/miscCal.html b/WebSimHelper/miscCal.html new file mode 100644 index 0000000..0bf4445 --- /dev/null +++ b/WebSimHelper/miscCal.html @@ -0,0 +1,287 @@ + + + + + + + + +

enA, pnA, pps convertor

+ + + + + + + + + + + + + + + + + + + + + + +
Charge state:
Attenuation:
enA pnA pps att. pps
+ +

Yield Calculator

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Integrated Xsec mb
Beam intensity pps
Target thickness ug/cm2
Target molar mass g/mol
Nucleus/molecule
Num of nucleus per area count/cm2
Yield per sec count/sec
Spectroscopic factor
Wanted Count
stat. uncertainty %
Beam Time required day
+ + + + + + + diff --git a/WebSimHelper/simpleSim.html b/WebSimHelper/simpleSim.html new file mode 100644 index 0000000..1995901 --- /dev/null +++ b/WebSimHelper/simpleSim.html @@ -0,0 +1,529 @@ + + + + DWBA and Monte Carlo Simulation + + + + + +

DWBA and Monte Carlo Simulation

+ + + +
+ +

Reaction:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + +
Beam EnergyMeV/u
Beam Jπ + Ex:MeV
Target
Light recoil
Number of events
+ +

Detector:

+ + + + HELIOS +
+ SOLARIS +
+ ISS + + + + + + + + + + + + + + + + +
B-field T (minus sign = field point to upstream)
Array Pos.mm (negative for upstream)
Recoil Pos.mm (negative for upstream)
+ +

DWBA and Ex List:

+ +

For 2-nucleon transfer, Orbital take the form NL=X, where N is number of node, X is momentum number. n and L are related by Σi (2ni+li) = 2N + X + 2n + l, where ni and li are the principle number and orbital angular momentum of the each transfered nucleon, and n and l are the internal quanta of the 2-nucleon. e.g. (t,p) reaction to 0f orbtial, the left-hand side would be ni = 0 and li = 3 and the sum is 3+3 = 6 = 2N + X + 2n+l. Assume n = l = 0, we have 6 = 2N+L. Thus, 3L=0, 2L=2,1L=4, 0L=6.

+ + TODO: guess the orbital +
+ Positive parity + Negative parity + Unknown parity +
+ + Isotope: + Max Ex: MeV + +

+ + + + + + + + + + + + + + + + + + + + + +
Ex [MeV]    Jπ Orbital
+ +

+ Cal. DWBA + + + + + + + + + +
Incoming Channel + + +
Outgoing Channel + + +
+ + + + Only DWBA and Don't Sim. Angle range (for only DWBA) + + + + + + + +
Min [deg]: Max [deg]:
+ +

Plot config:

+ + + + + + + + + + + + + + + + + +
E vs ZEx (cal.)ThetaCM
ThetaCM vs ZRecoil X vs YRecoil-R vs ThetaCM
Recoil R vs ZTime diff vs ZArray X vs Y
+ +

+ + Array Hit
+ Loop = 1
+ ThetaCM > 10 deg
+ +

+ + +
+ +
+ +

Advanced control

+ + + + + + + + + + + + + + + + + + + + + + + + + + + +
Download Sample files: + Reaction File +
+ DetectorGeo (SOLARIS) File +
+ DetectorGeo (HELIOS) File +
+ Ex File +
+ DWBA File +
+ Plot Config File +
+ +

+
+ + + + + + + + + + + + + + + + + + + + + + + +
Reaction File
DetectorGeo File +
Ex File +
DWBA File ^ +
*.in File ^ +
Plot Config File # +
^ can be alone
# can be omitted +
+
+ +
    +
  • File name can be customized.
    • +
  • For kinematic simulation, only the reactionConfig.txt, detectorGeo.txt, and Ex.txt are needed.
    • +
  • For DWBA calculation, only the DWBA file is needed.
    • +
  • If reactionConfig.txt, detectorGeo.txt, and DWBA file are presented, will do DWBA and use the result for simulation.
    • +
  • When the DWBA file is presented, the kinematic simulation will use the DWBA result and also the excitation energy. i.e. the user provide Ex.txt will not be used.
    • +
  • User can use a customs in File for DWBA calculation. Once the in File exist, it ignores the DWBA file.
    • +
  • To change DWBA angular range, download the in file, edit it. But becareful, DWBA for kinematic simulation must be 0 -180 deg.
    • +
+ +
+The source code for calculation can be found in Here + + + + + + +