* * $Id$ * * $Log$ * Revision 1.22 2005/12/27 01:14:28 jonesrt * gustep.F * - removed FDC cathode strips from the list of sensitive volumes where * hits can occur -- the cathode strips themselves are not sensitive, * they just pick up signals from the anode wires, as in the statement * no anode hit => no cathode hit! * --rtj-- * * Revision 1.21 2005/12/12 15:38:46 jonesrt * hitutil.F * - removed from top-level, moved into new subdir hitutil [rtj] * hitBCal.c * - increased the segmentation of the BCal readout [rtj] * cdcdump.c * - added a check for the existence of certain groups [rtj] * gustep.F * - added dispatch for additional BCal readout sections [rtj] * Makefile.orig * - removed compilation of hitutil.F, added libhitutil.a linkage [rtj] * * Revision 1.20 2005/06/22 15:32:16 zisis * * First draft of this routine to work with the new BCAL geometry and hits. CX. * * Revision 1.19 2005/04/26 14:51:52 ostrov * Addition of UPV * * Revision 1.19 2005/03/20 ao * hitUPV.c * Makefile.orig * gustep.F * - support for hits in UPV * * Revision 1.18 2005/01/29 19:25:28 jonesrt * hitLGD.c * - renamed to hitFCal.c [rtj] * Makefile.orig * - modified to reflect the name change for hitLGD.c [rtj] * control.in * - it seems that I always have touched this file at some point! [rtj] * gustep.F * - added argument ISTAK to argument list for hitXXX functions, so that * they can determine whether a given track is the primary or not [rtj] * - some of the names of volumes have been changed in the recent geometry * update, reflect that fact in the Makefile [rtj] * hddm_s.c, hddm_s.h * - updated from hddm (you should generate these using hddm-c and then * copy them over from hddm to this folder) [rtj] * hddsGeant3.F * - updated from hdds (you should generate this using hdds-geant and * then copy it over from hdds to this folder) [rtj] * uginit.F * - added a line to switch HBOOK from //LUN3 (closed after return from * GRFILE) back to the geant.hbook output file on unit 50 [rtj] * hddmInput.c * - modified to store the actual coordinates of the event vertex in the * Monte Carlo section of the output record, in case that the vertex * was generated by the simulator instead of the generator [rtj] * hitXXX.c * - modified to accommodate an extra tag primary="boolean" on all of * the cheat tags, to tell whether the hit was produced by one of the * original primaries, or by a secondary produced by one of them. * - hitFTOF.c modified to accommodate two layers instead of one [rtj] * - hitStart.c modified to accommodate the segmented readout [rtj] * - hitCerenkov.c - modified to accommodate the segmented readout [rtj] * - hitCerenkov.c - added a cheat tag to the Cerenkov readout [rtj] * - all cheat tags have been modified to report all three coordinates * (in the global reference system) instead of only two [rtj] * * Revision 1.17 2005/01/21 09:34:05 davidl * If ff card NOSECONDARIES set, then don't push any secondaries onto the stack * * Revision 1.16 2004/06/17 18:32:55 davidl * Fixed typo that caused comment to be misleading * * Revision 1.15 2004/06/07 19:05:53 jonesrt * Makefile, gustep.F * - added option CERENKOV_PID_NTUPLE to save information from a bg * simulation to an ntuple stored in geant.hbook [rjt,rem] * * Revision 1.14 2004/05/18 12:58:54 jonesrt * Makefile * - created a section at the top for global defines that are used to * build custom versions of the simulation [rtj] * hddm_s.c, hddm_s.h * - default i/o library modules (generated by hddm package) [rtj] * hddsGeant3.F * - default geometry module (generated by hdds package) [rtj] * hitStart.c * - changes to accomodate new vertex counter cylinder+plane structure [rtj] * guout.F, gustep.F * - defines for custom builds moved from sources to Makefile * gustep.F * - new conditional WERNERS_VTX_NTUPLE sections added for background * studies in the region of the vertex counter [rtj] * * Revision 1.13 2004/01/14 16:34:48 brash * Fixed bug in gustep.F regarding placement of certain assignment statements with * respect to certain #ifdef statements. Should work now with BACKGROUND_* defines * turned either on or off. (EJB) * * Revision 1.12 2004/01/14 16:28:10 brash * Updates in order to analyze different readout modules of the barrel calorimeter. (EJB) * * Revision 1.11 2003/12/10 15:32:57 jonesrt * -control.in : never mind [rtj] * -gukine.F : fixed a bug in the setting of polarization ppol used * for polarization studies of the coherent bremsstrahlung source [rtj] * -gustep.F : changed background studies facility to split data across * separate ntuples, one for each virtual detector [rtj] * -hit*.F : modified behaviour from "quit" to "truncate" in the case where * the number of hits exceeds the maximum allowed for that counter [rtj] * * Revision 1.10 2003/01/08 19:17:34 jonesrt * - gustep.F : added collection of information in a backgrounds ntuple - rtj * - gufld.F : enabled magnetic field in sweep magnets, was off before - rtj * * Revision 1.9 2003/01/02 23:49:33 jonesrt * - included updates in gustep.F with conditional code for background * studies, by R. Jones and C. Gauthier [rtj] * - moved the beam origin a meter upstream in gukine.F to make room for * additional shielding in the collimator cave [rtj] * * Revision 1.8 2001/12/19 02:34:55 jonesrt * Fixed the names of sensitive volumes in the save-hits section of gustep.F, * also added support for the r="float" parameter of in hitCDC.c. * -rtj- * * Revision 1.7 2001/10/29 17:39:23 jonesrt * - added mc truth info to output event for internal track/photon generators * - added special code for background studies, selected by the conditional * #define BACKGROUND_STUDIES (in gustep.F) * - added conditional code to disable normal event output for bg studies, using * #define DISABLE_OUTPUT (in guout.F) * Both of the above defines are disabled in the distribution code by default. * -rtj- * * Revision 1.6 2001/07/27 21:04:09 jonesrt * With this release, HDGeant version 1.0 is now in beta. -rtj * * Revision 1.5 2001/07/24 05:37:16 jonesrt * First working prototype of hits package -rtj * * Revision 1.4 2001/07/19 23:25:49 jonesrt * numerous new files as I develop the prototype hits libraries -rtj * * Revision 1.3 2001/07/17 22:38:40 jonesrt * Adding hits registry in gustep -rtj * * Revision 1.2 2001/07/15 07:31:37 jonesrt * HDGeant now supports kinematic input from Monte Carlo generators * via the routines in hddmInput.c -rtj * * Revision 1.1 2001/07/10 18:05:47 jonesrt * imported several of the gu*.F user subroutines for Hall D customization -rtj * * Revision 1.1.1.1 1995/10/24 10:21:52 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.35 by S.Giani *-- Author : SUBROUTINE GUSTEP * ************************************************************************ * * * User routine called at the end of each tracking step * * MEC is the mechanism origin of the step * * INWVOL is different from 0 when the track has reached * * a volume boundary * * ISTOP is different from 0 if the track has stopped * * * ************************************************************************ * #include "geant321/gckine.inc" #include "geant321/gcking.inc" #include "geant321/gcomis.inc" #include "geant321/gcvolu.inc" #include "geant321/gctmed.inc" #include "geant321/gctrak.inc" #include "geant321/gcflag.inc" #include "geant321/gcphys.inc" #include "geant321/gcbank.inc" #include "geant321/gcmate.inc" #include "hdtrackparams.inc" #include "controlparams.inc" * * ----------------------------------------------------------------- * #define TOP_CERENKOV_EFFICIENCY 0.25 character*80 title character*4 cnames(15) equivalence (NAMES(1),cnames(1)) real rnd(100) integer*4 cint character*4 cchar equivalence(cint,cchar) #ifdef BACKGROUND_STUDIES integer type real xv(4),Etot common /bgNtuple/type,xv,Etot character*80 bgntdef parameter (bgntdef='type:I,xv(4):R,Etot:R') integer bgnt parameter (bgnt=10) #endif #ifdef BACKGROUND_PROFILING integer det real vertx(3),tofgx,ubuf(99),xint(3,999) integer ntbeamx,nttargx,nubuf,mint save ubuf,xint,det,mint integer iorder(999) save iorder #endif #ifdef ACTIVE_COLLIMATOR_SIMS integer evno real vertx1(3),tofgx1 integer ntbeamx1,nttargx1,nubuf1 real qsum(8),qsumb(8),qsump(8),ubuf1(99) common /hdacol/evno,Egam,qsum,qsumb,qsump,ubuf1 save /hdacol/ character*100 ntacoldef parameter (ntacoldef='ev:I,Egam:R,qsum(8):R,qsumb(8):R,' + //'qsump(8):R,ppol:R,offset(2):R') integer ntacol parameter (ntacol=25) #endif #if defined WERNERS_VTX_NTUPLE integer evno,part,vid real xvtx(3),xdet(3),pdet(3) common /wernerNtuple/evno,xvtx,xdet,pdet,part,vid save /wernerNtuple/ character*80 ntwernerdef parameter (ntwernerdef='ev:I,xv(3):R,xt(3):R,p(3):R,part:I,vid:I') integer ntwerner parameter (ntwerner=10) #elif defined CERENKOV_PID_NTUPLE integer evno,part real xvtx(3),xdet(3),pdet(3) common /ckovNtuple/evno,xvtx,xdet,pdet,part save /ckovNtuple/ character*80 ntckovdef parameter (ntckovdef='ev:I,xv(3):R,xt(3):R,p(3):R,part:I,vid:I') integer ntckov parameter (ntckov=10) #elif defined FCAL_SPLASH_NTUPLE integer evno,part real xconv(6),pconv integer nfcal real xfcal(999),yfcal(999),Efcal(999) common /fsplNtuple/evno,part,xconv,pconv,nfcal,xfcal,yfcal,Efcal save /fsplNtuple/ character*160 ntfspldef parameter (ntfspldef='ev:I,part:I,xconv(6):R,pconv:R,' + //'nfcal[0,999]:I,xfcal(nfcal):R,yfcal(nfcal):R,Efcal(nfcal):R') integer ntfspl parameter (ntfspl=10) integer getrow,getcolumn external getrow,getcolumn #elif defined FCAL_SHOWER_PROFILE_NTUPLE c c NTUPLE CONTENTS: c evno: event number c part: initial particle type c xinit(3): initial particle location c pinit(3): initial particle momentum c xstop(3): location where initial particle stops c ndiv: number of fcal divisions in the following arrays... c xdiv(ndiv): array of x positions c ydiv(ndiv): array of y positions c zdiv(ndiv): array of z positions c ediv(ndiv): array of energy depositions at the above positions c #define NXDIVISIONS 400 #define NYDIVISIONS 400 #define NZDIVISIONS 60 #define XDIVISIONSIZE 1.0 #define YDIVISIONSIZE 1.0 #define ZDIVISIONSIZE 1.0 #define XDIVISIONSTART -200.0 #define YDIVISIONSTART -200.0 #define ZDIVISIONSTART 620.0 #define MINIMUMENERGYPERDIVISION 1.0e-3 c real edivisions(NXDIVISIONS,NYDIVISIONS,NZDIVISIONS) save edivisions integer ixdivision, iydivision, izdivision integer evno,part real xstop(3),xinit(3),pinit(3) integer ndiv real xdiv(10000),ydiv(10000),zdiv(10000),ediv(10000) common /fspNtuple/evno,part,xinit,pinit,xstop, + ndiv,xdiv,ydiv,zdiv,ediv save /fspNtuple/ character*160 ntfspdef parameter (ntfspdef='evno:I,part:I,xinit(3):R,pinit(3):R,' + //'xstop(3):R,ndiv[0,10000]:I,xdiv(ndiv):R,ydiv(ndiv):R,' + //'zdiv(ndiv):R,ediv(ndiv):R') integer ntfsp parameter (ntfsp=10) c #elif defined BCAL_RESOLUTION_NTUPLE integer evno real pphot(4),dE,xdep(3) common /bcalresnt/evno,pphot,dE,xdep save /bcalresnt/ character*80 ntbcaldef parameter (ntbcaldef='evno:I,pphot(4):R,dE:R,xdep(3):R') integer ntbcal parameter (ntbcal=33) #endif logical hexist external hexist #ifdef BACKGROUND_PROFILING ******************************************************************************** * The following defines an ntuple containing information on particle type * energy, position, polarization, and at what virtual detector the * particle passes through. the last colume entry is defined by a integer * 'det' whose value is the number of the detector the particle passes * through. The xint(icomp,iint) array records the vertex history of the * interaction sequence leading to the detected particle, iint=1...mint. ******************************************************************************** ********** assignment of the value of 'det' to a particle ****************** if (cnames(nlevel).eq.'DET1') then det = 1 elseif (cnames(nlevel).eq.'DET2') then det = 2 elseif (cnames(nlevel).eq.'DET3') then det = 3 elseif (cnames(nlevel).eq.'DET4') then det = 4 elseif (cnames(nlevel).eq.'DET5') then det = 5 elseif (cnames(nlevel).eq.'DET6') then det = 6 elseif (cnames(nlevel).eq.'DET7') then det = 7 elseif (cnames(nlevel).eq.'CONV') then det = 8 elseif (cnames(nlevel).eq.'LIH2') then det = 9 else det = 0 endif ********************** Defintion of ntuple **************************** if (inwvol.eq.1.and.det.ge.1 c + .and.((ipart.eq.2).or.(ipart.eq.3)) c + .and.(vect(1)**2+vect(2)**2).gt.25.0) then + ) then call gfvert(1,vertx,ntbeamx,nttargx,tofgx,ubuf,nubuf) if (.not.hexist(9+det)) then write(title,"('hits in virtual detector',i3)") det call hbnt(9+det,title,' ') call hbname(9+det,'hits',gekin,'totE:r') call hbname(9+det,'hits',vect,'x(7):r') call hbname(9+det,'hits',ubuf,'ppol:r') call hbname(9+det,'hits',ubuf(2),'xspot(2):r') call hbname(9+det,'hits',ipart,'ptype:i') call hbname(9+det,'hits',det,'det:i') call hbname(9+det,'hits',mint,'mint[0,999]:i') call hbname(9+det,'hits',xint,'xint(3,mint):r') endif call hfnt(9+det) endif if (nstep.eq.0) then if (istak.eq.0) then mint = 0 else do while (mint.gt.0.and.iorder(mint).gt.istak) mint = mint-1 enddo if (mint.lt.999) then mint = mint+1 iorder(mint) = istak xint(1,mint) = vect(1) xint(2,mint) = vect(2) xint(3,mint) = vect(3) endif endif endif ******************************************************************************** #endif * Stop beam photons which are initially generated outside the 12.4 cm * radius beam pipe at the entrance to the collimator cave IF(vect(3).eq.-2200.) THEN ph_radius = sqrt(vect(1)*vect(1)+vect(2)*vect(2)) IF(ph_radius.ge.12.4) THEN ISTOP = 1 ENDIF ENDIF * Stop wimpy charged particles that are taking forever to range out if ((NSTEP.ge.99999).and.(CHARGE.ne.0).or. + (ILOSS.eq.0).and.(GEKIN.lt.1e-3)) then DESTEP = GEKIN ISTOP = 98 NMEC = NMEC+1 LMEC(NMEC) = 30 endif * Place any secondaries generated during this step onto the stack if (nosecondaries.eq.0) then do i=1,NGKINE itypa = GKIN(5,i) if (itypa.ne.4) then c make primary except if in calorimeter volume and not hadronic interaction iflgk(i) = 1 cint = KCASE c print *, cchar rx = sqrt(VERT(1)**2+VERT(2)**2) if ( (((rx>65.) .and.((VERT(3)>-17.) .and. > (VERT(3)<390.))).or.(VERT(3)>625.)).and. > (cchar.ne.'HADR')) then iflgk(i) = 0 endif call GSKING(i) endif enddo endif * For explicit Cerenkov generation, apply an inefficiency factor if (NGPHOT.gt.100) then call GRANOR(rnd(1),rnd(2)) sigma = sqrt(NGPHOT * (TOP_CERENKOV_EFFICIENCY + * (1 - TOP_CERENKOV_EFFICIENCY))) NGPHOT = TOP_CERENKOV_EFFICIENCY*NGPHOT + rnd(1)*sigma + 0.5 call GSKPHO(0) elseif (NGPHOT.gt.0) then call GRNDM(rnd,NGPHOT) do i=1,NGPHOT if (rnd(i).le.TOP_CERENKOV_EFFICIENCY) then call GSKPHO(i) endif enddo endif #ifdef ACTIVE_COLLIMATOR_SIMS if (IDEVT.ne.evno) then if (IDEVT.eq.1) then call hbnt(ntacol,'active collimator ntuple, all events',' ') call hbname(ntacol,'charge',evno,ntacoldef) call hbnt(ntacol+1,'active collimator ntuple, charged events', + ' ') call hbname(ntacol+1,'charge',evno,ntacoldef) else call hfnt(ntacol) if (qsum(1).ne.0 .or. qsum(2).ne.0 .or. + qsum(3).ne.0 .or. qsum(4).ne.0 .or. + qsum(5).ne.0 .or. qsum(6).ne.0 .or. + qsum(7).ne.0 .or. qsum(8).ne.0 .or. + qsumb(1).ne.0 .or. qsumb(2).ne.0 .or. + qsumb(3).ne.0 .or. qsumb(4).ne.0 .or. + qsumb(5).ne.0 .or. qsumb(6).ne.0 .or. + qsumb(7).ne.0 .or. qsumb(8).ne.0 .or. + qsump(1).ne.0 .or. qsump(2).ne.0 .or. + qsump(3).ne.0 .or. qsump(4).ne.0 .or. + qsump(5).ne.0 .or. qsump(6).ne.0 .or. + qsump(7).ne.0 .or. qsump(8).ne.0) then call hfnt(ntacol+1) endif endif evno = IDEVT Egam = GETOT do i=1,8 qsum(i) = 0 qsumb(i) = 0 qsump(i) = 0 enddo call gfvert(1,vertx1,ntbeamx1,nttargx1,tofgx1,ubuf1,nubuf1) endif if (cnames(NLEVEL).eq.'PCTT' .or. + cnames(NLEVEL).eq.'PCPB') then ISTOP = 3 NMEC = NMEC+1 LMEC(NMEC) = 30 elseif (NLEVEL.ge.7.and.cnames(7).eq.'ACWI') then if (INWVOL.eq.2) then qsum(NUMBER(7)) = qsum(NUMBER(7))-CHARGE if (cnames(NLEVEL).eq.'ACBI') then c print *, '*-* leaving ACBI' qsumb(NUMBER(7)) = qsumb(NUMBER(7))-CHARGE elseif (cnames(NLEVEL).eq.'PIN1') then c print *, '*-* leaving ACWI/PIN1' qsump(NUMBER(7)) = qsump(NUMBER(7))-CHARGE else c print *, '*-* leaving ACWI' endif elseif ((INWVOL.eq.1).and.(NSTEP.gt.0)) then qsum(NUMBER(7)) = qsum(NUMBER(7))+CHARGE if (cnames(NLEVEL).eq.'ACBI') then c print *, '*-* entering ACBI' qsumb(NUMBER(7)) = qsumb(NUMBER(7))+CHARGE elseif (cnames(NLEVEL).eq.'PIN1') then c print *, '*-* entering ACWI/PIN1' qsump(NUMBER(7)) = qsump(NUMBER(7))+CHARGE else c print *, '*-* entering ACWI' endif endif elseif (NLEVEL.ge.7.and.cnames(7).eq.'ACWO') then if (INWVOL.eq.2) then qsum(NUMBER(7)+4) = qsum(NUMBER(7)+4)-CHARGE if (cnames(NLEVEL).eq.'ACBO') then c print *, '*-* leaving ACBO' qsumb(NUMBER(7)+4) = qsumb(NUMBER(7)+4)-CHARGE elseif (cnames(NLEVEL).eq.'PIN1') then c print *, '*-* leaving ACWO/PIN1' qsump(NUMBER(7)+4) = qsump(NUMBER(7)+4)-CHARGE else c print *, '*-* leaving ACWO' endif elseif ((INWVOL.eq.1).and.(NSTEP.gt.0)) then qsum(NUMBER(7)+4) = qsum(NUMBER(7)+4)+CHARGE if (cnames(NLEVEL).eq.'ACBO') then c print *, '*-* entering ACBO' qsumb(NUMBER(7)+4) = qsumb(NUMBER(7)+4)+CHARGE elseif (cnames(NLEVEL).eq.'PIN1') then c print *, '*-* entering ACWO/PIN1' qsump(NUMBER(7)+4) = qsump(NUMBER(7)+4)+CHARGE else c print *, '*-* entering ACWO' endif endif endif #endif #ifdef BACKGROUND_STUDIES if (.not.HEXIST(bgnt)) then call HBNT(bgnt,'background particles','') call HBNAME(bgnt,'tracks',type,bgntdef) endif if (INWVOL.eq.1) then z = VECT(3) r = sqrt(VECT(1)**2 + VECT(2)**2) if (cnames(NLEVEL).eq.'PCTT' .or. + cnames(NLEVEL).eq.'PCPB') then ISTOP = 3 NMEC = NMEC+1 LMEC(NMEC) = 30 elseif ((cnames(NLEVEL).eq.'UWIT' .and. abs(z-0.02).lt.0.1) + .or. (cnames(NLEVEL).eq.'VRTX' .and. abs(r-4.95).lt.0.1) + .or. (cnames(NLEVEL).eq.'CDCI' .and. abs(r-15.0).lt.0.1) + .or. (cnames(NLEVEL).eq.'DC12' .and. abs(r-37.0).lt.0.1) + .or. (cnames(NLEVEL).eq.'CDCO' .and. abs(r-59.0).lt.0.1) + .or. (cnames(NLEVEL).eq.'BCAL' .and. abs(r-65.0).lt.0.1) + .or. (cnames(NLEVEL).eq.'FDC ' .and. abs(z-224.).lt.1.0) + .or. (cnames(NLEVEL).eq.'CERE' .and. abs(z-410.).lt.1.0) + .or. (cnames(NLEVEL).eq.'FCAL' .and. abs(z-575.).lt.1.0)) + then xv(1) = VECT(1) xv(2) = VECT(2) xv(3) = VECT(3) xv(4) = r type = IPART Etot = GETOT call HFNT(bgnt) endif endif #endif #if defined WERNERS_VTX_NTUPLE if (.not.HEXIST(ntwerner)) then call HBNT(ntwerner,'vertex counter hits','') call HBNAME(ntwerner,'hits',evno,ntwernerdef) endif if (NSTEP.eq.0) then xvtx(1) = vect(1) xvtx(2) = vect(2) xvtx(3) = vect(3) elseif (INWVOL.eq.1) then if (cnames(NLEVEL).eq.'STRC') then vid = 1 elseif (cnames(NLEVEL).eq.'STRP') then vid = 2 else vid = 0 endif if (vid.gt.0) then xdet(1) = VECT(1) xdet(2) = VECT(2) xdet(3) = VECT(3) pdet(1) = VECT(4)*VECT(7) pdet(2) = VECT(5)*VECT(7) pdet(3) = VECT(6)*VECT(7) part = IPART evno = IDEVT call HFNT(ntwerner) endif endif #elif defined CERENKOV_PID_NTUPLE if (.not.HEXIST(ntckov)) then call HBNT(ntckov,'cerenkov counter hits','') call HBNAME(ntckov,'hits',evno,ntckovdef) endif if (NSTEP.eq.0) then xvtx(1) = vect(1) xvtx(2) = vect(2) xvtx(3) = vect(3) elseif (INWVOL.eq.1) then if (cnames(NLEVEL).eq.'CGAS') then xdet(1) = VECT(1) xdet(2) = VECT(2) xdet(3) = VECT(3) pdet(1) = VECT(4)*VECT(7) pdet(2) = VECT(5)*VECT(7) pdet(3) = VECT(6)*VECT(7) part = IPART evno = IDEVT call HFNT(ntckov) endif endif #elif defined FCAL_SPLASH_NTUPLE if (.not.HEXIST(ntfspl)) then call HBNT(ntfspl,'FCal splash hits','') call HBNAME(ntfspl,'hits',evno,ntfspldef) endif if (ISTAK.eq.0.and.NSTEP.eq.0) then part = IPART evno = IDEVT nfcal = 0 elseif (ISTAK.eq.0.and.ISTOP.ne.0) then xconv(1) = vect(1) xconv(2) = vect(2) xconv(3) = vect(3) xconv(4) = vect(4) xconv(5) = vect(5) xconv(6) = vect(6) pconv = vect(7) elseif (cnames(NLEVEL).eq.'LGBL'.and.DESTEP.gt.1e-3) then row = getrow() col = getcolumn() do i=1,nfcal if (col.eq.xfcal(i).and.row.eq.yfcal(i)) then goto 810 endif enddo nfcal = nfcal+1 i = nfcal xfcal(i) = col yfcal(i) = row Efcal(i) = 0 810 continue Efcal(i) = Efcal(i)+DESTEP endif if (ISTOP.ne.0.and.IQ(JSTAK+1).eq.0) then call HFNT(ntfspl) endif #elif defined FCAL_SHOWER_PROFILE_NTUPLE c c create ntuple c if (.not.HEXIST(ntfsp)) then call HBNT(ntfsp,'FCal shower profile','') call HBNAME(ntfsp,'hits',evno,ntfspdef) endif c c get information about primary and c reset division variables c if (ISTAK.eq.0.and.NSTEP.eq.0) then evno = IDEVT part = IPART ndiv = 0 do ix = 1, NXDIVISIONS do iy = 1, NYDIVISIONS do iz = 1, NZDIVISIONS edivisions(ix,iy,iz) = 0.0 enddo enddo enddo xinit(1) = vect(1) xinit(2) = vect(2) xinit(3) = vect(3) pinit(1) = vect(4)*vect(7) pinit(2) = vect(5)*vect(7) pinit(3) = vect(6)*vect(7) c c find the location where the primary stops c elseif (ISTAK.eq.0.and.ISTOP.ne.0) then xstop(1) = vect(1) xstop(2) = vect(2) xstop(3) = vect(3) c c record energy depositions in fcal divisions c elseif (cnames(NLEVEL).eq.'LGBL'.and.DESTEP.gt.0.0) then ixdivision = int((vect(1)-XDIVISIONSTART)/XDIVISIONSIZE) + 1 iydivision = int((vect(2)-YDIVISIONSTART)/YDIVISIONSIZE) + 1 izdivision = int((vect(3)-ZDIVISIONSTART)/ZDIVISIONSIZE) + 1 if ((ixdivision.ge.1).and.(ixdivision.le.NXDIVISIONS).and. + (iydivision.ge.1).and.(iydivision.le.NYDIVISIONS).and. + (izdivision.ge.1).and.(izdivision.le.NZDIVISIONS)) then edivisions(ixdivision,iydivision,izdivision) = + edivisions(ixdivision,iydivision,izdivision) + DESTEP endif endif c c rearrange division information and record ntuples c if (ISTOP.ne.0.and.IQ(JSTAK+1).eq.0) then do ix = 1, NXDIVISIONS do iy = 1, NYDIVISIONS do iz = 1, NZDIVISIONS if (edivisions(ix,iy,iz).ge.MINIMUMENERGYPERDIVISION) then if (ndiv.lt.10000) then ndiv = ndiv + 1 xdiv(ndiv) = ix*XDIVISIONSIZE + XDIVISIONSTART + - XDIVISIONSIZE/2.0 ydiv(ndiv) = iy*YDIVISIONSIZE + YDIVISIONSTART + - YDIVISIONSIZE/2.0 zdiv(ndiv) = iz*ZDIVISIONSIZE + ZDIVISIONSTART + - ZDIVISIONSIZE/2.0 ediv(ndiv) = edivisions(ix,iy,iz) endif endif enddo enddo enddo call HFNT(ntfsp) endif c c #elif defined FCAL_RADIATION_HISTOGRAMS c c book histograms c if (.not.HEXIST(1000)) then call HBOOK2(1000,'rad',240,-120.,120.,240,-120.,120.,0) do iz = 1, 60 call HBOOK2(1620+iz-1,'rad',240,-120.,120.,240,-120.,120.,0) enddo endif c c record energy depositions c if (cnames(NLEVEL).eq.'LGBL'.and.DESTEP.gt.0.0) then call HFILL(1000, vect(1),vect(2),DESTEP) call HFILL(1000+int(vect(3)),vect(1),vect(2),DESTEP) endif c c #elif defined BCAL_RESOLUTION_NTUPLE if (.not.HEXIST(ntbcal)) then call HBNT(ntbcal,'BCal shower energy deposition','') call HBNAME(ntbcal,'deposits',evno,ntbcaldef) endif if (ITRA.eq.1.and.ISTAK.eq.0.and.NSTEP.eq.0) then evno=IDEVT pphot(1) = VECT(4)*VECT(7) pphot(2) = VECT(5)*VECT(7) pphot(3) = VECT(6)*VECT(7) pphot(4) = VECT(7) elseif (DESTEP.gt.0.and. + NLEVEL.gt.6.and.cnames(6).eq.'BCAM') then dE=DESTEP xdep(1)=vect(1) xdep(2)=vect(2) xdep(3)=vect(3) call HFNT(ntbcal) endif #endif * Optionally store particle trajectory if (storetraj.ne.0) then call addtrajectorypoint(VECT,TOFG,DESTEP,GETOT,ITRA,ISTAK + ,IPART, RADL, STEP, NMEC, LMEC, storetraj) endif c #ifndef TRACK_SHOWERS_IN_COLLIMATOR if(showersincol.eq.0) then if ((cnames(NLEVEL).eq.'INSU').or. + (cnames(NLEVEL).eq.'PCTT').or. + (cnames(NLEVEL).eq.'PCPB')) then ISTOP = 3 NMEC = NMEC+1 LMEC(NMEC) = 30 endif endif * #endif #ifdef HISTOGRAM_MATERIAL_SEEN_BY_FIRST_TRACK if (ITRA.eq.1.and.ISTAK.eq.0) then if (NSTEP.eq.0) then title = 'material in g/cm2 seen by first track' call HBOOK1(IDEVT,title,1000,0.,200.,0) title = 'radiation lengths seen by first track' call HBOOK1(IDEVT+1000000,title,700,0.,700.,0) endif call HFILL(IDEVT,SLENG,0.,STEP*DENS) call HFILL(IDEVT+1000000,VECT(3),0.,STEP/RADL) endif #endif * Register hits in sensitive detectors here if (ISVOL.ne.0.or.ISTOP.ne.0.or.INWVOL.ge.2) then call savehits endif CALL GDEBUG END