Merge branch 'new_ecal_hists' into 'master'

Adding histograms and instructions for ECAL simulation

See merge request !143

Control/CalypsoExample/SimHitExample/CaloREADME.md

+
+# Running the ECAL simulation  
+
+Use runEcal.py to specify particle type, energy, mass, z position in calo, theta and phi ranges etc
+runEcal.py location: calypso/Simulation/G4Faser/G4FaserAlg/test/runEcal.py
+ - ConfigFlags.Output.HITSFileName = ""
+specifies the output name of the HITS file
+ - pg.sampler.pid = -11
+particle ID (-11 = electron)
+ - pg.sampler.mom = PG.EThetaMPhiSampler(energy=100*GeV, theta=[0,pi/20], phi=[0,2*pi], mass=0.511)
+set particle energy, theta, phi and mass
+ - pg.sampler.pos = PG.PosSampler(x=[-5.0, 5.0], y=[-5.0, 5.0], z=2744.1, t=0.0)
+set position of beam in terms of x and y coordinates 
+z = 2744.1 (mm) corresponds to calorimeter face 
+x=60.9 and y=49.6 corresponds to centre of top left module
+ - sys.exit(int(acc.run(maxEvents=1000).isFailure()))
+set the number of events to simulate
+1000 event simulations at 1 TeV for electrons and photons need to be split
+change the random seed in this case
+ - pg.randomSeed = 123456
+
+Run this file in Calypso's run directory to genetate HITS file of simulation.
+
+# Plotting Histograms in CaloSimHitAlg
+
+Once the HITS file has been generated from runEcal.py:
+ - use CaloSimHitExample_jobOptions.py to loop over HITS and generate HistoFile
+ - the histogram plotting is specified in CaloSimHitAlg.h and CaloSimHitAlg.cxx
+
+CaloSimHitAlg.cxx plots variations of 6 histograms:
+ - eLoss - energy deposited per hit in the 4 ECAL modules (GeV)
+    - eLossBR - eLoss in bottom right ECAL module
+    - eLossTR - eLoss in top right ECAL module
+    - eLossBL - eLoss in bottom left ECAL module
+    - eLossTL - eLoss in top left ECAL module
+ - eLossTot - total energy deposited by particle across all 4 ECAL modules (GeV)
+    - eLossBRTot - eLossTot in bottom right ECAL module
+    - eLossTRTot - eLossTot in top right ECAL module
+    - eLossBLTot - eLossTot in bottom left ECAL module
+    - eLossTLTot - eLossTot in top left ECAL module
+ - module - cross section of hits in the calo, where (0,0), (1,0), (0,1) and (1,1) correspond to ECAL module centres
+ - modulePos - cross section of hits in the calo, showing their position in mm
+ - meanTime - events in the calo as a function of mean time (ns)
+ - weightedTime - energy weighted events in the calo as a function of mean time (ns)
+
+There are also fractional histograms for eLoss and eLossTot that give the deposited energy as a fraction of injected energy.
+Due to the histograms not depending on injected energy of simulated particle, may need to Rebin() when drawing eLoss and eLossTot histograms.

Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx

@@ -14,22 +14,46 @@ StatusCode CaloSimHitAlg::initialize()
     // initialize a histogram 
     // letter at end of TH1 indicated variable type (D double, F float etc)
     // first string is root object name, second is histogram title
+ 
+
+  m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss;Energy Loss [GeV];Events", 1000000, 0, 100); // 100000 bins, can rebin if needed - 1000 GeV inj - max 100 per event 
+    m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Top Left Module;Energy Loss [GeV];Events", 100000, 0, 100);
+    m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Top Right Module;Energy Loss [GeV];Events", 100000, 0, 100);
+    m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Bottom Right Module;Energy Loss [GeV];Events", 100000, 0, 100);
+    m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Bottom Left Module;Energy Loss [GeV];Events", 100000, 0, 100);
+// 1 million bins
+    m_elossTot = new TH1D("eLossTot", "Total Energy Deposited in Calorimeter;Deposited Energy [GeV];Events", 1000000, 0, 1000);   
+    m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Deposited in Top Left Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+    m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Deposited in Top Right Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+    m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Deposited in Bottom Right Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+    m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Deposited in Bottom Left Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);    
 
-    m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss Fraction;Fraction;Events", 1000, 0, 0.0001); 
-    m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Fraction Top Left Module;Fraction;Events", 1000, 0, 0.0001);
-    m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Fraction Top Right Module;Fraction;Events", 1000, 0, 0.0001);
-    m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Fraction Bottom Right Module;Fraction;Events", 1000, 0, 0.0001);
-    m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Fraction Bottom Left Module;Fraction;Events", 1000, 0, 0.0001);
+    m_module = new TH2D("module;x;y", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
+    // m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
+    // x4 original number of bins
+    m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 52, -130, 130, 52, -130, 130 );
 
