#include <GoodSeedProducer.h>

Inheritance diagram for GoodSeedProducer:

Public Member Functions
	GoodSeedProducer (const edm::ParameterSet &)
	~GoodSeedProducer ()
Private Types
typedef TrajectoryStateOnSurface	TSOS
Private Member Functions
virtual void	beginRun (edm::Run &run, const edm::EventSetup &)
virtual void	endRun ()
void	fillPreIdRefValueMap (edm::Handle< reco::TrackCollection > tkhandle, const edm::OrphanHandle< reco::PreIdCollection > &, edm::ValueMap< reco::PreIdRef >::Filler &filler)
int	getBin (float, float)
	Find the bin in pt and eta.
int	getBin (float)
bool	IsIsolated (float charge, float P, math::XYZPointF, const reco::PFClusterCollection &ecalColl, const reco::PFClusterCollection &hcalColl)
virtual void	produce (edm::Event &, const edm::EventSetup &)
void	PSforTMVA (math::XYZTLorentzVector mom, math::XYZTLorentzVector pos)
Private Attributes
bool	applyIsolation_
	ISOLATION REQUEST AS DONE IN THE TAU GROUP.
math::XYZVector	B_
	B field.
float	chi
float	chiRatio
float	chired
double	clusThreshold_
	Cut on the energy of the clusters.
edm::ParameterSet	conf_
bool	disablePreId_
	switch to disable the pre-id
float	dpt
double	EcalStripSumE_deltaEta_
double	EcalStripSumE_deltaPhiOverQ_maxValue_
double	EcalStripSumE_deltaPhiOverQ_minValue_
double	EcalStripSumE_minClusEnergy_
float	eP
	VARIABLES NEEDED FOR TMVA.
float	eta
edm::ESHandle< TrajectoryFitter >	fitter_
	Fitter.
std::string	fitterName_
double	HcalIsolWindow_
double	maxEp_
double	maxEta_
double	maxHoverP_
double	maxPt_
std::string	method_
	TMVA method.
double	minEoverP_
double	minEp_
	Min and MAx allowed values forEoverP.
double	minPt_
	Minimum transverse momentum and maximum pseudorapidity.
float	nhit
int	nHitsInSeed_
	Number of hits in the seed;.
edm::InputTag	pfCLusTagECLabel_
edm::InputTag	pfCLusTagHCLabel_
edm::InputTag	pfCLusTagPSLabel_
PFTrackTransformer *	pfTransformer_
	PFTrackTransformer.
std::string	preidckf_
	Name of the Seed(Ckf) Collection.
std::string	preidgsf_
	Name of the Seed(Gsf) Collection.
std::string	preidname_
	Name of the preid Collection (FB)
bool	produceCkfseed_
	Produce the Seed for Ckf tracks?
bool	producePreId_
	Produce the pre-id debugging collection.
std::string	propagatorName_
float	ps1chi
std::vector< reco::PFCluster >	ps1Clus
	Vector of clusters of the PreShower.
float	ps1En
float	ps2chi
std::vector< reco::PFCluster >	ps2Clus
float	ps2En
float	pt
double	PtThresholdSavePredId_
	Threshold to save Pre Idinfo.
TMVA::Reader *	reader
	READER FOR TMVA.
std::map< reco::TrackRef, unsigned >	refMap_
	Map used to create the TrackRef, PreIdRef value map.
edm::ESHandle< TrajectorySmoother >	smoother_
	Smoother.
std::string	smootherName_
float	thr [150]
	vector of thresholds for different bins of eta and pt
float	thrPS [20]
reco::TrackBase::TrackQuality	trackQuality_
std::vector< edm::InputTag >	tracksContainers_
bool	usePreshower_
	Use of Preshower clusters.
bool	useQuality_
	TRACK QUALITY.
bool	useTmva_
	USE OF TMVA.
Static Private Attributes
static PFResolutionMap *	resMapEtaECAL_ = 0
static PFResolutionMap *	resMapPhiECAL_ = 0

Detailed Description

Definition at line 46 of file GoodSeedProducer.h.

Member Typedef Documentation

typedef TrajectoryStateOnSurface GoodSeedProducer::TSOS [private]

Definition at line 47 of file GoodSeedProducer.h.

Constructor & Destructor Documentation

GoodSeedProducer::GoodSeedProducer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 40 of file GoodSeedProducer.cc.

