[5] | 1 | #include <iostream> |
---|
| 2 | #include <fstream> |
---|
| 3 | #include <cmath> |
---|
| 4 | #include <iomanip> |
---|
[15] | 5 | #include <string.h> |
---|
| 6 | #include <stdio.h> |
---|
[5] | 7 | #include "simcrys.h" |
---|
| 8 | |
---|
[13] | 9 | //crystal routine of the proton beam |
---|
| 10 | #include "crystalproton.h" |
---|
| 11 | #include "crystalprotonfuns.h" |
---|
[5] | 12 | |
---|
[13] | 13 | extern "C" |
---|
| 14 | { |
---|
| 15 | struct { |
---|
| 16 | double Cry_length, Rcurv, C_xmax, C_ymax, Alayer, C_orient; |
---|
| 17 | }Par_Cry1_; |
---|
| 18 | }; |
---|
| 19 | |
---|
| 20 | |
---|
[5] | 21 | using namespace std; |
---|
| 22 | |
---|
[13] | 23 | //Fortran crystal routine for the proton beam |
---|
| 24 | // By Jianfeng Zhang @ LAL, 05/2014 |
---|
| 25 | // extern "C" { |
---|
| 26 | // // void cryst_(int,double,double,double,double,double,double); |
---|
| 27 | // void cryst_(int IS, double x, double xp, double y,double yp,double PC,double Length); |
---|
| 28 | // } |
---|
[5] | 29 | |
---|
| 30 | |
---|
[13] | 31 | |
---|
[5] | 32 | SimCrys::SimCrys(Crystal crys, Partcrys part) |
---|
| 33 | : crys(crys), part(part) //,moy(moy) //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 34 | { |
---|
| 35 | } |
---|
| 36 | |
---|
| 37 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 38 | ////////////////////////// FUNCTION BASES THAT MAKE SOME BASIC CALCULATION |
---|
| 39 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 40 | |
---|
| 41 | void SimCrys::bases(SimCrys& sim) |
---|
| 42 | { |
---|
| 43 | if ( (crys.miscut < 0 ) && ( part.x > 0) && (part.x < -crys.Cry_length * tan(crys.miscut))) { |
---|
| 44 | crys.L_chan = -part.x * tan(crys.miscut); |
---|
| 45 | } |
---|
[15] | 46 | crys.tchan = crys.L_chan / crys.Rcurv; |
---|
[5] | 47 | part.xp_rel = part.xp - crys.miscut; |
---|
| 48 | } |
---|
| 49 | |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 53 | ////////////////////////FUNCTIN CROSS TO KNOW IF THE PARTICULE CROSS THE CRYSTAL |
---|
| 54 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | bool SimCrys::cross(SimCrys& sim) |
---|
| 60 | { |
---|
[15] | 61 | bool crossing(true); |
---|
[5] | 62 | if ((part.y < crys.ymin) || (part.y > crys.ymax) || (part.x > crys.xmax)) { |
---|
| 63 | //cout << "OUT of the Crystal !!" << part.x << " , "<<part.xp << endl; |
---|
[15] | 64 | crossing = false; |
---|
[5] | 65 | part.x = part.x + part.xp * crys.s_length; |
---|
| 66 | |
---|
| 67 | part.y = part.y + part.yp * crys.s_length; |
---|
| 68 | |
---|
| 69 | part.nout = part.nout + 1; /////////////COUNT///////////// |
---|
| 70 | part.effet = 1; |
---|
| 71 | //cout<<" BLAAAAAAAAa"<<endl; |
---|
| 72 | |
---|
| 73 | } |
---|
| 74 | |
---|
[15] | 75 | return crossing; |
---|
[5] | 76 | } |
---|
| 77 | |
---|
| 78 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 79 | ////////////////////////////FUNCTION LAYER THAT DETERMIN IF THE PARTICULE IS IN THE AMORPHOUS LAYER///////////////////////////// |
---|
| 80 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 81 | |
---|
| 82 | bool SimCrys::layer(SimCrys& sim) |
---|
| 83 | { |
---|
| 84 | //cout<<"La fonction layer est utilisee>>>>>>>>>>>>>>"<<endl; |
---|
| 85 | bool out(true); |
---|
| 86 | if ((part.x < crys.Alayer) || ((part.y - crys.ymin) < crys.Alayer) || ((crys.ymax - part.y) < crys.Alayer)) { |
---|
[15] | 87 | cout<<"Amorphous layer! " << part.x <<" , " <<part.xp<<endl; |
---|
[5] | 88 | double a_eq; |
---|
| 89 | double b_eq; |
---|
| 90 | double c_eq; |
---|
| 91 | double delta; //a, b, c, Delta of the second order eq. |
---|
| 92 | double x0; |
---|
| 93 | double y0; |
---|
| 94 | double length_xs; //!Amorphous length |
---|
| 95 | double length_ys; //!Amorphous length |
---|
| 96 | double am_length; //!Amorphous length |
---|
| 97 | out = false; |
---|
| 98 | x0 = part.x; |
---|
| 99 | y0 = part.y; |
---|
| 100 | a_eq = (1 + part.xp * part.xp); |
---|
| 101 | b_eq = (2 * (part.x) * (part.xp) - 2 * (part.xp) * crys.Rcurv); |
---|
| 102 | c_eq = (part.x) * part.x - 2 * abs(part.x) * crys.Rcurv; |
---|
| 103 | delta = b_eq * b_eq - 4 * a_eq * c_eq; |
---|
| 104 | //cout<<"a : " << a_eq<<endl; //POUR TESTER FONCTION |
---|
| 105 | //cout<<"b: " << b_eq <<endl; //POUR TESTER FONCTION |
---|
| 106 | //cout<<"c : " << c_eq <<endl; //POUR TESTER FONCTION |
---|
| 107 | //cout<<"Delta : " << delta <<endl; //POUR TESTER FONCTION |
---|
| 108 | part.s = (-b_eq + sqrt(delta)) / ((2 * a_eq)) ; //in [m] |
---|
| 109 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
| 110 | if (part.s >= crys.s_length) { |
---|
| 111 | part.s = crys.s_length; |
---|
| 112 | } |
---|
| 113 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
| 114 | part.x = part.xp * part.s + x0 ; //in [m] |
---|
| 115 | //cout<<"x0 : " << x0 <<endl; //POUR TESTER FONCTION |
---|
| 116 | //cout<<"x : " << part.x <<endl; //POUR TESTER FONCTION |
---|
| 117 | |
---|
| 118 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
| 119 | //cout<<"Length_XS : " << length_xs <<endl; //POUR TESTER FONCTION |
---|
| 120 | |
---|
| 121 | if ((((part.yp >= 0) && ((part.y + part.yp * part.s) <= crys.ymax))) || (((part.yp < 0) && ((part.y + part.yp * part.s) >= crys.ymin)))) { |
---|
| 122 | length_ys = part.yp * length_xs; |
---|
| 123 | } else { |
---|
| 124 | part.s = (crys.ymax - part.y) / part.yp; |
---|
| 125 | length_ys = sqrt((crys.ymax - part.y) * (crys.ymax - part.y) + part.s * part.s); |
---|
| 126 | part.x = x0 + part.xp * part.s; |
---|
| 127 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | am_length = sqrt(length_xs * length_xs + length_ys * length_ys); |
---|
| 131 | //cout << "AM_Length :" << am_length << endl; // TESTER FONCTION |
---|
| 132 | part.s = part.s / 2; |
---|
| 133 | part.x = x0 + part.xp * part.s; |
---|
| 134 | part.y = y0 + part.yp * part.s; |
---|
| 135 | //cout<<" am_length "<<am_length<<endl; |
---|
| 136 | move_am(crys.IS, crys.NAM, am_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
| 137 | //cout<<"AM LAYER !!!! "<<part.PC<<endl; |
---|
| 138 | //cout<< "Amorphous"<< endl; |
---|
| 139 | part.x = part.x + part.xp * (crys.Cry_length - part.s); |
---|
| 140 | part.y = part.y + part.yp * (crys.Cry_length - part.s); |
---|
| 141 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 142 | part.effet = 2; |
---|
| 143 | } else if ((part.x > (crys.xmax - crys.Alayer)) && (part.x < crys.xmax)) { |
---|
| 144 | out = false; |
---|
| 145 | |
---|
| 146 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
| 147 | |
---|
| 148 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 149 | part.effet = 2; |
---|
| 150 | //cout<<"Fix HERE"<<endl; |
---|
| 151 | } |
---|
| 152 | return out; |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 157 | //SI LES DEUX FONCTIONS D AVT RETURNE TRUE, ALORS ON CALCULE LES PARA CRITIQUE ET D AUTRE PARA (POINT C - E) |
---|
| 158 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 159 | |
---|
| 160 | |
---|
| 161 | |
---|
| 162 | void SimCrys::parameters(SimCrys& sim) |
---|
| 163 | { |
---|
| 164 | //THE FIRST PART IS FOR THE (110) ORIENTATION |
---|
| 165 | //cout << "The parameters function is use !!!!! "<< endl; |
---|
| 166 | double c_v1; |
---|
| 167 | double c_v2; |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | crys.xpcrit0 = sqrt(2.10e-9 * crys.eUm[crys.IS] / part.PC); |
---|
| 171 | crys.Rcrit = part.PC / (2.e-6 * crys.eUm[crys.IS]) * crys.AI[crys.IS]; //if R>Rcritical=>no channeling is possible |
---|
| 172 | crys.ratio = crys.Rcurv / crys.Rcrit; |
---|
| 173 | crys.xpcrit = crys.xpcrit0 * (crys.Rcurv - crys.Rcrit) / crys.Rcurv; |
---|
| 174 | c_v1 = 1.7; |
---|
| 175 | c_v2 = -1.5; |
---|
| 176 | if (crys.ratio <= 1) { |
---|
| 177 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.4 * crys.ratio); //rms scattering |
---|
| 178 | part.Ang_avr = c_v2 * crys.xpcrit0 * 0.05 * crys.ratio; //average angle reflection |
---|
| 179 | part.Vcapt = 0.0; |
---|
| 180 | } else if (crys.ratio <= 3) { |
---|
| 181 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.571 * 0.3 * crys.ratio + 0.85); //rms scattering |
---|
| 182 | part.Ang_avr = c_v2 * crys.xpcrit0 * (0.1972 * crys.ratio - 0.1472); //average angle reflection |
---|
| 183 | part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2); //K=0.00070 is taken based on simulations using CATCH.f (V.Biryukov) |
---|
| 184 | } else { |
---|
| 185 | part.Ang_rms = c_v1 * crys.xpcrit0 * (1 / crys.ratio); |
---|
| 186 | part.Ang_avr = c_v2 * crys.xpcrit0 * (1 - 1.6667 / crys.ratio); |
---|
| 187 | part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2); |
---|
| 188 | } |
---|
| 189 | //cout <<crys.xpcrit<<endl; |
---|
| 190 | if (crys.C_orient == 2) { |
---|
| 191 | part.Ang_avr = part.Ang_avr * 0.93; |
---|
| 192 | part.Ang_rms = part.Ang_rms * 1.05; // FOR (111) ORIENTATION |
---|
| 193 | crys.xpcrit = crys.xpcrit * 0.98; |
---|
| 194 | } |
---|
| 195 | //cout <<crys.xpcrit<<endl; |
---|
| 196 | } |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 202 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE I CF SHEMA!!! C EST LA PARTIE CHANNELING D OU LE NOM |
---|
| 203 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | |
---|
| 208 | void SimCrys::channel(SimCrys& sim) |
---|
| 209 | { |
---|
| 210 | //cout << "La fonction channel est utilisee!!!! "<<endl; |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | part.Chann = sqrt(crys.xpcrit * crys.xpcrit - part.xp_rel * part.xp_rel) / crys.xpcrit0; //probability of channeling/volume capture |
---|
| 214 | |
---|
| 215 | //cout <<"Channeling Prob :" << part.Chann <<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
| 216 | //cout <<"Channeling ANGLE :" << crys.tchan<<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | double n_atom(0.1); //probability for entering channeling near atomic planes |
---|
| 220 | double Dxp(0); //change angle from channeling [mrad] |
---|
| 221 | double TLdech1(0); //tipical dechanneling length(1) [m] |
---|
| 222 | double tdech(0); |
---|
| 223 | double Ldech(0); //Dechanneling. length |
---|
| 224 | double Sdech(0); |
---|
| 225 | //cout<<"Chann : "<<part.Chann<<endl; |
---|
| 226 | if (random() <= part.Chann) { |
---|
| 227 | TLdech1 = 0.0005 * part.PC * pow(1 - 1 / crys.ratio, 2); // calculate tipical dechanneling length(m) |
---|
| 228 | if (random() <= n_atom) { |
---|
| 229 | TLdech1 = 0.000002 * part.PC * pow(1. - 1. / crys.ratio, 2); // dechanneling length (m) for short crystal for high amplitude particles |
---|
| 230 | } |
---|
| 231 | //cout<<"@ energy " << part.PC <<" dechanneling length "<< TLdech1 << endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
| 232 | part.Dechan = log(random()); |
---|
| 233 | //cout<<"dechannn " << part.Dechan <<endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
| 234 | Ldech = -TLdech1 * part.Dechan; //careful: the dechanneling lentgh is along the trajectory of the particle -not along the longitudinal coordinate... |
---|
| 235 | tdech = Ldech / crys.Rcurv; |
---|
| 236 | Sdech = Ldech * cos(crys.miscut + 0.5 * tdech); |
---|
| 237 | // 2 option if the particule channeling !!! |
---|
| 238 | // (1) dechanneling |
---|
| 239 | // (2) full channeling |
---|
| 240 | |
---|
| 241 | // (1) : |
---|
| 242 | |
---|
| 243 | if (Ldech < crys.L_chan) { |
---|
| 244 | //cout <<"Dechanneling ! "<< endl; |
---|
| 245 | Dxp = Ldech / crys.Rcurv; //change angle from channeling [mrad] |
---|
| 246 | |
---|
| 247 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
| 248 | part.effet = 4; |
---|
| 249 | part.x = part.x + Ldech * (sin(0.5 * Dxp + crys.miscut)); // trajectory at channeling exit |
---|
| 250 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
| 251 | part.y = part.y + part.yp * Sdech; |
---|
| 252 | |
---|
| 253 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
| 254 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | move_am(crys.IS, crys.NAM, crys.s_length - Sdech, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
| 258 | |
---|
| 259 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * part.y; //energy loss to ionization [GeV] |
---|
| 260 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
| 261 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | // (2) : |
---|
| 265 | |
---|
| 266 | else { |
---|
| 267 | //cout <<"Channeling ! "<< endl; |
---|
| 268 | Dxp = crys.L_chan / crys.Rcurv; //change angle [mrad] |
---|
| 269 | part.x = part.x + crys.L_chan * (sin(0.5 * Dxp + crys.miscut)); //trajectory at channeling exit |
---|
| 270 | |
---|
| 271 | |
---|
| 272 | // Removed dependence on incoming angle according to discussion with I.Yazynin on 24/10/2012 |
---|
| 273 | //part.xp = part.xp + Dxp + (random_gauss()-0.5)*crys.xpcrit*0.5; // [mrad] |
---|
| 274 | part.xp = /*part.xp*/ + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; // [mrad] |
---|
| 275 | |
---|
| 276 | part.y = part.y + crys.s_length * part.yp; |
---|
| 277 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * crys.Cry_length; // energy loss to ionization [GeV] |
---|
| 278 | part.nchann = part.nchann + 1; /////////////COUNT///////////// |
---|
| 279 | part.effet = 3; |
---|
| 280 | } |
---|
| 281 | } else { |
---|
| 282 | //cout <<"Volume Reflection ! "<< endl; |
---|
| 283 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
| 284 | part.effet = 5; |
---|
| 285 | Dxp = 0.5 * (part.xp_rel / crys.xpcrit + 1) * part.Ang_avr; |
---|
| 286 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
| 287 | /* |
---|
| 288 | next line new from sasha |
---|
| 289 | part.xp=part.xp+0.45*(part.xp/crys.xpcrit+1)*part.Ang_avr; |
---|
| 290 | valentina fix here |
---|
| 291 | */ |
---|
| 292 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 293 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 294 | |
---|
| 295 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
| 296 | |
---|
| 297 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 298 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 299 | |
---|
| 300 | } |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | |
---|
| 305 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 306 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE II CF SHEMA!!! C EST LA PARTIE VOLUME REFLEXION D OU LE NOM |
---|
| 307 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 308 | |
---|
| 309 | |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | void SimCrys::reflection(SimCrys& sim) |
---|
| 313 | { |
---|
| 314 | //cout<<"The function reflection is used" <<endl; |
---|
| 315 | double Dxp(0); //change angle from channeling [mrad] |
---|
| 316 | double Lrefl(0); //distance of the reflection point [m] |
---|
| 317 | double Srefl(0); //distance of the reflection point [m] |
---|
| 318 | double TLdech2(0); //tipical dechanneling length(2) [m] |
---|
| 319 | double tdech(0); |
---|
| 320 | double Ldech(0); //Dechanneling. length |
---|
| 321 | double Sdech(0); |
---|
| 322 | double Red_S(0); //Reduced length (in case of dechanneling) |
---|
| 323 | double Rlength(0); //Reduced length |
---|
| 324 | |
---|
| 325 | Lrefl = (part.xp_rel) * crys.Rcurv; |
---|
| 326 | Srefl = sin(part.xp) * Lrefl; |
---|
| 327 | |
---|
| 328 | |
---|
| 329 | if ((Lrefl > 0) && (Lrefl < crys.Cry_length)) { // IE : If the VR point is inside the crystal!!!! |
---|
| 330 | |
---|
| 331 | //cout<<"VRCapt : "<<part.Vcapt<<endl; |
---|
| 332 | if ((random() > part.Vcapt) || (crys.ZN == 0)) { |
---|
| 333 | //cout<< "Volume Reflexion! "<<endl; |
---|
| 334 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
| 335 | part.effet = 5; |
---|
| 336 | part.x = part.x + part.xp * Srefl; |
---|
| 337 | part.y = part.y + part.yp * Srefl; |
---|
| 338 | Dxp = part.Ang_avr; |
---|
| 339 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
| 340 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
| 341 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
| 342 | |
---|
| 343 | if (crys.ZN > 0) { |
---|
| 344 | move_am(crys.IS, crys.NAM, crys.s_length - Srefl, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
| 345 | } |
---|
| 346 | |
---|
| 347 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
| 348 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
| 349 | } else { |
---|
| 350 | //cout<< "Volume Capture! "; |
---|
| 351 | part.nvcapt = part.nvcapt + 1; /////////////COUNT///////////// |
---|
| 352 | part.effet = 6; |
---|
| 353 | part.x = part.x + part.xp * Srefl; |
---|
| 354 | part.y = part.y + part.yp * Srefl; |
---|
| 355 | /* |
---|
| 356 | TLdech2= 0.00011*pow(part.PC,0.25)*pow((1.-1./crys.ratio),2) //dechanneling length [m] |
---|
| 357 | part.Dechan = log(1.-random()) |
---|
| 358 | Ldech = -TLdech2*part.Dechan |
---|
| 359 | next 2 lines new from sasha - different dechanneling probability |
---|
| 360 | */ |
---|
| 361 | TLdech2 = 0.01 * part.PC * pow((1. - 1. / crys.ratio), 2); // typical dechanneling length [m] |
---|
| 362 | Ldech = 0.005 * TLdech2 * pow((sqrt(0.01 - log(random())) - 0.1), 2); // (dechanneling length) [m] |
---|
| 363 | tdech = Ldech / crys.Rcurv; |
---|
| 364 | Sdech = Ldech * cos(part.xp + 0.5 * tdech); |
---|
| 365 | |
---|
| 366 | |
---|
| 367 | |
---|
| 368 | if (Ldech < (crys.Cry_length - Lrefl)) { // for dechanneling part |
---|
| 369 | //cout<< "Dechanneling! "; |
---|
| 370 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
| 371 | part.effet = 4; |
---|
| 372 | Dxp = Ldech / crys.Rcurv; |
---|
| 373 | part.x = part.x + Ldech * (sin(0.5 * Dxp + part.xp)); //trajectory at channeling exit |
---|
| 374 | part.y = part.y + Sdech * part.yp; |
---|
| 375 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
| 376 | Red_S = crys.s_length - Srefl - Sdech; |
---|
| 377 | // move in amorphous substance---------------------------- |
---|
| 378 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
| 379 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
| 380 | |
---|
| 381 | if (crys.ZN > 0) { |
---|
| 382 | //cout<<"Amorphous! "<<endl; /////////////COUNT///////////// |
---|
| 383 | move_am(crys.IS, crys.NAM, Red_S, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
| 387 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
| 388 | //cout<<endl; |
---|
| 389 | } else { |
---|
| 390 | //cout<< "Volume Capture! "<<endl; |
---|
| 391 | //Pas besoin d additionner 1 dans le compte de vol capt car deja fait!!! |
---|
| 392 | Rlength = crys.Cry_length - Lrefl; |
---|
| 393 | crys.tchan = Rlength / crys.Rcurv; |
---|
| 394 | Red_S = Rlength * cos(part.xp + 0.5 * crys.tchan); |
---|
| 395 | Dxp = (crys.Cry_length - Lrefl) / crys.Rcurv; |
---|
| 396 | part.x = part.x + sin(0.5 * Dxp + part.xp) * Rlength; // trajectory at channeling exit |
---|
| 397 | part.y = part.y + Red_S * part.yp; |
---|
| 398 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; //[mrad] |
---|
| 399 | } |
---|
| 400 | } |
---|
| 401 | } |
---|
| 402 | // move in amorphous substance for big input angles--------------- |
---|
| 403 | else { |
---|
| 404 | //cout<<"Amorphous! "<<endl; |
---|
| 405 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 406 | part.effet = 2; |
---|
| 407 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 408 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 409 | |
---|
| 410 | if (crys.ZN > 0) { |
---|
| 411 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
| 412 | } |
---|
| 413 | |
---|
| 414 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 415 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 416 | } |
---|
| 417 | |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | |
---|
| 422 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 423 | //THE FUNCTION THAT DO ALL THE SCHEMA USING THE PREVIOUS FUNCTION !!!!!!!!!!111 |
---|
| 424 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 425 | |
---|
| 426 | |
---|
| 427 | |
---|
| 428 | |
---|
[15] | 429 | //called the FORTRAN cyrstal-particle routine from SixTrack |
---|
| 430 | /* |
---|
| 431 | void SimCrys::general(SimCrys& sim, const int& pas, char name, double emitx0, double emity0, double enom, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt) |
---|
| 432 | { |
---|
| 433 | |
---|
| 434 | long int max_npart = 1500000; |
---|
| 435 | |
---|
| 436 | char name_coll[50] = name; |
---|
| 437 | |
---|
| 438 | int NP = ; |
---|
| 439 | double ENOM = enom*1e-3; // [GeV] |
---|
| 440 | double EMITX0 = emitx0; |
---|
| 441 | double EMITY0 = emity0; |
---|
| 442 | double AX = , AY = , BX = , BY = ; |
---|
| 443 | |
---|
| 444 | double X_IN[max_npart] ; |
---|
| 445 | double XP_IN[max_npart]; |
---|
| 446 | double Y_IN[max_npart]; |
---|
| 447 | double YP_IN[max_npart]; |
---|
| 448 | double P_IN[max_npart]; |
---|
| 449 | double S_IN[max_npart] = ; |
---|
| 450 | |
---|
| 451 | int flagsec[max_npart] = ; |
---|
| 452 | bool dowrite_impact = true; |
---|
| 453 | int name[max_npart] = ; |
---|
| 454 | |
---|
| 455 | char C_MATERIAL[6]; |
---|
| 456 | double C_LENGTH = crys.s_length; |
---|
| 457 | double C_ROTATION = C_rotation; |
---|
| 458 | double C_APERTURE = C_aperture; |
---|
| 459 | double C_OFFSET = C_offset; |
---|
| 460 | double C_TILT[2] = {C_tilt, 0} ; |
---|
| 461 | |
---|
| 462 | int LHIT[max_npart] = ; //position of the particle been hit |
---|
| 463 | int PART_ABS[max_npart] = ; |
---|
| 464 | |
---|
| 465 | double IMPACT[max_npart] = ; |
---|
| 466 | double INDIV[max_npart] = ; |
---|
| 467 | double LINT[max_npart] = ; |
---|
| 468 | |
---|
| 469 | */ |
---|
| 470 | /* ----------------- |
---|
| 471 | !++ Copy particle data to 1-dim array and go back to meters |
---|
| 472 | */ |
---|
| 473 | /* |
---|
| 474 | for(j = 1, napx) |
---|
| 475 | rcx(j) = (xv(1,j)-torbx(ie))/1d3 |
---|
| 476 | rcxp(j) = (yv(1,j)-torbxp(ie))/1d3 |
---|
| 477 | rcy(j) = (xv(2,j)-torby(ie))/1d3 |
---|
| 478 | rcyp(j) = (yv(2,j)-torbyp(ie))/1d3 |
---|
| 479 | rcp(j) = ejv(j)/1d3 |
---|
| 480 | rcs(j) = 0d0 |
---|
| 481 | part_hit_before(j) = part_hit(j) |
---|
| 482 | rcx0(j) = rcx(j) |
---|
| 483 | rcxp0(j) = rcxp(j) |
---|
| 484 | rcy0(j) = rcy(j) |
---|
| 485 | rcyp0(j) = rcyp(j) |
---|
| 486 | rcp0(j) = rcp(j) |
---|
| 487 | ejf0v(j) = ejfv(j) |
---|
| 488 | ! |
---|
| 489 | !++ For zero length element track back half collimator length |
---|
| 490 | ! |
---|
| 491 | if (stracki.eq.0.) then |
---|
| 492 | rcx(j) = rcx(j) - 0.5d0*c_length*rcxp(j) |
---|
| 493 | rcy(j)= rcy(j) - 0.5d0*c_length*rcyp(j) |
---|
| 494 | else |
---|
| 495 | Write(*,*) |
---|
| 496 | & "ERROR: Non-zero length collimator!" |
---|
| 497 | STOP |
---|
| 498 | endif |
---|
| 499 | flukaname(j) = ipart(j)+100*samplenumber |
---|
| 500 | ! |
---|
| 501 | end do |
---|
| 502 | ! |
---|
| 503 | !++ Do the collimation tracking |
---|
| 504 | ! |
---|
| 505 | enom_gev = myenom*1d-3 |
---|
| 506 | ! |
---|
| 507 | !++ Allow |
---|
[5] | 508 | |
---|
[15] | 509 | |
---|
| 510 | |
---|
| 511 | |
---|
| 512 | if (crys.IS == 0) { |
---|
| 513 | C_MATERIAL = "CRY-Si"; |
---|
| 514 | } else if (crys.IS == 1) { |
---|
| 515 | C_MATERIAL = "CRY-W"; |
---|
| 516 | } else if (crys.IS == 2) { |
---|
| 517 | C_MATERIAL = "CRY-C"; |
---|
| 518 | } else if (crys.IS == 3) { |
---|
| 519 | C_MATERIAL = "CRY-Ge"; |
---|
| 520 | }else{ |
---|
| 521 | cerr << "Material of the Crystal not fund !!!" << endl; |
---|
| 522 | } |
---|
| 523 | |
---|
| 524 | |
---|
| 525 | |
---|
| 526 | int count = 0; |
---|
| 527 | string rest; |
---|
| 528 | ifstream entree; |
---|
| 529 | string name2; |
---|
| 530 | name2 = outputpath + "/ascii_before.csv"; |
---|
| 531 | entree.open(name2.c_str()); |
---|
| 532 | if (entree) { |
---|
| 533 | while (!entree.eof()) { |
---|
| 534 | count = count + 1; |
---|
| 535 | |
---|
| 536 | getline(entree, rest, ','); |
---|
| 537 | X_IN[count] = atof(rest.c_str()); |
---|
| 538 | //cout<<" X " << part.x<<endl; |
---|
| 539 | |
---|
| 540 | getline(entree, rest, ','); |
---|
| 541 | Y_IN[count] = atof(rest.c_str()); |
---|
| 542 | //cout<<" Y " << part.y<<endl; |
---|
| 543 | |
---|
| 544 | getline(entree, rest, ','); |
---|
| 545 | XP_IN[count] = atof(rest.c_str()); |
---|
| 546 | //cout<<" XP " << part.xp<<endl; |
---|
| 547 | |
---|
| 548 | getline(entree, rest, ','); |
---|
| 549 | YP_IN[count] = atof(rest.c_str()); |
---|
| 550 | //cout<<" YP " << part.yp<<endl; |
---|
| 551 | |
---|
| 552 | getline(entree, rest); |
---|
| 553 | P_IN[count] = atof(rest.c_str()); |
---|
| 554 | |
---|
| 555 | } |
---|
| 556 | } |
---|
| 557 | entree.close(); |
---|
| 558 | |
---|
| 559 | |
---|
| 560 | //from the file "crystaltrack_dan.f", 08/2014, by Jianfeng Zhang @ LAL |
---|
| 561 | COLLIMATE_CRY_(name_coll,C_MATERIAL, C_LENGTH, |
---|
| 562 | C_ROTATION, |
---|
| 563 | C_APERTURE, C_OFFSET, C_TILT, |
---|
| 564 | X_IN, XP_IN, Y_IN, |
---|
| 565 | YP_IN, P_IN, S_IN, NP, ENOM, LHIT, |
---|
| 566 | PART_ABS, IMPACT, INDIV, LINT, |
---|
| 567 | BX,BY,AX, |
---|
| 568 | AY,EMITX0,EMITY0, |
---|
| 569 | name,flagsec,dowrite_impact); |
---|
| 570 | } |
---|
| 571 | |
---|
| 572 | */ |
---|
| 573 | |
---|
| 574 | /******************************************************************************** |
---|
| 575 | * |
---|
| 576 | * original crystal routine |
---|
| 577 | * *******************************************************************************/ |
---|
[17] | 578 | void SimCrys::general(SimCrys& sim, const int& pas, long int& nhit, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt) |
---|
[5] | 579 | { |
---|
| 580 | srand ( time(NULL) ); |
---|
| 581 | double AV_xp0(0); |
---|
| 582 | double AV_yp0(0); |
---|
| 583 | double RMSxp0(0); |
---|
| 584 | double RMSyp0(0); |
---|
| 585 | double AV_xp1(0); |
---|
| 586 | double AV_yp1(0); |
---|
| 587 | double long RMSxp1(0); |
---|
| 588 | double long RMSyp1(0); |
---|
| 589 | double av_pc1(0); |
---|
| 590 | double av_pc0(0); |
---|
| 591 | double SIGxp0(0); |
---|
| 592 | double SIGyp0(0); |
---|
| 593 | double SIGxp1(0); |
---|
| 594 | double SIGyp1(0); |
---|
| 595 | double xp0(0); |
---|
| 596 | double xp1(0); |
---|
| 597 | double count(0); |
---|
| 598 | double zpos(-0.1); |
---|
| 599 | double xfin(0); |
---|
| 600 | double yfin(0); |
---|
| 601 | int const col(15); |
---|
| 602 | |
---|
[15] | 603 | |
---|
| 604 | //parameter to check whether the particle hit the crystal |
---|
| 605 | char Proc[50] = " "; // process of the particle-crystal interaction |
---|
| 606 | bool crossing = false; //the particle hit the collimator or not |
---|
| 607 | int layerflag = 1; //check the particle exit the layer in "AM" |
---|
[5] | 608 | //-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------ |
---|
[8] | 609 | int mat=0; |
---|
[5] | 610 | |
---|
[8] | 611 | double p0=0.0; |
---|
| 612 | double zlm=0.0; |
---|
| 613 | double shift=0.0; |
---|
| 614 | double x_shift=0.0, xp_shift=0.0, s_shift=0.0; //coordinates after shift/rotation |
---|
| 615 | double x_rot=0.0, xp_rot=0.0, s_rot=0.0; |
---|
| 616 | 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 |
---|
| 617 | 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) |
---|
| 618 | double x00=0.0; |
---|
| 619 | double z=0.0; |
---|
| 620 | double z00=0.0; |
---|
| 621 | double zp=0.0; |
---|
| 622 | double dpop=0.0; |
---|
| 623 | double s=0.0; |
---|
| 624 | double a_eq=0.0, b_eq=0.0, c_eq=0.0, Delta=0.0; |
---|
[15] | 625 | int part_abs=-1; |
---|
[8] | 626 | double x=0.0; |
---|
| 627 | double xp=0.0; |
---|
| 628 | double y=0.0; |
---|
| 629 | double yp=0.0; |
---|
| 630 | double p=0.0; //be careful: [Gev] |
---|
| 631 | double s_in=0.0; |
---|
[5] | 632 | |
---|
| 633 | // adding variables for the pencil beam. Variables in the absolute reference frame. |
---|
| 634 | |
---|
[8] | 635 | double x_in0=0.0; |
---|
| 636 | double xp_in0=0.0; |
---|
| 637 | double y_in0=0.0; |
---|
| 638 | double yp_in0=0.0; |
---|
| 639 | double p_in0=0.0; //be careful: [Gev] |
---|
| 640 | double s_in0=0.0; |
---|
| 641 | double s_impact=0.0; |
---|
[5] | 642 | double totcount(0); |
---|
| 643 | |
---|
[9] | 644 | double mirror=Mirror; |
---|
[8] | 645 | double tiltangle=0.0; |
---|
[5] | 646 | |
---|
| 647 | |
---|
[8] | 648 | double c_length=0.0; //Length in m |
---|
[5] | 649 | double c_rotation(C_rotation) ; //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED |
---|
| 650 | double c_aperture(C_aperture) ; //Aperture in m |
---|
| 651 | double c_offset(C_offset) ; //Offset in m |
---|
| 652 | double c_tilt[2] = {C_tilt, 0} ; //Tilt in radian |
---|
| 653 | double c_tilt0[2] ; //Tilt in radian |
---|
[15] | 654 | |
---|
[5] | 655 | |
---|
[8] | 656 | double xp_tangent=0.0; |
---|
[5] | 657 | |
---|
| 658 | double Cry_tilt(Crystal_tilt); //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++//////////////// |
---|
| 659 | |
---|
| 660 | |
---|
[17] | 661 | |
---|
| 662 | |
---|
| 663 | |
---|
[5] | 664 | //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS------------------------------------------------------------------------- |
---|
| 665 | c_length = crys.Cry_length; |
---|
| 666 | |
---|
| 667 | if (crys.IS == 0) { |
---|
[15] | 668 | mat = 8; |
---|
[5] | 669 | } else if (crys.IS == 1) { |
---|
[15] | 670 | mat = 9; |
---|
[5] | 671 | } else if (crys.IS == 2) { |
---|
[15] | 672 | mat = 10; |
---|
[5] | 673 | } else if (crys.IS == 3) { |
---|
[15] | 674 | mat = 11; |
---|
[5] | 675 | } else { |
---|
[15] | 676 | cout << "Material of the Crystal not fund !!!" << endl; |
---|
[5] | 677 | } |
---|
| 678 | |
---|
[15] | 679 | |
---|
[5] | 680 | if (c_length > 0 ) { |
---|
[15] | 681 | p0 = part.PC; |
---|
| 682 | c_tilt0[0] = c_tilt[0]; |
---|
| 683 | c_tilt0[1] = c_tilt[1]; |
---|
| 684 | tiltangle = c_tilt0[0]; |
---|
[5] | 685 | } |
---|
[15] | 686 | |
---|
| 687 | |
---|
| 688 | //new.......................... |
---|
| 689 | // call scatin(p0) ?????????????? |
---|
[5] | 690 | |
---|
[15] | 691 | // scatin_(p0); |
---|
| 692 | |
---|
| 693 | /* EVENT LOOP, initial distribution is here a flat distribution with |
---|
| 694 | * xmin=x-, xmax=x+, etc. from the input file |
---|
| 695 | */ |
---|
[17] | 696 | |
---|
| 697 | |
---|
| 698 | //initialize particle variables |
---|
| 699 | part.n_part = 0; |
---|
| 700 | part.nabsorbed = 0; |
---|
| 701 | part.namor = 0; |
---|
| 702 | part.nchann = 0; |
---|
| 703 | part.ndechann = 0; |
---|
| 704 | part.nout =0; |
---|
| 705 | part.nvcapt = 0; |
---|
| 706 | part.nvrefl = 0; |
---|
| 707 | part.nhit = 0; |
---|
| 708 | nhit = 0; |
---|
[15] | 709 | double fracab = 0; |
---|
| 710 | // int n_chan = 0; // valentina :initialize to zero the counters for crystal effects |
---|
| 711 | // int n_VR = 0; |
---|
| 712 | // int n_amorphous = 0; |
---|
| 713 | |
---|
| 714 | |
---|
| 715 | double impact = 0; |
---|
| 716 | double indiv = 0; |
---|
| 717 | double lint = 0; |
---|
| 718 | |
---|
| 719 | |
---|
| 720 | // c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
---|
| 721 | //c WRITE(*,*) ITURN,ICOLL |
---|
| 722 | /* do j = 1, nev |
---|
| 723 | //c |
---|
| 724 | impact(j) = -1.0; |
---|
| 725 | lint(j) = -1.0; |
---|
| 726 | indiv(j) = -1.0; |
---|
| 727 | |
---|
[5] | 728 | |
---|
[15] | 729 | idx_proc = name(j)-100*samplenumber; |
---|
[5] | 730 | |
---|
[15] | 731 | if (ITURN == 1){ // !daniele |
---|
| 732 | bool_proc_old(idx_proc)=-1;} // !daniele |
---|
| 733 | // elseif (bool_proc(idx_proc) != -1){ |
---|
| 734 | else{ |
---|
| 735 | bool_proc_old(idx_proc)=bool_proc(idx_proc); // !daniele |
---|
| 736 | } // !daniele |
---|
| 737 | |
---|
| 738 | PROC= "out"; // !the default process is 'out' !daniele |
---|
| 739 | bool_proc(idx_proc)=-1; // !daniele |
---|
| 740 | */ |
---|
| 741 | |
---|
| 742 | int nabs = 0; |
---|
| 743 | |
---|
[5] | 744 | //ofstream sortie("results.csv"); |
---|
| 745 | //sortie << "xp1" <<","<< "xp0"<<","<< "x"<<","<< "y"<<endl; |
---|
| 746 | ofstream sortieIco; |
---|
| 747 | string name1; |
---|
[15] | 748 | |
---|
| 749 | //file with the particle coordinate after the crystal |
---|
[5] | 750 | name1 = outputpath + "/ascii_after.csv"; |
---|
| 751 | |
---|
| 752 | sortieIco.open(name1.c_str()); |
---|
| 753 | |
---|
| 754 | string rest; |
---|
| 755 | |
---|
| 756 | ifstream entree; |
---|
| 757 | string name2; |
---|
| 758 | name2 = outputpath + "/ascii_before.csv"; |
---|
| 759 | |
---|
| 760 | entree.open(name2.c_str()); |
---|
| 761 | if (entree) { |
---|
[15] | 762 | while (!entree.eof()) { |
---|
| 763 | count = count + 1; |
---|
| 764 | |
---|
[17] | 765 | //particle hit the crystal or not |
---|
| 766 | bool hitflag = false; |
---|
| 767 | |
---|
| 768 | |
---|
[15] | 769 | getline(entree, rest, ','); |
---|
| 770 | part.x = atof(rest.c_str()); |
---|
| 771 | //cout<<" X " << part.x<<endl; |
---|
[5] | 772 | |
---|
[15] | 773 | getline(entree, rest, ','); |
---|
| 774 | part.y = atof(rest.c_str()); |
---|
| 775 | //cout<<" Y " << part.y<<endl; |
---|
[5] | 776 | |
---|
[15] | 777 | getline(entree, rest, ','); |
---|
| 778 | part.xp = atof(rest.c_str()); |
---|
| 779 | //cout<<" XP " << part.xp<<endl; |
---|
[5] | 780 | |
---|
[15] | 781 | getline(entree, rest, ','); |
---|
| 782 | part.yp = atof(rest.c_str()); |
---|
| 783 | //cout<<" YP " << part.yp<<endl; |
---|
[5] | 784 | |
---|
[15] | 785 | getline(entree, rest); |
---|
| 786 | part.PC = atof(rest.c_str()); |
---|
| 787 | |
---|
| 788 | if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT |
---|
[5] | 789 | |
---|
[15] | 790 | av_pc0 = av_pc0 + part.PC; |
---|
| 791 | |
---|
| 792 | //part.x=random_gauss()*0.00115; |
---|
| 793 | //part.xp=random_gauss()*sqrt(5*1e-6); |
---|
| 794 | //cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl; |
---|
[5] | 795 | |
---|
[15] | 796 | // IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 |
---|
| 797 | //bdelta=tan(part.xp)*10; |
---|
| 798 | part.x = (part.x); //-0.0045; |
---|
| 799 | |
---|
| 800 | |
---|
| 801 | //bdelta=tan(part.yp)*10; |
---|
| 802 | part.y = (part.y); //-0.0045; |
---|
[5] | 803 | |
---|
[15] | 804 | //cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl; |
---|
| 805 | xp0 = part.xp; |
---|
[5] | 806 | |
---|
[15] | 807 | //part.y=random_gauss()*sqrt(0.00095); |
---|
| 808 | //part.yp=5*1e-6; |
---|
[5] | 809 | |
---|
[15] | 810 | //part.PC=120.976; |
---|
| 811 | //cout<<" X :"<<part.x<<" XP :"<<part.xp<<" Y :"<<part.y<<" YP :"<<part.yp<<" PC :"<<part.PC<<endl; |
---|
| 812 | AV_xp0 = AV_xp0 + part.xp; //average x angle before impact |
---|
| 813 | AV_yp0 = AV_yp0 + part.yp; //average y angle before impact |
---|
| 814 | RMSxp0 = RMSxp0 + part.xp * part.xp; //rms x angle before impact |
---|
| 815 | RMSyp0 = RMSyp0 + part.yp * part.yp; //rms y angle before impact |
---|
[5] | 816 | |
---|
[15] | 817 | // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 818 | //change the beam frame to the crystal frame !!!!!!!!! |
---|
| 819 | // since the interaction between the particle and the crystal is observed in the |
---|
| 820 | // crystal reference system |
---|
[5] | 821 | |
---|
[15] | 822 | //initial coordinates in the beam frame |
---|
| 823 | s = 0; |
---|
| 824 | x = part.x; |
---|
| 825 | xp = part.xp; |
---|
| 826 | z = part.y; |
---|
| 827 | zp = part.yp; |
---|
| 828 | p = part.PC; |
---|
[5] | 829 | |
---|
[15] | 830 | x_temp = 0; |
---|
| 831 | x_int = 0; |
---|
| 832 | x_rot = 0; |
---|
| 833 | x_shift = 0; |
---|
| 834 | s_temp = 0; |
---|
| 835 | s_int = 0; |
---|
| 836 | s_rot = 0; |
---|
| 837 | s_shift = 0; |
---|
| 838 | xp_int = 0; |
---|
| 839 | xp_temp = 0; |
---|
| 840 | xp_rot = 0; |
---|
| 841 | xp_shift = 0; |
---|
| 842 | shift = 0; |
---|
| 843 | tilt_int = 0; |
---|
| 844 | dpop = (p - p0) / p0; //momentum spread |
---|
| 845 | |
---|
| 846 | |
---|
| 847 | /*---------------DANIELE------------------------- |
---|
| 848 | ! corrected position for variable association for FirstImpact.dat |
---|
| 849 | ! also for a vertical crystal case. |
---|
| 850 | ! Only x_in0 (i.e. b) have to be assigned after the change of reference frame |
---|
| 851 | !-----------------------------------------------*/ |
---|
[5] | 852 | |
---|
[15] | 853 | s_in0 = s_in; //daniele |
---|
| 854 | // x_in0(j) = x //!daniele |
---|
| 855 | xp_in0 = xp; //daniele |
---|
| 856 | y_in0 = z; //daniele |
---|
| 857 | yp_in0 = zp; //daniele |
---|
| 858 | p_in0 = p; //daniele |
---|
| 859 | |
---|
| 860 | /*----------------DANIELE------------------*/ |
---|
[5] | 861 | |
---|
[13] | 862 | |
---|
[15] | 863 | //********************************************************************************** |
---|
| 864 | |
---|
| 865 | // ------------------------------transform particle coordinates to get into collimator coordinate system ----------------------------------------- |
---|
| 866 | |
---|
[5] | 867 | |
---|
[15] | 868 | // first check whether particle was lost before |
---|
| 869 | totcount = totcount + 1; |
---|
| 870 | if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) { ////////////////////////// IF NUMERO 1////////////////// |
---|
| 871 | |
---|
| 872 | // first do rotation into collimator frame |
---|
| 873 | |
---|
| 874 | x = part.x * cos(c_rotation) + sin(c_rotation) * part.y; |
---|
| 875 | z = part.y * cos(c_rotation) - sin(c_rotation) * part.x; |
---|
| 876 | xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp; |
---|
| 877 | zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp; |
---|
| 878 | |
---|
| 879 | //cout<<"PARTICULE NON PERDU>>><<"<<endl; |
---|
| 880 | |
---|
| 881 | //To define the correct value of the mirror for the install locations |
---|
| 882 | // of the crystal in SPS & LHC, inner side or outside of the vacuum chamber???? |
---|
| 883 | // by Zhang @ CERN, 23/04/2014..... |
---|
| 884 | //mirror = 1; |
---|
[8] | 885 | |
---|
[15] | 886 | |
---|
| 887 | x = mirror * x; |
---|
| 888 | xp = mirror * xp; |
---|
[5] | 889 | |
---|
[15] | 890 | //Shift with opening and offset |
---|
[5] | 891 | |
---|
[15] | 892 | x = x - c_aperture / 2 - c_offset; |
---|
[5] | 893 | |
---|
[15] | 894 | // Include collimator tilt |
---|
[5] | 895 | |
---|
[15] | 896 | if (tiltangle > 0) { ////////////////////////// IF NUMERO 2////////////////// |
---|
| 897 | xp = xp - tiltangle; |
---|
| 898 | } //////////////////////////FIN IF NUMERO 2////////////////// |
---|
| 899 | else if (tiltangle < 0) { ////////////////////////// IF NUMERO 3////////////////// |
---|
| 900 | x = x + sin(tiltangle) * c_length; |
---|
| 901 | xp = xp - tiltangle; |
---|
| 902 | } //////////////////////////FIN IF NUMERO 3////////////////// |
---|
[5] | 903 | |
---|
[15] | 904 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 4////////////////// |
---|
| 905 | |
---|
| 906 | s_shift = s; |
---|
| 907 | shift = crys.Rcurv * (1 - cos(Cry_tilt)); |
---|
| 908 | if (Cry_tilt < (-crys.cry_bend) ) { ////////////////////////// IF NUMERO 5////////////////// |
---|
| 909 | shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) ); |
---|
| 910 | } //////////////////////////FIN IF NUMERO 5////////////////// |
---|
| 911 | x_shift = x - shift; |
---|
| 912 | } //////////////////////////fin IF NUMERO 4////////////////// |
---|
| 913 | else { //////////////////////////IF NUMERO 6////////////////// |
---|
| 914 | //cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl; |
---|
| 915 | s_shift = s; |
---|
| 916 | x_shift = x; |
---|
| 917 | } //////////////////////////fin IF NUMERO 6////////////////// |
---|
| 918 | // cout<<"debug - S shift" << s_shift<<endl; |
---|
| 919 | //cout<<"debug - X shift" << x_shift <<endl; |
---|
| 920 | // |
---|
| 921 | // 2nd transformation: rotation |
---|
| 922 | s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt); |
---|
| 923 | x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt); |
---|
| 924 | xp_rot = xp - Cry_tilt; |
---|
| 925 | //cout<< "debug - S rot" << s_rot <<endl; |
---|
| 926 | //cout<<"debug - X rot" << x_rot <<endl; |
---|
| 927 | //cout<<"debug - XP rot" << xp_rot<<endl; |
---|
| 928 | |
---|
| 929 | // 3rd transformation: drift to the new coordinate s=0 |
---|
| 930 | xp = xp_rot; |
---|
| 931 | x = x_rot - xp_rot * s_rot; |
---|
| 932 | z = z - zp * s_rot; |
---|
| 933 | s = 0; |
---|
| 934 | //cout<<"debug - S cryRF" << s_rot <<endl; |
---|
| 935 | //cout<<"debug - X cryRF" << x_rot <<endl; |
---|
| 936 | //cout<<"debug - XP cryRF" << xp_rot <<endl; |
---|
[5] | 937 | |
---|
| 938 | |
---|
[15] | 939 | |
---|
| 940 | //******************************************************************************************************************** |
---|
| 941 | //----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal--------------------------------------------------------------- |
---|
| 942 | //cout<<" X " << x<<endl; |
---|
| 943 | //Particle hit the crystal |
---|
| 944 | if (x >= 0) { ////////////////////////// IF NUMERO 7////////////////// |
---|
| 945 | s_impact = s_in0; //!(for the first impact) |
---|
| 946 | |
---|
[17] | 947 | if(0){ |
---|
[15] | 948 | cout<<"hit the cry entrance face"<<endl; |
---|
| 949 | cout<<"impact at s = "<< s_impact<< ", x = "<<x_in0<<endl; |
---|
| 950 | cout<<"with angle xp = "<<xp<<endl; |
---|
| 951 | cout<<"s before: "<< s<<endl; |
---|
| 952 | } |
---|
| 953 | |
---|
| 954 | part.x = x; |
---|
| 955 | part.xp = xp; |
---|
| 956 | part.y = z; |
---|
| 957 | part.yp = zp; |
---|
| 958 | part.PC = p; |
---|
[5] | 959 | |
---|
[15] | 960 | //------------------------------------------------------------------------------ |
---|
| 961 | // The proton crystal routine in Fortran version |
---|
| 962 | // |
---|
| 963 | //------------------------------------------------------------------------------ |
---|
| 964 | |
---|
| 965 | //Fortran crystal routine in SixTrack |
---|
| 966 | // Modified by Jianfeng Zhang @ LAL, 30/04/2014 |
---|
| 967 | double crysCry_length=crys.Cry_length; |
---|
| 968 | double crysRcurv=crys.Rcurv; |
---|
| 969 | double crysC_xmax=crys.C_xmax; |
---|
| 970 | double crysC_ymax=crys.C_ymax; |
---|
| 971 | double crysAlayer=crys.Alayer; |
---|
| 972 | int crysC_orient=crys.C_orient; |
---|
| 973 | double crysmiscut = crys.miscut; |
---|
| 974 | int IS =0; |
---|
| 975 | double length =0.0; |
---|
| 976 | IS = mat -7 ; |
---|
| 977 | length = crys.Cry_length; |
---|
| 978 | double eUmIS = crys.eUm[crys.IS]; |
---|
| 979 | double AIIS = crys.AI[crys.IS]; |
---|
| 980 | |
---|
| 981 | |
---|
| 982 | |
---|
| 983 | double crysparaDes[4]; |
---|
| 984 | crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3]; |
---|
| 985 | double crysparaDlyi[4]; |
---|
| 986 | crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3]; |
---|
| 987 | double crysparaDlri[4]; |
---|
| 988 | crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3]; |
---|
| 989 | |
---|
| 990 | double crysparaeUmIS[4]; |
---|
| 991 | crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3]; |
---|
| 992 | double crysparaAIIS[4]; |
---|
| 993 | crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3]; |
---|
| 994 | |
---|
| 995 | // cout << "before enter crystal routine...."<<endl; |
---|
[13] | 996 | |
---|
[15] | 997 | // cryst_(&IS,&x,&xp,&z,&zp,&p,&length,&crysCry_length,&crysRcurv,&crysC_xmax,&crysC_ymax,&crysAlayer,&crysC_orient,&crysmiscut, &Proc,&layerflag, &part.Chann); |
---|
| 998 | |
---|
| 999 | // cryst_(&length,&layerflag, crypara,partpara,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS,Proc); |
---|
| 1000 | |
---|
| 1001 | cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv, |
---|
| 1002 | &crys.s_length,&crys.IS,&crys.NAM, |
---|
| 1003 | &crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length, |
---|
| 1004 | &crys.ZN,&crys.C_xmax,&crys.C_ymax, |
---|
| 1005 | &part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms, |
---|
| 1006 | & part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS); |
---|
| 1007 | |
---|
[13] | 1008 | |
---|
[15] | 1009 | |
---|
| 1010 | //convert back the Fortran variables to C++ variables |
---|
| 1011 | //mat = IS +7; |
---|
| 1012 | |
---|
| 1013 | // cout << "IS in Fortran = ..." << crys.IS << endl; |
---|
[5] | 1014 | |
---|
[15] | 1015 | crys.IS = crys.IS - 1; |
---|
| 1016 | |
---|
| 1017 | // cout << "IS in C++ = ..." << crys.IS << endl; |
---|
| 1018 | |
---|
| 1019 | crys.Cry_length = length; |
---|
[16] | 1020 | |
---|
| 1021 | x=part.x; |
---|
| 1022 | xp=part.xp; |
---|
| 1023 | z=part.y; |
---|
| 1024 | zp=part.yp; |
---|
| 1025 | p=part.PC; |
---|
| 1026 | s=part.s; |
---|
[15] | 1027 | |
---|
| 1028 | |
---|
| 1029 | /* |
---|
| 1030 | crys.Rcurv = crysRcurv; |
---|
| 1031 | crys.C_xmax = crysC_xmax; |
---|
| 1032 | crys.C_ymax = crysC_ymax; |
---|
| 1033 | crys.Alayer = crysAlayer; |
---|
| 1034 | crys.C_orient = crysC_orient; |
---|
| 1035 | crys.miscut = crysmiscut; |
---|
| 1036 | */ |
---|
| 1037 | |
---|
| 1038 | // cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl; |
---|
| 1039 | // //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
| 1040 | // bases(sim); |
---|
| 1041 | // //cout << part.xp<< " XP "<<endl; |
---|
| 1042 | // bool crossing(cross(sim)); |
---|
| 1043 | // //cout<<part.xp-xp0<<endl; |
---|
| 1044 | // //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
| 1045 | // bool layering; |
---|
| 1046 | // if (crossing != 0) { ////////////////////////// IF NUMERO 8////////////////// |
---|
| 1047 | // layering = layer(sim); |
---|
| 1048 | // //cout<<"Crossing !!!!"<<endl; |
---|
| 1049 | // |
---|
| 1050 | // if (layering == true) { ////////////////////////// IF NUMERO 9////////////////// |
---|
| 1051 | // parameters(sim); |
---|
| 1052 | // //cout<<" Layering =true " <<endl; |
---|
| 1053 | // |
---|
| 1054 | // if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10////////////////// |
---|
| 1055 | // channel(sim); |
---|
| 1056 | // //cout<< " CHANNELING " << endl; |
---|
| 1057 | // } ////////////////////////// fin IF NUMERO 10////////////////// |
---|
| 1058 | // else { ////////////////////////// IF NUMERO 11////////////////// |
---|
| 1059 | // reflection(sim); |
---|
| 1060 | // //cout<<" REFLECTION "<<endl; |
---|
| 1061 | // } ////////////////////////// fin IF NUMERO 11////////////////// |
---|
| 1062 | // |
---|
| 1063 | // |
---|
| 1064 | // } ////////////////////////// fin IF NUMERO 9////////////////// |
---|
| 1065 | // } ////////////////////////// fin IF NUMERO 8////////////////// |
---|
| 1066 | |
---|
| 1067 | //-------------------------------------------------------------------------------------------------- |
---|
| 1068 | |
---|
[13] | 1069 | |
---|
[15] | 1070 | // xp = part.xp; |
---|
| 1071 | |
---|
| 1072 | //} |
---|
[5] | 1073 | |
---|
[15] | 1074 | s = crys.Rcurv * sin(crys.cry_bend); |
---|
| 1075 | zlm = crys.Rcurv * sin(crys.cry_bend); |
---|
| 1076 | //cout<<"process:"<<PROC<<endl; |
---|
| 1077 | //cout<<"s after"<< s<<endl; |
---|
| 1078 | //cout<<"part.xp"<<part.xp<<endl; |
---|
| 1079 | |
---|
| 1080 | //cout<<"----------------------1111111111111111111111111-----------------"<<endl; |
---|
| 1081 | |
---|
[17] | 1082 | if(layerflag == 0){ |
---|
| 1083 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 1084 | part.effet = 2; |
---|
| 1085 | } |
---|
| 1086 | |
---|
[15] | 1087 | if (strncmp(Proc,"out",3)==0) |
---|
| 1088 | crossing = false; |
---|
| 1089 | else |
---|
| 1090 | crossing = true; |
---|
| 1091 | |
---|
| 1092 | if (crossing == true){ |
---|
[17] | 1093 | hitflag = true; |
---|
| 1094 | nhit = nhit + 1; |
---|
| 1095 | part.nhit = nhit; |
---|
[15] | 1096 | // lhit = 100000000*ie + ITURN; |
---|
| 1097 | impact = x_in0; |
---|
| 1098 | indiv = xp_in0; |
---|
| 1099 | } |
---|
| 1100 | } |
---|
| 1101 | // Particle not hit the crystal |
---|
| 1102 | ////////////////////////// fin IF NUMERO 7////////////////// |
---|
| 1103 | else { ////////////////////////// IF NUMERO 12////////////////// |
---|
| 1104 | |
---|
| 1105 | //cout<<"OU LA ABBBBBBBBBBBB" <<endl; |
---|
| 1106 | xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv)); |
---|
| 1107 | //cout<<"RCURV "<<crys.Rcurv<<endl; |
---|
| 1108 | //cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl; |
---|
| 1109 | |
---|
| 1110 | //cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl; |
---|
| 1111 | //cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl; |
---|
| 1112 | if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13////////////////// |
---|
| 1113 | |
---|
| 1114 | //if it hits the crystal, calculate in which point and apply the |
---|
| 1115 | //transformation and drift to that point |
---|
| 1116 | a_eq = (1 + xp * xp); |
---|
| 1117 | b_eq = 2 * xp * (x - crys.Rcurv); |
---|
| 1118 | c_eq = -2 * x * crys.Rcurv + x * x; |
---|
| 1119 | Delta = b_eq * b_eq - 4 * (a_eq * c_eq); |
---|
| 1120 | s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq); |
---|
| 1121 | //cout<< " s int "<<s_int <<endl; |
---|
| 1122 | if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14////////////////// |
---|
| 1123 | // transform to a new ref system:shift and rotate |
---|
| 1124 | |
---|
| 1125 | x_int = xp * s_int + x; |
---|
| 1126 | xp_int = xp; |
---|
| 1127 | z = z + zp * s_int; |
---|
| 1128 | x = 0; |
---|
| 1129 | s = 0; |
---|
| 1130 | |
---|
[5] | 1131 | |
---|
[15] | 1132 | // tilt_int=2*X_int/S_int |
---|
| 1133 | tilt_int = s_int / crys.Rcurv; |
---|
| 1134 | xp = xp - tilt_int; |
---|
| 1135 | |
---|
| 1136 | //--------------------------------------------------------------------------------------- |
---|
| 1137 | // Fortran Cyrstal routine for the protons |
---|
| 1138 | // |
---|
| 1139 | //--------------------------------------------------------------------------------------- |
---|
| 1140 | |
---|
| 1141 | part.x = x; |
---|
| 1142 | part.xp = xp; |
---|
| 1143 | part.y = z; |
---|
| 1144 | part.yp = zp; |
---|
| 1145 | part.PC = p; |
---|
| 1146 | crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv); |
---|
| 1147 | |
---|
[13] | 1148 | |
---|
| 1149 | |
---|
[15] | 1150 | /* |
---|
[13] | 1151 | |
---|
[15] | 1152 | //Fortran crystal routine in SixTrack |
---|
| 1153 | // Modified by Jianfeng Zhang @ LAL, 30/04/2014 |
---|
| 1154 | CRYST_(mat-7,x,xp,z,zp,p,cry_length); |
---|
| 1155 | */ |
---|
| 1156 | |
---|
| 1157 | |
---|
| 1158 | //------------------------------------------------------------------------------ |
---|
| 1159 | // The proton crystal routine in Fortran version |
---|
| 1160 | // |
---|
| 1161 | //------------------------------------------------------------------------------ |
---|
| 1162 | |
---|
| 1163 | //Fortran crystal routine in SixTrack |
---|
| 1164 | // Modified by Jianfeng Zhang @ LAL, 30/04/2014 |
---|
| 1165 | double crysCry_length=crys.Cry_length; |
---|
| 1166 | double crysRcurv=crys.Rcurv; |
---|
| 1167 | double crysC_xmax=crys.C_xmax; |
---|
| 1168 | double crysC_ymax=crys.C_ymax; |
---|
| 1169 | double crysAlayer=crys.Alayer; |
---|
| 1170 | int crysC_orient=crys.C_orient; |
---|
| 1171 | double crysmiscut = crys.miscut; |
---|
| 1172 | int IS =0; |
---|
| 1173 | double length =0.0; |
---|
| 1174 | IS = mat -7 ; |
---|
| 1175 | length = crys.Cry_length; |
---|
| 1176 | double eUmIS = crys.eUm[crys.IS]; |
---|
| 1177 | double AIIS = crys.AI[crys.IS]; |
---|
| 1178 | |
---|
| 1179 | |
---|
| 1180 | double crysparaDes[4]; |
---|
| 1181 | crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3]; |
---|
| 1182 | double crysparaDlyi[4]; |
---|
| 1183 | crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3]; |
---|
| 1184 | double crysparaDlri[4]; |
---|
| 1185 | crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3]; |
---|
| 1186 | |
---|
| 1187 | double crysparaeUmIS[4]; |
---|
| 1188 | crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3]; |
---|
| 1189 | double crysparaAIIS[4]; |
---|
| 1190 | crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3]; |
---|
| 1191 | |
---|
| 1192 | |
---|
| 1193 | |
---|
| 1194 | |
---|
| 1195 | |
---|
| 1196 | // cout << "before enter crystal routine...."<<endl; |
---|
| 1197 | |
---|
| 1198 | cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv, |
---|
| 1199 | &crys.s_length,&crys.IS,&crys.NAM, |
---|
| 1200 | &crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length, |
---|
| 1201 | &crys.ZN,&crys.C_xmax,&crys.C_ymax, |
---|
| 1202 | &part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms, |
---|
| 1203 | & part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS); |
---|
| 1204 | |
---|
| 1205 | //convert back the Fortran variables to C++ variables |
---|
| 1206 | //mat = IS +7; |
---|
| 1207 | crys.Cry_length = length; |
---|
| 1208 | crys.IS = crys.IS - 1; |
---|
| 1209 | |
---|
[5] | 1210 | |
---|
[16] | 1211 | |
---|
| 1212 | x=part.x; |
---|
| 1213 | xp=part.xp; |
---|
| 1214 | z=part.y; |
---|
| 1215 | zp=part.yp; |
---|
| 1216 | p=part.PC; |
---|
| 1217 | s=part.s; |
---|
[15] | 1218 | |
---|
| 1219 | /* |
---|
| 1220 | crys.Rcurv = crysRcurv; |
---|
| 1221 | crys.C_xmax = crysC_xmax; |
---|
| 1222 | crys.C_ymax = crysC_ymax; |
---|
| 1223 | crys.Alayer = crysAlayer; |
---|
| 1224 | crys.C_orient = crysC_orient; |
---|
| 1225 | crys.miscut = crysmiscut; |
---|
| 1226 | */ |
---|
[5] | 1227 | |
---|
[15] | 1228 | /* |
---|
| 1229 | bases(sim); |
---|
| 1230 | |
---|
| 1231 | //cout<<"part.xp"<<part.xp<<endl; |
---|
[5] | 1232 | |
---|
[15] | 1233 | bool crossing(cross(sim)); |
---|
| 1234 | //cout<<"Crossing !!!!"<<endl; |
---|
| 1235 | bool layering(layer(sim)); |
---|
[5] | 1236 | |
---|
[15] | 1237 | if (layering == true) { |
---|
| 1238 | //cout<<" Layering =true " <<endl; |
---|
| 1239 | parameters(sim); |
---|
| 1240 | if ((abs(part.xp_rel) < crys.xpcrit)) { |
---|
| 1241 | channel(sim); |
---|
| 1242 | //cout<< " CHANNELING " << endl; |
---|
| 1243 | } else { |
---|
| 1244 | reflection(sim); |
---|
| 1245 | //cout<< " REFLECTION " << endl; |
---|
| 1246 | } |
---|
| 1247 | |
---|
| 1248 | } |
---|
| 1249 | */ |
---|
| 1250 | //--------------------------------------------------------------------------------------- |
---|
[5] | 1251 | |
---|
| 1252 | |
---|
[15] | 1253 | |
---|
[5] | 1254 | //cout<<"part.xp"<<part.xp<<endl; |
---|
| 1255 | //cout<<"----------------------2222222222222222222-----------------"<<endl; |
---|
[15] | 1256 | |
---|
| 1257 | s = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
| 1258 | zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
| 1259 | |
---|
[5] | 1260 | |
---|
[15] | 1261 | if(layerflag == 0){ |
---|
| 1262 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 1263 | part.effet = 2; |
---|
| 1264 | } |
---|
| 1265 | if (strncmp(Proc,"out",3)==0) |
---|
| 1266 | crossing = false; |
---|
| 1267 | else |
---|
| 1268 | crossing = true; |
---|
| 1269 | |
---|
| 1270 | if (crossing == true) { ////////////////////////// IF NUMERO 15////////////////// |
---|
| 1271 | x_rot = x_int; |
---|
| 1272 | s_rot = s_int; |
---|
| 1273 | xp_rot = xp_int; |
---|
| 1274 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
| 1275 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
| 1276 | xp_shift = xp_rot + Cry_tilt; |
---|
| 1277 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16////////////////// |
---|
| 1278 | s_impact = s_shift; |
---|
| 1279 | x_in0 = x_shift + shift; |
---|
| 1280 | xp_in0 = xp_shift; |
---|
| 1281 | } ////////////////////////// fin IF NUMERO 16////////////////// |
---|
| 1282 | else { ////////////////////////// IF NUMERO 17////////////////// |
---|
| 1283 | x_in0 = x_shift; |
---|
| 1284 | s_impact = s_shift; |
---|
| 1285 | xp_in0 = xp_shift; |
---|
| 1286 | } ////////////////////////// fin IF NUMERO 17////////////////// |
---|
[17] | 1287 | hitflag = true; |
---|
[15] | 1288 | nhit = nhit + 1; |
---|
[17] | 1289 | part.nhit = nhit; |
---|
[15] | 1290 | //lhit = 100000000*ie + ITURN; |
---|
| 1291 | impact = x_in0; |
---|
| 1292 | indiv = xp_in0; |
---|
| 1293 | } ////////////////////////// fin IF NUMERO 15////////////////// |
---|
| 1294 | |
---|
| 1295 | |
---|
[5] | 1296 | |
---|
[15] | 1297 | x_temp = x; |
---|
| 1298 | s_temp = s; |
---|
| 1299 | xp_temp = xp; |
---|
| 1300 | s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int); |
---|
| 1301 | x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int); |
---|
| 1302 | xp = xp_temp + tilt_int; |
---|
| 1303 | |
---|
| 1304 | // 2nd: shift back the 2 axis |
---|
| 1305 | x = x + x_int; |
---|
| 1306 | s = s + s_int; |
---|
| 1307 | |
---|
| 1308 | } ////////////////////////// fin IF NUMERO 14////////////////// |
---|
| 1309 | // Finish if (s_int < crys.Rcurv * sin(crys.cry_bend)) |
---|
| 1310 | else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"///////////////////////// |
---|
| 1311 | |
---|
| 1312 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
| 1313 | x = x + s * xp; |
---|
| 1314 | z = z + s * zp; |
---|
| 1315 | ++part.nout; |
---|
| 1316 | //cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl; |
---|
| 1317 | }// Finish else if (s_int < crys.Rcurv * sin(crys.cry_bend)) |
---|
| 1318 | |
---|
| 1319 | }else { |
---|
| 1320 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
| 1321 | x = x + s * xp; |
---|
| 1322 | z = z + s * zp; |
---|
| 1323 | ++part.nout; |
---|
| 1324 | //cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl; |
---|
| 1325 | } |
---|
| 1326 | //} ////////////////////////// fin IF NUMERO 13////////////////// |
---|
| 1327 | // Finish else if ( xp >= xp_tangent) |
---|
| 1328 | |
---|
| 1329 | }////////////////////// 12 |
---|
| 1330 | // Finish particle not hit the crystal else (x < 0) |
---|
| 1331 | |
---|
[13] | 1332 | |
---|
| 1333 | |
---|
[15] | 1334 | |
---|
| 1335 | //************************************************************************************************************************** |
---|
| 1336 | //trasform back from the crystal to the collimator reference system |
---|
| 1337 | // 1st: un-rotate the coordinates |
---|
| 1338 | |
---|
[5] | 1339 | |
---|
[15] | 1340 | x_rot = x; |
---|
| 1341 | s_rot = s; |
---|
| 1342 | xp_rot = xp; |
---|
| 1343 | //cout<< "debug - S Cry RF 2" , s_rot <<endl; |
---|
| 1344 | //cout<<"debug - X Cry RF 2" , x_rot <<endl; |
---|
| 1345 | //cout<<"debug - XP Cry RF 2" , xp_rot <<endl; |
---|
| 1346 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
| 1347 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
| 1348 | xp_shift = xp_rot + Cry_tilt; |
---|
| 1349 | // 2nd: shift back the reference frame |
---|
| 1350 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18////////////////// |
---|
| 1351 | s = s_shift; |
---|
| 1352 | x = x_shift + shift; |
---|
| 1353 | xp = xp_shift; |
---|
| 1354 | } ////////////////////////// fin IF NUMERO 18////////////////// |
---|
| 1355 | else { |
---|
| 1356 | x = x_shift; |
---|
| 1357 | s = s_shift; |
---|
| 1358 | xp = xp_shift; |
---|
| 1359 | } |
---|
| 1360 | // 3rd: shift to new S=Length position |
---|
| 1361 | x = xp * (c_length - s) + x; |
---|
| 1362 | z = zp * (c_length - s) + z; |
---|
| 1363 | s = c_length; |
---|
[5] | 1364 | |
---|
| 1365 | |
---|
[17] | 1366 | nabs = 0; //particle lost or not lost..... |
---|
[15] | 1367 | part.effet = nabs; |
---|
| 1368 | //cout<< "debug - S Coll RF 2" , s_rot <<endl; |
---|
| 1369 | //cout<<"debug - X Coll RF 2" , x_rot <<endl; |
---|
| 1370 | //cout<<"debug - XP Coll RF 2" , xp_rot <<endl; |
---|
[5] | 1371 | |
---|
[15] | 1372 | //cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl; |
---|
[5] | 1373 | |
---|
[15] | 1374 | |
---|
| 1375 | if(strncmp(Proc,"out",3)==0){ |
---|
| 1376 | part.nout = part.nout + 1; /////////////COUNT///////////// |
---|
| 1377 | part.effet = 1; |
---|
| 1378 | } |
---|
| 1379 | if(layerflag == 0){ |
---|
| 1380 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 1381 | part.effet = 2; |
---|
| 1382 | } |
---|
| 1383 | if(strncmp(Proc,"AM",2)==0){ |
---|
| 1384 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 1385 | part.effet = 2; |
---|
| 1386 | } |
---|
| 1387 | if(strncmp(Proc,"CH",2)==0){ |
---|
| 1388 | part.nchann = part.nchann + 1; |
---|
| 1389 | part.effet = 3; |
---|
| 1390 | } |
---|
| 1391 | if(strncmp(Proc,"DC",2)==0){ |
---|
| 1392 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
| 1393 | part.effet = 4;; |
---|
| 1394 | } |
---|
| 1395 | if(strncmp(Proc,"VR",2)==0){ |
---|
| 1396 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
| 1397 | part.effet = 5; |
---|
| 1398 | } |
---|
| 1399 | if(strncmp(Proc,"VC",2)==0){ |
---|
| 1400 | part.nvcapt = part.nvcapt + 1; /////////////COUNT///////////// |
---|
| 1401 | part.effet = 6;} |
---|
[17] | 1402 | if(strncmp(Proc,"absorbed",8==0)){ |
---|
| 1403 | part.nabsorbed = part.nabsorbed +1; |
---|
| 1404 | part.effet = 7; |
---|
| 1405 | } |
---|
| 1406 | /* |
---|
| 1407 | if (PROC(1:3).eq.'pne') |
---|
| 1408 | PROC_dan=7 |
---|
| 1409 | if (PROC(1:3).eq.'ppe') |
---|
| 1410 | PROC_dan=8 |
---|
| 1411 | if (PROC(1:4).eq.'diff') |
---|
| 1412 | PROC_dan=9 |
---|
| 1413 | if (PROC(1:4).eq.'ruth') |
---|
| 1414 | PROC_dan=10 |
---|
| 1415 | */ |
---|
| 1416 | |
---|
[15] | 1417 | |
---|
| 1418 | // ????????????????????????????????????? |
---|
| 1419 | if(part.effet == 2) { |
---|
| 1420 | part_abs = 0; |
---|
| 1421 | }else { |
---|
| 1422 | part_abs = 1; |
---|
| 1423 | } |
---|
[5] | 1424 | |
---|
| 1425 | |
---|
[15] | 1426 | //========================================================================================================================= |
---|
| 1427 | // Transform back to particle coordinates with opening and offset |
---|
[5] | 1428 | |
---|
[15] | 1429 | if (part_abs == 0) { ////////////////////////// IF NUMERO 19////////////////// |
---|
| 1430 | // |
---|
| 1431 | // Include collimator tilt |
---|
[5] | 1432 | |
---|
[15] | 1433 | if (tiltangle > 0) { |
---|
| 1434 | x = x + tiltangle * c_length; |
---|
| 1435 | xp = xp + tiltangle; |
---|
| 1436 | } else if (tiltangle < 0) { |
---|
| 1437 | x = x + tiltangle * c_length; |
---|
| 1438 | xp = xp + tiltangle; |
---|
| 1439 | |
---|
| 1440 | x = x - sin(tiltangle) * c_length; |
---|
| 1441 | } |
---|
[5] | 1442 | |
---|
[15] | 1443 | // Transform back to particle coordinates with opening and offset |
---|
| 1444 | |
---|
| 1445 | z00 = z; |
---|
| 1446 | x00 = x + mirror * c_offset; |
---|
| 1447 | x = x + c_aperture / 2 + mirror * c_offset; |
---|
[5] | 1448 | |
---|
[15] | 1449 | //+ Now mirror at the horizontal axis for negative X offset |
---|
[5] | 1450 | |
---|
[15] | 1451 | x = mirror * x; |
---|
| 1452 | xp = mirror * xp; |
---|
[5] | 1453 | |
---|
[15] | 1454 | // Last do rotation into collimator frame |
---|
| 1455 | |
---|
| 1456 | part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation); |
---|
| 1457 | part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation); |
---|
| 1458 | part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation); |
---|
| 1459 | part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation); |
---|
[5] | 1460 | |
---|
[15] | 1461 | |
---|
| 1462 | |
---|
| 1463 | //********************************************************************************** |
---|
| 1464 | //For output.... by Zhang @ CERN, 24/04/2014 |
---|
| 1465 | //p_in = (1 + dpop) * p0; |
---|
| 1466 | //From Dianiele |
---|
| 1467 | //p_in = p; |
---|
| 1468 | //s_in = s_in + s; |
---|
| 1469 | part.PC = p; |
---|
| 1470 | part.s = part.s + s; |
---|
| 1471 | |
---|
| 1472 | if (part.effet == 1) { ////////////////////////// IF NUMERO 21////////////////// |
---|
[5] | 1473 | |
---|
[15] | 1474 | part.x = 99.99e-3; |
---|
| 1475 | part.y = 99.99e-3; |
---|
| 1476 | cout<< "Mean that the particules is lost"<<endl; |
---|
| 1477 | } ////////////////////////// fin IF NUMERO 21////////////////// |
---|
| 1478 | } ////////////////////////// fin IF NUMERO 19////////////////// |
---|
[5] | 1479 | |
---|
| 1480 | |
---|
[15] | 1481 | // End of check for particles not being lost before (see @330) |
---|
[5] | 1482 | |
---|
[15] | 1483 | //} ////////////////////////// fin IF NUMERO 12////////////////// |
---|
| 1484 | // End of loop over all particles |
---|
[5] | 1485 | |
---|
[15] | 1486 | |
---|
| 1487 | |
---|
| 1488 | } //collimator with length = 0 ////////////////////////// fin IF NUMERO 1////////////////// |
---|
[13] | 1489 | |
---|
[15] | 1490 | // =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM==================================================================== |
---|
[16] | 1491 | |
---|
[15] | 1492 | part.xp=xp; |
---|
| 1493 | part.yp=zp; |
---|
| 1494 | part.y=z; |
---|
| 1495 | part.x=x; |
---|
| 1496 | part.PC=p; |
---|
[16] | 1497 | |
---|
[15] | 1498 | xp1 = part.xp; |
---|
[5] | 1499 | |
---|
[15] | 1500 | |
---|
| 1501 | //if(xp1-xp0!=0){ |
---|
| 1502 | //sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl; |
---|
| 1503 | //} |
---|
| 1504 | //cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl; |
---|
[5] | 1505 | |
---|
[15] | 1506 | |
---|
| 1507 | AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction |
---|
| 1508 | AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction |
---|
| 1509 | RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction |
---|
| 1510 | RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction |
---|
| 1511 | av_pc1 = av_pc1 + part.PC; |
---|
[5] | 1512 | |
---|
[17] | 1513 | if(0){ |
---|
[15] | 1514 | cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl; |
---|
| 1515 | //cout<<endl; |
---|
[17] | 1516 | } |
---|
| 1517 | //only record the particles that hit the collimator. |
---|
| 1518 | if(hitflag != true){ |
---|
[15] | 1519 | //sortieIco.setf(ios::scientific); |
---|
[17] | 1520 | sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl; |
---|
[15] | 1521 | //}//end of for.... |
---|
[17] | 1522 | } |
---|
| 1523 | |
---|
[15] | 1524 | } |
---|
| 1525 | } //End of while |
---|
[5] | 1526 | |
---|
| 1527 | |
---|
| 1528 | } else { |
---|
[15] | 1529 | cerr<< "We can not open the file !" << endl; |
---|
[5] | 1530 | } |
---|
| 1531 | entree.close(); |
---|
| 1532 | |
---|
| 1533 | |
---|
| 1534 | |
---|
| 1535 | AV_xp0 = AV_xp0 / count; |
---|
| 1536 | AV_yp0 = AV_yp0 / count; |
---|
[15] | 1537 | |
---|
[5] | 1538 | RMSxp0 = sqrt(RMSxp0 / (count)); |
---|
| 1539 | RMSyp0 = sqrt(RMSyp0 / (count)); |
---|
[15] | 1540 | |
---|
[5] | 1541 | SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0)); |
---|
| 1542 | SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0)); |
---|
| 1543 | |
---|
| 1544 | AV_xp1 = AV_xp1 / count; |
---|
| 1545 | AV_yp1 = AV_yp1 / count; |
---|
| 1546 | RMSxp1 = sqrt(RMSxp1 / (count)); |
---|
| 1547 | RMSyp1 = sqrt(RMSyp1 / (count)); |
---|
[15] | 1548 | |
---|
[5] | 1549 | SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1)); |
---|
| 1550 | SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1)); |
---|
| 1551 | av_pc0 = av_pc0 / count; |
---|
| 1552 | av_pc1 = av_pc1 / count; |
---|
| 1553 | //-------------------------------------------------------------------------- |
---|
| 1554 | //-------------write a summary file----------------------------------------- |
---|
| 1555 | //cout<<endl; |
---|
| 1556 | //cout<<endl; |
---|
| 1557 | //cout<<endl; |
---|
[15] | 1558 | |
---|
[5] | 1559 | //cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl; |
---|
| 1560 | //cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl; |
---|
| 1561 | //cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl; |
---|
| 1562 | //cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl; |
---|
[15] | 1563 | |
---|
[5] | 1564 | |
---|
| 1565 | |
---|
| 1566 | sortieIco.close(); |
---|
[15] | 1567 | // cout<<endl; |
---|
| 1568 | // affiche(); |
---|
[5] | 1569 | file_out(pas, outputpath); |
---|
[15] | 1570 | |
---|
[17] | 1571 | if(1){ |
---|
| 1572 | cout << "number of particles before the crystal: " << count<< endl; |
---|
| 1573 | cout << "number of particle hit the crystal: " << part.nhit << endl; |
---|
| 1574 | cout << "number of particles after the crystal: " << count - part.nhit << endl; |
---|
| 1575 | } |
---|
[5] | 1576 | //cout<<endl; |
---|
| 1577 | //cout<<endl; |
---|
| 1578 | //cout<<endl; |
---|
| 1579 | } |
---|
| 1580 | |
---|
| 1581 | |
---|
| 1582 | |
---|
| 1583 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1584 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1585 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1586 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1587 | |
---|
| 1588 | //WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW |
---|
| 1589 | |
---|
| 1590 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1591 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1592 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1593 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1594 | |
---|
| 1595 | |
---|
| 1596 | void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC) |
---|
| 1597 | { |
---|
| 1598 | //cout<<"The move_am function is use !!!!! "<<endl; |
---|
| 1599 | double DYA(0); |
---|
| 1600 | double Wp(0); |
---|
| 1601 | if (NAM == 0) { |
---|
| 1602 | return; |
---|
| 1603 | } |
---|
| 1604 | xp = xp * 1000; |
---|
| 1605 | yp = yp * 1000; |
---|
| 1606 | //cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl; |
---|
| 1607 | //DEI ---> dE/dx (stoping energy) |
---|
| 1608 | PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV] |
---|
| 1609 | |
---|
| 1610 | // Coulomb Scattering |
---|
| 1611 | //DYA ---> rms of coloumb scattering |
---|
| 1612 | DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad) |
---|
| 1613 | xp = xp + DYA * random_gauss_cut(); |
---|
| 1614 | yp = yp + DYA * random_gauss_cut(); |
---|
| 1615 | |
---|
| 1616 | //Elastic scattering |
---|
| 1617 | //cout<<"Plusieurs possibilitees : "<<endl; |
---|
| 1618 | if (random() <= (DZ / crys.DLAI[IS])) { |
---|
| 1619 | //cout<<"Poss. 1 !!!! "<<endl; |
---|
| 1620 | xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane |
---|
| 1621 | yp = yp + 196 * random_gauss_cut() / PC; |
---|
| 1622 | } |
---|
| 1623 | |
---|
| 1624 | //Diffraction interaction |
---|
| 1625 | if (random() <= (DZ / (DLY * 6.143))) { |
---|
| 1626 | //cout<<"Poss. 2 !!!! "<<endl; |
---|
| 1627 | xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr] |
---|
| 1628 | yp = yp + 257 * random_gauss_cut() / PC; |
---|
| 1629 | Wp = random(); |
---|
| 1630 | PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle |
---|
| 1631 | |
---|
| 1632 | } |
---|
| 1633 | |
---|
| 1634 | //Inelastic interaction |
---|
| 1635 | if (random() <= DZ / DLY) { |
---|
| 1636 | //cout<<" Absorbed !! "<< endl; |
---|
| 1637 | } |
---|
| 1638 | //cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl; |
---|
| 1639 | xp = xp / 1000; |
---|
| 1640 | yp = yp / 1000; |
---|
| 1641 | } |
---|
| 1642 | |
---|
| 1643 | |
---|
| 1644 | double SimCrys::random() |
---|
| 1645 | { |
---|
| 1646 | double scale = RAND_MAX; |
---|
| 1647 | double base = rand() / scale; |
---|
| 1648 | double fine = rand() / scale; |
---|
| 1649 | return base + fine / scale; |
---|
| 1650 | } |
---|
| 1651 | |
---|
| 1652 | double SimCrys::random_dist() |
---|
| 1653 | { |
---|
| 1654 | double scale = RAND_MAX; |
---|
| 1655 | double base = rand() / scale; |
---|
| 1656 | double fine = rand() / scale; |
---|
| 1657 | return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale)); |
---|
| 1658 | } |
---|
| 1659 | |
---|
| 1660 | double SimCrys::random_gauss() |
---|
| 1661 | { |
---|
| 1662 | const double PI (4.0 * atan(1.0)); |
---|
| 1663 | double U(random()); |
---|
| 1664 | double V(random()); |
---|
| 1665 | return (sqrt(-2 * log(U)) * cos(2 * PI * V)); |
---|
| 1666 | } |
---|
| 1667 | |
---|
| 1668 | double SimCrys::random_gauss_cut() |
---|
| 1669 | { |
---|
| 1670 | double resu; |
---|
| 1671 | do { |
---|
| 1672 | resu = random_gauss(); |
---|
| 1673 | } while (abs(resu) >= 1); |
---|
| 1674 | return resu; |
---|
| 1675 | } |
---|
| 1676 | |
---|
| 1677 | double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
| 1678 | { |
---|
| 1679 | const double PI (4.0 * atan(1.0)); |
---|
| 1680 | double U(random()); |
---|
| 1681 | double V(random()); |
---|
| 1682 | //moy=random(); |
---|
| 1683 | return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 1684 | } |
---|
| 1685 | |
---|
| 1686 | /////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1687 | //POUR 409 : |
---|
| 1688 | |
---|
| 1689 | double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
| 1690 | { |
---|
| 1691 | const double PI (4.0 * atan(1.0)); |
---|
| 1692 | double U(random()); |
---|
| 1693 | double V(random()); |
---|
| 1694 | //moy=random(); |
---|
| 1695 | return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 1696 | } |
---|
| 1697 | |
---|
| 1698 | ///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1699 | |
---|
| 1700 | void SimCrys::affiche() |
---|
| 1701 | { |
---|
| 1702 | cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
| 1703 | cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
| 1704 | cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1705 | cout << "Critical Radius :" << crys.Rcrit << endl; |
---|
| 1706 | cout << "Ratio :" << crys.ratio << endl; |
---|
| 1707 | cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
| 1708 | cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
| 1709 | cout << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
| 1710 | cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1711 | cout << endl; |
---|
| 1712 | cout << endl; |
---|
| 1713 | cout << "N.AMORPH :" << part.namor << endl; |
---|
| 1714 | cout << "N.DECHANNELING:" << part.ndechann << endl; |
---|
| 1715 | cout << "N.CHANNELING :" << part.nchann << endl; |
---|
| 1716 | cout << "N.OUT :" << part.nout << endl; |
---|
| 1717 | cout << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
[17] | 1718 | cout << "N.VCAPTURE:" << part.nvcapt << endl; |
---|
| 1719 | cout << "N.ABSORBED: "<< part.nabsorbed << endl; |
---|
[5] | 1720 | |
---|
| 1721 | //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!"; |
---|
| 1722 | } |
---|
| 1723 | |
---|
| 1724 | void SimCrys::file_out(const int& pas, string outputpath) |
---|
| 1725 | { |
---|
| 1726 | |
---|
| 1727 | ofstream out; |
---|
| 1728 | string file; |
---|
| 1729 | file = outputpath + "/crystal_output.out"; |
---|
| 1730 | |
---|
| 1731 | out.open(file.c_str(), ios::out | ios::app); |
---|
| 1732 | |
---|
| 1733 | if (out.fail()) { |
---|
| 1734 | cerr << "Warning: problem writing in the file " << file << "!" << endl; |
---|
| 1735 | } else { |
---|
| 1736 | |
---|
| 1737 | out << endl; |
---|
| 1738 | out << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
| 1739 | out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
[17] | 1740 | out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << 1e6*crys.Cry_length / crys.Rcurv << endl; |
---|
[5] | 1741 | out << "Critical Radius :" << crys.Rcrit << endl; |
---|
| 1742 | out << "Ratio :" << crys.ratio << endl; |
---|
| 1743 | out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
| 1744 | out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
| 1745 | out << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
| 1746 | out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1747 | out << endl; |
---|
| 1748 | out << endl; |
---|
| 1749 | out << "N.AMORPH :" << part.namor << endl; |
---|
| 1750 | out << "N.DECHANNELING:" << part.ndechann << endl; |
---|
| 1751 | out << "N.CHANNELING :" << part.nchann << endl; |
---|
| 1752 | out << "N.OUT :" << part.nout << endl; |
---|
| 1753 | out << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
[17] | 1754 | out << "N.VCAPTURE :" << part.nvcapt << endl ; |
---|
| 1755 | out << "N.ABSORBED: "<< part.nabsorbed << endl; |
---|
| 1756 | out << "N. of particles that hit the crystal collimator: " << part.nhit << endl; |
---|
| 1757 | cout << endl; |
---|
| 1758 | cout << endl; |
---|
| 1759 | |
---|
[5] | 1760 | out.close(); |
---|
| 1761 | } |
---|
| 1762 | |
---|
| 1763 | } |
---|
| 1764 | |
---|