/
molmc10.f
4031 lines (3723 loc) · 122 KB
/
molmc10.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c-------------------------------------------------------------------------------
c Original Monte Carlo code of Swartz originally written for the
c SLAC Moller polarimeter (NEWMOL); the corresponding paper:
c M. Swartz, NIM Nucl. Instr. and Meth. A 363 (1995) 526
c
c This is a modified version of NEWMOL for a two arm polarimeter
c Include real momentum distributions for Fe+Co electrons
c This version is modified for use with the CEBAF Moller Polarimeter
c af. 10.08.94.
c
c All subroutines for the CEBAF optics programmed by Laurens de Bever (L.B.)
c
ccc
ccc Modifications to get the information available in ntuples for
ccc reanalyzing with PAW, by Matthias Loppacher 12.94.
ccc Comments and additional modifications, by Matthias Loppacher 11.96. (M.L.)
ccc !Histo specific for histo analysis
ccc !Ntuple specific for ntuple analysis
ccc
c
c-------------------------------------------------------------------------------
C
PROGRAM TWOMOL
c
IMPLICIT NONE
C
character*128 argin(10)
integer nrarg, iargc, inunit
real*8 quad1x, quad1y, quad3x, quad3y !offset of quads
integer npads !can be used for #hodoscope channels (set also in
parameter (npads=1) !routine Detectors)
REAL*8 Lsol, LQ1, LQ2, LQ3, GdotL1,GdotL2, GdotL3
REAL*8 ZDIST,VEC(6,2)
real*8 quad1exit(2,2)
real*8 quad2entrance(2,2),quad2exit(2,2)
real*8 quad3entrance(2,2),quad3exit(2,2),q3mechexit(2,2)
real*8 pipeexit(2,2),armentr(2,2),armexit(2,2)
REAL*8 PEL,PBEAM,CST,TGEN,PHIGEN,WEIGHT,PWT(2,3)
real*4 p_up
real*8 pmoll(2) !lab momentum of electrons
real*8 mytheta,myphi
REAL*8 XTGT,DEPTH,XDEPTH,PB
REAL*8 WGTBM,WGTML
real*8 WGT,ASX2,ASY2,ASZ2,wgtcx
REAL*8 SIGNAL,SIGNL2,ASZ,ASX,ASY,DIFF,WEITOT,WEITO2
REAL*8 siSIGNAL, siSIGNL2, DBsignal, DBsignl2
real*8 WEITOTs(2),WEITO2s(2)
REAL*8 BLKTL(2,3),BLKTL2(2,3)
real*8 avgA,avgA2,davgA
real*8 weights(2)
real*8 TC, TM1, M1M2, M2M3,M1C, CM2, M3D(2)
real*8 BSolLong, MinSolFld, SM1, soldx, soldy, soldxp, soldyp
parameter (MinSolFld=0.1d0)
real*4 Azz(npads,2), Azze(npads,2)
real*4 DbAzz, DbAzze, Dbsig, Dbsige
INTEGER I,J,K,L,NUMEVENTS,ITRIAL,IHD
INTEGER solsim
integer nplt, nplt2, nplt3, nplt4
integer ienergy
real*8 padtl(2,3,npads,2), padtl2(2,3,npads,2)
INTEGER II, JJ
LOGICAL TCKBRM,TAILS,BINDE,Q1CLIP,QFRINGE,THREEQUAD
logical twomols, singlemol(2), scflag(3)
CHARACTER*12 TITLE(2)
CHARACTER*55 HEADING
CHARACTER*80 FILENAME
real*8 b_energy, mcur1, mcur2, r1, r2, r3, tip1, tip2, tip3
real*8 grad1, grad2, grad3
real*8 q1tip, q2tip, q3tip
real*8 q1rad, rq1_1, rq1_2 !Q1 clearance, and electron 1 & 2 radii at Q1 exit
data q1rad / 39.0/ !Q1 exit clearance (radius, in mm)
common /collimators/ pcol,colxmin,colxmax,colymin,colymax
real*8 pcol(7)
real*8 colxmin(7),colxmax(7),colymin(7),colymax(7)
real*8 col(2,7,2)
real*8 col1z, col3z, col5z, col6z
real*8 PPOL, PUNPOL, TGTPOL
real*8 targetpol, twom, efficiency
real*8 XQ1, XQ2, XQ3,XDET(2)
real*8 ran3
real*8 quad2offset, quad3offset
real*8 polfact, polfacte
real*8 dgpolfact,dgpolfacte
logical all_coll !function declaration
ccc
INTEGER NWPAWC, H, iquest
PARAMETER(NWPAWC= 5000000) ! OK for 1'100'000 events (6 effects on)
COMMON /PAWC/ H(NWPAWC)
c common/quest/iquest(100) ! default 16000
common/quest/iquest(16000) ! default 16000
real*8 climit, crange
common /anglerange/ climit,crange
c
integer*4 istat,icycle ! M.L. L.B. 94
c real*4 ntuple(20)
c common/cern/ ntuple
ccc
COMMON / RANSEED / II
REAL*8 SIGX,SIGXP,SIGY,SIGYP,X0,Y0,TX0,TY0
common /beamparam/ SIGX,SIGY,SIGXP,SIGYP,X0,Y0,TX0,TY0
character*1 answer
C
c gradient parameters for quadrupoles
c data Q1/-0.00438,-0.00578,-0.00608,-0.00572,-0.00502,-0.00408/
c data Q2/0.001229,0.002621,0.004308,0.006221,0.008336,0.010653/
C
C TCKBRM and TAILS switch on thick target bremsstrahlung and
C multiple scattering with Moliere tails, respectively
C BINDE switches on the electron binding effects
C
DATA TITLE/' 10-UM MOLIE',' 10-UM GAUSS'/
C
CALL HLIMIT(NWPAWC)
iquest(10)=64000 ! 16000 -> 64000 events times 4
! has to be placed after hlimit!
solsim = 0
twom = 0.
C
C this may need to be updated depending on field
targetpol = 0.08043
c targetpol = 1.0
C quad magnet aperture (radius)
c r1=47.625
r1=47.6
r2=127.
r3=127. !3rd quad, same as 2nd
C
C MAGNET LENGTH
C Lsol= 125. !This gets set later now
LQ1 = 361.92 !250.
c LQ2 = 965.2 ! this is 914.4+0.6*r2
LQ2 = 990.6 ! this is 39" eff. length
LQ3 = 990.6
C
C THE POLARIMETER GEOMETRY: TC = TARGET to COLLIMATOR 1 DISTANCE (MM)
C TM1 = TARGET TO 1st quad MAGNET ENTRANCE
C SM1 = Solenoid exit to 1st Quad Entrance
C M1M2= 1st quad MAGNET exit to 2nd quad entrance
C M2M3= 2nd quad MAGNET exit to 3rd quad entrance
c M3D = 3rd quad MAG exit to detectors
c M1C = 1st quad MAG exit to collimator
c CM2 = collimator 6 to 2nd quad MAG entrance
C DJG Original distances before 2010
c col1z = 2000.4 ! collimator 1 distance form target
c col3z = 2090.9 ! collimator 3
c col5z = 2203.6 ! collimator 5
c col6z = 2295.7 ! collimator 6
c TC = 2000.4
C DJG New distances for Qweak and after
C DJG Q2 has been shifted upstream 123.604 mm, so I assume the collimator box moves with it.
C DJG An additional 30.48 cm was needed to make more space for the Compton, so I took it up
C DJG between Q1 and Q2. Not ideal, but works ok.
c XDET = (11.592-0.3048)*1.0E3
c XDET = (11.592-0.3048-0.09)*1.0E3 ! rough number from Kelly T.
c XDET = (11190.45 - 0.06) !~11m is average targ to det. distance,
c !measured 7/10 JAM
XDET(1) = 11191.77 !left det survey Fall 2011
XDET(2) = 11188.23 !right det survey Fall 2011
c col1z = 2000.4-123.604-304.8
c col3z = 2090.9-123.604-304.8
c col5z = 2203.6-123.604-304.8
c col6z = 2295.7-123.604-304.8
c TC = 2000.4-123.604-304.8
C DJG: updated for survey results from July 2010
C Collimator box is fiducialized relative to center of collimator 5
C Ideal distance from target is 1757.09 mm, with 4.32 mm offset downstream
col5z = 1757.09+4.32
col3z = col5z-112.7
col1z = col5z-203.2
col6z = col5z+92.1 !note, this value of 92.1 mm is apparently 4.3 mm larger than survey
TC=col1z
col1z = 2000.4-123.604-304.8
col3z = 2090.9-123.604-304.8
col5z = 2203.6-123.604-304.8
col6z = 2295.7-123.604-304.8
TC = 2000.4-123.604-304.8
C DJG GEP/SANE position
c XQ1 = 1.0008e3
C 2010 position
c XQ1 = 0.8484e3
XQ1 = 848.4 - 0.06 + 2.6 !1024.6 is nominal value from dimad. .06 is offset from targ, 2.6 from Q1 JAM
C DJG GEP/SANE position
c XQ2 = 3.1195e3
C 2010 position
c XQ2 = 3.0016196e3-304.8
XQ2 = 2696.75 - 0.06 -2.84 !z-offset for Q2 is -2.84 JAM
c XQ3 = 4.334612e3-304.8
XQ3 = 4029.83 - 0.06 - 0.12 !z-offset for Q3 is -0.12 JAM
TM1 = XQ1 - LQ1/2.
C SM1 = TM1 - LSol
M1M2 = XQ2 - LQ2/2. - XQ1 - LQ1/2.
M2M3 = XQ3 - LQ3/2. - XQ2 - LQ2/2.
c M3D = XDET - XQ3 - LQ3/2.
M3D(1) = XDET(1) - XQ3 - LQ3/2.
M3D(2) = XDET(2) - XQ3 - LQ3/2.
M1C = TC - Tm1 - LQ1
CM2 = XQ2 - LQ2/2. - col6z
c-- Here we read in the parameters for this run..
c--
c-----------------
c-- If there is an input file specified as a command argument, use it.
c-- If there is no command line argumenet, try to use file mc_input.dat
c-- If that file does not exist, read input from stdin.
nrarg = iargc()
if (nrarg .gt. 0) then
inunit= 31
do j=1,nrarg
call getarg(j,argin(j))
enddo
write(filename,'("infiles/",(a))') argin(1)
i=index(filename,' ')
write(filename(i:),'(".inp")')
write(6,*) 'filename is', filename
open(31,file=filename,status='old')
else
inunit= 32
open(unit=32,file='infiles/mc_input.dat',status='old',err=1098)
endif
goto 1099
1098 inunit=5 !if all else failed, try to read from stdin
1099 continue
c-----------------
C RANDOM NUMBER GENERATOR SEED
if(inunit.eq.5) write(6,'(a)')' Random number generator seed: '
read(inunit,*) ii
print*, 'Random number generator seed: ',ii
C THE BEAM ENERGY
C
if(inunit.eq.5) write(6,'(a)')' BEAM Energy: '
read(inunit,*) b_energy
PBEAM= b_energy !beam momentem in GeV/c (E=p for e-'s)
p_up= int(pbeam/2.0)+1.0 !nice upper limit for moller momenta histograms
C
C THE QUADRUPOLE MAGNET CURRENTS
C
if(inunit.eq.5) write(6,'(a)')' QUAD currents q1 & q2: '
read(inunit,*)mcur1, mcur2
c tip1 = q1tip(mcur1) !pole tip fields
c tip2 = q2tip(mcur2)
c tip3 = q2tip(mcur2)
C DJG Let MC handle negative currents
C DJG note that these sign conventions are the opposite
C DJG what is usually done in TRANSPORT
C DJG Normally, horizontally focusing (Q1) = positive field
C DJG vertically focusing (Q2) = negative field
if(mcur1.lt.0.0) then
tip1 = q1tip(abs(mcur1)) !pole tip fields
else
tip1 = -q1tip(mcur1)
endif
if(mcur2.lt.0.0) then
tip2 = -q2tip(abs(mcur2))
else
tip2 = q2tip(mcur2)
endif
tip3 = tip2
c grad1=-tip1/r1 ! field gradients
grad1=tip1/r1 ! field gradients (already built into sign of field)
grad2=tip2/r2
grad3=tip3/r3
if(inunit.eq.5) print *,tip1,tip2,mcur1,mcur2,grad1,grad2,grad3
C
if(inunit.eq.5) write(6,'(a)')' Number of events: '
read(inunit,*)NUMEVENTS
if(inunit.eq.5) write(6,'(a)')' scatt. Angle range (from,howmuch) [rad]:'
read(inunit,*) climit,crange
if(inunit.eq.5) write(6,'(a)')' BEAM OFFSET (X/Y) and ANGLE (X/Y): '
read(inunit,*) X0,Y0,TX0,TY0
C SIGMAS FOR THE BEAM SIZE AND DIVERGENCE AT THE M0LLER TARGET
C
C DATA SIGX,SIGY,SIGXP,SIGYP/0.500,0.200,0.48D-6,0.55D-6/
if(inunit.eq.5) then
write(6,'(a)')' Beam size(X/Y) and divergence(X/Y): '
endif
read(inunit,*) SIGX,SIGY,SIGXP,SIGYP
if(inunit.eq.5) then
write(6,'(a)')' Radiative Corrections ?(y/n): '
endif
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
scflag(2) = .true.
TCKBRM = .true.
else
scflag(2) = .false.
TCKBRM = .false.
endif
if(inunit.eq.5) write(6,'(a)')' Multiple Scattering ?(y/n): '
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
TAILS = .true. !Multiple scattering with moliere radius
scflag(1) = .true. !multiple satt. before interation
scflag(3) = .true. !multiple satt. after interation
else
TAILS = .false.
scflag(1) = .false.
scflag(3) = .false.
endif
if(inunit.eq.5) write(6,'(a)')' Levchuk Effect ?(y/n): '
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
BINDE = .TRUE.
else
BINDE = .false.
endif
C THE TARGET THICKNESS (MM)
if(inunit.eq.5) write(6,'(a)')' Target Thickness [mm]: '
read(inunit,*) xtgt
* XTGT = 0.010
ccc XTGT = 0.020 ! 20 UM LONGITUDINAL
if(inunit.eq.5) write(6,'(a)')' Q2 position offset [mm]: '
read(inunit,*) quad2offset
quad3offset=quad2offset
C THE COLLIMATOR POSITIONS
if(inunit.eq.5) then
write(6,'(a)') 'Settings of Collimators 1,2,3,4,6,7 (cm) //
>[all .lt. 0.0]:'
endif
read(inunit,*) pcol(1),pcol(2),pcol(3),pcol(4),pcol(6),pcol(7)
pcol(5)= 0.0 !fixed -- centered on beamline and irrelevant.
do i=1,7
colxmin(i)= -999.9
colxmax(i)= 999.9
colymin(i)= -999.9
colymax(i)= 999.9
end do
colymax(1)= -10.0*pcol(1)
colymin(2)= 10.0*pcol(2)
colxmin(3)= 10.0*pcol(3)
colxmax(4)= -10.0*pcol(4)
colxmax(6)= 10.0*pcol(6)
colxmin(7)= -10.0*pcol(7)
if(inunit.eq.5) then
write(6,'(a)')' Impose Q1 Exit Aperture Restriction ?(y/n): '
endif
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
Q1CLIP = .TRUE.
else
Q1CLIP = .false.
endif
if(inunit.eq.5) then
write(6,'(a)')'Use 3rd order fringe eff. in quads ?(y/n): '
endif
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
QFRINGE = .TRUE.
else
QFRINGE = .false.
endif
if(inunit.eq.5) then
write(6,'(a)')'Use 3 quad tune (high energy) ?(y/n): '
endif
read(inunit,'(a)')answer
if ((answer.eq.'y').or.(answer.eq.'Y')) then
THREEQUAD= .TRUE.
else
THREEQUAD = .false.
endif
if(.not.THREEQUAD) then
tip2=0.0
grad2=0.0
endif
if(inunit.eq.5) write(6,'(a)')'Which solenoid simulation?
>F=full,S=simple,N=none'
read(inunit,'(a)')answer
if ((answer.eq.'f').or.(answer.eq.'F')) then
SOLSIM = 1
write(6,'(a)')'Complex field will be used'
read(inunit,*) BSolLong
write(6,'(a)')' Enter solenoid
>dx,dy,dxp,dyp (mm,radians)'
read(inunit,*) soldx, soldy, soldxp, soldyp
c LSol=500. !propagate further through field
LSol=450. !propagate further through field
elseif ((answer.eq.'s').or.(answer.eq.'S')) then
write(6,*) 'Simple solenoid will be used'
SOLSIM = -1
write(6,'(a)')' Enter Solenoid Field (T)'
read(inunit,*) BSolLong
write(6,'(a)')' Enter solenoid
>dx,dy,dxp,dyp (mm,radians)'
read(inunit,*) soldx, soldy, soldxp, soldyp
c LSol=125. !DJG: note that Lsol is Leff/2 in this context
LSol=155. !DJG: Leff is 310 mm (from field map).
!But we only go through 1/2 of solenoid
else
write(6,'(a)') 'No solenoid effect applied'
solsim=0
BSolLong=0.0
soldx=0
soldy=0
soldxp=0
soldyp=0
LSol=125.
endif
C DJG calculate input file-dep. distances
if(QFRINGE) then
TM1 = XQ1 - (LQ1-r1)/2. - r1
M1M2 = XQ2 - (LQ2-r2)/2. - r2 - XQ1 - (LQ1-r1)/2.- r1
M2M3 = XQ3 - (LQ3-r3)/2. - r3 - XQ2 - (LQ2-r2)/2.- r2
c M3D = XDET - XQ3 - (LQ3-r3)/2. - r3
M3D(1) = XDET(1) - XQ3 - (LQ3-r3)/2. - r3
M3D(2) = XDET(2) - XQ3 - (LQ3-r3)/2. - r3
M1C = TC - Tm1 - (LQ1+r1)
CM2 = XQ2 - (LQ2-r2)/2. - r2 - col6z
endif
SM1 = TM1 - LSol
C if(inunit.eq.5) write(6,'(a)')' Enter Solenoid Field (T)'
C read(inunit,*) BSolLong
C
C if(inunit.eq.5) write(6,'(a)')' Enter solenoid dx,dy,dxp,dyp (mm,radians)'
C read(inunit,*) soldx, soldy, soldxp, soldyp
c---------------------------------------------------------------
write(6,811) b_energy, PBEAM
811 format(' Beam Energy: ',f6.3,' Momentum: ',f6.3)
write(6,812) mcur1,mcur2
812 format(' Quad Currents: ',2f8.2, ' Amps')
write(6,813) tip1,tip2,tip3
813 format(' Quad Tip Fields: ',3f10.4, ' Tesla')
write(6,814) grad1,grad2,grad3
814 format(' Quad Gradients: ',3f10.5, ' T/mm')
GdotL1= (grad1*10.0*10.0)*(LQ1/10.0)
GdotL2= (grad2*10.0*10.0)*(LQ2/10.0)
write(6,8141) GdotL1,GdotL2
8141 format(' Quad GL ',2f10.4, ' (kG/cm)*cm')
write(6,815) NUMEVENTS
815 format(' Generating ',i10,' events.')
write(6,816) climit,crange
816 format(' Scattering Angle Range: ',2f10.4)
write(6,817) X0,Y0,TX0,TY0
817 format(' Beam Offsets: Position x/y and Angle x/y',4f10.4)
write(6,818) SIGX,SIGY,SIGXP,SIGYP
818 format(' Beam size(X/Y) and divergence(X/Y):',4f10.4)
write(6,819) scflag(2)
819 format(' Radiative Corrections?: ',L1)
write(6,820) scflag(1)
820 format(' MCS Before Interaction?: ',L1)
write(6,821) scflag(3)
821 format(' MCS After Interaction?: ',L1)
write(6,822) BINDE
822 format(' Levchuk Effect??: ',L1)
write(6,823) xtgt
823 format(' Target Thickness: ',f10.3)
write(6,824) pcol
824 format(' Collimators at:',7f6.1)
write(6,825) Q1CLIP
825 format(' Q1 Exit Aperture Clipping?? ',L1)
write(6,826) QFRINGE
826 format(' Quad fringe field effects?? ',L1)
write(6,827) THREEQUAD
827 format(' 3 QUAD tune?? ',L1)
write(6,828) BSolLong
828 format(' Solenoid Field (Tesla): ',f4.2)
write(6,829) soldx, soldy, soldxp, soldyp
829 format(' Solenoid offsets (x,y in mm, dxp,dyp in radians):',2f6.2,2f8.5)
c if(abs(BSolLong).le.MinSolFld) write(6,828)
c 828 format(' ... Solenoid field below cut. No solenoid effect applied.')
c---------------------------------------------------------------
IHD=0
C
C IHD DETERMINES THE HEADING TO BE PRINTED
C
IF(TAILS) THEN
IHD=IHD+1
ELSE
IHD=IHD+2
ENDIF
C
C PUT THE DETECTOR POSITION INTO THE HEADING
C
IF(BINDE) THEN
WRITE(HEADING,50) TITLE(IHD),Y0,X0,PBEAM
50 FORMAT('H2 P DIST',A12,1X,'Y',F5.2,' MM',' X',
> F5.2,' MM',F7.3,' GEV')
ELSE
WRITE(HEADING,51) TITLE(IHD),Y0,X0,PBEAM
51 FORMAT('FREE E: ',A12,1X,'Y', F5.2,' MM', ' X',
> F5.2,' MM',F7.3,' GEV')
ENDIF
C
C INITIALIZE EVERYTHING
ccc
c call hropen(16,'ntuple','scr2:cebaf.ntuple','nq',4096,istat) !Ntuple
if(nrarg.gt.0) then
write(filename,'("hbook/",(a))') argin(1)
i=index(filename,' ')
write(filename(i:),'(".hbook")')
call hropen(17,'histo',filename,'nq',4096,istat) !Histo
else
call hropen(17,'histo',
& 'hbook/mollermcb.hbook','nq',4096,istat) !Histo
endif
ccc
CALL HBOOKM(npads,p_up)
ccc
ccc initialize ntuple stored in file cebaf.ntuple M.L. 12.94
ccc
c call hcdir('//ntuple',' ') !Ntuple
c call hbnt(10,'cebaf_ntuple',' ') !Ntuple
c call hbset('bsize',4096*20,ierr) !Ntuple
c call hbname(10,' ', 0, '$clear') !Ntuple
c call hbname(10,'events',ntuple,'collx1,colly1,collx2
c * ,colly2,detx1,dety1,detx2,dety2,Axxp,Axxm,Axyp,Axym
c * ,Azzp,Azzm,Axxp2,Axxm2,Axyp2,Axym2,Azzp2,Azzm2') !Ntuple
c call hcdir('//histo',' ') !Histo
ccc
DO J=1,3
DO I=1,2
BLKTL(I,J)=0.
BLKTL2(I,J)=0.
ENDDO
ENDDO
avgA=0.0
avgA2=0.0
WEITOT=0.
WEITO2=0.
nplt= 0
nplt2= 0
nplt3= 0
nplt4= 0
C
ccc
ccc Define the number of event trials (need many to simulate bound e- tgt)
ccc Code works up to 4'000'000 events (with ntuples). There was a problem
ccc at the beginning with the CERN library that it crashed after 1'000'000
ccc which was solved by increasing iquest(10) from 16000 to 64000 and
ccc call hropen(16,'ntuple','scr2:cebaf.ntuple','nq',4096,istat)
ccc instead of ...'nq',1024,istat) the default. M.L. 11.96
ccc
C
c
c________________________________________________________________________
c________________________________________________________________________
c________________________________________________________________________
C
C LOOP THROUGH EVENTS
c________________________________________________________________________
C
DO 900 ITRIAL=1,NUMEVENTS
ccc
ccc Get some feed-back how many event are already done M.L. 11.96
ccc
if(mod(itrial,100000) .eq. 0) then
write(6,*) itrial
endif
ccc
twomols= .true.
singlemol(1)= .true.
singlemol(2)= .true.
C
C GENERATE BEAM PHASE SPACE FOR THE INCIDENT VECTOR
C
CALL BPHASE(VEC(1,1),PBEAM)
C
C ACCOUNT FOR VERTICAL AND HORIZONTAL BEAM OFFSET
C
call hfdp1(215,vec(1,1),1.0D0) !histo beamspot
call hfdp1(216,vec(2,1),1.0D0) !histo
C
c________________________________________________________________________
call hfill(10,1.,0.0,1.0) !count- after moller scattering
do j=1,2
if (singlemol(j)) call hfill(10,1.0+j,0.,1.)
end do
c________________________________________________________________________
c________________________________________________________________________
C
C CHOOSE A RANDOM INTERACTION DEPTH IN THE TARGET
C
DEPTH = 0.999*RAN3(II)+0.0005
XDEPTH=XTGT*DEPTH
C
C MULTIPLE SCATTER AND BREM THE E- TO THE INTERACTION POINT
C
if (scflag(1)) then
CALL BRMSCT(VEC(1,1),XDEPTH,1,TAILS,TCKBRM,WGTBM)
PB=SQRT(VEC(4,1)**2+VEC(5,1)**2+VEC(6,1)**2)
if (PB .lt. 0.1*b_energy) goto 900 !low energy cutoff
else
call zdrift(vec(1,1),xdepth)
wgtbm= 1.0
endif
C
c________________________________________________________________________
call hfill(10,4.,0.0,1.0) !count- after moller scattering
do j=1,2
if (singlemol(j)) call hfill(10,4.0+j,0.,1.)
end do
c________________________________________________________________________
c________________________________________________________________________
C
C GENERATE MOLLER SCATTERING
C
CALL VECGEN(VEC,XTGT,CST,TGEN,PHIGEN,PEL,
> WEIGHT,PWT,BINDE,scflag(2), PPOL,
> PUNPOL, TGTPOL, targetpol)
call hfdp1(219,cst,0.65D0) !Histo cosine of CM scat. angle
mytheta = acos(cst)*57.2958d0
myphi = phigen*57.2958d0
call hfdp1(401,mytheta,weight)
call hfdp1(4030,mytheta,1.0d0)
call hfdp1(402,myphi,weight)
wgtcx= weight !save unpol'd x-section weight
DO J=1,2 !weights for individual electrons
weights(j) = (WGTBM*WEIGHT)/(1.D0*NUMEVENTS)
PB=SQRT(VEC(4,j)**2+VEC(5,j)**2+VEC(6,j)**2)
if (PB .lt. 0.1*b_energy) then !low energy cutoff
twomols= .false. ! should be false
singlemol(j)= .false. ! should be false
endif
ENDDO
c weight for coincidence electron events
WEIGHT = (WGTBM*WEIGHT)/(1.D0*NUMEVENTS)
C
c________________________________________________________________________
call hfill(10,7.,0.0,1.0) !count- after moller scattering
do j=1,2
if (singlemol(j)) call hfill(10,7.0+j,0.,1.)
end do
c________________________________________________________________________
C
C MULTIPLE SCATTER AND BREM THE OUTGOING ELECTRONS
C
XDEPTH=XTGT*(1.-DEPTH)
WGTML=1.D0
DO J=1,2
if (singlemol(j)) then
if (scflag(3)) then
CALL BRMSCT(VEC(1,J),XDEPTH,1,TAILS,TCKBRM,WGT)
else
call zdrift(vec(1,j),xdepth)
wgt= 1.0
endif
else
endif
IF(WGT.LE.0.) then
twomols= .false.
singlemol(j)= .false.
weights(j)= 0.0
weight= 0.0
else
WGTML=WGTML*WGT
weights(j)= weights(j)*wgt
PB=SQRT(VEC(4,j)**2+VEC(5,j)**2+VEC(6,j)**2)
if (PB .lt. 0.1*b_energy) then !low energy cutoff
twomols= .false.
singlemol(j)= .false.
endif
endif
ENDDO
if (twomols) then
WEIGHT=WEIGHT*WGTML
call hfdp1(4031,mytheta,1.0d0)
endif
C
call hfill(10,10.,0.0,1.0) !count- begin tracking
do j=1,2
if (singlemol(j)) call hfill(10,10.0+j,0.,1.)
end do
c________________________________________________________________________
C
C PROPAGATE FROM TARGET, AT CENTER OF SOLENOID, TO END OF SOLENOID
C
if(solsim.lt.0) then
do j=1,2
call misalign(vec(1,j), soldx, soldy, soldxp, soldyp)
call solenoid(vec(1,j),Lsol,BSolLong,10)
call misalign(vec(1,j),-soldx,-soldy,-soldxp,-soldyp)
enddo
elseif(solsim.gt.0) then
do j=1,2
call misalign(vec(1,j),soldx,soldy,soldxp,soldyp)
call solswim(BsolLong,vec(1,j),vec(2,j),vec(3,j),vec(4,j),
> vec(5,j),vec(6,j))
c write(6,*) 'cheesy poofs',j,vec(1,j),vec(2,j),vec(3,j)
call misalign(vec(1,j),-soldx,-soldy,-soldxp,-soldyp)
enddo
endif
c________________________________________________________________________
C
C PROPAGATE TO FIRST QUAD MAGNET
C
if(solsim.ne.0) then
ZDIST = SM1
else
ZDIST = TM1
endif
c write(6,*) 'cheesy poofs',ZDIST,SM1
DO J=1,2
if (singlemol(j)) then
if(solsim.ne.0) then
CALL ZDRIFT(VEC(1,J),ZDIST+(Lsol-VEC(3,J)))
else
CALL ZDRIFT(VEC(1,J),ZDIST)
endif
CALL HFDP2(91,VEC(1,J),VEC(2,J),1.D0) !histo Q1 dist.
endif
ENDDO
C
c________________________________________________________________________
C
C PROPAGATE THROUGH THE 1ST QUAD MAGNET
C
c iii= iii +1
DO J=1,2
if (singlemol(j)) then
ccc
ccc Misalignment of the quad1, done by shifting beam forth and back
ccc Not the most elegant way but it works M.L. 9.96
ccc
ccc vec(1,2)=vec(1,2)-2 ! M.L. 9.96 to simulate y shift of quad
ccc vec(1,1)=vec(1,1)-2 ! M.L. 9.96 to simulate x shift of quad
if(QFRINGE) then
CALL QUADRUPOLE_FRINGE(VEC(1,J),LQ1,r1,tip1)
else
quad1x=0.01 !Mis-alignment of electron = -quad shift JAM
quad1y=0.12 !
CALL quad_misalign(vec(1,j),quad1x,quad1y) !Misaligns electrons JAM
CALL QUADRUPOLE(VEC(1,J),LQ1,grad1,10)
CALL quad_misalign(vec(1,j),-quad1x,-quad1y) !re-aligns electrons JAM
endif
ccc vec(1,1)=vec(1,1)+2 ! M.L. 9.96 to simulate x shift of quad
ccc vec(1,2)=vec(1,2)+2 ! M.L. 9.96 to simulate y shift of quad
quad1exit(1,j)= vec(1,j)
quad1exit(2,j)= vec(2,j)
CALL HFDP2(92,VEC(1,J),VEC(2,J),1.D0) !histo post-Q1 dist.
if( dabs(mytheta - 90.0) .lt. 0.1) then
CALL HFDP2(921,VEC(1,J),VEC(2,J),1.D0) !histo Q1 exit for 90 degrees
endif
endif
ENDDO
c________________________________________________________________________
C
c-- Account for any clipping at Q1 exit
if (Q1CLIP) then
DO J=1,2
if (singlemol(j)) then
rq1_1= sqrt(quad1exit(1,j)**2 + quad1exit(2,j)**2)
if( rq1_1 .ge. q1rad ) then
singlemol(j)=.false.
twomols= .false.
endif
endif
enddo
endif
if (twomols) then
call hfdp1(4032,mytheta,1.0d0)
endif
c________________________________________________________________________
C
C PROPAGATE TO THE COLLIMATOR 1
C
ZDIST = M1C
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
IF (singlemol(j)) THEN
CALL HFDP2(100,VEC(1,J),VEC(2,J),1.D0) !histo pre-coll dist.
CALL HFDP2(1001,VEC(1,J),VEC(2,J),1.D0)
col(1,1,j)=VEC(1,J)
col(2,1,j)=VEC(2,J)
NPLT=NPLT+1
c ntuple(2*j-1)=vec(1,j) !Ntuple coll x1 x2
c ntuple(2*j)=vec(2,j) !Ntuple coll y1 y2
ENDIF
ENDDO
c________________________________________________________________________
C
C PROPAGATE TO THE COLLIMATOR 3
C
ZDIST = col3z-col1z
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
IF (singlemol(j)) THEN
CALL HFDP2(100,VEC(1,J),VEC(2,J),1.D0) !histo pre-coll dist.
CALL HFDP2(1003,VEC(1,J),VEC(2,J),1.D0)
col(1,3,j)=VEC(1,J)
col(2,3,j)=VEC(2,J)
NPLT=NPLT+1
c ntuple(2*j-1)=vec(1,j) !Ntuple coll x1 x2
c ntuple(2*j)=vec(2,j) !Ntuple coll y1 y2
ENDIF
ENDDO
c________________________________________________________________________
C
C PROPAGATE TO THE COLLIMATOR 5
C
ZDIST = col5z-col3z
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
IF (singlemol(j)) THEN
CALL HFDP2(100,VEC(1,J),VEC(2,J),1.D0) !histo pre-coll dist.
CALL HFDP2(1005,VEC(1,J),VEC(2,J),1.D0)
col(1,5,j)=VEC(1,J)
col(2,5,j)=VEC(2,J)
NPLT=NPLT+1
c ntuple(2*j-1)=vec(1,j) !Ntuple coll x1 x2
c ntuple(2*j)=vec(2,j) !Ntuple coll y1 y2
ENDIF
ENDDO
c________________________________________________________________________
C
C PROPAGATE TO THE COLLIMATOR 6
C
ZDIST = col6z-col5z
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
IF (singlemol(j)) THEN
CALL HFDP2(100,VEC(1,J),VEC(2,J),1.D0) !histo pre-coll dist.
CALL HFDP2(1006,VEC(1,J),VEC(2,J),1.D0)
col(1,6,j)=VEC(1,J)
col(2,6,j)=VEC(2,J)
NPLT=NPLT+1
c ntuple(2*j-1)=vec(1,j) !Ntuple coll x1 x2
c ntuple(2*j)=vec(2,j) !Ntuple coll y1 y2
ENDIF
ENDDO
c________________________________________________________________________
C
C SIMULATE THE COLLIMATOR
C
ccc
ccc No Collimator is used when reanalyzing with ntuples is done. We can
ccc then simulate the collimator by cuts in the ntuple analysis M.L. 11.96
ccc
ienergy = 1
DO J=1,2
if (singlemol(j)) then
singlemol(j)= all_coll(col(1,1,j))
twomols= ( twomols .and. singlemol(j) )
endif
c if (singlemol(j)) then
c + call collimate(vec(1,j),singlemol(j),twomols,ienergy)
IF ((NPLT2.LT.5000).and.(singlemol(j))) THEN
CALL HFDP2(101,VEC(1,J),VEC(2,J),1.D0) !histo post-coll dist.
NPLT2=NPLT2+1
ENDIF
ENDDO
if (twomols) then
call hfdp1(4033,mytheta,1.0d0)
endif
ccc
c________________________________________________________________________
C
C PROPAGATE TO THE 2ND QUAD MAGNET
C
ZDIST = CM2
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
ENDDO
c________________________________________________________________________
C
C PROPAGATE THROUGH THE 2ND QUAD MAGNET
C
DO J=1,2
if (singlemol(j)) then
CALL HFDP2(1011,VEC(1,J),VEC(2,J),1.D0) !at Q2 aperture.
endif
enddo
if(dabs(VEC(1,1)).gt.100..or.dabs(vec(1,2)).gt.100.) go to 900
DO J=1,2
if (singlemol(j)) then
ccc
ccc Missalignment of the quad2, done by shifting beam forth and back
ccc Not the most elegant way but it works M.L. 9.96
ccc
ccc vec(1,2)=vec(1,2)-5 ! M.L. 9.96 to simulate y shift of quad
vec(1,j)=vec(1,j) + quad2offset ! simulate x shift of quad
quad2entrance(1,j)= vec(1,j)
quad2entrance(2,j)= vec(2,j)
c print *,vec(1,j),j
if(tip2.gt.0.0) then
if(QFRINGE) then
CALL QUADRUPOLE_FRINGE(VEC(1,J),LQ2,r2,tip2)
else
CALL QUADRUPOLE(VEC(1,J),LQ2,grad2,10)
endif
else
if(QFRINGE) then
zdist=LQ2+r2
else
zdist = LQ2
endif
call zdrift(vec(1,J),zdist)
endif
vec(1,j)=vec(1,j) - quad2offset ! simulate x shift of quad
quad2exit(1,j)= vec(1,j)
quad2exit(2,j)= vec(2,j)
ccc vec(1,2)=vec(1,2)+5 ! M.L. 9.96 to simulate y shift of quad
c write(6,*) vec(1,j),vec(2,j),vec(3,j),vec(4,j),vec(5,j),vec(6,j),j,'l'
endif
ENDDO
c________________________________________________________________________
C
C PROPAGATE TO THE 3rd QUAD MAGNET
C
ZDIST = M2M3
DO J=1,2
if (singlemol(j)) call zdrift(vec(1,j),zdist)
ENDDO
c________________________________________________________________________
C