References applyIsolation_, clusThreshold_, disablePreId_, EcalStripSumE_deltaEta_, EcalStripSumE_deltaPhiOverQ_maxValue_, EcalStripSumE_deltaPhiOverQ_minValue_, EcalStripSumE_minClusEnergy_, fitterName_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), HcalIsolWindow_, LogDebug, maxEp_, maxEta_, maxHoverP_, maxPt_, minEoverP_, minEp_, minPt_, nHitsInSeed_, pfCLusTagECLabel_, pfCLusTagHCLabel_, pfCLusTagPSLabel_, preidckf_, preidgsf_, preidname_, produceCkfseed_, producePreId_, edm::ProductRegistryHelper::produces(), PtThresholdSavePredId_, smootherName_, trackQuality_, tracksContainers_, usePreshower_, useQuality_, and useTmva_.

                                                             :
  pfTransformer_(0),
  conf_(iConfig)
{
  LogInfo("GoodSeedProducer")<<"Electron PreIdentification started  ";
  
  //now do what ever initialization is needed
 
  tracksContainers_ = 
    iConfig.getParameter< vector < InputTag > >("TkColList");
  
  minPt_=iConfig.getParameter<double>("MinPt");
  maxPt_=iConfig.getParameter<double>("MaxPt");
  maxEta_=iConfig.getParameter<double>("MaxEta");


  //ISOLATION REQUEST AS DONE IN THE TAU GROUP
  applyIsolation_ =iConfig.getParameter<bool>("ApplyIsolation");
  HcalIsolWindow_                       =iConfig.getParameter<double>("HcalWindow");
  EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
  EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
  EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
  EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
   minEoverP_= iConfig.getParameter<double>("EOverPLead_minValue");
   maxHoverP_= iConfig.getParameter<double>("HOverPLead_maxValue");
 
  //
  pfCLusTagECLabel_=
    iConfig.getParameter<InputTag>("PFEcalClusterLabel");

  pfCLusTagHCLabel_=
   iConfig.getParameter<InputTag>("PFHcalClusterLabel");  

  pfCLusTagPSLabel_=
    iConfig.getParameter<InputTag>("PFPSClusterLabel");
  
  preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
  preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
  preidname_= iConfig.getParameter<string>("PreIdLabel");
  
  
  fitterName_ = iConfig.getParameter<string>("Fitter");
  smootherName_ = iConfig.getParameter<string>("Smoother");
  
  
  nHitsInSeed_=iConfig.getParameter<int>("NHitsInSeed");

  clusThreshold_=iConfig.getParameter<double>("ClusterThreshold");
  
  minEp_=iConfig.getParameter<double>("MinEOverP");
  maxEp_=iConfig.getParameter<double>("MaxEOverP");

  //collection to produce
  produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed",false);

  // to disable the electron part (for HI collisions for examples) 
  disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId",false);  

  producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId",true);  
  //  if no electron, cannot produce the preid
  if(disablePreId_) 
    producePreId_=false;
  PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId",1.);  

  LogDebug("GoodSeedProducer")<<"Seeds for GSF will be produced ";

  // no disablePreId_ switch here. The collection will be empty if it is true
  produces<ElectronSeedCollection>(preidgsf_);

  if(produceCkfseed_){
    LogDebug("GoodSeedProducer")<<"Seeds for CKF will be produced ";
    produces<TrajectorySeedCollection>(preidckf_);
  }
  

  if(producePreId_){
    LogDebug("GoodSeedProducer")<<"PreId debugging information will be produced ";

    produces<PreIdCollection>(preidname_);
    if(tracksContainers_.size()==1) // do not make a value map if more than one input track collection
      produces<edm::ValueMap<reco::PreIdRef> >(preidname_);
  } 
  
  useQuality_   = iConfig.getParameter<bool>("UseQuality");
  trackQuality_=TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));

  useTmva_= iConfig.getUntrackedParameter<bool>("UseTMVA",false);
  
  usePreshower_ = iConfig.getParameter<bool>("UsePreShower");

}

GoodSeedProducer::~GoodSeedProducer ( )

Definition at line 133 of file GoodSeedProducer.cc.

References pfTransformer_.

{
  
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.) 

  delete pfTransformer_;
}

Member Function Documentation

