CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/RecoParticleFlow/PFProducer/plugins/PFElectronTranslator.cc

Go to the documentation of this file.
00001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00002 #include "FWCore/Framework/interface/ESHandle.h"
00003 #include "FWCore/Framework/interface/Event.h"
00004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00005 #include "RecoParticleFlow/PFProducer/plugins/PFElectronTranslator.h"
00006 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
00007 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
00008 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateElectronExtra.h"
00010 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00011 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
00012 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00013 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00014 #include "DataFormats/EgammaCandidates/interface/GsfElectronCore.h"
00015 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00016 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
00017 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
00018 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
00019 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
00020 
00021 
00022 PFElectronTranslator::PFElectronTranslator(const edm::ParameterSet & iConfig) {
00023   inputTagPFCandidates_ 
00024     = iConfig.getParameter<edm::InputTag>("PFCandidate");
00025   inputTagPFCandidateElectrons_ 
00026     = iConfig.getParameter<edm::InputTag>("PFCandidateElectron");
00027   inputTagGSFTracks_
00028     = iConfig.getParameter<edm::InputTag>("GSFTracks");
00029 
00030   edm::ParameterSet isoVals  = iConfig.getParameter<edm::ParameterSet> ("isolationValues");
00031   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfChargedHadrons"));
00032   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfPhotons"));
00033   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfNeutralHadrons"));
00034 
00035   PFBasicClusterCollection_ = iConfig.getParameter<std::string>("PFBasicClusters");
00036   PFPreshowerClusterCollection_ = iConfig.getParameter<std::string>("PFPreshowerClusters");
00037   PFSuperClusterCollection_ = iConfig.getParameter<std::string>("PFSuperClusters");
00038   GsfElectronCoreCollection_ = iConfig.getParameter<std::string>("PFGsfElectronCore");
00039   GsfElectronCollection_ = iConfig.getParameter<std::string>("PFGsfElectron");
00040   
00041   PFMVAValueMap_ = iConfig.getParameter<std::string>("ElectronMVA");
00042   PFSCValueMap_ = iConfig.getParameter<std::string>("ElectronSC");
00043   MVACut_ = (iConfig.getParameter<edm::ParameterSet>("MVACutBlock")).getParameter<double>("MVACut");
00044   checkStatusFlag_ = iConfig.getParameter<bool>("CheckStatusFlag");
00045   
00046   if (iConfig.exists("emptyIsOk")) emptyIsOk_ = iConfig.getParameter<bool>("emptyIsOk");
00047   else emptyIsOk_=false;
00048 
00049   produces<reco::BasicClusterCollection>(PFBasicClusterCollection_); 
00050   produces<reco::PreshowerClusterCollection>(PFPreshowerClusterCollection_); 
00051   produces<reco::SuperClusterCollection>(PFSuperClusterCollection_); 
00052   produces<reco::GsfElectronCoreCollection>(GsfElectronCoreCollection_);
00053   produces<reco::GsfElectronCollection>(GsfElectronCollection_);
00054   produces<edm::ValueMap<float> >(PFMVAValueMap_);
00055   produces<edm::ValueMap<reco::SuperClusterRef> >(PFSCValueMap_);
00056 }
00057 
00058 PFElectronTranslator::~PFElectronTranslator() {}
00059 
00060 void PFElectronTranslator::beginRun(edm::Run& run,const edm::EventSetup & es) {}
00061 
00062 void PFElectronTranslator::produce(edm::Event& iEvent,  
00063                                     const edm::EventSetup& iSetup) { 
00064   
00065   std::auto_ptr<reco::GsfElectronCoreCollection>
00066     gsfElectronCores_p(new reco::GsfElectronCoreCollection);
00067 
00068   std::auto_ptr<reco::GsfElectronCollection>
00069     gsfElectrons_p(new reco::GsfElectronCollection);
00070 
00071   std::auto_ptr<reco::SuperClusterCollection> 
00072     superClusters_p(new reco::SuperClusterCollection);
00073 
00074   std::auto_ptr<reco::BasicClusterCollection> 
00075     basicClusters_p(new reco::BasicClusterCollection);
00076 
00077   std::auto_ptr<reco::PreshowerClusterCollection>
00078     psClusters_p(new reco::PreshowerClusterCollection);
00079   
00080   std::auto_ptr<edm::ValueMap<float> > mvaMap_p(new edm::ValueMap<float>());
00081   edm::ValueMap<float>::Filler mvaFiller(*mvaMap_p);
00082 
00083   std::auto_ptr<edm::ValueMap<reco::SuperClusterRef> > 
00084     scMap_p(new edm::ValueMap<reco::SuperClusterRef>());
00085   edm::ValueMap<reco::SuperClusterRef>::Filler scRefFiller(*scMap_p);
00086   
00087 
00088   edm::Handle<reco::PFCandidateCollection> pfCandidates;
00089   bool status=fetchCandidateCollection(pfCandidates, 
00090                                        inputTagPFCandidates_, 
00091                                        iEvent );
00092 
00093   IsolationValueMaps isolationValues(inputTagIsoVals_.size());
00094   for (size_t j = 0; j<inputTagIsoVals_.size(); ++j) {
00095     iEvent.getByLabel(inputTagIsoVals_[j], isolationValues[j]);
00096   }
00097 
00098 
00099   // clear the vectors
00100   GsfTrackRef_.clear();
00101   CandidatePtr_.clear();
00102   ambiguousGsfTracks_.clear();
00103   kfTrackRef_.clear();
00104   basicClusters_.clear();
00105   pfClusters_.clear();
00106   preshowerClusters_.clear();
00107   superClusters_.clear();
00108   basicClusterPtr_.clear();
00109   preshowerClusterPtr_.clear();
00110   gsfPFCandidateIndex_.clear();
00111   gsfElectronCoreRefs_.clear();
00112   scMap_.clear();
00113 
00114  
00115   // loop on the candidates 
00116   //CC@@
00117   // we need first to create AND put the SuperCluster, 
00118   // basic clusters and presh clusters collection 
00119   // in order to get a working Handle
00120   unsigned ncand=(status)?pfCandidates->size():0;
00121   unsigned iGSF=0;
00122   for( unsigned i=0; i<ncand; ++i ) {
00123 
00124     const reco::PFCandidate& cand = (*pfCandidates)[i];    
00125     if(cand.particleId()!=reco::PFCandidate::e) continue; 
00126     if(cand.gsfTrackRef().isNull()) continue;
00127     // Note that -1 will still cut some total garbage candidates 
00128     // Fill the MVA map
00129     if(cand.mva_e_pi()<MVACut_) continue;
00130     
00131     // Check the status flag
00132     if(checkStatusFlag_ && !cand.electronExtraRef()->electronStatus(reco::PFCandidateElectronExtra::Selected)) {
00133       continue;
00134     }
00135 
00136     GsfTrackRef_.push_back(cand.gsfTrackRef());
00137     kfTrackRef_.push_back(cand.trackRef());
00138     gsfPFCandidateIndex_.push_back(i);
00139 
00140     reco::PFCandidatePtr ptrToPFElectron(pfCandidates,i);
00141     //CandidatePtr_.push_back(ptrToPFElectron->sourceCandidatePtr(0));
00142     CandidatePtr_.push_back(ptrToPFElectron);    
00143 
00144     basicClusters_.push_back(reco::BasicClusterCollection());
00145     pfClusters_.push_back(std::vector<const reco::PFCluster *>());
00146     preshowerClusters_.push_back(reco::PreshowerClusterCollection());
00147     ambiguousGsfTracks_.push_back(std::vector<reco::GsfTrackRef>());
00148 
00149     for(unsigned iele=0; iele<cand.elementsInBlocks().size(); ++iele) {
00150       // first get the block 
00151       reco::PFBlockRef blockRef = cand.elementsInBlocks()[iele].first;
00152       //
00153       unsigned elementIndex = cand.elementsInBlocks()[iele].second;
00154       // check it actually exists 
00155       if(blockRef.isNull()) continue;
00156       
00157       // then get the elements of the block
00158       const edm::OwnVector< reco::PFBlockElement >&  elements = (*blockRef).elements();
00159       
00160       const reco::PFBlockElement & pfbe (elements[elementIndex]); 
00161       // The first ECAL element should be the cluster associated to the GSF; defined as the seed
00162       if(pfbe.type()==reco::PFBlockElement::ECAL)
00163         {         
00164           //      const reco::PFCandidate * coCandidate = &cand;
00165           // the Brem photons are saved as daughter PFCandidate; this 
00166           // is convenient to access the corrected energy
00167           //      std::cout << " Found candidate "  << correspondingDaughterCandidate(coCandidate,pfbe) << " " << coCandidate << std::endl;
00168           createBasicCluster(pfbe,basicClusters_[iGSF],pfClusters_[iGSF],correspondingDaughterCandidate(cand,pfbe));
00169         }
00170       if(pfbe.type()==reco::PFBlockElement::PS1)
00171         {
00172           createPreshowerCluster(pfbe,preshowerClusters_[iGSF],1);
00173         }
00174       if(pfbe.type()==reco::PFBlockElement::PS2)
00175         {
00176           createPreshowerCluster(pfbe,preshowerClusters_[iGSF],2);
00177         }      
00178       if(pfbe.type()==reco::PFBlockElement::GSF)
00179         {
00180           getAmbiguousGsfTracks(pfbe,ambiguousGsfTracks_[iGSF]);
00181         }
00182 
00183     }   // loop on the elements
00184 
00185     // save the basic clusters
00186     basicClusters_p->insert(basicClusters_p->end(),basicClusters_[iGSF].begin(), basicClusters_[iGSF].end());
00187     // save the preshower clusters
00188     psClusters_p->insert(psClusters_p->end(),preshowerClusters_[iGSF].begin(),preshowerClusters_[iGSF].end());
00189 
00190     ++iGSF;
00191   } // loop on PFCandidates
00192 
00193   
00194    //Save the basic clusters and get an handle as to be able to create valid Refs (thanks to Claude)
00195   //  std::cout << " Number of basic clusters " << basicClusters_p->size() << std::endl;
00196   const edm::OrphanHandle<reco::BasicClusterCollection> bcRefProd = 
00197     iEvent.put(basicClusters_p,PFBasicClusterCollection_);
00198 
00199   //preshower clusters
00200   const edm::OrphanHandle<reco::PreshowerClusterCollection> psRefProd = 
00201     iEvent.put(psClusters_p,PFPreshowerClusterCollection_);
00202 
00203   // now that the Basic clusters are in the event, the Ref can be created
00204   createBasicClusterPtrs(bcRefProd);
00205   // now that the preshower clusters are in the event, the Ref can be created
00206   createPreshowerClusterPtrs(psRefProd);
00207   
00208   // and now the Super cluster can be created with valid references  
00209   if(status) createSuperClusters(*pfCandidates,*superClusters_p);
00210   
00211   // Let's put the super clusters in the event
00212   const edm::OrphanHandle<reco::SuperClusterCollection> scRefProd = iEvent.put(superClusters_p,PFSuperClusterCollection_); 
00213   // create the super cluster Ref
00214   createSuperClusterGsfMapRefs(scRefProd);
00215 
00216   // Now create the GsfElectronCoers
00217   createGsfElectronCores(*gsfElectronCores_p);
00218   // Put them in the as to get to be able to build a Ref
00219   const edm::OrphanHandle<reco::GsfElectronCoreCollection> gsfElectronCoreRefProd = 
00220     iEvent.put(gsfElectronCores_p,GsfElectronCoreCollection_);
00221 
00222   // now create the Refs 
00223   createGsfElectronCoreRefs(gsfElectronCoreRefProd);
00224 
00225   // now make the GsfElectron
00226   createGsfElectrons(*pfCandidates,isolationValues,*gsfElectrons_p);
00227   iEvent.put(gsfElectrons_p,GsfElectronCollection_);
00228 
00229   fillMVAValueMap(iEvent,mvaFiller);
00230   mvaFiller.fill();
00231 
00232   fillSCRefValueMap(iEvent,scRefFiller);
00233   scRefFiller.fill();
00234 
00235   // MVA map
00236   iEvent.put(mvaMap_p,PFMVAValueMap_);
00237   // Gsf-SC map
00238   iEvent.put(scMap_p,PFSCValueMap_);
00239 
00240 
00241   
00242 }
00243 
00244 
00245 
00246 bool PFElectronTranslator::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c, 
00247                                               const edm::InputTag& tag, 
00248                                               const edm::Event& iEvent) const {  
00249   bool found = iEvent.getByLabel(tag, c);
00250 
00251   if(!found && !emptyIsOk_)
00252     {
00253       std::ostringstream  err;
00254       err<<" cannot get PFCandidates: "
00255          <<tag<<std::endl;
00256       edm::LogError("PFElectronTranslator")<<err.str();
00257     }
00258   return found;
00259       
00260 }
00261 
00262 void PFElectronTranslator::fetchGsfCollection(edm::Handle<reco::GsfTrackCollection>& c, 
00263                                               const edm::InputTag& tag, 
00264                                               const edm::Event& iEvent) const {  
00265   bool found = iEvent.getByLabel(tag, c);
00266   
00267   if(!found ) {
00268     std::ostringstream  err;
00269     err<<" cannot get GSFTracks: "
00270        <<tag<<std::endl;
00271     edm::LogError("PFElectronTranslator")<<err.str();
00272     throw cms::Exception( "MissingProduct", err.str());
00273   }  
00274 }
00275 
00276 // The basic cluster is a copy of the PFCluster -> the energy is not corrected 
00277 // It should be possible to get the corrected energy (including the associated PS energy)
00278 // from the PFCandidate daugthers ; Needs some work 
00279 void PFElectronTranslator::createBasicCluster(const reco::PFBlockElement & PFBE, 
00280                                               reco::BasicClusterCollection & basicClusters, 
00281                                               std::vector<const reco::PFCluster *> & pfClusters,
00282                                               const reco::PFCandidate & coCandidate) const
00283 {
00284   reco::PFClusterRef myPFClusterRef= PFBE.clusterRef();
00285   if(myPFClusterRef.isNull()) return;  
00286 
00287   const reco::PFCluster & myPFCluster (*myPFClusterRef);
00288   pfClusters.push_back(&myPFCluster);
00289 //  std::cout << " Creating BC " << myPFCluster.energy() << " " << coCandidate.ecalEnergy() <<" "<<  coCandidate.rawEcalEnergy() <<std::endl;
00290 //  std::cout << " # hits " << myPFCluster.hitsAndFractions().size() << std::endl;
00291 
00292 //  basicClusters.push_back(reco::CaloCluster(myPFCluster.energy(),
00293   basicClusters.push_back(reco::CaloCluster(
00294                                             //      myPFCluster.energy(),
00295                                             coCandidate.rawEcalEnergy(),
00296                                             myPFCluster.position(),
00297                                             myPFCluster.caloID(),
00298                                             myPFCluster.hitsAndFractions(),
00299                                             myPFCluster.algo(),
00300                                             myPFCluster.seed()));
00301 }
00302 
00303 
00304 void PFElectronTranslator::createPreshowerCluster(const reco::PFBlockElement & PFBE, reco::PreshowerClusterCollection& preshowerClusters,unsigned plane) const
00305 {
00306   reco::PFClusterRef  myPFClusterRef= PFBE.clusterRef();
00307   preshowerClusters.push_back(reco::PreshowerCluster(myPFClusterRef->energy(),myPFClusterRef->position(),
00308                                                myPFClusterRef->hitsAndFractions(),plane));
00309 }
00310 
00311 void PFElectronTranslator::createBasicClusterPtrs(const edm::OrphanHandle<reco::BasicClusterCollection> & basicClustersHandle )
00312 {
00313   unsigned size=GsfTrackRef_.size();
00314   unsigned basicClusterCounter=0;
00315   basicClusterPtr_.resize(size);
00316 
00317   for(unsigned iGSF=0;iGSF<size;++iGSF) // loop on tracks
00318     {
00319       unsigned nbc=basicClusters_[iGSF].size();
00320       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00321         {
00322           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00323           reco::CaloClusterPtr bcPtr(basicClustersHandle,basicClusterCounter);
00324           basicClusterPtr_[iGSF].push_back(bcPtr);
00325           ++basicClusterCounter;
00326         }
00327     }
00328 }
00329 
00330 void PFElectronTranslator::createPreshowerClusterPtrs(const edm::OrphanHandle<reco::PreshowerClusterCollection> & preshowerClustersHandle )
00331 {
00332   unsigned size=GsfTrackRef_.size();
00333   unsigned psClusterCounter=0;
00334   preshowerClusterPtr_.resize(size);
00335 
00336   for(unsigned iGSF=0;iGSF<size;++iGSF) // loop on tracks
00337     {
00338       unsigned nbc=preshowerClusters_[iGSF].size();
00339       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00340         {
00341           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00342           reco::CaloClusterPtr psPtr(preshowerClustersHandle,psClusterCounter);
00343           preshowerClusterPtr_[iGSF].push_back(psPtr);
00344           ++psClusterCounter;
00345         }
00346     }
00347 }
00348 
00349 void PFElectronTranslator::createSuperClusterGsfMapRefs(const edm::OrphanHandle<reco::SuperClusterCollection> & superClustersHandle )
00350 {
00351   unsigned size=GsfTrackRef_.size();
00352 
00353   for(unsigned iGSF=0;iGSF<size;++iGSF) // loop on tracks
00354     {
00355       edm::Ref<reco::SuperClusterCollection> scRef(superClustersHandle,iGSF);
00356       scMap_[GsfTrackRef_[iGSF]]=scRef;
00357     }
00358 }
00359 
00360 
00361 void PFElectronTranslator::fillMVAValueMap(edm::Event& iEvent, edm::ValueMap<float>::Filler & filler) 
00362 {
00363   gsfMvaMap_.clear();
00364   edm::Handle<reco::PFCandidateCollection> pfCandidates;
00365   bool status=fetchCandidateCollection(pfCandidates, 
00366                                        inputTagPFCandidateElectrons_, 
00367                                        iEvent );
00368   
00369   unsigned ncand=(status)?pfCandidates->size():0;
00370   for( unsigned i=0; i<ncand; ++i ) {
00371     
00372     const reco::PFCandidate& cand = (*pfCandidates)[i];    
00373     if(cand.particleId()!=reco::PFCandidate::e) continue; 
00374     if(cand.gsfTrackRef().isNull()) continue;
00375     // Fill the MVA map
00376     gsfMvaMap_[cand.gsfTrackRef()]=cand.mva_e_pi();       
00377   }
00378   
00379   edm::Handle<reco::GsfTrackCollection> gsfTracks;
00380   fetchGsfCollection(gsfTracks,
00381                      inputTagGSFTracks_,
00382                      iEvent);
00383   unsigned ngsf=gsfTracks->size();
00384   std::vector<float> values;
00385   for(unsigned igsf=0;igsf<ngsf;++igsf)
00386     {
00387       reco::GsfTrackRef theTrackRef(gsfTracks, igsf);
00388       std::map<reco::GsfTrackRef,float>::const_iterator itcheck=gsfMvaMap_.find(theTrackRef);
00389       if(itcheck==gsfMvaMap_.end())
00390         {
00391           //      edm::LogWarning("PFElectronTranslator") << "MVA Map, missing GSF track ref " << std::endl;
00392           values.push_back(-99.);
00393           //      std::cout << " Push_back -99. " << std::endl;
00394         }
00395       else
00396         {
00397           //      std::cout <<  " Value " << itcheck->second << std::endl;
00398           values.push_back(itcheck->second);      
00399         }
00400     }
00401   filler.insert(gsfTracks,values.begin(),values.end());
00402 }
00403 
00404 
00405 void PFElectronTranslator::fillSCRefValueMap(edm::Event& iEvent, 
00406                                              edm::ValueMap<reco::SuperClusterRef>::Filler & filler) const
00407 {
00408   edm::Handle<reco::GsfTrackCollection> gsfTracks;
00409   fetchGsfCollection(gsfTracks,
00410                      inputTagGSFTracks_,
00411                      iEvent);
00412   unsigned ngsf=gsfTracks->size();
00413   std::vector<reco::SuperClusterRef> values;
00414   for(unsigned igsf=0;igsf<ngsf;++igsf)
00415     {
00416       reco::GsfTrackRef theTrackRef(gsfTracks, igsf);
00417       std::map<reco::GsfTrackRef,reco::SuperClusterRef>::const_iterator itcheck=scMap_.find(theTrackRef);
00418       if(itcheck==scMap_.end())
00419         {
00420           //      edm::LogWarning("PFElectronTranslator") << "SCRef Map, missing GSF track ref" << std::endl;
00421           values.push_back(reco::SuperClusterRef());
00422         }
00423       else
00424         {
00425           values.push_back(itcheck->second);      
00426         }
00427     }
00428   filler.insert(gsfTracks,values.begin(),values.end());
00429 }
00430 
00431 
00432 void PFElectronTranslator::createSuperClusters(const reco::PFCandidateCollection & pfCand,
00433                                                reco::SuperClusterCollection &superClusters) const
00434 {
00435   unsigned nGSF=GsfTrackRef_.size();
00436   for(unsigned iGSF=0;iGSF<nGSF;++iGSF)
00437     {
00438 
00439       // Computes energy position a la e/gamma 
00440       double sclusterE=0;
00441       double posX=0.;
00442       double posY=0.;
00443       double posZ=0.;
00444       
00445       unsigned nbasics=basicClusters_[iGSF].size();
00446       for(unsigned ibc=0;ibc<nbasics;++ibc)
00447         {
00448           double e = basicClusters_[iGSF][ibc].energy();
00449           sclusterE += e;
00450           posX += e * basicClusters_[iGSF][ibc].position().X();
00451           posY += e * basicClusters_[iGSF][ibc].position().Y();
00452           posZ += e * basicClusters_[iGSF][ibc].position().Z();   
00453         }
00454       posX /=sclusterE;
00455       posY /=sclusterE;
00456       posZ /=sclusterE;
00457       
00458       if(pfCand[gsfPFCandidateIndex_[iGSF]].gsfTrackRef()!=GsfTrackRef_[iGSF])
00459         {
00460           edm::LogError("PFElectronTranslator") << " Major problem in PFElectron Translator" << std::endl;
00461         }
00462       
00463       // compute the width
00464       PFClusterWidthAlgo pfwidth(pfClusters_[iGSF]);
00465       
00466       double correctedEnergy=pfCand[gsfPFCandidateIndex_[iGSF]].ecalEnergy();
00467       reco::SuperCluster mySuperCluster(correctedEnergy,math::XYZPoint(posX,posY,posZ));
00468       // protection against empty basic cluster collection ; the value is -2 in this case
00469       if(nbasics)
00470         {
00471 //        std::cout << "SuperCluster creation; energy " << pfCand[gsfPFCandidateIndex_[iGSF]].ecalEnergy();
00472 //        std::cout << " " <<   pfCand[gsfPFCandidateIndex_[iGSF]].rawEcalEnergy() << std::endl;
00473 //        std::cout << "Seed energy from basic " << basicClusters_[iGSF][0].energy() << std::endl;
00474           mySuperCluster.setSeed(basicClusterPtr_[iGSF][0]);
00475         }
00476       else
00477         {
00478           //      std::cout << "SuperCluster creation ; seed energy " << 0 << std::endl;
00479 //        std::cout << "SuperCluster creation ; energy " << pfCand[gsfPFCandidateIndex_[iGSF]].ecalEnergy();
00480 //        std::cout << " " <<   pfCand[gsfPFCandidateIndex_[iGSF]].rawEcalEnergy() << std::endl;
00481 //        std::cout << " No seed found " << 0 << std::endl;       
00482 //        std::cout << " MVA " << pfCand[gsfPFCandidateIndex_[iGSF]].mva_e_pi() << std::endl;
00483           mySuperCluster.setSeed(reco::CaloClusterPtr());
00484         }
00485       // the seed should be the first basic cluster
00486 
00487       for(unsigned ibc=0;ibc<nbasics;++ibc)
00488         {
00489           mySuperCluster.addCluster(basicClusterPtr_[iGSF][ibc]);
00490           //      std::cout <<"Adding Ref to SC " << basicClusterPtr_[iGSF][ibc].index() << std::endl;
00491           const std::vector< std::pair<DetId, float> > & v1 = basicClusters_[iGSF][ibc].hitsAndFractions();
00492           //      std::cout << " Number of cells " << v1.size() << std::endl;
00493           for( std::vector< std::pair<DetId, float> >::const_iterator diIt = v1.begin();
00494                diIt != v1.end();
00495                ++diIt ) {
00496             //      std::cout << " Adding DetId " << (diIt->first).rawId() << " " << diIt->second << std::endl;
00497             mySuperCluster.addHitAndFraction(diIt->first,diIt->second);
00498           } // loop over rechits      
00499         }      
00500 
00501       unsigned nps=preshowerClusterPtr_[iGSF].size();
00502       for(unsigned ips=0;ips<nps;++ips)
00503         {
00504           mySuperCluster.addPreshowerCluster(preshowerClusterPtr_[iGSF][ips]);
00505         }
00506       
00507 
00508       // Set the preshower energy
00509       mySuperCluster.setPreshowerEnergy(pfCand[gsfPFCandidateIndex_[iGSF]].pS1Energy()+
00510                                         pfCand[gsfPFCandidateIndex_[iGSF]].pS2Energy());
00511 
00512       // Set the cluster width
00513       mySuperCluster.setEtaWidth(pfwidth.pflowEtaWidth());
00514       mySuperCluster.setPhiWidth(pfwidth.pflowPhiWidth());
00515       // Force the computation of rawEnergy_ of the reco::SuperCluster
00516       mySuperCluster.rawEnergy();
00517       superClusters.push_back(mySuperCluster);
00518    }
00519 }
00520 
00521 
00522 const reco::PFCandidate & PFElectronTranslator::correspondingDaughterCandidate(const reco::PFCandidate & cand, const reco::PFBlockElement & pfbe) const
00523 {
00524   unsigned refindex=pfbe.index();
00525   //  std::cout << " N daughters " << cand.numberOfDaughters() << std::endl;
00526   reco::PFCandidate::const_iterator myDaughterCandidate=cand.begin();
00527   reco::PFCandidate::const_iterator itend=cand.end();
00528 
00529   for(;myDaughterCandidate!=itend;++myDaughterCandidate)
00530     {
00531       const reco::PFCandidate * myPFCandidate = (const reco::PFCandidate*)&*myDaughterCandidate;
00532       if(myPFCandidate->elementsInBlocks().size()!=1)
00533         {
00534           //      std::cout << " Daughter with " << myPFCandidate.elementsInBlocks().size()<< " element in block " << std::endl;
00535           return cand;
00536         }
00537       if(myPFCandidate->elementsInBlocks()[0].second==refindex) 
00538         {
00539           //      std::cout << " Found it " << cand << std::endl;
00540           return *myPFCandidate;
00541         }      
00542     }
00543   return cand;
00544 }
00545 
00546 void PFElectronTranslator::createGsfElectronCores(reco::GsfElectronCoreCollection & gsfElectronCores) const {
00547   unsigned nGSF=GsfTrackRef_.size();
00548   for(unsigned iGSF=0;iGSF<nGSF;++iGSF)
00549     {
00550       reco::GsfElectronCore myElectronCore(GsfTrackRef_[iGSF]);
00551       myElectronCore.setCtfTrack(kfTrackRef_[iGSF],-1.);
00552       std::map<reco::GsfTrackRef,reco::SuperClusterRef>::const_iterator 
00553         itcheck=scMap_.find(GsfTrackRef_[iGSF]);
00554       if(itcheck!=scMap_.end())
00555         myElectronCore.setPflowSuperCluster(itcheck->second);
00556       gsfElectronCores.push_back(myElectronCore);
00557     }
00558 }
00559 
00560 void PFElectronTranslator::createGsfElectronCoreRefs(const edm::OrphanHandle<reco::GsfElectronCoreCollection> & gsfElectronCoreHandle) {
00561   unsigned size=GsfTrackRef_.size();
00562   
00563   for(unsigned iGSF=0;iGSF<size;++iGSF) // loop on tracks
00564     {
00565       edm::Ref<reco::GsfElectronCoreCollection> elecCoreRef(gsfElectronCoreHandle,iGSF);
00566       gsfElectronCoreRefs_.push_back(elecCoreRef);
00567     }  
00568 }
00569 
00570 void PFElectronTranslator::getAmbiguousGsfTracks(const reco::PFBlockElement & PFBE, std::vector<reco::GsfTrackRef>& tracks) const {
00571   const reco::PFBlockElementGsfTrack *  GsfEl =  dynamic_cast<const reco::PFBlockElementGsfTrack*>(&PFBE);
00572   if(GsfEl==0) return;
00573   const std::vector<reco::GsfPFRecTrackRef>& ambPFRecTracks(GsfEl->GsftrackRefPF()->convBremGsfPFRecTrackRef());
00574   unsigned ntracks=ambPFRecTracks.size();
00575   for(unsigned it=0;it<ntracks;++it) {
00576     tracks.push_back(ambPFRecTracks[it]->gsfTrackRef());
00577   }
00578 }
00579 
00580 
00581 void PFElectronTranslator::createGsfElectrons(const reco::PFCandidateCollection & pfcand, 
00582                                               const IsolationValueMaps& isolationValues,
00583                                               reco::GsfElectronCollection &gsfelectrons) {
00584   unsigned size=GsfTrackRef_.size();
00585   
00586   for(unsigned iGSF=0;iGSF<size;++iGSF) // loop on tracks
00587     {
00588       const reco::PFCandidate& pfCandidate(pfcand[gsfPFCandidateIndex_[iGSF]]);
00589       // Electron
00590       reco::GsfElectron myElectron(gsfElectronCoreRefs_[iGSF]);
00591       // Warning set p4 error ! 
00592       myElectron.setP4(reco::GsfElectron::P4_PFLOW_COMBINATION, pfCandidate.p4(),pfCandidate.deltaP(), true);
00593       
00594       // MVA inputs
00595       reco::GsfElectron::MvaInput myMvaInput;
00596       myMvaInput.earlyBrem = pfCandidate.electronExtraRef()->mvaVariable(reco::PFCandidateElectronExtra::MVA_FirstBrem);
00597       myMvaInput.lateBrem = pfCandidate.electronExtraRef()->mvaVariable(reco::PFCandidateElectronExtra::MVA_LateBrem);
00598       myMvaInput.deltaEta = pfCandidate.electronExtraRef()->mvaVariable(reco::PFCandidateElectronExtra::MVA_DeltaEtaTrackCluster);
00599       myMvaInput.sigmaEtaEta = pfCandidate.electronExtraRef()->sigmaEtaEta();
00600       myMvaInput.hadEnergy = pfCandidate.electronExtraRef()->hadEnergy();
00601 
00602       // Mustache
00603       reco::Mustache myMustache;
00604       myMustache.MustacheID(*(myElectron. pflowSuperCluster()), myMvaInput.nClusterOutsideMustache, myMvaInput.etOutsideMustache );
00605 
00606       myElectron.setMvaInput(myMvaInput);
00607 
00608       // MVA output
00609       reco::GsfElectron::MvaOutput myMvaOutput;
00610       myMvaOutput.status = pfCandidate.electronExtraRef()->electronStatus();
00611       myMvaOutput.mva = pfCandidate.mva_e_pi();
00612       myElectron.setMvaOutput(myMvaOutput);
00613       
00614       // ambiguous tracks
00615       unsigned ntracks=ambiguousGsfTracks_[iGSF].size();
00616       for(unsigned it=0;it<ntracks;++it) {
00617         myElectron.addAmbiguousGsfTrack(ambiguousGsfTracks_[iGSF][it]);
00618       }
00619 
00620       // isolation
00621       reco::GsfElectron::PflowIsolationVariables myPFIso;
00622       myPFIso.chargedHadronIso=(*isolationValues[0])[CandidatePtr_[iGSF]];
00623       myPFIso.photonIso=(*isolationValues[1])[CandidatePtr_[iGSF]];
00624       myPFIso.neutralHadronIso=(*isolationValues[2])[CandidatePtr_[iGSF]];      
00625       myElectron.setPfIsolationVariables(myPFIso);
00626       
00627       gsfelectrons.push_back(myElectron);
00628     }
00629 
00630 }
00631 
00632