PrimaryElectronIndex corrected.
[u/mrichter/AliRoot.git] / TFluka / mgdraw.cxx
1 #include <Riostream.h>
2 #include "TVirtualMCApplication.h"
3 #include "TVirtualMCStack.h"
4
5 #include "TFluka.h"
6 #include "TFlukaCodes.h"
7 // Fluka include
8 #include "Fdimpar.h"  //(DIMPAR) fluka include
9 #include "Fdblprc.h"  //(DBLPRC) fluka common
10 #include "Ftrackr.h"  //(TRACKR) fluka common
11 #include "Fopphst.h"  //(OPPHST) fluka common
12 #include "Fflkstk.h"  //(FLKSTK) fluka common
13 #include "Fltclcm.h"  //(LTCLCM) fluka common
14 #include "Fpaprop.h"  //(PAPROP) fluka common
15 #include "Falldlt.h"  //(ALLDLT) fluka common
16
17 #ifndef WIN32
18 # define mgdraw mgdraw_
19 #else
20 # define mgdraw MGDRAW
21 #endif
22
23
24 #include "TGeoManager.h" // <- delete
25
26 extern "C" {
27 void mgdraw(Int_t& icode, Int_t& mreg)
28 {
29     TFluka* fluka =  (TFluka*) gMC;
30     if (mreg == fluka->GetDummyRegion()) return;
31 //
32 //  Make sure that stack has currrent track Id
33 //
34     Int_t trackId = -1;
35     TVirtualMCStack* cppstack = fluka->GetStack();
36     
37     if (TRACKR.jtrack == -1) {
38         trackId = OPPHST.louopp[OPPHST.lstopp];
39         if (trackId == 0) {
40             trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
41         }
42     } else {
43         trackId = TRACKR.ispusr[mkbmx2-1];
44     }
45     
46     Int_t verbosityLevel = fluka->GetVerbosityLevel();
47
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( 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     if (!TRACKR.ispusr[mkbmx2 - 2]) {
87         if (verbosityLevel >= 3) {
88             cout << endl << "mgdraw: energy deposition for:" << trackId
89                  << " icode=" << icode
90                  << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
91                  << " flukaid="<< TRACKR.jtrack
92                  << " mreg=" << mreg
93                  << " np  =" << ALLDLT.nalldl
94                  << endl;
95                  
96         }
97         Int_t nprim = 0;
98 //      printf("mgdraw %5d %5d \n", ALLDLT.nalldl, ALLDLT.lalldl );
99         
100         if ((nprim = ALLDLT.nalldl) > 0) {
101             //
102             // Multiple steps (primary ionisation)
103             if (nprim >= mxalld) {
104                 nprim = mxalld;
105                 Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
106                        ALLDLT.nalldl, 
107                        fluka->PDGFromId(TRACKR.jtrack), 
108                        mreg, 
109                        TRACKR.ptrack, 
110                        TRACKR.ctrack);
111                 
112             }
113             for (Int_t i = 0; i < nprim; i++) {
114                 fluka->SetCurrentPrimaryElectronIndex(i);
115                 (TVirtualMCApplication::Instance())->Stepping();
116                 if (i == 0) fluka->SetTrackIsNew(kFALSE);
117             }
118             fluka->SetCurrentPrimaryElectronIndex(-1);
119         } else {
120             // Single step
121             (TVirtualMCApplication::Instance())->Stepping();
122             fluka->SetTrackIsNew(kFALSE);
123         }
124         
125     } else {
126         //
127         // Tracking is being resumed after secondary tracking
128         //
129         if (verbosityLevel >= 3) {
130             cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
131         }
132
133         fluka->SetTrackIsNew(kTRUE);
134         fluka->SetCaller(kMGResumedTrack);
135         (TVirtualMCApplication::Instance())->Stepping();
136
137         // Reset flag and stored values
138         TRACKR.ispusr[mkbmx2 - 2] = 0;
139         for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
140
141
142         if (verbosityLevel >= 3) {
143             cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
144             cout << " Track= " << trackId << " region = " << mreg << endl;
145         }
146
147         fluka->SetTrackIsNew(kFALSE);
148         fluka->SetCaller(kMGDRAW);
149         (TVirtualMCApplication::Instance())->Stepping();
150     }
151 } // end of mgdraw
152 } // end of extern "C"
153