void GoodSeedProducer::beginRun	(	edm::Run &	run,
		const edm::EventSetup &	es
	)		`[private, virtual]`

Reimplemented from edm::EDProducer.

Definition at line 493 of file GoodSeedProducer.cc.

References B_, chi, chiRatio, chired, conf_, dpt, eP, eta, edm::FileInPath::fullPath(), edm::EventSetup::get(), edm::ParameterSet::getParameter(), method_, nhit, PFTrackTransformer::OnlyProp(), pfTransformer_, pt, reader, resMapEtaECAL_, resMapPhiECAL_, thr, thrPS, and useTmva_.

{
  //Magnetic Field
  ESHandle<MagneticField> magneticField;
  es.get<IdealMagneticFieldRecord>().get(magneticField);
  B_=magneticField->inTesla(GlobalPoint(0,0,0));
  
  pfTransformer_= new PFTrackTransformer(B_);
  pfTransformer_->OnlyProp();


  
  //Resolution maps
  FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
  FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
  resMapEtaECAL_ = new PFResolutionMap("ECAL_eta",ecalEtaMap.fullPath().c_str());
  resMapPhiECAL_ = new PFResolutionMap("ECAL_phi",ecalPhiMap.fullPath().c_str());

  if(useTmva_){
    reader = new TMVA::Reader("!Color:Silent");
    method_ = conf_.getParameter<string>("TMVAMethod");
    
    reader->AddVariable("eP",&eP);
    reader->AddVariable("chi",&chi);
    reader->AddVariable("chired",&chired);
    reader->AddVariable("chiRatio",&chiRatio);
    reader->AddVariable("dpt",&dpt);
    reader->AddVariable("nhit",&nhit);
    reader->AddVariable("eta",&eta);
    reader->AddVariable("pt",&pt);
    FileInPath Weigths(conf_.getParameter<string>("Weights"));
    reader->BookMVA( method_, Weigths.fullPath().c_str()  );    
    }

    
    //read threshold
    FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
    ifstream ifs(parFile.fullPath().c_str());
    for (int iy=0;iy<72;iy++) ifs >> thr[iy];
    
    //read PS threshold
    FileInPath parPSFile(conf_.getParameter<string>("PSThresholdFile"));
    ifstream ifsPS(parPSFile.fullPath().c_str());
    for (int iy=0;iy<12;iy++) ifsPS >> thrPS[iy];

}

void GoodSeedProducer::endRun ( void ) [private, virtual]

Definition at line 542 of file GoodSeedProducer.cc.

References pfTransformer_, reader, resMapEtaECAL_, resMapPhiECAL_, and useTmva_.

                         {
  delete pfTransformer_;
  delete resMapEtaECAL_;
  delete resMapPhiECAL_;
  if(useTmva_) delete reader;
}

void GoodSeedProducer::fillPreIdRefValueMap	(	edm::Handle< reco::TrackCollection >	tkhandle,
		const edm::OrphanHandle< reco::PreIdCollection > &	preidhandle,
		edm::ValueMap< reco::PreIdRef >::Filler &	filler
	)		`[private]`

Definition at line 689 of file GoodSeedProducer.cc.

References edm::helper::Filler< Map >::insert(), refMap_, and makeHLTPrescaleTable::values.

Referenced by produce().

{
  std::vector<reco::PreIdRef> values;

  unsigned ntracks=tracks->size();
  for(unsigned itrack=0;itrack<ntracks;++itrack)
   {
     reco::TrackRef theTrackRef(tracks,itrack);
     std::map<reco::TrackRef,unsigned>::const_iterator itcheck=refMap_.find(theTrackRef);
     if(itcheck==refMap_.end()) 
       {
         // the track has been early discarded
         values.push_back(reco::PreIdRef());
       }
     else
       {
         edm::Ref<reco::PreIdCollection> preIdRef(preidhandle,itcheck->second);
         values.push_back(preIdRef);
         //      std::cout << " Checking Refs " << (theTrackRef==preIdRef->trackRef()) << std::endl;
       }
   }
  filler.insert(tracks,values.begin(),values.end());
}

int GoodSeedProducer::getBin	(	float	eta,
		float	pt
	)		`[private]`

Find the bin in pt and eta.

Definition at line 560 of file GoodSeedProducer.cc.

References LogDebug.