-    m_elossTot = new TH1D("eLossTot", "Total Energy Fraction Deposited in Calorimeter;Fraction;Events", 1000, 0, 1);   
-    m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Fraction Deposited in Top Left Module;Fraction;Events", 1000, 0, 1);
-    m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Fraction Deposited in Top Right Module;Fraction;Events", 1000, 0, 1);
-    m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Fraction Deposited in Bottom Right Module;Fraction;Events", 1000, 0, 1);
-    m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Fraction Deposited in Bottom Left Module;Fraction;Events", 1000, 0, 1);    
+    // as a fraction of injected energy (when known)
 
-    m_module = new TH2D("module;x;y", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
-    m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
-        
+    m_eloss_frac = new TH1D("eLoss_frac", "Calo Hit Energy Loss Fraction;Fraction;Events", 100000, 0, 0.1); 
+    m_elossTL_frac = new TH1D("eLossTL_frac", "Calo Hit Energy Loss Fraction Top Left Module;Fraction;Events", 100000, 0, 0.01);
+    m_elossTR_frac = new TH1D("eLossTR_frac", "Calo Hit Energy Loss Fraction Top Right Module;Fraction;Events", 100000, 0, 0.01);
+    m_elossBR_frac = new TH1D("eLossBR_frac", "Calo Hit Energy Loss Fraction Bottom Right Module;Fraction;Events", 100000, 0, 0.01);
+    m_elossBL_frac = new TH1D("eLossBL_frac", "Calo Hit Energy Loss Fraction Bottom Left Module;Fraction;Events", 100000, 0, 0.01);
+
+    m_elossTot_frac = new TH1D("eLossTot_frac", "Total Energy Fraction Deposited in Calorimeter;Fraction;Events", 1000000, 0, 1);   
+    m_elossTLTot_frac = new TH1D("eLossTLTot_frac", "Total Energy Fraction Deposited in Top Left Module;Fraction;Events", 1000000, 0, 1);
+    m_elossTRTot_frac = new TH1D("eLossTRTot_frac", "Total Energy Fraction Deposited in Top Right Module;Fraction;Events", 1000000, 0, 1);
+    m_elossBRTot_frac = new TH1D("eLossBRTot_frac", "Total Energy Fraction Deposited in Bottom Right Module;Fraction;Events", 1000000, 0, 1);
+    m_elossBLTot_frac = new TH1D("eLossBLTot_frac", "Total Energy Fraction Deposited in Bottom Left Module;Fraction;Events", 1000000, 0, 1);  
+
+    // new histograms hits vs time 29/6/21
+
+    m_meanTime = new TH1D("meanTime", "Events vs Mean Time;Time [ns];Events",4000, 0, 2);
+    m_weightedTime = new TH1D("weightedTime", "Weighted Events vs Time;Time [ns];Events",4000, 0, 2); // weighted energy should be in GeV
+
+    ATH_CHECK(histSvc()->regHist("/HIST/meanTime", m_meanTime));
+    ATH_CHECK(histSvc()->regHist("/HIST/weightedTime", m_weightedTime));     
 
     ATH_CHECK(histSvc()->regHist("/HIST/eloss", m_eloss));
     ATH_CHECK(histSvc()->regHist("/HIST/elossTL", m_elossTL)); 
@@ -42,10 +66,24 @@ StatusCode CaloSimHitAlg::initialize()
     ATH_CHECK(histSvc()->regHist("/HIST/elossTRTot", m_elossTRTot));
     ATH_CHECK(histSvc()->regHist("/HIST/elossBRTot", m_elossBRTot));
     ATH_CHECK(histSvc()->regHist("/HIST/elossBLTot", m_elossBLTot));
-
+ 
     ATH_CHECK(histSvc()->regHist("/HIST/modules", m_module));
     ATH_CHECK(histSvc()->regHist("/HIST/modulePos", m_modulePos));
 
+    // as a fraction of injected energy 
+
+    ATH_CHECK(histSvc()->regHist("/HIST/eloss_frac", m_eloss_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossTL_frac", m_elossTL_frac)); 
+    ATH_CHECK(histSvc()->regHist("/HIST/elossTR_frac", m_elossTR_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossBR_frac", m_elossBR_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossBL_frac", m_elossBL_frac));
+
+    ATH_CHECK(histSvc()->regHist("/HIST/elossTot_frac", m_elossTot_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossTLTot_frac", m_elossTLTot_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossTRTot_frac", m_elossTRTot_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossBRTot_frac", m_elossBRTot_frac));
+    ATH_CHECK(histSvc()->regHist("/HIST/elossBLTot_frac", m_elossBLTot_frac));
+
     // initialize data handle keys
     ATH_CHECK( m_mcEventKey.initialize() );
     ATH_CHECK( m_faserCaloHitKey.initialize() );
@@ -78,8 +116,18 @@ StatusCode CaloSimHitAlg::execute()
     if (ev->particles_size() > 0)
     {
         HepMC::GenEvent::particle_const_iterator p = ev->particles_begin();
-        ePrimary = (*p)->momentum().e();
+        //ePrimary = (*p)->momentum().e();
+
+	for ( ; p !=  ev->particles_end(); ++p) {
+	  // Particle gun primaries have BCs starting at 10000, while others have 3 digit values
+	  if ((*p)->barcode() > 11000) continue;
+	  ePrimary += (*p)->momentum().e();
+	}
+	std::cout << "E Primary " << ePrimary << std::endl;
     } 
+    
+    //ev->print();
+    //ePrimary = 1000000;
 
     // The hit container might be empty because particles missed the modules
     if (h_caloHits->size() == 0) return StatusCode::SUCCESS;
@@ -91,18 +139,38 @@ StatusCode CaloSimHitAlg::execute()
     // BL = bottom left calorimeter module
     // BR = bottom right calorimeter module 
 
-    float e(0);;
+    float e(0);
     float tl(0);
     float tr(0);
     float br(0);
     float bl(0);
 
+    float e_frac(0);
+    float tl_frac(0);
+    float tr_frac(0);
+    float br_frac(0);
+    float bl_frac(0);
+
+  
+
     for (const CaloHit& hit : *h_caloHits)  {
       //hit.print();
-      float ehit = hit.energyLoss()/ePrimary;
+      float ehit_frac = hit.energyLoss()/ePrimary;
+      float ehit = hit.energyLoss()/1000; // to convert from MeV to GeV
+      m_eloss_frac->Fill(ehit_frac);
       m_eloss->Fill(ehit);
+      e_frac += ehit_frac;
       e += ehit;
 
+      // add extra time hist
+      float time = hit.meanTime();
+      m_meanTime->Fill(time);
+
+      // add weighted energy hist
+      m_weightedTime->Fill(time,ehit);
+
+      //
+
       m_module->Fill(hit.getModule(), hit.getRow());	
       m_modulePos->Fill(hit.getModule()==0 ? -121.9/2 + hit.localStartPosition().x() : 121.9/2 + hit.localStartPosition().x(), hit.getRow()==0 ? -121.9/2 + hit.localStartPosition().y() : 121.9/2 + hit.localStartPosition().y());
 
@@ -110,21 +178,29 @@ StatusCode CaloSimHitAlg::execute()
       if (hit.getModule() == 0 && hit.getRow() == 1) // 0 1 TL
 	{
 	  m_elossTL->Fill(ehit);
+	  m_elossTL_frac->Fill(ehit_frac);
+	  tl_frac += ehit_frac;
 	  tl += ehit;
 	}
       if (hit.getModule() == 1 && hit.getRow() == 1) // 1 1 TR
 	{
 	  m_elossTR->Fill(ehit);
+	  m_elossTR_frac->Fill(ehit_frac);
+	  tr_frac += ehit_frac;
 	  tr += ehit;
 	}
       if (hit.getModule() == 1 && hit.getRow() == 0) // 1 0 BR
 	{
 	  m_elossBR->Fill(ehit);
+	  m_elossBR_frac->Fill(ehit_frac);
+	  br_frac += ehit_frac;
 	  br += ehit;
 	}
       if (hit.getModule() == 0 && hit.getRow() == 0) // 0 0 BL
 	{
 	  m_elossBL->Fill(ehit);
+	  m_elossBL_frac->Fill(ehit_frac);
+	  bl_frac += ehit_frac;
 	  bl += ehit;
 	}
 
@@ -136,6 +212,14 @@ StatusCode CaloSimHitAlg::execute()
     m_elossBRTot->Fill(br);
     m_elossBLTot->Fill(bl); 
 
+    // as a fraction of injected energy 
+
+    m_elossTot_frac->Fill(e_frac);
+    m_elossTLTot_frac->Fill(tl_frac);
+    m_elossTRTot_frac->Fill(tr_frac);
+    m_elossBRTot_frac->Fill(br_frac);
+    m_elossBLTot_frac->Fill(bl_frac);
+
     return StatusCode::SUCCESS;
 }
 

Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h

@@ -21,6 +21,9 @@ class CaloSimHitAlg : public AthHistogramAlgorithm
     TH2* m_module;
     TH2* m_modulePos;
     TH1* m_elossTot; 
+ 
+    TH1* m_meanTime;
+    TH1* m_weightedTime;
 
     // add new histograms
     // TL = top left calorimeter module
@@ -35,7 +38,20 @@ class CaloSimHitAlg : public AthHistogramAlgorithm
     TH1* m_elossTLTot;
     TH1* m_elossTRTot;
     TH1* m_elossBRTot;
-    TH1* m_elossBLTot;   
+    TH1* m_elossBLTot;    
+
+    // add in fractional histograms
+
+    TH1* m_eloss_frac;
+    TH1* m_elossTL_frac;
+    TH1* m_elossTR_frac;
+    TH1* m_elossBR_frac; 
+    TH1* m_elossBL_frac; 
+    TH1* m_elossTot_frac;
+    TH1* m_elossTLTot_frac;
+    TH1* m_elossTRTot_frac;
+    TH1* m_elossBRTot_frac;
+    TH1* m_elossBLTot_frac;  
 
     // Read handle keys for data containers
     // Any other event data can be accessed identically