diff --git a/src/core/base2/CMakeLists.txt b/src/core/base2/CMakeLists.txt
index 6be16be07dc49b64105d314291a1ac1d29d49a04..ee41d305af73bef6aabb463a4df483022650d5d9 100644
--- a/src/core/base2/CMakeLists.txt
+++ b/src/core/base2/CMakeLists.txt
@@ -22,7 +22,7 @@ set(lib_files
PMTResponse
Target
TargetMaterial
- SANETargets
+ #SANETargets
#Luminosity
BremsstrahlungRadiator
)
diff --git a/src/core/base2/include/LinkDef.h b/src/core/base2/include/LinkDef.h
index 11bec40a8954230f5de05a9ac0c1775dacd98a18..dfb3082ac202bcddd681288f306a0ed5ecfda6b7 100644
--- a/src/core/base2/include/LinkDef.h
+++ b/src/core/base2/include/LinkDef.h
@@ -35,10 +35,11 @@
#pragma link C++ class insane::Target+;
#pragma link C++ class insane::SimpleTarget+;
#pragma link C++ class insane::SimpleTargetWithWindows+;
-#pragma link C++ class insane::UVAPolarizedAmmoniaTarget+;
-#pragma link C++ class insane::UVACarbonTarget+;
-#pragma link C++ class insane::UVACrossHairTarget+;
-#pragma link C++ class insane::UVAPureHeliumTarget+;
+
+//#pragma link C++ class insane::UVAPolarizedAmmoniaTarget+;
+//#pragma link C++ class insane::UVACarbonTarget+;
+//#pragma link C++ class insane::UVACrossHairTarget+;
+//#pragma link C++ class insane::UVAPureHeliumTarget+;
#pragma link C++ class insane::PMTResponse+;
#pragma link C++ class insane::PMTResponse2+;
diff --git a/src/core/base2/include/SANETargets.h b/src/core/base2/include/SANETargets.h
deleted file mode 100644
index 981d28aa4f2dbb895b60c33be11ad2dd818098bf..0000000000000000000000000000000000000000
--- a/src/core/base2/include/SANETargets.h
+++ /dev/null
@@ -1,75 +0,0 @@
-#ifndef SANETargets_HH
-#define SANETargets_HH 1
-
-#include "Target.h"
-
-
-namespace insane {
-
-/** UVA's polarized Ammonia target.
- *
- * \ingroup Apparatus
- */
-class UVAPolarizedAmmoniaTarget : public Target {
-
- public:
- UVAPolarizedAmmoniaTarget(const char * name = "UVA-polarized-ammonia-target",
- const char * title = "UVA Polarized Ammonia Target",
- Double_t pf = 0.6);
- virtual ~UVAPolarizedAmmoniaTarget();
-
- virtual void DefineMaterials();
-
-
- ClassDef(UVAPolarizedAmmoniaTarget,7)
-};
-
-/** UVA's Carbon target.
- *
- * \ingroup Apparatus
- */
-class UVACarbonTarget : public Target {
-
- public:
- UVACarbonTarget(const char * name = "UVA-Carbon-target", const char * title = "UVA Carbon Target");
- virtual ~UVACarbonTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(UVACarbonTarget,7)
-};
-
-/** Cross Hair.
- *
- * \ingroup Apparatus
- */
-class UVACrossHairTarget : public Target {
-
- public:
- UVACrossHairTarget(const char * name = "UVA-Cross-Hair-target", const char * title = "UVA Cross-Hair Target");
- virtual ~UVACrossHairTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(UVACrossHairTarget,7)
-};
-
-/** UVA's Carbon target.
- *
- * \ingroup Apparatus
- */
-class UVAPureHeliumTarget : public Target {
-
- public:
- UVAPureHeliumTarget(const char * name = "UVA-PureHelium-target", const char * title = "UVA PureHelium Target");
- virtual ~UVAPureHeliumTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(UVAPureHeliumTarget,7)
-};
-
-}
-
-#endif
-
diff --git a/src/core/base2/src/SANETargets.cxx b/src/core/base2/src/SANETargets.cxx
deleted file mode 100644
index 4ffffccec7f4f8ed32a619d62555df0753f7bfdd..0000000000000000000000000000000000000000
--- a/src/core/base2/src/SANETargets.cxx
+++ /dev/null
@@ -1,364 +0,0 @@
-#include "SANETargets.h"
-
-
-namespace insane {
-
-
-//_______________________________________________________________________________
-UVAPolarizedAmmoniaTarget::UVAPolarizedAmmoniaTarget(const char * name, const char * title,Double_t pf): Target(name, title)
-{
-/* TSQLServer * db = TSQLServer::Connect("mysql://quarks.temple.edu","sane","secret");
- TString sql = Form(" SELECT beam_energy,target_angle FROM SANE.run_info WHERE run_number=%d ", run_number);
- SELECT packing_fraction FROM `packing_fractions` WHERE first_run < 72600 AND last_run > 72600*/
-
- fPackingFraction = pf;
- DefineMaterials();
-}
-//_______________________________________________________________________________
-UVAPolarizedAmmoniaTarget::~UVAPolarizedAmmoniaTarget() {
-}
-//_______________________________________________________________________________
-
-void UVAPolarizedAmmoniaTarget::DefineMaterials() {
- // From https://hallcweb.jlab.org/experiments/sane/wiki/index.php/UVa_Polarized_Target#Materials_in_Target
- auto * matH3 = new TargetMaterial("H3", "H3 in Ammonia", 1, 1);
- matH3->fLength = 3.0; //cm
- matH3->fZposition = 0.15; //cm
- matH3->fDensity = 2.0 * 0.0770; // g/cm3, mutiplied by 2 because the denisty assumes a 50% packing fraction
- matH3->fIsPolarized = true;
- matH3->fPackingFraction = fPackingFraction;
-
- auto * matN14 = new TargetMaterial("14N", "14N in Ammonia", 7, 14);
- matN14->fLength = 3.0; //cm
- matN14->fZposition = 0.15; //cm
- matN14->fDensity = 2.0 * 0.3576; //g/cm3 mutiplied by 2 because the denisty assumes a 50% packing fraction
- matN14->fPackingFraction = fPackingFraction;
-
- auto * matHe = new TargetMaterial("He", "LHe around Ammonia beads", 2, 4);
- matHe->fLength = 3.0; //cm
- matHe->fZposition = 0.15; //cm
- matHe->fDensity = 0.1450; // g/cm3
- matHe->fPackingFraction = 1.0 - fPackingFraction;
-
- //TargetMaterial * matLHe = new TargetMaterial("LHe", "Nose Liquid Helium ", 2, 4);
- //matLHe->fLength = 0.50000 * 2.0; //cm
- //matLHe->fDensity = 0.1450; // g/cm3
-
- auto * matLHe1 = new TargetMaterial("LHe1", "Nose Liquid Helium upstream", 2, 4);
- matLHe1->fLength = 0.50000; //cm
- matLHe1->fDensity = 0.1450; //g/cm3
- matLHe1->fZposition = -1.60; //cm
-
- auto * matLHe2 = new TargetMaterial("LHe2", "Nose Liquid Helium downstream", 2, 4);
- matLHe2->fLength = 0.50000; //cm
- matLHe2->fDensity = 0.1450; // g/cm3
- matLHe2->fZposition = 1.90; //cm
-
- //TargetMaterial * matAlEndcap = new TargetMaterial("Al-endcap", "Aluminum target endcap", 13, 27);
- //matAlEndcap->fLength = 2.0 * 0.00381 ; //cm
- //matAlEndcap->fDensity = 2.7; // g/cm3
-
- auto * matAlEndcap1 = new TargetMaterial("Al-endcap1", "Aluminum target endcap upstream", 13, 27);
- matAlEndcap1->fLength = 0.00381 ; //cm
- matAlEndcap1->fDensity = 2.7; // g/cm3
- matAlEndcap1->fZposition = -1.35 ; //cm
-
- auto * matAlEndcap2 = new TargetMaterial("Al-endcap2", "Aluminum target endcap downstream", 13, 27);
- matAlEndcap2->fLength = 0.00381 ; //cm
- matAlEndcap2->fDensity = 2.7; // g/cm3
- matAlEndcap2->fZposition = 1.65; //cm
-
- auto * matAlTail1 = new TargetMaterial("Al-tail1", "Aluminum target tail upstream", 13, 27);
- matAlTail1->fLength = 0.01016 ; //cm
- matAlTail1->fDensity = 2.7; // g/cm3
- matAlTail1->fZposition = -2.05 ; //cm
-
- auto * matAlTail2 = new TargetMaterial("Al-tail2", "Aluminum target tail downstream", 13, 27);
- matAlTail2->fLength = 0.01016 ; //cm
- matAlTail2->fDensity = 2.7; // g/cm3
- matAlTail2->fZposition = 2.35 ; //cm
-
- auto * matAl4kShield1 = new TargetMaterial("Al-4kshield1", "Aluminum target 4K shield upstream", 13, 27);
- matAl4kShield1->fLength = 0.00254 ; //cm
- matAl4kShield1->fDensity = 2.7; // g/cm3
- matAl4kShield1->fZposition = -3.80 ; //cm
-
- auto * matAl4kShield2 = new TargetMaterial("Al-4kshield2", "Aluminum target 4K shield upstream", 13, 27);
- matAl4kShield2->fLength = 0.00254 ; //cm
- matAl4kShield2->fDensity = 2.7; // g/cm3
- matAl4kShield2->fZposition = 3.80 ; //cm
-
- auto * matCu = new TargetMaterial("Cu", "NMR-Cu ", 29, 64);
- matCu->fLength = 0.00050 ; //cm
- matCu->fDensity = 8.9600 ; // g/cm3
- matCu->fZposition = 0.15 ; //cm
-
- auto * matNi = new TargetMaterial("Ni", "NMR-Ni ", 28, 59);
- matNi->fLength = 0.00022 ; //cm
- matNi->fDensity = 8.9020 ; // g/cm3
- matNi->fZposition = 0.15 ; //cm
-
- this->AddMaterial(matH3);
- this->AddMaterial(matN14);
- this->AddMaterial(matHe);
- //this->AddMaterial(matLHe);
- this->AddMaterial(matLHe1);
- this->AddMaterial(matLHe2);
- this->AddMaterial(matAlEndcap1);
- this->AddMaterial(matAlEndcap2);
- this->AddMaterial(matAlTail1);
- this->AddMaterial(matAlTail2);
- this->AddMaterial(matAl4kShield1);
- this->AddMaterial(matAl4kShield2);
- this->AddMaterial(matCu);
- this->AddMaterial(matNi);
-
- matH3->fTGeoVolume->SetTransparency(50);
- matH3->fTGeoVolume->SetFillColor(4);
- matH3->fTGeoVolume->SetLineColor(4);
-
- matN14->fTGeoVolume->SetTransparency(100);
- matHe->fTGeoVolume->SetTransparency(100);
-
- matLHe1->fTGeoVolume->SetTransparency(50);
- matLHe1->fTGeoVolume->SetFillColor(6);
- matLHe1->fTGeoVolume->SetLineColor(6);
-
- matLHe2->fTGeoVolume->SetTransparency(50);
- matLHe2->fTGeoVolume->SetFillColor(6);
- matLHe2->fTGeoVolume->SetLineColor(6);
-
- matAlEndcap1->fTGeoVolume->SetTransparency(30);
- matAlEndcap1->fTGeoVolume->SetFillColor(2);
- matAlEndcap1->fTGeoVolume->SetLineColor(2);
-
- matAlEndcap2->fTGeoVolume->SetTransparency(30);
- matAlEndcap2->fTGeoVolume->SetFillColor(2);
- matAlEndcap2->fTGeoVolume->SetLineColor(2);
-
- matAlTail1->fTGeoVolume->SetTransparency(30);
- matAlTail1->fTGeoVolume->SetFillColor(3);
- matAlTail1->fTGeoVolume->SetLineColor(3);
-
- matAlTail2->fTGeoVolume->SetTransparency(30);
- matAlTail2->fTGeoVolume->SetFillColor(3);
- matAlTail2->fTGeoVolume->SetLineColor(3);
-
-}
-//_______________________________________________________________________________
-
-
-UVACarbonTarget::UVACarbonTarget(const char * name, const char * title): Target(name, title)
-{
- fPackingFraction = 0.6;
- DefineMaterials();
-}
-//_______________________________________________________________________________
-
-UVACarbonTarget::~UVACarbonTarget()
-{
-
-}
-//_______________________________________________________________________________
-
-void UVACarbonTarget::DefineMaterials()
-{
- auto * matH3 = new TargetMaterial("H3", "H3 in Ammonia", 1, 1);
- matH3->fLength = 3.0; //cm
- matH3->fDensity = 2.0 * 0.0770/*/3.0*/; // g/cm3
- matH3->fIsPolarized = true;
- matH3->fPackingFraction = fPackingFraction;
-
- auto * matN14 = new TargetMaterial("14N", "14N in Ammonia", 7, 14);
- matN14->fLength = 3.0; //cm
- matN14->fDensity = 2.0 * 0.3576; // g/cm3
- matN14->fPackingFraction = fPackingFraction;
-
- auto * matHe = new TargetMaterial("He", "LHe around Ammonia beads", 2, 4);
- matHe->fLength = 3.0; //cm
- matHe->fDensity = 0.1450; // g/cm3
- matHe->fPackingFraction = 1.0 - fPackingFraction;
-
- auto * matLHe = new TargetMaterial("LHe", "Liquid Helium", 2, 4);
- matLHe->fLength = 0.50000 * 2.0; //cm
- matLHe->fDensity = 0.1450; // g/cm3
-
- auto * matAlEndcap = new TargetMaterial("Al-endcap", "Aluminum target endcp", 13, 27);
- matAlEndcap->fLength = 2.0 * 0.00381 ; //cm
- matAlEndcap->fDensity = 2.7; // g/cm3
-
- auto * matAlTail = new TargetMaterial("Al-tail", "Aluminum target tail", 13, 27);
- matAlTail->fLength = 2.0 * 0.01016 ; //cm
- matAlTail->fDensity = 2.7; // g/cm3
-
- auto * matAl4kShield = new TargetMaterial("Al-4kshield", "Aluminum target 4K shield", 13, 27);
- matAl4kShield->fLength = 2.0 * 0.00254 ; //cm
- matAl4kShield->fDensity = 2.7; // g/cm3
-
- auto * matCu = new TargetMaterial("Cu", "NMR-Cu ", 29, 64);
- matCu->fLength = 0.00050 ; //cm
- matCu->fDensity = 8.9600 ; // g/cm3
-
- auto * matNi = new TargetMaterial("Ni", "NMR-Ni ", 28, 59);
- matCu->fLength = 0.00022 ; //cm
- matCu->fDensity = 8.9020 ; // g/cm3
-
- this->AddMaterial(matH3);
- this->AddMaterial(matN14);
- this->AddMaterial(matLHe);
- this->AddMaterial(matHe);
- this->AddMaterial(matAlEndcap);
- this->AddMaterial(matAlTail);
- this->AddMaterial(matAl4kShield);
- this->AddMaterial(matCu);
- this->AddMaterial(matNi);
-
-}
-//_______________________________________________________________________________
-
-
-UVACrossHairTarget::UVACrossHairTarget(const char * name, const char * title): Target(name, title)
-{
-
- fPackingFraction = 0.6;
- DefineMaterials();
-}
-//_______________________________________________________________________________
-
-UVACrossHairTarget::~UVACrossHairTarget()
-{
-
-}
-//_______________________________________________________________________________
-
-void UVACrossHairTarget::DefineMaterials()
-{
- auto * matH3 = new TargetMaterial("H3", "H3 in Ammonia", 1, 1);
- matH3->fLength = 3.0; //cm
- matH3->fDensity = 2.0 * 0.0770/3.0; // g/cm3
- matH3->fIsPolarized = true;
- matH3->fPackingFraction = fPackingFraction;
-
- auto * matN14 = new TargetMaterial("14N", "14N in Ammonia", 7, 14);
- matN14->fLength = 3.0; //cm
- matN14->fDensity = 2.0 * 0.3576; // g/cm3
- matN14->fPackingFraction = fPackingFraction;
-
- auto * matHe = new TargetMaterial("He", "LHe around Ammonia beads", 2, 4);
- matHe->fLength = 3.0; //cm
- matHe->fDensity = 0.1450; // g/cm3
- matHe->fPackingFraction = 1.0 - fPackingFraction;
-
- auto * matLHe = new TargetMaterial("LHe", "Liquid Helium", 2, 4);
- matLHe->fLength = 0.50000 * 2.0; //cm
- matLHe->fDensity = 0.1450; // g/cm3
-
- auto * matAlEndcap = new TargetMaterial("Al-endcap", "Aluminum target endcp", 13, 27);
- matAlEndcap->fLength = 2.0 * 0.00381 ; //cm
- matAlEndcap->fDensity = 2.7; // g/cm3
-
- auto * matAlTail = new TargetMaterial("Al-tail", "Aluminum target tail", 13, 27);
- matAlTail->fLength = 2.0 * 0.01016 ; //cm
- matAlTail->fDensity = 2.7; // g/cm3
-
- auto * matAl4kShield = new TargetMaterial("Al-4kshield", "Aluminum target 4K shield", 13, 27);
- matAl4kShield->fLength = 2.0 * 0.00254 ; //cm
- matAl4kShield->fDensity = 2.7; // g/cm3
-
- auto * matCu = new TargetMaterial("Cu", "NMR-Cu ", 29, 64);
- matCu->fLength = 0.00050 ; //cm
- matCu->fDensity = 8.9600 ; // g/cm3
-
- auto * matNi = new TargetMaterial("Ni", "NMR-Ni ", 28, 59);
- matCu->fLength = 0.00022 ; //cm
- matCu->fDensity = 8.9020 ; // g/cm3
-
- this->AddMaterial(matH3);
- this->AddMaterial(matN14);
- this->AddMaterial(matLHe);
- this->AddMaterial(matHe);
- this->AddMaterial(matAlEndcap);
- this->AddMaterial(matAlTail);
- this->AddMaterial(matAl4kShield);
- this->AddMaterial(matCu);
- this->AddMaterial(matNi);
-
-}
-//_______________________________________________________________________________
-
-
-
-//_______________________________________________________________________________
-
-
-UVAPureHeliumTarget::UVAPureHeliumTarget(const char * name, const char * title): Target(name, title)
-{
- fPackingFraction = 0.0;
- DefineMaterials();
-}
-//_______________________________________________________________________________
-
-UVAPureHeliumTarget::~UVAPureHeliumTarget()
-{
-
-}
-//_______________________________________________________________________________
-
-void UVAPureHeliumTarget::DefineMaterials()
-{
- //TargetMaterial * matH3 = new TargetMaterial("H3", "H3 in Ammonia", 1, 1);
- //matH3->fLength = 3.0; //cm
- //matH3->fDensity = 2.0 * 0.0770/*/3.0*/; // g/cm3
- //matH3->fIsPolarized = true;
- //matH3->fPackingFraction = fPackingFraction;
-
- //TargetMaterial * matN14 = new TargetMaterial("14N", "14N in Ammonia", 7, 14);
- //matN14->fLength = 3.0; //cm
- //matN14->fDensity = 2.0 * 0.3576; // g/cm3
- //matN14->fPackingFraction = fPackingFraction;
-
- auto * matHe = new TargetMaterial("He", "LHe around Ammonia beads", 2, 4);
- matHe->fLength = 3.0; //cm
- matHe->fZposition = 0.15; //cm
- matHe->fDensity = 0.1450; // g/cm3
- matHe->fPackingFraction = 1.0 - fPackingFraction;
-
- //TargetMaterial * matLHe = new TargetMaterial("LHe", "Liquid Helium", 2, 4);
- //matLHe->fLength = 0.50000 * 2.0; //cm
- //matLHe->fDensity = 0.1450; // g/cm3
-
- //TargetMaterial * matAlEndcap = new TargetMaterial("Al-endcap", "Aluminum target endcp", 13, 27);
- //matAlEndcap->fLength = 2.0 * 0.00381 ; //cm
- //matAlEndcap->fDensity = 2.7; // g/cm3
-
- //TargetMaterial * matAlTail = new TargetMaterial("Al-tail", "Aluminum target tail", 13, 27);
- //matAlTail->fLength = 2.0 * 0.01016 ; //cm
- //matAlTail->fDensity = 2.7; // g/cm3
-
- //TargetMaterial * matAl4kShield = new TargetMaterial("Al-4kshield", "Aluminum target 4K shield", 13, 27);
- //matAl4kShield->fLength = 2.0 * 0.00254 ; //cm
- //matAl4kShield->fDensity = 2.7; // g/cm3
-
- //TargetMaterial * matCu = new TargetMaterial("Cu", "NMR-Cu ", 29, 64);
- //matCu->fLength = 0.00050 ; //cm
- //matCu->fDensity = 8.9600 ; // g/cm3
-
- //TargetMaterial * matNi = new TargetMaterial("Ni", "NMR-Ni ", 28, 59);
- //matCu->fLength = 0.00022 ; //cm
- //matCu->fDensity = 8.9020 ; // g/cm3
-
- //this->AddMaterial(matH3);
- //this->AddMaterial(matN14);
- this->AddMaterial(matHe);
- //this->AddMaterial(matLHe);
- //this->AddMaterial(matAlEndcap);
- //this->AddMaterial(matAlTail);
- //this->AddMaterial(matAl4kShield);
- //this->AddMaterial(matCu);
- //this->AddMaterial(matNi);
-
-}
-//_______________________________________________________________________________
-
-}
-
diff --git a/src/core/physics/CMakeLists.txt b/src/core/physics/CMakeLists.txt
index f9abf1c861e758c539aba5a292fd7b409b9cc102..2d8bd7a68752bc9a6ec7b4f202bdaf40e21f90dc 100644
--- a/src/core/physics/CMakeLists.txt
+++ b/src/core/physics/CMakeLists.txt
@@ -34,7 +34,6 @@ set(lib_files
#Nucleus
#DilutionFactor
#Acceptance
- EventGenerator
_EG_Event
TargetEventGenerator
#TargetMaterial
@@ -42,10 +41,9 @@ set(lib_files
#BremsstrahlungRadiator
#Target
- CLASTargets
+ #CLASTargets
CrossSectionInfo
FunctionManager
- DiffXSec
CompositeDiffXSec
InclusiveDISXSec
DiffXSecKinematicKey
@@ -57,9 +55,7 @@ set(lib_files
eInclusiveElasticDiffXSec
pInclusiveElasticDiffXSec
XSections
- InclusiveDiffXSec
InclusiveBornDISXSec
- ExclusiveDiffXSec
MAIDKinematicKey
MAIDInclusiveDiffXSec
MAIDNucleusInclusiveDiffXSec
@@ -112,9 +108,6 @@ set(lib_files
AmrounFormFactors
BilenkayaFormFactors
BeamSpinAsymmetry
- PhaseSpaceSampler
- PhaseSpaceVariable
- PhaseSpace
PartonDistributionFunctions
PolarizedPartonDistributionFunctions
StatisticalQuarkFits
diff --git a/src/core/physics/include/CLASTargets.h b/src/core/physics/include/CLASTargets.h
deleted file mode 100644
index b50a01f010d35c187333b63e6651570dd9f2ee09..0000000000000000000000000000000000000000
--- a/src/core/physics/include/CLASTargets.h
+++ /dev/null
@@ -1,83 +0,0 @@
-#ifndef InSANE_CLASTargets_HH
-#define InSANE_CLASTargets_HH
-
-#include "Target.h"
-
-namespace insane {
- namespace physics {
-
- class Helium4GasTarget : public Target {
- public:
- Helium4GasTarget(const char * name = "Helium4GasTarget",
- const char * title = "Helium4 Gas Target" );
- virtual ~Helium4GasTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(Helium4GasTarget,1)
- };
-
- class DeuteriumGasTarget : public Target {
- public:
- DeuteriumGasTarget(const char * name = "DeuteriumGasTarget",
- const char * title = "Deuterium Gas Target" );
- virtual ~DeuteriumGasTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(DeuteriumGasTarget,1)
- };
-
- class HydrogenGasTarget : public Target {
- public:
- HydrogenGasTarget(const char * name = "HydrogenGasTarget",
- const char * title = "Hydrogen Gas Target" );
- virtual ~HydrogenGasTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(HydrogenGasTarget,1)
- };
-
- class LH2Target : public Target {
- public:
- LH2Target(const char * name = "LH2Target",
- const char * title = "Liquid H2 Target" );
- virtual ~LH2Target();
-
- virtual void DefineMaterials();
-
- ClassDef(LH2Target,1)
- };
-
-
- // TDIS targets (not really CLAS)
-
- class TDISDeuteriumGasTarget : public Target {
- public:
- TDISDeuteriumGasTarget(const char * name = "DeuteriumGasTarget",
- const char * title = "Deuterium Gas Target" );
- virtual ~TDISDeuteriumGasTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(TDISDeuteriumGasTarget,1)
- };
-
- class TDISHydrogenGasTarget : public Target {
- public:
- TDISHydrogenGasTarget(const char * name = "TDISHydrogenGasTarget",
- const char * title = "TDISHydrogen Gas Target" );
- virtual ~TDISHydrogenGasTarget();
-
- virtual void DefineMaterials();
-
- ClassDef(TDISHydrogenGasTarget,1)
- };
-
-
- }
-}
-
-#endif
-
diff --git a/src/core/physics/include/DiffXSec.h b/src/core/physics/include/DiffXSec.h
deleted file mode 100644
index 73743e27ac2ead4c3ac6671df365578e8670f2ba..0000000000000000000000000000000000000000
--- a/src/core/physics/include/DiffXSec.h
+++ /dev/null
@@ -1,378 +0,0 @@
-#ifndef InSANEDiffXSec_H
-#define InSANEDiffXSec_H 1
-
-#include <ostream>
-#include "TFoam.h"
-#include "TFoamIntegrand.h"
-#include "TMath.h"
-#include "TVector3.h"
-#include "Math/IFunction.h"
-#include "Math/Functor.h"
-#include "TString.h"
-#include "Math/Integrator.h"
-#include "Math/IntegratorMultiDim.h"
-#include "Math/AllIntegrationTypes.h"
-#include "Math/GaussIntegrator.h"
-
-#include "Nucleus.h"
-#include "PhaseSpace.h"
-#include "Particle.h"
-//#include "BETAG4MonteCarloEvent.h"
-#include "TargetMaterial.h"
-#include "FormFactors.h"
-#include "StructureFunctions.h"
-#include "PolarizedStructureFunctions.h"
-
-
-namespace insane {
-namespace physics {
-
- //class StructureFunctions;
- //class PolarizedStructureFunctions;
-
- /** Abstract base class for a Differential Cross Section.
- * \f$ \frac{d^N\sigma}{dx_1 dx_2 ... dx_N} \f$
- *
- * Note, that the solid angle \f$ d\Omega = sin(\theta)d\theta d\phi \f$
- * requires that the flagg fIncludeJacobian be set to true to include the
- * \f$sin(\theta)\f$ factor in calculated result.
- *
- * Units used:
- * - GeV
- * - radian
- * - steradian
- * - nanobarn
- *
- * Methods to override:
- * - EvaluateXSec
- * - InitPhaseSpace
- * - GetDependentVariables
- *
- * See xsec_ids.md for XS ID numbering convetion.
- *
- * Since we are often using natural units to do a calculation, i.e. h=1 and c=1,
- * we need to multiply cross sections by
- * (hc)^2=(197.33 MeV (-fm))^2 = 0.038939129 GeV^2 fm^2 =
- * 0.00038939129 GeV^2 barn = 0.38939129 GeV^2 mbarn
- *
- * \ingroup EventGen
- * \ingroup xsections
- *
- */
- class DiffXSec : public TFoamIntegrand, public ROOT::Math::IBaseFunctionMultiDim {
-
- private:
- Bool_t fIncludeJacobian;
- Bool_t fIsModified;
- Bool_t fUsePhaseSpace; // When false xsec evaluation should ignore values of phase space limits
- PhaseSpace * fPhaseSpace; //->
- mutable Double_t fVars[30]; //->
-
- protected:
- Int_t fID; ///< ID unique to the cross section. Used in e.g. for identifying source of thrown event.
- Int_t fBaseID; ///< ID which is used to set fID to differentiate cross sections for final state particles
- std::vector<Int_t> fPIDs; ///< PIDs of final state particles, should have fnParticles entries
- Int_t fnDim; ///< Dimension of differential.
- Int_t fnParticles; ///< Number of particles in the final state.
- Int_t fnPrimaryParticles; ///< Number of particles printed to lund file.
- std::vector<Int_t> fnParticleVars; ///< Number of PS variables associated with particle
- TList fAddedXSections; ///< Currently not used!!
- TList fParticles; ///< List of particles that are thrown for reaction.
- TList fVariableNames; ///< List of variable names
- Double_t fBeamEnergy; ///< Incoming beam energy
- Double_t fTotalXSec; ///< Total cross-section over phase space.
- Double_t fHelicity; ///< Beam helicity, 0 if unpolarized
- TVector3 fTargetPol; ///< Target Polariztion vector (in the lab)
-
- Nucleus fTargetNucleus; ///< Isotope used to calculate cross section
- TargetMaterial fTargetMaterial; ///< Target material associated with cross section
- Int_t fMaterialIndex; ///< Target material index associated with cross section
-
- std::vector<Double_t> fScaleFactors; ///< Scale factors are set so that the random variables fall between 0 and 1
- ///< Scale factor for multilying random variable such that it is less than 1.
- std::vector<Double_t> fOffsets; ///< Off sets for negative numbers
-
- mutable Double_t fNormalizedVariables[30]; //-> should be more clever here
- mutable Double_t fUnnormalizedVariables[30]; //-> should be more clever here
- mutable Double_t fDependentVariables[30]; //-> should be more clever here
-
- TString fTitle;
- TString fPlotTitle;
- TString fLabel;
- TString fUnits;
-
- StructureFunctions * fStructureFunctions; //!
- PolarizedStructureFunctions * fPolarizedStructureFunctions; //!
- FormFactors * fFormFactors; //!
-
- public:
- DiffXSec();
- DiffXSec(const DiffXSec& old);
- DiffXSec& operator=(const DiffXSec& old);
-
- //DiffXSec(const DiffXSec&) = default; // Copy constructor
- DiffXSec(DiffXSec&&) = default; // Move constructor
- //DiffXSec& operator=(const DiffXSec&) & = default; // Copy assignment operator
- DiffXSec& operator=(DiffXSec&&) & = default; // Move assignment operato
-
- virtual ~DiffXSec();
-
- virtual void SetTitle(const char * t){ fTitle = t;}
- const char * GetTitle() const {return( fTitle.Data());}
- virtual void SetPlotTitle(const char * t){ fPlotTitle = t;}
- const char * GetPlotTitle() const {return( fPlotTitle.Data());}
- virtual void SetLabel(const char * t){ fLabel = t;}
- const char * GetLabel(){ return(Form("%s (%s)",fLabel.Data(),fUnits.Data() )); }
- virtual void SetUnits(const char * t){ fUnits = t;}
- const char * GetUnits() const {return( fUnits.Data());}
-
- Int_t GetID() const { return fID; }
-
- virtual void UsePhaseSpace(Bool_t val = true) {fUsePhaseSpace = val; }
-
- virtual void SetTargetMaterial(const TargetMaterial& mat){ fTargetMaterial = mat; fTargetNucleus = fTargetMaterial.GetNucleus(); }
- const TargetMaterial& GetTargetMaterial() const {return fTargetMaterial;}
-
- virtual void SetTargetMaterialIndex(Int_t i){ fMaterialIndex = i ; }
- Int_t GetTargetMaterialIndex() const {return fMaterialIndex;}
-
- /** A method that should be overriden which creates a phase space and phase
- * space variables.
- * It should end with SetPhaseSpace() which sets the member
- * fPhaseSpace and does some checks.
- */
- virtual void InitializePhaseSpaceVariables() {}
-
- /** Virtual method to define the random event from the variables provided.
- * This is the transition point between the phase space variables and
- * the particles. The argument should be the full list of variables returned from
- * GetDependentVariables.
- */
- virtual void DefineEvent(Double_t * vars ){}
-
- /** Clone the cross section. Note: you should provide both versions of clone
- * with the 0 argument version that calls the 1 arg version.
- */
- virtual DiffXSec* Clone(const char * newname) const ;
- virtual DiffXSec* Clone() const { return( Clone("") ); }
-
- /** Creates an instance of Particle for each final state particle.
- * Using the array fPIDs it sets the appropriate PDG encoding for
- * each final state particle.
- */
- virtual void InitializeFinalStateParticles();
-
- /** The list of final state particles. */
- TList * GetParticles();
- TParticlePDG * GetParticlePDG(int i = 0) const;
-
- /** Returns the number of PS vars for the ith particle. */
- Int_t GetNParticleVars(unsigned int i=0){ if( i>=fnParticleVars.size() ) return(0); return(fnParticleVars[i]);}
-
- const Nucleus& GetTargetNucleus() const { return fTargetNucleus; }
- virtual void SetTargetNucleus(const Nucleus & targ){ fTargetNucleus = targ; }
-
- Double_t GetZ() const { return( Double_t( fTargetNucleus.GetZ() )) ;}
- Double_t GetA() const { return( Double_t( fTargetNucleus.GetA() )) ;}
- Double_t GetN() const { return( Double_t( fTargetNucleus.GetN() )) ;}
-
-
- Double_t GetHelicity() const { return fHelicity; }
- virtual void SetHelicity(Double_t h) {
- // make sure to add call to this in override
- if(TMath::Abs(h) > 1.0){
- Error("SetHelicity","Argument too big: |h|<=1.");
- }else{
- fHelicity = h;
- }
- }
- const TVector3& GetTargetPolarization(){ return fTargetPol; }
- virtual void SetTargetPolarization(const TVector3 &P){
- // make sure to add call to this in override
- fTargetPol = P;
- if(P.Mag() > 1.0) {
- Error("SetTargetPolarization","Magnitude too big. Setting to 1.");
- fTargetPol.SetMag(1.0);
- }
- }
-
- virtual void SetTargetPolDirection(Double_t theta, Double_t phi = 0.0){
- TVector3 pol = fTargetPol;
- if(pol.Mag() == 0.0) {
- Warning("SetTargetPolDirection","Target Polarization was 0. Setting to 1.");
- pol.SetMag(1.0);
- }
- pol.SetTheta(theta);
- pol.SetTheta(phi);
- SetTargetPolarization(pol);
- }
- virtual void SetPolarizations(Double_t pe, Double_t pt){
- // make sure to add call to this in override or repeat this
- TVector3 pol = fTargetPol;
- pol.SetMag(pt);
- //fHelicity = pe;
- SetHelicity(pe);
- SetTargetPolarization(pol);
- }
-
- /** Includes sin(theta) term for solid angle integration.
- * This is set automatically when the cross section is added to a phase space sampler,
- * otherwise it should be ignored.
- */
- Bool_t IncludeJacobian() const { return fIncludeJacobian ; }
- void SetIncludeJacobian(Bool_t v = true){ fIncludeJacobian = v; }
-
- /** Virtual method that returns the calculated values of dependent variables. This method
- * should be overridden for exclusive cross sections in order to conserve momentum and energy!
- * This is called from Density(const Double_t *y), where argument is an array
- * independent random variables only.
- *
- * For example, in mott scattering (elastic) there is really only two random variables,
- * the the scattered angles theta and phi. The rest can be calculated from these two angles
- * (assuming the beam energy is known).
- *
- * This implementation assumes that all variables are independent.
- */
- virtual Double_t * GetDependentVariables(const Double_t * x) const;
-
- const Double_t* DependentVariables() const { return fDependentVariables; }
-
- /** @name Evaluating the cross section
- * Various methods used to evaluate the cross sections.
- * Some are required by base classes some are just convenient
- * @{
- */
-
- /** The main virtual method which implements the cross section.
- *
- * The argument array should be passed as the results of GetDependentVariables() (especially for
- * exclusive processes.
- */
- virtual Double_t EvaluateXSec(const Double_t * x) const ;
-
- /** Calculate the Jacobian determinant.
- * Returns the jacobian needed to convert the cross section
- * into the right form for integrating over the independent
- * variables for the particles defined in the phase space.
- * Note: make sure the jacobian is used if IncludeJacobian
- * is set, otherwise, using the cross section with an event
- * generator will produce the wrong total cross section.
- * Takes the same arguments as EvaluateXSec which is the
- * result of GetDependentVariables.
- */
- virtual Double_t Jacobian(const Double_t * x) const;
-
- virtual Double_t GetBeamSpinAsymmetry(const Double_t * x) const { return 0.0; }
-
- /** Used to evaluate the base cross section.
- * Useful for radiative corrections where EvaluteXSec implements
- * methods that want to call the base. Override this so RADCOR and POLRAD
- * Use the base class.
- */
- virtual Double_t EvaluateBaseXSec(const Double_t * x) const { return EvaluateXSec(x); }
-
- /** Pure virtual method from TFoamIntegrand method used by TFoam.
- * \f$ \frac{d^3\sigma}{d\theta d\phi d\omega} \f$
- * DiffXSec::Density() is sampled by TFoam with 0 < x < 1
- * and dimension equal to the differntial dimension (not the number of phase space variables).
- */
- virtual Double_t Density(Int_t ndim, Double_t * x) ;
-
- /** Evaluates the cross section with the current values of
- * the phase space variables.
- */
- virtual Double_t EvaluateCurrent() const ;
-
- double DoEval(const double* x) const { return (double)EvaluateXSec((Double_t*) x) ; }
- double operator()(double *x, double *p) { return (double)EvaluateXSec((Double_t*) x); }
- unsigned int NDim() const { return fnDim; }
-
- //@}
-
- /** Refresh the cross section.
- * Recalculates the scale factors for each variable, which puts the
- * values between 0 and 1 (as required for TFoam).
- * Called from PhaseSpaceSampler::Refresh().
- */
- void Refresh();
-
- /** Used when object is read from disk.
- * If there are special initializations that need to be done for persistency, this method should
- * be implemented.
- */
- virtual Int_t InitFromDisk(){ return 0; }
-
- /** Sets the phase space and scale factors.
- * The argument PS object is owned by this class.
- * Make sure phase space has had all variables added prior to adding.
- * Initializes the final state particles ( see InitializeFinalStateParticles() )
- */
- virtual void SetPhaseSpace(PhaseSpace * ps);
- PhaseSpace * GetPhaseSpace() const { return(fPhaseSpace); }
-
- /** Sets the integral of the XSection over the phase space*/
- void SetTotalXSec(Double_t norm) { fTotalXSec = norm; }
- Double_t GetTotalXSec(){ return fTotalXSec;}
-
- Double_t GetPDFNormalization() const { return(fTotalXSec); } ///< Deprecated
- Double_t GetIntegratedCrossSection() const { return(fTotalXSec); }
- Double_t CalculateTotalCrossSection();
-
- /** Set Beam Energy in units of GeV */
- virtual void SetBeamEnergy(Double_t ebeam) { fBeamEnergy = ebeam; fIsModified = true; }
- virtual Double_t GetBeamEnergy() const { return(fBeamEnergy); }
-
- /** Returns the variables with ther physical values. */
- Double_t * GetUnnormalizedVariables(const Double_t * x);
-
- /** Returns the variables scaled between 0 and 1 (for TFoam). */
- Double_t * GetNormalizedVariables(const Double_t * x);
-
- /** Returns true if the set of variables fall within the phase space.
- * See DiffXSec::SetPhaseSpace().
- * Calls SetEventValues() which takes the raw array of values and sets the
- * variable's datamember, fCurrentValue.
- */
- bool VariablesInPhaseSpace(const Int_t ndim, const Double_t * x) const;
-
- virtual void Print(std::ostream &stream) const ;
- virtual void Print(const Option_t * opt = "") const ; // *MENU*
- virtual void PrintTable(); // *MENU*
- virtual void PrintTable(std::ostream &out);
- virtual void PrintGrid(std::ostream &out,Int_t N = 20);
-
- /** Returns the PDG particle encoding.
- * Note: For multi particle production cross sections you should use
- * the argument to get that particles PDG encoding.
- */
- Int_t GetParticleType(Int_t part = 0) const ;
-
- virtual void SetParticleType( Int_t pdgcode, Int_t part = 0) ;
- virtual void SetProductionParticleType( Int_t PDGcode, Int_t part = 0);
-
- void PrintFinalStatePIDs() const ;
-
- bool IsModified() const { return fIsModified; }
- void SetModified(bool val = false) { fIsModified = val; }
-
-
- /** Set the unpolarized structure functions to be used to calculate W1,W2,F1,F2, etc... */
- virtual void SetUnpolarizedStructureFunctions(StructureFunctions * sf) { fStructureFunctions = sf; }
- StructureFunctions * GetUnpolarizedStructureFunctions() const { return(fStructureFunctions); }
-
- /** Set the polarized structure functions to be used to calculate G1,G2,g1,g2, etc... */
- virtual void SetPolarizedStructureFunctions(PolarizedStructureFunctions * sf) { fPolarizedStructureFunctions = sf; }
- PolarizedStructureFunctions * GetPolarizedStructureFunctions() const { return(fPolarizedStructureFunctions); }
-
- /** Set the form factors */
- virtual void SetFormFactors(FormFactors * ff) { fFormFactors = ff; }
- FormFactors * GetFormFactors() const { return(fFormFactors); }
-
- ClassDef(DiffXSec,6)
- };
-
-}
-}
-
-#endif
-
diff --git a/src/core/physics/include/EventGenerator.h b/src/core/physics/include/EventGenerator.h
deleted file mode 100644
index dd48754e9ba5c16cf09210c5b9655eec615dd9ae..0000000000000000000000000000000000000000
--- a/src/core/physics/include/EventGenerator.h
+++ /dev/null
@@ -1,184 +0,0 @@
-#ifndef EventGenerator_H
-#define EventGenerator_H 1
-
-#include "TNamed.h"
-#include "PhaseSpaceSampler.h"
-#include "TList.h"
-//#include "BETAG4MonteCarloEvent.h"
-#include "TRandom3.h"
-#include "InSANEMathFunc.h"
-//#include "RunManager.h"
-#include "TLorentzVector.h"
-#include <iostream>
-#include "Particle.h"
-#include "TParticlePDG.h"
-#include "TClonesArray.h"
-#include <vector>
-#include "_EG_Event.h"
-
-/*! \page simulation Simulation Tools
-
- \section montecarlo Monte Carlo Event Generation
- An example that creates a few different cross sections and puts them together in an event generator.
- @include event_generator.cxx
-
- */
-
-
-namespace insane {
-namespace physics {
-
- class DiffXSec;
-
-/** Event generator for Montecarlo.
- * Add phase space samplers which have mutually exclusive phase spaces.
- * For example, you can add inclusive pi0 and inclusive e-, but NOT,
- * inclusive e- and a helicity dependent cross section.
- *
- * The basic operation of an event generator goes like this...
- * First the phase space samplers are added (see PhaseSpaceSampler for more).
- * Then the total cross-section, \f$ \sigma_{T} \f$, is calculated in order to normalize the total probability.
- * Also, the total cross-section for each sampler, \f$ \sigma_{T}^{i} \f$ added.
- * When an event is generated using GenerateEvent() a random sampler is selected
- * with probability \f$ P(i) = \sigma_{T}^{i}/\sigma_{T} \f$. If there is one sampler the probabilty is 1.
- * The event is then generated using PhaseSpaceSampler::GenerateEvent() which returns a TList of Particle.
- * For details of how the kinematic variables are sampled, see PhaseSpaceSampler.
- *
- * The volume for uniformly sampling the vertex can be either a box or
- * a cylinder (rastered beam).
- *
- * When creating an event generator, the Initialize() method should be implemented.
- * In Initialize() you should attach the appropriate phase space samplers. Here you
- * can also create cross sections and phase spaces for said phase space samplers.
- * There is no need to call Refresh() here. This is done just before using the event generator.
- *
- *
- * \ingroup EventGen
- */
-class EventGenerator : public TNamed {
-
- private:
- bool fIsModified;
-
- protected:
- Double_t fTotalXSection;
- TList fSamplers;
- Double_t fBeamEnergy;
-
- TRandom * fRandom; //!
- Double_t * fEventArray; //!
- TList * fParticles; //!
- PhaseSpaceSampler * fSampler; //! currently randomly selected sampler
-
- public :
- bool fIsBeamRastered;
- Double_t fSlowRasterSize;
- Double_t fdxVertex;
- Double_t fdyVertex;
- Double_t fdzVertex;
-
- _EG_Event fEG_Event;
-
- public:
- EventGenerator(const char * name = "eventGen", const char * title = "Event Gen") ;
- EventGenerator( const EventGenerator & eg );
- EventGenerator& operator=(const EventGenerator& v);
-
- virtual ~EventGenerator() ;
-
- /** Set the beam energy for EVERY cross section. */
- void SetBeamEnergy(Double_t E0);
- Double_t GetBeamEnergy() const { return fBeamEnergy; }
-
- /** Virtual method which should create the specific samplers, cross-sections, etc...
- * for the specific event generator implementation.
- * Otherwise it is the same as refresh.
- */
- virtual void Initialize() {
- //std::cout << "EventGenerator::Initialize" << std::endl;
- /// Add here.
- Refresh();
- }
- void InitFromDisk();
-
- void Refresh(bool print = false); ///< Called after phase spaces have been modified.
-
- void ListSamplers(); ///< Lists each phase space sampler
-
- void ListPhaseSpaceVariables(); ///< Prints the phase space variables associated with each phase space sampler.
-
- void Print(const Option_t * opt = "") const ;
- void Print(std::ostream& stream) const ;
-
- Double_t CalculateTotalCrossSection(); ///< Adds all the sampler's total Xsections together
- Double_t GetTotalCrossSection() const { return fTotalXSection; }
- Double_t GetTotalCrossSectionForParticle(Int_t pid =0) ;
-
- /** Add a configured phase space sampler to the event generator.
- */
- void AddSampler(PhaseSpaceSampler * ps) ;
-
- /** Clear the Event generator. Removes all samplers and resets the total Xsection to zero.
- */
- void Clear(Option_t * opt = "") ;
-
- /** Get a random sampler with a probablity proportional to it's total cross section.
- * Used to select the phase space sampler from which the event will generated.
- */
- PhaseSpaceSampler * GetRandomSampler() ;
-
- TList * GetSamplers(); ///< Returns pointer to the list of PS samplers
-
- /** Returns a list with the variable names that match.
- * Looks through all sameplers/phase spaces.
- */
- TList * GetPSVariables(const char * var );
-
- /** Generates a random vertex.
- * The argument is
- */
- TLorentzVector GetRandomVertex(Int_t i=-1) ;
-
- /** Checks each cross-section their phase-space and its associated variables by using
- * their methods IsModified().
- * Use as a check, but not at the event level!
- */
- bool NeedsRefreshed() ;
-
- bool IsModified() const { return fIsModified; }
- void SetModified(bool val = true) { fIsModified = val; }
-
- void PrintXSecSummary(std::ostream& s = std::cout);
-
- /** Generate an event.
- * This method retuns a list of Particles to be thrown.
- * Get a random sampler and use it to generate a random event
- *
- * \todo Need to translate P.S. variable (which are arbitrary quantities)
- * into well defined Particle
- */
- TList * GenerateEvent() ;
-
- Double_t * GetEventArray() { return(fEventArray); }
-
- void LundFormat( std::ostream& s = std::cout );
-
- void LundHeaderFormat(
- _EG_Event& eg_event,
- std::ostream& s = std::cout) ;
-
- void LundEventFormat(
- int i,
- _EG_Event& eg_event,
- std::ostream& s = std::cout) ;
-
-
- ClassDef(EventGenerator, 4)
-};
-
-}
-}
-
-
-#endif
-
diff --git a/src/core/physics/include/ExclusiveDiffXSec.h b/src/core/physics/include/ExclusiveDiffXSec.h
deleted file mode 100644
index 21d68aa24be8be9fd9a4ac6a89636e98df27f442..0000000000000000000000000000000000000000
--- a/src/core/physics/include/ExclusiveDiffXSec.h
+++ /dev/null
@@ -1,293 +0,0 @@
-#ifndef ExclusiveDiffXSec_H
-#define ExclusiveDiffXSec_H 1
-#include "DiffXSec.h"
-#include "TFoamIntegrand.h"
-#include "TMath.h"
-#include "TCanvas.h"
-#include "TGraph.h"
-#include "TGraph2D.h"
-#include "FunctionManager.h"
-#include "TParticlePDG.h"
-#include "TDatabasePDG.h"
-#include "TVector3.h"
-
-
-namespace insane {
-namespace physics {
-/** Base class for an Exclusive Differential Cross Section.
- * Examples include the Mott cross section, (e,e'p) elastic scattering, exclusive meson productions etc...
- * Uses units of GeV, radians, and nanobarns.
- *
- * \ingroup exclusiveXSec
- */
-class ExclusiveDiffXSec : public DiffXSec {
- protected:
-
- //Double_t * fTheta_p; //!
- //Double_t * fPhi_p; //!
- //Double_t * fEnergy_p; //!
-
- public:
- ExclusiveDiffXSec() ;
- virtual ~ExclusiveDiffXSec();
-
- ExclusiveDiffXSec(const ExclusiveDiffXSec& old) : DiffXSec(old) { (*this) = old; }
-
- ExclusiveDiffXSec& operator=(const ExclusiveDiffXSec& old) {
- if (this != &old) { // make sure not same object
- this->DiffXSec::operator=(old);
- }
- return *this; // Return ref for multiple assignment
- }
- //ExclusiveDiffXSec(const ExclusiveDiffXSec&) = default; // Copy constructor
- ExclusiveDiffXSec(ExclusiveDiffXSec&&) = default; // Move constructor
- //ExclusiveDiffXSec& operator=(const ExclusiveDiffXSec&) & = default; // Copy assignment operator
- ExclusiveDiffXSec& operator=(ExclusiveDiffXSec&&) & = default; // Move assignment operato
- //virtual ~ExclusiveDiffXSec() { }
-
- /** \todo how to add cross sections ? */
- ExclusiveDiffXSec operator+ (ExclusiveDiffXSec& right) {
- ExclusiveDiffXSec temp;
- return (temp);
- }
-
- /** Virtual method that returns the calculated values of dependent variables. This method
- * should be overridden for exclusive cross sections in order to conserve momentum and energy!
- * This is called from Density(const Double_t *y), where argument is an array of only the
- * independent random variables.
- *
- * For example, in mott scattering (elastic) there is really only two random variables,
- * the the scattered angles theta and phi. The rest can be calculated from these two angles
- * (assuming the beam energy is known).
- *
- */
- virtual Double_t * GetDependentVariables(const Double_t * x) const {
- std::cout << "ExclusiveDiffXSec::GetDependentVariables" << std::endl;
- TVector3 k1(0, 0, fBeamEnergy); // incident electron
- TVector3 k2(0, 0, 0); // scattered electron
- double eprime = fBeamEnergy / (1.0 + 2.0 * fBeamEnergy / M_p/GeV * TMath::Power(TMath::Sin(x[0] / 2.0), 2)) ;
- k2.SetMagThetaPhi(eprime, x[0], x[1]);
- TVector3 p2 = k1 - k2; // recoil proton
- fDependentVariables[0] = eprime;
- fDependentVariables[1] = x[0];
- fDependentVariables[2] = x[1];
- fDependentVariables[3] = p2.Mag();
- fDependentVariables[4] = p2.Theta();
- fDependentVariables[5] = p2.Phi();
- return(fDependentVariables);
- }
-
- virtual void DefineEvent(Double_t * vars);
-
- /** TFoamIntegrand method used by TFoam.
- * Density() is sampled by TFoam with 0 < x < 1
- * \frac{d^3\sigma}{d\theta d\phi d\omega}
- */
- virtual Double_t Density(Int_t ndim, Double_t * x) {
- return(EvaluateXSec(GetDependentVariables(GetUnnormalizedVariables(x))));
- }
-
- /** Initialize a phase space for (e,e'p) */
- virtual void InitializePhaseSpaceVariables();
-
- /** Virtual Method used to evaluate cross section */
- virtual Double_t EvaluateXSec(const Double_t * x) const;
-
- /** Needed by ROOT::Math::IBaseFunctionMultiDim. */
- unsigned int NDim() const { return fnDim; }
-
- /** Needed by ROOT::Math::IBaseFunctionMultiDim. */
- double DoEval(const double * x) const { return (double)EvaluateXSec((Double_t*) x); }
-
- virtual ExclusiveDiffXSec* Clone(const char * newname) const {
- std::cout << "ExclusiveDiffXSec::Clone()\n";
- auto * copy = new ExclusiveDiffXSec();
- (*copy) = (*this);
- return copy;
- }
- virtual ExclusiveDiffXSec* Clone() const { return( Clone("") ); }
-
- /** Set the unpolarized structure functions to be used to calculate W1,W2,F1,F2, etc...
- */
- void SetUnpolarizedStructureFunctions(StructureFunctions * sf) { fStructureFunctions = sf; }
- StructureFunctions * GetUnpolarizedStructureFunctions() { return(fStructureFunctions); }
-
- /** Set the polarized structure functions to be used to calculate G1,G2,g1,g2, etc...
- */
- void SetPolarizedStructureFunctions(PolarizedStructureFunctions * sf) { fPolarizedStructureFunctions = sf; }
- PolarizedStructureFunctions * GetPolarizedStructureFunctions() { return(fPolarizedStructureFunctions); }
-
-
- /** @name Useful functions for plotting...
- * @{
- */
- /** For using with TF1 functions */
- virtual double EnergyDepXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = x[0];
- y[1] = p[0];
- y[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- /** For using with TF1 functions. */
- virtual double EnergyFractionDepXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = x[0] * (fBeamEnergy);
- y[1] = p[0];
- y[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- /** Cross section as a function of energy.
- * x[0] = energy
- * p[0] = theta
- * p[1] = phi
- */
- virtual double EnergyDependentXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = x[0];
- y[1] = p[0];
- y[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- /** Cross section as a function of momentum.
- * x[0] = momentum
- * p[0] = theta
- * p[1] = phi
- */
- virtual double MomentumDependentXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = x[0];
- y[1] = p[0];
- y[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- /** Cross section as a function of polar angle .
- * x[0] = theta
- * p[0] = energy
- * p[1] = phi
- */
- virtual double PolarAngleDependentXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = x[0];
- y[1] = p[0];
- y[2] = p[0];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- /** Cross section as a function of polar angle .
- * x[0] = phi
- * p[0] = energy
- * p[1] = theta
- */
- virtual double AzimuthalAngleDependentXSec(double *x, double *p) {
- Double_t y[3];
- y[0] = p[0];
- y[1] = p[1];
- y[2] = x[0];
- return(EvaluateXSec(GetDependentVariables(y)));
- }
-
- //@}
-
- ClassDef(ExclusiveDiffXSec, 1)
-};
-
-
-
-/** A flat exclusive cross section for (e,e'p)
- *
- * \ingroup exclusiveXSec
- */
-class FlatExclusiveDiffXSec : public ExclusiveDiffXSec {
- public:
- FlatExclusiveDiffXSec() { }
- virtual ~FlatExclusiveDiffXSec() { }
- virtual FlatExclusiveDiffXSec* Clone(const char * newname) const {
- std::cout << "FlatExclusiveDiffXSec::Clone()\n";
- auto * copy = new FlatExclusiveDiffXSec();
- (*copy) = (*this);
- return copy;
- }
- virtual FlatExclusiveDiffXSec* Clone() const { return( Clone("") ); }
-
-
- /** Virtual method that returns the calculated values of dependent variables. This method
- * should be overridden for exclusive cross sections in order to conserve momentum and energy!
- *
- * For example, in mott scattering (elastic) there is really only two random variables,
- * the the scattered angles theta and phi. The rest can be calculated from these two angles
- * (assuming the beam energy is known).
- *
- */
- virtual Double_t * GetDependentVariables(const Double_t * x) const {
- // A faster way would be to solve the kinematics and plug the solutions in explicitly.
- TVector3 k1(0, 0, fBeamEnergy); // incident electron
- TVector3 k2(0, 0, 0); // scattered electron
- double eprime = fBeamEnergy / (1.0 + 2.0 * fBeamEnergy / M_p/GeV * TMath::Power(TMath::Sin(x[0] / 2.0), 2)) ;
- k2.SetMagThetaPhi(eprime, x[0], x[1]);
- TVector3 p2 = k1 - k2; // recoil proton
- fDependentVariables[0] = eprime;
- fDependentVariables[1] = x[0];
- fDependentVariables[2] = x[1];
- fDependentVariables[3] = p2.Mag();
- fDependentVariables[4] = p2.Theta();
- fDependentVariables[5] = p2.Phi();
- p2.Print();
- k2.Print();
- return(fDependentVariables);
- }
-
-
- /** Evaluate Cross Section. Flat cross section returns 1 */
- virtual Double_t EvaluateXSec(const Double_t * x) const;
-
- ClassDef(FlatExclusiveDiffXSec, 1)
-};
-
-
-
-/** Exclusive Mott Cross Section for scattering off a particle with proton mass.
- * The mott cross section describes scattering from a point-like spinless particle.
- *
- * \ingroup exclusiveXSec
- */
-class ExclusiveMottXSec : public ExclusiveDiffXSec {
- public:
- ExclusiveMottXSec() { }
-
- virtual ~ExclusiveMottXSec() { }
- virtual ExclusiveMottXSec* Clone(const char * newname) const {
- std::cout << "ExclusiveMottXSec::Clone()\n";
- auto * copy = new ExclusiveMottXSec();
- (*copy) = (*this);
- return copy;
- }
- virtual ExclusiveMottXSec* Clone() const { return( Clone("") ); }
-
- /** Needed by ROOT::Math::IBaseFunctionMultiDim */
- unsigned int NDim() const {
- return fnDim;
- }
-
- /** Returns the scattered electron energy using the angle.
- */
- Double_t GetEPrime(const Double_t theta)const {
- return (fBeamEnergy / (1.0 + 2.0 * fBeamEnergy / M_p/GeV * TMath::Power(TMath::Sin(theta / 2.0), 2)));
- }
-
- /** Evaluate Cross Section (mbarn/sr) */
- virtual Double_t EvaluateXSec(const Double_t * x) const;
-
-
- ClassDef(ExclusiveMottXSec, 1)
-};
-
-}
-}
-
-#endif
-
diff --git a/src/core/physics/include/InclusiveDiffXSec.h b/src/core/physics/include/InclusiveDiffXSec.h
deleted file mode 100644
index 8192705a140b99ce6d4074f32a320872a7d507b6..0000000000000000000000000000000000000000
--- a/src/core/physics/include/InclusiveDiffXSec.h
+++ /dev/null
@@ -1,192 +0,0 @@
-#ifndef InclusiveDiffXSec_H
-#define InclusiveDiffXSec_H 1
-#include "TROOT.h"
-#include "TObject.h"
-#include "DiffXSec.h"
-#include "TFoamIntegrand.h"
-#include "TMath.h"
-#include "StructureFunctions.h"
-#include "TCanvas.h"
-#include "TGraph.h"
-#include "TGraph2D.h"
-#include "TParticlePDG.h"
-#include "TDatabasePDG.h"
-
-namespace insane {
-namespace physics {
-
-/** Base class for an Inclusive Differential Cross Section.
- *
- * Uses units of GeV and radians. The Default unit for cross sections is
- *
- *
- * \ingroup xsections
- * \ingroup inclusiveXSec
- */
-class InclusiveDiffXSec : public DiffXSec {
-
- protected:
-
- Int_t fPolType; ///< what does this do?
- Double_t fFuncArgs[9]; // used in plotting functions
-
- public:
- InclusiveDiffXSec();
- InclusiveDiffXSec(const InclusiveDiffXSec& old);
- virtual ~InclusiveDiffXSec();
- InclusiveDiffXSec& operator=(const InclusiveDiffXSec& old);
- virtual InclusiveDiffXSec* Clone(const char * newname) const ;
- virtual InclusiveDiffXSec* Clone() const ;
-
- /** \todo how to add cross sections ? */
- InclusiveDiffXSec operator+ (InclusiveDiffXSec& right) {
- InclusiveDiffXSec temp;
- /* temp.fL*/
- return (temp);
- }
-
- virtual void DefineEvent(Double_t*);
-
- /// \todo A better way to implement polarization of cross section difference
- /// NEEDS DOCUMENTED. When and how is this used?
- virtual void SetPolarizationType(Int_t t){fPolType = t;}
-
- virtual void InitializePhaseSpaceVariables();
-
- /** TFoamIntegrand method used by TFoam.
- * Density() is sampled by TFoam with 0 < x < 1
- * \frac{d^3\sigma}{d\theta d\phi d\omega}
- */
- virtual Double_t Density(Int_t ndim, Double_t * x);
-
- /** Virtual Method used to evaluate cross section. */
- virtual Double_t EvaluateXSec(const Double_t * x) const;
-
-
-
- /** Needed by ROOT::Math::IBaseFunctionMultiDim */
- unsigned int NDim() const { return fnDim; }
-
- /** Needed by ROOT::Math::IBaseFunctionMultiDim */
- double DoEval(const double * x) const ;
-
- /** @name Useful functions for plotting.
- * For using with TF1 functions.
- * ------------------------------------------------------------
- * @{
- */
- /** Cross section as a function of W.
- * Calulates E'.
- * \param x[0] = W
- * \param p[0] = theta
- * \param p[1] = phi
- */
- virtual double Vs_W(double *x, double *p) ;
-
- /** sig(E). E=x[0], theta=p[0], phi=p[1] */
- virtual double EnergyDepXSec(double *x, double *p);
-
- /** sig(E). y=nu/Ebeam=x[0] */
- virtual double EnergyFractionDepXSec(double *x, double *p);
-
- /** Cross section as a function of energy.
- * \param x[0] = energy
- * \param p[0] = theta
- * \param p[1] = phi
- */
- virtual double EnergyDependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of energy.
- * \param x[0] = W2
- * \param p[0] = theta
- * \param p[1] = phi
- * \todo This function is incomplte.
- */
- virtual double W2DependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of energy.
- * \param x[0] = W
- * \param p[0] = Q2
- * \param p[1] = phi
- * \todo This function is incomplte.
- */
- virtual double WDependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of Bjorken x.
- * \param x[0] = x_bjorken
- * \param p[0] = Q2
- * \param p[1] = phi
- */
- virtual double xDependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of Photon energy.
- * \param x[0] = photon energy
- * \param p[0] = theta
- * \param p[1] = phi
- */
- virtual double PhotonEnergyDependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of momentum.
- * \param x[0] = momentum
- * \param p[0] = theta
- * \param p[1] = phi
- */
- virtual double MomentumDependentXSec(double *x, double *p) ;
-
- /** Cross section as a function of polar angle .
- * \param x[0] = theta
- * \param p[0] = energy
- * \param p[1] = phi
- */
- virtual double PolarAngleDependentXSec(double *x, double *p) ;
-
- virtual double PolarAngleDependentXSec_deg(double *x, double *p) ;
-
- /** Cross section as a function of polar angle .
- * \param x[0] = phi
- * \param p[0] = energy
- * \param p[1] = theta
- */
- virtual double AzimuthalAngleDependentXSec(double *x, double *p) ;
-
- virtual double Evaluate2D_p_theta(double *x, double *p){
- // using pion mass... add as parameter?...
- using namespace CLHEP;
- double pp = x[0];
- double Ep = TMath::Sqrt(x[0]*x[0]+(M_pi0/GeV)*(M_pi0/GeV));
- fFuncArgs[0] = Ep;
- fFuncArgs[1] = x[1];
- fFuncArgs[2] = p[0];
- return( (pp/Ep)*EvaluateXSec(fFuncArgs));
- }
-
- //@}
-
-
- ClassDef(InclusiveDiffXSec,2)
-};
-
-
-
-/** A flat inclusive cross section
- *
- *
- * \ingroup inclusiveXSec
- */
-class FlatInclusiveDiffXSec : public InclusiveDiffXSec {
- public:
- FlatInclusiveDiffXSec();
- virtual ~FlatInclusiveDiffXSec();
-
- /** Evaluate Cross Section. Flat cross section returns 1 */
- virtual Double_t EvaluateXSec(const Double_t * x) const;
-
- ClassDef(FlatInclusiveDiffXSec, 1)
-};
-
-}
-}
-
-#endif
-
-
diff --git a/src/core/physics/include/LinkDef.h b/src/core/physics/include/LinkDef.h
index 46cfc4e6ab13e6fd915d678977883db36dc5852c..34f14b8506b24c6c9e986e812c3fe6a35ef335a3 100644
--- a/src/core/physics/include/LinkDef.h
+++ b/src/core/physics/include/LinkDef.h
@@ -96,7 +96,6 @@
//
#pragma link C++ class insane::physics::_EG_Event+;
-#pragma link C++ class insane::physics::EventGenerator+;
#pragma link C++ class insane::physics::TargetEventGenerator+;
#pragma link C++ class insane::physics::BETAG4SavedEventGenerator+;
@@ -105,12 +104,12 @@
//#pragma link C++ class UVACarbonTarget+;
//#pragma link C++ class UVACrossHairTarget+;
//#pragma link C++ class UVAPureHeliumTarget+;
-#pragma link C++ class insane::physics::Helium4GasTarget+;
-#pragma link C++ class insane::physics::TDISDeuteriumGasTarget+;
-#pragma link C++ class insane::physics::DeuteriumGasTarget+;
-#pragma link C++ class insane::physics::HydrogenGasTarget+;
-#pragma link C++ class insane::physics::TDISHydrogenGasTarget+;
-#pragma link C++ class insane::physics::LH2Target+;
+//#pragma link C++ class insane::physics::Helium4GasTarget+;
+//#pragma link C++ class insane::physics::TDISDeuteriumGasTarget+;
+//#pragma link C++ class insane::physics::DeuteriumGasTarget+;
+//#pragma link C++ class insane::physics::HydrogenGasTarget+;
+//#pragma link C++ class insane::physics::TDISHydrogenGasTarget+;
+//#pragma link C++ class insane::physics::LH2Target+;
#pragma link C++ class std::vector<TString*>+;
#pragma link C++ class std::map<int, double >+;
@@ -228,20 +227,20 @@
#pragma link C++ class insane::physics::LSS2010PolarizedPDFs+;
#pragma link C++ class insane::physics::DSSVPolarizedPDFs+;
-#pragma link C++ class insane::physics::PhaseSpace+;
-#pragma link C++ class insane::physics::PhaseSpaceVariable+;
-#pragma link C++ class insane::physics::DiscretePhaseSpaceVariable+;
-#pragma link C++ class std::vector<insane::physics::PhaseSpaceVariable*>+;
+//#pragma link C++ class insane::physics::PhaseSpace+;
+//#pragma link C++ class insane::physics::PhaseSpaceVariable+;
+//#pragma link C++ class insane::physics::DiscretePhaseSpaceVariable+;
+//#pragma link C++ class std::vector<insane::physics::PhaseSpaceVariable*>+;
-#pragma link C++ class insane::physics::DiffXSec+;
+//#pragma link C++ class insane::physics::DiffXSec+;
#pragma link C++ class insane::physics::CompositeDiffXSec+;
#pragma link C++ class insane::physics::DiffXSecKinematicKey+;
#pragma link C++ class insane::physics::GridDiffXSec+;
#pragma link C++ class insane::physics::GridXSecValue+;
-#pragma link C++ class insane::physics::InclusiveDiffXSec+;
+//#pragma link C++ class insane::physics::InclusiveDiffXSec+;
#pragma link C++ class insane::physics::InclusiveDISXSec+;
#pragma link C++ class insane::physics::InclusiveBornDISXSec+;
-#pragma link C++ class insane::physics::ExclusiveDiffXSec+;
+//#pragma link C++ class insane::physics::ExclusiveDiffXSec+;
#pragma link C++ class insane::physics::MAIDInclusiveDiffXSec+;
#pragma link C++ class insane::physics::MAIDNucleusInclusiveDiffXSec+;
diff --git a/src/core/physics/include/PhaseSpace.h b/src/core/physics/include/PhaseSpace.h
deleted file mode 100644
index 69d677983f5144fa93069ded8e43adbeda564e94..0000000000000000000000000000000000000000
--- a/src/core/physics/include/PhaseSpace.h
+++ /dev/null
@@ -1,241 +0,0 @@
-#ifndef PhaseSpace_H
-#define PhaseSpace_H 1
-
-#include "PhaseSpaceVariable.h"
-
-#include <vector>
-#include <string>
-#include <iostream>
-#include <iomanip>
-
-#include "TMath.h"
-#include "TObject.h"
-#include "TNamed.h"
-#include "TString.h"
-#include "TFoamIntegrand.h"
-#include "TList.h"
-
-namespace insane {
-namespace physics {
-
-/** Base for phase space of a differential cross section.
- *
- * \f$ \frac{d^N\sigma}{dx_1 dx_2 ... dx_N} \f$
- *
- * \ingroup xsections
- * \ingroup EventGen
- */
-class PhaseSpace : public TObject {
-
- protected:
- TList fVariables;
- Int_t fnDim;
- Int_t fnIndependentVars;
- Bool_t fIsModified;
- Double_t fScaledIntegralJacobian; // For multiplying the foam integral result to get the actual integral
- // this is because foam uses variables between 0 and 1
-
- public:
- PhaseSpace() {
- fVariables.SetOwner(true);
- fVariables.Clear();
- fnDim = 0;
- fnIndependentVars = 0;
- fIsModified = false;
- fScaledIntegralJacobian = 1.0;
- }
- PhaseSpace(const PhaseSpace& rhs) :
- TObject(rhs), fnDim(rhs.fnDim), fnIndependentVars(rhs.fnIndependentVars), fIsModified(rhs.fIsModified)
- {
- fVariables.AddAll(&(rhs.fVariables));
- }
-
- virtual ~PhaseSpace() {
- // Delete the PS variables. The TList fVariables is the owner so they are deleted when cleared.
- fVariables.Clear();
- }
-
- //PhaseSpace(const PhaseSpace&) = default; // Copy constructor
- PhaseSpace(PhaseSpace&&) = default; // Move constructor
- PhaseSpace& operator=(const PhaseSpace&) & = default; // Copy assignment operator
- PhaseSpace& operator=(PhaseSpace&&) & = default; // Move assignment operato
- //virtual ~PhaseSpace() { }
-
- void ClearVariables() {
- fVariables.Clear();
- SetModified(true);
- }
-
- Double_t GetScaledIntegralJacobian() const { return fScaledIntegralJacobian; }
-
- Int_t GetNIndependentVariables() const { return fnIndependentVars; }
-
- Int_t GetDimension() const {
- return(fnDim) ;
- }
- void SetDimension(Int_t dim) {
- fnDim = dim;
- SetModified(true);
- }
-
- void AddVariable(PhaseSpaceVariable * var) {
- fVariables.Add(var);
- fnDim++;
- if(!var->IsDependent()) fnIndependentVars++;
- SetModified(true);
- }
-
-
- void Print(Option_t * opt = "") const ;
- void Print(std::ostream& stream) const ;
-
- /** Returns a list of PhaseSpaceVariable pointers.
- * Use pIndex to select a particle's variables.
- * By default it returns a list with all phase space variables
- */
- TList * GetVariables(Int_t pIndex = -1) {
- TList * thelist = nullptr;
- const int N = fVariables.GetEntries();
- if (pIndex == -1) thelist = &fVariables;
- else {
- thelist = new TList();
- for (int i = 0; i < N; i++) {
- if (GetVariable(i)->GetParticleIndex() == pIndex) thelist->Add(GetVariable(i));
- }
- }
- return(thelist);
- }
-
- /** get the ith variable added */
- PhaseSpaceVariable* GetVariable(Int_t i) const {
- if (i < fnDim) return((PhaseSpaceVariable*)fVariables.At(i));
- else {
- std::cout << " Index " << i << " is larger than phase space dimension, " << fnDim << ". \n";
- return(nullptr);
- }
- }
-
- /** Get the PhaseSpaceVariable by name
- * \note You should use theta,phi,energy,momentum,... etc when naming Phase space variables.
- * Usually there is one particle so if it has a PS variable named "theta_e", you can
- * just ask for "theta" and it will return the first variable to have "theta" in it.
- * If there are two thetas then it returns the first one,
- * and therefore you should be more explicit in your search.
- *
- */
- PhaseSpaceVariable* GetVariableWithName(const char * name) {
- const int N = fVariables.GetEntries();
- for (int i = 0; i < N; i++) {
- TString aVarName = GetVariable(i)->GetName();
- if (aVarName.Contains(name)) return(GetVariable(i));
- /* if( !strcmp(name,GetVariable(i)->GetName()) ) return(GetVariable(i));*/
- }
- std::cout << " PhaseSpaceVariable named " << name << " was not found!\n";
- return(nullptr);
- }
-
- PhaseSpaceVariable* GetIndependentVariable(int i_ind) {
- const int N = fVariables.GetEntries();
- int n_ind = 0;
- for (int i = 0; i < N; i++) {
- auto var = GetVariable(i);
- if( !(var->IsDependent()) ) {
- if( n_ind == i_ind ){
- return var;
- }
- n_ind++;
- }
- }
- return(nullptr);
- }
-
- /** Called by DiffXSec::Refresh() */
- void Refresh() {
- fScaledIntegralJacobian = 1.0;
- const int N = fVariables.GetEntries();
- for (int i = 0; i < N; i++) {
- auto * aVar = (PhaseSpaceVariable*) GetVariable(i);
- fScaledIntegralJacobian *= aVar->GetNormScale();
- aVar->SetCurrentValue( aVar->GetCentralValue() );
- aVar->SetModified(false);
- }
- }
-
- Double_t GetSolidAngle(Int_t ipart = 0){
- PhaseSpaceVariable * phi = GetVariableWithName("phi");
- PhaseSpaceVariable * theta = GetVariableWithName("theta");
- if(!phi) return 0.0;
- if(!theta) return 0.0;
- Double_t deltaPhi = phi->GetMaximum() - phi->GetMinimum();
- Double_t deltaTheta = -1.0*(TMath::Cos(theta->GetMaximum()) - TMath::Cos(theta->GetMinimum()));
- return deltaPhi*deltaTheta;
- }
- Double_t GetEnergyBite(Int_t ipart = 0){
- PhaseSpaceVariable * energy = GetVariableWithName("energy");
- Double_t deltaE = energy->GetMaximum() - energy->GetMinimum();
- return deltaE;
- }
-
- void ListVariables() {
- std::cout << " Variables in phase space are:\n";
- const int N = fVariables.GetEntries();
- for (int i = 0; i < N; i++) {
- std::cout << GetVariable(i)->GetName() << "\n";
- }
- std::cout << "\n";
- }
-
- /** Sets the values of each PS variable. Is called from DiffXSec::VariablesInPhaseSpace().
- * The dimension vals should be of dimension fNParticles*3
- */
- void SetEventValues(const Double_t * vals, int nVars=0) const {
- int N = fVariables.GetEntries();
- if(nVars > 0) {
- N= nVars;
- }
- for (int i = 0; i < N; i++) {
- GetVariable(i)->SetCurrentValue(vals[i]);
- }
- }
-
- /** Returns the array filled with the current values of each phase space
- * variable.
- */
- void GetEventValues(Double_t * vals) {
- const int N = fVariables.GetEntries();
- if(N>0){
- for (int i=0;i<N;i++) {
- vals[i] = GetVariable(i)->GetCurrentValue();
- }
- }else{
- std::cout << "[PhaseSpace::GetEventValues]: WARNING! ";
- std::cout << "There seems to be no phase space variables..." << std::endl;
- }
-
- }
-
- void PrintTableVariables() {
- const int N = fVariables.GetEntries();
- for (int i = 0; i < N; i++) {
- GetVariable(i)->PrintTable();
- }
- }
-
- bool IsModified() const {
- return fIsModified;
- }
- void SetModified(bool val = true) {
- fIsModified = val;
- }
-
-
-
- ClassDef(PhaseSpace,3)
-};
-
-}
-}
-
-#endif
-
-
diff --git a/src/core/physics/include/PhaseSpaceSampler.h b/src/core/physics/include/PhaseSpaceSampler.h
deleted file mode 100644
index 3217e09fa70d0585c8173ba92740c485bfab44ca..0000000000000000000000000000000000000000
--- a/src/core/physics/include/PhaseSpaceSampler.h
+++ /dev/null
@@ -1,152 +0,0 @@
-#ifndef PhaseSpaceSampler_H
-#define PhaseSpaceSampler_H 1
-#include <TROOT.h>
-#include <TObject.h>
-#include "TFoam.h"
-#include "TFoamIntegrand.h"
-#include "TMath.h"
-#include "DiffXSec.h"
-#include "TRandom3.h"
-#include "TList.h"
-#include <iostream>
-#include <iomanip>
-
-namespace insane {
-namespace physics {
-
-/** Base class for a random sampling of a phase space.
- * This class uses ROOT's TFoam classes for sampling the phase space.
- *
- * The phase-space variables are scaled to be between 0 and 1 in order to work with TFoam.
- * This is why the result of foam is pass through DiffXSec::GetUnnormalizedVariables();
- * The final result of PhaseSpaceSampler::GenerateEvent() comes from passing the random
- * (unnormalized variables) through the virutal method DiffXSec::GetDependentVariables(),
- * which appends any possible variables which are kinematically constrained (and thus not random).
- * For example, in exclusive elastic scattering where the electron variables are random, the
- * nuclei's variables are calculated and appended to the result.
- *
- * Note: The sampler owns the cross section.
- *
- * \ingroup EventGen
- */
-class PhaseSpaceSampler : public TObject {
- protected:
- Double_t fWeight; /// Weight used by event generator
- Bool_t fIsModified;
- Double_t fTotalXSection;
- Double_t fMCvect[30]; //
- DiffXSec * fDiffXSec; //!
- TFoam * fFoam; //->
- TRandom * fRandomNumberGenerator; //!
- TList fXSectionList;
- Int_t fFoamCells; // number of cells
- Int_t fFoamSample; // number of MC events per cell in build-up
- Int_t fFoamBins; // number of bins in build-up
- Int_t fFoamOptRej ; // Wted events for OptRej=0; wt=1 for OptRej=1 (default)
- Int_t fFoamOptDrive; // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax
- Int_t fFoamEvPerBin; // Maximum events (equiv.) per bin in buid-up
- Int_t fFoamChat ; // Chat level };
- Double_t fFoamMaxWtRej; // 1.0
-
- TList fCrossSection; // This exists to store the cross section which doesn't stream by itself
- // but does when it is in a container. When the sampler is read back from
- // a file it should initialized with this cross section by calling InitFromDisk
-
- public:
- PhaseSpaceSampler(DiffXSec * xsec = nullptr);
-
- PhaseSpaceSampler(const PhaseSpaceSampler& rhs);
- PhaseSpaceSampler& operator=(const PhaseSpaceSampler& rhs);
-
- //PhaseSpaceSampler(const PhaseSpaceSampler&) = default; // Copy constructor
- PhaseSpaceSampler(PhaseSpaceSampler&&) = default; // Move constructor
- //PhaseSpaceSampler& operator=(const PhaseSpaceSampler&) & = default; // Copy assignment operator
- PhaseSpaceSampler& operator=(PhaseSpaceSampler&&) & = default; // Move assignment operato
- //virtual ~PhaseSpaceSampler() { }
-
- /** The sampler owns the XSection so it is deleted. */
- virtual ~PhaseSpaceSampler();
-
- void InitFromDisk() ;
-
- TFoam * GetFoam(){ return(fFoam); }
-
- void SetFoamCells(Int_t v) { fFoamCells = v; if(fFoam) fFoam->SetnCells(fFoamCells); }
- void SetFoamSample(Int_t v) { fFoamSample = v; if(fFoam) fFoam->SetnSampl(fFoamSample); }
- void SetFoamBins(Int_t v) { fFoamBins = v; if(fFoam) fFoam->SetnBin(fFoamBins); }
- void SetFoamOptRej(Int_t v) { fFoamOptRej = v; if(fFoam) fFoam->SetOptRej(fFoamOptRej); }
- void SetFoamOptDrive(Int_t v) { fFoamOptDrive = v; if(fFoam) fFoam->SetOptDrive(fFoamOptDrive); }
- void SetFoamEvPerBin(Int_t v) { fFoamEvPerBin = v; if(fFoam) fFoam->SetEvPerBin(fFoamEvPerBin); }
- void SetFoamChat(Int_t v) { fFoamChat = v; if(fFoam) fFoam->SetChat(fFoamChat); }
- void SetFoamMaxWtRej(Double_t v) { fFoamMaxWtRej = v; if(fFoam) fFoam->SetMaxWtRej(fFoamMaxWtRej); }
-
- void PrintFoamConfig(std::ostream& c = std::cout) const;
-
- /** Refresh creates a new TFoam object after changing the parameters
- * of the phase space or cross section.
- * If a cross section argument is provided the old one is deleted and the new one is used
- * (creating a new phase space space too).
- */
- void Refresh(DiffXSec * xsec = nullptr) ;
-
- /** Generates a new event and returns an array of random variables.
- * Returns array fiilled with variables.
- */
- Double_t * GenerateEvent();
-
- /** Calculates the cross section normalization (ie total cross section).
- *
- * The cross section normalization is just the total cross section
- * which is just the integral of the differential cross section over the
- * entire phase space.
- */
- virtual Double_t NormalizePDF();
-
- /** Set Cross Section whos weighted phase-space is used to generate random variables.
- */
- void SetXSec(DiffXSec * xsec, Bool_t mod=true) ;
- DiffXSec * GetXSec() const { return(fDiffXSec); }
-
- void Print(const Option_t* option = "") const ;
- void Print(std::ostream& stream) const ;
-
- Double_t GetTotalXSection() const { return(fTotalXSection); }
- Double_t CalculateTotalXSection();
-
- bool IsModified() const { return fIsModified; }
- void SetModified(bool val = false) { fIsModified = val; }
-
- /** Set the weight used by the event generator. Must be between 0 and 1.0 */
- void SetWeight(Double_t w) {
- if( w == 0.0 ) Warning("SetWeight","Weight set to zero!");
- else if(w < 0.0 ) Error("SetWeight","Weight less than zero!");
- else if(w > 1.0 ) Error("SetWeight","Weight bigger than 1");
- else fWeight = w;
- }
- Double_t GetWeight(){return fWeight;}
-
-
-
- /** Work around to load the first cross section saved in list when the
- * cross section is saved to a file. F
- * INCOMPLETE
- */
- void Load() {
- if (fXSectionList.GetEntries() == 0) {
- std::cout << " NO SAVED CROSS SECTIONS!!!!\n";
- } else {
- SetXSec((DiffXSec*)fXSectionList.At(0));
- /* GetXSec()->InitializePhaseSpaceVariables(); */
- // fMomentum_pi = varEnergy2->GetCurrentValueAddress();
-
- }
- }
-
- ClassDef(PhaseSpaceSampler, 3)
-};
-
-}
-}
-
-#endif
-
diff --git a/src/core/physics/include/PhaseSpaceVariable.h b/src/core/physics/include/PhaseSpaceVariable.h
deleted file mode 100644
index 7a1dad7fb6849a7ba91de3871b767ade09aa754c..0000000000000000000000000000000000000000
--- a/src/core/physics/include/PhaseSpaceVariable.h
+++ /dev/null
@@ -1,203 +0,0 @@
-#ifndef PhaseSpaceVariable_HH
-#define PhaseSpaceVariable_HH 1
-
-#include <iostream>
-#include <iomanip>
-#include <vector>
-#include <string>
-
-#include "TMath.h"
-#include "TNamed.h"
-#include "TString.h"
-#include "TList.h"
-
-namespace insane {
-namespace physics {
-
-/** Base class for a phase space variable
- *
- * The variable normalization is for FOAM. Each variable is scaled using the
- * minima and maxima such that its scaled value falls between 0 and 1 (as required
- * for FOAM random variables)
- *
- * @param name The name given to refer to the variable.
- * @param varexp The displayed expression for the varaible.
- * This can include the ROOT style latex, eg #theta
- *
- * \ingroup EventGen
- * \ingroup xsections
- */
-class PhaseSpaceVariable : public TNamed {
- protected:
- Bool_t fIsModified = true;
- Bool_t fIsDependent = false;
- Bool_t fInverted = false;
- Bool_t fUniform = false;
- Int_t fiParticle = -1; // associated particle number (for events with multiple particles).
- Double_t fVariableMinima = 0.0;
- Double_t fVariableMaxima = 1.0;
- TString fVariableUnits = "";
- mutable Double_t fCurrentValue = 0.1;
-
- public:
- PhaseSpaceVariable(const char * name = "" , const char * varexp = "" , Double_t min = 0.0, Double_t max = 1.0);
- PhaseSpaceVariable(const PhaseSpaceVariable& rhs);
-
- //PhaseSpaceVariable(const PhaseSpaceVariable&) = default; // Copy constructor
- PhaseSpaceVariable(PhaseSpaceVariable&&) = default; // Move constructor
- PhaseSpaceVariable& operator=(const PhaseSpaceVariable&) & = default; // Copy assignment operator
- PhaseSpaceVariable& operator=(PhaseSpaceVariable&&) & = default; // Move assignment operato
- virtual ~PhaseSpaceVariable();
-
- Double_t GetCurrentValue(){ return fCurrentValue; }
- void SetCurrentValue(Double_t v){ fCurrentValue = v;}
-
- Double_t * GetCurrentValueAddress() { return &fCurrentValue ; }
-
- void SetParticleIndex(Int_t i) { fiParticle = i; }
- Int_t GetParticleIndex() const { return fiParticle; }
-
- Double_t GetCentralValue() const { return((fVariableMinima + fVariableMaxima) / 2.0); }
-
- /** Use for azimuthal angles near the branch cut at -pi/pi */
- void SetInverted(Bool_t inv = true) { fInverted = inv; fIsModified = true; }
-
- void SetUniform(Bool_t u) { fUniform = u; fIsModified = true; }
- Bool_t IsUniform() const { return fUniform; }
-
- void SetDependent(Bool_t val = true) { fIsDependent = val; fIsModified = true; }
- Bool_t IsDependent() { return fIsDependent; }
-
- Bool_t IsModified() const { return fIsModified; }
- void SetModified(bool val = false) { fIsModified = val; }
-
- /** Note that we invert the cut for azimuthal angles centered around the
- * branch cut near phi ~ -pi or pi .
- *
- * Here we use 0 < phi < 2pi
- */
- Bool_t IsInVariableRange() const {
- //Double_t min = GetMinimum();
- //Double_t max = GetMaximum();
- //if( !( (fCurrentValue >= min)&&(fCurrentValue <= max) ) ){
- // std::cout << "[PhaseSpace::IsInVariableRange]: Current value = " << fCurrentValue << std::endl;
- // std::cout << "min: " << min << std::endl;
- // std::cout << "max: " << max << std::endl;
- //}
- if (!fInverted) {
- return((fCurrentValue <= GetMaximum()) &&
- (fCurrentValue >= GetMinimum()));
- } /* else it is inverted (e.g. azimuthal angle near -pi/pi cut!) */
- if (fCurrentValue > 0.0) {
- return((fCurrentValue <= GetMaximum()) &&
- (fCurrentValue >= GetMinimum()));
- }
- if (fCurrentValue < 0.0) {
- return((fCurrentValue + 2.0 * TMath::Pi() <= GetMaximum()) &&
- (fCurrentValue + 2.0 * TMath::Pi() >= GetMinimum())) ;
- }
- return false;
- }
-
- /** Returns offset needed such that 0 < (variable-offset)/scale < 1
- */
- Double_t GetNormScale() const {
- return(fVariableMaxima - GetNormOffset());
- }
-
- /** Returns offset needed such that 0 < (variable-offset)/scale < 1
- */
- Double_t GetNormOffset() const {
- return(fVariableMinima);
- }
-
- void SetRange(double min, double max){
- if(min>max) Warning("SetRange","min is greater than max");
- SetMinimum(min);
- SetMaximum(max);
- }
-
- void SetMaximum(double val) {
- fVariableMaxima = val;
- fIsModified = true;
- }
- void SetMinimum(double val) {
- fVariableMinima = val;
- fIsModified = true;
- }
- double GetMaximum() const {
- return(fVariableMaxima);
- }
- double GetMinimum() const {
- return(fVariableMinima);
- }
- void SetVariableUnits(const char * unit) {
- fVariableUnits = unit;
- fIsModified = true;
- }
-
- const char * GetVariableUnits() const {
- return(fVariableUnits.Data());
- }
-
- void Print(const Option_t * opt = "") const ;
-
- void Print(std::ostream& stream) const ;
-
- void PrintTable() const {
- std::cout << " " << GetName() << " = " << fCurrentValue << "\n";
- }
-
-
- ClassDef(PhaseSpaceVariable, 3)
-};
-
-
-/** A discrete phase space variable such as helicity.
- *
- * \ingroup EventGen
- * \ingroup xsections
- */
-class DiscretePhaseSpaceVariable : public PhaseSpaceVariable {
- private:
- Int_t fNDivisions;
-
- public:
- DiscretePhaseSpaceVariable(const char * n = "", const char * t = "") : PhaseSpaceVariable(n,t) {
- fNDivisions = 2;
- }
- DiscretePhaseSpaceVariable(const DiscretePhaseSpaceVariable& rhs) : PhaseSpaceVariable(rhs),
- fNDivisions(rhs.fNDivisions) {}
-
- //PhaseSpaceVariable(const PhaseSpaceVariable&) = default; // Copy constructor
- DiscretePhaseSpaceVariable(DiscretePhaseSpaceVariable&&) = default; // Move constructor
- DiscretePhaseSpaceVariable& operator=(const DiscretePhaseSpaceVariable&) & = default; // Copy assignment operator
- DiscretePhaseSpaceVariable& operator=(DiscretePhaseSpaceVariable&&) & = default; // Move assignment operato
- virtual ~DiscretePhaseSpaceVariable() { }
-
- void SetNumberOfValues(int i) {
- fNDivisions = i;
- fIsModified = true;
- }
- Int_t GetNumberOfValues() const {
- return(fNDivisions);
- }
-
- /** Returns an integer value less than the set number of possible values
- * where the range 0 to 1 is divided equally.
- *
- * @param randomValue should be a value between 0 and 1
- */
- Int_t GetDiscreteVariable(Double_t randomValue) const {
- Double_t val = randomValue * ((Double_t)fNDivisions);
- return((Int_t)val);
- }
-
- ClassDef(DiscretePhaseSpaceVariable, 1)
-};
-
-}
-}
-
-#endif
-
diff --git a/src/core/physics/src/CLASTargets.cxx b/src/core/physics/src/CLASTargets.cxx
deleted file mode 100644
index a9f941a4f05eb87c02f857a119af3eedafd90371..0000000000000000000000000000000000000000
--- a/src/core/physics/src/CLASTargets.cxx
+++ /dev/null
@@ -1,294 +0,0 @@
-#include "CLASTargets.h"
-
-namespace insane {
- namespace physics {
-
- Helium4GasTarget::Helium4GasTarget(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- Helium4GasTarget::~Helium4GasTarget()
- { }
- //______________________________________________________________________________
-
- void Helium4GasTarget::DefineMaterials()
- {
- using namespace CLHEP;
- double target_length = 35.0;
- auto * matHe4 = new TargetMaterial("He4","He4",2,4);
-
- double density_at_1atm_300K = 0.164*kg/m3; // kg/m3
- double a_noUnit = 4.003;
- double pre_noUnit = 3.0;
- double tempe_noUnit = 298.;
- double pressure = 3.0*atmosphere;
- double temperature = 298.*kelvin;
- double He4_density = density_at_1atm_300K*pre_noUnit;
- // density is 0.178 at 0C, 0.169 at 15C, 0.165 at 25C (all in units of kg/m3)
- //
- //(a_noUnit*pre_noUnit)/(0.0821*tempe_noUnit)*kg/m3; //0.164*kg/m3 at 1 atm;
-
- matHe4->fLength = target_length; //cm
- matHe4->fDensity = He4_density/(g/cm3); // g/cm3
- matHe4->fZposition = 0.0; //cm
-
- double window_thickness = 15.0/10000.0 ; //15 um
- auto * matAlEndcap1 = new TargetMaterial("Al-endcap1", "Aluminum target endcap upstream", 13, 27);
- matAlEndcap1->fLength = window_thickness ; //15 um
- matAlEndcap1->fDensity = 2.7; // g/cm3
- matAlEndcap1->fZposition = -(target_length/2.0+window_thickness/2.0); //cm
-
- auto * matAlEndcap2 = new TargetMaterial("Al-endcap2", "Aluminum target endcap downstream", 13, 27);
- matAlEndcap2->fLength = window_thickness ; //cm
- matAlEndcap2->fDensity = 2.7; // g/cm3
- matAlEndcap2->fZposition = target_length/2.0; //cm
- matAlEndcap2->fZposition = (target_length/2.0+window_thickness/2.0); //cm
-
- this->AddMaterial(matHe4);
- this->AddMaterial(matAlEndcap1);
- this->AddMaterial(matAlEndcap2);
- }
- //______________________________________________________________________________
-
-
-
- DeuteriumGasTarget::DeuteriumGasTarget(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- DeuteriumGasTarget::~DeuteriumGasTarget()
- { }
- //______________________________________________________________________________
-
- void DeuteriumGasTarget::DefineMaterials()
- {
- using namespace CLHEP;
- double target_length = 35.0;
- auto * matD = new TargetMaterial("D","D",1,2);
-
- double density_at_1atm_300K = 0.165*kg/m3;
- // density is 0.179 at 0C, to 0.170 at 15C, to 0.165 at 25C.
- double a_noUnit = 2.01410178;
- double pre_noUnit = 3.0;
- double tempe_noUnit = 298.;
- double pressure = 3.0*atmosphere;
- double temperature = 298.*kelvin;
- double He4_density = density_at_1atm_300K*pre_noUnit;//(a_noUnit*pre_noUnit)/(0.0821*tempe_noUnit)*kg/m3; //0.164*kg/m3 at 1 atm;
-
- matD->fLength = target_length; //cm
- matD->fDensity = He4_density/(g/cm3); // g/cm3
- matD->fZposition = 0.0; //cm
-
- double window_thickness = 15.0/10000.0 ; //15 um
- auto * matAlEndcap1 = new TargetMaterial("Al-endcap1", "Aluminum target endcap upstream", 13, 27);
- matAlEndcap1->fLength = window_thickness ; //15 um
- matAlEndcap1->fDensity = 2.7; // g/cm3
- matAlEndcap1->fZposition = -(target_length/2.0+window_thickness/2.0); //cm
-
- auto * matAlEndcap2 = new TargetMaterial("Al-endcap2", "Aluminum target endcap downstream", 13, 27);
- matAlEndcap2->fLength = window_thickness ; //cm
- matAlEndcap2->fDensity = 2.7; // g/cm3
- matAlEndcap2->fZposition = target_length/2.0; //cm
- matAlEndcap2->fZposition = (target_length/2.0+window_thickness/2.0); //cm
-
- this->AddMaterial(matD);
- this->AddMaterial(matAlEndcap1);
- this->AddMaterial(matAlEndcap2);
- }
- //______________________________________________________________________________
-
-
-
- HydrogenGasTarget::HydrogenGasTarget(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- HydrogenGasTarget::~HydrogenGasTarget()
- { }
- //______________________________________________________________________________
-
- void HydrogenGasTarget::DefineMaterials()
- {
- using namespace CLHEP;
- double target_length = 35.0;
- auto * matD = new TargetMaterial("H","H",1,1);
-
- double density_at_1atm_300K = 0.0823*kg/m3;
- // density is 0.0899 at 0C, to 0.0852 at 15C, to 0.0823 at 25C.
- double a_noUnit = 2.01410178;
- double pre_noUnit = 3.0;
- double tempe_noUnit = 298.;
- double pressure = 3.0*atmosphere;
- double temperature = 298.*kelvin;
- double He4_density = density_at_1atm_300K*pre_noUnit;//(a_noUnit*pre_noUnit)/(0.0821*tempe_noUnit)*kg/m3; //0.164*kg/m3 at 1 atm;
-
- matD->fLength = target_length; //cm
- matD->fDensity = He4_density/(g/cm3); // g/cm3
- matD->fZposition = 0.0; //cm
-
- double window_thickness = 15.0/10000.0 ; //15 um
- auto * matAlEndcap1 = new TargetMaterial("Al-endcap1", "Aluminum target endcap upstream", 13, 27);
- matAlEndcap1->fLength = window_thickness ; //15 um
- matAlEndcap1->fDensity = 2.7; // g/cm3
- matAlEndcap1->fZposition = -(target_length/2.0+window_thickness/2.0); //cm
-
- auto * matAlEndcap2 = new TargetMaterial("Al-endcap2", "Aluminum target endcap downstream", 13, 27);
- matAlEndcap2->fLength = window_thickness ; //cm
- matAlEndcap2->fDensity = 2.7; // g/cm3
- matAlEndcap2->fZposition = target_length/2.0; //cm
- matAlEndcap2->fZposition = (target_length/2.0+window_thickness/2.0); //cm
-
- this->AddMaterial(matD);
- this->AddMaterial(matAlEndcap1);
- this->AddMaterial(matAlEndcap2);
- }
- //______________________________________________________________________________
-
-
- LH2Target::LH2Target(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- LH2Target::~LH2Target()
- { }
- //______________________________________________________________________________
-
- void LH2Target::DefineMaterials()
- {
- Double_t LH2_density = 0.07085; //g/cm^3
- Double_t target_length = 5.0;
- auto * matLH2 = new TargetMaterial("LH2","LH2",1,1);
- matLH2->fLength = target_length; //cm
- matLH2->fDensity = LH2_density; // g/cm3
- matLH2->fZposition = 0.0; //cm
-
- auto * matAlEndcap1 = new TargetMaterial("Al-endcap1", "Aluminum target endcap upstream", 13, 27);
- matAlEndcap1->fLength = 0.00381 ; //cm
- matAlEndcap1->fDensity = 2.7; // g/cm3
- matAlEndcap1->fZposition = -target_length/2.0; //cm
-
- auto * matAlEndcap2 = new TargetMaterial("Al-endcap2", "Aluminum target endcap downstream", 13, 27);
- matAlEndcap2->fLength = 0.00381 ; //cm
- matAlEndcap2->fDensity = 2.7; // g/cm3
- matAlEndcap2->fZposition = target_length/2.0; //cm
-
- this->AddMaterial(matLH2);
- this->AddMaterial(matAlEndcap1);
- this->AddMaterial(matAlEndcap2);
-
- }
- //______________________________________________________________________________
-
- TDISDeuteriumGasTarget::TDISDeuteriumGasTarget(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- TDISDeuteriumGasTarget::~TDISDeuteriumGasTarget()
- { }
- //______________________________________________________________________________
-
- void TDISDeuteriumGasTarget::DefineMaterials()
- {
- using namespace CLHEP;
- double target_length = 40.0;
- auto * matD = new TargetMaterial("D","D",1,2);
-
- double density_at_1atm_300K = 0.165*kg/m3;
- // density is 0.179 at 0C, to 0.170 at 15C, to 0.165 at 25C.
- double a_noUnit = 2.01410178;
- double pre_noUnit = 1.0;
- double tempe_noUnit = 298.;
- double pressure = 1.0*atmosphere;
- double temperature = 298.*kelvin;
- double He4_density = density_at_1atm_300K*pre_noUnit;//(a_noUnit*pre_noUnit)/(0.0821*tempe_noUnit)*kg/m3; //0.164*kg/m3 at 1 atm;
-
- matD->fLength = target_length; //cm
- matD->fDensity = He4_density/(g/cm3); // g/cm3
- matD->fZposition = 0.0; //cm
-
- double window_thickness = 20.0/10000.0 ; //20 um
- auto * matBeEndcap1 = new TargetMaterial("Be-endcap1", "Beryllium target endcap upstream", 4, 9);
- matBeEndcap1->fLength = window_thickness ; //15 um
- matBeEndcap1->fDensity = 1.85; // g/cm3
- matBeEndcap1->fZposition = -(target_length/2.0+window_thickness/2.0); //cm
-
- auto * matBeEndcap2 = new TargetMaterial("Be-endcap2", "Beryllium target endcap downstream", 4, 9);
- matBeEndcap2->fLength = window_thickness ; //cm
- matBeEndcap2->fDensity = 1.85; // g/cm3
- matBeEndcap2->fZposition = target_length/2.0; //cm
- matBeEndcap2->fZposition = (target_length/2.0+window_thickness/2.0); //cm
-
- this->AddMaterial(matD);
- this->AddMaterial(matBeEndcap1);
- this->AddMaterial(matBeEndcap2);
- }
- //______________________________________________________________________________
-
- TDISHydrogenGasTarget::TDISHydrogenGasTarget(const char * name, const char * title) :
- Target(name,title)
- {
- DefineMaterials();
- }
- //______________________________________________________________________________
-
- TDISHydrogenGasTarget::~TDISHydrogenGasTarget()
- { }
- //______________________________________________________________________________
-
- void TDISHydrogenGasTarget::DefineMaterials()
- {
- using namespace CLHEP;
- double target_length = 40.0;
- auto * matD = new TargetMaterial("H","H",1,1);
-
- double density_at_1atm_300K = 0.0823*kg/m3;
- // density is 0.0899 at 0C, to 0.0852 at 15C, to 0.0823 at 25C.
- double a_noUnit = 2.01410178;
- double pre_noUnit = 1.0;
- double tempe_noUnit = 298.;
- double pressure = 1.0*atmosphere;
- double temperature = 298.*kelvin;
- double He4_density = density_at_1atm_300K*pre_noUnit;//(a_noUnit*pre_noUnit)/(0.0821*tempe_noUnit)*kg/m3; //0.164*kg/m3 at 1 atm;
-
- matD->fLength = target_length; //cm
- matD->fDensity = He4_density/(g/cm3); // g/cm3
- matD->fZposition = 0.0; //cm
-
- double window_thickness = 20.0/10000.0 ; //20 um
- auto * matBeEndcap1 = new TargetMaterial("Be-endcap1", "Beryllium target endcap upstream", 4, 9);
- matBeEndcap1->fLength = window_thickness ; //15 um
- matBeEndcap1->fDensity = 1.85; // g/cm3
- matBeEndcap1->fZposition = -(target_length/2.0+window_thickness/2.0); //cm
-
- auto * matBeEndcap2 = new TargetMaterial("Be-endcap2", "Beryllium target endcap downstream", 4, 9);
- matBeEndcap2->fLength = window_thickness ; //cm
- matBeEndcap2->fDensity = 1.85; // g/cm3
- matBeEndcap2->fZposition = target_length/2.0; //cm
- matBeEndcap2->fZposition = (target_length/2.0+window_thickness/2.0); //cm
-
- this->AddMaterial(matD);
- this->AddMaterial(matBeEndcap1);
- this->AddMaterial(matBeEndcap2);
- }
- //______________________________________________________________________________
-
-
- }
-}
-//______________________________________________________________________________
-
diff --git a/src/core/physics/src/DiffXSec.cxx b/src/core/physics/src/DiffXSec.cxx
deleted file mode 100644
index c2005d9db65543662e501394b4583c892fb73161..0000000000000000000000000000000000000000
--- a/src/core/physics/src/DiffXSec.cxx
+++ /dev/null
@@ -1,429 +0,0 @@
-#include "DiffXSec.h"
-#include "FunctionManager.h"
-
-namespace insane {
-namespace physics {
-
-DiffXSec::DiffXSec() :
- fIncludeJacobian(false),fIsModified(true),fUsePhaseSpace(true),
- fPhaseSpace(nullptr), fID(-1), fBaseID(0),fnDim(3)
-{
- // By default a single particle.
- fMaterialIndex = -1;
- fnDim = 3;
- fPhaseSpace = nullptr;
- fnParticles = 1;
- fnPrimaryParticles = fnParticles;
- fPIDs.push_back(11); //electron
- fTitle = "";
- fPlotTitle = "";
- fLabel = "#frac{d#sigma}{dE'd#Omega}";
- fUnits = "nb/GeV/sr";
- fBeamEnergy = 5.9;
-
- fVariableNames.SetOwner(true);
- fVariableNames.Clear();
-
- fAddedXSections.SetOwner(true);
- fAddedXSections.Clear();
-
- fTargetNucleus = Nucleus::Proton();
-
- //fParticles.SetOwner(true);
- fParticles.Clear();
-
- fHelicity = 0.0;
- fTargetPol.SetMagThetaPhi(0.0,180*degree,0.0);
- fTotalXSec = 0.0;
-
- fStructureFunctions = FunctionManager::GetInstance()->GetStructureFunctions();
- fPolarizedStructureFunctions = FunctionManager::GetInstance()->GetPolarizedStructureFunctions();
- fFormFactors = FunctionManager::GetInstance()->GetFormFactors();
-}
-//_______________________________________________________________________
-
-DiffXSec::~DiffXSec()
-{
- fPIDs.clear();
- fScaleFactors.clear();
- fOffsets.clear();
- fVariableNames.Clear();
- fAddedXSections.Clear();
- fParticles.Clear();
- if (fPhaseSpace) delete fPhaseSpace;
- fPhaseSpace = nullptr;
-}
-//_______________________________________________________________________
-
-DiffXSec::DiffXSec(const DiffXSec& old)
-{
- (*this) = old;
-}
-//_______________________________________________________________________
-
-DiffXSec& DiffXSec::operator=(const DiffXSec& rhs)
-{
- if (this != &rhs) { // make sure not same object
- fIncludeJacobian = rhs.fIncludeJacobian ;
- fIsModified = rhs.fIsModified ;
- fUsePhaseSpace = rhs.fUsePhaseSpace ;
- fPhaseSpace = rhs.fPhaseSpace ;
- //mutable Double_t fVars[15] ;
- fPIDs = rhs.fPIDs ;
- fnDim = rhs.fnDim ;
- fnParticles = rhs.fnParticles ;
- fnParticleVars = rhs.fnParticleVars ;
- //fAddedXSections = rhs.fAddedXSections ;
- fParticles.Clear();
- fParticles.AddAll(&rhs.fParticles) ;
- fVariableNames.Clear();
- fVariableNames.AddAll(&rhs.fVariableNames ) ;
- fBeamEnergy = rhs.fBeamEnergy ;
- fTotalXSec = rhs.fTotalXSec ;
- fHelicity = rhs.fHelicity ;
- fTargetPol = rhs.fTargetPol ;
- fTargetNucleus = rhs.fTargetNucleus ;
- fTargetMaterial = rhs.fTargetMaterial ;
- fMaterialIndex = rhs.fMaterialIndex ;
-
- fScaleFactors = rhs.fScaleFactors ;
- fOffsets = rhs.fOffsets ;
-
- //mutable Double_t fNormalizedVariables[15] ;
- //mutable Double_t fUnnormalizedVariables[15] ;
- //mutable Double_t fDependentVariables[15] ;
-
- fTitle = rhs.fTitle ;
- fPlotTitle = rhs.fPlotTitle ;
- fLabel = rhs.fLabel ;
- fUnits = rhs.fUnits ;
- fStructureFunctions = rhs.fStructureFunctions ;
- fPolarizedStructureFunctions = rhs.fPolarizedStructureFunctions ;
- fFormFactors = rhs.fFormFactors ;
-
- //fnDim = old.fnDim;
- //fnParticles = old.fnParticles;
- //fPIDs = old.fPIDs;
- //fTitle = old.fTitle;
- //fPlotTitle = old.fPlotTitle;
- //fBeamEnergy = old.fBeamEnergy;
- }
- return *this; // Return ref for multiple assignment
-}
-//_______________________________________________________________________
-
-DiffXSec* DiffXSec::Clone(const char * ) const
-{
- std::cout << "DiffXSec::Clone()\n";
- auto * copy = new DiffXSec();
- (*copy) = (*this);
- return copy;
-}
-//______________________________________________________________________________
-
-Double_t DiffXSec::Density(Int_t ndim, Double_t * x) {
- //std::cout << " called Density in dummy base class!!\n;";
- return(EvaluateXSec(GetDependentVariables(GetUnnormalizedVariables(x))));
-}
-
-//_______________________________________________________________________
-Double_t DiffXSec::EvaluateCurrent() const
-{
- for(int i=0;i<fPhaseSpace->GetDimension();i++) fVars[i] = 0;
- /* std::cout << vars[0] << " , " << vars[1] << "\n";*/
- fPhaseSpace->GetEventValues(fVars);
- return(EvaluateXSec(GetDependentVariables(fVars)));
-}
-
-//_______________________________________________________________________
-
-Double_t DiffXSec::EvaluateXSec(const Double_t * x) const
-{
- std::cout << "You need to define EvaluateXSec in the derived class!!" << std::endl;
- return(0);
-}
-//______________________________________________________________________________
-
-Double_t DiffXSec::Jacobian(const Double_t * x) const
-{
- return( TMath::Sin(x[1]));
-}
-//_______________________________________________________________________
-
-bool DiffXSec::VariablesInPhaseSpace(const Int_t ndim, const Double_t * x) const {
- // also checks that the variables are not NaN
- for(int i = 0; i < ndim; i++) {
- if( std::isnan(x[i]) ) {
- Warning("VariablesInPhaseSpace","Variable %d is NaN",i);
- return false;
- }
- }
- if (!fPhaseSpace) {
- std::cout << "[DiffXSec::VariablesInPhaseSpace]: No phase space defined \n";
- return(false);
- }
- fPhaseSpace->SetEventValues(x,fPhaseSpace->GetDimension());//ndim); // allowed because PhaseSpaceVariable::fCurrentValue is mutable
- // ndim is the total number of phase space variables.
- if( !fUsePhaseSpace ) return true;
- bool result = true;
- if (ndim < fnDim) {
- // phase space must have the same or more variables than the differential
- std::cout << " Dimensions Do Not Match! \n";
- std::cout << " VariablesInPhaseSpace argument ndim = " << ndim << " \n";
- std::cout << " DiffXSec::fnDim = " << fnDim << " \n";
- return(false);
- }
- for (int i = 0; i < ndim /*fPhaseSpace->GetDimension()*/; i++) {
- if (!(fPhaseSpace->GetVariable(i)->IsInVariableRange())) {
- result = false;
- //std::cout << "DiffXSec::VariablesInPhaseSpace " << i << " out of range: " << x[i] << std::endl;
- //if(!(fPhaseSpace->GetVariable(i)->IsDependent())) {
- // std::cout << "DiffXSec::VariablesInPhaseSpace - independent variable " << i << " out of range: " << x[i] << std::endl;
- //}
- break;
- }
- }
- return(result);
-}
-//_______________________________________________________________________
-
-Double_t DiffXSec::CalculateTotalCrossSection()
-{
- Double_t result = 0.0;
- Double_t a[30], b[30];
- for (int i = 0; i < GetPhaseSpace()->GetDimension() && i < 30; i++) {
- a[i] = GetPhaseSpace()->GetVariable(i)->GetMinimum();
- b[i] = GetPhaseSpace()->GetVariable(i)->GetMaximum();
- }
- // ROOT::Math::IntegratorMultiDim ig2(ROOT::Math::IntegrationMultiDim::kPLAIN);
- // ROOT::Math::GSLMCIntegrator ig2(ROOT::Math::IntegrationMultiDim::kVEGAS);
- ROOT::Math::IntegratorMultiDim ig2(ROOT::Math::IntegrationMultiDim::kVEGAS);
- ig2.SetFunction((*this));
- /* ig2.SetRelTolerance(0.05);*/
- result = ig2.Integral(a, b);
- std::cout << " DiffXSec::CalculateTotalCrossSection() integral result is " << result << std::endl;
- fTotalXSec = result;
- return(result);
-}
-//_______________________________________________________________________
-void DiffXSec::Print(std::ostream& stream) const {
- //std::cout << " --------------------------------------------------------- " << std::endl;
- stream << " - cross section : " << fTitle.Data() << "\n";
- stream << " title : " << fPlotTitle.Data() << "\n";
- /* std::cout << " Cross Section address : " << this << "\n";*/
- stream << " N Particles : " << fnParticles << "\n";
- stream << " Beam Energy : " << fBeamEnergy << std::endl;
- stream << " Jacobian : " << fIncludeJacobian << std::endl;
- fTargetNucleus.Print(stream);
- stream << " Phase Space : " << "\n";
- if (fPhaseSpace) fPhaseSpace->Print(stream);
- stream << " PIDs : " ;
- for(unsigned int i=0; i< fPIDs.size(); i++) { std::cout << fPIDs[i] << " " ; }
- stream << std::endl;
- stream << " sigma : " << EvaluateCurrent() << " nb/GeV/sr" << std::endl;
- //PrintTable();
- stream << " --------------------------------------------------------- " << std::endl;
-}
-//_______________________________________________________________________
-void DiffXSec::Print(const Option_t * opt) const {
- Print(std::cout);
- ////std::cout << " --------------------------------------------------------- " << std::endl;
- //std::cout << " XSection title : " << fTitle.Data() << "\n";
- //std::cout << " plot title : " << fPlotTitle.Data() << "\n";
- ///* std::cout << " Cross Section address : " << this << "\n";*/
- //std::cout << " fnParticles : " << fnParticles << "\n";
- //std::cout << " Beam Energy : " << fBeamEnergy << std::endl;
- //std::cout << " Jacobian : " << fIncludeJacobian << std::endl;
- //fTargetNucleus.Print();
- //std::cout << " Phase Space : " << "\n";
- //if (fPhaseSpace){
- // fPhaseSpace->Print();
- //}
- //std::cout << " PIDs : " ;
- //for(unsigned int i=0; i< fPIDs.size(); i++) { std::cout << fPIDs[i] << " " ; }
- //std::cout << std::endl;
- //std::cout << " sigma : " << EvaluateCurrent() << " nb/GeV/sr" << std::endl;
- ////PrintTable();
- //std::cout << " --------------------------------------------------------- " << std::endl;
-}
-//_______________________________________________________________________
-void DiffXSec::PrintTable(std::ostream &out){
- if(fPhaseSpace){
- fPhaseSpace->GetEventValues(fVars);
- for(int i = 0; i<fPhaseSpace->GetDimension(); i++) out << fVars[i] << " ";
- out << EvaluateXSec(fVars) << std::endl;
- }
-}
-//_______________________________________________________________________
-void DiffXSec::PrintTable() {
- PrintTable(std::cout);
-}
-//_______________________________________________________________________
-void DiffXSec::PrintGrid(std::ostream &out, Int_t N){
- if(fPhaseSpace){
- std::vector<Double_t> minima;
- std::vector<Double_t> maxima;
- std::vector<Double_t> step;
- std::vector<Int_t> counter;
- for(int i = 0;i<fPhaseSpace->GetDimension(); i++){
- PhaseSpaceVariable * var = fPhaseSpace->GetVariable(i);
- minima.push_back(var->GetMinimum());
- maxima.push_back(var->GetMaximum());
- step.push_back( (maxima[i] - minima[i])/(double(N-1)) );
- counter.push_back(0);
- }
- auto Neval = (unsigned int) TMath::Power(N,fPhaseSpace->GetDimension());
-
- for(unsigned int n = 0; n<Neval ; n++){
- for(unsigned int i = 0; i<counter.size(); i++){
- fVars[i] = minima[i] + double(counter[i])*step[i];
- }
- fPhaseSpace->SetEventValues(fVars);
- PrintTable(out);
- counter[0]++;
- for(unsigned int i = 0; i<counter.size(); i++){
- if(counter[i] == N) {
- counter[i] = 0;
- if(i+1 < N)counter[i+1]++;
- }
- }
- }
- }
-}
-//______________________________________________________________________________
-
-Int_t DiffXSec::GetParticleType(Int_t part) const {
- Int_t N = fPIDs.size();
- if (N > part) return fPIDs.at(part);
- /* else */
- return(0);
-}
-//______________________________________________________________________________
-
-void DiffXSec::SetParticleType(Int_t pdgcode, Int_t part) {
- Int_t N = fPIDs.size();
- if (N > part) fPIDs[part] = pdgcode;
- else fPIDs.push_back(pdgcode);
- fIsModified = true;
-}
-//______________________________________________________________________________
-
-void DiffXSec::SetProductionParticleType( Int_t PDGcode, Int_t part )
-{
- SetParticleType(PDGcode,part);
-}
-//______________________________________________________________________________
-
-void DiffXSec::PrintFinalStatePIDs() const {
- Int_t N = fPIDs.size();
- for (int i = 0; i < N ; i++)
- std::cout << " particle: " << i << " , PID = " << fPIDs.at(i) << "\n";
-}
-//_______________________________________________________________________
-
-Double_t * DiffXSec::GetUnnormalizedVariables(const Double_t * x)
-{
- for (unsigned int i = 0; i < NDim() && i < 30 ; i++) {
- fUnnormalizedVariables[i] = x[i] * fScaleFactors[i] + fOffsets[i];
- }
- return(fUnnormalizedVariables);
-}
-//_______________________________________________________________________
-
-Double_t * DiffXSec::GetNormalizedVariables(const Double_t * x)
-{
- for (unsigned int i = 0; i < NDim() && i < 30 ; i++) {
- fNormalizedVariables[i] = (x[i] - fOffsets[i]) / fScaleFactors[i];
- }
- return(fNormalizedVariables);
-}
-//_______________________________________________________________________
-
-Double_t * DiffXSec::GetDependentVariables(const Double_t * x) const {
- for (unsigned int i = 0; i < NDim() && i < 30 ; i++) {
- fDependentVariables[i] = x[i];
- }
- return(fDependentVariables);
-}
-//_______________________________________________________________________
-void DiffXSec::Refresh() {
-
- if (GetPhaseSpace()) {
- fScaleFactors.clear();
- fOffsets.clear();
- for (int i = 0 ; i < GetPhaseSpace()->GetDimension() ; i++) {
- if(GetPhaseSpace()->GetVariable(i)->IsDependent()) continue;
- fScaleFactors.push_back(GetPhaseSpace()->GetVariable(i)->GetNormScale());
- fOffsets.push_back(GetPhaseSpace()->GetVariable(i)->GetNormOffset());
- }
- GetPhaseSpace()->Refresh();
- SetModified(false);
- } else {
- Error("Refresh","No phase space set!");
- }
-
- InitializeFinalStateParticles();
-}
-//_______________________________________________________________________
-
-void DiffXSec::SetPhaseSpace(PhaseSpace * ps)
-{
- //if( ps != fPhaseSpace) if( ps!=0 ) if( fPhaseSpace != 0) delete fPhaseSpace;
- fPhaseSpace = ps;
- if (fPhaseSpace->GetDimension() < (Int_t)NDim()) Error("SetPhaseSpace", "Phase space dimension less than XSec dimension");
- if (fPhaseSpace->GetNIndependentVariables() > (Int_t)NDim()) Error("SetPhaseSpace", "Phase space dimension greater than XSec dimension");
- fIsModified = true;
-}
-//______________________________________________________________________________
-
-void DiffXSec::InitializeFinalStateParticles()
-{
- // Sets the list of final state particles
- fParticles.Clear();
- if( fnParticles != (Int_t)fPIDs.size() ) {
- Warning("InitializeFinalStateParticles", "warning fnParticles do not have their PIDs defined");
- }
- fnParticleVars.clear();
- for (int i = 0; i < fnParticles; i++) {
-
- Int_t nVars = 0;
- for (int j = 0; j < GetPhaseSpace()->GetDimension(); j++) {
- if( GetPhaseSpace()->GetVariable(j)->GetParticleIndex() == i ) nVars++;
- }
- fnParticleVars.push_back(nVars);
- auto aParticle = new Particle();
- aParticle->SetPdgCode( fPIDs[i] );
- fParticles.Add( aParticle );
- }
-}
-//______________________________________________________________________________
-
-TList * DiffXSec::GetParticles()
-{
- return(&fParticles);
-}
-//______________________________________________________________________________
-
-TParticlePDG * DiffXSec::GetParticlePDG(int i) const
-{
- if( fParticles.GetEntries() > 0 ){
- auto * part = (TParticle*)fParticles.At(i);
- return(part->GetPDG(1)); // mode=1 does not return a new TParticlePDG, it reuses it (don't use with TclonesArray)
- }
- return nullptr;
-}
-//______________________________________________________________________________
-
-//Int_t DiffXSec::GetParticleType(int i = 0) const
-//{
-// if( fParticles.GetEntries() > 0 ){
-// TParticle * part = (TParticle*)fParticles.At(i);
-// return(part->GetPdgCode());
-// }
-// return 0;
-//}
-//______________________________________________________________________________
-
-}
-}
-
diff --git a/src/core/physics/src/EventGenerator.cxx b/src/core/physics/src/EventGenerator.cxx
deleted file mode 100644
index c9c8715c25ba618c7b8a4d8398cc357f195beb8d..0000000000000000000000000000000000000000
--- a/src/core/physics/src/EventGenerator.cxx
+++ /dev/null
@@ -1,487 +0,0 @@
-#include "EventGenerator.h"
-#include <sstream>
-
-namespace insane {
-namespace physics {
-
-EventGenerator::EventGenerator(const char * name, const char * title) : TNamed(name, title) {
-
- fIsModified = true;
- fTotalXSection = 0.0;
- fSamplers.Clear();
- fSamplers.SetOwner(true);
- fBeamEnergy = 5.9;//GeV
- fRandom = new TRandom3(0);
- fEventArray = nullptr;
- fParticles = nullptr;
- fSampler = nullptr;
- fIsBeamRastered = true;
- fSlowRasterSize = 3.0; //cm
- fdzVertex = 2.0; //cm
- fdxVertex = fSlowRasterSize/2.0; //cm
- fdyVertex = fSlowRasterSize/2.0; //cm
-}
-//________________________________________________________________________________
-
-EventGenerator::EventGenerator( const EventGenerator & eg ) : TNamed(eg) {
- fIsModified = eg.fIsModified;
- fTotalXSection = eg.fTotalXSection;
- fSamplers.AddAll( &(eg.fSamplers) );
- fSamplers.SetOwner(true); // not sure about this
- fBeamEnergy = eg.fBeamEnergy;//GeV
- fRandom = new TRandom3(0);
- fIsBeamRastered = eg.fIsBeamRastered;
- fSlowRasterSize = eg.fSlowRasterSize; //cm
- fdxVertex = eg.fdxVertex; //cm
- fdyVertex = eg.fdyVertex; //cm
- fdzVertex = eg.fdzVertex; //cm
-}
-//________________________________________________________________________________
-
-EventGenerator& EventGenerator::operator=(const EventGenerator& eg){
- TNamed::operator=(eg);
- if (this != &eg) {
- fIsModified = eg.fIsModified;
- fTotalXSection = eg.fTotalXSection;
- fSamplers.AddAll( &(eg.fSamplers) );
- fSamplers.SetOwner(true); // not sure about this
- fBeamEnergy = eg.fBeamEnergy;//GeV
- fRandom = new TRandom3(0);
- fIsBeamRastered = eg.fIsBeamRastered;
- fSlowRasterSize = eg.fSlowRasterSize; //cm
- fdxVertex = eg.fdxVertex; //cm
- fdyVertex = eg.fdyVertex; //cm
- fdzVertex = eg.fdzVertex; //cm
- }
- return *this;
-}
-//______________________________________________________________________________
-EventGenerator::~EventGenerator() {
-
-}
-//________________________________________________________________________________
-void EventGenerator::SetBeamEnergy(Double_t E0) {
- fBeamEnergy = E0;
- for(int i = 0; i< fSamplers.GetEntries(); i++) {
- auto * samp = (PhaseSpaceSampler*)fSamplers.At(i);
- DiffXSec * xsec = samp->GetXSec();
- if(xsec) xsec->SetBeamEnergy(E0);
- }
- SetModified(true);
-}
-
-//________________________________________________________________________________
-void EventGenerator::Refresh(bool print) {
- //std::cout << " Refreshing Event Generator\n";
- fTotalXSection = 0.0;
- PhaseSpaceSampler * aSampler = nullptr;
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- /* if( !(aSampler->GetXSec() ) ) aSampler->Load();*/
- //std::cout << "refreshing sampler of " << aSampler->GetXSec()->GetName() << " ... " << std::endl;;
- //std::cout << "solid angle: " << aSampler->GetXSec()->GetPhaseSpace()->GetSolidAngle() << " sr " << std::endl;
- //std::cout << "delta E : " << aSampler->GetXSec()->GetPhaseSpace()->GetEnergyBite() << " GeV " << std::endl;
-
- aSampler->Refresh();
- if( print ) {
- const TargetMaterial& mat = aSampler->GetXSec()->GetTargetMaterial();
- //if(mat)
- std::cout << "material (" << std::setw(10) << mat.GetName() << ") ";
- //std::cout << " Done refreshing PS Sampler " << aSampler->GetName() << std::endl;;
- std::cout << "XS_id: " << aSampler->GetXSec()->GetID() << " (" << std::setw(30) << aSampler->GetXSec()->GetName() << ") = " << aSampler->GetTotalXSection() << " nb" << std::endl;;
- }
- }
- //std::cout << " Calculating the total cross-section ... " << std::endl;
- CalculateTotalCrossSection();
- SetModified(false);
-}
-//________________________________________________________________________________
-
-void EventGenerator::ListSamplers()
-{
- PhaseSpaceSampler * aSampler = nullptr;
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- aSampler->Print();
- }
-}
-//________________________________________________________________________________
-
-void EventGenerator::PrintXSecSummary(std::ostream& s)
-{
- if(IsModified()) Refresh();
- PhaseSpaceSampler * aSampler = nullptr;
- Double_t total = 0;
- TListIter samplerIter(&fSamplers);
- s
- << std::setw(11) << "XS_ID"
- << std::setw(10) << "PDG_codes"
- << std::setw(20) << "XS (nb)"
- << std::endl;;
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
-
- Double_t tot = aSampler->GetTotalXSection();
- total += tot;
- //Int_t part = aSampler->GetXSec()->GetParticleType();
- std::stringstream ss_pdg;
- int npart = aSampler->GetXSec()->GetParticles()->GetEntries();
- for(int ipdg = 0; ipdg < npart; ipdg++) {
- if(ipdg != 0 ) ss_pdg << ", ";
- ss_pdg << aSampler->GetXSec()->GetParticleType(ipdg) ;
- }
- std::string part = ss_pdg.str();
-
- s
- << "{"
- << std::setw(11) << aSampler->GetXSec()->GetID()
- << ", {"
- << std::setw(10) << part
- << "}, "
- << std::setw(20) << std::right << std::setprecision(6) << std::scientific << tot << std::defaultfloat << std::left
- << std::endl;
- //std::cout << part << " " << tot << " nb" << std::endl;
- //aSampler->Refresh();
- }
- s << "-------------------------------" << std::endl;
- s
- << std::setw(21) << " "
- << std::setw(20) << total
- << std::endl;
- //std::cout << " total : " << total << " nb" << std::endl;
-
-}
-//________________________________________________________________________________
-
-Double_t EventGenerator::GetTotalCrossSectionForParticle(Int_t pid )
-{
- if(IsModified()) Refresh();
- if( pid==0 ) return fTotalXSection;
- PhaseSpaceSampler * aSampler = nullptr;
- Double_t total = 0;
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
-
- Double_t tot = aSampler->GetTotalXSection();
- Int_t part = aSampler->GetXSec()->GetParticleType();
- if(pid == part) total += tot;
- //aSampler->Refresh();
- }
- return total;
-}
-//________________________________________________________________________________
-
-void EventGenerator::ListPhaseSpaceVariables()
-{
- PhaseSpaceSampler * aSampler = nullptr;
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- aSampler->GetXSec()->GetPhaseSpace()->Print();
- }
-}
-//______________________________________________________________________________
-
-void EventGenerator::InitFromDisk()
-{
- TListIter samplerIter(&fSamplers);
- std::cout << "init from disk\n";
- PhaseSpaceSampler * aSampler = nullptr;
- while (( aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- aSampler->InitFromDisk();
- std::cout << "done init sampler\n";
- }
-}
-//________________________________________________________________________________
-
-void EventGenerator::Print(std::ostream& stream) const
-{
- stream << "================================================================================" << std::endl;
- stream << "EG: " << GetName() << std::endl;
- stream << "================================================================================" << std::endl;
- stream << " Total cross section = " <<
- std::scientific << fTotalXSection << " [nb] = " <<
- std::scientific << fTotalXSection*1.0e-9 << " [b] = " <<
- std::scientific << fTotalXSection*1.0e-33 << " [cm^2] " << std::endl;
- PhaseSpaceSampler * aSampler = nullptr;
- TListIter samplerIter(&fSamplers);
- int i = 0;
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- stream << " ----------------------------------------------------------------------------" << std::endl;
- stream << " PS Sampler #" << i << std::endl;
- aSampler->Print(stream);
- i++;
- }
-}
-//________________________________________________________________________________
-
-void EventGenerator::Print(const Option_t * opt) const
-{
- EventGenerator::Print(std::cout);
- //std::cout << "--------------------------------------------------------------------------------" << std::endl;
- //std::cout << "Event Generator : " << GetName() << std::endl;
- //std::cout << "Total Cross Section : " << fTotalXSection << " nb" << std::endl;
- //PhaseSpaceSampler * aSampler = 0;
- //TListIter samplerIter(&fSamplers);
- //int i = 1;
- //while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- // std::cout << " o----------------------------------------------------------------------------" << std::endl;
- // std::cout << " o phase space sampler " << i << std::endl;
- // aSampler->Print();
- // i++;
- //}
-}
-//________________________________________________________________________________
-
-Double_t EventGenerator::CalculateTotalCrossSection()
-{
- PhaseSpaceSampler * aSampler = nullptr;
- fTotalXSection = 0.0;
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- fTotalXSection += aSampler->CalculateTotalXSection();
- }
- //printf(" EventGenerator total Xsection is %f nb.\n", fTotalXSection);
- std::cout << " material (" << std::setw(20) << "All" << "), ";
- std::cout << " cross section (" << std::setw(30) << "Total" << ") = " << fTotalXSection << " nb" << std::endl;
- return(fTotalXSection);
-}
-
-//________________________________________________________________________________
-
-void EventGenerator::AddSampler(PhaseSpaceSampler * ps)
-{
- if(ps) fSamplers.Add(ps);
- else Error("AddSampler","Sampler is null!");
-}
-
-//________________________________________________________________________________
-
-void EventGenerator::Clear(Option_t * opt)
-{
- fSamplers.Clear(opt);
- fTotalXSection = 0.0;
-}
-//________________________________________________________________________________
-
-PhaseSpaceSampler * EventGenerator::GetRandomSampler()
-{
- if (fSamplers.GetEntries() == 0) {
- Error("GetRandomSampler","Event generator has no samplers.");
- }
- PhaseSpaceSampler * aSampler = nullptr;
- PhaseSpaceSampler * selectedSampler = nullptr;
- Double_t rand = fRandom->Uniform(fTotalXSection);
- TListIter samplerIter(&fSamplers);
- while ((aSampler = (PhaseSpaceSampler*)samplerIter.Next())) {
- rand = rand - aSampler->GetTotalXSection();
- if (rand < 0.0) {
- selectedSampler = aSampler;
- break;
- }
- }
- return(selectedSampler);
-}
-//________________________________________________________________________________
-
-TList * EventGenerator::GetSamplers()
-{
- return (&fSamplers);
-}
-//________________________________________________________________________________
-
-TList * EventGenerator::GetPSVariables(const char * var )
-{
- auto * vars = new TList();
- TList * samps = GetSamplers();
- for(int i = 0;i<samps->GetEntries(); i++) {
- auto * aSamp = (PhaseSpaceSampler*)samps->At(i);
- DiffXSec * xsec = aSamp->GetXSec();
- if(xsec){
- PhaseSpace * ps = xsec->GetPhaseSpace();
- if( ps ) {
- PhaseSpaceVariable * aVar = ps->GetVariableWithName( var );
- if(aVar)
- vars->Add(aVar);
- }
- }
- }
- return(vars);
-}
-//________________________________________________________________________________
-
-TLorentzVector EventGenerator::GetRandomVertex(Int_t i) {
-
- Double_t x = 0;
- Double_t y = 0;
- Double_t z = 0;
- if( i >= 0 ){
- // Use cross section's target material z position
- //fSampler->GetXSec()->GetTargetMaterial()->Dump();
- Double_t zpos = fSampler->GetXSec()->GetTargetMaterial().fZposition;
- Double_t l = fSampler->GetXSec()->GetTargetMaterial().fLength;
- z = fRandom->Uniform(zpos-l/2.0, zpos+l/2.0);
- } else {
- //
- z = fRandom->Uniform(-fdzVertex, fdzVertex);
- }
- if( fIsBeamRastered ) {
- fdxVertex = fSlowRasterSize/2.0; //cm
- fdyVertex = fSlowRasterSize/2.0; //cm
- do {
- //fRandom->Circle(x,y,fSlowRasterSize/2.0);
- x = fRandom->Uniform(-fdxVertex, fdxVertex);
- y = fRandom->Uniform(-fdyVertex, fdyVertex);
-
- } while ( x*x + y*y > (fSlowRasterSize/2.0)*(fSlowRasterSize/2.0)) ;
- } else {
- x = fRandom->Uniform(-fdxVertex, fdxVertex);
- y = fRandom->Uniform(-fdyVertex, fdyVertex);
- }
- TLorentzVector vertex(x,y,z,0.0);
- return(vertex);
-}
-//________________________________________________________________________________
-
-bool EventGenerator::NeedsRefreshed()
-{
- int res = 0;
- PhaseSpaceSampler * aSampler = nullptr;
- for (int i = 0; i < fSamplers.GetEntries(); i++) {
- aSampler = (PhaseSpaceSampler*)fSamplers.At(i);
- res += aSampler->IsModified();
- DiffXSec * aXSec = aSampler->GetXSec();
- if (aXSec) res += aXSec->IsModified();
- PhaseSpace * aPS = aXSec->GetPhaseSpace();
- if (aPS) res += aPS->IsModified();
- PhaseSpaceVariable * aPSVar = nullptr;
- for (int j = 0; j < aPS->GetDimension(); j++) {
- aPSVar = aPS->GetVariable(j);
- if (aPSVar) res += aPSVar->IsModified();
- }
- }
- fIsModified = res;
- return(res);
-}
-//________________________________________________________________________________
-
-TList * EventGenerator::GenerateEvent()
-{
- if( IsModified() ) Refresh();
-
- fSampler = GetRandomSampler();
-
- if(!fSampler) {
- Error("GenerateEvent","Null Sampler");
- return(nullptr);
- }
-
- fEventArray = fSampler->GenerateEvent();
- fParticles = fSampler->GetXSec()->GetParticles();
- Int_t matIndex = fSampler->GetXSec()->GetTargetMaterialIndex(); // index of associated target material, -1 if none
- Int_t xsec_index = fSamplers.IndexOf(fSampler);
-
- auto rand_vertex = GetRandomVertex(matIndex);
-
- // this is now done by a mandatory class DiffXSec::DefineEvent() which
- // is called from PhaseSpaceSampler::GenerateEvent();
- for (int i = 0; i < fParticles->GetEntries(); i++) {
- TParticle * apart = ((TParticle*)fParticles->At(i));
- apart->SetProductionVertex( rand_vertex );
- //anpart->SetStatusCode(matIndex);
- apart->SetFirstMother(matIndex);
- apart->SetStatusCode(xsec_index);
- }
-
- fEG_Event.SetPSVariables(fEventArray, fSampler->GetXSec()->GetPhaseSpace()->GetDimension());
- fEG_Event.SetParticles(fParticles);
- fEG_Event.fXSec = fSampler->GetXSec();
- fEG_Event.fXS_id = fSampler->GetXSec()->GetID();
- fEG_Event.fBeamPol = fSampler->GetXSec()->GetHelicity();
-
- return fParticles;
-}
-//________________________________________________________________________________
-
-void EventGenerator::LundFormat( std::ostream& s )
-{
- LundHeaderFormat(fEG_Event,s);
- for(int i = 0; i<fEG_Event.fParticles->GetEntries(); i++ ){
- LundEventFormat(i,fEG_Event,s);
- }
-}
-//________________________________________________________________________________
-
-void EventGenerator::LundHeaderFormat(
- _EG_Event& eg_event,
- std::ostream& s)
-{
- // Lund format reference: https://gemc.jlab.org/gemc/Documentation/Entries/2011/3/18_The_LUND_Format.html
- // Header:
- // Column Quantity
- // 1 Number of particles
- // 2* Number of target nucleons
- // 3* Number of target protons
- // 4* Target Polarization
- // 5 Beam Polarization
- // 6* x
- // 7* y
- // 8* W
- // 9* Q2
- // 10* nu
- s << eg_event.fParticles->GetEntries() << " "
- << eg_event.fXSec->GetA() << " "
- << eg_event.fXSec->GetZ() << " "
- << eg_event.fTargetPol << " "
- << eg_event.fBeamPol << " "
- << 0.0 /* x */ << " "
- << 0.0 /* y */ << " "
- << 0.0 /* W */ << " "
- << 0.0 /* Q2 */ << " "
- << 0.0 /* nu */ << "\n";
-}
-//________________________________________________________________________________
-void EventGenerator::LundEventFormat(
- int i,
- _EG_Event& eg_event,
- std::ostream& s)
-{
- // Lund format reference: https://gemc.jlab.org/gemc/Documentation/Entries/2011/3/18_The_LUND_Format.html
- // Particles:
- // Column Quantity
- // 1 index
- // 2 charge
- // 3 type(=1 is active)
- // 4 particle id
- // 5 parent id (decay bookkeeping)
- // 6 daughter (decay bookkeeping)
- // 7 px [GeV]
- // 8 py [GeV]
- // 9 pz [GeV]
- // 10 E [GeV]
- // 11 mass (not used)
- // 12 x vertex [m]
- // 13 y vertex [m]
- // 14 z vertex [m]
- // NOTE the vertex is in METERS unlike the website above
- int iLund = i+1;
- auto * part = (Particle*)eg_event.fParticles->At(i);
- TParticlePDG * part_pdg = part->GetPDG();
- s << " "
- << iLund << " "
- << part_pdg->Charge()/3.0 << " "
- << 1 << " "
- << part_pdg->PdgCode() << " "
- << 0 << " "
- << 0 << " "
- << part->Px() << " "
- << part->Py() << " "
- << part->Pz() << " "
- << part->Energy() << " "
- << part_pdg->Mass() << " "
- << part->Vx()/100.0 << " "
- << part->Vy()/100.0 << " "
- << part->Vz()/100.0 << "\n";
-}
-//______________________________________________________________________________
-
-}
-}
diff --git a/src/core/physics/src/ExclusiveDiffXSec.cxx b/src/core/physics/src/ExclusiveDiffXSec.cxx
deleted file mode 100644
index 890b0637a45437797a39885f4e113ac5007fbeaf..0000000000000000000000000000000000000000
--- a/src/core/physics/src/ExclusiveDiffXSec.cxx
+++ /dev/null
@@ -1,183 +0,0 @@
-#include "ExclusiveDiffXSec.h"
-#include "TMath.h"
-#include "Math/Integrator.h"
-#include "Math/IntegratorMultiDim.h"
-#include "Math/AllIntegrationTypes.h"
-#include "Math/Functor.h"
-#include "Math/GaussIntegrator.h"
-#include "StructureFunctions.h"
-#include "TCanvas.h"
-#include "TGraph2D.h"
-#include "TGraph.h"
-#include "TStyle.h"
-#include "TVirtualPad.h"
-#include "TAxis.h"
-
-namespace insane {
-namespace physics {
-
-
-//________________________________________________________________________
-ExclusiveDiffXSec::ExclusiveDiffXSec()
-{
- fID = 2000;
- fnDim = 2;
- fnParticles = 2;
- fPIDs.push_back(2212);//proton
- //InitializePhaseSpaceVariables();
-}
-
-//________________________________________________________________________
-ExclusiveDiffXSec::~ExclusiveDiffXSec() {
-}
-//______________________________________________________________________________
-void ExclusiveDiffXSec::DefineEvent(Double_t * vars) {
-
- Int_t totvars = 0;
- for (int i = 0; i < fParticles.GetEntries(); i++) {
- /// \todo fix this hard coding of 3 variables per event.
- /// here we are assuming the order E,theta,phi,then others
- /// \todo figure out how to handle vertex.
- insane::Kine::SetMomFromEThetaPhi((TParticle*)(fParticles.At(i)), &vars[totvars]);
- //((TParticle*)fParticles.At(i))->SetProductionVertex(GetRandomVertex());
- totvars += GetNParticleVars(i);
- }
-
-}
-
-//________________________________________________________________________
-Double_t ExclusiveDiffXSec::EvaluateXSec(const Double_t * x) const
-{
-
- /// Return zero if the variables are not in the phase space
- if (!VariablesInPhaseSpace(fnDim, x)) return(0.0);
-
- Double_t Eprime = x[0];
- Double_t theta = x[1];
- //Double_t phi = x[2];
- Double_t Qsquared = 4.0 * GetBeamEnergy() * Eprime * TMath::Power(TMath::Sin(theta / 2.0), 2);
- Double_t xbjorken = Qsquared / (2.0 * 0.938 * (GetBeamEnergy() - Eprime));
- // MeV and degrees
- /*Double_t hbarcSquaredperbarn = 6.5822e-22*6.5822e-22*9.0e16/(1.0e-28); // GeV^2 ubarn */
- Double_t hbarc2 = 0.38939129; /*(hbar*c)^2 = 0.38939129 GeV^2 mbarn */
- Double_t mottXSec = hbarc2 * (1. / 137.) * (1. / 137.) *
- TMath::Power(TMath::Cos(theta / 2.0), 2) /
- (4.0 * GetBeamEnergy() * GetBeamEnergy() * TMath::Power(TMath::Sin(theta / 2.0), 4))
- * (1.0 / 0.938);
- return(mottXSec);
-}
-//______________________________________________________________________________
-
-void ExclusiveDiffXSec::InitializePhaseSpaceVariables() {
-
- PhaseSpace * ps = GetPhaseSpace();
- if(ps) delete ps;
- ps = nullptr;
-
- if (!ps) {
- std::cout << " Creating NEW PhaseSpace for ExclusiveDiffXSec\n";
- ps = new PhaseSpace();
-
- // ------------------------------
- // Electron
- auto * varEnergy = new PhaseSpaceVariable();
- varEnergy = new PhaseSpaceVariable();
- varEnergy->SetNameTitle("energy_e", "E_{e'}");
- varEnergy->SetMinimum(0.50);
- varEnergy->SetMaximum(GetBeamEnergy());
- varEnergy->SetDependent(true);
- ps->AddVariable(varEnergy);
-
- auto * varTheta = new PhaseSpaceVariable();
- varTheta->SetNameTitle("theta_e", "#theta_{e'}");
- varTheta->SetMinimum(20.0 * TMath::Pi() / 180.0);
- varTheta->SetMaximum(160.0 * TMath::Pi() / 180.0);
- ps->AddVariable(varTheta);
-
- auto * varPhi = new PhaseSpaceVariable();
- varPhi->SetNameTitle("phi_e", "#phi_{e'}");
- varPhi->SetMinimum(-10.0 * TMath::Pi()/180.0 );
- varPhi->SetMaximum( 10.0 * TMath::Pi()/180.0 );
- varPhi->SetUniform(true);
- ps->AddVariable(varPhi);
-
-
- // ------------------------------
- // Proton
- auto * varEnergy_p = new PhaseSpaceVariable();
- varEnergy_p = new PhaseSpaceVariable();
- varEnergy_p->SetNameTitle("energy_p", "E_{p}");
- varEnergy_p->SetMinimum(0.9);
- varEnergy_p->SetMaximum(5.0);
- varEnergy_p->SetParticleIndex(1);
- varEnergy_p->SetDependent(true);
- ps->AddVariable(varEnergy_p);
-
- auto * varTheta_p = new PhaseSpaceVariable();
- varTheta_p->SetNameTitle("theta_p", "#theta_{p}");
- varTheta_p->SetMinimum(0.5*degree);
- varTheta_p->SetMaximum(170*degree);
- varTheta_p->SetParticleIndex(1);
- varTheta_p->SetDependent(true);
- ps->AddVariable(varTheta_p);
-
- auto * varPhi_p = new PhaseSpaceVariable();
- varPhi_p->SetNameTitle("phi_p", "#phi_{p}");
- varPhi_p->SetMinimum(- 1.0*TMath::Pi() );
- varPhi_p->SetMaximum( 2.0*TMath::Pi() );
- varPhi_p->SetParticleIndex(1);
- varPhi_p->SetDependent(true);
- varPhi_p->SetInverted(true);
- ps->AddVariable(varPhi_p);
-
- SetPhaseSpace(ps);
- } else {
- std::cout << " Using existing phase space variables !\n";
- }
-}
-//______________________________________________________________________________
-
-
-
-
-//______________________________________________________________________________
-Double_t FlatExclusiveDiffXSec::EvaluateXSec(const Double_t * x) const
-{
- if (!VariablesInPhaseSpace(fnDim, x)) return(0.0);
- /// Returns a constant
- return(1.0);
-}
-//______________________________________________________________________________
-
-
-
-
-//______________________________________________________________________________
-Double_t ExclusiveMottXSec::EvaluateXSec(const Double_t * x) const
-{
- /// Get the recoiling proton momentum
- if (GetBeamEnergy() < GetEPrime(x[1])) return(0.0);
-
-// TVector3 k1(0,0,GetBeamEnergy()); // incident electron
-// TVector3 k2(0,0,0); // scattered electron
-// k2.SetMagThetaPhi(GetEPrime(x[1]),x[1],x[2]);
-// TVector3 p2 = k1-k2; // recoil proton
-// Double_t y[6] = {x[0],x[1],x[2],x[3],x[4],x[5]};
-//
-// Double_t * z = GetDependentVariables(y);
-
- ///
- if (!VariablesInPhaseSpace(fnDim, x)) return(0.0);
-
- Double_t Eprime = x[0];
- Double_t theta = x[1];
- Double_t hbarc2 = 0.38939129; /*(hbar*c)^2 = 0.38939129 GeV^2 mbarn */
- Double_t mottXSec = hbarc2 * (1. / 137.) * (1. / 137.) *
- TMath::Power(TMath::Cos(theta / 2.0), 2) /
- (4.0 * GetBeamEnergy() * GetBeamEnergy() * TMath::Power(TMath::Sin(theta / 2.0), 4))
- * (Eprime / GetBeamEnergy());
- return(mottXSec);
-}
-
-
-}}
diff --git a/src/core/physics/src/InclusiveDiffXSec.cxx b/src/core/physics/src/InclusiveDiffXSec.cxx
deleted file mode 100644
index d748973a04d9748deeda2645bf3e31764907ed7c..0000000000000000000000000000000000000000
--- a/src/core/physics/src/InclusiveDiffXSec.cxx
+++ /dev/null
@@ -1,333 +0,0 @@
-#include "InclusiveDiffXSec.h"
-#include "TMath.h"
-#include "Math/Integrator.h"
-#include "Math/IntegratorMultiDim.h"
-#include "Math/AllIntegrationTypes.h"
-#include "Math/Functor.h"
-#include "Math/GaussIntegrator.h"
-//#include "StructureFunctions.h"
-#include "TCanvas.h"
-#include "TGraph2D.h"
-#include "TGraph.h"
-#include "TStyle.h"
-#include "TVirtualPad.h"
-#include "TAxis.h"
-
-namespace insane {
-namespace physics {
-
-
-InclusiveDiffXSec::InclusiveDiffXSec() : fPolType(0)
-{
- fID = 100000000;
- SetIncludeJacobian(false);
-}
-//______________________________________________________________________________
-
-InclusiveDiffXSec::InclusiveDiffXSec(const InclusiveDiffXSec& old) :
- DiffXSec(old)
-{
- (*this) = old;
-}
-//______________________________________________________________________________
-
-InclusiveDiffXSec::~InclusiveDiffXSec()
-{
-}
-//______________________________________________________________________________
-
-InclusiveDiffXSec& InclusiveDiffXSec::operator=(const InclusiveDiffXSec& old)
-{
- if (this != &old) { // make sure not same object
- DiffXSec::operator=(old);
- fPolType = old.fPolType;
- }
- return *this; // Return ref for multiple assignment
-}
-//______________________________________________________________________________
-
-InclusiveDiffXSec* InclusiveDiffXSec::Clone(const char * newname) const
-{
- std::cout << "InclusiveDiffXSec::Clone()\n";
- auto * copy = new InclusiveDiffXSec();
- (*copy) = (*this);
- return copy;
-}
-//______________________________________________________________________________
-
-InclusiveDiffXSec* InclusiveDiffXSec::Clone() const
-{
- return( Clone("") );
-}
-//______________________________________________________________________________
-
-double InclusiveDiffXSec::DoEval(const double * x) const
-{
- return (double)EvaluateXSec((Double_t*) x);
-}
-//______________________________________________________________________________
-
-Double_t InclusiveDiffXSec::EvaluateXSec(const Double_t * x) const
-{
-
- if (!VariablesInPhaseSpace(fnDim, x)) return(0.0);
-
- // for inelastic scattering
- Double_t M = M_p/GeV; // in GeV
- Double_t Ep = x[0]; // scattered electron energy in GeV
- Double_t th = x[1]; // electron scattering angle in radians
- Double_t Nu = fBeamEnergy - Ep;
- Double_t SIN = TMath::Sin(th/2.);
- Double_t SIN2 = SIN*SIN;
- Double_t COS2 = 1. - SIN2;
- Double_t TAN2 = SIN2/COS2;
- Double_t Q2 = insane::Kine::Qsquared(fBeamEnergy,Ep,th);
- Double_t xBj = insane::Kine::xBjorken_EEprimeTheta(fBeamEnergy,Ep,th);;
-
- // Calculate structure functions
- Double_t F1,F2;
-
- //Nucleus::NucleusType Target = fTargetNucleus.GetType();
- //fTargetNucleus.Dump();
-
- if( fTargetNucleus == Nucleus::Proton() ){
- F1 = fStructureFunctions->F1p(xBj,Q2);
- F2 = fStructureFunctions->F2p(xBj,Q2);
- } else if(fTargetNucleus == Nucleus::Neutron() ) {
- F1 = fStructureFunctions->F1n(xBj,Q2);
- F2 = fStructureFunctions->F2n(xBj,Q2);
- } else {
- Error("EvaluateXSec","Invalid target (or not functional yet). ");
- std::cout << " Use CompositeDiffXSec Instead ... using proton for now " << std::endl;
- F1 = fStructureFunctions->F1p(xBj,Q2);
- F2 = fStructureFunctions->F2p(xBj,Q2);
- }
- Double_t W1 = (1./M)*F1;
- Double_t W2 = (1./Nu)*F2;
- // compute the Mott cross section (units = nb):
- Double_t alpha = fine_structure_const;
- Double_t num = alpha*alpha*COS2;
- Double_t den = 4.*fBeamEnergy*fBeamEnergy*SIN2*SIN2;
- Double_t MottXS = num/den;
- // compute the full cross section (units = nb/GeV/sr)
- Double_t fullXsec = MottXS*(W2 + 2.0*TAN2*W1)*hbarc2_gev_nb;
-
- if ( IncludeJacobian() ) return fullXsec*TMath::Sin(th);
- return fullXsec;
-}
-//______________________________________________________________________________
-
-Double_t InclusiveDiffXSec::Density(Int_t ndim, Double_t * x)
-{
- // std::cout << "Density " << x[0] << " " << x[1] << " " << x[2] << " \n";
- //Double_t y[3]={x[0]*fScaleFactors[0]+fOffsets[0],x[1]*fScaleFactors[1]+fOffsets[1],x[2]*fScaleFactors[2]+fOffsets[2]};
- return(EvaluateXSec(GetDependentVariables(GetUnnormalizedVariables(x))));
-}
-//______________________________________________________________________________
-
-void InclusiveDiffXSec::DefineEvent(Double_t * vars)
-{
-
- Int_t totvars = 0;
- for (int i = 0; i < fParticles.GetEntries(); i++) {
- /// \todo fix this hard coding of 3 variables per event.
- /// here we are assuming the order E,theta,phi,then others
- /// \todo figure out how to handle vertex.
- insane::Kine::SetMomFromEThetaPhi((TParticle*)(fParticles.At(i)), &vars[totvars]);
- //((TParticle*)fParticles.At(i))->SetProductionVertex(GetRandomVertex());
- totvars += GetNParticleVars(i);
- }
-
-}
-//______________________________________________________________________________
-
-void InclusiveDiffXSec::InitializePhaseSpaceVariables()
-{
- //std::cout << " o InclusiveDiffXSec::InitializePhaseSpaceVariables() \n";
-
- PhaseSpace * ps = GetPhaseSpace();
- if (ps) delete ps;
- ps = nullptr;
- if (!ps) {
- // std::cout << " o creating new PhaseSpace\n";
-
- ps = new PhaseSpace();
-
- auto * varEnergy = new PhaseSpaceVariable("energy_e", "E_{e'}", 0.5, 5.0);
- ps->AddVariable(varEnergy);
-
- auto * varTheta = new PhaseSpaceVariable("theta_e", "#theta_{e'}", 25.0*degree, 60.0*degree );
- ps->AddVariable(varTheta);
-
- auto * varPhi = new PhaseSpaceVariable("phi_e", "#phi_{e'}",-180.0*degree,180.0*degree );
- varPhi->SetUniform(true);
- ps->AddVariable(varPhi);
-
- SetPhaseSpace(ps);
-
- }
-}
-//______________________________________________________________________________
-
-double InclusiveDiffXSec::Vs_W(double *x, double *p)
-{
- Double_t E0 = GetBeamEnergy();
- Double_t W = x[0];
- Double_t th = p[0];
- fFuncArgs[0] = insane::Kine::Eprime_W2theta(W*W,th,E0);
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-
-double InclusiveDiffXSec::EnergyDepXSec(double *x, double *p)
-{
- fFuncArgs[0] = x[0];
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-
-double InclusiveDiffXSec::EnergyFractionDepXSec(double *x, double *p)
-{
- fFuncArgs[0] = (1.0 - x[0])*GetBeamEnergy();
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-
-double InclusiveDiffXSec::EnergyDependentXSec(double *x, double *p)
-{
- fFuncArgs[0] = x[0];
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::W2DependentXSec(double *x, double *p) {
- Double_t W2 = x[0];
- Double_t Eprime = insane::Kine::Eprime_W2theta(W2,p[0],GetBeamEnergy());
- fFuncArgs[0] = Eprime;
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::WDependentXSec(double *x, double *p) {
- Double_t W = x[0];
- Double_t Q2 = p[0];
- Double_t xbj = insane::Kine::xBjorken_WQsq(W, Q2);
- Double_t E0 = GetBeamEnergy();
- Double_t y_frac = Q2/(2.0*(M_p/GeV)*xbj*E0);
- if( y_frac >1.0 ) return 0.0;
- Double_t Eprime = insane::Kine::Eprime_xQ2y(xbj,Q2,y_frac);
- Double_t theta = insane::Kine::Theta_xQ2y(xbj,Q2,y_frac);
- if( TMath::IsNaN(theta) ) return 0.0;
- if( TMath::IsNaN(Eprime) ) return 0.0;
- //std::cout << "W: " << W << ", Q2: " << Q2 << std::endl;
- //std::cout << "theta: " << theta << ", Eprime: " << Eprime ;
- //std::cout << ", y_frac: " << y_frac << ", x: " << xbj << std::endl;
- fFuncArgs[0] = Eprime;
- fFuncArgs[1] = theta;
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::xDependentXSec(double *x, double *p) {
- Double_t x_bj = x[0];
- Double_t Q2 = p[0];
- Double_t E0 = GetBeamEnergy();
- //Double_t phi = p[1];
- Double_t xmin = Q2/(2.0*(M_p/GeV)*E0);
- if(x_bj <= xmin) return 0.0;
- Double_t y_frac = Q2/(2.0*(M_p/GeV)*x_bj*E0);
- Double_t Eprime = insane::Kine::Eprime_xQ2y(x_bj,Q2,y_frac);
- Double_t theta = insane::Kine::Theta_xQ2y(x_bj,Q2,y_frac);
- //std::cout << "theta: " << theta << ", Eprime: " << Eprime << std::endl;
- //std::cout << "y_frac: " << y_frac << ", x: " << x_bj << std::endl;
- if(Eprime < 0.0) return 0.0;
- if(y_frac < 0.0) return 0.0;
- if(y_frac >= 1.0) return 0.0;
- fFuncArgs[0] = Eprime;
- fFuncArgs[1] = theta;
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::PhotonEnergyDependentXSec(double *x, double *p) {
- Double_t pEnergy = GetBeamEnergy() - x[0];
- fFuncArgs[0] = pEnergy;
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::MomentumDependentXSec(double *x, double *p) {
- TParticle * aParticle = nullptr;
- if(fParticles.GetEntries() > 0) aParticle = (TParticle*)fParticles.At(0);
- Double_t Eprime = x[0];
- if(aParticle) Eprime = TMath::Sqrt( x[0]*x[0] + TMath::Power(aParticle->GetMass(),2.0));
- fFuncArgs[0] = Eprime;
- fFuncArgs[1] = p[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::PolarAngleDependentXSec(double *x, double *p) {
- fFuncArgs[0] = p[0];
- fFuncArgs[1] = x[0];
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::PolarAngleDependentXSec_deg(double *x, double *p) {
- fFuncArgs[0] = p[0];
- fFuncArgs[1] = x[0]*degree;
- fFuncArgs[2] = p[1];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-double InclusiveDiffXSec::AzimuthalAngleDependentXSec(double *x, double *p) {
- fFuncArgs[0] = p[0];
- fFuncArgs[1] = p[1];
- fFuncArgs[2] = x[0];
- return(EvaluateXSec(GetDependentVariables(fFuncArgs)));
-}
-//______________________________________________________________________________
-
-
-
-
-
-//______________________________________________________________________________
-
-FlatInclusiveDiffXSec::FlatInclusiveDiffXSec()
-{
- fID = 100000001;
-}
-//______________________________________________________________________________
-
-FlatInclusiveDiffXSec::~FlatInclusiveDiffXSec()
-{ }
-//______________________________________________________________________________
-
-Double_t FlatInclusiveDiffXSec::EvaluateXSec(const Double_t * x) const {
-
- if (!VariablesInPhaseSpace(fnDim, x)) return(0.0);
-
- Double_t th = x[1];
-
- Double_t fullXsec = 1.0*hbarc2_gev_nb;
-
- if ( IncludeJacobian() ) {
- //std::cout << " using jacobian\n";
- return fullXsec*TMath::Sin(th);
- }
-
- return(fullXsec);
-}
-//______________________________________________________________________________
-}}
diff --git a/src/core/physics/src/PhaseSpace.cxx b/src/core/physics/src/PhaseSpace.cxx
deleted file mode 100644
index 6a4ddf456e1c1b3a20cc90108fe1fa4c41524942..0000000000000000000000000000000000000000
--- a/src/core/physics/src/PhaseSpace.cxx
+++ /dev/null
@@ -1,29 +0,0 @@
-#include "PhaseSpace.h"
-
-#include "TString.h"
-#include "TClonesArray.h"
-#include <iostream>
-#include <vector>
-
-
-namespace insane {
-namespace physics {
-
-//______________________________________________________________________________
-void PhaseSpace::Print(Option_t * opt ) const {
- //std::cout << " Phase Space Address : " << this << "\n";
- std::cout << std::setw(25) << std::right << "Phase Space Dimension : " << fnDim << "\n";
- const int N = fVariables.GetEntries();
- if (N != fnDim) std::cout << " ARRAY SIZE AND DIMENSION DO NOT MATCH!!!!\n";
- for (int i = 0; i < N; i++)((PhaseSpaceVariable*)fVariables.At(i))->Print();
-}
-//______________________________________________________________________________
-void PhaseSpace::Print(std::ostream& stream) const {
- //std::cout << " Phase Space Address : " << this << "\n";
- stream << std::setw(25) << std::right << "Phase Space Dimension : " << fnDim << "\n";
- const int N = fVariables.GetEntries();
- if (N != fnDim) stream << " ARRAY SIZE AND DIMENSION DO NOT MATCH!!!!\n";
- for (int i = 0; i < N; i++)((PhaseSpaceVariable*)fVariables.At(i))->Print(stream);
-}
-//______________________________________________________________________________
-}}
diff --git a/src/core/physics/src/PhaseSpaceSampler.cxx b/src/core/physics/src/PhaseSpaceSampler.cxx
deleted file mode 100644
index b5fc9e0ca6edbff6b74e39a2d0e86ab65fc71098..0000000000000000000000000000000000000000
--- a/src/core/physics/src/PhaseSpaceSampler.cxx
+++ /dev/null
@@ -1,233 +0,0 @@
-#include "PhaseSpaceSampler.h"
-
-#include "PhaseSpace.h"
-#include "TH2D.h"
-#include "TRandom3.h"
-#include "TFoam.h"
-#include "TCanvas.h"
-#include "PhaseSpace.h"
-#include "InclusiveDiffXSec.h"
-#include "Math/Integrator.h"
-#include "Math/IntegratorMultiDim.h"
-#include "Math/AllIntegrationTypes.h"
-#include "Math/Functor.h"
-#include "Math/GaussIntegrator.h"
-#include "time.h"
-
-namespace insane {
-namespace physics {
-
-PhaseSpaceSampler::PhaseSpaceSampler(DiffXSec * xsec)
-{
- fWeight = 1.0;
- fIsModified = true;
- fTotalXSection = 0.0;
- for(int i = 0;i<30;i++) fMCvect[i] = 0;
- fDiffXSec = nullptr;
- fFoam = nullptr;
- fRandomNumberGenerator = nullptr;
-
- fFoamCells = 2000; // default(1000) No of allocated number of cells,
- fFoamSample = 500; // default(200) No. of MC events in the cell MC exploration
- fFoamBins = 8; // default(8) No. of bins in edge-histogram in cell exploration
-
- fFoamOptRej = 1; // default(1) OptRej = 0, weighted; OptRej=1, wt=1 MC events
- fFoamOptDrive = 2; // default(2) Maximum weight reduction, =1 for variance reduction
- fFoamEvPerBin = 50; // default(25) Maximum number of the effective wt=1 events/bin,
- fFoamMaxWtRej = 1.1; //
- fFoamChat = 0; // default(1) 0,1,2 is the ``chat level'' in the standard output
- if (xsec) {
- SetXSec(xsec);
- }
- if (!fRandomNumberGenerator) {
- fRandomNumberGenerator = new TRandom3(0);//gRandom;//RunManager::GetRunManager()->GetRandom(); // Create random number generator
- }
-}
-//________________________________________________________________________________
-
-PhaseSpaceSampler::PhaseSpaceSampler(const PhaseSpaceSampler& rhs) :
- TObject(rhs),
- fWeight(rhs.fWeight) ,
- fIsModified(rhs.fIsModified) ,
- fTotalXSection(rhs.fTotalXSection) ,
- fDiffXSec(rhs.fDiffXSec->Clone()) ,
- fFoam(nullptr) ,
- fRandomNumberGenerator(gRandom) ,
- fFoamCells(rhs.fFoamCells) ,
- fFoamSample(rhs.fFoamSample) ,
- fFoamBins(rhs.fFoamBins) ,
- fFoamOptRej(rhs.fFoamOptRej) ,
- fFoamOptDrive(rhs.fFoamOptDrive) ,
- fFoamEvPerBin(rhs.fFoamEvPerBin) ,
- fFoamChat(rhs.fFoamChat)
-{
- fCrossSection.Clear();
- fXSectionList.AddAll(&(rhs.fXSectionList));
- fCrossSection.Clear();
- fCrossSection.AddAll(&(rhs.fCrossSection));
-}
-//______________________________________________________________________________
-
-PhaseSpaceSampler& PhaseSpaceSampler::operator=(const PhaseSpaceSampler& rhs)
-{
- if (this != &rhs) { // make sure not same object
- TObject::operator=(rhs);
- fWeight = rhs.fWeight ;
- fIsModified = rhs.fIsModified ;
- fTotalXSection = rhs.fTotalXSection ;
- fDiffXSec = rhs.fDiffXSec->Clone();
- fFoam = nullptr;
- fRandomNumberGenerator = gRandom;
- fFoamCells = rhs.fFoamCells ;
- fFoamSample = rhs.fFoamSample ;
- fFoamBins = rhs.fFoamBins ;
- fFoamOptRej = rhs.fFoamOptRej ;
- fFoamOptDrive = rhs.fFoamOptDrive ;
- fFoamEvPerBin = rhs.fFoamEvPerBin ;
- fFoamChat = rhs.fFoamChat ;
- fFoamMaxWtRej = rhs.fFoamMaxWtRej ;
-
- fXSectionList.Clear();
- fXSectionList.AddAll(&(rhs.fXSectionList));
-
- fCrossSection.Clear();
- fCrossSection.AddAll(&(rhs.fCrossSection));
- //fCrossSection; // This exists to store the cross section which doesn't stream by itself
- // // but does when it is in a container. When the sampler is read back from
- // // a file it should initialized with this cross section by calling InitFromDisk
- }
- return *this;
-}
-//______________________________________________________________________________
-
-PhaseSpaceSampler::~PhaseSpaceSampler()
-{
- if(fDiffXSec) delete fDiffXSec;
- fDiffXSec = nullptr;
- //if (fFoam) delete fFoam;
- //fFoam = 0;
- //delete fMCvect;
- //fMCvect = 0;
-}
-//______________________________________________________________________________
-void PhaseSpaceSampler::InitFromDisk() {
- if(gRandom){
- fRandomNumberGenerator = gRandom;//RunManager::GetRunManager()->GetRandom(); // Create random number generator
- }
- fFoam->SetPseRan(fRandomNumberGenerator); // Set random number generator
- if(fCrossSection.GetEntries()>0){
- SetXSec( (DiffXSec*)fCrossSection.At(0) , false);
- } else {
- Error("InitFromDisk","No Cross Section stored in list.");
- }
- std::cout << "going to init cross seciton...\n";
- GetXSec()->InitFromDisk();
-}
-//________________________________________________________________________________
-void PhaseSpaceSampler::SetXSec(DiffXSec * xsec, Bool_t mod) {
- SetModified(mod);
- if(xsec)xsec->SetIncludeJacobian(true);
- fDiffXSec = xsec;
- fCrossSection.Clear("nodelete");
- fCrossSection.Add(fDiffXSec);
- if(fFoam) fFoam->SetRho(xsec);
-}
-//________________________________________________________________________________
-void PhaseSpaceSampler::Print(const Option_t* option) const {
- // option s for short print
- std::cout << " | Weight : " << fWeight << std::endl;
- std::cout << " | Total xsec : " << std::scientific << fTotalXSection << " [nb]" << std::endl;
- std::cout << " | " << std::scientific << fTotalXSection*1.0e-9 << " [b] " << std::endl;
- std::cout << " | " << std::scientific << fTotalXSection*1.0e-33 << " [cm^2] " << std::endl;
- // don't print xsec if short
- if(strcmp("S",option) || strcmp("s",option))if (fDiffXSec) fDiffXSec->Print();
-}
-//________________________________________________________________________________
-void PhaseSpaceSampler::Print(std::ostream& stream) const {
- // option s for short print
- stream << " | Weight : " << fWeight << std::endl;
- stream << " | Total xsec : " << std::scientific << fTotalXSection << " [nb]" << std::endl;
- stream << " | " << std::scientific << fTotalXSection*1.0e-9 << " [b] " << std::endl;
- stream << " | " << std::scientific << fTotalXSection*1.0e-33 << " [cm^2] " << std::endl;
- // don't print xsec if short
- if(fDiffXSec) fDiffXSec->Print(stream);
-}
-//________________________________________________________________________________
-void PhaseSpaceSampler::PrintFoamConfig(std::ostream& c) const {
-}
-//________________________________________________________________________________
-void PhaseSpaceSampler::Refresh(DiffXSec * xsec ) {
- //std::cout << "PhaseSpaceSampler::Refresh" << std::cout ;
- if (xsec) {
- if (fDiffXSec) delete fDiffXSec;
- fDiffXSec = nullptr;
- SetXSec(xsec);
- }
- if (fDiffXSec) {
- fDiffXSec->Refresh();
- if(!fDiffXSec->GetPhaseSpace() ) Error("Refresh()","XSec has Null P.S.");
-
- //fDiffXSec->SetTotalXSec(NormalizePDF());
-
- if (fFoam) delete fFoam;
- fFoam = nullptr;
- fFoam = new TFoam("FoamX"); // Create Simulators
- fFoam->SetkDim( fDiffXSec->NDim()); // No. of dimensions, obligatory!
- fFoam->SetnCells( fFoamCells);
- fFoam->SetnSampl( fFoamSample); // optional
- fFoam->SetChat( fFoamChat); // optional
- fFoam->SetnBin( fFoamBins); // optional
- fFoam->SetOptRej( fFoamOptRej); // optional
- fFoam->SetOptDrive( fFoamOptDrive); // optional
- fFoam->SetMaxWtRej( fFoamMaxWtRej); // optional
- fFoam->SetEvPerBin( fFoamEvPerBin); // optional
- fFoam->SetRho( fDiffXSec); // Set 2-dim distribution,
- //fDiffXSec->Print();
- //std::cout << fDiffXSec << std::endl;
- for(int i = 0; i<fDiffXSec->GetPhaseSpace()->GetNIndependentVariables() ;i++ ) {
- bool is_uniform = fDiffXSec->GetPhaseSpace()->GetIndependentVariable(i)->IsUniform();
- if( is_uniform ){
- fFoam->SetInhiDiv(i,1);
- }
- }
- fFoam->SetPseRan(fRandomNumberGenerator); // Set random number generator
- fFoam->Initialize(); // Initialize simulator, takes a few seconds...
-
- // Note here we are using the weight in calculating the total cross section
- // This could be passed in as luminosity which means that totalXSection is
- // actually a rate.
- fTotalXSection = fWeight*fFoam->GetPrimary()*fDiffXSec->GetPhaseSpace()->GetScaledIntegralJacobian();
-
- fDiffXSec->SetTotalXSec(fTotalXSection);
- } else {
- Error("Refresh()", "Null fDiffXSec!");
- }
- SetModified(false);
-}
-//________________________________________________________________________________
-Double_t PhaseSpaceSampler::NormalizePDF() {
- /// Deprecated
- //fTotalXSection = fFoam->GetPrimary();
- return(fTotalXSection);
-}
-//________________________________________________________________________________
-Double_t PhaseSpaceSampler::CalculateTotalXSection() {
- /// Deprecated
- return(NormalizePDF());
-}
-//________________________________________________________________________________
-Double_t * PhaseSpaceSampler::GenerateEvent() {
- fFoam->MakeEvent(); // generate MC event
- fFoam->GetMCvect(fMCvect); // get generated vector (x,y)
- // for(int i = 0; i<fDiffXSec->GetPhaseSpace()->GetDimension() ;i++ )
- Double_t * vars = fDiffXSec->GetUnnormalizedVariables(fMCvect);
- Double_t * allVars = fDiffXSec->GetDependentVariables(vars);
- // is this needed?
- fDiffXSec->GetPhaseSpace()->SetEventValues(allVars);
- fDiffXSec->DefineEvent(allVars); // Defines all the particle's values
- return(allVars);
-}
-//________________________________________________________________________________
-
-}
-}
diff --git a/src/core/physics/src/PhaseSpaceVariable.cxx b/src/core/physics/src/PhaseSpaceVariable.cxx
deleted file mode 100644
index 3d0302215dc0edb756dcbd80f19cc0d34d4132a5..0000000000000000000000000000000000000000
--- a/src/core/physics/src/PhaseSpaceVariable.cxx
+++ /dev/null
@@ -1,63 +0,0 @@
-#include "PhaseSpaceVariable.h"
-#include <iostream>
-
-
-namespace insane {
-namespace physics {
-
-PhaseSpaceVariable::PhaseSpaceVariable(const char * name, const char * varexp, Double_t min, Double_t max)
- : TNamed(name, varexp)
-{
- fVariableMinima = min;
- fVariableMaxima = max;
- fVariableUnits = "";
- fCurrentValue = min + (max-min)/2.0;
- fiParticle = 0;
- fInverted = false;
- fIsDependent = false;
- fIsModified = true;
-}
-//______________________________________________________________________________
-
-PhaseSpaceVariable::PhaseSpaceVariable(const PhaseSpaceVariable& rhs) : TNamed(rhs),
- fIsModified(rhs.fIsModified),
- fIsDependent(rhs.fIsDependent),
- fInverted(rhs.fInverted),
- fiParticle(rhs.fiParticle),
- fVariableMinima(rhs.fVariableMinima),
- fVariableMaxima(rhs.fVariableMaxima),
- fVariableUnits(rhs.fVariableUnits),
- fCurrentValue(rhs.fCurrentValue)
-{ }
-//______________________________________________________________________________
-
-PhaseSpaceVariable::~PhaseSpaceVariable()
-{ }
-//______________________________________________________________________________
-
-void PhaseSpaceVariable::Print(std::ostream& stream) const
-{
- stream << " "
- << std::setw(10) << std::left << GetName() << " "
- << "part. " << std::left << std::setw(2) << fiParticle << " : "
- << std::setw(10) << std::right << GetMinimum() << " "
- << " < " << std::setw(10) << std::left << GetTitle() << " < "
- << std::setw(10) << std::left << GetMaximum() << " [" << std::setw(5) << GetVariableUnits() << "] with ";
- stream << "value = " << std::setw(4) << fCurrentValue << "\n";
-}
-//______________________________________________________________________________
-
-void PhaseSpaceVariable::Print(const Option_t * opt ) const
-{
- Print(std::cout);
- this->Dump();
- //std::cout << "Variable, "
- // << std::setw(10) << std::left << GetName() << ", "
- // << "particle " << std::left << std::setw(2) << fiParticle << ", limits: "
- // << std::setw(10) << std::right << GetMinimum() << " "
- // << " < " << std::setw(10) << std::left << GetTitle() << " < "
- // << std::setw(10) << std::left << GetMaximum() << " [" << std::setw(10) << GetVariableUnits() << "] with ";
- //std::cout << "current value = " << std::setw(4) << fCurrentValue << "\n";
-}
-//______________________________________________________________________________
-}}
diff --git a/src/physics2/xsections/CMakeLists.txt b/src/physics2/xsections/CMakeLists.txt
index dfd48842277968f2867072035c62035f4d61d485..2512bd6a81d8bdba3f4967b664016f858192066e 100644
--- a/src/physics2/xsections/CMakeLists.txt
+++ b/src/physics2/xsections/CMakeLists.txt
@@ -1,6 +1,6 @@
set(aname xsections)
set(subdirname physics2)
-set(needs_libs InSANEbase2 InSANEphysics InSANExsecs)
+set(needs_libs InSANEbase2 InSANEphysics InSANExsec)
set(libname "InSANE${aname}")
set(dictname "${libname}Dict")
set(lib_LINKDEF "${PROJECT_SOURCE_DIR}/src/${subdirname}/${aname}/include/LinkDef.h")