* * $Id: gustep.F 4225 2008-09-17 22:27:02Z jonesrt $ * * $Log$ * Revision 1.23 2008/12/16 19:36:04 somov * gustep.F_review2008 * - Stop beam photons which are initially generated * outside the 12.4 cm radius beam pipe at the entrance * to the collimator cave. * - The primary collimator volume COL1 was replaced by * volumes PCTT, PCTB, and PCSD (these volumes are used * if one wants to stop showers in the 1st collimator). * * 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" * * ----------------------------------------------------------------- * #define TOP_CERENKOV_EFFICIENCY 0.25 character*80 title character*4 cnames(15) equivalence (NAMES(1),cnames(1)) real rnd(100) #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 #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 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 real ph_radius #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 else det = 0 endif ********************** Defintion of ntuple **************************** if (inwvol.eq.1.and.det.gt.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',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 CALL GDEBUG * Register hits in sensitive detectors here if (ISVOL.ne.0) then call savehits 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 * 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) call GSKING(i) 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 * Stop wimpy charged particles that are taking forever to range out if ((NSTEP.ge.9999).and.(CHARGE.ne.0).or. + (ILOSS.eq.0).and.(GEKIN.lt.1e-3)) then DESTEP = GEKIN ISTOP = 1 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.'PCTB').or. + (cnames(NLEVEL).eq.'PCSD')) then ISTOP = 1 return 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 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 #ifndef TRACK_SHOWERS_IN_COLLIMATOR if ((cnames(NLEVEL).eq.'INSU').or. + (cnames(NLEVEL).eq.'PCTT').or. + (cnames(NLEVEL).eq.'PCTB').or. + (cnames(NLEVEL).eq.'PCSD')) then ISTOP=1 endif if (cnames(NLEVEL).eq.'OCOL') then ISTOP=1 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 END