- Improved support for heavy fragment transport
[u/mrichter/AliRoot.git] / TFluka / mgdraw.cxx
1 #include <Riostream.h>
2 #include "TVirtualMCApplication.h"
3 #include "TVirtualMCStack.h"
4 #include "TFluka.h"
5 #include "TFlukaCodes.h"
6
7 // Fluka includes
8 #include "Fdimpar.h"  //(DIMPAR) 
9 #include "Fdblprc.h"  //(DBLPRC) 
10 #include "Ftrackr.h"  //(TRACKR) 
11 #include "Fopphst.h"  //(OPPHST) 
12 #include "Fflkstk.h"  //(FLKSTK) 
13 #include "Fltclcm.h"  //(LTCLCM) 
14 #include "Fpaprop.h"  //(PAPROP) 
15 #include "Falldlt.h"  //(ALLDLT) 
16 #include "Fdpdxcm.h"  //(DPDXCM) 
17 #include "Fflkmat.h"  //(FLKMAT) 
18
19 #include "TGeoManager.h"
20
21 #ifndef WIN32
22 # define mgdraw mgdraw_
23 #else
24 # define mgdraw MGDRAW
25 #endif
26
27 extern "C" {
28 void mgdraw(Int_t& icode, Int_t& mreg)
29 {
30     TFluka* fluka =  (TFluka*) gMC;
31     if (mreg == fluka->GetDummyRegion()) return;
32 //
33 //  Make sure that stack has currrent track Id
34 //
35     Int_t trackId = -1;
36     TVirtualMCStack* cppstack = fluka->GetStack();
37     
38     if (TRACKR.jtrack == -1) {
39         trackId = OPPHST.louopp[OPPHST.lstopp];
40         if (trackId == 0) {
41             trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
42         }
43     } else {
44         trackId = TRACKR.ispusr[mkbmx2-1];
45     }
46     
47     Int_t verbosityLevel = fluka->GetVerbosityLevel();
48     if (TRACKR.jtrack < -6) {
49        // from -7 to -12 = "heavy" fragment
50        // assing parent id
51        // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
52        if (verbosityLevel >= 3) {
53           cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
54                << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
55        }
56 //       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
57     }
58     
59     cppstack->SetCurrentTrack(trackId);
60 //
61 //    
62     Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
63     fluka->SetMreg(mreg, mlttc);
64     fluka->SetIcode((FlukaProcessCode_t) icode);
65     fluka->SetCaller(kMGDRAW);
66
67     Int_t nodeId;
68     Int_t volId   = fluka->CurrentVolID(nodeId);
69     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
70
71     // check region lattice consistency (debug Ernesto)
72     // *****************************************************
73     if(verbosityLevel>=3 && mreg != volId  && !gGeoManager->IsOutside() ) {
74        cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
75             << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
76             << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
77             << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
78             << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
79             << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
80             << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
81             << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
82         if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
83     }
84     // *****************************************************
85
86     
87     Int_t med   = FLKMAT.medium[mreg - 1];  // Medium
88     Int_t msd   = DPDXCM.msdpdx[med  - 1];  // Iionisation model
89
90     if (!TRACKR.ispusr[mkbmx2 - 2]) {
91         if (verbosityLevel >= 3) {
92             cout << endl << "mgdraw: energy deposition for:" << trackId
93                  << " icode=" << icode
94                  << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
95                  << " flukaid="<< TRACKR.jtrack
96                  << " mreg=" << mreg
97                  << " np  =" << ALLDLT.nalldl
98                  << endl;
99                  
100         }
101
102         if (msd > 0) {
103 //
104 // Primary ionsiations
105             Int_t nprim = ALLDLT.nalldl;
106             if (nprim >= mxalld) {
107                 nprim = mxalld;
108                 Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
109                         ALLDLT.nalldl, 
110                         fluka->PDGFromId(TRACKR.jtrack), 
111                         mreg, 
112                         TRACKR.ptrack, 
113                         TRACKR.ctrack);
114             }
115             fluka->PrimaryIonisationStepping(nprim);
116         } else {
117             // Single step
118             (TVirtualMCApplication::Instance())->Stepping();
119             fluka->SetTrackIsNew(kFALSE);
120         }
121         
122     } else {
123         //
124         // Tracking is being resumed after secondary tracking
125         //
126         if (verbosityLevel >= 3) {
127             cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
128         }
129
130         fluka->SetTrackIsNew(kTRUE);
131         fluka->SetCaller(kMGResumedTrack);
132         (TVirtualMCApplication::Instance())->Stepping();
133
134         // Reset flag and stored values
135         TRACKR.ispusr[mkbmx2 - 2] = 0;
136         for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
137
138
139         if (verbosityLevel >= 3) {
140             cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
141             cout << " Track= " << trackId << " region = " << mreg << endl;
142         }
143
144         fluka->SetTrackIsNew(kFALSE);
145         fluka->SetCaller(kMGDRAW);
146         if (msd == 0) {
147             (TVirtualMCApplication::Instance())->Stepping();
148         } else {
149             Int_t nprim = ALLDLT.nalldl;
150 // Protection against nprim > mxalld
151             if (nprim >= mxalld) {
152                 nprim = mxalld;
153                 Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
154                         ALLDLT.nalldl, 
155                         fluka->PDGFromId(TRACKR.jtrack), 
156                         mreg, 
157                         TRACKR.ptrack, 
158                         TRACKR.ctrack);
159             }
160             fluka->PrimaryIonisationStepping(nprim);
161         } // primary ionisation switched on
162     } // tracking resumed
163 } // end of mgdraw
164 } // end of extern "C"
165