From 6b6375fd1e237583c528d6ed73ddd986837da279 Mon Sep 17 00:00:00 2001
From: Jochen Meyer <Jochen.Meyer@cern.ch>
Date: Mon, 15 Sep 2014 15:59:17 +0200
Subject: [PATCH] fixing compiler warning (CscCalibAlgs-00-12-07)
---
.../CscCalib/CscCalibAlgs/cmt/cdb.log | 0
.../CscCalib/CscCalibAlgs/cmt/requirements | 65 +
.../CscCalib/CscCalibAlgs/python/BulkRun.py | 203 +++
.../CscCalibAlgs/python/BulkRunFollowup.py | 221 ++++
.../CscCalibAlgs/python/CscCalibQuery.py | 421 ++++++
.../CscCalibAlgs/share/CscCalcPedMon.py | 227 ++++
.../CscCalibAlgs/share/CscCalcSlopeMon.py | 198 +++
.../share/CscCalibCommonOptions.py | 119 ++
.../CscCalibAlgs/share/CscExtractPed.py | 57 +
.../share/cscPulserLinearityCorrection.txt | 128 ++
.../CscCalibAlgs/share/cscWritePedCool.py | 185 +++
.../CscCalib/CscCalibAlgs/src/BipolarFit.cxx | 603 +++++++++
.../CscCalib/CscCalibAlgs/src/BipolarFit.h | 48 +
.../CscCalib/CscCalibAlgs/src/CscCalcPed.cxx | 1111 ++++++++++++++++
.../CscCalib/CscCalibAlgs/src/CscCalcPed.h | 203 +++
.../CscCalibAlgs/src/CscCalcSlope.cxx | 1124 +++++++++++++++++
.../CscCalib/CscCalibAlgs/src/CscCalcSlope.h | 184 +++
.../src/components/CscCalibAlgs_entries.cxx | 14 +
.../src/components/CscCalibAlgs_load.cxx | 3 +
19 files changed, 5114 insertions(+)
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/cdb.log
create mode 100755 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/requirements
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRun.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRunFollowup.py
create mode 100755 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/CscCalibQuery.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcPedMon.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcSlopeMon.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalibCommonOptions.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscExtractPed.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscPulserLinearityCorrection.txt
create mode 100755 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscWritePedCool.py
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.cxx
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.h
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.cxx
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.h
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.cxx
create mode 100644 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.h
create mode 100755 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_entries.cxx
create mode 100755 MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_load.cxx
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/cdb.log b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/cdb.log
new file mode 100644
index 0000000000000..e69de29bb2d1d
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/requirements b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/requirements
new file mode 100755
index 0000000000000..5103ea96159d9
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/cmt/requirements
@@ -0,0 +1,65 @@
+package CscCalibAlgs
+########
+########
+
+private
+use AtlasPolicy AtlasPolicy-*
+
+use AtlasROOT AtlasROOT-* External
+
+use GaudiInterface GaudiInterface-* External
+#use AtlasCORAL AtlasCORAL-* External
+#use AtlasCOOL AtlasCOOL-* External
+#Not used anymore: use AtlasSEAL AtlasSEAL-* External
+
+#use CondDBObjects CondDBObjects-* Database
+#use AthenaPoolUtilities AthenaPoolUtilities-* Database/AthenaPOOL
+
+#use Identifier Identifier-* DetectorDescription
+use MuonIdHelpers MuonIdHelpers-* MuonSpectrometer
+use MuonRDO MuonRDO-* MuonSpectrometer
+#use MuonGeoModel MuonGeoModel-* MuonSpectrometer
+#use MuonByteStream MuonByteStream-* MuonSpectrometer/MuonCnv
+use MuonReadoutGeometry MuonReadoutGeometry-* MuonSpectrometer/MuonDetDescr
+
+
+use StoreGate StoreGate-* Control
+use DataModel DataModel-* Control
+#use CLIDSvc CLIDSvc-* Control
+
+use CscCalibTools CscCalibTools-* MuonSpectrometer/MuonCalib/CscCalib
+use CscCalibData CscCalibData-* MuonSpectrometer/MuonCalib/CscCalib
+
+use MuonCSC_CnvTools MuonCSC_CnvTools-* MuonSpectrometer/MuonCnv
+
+#use MuonCondSvc MuonCondSvc-* MuonSpectrometer/MuonConditions/MuonCondGeneral
+use MuonCondInterface MuonCondInterface-* MuonSpectrometer/MuonConditions/MuonCondGeneral
+#use MuonCondData MuonCondData-* MuonSpectrometer/MuonConditions/MuonCondGeneral
+
+
+end_private
+
+#Monitoring
+#use MuonCalibMonitoring MuonCalibMonitoring-* MuonSpectrometer/MuonValidation/MuonDQA
+
+# Test infrastructure
+use TestPolicy TestPolicy-*
+use MinimalRunTime MinimalRunTime-* Control -no_auto_imports
+
+apply_pattern component_library
+library CscCalibAlgs *.cxx components/*.cxx
+#apply_tag ROOTBasicLibs
+##apply_tag ROOTMathLibs
+#apply_tag ROOTSTLDictLibs
+#apply_tag ROOTGraphicsLibs
+#apply_tag ROOTTableLibs
+
+#apply_pattern dual_use_library files = *.cxx
+
+apply_pattern declare_joboptions files="*.py"
+apply_pattern declare_python_modules files="*.py"
+
+#To compile as debug separately from all else
+private
+#macro cppdebugflags '$(cppdebugflags_s)'
+#macro_remove componentshr_linkopts "-Wl, -s"
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRun.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRun.py
new file mode 100644
index 0000000000000..539dbb52d79ac
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRun.py
@@ -0,0 +1,203 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+from glob import glob
+import subprocess
+import os
+import pickle
+import sys
+
+class BulkRun :
+ """Run on multiple pedestal files
+
+ Used to run on multiple pedestal bytestream files located somewhere on disk.
+ Puts results into seperate directories by run-number. It extracts pertinent
+ information from the file names, which are expected to be of the form:
+
+ data10_calib.00157081.calibration_pedCSC.daq.RAW._lb0000._CSC-EB._0001.data
+
+ run(numToRun) - process NumToRun (at most) runs out of those not yet run out of. This is the primary interface.
+ ReadProcessedFilesList() - find all files already run on from a local text file
+ AddProcessedFiles() - save newly run on files to disk
+ FindFiles() - Find pattern and run number of a set of files that have not yet been run on
+ RunAthena() - run athena job for a set of runs
+
+ """
+ def __init__(
+ self,
+ inputPattern = "/raid02/schernau/ped/csc/csc/*.data" ,
+ processedFilesList = "ProcessedFiles.list",
+ outputDirBase = "/raid02/lampen/datasets/csc/PedProcessing2",
+ debug = False,
+ allowDirOverwrite = False,
+ ):
+ """initialize internal variables"""
+ self.InputPattern = inputPattern
+ self.ProcessedFilesList = processedFilesList
+ self.OutputDirBase = outputDirBase
+ self.debug = debug
+ self.AllowDirOverwrite = allowDirOverwrite
+ print self.InputPattern
+
+ def run(self, numToRun = 10) :
+ """
+ Run over all run numbers, find files for each, and submit each set to
+ CscCalcPedMon.py
+ """
+
+ print "Running on " + str(numToRun) + " runsSet."
+
+
+ runNumbers = []
+ for runCnt in range(numToRun) :
+ print ">>>>Running on runSet " + str(runCnt+1) + " of " + str(numToRun)
+ pattern,runNumber = self.FindFiles()
+ if(pattern != ""):
+ #have files to run on
+ runNumbers += [runNumber]
+ try:
+ self.RunAthena(pattern,runNumber)
+ except:
+ "Failure during RunAthena!"
+ raise
+ else:
+ print "No more unprocessed files. Congrats!"
+ print "N runs done: " + str(runCnt +1)
+ print runNumbers
+ return
+ print "finished all " + str(numToRun)
+ print "Run numbers include:"
+ print runNumbers
+ print
+ print "All Processed files:"
+ print self.ReadProcessedFilesList()
+
+ #Read list of previously processed files
+ def ReadProcessedFilesList(self):
+
+ ProcessedFiles = []
+
+ try:
+ #Get processed files
+ f = open(self.ProcessedFilesList,"r")
+ ProcessedFiles = pickle.load(f)
+ f.close()
+ except:
+ print "No processed file list yet..."
+
+ #Remove newline character from each filename
+
+ #for index in range(len(ProcessedFiles)) :
+ # file = ProcessedFiles[index]
+ # ProcessedFiles[index] = file[0:len(file)-1]
+
+ print 'Processed String: '
+ print ProcessedFiles
+
+ return ProcessedFiles
+
+ def AddProcessedFiles(self,newFiles):
+ """Save new processed files to disk"""
+ ProcessedFiles = self.ReadProcessedFilesList()
+ ProcessedFiles += newFiles
+ f = open(self.ProcessedFilesList,"w")
+ pickle.dump(ProcessedFiles,f)
+ f.close()
+
+ return True
+
+
+ #Read list of previously processed files
+ def FindFiles(self):
+
+ #Initial setup
+ FoundUnprocessedFile = False
+
+ #Get processed file list
+ ProcessedFiles = self.ReadProcessedFilesList()
+
+ #Get list of files in input dir
+ print "Input pattern: " + self.InputPattern
+ inputFiles = glob(self.InputPattern)
+
+ if(self.debug):
+ print
+ print "Searching for file"
+ print "InputPattern: " + self.InputPattern
+ #print "inputFiles: "
+ #print inputFiles
+
+ pattern = ""
+ runNumber = ""
+
+ #Loop through list until find file that is
+ for file in inputFiles:
+ if not ProcessedFiles.count(file):
+
+ index = file.find("data10")
+ if(index == -1):
+ index = file.find("data11")
+ if(index == -1):
+ print "ERROR! Index of -1!"
+ raise Exception("Index error")
+ FoundUnprocessedFile = True
+ pattern = file[0:index + 22] + "*" #includes run number
+ runNumber = file[index + 13: index + 21]
+
+ if(not FoundUnprocessedFile):
+ return "",""
+
+ if(self.debug) :
+ print "Found unprocessed file with pattern: " + pattern
+ print "This includes files:"
+ print glob(pattern)
+
+ return pattern, runNumber
+
+ def RunAthena(self, pattern, runNumber):
+ """Run athena on a particular set of files matching pattern"""
+ outputDirPath = self.OutputDirBase + "/" + runNumber
+
+ if(self.AllowDirOverwrite):
+ try:
+ subprocess.call("rm -rf " + outputDirPath)
+ except:
+ pass
+
+ print "Making directory" + outputDirPath
+ #Create output directory
+ os.mkdir(outputDirPath,0755)
+
+ #Change directory to output directory
+ #os.chdir(outputDirPath)
+
+ #Athena options
+ athOpt = "outputPre='" + outputDirPath +"/" + runNumber
+ athOpt += "';inputPat='" + pattern
+ athOpt += "';reportPrefix='runNumber = "
+ athOpt += runNumber + "'"
+
+ #athena arguments
+ athArgs = ["athena.py", "-c", athOpt, "CscCalcPedMon.py"]
+
+ #Output log file
+ logFile = open(outputDirPath + "/run.log","w")
+
+ try:
+ print
+ print"**************************************"
+ print"Starting running on run " + str(runNumber)
+ sys.stdout.flush()
+ subprocess.Popen(athArgs,stdout=logFile,stderr=subprocess.STDOUT).wait()
+ print"Finished run " + str(runNumber)
+ print"**************************************"
+ print
+ except:
+ print "Error while running athena!"
+ raise
+
+ logFile.close()
+
+
+ #Add files we just ran on to file list
+ newProcessedFiles = glob(pattern);
+ self.AddProcessedFiles(newProcessedFiles)
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRunFollowup.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRunFollowup.py
new file mode 100644
index 0000000000000..89e91fa511bcf
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/BulkRunFollowup.py
@@ -0,0 +1,221 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+from glob import glob
+import subprocess
+import os
+import pickle
+import sys
+
+class BulkRunFollowup :
+ """Run on multiple pedestal files
+
+ #Used to run on multiple pedestal bytestream files located somewhere on disk.
+ #Puts results into seperate directories by run-number. It extracts pertinent
+ #information from the file names, which are expected to be of the form:
+
+ Quickie folowup when BulkRun failed on a few. This will re-enter existing directories
+ of particular run numbers. It should be given a sepearte processedFiles.list
+
+ data10_calib.00157081.calibration_pedCSC.daq.RAW._lb0000._CSC-EB._0001.data
+
+ run() - run over all runs
+ ReadProcessedFilesList() - find all files already run on from a local text file
+ AddProcessedFiles() - save newly run on files to disk
+ FindFiles() - Find pattern and run number of a set of files that have not yet been run on
+ RunAthena() - run athena job for a set of runs
+
+ """
+ def __init__(
+ self,
+ allowedRuns,
+ inputPattern = "/raid02/schernau/ped/csc/*.data",
+ processedFilesList = "ProcessedFiles_Followup.list",
+ outputDirBase = "/raid02/lampen/datasets/csc/PedProcessing",
+ debug = False,
+ allowDirOverwrite = False,
+ ):
+ """initialize internal variables"""
+ self.AllowedRuns = allowedRuns
+ self.InputPattern = inputPattern
+ self.ProcessedFilesList = processedFilesList
+ self.OutputDirBase = outputDirBase
+ self.debug = debug
+ self.AllowDirOverwrite = allowDirOverwrite
+
+ def run(self, numToRun = 10) :
+ """
+ Run over all run numbers, find files for each, and submit each set to
+ CscCalcPedMon.py
+ """
+
+ print "Running on " + str(numToRun) + " runsSet."
+
+
+ runNumbers = []
+ for runCnt in range(numToRun) :
+ print ">>>>Running on runSet " + str(runCnt+1) + " of " + str(numToRun)
+ pattern,runNumber = self.FindFiles()
+ if(pattern != ""):
+ #have files to run on
+ if runNumber in self.AllowedRuns:
+ runNumbers += [runNumber]
+ try:
+ self.RunAthena(pattern,runNumber)
+ except:
+ "Failure during RunAthena!"
+ raise
+ else:
+ if(self.debug):
+ print "Skipping file from runNumber " + str(runNumber) + " (Not allowed)"
+ #Add files we're skipping into file list
+ newProcessedFiles = glob(pattern);
+ self.AddProcessedFiles(newProcessedFiles)
+
+ else:
+ print "No more unprocessed files. Congrats!"
+ print "N runs done: " + str(runCnt +1)
+ print runNumbers
+ return
+ print "finished all " + str(numToRun)
+ print "Run numbers include:"
+ print runNumbers
+ print
+ print "All Processed files:"
+ print self.ReadProcessedFilesList()
+
+ #Read list of previously processed files
+ def ReadProcessedFilesList(self):
+
+ ProcessedFiles = []
+
+ try:
+ #Get processed files
+ f = open(self.ProcessedFilesList,"r")
+ ProcessedFiles = pickle.load(f)
+ f.close()
+ except:
+ print "No processed file list yet..."
+
+ #Remove newline character from each filename
+
+ #for index in range(len(ProcessedFiles)) :
+ # file = ProcessedFiles[index]
+ # ProcessedFiles[index] = file[0:len(file)-1]
+
+ print 'Processed String: '
+ print ProcessedFiles
+
+ return ProcessedFiles
+
+ def AddProcessedFiles(self,newFiles):
+ """Save new processed files to disk"""
+ ProcessedFiles = self.ReadProcessedFilesList()
+ print "Got old processed files"
+ print ProcessedFiles
+ ProcessedFiles += newFiles
+ print "New version:"
+ print ProcessedFiles
+ f = open(self.ProcessedFilesList,"w")
+ pickle.dump(ProcessedFiles,f)
+ f.close()
+
+ return True
+
+
+ #Read list of previously processed files
+ def FindFiles(self):
+
+ #Initial setup
+ FoundUnprocessedFile = False
+
+ #Get processed file list
+ ProcessedFiles = self.ReadProcessedFilesList()
+
+ #Get list of files in input dir
+ inputFiles = glob(self.InputPattern)
+
+ if(self.debug):
+ print
+ print "Searching for file"
+ print "InputPattern: " + self.InputPattern
+ #print "inputFiles: "
+ #print inputFiles
+
+ pattern = ""
+ runNumber = ""
+
+ #Loop through list until find file that is
+ for file in inputFiles:
+ if not ProcessedFiles.count(file):
+
+ index = file.find("data10")
+ if(index ==-1):
+ index = file.find("data11")
+ if(index == -1):
+ print "ERROR! Index of -1 searching for data1X prefix in filename:"
+ print file
+ raise Exception("Index error")
+ FoundUnprocessedFile = True
+ pattern = file[0:index + 22] + "*" #includes run number
+ runNumber = file[index + 13: index + 21]
+
+ if(not FoundUnprocessedFile):
+ if(self.debug): print "Did not find unprocessed file"
+ return "",""
+
+ if(self.debug) :
+ print "Found unprocessed file with pattern: " + pattern
+ print "This includes files:"
+ print glob(pattern)
+
+ return pattern, runNumber
+
+ def RunAthena(self, pattern, runNumber):
+ """Run athena on a particular set of files matching pattern"""
+ outputDirPath = self.OutputDirBase + "/" + runNumber
+
+ if(self.AllowDirOverwrite):
+ try:
+ subprocess.call("rm -rf " + outputDirPath)
+ except:
+ pass
+
+ print "Making directory" + outputDirPath
+ #Create output directory
+ os.mkdir(outputDirPath,0755)
+
+ #Change directory to output directory
+ #(Causes problems)
+ #os.chdir(outputDirPath)
+
+ #Athena options
+ athOpt = "outputPre='" + outputDirPath +"/" + runNumber
+ athOpt += "';inputPat='" + pattern
+ athOpt += "';reportPrefix='runNumber = "
+ athOpt += runNumber + "'"
+
+ #athena arguments
+ athArgs = ["athena.py", "-c", athOpt, "CscCalcPedMon.py"]
+
+ #Output log file
+ logFile = open(outputDirPath + "/run2.log","w")
+
+ try:
+ print
+ print"**************************************"
+ print"Starting running on run " + str(runNumber)
+ sys.stdout.flush()
+ subprocess.Popen(athArgs,stdout=logFile,stderr=subprocess.STDOUT).wait()
+ print"Finished run " + str(runNumber)
+ print"**************************************"
+ print
+ except:
+ print "Error while running athena!"
+ raise
+
+ logFile.close()
+
+
+ #Add files we just ran on to file list
+ newProcessedFiles = glob(pattern);
+ self.AddProcessedFiles(newProcessedFiles)
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/CscCalibQuery.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/CscCalibQuery.py
new file mode 100755
index 0000000000000..e14f6a452d910
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/CscCalibQuery.py
@@ -0,0 +1,421 @@
+#!/afs/cern.ch/sw/lcg/external/Python/2.5.4p2/slc4_ia32_gcc34/bin/python
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+#/usr/bin/python
+import os
+import sys
+import re
+import glob
+import subprocess
+import pickle
+
+
+#This script looks for new CSC calibration files on castor and runs on them. It:
+#1: Checks to see if a job is currently running by looking for a calibrationRunning### file in the script directory. Quits if anything is running.
+#2: If no job is currently running, sees if there are any new calibration directories that we haven't yet processed. Castor is checked and its contents are checked against a Castor Contents local text file. This file isn't nesscary (there is a redundency later on to prevent duplicate calib jobs), but speeds things up.
+#3: If there are new directories, we check to see if any have at least numFilesToRun files. If they don't, we don't update Castor Contents, but we don't run on them either (its assumed they are being staged and we'll pick them up next time CscCalibQuery is run)
+#4: If we find a valid calibration directory with enough files, we check to see if a directory with calibration output for the new calibration run exists. If it does, we don't run on it again.
+#5: If it is indeed a new calibration directory, a bash script is created and submitted to a batch queue. The bash script:
+# a) Runs athena calibration algorithm
+# b) Produces web page with calibration monitoring plots
+# c) Emails mailing list about state of calibration (any problems)
+# d) Generates mysql database file for merging into COOL
+# e) (Disabled) automatically adds to database
+#...
+
+#minimum number of calibration file in castor to start calibration
+numFilesToRun=3
+responsiblePerson = "lampen@physics.arizona.edu"
+#maillist = 'lampen@physics.arizona.edu schernau@positron.ps.uci.edu elliott@physics.arizona.edu prolay@physics.arizona.edu'
+maillist = responsiblePerson
+
+CoolMergeByDefault = False
+
+#Utility functions################################
+
+
+def updateList(oldListPath,newList):
+ print 'updating file list'
+ #update the old file list
+ outFile = open(oldListPath, 'w')
+ outFile.write(newList)
+ outFile.close()
+
+#runs calibration on input file.
+#It creates a bash script which can setup the job, run it, and do post job
+#processing. The bash script is submitted to lxbatch
+def runCalib(calType, runNumber,workDir,castorCopyCmd):
+ #print 'running calib'
+ #initialize based on calibration type
+ scriptDir = '${HOME}/CSC/run/'
+ if calType == 'pulser':
+ calibScript = scriptDir + 'CscCalcSlopeMon.py'
+ dbScript = scriptDir + 'cscWritePSlopeCool.py'
+ webDir = '${HOME}/www/csc/pulser'
+ runListFile = "pulserRunList.pickle"
+ #extractScript = scriptDir +'CscExtractPulser.py'
+ elif( calType == 'ped'):
+ calibScript = scriptDir + 'CscCalcPedMon.py'
+ dbScript = scriptDir + 'cscWritePedRefCool.py'
+ onlDbFile = scriptDir + 'online.cal'
+ webDir = '${HOME}/www/csc/ped'
+ webPageUrl = 'https://atlas-csc-calib.web.cern.ch/atlas-csc-calib/ped/pedRun_' +runNumber
+ runListFile = "pedRunList.pickle"
+ #extractScript = scriptDir + 'CscExtractPed.py'
+
+ outputDir = workDir +"/CalibResults"
+ bsDir = workDir +"/Bytestream"
+
+ print 'outputDir = ' + outputDir
+
+
+ #Setup finished email message
+ emailMessage = 'Finished ' + calType + ' calibration on run number ' + runNumber
+ emailMessage += '. Output in ' + outputDir + '.'
+ emailMessage += '\\nAnd website at:\\n'
+ emailMessage += webPageUrl
+
+ #Setup finished email subjects
+ goodEmailSubject = '[CSC CALIB PROC]: SUCCESS with ' + calType + 'calib run' + runNumber
+ badEmailSubject = '[CSC CALIB PROC]: PROBLEMS with ' + calType + 'calib run' + runNumber
+
+ #Prepare bash script for batch system
+ bashFilePath = workDir+'/CscCalib_' + calType + '_' + runNumber + '.sh'
+
+ bsubCmd ='cd ' + outputDir + ';bsub -q 2nd -R "type==SLC5_64&&mem>420" ' + bashFilePath
+
+ bashFileContents = "#!/bin/bash\n"
+ bashFileContents += "#To resubmit this job, submit it to the atlasmuonqueu like so:\n#"
+ bashFileContents += bsubCmd +"\n"
+ bashFileContents += "source ~/CSC/CscSetup.sh\n"
+ bashFileContents += "\n"
+ bashFileContents += "resultDir=\"" + outputDir + "\"\n"
+ bashFileContents += 'bytestreamDir="' + bsDir + '"\n'
+ bashFileContents += 'maillist="' + maillist + '"\n'
+ bashFileContents += 'webSiteDir="' + webDir + '"\n'
+
+ calFilePrefix= "${resultDir}/" + runNumber
+ inputPattern = "${bytestreamDir}/*.data"
+
+ calibCommand = 'echo "Running calibration"\n' \
+ + 'mkdir ${resultDir}\n' \
+ + 'athena.py -c "outputPre=\'' + calFilePrefix \
+ + '\';inputOnlCalibFile=\'' +onlDbFile \
+ + '\';inputPat=\'' + inputPattern \
+ + '\';reportPrefix=\'' + emailMessage \
+ + '\';" ' \
+ + calibScript
+
+
+ goodEmailCommand = ' mail -s "' + goodEmailSubject + '" $maillist < ' + calFilePrefix + "CalibReport.txt"
+ badEmailCommand = ' mail -s "' + badEmailSubject + '" $maillist < ' + calFilePrefix + "CalibReport.txt"
+
+ #For reference tag, we actually want the IOV to start just after the LAST run number
+ #Get old run numbers
+ infile = open(runListFile,"r")
+ runList = pickle.load(infile)
+ runList.sort()
+ print "got runs"
+ print runList
+ infile.close()
+
+
+ if(runNumber in runList):
+ print "Mailing message"
+ message =["mail","-s",\
+ '"New castor run directory found for previously processed run ' + str(runNumber) + '"',\
+ responsiblePerson,\
+ "<",\
+ "runAlreadyProcessed.mail"]
+ print message
+ subprocess.call(message)
+ sys.exit()
+
+
+ highestRun = runList[-1]
+
+ isRunNumberConflict = False
+
+
+ if(highestRun > runNumber):
+ #Something odd happening. The Input run number is lower than the last run number
+ #this script (thinks it) processed
+ #Notify someone important, and don't add to cool when done...
+ subprocess.call(["mail","-s",\
+ "[CSC CALIB PROC] Wrong run number ordering on run " + str(runNumber) \
+ + "! Human intervension required!",\
+ responsiblePerson,\
+ "<",\
+ "runNumberConflict.mail"]\
+ )
+ isRunNumberConflict = True
+ else:
+ #No problem, update run list
+ runList += [runNumber]
+ outfile = open(runListFile,"w")
+ pickle.dump(runList,outfile)
+ outfile.close()
+
+
+ #Label that we're working, so that other instances of CscCalibQuery won't run
+ subprocess.call(['touch','/afs/cern.ch/user/m/muoncali/CSC/run/runningCalibration' + runNumber])
+
+ #Command to create .db file
+ DbCommand = 'athena.py -c "input=\'' + calFilePrefix+'.cal\';output=\'' \
+ + calFilePrefix + '.db\';IOVRunStart=int(\'' + highestRun + '\')" ' + dbScript
+
+ UploadCommand = ''
+ #Noexec prevents actual execution of cool
+ #UploadCommand = '/afs/cern.ch/user/a/atlcond/utils/AtlCoolMerge.py ' \
+ # + "--batch --noexec --comment=\'Automated UPD1 update from " + calType \
+ # + ' run ' + runNumber +'\' ' \
+ # + calFilePrefix + '.db' \
+ # + ' COMP200 ATONR_COOL ATLAS_COOLONL_CSC_W WCOOLONL4CSC17 '
+ #UploadCommand += "\n"
+
+
+
+
+
+ #Command to upload .db file to database
+ UploadCommand += '/afs/cern.ch/user/a/atlcond/utils/AtlCoolMerge.py ' \
+ + "--batch "
+ if(not CoolMergeByDefault or isRunNumberConflict):
+ UploadCommand += " --noexec "
+
+ UploadCommand += "--comment=\'Automated reference update from " + calType \
+ + ' run ' + runNumber + ' to IOV starting at ' + highestRun + "' " \
+ + calFilePrefix + '.db' \
+ + ' COMP200 ATLAS_COOLWRITE ATLAS_COOLOFL_CSC_W WCOOLOFL4CSC17 '
+
+ WebSiteCommand = '\nfs sa ${resultDir} webserver:afs read\n'
+ WebSiteCommand += 'cd $resultDir\n'
+ WebSiteCommand += 'ln -s ${resultDir} ${webSiteDir}/pedRun_' + runNumber + '\n'
+ WebSiteCommand += 'MakeCscPedSite.exe ' + runNumber + '.root ' + runNumber + '.cal_online\n'
+
+ t1 = '\t'
+
+ #Run calibration. If no problems detected, go ahead and upload.
+ bashFileContents += '\ncd ' + workDir + '\n' \
+ + "#Copying files from castor#\n" \
+ + castorCopyCmd \
+ + '\n'\
+ + "#Running calibration (and calib monitoring)\n"\
+ + "#, resulting in .cal file and status report\n"\
+ + calibCommand +'\n' \
+ + "\n"\
+ + "#Athena job to transform *.cal file to *.db file.\n"\
+ + DbCommand + '\n' \
+ + "#Python utility to upload *.db file to database. When entering by\n"\
+ + "#hand, I recomend removing '--batch' flag so you can check puposed\n"\
+ + "#operations before submision to database\n"\
+ + UploadCommand + '\n' \
+ + '#Check if AllCalibMonGood file was created, which means that\n'\
+ + '#this run passed acceptable criteria in the calibration monitoring\n'\
+ + 'if [ -a AllCalibMonGood ]; then\n' \
+ + t1 + "#Email list that the calibration looks good\n"\
+ + t1 + goodEmailCommand +'\n' \
+ + t1 + "################################################################\n"\
+ + t1 + "#Execute next two commands if you want to submit database entry#\n"\
+ + t1 + "#Useful if these steps were skipped due to suspicious behaviour#\n"\
+ + t1 + "#during calibration. #\n"\
+ + t1 + "###############################################################\n"\
+ + t1 + "\n"\
+ + t1 + '\n'\
+ + 'else\n' \
+ + t1 + "#Suspicious behaviour in calibration. Notify mail list of this fact\n"\
+ + t1 + badEmailCommand + '\n' \
+ + 'fi\n'\
+ + '\n'\
+ + '#Always create website'\
+ + WebSiteCommand + '\n' \
+ + 'rm -rf $bytestreamDir\n' \
+ + 'rm -rf ' + scriptDir +"runningCalibration" + runNumber +'\n'
+
+ #Write bashfile
+ print "Printing bash file to: " +bashFilePath
+ bashFile = open(bashFilePath, 'w')
+ bashFile.write(bashFileContents)
+ bashFile.close()
+
+ #Submit script
+ os.system('chmod +x ' + bashFilePath)
+ bsubsMessage = os.popen(bsubCmd).read()
+
+ #Send alert email
+ emailMessage = 'Starting ' + calType + ' calibration on run number ' + runNumber
+ emailSubject = '[CSC CALIB PROC]: ' + emailMessage
+ emailMessage += '\nbsubs output:\n' + bsubsMessage
+ os.system('echo "' + emailMessage + '" | mail -s "' + emailSubject + '" ' + maillist)
+
+#########################################
+#Main program#############################
+#########################################
+
+print 'running'
+#first command line argument should be the calibration type, pulser or ped.
+try:
+ calType = sys.argv[1]
+ if calType == 'pulser':
+ calibFileDir = '/castor/cern.ch/user/l/lampen/CalibRunTest/slope/'
+ oldListFilePath = '/afs/cern.ch/user/m/muoncali/CSC/run/pulserList.txt'
+ outputDir = '/afs/cern.ch/user/m/muoncali/w0/CSC/runs/pulser/pulser'
+ elif calType == 'ped':
+ calibFileDir = '/castor/cern.ch/grid/atlas/DAQ/muon/csc/'
+ oldListFilePath = '/afs/cern.ch/user/m/muoncali/CSC/run/pedDirList.txt'
+
+ #oldListFilePath = '/afs/cern.ch/user/l/lampen/Calib/dev/rel_0/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/pedList.txt'
+ #outputDir = '/afs/cern.ch/user/l/lampen/Calib/dev/rel_0/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/python/'
+ outputDir = '/afs/cern.ch/user/m/muoncali/w0/CSC/runs/ped/ped'
+ #data11_calib.00172807.calibration_pedCSC.daq.RAW._lb0000._CSC-EB._0001.data
+ #calibRe = re.compile('data10_calib.*calibration_ped\.daq\.RAQ.*\.data')
+ #calibRe = re.compile('data10_.*pedCSC.*\.data')
+ calibRe = re.compile('data1.*_calib.*calibration_pedCSC\.daq\.RAW.*\.data')
+ else:
+ print 'Need to specify pulser or ped'
+ os._exit(0)
+except:
+ print 'Need to specify pulser or ped'
+ os._exit(0)
+
+
+#First, see if a calibration is already running
+#If already running, stop
+testFile = glob.glob("/afs/cern.ch/user/m/muoncali/CSC/run/runningCalibration*")
+if(testFile):
+ print 'already running: ' + str(testFile)
+ sys.exit()
+
+
+#Get the current list of directories in castor
+print 'rfdir ' + calibFileDir
+currentLs = os.popen('rfdir ' + calibFileDir).read()
+
+os.popen('touch testing')
+
+#Get the old list of files from castor
+inFile = open(oldListFilePath, 'r')
+oldLs = inFile.read()
+
+inFile.close()
+
+print 'checking for changes'
+
+import datetime
+now = datetime.datetime.now()
+today = now.day
+
+#Check if there are any changes between the lists
+if(oldLs != currentLs):
+ print 'There has been a change'
+ currentDirList = currentLs.split('\n')
+ oldDirList = oldLs.split('\n')
+
+ updateRunList = True
+ runningDir = ""
+ #Run on any new directories
+ for Dir in currentDirList:
+ #print 'Checking ' + Dir
+ if Dir not in oldDirList:
+ print "**********************************************************"
+ splitDir = Dir.split()
+ day = int(splitDir[-3])
+ DirName =splitDir[-1] #last ' ' delimited word in line is Dirname
+
+ print splitDir
+ timediff = today - day
+ if(timediff > 7 or timediff < -23):
+ print timediff
+
+ print "Found old dir " + DirName + ", but its over at least a week old, so we're skipping it"
+
+ continue
+
+
+
+ print "day is " + str(day)
+
+ #print 'Dirname is ' + DirName
+ cmd = 'rfdir ' + calibFileDir + DirName
+ fileList = os.popen(cmd).read().split('\n')
+
+ nFiles = len(fileList) -1
+ print "found " + str(nFiles) + " files"
+
+ runNumber = DirName
+
+ #prepare output directory
+ outputDirFull = outputDir + 'Run_' + runNumber
+ print 'outputDirFull is ' + outputDirFull
+
+ #Loop through files in directory.
+ #Start building castor copy cmd
+ #If any files don't match the calibration file
+ #requirement, mark this directory as something not interested in.
+ ThisCalibDir = False
+ castorCopyCmd = "mkdir ${bytestreamDir}"
+ for file in fileList:
+ fileName = (file.split(' '))[-1] #last ' ' delimited word in line is fileName
+ if(fileName != ""):
+ print "fileName: " +fileName
+ ThisCalibDir = re.match(calibRe,fileName)
+ if(nFiles < numFilesToRun):
+ print "only " +str(nFiles) + " files. Breaking."
+ if(ThisCalibDir):
+ print "There is a calib dir, but it only has " + str(nFiles) \
+ + " file. Will not process."
+ updateRunList = False
+ break;
+ if(ThisCalibDir):
+ fullPath = calibFileDir + DirName + '/' + fileName
+ castorCopyCmd += "\nxrdcp root://castoratlas/" + fullPath + " ${bytestreamDir}";
+ #print 'found ' + fullPath
+ else:
+ break
+
+ print "found " + str(nFiles) + " files"
+
+ #only run if have enough files
+ if(nFiles >= numFilesToRun):
+ #Comment next 4 lines out if you want to run on multiple runs at once
+ if(runningDir):
+ updateRunList = False
+ print "Found new calibration directory to run on (" +DirName + "), but already running on " + runningDir
+ continue
+ print 'running on ' + DirName
+
+ if(nFiles > numFilesToRun):
+ print 'WARNING! Found ' + str(nFiles) + ' calib files, when only ' + str(numFilesToRun) + ' were expected!'
+
+
+ if( os.path.exists(outputDirFull)):
+ print "There is already a directory for run number " + runNumber
+ print "Will not run on this pedestal run"
+ continue
+
+
+ if(ThisCalibDir):
+ print castorCopyCmd
+ #updateList(oldListFilePath,currentLs)
+ #os.makedirs(outputDirFull+"/Bytestream")
+ os.makedirs(outputDirFull+"/CalibResults")
+
+ runCalib(calType, runNumber, outputDirFull, castorCopyCmd)
+
+ print('Launched job, but will continue to check if we have more runs to run on later')
+ runningDir = DirName
+ updateRunList = False
+ continue
+
+
+ if( updateRunList):
+ #We update the run list so long as there are no runs on it that we expect to run on in the future.
+ #The safest thing to do is to wait until all runs are processed.
+ print "No unprocessed or currently being processed calibration runs. Will update run list."
+ updateList(oldListFilePath,currentLs)
+ else:
+ print "NOT updating run list"
+
+else:
+ print 'No changes between old and new runs'
+ pass
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcPedMon.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcPedMon.py
new file mode 100644
index 0000000000000..d51b03feaeeb7
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcPedMon.py
@@ -0,0 +1,227 @@
+#**************************************************************
+#
+# Csc Pedestal Calibration
+#
+#==============================================================
+
+runNumber = 2147483647
+#Input files:
+#prelimDir = "/raid01/lampen/testarea/Calib/August09CosmicRun/data/preliminaryData/"
+
+#AugustRunDir = "rfio:/castor/cern.ch/grid/atlas/DAQ/muon/csc/00124401/"
+#myInputFiles = [
+# AugustRunDir + 'data_test.00124401.calibration_ped.daq.RAW._lb0000._CSC-EB._0001.data',
+# AugustRunDir + 'data_test.00124401.calibration_ped.daq.RAW._lb0000._CSC-EB._0002.data',
+# AugustRunDir + 'data_test.00124401.calibration_ped.daq.RAW._lb0000._CSC-EB._0003.data',
+# ]
+
+#SeptRunDir = 'rfio:/castor/cern.ch/grid/atlas/DAQ/muon/csc/00129083/'
+#myInputFiles = [
+# SeptRunDir + 'data_test.00129083.calibration_ped.daq.RAW._lb0000._CSC-EB._0001.data',
+# SeptRunDir + 'data_test.00129083.calibration_ped.daq.RAW._lb0000._CSC-EB._0002.data',
+# SeptRunDir + 'data_test.00129083.calibration_ped.daq.RAW._lb0000._CSC-EB._0003.data',
+# ]
+
+#OctRunDir = '/raid02/lampen/datasets/csc/pedRuns/oct/raw/'
+#OctRUn1
+#myInputFiles = [
+# OctRunDir + 'data_test.00132435.calibration_ped.daq.RAW._lb0000._CSC-EB._0001.data',
+# OctRunDir + 'data_test.00132435.calibration_ped.daq.RAW._lb0000._CSC-EB._0002.data',
+# OctRunDir + 'data_test.00132435.calibration_ped.daq.RAW._lb0000._CSC-EB._0003.data'
+# ]
+
+#OctRun2
+#myInputFiles = [
+# OctRunDir + 'data_test.00133116.calibration_ped.daq.RAW._lb0000._CSC-EB._0001.data',
+#OctRunDir + 'data_test.00133116.calibration_ped.daq.RAW._lb0000._CSC-EB._0002.data',
+#OctRunDir + 'data_test.00133116.calibration_ped.daq.RAW._lb0000._CSC-EB._0003.data'
+#]
+NovRunDir="/raid02/lampen/datasets/csc/pedRuns/nov/raw/*.data"
+
+import glob
+
+if('inputPat' not in dir()):
+ inputPat = NovRunDir
+
+myInputFiles = glob.glob(inputPat)
+#ServiceMgr.EventSelector.InputCollections=ListFromCastor(CastorDir)
+
+
+if('outputPre' not in dir()):
+ outputPre = "pedestal"
+
+if('outputFile' not in dir()):
+ outputFile = outputPre + '.root'
+
+if( 'inputOnlCalibFile' not in dir()):
+ inputOnlCalibFile = ""
+
+#myInputFiles = [
+# prelimDir + 'data_test.00124166.calibration_ped.daq.RAW._lb0000._CSC-EB._0002.data'
+#]
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+
+from AthenaCommon.AppMgr import ServiceMgr
+from AthenaCommon.AppMgr import ToolSvc
+
+#---------------------------------------------------------------
+# Detector Initialization #---------------------------------------------------------------
+#Set to read from real data COOL data
+
+from AthenaCommon.GlobalFlags import GlobalFlags
+GlobalFlags.DetGeo.set_atlas()
+GlobalFlags.DataSource.set_data()
+from IOVDbSvc.CondDB import conddb
+
+DetDescrVersion = "ATLAS-GEONF-08-00-00"
+include ("RecExCond/RecExCommon_flags.py")
+DetFlags.ID_setOff()
+DetFlags.Calo_setOff()
+include ("RecExCond/AllDet_detDescr.py")
+
+#Set your conditions global tag
+conddb.setGlobalTag('COMCOND-ES1P-002-00')
+
+#August 09 cosmics
+#conddb.setGlobalTag('COMCOND-MONC-002-00')
+
+#Sept 08 cosmics
+#conddb.setGlobalTag('COMCOND-REPC-002-02')
+
+conddb.addOverride("/CSC","Csc-REF-001-00")
+
+Service("IOVDbSvc").forceRunNumber = runNumber;
+Service("IOVDbSvc").forceLumiblockNumber = 1
+
+#---------------------------------------------------------------
+# Event Loop
+#---------------------------------------------------------------
+include( "ByteStreamCnvSvc/BSEventStorageEventSelector_jobOptions.py" )
+
+
+
+#---------------------------------------------------------------
+# CSC ROD Decoder
+#---------------------------------------------------------------
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import Muon__CscROD_Decoder
+CscRodDecoder = Muon__CscROD_Decoder( name = "CscROD_Decoder",
+ IsCosmics = False, #False for collision
+ #IsCosmics = True, #False for collision
+ IsOldCosmics = False)
+CscRodDecoder.OutputLevel = WARNING
+ToolSvc += CscRodDecoder
+
+#---------------------------------------------------------------
+# CSC BS -> RDO
+#---------------------------------------------------------------
+# --- ROD Decoder
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import Muon__CSC_RawDataProviderTool
+MuonCscRawDataProviderTool = Muon__CSC_RawDataProviderTool(name = "MuonCscRawDataProviderTool",
+ Decoder = CscRodDecoder)
+ToolSvc += MuonCscRawDataProviderTool
+print MuonCscRawDataProviderTool
+
+# --- RawData Provider
+from MuonByteStream.MuonByteStreamConf import Muon__CscRawDataProvider
+topSequence += Muon__CscRawDataProvider(name = "MuonCscRawDataProvider",
+ ProviderTool = ToolSvc.MuonCscRawDataProviderTool)
+print topSequence.MuonCscRawDataProvider
+
+# --- BS Converter
+theApp.Dlls += [ "ByteStreamCnvSvc", "ByteStreamCnvSvcBase"]
+
+# --- Services
+theApp.ExtSvc += ["ByteStreamEventStorageInputSvc/ByteStreamInputSvc"]
+theApp.ExtSvc += [ "ROBDataProviderSvc/ROBDataProviderSvc" ]
+theApp.ExtSvc += [ "EventSelectorByteStream/EventSelector"]
+theApp.EvtSel = "EventSelector"
+
+#---------------------------------------------------------------
+# Muon Converters
+#---------------------------------------------------------------
+include ( "MuonEventCnvTools/MuonEventCnvTools_jobOptions.py" )
+include ( "MuonEventAthenaPool/MuonEventAthenaPool_joboptions.py" )
+#include ("MuonRdoToPrepData/MuonRdoToMuonPrepData_jobOptions.py")
+
+
+#---------------------------------------------------------------
+# BS input files
+#---------------------------------------------------------------
+
+print 'File list is'
+print myInputFiles
+
+ServiceMgr.ByteStreamInputSvc.FullFileName = myInputFiles
+#ServiceMgr.ByteStreamInputSvc.FullFileName = myTestFiles
+
+#---------------------------------------------------------------
+# CSC Calibration from DB
+#---------------------------------------------------------------
+include ("CscCalibTools/CscCalibTool_jobOptions.py")
+
+theApp.EvtMax = -1
+
+#load main algorithms
+from CscCalibAlgs.CscCalibAlgsConf import MuonCalib__CscCalcPed
+CscCalcPed = MuonCalib__CscCalcPed( "CscCalcPed" )
+topSequence += CscCalcPed
+
+#Specify output calibration file
+CscCalcPed.OutputFile = outputPre + ".cal"
+#CscCalcPed.CalOutputVersion = "00-00" #old hash based
+CscCalcPed.CalOutputVersion = "03-00" #New, sca address based
+CscCalcPed.OnlineCalibFile = inputOnlCalibFile
+CscCalcPed.DoBitHists = False
+CscCalcPed.DoSampHists = False
+CscCalcPed.DoCorrelation = False;
+#CscCalcPed.OutputLevel =VERBOSE
+#Define calculations to do. PedestalRun is contradictory to the other options.
+
+from AthenaMonitoring.AthenaMonitoringConf import AthenaMonManager
+topSequence += AthenaMonManager( "MuonMonManager" )
+
+# ******************************
+# Muon Monitoring configuration
+# ******************************
+monMan = topSequence.MuonMonManager
+#monMan = topSequence.PrimaryManager
+monMan.FileKey = "CscCalibMon" #"muoncosmics" ## FileKey must match that given to THistSvc
+
+## Set global monitoring parameters: see the AthenaMonManager class
+# in the Control/AthenaMonitoring package
+monMan.ManualDataTypeSetup = True
+monMan.DataType = "monteCarlo"
+monMan.Environment = "user"
+monMan.ManualRunLBSetup = True
+monMan.Run = 1
+monMan.LumiBlock = 1
+
+from MuonCalibMonitoring.MuonCalibMonitoringConf import CscCalibMonToolPed
+CscCalibMonTool = CscCalibMonToolPed ("CscCalibMonToolPed" )
+ToolSvc += CscCalibMonTool
+monMan.AthenaMonTools += [CscCalibMonTool]
+CscCalibMonTool.GetDebugForAllChannels = False
+CscCalibMonTool.StatusReportFileName = outputPre + "CalibReport.txt"
+if('reportPrefix' in dir()):
+ CscCalibMonTool.StatusReportPrefix = reportPrefix
+
+#CscCalibMonTool.OutputLevel = VERBOSE
+
+#Ensure what we need from cond database is loaded
+from MuonCondSvc.CscCondDB import cscCondDB
+cscCondDB.addPedFolders() #Set to read pedestal and noise folders from the database
+cscCondDB.addRmsFolder()
+cscCondDB.addF001Folder()
+
+#--------------------------------------------------------------
+# Output files.
+#--------------------------------------------------------------
+
+from GaudiSvc.GaudiSvcConf import THistSvc
+ServiceMgr += THistSvc()
+ServiceMgr.THistSvc.Output = ["CscCalibMon DATAFILE='" + outputFile + "' OPT='RECREATE'"]
+
+ServiceMgr.MessageSvc.infoLimit = 0
+ServiceMgr.MessageSvc.debugLimit = 0
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcSlopeMon.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcSlopeMon.py
new file mode 100644
index 0000000000000..4fcc24014f740
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalcSlopeMon.py
@@ -0,0 +1,198 @@
+runNumber = 2147483647
+
+if('inputPat' not in dir()):
+ inputPat = NovRunDir
+
+import glob
+myInputFiles = glob.glob(inputPat)
+
+#Lots of pulses
+#defInputFile = "/afs/cern.ch/user/l/lampen/public/bytestreamTest/daq.cal.0079571.No.Streaming.LB0000.CSC-EB._0001.data"
+#myGeometryTag = 'ATLAS-GEO-01-00-00'
+#myIsCosmics = False
+#Old Cosmics
+#defInputFile = "rfio:/castor/cern.ch/atlas/muon/CSC/daq_CSC-EB-RCD__0001209_file01.data"
+#myGeometryTag = "MuonSpectrometer-CSC-CS-00"
+#myIsCosmics = True
+#Cosmics
+#defInputFile = "rfio:/castor/cern.ch/atlas/muon/CSC/daq_CSC-EB-RCD__0001853_file01.data"
+#myGeometryTag = "MuonSpectrometer-CSC-CS-00"
+#myIsCosmics = True
+############################
+
+#######################################################################
+#Input/Ouptut job settings.
+#Command line options can give an input file as well as an output prefix
+#######################################################################
+if('outputPre' not in dir()):
+ outputPre = "pulser"
+
+if('outputTextFile' not in dir()):
+ outputTextFile = outputPre + ".cal"
+
+if('outputFile' not in dir()):
+ outputFile = outputPre + '.root'
+
+#include('CscCalibAlgs/CscCalibCommonOptions.py')
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+
+from AthenaCommon.AppMgr import ServiceMgr
+from AthenaCommon.AppMgr import ToolSvc
+
+#---------------------------------------------------------------
+# Detector Initialization #---------------------------------------------------------------
+#Set to read from real data COOL data
+
+from AthenaCommon.GlobalFlags import GlobalFlags
+GlobalFlags.DetGeo.set_atlas()
+GlobalFlags.DataSource.set_data()
+from IOVDbSvc.CondDB import conddb
+
+DetDescrVersion = "ATLAS-GEONF-08-00-00"
+include ("RecExCond/RecExCommon_flags.py")
+DetFlags.ID_setOff()
+DetFlags.Calo_setOff()
+include ("RecExCond/AllDet_detDescr.py")
+
+#Setup conditions global tag
+conddb.setGlobalTag('COMCOND-ES1P-002-00')
+
+#Want to view most recent values
+Service("IOVDbSvc").forceRunNumber = runNumber;
+Service("IOVDbSvc").forceLumiblockNumber = 1
+
+#---------------------------------------------------------------
+# Event Loop
+#---------------------------------------------------------------
+include( "ByteStreamCnvSvc/BSEventStorageEventSelector_jobOptions.py" )
+
+#---------------------------------------------------------------
+# CSC ROD Decoder
+#---------------------------------------------------------------
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import Muon__CscROD_Decoder
+CscRodDecoder = Muon__CscROD_Decoder( name = "CscROD_Decoder",
+ IsCosmics = False, #False for collision
+ #IsCosmics = True, #False for collision
+ IsOldCosmics = False)
+CscRodDecoder.OutputLevel = WARNING
+ToolSvc += CscRodDecoder
+
+#---------------------------------------------------------------
+# CSC BS -> RDO
+#---------------------------------------------------------------
+# --- ROD Decoder
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import Muon__CSC_RawDataProviderTool
+MuonCscRawDataProviderTool = Muon__CSC_RawDataProviderTool(name = "MuonCscRawDataProviderTool",
+ Decoder = CscRodDecoder)
+ToolSvc += MuonCscRawDataProviderTool
+print MuonCscRawDataProviderTool
+
+# --- RawData Provider
+from MuonByteStream.MuonByteStreamConf import Muon__CscRawDataProvider
+topSequence += Muon__CscRawDataProvider(name = "MuonCscRawDataProvider",
+ ProviderTool = ToolSvc.MuonCscRawDataProviderTool)
+print topSequence.MuonCscRawDataProvider
+
+# --- BS Converter
+theApp.Dlls += [ "ByteStreamCnvSvc", "ByteStreamCnvSvcBase"]
+
+# --- Services
+theApp.ExtSvc += ["ByteStreamEventStorageInputSvc/ByteStreamInputSvc"]
+theApp.ExtSvc += [ "ROBDataProviderSvc/ROBDataProviderSvc" ]
+theApp.ExtSvc += [ "EventSelectorByteStream/EventSelector"]
+theApp.EvtSel = "EventSelector"
+
+#---------------------------------------------------------------
+# Muon Converters
+#---------------------------------------------------------------
+include ( "MuonEventCnvTools/MuonEventCnvTools_jobOptions.py" )
+include ( "MuonEventAthenaPool/MuonEventAthenaPool_joboptions.py" )
+#include ("MuonRdoToPrepData/MuonRdoToMuonPrepData_jobOptions.py")
+
+
+#Set input files
+ServiceMgr.ByteStreamInputSvc.FullFileName = myInputFiles
+
+#---------------------------------------------------------------
+# CSC Calibration from DB
+#---------------------------------------------------------------
+include ("CscCalibTools/CscCalibTool_jobOptions.py")
+
+theApp.EvtMax = -1
+
+#load main algorithms
+
+from CscCalibAlgs.CscCalibAlgsConf import MuonCalib__CscCalcSlope
+CscCalcSlope = MuonCalib__CscCalcSlope( "CscCalcSlope" )
+topSequence += CscCalcSlope
+#Specify output calibration file
+CscCalcSlope.OutputFile =outputTextFile
+
+#Define calculations to do. PedestalRun is contradictory to the other options.
+CscCalcSlope.DoCrossTalkFix = True
+CscCalcSlope.DoBipolarFit = True
+CscCalcSlope.GetPedFromFile = False
+#CscCalcSlope.OutputLevel = DEBUG
+
+#if(ToolSvc.CscRdoContByteStreamTool.IsCosmicData) :
+# CscCalcSlope.ExpectedChamberLayer = 1
+
+from AthenaMonitoring.AthenaMonitoringConf import AthenaMonManager
+topSequence += AthenaMonManager( "MuonMonManager" )
+
+# ******************************
+# Muon Monitoring configuration
+# ******************************
+monMan = topSequence.MuonMonManager
+#monMan = topSequence.PrimaryManager
+monMan.FileKey = "CscCalibMon" #"muoncosmics" ## FileKey must match that given to THistSvc
+
+## Set global monitoring parameters: see the AthenaMonManager class
+# in the Control/AthenaMonitoring package monMan.ManualDataTypeSetup = True
+monMan.DataType = "monteCarlo"
+monMan.Environment = "user"
+monMan.ManualRunLBSetup = True
+monMan.Run = 1
+monMan.LumiBlock = 1
+
+from MuonCalibMonitoring.MuonCalibMonitoringConf import CscCalibMonToolSlope
+CscCalibMonTool = CscCalibMonToolSlope (
+ name = "CscCalibMonToolSlope"
+ )
+ToolSvc += CscCalibMonTool
+##CscCalibMonTool.OutputLevel = DEBUG
+monMan.AthenaMonTools += [CscCalibMonTool]
+#Conditions DB service
+#from IOVDbSvc.CondDB import conddb
+#conddb.setGlobalTag("COMCOND-003-00")
+from MuonCondSvc.CscCondDB import cscCondDB
+cscCondDB.addPedFolders() #Set to read pedestal and noise folders from the database
+cscCondDB.addPSlopeFolder() #Set to read slope folder from database
+#if(ToolSvc.CscRdoContByteStreamTool.IsCosmicData) :
+# cscCondDB.overrideBaseTag('Csc-cosmic-DEFAULTCOND')
+
+
+svcMgr.MessageSvc.Format = "% F%30W%S%7W%R%T %0W%M"
+svcMgr.MessageSvc.defaultLimit=1000000;
+
+#--------------------------------------------------------------
+# Set output level threshold (1=VERBOSE, 2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL )
+#--------------------------------------------------------------
+
+#---------------------------------------------------------------------------
+print "List all DLL"
+print theApp.Dlls
+print "List all ExtSvc"
+print theApp.ExtSvc
+print "List of all top algorithms"
+print theApp.TopAlg
+#---------------------------------------------------------------------------
+
+
+from GaudiSvc.GaudiSvcConf import THistSvc
+ServiceMgr += THistSvc()
+ServiceMgr.THistSvc.Output = ["CscCalibMon DATAFILE='" + outputFile + "' OPT='RECREATE'"]
+
+
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalibCommonOptions.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalibCommonOptions.py
new file mode 100644
index 0000000000000..433c7fdcc3b29
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscCalibCommonOptions.py
@@ -0,0 +1,119 @@
+#**************************************************************
+#
+# Read in ByteStream events, and Access CSC RDO
+#
+#==============================================================
+
+from AthenaCommon.AppMgr import AlgSequence
+topSequence = AlgSequence()
+
+from AthenaCommon.AppMgr import ServiceMgr
+from AthenaCommon.AppMgr import ToolSvc
+
+from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
+
+from AthenaCommon.GlobalFlags import globalflags
+globalflags.DetDescrVersion.set_Value_and_Lock('ATLAS-GEO-04-00-00')
+globalflags.ConditionsTag.set_Value_and_Lock('COMCOND-ES1C-000-00')
+globalflags.InputFormat.set_Value_and_Lock('bytestream')
+globalflags.DetGeo = 'commis'
+globalflags.DataSource.set_Value_and_Lock('data')
+
+
+from AthenaCommon.BeamFlags import jobproperties
+jobproperties.Beam.beamType.set_Value_and_Lock("cosmics")
+
+from AthenaCommon.BFieldFlags import jobproperties
+jobproperties.BField.solenoidOn.set_Value_and_Lock(True)
+jobproperties.BField.barrelToroidOn.set_Value_and_Lock(True)
+jobproperties.BField.endcapToroidOn.set_Value_and_Lock(True)
+
+
+from MuonRecExample.MuonRecFlags import muonRecFlags
+# Run on CSC data
+muonRecFlags.doCSCs.set_Value_and_Lock(True)
+
+from RecExConfig.RecFlags import rec
+rec.Commissioning.set_Value_and_Lock(True)
+rec.oldRecExCommissionConfig.set_Value_and_Lock(False)
+rec.oldFlagTopSteering.set_Value_and_Lock(False)
+rec.ScopingLevel=4
+
+
+# Detector Initialization
+include ("RecExCond/RecExCommon_flags.py")
+DetFlags.ID_setOff()
+DetFlags.Tile_setOff()
+DetFlags.LAr_setOff()
+DetFlags.Calo_setOff()
+include ("RecExCond/AllDet_detDescr.py")
+
+
+include( "ByteStreamCnvSvc/BSEventStorageEventSelector_jobOptions.py" )
+
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import CscROD_Decoder
+CscRodDecoder = CscROD_Decoder(name = "CscROD_Decoder",
+ IsCosmics = False,
+ IsOldCosmics = False)
+CscRodDecoder.OutputLevel = VERBOSE
+ToolSvc += CscRodDecoder
+
+from MuonCSC_CnvTools.MuonCSC_CnvToolsConf import Muon__CSC_RawDataProviderTool
+MuonCscRawDataProviderTool = Muon__CSC_RawDataProviderTool(name = "MuonCscRawDataProviderTool",
+ Decoder = CscRodDecoder)
+ToolSvc += MuonCscRawDataProviderTool
+print MuonCscRawDataProviderTool
+
+# load the CscRawDataProvider
+from MuonByteStream.MuonByteStreamConf import Muon__CscRawDataProvider
+topSequence += Muon__CscRawDataProvider(name = "MuonCscRawDataProvider",
+ ProviderTool = ToolSvc.MuonCscRawDataProviderTool)
+print topSequence.MuonCscRawDataProvider
+
+
+#from MuonMDT_CnvTools.MuonMDT_CnvToolsConf import MdtROD_Decoder
+#MdtRodDecoder = MdtROD_Decoder(name = "MdtROD_Decoder")
+#ToolSvc += MdtRodDecoder
+
+#from MuonMDT_CnvTools.MuonMDT_CnvToolsConf import Muon__MDT_RawDataProviderTool
+#MuonMdtRawDataProviderTool = Muon__MDT_RawDataProviderTool(name = "MuonMdtRawDataProviderTool",
+# Decoder = MdtRodDecoder)
+#ToolSvc += MuonMdtRawDataProviderTool
+
+#from MuonByteStream.MuonByteStreamConf import Muon__MdtRawDataProvider
+#topSequence += Muon__MdtRawDataProvider(name = "MuonMdtRawDataProvider",
+# ProviderTool = ToolSvc.MuonMdtRawDataProviderTool)
+
+
+# --- ByteStream
+theApp.Dlls += [ "ByteStreamCnvSvc", "ByteStreamCnvSvcBase"]
+# --- Services
+theApp.ExtSvc += ["ByteStreamEventStorageInputSvc/ByteStreamInputSvc"]
+theApp.ExtSvc += [ "ROBDataProviderSvc/ROBDataProviderSvc" ]
+theApp.ExtSvc += [ "EventSelectorByteStream/EventSelector"]
+theApp.EvtSel = "EventSelector"
+
+include ( "MuonEventCnvTools/MuonEventCnvTools_jobOptions.py" )
+include ( "MuonEventAthenaPool/MuonEventAthenaPool_joboptions.py" )
+#include ("MuonRdoToPrepData/MuonRdoToMuonPrepData_jobOptions.py")
+
+if('inputFiles' not in dir()):
+ inputFiles = defInputFiles
+
+print 'the input files are: '
+print inputFiles
+ServiceMgr.ByteStreamInputSvc.FullFileName= inputFiles
+
+include ("CscCalibTools/CscCalibTool_jobOptions.py")
+ToolSvc.CscCalibTool.ReadFromDatabase = True
+
+#--------------------------------------------------------------
+# Output files.
+#--------------------------------------------------------------
+if('outputFile' not in dir()):
+ outputFile = 'csc_calib_mon.root'
+
+from GaudiSvc.GaudiSvcConf import THistSvc
+ServiceMgr += THistSvc()
+ServiceMgr.THistSvc.Output = ["CscCalibMon DATAFILE='" + outputFile + "' OPT='RECREATE'"]
+
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscExtractPed.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscExtractPed.py
new file mode 100644
index 0000000000000..d81821e44f23f
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/CscExtractPed.py
@@ -0,0 +1,57 @@
+from ROOT import TCanvas, TH1D, TFile, TDirectory
+import sys
+
+
+inFileName = sys.argv[1]
+outFileName= sys.argv[2]
+outFileType = sys.argv[3]
+
+print "Input " + inFileName
+print "Output " + outFileName
+print "Ouput File Type " + outFileType
+
+cscRootName ="MUON_CSC_PED"
+overviewName = "Overview"
+
+inFile = TFile(inFileName);
+
+cnv = TCanvas("c1","c1");
+
+rootDir = TDirectory()
+overDir = TDirectory()
+
+rootDir = inFile.FindObjectAny(cscRootName)
+overDir = rootDir.FindObjectAny(overviewName)
+
+hist = TH1D()
+
+#Page 1
+hist = overDir.Get("h_csc_calib_numSignificant")
+hist.Draw()
+cnv.Print(outFileName+"(",outFileType)
+
+#Page 2
+cnv.Clear()
+cnv.Divide(2,2)
+
+i = 1
+cnv.cd(i)
+hist = overDir.Get("h_csc_calib_pedCompareOverview")
+hist.Draw()
+
+i += 1
+cnv.cd(i)
+hist = overDir.Get("h_csc_calib_noiseCompareOverview")
+hist.Draw()
+
+i += 1
+cnv.cd(i)
+hist = overDir.Get("h_csc_calib_pedChi2Overview")
+hist.Draw()
+
+i += 1
+cnv.cd(i)
+hist = overDir.Get("h_csc_calib_pedMissingChannels")
+hist.Draw()
+
+cnv.Print(outFileName+")",outFileType);
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscPulserLinearityCorrection.txt b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscPulserLinearityCorrection.txt
new file mode 100644
index 0000000000000..9ea29d26de45e
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscPulserLinearityCorrection.txt
@@ -0,0 +1,128 @@
+0 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+0 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+0 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+0 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+1 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+1 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+1 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+1 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+2 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+2 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+2 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+2 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+3 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+3 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+3 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+3 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+4 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+4 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+4 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+4 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+5 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+5 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+5 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+5 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+6 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+6 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+6 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+6 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+7 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+7 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+7 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+7 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+8 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+8 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+8 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+8 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+9 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+9 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+9 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+9 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+10 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+10 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+10 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+10 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+11 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+11 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+11 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+11 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+12 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+12 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+12 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+12 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+13 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+13 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+13 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+13 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+14 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+14 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+14 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+14 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+15 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+15 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+15 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+15 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+16 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+16 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+16 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+16 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+17 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+17 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+17 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+17 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+18 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+18 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+18 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+18 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+19 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+19 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+19 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+19 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+20 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+20 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+20 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+20 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+21 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+21 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+21 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+21 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+22 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+22 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+22 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+22 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+23 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+23 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+23 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+23 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+24 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+24 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+24 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+24 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+25 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+25 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+25 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+25 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+26 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+26 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+26 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+26 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+27 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+27 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+27 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+27 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+28 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+28 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+28 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+28 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+29 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+29 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+29 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+29 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+30 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+30 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+30 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+30 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+31 0 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+31 1 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+31 2 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
+31 3 1.0521 1.0141 1.0329 1.0169 1.0044 1.0033
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscWritePedCool.py b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscWritePedCool.py
new file mode 100755
index 0000000000000..7fd30c485b045
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/share/cscWritePedCool.py
@@ -0,0 +1,185 @@
+###############################################################
+# author Caleb Lampen lampen@physics.arizona.edu
+#
+# cscWritePedCool.py
+# Standard script to update CSC UPD1 pedestal tags
+#==============================================================
+#Input calibration text file
+if("input" not in dir()):
+ input = "/raid02/lampen/datasets/csc/pedRuns/nov/pedestal.cal"
+#Output db file name
+if("output" not in dir()):
+ output = "cscPed.db"
+#Should be "COMP200" or "OFLP200"
+if("dbname" not in dir()):
+ dbname = "COMP200"
+#Csc tag that will be used for any merging
+if("mergeCscTag" not in dir()):
+ mergeCscTag = "Csc-ES1-002-00"
+
+IOVRunStart=0
+
+#use McEventSelector
+DetDescrVersion = "ATLAS-GEONTF-08-00-00"
+
+import AthenaCommon.AtlasUnixStandardJob
+from AthenaCommon.AppMgr import theApp
+theApp.EvtMax =1
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+
+#switch Inner Detector and/or Calorimeter if needed
+# DetFlags modifications are best set here (uncomment RecExCommon_flags first)
+#from AthenaCommon.GlobalFlags import GlobalFlags
+from AthenaCommon.GlobalFlags import GlobalFlags
+#GlobalFlags.DetGeo.set_atlas()
+#GlobalFlags.DataSource.set_geant4()
+GlobalFlags.DetGeo.set_commis();
+GlobalFlags.DataSource.set_data();
+#inc ("RecExCommon/RecExCommon_flags.py")
+include ("RecExCond/RecExCommon_flags.py")
+
+from AthenaCommon.DetFlags import DetFlags
+DetFlags.detdescr.Muon_setOn()
+DetFlags.detdescr.ID_setOff()
+DetFlags.detdescr.LAr_setOff()
+DetFlags.detdescr.Tile_setOff()
+DetFlags.Print()
+
+#GlobalFlags.InputFormat.set_bytestream()
+
+include ("RecExCond/AllDet_detDescr.py")
+
+# Detector Initialization
+from AtlasGeoModel import SetGeometryVersion
+from GeoModelSvc.GeoModelSvcConf import GeoModelSvc
+from AtlasGeoModel import GeoModelInit
+
+
+
+#CscReadWRiteCoolStr is a wrapper for the CscCoolStrSvc. It will have the
+#flat file coppied to the database
+from MuonCondCool.MuonCondCoolConf import MuonCalib__CscReadWriteCoolStr
+topSequence += MuonCalib__CscReadWriteCoolStr()
+CscReadWriteDatabase= MuonCalib__CscReadWriteCoolStr()
+CscReadWriteDatabase.Write=True
+CscReadWriteDatabase.iFiles=[input]
+
+#Define cool foldernames. These need to be given to both
+#CoolStrSvc and IovDbSvc
+gainFolder = "/CSC/GAIN"
+pslopeFolder = "/CSC/PSLOPE"
+rslopeFolder = "/CSC/RSLOPE"
+pedFolder = "/CSC/PED"
+noiseFolder = "/CSC/NOISE"
+rmsFolder = "/CSC/RMS"
+f001Folder = "/CSC/FTHOLD"
+tholdFolder = "/CSC/THOLD"
+peaktFolder = "/CSC/PEAKT"
+widthFolder = "/CSC/WIDTH"
+sat1Folder = "/CSC/SAT1"
+sat2Folder = "/CSC/SAT2"
+peakcFolder = "/CSC/PEAKC"
+sampleTimeRatioFolder = "/CSC/SAMPLERATIO"
+occupancyFolder = "/CSC/OCCUPANCY"
+statFolder = "/CSC/STAT"
+
+#suspendCaching()
+
+# data is written to conditions database using OutputConditionsAlg
+from RegistrationServices.OutputConditionsAlg import OutputConditionsAlg
+OutCondAlg = OutputConditionsAlg("OutCondAlg","dummy.root")
+
+# List of things to be written.
+# Make sure that only the folders you actually want to write are
+# in this list. If something is in this list, and isn't in
+# the calibration file, the entire job will fail.
+prefix = "CondAttrListCollection#"
+OutCondAlg.ObjectList=[
+#prefix + gainFolder,
+#prefix + pslopeFolder,
+#prefix + rslopeFolder,
+prefix + pedFolder,
+prefix + noiseFolder,
+prefix + rmsFolder,
+prefix + f001Folder,
+#prefix + f001Folder,
+#prefix + tholdFolder,
+#prefix + peaktFolder,
+#prefix + widthFolder,
+#prefix + sat1Folder ,
+#prefix +sat2Folder,
+#prefix + peakcFolder,
+#prefix + sampleTimeRatioFolder,
+#prefix + occupancyFolder,
+#prefix + statFolder
+]
+
+#Taglist must be in same order as folder list!
+OutCondAlg.IOVTagList = [
+#"CscGain-sim-100-00",
+#"CscPslope-sim-000-01",
+#"CscRslope-temp-000-00",
+"CscPed-ES1-UPD1-002-00",
+"CscNoise-ES1-UPD1-002-00",
+"CscRms-ES1-UPD1-002-00",
+"CscFthold-ES1-UPD1-002-00",
+#"CscFthold-comm-ES1-UPD1-002-03",
+#"CscThold-test-000-00",
+#"CscPeakt-test-000-00",
+#"CscWidth-sim-100-00",
+#"CscSat1-sim-100-00",
+#"CscSat2-sim-100-00",
+#"CscPeakc-sim-100-00",
+#"CscSampleratio-sim-100-00",
+#"CscOccupancy-sim-100-00",
+#"CscStat-comm-002-01"
+]
+
+OutCondAlg.WriteIOV=True
+# set the interval of validity for the file here
+# putting nothing for Run2 (uppser limit) falls back on default
+#which is to be valid for all run/event
+OutCondAlg.Run1=IOVRunStart
+#OutputConditionsAlg.Event1=0
+#OutputConditionsAlg.Run2=9999
+#OutputConditionsAlg.Event2=9999
+
+##############
+from IOVDbSvc.CondDB import conddb
+conddb.setGlobalTag('COMCOND-HLTP-002-00')
+#conddb.setGlobalTag("COMCOND-006-00")
+include("RegistrationServices/IOVRegistrationSvc_jobOptions.py")
+#conddb.setGlobalTag('DEFAULTCOND')
+from AthenaCommon.AppMgr import ServiceMgr
+ServiceMgr.IOVRegistrationSvc.OverrideNames+=["Data"]
+ServiceMgr.IOVRegistrationSvc.OverrideTypes+=["String64k"]
+#WARNING! This will APPEND to the end of the old database file! Not overwrite!
+ServiceMgr.IOVDbSvc.dbConnection ="sqlite://;schema=" + output + ";dbname=" + dbname
+#Should never be changed to True unless absolutely sure!!!
+#ServiceMgr.IOVRegistrationSvc.RecreateFolders = False
+
+
+###############
+#Set detector description
+#DetDescrVersion = "ATLAS-Comm-01-00-00"
+include ("RecExCond/AllDet_detDescr.py")
+
+# CscCoolStrSvc preps data to be sent to cool database
+from MuonCondSvc.CscCondDB import cscCondDB
+#Stop caching since we aren't interested in reading it out right now
+#cscCondDB.useLocalFile()
+#cscCondDB.SetupForNewFolder()
+cscCondDB.CscCoolStrSvc.DoCaching = True
+cscCondDB.CscCoolStrSvc.DoMerge = True
+cscCondDB.addPedFolders()
+cscCondDB.addRmsFolder()
+cscCondDB.addF001Folder()
+conddb.addOverride("/CSC",mergeCscTag)
+#cscCondDB.addStatusFolder()
+#conddb.addOverride("/CSC/STAT","CscStat-comm-002-00")
+#cscCondDB.addStatusFolder()
+#ServiceMgr.MessageSvc.OutputLevel = DEBUG
+ServiceMgr.MessageSvc.debugLimit = 0
+#ServiceMgr.CscCoolStrSvc.OutputLevel = VERBOSE
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.cxx b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.cxx
new file mode 100644
index 0000000000000..2d13ef5a4598a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.cxx
@@ -0,0 +1,603 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//////////////////////////////////////////////////////
+// A programme to determine the parameters of
+// a bipolar pulse by chi2 minimization
+//
+// http://knikolop.home.cern.ch/knikolop/csc/bipolarfit.c
+//
+// for the function of the bipolar pulse look
+// http://positron.ps.uci.edu/~schernau/ROD/SIT/results/cluster1.html
+//
+// Konstantinos Nikolopoulos
+// 23/5/2007
+//
+// This is temporary
+
+// this version : 27/7/2007
+//////////////////////////////////////////////////////
+#include "BipolarFit.h"
+#include <stdio.h>
+#include <iostream>
+BipolarFit::BipolarFit()
+{
+ n=12.;
+ powcachez = -9999.;
+ powcachezn= -9999.;
+ zmax= n+1 - sqrt(n+1.0);
+ bipolarNormalization = FindPow(zmax)*(1-zmax/(n+1))*exp(-zmax);
+ tsampling = 25.;
+
+}
+
+BipolarFit::~BipolarFit()
+{}
+
+ double
+BipolarFit::FindPow(double z)
+{
+ if(fabs(powcachez-z)<1.e-4)
+ return powcachezn;
+
+ double zpower = z*z*z;
+ zpower *= zpower;
+ zpower *= zpower;
+
+ powcachez = z;
+ powcachezn= zpower;
+ return zpower;
+}
+
+ double
+BipolarFit::FindInitValues(double*x,double *initValues,int *maxsample)
+{
+ // find maximum sample imax:
+ double peakingTime=-99.; // interpolated peaking time in samples
+ double amplitude=-99.; // interpolated amplitude
+ double amin, amax;
+ int imax = 0;
+ const int numSample=4;
+ amax = x[0];
+ amin = x[0];
+ for(int j=1; j<numSample; j++)
+ {
+ if(amax<x[j])
+ {
+ amax = x[j];
+ imax = j;
+ }
+ if(amin>x[j])
+ amin = x[j];
+ }
+
+ // calculate peak and amplitude:
+ (*maxsample)=imax;
+ if((imax==0) || (imax==numSample-1)) // no interpolation possible
+ {
+ peakingTime = imax;
+ amplitude = amax;
+ }
+ else // do parabolic interpolation
+ {
+ double a, b, c; // coeffients of parabola y=a*x*x+b*x+c
+ a = 0.5*(x[imax+1]+x[imax-1]-2.0*x[imax]);
+ b = 0.5*(x[imax+1]-x[imax-1]);
+ c = x[imax];
+
+ peakingTime = -0.5*b/a;
+ amplitude = a*peakingTime*peakingTime + b*peakingTime + c;
+ peakingTime += imax; // it was relative to imax
+ }
+
+ initValues[0] = amplitude;
+ initValues[1] = peakingTime - zmax*initValues[2]/tsampling;
+ return x[imax];
+}
+
+ double
+BipolarFit::bipolar(double *x, double *parm) // the bipolar pulse function
+{
+ // the bipolar pulse function is
+ //
+ // z = (x-parm[0])*tsampling/parm[3]
+ // zmax = n+1 - sqrt(n+1)
+ // aa = exp(n*log(zmax))*(1-zmax/(n+1))*exp(-zmax)
+ // parm[0]*exp(n*log(z))*(1-z/(n+1))*exp(-z)/aa
+ //
+ // for timing reasons instead of x (ie # of sample)
+ // the z is given
+ double z = x[0];
+
+
+ //const double tsampling = 25.;// nsec
+ //z=(x[0]-parm[1])*tsampling/(parm[2]);
+ if(z<0.)
+ return 0.;
+
+ return parm[0]*FindPow(z)*(1-z/(n+1))*exp(-z)/bipolarNormalization;
+}
+
+ void
+BipolarFit::Derivative(double A[][3],double fp[][1], double p0[][1],int imeas, int*meas)
+{
+ //calculate the derivatives and the 0th order approximation
+ //around the ADC samplings
+ double norm = p0[0][0];
+ //double parm[3]={1.,p0[1][0],p0[2][0]};
+ for(int i=0;i<imeas;i++)
+ {
+ int ii = meas[i];
+ double z = (ii-p0[1][0])*tsampling/p0[2][0];
+ double repquant = 0.;
+ double dFdzNormalized = 0.;
+ if(z>0.)
+ {
+ repquant = FindPow(z)*exp(-z)/bipolarNormalization;
+ dFdzNormalized= repquant*(n/z+z/13.-2.);
+ }
+
+ A[ii][0] = repquant*(1.-z/13.);
+ //A[ii][0] = bipolar(&z,parm);
+ fp[ii][0] = norm * A[ii][0];
+
+ //double normOverZmax = norm/bipolarNormalization;
+ double commonpart = norm* dFdzNormalized;//(z,parm);
+ A[ii][1] = commonpart * (-tsampling/p0[2][0]);
+ A[ii][2] = commonpart * (-z/p0[2][0]);
+ }
+ // end of derivative/zeroth order calculations
+}
+
+//void BipolarFit::InvertMatrix(HepMatrix& matrix,const int dim,int*correspdim)
+void BipolarFit::InvertMatrix(double matrix[][3],const int dim,int*correspdim)
+{
+ // invert 2x2 or 3x3 symmetric matrix
+ if(dim==2)
+ {
+ int ii=correspdim[0];
+ int jj=correspdim[1];
+ double determinant= -matrix[jj][ii]*matrix[ii][jj] +matrix[ii][ii]*matrix[jj][jj];
+ //if(fabs(determinant)<1.e-13)
+ //std::cout<<" zero determinant "<<std::endl;
+ double i00 = matrix[ii][ii];
+ matrix[ii][ii] = matrix[jj][jj]/determinant;
+ matrix[jj][jj] = i00/determinant;
+ matrix[ii][jj] = -matrix[ii][jj]/determinant;
+ matrix[jj][ii] = matrix[ii][jj];
+
+ }
+ else if(dim==3)
+ {
+ double sm12 = matrix[0][1]*matrix[0][1];
+ double sm13 = matrix[0][2]*matrix[0][2];
+ double sm23 = matrix[1][2]*matrix[1][2];
+ double determinant = matrix[0][0]*matrix[1][1]*matrix[2][2]
+ -matrix[0][0]*sm23-sm12*matrix[2][2]
+ +2.*matrix[0][1]*matrix[0][2]*matrix[1][2]
+ -matrix[1][1]*sm13;
+
+ //if(fabs(determinant)<1.e-13)
+ //std::cout << "zero determinant"<<std::endl;
+ double i00 = matrix[1][1]*matrix[2][2]-sm23;
+ double i11 = matrix[0][0]*matrix[2][2]-sm13;
+ double i22 = matrix[0][0]*matrix[1][1]-sm12;
+
+ matrix[1][0] = (matrix[0][2]*matrix[1][2]-matrix[2][2]*matrix[0][1])/determinant;
+ matrix[2][0] = (matrix[0][1]*matrix[1][2]-matrix[0][2]*matrix[1][1])/determinant;
+ matrix[2][1] = (matrix[0][1]*matrix[0][2]-matrix[0][0]*matrix[1][2])/determinant;
+ matrix[0][1] = matrix[1][0];
+ matrix[0][2] = matrix[2][0];
+ matrix[1][2] = matrix[2][1];
+ matrix[0][0] = i00/determinant;
+ matrix[1][1] = i11/determinant;
+ matrix[2][2] = i22/determinant;
+
+ }
+ else
+ {
+ //int ierr;
+ //matrix.invert(ierr);
+ printf("this is not a 2x2 or 3x3 matrix\n");
+ }
+}
+
+ int
+BipolarFit::Fit(double *x,const double ex,const double pedestal, const double predefinedwidth,double *result,double * /**errors*/,double *chi2)
+{
+ /**
+Input :
+x -> the ADC measurements in the 4 samplings
+ex -> the electronic noise in ADC counts
+Output:
+result -> the result of the fit
+errors -> the errors of the fit parameters
+chi2 -> the chi2 of the fit
+Returns:
+status of the fit
+
+The function return an integer representing different exit status.
+--- Bipolar fit tried ---
+0 --> Bipolar fit converged OK
+1 --> The fit converged but its chi2 is bad
+2 --> The fit converged but estimated shaping time very different from expected
+3 --> Maximum number of iterations reached OK
+4 --> Diverging fit. Return parabola interpolation OK
+5 --> One sample saturated. Bipolar fit tried nevertheless
+ //6 --> Late pulse the. The constant shaping time fit converged but the full fit didn't.
+ //Return fix shaping time fit results. OK
+ --- Bipolar fit not tried ---
+ 6 --> More samples saturated. Return estimate based on parabola interpolation
+ 7 --> Parabola interpolation indicates amplitude <5sigmaNoise.
+ Return parabola estimated parameters. OK
+ 8 --> Distorted pulse shape. Fit likely to fail if attempted.parabola interpolation returned
+ 9 --> Too late pulse. Bipolar fit not attempted, return rough estimate OK
+ 10 --> It seems that a reasonable pulse exist, but the fit gave substantially lower value than the parabola. return the parabola (which is kind of lower bound for normal pulses).
+ **/
+
+ // initial parameter estimates using parabola interpolation
+ double initValues[3]={0.};
+ //double FWHM=0.;
+ int imax = -1;
+ initValues[2]=predefinedwidth;
+ double samplemax = FindInitValues(x,initValues,&imax);
+ //std::cout << " Init " << initValues[0] << " "<<initValues[1] << " "<<initValues[2]<<std::endl;
+ result[0] = initValues[0];
+ result[1] = initValues[1];
+ result[2] = initValues[2];
+
+ // do not fit noise
+ if(samplemax < 5.*ex)
+ return 7;
+
+ bool onesaturated = false;
+ if(imax==0)
+ {
+ if(x[imax+1]<0. && x[imax+2]<0.)
+ return 8;
+ }
+ else if(imax==3)
+ {
+ // don't fit too late pulses
+ if(initValues[1]+zmax*initValues[2]/tsampling>2.75)
+ return 9;
+ }
+ else
+ {
+ if(x[imax-1]<0. && (x[imax+1]<2. || -x[imax-1]-x[imax]>0.))
+ return 8;
+ if(x[imax]+pedestal>4000.)
+ {
+ if(x[imax-1]+pedestal>4000. || x[imax+1]+pedestal>4000.)
+ return 6;
+ else
+ {
+ onesaturated=true;
+ }
+ }
+ }
+ //always fix width and fit two parameters
+ bool fitpar[3] = {true,true,false};
+ bool usemeas[4] = {true,true,true,true};
+ if(initValues[1]+zmax*initValues[2]/tsampling<2.0)
+ {
+ fitpar[2] = false;
+ usemeas[3]= false;
+ }
+ //if(initValues[1]>0.)
+ //{
+ // usemeas[0]= false;
+ // fitpar[2] = false;
+ //}
+
+ if(onesaturated)
+ {
+ usemeas[imax]= false;
+ fitpar[2] = false;
+ }
+ int imeas =0;
+ int meas[4] = {0};
+ for(int i=0;i<4;i++)
+ if(usemeas[i])
+ {
+ meas[imeas]=i;
+ imeas++;
+ }
+ int ipar =0;
+ int par[3] = {0};
+ for(int i=0;i<3;i++)
+ if(fitpar[i])
+ {
+ par[ipar]=i;
+ ipar++;
+ }
+
+ /* std::cout << " Use Meas " <<std::endl;
+ for(int i =0;i<4;i++)
+ std::cout << usemeas[i]<< " ";
+ std::cout<<std::endl << " fit par " <<std::endl;
+ for(int i=0;i<3;i++)
+ std::cout << fitpar[i] << " ";
+ std::cout<<std::endl;
+ */
+ int FitStatus = TheFitter(x,ex,initValues,imeas,meas,ipar,par,chi2,result);
+ //std::cout << "chi2 " << (*chi2) <<std::endl;
+ //std::cout << result[0] << " " << result[1] << " "<<result[2] <<std::endl;
+ //std::cout<< initValues[0] << " " << initValues[1] << " "<< initValues[2]<<std::endl;
+ // the parabola interpolated estimate is most of the time a lower bound (for high pulses)!
+ if(result[0]> 10.*ex && result[0]<0.90*initValues[0])
+ {
+ result[0] = initValues[0];
+ // std::cout << "using the lower bound "<<std::endl;
+ return 10;
+ }
+
+ if(onesaturated)
+ return 5;
+
+
+ return FitStatus;
+}
+
+ void
+BipolarFit::InvertSymmetric4x4(double W[][4])
+{
+ double Determinant = W[0][3]*W[0][3]*W[1][2]*W[1][2] - 2.*W[0][2]*W[0][3]*W[1][2]*W[1][3]
+ +W[0][2]*W[0][2]*W[1][3]*W[1][3]-W[0][3]*W[0][3]*W[1][1]*W[2][2] + 2.*W[0][1]*W[0][3]*W[1][3]*W[2][2] - W[0][0]*W[1][3]*W[1][3]*W[2][2]+2.*W[0][2]*W[0][3]*W[1][1]*W[2][3]
+ - 2.*W[0][1]*W[0][3]*W[1][2]*W[2][3] - 2.*W[0][1]*W[0][2]*W[1][3]*W[2][3]+2.*W[0][0]*W[1][2]*W[1][3]*W[2][3]+W[0][1]*W[0][1]*W[2][3]*W[2][3]-W[0][0]*W[1][1]*W[2][3]*W[2][3]-W[0][2]*W[0][2]*W[1][1]*W[3][3]+2.*W[0][1]*W[0][2]*W[1][2]*W[3][3]-W[0][0]*W[1][2]*W[1][2]*W[3][3]-W[0][1]*W[0][1]*W[2][2]*W[3][3]+W[0][0]*W[1][1]*W[2][2]*W[3][3];
+
+ W[0][1] = W[3][0]*W[3][1]*W[2][2]-W[3][0]*W[2][1]*W[3][2] - W[2][0]*W[3][1]*W[3][2]+W[1][0]*W[3][2]*W[3][2]+W[2][0]*W[2][1]*W[3][3]-W[1][0]*W[2][2]*W[3][3];
+ W[0][2] = -W[3][0]*W[2][1]*W[3][1]+W[2][0]*W[3][1]*W[3][1]+W[3][0]*W[1][1]*W[3][2]-W[1][0]*W[3][1]*W[3][2]-W[2][0]*W[1][1]*W[3][3]+W[1][0]*W[2][1]*W[3][3];
+ W[0][3] = W[3][0]*W[2][1]*W[2][1]-W[2][0]*W[2][1]*W[3][1]-W[3][0]*W[1][1]*W[2][2]+W[1][0]*W[3][1]*W[2][2]+W[2][0]*W[1][1]*W[3][2]-W[1][0]*W[2][1]*W[3][2];
+ W[1][2] = W[3][0]*W[3][0]*W[2][1]-W[2][0]*W[3][0]*W[3][1]-W[1][0]*W[3][0]*W[3][2]+W[0][0]*W[3][1]*W[3][2]+W[1][0]*W[2][0]*W[3][3]-W[0][0]*W[2][1]*W[3][3];
+
+ W[1][3] = -W[2][0]*W[3][0]*W[2][1]+W[2][0]*W[2][0]*W[3][1]+W[1][0]*W[3][0]*W[2][2]-W[0][0]*W[3][1]*W[2][2]-W[1][0]*W[2][0]*W[3][2]+W[0][0]*W[2][1]*W[3][2];
+ W[2][3] = W[2][0]*W[3][0]*W[1][1]-W[1][0]*W[3][0]*W[2][1]-W[1][0]*W[2][0]*W[3][1]+W[0][0]*W[2][1]*W[3][1]+W[1][0]*W[1][0]*W[3][2]-W[0][0]*W[1][1]*W[3][2];
+
+ double W00 = -W[3][1]*W[3][1]*W[2][2]+2.*W[2][1]*W[3][1]*W[3][2]-W[1][1]*W[3][2]*W[3][2]-W[2][1]*W[2][1]*W[3][3]+W[1][1]*W[2][2]*W[3][3];
+ double W11 = -W[3][0]*W[3][0]*W[2][2]+2.*W[2][0]*W[3][0]*W[3][2]-W[0][0]*W[3][2]*W[3][2]-W[2][0]*W[2][0]*W[3][3]+W[0][0]*W[2][2]*W[3][3];
+ double W22 = -W[3][0]*W[3][0]*W[1][1]+2.*W[1][0]*W[3][0]*W[3][1]-W[0][0]*W[3][1]*W[3][1]-W[1][0]*W[1][0]*W[3][3]+W[0][0]*W[1][1]*W[3][3];
+ double W33 = -W[2][0]*W[2][0]*W[1][1]+2.*W[1][0]*W[2][0]*W[2][1]-W[0][0]*W[2][1]*W[2][1]-W[1][0]*W[1][0]*W[2][2]+W[0][0]*W[1][1]*W[2][2];
+
+ for(int i=0;i<3;i++)
+ {
+ for(int j=1;j<4;j++)
+ {
+ if(i>=j)
+ continue;
+ W[i][j] = W[i][j]/Determinant;
+ W[j][i] = W[i][j];
+ }
+ }
+ W[0][0] = W00/Determinant;
+ W[1][1] = W11/Determinant;
+ W[2][2] = W22/Determinant;
+ W[3][3] = W33/Determinant;
+}
+
+ int
+BipolarFit::TheFitter(double*x,const double ex,double *initValues, int imeas, int *meas, int ipar, int *par,double *chi2,double *result)
+{
+ // maximum iterations
+ const int maxIter = 7;
+ // tolerances
+ double fitTolerance0 = 0.10;
+ double fitTolerance1 = 0.01;
+ //HepMatrix p0(3,1,0); // the matrix of the initial fit parameters
+ double p0[3][1]={{0.}};
+ for(int j=0;j<3;j++)
+ p0[j][0] = initValues[j];
+
+ //HepMatrix m(4,1,0); // the matrix of ADC measurements (samples: 0,1,2,3)
+ double m[4][1]={{0.}};
+ //HepMatrix W(4,4,0); // the error matrix of the ADC measurements
+ double W[4][4]={{0.}};
+ for(int i=0;i<4;i++)
+ {
+ m[i][0] = x[i];
+ W[i][i] = ex*ex;
+ }
+ // covariances
+ W[0][1] = 0.03*ex*ex;
+ W[0][2] = -0.411*ex*ex;
+ W[0][3] = -0.188*ex*ex;
+ W[1][2] = 0.0275*ex*ex;
+ W[1][3] = -0.4303*ex*ex;
+ W[2][3] = 0.*ex*ex;
+ W[1][0] = W[0][1];
+ W[2][0] = W[0][2];
+ W[3][0] = W[0][3];
+ W[2][1] = W[1][2];
+ W[3][1] = W[1][3];
+ W[3][2] = W[2][3];
+
+ //WW.invert(ierr);
+ InvertSymmetric4x4(W);
+
+ // Taylor expansion of the bipolar pulse model around the
+ // samplings : F(x) = F(p0) + A *(p-p0) + higher.order
+ //HepMatrix fp(4,1,0); // the matrix of 0th order approximation
+ double fp[4][1]={{0.}};
+ //HepMatrix A(4,3,0); // the matrix of derivatives
+ double A[4][3]={{0.}};
+ // remarks :
+ // if the pulse peaks in the last sampling fit with a constant shaping time
+ // if the pulse peaks in the first sampling fit without using the last sampling
+ // (too large contribution to the chi2
+ int counter=0;
+ bool converged=false;
+ double amplitudeChangeOld = 0.;
+ bool diverganceCandidate = false;
+ //HepMatrix weight(3,3,1); // weight matrix allocated once
+ // the non-fitted parts are taken care appropriately
+ // at least if the fitting parameters or measurements
+ // don't change during the fitting procedure
+ double weight[3][3]={{0.}};
+ weight[0][0]=1.;
+ weight[1][1]=1.;
+ weight[2][2]=1.;
+
+ //HepMatrix residFactor(3,4,0); // residFactor allocated once
+ double residFactor[3][4]={{0.}};
+ while(!converged && counter<maxIter) // fit loop
+ {
+ Derivative(A,fp,p0,imeas,meas);// calculate the matrix of derivatives and 0th order approximation
+ // matrix multiplication
+ // the weight matrix is symmetric
+ // weight= A.T()*W*A;//.assign(A.T()*W*A);
+
+ double helpmatrix[4][3]={{0.}};
+ for(int i=0;i<imeas;i++)
+ {
+ int ii=meas[i];
+ for(int j=0;j<ipar;j++)
+ {
+ int jj=par[j];
+ for(int k=0;k<imeas;k++)
+ {
+ int kk=meas[k];
+ helpmatrix[ii][jj] += W[ii][kk]*A[kk][jj];
+ }
+ }
+ }
+ for(int i=0;i<ipar;i++)
+ {
+ int ii=par[i];
+ for(int j=i;j<ipar;j++)
+ {
+ int jj=par[j];
+ weight[ii][jj] = 0.;
+ for(int k=0;k<imeas;k++)
+ {
+ int kk=meas[k];
+ weight[ii][jj] += A[kk][ii]*helpmatrix[kk][jj];//A[kk][ii]*A[kk][jj];
+ }
+ //weight[ii][jj]*=W[0][0];
+ weight[jj][ii] =weight[ii][jj];
+ }
+ }
+ //weight.invert(ierr); // inversion of weight matrix
+ // hand-made inversion of 2x2 or 3x3 symmetric matrix
+ InvertMatrix(weight,ipar,par);
+
+ //calculate W*(A.T()*W)
+ //residFactor = weight*(A.T()*W);
+ double helpmatrix2[3][4]={{0.}};
+ for(int i=0;i<ipar;i++)
+ {
+ int ii=par[i];
+ for(int j=0;j<imeas;j++)
+ {
+ int jj=meas[j];
+ for(int k=0;k<imeas;k++)
+ {
+ int kk=meas[k];
+ helpmatrix2[ii][jj] += A[kk][ii] * W[kk][jj];
+ }
+ }
+ }
+
+ for(int i=0;i<ipar;i++)
+ {
+ int ii = par[i];
+ for(int j=0;j<imeas;j++)
+ {
+ int jj=meas[j];
+ residFactor[ii][jj]=0.;
+ for(int k=0;k<ipar;k++)
+ {
+ int kk=par[k];
+ residFactor[ii][jj] += weight[ii][kk]*helpmatrix2[kk][jj];
+ }
+ //residFactor[ii][jj]*=W[0][0];
+ }
+ }
+
+ double paramDiff[3] = {0.};
+ for(int i=0;i<ipar;i++)
+ {
+ int ii=par[i];
+ //estimation of new parameters
+ //paramDiff[i][0] += (weight*(A.T()*W)*(m-fp))[i][0];
+ for(int j=0;j<imeas;j++)
+ {
+ int jj = meas[j];
+ paramDiff[ii] += residFactor[ii][jj]*(m[jj][0]-fp[jj][0]);
+ }
+ p0[ii][0] += paramDiff[ii];
+ }
+ //std::cout << "##### "<<p0[0][0]<< " "<<p0[1][0] << " "<<p0[2][0]<<std::endl;
+ // if the parameters are going nuts keep them sensible
+ if(p0[1][0]>1. || p0[1][0]<-3.)
+ p0[1][0] = initValues[1];
+
+ double amplitudeChangeNew = fabs(paramDiff[0]);
+ if(fabs(paramDiff[0])<fitTolerance0 && fabs(paramDiff[1])<fitTolerance1)
+ {
+ converged = true;
+ // calculate chi2
+ // (m-fp).T()*W*(m-fp)
+ double residual[4]= {0.};
+ for(int i=0;i<imeas;i++)
+ {
+ int ii=meas[i];
+ residual[i] = m[ii][0]-fp[ii][0];
+ }
+
+ double tmpChi2 = 0.;
+ double helpmatrixchi2[4][1]={{0.}};
+ for(int i=0;i<imeas;i++)
+ {
+ int ii=meas[i];
+ for(int k=0;k<imeas;k++)
+ {
+ int kk=meas[k];
+ helpmatrixchi2[ii][0] += W[ii][kk]*residual[kk];
+ }
+ }
+ for(int k=0;k<imeas;k++)
+ {
+ int kk=meas[k];
+ tmpChi2 += residual[kk]*helpmatrixchi2[kk][0];
+ }
+ (*chi2) = tmpChi2;
+ }
+ else if(counter>4 && (amplitudeChangeNew>4.*amplitudeChangeOld))
+ {
+ if(diverganceCandidate)
+ {
+ //diverging fit
+ //return parabola interpolation
+ printf("Diverging fit\n");
+ return 4;
+ }
+ else
+ diverganceCandidate = true;
+ }
+ if(p0[0][0]<0.)
+ {
+ //negative amplitude
+ //fit diverged
+ // return parabola
+ return 4;
+ }
+ //if after a couple of iterations the amplitude is low
+ // reduce the tolerances (or the maximum iterations)
+ // low amplitude pulses tend to oscillate and exhaust all iterations
+ if(p0[0][0]<20.)
+ {
+ fitTolerance0 = 0.1;
+ fitTolerance1 = 0.05;
+ }
+ amplitudeChangeOld = amplitudeChangeNew;
+ counter++;
+ }
+ result[0]=p0[0][0];
+ result[1]=zmax*p0[2][0]/tsampling+p0[1][0];
+ result[2]=p0[2][0];
+
+ if(counter==maxIter)
+ return 3;
+ return 0;
+}
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.h b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.h
new file mode 100644
index 0000000000000..c11491ff09abb
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/BipolarFit.h
@@ -0,0 +1,48 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef BIPOLARFIT_H
+#define BIPOLARFIT_H
+//////////////////////////////////////////////////////
+// A programme to determine the parameters of
+// a bipolar pulse by chi2 minimization
+//
+// http://knikolop.home.cern.ch/knikolop/csc/bipolarfit.c
+//
+// for the function of the bipolar pulse look
+// http://positron.ps.uci.edu/~schernau/ROD/SIT/results/cluster1.html
+//
+// Konstantinos Nikolopoulos
+// 23/5/2007
+// this version : 27/7/2007
+//////////////////////////////////////////////////////
+#include<math.h>
+//#include"CLHEP/Matrix/Matrix.h"
+//#include"CLHEP/Matrix/SymMatrix.h"
+class BipolarFit
+{
+ public:
+ BipolarFit();
+ ~BipolarFit();
+ int Fit(double *x,const double ex,const double pedestal, const double predefinedwidth,double *result,double *errors,double *chi2);
+
+ private:
+ void InvertMatrix(double matrix[][3],const int dim,int*);
+ void InvertSymmetric4x4(double W[][4]);
+ double FindInitValues(double*x,double *initValues,int *maxsample);
+ int TheFitter(double*x,const double ex,double *initValues, int imeas, int *meas, int ipar, int *par,double *chi2,double *result);
+ //void Derivative(HepMatrix&,HepMatrix&,HepMatrix&,bool*usemeas);
+ void Derivative(double A[][3],double fp[][1], double p0[][1],int imeas, int*meas);
+ //double dFdz(double,double *);
+ double bipolar(double*, double*);
+ double FindPow(double z);
+
+ double n;
+ double powcachez;
+ double powcachezn;
+ double zmax;
+ double bipolarNormalization;
+ double tsampling;
+};
+#endif
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.cxx b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.cxx
new file mode 100644
index 0000000000000..3310b8296b138
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.cxx
@@ -0,0 +1,1111 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "GaudiKernel/Chrono.h"
+
+#include "TH1.h"
+#include "TF1.h"
+#include "TFile.h"
+
+#include "StoreGate/StoreGate.h"
+
+#include "MuonIdHelpers/CscIdHelper.h"
+#include "MuonCSC_CnvTools/ICSC_RDO_Decoder.h"
+
+#include "MuonRDO/CscRawData.h"
+#include "MuonRDO/CscRawDataCollection.h"
+#include "MuonRDO/CscRawDataContainer.h"
+
+#include "CscCalibData/CscCalibResultContainer.h"
+#include "CscCalibData/CscCalibReportContainer.h"
+#include "CscCalibData/CscCalibReportPed.h"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <bitset>
+#include <inttypes.h>
+//#include <math>
+
+#include "CscCalcPed.h"
+
+using namespace std;
+
+namespace MuonCalib {
+
+ CscCalcPed::CscCalcPed(const string& name, ISvcLocator* pSvcLocator) :
+ Algorithm(name,pSvcLocator),
+ m_chronoSvc(0),
+ m_cscCoolStrSvc("MuonCalib::CscCoolStrSvc",name),
+ m_cscRdoDecoderTool ("Muon::CscRDO_Decoder"),
+ m_numBits(12),
+ m_bitHists(NULL),
+ m_bitProds(NULL),
+ m_eventCnt(0),
+ m_debug(0),
+ m_verbose(0)
+ {
+ declareProperty("OutputFile", m_outputFileName = "output.cal");
+ declareProperty("TitlePrefix",m_titlePrefix = ""); //Prefix appended to title of histograms and graphs
+ declareProperty("TitlePostfix",m_titlePostfix = ""); //Prefix appended to title of histograms and graphs
+ declareProperty("CscRdoDecoderTool", m_cscRdoDecoderTool );
+ declareProperty("PedAmpHighBound", m_ampHistHighBound = 2600);
+ declareProperty("PedAmpLowBound", m_ampHistLowBound = 1800);
+ declareProperty("PedAmpNumBins", m_ampHistNumBins = 800);
+ //declareProperty("PedAmpHighBound", m_ampHistHighBound = 2300);
+ //declareProperty("PedAmpLowBound", m_ampHistLowBound = 1900);
+ //declareProperty("PedAmpNumBins", m_ampHistNumBins = 400);
+ declareProperty("ExpectedChamberLayer", m_expectedChamberLayer = 2);
+ declareProperty("ThresholdMultiplier", m_thresholdMultiplier = 3.1);
+ declareProperty("DoBitHists", m_doBitHists = false);
+ declareProperty("CalOutputVersion", m_calOutputVersion="03-00");
+ declareProperty("DoCorrelation", m_doCorrelation = false);
+
+ declareProperty("DoSampHists", m_doSampleHists = false);
+ declareProperty("NumSamplesExpected", m_numSamplesExpected = 4);
+ declareProperty("DoF001", m_doF001 = true);
+
+ declareProperty("OnlineCalibFile", m_onlineDbFile = "", "File with data currently stored in online configuration database");
+ declareProperty("CompareOnlineCalibFile",m_doOnlineDbFile = true, "Compare new to online data?");
+
+
+ //Can't do correlation without bitHists
+ if(!m_doBitHists){
+ m_doCorrelation = false;
+ //MsgStream mLog( msgSvc(), name() );
+ // mLog << MSG::WARNING
+ // << "Tried to do correlation without bitHists. This isn't possible." << endreq;
+ }
+
+ }
+
+ CscCalcPed::~CscCalcPed() {}
+
+ StatusCode CscCalcPed::initialize()
+ {
+ MsgStream mLog( msgSvc(), name() );
+
+ //Set debug level
+ m_debug = (mLog.level() <= MSG::DEBUG);
+ m_verbose = (mLog.level() <= MSG::VERBOSE);
+
+ mLog << MSG::INFO << "CscCalcPed::initialize() called" << endreq;
+
+ if(m_doOnlineDbFile && m_onlineDbFile ==""){
+ mLog << MSG::FATAL <<"Either specify an OnlineCalibFile or set CompareOnlineCalibFile to false"<< endreq;
+ return StatusCode::FAILURE;
+ }
+
+ //*******Register services and tools *********/
+ /// Store Gate active store
+ StatusCode sc = serviceLocator()->service("StoreGateSvc", m_storeGate);
+ if (sc != StatusCode::SUCCESS )
+ {
+ mLog << MSG::ERROR << " Cannot get StoreGateSvc " << endreq;
+ return sc ;
+ }
+
+ sc = serviceLocator()->service("DetectorStore", m_detStore);
+
+ if(sc.isSuccess())
+ {
+ sc = m_detStore->retrieve(m_cscId,"CSCIDHELPER");
+ if( sc.isFailure())
+ {
+ mLog << MSG::ERROR << " Cannot retrieve CscIdHelper " << endreq;
+ return sc;
+ }
+ }
+ else
+ {
+ mLog << MSG::ERROR << "couldn't retrieve CscIdHelper" << endreq;
+
+ return StatusCode::FAILURE;
+ }
+
+ if(m_cscCoolStrSvc.retrieve().isFailure())
+ {
+ mLog << MSG::FATAL << "Unable to retrieve CscCoolStrSvc" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ sc = service("ChronoStatSvc",m_chronoSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Cannot retrieve ChronoStatSvc!" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ if (m_cscRdoDecoderTool.retrieve().isFailure()){
+ mLog << MSG::FATAL << "Cannot retrieve Csc_Rdo_Decoder Tool!" << endreq;
+ return StatusCode::FAILURE;
+ }
+ /*sc = service("THistSvc", m_thistSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Cannot retrieve THistSvc!" << endreq;
+ return StatusCode::FAILURE;
+ }*/
+
+ IToolSvc* toolSvc;
+ sc = service("ToolSvc",toolSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::ERROR << "Unable to retrieve ToolSvc" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ //Set to SG::VIEW_ELEMENTS, since we want root to do its own memory
+ //management.
+ m_ampHists = new DataVector<TH1I>(SG::VIEW_ELEMENTS);
+ if(m_doSampleHists) m_sampHists = new DataVector< DataVector<TH1I> >;
+ if(m_doBitHists) m_bitHists = new DataVector<TH1I>(SG::VIEW_ELEMENTS);
+
+ //Loop through ids to find out what hash range we're working on (in case we're using some
+ //unusual geometry)
+ vector<Identifier> ids = m_cscId->idVector();
+ vector<Identifier>::const_iterator chamItr = ids.begin();
+ vector<Identifier>::const_iterator chamEnd = ids.end();
+ m_maxStripHash = 0;
+
+ for(; chamItr != chamEnd; chamItr++)
+ {
+ vector<Identifier> stripVect;
+ m_cscId->idChannels(*chamItr,stripVect);
+
+ vector<Identifier>::const_iterator stripItr = stripVect.begin();
+ vector<Identifier>::const_iterator stripEnd = stripVect.end();
+ for(;stripItr != stripEnd; stripItr++)
+ {
+ IdentifierHash stripHash;
+ m_cscId->get_channel_hash(*stripItr,stripHash);
+ if((unsigned int)stripHash > m_maxStripHash)
+ m_maxStripHash = (int)stripHash;
+ }//End strip loop
+ }//end chamber loop
+
+
+ //Now creating ampHists. It wasn't done in last loop since there
+ //is no gaurantee that it will come out in strip order, and we assume
+ //later that m_ampHists's index = stripHash
+ if(m_debug) mLog << MSG::DEBUG << "Preparing ampHists. Only allowing those for "
+ << " chamberLayer " << m_expectedChamberLayer << endreq;
+ for(unsigned int stripItr = 0 ; stripItr <=m_maxStripHash; stripItr++)
+ {
+
+ IdentifierHash stripHash =stripItr;
+ Identifier stripId;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id(stripHash, stripId, &channelContext);
+
+ int chamLayer = m_cscId->chamberLayer(stripId);
+ if(chamLayer == m_expectedChamberLayer) //Only second chamber layer exists
+ {
+ int stationEta = m_cscId->stationEta(stripId);
+ int stationPhi = m_cscId->stationPhi(stripId);
+ int stripNumber = m_cscId->strip(stripId);
+ int wireLayer = m_cscId->wireLayer(stripId);
+ char orientation = m_cscId->measuresPhi(stripId) ? 'Y':'X';
+
+
+ char name[30],titleSeed[600];
+ TH1I* hist = NULL;
+
+ int stationName = m_cscId->stationName(stripId);
+ //Amplitude histogram
+ sprintf(name, "ampHist%d",stripItr);
+ sprintf(titleSeed, "Amplitude Histogram for eta %d, sector %d, layer %d%c, strip %d",
+ stationEta,(2*stationPhi+50 - stationName),wireLayer,orientation,stripNumber);
+ string title = m_titlePrefix + titleSeed + m_titlePostfix;
+
+ hist = new TH1I(name,title.c_str(),m_ampHistNumBins,m_ampHistLowBound,
+ m_ampHistHighBound);
+ hist->GetXaxis()->SetTitle("Amplitude (ADC value)");
+ hist->GetYaxis()->SetTitle("Counts");
+ m_ampHists->push_back(hist);
+
+ if(m_doSampleHists) {
+ DataVector<TH1I>* tempVect = new DataVector<TH1I>(SG::VIEW_ELEMENTS);
+ for(int cnt = 0; cnt < m_numSamplesExpected ; cnt++) {
+ sprintf(name, "sampHist%d_%d",stripItr,cnt);
+ sprintf(titleSeed, "Amplitude Histogram for eta %d, sector %d, layer %d%c, strip %d, sample %d",
+ stationEta,(2*stationPhi+50 - stationName),wireLayer,orientation,stripNumber,cnt);
+
+ hist = new TH1I(name,title.c_str(),m_ampHistNumBins,m_ampHistLowBound,
+ m_ampHistHighBound);
+ tempVect->push_back(hist);
+
+ }
+ m_sampHists->push_back(tempVect);
+ }
+
+ if(m_bitHists)
+ {
+ //Bit histogram (for looking for stuck-bits)
+ sprintf(name, "bitHist%d",stripItr);
+ sprintf(titleSeed, "Bit histogram for eta %d, sector %d, layer %d%c strip %d",
+ stationEta,(2*stationPhi+50 - stationName),wireLayer,orientation,stripNumber);
+ title = m_titlePrefix + titleSeed + m_titlePostfix;
+ hist = new TH1I(name, title.c_str(), m_numBits, 0, m_numBits); //12 bits
+ hist->GetXaxis()->SetTitle("Bit");
+ hist->GetYaxis()->SetTitle("Counts");
+ m_bitHists->push_back(hist);
+ if(m_doCorrelation) {
+ (*m_bitProds)[stripItr] = new TH2F(TString::Format("bitProds%d",stripItr),"Bit products", m_numBits,0,m_numBits, m_numBits, 0, m_numBits);
+
+ }
+ }
+ }
+ else
+ {
+ m_ampHists->push_back(NULL);
+ if(m_bitHists) m_bitHists->push_back(NULL);
+ if(m_doSampleHists) m_sampHists->push_back(NULL);
+ }
+ }//end strip loop
+
+ //Setup result collections. These will be passed to storegate for any monitoring later
+ m_peds = new CscCalibResultCollection("ped");
+ m_noises = new CscCalibResultCollection("noise");
+ m_rmses = new CscCalibResultCollection("rms");
+ if(m_doF001){
+ if(m_debug) mLog << MSG::DEBUG << "Doing f001" << endreq;
+ //For f001 values
+ m_f001s = new CscCalibResultCollection("f001");
+
+ //Initializing for comparing new values to an online database file
+ if(m_doOnlineDbFile){
+
+ //How many samples failed the online threshold test of f001 +2*RMS
+ //(f001 and RMS read from a file from online configuration db)
+ m_onlineTHoldBreaches = new CscCalibResultCollection("OnlTHoldBreaches");
+
+
+
+ //Vector we use to count the online thresholds
+ m_onlineThresholdFailureCount.resize(m_maxStripHash+1,0);
+
+ //Retrieve current online thresholds
+ m_onlineThresholds.resize(m_maxStripHash+1);
+ ifstream ifile; ifile.open(m_onlineDbFile.c_str());
+ if(!ifile.is_open()){
+ mLog << MSG::FATAL << "Failed to open online database file " << m_onlineDbFile << endreq;
+ return StatusCode::FAILURE;
+ }
+ std::string buf;
+ ifile >> buf; // skip 32 at start of file
+ unsigned int onlineId;
+ unsigned int hashId;
+ double rms;
+ double f001;
+
+ if(!ifile){
+ mLog << MSG::FATAL << "Problem with file after one word read in." << endreq;
+ return StatusCode::FAILURE;
+ }
+
+
+ if(m_debug) mLog << MSG::DEBUG << "Reading in online thresholds from file " << m_onlineDbFile << endreq;
+ if(m_debug) mLog << MSG::DEBUG << "First (junk) word: " << buf << endreq;
+ int chanCnt = 0;
+ while(ifile >> hex >> onlineId >> dec) {
+ chanCnt++;
+ onlineToOfflineHashId(onlineId,hashId);
+
+ ifile >> buf >> buf >> buf >> buf >> rms >> buf >> f001;
+ double thold = f001 + 2*rms;
+ if(m_verbose) mLog << MSG::VERBOSE
+ << "onlid: " << hex << onlineId << dec << " hash: " << hashId << " rms: " << rms << " f001: " << f001 << " thold: " << thold
+ << endreq;
+ m_onlineThresholds.at(hashId) = thold;
+ if(m_verbose){
+ if(!ifile)
+ mLog << MSG::VERBOSE << "input file is done, ready to close!"<< std::endl;
+ }
+ else
+ mLog << MSG::VERBOSE << "Input file still good!" <<std::endl;
+
+
+ }
+ if(chanCnt != 30720){
+ mLog << MSG::FATAL << "Did not retrieve expected 30720 channels from online database! Retrieved: " << chanCnt << endreq;
+ mLog << MSG::FATAL << "Last onlineId read: " << hex << onlineId << dec << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ }//if m_doOnlineDBFile
+ }//db file
+
+ mLog << MSG::INFO << "highest strip hash id is " << m_maxStripHash << endreq;
+
+ //If we're doing correlation plots, set up the product histogram array
+ if(m_doCorrelation){
+ m_bitProds = new DataVector<TH2F>(SG::VIEW_ELEMENTS);
+ m_bitProds->resize(m_maxStripHash+1);
+ }
+
+
+
+
+
+ cout << "m_prods value: " << m_bitProds << "\ndoCorrelation" << m_doCorrelation << endl;
+ return StatusCode::SUCCESS;
+ }
+
+ //Execute loops through all strips and fills histograms
+ StatusCode CscCalcPed::execute()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ if(m_debug) mLog << MSG::DEBUG << "Begin execute" << endreq;
+ //collectEventInfo collects infomation about each event by filling ampHistCollection and peaktHist.
+ StatusCode sc = collectEventInfo();
+
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::ERROR << "There was an error collecting information from the RDO this event." << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+ mLog << MSG::INFO << "End execute" << endreq;
+ return StatusCode::SUCCESS;
+ } //end execute()
+
+ StatusCode CscCalcPed::finalize() {
+
+ MsgStream mLog( msgSvc(), name() );
+
+ if(m_eventCnt ==0)
+ {
+ mLog <<MSG::FATAL << "No events processed!" << endreq;
+ return StatusCode::FAILURE;
+ }
+ else
+ mLog << MSG::INFO << "In finalize() after analyzing " << m_eventCnt << " events " <<endreq;
+
+ StatusCode sc;
+
+ mLog << MSG::INFO << "not dump all hists!" << endreq;
+
+ //calculateParameters() finds means and fits gain curves from the data in
+ //m_ampHistCollection and/or m_peaktHist
+ sc =calculateParameters();
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::FATAL << "Calculation of parameters failed!" << endreq;
+ return StatusCode::FAILURE;
+ }
+ if(m_debug) mLog << MSG::DEBUG << "Finished calculating parameters" << endreq;
+
+
+ //writeCalibrationFile() writes the calculated parameters into a calibration fie.
+ sc = writeCalibrationFile();
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::FATAL << "Calculation of parameters failed!" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ //Record report (ampHists) and results to storegate for future algorithms
+ //like COOL database storage or monitoring to collect
+ sc =storeGateRecord();
+ if (!sc.isSuccess())
+ {
+ mLog << MSG::FATAL << "Failed recording data in storegate" <<endreq;
+ }
+
+ mLog << MSG::INFO << "Finished finalize" << endreq;
+ return StatusCode::SUCCESS;
+ }//end finalize()
+
+
+ //Collect event info is run every event, gathering amplitutdes and peakting times.
+ StatusCode CscCalcPed::collectEventInfo()
+ {
+ MsgStream mLog( msgSvc(), name() );
+
+
+ //start the timer
+ Chrono chrono(m_chronoSvc,"collectEventInfo");
+
+
+ m_eventCnt++;
+ if(m_debug) mLog << MSG::DEBUG <<"Collecting event info for event " << m_eventCnt << endreq;
+ //Below might need to be changed depending on how we get data
+ const CscRawDataContainer* rawDataContainter;
+ StatusCode sc_read = m_storeGate->retrieve(rawDataContainter, "CSCRDO");
+ if (sc_read != StatusCode::SUCCESS)
+ {
+ mLog << MSG::FATAL << "Could not find event" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ if(m_verbose) mLog << MSG::VERBOSE <<"Retrieved RDO from storegate " << endreq;
+
+ if(rawDataContainter->size() == 0)
+ {
+ mLog << MSG::FATAL << "no rods in RDO!" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ if(m_verbose) mLog << MSG::VERBOSE <<"There are " << rawDataContainter->size() << " rods in the RDO" << endreq;
+
+ IdContext channelContext = m_cscId->channel_context();
+
+ //Loop over RODs (data from 2 chambers), each of which is in
+ //a single CscRawaData collection
+ CscRawDataContainer::const_iterator rodItr = rawDataContainter->begin();
+ CscRawDataContainer::const_iterator rodEnd = rawDataContainter->end();
+ if(m_verbose) mLog << MSG::VERBOSE <<"does rodItr == rodEnd? " << (rodItr == rodEnd) << endreq;
+ // mLog << MSG::INFO<< " going to read all " << rawDataContainter->size() << " rod's" << endreq;
+ for(;rodItr != rodEnd; rodItr++)
+ {
+ Chrono chronoRod(m_chronoSvc,"RodItr");
+ if(m_verbose) mLog << MSG::VERBOSE <<"Examining a ROD" << endreq;
+ const CscRawDataCollection * rod = (*rodItr); //Removing another "pointer layer" to make
+ // syntax simpler
+
+ if(m_verbose) mLog << MSG::VERBOSE <<"There are " << rod->size() << " clusters in the ROD" << endreq;
+ if(rod->size() >0)
+ {
+ //Loop over strips in rod
+ CscRawDataCollection::const_iterator clusItr = rod->begin();
+ CscRawDataCollection::const_iterator clusEnd = rod->end();
+ for(; clusItr!=clusEnd ; clusItr++)
+ {
+ Chrono chronoClus(m_chronoSvc,"ClusterItr");
+ const CscRawData * cluster = (*clusItr); //Note: For a ped run, the "cluster"
+ // is the size of an entire layer.
+ int numStrips = cluster->width();
+ int samplesPerStrip = (cluster->samples()).size()/numStrips;
+
+ if(m_verbose) mLog << MSG::VERBOSE << "About to collect info from " << numStrips << " strips" << endl;
+ for(int stripItr = 0; stripItr <numStrips; stripItr++)
+ {
+
+ Chrono chronoStrip(m_chronoSvc,"stripItr");
+ // WP Added
+ Identifier channelId =m_cscRdoDecoderTool->channelIdentifier(cluster, stripItr);
+ IdentifierHash cscChannelHashId;
+ m_cscId->get_channel_hash(channelId, cscChannelHashId);
+
+ ///////////
+ int stripHash = cscChannelHashId;
+
+ Identifier stripId;
+ m_cscId->get_id(stripHash, stripId, &channelContext);
+
+
+ Chrono chronoAfterId(m_chronoSvc,"afterID1");
+
+ if( m_cscId->chamberLayer(channelId) != m_expectedChamberLayer)
+ {
+ // cout << "Wrong chamber layer" << endl;
+ mLog << MSG::WARNING << "Wrong chamber layer a hash ("
+ << stripHash << ") from the wrong multilayer has appeared in the data. Its string id is " << m_cscId->show_to_string(stripId)
+ << " " << m_cscId->show_to_string(channelId)
+ << endreq;
+
+ mLog << MSG::INFO << "WP added (1) "
+ << m_cscId->stationEta(stripId) << " " << m_cscId->measuresPhi(stripId) << " "
+ << stripHash << " " << cscChannelHashId
+ << endreq;
+
+ mLog << MSG::INFO << "WP added (2) "
+ << m_cscId->stationEta(stripId) << " " << m_cscId->measuresPhi(stripId) << " "
+ << stripId << " " << channelId
+ << endreq;
+
+
+
+ // if(m_cscId->measuresPhi(stripId))
+ // mLog << MSG::DEBUG <<" bad id Measures Phi" << endreq;
+ // else
+ // mLog << MSG::DEBUG <<" bad id is eta" << endreq;
+ //continue;
+ stripId = m_cscId->channelID(
+ m_cscId->stationName(stripId),
+ m_cscId->stationEta(stripId),
+ m_cscId->stationPhi(stripId),
+ 2,
+ m_cscId->wireLayer(stripId),
+ m_cscId->measuresPhi(stripId),
+ m_cscId->strip(stripId)
+ );
+ IdentifierHash newHash;
+ m_cscId->get_channel_hash(stripId, newHash );
+ stripHash = newHash;
+ if(m_debug) mLog << MSG::DEBUG << "New hash " << stripHash << endreq;
+ }
+ else{
+ if(m_cscId->measuresPhi(stripId))
+ mLog << MSG::VERBOSE <<" good id Measures Phi" << endreq;
+ else
+ mLog << MSG::VERBOSE <<" good id is eta" << endreq;
+ }
+
+ Chrono chronoAfterId2(m_chronoSvc,"afterID2");
+
+ //Get samples. Each shows amplitude of pulse at different
+ //time slice.
+ vector<uint16_t> samples;
+ cluster->samples(stripItr,samplesPerStrip,samples);
+
+ //Test for threshold breach...
+
+ //int minMax = GetMinMax(samples);
+ //if(minMax < (m_cscId->measuresPhi(stripId) ? 50 :30)) {
+ size_t sampCnt = 0;
+ std::vector<uint16_t>::const_iterator sampItr = samples.begin();
+ std::vector<uint16_t>::const_iterator sampEnd = samples.end();
+ for(;sampItr != sampEnd; sampItr++)
+ {
+ //if(m_debug) mLog << MSG::DEBUG << "Filling amplitude histogram" << endreq;
+ (*m_ampHists)[stripHash]->Fill(*sampItr);
+ //cerr << "Filling sample hist " << sampCnt << endl;
+ if(m_doSampleHists)
+ (*((*m_sampHists)[stripHash]))[sampCnt]->Fill(*sampItr);
+ //cerr << "Filled" << sampCnt << endl;
+ if(m_bitHists && sampCnt==1)
+ {
+ TH2F* prodHist = NULL;
+ if(m_bitProds!=NULL)
+ prodHist = (*m_bitProds)[stripHash];
+ if(!fillBitHist((*m_bitHists)[stripHash],*sampItr, prodHist).isSuccess())
+ mLog << MSG::WARNING << "Failed recording bits for strip " << stripHash << endreq;
+ }//end if(m_bitHists)
+
+ if(m_doOnlineDbFile){//m_doF001){
+ //test if any samples are obvoe the online threshold
+ if (*sampItr > m_onlineThresholds[stripHash] ){
+ m_onlineThresholdFailureCount[stripHash]++;
+ if(m_verbose) mLog << MSG::VERBOSE << "StripHash: " << stripHash <<
+ " has online threshold breach. Sample: " << *sampItr << " Thold: "
+ << m_onlineThresholds[stripHash] << endreq;
+ }
+ }
+ sampCnt++;
+ }//end Sample loop
+ // }//end if GetMinMax
+ // else{
+ // mLog << MSG::DEBUG << "Didn't pass GetMinMax() on hash " << stripHash << std::endl;
+ // }
+ }//end strip loop
+ }//end cluster loop
+ }
+ else
+ if(m_debug) mLog << MSG::DEBUG << "There is an empty rod (CscRawDataContainer)." <<endreq;
+ }//end rod loop
+ if(m_debug) mLog << MSG::DEBUG << "end collectEventInfo()" << endreq;
+ return StatusCode::SUCCESS;
+ }// end collectEventInfo()
+
+
+ //Calculate parameters is called during finalize,
+ //calculates the parameter values.
+ StatusCode CscCalcPed::calculateParameters()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ Chrono chrono(m_chronoSvc,"calculateParameters");
+
+
+ //**********Calculate all the parameters from data collected during execute*****************//
+
+ for(unsigned int stripHash = 0 ;stripHash <= m_maxStripHash; stripHash++)
+ {
+ if(stripHash < 50 || stripHash%1000 == 0)
+ {
+ mLog << MSG::INFO << "Analyzing strip with hash " << stripHash << " out of " << m_maxStripHash << endreq;
+ if(m_verbose) mLog <<MSG::VERBOSE << (float)clock()/((float)CLOCKS_PER_SEC) << " is the time" << endreq;
+ }
+
+ TH1I * ampHist = (*m_ampHists)[stripHash];
+ if(ampHist != NULL)
+ {
+ if(m_verbose) mLog << MSG::VERBOSE << "Have data for strip hash " << stripHash << endreq;
+ if(ampHist->GetEntries() >0) //If strip wasn't tested, it won't have entries
+ {
+ //Following Schernau's work
+ float histMean = ampHist->GetMean();
+ float histRMS = ampHist->GetRMS();
+ float histRMSError = ampHist->GetRMSError();
+
+ float lowbound = histMean - 3*histRMS;
+ float highbound = histMean + 3*histRMS;
+ if(m_verbose) mLog << MSG::VERBOSE << "About to fit..." << endreq;
+
+ int result = ampHist->Fit("gaus","QL","",lowbound,highbound);
+ if(m_verbose) mLog << MSG::VERBOSE << "Result is " << result << endreq;
+ TF1 * fittedFunction = ampHist->GetFunction("gaus");
+ //double mean = fittedFunction->GetParameter(1);
+ double meanError = fittedFunction->GetParError(1);
+ double sigma = fittedFunction->GetParameter(2);
+ double sigmaError = fittedFunction->GetParError(2);
+ double chi2 = fittedFunction->GetChisquare();
+ int ndf = fittedFunction->GetNDF();
+
+ m_peds->push_back(new CscCalibResult(stripHash,histMean,meanError,chi2,ndf));
+ m_noises->push_back(new CscCalibResult(stripHash,sigma,sigmaError,chi2,ndf));
+ m_rmses->push_back(new CscCalibResult(stripHash,histRMS,histRMSError,0,0));
+
+
+
+ //Integrated threshold (f001) calculation
+ if(m_doF001){
+ int num = (int)ampHist->GetEntries();
+ int thr = ampHist->GetNbinsX() + 1; // start at overflow bin
+ double maxSum = 0.001*num; // 99.90 of pedestals under thr
+ //double maxSum = 0.0001*num; // 99.99 of pedestals under thr
+
+ double sum = 0;
+ do{
+ sum += ampHist->GetBinContent(thr);
+ thr--;
+ } while ((thr>0)&&(sum<maxSum));
+
+ //double threshold = ampHist->GetXaxis()->GetBinLowEdge(thr) +2; //For some reason +2 here matches Michael Schernau's +1
+ double threshold = ampHist->GetXaxis()->GetBinLowEdge(thr) +1; //For some reason +2 here matches Michael Schernau's +1
+ m_f001s->push_back(new CscCalibResult(stripHash,threshold,0,0,0));
+
+ if(m_doOnlineDbFile){
+ m_onlineTHoldBreaches->push_back(new CscCalibResult(stripHash,m_onlineThresholdFailureCount[stripHash],0));
+ }
+ }
+
+ //Threshold calculation
+ //double threshold = mean + m_thresholdMultiplier*sigma;
+ //double thresholdError = sqrt((meanError*meanError +
+ // m_thresholdMultiplier*sigmaError*m_thresholdMultiplier*sigmaError)/2);
+
+ //if(m_verbose) mLog << MSG::VERBOSE << "===> Ped: " << mean << "\tNoise: " << sigma
+ // << "\tthreshold: " << threshold << endreq;
+
+ }
+ }
+ else
+ if(m_verbose) mLog << MSG::VERBOSE << "Don't have data for strip hash " << stripHash << endreq;
+ }//end loop over strips
+
+
+ //don't need it anymore, clear ram taken by m_failure tests
+ mLog << MSG::DEBUG << "Clearing m_onlineThresholdFailureCount" << endreq;
+ m_onlineThresholdFailureCount.resize(0);
+
+
+ mLog << MSG::INFO << "Completed calculating parameters." << endreq;
+
+ return StatusCode::SUCCESS;
+ }//End calculateParameters()
+
+ //writeCalibrationFile() dumps the parameters to disk
+ StatusCode CscCalcPed::writeCalibrationFile()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ Chrono chrono(m_chronoSvc,"writeFile");
+ //***Take conditions data held in summary histograms and print to the calibration file***//
+ mLog << MSG::INFO << "Parameters calculated, preparing to output to file: " << m_outputFileName << " Types 1 and " << m_calOutputVersion << endreq;
+
+ calOutput1();
+
+ if(m_calOutputVersion == "00-00"){
+ return calOutput0();
+ }
+ else if(m_calOutputVersion == "03-00") {
+ return calOutput3();
+ }
+ else{
+ mLog << "Don't know how to write calibration file version " << m_calOutputVersion << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode CscCalcPed::calOutput0() {
+ MsgStream mLog( msgSvc(), name() );
+
+ ofstream out;
+ out.open(m_outputFileName.c_str());
+ if(!out.is_open())
+ {
+ mLog << MSG::ERROR << "Can't open file " << m_outputFileName.c_str() << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+
+ out << "00-00 ";
+ out << m_peds->size() << " ";
+ out << "ped ";
+ out << "noise ";
+ cout << "rms ";
+ //out << "thold ";
+ out << "END_HEADER\n";
+
+ if(m_debug) mLog << MSG::DEBUG << "Begining loop over all " << m_peds->size()
+ << " channels data was collected for." << endreq;
+
+ //form is:hashID chamber LayerOrientationStrip parametervalue parametervalue
+ CscCalibResultCollection::iterator pedItr = m_peds->begin();
+ CscCalibResultCollection::iterator pedEnd = m_peds->end();
+ CscCalibResultCollection::iterator noiseItr = m_noises->begin();
+ CscCalibResultCollection::iterator rmsItr = m_rmses->begin();
+ for(;pedItr!= pedEnd;pedItr++,noiseItr++, rmsItr++)//,tholdItr++)
+ {
+ int hashId = (*pedItr)->hashId();
+ double ped = (*pedItr)->value();
+ double noise = (*noiseItr)->value();
+ double rms = (*rmsItr)->value();
+ //double thold = (*tholdItr)->value();
+
+ if(m_debug) mLog << MSG::DEBUG << "we're on hash " << hashId << " with pedestal " << ped
+ << "and noise " << noise << endreq;//<< " and threshold " << thold << endreq;
+ Identifier id;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id(hashId,id, &channelContext);
+
+ Identifier chamberId = m_cscId->elementID(id);
+ if(!m_cscId->valid(chamberId))
+ {
+ mLog << MSG::WARNING << chamberId.getString() << " is not a valid id!" << endreq;
+ mLog << MSG::WARNING << "identifier is: " << m_cscId->show_to_string(chamberId) << endreq;
+ //id.show();
+ //in.ignore(10000,'\n');
+ }
+
+ IdentifierHash chamberHash;
+ m_cscId->get_module_hash(id,chamberHash);
+
+ //print out values.
+ out << hashId;
+ out <<" " << chamberHash;
+ out << " " << m_cscId->show_to_string(id) << " ";
+ out << " " << ped;
+ out << " " << noise;
+ out << " " << rms;
+ //out << " " << thold;
+ out << "\n" ;
+ } //end loop over hash Ids
+
+ out.close(); //done writing
+ mLog <<MSG::INFO <<"File written" << endreq;
+ return StatusCode::SUCCESS;
+ }//end calOutput0
+
+ //calOutput1 prints out calibration output file in format Michael Schernau's
+ //online software likes to read
+ StatusCode CscCalcPed::calOutput1() {
+ MsgStream mLog( msgSvc(), name() );
+
+ ofstream out;
+ string onlineFileName = m_outputFileName + "_online";
+
+ out.open(onlineFileName.c_str());
+ if(!out.is_open())
+ {
+ mLog << MSG::ERROR << "Can't open online file " << m_outputFileName.c_str() << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+
+ out << "32\n";
+
+ if(m_debug) mLog << MSG::DEBUG << "Begining loop over all " << m_peds->size()
+ << " channels data was collected for." << endreq;
+
+ //form is:hashID chamber LayerOrientationStrip parametervalue parametervalue
+ CscCalibResultCollection::iterator pedItr = m_peds->begin();
+ CscCalibResultCollection::iterator pedEnd = m_peds->end();
+ CscCalibResultCollection::iterator noiseItr = m_noises->begin();
+ CscCalibResultCollection::iterator rmsItr = m_rmses->begin();
+ CscCalibResultCollection::iterator f001Itr = m_f001s->begin();
+ for(;pedItr!= pedEnd;pedItr++,noiseItr++, rmsItr++, f001Itr++)//,tholdItr++)
+ {
+ int hashId = (*pedItr)->hashId();
+ double ped = (*pedItr)->value();
+ double noise = (*noiseItr)->value();
+ //double rms = (*rmsItr)->value();
+ double f001 = (*f001Itr)->value();
+
+ //double thold = (*tholdItr)->value();
+ string onlineHexId;
+
+ //Online ids are same as "string ids" used internally in COOL db.
+ m_cscCoolStrSvc->indexToStringId(hashId, "CHANNEL", onlineHexId);
+
+ if(m_debug) mLog << MSG::DEBUG << "we're on hash " << hashId << " with pedestal " << ped
+ << "and noise " << noise << endreq;//<< " and threshold " << thold << endreq;
+ Identifier id;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id(hashId,id, &channelContext);
+
+ Identifier chamberId = m_cscId->elementID(id);
+ if(!m_cscId->valid(chamberId))
+ {
+ mLog << MSG::WARNING << chamberId.getString() << " is not a valid id!" << endreq;
+ mLog << MSG::WARNING << "identifier is: " << m_cscId->show_to_string(chamberId) << endreq;
+ //id.show();
+ //in.ignore(10000,'\n');
+ }
+
+ char orientationChar = (m_cscId->measuresPhi(id) ? 'Y':'X');
+
+
+ IdentifierHash chamberHash;
+ m_cscId->get_module_hash(id,chamberHash);
+
+ //print out values.
+ out.setf(ios::right);//right aligned columns
+ out << setfill('0') << setw(8) << onlineHexId;
+ out <<" "
+ << setw(2) << chamberHash << orientationChar << (m_cscId->wireLayer(id)-1)
+ <<" "
+ << setw(3) << m_cscId->strip(id) -1 << " " ;
+ out.setf(ios::fixed);
+
+
+ out << " " << setprecision(3) << setw(8) << ped << " 0000.00";
+ out << " " << setprecision(3) << setw(8) << noise << " 0000.000";
+ out << " " << f001;
+ out << "\n" ;
+ } //end loop over hash Ids
+
+ out.close(); //done writing
+ mLog <<MSG::INFO <<"File written" << endreq;
+ return StatusCode::SUCCESS;
+ }//end calOutput1
+
+ //calOutput3 outputs version 03-00 calibration files, which are what the most recent version
+ //of CscCoolReadWrite takes for input
+ StatusCode CscCalcPed::calOutput3() {
+ MsgStream mLog( msgSvc(), name() );
+
+ ofstream out;
+ out.open(m_outputFileName.c_str());
+ if(!out.is_open())
+ {
+ mLog << MSG::ERROR << "Can't open output 3 type file " << m_outputFileName.c_str() << " for writing " << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+ out << "03-00 <END_HEADER>";
+
+ outputParameter3(*m_peds, out);
+ outputParameter3(*m_noises,out);
+ outputParameter3(*m_rmses,out);
+ if(m_doF001)
+ outputParameter3(*m_f001s,out);
+ out << "\n<END_FILE>";
+ out.close();
+
+ return StatusCode::SUCCESS;
+ }
+
+
+ //Outputs a single parameter in version 03-00
+ void CscCalcPed::outputParameter3(const CscCalibResultCollection & results, ofstream & out){
+
+ out << "\n";
+ out << "<NEW_PAR> " << results.parName() << "\n";
+ std::string idString;
+
+ CscCalibResultCollection::const_iterator resItr = results.begin();
+ CscCalibResultCollection::const_iterator resEnd = results.end();
+ for(; resItr != resEnd; resItr++){
+ unsigned int hashId = (*resItr)->hashId();
+ double value = (*resItr)->value();
+ std::string idString;
+
+ m_cscCoolStrSvc->indexToStringId(hashId, "CHANNEL", idString);
+
+ out << idString << " " << value << "\n";
+ }
+
+ out << "<END_PAR>" ;
+ }
+
+
+
+ //Record ampHist and calib results to storegate for monitoring and maybe other
+ //programs
+ StatusCode CscCalcPed::storeGateRecord()
+ {
+ MsgStream mLog( msgSvc(), name() );
+
+ StatusCode sc = StatusCode::SUCCESS;
+
+ bool thereIsAnError = false;
+
+ string histKey = "cscPedCalibReport";
+ if(m_debug) mLog << MSG::DEBUG << "Recording pedestal amplitude histograms to TDS with key "
+ << histKey << endreq;
+
+ //CscCalibReport has extraraneous monitoring information
+ CscCalibReportPed * report = new CscCalibReportPed("pedAmps");
+ report->setPedAmpHists(m_ampHists); //report now has ownership of the DataVector
+ report->setBitHists(m_bitHists); //report now has ownership of the DataVector
+ if(m_doSampleHists){
+ report->setSampHists(m_sampHists);
+ }
+ if(m_bitProds){
+ report->setBitCorrelation( makeBitCorrelation());
+ }
+
+ CscCalibReportContainer * repCont = new CscCalibReportContainer(histKey);
+ repCont->push_back(report);
+
+ sc = m_storeGate->record(repCont, histKey);
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::ERROR << "Failed to record CscCalibReportPed to storegate" << endreq;
+ thereIsAnError = true;
+ delete repCont;
+ }
+
+
+ //CscCalibResult contains the actual parameters that we recorded, mostly things that should be entered
+ //into cool
+ string key = "CscCalibResultPed";
+ if(m_debug) mLog << MSG::DEBUG << "Recording calibration results to TDS with key " << key << endreq;
+
+ CscCalibResultContainer * calibResults
+ = new CscCalibResultContainer("CscCalibResultPed");
+ calibResults->push_back(m_peds);
+ calibResults->push_back(m_noises);
+ calibResults->push_back(m_rmses);
+ calibResults->push_back(m_f001s);
+ if(m_doOnlineDbFile)
+ calibResults->push_back(m_onlineTHoldBreaches);
+
+ sc = m_storeGate->record(calibResults,key);
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::ERROR << "Failed to record data to storegate" << endreq;
+ thereIsAnError = true;
+ delete calibResults;
+ }
+
+ if(thereIsAnError)
+ return StatusCode::RECOVERABLE;
+
+ return StatusCode::SUCCESS;
+ }
+
+ //Does bit correlation study. Likely will be dropped
+ DataVector<TH2F> * CscCalcPed::makeBitCorrelation() {
+
+ if(!m_bitProds || !m_bitHists)
+ return NULL;
+
+ DataVector<TH2F> * correlations = new DataVector<TH2F>(SG::VIEW_ELEMENTS);
+ correlations->resize(m_maxStripHash +1);
+
+ for(unsigned int hashItr =0; hashItr <= m_maxStripHash; hashItr++) {
+ IdentifierHash stripHash =hashItr;
+ Identifier stripId;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id(stripHash, stripId, &channelContext);
+
+ int chamLayer = m_cscId->chamberLayer(stripId);
+ if(chamLayer == m_expectedChamberLayer) //Only second chamber layer exists
+ {
+ int stationName = m_cscId->stationName(stripId);
+ //int stationEta = m_cscId->stationEta(stripId);
+ int stationPhi = m_cscId->stationPhi(stripId);
+ int stripNumber = m_cscId->strip(stripId);
+ int wireLayer = m_cscId->wireLayer(stripId);
+ char orientation = m_cscId->measuresPhi(stripId) ? 'Y':'X';
+
+ int sector = 2*stationPhi + 50 - stationName;
+
+ stringstream name;
+ name << "h_bitCorr_sector_" << setfill('0') << setw(2) << sector <<
+ "_lay_" << wireLayer << orientation << "_strip_" << setw(3) << stripNumber;
+ stringstream title;
+ title << "h_bitCorr_sector_" << setfill('0') << setw(2) << sector <<
+ "_lay_" << wireLayer << orientation << "_strip_" << setw(3) << stripNumber;
+
+ TH2F* correlationHist = new TH2F(
+ name.str().c_str(),
+ title.str().c_str(),
+ m_numBits, 0, m_numBits,
+ m_numBits, 0, m_numBits
+ );
+ (*correlations)[hashItr] = correlationHist;
+
+ //each amphis is filled exaclty the number of times the bits are sampled,
+ //so its a good place to get n
+ double n = (*m_ampHists)[hashItr]->GetEntries();
+ TH1I * bitHist = (*m_bitHists)[hashItr];
+ TH2F * bitProds = (*m_bitProds)[hashItr];
+ for(unsigned int bit1 = 1; bit1 <=m_numBits; bit1++){
+ for(unsigned int bit2 = 1; bit2 <=bit1; bit2++){
+
+ /*
+ //Prolay's method
+
+ float xyn = bitProds->GetBinContent(bit1,bit2)/n; //Avg product
+ float x = bitHist->GetBinContent(bit1);
+ float xn = x/n; //avg bit1
+ float y = bitHist->GetBinContent(bit2);
+ float yn = y/n; //avg bit2
+ float sigx = sqrt(x*(1-xn)); //sigma 1
+ float sigy = sqrt(y*(1-yn)); //sigma 2
+
+ float pxy = (xyn - xn*yn)/sigx/sigy;
+
+ correlationHist->SetBinContent(bit1,bit2,pxy);
+ if(bit1!=bit2)
+ correlationHist->SetBinContent(bit2,bit1,pxy);
+ */
+ float xy = bitProds->GetBinContent(bit1,bit2);
+ float x = bitHist->GetBinContent(bit1);
+ float y = bitHist->GetBinContent(bit2);
+
+ float r;
+ float denom = (n*x-x*x)*(n*y-y*y);
+ if(denom <= 0 )
+ r= 0;
+ else
+ r = (n*xy - x*y)/sqrt(denom);
+
+ //Pearson r
+ correlationHist->SetBinContent(bit1,bit2,r);
+ if(bit1!=bit2)
+ correlationHist->SetBinContent(bit2,bit1,r);
+
+
+ }
+ }
+ }
+ }
+ return correlations;
+ }//end makeBitCorrelations
+
+
+ StatusCode CscCalcPed::onlineToOfflineHashId(const unsigned int & onlineId, unsigned int &hashId) const
+ {
+ int stationName = ((onlineId >> 16)&0x1) + 50;
+ int phi = ((onlineId >> 13)&0x7)+1;
+ int eta = ((((onlineId >> 12)&0x1) == 1) ? 1:-1);
+ int chamLay = ((onlineId>>11)&0x1) +1;
+ int wireLay = ((onlineId>>9)&0x3) +1;
+ int measuresPhi = ((onlineId >> 8)&0x1);
+ int strip;
+
+ // Online and offline phi ids are flipped on A wheel
+ if( measuresPhi && eta == 1){
+ strip = 48 - ((onlineId)&0xff) ; //equivalent: 49 -( onlineId&0xff +1)
+ }
+ else {
+ strip = ((onlineId)&0xff) +1;
+ }
+
+ Identifier chanId = m_cscId->channelID(stationName,eta,phi,chamLay,wireLay,measuresPhi,strip);
+
+ IdentifierHash chanHash;
+ m_cscId->get_channel_hash(chanId, chanHash);
+
+ hashId = (unsigned int)chanHash;
+
+ return StatusCode::SUCCESS;
+ }
+
+ }//end namespace MuonCalib
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.h b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.h
new file mode 100644
index 0000000000000..a8def2813b701
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcPed.h
@@ -0,0 +1,203 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef CSCCALCPED_H
+#define CSCCALCPED_H_
+/**CscCalcPed - algorithm that finds the Cathode Strip Chamber pedestals from
+an RDO
+*/
+#include "GaudiKernel/Algorithm.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "StoreGate/DataHandle.h"
+
+#include "GaudiKernel/ToolHandle.h"
+
+#include "MuonReadoutGeometry/MuonDetectorManager.h"
+#include "DataModel/DataVector.h"
+
+#include <vector>
+#include <bitset>
+
+#include "CscCalibData/CscCalibResultCollection.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH2F.h"
+#include "MuonCondInterface/CscICoolStrSvc.h"
+
+class cscIdHelper;
+class TFile;
+class IdentifierHash;
+//class ICscCalibTool;
+
+namespace Muon {
+ class ICSC_RDO_Decoder;
+}
+
+namespace MuonCalib{
+ /**
+ @class CscCalcPed
+
+ @brief does calibration of the CSC chambers
+
+ @author lampen@physics.arizona.edu
+
+ @section Description
+ CscCalcPed is an algorithm that cycles through calibration events and generates
+ the pedestals. A root file is also generated where the user can view the validity of the constants.
+ */
+
+ class CscCalcPed: public Algorithm
+ {
+ public:
+ CscCalcPed(const std::string& name, ISvcLocator* pSvcLocator);
+ ~CscCalcPed(void);
+
+ /**basic required functions*/
+ StatusCode initialize(void);
+ StatusCode execute(void);
+ StatusCode finalize(void);
+
+ private:
+
+ /***********Private member functions*/
+ /**event loop functions*/
+ StatusCode collectEventInfo();
+ /**Finalize functions*/
+ StatusCode calculateParameters();
+ StatusCode writeCalibrationFile();
+ StatusCode storeGateRecord();
+ DataVector<TH2F> * makeBitCorrelation();
+
+ StatusCode calOutput0();
+ StatusCode calOutput1();
+ StatusCode calOutput3();
+ void outputParameter3(const CscCalibResultCollection & results, ofstream & out);
+
+ /**Utility functions*/
+ StatusCode hashToChamberName(IdentifierHash,std::string);
+ StatusCode fillBitHist(TH1I * bitHist, const uint16_t & val, TH2F* bitProds);
+ template <typename dataType> dataType GetMinMax(std::vector<dataType> & vec) {
+
+ typename std::vector<dataType>::const_iterator itr = vec.begin();
+ typename std::vector<dataType>::const_iterator end = vec.end();
+ if(itr == end)
+ return 0;
+ dataType max =*itr;
+ dataType min =*itr;
+ itr++;
+ for(; itr != end ; ++itr) {
+ if(*itr < min)
+ min = *itr;
+ if(*itr > max)
+ max = *itr;
+ }
+ return max - min;
+ }
+
+ StatusCode onlineToOfflineHashId(const unsigned int & onlineId, unsigned int &hashId) const;
+
+ /*********Private member variables*/
+ /**Services and tools*/
+ StoreGateSvc * m_storeGate;
+ StoreGateSvc* m_detStore;
+
+ // ITHistSvc * m_thistSvc;
+ // ICscCalibTool * m_cscCalibTool;
+ const CscIdHelper *m_cscId;
+ const MuonGM::MuonDetectorManager * m_muon_mgr;
+ IChronoStatSvc* m_chronoSvc;
+ ServiceHandle<CscICoolStrSvc> m_cscCoolStrSvc;
+ ToolHandle<Muon::ICSC_RDO_Decoder> m_cscRdoDecoderTool;
+
+
+ /**Parameters input through joboptions*/
+ std::string m_outputFileName;
+ std::string m_titlePrefix, m_titlePostfix;
+ std::string m_calOutputVersion;
+ bool m_doCorrelation;
+
+
+ //bool m_makeHists, m_dumpAllHists;
+ float m_thresholdMultiplier;
+ int m_expectedChamberLayer;
+
+ std::string m_onlineDbFile; //!< filename for file with online database information
+
+ /**Internally global variables*/
+ unsigned int m_maxStripHash;
+
+ unsigned int m_ampHistLowBound, m_ampHistHighBound, m_ampHistNumBins;
+
+ const unsigned int m_numBits;
+
+ bool m_doBitHists;
+ bool m_doSampleHists;
+
+ int m_numSamplesExpected;
+
+ bool m_doF001;
+
+ DataVector<TH1I> * m_ampHists;
+ DataVector< DataVector<TH1I> > * m_sampHists;
+ DataVector<TH1I> * m_bitHists;
+ DataVector<TH2F> * m_bitProds;
+ DataVector<TH1F> * m_bitCorrelation;
+ std::vector<int> m_onlineThresholds;
+ std::vector<int> m_onlineThresholdFailureCount;
+ CscCalibResultCollection * m_peds;
+ CscCalibResultCollection * m_noises;
+ CscCalibResultCollection * m_rmses;
+ CscCalibResultCollection * m_f001s;
+ CscCalibResultCollection * m_onlineTHoldBreaches;
+ int m_eventCnt;
+ double * crossTalkFix;
+
+ //String for interface to patch changes
+ std::string cmt_parameter;
+
+ //debug level
+ bool m_debug, m_verbose;
+
+ //Temp
+ int m_tempCnt;
+
+
+ bool m_doOnlineDbFile;
+
+ };//end class CscCalcPed
+
+
+ inline StatusCode CscCalcPed::fillBitHist(TH1I * bitHist, const uint16_t & val, TH2F* bitProds)
+ {
+ if(!bitHist)
+ return StatusCode::RECOVERABLE;
+
+ //Num bits should always be m_numBits
+ std::bitset<12> bitVal(val);
+
+ for(unsigned int bitIndex = 0; bitIndex < m_numBits; bitIndex++){
+ if(bitVal[bitIndex]){
+ bitHist->Fill(bitIndex);
+ if(bitProds){
+ for(unsigned int bitIndex2 = 0 ; bitIndex2 <= bitIndex ; bitIndex2++) {
+ if(bitVal[bitIndex2]){
+ bitProds->Fill(bitIndex,bitIndex2);
+ if(bitIndex != bitIndex2)
+ bitProds->Fill(bitIndex2,bitIndex);
+ }
+
+ }
+ }
+ }
+
+ }
+
+
+ return StatusCode::SUCCESS;
+ }
+
+}//end namespace:
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.cxx b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.cxx
new file mode 100644
index 0000000000000..5f0bebea72548
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.cxx
@@ -0,0 +1,1124 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "GaudiKernel/Chrono.h"
+
+#include "TF1.h"
+#include "TGraphErrors.h"
+#include "TFile.h"
+#include "TProfile.h"
+
+#include "StoreGate/StoreGate.h"
+
+#include "MuonIdHelpers/CscIdHelper.h"
+
+#include "MuonRDO/CscRawData.h"
+#include "MuonRDO/CscRawDataCollection.h"
+#include "MuonRDO/CscRawDataContainer.h"
+#include "CscCalibTools/ICscCalibTool.h"
+
+#include "MuonCondInterface/CscICoolStrSvc.h"
+#include "MuonCSC_CnvTools/ICSC_RDO_Decoder.h"
+//#include "MuonCondSvc/CscCoolStrSvc.h"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <bitset>
+#include <inttypes.h>
+//#include <math>
+
+#include "CscCalcSlope.h"
+#include "BipolarFit.h"
+
+#include "CscCalibData/CscCalibResultContainer.h"
+#include "CscCalibData/CscCalibReportContainer.h"
+#include "CscCalibData/CscCalibReportSlope.h"
+
+
+using namespace std;
+
+namespace MuonCalib {
+
+ CscCalcSlope::CscCalcSlope(const string& name, ISvcLocator* pSvcLocator) :
+ Algorithm(name,pSvcLocator),
+ m_cscCoolStrSvc("MuonCalib::CscCoolStrSvc",name),
+ m_cscRdoDecoderTool ("Muon::CscRDO_Decoder"),
+ m_outputFileName("output.cal"),
+ m_lastPulserLevel(-999),
+ m_currentAmpProf(NULL),
+ m_ampProfs(NULL),
+ m_pulsedChambers(NULL),
+ m_eventCnt(0),
+ m_peakTimeProf(NULL),
+ m_numBits(12)
+ {
+ declareProperty("OutputFile", m_outputFileName = "");
+ //declareProperty("FitCutoff",m_fitCutoff=0);
+ declareProperty("IgnoreDatabaseError",m_ignoreDatabaseError = false); //Set to true to ignore database errors
+ declareProperty("TitlePrefix",m_titlePrefix = ""); //Prefix appended to title of histograms and graphs
+ declareProperty("TitlePostfix",m_titlePostfix = ""); //Postfix appended to title of histograms and graphs
+
+ //Don't put the min at 0! At least some offset!
+ //declareProperty("FitMin",m_calFuncMin=2);
+ //declareProperty("EdgeChanFitMin", m_calFuncEdgeChanMin = 6);
+ //declareProperty("FitMax",m_calFuncMax=31.5);
+
+ //test parameters
+ declareProperty("DoBipolarFit", m_doBipolarFit = true);
+ declareProperty("DoCrossTalkFix",m_doCrossTalkFix = true);
+
+ //declareProperty("DoPulserLinearityCorrection",m_doPulserFactor = true);
+ //declareProperty("PulserLinearityFile", m_pulserFactorFile = "pulserFactors.txt");
+
+ declareProperty("GetPedFromFile",m_pedFile = false);
+ declareProperty("PedFileName",m_pedFileName = "");
+
+ declareProperty("ExpectedChamberLayer", m_expectedChamberLayer = 2);
+
+ declareProperty("DoLinPlot" , m_doLinPlot = false);
+ declareProperty("CalibFitFunction" , m_calFitFunc = "[0] + [1]*10^(x/-20)");
+ declareProperty("MinDeltaAdc",m_minDeltaAdc = 10, "Minimum change in ADC a calgraph needs to drop for a fit lower bound to be set");
+
+ /*
+ ADC = mC(db) = m*C_max*10^(db/20) + intercept // C_max = maximum charge
+
+ we divide by c_max later to get m (the function as in m_calFitFunc has
+ [1] = max ADC, so we divide max charge to get
+
+ since attenuation is negative
+
+ db = -20 log(ADC+intercept/mC_max)
+
+ db = -20 log(ADC - intercept / mC_max)
+
+ (-20 since V^2 gives power disipation)
+
+ */
+ //declareProperty("CalibFitFunction", m_calFitFunc = "-20*TMath::Log10((x-[0])/[1])");
+ //declareProperty("CalibFitFunction", m_calFitFunc = "-20*TMath::Log10((x-[0])/[1])");
+
+ declareProperty("FindPeakTime", m_findPeakTime = true);
+ declareProperty("DoBitHists", m_doBitHists = true);
+
+ declareProperty("CalOutputVersion", m_calOutputVersion="03-00");
+
+ m_crossTalkFix = new double[24];
+ m_crossTalkFix[0] = 1.0322840930;
+ m_crossTalkFix[1] = 1.0422690324;
+ m_crossTalkFix[2] = 1.0235384586;
+ m_crossTalkFix[3] = 1.0183445962;
+ m_crossTalkFix[4] = 1.0151409212;
+ m_crossTalkFix[5] = 1.0152511102;
+ m_crossTalkFix[6] = 1.0103618910;
+ m_crossTalkFix[7] = 1.0113985580;
+ m_crossTalkFix[8] = 1.0040464232;
+ m_crossTalkFix[9] = 1.0049431193;
+ m_crossTalkFix[10] = 0.9997829589;
+ m_crossTalkFix[11] = 1.0003994005;
+ m_crossTalkFix[12] = 0.9826108255;
+ m_crossTalkFix[13] = 0.9850002836;
+ m_crossTalkFix[14] = 0.9831852065;
+ m_crossTalkFix[15] = 0.9826508145;
+ m_crossTalkFix[16] = 0.9804885017;
+ m_crossTalkFix[17] = 0.9811262196;
+ m_crossTalkFix[18] = 0.9784119832;
+ m_crossTalkFix[19] = 0.9777689757;
+ m_crossTalkFix[20] = 0.9704773978;
+ m_crossTalkFix[21] = 0.9738781078;
+ m_crossTalkFix[22] = 0.9710430303;
+ m_crossTalkFix[23] = 0.9743144079;
+ }
+
+ CscCalcSlope::~CscCalcSlope() {}
+
+ StatusCode CscCalcSlope::initialize()
+ {
+ MsgStream mLog( msgSvc(), name() );
+
+ mLog << MSG::INFO << "CscCalcSlope::initialize() called" << endreq;
+
+ //*******Register services and tools *********/
+ // Store Gate active store
+ StatusCode sc = serviceLocator()->service("StoreGateSvc", m_storeGate);
+ if (sc != StatusCode::SUCCESS )
+ {
+ mLog << MSG::FATAL << " Cannot get StoreGateSvc " << endreq;
+ return StatusCode::FAILURE ;
+ }
+
+cerr << "detstore" << endl;
+ StoreGateSvc* detStore = 0;
+ sc = serviceLocator()->service("DetectorStore", detStore);
+
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::FATAL << " DetectorStore not found" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+cerr << "idhelper" << endl;
+ sc = detStore->retrieve(m_cscId,"CSCIDHELPER");
+ if( sc.isFailure())
+ {
+ mLog << MSG::ERROR << " Cannot retrieve CscIdHelper " << endreq;
+ return sc;
+ }
+
+cerr << "chorno" << endl;
+ sc = service("ChronoStatSvc",m_chronoSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Cannot retrieve ChronoStatSvc!" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ /*sc = service("THistSvc", m_thistSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Cannot retrieve THistSvc!" << endreq;
+ return StatusCode::FAILURE;
+ }*/
+
+cerr << "toolsvc" << endl;
+ IToolSvc* toolSvc;
+ sc = service("ToolSvc",toolSvc);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Unable to retrieve ToolSvc" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+cerr << "calibtool" << endl;
+ sc = toolSvc->retrieveTool("CscCalibTool",m_cscCalibTool);
+ if(sc.isFailure())
+ {
+ mLog << MSG::FATAL << "Unable to retrieve CscCalibTool" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+
+
+cerr << "decoder tool" << endl;
+ if (m_cscRdoDecoderTool.retrieve().isFailure()){
+ mLog << MSG::FATAL << "Cannot retrieve Csc_Rdo_Decoder Tool!" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ mLog << MSG::INFO <<"Finished initializing services. " << endreq;
+cerr << "done init services" << endl;
+ //*****Initialize internal variables and histograms*******/
+ m_ampProfs = new std::map<int, TProfile* >();
+ //m_ampHists = new DataVector<TH1F>();
+ //m_pulserLevels = new DataVector< std::vector<int> >;
+ //m_adcValues = new DataVector< std::vector<double> >;
+ //m_adcErrors = new DataVector< std::vector<double> >;
+
+ //Setup lookup table for pulser levels
+ m_dbLevels.resize(64);
+ for(unsigned int pulserLevel=0; pulserLevel < 64; pulserLevel++)
+ m_dbLevels[pulserLevel] = pulserLevel*.5;
+
+ IdContext channelContext = m_cscId->channel_context();
+
+ if(m_doBitHists) m_bitHists = new DataVector<TH1I>(SG::VIEW_ELEMENTS);
+ //Loop through ids to find out what hash range we're working on, and to
+ //initialize histograms.
+ vector<Identifier> ids = m_cscId->idVector();
+ vector<Identifier>::const_iterator chamItr = ids.begin();
+ vector<Identifier>::const_iterator chamEnd = ids.end();
+ m_maxStripHash = 0;
+ for(; chamItr != chamEnd; chamItr++)
+ {
+ std::vector<Identifier> stripVect;
+ m_cscId->idChannels(*chamItr,stripVect);
+
+ /*
+ vector<vector<float> > chamberDbs;
+ for(unsigned int layItr = 0; layItr <4; layItr++)
+ {
+ vector<float> layerDbs;
+ for(unsigned int pulserLevel=0;pulserLevel< 64;pulserLevel++)
+ {
+ layerDbs.push_back(pulserLevel*.5);
+ }
+ chamberDbs.push_back(layerDbs);
+ }
+ m_dbLevels.push_back(chamberDbs);*/
+
+ vector<Identifier>::const_iterator stripItr = stripVect.begin();
+ vector<Identifier>::const_iterator stripEnd = stripVect.end();
+ for(;stripItr != stripEnd; stripItr++)
+ {
+ IdentifierHash stripHash;
+ m_cscId->get_channel_hash(*stripItr,stripHash);
+
+ if(m_maxStripHash < (unsigned int)stripHash)
+ m_maxStripHash = (unsigned int)stripHash;
+
+ if(m_bitHists)
+ {
+ Identifier id;
+ m_cscId->get_id((IdentifierHash)stripHash,id,&channelContext);
+ int wireLayer = m_cscId->wireLayer(id);
+ char orientation = (m_cscId->measuresPhi(id) ? 'Y':'X');
+
+ int stationName = m_cscId->stationName(id);
+ int stationPhi = m_cscId->stationPhi(id);
+ int stationEta = m_cscId->stationEta(id);
+
+ int stripNumber = m_cscId->strip(id);
+
+ char bitName[200], titleSeed[500];
+ //Bit histogram (for looking for stuck-bits)
+ sprintf(bitName, "bitHist%d",(int)stripHash);
+ sprintf(titleSeed, "Bit histogram for eta %d, sector %d, layer %d%c strip %d",
+ stationEta,(2*stationPhi+50 - stationName),wireLayer,orientation,stripNumber);
+ string title = m_titlePrefix + titleSeed + m_titlePostfix;
+ TH1I* hist = new TH1I(bitName, title.c_str(), m_numBits, 0, m_numBits); //12 bits
+ hist->GetXaxis()->SetTitle("Bit");
+ hist->GetYaxis()->SetTitle("Counts");
+ m_bitHists->push_back(hist);
+ }
+ }
+ }//end chamber loop
+
+ m_fitReturns = new std::vector<float>;
+ m_fitReturns->resize(m_maxStripHash+1,0);
+ /*m_allowedStrips = new bool[m_maxStripHash+1]*/;
+
+ m_calGraphs = new DataVector<TGraphErrors>(SG::VIEW_ELEMENTS);
+ for(unsigned int chanItr =0; chanItr <= m_maxStripHash; chanItr++)
+ {
+ m_calGraphs->push_back(NULL);
+ /*IdentifierHash stripHash =chanItr;
+ Identifier stripId;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id(stripHash, stripId, &channelContext);
+ m_allowedStrips[chanItr] = (m_cscId->chamberLayer(stripId) == m_expectedChamberLayer);*/
+ }
+
+
+ if(m_pedFile)
+ {
+ mLog << MSG::INFO << "Opening pedestal file" << endreq;
+ ifstream in(m_pedFileName.c_str());
+ int stripHash;
+ double ped,noise;//,pedError,noiseError;
+ string buff;
+ in >> buff >> buff >> buff >> buff >> buff ;// skip header
+ m_peds = new float[m_maxStripHash+1];
+ m_noises = new float[m_maxStripHash+1];
+
+ while(!in.eof())
+ {
+ in >> stripHash >> buff >> buff >> ped >> noise;
+ mLog << MSG::INFO << stripHash << "\t" << ped << "\t" << noise << endreq;
+ m_peds[stripHash] = ped;
+ m_noises[stripHash] = noise;
+ }
+ }
+ else
+ {
+ if(m_cscCoolStrSvc.retrieve().isFailure())
+ {
+ mLog << MSG::FATAL << "Unable to retrieve CscCoolStrSvc" << endreq;
+ return StatusCode::FAILURE;
+ }
+ }
+
+ /*for(unsigned int stripItr = 0 ; stripItr < m_maxStripHash+1; stripItr++)
+ {
+ m_pulserLevels->push_back(new std::vector<int> );
+ m_adcValues->push_back(new std::vector<double> );
+ m_adcErrors->push_back(new std::vector<double> );
+ }//End strip loop*/
+ mLog << MSG::INFO << "Counted " << m_maxStripHash +1 << " strips." << endreq;
+
+ //m_fracProfs = new DataVector<DataVector<TProfile> >();
+ //m_fracGraphs = new DataVector<DataVector<TGraph> >();
+ m_slopes = new CscCalibResultCollection("pslope");
+ m_intercepts = new CscCalibResultCollection("pinter");
+
+ if(m_findPeakTime)
+ {
+ m_peakTimeProf = new TProfile("PeakTimes","Peaking Time for each channel",m_maxStripHash+1,
+ 0,m_maxStripHash+1);
+ m_peakTimes = new CscCalibResultCollection("peakt");
+ }
+
+ m_pulsedChambers = new std::set<int>;
+
+
+ mLog <<MSG::DEBUG << "End initialize" << endreq;
+ return StatusCode::SUCCESS;
+ }//end initialize
+
+ //Execute loops through all strips and fills histograms
+ StatusCode CscCalcSlope::execute()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ mLog << MSG::INFO << "Begin execute" << endreq;
+ //collectEventInfo collects infomation about each event by filling ampHistCollection and peaktHist.
+ StatusCode sc = collectEventInfo();
+
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::WARNING << "There was an error collecting information from the RDO this event." << endreq;
+ return sc;
+ }
+ mLog << MSG::INFO << "End execute" << endreq;
+ return StatusCode::SUCCESS;
+ } //end execute()
+
+ StatusCode CscCalcSlope::finalize()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ mLog << MSG::INFO << "In finalize()" << endreq;
+
+ StatusCode sc;
+
+ bool thereIsAFatal=false;
+
+ //makeCalibPoints();
+
+ //calculateParameters() finds means and fits gain curves from the data in
+ ///m_adcValues and/or m_allPeaktsHist
+ sc =calculateParameters();
+ if(sc.isFailure())
+ {
+ mLog << MSG::WARNING << "Calculation of parameters failed!" << endreq;
+ }
+ mLog << MSG::DEBUG << "Finished calculating parameters" << endreq;
+
+ //writeCalibrationFile() writes the calculated parameters into a calibration fie.
+ sc = writeCalibrationFile();
+ if(!sc.isSuccess())
+ {
+ mLog << MSG::FATAL << "Failed to write parameters to disk!" << endreq;
+ thereIsAFatal = true; //Not quiting yet to ensure memory is properly deleted
+ }
+
+ sc = storeGateRecord();
+ if(sc.isFailure())
+ {
+ mLog <<MSG::FATAL << "Failed to record parameters in StoreGate " << endreq;
+ thereIsAFatal = true;
+ }
+
+ //delete m_ampHists;
+ //delete m_pulserLevels;
+ //delete m_adcValues;
+ //delete m_adcErrors;
+
+ delete m_peakTimeProf;
+
+ delete [] m_crossTalkFix;
+ mLog <<MSG::DEBUG << "Finished finalize()" << endl;
+
+ if(thereIsAFatal)
+ return StatusCode::FAILURE;
+
+ return StatusCode::SUCCESS;
+ }//end finalize()
+
+
+ //Collect event info is run every event, gathering amplitutdes and peakting times.
+ StatusCode CscCalcSlope::collectEventInfo()
+ {
+ MsgStream mLog( msgSvc(), name() );
+
+ bool thereIsAnError = false;
+
+ Chrono chrono(m_chronoSvc,"collectEventInfo");
+ mLog << MSG::DEBUG <<"Collecting event info for event " << m_eventCnt << endreq;
+ //Below might need to be changed depending on how we get data
+ const CscRawDataContainer* fullRDO;
+ StatusCode sc_read = m_storeGate->retrieve(fullRDO, "CSCRDO");
+ if (sc_read != StatusCode::SUCCESS)
+ {
+ mLog << MSG::FATAL << "Could not find event" << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ // mLog << MSG::INFO <<"Got raw data " << endreq;
+ //Loop over RODs (data from 2 chambers), each of which is in
+ //a single CscRawaData collection
+ CscRawDataContainer::const_iterator rodItr = fullRDO->begin();
+ CscRawDataContainer::const_iterator rodEnd = fullRDO->end();
+ for(;rodItr != rodEnd; rodItr++)
+ {
+ const CscRawDataCollection * rod = (*rodItr); //Removing another "pointer layer" to make
+ if(rod->size() >0)
+ {
+
+ uint16_t pulsedWireLayer = rod->calLayer();
+
+ int pulserLevel = rod->calAmplitude();
+ mLog << MSG::VERBOSE << "Pulser level is " << pulserLevel << endreq;
+ if( pulserLevel != m_lastPulserLevel)
+ {
+ mLog <<MSG::INFO << "New pulser level found. (" << pulserLevel <<")." << endreq;
+
+ //m_ampHists->clear();
+ map<int,TProfile*>::iterator alreadyExistingProfile = m_ampProfs->find(pulserLevel);
+
+ if(alreadyExistingProfile == m_ampProfs->end())
+ {//No previous profile for this amplitude exists
+
+ mLog << MSG::DEBUG << " creating new amplitude profile" << endreq;
+ stringstream name, title;
+ name << "ampProf_" << pulserLevel;
+ title << m_titlePrefix << "Amplitudes For Pulser Level " << pulserLevel << m_titlePostfix;
+ m_currentAmpProf = new TProfile(name.str().c_str(), title.str().c_str(),
+ m_maxStripHash+1, 0, m_maxStripHash);
+ m_currentAmpProf->GetXaxis()->SetTitle("Channel (Hash Id)");
+ m_currentAmpProf->GetYaxis()->SetTitle("Amplitude (ADC value)");
+ mLog << MSG::DEBUG << "Adding new amplitude profile" << endreq;
+ m_ampProfs->insert(pair<int, TProfile*>( pulserLevel, m_currentAmpProf));
+ }
+ else
+ {
+ mLog << MSG::DEBUG << " using existing amplitude profile" << endreq;
+ m_currentAmpProf = alreadyExistingProfile->second;
+ }
+
+ m_lastPulserLevel = pulserLevel;
+ }
+
+ unsigned int samplingPhase = rod->samplingPhase();
+ uint8_t samplingPeriod = rod->rate(); //sampling period in ns
+
+
+ //Loop over strips in rod
+ CscRawDataCollection::const_iterator clusItr = rod->begin();
+ CscRawDataCollection::const_iterator clusEnd = rod->end();
+ for(; clusItr!=clusEnd ; clusItr++)
+ {
+ //mLog << MSG::INFO<< " converting cluster" << endreq;
+ const CscRawData * cluster = (*clusItr); //Note: For a pulser or ped run, the "cluster"
+ //is the size of an entire layer
+ int numStrips = cluster->width();
+ int samplesPerStrip = (cluster->samples()).size()/numStrips;
+
+ IdContext channelContext = m_cscId->channel_context();
+ /*
+ Identifier startId;
+ m_cscId->get_id(hashStart, startId, &channelContex);
+
+ //check the chamber layer.
+ int chamberLayer = m_cscId->chamberLayer(startId) << endreq;
+ if(chamberLayer != m_expectedChamberLayer)
+ {
+ m_log << MSG::ERROR << "A chamber with a chamber layer of " << chamberLayer << " was found. "
+ << m_expectedChamberLayer << " was expected. Something is probably wrong." << endreq;
+ continue;
+ } */
+
+ for(int stripItr = 0; stripItr <numStrips; stripItr++)
+ {
+ Identifier stripId =m_cscRdoDecoderTool->channelIdentifier(cluster, stripItr);
+ IdentifierHash cscChannelHashId;
+ m_cscId->get_channel_hash(stripId, cscChannelHashId);
+ int stripHash = cscChannelHashId;
+ mLog << MSG::VERBOSE << "The eta of this strip is: " << m_cscId->stationEta(stripId) << endreq;
+
+ int chamberLayer = m_cscId->chamberLayer(stripId);
+ if(chamberLayer != m_expectedChamberLayer)
+ {
+ mLog << MSG::FATAL << "Cluster includes strip in chamber layer "
+ << chamberLayer << ". Only " << m_expectedChamberLayer
+ << " is valid." << endreq;
+ thereIsAnError = true;
+ return StatusCode::FAILURE;
+ }
+
+ int currentWireLayer = m_cscId->wireLayer(stripId) - 1 ;
+ bool isThisLayerPulsed = (pulsedWireLayer >> currentWireLayer)&0x1;
+ if(isThisLayerPulsed)
+ {
+ /*Usefull for debug, but slows things down a surprising amount
+ if(!m_cscId->valid(stripId))
+ {
+ mLog << MSG::ERROR << stripId.getString() << " is not a valid id!" << endreq;
+ }*/
+
+ vector<uint16_t> samples;
+ cluster->samples(stripItr,samplesPerStrip,samples); //Retrieve samples for a single strip
+
+ float ped;
+ float noise;
+ if(m_pedFile)
+ {
+ ped = m_peds[stripHash];
+ noise = m_noises[stripHash];
+ }
+ else{
+
+ if(StatusCode::SUCCESS != m_cscCoolStrSvc->getParameter(ped,"ped",stripHash))
+ {
+ ped = 2054;
+ mLog << (m_ignoreDatabaseError ? MSG::WARNING : MSG::ERROR)
+ << "Failed at getting pedestal from COOL for hash " << stripHash << endreq;
+ if(!m_ignoreDatabaseError)
+ return StatusCode::RECOVERABLE;
+ mLog << MSG::WARNING << "Setting to " << ped << endreq;
+ }
+ else
+ mLog << MSG::VERBOSE << "Got pedestal of " << ped << endreq;
+ if(StatusCode::SUCCESS != m_cscCoolStrSvc->getParameter(
+ noise, "noise", stripHash))
+ {
+ noise = .001;
+ mLog << (m_ignoreDatabaseError ? MSG::WARNING : MSG::ERROR) << "Failed at getting noise from COOL for hash " << stripHash << endreq;
+ if(!m_ignoreDatabaseError)
+ return StatusCode::FAILURE;
+ mLog << MSG::WARNING << "Setting to " << noise << endreq;
+ }
+
+ }
+
+ double peakAmp, peakTime;
+ int success;
+ if(!m_doBipolarFit)
+ {
+ //Need to convert vector from ints to floats to pass to findCharge
+ std::vector<float> floatSamples;
+ std::vector<uint16_t>::const_iterator sampItr = samples.begin();
+ std::vector<uint16_t>::const_iterator sampEnd = samples.end();
+ for(;sampItr != sampEnd; sampItr++){
+
+ floatSamples.push_back((*sampItr)-ped);
+ if(m_bitHists){
+ if(!fillBitHist((*m_bitHists)[stripHash],*sampItr)){
+ mLog << MSG::WARNING << "Failed recording bits for strip " << stripHash << endreq;
+ }
+
+ }
+ }
+
+ success = m_cscCalibTool->findCharge((float)samplingPeriod, samplingPhase,floatSamples,peakAmp,peakTime);
+
+ }
+ else
+ {
+ //Need to convert vector from ints to doubles to pass to bipolar fit
+ double adcSamples[4];
+ for(int i = 0; i < 4; i++) adcSamples[i] = samples[i] -ped;
+ double fitResult[3],fitErrors[3], chi2;
+ double width = samplingPeriod == 50 ? 7.2:14.4; //check if 50 or 25ns period
+
+ m_bipolarFit.Fit(adcSamples,noise,ped,width,fitResult,fitErrors,&chi2);
+ success = true;
+ peakAmp = fitResult[0];
+ peakTime = fitResult[1] - (samplingPhase ? 25 : 0);
+ }//end if m_doBipolarFit
+
+
+ if(success)
+ {
+ m_currentAmpProf->Fill(stripHash,peakAmp);
+ //((*m_ampHists)[stripHash])->Fill(peakAmp);
+
+ if(m_findPeakTime)
+ m_peakTimeProf->Fill(stripHash,peakTime);
+ }
+ else
+ {
+ mLog << MSG::WARNING << "Failed at fitting pulse shape. Debug info: " <<endreq;
+ mLog << MSG::WARNING << "stripHash " << stripHash << endreq;
+ mLog << MSG::WARNING << "strip in chamber " << stripItr << endreq;
+ mLog << MSG::WARNING
+ << " and detailed id " << m_cscId->show_to_string(stripId,&channelContext)
+ << endreq;
+ mLog << "Pulsed layer " << pulsedWireLayer <<endreq;
+ mLog << ", Samples: " << samples[0] <<", " << samples[1] << ", "
+ << samples[2] << ", " << samples[3] << endreq;
+ }
+ }//end if (islayerPulsedand and is precision layer)
+ }//end strip loop
+
+ }//end cluster loop
+ }//end if rod >1
+ }//end rod loop
+
+
+ mLog << MSG::DEBUG << "end collectEventInfo()" << endreq;
+ m_eventCnt++;
+
+ if(thereIsAnError)
+ return StatusCode::RECOVERABLE;
+
+ return StatusCode::SUCCESS;
+ }// end collectEventInfo()
+
+
+ //Calculate parameters is called during finalize, and calculates the parameter values from the
+ //data taken by collectEventData()
+ StatusCode CscCalcSlope::calculateParameters()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ Chrono chrono(m_chronoSvc,"calculateParameters");
+ StatusCode sc;
+ mLog << MSG::INFO << "Calculating calibration constants." << endreq;
+
+ unsigned int numCalibPoints = m_ampProfs->size();
+ mLog << MSG::INFO << "There are " << numCalibPoints << " pulser levels to evaluate." << endreq;
+
+
+
+
+ IdContext channelContext = m_cscId->channel_context();
+
+ float chargeMax = 530.88; //in fC
+
+ int crossTalkCnt = 0;
+
+ for(unsigned int stripHash = 0 ;stripHash <= m_maxStripHash; stripHash++)
+ {
+
+ if(true)//stripHash < 50 || stripHash%1000 == 0)
+ {
+ mLog << MSG::INFO << "Analyzing strip with hash " << stripHash << " out of " << m_maxStripHash << endreq;
+ //mLog <<MSG::DEBUG << (float)clock()/((float)CLOCKS_PER_SEC) << " is the time " << endreq;
+ }
+
+ //**Now tackle slope calculation
+
+
+ Identifier id;
+ m_cscId->get_id((IdentifierHash)stripHash,id,&channelContext);
+ int chamberLayer = m_cscId->chamberLayer(id);
+ char orientation = (m_cscId->measuresPhi(id) ? 'Y':'X');
+
+ int wireLayer = m_cscId->wireLayer(id);
+
+
+ int stationName = m_cscId->stationName(id);
+ int stationPhi = m_cscId->stationPhi(id);
+ int stationEta = m_cscId->stationEta(id);
+ int stripNumber = m_cscId->strip(id);
+
+
+ //Decide if we're fitting to an edge strip or not, which
+ //use different TF1s
+ /*
+ TF1 * myFunc = NULL;
+ if(stripNumber <= 2 || stripNumber >= 191)
+ myFunc = edgeFunc;
+ else
+ myFunc = simpleFunc;
+ */
+
+
+ IdentifierHash chamHash;
+ m_cscId->get_module_hash(id,chamHash);
+
+ if(chamberLayer != m_expectedChamberLayer)
+ continue;
+
+ if(m_findPeakTime)
+ {
+ if(m_peakTimeProf->GetBinEntries(stripHash+1)) //See if any peaking times were recorded for strip
+ {
+ float peakt = m_peakTimeProf->GetBinContent(stripHash+1);
+ float peaktError = m_peakTimeProf->GetBinError(stripHash+1);
+ CscCalibResult * peaktResult = new CscCalibResult(stripHash,peakt, peaktError);
+ m_peakTimes->push_back(peaktResult);
+ }
+ }//end if(m_findPeakTime)
+
+ //Don't find slope for this strip if it is a transverse strip
+ if(orientation != 'X')
+ continue;
+
+ //For removing plateau's from fit
+ bool foundMin(false);
+ double fitMinX = 0;
+ double fitMaxX = 0;
+ double lastVal = -1;
+ double lastDrop=0;
+ double thisDrop=0;
+
+
+ TGraphErrors * calGraph = new TGraphErrors(numCalibPoints); //calGraph will be what the gain will be found on
+ char calName[20],titleSeed[500];
+ sprintf(calName, "calGraph%d",stripHash);
+ sprintf(titleSeed, "Calgraph for eta %d, sector %d, layer %d%c, strip %d",stationEta,(2*stationPhi+50 - stationName),wireLayer,orientation, stripNumber);
+ calGraph->SetName(calName);
+ string title = m_titlePrefix + titleSeed + m_titlePostfix;
+ calGraph->SetTitle(title.c_str());
+ calGraph->GetYaxis()->SetTitle("ADC counts");
+ calGraph->GetXaxis()->SetTitle("Attenuation (-db)");
+
+ //
+ mLog << MSG::DEBUG << " Generating " << title << endreq;
+
+ bool isGoodStrip = false;
+
+ //Loop over all attenuation levels, filling the calGraph with the amplitudes
+ //for this strip
+ mLog << MSG::DEBUG << "Number of ampProfs " << m_ampProfs->size() << endreq;
+ if(!m_ampProfs){
+ mLog << MSG::FATAL << "m_ampProfs empty!" << endreq;
+ return StatusCode::FAILURE;
+ }
+ int calPointItr = 0;
+ map<int, TProfile*>::const_iterator ampProfItr = m_ampProfs->begin();
+ map<int, TProfile*>::const_iterator ampProfEnd = m_ampProfs->end();
+ for(; ampProfItr != ampProfEnd; ampProfItr++)
+ {
+ if(!ampProfItr->second){
+ mLog << MSG::FATAL << "Failed at accessing ampProf!" << endreq;
+ return StatusCode::FAILURE;
+ }
+ mLog << MSG::DEBUG << "\tLooking for data for pulser level "
+ << ampProfItr->first << endreq;
+
+ if(ampProfItr->second->GetBinEntries(stripHash+1))
+ {
+
+ mLog << MSG::VERBOSE << "\nHave data for strip " << stripHash<< endreq;
+
+ isGoodStrip = true;
+
+ int pulserLevel = ampProfItr->first;
+ float adcValue = ampProfItr->second->GetBinContent(stripHash+1);
+ float adcError = ampProfItr->second->GetBinError(stripHash+1);
+ if(m_doCrossTalkFix)
+ {
+ mLog <<MSG::VERBOSE << "\tCrosstalk fix " << m_crossTalkFix[crossTalkCnt] << endreq;
+ adcValue /= m_crossTalkFix[crossTalkCnt];
+ adcError /= m_crossTalkFix[crossTalkCnt];
+ }
+ if(adcError != adcError)
+ adcError = 0.01;
+
+ float db = m_dbLevels[pulserLevel];
+
+
+ float attenValue =0;
+ if(m_doLinPlot)
+ attenValue = 300*pow(10,db/20);
+ else
+ attenValue = db;
+
+ mLog << MSG::DEBUG << "\tStoring at db of " << db << " with attenValue " << attenValue << " from pulser level of " << pulserLevel << " and adcValue " << adcValue << endreq;
+
+
+
+ //See if the last two drops were far down enough
+ if(!foundMin){
+ thisDrop = lastVal - adcValue;
+ mLog << MSG::DEBUG << "\tFinding fit min:"
+ << "\tlastVal = " << lastVal
+ << ";lastDrop " << lastDrop << "; thisDrop " << thisDrop << endreq;
+ if(thisDrop > m_minDeltaAdc && lastDrop > m_minDeltaAdc){
+ mLog << MSG::DEBUG << "Found fitMin!" << std::endl;
+ foundMin = true;
+ fitMinX = attenValue;
+ }
+ else{
+ //Not enough deltaADC, store this iterations values for the next loop
+ lastDrop = thisDrop;
+ lastVal = adcValue;
+ }
+ }
+
+ //Use highest attenuation level as fitMaxX
+ if(attenValue > fitMaxX)
+ fitMaxX = attenValue;
+
+ calGraph->SetPoint(calPointItr,attenValue,adcValue);
+ calGraph->SetPointError(calPointItr,0.01,adcError);
+ //calGraph->SetPoint(calPointItr,adcValue,attenValue);
+ //calGraph->SetPointError(calPointItr,adcError,0.01);
+ calPointItr++;
+ //fittedFunction->GetParError(1)); //WARNING: Assumed small charge error
+ }//done if(entries >0)
+
+ }//Done ampProfItr loop
+
+ if(!foundMin && isGoodStrip){
+ mLog << MSG::WARNING << "Failed to find minium for " << title << endreq;
+ }
+
+ //***Do a simple fit to calGraph***
+ //Here we only fit the linear part of the plot. m_fitCutoff can be set by user.
+ if(isGoodStrip)
+ {
+ mLog << MSG::INFO << "we have a good stripHash at " << stripHash << endreq;
+
+ m_pulsedChambers->insert(chamHash); //Programer note: Only gets filled on x-axis. Probably OK.
+
+ float slope, slopeError, intercept, interceptError, chiSquared;
+ int ndf;
+ int fitRet=0;
+
+ //Setup our gain fit function
+ TF1 myFunc("myFunction", m_calFitFunc.c_str(), fitMinX, fitMaxX);
+ myFunc.SetLineColor(kRed);
+ if(m_doLinPlot)
+ {
+ myFunc.SetParameters(0,5);
+ slope = myFunc.GetParameter(1);
+ slopeError = myFunc.GetParError(1);
+ intercept = myFunc.GetParameter(0);
+ interceptError = myFunc.GetParError(0);
+ chiSquared = myFunc.GetChisquare();
+ ndf = myFunc.GetNDF();
+ }
+ else
+ {
+ myFunc.SetParameters(0.1,2000);
+
+ fitRet = calGraph->Fit(&myFunc,"RV");
+
+ slope = myFunc.GetParameter(1)/chargeMax;
+ slopeError = myFunc.GetParError(1);
+ intercept = myFunc.GetParameter(0);
+ interceptError = myFunc.GetParError(0);
+ chiSquared = myFunc.GetChisquare();
+ ndf = myFunc.GetNDF();
+ }
+
+ float invertedSlope;
+ if(abs(slope) < 0.00001 || slope == -999) //watch out for slope==0
+ {
+ mLog << MSG::WARNING << "Slope invalid " << endreq;
+ continue;
+ }
+
+ invertedSlope = 1/slope;
+
+ cerr << "Inserting calgraph in for hash " << stripHash << endl;
+ (*m_calGraphs)[stripHash] = calGraph;
+
+ mLog << MSG::DEBUG << "StripHash: " << stripHash << "; slope: " <<slope
+ << "; intercept: " << intercept
+ << "; chi^2/ndf: " << chiSquared << "/" << ndf << endreq;
+ CscCalibResult * slopeResult = new CscCalibResult(stripHash,invertedSlope,slopeError,chiSquared,ndf);
+ CscCalibResult * interceptResult = new CscCalibResult(stripHash, intercept, interceptError, chiSquared, ndf);
+
+ m_slopes->push_back(slopeResult);
+ m_intercepts->push_back(interceptResult);
+ (*m_fitReturns)[stripHash] = fitRet;
+
+ }//end if(isGoodStrip)
+
+
+ if(crossTalkCnt == 23)
+ crossTalkCnt = 0;
+ else
+ crossTalkCnt++;
+ mLog << MSG::DEBUG << "Looping over next strip..." << endreq;
+ }//end loop over strips
+ mLog << MSG::INFO << "Completed calculating parameters for each strip" << endreq;
+ return StatusCode::SUCCESS;
+ }//End calculateParameters()
+
+ //writeCalibrationFile() dumps the parameters to disk
+ StatusCode CscCalcSlope::writeCalibrationFile()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ Chrono chrono(m_chronoSvc,"writeCalibrationFile");
+ if(m_calOutputVersion == "00-00"){
+ mLog << MSG::INFO << "Printing output file version 00-00" << endreq;
+ return calOutput0();
+ }
+ else if(m_calOutputVersion == "03-00") {
+ mLog << MSG::INFO << "Printing output file version 03-00" << endreq;
+ return calOutput3();
+ }
+ else{
+ mLog << "Don't know how to write calibration file version " << m_calOutputVersion << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+ }
+
+ StatusCode CscCalcSlope::calOutput0(){
+ MsgStream mLog( msgSvc(), name() );
+ //***Take conditions data held in summary histograms and print to the calibration file***//
+ mLog << MSG::INFO << "Parameters calculated, preparing to outputing to file: " << m_outputFileName << endreq;
+ ofstream out;
+ out.open(m_outputFileName.c_str());
+ if(!out.is_open())
+ {
+ mLog << MSG::FATAL << "Can't open file " << m_outputFileName.c_str() << "for writing" << endreq;
+ return StatusCode::FAILURE;
+ }
+ //Start by writing file version number (mainly for COOL program to read)
+ out << "00-00 ";
+ //Number of strips we have info for:
+ out << m_slopes->size() << " ";
+ //print out header
+ out << "pslope ";
+ if(m_findPeakTime) out << "peakt ";
+ out << "END_HEADER\n";
+ //Now we loop over each strip's parameters and print them out
+ mLog << MSG::DEBUG << "Begining loop over all " << m_maxStripHash << " hash ids." << endreq;
+
+ //form is:
+ //hashID chamber LayerOrientationStrip parametervalue parametervalue
+ CscCalibResultCollection::iterator slopeItr = m_slopes->begin();
+ CscCalibResultCollection::iterator slopeEnd = m_slopes->end();
+ CscCalibResultCollection::iterator peaktItr;
+ CscCalibResultCollection::iterator peaktEnd;
+ if(m_findPeakTime)
+ {
+ peaktItr = m_peakTimes->begin();
+ peaktEnd = m_peakTimes->end();
+ }
+ for(; slopeItr != slopeEnd; slopeItr++)
+ {
+ if(m_findPeakTime && (peaktItr == peaktEnd) )
+ {
+ mLog << MSG::FATAL << "Peaktimes out of sync with slopes. Quiting write." << endreq;
+
+ return StatusCode::FAILURE;
+ }
+
+ //Output hashId
+ out << (*slopeItr)->hashId();
+
+ //get id for next few outputs
+ Identifier id;
+ IdContext channelContext = m_cscId->channel_context();
+ m_cscId->get_id((*slopeItr)->hashId(),id, &channelContext);
+
+ //output chamber #
+ IdentifierHash chamberHash;
+ Identifier chamberId = m_cscId->elementID(id);
+ if(!m_cscId->valid(chamberId))
+ {
+ mLog << MSG::FATAL << chamberId.getString() << " is not a valid id!" << endreq;
+ mLog << MSG::FATAL << "identifier is: " << m_cscId->show_to_string(chamberId) << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ m_cscId->get_module_hash(id,chamberHash);
+ out <<" " << chamberHash;
+
+ //output strip details
+ out << " " << m_cscId->show_to_string(id) << " ";
+
+ //output parameter values
+ out << " " << (*slopeItr)->value();
+ if(m_findPeakTime) out << " " << (*peaktItr)->value();
+ out << "\n" ; //to improve readability
+ } //end loop over hash Ids
+
+ out.close(); //done writing
+
+ return StatusCode::SUCCESS;
+ }//end calOutput0
+
+ StatusCode CscCalcSlope::storeGateRecord()
+ {
+ MsgStream mLog( msgSvc(), name() );
+ mLog << MSG::INFO << "Recording csc calibration report." << endreq;
+
+ StatusCode sc = StatusCode::SUCCESS;
+
+ bool thereIsAnError = false;
+
+ string histKey = "cscSlopeCalibReport";
+ mLog <<MSG::DEBUG << "Recording calibration graphs to TDS with key "
+ << histKey << endreq;
+
+ CscCalibReportSlope * report = new CscCalibReportSlope("calGraphs");
+
+ report->setCalGraphs(m_calGraphs);
+ report->setAmpProfs(m_ampProfs);
+ report->setPulsedChambers(m_pulsedChambers);
+ report->setBitHists(m_bitHists);
+ report->setFitResults(m_fitReturns);
+
+ CscCalibReportContainer * repCont = new CscCalibReportContainer(histKey);
+ repCont->push_back(report);
+
+ sc = m_storeGate->record(repCont, histKey);
+ if(sc.isFailure())
+ {
+ mLog << MSG::ERROR << "Failed to record CscCalibReportSlope to storegate" << endreq;
+ thereIsAnError = true;
+ //Since storegate isn't taking ownership, we'll delete it:
+ delete repCont;
+ }
+
+ CscCalibResultContainer * calibResults
+ = new CscCalibResultContainer("CscCalibResultSlope");
+ calibResults->push_back(m_slopes);
+ calibResults->push_back(m_intercepts);
+ if(m_findPeakTime)
+ calibResults->push_back(m_peakTimes);
+ sc = m_storeGate->record(calibResults,"CscCalibResultSlope");
+ if(sc.isFailure())
+ {
+ mLog << MSG::ERROR << "Failed to record results to storegate" << endreq;
+ thereIsAnError = true;
+ //Since storegate isn't taking ownership, we'll delete it
+ delete calibResults;
+ }
+
+ if(thereIsAnError)
+ return StatusCode::RECOVERABLE;
+
+ return StatusCode::SUCCESS;
+ }
+
+
+
+ StatusCode CscCalcSlope::calOutput3() {
+ MsgStream mLog( msgSvc(), name() );
+
+ ofstream out;
+ out.open(m_outputFileName.c_str());
+ if(!out.is_open())
+ {
+ mLog << MSG::ERROR << "Can't open file " << m_outputFileName.c_str() << endreq;
+ return StatusCode::RECOVERABLE;
+ }
+ out << "03-00 <END_HEADER>";
+
+ outputParameter3(*m_slopes, out);
+ //outputParameter3(*m_intercepts, out);
+ //if(m_peakTimes)
+ // outputParameter3(*m_peakTimes, out);
+ out << "\n<END_FILE>";
+ out.close();
+
+ mLog << MSG::INFO << "Successfully opened file " << m_outputFileName << endreq;
+
+ return StatusCode::SUCCESS;
+ }
+
+
+ void CscCalcSlope::outputParameter3(const CscCalibResultCollection & results, ofstream & out){
+ MsgStream mLog( msgSvc(), name() );
+
+ mLog << MSG::INFO << "Printing out parameter " << results.parName() << endreq;
+
+ out << "\n";
+ out << "<NEW_PAR> " << results.parName() << "\n";
+ std::string idString;
+
+ CscCalibResultCollection::const_iterator resItr = results.begin();
+ CscCalibResultCollection::const_iterator resEnd = results.end();
+ for(; resItr != resEnd; resItr++){
+ unsigned int hashId = (*resItr)->hashId();
+ double value = (*resItr)->value();
+ std::string idString;
+
+ m_cscCoolStrSvc->indexToStringId(hashId, "CHANNEL", idString);
+
+ out << idString << " " << value << "\n";
+ }
+
+ out << "<END_PAR>" ;
+ }
+}//end namespace MuonCalib
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.h b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.h
new file mode 100644
index 0000000000000..d120e04d5c3ee
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/CscCalcSlope.h
@@ -0,0 +1,184 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef CSCCALCSLOPE_H
+#define CSCCALCSLOPE_H
+/**CscCalcSlope - algorithm which cycles through all calibration data and
+generates a flat file with the calibration constants. It also generates
+a root file where the user can view the histograms used to find the constants,
+so that he can determine the validity of the constants
+*/
+#include "GaudiKernel/Algorithm.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "StoreGate/DataHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+//#include "MuonCondData/CscCalibData.h"
+#include "DataModel/DataVector.h"
+#include "CscCalibData/CscCalibResultCollection.h"
+#include "MuonCSC_CnvTools/ICSC_RDO_Decoder.h"
+
+//temporary for tests
+#include "BipolarFit.h"
+
+#include <vector>
+#include <map>
+#include <set>
+
+#include "TH1.h"
+
+
+class TProfile;
+class CscICoolStrSvc;
+class ICscCalibTool;
+class cscIdHelper;
+class CscCalibResultCollection;
+class TGraphErrors;
+
+namespace MuonCalib{
+ /**
+ @class CscCalcSlope
+
+ @brief does calibration of the CSC chambers
+
+ @author lampen@physics.arizona.edu
+
+ @section Description
+ CscCalcSlope is an algorithm that cycles through calibration events and generates
+ the calibration constants. A root file is also generated where the user can
+ view the validity of the constants.
+ */
+
+ class CscCalcSlope: public Algorithm
+ {
+ public:
+ CscCalcSlope(const std::string& name, ISvcLocator* pSvcLocator);
+ ~CscCalcSlope(void);
+
+ /**basic required functions*/
+ StatusCode initialize(void);
+ StatusCode execute(void);
+ StatusCode finalize(void);
+
+ private:
+
+ /***********Private member functions*/
+ /**event loop functions*/
+ StatusCode collectEventInfo();
+ /**Finalize functions*/
+ StatusCode calculateParameters();
+ StatusCode writeCalibrationFile();
+ StatusCode storeGateRecord();
+ /**Utility functions*/
+ StatusCode makeCalibPoints();
+ double calShape(double *x, double * par);
+ StatusCode calOutput0();
+ StatusCode calOutput3();
+ void outputParameter3(const CscCalibResultCollection & results, ofstream & out);
+
+ /*********Private member variables*/
+ /**Services and tools*/
+ StoreGateSvc * m_storeGate;
+ //ITHistSvc * m_thistSvc;
+ ICscCalibTool * m_cscCalibTool;
+ ServiceHandle<CscICoolStrSvc> m_cscCoolStrSvc;
+ ToolHandle<Muon::ICSC_RDO_Decoder> m_cscRdoDecoderTool;
+ const CscIdHelper *m_cscId;
+ IChronoStatSvc* m_chronoSvc;
+
+ /**Parameters input through joboptions*/
+ std::string m_outputFileName;
+ std::string m_calOutputVersion;
+
+ //double m_calFuncMin, m_calFuncMax;
+ double m_minDeltaAdc;
+ //double m_calFuncEdgeChanMin;
+
+ bool m_dumpAllHists;
+ bool m_ignoreDatabaseError;
+ bool m_doBitHists;
+
+ std::string m_titlePrefix, m_titlePostfix;
+ /**Internally global variables*/
+ unsigned int m_maxStripHash;
+ int m_lastPulserLevel;
+
+ //DataVector< std::vector<double> > * m_adcValues;
+ //DataVector< std::vector<double> > * m_adcErrors;
+ //DataVector< std::vector<int> > * m_pulserLevels;
+ DataVector<DataVector<TProfile> >* m_fracProfs;
+ DataVector<DataVector<TGraph> >* m_fracGraphs;
+ DataVector<TH1I> * m_bitHists;
+ std::vector<float> * m_fitReturns;
+
+ TGraph * m_resGraph;
+ DataVector<TGraphErrors> * m_calGraphs;
+ //DataVector<TH1F> * m_ampHists;
+ TProfile * m_currentAmpProf;
+ std::map<int, TProfile*> * m_ampProfs;
+ std::set<int> * m_pulsedChambers;
+
+ /**coherent correction array has the corrections to the coherently pulsed channels to get the basic channels**/
+ int m_eventCnt;
+
+ CscCalibResultCollection *m_slopes, *m_intercepts;
+
+
+ //Temporary for testing:
+ BipolarFit m_bipolarFit;
+ bool m_doBipolarFit;
+ double * m_crossTalkFix;
+ bool m_doCrossTalkFix;
+ //bool m_doPulserFactor;
+ std::vector<float> m_dbLevels;
+
+ float *m_peds, *m_noises;
+ //StatusCode onlineToOfflineId(int onlineId, Identifier &chamberId, Identifier &stripId)const;
+ StatusCode fillBitHist(TH1I * bitHist, const uint16_t & val);
+ bool m_pedFile;
+ std::string m_pedFileName;
+ //std::string m_pulserFactorFile;
+
+ int m_expectedChamberLayer;
+ //bool * m_allowedStrips;
+
+ std::string m_calFitFunc;
+
+ //For peaking time
+ bool m_findPeakTime;
+ TProfile * m_peakTimeProf;
+ CscCalibResultCollection * m_peakTimes;
+
+ //allows linear ploting
+ bool m_doLinPlot;
+
+ //String command to allow interface to patched version
+ std::string cmd_parameters;
+
+ //Number of bits
+ unsigned int m_numBits;
+ };//end class CscCalcSlope
+
+ inline StatusCode CscCalcSlope::fillBitHist(TH1I * bitHist, const uint16_t & val)
+ {
+ if(!bitHist)
+ return StatusCode::RECOVERABLE;
+
+ //Num bits should always be m_numBits
+ std::bitset<12> bitVal(val);
+
+ for(unsigned int bitIndex = 0; bitIndex < m_numBits; bitIndex++){
+ if(bitVal[bitIndex]){
+ bitHist->Fill(bitIndex);
+ }
+ }
+
+
+ return StatusCode::SUCCESS;
+ }
+
+}//end namespace:
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_entries.cxx b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_entries.cxx
new file mode 100755
index 0000000000000..3a782a0c6cd1b
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_entries.cxx
@@ -0,0 +1,14 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "../CscCalcPed.h"
+#include "../CscCalcSlope.h"
+
+using namespace MuonCalib;
+
+DECLARE_ALGORITHM_FACTORY( CscCalcPed )
+DECLARE_ALGORITHM_FACTORY( CscCalcSlope )
+
+DECLARE_FACTORY_ENTRIES( CscCalibAlgs )
+{
+ DECLARE_ALGORITHM( CscCalcPed )
+ DECLARE_ALGORITHM( CscCalcSlope )
+}
diff --git a/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_load.cxx b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_load.cxx
new file mode 100755
index 0000000000000..0cbdda671ffcf
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/CscCalib/CscCalibAlgs/src/components/CscCalibAlgs_load.cxx
@@ -0,0 +1,3 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( CscCalibAlgs )
--
GitLab