c07ec576313a9365024e8efb1fd77d1af53f068b
[u/mrichter/AliRoot.git] / TFluka / source.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 "Fsourcm.h"  //(EPISOR) fluka common
8 #include "Fflkstk.h"  //(FLKSTK) fluka common
9 #include "Fsumcou.h"  //(SUMCOU) fluka common
10 #include "Fpaprop.h"  //(PAPROP) fluka common
11 #include "Fltclcm.h"  //(LTCLCM) fluka common
12 #include "Fopphst.h"  //(OPPHST) fluka common
13 #include "Fioiocm.h"  //(IOIOCM) fluka common
14 #include "Fbeamcm.h"  //(BEAMCM) fluka common
15 //Virutal MC
16 #include "TFluka.h"
17 #include "TFlukaIon.h"
18
19 #include "TVirtualMCStack.h"
20 //#include "TVirtualMCApplication.h"
21
22 #include "TParticle.h"
23 #include "TVector3.h"
24
25 //Other
26 #include <Riostream.h>
27
28 #ifndef WIN32
29 # define source source_
30 # define geocrs geocrs_
31 # define georeg georeg_
32 # define geohsm geohsm_
33 # define soevsv soevsv_
34 # define dcdion dcdion_
35 # define setion setion_
36 # define stisbm stisbm_
37 #else
38 # define source SOURCE
39 # define geocrs GEOCRS
40 # define georeg GEOREG
41 # define geohsm GEOHSM
42 # define soevsv SOEVSV
43 # define dcdion DCDION
44 # define setion SETION
45 # define stisbm STISBM
46 #endif
47 extern "C" {
48   //
49   // Prototypes for FLUKA functions
50   //
51   void type_of_call geocrs(Double_t &, Double_t &, Double_t &);
52   void type_of_call georeg(Double_t &, Double_t &, Double_t &, 
53                            Int_t &, Int_t &);
54   void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
55   void type_of_call soevsv();
56   void type_of_call dcdion(Int_t &);
57   void type_of_call setion(Int_t &);
58   void type_of_call stisbm(Int_t &, Int_t &, Int_t &);
59
60  /*
61    *----------------------------------------------------------------------*
62    *                                                                      *
63    *     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
64    *                                                   Infn - Milan       *
65    *                                                                      *
66    *     Last change on 21-jun-98     by    Alfredo Ferrari               *
67    *                                                                      *
68    *     C++ version on 27-sep-02     by    Isidro Gonzalez               *
69    *                                                                      *
70    *  This is just an example of a possible user written source routine.  *
71    *  note that the beam card still has some meaning - in the scoring the *
72    *  maximum momentum used in deciding the binning is taken from the     *
73    *  beam momentum.  Other beam card parameters are obsolete.            *
74    *                                                                      *
75    *----------------------------------------------------------------------*/
76
77   void source(Int_t& nomore) {
78 // Get the pointer to TFluka
79     TFluka* fluka = (TFluka*)gMC;
80
81     Int_t verbosityLevel = fluka->GetVerbosityLevel();
82     Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
83     if (debug) {
84       cout << "==> source(" << nomore << ")" << endl;
85       cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
86     }  
87
88     static Bool_t lfirst = true;
89     static Bool_t particleIsPrimary = true;
90     static Bool_t lastParticleWasPrimary = true;
91
92     nomore = 0;
93
94     
95 //  Get the stack 
96     TVirtualMCStack* cppstack = fluka->GetStack();
97
98     TParticle* particle;
99     Int_t itrack = -1;
100     Int_t  nprim  = cppstack->GetNprimary();
101 //  Get the next particle from the stack
102     particle  = cppstack->PopNextTrack(itrack);
103     fluka->SetTrackIsNew(kTRUE);
104     if (itrack == (nprim - 1)) lfirst = true;
105 //  Is this a secondary not handled by Fluka, i.e. a particle added by user action ?
106     lastParticleWasPrimary = particleIsPrimary;
107     
108     if (itrack >= nprim) {
109         particleIsPrimary = kFALSE;
110     } else {
111         particleIsPrimary = kTRUE;
112     }
113
114     if (lfirst) {
115         SOURCM.tkesum = zerzer;
116         lfirst = false;
117         SOURCM.lussrc = true;
118     } else {
119 //
120 // Post-track actions for primary track
121 //
122         if (particleIsPrimary) {
123             TVirtualMCApplication::Instance()->PostTrack();
124             TVirtualMCApplication::Instance()->FinishPrimary();
125             if ((itrack%10)==0)
126                 cout << "=== TRACKING PRIMARY "<< itrack <<" ===" << endl;
127             //printf("=== TRACKING PRIMARY %d ===\n", itrack);
128         }
129     }
130
131     // Exit if itrack is negative (-1). Set lsouit to false to mark last track for this event
132
133     if (itrack<0) {
134       nomore = 1;
135       SOURCM.lsouit = false;
136       if (debug) {
137          cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
138          cout << "\t* No more particles. Exiting..." << endl;
139          cout << "<== source(" << nomore << ")" << endl;
140       }   
141       return;
142     }
143     
144     //
145     // Handle user event abortion
146     if (fluka->EventIsStopped()) {
147         printf("Event has been stopped by user !");
148         fluka->SetStopEvent(kFALSE);
149         nomore = 1;
150         SOURCM.lsouit = false;
151         return;
152     }
153
154     //Get some info about the particle and print it
155     //
156     //pdg code
157     Int_t pdg = particle->GetPdgCode();
158     TVector3 polarisation;
159     particle->GetPolarisation(polarisation);
160     if (debug) {
161        cout << "\t* Particle " << itrack << " retrieved..." << endl;
162        cout << "\t\t+ Name = " << particle->GetName() << endl;
163        cout << "\t\t+ PDG/Fluka code = " << pdg 
164             << " / " << fluka->IdFromPDG(pdg) << endl;
165        cout << "\t\t+ P = (" 
166             << particle->Px() << " , "
167             << particle->Py() << " , "
168             << particle->Pz() << " ) --> "
169             << particle->P() << " GeV "
170             << particle->Energy() << " GeV "
171             << particle->GetMass() << " GeV " << endl;
172     }   
173     /* Npflka is the stack counter: of course any time source is called it
174      * must be =0
175      */
176     /* Cosines (tx,ty,tz)*/
177     Double_t cosx = particle->Px()/particle->P();
178     Double_t cosy = particle->Py()/particle->P();
179     Double_t cosxy = cosx * cosx + cosy * cosy;
180     Double_t cosz;
181
182     if (cosxy < 1.) {
183         cosz = TMath::Sqrt(oneone - cosxy);
184     } else {
185         cosx /= TMath::Sqrt(cosxy);
186         cosy /= TMath::Sqrt(cosxy);     
187         cosz = 0.;
188     }
189     
190     
191     
192     if (particle->Pz() < 0.) cosz = -cosz;    
193
194     if (pdg != 50000050 &&  pdg !=  50000051) {
195         FLKSTK.npflka++;
196         Int_t ifl =  fluka-> IdFromPDG(pdg);
197         Int_t ionid;
198         
199         if (ifl == -2) {
200             Int_t ia = TFlukaIon::GetA(pdg);
201             Int_t iz = TFlukaIon::GetZ(pdg);
202             Int_t is = TFlukaIon::GetIsomerNumber(pdg);
203             printf("Isomer %5d \n", is);
204             
205             if (is == 0) {
206                 IOIOCM.iproa = ia;
207                 IOIOCM.iproz = iz;
208                 BEAMCM.ijhion = iz * 1000 + ia;
209                 BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
210                 ionid = BEAMCM.ijhion;
211                 dcdion(ionid);
212                 setion(ionid);
213                 FLKSTK.iloflk[FLKSTK.npflka] = ionid;
214                 FLKSTK.lraddc[FLKSTK.npflka] = 0;
215             } else {
216                 BEAMCM.lrdbea = 1;
217                 stisbm(ia, iz, is);
218                 BEAMCM.ijhion = iz * 1000 + ia;
219                 BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
220                 ionid = BEAMCM.ijhion;
221                 dcdion(ionid);
222                 setion(ionid);
223                 FLKSTK.iloflk[FLKSTK.npflka] = ionid;
224                 FLKSTK.lraddc[FLKSTK.npflka] = 1;
225             }
226         } else {
227             FLKSTK.iloflk[FLKSTK.npflka] = ifl;
228             FLKSTK.lraddc[FLKSTK.npflka] = 0;
229             ionid = ifl;
230         }
231         
232         /* Wtflk is the weight of the particle*/
233         FLKSTK.wtflk[FLKSTK.npflka] = oneone;
234         SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
235         
236         FLKSTK.loflk[FLKSTK.npflka] = 1;
237         
238         /* User dependent flag:*/
239         FLKSTK.louse[FLKSTK.npflka] = 0;
240
241         /* User dependent spare variables:*/
242         Int_t ispr = 0;
243         for (ispr = 0; ispr < mkbmx1; ispr++)
244             FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
245
246         /* User dependent spare flags:*/
247         for (ispr = 0; ispr < mkbmx2; ispr++)
248             FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
249
250         /* Save the track number of the stack particle:*/
251         FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
252         FLKSTK.nparma++;
253         FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
254         FLKSTK.nevent[FLKSTK.npflka] = 0;
255         FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
256         
257         /* Particle age (s)*/
258         FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
259         FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
260         FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
261
262         /* Group number for "low" energy neutrons, set to 0 anyway*/
263         FLKSTK.igroup[FLKSTK.npflka] = 0;
264         
265         /* Kinetic energy */
266         Double_t p    = particle->P();
267 //        Double_t mass = PAPROP.am[ifl + 6];
268         Double_t mass = PAPROP.am[ionid + 6];
269         FLKSTK.tkeflk[FLKSTK.npflka] =  TMath::Sqrt( p * p + mass * mass) - mass;
270         /* Particle momentum*/
271         FLKSTK.pmoflk [FLKSTK.npflka] = p;
272
273         FLKSTK.txflk [FLKSTK.npflka] = cosx;
274         FLKSTK.tyflk [FLKSTK.npflka] = cosy;
275         FLKSTK.tzflk [FLKSTK.npflka] = cosz;
276     
277         /* Polarization cosines:*/
278         if (polarisation.Mag()) {
279             Double_t cospolx = polarisation.Px() / polarisation.Mag();
280             Double_t cospoly = polarisation.Py() / polarisation.Mag();
281             Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
282             FLKSTK.txpol [FLKSTK.npflka] = cospolx;
283             FLKSTK.typol [FLKSTK.npflka] = cospoly;
284             FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
285         }
286         else {
287             FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
288             FLKSTK.typol [FLKSTK.npflka] = +zerzer;
289             FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
290         }
291
292         /* Particle coordinates*/
293         // Vertext coordinates;
294         FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
295         FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
296         FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
297
298         /*  Calculate the total kinetic energy of the primaries: don't change*/
299         Int_t st_ilo =  FLKSTK.iloflk[FLKSTK.npflka];
300         if ( st_ilo != 0 )
301             SOURCM.tkesum +=
302                 ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
303                  * FLKSTK.wtflk[FLKSTK.npflka]);
304         else
305             SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
306
307         /*  Here we ask for the region number of the hitting point.
308          *     NRGFLK (LFLKSTK) = ...
309          *  The following line makes the starting region search much more
310          *  robust if particles are starting very close to a boundary:
311          */
312         geocrs( FLKSTK.txflk[FLKSTK.npflka],
313                 FLKSTK.tyflk[FLKSTK.npflka],
314                 FLKSTK.tzflk[FLKSTK.npflka] );
315     
316         Int_t idisc;
317
318         georeg ( FLKSTK.xflk[FLKSTK.npflka],
319                  FLKSTK.yflk[FLKSTK.npflka],
320                  FLKSTK.zflk[FLKSTK.npflka],
321                  FLKSTK.nrgflk[FLKSTK.npflka],
322                  idisc);//<-- dummy return variable not used
323         /*  Do not change these cards:*/
324         Int_t igeohsm1 = 1;
325         Int_t igeohsm2 = -11;
326         geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
327         FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
328         soevsv();
329     } else {
330         //
331         // Next particle is optical photon
332         //
333         OPPHST.lstopp++;
334         OPPHST.donear [OPPHST.lstopp - 1] = 0.;
335
336         OPPHST.xoptph [OPPHST.lstopp - 1] = particle->Vx();
337         OPPHST.yoptph [OPPHST.lstopp - 1] = particle->Vy();
338         OPPHST.zoptph [OPPHST.lstopp - 1] = particle->Vz();
339
340         OPPHST.txopph [OPPHST.lstopp - 1] = cosx;
341         OPPHST.tyopph [OPPHST.lstopp - 1] = cosy;
342         OPPHST.tzopph [OPPHST.lstopp - 1] = cosz;
343
344
345         if (polarisation.Mag()) {
346             Double_t cospolx = polarisation.Px() / polarisation.Mag();
347             Double_t cospoly = polarisation.Py() / polarisation.Mag();
348             Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
349             OPPHST.txpopp [OPPHST.lstopp - 1] = cospolx;
350             OPPHST.typopp [OPPHST.lstopp - 1] = cospoly;
351             OPPHST.tzpopp [OPPHST.lstopp - 1] = cospolz;
352         }
353         else {
354             OPPHST.txpopp [OPPHST.lstopp - 1] = -twotwo;
355             OPPHST.typopp [OPPHST.lstopp - 1] = +zerzer;
356             OPPHST.tzpopp [OPPHST.lstopp - 1] = +zerzer;
357         }
358
359         geocrs( OPPHST.txopph[OPPHST.lstopp - 1],
360                 OPPHST.tyopph[OPPHST.lstopp - 1],
361                 OPPHST.tzopph[OPPHST.lstopp - 1] );
362
363         Int_t idisc;
364
365         georeg ( OPPHST.xoptph[OPPHST.lstopp - 1],
366                  OPPHST.yoptph[OPPHST.lstopp - 1],
367                  OPPHST.zoptph[OPPHST.lstopp - 1],
368                  OPPHST.nregop[OPPHST.lstopp - 1],
369                  idisc);//<-- dummy return variable not used
370
371         OPPHST.wtopph [OPPHST.lstopp - 1] = particle->GetWeight();
372         OPPHST.poptph [OPPHST.lstopp - 1] = particle->P();
373         OPPHST.agopph [OPPHST.lstopp - 1] = particle->T();
374         OPPHST.cmpopp [OPPHST.lstopp - 1] = +zerzer;
375         OPPHST.loopph [OPPHST.lstopp - 1] = 0;
376         OPPHST.louopp [OPPHST.lstopp - 1] = itrack;
377         OPPHST.nlatop [OPPHST.lstopp - 1] = LTCLCM.mlattc;
378      }
379
380 //
381 //  Pre-track actions at for primary tracks
382 //
383     if (particleIsPrimary) {
384         fluka->SetCaller(kSODRAW);
385         TVirtualMCApplication::Instance()->BeginPrimary();
386         TVirtualMCApplication::Instance()->PreTrack();
387     }
388 //
389     if (debug) cout << "<== source(" << nomore << ")" << endl;
390   }
391 }