#include <ElectronSeedAnalyzer.h>

Inheritance diagram for ElectronSeedAnalyzer:

Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
	ElectronSeedAnalyzer (const edm::ParameterSet &conf)
virtual void	endJob ()
virtual	~ElectronSeedAnalyzer ()
Private Attributes
TH1F *	histecaldriveneMCmatched_
TH1F *	histecaldrivenetaMCmatched_
TH1F *	histecaldrivenptMCmatched_
TH1F *	histeclu_
TH1F *	histeffeta_
TH1F *	histeffpt_
TH1F *	histeMC_
TH1F *	histeMCmatched_
TH1F *	histeoverp_
TH1F *	histeta_
TH1F *	histetaclu_
TH1F *	histetaMC_
TH1F *	histetaMCmatched_
TH1F *	histetclu_
TFile *	histfile_
TH1I *	histnbclus_
TH1I *	histnbseeds_
TH1I *	histnrseeds_
TH1F *	histp_
TH1F *	histpt_
TH1F *	histptMC_
TH1F *	histptMCmatched_
TH1F *	histq_
TH1F *	histtrackerdriveneMCmatched_
TH1F *	histtrackerdrivenetaMCmatched_
TH1F *	histtrackerdrivenptMCmatched_
edm::InputTag	inputCollection_
float	mcEnergy [10]
float	mcEta [10]
float	mcPhi [10]
float	mcPt [10]
float	mcQ [10]
float	seedDphi1 [10]
float	seedDphi2 [10]
float	seedDrz1 [10]
float	seedDrz2 [10]
float	seedEta [10]
int	seedLayer1 [10]
int	seedLayer2 [10]
float	seedMomentum [10]
float	seedPhi [10]
float	seedPhi1 [10]
float	seedPhi2 [10]
float	seedPt [10]
float	seedQ [10]
float	seedRz1 [10]
float	seedRz2 [10]
int	seedSide1 [10]
int	seedSide2 [10]
int	seedSubdet1 [10]
int	seedSubdet2 [10]
float	superclusterEnergy [10]
float	superclusterEt [10]
float	superclusterEta [10]
float	superclusterPhi [10]
TrajectoryStateTransform	transformer_
TTree *	tree_

Detailed Description

Definition at line 35 of file ElectronSeedAnalyzer.h.

Constructor & Destructor Documentation

ElectronSeedAnalyzer::ElectronSeedAnalyzer ( const edm::ParameterSet & conf ) [explicit]

Definition at line 62 of file ElectronSeedAnalyzer.cc.

References edm::ParameterSet::getParameter().

{
  inputCollection_=conf.getParameter<edm::InputTag>("inputCollection") ;
  histfile_ = new TFile("electronpixelseeds.root","RECREATE");
}

ElectronSeedAnalyzer::~ElectronSeedAnalyzer ( ) [virtual]

Definition at line 176 of file ElectronSeedAnalyzer.cc.

{
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  //tree_->Print();
  histfile_->Write();
  histeMC_->Write();
  histfile_->Close();
}

Member Function Documentation

void ElectronSeedAnalyzer::analyze	(	const edm::Event &	e,
		const edm::EventSetup &	iSetup
	)		`[virtual]`

Implements edm::EDAnalyzer.

Definition at line 186 of file ElectronSeedAnalyzer.cc.

