Upgrade

CMGTools/CaloUpgrade/BuildFile.xml

@@ -0,0 +1,12 @@
+<use   name="root"/>
+<use   name="rootcore"/>
+<use   name="rootinteractive"/>
+<use   name="rootrflx"/>
+<use   name="FWCore/Framework"/>
+<use   name="FWCore/ParameterSet"/>
+<use   name="FWCore/Sources"/>
+<use   name="PhysicsTools/Utilities"/>
+<use   name="boost"/>
+<export>
+   <lib name="1"/>
+</export>

CMGTools/CaloUpgrade/python/analyzers/HCALShowerAnalyzer.py

@@ -0,0 +1,75 @@
+import operator 
+import itertools
+import copy
+from math import fabs,sqrt
+from sets import Set
+
+from CMGTools.RootTools.fwlite.Analyzer import Analyzer
+from CMGTools.RootTools.fwlite.Event import Event
+from CMGTools.RootTools.fwlite.AutoHandle import AutoHandle
+from CMGTools.CaloUpgrade.tools.DataFormats import ShowerFromChargedPion 
+import ROOT
+
+
+
+class HCALShowerAnalyzer( Analyzer ):
+
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        self.doVis=True
+        self.visFile=ROOT.TFile("visInput.root","RECREATE")
+        self.doVis=True
+        self.counter=0
+        super(HCALShowerAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+
+
+
+
+    def declareHandles(self):
+        ''' Here declare handles of all objects we possibly need
+        '''
+        super(HCALShowerAnalyzer, self).declareHandles()
+
+        self.handles['hcalHits'] = AutoHandle( ('particleFlowRecHitHCAL','',''),'std::vector<reco::PFRecHit>')
+        self.handles['hoHits'] = AutoHandle( ('particleFlowRecHitHO',''),'std::vector<reco::PFRecHit>')
+        self.handles['ecalHits'] = AutoHandle( ('particleFlowRecHitECAL',''),'std::vector<reco::PFRecHit>')
+        self.handles['puInfo'] = AutoHandle( ('addPileupInfo',''),'std::vector<PileupSummaryInfo>')
+        self.handles['tracks'] = AutoHandle( ('pfTrack',''),'std::vector<reco::PFRecTrack>')
+
+    def beginLoop(self):
+        super(HCALShowerAnalyzer,self).beginLoop()
+
+
+    def process(self, iEvent, event):
+        self.event = iEvent.eventAuxiliary().id().event()
+        self.readCollections( iEvent )
+        event.pu = self.handles['puInfo'].product()
+        event.puInteractions=event.pu[0].getTrueNumInteractions()
+
+        event.tracks = self.handles['tracks'].product()
+        event.hcalHits = self.handles['hcalHits'].product()
+        event.ecalHits = self.handles['ecalHits'].product()
+        event.hoHits = self.handles['hoHits'].product()
+
+        event.showers=[]
+
+        #loop on tracks and match them with clusters
+        for track in event.tracks:
+            shower = ShowerFromChargedPion(track,0.5)
+            for hit in event.ecalHits:
+                shower.addConstituent(hit,True)
+            for hit in (event.hcalHits):
+                shower.addConstituent(hit)
+            for hit in (event.hoHits):
+                shower.addConstituent(hit)
+            event.showers.append(shower)    
+            print shower
+
+
+        if self.doVis:
+            for i,shower in enumerate(event.showers):
+                shower.makeVisTree(self.visFile,"t_"+str(self.counter)+"_"+str(i))
+
+        self.counter=self.counter+1
+        return True
+
+

CMGTools/CaloUpgrade/python/analyzers/HCALShowerTree.py

@@ -0,0 +1,35 @@
+import operator 
+import itertools
+import copy
+from math import fabs
+from ROOT import TLorentzVector
+
+from CMGTools.RootTools.analyzers.TreeAnalyzerNumpy import TreeAnalyzerNumpy
+
+class HCALShowerTree( TreeAnalyzerNumpy ):
+
+    def declareVariables(self):
+
+        self.branch('nPU')
+        self.branch('trackPt')
+        self.branch('trackP')
+        self.branch('trackEta')
+        self.branch('NECAL')
+        self.branch('NHCAL')
+
+        super(HCALShowerTree, self).declareVariables()
+
+
+
+    def process(self, iEvent, event):
+        for shower in event.showers:
+            self.reset()
+            self.set('nPU',event.puInteractions)
+            self.set('trackPt',shower.trackVector.Pt())
+            self.set('trackEta',shower.trackVector.Eta())
+            self.set('trackP',shower.trackMomentum)
+            self.set('NECAL',float(len(shower.ecalConstituents)))
+            self.set('NHCAL',float(len(shower.hcalConstituents)))
+            self.fill()
+
+

CMGTools/CaloUpgrade/python/tools/DataFormats.py

@@ -0,0 +1,160 @@
+from ROOT import TLorentzVector
+import ROOT
+from math import fabs,sqrt
+from CMGTools.RootTools.utils.DeltaR import deltaR,deltaPhi
+from array import array
+import numpy
+
+
+class ShowerFromChargedPion(object):
+    def __init__(self,track,dr = 0.5):
+        self.track=track
+        self.dr = dr
+        self.ecalEntrance =track.extrapolatedPoint(4).position()
+        self.hcalEntrance =track.extrapolatedPoint(6).position()
+        self.trackMomentum = self.track.trajectoryPoints()[0].momentum().energy()
+        self.trackVector = self.track.trajectoryPoints()[0].momentum()
+        self.ecalVector=None
+        self.hcalVector=None
+        self.ecalConstituents=[] 
+        self.hcalConstituents=[] 
+        self.hcalTiming=[]
+
+        #for visualization
+        #type = 0 track 1 ECAL , 2 HCAL
+
+        self.typeO=[-1,0]
+        self.eta=[self.ecalEntrance.Eta(),self.hcalEntrance.Eta()]
+        self.phi=[self.ecalEntrance.Phi(),self.hcalEntrance.Phi()]
+        self.ieta=[0,0]
+        self.iphi=[0,0]
+
+        self.rho=[self.ecalEntrance.Rho(),self.ecalEntrance.Rho()]
+        self.energy=[self.trackMomentum,self.trackMomentum]
+        self.depth=[-1,-1]
+
+
+    def vectorFromConstituent(self,cluster):
+        vec = ROOT.TVector3(cluster.position().x(),cluster.position().y(),cluster.position().z())
+        vec= vec.Unit()
+        vec*=cluster.energy()
+        return ROOT.TLorentzVector(vec.x(),vec.y(),vec.z(),cluster.energy())
+
+    def addConstituent(self,constituent,ecal = False):
+        vec = self.vectorFromConstituent(constituent)
+
+        if ecal and deltaR(vec.Eta(),vec.Phi(),self.trackVector.Eta(),self.trackVector.Phi())<self.dr:
+            self.ecalConstituents.append(constituent)
+            if self.ecalVector is None:
+                self.ecalVector = vec
+            else:    
+                self.ecalVector += vec
+            self.typeO.append(1)
+            self.depth.append(0)
+            self.eta.append(constituent.position().Eta()-self.ecalEntrance.Eta())
+            self.phi.append(deltaPhi(constituent.position().Phi(),self.ecalEntrance.Phi()))
+            self.rho.append(constituent.position().Rho()-self.ecalEntrance.Rho())
+            self.energy.append(constituent.energy())
+            self.ieta.append(0)
+            self.iphi.append(0)
+
+            return True
+
+        elif  deltaR(vec.Eta(),vec.Phi(),self.trackVector.Eta(),self.trackVector.Phi())<self.dr:
+            self.hcalConstituents.append(constituent)
+            self.hcalTiming.append(constituent.time())
+            if self.hcalVector is None:
+                self.hcalVector = vec
+            else:    
+                self.hcalVector += vec
+            self.typeO.append(2)
+
+            hcalid = ROOT.HcalDetId(constituent.detId())
+            self.depth.append(hcalid.depth())
+            self.ieta.append(hcalid.ieta())
+            self.iphi.append(hcalid.iphi())
+            self.eta.append(constituent.position().Eta()-self.ecalEntrance.Eta())
+            self.phi.append(deltaPhi(constituent.position().Phi(),self.ecalEntrance.Phi()))
+            self.rho.append(constituent.position().Rho()-self.ecalEntrance.Rho())
+            self.energy.append(constituent.energy())
+            return True
+        else:
+            return False
+
+
+
+    def calculateTrackDetId(self):
+        if len(self.hcalConstituents)>0:
+            constituents=sorted(self.hcalConstituents,key = lambda x: deltaR(x.position().Eta(),x.position().Phi(),self.ecalEntrance.Eta(),self.ecalEntrance.Phi()))
+            hcalid = ROOT.HcalDetId(constituents[0].detId())
+            self.ieta[1] = hcalid.ieta()
+            self.iphi[1] = hcalid.iphi()
+
+    def makeVisTree(self,tfile,treename):
+        tfile.cd()
+        self.calculateTrackDetId()
+
+        tree =ROOT.TTree(treename,treename)
+
+        typeO=numpy.zeros(1,float)
+        tree.Branch("type",typeO,'type/D')
+
+        eta=numpy.zeros(1,float)
+        tree.Branch("eta",eta,'eta/D')
+
+        phi=numpy.zeros(1,float)
+        tree.Branch("phi",phi,'phi/D')
+
+        ieta=numpy.zeros(1,float)
+        tree.Branch("ieta",ieta,'ieta/D')
+
+        iphi=numpy.zeros(1,float)
+        tree.Branch("iphi",iphi,'iphi/D')
+
+        rho=numpy.zeros(1,float)
+        tree.Branch("rho",rho,'rho/D')
+
+        energy=numpy.zeros(1,float)
+        tree.Branch("energy",energy,'energy/D')
+
+        depth=numpy.zeros(1,float)
+        tree.Branch("depth",depth,'depth/D')
+
+        for t,e,p,r,en,d,ie,ip in zip(self.typeO,self.eta,self.phi,self.rho,self.energy,self.depth,self.ieta,self.iphi):
+            typeO[0]=t
+            eta[0]=e
+            phi[0]=p
+            rho[0]=r
+            energy[0]=en
+            depth[0]=d
+            ieta[0]=ie
+            iphi[0]=ip
+            tree.Fill()
+        tree.Write()
+
+
+
+
+    def __str__(self):
+        str = """
+        New Shower from Pion
+        --------------------
+        """
+        str=str+ 'track pt,eta,phi,p,ecal eta, ecal phi, hcal eta, hcal phi  {pt},{eta},{phi},{p},{eeta},{ephi},{heta},{hphi} '.format(pt=self.trackVector.Pt(),eta=self.trackVector.Eta(),phi=self.trackVector.Phi(),p=self.trackMomentum,eeta=self.ecalEntrance.Eta(),ephi=self.ecalEntrance.Phi(),heta=self.hcalEntrance.Eta(),hphi=self.hcalEntrance.Phi())+'\n'
+        str=str+  """ECAL constituents
+                     -----------------
+                     """
+        for c in self.ecalConstituents:
+            str=str+ 'ID={id},Drho={rho},eta={eta},Dphi={phi},z={z},energy={energy} \n'.format(id= c.detId(),rho=c.position().rho(),eta=c.position().eta(),phi=c.position().phi(),z=c.position().z(),energy=c.energy())
+
+        str=str+  """HCAL constituents
+                     -----------------
+                     """
+        for c in self.hcalConstituents:
+            hcalid = ROOT.HcalDetId(c.detId())
+
+            str=str+ 'ID={id},Drho={rho},Deta={eta},Dphi={phi},z={z},energy={energy},depth={depth}\n'.format(id= c.detId(),rho=c.position().rho(),eta=c.position().eta(),phi=c.position().phi(),z=c.position().z(),energy=c.energy(),depth=hcalid.depth())
+
+
+        return str
+

CMGTools/CaloUpgrade/python/tools/Visualizer.py

@@ -0,0 +1,100 @@
+import ROOT
+from array import array
+class Visualizer(object):
+    def __init__(self,filename):
+        self.f = ROOT.TFile(filename)
+        self.hcal3D = self.createHCALHistogram()
+        self.legend,self.hcalLayers=self.createHCALLayers()
+        self.ecal = self.createECALHistogram()
+
+    def exit(self):
+        self.f.Close()
+
+
+    def createECALHistogram(self):
+        h = ROOT.TH2D("ecal","ecal",60,-0.522,0.522,60,-0.522,0.522)
+        h.GetXaxis().SetTitle(" #eta")
+        h.GetYaxis().SetTitle(" #phi")
+        h.SetLineColor(ROOT.kOrange)    
+
+        return h
+
+    def createHCALHistogram(self):
+        h = ROOT.TH3D("hcal","hcal",16,-8,8,16,-8,8,5,0,5)
+        h.GetXaxis().SetTitle("i #eta")
+        h.GetYaxis().SetTitle("i #phi")
+        h.GetZaxis().SetTitle("depth")
+
+        return h
+
+    def createHCALLayers(self):
+        layers=[]
+        legend =ROOT.TLegend(0.8,0.8,0.9,0.9)
+        legend.SetBorderSize(0)
+        legend.SetFillStyle(0)
+        for i in range(1,5+1):
+            layers.append(ROOT.TH2D("hcal_"+str(i),"hcal",16,-8,8,16,-8,8))
+            layers[-1].GetXaxis().SetTitle("i #eta")
+            layers[-1].GetYaxis().SetTitle("i #phi")
+            legend.AddEntry(layers[-1],"Layer "+str(i),"l")
+        layers[0].SetLineColor(ROOT.kRed)    
+        layers[1].SetLineColor(ROOT.kGreen)    
+        layers[2].SetLineColor(ROOT.kMagenta)    
+        layers[3].SetLineColor(ROOT.kYellow)    
+        layers[3].SetLineColor(ROOT.kOrange)    
+        layers[3].SetLineColor(ROOT.kAzure)    
+        return legend,layers
+
+
+    def process(self,event,shower,ecalThreshold=0.0,hcalThreshold=0.0):
+        self.hcal3D.Reset()
+        for layer in self.hcalLayers:
+            layer.Reset()
+
+        tree = self.f.Get("t_"+str(event)+"_"+str(shower))
+        iEta=0
+        iPhi=0
+        for event in tree:
+            if event.type>-0.5 and event.type<0.5:
+                iEta=event.ieta
+                iPhi=event.iphi
+                print 'Track energy',event.energy
+
+
+            if event.type>0.5 and event.type<1.5 and event.energy>ecalThreshold:
+                self.ecal.Fill(event.eta,event.phi,event.energy)
+                print 'ECAL',event.eta,event.phi,event.energy
+
+
+            if event.type>1.5 and event.type<2.5 and event.energy>hcalThreshold:
+                deltaPhi = event.iphi-iPhi
+                if deltaPhi>40:
+                    deltaPhi = deltaPhi-72
+                elif deltaPhi<-40:    
+                    deltaPhi = deltaPhi+72
+
+                self.hcal3D.Fill(event.ieta-iEta,deltaPhi,event.depth,event.energy)
+                self.hcalLayers[int(event.depth)-1].Fill(event.ieta-iEta,deltaPhi,event.energy)
+                print 'HCAL',event.ieta,event.iphi,event.depth,event.energy
+
+        self.canvases=[(ROOT.TCanvas("c3d",""))]
+        self.canvases[-1].cd()
+        self.hcal3D.Draw("box")
+
+        self.canvases.append(ROOT.TCanvas("c32d",""))
+        self.canvases[-1].cd()
+        sortedLayers=sorted(self.hcalLayers,key=lambda x : x.GetMaximum(),reverse=True)
+        for i,layer in enumerate(sortedLayers):
+            if i==0:
+                layer.Draw("box")
+            else:    
+                layer.Draw("box,same")
+
+        self.legend.Draw()        
+
+        self.canvases.append(ROOT.TCanvas("ecal",""))
+        self.canvases[-1].cd()
+        self.ecal.Draw("box")
+
+
+