Referenced by produce().

                                           {
  int ie=0;
  int ip=0;
  if (fabs(eta)<1.2) ie=0;
  else{ if (fabs(eta)<1.68) ie=1;
    else ie=2;
  }
  if (pt<6) ip=0;
  else {  if (pt<12) ip=1;     
        else ip=2;
  }
  int iep= ie*3+ip;
  LogDebug("GoodSeedProducer")<<"Track pt ="<<pt<<" eta="<<eta<<" bin="<<iep;
  return iep;
}

int GoodSeedProducer::getBin ( float pt ) [private]

Definition at line 550 of file GoodSeedProducer.cc.

                                {
int ip=0;
  if (pt<6) ip=0;
  else {  if (pt<12) ip=1;
        else ip=2;
  }
return ip;
}

bool GoodSeedProducer::IsIsolated	(	float	charge,
		float	P,
		math::XYZPointF	myElecTrkEcalPos,
		const reco::PFClusterCollection &	ecalColl,
		const reco::PFClusterCollection &	hcalColl
	)		`[private]`

Definition at line 637 of file GoodSeedProducer.cc.

References abs, DeDxDiscriminatorTools::charge(), Geom::deltaPhi(), deltaR(), EcalStripSumE_deltaPhiOverQ_maxValue_, EcalStripSumE_deltaPhiOverQ_minValue_, EcalStripSumE_minClusEnergy_, HcalIsolWindow_, maxHoverP_, minEoverP_, and P.

Referenced by produce().

                                                                 {


  double myHCALenergy3x3=0.;
  double myStripClusterE=0.;
 

  //  reco::TrackRef myElecTrk;
  
  if (fabs(myElecTrkEcalPos.z())<1. && myElecTrkEcalPos.x()<1. && myElecTrkEcalPos.y()<1. ) return false; 

  
  
  PFClusterCollection::const_iterator hc=hcalColl.begin();
  PFClusterCollection::const_iterator hcend=hcalColl.end();
  for (;hc!=hcend;++hc){
    math::XYZPoint clusPos = hc->position();
    double en = hc->energy();
    double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos,clusPos);
    if (deltaR<HcalIsolWindow_) {
      myHCALenergy3x3 += en;
      
    }
  }



  PFClusterCollection::const_iterator ec=ecalColl.begin();
  PFClusterCollection::const_iterator ecend=ecalColl.end();
  for (;ec!=ecend;++ec){
    math::XYZPoint clusPos = ec->position();
    double en = ec->energy();


    double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos,clusPos);
    double deltaEta = abs(myElecTrkEcalPos.eta()-clusPos.eta());
    double deltaPhiOverQ = deltaPhi/charge;
    if (en >= EcalStripSumE_minClusEnergy_ && deltaEta<EcalStripSumE_deltaEta_ && deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ && deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) { 
      myStripClusterE += en;
    }
  }       
  
  double EoP=myStripClusterE/P;
  double HoP=myHCALenergy3x3/P;


  return ((EoP>minEoverP_)&&(EoP<2.5) && (HoP<maxHoverP_))?true:false;
}

void GoodSeedProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[private, virtual]`

Implements edm::EDProducer.

Definition at line 149 of file GoodSeedProducer.cc.

References applyIsolation_, B_, DeDxDiscriminatorTools::charge(), chi, chiRatio, chired, clusThreshold_, disablePreId_, dpt, eP, eta, edm::EventID::event(), edm::helper::Filler< Map >::fill(), fillPreIdRefValueMap(), PFResolutionMap::FindBin(), fitter_, fitterName_, edm::EventSetup::get(), getBin(), edm::Event::getByLabel(), BaseParticlePropagator::getSuccess(), i, edm::EventBase::id(), IsIsolated(), LogDebug, maxEp_, maxEta_, maxPt_, method_, minEp_, minPt_, RawParticle::momentum(), nhit, L1TEmulatorMonitor_cff::p, pfCLusTagECLabel_, pfCLusTagHCLabel_, pfCLusTagPSLabel_, pos, funct::pow(), preidckf_, preidgsf_, preidname_, produceCkfseed_, producePreId_, edm::ESHandle< T >::product(), edm::Handle< T >::product(), BaseParticlePropagator::propagateToEcalEntrance(), ps1chi, ps1Clus, ps1En, ps2chi, ps2Clus, ps2En, PSforTMVA(), pt, reco::PreId::pt(), PtThresholdSavePredId_, edm::Event::put(), reader, refMap_, resMapEtaECAL_, resMapPhiECAL_, edm::EventID::run(), RawParticle::setCharge(), reco::ElectronSeed::setCtfTrack(), reco::PreId::setECALMatching(), reco::PreId::setECALMatchingProperties(), reco::PreId::setESMatching(), reco::PreId::setFinalDecision(), reco::PreId::setMVA(), reco::PreId::setTrack(), reco::PreId::setTrackFiltering(), reco::PreId::setTrackProperties(), smoother_, smootherName_, mathSSE::sqrt(), thr, thrPS, tmp, trackQuality_, tracksContainers_, TwoPi, usePreshower_, useQuality_, useTmva_, RawParticle::vertex(), ExpressReco_HICollisions_FallBack::x, ExpressReco_HICollisions_FallBack::y, and z.

{

  LogDebug("GoodSeedProducer")<<"START event: "<<iEvent.id().event()
                              <<" in run "<<iEvent.id().run();
  //Create empty output collections
  auto_ptr<ElectronSeedCollection> output_preid(new ElectronSeedCollection);
  auto_ptr<TrajectorySeedCollection> output_nopre(new TrajectorySeedCollection);
  auto_ptr<PreIdCollection> output_preidinfo(new PreIdCollection);
  auto_ptr<edm::ValueMap<reco::PreIdRef> > preIdMap_p(new edm::ValueMap<reco::PreIdRef>);
  edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);

  //Tracking Tools
  if(!disablePreId_)
    {
      iSetup.get<TrajectoryFitter::Record>().get(fitterName_, fitter_);
      iSetup.get<TrajectoryFitter::Record>().get(smootherName_, smoother_);
    }

  // clear temporary maps
  refMap_.clear();

  //Handle input collections
  //ECAL clusters             
  Handle<PFClusterCollection> theECPfClustCollection;
  iEvent.getByLabel(pfCLusTagECLabel_,theECPfClustCollection);
  
  vector<PFCluster> basClus;
  vector<PFCluster>::const_iterator iklus;
  for (iklus=theECPfClustCollection.product()->begin();
       iklus!=theECPfClustCollection.product()->end();
       iklus++){
    if((*iklus).energy()>clusThreshold_) basClus.push_back(*iklus);
  }

  //HCAL clusters
  Handle<PFClusterCollection> theHCPfClustCollection;
  iEvent.getByLabel(pfCLusTagHCLabel_,theHCPfClustCollection);
  
  //PS clusters
  Handle<PFClusterCollection> thePSPfClustCollection;
  iEvent.getByLabel(pfCLusTagPSLabel_,thePSPfClustCollection);
  
  ps1Clus.clear();
  ps2Clus.clear();
  
  for (iklus=thePSPfClustCollection.product()->begin();
       iklus!=thePSPfClustCollection.product()->end();
       iklus++){
    //layer==-11 first layer of PS
    //layer==-12 secon layer of PS
    if ((*iklus).layer()==-11) ps1Clus.push_back(*iklus);
    if ((*iklus).layer()==-12) ps2Clus.push_back(*iklus);
  }
  
  //Vector of track collections
  for (unsigned int istr=0; istr<tracksContainers_.size();istr++){
    
    //Track collection
    Handle<TrackCollection> tkRefCollection;
    iEvent.getByLabel(tracksContainers_[istr], tkRefCollection);
    const TrackCollection&  Tk=*(tkRefCollection.product());
    
    //Trajectory collection
    Handle<vector<Trajectory> > tjCollection;
    iEvent.getByLabel(tracksContainers_[istr], tjCollection);
    vector<Trajectory> Tj=*(tjCollection.product());
    
    
    LogDebug("GoodSeedProducer")<<"Number of tracks in collection "
                                <<tracksContainers_[istr] <<" to be analyzed "
                                <<Tj.size();
    

    //loop over the track collection
    for(unsigned int i=0;i<Tk.size();i++){              
      if (useQuality_ &&
          (!(Tk[i].quality(trackQuality_)))) continue;
      
      reco::PreId myPreId;
      bool GoodPreId=false;

      TrackRef trackRef(tkRefCollection, i);
      // TrajectorySeed Seed=Tj[i].seed();
      TrajectorySeed Seed=(*trackRef->seedRef());
      if(!disablePreId_)
        {
          int ipteta=getBin(Tk[i].eta(),Tk[i].pt());
          int ibin=ipteta*8;
          
          float PTOB=Tj[i].lastMeasurement().updatedState().globalMomentum().mag();
          float chikfred=Tk[i].normalizedChi2();
          int nhitpi=Tj[i].foundHits();
          float EP=0;
      
          // set track info
          myPreId.setTrack(trackRef);
          //CLUSTERS - TRACK matching
      
          float pfmass=  0.0005;
          float pfoutenergy=sqrt((pfmass*pfmass)+Tk[i].outerMomentum().Mag2());
          XYZTLorentzVector mom =XYZTLorentzVector(Tk[i].outerMomentum().x(),
                                                   Tk[i].outerMomentum().y(),
                                                   Tk[i].outerMomentum().z(),
                                                   pfoutenergy);
          XYZTLorentzVector pos =   XYZTLorentzVector(Tk[i].outerPosition().x(),
                                                      Tk[i].outerPosition().y(),
                                                      Tk[i].outerPosition().z(),
                                                      0.);

          BaseParticlePropagator theOutParticle = 
            BaseParticlePropagator( RawParticle(mom,pos),
                                    0,0,B_.z());
          theOutParticle.setCharge(Tk[i].charge());
      
          theOutParticle.propagateToEcalEntrance(false);
      
      
          float toteta=1000;
          float totphi=1000;
          float dr=1000;
          float EE=0;
          float feta=0;
          math::XYZPointF ElecTrkEcalPos(0,0,0);
          PFClusterRef clusterRef;
          math::XYZPoint meanShowerSaved;
          if(theOutParticle.getSuccess()!=0){
            ElecTrkEcalPos=math::XYZPointF(theOutParticle.vertex().x(),
                                           theOutParticle.vertex().y(),
                                           theOutParticle.vertex().z());
            bool isBelowPS=(fabs(theOutParticle.vertex().eta())>1.65) ? true :false;    
        
            unsigned clusCounter=0;

            for(vector<PFCluster>::const_iterator aClus = basClus.begin();
                aClus != basClus.end(); aClus++,++clusCounter) {
          
              double ecalShowerDepth
                = PFCluster::getDepthCorrection(aClus->energy(),
                                                isBelowPS,
                                                false);
        
              math::XYZPoint meanShower=math::XYZPoint(theOutParticle.vertex())+
                math::XYZTLorentzVector(theOutParticle.momentum()).Vect().Unit()*ecalShowerDepth;       
          
              float etarec=meanShower.eta();
              float phirec=meanShower.phi();
              float tmp_ep=aClus->energy()/PTOB;
              float tmp_phi=fabs(aClus->position().phi()-phirec);
              if (tmp_phi>TMath::TwoPi()) tmp_phi-= TMath::TwoPi();
              float tmp_dr=sqrt(pow(tmp_phi,2)+
                                pow((aClus->position().eta()-etarec),2));
          
              if ((tmp_dr<dr)&&(tmp_ep>minEp_)&&(tmp_ep<maxEp_)){
                dr=tmp_dr;
                toteta=aClus->position().eta()-etarec;
                totphi=tmp_phi;
                EP=tmp_ep;
                EE=aClus->energy();
                feta= aClus->position().eta();
                clusterRef = PFClusterRef(theECPfClustCollection,clusCounter);
                meanShowerSaved = meanShower;
              }
            }
          }

          //Resolution maps
          double ecaletares 
            = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta,EE));
          double ecalphires 
            = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta,EE)); 
      
          //geomatrical compatibility
          float chieta=(toteta!=1000)? toteta/ecaletares : toteta;
          float chiphi=(totphi!=1000)? totphi/ecalphires : totphi;
          float chichi= sqrt(chieta*chieta + chiphi*chiphi);
      
      
          //Matching criteria
          float chi2cut=thr[ibin+0];
          float ep_cutmin=thr[ibin+1];
          bool GoodMatching= ((chichi<chi2cut) &&(EP>ep_cutmin) && (nhitpi>10));
          bool EcalMatching=GoodMatching;
      
          if (Tk[i].pt()>maxPt_) GoodMatching=true;
          if (Tk[i].pt()<minPt_) GoodMatching=false;

          //ENDCAP
          //USE OF PRESHOWER 
          bool GoodPSMatching=false;
          if ((fabs(Tk[i].eta())>1.68)&&(usePreshower_)){
            int iptbin =4*getBin(Tk[i].pt());
            ps2En=0;ps1En=0;
            ps2chi=100.; ps1chi=100.;
            PSforTMVA(mom,pos);
            float p1e=thrPS[iptbin];
            float p2e=thrPS[iptbin+1];
            float p1c=thrPS[iptbin+2];
            float p2c=thrPS[iptbin+3];
            GoodPSMatching= 
              ((ps2En>p2e)
               &&(ps1En>p1e)
               &&(ps1chi<p1c)
               &&(ps2chi<p2c));
            GoodMatching = (GoodMatching && GoodPSMatching);
          }
  
          math::XYZPoint myPoint(ElecTrkEcalPos.X(),ElecTrkEcalPos.Y(),ElecTrkEcalPos.Z());
          myPreId.setECALMatchingProperties(clusterRef,myPoint,meanShowerSaved,toteta,totphi,chieta,
                                            chiphi,chichi,EP);
          myPreId.setECALMatching(EcalMatching);
          myPreId.setESMatching(GoodPSMatching);

          if(applyIsolation_){
            if(IsIsolated(float(Tk[i].charge()),Tk[i].p(),
                          ElecTrkEcalPos,*theECPfClustCollection,*theHCPfClustCollection)) 
              GoodMatching=true;
          }
          bool GoodRange= ((fabs(Tk[i].eta())<maxEta_) && 
                           (Tk[i].pt()>minPt_));
          //KF FILTERING FOR UNMATCHED EVENTS
          int hit1max=int(thr[ibin+2]);
          float chiredmin=thr[ibin+3];
          bool GoodKFFiltering =
            ((chikfred>chiredmin) || (nhitpi<hit1max));
      
          myPreId.setTrackFiltering(GoodKFFiltering);

          bool GoodTkId= false;
      
          if((!GoodMatching) &&(GoodKFFiltering) &&(GoodRange)){
            chi=chichi;
            chired=1000;
            chiRatio=1000;
            dpt=0;
            nhit=nhitpi;
      
            Trajectory::ConstRecHitContainer tmp;
            Trajectory::ConstRecHitContainer hits=Tj[i].recHits();
            for (int ih=hits.size()-1; ih>=0; ih--)  tmp.push_back(hits[ih]);
            vector<Trajectory> FitTjs=(fitter_.product())->fit(Seed,tmp,Tj[i].lastMeasurement().updatedState());
        
            if(FitTjs.size()>0){
              if(FitTjs[0].isValid()){
                vector<Trajectory> SmooTjs=(smoother_.product())->trajectories(FitTjs[0]);
                if(SmooTjs.size()>0){
                  if(SmooTjs[0].isValid()){
                
                    //Track refitted with electron hypothesis
                
                    float pt_out=SmooTjs[0].firstMeasurement().
                      updatedState().globalMomentum().perp();
                    float pt_in=SmooTjs[0].lastMeasurement().
                      updatedState().globalMomentum().perp();
                    dpt=(pt_in>0) ? fabs(pt_out-pt_in)/pt_in : 0.;
                    // the following is simply the number of degrees of freedom
                    chiRatio=SmooTjs[0].chiSquared()/Tj[i].chiSquared();
                    chired=chiRatio*chikfred;
                  }
                }
              }
            }
        
            //TMVA Analysis
            if(useTmva_){
        
              eta=Tk[i].eta();
              pt=Tk[i].pt();
              eP=EP;
          
              float Ytmva=reader->EvaluateMVA( method_ );
              
              float BDTcut=thr[ibin+4]; 
              if ( Ytmva>BDTcut) GoodTkId=true;
              myPreId.setMVA(GoodTkId,Ytmva);
              myPreId.setTrackProperties(chired,chiRatio,dpt);
            }else{ 
                          
              //
              float chiratiocut=thr[ibin+5]; 
              float gschicut=thr[ibin+6]; 
              float gsptmin=thr[ibin+7];

              GoodTkId=((dpt>gsptmin)&&(chired<gschicut)&&(chiRatio<chiratiocut));      
       
            }
          }
    
          GoodPreId= (GoodTkId || GoodMatching); 

          myPreId.setFinalDecision(GoodPreId);
      
          if(GoodPreId)
            LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
              "GeV/c, eta= "<<Tk[i].eta() <<
              ") preidentified for agreement between  track and ECAL cluster";
          if(GoodPreId &&(!GoodMatching))
            LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
              "GeV/c, eta= "<<Tk[i].eta() <<
              ") preidentified only for track properties";
        
        } // end of !disablePreId_
      
      if (GoodPreId){

        //NEW SEED with n hits  
        ElectronSeed NewSeed(Seed);
        NewSeed.setCtfTrack(trackRef);
        output_preid->push_back(NewSeed);
      }else{
        if (produceCkfseed_){
          output_nopre->push_back(Seed);
        }
      }
      if(producePreId_ && myPreId.pt()>PtThresholdSavePredId_)
        {
          // save the index of the PreId object as to be able to create a Ref later
          refMap_[trackRef] = output_preidinfo->size();
          output_preidinfo->push_back(myPreId);   
        }
    } //end loop on track collection
  } //end loop on the vector of track collections
  
  // no disablePreId_ switch, it is simpler to have an empty collection rather than no collection
  iEvent.put(output_preid,preidgsf_);
  if (produceCkfseed_)
    iEvent.put(output_nopre,preidckf_);
  if(producePreId_)
    {
      const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(output_preidinfo,preidname_);
      // now make the Value Map, but only if one input collection
      if(tracksContainers_.size()==1)
        {
          Handle<TrackCollection> tkRefCollection ;
          iEvent.getByLabel(tracksContainers_[0],tkRefCollection);
          fillPreIdRefValueMap(tkRefCollection,preIdRefProd,mapFiller);
          mapFiller.fill();
          iEvent.put(preIdMap_p,preidname_);
        }
    }
  
}

void GoodSeedProducer::PSforTMVA	(	math::XYZTLorentzVector	mom,
		math::XYZTLorentzVector	pos
	)		`[private]`

Definition at line 577 of file GoodSeedProducer.cc.

References B_, ExpressReco_HICollisions_FallBack::chi2, PFGeometry::innerRadius(), PFGeometry::outerRadius(), BaseParticlePropagator::propagateToPreshowerLayer1(), PFGeometry::PS1, ps1chi, ps1Clus, ps1En, PFGeometry::PS2, ps2chi, ps2Clus, ps2En, and mathSSE::sqrt().

Referenced by produce().

                                                                       {

  BaseParticlePropagator OutParticle(RawParticle(mom,pos)
                                     ,0.,0.,B_.z()) ;

  OutParticle.propagateToPreshowerLayer1(false);
  if (OutParticle.getSuccess()!=0){
    //   GlobalPoint v1=ps1TSOS.globalPosition();
    math::XYZPoint v1=math::XYZPoint(OutParticle.vertex());
    if ((v1.Rho() >=
         PFGeometry::innerRadius(PFGeometry::PS1)) &&
        (v1.Rho() <=
         PFGeometry::outerRadius(PFGeometry::PS1))) {
      float enPScl1=0;
      float chi1=100;
      vector<PFCluster>::const_iterator ips;
      for (ips=ps1Clus.begin(); ips!=ps1Clus.end();ips++){
        float ax=((*ips).position().x()-v1.x())/0.114;
        float ay=((*ips).position().y()-v1.y())/2.43;
        float pschi= sqrt(ax*ax+ay*ay);
        if (pschi<chi1){
          chi1=pschi;
          enPScl1=(*ips).energy();
        }
      }
      ps1En=enPScl1;
      ps1chi=chi1;


      OutParticle.propagateToPreshowerLayer2(false);
      if (OutParticle.getSuccess()!=0){
        math::XYZPoint v2=math::XYZPoint(OutParticle.vertex());
        if ((v2.Rho() >=
             PFGeometry::innerRadius(PFGeometry::PS2)) &&
            (v2.Rho() <=
             PFGeometry::outerRadius(PFGeometry::PS2))){
          float enPScl2=0;
          float chi2=100;
          for (ips=ps2Clus.begin(); ips!=ps2Clus.end();ips++){
            float ax=((*ips).position().x()-v2.x())/1.88;
            float ay=((*ips).position().y()-v2.y())/0.1449;
            float pschi= sqrt(ax*ax+ay*ay);
            if (pschi<chi2){
              chi2=pschi;
              enPScl2=(*ips).energy();
            }
          }

          ps2En=enPScl2;
          ps2chi=chi2;

        }
      }
      
    }
    
  }
}