References abs, alongMomentum, edm::RefToBase< T >::castTo(), DeDxDiscriminatorTools::charge(), deltaR(), DetId::det(), PXFDetId::disk(), relval_parameters_module::energy, eta(), PV3DBase< T, PVType, FrameType >::eta(), funct::exp(), configurableAnalysis::GenParticle, edm::EventSetup::get(), edm::Event::getByLabel(), edm::Event::getByType(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), edm::EventBase::id(), edm::RefToBase< T >::isNull(), TrajectoryStateOnSurface::isValid(), TIBDetId::layer(), TOBDetId::layer(), PXBDetId::layer(), LogDebug, PV3DBase< T, PVType, FrameType >::mag(), scaleCards::mass, FreeTrajectoryState::momentum(), oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, parents, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, pi, FreeTrajectoryState::position(), funct::pow(), edm::Handle< T >::product(), PropagatorWithMaterial::propagate(), alignCSCRings::r, DetId::rawId(), funct::sin(), mathSSE::sqrt(), DetId::subdetId(), GeomDet::surface(), lumiQTWidget::t, Surface::toGlobal(), GeomDet::toGlobal(), trajectoryStateTransform::transientState(), TIDDetId::wheel(), TECDetId::wheel(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

{

  edm::ESHandle<TrackerGeometry> pDD ;
  edm::ESHandle<MagneticField> theMagField ;
  iSetup.get<TrackerDigiGeometryRecord> ().get(pDD);
  iSetup.get<IdealMagneticFieldRecord>().get(theMagField);

  // rereads the seeds for test purposes
  typedef edm::OwnVector<TrackingRecHit> recHitContainer;
  typedef recHitContainer::const_iterator const_iterator;
  typedef std::pair<const_iterator,const_iterator> range;


  // get beam spot
  edm::Handle<reco::BeamSpot> theBeamSpot;
  e.getByType(theBeamSpot);

  // get seeds

  edm::Handle<ElectronSeedCollection> elSeeds;
  e.getByLabel(inputCollection_,elSeeds);
  edm::LogInfo("")<<"\n\n =================> Treating event "<<e.id()<<" Number of seeds "<<elSeeds.product()->size();
  int is=0;

  FTSFromVertexToPointFactory   myFTS;
  float mass=.000511; // electron propagation
  PropagatorWithMaterial* prop1stLayer = new PropagatorWithMaterial(oppositeToMomentum,mass,&(*theMagField));
  PropagatorWithMaterial* prop2ndLayer = new PropagatorWithMaterial(alongMomentum,mass,&(*theMagField));

  float dphi1=0., dphi2=0., drz1=0., drz2=0.;
  float phi1=0., phi2=0., rz1=0., rz2=0.;

  for( ElectronSeedCollection::const_iterator MyS= (*elSeeds).begin(); MyS != (*elSeeds).end(); ++MyS) {

    LogDebug("") <<"\nSeed nr "<<is<<": ";
    range r=(*MyS).recHits();
     LogDebug("")<<" Number of RecHits= "<<(*MyS).nHits();
    const GeomDet *det1=0;const GeomDet *det2=0;

    TrajectorySeed::const_iterator it=r.first;
    DetId id1 = (*it).geographicalId();
    det1 = pDD->idToDet(id1);
    LogDebug("") <<" First hit local x,y,z "<<(*it).localPosition()<<" det "<<id1.det()<<" subdet "<<id1.subdetId();
    LogDebug("") <<" First hit global  "<<det1->toGlobal((*it).localPosition());
    //std::cout <<" First hit local x,y,z "<<(*it).localPosition()<<" det "<<id1.det()<<" subdet "<<id1.subdetId()<< std::endl;
    //std::cout <<" First hit global  "<<det1->toGlobal((*it).localPosition())<< std::endl;
    it++;
    DetId id2 = (*it).geographicalId();
    det2 = pDD->idToDet(id2);
    LogDebug("") <<" Second hit local x,y,z "<<(*it).localPosition()<<" det "<<id2.det()<<" subdet "<<id2.subdetId();
    LogDebug("") <<" Second hit global  "<<det2->toGlobal((*it).localPosition());
    //std::cout <<" Second hit local x,y,z "<<(*it).localPosition()<<" det "<<id2.det()<<" subdet "<<id2.subdetId()<< std::endl;
    //std::cout <<" Second hit global  "<<det2->toGlobal((*it).localPosition()) << std::endl;

    // state on last det
    const GeomDet *det=0;
    for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
    TrajectoryStateOnSurface t=  trajectoryStateTransform::transientState((*MyS).startingState(), &(det->surface()), &(*theMagField));

    // debug

    LogDebug("")<<" ElectronSeed outermost state position: "<<t.globalPosition();
    LogDebug("")<<" ElectronSeed outermost state momentum: "<<t.globalMomentum();
    edm::RefToBase<CaloCluster> caloCluster = (*MyS).caloCluster() ;
    if (caloCluster.isNull()) continue;
    edm::Ref<SuperClusterCollection> theClus = caloCluster.castTo<SuperClusterRef>() ;
    LogDebug("")<<" ElectronSeed superCluster energy: "<<theClus->energy()<<", position: "<<theClus->position();
    LogDebug("")<<" ElectronSeed outermost state Pt: "<<t.globalMomentum().perp();
    LogDebug("")<<" ElectronSeed supercluster Et: "<<theClus->energy()*sin(2.*atan(exp(-theClus->position().eta())));
    LogDebug("")<<" ElectronSeed outermost momentum direction eta: "<<t.globalMomentum().eta();
    LogDebug("")<<" ElectronSeed supercluster eta: "<<theClus->position().eta();
    LogDebug("")<<" ElectronSeed seed charge: "<<(*MyS).getCharge();
    LogDebug("")<<" ElectronSeed E/p: "<<theClus->energy()/t.globalMomentum().mag();

    // retreive SC and compute distances between hit position and prediction the same
    // way as in the PixelHitMatcher
    
    // inputs are charge, cluster position, vertex position, cluster energy and B field
    int charge = int((*MyS).getCharge());
    GlobalPoint xmeas(theClus->position().x(),theClus->position().y(),theClus->position().z());
    GlobalPoint vprim(theBeamSpot->position().x(),theBeamSpot->position().y(),theBeamSpot->position().z());
    float energy = theClus->energy();

    FreeTrajectoryState fts = myFTS(&(*theMagField),xmeas, vprim,
                                 energy, charge);
    //std::cout << "[PixelHitMatcher::compatibleSeeds] fts position, momentum " <<
    // fts.parameters().position() << " " << fts.parameters().momentum() << std::endl;

    PerpendicularBoundPlaneBuilder bpb;
    TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));

    //      TrajectorySeed::range r=(*seeds.product())[i].recHits();
   // TrajectorySeed::range r=(*seeds)[i].recHits();

    // first Hit
    it=r.first;
    DetId id=(*it).geographicalId();
    const GeomDet *geomdet=pDD->idToDet((*it).geographicalId());
    LocalPoint lp=(*it).localPosition();
    GlobalPoint hitPos=geomdet->surface().toGlobal(lp);

    TrajectoryStateOnSurface tsos1;
    tsos1 = prop1stLayer->propagate(tsos,geomdet->surface()) ;

    if (tsos1.isValid()) {

      std::pair<bool,double> est;

      //UB add test on phidiff
      float SCl_phi = xmeas.phi();
      float localDphi = SCl_phi-hitPos.phi();
      if(localDphi>CLHEP::pi)localDphi-=(2*CLHEP::pi);
      if(localDphi<-CLHEP::pi)localDphi+=(2*CLHEP::pi);
      if(std::abs(localDphi)>2.5)continue;

      phi1 = hitPos.phi();
      dphi1 = hitPos.phi() - tsos1.globalPosition().phi();
      rz1 = hitPos.perp();
      drz1 = hitPos.perp() - tsos1.globalPosition().perp();
      if (id.subdetId()%2==1) {
        drz1 = hitPos.z() - tsos1.globalPosition().z();
        rz1 = hitPos.z();
      }

      // now second Hit
      it++;
      DetId id2=(*it).geographicalId();
      const GeomDet *geomdet2=pDD->idToDet((*it).geographicalId());
      TrajectoryStateOnSurface tsos2;

      // compute the z vertex from the cluster point and the found pixel hit
      double pxHit1z = hitPos.z();
      double pxHit1x = hitPos.x();
      double pxHit1y = hitPos.y();
      double r1diff = (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) + (pxHit1y-vprim.y())*(pxHit1y-vprim.y());
      r1diff=sqrt(r1diff);
      double r2diff = (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) + (xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y);
      r2diff=sqrt(r2diff);
      double zVertexPred = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;

      GlobalPoint vertexPred(vprim.x(),vprim.y(),zVertexPred);

      FreeTrajectoryState fts2 = myFTS(&(*theMagField),hitPos,vertexPred,energy, charge);
      tsos2 = prop2ndLayer->propagate(fts2,geomdet2->surface()) ;

      if (tsos2.isValid()) {
        LocalPoint lp2=(*it).localPosition();
        GlobalPoint hitPos2=geomdet2->surface().toGlobal(lp2);
        phi2 = hitPos2.phi();
        dphi2 = hitPos2.phi() - tsos2.globalPosition().phi();
        rz2 = hitPos2.perp();
        drz2 = hitPos2.perp() - tsos2.globalPosition().perp();
        if (id2.subdetId()%2==1) {
          rz2 = hitPos2.z();
          drz2 = hitPos2.z() - tsos2.globalPosition().z();
        }
      }

    }

    // fill the tree and histos

    histpt_->Fill(t.globalMomentum().perp());
    histetclu_->Fill(theClus->energy()*sin(2.*atan(exp(-theClus->position().eta()))));
    histeta_->Fill(t.globalMomentum().eta());
    histetaclu_->Fill(theClus->position().eta());
    histq_->Fill((*MyS).getCharge());
    histeoverp_->Fill(theClus->energy()/t.globalMomentum().mag());

    if (is<10) {
      superclusterEnergy[is] = theClus->energy();
      superclusterEta[is] = theClus->position().eta();
      superclusterPhi[is] = theClus->position().phi();
      superclusterEt[is] = theClus->energy()*sin(2.*atan(exp(-theClus->position().eta())));
      seedMomentum[is] = t.globalMomentum().mag();
      seedEta[is] = t.globalMomentum().eta();
      seedPhi[is] = t.globalMomentum().phi();
      seedPt[is] = t.globalMomentum().perp();
      seedQ[is] = (*MyS).getCharge();
      seedSubdet1[is] = id1.subdetId();
      switch (seedSubdet1[is]) {
        case 1:
          {
          PXBDetId pxbid1( id1.rawId() );
          seedLayer1[is] = pxbid1.layer();
          seedSide1[is] = 0;
          break;
          }
        case 2:
        {
          PXFDetId pxfid1( id1.rawId() );
          seedLayer1[is] = pxfid1.disk();
          seedSide1[is] = pxfid1.side();
          break;
          }
        case 3:
          {
          TIBDetId tibid1( id1.rawId() );
          seedLayer1[is] = tibid1.layer();
          seedSide1[is] = 0;
          break;
          }
        case 4:
          {
          TIDDetId tidid1( id1.rawId() );
          seedLayer1[is] = tidid1.wheel();
          seedSide1[is] = tidid1.side();
          break;
          }
        case 5:
          {
          TOBDetId tobid1( id1.rawId() );
          seedLayer1[is] = tobid1.layer();
          seedSide1[is] = 0;
          break;
          }
        case 6:
          {
          TECDetId tecid1( id1.rawId() );
          seedLayer1[is] = tecid1.wheel();
          seedSide1[is] = tecid1.side();
          break;
          }
      }
      seedPhi1[is] = phi1;
      seedRz1[is] = rz1;
      seedDphi1[is] = dphi1;
      seedDrz1[is] = drz1;
      seedSubdet2[is] = id2.subdetId();
      switch (seedSubdet2[is]) {
        case 1:
          {
          PXBDetId pxbid2( id2.rawId() );
          seedLayer2[is] = pxbid2.layer();
          seedSide2[is] = 0;
          break;
          }
        case 2:
          {
          PXFDetId pxfid2( id2.rawId() );
          seedLayer2[is] = pxfid2.disk();
          seedSide2[is] = pxfid2.side();
          break;
          }
        case 3:
          {
          TIBDetId tibid2( id2.rawId() );
          seedLayer2[is] = tibid2.layer();
          seedSide2[is] = 0;
          break;
          }
        case 4:
          {
          TIDDetId tidid2( id2.rawId() );
          seedLayer2[is] = tidid2.wheel();
          seedSide2[is] = tidid2.side();
          break;
          }
        case 5:
          {
          TOBDetId tobid2( id2.rawId() );
          seedLayer2[is] = tobid2.layer();
          seedSide2[is] = 0;
          break;
          }
        case 6:
          {
          TECDetId tecid2( id2.rawId() );
          seedLayer2[is] = tecid2.wheel();
          seedSide2[is] = tecid2.side();
          break;
          }
      }
      seedDphi2[is] = dphi2;
      seedDrz2[is] = drz2;
      seedPhi2[is] = phi2;
      seedRz2[is] = rz2;
    }

    is++;

  }

  histnbseeds_->Fill(elSeeds.product()->size());

  // get input clusters

  edm::Handle<SuperClusterCollection> clusters;
  //CC to be changed according to supercluster input
  e.getByLabel("correctedHybridSuperClusters", clusters);
  histnbclus_->Fill(clusters.product()->size());
  if (clusters.product()->size()>0) histnrseeds_->Fill(elSeeds.product()->size());
  // get MC information

  edm::Handle<edm::HepMCProduct> HepMCEvt;
  // this one is empty branch in current test files
  //e.getByLabel("VtxSmeared", "", HepMCEvt);
  //e.getByLabel("source", "", HepMCEvt);
  e.getByLabel("generator", "", HepMCEvt);

  const HepMC::GenEvent* MCEvt = HepMCEvt->GetEvent();
  HepMC::GenParticle* genPc=0;
  HepMC::FourVector pAssSim;
  int ip=0;
  for (HepMC::GenEvent::particle_const_iterator partIter = MCEvt->particles_begin();
   partIter != MCEvt->particles_end(); ++partIter) {

    for (HepMC::GenEvent::vertex_const_iterator vertIter = MCEvt->vertices_begin();
     vertIter != MCEvt->vertices_end(); ++vertIter) {

      //      CLHEP::HepLorentzVector creation = (*partIter)->CreationVertex();
      HepMC::GenVertex * creation = (*partIter)->production_vertex();
      //      CLHEP::HepLorentzVector momentum = (*partIter)->Momentum();
      HepMC::FourVector momentum = (*partIter)->momentum();
      //      HepPDT::ParticleID id = (*partIter)->particleID();  // electrons and positrons are 11 and -11
       int id = (*partIter)->pdg_id();  // electrons and positrons are 11 and -11
     LogDebug("")  << "MC particle id " << id << ", creationVertex " << (*creation) << " cm, initialMomentum " << momentum.rho() << " GeV/c" << std::endl;

      if (id == 11 || id == -11) {

      // single primary electrons or electrons from Zs or Ws
      HepMC::GenParticle* mother = 0;
      if ( (*partIter)->production_vertex() )  {
       if ( (*partIter)->production_vertex()->particles_begin(HepMC::parents) !=
           (*partIter)->production_vertex()->particles_end(HepMC::parents))
            mother = *((*partIter)->production_vertex()->particles_begin(HepMC::parents));
      }
      if ( ((mother == 0) || ((mother != 0) && (mother->pdg_id() == 23))
                          || ((mother != 0) && (mother->pdg_id() == 32))
                          || ((mother != 0) && (std::abs(mother->pdg_id()) == 24)))) {
      genPc=(*partIter);
      pAssSim = genPc->momentum();

      // EWK fiducial
      //if (pAssSim.perp()> 100. || std::abs(pAssSim.eta())> 2.5) continue;     
      //if (pAssSim.perp()< 20. || (std::abs(pAssSim.eta())> 1.4442 && std::abs(pAssSim.eta())< 1.56) || std::abs(pAssSim.eta())> 2.5) continue;
      // reconstruction fiducial
      //if (pAssSim.perp()< 5. || std::abs(pAssSim.eta())> 2.5) continue;
      if (std::abs(pAssSim.eta())> 2.5) continue;

      histptMC_->Fill(pAssSim.perp());
      histetaMC_->Fill(pAssSim.eta());
      histeMC_->Fill(pAssSim.rho());

     // looking for the best matching gsf electron
      bool okSeedFound = false;
      double seedOkRatio = 999999.;

      // find best matched seed
      reco::ElectronSeed bestElectronSeed;
      for( ElectronSeedCollection::const_iterator gsfIter= (*elSeeds).begin(); gsfIter != (*elSeeds).end(); ++gsfIter) {

        range r=gsfIter->recHits();
        const GeomDet *det=0;
        for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
         TrajectoryStateOnSurface t= trajectoryStateTransform::transientState(gsfIter->startingState(), &(det->surface()), &(*theMagField));

        float eta = t.globalMomentum().eta();
        float phi = t.globalMomentum().phi();
        float p = t.globalMomentum().mag();
        double dphi = phi-pAssSim.phi();
        if (std::abs(dphi)>CLHEP::pi)
         dphi = dphi < 0? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
        double deltaR = sqrt(std::pow((eta-pAssSim.eta()),2) + std::pow(dphi,2));
        if ( deltaR < 0.15 ){
//      if ( deltaR < 0.3 ){
        //if ( (genPc->pdg_id() == 11) && (gsfIter->charge() < 0.) || (genPc->pdg_id() == -11) &&
        //(gsfIter->charge() > 0.) ){
          double tmpSeedRatio = p/pAssSim.t();
          if ( std::abs(tmpSeedRatio-1) < std::abs(seedOkRatio-1) ) {
            seedOkRatio = tmpSeedRatio;
            bestElectronSeed=*gsfIter;
            okSeedFound = true;
          }
        //}
        }
      } // loop over rec ele to look for the best one

      // analysis when the mc track is found
     if (okSeedFound){

        histptMCmatched_->Fill(pAssSim.perp());
        histetaMCmatched_->Fill(pAssSim.eta());
        histeMCmatched_->Fill(pAssSim.rho());
        if (ip<10) {
          mcEnergy[ip] = pAssSim.rho();
          mcEta[ip] = pAssSim.eta();
          mcPhi[ip] = pAssSim.phi();
          mcPt[ip] = pAssSim.perp();
          mcQ[ip] = ((id == 11) ? -1.: +1.);
        }
      }
      
     // efficiency for ecal driven only
      okSeedFound = false;
      seedOkRatio = 999999.;

      // find best matched seed
      for( ElectronSeedCollection::const_iterator gsfIter= (*elSeeds).begin(); gsfIter != (*elSeeds).end(); ++gsfIter) {

        range r=gsfIter->recHits();
        const GeomDet *det=0;
        for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
         TrajectoryStateOnSurface t= trajectoryStateTransform::transientState(gsfIter->startingState(), &(det->surface()), &(*theMagField));

        float eta = t.globalMomentum().eta();
        float phi = t.globalMomentum().phi();
        float p = t.globalMomentum().mag();
        double dphi = phi-pAssSim.phi();
        if (std::abs(dphi)>CLHEP::pi)
         dphi = dphi < 0? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
        double deltaR = sqrt(std::pow((eta-pAssSim.eta()),2) + std::pow(dphi,2));
        if (gsfIter->isEcalDriven()) {
        if ( deltaR < 0.15 ){
//      if ( deltaR < 0.3 ){
        //if ( (genPc->pdg_id() == 11) && (gsfIter->charge() < 0.) || (genPc->pdg_id() == -11) &&
        //(gsfIter->charge() > 0.) ){
          double tmpSeedRatio = p/pAssSim.t();
          if ( std::abs(tmpSeedRatio-1) < std::abs(seedOkRatio-1) ) {
            seedOkRatio = tmpSeedRatio;
            bestElectronSeed=*gsfIter;
            okSeedFound = true;
          }
        //}
        }
        } // end if ecal driven
      } // loop over rec ele to look for the best one

      // analysis when the mc track is found
     if (okSeedFound){

        histecaldrivenptMCmatched_->Fill(pAssSim.perp());
        histecaldrivenetaMCmatched_->Fill(pAssSim.eta());
        histecaldriveneMCmatched_->Fill(pAssSim.rho());
      }

      // efficiency for tracker driven only
      okSeedFound = false;
      seedOkRatio = 999999.;

      // find best matched seed
      for( ElectronSeedCollection::const_iterator gsfIter= (*elSeeds).begin(); gsfIter != (*elSeeds).end(); ++gsfIter) {

        range r=gsfIter->recHits();
        const GeomDet *det=0;
        for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
         TrajectoryStateOnSurface t= trajectoryStateTransform::transientState(gsfIter->startingState(), &(det->surface()), &(*theMagField));

        float eta = t.globalMomentum().eta();
        float phi = t.globalMomentum().phi();
        float p = t.globalMomentum().mag();
        double dphi = phi-pAssSim.phi();
        if (std::abs(dphi)>CLHEP::pi)
         dphi = dphi < 0? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
        double deltaR = sqrt(std::pow((eta-pAssSim.eta()),2) + std::pow(dphi,2));
        if (gsfIter->isTrackerDriven()) {
        if ( deltaR < 0.15 ){
//      if ( deltaR < 0.3 ){
        //if ( (genPc->pdg_id() == 11) && (gsfIter->charge() < 0.) || (genPc->pdg_id() == -11) &&
        //(gsfIter->charge() > 0.) ){
          double tmpSeedRatio = p/pAssSim.t();
          if ( std::abs(tmpSeedRatio-1) < std::abs(seedOkRatio-1) ) {
            seedOkRatio = tmpSeedRatio;
            bestElectronSeed=*gsfIter;
            okSeedFound = true;
          }
        //}
        }
        } // end if ecal driven
      } // loop over rec ele to look for the best one

      // analysis when the mc track is found
     if (okSeedFound){

        histtrackerdrivenptMCmatched_->Fill(pAssSim.perp());
        histtrackerdrivenetaMCmatched_->Fill(pAssSim.eta());
        histtrackerdriveneMCmatched_->Fill(pAssSim.rho());
      }

     } // end if mother W or Z

    } // end if gen part is electron

    } // end loop on vertices

    ip++;

  } // end loop on gen particles

  //tree_->Fill();

}

