diff --git a/HOWTO b/HOWTO
index 4eb80a248efe7c050a16c7b25a356c00ab39990a..a4a3a4a4adb950cd36aad437ea88907b3ca7a352 100644
--- a/HOWTO
+++ b/HOWTO
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/INSTALL b/INSTALL
index eb99292569a1add4c4a22b90fb3b9ddc7ca08b08..f712ec176f846e67e5bf78f5fdb42feb40503747 100644
--- a/INSTALL
+++ b/INSTALL
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector
index bee456b3d8bd5104bcb017232ee3a0cc5bc03615..2cace566d03d4546598209b313a2100903f7e81a 100644
--- a/Inputs/DetectorConfiguration/Riken_65mm.detector
+++ b/Inputs/DetectorConfiguration/Riken_65mm.detector
@@ -17,12 +17,12 @@ GeneralTarget
Z= 0
CryoTarget
- THICKNESS= 3000
+ THICKNESS= 3
RADIUS= 45
TEMPERATURE= 26
PRESSURE= 1
MATERIAL= D2
- WINDOWSTHICKNESS= 15
+ WINDOWSTHICKNESS= 0
WINDOWSMATERIAL= Mylar
X= 0
Y= 0
@@ -149,6 +149,28 @@ D= 17.61 -9.85 104.11
ScintillatorPlastic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Plastic
+ X= 0
+ Y= -15
+ Z= 318
+ Shape= Square
+ Height= 70
+ Width= 40
+ Thickness= 20
+ Scintillator= BC400
+ LeadThickness= 0
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Plastic
+ X= 0
+ Y= -15
+ Z= 343
+ Thickness= 30
+ Shape= Square
+ Height= 70
+ Width= 40
+ Scintillator= BC400
+ LeadThickness= 2
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Plastic
THETA= 0
PHI= 0
R= 318
@@ -157,7 +179,7 @@ Plastic
Scintillator= BC400
LeadThickness= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Plastic
+%Plastic
THETA= 0
PHI= 0
R= 343
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 8567e32e3e2aad009e6f1f14e403149c23216e86..9b71622344f5ef1e607ba911c66cf308cd444332 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -7,7 +7,7 @@ TransfertToResonance
Target= 2H
Light= 3He
Heavy= 10He
- ExcitationEnergy= 0.0
+ ExcitationEnergy= 1.0
BeamEnergy= 550
BeamEnergySpread= 0
SigmaThetaX= 0.6921330164
@@ -19,7 +19,7 @@ TransfertToResonance
ResonanceDecayA= 8
CrossSectionPath= 11Li(d,3He)10He.txt
ShootLight= 1
- ShootHeavy= 1
+ ShootHeavy= 0
ShootDecayProduct= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/11Li.beam b/Inputs/EventGenerator/11Li.beam
index 5a06c27989f836364fc4ee579144f4e531185245..c32a2702bc904b44898f856c982c038e5fdda9a2 100644
--- a/Inputs/EventGenerator/11Li.beam
+++ b/Inputs/EventGenerator/11Li.beam
@@ -6,8 +6,8 @@
Beam
ParticleZ= 3
ParticleA= 11
- BeamEnergy= 550
- BeamEnergySpread= 10
+ BeamEnergy= 500
+ BeamEnergySpread= 50
SigmaX= 6.232
SigmaY= 9.069
SigmaThetaX= 0.6921330164
diff --git a/Inputs/EventGenerator/34Si.reaction b/Inputs/EventGenerator/34Si.reaction
index edecfe2aa56337481e8871aedf2e9e364b0ce9cc..41a35ec8e0a24a92782a3bf5202d1700d4c628ac 100644
--- a/Inputs/EventGenerator/34Si.reaction
+++ b/Inputs/EventGenerator/34Si.reaction
@@ -1,21 +1,37 @@
+Transfert
+ Beam= 19C
+ Target= 208Pb
+ Light= 208Pb
+ Heavy= 19C
+ ExcitationEnergy= 5.9
+ BeamEnergy= 722
+ BeamEnergySpread= 0
+ SigmaThetaX= 0.6921330164
+ SigmaPhiY= 0.963142053
+ SigmaX= 6.232
+ SigmaY= 9.069
+ CrossSectionPath= flat.txt
+ ShootLight= 0
+ ShootHeavy= 1
+
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Beam energy given in MeV ; Excitation in MeV 748 612
-Transfert
- Beam= 34Si
- Target= 2H
- Light= 1H
- Heavy= 35Si
- ExcitationEnergy= 2.0
- BeamEnergy= 748
- BeamEnergySpread= 0
- SigmaX= 0.5
- SigmaY= 0.5
- SigmaThetaX= 0
- SigmaPhiY= 0
- BeamSpreadX= 0
- CrossSectionPath= ni69_g7_01.n
- ShootLight= 1
- ShootHeavy= 0
+%Transfert
+% Beam= 34Si
+% Target= 2H
+% Light= 1H
+% Heavy= 35Si
+% ExcitationEnergy= 2.0
+% BeamEnergy= 748
+% BeamEnergySpread= 0
+% SigmaX= 0.5
+% SigmaY= 0.5
+% SigmaThetaX= 0
+% SigmaPhiY= 0
+% CrossSectionPath= flat.txt
+% ShootLight= 1
+% ShootHeavy= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index 8fa0b8698412b53f993a685d44e12862aabcb8ff..22957a57d0c16cf6c69813e071d726986534ad59 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 21.6
+ EnergyHigh= 300
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/Inputs/EventGenerator/60Fe.reaction b/Inputs/EventGenerator/60Fe.reaction
index 467bd51d9509590ead464e0f56adc36d461b14b9..ea0874438d3a70b841bdbd8852d269b7e4d3e949 100644
--- a/Inputs/EventGenerator/60Fe.reaction
+++ b/Inputs/EventGenerator/60Fe.reaction
@@ -10,10 +10,10 @@ Transfert
ExcitationEnergy= 2.0
BeamEnergy= 1200
BeamEnergySpread= 0
- SigmaX= 1
- SigmaY= 1
- SigmaThetaX= 2
- SigmaPhiY= 2
+ SigmaX= 3
+ SigmaY= 3
+ SigmaThetaX= 0.5
+ SigmaPhiY= 0.5
CrossSectionPath= ni69_g7_01.n
ShootLight= 1
ShootHeavy= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index fd5f19e5259addc2b6bbccaeb2e2bcb97b08b3e1..343524fe6bf0f84f05874a7e5090b5475ba93b40 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 21.6
+ EnergyHigh= 25
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index 336fb8c7142848ae3f053a591911a9684948e3f2..47e09a3cb2df16b98579a2ef3f7a8f651ef9c268 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh
@@ -31,6 +31,7 @@
#include "RootInput.h"
#include "RootOutput.h"
#include "TInitialConditions.h"
+#include "TInteractionCoordinates.h"
#include "TMust2Physics.h"
#include "TSSSDPhysics.h"
#include "TPlasticPhysics.h"
@@ -74,7 +75,7 @@ namespace GRAPH
{
// Declare your Spectra here:
- TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
+ TH1F* myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 100 , -40 , 40 ) ;
TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
@@ -110,58 +111,56 @@ namespace ENERGYLOSS
// 3He Energy Loss
EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" ,
"G4Table",
- 10000 );
+ 1000,
+ 3 );
EnergyLoss He3TargetGaz = EnergyLoss ( "He3_D2.G4table" ,
"G4Table",
- 10000 );
+ 1000,
+ 3 );
- EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
- "LISE",
- 10000 ,
- 1 ,
- 3 );
-
- EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
- "LISE",
- 10000 ,
- 1 ,
+ EnergyLoss He3StripAl = EnergyLoss ( "He3_Aluminium.G4table" ,
+ "G4Table",
+ 1000,
3 );
-
-
-// // 3He Energy Loss
// EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
-// 10000 ,
-// 1 ,
-// 3 );
+// "LISE",
+// 10000 ,
+// 1 ,
+// 3 );
//
-// EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
-// 10000 ,
-// 1 ,
-// 3 );
+// EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
+// "LISE",
+// 10000 ,
+// 1 ,
+// 3 );
//
-// EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
-// 10000 ,
-// 1 ,
-// 3 );
-//
-// EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
-// 10000 ,
-// 1 ,
-// 3 );
+// EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
+// "LISE",
+// 10000 ,
+// 1 ,
+// 3 );
+
+
+
+ EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
+ "LISE",
+ 100 ,
+ 3 ,
+ 1 );
// proton Energy Loss
EnergyLoss protonTargetWind = EnergyLoss ( "proton_Mylar.txt" ,
"LISE",
- 1000 ,
+ 100 ,
1 ,
1 );
EnergyLoss protonTargetGaz = EnergyLoss ( "proton_D2gaz_1b_26K.txt" ,
"LISE",
- 1000 ,
+ 100 ,
1 ,
1 );
diff --git a/NPAnalysis/10He_Riken/macro/CrossSection.c b/NPAnalysis/10He_Riken/macro/CrossSection.c
index 5640efb446f5464a98882b5a40ea015ef2745171..56e071c820cdedaccfe71c6e3a1a7894070474cf 100644
--- a/NPAnalysis/10He_Riken/macro/CrossSection.c
+++ b/NPAnalysis/10He_Riken/macro/CrossSection.c
@@ -11,16 +11,8 @@
double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
double RadToDeg = 180./Pi ; // 360/2Pi = 180/Pi
-TFile *file0 = TFile::Open("./Result/myResult.root");
-
- cEA = new TCanvas("cEA","Kinematic Line" ,100,100,900,900);
- hEA->Draw("COLZ");
- cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600);
- hEx->Draw();
-
- cEHexa = new TCanvas("cEHexa","Hexaneutron bound Energy" ,100,100,600,600);
- hEHexa->Draw();
-
+TFile *file0 = TFile::Open("./Result/thetaCM.root");
+
cCM = new TCanvas("cCm" , "Cross Section (CM)" , 100 , 100 , 900, 900) ;
hThetaCM->Draw();
@@ -79,8 +71,7 @@ TFile *file0 = TFile::Open("./Result/myResult.root");
//Normalisation:
- //Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
- Int_t Maximum_Bin = 3 ;
+ Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ;
Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ;
Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ;
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index ed1a372c4b79fa6f9a43e552785d8beb9bd4eb98..bebe502b09a365729d464274121a6f8c7965f70a 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -34,13 +34,14 @@ int main(int argc,char** argv)
myDetector -> ReadConfigurationFile(detectorfileName) ;
// Instantiate the Calibration Manger using a file
- CalibrationManager* myCalibration = new CalibrationManager(calibrationfileName) ;
+ CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
// Attach more branch to the output
double ELab[2], ExcitationEnergy[2] ;
double ThetaLab[2] , ThetaCM[2] ;
double X[2] , Y[2] ;
+ double EventWeight=0;
// Exitation Energy
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",ExcitationEnergy,"ExcitationEnergy[2]/D") ;
@@ -52,38 +53,76 @@ int main(int argc,char** argv)
RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") ;
RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]") ;
+ // For phasespace event
+ RootOutput::getInstance()->GetTree()->Branch("EventWeight",&EventWeight,"EventWeight/D") ;
// Get the formed Chained Tree and Treat it
TChain* Chain = RootInput:: getInstance() -> GetChain() ;
// Open the ThinSi Branch
+ Chain->SetBranchStatus("EventWeight",true);
+ Chain->SetBranchAddress("EventWeight" ,&EventWeight );
+
Chain->SetBranchStatus("InitialConditions",true) ;
Chain->SetBranchStatus("fIC_*",true) ;
TInitialConditions* Init = new TInitialConditions();
Chain->SetBranchAddress("InitialConditions" ,&Init );
- double XTarget=0 ; double YTarget=0; double BeamTheta = 0 ; double BeamPhi = 0 ; double E=-1000;
+ TInteractionCoordinates* ICoord = new TInteractionCoordinates();
+ Chain->SetBranchAddress("InteractionCoordinates" ,&ICoord );
+
+
+
+ double XTarget=0 ; double YTarget=0; double ZTarget = 0 ; double BeamTheta = 0 ; double BeamPhi = 0 ; double E=-1000;
// Get Detector Pointer:
TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ;
TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ;
TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ;
- cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
-
+
+
+RootOutput::getInstance()->GetList()->Add(myHist1D);
+
+ TCutG *cut3He_MUST2 = new TCutG("Cut3HeMUST2",11);
+ cut3He_MUST2->SetPoint(0,-1.49426,24.2781);
+ cut3He_MUST2->SetPoint(1,15.3161,14.3151);
+ cut3He_MUST2->SetPoint(2,47.069,7.6732);
+ cut3He_MUST2->SetPoint(3,110.575,4.35222);
+ cut3He_MUST2->SetPoint(4,308.563,2.25477);
+ cut3He_MUST2->SetPoint(5,310.431,1.20604);
+ cut3He_MUST2->SetPoint(6,232.917,1.11864);
+ cut3He_MUST2->SetPoint(7,89.0948,2.42955);
+ cut3He_MUST2->SetPoint(8,32.1264,6.71186);
+ cut3He_MUST2->SetPoint(9,1.30747,11.9555);
+ cut3He_MUST2->SetPoint(10,-1.49426,24.2781);
+
+ TCutG *cut3He_M2_SSSD = new TCutG("Cut3HeM2SSSD",11);
+ cut3He_M2_SSSD->SetPoint(0,7.44252,1.45432);
+ cut3He_M2_SSSD->SetPoint(1,40.6322,0.684984);
+ cut3He_M2_SSSD->SetPoint(2,95.9483,0.406912);
+ cut3He_M2_SSSD->SetPoint(3,218.649,0.216896);
+ cut3He_M2_SSSD->SetPoint(4,362.471,0.101033);
+ cut3He_M2_SSSD->SetPoint(5,355.431,0.0315148);
+ cut3He_M2_SSSD->SetPoint(6,100.977,0.0315148);
+ cut3He_M2_SSSD->SetPoint(7,25.546,0.314221);
+ cut3He_M2_SSSD->SetPoint(8,3.41954,0.995498);
+ cut3He_M2_SSSD->SetPoint(9,3.41954,1.40797);
+ cut3He_M2_SSSD->SetPoint(10,7.44252,1.45432);
+
+ cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
int i ,N=Chain -> GetEntries();
cout << " Number of Event to be treated : " << N << endl ;
-
clock_t begin=clock();
clock_t end=begin;
- for ( i = 0 ; i < N ; i ++ )
+ for ( i = 0 ; i < N ; i ++ )
{
// Clear local branch
for(int hh = 0 ; hh <2 ; hh++)
{
- ELab[hh] = -1 ; ExcitationEnergy[hh] = -1 ; ThetaLab[hh] = -1 ;
- X[hh] = -1000 ; Y[hh] = -1000 ; ThetaCM[hh] = -1 ;
+ ELab[hh] = -10000 ; ExcitationEnergy[hh] = -10000 ; ThetaLab[hh] = -10000 ;
+ X[hh] = -10000 ; Y[hh] = -10000 ; ThetaCM[hh] = -10000 ;
}
// Minimum code
@@ -107,65 +146,73 @@ int main(int argc,char** argv)
// Build the new one
myDetector -> BuildPhysicalEvent() ;
////
-
-
-
-
// Target (from initial condition)
XTarget = Init->GetICPositionX(0);
YTarget = Init->GetICPositionY(0);
+ ZTarget = Init->GetICPositionZ(0);
// XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
// YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
- ////
+ ////
+
// Must 2 And ThinSi //
- //for(int hit = 0; hit < M2 -> GetEventMultiplicity() ; hit ++)
for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
{
- ELab[hit] = -1 ; ThetaLab[hit] = -1;
// Get Hit Direction
TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0);
-
// Angle between beam and particle
- ThetaLab[hit] = ThetaCalculation ( HitDirection , BeamDirection ) ;
+ ThetaLab[hit] = ThetaCalculation ( HitDirection , BeamDirection ) ;
+
+ double Xint = ICoord-> GetDetectedPositionX(hit) - M2 -> GetPositionOfInteraction(hit).X();
+ double Yint = ICoord-> GetDetectedPositionY(hit) - M2 -> GetPositionOfInteraction(hit).Y();
+ double Zint = ICoord-> GetDetectedPositionZ(hit) - M2 -> GetPositionOfInteraction(hit).Z();
+
+
+// cout << Xint << " " << Yint << " " << Zint << endl;
+
// Angle between particule and z axis (target Normal)
double ThetaN = ThetaCalculation ( HitDirection , TVector3(0,0,1) ) ;
- // ANgle between particule and Must2 Si surface
+ // Angle between particule and Must2 Si surface
double ThetaMM2Surface = ThetaCalculation ( HitDirection , M2 -> GetTelescopeNormal(hit) );
if(M2 -> GetPositionOfInteraction(hit).Z() > 0)
{
- if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0)
+ if( M2 -> CsI_E[hit] < 0 && M2 -> Si_E[hit] > 0 )
{
ELab[hit] = M2 -> Si_E[hit] ;
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 2*0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
-
-// ELab[hit]= He3StripSi.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
-// 20*micrometer , // Target Thickness at 0 degree
-// ThetaMM2Surface );
-
- if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit];
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
+ if(ThinSi -> Energy.size() > 0)
+ {
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
+ ELab[hit] += ThinSi-> Energy[hit];
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
+ }
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
+ ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0*micrometer , // Target Thickness at 0 degree
+ ThetaN );
- ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
+ ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 3*micrometer , // Target Thickness at 0 degree
+ ThetaN );
ThetaCM[hit] = He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ;
+
+ if(ThinSi -> Energy.size() > 0 )
+ if(cut3He_M2_SSSD->IsInside(ThinSi -> Energy[hit], M2 -> Si_E[hit]) )
+ myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+
X[hit] = HitDirection . X();
Y[hit] = HitDirection . Y();
ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -176,36 +223,44 @@ int main(int argc,char** argv)
ELab[hit]= M2 ->CsI_E[hit] ;
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 3*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
+ ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 3*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
ELab[hit]+= M2 ->Si_E[hit];
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
- if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit];
+ if(ThinSi -> Energy.size() > 0)
+ {
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
+ ELab[hit] += ThinSi-> Energy[hit];
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
+ }
+
+ ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0*micrometer , // Target Thickness at 0 degree
+ ThetaN );
-
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
-
- ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
-
+ ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 3*micrometer , // Target Thickness at 0 degree
+ ThetaN );
+
ThetaCM[hit]= He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
- ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ;
+ ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit]) ;
+
+ if( cut3He_MUST2->IsInside(M2 -> Si_E[hit], M2 -> CsI_E[hit]) )
+ myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+
X[hit] = HitDirection . X();
Y[hit] = HitDirection . Y();
ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -216,47 +271,15 @@ int main(int argc,char** argv)
}
/*else if(M2 -> GetPositionOfInteraction(hit).Z()<0)
- {
-
- if(ELab[hit]>-1000 )
- {
- if(ELab[hit]>18)
- {
- ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
- }
-
- ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= protonTargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
-
- ELab[hit]= protonTargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
- ThetaCM[hit] = myReaction -> EnergyLabToThetaCM( ELab[hit] , 1 ) /deg ;
- ExcitationEnergy[hit] = myReaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ;
- X[hit] = HitDirection . X();
- Y[hit] = HitDirection . Y();
- ThetaLab[hit] = ThetaLab[hit] / deg ;
- }
-
- else {ExcitationEnergy[hit]=-100 ; X[hit] = -100 ; Y[hit] = -100 ;ThetaLab[hit]=-100; ThetaCM[hit]=-100 ;}
-
- } */
-
+ {} */
}
-
RootOutput::getInstance()->GetTree()->Fill() ;
}
cout << " A total of " << i << " event has been annalysed " << endl ;
- cout << endl << " ///////////////////////////////////// "<< endl<< endl ;
+ cout << endl << "/////////////////////////////////"<< endl<< endl ;
+ myHist1D->Write();
RootOutput::getInstance()->Destroy();
return 0 ;
}
@@ -264,6 +287,6 @@ int main(int argc,char** argv)
double ThetaCalculation (TVector3 A , TVector3 B)
{
double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
- return Theta*rad ;
+ return Theta ;
}
diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc
index 8162feff58c6a792513052f34f76118579d00000..6d46e9a27a3575d3593fc0e69e68e19848ca9989 100644
--- a/NPAnalysis/Gaspard/src/Analysis.cc
+++ b/NPAnalysis/Gaspard/src/Analysis.cc
@@ -41,8 +41,14 @@ int main(int argc,char** argv)
// Set energy beam at target middle
myReaction->SetBeamEnergy(BeamEnergy);
- // Print target thickness
- cout << myDetector->GetTargetThickness() << endl;
+ // nominal beam energy
+ Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV;
+ cout << BeamEnergyNominal << endl;
+ // slow beam at target middle
+ Double_t BeamEnergy = BeamEnergyNominal - BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
+ cout << BeamEnergy << endl;
+ // set energy beam at target middle
+ myReaction->SetBeamEnergy(BeamEnergy);
// Attach more branch to the output
double Ex = 0 ; double ExNoStrips = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
diff --git a/NPAnalysis/e530/Makefile b/NPAnalysis/e530/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..f5bfeb69ffd25b28c9ccb8a81eebb1b66664dfa2
--- /dev/null
+++ b/NPAnalysis/e530/Makefile
@@ -0,0 +1,6 @@
+
+Analyse:
+ make -C ./src
+
+clean:
+ make clean -C ./src
diff --git a/NPAnalysis/e530/RunToTreat.txt b/NPAnalysis/e530/RunToTreat.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c56d1ce44cd471573db652b261e0f3be3adc4cd8
--- /dev/null
+++ b/NPAnalysis/e530/RunToTreat.txt
@@ -0,0 +1,7 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ %../../Outputs/Simulation/3He_source.root
+ %../../Outputs/Simulation/alpha_source.root
+ %../../Outputs/Simulation/REALIST_DATA.root
+ ../../Outputs/Simulation/mySimul.root
diff --git a/NPAnalysis/e530/include/ObjectManager.hh b/NPAnalysis/e530/include/ObjectManager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..036e1a34f65df03c02e04424554c6185972b913c
--- /dev/null
+++ b/NPAnalysis/e530/include/ObjectManager.hh
@@ -0,0 +1,194 @@
+// You can use this file to declare your spectra, file, energy loss , ... and whatever you want.
+// This way you can remove all unnecessary declaration in the main programm.
+// In order to help debugging and organizing we use Name Space.
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
+
+// NPA
+#include "DetectorManager.h"
+
+// STL C++
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+#include <time.h>
+
+// ROOT
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TLeaf.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+// NPL
+#include "TPlasticData.h"
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TInitialConditions.h"
+#include "TMust2Physics.h"
+#include "TSSSDPhysics.h"
+#include "TPlasticPhysics.h"
+
+// Use CLHEP System of unit and Physical Constant
+#include "CLHEP/Units/GlobalSystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+
+// ----------------------------------------------------------------------------------------------
+double ThetaCalculation (TVector3 A , TVector3 B) ;
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
+namespace VARIABLE
+ {
+ // Declare your Variable here:
+
+ double X1,Y1,Z1 ;
+ int N1,N2 = 0 ;
+ bool check= false ;
+
+ // A Usefull Simple Random Generator
+ TRandom Rand;
+ }
+
+using namespace VARIABLE ;
+// ----------------------------------------------------------------------------------------------
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------GRAPH------------------------------------------------------
+#include <TObject.h>
+#include <TH1.h>
+#include <TH1F.h>
+#include <TH2.h>
+#include <TH2F.h>
+#include <TGraph2D.h>
+
+namespace GRAPH
+ {
+ // Declare your Spectra here:
+
+ TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
+
+ TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
+
+ }
+
+using namespace GRAPH ;
+// --------------------------------------------------------------------------------------------
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------CUT------------------------------------------------------
+#include <TCutG.h>
+namespace CUT
+ {
+ // Declare your Cut here:
+
+ }
+
+using namespace CUT ;
+// --------------------------------------------------------------------------------------------
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------ENERGY LOSS----------------------------------------------
+#include "NPEnergyLoss.h"
+using namespace NPL ;
+namespace ENERGYLOSS
+ {
+
+
+ // 3He Energy Loss
+ /*
+ EnergyLoss He3TargetWind = EnergyLoss ("He3_Mylar.G4table",
+ "G4Table" ,
+ 10000 );
+
+ EnergyLoss He3TargetGaz = EnergyLoss ("He3_D2.G4table" ,
+ "G4Table" ,
+ 10000 );
+
+ EnergyLoss He3StripAl = EnergyLoss ("3He_Al.txt" ,
+ "LISE" ,
+ 10000 ,
+ 1 ,
+ 3 );
+
+ EnergyLoss He3StripSi = EnergyLoss ("3He_Si.txt" ,
+ "LISE" ,
+ 10000 ,
+ 1 ,
+ 3 );
+
+ */
+
+ // // 3He Energy Loss
+ // EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
+ // 10000 ,
+ // 1 ,
+ // 3 );
+ //
+ // EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
+ // 10000 ,
+ // 1 ,
+ // 3 );
+ //
+ // EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
+ // 10000 ,
+ // 1 ,
+ // 3 );
+ //
+ // EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
+ // 10000 ,
+ // 1 ,
+ // 3 );
+
+
+ // proton Energy Loss
+ EnergyLoss protonTargetWind = EnergyLoss ("proton_Mylar.txt",
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 1 );
+
+ EnergyLoss protonTargetGaz = EnergyLoss ("proton_D2gaz_1b_26K.txt" ,
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 1 );
+
+ EnergyLoss protonStripAl = EnergyLoss ("proton_Al.txt" ,
+ "LISE" ,
+ 100 ,
+ 1 ,
+ 1 );
+
+
+ }
+
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------Random Engine----------------------------------------------
+#include "TRandom3.h"
+namespace RANDOMENGINE
+ {
+
+ TRandom3 RandomEngine = TRandom3();
+
+ }
+
+using namespace RANDOMENGINE ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
diff --git a/NPAnalysis/e530/macro/CrossSection.c b/NPAnalysis/e530/macro/CrossSection.c
new file mode 100644
index 0000000000000000000000000000000000000000..5640efb446f5464a98882b5a40ea015ef2745171
--- /dev/null
+++ b/NPAnalysis/e530/macro/CrossSection.c
@@ -0,0 +1,121 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+
+ double Pi = 3.141592654 ;
+
+///Usefull myAngleInRad * RadToDeg = myAngle In deg :p
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+ double RadToDeg = 180./Pi ; // 360/2Pi = 180/Pi
+
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ cEA = new TCanvas("cEA","Kinematic Line" ,100,100,900,900);
+ hEA->Draw("COLZ");
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600);
+ hEx->Draw();
+
+ cEHexa = new TCanvas("cEHexa","Hexaneutron bound Energy" ,100,100,600,600);
+ hEHexa->Draw();
+
+ cCM = new TCanvas("cCm" , "Cross Section (CM)" , 100 , 100 , 900, 900) ;
+ hThetaCM->Draw();
+
+ vector<double> relative_error;
+ double NumberOfHitX=0;
+ ifstream efficFile;
+ efficFile.open("./17cm.efficiency");
+ string EffString, ThetaString;
+ double xxx = hThetaLabCM->GetNbinsY();
+
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ // Calculate Relative error
+ NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx );
+ if(NumberOfHitX!=0)relative_error.push_back( 1. / ( sqrt(NumberOfHitX) ) );
+ else {relative_error.push_back(0) ;}
+ }
+
+ for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++)
+ {
+ //Efficacity correction
+ efficFile >> ThetaString >> EffString;
+ double Eff = atoi(EffString.c_str());
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ if (Eff!=0)
+ {
+ double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff) ;
+ hThetaLabCM->SetCellContent( i, j, pipo );
+ }
+ else hThetaLabCM->SetCellContent( i, j, 0 );
+ }
+ }
+
+ TH1D* hCrossSection = new TH1D();
+ hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ;
+
+ for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++)
+ {
+ // Calculate Error
+ hCrossSection->SetBinContent(i,
+ NumberOfHitX=hCrossSection->GetBinContent(i)/(sin(hCrossSection->GetBinCenter(i)*DegToRad) ));
+ double error=NumberOfHitX*relative_error[i];
+ hCrossSection->SetBinError(i,error);
+ }
+
+ hCrossSection->Draw();
+
+ string Path="../Data/CrossSection/11Li(d,3He)10He.txt";
+ ifstream CrossSectionFile;
+ CrossSectionFile.open(Path.c_str());
+ if( CrossSectionFile.is_open() )
+ cout << " file " << Path << " loading " << endl;
+ else{
+ cout << " Error, no file " << Path << " found" << endl; return;}
+
+
+ //Normalisation:
+ //Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
+ Int_t Maximum_Bin = 3 ;
+ Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ;
+ Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ;
+ Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor = 0 ;
+ Double_t ThetaVal = 0 ;
+ Double_t CSVal = 0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+ CrossSectionFile >> theta >> CS ;
+
+ ThetaVal=atof(theta.c_str()) ;
+
+ if(ThetaVal>(Maximum_theta-Bin_Width/1000) && ThetaVal<(Maximum_theta+Bin_Width/1000) )
+ {
+ CSVal=atof(CS.c_str()) ;
+ Normalisation_Factor = Maximum/CSVal;
+ }
+ }
+ CrossSectionFile.close();
+ CrossSectionFile.open(Path.c_str());
+
+ //Reading Cross Section
+
+ ThetaVal=0;
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor ) ;
+ }
+
+}
diff --git a/NPAnalysis/e530/macro/TimeOfFlight.c b/NPAnalysis/e530/macro/TimeOfFlight.c
new file mode 100644
index 0000000000000000000000000000000000000000..916932fc04c540832f14c0c4cf14b4889d8972be
--- /dev/null
+++ b/NPAnalysis/e530/macro/TimeOfFlight.c
@@ -0,0 +1,14 @@
+{
+
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ //gStyle->SetPalette(1);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800);
+ hTOF1234->Draw("COLZ");
+
+}
diff --git a/NPAnalysis/e530/macro/affich.c b/NPAnalysis/e530/macro/affich.c
new file mode 100644
index 0000000000000000000000000000000000000000..2c8537ea62be3b91ec797d490b630dd8f4ecc6eb
--- /dev/null
+++ b/NPAnalysis/e530/macro/affich.c
@@ -0,0 +1,291 @@
+{
+
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ //gStyle->SetPalette(51,0);
+ gStyle->SetPalette(1);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+/* TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800);
+ cTOF->Divide(3,2);
+
+ cTOF->cd(1);
+ hTOF1->Draw("COLZ");
+
+ cTOF->cd(2);
+ hTOF2->Draw("COLZ");
+
+ cTOF->cd(3);
+ hTOF3->Draw("COLZ");
+
+ cTOF->cd(4);
+ hTOF4->Draw("COLZ");
+
+ cTOF->cd(5);
+ hTOF5->Draw("COLZ"); */
+
+ double Pi = 3.141592654 ;
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+
+/* vector<double> relative_error;
+ double NumberOfHitX=0;
+ cTheta2D = new TCanvas("cTheta2d","Cross Section" ,100,100,900,900);
+ hThetaLabCM->Draw("COLZ");*/
+ cTheta = new TCanvas("cTheta","Cross Section" ,100,100,900,900);
+ hTheta->Draw();
+/* ifstream efficFile;
+
+ efficFile.open("./efficiency/20cm2.efficiency");
+ string EffString, ThetaString;
+ double xxx = hThetaLabCM->GetNbinsY();
+
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ // Calculate Relative error
+ NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx );
+ if(NumberOfHitX!=0)relative_error.push_back( ( sqrt(NumberOfHitX) ) / (NumberOfHitX) );
+ else {relative_error.push_back(0) ;}
+ }
+
+ for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++)
+ {
+ //Efficacity correction
+ efficFile >> ThetaString >> EffString;
+ double Eff = atoi(EffString.c_str());
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ if (Eff!=0)
+ {
+ double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff/100) ;
+ hThetaLabCM->SetCellContent( i, j, pipo );
+ }
+ else hThetaLabCM->SetCellContent( i, j, 0 );
+ }
+ }
+
+
+ TH1D* hCrossSection = new TH1D();
+ hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ;
+
+
+ for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++)
+ {
+ // Calculate Error
+ NumberOfHitX=hCrossSection->GetBinContent(i);
+ double error=NumberOfHitX*relative_error[i];
+ hCrossSection->SetBinError(i,error);
+ }
+
+ hCrossSection->Draw();
+
+/* ofstream efficFile;
+ efficFile.open("20cm.efficiency");
+
+ for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++)
+ {
+
+ Double_t Flux = 100000*sin(hTheta->GetBinCenter(i)*DegToRad)*2*DegToRad/2;
+
+ if(Flux!=0)
+ {
+ Double_t Efficiency = ( ( Flux - hTheta->GetBinContent(i) ) / Flux ) * 100 ;
+
+ hTheta->SetBinContent(i, 100-Efficiency ) ;
+
+ efficFile << hTheta->GetBinCenter(i) << " " << 100-Efficiency << endl;
+ }
+
+ else hTheta->SetBinContent(i, 0 );
+ }
+
+
+ //hTheta->Draw();*/
+
+ //efficiency correction and error bar calculation
+/* ifstream efficFile;
+ efficFile.open("./efficiency/20cmCM.efficiency");
+ string EffString, ThetaString;
+ double Eff, relative_error, error, toto ;
+
+ for(Int_t i = 1 ; i < hThetaCM->GetNbinsX() ; i++)
+ {
+
+ efficFile >> ThetaString >> EffString;
+ Eff = atoi(EffString.c_str());
+ if(hThetaCM->GetBinContent(i)!=0)relative_error = sqrt(hThetaCM->GetBinContent(i))/(hThetaCM->GetBinContent(i));
+
+ toto=hThetaCM->GetBinContent(i)*(100-Eff) ;
+ hThetaCM->SetBinContent(i, toto) ;
+
+ error=hThetaCM->GetBinContent(i)*relative_error;
+ hThetaCM->SetBinError(i,error);
+ }*/
+
+
+
+
+ //hThetaCM->Draw();
+
+/* string Path="/home/Adrien/Desktop/geant/Simulation/Data/CrossSection/11Li(d,3He)10He.txt";
+ ifstream CrossSectionFile;
+ CrossSectionFile.open(Path.c_str());
+ if( CrossSectionFile.is_open() )
+ cout << " file " << Path << " loading " << endl;
+ else{
+ cout << " Error, no file " << Path << " found" << endl; return;}
+
+
+ //Normalisation:
+ Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
+ Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ;
+ Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ;
+ Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor=0;
+
+
+
+ Double_t ThetaVal=0 ;
+ Double_t CSVal=0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+ CrossSectionFile >> theta >> CS ;
+
+ ThetaVal=atof(theta.c_str()) ;
+
+ if(ThetaVal>(Maximum_theta-Bin_Width/4) && ThetaVal<(Maximum_theta+Bin_Width/4) )
+ {
+ CSVal=atof(CS.c_str()) ;
+ Normalisation_Factor = Maximum/CSVal;
+ }
+ }
+ CrossSectionFile.close();
+ CrossSectionFile.open(Path.c_str());
+
+ //Reading Cross Section
+
+ ThetaVal=0;
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ;
+ }
+
+
+
+/*
+ //Normalisation:
+ Double_t Maximum_Bin = hTheta->GetMaximumBin() ;
+ Double_t Maximum = hTheta->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor = Maximum/3.;
+
+
+
+ //Reading Cross Section
+ Double_t ThetaVal=0 ;
+ Double_t CSVal=0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+
+ ThetaVal=0.000479597*pow(ThetaVal,3)-0.0670771*pow(ThetaVal,2)+3.08149*ThetaVal-0.339958 ;
+
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ;
+ }
+// cTheta->SetLogy(1);
+
+
+
+
+ /*
+ cKine = new TCanvas("cKine","Kinematics lines" ,10,10,800,800);
+ hKine->Draw("COLZ");*/
+
+/* cKineInit = new TCanvas("cKineInit","Initial Kinematics lines" ,100,100,600,600);
+ hKineInit->Draw("COLZ");
+
+ cEDE = new TCanvas("EDE","EDE add Strip" ,100,100,600,600);
+ hEDEAddStrip->Draw("COLZ");
+/* cG = new TCanvas("cG","Strip position",500,100,800,600);
+ cG->Divide(3,2);
+
+ cG->cd(1);
+ Agraph2D->Draw("P0") ;
+ cG->cd(2);
+ Bgraph2D->Draw("P0") ;
+ cG->cd(3);
+ Cgraph2D->Draw("P0") ;
+ cG->cd(4);
+ Dgraph2D->Draw("P0") ;
+ cG->cd(5);
+ Egraph2D->Draw("P0") ;
+
+ cG2 = new TCanvas("cG2","all Strip position",500,100,800,600);
+ TOTgraph2D->Draw("P0") ;
+ */
+ cH = new TCanvas("cH","Hit density",500,100,1000,800);
+ cH->Divide(3,2);
+
+ cH->cd(1);
+ hHIT4->Draw("COLZ");
+ cH->cd(2);
+ hHIT2->Draw("COLZ");
+ cH->cd(3);
+ hHIT5->Draw("COLZ");
+ cH->cd(4);
+ hHIT1->Draw("COLZ");
+ cH->cd(5);
+ hHIT3->Draw("COLZ");
+ cH->cd(6);
+ hThetaHeavy->Draw();
+
+/* cH->cd(5);
+ hXZ->Draw("COLZ");
+ cH->cd(6);
+ hXY->Draw("COLZ");
+
+ cEDE = new TCanvas("cEDE","EDE indentification Spectra",500,100,1000,800);
+ cEDE->Divide(3,2);
+
+ cEDE->cd(1);
+ hEDE1->Draw("COLZ");
+ cEDE->cd(2);
+ hEDE2->Draw("COLZ");
+ cEDE->cd(3);
+ hEDE3->Draw("COLZ");
+ cEDE->cd(4);
+ hEDE4->Draw("COLZ");
+ cEDE->cd(5);
+ hEDE5->Draw("COLZ");
+ cEDE->cd(6);
+ TH2F* hEDET= new TH2F("hEDET","",4000,-1,600, 400, -1, 25) ;
+ hEDET->Add(hEDE1);
+ hEDET->Add(hEDE2);
+ hEDET->Add(hEDE3);
+ hEDET->Add(hEDE4);
+ hEDET->Draw("CLOZ");*/
+
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600);
+ hEx->Draw();
+ /*cE = new TCanvas("cE","Light Energy" ,500,100,1000,800);
+ cE->Divide(2);
+ cE->cd(1);
+ hE1234->Draw();
+ cE->cd(2);
+ hE5->Draw();*/
+}
diff --git a/NPAnalysis/e530/macro/efficiency.c b/NPAnalysis/e530/macro/efficiency.c
new file mode 100644
index 0000000000000000000000000000000000000000..9a04af185f75853d0bce5bae2009cd8cd01cae1a
--- /dev/null
+++ b/NPAnalysis/e530/macro/efficiency.c
@@ -0,0 +1,58 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ double Pi = 3.141592654 ;
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+
+
+ ofstream efficFile;
+ efficFile.open("12cm.efficiency");
+
+ for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++)
+ {
+ Double_t Flux = 1000000./2. * sin(hTheta->GetBinCenter(i)*DegToRad) * (2*DegToRad) ;
+
+ if(Flux!=0)
+ {
+ Double_t Efficiency = 100*( (hTheta->GetBinContent(i)) / Flux );
+ hTheta->SetBinContent(i, Efficiency ) ;
+ }
+
+ else hTheta->SetBinContent(i, 0 );
+
+ efficFile << hTheta->GetBinCenter(i) << " " << Efficiency << endl;
+ }
+
+ cEff = new TCanvas("cEff","Efficiency" ,100,100,600,600);
+ hTheta->Draw();
+
+ cHit = new TCanvas("cHit","Hit" ,100,100,600,600);
+ hXY->Draw("COLZ");
+
+ cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600);
+ cHit2->Divide(3,2);
+
+ cHit2->cd(1);
+ hXY1->Draw("COLZ");
+
+ cHit2->cd(2);
+ hXY2->Draw("COLZ");
+
+ cHit2->cd(3);
+ hXY3->Draw("COLZ");
+
+ cHit2->cd(4);
+ hXY4->Draw("COLZ");
+
+ cHit2->cd(5);
+ hXY5->Draw("COLZ");
+
+ cHit2->cd(6);
+ hXY6->Draw("COLZ");
+
+}
diff --git a/NPAnalysis/e530/macro/hit.c b/NPAnalysis/e530/macro/hit.c
new file mode 100644
index 0000000000000000000000000000000000000000..50db19da5e7c54b1cba977ffe34232fb328dcdd0
--- /dev/null
+++ b/NPAnalysis/e530/macro/hit.c
@@ -0,0 +1,44 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ cEff = new TCanvas("cEff","Theta Distribution" ,100,100,600,600);
+ hTheta->Draw();
+
+
+ cCM = new TCanvas("cCm" , "ThetaCM" , 100 , 100 , 600, 600) ;
+ hThetaCM->Draw();
+
+ cHit = new TCanvas("cHit","Hit" ,100,100,600,600);
+ hXY->Draw("COLZ");
+
+
+ cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600);
+ cHit2->Divide(3,2);
+
+ cHit2->cd(1);
+ hXY1->Draw("COLZ");
+
+ cHit2->cd(2);
+ hXY2->Draw("COLZ");
+
+ cHit2->cd(3);
+ hXY3->Draw("COLZ");
+
+ cHit2->cd(4);
+ hXY4->Draw("COLZ");
+
+ cHit2->cd(5);
+ hXY5->Draw("COLZ");
+
+ cHit2->cd(6);
+ hXY6->Draw("COLZ");
+
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,300,300);
+ hEx->Draw();
+
+}
diff --git a/NPAnalysis/e530/pipo.txt b/NPAnalysis/e530/pipo.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a6ffecd1bfe9a754e329aa65d0f5fedaec9284f8
--- /dev/null
+++ b/NPAnalysis/e530/pipo.txt
@@ -0,0 +1,2 @@
+CalibrationFilePath
+ /home/Adrien/Desktop/NPTool/NPTool.dev.CalibMust2/NPAnalysis/e530/calibrrrPipo.txt
diff --git a/NPAnalysis/e530/src/Analysis.cc b/NPAnalysis/e530/src/Analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2561d14e48382f007406e0f1693d80867e5e0aaa
--- /dev/null
+++ b/NPAnalysis/e530/src/Analysis.cc
@@ -0,0 +1,111 @@
+#include "ObjectManager.hh"
+
+using namespace std;
+
+int main(int argc,char** argv)
+{
+
+ if(argc!=5)
+ {
+ cout <<
+ "you need to specify both a Reaction file and a Detector file such as : "<< endl;
+ cout <<
+ "Analysis myReaction.reaction myDetector.detector runToRead.run" << endl;
+ return 0;
+ }
+
+ string reactionfileName = argv[1] ;
+ string detectorfileName = argv[2] ;
+ string calibrationfileName = argv[3] ;
+ string runToTreatFileName = argv[4] ;
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
+ // First of All instantiate RootInput and Output
+ // Detector will be attached later
+ RootInput:: getInstance(runToTreatFileName) ;
+ RootOutput::getInstance("Analysis/60Fe_AnalyzedData", "AnalyzedTree") ;
+
+ // Instantiate some Reaction
+ NPL::Reaction* e530Reaction = new Reaction ;
+ e530Reaction -> ReadConfigurationFile(reactionfileName) ;
+
+ // Instantiate the Calibration Manger using a file (WARNING:prior to the detector instantiation)
+ CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
+
+ // Instantiate the detector using a file
+ NPA::DetectorManager* myDetector = new DetectorManager ;
+ myDetector -> ReadConfigurationFile(detectorfileName) ;
+
+ // Ask the detector manager to load the parameter added by the detector in the calibrationfileName
+ myCalibration->LoadParameterFromFile() ;
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // Attach more branch to the output
+ double ELab=0,ThetaLab=0,ExcitationEnergy=0;
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&ExcitationEnergy,"ExcitationEnergy/D") ;
+ // Get the formed Chained Tree and Treat it
+ TChain* Chain = RootInput:: getInstance() -> GetChain() ;
+
+ TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ;
+ cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
+
+ int i ,N=Chain -> GetEntries();
+
+ //N = 1000;
+ cout << " Number of Event to be treated : " << N << endl ;
+
+ clock_t begin=clock();
+ clock_t end=begin;
+ for ( i = 0 ; i < N ; i ++ )
+ {
+ // Minimum code
+ if( i%10000 == 0 && i!=0) {
+ cout.precision(5);
+ end=clock();
+ double TimeElapsed = (end-begin)/CLOCKS_PER_SEC;
+ double percent = (double)i/N ;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed ;
+ cout << "\r Progression:" << percent*100
+ << " % \t | \t Remaining time : ~"
+ << TimeToWait <<"s"<< flush;
+ }
+
+ else if (i==N-1) cout << "\r Progression:"
+ << " 100% " <<endl;
+
+ Chain -> GetEntry(i);
+ // Clear Previous Event
+ myDetector -> ClearEventPhysics() ;
+ // Build the new one
+ myDetector -> BuildPhysicalEvent() ;
+ ////
+
+
+ // Must 2
+ for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
+ {
+ ELab = -1 ; ThetaLab = -1;
+ // Get Hit Direction
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(0,0,-40);
+ // Angle between beam and particle
+ ThetaLab = ThetaCalculation ( HitDirection , TVector3(0,0,1) ) ;
+ ELab = M2 -> Si_E[hit] + M2 -> SiLi_E[hit] ;
+ }
+
+ RootOutput::getInstance()->GetTree()->Fill() ;
+ }
+
+ cout << " A total of " << i << " event has been annalysed " << endl ;
+ cout << endl << " ///////////////////////////////////// "<< endl<< endl ;
+ RootOutput::getInstance()->Destroy();
+ return 0 ;
+}
+
+double ThetaCalculation (TVector3 A , TVector3 B)
+{
+ double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
+ return Theta*rad ;
+}
+
diff --git a/NPAnalysis/e530/src/GNUmakefile b/NPAnalysis/e530/src/GNUmakefile
new file mode 100644
index 0000000000000000000000000000000000000000..28c404622cec4d40b2e93aa0009fc4fc4fa97d26
--- /dev/null
+++ b/NPAnalysis/e530/src/GNUmakefile
@@ -0,0 +1,42 @@
+###Make file made by Adrien MATTA/ Institut de Physique Nucleaire d'Orsay IPNO###
+# Made to compile the ROOT Analyser for MUST2 experiment
+
+CPP=g++
+EXEC=Analysis
+
+# local includes
+NPAINCLUDES = ../include
+
+# ROOT includes
+CXXFLAGS += `root-config --cflags`
+
+# CLHEP includes
+CXXFLAGS += -I$(CLHEP_INCLUDE_DIR)
+CXXFLAGS += -I$(NPAINCLUDES)
+CXXFLAGS += -I$(NPLIB)/include
+
+LDFLAGS = `root-config --libs` -lMathMore
+LDFLAGS+= `$(NPLIB)/liblist`
+LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
+
+SRC= $(wildcard *.cc)
+INC= $(wildcard $(NPAINCLUDES)/*.hh)
+OBJ=$(SRC:.cc=.o)
+
+#all:$(EXEC)
+# @$(CPP) -o $@ -c $< $(CXXFLAGS)
+
+Analysis:$(OBJ) $(INC)
+ @$(CPP) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)
+ mv Analysis ../Analysis
+
+%.o: %.cc
+ @$(CPP) $(CXXFLAGS) -o $@ -c $<
+
+.PHONY: clean mrproper
+
+clean:
+ rm -rf *.o
+
+mrproper: clean
+ rm -rf $(EXEC)
diff --git a/NPLib/AnnularS1/TS1Data.cxx b/NPLib/AnnularS1/TS1Data.cxx
index e5224851600abf587d58a78e8cbce7909926d20d..6ac39c0919556ad3ff100df324b51ee6d6a8253f 100644
--- a/NPLib/AnnularS1/TS1Data.cxx
+++ b/NPLib/AnnularS1/TS1Data.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/AnnularS1/TS1Data.h b/NPLib/AnnularS1/TS1Data.h
index 225f6e714cbbb9ec411a9c56f71b76378bbf0eda..934581f4f74de532b9e16f6a8dcf81926fbf6d93 100644
--- a/NPLib/AnnularS1/TS1Data.h
+++ b/NPLib/AnnularS1/TS1Data.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/CalibrationManager/CalibrationManager.cxx b/NPLib/CalibrationManager/CalibrationManager.cxx
index 085c08285c60962b323ad2a89135c9683dd309bb..49fe4db9f20e8711784d7230d9e148132dfe9eea 100644
--- a/NPLib/CalibrationManager/CalibrationManager.cxx
+++ b/NPLib/CalibrationManager/CalibrationManager.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -25,6 +25,7 @@
#include <cstdlib>
#include <limits>
#include <cmath>
+#include <sstream>
//////////////////////////////////////////////////////////////////
CalibrationManager* CalibrationManager::instance = 0;
@@ -58,15 +59,17 @@ CalibrationManager::CalibrationManager(string configFileName)
return;
}
- else {
+ else
+ {
+ cout << "Reading list of file from :" << configFileName << endl;
while (!inputConfigFile.eof()) {
getline(inputConfigFile, lineBuffer);
// search for token giving the list of Root files to treat
- if ( lineBuffer.compare(0, 12, "CalibrationFilePath") == 0 ) {
+ if ( lineBuffer.compare(0, 19, "CalibrationFilePath") == 0 ) {
while (!inputConfigFile.eof()) {
inputConfigFile >> dataBuffer;
-
+
// ignore comment Line
if (dataBuffer.compare(0, 1, "%") == 0) {
inputConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
@@ -76,11 +79,12 @@ CalibrationManager::CalibrationManager(string configFileName)
AddFile(dataBuffer);
cout << "Adding file " << dataBuffer << " to Calibration" << endl;
}
- }
- }
- }
+ }
+ }
+ }
+ }
+ cout << "/////////////////////////////////" << endl;
}
- cout << "/////////////////////////////////" << endl;}
//////////////////////////////////////////////////////////////////
CalibrationManager::~CalibrationManager()
@@ -99,49 +103,56 @@ void CalibrationManager::LoadParameterFromFile()
{
ifstream CalibFile ;
string DataBuffer ;
+ string LineBuffer ;
for(unsigned int i = 0 ; i < fFileList.size() ; i++)
{
CalibFile.open( fFileList[i].c_str() );
- vector<double> Coeff ;
map<string,string>::iterator it ;
if(!CalibFile)
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
- cout << " WARNING: FILE " << fFileList[i] << " IS MISSING " << endl ;
+ cout << " WARNING: FILE " << fFileList[i] << " IS MISSING " << endl ;
cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
-
}
else while( !CalibFile.eof() )
+ {
+ // Read the file Line by line
+ getline(CalibFile, LineBuffer);
+
+ // Create a istringstream to manipulate the line easely
+ istringstream theLine (LineBuffer,istringstream::in);
+ theLine >> DataBuffer ;
+
+ // Comment support, comment symbole is %
+ if(DataBuffer.compare(0, 1, "%") == 0) {
+ CalibFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ // Search word in the token list
+ it=fToken.find(DataBuffer);
+
+ // if the word is find, values are read
+ if( it!=fToken.end() )
{
- CalibFile >> DataBuffer ;
-
- // Search word in the token list
- it=fToken.find(DataBuffer);
-
- // if the word is find, values are read
- if( it!=fToken.end() )
+ vector<double> Coeff ;
+ while( !theLine.eof() )
{
-
- Coeff.clear();
- while(DataBuffer!="\n")
- {
- CalibFile >> DataBuffer ; Coeff.push_back( atof(DataBuffer.c_str()) ) ;
- }
+ theLine >> DataBuffer ; Coeff.push_back( atof(DataBuffer.c_str()) ) ;
+ }
- // Check this parameter is not already define
- if( fCalibrationCoeff.find(it->second) != fCalibrationCoeff.end() ) cout << "WARNING: Parameter " << it->second << " Already found. It will be rewritted " << endl;
+ // Check this parameter is not already define
+ if( fCalibrationCoeff.find(it->second) != fCalibrationCoeff.end() )
+ cout << "WARNING: Parameter " << it->second << " Already found. It will be rewritted " << endl;
- // Add the list of Coeff to the Coeff map using Parameter Path as index
- fCalibrationCoeff[ it->second ] = Coeff ;
- }
-
+ // Add the list of Coeff to the Coeff map using Parameter Path as index
+ fCalibrationCoeff[ it->second ] = Coeff ;
}
- CalibFile.close() ;
+
+ }
+ CalibFile.close() ;
}
-
}
//////////////////////////////////////////////////////////////////
@@ -149,8 +160,12 @@ double CalibrationManager::ApplyCalibration(string ParameterPath , double RawVal
{
double CalibratedValue = 0 ;
map< string , vector<double> >::iterator it ;
+
+ // Find the good parameter in the Map
+ // Using Find method of stl is the fastest way
it = fCalibrationCoeff.find(ParameterPath) ;
+ // If the find methods return the end iterator it's mean the parameter was not found
if(it == fCalibrationCoeff.end() )
{
/* cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
@@ -160,8 +175,13 @@ double CalibrationManager::ApplyCalibration(string ParameterPath , double RawVal
return RawValue ;
}
+ // Else we take the second part of the element (first is index, ie: parameter path)
+ // Second is the vector of Coeff
vector<double> Coeff = it->second ;
+ // The vector size give the degree of calibration
+ // We just apply the coeff and returned the calibrated value
+
for(unsigned int i = 0 ; i < Coeff.size() ; i++)
{
CalibratedValue += Coeff[i]*pow(RawValue, (double)i);
diff --git a/NPLib/CalibrationManager/CalibrationManager.h b/NPLib/CalibrationManager/CalibrationManager.h
index 69df3a9f7f5e1ed8ee4519c1dbeaaa488cdd128d..9472a3123a1b70717593901ffe294a8c82cb853e 100644
--- a/NPLib/CalibrationManager/CalibrationManager.h
+++ b/NPLib/CalibrationManager/CalibrationManager.h
@@ -1,7 +1,7 @@
#ifndef ConfigManager_h
#define ConfigManager_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -33,7 +33,7 @@ using namespace std ;
class CalibrationManager
{
- public: // Constructor and Destructor
+ protected: // Constructor and Destructor are protected because the class is a singleton
CalibrationManager(string configFileName);
~CalibrationManager();
@@ -50,20 +50,20 @@ class CalibrationManager
inline void AddFile(string Path) { fFileList.push_back(Path) ;} ;
- public: // Declaration of Calibration
-
+ public: // Calibration Parameter Related
+
// call like : myCalibrationManager->AddParameter( "MUST2" ,"Telescope5_Si_X38_E", "T5_Si_X38_E" )
// return false if the token is not found in the file list
- bool AddParameter(string DetectorName , string ParameterName , string Token ) ;
+ bool AddParameter(string DetectorName , string ParameterName , string Token) ;
// call like : myCalibrationManager->ApplyCalibration( "MUST2/Telescope5_Si_X38_E" , RawEnergy )
// return the Calibrated value
double ApplyCalibration(string ParameterPath , double RawValue);
- public: // To be called after initialisation
- // Loop over the file list and catch the file used for calibration
- void LoadParameterFromFile();
+ public: // To be called after initialisation
+ // Loop over the file list and catch the file used for calibration
+ void LoadParameterFromFile();
private:
diff --git a/NPLib/DummyDetector/TDUMMYDetectorData.cxx b/NPLib/DummyDetector/TDUMMYDetectorData.cxx
index 1c15b58bdaa04b846a331cac3de2a18b26d11ee0..e1b58093152f087ab5427260393a7e89481f4045 100755
--- a/NPLib/DummyDetector/TDUMMYDetectorData.cxx
+++ b/NPLib/DummyDetector/TDUMMYDetectorData.cxx
@@ -1,7 +1,7 @@
#include <iostream>
#include "TDUMMYDetectorData.h"
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/DummyDetector/TDUMMYDetectorData.h b/NPLib/DummyDetector/TDUMMYDetectorData.h
index cc5c29fa18c8300f629b2c221e06e6d60d6e5d0e..e340a207342c71299d248d65c604b1a8b3d9124c 100755
--- a/NPLib/DummyDetector/TDUMMYDetectorData.h
+++ b/NPLib/DummyDetector/TDUMMYDetectorData.h
@@ -1,7 +1,7 @@
#ifndef __DUMMYDetectorDATA__
#define __DUMMYDetectorDATA__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/GaspardTracker.cxx b/NPLib/GASPARD/GaspardTracker.cxx
index 7c8aa36ebd1cf849dae9f0fbce6216dbe06987d6..b43e899f317ddd4535482fa5076b8b85929bdece 100644
--- a/NPLib/GASPARD/GaspardTracker.cxx
+++ b/NPLib/GASPARD/GaspardTracker.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/GaspardTracker.h b/NPLib/GASPARD/GaspardTracker.h
index 5de55bd49c53f3840d85969093c612f39452aa8c..346728d1b9858cac48786d5b56d03760680ad5a4 100644
--- a/NPLib/GASPARD/GaspardTracker.h
+++ b/NPLib/GASPARD/GaspardTracker.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/TGaspardTrackerData.cxx b/NPLib/GASPARD/TGaspardTrackerData.cxx
index 97ae3e74d203491dc7d409247d6830ea36f23385..dc88e6ab2f47cda90860540082cfc499d91250b8 100644
--- a/NPLib/GASPARD/TGaspardTrackerData.cxx
+++ b/NPLib/GASPARD/TGaspardTrackerData.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/TGaspardTrackerData.h b/NPLib/GASPARD/TGaspardTrackerData.h
index 3be1dc39375ce6418a552254a29d82d12df23bbc..bc45afa51b5a681629536f305c28da662e9a9ff3 100644
--- a/NPLib/GASPARD/TGaspardTrackerData.h
+++ b/NPLib/GASPARD/TGaspardTrackerData.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
index 8487f87751da5c57a52d2567d234347d9f7714bd..09c69515c360efc13d32aab0106566ece3f2976d 100644
--- a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
+++ b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/GASPARD/TGaspardTrackerPhysics.h b/NPLib/GASPARD/TGaspardTrackerPhysics.h
index 425d5e601e76a95d638d894893fa6492fc23cabd..0e9d115c39be1c0554ca118ef509fae0f20fac44 100644
--- a/NPLib/GASPARD/TGaspardTrackerPhysics.h
+++ b/NPLib/GASPARD/TGaspardTrackerPhysics.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/INSTALL b/NPLib/INSTALL
index b93077b70d97a8fa963e57928a130d76561fd0b8..ab085fec55c70335fa113f662cd9deb55755fad2 100644
--- a/NPLib/INSTALL
+++ b/NPLib/INSTALL
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/IORoot/RootInput.cxx b/NPLib/IORoot/RootInput.cxx
index 327343c5ff756da039a2ace4878c0dd9c14f4b1e..146c604482d17d2aa261040d91e09dc6e3ef0b1d 100644
--- a/NPLib/IORoot/RootInput.cxx
+++ b/NPLib/IORoot/RootInput.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/IORoot/RootInput.h b/NPLib/IORoot/RootInput.h
index 4a05b604704f8a417bb0d250573ad0b4887749b9..776b0d3b1879f3f7533583b2a5892e7a5ebafe20 100644
--- a/NPLib/IORoot/RootInput.h
+++ b/NPLib/IORoot/RootInput.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/IORoot/RootOutput.cxx b/NPLib/IORoot/RootOutput.cxx
index 625edb299312ab1829bc9feebe92a41f891850e2..f3b03096f2b656936164741e6b5bc084f520f789 100644
--- a/NPLib/IORoot/RootOutput.cxx
+++ b/NPLib/IORoot/RootOutput.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/IORoot/RootOutput.h b/NPLib/IORoot/RootOutput.h
index 0a56bf7cd93497b8c3c94acd8266debae152bf98..4d25455b13c31ce31af69276cfaa06ca4218f7da 100644
--- a/NPLib/IORoot/RootOutput.h
+++ b/NPLib/IORoot/RootOutput.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/InitialConditions/TInitialConditions.cxx b/NPLib/InitialConditions/TInitialConditions.cxx
index be597460ded1fc13e750111c15194a477cbf6d22..6dbd3411f69fa17165421413b996a5899f6dcfe0 100644
--- a/NPLib/InitialConditions/TInitialConditions.cxx
+++ b/NPLib/InitialConditions/TInitialConditions.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 468e941d2ea92152d55524ceccf0cda2e1089a47..67430e715990e4f61e8b006a8088f4c600e0f03a 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/InteractionCoordinates/TInteractionCoordinates.cxx b/NPLib/InteractionCoordinates/TInteractionCoordinates.cxx
index 624060e12b144eb54f64a5aeffbd5adcbba363eb..df7705ad5a909a40d2de9e7cd0706c5f283e1439 100644
--- a/NPLib/InteractionCoordinates/TInteractionCoordinates.cxx
+++ b/NPLib/InteractionCoordinates/TInteractionCoordinates.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/InteractionCoordinates/TInteractionCoordinates.h b/NPLib/InteractionCoordinates/TInteractionCoordinates.h
index ffd1fcc3c2eb4b0c10da88bdf23081d88802e64f..c3668039d9b5146942689c92086f95b0ff820aa6 100644
--- a/NPLib/InteractionCoordinates/TInteractionCoordinates.h
+++ b/NPLib/InteractionCoordinates/TInteractionCoordinates.h
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/MUST2/TMust2Data.cxx b/NPLib/MUST2/TMust2Data.cxx
index 308e52a0b8d581e4b55d819dd2cd584245e767f3..50fc48303522d333d30489dff7fe3db2cade995d 100644
--- a/NPLib/MUST2/TMust2Data.cxx
+++ b/NPLib/MUST2/TMust2Data.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/MUST2/TMust2Data.h b/NPLib/MUST2/TMust2Data.h
index 654e999cdb441c7e62ae992f3b80cef83acd05fd..f8540176b29b00d75b3a97567d3f8114a78f0bed 100644
--- a/NPLib/MUST2/TMust2Data.h
+++ b/NPLib/MUST2/TMust2Data.h
@@ -1,7 +1,7 @@
#ifndef __MUST2DATA__
#define __MUST2DATA__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -151,7 +151,7 @@ class TMust2Data : public TObject {
// CsI
//(E)
- UShort_t GetMMCsIEMult() {return fMM_CsIE_DetectorNbr.size();}
+ UShort_t GetMMCsIEMult() {return fMM_CsIE_DetectorNbr.size();}
UShort_t GetMMCsIEDetectorNbr(Int_t i) {return fMM_CsIE_DetectorNbr.at(i);}
UShort_t GetMMCsIECristalNbr(Int_t i) {return fMM_CsIE_CristalNbr.at(i);}
Double_t GetMMCsIEEnergy(Int_t i) {return fMM_CsIE_Energy.at(i);}
diff --git a/NPLib/MUST2/TMust2Physics.cxx b/NPLib/MUST2/TMust2Physics.cxx
index 4cd918223b2ae842225aa0c02b6217b5b48634cd..e64af68be2201350914283bb6ff10005ced77d27 100644
--- a/NPLib/MUST2/TMust2Physics.cxx
+++ b/NPLib/MUST2/TMust2Physics.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -16,8 +16,6 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
- * Only multiplicity one and multiplicity 2 are down. *
- * Improvment needed *
* *
*****************************************************************************/
#include "TMust2Physics.h"
@@ -40,17 +38,17 @@ using namespace LOCAL;
ClassImp(TMust2Physics)
///////////////////////////////////////////////////////////////////////////
TMust2Physics::TMust2Physics()
- {
- EventMultiplicity = 0 ;
- EventData = new TMust2Data ;
- EventPhysics = this ;
+ {
+ EventMultiplicity = 0 ;
+ EventData = new TMust2Data ;
+ EventPhysics = this ;
+ NumberOfTelescope = 0 ;
}
///////////////////////////////////////////////////////////////////////////
void TMust2Physics::BuildSimplePhysicalEvent()
{
BuildPhysicalEvent();
-
}
///////////////////////////////////////////////////////////////////////////
@@ -95,14 +93,24 @@ void TMust2Physics::BuildPhysicalEvent()
{
if(EventData->GetMMSiLiEDetectorNbr(j)==N)
{
- // if SiLi energy is above threshold check the compatibility
+ // SiLi energy is above threshold check the compatibility
if( fSiLi_E(EventData , j)>SiLi_E_Threshold )
{
+ // pad vs strip number match
if( Match_Si_SiLi( X, Y , EventData->GetMMSiLiEPadNbr(j) ) )
{
SiLi_N.push_back(EventData->GetMMSiLiEPadNbr(j)) ;
- SiLi_E.push_back(fSiLi_E(EventData , j)) ;
- SiLi_T.push_back(fSiLi_T(EventData , j)) ;
+ SiLi_E.push_back(fSiLi_E(EventData , j)) ;
+
+ // Look for associate energy
+ // Note: in case of use of SiLi "Orsay" time is not coded.
+ for(int k =0 ; k < EventData->GetMMSiLiTMult() ; k ++)
+ {
+ // Same Pad, same Detector
+ if( EventData->GetMMSiLiEPadNbr(j)==EventData->GetMMSiLiEPadNbr(k) && EventData->GetMMSiLiEDetectorNbr(j)==EventData->GetMMSiLiTDetectorNbr(k) )
+ {SiLi_T.push_back(fSiLi_T(EventData , k)) ; break ;}
+ }
+
check_SILI = true ;
}
@@ -114,14 +122,20 @@ void TMust2Physics::BuildPhysicalEvent()
{
if(EventData->GetMMCsIEDetectorNbr(j)==N)
{
- // ifCsI energy is above threshold check the compatibility
+ // CsI energy is above threshold check the compatibility
if( fCsI_T(EventData , j)>CsI_E_Threshold )
{
- if( Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
+ if(Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
{
CsI_N.push_back(EventData->GetMMCsIECristalNbr(j)) ;
- CsI_E.push_back(fCsI_E(EventData , j)) ;
- CsI_T.push_back(fCsI_T(EventData , j)) ;
+ CsI_E.push_back(fCsI_E(EventData , j)) ;
+
+ for(int k =0 ; k < EventData->GetMMCsITMult() ; k ++)
+ {
+ if( EventData->GetMMCsIECristalNbr(j)==EventData->GetMMCsITCristalNbr(k) && EventData->GetMMCsIEDetectorNbr(j)==EventData->GetMMCsITDetectorNbr(k) )
+ {CsI_T.push_back(fCsI_T(EventData , k)) ; break ;}
+ }
+
check_CSI = true ;
}
}
@@ -132,15 +146,15 @@ void TMust2Physics::BuildPhysicalEvent()
if(!check_SILI)
{
SiLi_N.push_back(0) ;
- SiLi_E.push_back(0) ;
- SiLi_T.push_back(0) ;
+ SiLi_E.push_back(-10000) ;
+ SiLi_T.push_back(-10000) ;
}
if(!check_CSI)
{
CsI_N.push_back(0) ;
- CsI_E.push_back(0) ;
- CsI_T.push_back(0) ;
+ CsI_E.push_back(-10000) ;
+ CsI_T.push_back(-10000) ;
}
}
@@ -177,6 +191,11 @@ vector < TVector2 > TMust2Physics :: Match_X_Y()
{
vector < TVector2 > ArrayOfGoodCouple ;
+ // Prevent code from treating very high multiplicity Event
+ // Those event are not physical anyway and that improve speed.
+ if( EventData->GetMMStripXEMult()>6 || EventData->GetMMStripYEMult()>6 )
+ return ArrayOfGoodCouple;
+
for(int i = 0 ; i < EventData->GetMMStripXEMult(); i++)
{
// if X value is above threshold, look at Y value
@@ -191,7 +210,7 @@ vector < TVector2 > TMust2Physics :: Match_X_Y()
// if same detector check energy
if ( EventData->GetMMStripXEDetectorNbr(i) == EventData->GetMMStripYEDetectorNbr(j) )
{
- // Look if energy match
+ // Look if energy match (within 10%)
if( ( fSi_X_E(EventData , i) - fSi_Y_E(EventData , j) ) / fSi_X_E(EventData , i) < 0.1 )
ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
}
@@ -744,10 +763,12 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 ,
NumberOfTelescope++;
// Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U = C_X1_Y1 - C_X128_Y1 ;
+ TVector3 U = C_X128_Y1 - C_X1_Y1 ;
+ double Ushift = (U.Mag()-98)/2. ;
U = U.Unit() ;
// Vector V on Telescope Face (parallele to X Strip)
- TVector3 V = C_X128_Y128 - C_X128_Y1 ;
+ TVector3 V = C_X1_Y128 - C_X1_Y1 ;
+ double Vshift = (V.Mag() -98)/2. ;
V = V.Unit() ;
// Position Vector of Strip Center
@@ -756,10 +777,9 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 ,
TVector3 Strip_1_1 ;
// Geometry Parameter
- double Face = 98 ; //mm
- double NumberOfStrip = 128 ;
+ double Face = 98 ; //mm
+ double NumberOfStrip = 128 ;
double StripPitch = Face/NumberOfStrip ; //mm
-
// Buffer object to fill Position Array
vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
@@ -768,8 +788,9 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 ,
vector< vector< double > > OneTelescopeStripPositionZ ;
// Moving StripCenter to 1.1 corner:
- Strip_1_1 = C_X128_Y1 + (U+V) * (StripPitch/2.) ;
-
+ Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.) ;
+ Strip_1_1+= U*Ushift+V*Vshift ;
+
for( int i = 0 ; i < 128 ; i++ )
{
lineX.clear() ;
@@ -780,7 +801,6 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 ,
{
StripCenter = Strip_1_1 + StripPitch*( i*U + j*V ) ;
//StripCenter += -TargetPosition ;
-
lineX.push_back( StripCenter.X() ) ;
lineY.push_back( StripCenter.Y() ) ;
lineZ.push_back( StripCenter.Z() ) ;
@@ -791,8 +811,6 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 ,
OneTelescopeStripPositionZ.push_back(lineZ) ;
}
-
-
StripPositionX.push_back( OneTelescopeStripPositionX ) ;
StripPositionY.push_back( OneTelescopeStripPositionY ) ;
StripPositionZ.push_back( OneTelescopeStripPositionZ ) ;
@@ -835,7 +853,7 @@ void TMust2Physics::AddTelescope( double theta ,
W = C.Unit() ;
U = W .Cross ( P ) ;
- V = W .Cross ( U );
+ V = W .Cross ( U );
U = U.Unit();
V = V.Unit();
diff --git a/NPLib/MUST2/TMust2Physics.h b/NPLib/MUST2/TMust2Physics.h
index 1f649c60c2ae4298b4a44380e7fdb5bf7a9d903e..67cb9c31633b521d3ce59a5f0852baf9a329c7ec 100644
--- a/NPLib/MUST2/TMust2Physics.h
+++ b/NPLib/MUST2/TMust2Physics.h
@@ -1,7 +1,7 @@
#ifndef TMUST2PHYSICS_H
#define TMUST2PHYSICS_H
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/Plastic/TPlasticPhysics.cxx b/NPLib/Plastic/TPlasticPhysics.cxx
index 093f1ee58f6e74d8d2d0b1353c2b662821bfcbe0..394c52311177561dfeb8c4a3c9909f4dcbe0c119 100644
--- a/NPLib/Plastic/TPlasticPhysics.cxx
+++ b/NPLib/Plastic/TPlasticPhysics.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -76,18 +76,20 @@ void TPlasticPhysics::ReadConfiguration(string Path)
string LineBuffer ;
string DataBuffer ;
- double Theta = 0 , Phi = 0 , R = 0 , Thickness = 0 , Radius = 0 , LeadThickness = 0;
- string Scintillator ;
-
- bool check_Theta = false ;
- bool check_Phi = false ;
- bool check_R = false ;
+ bool check_Theta = false ;
+ bool check_Phi = false ;
+ bool check_R = false ;
bool check_Thickness = false ;
- bool check_Radius = false ;
- bool check_LeadThickness = false ;
+ bool check_Radius = false ;
+ bool check_LeadThickness = false ;
bool check_Scintillator = false ;
- bool ReadingStatus = false ;
-
+ bool check_Height = false ;
+ bool check_Width = false ;
+ bool check_Shape = false ;
+ bool check_X = false ;
+ bool check_Y = false ;
+ bool check_Z = false ;
+ bool ReadingStatus = false ;
while (!ConfigFile.eof())
{
@@ -117,54 +119,96 @@ void TPlasticPhysics::ReadConfiguration(string Path)
// Finding another telescope (safety), toggle out
else if (DataBuffer.compare(0, 6, "Plastic") == 0) {
- cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
+ cout << "WARNING: Another Detector is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
ReadingStatus = false ;
}
- //Angle method
+ //Angle method
else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
check_Theta = true;
ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
+ cout << "Theta: " << atof(DataBuffer.c_str()) << "deg" << endl;
}
else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
check_Phi = true;
ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
+ cout << "Phi: " << atof( DataBuffer.c_str() ) << "deg" << endl;
}
else if (DataBuffer.compare(0, 2, "R=") == 0) {
check_R = true;
ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
+ cout << "R: " << atof( DataBuffer.c_str() ) << "mm" << endl;
+ }
+
+ //Position method
+ else if (DataBuffer.compare(0, 2, "X=") == 0) {
+ check_X = true;
+ ConfigFile >> DataBuffer ;
+ cout << "X: " << atof( DataBuffer.c_str() ) << "mm" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 2, "Y=") == 0) {
+ check_Y = true;
+ ConfigFile >> DataBuffer ;
+ cout << "Y: " << atof( DataBuffer.c_str() ) << "mm"<< endl;
+ }
+
+ else if (DataBuffer.compare(0, 2, "Z=") == 0) {
+ check_Z = true;
+ ConfigFile >> DataBuffer ;
+ cout << "Z: " << atof( DataBuffer.c_str() ) << "mm" << endl;
}
+
+ //General
+ else if (DataBuffer.compare(0, 6, "Shape=") == 0) {
+ check_Shape = true;
+ ConfigFile >> DataBuffer ;
+ cout << "Shape: " << DataBuffer << endl;
+ }
+
+ // Cylindrical shape
else if (DataBuffer.compare(0, 7, "Radius=") == 0) {
check_Radius = true;
ConfigFile >> DataBuffer ;
- Radius = atof(DataBuffer.c_str()) ;
+ cout << "Plastic Radius: " << atof( DataBuffer.c_str() ) << "mm" << endl;
+ }
+
+ // Squared shape
+ else if (DataBuffer.compare(0, 7, "Width=") == 0) {
+ check_Width = true;
+ ConfigFile >> DataBuffer ;
+ cout << "Plastic Width: " <<atof( DataBuffer.c_str() ) << "mm" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "Height=") == 0) {
+ check_Height = true;
+ ConfigFile >> DataBuffer ;
+ cout << "Plastic Height: " << atof( DataBuffer.c_str() ) << "mm" << endl;
}
+ // Common
else if (DataBuffer.compare(0, 10, "Thickness=") == 0) {
check_Thickness = true;
ConfigFile >> DataBuffer ;
- Thickness = atof(DataBuffer.c_str()) ;
+ cout << "Plastic Thickness: " << atof( DataBuffer.c_str() ) << "mm" << endl;
}
else if (DataBuffer.compare(0, 13, "Scintillator=") == 0) {
check_Scintillator = true ;
ConfigFile >> DataBuffer ;
- Scintillator = DataBuffer ;
+ cout << "Plastic Scintillator type: " << DataBuffer << endl;
}
else if (DataBuffer.compare(0, 14, "LeadThickness=") == 0) {
check_LeadThickness = true;
ConfigFile >> DataBuffer ;
- LeadThickness = atof(DataBuffer.c_str()) ;
+ cout << "Lead Thickness : " << atof( DataBuffer.c_str() ) << "mm" << endl;
}
-
- ///////////////////////////////////////////////////
+
+ ///////////////////////////////////////////////////
// If no Detector Token and no comment, toggle out
else
{ReadingStatus = false; cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl ;}
@@ -172,21 +216,26 @@ void TPlasticPhysics::ReadConfiguration(string Path)
/////////////////////////////////////////////////
// If All necessary information there, toggle out
- if ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator)
+ if ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator && check_Height && check_Width && check_Shape && check_X && check_Y && check_Z )
{
NumberOfDetector++;
- // Reinitialisation of Check Boolean
-
- check_Theta = false ;
- check_Phi = false ;
- check_R = false ;
- check_Thickness = false ;
- check_Radius = false ;
- check_LeadThickness = false ;
- check_Scintillator = false ;
- ReadingStatus = false ;
- cout << "///"<< endl ;
+ // Reinitialisation of Check Boolean
+ check_Theta = false ;
+ check_Phi = false ;
+ check_R = false ;
+ check_Thickness = false ;
+ check_Radius = false ;
+ check_LeadThickness = false ;
+ check_Scintillator = false ;
+ check_Height = false ;
+ check_Width = false ;
+ check_Shape = false ;
+ check_X = false ;
+ check_Y = false ;
+ check_Z = false ;
+ ReadingStatus = false ;
+ cout << "///"<< endl ;
}
}
diff --git a/NPLib/Plastic/TPlasticPhysics.h b/NPLib/Plastic/TPlasticPhysics.h
index a00c34437ad3ac4caa7aee9195a241ec42fd6fdf..9bf4f9af0a4a093ca78c4284e44a84b352f1a2ca 100644
--- a/NPLib/Plastic/TPlasticPhysics.h
+++ b/NPLib/Plastic/TPlasticPhysics.h
@@ -1,7 +1,7 @@
#ifndef __PlasticPhysics__
#define __PlasticPhysics__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/README b/NPLib/README
index 5723d77390845d9aa5512a46bb0bdda223619ffb..03dd1b66720c58bef45d130bbf592f3f03015d67 100644
--- a/NPLib/README
+++ b/NPLib/README
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/SSSD/TSSSDData.cxx b/NPLib/SSSD/TSSSDData.cxx
index 66d7b22393986cb749027650c8e7f102f822905d..f672bae35865888cdc0de6885006fd4ae62e91dd 100644
--- a/NPLib/SSSD/TSSSDData.cxx
+++ b/NPLib/SSSD/TSSSDData.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/SSSD/TSSSDData.h b/NPLib/SSSD/TSSSDData.h
index 7e90e7f0515b029ea4799ac39516cd4d5c6980c5..85de5af28b45547d82b15aa7b8d1cee23046bc65 100644
--- a/NPLib/SSSD/TSSSDData.h
+++ b/NPLib/SSSD/TSSSDData.h
@@ -1,7 +1,7 @@
#ifndef __SSSDDATA__
#define __SSSDDATA__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/SSSD/TSSSDPhysics.cxx b/NPLib/SSSD/TSSSDPhysics.cxx
index 63be1855e636503767ddc451f875d8821984c922..9ba970a25192fdd0ecb925002a526f0fee21e616 100644
--- a/NPLib/SSSD/TSSSDPhysics.cxx
+++ b/NPLib/SSSD/TSSSDPhysics.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/SSSD/TSSSDPhysics.h b/NPLib/SSSD/TSSSDPhysics.h
index 21bc790fb7067a44b6beef2e3f1da94890a10b97..fa341dc97e52f53da1a3e3d424a379b76a656cc0 100644
--- a/NPLib/SSSD/TSSSDPhysics.h
+++ b/NPLib/SSSD/TSSSDPhysics.h
@@ -1,7 +1,7 @@
#ifndef __SSSDPhysics__
#define __SSSDPhysics__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/Tools/NPNucleus.cxx b/NPLib/Tools/NPNucleus.cxx
index 6e59692165907d6b360039e897b1cea54318fee9..3deac59b33a5cf9a866b451374397f7532d44757 100644
--- a/NPLib/Tools/NPNucleus.cxx
+++ b/NPLib/Tools/NPNucleus.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/Tools/NPNucleus.h b/NPLib/Tools/NPNucleus.h
index a56d946777fbb02dfb09f251b2831b71268b45da..79142544abe6e77ad0d78e706f35e762d72541fe 100644
--- a/NPLib/Tools/NPNucleus.h
+++ b/NPLib/Tools/NPNucleus.h
@@ -1,7 +1,7 @@
#ifndef __NUCLEUS__
#define __NUCLEUS__
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPLib/Tools/NPReaction.cxx b/NPLib/Tools/NPReaction.cxx
index 13121b3ec5736bede26ffdb8bb3803707e436e9b..140776985584bc51f8edf0efd4e200a838f9108f 100644
--- a/NPLib/Tools/NPReaction.cxx
+++ b/NPLib/Tools/NPReaction.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -54,15 +54,15 @@ Reaction::Reaction()
{
//------------- Default Constructor -------------
- fNoy1 = new Nucleus();
- fNoy2 = new Nucleus();
- fNoy3 = new Nucleus();
- fNoy4 = new Nucleus();
- fBeamEnergy = 0;
- fThetaCM = 0;
- fExcitation = 0;
- fQValue = 0;
- initializePrecomputeVariable();
+ fNuclei1 = new Nucleus() ;
+ fNuclei2 = new Nucleus() ;
+ fNuclei3 = new Nucleus() ;
+ fNuclei4 = new Nucleus() ;
+ fBeamEnergy = 0 ;
+ fThetaCM = 0 ;
+ fExcitation = 0 ;
+ fQValue = 0 ;
+ initializePrecomputeVariable() ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -76,18 +76,18 @@ void Reaction::SetEveryThing(string name1, string name2, string name3, string na
{
//------------- Constructor with nuclei names and beam energy ------------
- fNoy1 = new Nucleus(name1);
- fNoy2 = new Nucleus(name2);
- fNoy3 = new Nucleus(name3);
- fNoy4 = new Nucleus(name4);
+ fNuclei1 = new Nucleus(name1);
+ fNuclei2 = new Nucleus(name2);
+ fNuclei3 = new Nucleus(name3);
+ fNuclei4 = new Nucleus(name4);
fBeamEnergy = BeamEnergy;
fThetaCM = 0;
fExcitation = ExcitationEnergy;
- fQValue = ( fNoy1->GetMassExcess() + fNoy2->GetMassExcess()
- - fNoy3->GetMassExcess() - fNoy4->GetMassExcess()) / 1000;
+ fQValue = ( fNuclei1->GetMassExcess() + fNuclei2->GetMassExcess()
+ - fNuclei3->GetMassExcess() - fNuclei4->GetMassExcess()) / 1000;
- int masse = fNoy1->GetA() + fNoy2->GetA() - fNoy3->GetA() - fNoy4->GetA();
- int charge = fNoy1->GetZ() + fNoy2->GetZ() - fNoy3->GetZ() - fNoy4->GetZ();
+ int masse = fNuclei1->GetA() + fNuclei2->GetA() - fNuclei3->GetA() - fNuclei4->GetA();
+ int charge = fNuclei1->GetZ() + fNuclei2->GetZ() - fNuclei3->GetZ() - fNuclei4->GetZ();
if (masse || charge) cout << "Problem with charge or mass conservation" << endl;
///Read the differential cross section
@@ -122,14 +122,31 @@ Reaction::~Reaction()
{
//------------- Default Destructor ------------
- delete fNoy1;
- delete fNoy2;
- delete fNoy3;
- delete fNoy4;
+ delete fNuclei1;
+ delete fNuclei2;
+ delete fNuclei3;
+ delete fNuclei4;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+bool Reaction::CheckKinematic()
+ {
+ // Check if kinematics is allowed
+
+ // case of inverse kinematics
+ double theta = fThetaCM;
+ if (m1 > m2) theta = M_PI - fThetaCM;
+
+ // total and kinetic energies in the lab
+ double W3lab = W3cm * G * (1 + B*beta3cm*cos(theta));
+ double W4lab = W4cm * G * (1 + B*beta4cm*cos(theta + M_PI));
+ // test for total energy conversion
+ if (fabs(WtotLab - (W3lab+W4lab)) > 1e-6)
+ return false ;
+
+ else return true ;
+ }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Reaction::KineRelativistic(double &ThetaLab3, double &EnergieLab3,
@@ -138,7 +155,7 @@ void Reaction::KineRelativistic(double &ThetaLab3, double &EnergieLab3,
// 2-body relativistic kinematics: direct + inverse
// EnergieLab3,4 : lab energy in MeV of the 2 ejectiles
// ThetaLab3,4 : angles in rad
-
+
// case of inverse kinematics
double theta = fThetaCM;
if (m1 > m2) theta = M_PI - fThetaCM;
@@ -168,36 +185,15 @@ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab) cons
{
// EnergyLab in MeV
// ThetaLab in rad
-
double P1 = sqrt(2*m1*fBeamEnergy+(fBeamEnergy*fBeamEnergy)) ;
double P3 = sqrt(2*m3*EnergyLab+(EnergyLab*EnergyLab)) ;
double P4 = sqrt(P1*P1+P3*P3-(2*P1*P3*cos(ThetaLab))) ;
double E4 = fBeamEnergy+m1+m2-(EnergyLab+m3) ;
double m4e = sqrt((E4*E4)-(P4*P4)) ;
- double Eex= m4e-m4 ;
+ double Eex= m4e-fNuclei4->Mass() ;
return Eex;
}
-
-/*double Reaction::ReconstructThetaHeavy(double EnergyLab, double ThetaLab)
- {
- // EnergyLab in MeV
- // ThetaLab in rad
-
- double m1 = fNoy1->Mass() ;
- double m2 = fNoy2->Mass() ;
- double m3 = fNoy3->Mass() ;
- double m4 = fNoy4->Mass() ;
-
- double P1 = sqrt(2*m1*fBeamEnergy+(fBeamEnergy*fBeamEnergy)) ;
- double P3 = sqrt(2*m3*EnergyLab+(EnergyLab*EnergyLab)) ;
- double P4 = sqrt(P1*P1+P3*P3-(2*P1*P3*cos(ThetaLab))) ;
-
-
- double cosThetaHeavy = (P1 - P4cos(hetaLab))/P3 ;
-
- }*/
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//Return ThetaCM
@@ -234,8 +230,8 @@ void Reaction::Print() const
{
// Print informations concerning the reaction
- cout << "Reaction : " << fNoy2->GetName() << "(" << fNoy1->GetName()
- << "," << fNoy3->GetName() << ")" << fNoy4->GetName() << " @ "
+ cout << "Reaction : " << fNuclei2->GetName() << "(" << fNuclei1->GetName()
+ << "," << fNuclei3->GetName() << ")" << fNuclei4->GetName() << " @ "
<< fBeamEnergy << " MeV"
<< endl ;
@@ -365,11 +361,11 @@ void Reaction::ReadConfigurationFile(string Path)
void Reaction::initializePrecomputeVariable()
- {
- m1 = fNoy1->Mass();
- m2 = fNoy2->Mass();
- m3 = fNoy3->Mass();
- m4 = fNoy4->Mass() + fExcitation;
+ {
+ m1 = fNuclei1->Mass();
+ m2 = fNuclei2->Mass();
+ m3 = fNuclei3->Mass();
+ m4 = fNuclei4->Mass() + fExcitation;
// center-of-mass velocity
WtotLab = (fBeamEnergy + m1) + m2;
@@ -390,7 +386,6 @@ void Reaction::initializePrecomputeVariable()
// Constants of the kinematics
K3 = B / beta3cm;
K4 = B / beta4cm;
-
}
diff --git a/NPLib/Tools/NPReaction.h b/NPLib/Tools/NPReaction.h
index bb8f8932e62016176f8bd959bb94d29374da07af..771aedbea5e67073f181fcc4fed4424199a0fede 100644
--- a/NPLib/Tools/NPReaction.h
+++ b/NPLib/Tools/NPReaction.h
@@ -54,10 +54,10 @@ namespace NPL
void ReadConfigurationFile(string Path) ;
private:
- Nucleus *fNoy1 ; // Beam
- Nucleus *fNoy2 ; // Target
- Nucleus *fNoy3 ; // Light ejectile
- Nucleus *fNoy4 ; // Heavy ejectile
+ Nucleus *fNuclei1 ; // Beam
+ Nucleus *fNuclei2 ; // Target
+ Nucleus *fNuclei3 ; // Light ejectile
+ Nucleus *fNuclei4 ; // Heavy ejectile
double fQValue ; // Q-value in MeV
double fBeamEnergy ; // Beam energy in MeV
double fThetaCM ; // Center-of-mass angle in radian
@@ -67,19 +67,19 @@ namespace NPL
public:
// Getters and Setters
- void SetBeamEnergy (double efais) {fBeamEnergy = efais;initializePrecomputeVariable();}
- void SetThetaCM (double angle) {fThetaCM = angle;}
- void SetExcitation (double exci) {fExcitation = exci;initializePrecomputeVariable();}
- double GetBeamEnergy() const {return fBeamEnergy;}
- double GetThetaCM() const {return fThetaCM;}
- double GetExcitation() const {return fExcitation;}
- double GetQValue() const {return fQValue;}
- Nucleus* GetNucleus1() const {return fNoy1;}
- Nucleus* GetNucleus2() const {return fNoy2;}
- Nucleus* GetNucleus3() const {return fNoy3;}
- Nucleus* GetNucleus4() const {return fNoy4;}
- double* GetCrossSection() const {return CrossSection;}
- int GetCrossSectionSize() const {return CrossSectionSize;}
+ void SetBeamEnergy (double eBeam) {fBeamEnergy = eBeam;initializePrecomputeVariable();}
+ void SetThetaCM (double angle) {fThetaCM = angle;initializePrecomputeVariable();}
+ void SetExcitation (double exci) {fExcitation = exci;initializePrecomputeVariable();}
+ double GetBeamEnergy() const {return fBeamEnergy;}
+ double GetThetaCM() const {return fThetaCM;}
+ double GetExcitation() const {return fExcitation;}
+ double GetQValue() const {return fQValue;}
+ Nucleus* GetNucleus1() const {return fNuclei1;}
+ Nucleus* GetNucleus2() const {return fNuclei2;}
+ Nucleus* GetNucleus3() const {return fNuclei3;}
+ Nucleus* GetNucleus4() const {return fNuclei4;}
+ double* GetCrossSection() const {return CrossSection;}
+ int GetCrossSectionSize() const {return CrossSectionSize;}
@@ -111,6 +111,10 @@ namespace NPL
public: // Kinematics
+ // Check that the reaction is alowed
+ bool CheckKinematic();
+
+
// Compute ThetaLab and EnergyLab for product of reaction
void KineRelativistic( double &ThetaLab3 ,
double &EnergieLab3 ,
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 398198c956ba99c9791339f38dd3d3291b79f9b0..92991cdee1f40dc73a2d35a2e8955d8452f7aac5 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -302,6 +302,7 @@ void DetectorManager::InitializeRootOutput()
void DetectorManager::AddDetector(string DetectorName , VDetector* newDetector)
{
m_Detector[DetectorName] = newDetector;
+ newDetector->AddParameterToCalibrationManager();
}
diff --git a/NPSimulation/include/AnnularS1.hh b/NPSimulation/include/AnnularS1.hh
index 2f5ffb7dd9fa31ac91cd13e32770daa10e164328..71e6c97b25b052f8d40163e76d27177f99d239f4 100644
--- a/NPSimulation/include/AnnularS1.hh
+++ b/NPSimulation/include/AnnularS1.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/AnnularS1Scorers.hh b/NPSimulation/include/AnnularS1Scorers.hh
index c38cdb2fd0ec21eb90bdc58162b3a78b8748bedf..d689a2c86dd83fca1a455275cbd47f785039c55e 100644
--- a/NPSimulation/include/AnnularS1Scorers.hh
+++ b/NPSimulation/include/AnnularS1Scorers.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/DetectorConstruction.hh b/NPSimulation/include/DetectorConstruction.hh
index 8a65ee3abfe055a6a3f392dbd1ef0de6398a9219..6f305e6f69708253eaae6d71d203422a54c582cd 100644
--- a/NPSimulation/include/DetectorConstruction.hh
+++ b/NPSimulation/include/DetectorConstruction.hh
@@ -1,7 +1,7 @@
#ifndef DetectorConstruction_H
#define DetectorConstruction_H 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/DummyDetector.hh b/NPSimulation/include/DummyDetector.hh
index 700cec3a1a6302d36d641aa99178eb9d128f8f49..21b01f77daed6aea31d55bf7b239ff780b4ceed7 100755
--- a/NPSimulation/include/DummyDetector.hh
+++ b/NPSimulation/include/DummyDetector.hh
@@ -1,7 +1,7 @@
#ifndef DUMMYDetector_h
#define DUMMYDetector_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/EventAction.hh b/NPSimulation/include/EventAction.hh
index cbcf323de4775ea4f09a745dfc09ec96b4140cbb..6cc16bd694b129b740b733a9d4a792945f27b21b 100644
--- a/NPSimulation/include/EventAction.hh
+++ b/NPSimulation/include/EventAction.hh
@@ -1,7 +1,7 @@
#ifndef EventAction_h
#define EventAction_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/EventGeneratorBeam.hh b/NPSimulation/include/EventGeneratorBeam.hh
index ed5de49b6ff184d4ec702a3442c7bda4e0d9fe09..e18dba96b3f3f5a5573acc982280f4011999cb8f 100644
--- a/NPSimulation/include/EventGeneratorBeam.hh
+++ b/NPSimulation/include/EventGeneratorBeam.hh
@@ -1,7 +1,7 @@
#ifndef EventGeneratorBeam_h
#define EventGeneratorBeam_h
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -57,6 +57,8 @@ private: // TTree to store initial value of beam and reaction
private: // Source parameter
G4ParticleDefinition* m_particle; // Kind of particle to shoot
+ G4int m_beamA;
+ G4int m_beamZ;
G4double m_BeamEnergy;
G4double m_BeamEnergySpread;
G4double m_SigmaX;
diff --git a/NPSimulation/include/EventGeneratorIsotropic.hh b/NPSimulation/include/EventGeneratorIsotropic.hh
index 03210f62d45df4ccc48126c8824eeff350d85ea0..cc6b644b8b1fd9d841b606da29243a38c332035d 100644
--- a/NPSimulation/include/EventGeneratorIsotropic.hh
+++ b/NPSimulation/include/EventGeneratorIsotropic.hh
@@ -1,7 +1,7 @@
#ifndef EventGeneratorIsotropic_h
#define EventGeneratorIsotropic_h
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/EventGeneratorPhaseSpace.hh b/NPSimulation/include/EventGeneratorPhaseSpace.hh
new file mode 100755
index 0000000000000000000000000000000000000000..600d6cae02c9b4293c595b42d054807c4fdad850
--- /dev/null
+++ b/NPSimulation/include/EventGeneratorPhaseSpace.hh
@@ -0,0 +1,106 @@
+#ifndef EventGeneratorPhaseSpace_H
+#define EventGeneratorPhaseSpace_H
+/*****************************************************************************
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : Feb 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated pure phase space event and *
+ * especially evaluate the background contribution of different phase space *
+ * channel during an experiment. Manage up to 18 body phase space *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+
+// G4 headers
+#include "G4ParticleDefinition.hh"
+
+// NPSimulation
+#include "VEventGenerator.hh"
+#include "Target.hh"
+
+// NPLib
+#include "TInitialConditions.h"
+
+using namespace std;
+
+class EventGeneratorPhaseSpace : public VEventGenerator
+{
+ public: // Constructors and Destructors
+
+ // Default constructor used to allocate memory
+ EventGeneratorPhaseSpace();
+
+ // This is the constructor to be used
+ EventGeneratorPhaseSpace( G4ParticleDefinition* BeamParticle ,
+ G4ParticleDefinition* TargetParticle ,
+ vector<G4ParticleDefinition*> ReactionProducts ,
+ double BeamEnergy ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY );
+
+ // Default Destructor
+ virtual ~EventGeneratorPhaseSpace();
+
+ public: // Inherit from VEventGenerator class
+ void ReadConfiguration(string) ;
+ void GenerateEvent(G4Event*, G4ParticleGun*) ;
+ void SetTarget(Target* Target);
+
+ private: // Particle Shoot Option
+ bool m_ShootLight ;
+ bool m_ShootHeavy ;
+ bool m_ShootDecayProduct ;
+
+ private: // Value to store in the output tree : initial value of beam and reaction and phase space event weight
+ TInitialConditions* m_InitConditions;
+ // When the Phase Space Generator is called, the weight of the current configuration is return and stored in this variable
+ // Spectrum then need to be weighted by this paramater to be realistic
+ // NB: This procedure avoid long calculation time of the rejection methods previously instantiate and therefore allow simulation of manybody phase space decay
+ double m_EventWeight ;
+
+ private: // Target Parameter
+ Target* m_Target;
+
+ private: // Experimental Parameter
+ G4ParticleDefinition* m_BeamParticle ;
+ G4ParticleDefinition* m_TargetParticle ;
+ double m_BeamEnergy ;
+ double m_BeamEnergySpread ;
+ double m_SigmaX ;
+ double m_SigmaY ;
+ double m_SigmaThetaX ;
+ double m_SigmaPhiY ;
+ vector<G4ParticleDefinition*> m_ReactionProducts ;
+
+ private: //Other
+ void InitializeRootOutput() ;
+
+ void SetEverything( G4ParticleDefinition* BeamParticle ,
+ G4ParticleDefinition* TargetParticle ,
+ vector<G4ParticleDefinition*> ReactionProducts ,
+ double BeamEnergy ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY );
+};
+
+#endif
diff --git a/NPSimulation/include/EventGeneratorTransfert.hh b/NPSimulation/include/EventGeneratorTransfert.hh
index 8c570880d17397018574fe2784c2d97d24adc9a5..1239a93946a42db217d393342365100505065c77 100644
--- a/NPSimulation/include/EventGeneratorTransfert.hh
+++ b/NPSimulation/include/EventGeneratorTransfert.hh
@@ -1,7 +1,7 @@
#ifndef EventGeneratorTransfert_H
#define EventGeneratorTransfert_H
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/EventGeneratorTransfertToResonance.hh b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
index 2a5f1c30f8c822ffae01140825a6fbd6c346b0fe..dbc1873691d164432e23386a1e73f77c03dc4b62 100644
--- a/NPSimulation/include/EventGeneratorTransfertToResonance.hh
+++ b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
@@ -1,7 +1,7 @@
#ifndef EventGeneratorTransfertToResonance_H
#define EventGeneratorTransfertToResonance_H
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -99,13 +99,17 @@ class EventGeneratorTransfertToResonance : public VEventGenerator
Reaction* m_Reaction ;
private: // Resonance decay channel
- double m_ResonanceWidth ;
- double m_ResonanceMean ;
- int m_ResonanceDecayZ ;
- int m_ResonanceDecayA ;
+ double m_ResonanceWidth ;
+ double m_ResonanceMean ;
+ int m_ResonanceDecayZ ;
+ int m_ResonanceDecayA ;
+
+ // When the Phase Space Generator is called, the weight of the current configuration is return and stored in this variable
+ // Spectrum then need to be weighted by this paramater to be realistic
+ // NB: This procedure avoid long calculation time of the rejection methods previously instantiate and therefore allow simulation of manybody phase space decay
+ double m_EventWeight ;
//Other
- void Print() const ;
void InitializeRootOutput() ;
void ResonanceDecay( G4double EnergyHeavy ,
G4double ThetaHeavy ,
diff --git a/NPSimulation/include/GaspardScorers.hh b/NPSimulation/include/GaspardScorers.hh
index e009528060643ba63a34595355439f353ea93703..00ea23df1dc2153c239dd368a86a3293837cea99 100644
--- a/NPSimulation/include/GaspardScorers.hh
+++ b/NPSimulation/include/GaspardScorers.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTracker.hh b/NPSimulation/include/GaspardTracker.hh
index 0f094da67b2d45360570f61e480a214a809dfb20..4ca87d728a57b715068469d799170ceb8ebe89d8 100644
--- a/NPSimulation/include/GaspardTracker.hh
+++ b/NPSimulation/include/GaspardTracker.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTrackerAnnular.hh b/NPSimulation/include/GaspardTrackerAnnular.hh
index 0437e4bad965c6ca04d1738e511046f5594fa643..a07513fd035fc4cb2c4dd51452778bb6f0e4289d 100644
--- a/NPSimulation/include/GaspardTrackerAnnular.hh
+++ b/NPSimulation/include/GaspardTrackerAnnular.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTrackerDummyShape.hh b/NPSimulation/include/GaspardTrackerDummyShape.hh
index aafc1c50503b4e01502ccfb5b20800400d428fc8..cb7c3501587c4b5ba885871a0e5421ad869c0cac 100644
--- a/NPSimulation/include/GaspardTrackerDummyShape.hh
+++ b/NPSimulation/include/GaspardTrackerDummyShape.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTrackerModule.hh b/NPSimulation/include/GaspardTrackerModule.hh
index 0281068851aaeece058d182987dbb7edebf66a97..489643a65c76c8eba02955a0c85c1e46d5245fbb 100644
--- a/NPSimulation/include/GaspardTrackerModule.hh
+++ b/NPSimulation/include/GaspardTrackerModule.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTrackerSquare.hh b/NPSimulation/include/GaspardTrackerSquare.hh
index 180b0d1b125923308a2a618e0db785290f8ca763..c8c33b929f84b82fc31d3e7c164d671aefdf3d33 100644
--- a/NPSimulation/include/GaspardTrackerSquare.hh
+++ b/NPSimulation/include/GaspardTrackerSquare.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GaspardTrackerTrapezoid.hh b/NPSimulation/include/GaspardTrackerTrapezoid.hh
index 05bf5bef1a4a9c3d220ad3be5a6b5d6cc497e2fc..ac0c728ca17de05491e96f08e4921253f7afe114 100644
--- a/NPSimulation/include/GaspardTrackerTrapezoid.hh
+++ b/NPSimulation/include/GaspardTrackerTrapezoid.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/GeneralScorers.hh b/NPSimulation/include/GeneralScorers.hh
index 70dcc65e51c9bfd69172aa7eb0e40f22a49492f0..063db646c24ce97ffbf86aa089dca1d42eef9d31 100644
--- a/NPSimulation/include/GeneralScorers.hh
+++ b/NPSimulation/include/GeneralScorers.hh
@@ -1,7 +1,7 @@
#ifndef GeneralScorers_h
#define GeneralScorers_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/MUST2Array.hh b/NPSimulation/include/MUST2Array.hh
index 71196c2f5cb5c84244ea51df37bf401b7d8fe128..b2d794b383625789c9748755ee72c6ff045a1e9a 100644
--- a/NPSimulation/include/MUST2Array.hh
+++ b/NPSimulation/include/MUST2Array.hh
@@ -1,7 +1,7 @@
#ifndef MUST2Array_h
#define MUST2Array_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/Must2Scorers.hh b/NPSimulation/include/Must2Scorers.hh
index f5867f6ad489637d7873d74f84b84cfcd55c66df..2eb5431d8dfceb91e36736a76959945e9ae986d9 100644
--- a/NPSimulation/include/Must2Scorers.hh
+++ b/NPSimulation/include/Must2Scorers.hh
@@ -1,7 +1,7 @@
#ifndef Must2Scorers_h
#define Must2Scorers_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/PhysicsList.hh b/NPSimulation/include/PhysicsList.hh
index 79ecf3f40f302ea681fafb9e65692a0652cd9454..9c8e3b55d398b21a7cddfc6879d3fcf3f497b20f 100644
--- a/NPSimulation/include/PhysicsList.hh
+++ b/NPSimulation/include/PhysicsList.hh
@@ -1,7 +1,7 @@
#ifndef PhysicsList_h
#define PhysicsList_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/Plastic.hh b/NPSimulation/include/Plastic.hh
index 46f3b24fa967e0deed01cfe842a230783326a32f..02f2020db831720ef16a52c84de051b2b506bf36 100644
--- a/NPSimulation/include/Plastic.hh
+++ b/NPSimulation/include/Plastic.hh
@@ -1,7 +1,7 @@
#ifndef Plastic_h
#define Plastic_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -54,14 +54,24 @@ public:
//////// Specific Function of this Class ///////////
////////////////////////////////////////////////////
public:
- // By Angle Method
- void AddPlastic( G4double R ,
- G4double Theta ,
- G4double Phi ,
- G4double PlasticThickness ,
- G4double PlasticRadius ,
- G4String Scintillator ,
- G4double LeadThickness );
+ // Cylindric plastic
+ void AddPlastic( G4double R ,
+ G4double Theta ,
+ G4double Phi ,
+ G4double PlasticThickness ,
+ G4double PlasticRadius ,
+ G4String Scintillator ,
+ G4double LeadThickness );
+
+ // Squared Plastic
+ void AddPlastic( G4double R ,
+ G4double Theta ,
+ G4double Phi ,
+ G4double Height ,
+ G4double Width ,
+ G4double Thickness ,
+ G4String Scintillator ,
+ G4double LeadThickness );
void VolumeMaker(G4ThreeVector Det_pos, int DetNumber,G4LogicalVolume* world) ;
////////////////////////////////////////////////////
@@ -111,12 +121,16 @@ private:
////////////////////////////////////////////////////
private:
- // if true a Lead plate is added in front or back of the detector
+
+
+ // Lead plate is added in front or back of the detector
vector<double> m_LeadThickness ;
vector<double> m_PlasticThickness ;
- vector<double> m_PlasticRadius ;
-
+ vector<double> m_PlasticRadius ; // cylindrical shape
+ vector<double> m_PlasticHeight ; // squared shape
+ vector<double> m_PlasticWidth ; // squared shape
+
// Used for "By Angle Definition"
vector<G4double> m_R ; // |
vector<G4double> m_Theta ; // > Spherical coordinate plastic volume center
diff --git a/NPSimulation/include/PrimaryGeneratorAction.hh b/NPSimulation/include/PrimaryGeneratorAction.hh
index 97475ec1b592513574e707656d59d8f2c53fe100..9638cea6a00e59c86cd7f1cab1e9870b70c9bdf6 100644
--- a/NPSimulation/include/PrimaryGeneratorAction.hh
+++ b/NPSimulation/include/PrimaryGeneratorAction.hh
@@ -1,7 +1,7 @@
#ifndef PrimaryGeneratorAction_h
#define PrimaryGeneratorAction_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/SteppingVerbose.hh b/NPSimulation/include/SteppingVerbose.hh
index c9ab7c53e5614b6c3810817c8ed213cbbe43b435..fc8ec37052dc529edb38ef0a50b8e61a4a29081f 100644
--- a/NPSimulation/include/SteppingVerbose.hh
+++ b/NPSimulation/include/SteppingVerbose.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/Target.hh b/NPSimulation/include/Target.hh
index 170f7394cb4769e9fc416af1668266d5d88841ba..285b175e29b80b444057b0976ad6d2c623c63f25 100644
--- a/NPSimulation/include/Target.hh
+++ b/NPSimulation/include/Target.hh
@@ -1,7 +1,7 @@
#ifndef Target2Array_h
#define Target2Array_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/ThinSi.hh b/NPSimulation/include/ThinSi.hh
index 48febf6c21ae00357bf7ab0ddbe10d0ca78a48f3..562c0d4d7577850786f28d23791c76bf3cdd7985 100644
--- a/NPSimulation/include/ThinSi.hh
+++ b/NPSimulation/include/ThinSi.hh
@@ -1,7 +1,7 @@
#ifndef ThinSi_h
#define ThinSi_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/ThinSiScorers.hh b/NPSimulation/include/ThinSiScorers.hh
index 9449cffa80c1396dc988e3899c951de766a70bba..a2dc8630b4fdcddcd38f27773435446aa9dd6974 100644
--- a/NPSimulation/include/ThinSiScorers.hh
+++ b/NPSimulation/include/ThinSiScorers.hh
@@ -1,7 +1,7 @@
#ifndef ThinSiScorers_h
#define ThinSiScorers_h 1
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/VDetector.hh b/NPSimulation/include/VDetector.hh
index dd98b31e4954b4c0f6a82c41319822d163db637d..4b05066eb9d6405e7ebe372e210238d115446900 100644
--- a/NPSimulation/include/VDetector.hh
+++ b/NPSimulation/include/VDetector.hh
@@ -1,7 +1,7 @@
#ifndef VDectector_H
#define VDectector_H
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index 59bbd3bca32a6e67823c58929591f519637b554c..b158b6ae8490ef4cefb120166a8605cd39d85004 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -1,7 +1,7 @@
#ifndef VEventGenerator_h
#define VEventGenerator_h
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/AnnularS1.cc b/NPSimulation/src/AnnularS1.cc
index c5fb036b41be93b6215c973a1288b2b59773d4c7..6b3426484162d020a23fb354701b3eca186b079c 100644
--- a/NPSimulation/src/AnnularS1.cc
+++ b/NPSimulation/src/AnnularS1.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/AnnularS1Scorers.cc b/NPSimulation/src/AnnularS1Scorers.cc
index 3e566c8c91a741be8b51fc267d66aafadcc51edd..c8a03e3b7b856c6df8726deafc6be4bbfe1c87b2 100644
--- a/NPSimulation/src/AnnularS1Scorers.cc
+++ b/NPSimulation/src/AnnularS1Scorers.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc
index 391bf2f7e364e189107955159c93a9c5b73312d9..56a612510834c0a6cd604414b8efb1766fbb8897 100644
--- a/NPSimulation/src/DetectorConstruction.cc
+++ b/NPSimulation/src/DetectorConstruction.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/DummyDetector.cc b/NPSimulation/src/DummyDetector.cc
index 0e4e17d8586e596d56a245c48bd8055fb29c3140..bf79671be489264c32293b2edd725569a39cf672 100644
--- a/NPSimulation/src/DummyDetector.cc
+++ b/NPSimulation/src/DummyDetector.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/EventAction.cc b/NPSimulation/src/EventAction.cc
index c370fd39479cd49197bf418c082be2bff02cb6c3..4692e969b094550dcd5ca17ee430f2f62f53e8d4 100644
--- a/NPSimulation/src/EventAction.cc
+++ b/NPSimulation/src/EventAction.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index a8c8db2784d1220bd6c98a6e5744e76339d8e86f..5d263936109277fc965f0ec36e28c2c0aeb9ebda 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -42,7 +42,10 @@ using namespace CLHEP;
EventGeneratorBeam::EventGeneratorBeam()
{
m_InitConditions = new TInitialConditions();
- m_Target = 0;
+ m_Target = NULL ;
+ m_beamA = 0 ;
+ m_beamZ = 0 ;
+ m_particle = NULL ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -57,7 +60,6 @@ void EventGeneratorBeam::SetTarget(Target* Target)
if(Target!=0)
{
m_Target = Target;
- m_Target->WriteDEDXTable(m_particle ,0, m_BeamEnergy);
}
}
@@ -68,8 +70,6 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
string LineBuffer;
string DataBuffer;
- ////////Reaction Setting needs///////
- G4double particleZ = 0 , particleA = 0 ;
//////////////////////////////////////////////////////////////////////////////////////////
ifstream ReactionFile;
ReactionFile.open(Path.c_str());
@@ -92,14 +92,11 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
//Pick-up next line
getline(ReactionFile, LineBuffer);
-
-
if (LineBuffer.compare(0, 4, "Beam") == 0) {
G4cout << "Beam Found" << G4endl ;
ReadingStatus = true ;
}
-
while(ReadingStatus){
ReactionFile >> DataBuffer;
@@ -110,16 +107,22 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 10, "ParticleZ=") == 0) {
check_Z = true ;
ReactionFile >> DataBuffer;
- particleZ = atof(DataBuffer.c_str());
+ m_beamZ = atof(DataBuffer.c_str());
+
+ if(check_A)
+ G4cout << "Beam Particle: Z:" << m_beamZ << " A:" << m_beamA << G4endl;
+
}
else if (DataBuffer.compare(0, 10, "ParticleA=") == 0) {
check_A = true ;
ReactionFile >> DataBuffer;
- particleA = atof(DataBuffer.c_str());
- G4cout << "Beam Particle: Z:" << particleZ << " A:" << particleA << G4endl;
- m_particle = G4ParticleTable::GetParticleTable()->GetIon(particleZ, particleA, 0.);
+ m_beamA = atof(DataBuffer.c_str());
+
+ if(check_Z)
+ G4cout << "Beam Particle: Z:" << m_beamZ << " A:" << m_beamA << G4endl;
+
}
else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
@@ -185,11 +188,15 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
{
- //--------------write the DeDx Table -------------------
- if(m_Target!=0)
- m_Target->WriteDEDXTable(m_particle ,0, m_BeamEnergy+4*m_BeamEnergySpread);
-
-
+ //--------------write the DeDx Table -------------------
+ if( anEvent->GetEventID()==0)
+ {
+ m_particle = G4ParticleTable::GetParticleTable()->GetIon(m_beamZ, m_beamA , 0.);
+
+ if(m_Target!=0 )
+ m_Target->WriteDEDXTable(m_particle ,0, m_BeamEnergy+4*m_BeamEnergySpread);
+ }
+
m_InitConditions->Clear();
///////////////////////////////////////////////////////////////////////
diff --git a/NPSimulation/src/EventGeneratorIsotropic.cc b/NPSimulation/src/EventGeneratorIsotropic.cc
index b2814747a18ab3633cd2e20d5b7071001c05f624..503766e24e3922396bd9a5a85220beeff9f8e4c7 100644
--- a/NPSimulation/src/EventGeneratorIsotropic.cc
+++ b/NPSimulation/src/EventGeneratorIsotropic.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/EventGeneratorPhaseSpace.cc b/NPSimulation/src/EventGeneratorPhaseSpace.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ed612b9adef5a5aeec3d7f4077247ca93788a707
--- /dev/null
+++ b/NPSimulation/src/EventGeneratorPhaseSpace.cc
@@ -0,0 +1,437 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : Feb 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated pure phase space event and *
+ * especially evaluate the background contribution of different phase space *
+ * channel during an experiment. Manage up to 18 body phase space *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <iostream>
+#include <fstream>
+#include <limits>
+
+// G4 header defining G4 types
+#include "globals.hh"
+
+// G4 headers
+#include "G4ParticleTable.hh"
+#include "G4ParticleGun.hh"
+#include "G4RotationMatrix.hh"
+
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// NPTool headers
+#include "EventGeneratorPhaseSpace.hh"
+#include "RootOutput.h"
+
+//Root Headers
+#include "TGenPhaseSpace.h"
+
+using namespace std;
+using namespace CLHEP;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorPhaseSpace::EventGeneratorPhaseSpace()
+{
+ //------------- Default Constructor -------------
+ m_InitConditions = new TInitialConditions() ;
+ m_Target = new Target() ;
+ m_SigmaX = 0 ;
+ m_SigmaY = 0 ;
+ m_SigmaThetaX = 0 ;
+ m_SigmaPhiY = 0 ;
+ m_EventWeight = 0 ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorPhaseSpace::~EventGeneratorPhaseSpace()
+{
+ //------------- Default Destructor ------------
+ delete m_InitConditions ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorPhaseSpace::SetTarget(Target* Target)
+ {
+ if(Target!=0)
+ m_Target = Target;
+ }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorPhaseSpace::EventGeneratorPhaseSpace( G4ParticleDefinition* BeamParticle ,
+ G4ParticleDefinition* TargetParticle ,
+ vector<G4ParticleDefinition*> ReactionProducts ,
+ double BeamEnergy ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY )
+{
+ //------------- Constructor with nuclei names and beam energy ------------
+
+ SetEverything( BeamParticle ,
+ TargetParticle ,
+ ReactionProducts ,
+ BeamEnergy ,
+ BeamEnergySpread ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY );
+
+ m_EventWeight = 0 ;
+
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorPhaseSpace::InitializeRootOutput()
+{
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
+ pTree->Branch("EventWeight",&m_EventWeight,"EventWeigt/D");
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Inherit from VEventGenerator
+
+void EventGeneratorPhaseSpace::ReadConfiguration(string Path)
+{
+////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+
+////////Reaction Setting needs///////
+ string Beam, Target, Heavy, Light, CrossSectionPath ;
+ G4double BeamEnergy = 0 , BeamEnergySpread = 0 , SigmaX = 0 , SigmaY = 0 , SigmaThetaX = 0 , SigmaPhiY=0;
+ vector<G4ParticleDefinition*> ReactionProducts;
+ G4ParticleDefinition* BeamParticle = NULL ;
+ G4ParticleDefinition* TargetParticle = NULL ;
+
+ bool ReadingStatus = false ;
+ bool check_Beam = false ;
+ bool check_Target = false ;
+ bool check_ReactionProducts = false ;
+ bool check_BeamEnergy = false ;
+ bool check_BeamEnergySpread = false ;
+ bool check_FWHMX = false ;
+ bool check_FWHMY = false ;
+ bool check_EmmitanceTheta = false ;
+ bool check_EmmitancePhi = false ;
+
+//////////////////////////////////////////////////////////////////////////////////////////
+ ifstream ReactionFile;
+ ReactionFile.open(Path.c_str());
+
+ if (ReactionFile.is_open()) {}
+ else {
+ return;
+ }
+
+ while (!ReactionFile.eof()) {
+ //Pick-up next line
+ getline(ReactionFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 10, "PhaseSpace") == 0) { ReadingStatus = true ;}
+
+
+while(ReadingStatus){
+
+ ReactionFile >> DataBuffer;
+
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ else if (DataBuffer.compare(0, 5, "Beam=") == 0) {
+ check_Beam = true ;
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int A = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int Z = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int EXX = atoi(DataBuffer.c_str())*MeV;
+
+ G4ParticleDefinition* ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
+ BeamParticle = ParticleBuffer ;
+
+ G4cout << "Beam : A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "Target=") == 0) {
+ check_Target = true ;
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int A = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int Z = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int EXX = atoi(DataBuffer.c_str())*MeV;
+
+ G4ParticleDefinition* ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
+ TargetParticle = ParticleBuffer ;
+
+ G4cout << "Target : A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
+ check_BeamEnergy = true ;
+ ReactionFile >> DataBuffer;
+ BeamEnergy = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "Beam Energy " << BeamEnergy / MeV << " MeV" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 17, "BeamEnergySpread=") == 0) {
+ check_BeamEnergySpread = true ;
+ ReactionFile >> DataBuffer;
+ BeamEnergySpread = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
+ check_FWHMX = true ;
+ ReactionFile >> DataBuffer;
+ SigmaX = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
+ check_FWHMY = true ;
+ ReactionFile >> DataBuffer;
+ SigmaY = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 19, "SigmaThetaX=") == 0) {
+ check_EmmitanceTheta = true ;
+ ReactionFile >> DataBuffer;
+ SigmaThetaX = atof(DataBuffer.c_str()) * deg;
+ G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 17, "SigmaPhiY=") == 0) {
+ check_EmmitancePhi = true ;
+ ReactionFile >> DataBuffer;
+ SigmaPhiY = atof(DataBuffer.c_str()) * deg;
+ G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 13, "DecayProduct=") == 0) {
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int A = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int Z = atoi(DataBuffer.c_str());
+ ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
+ int EXX = atoi(DataBuffer.c_str())*MeV;
+
+ G4ParticleDefinition* ParticleBuffer ;
+
+ if ( A==1 && Z==0 )
+ ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+ else if ( A==1 && Z==1)
+ ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("proton");
+// else if ( A==2 && Z==0)
+// ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+// ParticleBuffer->SetPDGMass( ParticleBuffer->GetPDGMass*2);
+ else
+ ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
+
+ ReactionProducts.push_back(ParticleBuffer);
+
+ G4cout << "Decay Product: A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
+ }
+
+ else if (DataBuffer.compare(0, 21, "EndOfDecayProductList") == 0) {
+ check_ReactionProducts = true ;
+ }
+
+
+ ///////////////////////////////////////////////////
+ // If no Transfert Token and no comment, toggle out
+ else
+ {ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
+
+ ///////////////////////////////////////////////////
+ // If all Token found toggle out
+ if( check_Beam && check_Target && check_BeamEnergy && check_BeamEnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta
+ && check_EmmitancePhi && check_ReactionProducts)
+ ReadingStatus = false ;
+
+ }
+
+
+ }
+
+ SetEverything( BeamParticle ,
+ TargetParticle ,
+ ReactionProducts ,
+ BeamEnergy ,
+ BeamEnergySpread ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY );
+
+ ReactionFile.close();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorPhaseSpace::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
+{
+
+ // If first time, write the DeDx table
+ if(anEvent->GetEventID()==0 && m_Target!=0)
+ {
+ m_Target->WriteDEDXTable(m_BeamParticle ,0, m_BeamEnergy+4*m_BeamEnergySpread);
+ for(unsigned int k = 0 ; k < m_ReactionProducts.size() ; k++)
+ {
+ m_Target->WriteDEDXTable(m_ReactionProducts[k] ,0, m_BeamEnergy+4*m_BeamEnergySpread);
+ }
+ }
+
+ // Clear contents of Precedent event (now stored in ROOTOutput)
+ m_InitConditions->Clear() ;
+ m_EventWeight = 0 ;
+
+ ///////////////////////////////////////////////////////////////////////
+ ///// Calculate the incident beam direction as well as the vertex /////
+ ///// of interaction in target and Energy Loss of the beam within /////
+ ///// the target. /////
+ ///////////////////////////////////////////////////////////////////////
+ G4ThreeVector InterCoord;
+
+ G4double Beam_thetaX = 0, Beam_phiY = 0;
+ G4double Beam_theta = 0, Beam_phi = 0;
+ G4double FinalBeamEnergy = 0 ;
+ G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
+
+ m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ InitialBeamEnergy,
+ m_BeamParticle,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi,
+ FinalBeamEnergy);
+
+ m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
+
+ // write vertex position to ROOT file
+ G4double x0 = InterCoord.x();
+ G4double y0 = InterCoord.y();
+ G4double z0 = InterCoord.z();
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
+
+ // write emittance angles to ROOT file
+ m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
+ m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
+
+ // write angles to ROOT file
+ m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
+ m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
+
+ //////////////////////////////////////////////////////////
+ //////////////// Phase Space Calculation /////////////////
+ //////////////////////////////////////////////////////////
+
+ // Masses array of ecay products in GeV (unit used by ROOT)
+ int NumberOfReactionProducts = m_ReactionProducts.size() ;
+ double* masses = new double[NumberOfReactionProducts];
+
+ for(int k = 0 ; k < NumberOfReactionProducts ; k++ )
+ {
+ masses[k] = m_ReactionProducts[k]-> GetPDGMass()/GeV ;
+ }
+
+ // Kinematics of reaction
+ G4double M = m_BeamParticle -> GetPDGMass() ;
+ G4double InitialE = FinalBeamEnergy + M ;
+ G4double InitialMomentumX = sqrt( InitialE*InitialE - M*M) * sin(Beam_theta) * cos(Beam_phi) ;
+ G4double InitialMomentumY = sqrt( InitialE*InitialE - M*M) * sin(Beam_theta) * sin(Beam_phi) ;
+ G4double InitialMomentumZ = sqrt( InitialE*InitialE - M*M) * cos(Beam_theta) ;
+
+ TLorentzVector Initial = TLorentzVector(InitialMomentumX/GeV, InitialMomentumY/GeV, InitialMomentumZ/GeV,InitialE/GeV) + TLorentzVector(0,0,0,m_TargetParticle -> GetPDGMass() / GeV);
+
+
+ // Instentiate a Phase Space Generator, with flat distrution
+ TGenPhaseSpace TPhaseSpace ;
+
+ if( !TPhaseSpace.SetDecay(Initial, NumberOfReactionProducts , masses) ) cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles "<<endl;
+
+ // Generate event and store the associate weight. Think to use this weigt to get correcte spectrum
+ m_EventWeight = TPhaseSpace.Generate() ;
+
+
+ TLorentzVector* ProductLV ;
+ for ( int u = 0; u < NumberOfReactionProducts ; u++)
+ {
+ ProductLV = TPhaseSpace.GetDecay(u);
+
+ G4ThreeVector Momentum ( ProductLV->X()*GeV ,
+ ProductLV->Y()*GeV ,
+ ProductLV->Z()*GeV );
+ Momentum.unit() ;
+
+ G4double Energy = ProductLV->E()*GeV-masses[u]*GeV ;
+
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(m_ReactionProducts[u]) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+
+ // get theta and phi in the world frame
+ G4double theta_world = Momentum.theta();
+ G4double phi_world = Momentum.phi();
+ if (phi_world < 1e-6) phi_world += 2*pi;
+ // write angles in ROOT file
+ m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
+ m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
+ m_InitConditions->SetICEmittedEnergy(Energy);
+
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorPhaseSpace::SetEverything( G4ParticleDefinition* BeamParticle ,
+ G4ParticleDefinition* TargetParticle ,
+ vector<G4ParticleDefinition*> ReactionProducts ,
+ double BeamEnergy ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY )
+{
+ //------------- Constructor with nuclei names and beam energy ------------
+ m_BeamParticle = BeamParticle ;
+ m_TargetParticle = TargetParticle ;
+ m_ReactionProducts = ReactionProducts ;
+ m_BeamEnergy = BeamEnergy ;
+ m_BeamEnergySpread = BeamEnergySpread ;
+ m_SigmaX = SigmaX ;
+ m_SigmaY = SigmaY ;
+ m_SigmaThetaX = SigmaThetaX ;
+ m_SigmaPhiY = SigmaPhiY ;
+}
+
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 0929619972a6176532cc38921cbb9188ded31200..1e251c2a2378a3d0585d49b99a596e303e008d7a 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index dab91dec0c35eed7f27f43f8df34e9893adca48a..f7338b3a15f0c470c67fe3d69df7e4e45ae31cea 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -46,32 +46,34 @@
//Root Headers
#include "TGenPhaseSpace.h"
+//CLHEP
+#include "CLHEP/Random/RandBreitWigner.h"
+
using namespace std;
using namespace CLHEP;
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance()
{
//------------- Default Constructor -------------
m_InitConditions = new TInitialConditions() ;
- m_Reaction = new Reaction() ;
- m_Target = new Target();
- m_SigmaX = 0 ;
- m_SigmaY = 0 ;
- m_SigmaThetaX = 0 ;
- m_SigmaPhiY = 0 ;
- m_ResonanceDecayZ = 0 ;
- m_ResonanceDecayA = 0 ;
+ m_Reaction = new Reaction() ;
+ m_Target = new Target() ;
+ m_SigmaX = 0 ;
+ m_SigmaY = 0 ;
+ m_SigmaThetaX = 0 ;
+ m_SigmaPhiY = 0 ;
+ m_ResonanceDecayZ = 0 ;
+ m_ResonanceDecayA = 0 ;
+ m_EventWeight = 0 ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance()
{
//------------- Default Destructor ------------
- delete m_InitConditions;
- delete m_Reaction ;
+ delete m_InitConditions ;
+ delete m_Reaction ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfertToResonance::SetTarget(Target* Target)
@@ -119,6 +121,8 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string
ShootHeavy ,
ShootDecayProduct ,
Path );
+
+ m_EventWeight = 0 ;
}
@@ -129,13 +133,7 @@ void EventGeneratorTransfertToResonance::InitializeRootOutput()
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::Print() const
-{
-
+ pTree->Branch("EventWeight",&m_EventWeight,"EventWeigt/D");
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -188,7 +186,7 @@ void EventGeneratorTransfertToResonance::ReadConfiguration(string Path)
- if (LineBuffer.compare(0, 9, "Transfert") == 0) { ReadingStatus = true ;}
+ if (LineBuffer.compare(0, 20, "TransfertToResonance") == 0) { ReadingStatus = true ;}
while(ReadingStatus){
@@ -345,24 +343,24 @@ while(ReadingStatus){
}
- SetEverything( Beam ,
- Target ,
- Light ,
- Heavy ,
- BeamEnergy ,
- ExcitationEnergy ,
- BeamEnergySpread ,
- SigmaX ,
- SigmaY ,
- SigmaThetaX ,
- SigmaPhiY ,
- ResonanceWidth,
- ResonanceDecayZ ,
- ResonanceDecayA ,
- ShootLight ,
- ShootHeavy ,
- ShootDecayProduct ,
- CrossSectionPath );
+ SetEverything( Beam ,
+ Target ,
+ Light ,
+ Heavy ,
+ BeamEnergy ,
+ ExcitationEnergy ,
+ BeamEnergySpread ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY ,
+ ResonanceWidth ,
+ ResonanceDecayZ ,
+ ResonanceDecayA ,
+ ShootLight ,
+ ShootHeavy ,
+ ShootDecayProduct ,
+ CrossSectionPath );
ReactionFile.close();
}
@@ -370,7 +368,9 @@ while(ReadingStatus){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
{
-
+ // Initialize event weight to one.
+ m_EventWeight = 1 ;
+
// If first time, write the DeDx table
if(anEvent->GetEventID()==0)
{
@@ -403,26 +403,23 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
G4ParticleDefinition* LightName
= G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, 0.);
- // Recoil
- G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
- G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
-
- G4ParticleDefinition* HeavyName
- = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
-
+ // Shoot the Resonance energy following the mean and width value
+ // EXX should always be more than specific heat of the reaction
+ // double EXX = RandBreitWigner::shoot(m_ResonanceMean,m_ResonanceWidth) ;
+ double EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
+ m_Reaction->SetExcitation( EXX );
+
+ while ( m_Reaction->CheckKinematic()==false )
+ {
+// EXX = RandBreitWigner::shoot(m_ResonanceMean,m_ResonanceWidth) ;
+ EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
+ m_Reaction->SetExcitation( EXX );
+ }
// Beam
G4int BeamZ = m_Reaction->GetNucleus1()->GetZ();
G4int BeamA = m_Reaction->GetNucleus1()->GetA();
G4ParticleDefinition* BeamName = G4ParticleTable::GetParticleTable()->GetIon(BeamZ, BeamA, 0);
- // Shoot the Resonance energy following the mean and width value
- double EXX = -10 ;
-
- while(EXX<0)
- EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
-
- m_Reaction->SetExcitation( EXX );
-
///////////////////////////////////////////////////////////////////////
///// Calculate the incident beam direction as well as the vertex /////
///// of interaction in target and Energy Loss of the beam within /////
@@ -482,11 +479,6 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
///// Angles for emitted particles following Cross Section //////
///// Angles are in the beam frame //////
/////////////////////////////////////////////////////////////////
- // Beam incident energy
- G4double NominalBeamEnergy = m_BeamEnergy;
- G4double IncidentBeamEnergy = RandGauss::shoot(NominalBeamEnergy, m_BeamEnergySpread);
- m_Reaction->SetBeamEnergy(IncidentBeamEnergy);
- m_InitConditions->SetICIncidentEnergy(IncidentBeamEnergy / MeV);
// Angles
RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
G4double ThetaCM = CrossSectionShoot.shoot() * (180*deg);
@@ -545,12 +537,6 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
}
// Case of recoil particle
- /* // Particle type
- particleGun->SetParticleDefinition(HeavyName);
- // Particle energy
- particleGun->SetParticleEnergy(EnergyHeavy);
- // Particle vertex position
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
// Particle direction
// Kinematical angles in the beam frame are transformed
// to the "world" frame*/
@@ -560,14 +546,15 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
G4double phi_world = momentum_kine_world.phi();
if (phi_world < 1e-6) phi_world += 2*pi;
+ if(m_ShootHeavy || m_ShootDecayProduct)
EventGeneratorTransfertToResonance::ResonanceDecay( EnergyHeavy ,
- theta_world ,
- phi_world ,
- x0 ,
- y0 ,
- z0 ,
- anEvent ,
- particleGun );
+ theta_world ,
+ phi_world ,
+ x0 ,
+ y0 ,
+ z0 ,
+ anEvent ,
+ particleGun );
}
@@ -587,8 +574,8 @@ void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy
G4double parentA = m_Reaction->GetNucleus4()->GetA() ;
G4int NumberOfNeutrons = (parentA - parentZ) - (m_ResonanceDecayA - m_ResonanceDecayZ) ;
- G4int NumberOfProtons = parentZ - m_ResonanceDecayZ ;
- G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;
+ G4int NumberOfProtons = parentZ - m_ResonanceDecayZ ;
+ G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;
if (NumberOfNeutrons < 0 || NumberOfProtons < 0) {
G4cout << "Error input for Resonance decay" << G4endl;
@@ -598,7 +585,7 @@ void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy
else {
//Obtain Mass of daughter Nuclei
G4ParticleDefinition* parent
- = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation()) ;
+ = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, 0 ) ;
G4ParticleDefinition* daughter
= G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
@@ -609,11 +596,20 @@ void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy
G4ParticleDefinition* proton
= G4ParticleTable::GetParticleTable()->FindParticle("proton");
- G4double M = parent -> GetPDGMass() ;
+ G4double M = parent -> GetPDGMass() + m_Reaction->GetExcitation()*MeV ;
G4double md = daughter -> GetPDGMass() ;
G4double mn = neutron -> GetPDGMass() ;
G4double mp = proton -> GetPDGMass() ;
+ // Check that we are above threshold:
+ // If the Resonnance go below the threshold, decay is forced at thereshold
+ if (M < md + NumberOfNeutrons*mn + NumberOfProtons*mp)
+ {
+ double NewExx = parent -> GetPDGMass() - md - NumberOfNeutrons*mn - NumberOfProtons*mp ;
+ M = parent -> GetPDGMass() + NewExx ;
+ }
+
+
G4double InitialE = EnergyHeavy + M ;
G4double InitialMomentumX = sqrt( InitialE*InitialE - M*M) * sin(ThetaHeavy) * cos(PhiHeavy) ;
G4double InitialMomentumY = sqrt( InitialE*InitialE - M*M) * sin(ThetaHeavy) * sin(PhiHeavy) ;
@@ -638,15 +634,9 @@ void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy
TGenPhaseSpace TPhaseSpace ;
if( !TPhaseSpace.SetDecay(Initial, NumberOfDecayProducts+1, masses) ) cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles "<<endl;
- double MaxWt=TPhaseSpace.GetWtMax() ;
- double Weight = 0 ;
- double Rand = 1 ;
-
- while( Rand > Weight )
- {
- Weight = TPhaseSpace.Generate() ;
- Rand = CLHEP::RandFlat::shoot()*MaxWt ;
- }
+
+ // Generate an event and store is weight in the Output Tree
+ m_EventWeight = TPhaseSpace.Generate() ;
TLorentzVector* daugterLV ;
diff --git a/NPSimulation/src/GaspardScorers.cc b/NPSimulation/src/GaspardScorers.cc
index 266edcfdfdb814c35fe5536c4f7e4d29639634b0..3897c87ebd874a10c61c58505aa7fcc8b7750860 100644
--- a/NPSimulation/src/GaspardScorers.cc
+++ b/NPSimulation/src/GaspardScorers.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTracker.cc b/NPSimulation/src/GaspardTracker.cc
index c822520055818fe99167a955f3ad8c184284ae8e..bd377f355efba4bfea1b508ac7aecfe954b2ec60 100644
--- a/NPSimulation/src/GaspardTracker.cc
+++ b/NPSimulation/src/GaspardTracker.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTrackerAnnular.cc b/NPSimulation/src/GaspardTrackerAnnular.cc
index 9e00cbe189133cd2c8b5bd81e467555bf3fa00cb..b6014d4d91620d02e7ae779b248432e91bdb000e 100644
--- a/NPSimulation/src/GaspardTrackerAnnular.cc
+++ b/NPSimulation/src/GaspardTrackerAnnular.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTrackerDummyShape.cc b/NPSimulation/src/GaspardTrackerDummyShape.cc
index 20b87009b32667c511016f36c0a1ea139a89c5cc..c604ce37f0e879bcb5b68f1794860f426f090569 100644
--- a/NPSimulation/src/GaspardTrackerDummyShape.cc
+++ b/NPSimulation/src/GaspardTrackerDummyShape.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTrackerModule.cc b/NPSimulation/src/GaspardTrackerModule.cc
index 890cd77702c06e4720d4ca5b54e55f513e67a61e..4fb280990d0bfcaae6913815f275046705972902 100644
--- a/NPSimulation/src/GaspardTrackerModule.cc
+++ b/NPSimulation/src/GaspardTrackerModule.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTrackerSquare.cc b/NPSimulation/src/GaspardTrackerSquare.cc
index b9eee508a5fdb67749a45c38b148e716b16b066c..5fbeb423a30c59a3550e82cee0a485547f59e547 100644
--- a/NPSimulation/src/GaspardTrackerSquare.cc
+++ b/NPSimulation/src/GaspardTrackerSquare.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GaspardTrackerTrapezoid.cc b/NPSimulation/src/GaspardTrackerTrapezoid.cc
index 2a38ba8c373085604a94e8fc14af6140fdc87ad3..a0b6abe10806a849275bf7efa84832b57ab9afef 100644
--- a/NPSimulation/src/GaspardTrackerTrapezoid.cc
+++ b/NPSimulation/src/GaspardTrackerTrapezoid.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/GeneralScorers.cc b/NPSimulation/src/GeneralScorers.cc
index c22470385ec78a2e41f3deaac458714dfd2bd302..50e8c41d047b0b0a59666ddc973ce85287ac94ef 100644
--- a/NPSimulation/src/GeneralScorers.cc
+++ b/NPSimulation/src/GeneralScorers.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/MUST2Array.cc b/NPSimulation/src/MUST2Array.cc
index 8b18d4b053b8d28e568af766f7d27e05afb65bbd..b8c5f43281a739a3ab1ec0d4bc16fcc661fbe14a 100644
--- a/NPSimulation/src/MUST2Array.cc
+++ b/NPSimulation/src/MUST2Array.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -163,12 +163,12 @@ void MUST2Array::VolumeMaker(G4int TelescopeNumber ,
G4String Name = "MUST2Telescope" + DetectorNumber ;
- PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos) ,
- logicMM ,
- Name ,
- world ,
- false ,
- 0);
+ PVPBuffer = new G4PVPlacement( G4Transform3D(*MMrot, MMpos) ,
+ logicMM ,
+ Name ,
+ world ,
+ false ,
+ 0);
@@ -879,21 +879,19 @@ void MUST2Array::ConstructDetector(G4LogicalVolume* world)
// (u,v) // to silicon plan
// w perpendicular to (u,v) plan and pointing CsI
MMu = m_X128_Y1[i] - m_X1_Y1[i] ;
- MMu = -MMu.unit() ;
+ MMu = MMu.unit() ;
MMv = m_X1_Y128[i] - m_X1_Y1[i] ;
MMv = MMv.unit() ;
-
- MMw = MMu.cross(MMv) ;
+
+ MMw = MMv.cross(MMu) ;
// if (MMw.z() > 0)MMw = MMv.cross(MMu) ;
MMw = MMw.unit() ;
-
+
MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ;
// Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // translation to place Telescope
+ MMrot = new G4RotationMatrix(MMv, MMu, MMw);
MMpos = MMw * Length * 0.5 + MMCenter ;
}
diff --git a/NPSimulation/src/Must2Scorers.cc b/NPSimulation/src/Must2Scorers.cc
index 2b057adf11486a6508e9239bc808b555e870b688..1d9d6bc813ab1b3a4b362744fefb6c2b6f525d20 100644
--- a/NPSimulation/src/Must2Scorers.cc
+++ b/NPSimulation/src/Must2Scorers.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/PhysicsList.cc b/NPSimulation/src/PhysicsList.cc
index 104f0673950be17dba416e4ab23c08eb5a4eef35..0bb173f6cdb65ef72f24e41b2e5afc9c3aea5ebd 100644
--- a/NPSimulation/src/PhysicsList.cc
+++ b/NPSimulation/src/PhysicsList.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/Plastic.cc b/NPSimulation/src/Plastic.cc
index 5e1e3fda47083a810e4cbdd617b9808a6d135c24..a8ea8618837d62f6d511dbd9b013a29da78da7c5 100644
--- a/NPSimulation/src/Plastic.cc
+++ b/NPSimulation/src/Plastic.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -27,6 +27,7 @@
#include <limits>
//G4 Geometry object
#include "G4Tubs.hh"
+#include "G4Box.hh"
//G4 sensitive
#include "G4SDManager.hh"
@@ -91,16 +92,40 @@ void Plastic::AddPlastic( G4double R ,
G4double PlasticRadius ,
G4String Scintillator ,
G4double LeadThickness )
-{
+ {
+
+ m_R.push_back(R) ;
+ m_Theta.push_back(Theta) ;
+ m_Phi.push_back(Phi) ;
+ m_PlasticThickness.push_back(PlasticThickness) ;
+ m_LeadThickness.push_back(LeadThickness) ;
+ m_Scintillator.push_back(Scintillator) ;
+ m_PlasticRadius.push_back(PlasticRadius) ; // cylindrical shape
+ m_PlasticHeight.push_back(-1) ; // squared shape
+ m_PlasticWidth.push_back(-1) ; // squared shape
+ }
+
+void Plastic::AddPlastic( G4double R ,
+ G4double Theta ,
+ G4double Phi ,
+ G4double Height ,
+ G4double Width ,
+ G4double PlasticThickness ,
+ G4String Scintillator ,
+ G4double LeadThickness )
+ {
+ m_R.push_back(R) ;
+ m_Theta.push_back(Theta) ;
+ m_Phi.push_back(Phi) ;
+ m_PlasticThickness.push_back(PlasticThickness) ;
+ m_LeadThickness.push_back(LeadThickness) ;
+ m_Scintillator.push_back(Scintillator) ;
+ m_PlasticRadius.push_back(-1) ; // cylindrical shape
+ m_PlasticHeight.push_back(Height) ; // squared shape
+ m_PlasticWidth.push_back(Width) ; // squared shape
+
+ }
- m_R.push_back(R) ;
- m_Theta.push_back(Theta) ;
- m_Phi.push_back(Phi) ;
- m_PlasticThickness.push_back(PlasticThickness) ;
- m_PlasticRadius.push_back(PlasticRadius) ;
- m_LeadThickness.push_back(LeadThickness) ;
- m_Scintillator.push_back(Scintillator) ;
-}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -117,18 +142,24 @@ void Plastic::ReadConfiguration(string Path)
string LineBuffer ;
string DataBuffer ;
- G4double Theta = 0 , Phi = 0 , R = 0 , Thickness = 0 , Radius = 0 , LeadThickness = 0;
- G4String Scintillator ;
+ G4double Theta = 0 , Phi = 0 , R = 0 , Thickness = 0 , Radius = 0 , LeadThickness = 0, X = 0 , Y = 0 , Z = 0 , Width = 0 , Height = 0 ;
+ G4String Scintillator, Shape ;
- bool check_Theta = false ;
- bool check_Phi = false ;
- bool check_R = false ;
+ bool check_Theta = false ;
+ bool check_Phi = false ;
+ bool check_R = false ;
bool check_Thickness = false ;
- bool check_Radius = false ;
- bool check_LeadThickness = false ;
+ bool check_Radius = false ;
+ bool check_LeadThickness = false ;
bool check_Scintillator = false ;
- bool ReadingStatus = false ;
-
+ bool check_Height = false ;
+ bool check_Width = false ;
+ bool check_Shape = false ;
+ bool check_X = false ;
+ bool check_Y = false ;
+ bool check_Z = false ;
+ bool ReadingStatus = false ;
+
while (!ConfigFile.eof())
{
@@ -187,6 +218,41 @@ void Plastic::ReadConfiguration(string Path)
cout << "R: " << R/mm << endl;
}
+ //Position method
+ else if (DataBuffer.compare(0, 2, "X=") == 0) {
+ check_X = true;
+ ConfigFile >> DataBuffer ;
+ X = atof(DataBuffer.c_str()) ;
+ X = X * mm;
+ cout << "X: " << X / mm << endl;
+ }
+
+ else if (DataBuffer.compare(0, 2, "Y=") == 0) {
+ check_Y = true;
+ ConfigFile >> DataBuffer ;
+ Y = atof(DataBuffer.c_str()) ;
+ Y = Y * mm;
+ cout << "Y: " << Y / mm << endl;
+ }
+
+ else if (DataBuffer.compare(0, 2, "Z=") == 0) {
+ check_Z = true;
+ ConfigFile >> DataBuffer ;
+ Z = atof(DataBuffer.c_str()) ;
+ Z = Z * mm;
+ cout << "Z: " << Z / mm << endl;
+ }
+
+
+ //General
+ else if (DataBuffer.compare(0, 6, "Shape=") == 0) {
+ check_Shape = true;
+ ConfigFile >> DataBuffer ;
+ Shape = DataBuffer ;
+ cout << "Shape: " << Shape << endl;
+ }
+
+ // Cylindrical shape
else if (DataBuffer.compare(0, 7, "Radius=") == 0) {
check_Radius = true;
ConfigFile >> DataBuffer ;
@@ -195,6 +261,24 @@ void Plastic::ReadConfiguration(string Path)
cout << "Plastic Radius: " << Radius/mm << endl;
}
+ // Squared shape
+ else if (DataBuffer.compare(0, 7, "Width=") == 0) {
+ check_Width = true;
+ ConfigFile >> DataBuffer ;
+ Width = atof(DataBuffer.c_str()) ;
+ Width = Width * mm;
+ cout << "Plastic Width: " << Width/mm << endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "Height=") == 0) {
+ check_Height = true;
+ ConfigFile >> DataBuffer ;
+ Height = atof(DataBuffer.c_str()) ;
+ Height = Height * mm;
+ cout << "Plastic Height: " << Height/mm << endl;
+ }
+
+ // Common
else if (DataBuffer.compare(0, 10, "Thickness=") == 0) {
check_Thickness = true;
ConfigFile >> DataBuffer ;
@@ -226,27 +310,57 @@ void Plastic::ReadConfiguration(string Path)
/////////////////////////////////////////////////
// If All necessary information there, toggle out
- if ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator)
+ if ( ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator && check_Shape) // Cylindrical case
+ || ( check_X && check_Y && check_Z && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator )
+
+ || ( check_Theta && check_Phi && check_R && check_Thickness && check_Width && check_Height && check_LeadThickness && check_Scintillator && check_Shape ) // Squared case
+ || ( check_X && check_Y && check_Z && check_Thickness && check_Width && check_Height && check_LeadThickness && check_Scintillator )
+ )
{
- AddPlastic( R ,
- Theta ,
- Phi ,
- Thickness ,
- Radius ,
- Scintillator ,
- LeadThickness );
+
+ if (check_X && check_Y && check_Z)
+ {
+ R = sqrt (X*X+Y*Y+Z*Z) ;
+ Theta = acos(Z / (R) ) ;
+ Phi = atan2(Y,X) ;
+ }
+
+ if (Shape == "Cylinder")
+ AddPlastic( R ,
+ Theta ,
+ Phi ,
+ Thickness ,
+ Radius ,
+ Scintillator ,
+ LeadThickness );
+
+ else if (Shape == "Square")
+ AddPlastic( R ,
+ Theta ,
+ Phi ,
+ Height ,
+ Width ,
+ Thickness ,
+ Scintillator ,
+ LeadThickness );
// Reinitialisation of Check Boolean
check_Theta = false ;
- check_Phi = false ;
+ check_Phi = false ;
check_R = false ;
check_Thickness = false ;
check_Radius = false ;
- check_LeadThickness = false ;
- check_Scintillator = false ;
+ check_LeadThickness = false ;
+ check_Scintillator = false ;
+ check_Height = false ;
+ check_Width = false ;
+ check_Shape = false ;
+ check_X = false ;
+ check_Y = false ;
+ check_Z = false ;
ReadingStatus = false ;
- cout << "///"<< endl ;
+ cout << "///"<< endl ;
}
}
@@ -303,58 +417,103 @@ void Plastic::VolumeMaker(G4ThreeVector Det_pos, int DetNumber, G4LogicalVolume*
// Definition of the volume containing the sensitive detector
- if(m_PlasticThickness[i]>0 && m_PlasticRadius[i]>0)
- {
- G4Tubs* solidPlastic = new G4Tubs( Name ,
- 0 ,
- m_PlasticRadius[i] ,
- m_PlasticThickness[i]/2 ,
- 0*deg ,
- 360*deg );
-
- G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0);
- logicPlastic->SetSensitiveDetector(m_PlasticScorer);
+
+
+ // Cylindrical Case
+ if(m_PlasticRadius[i]!=-1)
+ {
+ if(m_PlasticThickness[i]>0 && m_PlasticRadius[i]>0)
+ {
+ G4Tubs* solidPlastic = new G4Tubs( Name ,
+ 0 ,
+ m_PlasticRadius[i] ,
+ m_PlasticThickness[i]/2 ,
+ 0*deg ,
+ 360*deg );
+
+ G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0);
+ logicPlastic->SetSensitiveDetector(m_PlasticScorer);
+
+ G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ;
+ logicPlastic->SetVisAttributes(PlastVisAtt) ;
+
+
+
+ PVPBuffer = new G4PVPlacement( 0 ,
+ Det_pos ,
+ logicPlastic ,
+ Name + "_Scintillator" ,
+ world ,
+ false ,
+ 0 );
+ }
- G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ;
- logicPlastic->SetVisAttributes(PlastVisAtt) ;
-
-
-
- PVPBuffer = new G4PVPlacement( 0 ,
- Det_pos ,
- logicPlastic ,
- Name + "_Scintillator" ,
- world ,
- false ,
- 0 );
-
-
+ if(m_LeadThickness[i]>0&& m_PlasticRadius[i]>0)
+ {
+ G4Tubs* solidLead = new G4Tubs( Name+"_Lead" ,
+ 0 ,
+ m_PlasticRadius[i] ,
+ m_LeadThickness[i]/2 ,
+ 0*deg ,
+ 360*deg );
+
+ G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0);
+ G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ;
+ logicLead->SetVisAttributes(LeadVisAtt) ;
+
+ PVPBuffer = new G4PVPlacement( 0 ,
+ Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() ,
+ logicLead ,
+ Name+"_Lead" ,
+ world ,
+ false ,
+ 0 );
+
+ }
+ }
+
+ // Squared case
+ if(m_PlasticHeight[i]!=-1)
+ {
+
+ if(m_PlasticThickness[i]>0 && m_PlasticHeight[i]>0 && m_PlasticWidth[i]>0)
+ {
+ G4Box* solidPlastic = new G4Box(Name, 0.5*m_PlasticWidth[i], 0.5*m_PlasticHeight[i], 0.5*m_PlasticThickness[i]);
+ G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0);
+ logicPlastic->SetSensitiveDetector(m_PlasticScorer);
+
+ G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ;
+ logicPlastic->SetVisAttributes(PlastVisAtt) ;
+
+ PVPBuffer = new G4PVPlacement( 0 ,
+ Det_pos ,
+ logicPlastic ,
+ Name + "_Scintillator" ,
+ world ,
+ false ,
+ 0 );
+ }
+
+ if(m_LeadThickness[i]>0&& m_PlasticHeight[i]>0 && m_PlasticWidth[i]>0)
+ {
+ G4Box* solidLead = new G4Box(Name+"_Lead", 0.5*m_PlasticWidth[i], 0.5*m_PlasticHeight[i], 0.5*m_LeadThickness[i]);
+
+ G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0);
+ G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ;
+ logicLead->SetVisAttributes(LeadVisAtt) ;
+
+ PVPBuffer = new G4PVPlacement( 0 ,
+ Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() ,
+ logicLead ,
+ Name+"_Lead" ,
+ world ,
+ false ,
+ 0 );
+
+ }
}
-
- if(m_LeadThickness[i]>0&& m_PlasticRadius[i]>0)
- {
- G4Tubs* solidLead = new G4Tubs( Name+"_Lead" ,
- 0 ,
- m_PlasticRadius[i] ,
- m_LeadThickness[i]/2 ,
- 0*deg ,
- 360*deg );
-
- G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0);
- G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ;
- logicLead->SetVisAttributes(LeadVisAtt) ;
-
- PVPBuffer = new G4PVPlacement( 0 ,
- Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() ,
- logicLead ,
- Name+"_Lead" ,
- world ,
- false ,
- 0 );
-
- }
}
// Add Detector branch to the EventTree.
diff --git a/NPSimulation/src/PrimaryGeneratorAction.cc b/NPSimulation/src/PrimaryGeneratorAction.cc
index a5d8a264e2eb54c33cb78404ac210f7ee0ab6ba4..c6a852cdb4fe334a5a68e3425e3bd3908571b835 100644
--- a/NPSimulation/src/PrimaryGeneratorAction.cc
+++ b/NPSimulation/src/PrimaryGeneratorAction.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -35,6 +35,7 @@
#include "EventGeneratorTransfertToResonance.hh"
#include "EventGeneratorIsotropic.hh"
#include "EventGeneratorBeam.hh"
+#include "EventGeneratorPhaseSpace.hh"
#include <cstdlib>
@@ -83,10 +84,11 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
ifstream EventGeneratorFile;
EventGeneratorFile.open(Path.c_str());
- bool check_TransfertToResonance = false ;
- bool check_Isotropic = false ;
- bool check_Transfert = false ;
- bool check_Beam = false ;
+ bool check_TransfertToResonance = false ;
+ bool check_PhaseSpace = false ;
+ bool check_Isotropic = false ;
+ bool check_Transfert = false ;
+ bool check_Beam = false ;
if (EventGeneratorFile.is_open()) cout << " Event Generator file " << Path << " loading " << endl ;
else {
@@ -151,6 +153,18 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
myEventGenerator->SetTarget(m_detector->GetTarget());
m_EventGenerator = myEventGenerator ;
}
+
+ //Search for Transfert To Resonance
+ else if (LineBuffer.compare(0, 10, "PhaseSpace") == 0 && !check_PhaseSpace) {
+ check_PhaseSpace = true ;
+ VEventGenerator* myEventGenerator = new EventGeneratorPhaseSpace() ;
+ EventGeneratorFile.close() ;
+ myEventGenerator->ReadConfiguration(Path) ;
+ EventGeneratorFile.open(Path.c_str()) ;
+ myEventGenerator->InitializeRootOutput() ;
+ myEventGenerator->SetTarget(m_detector->GetTarget()) ;
+ m_EventGenerator = myEventGenerator ;
+ }
}
EventGeneratorFile.close();
}
diff --git a/NPSimulation/src/SteppingVerbose.cc b/NPSimulation/src/SteppingVerbose.cc
index 73018998e694a372f53a664bea20aa5e4680b77b..eb8acb7d353702b254f71db0961605f3360fe23b 100644
--- a/NPSimulation/src/SteppingVerbose.cc
+++ b/NPSimulation/src/SteppingVerbose.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index b8d5ef2c33790050ca93a519616e65bae609935b..6e8d1d6b4b46bb7ab898c5ddf8309e0c8ab16f80 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -77,6 +77,7 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
if (Pressure == 1) {
G4cout << "CryoTarget pressure set to 1 bar" << G4endl;
+
if (Temperature == 24) {
density = 0.0020182 * g / cm3;
G4cout << "CryoTarget temp set to 24K" << G4endl;
@@ -91,14 +92,17 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
density = 0.001863 * g / cm3;
G4cout << "CryoTarget temp set to 26K" << G4endl;
}
+
+ else if (Temperature == 30) {
+ density = 0.00020475 * g / cm3;
+ G4cout << "CryoTarget temp set to 30K" << G4endl;
+ }
else {
G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
}
}
-
-
else if (Pressure == 0.5) {
G4cout << "CryoTarget pressure set to 0.5 bar" << G4endl;
@@ -117,6 +121,17 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
G4cout << "CryoTarget temp set to 26K" << G4endl;
}
+
+ else if (Pressure == 0.7) {
+ G4cout << "CryoTarget pressure set to 0.7 bar" << G4endl;
+
+ if (Temperature == 26) {
+ density = 0.0013125 * g / cm3;
+ G4cout << "CryoTarget temp set to 26K" << G4endl;
+ }
+ }
+
+
else {
G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
}
@@ -294,21 +309,21 @@ void Target::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 2, "X=") == 0) {
check_X = true ;
ConfigFile >> DataBuffer;
- m_TargetX = atoi(DataBuffer.c_str()) * mm;
+ m_TargetX = atof(DataBuffer.c_str()) * mm;
cout << "Target coordinate (mm): ( " << m_TargetX / mm << " ; ";
}
else if (DataBuffer.compare(0, 2, "Y=") == 0) {
check_Y = true ;
ConfigFile >> DataBuffer;
- m_TargetY = atoi(DataBuffer.c_str()) * mm;
+ m_TargetY = atof(DataBuffer.c_str()) * mm;
cout << m_TargetY / mm << " ; ";
}
else if (DataBuffer.compare(0, 2, "Z=") == 0) {
check_Z = true ;
ConfigFile >> DataBuffer;
- m_TargetZ = atoi(DataBuffer.c_str()) * mm;
+ m_TargetZ = atof(DataBuffer.c_str()) * mm;
cout << m_TargetZ / mm << " )" << endl ;
}
@@ -355,13 +370,13 @@ void Target::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 12, "TEMPERATURE=") == 0) {
check_Temperature = true ;
ConfigFile >> DataBuffer;
- m_TargetTemperature = atoi(DataBuffer.c_str());
+ m_TargetTemperature = atof(DataBuffer.c_str());
}
else if (DataBuffer.compare(0, 9, "PRESSURE=") == 0) {
check_Pressure = true ;
ConfigFile >> DataBuffer;
- m_TargetPressure = atoi(DataBuffer.c_str());
+ m_TargetPressure = atof(DataBuffer.c_str());
}
else if (DataBuffer.compare(0, 9, "MATERIAL=") == 0) {
@@ -374,7 +389,7 @@ void Target::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 17, "WINDOWSTHICKNESS=") == 0) {
check_WinThickness = true ;
ConfigFile >> DataBuffer;
- m_WindowsThickness = atoi(DataBuffer.c_str()) * micrometer;
+ m_WindowsThickness = atof(DataBuffer.c_str()) * micrometer;
cout << "Windows Thickness: " << m_WindowsThickness / micrometer << "um" << endl ;
}
@@ -388,21 +403,21 @@ void Target::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 2, "X=") == 0) {
check_X = true ;
ConfigFile >> DataBuffer;
- m_TargetX = atoi(DataBuffer.c_str()) * mm;
+ m_TargetX = atof(DataBuffer.c_str()) * mm;
cout << "Target coordinate (mm): ( " << m_TargetX / mm << " ; ";
}
else if (DataBuffer.compare(0, 2, "Y=") == 0) {
check_Y = true ;
ConfigFile >> DataBuffer;
- m_TargetY = atoi(DataBuffer.c_str()) * mm;
+ m_TargetY = atof(DataBuffer.c_str()) * mm;
cout << m_TargetY / mm << " ; ";
}
else if (DataBuffer.compare(0, 2, "Z=") == 0) {
check_Z = true ;
ConfigFile >> DataBuffer;
- m_TargetZ = atoi(DataBuffer.c_str()) * mm;
+ m_TargetZ = atof(DataBuffer.c_str()) * mm;
cout << m_TargetZ / mm << " )" << endl ;
}
@@ -561,9 +576,8 @@ void Target::CalculateBeamInteraction( double MeanPosX, double SigmaPosX, double
if (AngleIncidentPhi < 0) AngleIncidentPhi += 2*pi ;
if (AngleIncidentTheta < 1e-6) AngleIncidentPhi = 0 ;
- // Calculation of effective target thickness and z-position of interaction
+ // effective target thickness is dz and z-position of interaction
// when the target is tilted wrt the beam axis
- double EffectiveThickness = m_TargetThickness / (BeamDir.unit()).dot(TargetNormal.unit());
double uniform = RandFlat::shoot();
z0 = dz + (-m_TargetThickness / 2 + uniform * m_TargetThickness);
diff --git a/NPSimulation/src/ThinSi.cc b/NPSimulation/src/ThinSi.cc
index 3af43f06ac2bb0ff45dc5cae53c8bd7b388273dd..9c12bc02eb1ff80120b7ccab2eef2f2f7ed79219 100644
--- a/NPSimulation/src/ThinSi.cc
+++ b/NPSimulation/src/ThinSi.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/ThinSiScorers.cc b/NPSimulation/src/ThinSiScorers.cc
index 53faa61b07698e34d45c0fd94c3e67aab908c55c..6ea27810b730fc892bf4bf3b4825e3bebba72496 100644
--- a/NPSimulation/src/ThinSiScorers.cc
+++ b/NPSimulation/src/ThinSiScorers.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/VDetector.cc b/NPSimulation/src/VDetector.cc
index df74f1f014ae9d3907e9fdc8783898f5add0516b..59ce7bc05e694cb3f7e97c62ddae95adb54e1719 100644
--- a/NPSimulation/src/VDetector.cc
+++ b/NPSimulation/src/VDetector.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index bf2b8e82b1f30a6be7f2e75b8f067898ca3d2347..2da902dfea0faf122ce31ff2ba6818b53dd781e1 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
diff --git a/README b/README
index a0cad80d2ede8ac15ab6d25d35f5e14b8aef0e64..63df8bcb12fa09deea23089a4c530d42b083922e 100644
--- a/README
+++ b/README
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *