CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    PFTracking
00004 // Class:      GoodSeedProducer
00005 // 
00006 // Original Author:  Michele Pioppi
00007 // March 2010. F. Beaudette. Produce PreId information
00008 
00009 
00010 #include "RecoParticleFlow/PFTracking/interface/GoodSeedProducer.h"
00011 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
00012 #include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"
00013 
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 #include "FWCore/ParameterSet/interface/FileInPath.h"
00016 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00017 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
00018 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00019 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
00020 #include "DataFormats/ParticleFlowReco/interface/PreId.h"
00021 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00022 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
00023 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"  
00024 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h" 
00025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00026 #include "FastSimulation/BaseParticlePropagator/interface/BaseParticlePropagator.h"
00027 #include "MagneticField/Engine/interface/MagneticField.h"
00028 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00029 #include <fstream>
00030 #include <string>
00031 #include "TMath.h"
00032 #include "Math/VectorUtil.h"
00033 
00034 using namespace edm;
00035 using namespace std;
00036 using namespace reco;
00037 PFResolutionMap* GoodSeedProducer::resMapEtaECAL_ = 0;                                        
00038 PFResolutionMap* GoodSeedProducer::resMapPhiECAL_ = 0;
00039 
00040 GoodSeedProducer::GoodSeedProducer(const ParameterSet& iConfig):
00041   pfTransformer_(0),
00042   conf_(iConfig)
00043 {
00044   LogInfo("GoodSeedProducer")<<"Electron PreIdentification started  ";
00045   
00046   //now do what ever initialization is needed
00047  
00048   tracksContainers_ = 
00049     iConfig.getParameter< vector < InputTag > >("TkColList");
00050   
00051   minPt_=iConfig.getParameter<double>("MinPt");
00052   maxPt_=iConfig.getParameter<double>("MaxPt");
00053   maxEta_=iConfig.getParameter<double>("MaxEta");
00054 
00055 
00056   //ISOLATION REQUEST AS DONE IN THE TAU GROUP
00057   applyIsolation_ =iConfig.getParameter<bool>("ApplyIsolation");
00058   HcalIsolWindow_                       =iConfig.getParameter<double>("HcalWindow");
00059   EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
00060   EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
00061   EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
00062   EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
00063    minEoverP_= iConfig.getParameter<double>("EOverPLead_minValue");
00064    maxHoverP_= iConfig.getParameter<double>("HOverPLead_maxValue");
00065  
00066   //
00067   pfCLusTagECLabel_=
00068     iConfig.getParameter<InputTag>("PFEcalClusterLabel");
00069 
00070   pfCLusTagHCLabel_=
00071    iConfig.getParameter<InputTag>("PFHcalClusterLabel");  
00072 
00073   pfCLusTagPSLabel_=
00074     iConfig.getParameter<InputTag>("PFPSClusterLabel");
00075   
00076   preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
00077   preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
00078   preidname_= iConfig.getParameter<string>("PreIdLabel");
00079   
00080   
00081   fitterName_ = iConfig.getParameter<string>("Fitter");
00082   smootherName_ = iConfig.getParameter<string>("Smoother");
00083   
00084   
00085   nHitsInSeed_=iConfig.getParameter<int>("NHitsInSeed");
00086 
00087   clusThreshold_=iConfig.getParameter<double>("ClusterThreshold");
00088   
00089   minEp_=iConfig.getParameter<double>("MinEOverP");
00090   maxEp_=iConfig.getParameter<double>("MaxEOverP");
00091 
00092   //collection to produce
00093   produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed",false);
00094 
00095   // to disable the electron part (for HI collisions for examples) 
00096   disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId",false);  
00097 
00098   producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId",true);  
00099   //  if no electron, cannot produce the preid
00100   if(disablePreId_) 
00101     producePreId_=false;
00102   PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId",1.);  
00103 
00104   LogDebug("GoodSeedProducer")<<"Seeds for GSF will be produced ";
00105 
00106   // no disablePreId_ switch here. The collection will be empty if it is true
00107   produces<ElectronSeedCollection>(preidgsf_);
00108 
00109   if(produceCkfseed_){
00110     LogDebug("GoodSeedProducer")<<"Seeds for CKF will be produced ";
00111     produces<TrajectorySeedCollection>(preidckf_);
00112   }
00113   
00114 
00115   if(producePreId_){
00116     LogDebug("GoodSeedProducer")<<"PreId debugging information will be produced ";
00117 
00118     produces<PreIdCollection>(preidname_);
00119     if(tracksContainers_.size()==1) // do not make a value map if more than one input track collection
00120       produces<edm::ValueMap<reco::PreIdRef> >(preidname_);
00121   } 
00122   
00123   useQuality_   = iConfig.getParameter<bool>("UseQuality");
00124   trackQuality_=TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
00125 
00126   useTmva_= iConfig.getUntrackedParameter<bool>("UseTMVA",false);
00127   
00128   usePreshower_ = iConfig.getParameter<bool>("UsePreShower");
00129 
00130 }
00131 
00132 
00133 GoodSeedProducer::~GoodSeedProducer()
00134 {
00135   
00136   // do anything here that needs to be done at desctruction time
00137   // (e.g. close files, deallocate resources etc.) 
00138 
00139   delete pfTransformer_;
00140 }
00141 
00142 
00143 //
00144 // member functions
00145 //
00146 
00147 // ------------ method called to produce the data  ------------
00148 void
00149 GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup)
00150 {
00151 
00152   LogDebug("GoodSeedProducer")<<"START event: "<<iEvent.id().event()
00153                               <<" in run "<<iEvent.id().run();
00154   //Create empty output collections
00155   auto_ptr<ElectronSeedCollection> output_preid(new ElectronSeedCollection);
00156   auto_ptr<TrajectorySeedCollection> output_nopre(new TrajectorySeedCollection);
00157   auto_ptr<PreIdCollection> output_preidinfo(new PreIdCollection);
00158   auto_ptr<edm::ValueMap<reco::PreIdRef> > preIdMap_p(new edm::ValueMap<reco::PreIdRef>);
00159   edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);
00160 
00161   //Tracking Tools
00162   if(!disablePreId_)
00163     {
00164       iSetup.get<TrajectoryFitter::Record>().get(fitterName_, fitter_);
00165       iSetup.get<TrajectoryFitter::Record>().get(smootherName_, smoother_);
00166     }
00167 
00168   // clear temporary maps
00169   refMap_.clear();
00170 
00171   //Handle input collections
00172   //ECAL clusters             
00173   Handle<PFClusterCollection> theECPfClustCollection;
00174   iEvent.getByLabel(pfCLusTagECLabel_,theECPfClustCollection);
00175   
00176   vector<PFCluster> basClus;
00177   vector<PFCluster>::const_iterator iklus;
00178   for (iklus=theECPfClustCollection.product()->begin();
00179        iklus!=theECPfClustCollection.product()->end();
00180        iklus++){
00181     if((*iklus).energy()>clusThreshold_) basClus.push_back(*iklus);
00182   }
00183 
00184   //HCAL clusters
00185   Handle<PFClusterCollection> theHCPfClustCollection;
00186   iEvent.getByLabel(pfCLusTagHCLabel_,theHCPfClustCollection);
00187   
00188   //PS clusters
00189   Handle<PFClusterCollection> thePSPfClustCollection;
00190   iEvent.getByLabel(pfCLusTagPSLabel_,thePSPfClustCollection);
00191   
00192   ps1Clus.clear();
00193   ps2Clus.clear();
00194   
00195   for (iklus=thePSPfClustCollection.product()->begin();
00196        iklus!=thePSPfClustCollection.product()->end();
00197        iklus++){
00198     //layer==-11 first layer of PS
00199     //layer==-12 secon layer of PS
00200     if ((*iklus).layer()==-11) ps1Clus.push_back(*iklus);
00201     if ((*iklus).layer()==-12) ps2Clus.push_back(*iklus);
00202   }
00203   
00204   //Vector of track collections
00205   for (unsigned int istr=0; istr<tracksContainers_.size();istr++){
00206     
00207     //Track collection
00208     Handle<TrackCollection> tkRefCollection;
00209     iEvent.getByLabel(tracksContainers_[istr], tkRefCollection);
00210     const TrackCollection&  Tk=*(tkRefCollection.product());
00211     
00212     //Trajectory collection
00213     Handle<vector<Trajectory> > tjCollection;
00214     iEvent.getByLabel(tracksContainers_[istr], tjCollection);
00215     vector<Trajectory> Tj=*(tjCollection.product());
00216     
00217     
00218     LogDebug("GoodSeedProducer")<<"Number of tracks in collection "
00219                                 <<tracksContainers_[istr] <<" to be analyzed "
00220                                 <<Tj.size();
00221     
00222 
00223     //loop over the track collection
00224     for(unsigned int i=0;i<Tk.size();i++){              
00225       if (useQuality_ &&
00226           (!(Tk[i].quality(trackQuality_)))) continue;
00227       
00228       reco::PreId myPreId;
00229       bool GoodPreId=false;
00230 
00231       TrackRef trackRef(tkRefCollection, i);
00232       // TrajectorySeed Seed=Tj[i].seed();
00233       TrajectorySeed Seed=(*trackRef->seedRef());
00234       if(!disablePreId_)
00235         {
00236           int ipteta=getBin(Tk[i].eta(),Tk[i].pt());
00237           int ibin=ipteta*8;
00238           
00239           float PTOB=Tj[i].lastMeasurement().updatedState().globalMomentum().mag();
00240           float chikfred=Tk[i].normalizedChi2();
00241           int nhitpi=Tj[i].foundHits();
00242           float EP=0;
00243       
00244           // set track info
00245           myPreId.setTrack(trackRef);
00246           //CLUSTERS - TRACK matching
00247       
00248           float pfmass=  0.0005;
00249           float pfoutenergy=sqrt((pfmass*pfmass)+Tk[i].outerMomentum().Mag2());
00250           XYZTLorentzVector mom =XYZTLorentzVector(Tk[i].outerMomentum().x(),
00251                                                    Tk[i].outerMomentum().y(),
00252                                                    Tk[i].outerMomentum().z(),
00253                                                    pfoutenergy);
00254           XYZTLorentzVector pos =   XYZTLorentzVector(Tk[i].outerPosition().x(),
00255                                                       Tk[i].outerPosition().y(),
00256                                                       Tk[i].outerPosition().z(),
00257                                                       0.);
00258 
00259           BaseParticlePropagator theOutParticle = 
00260             BaseParticlePropagator( RawParticle(mom,pos),
00261                                     0,0,B_.z());
00262           theOutParticle.setCharge(Tk[i].charge());
00263       
00264           theOutParticle.propagateToEcalEntrance(false);
00265       
00266       
00267           float toteta=1000;
00268           float totphi=1000;
00269           float dr=1000;
00270           float EE=0;
00271           float feta=0;
00272           math::XYZPointF ElecTrkEcalPos(0,0,0);
00273           PFClusterRef clusterRef;
00274           math::XYZPoint meanShowerSaved;
00275           if(theOutParticle.getSuccess()!=0){
00276             ElecTrkEcalPos=math::XYZPointF(theOutParticle.vertex().x(),
00277                                            theOutParticle.vertex().y(),
00278                                            theOutParticle.vertex().z());
00279             bool isBelowPS=(fabs(theOutParticle.vertex().eta())>1.65) ? true :false;    
00280         
00281             unsigned clusCounter=0;
00282 
00283             for(vector<PFCluster>::const_iterator aClus = basClus.begin();
00284                 aClus != basClus.end(); aClus++,++clusCounter) {
00285           
00286               double ecalShowerDepth
00287                 = PFCluster::getDepthCorrection(aClus->energy(),
00288                                                 isBelowPS,
00289                                                 false);
00290         
00291               math::XYZPoint meanShower=math::XYZPoint(theOutParticle.vertex())+
00292                 math::XYZTLorentzVector(theOutParticle.momentum()).Vect().Unit()*ecalShowerDepth;       
00293           
00294               float etarec=meanShower.eta();
00295               float phirec=meanShower.phi();
00296               float tmp_ep=aClus->energy()/PTOB;
00297               float tmp_phi=fabs(aClus->position().phi()-phirec);
00298               if (tmp_phi>TMath::TwoPi()) tmp_phi-= TMath::TwoPi();
00299               float tmp_dr=sqrt(pow(tmp_phi,2)+
00300                                 pow((aClus->position().eta()-etarec),2));
00301           
00302               if ((tmp_dr<dr)&&(tmp_ep>minEp_)&&(tmp_ep<maxEp_)){
00303                 dr=tmp_dr;
00304                 toteta=aClus->position().eta()-etarec;
00305                 totphi=tmp_phi;
00306                 EP=tmp_ep;
00307                 EE=aClus->energy();
00308                 feta= aClus->position().eta();
00309                 clusterRef = PFClusterRef(theECPfClustCollection,clusCounter);
00310                 meanShowerSaved = meanShower;
00311               }
00312             }
00313           }
00314 
00315           //Resolution maps
00316           double ecaletares 
00317             = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta,EE));
00318           double ecalphires 
00319             = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta,EE)); 
00320       
00321           //geomatrical compatibility
00322           float chieta=(toteta!=1000)? toteta/ecaletares : toteta;
00323           float chiphi=(totphi!=1000)? totphi/ecalphires : totphi;
00324           float chichi= sqrt(chieta*chieta + chiphi*chiphi);
00325       
00326       
00327           //Matching criteria
00328           float chi2cut=thr[ibin+0];
00329           float ep_cutmin=thr[ibin+1];
00330           bool GoodMatching= ((chichi<chi2cut) &&(EP>ep_cutmin) && (nhitpi>10));
00331           bool EcalMatching=GoodMatching;
00332       
00333           if (Tk[i].pt()>maxPt_) GoodMatching=true;
00334           if (Tk[i].pt()<minPt_) GoodMatching=false;
00335 
00336           //ENDCAP
00337           //USE OF PRESHOWER 
00338           bool GoodPSMatching=false;
00339           if ((fabs(Tk[i].eta())>1.68)&&(usePreshower_)){
00340             int iptbin =4*getBin(Tk[i].pt());
00341             ps2En=0;ps1En=0;
00342             ps2chi=100.; ps1chi=100.;
00343             PSforTMVA(mom,pos);
00344             float p1e=thrPS[iptbin];
00345             float p2e=thrPS[iptbin+1];
00346             float p1c=thrPS[iptbin+2];
00347             float p2c=thrPS[iptbin+3];
00348             GoodPSMatching= 
00349               ((ps2En>p2e)
00350                &&(ps1En>p1e)
00351                &&(ps1chi<p1c)
00352                &&(ps2chi<p2c));
00353             GoodMatching = (GoodMatching && GoodPSMatching);
00354           }
00355   
00356           math::XYZPoint myPoint(ElecTrkEcalPos.X(),ElecTrkEcalPos.Y(),ElecTrkEcalPos.Z());
00357           myPreId.setECALMatchingProperties(clusterRef,myPoint,meanShowerSaved,toteta,totphi,chieta,
00358                                             chiphi,chichi,EP);
00359           myPreId.setECALMatching(EcalMatching);
00360           myPreId.setESMatching(GoodPSMatching);
00361 
00362           if(applyIsolation_){
00363             if(IsIsolated(float(Tk[i].charge()),Tk[i].p(),
00364                           ElecTrkEcalPos,*theECPfClustCollection,*theHCPfClustCollection)) 
00365               GoodMatching=true;
00366           }
00367           bool GoodRange= ((fabs(Tk[i].eta())<maxEta_) && 
00368                            (Tk[i].pt()>minPt_));
00369           //KF FILTERING FOR UNMATCHED EVENTS
00370           int hit1max=int(thr[ibin+2]);
00371           float chiredmin=thr[ibin+3];
00372           bool GoodKFFiltering =
00373             ((chikfred>chiredmin) || (nhitpi<hit1max));
00374       
00375           myPreId.setTrackFiltering(GoodKFFiltering);
00376 
00377           bool GoodTkId= false;
00378       
00379           if((!GoodMatching) &&(GoodKFFiltering) &&(GoodRange)){
00380             chi=chichi;
00381             chired=1000;
00382             chiRatio=1000;
00383             dpt=0;
00384             nhit=nhitpi;
00385       
00386             Trajectory::ConstRecHitContainer tmp;
00387             Trajectory::ConstRecHitContainer hits=Tj[i].recHits();
00388             for (int ih=hits.size()-1; ih>=0; ih--)  tmp.push_back(hits[ih]);
00389             vector<Trajectory> FitTjs=(fitter_.product())->fit(Seed,tmp,Tj[i].lastMeasurement().updatedState());
00390         
00391             if(FitTjs.size()>0){
00392               if(FitTjs[0].isValid()){
00393                 vector<Trajectory> SmooTjs=(smoother_.product())->trajectories(FitTjs[0]);
00394                 if(SmooTjs.size()>0){
00395                   if(SmooTjs[0].isValid()){
00396                 
00397                     //Track refitted with electron hypothesis
00398                 
00399                     float pt_out=SmooTjs[0].firstMeasurement().
00400                       updatedState().globalMomentum().perp();
00401                     float pt_in=SmooTjs[0].lastMeasurement().
00402                       updatedState().globalMomentum().perp();
00403                     dpt=(pt_in>0) ? fabs(pt_out-pt_in)/pt_in : 0.;
00404                     // the following is simply the number of degrees of freedom
00405                     chiRatio=SmooTjs[0].chiSquared()/Tj[i].chiSquared();
00406                     chired=chiRatio*chikfred;
00407                   }
00408                 }
00409               }
00410             }
00411         
00412             //TMVA Analysis
00413             if(useTmva_){
00414         
00415               eta=Tk[i].eta();
00416               pt=Tk[i].pt();
00417               eP=EP;
00418           
00419               float Ytmva=reader->EvaluateMVA( method_ );
00420               
00421               float BDTcut=thr[ibin+4]; 
00422               if ( Ytmva>BDTcut) GoodTkId=true;
00423               myPreId.setMVA(GoodTkId,Ytmva);
00424               myPreId.setTrackProperties(chired,chiRatio,dpt);
00425             }else{ 
00426                           
00427               //
00428               float chiratiocut=thr[ibin+5]; 
00429               float gschicut=thr[ibin+6]; 
00430               float gsptmin=thr[ibin+7];
00431 
00432               GoodTkId=((dpt>gsptmin)&&(chired<gschicut)&&(chiRatio<chiratiocut));      
00433        
00434             }
00435           }
00436     
00437           GoodPreId= (GoodTkId || GoodMatching); 
00438 
00439           myPreId.setFinalDecision(GoodPreId);
00440       
00441           if(GoodPreId)
00442             LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
00443               "GeV/c, eta= "<<Tk[i].eta() <<
00444               ") preidentified for agreement between  track and ECAL cluster";
00445           if(GoodPreId &&(!GoodMatching))
00446             LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
00447               "GeV/c, eta= "<<Tk[i].eta() <<
00448               ") preidentified only for track properties";
00449         
00450         } // end of !disablePreId_
00451       
00452       if (GoodPreId){
00453 
00454         //NEW SEED with n hits  
00455         ElectronSeed NewSeed(Seed);
00456         NewSeed.setCtfTrack(trackRef);
00457         output_preid->push_back(NewSeed);
00458       }else{
00459         if (produceCkfseed_){
00460           output_nopre->push_back(Seed);
00461         }
00462       }
00463       if(producePreId_ && myPreId.pt()>PtThresholdSavePredId_)
00464         {
00465           // save the index of the PreId object as to be able to create a Ref later
00466           refMap_[trackRef] = output_preidinfo->size();
00467           output_preidinfo->push_back(myPreId);   
00468         }
00469     } //end loop on track collection
00470   } //end loop on the vector of track collections
00471   
00472   // no disablePreId_ switch, it is simpler to have an empty collection rather than no collection
00473   iEvent.put(output_preid,preidgsf_);
00474   if (produceCkfseed_)
00475     iEvent.put(output_nopre,preidckf_);
00476   if(producePreId_)
00477     {
00478       const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(output_preidinfo,preidname_);
00479       // now make the Value Map, but only if one input collection
00480       if(tracksContainers_.size()==1)
00481         {
00482           Handle<TrackCollection> tkRefCollection ;
00483           iEvent.getByLabel(tracksContainers_[0],tkRefCollection);
00484           fillPreIdRefValueMap(tkRefCollection,preIdRefProd,mapFiller);
00485           mapFiller.fill();
00486           iEvent.put(preIdMap_p,preidname_);
00487         }
00488     }
00489 
00490    // clear temporary maps
00491   refMap_.clear();
00492 
00493 }
00494 // ------------ method called once each job just before starting event loop  ------------
00495 void 
00496 GoodSeedProducer::beginRun(edm::Run & run,
00497                            const EventSetup& es)
00498 {
00499   //Magnetic Field
00500   ESHandle<MagneticField> magneticField;
00501   es.get<IdealMagneticFieldRecord>().get(magneticField);
00502   B_=magneticField->inTesla(GlobalPoint(0,0,0));
00503   
00504   pfTransformer_= new PFTrackTransformer(B_);
00505   pfTransformer_->OnlyProp();
00506 
00507 
00508   
00509   //Resolution maps
00510   FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
00511   FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
00512   resMapEtaECAL_ = new PFResolutionMap("ECAL_eta",ecalEtaMap.fullPath().c_str());
00513   resMapPhiECAL_ = new PFResolutionMap("ECAL_phi",ecalPhiMap.fullPath().c_str());
00514 
00515   if(useTmva_){
00516     reader = new TMVA::Reader("!Color:Silent");
00517     method_ = conf_.getParameter<string>("TMVAMethod");
00518     
00519     reader->AddVariable("eP",&eP);
00520     reader->AddVariable("chi",&chi);
00521     reader->AddVariable("chired",&chired);
00522     reader->AddVariable("chiRatio",&chiRatio);
00523     reader->AddVariable("dpt",&dpt);
00524     reader->AddVariable("nhit",&nhit);
00525     reader->AddVariable("eta",&eta);
00526     reader->AddVariable("pt",&pt);
00527     FileInPath Weigths(conf_.getParameter<string>("Weights"));
00528     reader->BookMVA( method_, Weigths.fullPath().c_str()  );    
00529     }
00530 
00531     
00532     //read threshold
00533     FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
00534     ifstream ifs(parFile.fullPath().c_str());
00535     for (int iy=0;iy<72;iy++) ifs >> thr[iy];
00536     
00537     //read PS threshold
00538     FileInPath parPSFile(conf_.getParameter<string>("PSThresholdFile"));
00539     ifstream ifsPS(parPSFile.fullPath().c_str());
00540     for (int iy=0;iy<12;iy++) ifsPS >> thrPS[iy];
00541 
00542 }
00543 
00544 void 
00545 GoodSeedProducer::endRun() {
00546   delete pfTransformer_;
00547   delete resMapEtaECAL_;
00548   delete resMapPhiECAL_;
00549   if(useTmva_) delete reader;
00550 }
00551 
00552 int 
00553 GoodSeedProducer::getBin(float pt){
00554 int ip=0;
00555   if (pt<6) ip=0;
00556   else {  if (pt<12) ip=1;
00557         else ip=2;
00558   }
00559 return ip;
00560 }
00561 
00562 int 
00563 GoodSeedProducer::getBin(float eta, float pt){
00564   int ie=0;
00565   int ip=0;
00566   if (fabs(eta)<1.2) ie=0;
00567   else{ if (fabs(eta)<1.68) ie=1;
00568     else ie=2;
00569   }
00570   if (pt<6) ip=0;
00571   else {  if (pt<12) ip=1;     
00572         else ip=2;
00573   }
00574   int iep= ie*3+ip;
00575   LogDebug("GoodSeedProducer")<<"Track pt ="<<pt<<" eta="<<eta<<" bin="<<iep;
00576   return iep;
00577 }
00578 
00579 void 
00580 GoodSeedProducer::PSforTMVA(XYZTLorentzVector mom,XYZTLorentzVector pos ){
00581 
00582   BaseParticlePropagator OutParticle(RawParticle(mom,pos)
00583                                      ,0.,0.,B_.z()) ;
00584 
00585   OutParticle.propagateToPreshowerLayer1(false);
00586   if (OutParticle.getSuccess()!=0){
00587     //   GlobalPoint v1=ps1TSOS.globalPosition();
00588     math::XYZPoint v1=math::XYZPoint(OutParticle.vertex());
00589     if ((v1.Rho() >=
00590          PFGeometry::innerRadius(PFGeometry::PS1)) &&
00591         (v1.Rho() <=
00592          PFGeometry::outerRadius(PFGeometry::PS1))) {
00593       float enPScl1=0;
00594       float chi1=100;
00595       vector<PFCluster>::const_iterator ips;
00596       for (ips=ps1Clus.begin(); ips!=ps1Clus.end();ips++){
00597         float ax=((*ips).position().x()-v1.x())/0.114;
00598         float ay=((*ips).position().y()-v1.y())/2.43;
00599         float pschi= sqrt(ax*ax+ay*ay);
00600         if (pschi<chi1){
00601           chi1=pschi;
00602           enPScl1=(*ips).energy();
00603         }
00604       }
00605       ps1En=enPScl1;
00606       ps1chi=chi1;
00607 
00608 
00609       OutParticle.propagateToPreshowerLayer2(false);
00610       if (OutParticle.getSuccess()!=0){
00611         math::XYZPoint v2=math::XYZPoint(OutParticle.vertex());
00612         if ((v2.Rho() >=
00613              PFGeometry::innerRadius(PFGeometry::PS2)) &&
00614             (v2.Rho() <=
00615              PFGeometry::outerRadius(PFGeometry::PS2))){
00616           float enPScl2=0;
00617           float chi2=100;
00618           for (ips=ps2Clus.begin(); ips!=ps2Clus.end();ips++){
00619             float ax=((*ips).position().x()-v2.x())/1.88;
00620             float ay=((*ips).position().y()-v2.y())/0.1449;
00621             float pschi= sqrt(ax*ax+ay*ay);
00622             if (pschi<chi2){
00623               chi2=pschi;
00624               enPScl2=(*ips).energy();
00625             }
00626           }
00627 
00628           ps2En=enPScl2;
00629           ps2chi=chi2;
00630 
00631         }
00632       }
00633       
00634     }
00635     
00636   }
00637 }
00638 
00639 bool 
00640 GoodSeedProducer::IsIsolated(float charge, float P,
00641                              math::XYZPointF myElecTrkEcalPos,
00642                              const PFClusterCollection &ecalColl,
00643                              const PFClusterCollection &hcalColl){
00644 
00645 
00646   double myHCALenergy3x3=0.;
00647   double myStripClusterE=0.;
00648  
00649 
00650   //  reco::TrackRef myElecTrk;
00651   
00652   if (fabs(myElecTrkEcalPos.z())<1. && myElecTrkEcalPos.x()<1. && myElecTrkEcalPos.y()<1. ) return false; 
00653 
00654   
00655   
00656   PFClusterCollection::const_iterator hc=hcalColl.begin();
00657   PFClusterCollection::const_iterator hcend=hcalColl.end();
00658   for (;hc!=hcend;++hc){
00659     math::XYZPoint clusPos = hc->position();
00660     double en = hc->energy();
00661     double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos,clusPos);
00662     if (deltaR<HcalIsolWindow_) {
00663       myHCALenergy3x3 += en;
00664       
00665     }
00666   }
00667 
00668 
00669 
00670   PFClusterCollection::const_iterator ec=ecalColl.begin();
00671   PFClusterCollection::const_iterator ecend=ecalColl.end();
00672   for (;ec!=ecend;++ec){
00673     math::XYZPoint clusPos = ec->position();
00674     double en = ec->energy();
00675 
00676 
00677     double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos,clusPos);
00678     double deltaEta = abs(myElecTrkEcalPos.eta()-clusPos.eta());
00679     double deltaPhiOverQ = deltaPhi/charge;
00680     if (en >= EcalStripSumE_minClusEnergy_ && deltaEta<EcalStripSumE_deltaEta_ && deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ && deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) { 
00681       myStripClusterE += en;
00682     }
00683   }       
00684   
00685   double EoP=myStripClusterE/P;
00686   double HoP=myHCALenergy3x3/P;
00687 
00688 
00689   return ((EoP>minEoverP_)&&(EoP<2.5) && (HoP<maxHoverP_))?true:false;
00690 }
00691 
00692 void GoodSeedProducer::fillPreIdRefValueMap( Handle<TrackCollection> tracks,
00693                                              const edm::OrphanHandle<reco::PreIdCollection>& preidhandle,
00694                                              edm::ValueMap<reco::PreIdRef>::Filler & filler)
00695 {
00696   std::vector<reco::PreIdRef> values;
00697 
00698   unsigned ntracks=tracks->size();
00699   for(unsigned itrack=0;itrack<ntracks;++itrack)
00700    {
00701      reco::TrackRef theTrackRef(tracks,itrack);
00702      std::map<reco::TrackRef,unsigned>::const_iterator itcheck=refMap_.find(theTrackRef);
00703      if(itcheck==refMap_.end()) 
00704        {
00705          // the track has been early discarded
00706          values.push_back(reco::PreIdRef());
00707        }
00708      else
00709        {
00710          edm::Ref<reco::PreIdCollection> preIdRef(preidhandle,itcheck->second);
00711          values.push_back(preIdRef);
00712          //      std::cout << " Checking Refs " << (theTrackRef==preIdRef->trackRef()) << std::endl;
00713        }
00714    }
00715   filler.insert(tracks,values.begin(),values.end());
00716 }