cvsroot/UserCode/test_16Jul_nav.py

# 15 July changes: moved shelve.close() statements out one loop level. Does it work?
# 15 July ~4pm changes: moved shelve.open() and shelve.update() statement out one level, all shelve statements now out one level
import numpy as n
from ROOT import TFile, TObject, gDirectory, TCanvas, TH1F, TProfile, TList, THashList, TIter, TAttMarker, TAttLine, TDirectory, TGraph, TString
from collections import defaultdict
import math
#import whichdb
import re
import os
from os import path
import glob
import subprocess
import time
import datetime
import shelve
#import cPickle as pickle

outputDictApv1Mean = defaultdict(dict)
titleDict = defaultdict(dict)
outputDictApv1Rms = defaultdict(dict)
outputDictApv2Mean = defaultdict(dict)
outputDictApv2Rms = defaultdict(dict)
fitInputFilesList = []
filePath = '/afs/cern.ch/work/e/eorcutt/private/noise/test15Jul_try2/'
#filePath = ''
#count = 0

#def makeFitOutputFiles():
#    """
#    This function creates the needed output files for the V_depl voltages.
#    There should be only one file per bias scan.
#    """
#    
#    fitOutputFile = TFile("Vdepl "+ str(time.time()) + ".root", "RECREATE")
#    fitOutputFile.Close()
#
##    for module in moduleDictionary:
#
#
def fittingFunction(x, par):
    """
    Piecewise function for fitting, in order to find HV at which the knee exists.
    The function being fit is noise = sqrt( A + B*sqrt( C / V ) ), for
    A : electronics-dependent parameter
    B : electronics-dependent parameter
    C : V_{depl}
    V : high voltage points (x in TF1)
    """

    if x[0] < par[2]:
        value =  par[0] + par[1] * sqrt( par[2] / x[0] )  # won't be a simple polynomial if there exist other noise source (need to take them in quadrature). 
    else:
        value = par[0] + par[1] 

    return value

def getDepletionVoltage(inputProfile):
#    print CurrentDirectory    
##    CurrentDirectory = gDirectory
##    print gDirectory
##    gDirectory.cd()
    tempProfile = TH1F("histo", "Test histo", 21, 0, 400)
##    inputName = inputProfile.GetName()
    profileTitle = str(inputProfile)[24:-17]
#    print profileTitle
    gDirectory.GetObject(profileTitle, tempProfile) 
    fit = TH1("fitIt", fittingFunction, 10, 400, 3)
#
#    # Much thanks to Christian Barth for the parameter ranges. Need to show how to get them myself, though.
    fit.SetParameter(0, 4.)
    fit.SetParameter(2, 100.)
    fit.SetParLimits(2, 40, 400)
    fit.SetParLimits(1, .25, .55)
    fit.SetParName(0, "ElecConst1")
    fit.SetParName(1, "ElecConst2")
    fit.SetParName(2, "V_depl")
#
    tempProfile.Fit(fit, "R")
##    print "Estimated V_depl: %f" % ( tempProfile.GetFunction("fitIt").GetParameter("V_depl") )
    return { tempProfile.GetFunction("fitIt").GetParameter("V_depl"): profileTitle }
#
##This function starts the navigation and analysis process
#def launch():
#    CurrentDirectory = gDirectory
#    DirectoryDigger(CurrentDirectory)
#
##==============================================================================
##This function will recursively traverse the ROOT file until it has reached a directory with no sub directories
##==============================================================================
def DirectoryDigger(CurrentDirectory, outputFile, dir):
#        print "DirectoryDigger run" 
        numberOfSubDirectories = 0
    
        ListOfKeys = CurrentDirectory.GetListOfKeys()
        Iterator = TIter(ListOfKeys)
        NextDirectory = TDirectory()
    
        for next in Iterator:

            LoopObject = next

            if (LoopObject.IsFolder()):
                numberOfSubDirectories += 1
                Title = LoopObject.GetTitle()
                CurrentDirectory.GetObject(Title, NextDirectory)
            
                currentFile = outputFile 
                DirectoryDigger(NextDirectory, currentFile, dir)
        currentFile = outputFile 
        if (numberOfSubDirectories == 0):
#            print "\nWe have reached a bottom directory.\nThe current directory is:\n",CurrentDirectory.pwd() #this was commented out to speed up output.
#            Count -= 1
#            print "There are ", Count, " modules remaining."
            ProfileLooper(CurrentDirectory, currentFile, dir)


def ProfileLooper(CurrentDirectory, currentFile, dir): #This function will identify which histograms to analyze
    Profile = TProfile()
    ListOfKeys = CurrentDirectory.GetListOfKeys()
    Iterator = TIter(ListOfKeys)
   
    
    for next in Iterator:
        Title = next.GetTitle()
        Profile = CurrentDirectory.Get(Title)
        
        if(Profile != None):
            #scanResult=sscanf(myProfile->GetTitle(), "ExpertHisto_Pedestals_FedKey0x%x_LldChannel%d_%s", &fedKey, &laserChannel, plotType);
            ScanExpression = re.compile(r'ExpertHisto_Pedestals_FedKey(?P<fedKey>0[xX][\dA-Fa-f]+)_LldChannel(?P<laserChannel>[-+]?\d+)_(?P<plotType>\S+)')
            MatchObject = ScanExpression.match(Profile.GetTitle())
            
            if(MatchObject != None):
                FedKey = MatchObject.group('fedKey')
                LaserChannel = MatchObject.group('laserChannel')
                PlotType = MatchObject.group('plotType')

                if(PlotType == 'Noise'):   
