diff --git a/Histos.py b/Histos.py index f347c8a..94a5923 100644 --- a/Histos.py +++ b/Histos.py @@ -16,44 +16,42 @@ # Draw Histograms # """ # -# def __init__(self): +class Histos(object): + def __init__(self): -file=ROOT.TFile("histos.root","read") -Gfile=ROOT.TFile("histos_good.root", "read") -histo=file.Get('h_pt') + file=ROOT.TFile("histos.root","read") + Gfile=ROOT.TFile("histos_good.root", "read") -ROOT.TCanvas.__init__._creates = False -canvas = ROOT.TCanvas() + #### bins and bounds????? + #def drawHisto(self, variable, bins, x_min, x_max): + def drawHisto(self, *variable) -ROOT.SetOwnership(canvas, False) + for i in variable: + + if variable="pt": + histo=file.Get('h_pt') + createCanvas(histo) + if variable="eta": + histo=file.Get('h_eta') + + def createCanvas(self, histo): + + ROOT.TCanvas.__init__._creates = False + canvas = ROOT.TCanvas() + + ROOT.SetOwnership(canvas, False) -canvas.cd() -canvas.Range(-68.75, -7.5, 856.25, 42.5) -canvas.SetFillColor(0) -canvas.SetBorderMode(0) -canvas.SetBorderSize(2) -canvas.SetTickx(1) -canvas.SetTicky(1) -canvas.SetLeftMargin(0.15) -canvas.SetRightMargin(0.05) -canvas.SetTopMargin(0.05) -canvas.SetBottomMargin(0.15) -canvas.SetFrameFillStyle(0) -canvas.SetFrameBorderMode(0) -canvas.SetFrameFillStyle(0) -canvas.SetFrameBorderMode(0) -canvas._showGuideLines = False - - -histo.Draw() -plot = {'canvas': canvas} - -canvas.Update() -canvas.Draw() - -canvas.SaveAs("h_pt.png") + canvas.cd() + + histo.Draw() + plot = {'canvas': canvas} + + canvas.Update() + canvas.Draw() + + canvas.SaveAs(histo.name()+".png") ROOT.gApplication.Run() diff --git a/Muon.py b/Muon.py index dd79653..59866b4 100644 --- a/Muon.py +++ b/Muon.py @@ -6,18 +6,28 @@ class Muon(object): - def __init__(self, event, position, pt, eta, energy, vertex): + def __init__(self, event, position, pt,px,py,pz, eta, energy,vertex_z,dB,edB,isolation_sumPt,isolation_emEt,isolation_hadEt,isGlobalMuon,isTrackerMuon,numberOfValidHits,normChi2,charge,Vertex_Z): self.event = event self.position = position self.pt = pt + self.px=px + self.py=py + self.pz=pz self.eta = eta self.energy = energy - self.vertex = vertex - - - def setPt(pt): - self.pt=pt + self.vertex_z = vertex_z + self.dB=dB + self.edB=edB + self.isolation_sumPt=isolation_sumPt + self.isolation_emEt=isolation_emEt + self.isolation_hadEt=isolation_hadEt + self.isGlobalMuon=isGlobalMuon + self.isTrackerMuon=isTrackerMuon + self.numberOfValidHits=numberOfValidHits + self.normChi2=normChi2 + self.charge=charge + self.Vertex_Z=Vertex_Z def setMuon(pt, eta, energy): self.pt=pt @@ -34,12 +44,57 @@ def getPosition (self): def getPt(self): return self.pt + def getPx(self): + return self.px + + def getPy(self): + return self.py + + def getPz(self): + return self.pz + def getEta(self): return self.eta def getEnergy(self): return self.energy + def getvertex_z(self): + return self.vertex_z + + def getdB(self): + return self.dB + + def getedB(self): + return self.edB + + def getIsolation_sumPt(self): + return self.isolation_sumPt + + def getIsolation_emEt(self): + return self.isolation_emEt + + def getIsolation_hadEt(self): + return self.isolation_hadEt + + def getIsGlobalMuon(self): + return self.isGlobalMuon + + def getIsTrackerMuon(self): + return self.isTrackerMuon + + def getNumberOfValidHits(self): + return self.numberOfValidHits + + def getNormChi2(self): + return self.normChi2 + + def getCharge(self): + return self.charge + + def getVertex_Z(self): + return self.Vertex_Z + def printMuon(self): print " Muon del evento: "+ repr(self.getEvent()) + " con posicion: " + repr(self.getPosition()) + " tiene como variables: \n "+ "pt: "+ repr(self.getPt())+ ", eta: "+ repr(self.getEta()) diff --git a/Selector.py b/Selector.py index 52bb9f1..27a5b6e 100644 --- a/Selector.py +++ b/Selector.py @@ -1,4 +1,5 @@ import ROOT as ROOT +from ROOT import * import sys import getopt from DataFormats.FWLite import Events, Handle @@ -6,62 +7,148 @@ from readTTree import readTTree -def main (): +class Selection(object): - # Max and min variables + def __init__(self,pt_min,eta_max,distance,dB_max,isolation,mass_min,normChi2_max,numValidHits ): +#Otros posibles constructores eleigiendo el numero de cortes y en cuales variables."Seian Muchos" + self.pt_min = pt_min + self.eta_max = eta_max + self.distance = distance + self.dB_max = dB_max + self.normChi2_max = normChi2_max + self.isolation = isolation + #dimensionless. (sumPt+emEnergy+hadEnergy)/muon.pt = maxima energia antes de considerarlo como un jet de particulas. + self.mass_min = mass_min + self.numValidHits = numValidHits - pt_min = 12 - eta_max = 2.4 - distance = 0.2 - dB_max = 0.02 # cm. dB=impact parameter - #normChi2_max = 10 - isolation = 0.15 - #dimensionless. (sumPt+emEnergy+hadEnergy)/muon.pt = maxima energia antes de considerarlo como un jet de particulas. - mass_min = 60 - chi2 = 10 - numValidHits = 10 + self.good_muons = [] + self.mass=[] + self.pt_1=[] + self.pt_2=[] - # Good muons list - good_muons = [] - +#####self read y self muons???????########Tambien puedo llamar a esta clase y su metodo desde el script principal(donde este el main, no??? read=readTTree() muons=read.process() - bin=[] - + self.bin=[0]*10 + self.bin[0]=len(muons) + # Histograms for good muons variables - h_pt = ROOT.TH1F('g_pt', 'good_pt', 50, -2, 300) + + self.G_pt = ROOT.TH1F('g_pt', 'Transverse Momentum good Muons', 50, -2, 300) + self.G_eta=ROOT.TH1F( 'g_eta', 'Angle Transverse good Muons', 50, -50 , 50 ) + + for iMuon in range(0, len(muons)): #muons[iMuon].printMuon() - # Good muons list - if selector(muons[iMuon], pt_min): - #muons[iMuon].printMuon() - good_muons.append(muons[iMuon]) - h_pt.Fill(muons[iMuon].getPt()) - # Print elements of good_muons + # The selector function evaluates if the inner muon is GlobalMuon and TrackerMuon returning a boolean type before starting the selection for the different variables cuts. + if selector(muons[iMuon],self.pt_min,self.eta_max,self.distance,self.dB_max,self.isolation,self.normChi2_max,self.numValidHits ): + #append this muon to the good muons list. + self.good_muons.append(muons[iMuon]) + + self.G_pt.Fill(muons[iMuon].getPt()) + self.G_eta.Fill(muons[iMuon].getEta()) + + self.Ghistos = ROOT.TFile("Ghistos.root", "RECREATE") + self.G_pt.Write() + self.G_eta.Write() + self.Ghistos.Close() + + ######??????????meter los histogramas en el bucle de good muos y comentar. + + for i in range (0,len(self.good_muons)): + j=i + self.good_muons[i].getEvent() + while self.good_muons[j+1].getEvent()== self.good_muons[j].getEvent(): + if (self.good_muons[j].getCharge()*self.good_muons[j+1].getCharge())<0: + #get its Lorentz vector through a ROOT function + + tlv1=ROOT.TLorentzVector() + tlv1.SetPxPyPzE(self.good_muons[j].getPx(), self.good_muons[j].getPy(), self.good_muons[j].getPz(), self.good_muons[j].getEnergy()) + tlv2=ROOT.TLorentzVector() + tlv2.SetPxPyPzE(self.good_muons[j+1].getPx(), self.good_muons[j+1].getPy(), self.good_muons[j+1].getPz(), self.good_muons[j+1].getEnergy()) + + #self.pt_1.append(self.good_muons[i].getPt()) + #self.pt_2.append(self.good_muons[i+1].getPt()) + mass=(tlv1+tlv2).M() + if self.mass_min pt_min : + + if muon.getIsGlobalMuon(): + + i=+1 + self.bin.addBin(i) + return True + + if muon.getIsTrackerMuon(): + + i=+1 + self.bin.addBin(i) + return True + + + if muon.getPt() >self.pt_min: + + i=+1 + self.bin.addBin(i) + return True + + if muon.getEta()< self.eta_max: + i=+1 + self.bin.addBin(i) + return True - return False + if muon.getdB()< self.dB_max: + + i=+1 + self.bin.addBin(i) + return True + + + if ((muon.getIsolation_sumPt+getIsolation_emEt()+getIsolation_hadEt)/muon.getPt())self.numValidHits: + + i=+1 + self.bin.addBin(i) + return True + + #Last bin for the efficiency.Just the number of good muons. + i=+1 + self.bin.addBin(i) + return True -def addBin(i): - bin[i]=0 - - -if __name__ == "__main__": - main() diff --git a/prueba.py b/prueba.py new file mode 100644 index 0000000..caa8b70 --- /dev/null +++ b/prueba.py @@ -0,0 +1,43 @@ +import ROOT +from ROOT import * +import numpy as n +import array + +x = array.array("i", [1,2,3,4]) +y = array.array("i",[67,45,33,11]) +n=4 ## Numero de puntos que vas a pintar. Parametro imprescindible de Tgraph. +f1=ROOT.TGraph(4,x,y) +#f1.SetLineColor(2) +#f1.SetLineWidth(1) +f1.SetMarkerColor(46) +f1.SetMarkerStyle(21) +f1.SetMarkerSize(6) +f1.SetTitle("Efficiency") +f1.GetXaxis().SetTitle("Number of Cuts") +f1.GetYaxis().SetTitle("muons") +#f1.Draw("A*H") +#ROOT.gApplication.Run() + +h_eff=ROOT.TH1F("eff","efficiency",4,0,6) +for i in range(0,len(y)): + h_eff.SetBinContent(i,y[i]) + +h_eff.Draw() +ROOT.gApplication.Run() + +##Creo una funcion con parametros. Relleno mi histograma con valores aleatorios de esta. +#fa1=ROOT.TF1("fa1","sin(x)/x",0,10); +#fa1=ROOT.TF1("fa1","TMath::DiLog(x)",0,10) +fa1=ROOT.TF1("fa1","[0]*exp(-0.5*((x-[1])/[2])**2)",0,10) +fa1.SetParameter(0,1) +fa1.SetParameter(1,1) +fa1.SetParameter(2,1) +h=ROOT.TH1F("histo","Histo_Fit_Prueba",50,0,10) +for i in range (0,300): + h.Fill(fa1.GetRandom()) +#Hago un fit con gauss y luego con Breit Wigner +p1=h.Fit("gaus") +#Breit_Wigner=ROOT.TF1("Breit_Wigner",....) +gStyle.SetOptFit() +#h.Draw() +#ROOT.gApplication.Run() diff --git a/readTTree.py b/readTTree.py index b601ec8..7db1b1b 100644 --- a/readTTree.py +++ b/readTTree.py @@ -21,69 +21,125 @@ def __init__(self): # Define and init the variables for each branch as ROOT vectors self.Muon_pt = ROOT.std.vector('float')() + self.Muon_px= ROOT.std.vector('float')() + self.Muon_py= ROOT.std.vector('float')() + self.Muon_pz= ROOT.std.vector('float')() self.Muon_eta = ROOT.std.vector('float')() self.Muon_energy = ROOT.std.vector('float')() self.Muon_vertex_z = ROOT.std.vector('float')() + self.dB = ROOT.std.vector('float')() + self.edB = ROOT.std.vector('float')() + self.Muon_isolation_sumPt = ROOT.std.vector('float')() + self.Muon_isolation_emEt = ROOT.std.vector('float')() + self.Muon_isolation_hadEt = ROOT.std.vector('float')() + self.Muon_isGlobalMuon = ROOT.std.vector('int')() + self.Muon_isTrackerMuon = ROOT.std.vector('int')() + self.Muon_numberOfValidHits = ROOT.std.vector('int')() + self.Muon_normChi2 = ROOT.std.vector('float')() + self.Muon_charge = ROOT.std.vector('int')() + self.Vertex_z = ROOT.std.vector('float')() + # List of all muons as a list of Muon object (class Muon) self.all_muons = [] - self.good_muons = [] - self.h_pt=ROOT.TH1F( 'h_pt', 'Muons Transverse Momentun', 50, -2, 300 ) + # Define and init the histograms for each branch using TH1F from ROOT + self.h_pt=ROOT.TH1F( 'h_pt', 'Muons Transverse Momentun', 50, -2, 20000 ) + self.h_px=ROOT.TH1F( 'h_px', 'Muons x- Momentun', 50, -300, 300 ) + self.h_py=ROOT.TH1F( 'h_py', 'Muons y- Momentun', 50, -300, 300 ) + self.h_pz=ROOT.TH1F( 'h_pz', 'Muons z- Momentun', 50, -300, 300 ) + self.h_eta=ROOT.TH1F( 'h_eta', 'Angle Transvese', 50, -50 , 50 ) + self.h_energy=ROOT.TH1F('h_energy','Muons Energy', 50, -300,300) + self.h_charge=ROOT.TH1F('h_charge','Muons Charge', 50,-2,2) + self.h_normChi2=ROOT.TH1F('h_normChi2', 'Muons Chi2', 50, -200,200) + self.h_numberOfValidHits=ROOT.TH1F('h_numberOfValidHits', 'Number of Valid Hits', 50, -200,200) + self.h_dB=ROOT.TH1F('h_dB','Impact Parameter',50,-1,200) + #self.h_edB=ROOT.TH1F('h_edB','Impact Parameter Error',50,-1,200) Pintar como barras de error en el histograma? + self.h_isolation_sumPt=ROOT.TH1F('h_isolation_sumPt','IsolationX',50, -300,300) + self.h_isolation_emEt=ROOT.TH1F('h_isolation_emEt','IsolationX',50, -300,300) + + + def process(self): - # Tell the tree to populate these given variables when reading an entry. - # First parameter is the branch name and second is the address of the variable where the branch data is placed. + ''' + Initialization of the reading process. This function address the data stored in the tree to the different branches previously defined above. For each variable, the + first parameter in the sentence is the name you want for the branch and the second one + is the name of the real variable where the data is adressed,as we said before. This + variable or branch name create a physical memory direction in your computer where store the information that it contains. + ''' + self.tree.SetBranchAddress("Muon_pt", self.Muon_pt) + self.tree.SetBranchAddress("Muon_px", self.Muon_px) + self.tree.SetBranchAddress("Muon_py", self.Muon_py) + self.tree.SetBranchAddress("Muon_pz", self.Muon_pz) self.tree.SetBranchAddress("Muon_eta", self.Muon_eta) self.tree.SetBranchAddress("Muon_energy", self.Muon_energy) self.tree.SetBranchAddress("Muon_vertex_z", self.Muon_vertex_z) - + self.tree.SetBranchAddress("Muon_dB", self.dB) + self.tree.SetBranchAddress("Muon_edB", self.edB) + self.tree.SetBranchAddress("Muon_isolation_sumPt",self.Muon_isolation_sumPt ) + self.tree.SetBranchAddress("Muon_isolation_emEt",self.Muon_isolation_emEt ) + self.tree.SetBranchAddress("Muon_isolation_hadEt",self.Muon_isolation_hadEt) + self.tree.SetBranchAddress("Muon_isGlobalMuon", self.Muon_isGlobalMuon) + self.tree.SetBranchAddress("Muon_isTrackerMuon", self.Muon_isTrackerMuon) + self.tree.SetBranchAddress("Muon_numberOfValidHits",self.Muon_numberOfValidHits) + self.tree.SetBranchAddress("Muon_normChi2",self.Muon_normChi2) + self.tree.SetBranchAddress("Muon_charge",self.Muon_charge) + self.tree.SetBranchAddress("Muon_Vertex_z",self.Vertex_z) + + + # numEntries: Number of entries(events) of the tree numEntries= self.tree.GetEntries() - #print numEntries - # For each event, populate the tree branches, create every muon and add it to all_muons list + + # For each event or entry,the following loop populates the tree branches, creates every muon and add it to all_muons list for event in range(0, numEntries): - # Populate the tree + # Address the data of each physical variable registed in this event or entry number to its branch associated listed above. self.tree.GetEntry(event) # Loop all muons in each entry = vector size for position in range(0,self.Muon_pt.size()): - muon=Muon(event, position, self.Muon_pt[position], self.Muon_eta[position], self.Muon_energy[position], self.Muon_vertex_z[position]) + muon=Muon(event, position, self.Muon_pt[position], self.Muon_px[position],self.Muon_py[position],self.Muon_pz[position],self.Muon_eta[position], self.Muon_energy[position], self.Muon_vertex_z[position], self.dB[position], self.edB[position],self.Muon_isolation_sumPt[position],self.Muon_isolation_emEt[position],self.Muon_isolation_hadEt[position],self.Muon_isGlobalMuon[position],self.Muon_isTrackerMuon[position], self.Muon_numberOfValidHits[position],self.Muon_normChi2[position],self.Muon_charge[position],self.Vertex_z[position]) # Add each muon to all_muon list self.all_muons.append(muon) + # print muon components #muon.printMuon() + # Fill the histogram for each variable self.h_pt.Fill(self.Muon_pt[position]) + self.h_px.Fill(self.Muon_px[position]) + self.h_py.Fill(self.Muon_py[position]) + self.h_pz.Fill(self.Muon_pz[position]) + self.h_eta.Fill(self.Muon_eta[position]) + self.h_energy.Fill(self.Muon_energy[position]) + self.h_charge.Fill(self.Muon_charge[position]) + self.h_normChi2.Fill(self.Muon_normChi2[position]) + self.h_numberOfValidHits.Fill(self.Muon_numberOfValidHits[position]) + self.h_dB.Fill(self.dB[position]) + self.h_isolation_sumPt.Fill(self.Muon_isolation_sumPt[position]) + self.h_isolation_emEt.Fill(self.Muon_isolation_emEt[position]) + return self.all_muons + # Create a rootfile for the histogramas self.fhistos = ROOT.TFile("histos.root", "RECREATE") # Write histograms in file self.h_pt.Write() + self.h_px.Write() + self.h_py.Write() + self.h_pz.Write() + self.h_eta.Write() + self.h_energy.Write() + self.h_normChi2.Write() + self.h_numberOfValidHits.Write() + self.h_dB.Write() + self.h_isolation_sumPt.Write() + self.h_isolation_emEt.Write() + #Close the file self.fhistos.Close() #return muon list - return self.all_muons - def plotter(self): - """ - Plots the histograms - """ -# fig1 = P.figure() - -# ax_1 = fig1.add_subplot(211) -# ax_1.hist(self.all_pt, bins = 60, alpha=0.5, label="All Muons pt", log = True) -# ax_1.set_xlabel("Transverse Momentum") -# ax_1.set_ylabel("frequency") -# ax_1.legend(loc='upper right') -# P.show() - - #c1 = ROOT.TCanvas( 'c1', 'Muons', 1) - h_pt= TH1F( 'self.Muon_pt', 'Muons Transverse Momentun', 50, -2, 300 ) - for i in range (0, self.Muon_pt.size()): - - h_pt.Fill(self.Muon_pt[i]) - h_pt.Draw() - #c1.Update()