diff --git a/LambaAnalyzer/src/LambdaAnalyzer.cc b/LambaAnalyzer/src/LambdaAnalyzer.cc index 5ee5dba..0f21b26 100755 --- a/LambaAnalyzer/src/LambdaAnalyzer.cc +++ b/LambaAnalyzer/src/LambdaAnalyzer.cc @@ -26,6 +26,14 @@ #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/TrackFitters/interface/TrajectoryFitterRecord.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" + #include "DataFormats/Math/interface/LorentzVector.h" #include "DataFormats/Math/interface/Vector.h" @@ -52,10 +60,17 @@ #include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "RecoTracker/TrackProducer/interface/KfTrackProducerBase.h" +#include "RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h" class LambdaAnalyzer : public edm::EDAnalyzer { + using Base = AlgoProductTraits; + //using TrackCollection = typename Base::TrackCollection; + using AlgoProductCollection = typename Base::AlgoProductCollection; + using TrackView = typename Base::TrackView; public: explicit LambdaAnalyzer(const edm::ParameterSet&); ~LambdaAnalyzer(); @@ -68,6 +83,8 @@ class LambdaAnalyzer : public edm::EDAnalyzer { void loop(const edm::Event& iEvent, const edm::EventSetup&, const reco::Vertex& RecVtx); void assignStableDaughters(const reco::Candidate* p, std::vector & pids); void initialize(); + int nSharedPixelLayerHits(reco::TransientTrack *track1, reco::TransientTrack *track2); +// void getFromEvt(edm::Event& theEvent,edm::Handle& theTCollection); // ----------member data --------------------------- bool doGen, doK3pi, doKpi; double m_pi, m_K, m_p; @@ -78,11 +95,11 @@ class LambdaAnalyzer : public edm::EDAnalyzer { std::vector t_tks; std::vector t_tks_noPXB1; TTree *tree1; + TrackProducerAlgorithm theAlgo; const std::vector *tPC; //std::vector trkParticle; - //ntuple variables int NKpiCand,NK3piCand,trigflag[160],NKpiMC,NK3piMC; @@ -96,7 +113,9 @@ class LambdaAnalyzer : public edm::EDAnalyzer { std::vector LambdaVtxPosx3,LambdaVtxPosy3,LambdaVtxPosz3,LambdaVtxerrx3,LambdaVtxerry3,LambdaVtxerrz3,LambdaVtx3cxy,LambdaVtx3cxz,LambdaVtx3cyz,LambdaetaK3pi,LambdaphiK3pi; std::vector DSetaK3pi,DSphiK3pi,LambdaMassK3proton,DSMassK3proton; - std::vector LambdaSharedHitLayer,LambdaSharedHitLadder,LambdaSharedHitModule; + std::vector LambdaSharedHitLayer,LambdaSharedHitLadder,LambdaSharedHitModule,LambdaSharedHitPixelHits_x,LambdaSharedHitPixelHits_y,LambdaSharedHitPixelHits_adc; + std::vector ProtonPixelHit_x,ProtonPixelHit_y,ProtonPixelHit_adc,PionPixelHit_x,PionPixelHit_y,PionPixelHit_adc; + int PionPixelHitLayer,ProtonPixelHitLayer; //Gen quantities std::vector GenLambdaVtxPosx, GenLambdaVtxPosy, GenLambdaVtxPosz, GenLambdaSourceVtxPosx, GenLambdaSourceVtxPosy, GenLambdaSourceVtxPosz,GenLambdaPt, GenLambdaP, GenLambdaPhi, GenLambdaEta, GenLambdaMass, GenLambdaMt, GenLambdaE, GenLambdaEt, GenLambdaPx, GenLambdaPy, GenLambdaPz, GenFlightLength, GenDeltaR; @@ -140,7 +159,9 @@ class LambdaAnalyzer : public edm::EDAnalyzer { edm::EDGetTokenT TrackCollT_noPXB1_; edm::EDGetTokenT VtxCollT_; edm::EDGetTokenT GenCollT_; - + edm::EDGetTokenT BSCollT_; + edm::EDGetTokenT> theTCCollection; + edm::Handle genParticles; typedef edm::AssociationMap > TrackVertexAssMap; @@ -156,7 +177,8 @@ class LambdaAnalyzer : public edm::EDAnalyzer { LambdaAnalyzer::LambdaAnalyzer(const edm::ParameterSet& iConfig): doGen(iConfig.getParameter("doGen")), doK3pi(iConfig.getParameter("doK3pi")), - doKpi(iConfig.getParameter("doKpi")) + doKpi(iConfig.getParameter("doKpi")), + theAlgo(iConfig) { //now do what ever initialization is needed edm::Service fs; @@ -168,6 +190,8 @@ LambdaAnalyzer::LambdaAnalyzer(const edm::ParameterSet& iConfig): VtxCollT_ = consumes(iConfig.getUntrackedParameter("vertices")); GenCollT_ = consumes(iConfig.getUntrackedParameter("genParticles")); T2VCollT_ = consumes(iConfig.getUntrackedParameter("T2V")); + BSCollT_ = consumes(iConfig.getUntrackedParameter("BeamSpot")); + theTCCollection = consumes>(iConfig.getUntrackedParameter("trackCandidates")); //MixCollT_ = consumes(iConfig.getUntrackedParameter("mix")); TrackingParticleCollT_ = consumes(iConfig.getParameter("trackingParticles") ); @@ -252,9 +276,9 @@ void LambdaAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iS BSerrz=vertexBeamSpot.z0Error();*/ - Handle generalTracks; + Handle generalTracks; iEvent.getByToken(TrackCollT_, generalTracks); - Handle generalTracks_noPXB1; + Handle generalTracks_noPXB1; iEvent.getByToken(TrackCollT_noPXB1_, generalTracks_noPXB1); //iEvent.getByLabel("generalTracks",generalTracks); @@ -567,19 +591,85 @@ void LambdaAnalyzer::loop(const edm::Event& iEvent, const edm::EventSetup& iSetu // } // } // std::cout<<"nSharedPBHits = "<recHitsBegin(); trackingRecHit_iterator hb_proton = proton->recHitsBegin(); TrackingRecHit const * h1[4] = { (*hb_pi1), (*(hb_pi1+1)), (*(hb_pi1+2)), (*(hb_pi1+3)) }; TrackingRecHit const * h2[4] = { (*hb_proton), (*(hb_proton+1)), (*(hb_proton+2)), (*(hb_proton+3)) }; + const SiPixelRecHit* pixelhit_p = dynamic_cast(h2[0]); + if(pixelhit_p!=nullptr && h2[0]->isValid() && LambdaMass.size()==0) { + // only fill for the first row in the event since we can't have a vector of vectors + std::vector pixels_p(pixelhit_p->cluster()->pixels()); + // std::cout<<"pixels.size() = "<clusterProbability(0) = "<clusterProbability(0)<clusterProbability(1) = "<clusterProbability(1)<hasFilledProb() = "<geographicalId(); + ProtonPixelHitLayer = pxb_id_p.layer(); + + } + else { + ProtonPixelHit_x.push_back(-99); + ProtonPixelHit_y.push_back(-99); + ProtonPixelHit_adc.push_back(-99); + ProtonPixelHitLayer = -1; + } + const SiPixelRecHit* pixelhit_pi = dynamic_cast(h1[0]); + if(pixelhit_pi!=nullptr && h1[0]->isValid() && LambdaMass.size()==0) { + // only fill for the first row in the event since we can't have a vector of vectors + std::vector pixels_pi(pixelhit_pi->cluster()->pixels()); + // std::cout<<"pixels.size() = "<clusterProbability(0) = "<clusterProbability(0)<clusterProbability(1) = "<clusterProbability(1)<hasFilledProb() = "<geographicalId(); + PionPixelHitLayer = pxb_id_pi.layer(); + + } + else { + PionPixelHit_x.push_back(-99); + PionPixelHit_y.push_back(-99); + PionPixelHit_adc.push_back(-99); + PionPixelHitLayer = -1; + } //std::cout<<"new lambda pairing"<isValid()) continue; for (int k2=0; k2<4; k2++) { if (!h2[k2]->isValid()) continue; bool shared = h1[k1]->sharesInput(h2[k2],TrackingRecHit::some); if (shared) { + const SiPixelRecHit* pixelhit = dynamic_cast(h1[k1]); + //pixelhit->cluster()->pixels(); + if(pixelhit!=nullptr) { + std::vector pixels(pixelhit->cluster()->pixels()); + if (LambdaSharedHitLayer.size() == 0) { + // can't have a vector of vectors so I don't have a better way to do this now then only deal with the first + // Lambda if there are more than one for this vertex + //std::cout<<"pixels.size() = "<geographicalId()< 6.) { + // // shared hit at layer 0 but clearly coming from layer 2 + // std::cout<<"found outlier shared hit at layer 0!"<geographicalId())<det().subDetector().< 0.25){ - if( (t_trk.track().numberOfValidHits() >= 6) && (t_trk.track().pt() > 0.35) && - fabs(t_trk.track().dz(RecVtx.position()))<2.0 && - fabs(t_trk.track().dxy(RecVtx.position()) / t_trk.track().d0Error()) > 2.0 ) { - double deta_pi1 = fabs(pi1->track().eta() - t_trk.track().eta()); - double dphi_pi1 = fabs(pi1->track().phi() - t_trk.track().phi()); - double dr_pi1 = TMath::Sqrt(TMath::Power(deta_pi1,2) + TMath::Power(dphi_pi1,2)); - if (dr_pi1 < minDR_pi1) { - minDR_pi1 = dr_pi1; - pi1_noPXB1 = t_trk; - } - double deta_proton = fabs(proton->track().eta() - t_trk.track().eta()); - double dphi_proton = fabs(proton->track().phi() - t_trk.track().phi()); - double dr_proton = TMath::Sqrt(TMath::Power(deta_proton,2) + TMath::Power(dphi_proton,2)); - if (dr_proton < minDR_proton && minDR_pi1 != dr_pi1) { - minDR_proton = dr_proton; - proton_noPXB1 = t_trk; - } + //if (foundSharedHitLayer0) { + double minDR_pi1 = 999.; + double minDR_proton = 999.; + TransientTrack pi1_noPXB1; + TransientTrack proton_noPXB1; + //std::cout<<"starting new loop"< 0.25){ + if( (t_trk.track().numberOfValidHits() >= 6) && (t_trk.track().pt() > 0.35) && + fabs(t_trk.track().dz(RecVtx.position()))<2.0 && + fabs(t_trk.track().dxy(RecVtx.position()) / t_trk.track().d0Error()) > 2.0 ) { + double deta_pi1 = fabs(pi1->track().eta() - t_trk.track().eta()); + double dphi_pi1 = fabs(pi1->track().phi() - t_trk.track().phi()); + double dr_pi1 = TMath::Sqrt(TMath::Power(deta_pi1,2) + TMath::Power(dphi_pi1,2)); + if (dr_pi1 < minDR_pi1) { + minDR_pi1 = dr_pi1; + pi1_noPXB1 = t_trk; + } + double deta_proton = fabs(proton->track().eta() - t_trk.track().eta()); + double dphi_proton = fabs(proton->track().phi() - t_trk.track().phi()); + double dr_proton = TMath::Sqrt(TMath::Power(deta_proton,2) + TMath::Power(dphi_proton,2)); + if (dr_proton < minDR_proton && minDR_pi1 != dr_pi1) { + minDR_proton = dr_proton; + proton_noPXB1 = t_trk; + } + int nSharedHits = nSharedPixelLayerHits(pi1,&t_trk); + //std::cout<<"going to run a sanity check on the no-PXB1 track"<hitPattern(); + for (int i = 0; i < p2.numberOfAllHits(HitPattern::TRACK_HITS); i++) { + uint32_t hit = p2.getHitPattern(HitPattern::TRACK_HITS, i); + // if the hit is valid and in pixel barrel, print out the layer + if (p2.validHitFilter(hit) && p2.pixelBarrelHitFilter(hit)){ + //std::cout << "valid hit found in pixel barrel layer " + // << p2.getLayer(hit) + // << std::endl; + } + } + //std::cout<<"nSharedHits = "<track().eta()<<" : "<track().phi()< tks_noPXB1; - tks_noPXB1.push_back(pi1_noPXB1); - tks_noPXB1.push_back(proton_noPXB1); - //KalmanVertexFitter kalman(true); - TransientVertex v_noPXB1 = kalman.vertex(tks_noPXB1); - if(v_noPXB1.isValid() && v_noPXB1.hasRefittedTracks()) { - //double vtxProb =TMath::Prob( (Double_t) v.totalChiSquared(), (Int_t) v.degreesOfFreedom()); - //if (vtxProb < 0.05) continue; - //TransientTrack pi1_f_noPXB1 = v_noPXB1.refittedTrack(*pi1_noPXB1); - //TransientTrack proton_f_noPXB1 = v_noPXB1.refittedTrack(*proton_noPXB1); - math::XYZVector Lambda_position_noPXB1 = math::XYZVector(v_noPXB1.position().x(), v_noPXB1.position().y(), v_noPXB1.position().z() ); - - math::XYZVector displacement_noPXB1 = Lambda_position_noPXB1 - PV_position; - flightlength_noPXB1 = sqrt( displacement_noPXB1.perp2()); - //std::cout<<"flight length with PXB1: "< tks_noPXB1; + tks_noPXB1.push_back(pi1_noPXB1); + tks_noPXB1.push_back(proton_noPXB1); + //KalmanVertexFitter kalman(true); + TransientVertex v_noPXB1 = kalman.vertex(tks_noPXB1); + if(v_noPXB1.isValid() && v_noPXB1.hasRefittedTracks()) { + //double vtxProb =TMath::Prob( (Double_t) v.totalChiSquared(), (Int_t) v.degreesOfFreedom()); + //if (vtxProb < 0.05) continue; + //TransientTrack pi1_f_noPXB1 = v_noPXB1.refittedTrack(*pi1_noPXB1); + //TransientTrack proton_f_noPXB1 = v_noPXB1.refittedTrack(*proton_noPXB1); + math::XYZVector Lambda_position_noPXB1 = math::XYZVector(v_noPXB1.position().x(), v_noPXB1.position().y(), v_noPXB1.position().z() ); + + math::XYZVector displacement_noPXB1 = Lambda_position_noPXB1 - PV_position; + flightlength_noPXB1 = sqrt( displacement_noPXB1.perp2()); + //std::cout<<"flight length with PXB1: "< theG; + //edm::ESHandle theMF; + //edm::ESHandle theFitter; + //edm::ESHandle thePropagator; + //edm::ESHandle theMeasTk; + //edm::ESHandle theBuilder; + //iSetup.get().get(theG); + ////getFromES(iSetup,theG,theMF,theFitter,thePropagator,theMeasTk,theBuilder); + //iSetup.get().get(theMF); + //iSetup.get().get("KFFittingSmootherWithOutliersRejectionAndRK",theFitter); + //iSetup.get().get("RungeKuttaTrackerPropagator",thePropagator); + //iSetup.get().get("",theMeasTk); + //iSetup.get().get("WithAngleAndTemplate",theBuilder); + + //edm::Handle recoBeamSpotHandle; + //iEvent.getByToken(BSCollT_,recoBeamSpotHandle); + //reco::BeamSpot bs = *recoBeamSpotHandle; + + //AlgoProductCollection algoResults; + //edm::Handle> theTCollection; + //AlgoProductTraits::AlgoProductCollection algoResults; + //AlgoProductTraits::AlgoProduct algoResults; + //edm::View theTCollection; + //reco::Track t_test; + // Track(double chi2, double ndof, const Point & referencePoint, + // const Vector & momentum, int charge, const CovarianceMatrix &, + // TrackAlgorithm = undefAlgorithm, TrackQuality quality = undefQuality, +// float t0 = 0, float beta = 0, +// float covt0t0 = -1., float covbetabeta = -1.); + //Handle test_generalTracks; + //iEvent.getByToken(TrackCollT_, test_generalTracks); + + //edm::Handle> theTCollection; + //edm::Handle theTCollection; + //edm::View theTCollection; + //edm::Handle theTCollection; + //iEvent.getByToken(theTCCollection,theTCollection); + //ckfTrackCandidates + //iEvent.getByLabel("generalTracks",theTCollection ); + + ///edm::Handle> theTCollection; + //getFromEvt(iEvent,theTCollection); + + //edm::Handle> theTCollection_skimmed; + //std::cout<<"theTCollection->size() = "<size()<at(0) = "<at(0)<at(0); + //std::cout<> to refit, i.e. just the + // p,pi tracks with layer 1 hits removed? + + //theAlgo.runWithTrack(theG.product(), theMF.product(), *theTCollection, + // theFitter.product(), thePropagator.product(), + // theBuilder.product(), bs, algoResults); + + //std::cout<<"algoResults.size() = "<Branch("LambdaSharedHitLayer",&LambdaSharedHitLayer); tree1->Branch("LambdaSharedHitLadder",&LambdaSharedHitLadder); tree1->Branch("LambdaSharedHitModule",&LambdaSharedHitModule); - +tree1->Branch("LambdaSharedHitPixelHits_x",&LambdaSharedHitPixelHits_x); +tree1->Branch("LambdaSharedHitPixelHits_y",&LambdaSharedHitPixelHits_y); +tree1->Branch("LambdaSharedHitPixelHits_adc",&LambdaSharedHitPixelHits_adc); +tree1->Branch("PionPixelHit_x",&PionPixelHit_x); +tree1->Branch("PionPixelHit_y",&PionPixelHit_y); +tree1->Branch("PionPixelHit_adc",&PionPixelHit_adc); +tree1->Branch("ProtonPixelHit_x",&ProtonPixelHit_x); +tree1->Branch("ProtonPixelHit_y",&ProtonPixelHit_y); +tree1->Branch("ProtonPixelHit_adc",&ProtonPixelHit_adc); +tree1->Branch("PionPixelHitLayer",&PionPixelHitLayer); +tree1->Branch("ProtonPixelHitLayer",&ProtonPixelHitLayer); //tracks tree1->Branch("ntracks",&ntracks,"ntracks/I"); @@ -1076,5 +1275,46 @@ void LambdaAnalyzer::endJob() { } + +int LambdaAnalyzer::nSharedPixelLayerHits(reco::TransientTrack *track1, reco::TransientTrack *track2) { + //std::cout<<"called nSharedPixelLayerHits"<recHitsBegin(); + trackingRecHit_iterator hb_t2 = track2->recHitsBegin(); + TrackingRecHit const * h1[7] = { (*hb_t1), (*(hb_t1+1)), (*(hb_t1+2)), (*(hb_t1+3)), (*(hb_t1+4)), (*(hb_t1+5)), (*(hb_t1+6)), }; + TrackingRecHit const * h2[7] = { (*hb_t2), (*(hb_t2+1)), (*(hb_t2+2)), (*(hb_t2+3)), (*(hb_t2+4)), (*(hb_t2+5)), (*(hb_t2+6)), }; + for (int k1=0; k1<7; k1++) { + if (!h1[k1]->isValid()) continue; + //std::cout<<"track2 # recHits = "<recHitsSize()<isValid()) continue; + //if (k1==0) { + //PXBDetId tmp = h2[k2]->geographicalId(); + //if (tmp.subdetId()== 1) { + //std::cout<<"k2 = "<sharesInput(h2[k2],TrackingRecHit::some); + if (shared) { + nShared += 1; + } + } + } + //if (nShared>=3) { + // std::cout<<"nShared >= 3"<isValid()) { + // PXBDetId pxb_id = h2[0]->geographicalId(); + // std::cout<<"layer = "<& theTCollection) { +// theEvent.getByLabel("generalTracks",theTCollection ); +//} + //define this as a plug-in DEFINE_FWK_MODULE(LambdaAnalyzer); diff --git a/LambaAnalyzer/test/addBranchesForNNTraining.py b/LambaAnalyzer/test/addBranchesForNNTraining.py new file mode 100644 index 0000000..72d4a33 --- /dev/null +++ b/LambaAnalyzer/test/addBranchesForNNTraining.py @@ -0,0 +1,114 @@ +import ROOT +import sys +from math import pow,sqrt +import numpy as np + +ifile = ROOT.TFile(sys.argv[1]) +ofile = ROOT.TFile("ofile2.root","RECREATE") +tree = ifile.Get("tree1") +otree = tree.CloneTree(0) + +nentries = tree.GetEntries() +print "nentries = ",nentries +gridSize = 16 +#nentries = 1200000 + +isSharedHit = np.zeros(1,dtype=int) +otree.Branch("isSharedHit",isSharedHit,"isSharedHit/I") +pixel_references = [0.]*gridSize*gridSize +for i in xrange(gridSize*gridSize): + pixel_references[i] = np.zeros(1,dtype=float) + otree.Branch("pixel_%i" % i,pixel_references[i],"pixel_%i/D" %i ) +trackPt = np.zeros(1,dtype=float) +trackEta = np.zeros(1,dtype=float) +trackPhi = np.zeros(1,dtype=float) +otree.Branch("trackPt",trackPt,"trackPt/D") +otree.Branch("trackEta",trackEta,"trackEta/D") +otree.Branch("trackPhi",trackPt,"trackPhi/D") + +def getPixelHist(pixels,gridSize): + xmin = -1 + xmax = -1 + ymin = -1 + ymax = -1 + xavg = 0. + yavg = 0. + tot_adc = 0. + for x,y,adc in pixels: + #print x,y,adc + if x < xmin or xmin == -1: + xmin = x + if y < ymin or ymin == -1: + ymin = y + xavg += x*adc + yavg += y*adc + tot_adc += adc + xavg = xavg / (tot_adc) + yavg = yavg / (tot_adc) + xavg_int = int(round(xavg)) + yavg_int = int(round(yavg)) + hist = ROOT.TH2F("hist_%i" % iEntry,"hist_%i" % iEntry,gridSize,0,gridSize,gridSize,0,gridSize) + for x,y,adc in pixels: + #print (x-xavg_int),(y-ymin),adc + hist.Fill(x-xavg_int+gridSize/2.,y-yavg_int+gridSize/2.,adc) + if hist.Integral() > 0: + hist.Scale(1./hist.Integral()) + else: + hist = ROOT.TH2F("hist_shared","hist_shared",gridSize,0,gridSize,gridSize,0,gridSize) + return hist + +#for iEntry in xrange(nentries): +for iEntry in xrange(1200000): +#for iEntry in xrange(1200000,nentries): + if (iEntry % 1000 == 0): + print "processing entry: ",iEntry + tree.GetEntry(iEntry) + pixels_shared = [] + pixels_pion = [] + pixels_proton = [] + if len(tree.PionPixelHit_x)>0 and tree.PionPixelHitLayer==0 and iEntry%100==0: + for i in xrange(len(tree.PionPixelHit_x)): + pixels_pion.append((tree.PionPixelHit_x[i],tree.PionPixelHit_y[i],tree.PionPixelHit_adc[i])) + hist = getPixelHist(pixels_pion,gridSize) + isSharedHit[0] = 0 + for i in xrange(gridSize): + for j in xrange(gridSize): + pixel_references[i+gridSize*j][0] = hist.GetBinContent(i+1,j+1) + trackPt[0] = tree.TrkPi1pt[0] + trackEta[0] = tree.TrkPi1eta[0] + trackPhi[0] = tree.TrkPi1phi[0] + otree.Fill() + if len(tree.ProtonPixelHit_x)>0 and tree.ProtonPixelHitLayer==0 and iEntry%100==0: + for i in xrange(len(tree.ProtonPixelHit_x)): + pixels_proton.append((tree.ProtonPixelHit_x[i],tree.ProtonPixelHit_y[i],tree.ProtonPixelHit_adc[i])) + hist = getPixelHist(pixels_proton,gridSize) + isSharedHit[0] = 0 + for i in xrange(gridSize): + for j in xrange(gridSize): + pixel_references[i+gridSize*j][0] = hist.GetBinContent(i+1,j+1) + trackPt[0] = tree.TrkProtonpt[0] + trackEta[0] = tree.TrkProtoneta[0] + trackPhi[0] = tree.TrkProtonphi[0] + otree.Fill() + if tree.LambdaMass[0] > 0 and len(tree.LambdaSharedHitPixelHits_x) > 0 and tree.LambdaSharedHitLayer[0]==0 and tree.flightLength[0]<4.: + for i in xrange(len(tree.LambdaSharedHitPixelHits_x)): + pixels_shared.append((tree.LambdaSharedHitPixelHits_x[i],tree.LambdaSharedHitPixelHits_y[i],tree.LambdaSharedHitPixelHits_adc[i])) + isSharedHit[0] = 1 + hist = getPixelHist(pixels_shared,gridSize) + for i in xrange(gridSize): + for j in xrange(gridSize): + pixel_references[i+gridSize*j][0] = hist.GetBinContent(i+1,j+1) + if tree.TrkPi1pt[0] > tree.TrkProtonpt[0]: + trackPt[0] = tree.TrkPi1pt[0] + trackEta[0] = tree.TrkPi1eta[0] + trackPhi[0] = tree.TrkPi1phi[0] + else: + trackPt[0] = tree.TrkProtonpt[0] + trackEta[0] = tree.TrkProtoneta[0] + trackPhi[0] = tree.TrkProtonphi[0] + + otree.Fill() + #dr = sqrt( pow((abs(tree.TrkProtonphi[0]-tree.TrkPi1phi[0])-(abs(tree.TrkProtonphi[0]-tree.TrkPi1phi[0])>3.14)*2*3.14),2) + pow((tree.TrkProtoneta[0]-tree.TrkPi1eta[0]),2) ) + +ofile.cd() +otree.Write() diff --git a/LambaAnalyzer/test/converRootToPandas.py b/LambaAnalyzer/test/converRootToPandas.py new file mode 100644 index 0000000..04fe515 --- /dev/null +++ b/LambaAnalyzer/test/converRootToPandas.py @@ -0,0 +1,11 @@ +from root_pandas import read_root +import sys + +cols = ['isSharedHit','trackPt','trackEta','trackPhi'] +for i in xrange(16*16): + cols.append('pixel_%i' % i) + +df = read_root(sys.argv[1], columns=cols) +print df +df.to_hdf("pixelTrain.h5",key='df',mode='w') +print "done" diff --git a/LambaAnalyzer/test/makePixelClusterPlots.py b/LambaAnalyzer/test/makePixelClusterPlots.py new file mode 100644 index 0000000..a8e4d32 --- /dev/null +++ b/LambaAnalyzer/test/makePixelClusterPlots.py @@ -0,0 +1,71 @@ +import ROOT +import sys +from math import pow,sqrt + +ifile = ROOT.TFile(sys.argv[1]) +ofile = ROOT.TFile("ofile.root","RECREATE") +tree = ifile.Get("tree1") + +nentries = tree.GetEntries() +print "nentries = ",nentries +gridSize = 16 +hist_shared = ROOT.TH2F("hist_shared","hist_shared",gridSize,0,gridSize,gridSize,0,gridSize) +hist_pion = ROOT.TH2F("hist_pion","hist_pion",gridSize,0,gridSize,gridSize,0,gridSize) +hist_proton = ROOT.TH2F("hist_proton","hist_proton",gridSize,0,gridSize,gridSize,0,gridSize) +#nentries = 1000 + +def getPixelHist(pixels,gridSize): + xmin = -1 + xmax = -1 + ymin = -1 + ymax = -1 + xavg = 0. + yavg = 0. + tot_adc = 0. + for x,y,adc in pixels: + #print x,y,adc + if x < xmin or xmin == -1: + xmin = x + if y < ymin or ymin == -1: + ymin = y + xavg += x*adc + yavg += y*adc + tot_adc += adc + xavg = xavg / (tot_adc) + yavg = yavg / (tot_adc) + xavg_int = int(round(xavg)) + yavg_int = int(round(yavg)) + hist = ROOT.TH2F("hist_%i" % iEntry,"hist_%i" % iEntry,gridSize,0,gridSize,gridSize,0,gridSize) + for x,y,adc in pixels: + #print (x-xavg_int),(y-ymin),adc + hist.Fill(x-xavg_int+gridSize/2.,y-yavg_int+gridSize/2.,adc) + if hist.Integral() > 0: + hist.Scale(1./hist.Integral()) + else: + hist = ROOT.TH2F("hist_shared","hist_shared",gridSize,0,gridSize,gridSize,0,gridSize) + return hist + +for iEntry in xrange(nentries): + if (iEntry % 1000 == 0): + print "processing entry: ",iEntry + tree.GetEntry(iEntry) + pixels_shared = [] + pixels_pion = [] + pixels_proton = [] + if len(tree.PionPixelHit_x)>0: + for i in xrange(len(tree.PionPixelHit_x)): + pixels_pion.append((tree.PionPixelHit_x[i],tree.PionPixelHit_y[i],tree.PionPixelHit_adc[i])) + hist_pion.Add(getPixelHist(pixels_pion,gridSize)) + if len(tree.ProtonPixelHit_x)>0: + for i in xrange(len(tree.ProtonPixelHit_x)): + pixels_proton.append((tree.ProtonPixelHit_x[i],tree.ProtonPixelHit_y[i],tree.ProtonPixelHit_adc[i])) + hist_proton.Add(getPixelHist(pixels_proton,gridSize)) + if tree.LambdaMass[0] > 0 and len(tree.LambdaSharedHitPixelHits_x) > 0 and tree.LambdaSharedHitLayer[0]==0 and tree.flightLength[0]<4.: + for i in xrange(len(tree.LambdaSharedHitPixelHits_x)): + pixels_shared.append((tree.LambdaSharedHitPixelHits_x[i],tree.LambdaSharedHitPixelHits_y[i],tree.LambdaSharedHitPixelHits_adc[i])) + hist_shared.Add(getPixelHist(pixels_shared,gridSize)) + #dr = sqrt( pow((abs(tree.TrkProtonphi[0]-tree.TrkPi1phi[0])-(abs(tree.TrkProtonphi[0]-tree.TrkPi1phi[0])>3.14)*2*3.14),2) + pow((tree.TrkProtoneta[0]-tree.TrkPi1eta[0]),2) ) +ofile.cd() +hist_shared.Write() +hist_pion.Write() +hist_proton.Write() diff --git a/LambaAnalyzer/test/trkeffanalyzer_Data_GeneralTracks_cfg.py b/LambaAnalyzer/test/trkeffanalyzer_Data_GeneralTracks_cfg.py index bea3839..8f10194 100755 --- a/LambaAnalyzer/test/trkeffanalyzer_Data_GeneralTracks_cfg.py +++ b/LambaAnalyzer/test/trkeffanalyzer_Data_GeneralTracks_cfg.py @@ -26,7 +26,8 @@ #process.GlobalTag.globaltag = 'GR_P_V56::All' process.GlobalTag.globaltag = '101X_dataRun2_Prompt_v11' -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10) ) +#process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000) ) #process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10000) ) #filelist = FileUtils.loadListFromFile("inputlist.list") @@ -159,6 +160,11 @@ vertices = cms.untracked.InputTag("offlinePrimaryVertices"), genParticles = cms.untracked.InputTag("genParticles"), T2V = cms.untracked.InputTag("Tracks2Vertex"), + + AlgorithmName = cms.string('undefAlgorithm'), + BeamSpot = cms.untracked.InputTag('offlineBeamSpot'), + trackCandidates = cms.untracked.InputTag('generalTracks'), + #trackCandidates = cms.untracked.InputTag('ckfTrackCandidates'), ) process.TFileService = cms.Service("TFileService", @@ -193,6 +199,7 @@ process.TrackerTrackHitFilter.src = 'TrackRefitter1' #process.TrackerTrackHitFilter.src = 'generalTracks' process.TrackerTrackHitFilter.commands = cms.vstring("drop PXB","keep PXB 2","keep PXB 3","keep PXB 4","keep PXE","keep TIB","keep TID","keep TOB","keep TEC") +#process.TrackerTrackHitFilter.commands = cms.vstring("drop PXB","keep PXB 2","keep PXB 3","keep PXB 4","keep PXE","keep TIB","keep TID","keep TOB","keep TEC") #process.TrackerTrackHitFilter.useTrajectories= True # this is needed only if you require some selections; but it will work even if you don't ask for them #process.TrackerTrackHitFilter.minimumHits = 6 diff --git a/LambaAnalyzer/test/trkeffanalyzer_MC_GeneralTracks_cfg.py b/LambaAnalyzer/test/trkeffanalyzer_MC_GeneralTracks_cfg.py index 178cf30..ddc151b 100755 --- a/LambaAnalyzer/test/trkeffanalyzer_MC_GeneralTracks_cfg.py +++ b/LambaAnalyzer/test/trkeffanalyzer_MC_GeneralTracks_cfg.py @@ -167,6 +167,9 @@ T2V = cms.untracked.InputTag("Tracks2Vertex"), trackingParticles = cms.InputTag("mix","MergedTrackTruth"), trackingVertices = cms.InputTag("mix","MergedTrackTruth"), + AlgorithmName = cms.string('undefAlgorithm'), + BeamSpot = cms.untracked.InputTag('offlineBeamSpot'), + trackCandidates = cms.untracked.InputTag('generalTracks'), ) process.TFileService = cms.Service("TFileService",