refactor: extract MeanMaterialBudget to standalone function

CHANGELOG.md

@@ -208,6 +208,7 @@ it in future.
   - Update Python code in shipDigiReco.py to use modern constructors
   - Bump ClassDef versions to 2 for schema evolution
   - Add TrackInfo to RNTuple I/O test suite
+  - Move MeanMaterialBudget to standalone function to reduce duplication and allow proper attribution of code to ALICE
 
 ### Removed
 

LICENSES/BSD-3-Clause.txt

@@ -0,0 +1,11 @@
+Copyright (c) <year> <owner>.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

REUSE.toml

@@ -3,6 +3,14 @@ SPDX-PackageName = "FairShip"
 SPDX-PackageSupplier = "CERN for the benefit of the SHiP Collaboration"
 SPDX-PackageDownloadLocation = "https://github.com/ShipSoft/FairShip"
 
+# MeanMaterialBudget function from ALICE AliRoot (BSD-3-Clause)
+[[annotations]]
+path = ["shipgen/MeanMaterialBudget.h", "shipgen/MeanMaterialBudget.cxx"]
+precedence = "aggregate"
+SPDX-FileCopyrightText = "ALICE Experiment at CERN, All rights reserved"
+SPDX-License-Identifier = "BSD-3-Clause"
+
+# Default FairShip licence for all other files
 [[annotations]]
 path = "**"
 precedence = "aggregate"

shipgen/CMakeLists.txt

@@ -76,7 +76,8 @@ set(SRCS
     FixedTargetGenerator.cxx
     EvtCalcGenerator.cxx
     TEvtGenDecayer.cxx
-    BeamSmearingUtils.cxx)
+    BeamSmearingUtils.cxx
+    MeanMaterialBudget.cxx)
 
 set(LINKDEF LinkDef.h)
 set(LIBRARY_NAME ShipGen)

shipgen/FixedTargetGenerator.cxx

@@ -9,6 +9,7 @@
 #include "Pythia8Plugins/EvtGen.h"
 #include "FixedTargetGenerator.h"
 #include "HNLPythia8Generator.h"
+#include "MeanMaterialBudget.h"
 #include "ShipUnit.h"
 #include "TGeoBBox.h"
 #include "TGeoNode.h"
@@ -265,7 +266,7 @@ Bool_t FixedTargetGenerator::Init()
    end[2]=endZ;
    LOG(info) << "FixedTargetGenerator: Using geometry-based target coordinates startZ=" << startZ << " endZ=" << endZ;
    //find maximum interaction length
