]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TFluka/mgdraw.cxx
Common blocks to access material properties and
[u/mrichter/AliRoot.git] / TFluka / mgdraw.cxx
index f2ef1cb1e44ae0e4ad1f1ca446b8537916ba3106..fa6540b570c2b6edfd96fc22c9fb0241309c1994 100644 (file)
@@ -3,13 +3,16 @@
 #include "TVirtualMCStack.h"
 
 #include "TFluka.h"
-
+#include "TFlukaCodes.h"
 // Fluka include
 #include "Fdimpar.h"  //(DIMPAR) fluka include
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Ftrackr.h"  //(TRACKR) fluka common
 #include "Fopphst.h"  //(OPPHST) fluka common
 #include "Fflkstk.h"  //(FLKSTK) fluka common
+#include "Fltclcm.h"  //(LTCLCM) fluka common
+#include "Fpaprop.h"  //(PAPROP) fluka common
+#include "Falldlt.h"  //(ALLDLT) fluka common
 
 #ifndef WIN32
 # define mgdraw mgdraw_
 # define mgdraw MGDRAW
 #endif
 
+
+#include "TGeoManager.h" // <- delete
+
 extern "C" {
 void mgdraw(Int_t& icode, Int_t& mreg)
 {
     TFluka* fluka =  (TFluka*) gMC;
-    Int_t verbosityLevel = fluka->GetVerbosityLevel();
+    if (mreg == fluka->GetDummyRegion()) return;
 //
 //  Make sure that stack has currrent track Id
 //
@@ -29,64 +35,119 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     TVirtualMCStack* cppstack = fluka->GetStack();
     
     if (TRACKR.jtrack == -1) {
-       trackId = OPPHST.louopp[OPPHST.lstopp];
-       if (trackId == 0) {
-           trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
-       }
+        trackId = OPPHST.louopp[OPPHST.lstopp];
+        if (trackId == 0) {
+            trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
+        }
     } else {
-       trackId = TRACKR.ispusr[mkbmx2-1];
+        trackId = TRACKR.ispusr[mkbmx2-1];
+    }
+    
+    Int_t verbosityLevel = fluka->GetVerbosityLevel();
+
+    if (TRACKR.jtrack < -6) {
+       // from -7 to -12 = "heavy" fragment
+       // assing parent id
+       // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
+       if (verbosityLevel >= 3) {
+          cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
+               << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
+       }
+       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
     }
     
     cppstack->SetCurrentTrack(trackId);
 //
 //    
-    fluka->SetMreg(mreg);
-    fluka->SetNewreg(mreg);
-    fluka->SetIcode(icode);
-    fluka->SetCaller(4);
-
-    if (!TRACKR.ispusr[mkbmx2 - 2]) {
-       //
-       // Single step
-       if (TRACKR.jtrack == -1 && trackId == 109340) {
-           cout << endl << " !!! I am in mgdraw - calling Stepping(): " << icode << endl;
-           cout << endl << " Track Id = " << trackId << " region = " << mreg << endl;
-           printf("Stepsize %13.5e \n", fluka->TrackStep());
-       }
+    Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
+    fluka->SetMreg(mreg, mlttc);
+    fluka->SetIcode((FlukaProcessCode_t) icode);
+    fluka->SetCaller(kMGDRAW);
 
+    Int_t nodeId;
+    Int_t volId   = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
 
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+       cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
+    }
+    // *****************************************************
 
-      
-       (TVirtualMCApplication::Instance())->Stepping();
-       fluka->SetTrackIsNew(kFALSE);
-    } else {
-       //
-       // Tracking is being resumed after secondary tracking
-       //
+    if (!TRACKR.ispusr[mkbmx2 - 2]) {
        if (verbosityLevel >= 3) {
-           cout << endl << " !!! I am in mgdraw - resuming Stepping(): " << trackId << endl;
+           cout << endl << "mgdraw: energy deposition for:" << trackId
+                << " icode=" << icode
+                << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+                << " flukaid="<< TRACKR.jtrack
+                << " mreg=" << mreg
+                << " np  =" << ALLDLT.nalldl
+                << endl;
+                
        }
+       Int_t nprim = 0;
+//     printf("mgdraw %5d %5d \n", ALLDLT.nalldl, ALLDLT.lalldl );
        
-       fluka->SetTrackIsNew(kTRUE);
-       fluka->SetCaller(40);
-       (TVirtualMCApplication::Instance())->Stepping();
+       if ((nprim = ALLDLT.nalldl) > 0) {
+           //
+           // Multiple steps (primary ionisation)
+           if (nprim >= mxalld) {
+               nprim = mxalld;
+               Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
+                      ALLDLT.nalldl, 
+                      fluka->PDGFromId(TRACKR.jtrack), 
+                      mreg, 
+                      TRACKR.ptrack, 
+                      TRACKR.ctrack);
+               
+           }
+           for (Int_t i = 0; i < nprim; i++) {
+               fluka->SetCurrentPrimaryElectronIndex(i);
+               (TVirtualMCApplication::Instance())->Stepping();
+               if (i == 0) fluka->SetTrackIsNew(kFALSE);
+           }
+           fluka->SetCurrentPrimaryElectronIndex(-1);
+       } else {
+           // Single step
+           (TVirtualMCApplication::Instance())->Stepping();
+           fluka->SetTrackIsNew(kFALSE);
+       }
+       
+    } else {
+        //
+        // Tracking is being resumed after secondary tracking
+        //
+        if (verbosityLevel >= 3) {
+            cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
+        }
 
-       // Reset flag and stored values
-       TRACKR.ispusr[mkbmx2 - 2] = 0;
-       for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
+        fluka->SetTrackIsNew(kTRUE);
+        fluka->SetCaller(kMGResumedTrack);
+        (TVirtualMCApplication::Instance())->Stepping();
 
+        // Reset flag and stored values
+        TRACKR.ispusr[mkbmx2 - 2] = 0;
+        for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
 
-       if (verbosityLevel >= 3) {
-           cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
-           cout << endl << " Track Id = " << trackId << " region = " << mreg << endl;
-       }
 
-       fluka->SetTrackIsNew(kFALSE);
-       fluka->SetCaller(4);
-       (TVirtualMCApplication::Instance())->Stepping();
+        if (verbosityLevel >= 3) {
+            cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
+            cout << " Track= " << trackId << " region = " << mreg << endl;
+        }
+
+        fluka->SetTrackIsNew(kFALSE);
+        fluka->SetCaller(kMGDRAW);
+        (TVirtualMCApplication::Instance())->Stepping();
     }
-    
-    
 } // end of mgdraw
 } // end of extern "C"