void ElectronSeedAnalyzer::beginJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 68 of file ElectronSeedAnalyzer.cc.

{
  histfile_->cd();
  tree_ = new TTree("ElectronSeeds","ElectronSeed validation ntuple");
  tree_->Branch("mcEnergy",mcEnergy,"mcEnergy[10]/F");
  tree_->Branch("mcEta",mcEta,"mcEta[10]/F");
  tree_->Branch("mcPhi",mcPhi,"mcPhi[10]/F");
  tree_->Branch("mcPt",mcPt,"mcPt[10]/F");
  tree_->Branch("mcQ",mcQ,"mcQ[10]/F");
  tree_->Branch("superclusterEnergy",superclusterEnergy,"superclusterEnergy[10]/F");
  tree_->Branch("superclusterEta",superclusterEta,"superclusterEta[10]/F");
  tree_->Branch("superclusterPhi",superclusterPhi,"superclusterPhi[10]/F");
  tree_->Branch("superclusterEt",superclusterEt,"superclusterEt[10]/F");
  tree_->Branch("seedMomentum",seedMomentum,"seedMomentum[10]/F");
  tree_->Branch("seedEta",seedEta,"seedEta[10]/F");
  tree_->Branch("seedPhi",seedPhi,"seedPhi[10]/F");
  tree_->Branch("seedPt",seedPt,"seedPt[10]/F");
  tree_->Branch("seedQ",seedQ,"seedQ[10]/F");
  tree_->Branch("seedSubdet1",seedSubdet1,"seedSubdet1[10]/I");
  tree_->Branch("seedLayer1",seedLayer1,"seedLayer1[10]/I");
  tree_->Branch("seedSide1",seedSide1,"seedSide1[10]/I");
  tree_->Branch("seedPhi1",seedPhi1,"seedPhi1[10]/F");
  tree_->Branch("seedDphi1",seedDphi1,"seedDphi1[10]/F");
  tree_->Branch("seedDrz1",seedDrz1,"seedDrz1[10]/F");
  tree_->Branch("seedRz1",seedRz1,"seedRz1[10]/F");
  tree_->Branch("seedSubdet2",seedSubdet2,"seedSubdet2[10]/I");
  tree_->Branch("seedLayer2",seedLayer2,"seedLayer2[10]/I");
  tree_->Branch("seedSide2",seedSide2,"seedSide2[10]/I");
  tree_->Branch("seedPhi2",seedPhi2,"seedPhi2[10]/F");
  tree_->Branch("seedDphi2",seedDphi2,"seedDphi2[10]/F");
  tree_->Branch("seedRz2",seedRz2,"seedRz2[10]/F");
  tree_->Branch("seedDrz2",seedDrz2,"seedDrz2[10]/F");
  histeMC_ = new TH1F("eMC","MC particle energy",100,0.,100.);
  histeMCmatched_ = new TH1F("eMCmatched","matched MC particle energy",100,0.,100.);
  histecaldriveneMCmatched_ = new TH1F("ecaldriveneMCmatched","matched MC particle energy, ecal driven",100,0.,100.);
  histtrackerdriveneMCmatched_ = new TH1F("trackerdriveneMCmatched","matched MC particle energy, tracker driven",100,0.,100.);
  histp_ = new TH1F("p","seed p",100,0.,100.);
  histeclu_ = new TH1F("clus energy","supercluster energy",100,0.,100.);
  histpt_ = new TH1F("pt","seed pt",100,0.,100.);
  histptMC_ = new TH1F("ptMC","MC particle pt",100,0.,100.);
  histptMCmatched_ = new TH1F("ptMCmatched","matched MC particle pt",100,0.,100.);
  histecaldrivenptMCmatched_ = new TH1F("ecaldrivenptMCmatched","matched MC particle pt, ecal driven",100,0.,100.);
  histtrackerdrivenptMCmatched_ = new TH1F("trackerdrivenptMCmatched","matched MC particle pt, tracker driven",100,0.,100.);
  histetclu_ = new TH1F("Et","supercluster Et",100,0.,100.);
  histeffpt_ = new TH1F("pt eff","seed effciency vs pt",100,0.,100.);
  histeta_ = new TH1F("seed eta","seed eta",100,-2.5,2.5);
  histetaMC_ = new TH1F("etaMC","MC particle eta",100,-2.5,2.5);
  histetaMCmatched_ = new TH1F("etaMCmatched","matched MC particle eta",100,-2.5,2.5);
  histecaldrivenetaMCmatched_ = new TH1F("ecaldrivenetaMCmatched","matched MC particle eta, ecal driven",100,-2.5,2.5);
  histtrackerdrivenetaMCmatched_ = new TH1F("trackerdrivenetaMCmatched","matched MC particle eta, tracker driven",100,-2.5,2.5);
  histetaclu_ = new TH1F("clus eta","supercluster eta",100,-2.5,2.5);
  histeffeta_ = new TH1F("eta eff","seed effciency vs eta",100,-2.5,2.5);
  histq_ = new TH1F("q","seed charge",100,-2.5,2.5);
  histeoverp_ = new TH1F("E/p","seed E/p",100,0.,10.);
  histnbseeds_ = new TH1I("nrs","Nr of seeds ",50,0.,25.);
  histnbclus_ = new TH1I("nrclus","Nr of superclusters ",50,0.,25.);
  histnrseeds_ = new TH1I("ns","Nr of seeds if clusters",50,0.,25.);
}

