TFluka/source_ray.cxx

   1 // Fortran
   2 #include "TCallf77.h"
   3
   4 // Fluka commons
   5 #include "Fdblprc.h"  //(DBLPRC) fluka common
   6 #include "Fdimpar.h"  //(DIMPAR) fluka parameters
   7 #include "Fepisor.h"  //(EPISOR) fluka common
   8 #include "Fstack.h"   //(STACK)  fluka common
   9 #include "Fstars.h"   //(STARS)  fluka common
  10 #include "Fbeam.h"    //(BEAM)   fluka common
  11 #include "Fpaprop.h"  //(PAPROP) fluka common
  12 #include "Fltclcm.h"  //(LTCLCM) fluka common
  13 #include "Fpart.h"
  14 //#include "Fcaslim.h"  //(CASLIM) fluka common
  15
  16 #include "TVector3.h"
  17 #include "TRandom.h"
  18 #include "TParticle.h"
  19 #include "../EVGEN/AliGenScan.h"
  20 #include "AliStack.h"
  21
  22 //Other
  23 #include <Riostream.h>
  24
  25 #ifndef WIN32
  26 # define source source_
  27 # define geocrs geocrs_
  28 # define georeg georeg_
  29 # define geohsm geohsm_
  30 # define soevsv soevsv_
  31 # define mcihad mcihad_
  32 # define source_ray source_ray__
  33 #else
  34 # define source SOURCE
  35 # define geocrs GEOCRS
  36 # define georeg GEOREG
  37 # define geohsm GEOHSM
  38 # define soevsv SOEVSV
  39 # define mcihad MCIHAD
  40 # define source_ray SOURCE_RAY
  41 #endif
  42
  43 extern "C" {
  44   //
  45   // Prototypes for FLUKA functions
  46   //
  47   void type_of_call geocrs(Double_t &, Double_t &, Double_t &);
  48   void type_of_call georeg(Double_t &, Double_t &, Double_t &,
  49                            Int_t &, Int_t &);
  50   void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
  51   void type_of_call soevsv();
  52   int  type_of_call mcihad(const int&);
  53  /*
  54    *----------------------------------------------------------------------*
  55    *                                                                      *
  56    *     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
  57    *                                                   Infn - Milan       *
  58    *                                                                      *
  59    *     Last change on 21-jun-98     by    Alfredo Ferrari               *
  60    *                                                                      *
  61    *     C++ version on 27-sep-02     by    Isidro Gonzalez               *
  62    *                                                                      *
  63    *  This is just an example of a possible user written source routine.  *
  64    *  note that the beam card still has some meaning - in the scoring the *
  65    *  maximum momentum used in deciding the binning is taken from the     *
  66    *  beam momentum.  Other beam card parameters are obsolete.            *
  67    *                                                                      *
  68    *----------------------------------------------------------------------*/
  69
  70   void source_ray(Int_t& nomore) {
  71
  72       static Bool_t lfirst       = true;
  73       static AliGenScan* gener   = 0;
  74       static AliStack* stack     = 0;
  75
  76       nomore = 0;
  77       TParticle* particle = 0;
  78       Int_t itrack        = -1;
  79       if (lfirst) {
  80           EPISOR.tkesum = zerzer;
  81           lfirst = false;
  82           EPISOR.lussrc = true;
  83 //
  84 // The generator
  85 //
  86           gener  = new AliGenScan(0);
  87           gener->SetRange(40, -200., 200., 40, -200., 200., 1, -700., -699.);
  88           gener->SetPart(0);
  89           gener->SetThetaRange(180., 180.);
  90           gener->SetPhiRange(0.,0.);
  91
  92
  93 // Stack
  94           stack = new AliStack(10000);
  95           gener->SetStack(stack);
  96           gener->Init();
  97
  98       } else {
  99           //
 100           // Generate event
 101           stack->Reset();
 102           gener->Generate();
 103           Int_t npart = stack->GetNprimary();
 104           for (Int_t part = 0; part < npart; part++) {
 105               STACK.lstack++;
 106               particle = stack->Particle(part);
 107
 108               printf("Particle %5d %10s %10.3f %10.3f %10.3f \n",
 109                      STACK.lstack, particle->GetName(),
 110                      particle->Px(), particle->Py(), particle->Pz());
 111
 112
 113
 114               /* Wt is the weight of the particle*/
 115               STACK.wt[STACK.lstack] = oneone;
 116               STARS.weipri += STACK.wt[STACK.lstack];
 117
 118               // ray
 119               STACK.ilo[STACK.lstack] = 0;
 120               /* From this point .....
 121                * Particle generation (1 for primaries)
 122                */
 123               STACK.lo[STACK.lstack] = 1;
 124
 125               /* User dependent flag:*/
 126               STACK.louse[STACK.lstack] = 0;
 127
 128               /* User dependent spare variables:*/
 129               Int_t ispr = 0;
 130               for (ispr = 0; ispr < mkbmx1; ispr++)
 131                   STACK.sparek[STACK.lstack][ispr] = zerzer;
 132
 133               /* User dependent spare flags:*/
 134               for (ispr = 0; ispr < mkbmx2; ispr++)
 135                   STACK.ispark[STACK.lstack][ispr] = 0;
 136
 137               /* Save the track number of the stack particle:*/
 138               STACK.ispark[STACK.lstack][mkbmx2-1] = itrack;
 139               STACK.nparma++;
 140               STACK.numpar[STACK.lstack] = STACK.nparma;
 141               STACK.nevent[STACK.lstack] = 0;
 142               STACK.dfnear[STACK.lstack] = +zerzer;
 143
 144               /* Particle age (s)*/
 145               STACK.agestk[STACK.lstack] = +zerzer;
 146               STACK.aknshr[STACK.lstack] = -twotwo;
 147
 148               /* Group number for "low" energy neutrons, set to 0 anyway*/
 149               STACK.igroup[STACK.lstack] = 0;
 150
 151               /* Kinetic energy */
 152               STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass();
 153
 154
 155               /* Particle momentum*/
 156               STACK.pmom [STACK.lstack] = particle->P();
 157
 158               /* Cosines (tx,ty,tz)*/
 159               Double_t cosx = particle->Px()/particle->P();
 160               Double_t cosy = particle->Py()/particle->P();
 161               Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy);
 162               if (particle->Pz() < 0.) cosz = -cosz;
 163
 164               STACK.tx [STACK.lstack] = cosx;
 165               STACK.ty [STACK.lstack] = cosy;
 166               STACK.tz [STACK.lstack] = cosz;
 167
 168               /* Polarization cosines:*/
 169               STACK.txpol [STACK.lstack] = -twotwo;
 170               STACK.typol [STACK.lstack] = +zerzer;
 171               STACK.tzpol [STACK.lstack] = +zerzer;
 172
 173               /* Particle coordinates*/
 174               STACK.xa [STACK.lstack] =   particle->Vx();
 175               STACK.ya [STACK.lstack] =   particle->Vy();
 176               STACK.za [STACK.lstack] =   particle->Vz();
 177
 178               printf("Particle Vertex %10.3f %10.3f %10.3f  \n",
 179                      STACK.xa [STACK.lstack],  STACK.ya [STACK.lstack],  STACK.za [STACK.lstack]);
 180
 181
 182
 183               /*  Calculate the total kinetic energy of the primaries: don't change*/
 184               Int_t st_ilo =  STACK.ilo[STACK.lstack];
 185               if ( st_ilo != 0 )
 186                   EPISOR.tkesum +=
 187                       ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6])
 188                        * STACK.wt[STACK.lstack]);
 189               else
 190                   EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]);
 191
 192               /*  Here we ask for the region number of the hitting point.
 193                *     NREG (LSTACK) = ...
 194                *  The following line makes the starting region search much more
 195                *  robust if particles are starting very close to a boundary:
 196                */
 197               geocrs( STACK.tx[STACK.lstack],
 198                       STACK.ty[STACK.lstack],
 199                       STACK.tz[STACK.lstack] );
 200
 201               Int_t idisc;
 202
 203               georeg ( STACK.xa[STACK.lstack],
 204                        STACK.ya[STACK.lstack],
 205                        STACK.za[STACK.lstack],
 206                        STACK.nreg[STACK.lstack],
 207                        idisc);//<-- dummy return variable not used
 208               /*  Do not change these cards:*/
 209               Int_t igeohsm1 = 1;
 210               Int_t igeohsm2 = -11;
 211               geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc );
 212               STACK.nlattc[STACK.lstack] = LTCLCM.mlattc;
 213               soevsv();
 214           }
 215       }
 216   }
 217 }
 218
 219