From 8155a97ea082a575972f12c73eb99bed8f6072a0 Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Sun, 16 Jan 2011 18:26:35 +0000
Subject: [PATCH] * Fix bugs in NPLib for W1 detector + Call to
AddDetector() method was missing
* Finish NPAnalysis for W1 detector
+ Cases with multiplicity greater than one are treated
+ Excitation energy is reconstructed OK without strip
+ Excitation energy taking into account strips should be checked
in more details.
* Cosmetic for screen printing in NPLib classes
---
Inputs/DetectorConfiguration/W1.detector | 12 +-
Inputs/EventGenerator/alpha.source | 4 +-
NPAnalysis/W1/RunToTreat.txt | 3 -
NPAnalysis/W1/configs/ConfigW1.dat | 3 +
NPAnalysis/W1/src/Analysis.cc | 110 +++++-
.../CalibrationManager/CalibrationManager.cxx | 3 +-
NPLib/IORoot/RootInput.cxx | 5 +-
NPLib/Tools/NPEnergyLoss.cxx | 6 +-
NPLib/VDetector/DetectorManager.cxx | 10 +-
NPLib/W1/TW1Physics.cxx | 338 ++++++++++--------
NPLib/W1/TW1Physics.h | 1 +
11 files changed, 301 insertions(+), 194 deletions(-)
create mode 100644 NPAnalysis/W1/configs/ConfigW1.dat
diff --git a/Inputs/DetectorConfiguration/W1.detector b/Inputs/DetectorConfiguration/W1.detector
index 7b2cf7f9c..bf207919a 100644
--- a/Inputs/DetectorConfiguration/W1.detector
+++ b/Inputs/DetectorConfiguration/W1.detector
@@ -1,4 +1,14 @@
-%%%%%%% Telescope 1 %%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GeneralTarget
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 0.001
+ RADIUS= 7.5
+ MATERIAL= CD2
+ X= 0
+ Y= 0
+ Z= 0
+%%%%%%% Detector 1 %%%%%%%
W1
THETA= 0
PHI= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index f840fa369..570339771 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -6,8 +6,8 @@
Isotropic
EnergyLow= 5
EnergyHigh= 5
- HalfOpenAngleMin= 0
- HalfOpenAngleMax= 20
+ HalfOpenAngleMin= 140
+ HalfOpenAngleMax= 160
x0= 0
y0= 0
z0= 0
diff --git a/NPAnalysis/W1/RunToTreat.txt b/NPAnalysis/W1/RunToTreat.txt
index c4154eda3..698f27aa2 100644
--- a/NPAnalysis/W1/RunToTreat.txt
+++ b/NPAnalysis/W1/RunToTreat.txt
@@ -2,6 +2,3 @@ TTreeName
SimulatedTree
RootFileName
../../Outputs/Simulation/myResult.root
-% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T0_B0_E0_S2mm.root
-% ../../Outputs/Simulation/132Sndp_3p3_10MeVA_T0_B1_E0_S05mm.root
-% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T1_B1_E0_S05mm.root
diff --git a/NPAnalysis/W1/configs/ConfigW1.dat b/NPAnalysis/W1/configs/ConfigW1.dat
new file mode 100644
index 000000000..a019d238b
--- /dev/null
+++ b/NPAnalysis/W1/configs/ConfigW1.dat
@@ -0,0 +1,3 @@
+ConfigW1
+ MAX_STRIP_MULTIPLICITY 1
+ STRIP_ENERGY_MATCHING_TOLERANCE 10
diff --git a/NPAnalysis/W1/src/Analysis.cc b/NPAnalysis/W1/src/Analysis.cc
index 152c3fc40..ff556925d 100644
--- a/NPAnalysis/W1/src/Analysis.cc
+++ b/NPAnalysis/W1/src/Analysis.cc
@@ -19,15 +19,18 @@ int main(int argc,char** argv)
RootOutput::getInstance("Analysis/W1_AnalyzedData", "AnalyzedTree");
// Initialize the reaction
+ cout << endl << "/////////// Event generator ///////////" << endl;
NPL::Reaction* myReaction = new Reaction();
myReaction->ReadConfigurationFile(reactionfileName);
// Initialize the detector
+// cout << endl << "/////////// Detector geometry ///////////" << endl;
NPA::DetectorManager* myDetector = new DetectorManager;
myDetector->ReadConfigurationFile(detectorfileName);
// Calculate beam energy at target middle
// Target informations
+ cout << endl;
cout << "/////////// Target information ///////////" << endl;
cout << "Thickness (um): " << myDetector->GetTargetThickness() << endl;
// Get nominal beam energy
@@ -35,7 +38,7 @@ int main(int argc,char** argv)
cout << "Nominal beam energy (MeV): " << BeamEnergyNominal << endl;
// Slow beam at target middle
Double_t BeamEnergy = BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
- cout << "Middle-target beam energy (MeV): " << BeamEnergy << endl << endl;
+ cout << "Middle-target beam energy (MeV): " << BeamEnergy << endl;
// Set energy beam at target middle
myReaction->SetBeamEnergy(BeamEnergy);
@@ -48,7 +51,7 @@ int main(int argc,char** argv)
RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
- // Get GaspardTracker pointer
+ // Get TW1Physics pointer
TW1Physics* W1 = (TW1Physics*) myDetector->m_Detector["W1"];
// Get the input TChain and treat it
@@ -64,10 +67,8 @@ int main(int argc,char** argv)
chain->SetBranchAddress("InteractionCoordinates", &interCoord);
chain->SetBranchStatus("InteractionCoordinates", 0);
- // Analysis is here!
- int nentries = chain->GetEntries();
- cout << "/////////// Loop information ///////////" << endl;
- cout << "Number of entries to be analysed: " << nentries << endl;
+ // random generator
+ TRandom3 *gene = new TRandom3();
// Default initialization
double XTarget = 0;
@@ -75,18 +76,98 @@ int main(int argc,char** argv)
double BeamTheta = 0;
double BeamPhi = 0;
- // random generator
- TRandom3 *gene = new TRandom3();
+ // Analysis is here!
+ int nentries = chain->GetEntries();
+ cout << endl;
+ cout << "/////////// Loop information ///////////" << endl;
+ cout << "Number of entries to be analysed: " << nentries << endl;
+ clock_t begin = clock();
+ clock_t end = begin;
// Loop on all events
for (int i = 0; i < nentries; i ++) {
- if (i%10000 == 0 && i!=0) cout << "\r" << i << " analyzed events" << flush;
+ if (i%10000 == 0 && i!=0) {
+ cout.precision(5);
+ end = clock();
+ double TimeElapsed = (end-begin) / CLOCKS_PER_SEC;
+ double percent = (double)i / nentries;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
+ cout << "\rProgression:" << percent*100 << " % \t | \t Remaining time : ~" << TimeToWait << "s" << flush;
+ }
+ else if (i == nentries-1) cout << "\rProgression:" << " 100%" << endl;
+
+ // Get raw data
chain -> GetEntry(i);
- // Treat Gaspard event
+ // Treat W1 event
myDetector->ClearEventPhysics();
myDetector->BuildPhysicalEvent();
+ // Get Target information from TInitialConditions
+ XTarget = initCond->GetICPositionX(0);
+ YTarget = initCond->GetICPositionY(0);
+ BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
+ BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
+ TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
+
+ // loop on multiplicity event
+ for (int hit = 0; hit < W1->GetEventMultiplicity(); hit++) {
+ // Get c.m. angle
+ double ThetaCM = initCond->GetICEmittedAngleThetaCM(0) * deg;
+
+ // Get exact scattering angle from TInteractionCoordinates object
+ double DetecX = interCoord->GetDetectedPositionX(hit);
+ double DetecY = interCoord->GetDetectedPositionY(hit);
+ double DetecZ = interCoord->GetDetectedPositionZ(hit);
+ TVector3 Detec(DetecX, DetecY, DetecZ);
+
+ // Get interaction position in detector
+ // This takes into account the strips
+ TVector3 A = W1->GetPositionOfInteraction(hit);
+
+ // Hit direction taking into account beam position on target
+ TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
+
+ // Calculate scattering angle w.r.t. optical beam axis (do not take into account beam position on target)
+ double ThetaStrip = ThetaCalculation(A, TVector3(0,0,1));
+ double Theta = ThetaCalculation(Detec, TVector3(0, 0, 1));
+
+ // Calculate scattering angle w.r.t. real beam direction (ideal case)
+// ThetaStrip = ThetaCalculation(HitDirection, BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection);
+ // Calculate scattering angle w.r.t. beam (finite spatial resolution)
+// double resol = 800; // in micrometer
+// double angle = gene->Rndm() * 2*3.14;
+// double r = fabs(gene->Gaus(0, resol)) * micrometer;
+// ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+//
+ // Correct for energy loss in the target
+ double E = W1->GetEnergy(hit);
+ E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip);
+
+ // Calculate excitation energy
+// if (Theta/deg > 150 && Theta/deg < 180) {
+// if (Theta/deg < 60 && ThetaCM/deg < 90) {
+// if (Theta/deg > 35 && Theta/deg < 45 && E/MeV < 17) {
+// if (Theta/deg < 45) {
+// if (E/MeV < 38) { // for (p,t) reaction
+ if (Theta/deg > 30) { // for (d,p) reaction
+ ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
+ Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
+ }
+ else {
+ Ex = -200;
+ ExNoStrips = -200;
+ }
+
+ // Fill output tree
+ RootOutput::getInstance()->GetTree()->Fill();
+ } // end loop on multiplicity event
+ } // end loop on number of events to treat
+
+
+/*
// Get total energy
double E = W1->GetEnergy(0);
@@ -111,13 +192,6 @@ int main(int argc,char** argv)
// This takes into account the strips
A = W1->GetPositionOfInteraction(0);
- // Get beam interaction coordinates on target (from initial condition)
- XTarget = initCond->GetICPositionX(0);
- YTarget = initCond->GetICPositionY(0);
- BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
- BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
- TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
-
// Hit direction taking into account beam position on target
TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
// cout << "A: " << A.X() << " " << A.Y() << " " << A.Z() << endl;
@@ -171,7 +245,7 @@ int main(int argc,char** argv)
// Fill output tree
RootOutput::getInstance()->GetTree()->Fill();
}
-
+*/
// delete singleton classes
RootOutput::getInstance()->Destroy();
RootInput::getInstance()->Destroy();
diff --git a/NPLib/CalibrationManager/CalibrationManager.cxx b/NPLib/CalibrationManager/CalibrationManager.cxx
index 49fe4db9f..564fa8a49 100644
--- a/NPLib/CalibrationManager/CalibrationManager.cxx
+++ b/NPLib/CalibrationManager/CalibrationManager.cxx
@@ -51,7 +51,8 @@ CalibrationManager::CalibrationManager(string configFileName)
ifstream inputConfigFile;
inputConfigFile.open(configFileName.c_str());
- cout << "/////////////////////////////////" << endl;
+ cout << endl;
+ cout << "/////////// Calibration Information ///////////" << endl;
cout << "Getting list of Calibration File" << endl;
if (!inputConfigFile) {
diff --git a/NPLib/IORoot/RootInput.cxx b/NPLib/IORoot/RootInput.cxx
index 4d00e61e5..61f5f2950 100644
--- a/NPLib/IORoot/RootInput.cxx
+++ b/NPLib/IORoot/RootInput.cxx
@@ -65,7 +65,8 @@ RootInput::RootInput(string configFileName)
pRootChain = new TChain();
- cout << "/////////////////////////////////" << endl;
+ cout << endl;
+ cout << "/////////// ROOT Input files ///////////" << endl;
cout << "Initializing input TChain" << endl;
if (!inputConfigFile) {
@@ -107,8 +108,6 @@ RootInput::RootInput(string configFileName)
if (!CheckRootFileName || !CheckTreeName)
cout << "WARNING: Token not found for InputTree Declaration : Input Tree may not be instantiate properly" << endl;
-
- cout << "/////////////////////////////////" << endl;
}
void RootInput::AddFriendChain( string RunToAdd)
diff --git a/NPLib/Tools/NPEnergyLoss.cxx b/NPLib/Tools/NPEnergyLoss.cxx
index 29fd40ba8..3bd1c02ad 100644
--- a/NPLib/Tools/NPEnergyLoss.cxx
+++ b/NPLib/Tools/NPEnergyLoss.cxx
@@ -66,8 +66,9 @@ EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int
string globalPath = getenv("NPTOOL");
string StandardPath = globalPath + "/Inputs/EnergyLoss/" + Path;
- cout << "///////////////////////////////// " << endl ;
- cout << "Initialising an EnergyLoss object " << endl ;
+ cout << endl;
+ cout << "/////////// Energy loss ///////////" << endl ;
+ cout << "Initializing an EnergyLoss object " << endl ;
ifstream TableFile ;
TableFile.open(StandardPath.c_str()) ;
@@ -157,7 +158,6 @@ EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int
}
fInter = new Interpolator( fEnergy , fdEdX_Total ) ;
- cout << "///////////////////////////////// " << endl ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 210ac26ee..41f5f05ab 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -57,8 +57,8 @@ void DetectorManager::ReadConfigurationFile(string Path)
if (ConfigFile.is_open())
{
- cout << "/////////////////////////////" << endl;
- cout << " Configuration file " << Path << " loading " << endl;
+ cout << endl << "/////////// Detector geometry ///////////" << endl;
+ cout << "Configuration file " << Path << " loading " << endl;
Path = StandardPath;
}
@@ -66,7 +66,7 @@ void DetectorManager::ReadConfigurationFile(string Path)
{
ConfigFile.open( Path.c_str() );
if(ConfigFile.is_open()) {
- cout << "/////////////////////////////" << endl;
+ cout << endl << "/////////// Detector geometry ///////////" << endl;
cout << " Configuration file " << Path << " loading " << endl;
}
@@ -157,7 +157,7 @@ void DetectorManager::ReadConfigurationFile(string Path)
////////////////////////////////////////////
else if (LineBuffer.compare(0, 2, "W1") == 0 && W1 == false) {
W1 = true;
- cout << "//////// W1 ////////" << endl << endl;
+ cout << "//////// W1 ////////" << endl;
// Instantiate the new array as a VDetector Object
VDetector* myDetector = new TW1Physics();
@@ -214,7 +214,7 @@ void DetectorManager::ReadConfigurationFile(string Path)
////////////////////////////////////////////
else if (LineBuffer.compare(0, 13, "GeneralTarget") == 0 && GeneralTarget == false) {
GeneralTarget = true ;
- cout << "////////// Target ///////////" << endl << endl;
+ cout << "////////// Target ///////////" << endl;
// jump one line
getline(ConfigFile, LineBuffer);
diff --git a/NPLib/W1/TW1Physics.cxx b/NPLib/W1/TW1Physics.cxx
index f53322512..9be84f76c 100644
--- a/NPLib/W1/TW1Physics.cxx
+++ b/NPLib/W1/TW1Physics.cxx
@@ -94,174 +94,177 @@ void TW1Physics::Clear()
///////////////////////////////////////////////////////////////////////////
void TW1Physics::ReadConfiguration(string Path)
- {
- ifstream ConfigFile ;
- ConfigFile.open(Path.c_str()) ;
- string LineBuffer ;
- string DataBuffer ;
-
- double TLX , BLX , BRX , TRX , TLY , BLY , BRY , TRY , TLZ , BLZ , BRZ , TRZ ;
- double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
- bool check_A = false ;
- bool check_B = false ;
- bool check_C = false ;
- bool check_D = false ;
-
- bool check_Theta = false ;
- bool check_Phi = false ;
- bool check_R = false ;
- bool check_beta = false ;
- bool ReadingStatus = false ;
-
- while (!ConfigFile.eof())
- {
-
- getline(ConfigFile, LineBuffer);
-
- // If line is a Start Up ThinSi bloc, Reading toggle to true
- if (LineBuffer.compare(0, 6, "ThinSi") == 0)
- {
- cout << "Detector found: " << endl ;
- ReadingStatus = true ;
- }
-
- // Else don't toggle to Reading Block Status
- else ReadingStatus = false ;
-
- // Reading Block
- while(ReadingStatus)
- {
- // Pickup Next Word
- ConfigFile >> DataBuffer ;
-
- // Comment Line
- if (DataBuffer.compare(0, 1, "%") == 0) { ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
-
- // Finding another telescope (safety), toggle out
- else if (DataBuffer.compare(0, 6, "ThinSi") == 0) {
- cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
- ReadingStatus = false ;
- }
+{
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer, DataBuffer;
+
+ double TLX, BLX, BRX, TRX, TLY, BLY, BRY, TRY, TLZ, BLZ, BRZ, TRZ;
+ TVector3 A, B, C, D;
+ double Theta = 0, Phi = 0, R = 0, beta_u = 0, beta_v = 0, beta_w = 0;
+ bool check_A = false;
+ bool check_B = false;
+ bool check_C = false;
+ bool check_D = false;
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+ bool check_beta = false;
+ bool ReadingStatus = false;
- //Position method
- else if (DataBuffer.compare(0, 3, "A=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- TLX = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- TLY = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- TLZ = atof(DataBuffer.c_str()) ;
-
- }
-
- else if (DataBuffer.compare(0, 3, "B=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- BLX = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- BLY = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- BLZ = atof(DataBuffer.c_str()) ;
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
- }
+ // If W1 detector found, toggle Reading Block Status
+ if (LineBuffer.compare(0, 2, "W1") == 0) {
+ cout << "Detector found: " << endl;
+ ReadingStatus = true;
+ }
- else if (DataBuffer.compare(0, 3, "C=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- BRX = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- BRY = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- BRZ = atof(DataBuffer.c_str()) ;
+ // else don't toggle to Reading Block Status
+ else ReadingStatus = false;
- }
+ // Reading Block
+ while (ReadingStatus) {
+ // Pickup Next Word
+ ConfigFile >> DataBuffer;
- else if (DataBuffer.compare(0, 3, "D=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- TRX = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- TRY = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- TRZ = atof(DataBuffer.c_str()) ;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
- }
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer.compare(0, 2, "W1") == 0) {
+ cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
+ ReadingStatus = false;
+ }
-
- //Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- }
+ // Position method
+ else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ check_A = true;
+ ConfigFile >> DataBuffer;
+ TLX = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ TLY = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ TLZ = atof(DataBuffer.c_str());
+
+ A = TVector3(TLX, TLY, TLZ);
+ cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+ }
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- }
+ else if (DataBuffer.compare(0, 7, "X16_Y1=") == 0) {
+ check_B = true;
+ ConfigFile >> DataBuffer;
+ BLX = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ BLY = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ BLZ = atof(DataBuffer.c_str());
+
+ B = TVector3(BLX, BLY, BLZ);
+ cout << "X16 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+ }
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
+ else if (DataBuffer.compare(0, 7, "X1_Y16=") == 0) {
+ check_C = true;
+ ConfigFile >> DataBuffer;
+ BRX = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ BRY = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ BRZ = atof(DataBuffer.c_str());
+
+ C = TVector3(BRX, BRY, BRZ);
+ cout << "X1 Y16 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
}
+ else if (DataBuffer.compare(0, 8, "X16_Y16=") == 0) {
+ check_D = true;
+ ConfigFile >> DataBuffer;
+ TRX = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ TRY = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ TRZ = atof(DataBuffer.c_str());
+
+ D = TVector3(TRX, TRY, TRZ);
+ cout << "X16 Y16 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+ }
+
+
+ //Angle method
+ else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
+ check_Theta = true;
+ ConfigFile >> DataBuffer;
+ Theta = atof(DataBuffer.c_str());
+ cout << "Theta: " << Theta << endl;
+ }
+
+ else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
+ check_Phi = true;
+ ConfigFile >> DataBuffer;
+ Phi = atof(DataBuffer.c_str());
+ cout << "Phi: " << Phi << endl;
+ }
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
- check_beta = true;
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- }
+ else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ check_R = true;
+ ConfigFile >> DataBuffer;
+ R = atof(DataBuffer.c_str());
+ cout << "R: " << R << endl;
+ }
+
+ else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
+ check_beta = true;
+ ConfigFile >> DataBuffer;
+ beta_u = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ beta_v = atof(DataBuffer.c_str());
+ ConfigFile >> DataBuffer;
+ beta_w = atof(DataBuffer.c_str());
+ cout << "Beta: " << beta_u << " " << beta_v << " " << beta_w << endl;
+ }
- ///////////////////////////////////////////////////
- // If no Detector Token and no comment, toggle out
- else
- {ReadingStatus = false; cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl ;}
-
- /////////////////////////////////////////////////
- // If All necessary information there, toggle out
-
- if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) )
- {
- ReadingStatus = false;
+ ///////////////////////////////////////////////////
+ // If no Detector Token and no comment, toggle out
+ else {
+ ReadingStatus = false;
+ cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl;
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ((check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta)) {
+ ReadingStatus = false;
- ///Add The previously define telescope
- //With position method
- if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R)) {
- m_NumberOfDetector++;
- }
-
- //with angle method
- else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D)) {
- m_NumberOfDetector++;
- }
-
- // Reinitialisation of Check Boolean
-
- check_A = false ;
- check_B = false ;
- check_C = false ;
- check_D = false ;
-
- check_Theta = false ;
- check_Phi = false ;
- check_R = false ;
- check_beta = false ;
- ReadingStatus = false ;
-
- }
-
+ // Add The previously define telescope
+ // With position method
+ if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R)) {
+ AddDetector(A, B, C, D);
}
- }
+
+ // with angle method
+ else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D)) {
+ AddDetector(Theta, Phi, R, beta_u, beta_v, beta_w);
+ }
+
+ // Reinitialisation of Check Boolean
+ check_A = false;
+ check_B = false;
+ check_C = false;
+ check_D = false;
+
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
+ check_beta = false;
+ ReadingStatus = false;
+ }
+ }
+ }
- InitializeStandardParameter() ;
- ReadAnalysisConfig() ;
+ InitializeStandardParameter();
+ ReadAnalysisConfig();
}
@@ -482,6 +485,23 @@ TVector3 TW1Physics::GetDetectorNormal(int i)
+void TW1Physics::DumpStrippingScheme(int detecNumber)
+{
+ cout << endl << "TW1Physics::DumpStrippingScheme()" << endl;
+ cout << "Detector number " << detecNumber << endl;
+
+ for (int i = 1; i < m_NumberOfStrips+1; i++) { // front part
+ for (int j = 1; j < m_NumberOfStrips+1; j++) { // back part
+ cout << "strips Front, Back: " << i << " " << j << "\t--->\t (X,Y,Z) mm: "
+ << GetStripPositionX(detecNumber, i, j) << "\t"
+ << GetStripPositionY(detecNumber, i, j) << "\t"
+ << GetStripPositionZ(detecNumber, i, j) << endl;
+ }
+ }
+}
+
+
+
///////////////////////////////////////////////////////////////////////////
void TW1Physics::BuildPhysicalEvent()
{
@@ -683,6 +703,8 @@ void TW1Physics::ReadAnalysisConfig()
bool check_mult = false;
bool check_match = false;
+ cout << "\t/////////// Reading ConfigW1.dat file ///////////" << endl;
+
// path to file
string FileName = "./configs/ConfigW1.dat";
@@ -691,10 +713,10 @@ void TW1Physics::ReadAnalysisConfig()
AnalysisConfigFile.open(FileName.c_str());
if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigW1.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ cout << "\tNo ConfigW1.dat found: Default parameter loaded for Analayis " << FileName << endl;
return;
}
- cout << " Loading user parameter for Analysis from ConfigW1.dat " << endl;
+ cout << "\tLoading user parameters from ConfigW1.dat " << endl;
// read analysis config file
string LineBuffer,DataBuffer;
@@ -718,14 +740,14 @@ void TW1Physics::ReadAnalysisConfig()
check_mult = true;
AnalysisConfigFile >> DataBuffer;
m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() );
- cout << "Maximun strip multiplicity= " << m_MaximumStripMultiplicityAllowed << endl;
+ cout << "\tMaximun strip multiplicity= " << m_MaximumStripMultiplicityAllowed << endl;
}
else if (DataBuffer.compare(0, 31, "STRIP_ENERGY_MATCHING_TOLERANCE") == 0) {
check_match = true;
AnalysisConfigFile >> DataBuffer;
m_StripEnergyMatchingTolerance = atoi(DataBuffer.c_str() );
- cout << "Strip energy matching tolerance= " << m_StripEnergyMatchingTolerance << endl;
+ cout << "\tStrip energy matching tolerance= " << m_StripEnergyMatchingTolerance << endl;
}
else if (DataBuffer.compare(0, 2, "W1") == 0) {
@@ -756,7 +778,7 @@ void TW1Physics::ReadAnalysisConfig()
}
else {
- cout << "Warning: detector type for W1 unknown!" << endl;
+ cout << "\tWarning: detector type for W1 unknown!" << endl;
}
}
@@ -785,7 +807,7 @@ void TW1Physics::ReadAnalysisConfig()
}
else {
- cout << "Warning: don't know what to do (lost in translation)" << endl;
+ cout << "\tWarning: don't know what to do (lost in translation)" << endl;
}
}
else {
diff --git a/NPLib/W1/TW1Physics.h b/NPLib/W1/TW1Physics.h
index 03fc0b7a3..28d62c54a 100644
--- a/NPLib/W1/TW1Physics.h
+++ b/NPLib/W1/TW1Physics.h
@@ -176,6 +176,7 @@ class TW1Physics : public TObject, public NPA::VDetector
double GetStripPositionZ(int N, int X, int Y) {return m_StripPositionZ[N-1][X-1][Y-1];};
TVector3 GetPositionOfInteraction(int i);
TVector3 GetDetectorNormal(int i);
+ void DumpStrippingScheme(int detecNumber);
private: // Geometry and strip number
--
GitLab