Changeset 7 in ZHANGProjects

Timestamp:

Mar 25, 2014, 11:22:11 AM (10 years ago)

Author:

zhangj

Message:

Fix the bugs to read the collimator files and optics files.

Location:

ICOSIM/CPP/trunk

Files:

• . (modified) (1 prop)
• sample (modified) (1 prop)
• sample/collimatorfile_standard_ions_nominal.csv (added)
• sample/opticsfile_standard_ions_nominal.csv (added)
• source (modified) (1 prop)
• source/.kdev4/source.kdev4 (modified) (1 diff)
• source/Particle.cc (modified) (2 diffs)
• source/StandardCollimator.cc (modified) (3 diffs)
• source/lattice.cc (modified) (4 diffs)
• source/simulation.cc (modified) (20 diffs)

ICOSIM/CPP/trunk/source/.kdev4/source.kdev4

@@ -1 +1 @@
 [Buildset]
-BuildItems=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x00\x01\x00\x00\x00\x10\x00I\x00C\x00O\x00S\x00Y\x00M\x00+\x00+)
+BuildItems=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x00)
+
+[Launch]
+Launch Configurations=Launch Configuration 0
+
+[Launch][Launch Configuration 0]
+Configured Launch Modes=execute
+Configured Launchers=nativeAppLauncher
+Name=New Native Application Configuration
+Type=Native Application
+
+[Launch][Launch Configuration 0][Data]
+Arguments=../sample/collimatorfile_standard_ions_nominal.csv   ../outputs
+Dependencies=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x00)
+Dependency Action=Nothing
+EnvironmentGroup=default
+Executable=icosim++
+External Terminal=konsole --noclose --workdir %workdir -e %exe
+Project Target=
+Use External Terminal=false
+Working Directory=file:///home/jfz/codes/ICOSIMPP/CPP/source
+isExecutable=true
+
+[Project]
+VersionControlSupport=kdevsubversion

ICOSIM/CPP/trunk/source/Particle.cc

@@ -9 +9 @@
 //==================================Constructors, destructor====================================
 
+/*
+   Initialize the 6-D particle coordinates.
+ */
 Particle::Particle(double x, double dx, double y, double dy, double deltaE, double t, double mass, double charge, double moment)
     : Ap0(mass), Zp0(charge), dpoporiginal(moment)
@@ -123 +126 @@
 }
 