#                    count++
#                    #print "This histogram is a noise histogram."
#                    print LaserChannel + " " + PlotType
                    
                    profTitle = Profile.GetTitle()

                    outputDictApv1Mean[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv1Mean']
                    outputDictApv1Rms[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv1Rms']
                    outputDictApv2Mean[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv2Mean']
                    outputDictApv2Rms[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv2Rms']
                    titleDict[profTitle][currentFile.GetName()] = ('TProfile: ' + profTitle )

#                    print outputDictApv1Mean[profTitle][currentFile.GetName()]

#### START SHELVE CHUNK ###
                    #startShelveTime = time.time()
#                    print "Shelving..."
                apv1MeanShelf = shelve.open(str(os.path.join(filePath,dir))+"apv1Mean", writeback=True)
#                    apv1MeanShelf = shelve.open("apv1Mean"+dir, writeback=True)
#                    if '"apv1Mean"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv1Mean"+dir')  
                apv1MeanShelf.update(outputDictApv1Mean)
#                    apv2MeanShelf = shelve.open("apv2Mean"+dir, writeback=True)

#                    if '"apv2Mean"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv2Mean"+dir')  
                apv2MeanShelf = shelve.open(str(os.path.join(filePath,dir))+"apv2Mean", writeback=True)
                apv2MeanShelf.update(outputDictApv2Mean)
#                    apv1RmsShelf = shelve.open("apv1Rms"+dir, writeback=True)

#                    if '"apv1Rms"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv1Rms"+dir')  
                apv1RmsShelf = shelve.open(str(os.path.join(filePath,dir))+"apv1Rms", writeback=True)
                apv1RmsShelf.update(outputDictApv1Rms)

#                    apv2RmsShelf = shelve.open("apv2Rms"+dir, writeback=True)
                apv2RmsShelf = shelve.open(str(os.path.join(filePath,dir))+"apv2Rms", writeback=True)
#                    if '"apv2Rms"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv2Rms"+dir')  
                apv2RmsShelf.update(outputDictApv2Rms)

                apv2RmsShelf.close() 
                apv2MeanShelf.close() 
                apv1RmsShelf.close() 
                apv1MeanShelf.close()
#### END SHELVE CHUNK ###
##                    print "\nDone shelving this TProfile.\n"
#
#### NEW SHELVE TEST: ###
#
#                apv1MeanShelf.sync()
#                apv1RmsShelf.sync()
#                apv2MeanShelf.sync()
#                apv2RmsShelf.sync()
                
#                endShelveTime = time.time()
#                print "Sync time: %2f sec" % (endShelveTime - startShelveTime) 
##                    startingTime = time.time()
##                    apv1MeanShelf = shelve.open("apv1Mean")
##                    apv2MeanShelf = shelve.open("apv2Mean")
##                    apv1RmsShelf = shelve.open("apv1Rms")
##                    apv2RmsShelf = shelve.open("apv2Rms")
##                    
##                    apv1Mtemp = apv1MeanShelf[outputDictApv1Mean[Profile.GetTitle()][file1.GetName()]] 
#
#            print "File " + dir + " done." 
#
#### END SHELVE TEST ###
#
#
#
#                    apv1meanPick = open('apv1mean.p','w')
#                    apv2meanPick = open('apv2mean.p','w')
#                    apv1rmsPick = open('apv1rms.p','w')
#                    apv2rmsPick = open('apv2rms.p','w')
#                    pickle.dump(outputDictApv1Mean, apv1meanPick) 
#                    pickle.dump(outputDictApv2Mean, apv2meanPick) 
#                    pickle.dump(outputDictApv1Rms, apv1rmsPick) 
#                    pickle.dump(outputDictApv2Rms, apv2rmsPick) 
#                    endShelveTime = time.time()
#                    print "Time: %3f sec" % (endShelveTime - startShelveTime) 
##                    else:
##                        with open('apv1mean.p','a') as file:
##                            file.write("{0}, {1}\n".format(name, mean))
##
##                        with open('apv2mean.p','a') as file:
##                            file.write("{0}, {1}\n".format(name, mean))
##                        with open('apv1mean.p','a') as file:
##                            file.write("{0}, {1}\n".format(name, mean))
##                        with open('apv1mean.p','a') as file:
##                            file.write("{0}, {1}\n".format(name, mean))
##                    print "%f" % getDepletionVoltage(Profile, CurrentDirectory)
#
## END ProfileLooper() function. ##
#
##==============================================================================

def findNewMeanAndSigma(inputProfile):
    sum = 0
    averageNoise = 0
    allBins = {}
    twoSigmaBins = {}
    binsToRemove = {}

    allBinsApv2 = {}
    twoSigmaBinsApv2 = {}
    binsToRemoveApv2 = {}

    for bin in range(128):
        noise = inputProfile.GetBinContent(bin+1)
        allBins[bin] = inputProfile.GetBinContent(bin+1)
        sum += noise

    averageNoise = sum / 128
    
    for bin in range(128):
        noise = inputProfile.GetBinContent(bin+1)
        if abs(noise - averageNoise) >= 2 * (inputProfile.GetRMS(2)): binsToRemove[bin] = noise

    for key, val in allBins.iteritems():
        if key not in binsToRemove: twoSigmaBins[key] = val

    bins = n.array([float(bin) for bin in allBins.keys()])
    noiseValues = n.array(allBins.values())

    removedValuesProfile = TGraph(len(bins), bins, noiseValues)
    apv1mean = removedValuesProfile.GetMean(2)
    apv1sigma = removedValuesProfile.GetRMS(2)

   #print "Mean (sans outliers): %f" % apv1mean
   #print "Sigma (sans outliers): %f" % apv1sigma

    averageNoise = 0
    sum = 0

    for bin in range(129,256):
       noise = inputProfile.GetBinContent(bin)
       allBinsApv2[bin] = inputProfile.GetBinContent(bin)
       sum += noise

    averageNoise = sum / 128
    inputProfile.GetXaxis().SetRangeUser(129,256)

    for bin in range(129,256):
       noise = inputProfile.GetBinContent(bin+1)
       if abs(noise - averageNoise) >= 2 * (inputProfile.GetRMS(2)): binsToRemoveApv2[bin] = noise

    for key, val in allBinsApv2.iteritems(): twoSigmaBinsApv2[key] = val

    bins = n.array([float(bin) for bin in allBins.keys()])
    noiseValues = n.array(allBinsApv2.values())

    removedValuesProfileApv2 = TGraph(len(bins), bins, noiseValues) # We're storing the extracted and fixed voltages in a TGraph, so we can get statistics on it.
    apv2mean = removedValuesProfileApv2.GetMean(2)
    apv2sigma = removedValuesProfileApv2.GetRMS(2)

    return { 'apv1Mean':apv1mean, 'apv1Rms':apv1sigma, 'apv2Mean':apv2mean, 'apv2Rms':apv2sigma }

#    ### END FindNewMeanAndSigma() function ###

def main():

    #pseudocode
    #loop over date
        #loop over subdetector
            # go through TProfiles, get mean rms(noise)

    eosCommand = "/afs/cern.ch/project/eos/installation/0.1.0-22d/bin/eos.select ls eos/cms/store/user/gbenelli/NoiseBiasScan/"
    commandOutput = subprocess.Popen(eosCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read()
    #dirCommand = listSpecificEosDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read(
#    outputFileKeeperList = []
#    outputFileDict = []
#    subDetectorDict = []

    dirs = [dir for dir in commandOutput.split('\n') if dir]
    inputStartTime = time.time()


    for dir in dirs:#[:1]: #NB: dir is the Jul2011, Dec2010, etc. directories on EOS.
       listSpecificEosDirCommand = eosCommand + dir
       subdirs = [subdir for subdir in subprocess.Popen(listSpecificEosDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read().split('\n') if subdir]
#       print dir 
       outputName = dir + "_" + subdir 
       outputRootFileApv1 = TFile(filePath + str(outputName)+"_apv1.root", "RECREATE")
       outputRootFileApv1.Close()
       outputRootFileApv2 = TFile(filePath + str(outputName)+ "_apv2.root", "RECREATE")
       outputRootFileApv2.Close()

       for subdir in subdirs:#[:1]: #NB: subdir is the subdetector (eg, TECM).
    #       print "Looking into dir %s subdir %s"%(dir,subdir)
           listSpecificEosSubDirCommand = listSpecificEosDirCommand + '/' + subdir
           eosFiles = [file for file in subprocess.Popen(listSpecificEosSubDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read().split('\n') if file]
           print subdir

           for file in eosFiles:
#               print "eosFiles loop run"
               oldFilePath = 'root://eoscms//eos/cms/store/user/gbenelli/NoiseBiasScan/' + dir + '/' + subdir + '/' + file
               outputName = dir + "_" + subdir + ".root"
#               print outputName
    #           makeOutFile = TFile(str(outputName), "RECREATE")
               if path.exists(outputName):

                   currentFile = TFile.Open(oldFilePath)
                   CurrentDirectory = gDirectory
                   print gDirectory
                   DirectoryDigger(CurrentDirectory, currentFile, dir)
#                   print str(currentFile) + " existing"
                   continue
               else:
#                   makeOutFile = TFile(str(outputName), "RECREATE")
#                   makeOutFile.Close()
                   currentFile = TFile.Open(oldFilePath)
                   CurrentDirectory = gDirectory
                   print gDirectory
                   DirectoryDigger(CurrentDirectory, currentFile, dir)
#                   print str(currentFile) + " making"
                
#               DirectoryDigger(CurrentDirectory, currentFile)

               #Here, put the beginning of the output-sorting and analysis stuff.

    #           print makeOutFile.GetName()
    #           print outputName

    #           print subdir
    #           print filePath
    #           file1 = TFile.Open(filePath)
    #           Count = 1
    #           launch(outputFileName)

#            inputEndTime = time.time()
#            outputStartTime = time.time()

#                   outputDictApv1Mean = shelve.open('"Apv1Mean"+dir')

    allowed1 = ['Apv1MeanJul2011', 'Apv1MeanDec2010']            
    for inputFile in allowed1:
        outputDictApv1Mean = shelve.open(filePath + inputFile)
        for key in outputDictApv1Mean: #This (long) loop outputs the rms(noise) points against HV histos into a ROOT file for APV1.
            listOfKeys = outputDictApv1Mean[key].keys()
                
            tempProfile = TH1F(key, key, 425, 0, 425)
                
            for subKey in listOfKeys:
                #This block of code determines the voltage of each ROOT file from the file name
                fileExpression = re.compile(r'(([0-9a-zA-Z/_.:]+[/])?(?P<fileName>[a-zA-Z]+\d+\.root))')
                fileMatchObject = fileExpression.match(subKey)
                
                print "The subkey being matched is ", subKey
            #           
                fileName = fileMatchObject.group('fileName')
            #            
                voltageExpression = re.compile(r'(?P<subDetector>[a-zA-Z]+)(?P<voltage>\d+)\.root')
                voltageMatchObject = voltageExpression.match(fileName)
            #        
                if (voltageMatchObject != None):
                     subDetector = voltageMatchObject.group('subDetector')
                     voltageString = voltageMatchObject.group('voltage')
                     voltage = int(voltageString)
                    
                     outputRootFileApv1 = TFile(filePath + str(outputName) + dir + "_apv1.root", "UPDATE")
                     tempProfile.SetBinContent(voltage, outputDictApv1Mean[key][subKey])
                     tempProfile.SetBinError(voltage, outputDictApv1Rms[key][subKey])
                     tempProfile.SetTitle(inputFile)
                #     getDepletionVoltage(tempProfile)
                #      global profileTitle
                #      profileTitle = titleDict[key][subKey]
                     print "TESTING TITLE: %s" % tempProfile.GetTitle()
                     print "TESTING PROFILE: %f" % tempProfile.GetBinContent(voltage)
            #            
            ##        else:
            ##            print "Something is wrong. The key is: ", key
            #
                tempProfile.Write()
            ##    outputDictApv1Mean[key][subKey]
                outputRootFileApv1.Close()
    #
    outputRootFileApv1.Close()

#    outputDictApv2Mean = shelve.open()

    allowed2 = ['Apv2MeanJul2011', 'Apv2MeanDec2010']            
    for inputFile in allowed2:
        outputDictApv2Mean = shelve.open(filePath + inputFile)

        for key in outputDictApv2Mean: #This (long) loop outputs the rms(noise) points against HV histos into a ROOT file for APV1.
            listOfKeys = outputDictApv2Mean[key].keys()
            tempProfileApv2 = TH1F(key, key, 425, 0, 425)
                
            for subKey in listOfKeys:
                #This block of code determines the voltage of each ROOT file from the file name
                fileExpression = re.compile(r'(([0-9a-zA-Z/_.:]+[/])?(?P<fileName>[a-zA-Z]+\d+\.root))')
                fileMatchObject = fileExpression.match(subKey)
                print "The subkey being matched is ", subKey
                fileName = fileMatchObject.group('fileName')
                voltageExpression = re.compile(r'(?P<subDetector>[a-zA-Z]+)(?P<voltage>\d+)\.root')
                voltageMatchObject = voltageExpression.match(fileName)
            #        
                if (voltageMatchObject != None):
                     subDetector = voltageMatchObject.group('subDetector')
                     voltageString = voltageMatchObject.group('voltage')
                     voltage = int(voltageString)
                     outputRootFileApv2 = TFile(filePath + str(outputName) + dir + +"_apv2" + ".root", "UPDATE")
                     tempProfileApv2.SetBinContent(voltage, outputDictApv1Mean[key][subKey])
                     tempProfileApv2.SetBinError(voltage, outputDictApv1Rms[key][subKey])
                     tempProfileApv2.SetTitle(inputFile)
                     print "TESTING TITLE: %s" % tempProfileApv2.GetTitle()
                     print "TESTING PROFILE: %f" % tempProfileApv2.GetBinContent(voltage)

                tempProfileApv2.Write()
                outputRootFileApv2.Close()
    outputRootFileApv2.Close()
    vdeplOutput = open(filePath + 'vdeplOutput.txt', 'w')

    for inputFile in fitInputFilesList:
        print getDepletionVoltage(inputFile)
        vdeplOutput.write(getDepletionVoltage(inputFile))

        #
#        print "Reading the input took ", ((inputEndTime - inputStartTime)/60), " minutes."
#        print "Outputting the results took ", ((outputEndTime - inputStartTime)/60), " minutes."
#        print "The total time taken was ", (((inputEndTime - inputStartTime)/60) + ((outputEndTime - inputStartTime)/60)), " minutes."
        #
    #END MAIN#

if __name__=="__main__":
    main()


Revision:	1.1
Committed:	Mon Jul 16 09:15:13 2012 UTC (12 years, 9 months ago) by eorcutt
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	JO_21_1, JO_19_12, JO_18_12, t_11_12_JO, JO_18_11, James_1_10, james_afternoon_19_9_2012, logger_myalvac_19Aug12-19h53m43s, James13_8, JO_3Aug, HEAD
Log Message:	Combination of looping and fitting. Runs very slowly.
#	User	Rev	Content
1	eorcutt	1.1	# 15 July changes: moved shelve.close() statements out one loop level. Does it work?
2			# 15 July ~4pm changes: moved shelve.open() and shelve.update() statement out one level, all shelve statements now out one level
3			import numpy as n
4			from ROOT import TFile, TObject, gDirectory, TCanvas, TH1F, TProfile, TList, THashList, TIter, TAttMarker, TAttLine, TDirectory, TGraph, TString
5			from collections import defaultdict
6			import math
7			#import whichdb
8			import re
9			import os
10			from os import path
11			import glob
12			import subprocess
13			import time
14			import datetime
15			import shelve
16			#import cPickle as pickle
17
18			outputDictApv1Mean = defaultdict(dict)
19			titleDict = defaultdict(dict)
20			outputDictApv1Rms = defaultdict(dict)
21			outputDictApv2Mean = defaultdict(dict)
22			outputDictApv2Rms = defaultdict(dict)
23			fitInputFilesList = []
24			filePath = '/afs/cern.ch/work/e/eorcutt/private/noise/test15Jul_try2/'
25			#filePath = ''
26			#count = 0
27
28			#def makeFitOutputFiles():
29			# """
30			# This function creates the needed output files for the V_depl voltages.
31			# There should be only one file per bias scan.
32			# """
33			#
34			# fitOutputFile = TFile("Vdepl "+ str(time.time()) + ".root", "RECREATE")
35			# fitOutputFile.Close()
36			#
37			## for module in moduleDictionary:
38			#
39			#
40			def fittingFunction(x, par):
41			"""
42			Piecewise function for fitting, in order to find HV at which the knee exists.
43			The function being fit is noise = sqrt( A + B*sqrt( C / V ) ), for
44			A : electronics-dependent parameter
45			B : electronics-dependent parameter
46			C : V_{depl}
47			V : high voltage points (x in TF1)
48			"""
49
50			if x[0] < par[2]:
51			value = par[0] + par[1] * sqrt( par[2] / x[0] ) # won't be a simple polynomial if there exist other noise source (need to take them in quadrature).
52			else:
53			value = par[0] + par[1]
54
55			return value
56
57			def getDepletionVoltage(inputProfile):
58			# print CurrentDirectory
59			## CurrentDirectory = gDirectory
60			## print gDirectory
61			## gDirectory.cd()
62			tempProfile = TH1F("histo", "Test histo", 21, 0, 400)
63			## inputName = inputProfile.GetName()
64			profileTitle = str(inputProfile)[24:-17]
65			# print profileTitle
66			gDirectory.GetObject(profileTitle, tempProfile)
67			fit = TH1("fitIt", fittingFunction, 10, 400, 3)
68			#
69			# # Much thanks to Christian Barth for the parameter ranges. Need to show how to get them myself, though.
70			fit.SetParameter(0, 4.)
71			fit.SetParameter(2, 100.)
72			fit.SetParLimits(2, 40, 400)
73			fit.SetParLimits(1, .25, .55)
74			fit.SetParName(0, "ElecConst1")
75			fit.SetParName(1, "ElecConst2")
76			fit.SetParName(2, "V_depl")
77			#
78			tempProfile.Fit(fit, "R")
79			## print "Estimated V_depl: %f" % ( tempProfile.GetFunction("fitIt").GetParameter("V_depl") )
80			return { tempProfile.GetFunction("fitIt").GetParameter("V_depl"): profileTitle }
81			#
82			##This function starts the navigation and analysis process
83			#def launch():
84			# CurrentDirectory = gDirectory
85			# DirectoryDigger(CurrentDirectory)
86			#
87			##==============================================================================
88			##This function will recursively traverse the ROOT file until it has reached a directory with no sub directories
89			##==============================================================================
90			def DirectoryDigger(CurrentDirectory, outputFile, dir):
91			# print "DirectoryDigger run"
92			numberOfSubDirectories = 0
93
94			ListOfKeys = CurrentDirectory.GetListOfKeys()
95			Iterator = TIter(ListOfKeys)
96			NextDirectory = TDirectory()
97
98			for next in Iterator:
99
100			LoopObject = next
101
102			if (LoopObject.IsFolder()):
103			numberOfSubDirectories += 1
104			Title = LoopObject.GetTitle()
105			CurrentDirectory.GetObject(Title, NextDirectory)
106
107			currentFile = outputFile
108			DirectoryDigger(NextDirectory, currentFile, dir)
109			currentFile = outputFile
110			if (numberOfSubDirectories == 0):
111			# print "\nWe have reached a bottom directory.\nThe current directory is:\n",CurrentDirectory.pwd() #this was commented out to speed up output.
112			# Count -= 1
113			# print "There are ", Count, " modules remaining."
114			ProfileLooper(CurrentDirectory, currentFile, dir)
115
116
117			def ProfileLooper(CurrentDirectory, currentFile, dir): #This function will identify which histograms to analyze
118			Profile = TProfile()
119			ListOfKeys = CurrentDirectory.GetListOfKeys()
120			Iterator = TIter(ListOfKeys)
121
122
123			for next in Iterator:
124			Title = next.GetTitle()
125			Profile = CurrentDirectory.Get(Title)
126
127			if(Profile != None):
128			#scanResult=sscanf(myProfile->GetTitle(), "ExpertHisto_Pedestals_FedKey0x%x_LldChannel%d_%s", &fedKey, &laserChannel, plotType);
129			ScanExpression = re.compile(r'ExpertHisto_Pedestals_FedKey(?P<fedKey>0[xX][\dA-Fa-f]+)_LldChannel(?P<laserChannel>[-+]?\d+)_(?P<plotType>\S+)')
130			MatchObject = ScanExpression.match(Profile.GetTitle())
131
132			if(MatchObject != None):
133			FedKey = MatchObject.group('fedKey')
134			LaserChannel = MatchObject.group('laserChannel')
135			PlotType = MatchObject.group('plotType')
136
137			if(PlotType == 'Noise'):
138			# count++
139			# #print "This histogram is a noise histogram."
140			# print LaserChannel + " " + PlotType
141
142			profTitle = Profile.GetTitle()
143
144			outputDictApv1Mean[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv1Mean']
145			outputDictApv1Rms[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv1Rms']
146			outputDictApv2Mean[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv2Mean']
147			outputDictApv2Rms[profTitle][currentFile.GetName()] = findNewMeanAndSigma(Profile)['apv2Rms']
148			titleDict[profTitle][currentFile.GetName()] = ('TProfile: ' + profTitle )
149
150			# print outputDictApv1Mean[profTitle][currentFile.GetName()]
151
152			#### START SHELVE CHUNK ###
153			#startShelveTime = time.time()
154			# print "Shelving..."
155			apv1MeanShelf = shelve.open(str(os.path.join(filePath,dir))+"apv1Mean", writeback=True)
156			# apv1MeanShelf = shelve.open("apv1Mean"+dir, writeback=True)
157			# if '"apv1Mean"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv1Mean"+dir')
158			apv1MeanShelf.update(outputDictApv1Mean)
159			# apv2MeanShelf = shelve.open("apv2Mean"+dir, writeback=True)
160
161			# if '"apv2Mean"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv2Mean"+dir')
162			apv2MeanShelf = shelve.open(str(os.path.join(filePath,dir))+"apv2Mean", writeback=True)
163			apv2MeanShelf.update(outputDictApv2Mean)
164			# apv1RmsShelf = shelve.open("apv1Rms"+dir, writeback=True)
165
166			# if '"apv1Rms"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv1Rms"+dir')
167			apv1RmsShelf = shelve.open(str(os.path.join(filePath,dir))+"apv1Rms", writeback=True)
168			apv1RmsShelf.update(outputDictApv1Rms)
169
170			# apv2RmsShelf = shelve.open("apv2Rms"+dir, writeback=True)
171			apv2RmsShelf = shelve.open(str(os.path.join(filePath,dir))+"apv2Rms", writeback=True)
172			# if '"apv2Rms"+dir' not in fitInputFilesList: fitInputFilesList.append('"apv2Rms"+dir')
173			apv2RmsShelf.update(outputDictApv2Rms)
174
175			apv2RmsShelf.close()
176			apv2MeanShelf.close()
177			apv1RmsShelf.close()
178			apv1MeanShelf.close()
179			#### END SHELVE CHUNK ###
180			## print "\nDone shelving this TProfile.\n"
181			#
182			#### NEW SHELVE TEST: ###
183			#
184			# apv1MeanShelf.sync()
185			# apv1RmsShelf.sync()
186			# apv2MeanShelf.sync()
187			# apv2RmsShelf.sync()
188
189			# endShelveTime = time.time()
190			# print "Sync time: %2f sec" % (endShelveTime - startShelveTime)
191			## startingTime = time.time()
192			## apv1MeanShelf = shelve.open("apv1Mean")
193			## apv2MeanShelf = shelve.open("apv2Mean")
194			## apv1RmsShelf = shelve.open("apv1Rms")
195			## apv2RmsShelf = shelve.open("apv2Rms")
196			##
197			## apv1Mtemp = apv1MeanShelf[outputDictApv1Mean[Profile.GetTitle()][file1.GetName()]]
198			#
199			# print "File " + dir + " done."
200			#
201			#### END SHELVE TEST ###
202			#
203			#
204			#
205			# apv1meanPick = open('apv1mean.p','w')
206			# apv2meanPick = open('apv2mean.p','w')
207			# apv1rmsPick = open('apv1rms.p','w')
208			# apv2rmsPick = open('apv2rms.p','w')
209			# pickle.dump(outputDictApv1Mean, apv1meanPick)
210			# pickle.dump(outputDictApv2Mean, apv2meanPick)
211			# pickle.dump(outputDictApv1Rms, apv1rmsPick)
212			# pickle.dump(outputDictApv2Rms, apv2rmsPick)
213			# endShelveTime = time.time()
214			# print "Time: %3f sec" % (endShelveTime - startShelveTime)
215			## else:
216			## with open('apv1mean.p','a') as file:
217			## file.write("{0}, {1}\n".format(name, mean))
218			##
219			## with open('apv2mean.p','a') as file:
220			## file.write("{0}, {1}\n".format(name, mean))
221			## with open('apv1mean.p','a') as file:
222			## file.write("{0}, {1}\n".format(name, mean))
223			## with open('apv1mean.p','a') as file:
224			## file.write("{0}, {1}\n".format(name, mean))
225			## print "%f" % getDepletionVoltage(Profile, CurrentDirectory)
226			#
227			## END ProfileLooper() function. ##
228			#
229			##==============================================================================
230
231			def findNewMeanAndSigma(inputProfile):
232			sum = 0
233			averageNoise = 0
234			allBins = {}
235			twoSigmaBins = {}
236			binsToRemove = {}
237
238			allBinsApv2 = {}
239			twoSigmaBinsApv2 = {}
240			binsToRemoveApv2 = {}
241
242			for bin in range(128):
243			noise = inputProfile.GetBinContent(bin+1)
244			allBins[bin] = inputProfile.GetBinContent(bin+1)
245			sum += noise
246
247			averageNoise = sum / 128
248
249			for bin in range(128):
250			noise = inputProfile.GetBinContent(bin+1)
251			if abs(noise - averageNoise) >= 2 * (inputProfile.GetRMS(2)): binsToRemove[bin] = noise
252
253			for key, val in allBins.iteritems():
254			if key not in binsToRemove: twoSigmaBins[key] = val
255
256			bins = n.array([float(bin) for bin in allBins.keys()])
257			noiseValues = n.array(allBins.values())
258
259			removedValuesProfile = TGraph(len(bins), bins, noiseValues)
260			apv1mean = removedValuesProfile.GetMean(2)
261			apv1sigma = removedValuesProfile.GetRMS(2)
262
263			#print "Mean (sans outliers): %f" % apv1mean
264			#print "Sigma (sans outliers): %f" % apv1sigma
265
266			averageNoise = 0
267			sum = 0
268
269			for bin in range(129,256):
270			noise = inputProfile.GetBinContent(bin)
271			allBinsApv2[bin] = inputProfile.GetBinContent(bin)
272			sum += noise
273
274			averageNoise = sum / 128
275			inputProfile.GetXaxis().SetRangeUser(129,256)
276
277			for bin in range(129,256):
278			noise = inputProfile.GetBinContent(bin+1)
279			if abs(noise - averageNoise) >= 2 * (inputProfile.GetRMS(2)): binsToRemoveApv2[bin] = noise
280
281			for key, val in allBinsApv2.iteritems(): twoSigmaBinsApv2[key] = val
282
283			bins = n.array([float(bin) for bin in allBins.keys()])
284			noiseValues = n.array(allBinsApv2.values())
285
286			removedValuesProfileApv2 = TGraph(len(bins), bins, noiseValues) # We're storing the extracted and fixed voltages in a TGraph, so we can get statistics on it.
287			apv2mean = removedValuesProfileApv2.GetMean(2)
288			apv2sigma = removedValuesProfileApv2.GetRMS(2)
289
290			return { 'apv1Mean':apv1mean, 'apv1Rms':apv1sigma, 'apv2Mean':apv2mean, 'apv2Rms':apv2sigma }
291
292			# ### END FindNewMeanAndSigma() function ###
293
294			def main():
295
296			#pseudocode
297			#loop over date
298			#loop over subdetector
299			# go through TProfiles, get mean rms(noise)
300
301			eosCommand = "/afs/cern.ch/project/eos/installation/0.1.0-22d/bin/eos.select ls eos/cms/store/user/gbenelli/NoiseBiasScan/"
302			commandOutput = subprocess.Popen(eosCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read()
303			#dirCommand = listSpecificEosDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read(
304			# outputFileKeeperList = []
305			# outputFileDict = []
306			# subDetectorDict = []
307
308			dirs = [dir for dir in commandOutput.split('\n') if dir]
309			inputStartTime = time.time()
310
311
312			for dir in dirs:#[:1]: #NB: dir is the Jul2011, Dec2010, etc. directories on EOS.
313			listSpecificEosDirCommand = eosCommand + dir
314			subdirs = [subdir for subdir in subprocess.Popen(listSpecificEosDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read().split('\n') if subdir]
315			# print dir
316			outputName = dir + "_" + subdir
317			outputRootFileApv1 = TFile(filePath + str(outputName)+"_apv1.root", "RECREATE")
318			outputRootFileApv1.Close()
319			outputRootFileApv2 = TFile(filePath + str(outputName)+ "_apv2.root", "RECREATE")
320			outputRootFileApv2.Close()
321
322			for subdir in subdirs:#[:1]: #NB: subdir is the subdetector (eg, TECM).
323			# print "Looking into dir %s subdir %s"%(dir,subdir)
324			listSpecificEosSubDirCommand = listSpecificEosDirCommand + '/' + subdir
325			eosFiles = [file for file in subprocess.Popen(listSpecificEosSubDirCommand, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT).stdout.read().split('\n') if file]
326			print subdir
327
328			for file in eosFiles:
329			# print "eosFiles loop run"
330			oldFilePath = 'root://eoscms//eos/cms/store/user/gbenelli/NoiseBiasScan/' + dir + '/' + subdir + '/' + file
331			outputName = dir + "_" + subdir + ".root"
332			# print outputName
333			# makeOutFile = TFile(str(outputName), "RECREATE")
334			if path.exists(outputName):
335
336			currentFile = TFile.Open(oldFilePath)
337			CurrentDirectory = gDirectory
338			print gDirectory
339			DirectoryDigger(CurrentDirectory, currentFile, dir)
340			# print str(currentFile) + " existing"
341			continue
342			else:
343			# makeOutFile = TFile(str(outputName), "RECREATE")
344			# makeOutFile.Close()
345			currentFile = TFile.Open(oldFilePath)
346			CurrentDirectory = gDirectory
347			print gDirectory
348			DirectoryDigger(CurrentDirectory, currentFile, dir)
349			# print str(currentFile) + " making"
350
351			# DirectoryDigger(CurrentDirectory, currentFile)
352
353			#Here, put the beginning of the output-sorting and analysis stuff.
354
355			# print makeOutFile.GetName()
356			# print outputName
357
358			# print subdir
359			# print filePath
360			# file1 = TFile.Open(filePath)
361			# Count = 1
362			# launch(outputFileName)
363
364			# inputEndTime = time.time()
365			# outputStartTime = time.time()
366
367			# outputDictApv1Mean = shelve.open('"Apv1Mean"+dir')
368
369			allowed1 = ['Apv1MeanJul2011', 'Apv1MeanDec2010']
370			for inputFile in allowed1:
371			outputDictApv1Mean = shelve.open(filePath + inputFile)
372			for key in outputDictApv1Mean: #This (long) loop outputs the rms(noise) points against HV histos into a ROOT file for APV1.
373			listOfKeys = outputDictApv1Mean[key].keys()
374
375			tempProfile = TH1F(key, key, 425, 0, 425)
376
377			for subKey in listOfKeys:
378			#This block of code determines the voltage of each ROOT file from the file name
379			fileExpression = re.compile(r'(([0-9a-zA-Z/_.:]+[/])?(?P<fileName>[a-zA-Z]+\d+\.root))')
380			fileMatchObject = fileExpression.match(subKey)
381
382			print "The subkey being matched is ", subKey
383			#
384			fileName = fileMatchObject.group('fileName')
385			#
386			voltageExpression = re.compile(r'(?P<subDetector>[a-zA-Z]+)(?P<voltage>\d+)\.root')
387			voltageMatchObject = voltageExpression.match(fileName)
388			#
389			if (voltageMatchObject != None):
390			subDetector = voltageMatchObject.group('subDetector')
391			voltageString = voltageMatchObject.group('voltage')
392			voltage = int(voltageString)
393
394			outputRootFileApv1 = TFile(filePath + str(outputName) + dir + "_apv1.root", "UPDATE")
395			tempProfile.SetBinContent(voltage, outputDictApv1Mean[key][subKey])
396			tempProfile.SetBinError(voltage, outputDictApv1Rms[key][subKey])
397			tempProfile.SetTitle(inputFile)
398			# getDepletionVoltage(tempProfile)
399			# global profileTitle
400			# profileTitle = titleDict[key][subKey]
401			print "TESTING TITLE: %s" % tempProfile.GetTitle()
402			print "TESTING PROFILE: %f" % tempProfile.GetBinContent(voltage)
403			#
404			## else:
405			## print "Something is wrong. The key is: ", key
406			#
407			tempProfile.Write()
408			## outputDictApv1Mean[key][subKey]
409			outputRootFileApv1.Close()
410			#
411			outputRootFileApv1.Close()
412
413			# outputDictApv2Mean = shelve.open()
414
415			allowed2 = ['Apv2MeanJul2011', 'Apv2MeanDec2010']
416			for inputFile in allowed2:
417			outputDictApv2Mean = shelve.open(filePath + inputFile)
418
419			for key in outputDictApv2Mean: #This (long) loop outputs the rms(noise) points against HV histos into a ROOT file for APV1.
420			listOfKeys = outputDictApv2Mean[key].keys()
421			tempProfileApv2 = TH1F(key, key, 425, 0, 425)
422
423			for subKey in listOfKeys:
424			#This block of code determines the voltage of each ROOT file from the file name
425			fileExpression = re.compile(r'(([0-9a-zA-Z/_.:]+[/])?(?P<fileName>[a-zA-Z]+\d+\.root))')
426			fileMatchObject = fileExpression.match(subKey)
427			print "The subkey being matched is ", subKey
428			fileName = fileMatchObject.group('fileName')
429			voltageExpression = re.compile(r'(?P<subDetector>[a-zA-Z]+)(?P<voltage>\d+)\.root')
430			voltageMatchObject = voltageExpression.match(fileName)
431			#
432			if (voltageMatchObject != None):
433			subDetector = voltageMatchObject.group('subDetector')
434			voltageString = voltageMatchObject.group('voltage')
435			voltage = int(voltageString)
436			outputRootFileApv2 = TFile(filePath + str(outputName) + dir + +"_apv2" + ".root", "UPDATE")
437			tempProfileApv2.SetBinContent(voltage, outputDictApv1Mean[key][subKey])
438			tempProfileApv2.SetBinError(voltage, outputDictApv1Rms[key][subKey])
439			tempProfileApv2.SetTitle(inputFile)
440			print "TESTING TITLE: %s" % tempProfileApv2.GetTitle()
441			print "TESTING PROFILE: %f" % tempProfileApv2.GetBinContent(voltage)
442
443			tempProfileApv2.Write()
444			outputRootFileApv2.Close()
445			outputRootFileApv2.Close()
446			vdeplOutput = open(filePath + 'vdeplOutput.txt', 'w')
447
448			for inputFile in fitInputFilesList:
449			print getDepletionVoltage(inputFile)
450			vdeplOutput.write(getDepletionVoltage(inputFile))
451
452			#
453			# print "Reading the input took ", ((inputEndTime - inputStartTime)/60), " minutes."
454			# print "Outputting the results took ", ((outputEndTime - inputStartTime)/60), " minutes."
455			# print "The total time taken was ", (((inputEndTime - inputStartTime)/60) + ((outputEndTime - inputStartTime)/60)), " minutes."
456			#
457			#END MAIN#
458
459			if __name__=="__main__":
460			main()
461
462
463
464
465