diff --git a/Inputs/EventGenerator/60Fe.reaction b/Inputs/EventGenerator/60Fe.reaction
index 6430739f1ca663828d0ae1c5f912642e5f2a9aa0..467bd51d9509590ead464e0f56adc36d461b14b9 100644
--- a/Inputs/EventGenerator/60Fe.reaction
+++ b/Inputs/EventGenerator/60Fe.reaction
@@ -12,8 +12,8 @@ Transfert
BeamEnergySpread= 0
SigmaX= 1
SigmaY= 1
- SigmaThetaX= 0.5
- SigmaPhiY= 0.5
+ SigmaThetaX= 2
+ SigmaPhiY= 2
CrossSectionPath= ni69_g7_01.n
ShootLight= 1
ShootHeavy= 0
diff --git a/NPAnalysis/Gaspard/Result/myResult.root b/NPAnalysis/Gaspard/Result/myResult.root
index 040f6ae05a7b10050c1f4674773a35c8f4b62b8a..eb7aa8c9045793132a6c1c12c390e3b06e31eeb7 100644
Binary files a/NPAnalysis/Gaspard/Result/myResult.root and b/NPAnalysis/Gaspard/Result/myResult.root differ
diff --git a/NPAnalysis/Gaspard/macros/ControlSimu.C b/NPAnalysis/Gaspard/macros/ControlSimu.C
new file mode 100644
index 0000000000000000000000000000000000000000..56f2ee46029f7451fcf69349a615f5cd3b5a17bc
--- /dev/null
+++ b/NPAnalysis/Gaspard/macros/ControlSimu.C
@@ -0,0 +1,112 @@
+#include <iostream>
+
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TSystem.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TTree.h"
+#include "TBranch.h"
+#include "TH1F.h"
+#include "TH2F.h"
+
+#include "TGaspardTrackerData.h"
+#include "TInitialConditions.h"
+#include "TInteractionCoordinates.h"
+
+void ControlSimu(const char * fname = "mySimul")
+{
+ // Open output ROOT file from NPTool simulation run
+ TString path = gSystem->Getenv("NPTOOL");
+ path += "/Outputs/Simulation/";
+ TString inFileName = fname;
+ if (!inFileName.Contains("root")) inFileName += ".root";
+ TFile *inFile = new TFile(path + inFileName);
+ TTree *tree = (TTree*) inFile->Get("SimulatedTree");
+
+ // Connect the branches of the TTree and activate then if necessary
+ // TInitialConditions branch
+ TInitialConditions *initCond = 0;
+ tree->SetBranchAddress("InitialConditions", &initCond);
+ tree->SetBranchStatus("InitialConditions", 1);
+ // TInteractionCoordinates branch
+ TInteractionCoordinates *interCoord = 0;
+ tree->SetBranchAddress("InteractionCoordinates", &interCoord);
+ tree->SetBranchStatus("InteractionCoordinates", 0);
+
+ // Declare histograms
+ // for incident beam
+ TH2F *hEmittanceXY = new TH2F("hEmittanceXY", "Position on Target (x, y)", 200, -10, 10, 200, -10, 10);
+ TH2F *hEmittanceXTheta = new TH2F("hEmittanceXTheta", "Emittance (x, #theta)", 200, -10, 10, 200, -10, 10);
+ TH2F *hEmittanceYPhi = new TH2F("hEmittanceYPhi", "Emittance (y, #phi)", 200, -10, 10, 200, -10, 10);
+ TH1F *hIncidentTheta = new TH1F("hIncidentTheta", "Incident Theta", 100, 0, 10);
+ TH1F *hIncidentPhi = new TH1F("hIncidentPhi", "Incident Phi", 360, 0, 360);
+ // for emitted particle
+ TH1F *hEmittedThetaCM = new TH1F("hEmittedThetaCM", "Emitted Theta CM", 180, 0, 180);
+ TH1F *hEmittedThetaIF = new TH1F("hEmittedThetaIF", "Emitted Theta (reaction frame)", 180, 0, 180);
+ TH1F *hEmittedPhiIF = new TH1F("hEmittedPhiIF", "Emitted Phi (reaction frame)", 360, 0, 360);
+ TH1F *hEmittedThetaWF = new TH1F("hEmittedThetaWF", "Emitted Theta (detector frame)", 180, 0, 180);
+ TH1F *hEmittedPhiWF = new TH1F("hEmittedPhiWF", "Emitted Phi (detector frame)", 360, 0, 360);
+ TH1F *hControlPhi = new TH1F("hControlPhi", "IncidentPhi + EmittedPhi", 360, 0, 360);
+ TH2F *hEmittedETheta = new TH2F("hEmittedETheta", "Kinematic line", 180, 0, 180, 200, 0, 150);
+
+ // Read the TTree
+ Int_t nentries = tree->GetEntries();
+ cout << "TTree contains " << nentries << " events" << endl;
+ for (Int_t i = 0; i < nentries; i++) {
+ if (i%10000 == 0) cout << "Entry " << i << endl;
+ tree->GetEntry(i);
+
+ // Fill histos
+ // incident beam
+ hEmittanceXY -> Fill(initCond->GetICPositionX(0), initCond->GetICPositionY(0));
+ hEmittanceXTheta -> Fill(initCond->GetICPositionX(0), initCond->GetICIncidentEmittanceTheta(0));
+ hEmittanceYPhi -> Fill(initCond->GetICPositionX(0), initCond->GetICIncidentEmittancePhi(0));
+ hIncidentTheta -> Fill(initCond->GetICIncidentAngleTheta(0));
+ hIncidentPhi -> Fill(initCond->GetICIncidentAnglePhi(0));
+ // ejected particle
+ hEmittedThetaCM -> Fill(initCond->GetICEmittedAngleThetaCM(0));
+ hEmittedThetaIF -> Fill(initCond->GetICEmittedAngleThetaLabIncidentFrame(0));
+ hEmittedPhiIF -> Fill(initCond->GetICEmittedAnglePhiIncidentFrame(0));
+ hEmittedThetaWF -> Fill(initCond->GetICEmittedAngleThetaLabWorldFrame(0));
+ hEmittedPhiWF -> Fill(initCond->GetICEmittedAnglePhiWorldFrame(0));
+ hEmittedETheta -> Fill(initCond->GetICEmittedAngleThetaLabIncidentFrame(0), initCond->GetICEmittedEnergy(0));
+
+ // Control histo
+ Double_t phi_control = initCond->GetICIncidentAnglePhi(0) + initCond->GetICEmittedAnglePhiIncidentFrame(0);
+ if (phi_control > 360) phi_control -= 360;
+ hControlPhi->Fill(phi_control);
+ }
+
+ // Display histograms
+ TCanvas *canvas1 = new TCanvas("canvas1", "Incident beam properties");
+ canvas1->Divide(2,3);
+ canvas1->cd(1);
+ hEmittanceXTheta->Draw("colz");
+ canvas1->cd(2);
+ hEmittanceYPhi->Draw("colz");
+ canvas1->cd(3);
+ hIncidentTheta->Draw();
+ canvas1->cd(4);
+ hIncidentPhi->Draw();
+ canvas1->cd(5);
+ hEmittanceXY->Draw("colz");
+
+ // Display histograms
+ TCanvas *canvas2 = new TCanvas("canvas2", "Emitted particle properties");
+ canvas2->Divide(2,3);
+ canvas2->cd(1);
+ hEmittedThetaCM->Draw();
+ canvas2->cd(2);
+ hEmittedETheta->Draw();
+ canvas2->cd(3);
+ hEmittedThetaIF->Draw();
+ canvas2->cd(4);
+ hEmittedPhiIF->Draw();
+ canvas2->cd(5);
+ hEmittedThetaWF->Draw();
+ canvas2->cd(6);
+ hEmittedPhiWF->Draw();
+// hControlPhi->SetLineColor(kRed);
+// hControlPhi->Draw("same");
+}
diff --git a/NPAnalysis/Gaspard/macros/ControlSimuGaspard.C b/NPAnalysis/Gaspard/macros/ControlSimuGaspard.C
deleted file mode 100644
index 9bb84c80f49e0fbdb3e8ea35c5361961e4f5cd57..0000000000000000000000000000000000000000
--- a/NPAnalysis/Gaspard/macros/ControlSimuGaspard.C
+++ /dev/null
@@ -1,84 +0,0 @@
-#include <iostream>
-
-#include "TROOT.h"
-#include "TSystem.h"
-#include "TFile.h"
-#include "TString.h"
-#include "TTree.h"
-#include "TBranch.h"
-#include "TH1F.h"
-
-#include "TGaspardTrackerData.h"
-#include "TInitialConditions.h"
-#include "TInteractionCoordinates.h"
-
-void ControlSimuGaspard(const char * fname = "mySimul")
-{
- // Open output ROOT file from NPTool simulation run
- TString path = gSystem->Getenv("NPTOOL");
- path += "/Outputs/";
- TString inFileName = fname;
- if (!inFileName.Contains("root")) inFileName += ".root";
- TFile *inFile = new TFile(path + inFileName);
- TTree *tree = (TTree*) inFile->Get("EventTree");
-
- // Connect the branches of the TTree and activate then if necessary
- // TGaspardTrackerData branch
- TGaspardTrackerData *gpdTrkData = 0;
- tree->SetBranchAddress("GASPARD", &gpdTrkData);
- tree->SetBranchStatus("GASPARD", 0);
- // TInitialConditions branch
- TInitialConditions *initCond = 0;
- tree->SetBranchAddress("InitialConditions", &initCond);
- tree->SetBranchStatus("InitialConditions", 1);
- // TInteractionCoordinates branch
- TInteractionCoordinates *interCoord = 0;
- tree->SetBranchAddress("InteractionCoordinates", &interCoord);
- tree->SetBranchStatus("InteractionCoordinates", 0);
-
- // Prepare histograms
- TH1F *hIncidentTheta = new TH1F("hIncidentTheta", "Incident Theta", 180, 0, 180);
- TH1F *hIncidentPhi = new TH1F("hIncidentPhi", "Incident Phi", 360, 0, 360);
- TH1F *hEmittedThetaIF = new TH1F("hEmittedThetaIF", "Emitted Theta (incident frame)", 180, 0, 180);
- TH1F *hEmittedPhiIF = new TH1F("hEmittedPhiIF", "Emitted Phi (incident frame)", 360, 0, 360);
- TH1F *hEmittedThetaWF = new TH1F("hEmittedThetaWF", "Emitted Theta (laboratory frame)", 180, 0, 180);
- TH1F *hEmittedPhiWF = new TH1F("hEmittedPhiWF", "Emitted Phi (laboratory frame)", 360, 0, 360);
- TH1F *hControlPhi = new TH1F("hControlPhi", "IncidentPhi + EmittedPhi", 360, 0, 360);
-
- // Read the TTree
- Int_t nentries = tree->GetEntries();
- cout << "TTree contains " << nentries << " events" << endl;
- for (Int_t i = 0; i < nentries; i++) {
- if (i%1000 == 0) cout << "Entry " << i << endl;
- tree->GetEntry(i);
- // Fill histos
- hIncidentTheta->Fill(initCond->GetICIncidentAngleTheta(0));
- hIncidentPhi->Fill(initCond->GetICIncidentAnglePhi(0));
- hEmittedThetaIF->Fill(initCond->GetICEmittedAngleThetaLabIncidentFrame(0));
- hEmittedPhiIF->Fill(initCond->GetICEmittedAnglePhiIncidentFrame(0));
- hEmittedThetaWF->Fill(initCond->GetICEmittedAngleThetaLabWorldFrame(0));
- hEmittedPhiWF->Fill(initCond->GetICEmittedAnglePhiWorldFrame(0));
- // Control histo
- Double_t phi_control = initCond->GetICIncidentAnglePhi(0) + initCond->GetICEmittedAnglePhiIncidentFrame(0);
- if (phi_control > 360) phi_control -= 360;
- hControlPhi->Fill(phi_control);
- }
-
- // Display histograms
- TCanvas *canvas = new TCanvas("canvas", "Angles");
- canvas->Divide(2,3);
- canvas->cd(1);
- hIncidentTheta->Draw();
- canvas->cd(2);
- hIncidentPhi->Draw();
- canvas->cd(3);
- hEmittedThetaIF->Draw();
- canvas->cd(4);
- hEmittedPhiIF->Draw();
- canvas->cd(5);
- hEmittedThetaWF->Draw();
- canvas->cd(6);
- hEmittedPhiWF->Draw();
- hControlPhi->SetLineColor(kRed);
- hControlPhi->Draw("same");
-}
diff --git a/NPAnalysis/Gaspard/macros/GeometricalEfficiency.C b/NPAnalysis/Gaspard/macros/GeometricalEfficiency.C
index ae87e6e7a49bc5ec59a8408033542d2d4c5f4f8f..90ca0e83b863b7b45116682d97a9fafb246068af 100644
--- a/NPAnalysis/Gaspard/macros/GeometricalEfficiency.C
+++ b/NPAnalysis/Gaspard/macros/GeometricalEfficiency.C
@@ -15,11 +15,11 @@ void GeometricalEfficiency(const char * fname = "Efficiency_10000")
{
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
- path += "/Outputs/";
+ path += "/Outputs/Simulation/";
TString inFileName = fname;
if (!inFileName.Contains("root")) inFileName += ".root";
TFile *inFile = new TFile(path + inFileName);
- TTree *tree = (TTree*) inFile->Get("EventTree");
+ TTree *tree = (TTree*) inFile->Get("SimulatedTree");
// Connect the branches of the TTree and activate then if necessary
// TInitialConditions branch
diff --git a/NPAnalysis/Gaspard/macros/TestAngleReconstruction.C b/NPAnalysis/Gaspard/macros/TestAngleReconstruction.C
index 7f193f0ff8b20865379a78b128e6095a8d08f5e8..217efe959dae85053ab602d6efdc92ff9b05df2b 100644
--- a/NPAnalysis/Gaspard/macros/TestAngleReconstruction.C
+++ b/NPAnalysis/Gaspard/macros/TestAngleReconstruction.C
@@ -19,7 +19,7 @@ void TestAngleReconstruction(const char * fname = "mySimul")
{
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
- path += "/Outputs/";
+ path += "/Outputs/Simulation/";
TString inFileName = fname;
if (!inFileName.Contains("root")) inFileName += ".root";
TFile *inFile = new TFile(path + inFileName);
diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc
index a1ad0575239021cad3339f443675fc41da450421..ad75bb183b1ccf0651bc89fa15d53c912c222f65 100644
--- a/NPAnalysis/Gaspard/src/Analysis.cc
+++ b/NPAnalysis/Gaspard/src/Analysis.cc
@@ -117,7 +117,7 @@ int main(int argc,char** argv)
Int_t detecXT = EventGPD->GetGPDTrkFirstStageFrontTDetectorNbr(0) / det_ref;
Int_t detecYE = EventGPD->GetGPDTrkFirstStageBackEDetectorNbr(0) / det_ref;
Int_t detecYT = EventGPD->GetGPDTrkFirstStageBackTDetectorNbr(0) / det_ref;
-// det_ref -= 1000; // for TGaspardTrackerDummyShape
+ det_ref -= 1000; // for TGaspardTrackerDummyShape
// case of same detector
if (detecXE*detecXT*detecYE*detecYT == 1) {
// calculate strip number
diff --git a/NPLib/GASPARD/GaspardTracker.cxx b/NPLib/GASPARD/GaspardTracker.cxx
index 86af36f63a08a8dea5102732dfedc14d3cd018f1..d9971e51b4cbcb4c8a5b24d90a3be91dc2b0131a 100644
--- a/NPLib/GASPARD/GaspardTracker.cxx
+++ b/NPLib/GASPARD/GaspardTracker.cxx
@@ -279,7 +279,7 @@ void GaspardTracker::ReadConfiguration(string Path)
ConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
}
// Finding another telescope (safety), toggle out
- else if (DataBuffer.compare(0, 9, "GPDSquare") == 0) {
+ else if (DataBuffer.compare(0, 13, "GPDDummyShape") == 0) {
cout << "WARNING: Another Module is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
ReadingStatus = false;
}
diff --git a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
index 5abd42f89bab02139d82512f8094d5fac1e2a478..6bca92aea82e27c4d075f453a5f9370a91a82093 100644
--- a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
+++ b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
@@ -79,6 +79,7 @@ void TGaspardTrackerPhysics::BuildPhysicalEvent(TGaspardTrackerData* Data)
int detecXT = Data->GetGPDTrkFirstStageFrontTDetectorNbr(0) / det_ref;
int detecYE = Data->GetGPDTrkFirstStageBackEDetectorNbr(0) / det_ref;
int detecYT = Data->GetGPDTrkFirstStageBackTDetectorNbr(0) / det_ref;
+ det_ref -= 1000; // only for GaspardDummyShape
// case of same detector
if (detecXE*detecXT*detecYE*detecYT == 1) {
diff --git a/NPSimulation/src/EventGeneratorIsotropic.cc b/NPSimulation/src/EventGeneratorIsotropic.cc
index 9f1faf0ff77e0e0ca2f01b195bbe65994dfa20a5..c60482c0471de284466993cd5090327195323136 100644
--- a/NPSimulation/src/EventGeneratorIsotropic.cc
+++ b/NPSimulation/src/EventGeneratorIsotropic.cc
@@ -210,6 +210,13 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent, G4ParticleGun* par
m_InitConditions->SetICPositionX(m_x0 / mm);
m_InitConditions->SetICPositionY(m_y0 / mm);
m_InitConditions->SetICPositionZ(m_z0 / mm);
+ // Incident "particles"
+ // Everything is zero for a source
+ m_InitConditions->SetICIncidentEmittanceTheta(0);
+ m_InitConditions->SetICIncidentEmittancePhi(0);
+ m_InitConditions->SetICIncidentAngleTheta(0);
+ m_InitConditions->SetICIncidentAnglePhi(0);
+ m_InitConditions->SetICIncidentEnergy(0);
// Emitted particle angles
m_InitConditions->SetICEmittedAngleThetaCM(theta / deg);
m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(theta / deg);
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 28521e20ea51b783988d3ad91cb42d0cbcb982d5..12a9205e554c38723eab62ef7e6bc12d8c834da3 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -354,12 +354,13 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
// Calculate Angle in spherical coordinate, passing by the direction vector dir
- G4double Xdir = cos( pi/2. - Beam_thetaX );
- G4double Ydir = cos( pi/2. - Beam_phiY );
- G4double Zdir = sin( pi/2. - Beam_thetaX ) + sin( pi/2. - Beam_phiY);
+ G4double Xdir = cos(pi/2. - Beam_thetaX);
+ G4double Ydir = cos(pi/2. - Beam_phiY );
+ G4double Zdir = sin(pi/2. - Beam_thetaX) + sin(pi/2. - Beam_phiY);
- G4double Beam_theta = acos ( Zdir / sqrt( Xdir*Xdir + Ydir*Ydir + Zdir*Zdir ) );
- G4double Beam_phi = atan2( Ydir , Xdir ) ;
+ G4double Beam_theta = acos(Zdir / sqrt(Xdir*Xdir + Ydir*Ydir + Zdir*Zdir)) * rad;
+ G4double Beam_phi = atan2(Ydir, Xdir) * rad;
+ if (Beam_phi < 0) Beam_phi += 2*pi;
// write angles to ROOT file
m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);