Update Tile Pulse Simulator to add continuous readout mode

Update Tile Pulse Simulator to add continuous readout mode

TileCalorimeter/TilePulseSimulator/TilePulseSimulator/TileSampleGenerator.h

@@ -23,7 +23,8 @@ public:
 		m_ps = ps;
 	}
 	void fillSamples(double t0, double pedestal, double amplitude1, double amplitude2, TF1* pdf, bool addNoise, double itOffset = 0, double otOffset = 50);
-	void fillNSamples(double t0, double pedestal, double amp_it, std::vector<float> amp_pu, TF1* pdf, bool addNoise, double itOffset = 0, int nSamples = 7, int nPul = 21);
+	void fillNSamples(double t0, double pedestal, double amp_it, const std::vector<float>& amp_pu, TF1* pdf, bool addNoise, double itOffset = 0, int nSamples = 7, int nPul = 21);
+        float fillSample(double t0, double pedestal, const std::vector<float>& amp_pu, TF1* pdf, bool addNoise, int nPul = 21, int gain = 1);
 	void fill7SamplesQIE(float amp_it, float *amp_pu); //The function calculates charges for 7 intervals with each interval of 25 ns
 
 private:

TileCalorimeter/TilePulseSimulator/src/TileSampleGenerator.cxx

@@ -61,9 +61,9 @@ void TileSampleGenerator::fillSamples(double t0, double pedestal, double amplitu
 
 //
 //________________________________________________________
-void TileSampleGenerator::fillNSamples(double t0, double pedestal, double amp_it, vector<float> amp_pu, TF1* pdf, bool addNoise, double itOffset, int nSamples, int nPul) {
+void TileSampleGenerator::fillNSamples(double t0, double pedestal, double amp_it, const std::vector<float>& amp_pu, TF1* pdf, bool addNoise, double itOffset, int nSamples, int nPul) {
 
-    std::unique_ptr<TileSampleBuffer> bufall(new TileSampleBuffer(nPul, -25*((nPul-1)/2), 25.));
+        std::unique_ptr<TileSampleBuffer> bufall = std::make_unique<TileSampleBuffer>(nPul, -25*((nPul-1)/2), 25.);
 
 	if(m_DEBUG){
 		cout << "Pileup pulses:" << std::endl;
@@ -127,6 +127,47 @@ void TileSampleGenerator::fillNSamples(double t0, double pedestal, double amp_it
 	return;
 }
 
+float TileSampleGenerator::fillSample(double t0, double pedestal, const vector<float>& amp_pu, TF1* pdf, bool addNoise, int nPul, int gain) {
+
+        std::unique_ptr<TileSampleBuffer> bufall = std::make_unique<TileSampleBuffer>(nPul, -25*((nPul-1)/2), 25.);
+
+	if(m_DEBUG){
+		std::cout << "Pileup pulses:" << std::endl;
+		for (std::vector<float>::const_iterator i = amp_pu.begin(); i != amp_pu.end(); ++i)
+		std::cout << *i << ' ';
+		std::cout << std::endl;
+	}
+
+	double amp_it_out = pedestal;
+	vector<int> t(nPul);
+
+	for (int pul=0; pul < nPul; pul++){
+
+		t[pul] = bufall->getTime(pul) - t0;
+
+		if(gain == 1){
+		  amp_it_out += m_ps->eval(t[pul], false, true) * amp_pu.at(pul);
+		}
+		else{
+		  amp_it_out += (m_ps->eval(t[pul], false, true) * amp_pu.at(pul)) / 40;
+		}
+
+		if(m_DEBUG){
+			std::cout << "nPul-1-pul: " << pul << " getTime(nPul-1-pul) " << bufall->getTime(pul) << " ps " << m_ps->eval(bufall->getTime(pul), false, true) << std::endl;
+			std::cout << "PU sample " << pul << ", pulse shape evaluated for t'=" << t[pul] << std::endl;
+			std::cout << "Contribution for in time amp " << m_ps->eval(t[pul], false, true) << " * " << amp_pu.at(pul) << std::endl;
+			std::cout << std::endl;
+		}
+	}
+
+	if (addNoise) {
+		amp_it_out += pdf->GetRandom();
+	}
+
+	return amp_it_out;
+
+}
+
 //
 //________________________________________________________
 // Simulation of the PMT pulse shapes for QIE FEB

TileCalorimeter/TileSimAlgs/TileSimAlgs/TileDigitsFromPulse.h

@@ -78,6 +78,7 @@ class IAthRNGSvc;
 
 class TH1F;
 class TFile;
+class TRandom3;
 
 /** 
 @class TileDigitsFromPulse
@@ -152,13 +153,14 @@ private:
     std::string m_ootADistHistName; //!< Name of  histogram for out-of-time amplitude distribution
     
     bool m_simQIE; //!<Raw PMT pulses are generated if the option is set to true. The option is intended to simulate the QIE FEB.
+    bool m_simPulseChain; //!< Simulate continous output of readout for HL-LHC paradigm
 
     int m_seed;
     int m_BunchSpacing; //!< Time between pulses in ms 25, 50 or 75
     int m_nSamples;  //!< number of read out samples
     int m_nPul;  //!< number of pileup pulses 
-    int m_nPul_eff;  //Used for symetrization of PU in computation
-    std::vector<float> m_PUAmp;
+    int m_nPul_eff{};  //Used for symetrization of PU in computation
+    std::vector<std::vector<std::vector<std::vector<float>>>> m_PUAmp; // Used to store PU amplitudes
     bool m_PhaseII; //Use parameters of TilePhaseII if the option is set to true
     bool m_bigain; //If true, save the two gains in the ntuples
 
@@ -186,6 +188,12 @@ private:
     bool makeDist(TFile*& file, TH1F*& hist, const std::string& fileName, const std::string& histName="h_Eopt_hi"); //!< Method to read distribution from file
     bool makeDist(TFile*& file, std::vector<std::vector<TH1F*>>& hists, const std::string& fileName);
 
+    void addPileUp(double &n_inTimeAmp, int gain, int ros, int drawer, int channel); //!< Fill vector with pile-up amplitudes
+    void addPileUpSample(int gain, int ros, int drawer, int channel); //!< Fill only a BC with pile-up amplitude
+
+    float m_sample_tru = 0;
+
+    std::unique_ptr<TRandom3> m_random;
 };
 
 #endif // TILESIMALGS_TILEDIGITSFROMPULSE_H

TileCalorimeter/TileSimAlgs/python/RunTilePulseSim.py

@@ -41,6 +41,15 @@ def TileDigitsFromPulseCfg(flags, **kwargs):
         AmpDistLowerLimit             -- Set all bins lower than this to zero. Default = 135
         InTimeAmpDistHistogramName    -- Name of the histogram to use for in-time amplitude distribution
         OutOfTimeAmpDistHistogramName -- Name of the histogram to use for out-of-time amplitude distribution
+        PedestalValueHG               -- Pedestal in HG if not taken from database
+        PedestalValueLG               -- Pedestal in LG if not taken from database
+        SimulatePileUpWithPoiss       -- Simulate pile-up overlaying signals from distribution
+        AvgMuForPileUpSimulation      -- Average number of pp collisions for pile-up simulation with SimulatePileUpWithPoiss
+        PileUpAmpDistFileName         -- Distribution to simulate pile-up with SimulatePileUpWithPoiss
+        RandomSeed                    -- Random seed for random number generator
+        SimulatePulseChain            -- Simulate continuous output from readout cosidering HL-LHC paradigm
+        Bigain                        -- Save two gains in ntuple
+        NPulses                       -- The number of neighboring bunch crossings (before and after the in-time crossing) whose signals are accounted for when simulating the total contribution to a given bunch crossing
     """
 
     kwargs.setdefault('InTimeAmp', 1000)
@@ -48,8 +57,13 @@ def TileDigitsFromPulseCfg(flags, **kwargs):
     kwargs.setdefault('ImperfectionRms', 0)
     kwargs.setdefault('TilePhaseII', False)
     kwargs.setdefault('NSamples', 7)
+    kwargs.setdefault('NPulses', 21)
+    kwargs.setdefault('Bigain', False)
+    kwargs.setdefault('SimulatePulseChain', False)
 
     PhaseII = kwargs['TilePhaseII']
+    PulseChain = kwargs['SimulatePulseChain']
+
     # PhaseII parameters
     if PhaseII:
         kwargs.setdefault('PedestalValueHG', 100)
@@ -62,6 +76,7 @@ def TileDigitsFromPulseCfg(flags, **kwargs):
 
     kwargs.setdefault('PileUpFraction', 0)
     kwargs.setdefault('AmpDistLowerLimit', 0)
+    kwargs.setdefault('SimulatePileUpWithPoiss', False)
     kwargs.setdefault('AvgMuForPileUpSimulation', 80)
 
     from TileGeoModel.TileGMConfig import TileGMCfg
@@ -77,7 +92,7 @@ def TileDigitsFromPulseCfg(flags, **kwargs):
     kwargs['RndmSvc'] = acc.getPrimaryAndMerge( AthRNGSvcCfg(flags) ).name
 
     # Configure TileInfoLoader to set up number of samples
-    nSamples = kwargs['NSamples']
+    nSamples = kwargs['NSamples'] if not PulseChain else 1
     ADCmax = 4095 if PhaseII else 1023
     ADCmaskValue = 4800 if PhaseII else 2047
     from TileConditions.TileInfoLoaderConfig import TileInfoLoaderCfg
@@ -118,16 +133,23 @@ if __name__ == '__main__':
     parser.add_argument('--phaseII', type=bool, default=False, help='Use parameters of TilePhaseII')
     parser.add_argument('--run', type=int, default=410000, help='Run number')
     parser.add_argument('--save-true-amplitude', action='store_true', help='Save true Tile raw channel amplitude into h2000')
+    parser.add_argument('--pulseChain', type=bool, default=False, help='Simulate continuous output of readout across bunch crossings')
     parser.add_argument('--acr-db', action='store_true', help='Use auto correlation matrix from DB')
 
     args, _ = parser.parse_known_args()
 
     flags.Tile.RunType = TileRunType.PHY
+
     # Set up Tile reconstuction methods
-    flags.Tile.doOpt2 = True
-    flags.Tile.doOptATLAS = True
-    if args.nsamples != 7:
+    if(args.pulseChain): # no reconstruction ran for continuous readout simulation
+        flags.Tile.doOpt2 = False
+        flags.Tile.doOptATLAS = False
         flags.Tile.OfcFromCOOL = False
+    else:
+        flags.Tile.doOpt2 = True
+        flags.Tile.doOptATLAS = True
+        if args.nsamples != 7:
+            flags.Tile.OfcFromCOOL = False
 
     flags.Input.isMC = True
     flags.Input.Files = []
@@ -162,7 +184,7 @@ if __name__ == '__main__':
         cfg.merge(PoolReadCfg(flags))
 
     # =======>>> Configure Tile digits from pulse algorithm
-    cfg.merge( TileDigitsFromPulseCfg(flags, NSamples=args.nsamples, TilePhaseII=args.phaseII) )
+    cfg.merge( TileDigitsFromPulseCfg(flags, NSamples=args.nsamples, TilePhaseII=args.phaseII, SimulatePulseChain=args.pulseChain) )
 
     # =======>>> Configure Tile raw channel maker
     from TileRecUtils.TileRawChannelMakerConfig import TileRawChannelMakerCfg
@@ -178,20 +200,27 @@ if __name__ == '__main__':
     # =======>>> Configure Tile h2000 ntuple production
     ntupleFile = f'{args.outputDirectory}/tile_{runNumber}_{args.outputVersion}.aan.root'
     from TileRec.TileAANtupleConfig import TileAANtupleCfg
+
     cfg.merge( TileAANtupleCfg(flags,
-                               saveTMDB=False,
-                               NSamples=args.nsamples,
-                               TileL2Cnt='',
-                               TileDigitsContainerFlt='',
-                               TileDigitsContainer='TileDigitsCnt',
-                               TileRawChannelContainer='TileRawChannelCnt',
-                               TileRawChannelContainerOpt='TileRawChannelOpt2',
-                               CalibrateEnergy=False,
-                               OfflineUnits=0,
-                               CalibMode=True,
-                               outputFile=ntupleFile) )
-    if args.save_true_amplitude:
-        cfg.getEventAlgo('TileNtuple').TileRawChannelContainerFit = 'TrueAmp'
+                            saveTMDB=False,
+                            TileL2Cnt='',
+                            TileDigitsContainerFlt='',
+                            TileDigitsContainer='TileDigitsCnt',
+                            CalibrateEnergy=False,
+                            OfflineUnits=0,
+                            CalibMode=True,
+                            outputFile=ntupleFile) )
+
+    if args.pulseChain:
+        cfg.getEventAlgo('TileNtuple').NSamples = 1
+        cfg.getEventAlgo('TileNtuple').TileRawChannelContainer = 'TrueAmp'
+    else:
+        cfg.getEventAlgo('TileNtuple').NSamples = args.nsamples
+        cfg.getEventAlgo('TileNtuple').TileRawChannelContainer='TileRawChannelCnt'
+        cfg.getEventAlgo('TileNtuple').TileRawChannelContainerOpt='TileRawChannelOpt2'
+
+        if args.save_true_amplitude:
+            cfg.getEventAlgo('TileNtuple').TileRawChannelContainerFit = 'TrueAmp'
 
     # =======>>> Any last things to do?
     if args.postExec: