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