Added 2D Histogram of Angle vs MM

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")

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()
 

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;   

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);
  
 }