void ElectronSeedAnalyzer::endJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 127 of file ElectronSeedAnalyzer.cc.

{
  histfile_->cd();
  tree_->Print();
  tree_->Write();

  // efficiency vs eta
  TH1F *histetaEff = (TH1F*)histetaMCmatched_->Clone("histetaEff");
  histetaEff->Reset();
  histetaEff->Divide(histetaMCmatched_,histeta_,1,1,"b");
  histetaEff->Print();
  histetaEff->GetXaxis()->SetTitle("#eta");
  histetaEff->GetYaxis()->SetTitle("Efficiency");

  // efficiency vs pt
  TH1F *histptEff = (TH1F*)histptMCmatched_->Clone("histotEff");
  histptEff->Reset();
  histptEff->Divide(histptMCmatched_,histpt_,1,1,"b");
  histptEff->Print();
  histptEff->GetXaxis()->SetTitle("p_{T}");
  histptEff->GetYaxis()->SetTitle("Efficiency");

  histeMCmatched_->Write();
  histecaldriveneMCmatched_->Write();
  histtrackerdriveneMCmatched_->Write();
  histeMC_->Write();
  histp_->Write();
  histeclu_->Write();
  histpt_->Write();
  histptMCmatched_->Write();
  histecaldrivenptMCmatched_->Write();
  histtrackerdrivenptMCmatched_->Write();
  histptMC_->Write();
  histetclu_->Write();
  histeffpt_->Write();
  histeta_->Write();
  histetaMCmatched_->Write();
  histecaldrivenetaMCmatched_->Write();
  histtrackerdrivenetaMCmatched_->Write();
  histetaMC_->Write();
  histetaclu_->Write();
  histeffeta_->Write();
  histq_->Write();
  histeoverp_->Write();
  histnbseeds_->Write();
  histnbclus_->Write();
  histnrseeds_->Write();
}