diff --git a/Projects/T40/Analysis.cxx b/Projects/T40/Analysis.cxx
index 77d2e4fbf08f3d11518c6ed910a4926eaed8ba0c..c64c6e8644ab594664930a49046a35455641f155 100644
--- a/Projects/T40/Analysis.cxx
+++ b/Projects/T40/Analysis.cxx
@@ -19,11 +19,13 @@
* *
* *
*****************************************************************************/
-#include<cassert>
-#include<iostream>
+#include <cassert>
+#include <iostream>
using namespace std;
+#include "TSystem.h"
+
#include "Analysis.h"
#include "NPAnalysisFactory.h"
#include "NPDetectorManager.h"
@@ -65,9 +67,9 @@ double calculate_fit_slope(int len, double* Aw_X, double* Aw_Z, double& R2)
/// TODO::: R2 doesn't seem to make sense... look into it!
return slope;
-} }
-
+}
+}
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -75,21 +77,76 @@ Analysis::Analysis(){
Analysis::~Analysis(){
}
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::Init(){
+namespace {
- TH = (TTiaraHyballPhysics*) m_DetectorManager -> GetDetector("HyballWedge");
- TB = (TTiaraBarrelPhysics*) m_DetectorManager -> GetDetector("Tiara");
- TF = (TFPDTamuPhysics*) m_DetectorManager -> GetDetector("FPDTamu");
- TG = (TGeTAMUPhysics*) m_DetectorManager -> GetDetector("GeTAMU");
+NPL::EnergyLoss check_energy_loss(const string& particle, const string& target){
+ string globalPath = getenv("NPTOOL");
+ string standardPath = globalPath + "/Inputs/EnergyLoss/";
+ string srimPath = particle+"_"+target+".SRIM";
+ string geantPath = particle+"_"+target+".G4table";
+
+ string pathList[4] = {
+ srimPath, // SRIM in PWD
+ standardPath + srimPath, // SRIM in standard location
+ geantPath, // GEANT in PWD
+ standardPath + geantPath // GEANT in standard location
+ };
+
+ int i=0;
+ for(;i<4;++i){
+ std::ifstream ifs(pathList[i].c_str());
+ if(ifs.good()) { break; }
+ }
+ int nstep;
+ string type;
+ switch(i){
+ case 0:
+ cout << "Using SRIM energy loss file in \"Projects/T40\"\n";
+ nstep = 10;
+ type = "SRIM";
+ break;
+ case 1:
+ cout << "Using SRIM energy loss file in \"Inputs/EnergyLoss\"\n";
+ nstep = 10;
+ type = "SRIM";
+ break;
+ case 2:
+ cout << "Using GEANT energy loss file in \"Projects/T40\"\n";
+ nstep = 100;
+ type = "G4table";
+ break;
+ case 3:
+ cout << "Using GEANT energy loss file in \"Inputs/EnergyLoss\"\n";
+ nstep = 100;
+ type = "G4table";
+ break;
+ default:
+ cerr << "FATAL ERROR:: Energy Loss file not found for particle: \"" << particle
+ << "\", and material: \"" << target << "\" (neither SRIM nor GEANT)\n";
+ exit(1);
+ break;
+ }
+
+ return NPL::EnergyLoss(pathList[i], type.c_str(), nstep);
+} }
+
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init(){
+ TH = (TTiaraHyballPhysics*) m_DetectorManager -> GetDetector("HyballWedge"); assert(TH);
+ TB = (TTiaraBarrelPhysics*) m_DetectorManager -> GetDetector("Tiara"); assert(TB);
+ TF = (TFPDTamuPhysics*) m_DetectorManager -> GetDetector("FPDTamu"); assert(TF);
+ TG = (TGeTAMUPhysics*) m_DetectorManager -> GetDetector("GeTAMU"); assert(TG);
+ MDM = (TMDMPhysics*) m_DetectorManager -> GetDetector("MDM"); assert(MDM);
- // get reaction information
+ // get reaction information
myReaction = new NPL::Reaction();
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
OriginalBeamEnergy = myReaction->GetBeamEnergy();
cout << "Original Beam energy (entrance of target): " << OriginalBeamEnergy << endl ;
-
+ if(MDM) { MDM->SetReaction(myReaction); }
+
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -109,20 +166,27 @@ void Analysis::Init(){
*/
//Copied from Momo's Slack 170222.
+//GAC 171003 - check for existence of SRIM file, use that if available; if not,
+//default to G4table
+//
string light=NPL::ChangeNameToG4Standard(myReaction->GetNucleus3().GetName());
string beam=NPL::ChangeNameToG4Standard(myReaction->GetNucleus1().GetName());
- LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".SRIM","SRIM",10 );
+
+ LightTarget = check_energy_loss(light,TargetMaterial);
+
//by Shuya 170505
//Note when you analyze the triple alpha calibration run, use He4_Al and He4_Si
- LightAl = NPL::EnergyLoss(light+"_Al.SRIM","SRIM",10);
+ LightAl = check_energy_loss(light,"Al");
+
//LightAl = NPL::EnergyLoss("He4_Al.SRIM","SRIM",10);
- LightSi = NPL::EnergyLoss(light+"_Si.SRIM","SRIM",10);
+ LightSi = check_energy_loss(light,"Si");
//LightSi = NPL::EnergyLoss("He4_Si.SRIM","SRIM",10);
//by Shuya 170530
//LightCBacking = NPL::EnergyLoss(light+"_C.SRIM","SRIM",10);
- BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".SRIM","SRIM",10);
+ BeamTarget = check_energy_loss(beam,TargetMaterial);
+
FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness*0.5, 0);
myReaction->SetBeamEnergy(FinalBeamEnergy);
cout << "Final Beam energy (middle of target): " << FinalBeamEnergy << endl;
@@ -141,8 +205,8 @@ void Analysis::Init(){
//Original_ELab=0;
//Original_ThetaLab=0;
- XTarget =-1.026126;
- YTarget =-2.3589;
+ XTarget =0;
+ YTarget =0;
BeamDirection = TVector3(0,0,1);
InitOutputBranch();
InitInputBranch();
@@ -228,15 +292,9 @@ void Analysis::TreatEvent(){
ThetaTHSurface = 0;
ThetaNormalTarget = 0;
if(XTarget>-1000 && YTarget>-1000){
- //TVector3 BeamImpact(XTarget,YTarget,0);
- //by Shuya 171020 from 22Ne(d,d) IB7,8 5 degrees adjustment
- //TVector3 BeamImpact(0.173098, 2.67341, 4.07383);
- //by Shuya 171020 from 22Ne(d,d) IB7,8 5 degrees adjustment, but X, Y=0
- TVector3 BeamImpact(XTarget, YTarget, 4.07383);
+ TVector3 BeamImpact(XTarget,YTarget,0);
- //by Shuya 171218 (because of T40 meeting's discussion)
- //TVector3 HitDirection = TH -> GetRandomisedPositionOfInteraction(countTiaraHyball) - BeamImpact ;
- TVector3 HitDirection = TH -> GetPositionOfInteraction(countTiaraHyball) - BeamImpact ;
+ TVector3 HitDirection = TH -> GetRandomisedPositionOfInteraction(countTiaraHyball) - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
ThetaTHSurface = HitDirection.Angle(TVector3(0,0,-1)); // vector Normal on Hyball
@@ -245,8 +303,14 @@ void Analysis::TreatEvent(){
//by Shuya 171019
PhiLab = HitDirection.Phi();
PhiLab = PhiLab/(TMath::Pi())*180.0;
- //by Shuya 171208
- PhiLab_Hyball = PhiLab;
+
+ // GAC 171020
+ {
+ TVector3 v;
+ v.SetMagThetaPhi(1,ThetaLab*deg,PhiLab*deg);
+ ThetaXLab = atan(v.X()/v.Z()) / deg;
+ ThetaYLab = atan(v.Y()/v.z()) / deg;
+ }
}
else{
BeamDirection = TVector3(-1000,-1000,-1000);
@@ -257,9 +321,6 @@ void Analysis::TreatEvent(){
/////////////////////////////
// Part 2 : Impact Energy
Energy = ELab = 0;
-//by Shuya 171206
- ELab_Hyball = 0;
-
Si_E_TH = TH->Strip_E[countTiaraHyball];
Energy = Si_E_TH; // calibration for hyball is in MeV
// Correct for energy loss using the thickness of the target and the dead layer
@@ -281,9 +342,6 @@ void Analysis::TreatEvent(){
ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
ThetaLab=ThetaLab/deg;
-//by Shuya 171206
- ELab_Hyball = ELab;
- ThetaLab_Hyball = ThetaLab;
//by Shuya 170703
Ex_Hyball = Ex;
@@ -303,16 +361,9 @@ void Analysis::TreatEvent(){
ThetaTBSurface = 0;
ThetaNormalTarget = 0;
if(XTarget>-1000 && YTarget>-1000){
- //TVector3 BeamImpact(XTarget,YTarget,0);
- //by Shuya 171020 from 22Ne(d,d) IB7,8 5 degrees adjustment
- //TVector3 BeamImpact(0.173098, 2.67341, 4.07383);
- //by Shuya 171020 from 22Ne(d,d) IB7,8 5 degrees adjustment, but X, Y=0
- TVector3 BeamImpact(XTarget, YTarget, 4.07383);
-
- //by Shuya 171218 (because of T40 meeting's discussion)
- //TVector3 HitDirection = TB -> GetRandomisedPositionOfInteraction(countTiaraBarrel) - BeamImpact ;
- TVector3 HitDirection = TB -> GetPositionOfInteraction(countTiaraBarrel) - BeamImpact ;
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = TB -> GetRandomisedPositionOfInteraction(countTiaraBarrel) - BeamImpact ;
//Angle of emission wrt to beam
ThetaLab = HitDirection.Angle( BeamDirection );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
@@ -325,8 +376,14 @@ void Analysis::TreatEvent(){
//by Shuya 171019
PhiLab = HitDirection.Phi();
PhiLab = PhiLab/(TMath::Pi())*180.0;
- //by Shuya 171208
- PhiLab_Barrel = PhiLab;
+ //GAC 171020
+ {
+ TVector3 v;
+ v.SetMagThetaPhi(1,ThetaLab*deg,PhiLab*deg);
+ ThetaXLab = atan(v.X()/v.Z()) / deg;
+ ThetaYLab = atan(v.Y()/v.z()) / deg;
+ }
+
}
else{
BeamDirection = TVector3(-1000,-1000,-1000);
@@ -337,9 +394,6 @@ void Analysis::TreatEvent(){
/////////////////////////////
// Part 2 : Impact Energy
Energy = ELab = 0;
-//by Shuya 171206
- ELab_Barrel = 0;
-
Si_E_InnerTB = TB->Strip_E[countTiaraBarrel];
Energy = Si_E_InnerTB*keV;// calibration for barrel is in keV
@@ -350,14 +404,6 @@ void Analysis::TreatEvent(){
Energy = Si_E_InnerTB*keV + Si_E_OuterTB*keV;
}
- //by Shuya 171208. If you need E+dE for Barrel.
- for(unsigned int countTiaraOuterBarrel = 0 ; countTiaraOuterBarrel < TB->Outer_Strip_E.size() ; countTiaraOuterBarrel++){
- if(TB->Outer_Detector_N[countTiaraOuterBarrel]==TB->Detector_N[countTiaraBarrel] && TB->Outer_Strip_E[countTiaraOuterBarrel]>0){
- Si_E_OuterTB = TB->Outer_Strip_E[countTiaraOuterBarrel];
- Energy = Si_E_InnerTB*keV + Si_E_OuterTB*keV;
- }
- }
-
// Evaluate energy using the thickness, Target and Si dead layer Correction
if(Energy > 0){
ELab = LightSi.EvaluateInitialEnergy( Energy ,0.3*micrometer, ThetaTBSurface);
@@ -379,10 +425,6 @@ void Analysis::TreatEvent(){
ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
ThetaLab=ThetaLab/deg;
-//by Shuya 171206
- ELab_Barrel = ELab;
- ThetaLab_Barrel = ThetaLab;
-
/////////////////////////////
// Part 5 : Implementing randomised position impact matrix for both the entire Barrel (all 8 strips) and each strip making up the octagonal Barrel individually
TVector3 BarrelRandomImpactPosition = TB -> GetRandomisedPositionOfInteraction(countTiaraBarrel);
@@ -516,6 +558,16 @@ void Analysis::TreatEvent(){
if(detNumber >=0 && detNumber< 4) {
Aw_X[detNumber] = TF->AWirePositionX[iw];
Aw_Z[detNumber] = TF->AWirePositionZ[iw];
+
+ // Fill MDM class with FPD data
+ // Only do this if it's not there already
+ // (e.g. data files, not similation)
+ if(MDM && RootInput::getInstance()->GetChain()->GetBranch("MDM") == 0) {
+ MDM->DetectorNumber.push_back(detNumber);
+ MDM->Xpos.push_back(Aw_X[detNumber]);
+ MDM->Ypos.push_back(0);
+ MDM->Zpos.push_back(Aw_Z[detNumber]);
+ }
}
else {
cerr << "WARNING:: Wire number not bewtween 0 and 4!\n";
@@ -596,6 +648,19 @@ void Analysis::TreatEvent(){
else
TacSiGeOR = -999;
+ /////////////////////////////////////////////////////////////
+ // MDM RECONSTRUCTION OF TARGET PARAMETERS //////////////////
+ /////////////////////////////////////////////////////////////
+ if(MDM){
+ if(Ex != -1000) MDM->SetEx4(Ex);
+ else MDM->SetEx4(0);
+ MDM->SetLightParticleAngles(ThetaLab, PhiLab);
+
+ // do target parameter minimization
+ MDM->MinimizeTarget();
+ }
+
+
//by Shuya 170524
//RunNumber = RootInput::getInstance()->GetChain()->GetFileNumber() + 1;
if(currentfilename != RootInput::getInstance()->GetChain()->GetCurrentFile()->GetName())
@@ -620,12 +685,6 @@ void Analysis::ReInitValue(){
//Silicon
Ex = -1000 ;
ELab = -1000;
-//by Shuya 171206
- ELab_Hyball = -1000;
- ELab_Barrel = -1000;
- ThetaLab_Hyball = -1000;
- ThetaLab_Barrel = -1000;
-
ThetaLab = -1000;
ThetaCM = -1000;
LightParticleDetected = false ;
@@ -634,9 +693,7 @@ void Analysis::ReInitValue(){
Ex_Barrel = -1000 ;
//by Shuya 171019
PhiLab = -1000;
-//by Shuya 171208
- PhiLab_Hyball = -1000;
- PhiLab_Barrel = -1000;
+ ThetaXLab = ThetaYLab = -1000;
//Simu
//Original_ELab = -1000;
@@ -714,20 +771,13 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex_Barrel",&Ex_Barrel,"Ex_Barrel/D");
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
-//by Shuya 171206
- RootOutput::getInstance()->GetTree()->Branch("ELab_Hyball",&ELab_Hyball,"ELab_Hyball/D");
- RootOutput::getInstance()->GetTree()->Branch("ELab_Barrel",&ELab_Barrel,"ELab_Barrel/D");
-
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
-//by Shuya 171206
- RootOutput::getInstance()->GetTree()->Branch("ThetaLab_Hyball",&ThetaLab_Hyball,"ThetaLab_Hyball/D");
- RootOutput::getInstance()->GetTree()->Branch("ThetaLab_Barrel",&ThetaLab_Barrel,"ThetaLab_Barrel/D");
-
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
//by Shuya 171019
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
- RootOutput::getInstance()->GetTree()->Branch("PhiLab_Hyball",&PhiLab_Hyball,"PhiLab_Hyball/D");
- RootOutput::getInstance()->GetTree()->Branch("PhiLab_Barrel",&PhiLab_Barrel,"PhiLab_Barrel/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaXLab",&ThetaXLab,"ThetaXLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaYLab",&ThetaYLab,"ThetaYLab/D");
+
RootOutput::getInstance()->GetTree()->Branch("TiaraImpactMatrixX",&TiaraIMX,"TiaraImpactMatrixX/D");
RootOutput::getInstance()->GetTree()->Branch("TiaraImpactMatrixY",&TiaraIMY,"TiaraImpactMatrixY/D");
@@ -792,8 +842,7 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("TacSiMicro_dE",&TacSiMicro_dE,"TacSiMicro_dE/D");
RootOutput::getInstance()->GetTree()->Branch("TacSiPlastLeft",&TacSiPlastLeft,"TacSiPlastLeft/D");
- RootOutput::getInstance()->GetTree()->Branch("TacSiPlastRight",&TacSiPlastRight,"TacSiPlastRight/D");
-
+ RootOutput::getInstance()->GetTree()->Branch("TacSiPlastRight",&TacSiPlastRight,"TacSiPlastRight/D");
// Other
RootOutput::getInstance()->GetTree()->Branch("RunNumber", &RunNumber,"RunNumber/I");
diff --git a/Projects/T40/Analysis.h b/Projects/T40/Analysis.h
index aae93b6640c5ae6dc882405d4c0846c5e3936c00..938e3fffcaf5d740e2ebfe36e8a2d2f1334d30a0 100644
--- a/Projects/T40/Analysis.h
+++ b/Projects/T40/Analysis.h
@@ -27,7 +27,7 @@
#include "TTiaraBarrelPhysics.h"
#include "TFPDTamuPhysics.h"
#include "TGeTAMUPhysics.h"
-//#include "TMDMPhysics.h"
+#include "TMDMPhysics.h"
#include "TInitialConditions.h"
#include "NPEnergyLoss.h"
#include "NPReaction.h"
@@ -61,15 +61,7 @@ class Analysis: public NPL::VAnalysis{
double Ex_Barrel;
double ELab;
-//by Shuya 171206
- double ELab_Hyball;
- double ELab_Barrel;
-
double ThetaLab;
-//by Shuya 171206
- double ThetaLab_Hyball;
- double ThetaLab_Barrel;
-
double ThetaCM;
double TiaraIMX;
double TiaraIMY;
@@ -79,10 +71,6 @@ class Analysis: public NPL::VAnalysis{
bool LightParticleDetected;
//by Shuya 171019
double PhiLab;
-//by Shuya 171206
- double PhiLab_Hyball;
- double PhiLab_Barrel;
-
// GAC 171020
double ThetaXLab;
double ThetaYLab;
@@ -103,8 +91,8 @@ class Analysis: public NPL::VAnalysis{
TTiaraBarrelPhysics* TB;
TFPDTamuPhysics* TF;
TGeTAMUPhysics* TG;
- //TMDMPhysics* MDM;
- //bool MDM_in_file; //!
+ TMDMPhysics* MDM;
+ bool MDM_in_file; //!
TRandom *Rand ;
double ThetaNormalTarget ;