[5] | 1 | #include <iostream> |
---|
| 2 | #include <fstream> |
---|
| 3 | #include <cmath> |
---|
| 4 | #include <iomanip> |
---|
| 5 | #include <string> |
---|
| 6 | #include "simcrys.h" |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | using namespace std; |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | |
---|
| 13 | SimCrys::SimCrys(Crystal crys, Partcrys part) |
---|
| 14 | : crys(crys), part(part) //,moy(moy) //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 15 | { |
---|
| 16 | } |
---|
| 17 | |
---|
| 18 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 19 | ////////////////////////// FUNCTION BASES THAT MAKE SOME BASIC CALCULATION |
---|
| 20 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 21 | |
---|
| 22 | void SimCrys::bases(SimCrys& sim) |
---|
| 23 | { |
---|
| 24 | if ( (crys.miscut < 0 ) && ( part.x > 0) && (part.x < -crys.Cry_length * tan(crys.miscut))) { |
---|
| 25 | crys.L_chan = -part.x * tan(crys.miscut); |
---|
| 26 | crys.tchan = crys.L_chan / crys.Rcurv; |
---|
| 27 | } |
---|
| 28 | part.xp_rel = part.xp - crys.miscut; |
---|
| 29 | } |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 34 | ////////////////////////FUNCTIN CROSS TO KNOW IF THE PARTICULE CROSS THE CRYSTAL |
---|
| 35 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 36 | |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | bool SimCrys::cross(SimCrys& sim) |
---|
| 41 | { |
---|
| 42 | bool out(true); |
---|
| 43 | if ((part.y < crys.ymin) || (part.y > crys.ymax) || (part.x > crys.xmax)) { |
---|
| 44 | //cout << "OUT of the Crystal !!" << part.x << " , "<<part.xp << endl; |
---|
| 45 | out = false; |
---|
| 46 | part.x = part.x + part.xp * crys.s_length; |
---|
| 47 | |
---|
| 48 | part.y = part.y + part.yp * crys.s_length; |
---|
| 49 | |
---|
| 50 | part.nout = part.nout + 1; /////////////COUNT///////////// |
---|
| 51 | part.effet = 1; |
---|
| 52 | //cout<<" BLAAAAAAAAa"<<endl; |
---|
| 53 | |
---|
| 54 | } |
---|
| 55 | |
---|
| 56 | return out; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 60 | ////////////////////////////FUNCTION LAYER THAT DETERMIN IF THE PARTICULE IS IN THE AMORPHOUS LAYER///////////////////////////// |
---|
| 61 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 62 | |
---|
| 63 | bool SimCrys::layer(SimCrys& sim) |
---|
| 64 | { |
---|
| 65 | //cout<<"La fonction layer est utilisee>>>>>>>>>>>>>>"<<endl; |
---|
| 66 | bool out(true); |
---|
| 67 | if ((part.x < crys.Alayer) || ((part.y - crys.ymin) < crys.Alayer) || ((crys.ymax - part.y) < crys.Alayer)) { |
---|
| 68 | //cout<<"Amorphous layer! " << part.x <<" , " <<part.xp<<endl; |
---|
| 69 | double a_eq; |
---|
| 70 | double b_eq; |
---|
| 71 | double c_eq; |
---|
| 72 | double delta; //a, b, c, Delta of the second order eq. |
---|
| 73 | double x0; |
---|
| 74 | double y0; |
---|
| 75 | double length_xs; //!Amorphous length |
---|
| 76 | double length_ys; //!Amorphous length |
---|
| 77 | double am_length; //!Amorphous length |
---|
| 78 | out = false; |
---|
| 79 | x0 = part.x; |
---|
| 80 | y0 = part.y; |
---|
| 81 | a_eq = (1 + part.xp * part.xp); |
---|
| 82 | b_eq = (2 * (part.x) * (part.xp) - 2 * (part.xp) * crys.Rcurv); |
---|
| 83 | c_eq = (part.x) * part.x - 2 * abs(part.x) * crys.Rcurv; |
---|
| 84 | delta = b_eq * b_eq - 4 * a_eq * c_eq; |
---|
| 85 | //cout<<"a : " << a_eq<<endl; //POUR TESTER FONCTION |
---|
| 86 | //cout<<"b: " << b_eq <<endl; //POUR TESTER FONCTION |
---|
| 87 | //cout<<"c : " << c_eq <<endl; //POUR TESTER FONCTION |
---|
| 88 | //cout<<"Delta : " << delta <<endl; //POUR TESTER FONCTION |
---|
| 89 | part.s = (-b_eq + sqrt(delta)) / ((2 * a_eq)) ; //in [m] |
---|
| 90 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
| 91 | if (part.s >= crys.s_length) { |
---|
| 92 | part.s = crys.s_length; |
---|
| 93 | } |
---|
| 94 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
| 95 | part.x = part.xp * part.s + x0 ; //in [m] |
---|
| 96 | //cout<<"x0 : " << x0 <<endl; //POUR TESTER FONCTION |
---|
| 97 | //cout<<"x : " << part.x <<endl; //POUR TESTER FONCTION |
---|
| 98 | |
---|
| 99 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
| 100 | //cout<<"Length_XS : " << length_xs <<endl; //POUR TESTER FONCTION |
---|
| 101 | |
---|
| 102 | if ((((part.yp >= 0) && ((part.y + part.yp * part.s) <= crys.ymax))) || (((part.yp < 0) && ((part.y + part.yp * part.s) >= crys.ymin)))) { |
---|
| 103 | length_ys = part.yp * length_xs; |
---|
| 104 | } else { |
---|
| 105 | part.s = (crys.ymax - part.y) / part.yp; |
---|
| 106 | length_ys = sqrt((crys.ymax - part.y) * (crys.ymax - part.y) + part.s * part.s); |
---|
| 107 | part.x = x0 + part.xp * part.s; |
---|
| 108 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
| 109 | } |
---|
| 110 | |
---|
| 111 | am_length = sqrt(length_xs * length_xs + length_ys * length_ys); |
---|
| 112 | //cout << "AM_Length :" << am_length << endl; // TESTER FONCTION |
---|
| 113 | part.s = part.s / 2; |
---|
| 114 | part.x = x0 + part.xp * part.s; |
---|
| 115 | part.y = y0 + part.yp * part.s; |
---|
| 116 | //cout<<" am_length "<<am_length<<endl; |
---|
| 117 | 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 |
---|
| 118 | //cout<<"AM LAYER !!!! "<<part.PC<<endl; |
---|
| 119 | //cout<< "Amorphous"<< endl; |
---|
| 120 | part.x = part.x + part.xp * (crys.Cry_length - part.s); |
---|
| 121 | part.y = part.y + part.yp * (crys.Cry_length - part.s); |
---|
| 122 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 123 | part.effet = 2; |
---|
| 124 | } else if ((part.x > (crys.xmax - crys.Alayer)) && (part.x < crys.xmax)) { |
---|
| 125 | out = false; |
---|
| 126 | |
---|
| 127 | 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 |
---|
| 128 | |
---|
| 129 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 130 | part.effet = 2; |
---|
| 131 | //cout<<"Fix HERE"<<endl; |
---|
| 132 | } |
---|
| 133 | return out; |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 138 | //SI LES DEUX FONCTIONS D AVT RETURNE TRUE, ALORS ON CALCULE LES PARA CRITIQUE ET D AUTRE PARA (POINT C - E) |
---|
| 139 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 140 | |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | void SimCrys::parameters(SimCrys& sim) |
---|
| 144 | { |
---|
| 145 | //THE FIRST PART IS FOR THE (110) ORIENTATION |
---|
| 146 | //cout << "The parameters function is use !!!!! "<< endl; |
---|
| 147 | double c_v1; |
---|
| 148 | double c_v2; |
---|
| 149 | |
---|
| 150 | |
---|
| 151 | crys.xpcrit0 = sqrt(2.10e-9 * crys.eUm[crys.IS] / part.PC); |
---|
| 152 | crys.Rcrit = part.PC / (2.e-6 * crys.eUm[crys.IS]) * crys.AI[crys.IS]; //if R>Rcritical=>no channeling is possible |
---|
| 153 | crys.ratio = crys.Rcurv / crys.Rcrit; |
---|
| 154 | crys.xpcrit = crys.xpcrit0 * (crys.Rcurv - crys.Rcrit) / crys.Rcurv; |
---|
| 155 | c_v1 = 1.7; |
---|
| 156 | c_v2 = -1.5; |
---|
| 157 | if (crys.ratio <= 1) { |
---|
| 158 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.4 * crys.ratio); //rms scattering |
---|
| 159 | part.Ang_avr = c_v2 * crys.xpcrit0 * 0.05 * crys.ratio; //average angle reflection |
---|
| 160 | part.Vcapt = 0.0; |
---|
| 161 | } else if (crys.ratio <= 3) { |
---|
| 162 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.571 * 0.3 * crys.ratio + 0.85); //rms scattering |
---|
| 163 | part.Ang_avr = c_v2 * crys.xpcrit0 * (0.1972 * crys.ratio - 0.1472); //average angle reflection |
---|
| 164 | 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) |
---|
| 165 | } else { |
---|
| 166 | part.Ang_rms = c_v1 * crys.xpcrit0 * (1 / crys.ratio); |
---|
| 167 | part.Ang_avr = c_v2 * crys.xpcrit0 * (1 - 1.6667 / crys.ratio); |
---|
| 168 | part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2); |
---|
| 169 | } |
---|
| 170 | //cout <<crys.xpcrit<<endl; |
---|
| 171 | if (crys.C_orient == 2) { |
---|
| 172 | part.Ang_avr = part.Ang_avr * 0.93; |
---|
| 173 | part.Ang_rms = part.Ang_rms * 1.05; // FOR (111) ORIENTATION |
---|
| 174 | crys.xpcrit = crys.xpcrit * 0.98; |
---|
| 175 | } |
---|
| 176 | //cout <<crys.xpcrit<<endl; |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 183 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE I CF SHEMA!!! C EST LA PARTIE CHANNELING D OU LE NOM |
---|
| 184 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 185 | |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | |
---|
| 189 | void SimCrys::channel(SimCrys& sim) |
---|
| 190 | { |
---|
| 191 | //cout << "La fonction channel est utilisee!!!! "<<endl; |
---|
| 192 | |
---|
| 193 | |
---|
| 194 | part.Chann = sqrt(crys.xpcrit * crys.xpcrit - part.xp_rel * part.xp_rel) / crys.xpcrit0; //probability of channeling/volume capture |
---|
| 195 | |
---|
| 196 | //cout <<"Channeling Prob :" << part.Chann <<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
| 197 | //cout <<"Channeling ANGLE :" << crys.tchan<<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | double n_atom(0.1); //probability for entering channeling near atomic planes |
---|
| 201 | double Dxp(0); //change angle from channeling [mrad] |
---|
| 202 | double TLdech1(0); //tipical dechanneling length(1) [m] |
---|
| 203 | double tdech(0); |
---|
| 204 | double Ldech(0); //Dechanneling. length |
---|
| 205 | double Sdech(0); |
---|
| 206 | //cout<<"Chann : "<<part.Chann<<endl; |
---|
| 207 | if (random() <= part.Chann) { |
---|
| 208 | TLdech1 = 0.0005 * part.PC * pow(1 - 1 / crys.ratio, 2); // calculate tipical dechanneling length(m) |
---|
| 209 | if (random() <= n_atom) { |
---|
| 210 | TLdech1 = 0.000002 * part.PC * pow(1. - 1. / crys.ratio, 2); // dechanneling length (m) for short crystal for high amplitude particles |
---|
| 211 | } |
---|
| 212 | //cout<<"@ energy " << part.PC <<" dechanneling length "<< TLdech1 << endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
| 213 | part.Dechan = log(random()); |
---|
| 214 | //cout<<"dechannn " << part.Dechan <<endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
| 215 | Ldech = -TLdech1 * part.Dechan; //careful: the dechanneling lentgh is along the trajectory of the particle -not along the longitudinal coordinate... |
---|
| 216 | tdech = Ldech / crys.Rcurv; |
---|
| 217 | Sdech = Ldech * cos(crys.miscut + 0.5 * tdech); |
---|
| 218 | // 2 option if the particule channeling !!! |
---|
| 219 | // (1) dechanneling |
---|
| 220 | // (2) full channeling |
---|
| 221 | |
---|
| 222 | // (1) : |
---|
| 223 | |
---|
| 224 | if (Ldech < crys.L_chan) { |
---|
| 225 | //cout <<"Dechanneling ! "<< endl; |
---|
| 226 | Dxp = Ldech / crys.Rcurv; //change angle from channeling [mrad] |
---|
| 227 | |
---|
| 228 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
| 229 | part.effet = 4; |
---|
| 230 | part.x = part.x + Ldech * (sin(0.5 * Dxp + crys.miscut)); // trajectory at channeling exit |
---|
| 231 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
| 232 | part.y = part.y + part.yp * Sdech; |
---|
| 233 | |
---|
| 234 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
| 235 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
| 236 | |
---|
| 237 | |
---|
| 238 | 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 |
---|
| 239 | |
---|
| 240 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * part.y; //energy loss to ionization [GeV] |
---|
| 241 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
| 242 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | // (2) : |
---|
| 246 | |
---|
| 247 | else { |
---|
| 248 | //cout <<"Channeling ! "<< endl; |
---|
| 249 | Dxp = crys.L_chan / crys.Rcurv; //change angle [mrad] |
---|
| 250 | part.x = part.x + crys.L_chan * (sin(0.5 * Dxp + crys.miscut)); //trajectory at channeling exit |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | // Removed dependence on incoming angle according to discussion with I.Yazynin on 24/10/2012 |
---|
| 254 | //part.xp = part.xp + Dxp + (random_gauss()-0.5)*crys.xpcrit*0.5; // [mrad] |
---|
| 255 | part.xp = /*part.xp*/ + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; // [mrad] |
---|
| 256 | |
---|
| 257 | part.y = part.y + crys.s_length * part.yp; |
---|
| 258 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * crys.Cry_length; // energy loss to ionization [GeV] |
---|
| 259 | part.nchann = part.nchann + 1; /////////////COUNT///////////// |
---|
| 260 | part.effet = 3; |
---|
| 261 | } |
---|
| 262 | } else { |
---|
| 263 | //cout <<"Volume Reflection ! "<< endl; |
---|
| 264 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
| 265 | part.effet = 5; |
---|
| 266 | Dxp = 0.5 * (part.xp_rel / crys.xpcrit + 1) * part.Ang_avr; |
---|
| 267 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
| 268 | /* |
---|
| 269 | next line new from sasha |
---|
| 270 | part.xp=part.xp+0.45*(part.xp/crys.xpcrit+1)*part.Ang_avr; |
---|
| 271 | valentina fix here |
---|
| 272 | */ |
---|
| 273 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 274 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 275 | |
---|
| 276 | 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 |
---|
| 277 | |
---|
| 278 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 279 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 280 | |
---|
| 281 | } |
---|
| 282 | } |
---|
| 283 | |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 287 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE II CF SHEMA!!! C EST LA PARTIE VOLUME REFLEXION D OU LE NOM |
---|
| 288 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 289 | |
---|
| 290 | |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | void SimCrys::reflection(SimCrys& sim) |
---|
| 294 | { |
---|
| 295 | //cout<<"The function reflection is used" <<endl; |
---|
| 296 | double Dxp(0); //change angle from channeling [mrad] |
---|
| 297 | double Lrefl(0); //distance of the reflection point [m] |
---|
| 298 | double Srefl(0); //distance of the reflection point [m] |
---|
| 299 | double TLdech2(0); //tipical dechanneling length(2) [m] |
---|
| 300 | double tdech(0); |
---|
| 301 | double Ldech(0); //Dechanneling. length |
---|
| 302 | double Sdech(0); |
---|
| 303 | double Red_S(0); //Reduced length (in case of dechanneling) |
---|
| 304 | double Rlength(0); //Reduced length |
---|
| 305 | |
---|
| 306 | Lrefl = (part.xp_rel) * crys.Rcurv; |
---|
| 307 | Srefl = sin(part.xp) * Lrefl; |
---|
| 308 | |
---|
| 309 | |
---|
| 310 | if ((Lrefl > 0) && (Lrefl < crys.Cry_length)) { // IE : If the VR point is inside the crystal!!!! |
---|
| 311 | |
---|
| 312 | //cout<<"VRCapt : "<<part.Vcapt<<endl; |
---|
| 313 | if ((random() > part.Vcapt) || (crys.ZN == 0)) { |
---|
| 314 | //cout<< "Volume Reflexion! "<<endl; |
---|
| 315 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
| 316 | part.effet = 5; |
---|
| 317 | part.x = part.x + part.xp * Srefl; |
---|
| 318 | part.y = part.y + part.yp * Srefl; |
---|
| 319 | Dxp = part.Ang_avr; |
---|
| 320 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
| 321 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
| 322 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
| 323 | |
---|
| 324 | if (crys.ZN > 0) { |
---|
| 325 | 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 |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
| 329 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
| 330 | } else { |
---|
| 331 | //cout<< "Volume Capture! "; |
---|
| 332 | part.nvcapt = part.nvcapt + 1; /////////////COUNT///////////// |
---|
| 333 | part.effet = 6; |
---|
| 334 | part.x = part.x + part.xp * Srefl; |
---|
| 335 | part.y = part.y + part.yp * Srefl; |
---|
| 336 | /* |
---|
| 337 | TLdech2= 0.00011*pow(part.PC,0.25)*pow((1.-1./crys.ratio),2) //dechanneling length [m] |
---|
| 338 | part.Dechan = log(1.-random()) |
---|
| 339 | Ldech = -TLdech2*part.Dechan |
---|
| 340 | next 2 lines new from sasha - different dechanneling probability |
---|
| 341 | */ |
---|
| 342 | TLdech2 = 0.01 * part.PC * pow((1. - 1. / crys.ratio), 2); // typical dechanneling length [m] |
---|
| 343 | Ldech = 0.005 * TLdech2 * pow((sqrt(0.01 - log(random())) - 0.1), 2); // (dechanneling length) [m] |
---|
| 344 | tdech = Ldech / crys.Rcurv; |
---|
| 345 | Sdech = Ldech * cos(part.xp + 0.5 * tdech); |
---|
| 346 | |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | if (Ldech < (crys.Cry_length - Lrefl)) { // for dechanneling part |
---|
| 350 | //cout<< "Dechanneling! "; |
---|
| 351 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
| 352 | part.effet = 4; |
---|
| 353 | Dxp = Ldech / crys.Rcurv; |
---|
| 354 | part.x = part.x + Ldech * (sin(0.5 * Dxp + part.xp)); //trajectory at channeling exit |
---|
| 355 | part.y = part.y + Sdech * part.yp; |
---|
| 356 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
| 357 | Red_S = crys.s_length - Srefl - Sdech; |
---|
| 358 | // move in amorphous substance---------------------------- |
---|
| 359 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
| 360 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
| 361 | |
---|
| 362 | if (crys.ZN > 0) { |
---|
| 363 | //cout<<"Amorphous! "<<endl; /////////////COUNT///////////// |
---|
| 364 | 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 |
---|
| 365 | } |
---|
| 366 | |
---|
| 367 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
| 368 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
| 369 | //cout<<endl; |
---|
| 370 | } else { |
---|
| 371 | //cout<< "Volume Capture! "<<endl; |
---|
| 372 | //Pas besoin d additionner 1 dans le compte de vol capt car deja fait!!! |
---|
| 373 | Rlength = crys.Cry_length - Lrefl; |
---|
| 374 | crys.tchan = Rlength / crys.Rcurv; |
---|
| 375 | Red_S = Rlength * cos(part.xp + 0.5 * crys.tchan); |
---|
| 376 | Dxp = (crys.Cry_length - Lrefl) / crys.Rcurv; |
---|
| 377 | part.x = part.x + sin(0.5 * Dxp + part.xp) * Rlength; // trajectory at channeling exit |
---|
| 378 | part.y = part.y + Red_S * part.yp; |
---|
| 379 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; //[mrad] |
---|
| 380 | } |
---|
| 381 | } |
---|
| 382 | } |
---|
| 383 | // move in amorphous substance for big input angles--------------- |
---|
| 384 | else { |
---|
| 385 | //cout<<"Amorphous! "<<endl; |
---|
| 386 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
| 387 | part.effet = 2; |
---|
| 388 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 389 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 390 | |
---|
| 391 | if (crys.ZN > 0) { |
---|
| 392 | 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 |
---|
| 393 | } |
---|
| 394 | |
---|
| 395 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
| 396 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | |
---|
| 402 | |
---|
| 403 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 404 | //THE FUNCTION THAT DO ALL THE SCHEMA USING THE PREVIOUS FUNCTION !!!!!!!!!!111 |
---|
| 405 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 406 | |
---|
| 407 | |
---|
| 408 | |
---|
| 409 | |
---|
| 410 | |
---|
| 411 | void SimCrys::general(SimCrys& sim, const int& pas, string outputpath, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt) |
---|
| 412 | { |
---|
| 413 | srand ( time(NULL) ); |
---|
| 414 | double AV_xp0(0); |
---|
| 415 | double AV_yp0(0); |
---|
| 416 | double RMSxp0(0); |
---|
| 417 | double RMSyp0(0); |
---|
| 418 | double AV_xp1(0); |
---|
| 419 | double AV_yp1(0); |
---|
| 420 | double long RMSxp1(0); |
---|
| 421 | double long RMSyp1(0); |
---|
| 422 | double av_pc1(0); |
---|
| 423 | double av_pc0(0); |
---|
| 424 | double SIGxp0(0); |
---|
| 425 | double SIGyp0(0); |
---|
| 426 | double SIGxp1(0); |
---|
| 427 | double SIGyp1(0); |
---|
| 428 | double xp0(0); |
---|
| 429 | double xp1(0); |
---|
| 430 | double count(0); |
---|
| 431 | double zpos(-0.1); |
---|
| 432 | double xfin(0); |
---|
| 433 | double yfin(0); |
---|
| 434 | int const col(15); |
---|
| 435 | |
---|
| 436 | |
---|
| 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/rotation |
---|
| 444 | 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 below |
---|
| 446 | 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; |
---|
| 461 | |
---|
| 462 | // adding variables for the pencil beam. Variables in the absolute reference frame. |
---|
| 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; |
---|
| 471 | double totcount(0); |
---|
| 472 | |
---|
| 473 | double mirror; |
---|
| 474 | double tiltangle; |
---|
| 475 | |
---|
| 476 | |
---|
| 477 | double c_length; //Length in m |
---|
| 478 | double c_rotation(C_rotation) ; //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED |
---|
| 479 | double c_aperture(C_aperture) ; //Aperture in m |
---|
| 480 | double c_offset(C_offset) ; //Offset in m |
---|
| 481 | double c_tilt[2] = {C_tilt, 0} ; //Tilt in radian |
---|
| 482 | double c_tilt0[2] ; //Tilt in radian |
---|
| 483 | |
---|
| 484 | |
---|
| 485 | double xp_tangent; |
---|
| 486 | |
---|
| 487 | double Cry_tilt(Crystal_tilt); //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++//////////////// |
---|
| 488 | |
---|
| 489 | |
---|
| 490 | |
---|
| 491 | //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS------------------------------------------------------------------------- |
---|
| 492 | c_length = crys.Cry_length; |
---|
| 493 | |
---|
| 494 | if (crys.IS == 0) { |
---|
| 495 | mat = 8; |
---|
| 496 | } else if (crys.IS == 1) { |
---|
| 497 | mat = 9; |
---|
| 498 | } else if (crys.IS == 2) { |
---|
| 499 | mat = 10; |
---|
| 500 | } else if (crys.IS == 3) { |
---|
| 501 | mat = 11; |
---|
| 502 | } else { |
---|
| 503 | cout << "Material of the Crystal not fund !!!" << endl; |
---|
| 504 | } |
---|
| 505 | |
---|
| 506 | |
---|
| 507 | if (c_length > 0 ) { |
---|
| 508 | p0 = part.PC; |
---|
| 509 | c_tilt0[0] = c_tilt[0]; |
---|
| 510 | c_tilt0[1] = c_tilt[1]; |
---|
| 511 | tiltangle = c_tilt0[0]; |
---|
| 512 | } |
---|
| 513 | |
---|
| 514 | |
---|
| 515 | |
---|
| 516 | //ofstream sortie("results.csv"); |
---|
| 517 | //sortie << "xp1" <<","<< "xp0"<<","<< "x"<<","<< "y"<<endl; |
---|
| 518 | ofstream sortieIco; |
---|
| 519 | string name1; |
---|
| 520 | name1 = outputpath + "/ascii_after.csv"; |
---|
| 521 | |
---|
| 522 | sortieIco.open(name1.c_str()); |
---|
| 523 | |
---|
| 524 | string rest; |
---|
| 525 | |
---|
| 526 | ifstream entree; |
---|
| 527 | string name2; |
---|
| 528 | name2 = outputpath + "/ascii_before.csv"; |
---|
| 529 | |
---|
| 530 | entree.open(name2.c_str()); |
---|
| 531 | if (entree) { |
---|
| 532 | while (!entree.eof()) { |
---|
| 533 | count = count + 1; |
---|
| 534 | |
---|
| 535 | getline(entree, rest, ','); |
---|
| 536 | part.x = atof(rest.c_str()); |
---|
| 537 | |
---|
| 538 | |
---|
| 539 | //cout<<" X " << part.x<<endl; |
---|
| 540 | |
---|
| 541 | getline(entree, rest, ','); |
---|
| 542 | part.y = atof(rest.c_str()); |
---|
| 543 | |
---|
| 544 | |
---|
| 545 | //cout<<" Y " << part.y<<endl; |
---|
| 546 | |
---|
| 547 | getline(entree, rest, ','); |
---|
| 548 | part.xp = atof(rest.c_str()); |
---|
| 549 | |
---|
| 550 | //cout<<" XP " << part.xp<<endl; |
---|
| 551 | |
---|
| 552 | |
---|
| 553 | |
---|
| 554 | getline(entree, rest, ','); |
---|
| 555 | part.yp = atof(rest.c_str()); |
---|
| 556 | |
---|
| 557 | |
---|
| 558 | //cout<<" YP " << part.yp<<endl; |
---|
| 559 | |
---|
| 560 | getline(entree, rest); |
---|
| 561 | part.PC = atof(rest.c_str()); |
---|
| 562 | |
---|
| 563 | if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT |
---|
| 564 | |
---|
| 565 | av_pc0 = av_pc0 + part.PC; |
---|
| 566 | |
---|
| 567 | //part.x=random_gauss()*0.00115; |
---|
| 568 | //part.xp=random_gauss()*sqrt(5*1e-6); |
---|
| 569 | //cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl; |
---|
| 570 | |
---|
| 571 | // IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 |
---|
| 572 | //bdelta=tan(part.xp)*10; |
---|
| 573 | part.x = (part.x); //-0.0045; |
---|
| 574 | |
---|
| 575 | |
---|
| 576 | //bdelta=tan(part.yp)*10; |
---|
| 577 | part.y = (part.y); //-0.0045; |
---|
| 578 | |
---|
| 579 | |
---|
| 580 | |
---|
| 581 | //cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl; |
---|
| 582 | xp0 = part.xp; |
---|
| 583 | |
---|
| 584 | //part.y=random_gauss()*sqrt(0.00095); |
---|
| 585 | //part.yp=5*1e-6; |
---|
| 586 | |
---|
| 587 | |
---|
| 588 | //part.PC=120.976; |
---|
| 589 | //cout<<" X :"<<part.x<<" XP :"<<part.xp<<" Y :"<<part.y<<" YP :"<<part.yp<<" PC :"<<part.PC<<endl; |
---|
| 590 | AV_xp0 = AV_xp0 + part.xp; //average x angle before impact |
---|
| 591 | AV_yp0 = AV_yp0 + part.yp; //average y angle before impact |
---|
| 592 | RMSxp0 = RMSxp0 + part.xp * part.xp; //rms x angle before impact |
---|
| 593 | RMSyp0 = RMSyp0 + part.yp * part.yp; //rms y angle before impact |
---|
| 594 | |
---|
| 595 | |
---|
| 596 | // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 597 | |
---|
| 598 | |
---|
| 599 | |
---|
| 600 | s = 0; |
---|
| 601 | x = part.x; |
---|
| 602 | xp = part.xp; |
---|
| 603 | z = part.y; |
---|
| 604 | zp = part.yp; |
---|
| 605 | p = part.PC; |
---|
| 606 | |
---|
| 607 | x_temp = 0; |
---|
| 608 | x_int = 0; |
---|
| 609 | x_rot = 0; |
---|
| 610 | x_shift = 0; |
---|
| 611 | s_temp = 0; |
---|
| 612 | s_int = 0; |
---|
| 613 | s_rot = 0; |
---|
| 614 | s_shift = 0; |
---|
| 615 | xp_int = 0; |
---|
| 616 | xp_temp = 0; |
---|
| 617 | xp_rot = 0; |
---|
| 618 | xp_shift = 0; |
---|
| 619 | shift = 0; |
---|
| 620 | tilt_int = 0; |
---|
| 621 | dpop = (p - p0) / p0; |
---|
| 622 | |
---|
| 623 | // ------------------------------transform particle coordinates to get into collimator coordinate system ----------------------------------------- |
---|
| 624 | |
---|
| 625 | |
---|
| 626 | // first check whether particle was lost before |
---|
| 627 | totcount = totcount + 1; |
---|
| 628 | if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) { ////////////////////////// IF NUMERO 1////////////////// |
---|
| 629 | |
---|
| 630 | // first do rotation into collimator frame |
---|
| 631 | |
---|
| 632 | x = part.x * cos(c_rotation) + sin(c_rotation) * part.y; |
---|
| 633 | z = part.y * cos(c_rotation) - sin(c_rotation) * part.x; |
---|
| 634 | xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp; |
---|
| 635 | zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp; |
---|
| 636 | |
---|
| 637 | //cout<<"PARTICULE NON PERDU>>><<"<<endl; |
---|
| 638 | |
---|
| 639 | mirror = 1; |
---|
| 640 | x = mirror * x; |
---|
| 641 | xp = mirror * xp; |
---|
| 642 | |
---|
| 643 | //Shift with opening and offset |
---|
| 644 | |
---|
| 645 | x = x - c_aperture / 2 - c_offset; |
---|
| 646 | |
---|
| 647 | // Include collimator tilt |
---|
| 648 | |
---|
| 649 | if (tiltangle > 0) { ////////////////////////// IF NUMERO 2////////////////// |
---|
| 650 | xp = xp - tiltangle; |
---|
| 651 | } //////////////////////////FIN IF NUMERO 2////////////////// |
---|
| 652 | else if (tiltangle < 0) { ////////////////////////// IF NUMERO 3////////////////// |
---|
| 653 | x = x + sin(tiltangle) * c_length; |
---|
| 654 | xp = xp - tiltangle; |
---|
| 655 | } //////////////////////////FIN IF NUMERO 3////////////////// |
---|
| 656 | |
---|
| 657 | |
---|
| 658 | |
---|
| 659 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 4////////////////// |
---|
| 660 | |
---|
| 661 | s_shift = s; |
---|
| 662 | shift = crys.Rcurv * (1 - cos(Cry_tilt)); |
---|
| 663 | if (Cry_tilt < (-crys.cry_bend) ) { ////////////////////////// IF NUMERO 5////////////////// |
---|
| 664 | shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) ); |
---|
| 665 | } //////////////////////////FIN IF NUMERO 5////////////////// |
---|
| 666 | x_shift = x - shift; |
---|
| 667 | } //////////////////////////fin IF NUMERO 4////////////////// |
---|
| 668 | else { //////////////////////////IF NUMERO 6////////////////// |
---|
| 669 | //cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl; |
---|
| 670 | s_shift = s; |
---|
| 671 | x_shift = x; |
---|
| 672 | } //////////////////////////fin IF NUMERO 6////////////////// |
---|
| 673 | // cout<<"debug - S shift" << s_shift<<endl; |
---|
| 674 | //cout<<"debug - X shift" << x_shift <<endl; |
---|
| 675 | // |
---|
| 676 | // 2nd transformation: rotation |
---|
| 677 | s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt); |
---|
| 678 | x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt); |
---|
| 679 | xp_rot = xp - Cry_tilt; |
---|
| 680 | //cout<< "debug - S rot" << s_rot <<endl; |
---|
| 681 | //cout<<"debug - X rot" << x_rot <<endl; |
---|
| 682 | //cout<<"debug - XP rot" << xp_rot<<endl; |
---|
| 683 | |
---|
| 684 | // 3rd transformation: drift to the new coordinate s=0 |
---|
| 685 | xp = xp_rot; |
---|
| 686 | x = x_rot - xp_rot * s_rot; |
---|
| 687 | z = z - zp * s_rot; |
---|
| 688 | s = 0; |
---|
| 689 | //cout<<"debug - S cryRF" << s_rot <<endl; |
---|
| 690 | //cout<<"debug - X cryRF" << x_rot <<endl; |
---|
| 691 | //cout<<"debug - XP cryRF" << xp_rot <<endl; |
---|
| 692 | |
---|
| 693 | |
---|
| 694 | //----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal--------------------------------------------------------------- |
---|
| 695 | //cout<<" X " << x<<endl; |
---|
| 696 | if (x >= 0) { ////////////////////////// IF NUMERO 7////////////////// |
---|
| 697 | s_impact = s_in0; //!(for the first impact) |
---|
| 698 | //cout<<"hit the cry entrance face"<<endl; |
---|
| 699 | //cout<<"impact at s,x = "<< s_impact<<x_in0<<endl; |
---|
| 700 | //cout<<"with angle xp = "<<xp<<endl; |
---|
| 701 | //cout<<"s before"<< s<<endl; |
---|
| 702 | |
---|
| 703 | part.x = x; |
---|
| 704 | part.xp = xp; |
---|
| 705 | part.y = z; |
---|
| 706 | part.yp = zp; |
---|
| 707 | part.PC = p; |
---|
| 708 | //cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl; |
---|
| 709 | //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
| 710 | bases(sim); |
---|
| 711 | //cout << part.xp<< " XP "<<endl; |
---|
| 712 | bool crossing(cross(sim)); |
---|
| 713 | //cout<<part.xp-xp0<<endl; |
---|
| 714 | //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
| 715 | bool layering; |
---|
| 716 | if (crossing != 0) { ////////////////////////// IF NUMERO 8////////////////// |
---|
| 717 | layering = layer(sim); |
---|
| 718 | //cout<<"Crossing !!!!"<<endl; |
---|
| 719 | |
---|
| 720 | if (layering == true) { ////////////////////////// IF NUMERO 9////////////////// |
---|
| 721 | parameters(sim); |
---|
| 722 | //cout<<" Layering =true " <<endl; |
---|
| 723 | |
---|
| 724 | if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10////////////////// |
---|
| 725 | channel(sim); |
---|
| 726 | //cout<< " CHANNELING " << endl; |
---|
| 727 | } ////////////////////////// fin IF NUMERO 10////////////////// |
---|
| 728 | else { ////////////////////////// IF NUMERO 11////////////////// |
---|
| 729 | reflection(sim); |
---|
| 730 | //cout<<" REFLECTION "<<endl; |
---|
| 731 | } ////////////////////////// fin IF NUMERO 11////////////////// |
---|
| 732 | |
---|
| 733 | |
---|
| 734 | } ////////////////////////// fin IF NUMERO 9////////////////// |
---|
| 735 | } ////////////////////////// fin IF NUMERO 8////////////////// |
---|
| 736 | |
---|
| 737 | xp = part.xp; |
---|
| 738 | |
---|
| 739 | //} |
---|
| 740 | |
---|
| 741 | s = crys.Rcurv * sin(crys.cry_bend); |
---|
| 742 | zlm = crys.Rcurv * sin(crys.cry_bend); |
---|
| 743 | //cout<<"process:"<<PROC<<endl; |
---|
| 744 | //cout<<"s after"<< s<<endl; |
---|
| 745 | //cout<<"part.xp"<<part.xp<<endl; |
---|
| 746 | |
---|
| 747 | //cout<<"----------------------1111111111111111111111111-----------------"<<endl; |
---|
| 748 | |
---|
| 749 | |
---|
| 750 | } ////////////////////////// fin IF NUMERO 7////////////////// |
---|
| 751 | else { ////////////////////////// IF NUMERO 12////////////////// |
---|
| 752 | |
---|
| 753 | //cout<<"OU LA ABBBBBBBBBBBB" <<endl; |
---|
| 754 | xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv)); |
---|
| 755 | //cout<<"RCURV "<<crys.Rcurv<<endl; |
---|
| 756 | //cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl; |
---|
| 757 | |
---|
| 758 | //cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl; |
---|
| 759 | //cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl; |
---|
| 760 | if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13////////////////// |
---|
| 761 | |
---|
| 762 | //if it hits the crystal, calculate in which point and apply the |
---|
| 763 | //transformation and drift to that point |
---|
| 764 | a_eq = (1 + xp * xp); |
---|
| 765 | b_eq = 2 * xp * (x - crys.Rcurv); |
---|
| 766 | c_eq = -2 * x * crys.Rcurv + x * x; |
---|
| 767 | Delta = b_eq * b_eq - 4 * (a_eq * c_eq); |
---|
| 768 | s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq); |
---|
| 769 | //cout<< " s int "<<s_int <<endl; |
---|
| 770 | if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14////////////////// |
---|
| 771 | // transform to a new ref system:shift and rotate |
---|
| 772 | |
---|
| 773 | x_int = xp * s_int + x; |
---|
| 774 | xp_int = xp; |
---|
| 775 | z = z + zp * s_int; |
---|
| 776 | x = 0; |
---|
| 777 | s = 0; |
---|
| 778 | |
---|
| 779 | |
---|
| 780 | // tilt_int=2*X_int/S_int |
---|
| 781 | tilt_int = s_int / crys.Rcurv; |
---|
| 782 | xp = xp - tilt_int; |
---|
| 783 | |
---|
| 784 | |
---|
| 785 | |
---|
| 786 | part.x = x; |
---|
| 787 | part.xp = xp; |
---|
| 788 | part.y = z; |
---|
| 789 | part.yp = zp; |
---|
| 790 | part.PC = p; |
---|
| 791 | crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv); |
---|
| 792 | bases(sim); |
---|
| 793 | |
---|
| 794 | //cout<<"part.xp"<<part.xp<<endl; |
---|
| 795 | |
---|
| 796 | bool crossing(cross(sim)); |
---|
| 797 | //cout<<"Crossing !!!!"<<endl; |
---|
| 798 | bool layering(layer(sim)); |
---|
| 799 | |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | if (layering == true) { |
---|
| 803 | //cout<<" Layering =true " <<endl; |
---|
| 804 | parameters(sim); |
---|
| 805 | if ((abs(part.xp_rel) < crys.xpcrit)) { |
---|
| 806 | channel(sim); |
---|
| 807 | //cout<< " CHANNELING " << endl; |
---|
| 808 | } else { |
---|
| 809 | reflection(sim); |
---|
| 810 | //cout<< " REFLECTION " << endl; |
---|
| 811 | } |
---|
| 812 | |
---|
| 813 | |
---|
| 814 | } |
---|
| 815 | |
---|
| 816 | //cout<<"part.xp"<<part.xp<<endl; |
---|
| 817 | //cout<<"----------------------2222222222222222222-----------------"<<endl; |
---|
| 818 | |
---|
| 819 | s = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
| 820 | zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | if (crossing != 0) { ////////////////////////// IF NUMERO 15////////////////// |
---|
| 824 | x_rot = x_int; |
---|
| 825 | s_rot = s_int; |
---|
| 826 | xp_rot = xp_int; |
---|
| 827 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
| 828 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
| 829 | xp_shift = xp_rot + Cry_tilt; |
---|
| 830 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16////////////////// |
---|
| 831 | s_impact = s_shift; |
---|
| 832 | x_in0 = x_shift + shift; |
---|
| 833 | xp_in0 = xp_shift; |
---|
| 834 | } ////////////////////////// fin IF NUMERO 16////////////////// |
---|
| 835 | else { ////////////////////////// IF NUMERO 17////////////////// |
---|
| 836 | x_in0 = x_shift; |
---|
| 837 | s_impact = s_shift; |
---|
| 838 | xp_in0 = xp_shift; |
---|
| 839 | } ////////////////////////// fin IF NUMERO 17////////////////// |
---|
| 840 | |
---|
| 841 | } ////////////////////////// fin IF NUMERO 15////////////////// |
---|
| 842 | |
---|
| 843 | |
---|
| 844 | x_temp = x; |
---|
| 845 | s_temp = s; |
---|
| 846 | xp_temp = xp; |
---|
| 847 | s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int); |
---|
| 848 | x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int); |
---|
| 849 | xp = xp_temp + tilt_int; |
---|
| 850 | |
---|
| 851 | // 2nd: shift back the 2 axis |
---|
| 852 | x = x + x_int; |
---|
| 853 | s = s + s_int; |
---|
| 854 | |
---|
| 855 | } ////////////////////////// fin IF NUMERO 14////////////////// |
---|
| 856 | else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"///////////////////////// |
---|
| 857 | |
---|
| 858 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
| 859 | x = x + s * xp; |
---|
| 860 | z = z + s * zp; |
---|
| 861 | ++part.nout; |
---|
| 862 | //cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl; |
---|
| 863 | } |
---|
| 864 | } else { |
---|
| 865 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
| 866 | x = x + s * xp; |
---|
| 867 | z = z + s * zp; |
---|
| 868 | ++part.nout; |
---|
| 869 | //cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl; |
---|
| 870 | } |
---|
| 871 | //} ////////////////////////// fin IF NUMERO 13////////////////// |
---|
| 872 | |
---|
| 873 | }////////////////////// 12 |
---|
| 874 | |
---|
| 875 | |
---|
| 876 | //trasform back from the crystal to the collimator reference system |
---|
| 877 | // 1st: un-rotate the coordinates |
---|
| 878 | |
---|
| 879 | |
---|
| 880 | x_rot = x; |
---|
| 881 | s_rot = s; |
---|
| 882 | xp_rot = xp; |
---|
| 883 | //cout<< "debug - S Cry RF 2" , s_rot <<endl; |
---|
| 884 | //cout<<"debug - X Cry RF 2" , x_rot <<endl; |
---|
| 885 | //cout<<"debug - XP Cry RF 2" , xp_rot <<endl; |
---|
| 886 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
| 887 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
| 888 | xp_shift = xp_rot + Cry_tilt; |
---|
| 889 | // 2nd: shift back the reference frame |
---|
| 890 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18////////////////// |
---|
| 891 | s = s_shift; |
---|
| 892 | x = x_shift + shift; |
---|
| 893 | xp = xp_shift; |
---|
| 894 | } ////////////////////////// fin IF NUMERO 18////////////////// |
---|
| 895 | else { |
---|
| 896 | x = x_shift; |
---|
| 897 | s = s_shift; |
---|
| 898 | xp = xp_shift; |
---|
| 899 | } |
---|
| 900 | // 3rd: shift to new S=Length position |
---|
| 901 | x = xp * (c_length - s) + x; |
---|
| 902 | z = zp * (c_length - s) + z; |
---|
| 903 | s = c_length; |
---|
| 904 | |
---|
| 905 | |
---|
| 906 | //cout<< "debug - S Coll RF 2" , s_rot <<endl; |
---|
| 907 | //cout<<"debug - X Coll RF 2" , x_rot <<endl; |
---|
| 908 | //cout<<"debug - XP Coll RF 2" , xp_rot <<endl; |
---|
| 909 | |
---|
| 910 | //cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl; |
---|
| 911 | |
---|
| 912 | if (part.effet == 2) { |
---|
| 913 | part_abs = 0; |
---|
| 914 | } else { |
---|
| 915 | part_abs = 1; |
---|
| 916 | } |
---|
| 917 | |
---|
| 918 | |
---|
| 919 | //========================================================================================================================= |
---|
| 920 | // Transform back to particle coordinates with opening and offset |
---|
| 921 | |
---|
| 922 | if (part_abs == 0) { ////////////////////////// IF NUMERO 19////////////////// |
---|
| 923 | // |
---|
| 924 | // Include collimator tilt |
---|
| 925 | |
---|
| 926 | if (tiltangle > 0) { |
---|
| 927 | x = x + tiltangle * c_length; |
---|
| 928 | xp = xp + tiltangle; |
---|
| 929 | } else if (tiltangle < 0) { |
---|
| 930 | x = x + tiltangle * c_length; |
---|
| 931 | xp = xp + tiltangle; |
---|
| 932 | |
---|
| 933 | x = x - sin(tiltangle) * c_length; |
---|
| 934 | } |
---|
| 935 | |
---|
| 936 | // Transform back to particle coordinates with opening and offset |
---|
| 937 | |
---|
| 938 | z00 = z; |
---|
| 939 | x00 = x + mirror * c_offset; |
---|
| 940 | x = x + c_aperture / 2 + mirror * c_offset; |
---|
| 941 | |
---|
| 942 | //+ Now mirror at the horizontal axis for negative X offset |
---|
| 943 | |
---|
| 944 | x = mirror * x; |
---|
| 945 | xp = mirror * xp; |
---|
| 946 | |
---|
| 947 | // Last do rotation into collimator frame |
---|
| 948 | |
---|
| 949 | part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation); |
---|
| 950 | part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation); |
---|
| 951 | part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation); |
---|
| 952 | part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation); |
---|
| 953 | |
---|
| 954 | |
---|
| 955 | //p_in = (1 + dpop) * p0; |
---|
| 956 | s_in = s_in + s; |
---|
| 957 | |
---|
| 958 | if (part.effet == 1) { ////////////////////////// IF NUMERO 21////////////////// |
---|
| 959 | |
---|
| 960 | part.x = 99.99e-3; |
---|
| 961 | part.y = 99.99e-3; |
---|
| 962 | //cout<< "Mean that the particules is lost"<<endl; |
---|
| 963 | } ////////////////////////// fin IF NUMERO 21////////////////// |
---|
| 964 | } ////////////////////////// fin IF NUMERO 19////////////////// |
---|
| 965 | |
---|
| 966 | |
---|
| 967 | // End of check for particles not being lost before (see @330) |
---|
| 968 | |
---|
| 969 | //} ////////////////////////// fin IF NUMERO 12////////////////// |
---|
| 970 | // End of loop over all particles |
---|
| 971 | |
---|
| 972 | |
---|
| 973 | |
---|
| 974 | } //collimator with length = 0 ////////////////////////// fin IF NUMERO 1////////////////// |
---|
| 975 | |
---|
| 976 | // =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM==================================================================== |
---|
| 977 | /* |
---|
| 978 | part.xp=xp; |
---|
| 979 | part.yp=zp; |
---|
| 980 | part.y=z; |
---|
| 981 | part.x=x; |
---|
| 982 | part.PC=p; |
---|
| 983 | */ |
---|
| 984 | xp1 = part.xp; |
---|
| 985 | |
---|
| 986 | |
---|
| 987 | //if(xp1-xp0!=0){ |
---|
| 988 | //sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl; |
---|
| 989 | //} |
---|
| 990 | //cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl; |
---|
| 991 | |
---|
| 992 | |
---|
| 993 | AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction |
---|
| 994 | AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction |
---|
| 995 | RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction |
---|
| 996 | RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction |
---|
| 997 | av_pc1 = av_pc1 + part.PC; |
---|
| 998 | |
---|
| 999 | |
---|
| 1000 | //cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl; |
---|
| 1001 | //cout<<endl; |
---|
| 1002 | //sortieIco.setf(ios::scientific); |
---|
| 1003 | sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl; |
---|
| 1004 | //}//end of for.... |
---|
| 1005 | |
---|
| 1006 | } |
---|
| 1007 | } //End of while |
---|
| 1008 | |
---|
| 1009 | |
---|
| 1010 | } else { |
---|
| 1011 | //cout<< "We can not open the file !" << endl; |
---|
| 1012 | } |
---|
| 1013 | entree.close(); |
---|
| 1014 | |
---|
| 1015 | |
---|
| 1016 | |
---|
| 1017 | AV_xp0 = AV_xp0 / count; |
---|
| 1018 | AV_yp0 = AV_yp0 / count; |
---|
| 1019 | |
---|
| 1020 | RMSxp0 = sqrt(RMSxp0 / (count)); |
---|
| 1021 | RMSyp0 = sqrt(RMSyp0 / (count)); |
---|
| 1022 | |
---|
| 1023 | SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0)); |
---|
| 1024 | SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0)); |
---|
| 1025 | |
---|
| 1026 | AV_xp1 = AV_xp1 / count; |
---|
| 1027 | AV_yp1 = AV_yp1 / count; |
---|
| 1028 | RMSxp1 = sqrt(RMSxp1 / (count)); |
---|
| 1029 | RMSyp1 = sqrt(RMSyp1 / (count)); |
---|
| 1030 | |
---|
| 1031 | SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1)); |
---|
| 1032 | SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1)); |
---|
| 1033 | av_pc0 = av_pc0 / count; |
---|
| 1034 | av_pc1 = av_pc1 / count; |
---|
| 1035 | //-------------------------------------------------------------------------- |
---|
| 1036 | //-------------write a summary file----------------------------------------- |
---|
| 1037 | //cout<<endl; |
---|
| 1038 | //cout<<endl; |
---|
| 1039 | //cout<<endl; |
---|
| 1040 | |
---|
| 1041 | //cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl; |
---|
| 1042 | //cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl; |
---|
| 1043 | //cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl; |
---|
| 1044 | //cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl; |
---|
| 1045 | |
---|
| 1046 | |
---|
| 1047 | |
---|
| 1048 | sortieIco.close(); |
---|
| 1049 | //cout<<endl; |
---|
| 1050 | //affiche(); |
---|
| 1051 | file_out(pas, outputpath); |
---|
| 1052 | |
---|
| 1053 | //cout<<endl; |
---|
| 1054 | //cout<<endl; |
---|
| 1055 | //cout<<endl; |
---|
| 1056 | } |
---|
| 1057 | |
---|
| 1058 | |
---|
| 1059 | |
---|
| 1060 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1061 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1062 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1063 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1064 | |
---|
| 1065 | //WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW |
---|
| 1066 | |
---|
| 1067 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1068 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1069 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1070 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1071 | |
---|
| 1072 | |
---|
| 1073 | void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC) |
---|
| 1074 | { |
---|
| 1075 | //cout<<"The move_am function is use !!!!! "<<endl; |
---|
| 1076 | double DYA(0); |
---|
| 1077 | double Wp(0); |
---|
| 1078 | if (NAM == 0) { |
---|
| 1079 | return; |
---|
| 1080 | } |
---|
| 1081 | xp = xp * 1000; |
---|
| 1082 | yp = yp * 1000; |
---|
| 1083 | //cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl; |
---|
| 1084 | //DEI ---> dE/dx (stoping energy) |
---|
| 1085 | PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV] |
---|
| 1086 | |
---|
| 1087 | // Coulomb Scattering |
---|
| 1088 | //DYA ---> rms of coloumb scattering |
---|
| 1089 | DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad) |
---|
| 1090 | xp = xp + DYA * random_gauss_cut(); |
---|
| 1091 | yp = yp + DYA * random_gauss_cut(); |
---|
| 1092 | |
---|
| 1093 | //Elastic scattering |
---|
| 1094 | //cout<<"Plusieurs possibilitees : "<<endl; |
---|
| 1095 | if (random() <= (DZ / crys.DLAI[IS])) { |
---|
| 1096 | //cout<<"Poss. 1 !!!! "<<endl; |
---|
| 1097 | xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane |
---|
| 1098 | yp = yp + 196 * random_gauss_cut() / PC; |
---|
| 1099 | } |
---|
| 1100 | |
---|
| 1101 | //Diffraction interaction |
---|
| 1102 | if (random() <= (DZ / (DLY * 6.143))) { |
---|
| 1103 | //cout<<"Poss. 2 !!!! "<<endl; |
---|
| 1104 | xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr] |
---|
| 1105 | yp = yp + 257 * random_gauss_cut() / PC; |
---|
| 1106 | Wp = random(); |
---|
| 1107 | PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle |
---|
| 1108 | |
---|
| 1109 | } |
---|
| 1110 | |
---|
| 1111 | //Inelastic interaction |
---|
| 1112 | if (random() <= DZ / DLY) { |
---|
| 1113 | //cout<<" Absorbed !! "<< endl; |
---|
| 1114 | } |
---|
| 1115 | //cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl; |
---|
| 1116 | xp = xp / 1000; |
---|
| 1117 | yp = yp / 1000; |
---|
| 1118 | } |
---|
| 1119 | |
---|
| 1120 | |
---|
| 1121 | double SimCrys::random() |
---|
| 1122 | { |
---|
| 1123 | double scale = RAND_MAX; |
---|
| 1124 | double base = rand() / scale; |
---|
| 1125 | double fine = rand() / scale; |
---|
| 1126 | return base + fine / scale; |
---|
| 1127 | } |
---|
| 1128 | |
---|
| 1129 | double SimCrys::random_dist() |
---|
| 1130 | { |
---|
| 1131 | double scale = RAND_MAX; |
---|
| 1132 | double base = rand() / scale; |
---|
| 1133 | double fine = rand() / scale; |
---|
| 1134 | return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale)); |
---|
| 1135 | } |
---|
| 1136 | |
---|
| 1137 | double SimCrys::random_gauss() |
---|
| 1138 | { |
---|
| 1139 | const double PI (4.0 * atan(1.0)); |
---|
| 1140 | double U(random()); |
---|
| 1141 | double V(random()); |
---|
| 1142 | return (sqrt(-2 * log(U)) * cos(2 * PI * V)); |
---|
| 1143 | } |
---|
| 1144 | |
---|
| 1145 | double SimCrys::random_gauss_cut() |
---|
| 1146 | { |
---|
| 1147 | double resu; |
---|
| 1148 | do { |
---|
| 1149 | resu = random_gauss(); |
---|
| 1150 | } while (abs(resu) >= 1); |
---|
| 1151 | return resu; |
---|
| 1152 | } |
---|
| 1153 | |
---|
| 1154 | double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
| 1155 | { |
---|
| 1156 | const double PI (4.0 * atan(1.0)); |
---|
| 1157 | double U(random()); |
---|
| 1158 | double V(random()); |
---|
| 1159 | //moy=random(); |
---|
| 1160 | return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 1161 | } |
---|
| 1162 | |
---|
| 1163 | /////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1164 | //POUR 409 : |
---|
| 1165 | |
---|
| 1166 | double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
| 1167 | { |
---|
| 1168 | const double PI (4.0 * atan(1.0)); |
---|
| 1169 | double U(random()); |
---|
| 1170 | double V(random()); |
---|
| 1171 | //moy=random(); |
---|
| 1172 | return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
| 1173 | } |
---|
| 1174 | |
---|
| 1175 | ///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1176 | |
---|
| 1177 | void SimCrys::affiche() |
---|
| 1178 | { |
---|
| 1179 | cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
| 1180 | cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
| 1181 | cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1182 | cout << "Critical Radius :" << crys.Rcrit << endl; |
---|
| 1183 | cout << "Ratio :" << crys.ratio << endl; |
---|
| 1184 | cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
| 1185 | cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
| 1186 | cout << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
| 1187 | cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1188 | cout << endl; |
---|
| 1189 | cout << endl; |
---|
| 1190 | cout << "N.AMORPH :" << part.namor << endl; |
---|
| 1191 | cout << "N.DECHANNELING:" << part.ndechann << endl; |
---|
| 1192 | cout << "N.CHANNELING :" << part.nchann << endl; |
---|
| 1193 | cout << "N.OUT :" << part.nout << endl; |
---|
| 1194 | cout << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
| 1195 | cout << "N.VCAPTURE :" << part.nvcapt << endl; |
---|
| 1196 | |
---|
| 1197 | //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!"; |
---|
| 1198 | } |
---|
| 1199 | |
---|
| 1200 | void SimCrys::file_out(const int& pas, string outputpath) |
---|
| 1201 | { |
---|
| 1202 | |
---|
| 1203 | ofstream out; |
---|
| 1204 | string file; |
---|
| 1205 | file = outputpath + "/crystal_output.out"; |
---|
| 1206 | |
---|
| 1207 | out.open(file.c_str(), ios::out | ios::app); |
---|
| 1208 | |
---|
| 1209 | if (out.fail()) { |
---|
| 1210 | cerr << "Warning: problem writing in the file " << file << "!" << endl; |
---|
| 1211 | } else { |
---|
| 1212 | |
---|
| 1213 | out << endl; |
---|
| 1214 | out << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
| 1215 | out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
| 1216 | out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1217 | out << "Critical Radius :" << crys.Rcrit << endl; |
---|
| 1218 | out << "Ratio :" << crys.ratio << endl; |
---|
| 1219 | out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
| 1220 | out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
| 1221 | out << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
| 1222 | out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
| 1223 | out << endl; |
---|
| 1224 | out << endl; |
---|
| 1225 | out << "N.AMORPH :" << part.namor << endl; |
---|
| 1226 | out << "N.DECHANNELING:" << part.ndechann << endl; |
---|
| 1227 | out << "N.CHANNELING :" << part.nchann << endl; |
---|
| 1228 | out << "N.OUT :" << part.nout << endl; |
---|
| 1229 | out << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
| 1230 | out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl; |
---|
| 1231 | |
---|
| 1232 | out.close(); |
---|
| 1233 | } |
---|
| 1234 | |
---|
| 1235 | } |
---|
| 1236 | |
---|