Pluto v6.3.0
Loading...
Searching...
No Matches
piony.C
Go to the documentation of this file.
1#include <cstdlib>
2#include <string>
3#include <cstring>
4#include <cstdio>
5#include <ctime>
6#include <iomanip>
7#include <cstdlib>
8#include <cstring>
9#include <cstdio>
10#include <ctime>
11#include <iostream>
12#include <math.h>
13#include <fstream>
14#include "iostream"
15#include "math.h"
16#include "stdio.h"
17#include "stdlib.h"
18#include "time.h"
19#include "string"
20#include "unistd.h"
21#include "string.h"
22#include "fstream"
23//#include "Math.h"
24#include "Math/Math.h"
25#include "TROOT.h"
26#include "TTree.h"
27#include <TFile.h>
28#include "TChain.h"
29#include "TCanvas.h"
30#include "TLine.h"
31#include "TArc.h"
32#include "TF1.h"
33#include "TMath.h"
34#include "TFile.h"
35#include "TMath.h"
36#include "TH1F.h"
37
38//#include "header.h"
39
40using namespace std;
41
42int piony()
43
44{
45
46 long double T11[51],T12[51],T13[51],T14[51],T16[51];
47 long double T111[51],T112[51],T113[51],T114[51],T116[51];
48 long double T122[51],T123[51],T124[51],T126[51];
49 long double T133[51],T134[51],T136[51];
50 long double T144[51],T146[51];
51 long double T166[51],along[51];
52
53 long double T21[51],T22[51],T23[51],T24[51],T26[51];
54 long double T211[51],T212[51],T213[51],T214[51],T216[51];
55 long double T222[51],T223[51],T224[51],T226[51];
56 long double T233[51],T234[51],T236[51];
57 long double T244[51],T246[51];
58 long double T266[51];
59
60 long double T31[51],T32[51],T33[51],T34[51],T36[51];
61 long double T311[51],T312[51],T313[51],T314[51],T316[51];
62 long double T322[51],T323[51],T324[51],T326[51];
63 long double T333[51],T334[51],T336[51];
64 long double T344[51],T346[51];
65 long double T366[51];
66
67 long double T41[51],T42[51],T43[51],T44[51],T46[51];
68 long double T411[51],T412[51],T413[51],T414[51],T416[51];
69 long double T422[51],T423[51],T424[51],T426[51];
70 long double T433[51],T434[51],T436[51];
71 long double T444[51],T446[51];
72 long double T466[51];
73
74 Int_t ipass[51];
75 Int_t kpass[51],rpass[51];
76 char tm[11];
77 char tc[5];
78 char tl[31];
79 char txt14[15];
80 char txt22[23];
81 char txt2[3];
82 char txt5[6];
83 char txt6[7];
84
85 long double rapport = 0;
86 long double offset = 0;
87 long double rx = 0;
88 Int_t itotpaass = 0;
89 Int_t itotq9 = 0;
90 Int_t idet1 = 0;
91 Int_t idet2 = 0;
92 long double dist1 = 0;
93 long double dist2 = 0;
94 long double ouvrescaleh = 0;
95 long double ouvrescalev = 0;
96 long double r1 = 0;
97 Int_t ii1 = 0;
98 Int_t ii2 = 0;
99 long double r2 = 0;
100 long double r3 = 0;
101 Int_t ii3 = 0;
102 Int_t ii4 = 0;
103 long double r4 = 0;
104 long double rh = 0;
105 long double rv = 0;
106 long double rres = 0;
107 long double rdm = 0;
108 long double shiftx = 0;
109 long double shifty = 0;
110 long double dx = 0;
111 long double dxp = 0;
112 long double dy = 0;
113 long double dyp = 0;
114 long double dpspmax = 0;
115 long double dpsppas = 0;
116 Int_t ntir = 0;
117 long double ppion = 0;
118 long double esi = 0;
119 long double resoldet = 0;
120 long double betapion = 0;
121 long double diffmul = 0;
122 long double resoldmx = 0;
123 long double resoldmy = 0;
124 Int_t idineg = 0;
125 long double xmax = 0;
126 long double ymax = 0;
127 long double xpmax = 0;
128 long double ypmax = 0;
129 long double dpmax = 0;
130 //Int_t k = 0;
131 long double relem = 0;
132 long double z1 = 0;
133 long double z2 = 0;
134 long double z3 = 0;
135 long double z4 = 0;
136 long double z5 = 0;
137 long double z6 = 0;
138 long double z7 = 0;
139 long double z8 = 0;
140 long double z9 = 0;
141 long double z10 = 0;
142 long double z11 = 0;
143 long double z12 = 0;
144 long double z13 = 0;
145 long double z14 = 0;
146 long double z15 = 0;
147 Int_t itir = 1;
148 long double yy = 0;
149 long double xi = 0;
150 long double xpi = 0;
151 long double yi = 0;
152 long double ypi = 0;
153 long double dp = 0;
154 long double dd = 0;
155 long double ndp = 0;
156 long double uns = 0;
157 Int_t kii = 0;
158 long double seuil = 0;
159 Int_t iok = 0;
160 long double xf = 0;
161 long double yf = 0;
162 long double xpf = 0;
163 long double ypf = 0;
164 long double xf1 = 0;
165 long double yf1 = 0;
166 long double xpf1 = 0;
167 long double ypf1 = 0;
168 long double xf2 = 0;
169 long double xpf2 = 0;
170 long double yf2 = 0;
171 long double ypf2 = 0;
172 long double xdet1 = 0;
173 long double ydet1 = 0;
174 long double xdet2 = 0;
175 long double ydet2 = 0;
176 Int_t kk = 0;
177 long double deltax = 0;
178 long double deltay = 0;
179 long double xhad = 0;
180 long double yhad = 0;
181 long double radius = 0;
182 Int_t ijk = 0;
183 Int_t itotpass = 0;
184 long double xdet1sav = 0;
185 long double xca = 0;
186 Int_t nca = 0;
187 long double xdet2sav = 0;
188 long double ydet1sav = 0;
189 long double ydet2sav = 0;
190 long double yca = 0;
191 long double alfa = 0;
192 long double beta = 0;
193 long double gamma = 0;
194 long double discri = 0;
195 long double galfa = 0;
196 long double gbeta = 0;
197 long double ggamma = 0;
198 long double gdelta = 0;
199 long double aa1 = 0;
200 long double aa2 = 0;
201 long double aa3 = 0;
202 long double aa4 = 0;
203 long double sol = 0;
204 Int_t iter = 1;
205 long double yyy1 = 0;
206 long double yyp1 = 0;
207 long double yyy2 = 0;
208 long double yyp2 = 0;
209 long double yyy3 = 0;
210 long double yyp3 = 0;
211 long double yyy4 = 0;
212 long double yyp4 = 0;
213 long double yyy5 = 0;
214 long double yyp5 = 0;
215 long double teta = 0;
216 long double bb = 0;
217 long double cc = 0;
218 long double aap = 0;
219 long double bbp = 0;
220 long double ccp = 0;
221 long double deter = 0;
222 long double yci = 0;
223 long double phi = 0;
224 long double poids = 0;
225 long double xfth = 0;
226 long double yfth = 0;
227 long double xfcal = 0;
228 long double yfcal = 0;
229 long double rr = 0;
230 Int_t ipi = 0;
231 Int_t ips2 = 0;
232 Int_t idif = 0;
233 long double xabs = 0;
234 long double yabs = 0;
235 Int_t ki = 0;
236 long double reste = 0;
237 long double ri = 0;
238 Int_t kij = 0;
239 long double alossfraction = 0;
240
241 // string txt2 = 0;
242 long double aa = 0;
243 long double sig11 = 0;
244// string txt5 = 0;
245 long double sig22 = 0;
246 long double r12 = 0;
247 long double sig33 = 0;
248 long double r13 = 0;
249 long double r23 = 0;
250 long double sig44 = 0;
251 long double r14 = 0;
252 long double r24 = 0;
253 long double r34 = 0;
254 long double sig55 = 0;
255 long double r15 = 0;
256 long double r25 = 0;
257 long double r35 = 0;
258 long double r45 = 0;
259 long double sig66 = 0;
260 long double r16 = 0;
261 long double r26 = 0;
262 long double r36 = 0;
263 long double r46 = 0;
264 long double r56 = 0;
265 // string txt14 = 0;
266 Int_t ii = 0;
267 Int_t ind1 = 0;
268 Int_t ind23 = 0;
269 double rreste = 0;
270 Int_t retape = 0;
271
272 Int_t ouvh[51] = {0,6.,6.,6.,6.,7.,7.,7.,7.,7.,6.,
273 6.,6.,6.,6.,6.,5.,5.,6.,6.,6.,
274 6.,6.,6.,6.,6.,6.,6.,6.,6.,6.,
275 6.,6.,0.6,5.,5.,5.,5.,5.,5.,5.,
276 5.,5.,5.,5.,5.,5.,5.,5.,5.,5. };
277 Int_t ouvv[51] = {0,6.,6.,6.,6.,3.5,3.5,3.5,3.5,3.5,3.5,
278 3.5,6.,6.,6.,6.,6.,4.,4.,6.,6.,
279 6.,3.5,3.5,6.,6.,6.,6.,6.,6.,6.,
280 6.,6.,0.6,5.,5.,5.,5.,5.,5.,5.,
281 5.,5.,5.,5.,5.,5.,5.,5.,5.,5.};
282 Int_t kind[11]= {0, 2, 4, 6, 10, 13, 17, 23, 26, 29, 33 };
283 // string str;
284 // std::cin.fixed;
285 // cin.fixed;
286 // cout.fixed;
287 std::ifstream plik1;
288 plik1.open("pibeam_set6.data");
289 // getline(plik1,str);
290 // std::cout<<"dupa"<<str<<endl;
291 //ifstream plik1("pibeam_set3.data");
292 std::ofstream plik2;
293 plik2.open("pibeam_set3.out");
294
295 // cout.precision(10);
296 // cin.precision(10);
297
298 //ios::precision();
299 // ifstream.precision(10);
300 // ofstream.precision (std::numeric_limits<double>::digits10 + 1);
301 Int_t nelem=34;
302 long double renorgauss=sqrt(3./5.);
303 Int_t x0si=22.;
304 long double rhosi=2.33;
305 long double amasspi=0.13958;
306 //Int_t dupa;
307 //char asdf[3];
308
309 for(Int_t iii=0; iii<51; iii++)
310 {
311 ipass[iii] = 0;
312 kpass[iii] = 0;
313 rpass[iii] = 0.;
314 }
315
316 //--TTree--- branch variables----------------------
317 Int_t sector =0;
318 double rT166 = 0,rT266 = 0,rT366 = 0,rT466 = 0,raa = 0,ralong = 0;
319
320 TFile *output = new TFile("piony.root", "RECREATE");
321
322 TTree *tree = new TTree("Pion","Pionbeam");
323 tree->Branch("sectors", &sector);
324 tree->Branch("T166", &rT166);
325 tree->Branch("T266", &rT266);
326 tree->Branch("T366", &rT366);
327 tree->Branch("T466", &rT466);
328 tree->Branch("aa", &raa);
329 tree->Branch("along", &ralong);
330
331 sector = 1;
332 string str;
333 for(Int_t i=1; i<nelem; i++)
334 {
335 //i= i + 1;
336 str.clear();
337 //getline(plik1,str);
338 //stringstream plik1(str);
339 //std::cout<<ouvh[i]<<endl;
340 //C
341 //C plik1.read (5,*) tm,typ,txt6,al,txt5
342 //C 100 format (A10)
343 //plik1.read(plik1), along[i], txt2, aa, sig11, txt5;
344
345
346 plik1>>along[i]>>txt2>>aa>>sig11>>txt5;
347 plik1>>aa>>sig22>>txt5>>r12;
348 plik1>>aa>>sig33>>txt5>>r13>>r23;
349 plik1>>aa>>sig44>>txt5>>r14>>r24>>r34;
350 plik1>>aa>>sig55>>txt5>>r15>>r25>>r35>>r45;
351 plik1>>aa>>sig66>>txt5>>r16>>r26>>r36>>r46>>r56;
352
353
354
355 getline(plik1,str,'\n');
356 getline(plik1,str,'\n');
357
358
359 plik1>>T11[i]>>T12[i]>>T13[i]>>T14[i]>>aa>>T16[i];
360
361 cout << T11[i] << endl;
362
363 plik1>>T21[i]>>T22[i]>>T23[i]>>T24[i]>>aa>>T26[i];
364 plik1>>T31[i]>>T32[i]>>T33[i]>>T34[i]>>aa>>T36[i]>>std::fixed;
365 // std::cout<<T31[i]<<endl;
366 plik1>>T41[i]>>T42[i]>>T43[i]>>T44[i]>>aa>>T46[i];
367 plik1>>aa>>aa>>aa>>aa>>aa>>aa;
368 plik1>>aa>>aa>>aa>>aa>>aa>>aa;
369 getline(plik1,str);
370
371
372 //C
373 //C lecture second ordre T1ij
374 //C
375 getline(plik1,str);
376 plik1>>ind1>>ind23>>T111[i];
377 plik1>>ind1>>ind23>>T112[i]>>ind1>>ind23>>T122[i];
378 plik1>>ind1>>ind23>>T113[i]>>ind1>>ind23>>T123[i]>>ind1>>ind23>>T133[i];
379 plik1>>ind1>>ind23>>T114[i]>>ind1>>ind23>>T124[i]>>ind1>>ind23>>T134[i]>>ind1>>ind23>>T144[i];
380 plik1>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa;
381 plik1>>ind1>>ind23>>T116[i]>>ind1>>ind23>>T126[i]>>ind1>>ind23>>T136[i]>>ind1>>ind23>>T146[i]>>ind1>>ind23>>aa>>ind1>>ind23>>T166[i];
382 getline(plik1,str);
383 plik2<<".."<<i<<".."<<T166[i]<<endl;
384
385 //C
386 //C lecture second ordre T2ij
387 //C
388 plik1>>ii>>ii>>T211[i];
389 plik1>>ii>>ii>>T212[i]>>ii>>ii>>T222[i];
390 plik1>>ii>>ii>>T213[i]>>ii>>ii>>T223[i]>>ii>>ii>>T233[i];
391 plik1>>ii>>ii>>T214[i]>>ii>>ii>>T224[i]>>ii>>ii>>T234[i]>>ii>>ii>>T244[i];
392 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
393 plik1>>ii>>ii>>T216[i]>>ii>>ii>>T226[i]>>ii>>ii>>T236[i]>>ii>>ii>>T246[i]>>ii>>ii>>aa>>ii>>ii>>T266[i];
394 getline(plik1,str);
395 plik2<<".."<<i<<".."<<T266[i]<<endl;
396 rT266 = T266[i];
397
398 //C
399 //C lecture second ordre T3ij
400 //C
401 plik1>>ind1>>ind23>>T311[i];
402 plik1>>ind1>>ind23>>T312[i]>>ind1>>ind23>>T322[i];
403 plik1>>ind1>>ind23>>T313[i]>>ind1>>ind23>>T323[i]>>ind1>>ind23>>T333[i];
404 plik1>>ind1>>ind23>>T314[i]>>ind1>>ind23>>T324[i]>>ind1>>ind23>>T334[i]>>ind1>>ind23>>T344[i];
405 plik1>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa;
406 plik1>>ind1>>ind23>>T316[i]>>ind1>>ind23>>T326[i]>>ind1>>ind23>>T336[i]>>ind1>>ind23>>T346[i]>>ind1>>ind23>>aa>>ind1>>ind23>>T366[i];
407 getline(plik1,str);
408 plik2<<".."<<i<<".."<<T366[i]<<endl;
409 rT366 = T366[i];
410
411 //C
412 //C lecture second ordre T4ij
413 //C
414
415 plik1>>ii>>ii>>T411[i];
416 plik1>>ii>>ii>>T412[i]>>ii>>ii>>T422[i];
417 plik1>>ii>>ii>>T413[i]>>ii>>ii>>T423[i]>>ii>>ii>>T433[i];
418 plik1>>ii>>ii>>T414[i]>>ii>>ii>>T424[i]>>ii>>ii>>T434[i]>>ii>>ii>>T444[i];
419 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
420 plik1>>ii>>ii>>T416[i]>>ii>>ii>>T426[i]>>ii>>ii>>T436[i]>>ii>>ii>>T446[i]>>ii>>ii>>aa>>ii>>ii>>T466[i];
421 getline(plik1,str);
422 plik2<<".."<<i<<".."<<T466[i]<<endl;
423 rT466 = T466[i];
424 //C
425 //C lecture second ordre T5ij
426 //C
427
428 plik1>>ii>>ii>>aa;
429 plik1>>ii>>ii>>aa>>ii>>ii>>aa;
430 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
431 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
432 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
433 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
434 getline(plik1,str);
435 plik2<<".."<<i<<".."<<aa<<endl;
436 raa = aa;
437
438 // cout<<"nowa petla"<<endl;
439 // getline(plik1,str);
440 // std::cout<<"pusto"<<asdf<<"???????????????"<<endl;
441
442 //C
443 //C fin de lecture des Tijk
444 //C
445 rapport = 0;
446 //C 10 cm :offset pour digitaliser la position
447 offset = 10;
448 //C sur les detecteurs
449 if (T126[i] != 0) {
450 rapport = -100 * T12[i] / T126[i];
451 }
452 if (rapport < - 15) {
453 rapport = -15;
454 }
455 if (rapport > 15) {
456 rapport = 15;
457 }
458 // rx = fem::ffloat(i) + 0.01f;
459 rx = float(i) + 0.01;
460 //C call hfill (40,rx,0.,rapport+20.)
461 //C write (6,*) i,rapport
462 plik2<<"..along ="<<along[i]<<"..T126 ="<<T126[i]<<endl;
463 ralong = along[i];
464 getline(plik1,str,'\n');
465 // getline(plik1,str);
466 rT166 = T166[i];
467 sector = i + 1;
468 //if (sector > 0)
469 tree->Fill();
470 }
471
472 plik1.close();
473
474 itotpaass = 0;
475 itotq9 = 0;
476 //C plan intermediaire
477 idet1 = 17;
478 //C idet2=20 ! entree Q5
479 //C idet2=22 ! entree dipole 2
480 //C idet2=24 ! entre Q7
481 //C point intermediaire entre Q7 et Q8
482 idet2 = 26;
483 //C idet2=27 ! entree Q8
484 //C idet2=28 ! sortie Q8
485 //C idet2=29 ! entree Q9
486 //C idet2=30 ! sortie Q9
487 //C idet2=31 ! 1 metre avant cible HADES
488 //C calcul des T1jk a une distance dist2 du point standard pour detecteur 2
489 //C distance a mettre en metre
490 dist1 = 0;
491 //C write (6,*) t11[idet2],t12[idet2],t13[idet2],t14[idet2],t16[idet2]
492 //std::setprecision(20);
493 // is it really necessary? dist1=0.0 !
494 /*
495 T11[idet1]=T11[idet1]+0.1*dist1*T21[idet1];
496 T12[idet1]=T12[idet1]+0.1*dist1*T22[idet1];
497 T13[idet1]=T13[idet1]+0.1*dist1*T23[idet1];
498 T14[idet1]=T14[idet1]+0.1*dist1*T24[idet1];
499 T16[idet1]=T16[idet1]+0.1*dist1*T26[idet1];
500 T111[idet1]=T111[idet1]+0.1*dist1*T211[idet1];
501 T112[idet1]=T112[idet1]+0.1*dist1*T212[idet1];
502 T113[idet1]=T113[idet1]+0.1*dist1*T213[idet1];
503 T114[idet1]=T114[idet1]+0.1*dist1*T214[idet1];
504 T116[idet1]=T116[idet1]+0.1*dist1*T216[idet1];
505 T122[idet1]=T122[idet1]+0.1*dist1*T222[idet1];
506 T123[idet1]=T123[idet1]+0.1*dist1*T223[idet1];
507 T124[idet1]=T124[idet1]+0.1*dist1*T224[idet1];
508 T126[idet1]=T126[idet1]+0.1*dist1*T226[idet1];
509 T133[idet1]=T133[idet1]+0.1*dist1*T233[idet1];
510 T134[idet1]=T134[idet1]+0.1*dist1*T234[idet1];
511 T136[idet1]=T136[idet1]+0.1*dist1*T236[idet1];
512 T144[idet1]=T144[idet1]+0.1*dist1*T244[idet1];
513 T146[idet1]=T146[idet1]+0.1*dist1*T246[idet1];
514 T166[idet1]=T166[idet1]+0.1*dist1*T266[idet1];
515 T31[idet1]=T31[idet1]+0.1*dist1*T41[idet1];
516 T32[idet1]=T32[idet1]+0.1*dist1*T42[idet1];
517 T33[idet1]=T33[idet1]+0.1*dist1*T43[idet1];
518 T34[idet1]=T34[idet1]+0.1*dist1*T44[idet1];
519 T36[idet1]=T36[idet1]+0.1*dist1*T46[idet1];
520 T311[idet1]=T311[idet1]+0.1*dist1*T411[idet1];
521 T312[idet1]=T312[idet1]+0.1*dist1*T412[idet1];
522 T313[idet1]=T313[idet1]+0.1*dist1*T413[idet1];
523 T314[idet1]=T314[idet1]+0.1*dist1*T414[idet1];
524 T316[idet1]=T316[idet1]+0.1*dist1*T416[idet1];
525 T322[idet1]=T322[idet1]+0.1*dist1*T422[idet1];
526 T323[idet1]=T323[idet1]+0.1*dist1*T423[idet1];
527 T324[idet1]=T324[idet1]+0.1*dist1*T424[idet1];
528 T326[idet1]=T326[idet1]+0.1*dist1*T426[idet1];
529 T333[idet1]=T333[idet1]+0.1*dist1*T433[idet1];
530 T334[idet1]=T334[idet1]+0.1*dist1*T434[idet1];
531 T336[idet1]=T336[idet1]+0.1*dist1*T436[idet1];
532 T344[idet1]=T344[idet1]+0.1*dist1*T444[idet1];
533 T346[idet1]=T346[idet1]+0.1*dist1*T446[idet1];
534 T366[idet1]=T366[idet1]+0.1*dist1*T466[idet1];
535 */
536 //----------------------------------------------
537 // dist2=-1.6; // ! distance a mettre en metre
538 // T11[idet2]=T11[idet2]+0.1*dist2*T21[idet2];
539 // T12[idet2]=T12[idet2]+0.1*dist2*T22[idet2];
540 // T13[idet2]=T13[idet2]+0.1*dist2*T23[idet2];
541 // T14[idet2]=T14[idet2]+0.1*dist2*T24[idet2];
542 // T16[idet2]=T16[idet2]+0.1*dist2*T26[idet2];
543 // T111[idet2]=T111[idet2]+0.1*dist2*T211[idet2];
544 // T112[idet2]=T112[idet2]+0.1*dist2*T212[idet2];
545 // T113[idet2]=T113[idet2]+0.1*dist2*T213[idet2];
546 // T114[idet2]=T114[idet2]+0.1*dist2*T214[idet2];
547 // T116[idet2]=T116[idet2]+0.1*dist2*T216[idet2];
548 // T122[idet2]=T122[idet2]+0.1*dist2*T222[idet2];
549 // T123[idet2]=T123[idet2]+0.1*dist2*T223[idet2];
550 // T124[idet2]=T124[idet2]+0.1*dist2*T224[idet2];
551 // T126[idet2]=T126[idet2]+0.1*dist2*T226[idet2];
552 // T133[idet2]=T133[idet2]+0.1*dist2*T233[idet2];
553 // T134[idet2]=T134[idet2]+0.1*dist2*T234[idet2];
554 // T136[idet2]=T136[idet2]+0.1*dist2*T236[idet2];
555 // T144[idet2]=T144[idet2]+0.1*dist2*T244[idet2];
556 // T146[idet2]=T146[idet2]+0.1*dist2*T246[idet2];
557 // T166[idet2]=T166[idet2]+0.1*dist2*T266[idet2];
558 // T31[idet2]=T31[idet2]+0.1*dist2*T41[idet2];
559 // T32[idet2]=T32[idet2]+0.1*dist2*T42[idet2];
560 // T33[idet2]=T33[idet2]+0.1*dist2*T43[idet2];
561 // T34[idet2]=T34[idet2]+0.1*dist2*T44[idet2];
562 // T36[idet2]=T36[idet2]+0.1*dist2*T46[idet2];
563 // T311[idet2]=T311[idet2]+0.1*dist2*T411[idet2];
564 // T312[idet2]=T312[idet2]+0.1*dist2*T412[idet2];
565 // T313[idet2]=T313[idet2]+0.1*dist2*T413[idet2];
566 // T314[idet2]=T314[idet2]+0.1*dist2*T414[idet2];
567 // T316[idet2]=T316[idet2]+0.1*dist2*T416[idet2];
568 // T322[idet2]=T322[idet2]+0.1*dist2*T422[idet2];
569 // T323[idet2]=T323[idet2]+0.1*dist2*T423[idet2];
570 // T324[idet2]=T324[idet2]+0.1*dist2*T424[idet2];
571 // T326[idet2]=T326[idet2]+0.1*dist2*T426[idet2];
572 // T333[idet2]=T333[idet2]+0.1*dist2*T433[idet2];
573 // T334[idet2]=T334[idet2]+0.1*dist2*T434[idet2];
574 // T336[idet2]=T336[idet2]+0.1*dist2*T436[idet2];
575 // T344[idet2]=T344[idet2]+0.1*dist2*T444[idet2];
576 // T346[idet2]=T346[idet2]+0.1*dist2*T446[idet2];
577 // T366[idet2]=T366[idet2]+0.1*dist2*T466[idet2];
578
579
580 plik2<<".."<<T11[idet2 - 1]<<".."<<T12[idet2 - 1]<<".."<<T13[idet2]<<".."<<T16[idet2 - 1]<<endl;
581 plik2<<".."<<T11[idet2]<<".."<<T12[idet2]<<".."<<T13[idet2]<<".."<<T14[idet2]<<".."<<T16[idet2]<<endl;;
582 plik2<<".."<<T31[idet2 - 1]<<".."<<T32[idet2 - 1]<<".."<<T33[idet2]<<".."<<T36[idet2 - 1]<<endl;
583 plik2<<".."<<T31[idet2]<<".."<<T32[idet2]<<".."<<T33[idet2]<<".."<<T34[idet2]<<".."<<T36[idet2]<<endl;
584
585
586
587 TFile *output1 = new TFile("piony1.root", "RECREATE");
588 TTree *treex = new TTree("Pion1","Pionbeam1");
589 treex->Branch("reste", &rreste);
590 treex->Branch("etape", &retape);
591 TH1F *random = new TH1F("random","random",3000,-1.5,1.5);
592 //C
593 //C zone utile detecteurs det1 et det2
594 ouvh[idet1] = 5.0;
595 //C zone utile detecteurs det1 et det2
596 ouvv[idet1] = 5.0;
597 //C zone utile detecteurs det1 et det2
598 ouvh[idet2] = 5.0;
599 //C zone utile detecteurs det1 et det2
600 ouvv[idet2] = 5.0;
601 //C
602 ouvh[33] = 4.6;
603 ouvv[33] = 4.6;
604 //C
605 ouvrescaleh = 1.;
606 ouvrescalev = 1.;
607 //float ide2,ii1,ii2,idet2,
608
609
610 r1 = float(idet2) + 0.001;
611 r1 = r1 / 10.;
612 ii1 = r1;
613 r1 = float(ii1) + 0.001;
614 ii2 = idet1 - 10 * ii1;
615 r2 = float(ii2) + 0.001;
616 r3 = float(idet2) + 0.001;
617 r3 = r3 / 10;
618 ii3 = r3;
619 r3 = float(ii3) + 0.001;
620 ii4 = idet2 - 10 * ii3;
621 r4 = float(ii4) + 0.001;
622
623 //std::cout<<r4<<endl;
624
625 //C
626 //C r1 a r4 quand le numero d'element ou sont mis les detecteurs
627 //C rh,rv,rres,rdm donne le facteur de reduction sur l'effet de
628 //C la taille h de l'objet
629 //C la taille v de l'objet
630 //C la resolution des detecteurs
631 //C la diffusion multiple
632 //C une valeur =0 indique que c'est pris en compte totalement
633 //C une valeur non nulle indique qu'on divise l'effet
634 //C correspondant par 10**rh, 10**rv, 10**rres ou 10**rdm
635 //C
636 //C on multiplie la taille hor de l'objet par 2
637 rh = -0.30103;
638 rh = 0.;
639 rv = 0.;
640 //C groupement des pistes par 4
641 rres = -0.60206;
642 //C groupement des pistes par 2
643 rres = -0.30103;
644 //C non groupement des pistes
645 rres = 0.;
646 //C rres=0.30103 ! multiplication du nombre de pistes par 2
647 //C multiplication du nombre de pistes par 2**3
648 rres = 0.30103 * 3.;
649 if (rres > 5.) {
650 rres = 5.;
651 }
652 rdm = 0.;
653
654 //C beam offset (in cm) at prod target in x
655 shiftx = 0.0;
656 //C beam offset (in cm) at prod target in y
657 shifty = 0.0;
658 dx = 0.05 / pow(10., rh);
659 //C change change change
660 dxp = 10.;
661 dy = 0.05 / pow(10., rv);
662 dyp = 50.;
663 //C doit etre une valeur entiere en %
664 dpspmax = 6.;
665 dpsppas = 1.0;
666
667
668 //ntir = 100000;
669 ntir = 1;
670
671
672
673
674 //C p en gev/c
675 ppion = 1.3;
676 //C ppion=0.7 ! p en gev/c
677 //C call hfill (41,9.001,0.,10.*ppion)
678 //C epaisseur 300 microns Si -> g.cm-2
679 esi = 300. / 10000. * rhosi;
680 //C =800 microns carre
681 resoldet = 800. / pow(10, rres) / 10000.;
682 //C
683 //C ********* fin des donnees definissant le probleme
684 //C
685 betapion = ppion / sqrt(pow(ppion,2) + pow(amasspi,2));
686 diffmul = 13.6 / ppion / betapion * sqrt(esi / x0si);
687 //Cdm raporte a det2
688 resoldmx = T12[idet2] * diffmul / T22[idet1];
689 resoldmx = resoldmx / pow(10., rdm);
690 //Cdm raporte a det2
691 resoldmy = T34[idet2] * diffmul / T44[idet1];
692 resoldmy = resoldmy / pow(10., rdm);
693 //write(plik2,"(' diff mult=',f5.2,' mrad resol equiv det2=',f5.2,' mm')"),
694 // write (6,777) diffmul,10.*resoldet
695 plik2<<"diff mul = "<<diffmul<<"resol="<<resoldet*10<<endl;
696 // diffmul, 10. * resoldet;
697 idineg = 0;
698 //C
699 //C ******** debut boucle sur les tirages
700 //C
701 xmax = 0.05;
702 ymax = 0.05;
703 xpmax = 10.;
704 ypmax = 50.;
705 //C utilise uniquement pour le calcul des contributions ma
706 dpmax = 5.;
707 //FEM_DO_SAFE(k, 1, nelem)
708
709 for(Int_t k = 1; k < nelem; k++){
710 // relem = fem::ffloat[k] + 0.1f;
711 relem = float(k) + 0.1;
712
713
714 z1 = 10. * fabs(T31[k] * xmax);
715 z2 = 10. * fabs(T32[k] * xpmax);
716 z3 = 10. * fabs(T33[k] * ymax);
717 z4 = 10. * fabs(T34[k] * ypmax);
718 z5 = 10. * fabs(T36[k] * dpmax);
719 z6 = 10. * fabs(T322[k] * xpmax * xpmax);
720 z7 = 10. * fabs(T323[k] * xpmax * ymax);
721 z8 = 10. * fabs(T324[k] * xpmax * ypmax);
722 z9 = 10. * fabs(T326[k] * xpmax * dpmax);
723 z10 = 10. * fabs(T333[k] * ymax * ymax);
724 z11 = 10. * fabs(T334[k] * ymax * ypmax);
725 z12 = 10. * fabs(T336[k] * ymax * dpmax);
726 z13 = 10. * fabs(T344[k] * ypmax * ypmax);
727 z14 = 10. * fabs(T346[k] * ypmax * dpmax);
728 z15 = 10. * fabs(T366[k] * dpmax * dpmax);
729
730
731
732 }
733
734 double rxf = 0,ryf = 0,rdp = 0,rxpf = 0,rypf = 0, relement = 0;
735
736 TFile *output_new = new TFile("piony_new.root", "RECREATE");
737 TTree *treec = new TTree("Pion5","Pionbeam5");
738 treec->Branch("xi", &rxf);
739 treec->Branch("yi", &ryf);
740 treec->Branch("dp", &rdp);
741 treec->Branch("xpi", &rxpf);
742 treec->Branch("ypi", &rypf);
743 treec->Branch("elemnt", &relement);
744
745
746 for(Int_t itir = 0; itir < ntir; itir++)
747 {
748 //C
749 //C tirage gaussien dans l'extension objet en x et y
750 //C
751 //C tirage uniforme en theta et phi
752 //C
753 //C dp par pas de 1% de +dpspmax a -dpspmax
754 //C
755 // yy = (2. * float(rand())/RAND_MAX - 1.);
756 // //cout<<yy<<endl;
757 // random->Fill(yy);
758 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
759 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
760 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
761 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
762 // yy = yy * renorgauss;
763 // xi = dx * yy + shiftx;
764 // xpi = dxp * (2. * float(rand())/RAND_MAX - 1.);
765 // yy = 2. * float(rand())/RAND_MAX - 1.;
766 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
767 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
768 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
769 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
770 // yy = yy * renorgauss;
771 // yi = dy * yy + shifty;
772 // ypi = dyp * (2. * float(rand())/RAND_MAX - 1.);
773
774 xi = 0;
775 yi = 0;
776 xpi = 0;
777 ypi = 0.05;
778
779
780
781 //C
782 //C on genere de facon equiprobable des delta(p)/p de dpspmax
783 //C
784 //C entre -dpspmax et +dpspmax (bornes comprises)
785 //C
786 //C par pas de dpsppas
787 //C
788 // dp = dpspmax * 2. * float(rand())/RAND_MAX - 1.;
789 // dd = float(rand())/RAND_MAX;
790 // ndp = 2 * int((dpspmax + 0.01) / dpsppas) + 1;
791 // dp = float((ndp - 1) / 2) * dpsppas;
792 // uns = 1. / float(ndp);
793 // dp = 1;
794
795 dp = -0.08;
796 // dd = 1;
797
798 // //FEM_DO_SAFE(kii, 1, ndp)
799 // for(Int_t kii = 1; kii < ndp; kii++){
800 // seuil = uns * float(ndp - kii);
801 // if (dd < (seuil)) {
802 // dp = dp - dpsppas;
803 // }
804 // }
805
806 //C xi=1.
807 //C xpi=0.yf=T31(k)*xi+T32(k)*xpi+T33(k)*yi+T34(k)*ypi+T36(k)*dp
808 //C yi=1.
809 //C ypi=20.
810 //C dp=-6.
811 //C
812 //C ******** debut boucle sur les elements magnetiques
813 //C
814 iok = 1;
815 //FEM_DO_SAFE(k, 1, nelem)
816 for(Int_t k = 1; k < nelem; k++){ //xf is OK!!!!!!!!!!!
817 xf =0;
818 // yf =0;
819 //std::setprecision(20);
820 //std::cin.scientific
821 z1=T11[k]*xi;
822 z2=T12[k]*xpi;
823 z3=T13[k]*yi;
824 z4=T14[k]*ypi;
825 z5=T16[k]*dp;
826 // cout << xf << endl;
827
828 //cout << z1 << ":" << z2 << ":" << z3 << ":" << z4 << ":" << z5 << endl;
829
830 // cout<<z5<<endl;
831 xf=z1+z2+z3+z4+z5;
832 yf=T31[k]*xi+T32[k]*xpi+T33[k]*yi+T34[k]*ypi+T36[k]*dp; //yf=T31(k)*xi+T32(k)*xpi+T33(k)*yi+T34(k)*ypi+T36(k)*dp
833 // std::cout<<"T31 = .. "<<T31[k]<<"..T32 = .. "<<T32[k]<<"..T33 = .. "<<T33[k]<<"..T34 = .. "<<T34[k]<<"..T36 = .. "<<T36[k]<<std::fixed<<endl;
834 xpf=T21[k]*xi+T22[k]*xpi+T23[k]*yi+T24[k]*ypi+T26[k]*dp;
835 ypf=T41[k]*xi+T42[k]*xpi+T43[k]*yi+T44[k]*ypi+T46[k]*dp;
836 xf1=xf;
837 yf1=yf;
838 xpf1=xpf;
839 ypf1=ypf;
840
841 //cout << xf << endl;
842
843 xf=xf1+T111[k]*xi*xi+T112[k]*xi*xpi+T113[k]*xi*yi;
844 xf=xf+T114[k]*xi*ypi+T116[k]*xi*dp;
845 xf=xf+T122[k]*xpi*xpi+T123[k]*xpi*yi+T124[k]*xpi*ypi;
846 xf=xf+T126[k]*xpi*dp;
847 xf=xf+T133[k]*yi*yi+T134[k]*yi*ypi+T136[k]*yi*dp;
848 xf=xf+T144[k]*ypi*ypi+T146[k]*ypi*dp;
849 xf=xf+T166[k]*dp*dp;
850
851 //cout << xf << endl;
852
853 // xf=xf1+T126[k]*xpi*dp+T166[k]*dp*dp;
854 // xf=xf1+T126[k]*xpi*dp;
855
856 xpf=xpf1+T211[k]*xi*xi+T212[k]*xi*xpi+T213[k]*xi*yi;
857 xpf=xpf+T214[k]*xi*ypi+T216[k]*xi*dp;
858 xpf=xpf+T222[k]*xpi*xpi+T223[k]*xpi*yi+T224[k]*xpi*ypi;
859 xpf=xpf+T226[k]*xpi*dp;
860 xpf=xpf+T233[k]*yi*yi+T234[k]*yi*ypi+T236[k]*yi*dp;
861 xpf=xpf+T244[k]*ypi*ypi+T246[k]*ypi*dp;
862 xpf=xpf+T266[k]*dp*dp;
863
864 yf=yf1+T311[k]*xi*xi+T312[k]*xi*xpi+T313[k]*xi*yi;
865 yf=yf+T314[k]*xi*ypi+T316[k]*xi*dp;
866 yf=yf+T322[k]*xpi*xpi+T323[k]*xpi*yi+T324[k]*xpi*ypi;
867 yf=yf+T326[k]*xpi*dp;
868 yf=yf+T333[k]*yi*yi+T334[k]*yi*ypi+T336[k]*yi*dp;
869 yf=yf+T344[k]*ypi*ypi+T346[k]*ypi*dp;
870 yf=yf+T366[k]*dp*dp;
871
872 ypf=ypf1+T411[k]*xi*xi+T412[k]*xi*xpi+T413[k]*xi*yi;
873 ypf=ypf+T414[k]*xi*ypi+T416[k]*xi*dp;
874 ypf=ypf+T422[k]*xpi*xpi+T423[k]*xpi*yi+T424[k]*xpi*ypi;
875 ypf=ypf+T426[k]*xpi*dp;
876 ypf=ypf+T433[k]*yi*yi+T434[k]*yi*ypi+T436[k]*yi*dp;
877 ypf=ypf+T444[k]*ypi*ypi+T446[k]*ypi*dp;
878 ypf=ypf+T466[k]*dp*dp;
879 //C
880 //C yf=yf1+T346[k]*ypi*dp+T366[k]*dp*dp
881 //C yf=yf1+T346[k]*ypi*dp
882 // rxf =
883 xf2 = xf;
884 xpf2 = xpf;
885 yf2 = yf;
886 ypf2 = ypf;
887
888 //xf is ok
889 // std::cout<<xf<<endl;
890 if (fabs(xf) > (ouvrescaleh * ouvh[k])) { iok = 0;}
891 if (fabs(yf) > (ouvrescalev * ouvv[k])) { iok = 0;}
892
893 if (iok > 0) {
894 ipass[k]++;
895 }
896 deltax = 0;
897 deltay = 0;
898 if (xf1 != 0) {
899 deltax = fabs((xf2 - xf1) / xf1 * 100);
900 }
901 if (yf1 != 0) {
902 deltay = fabs((yf2 - yf1) / yf1 * 100);
903 }
904 xhad = xf2;
905 yhad = yf2;
906
907 if ((k == 30) && (iok == 1)) {
908 itotq9++;
909
910 xf = xf1 + T111[33] * xi * xi + T112[33] * xi * xpi + T113[33] * xi * yi;
911 xf = xf + T114[33] * xi * ypi + T116[33] * xi * dp;
912 xf = xf + T122[33] * xpi * xpi + T123[33] * xpi * yi + T124[33] * xpi * ypi;
913 xf = xf + T126[33] * xpi * dp;
914 xf = xf + T133[33] * yi * yi + T134[33] * yi * ypi + T136[33] * yi * dp;
915 xf = xf +T144[33] * ypi * ypi + T146[33] * ypi * dp;
916 xf = xf +T166[33] * dp * dp;
917 yf = yf1 + T311[33] * xi * xi + T312[33] * xi * xpi + T313[33] * xi * yi;
918 yf = yf + T314[33] * xi * ypi + T316[33] * xi * dp;
919 yf = yf + T322[33] * xpi * xpi + T323[33] * xpi * yi + T324[33] * xpi * ypi;
920 yf = yf + T326[33] * xpi * dp;
921 yf = yf + T333[33] * yi * yi + T334[33] * yi * ypi + T336[33] * yi * dp;
922 yf = yf + T344[33] * ypi * ypi + T346[33] * ypi * dp;
923 yf = yf + T366[33] * dp * dp;
924 radius = sqrt(pow(xf,2) + pow(yf,2));
925 //C call hfill (42,radius,0.,1./282450.)
926 //C call hfill (42,radius,0.,1./240539.)
927 }
928 rxf = xhad;
929 ryf =yhad;
930 rdp = dp;
931 rxpf = xpf2;
932 rypf = ypf2;
933 relement = k + 1;
934 treec ->Fill();
935
936 cout << xf2 << ":" << yf2 << ":" << xpf2 << ":" << ypf2 << endl;
937
938
939 }
940 //C
941 //C ******** fin boucle sur les elements magnetiques
942 //C
943 ijk = (dp + 10.01);
944
945 if (iok > 0) {
946 kpass[ijk]++;
947 }
948 // if (iok == 1) {
949 // itotpass++;
950
951 // //C call hfill (36,ydet1,ydet2,1.)
952 // //C
953 // //C discretisation du detecteur 1
954 // xdet1sav = xdet1;
955 // xca = (xdet1 + 0.5 * resoldet + offset) / resoldet;
956 // nca = xca;
957 // xdet1 = float(nca) * resoldet - offset;
958 // //C
959 // //C discretisation du detecteur 2
960 // xca = (xdet2 + 0.5 * resoldet + offset) / resoldet;
961 // nca = xca;
962 // xdet2 = float(nca) * resoldet - offset;
963 // xdet2sav = xdet2;
964 // //C
965 // //C diffusion multiple
966 // yy = (2. * float(rand())/RAND_MAX - 1.);
967 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
968 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
969 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
970 // yy = yy + (2. * float(rand())/RAND_MAX - 1.);
971 // yy = yy * renorgauss;
972 // //C diff mult gaussien
973 // xdet2 = xdet2 + resoldmx * yy;
974 // //C
975 // //C discretisation du detecteur 1 en y
976 // //C
977 // ydet1sav = ydet1;
978 // ydet2sav = ydet2;
979
980 // alfa = -(T16[idet1] * T126[idet2] - T16[idet2] * T126[idet1]);
981 // beta = T16[idet2] * T12[idet1] - xdet2 * T126[idet1];
982 // beta = beta - T16[idet1] * T12[idet2] + xdet1 * T126[idet2];
983 // gamma = xdet1 * T12[idet2] - xdet2 * T12[idet1];
984 // discri = pow(beta,2) - 4. * alfa * gamma;
985 // galfa = T126[idet1] * T166[idet2] - T126[idet2] * T166[idet1];
986 // gbeta = T12[idet1] * T166[idet2] - T12[idet2] * T166[idet1];
987 // gbeta = gbeta - T16[idet1] * T126[idet2] + T16[idet2] * T126[idet1];
988 // ggamma = T12[idet1] * T16[idet2] - T12[idet2] * T16[idet1];
989 // ggamma = ggamma + xdet1 * T126[idet2] - xdet2 * T126[idet1];
990 // gdelta = xdet1 * T12[idet2] - xdet2 * T12[idet1];
991 // aa1 = 1.;
992 // aa2 = gbeta / galfa;
993 // aa3 = ggamma / galfa;
994 // aa4 = gdelta / galfa;
995 // sol = 3.;
996 // iter = 0;
997 // //C write (6,755) iter,sol,dp
998 // yyy1 = galfa * pow(sol,3) + gbeta * pow(sol,2) +
999 // ggamma * sol + gdelta;
1000 // yyp1 = 3. * galfa * pow(sol,2) + 2. * gbeta * sol + ggamma;
1001 // sol = sol - yyy1 / yyp1;
1002 // iter = iter + 1;
1003 // //C write (6,755) iter,sol,dp
1004 // yyy2 = galfa * pow(sol,3) + gbeta * pow(sol,2) +
1005 // ggamma * sol + gdelta;
1006 // yyp2 = 3. * galfa * pow(sol,2) + 2. * gbeta * sol + ggamma;
1007 // sol = sol - yyy2 / yyp2;
1008 // iter = iter + 1;
1009 // //C write (6,755) iter,sol,dp
1010 // yyy3 = galfa * pow(sol,3) + gbeta * pow(sol,2) +
1011 // ggamma * sol + gdelta;
1012 // yyp3 = 3. * galfa * pow(sol,2) + 2. * gbeta * sol + ggamma;
1013 // sol = sol - yyy3 / yyp3;
1014 // iter = iter + 1;
1015 // //C write (6,755) iter,sol,dp
1016 // yyy4 = galfa * pow(sol,3) + gbeta * pow(sol,2) +
1017 // ggamma * sol + gdelta;
1018 // yyp4 = 3. * galfa * pow(sol,2) + 2. * gbeta * sol + ggamma;
1019 // sol = sol - yyy4 / yyp4;
1020 // iter = iter + 1;
1021 // //C write (6,755) iter,sol,dp
1022 // yyy5 = galfa * pow(sol,3) + gbeta * pow(sol,2) +
1023 // ggamma * sol + gdelta;
1024 // yyp5 = 3. * galfa * pow(sol,2) + 2. * gbeta * sol + ggamma;
1025 // sol = sol - yyy5 / yyp5;
1026 // iter = iter + 1;
1027 // //C write (6,755) iter,sol,dp
1028 // teta = (xdet1 - T16[idet1] * sol - T166[idet1] * sol * sol);
1029 // teta = teta / (T12[idet1] + T126[idet1] * sol);
1030
1031 // // goto statment_999;
1032 // aa = T33[idet1];
1033 // bb = T34[idet1] + T346[idet1] * sol + T324[idet1] * teta;
1034 // cc = ydet1 - T36[idet1] * sol;
1035 // cc = cc - T326[idet1] * teta * sol - T366[idet1] * sol * sol;
1036 // aap = T33[idet2];
1037 // bbp = T34[idet2] + T346[idet2] * sol + T324[idet2] * teta;
1038 // ccp = ydet2 - T36[idet2] * sol;
1039 // ccp = ccp - T326[idet2] * teta * sol - T366[idet2] * sol * sol;
1040 // deter = aa * bbp - aap * bb;
1041 // yci = (cc * bbp - bb * ccp) / deter;
1042 // phi = (aa * ccp - cc * aap) / deter;
1043 // // statment_999:
1044
1045 // poids = 100. / float(ntir) * float(ndp);
1046 // poids = 100. / 102962.;
1047
1048 // xf = xf1 + T111[33] * xi * xi + T112[33] * xi * xpi + T113[33] * xi * yi;
1049 // xf = xf + T114[33] * xi * ypi + T116[33] * xi * dp;
1050 // xf = xf + T122[33] * xpi * xpi + T123[33] * xpi * yi + T124[33] * xpi * ypi;
1051 // xf = xf + T126[33] * xpi * dp;
1052 // xf = xf + T133[33] * yi * yi + T134[33] * yi * ypi + T136[33] * yi * dp;
1053 // xf = xf + T144[33] * ypi * ypi + T146[33] * ypi * dp;
1054 // xfth = xf + T166[33] * dp * dp;
1055 // yf = yf1 + T311[33] * xi * xi + T312[33] * xi * xpi + T313[33] * xi * yi;
1056 // yf = yf + T314[33] * xi * ypi + T316[33] * xi * dp;
1057 // yf = yf + T322[33] * xpi * xpi + T323[33] * xpi * yi + T324[33] * xpi * ypi;
1058 // yf = yf + T326[33] * xpi * dp;
1059 // yf = yf + T333[33] * yi * yi + T334[33] * yi * ypi + T336[33] * yi * dp;
1060 // yf = yf + T344[33] * ypi * ypi + T346[33] * ypi * dp;
1061 // yfth = yf + T366[33] * dp * dp;
1062 // //C
1063 // xf = xf1 + T111[33] * 0. * 0. + T112[33] * 0. * teta + T113[33] * 0. * yci;
1064 // xf = xf + T114[33] * 0. * phi + T116[33] * 0. * sol;
1065 // xf = xf + T122[33] * teta * teta + T123[33] * teta * yci + T124[33] * teta * phi;
1066 // xf = xf + T126[33] * teta * sol;
1067 // xf = xf + T133[33] * yci * yci + T134[33] * yci * phi + T136[33] * yci * sol;
1068 // xf = xf + T144[33] * phi * phi + T146[33] * phi * sol;
1069 // xfcal = xf + T166[33] * sol * sol;
1070 // yf = yf1 + T311[33] * 0. * 0. + T312[33] * 0. * teta + T313[33] * 0. * yci;
1071 // yf = yf + T314[33] * 0. * phi + T316[33] * 0. * sol;
1072 // yf = yf + T322[33] * teta * teta + T323[33] * teta * yci + T324[33] * teta * phi;
1073 // yf = yf + T326[33] * teta * sol;
1074 // yf = yf + T333[33] * yci * yci + T334[33] * yci * phi + T336[33] * yci * sol;
1075 // yf = yf + T344[33] * phi * phi + T346[33] * phi * sol;
1076 // yfcal = yf + T366[33] * sol * sol;
1077 // rr = sqrt(pow((xfth - xfcal),2) + pow((yfth - yfcal),2)); // matlab.hpp
1078 // ipi = dp + 6.001;
1079 // ips2 = ipi / 1;
1080 // idif = ipi - ips2 * 1;
1081 // //C if (idif.eq.0) call hfill (230+ips2,10.*yfth,0.,1.)
1082 // xabs = 10. * xfth;
1083 // yabs = 10. * (xfcal - xfth) + dp * 10.;
1084
1085 // xabs = 10. * yfth;
1086 // yabs = 10. * (yfcal - yfth) + dp * 5.;
1087
1088 // }
1089
1090 }
1091 random->Write();
1092 double rrpass = 0;
1093 Int_t dpsp = 0;
1094
1095 TFile *output2 = new TFile("piony2.root", "RECREATE");
1096 TTree *treez = new TTree("Pion2","Pionbeam2");
1097 treez->Branch("alossfraction", &rreste);
1098 treez->Branch("dpsp", &retape);
1099 // output1->Write();
1101
1102
1103 for(Int_t ki=1; ki<nelem; ki++){
1104 reste = 100. * float(ipass[ki]) / float(ntir);
1105 plik2<<"etape.."<<ki<<"..reste.."<<reste<<endl;
1106 ri = float(ki) + 0.01;
1107 rreste = reste;
1108 retape = ki+1;
1109 if (retape > 0) treex->Fill();
1110
1111
1112 }
1113
1114 for(Int_t kij = 1; kij<16; kij++){
1115 rpass[kij] = float(kpass[kij]) * 15. / float(ntir) * 100.;
1116 plik2<<"dpsp.."<<kij - 7<<"..transsmision.."<<rpass[kij]<<endl;
1117 rrpass = rpass[kij];
1118 dpsp = (kij - 7)+1;
1119 if (dpsp > 0) treez->Fill();
1120 }
1121 plik2<<"idineg...."<<idineg<<endl;
1122 alossfraction = float(itotq9 - itotpass) / float(itotq9);
1123 plik2<<"fraction hors cible.."<<alossfraction * 100.;
1124
1125
1126 plik2.close();
1127 output->Write();
1128 output1->Write();
1129 output2->Write();
1130 output_new->Write();
1131}
cc
PChannel * cc[]
Definition dp_2delta_dalitz.C:36
piony
int piony()
Definition piony.C:42