- Timestamp:
- Apr 24, 2014, 4:01:58 PM (10 years ago)
- Location:
- ICOSIM/CPP/trunk/source
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
ICOSIM/CPP/trunk/source/.kdev4/source.kdev4
r7 r8 12 12 13 13 [Launch][Launch Configuration 0][Data] 14 Arguments=../sample/collimatorfile_ standard_ions_nominal.csv ../outputs14 Arguments=../sample/collimatorfile_crystal_matlab.csv ../outputs/crystal 15 15 Dependencies=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x00) 16 16 Dependency Action=Nothing -
ICOSIM/CPP/trunk/source/Particle.cc
r7 r8 165 165 { 166 166 167 int nh ;168 double nr , sum, Rx, Ry;167 int nh=0; 168 double nr=0.0, sum=0.0, Rx=0.0, Ry=0.0; 169 169 vector < vector < double > > h1, h; 170 170 vector < double > hsum1, hsum, xh, yh, xsh, ysh; -
ICOSIM/CPP/trunk/source/Particle.h
r5 r8 26 26 //==================Constructors, destructor========================================= 27 27 28 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);28 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); 29 29 30 30 Particle(const Particle& obj); … … 72 72 73 73 74 double Ap0; // massnumber of particle75 double Zp0; //charge state of particle74 double Ap0; //atomic number of particle 75 double Zp0; //charge state of particle 76 76 double dpoporiginal; //momentum (normalement on n en a pas besoin, = p1.coordonnees[][4] 77 77 int inabs; //0 if gets lost, 1 if continue -
ICOSIM/CPP/trunk/source/crystal.cc
r5 r8 17 17 Alayer = 0; 18 18 19 cry_bend = (Cry_length / Rcurv); //total bending for this crystal;20 s_length = Rcurv * (sin(cry_bend)); 19 cry_bend = (Cry_length / Rcurv); //total bending angle for this crystal; 20 s_length = Rcurv * (sin(cry_bend)); // projected length [m] 21 21 xmax = C_xmax ; 22 22 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
r7 r8 346 346 sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle 347 347 lcoll = resColli[ipcoll[i]]->L;//length of the collimator 348 348 349 349 350 if (sa == 0) { … … 393 394 }//end loop over i 394 395 395 cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;396 // cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl; 396 397 397 398 if (total - count == 0) { //we return if all the particles are gone … … 419 420 if ((k + 1) % blowupperiod == 0) { //blowup/diffusion 420 421 if (bunch[p].in == 1) { 421 //cout << "Blow-up!!" << k+1 << endl;422 int im ;422 // cout << "Blow-up!!" << k+1 << endl; 423 int im=0; 423 424 im = reseau.size() - 1; 424 425 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; 491 492 492 double K ;493 double cx , sx, cy, sy;493 double K=0.0; 494 double cx=0.0, sx=0.0, cy=0.0, sy=0.0; 494 495 495 496 R11X.clear(); … … 524 525 525 526 Lelem.push_back(reseau[i]->S - reseau[i - 1]->S); 527 526 528 if (reseau[i]->K1L > 0) { 527 529 if (Lelem[i] == 0) { … … 617 619 if (icop != -1) { //the element we consider is a collimator 618 620 621 //cout << "icop = "<<icop << endl; 622 619 623 if (resColli[icop]->method == "standard") { 620 624 -
ICOSIM/CPP/trunk/source/simcrys.cc
r5 r8 436 436 437 437 //-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------ 438 int mat ;439 440 double p0 ;441 double zlm ;442 double shift ;443 double x_shift , xp_shift, s_shift; //coordinates after shift/rotation444 double x_rot , xp_rot, s_rot;445 double x_temp , xp_temp, s_temp; //all the _temp variables are used when you hit the cry from below446 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)447 double x00 ;448 double z ;449 double z00 ;450 double zp ;451 double dpop ;452 double s ;453 double a_eq , b_eq, c_eq, Delta;454 int part_abs ;455 double x ;456 double xp ;457 double y ;458 double yp ;459 double p ; //be careful: [Gev]460 double s_in ;438 int mat=0; 439 440 double p0=0.0; 441 double zlm=0.0; 442 double shift=0.0; 443 double x_shift=0.0, xp_shift=0.0, s_shift=0.0; //coordinates after shift/rotation 444 double x_rot=0.0, xp_rot=0.0, s_rot=0.0; 445 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 446 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) 447 double x00=0.0; 448 double z=0.0; 449 double z00=0.0; 450 double zp=0.0; 451 double dpop=0.0; 452 double s=0.0; 453 double a_eq=0.0, b_eq=0.0, c_eq=0.0, Delta=0.0; 454 int part_abs=0.0; 455 double x=0.0; 456 double xp=0.0; 457 double y=0.0; 458 double yp=0.0; 459 double p=0.0; //be careful: [Gev] 460 double s_in=0.0; 461 461 462 462 // adding variables for the pencil beam. Variables in the absolute reference frame. 463 463 464 double x_in0 ;465 double xp_in0 ;466 double y_in0 ;467 double yp_in0 ;468 double p_in0 ; //be careful: [Gev]469 double s_in0 ;470 double s_impact ;464 double x_in0=0.0; 465 double xp_in0=0.0; 466 double y_in0=0.0; 467 double yp_in0=0.0; 468 double p_in0=0.0; //be careful: [Gev] 469 double s_in0=0.0; 470 double s_impact=0.0; 471 471 double totcount(0); 472 472 473 double mirror ;474 double tiltangle ;475 476 477 double c_length ; //Length in m473 double mirror=0.0; 474 double tiltangle=0.0; 475 476 477 double c_length=0.0; //Length in m 478 478 double c_rotation(C_rotation) ; //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED 479 479 double c_aperture(C_aperture) ; //Aperture in m … … 483 483 484 484 485 double xp_tangent ;485 double xp_tangent=0.0; 486 486 487 487 double Cry_tilt(Crystal_tilt); //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++//////////////// … … 595 595 596 596 // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 597 598 599 597 //change the beam frame to the crystal frame !!!!!!!!! 598 // since the interaction between the particle and the crystal is observed in the 599 // crystal reference system 600 601 //initial coordinates in the beam frame 600 602 s = 0; 601 603 x = part.x; … … 619 621 shift = 0; 620 622 tilt_int = 0; 621 dpop = (p - p0) / p0; 623 dpop = (p - p0) / p0; //momentum spread 622 624 623 625 // ------------------------------transform particle coordinates to get into collimator coordinate system ----------------------------------------- … … 637 639 //cout<<"PARTICULE NON PERDU>>><<"<<endl; 638 640 641 //?????????????? Need to define the correct value of the mirror for the install locations 642 // of the crystal in SPS & LHC, inner side or outside of the vacuum chamber???? 643 // by Zhang @ CERN, 23/04/2014..... 639 644 mirror = 1; 645 640 646 x = mirror * x; 641 647 xp = mirror * xp; … … 953 959 954 960 961 //For output.... by Zhang @ CERN, 24/04/2014 955 962 //p_in = (1 + dpop) * p0; 956 963 s_in = s_in + s; -
ICOSIM/CPP/trunk/source/simulation.cc
r7 r8 92 92 } else { 93 93 94 double temp2; 95 96 enterPart >> temp2; 97 98 while (!enterPart.eof()) { 94 double temp2=0.0; 95 //size_t pos=0; 96 // string input; 97 98 /*read one cubic number of the file*/ 99 while (enterPart >> temp2) { 99 100 100 101 Particle p; … … 113 114 114 115 bunch.push_back(p); 115 116 enterPart >> temp2;117 }118 }116 } 117 } 118 119 enterPart.close(); 119 120 } 120 121 121 122 //read information from collimator file 122 123 123 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)124 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) 124 125 { 125 126 … … 134 135 //reading the parameter values in the top part of the collimator file 135 136 136 string name, comment, date, strtemp; 137 entree >> ws; 137 string name, comment, date; 138 string strtemp,ws; 139 size_t pos=0; 140 string delimiter=","; 141 142 entree >> ws; 138 143 getline(entree, comment); 139 144 getline(entree, date); 140 145 141 while (!entree.eof()) { 142 143 entree >> ws; 144 145 getline(entree, name, ','); 146 147 if (name == "MASSNUMBER") { 148 getline(entree, strtemp); 146 while (entree >> ws) { 147 148 pos=ws.find(delimiter); 149 name=ws.substr(0,pos); 150 ws.erase(0,pos+delimiter.size()); 151 152 strtemp=ws; 153 154 // Fix the bugs to read long particle numbers & set the RF flag, by Jianfeng Zhang @ LAL, 15/04/2014 155 if (name == "MASSNUMBER" && !strtemp.empty()) { 149 156 Apr = atof(strtemp.c_str()); 150 } else if (name == "CHARGESTATE") { 151 getline(entree, strtemp); 157 } else if (name == "CHARGESTATE" && !strtemp.empty()) { 152 158 Zpr = atof(strtemp.c_str()); 153 } else if (name == "MASS") { 154 getline(entree, strtemp); 159 } else if (name == "MASS" && !strtemp.empty()) { 155 160 mass = atof(strtemp.c_str()); 156 } else if (name == "ENERGY") { 157 getline(entree, strtemp); 161 } else if (name == "ENERGY" && !strtemp.empty()) { 158 162 energyPerIon = atof(strtemp.c_str()); 159 } else if (name == "EX") { 160 getline(entree, strtemp); 163 } else if (name == "EX" && !strtemp.empty()) { 161 164 emix = atof(strtemp.c_str()); 162 } else if (name == "EY") { 163 getline(entree, strtemp); 165 } else if (name == "EY" && !strtemp.empty()) { 164 166 emiy = atof(strtemp.c_str()); 165 } else if (name == "NPART") { 166 getline(entree, strtemp); 167 npart = atoi(strtemp.c_str()); 168 } else if (name == "KBUNCH") { 169 getline(entree, strtemp); 167 } else if (name == "NPART" && !strtemp.empty()) { 168 npart = atol(strtemp.c_str()); 169 } else if (name == "KBUNCH" && !strtemp.empty()) { 170 170 kbunch = atoi(strtemp.c_str()); 171 } else if (name == "SIGDPP") { 172 getline(entree, strtemp); 171 } else if (name == "SIGDPP" && !strtemp.empty()) { 173 172 sigdpp = atof(strtemp.c_str()); 174 } else if (name == "TAUBEAM") { 175 getline(entree, strtemp); 173 } else if (name == "TAUBEAM" && !strtemp.empty()) { 176 174 taubeam = atoi(strtemp.c_str()); 177 } else if (name == "NPARTI") { 178 getline(entree, strtemp); 179 nparti = atoi(strtemp.c_str()); 180 } else if (name == "PARTDISTR") { 181 getline(entree, partdistr); 182 } else if (name == "R1R2SKIN") { 183 if (partdistr == "r1r2") { 184 getline(entree, strtemp); 175 } else if (name == "NPARTI" && !strtemp.empty()) { 176 nparti = atol(strtemp.c_str()); 177 } else if (name == "PARTDISTR" && !strtemp.empty()) { 178 partdistr = strtemp; 179 } else if (name == "R1R2SKIN" && !strtemp.empty()) { 180 if (partdistr == "r1r2" && !strtemp.empty()) { 185 181 r1r2skin = atof(strtemp.c_str()); 186 182 } else { 187 getline(entree, strtemp); 188 } 189 } else if (name == "NREV1") { 190 getline(entree, strtemp); 191 nrev1 = atoi(strtemp.c_str()); 192 } else if (name == "NREV2") { 193 getline(entree, strtemp); 183 // getline(entree, strtemp); 184 } 185 } else if (name == "NREV1" && !strtemp.empty()){ 186 nrev1 = atol(strtemp.c_str()); 187 } else if (name == "NREV2" && !strtemp.empty()) { 194 188 nrev2 = atoi(strtemp.c_str()); 195 } else if (name == "BLOWUP2") { 196 getline(entree, strtemp); 189 } else if (name == "BLOWUP2" && !strtemp.empty()) { 197 190 blowup2 = atof(strtemp.c_str()); 198 } else if (name == "BLOWUPPERIOD") { 199 getline(entree, strtemp); 191 } else if (name == "BLOWUPPERIOD" && !strtemp.empty()) { 200 192 blowupperiod = atoi(strtemp.c_str()); 201 } else if (name == "STOPLINEARTRACKING") { 202 getline(entree, stoplineartracking); 203 } else if (name == "SCALEORBIT") { 204 getline(entree, strtemp); 193 } else if (name == "STOPLINEARTRACKING" && !strtemp.empty()) { 194 stoplineartracking = strtemp; 195 } else if (name == "SCALEORBIT" && !strtemp.empty()) { 205 196 scaleorbit = atoi(strtemp.c_str()); 206 } else if (name == "WECOLLI") { 207 getline(entree, strtemp); 197 } else if (name == "WECOLLI" && !strtemp.empty()) { 208 198 wecolli = atof(strtemp.c_str()); 209 } else if (name == "NSIGI") { 210 getline(entree, strtemp); 199 } else if (name == "NSIGI" && !strtemp.empty()) { 211 200 nsigi = atof(strtemp.c_str()); 212 } else if (name == "CROSSSECTIONSPATH") { 213 getline(entree, crosssectionspath); 214 } else if (name == "OPTICSPATH") { 215 getline(entree, opticsFile); 216 } else if (name == "INPUTPATH") { 217 getline(entree, inputpath); 218 } else if (name == "THICKNESSMAGNETICFIELD") { 219 getline(entree, strtemp); 201 } else if (name == "CROSSSECTIONSPATH" && !strtemp.empty()) { 202 crosssectionspath = strtemp; 203 } else if (name == "OPTICSPATH" && !strtemp.empty()) { 204 opticsFile = strtemp; 205 } else if (name == "INPUTPATH" && !strtemp.empty()) { 206 inputpath = strtemp; 207 } else if (name == "THICKNESSMAGNETICFIELD" && !strtemp.empty()) { 220 208 thicknessmagneticfield = atof(strtemp.c_str()); 221 } else if (name == "BMAX") { 222 getline(entree, strtemp); 209 } else if (name == "BMAX" && !strtemp.empty()) { 223 210 Bmax = atof(strtemp.c_str()); 224 } else if (name == "DELTAGAP") { 225 getline(entree, strtemp); 211 } else if (name == "DELTAGAP" && !strtemp.empty()) { 226 212 deltaGap = atof(strtemp.c_str()); 227 } else if (name == "BEAMFLAG") { 228 getline(entree, strtemp); 213 } else if (name == "BEAMFLAG" && !strtemp.empty()) { 229 214 beamflag = atoi(strtemp.c_str()); 230 } else if (name == "RUNNINGFLAG") { 231 getline(entree, strtemp); 215 } else if (name == "RUNNINGFLAG" && !strtemp.empty()) { 232 216 RunningFlag = atoi(strtemp.c_str()); 233 } else if (name == "IDPART") { 234 getline(entree, strtemp); 217 } else if (name == "IDPART" && !strtemp.empty()) { 235 218 idpart = atoi(strtemp.c_str()); 236 } else if (name == "OUTCOORD") { 237 getline(entree, strtemp); 219 } else if (name == "OUTCOORD" && !strtemp.empty()) { 238 220 outcoord = atoi(strtemp.c_str()); 239 } else if (name == "IDELT") { 240 getline(entree, strtemp); 221 } else if (name == "IDELT" && !strtemp.empty()) { 241 222 idelt = atoi(strtemp.c_str()); 242 } else if (name == "FREQRF") { 243 getline(entree, strtemp); 223 } else if (name == "FREQRF" && !strtemp.empty()) { 244 224 freqrf = atof(strtemp.c_str()); 245 225 RFflag = 1; 246 } else if (name == "RFVOLTAGE") { 247 getline(entree, strtemp); 226 } else if (name == "RFVOLTAGE" && !strtemp.empty()) { 248 227 rfvoltage = atof(strtemp.c_str()); 249 228 RFflag = 1; 250 } else if (name == "RFHARMONBER") { 251 getline(entree, strtemp); 229 } else if (name == "RFHARMONBER" && !strtemp.empty()) { 252 230 rfharmonic = atof(strtemp.c_str()); 253 231 RFflag = 1; 254 } else if (name == "MOMENTUM") { 255 getline(entree, strtemp); 232 } else if (name == "MOMENTUM" && !strtemp.empty()) { 256 233 momentum = atof(strtemp.c_str()); 257 } else if (name == "PLOTFLAG" ) {258 getline(entree, plotflag);234 } else if (name == "PLOTFLAG" && !strtemp.empty()) { 235 plotflag = strtemp; 259 236 } else if (name == "NAME") { 260 237 break; 261 238 } else { 262 getline(entree, strtemp);239 // getline(entree, strtemp); 263 240 } 264 241 … … 266 243 } 267 244 entree.close(); 245 246 268 247 269 248 //reads optics file … … 282 261 } else { 283 262 284 string NAME, MATERIAL, METHOD, strtemp; 285 double ANGLE, LENGTH; 286 int PHASE; 287 double NSIG2; 263 string NAME, MATERIAL, METHOD, strtemp1,ws1; 264 string delimiter1 = ","; 265 size_t pos=0; 266 267 double ANGLE=0.0, LENGTH=0.0; 268 int PHASE=1; 269 double NSIG2=0.0; 288 270 289 271 290 272 getline(entree1, NAME, ','); 291 273 while (NAME != "NAME") { 292 getline(entree1, strtemp );274 getline(entree1, strtemp1); 293 275 getline(entree1, NAME, ','); 294 276 } 295 getline(entree1, strtemp); 296 297 while (!entree1.eof()) { 298 getline(entree1, NAME, ','); 299 getline(entree1, strtemp, ','); 300 ANGLE = atof(strtemp.c_str()); 301 getline(entree1, strtemp, ','); 302 LENGTH = atof(strtemp.c_str()); 303 getline(entree1, strtemp, ','); 304 PHASE = atoi(strtemp.c_str()); 305 getline(entree1, strtemp, ','); 306 NSIG2 = atof(strtemp.c_str()); 307 getline(entree1, MATERIAL, ','); 308 getline(entree1, METHOD); 309 310 //we construct some special vectors related to collimators 311 277 getline(entree1, strtemp1); 278 279 // while (!entree1.eof()) { 280 while (entree1 >> ws1) { 281 282 pos = ws1.find(delimiter1); 283 NAME = ws1.substr(0,pos); 284 ws1.erase(0,pos+delimiter1.size()); 285 286 287 pos = ws1.find(delimiter1); 288 strtemp1 = ws1.substr(0,pos); 289 ANGLE = atof(strtemp1.c_str()); 290 ws1.erase(0,pos+delimiter1.size()); 291 292 293 pos = ws1.find(delimiter1); 294 strtemp1 = ws1.substr(0,pos); 295 LENGTH = atof(strtemp1.c_str()); 296 ws1.erase(0,pos+delimiter1.size()); 297 298 pos = ws1.find(delimiter1); 299 strtemp1 = ws1.substr(0,pos); 300 PHASE = atoi(strtemp1.c_str()); 301 ws1.erase(0,pos+delimiter1.size()); 302 303 pos = ws1.find(delimiter1); 304 strtemp1 = ws1.substr(0,pos); 305 NSIG2 = atof(strtemp1.c_str()); 306 ws1.erase(0,pos+delimiter1.size()); 307 308 309 pos = ws1.find(delimiter1); 310 MATERIAL = ws1.substr(0,pos); 311 ws1.erase(0,pos+delimiter1.size()); 312 313 METHOD = ws1; 314 ws1.erase(0,ws1.size()); 315 316 317 318 // cout << NAME << endl; 319 // cout << ANGLE << endl; 320 // cout << LENGTH << endl; 321 // cout << PHASE << endl; 322 // cout << NSIG2 << endl; 323 // cout << MATERIAL << endl; 324 // cout << METHOD << endl; 325 // 312 326 313 327 for (int i(0); i < nbre_elts; ++i) { 328 329 // if(i==nbre_elts-1) 330 // cout << nbre_elts - 1 << endl; 331 314 332 if (lat.reseau[i]->NAME == NAME) { 333 334 // cout << "NAME is " << NAME << "i = " << i << endl; 315 335 //we skip collimators that are not found in the optics file and collimators with names ending in the wrong letter 316 336 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 317 337 318 338 lat.ips.push_back(i); 319 339 … … 328 348 } else if (NAME.substr(0, 4) == "TCLA") { 329 349 lat.itcla.push_back(i); 330 } 350 }else 351 cout << "Simulation:readInputFile(): Warning! The name of the collimator ("<< 352 NAME <<") does not start with 'TCI', 'TCP', 'TCRYO', 'CRYSTAL', 'TCS', 'TCT' or 'TCLA'"<<endl;; 331 353 332 354 … … 353 375 CrystalCollimator crystalcolli(*lat.reseau[i], ANGLE, NSIG2, METHOD); 354 376 temp.push_back(new CrystalCollimator(crystalcolli)); 355 356 377 } else { 357 378 … … 364 385 } 365 386 387 entree1.close(); 366 388 367 389 sort(lat.ip.begin(), lat.ip.end()); … … 456 478 getline(readinfocrystal, temp1); 457 479 lat.resColli[j]->Cry_tilt = atof(temp1.c_str()); 458 } else { 459 getline(readinfocrystal, temp1); 460 } 480 } //else { 481 //cout << "Simulation::readInputFile(): Warning: the crystal (" 482 // << temp1 <<") is not defined in the collimator file" <<endl; 483 // getline(readinfocrystal, temp1); 484 //} 461 485 getline(readinfocrystal, temp1, ','); 462 486 } … … 510 534 511 535 for (int k(0); k < lat.resColli.size(); ++k) { 512 lat.resColli[k]->phi = (lat.resColli[k]->hgap2 - lat.resColli[k]->hgap) / lat.resColli[k]->L; //Calculate absolute angle of collimator 536 lat.resColli[k]->phi = (lat.resColli[k]->hgap2 - lat.resColli[k]->hgap) / lat.resColli[k]->L; //Calculate absolute angle of collimator; installation angle 513 537 } 514 538 … … 523 547 524 548 549 for(int k(0); k < lat.resColli.size(); ++k){ 550 cout << "test....."; 551 cout << lat.resColli[k]->hgap <<" " <<lat.resColli[k]->hgap2 <<endl; 552 553 } 554 525 555 lat.freqrf = freqrf; 526 556 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, 543 573 r1r2skin=0.0, blowup2=0.0, thicknessmagneticfield=0.0, Bmax=0.0, deltaGap=0.0, 544 gamma=0.0, betgam=0.0, W , PlossPb, freqrf=0.0, rfvoltage=0.0, rfharmonic=0.0,574 gamma=0.0, betgam=0.0, W=0.0, PlossPb=0.0, freqrf=0.0, rfvoltage=0.0, rfharmonic=0.0, 545 575 wecolli=0.0, nsigi=0.0, momentum=0.0; 546 int npart=0, kbunch=0, taubeam=0, nparti=0, nrev1=0, nrev2=0, blowupperiod=0, scaleorbit=0, 576 long int npart=0L, nparti=0L, nrev1=0L; 577 int kbunch=0, taubeam=0, nrev2=0, blowupperiod=0, scaleorbit=0, 547 578 beamflag=1, stopLinearTrackingNo=0, RunWithFluka=0, RunWithCrystal=0, RunningFlag=1, 548 579 outcoord=0, idpart=0, idelt=0, RFflag=0; … … 672 703 673 704 //================================================================= END OF TEMP. MODIFICATIONS =========================================================================================== 674 int particlehit ;705 int particlehit=0; 675 706 676 707 if ((RunningFlag == 1) || (RunningFlag == 2)) { 677 708 cout << "FIRST TRACKING." << endl; 678 709 679 //We track the particles linearly between the primary collimators 710 //We track the particles linearly between the primary collimators, to find the particles hit on the collimators in "nrev1" turns 711 // to save the simulation time in the following nonlinear tracking... 680 712 lat.trackensemblelinearnew(bunch, bunchhit, nrev1, blowup2, blowupperiod); 681 713 particlehit = bunchhit.size(); … … 772 804 vector <vector <double> > xco, yco; 773 805 int irev(0); 774 int nonlinflag ;806 int nonlinflag(0); 775 807 double attr1(0); 776 808 nonlinflag = 0;//the second tracking is linear … … 1459 1491 } 1460 1492 1461 double mahx ;1493 double mahx=0.0; 1462 1494 1463 1495 mahx = *max_element(aperx2.begin(), aperx2.end()); … … 1476 1508 } 1477 1509 1478 double mahy ;1510 double mahy=0.0; 1479 1511 1480 1512 mahy = *max_element(apery2.begin(), apery2.end()); … … 1569 1601 * 1570 1602 * 1571 * Jianfeng Zhang @ LAL, 21/03/2014 1572 * Fix the bug to read the lattice element parameters. 1603 * Jianfeng Zhang @ LAL, 15/04/2014 1604 * Fix the bug to read the optics files of the lattice element parameters. 1605 * Now the code can correctly the external file until the end of the file is 1606 * arrived. 1573 1607 * 1574 1608 ***************************************************************************/ … … 1603 1637 vector <int> order; 1604 1638 1639 size_t pos=0; 1640 std::string input; 1641 std::string delimiter = ","; 1642 1605 1643 //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 1606 1644 //we can have any order of the headers in the opticsfile … … 1659 1697 1660 1698 /*start read the lattice elements parameters*/ 1661 while (!entree.eof()) { 1662 1699 while (entree >> input) { 1663 1700 1664 1701 ++taille; 1665 1702 /* lattice elment*/ 1666 for (int n(0); n < 31; ++n) { 1667 getline(entree, strtemp, ','); 1668 vectemp.push_back(strtemp); 1703 while ((pos=input.find(delimiter)) !=std::string::npos) { 1704 strtemp = input.substr(0,pos); 1705 vectemp.push_back(strtemp); 1706 input.erase(0,pos+delimiter.size()); 1669 1707 } 1670 1708 /* last variable of the lattice element */ 1671 getline(entree, strtemp); 1672 /* remove the "\r" symbol from the last variable in the "*.csv file", 1673 maybe this is a potential bug for other file format... By jianfeng Zhang @ LAL, 21/03/2014*/ 1674 strtemp = strtemp.substr(0,strtemp.size()-1); 1675 vectemp.push_back(strtemp); 1676 1709 strtemp=input; 1710 vectemp.push_back(strtemp); 1711 1712 1677 1713 1678 1714 for (int h(0); h < order.size(); ++h) { … … 1745 1781 } 1746 1782 1747 /*to fix the bug that: the c++ code can't stop when the end of fine is reached in the *.csv files1748 by Jianfeng zhang @ LAL, 21/03/20141749 */1750 if(name.substr(0, 1) =="\r" || name.substr(0, 1) =="\n")1751 {1752 vectemp.clear();1753 taille = taille - 1;1754 break;1755 }1756 1757 1758 1783 1759 1784 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(); 1764 1789 1765 if (name == "LHCB1$END") { 1766 cout << name << " "<<taille<<endl; 1767 } 1768 1769 } 1770 1790 1791 1792 } 1793 // cout << name << " "<<taille<<endl; 1771 1794 /*number of lattice elements*/ 1772 1795 this ->nbre_elts = taille; … … 1838 1861 } 1839 1862 1840 1863 /************************************************************************************** 1864 * 1865 * 1866 * 1867 * 1868 * 1869 ***************************************************************************************/ 1841 1870 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) 1842 1871 { … … 1851 1880 1852 1881 for (int i(0); i < lat.ip.size(); ++i) { 1853 sortie << lat.ip[i] << ", ";1882 sortie << lat.ip[i]+1 << ", "; 1854 1883 } 1855 1884 … … 1859 1888 1860 1889 for (int i(0); i < lat.is.size(); ++i) { 1861 sortie << lat.is[i] << ", ";1890 sortie << lat.is[i]+1 << ", "; 1862 1891 } 1863 1892 … … 1867 1896 1868 1897 for (int i(0); i < lat.ips.size(); ++i) { 1869 sortie << lat.ips[i] << ", ";1898 sortie << lat.ips[i]+1 << ", "; 1870 1899 } 1871 1900 … … 1884 1913 sortie << "Charge state of the reference particle (Zpr): " << Zpr << endl; 1885 1914 sortie << "Number of particles in the simulation (nparti): " << nparti << endl; 1886 sortie << "Average power loss of the beam (PlossPb ): " << PlossPb << endl;1915 sortie << "Average power loss of the beam (PlossPb [J/s]): " << PlossPb << endl; 1887 1916 1888 1917 sortie << "Mass for the particles hitting the aperture (Aphit):"; … … 1944 1973 1945 1974 1946 sortie << "Expected lifetime of the beam (taubeam ): " << taubeam;1975 sortie << "Expected lifetime of the beam (taubeam [second]): " << taubeam; 1947 1976 1948 1977 sortie.close(); … … 1960 1989 1961 1990 1962 1991 /************************************************************************************** 1992 * 1993 * 1994 * 1995 * Fix the bugs to save ips, ip, is, it, itcla. By Jianfeng Zhang @ LAL, 16/04/2014. 1996 * 1997 ***************************************************************************************/ 1963 1998 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) 1964 1999 { … … 1967 2002 1968 2003 int NLostTotal(0); 1969 int irevloc , i0loc, imloc;2004 int irevloc(0), i0loc(0), imloc(0); 1970 2005 1971 2006 NLostTotal = asumhitcolli[asumhitcolli.size() - 1] + asumhits[asumhits.size() - 1]; … … 2100 2135 } 2101 2136 2137 //plot the particle hit the collimators with nsig <= 100. 2102 2138 ofstream outplot3; 2103 2139 string plot3; … … 2127 2163 2128 2164 for (int i(0); i < newnhitcolli.size(); ++i) { 2129 gnuout2 << "set label " << i + 1 << " " << '"' << lat.reseau[ipsnew[i]]->NAME << '"' << """ at " << i << ",13 rotate front" << endl;2165 gnuout2 << "set label " << i + 1 << " " << '"' << lat.reseau[ipsnew[i]]->NAME << '"' << """ at " << i+1 << ",13 rotate front" << endl; 2130 2166 } 2131 2167 gnuout2.close(); … … 2149 2185 2150 2186 //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 2187 //asumrem[i]: particle rest in the beam 2188 //asumhitcolli[i]: particle lost on the collimators 2189 //asumhits[i]: particle lost elsewhere 2151 2190 2152 2191 ofstream outplot4; … … 2554 2593 } 2555 2594 2556 2595 /*********************************************************************************************************** 2596 * 2597 * 2598 * 2599 * ********************************************************************************************************/ 2557 2600 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) 2558 2601 { -
ICOSIM/CPP/trunk/source/simulation.h
r5 r8 53 53 //method used to read the parameters in the collimator file 54 54 55 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);55 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); 56 56 57 57
Note: See TracChangeset
for help on using the changeset viewer.