Changeset 8 in ZHANGProjects

Timestamp:

Apr 24, 2014, 4:01:58 PM (10 years ago)

Author:

zhangj

Message:

Fix the bug to read files.

Location:

ICOSIM/CPP/trunk/source

Files:

• .kdev4/source.kdev4 (modified) (1 diff)
• Particle.cc (modified) (1 diff)
• Particle.h (modified) (2 diffs)
• crystal.cc (modified) (1 diff)
• lattice.cc (modified) (6 diffs)
• simcrys.cc (modified) (6 diffs)
• simulation.cc (modified) (33 diffs)
• simulation.h (modified) (1 diff)

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

@@ -12 +12 @@
 
 [Launch][Launch Configuration 0][Data]
-Arguments=../sample/collimatorfile_standard_ions_nominal.csv   ../outputs
+Arguments=../sample/collimatorfile_crystal_matlab.csv    ../outputs/crystal
 Dependencies=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x00)
 Dependency Action=Nothing

ICOSIM/CPP/trunk/source/Particle.cc

@@ -165 +165 @@
 {
 
-    int nh;
-    double nr, sum, Rx, Ry;
+    int nh=0;
+    double nr=0.0, sum=0.0, Rx=0.0, Ry=0.0;
     vector < vector < double > > h1, h;
     vector < double > hsum1, hsum, xh, yh, xsh, ysh;

ICOSIM/CPP/trunk/source/Particle.h

@@ -26 +26 @@
     //==================Constructors, destructor=========================================
 
-    Particle(double x = 0, double dx = 0, double y = 0, double dy = 0, double deltaE = 0, double t = 0, double massnumber = 208, double chargestate = 82, double moment = 1.2);
+    Particle(double x = 0.0, double dx = 0.0, double y = 0.0, double dy = 0.0, double deltaE = 0.0, double t = 0.0, double massnumber = 208, double chargestate = 82, double moment = 1.2);
 
     Particle(const Particle& obj);
@@ -72 +72 @@
 
 
-    double Ap0; //massnumber of particle
-    double Zp0; //chargestate of particle
+    double Ap0; //atomic number of particle
+    double Zp0; //charge state of particle
     double dpoporiginal; //momentum (normalement on n en a pas besoin, = p1.coordonnees[][4]
     int inabs; //0 if gets lost, 1 if continue

ICOSIM/CPP/trunk/source/crystal.cc

@@ -17 +17 @@
     Alayer = 0;
 
-    cry_bend = (Cry_length / Rcurv); //total bending for this crystal;
-    s_length = Rcurv * (sin(cry_bend));
+    cry_bend = (Cry_length / Rcurv); //total bending angle for this crystal;
+    s_length = Rcurv * (sin(cry_bend)); // projected length [m]
     xmax = C_xmax ;
     L_chan = Cry_length; //Careful : the channeling length is the length of the particle trajectory, not the long. one

ICOSIM/CPP/trunk/source/lattice.cc

@@ -346 +346 @@
                         sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle
                         lcoll = resColli[ipcoll[i]]->L;//length of the collimator
+                        
 
                         if (sa == 0) {
@@ -393 +394 @@
         }//end loop over i
 
-        cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;
+       // cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;
 
         if (total - count == 0) { //we return if all the particles are gone
@@ -419 +420 @@
             if ((k + 1) % blowupperiod == 0) { //blowup/diffusion
                 if (bunch[p].in == 1) {
-                  //  cout << "Blow-up!!" << k+1 << endl;
-                    int im;
+                //    cout << "Blow-up!!" << k+1 << endl;
+                    int im=0;
                     im = reseau.size() - 1;
                     bunch[p].coordonnees[0][0] = (bunch[p].coordonnees[0][0] - reseau[im]->DX * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DX * bunch[p].coordonnees[0][4];
@@ -490 +491 @@
         cout << "Initialising R-matrix." << endl;
 
-        double K;
-        double cx, sx, cy, sy;
+        double K=0.0;
+        double cx=0.0, sx=0.0, cy=0.0, sy=0.0;
 
         R11X.clear();
@@ -524 +525 @@
 
             Lelem.push_back(reseau[i]->S - reseau[i - 1]->S);
+            
             if (reseau[i]->K1L > 0) {
                 if (Lelem[i] == 0) {
@@ -617 +619 @@
                 if (icop != -1) { //the element we consider is a collimator
 
+//cout << "icop = "<<icop << endl;
+
                     if (resColli[icop]->method == "standard") {
 

ICOSIM/CPP/trunk/source/simcrys.cc

@@ -436 +436 @@
 
     //-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------
-    int mat;
-
-    double p0;
-    double zlm;
-    double shift;
-    double x_shift, xp_shift, s_shift; //coordinates after shift/rotation
-    double x_rot, xp_rot, s_rot;
-    double x_temp, xp_temp, s_temp;  //all the _temp variables are used when you hit the cry from below
-    double tilt_int, x_int, xp_int, s_int;     //all the _int variables are used when you hit the cry from below (int=interaction point)
-    double x00;
-    double z;
-    double z00;
-    double zp;
-    double dpop;
-    double s;
-    double a_eq, b_eq, c_eq, Delta;
-    int part_abs;
-    double x;
-    double xp;
-    double y;
-    double yp;
-    double p;    //be careful: [Gev]
-    double s_in;
+    int mat=0;
+
+    double p0=0.0;
+    double zlm=0.0;
+    double shift=0.0;
+    double x_shift=0.0, xp_shift=0.0, s_shift=0.0; //coordinates after shift/rotation
+    double x_rot=0.0, xp_rot=0.0, s_rot=0.0;
+    double x_temp=0.0, xp_temp=0.0, s_temp=0.0;  //all the _temp variables are used when you hit the cry from below
+    double tilt_int=0.0, x_int=0.0, xp_int=0.0, s_int=0.0;//all the _int variables are used when you hit the cry from below (int=interaction point)
+    double x00=0.0;
+    double z=0.0;
+    double z00=0.0;
+    double zp=0.0;
+    double dpop=0.0;
+    double s=0.0;
+    double a_eq=0.0, b_eq=0.0, c_eq=0.0, Delta=0.0;
+    int part_abs=0.0;
+    double x=0.0;
+    double xp=0.0;
+    double y=0.0;
+    double yp=0.0;
+    double p=0.0;    //be careful: [Gev]
+    double s_in=0.0;
 
     //        adding variables for the pencil beam. Variables in the absolute reference frame.
 
-    double x_in0;
-    double xp_in0;
-    double y_in0;
-    double yp_in0;
-    double p_in0;    //be careful: [Gev]
-    double s_in0;
-    double s_impact;
+    double x_in0=0.0;
+    double xp_in0=0.0;
+    double y_in0=0.0;
+    double yp_in0=0.0;
+    double p_in0=0.0;    //be careful: [Gev]
+    double s_in0=0.0;
+    double s_impact=0.0;
     double totcount(0);
 
-    double mirror;
-    double tiltangle;
-
-
-    double c_length;      //Length in m
+    double mirror=0.0;
+    double tiltangle=0.0;
+
+
+    double c_length=0.0;      //Length in m
     double c_rotation(C_rotation) ;    //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED
     double c_aperture(C_aperture) ;   //Aperture in m
@@ -483 +483 @@
 
 
-    double xp_tangent;
+    double xp_tangent=0.0;
 
     double Cry_tilt(Crystal_tilt);     //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++////////////////
@@ -595 +595 @@
 
                 // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-
+               //change the beam frame to the crystal frame !!!!!!!!!
+               // since the interaction between the particle and the crystal is observed in the 
+               // crystal reference system
+
+               //initial coordinates in the beam frame
                 s   = 0;
                 x   = part.x;
@@ -619 +621 @@
                 shift = 0;
                 tilt_int = 0;
-                dpop = (p - p0) / p0;
+                dpop = (p - p0) / p0;  //momentum spread
 
                 // ------------------------------transform particle coordinates to get into collimator coordinate system -----------------------------------------
@@ -637 +639 @@
                     //cout<<"PARTICULE NON PERDU>>><<"<<endl;
 
+                    //?????????????? Need to define the correct value of the mirror for the install locations
+                    // of the crystal in SPS & LHC, inner side or outside of the vacuum chamber????
+                    //  by Zhang @ CERN, 23/04/2014.....
                     mirror = 1;
+                    
                     x  = mirror * x;
                     xp = mirror * xp;
@@ -953 +959 @@
 
 
+                        //For output.... by Zhang @ CERN, 24/04/2014
                         //p_in = (1 + dpop) * p0;
                         s_in = s_in + s;

ICOSIM/CPP/trunk/source/simulation.cc

@@ -92 +92 @@
     } else {
 
-        double temp2;
-
-        enterPart >> temp2;
-
-        while (!enterPart.eof()) {
+        double temp2=0.0;
+         //size_t  pos=0;
+        // string input;
+        
+     /*read one cubic number of the file*/
+        while (enterPart >> temp2) {
 
             Particle p;
@@ -113 +114 @@
 
             bunch.push_back(p);
-
-            enterPart >> temp2;
-        }
-    }
+        }
+    }
+    
+    enterPart.close();
 }
 
 //read information from collimator file
 
-void Simulation::readInputFile(string inputfile, double& Apr, double& Zpr, double& mass, double& energyPerIon, double& emix, double& emiy, int& npart, int& kbunch, double& sigdpp, int& taubeam, int& nparti, string& partdistr, double& r1r2skin, int& nrev1, int& nrev2, double& blowup2, int& blowupperiod, string& stoplineartracking, int& scaleorbit, double& wecolli, double& nsigi, string& opticsFile, double& thicknessmagneticfield, double& Bmax, double& deltaGap, int& beamflag, int& RunWithFluka, int& RunWithCrystal, double& freqrf, double& rfvoltage, double& rfharmonic, int& RunningFlag, int& idpart, int& idelt, int& outcoord, double& momentum, string& plotflag, string& inputpath, int& RFflag)
+void Simulation::readInputFile(string inputfile, double& Apr, double& Zpr, double& mass, double& energyPerIon, double& emix, double& emiy, long int& npart, int& kbunch, double& sigdpp, int& taubeam, long int& nparti, string& partdistr, double& r1r2skin, long int& nrev1, int& nrev2, double& blowup2, int& blowupperiod, string& stoplineartracking, int& scaleorbit, double& wecolli, double& nsigi, string& opticsFile, double& thicknessmagneticfield, double& Bmax, double& deltaGap, int& beamflag, int& RunWithFluka, int& RunWithCrystal, double& freqrf, double& rfvoltage, double& rfharmonic, int& RunningFlag, int& idpart, int& idelt, int& outcoord, double& momentum, string& plotflag, string& inputpath, int& RFflag)
 {
 
@@ -134 +135 @@
         //reading the parameter values in the top part of the collimator file
 
-        string name, comment, date, strtemp;
-        entree >> ws;
+        string name, comment, date;
+        string strtemp,ws;
+        size_t  pos=0;
+        string  delimiter=",";
+        
+        entree >> ws;
         getline(entree, comment);
         getline(entree, date);
 
-        while (!entree.eof()) {
-
-            entree >> ws;
-
-            getline(entree, name, ',');
-
-            if (name == "MASSNUMBER") {
-                getline(entree, strtemp);
+        while (entree >> ws) {
+ 
+             pos=ws.find(delimiter);
+            name=ws.substr(0,pos);
+            ws.erase(0,pos+delimiter.size());
+            
+            strtemp=ws;
+            
+           // Fix the bugs to read long particle numbers & set the RF flag, by Jianfeng Zhang @ LAL, 15/04/2014
+            if (name == "MASSNUMBER" && !strtemp.empty()) {
                 Apr = atof(strtemp.c_str());
-            } else if (name == "CHARGESTATE") {
-                getline(entree, strtemp);
+            } else if (name == "CHARGESTATE" && !strtemp.empty()) {
                 Zpr = atof(strtemp.c_str());
-            } else if (name == "MASS") {
-                getline(entree, strtemp);
+            } else if (name == "MASS" && !strtemp.empty()) {
                 mass = atof(strtemp.c_str());
-            } else if (name == "ENERGY") {
-                getline(entree, strtemp);
+            } else if (name == "ENERGY" && !strtemp.empty()) {
                 energyPerIon = atof(strtemp.c_str());
-            } else if (name == "EX") {
-                getline(entree, strtemp);
+            } else if (name == "EX" && !strtemp.empty()) {
                 emix = atof(strtemp.c_str());
-            } else if (name == "EY") {
-                getline(entree, strtemp);
+            } else if (name == "EY" && !strtemp.empty()) {
                 emiy = atof(strtemp.c_str());
-            } else if (name == "NPART") {
-                getline(entree, strtemp);
-                npart = atoi(strtemp.c_str());
-            } else if (name == "KBUNCH") {
-                getline(entree, strtemp);
+            } else if (name == "NPART" && !strtemp.empty()) {
+                npart = atol(strtemp.c_str());
+            } else if (name == "KBUNCH" && !strtemp.empty()) {
                 kbunch = atoi(strtemp.c_str());
-            } else if (name == "SIGDPP") {
-                getline(entree, strtemp);
+            } else if (name == "SIGDPP" && !strtemp.empty()) {
                 sigdpp = atof(strtemp.c_str());
-            } else if (name == "TAUBEAM") {
-                getline(entree, strtemp);
+            } else if (name == "TAUBEAM" && !strtemp.empty()) {
                 taubeam = atoi(strtemp.c_str());
-            } else if (name == "NPARTI") {
-                getline(entree, strtemp);
-                nparti = atoi(strtemp.c_str());
-            } else if (name == "PARTDISTR") {
-                getline(entree, partdistr);
-            } else if (name == "R1R2SKIN") {
-                if (partdistr == "r1r2") {
-                    getline(entree, strtemp);
+            } else if (name == "NPARTI" && !strtemp.empty()) {
+                nparti = atol(strtemp.c_str());
+            } else if (name == "PARTDISTR" && !strtemp.empty()) {
+                partdistr = strtemp;
+            } else if (name == "R1R2SKIN" && !strtemp.empty()) {
+                if (partdistr == "r1r2" && !strtemp.empty()) {
                     r1r2skin = atof(strtemp.c_str());
                 } else {
-                    getline(entree, strtemp);
-                }
-            } else if (name == "NREV1") {
-                getline(entree, strtemp);
-                nrev1 = atoi(strtemp.c_str());
-            } else if (name == "NREV2") {
-                getline(entree, strtemp);
+              //      getline(entree, strtemp);
+                }
+            } else if (name == "NREV1" && !strtemp.empty()){ 
+                nrev1 = atol(strtemp.c_str());
+            } else if (name == "NREV2" && !strtemp.empty()) {
                 nrev2 = atoi(strtemp.c_str());
-            } else if (name == "BLOWUP2") {
-                getline(entree, strtemp);
+            } else if (name == "BLOWUP2" && !strtemp.empty()) {
                 blowup2 = atof(strtemp.c_str());
-            } else if (name == "BLOWUPPERIOD") {
-                getline(entree, strtemp);
+            } else if (name == "BLOWUPPERIOD" && !strtemp.empty()) {
                 blowupperiod = atoi(strtemp.c_str());
-            } else if (name == "STOPLINEARTRACKING") {
-                getline(entree, stoplineartracking);
-            } else if (name == "SCALEORBIT") {
-                getline(entree, strtemp);
+            } else if (name == "STOPLINEARTRACKING" && !strtemp.empty()) {
+             stoplineartracking = strtemp;
+            } else if (name == "SCALEORBIT" && !strtemp.empty()) {
                 scaleorbit = atoi(strtemp.c_str());
-            } else if (name == "WECOLLI") {
-                getline(entree, strtemp);
+            } else if (name == "WECOLLI" && !strtemp.empty()) {
                 wecolli = atof(strtemp.c_str());
-            } else if (name == "NSIGI") {
-                getline(entree, strtemp);
+            } else if (name == "NSIGI" && !strtemp.empty()) {
                 nsigi = atof(strtemp.c_str());
-            } else if (name == "CROSSSECTIONSPATH") {
-                getline(entree, crosssectionspath);
-            } else if (name == "OPTICSPATH") {
-                getline(entree, opticsFile);
-            } else if (name == "INPUTPATH") {
-                getline(entree, inputpath);
-            } else if (name == "THICKNESSMAGNETICFIELD") {
-                getline(entree, strtemp);
+            } else if (name == "CROSSSECTIONSPATH" && !strtemp.empty()) {
+              crosssectionspath = strtemp;
+            } else if (name == "OPTICSPATH" && !strtemp.empty()) {
+              opticsFile = strtemp;
+            } else if (name == "INPUTPATH" && !strtemp.empty()) {
+             inputpath = strtemp;
+            } else if (name == "THICKNESSMAGNETICFIELD" && !strtemp.empty()) {
                 thicknessmagneticfield = atof(strtemp.c_str());
-            } else if (name == "BMAX") {
-                getline(entree, strtemp);
+            } else if (name == "BMAX" && !strtemp.empty()) {
                 Bmax = atof(strtemp.c_str());
-            } else if (name == "DELTAGAP") {
-                getline(entree, strtemp);
+            } else if (name == "DELTAGAP" && !strtemp.empty()) {
                 deltaGap = atof(strtemp.c_str());
-            } else if (name == "BEAMFLAG") {
-                getline(entree, strtemp);
+            } else if (name == "BEAMFLAG" && !strtemp.empty()) {
                 beamflag = atoi(strtemp.c_str());
-            } else if (name == "RUNNINGFLAG") {
-                getline(entree, strtemp);
+            } else if (name == "RUNNINGFLAG" && !strtemp.empty()) {
                 RunningFlag = atoi(strtemp.c_str());
-            } else if (name == "IDPART") {
-                getline(entree, strtemp);
+            } else if (name == "IDPART" && !strtemp.empty()) {
                 idpart = atoi(strtemp.c_str());
-            } else if (name == "OUTCOORD") {
-                getline(entree, strtemp);
+            } else if (name == "OUTCOORD" && !strtemp.empty()) {
                 outcoord = atoi(strtemp.c_str());
-            } else if (name == "IDELT") {
-                getline(entree, strtemp);
+            } else if (name == "IDELT" && !strtemp.empty()) {
                 idelt = atoi(strtemp.c_str());
-            } else if (name == "FREQRF") {
-                getline(entree, strtemp);
+            } else if (name == "FREQRF" && !strtemp.empty()) {
                 freqrf = atof(strtemp.c_str());
                 RFflag = 1;
-            } else if (name == "RFVOLTAGE") {
-                getline(entree, strtemp);
+            } else if (name == "RFVOLTAGE" && !strtemp.empty()) {
                 rfvoltage = atof(strtemp.c_str());
                 RFflag = 1;
-            } else if (name == "RFHARMONBER") {
-                getline(entree, strtemp);
+            } else if (name == "RFHARMONBER" && !strtemp.empty()) {
                 rfharmonic = atof(strtemp.c_str());
                 RFflag = 1;
-            } else if (name == "MOMENTUM") {
-                getline(entree, strtemp);
+            } else if (name == "MOMENTUM" && !strtemp.empty()) {
                 momentum = atof(strtemp.c_str());
-            } else if (name == "PLOTFLAG") {
-                getline(entree, plotflag);
+            } else if (name == "PLOTFLAG" && !strtemp.empty()) {
+                plotflag = strtemp;
             } else if (name == "NAME") {
                 break;
             } else {
-                getline(entree, strtemp);
+              //  getline(entree, strtemp);
             }
 
@@ -266 +243 @@
     }
     entree.close();
+    
+   
 
     //reads optics file
@@ -282 +261 @@
     } else {
 
-        string NAME, MATERIAL, METHOD, strtemp;
-        double ANGLE, LENGTH;
-        int PHASE;
-        double NSIG2;
+        string NAME, MATERIAL, METHOD, strtemp1,ws1;
+        string delimiter1 = ",";
+        size_t  pos=0;
+        
+        double ANGLE=0.0, LENGTH=0.0;
+        int PHASE=1;
+        double NSIG2=0.0;
 
 
         getline(entree1, NAME, ',');
         while (NAME != "NAME") {
-            getline(entree1, strtemp);
+            getline(entree1, strtemp1);
             getline(entree1, NAME, ',');
         }
-        getline(entree1, strtemp);
-
-        while (!entree1.eof()) {
-            getline(entree1, NAME, ',');
-            getline(entree1, strtemp, ',');
-            ANGLE = atof(strtemp.c_str());
-            getline(entree1, strtemp, ',');
-            LENGTH = atof(strtemp.c_str());
-            getline(entree1, strtemp, ',');
-            PHASE = atoi(strtemp.c_str());
-            getline(entree1, strtemp, ',');
-            NSIG2 = atof(strtemp.c_str());
-            getline(entree1, MATERIAL, ',');
-            getline(entree1, METHOD);
-
-            //we construct some special vectors related to collimators
-
+        getline(entree1, strtemp1);
+
+       // while (!entree1.eof()) {
+          while (entree1 >> ws1) {
+            
+            pos = ws1.find(delimiter1);
+            NAME = ws1.substr(0,pos);
+            ws1.erase(0,pos+delimiter1.size());
+            
+            
+            pos = ws1.find(delimiter1);
+            strtemp1 = ws1.substr(0,pos);
+            ANGLE = atof(strtemp1.c_str());
+            ws1.erase(0,pos+delimiter1.size());
+            
+           
+            pos = ws1.find(delimiter1);
+            strtemp1 = ws1.substr(0,pos);
+            LENGTH = atof(strtemp1.c_str());
+            ws1.erase(0,pos+delimiter1.size());
+            
+            pos = ws1.find(delimiter1);
+            strtemp1 = ws1.substr(0,pos);
+             PHASE = atoi(strtemp1.c_str());
+            ws1.erase(0,pos+delimiter1.size());
+            
+            pos = ws1.find(delimiter1);
+            strtemp1 = ws1.substr(0,pos);
+            NSIG2 = atof(strtemp1.c_str());
+            ws1.erase(0,pos+delimiter1.size());
+            
+            
+            pos = ws1.find(delimiter1);
+            MATERIAL = ws1.substr(0,pos);
+            ws1.erase(0,pos+delimiter1.size());
+            
+            METHOD = ws1;
+            ws1.erase(0,ws1.size());
+        
+         
+
+//          cout << NAME << endl;
+//          cout << ANGLE << endl;
+//          cout << LENGTH << endl;
+//          cout << PHASE << endl;
+//          cout << NSIG2 << endl;
+//          cout << MATERIAL << endl;
+//          cout << METHOD << endl;
+//          
 
             for (int i(0); i < nbre_elts; ++i) {
+              
+//            if(i==nbre_elts-1)
+//            cout << nbre_elts - 1 << endl;
+              
                 if (lat.reseau[i]->NAME == NAME) {
+                  
+                 // cout << "NAME is " << NAME << "i = " << i << endl; 
                     //we skip collimators that are not found in the optics file and collimators with names ending in the wrong letter
                     if ((atoi(&NAME[NAME.size() - 1]) == beamflag) || (NAME == "TCTV.4R2.B2") || (NAME == "TCTV.4R8.B2") || (NAME == "TCTV.4L2.B1") || (NAME == "TCTV.4L8.B1")) { //attention, il faut peut etre egalement supprimer ces elements dans le reseau
-
+                       
                         lat.ips.push_back(i);
 
@@ -328 +348 @@
                         } else if (NAME.substr(0, 4) == "TCLA") {
                             lat.itcla.push_back(i);
-                        }
+                        }else
+                          cout << "Simulation:readInputFile(): Warning! The name of the collimator ("<<
+                                   NAME <<") does not start with 'TCI', 'TCP', 'TCRYO', 'CRYSTAL', 'TCS', 'TCT'  or 'TCLA'"<<endl;;
 
 
@@ -353 +375 @@
                             CrystalCollimator crystalcolli(*lat.reseau[i], ANGLE, NSIG2, METHOD);
                             temp.push_back(new CrystalCollimator(crystalcolli));
-
                         } else {
 
@@ -364 +385 @@
     }
 
+    entree1.close();
 
     sort(lat.ip.begin(), lat.ip.end());
@@ -456 +478 @@
                         getline(readinfocrystal, temp1);
                         lat.resColli[j]->Cry_tilt = atof(temp1.c_str());
-                    } else {
-                        getline(readinfocrystal, temp1);
-                    }
+                    } //else {
+                      //cout << "Simulation::readInputFile(): Warning: the crystal ("
+                       //    << temp1 <<") is not defined in the collimator file" <<endl;
+                      //  getline(readinfocrystal, temp1);
+                    //}
                     getline(readinfocrystal, temp1, ',');
                 }
@@ -510 +534 @@
 
     for (int k(0); k < lat.resColli.size(); ++k) {
-        lat.resColli[k]->phi = (lat.resColli[k]->hgap2 - lat.resColli[k]->hgap) / lat.resColli[k]->L; //Calculate absolute angle of collimator
+        lat.resColli[k]->phi = (lat.resColli[k]->hgap2 - lat.resColli[k]->hgap) / lat.resColli[k]->L; //Calculate absolute angle of collimator; installation angle
     }
 
@@ -523 +547 @@
 
 
+    for(int k(0); k < lat.resColli.size(); ++k){
+     cout << "test.....";
+      cout << lat.resColli[k]->hgap <<"    " <<lat.resColli[k]->hgap2 <<endl;  
+    
+    }
+      
     lat.freqrf = freqrf;
     lat.rfvoltage = rfvoltage;
@@ -542 +572 @@
     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, 
+    gamma=0.0, betgam=0.0, W=0.0, PlossPb=0.0, 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,
+    long int npart=0L, nparti=0L, nrev1=0L;
+    int kbunch=0, taubeam=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;
@@ -672 +703 @@
 
     //================================================================= END OF TEMP. MODIFICATIONS ===========================================================================================
-    int particlehit;
+    int particlehit=0;
 
     if ((RunningFlag == 1) || (RunningFlag == 2)) {
         cout << "FIRST TRACKING." << endl;
 
-        //We track the particles linearly between the primary collimators
+        //We track the particles linearly between the primary collimators, to find the particles hit on the collimators in "nrev1" turns
+        // to save the simulation time in the following nonlinear tracking...
         lat.trackensemblelinearnew(bunch, bunchhit, nrev1, blowup2, blowupperiod);
         particlehit = bunchhit.size();
@@ -772 +804 @@
     vector <vector <double> > xco, yco;
     int irev(0);
-    int nonlinflag;
+    int nonlinflag(0);
     double attr1(0);
     nonlinflag = 0;//the second tracking is linear
@@ -1459 +1491 @@
         }
 
-        double mahx;
+        double mahx=0.0;
 
         mahx = *max_element(aperx2.begin(), aperx2.end());
@@ -1476 +1508 @@
         }
 
-        double mahy;
+        double mahy=0.0;
 
         mahy = *max_element(apery2.begin(), apery2.end());
@@ -1569 +1601 @@
  * 
  * 
- *  Jianfeng Zhang @ LAL, 21/03/2014
- *   Fix the bug to read the lattice element parameters.
+ *  Jianfeng Zhang @ LAL, 15/04/2014
+ *   Fix the bug to read the optics files of the lattice element parameters.
+ *   Now the code can correctly the external file until the end of the file is 
+ *   arrived.
  * 
  ***************************************************************************/
@@ -1603 +1637 @@
         vector <int> order;
 
+        size_t pos=0;
+        std::string input;
+        std::string delimiter = ",";
+        
         //header_base is: NAME,KEYWORD,PARENT,S,L,K0L,K0SL,K1L,K1SL,K2L,K2SL,XC,PXC,YC,PYC,TC,PTC,BETX,ALFX,MUX,DX,DPX,BETY,ALFY,MUY,DY,DPY,APERTYPE,APER_1,APER_2,APER_3,APER_4
         //we can have any order of the headers in the opticsfile
@@ -1659 +1697 @@
 
         /*start read the lattice elements parameters*/
-        while (!entree.eof()) {
-           
+        while (entree >> input) {
           
                ++taille;
             /* lattice elment*/
-            for (int n(0); n < 31; ++n) {
-                getline(entree, strtemp, ',');
-                vectemp.push_back(strtemp);
+            while ((pos=input.find(delimiter)) !=std::string::npos) {
+              strtemp = input.substr(0,pos);  
+               vectemp.push_back(strtemp);
+                input.erase(0,pos+delimiter.size());
             }
             /* 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);
-
+            strtemp=input;
+            vectemp.push_back(strtemp);
+           
+           
            
             for (int h(0); h < order.size(); ++h) {
@@ -1745 +1781 @@
             }
 
-             /*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);
@@ -1763 +1788 @@
             vectemp.clear();
 
-             if (name == "LHCB1$END") {
-                 cout << name << "         "<<taille<<endl;
-             }
-
-        }
-        
+            
+
+        }
+    //    cout << name << "         "<<taille<<endl;
         /*number of lattice elements*/
         this ->nbre_elts = taille;
@@ -1838 +1861 @@
 }
 
-
+/**************************************************************************************
+ * 
+ * 
+ * 
+ * 
+ * 
+ ***************************************************************************************/
 void Simulation::writeoutput(string outputfile, const double& emix, const double& emiy, const double& Apr, const double& Zpr, const double& nparti, const double& beamflag, const double& PlossPb, const int& i0, const int& im, const int& npart, const int& kbunch, const vector <int>& asumrem, const vector <int>& asumhitcolli, const vector <int>& asumhits, const int& taubeam, vector <double>& hitso, const vector <int>& nhitcollio, vector <double>& Aphito, vector <double>& Zphito, const int& ibeg, const int& iend, vector <vector <double> >& xco, vector <vector <double> >& yco, const string& plotflag, string outputpath, string inputpath)
 {
@@ -1851 +1880 @@
 
         for (int i(0); i < lat.ip.size(); ++i) {
-            sortie << lat.ip[i] << ", ";
+            sortie << lat.ip[i]+1 << ", ";
         }
 
@@ -1859 +1888 @@
 
         for (int i(0); i < lat.is.size(); ++i) {
-            sortie << lat.is[i] << ", ";
+            sortie << lat.is[i]+1 << ", ";
         }
 
@@ -1867 +1896 @@
 
         for (int i(0); i < lat.ips.size(); ++i) {
-            sortie << lat.ips[i] << ", ";
+            sortie << lat.ips[i]+1 << ", ";
         }
 
@@ -1884 +1913 @@
         sortie << "Charge state of the reference particle (Zpr): " << Zpr << endl;
         sortie << "Number of particles in the simulation (nparti): " << nparti << endl;
-        sortie << "Average power loss of the beam (PlossPb): " << PlossPb << endl;
+        sortie << "Average power loss of the beam (PlossPb [J/s]): " << PlossPb << endl;
 
         sortie << "Mass for the particles hitting the aperture (Aphit):";
@@ -1944 +1973 @@
 
 
-        sortie << "Expected lifetime of the beam (taubeam): " << taubeam;
+        sortie << "Expected lifetime of the beam (taubeam [second]): " << taubeam;
 
         sortie.close();
@@ -1960 +1989 @@
 
 
-
+/**************************************************************************************
+ * 
+ * 
+ * 
+ * Fix the bugs to save ips, ip, is, it, itcla. By Jianfeng Zhang @ LAL, 16/04/2014.
+ * 
+ ***************************************************************************************/
 void Simulation::PlotRunSummary(const vector <int>& asumrem, const vector <int>& asumhitcolli, const vector <int>& asumhits, const vector <int>& nhitcollio, vector <double>& hitso, const double& PlossPb, string outputpath)
 {
@@ -1967 +2002 @@
 
     int NLostTotal(0);
-    int irevloc, i0loc, imloc;
+    int irevloc(0), i0loc(0), imloc(0);
 
     NLostTotal = asumhitcolli[asumhitcolli.size() - 1] + asumhits[asumhits.size() - 1];
@@ -2100 +2135 @@
     }
 
+    //plot the particle hit the collimators with nsig <= 100.
     ofstream outplot3;
     string plot3;
@@ -2127 +2163 @@
 
         for (int i(0); i < newnhitcolli.size(); ++i) {
-            gnuout2 << "set label " << i + 1 << " " << '"' << lat.reseau[ipsnew[i]]->NAME << '"' << """ at " << i << ",13 rotate front" << endl;
+            gnuout2 << "set label " << i + 1 << " " << '"' << lat.reseau[ipsnew[i]]->NAME << '"' << """ at " << i+1 << ",13 rotate front" << endl;
         }
         gnuout2.close();
@@ -2149 +2185 @@
 
     //Preparing files for the plotting of the number of particles left in the beam, the number of particles that have hit a collimator and the number of particles that have hit the aperture
+//asumrem[i]: particle rest in the beam
+//asumhitcolli[i]: particle lost on the collimators
+//asumhits[i]: particle lost elsewhere
 
     ofstream outplot4;
@@ -2554 +2593 @@
 }
 
-
+/***********************************************************************************************************
+ *
+ * 
+ * 
+ * ********************************************************************************************************/
 void Simulation::PlotLossSpectra(const double& beamflag, const double& nparti, const vector <int>& asumrem, const double& Apr, const double& Zpr, vector <double>& hitso, vector <double>& Aphito, vector <double>& Zphito, const double& PlossPb, const int& npart, const int& kbunch, const int& taubeam, const vector <int>& nhitcollio, string outputpath, string inputpath)
 {

ICOSIM/CPP/trunk/source/simulation.h

@@ -53 +53 @@
     //method used to read the parameters in the collimator file
 
-    void readInputFile(string inputfile, double& Apr, double& Zpr, double& mass, double& energyPerIon, double& emix, double& emiy, int& npart, int& kbunch, double& sigdpp, int& taubeam, int& nparti, string& partdistr, double& r1r2skin, int& nrev1, int& nrev2, double& blowup2, int& blowupperiod, string& stoplineartracking, int& scaleorbit, double& wecolli, double& nsigi, string& opticsFile, double& thicknessmagneticfield, double& Bmax, double& deltaGap, int& beamflag, int& RunWithFluka, int& RunWithCrystal, double& freqrf, double& rfvoltage, double& rfharmonic, int& RunningFlag, int& idpart, int& idelt, int& outcoord, double& momentum, string& plotflag, string& inputpath, int& RFflag);
+    void readInputFile(string inputfile, double& Apr, double& Zpr, double& mass, double& energyPerIon, double& emix, double& emiy, long int& npart, int& kbunch, double& sigdpp, int& taubeam, long int& nparti, string& partdistr, double& r1r2skin, long int& nrev1, int& nrev2, double& blowup2, int& blowupperiod, string& stoplineartracking, int& scaleorbit, double& wecolli, double& nsigi, string& opticsFile, double& thicknessmagneticfield, double& Bmax, double& deltaGap, int& beamflag, int& RunWithFluka, int& RunWithCrystal, double& freqrf, double& rfvoltage, double& rfharmonic, int& RunningFlag, int& idpart, int& idelt, int& outcoord, double& momentum, string& plotflag, string& inputpath, int& RFflag);