diff --git a/NPLib/Detectors/Epic/TEpicData.h b/NPLib/Detectors/Epic/TEpicData.h
index 13d0ad82dd25056c8778047fde251de5e6b2b227..6e63f6a435e58647287cddd594254bae36897403 100644
--- a/NPLib/Detectors/Epic/TEpicData.h
+++ b/NPLib/Detectors/Epic/TEpicData.h
@@ -42,8 +42,10 @@ class TEpicData : public TObject {
Double_t Time_HF;
Double_t ToF;
Bool_t isFakeFission;
+ // vector of {Zstep, DEstep, DTstep}
std::vector<Double_t> PosZ;
std::vector<Double_t> DE;
+ std::vector<Double_t> DT; // DT=Ti-T0, with T0=Tfission or Talpha
// Getters
UShort_t GetAnodeNbr() const { return AnodeNbr; }
@@ -56,6 +58,7 @@ class TEpicData : public TObject {
Bool_t IsFakeFission() const { return isFakeFission; }
const std::vector<Double_t>& GetPosZ() const { return PosZ; }
const std::vector<Double_t>& GetDE() const { return DE; }
+ const std::vector<Double_t>& GetDT() const { return DT; }
};
//////////////////////////////////////////////////////////////
@@ -86,10 +89,10 @@ class TEpicData : public TObject {
// add //! to avoid ROOT creating dictionnary for the methods
public:
////////////////////// SETTERS ////////////////////////
- void Set(const UShort_t& n, const Double_t& q1, const Double_t& q2, const Double_t& qmax,
+ void Set(const UShort_t& n, const Double_t& q1, const Double_t& q2, const Double_t& qmax,
const Double_t& t, const Double_t& thf, const Double_t& tof, const Bool_t& fake,
- const std::vector<double>& z, const std::vector<double>& de) {
- fEpic_Data.push_back({n, q1, q2, qmax, t, thf, tof, fake, z, de});
+ const std::vector<double>& z, const std::vector<double>& de, const std::vector<double>& dt) {
+ fEpic_Data.push_back({n, q1, q2, qmax, t, thf, tof, fake, z, de, dt});
}
const EpicAnodeData& operator[](const unsigned int& i) const {return fEpic_Data[i];}
@@ -105,6 +108,7 @@ class TEpicData : public TObject {
Bool_t GetFakeFissionStatus(const unsigned int& i) const { return fEpic_Data[i].isFakeFission; }
const std::vector<Double_t>& GetPosZ(const unsigned int& i) const { return fEpic_Data[i].PosZ; }
const std::vector<Double_t>& GetDE(const unsigned int& i) const { return fEpic_Data[i].DE; }
+ const std::vector<Double_t>& GetDT(const unsigned int& i) const { return fEpic_Data[i].DT; }
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
ClassDef(TEpicData,1) // EpicData structure
diff --git a/NPLib/Detectors/Epic/TEpicPhysics.cxx b/NPLib/Detectors/Epic/TEpicPhysics.cxx
index 0085c780c9ba0739759b71f866e0e9da82d89cc2..983188d38a25924e4bcbee30d7c450d09daee411 100644
--- a/NPLib/Detectors/Epic/TEpicPhysics.cxx
+++ b/NPLib/Detectors/Epic/TEpicPhysics.cxx
@@ -113,23 +113,24 @@ void TEpicPhysics::PreTreat() {
unsigned int mysize = m_EventData->GetMultiplicity();
for (UShort_t i = 0; i < mysize ; ++i) {
Double_t Q1 = m_EventData->GetQ1(i);
- Double_t Q2 = m_EventData->GetQ2(i);
- Double_t Qmax = m_EventData->GetQmax(i);
if (Q1 > m_E_Threshold) {
-
unsigned short AnodeNumber = m_EventData->GetAnodeNbr(i);
double TimeOffset = 0;
TimeOffset = Cal->GetValue("Epic/ANODE"+NPL::itoa(AnodeNumber)+"_TIMEOFFSET",0);
double Time = m_EventData->GetTime(i);
double TimeHF = m_EventData->GetTimeHF(i);
- double tof = Time - m_EventData->GetTimeHF(i) - TimeOffset;
- bool kFakeStatus = m_EventData->GetFakeFissionStatus(i);
- vector<double> posz = m_EventData->GetPosZ(i);
- vector<double> de = m_EventData->GetDE(i);
+ double tof = Time - TimeHF - TimeOffset;
if(tof < 0){
tof += Cal->GetValue("Epic/PULSE_TIMEOFFSET",0) ;
}
- m_PreTreatedData->Set(AnodeNumber, Q1, Q2, Qmax, Time, TimeHF, tof, kFakeStatus, posz, de);
+ m_PreTreatedData->Set(AnodeNumber, Q1,
+ m_EventData->GetQ2(i),
+ m_EventData->GetQmax(i),
+ Time, TimeHF, tof,
+ m_EventData->GetFakeFissionStatus(i),
+ m_EventData->GetPosZ(i),
+ m_EventData->GetDE(i),
+ m_EventData->GetDT(i));
}
}
diff --git a/NPSimulation/Detectors/Coaxial_Germanium/Coaxial_Germanium.cc b/NPSimulation/Detectors/Coaxial_Germanium/Coaxial_Germanium.cc
index c3c46abed99ff36edefdc1622f32cc4525f46829..0628edaa5935aa290c228fae2f19756229a58051 100644
--- a/NPSimulation/Detectors/Coaxial_Germanium/Coaxial_Germanium.cc
+++ b/NPSimulation/Detectors/Coaxial_Germanium/Coaxial_Germanium.cc
@@ -120,9 +120,9 @@ G4LogicalVolume* Coaxial_Germanium::BuildDetector(G4int DetNumber, G4ThreeVector
vol_Coaxial_Germanium->SetVisAttributes(light_GreyAtt);
// Germanium crystal
- G4Tubs* Coaxial_Germanium_crys = new G4Tubs("Coaxial_Germanium_crys",0 , 13*mm, 10.2*0.5*mm, 0, 360*deg);
+ G4Tubs* Coaxial_Germanium_crys = new G4Tubs("Coaxial_Germanium_crys",0 , 25.*mm, 35.*mm, 0, 360*deg);
G4LogicalVolume* vol_crys = new G4LogicalVolume(Coaxial_Germanium_crys, m_MaterialGermanium, "logic_Coaxial_Germanium_crys", 0, 0, 0);
- G4ThreeVector crys_Pos = G4ThreeVector(0, 0, 0.6*1.5*mm+11*0.5*mm);
+ G4ThreeVector crys_Pos = G4ThreeVector(0, 0, 50.*mm);
new G4PVPlacement(0, crys_Pos, vol_crys, "Coaxial_Germanium_crys", logicCoaxial_Germanium, false, DetNumber);
vol_crys->SetVisAttributes(RedAtt);
diff --git a/NPSimulation/Detectors/Epic/Epic.cc b/NPSimulation/Detectors/Epic/Epic.cc
index 04d03022982e2417a1c2dbb59a864b8296f6903c..bd3c520651e38c99a53415501c873be553517bdd 100644
--- a/NPSimulation/Detectors/Epic/Epic.cc
+++ b/NPSimulation/Detectors/Epic/Epic.cc
@@ -523,7 +523,7 @@ void Epic::InitializeRootOutput(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAvtion
+// Called at in the EventAction::EndOfEventAction
void Epic::ReadSensitive(const G4Event* ){
m_Event->Clear();
@@ -531,18 +531,15 @@ void Epic::ReadSensitive(const G4Event* ){
// GaseousDetector scorer
GaseousDetectorScorers::PS_GaseousDetector* Scorer= (GaseousDetectorScorers::PS_GaseousDetector*) m_EpicScorer->GetPrimitive(0);
- unsigned int size = Scorer->GetMult();
+ unsigned int size = Scorer->GetMult();
for(unsigned int i = 0 ; i < size ; i++){
vector<unsigned int> level = Scorer->GetLevel(i);
//double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Epic_NS::ResoEnergy);
- double Energy = Scorer->GetEnergy(i);
- vector<string> slice_name = Scorer->GetParticleName(i);
- vector<double> slice_posZ = Scorer->GetSlicePosZ(i);
- vector<double> slice_DE = Scorer->GetEnergyLossPerSlice(i);
- //cout << "level[0] = " << level[0] << endl;
- //cout << "size slice_name : " << slice_name.size() << endl;
- //cout << "size slice_posZ : " << slice_posZ.size() << endl;
- //cout << "size slice_DE : " << slice_DE.size() << endl;
+ double Energy = Scorer->GetEnergy(i); // Attention Energy of all particles (FF, alpha, e-, ...)
+ vector<string> step_name = Scorer->GetParticleName(i);
+ vector<double> step_posZ = Scorer->GetStepPosZ(i);
+ vector<double> step_DE = Scorer->GetEnergyLossPerStep(i);
+ vector<double> step_time = Scorer->GetStepTime(i);
double influence = 0 ;
if(Energy>Epic_NS::EnergyThreshold){
@@ -554,26 +551,30 @@ void Epic::ReadSensitive(const G4Event* ){
posZ_anode += m_mapping_A[Anode] * (2. * (m_Distance_AK*mm + m_Thickness_K*mm) + m_InterDistance_KK*mm + thickness_A) ;
vector<double> z;
vector<double> de;
+ vector<double> dt;
- for(int j=0; j<slice_name.size(); j++){
- if(slice_name.at(j)=="e-") {
+ for(int j=0; j<step_name.size(); j++){
+ //cout << "step #" << j << " , particle name = " << step_name.at(j) << " , z = " << step_posZ.at(j) << " , dE = " << step_DE.at(j) << " , t = " << step_time.at(j) << endl;
+ if(step_name.at(j)=="e-" || step_name.at(j)=="anti_nu_e") {
continue;
}
- else{
- influence += slice_DE.at(j) * TMath::Abs(slice_posZ.at(j) - posZ_anode);
- z.push_back(slice_posZ.at(j));
- de.push_back(slice_DE.at(j));
+ else{ // FF or alpha
+ //cout << "[FF or alpha]step #" << j << " , particle name = " << step_name.at(j) << " , z = " << step_posZ.at(j) << " , dE = " << step_DE.at(j) << " , t = " << step_time.at(j) << endl;
+ influence += step_DE.at(j) * TMath::Abs(step_posZ.at(j) - posZ_anode);
+ z.push_back(step_posZ.at(j));
+ de.push_back(step_DE.at(j));
+ dt.push_back(step_time.at(j) - Scorer->GetTime(i)); // dt = time from the first step in SensitiveVolume
}
}
m_Event->Set(m_mapping_A[Anode], // set anode number
- influence / m_Distance_AK*mm, // set Q1
+ influence / m_Distance_AK*mm, // set Q1 FF or alpha only
Energy, // set Q2
Energy, // set Qmax
Time, // set Time
0., // set TimeHF
0., // set ToF
0, // status FakeFission
- z,de);
+ z, de, dt);
}
}
}
diff --git a/NPSimulation/Detectors/Plastic_BEDO/Plastic_BEDO.cc b/NPSimulation/Detectors/Plastic_BEDO/Plastic_BEDO.cc
index 7b8356f8456f166b03f354b8fa6149eae4bbcd56..bb677b19678f1de043374c48ddbefaa9abcee266 100644
--- a/NPSimulation/Detectors/Plastic_BEDO/Plastic_BEDO.cc
+++ b/NPSimulation/Detectors/Plastic_BEDO/Plastic_BEDO.cc
@@ -26,6 +26,7 @@
//G4 Geometry object
#include "G4Tubs.hh"
#include "G4Box.hh"
+#include "G4SubtractionSolid.hh"
//G4 sensitive
#include "G4SDManager.hh"
@@ -103,7 +104,13 @@ void Plastic_BEDO::AddDetector(double R, double Theta, double Phi){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Plastic_BEDO::BuildCylindricalDetector(){
if(!m_CylindricalDetector){
- G4Tubs* tub = new G4Tubs("Plastic_BEDO",Plastic_BEDO_NS::Radius-Plastic_BEDO_NS::Thickness,Plastic_BEDO_NS::Radius,Plastic_BEDO_NS::Length*0.5,0,360*deg);
+
+ G4ThreeVector zTrans(0., 0., 3.*mm);
+ G4Tubs* tubOut = new G4Tubs("Plastic_BEDO",0,Plastic_BEDO_NS::Radius,Plastic_BEDO_NS::Length*0.5,0,360*deg);
+
+ G4Tubs* tubIn = new G4Tubs("Plastic_BEDO",0,Plastic_BEDO_NS::Radius-Plastic_BEDO_NS::Thickness,Plastic_BEDO_NS::Length*0.5,0,360*deg);
+
+ G4SubtractionSolid* tub = new G4SubtractionSolid("tub", tubOut, tubIn, 0, zTrans);
G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(Plastic_BEDO_NS::Material);
m_CylindricalDetector = new G4LogicalVolume(tub,DetectorMaterial,"logic_Plastic_BEDO_tub",0,0,0);
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
index 68efeba6eb0250b246e035b466b9d0900f8eaa6d..a37fd82b8102064a120ca6f322dd2669afa555bd 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
@@ -62,7 +62,6 @@ EventGeneratorGEFReader::~EventGeneratorGEFReader(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorGEFReader::SourceParameters::SourceParameters(){
-
m_x0 = 0 ;
m_y0 = 0 ;
m_z0 = 0 ;
@@ -71,6 +70,7 @@ EventGeneratorGEFReader::SourceParameters::SourceParameters(){
m_SigmaZ = 0 ;
m_Boost = 0 ;
m_direction = 'z' ;
+ m_BeamProfile = "Gauss" ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
index 1c83fad481113ff7bf58fd93519eb4f53df20518..0788d3c160ea84d6d10c2c9c9f1522f97eda9a38 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
@@ -68,6 +68,9 @@ EventGeneratorIsotropic::SourceParameters::SourceParameters(){
m_SigmaZ = 0 ;
m_particle = NULL;
m_direction = 'z' ;
+ m_SourceProfile = "Gauss" ;
+ m_DecayLaw = "off";
+ m_ActivityBq = 1. ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -117,6 +120,8 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
if(blocks[i]->HasToken("ExcitationEnergy"))
it->m_ExcitationEnergy =blocks[i]->GetVectorDouble("ExcitationEnergy","MeV");
+ if(blocks[i]->HasToken("SourceProfile"))
+ it->m_SourceProfile=blocks[i]->GetString("SourceProfile");
if(blocks[i]->HasToken("SigmaX"))
it->m_SigmaX=blocks[i]->GetDouble("SigmaX","mm");
if(blocks[i]->HasToken("SigmaY"))
@@ -125,6 +130,10 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
it->m_SigmaZ=blocks[i]->GetDouble("SigmaZ","mm");
if(blocks[i]->HasToken("Multiplicity"))
it->m_Multiplicty=blocks[i]->GetVectorInt("Multiplicity");
+ if(blocks[i]->HasToken("DecayLaw"))
+ it->m_DecayLaw=blocks[i]->GetString("DecayLaw");
+ if(blocks[i]->HasToken("Activity_Bq"))
+ it->m_ActivityBq=blocks[i]->GetDouble("ActivityBq","Bq");
}
else{
cout << "ERROR: check your input file formatting \033[0m" << endl;
@@ -153,22 +162,27 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
void EventGeneratorIsotropic::GenerateEvent(G4Event*){
for(auto& par : m_Parameters) {
+ double delta_t = 0;
for(unsigned int p=0; p<par.m_particleName.size(); p++){
for(int i=0; i<par.m_Multiplicty[p]; i++){
par.m_particle=NULL;
if(par.m_particle==NULL){
- if(par.m_particleName[p]=="gamma" || par.m_particleName[p]=="neutron" || par.m_particleName[p]=="opticalphoton" || par.m_particleName[p]=="mu+" || par.m_particleName[p]=="mu-" || par.m_particleName[p]=="e-" || par.m_particleName[p]=="e+"){
+ if(par.m_particleName[p]=="gamma" || par.m_particleName[p]=="neutron" ||
+ par.m_particleName[p]=="opticalphoton" || par.m_particleName[p]=="mu+" ||
+ par.m_particleName[p]=="mu-" || par.m_particleName[p]=="e-" || par.m_particleName[p]=="e+"){
par.m_particle = G4ParticleTable::GetParticleTable()->FindParticle(par.m_particleName[p].c_str());
}
else{
NPL::Nucleus* N = new NPL::Nucleus(par.m_particleName[p]);
par.m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),par.m_ExcitationEnergy[p]);
delete N;
+ if(par.m_DecayLaw == "on" && i>0){
+ delta_t += G4RandExponential::shoot(1./par.m_ActivityBq) ;
+ }
}
-
}
-
+ if(delta_t > 50.e-9) continue; //to do to parametrize (time_window)
G4double cos_theta_min = cos(par.m_HalfOpenAngleMin);
G4double cos_theta_max = cos(par.m_HalfOpenAngleMax);
G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
@@ -191,9 +205,18 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
G4double x0, y0, z0;
G4double momentum_x, momentum_y, momentum_z;
- x0 = RandGauss::shoot(par.m_x0,par.m_SigmaX);
- y0 = RandGauss::shoot(par.m_y0,par.m_SigmaY);
- z0 = RandGauss::shoot(par.m_z0,par.m_SigmaZ);
+ if(par.m_SourceProfile=="Flat"){
+ double rand_r = par.m_SigmaX*sqrt(RandFlat::shoot());
+ double rand_theta = RandFlat::shoot() * 2. * pi;
+ x0 = par.m_x0 + rand_r * cos(rand_theta);
+ y0 = par.m_y0 + rand_r * sin(rand_theta);
+ z0 = RandFlat::shoot(par.m_z0*mm-par.m_SigmaZ*mm, par.m_z0*mm+par.m_SigmaZ*mm);
+ }
+ else{ // by default gaussian source profile is considered
+ x0 = RandGauss::shoot(par.m_x0,par.m_SigmaX);
+ y0 = RandGauss::shoot(par.m_y0,par.m_SigmaY);
+ z0 = RandGauss::shoot(par.m_z0,par.m_SigmaZ);
+ }
if(par.m_direction == 'z')
{
// Direction of particle, energy and laboratory angle
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
index b60d01b594c0567f50c6828785c14e0b24f0a88c..e15e00f2d60b71e5faadd334951a9107691e74f3 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
@@ -63,15 +63,18 @@ private: // Source parameter from input file
G4double m_x0 ; // Vertex Position X
G4double m_y0 ; // Vertex Position Y
G4double m_z0 ; // Vertex Position Z
- G4double m_SigmaX ;
- G4double m_SigmaY ;
+ G4double m_SigmaX ; // if m_SourceProfile==Flat, m_SigmaX=radius
+ G4double m_SigmaY ; // if m_SourceProfile==Flat, m_SigmaY=0
G4double m_SigmaZ ;
TString m_direction ;
+ TString m_SourceProfile ; // either "Gauss" (by default) or "Flat" on disk with radius of m_SigmaX
G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
vector<G4double> m_ExcitationEnergy ; // Excitation energy of the emitted particle
vector<string> m_particleName ;
vector<G4int> m_Multiplicty ;
- };
+ TString m_DecayLaw ; // [on]: mult follows the decay law controlled by T1/2 [off]: fix mult
+ G4double m_ActivityBq ; // if m_DecayLaw = "on" gives here the Activity in Bq of the sample
+ };
vector<SourceParameters> m_Parameters ;
ParticleStack* m_ParticleStack ;
TF1* fEnergyDist;
diff --git a/NPSimulation/Scorers/GaseousDetectorScorers.cc b/NPSimulation/Scorers/GaseousDetectorScorers.cc
index 9de4d195b4fccf8c18490862cd879a6f421c4deb..deb5c0563908c320855e20236332f5be623a8ad1 100644
--- a/NPSimulation/Scorers/GaseousDetectorScorers.cc
+++ b/NPSimulation/Scorers/GaseousDetectorScorers.cc
@@ -73,38 +73,46 @@ G4bool PS_GaseousDetector::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
// Contain Particle Id, Energy, Time and as many copy number as nested volume
unsigned int mysize = m_NestingLevel.size();
string particlename = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
- double slice_posZ = aStep->GetPreStepPoint()->GetPosition().z();
-
+ double step_posZ = aStep->GetPreStepPoint()->GetPosition().z();
+
t_Energy = aStep->GetTotalEnergyDeposit();
- t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
+ t_Time = aStep->GetPreStepPoint()->GetGlobalTime(); // DeltaT [ns] since the begining of the simulated event up to the begining of the step
t_ParticleName.push_back(particlename);
- t_SlicePosZ.push_back(slice_posZ);
- t_EnergyLossPerSlice.push_back(t_Energy);
+ t_StepPosZ.push_back(step_posZ);
+ t_EnergyLossPerStep.push_back(t_Energy);
+ t_StepTime.push_back(t_Time);
//cout << "mysize = " << mysize << endl;
//cout << "t_ParticleName : " << particlename << endl;
+ //cout << "step_posZ : " << step_posZ << endl;
//cout << "t_Energy : " << t_Energy << endl;
- //cout << "slice_posZ : " << slice_posZ << endl; // it seems that step_size is not uniform
- //cout << "t_ParticleName.size() = " << t_ParticleName.size() << endl;
- //cout << "t_SlicePosZ.size() = " << t_SlicePosZ.size() << endl;
- //cout << "t_EnergyLossPerSlice.size() = " << t_EnergyLossPerSlice.size() << endl;
+ //cout << "step_time : " << t_Time << endl;
+ //cout << "t_ParticleName.size() = " << t_ParticleName.size() << endl;
+ //cout << "t_StepPosZ.size() = " << t_StepPosZ.size() << endl;
+ //cout << "t_EnergyLossPerStep.size() = " << t_EnergyLossPerStep.size() << endl;
+ //cout << "t_StepTime.size() = " << t_StepTime.size() << endl;
t_Level.clear();
for (unsigned int i = 0; i < mysize; i++) {
t_Level.push_back(
aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(
m_NestingLevel[i]));
}
- // Check if the particle has interact before, if yes, add up the energies.
+ // Check if the particle has interact before, if yes, add up the energies and update the track
vector<GaseousDetectorData>::iterator it;
it = m_Data.find(GaseousDetectorData::CalculateIndex(t_Level));
if (it != m_Data.end()) {
+ // update the m_Data vector at each step
it->Add(t_Energy);
it->SetParticleName(particlename);
- it->SetSlicePosZ(slice_posZ);
- it->SetEnergyLossPerSlice(t_Energy);
- } else {
- //m_Data.Set(t_Energy, t_Time, t_Level, t_ParticleName, t_SliceNbr, t_SlicePosZ, t_EnergyLossPerSlice);
- m_Data.Set(t_Energy, t_Time, t_Level, t_ParticleName, t_SlicePosZ, t_EnergyLossPerSlice);
+ it->SetStepPosZ(step_posZ);
+ it->SetEnergyLossPerStep(t_Energy);
+ it->SetStepTime(t_Time);
+ }
+ else {
+ // first step: initialize a new m_Data vector
+ // t_Time [ns] is the time at first step
+ m_Data.Set(t_Energy, t_Time, t_Level, t_ParticleName, t_StepPosZ, t_EnergyLossPerStep, t_StepTime);
}
+
return TRUE;
}
@@ -119,9 +127,9 @@ void PS_GaseousDetector::clear() {
m_Data.clear();
t_Level.clear();
t_ParticleName.clear();
- //t_SliceNbr.clear();
- t_SlicePosZ.clear();
- t_EnergyLossPerSlice.clear();
+ t_StepPosZ.clear();
+ t_EnergyLossPerStep.clear();
+ t_StepTime.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Scorers/GaseousDetectorScorers.hh b/NPSimulation/Scorers/GaseousDetectorScorers.hh
index 874613dc619c7ce6362bc2ce90a63b1fb0be2a5d..fa5a8d2e6babc7f345cacc24fc55b97426f87b70 100644
--- a/NPSimulation/Scorers/GaseousDetectorScorers.hh
+++ b/NPSimulation/Scorers/GaseousDetectorScorers.hh
@@ -8,55 +8,60 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * Original Author: Audrey Chatillon, contact address: audrey.chatillon@cea.fr *
* *
- * Creation Date : February 2013 *
+ * Creation Date : January 2025 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* *
*---------------------------------------------------------------------------*
- * Comment: *** BASED ON CalorimeterScorers *** *
+ * Comment: *** COMPLETELY BASED ON AND COPIED FROM CalorimeterScorers *** *
+ * *** from Adrien MATTA, contact adress matta@lpccaen.in2p3.f *** *
*****************************************************************************/
#include "G4VPrimitiveScorer.hh"
#include "NPSHitsMap.hh"
-//#include "NPSecondaries.hh"
#include <map>
using namespace std;
using namespace CLHEP;
namespace GaseousDetectorScorers {
- // Hold One hit info
+
class GaseousDetectorData{
+ // --- ------------------------------------------- --- //
+ // --- Hold info of one hit ---------------------- --- //
+ // --- one hit = all the interaction steps of the --- //
+ // --- particle in the SensitiveVolume --- //
+ // --- ------------------------------------------- --- //
public:
GaseousDetectorData(const double& Energy,
const double& Time,
const vector<unsigned int>& Nesting,
const vector<string>& name,
- //const vector<int>& Slice,
const vector<double>& posZ,
- const vector<double>& DE){
+ const vector<double>& DE,
+ const vector<double>& Tstep){
m_Index=CalculateIndex(Nesting);
m_Level=Nesting;
m_Energy=Energy;
m_Time=Time;
m_ParticleName=name;
- //m_SliceNbr=Slice;
- m_SlicePosZ=posZ;
- m_EnergyLossPerSlice=DE;
+ m_StepPosZ=posZ;
+ m_EnergyLossPerStep=DE;
+ m_StepTime=Tstep;
};
~GaseousDetectorData(){};
private:
- unsigned int m_Index;
+ unsigned int m_Index;
vector<unsigned int> m_Level;
- double m_Energy;
- double m_Time;
+ double m_Energy;
+ double m_Time;
vector<string> m_ParticleName;
- //vector<int> m_SliceNbr;
- vector<double> m_SlicePosZ;
- vector<double> m_EnergyLossPerSlice;
+ vector<double> m_StepPosZ;
+ vector<double> m_EnergyLossPerStep;
+ vector<double> m_StepTime;
public:
static unsigned int CalculateIndex(const vector<unsigned int>& Nesting);
@@ -67,21 +72,26 @@ namespace GaseousDetectorScorers {
inline double GetEnergy() const {return m_Energy;}
inline double GetTime() const {return m_Time;}
inline vector<string> GetParticleName() const {return m_ParticleName;}
- //inline vector<int> GetSliceNbr() const {return m_SliceNbr;}
- inline vector<double> GetSlicePosZ() const {return m_SlicePosZ;}
- inline vector<double> GetEnergyLossPerSlice() const {return m_EnergyLossPerSlice;}
+ inline vector<double> GetStepPosZ() const {return m_StepPosZ;}
+ inline vector<double> GetEnergyLossPerStep() const {return m_EnergyLossPerStep;}
+ inline vector<double> GetStepTime() const {return m_StepTime;}
public:
- void Add(const double& Energy) {m_Energy+=Energy;};
+ // to be used in the PS_GaseousDetector::ProcessHits(G4Step*, G4TouchableHistory*)
+ // add the current step energy loss to the total energy of the particle
+ void Add(const double& DE) {m_Energy+=DE;};
+ // add specific information at each step
inline void SetParticleName(const string& name){m_ParticleName.push_back(name);}
- //inline void SetSliceNbr(const int& SliceNbr){m_SliceNbr.push_back(SliceNbr);}
- inline void SetSlicePosZ(const double& z){m_SlicePosZ.push_back(z);}
- inline void SetEnergyLossPerSlice(const double& DE){m_EnergyLossPerSlice.push_back(DE);}
- };
+ inline void SetStepPosZ(const double& z){m_StepPosZ.push_back(z);}
+ inline void SetEnergyLossPerStep(const double& DE){m_EnergyLossPerStep.push_back(DE);}
+ inline void SetStepTime(const double& z){m_StepTime.push_back(z);}
+ };
- // Manage a vector of DSSD hit
class GaseousDetectorDataVector{
+ // --- ----------------------- --- //
+ // --- Manage a vector of hits --- //
+ // --- ----------------------- --- //
public:
GaseousDetectorDataVector(){};
~GaseousDetectorDataVector(){};
@@ -100,18 +110,17 @@ namespace GaseousDetectorScorers {
const double& Time,
const vector<unsigned int>& Nesting,
const vector<string>& name,
- //const vector<int>& Slice,
const vector<double>& z,
- const vector<double>& DE) {
- //m_Data.push_back(GaseousDetectorData(Energy,Time,Nesting,name,Slice,z,DE));
- m_Data.push_back(GaseousDetectorData(Energy,Time,Nesting,name,z,DE));
+ const vector<double>& DE,
+ const vector<double>& t) {
+ m_Data.push_back(GaseousDetectorData(Energy,Time,Nesting,name,z,DE,t));
};
const GaseousDetectorData* operator[](const unsigned int& i) const {return &m_Data[i];};
};
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- class PS_GaseousDetector : public G4VPrimitiveScorer{
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ class PS_GaseousDetector : public G4VPrimitiveScorer{
public: // with description
PS_GaseousDetector(G4String name, vector<G4int> NestingLevel,G4int depth=0);
@@ -119,6 +128,7 @@ namespace GaseousDetectorScorers {
protected: // with description
+ // ProcessHists is calles at each step of the simulation
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
public:
@@ -135,20 +145,23 @@ namespace GaseousDetectorScorers {
private:
GaseousDetectorDataVector m_Data;
- double t_Energy;
- double t_Time;
- vector<unsigned int> t_Level;
+ double t_Energy; // t_Energy is incremented at each step
+ double t_Time; // time at first step
+ // all these vectors are used to accumulate info at each step
vector<string> t_ParticleName;
- //vector<int> t_SliceNbr;
- vector<double> t_SlicePosZ;
- vector<double> t_EnergyLossPerSlice;
+ vector<double> t_StepPosZ;
+ vector<double> t_EnergyLossPerStep;
+ vector<double> t_StepTime;
+ // unique vector which is cleared at each step
+ vector<unsigned int> t_Level;
public:
inline unsigned int GetMult() {return m_Data.size();};
inline double GetEnergy(const unsigned int& i) {return m_Data[i]->GetEnergy();};
inline double GetTime(const unsigned int& i) {return m_Data[i]->GetTime();};
inline vector<string> GetParticleName(const unsigned int& i) const {return m_Data[i]->GetParticleName();}
- inline vector<double> GetSlicePosZ(const unsigned int& i) const {return m_Data[i]->GetSlicePosZ();}
- inline vector<double> GetEnergyLossPerSlice(const unsigned int& i) const {return m_Data[i]->GetEnergyLossPerSlice();}
+ inline vector<double> GetStepPosZ(const unsigned int& i) const {return m_Data[i]->GetStepPosZ();}
+ inline vector<double> GetEnergyLossPerStep(const unsigned int& i) const {return m_Data[i]->GetEnergyLossPerStep();}
+ inline vector<double> GetStepTime(const unsigned int& i) const {return m_Data[i]->GetStepTime();}
inline vector<unsigned int> GetLevel(const unsigned int& i) {return m_Data[i]->GetLevel();};
};
}
diff --git a/Projects/AlPhaPha/2024/Snakefile b/Projects/AlPhaPha/2024/Snakefile
index c2b42c1666d746178afc26f844ea7121daa18566..374a37903e8b95991d85c8577260d2d86afebc5d 100644
--- a/Projects/AlPhaPha/2024/Snakefile
+++ b/Projects/AlPhaPha/2024/Snakefile
@@ -3,7 +3,7 @@ import subprocess
# Lire le répertoire d'entrée depuis les arguments de configuration
#input_directory = config["folder"]
-input_directory = os.getcwd() + "/../DataMacro/output/run_247"
+input_directory = os.getcwd() + "/../DataMacro/output/run_246"
origin = []
# Iterate over files in input_directory
for filename in os.listdir(input_directory):
@@ -12,14 +12,14 @@ for filename in os.listdir(input_directory):
origin.append(filename)
# Définir le répertoire de sortie pour les fichiers convertis
-phy_directory = os.getcwd() + "/../DataMacro/output/analysis/run_247"
+phy_directory = os.getcwd() + "/../DataMacro/output/analysis/run_246"
#phy_directory = "./"
# define target files directory
analysedfile = []
for inputfile in origin:
#analysedfile.append("/home/morfouacep/Physics/NPTool/nptool/Projects/ana_e850/root/analysis/"+inputfile.replace("_raw_","_"))
- analysedfile.append( os.getcwd() + "/../DataMacro/output/analysis/run_247/"+inputfile)
+ analysedfile.append( os.getcwd() + "/../DataMacro/output/analysis/run_246/"+inputfile)
## batch rules
rule all:
diff --git a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
index 7ead632c0d8a841894512a0d7f58ade995a08f5e..526de6830dad5f03f9665d95bbb4ff429da0db62 100755
--- a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
+++ b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
@@ -3,4 +3,4 @@ echo "- executing snakemake file for npanalysis..."
snakemake --cores 30 --forceall --keep-incomplete --keep-going --rerun-incomplete
echo "- snakemake executed successfully!"
echo "- Merging file..."
-root -q '../DataMacro/Merger.C(30,"root/analysis","VamosCalib247","../DataMacro/output/analysis/run_247/run_raw_247_")'
+root -q '../DataMacro/Merger.C(30,"root/analysis","Run246","../DataMacro/output/analysis/run_246/run_raw_246_")'
diff --git a/Projects/AlPhaPha/2024/macro/PLOT/PlotPres26JUIN.cxx b/Projects/AlPhaPha/2024/macro/PLOT/PlotPres26JUIN.cxx
index 0b83e3751177bd1b673d5fb4cdc0d3dcc1916a0e..822aacfa65a22cd55042c37e9f0eaf680648942b 100644
--- a/Projects/AlPhaPha/2024/macro/PLOT/PlotPres26JUIN.cxx
+++ b/Projects/AlPhaPha/2024/macro/PLOT/PlotPres26JUIN.cxx
@@ -6,7 +6,7 @@
void PlotPres26JUIN(){
- TFile *f = new TFile("../../root/analysis/Exc240Linea.root");
+ TFile *f = new TFile("../../root/analysis/VamosCalib247.root");
TFile *cutfile = TFile::Open("cut12C.root");
TCut *C12 =(TCut*)cutfile->Get("C12");
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
index 2f5609df9ad4428c9063c1e4e9bff464af363e9c..4225a9fd4bbfc0e62eff7e64bb9614221384114f 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
@@ -5,8 +5,8 @@
#include <unordered_map>
// Important variable
-#define STEP 3; // Step to parse through the histogram in X
-#define BINSIZE 10; // Size of the bin in X on which to projectY the plot
+#define STEP 2; // Step to parse through the histogram in X
+#define BINSIZE 15; // Size of the bin in X on which to projectY the plot
//Main
@@ -59,7 +59,13 @@ void CutAutoDEE(){
for (auto cut : CutLine){
cut->Draw("SAME");
}
-
+
+ //Save the cut !
+ TFile *ofile = new TFile("Output/CutAutoZ.root","RECREATE");
+ for (auto cut : CutLine){
+ cut->Write();
+ }
+ ofile->Close();
}
/**
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
index 860210dd763c14d12d60efe98dec99dde0c57f72..461cd1e6158683dc4571d3bebfe75565f6ec6d51 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
@@ -27,6 +27,7 @@ void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
TChain* chain = new TChain("PhysicsTree");
chain->Add("../../../root/analysis/VamosCalib247.root");
+ chain->Add("../../../root/analysis/Run246.root");
TICPhysics* IC = new TICPhysics() ;
TTimeData *Time = new TTimeData() ;
@@ -440,13 +441,62 @@ void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
gStyle->SetPalette(kRainBow);
+ // =========== Poster plot=============================
+
+ TCanvas* c8 = new TCanvas("DE Before correction","DE Before correction",4000,4000);
+
+ hDE_E->GetXaxis()->SetTitle("E (A.U)");
+ hDE_E->GetXaxis()->SetRangeUser(3000,35000);
+ hDE_E->GetXaxis()->CenterTitle();
+
+ hDE_E->GetYaxis()->SetTitle("#Delta E (A.U)");
+ hDE_E->GetYaxis()->SetRangeUser(6000,26000);
+ hDE_E->GetYaxis()->CenterTitle();
+
+ hDE_E->GetZaxis()->SetRangeUser(0,140);
+ hDE_E->GetZaxis()->SetTitle("Number of events");
+ hDE_E->GetZaxis()->CenterTitle();
+
+ hDE_E->SetTitle("PID VAMOS before correction");
+ hDE_E->SetStats(0);
+ hDE_E->Draw("COLZ");
+
+ c8->SetCanvasSize(4000,3000);
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.png");
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.pdf");
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.jpg");
+
+ TCanvas* c9 = new TCanvas("DE after correction","DE after correction",1000,1000);
+
+ hDE0234_E_corr->GetXaxis()->SetTitle("E (A.U)");
+ hDE0234_E_corr->GetXaxis()->SetRangeUser(3000,35000);
+ hDE0234_E_corr->GetXaxis()->CenterTitle();
+
+ hDE0234_E_corr->GetYaxis()->SetTitle("#Delta E (A.U)");
+ hDE0234_E_corr->GetYaxis()->SetRangeUser(6000,26000);
+ hDE0234_E_corr->GetYaxis()->CenterTitle();
+
+ hDE0234_E_corr->GetZaxis()->SetRangeUser(0,140);
+ hDE0234_E_corr->GetZaxis()->SetTitle("Number of events");
+ hDE0234_E_corr->GetZaxis()->CenterTitle();
+
+ hDE0234_E_corr->SetTitle("PID VAMOS after correction");
+ hDE0234_E_corr->SetStats(0);
+ hDE0234_E_corr->Draw("COLZ");
+
+ c9->SetCanvasSize(4000,3000);
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.png");
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.pdf");
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.jpg");
+ //gStyle->SetOptStat(1);
+ // =========== Poster plot=============================
+
+
TCanvas* c7 = new TCanvas("correction succesive","correction succesive",1500,1000);
c7->Divide(6);
c7->cd(1);
- hDE_E->GetXaxis()->SetTitle("E");
- hDE_E->GetYaxis()->SetTitle("DE");
hDE_E->Draw();
c7->cd(2);
@@ -466,7 +516,7 @@ void DECorr(bool cut = false, bool spline = false) {
c7->cd(5);
hDE_E0234->GetXaxis()->SetTitle("E");
- hDE_E0234->GetYaxis()->SetTitle("DE");
+ //hDE_E0234->GetYaxis()->SetTitle("DE");
hDE_E0234->Draw();
c7->cd(6);
@@ -548,8 +598,6 @@ void DECorr(bool cut = false, bool spline = false) {
gPad->Modified();
gPad->Update();
hDE0234_E_corr->Draw("colz");
- hDE0234_E_corr->GetXaxis()->SetTitle("E");
- hDE0234_E_corr->GetYaxis()->SetTitle("DE");
c2->cd(6);
hDE_Ebis->SetTitle("online");
@@ -565,7 +613,6 @@ void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
TFile *oSaved = new TFile("Output/DE_E.root","recreate");
- hDE0234_E_corr->GetXaxis()->SetRangeUser(2000,34000);
hDE0234_E_corr->Write();
TH2toCSV(hDE0234_E_corr,"Output/DE_E.csv");
diff --git a/Projects/Epic_sim_proto/K1_238Pu.source b/Projects/Epic_sim_proto/K1_238Pu.source
new file mode 100755
index 0000000000000000000000000000000000000000..c8988161b1dc8ee3357c9159a1efe9d449bcdacf
--- /dev/null
+++ b/Projects/Epic_sim_proto/K1_238Pu.source
@@ -0,0 +1,26 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy are given in MeV , Position in mm %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Isotropic
+ EnergyLow= 0 MeV
+ EnergyHigh= 0 MeV
+ HalfOpenAngleMin= 0 deg
+ HalfOpenAngleMax= 180 deg
+ x0= 0
+ y0= 0
+ z0= -14.1929 mm
+ SourceProfile= Flat
+ SigmaX= 16.5 mm
+ SigmaY= 0 mm
+ SigmaZ= 0 mm
+ particle= 238Pu
+ ExcitationEnergy= 0 MeV
+ Direction= z
+ Multiplicity= 10
+% DecayLaw= on: multiplicity following the decay law controlled by the halflife
+% off: fix multiplcitity
+ DecayLaw= on
+ ActivityBq= 1.e7
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/Epic_sim_proto/K1_238U.source b/Projects/Epic_sim_proto/K1_238U.source
index 00b7a26dd864bda460d89e93809680ff2f8c1f9e..0edffac973f560f8bbfe5caadf482bfd8ca95ee7 100755
--- a/Projects/Epic_sim_proto/K1_238U.source
+++ b/Projects/Epic_sim_proto/K1_238U.source
@@ -11,7 +11,16 @@ Isotropic
x0= 0
y0= 0
z0= -14.1929 mm
+ SourceProfile= Flat
+ SigmaX= 16.5 mm
+ SigmaY= 0 mm
+ SigmaZ= 0 mm
particle= 238U
ExcitationEnergy= 0 MeV
Direction= z
+ Multiplicity= 10
+% DecayLaw= on: multiplicity following the decay law controlled by the halflife
+% off: fix multiplcitity
+ DecayLaw= on
+ ActivityBq= 50
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/Epic_sim_proto/info.dat b/Projects/Epic_sim_proto/info.dat
new file mode 100644
index 0000000000000000000000000000000000000000..1287e8e878d5a7cafb2e2402cbb5ac1a745ecca5
--- /dev/null
+++ b/Projects/Epic_sim_proto/info.dat
@@ -0,0 +1,38 @@
+%%%%%%%%%%%%%%%%
+% 238U samples %
+%%%%%%%%%%%%%%%%
+ InputDataFile= /home/audrey/EPIC/GEF_lmd_238Unf/En_10MeV_CN_239U.lmd
+ K1: z0= -14.1929 mm
+ K2: z0= -9.17507 mm
+ K5: z0= -2.50894 mm
+ K6: z0= 2.50894 mm
+ K9: z0= 9.17506 mm
+ K10: z0= 14.1929 mm
+
+%%%%%%%%%%%%%%%%
+% 238U samples %
+%%%%%%%%%%%%%%%%
+ InputDataFile= /home/audrey/EPIC/GEF_lmd_235Unf/En_10MeV_CN_236U.lmd
+ K3: z0= -8.35094 mm
+ K4: z0= -3.33307 mm
+ K7: z0= 3.33307 mm
+ K8: z0= 3.33307 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% COMMANDS FOR SIMULATION %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ npsimulation -D Epic.detector -E K1_238U.source -M U238.mac --cut-parent-id 1 -O K1_alpha238U.root
+ npsimulation -D Epic.detector -E GEF.source --cut-parent-id 0 -O K1_GEF.root
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% COMMANDS TO READ EPIC DATA %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ TFile * falpha = new TFile("K1_alpha238U.root","read")
+ TTree * talpha = (TTree*)falpha->Get("SimulatedTree")
+ talpha->Draw("Epic.fEpic_Data.Q1>>h_alpha_Q1(3000,0,150)")
+
+ TFile * fGEF = new TFile("K1_GEF.root","read")
+ TTree * tGEF = (TTree*)fGEF->Get("SimulatedTree")
+ tGEF->Draw("Epic.fEpic_Data.Q1>>h_FF_Q1(3000,0,150)","","sames")
+
+