D0 analysis moved to PWG3 (Andrea)

[u/mrichter/AliRoot.git] / TFluka / mgdraw.cxx

diff --git a/TFluka/mgdraw.cxx b/TFluka/mgdraw.cxx
index 5e89412aae5e29384699e26898531f56fe90ac58..33accf8bf6a0155efe12321d9589a891eb092ba9 100644 (file)
--- a/TFluka/mgdraw.cxx
+++ b/TFluka/mgdraw.cxx
@@ -1,17 +1,22 @@
 #include <Riostream.h>
 #include "TVirtualMCApplication.h"
 #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
+
+// Fluka includes
+#include "Fdimpar.h"  //(DIMPAR) 
+#include "Fdblprc.h"  //(DBLPRC) 
+#include "Ftrackr.h"  //(TRACKR) 
+#include "Fopphst.h"  //(OPPHST) 
+#include "Fflkstk.h"  //(FLKSTK) 
+#include "Fltclcm.h"  //(LTCLCM) 
+#include "Fpaprop.h"  //(PAPROP) 
+#include "Falldlt.h"  //(ALLDLT) 
+#include "Fdpdxcm.h"  //(DPDXCM) 
+#include "Fflkmat.h"  //(FLKMAT) 
+
+#include "TGeoManager.h"
 
 #ifndef WIN32
 # define mgdraw mgdraw_
@@ -42,11 +47,11 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     Int_t verbosityLevel = fluka->GetVerbosityLevel();
 
     if (TRACKR.jtrack < -6) {
-       // (Unknow heavy Ion???)
-       // assing parent id ???
+       // from -7 to -12 = "heavy" fragment
+       // assing parent id
        // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
        if (verbosityLevel >= 3) {
-          cout << "mgdraw: jtrack < -6 =" << TRACKR.jtrack
+          cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
                << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
        }
        TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
@@ -55,22 +60,82 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     cppstack->SetCurrentTrack(trackId);
 //
 //    
-    Int_t mlttc = LTCLCM.mlatm1;
+    Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
     fluka->SetMreg(mreg, mlttc);
-    fluka->SetNewreg(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(verbosityLevel>=3 && 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;
+    }
+    // *****************************************************
+
+    
+    Int_t med   = FLKMAT.medium[mreg - 1];  // Medium
+    Int_t msd   = DPDXCM.msdpdx[med  - 1];  // Iionisation model
+
     if (!TRACKR.ispusr[mkbmx2 - 2]) {
-        //
-        // Single step
-        if (verbosityLevel >= 3) {
-           cout << endl << "mgdraw: energy deposition for:" << trackId
-                 << " icode=" << icode
-                 << " pdg=" << fluka->PDGFromId(TRACKR.jtrack) << " flukaid="<< TRACKR.jtrack << endl;
-        }
-        (TVirtualMCApplication::Instance())->Stepping();
-        fluka->SetTrackIsNew(kFALSE);
+       if (verbosityLevel >= 3) {
+           cout << endl << "mgdraw: energy deposition for:" << trackId
+                << " icode=" << icode
+                << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+                << " flukaid="<< TRACKR.jtrack
+                << " mreg=" << mreg
+                << " np  =" << ALLDLT.nalldl
+                << endl;
+                
+       }
+
+       if (msd > 0) {
+//
+// Primary ionsiations
+           Int_t nprim = ALLDLT.nalldl;
+// Protection against nprim > mxalld
+           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);
+               
+           }
+// Multiple steps for nprim > 0
+           if (nprim > 0) {
+               for (Int_t i = 0; i < nprim; i++) {
+                   fluka->SetCurrentPrimaryElectronIndex(i);
+                   (TVirtualMCApplication::Instance())->Stepping();
+                   if (i == 0) fluka->SetTrackIsNew(kFALSE);
+               }
+           } else {
+               // No primary electron ionisation
+               // Call Stepping anyway but flag nprim = 0 as index = -2
+               fluka->SetCurrentPrimaryElectronIndex(-2);
+               (TVirtualMCApplication::Instance())->Stepping();
+           }
+           // Reset the index
+           fluka->SetCurrentPrimaryElectronIndex(-1);
+       } else {
+           // Single step
+           (TVirtualMCApplication::Instance())->Stepping();
+           fluka->SetTrackIsNew(kFALSE);
+       }
+       
     } else {
         //
         // Tracking is being resumed after secondary tracking
@@ -95,7 +160,9 @@ void mgdraw(Int_t& icode, Int_t& mreg)
 
         fluka->SetTrackIsNew(kFALSE);
         fluka->SetCaller(kMGDRAW);
+       if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-2);
         (TVirtualMCApplication::Instance())->Stepping();
+       if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-1);
     }
 } // end of mgdraw
 } // end of extern "C"