-   bparam = fMaterialInvestigator->MeanMaterialBudget(start, end, mparam);
+   bparam = shipgen::MeanMaterialBudget(start, end, mparam);
    maxCrossSection =  mparam[9];
   } else if (targetName!=""){
    // Fallback to fragile TGeo navigation for backward compatibility
@@ -301,7 +302,7 @@ Bool_t FixedTargetGenerator::Init()
    end[1]=yOff;
    end[2]=endZ;
 //find maximum interaction length
-   bparam = fMaterialInvestigator->MeanMaterialBudget(start, end, mparam);
+   bparam = shipgen::MeanMaterialBudget(start, end, mparam);
    maxCrossSection =  mparam[9];
   } else {
       LOG(fatal) << "No target set.";
@@ -344,7 +345,7 @@ Bool_t FixedTargetGenerator::ReadEvent(FairPrimaryGenerator* cpg)
  //place x,y,z uniform along path
       zinter = gRandom->Uniform(zinterStart,end[2]);
       Double_t point[3]={xOff,yOff,zinter};
-      bparam = fMaterialInvestigator->MeanMaterialBudget(start, point, mparam);
+      bparam = shipgen::MeanMaterialBudget(start, point, mparam);
       Double_t interLength = mparam[8];
       TGeoNode *node = gGeoManager->FindNode(point[0],point[1],point[2]);
       TGeoMaterial *mat = 0;

shipgen/GenieGenerator.cxx

@@ -9,6 +9,7 @@
 #include "TRandom.h"
 #include "FairPrimaryGenerator.h"
 #include "GenieGenerator.h"
+#include "MeanMaterialBudget.h"
 #include "TGeoVolume.h"
 #include "TGeoNode.h"
 #include "TGeoManager.h"
@@ -67,172 +68,6 @@ Bool_t GenieGenerator::Init(const char* fileName, const int firstEvent) {
   return kTRUE;
 }
 // -------------------------------------------------------------------------
-Double_t GenieGenerator::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
-{
-  //
-  // Calculate mean material budget and material properties between
-  //    the points "start" and "end".
-  //
-  // "mparam" - parameters used for the energy and multiple scattering
-  //  corrections:
-  //
-  // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
-  // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
-  // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
-  // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
-  // mparam[4] - length: sum(x_i) [cm]
-  // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
-  // mparam[6] - number of boundary crosses
-  // mparam[7] - maximum density encountered (g/cm^3)
-  // mparam[8] - equivalent interaction length fraction: sum(x_i/I0_i) [adimensional]
-  // mparam[9] - maximum cross section encountered (mbarn)
-  //
-  //  Origin:  Marian Ivanov, [email protected]
-  //
-  //  Corrections and improvements by
-  //        Andrea Dainese, [email protected],
-  //        Andrei Gheata,  [email protected]
-  //        Thomas Ruf,  [email protected]
-  //
-
-  mparam[0]=0; mparam[1]=1; mparam[2]=0; mparam[3]=0; mparam[4]=0;
-  mparam[5]=0; mparam[6]=0; mparam[7]=0; mparam[8]=0; mparam[9]=0;
-  //
-  Double_t bparam[7]; // total parameters
-  Double_t lparam[7]; // local parameters
-  Double_t mbarn = 1E-3*1E-24*TMath::Na(); // cm^2 * Avogadro
-
-  for (Int_t i=0;i<7;i++) bparam[i]=0;
-
-  if (!gGeoManager) {
-    //AliFatalClass("No TGeo\n");
-    return 0.;
-  }
-  //
-  Double_t length;
-  Double_t dir[3];
-  length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
-                       (end[1]-start[1])*(end[1]-start[1])+
-                       (end[2]-start[2])*(end[2]-start[2]));
-  mparam[4]=length;
-  if (length<TGeoShape::Tolerance()) return 0.0;
-  Double_t invlen = 1./length;
-  dir[0] = (end[0]-start[0])*invlen;
-  dir[1] = (end[1]-start[1])*invlen;
-  dir[2] = (end[2]-start[2])*invlen;
-
-  // Initialize start point and direction
-  TGeoNode *currentnode = 0;
-  TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
-  if (!startnode) {
-    //AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
-    //		  start[0],start[1],start[2]));
-    return 0.0;
-  }
-  TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
-  lparam[0]   = material->GetDensity();
-  if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
-  lparam[1]   = material->GetRadLen();
-  lparam[2]   = material->GetA();
-  lparam[3]   = material->GetZ();
-  lparam[4]   = length;
-  lparam[5]   = lparam[3]/lparam[2];
-  lparam[6]   = material->GetIntLen();
-  Double_t  n = lparam[0]/lparam[2];
-  Double_t  sigma = 1./(n*lparam[6])/mbarn;
-  if (sigma > mparam[9]) mparam[9]=sigma;
-  if (material->IsMixture()) {
-    TGeoMixture * mixture = dynamic_cast<TGeoMixture*>(material);
-    lparam[5] =0;
-    Double_t sum =0;
-    for (Int_t iel=0;iel<mixture->GetNelements();iel++){
-      sum  += mixture->GetWmixt()[iel];
-      lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
-    }
-    lparam[5]/=sum;
-  }
-
-  // Locate next boundary within length without computing safety.
-  // Propagate either with length (if no boundary found) or just cross boundary
-  gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
-  Double_t step = 0.0; // Step made
-  Double_t snext = gGeoManager->GetStep();
-  // If no boundary within proposed length, return current density
-  if (!gGeoManager->IsOnBoundary()) {
-    mparam[0] = lparam[0];
-    mparam[1] = lparam[4]/lparam[1];
-    mparam[2] = lparam[2];
-    mparam[3] = lparam[3];
-    mparam[4] = lparam[4];
-    mparam[8] = lparam[4]/lparam[6];
-    return lparam[0];
-  }
-  // Try to cross the boundary and see what is next
-  Int_t nzero = 0;
-  while (length>TGeoShape::Tolerance()) {
-    currentnode = gGeoManager->GetCurrentNode();
-    if (snext<2.*TGeoShape::Tolerance()) nzero++;
-    else nzero = 0;
-    if (nzero>3) {
-      // This means navigation has problems on one boundary
-      // Try to cross by making a small step
-      //AliErrorClass("Cannot cross boundary\n");
-      mparam[0] = bparam[0]/step;
-      mparam[1] = bparam[1];
-      mparam[2] = bparam[2]/step;
-      mparam[3] = bparam[3]/step;
-      mparam[5] = bparam[5]/step;
-      mparam[8] = bparam[6];
-      mparam[4] = step;
-      mparam[0] = 0.;             // if crash of navigation take mean density 0
-      mparam[1] = 1000000;        // and infinite rad length
-      return bparam[0]/step;
-    }
-    mparam[6]+=1.;
-    step += snext;
-    bparam[1]    += snext/lparam[1];
-    bparam[2]    += snext*lparam[2];
-    bparam[3]    += snext*lparam[3];
-    bparam[5]    += snext*lparam[5];
-    bparam[6]    += snext/lparam[6];
-    bparam[0]    += snext*lparam[0];
-
-    if (snext>=length) break;
-    if (!currentnode) break;
-    length -= snext;
-    material = currentnode->GetVolume()->GetMedium()->GetMaterial();
-    lparam[0] = material->GetDensity();
-    if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
-    lparam[1]  = material->GetRadLen();
-    lparam[2]  = material->GetA();
-    lparam[3]  = material->GetZ();
-    lparam[5]   = lparam[3]/lparam[2];
-    lparam[6]   = material->GetIntLen();
-    n = lparam[0]/lparam[2];
-    sigma = 1./(n*lparam[6])/mbarn;
-    if (sigma > mparam[9]) mparam[9]=sigma;
-    if (material->IsMixture()) {
-      TGeoMixture * mixture = dynamic_cast<TGeoMixture*>(material);
-      lparam[5]=0;
-      Double_t sum =0;
-      for (Int_t iel=0;iel<mixture->GetNelements();iel++){
-        sum+= mixture->GetWmixt()[iel];
-        lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
-      }
-      lparam[5]/=sum;
-    }
-    gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
-    snext = gGeoManager->GetStep();
-  }
-  mparam[0] = bparam[0]/step;
-  mparam[1] = bparam[1];
-  mparam[2] = bparam[2]/step;
-  mparam[3] = bparam[3]/step;
-  mparam[5] = bparam[5]/step;
-  mparam[8] = bparam[6];
-  return bparam[0]/step;
-}
-
 std::vector<double> GenieGenerator::Rotate(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
 {
   //rotate vector px,py,pz to point at x,y,z at origin.
@@ -412,7 +247,7 @@ Bool_t GenieGenerator::ReadEvent(FairPrimaryGenerator* cpg)
     if (fFirst){
       Double_t bparam=0.;
       Double_t mparam[10];
-      bparam=MeanMaterialBudget(start, end, mparam);
+      bparam=shipgen::MeanMaterialBudget(start, end, mparam);
       cout << "Info GenieGenerator: MaterialBudget " << start[2] << " - "<< end[2] <<  endl;
       cout << "Info GenieGenerator: MaterialBudget " << bparam <<  endl;
       cout << "Info GenieGenerator: MaterialBudget 0 " << mparam[0] <<  endl;
@@ -511,7 +346,7 @@ Bool_t GenieGenerator::ReadEvent(FairPrimaryGenerator* cpg)
         end[1]=tynu*(end[2]-ztarget);
         //cout << "Info GenieGenerator: neutrino xyz-end " << end[0] << "-" << end[1] << "-" << end[2] << endl;
         //get material density between these two points
-        bparam=MeanMaterialBudget(start, end, mparam);
+        bparam=shipgen::MeanMaterialBudget(start, end, mparam);
         //printf("param %e %e %e \n",bparam,mparam[6],mparam[7]);
        }
     }

shipgen/GenieGenerator.h

@@ -39,7 +39,6 @@ class GenieGenerator : public FairGenerator
     endZ        = zE;
   }
   void AddBox(TVector3 dVec, TVector3 box);
-  Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam);
  private:
   std::vector<double> Rotate(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz);