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
  49     if (TRACKR.jtrack < -6) {
  50        // from -7 to -12 = "heavy" fragment
  51        // assing parent id
  52        // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
  53        if (verbosityLevel >= 3) {
  54           cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
  55                << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
  56        }
  57        TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
  58     }
  59
  60     cppstack->SetCurrentTrack(trackId);
  61 //
  62 //
  63     Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
  64     fluka->SetMreg(mreg, mlttc);
  65     fluka->SetIcode((FlukaProcessCode_t) icode);
  66     fluka->SetCaller(kMGDRAW);
  67
  68     Int_t nodeId;
  69     Int_t volId   = fluka->CurrentVolID(nodeId);
  70     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
  71
  72     // check region lattice consistency (debug Ernesto)
  73     // *****************************************************
  74     if( mreg != volId  && !gGeoManager->IsOutside() ) {
  75        cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
  76             << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
  77             << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
  78             << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
  79             << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
  80             << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
  81             << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
  82             << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
  83         if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
  84     }
  85     // *****************************************************
  86
  87
  88     Int_t med   = FLKMAT.medium[mreg - 1];  // Medium
  89     Int_t msd   = DPDXCM.msdpdx[med  - 1];  // Iionisation model
  90
  91     if (!TRACKR.ispusr[mkbmx2 - 2]) {
  92         if (verbosityLevel >= 3) {
  93             cout << endl << "mgdraw: energy deposition for:" << trackId
  94                  << " icode=" << icode
  95                  << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
  96                  << " flukaid="<< TRACKR.jtrack
  97                  << " mreg=" << mreg
  98                  << " np  =" << ALLDLT.nalldl
  99                  << endl;
 100
 101         }
 102
 103         if (msd > 0) {
 104 //
 105 // Primary ionsiations
 106             Int_t nprim = ALLDLT.nalldl;
 107 // Protection against nprim > mxalld
 108             if (nprim >= mxalld) {
 109                 nprim = mxalld;
 110                 Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n",
 111                        ALLDLT.nalldl,
 112                        fluka->PDGFromId(TRACKR.jtrack),
 113                        mreg,
 114                        TRACKR.ptrack,
 115                        TRACKR.ctrack);
 116
 117             }
 118 // Multiple steps for nprim > 0
 119             if (nprim > 0) {
 120                 for (Int_t i = 0; i < nprim; i++) {
 121                     fluka->SetCurrentPrimaryElectronIndex(i);
 122                     (TVirtualMCApplication::Instance())->Stepping();
 123                     if (i == 0) fluka->SetTrackIsNew(kFALSE);
 124                 }
 125             } else {
 126                 // No primary electron ionisation
 127                 // Call Stepping anyway but flag nprim = 0 as index = -2
 128                 fluka->SetCurrentPrimaryElectronIndex(-2);
 129                 (TVirtualMCApplication::Instance())->Stepping();
 130             }
 131             // Reset the index
 132             fluka->SetCurrentPrimaryElectronIndex(-1);
 133         } else {
 134             // Single step
 135             (TVirtualMCApplication::Instance())->Stepping();
 136             fluka->SetTrackIsNew(kFALSE);
 137         }
 138
 139     } else {
 140         //
 141         // Tracking is being resumed after secondary tracking
 142         //
 143         if (verbosityLevel >= 3) {
 144             cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
 145         }
 146
 147         fluka->SetTrackIsNew(kTRUE);
 148         fluka->SetCaller(kMGResumedTrack);
 149         (TVirtualMCApplication::Instance())->Stepping();
 150
 151         // Reset flag and stored values
 152         TRACKR.ispusr[mkbmx2 - 2] = 0;
 153         for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
 154
 155
 156         if (verbosityLevel >= 3) {
 157             cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
 158             cout << " Track= " << trackId << " region = " << mreg << endl;
 159         }
 160
 161         fluka->SetTrackIsNew(kFALSE);
 162         fluka->SetCaller(kMGDRAW);
 163         if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-2);
 164         (TVirtualMCApplication::Instance())->Stepping();
 165         if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-1);
 166     }
 167 } // end of mgdraw
 168 } // end of extern "C"
 169