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             for (Int_t i = 0; i < nprim; i++) {
 104                 (TVirtualMCApplication::Instance())->Stepping();
 105                 fluka->SetCurrentPrimaryElectronIndex(i);
 106                 if (i == 0) fluka->SetTrackIsNew(kFALSE);
 107             }
 108             fluka->SetCurrentPrimaryElectronIndex(-1);
 109         } else {
 110             // Single step
 111             (TVirtualMCApplication::Instance())->Stepping();
 112             fluka->SetTrackIsNew(kFALSE);
 113         }
 114
 115     } else {
 116         //
 117         // Tracking is being resumed after secondary tracking
 118         //
 119         if (verbosityLevel >= 3) {
 120             cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
 121         }
 122
 123         fluka->SetTrackIsNew(kTRUE);
 124         fluka->SetCaller(kMGResumedTrack);
 125         (TVirtualMCApplication::Instance())->Stepping();
 126
 127         // Reset flag and stored values
 128         TRACKR.ispusr[mkbmx2 - 2] = 0;
 129         for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
 130
 131
 132         if (verbosityLevel >= 3) {
 133             cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
 134             cout << " Track= " << trackId << " region = " << mreg << endl;
 135         }
 136
 137         fluka->SetTrackIsNew(kFALSE);
 138         fluka->SetCaller(kMGDRAW);
 139         (TVirtualMCApplication::Instance())->Stepping();
 140     }
 141 } // end of mgdraw
 142 } // end of extern "C"
 143