diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
index 527d89ef1d0c1f8aaf5bbc40ee893882ebd2dcc1..e7dde52b04aa9b1f5c9bb0fdde95fad2b2a78f00 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
@@ -69,7 +69,6 @@ G4String dimu_xsec_fac = "1000";
int main(int argc,char** argv) {
-
for(int i{0}; i<argc; ++i)
{
if(G4String(argv[i]) == "--nevts")
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
index 9f30f0c97ba2e9c66fe5e8f4854d0e479c012597..15e4f52fb71e2eb1d38fe79ea89219a37e29c53b 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
@@ -59,13 +59,15 @@ import ROOT
_output = ROOT.TFile.Open("G4GammaToDiLeptonConversionTest.root", "NEW")
_histo_dim_dict = {'DiElectron' : {'Hist' : {'InvMass' : (100, 0, 100),
- 'ThetaPlus' : (100, 0, 0.02),
- 'ThetaMinus' : (100, 0, 0.02),
+ 'ThetaSep' : (50, 0, 6E-4),
+ 'InvMass2D' : (100, 1, 3),
+ 'Theta' : (100, 0, 0.02),
'EFrac' : (100,0,1)},
'Symbol' : 'e'},
'DiMuon' : {'Hist' : {'InvMass' : (150, 250, 550),
- 'ThetaPlus' : (100, 0, 0.07),
- 'ThetaMinus' : (100, 0, 0.07),
+ 'Theta' : (100, 0, 0.07),
+ 'ThetaSep' : (50, 0, 1E-1),
+ 'InvMass2D' : (50, 200, 500),
'EFrac' : (100,0,1)},
'Symbol' : '#mu' } }
@@ -85,48 +87,49 @@ for energy in energies:
hist_mm.GetXaxis().SetTitle("M_{%s}/MeV" % '{s}{s}'.format(s=_histo_dim_dict[particle]['Symbol']))
- hist_tp = ROOT.TH1F("{}_AnglePlus_{}GeV_{}mm_Al".format(particle, energy, t),
- "#gamma-{} Angle for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{+}', energy, t),
- *_histo_dim_dict[particle]['Hist']['ThetaPlus'])
- hist_tp.GetXaxis().SetTitle("#theta(p_{%s}, p_{%s})" % ('#gamma', _histo_dim_dict[particle]['Symbol']+'^{+}'))
+ hist_tp = ROOT.TH1F("{}_Angle_{}GeV_{}mm_Al".format(particle, energy, t),
+ "#gamma-{} Angle for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{#pm}', energy, t),
+ *_histo_dim_dict[particle]['Hist']['Theta'])
+ hist_tp.GetXaxis().SetTitle("#theta(p_{%s}, p_{%s})" % ('#gamma', _histo_dim_dict[particle]['Symbol']+'^{#pm}'))
-
- hist_tm = ROOT.TH1F("{}_AngleMinus_{}GeV_{}mm_Al".format(particle, energy, t),
- "#gamma-{} Angle for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{+}', energy, t),
- *_histo_dim_dict[particle]['Hist']['ThetaMinus'])
- hist_tm.GetXaxis().SetTitle("#theta(p_{%s}, p_{%s})" % ('#gamma', _histo_dim_dict[particle]['Symbol']+'^{-}'))
-
- hist_efp = ROOT.TH1F("{}_EnergyFrac_Plus_{}GeV_{}mm_Al".format(particle, energy, t),
- "Energy Fraction of {} for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{+}', energy, t),
+ hist_efp = ROOT.TH1F("{}_EnergyFrac_{}GeV_{}mm_Al".format(particle, energy, t),
+ "Energy Fraction of {} for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{#pm}', energy, t),
*_histo_dim_dict[particle]['Hist']['EFrac'])
- hist_efp.GetXaxis().SetTitle("E_{%s}/E_{%s}" % (_histo_dim_dict[particle]['Symbol']+'^{+}', "#gamma"))
+ hist_efp.GetXaxis().SetTitle("E_{%s}/E_{%s}" % (_histo_dim_dict[particle]['Symbol']+'^{#pm}', "#gamma"))
- hist_efm = ROOT.TH1F("{}_EnergyFrac_Minus_{}GeV_{}mm_Al".format(particle, energy, t),
- "Energy Fraction of {} for {}GeV Photons in {}mm of Aluminium".format(_histo_dim_dict[particle]['Symbol']+'^{-}', energy, t),
- *_histo_dim_dict[particle]['Hist']['EFrac'])
- hist_efm.GetXaxis().SetTitle("E_{%s}/E_{%s}" % (_histo_dim_dict[particle]['Symbol']+'^{-}', "#gamma"))
+ hist_2d_sep_ang_mm = ROOT.TH2F("{}_Separation_Angle_vs_InvMass_{}GeV_{}mm_Al".format(particle, energy, t),
+ "Separation Angle Between {} Pair vs Invariant Mass for {}GeV Photons in {}mm of Aluminium".format(particle, energy, t),
+ _histo_dim_dict[particle]['Hist']['InvMass2D'][0],
+ _histo_dim_dict[particle]['Hist']['InvMass2D'][1],
+ _histo_dim_dict[particle]['Hist']['InvMass2D'][2],
+ _histo_dim_dict[particle]['Hist']['ThetaSep'][0],
+ _histo_dim_dict[particle]['Hist']['ThetaSep'][1],
+ 10*_histo_dim_dict[particle]['Hist']['ThetaSep'][2]/float(energy))
+
+ hist_2d_sep_ang_mm.GetXaxis().SetTitle("M_{%s}/MeV" % '{s}{s}'.format(s=_histo_dim_dict[particle]['Symbol']))
+ hist_2d_sep_ang_mm.GetYaxis().SetTitle("#theta(p_{%s},p_{%s})" % (_histo_dim_dict[particle]['Symbol']+'^{+}',
+ _histo_dim_dict[particle]['Symbol']+'^{-}'))
+ hist_2d_sep_ang_mm.GetYaxis().SetTitleOffset(1.2)
for event in _tree: # ROOT broke when using Draw(hist >> hist_tmp(i,j,k)) method
hist_mm.Fill(getattr(event, '{}_InvMass'.format(particle)))
hist_tp.Fill(getattr(event, '{}_ThetaPlus'.format(particle)))
- hist_tm.Fill(getattr(event, '{}_ThetaMinus'.format(particle)))
+ hist_tp.Fill(getattr(event, '{}_ThetaMinus'.format(particle)))
hist_efp.Fill(getattr(event, '{}_EnergyFracPlus'.format(particle)))
- hist_efm.Fill(getattr(event, '{}_EnergyFracMinus'.format(particle)))
-
+ hist_efp.Fill(getattr(event, '{}_EnergyFracMinus'.format(particle)))
+ hist_2d_sep_ang_mm.Fill(getattr(event, '{}_InvMass'.format(particle)), getattr(event, '{}_ThetaSep'.format(particle)))
_output.cd()
- hist_mm.Write()
- hist_tp.Write()
- hist_tm.Write()
- hist_efp.Write()
- hist_efm.Write()
+ hist_mm.Write(hist_mm.GetName(),ROOT.TObject.kOverwrite)
+ hist_tp.Write(hist_tp.GetName(),ROOT.TObject.kOverwrite)
+ hist_efp.Write(hist_efp.GetName(),ROOT.TObject.kOverwrite)
+ hist_2d_sep_ang_mm.Write(hist_2d_sep_ang_mm.GetName(),ROOT.TObject.kOverwrite)
hist_mm.Delete()
hist_tp.Delete()
- hist_tm.Delete()
hist_efp.Delete()
- hist_efm.Delete()
+ hist_2d_sep_ang_mm.Delete()
_tmp.Close()
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
index a76c886147c84e74113921dd41a1f5a09b848d7e..dfc52a0c86a6bdaf00cf36bb17e66d00ca53a037 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
@@ -111,6 +111,7 @@ void RunAction::BeginOfRunAction(const G4Run* aRun)
fAnalysis->CreateNtupleDColumn("DiMuon_PX_Minus");
fAnalysis->CreateNtupleDColumn("DiMuon_PY_Minus");
fAnalysis->CreateNtupleDColumn("DiMuon_PZ_Minus");
+ fAnalysis->CreateNtupleDColumn("DiMuon_ThetaSep"); //Separation Angle Between DiLeptons
fAnalysis->FinishNtuple();
fAnalysis->CreateNtuple("DiElectron", "DiElectron Production");
@@ -125,6 +126,7 @@ void RunAction::BeginOfRunAction(const G4Run* aRun)
fAnalysis->CreateNtupleDColumn("DiElectron_PX_Minus");
fAnalysis->CreateNtupleDColumn("DiElectron_PY_Minus");
fAnalysis->CreateNtupleDColumn("DiElectron_PZ_Minus");
+ fAnalysis->CreateNtupleDColumn("DiElectron_ThetaSep"); //Separation Angle Between DiLeptons
fAnalysis->FinishNtuple();
G4cout << "\n----> NTuple file is opened in " << fileName << G4endl;
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
index 3ad44d9545de01f1830a41cbadf5340a163c019f..906dbb61f6176bebb991c2e86f67c19411c071ee 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
@@ -94,9 +94,9 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep)
ETot += (*secondary)[lp]->GetTotalEnergy();
PTot += (*secondary)[lp]->GetMomentum();
-
}
-
+
+ G4double thetaSep = Pplus.angle(Pminus);
G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma;
G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus);
G4double GammaPlus = Eplus/particle_mass;
@@ -122,6 +122,7 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep)
analysisManager->FillNtupleDColumn(tuple_index, 8, Pminus.x());
analysisManager->FillNtupleDColumn(tuple_index, 9, Pminus.y());
analysisManager->FillNtupleDColumn(tuple_index, 10, Pminus.z());
+ analysisManager->FillNtupleDColumn(tuple_index, 11, thetaSep);
analysisManager->AddNtupleRow(tuple_index);
}