Reference Manual

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/PFPhotonTranslator.h"
00006 
00007 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
00008 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00009 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00010 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
00011 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00012 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
00013 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00014 #include "DataFormats/EgammaCandidates/interface/PhotonCoreFwd.h"
00015 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00016 #include "DataFormats/VertexReco/interface/Vertex.h"
00017 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00018 //#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00019 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
00020 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
00021 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
00022 
00023 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00024 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
00025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00026 
00027 //#include "Geometry/Records/interface/CaloGeometryRecord.h"
00028 //#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00029 //#include "Geometry/CaloTopology/interface/CaloTopology.h"
00030 //#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00031 //#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00032 //#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00033 //#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
00034 //#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00035 //#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
00036 //#include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
00037 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
00038 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h" 
00039 //#include "CondFormats/EcalObjects/interface/EcalFunctionParameters.h" 
00040 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00041 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidatePhotonExtra.h"
00042 
00043 #include "DataFormats/Math/interface/Vector3D.h"
00044 #include "DataFormats/Math/interface/LorentzVector.h"
00045 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00046 #include "DataFormats/Common/interface/AssociationVector.h"
00047 
00048 #include <Math/VectorUtil.h>
00049 #include <vector>
00050 #include "TLorentzVector.h"
00051 #include "TMath.h"
00052 
00053 using namespace edm;
00054 using namespace std;
00055 using namespace reco;
00056 
00057 using namespace ROOT::Math::VectorUtil;
00058 typedef math::XYZTLorentzVector LorentzVector;
00059 typedef math::XYZPoint Point;
00060 typedef math::XYZVector Vector;
00061 
00062 
00063 PFPhotonTranslator::PFPhotonTranslator(const edm::ParameterSet & iConfig) {
00064 
00065   //std::cout << "PFPhotonTranslator" << std::endl;
00066 
00067   inputTagPFCandidates_ 
00068     = iConfig.getParameter<edm::InputTag>("PFCandidate");
00069   
00070   edm::ParameterSet isoVals  = iConfig.getParameter<edm::ParameterSet> ("isolationValues");
00071   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfChargedHadrons"));
00072   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfPhotons"));
00073   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfNeutralHadrons"));
00074   
00075 
00076   PFBasicClusterCollection_ = iConfig.getParameter<std::string>("PFBasicClusters");
00077   PFPreshowerClusterCollection_ = iConfig.getParameter<std::string>("PFPreshowerClusters");
00078   PFSuperClusterCollection_ = iConfig.getParameter<std::string>("PFSuperClusters");
00079   PFConversionCollection_ = iConfig.getParameter<std::string>("PFConversionCollection");
00080   PFPhotonCoreCollection_ = iConfig.getParameter<std::string>("PFPhotonCores");
00081   PFPhotonCollection_ = iConfig.getParameter<std::string>("PFPhotons");
00082 
00083   EGPhotonCollection_ = iConfig.getParameter<std::string>("EGPhotons");
00084 
00085   vertexProducer_   = iConfig.getParameter<std::string>("primaryVertexProducer");
00086 
00087   barrelEcalHits_   = iConfig.getParameter<edm::InputTag>("barrelEcalHits");
00088   endcapEcalHits_   = iConfig.getParameter<edm::InputTag>("endcapEcalHits");
00089 
00090   hcalTowers_ = iConfig.getParameter<edm::InputTag>("hcalTowers");
00091   hOverEConeSize_   = iConfig.getParameter<double>("hOverEConeSize");
00092 
00093   if (iConfig.exists("emptyIsOk")) emptyIsOk_ = iConfig.getParameter<bool>("emptyIsOk");
00094   else emptyIsOk_=false;
00095 
00096   produces<reco::BasicClusterCollection>(PFBasicClusterCollection_); 
00097   produces<reco::PreshowerClusterCollection>(PFPreshowerClusterCollection_); 
00098   produces<reco::SuperClusterCollection>(PFSuperClusterCollection_); 
00099   produces<reco::PhotonCoreCollection>(PFPhotonCoreCollection_);
00100   produces<reco::PhotonCollection>(PFPhotonCollection_); 
00101   produces<reco::ConversionCollection>(PFConversionCollection_);
00102 }
00103 
00104 PFPhotonTranslator::~PFPhotonTranslator() {}
00105 
00106 void PFPhotonTranslator::beginRun(edm::Run& run,const edm::EventSetup & es) {
00107 
00108 }
00109 
00110 void PFPhotonTranslator::produce(edm::Event& iEvent,  
00111                                     const edm::EventSetup& iSetup) { 
00112 
00113   //cout << "NEW EVENT"<<endl;
00114 
00115   std::auto_ptr<reco::BasicClusterCollection> 
00116     basicClusters_p(new reco::BasicClusterCollection);
00117 
00118   std::auto_ptr<reco::PreshowerClusterCollection>
00119     psClusters_p(new reco::PreshowerClusterCollection);
00120 
00121   /*
00122   std::auto_ptr<reco::ConversionCollection>
00123     SingleLeg_p(new reco::ConversionCollection);
00124   */
00125 
00126   reco::SuperClusterCollection outputSuperClusterCollection;
00127   reco::ConversionCollection outputOneLegConversionCollection;
00128   reco::PhotonCoreCollection outputPhotonCoreCollection;
00129   reco::PhotonCollection outputPhotonCollection;
00130 
00131   outputSuperClusterCollection.clear();
00132   outputOneLegConversionCollection.clear();
00133   outputPhotonCoreCollection.clear();
00134   outputPhotonCollection.clear();
00135 
00136 
00137   edm::Handle<reco::PFCandidateCollection> pfCandidates;
00138   bool status=fetchCandidateCollection(pfCandidates, 
00139                                        inputTagPFCandidates_, 
00140                                        iEvent );
00141 
00142   edm::Handle<reco::PhotonCollection> egPhotons;
00143   iEvent.getByLabel(EGPhotonCollection_, egPhotons);
00144   
00145 
00146   Handle<reco::VertexCollection> vertexHandle;
00147 
00148   
00149   IsolationValueMaps isolationValues(inputTagIsoVals_.size());
00150   for (size_t j = 0; j<inputTagIsoVals_.size(); ++j) {
00151     iEvent.getByLabel(inputTagIsoVals_[j], isolationValues[j]);
00152   }
00153   
00154 
00155   // clear the vectors
00156   photPFCandidateIndex_.clear();
00157   basicClusters_.clear();
00158   pfClusters_.clear();
00159   preshowerClusters_.clear();
00160   superClusters_.clear();
00161   basicClusterPtr_.clear();
00162   preshowerClusterPtr_.clear();
00163   CandidatePtr_.clear();
00164   egSCRef_.clear();
00165   egPhotonRef_.clear();
00166 
00167   pfConv_.clear();
00168   pfSingleLegConv_.clear();
00169   pfSingleLegConvMva_.clear();
00170   conv1legPFCandidateIndex_.clear();
00171   conv2legPFCandidateIndex_.clear();
00172 
00173   // loop on the candidates 
00174   //CC@@
00175   // we need first to create AND put the SuperCluster, 
00176   // basic clusters and presh clusters collection 
00177   // in order to get a working Handle
00178   unsigned ncand=(status)?pfCandidates->size():0;
00179 
00180   unsigned iphot=0;
00181   unsigned iconv1leg=0;
00182   unsigned iconv2leg=0;
00183 
00184   for( unsigned i=0; i<ncand; ++i ) {
00185 
00186     const reco::PFCandidate& cand = (*pfCandidates)[i];    
00187     if(cand.particleId()!=reco::PFCandidate::gamma) continue;
00188     //cout << "cand.mva_nothing_gamma()="<<cand. mva_nothing_gamma()<<endl;
00189     if(cand. mva_nothing_gamma()>0.001)//Found PFPhoton with PFPhoton Extras saved
00190       {
00191 
00192         //cout << "NEW PHOTON" << endl;
00193 
00194         //std::cout << "nDoubleLegConv="<<cand.photonExtraRef()->conversionRef().size()<<std::endl;
00195 
00196         if (cand.photonExtraRef()->conversionRef().size()>0){
00197 
00198           pfConv_.push_back(reco::ConversionRefVector());
00199 
00200           const reco::ConversionRefVector & doubleLegConvColl = cand.photonExtraRef()->conversionRef();
00201           for (unsigned int iconv=0; iconv<doubleLegConvColl.size(); iconv++){
00202             pfConv_[iconv2leg].push_back(doubleLegConvColl[iconv]);
00203           }
00204 
00205           conv2legPFCandidateIndex_.push_back(iconv2leg);
00206           iconv2leg++;
00207         }
00208         else conv2legPFCandidateIndex_.push_back(-1);
00209 
00210         const std::vector<reco::TrackRef> & singleLegConvColl = cand.photonExtraRef()->singleLegConvTrackRef();
00211         const std::vector<float>& singleLegConvCollMva = cand.photonExtraRef()->singleLegConvMva();
00212         
00213         //std::cout << "nSingleLegConv=" <<singleLegConvColl.size() << std::endl;
00214 
00215         if (singleLegConvColl.size()>0){
00216 
00217           pfSingleLegConv_.push_back(std::vector<reco::TrackRef>());
00218           pfSingleLegConvMva_.push_back(std::vector<float>());
00219 
00220           //cout << "nTracks="<< singleLegConvColl.size()<<endl;
00221           for (unsigned int itk=0; itk<singleLegConvColl.size(); itk++){
00222             //cout << "Track pt="<< singleLegConvColl[itk]->pt() <<endl;
00223 
00224             pfSingleLegConv_[iconv1leg].push_back(singleLegConvColl[itk]);
00225             pfSingleLegConvMva_[iconv1leg].push_back(singleLegConvCollMva[itk]);
00226           }
00227 
00228 
00229           conv1legPFCandidateIndex_.push_back(iconv1leg);
00230         
00231           iconv1leg++;
00232         }
00233         else conv1legPFCandidateIndex_.push_back(-1);
00234         
00235       }
00236 
00237     photPFCandidateIndex_.push_back(i);
00238 
00239     basicClusters_.push_back(reco::BasicClusterCollection());
00240     pfClusters_.push_back(std::vector<const reco::PFCluster *>());
00241     preshowerClusters_.push_back(reco::PreshowerClusterCollection());
00242     superClusters_.push_back(reco::SuperClusterCollection());
00243 
00244     reco::PFCandidatePtr ptrToPFPhoton(pfCandidates,i);
00245     CandidatePtr_.push_back(ptrToPFPhoton);  
00246     egSCRef_.push_back(cand.superClusterRef());
00247     //std::cout << "PFPhoton cand " << iphot << std::endl;
00248 
00249     int iegphot=0;
00250     for (reco::PhotonCollection::const_iterator gamIter = egPhotons->begin(); gamIter != egPhotons->end(); ++gamIter){
00251       if (cand.superClusterRef()==gamIter->superCluster()){
00252         reco::PhotonRef PhotRef(reco::PhotonRef(egPhotons, iegphot));
00253         egPhotonRef_.push_back(PhotRef);
00254       }
00255       iegphot++;
00256     }
00257 
00258 
00259     //std::cout << "Cand elements in blocks : " << cand.elementsInBlocks().size() << std::endl;
00260 
00261     for(unsigned iele=0; iele<cand.elementsInBlocks().size(); ++iele) {
00262       // first get the block 
00263       reco::PFBlockRef blockRef = cand.elementsInBlocks()[iele].first;
00264       //
00265       unsigned elementIndex = cand.elementsInBlocks()[iele].second;
00266       // check it actually exists 
00267       if(blockRef.isNull()) continue;
00268       
00269       // then get the elements of the block
00270       const edm::OwnVector< reco::PFBlockElement >&  elements = (*blockRef).elements();
00271       
00272       const reco::PFBlockElement & pfbe (elements[elementIndex]); 
00273       // The first ECAL element should be the cluster associated to the GSF; defined as the seed
00274       if(pfbe.type()==reco::PFBlockElement::ECAL)
00275         {         
00276 
00277           //std::cout << "BlockElement ECAL" << std::endl;
00278           // the Brem photons are saved as daughter PFCandidate; this 
00279           // is convenient to access the corrected energy
00280           //      std::cout << " Found candidate "  << correspondingDaughterCandidate(coCandidate,pfbe) << " " << coCandidate << std::endl;
00281           createBasicCluster(pfbe,basicClusters_[iphot],pfClusters_[iphot],correspondingDaughterCandidate(cand,pfbe));
00282         }
00283       if(pfbe.type()==reco::PFBlockElement::PS1)
00284         {
00285           //std::cout << "BlockElement PS1" << std::endl;
00286           createPreshowerCluster(pfbe,preshowerClusters_[iphot],1);
00287         }
00288       if(pfbe.type()==reco::PFBlockElement::PS2)
00289         {
00290           //std::cout << "BlockElement PS2" << std::endl;
00291           createPreshowerCluster(pfbe,preshowerClusters_[iphot],2);
00292         }    
00293 
00294 
00295     }   // loop on the elements
00296 
00297         // save the basic clusters
00298     basicClusters_p->insert(basicClusters_p->end(),basicClusters_[iphot].begin(), basicClusters_[iphot].end());
00299     // save the preshower clusters
00300     psClusters_p->insert(psClusters_p->end(),preshowerClusters_[iphot].begin(),preshowerClusters_[iphot].end());
00301 
00302     ++iphot;
00303 
00304   } // loop on PFCandidates
00305 
00306 
00307    //Save the basic clusters and get an handle as to be able to create valid Refs (thanks to Claude)
00308   //  std::cout << " Number of basic clusters " << basicClusters_p->size() << std::endl;
00309   const edm::OrphanHandle<reco::BasicClusterCollection> bcRefProd = 
00310     iEvent.put(basicClusters_p,PFBasicClusterCollection_);
00311 
00312   //preshower clusters
00313   const edm::OrphanHandle<reco::PreshowerClusterCollection> psRefProd = 
00314     iEvent.put(psClusters_p,PFPreshowerClusterCollection_);
00315   
00316   // now that the Basic clusters are in the event, the Ref can be created
00317   createBasicClusterPtrs(bcRefProd);
00318   // now that the preshower clusters are in the event, the Ref can be created
00319   createPreshowerClusterPtrs(psRefProd);
00320 
00321   // and now the Super cluster can be created with valid references  
00322   //if(status) createSuperClusters(*pfCandidates,*superClusters_p);
00323   if(status) createSuperClusters(*pfCandidates,outputSuperClusterCollection);
00324   
00325   //std::cout << "nb superclusters in collection : "<<outputSuperClusterCollection.size()<<std::endl;
00326 
00327   // Let's put the super clusters in the event
00328   std::auto_ptr<reco::SuperClusterCollection> superClusters_p(new reco::SuperClusterCollection(outputSuperClusterCollection));  
00329   const edm::OrphanHandle<reco::SuperClusterCollection> scRefProd = iEvent.put(superClusters_p,PFSuperClusterCollection_); 
00330 
00331 
00332   /*
00333   int ipho=0;
00334   for (reco::SuperClusterCollection::const_iterator gamIter = scRefProd->begin(); gamIter != scRefProd->end(); ++gamIter){
00335     std::cout << "SC i="<<ipho<<" energy="<<gamIter->energy()<<std::endl;
00336     ipho++;
00337   }
00338   */
00339 
00340 
00341   //1-leg conversions
00342 
00343 
00344   if (status) createOneLegConversions(scRefProd, outputOneLegConversionCollection);
00345 
00346 
00347   std::auto_ptr<reco::ConversionCollection> SingleLeg_p(new reco::ConversionCollection(outputOneLegConversionCollection));  
00348   const edm::OrphanHandle<reco::ConversionCollection> ConvRefProd = iEvent.put(SingleLeg_p,PFConversionCollection_);
00349   /*
00350   int iconv = 0;
00351   for (reco::ConversionCollection::const_iterator convIter = ConvRefProd->begin(); convIter != ConvRefProd->end(); ++convIter){
00352 
00353     std::cout << "OneLegConv i="<<iconv<<" nTracks="<<convIter->nTracks()<<" EoverP="<<convIter->EoverP() <<std::endl;
00354     std::vector<edm::RefToBase<reco::Track> > convtracks = convIter->tracks();
00355     for (unsigned int itk=0; itk<convtracks.size(); itk++){
00356       std::cout << "Track pt="<< convtracks[itk]->pt() << std::endl;
00357     }  
00358 
00359     iconv++;
00360   }
00361   */
00362 
00363   //create photon cores
00364   //if(status) createPhotonCores(pfCandidates, scRefProd, *photonCores_p);
00365   if(status) createPhotonCores(scRefProd, ConvRefProd, outputPhotonCoreCollection);
00366   
00367   //std::cout << "nb photoncores in collection : "<<outputPhotonCoreCollection.size()<<std::endl;
00368 
00369   // Put the photon cores in the event
00370   std::auto_ptr<reco::PhotonCoreCollection> photonCores_p(new reco::PhotonCoreCollection(outputPhotonCoreCollection));  
00371   //std::cout << "photon core collection put in auto_ptr"<<std::endl;
00372   const edm::OrphanHandle<reco::PhotonCoreCollection> pcRefProd = iEvent.put(photonCores_p,PFPhotonCoreCollection_); 
00373   
00374   //std::cout << "photon core have been put in the event"<<std::endl;
00375   /*
00376   int ipho=0;
00377   for (reco::PhotonCoreCollection::const_iterator gamIter = pcRefProd->begin(); gamIter != pcRefProd->end(); ++gamIter){
00378     std::cout << "PhotonCore i="<<ipho<<" energy="<<gamIter->pfSuperCluster()->energy()<<std::endl;
00379     //for (unsigned int i=0; i<)
00380 
00381     std::cout << "PhotonCore i="<<ipho<<" nconv2leg="<<gamIter->conversions().size()<<" nconv1leg="<<gamIter->conversionsOneLeg().size()<<std::endl;
00382 
00383     const reco::ConversionRefVector & conv = gamIter->conversions();
00384     for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00385       cout << "2-leg iconv="<<iconv<<endl;
00386       cout << "2-leg nTracks="<<conv[iconv]->nTracks()<<endl;
00387       cout << "2-leg EoverP="<<conv[iconv]->EoverP()<<endl;
00388       cout << "2-leg ConvAlgorithm="<<conv[iconv]->algo()<<endl;
00389     }
00390 
00391     const reco::ConversionRefVector & convbis = gamIter->conversionsOneLeg();
00392     for (unsigned int iconv=0; iconv<convbis.size(); iconv++){
00393       cout << "1-leg iconv="<<iconv<<endl;
00394       cout << "1-leg nTracks="<<convbis[iconv]->nTracks()<<endl;
00395       cout << "1-leg EoverP="<<convbis[iconv]->EoverP()<<endl;
00396       cout << "1-leg ConvAlgorithm="<<convbis[iconv]->algo()<<endl;
00397     }
00398 
00399     ipho++;
00400   }
00401   */
00402 
00403   //load vertices
00404   reco::VertexCollection vertexCollection;
00405   bool validVertex=true;
00406   iEvent.getByLabel(vertexProducer_, vertexHandle);
00407   if (!vertexHandle.isValid()) {
00408     edm::LogWarning("PhotonProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
00409     validVertex=false;
00410   }
00411   if (validVertex) vertexCollection = *(vertexHandle.product());
00412 
00413   /*
00414   //load Ecal rechits
00415   bool validEcalRecHits=true;
00416   Handle<EcalRecHitCollection> barrelHitHandle;
00417   EcalRecHitCollection barrelRecHits;
00418   iEvent.getByLabel(barrelEcalHits_, barrelHitHandle);
00419   if (!barrelHitHandle.isValid()) {
00420     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
00421     validEcalRecHits=false; 
00422   }
00423   if (  validEcalRecHits)  barrelRecHits = *(barrelHitHandle.product());
00424   
00425   Handle<EcalRecHitCollection> endcapHitHandle;
00426   iEvent.getByLabel(endcapEcalHits_, endcapHitHandle);
00427   EcalRecHitCollection endcapRecHits;
00428   if (!endcapHitHandle.isValid()) {
00429     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
00430     validEcalRecHits=false; 
00431   }
00432   if( validEcalRecHits) endcapRecHits = *(endcapHitHandle.product());
00433 
00434   //load detector topology & geometry
00435   iSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
00436 
00437   edm::ESHandle<CaloTopology> pTopology;
00438   iSetup.get<CaloTopologyRecord>().get(theCaloTopo_);
00439   const CaloTopology *topology = theCaloTopo_.product();
00440 
00441   // get Hcal towers collection 
00442   Handle<CaloTowerCollection> hcalTowersHandle;
00443   iEvent.getByLabel(hcalTowers_, hcalTowersHandle);
00444   */
00445 
00446   //create photon collection
00447   //if(status) createPhotons(vertexCollection, pcRefProd, topology, &barrelRecHits, &endcapRecHits, hcalTowersHandle, isolationValues, outputPhotonCollection);
00448   if(status) createPhotons(vertexCollection, egPhotons, pcRefProd, isolationValues, outputPhotonCollection);
00449 
00450   // Put the photons in the event
00451   std::auto_ptr<reco::PhotonCollection> photons_p(new reco::PhotonCollection(outputPhotonCollection));  
00452   //std::cout << "photon collection put in auto_ptr"<<std::endl;
00453   const edm::OrphanHandle<reco::PhotonCollection> photonRefProd = iEvent.put(photons_p,PFPhotonCollection_); 
00454   //std::cout << "photons have been put in the event"<<std::endl;
00455   
00456   /*
00457   ipho=0;
00458   for (reco::PhotonCollection::const_iterator gamIter = photonRefProd->begin(); gamIter != photonRefProd->end(); ++gamIter){
00459     std::cout << "Photon i="<<ipho<<" pfEnergy="<<gamIter->pfSuperCluster()->energy()<<std::endl;
00460     
00461     const reco::ConversionRefVector & conv = gamIter->conversions();
00462     cout << "conversions obtained : conv.size()="<< conv.size()<<endl;
00463     for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00464       cout << "iconv="<<iconv<<endl;
00465       cout << "nTracks="<<conv[iconv]->nTracks()<<endl;
00466       cout << "EoverP="<<conv[iconv]->EoverP()<<endl;
00467       
00468       cout << "Vtx x="<<conv[iconv]->conversionVertex().x() << " y="<< conv[iconv]->conversionVertex().y()<<" z="<<conv[iconv]->conversionVertex().z()<< endl;
00469       cout << "VtxError x=" << conv[iconv]->conversionVertex().xError() << endl;
00470       
00471       std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
00472       //cout << "nTracks="<< convtracks.size()<<endl;
00473       for (unsigned int itk=0; itk<convtracks.size(); itk++){
00474         double convtrackpt = convtracks[itk]->pt();
00475         const edm::RefToBase<reco::Track> & mytrack = convtracks[itk];
00476         cout << "Track pt="<< convtrackpt <<endl;
00477         cout << "Track origin="<<gamIter->conversionTrackOrigin(mytrack)<<endl;
00478       }
00479     }
00480     
00481     //1-leg
00482     const reco::ConversionRefVector & convbis = gamIter->conversionsOneLeg();
00483     cout << "conversions obtained : conv.size()="<< convbis.size()<<endl;
00484     for (unsigned int iconv=0; iconv<convbis.size(); iconv++){
00485       cout << "iconv="<<iconv<<endl;
00486       cout << "nTracks="<<convbis[iconv]->nTracks()<<endl;
00487       cout << "EoverP="<<convbis[iconv]->EoverP()<<endl;
00488       
00489       cout << "Vtx x="<<convbis[iconv]->conversionVertex().x() << " y="<< convbis[iconv]->conversionVertex().y()<<" z="<<convbis[iconv]->conversionVertex().z()<< endl;
00490       cout << "VtxError x=" << convbis[iconv]->conversionVertex().xError() << endl;
00491        
00492       std::vector<edm::RefToBase<reco::Track> > convtracks = convbis[iconv]->tracks();
00493       //cout << "nTracks="<< convtracks.size()<<endl;
00494       for (unsigned int itk=0; itk<convtracks.size(); itk++){
00495         double convtrackpt = convtracks[itk]->pt();
00496         cout << "Track pt="<< convtrackpt <<endl;
00497         cout << "Track origin="<<gamIter->conversionTrackOrigin((convtracks[itk]))<<endl;
00498       }
00499     }
00500     ipho++;
00501   }
00502   */
00503   
00504 }
00505 
00506 bool PFPhotonTranslator::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c, 
00507                                               const edm::InputTag& tag, 
00508                                               const edm::Event& iEvent) const {  
00509   bool found = iEvent.getByLabel(tag, c);
00510 
00511   if(!found && !emptyIsOk_)
00512     {
00513       std::ostringstream  err;
00514       err<<" cannot get PFCandidates: "
00515          <<tag<<std::endl;
00516       edm::LogError("PFPhotonTranslator")<<err.str();
00517     }
00518   return found;
00519       
00520 }
00521 
00522 // The basic cluster is a copy of the PFCluster -> the energy is not corrected 
00523 // It should be possible to get the corrected energy (including the associated PS energy)
00524 // from the PFCandidate daugthers ; Needs some work 
00525 void PFPhotonTranslator::createBasicCluster(const reco::PFBlockElement & PFBE, 
00526                                               reco::BasicClusterCollection & basicClusters, 
00527                                               std::vector<const reco::PFCluster *> & pfClusters,
00528                                               const reco::PFCandidate & coCandidate) const
00529 {
00530   reco::PFClusterRef myPFClusterRef = PFBE.clusterRef();
00531   if(myPFClusterRef.isNull()) return;  
00532 
00533   const reco::PFCluster & myPFCluster (*myPFClusterRef);
00534   pfClusters.push_back(&myPFCluster);
00535   //std::cout << " Creating BC " << myPFCluster.energy() << " " << coCandidate.ecalEnergy() <<" "<<  coCandidate.rawEcalEnergy() <<std::endl;
00536   //std::cout << " # hits " << myPFCluster.hitsAndFractions().size() << std::endl;
00537 
00538 //  basicClusters.push_back(reco::CaloCluster(myPFCluster.energy(),
00539   basicClusters.push_back(reco::CaloCluster(coCandidate.rawEcalEnergy(),
00540                                             myPFCluster.position(),
00541                                             myPFCluster.caloID(),
00542                                             myPFCluster.hitsAndFractions(),
00543                                             myPFCluster.algo(),
00544                                             myPFCluster.seed()));
00545 }
00546 
00547 void PFPhotonTranslator::createPreshowerCluster(const reco::PFBlockElement & PFBE, reco::PreshowerClusterCollection& preshowerClusters,unsigned plane) const
00548 {
00549   reco::PFClusterRef  myPFClusterRef= PFBE.clusterRef();
00550   preshowerClusters.push_back(reco::PreshowerCluster(myPFClusterRef->energy(),myPFClusterRef->position(),
00551                                                myPFClusterRef->hitsAndFractions(),plane));
00552 }
00553 
00554 void PFPhotonTranslator::createBasicClusterPtrs(const edm::OrphanHandle<reco::BasicClusterCollection> & basicClustersHandle )
00555 {
00556   unsigned size=photPFCandidateIndex_.size();
00557   unsigned basicClusterCounter=0;
00558   basicClusterPtr_.resize(size);
00559 
00560   for(unsigned iphot=0;iphot<size;++iphot) // loop on tracks
00561     {
00562       unsigned nbc=basicClusters_[iphot].size();
00563       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00564         {
00565           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00566           reco::CaloClusterPtr bcPtr(basicClustersHandle,basicClusterCounter);
00567           basicClusterPtr_[iphot].push_back(bcPtr);
00568           ++basicClusterCounter;
00569         }
00570     }
00571 }
00572 
00573 void PFPhotonTranslator::createPreshowerClusterPtrs(const edm::OrphanHandle<reco::PreshowerClusterCollection> & preshowerClustersHandle )
00574 {
00575   unsigned size=photPFCandidateIndex_.size();
00576   unsigned psClusterCounter=0;
00577   preshowerClusterPtr_.resize(size);
00578 
00579   for(unsigned iphot=0;iphot<size;++iphot) // loop on tracks
00580     {
00581       unsigned nbc=preshowerClusters_[iphot].size();
00582       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00583         {
00584           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00585           reco::CaloClusterPtr psPtr(preshowerClustersHandle,psClusterCounter);
00586           preshowerClusterPtr_[iphot].push_back(psPtr);
00587           ++psClusterCounter;
00588         }
00589     }
00590 }
00591 
00592 void PFPhotonTranslator::createSuperClusters(const reco::PFCandidateCollection & pfCand,
00593                                                reco::SuperClusterCollection &superClusters) const
00594 {
00595   unsigned nphot=photPFCandidateIndex_.size();
00596   for(unsigned iphot=0;iphot<nphot;++iphot)
00597     {
00598 
00599       //cout << "SC iphot=" << iphot << endl;
00600 
00601       // Computes energy position a la e/gamma 
00602       double sclusterE=0;
00603       double posX=0.;
00604       double posY=0.;
00605       double posZ=0.;
00606       
00607       unsigned nbasics=basicClusters_[iphot].size();
00608       for(unsigned ibc=0;ibc<nbasics;++ibc)
00609         {
00610           //cout << "BC in SC : iphot="<<iphot<<endl;
00611           
00612           double e = basicClusters_[iphot][ibc].energy();
00613           sclusterE += e;
00614           posX += e * basicClusters_[iphot][ibc].position().X();
00615           posY += e * basicClusters_[iphot][ibc].position().Y();
00616           posZ += e * basicClusters_[iphot][ibc].position().Z();          
00617         }
00618       posX /=sclusterE;
00619       posY /=sclusterE;
00620       posZ /=sclusterE;
00621       
00622       /*
00623       if(pfCand[gsfPFCandidateIndex_[iphot]].gsfTrackRef()!=GsfTrackRef_[iphot])
00624         {
00625           edm::LogError("PFElectronTranslator") << " Major problem in PFElectron Translator" << std::endl;
00626         }
00627       */      
00628 
00629       // compute the width
00630       PFClusterWidthAlgo pfwidth(pfClusters_[iphot]);
00631       
00632       double correctedEnergy=pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
00633       reco::SuperCluster mySuperCluster(correctedEnergy,math::XYZPoint(posX,posY,posZ));
00634       // protection against empty basic cluster collection ; the value is -2 in this case
00635       if(nbasics)
00636         {
00637 //        std::cout << "SuperCluster creation; energy " << pfCand[gsfPFCandidateIndex_[iphot]].ecalEnergy();
00638 //        std::cout << " " <<   pfCand[gsfPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
00639 //        std::cout << "Seed energy from basic " << basicClusters_[iphot][0].energy() << std::endl;
00640           mySuperCluster.setSeed(basicClusterPtr_[iphot][0]);
00641         }
00642       else
00643         {
00644           //      std::cout << "SuperCluster creation ; seed energy " << 0 << std::endl;
00645           //std::cout << "SuperCluster creation ; energy " << pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
00646           //std::cout << " " <<   pfCand[photPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
00647 //        std::cout << " No seed found " << 0 << std::endl;       
00648 //        std::cout << " MVA " << pfCand[gsfPFCandidateIndex_[iphot]].mva_e_pi() << std::endl;
00649           mySuperCluster.setSeed(reco::CaloClusterPtr());
00650         }
00651       // the seed should be the first basic cluster
00652 
00653       for(unsigned ibc=0;ibc<nbasics;++ibc)
00654         {
00655           mySuperCluster.addCluster(basicClusterPtr_[iphot][ibc]);
00656           //      std::cout <<"Adding Ref to SC " << basicClusterPtr_[iphot][ibc].index() << std::endl;
00657           const std::vector< std::pair<DetId, float> > & v1 = basicClusters_[iphot][ibc].hitsAndFractions();
00658           //      std::cout << " Number of cells " << v1.size() << std::endl;
00659           for( std::vector< std::pair<DetId, float> >::const_iterator diIt = v1.begin();
00660                diIt != v1.end();
00661                ++diIt ) {
00662             //      std::cout << " Adding DetId " << (diIt->first).rawId() << " " << diIt->second << std::endl;
00663             mySuperCluster.addHitAndFraction(diIt->first,diIt->second);
00664           } // loop over rechits      
00665         }      
00666 
00667       unsigned nps=preshowerClusterPtr_[iphot].size();
00668       for(unsigned ips=0;ips<nps;++ips)
00669         {
00670           mySuperCluster.addPreshowerCluster(preshowerClusterPtr_[iphot][ips]);
00671         }
00672       
00673 
00674       // Set the preshower energy
00675       mySuperCluster.setPreshowerEnergy(pfCand[photPFCandidateIndex_[iphot]].pS1Energy()+
00676                                         pfCand[photPFCandidateIndex_[iphot]].pS2Energy());
00677 
00678       // Set the cluster width
00679       mySuperCluster.setEtaWidth(pfwidth.pflowEtaWidth());
00680       mySuperCluster.setPhiWidth(pfwidth.pflowPhiWidth());
00681       // Force the computation of rawEnergy_ of the reco::SuperCluster
00682       mySuperCluster.rawEnergy();
00683 
00684       //cout << "SC energy="<< mySuperCluster.energy()<<endl;
00685 
00686       superClusters.push_back(mySuperCluster);
00687       //std::cout << "nb super clusters in collection : "<<superClusters.size()<<std::endl;
00688     }
00689 }
00690 
00691 void PFPhotonTranslator::createOneLegConversions(const edm::OrphanHandle<reco::SuperClusterCollection> & superClustersHandle, reco::ConversionCollection &oneLegConversions)
00692 {
00693 
00694   //std::cout << "createOneLegConversions" << std::endl;
00695 
00696     unsigned nphot=photPFCandidateIndex_.size();
00697     for(unsigned iphot=0;iphot<nphot;++iphot)
00698       {
00699 
00700         //if (conv1legPFCandidateIndex_[iphot]==-1) cout << "No OneLegConversions to add"<<endl;
00701         //else std::cout << "Phot "<<iphot<< " nOneLegConversions to add : "<<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size()<<endl;
00702 
00703 
00704         if (conv1legPFCandidateIndex_[iphot]>-1){
00705 
00706           for (unsigned iConv=0; iConv<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++){
00707 
00708             reco::CaloClusterPtrVector scPtrVec;
00709 
00710             math::Error<3>::type error;
00711             const reco::Vertex  * convVtx = new reco::Vertex(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition(), error);
00712 
00713             //cout << "Vtx x="<<convVtx->x() << " y="<< convVtx->y()<<" z="<<convVtx->z()<< endl;
00714             //cout << "VtxError x=" << convVtx->xError() << endl;
00715 
00716             std::vector<reco::TrackRef> OneLegConvVector;
00717             OneLegConvVector.push_back(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]);
00718         
00719             reco::Conversion myOneLegConversion(scPtrVec, OneLegConvVector, *convVtx, reco::Conversion::pflow);
00720         
00721 
00722             std::vector< float > OneLegMvaVector;
00723             OneLegMvaVector.push_back(pfSingleLegConvMva_[conv1legPFCandidateIndex_[iphot]][iConv]);
00724             myOneLegConversion.setOneLegMVA(OneLegMvaVector);
00725             
00726             /*
00727             std::cout << "One leg conversion created" << endl;
00728             std::vector<edm::RefToBase<reco::Track> > convtracks = myOneLegConversion.tracks();
00729             const std::vector<float> mvalist = myOneLegConversion.oneLegMVA();
00730 
00731             cout << "nTracks="<< convtracks.size()<<endl;
00732             for (unsigned int itk=0; itk<convtracks.size(); itk++){
00733               //double convtrackpt = convtracks[itk]->pt();
00734               std::cout << "Track pt="<< convtracks[itk]->pt() << std::endl;
00735               std::cout << "Track mva="<< mvalist[itk] << std::endl;
00736             }   
00737             */
00738             oneLegConversions.push_back(myOneLegConversion);
00739             
00740             //cout << "OneLegConv added"<<endl;
00741             
00742           }
00743         }
00744       }
00745 }
00746 
00747 
00748 void PFPhotonTranslator::createPhotonCores(const edm::OrphanHandle<reco::SuperClusterCollection> & superClustersHandle, const edm::OrphanHandle<reco::ConversionCollection> & oneLegConversionHandle, reco::PhotonCoreCollection &photonCores)
00749 {
00750   
00751   //std::cout << "createPhotonCores" << std::endl;
00752 
00753   unsigned nphot=photPFCandidateIndex_.size();
00754 
00755   unsigned i1legtot = 0;
00756 
00757   for(unsigned iphot=0;iphot<nphot;++iphot)
00758     {
00759       //std::cout << "iphot="<<iphot<<std::endl;
00760 
00761       reco::PhotonCore myPhotonCore;
00762 
00763       reco::SuperClusterRef SCref(reco::SuperClusterRef(superClustersHandle, iphot));
00764       
00765       myPhotonCore.setPFlowPhoton(true);
00766       myPhotonCore.setStandardPhoton(false);
00767       myPhotonCore.setPflowSuperCluster(SCref);
00768       myPhotonCore.setSuperCluster(egSCRef_[iphot]);
00769 
00770       reco::ElectronSeedRefVector pixelSeeds = egPhotonRef_[iphot]->electronPixelSeeds();
00771       for (unsigned iseed=0; iseed<pixelSeeds.size(); iseed++){
00772         myPhotonCore.addElectronPixelSeed(pixelSeeds[iseed]);
00773       }
00774 
00775       //cout << "PhotonCores : SC OK" << endl;
00776 
00777       //cout << "conv1legPFCandidateIndex_[iphot]="<<conv1legPFCandidateIndex_[iphot]<<endl;
00778       //cout << "conv2legPFCandidateIndex_[iphot]="<<conv2legPFCandidateIndex_[iphot]<<endl;
00779 
00780       if (conv1legPFCandidateIndex_[iphot]>-1){
00781         for (unsigned int iConv=0; iConv<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++){
00782           
00783           const reco::ConversionRef & OneLegRef(reco::ConversionRef(oneLegConversionHandle, i1legtot));
00784           myPhotonCore.addOneLegConversion(OneLegRef);
00785           
00786           //cout << "PhotonCores : 1-leg OK" << endl;
00787           /*
00788           cout << "Testing 1-leg :"<<endl;
00789           const reco::ConversionRefVector & conv = myPhotonCore.conversionsOneLeg();
00790           for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00791             cout << "Testing 1-leg : iconv="<<iconv<<endl;
00792             cout << "Testing 1-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
00793             cout << "Testing 1-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
00794             std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
00795             for (unsigned int itk=0; itk<convtracks.size(); itk++){
00796               //double convtrackpt = convtracks[itk]->pt();
00797               std::cout << "Testing 1-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
00798               // std::cout << "Track mva="<< mvalist[itk] << std::endl;
00799             }  
00800           }
00801           */
00802 
00803           i1legtot++;
00804         }
00805       }
00806 
00807       if (conv2legPFCandidateIndex_[iphot]>-1){
00808         for(unsigned int iConv=0; iConv<pfConv_[conv2legPFCandidateIndex_[iphot]].size(); iConv++) {
00809 
00810           const reco::ConversionRef & TwoLegRef(pfConv_[conv2legPFCandidateIndex_[iphot]][iConv]);
00811           myPhotonCore.addConversion(TwoLegRef);
00812 
00813         }
00814         //cout << "PhotonCores : 2-leg OK" << endl;
00815 
00816         /*
00817         cout << "Testing 2-leg :"<<endl;
00818         const reco::ConversionRefVector & conv = myPhotonCore.conversions();
00819         for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00820           cout << "Testing 2-leg : iconv="<<iconv<<endl;
00821           cout << "Testing 2-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
00822           cout << "Testing 2-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
00823           std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
00824           for (unsigned int itk=0; itk<convtracks.size(); itk++){
00825             //double convtrackpt = convtracks[itk]->pt();
00826             std::cout << "Testing 2-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
00827             // std::cout << "Track mva="<< mvalist[itk] << std::endl;
00828           }  
00829         }
00830         */
00831       }
00832 
00833       photonCores.push_back(myPhotonCore);
00834       
00835     }
00836 
00837   //std::cout << "end of createPhotonCores"<<std::endl;
00838 }
00839 
00840 
00841 void PFPhotonTranslator::createPhotons(reco::VertexCollection &vertexCollection, edm::Handle<reco::PhotonCollection> &egPhotons, const edm::OrphanHandle<reco::PhotonCoreCollection> & photonCoresHandle, const IsolationValueMaps& isolationValues, reco::PhotonCollection &photons) 
00842 {
00843 
00844   //cout << "createPhotons" << endl;
00845   
00846   unsigned nphot=photPFCandidateIndex_.size();
00847 
00848   for(unsigned iphot=0;iphot<nphot;++iphot)
00849     {
00850       //std::cout << "iphot="<<iphot<<std::endl;
00851 
00852       reco::PhotonCoreRef PCref(reco::PhotonCoreRef(photonCoresHandle, iphot));
00853 
00854       math::XYZPoint vtx(0.,0.,0.);
00855       if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
00856       //std::cout << "vtx made" << std::endl;
00857 
00858       math::XYZVector direction =  PCref->pfSuperCluster()->position() - vtx;
00859 
00860       //It could be that pfSC energy gives not the best resolution : use smaller agregates for some cases ?
00861       math::XYZVector P3 = direction.unit() * PCref->pfSuperCluster()->energy();
00862       LorentzVector P4(P3.x(), P3.y(), P3.z(), PCref->pfSuperCluster()->energy());
00863 
00864       reco::Photon myPhoton(P4, PCref->pfSuperCluster()->position(), PCref, vtx);
00865       //cout << "photon created"<<endl;
00866 
00867 
00868 
00869       reco::Photon::ShowerShape  showerShape;
00870       reco::Photon::FiducialFlags fiducialFlags;
00871       reco::Photon::IsolationVariables isolationVariables03;
00872       reco::Photon::IsolationVariables isolationVariables04;
00873 
00874       showerShape.e1x5= egPhotonRef_[iphot]->e1x5();
00875       showerShape.e2x5= egPhotonRef_[iphot]->e2x5();
00876       showerShape.e3x3= egPhotonRef_[iphot]->e3x3();
00877       showerShape.e5x5= egPhotonRef_[iphot]->e5x5();
00878       showerShape.maxEnergyXtal =  egPhotonRef_[iphot]->maxEnergyXtal();
00879       showerShape.sigmaEtaEta =    egPhotonRef_[iphot]->sigmaEtaEta();
00880       showerShape.sigmaIetaIeta =  egPhotonRef_[iphot]->sigmaIetaIeta();
00881       showerShape.hcalDepth1OverEcal = egPhotonRef_[iphot]->hadronicDepth1OverEm();
00882       showerShape.hcalDepth2OverEcal = egPhotonRef_[iphot]->hadronicDepth2OverEm();
00883       myPhoton.setShowerShapeVariables ( showerShape ); 
00884           
00885       fiducialFlags.isEB = egPhotonRef_[iphot]->isEB();
00886       fiducialFlags.isEE = egPhotonRef_[iphot]->isEE();
00887       fiducialFlags.isEBEtaGap = egPhotonRef_[iphot]->isEBEtaGap();
00888       fiducialFlags.isEBPhiGap = egPhotonRef_[iphot]->isEBPhiGap();
00889       fiducialFlags.isEERingGap = egPhotonRef_[iphot]->isEERingGap();
00890       fiducialFlags.isEEDeeGap = egPhotonRef_[iphot]->isEEDeeGap();
00891       fiducialFlags.isEBEEGap = egPhotonRef_[iphot]->isEBEEGap();
00892       myPhoton.setFiducialVolumeFlags ( fiducialFlags );
00893 
00894       isolationVariables03.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR03();
00895       isolationVariables03.hcalTowerSumEt = egPhotonRef_[iphot]->hcalTowerSumEtConeDR03();
00896       isolationVariables03.hcalDepth1TowerSumEt = egPhotonRef_[iphot]->hcalDepth1TowerSumEtConeDR03();
00897       isolationVariables03.hcalDepth2TowerSumEt = egPhotonRef_[iphot]->hcalDepth2TowerSumEtConeDR03();
00898       isolationVariables03.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR03();
00899       isolationVariables03.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR03();
00900       isolationVariables03.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR03();
00901       isolationVariables03.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR03();
00902       isolationVariables04.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR04();
00903       isolationVariables04.hcalTowerSumEt = egPhotonRef_[iphot]->hcalTowerSumEtConeDR04();
00904       isolationVariables04.hcalDepth1TowerSumEt = egPhotonRef_[iphot]->hcalDepth1TowerSumEtConeDR04();
00905       isolationVariables04.hcalDepth2TowerSumEt = egPhotonRef_[iphot]->hcalDepth2TowerSumEtConeDR04();
00906       isolationVariables04.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR04();
00907       isolationVariables04.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR04();
00908       isolationVariables04.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR04();
00909       isolationVariables04.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR04();
00910       myPhoton.setIsolationVariables(isolationVariables04, isolationVariables03);
00911      
00912           
00913       
00914 
00915       reco::Photon::PflowIsolationVariables myPFIso;
00916       myPFIso.chargedHadronIso=(*isolationValues[0])[CandidatePtr_[iphot]];
00917       myPFIso.photonIso=(*isolationValues[1])[CandidatePtr_[iphot]];
00918       myPFIso.neutralHadronIso=(*isolationValues[2])[CandidatePtr_[iphot]];   
00919       myPhoton.setPflowIsolationVariables(myPFIso);
00920       
00921       //cout << "chargedHadronIso="<<myPhoton.chargedHadronIso()<<" photonIso="<<myPhoton.photonIso()<<" neutralHadronIso="<<myPhoton.neutralHadronIso()<<endl;
00922       
00923       /*
00924       if (basicClusters_[iphot].size()>0){
00925       // Cluster shape variables
00926       //Algorithms from EcalClusterTools could be adapted to PF photons ? (not based on 5x5 BC)
00927       //It happens that energy computed in eg e5x5 is greater than pfSC energy (EcalClusterTools gathering energies from adjacent crystals even if not belonging to the SC)
00928       const EcalRecHitCollection* hits = 0 ;
00929       int subdet = PCref->pfSuperCluster()->seed()->hitsAndFractions()[0].first.subdetId();
00930       if (subdet==EcalBarrel) hits = barrelRecHits;
00931       else if  (subdet==EcalEndcap) hits = endcapRecHits;
00932       const CaloGeometry* geometry = theCaloGeom_.product();
00933 
00934       float maxXtal =   EcalClusterTools::eMax( *(PCref->pfSuperCluster()->seed()), &(*hits) );
00935       //cout << "maxXtal="<<maxXtal<<endl;
00936       float e1x5    =   EcalClusterTools::e1x5(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00937       //cout << "e1x5="<<e1x5<<endl;
00938       float e2x5    =   EcalClusterTools::e2x5Max(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00939       //cout << "e2x5="<<e2x5<<endl;
00940       float e3x3    =   EcalClusterTools::e3x3(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00941       //cout << "e3x3="<<e3x3<<endl;
00942       float e5x5    =   EcalClusterTools::e5x5( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00943       //cout << "e5x5="<<e5x5<<endl;
00944       std::vector<float> cov =  EcalClusterTools::covariances( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology), geometry); 
00945       float sigmaEtaEta = sqrt(cov[0]);
00946       //cout << "sigmaEtaEta="<<sigmaEtaEta<<endl;
00947       std::vector<float> locCov =  EcalClusterTools::localCovariances( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00948       float sigmaIetaIeta = sqrt(locCov[0]);
00949       //cout << "sigmaIetaIeta="<<sigmaIetaIeta<<endl;
00950       //float r9 =e3x3/(PCref->pfSuperCluster()->rawEnergy());
00951 
00952 
00953       // calculate HoE
00954       const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
00955       EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;  
00956       EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;  
00957       double HoE1=towerIso1.getTowerESum(&(*PCref->pfSuperCluster()))/PCref->pfSuperCluster()->energy();
00958       double HoE2=towerIso2.getTowerESum(&(*PCref->pfSuperCluster()))/PCref->pfSuperCluster()->energy(); 
00959       //cout << "HoE1="<<HoE1<<endl;
00960       //cout << "HoE2="<<HoE2<<endl;  
00961 
00962       reco::Photon::ShowerShape  showerShape;
00963       showerShape.e1x5= e1x5;
00964       showerShape.e2x5= e2x5;
00965       showerShape.e3x3= e3x3;
00966       showerShape.e5x5= e5x5;
00967       showerShape.maxEnergyXtal =  maxXtal;
00968       showerShape.sigmaEtaEta =    sigmaEtaEta;
00969       showerShape.sigmaIetaIeta =  sigmaIetaIeta;
00970       showerShape.hcalDepth1OverEcal = HoE1;
00971       showerShape.hcalDepth2OverEcal = HoE2;
00972       myPhoton.setShowerShapeVariables ( showerShape ); 
00973       //cout << "shower shape variables filled"<<endl;
00974       }
00975       */
00976       
00977 
00978       photons.push_back(myPhoton);
00979 
00980     }
00981 
00982   //std::cout << "end of createPhotons"<<std::endl;
00983 }
00984 
00985 
00986 const reco::PFCandidate & PFPhotonTranslator::correspondingDaughterCandidate(const reco::PFCandidate & cand, const reco::PFBlockElement & pfbe) const
00987 {
00988   unsigned refindex=pfbe.index();
00989   //  std::cout << " N daughters " << cand.numberOfDaughters() << std::endl;
00990   reco::PFCandidate::const_iterator myDaughterCandidate=cand.begin();
00991   reco::PFCandidate::const_iterator itend=cand.end();
00992 
00993   for(;myDaughterCandidate!=itend;++myDaughterCandidate)
00994     {
00995       const reco::PFCandidate * myPFCandidate = (const reco::PFCandidate*)&*myDaughterCandidate;
00996       if(myPFCandidate->elementsInBlocks().size()!=1)
00997         {
00998           //      std::cout << " Daughter with " << myPFCandidate.elementsInBlocks().size()<< " element in block " << std::endl;
00999           return cand;
01000         }
01001       if(myPFCandidate->elementsInBlocks()[0].second==refindex) 
01002         {
01003           //      std::cout << " Found it " << cand << std::endl;
01004           return *myPFCandidate;
01005         }      
01006     }
01007   return cand;
01008 }
01009
CMSSW_4_4_3_patch1/src/RecoParticleFlow/PFProducer/plugins/PFPhotonTranslator.cc
