commit a8b83521908d8cc37dffe3315248b29d3384dd91 Author: Thomas Prokscha Date: Thu Jun 17 13:05:47 2004 +0000 Added to repository. diff --git a/geant3/inp/geant_lemsr.input b/geant3/inp/geant_lemsr.input new file mode 100644 index 0000000..5385139 --- /dev/null +++ b/geant3/inp/geant_lemsr.input @@ -0,0 +1,114 @@ +plot 1 +rung 3022 +trig 10 +rndm 9871 59367 +COMM 1 2 3 4 5 6 7 8 9 +sets 5 -74 1 9 0 0 0 0 0 +bfie 0. 0. 0. 0. 0. +spar ' e+' ' e-' ' gam' +cuts 0.0001 0.0001 0.01 0.01 0.00001 +COMM CutGam CutEle CutNeu CutHad CutMuo +ntva 'init' 'time' 'flag' 'scie' 'scoe' 'ltrb' +geom 'mcps' 'samp' 'efla' 'cryo' 'crsh' 'saal' 'ru11' +COMM geom 'mcps' 'efla' 'ru11' 'mcp2' 'mcpa' +view 1 60 40 0 10 10 4 4 +stop + + +This is the Format Free Input File for GEANT_LEMSR.FOR: +======================================================= + +how to use the data cards: + +COMM ---> comment line, no effect when reading + +plot 0 ---> no graphics. +plot 1 ---> graphics on screen, and also information on detector hits. +plot 2 ---> graphics on screen and PS file. + +rung 12 ---> define GEANT run number; will be appended to the NTP filename + +trig 100000 ---> number of events to be processed; on PSICL1 10^5 takes + up to 3:30 min CPU time. + + 1 2 3 4 5 6 7 8 9 +sets 5 -74 2 5 0 1 0 2 5 ---> first number : particle ID, 5 = mu+ (2 = e+) + ID = 500; user defined mu+ + with anistropic decay prob. + second number: if == 0, michel-distributed at + MCP2 position + if ne 0, start position is at z=sets(2) + (here 74 cm upstream of MCP tube center) + with beam momentum + sets(3).sets(4)sets(5), + this is 2.50 MeV/c in this example + sets(6): gaussian beam divergence (in degree) + sets(7): if = 1, read initial values from file + MC$INP:GEANT_LEMSR_INPUT.DAT + sets(8): gaussian dp/p distribution, + sigma_p = sets(8).sets(9)% +spar ' e+' ' e-' ' gam' --> select secondary particles to be tracked. + +BFIE 0. 0. 100. 90. 50. ---> Sample B-field (Gauss), Bx, By, Bz, and muon polarization=90% + and muon polarization=50% in sample. + +ntva 'time' 'flag' ---> see comments below + + +view 1 40 40 0 10 10 3 3 ---> in case of graphical output define view angle, + (60)(40)(0)(10)(10)(2)(2) position and scaling of graph. + default values view(1): if == 1 then use the command lines + parameters for graphical output; + otherwise use default parameter + view(2): theta angle (degree) + view(3): phi angle + view(4): psi angle + view(5): x0, x-position in HIGZ window + view(6): y0, y-position in HIGZ window + view(7): 10*x-scaling, view(7)=3 <==> xscal=0.3 + view(8): 10*y-scaling, view(8)=3 <==> yscal=0.3 +geom 'mcp2' 'efla' ---> see comment below + +stop ---> stops reading this file and continues the program. + +c------------------------------------------------------------------------------ +c +c flags for geometrical setup; it is possible to mount or remove the +c +c 'mcp2' MCP2 plates +c 'mcpa' MCP2 anode +c 'mcps' MCP2 stainless steel vacuum tube +c 'samp' sample (this is at the moment the sample holder) +c 'cryo' material of the cryostat without shielding +c 'crsh' cryo shield +c 'efla' the 100 CF flange at the end of the sample/MCP2 tube +c 'ru11' use MCP setup of Run 11 (delay line anode) +c 'ru10' use MCP setup of Run 10 and before (WSZ detector) +c has no influence for cryo geometry +c 'saal' uses Al instead of Cu for sample holder plates +c +c the volumes are activated when printing +c +c geom 'mcp2' 'efla' for example when data cards are requested +c +c--------------------------------------- +c +c flags to define NTuple variables for CWN +c +c 'init' write initial parameter to NT +c 'code' initial particle ID to NT +c 'endp' end position and energy of e+ +c 'time' t-mcp and t-sci +c 'scip' position when particle hits inner Sc. +c 'scie' energy loss of particles in inner sc. +c 'scop' position when particle hits outer Sc. +c 'scoe' energy loss of particles in outer sc. +c 'ltrb' energy losses in single detectors (inner and outer sc.) +c 'flag' Vol. numbers and particle codes; needed to get the single +c decay spectra +c 'nhit' Number of Hits in detectors MCP, SCI and SCO +c 'posm' Generated momentum of decay positron +c 'mcpe' energy loss in MCP detector +c 'allv' use all possible NT variables +c +c-------------------------------------------------------------------------------- diff --git a/geant3/paw/fit_prec.for b/geant3/paw/fit_prec.for new file mode 100644 index 0000000..1404737 --- /dev/null +++ b/geant3/paw/fit_prec.for @@ -0,0 +1,35 @@ + function fit_prec(t) + + real t, t0 + vector tzero(1) + common/pawpar/par(6) + + +c t: time [ns] + +c par(1): Background [1] +c par(2): normalization [1] +c par(3): asymmetry [1] +c par(4): frequency [MHz] +c par(5): phase [degree] +c par(6): muon lifetime [microsec] +c par(7): relaxation [1/microsec] + +c um richtige Fitergebnisse zu erhalten, muss der Fitroutine der Zeitnullpunkt +c des Spektrums mitgeteilt werden: + + t0 = tzero(1) + t = t-t0 + + + twopi = 6.283185 + fit_prec = par(1)+ ! Background + + par(2)* ! normalization + + exp(-t/(1000.*par(6)))* ! decay + + ( 1.+par(3)* ! asymmetry +c + exp(-t*(par(7)/1000.))* ! relaxation + + cos( twopi*(par(4)/1000. *t + par(5)/360.)) ! oscillation + + ) + + end + diff --git a/geant3/paw/fit_prec2_ge.kumac b/geant3/paw/fit_prec2_ge.kumac new file mode 100644 index 0000000..b9bb649 --- /dev/null +++ b/geant3/paw/fit_prec2_ge.kumac @@ -0,0 +1,383 @@ +if [1].EQ.'?' .OR. [1].EQ.' ' then + mess '===============================================================================' + mess 'FIT_PREC IDLeft B fromFit toFit t0' + mess + mess 'Kumac to fit the GEANT simulated decay spectra of L and R.' + mess 'If IDleft is given, then it is assumed that' + mess ' IDright = IDleft+3,' + mess + mess 'B : magnetic field in gauss' + mess 'fromFit : lower boundary for fit' + mess 'toFit : upper boundary for fit' + mess 't0 : time zero' + mess + mess '===============================================================================' + exitm +endif +* +*=========================================================================== +current = $hcdir() + +* define Histo ID's + + idl = [1] + idt = [idl]+2 + idr = [idl]+3 + idb = [idl]+4 +* +* check if histograms exist in PAW memory +* + cd //pawc + + if ( $hexist([idl]) .eq. 0 ) then + mess + mess Did not found Left histogram ID [idl] + mess + exitm + endif + if ( $hexist([idt]) .eq. 0 ) then + mess + mess Did not found Top histogram ID [idt] + mess + exitm + endif + if ( $hexist([idr]) .eq. 0 ) then + mess + mess Did not found Right histogram ID [idr] + mess + exitm + endif + if ( $hexist([idb]) .eq. 0 ) then + mess + mess Did not found Bottom histogram ID [idb] + mess + exitm + endif + + shift + +* transfer magnetic field strength: + + B = [1] + +* get region of histo to be used for fit: + + if ($index([2],'.').NE.0) then + von = $rsigma([2]) + else + von = [2] + endif + if ($index([3],'.').NE.0) then + bis = $rsigma([3]) + else + bis = [3] + endif + +* get rebinning factor and fit options: + + t0 = [4] + + if ( [t0] .eq. ' ' .or. [t0] .eq. '!' ) then + mess + mess Missing time zero ! + mess + exitm + endif + if ( $vexist(tzero) ) then + v/de tzero + endif + v/cre tzero(1) + v/inp tzero(1) [t0] + +* check if added histos are to be used and read in histos: + +* do the fits: + + freq = [B] * 0.013554 | precession frequency in MHz + freqmin = 0.8*[freq] + freqmax = 1.2*[freq] + freqstep = 0.005*[freq] +* + if ( $vexist(par) ) then + v/de par + endif + if ( $vexist(pmin) ) then + v/de pmin + endif + if ( $vexist(pmax) ) then + v/de pmax + endif + if ( $vexist(step) ) then + v/de step + endif + if ( $vexist(errpar) ) then + v/de errpar + endif +* +* +* +* BKG N0 Asym mu+ freq. Phase tau_mu lambda +* ======================================================= +* +* v/cre par(7) r 0. 500. 0.25 0. 0. 2.19703 0. +* v/cre pmin(7) r -10. 0. 0. 0. -400. 0. 0. +* v/cre pmax(7) r 100. 10000. 0.60 0. +400. 3. 100. +* v/cre step(7) r 0. 10. 0.01 0. 5. 0. 0. +* v/cre errpar(7) r 1. 1. 1. 1. 1. 1. 1. + v/cre par(6) r 0. 500. 0.25 [freq] 0. 2.19703 + v/cre pmin(6) r -10. 0. 0. [freqmin] -400. 0. + v/cre pmax(6) r 100. 10000. 0.40 [freqmax] +400. 3. + v/cre step(6) r 0. 10. 0.01 [freqstep] 5. 0. + v/cre errpar(6) r 0. 1. 1. 1. 1. 0. + +* + tdec0 = tzero(1) | vector tzero filled in FIT_PREC.KUMAC + ttt = Fit from [von] to [bis], t0 = [tdec0] + title [ttt] +* +* + hi/fit [idl]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar +* +* store fit result and error for A_Left, frequency and relaxation rate +* + a_l = par(3) + a_lerr = errpar(3) + a_l_weight = [a_l]/([a_lerr]*[a_lerr]) + inv_a_lerr = 1./([a_lerr]*[a_lerr]) + f_l = par(4) + f_lerr = errpar(4) + f_l_weight = [f_l]/([f_lerr]*[f_lerr]) + inv_f_lerr = 1./([f_lerr]*[f_lerr]) +* r_l = par(7) +* r_lerr = errpar(7) +* r_l_weight = [r_l]/([r_lerr]*[r_lerr]) +* inv_r_lerr = 1./([r_lerr]*[r_lerr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_l = par(8) +* beta_lerr = errpar(8) +* beta_l_weight = [beta_l]/([beta_lerr]*[beta_lerr]) +* inv_beta_lerr = 1./([beta_lerr]*[beta_lerr]) +* endif +* +* +* initialize the phase for the top detector +* + goto right +* + phase = par(5) + phase = [phase]+90. + v/inp par(5) [phase] + hi/fit [idt]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar +* +* store fit result and error for A_top +* + a_t = par(3) + a_terr = errpar(3) + a_t_weight = [a_t]/([a_terr]*[a_terr]) + inv_a_terr = 1./([a_terr]*[a_terr]) + f_t = par(4) + f_terr = errpar(4) + f_t_weight = [f_t]/([f_terr]*[f_terr]) + inv_f_terr = 1./([f_terr]*[f_terr]) +* r_t = par(7) +* r_terr = errpar(7) +* r_t_weight = [r_t]/([r_terr]*[r_terr]) +* inv_r_terr = 1./([r_terr]*[r_terr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_t = par(8) +* beta_terr = errpar(8) +* beta_t_weight = [beta_t]/([beta_terr]*[beta_terr]) +* inv_beta_terr = 1./([beta_terr]*[beta_terr]) +* endif +* +* +* +* initialize the phase for the right detector +* +right: + phase = par(5) + phase = [phase]+180. + v/inp par(5) [phase] + hi/fit [idr]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar + +* store fit result and error for A_right +* + a_r = par(3) + a_rerr = errpar(3) + a_r_weight = [a_r]/([a_rerr]*[a_rerr]) + inv_a_rerr = 1./([a_rerr]*[a_rerr]) + f_r = par(4) + f_rerr = errpar(4) + f_r_weight = [f_r]/([f_rerr]*[f_rerr]) + inv_f_rerr = 1./([f_rerr]*[f_rerr]) +* r_r = par(7) +* r_rerr = errpar(7) +* r_r_weight = [r_r]/([r_rerr]*[r_rerr]) +* inv_r_rerr = 1./([r_rerr]*[r_rerr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_r = par(8) +* beta_rerr = errpar(8) +* beta_r_weight = [beta_r]/([beta_rerr]*[beta_rerr]) +* inv_beta_rerr = 1./([beta_rerr]*[beta_rerr]) +* endif +* +* +* +* initialize the phase for the bottom detector +* + goto end +* + phase = par(5) + phase = [phase]+90. + v/inp par(5) [phase] + hi/fit [idb]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar +* +* store fit result and error for A_bottom +* + a_b = par(3) + a_berr = errpar(3) + a_b_weight = [a_b]/([a_berr]*[a_berr]) + inv_a_berr = 1./([a_berr]*[a_berr]) + f_b = par(4) + f_berr = errpar(4) + f_b_weight = [f_b]/([f_berr]*[f_berr]) + inv_f_berr = 1./([f_berr]*[f_berr]) +* r_b = par(7) +* r_berr = errpar(7) +* r_b_weight = [r_b]/([r_berr]*[r_berr]) +* inv_r_berr = 1./([r_berr]*[r_berr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_b = par(8) +* beta_berr = errpar(8) +* beta_b_weight = [beta_b]/([beta_berr]*[beta_berr]) +* inv_beta_berr = 1./([beta_berr]*[beta_berr]) +* endif +* +* +* plot result: + +end: + zone 2 2 + + hi/pl [idl] e + hi/pl [idt] e + hi/pl [idr] e + hi/pl [idb] e +* +*-------------------------------- +* +* calculate the weighted mean and the "inner" and "outer" standard deviation +* fo asymmetry +* + sum1 = [a_l_weight] + [a_r_weight] + sum2 = [inv_a_lerr] + [inv_a_rerr] + a_mean = [sum1]/[sum2] + in_a_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* now calculate the outer standard deviation +* + dev_l = ([a_mean]-[a_l])*([a_mean]-[a_l])*[inv_a_lerr] + dev_r = ([a_mean]-[a_r])*([a_mean]-[a_r])*[inv_a_rerr] + sum3 = [dev_l] + [dev_r] + out_a_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +*---------------------------------- +* +* mean of frequency +* + sum1 = [f_l_weight] + [f_r_weight] + sum2 = [inv_f_lerr] + [inv_f_rerr] + f_mean = [sum1]/[sum2] + in_f_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* + dev_l = ([f_mean]-[f_l])*([f_mean]-[f_l])*[inv_f_lerr] + dev_r = ([f_mean]-[f_r])*([f_mean]-[f_r])*[inv_f_rerr] + sum3 = [dev_l] + [dev_r] + out_f_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +*-------------------------------------- +* +* mean of relaxation +* +* sum1 = [r_l_weight] + [r_t_weight] + [r_r_weight] + [r_b_weight] +* sum2 = [inv_r_lerr] + [inv_r_terr] + [inv_r_rerr] + [inv_r_berr] +* r_mean = [sum1]/[sum2] +* in_r_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* +* dev_l = ([r_mean]-[r_l])*([r_mean]-[r_l])*[inv_r_lerr] +* dev_t = ([r_mean]-[r_t])*([r_mean]-[r_t])*[inv_r_terr] +* dev_r = ([r_mean]-[r_r])*([r_mean]-[r_r])*[inv_r_rerr] +* dev_b = ([r_mean]-[r_b])*([r_mean]-[r_b])*[inv_r_berr] +* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] +* out_r_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +* mean of stretched exponent +* +* +* if ( [opt] .eq. 'strexp' ) then +* sum1 = [beta_l_weight] + [beta_t_weight] + [beta_r_weight] + [beta_b_weight] +* sum2 = [inv_beta_lerr] + [inv_beta_terr] + [inv_beta_rerr] + [inv_beta_berr] +* beta_mean = [sum1]/[sum2] +* in_beta_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* +* dev_l = ([beta_mean]-[beta_l])*([beta_mean]-[beta_l])*[inv_beta_lerr] +* dev_t = ([beta_mean]-[beta_t])*([beta_mean]-[beta_t])*[inv_beta_terr] +* dev_r = ([beta_mean]-[beta_r])*([beta_mean]-[beta_r])*[inv_beta_rerr] +* dev_b = ([beta_mean]-[beta_b])*([beta_mean]-[beta_b])*[inv_beta_berr] +* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] +* out_beta_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +* endif +* +*---------------------------------------- +* +* output fit parameter definitions to text window: + + mess + mess 'par(1): Background [1]' + mess 'par(2): normalization [1]' + mess 'par(3): asymmetry [1]' + mess 'par(4): frequency [MHz]' + mess 'par(5): phase [degree]' + mess 'par(6): muon lifetime [microSec]' + mess 'par(7): relaxation [1/microSec]' + mess + mess + mess A = [a_mean] + mess A_err_in = [in_a_mean] + mess A_err_out = [out_a_mean] + mess + mess Freq. = [f_mean] + mess F_err_in = [in_f_mean] + mess F_err_out = [out_f_mean] + mess +* mess Relax. = [r_mean] +* mess r_err_in = [in_r_mean] +* mess r_err_out = [out_r_mean] +* mess +* if ( [opt] .eq. 'strexp' ) then +* mess +* mess Beta = [beta_mean] +* mess Beta_err_in = [in_beta_mean] +* mess Beta_err_out = [out_beta_mean] +* mess +* endif + +* restore previous settings: + + cd [current] + zone + title ' ' + set ysiz + set xsiz + diff --git a/geant3/paw/fit_prec_ge.kumac b/geant3/paw/fit_prec_ge.kumac new file mode 100644 index 0000000..d383266 --- /dev/null +++ b/geant3/paw/fit_prec_ge.kumac @@ -0,0 +1,415 @@ +if [1].EQ.'?' .OR. [1].EQ.' ' then + mess '===============================================================================' + mess 'FIT_PREC IDLeft B fromFit toFit t0' + mess + mess 'Kumac to fit the GEANT simulated decay spectra of L,T,R,B.' + mess 'If IDleft is given, then it is assumed that' + mess ' IDtop = IDleft+2,' + mess ' IDright = IDleft+3,' + mess ' IDbottom = IDleft+4.' + mess + mess 'B : magnetic field in gauss' + mess 'fromFit : lower boundary for fit' + mess 'toFit : upper boundary for fit' + mess 't0 : time zero' + mess + mess '===============================================================================' + exitm +endif +* +*=========================================================================== +current = $hcdir() + +fitopt = 'bl0' | fit option b=parameters with boundaries + | l=log likelihood fit (Poisson stat.) + | 0=do not plot fit result + +* define Histo ID's + + idl = [1] + idt = [idl]+2 + idr = [idl]+3 + idb = [idl]+4 +* +* check if histograms exist in PAW memory +* + cd //pawc + + if ( $hexist([idl]) .eq. 0 ) then + mess + mess Did not found Left histogram ID [idl] + mess + exitm + endif + if ( $hexist([idt]) .eq. 0 ) then + mess + mess Did not found Top histogram ID [idt] + mess + exitm + endif + if ( $hexist([idr]) .eq. 0 ) then + mess + mess Did not found Right histogram ID [idr] + mess + exitm + endif + if ( $hexist([idb]) .eq. 0 ) then + mess + mess Did not found Bottom histogram ID [idb] + mess + exitm + endif + + shift + +* transfer magnetic field strength: + + B = [1] + +* get region of histo to be used for fit: + + if ($index([2],'.').NE.0) then + von = $rsigma([2]) + else + von = [2] + endif + if ($index([3],'.').NE.0) then + bis = $rsigma([3]) + else + bis = [3] + endif + +* get rebinning factor and fit options: + + t0 = [4] + + if ( [t0] .eq. ' ' .or. [t0] .eq. '!' ) then + mess + mess Missing time zero ! + mess + exitm + endif + if ( $vexist(tzero) ) then + v/de tzero + endif + v/cre tzero(1) + v/inp tzero(1) [t0] + +* check if added histos are to be used and read in histos: + +* do the fits: + + freq = [B] * 0.013554 | precession frequency in MHz + freqmin = 0.8*[freq] + freqmax = 1.2*[freq] + freqstep = 0.005*[freq] +* + if ( $vexist(par) ) then + v/de par + endif + if ( $vexist(pmin) ) then + v/de pmin + endif + if ( $vexist(pmax) ) then + v/de pmax + endif + if ( $vexist(step) ) then + v/de step + endif + if ( $vexist(errpar) ) then + v/de errpar + endif +* +* +* +* BKG N0 Asym mu+ freq. Phase tau_mu lambda +* ======================================================= +* +* v/cre par(7) r 0. 500. 0.25 0. 0. 2.19703 0. +* v/cre pmin(7) r -10. 0. 0. 0. -400. 0. 0. +* v/cre pmax(7) r 100. 10000. 0.60 0. +400. 3. 100. +* v/cre step(7) r 0. 10. 0.01 0. 5. 0. 0. +* v/cre errpar(7) r 1. 1. 1. 1. 1. 1. 1. + v/cre par(6) r 0. 500. 0.25 [freq] 90. 2.19703 + v/cre pmin(6) r -10. 0. 0. [freqmin] -400. 0. + v/cre pmax(6) r 100. 10000. 0.40 [freqmax] +400. 3. + v/cre step(6) r 0. 10. 0.01 [freqstep] 5. 0. + v/cre errpar(6) r 0. 1. 1. 1. 1. 0. + +* + tdec0 = tzero(1) | vector tzero filled in FIT_PREC.KUMAC + ttt = Fit from [von] to [bis], t0 = [tdec0] +* title [ttt] +* +* + hi/fit [idl]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar +* +* store fit result and error for A_Left, frequency and relaxation rate +* + a_l = par(3) + a_lerr = errpar(3) + a_l_weight = [a_l]/([a_lerr]*[a_lerr]) + inv_a_lerr = 1./([a_lerr]*[a_lerr]) + f_l = par(4) + f_lerr = errpar(4) + f_l_weight = [f_l]/([f_lerr]*[f_lerr]) + inv_f_lerr = 1./([f_lerr]*[f_lerr]) +* r_l = par(7) +* r_lerr = errpar(7) +* r_l_weight = [r_l]/([r_lerr]*[r_lerr]) +* inv_r_lerr = 1./([r_lerr]*[r_lerr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_l = par(8) +* beta_lerr = errpar(8) +* beta_l_weight = [beta_l]/([beta_lerr]*[beta_lerr]) +* inv_beta_lerr = 1./([beta_lerr]*[beta_lerr]) +* endif +* +* +* initialize the phase for the top detector +* + phase = par(5) + phase = [phase]+90. + v/inp par(5) [phase] + hi/fit [idt]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar +* +* store fit result and error for A_top +* + a_t = par(3) + a_terr = errpar(3) + a_t_weight = [a_t]/([a_terr]*[a_terr]) + inv_a_terr = 1./([a_terr]*[a_terr]) + f_t = par(4) + f_terr = errpar(4) + f_t_weight = [f_t]/([f_terr]*[f_terr]) + inv_f_terr = 1./([f_terr]*[f_terr]) +* r_t = par(7) +* r_terr = errpar(7) +* r_t_weight = [r_t]/([r_terr]*[r_terr]) +* inv_r_terr = 1./([r_terr]*[r_terr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_t = par(8) +* beta_terr = errpar(8) +* beta_t_weight = [beta_t]/([beta_terr]*[beta_terr]) +* inv_beta_terr = 1./([beta_terr]*[beta_terr]) +* endif +* +* +* +* initialize the phase for the right detector +* + phase = par(5) + phase = [phase]+90. + v/inp par(5) [phase] + hi/fit [idr]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar + +* store fit result and error for A_right +* + a_r = par(3) + a_rerr = errpar(3) + a_r_weight = [a_r]/([a_rerr]*[a_rerr]) + inv_a_rerr = 1./([a_rerr]*[a_rerr]) + f_r = par(4) + f_rerr = errpar(4) + f_r_weight = [f_r]/([f_rerr]*[f_rerr]) + inv_f_rerr = 1./([f_rerr]*[f_rerr]) +* r_r = par(7) +* r_rerr = errpar(7) +* r_r_weight = [r_r]/([r_rerr]*[r_rerr]) +* inv_r_rerr = 1./([r_rerr]*[r_rerr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_r = par(8) +* beta_rerr = errpar(8) +* beta_r_weight = [beta_r]/([beta_rerr]*[beta_rerr]) +* inv_beta_rerr = 1./([beta_rerr]*[beta_rerr]) +* endif +* +* +* +* initialize the phase for the bottom detector +* + phase = par(5) + phase = [phase]+90. + v/inp par(5) [phase] + hi/fit [idb]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar +* +* store fit result and error for A_bottom +* + a_b = par(3) + a_berr = errpar(3) + a_b_weight = [a_b]/([a_berr]*[a_berr]) + inv_a_berr = 1./([a_berr]*[a_berr]) + f_b = par(4) + f_berr = errpar(4) + f_b_weight = [f_b]/([f_berr]*[f_berr]) + inv_f_berr = 1./([f_berr]*[f_berr]) +* r_b = par(7) +* r_berr = errpar(7) +* r_b_weight = [r_b]/([r_berr]*[r_berr]) +* inv_r_berr = 1./([r_berr]*[r_berr]) +* +* if ( [opt] .eq. 'strexp' ) then +* beta_b = par(8) +* beta_berr = errpar(8) +* beta_b_weight = [beta_b]/([beta_berr]*[beta_berr]) +* inv_beta_berr = 1./([beta_berr]*[beta_berr]) +* endif +* +* +* plot result: + + pics kum + set ysiz 25 + zone 2 3 + + hi/pl [idl] e + hi/pl [idt] e + hi/pl [idr] e + hi/pl [idb] e +* +*-------------------------------- +* +* calculate the weighted mean and the "inner" and "outer" standard deviation +* fo asymmetry +* + sum1 = [a_l_weight] + [a_t_weight] + [a_r_weight] + [a_b_weight] + sum2 = [inv_a_lerr] + [inv_a_terr] + [inv_a_rerr] + [inv_a_berr] + a_mean = [sum1]/[sum2] + in_a_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* now calculate the outer standard deviation +* + dev_l = ([a_mean]-[a_l])*([a_mean]-[a_l])*[inv_a_lerr] + dev_t = ([a_mean]-[a_t])*([a_mean]-[a_t])*[inv_a_terr] + dev_r = ([a_mean]-[a_r])*([a_mean]-[a_r])*[inv_a_rerr] + dev_b = ([a_mean]-[a_b])*([a_mean]-[a_b])*[inv_a_berr] + sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] + out_a_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +*---------------------------------- +* +* mean of frequency +* + sum1 = [f_l_weight] + [f_t_weight] + [f_r_weight] + [f_b_weight] + sum2 = [inv_f_lerr] + [inv_f_terr] + [inv_f_rerr] + [inv_f_berr] + f_mean = [sum1]/[sum2] + in_f_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* + dev_l = ([f_mean]-[f_l])*([f_mean]-[f_l])*[inv_f_lerr] + dev_t = ([f_mean]-[f_t])*([f_mean]-[f_t])*[inv_f_terr] + dev_r = ([f_mean]-[f_r])*([f_mean]-[f_r])*[inv_f_rerr] + dev_b = ([f_mean]-[f_b])*([f_mean]-[f_b])*[inv_f_berr] + sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] + out_f_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +*-------------------------------------- +* +* mean of relaxation +* +* sum1 = [r_l_weight] + [r_t_weight] + [r_r_weight] + [r_b_weight] +* sum2 = [inv_r_lerr] + [inv_r_terr] + [inv_r_rerr] + [inv_r_berr] +* r_mean = [sum1]/[sum2] +* in_r_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* +* dev_l = ([r_mean]-[r_l])*([r_mean]-[r_l])*[inv_r_lerr] +* dev_t = ([r_mean]-[r_t])*([r_mean]-[r_t])*[inv_r_terr] +* dev_r = ([r_mean]-[r_r])*([r_mean]-[r_r])*[inv_r_rerr] +* dev_b = ([r_mean]-[r_b])*([r_mean]-[r_b])*[inv_r_berr] +* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] +* out_r_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +* mean of stretched exponent +* +* +* if ( [opt] .eq. 'strexp' ) then +* sum1 = [beta_l_weight] + [beta_t_weight] + [beta_r_weight] + [beta_b_weight] +* sum2 = [inv_beta_lerr] + [inv_beta_terr] + [inv_beta_rerr] + [inv_beta_berr] +* beta_mean = [sum1]/[sum2] +* in_beta_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation +* +* no calculate the outer standard deviation +* +* dev_l = ([beta_mean]-[beta_l])*([beta_mean]-[beta_l])*[inv_beta_lerr] +* dev_t = ([beta_mean]-[beta_t])*([beta_mean]-[beta_t])*[inv_beta_terr] +* dev_r = ([beta_mean]-[beta_r])*([beta_mean]-[beta_r])*[inv_beta_rerr] +* dev_b = ([beta_mean]-[beta_b])*([beta_mean]-[beta_b])*[inv_beta_berr] +* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b] +* out_beta_mean = $sigma(sqrt(1./3.*[sum3]/[sum2])) +* +* endif +* +*---------------------------------------- +* +* output fit parameter definitions to text window: + + mess + mess 'par(1): Background [1]' + mess 'par(2): normalization [1]' + mess 'par(3): asymmetry [1]' + mess 'par(4): frequency [MHz]' + mess 'par(5): phase [degree]' + mess 'par(6): muon lifetime [microSec]' + mess 'par(7): relaxation [1/microSec]' + mess + mess + mess A = [a_mean] + mess A_err_in = [in_a_mean] + mess A_err_out = [out_a_mean] + mess + mess Freq. = [f_mean] + mess F_err_in = [in_f_mean] + mess F_err_out = [out_f_mean] + mess +* mess Relax. = [r_mean] +* mess r_err_in = [in_r_mean] +* mess r_err_out = [out_r_mean] +* mess +* if ( [opt] .eq. 'strexp' ) then +* mess +* mess Beta = [beta_mean] +* mess Beta_err_in = [in_beta_mean] +* mess Beta_err_out = [out_beta_mean] +* mess +* endif + +* output to graphics window (quick and dirty) + v/cre x(2) r 0 10 + v/cre y(2) r 0 10 + set xtic 0.0001 + set ytic 0.0001 + set vsiz 0.0001 + zone 1 3 3 s + graph 2 x y w +* + asym = $format([a_mean],F8.4) + in_err = $format([in_a_mean],F7.4) + out_err = $format([out_a_mean],F7.4) +* + text 1 8 'A = '//[asym]//'+-'//[in_err]//'+-'//[out_err] 0.4 +* + freq = $format([f_mean],F8.5) + in_err = $format([in_f_mean],F7.5) + out_err = $format([out_f_mean],F7.5) +* + app = ' MHz' + text 1 5 '[n] = '//[freq]//'+-'//[in_err]//'+-'//[out_err]//[app] 0.4 + + +* restore previous settings: + + cd [current] + zone +* title ' ' + set xtic + set ytic + set vsiz + set ysiz + set xsiz + diff --git a/geant3/paw/geant_dohist.kumac b/geant3/paw/geant_dohist.kumac new file mode 100644 index 0000000..c28218f --- /dev/null +++ b/geant3/paw/geant_dohist.kumac @@ -0,0 +1,158 @@ +* +* quicky KUMAC to get decay spectra from GEANT simulations +* (GEANT_LEMSR.KUMAC) +* +* TP, 03-Feb-1999, PSI +* +* [1] Geant Runnumber; for histogram titles +* +*---------------------------------------------------------------------- +* + if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then + mess + mess Missing GEANT RunNumber as parameter. + mess Needed for Histogram titles. + mess + mess Syntax: + mess + mess GEANT_DOHIST RunNo + mess + exitm + endif + runno = [1] +* +* cuts for e+ and e- and gammas +* + cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec. + cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec. + cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec. +* + cut $11 p0>30. + cut $10 desci>0.and.desco>0 + cut $1 de_left_i>0.and.de_left_o>0 + cut $2 de_top_i>0.and.de_top_o>0 + cut $3 de_rite_i>0.and.de_rite_o>0 + cut $4 de_bot_i>0.and.de_bot_o>0 +* +cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out +cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out +cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out +cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out +* + 1d 1 'GEANT '//[runno]//' volno, dE gt 0' 256 -.5 255.5 + 1d 2 'GEANT '//[runno]//' volno, dE gt 0, e+ bit or e- bit required' 256 -.5 255.5 + 1d 3 'GEANT '//[runno]//' volno, dE ge 0, only gammas' 256 -.5 255.5 + 1d 4 'GEANT '//[runno]//' volno, dE gt 0, only gammas' 256 -.5 255.5 + 1d 5 'GEANT '//[runno]//' volno, volno gt 0' 256 -.5 255.5 +*goto fine + 1d 404 'GEANT '//[runno]//' only left, de gt 0, 50ns bins' 220 0. 11000. + 1d 406 'GEANT '//[runno]//' only top, de gt 0, 50ns bins' 220 0. 11000. + 1d 407 'GEANT '//[runno]//' only right, de gt 0, 50ns bins' 220 0. 11000. + 1d 408 'GEANT '//[runno]//' only bottom, de gt 0, 50ns bins' 220 0. 11000. + 1d 504 'GEANT '//[runno]//' left , de gt 0, 50ns bins' 220 0. 11000. + 1d 506 'GEANT '//[runno]//' top , de gt 0, 50ns bins' 220 0. 11000. + 1d 507 'GEANT '//[runno]//' right, de gt 0, 50ns bins' 220 0. 11000. + 1d 508 'GEANT '//[runno]//' bottom, de gt 0, 50ns bins' 220 0. 11000. + 1d 604 'GEANT '//[runno]//' left , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 606 'GEANT '//[runno]//' top , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 607 'GEANT '//[runno]//' right , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 608 'GEANT '//[runno]//' bottom, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 704 'GEANT '//[runno]//' left , only gamma, de ge 0, 50ns bins' 220 0. 11000. + 1d 706 'GEANT '//[runno]//' top , only gamma, de ge 0, 50ns bins' 220 0. 11000. + 1d 707 'GEANT '//[runno]//' right , only gamma, de ge 0, 50ns bins' 220 0. 11000. + 1d 708 'GEANT '//[runno]//' bottom, only gamma, de ge 0, 50ns bins' 220 0. 11000. +* goto cont +* +fine: + 1d 1404 'GEANT '//[runno]//' only left, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1406 'GEANT '//[runno]//' only top, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1407 'GEANT '//[runno]//' only right, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1408 'GEANT '//[runno]//' only bottom, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1504 'GEANT '//[runno]//' left , de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1506 'GEANT '//[runno]//' top , de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1507 'GEANT '//[runno]//' right, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1508 'GEANT '//[runno]//' bottom, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1604 'GEANT '//[runno]//' left , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1606 'GEANT '//[runno]//' top , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1607 'GEANT '//[runno]//' right , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 + 1d 1608 'GEANT '//[runno]//' bottom, e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 +* + +mess Plot 1=Volno dE>0 +nt/pl 111.volno volno>0.and.$10 -1 + +mess Plot 2=Volno dE>0 and e+ or e- bit reuqired +nt/pl 111.volno volno>0.and.$10.and.($21.or.$22) -2 + +mess Plot 3=Volno dE>=0, only gammas +nt/pl 111.volno volno>0.and.partcode=17 -3 + +mess Plot 4=Volno dE>0, only gammas +nt/pl 111.volno volno>0.and.partcode=17.and.$10 -4 + +mess Plot 5=Volno volno>0 +nt/pl 111.volno volno>0 -5 + +mess NTPlot left=404 +nt/pl 111.tsci volno=17.and.$10 -404 +mess NTPlot top =406 +nt/pl 111.tsci volno=34.and.$10 -406 +mess NTPlot right=407 +nt/pl 111.tsci volno=68.and.$10 -407 +mess NTPlot bottom=408 +nt/pl 111.tsci volno=136.and.$10 -408 + +mess NTPlot left=1404 +nt/pl 111.tsci volno=17.and.$10 -1404 +mess NTPlot top =1406 +nt/pl 111.tsci volno=34.and.$10 -1406 +mess NTPlot right=1407 +nt/pl 111.tsci volno=68.and.$10 -1407 +mess NTPlot bottom=1408 +nt/pl 111.tsci volno=136.and.$10 -1408 + +mess NTPlot left id=504 +nt/pl 111.tsci $1.and.$5 -504 +mess NTPlot top id=506 +nt/pl 111.tsci $2.and.$6 -506 +mess NTPlot rite id=507 +nt/pl 111.tsci $3.and.$7 -507 +mess NTPlot bot id=508 +nt/pl 111.tsci $4.and.$8 -508 + +mess NTPlot left id=604 +nt/pl 111.tsci $1.and.$5.and.($21.or.$22) -604 +mess NTPlot top id=606 +nt/pl 111.tsci $2.and.$6.and.($21.or.$22) -606 +mess NTPlot rite id=607 +nt/pl 111.tsci $3.and.$7.and.($21.or.$22) -607 +mess NTPlot bot id=608 +nt/pl 111.tsci $4.and.$8.and.($21.or.$22) -608 + +mess NTPlot left id=704 +nt/pl 111.tsci partcode=17.and.$5 -704 +mess NTPlot top id=706 +nt/pl 111.tsci partcode=17.and.$6 -706 +mess NTPlot rite id=707 +nt/pl 111.tsci partcode=17.and.$7 -707 +mess NTPlot bot id=708 +nt/pl 111.tsci partcode=17.and.$8 -708 + +mess NTPlot left id=1504 +nt/pl 111.tsci $1.and.$5 -1504 +mess NTPlot top id=1506 +nt/pl 111.tsci $2.and.$6 -1506 +mess NTPlot rite id=1507 +nt/pl 111.tsci $3.and.$7 -1507 +mess NTPlot bot id=1508 +nt/pl 111.tsci $4.and.$8 -1508 + +mess NTPlot left id=1604 +nt/pl 111.tsci $1.and.$5.and.($21.or.$22) -1604 +mess NTPlot top id=1606 +nt/pl 111.tsci $2.and.$6.and.($21.or.$22) -1606 +mess NTPlot rite id=1607 +nt/pl 111.tsci $3.and.$7.and.($21.or.$22) -1607 +mess NTPlot bot id=1608 +nt/pl 111.tsci $4.and.$8.and.($21.or.$22) -1608 + diff --git a/geant3/paw/geant_dohist1.kumac b/geant3/paw/geant_dohist1.kumac new file mode 100644 index 0000000..71db0b8 --- /dev/null +++ b/geant3/paw/geant_dohist1.kumac @@ -0,0 +1,58 @@ +* +* quicky KUMAC to get decay spectra from GEANT simulations +* (GEANT_LEMSR.KUMAC) +* +* TP, 03-Feb-1999, PSI +* +* [1] Geant Runnumber; for histogram titles +* +* +* decay spectra, only 1 detector has hit, e+ e- bit required +*---------------------------------------------------------------------- +* + if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then + mess + mess Missing GEANT RunNumber as parameter. + mess Needed for Histogram titles. + mess + mess Syntax: + mess + mess GEANT_DOHIST RunNo + mess + exitm + endif + runno = [1] +* +* cuts for e+ and e- and gammas +* + cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec. + cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec. + cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec. +* + cut $11 p0>30. + cut $10 desci>0.and.desco>0 + cut $1 de_left_i>0.and.de_left_o>0 + cut $2 de_top_i>0.and.de_top_o>0 + cut $3 de_rite_i>0.and.de_rite_o>0 + cut $4 de_bot_i>0.and.de_bot_o>0 +* +cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out +cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out +cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out +cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out +* + 1d 614 'GEANT '//[runno]//' only left , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 616 'GEANT '//[runno]//' only top , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 617 'GEANT '//[runno]//' only right , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + 1d 618 'GEANT '//[runno]//' only bottom, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + + +mess NTPlot left id=614 +nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and.volno.eq.17 -614 +mess NTPlot top id=616 +nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and.volno.eq.34 -616 +mess NTPlot rite id=617 +nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and.volno.eq.68 -617 +mess NTPlot bot id=618 +nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and.volno.eq.136 -618 + diff --git a/geant3/paw/geant_dohist2.kumac b/geant3/paw/geant_dohist2.kumac new file mode 100644 index 0000000..ff01d71 --- /dev/null +++ b/geant3/paw/geant_dohist2.kumac @@ -0,0 +1,72 @@ +* +* quicky KUMAC to get decay spectra from GEANT simulations +* (GEANT_LEMSR.KUMAC) +* +* TP, 03-Feb-1999, PSI +* +* [1] Geant Runnumber; for histogram titles +* +* +* decay spectra, only 1 detector has hit, e+ e- bit required +*---------------------------------------------------------------------- +* + if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then + mess + mess Missing GEANT RunNumber as parameter. + mess Needed for Histogram titles. + mess + mess Syntax: + mess + mess GEANT_DOHIST RunNo + mess + exitm + endif + runno = [1] +* +* cuts for e+ and e- and gammas +* + cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec. + cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec. + cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec. +* + cut $11 p0>30. + cut $10 desci>0.and.desco>0 + cut $1 de_left_i>0.and.de_left_o>0 + cut $2 de_top_i>0.and.de_top_o>0 + cut $3 de_rite_i>0.and.de_rite_o>0 + cut $4 de_bot_i>0.and.de_bot_o>0 +* +cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out +cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out +cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out +cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out +* +1d 624 'GEANT '//[runno]//' left, dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. +1d 626 'GEANT '//[runno]//' top , dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. +1d 627 'GEANT '//[runno]//' right, dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. +1d 628 'GEANT '//[runno]//' bottom, dE others=0, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000. + +1d 1624 'GEANT '//[runno]//' left, dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 +1d 1626 'GEANT '//[runno]//' top , dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 +1d 1627 'GEANT '//[runno]//' right, dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 +1d 1628 'GEANT '//[runno]//' bottom, dE others=0, e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5 + + +mess NTPlot left id=624 +nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and..not.($2.or.$3.or.$4) -624 +mess NTPlot top id=626 +nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and..not.($1.or.$3.or.$4) -626 +mess NTPlot rite id=627 +nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and..not.($1.or.$2.or.$4) -627 +mess NTPlot bot id=628 +nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and..not.($1.or.$2.or.$3) -628 + +mess NTPlot left id=1624 +nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and..not.($2.or.$3.or.$4) -1624 +mess NTPlot top id=1626 +nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and..not.($1.or.$3.or.$4) -1626 +mess NTPlot rite id=1627 +nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and..not.($1.or.$2.or.$4) -1627 +mess NTPlot bot id=1628 +nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and..not.($1.or.$2.or.$3) -1628 + diff --git a/geant3/paw/geant_eloss.kumac b/geant3/paw/geant_eloss.kumac new file mode 100644 index 0000000..c142b5d --- /dev/null +++ b/geant3/paw/geant_eloss.kumac @@ -0,0 +1,92 @@ +* +* quicky KUMAC to get energy loss spectra from GEANT simulations +* (GEANT_ELOSS.KUMAC) +* +* TP, 03-Feb-1999, PSI +* +* [1] Geant Runnumber; for histogram titles +* +*---------------------------------------------------------------------- +* + if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then + mess + mess Missing GEANT RunNumber as parameter. + mess Needed for Histogram titles. + mess + mess Syntax: + mess + mess GEANT_DOHIST RunNo + mess + exitm + endif + runno = [1] +* +* cuts for e+ and e- and gammas +* + cut $20 btest(int(partcode),0).and.btest(int(partcode),4) | bits for gamma's + cut $21 btest(int(partcode),1).and.btest(int(partcode),5) | bits for e+ + cut $22 btest(int(partcode),2).and.btest(int(partcode),6) | bits for e- +* cut $30 partcode=17 | only gammas +* +* cut $11 p0>30. +* cut $10 desci>0.and.desco>0 +* cut $1 de_left_i>0.and.de_left_o>0 +* cut $2 de_top_i>0.and.de_top_o>0 +* cut $3 de_rite_i>0.and.de_rite_o>0 +* cut $4 de_bot_i>0.and.de_bot_o>0 +* +* cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out +* cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out +* cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out +* cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out +* + 1d 100 'GEANT '//[runno]//' desci gt 0' 200 0. 10. + 1d 101 'GEANT '//[runno]//' descie+ gt 0' 200 0. 10. + 1d 102 'GEANT '//[runno]//' descie- gt 0' 200 0. 10. + 1d 103 'GEANT '//[runno]//' desciga gt 0' 200 0. 1. + 1d 110 'GEANT '//[runno]//' desci gt 0, require e+ or e- bits' 200 0. 10. + + 1d 200 'GEANT '//[runno]//' desco gt 0' 200 0. 10. + 1d 201 'GEANT '//[runno]//' descoe+ gt 0' 200 0. 10. + 1d 202 'GEANT '//[runno]//' descoe- gt 0' 200 0. 10. + 1d 203 'GEANT '//[runno]//' descoga gt 0' 200 0. 1. + 1d 210 'GEANT '//[runno]//' desco gt 0, require e+ or e- bits' 200 0. 10. +* +* inner detectors +* + +mess Plot 100=Desci gt 0 +nt/pl 111.desci desci>0 -100 + +mess Plot 101=Descie+ gt 0 +nt/pl 111.descipos descipos>0 -101 + +mess Plot 102=Descie- gt 0 +nt/pl 111.desciele desciele>0 -102 + +mess Plot 103=Descigamma gt 0 +nt/pl 111.descigam descigam>0 -103 + +mess Plot 110=Desci gt 0, e+ or e- bits required +nt/pl 111.desci desci>0.and.($21.or.$22) -110 + +* +* outer detectors +* + +mess Plot 200=Desco gt 0 +nt/pl 111.desco desco>0 -200 + +mess Plot 201=Descoe+ gt 0 +nt/pl 111.descopos descopos>0 -201 + +mess Plot 202=Descoe- gt 0 +nt/pl 111.descoele descoele>0 -202 + +mess Plot 203=Descogamma gt 0 +nt/pl 111.descogam descogam>0 -203 + +mess Plot 210=Desco gt 0, e+ or e- bits required +nt/pl 111.desco desco>0.and.($21.or.$22) -210 + + diff --git a/geant3/paw/geant_part.kumac b/geant3/paw/geant_part.kumac new file mode 100644 index 0000000..1b5c5f3 --- /dev/null +++ b/geant3/paw/geant_part.kumac @@ -0,0 +1,196 @@ +* +* paw$dir:geant_part.kumac +* +* KUMAC to copy first 400 ns of GEANT decay spectra to +* histogram. +* +* [1] ID of left histogram; other ID's will be reconstructed +* +* TP, 23-feb-1999 +* +*---------------------------------------------------------------------------- +* + if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!' ) then + mess + mess GEANT_PART IDleft + mess + mess Copies first 400ns of existing GEANT decay spectra + mess to new IDs; NewID = OldID+1000 + mess + exitm + endif + idl = [1] + idt = [idl]+2 + idr = [idt]+1 + idb = [idr]+1 + idln = [idl]+1000 + idtn = [idt]+1000 + idrn = [idr]+1000 + idbn = [idb]+1000 +* + upp = 400. | upper bin center of new histograms +* + if ( $vexist(data) ) then + v/de data + endif +* +* get histo information +* + if ( $hexist([idl]) .eq. 0 ) then + mess + mess Left Histogram ID = [idl] does not exist ! + mess + exitm + endif + if ( $hexist([idt]) .eq. 0 ) then + mess + mess Top Histogram ID = [idt] does not exist ! + mess + exitm + endif + if ( $hexist([idr]) .eq. 0 ) then + mess + mess Right Histogram ID = [idr] does not exist ! + mess + exitm + endif + if ( $hexist([idb]) .eq. 0 ) then + mess + mess Bottom Histogram ID = [idb] does not exist ! + mess + exitm + endif +* +*---------------------------------------------------------------------------- +* +* L E F T +* + xmin = $hinfo([idl], 'XMIN') + xmax = $hinfo([idl], 'XMAX') + nbin = $hinfo([idl], 'XBINS') + bin = ([xmax] - [xmin])/[nbin] +* + mess Left ID = [idl] with bin size [bin] +* +* copy histogram to vector +* + v/cre data([nbin] + get_vec/con [idl] data +* +* construct new upper limit +* + nbin = [upp]/[bin] + 1 + xmax = [xmin] + [nbin]*[bin] + mess New LeftId = [idln] with xmax = [xmax] and nbin = [nbin] +* +* book histogram +* + if ( $hexist([idln]) ) then + hi/de [idln] + endif + 1d [idln] 'part of left, id = '//[idl] [nbin] [xmin] [xmax] + put_vec/con [idln] data(1:[nbin]) + v/de data +* +*--------------------------------------------------------------------------- +* +* T O P +* + xmin = $hinfo([idt], 'XMIN') + xmax = $hinfo([idt], 'XMAX') + nbin = $hinfo([idt], 'XBINS') + bin = ([xmax] - [xmin])/[nbin] +* + mess Top ID = [idt] with bin size [bin] +* +* copy histogram to vector +* + v/cre data([nbin] + get_vec/con [idt] data +* +* construct new upper limit +* + nbin = [upp]/[bin] + 1 + xmax = [xmin] + [nbin]*[bin] + mess New TopId = [idtn] with xmax = [xmax] and nbin = [nbin] +* +* book histogram +* + if ( $hexist([idtn]) ) then + hi/de [idtn] + endif + 1d [idtn] 'part of top, id = '//[idt] [nbin] [xmin] [xmax] + put_vec/con [idtn] data(1:[nbin]) + v/de data +* +*--------------------------------------------------------------------------- +* +* Right +* + xmin = $hinfo([idr], 'XMIN') + xmax = $hinfo([idr], 'XMAX') + nbin = $hinfo([idr], 'XBINS') + bin = ([xmax] - [xmin])/[nbin] +* + mess Right ID = [idr] with bin size [bin] +* +* copy histogram to vector +* + v/cre data([nbin] + get_vec/con [idr] data +* +* construct new upper limit +* + nbin = [upp]/[bin] + 1 + xmax = [xmin] + [nbin]*[bin] + mess New RightId = [idrn] with xmax = [xmax] and nbin = [nbin] +* +* book histogram +* + if ( $hexist([idrn]) ) then + hi/de [idrn] + endif + 1d [idrn] 'part of Right, id = '//[idr] [nbin] [xmin] [xmax] + put_vec/con [idrn] data(1:[nbin]) + v/de data +* +*------------------------------------------------------------------------------ +* +* Bottom +* + xmin = $hinfo([idb], 'XMIN') + xmax = $hinfo([idb], 'XMAX') + nbin = $hinfo([idb], 'XBINS') + bin = ([xmax] - [xmin])/[nbin] +* + mess Bottom ID = [idb] with bin size [bin] +* +* copy histogram to vector +* + v/cre data([nbin] + get_vec/con [idb] data +* +* construct new upper limit +* + nbin = [upp]/[bin] + 1 + xmax = [xmin] + [nbin]*[bin] + mess New BottomId = [idbn] with xmax = [xmax] and nbin = [nbin] +* +* book histogram +* + if ( $hexist([idbn]) ) then + hi/de [idbn] + endif + 1d [idbn] 'part of Bottom, id = '//[idb] [nbin] [xmin] [xmax] + put_vec/con [idbn] data(1:[nbin]) + v/de data +* +*------------------------------------------------------------------------------ +* + zone 2 2 + hi/pl [idln] + hi/pl [idtn] + hi/pl [idrn] + hi/pl [idbn] + zone +* diff --git a/geant3/paw/readge.kumac b/geant3/paw/readge.kumac new file mode 100644 index 0000000..1f4ebf2 --- /dev/null +++ b/geant3/paw/readge.kumac @@ -0,0 +1,62 @@ +* +* READGE.KUMAC +* +* KUMAC to open NT file from GEANT_LEMSR simulation on next free +* LUN. +* +* Input parameter: +* +* [1]: GEant Runnumber +* +* TP, 10-feb-1999, PSI +* +* TP, 06-Apr-2001, PSI: Unix, Linux and Windows version +* +*--------------------------------------------------------------------------- +* + if ([1] .eq. '?' .or. [1] .eq. ' ' .or. [1] .eq. '!') then + mess + mess READGE kumac to open GEANT_LEMSR NTUPLE file. + mess + mess Syntax: + mess + mess READGE runNr + mess + exitm + endif + runnr = [1] +* +* determine operating system to build the filename; +* only Windows or Unix/Linux +* + os = $OS +* + if ( [os] .eq. 'Windows') then + file = 'c:\users\thomas\mc\geant\data\geant_lemsr_'//[runnr]//'.nt' + goto check_file + endif + if ( [os] .eq. 'UNIX') then +* if ( $OSTYPE .eq. 'Linux' ) then + file = /scratch/$USER/geant/geant_lemsr_[runnr].nt + goto check_file +* else +* file = $HOME//'/mc/geant/data/geant_lemsr_'//[runnr]//'.nt' +* goto check_file +* endif + endif +* + mess Operating system [os] not supported + exitm +* +check_file: + + if ( $fexist([file]) ) then + hi/file 0 [file] 4096 + else + mess + mess File [file] does not exist... + mess + endif +* +* + diff --git a/geant3/src/lemsr/Makefile b/geant3/src/lemsr/Makefile new file mode 100644 index 0000000..7ec9b05 --- /dev/null +++ b/geant3/src/lemsr/Makefile @@ -0,0 +1,50 @@ +# +# $HOME/mc/geant/src/lemsr/Makefile +# +# makefile to compile and link the program geant_lemsr to +# simulate decaying muons at the LE-muSR spectrometer at PSI +# NTuple file to be analysed by PAW. +# +# TP, 15-Sep-2000, PSI +# 09-Apr-2001 +# +#-------------------------------------------------------------------------- +# +# The following lines contain specific switches for different +# Unix systems. Outcomment the unwanted. +# +#-------------------------------------------------------------------------- +# +# For Digital Unix +# +#FF = f77 +# +#-------------------------------------------------------------------------- +# +# For Linux +# +FF = g77 +CFLAGS = -c -I$(CERN_ROOT)/include/ -Wno-globals -finit-local-zero +#CFLAGS += -fno-automatic # not really necessary +# +#-------------------------------------------------------------------------- +# +OBJ = geant_lemsr_main.o get_direction.o gudcay.o gukine.o \ + guout.o gustep.o gutrack.o gutreve.o michel.o rotate.o ugeom.o +.F.o: + $(FF) $(CFLAGS) $*.F + +geant_lemsr: $(OBJ) + $(FF) -o $(HOME)/mc/geant/bin/geant_lemsr $(OBJ) \ + `cernlib geant321 pawlib graflib packlib mathlib` +# `cernlib geant321 pawlib graflib/Motif packlib mathlib` +# +# the program cernlib searches the cernlib libraries as well as +# needed system libraries +# +clean: + rm -f *.o *~ \#* .#* +#-------------------------------------------------------------------------- +# +# EOF +# diff --git a/geant3/src/lemsr/common.inc b/geant3/src/lemsr/common.inc new file mode 100644 index 0000000..e04f079 --- /dev/null +++ b/geant3/src/lemsr/common.inc @@ -0,0 +1,124 @@ +c +c ~/mc/geant/lemsr/src/common.inc +c +c GEANT_LEMSR program, +c include common blocks other than NTuple variables +c +c TP, 05-Apr-2001, PSI +c +c 17-Dec-2002 TP extended bfield array from 4 to 5, +c bfield(5)=polarisation in sample +c +c----------------------------------------------------------------------- +c +c check OS +c +#if defined( _WIN32 ) +#define OS_WIN +#else +#define OS_UNIX +#endif +c +c parameter +c + character*(*) IN_DIR, OUT_DIR +#if defined(OS_WIN) + parameter IN_DIR = 'c:\users\thomas\mc\geant\inp\' + parameter OUT_DIR = 'c:\users\thomas\mc\geant\data\' +#else + parameter (IN_DIR = '/afs/psi.ch/user/p/prokscha/mc/geant/inp/') +c parameter (OUT_DIR = '/afs/psi.ch/user/p/prokscha/mc/geant/data/') + parameter (OUT_DIR = '/scratch/prokscha/geant/') +#endif +c + character*(*) INFILE, PARFILE + parameter (INFILE = IN_DIR//'geant_lemsr_input.dat') + parameter (PARFILE = IN_DIR//'geant_lemsr.input') +c +c +c INFILE: read input data from file instead of throwing them +c in the program +c +c PARFILE: read input parameter from this file +c + real*4 MU_MASS, TWO_PI, POS_MASS, SQR_POS_MASS, + 1 MAX_TOTAL_ENERGY, TAU_MUON, MU_GYRO, R_MCP + parameter (MU_MASS = .105658389 ) !muon mass, GeV/c^2 + parameter (POS_MASS = 0.51099906 ) ! e+ mass, MeV/c^2 + parameter (SQR_POS_MASS = 0.26112 ) ! squared positron mass + parameter (MAX_TOTAL_ENERGY = 52.8304302) ! MeV +c + parameter (TAU_MUON = 2.19703E-6 ) ! muon life time in s + parameter (MU_GYRO = 13554.) ! Hz/G +c + parameter (R_MCP = 1.00 ) ! radius of MCP in cm + parameter (TWO_PI = 6.283185307) +c +c------------------------------------------------------------------- +c + real*4 bfield(5) ! new data card with B-field components +c ! bfield(4) = polarization (%) +c ! bfield(5) = polarization (%) in sample + integer*4 ix1 !random generator seed +c +c count successfull and rejected MICHEL-energy throws +c + integer*4 cou, rejcou +c +c------------------------------------------------------------------------------ +c +c flags for geometrical setup; it is possible to mount or remove the +c +c 'mcp2' MCP2 plates +c 'mcpa' MCP2 anode +c 'mcps' MCP2 stainless steel vacuum tube +c 'samp' sample (this is at the moment the sample holder) +c 'cryo' material of the cryostat without shielding +c 'crsh' cryo shield +c 'efla' the 100 CF flange at the end of the sample/MCP2 tube +c +c the volumes are activated when printing +c +c geom 'mcp2' 'efla' for example when data cards are requested +c + logical*1 l_mcp2 + logical*1 l_mcpa + logical*1 l_mcps + logical*1 l_samp + logical*1 l_cryo + logical*1 l_crsh + logical*1 l_efla + logical*1 l_run11 + character*4 samp_mat +c +c------------------------------------------------------------------------ +c +c flags for secondary particles +c + integer*4 lspar(10) ! new data card defined by CALL FFKEY to +c ! select seoncdary particles +c ! ' e+', ' e-' and ' gam' possible. +c + logical*1 l_pos + logical*1 l_ele + logical*1 l_gam +c +c------------------------------------------------------------------------- +c +c flags for output of percentage of processed events +c + integer*4 percent(9) +c +c------------------------------------------------------------------------- +c +c--- the common blocks +c + common /magnetic_field/ bfield + common /seeds/ ix1 + common /counter/ cou, rejcou + common /geo_flags/ l_mcp2, l_mcpa, l_mcps, l_samp, l_cryo, + 1 l_crsh, l_efla, l_run11, samp_mat + common /part_flag/ l_pos, l_ele, l_gam + common /proc_event/ percent + +c----------------------------------------------------------------------- diff --git a/geant3/src/lemsr/cwn.inc b/geant3/src/lemsr/cwn.inc new file mode 100644 index 0000000..1845c43 --- /dev/null +++ b/geant3/src/lemsr/cwn.inc @@ -0,0 +1,45 @@ +c +c cwn.inc +c +c include variable declaration for cwn Ntuples +c +c TP, 15-Sep-2000, PSI +c +c----------------------------------------------------------------------------- +c +c CWN variables +c + real*4 nt_x0, nt_y0, nt_z0, nt_px0, nt_py0, nt_pz0, nt_p0 + integer*4 nt_ipart + real*4 nt_xe, nt_ye, nt_ze, nt_eposend + real*4 nt_tmcp, nt_tsci ! times are integer, use 1ns bins + real*4 nt_xsci, nt_ysci, nt_zsci + real*4 nt_desci, nt_descipos, nt_descigam, nt_desciele + real*4 nt_xsco, nt_ysco, nt_zsco + real*4 nt_desco, nt_descopos, nt_descogam, nt_descoele + integer*4 nt_volno, nt_partcode + integer*4 nt_mcpnhits, nt_scinhits, nt_sconhits + real*4 nt_pxpos, nt_pypos, nt_pzpos + real*4 nt_demcp + real*4 nt_descil, nt_descit, nt_descir, nt_descib, + 1 nt_descol, nt_descot, nt_descor, nt_descob +c +c the common blocks +c + common /de_ltrb/ nt_descil, nt_descit, nt_descir, nt_descib, + 1 nt_descol, nt_descot, nt_descor, nt_descob + common /init_par/ nt_x0, nt_y0, nt_z0, nt_px0, nt_py0, nt_pz0, + 1 nt_p0 + common /partcode/ nt_ipart + common /end_pos/ nt_xe, nt_ye, nt_ze, nt_eposend + common /times/ nt_tmcp, nt_tsci + common /sci_pos/ nt_xsci, nt_ysci, nt_zsci + common /sci_elos/ nt_desci, nt_descipos, nt_descigam, nt_desciele + common /sco_pos/ nt_xsco, nt_ysco, nt_zsco + common /sco_elos/ nt_desco , nt_descopos, nt_descogam, nt_descoele + common /flags/ nt_volno, nt_partcode + common /nhits/ nt_mcpnhits, nt_scinhits, nt_sconhits + common /pos_mom/ nt_pxpos, nt_pypos, nt_pzpos + common /de_mcp/ nt_demcp +c +c----------------------------------------------------------------------------- diff --git a/geant3/src/lemsr/geant_lemsr_main.F b/geant3/src/lemsr/geant_lemsr_main.F new file mode 100644 index 0000000..faa5809 --- /dev/null +++ b/geant3/src/lemsr/geant_lemsr_main.F @@ -0,0 +1,816 @@ +#include "geant321/pilot.h" +c----------------------------------------------- +c ~/mc/geant/src/lemsr/geant_lemsr_main.F +c----------------------------------------------- +c + program geant_lemsr +c +c Program very much "grown", so many things could be programmed in +c a better way... +c +c program to simulate the response of the low energy muon muSR +c detector setup at PSI. Some parts are taken from Nai_simul.for. +c It is intended to steer the simulation by format free inputs +c (routine gffgo); this can be done interactively or data sets +c can be read from file. +c +c Uses Run11 setup geometry. +c +c For precession: at the moment use the total TOF for the calculation +c of the precession angle. When the muon is started as particle code=500 +c and upstream of the target, the precession is calculated as if the +c muon is flying all the time through the B-field. +c +c TP 18-Jan-1999 PSI +c +c TP, 4-Feb-1999: CW Ntuples +c TP, 9-Feb-1999: store energy loss per detector in NTuple +c TP,16-Feb-1999: time in NTUPLE as REAL*4 instead of INTEGER +c TP,26-Feb-1999: If particles start upstream gaussian divergence of +c beam now possible. +c +c***************************************************************************** +c +c TP, 11-May-2000 PSI +c Windows NT 4.0 version; uses include files which +c I created manually from the VMS GEANT321.TLB text library. +c +c TP, 15-Sep-2000 PSI +c Unix version Used fsplit to separate program into +c subroutines and functions. +c +c +c TP, 05-Apr-2001 PSI +c Unix, NT, and Linux (to be checked) +c version; use the Compaq Fortran fpp +c precompiler on NT +c +c 17-Dec-2002 TP extend bfield array to account for different +c polarization in sample (=bfield(5)) +c gaussian muon momentum distribution possible +c via sets(8) and sets(9): +c sigma_p = p0 * (sets(8)/100.+sets(9)/1000.) +c +c------------------------------------------------------------------ +c + implicit none +c +c include some files from $CERN_ROOT/include/ +c +c to avoid explicitly including the gt***.inc file +c define CERNLIB_TYPE; the gt***.inc files are then +c automatically included in the gc***.inc files. +c +#define CERNLIB_TYPE +#include "geant321/gclist.inc" +#include "geant321/gcflag.inc" +#include "geant321/gcphys.inc" +#include "geant321/gccuts.inc" +#include "geant321/gckine.inc" +c +#include "cwn.inc" +#include "common.inc" +c +c +c the particle definition (variable IPART): see GEANT manual CONS300-1 +c +c particle IPART +c +c gamma 1 +c positron 2 +c electron 3 +c neutrino 4 +c muon+ 5 +c muon- 6 +c pion0 7 +c pion+ 8 +c pion- 9 +c ... +c + integer*4 NG, NH + parameter (NG = 2000000) + parameter (NH = 1000000) + integer*4 ipawc,q + integer*4 iwktyp + integer*4 iquest(100) +c + common /quest/iquest !needed to get enough space for NTUPLE + common /pawc/ipawc(NH) + common /gcbank/q(NG) +c +c parameter for drawing (used in GEANT routine GDRAW) +c + real*4 gdtheta,gdphi,gdpsi + real*4 gdx0,gdy0,gdsx,gdsy +c + real*4 xn,yn,scal,r(2) + integer*4 i +c +c for opening NTUPLE file and NTUPLE +c + integer*4 istat + character*72 ntfile +c +c the run number copied from GEAN variable IDRUN +c + character*4 runno +c +c for graphic output to PostScript file +c + character*72 psfile +c + integer*4 lunin /4/ + integer*4 ipres + real*4 ub(1) +c +c----------------------------------------------------------------- +c +c flags for NTUPLE variable definitions +c +c 'init' write initial parameter to NT +c 'code' initial particle ID to NT +c 'endp' end position and energy of e+ +c 'time' t-mcp and t-sci +c 'scip' position when particle hits inner Sc. +c 'scie' energy loss of particles in inner sc. +c 'scop' position when particle hits outer Sc. +c 'scoe' energy loss of particles in outer sc. +c 'flag' Vol. numbers and particle codes; needed to get the single +c decay spectra +c 'nhit' Number of Hits in detectors MCP, SCI and SCO +c 'posm' Generated momentum of decay positron +c 'mcpe' energy loss in MCP detector +c 'allv' use all possible NT variables +c 'ltrb' energy loss in single detectors (inner and outer separate) +c + integer*4 ntvar(20) ! new data card defined by CALL FFKEY to +c ! select NTuple variables +c + logical*1 l_init_par /.false./ + logical*1 l_code /.false./ + logical*1 l_endp /.false./ + logical*1 l_time /.false./ + logical*1 l_scip /.false./ + logical*1 l_scie /.false./ + logical*1 l_scop /.false./ + logical*1 l_scoe /.false./ + logical*1 l_flag /.false./ + logical*1 l_nhit /.false./ + logical*1 l_posm /.false./ + logical*1 l_mcpe /.false./ + logical*1 l_allv /.false./ + logical*1 l_ltrb /.false./ +c +c------------------------------------------------------------------ +c +c start of program code +c +c------------------------------------------------------------------ +c +c general initializations +c +c initial random generator seeds +c + ix1 = 343156821 +c + do i = 1, 9 + percent(i) = 0 + enddo +c + l_mcp2 = .false. + l_mcpa = .false. + l_mcps = .false. + l_samp = .false. + l_cryo = .false. + l_crsh = .false. + l_efla = .false. + l_run11 = .false. + samp_mat = ' ' +c + l_pos = .false. + l_ele = .false. + l_gam = .false. +c +c----------------------------------------- +c + call gzebra(NG) + call hlimit(-NH) + call ginit ! initialize GEANT variables + call gzinit ! initialize GEANT data structure + call gpart ! GEANT standard particles + call gspart(500, 'muon+', 5, 0.1056584, 1., 2.19703E-6, ub, 1) +c +c------------------------------------------------------------------ +c +c open parfile in any case, the user can decide if data should +c be read from parfile or entered by terminal +c + open(lunin, file=PARFILE, status='OLD',err=1003) + write(6,*) ' ' + write(6,*) '%GEANT-LEMSR-I- Successfully opened data card file ' + write(6,*) + write(6,*) PARFILE +c +c read data cards +c + write(6,*) ' ' + write(6,*) 'Type in the data cards now... finish with ''STOP''' + write(6,*) '--------------------------------------------------' + write(6,*) ' ' + write(6,*) 'READ 4 : reads input data from '// + 1 'LUN=4--> file ' + write(6,*) ' ' + write(6,*) PARFILE + write(6,*) ' ' +c write(6,*) 'RUNG 24 : user run number' +c write(6,*) ' ' +c write(6,*) 'PLOT 0 : no graphical output' +c write(6,*) 'PLOT 1 : output on workstation only' +c write(6,*) 'PLOT 2 : output on workstation and metafile' +c write(6,*) ' ' +c write(6,*) 'SETS 2 0 : particle is e+(=2), michel '// +c 1 'distributed' +c write(6,*) ' start MCP2 detector' +c write(6,*) 'SETS 2 -74 27 5 1 : beam positron with 27.51 MeV/c '// +c 1 'at z=-74cm' +c write(6,*) 'SETS 5 -74 2 0 0 : beam muon(=5) with 2.00 MeV/c' +c write(6,*) ' particle started 74cm upstream of' +c write(6,*) ' of center of MCP2 tube.' +c write(6,*) ' ' +c write(6,*) 'TRIG 100 : 100 events to be generated' +c write(6,*) ' ' +c write(6,*) 'GEOM ''mcp2'' ''samp'' : mcp2 and cryo are mounted.' +c write(6,*) ' ' +c write(6,*) 'SPAR '' e+'' '' e-'' '' gam'' : select sec. '// +c 1 'particles.' +c write(6,*) ' ' +c write(6,*) 'BFIE 0. 0. 100. 90. : Sample B-field (Gauss), '// +c 1 'Bx, By, Bz,' +c write(6,*) ' and Polarization = 90 % here.' +c write(6,*) +c write(6,*) 'NTVA ''init'' ''time'' ''flag'' : define NT '// +c 1 'variables.' +c write(6,*) 'See geant_lemsr.input for more information.' +c write(6,*) +c write(6,*) 'View 1 Thet Phi Psi u0 v0 10*scal_x 10*scal_y : '// +c 1 'for graphical output' +c write(6,*) 'Def. 60 40 0 10 10 2 2' +c write(6,*) ' ' +c +c--------------------------------------------------------------- +c +c define new data card for particles +c + call ffkey('SPAR', lspar, 10, 'M') +c +c define new data card for magnetic field components +c + call ffkey('BFIE', bfield, 5, 'R') +c +c define new data card for NTUPLE variables +c + call ffkey('NTVA', ntvar, 20, 'M') +c +c read data card +c + call gffgo +c + close(lunin) +c + write(6,*) ' ' + write(6,*) '%GEANT-LEMSR-I- Closed file ', PARFILE + write(6,*) +c + write(6,*) 'B-field components Bx, By, Bz = ', bfield(1), + 1 bfield(2), bfield(3) + write(6,*) ' ' + write(6,*) 'Initial muon polarization = ', bfield(4),' %.' + write(6,*) +c + if ( bfield(1) .ne. 0. ) then + write(6,*) '******************************************' + write(6,*) ' ' + write(6,*) 'GEANT-LEMSR-E- B_x not yet implemented...' + write(6,*) 'Sorry, - - - - > E X I T' + write(6,*) ' ' + write(6,*) '******************************************' + call exit + endif + if ( bfield(2) .ne. 0. .and. bfield(3) .ne. 0. ) then + write(6,*) '******************************************' + write(6,*) ' ' + write(6,*) 'GEANT-LEMSR-E- B_y >0 and B_z >0 not yet '// + 1 'implemented...' + write(6,*) 'Sorry, - - - - > E X I T' + write(6,*) ' ' + write(6,*) '******************************************' + call exit + endif +c + if ( lsets(2) .ne. 0 .and. lsets(7) .ne. 1) then + write(6,*) ' ' + write(6,*) 'particles of type ',lsets(1),' start at z0 = ', + 1 lsets(2) + write(6,*) ' ' + else if ( lsets(7) .eq. 1 ) then + write(6,*) ' ' + write(6,*) 'particles of type ',lsets(1),' start at z0 '// + 1 'read from file' + write(6,*) ' ' + else + write(6,*) ' ' + write(6,*) 'particles of type ',lsets(1),' start at z0 = '// + 1 '1.40002 cm.' + write(6,*) ' ' + endif +c +c------------------------------------- +c +c look for particles to be created during tracking +c + ipres = 0 + call glook(' e+',lspar,10,ipres) + if ( ipres .ne. 0 ) then + l_pos = .true. + write(6,*) ' ' + write(6,*) 'secondary particle: e+...' + endif + call glook(' e-',lspar,10,ipres) + if ( ipres .ne. 0 ) then + l_ele = .true. + write(6,*) 'secondary particle: e-...' + endif + call glook(' gam',lspar,10,ipres) + if ( ipres .ne. 0 ) then + l_gam = .true. + write(6,*) 'secondary particle: gamma...' + endif +c +c------------------------------------------------------------ +c +c now look for geometrical setup variables +c + call glook('mcp2',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_mcp2 = .true. + write(6,*) 'Setup with: MCP2' + endif +c + call glook('mcpa',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_mcpa = .true. + write(6,*) 'Setup with: MCP2 anode and support ring' + endif +c + call glook('ru11',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_run11 = .true. + if ( l_mcp2 .or. l_mcpa ) then + write(6,*) 'Setup with: MCP2 Run 11 geometry, delay '// + 1 'line anode' + endif + else + if ( l_mcp2 .or. l_mcpa ) then + write(6,*) 'Setup with: MCP2 Run 10 geometry, WSZ anode' + endif + endif +c + call glook('mcps',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_mcps = .true. + write(6,*) 'Setup with: MCP2 stainless steel vacuum tube' + endif +c + call glook('samp',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_samp = .true. + write(6,*) 'Setup with: sample holder' + endif +c + call glook('saal',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + samp_mat = 'saal' + if ( l_samp ) then + write(6,*) 'Setup with: Al - sample holder' + endif + else + if ( l_samp ) then + write(6,*) 'Setup with: Cu - sample holder' + endif + endif + call glook('cryo',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_cryo = .true. + write(6,*) 'Setup with: cryostat' + endif +c + call glook('crsh',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_crsh = .true. + write(6,*) 'Setup with: He shield' + endif +c + call glook('efla',lgeom,20,ipres) + if ( ipres .ne. 0 ) then + l_efla = .true. + write(6,*) 'Setup with: 100 CF flange at the end' + endif +c +c---------------------------------------------------- +c +c look for NT variables +c + call glook('init', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_init_par = .true. + write(6,*) 'NT variable: initial position and momenta.' + endif +c + call glook('code', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_code = .true. + write(6,*) 'NT variable: Code of primary particle.' + endif +c + call glook('endp', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_endp = .true. + write(6,*) 'NT variable: End position and end energy of e+.' + endif +c + call glook('time', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_time = .true. + write(6,*) 'NT variable: t-MCP2 and t-SCI.' + endif + call glook('scip', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_scip = .true. + write(6,*) 'NT variable: Position of particle when '// + 1 'hitting inner Sc..' + endif +c + call glook('scie', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_scie = .true. + write(6,*) 'NT variable: Energy deposited in inner Sc.' + endif +c + call glook('scop', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_scop = .true. + write(6,*) 'NT variable: Position of particle when '// + 1 'hitting outer Sc..' + endif +c + call glook('scoe', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_scoe = .true. + write(6,*) 'NT variable: Energy deposited in outer Sc.' + endif +c + call glook('flag', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_flag = .true. + write(6,*) 'NT variable: Volume number and particle codes.' + write(6,*) 'Note: the Vol. No. is necessary to get '// + 1 'single decay spectra.' + endif +c + call glook('nhit', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_nhit = .true. + write(6,*) 'NT variable: Hits per Track in MCP, Sci Sco.' + endif +c + call glook('posm', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_posm = .true. + write(6,*) 'NT variable: inital e+ momentum.' + endif +c + call glook('mcpe', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_mcpe = .true. + write(6,*) 'NT variable: Energy deposited in MCP2.' + endif +c + call glook('allv', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_allv = .true. + write(6,*) 'NT variable: all variables.' + endif +c + call glook('ltrb', ntvar,20,ipres) + if ( ipres .ne. 0 ) then + l_ltrb = .true. + write(6,*) 'NT variable: energies deposited in single '// + 1 'scintill.' + endif +c +c---------------------------------------------------- +c + write(runno,'(I4)',err=1002) idrun ! copy run number to string +c + if ( lsets(2) .ne. 0 .and. lsets(7) .ne. 1) then + if ( lsets(3) .gt. 0 ) then + pkine(1) = float(lsets(3)) + float(lsets(4))/10. + + 1 float(lsets(5))/100. + else + pkine(1) = float(lsets(3)) - float(lsets(4))/10. - + 1 float(lsets(5))/100. + endif + write(6,'(a,1x,f10.4)') ' %GEANT-LEMSR-I- beam momentum '// + 1 'set to ', pkine(1) + endif + if ( lsets(6) .ne. 0 ) then + write(6,*) ' ' + write(6,*) 'Gaussian beam divergence set to (degree) ', + 1 lsets(6) + write(6,*) ' ' + endif +c +c define volumes and detectors +c + call ugeom +c +c compute the cross-sections for the different tracking media +c + call gphysi +c +c make MARS invisible +c + call gsatt('MARS','SEEN',0) +c + call ggclos ! end of the geometry, must be called +c +c print some information about volumes and tracking media on screen +c +c call gprint('VOLU',0) +c call gptmed(0) +c call gpmate(17) +c call gpmate(18) +c +c----------------------------------------------------------------------------- +c +c from GEXAM.FOR example : +c +c the card PLOT is used in the following way: +c no PLOT card --> LPLOT(1)=0 --> no graphic output +c PLOT 1 --> LPLOT(1)=1 --> plotting onto workstation +c PLOT 2 --> LPLOT(1)=2 --> plottinig onto workstation +c and PostScript file + if (lplot(1) .gt. 0) then +c + if ( lview(1) .eq. 1 ) then + gdtheta = lview(2) + gdphi = lview(3) + gdpsi = lview(4) + gdx0 = lview(5) + gdy0 = lview(6) + gdsx = lview(7)/10. + gdsy = lview(8)/10. + else + gdtheta = 60. + gdphi = 40. + gdpsi = 0. + gdx0 = 10. + gdy0 = 10. + gdsx = 0.2 + gdsy = 0.2 + endif +c +c + call iginit(0) ! init HIGZ +c +c ask for worktsation type +c this defines the geometry of the X11 window, according to +c HIGZ_WINDOWS.DAT +c + write(6,'(''0Enter the workstation type : '')') + read(5,*) iwktyp + call hplint(iwktyp) + call igset('2BUF',11.) +c + if (lplot(1) .gt. 1) then +c +c open graphics output file +c + psfile = OUT_DIR//'geant_lemsr_'//runno//'.ps' + open(unit=10, file=psfile, status='NEW', + 1 form='FORMATTED', + 2 access='SEQUENTIAL', err=1001) + + write(6,*) '%GEANT-LEMSR-I- succesfully opened '// + 1 'metafile ' + write(6,*) + write(6,*) psfile +c +c open and activate workstation for metafile output +c workstation ID = 5 +c goes to file LUN = 10 +c type = -4112 = PostScript, Landscape, A4 (--> HIGZ manual) +c + call iopwk(5, 10, -4112) + call iacwk(5) +c +c find out and select the appropriate transformation and viewport +c + call igqwk(5,'MXDS',r) + xn = min( 1., r(1)/r(2) ) + yn = min( 1., r(2)/r(1) ) + call isvp( 5, 0., xn, 0., yn) + scal = 20. / min( xn, yn ) + call iswn( 5, 0., scal, 0., scal) + call iselnt(5) + endif ! if LPLOT(1) .GT. 1 +c +c initialize the GEANT's drawing package before the first plot +c + call gdinit +c +c draw the chamber +c + call gdraw('MARS',gdtheta,gdphi,gdpsi,gdx0,gdy0,gdsx,gdsy) + call gdaxis(-10.,15.,-20.,8.) + call gdscal(gdx0, gdy0) + endif ! IF LPLOT(1) .GT. 0 +c +c------------------------------------------------------------------------------ +c +c open NTUPLE file and book NTUPLE +c + ntfile = OUT_DIR//'geant_lemsr_'//runno//'.nt' +c + iquest(10) = 65000 ! to allocate enough space for NTUPLE + call hropen(80,'GEANT',ntfile,'NQ',4096,istat) + if (istat.eq.0) then + write(6,*) '%GEANT-LEMSR-I- succesfully opened NTUPLE file ' + write(6,*) ' ' + write(6,*) ntfile + else + write(6,'('' %GEANT-LEMSR-E- error opening NTUPLE file '',/,a)') + 1 ntfile + call exit + endif +c +c book CWN +c + call hbnt(111,'GEANT_LEMSR', 'D') ! 'D' disk resident NTUPLE +c +c now define the variable blocks for CWN +c + if ( l_init_par .or. l_allv ) then + call hbname(111,'init_par', nt_x0, + 1 'x0:r*4, y0:r*4, z0:r*4,'// + 2 'px0:r*4, py0:r*4, pz0:r*4, p0:r*4') + endif +c + if ( l_code .or. l_allv ) then + call hbname(111,'partcode', nt_ipart, 'Ipart:u*4:9') + endif + if ( l_endp .or. l_allv ) then + call hbname(111,'endpos' , nt_xe, + 1 'xe:r*4, ye:r*4, ze:r*4, epos:r*4') + endif + if ( l_time .or. l_allv ) then + call hbname(111,'times' , nt_tmcp, 'tmcp:r*4, tsci:r*4') + endif + if ( l_scip .or. l_allv ) then + call hbname(111,'sci_pos' , nt_xsci, + 1 'xsci:r*4, ysci:r*4, zsci:r*4') + endif + if ( l_scie .or. l_allv ) then + call hbname(111,'sci_elos', nt_desci, + 1 'desci:r*4, descipos:r*4, descigam:r*4, desciele:r*4') + endif +c + if ( l_scop .or. l_allv ) then + call hbname(111,'sco_pos' , nt_xsco, + 1 'xsco:r*4, ysco:r*4, zsco:r*4') + endif + if ( l_scoe .or. l_allv ) then + call hbname(111,'sco_elos', nt_desco, + 1 'desco:r*4, descopos:r*4, descogam:r*4, descoele:r*4') + endif +c + if ( l_flag .or. l_allv ) then + call hbname(111,'flags' , nt_volno, 'volno:u*4:8,'// + 1 'partcode:u*4:7') + endif + if ( l_nhit .or. l_allv ) then + call hbname(111,'nhits' , nt_mcpnhits, + 1 'mcpnhits:u*4:8, scinhits:u*4:8, sconhits:u*4:8') + endif + if ( l_posm .or. l_allv ) then + call hbname(111,'pos_mom' , nt_pxpos, + 1 'pxpos:r*4, pypos:r*4, pzpos:r*4') + endif + if ( l_mcpe .or. l_allv ) then + call hbname(111,'de_mcp' , nt_demcp, 'demcp:r*4') + endif +c + if ( l_ltrb .or. l_allv ) then + call hbname(111,'de_ltrb' , nt_descil, 'de_left_i:r*4,'// + 1 'de_top_i:r*4, de_rite_i:r*4, de_bot_i:r*4,'// + 2 'de_left_o:r*4,'// + 3 'de_top_o:r*4, de_rite_o:r*4, de_bot_o:r*4') + endif +c + if ( lsets(7) .eq. 1 ) then ! read initial data from file +c + write(6,*) ' ' + write(6,*) ' lsets(7) = 1 ---> I will read data from file ', + 1 INFILE + open(70, file=INFILE, form='UNFORMATTED', status='OLD', + 1 err=2000) + write(6,*) ' ' + write(6,*) '%GEANT_LEMSR-I- Successfully opened file ', INFILE + read(70) nevent + write(6,*) ' ' + write(6,*) ' Number of events in file = ', nevent + write(6,*) ' ' + endif +c +c----------------------------------------------------------------- c +c +c E N D OF I N I T I A L I S A T I O N +c c +c----------------------------------------------------------------- +c +c now do the tracking +c + call grun +c +c update all workstations +c + if ( lplot(1) .gt. 0 ) call iuwk(0,1) +c +c accept some input to have time to look at the picture +c + if ( lplot(1) .gt. 0 ) then + write(6,*) ' ' + write(6,*) ' type anything to continue ' + read(5,*) + endif + if ( lplot(1) .gt. 1 ) then + call iclwk(5) ! advance paper + close(10) + endif +c +c end GEANT and HPLOT +c + call glast +c + if ( lplot(1) .gt. 0 ) call igend +c + if ( lsets(7) .eq. 1) then + close(70) + write(6,*) ' ' + write(6,*) '%GEANT-LEMSR-I- closed input file ', INFILE + write(6,*) ' ' + endif +c +c write data and close NTUPLE file +c + call hrout(111,i,' ') + call hrend('GEANT') + close(80) + write(6,*) '%GEANT-LEMSR-I- closed NTUPLE file ', ntfile + write(6,*) '%GEANT-LEMSR-I- results of NEUMANN rejection for '// + 1 'MICHEL spectrum : ' + write(6,'(t10,'' rejected trials : '',i10)') rejcou + write(6,'(t10,'' succes. trials : '',i10)') cou + write(6,*) ' ' + write(6,'(t5,'' RNDM generator seed at the end : '')') + write(6,'(t10,'' ix1 = '',i15)') ix1 + write(6,*) ' ' + call exit +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c error messages +c +1000 continue + write(6,'(''0%GEANT-LEMSR-E- error reading GDRAW parameter...'')') + call exit +c +1001 continue + write(6,'(''0%GEANT-LEMSR-E- error opening '',a)') psfile + call exit +c +1002 continue + write(6,'(''0%GEANT-LEMSR-E- error write RunNo to string...'')') + write(6,*) ' may be, you entered a run number with more than ' + write(6,*) ' 4 digits...' + call exit +c +1003 continue + write(6,'(''0%GEANT-LEMSR-E- error opening input file '',a)') + 1 PARFILE + call exit +c +2000 continue + write(6,*) '%GEANT-LEMSR-E- error opening input file ', INFILE + call exit + end diff --git a/geant3/src/lemsr/get_direction.F b/geant3/src/lemsr/get_direction.F new file mode 100644 index 0000000..2e39ebd --- /dev/null +++ b/geant3/src/lemsr/get_direction.F @@ -0,0 +1,77 @@ +c +c---------------------------------------------------------------------------- +c + subroutine get_direction(costheta, sintheta, cosphi, sinphi, + 1 pos_tot) +c +c +c T. Prokscha, 15-Sep-2000 PSI +c +c Unix Version +c +c------------------------------------------------------------------------------ +c +c generates the initial direction of the decay-e+ +c + implicit none +c +#define CERNLIB_TYPE +#include "geant321/gctrak.inc" +#include "common.inc" +c + real*4 pos_tot ! total e+ energy including mass + real*4 p ! polarization == bfield(4) + real*4 costheta,phi ! decay direction in rotated system + real*4 sintheta + real*4 cosphi, sinphi + real*4 eps ! pos_tot/MAX_TOTAL_ENERGY + real*4 Dw ! asymmetry factor + real*4 random(3) + real*4 x,y,z,r +c +c----------------------------------------------------------------------------- +c + x = vect(1) + y = vect(2) + z = vect(3) + r = sqrt(x**2 + y**2) + if (1.4 .le. z .and. z .le. 1.41 .and. r .le. 3.5) then + p = bfield(5) / 100. + else + p = bfield(4) / 100. + endif +c +#if defined (OS_UNIX) + call ranlux(random, 3) ! random generator from Mathlib +#else + random(1) = ran(ix1) + random(2) = ran(ix1) + random(3) = ran(ix1) +#endif + if (p.eq.0.0) then +c + costheta = 1. - 2.*random(1) +c + else +c +c see TP Logbook 7, p. 167 for getting the equations +c + eps = pos_tot / MAX_TOTAL_ENERGY + Dw = (2. * eps - 1.) / (3. - 2. * eps) * p + + costheta = 1/Dw*(-1. + + 1 sqrt( 1. - 2*Dw*( 2.*random(2) - 1. ) + Dw**2)) + endif +c + sintheta = sqrt(1. - costheta**2) +c + phi=random(3) * TWO_PI +c + cosphi = cos(phi) + sinphi = sin(phi) +c +c now we have the decay direction in the rest system of the muon +c return. +c + return + end diff --git a/geant3/src/lemsr/gudcay.F b/geant3/src/lemsr/gudcay.F new file mode 100644 index 0000000..174c43c --- /dev/null +++ b/geant3/src/lemsr/gudcay.F @@ -0,0 +1,110 @@ +c +#include "geant321/pilot.h" +c--------------------------------------------------------------------------- +c + subroutine gudcay +c +c in case of muon started with partcode = 500, throw michel spectrum +c and anisotropic angular distribution +c +c T. Prokscha, 15-Sep-2000 PSI +c +c Unix version +c +c TP, 05-Apr-2001 PSI Unix, NT, Linux +c +c--------------------------------------------------------------------------- +c + implicit none +c +#define CERNLIB_TYPE +#include "geant321/gctrak.inc" +#include "geant321/gckine.inc" +#include "geant321/gcking.inc" +c +#include "common.inc" +c +c Michel spectrum +c + real*4 pos_tot ! total positron energy (including rest mass) + real*4 p0 ! total positron momentum (MeV) + real*4 px0, py0, pz0 ! components of p0 +c +c decay direction in muon rest system (Spin parallel to x-axis !!!) +c + real*4 costheta, sintheta, cosphi, sinphi +c +c 4-momentum of muon particle (if decay in flight); if decay in flight +c do a Lorentz-transformation to get momenta in real space +c + real*4 decay_rest(4), beta(4), decay_lab(4), gamma, beta_tot + real*4 decay_rot(3) ! copy of momentum comp. of decay_lab(4) +c +c---------------------------------------------------------------------------- +c +c get the Michel spectrum +c + call michel(pos_tot, p0) +c + call get_direction(costheta, sintheta, cosphi, sinphi, pos_tot) +c +c get the single momentum components in the laboratory system +c the muon spin is directed along the x-axis which points to the +c LEFT detector.) +c + px0 = p0 * costheta + py0 = p0 * sintheta * cosphi + pz0 = p0 * sintheta * sinphi +c + decay_rest(1) = px0/1000. ! GEANT uses GeV/c instead of MeV/c + decay_rest(2) = py0/1000. + decay_rest(3) = pz0/1000. + decay_rest(4) = pos_tot/1000. +c +c for the Lorentz transformation, maybe not important, because the +c mu+ are slow +c + gamma = getot / MU_MASS + beta_tot = vect(7) / getot ! this is p/E, vect(7) = p + beta(1) = -vect(4) * beta_tot ! vect(4) = px/p + beta(2) = -vect(5) * beta_tot ! vect(5) = py/p + beta(3) = -vect(6) * beta_tot ! vect(6) = pz/p +c +c negative sign because we need the velocity of +c the lab-system with respect to the rest-system +c + beta(4) = gamma +c + call gloren(beta, decay_rest, decay_lab) +c +c if there is a B-field present we have to rotate the momentum vector +c at the moment only Bz implemented; exit program if Bx or By are +c given. +c +c + if ( bfield(1) .ne. 0. .or. bfield(2) .ne. 0. or. bfield(3) .ne. + 1 0.) then + decay_rot(1) = decay_lab(1) + decay_rot(2) = decay_lab(2) + decay_rot(3) = decay_lab(3) +c + call rotate(decay_rot, tofg) +c + decay_lab(1) = decay_rot(1) + decay_lab(2) = decay_rot(2) + decay_lab(3) = decay_rot(3) + endif +c + ngkine = 1 + tofd(1) = 0. + gkin(1,1) = decay_lab(1) + gkin(2,1) = decay_lab(2) + gkin(3,1) = decay_lab(3) + gkin(4,1) = decay_lab(4) + gkin(5,1) = 2 + gpos(1,1) = vect(1) + gpos(2,1) = vect(2) + gpos(3,1) = vect(3) +c + return + end diff --git a/geant3/src/lemsr/gukine.F b/geant3/src/lemsr/gukine.F new file mode 100644 index 0000000..839d93e --- /dev/null +++ b/geant3/src/lemsr/gukine.F @@ -0,0 +1,229 @@ +* +#include "geant321/pilot.h" +* +c------------------------------------------------------------------------------ +c + subroutine gukine +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version +c +c TP, 05-Apr-2001 PSI: Unix, NT, Linux +c +c------------------------------------------------------------------------------ +c + implicit none +c +#define CERNLIB_TYPE +#include "geant321/gckine.inc" +#include "geant321/gclist.inc" +#include "geant321/gctrak.inc" +c +#include "cwn.inc" +#include "common.inc" +c +c----------------------------------------------------------------------------- +c +c local variables +c + real*4 decay_rot(3) ! rotation of muon spin, +c ! rotated e+-momenta + real*4 x0, y0, z0, px0, py0, pz0, p0, pos_tot + real*4 radius, phi0, cosphi, sinphi, costheta, sintheta + real*4 theta2d, theta_beam, sintheta_beam, costheta_beam + real*4 proj, phi_beam, sinphi_beam, cosphi_beam, sigma_beam + integer*4 ntbeam,nttarg,nubuf,i + parameter (nubuf = 1) + real*4 ubuf(nubuf) + real*4 random_1(3), random_2(7) + real*4 sigma_p +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + ipart = lsets(1) +c +c----------------------------------------------------------------------------- +c +c throw decay point at MCP, isotropic distribution and MICHEL-spectrum +c + if ( lsets(2) .eq. 0 .and. lsets(7) .ne. 1) then +c +#if defined (OS_UNIX) + call ranlux(random_1,3) +#else + do i = 1, 3 + random_1(i) = ran(ix1) + enddo +#endif +c +c------------------ +c +c the decay point at MCP or sample +c +c z0 = +1.6 ! distance of MCP center to center of MARS +c z0 = +1.48002 ! the Hamamatsu MCP's end is at z=1.48cm +c +c set MCP2 to 1.52cm; doing this we have the same initial z-position +c for the e+ with MCP mounting and sample mounting: +c +c Cu ring of sample (z/2 = 0.2cm) at z = 1.6 cm +c MCP Hamamatsu Z-stack (z/2 = 0.12cm) at z = 1.52 cm +c + z0 = +1.40002 ! the Hamamatsu, Galileo MCP's end and sample end is at z=1.40cm now +c ! z=1.40002 corresponds to an implatation depth of the +c ! mu+ of 200 nm. +c + radius = R_MCP * sqrt(random_1(1)) ! uniform distribution on MCP +c radius = 4.7 + phi0 = TWO_PI * random_1(2) + x0 = radius * cos(phi0) + y0 = radius * sin(phi0) +c +c--------------------------- +c +c the decay time +c + tofg = -TAU_MUON*alog(1.-random_1(3)) +c +c the michel spectrum +c + call michel(pos_tot, p0) +c +c anisotropic angular distribution, P=100% +c + call get_direction(costheta, sintheta, cosphi, sinphi, + 1 pos_tot) +c +c get the single momentum components in the laboratory system +c (the muon spin is directed along the x-axis which points to the +c LEFT detector.) +c + if ( bfield(1) .ne. 0. or. bfield(2) .ne. 0. or. bfield(3) + 1 .ne. 0.) then +c +c rotate when there is a B-field present; the subroutine +c ROTATE(...) will return without any action if bfield(4)=pol.=0. +c + decay_rot(1) = p0 * costheta + decay_rot(2) = p0 * sintheta * cosphi + decay_rot(3) = p0 * sintheta * sinphi +c + call rotate(decay_rot, tofg) +c + px0 = decay_rot(1) + py0 = decay_rot(2) + pz0 = decay_rot(3) +c + else +c + px0 = p0 * costheta + py0 = p0 * sintheta * cosphi + pz0 = p0 * sintheta * sinphi +c + endif +c +c------------------------- +c +c read input data from file +c + else if ( lsets(7) .eq. 1 ) then +c + read(70) x0, y0, z0, px0, py0, pz0 + z0 = z0-0.3 +c +c------------------------- +c +c for beam particles hitting MCP with momentum defined by +c lsets(3), lsets(4), lsets(5) +c + else +#if defined (OS_UNIX) + call ranlux(random_2,7) +#else + do i = 1, 5 + random_2(i) = ran(ix1) + enddo +#endif +c + +c z0 = -74. ! start in vacuum upstream from MCP + z0 = float(lsets(2)) ! start in vacuum at z=z0; + radius = R_MCP * sqrt(random_2(1)) ! uniform distribution on MCP + phi0 = TWO_PI * random_2(2) + x0 = radius * cos(phi0) + y0 = radius * sin(phi0)+0.25 +c +c if sets(8) > 0 throw gaussian momentum distribution +c + p0 = pkine(1) + if ( lsets(8) .gt. 0 ) then + sigma_p = p0 * (lsets(8)/100. + lsets(9)/1000.) + theta2d = sigma_p * + 1 sqrt( -2. * alog(1.-random_2(6))) + proj = sin( TWO_PI * random_2(7)) + p0 = p0 + theta2d*proj + endif +c +c throw THETA gaussian distributed to get some divergence of +c the beam +c + sigma_beam = float(lsets(6)) ! in degree +c + if ( sigma_beam .eq. 0. ) then + px0 = 0. + py0 = 0. + pz0 = p0 + else + sigma_beam = sigma_beam / 57.29577951 ! in rad + theta2d = sigma_beam * + 1 sqrt( -2. * alog(1.-random_2(3))) + proj = sin( TWO_PI * random_2(4)) + theta_beam = proj * theta2d + sintheta_beam= sin(theta_beam) + costheta_beam= cos(theta_beam) +c + phi_beam = TWO_PI * random_2(5) + sinphi_beam = sin(phi_beam) + cosphi_beam = cos(phi_beam) +c + px0 = p0 * sintheta_beam * cosphi_beam + py0 = p0 * sintheta_beam * sinphi_beam + pz0 = p0 * costheta_beam +c + endif +c + endif +c +c----------------------------------------------------------------------------- +c +c type*,' ivert = ',ivert +c + vert(1) = x0 + vert(2) = y0 + vert(3) = z0 + pvert(1) = px0/1000. ! GEANT uses GeV/c + pvert(2) = py0/1000. + pvert(3) = pz0/1000. + nt_x0 = x0 + nt_y0 = y0 + nt_z0 = z0 + nt_px0 = px0 + nt_py0 = py0 + nt_pz0 = pz0 + nt_p0 = p0 + nt_ipart = ipart +c +c write(6,*) x0,y0,z0,px0,py0,pz0 +c +c store "vertex" (this means in this case the initial position and +c momentum of the paricle) data to the GEANT data structures +c JVERTX and JKINE. It has to be done ! +c + call gsvert(vert,ntbeam,nttarg,ubuf,nubuf,ivert) + call gskine(pvert,ipart,ivert,ubuf,nubuf,itra) +c call gpvert(ivert) +c call gpkine(itra) +c + return + end diff --git a/geant3/src/lemsr/guout.F b/geant3/src/lemsr/guout.F new file mode 100644 index 0000000..2639c29 --- /dev/null +++ b/geant3/src/lemsr/guout.F @@ -0,0 +1,439 @@ +#include "geant321/pilot.h" +c +c----------------------------------------------------------------------------- +c + subroutine guout +c +c this routine writes the data at the end of one event to NTUPLE. +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix Version, part of the lemsr_geant program +c +c TP, 05-Apr-2001, PSI Unix, NT, Linux version +c +c----------------------------------------------------------------------------- +c + implicit none +c +#define CERNLIB_TYPE +#include "geant321/gcsets.inc" +#include "geant321/gcflag.inc" +#include "geant321/gclist.inc" +c +c ntuple variables and other common variablea +c +#include "cwn.inc" +#include "common.inc" +c +c------------------------------------------------------------------------- +c + integer*4 i, j, NVDIM, NHMAX, NHDIM +c + parameter (NVDIM = 1) + parameter (NHMAX = 100) + parameter (NHDIM = 7) + integer*4 numvs(NVDIM) + integer*4 itra(NHMAX) + integer*4 nbv(NVDIM,NHMAX) + real*4 hits(NHDIM,NHMAX) + integer*4 nhits, maxind +c + real*4 detot, descipos, descigam, desciele ! energy loss total, +c ! and in Sc. due +c ! to e+, gammas, e- +c + real*4 descil, descit, descir, descib ! energy loss in one detec. +c + integer*4 bit + integer*4 volno ! volume number; this is to decide which +c ! scintillator was hit: +c ! volno = 1 (LEFT) ----> write 1 to NTuple variable +c ! = 2 (TOP) ----> write 2 +c ! = 3 (RIGHT) ----> write 4 +c ! = 4 (BOTTOM)----> write 8 +c ---> sum variable entry if more than one detector +c was hit. +c ! for outer scintillators, start LEFT ---> 16 +c and so on. +c + integer*4 partcode ! store in NTuple which particles had contributed +c ! partcode = 1 (Gamma) ---> write 1 to NT variable +c ! partcode = 2 (e+) ---> write 2 to NT +c ! partcode = 3 (e-) ---> write 4 to NT +c ---> sum NT variable values if two or three +c different particles were generated. +c for outer scintillators start with +c GAMMA ---> 8 +c + real*4 ubuf(1) + integer*4 nubuf /1/ + integer*4 nvert, part + real*4 vert(3), pvert(3) ! e+ vertex values +c +c----------------------------------------------------------------------------- +c + if ( lsets(1) .eq. 5 .or. lsets(1) .eq. 500) then +c ! if initial particle was a muon, +c ! then, look for e+ momentum which +c ! is created as 2. vertex +c ! "2" is therefore the 1. argument +c ! of gfkine +c + call gfkine(2, vert, pvert, part, nvert, ubuf, nubuf) +c + nt_pxpos = 1000. * pvert(1) ! MeV/c + nt_pypos = 1000. * pvert(2) ! MeV/c + nt_pzpos = 1000. * pvert(3) ! MeV/c +c + endif +c +c now look for detector hits to get information on energy loss, +c particle and which detector has a hit, and fill NTUPLE variuables +c +c------------------------------ +c +c the MCP2 detector +c + do i = 1, NHDIM ! reset first HITS(,) array + do j = 1, NHMAX + hits(i, j) = 0. + enddo + enddo + detot = 0. +c + call gfhits('MDET','DMCP',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv, + 1 hits, nhits) +c + if ( nhits .gt. 0 ) then ! ~MCP2 Nhits>0 + if ( nhits .le. NHMAX ) then + maxind = nhits + else + maxind = NHMAX + endif + do j = 1, maxind +c + detot = detot + hits(6,j) ! sum the energy loss per step +c +c write(6,'(1x,7E10.3)') hits(1,j),hits(2,j),hits(3,j),hits(4,j), +c + hits(5,j),hits(6,j),hits(7,j) + enddo +c +c now fill MCP relevant data to NTUPLE variables +c + nt_tmcp = int(hits(4,1)*1.e9) ! decay time at MCP in ns + nt_demcp = 1000. * detot ! total energy loss in MCP in MeV +c + endif ! ~ endif MCP2 hits >0 +c + nt_mcpnhits = nhits +c +c------------------------------------- +c +c the inner scintillators +c + do i = 1, NHDIM ! reset first HITS(,) array + do j = 1, NHMAX + hits(i, j) = 0. + enddo + enddo +c + nhits = 0 + detot = 0. + descipos = 0. + descigam = 0. + desciele = 0. + volno = 0 + partcode = 0 + descil = 0. + descit = 0. + descir = 0. + descib = 0. +c + call gfhits('ISSC','SCIS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits, + 1 nhits) +c +c write(6,*) 'Number of Hits in inner Sc.: ',nhits +c + if ( nhits .gt. 0 ) then ! ~ hit in inner Sc. +c + if ( nhits .le. NHMAX ) then + maxind = nhits + else + maxind = NHMAX + endif + do j = 1, maxind + detot = detot + hits(6,j) ! sum the energy loss per step + if ( hits(7,j) .eq. 1. ) then ! particle is gamma + descigam = descigam + hits(6,j) + endif + if ( hits(7,j) .eq. 2. ) then ! particle is positron + descipos = descipos + hits(6,j) + endif + if ( hits(7,j) .eq. 3. ) then ! particle is electron + desciele = desciele + hits(6,j) + endif +c +c now look for volume +c + bit = nbv(1,j) - 1 + volno = ibset(volno, bit) +c +c sum up the energy deposited in one detector +c + if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector + descil = descil + hits(6,j) + endif + if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector + descit = descit + hits(6,j) + endif + if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector + descir = descir + hits(6,j) + endif + if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector + descib = descib + hits(6,j) + endif +c +c look for particles +c + bit = int(hits(7,j)) - 1 + partcode = ibset(partcode, bit) +c +c write(6,'(1x,7E10.3,2I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j), +c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno +c + enddo +c + nt_tsci = hits(4,1)*1.e9 ! decay time at Sci. in ns + nt_xsci = hits(1,1) + nt_ysci = hits(2,1) + nt_zsci = hits(3,1) + nt_desci = detot * 1000. + nt_descipos = descipos * 1000. + nt_descigam = descigam * 1000. + nt_desciele = desciele * 1000. + nt_descil = descil * 1000. + nt_descit = descit * 1000. + nt_descir = descir * 1000. + nt_descib = descib * 1000. + nt_volno = volno + nt_partcode = partcode +c + endif ! ~ hit in inner Sc. +c + nt_scinhits = nhits +c +c end of inner sc. +c +c------------------------------- +c +c the outer scintillators +c + do i = 1, NHDIM ! reset first HITS(,) array + do j = 1, NHMAX + hits(i, j) = 0. + enddo + enddo +c + nhits = 0 + detot = 0. + descipos = 0. + descigam = 0. + desciele = 0. + descil = 0. + descit = 0. + descir = 0. + descib = 0. +c +c + call gfhits('OSSC','SCOS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits, + 1 nhits) +c +c write(6,*) 'Number of Hits in outer Sc.: ',nhits +c + if ( nhits .gt. 0 ) then ! ~ hit in outer Sc. +c + if ( nhits .le. NHMAX ) then + maxind = nhits + else + maxind = NHMAX + endif + do j = 1, maxind + detot = detot + hits(6,j) ! sum the energy loss per step + if ( hits(7,j) .eq. 1. ) then ! particle is gamma + descigam = descigam + hits(6,j) + endif + if ( hits(7,j) .eq. 2. ) then ! particle is positron + descipos = descipos + hits(6,j) + endif + if ( hits(7,j) .eq. 3. ) then ! particle is electron + desciele = desciele + hits(6,j) + endif +c +c now look for volume +c + bit = nbv(1,j) - 1 + 4 + volno = ibset(volno, bit) +c +c sum up the energy deposited in one detector +c + if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector + descil = descil + hits(6,j) + endif + if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector + descit = descit + hits(6,j) + endif + if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector + descir = descir + hits(6,j) + endif + if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector + descib = descib + hits(6,j) + endif +c +c look for particles +c + bit = int(hits(7,j)) - 1 + 4 + partcode = ibset(partcode, bit) +c +c write(6,'(1x,7E10.3,3I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j), +c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno,idtype + enddo +c +c +c nt_tsco = hits(4,1)*1.e9 ! decay time at Sci. in ns + nt_xsco = hits(1,1) + nt_ysco = hits(2,1) + nt_zsco = hits(3,1) + nt_desco = detot * 1000. + nt_descopos = descipos * 1000. + nt_descogam = descigam * 1000. + nt_descoele = desciele * 1000. + nt_descol = descil * 1000. + nt_descot = descit * 1000. + nt_descor = descir * 1000. + nt_descob = descib * 1000. + nt_volno = volno + nt_partcode = partcode +c + endif ! ~ hit in inner Sc. +c + nt_sconhits = nhits +c +c end out. scintillators +c +c--------------------------------- +c +c fill NTUPLE +c + call hfnt(111) +c +c reset NTUPLE variables +c + nt_x0 = 0. + nt_y0 = 0. + nt_z0 = 0. + nt_px0 = 0. + nt_py0 = 0. + nt_pz0 = 0. + nt_p0 = 0. + nt_ipart = 0 + nt_xe = 0. + nt_ye = 0. + nt_ze = 0. + nt_eposend = 0. + nt_tmcp = 0 + nt_tsci = 0 + nt_xsci = 0. + nt_ysci = 0. + nt_zsci = 0. + nt_desci = 0. + nt_descipos = 0. + nt_descigam = 0. + nt_desciele = 0. +c + nt_xsco = 0. + nt_ysco = 0. + nt_zsco = 0. + nt_desco = 0. + nt_descopos = 0. + nt_descogam = 0. + nt_descoele = 0. +c + nt_volno = 0 + nt_partcode = 0 + nt_mcpnhits = 0 + nt_scinhits = 0 + nt_sconhits = 0 +c + nt_pxpos = 0. + nt_pypos = 0. + nt_pzpos = 0. + nt_demcp = 0. +c + nt_descil = 0. + nt_descit = 0. + nt_descir = 0. + nt_descib = 0. +c + nt_descol = 0. + nt_descot = 0. + nt_descor = 0. + nt_descob = 0. +c +c----------------------------------- +c +c some information on progress of simulation +c + if ( float(idevt)/float(nevent) .ge. 0.1 .and. percent(1) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 10%' + percent(1) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.2 .and. percent(2) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 20%' + percent(2) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.3 .and. percent(3) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 30%' + percent(3) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.4 .and. percent(4) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 40%' + percent(4) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.5 .and. percent(5) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 50%' + percent(5) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.6 .and. percent(6) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 60%' + percent(6) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.7 .and. percent(7) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 70%' + percent(7) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.8 .and. percent(8) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 80%' + percent(8) = 1 + endif + if ( float(idevt)/float(nevent) .ge. 0.9 .and. percent(9) .eq. 0) + 1 then + write(6,*) ' fraction of processed events: 90%' + percent(9) = 1 + endif + if ( idevt .eq. nevent) then + write(6,*) ' fraction of processed events: 100%' + endif +c +c----------------------------- +c + return + end diff --git a/geant3/src/lemsr/gustep.F b/geant3/src/lemsr/gustep.F new file mode 100644 index 0000000..fbc5123 --- /dev/null +++ b/geant3/src/lemsr/gustep.F @@ -0,0 +1,161 @@ +#include "geant321/pilot.h" +c +c---------------------------------------------------------------------- +c + subroutine gustep +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix Version, part of the geant_lemsr program +c +c TP, 05-Apr-2001, PSI Unix, NT, Linux +c +c---------------------------------------------------------------------- +c + implicit none +c +c include again some GEANT variables +c +#define CERNLIB_TYPE +#include "geant321/gckine.inc" +#include "geant321/gctrak.inc" +#include "geant321/gctmed.inc" +#include "geant321/gcking.inc" +#include "geant321/gcsets.inc" +#include "geant321/gcvolu.inc" +#include "geant321/gclist.inc" +#include "geant321/gcflag.inc" +c +c NTuple variable and other common variables +c +#include "cwn.inc" +#include "common.inc" +c + integer*4 j, part + integer*4 ihit +c + real*4 hits(7) +c +c------------------------------------------------------------------ +c +c store secondary particle data to GEANT structure JSTAK +c + do j = 1,ngkine + part = gkin(5,j) + if ( part .eq. 4 ) then ! neutrinos, discard them + iflgk(j) = -1 + elseif ( part .eq. 2) then ! save e+ in JKINE and transport + if ( l_pos) then + iflgk(j) = 1 + else + iflgk(j) = -1 + endif + elseif ( part .eq. 1) then + if ( l_gam ) then + iflgk(j) = 0 ! simply transport the rest + else + iflgk(j) = -1 + endif + elseif ( part .eq. 3) then + if ( l_ele ) then + iflgk(j) = 0 ! simply transport the rest + else + iflgk(j) = -1 + endif + else ! discard other part. in any + iflgk(j) = -1 ! case + endif + enddo + call gsking(0) ! store secondary particle data to JSTAK +c +c +c print*,' particle type = ',ipart +c print*,' tracking medium number = ',numed +c print*,' isvol = ',isvol +c print*,' current track position = ',vect(1),vect(2),vect(3) +c print*,' current momentum = ',vect(4),vect(5),vect(6) +c print*,' total momentum = ',vect(7) +c print*,' current kinetic energy = ',gekin +c print*,' number of steps for = ',nstep +c print*,' total energy lost = ',destep +c print*,' current TOF = ',tofg +c print*,' istop = ',istop +c print*,' igauto = ',idecad +c print*,' inwvol = ',inwvol +c print*,'---------------------' +c +c + if ( lplot(1) .gt. 0) call gdcxyz +c +c store spatial,momentum and energy loss data to GEANT data structure +c JXYZ +c +c call gsxyz +c call gpcxyz +c +c look where the e+ stops +c + if (ipart.eq.2.and.(inwvol.eq.3.or.istop.ge.1)) then +c + nt_xe = vect(1) + nt_ye = vect(2) + nt_ze = vect(3) + nt_eposend = 1000. * gekin +c + endif +c +c look now for detector hits +c + if ( idtype .eq. 1) then +c + hits(1) = vect(1) + hits(2) = vect(2) + hits(3) = vect(3) + hits(4) = tofg + hits(5) = gekin + hits(6) = destep + hits(7) = float(ipart) + call gsahit(iset,idet,itra,numbv,hits,ihit) ! MCP + if ( lplot(1) .gt. 0 ) then ! draw and plot hit + call gphits('MDET','DMCP') + call gdhits('MDET','DMCP',0,854,.3) + endif +c + endif +c + if ( idtype .eq. 2) then +c + hits(1) = vect(1) + hits(2) = vect(2) + hits(3) = vect(3) + hits(4) = tofg + hits(5) = gekin + hits(6) = destep + hits(7) = float(ipart) + call gsahit(iset,idet,itra,numbv,hits,ihit) ! inner Scint. + if ( lplot(1) .gt. 0 ) then ! plot hit + call gphits('ISSC','SCIS') + call gdhits('ISSC','SCIS',0,852,.3) + endif +c + endif +c + if ( idtype .eq. 3) then +c + hits(1) = vect(1) + hits(2) = vect(2) + hits(3) = vect(3) + hits(4) = tofg + hits(5) = gekin + hits(6) = destep + hits(7) = float(ipart) + call gsahit(iset,idet,itra,numbv,hits,ihit) ! outer Scint. +c + if ( lplot(1) .gt. 0 ) then ! plot hit + call gphits('OSSC','SCOS') + call gdhits('OSSC','SCOS',0,842,.3) + endif +c + endif + return + end diff --git a/geant3/src/lemsr/gutrack.F b/geant3/src/lemsr/gutrack.F new file mode 100644 index 0000000..4e9751f --- /dev/null +++ b/geant3/src/lemsr/gutrack.F @@ -0,0 +1,19 @@ +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + subroutine gutrack +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version, part of the geant_lemsr program +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + implicit none +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c type*,' routine GUTRACK' + call gtrack + return + end diff --git a/geant3/src/lemsr/gutreve.F b/geant3/src/lemsr/gutreve.F new file mode 100644 index 0000000..ac31409 --- /dev/null +++ b/geant3/src/lemsr/gutreve.F @@ -0,0 +1,29 @@ +c#include "geant321/pilot.h" +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + subroutine gutreve +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version, part of the geant_lemsr program +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + implicit none +c +c include again some GEANT variables +c +c INTEGER IKINE,ITRA,ISTAK,IVERT,IPART,ITRTYP,NAPART,IPAOLD +c REAL PKINE,AMASS,CHARGE,TLIFE,VERT,PVERT +c +c#include "geant321/gtkine.inc" +c#include "geant321/gckine.inc" +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c type*,' routine GUTREVE' + call gtreve +c + return + end diff --git a/geant3/src/lemsr/michel.F b/geant3/src/lemsr/michel.F new file mode 100644 index 0000000..de49295 --- /dev/null +++ b/geant3/src/lemsr/michel.F @@ -0,0 +1,56 @@ +c +c---------------------------------------------------------------------------- +c + subroutine michel(pos_tot, p0) +c +c routine to throw the e+ energy spectrum in muon decay. +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version +c +c TP, 05-Apr-2001, PSI Unix, NT, Linux +c +c---------------------------------------------------------------------------- +c + implicit none +c +#include "common.inc" +c + real*4 pos_tot, p0 + real*4 eta, temp, rho, ratio + real*4 random(2) +c +c--------------------------- +c +c MICHEL distributed e+ energy by using NEUMANN rejection +c "rejection function" is 2.25 * eta +c +10 continue +#if defined (OS_UNIX) + call ranlux(random, 2) +#else + random(1) = ran(ix1) + random(2) = ran(ix2) +#endif + eta = sqrt(random(1)) + temp = 2. * eta**2 * (3. - 2. * eta) + ratio = temp / (2.25 * eta) + rho = random(2) + if ( rho .gt. ratio ) then + rejcou = rejcou + 1 + goto 10 + else + cou = cou + 1 + endif +c + pos_tot = MAX_TOTAL_ENERGY * eta +c + if ( pos_tot .ge. POS_MASS ) then + p0 = sqrt(pos_tot**2 - SQR_POS_MASS) + else + p0 = 0. + endif +c + return + end diff --git a/geant3/src/lemsr/rotate.F b/geant3/src/lemsr/rotate.F new file mode 100644 index 0000000..941f24c --- /dev/null +++ b/geant3/src/lemsr/rotate.F @@ -0,0 +1,52 @@ +c +c---------------------------------------------------------------------------- +c + subroutine rotate(decay_rot, time) +c +c routine to rotate the e+-momentum vector in case of an applied +c bfield; BFIELD is a 4-component array with bfield(4) = polarization. +c If Polarization = 0 return. +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version, part of the geant_lemsr program +c +c----------------------------------------------------------------------------- +c + implicit none +c +#include "common.inc" +c + real*4 decay_rot(3), time + real*4 omega, prec_angle, costh, sinth, cosph, sinph +c +c------------------------------------------------ +c + if ( bfield(4) .eq. 0. ) return +c + if ( bfield(2) .ne. 0. ) then + omega = bfield(2) * TWO_PI * MU_GYRO ! cycle/nsec + prec_angle = omega * time + costh = cos(prec_angle) + sinth = -sin(prec_angle) ! clockwise rotation + cosph = 1. + sinph = 0. +c + call gdrot(decay_rot, costh, sinth, cosph, sinph) +c + endif +c + if ( bfield(3) .ne. 0. ) then + omega = bfield(3) * TWO_PI * MU_GYRO ! cycle/nsec + prec_angle = omega * time + cosph = cos(prec_angle) + sinph = -sin(prec_angle) ! clockwise rotation + costh = 1. + sinth = 0. +c + call gdrot(decay_rot, costh, sinth, cosph, sinph) +c + endif +c + return + end diff --git a/geant3/src/lemsr/ugeom.F b/geant3/src/lemsr/ugeom.F new file mode 100644 index 0000000..fee2751 --- /dev/null +++ b/geant3/src/lemsr/ugeom.F @@ -0,0 +1,903 @@ +#include "geant321/pilot.h" +c +c----------------------------------------------------------------------- +c + subroutine ugeom +c +c this subroutine defines the geometrical set-up. +c The z-axis is defines as the low-energy mu+ beam axis. +c This means x-axis --> horizontal, pointing to the left +c y-axis --> vertical, upwards. +c +c T. Prokscha, 15-Sep-2000, PSI +c +c Unix version, part of the geant_lemsr program +c +c TP, 05-Apr-2001, PSI Unix, NT, Linux +c +c------------------------------------------------------------------------ +c + implicit none +c +#include "common.inc" +c +#define CERNLIB_TYPE +#include "geant321/gctmed.inc" +#include "geant321/gcsets.inc" +#include "geant321/gclist.inc" +#include "geant321/gcflag.inc" +#include "geant321/gcphys.inc" +#include "geant321/gccuts.inc" +#include "geant321/gckine.inc" +c +c----------------------------------------------------------------------------- +c + real*4 par_box(3) + real*4 par_tube(3) + real*4 par_tubs(5) ! segment of a tube + real*4 par_cone(5) +c + real*4 theta1,theta2,theta3,phi1,phi2,phi3 !rotation angles + integer*4 ivol !for the side arm +c + character*4 namesv(4) !for detector descr. + character*4 namesh(7) ! to define hit + integer*4 nbitsh(7) ! paramters on MCP + real*4 orig(7),fact(7) +c + real*4 atomw(9),zatom(9),wmat(9),dens ! to define MCP material + +c + integer*4 mat_type ! material type for sample plate + real*4 off + real*4 he_shield_thickness + integer*4 Air, Vacuum, Al, Fe, Cu, Steel, Brass, Sapphire, + 1 Macor, MCP, NaI, Scint ! medium numbers +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + data namesh /'x ','y ','z ','t ','e ','elos', 'part'/ + data nbitsh /7*32/ + data orig /3*100.,0.,0.,0.,0./ + data fact /3*100.,1.e+9,1000000.,1000000.,1./ +c +c the values for x,y,z,t,e and eloss will be internally converted +c by GEANT to integers; +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c MCP material data are taken from J.L.WIZA,Microchannel Plate Detectors, +c Nucl.Instr.Meth. 162 (1979), 587-601 +c +c data dens /4./ ! density of MCP glass + data dens /2./ ! about half of the MCP is "empty" due to +c ! the MCP channels + data zatom /82.,8.,14.,19.,37.,56.,33.,55.,11./ + data wmat /.479,.258,.182,.042,.018,.013,.004,.002,.001/ + data atomw /207.19,15.999,28.086,39.10,85.47,137.33,74.92, + 1 132.91,22.99/ +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define materials +c + call gmate ! standard GEANT materials +c +c some additional materials +c +c define LEAD GLASS (MCP), NaI +c + call gsmixt(17,'MCPGLASS',atomw,zatom,dens,9,wmat) +c +c define NaI +c + zatom(1) = 11. + atomw(1) = 22.99 + zatom(2) = 53 + atomw(2) = 126.9 + wmat(1) = 1. + wmat(2) = 1. + dens = 3.67 + call gsmixt(18,'NaICALOR',atomw,zatom,dens,-2,wmat) +c +c define Scintillator NE102A +c + zatom(1) = 1. + atomw(1) = 1.0079 + zatom(2) = 6. + atomw(2) = 12.011 + wmat(1) = 1.104 + wmat(2) = 1. + dens = 1.032 + call gsmixt(19,'SCINT',atomw,zatom,dens,-2,wmat) +c +c stainless steel +c + dens = 7.930 + zatom(1) = 26. + zatom(2) = 24. + zatom(3) = 28. + atomw(1) = 55.85 + atomw(2) = 52.00 + atomw(3) = 58.70 + wmat(1) = 0.71 + wmat(2) = 0.18 + wmat(3) = 0.11 + call gsmixt(20,'STEEL',atomw,zatom,dens,3,wmat) +c +c brass +c + dens = 8.67 + zatom(1) = 29. + zatom(2) = 30. + atomw(1) = 63.55 + atomw(2) = 65.38 + wmat(1) = .80 + wmat(2) = .19 +c + call gsmixt(21,'BRASS',atomw,zatom,dens,2,wmat) +c +c define Sapphire == AL2O3 +c + zatom(1) = 13. + atomw(1) = 26.98 + zatom(2) = 8. + atomw(2) = 16.00 + wmat(1) = 2. + wmat(2) = 3. + dens = 3.985 +c + call gsmixt(22,'SAPPHIRE',atomw,zatom,dens,-2,wmat) +c +c define MACOR as SiO2 with higher density; the correct compound is +c +c (SiO2)46 (Al2O3)16 (MgO)17 (K2O)10 (B2O3)7 +c + zatom(1) = 14. + atomw(1) = 28.0855 + zatom(2) = 8. + atomw(2) = 16.00 + wmat(1) = 1. + wmat(2) = 2. + dens = 2.52 +c + call gsmixt(23,'MACOR',atomw,zatom,dens,-2,wmat) +c +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define tracking medium numbers corresponding to the material numbers +c stored in JMATE by calling GMATE +c +c 9 = Aluminum +c 10 = Iron +c 11 = Copper +c 15 = Air +c 16 = Vacuum +c 17 = MCP +c 18 = NaI <--- not used for present setup +c 19 = Scintillator +c 20 = Stainless steel +c 21 = brass <--- not used +c 22 = sapphire +c 23 = macor +c + isvol = 0 + ifield = 0 +c epsil = 0.0001 ! boundary crossing precision set to 1 micron + epsil = 0.001 ! boundary crossing precision set to 10 micron +c + call gstmed( 1,'Air $', 15, isvol, ifield, fieldm, + 1 tmaxfd, STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 2,'Vacuum $', 16, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 3,'Aluminum $', 9, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 4,'Iron $', 10, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 5,'Copper $', 11, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 6,'Steel $', 20, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 7,'Brass $', 21, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 8,'Sapphire $', 22, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed( 9,'Macor $', 23, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) +c + isvol = 1 ! sensitive volume; this means this volume is a detector + call gstmed(10,'MCPDetector$', 17, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed(11,'NaIDetector$', 18, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) + call gstmed(12,'Scintillato$', 19, isvol, ifield, fieldm, tmaxfd, + 1 STEMAX, deemax, epsil, stmin, 0, 0) +c + Air = 1 + Vacuum = 2 + Al = 3 + Fe = 4 + Cu = 5 + Steel = 6 + Brass = 7 + Sapphire = 8 + Macor = 9 + MCP = 10 + NaI = 11 + Scint = 12 +c +c redefine kinetic cut energies for muons in all media from its +c default value = 10 MeV to 10keV =0.00001 GeV +c +c ---> define cut energies via CUTS input cards +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define MOTHER VOLUME MARS +c + par_box(1) = 100. + par_box(2) = 100. + par_box(3) = 100. +c + call gsvolu('MARS','BOX ', Air, par_box, 3, ivol) + if (ivol.lt.0) goto 900 !this means an error during +c !call to GSVOLU +c------------------------------------ +c +c vacuum tube of sample or MCP2, data from Balzers offer +c + par_tube(1) = 0. + par_tube(2) = 7.65 + par_tube(3) = 16.2 +c + call gsvolu('MCPV','TUBE', Vacuum, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c now stainless steel part of tube +c + par_tube(1) = 7.65 + par_tube(2) = 7.95 + par_tube(3) = 16.2 +c + call gsvolu('MCPS','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c 160 CF flange upstream of MCP2 tube +c + par_tube(1) = 7.95 + par_tube(2) = 10.125 + par_tube(3) = 1.1 +c + call gsvolu('F160','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c 100 CF flange at the end +c + par_tube(1) = 0. + par_tube(2) = 7.7 + par_tube(3) = 1.0 +c + call gsvolu('F100','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c 200 CF flange upstream of MCP2 tube to connect to gate valve chamber +c + par_tube(1) = 7.65 + par_tube(2) = 10.325 + par_tube(3) = 1.2 +c + call gsvolu('F200','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c par_tube(1) = 0. +c par_tube(2) = 7.95 +c par_tube(3) = 1.2 +c +c call gsvolu('T200','TUBE', Vacuum, par_tube, 3, ivol) +c if (ivol.lt.0) goto 900 +c +c the gate valve chamber +c +c + par_tube(1) = 0. + par_tube(2) = 10.325 + par_tube(3) = 9.25 +c + call gsvolu('GATV','TUBE', Vacuum, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c for simplicity choose outer diameter of the Gate valve chamber equal +c to the outer diameter of 200 CF flange +c + par_tube(1) = 10.325 + par_tube(2) = 12.65 + par_tube(3) = 9.25 +c + call gsvolu('GATS','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c the L3 tube, vacuum +c + par_tube(1) = 0. + par_tube(2) = 10.0 + par_tube(3) = 23.2 +c + call gsvolu('L3VA','TUBE', Vacuum, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c the L3 tube, stainless steel +c + par_tube(1) = 10.0 + par_tube(2) = 10.3 + par_tube(3) = 23.2 +c + call gsvolu('L3ST','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c conical anode +c +c the downstream part of HV piece +c + par_cone(1) = 2.5 ! half length dz + par_cone(2) = 4.8 ! inside radius at -dz + par_cone(3) = 6.2 ! outside radius at -dz + par_cone(4) = 3.3 ! inside radius at dz + par_cone(5) = 3.6 ! outside radius at dz +c + call gsvolu('RA-E','CONE', Steel, par_cone, 5, ivol) + if (ivol.lt.0) goto 900 +c +c the middle part, still on HV +c + par_cone(1) = 2.0 ! half length dz + par_cone(2) = 5.8 ! inside radius at -dz + par_cone(3) = 6.2 ! outside radius at -dz + par_cone(4) = 4.8 ! inside radius at dz + par_cone(5) = 6.2 ! outside radius at dz +c + call gsvolu('RA-M','CONE', Steel, par_cone, 5, ivol) + if (ivol.lt.0) goto 900 +c +c ground cylinder, stainless steel +c + par_tube(1) = 5.8 + par_tube(2) = 6.2 + par_tube(3) = 5.8 +c + call gsvolu('RA-G','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c------------------------------------------------------ +c +c sample cryo things +c +czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz +c +c sample holder +c +c 1. Cu plate (sample holder) on Cold finger, 0.5cm +c + par_tube(1) = 0. + par_tube(2) = 3.5 + par_tube(3) = 0.25 +c + if ( samp_mat .eq. 'saal' ) then + mat_type = Al ! Al + else + mat_type = Cu ! CU + endif +c + call gsvolu('SAH1','TUBE',mat_type,par_tube,3,ivol) +c print*,'SAH1 ivol = ',ivol + if (ivol.lt.0) goto 900 +c +c 2. Cu plate (0.4cm), sample is mounted on this plate +c + par_tube(1) = 0. + par_tube(2) = 3.5 + par_tube(3) = 0.2 +c + call gsvolu('SAH2','TUBE',mat_type,par_tube,3,ivol) +c print*,'SAH2 ivol = ',ivol + if (ivol.lt.0) goto 900 +c +c now select new color if sample plates are Al +c + if ( mat_type .eq. Al ) then + call gsatt('SAH1', 'COLO', 4) + call gsatt('SAH2', 'COLO', 4) + endif +c +c Sample mounting ring (0.1cm), 4cm inner diameter +c + par_tube(1) = 2. + par_tube(2) = 3.5 + par_tube(3) = 0.05 +c + call gsvolu('SAH3','TUBE',mat_type,par_tube,3,ivol) +c print*,'SAH3 ivol = ',ivol + if ( mat_type .eq. Al ) then + call gsatt('SAH3', 'COLO', 4) ! set color for aluminum + endif + if (ivol.lt.0) goto 900 + if ( l_samp) then + write(6,*) ' ' + if ( mat_type .eq. Cu ) then + write(6,*) ' Sample holder ring: Cu.' + else if ( mat_type .eq. Al ) then + write(6,*) ' Sample holder ring: Al.' + else + write(6,*) ' Sample holder ring: ', mat_type + endif + write(6,*) ' ' + endif +c +c +c Sapphire mounted between 1. and 2. Cu plate, 0.6cm thick, 6cm diameter +c + par_tube(1) = 0. + par_tube(2) = 3.0 + par_tube(3) = 0.3 +c + call gsvolu('SAPH','TUBE', Sapphire, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 + call gsatt('SAPH', 'COLO', 3) ! green color for saphhire +c +czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz +c +c cryo +c +c 1 cm plate, "cold finger" +c + par_tube(1) = 0. + par_tube(2) = 2.75 + par_tube(3) = 0.5 +c + call gsvolu('COFI','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c end plate of cryo (3cm diameter, 0.7cm thick) +c + par_tube(1) = 0. + par_tube(2) = 1.5 + par_tube(3) = 0.35 +c + call gsvolu('CRY1','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c "Waermetauscher", assume that it has 1cm open cylinder ? +c 5 cm long +c + par_tube(1) = 0.5 + par_tube(2) = 1.5 + par_tube(3) = 2.5 +c + call gsvolu('CRY2','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c the mounting ring for He-shield +c + par_tube(1) = 3.8 + par_tube(2) = 4.7 + par_tube(3) = 0.55 +c + call gsvolu('CRY3','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c add 2mm thick plate for mounting ring +c this is just to close the downstream side, I don't know if +c dimensions and position are right. +c + par_tube(1) = 1.5 + par_tube(2) = 3.8 + par_tube(3) = 0.1 +c + call gsvolu('CRY4','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c call gsatt('CRY4', 'COLO', 3) +czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz +c +c He shield +c + par_tube(1) = 4.7 + par_tube(2) = 4.8 + par_tube(3) = 4.5 +c + he_shield_thickness = par_tube(2) - par_tube(1) + if ( l_crsh ) then + write(6,*) ' ' + write(6,*) 'Thickness of Cu-He-shield: ',he_shield_thickness, + 1 ' cm.' + write(6,*) ' ' + endif +c + call gsvolu('CRSH','TUBE', Cu, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c------------------------------------------------------ +c +c MCP detector +c + par_tube(1) = 0. + par_tube(2) = 2.5 +c + if ( l_run11 ) then + par_tube(3) = 0.15 ! Galileo Z-stack, 2.4mm totallength + ! Galileo would have 3.0mm totallength + else + par_tube(3) = 0.12 ! hamamatsu Z-stack, 2.4mm totallength + endif +c + if ( l_mcp2 ) then + call gsvolu('DMCP','TUBE', MCP, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 + endif +c +c MCP2 mounting and anode up to Run11 +c + if ( .not. l_run11) then +c +c macor ring +c + par_tube(1) = 2.5 + par_tube(2) = 3.5 + par_tube(3) = 0.55 + call gsvolu('MCPM', 'TUBE', Macor, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c anode, use MACOR +c + par_tube(1) = 0. + par_tube(2) = 3.5 + par_tube(3) = 0.05 +c + call gsvolu('MCPA', 'TUBE', Macor, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c stainless steel ring +c + par_tube(1) = 2.5 + par_tube(2) = 3.5 + par_tube(3) = 0.05 + call gsvolu('MCSR', 'TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c stainless steel support ring +c + par_tube(1) = 2.75 + par_tube(2) = 4.25 + par_tube(3) = 0.2 + call gsvolu('MCSS', 'TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c-------------------------------------------- +c + else +c +c Run 11 MCP detector +c +c ceramic rings, use MACOR as material +c + par_tube(1) = 2.4 + par_tube(2) = 3.25 + par_tube(3) = 0.075 + call gsvolu('MCPM', 'TUBE', Macor, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c anode in Run11, +c + par_box(1) = 3.65 + par_box(2) = 3.65 + par_box(3) = 0.4 + call gsvolu('MCPA','BOX ', Steel, par_box, 3, ivol) + if (ivol.lt.0) goto 900 +c +c now the anode is not completely made of steel but there +c the are "vacuum cylinders". +c + par_tube(1) = 0. + par_tube(2) = 2.75 + par_tube(3) = 0.15 +c + call gsvolu('ANVA','TUBE', Vacuum, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c stainless steel ring for MCP2 mounting +c + par_box(1) = 3.65 + par_box(2) = 3.65 + par_box(3) = 0.1 + call gsvolu('MCSR','BOX ', Steel, par_box, 3, ivol) + if (ivol.lt.0) goto 900 +c + par_tube(1) = 0 + par_tube(2) = 2.75 + par_tube(3) = 0.1 + call gsvolu('MCVR','TUBE', Vacuum, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c +c stainless steel support ring +c + par_tube(1) = 4.0 + par_tube(2) = 4.8 + par_tube(3) = 0.25 + call gsvolu('MCSS','TUBE', Steel, par_tube, 3, ivol) + if (ivol.lt.0) goto 900 +c + endif ! if .not. l_run11 +c +c------------------------------------------------------ +c +c inner scintillator +c + par_tubs(1) = 9.0 ! inner radius + par_tubs(2) = 9.5 ! outer radius + par_tubs(3) = 13.0 ! half length in z + par_tubs(4) = -45. ! starting angle + par_tubs(5) = +45. ! ending angle + call gsvolu('SCIS','TUBS', Scint, par_tubs, 5, ivol) + if (ivol.lt.0) goto 900 +c +c outer scintillator +c + par_tubs(1) = 9.6 ! inner radius + par_tubs(2) = 10.1 ! outer radius + par_tubs(3) = 13.0 ! half length in z + par_tubs(4) = -45. ! starting angle + par_tubs(5) = +45. ! ending angle + call gsvolu('SCOS','TUBS', Scint, par_tubs, 5, ivol) + if (ivol.lt.0) goto 900 +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define now the positions of vacuum chambers and detectors +c within the MARS +c +c z defines the beam direction +c + if ( l_mcps) then + call gspos('MCPS',1,'MARS',0.,0.,0.,0,'ONLY') + endif + if ( l_efla ) then + call gspos('F100',1,'MARS',0.,0.,17.2,0,'ONLY') + endif +c + call gspos('MCPV',1,'MARS',0.,0.,0.,0,'ONLY') + call gspos('F160',1,'MARS',0.,0.,-15.10,0,'ONLY') +c call gspos('T200',1,'MARS',0.,0.,-17.40,0,'ONLY') + call gspos('GATV',1,'MARS',0.,0.,-25.45,0,'ONLY') + call gspos('GATS',1,'MARS',0.,0.,-25.45,0,'ONLY') + call gspos('F200',1,'GATV',0.,0., 8.05,0,'ONLY') + call gspos('L3VA',1,'MARS',0.,0.,-57.90,0,'ONLY') + call gspos('L3ST',1,'MARS',0.,0.,-57.90,0,'ONLY') +c +c Ring anode +c + if ( l_run11) then + call gspos('RA-E',1,'MCPV',0.,0., -8.6,0,'ONLY') + call gspos('RA-M',1,'MCPV',0.,0.,-13.1,0,'ONLY') + call gspos('RA-G',1,'GATV',0.,0.,+3.35,0,'ONLY') + endif +c +c position of sample cryo things +c +c off = 0.25 + off = 0.0 + if ( l_samp ) then + call gspos('SAH3',1,'MCPV',0.,0.,+1.35,0,'ONLY') + call gspos('SAH2',1,'MCPV',0.,0.,+1.6,0,'ONLY') + call gspos('SAPH',1,'MCPV',0.,0.,+2.1,0,'ONLY') + call gspos('SAH1',1,'MCPV',0.,0.,+2.65+off,0,'ONLY') + endif +c +c off = 0.5 + off = 0.0 + if ( l_cryo ) then ! parts of the CryoVac cryostat + call gspos('COFI',1,'MCPV',0.,0.,+3.40+off,0,'ONLY') + call gspos('CRY1',1,'MCPV',0.,0.,+4.25+off,0,'ONLY') + call gspos('CRY2',1,'MCPV',0.,0.,+6.75+off,0,'ONLY') + call gspos('CRY3',1,'MCPV',0.,0.,+5.10+off,0,'ONLY') + call gspos('CRY4',1,'MCPV',0.,0.,+5.10+off,0,'ONLY') + endif +c + if ( l_crsh ) then + call gspos('CRSH',1,'MCPV',0.,0.,+1.6,0,'ONLY') + endif +c +c position of MCP +c + if ( l_mcp2 ) then + if ( l_run11) then + call gspos('DMCP',1,'MCPV',0.,0.,+1.55,0,'ONLY') ! Galileo MCP + else + call gspos('DMCP',1,'MCPV',0.,0.,+1.52,0,'ONLY') ! Hamamatsu MCP + endif + endif +c +c anode and stainless steel ring +c + if ( l_mcpa ) then +c + if ( .not. l_run11 ) then +c +c Run 10 +c + call gspos('MCPM', 1, 'MCPV', 0., 0., +1.85, 0, 'ONLY') ! MACOR ring + call gspos('MCPA', 1, 'MCPV', 0., 0., +3.45, 0, 'ONLY') ! MACOR anodee + call gspos('MCSS', 1, 'MCPV', 0., 0., +3.90, 0, 'ONLY') ! steel support ring + call gspos('MCSR', 1, 'MCPV', 0., 0., +2.45, 0, 'ONLY') ! guard ring + call gspos('MCSR', 2, 'MCPV', 0., 0., +2.95, 0, 'ONLY') ! guard ring + else +c +c Run 11 +c + call gspos('MCPM',1,'MCPV',0.,0.,+1.325,0,'ONLY') ! ceramic ring + call gspos('MCPM',2,'MCPV',0.,0.,+1.775,0,'ONLY') ! ceramic ring + call gspos('MCSR',1,'MCPV',0.,0.,+2.00,0,'ONLY') ! steel plate + call gspos('MCVR',1,'MCSR',0.,0.,+0.00,0,'ONLY') ! hole in plate + call gspos('MCPA',1,'MCPV',0.,0.,+3.10,0,'ONLY') ! the anode + call gspos('ANVA',1,'MCPA',0.,0.,-0.25,0,'ONLY') ! part of anode + call gspos('ANVA',2,'MCPA',0.,0.,+0.25,0,'ONLY') ! part of anode + call gspos('MCSS',1,'MCPV',0.,0.,+5.75,0,'ONLY') ! support ring +c + endif +c + endif +c +c position of the four scintillators +c +c-------------------------- +c +c the LEFT detector +c + call gspos('SCIS',1,'MARS',0.,0.,0.,0,'ONLY') + call gspos('SCOS',1,'MARS',0.,0.,0.,0,'ONLY') +c +c-------------------------- +c +c the TOP detector; for comments in rotation matrix, see RIGHT detector +c below +c + theta1 = 90. + phi1 = 90. + theta2 = 90. + phi2 = 180. + theta3 = 0. + phi3 = 0. +c +c GSROTM(97...) defines rotation matrix 97 +c + call gsrotm(97,theta1,phi1,theta2,phi2,theta3,phi3) +c + call gspos('SCIS',2,'MARS',0.,0.,0.,97,'ONLY') + call gspos('SCOS',2,'MARS',0.,0.,0.,97,'ONLY') +c +c-------------------------------- +c +c the RIGHT detector +c +c Now, I can get 3 copies of the SCIS and SCOS by defining a rotation +c matrix to place the other detectors in the MARS +c +c THETA1 = polar angle of new x axis +c PHI1 = azimuth of new x axis +c THETA2 = polar angle of new y axis +c PHI2 = azimuth of new y axis +c THETA3 = polar angle of new z axis +c PHI3 = azimuth of new z axis +c +c SCIS and SCOS above correspond to the LEFT detector of the LEM +c setup; now I make a copy to get the RIGHT detector +c +c z remains unchanged ==> THETA3 = 0, PHI3 = 0 +c the x axis becomes -x ==> THETA1 = 90, PHI1 = 180 +c the y axis becomes -y ==> THETA2 = 90, PHI2 = 270 +c +c the unity matrix is +c +c theta1 = 90, phi1 = 0 +c theta2 = 90, phi2 = 90 +c theta3 = 0 , phi3 = 0 +c + theta1 = 90. + phi1 = 180. + theta2 = 90. + phi2 = 270. + theta3 = 0. + phi3 = 0. + +c GSROTM(98...) defines rotation matrix 98 +c + call gsrotm(98,theta1,phi1,theta2,phi2,theta3,phi3) +c + call gspos('SCIS',3,'MARS',0.,0.,0.,98,'ONLY') + call gspos('SCOS',3,'MARS',0.,0.,0.,98,'ONLY') +c +c-------------------------- +c +c the BOTTOM detector; for comments in rotation matrix, see RIGHT detector +c above +c + theta1 = 90. + phi1 = 270. + theta2 = 90. + phi2 = 0. ! == 360 + theta3 = 0. + phi3 = 0. +c +c GSROTM(99...) defines rotation matrix 99 +c + call gsrotm(99,theta1,phi1,theta2,phi2,theta3,phi3) +c + call gspos('SCIS',4,'MARS',0.,0.,0.,99,'ONLY') + call gspos('SCOS',4,'MARS',0.,0.,0.,99,'ONLY') +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define basic detector parameter and allocation of the detector +c volume to a SET using GSDET; a SET becomes useful in large +c detectors where it is possible to group together several detetectors +c like the different components of a calorimeter; +c in this simple case it isn't necessary, but GEANT aquires it +c + namesv(1) = 'DMCP' +c +c the following variables have been "included" +c + nvname = 1 + numbv(1) = 1 + idtype = 1 ! user defined detector ID + if ( l_mcp2 ) then + call gsdet('MDET','DMCP',nvname,namesv,numbv,idtype,100, + 1 100,iset,idet) +c +c define now the hit parameters using GSDETH +c MCP : the interesting variables I want to store are +c +c - position x,y,z +c - time t +c - energy e +c - energyloss eloss +c - particle type ipart +c + call gsdeth('MDET','DMCP',7,namesh,nbitsh,orig,fact) + endif +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c define Scintillator detector +c + nvname = 1 + namesv(1) = 'SCIS' + numbv(1) = 3 + idtype = 2 ! user defined detector ID + call gsdet('ISSC','SCIS',nvname,namesv,numbv,idtype,100,100, + 1 iset,idet) + call gsdeth('ISSC','SCIS',7,namesh,nbitsh,orig,fact) + nvname = 1 + namesv(1) = 'SCOS' + numbv(1) = 3 + idtype = 3 ! user defined detector ID +c + call gsdet('OSSC','SCOS',nvname,namesv,numbv,idtype,100,100, + 1 iset, idet) + call gsdeth('OSSC','SCOS',7,namesh,nbitsh,orig,fact) +c +c call gpsets('MDET','DMCP') +c call gpsets('ISSC','SCIS') +c call gpsets('OSSC','SCOS') +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + return +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c error handling +c +900 write(6,'(''0GEANT_LEMSR-E- Error in routine GSVOLU...'',/)') + call exit +c + end