-
+/*
+   set the ID of the registed particle ID.
+ */
 void Particle::setidentification(int id)
 {

ICOSIM/CPP/trunk/source/StandardCollimator.cc

@@ -25 +25 @@
 {
 
-    long double ca, sa;
+    long double ca=0.0, sa=0.0;
 
     sa = sin(this->tcang);
@@ -34 +34 @@
     p1.coordonnees[0][2] = p1.coordonnees[0][2] - scaleorbit * this->YC;
 
-    long double xl, xsl, alf, lcolleff;
+    long double xl=0.0, xsl=0.0, alf=0.0, lcolleff=0.0;
 
 
@@ -109 +109 @@
 
         //Apply thin collimator at impact position
-
+          /* Initialize the 6-D particle coordinates*/
         Particle ptemp(xint, xsint, yint, ysint, p1.coordonnees[0][4], 0, p1.Ap0, p1.Zp0, p1.coordonnees[0][4]);
 

ICOSIM/CPP/trunk/source/lattice.cc

@@ -269 +269 @@
 
 
-
+/*
+ * 
+ */
 void Lattice::trackensemblelinearnew(vector <Particle>& bunch, vector <Particle>& bunchhit, const int& nrev, const double& blowup2, const int& blowupperiod)
 {
 
-    int nip, niph, nrevhitp;
+  /*nip: number of primary collimators in the accelerator; niph: number of primary collimators hit by the particles; nrevhitp: number of turns when the particle hit the primary collimators */
+    int nip=0, niph=0, nrevhitp=0; 
     vector <int> iph;
-    double blowup;
+    double blowup=0.0; /* sqrt(blowup2), beam blow up strength*/
 
     //we only take the primary collimators into account here, as well as the first and the last elements of the lattice
@@ -331 +334 @@
                 if (bunch[p].in == 1) { //to assure the particle is still remainding in the accelerator
 
-                    long double pdepth, pdepth2;
+                    long double pdepth=0.0, pdepth2=0.0;
 
                     bunch[p].coordonnees[1][0] = R11X[i] * bunch[p].coordonnees[0][0] + R12X[i] * bunch[p].coordonnees[0][1] + (reseau[iph[i + 1]]->DX - R11X[i] * reseau[iph[i]]->DX - R12X[i] * reseau[iph[i]]->DPX) * bunch[p].coordonnees[0][4];
@@ -339 +342 @@
 
                     if (i < niph - 1) {
-                        long double lcoll, sa, ca;
+                        long double lcoll=0.0, sa=0.0, ca=0.0;
 
                         sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle
@@ -416 +419 @@
             if ((k + 1) % blowupperiod == 0) { //blowup/diffusion
                 if (bunch[p].in == 1) {
-                    //cout << "Blow-up!!" << k+1 << endl;
+                  //  cout << "Blow-up!!" << k+1 << endl;
                     int im;
                     im = reseau.size() - 1;

ICOSIM/CPP/trunk/source/simulation.cc

@@ -35 +35 @@
 
     for (int l(0); l < sizeBunch; ++l) {
-
+        /* generate the particle distributions */
         p1.genpartdist(1, 1, genere, 1.00095, 1.2e-8, 1.2e-8, 0.00011, bx, ax, dx, dpx, by, ay, dy, dpy, 6);
         p1.afficheFirstCoordonnees();
@@ -67 +67 @@
 
 
+/*
+ void Simulation::readParticle(const double& Apr, const double& Zpr, string inputfile)
+ 
+ */
 void Simulation::readParticle(const double& Apr, const double& Zpr, string inputfile)
 {
@@ -536 +540 @@
 
 
-    double Apr, Zpr, mass, energyPerIon, emix, emiy, sigdpp, r1r2skin, blowup2, thicknessmagneticfield, Bmax, deltaGap, gamma, betgam, W, PlossPb, freqrf, rfvoltage, rfharmonic, wecolli, nsigi, momentum;
-    int npart, kbunch, taubeam, nparti, nrev1, nrev2, blowupperiod, scaleorbit, beamflag, stopLinearTrackingNo, RunWithFluka, RunWithCrystal, RunningFlag, outcoord, idpart, idelt, RFflag;
-    string partdistr, stoplineartracking, crosssectionspath, opticsFile, plotflag, inputpath;
+    double Apr=0.0, Zpr=0.0, mass=0.0, energyPerIon=0.0, emix=0.0, emiy=0.0, sigdpp=0.0,
+    r1r2skin=0.0, blowup2=0.0, thicknessmagneticfield=0.0, Bmax=0.0, deltaGap=0.0, 
+    gamma=0.0, betgam=0.0, W, PlossPb, freqrf=0.0, rfvoltage=0.0, rfharmonic=0.0, 
+    wecolli=0.0, nsigi=0.0, momentum=0.0;
+    int npart=0, kbunch=0, taubeam=0, nparti=0, nrev1=0, nrev2=0, blowupperiod=0, scaleorbit=0,
+    beamflag=1, stopLinearTrackingNo=0, RunWithFluka=0, RunWithCrystal=0, RunningFlag=1,
+    outcoord=0, idpart=0, idelt=0, RFflag=0;
+    string partdistr="", stoplineartracking="", crosssectionspath="", opticsFile="", plotflag="", inputpath="";
 
     int indication(1);
@@ -1173 +1182 @@
 {
 
-
-    int JohnJaperFlag;
+//for the moment, ICOSIM++ cannot run without aperture parameters, so JohnJaperflag is fixed to 1; the variable is still here to facilitate further developpment
+    int JohnJaperFlag=1;
 
     defineParameters(opticsFile, beamflag, JohnJaperFlag);
 
     bool apertureIsZero(false);
-    int temp1(0);
-    int temp2(0);
-    int temp3;
-
+    int temp1 = 0;
+    int temp2 = 0;
+    int temp3 = 0;
+     
+    
     if (JohnJaperFlag == 1) { //JohnJaperflag=1 if we have information about the aperture
 
         for (int w(0); w < nbre_elts; ++w) {
             if (lat.reseau[w]->APERTYPE == "RECTELLIPSE") {
-                ++temp1;
+                temp1 = temp1 + 1;
                 if ((lat.reseau[w]->APER_1 == 0) && (lat.reseau[w]->APER_2 == 0) && (lat.reseau[w]->APER_3 == 0) && (lat.reseau[w]->APER_4 == 0)) {
-                    ++temp2;
-                } else {
-                    break;
-                }
-            } else {
-                break;
-            }
+                    temp2 = temp2 + 1;
+                } 
+//                 else {
+//                     break;
+//                 }
+            } 
+//             else {
+//                 break;
+//             }
         }
 
         //If all apertures are RECTELLIPSES:
-        //If an apeture variable is 0, it takes the value of the preceeding aperture variable.
-        //If all the apertures are 0, we define them to be 1 instead.
+           //If all the aperture variables of all lattice elements are 0, we define them to be 1 instead.
+        //  If an apeture variable is 0, it takes the value of the preceeding aperture variable.
 
         if (temp1 == nbre_elts) {
@@ -1214 +1226 @@
             if (lat.reseau[0]->APER_1 == 0) {
                 apertureIsZero = true;
-                temp3 = nbre_elts;
-                while (lat.reseau[temp3]->APER_1 == 0) {
+                temp3 = nbre_elts -1;
+                if (lat.reseau[temp3]->APER_1 == 0) {
                     temp3 = temp3 - 1;
                 }
@@ -1222 +1234 @@
             if (lat.reseau[0]->APER_2 == 0) {
                 apertureIsZero = true;
-                temp3 = nbre_elts;
-                while (lat.reseau[temp3]->APER_2 == 0) {
+                temp3 = nbre_elts -1;
+                if (lat.reseau[temp3]->APER_2 == 0) {
                     temp3 = temp3 - 1;
                 }
@@ -1230 +1242 @@
             if (lat.reseau[0]->APER_3 == 0) {
                 apertureIsZero = true;
-                temp3 = nbre_elts;
-                while (lat.reseau[temp3]->APER_3 == 0) {
+                temp3 = nbre_elts -1;
+                if (lat.reseau[temp3]->APER_3 == 0) {
                     temp3 = temp3 - 1;
                 }
@@ -1238 +1250 @@
             if (lat.reseau[0]->APER_4 == 0) {
                 apertureIsZero = true;
-                temp3 = nbre_elts;
-                while (lat.reseau[temp3]->APER_4 == 0) {
+                temp3 = nbre_elts - 1;
+                if (lat.reseau[temp3]->APER_4 == 0) {
                     temp3 = temp3 - 1;
                 }
@@ -1279 +1291 @@
                 }
             } else if (lat.reseau[r]->APERTYPE == "NONE") {
+              
+                if(r==0){
+                  lat.reseau[r]->APER_1 = lat.reseau[nbre_elts]->APER_1;
+                  lat.reseau[r]->APER_2 = lat.reseau[nbre_elts]->APER_2;
+                  lat.reseau[r]->APER_3 = lat.reseau[nbre_elts]->APER_3;
+                  lat.reseau[r]->APER_4 = lat.reseau[nbre_elts]->APER_4;
+                  lat.reseau[r]->APERTYPE = lat.reseau[nbre_elts]->APERTYPE;
+                }else{
                 lat.reseau[r]->APER_1 = lat.reseau[r - 1]->APER_1;
                 lat.reseau[r]->APER_2 = lat.reseau[r - 1]->APER_2;
@@ -1284 +1304 @@
                 lat.reseau[r]->APER_4 = lat.reseau[r - 1]->APER_4;
                 lat.reseau[r]->APERTYPE = lat.reseau[r - 1]->APERTYPE;
+                }
             }
         }
@@ -1291 +1312 @@
         }
 
-        //we compute the final apertures
+        //we compute the final apertures; different aperture type
         for (int i(0); i < lat.reseau.size(); ++i) {
             if (lat.reseau[i]->APERTYPE == "CIRCLE") {
@@ -1504 +1525 @@
     }
 
-    double hk;
+    double hk=0.0;
 
     for (int k(0); k < indice.size(); ++k) {
@@ -1529 +1550 @@
     }
 
-    double vk;
+    double vk=0.0;
 
     for (int k(0); k < indice.size(); ++k) {
@@ -1541 +1562 @@
 
 }
-
+/***************************************************************************
+ * 
+ * Read the lattice elements parameters
+ * 
+ * 
+ * 
+ * 
+ *  Jianfeng Zhang @ LAL, 21/03/2014
+ *   Fix the bug to read the lattice element parameters.
+ * 
+ ***************************************************************************/
 void Simulation::defineParameters(string opticsFile, const int& beamflag, int& JohnJaperFlag)
 {
 
     ifstream entree;
-    string APERTYPE, name, KEYWORD, PARENT;
+    string APERTYPE="", name="", KEYWORD="", PARENT="";
     entree.open(opticsFile.c_str());
 
@@ -1553 +1584 @@
     } else {
 
-        string name, temp;
-        string strtemp, strtemp1;
-
+        string name="", temp="";
+        string strtemp="", strtemp1="";
+         /* read the parameters at the beginning of the MADX twiss file*/
         getline(entree, strtemp, ',');
         while ((strtemp != "NAME") && (strtemp != "ALFX") && (strtemp != "S") && (strtemp != "L") && (strtemp != "K0L") && (strtemp != "K0SL") && (strtemp != "K1L") && (strtemp != "K1SL") && (strtemp != "K2L") && (strtemp != "K2SL") && (strtemp != "XC") && (strtemp != "PXC") && (strtemp != "KEYWORD") && (strtemp != "PARENT") && (strtemp != "YC") && (strtemp != "PYC") && (strtemp != "TC") && (strtemp != "PTC") && (strtemp != "BETX") && (strtemp != "MUX") && (strtemp != "DX") && (strtemp != "DPX") && (strtemp != "BETY") && (strtemp != "ALFY") && (strtemp != "MUY") && (strtemp != "DY") && (strtemp != "DPY") && (strtemp != "APERTYPE") && (strtemp != "APER_1") && (strtemp != "APER_2") && (strtemp != "APER_3") && (strtemp != "APER_4")) {
@@ -1561 +1592 @@
             getline(entree, strtemp, ',');
         }
-
-        double S, L, K0L, K0SL, K1L, K1SL, K2L, K2SL, XC, PXC, YC, PYC, TC, PTC, BETX, ALFX, MUX, DX, DPX, BETY, ALFY, MUY, DY, DPY, APER_1, APER_2, APER_3, APER_4;
+         /* read the parameters of the lattice elements*/
+        double S=0.0, L=0.0, K0L=0.0, K0SL=0.0, K1L=0.0, K1SL=0.0, K2L=0.0, K2SL=0.0, XC=0.0, PXC=0.0, YC=0.0, PYC=0.0, TC=0.0, PTC=0.0, 
+               BETX=0.0, ALFX=0.0, MUX=0.0, DX=0.0, DPX=0.0, BETY=0.0, ALFY=0.0, MUY=0.0, DY=0.0, DPY=0.0, APER_1=0.0, APER_2=0.0, APER_3=0.0, APER_4=0.0;
 
         int taille(0);
@@ -1612 +1644 @@
         }
         getline(entree, strtemp);
+        /* remove the "\r" symbol from the last variable in the "*.csv file", 
+               maybe this is a potential bug for other file format... By jianfeng Zhang @ LAL, 21/03/2014*/
+            strtemp = strtemp.substr(0,strtemp.size()-1);
         header.push_back(strtemp);
 
@@ -1623 +1658 @@
         }
 
+        /*start read the lattice elements parameters*/
         while (!entree.eof()) {
-            ++taille;
+           
+          
+               ++taille;
+            /* lattice elment*/
             for (int n(0); n < 31; ++n) {
                 getline(entree, strtemp, ',');
                 vectemp.push_back(strtemp);
             }
+            /* last variable of the lattice element */
             getline(entree, strtemp);
+            /* remove the "\r" symbol from the last variable in the "*.csv file", 
+               maybe this is a potential bug for other file format... By jianfeng Zhang @ LAL, 21/03/2014*/
+            strtemp = strtemp.substr(0,strtemp.size()-1);
             vectemp.push_back(strtemp);
 
+           
             for (int h(0); h < order.size(); ++h) {
 
@@ -1701 +1745 @@
             }
 
+             /*to fix the bug that: the c++ code can't stop when the end of fine is reached in the *.csv files
+               by Jianfeng zhang @ LAL, 21/03/2014
+             */
+            if(name.substr(0, 1) =="\r" || name.substr(0, 1) =="\n")
+            {
+              vectemp.clear();
+              taille = taille - 1; 
+              break;
+            }
+           
+             
+            
             Element elt(ALFX, ALFY, APER_1, APER_2, APER_3, APER_4, APERTYPE, BETX, BETY, DPX, DPY, DX, DY, KEYWORD, L, MUX, MUY, name, PTC, PXC, PYC, S, TC, XC, YC, K0L, K0SL, K1L, K1SL, K2L, K2SL, PARENT);
 
@@ -1707 +1763 @@
             vectemp.clear();
 
-            if ((name.substr(0, 3) == "END") || (name.substr(name.size() - 3, 3) == "END")) {
-                break;
-            }
-
-        }
+             if (name == "LHCB1$END") {
+                 cout << name << "         "<<taille<<endl;
+             }
+
+        }
+        
+        /*number of lattice elements*/
         this ->nbre_elts = taille;
     }