Primary ionisation electrons come now correctly ordered from FLUKA.

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

diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx
index a192431253b0a9093aa2a1bb51c78a5fbcdc10d7..1c3a90f0fda5a30f7bffaf531116b0048eb88010 100644 (file)
--- a/TFluka/TFluka.cxx
+++ b/TFluka/TFluka.cxx
@@ -255,6 +255,7 @@ void TFluka::Init() {
     // Add ions to PDG Data base
     //
      AddParticlesToPdgDataBase();
+     //
 }
 
 
@@ -265,9 +266,10 @@ void TFluka::FinishGeometry() {
 //
   if (fVerbosityLevel >=3) {
     cout << "==> TFluka::FinishGeometry() called." << endl;
-    printf("----FinishGeometry - nothing to do with TGeo\n");
+    printf("----FinishGeometry - applying misalignment if any\n");
     cout << "<== TFluka::FinishGeometry() called." << endl;
   }  
+  TVirtualMCApplication::Instance()->MisalignGeometry();
 } 
 
 //______________________________________________________________________________ 
@@ -986,7 +988,7 @@ Int_t TFluka::PDGFromId(Int_t id) const
   //
   // Return PDG code and pseudo ENDF code from Fluka code
   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon   
-    Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
+    Int_t idSpecial[6] = {GetIonPdg(2,4), GetIonPdg(2, 3), GetIonPdg(1,3), GetIonPdg(1,2), GetIonPdg(0,0), 50000050};
   // IPTOKP array goes from official to internal
 
     if (id == kFLUKAoptical) {
@@ -998,7 +1000,7 @@ Int_t TFluka::PDGFromId(Int_t id) const
 // Error id    
     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
         if (fVerbosityLevel >= 3)
-            printf("PDGFromId: Error id = 0\n");
+            printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
         return -1;
     }
 // Good id    
@@ -1216,7 +1218,7 @@ void TFluka::InitPhysics()
 // Process Fluka specific scoring options
 //
     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
-    Float_t loginp        = 49.0;
+    Float_t loginp        = -49.0;
     Int_t inp             = 0;
     Int_t nscore          = fUserScore->GetEntries();
     
@@ -1270,7 +1272,7 @@ void TFluka::InitPhysics()
 // Initialisation needed for Cerenkov photon production and transport
     TObjArray *matList = GetFlukaMaterials();
     Int_t nmaterial =  matList->GetEntriesFast();
-    fMaterials = new Int_t[nmaterial+3];
+    fMaterials = new Int_t[nmaterial+25];
     
     for (Int_t im = 0; im < nmaterial; im++)
     {
@@ -1608,7 +1610,8 @@ Int_t TFluka::CorrectFlukaId() const
    // and there is a call to endraw for energy deposition for each of them
    // and they have the track number of their parent, but different identity (pdg)
    // so we want to assign also their parent identity.
-   if( (IsTrackStop() )
+
+   if( (IsTrackStop())
         && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
         && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
       if (fVerbosityLevel >=3)
@@ -1618,7 +1621,11 @@ Int_t TFluka::CorrectFlukaId() const
                << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
       return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
    }
-   return TRACKR.jtrack;
+   if (TRACKR.jtrack <= 64){
+       return TRACKR.jtrack;
+   } else {
+       return TRACKR.j0trck;
+   }
 }
 
 
@@ -1641,6 +1648,7 @@ Double_t TFluka::TrackCharge() const
 // Return charge of the track currently transported
 // PAPROP.ichrge = electric charge of the particle
 // TRACKR.jtrack = identity number of the particle
+    
   FlukaCallerCode_t caller = GetCaller();
   if (caller != kEEDRAW) 
      return PAPROP.ichrge[CorrectFlukaId()+6];
@@ -1936,7 +1944,7 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
         break;
     case kKASOPHrefraction:
         iproc = kPLightRefraction;
-    case kEMSCOlocaledep : 
+    case kEMFSCOlocaldep : 
         iproc = kPPhotoelectric;
         break;
     default:
@@ -2292,8 +2300,6 @@ void TFluka::AddParticlesToPdgDataBase() const
 
     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
 
-    const Int_t kion=10000000;
-
     const Double_t kAu2Gev   = 0.9314943228;
     const Double_t khSlash   = 1.0545726663e-27;
     const Double_t kErg2Gev  = 1/1.6021773349e-3;
@@ -2302,15 +2308,14 @@ void TFluka::AddParticlesToPdgDataBase() const
 //
 // Ions
 //
-
   pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
-                     0,3,"Ion",kion+10020);
+                     0,3,"Ion",GetIonPdg(1,2));
   pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
-                     khShGev/(12.33*kYear2Sec),3,"Ion",kion+10030);
+                     khShGev/(12.33*kYear2Sec),3,"Ion",GetIonPdg(1,3));
   pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
-                     khShGev/(12.33*kYear2Sec),6,"Ion",kion+20040);
+                     khShGev/(12.33*kYear2Sec),6,"Ion",GetIonPdg(2,4));
   pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
-                     0,6,"Ion",kion+20030);
+                     0,6,"Ion",GetIonPdg(2,3));
 }
 
 //
@@ -2330,6 +2335,7 @@ Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
     // Returns kinetic energy of primary electron i
 
     Double_t ekin = -1.;
+    
     if (i >= 0 && i < ALLDLT.nalldl) {
         ekin =  ALLDLT.talldl[i];
     } else {
@@ -2357,4 +2363,33 @@ void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Doubl
        return;
 }
 
+Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
+{
+// Acording to
+// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
 
+  return 1000000000 + 10*1000*z + 10*a + i;
+}  
+     
+void  TFluka::PrimaryIonisationStepping(Int_t nprim)
+{
+// Call Stepping for primary ionisation electrons
+    Int_t i;
+// Protection against nprim > mxalld
+
+// Multiple steps for nprim > 0
+    if (nprim > 0) {
+       for (i = 0; i < nprim; i++) {
+           SetCurrentPrimaryElectronIndex(i);
+           (TVirtualMCApplication::Instance())->Stepping();
+           if (i == 0) SetTrackIsNew(kFALSE);
+       }       
+    } else {
+       // No primary electron ionisation
+       // Call Stepping anyway but flag nprim = 0 as index = -2
+       SetCurrentPrimaryElectronIndex(-2);
+       (TVirtualMCApplication::Instance())->Stepping();
+    }
+    // Reset the index
+    SetCurrentPrimaryElectronIndex(-1);
+}