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.energy(),
00546                                             myPFCluster.position(),
00547                                             myPFCluster.caloID(),
00548                                             myPFCluster.hitsAndFractions(),
00549                                             myPFCluster.algo(),
00550                                             myPFCluster.seed()));
00551 }
00552 
00553 void PFPhotonTranslator::createPreshowerCluster(const reco::PFBlockElement & PFBE, reco::PreshowerClusterCollection& preshowerClusters,unsigned plane) const
00554 {
00555   reco::PFClusterRef  myPFClusterRef= PFBE.clusterRef();
00556   preshowerClusters.push_back(reco::PreshowerCluster(myPFClusterRef->energy(),myPFClusterRef->position(),
00557                                                myPFClusterRef->hitsAndFractions(),plane));
00558 }
00559 
00560 void PFPhotonTranslator::createBasicClusterPtrs(const edm::OrphanHandle<reco::BasicClusterCollection> & basicClustersHandle )
00561 {
00562   unsigned size=photPFCandidateIndex_.size();
00563   unsigned basicClusterCounter=0;
00564   basicClusterPtr_.resize(size);
00565 
00566   for(unsigned iphot=0;iphot<size;++iphot) // loop on tracks
00567     {
00568       unsigned nbc=basicClusters_[iphot].size();
00569       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00570         {
00571           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00572           reco::CaloClusterPtr bcPtr(basicClustersHandle,basicClusterCounter);
00573           basicClusterPtr_[iphot].push_back(bcPtr);
00574           ++basicClusterCounter;
00575         }
00576     }
00577 }
00578 
00579 void PFPhotonTranslator::createPreshowerClusterPtrs(const edm::OrphanHandle<reco::PreshowerClusterCollection> & preshowerClustersHandle )
00580 {
00581   unsigned size=photPFCandidateIndex_.size();
00582   unsigned psClusterCounter=0;
00583   preshowerClusterPtr_.resize(size);
00584 
00585   for(unsigned iphot=0;iphot<size;++iphot) // loop on tracks
00586     {
00587       unsigned nbc=preshowerClusters_[iphot].size();
00588       for(unsigned ibc=0;ibc<nbc;++ibc) // loop on basic clusters
00589         {
00590           //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
00591           reco::CaloClusterPtr psPtr(preshowerClustersHandle,psClusterCounter);
00592           preshowerClusterPtr_[iphot].push_back(psPtr);
00593           ++psClusterCounter;
00594         }
00595     }
00596 }
00597 
00598 void PFPhotonTranslator::createSuperClusters(const reco::PFCandidateCollection & pfCand,
00599                                                reco::SuperClusterCollection &superClusters) const
00600 {
00601   unsigned nphot=photPFCandidateIndex_.size();
00602   for(unsigned iphot=0;iphot<nphot;++iphot)
00603     {
00604 
00605       //cout << "SC iphot=" << iphot << endl;
00606 
00607       // Computes energy position a la e/gamma 
00608       double sclusterE=0;
00609       double posX=0.;
00610       double posY=0.;
00611       double posZ=0.;
00612       
00613       unsigned nbasics=basicClusters_[iphot].size();
00614       for(unsigned ibc=0;ibc<nbasics;++ibc)
00615         {
00616           //cout << "BC in SC : iphot="<<iphot<<endl;
00617           
00618           double e = basicClusters_[iphot][ibc].energy();
00619           sclusterE += e;
00620           posX += e * basicClusters_[iphot][ibc].position().X();
00621           posY += e * basicClusters_[iphot][ibc].position().Y();
00622           posZ += e * basicClusters_[iphot][ibc].position().Z();          
00623         }
00624       posX /=sclusterE;
00625       posY /=sclusterE;
00626       posZ /=sclusterE;
00627       
00628       /*
00629       if(pfCand[gsfPFCandidateIndex_[iphot]].gsfTrackRef()!=GsfTrackRef_[iphot])
00630         {
00631           edm::LogError("PFElectronTranslator") << " Major problem in PFElectron Translator" << std::endl;
00632         }
00633       */      
00634 
00635       // compute the width
00636       PFClusterWidthAlgo pfwidth(pfClusters_[iphot]);
00637       
00638       double correctedEnergy=pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
00639       reco::SuperCluster mySuperCluster(correctedEnergy,math::XYZPoint(posX,posY,posZ));
00640       // protection against empty basic cluster collection ; the value is -2 in this case
00641       if(nbasics)
00642         {
00643 //        std::cout << "SuperCluster creation; energy " << pfCand[gsfPFCandidateIndex_[iphot]].ecalEnergy();
00644 //        std::cout << " " <<   pfCand[gsfPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
00645 //        std::cout << "Seed energy from basic " << basicClusters_[iphot][0].energy() << std::endl;
00646           mySuperCluster.setSeed(basicClusterPtr_[iphot][0]);
00647         }
00648       else
00649         {
00650           //      std::cout << "SuperCluster creation ; seed energy " << 0 << std::endl;
00651           //std::cout << "SuperCluster creation ; energy " << pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
00652           //std::cout << " " <<   pfCand[photPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
00653 //        std::cout << " No seed found " << 0 << std::endl;       
00654 //        std::cout << " MVA " << pfCand[gsfPFCandidateIndex_[iphot]].mva_e_pi() << std::endl;
00655           mySuperCluster.setSeed(reco::CaloClusterPtr());
00656         }
00657       // the seed should be the first basic cluster
00658 
00659       for(unsigned ibc=0;ibc<nbasics;++ibc)
00660         {
00661           mySuperCluster.addCluster(basicClusterPtr_[iphot][ibc]);
00662           //      std::cout <<"Adding Ref to SC " << basicClusterPtr_[iphot][ibc].index() << std::endl;
00663           const std::vector< std::pair<DetId, float> > & v1 = basicClusters_[iphot][ibc].hitsAndFractions();
00664           //      std::cout << " Number of cells " << v1.size() << std::endl;
00665           for( std::vector< std::pair<DetId, float> >::const_iterator diIt = v1.begin();
00666                diIt != v1.end();
00667                ++diIt ) {
00668             //      std::cout << " Adding DetId " << (diIt->first).rawId() << " " << diIt->second << std::endl;
00669             mySuperCluster.addHitAndFraction(diIt->first,diIt->second);
00670           } // loop over rechits      
00671         }      
00672 
00673       unsigned nps=preshowerClusterPtr_[iphot].size();
00674       for(unsigned ips=0;ips<nps;++ips)
00675         {
00676           mySuperCluster.addPreshowerCluster(preshowerClusterPtr_[iphot][ips]);
00677         }
00678       
00679 
00680       // Set the preshower energy
00681       mySuperCluster.setPreshowerEnergy(pfCand[photPFCandidateIndex_[iphot]].pS1Energy()+
00682                                         pfCand[photPFCandidateIndex_[iphot]].pS2Energy());
00683 
00684       // Set the cluster width
00685       mySuperCluster.setEtaWidth(pfwidth.pflowEtaWidth());
00686       mySuperCluster.setPhiWidth(pfwidth.pflowPhiWidth());
00687       // Force the computation of rawEnergy_ of the reco::SuperCluster
00688       mySuperCluster.rawEnergy();
00689 
00690       //cout << "SC energy="<< mySuperCluster.energy()<<endl;
00691 
00692       superClusters.push_back(mySuperCluster);
00693       //std::cout << "nb super clusters in collection : "<<superClusters.size()<<std::endl;
00694     }
00695 }
00696 
00697 void PFPhotonTranslator::createOneLegConversions(const edm::OrphanHandle<reco::SuperClusterCollection> & superClustersHandle, reco::ConversionCollection &oneLegConversions)
00698 {
00699 
00700   //std::cout << "createOneLegConversions" << std::endl;
00701 
00702     unsigned nphot=photPFCandidateIndex_.size();
00703     for(unsigned iphot=0;iphot<nphot;++iphot)
00704       {
00705 
00706         //if (conv1legPFCandidateIndex_[iphot]==-1) cout << "No OneLegConversions to add"<<endl;
00707         //else std::cout << "Phot "<<iphot<< " nOneLegConversions to add : "<<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size()<<endl;
00708 
00709 
00710         if (conv1legPFCandidateIndex_[iphot]>-1){
00711 
00712           for (unsigned iConv=0; iConv<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++){
00713 
00714             reco::CaloClusterPtrVector scPtrVec;
00715             std::vector<reco::CaloClusterPtr>matchingBC;
00716             math::Error<3>::type error;
00717             const reco::Vertex  * convVtx = new reco::Vertex(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition(), error);
00718             
00719             //cout << "Vtx x="<<convVtx->x() << " y="<< convVtx->y()<<" z="<<convVtx->z()<< endl;
00720             //cout << "VtxError x=" << convVtx->xError() << endl;
00721 
00722             std::vector<reco::TrackRef> OneLegConvVector;
00723             OneLegConvVector.push_back(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]);
00724             
00725             reco::CaloClusterPtrVector clu=scPtrVec;
00726             std::vector<reco::TrackRef> tr=OneLegConvVector;
00727             std::vector<math::XYZPointF>trackPositionAtEcalVec;
00728             std::vector<math::XYZPointF>innPointVec;
00729             std::vector<math::XYZVectorF>trackPinVec;
00730             std::vector<math::XYZVectorF>trackPoutVec;
00731             math::XYZPointF trackPositionAtEcal(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00732                                                 outerPosition().X(), 
00733                                                 pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00734                                                 outerPosition().Y(),
00735                                                 pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00736                                                 outerPosition().Z());
00737             math::XYZPointF innPoint(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00738                                      innerPosition().X(), 
00739                                      pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00740                                      innerPosition().Y(),
00741                                      pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00742                                      innerPosition().Z());
00743             math::XYZVectorF trackPin(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00744                                      innerMomentum().X(), 
00745                                      pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00746                                      innerMomentum().Y(),
00747                                      pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00748                                      innerMomentum().Z());
00749             math::XYZVectorF trackPout(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00750                                       outerMomentum().X(), 
00751                                       pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00752                                       outerMomentum().Y(),
00753                                       pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->
00754                                       outerMomentum().Z());
00755             float DCA=pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->d0();
00756             trackPositionAtEcalVec.push_back(trackPositionAtEcal);
00757             innPointVec.push_back(innPoint);
00758             trackPinVec.push_back(trackPin);
00759             trackPoutVec.push_back(trackPout);
00760             std::vector< float > OneLegMvaVector;
00761             reco::Conversion myOneLegConversion(scPtrVec, 
00762                                                 OneLegConvVector,
00763                                                 trackPositionAtEcalVec,
00764                                                 *convVtx,
00765                                                 matchingBC,
00766                                                 DCA,
00767                                                 innPointVec,
00768                                                 trackPinVec,
00769                                                 trackPoutVec,
00770                                                 pfSingleLegConvMva_[conv1legPFCandidateIndex_[iphot]][iConv],                     
00771                                                 reco::Conversion::pflow);
00772             OneLegMvaVector.push_back(pfSingleLegConvMva_[conv1legPFCandidateIndex_[iphot]][iConv]);
00773             myOneLegConversion.setOneLegMVA(OneLegMvaVector);
00774             //reco::Conversion myOneLegConversion(scPtrVec, 
00775             //OneLegConvVector, *convVtx, reco::Conversion::pflow);
00776             
00777             /*
00778             std::cout << "One leg conversion created" << endl;
00779             std::vector<edm::RefToBase<reco::Track> > convtracks = myOneLegConversion.tracks();
00780             const std::vector<float> mvalist = myOneLegConversion.oneLegMVA();
00781 
00782             cout << "nTracks="<< convtracks.size()<<endl;
00783             for (unsigned int itk=0; itk<convtracks.size(); itk++){
00784               //double convtrackpt = convtracks[itk]->pt();
00785               std::cout << "Track pt="<< convtracks[itk]->pt() << std::endl;
00786               std::cout << "Track mva="<< mvalist[itk] << std::endl;
00787             }   
00788             */
00789             oneLegConversions.push_back(myOneLegConversion);
00790             
00791             //cout << "OneLegConv added"<<endl;
00792             
00793           }
00794         }
00795       }
00796 }
00797 
00798 
00799 void PFPhotonTranslator::createPhotonCores(const edm::OrphanHandle<reco::SuperClusterCollection> & superClustersHandle, const edm::OrphanHandle<reco::ConversionCollection> & oneLegConversionHandle, reco::PhotonCoreCollection &photonCores)
00800 {
00801   
00802   //std::cout << "createPhotonCores" << std::endl;
00803 
00804   unsigned nphot=photPFCandidateIndex_.size();
00805 
00806   unsigned i1legtot = 0;
00807 
00808   for(unsigned iphot=0;iphot<nphot;++iphot)
00809     {
00810       //std::cout << "iphot="<<iphot<<std::endl;
00811 
00812       reco::PhotonCore myPhotonCore;
00813 
00814       reco::SuperClusterRef SCref(reco::SuperClusterRef(superClustersHandle, iphot));
00815       
00816       myPhotonCore.setPFlowPhoton(true);
00817       myPhotonCore.setStandardPhoton(false);
00818       myPhotonCore.setPflowSuperCluster(SCref);
00819       myPhotonCore.setSuperCluster(egSCRef_[iphot]);
00820 
00821       reco::ElectronSeedRefVector pixelSeeds = egPhotonRef_[iphot]->electronPixelSeeds();
00822       for (unsigned iseed=0; iseed<pixelSeeds.size(); iseed++){
00823         myPhotonCore.addElectronPixelSeed(pixelSeeds[iseed]);
00824       }
00825 
00826       //cout << "PhotonCores : SC OK" << endl;
00827 
00828       //cout << "conv1legPFCandidateIndex_[iphot]="<<conv1legPFCandidateIndex_[iphot]<<endl;
00829       //cout << "conv2legPFCandidateIndex_[iphot]="<<conv2legPFCandidateIndex_[iphot]<<endl;
00830 
00831       if (conv1legPFCandidateIndex_[iphot]>-1){
00832         for (unsigned int iConv=0; iConv<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++){
00833           
00834           const reco::ConversionRef & OneLegRef(reco::ConversionRef(oneLegConversionHandle, i1legtot));
00835           myPhotonCore.addOneLegConversion(OneLegRef);
00836           
00837           //cout << "PhotonCores : 1-leg OK" << endl;
00838           /*
00839           cout << "Testing 1-leg :"<<endl;
00840           const reco::ConversionRefVector & conv = myPhotonCore.conversionsOneLeg();
00841           for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00842             cout << "Testing 1-leg : iconv="<<iconv<<endl;
00843             cout << "Testing 1-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
00844             cout << "Testing 1-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
00845             std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
00846             for (unsigned int itk=0; itk<convtracks.size(); itk++){
00847               //double convtrackpt = convtracks[itk]->pt();
00848               std::cout << "Testing 1-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
00849               // std::cout << "Track mva="<< mvalist[itk] << std::endl;
00850             }  
00851           }
00852           */
00853 
00854           i1legtot++;
00855         }
00856       }
00857 
00858       if (conv2legPFCandidateIndex_[iphot]>-1){
00859         for(unsigned int iConv=0; iConv<pfConv_[conv2legPFCandidateIndex_[iphot]].size(); iConv++) {
00860 
00861           const reco::ConversionRef & TwoLegRef(pfConv_[conv2legPFCandidateIndex_[iphot]][iConv]);
00862           myPhotonCore.addConversion(TwoLegRef);
00863 
00864         }
00865         //cout << "PhotonCores : 2-leg OK" << endl;
00866 
00867         /*
00868         cout << "Testing 2-leg :"<<endl;
00869         const reco::ConversionRefVector & conv = myPhotonCore.conversions();
00870         for (unsigned int iconv=0; iconv<conv.size(); iconv++){
00871           cout << "Testing 2-leg : iconv="<<iconv<<endl;
00872           cout << "Testing 2-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
00873           cout << "Testing 2-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
00874           std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
00875           for (unsigned int itk=0; itk<convtracks.size(); itk++){
00876             //double convtrackpt = convtracks[itk]->pt();
00877             std::cout << "Testing 2-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
00878             // std::cout << "Track mva="<< mvalist[itk] << std::endl;
00879           }  
00880         }
00881         */
00882       }
00883 
00884       photonCores.push_back(myPhotonCore);
00885       
00886     }
00887 
00888   //std::cout << "end of createPhotonCores"<<std::endl;
00889 }
00890 
00891 
00892 void PFPhotonTranslator::createPhotons(reco::VertexCollection &vertexCollection, edm::Handle<reco::PhotonCollection> &egPhotons, const edm::OrphanHandle<reco::PhotonCoreCollection> & photonCoresHandle, const IsolationValueMaps& isolationValues, reco::PhotonCollection &photons) 
00893 {
00894 
00895   //cout << "createPhotons" << endl;
00896   
00897   unsigned nphot=photPFCandidateIndex_.size();
00898 
00899   for(unsigned iphot=0;iphot<nphot;++iphot)
00900     {
00901       //std::cout << "iphot="<<iphot<<std::endl;
00902 
00903       reco::PhotonCoreRef PCref(reco::PhotonCoreRef(photonCoresHandle, iphot));
00904 
00905       math::XYZPoint vtx(0.,0.,0.);
00906       if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
00907       //std::cout << "vtx made" << std::endl;
00908 
00909       math::XYZVector direction =  PCref->pfSuperCluster()->position() - vtx;
00910 
00911       //It could be that pfSC energy gives not the best resolution : use smaller agregates for some cases ?
00912       math::XYZVector P3 = direction.unit() * PCref->pfSuperCluster()->energy();
00913       LorentzVector P4(P3.x(), P3.y(), P3.z(), PCref->pfSuperCluster()->energy());
00914 
00915       reco::Photon myPhoton(P4, PCref->pfSuperCluster()->position(), PCref, vtx);
00916       //cout << "photon created"<<endl;
00917 
00918 
00919 
00920       reco::Photon::ShowerShape  showerShape;
00921       reco::Photon::FiducialFlags fiducialFlags;
00922       reco::Photon::IsolationVariables isolationVariables03;
00923       reco::Photon::IsolationVariables isolationVariables04;
00924 
00925       showerShape.e1x5= egPhotonRef_[iphot]->e1x5();
00926       showerShape.e2x5= egPhotonRef_[iphot]->e2x5();
00927       showerShape.e3x3= egPhotonRef_[iphot]->e3x3();
00928       showerShape.e5x5= egPhotonRef_[iphot]->e5x5();
00929       showerShape.maxEnergyXtal =  egPhotonRef_[iphot]->maxEnergyXtal();
00930       showerShape.sigmaEtaEta =    egPhotonRef_[iphot]->sigmaEtaEta();
00931       showerShape.sigmaIetaIeta =  egPhotonRef_[iphot]->sigmaIetaIeta();
00932       showerShape.hcalDepth1OverEcal = egPhotonRef_[iphot]->hadronicDepth1OverEm();
00933       showerShape.hcalDepth2OverEcal = egPhotonRef_[iphot]->hadronicDepth2OverEm();
00934       myPhoton.setShowerShapeVariables ( showerShape ); 
00935           
00936       fiducialFlags.isEB = egPhotonRef_[iphot]->isEB();
00937       fiducialFlags.isEE = egPhotonRef_[iphot]->isEE();
00938       fiducialFlags.isEBEtaGap = egPhotonRef_[iphot]->isEBEtaGap();
00939       fiducialFlags.isEBPhiGap = egPhotonRef_[iphot]->isEBPhiGap();
00940       fiducialFlags.isEERingGap = egPhotonRef_[iphot]->isEERingGap();
00941       fiducialFlags.isEEDeeGap = egPhotonRef_[iphot]->isEEDeeGap();
00942       fiducialFlags.isEBEEGap = egPhotonRef_[iphot]->isEBEEGap();
00943       myPhoton.setFiducialVolumeFlags ( fiducialFlags );
00944 
00945       isolationVariables03.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR03();
00946       isolationVariables03.hcalTowerSumEt = egPhotonRef_[iphot]->hcalTowerSumEtConeDR03();
00947       isolationVariables03.hcalDepth1TowerSumEt = egPhotonRef_[iphot]->hcalDepth1TowerSumEtConeDR03();
00948       isolationVariables03.hcalDepth2TowerSumEt = egPhotonRef_[iphot]->hcalDepth2TowerSumEtConeDR03();
00949       isolationVariables03.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR03();
00950       isolationVariables03.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR03();
00951       isolationVariables03.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR03();
00952       isolationVariables03.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR03();
00953       isolationVariables04.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR04();
00954       isolationVariables04.hcalTowerSumEt = egPhotonRef_[iphot]->hcalTowerSumEtConeDR04();
00955       isolationVariables04.hcalDepth1TowerSumEt = egPhotonRef_[iphot]->hcalDepth1TowerSumEtConeDR04();
00956       isolationVariables04.hcalDepth2TowerSumEt = egPhotonRef_[iphot]->hcalDepth2TowerSumEtConeDR04();
00957       isolationVariables04.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR04();
00958       isolationVariables04.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR04();
00959       isolationVariables04.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR04();
00960       isolationVariables04.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR04();
00961       myPhoton.setIsolationVariables(isolationVariables04, isolationVariables03);
00962      
00963           
00964       
00965 
00966       reco::Photon::PflowIsolationVariables myPFIso;
00967       myPFIso.chargedHadronIso=(*isolationValues[0])[CandidatePtr_[iphot]];
00968       myPFIso.photonIso=(*isolationValues[1])[CandidatePtr_[iphot]];
00969       myPFIso.neutralHadronIso=(*isolationValues[2])[CandidatePtr_[iphot]];   
00970       myPhoton.setPflowIsolationVariables(myPFIso);
00971       
00972       reco::Photon::PflowIDVariables myPFVariables;
00973 
00974       reco::Mustache myMustache;
00975       myMustache.MustacheID(*(myPhoton.pfSuperCluster()), myPFVariables.nClusterOutsideMustache, myPFVariables.etOutsideMustache );
00976       myPFVariables.mva = pfPhotonMva_[iphot];
00977       myPhoton.setPflowIDVariables(myPFVariables);
00978 
00979       //cout << "chargedHadronIso="<<myPhoton.chargedHadronIso()<<" photonIso="<<myPhoton.photonIso()<<" neutralHadronIso="<<myPhoton.neutralHadronIso()<<endl;
00980       
00981       // set PF-regression energy
00982       myPhoton.setCorrectedEnergy(reco::Photon::regression2,energyRegression_[iphot],energyRegressionError_[iphot],false);
00983       
00984 
00985       /*
00986       if (basicClusters_[iphot].size()>0){
00987       // Cluster shape variables
00988       //Algorithms from EcalClusterTools could be adapted to PF photons ? (not based on 5x5 BC)
00989       //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)
00990       const EcalRecHitCollection* hits = 0 ;
00991       int subdet = PCref->pfSuperCluster()->seed()->hitsAndFractions()[0].first.subdetId();
00992       if (subdet==EcalBarrel) hits = barrelRecHits;
00993       else if  (subdet==EcalEndcap) hits = endcapRecHits;
00994       const CaloGeometry* geometry = theCaloGeom_.product();
00995 
00996       float maxXtal =   EcalClusterTools::eMax( *(PCref->pfSuperCluster()->seed()), &(*hits) );
00997       //cout << "maxXtal="<<maxXtal<<endl;
00998       float e1x5    =   EcalClusterTools::e1x5(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
00999       //cout << "e1x5="<<e1x5<<endl;
01000       float e2x5    =   EcalClusterTools::e2x5Max(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
01001       //cout << "e2x5="<<e2x5<<endl;
01002       float e3x3    =   EcalClusterTools::e3x3(  *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
01003       //cout << "e3x3="<<e3x3<<endl;
01004       float e5x5    =   EcalClusterTools::e5x5( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
01005       //cout << "e5x5="<<e5x5<<endl;
01006       std::vector<float> cov =  EcalClusterTools::covariances( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology), geometry); 
01007       float sigmaEtaEta = sqrt(cov[0]);
01008       //cout << "sigmaEtaEta="<<sigmaEtaEta<<endl;
01009       std::vector<float> locCov =  EcalClusterTools::localCovariances( *(PCref->pfSuperCluster()->seed()), &(*hits), &(*topology)); 
01010       float sigmaIetaIeta = sqrt(locCov[0]);
01011       //cout << "sigmaIetaIeta="<<sigmaIetaIeta<<endl;
01012       //float r9 =e3x3/(PCref->pfSuperCluster()->rawEnergy());
01013 
01014 
01015       // calculate HoE
01016       const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
01017       EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;  
01018       EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;  
01019       double HoE1=towerIso1.getTowerESum(&(*PCref->pfSuperCluster()))/PCref->pfSuperCluster()->energy();
01020       double HoE2=towerIso2.getTowerESum(&(*PCref->pfSuperCluster()))/PCref->pfSuperCluster()->energy(); 
01021       //cout << "HoE1="<<HoE1<<endl;
01022       //cout << "HoE2="<<HoE2<<endl;  
01023 
01024       reco::Photon::ShowerShape  showerShape;
01025       showerShape.e1x5= e1x5;
01026       showerShape.e2x5= e2x5;
01027       showerShape.e3x3= e3x3;
01028       showerShape.e5x5= e5x5;
01029       showerShape.maxEnergyXtal =  maxXtal;
01030       showerShape.sigmaEtaEta =    sigmaEtaEta;
01031       showerShape.sigmaIetaIeta =  sigmaIetaIeta;
01032       showerShape.hcalDepth1OverEcal = HoE1;
01033       showerShape.hcalDepth2OverEcal = HoE2;
01034       myPhoton.setShowerShapeVariables ( showerShape ); 
01035       //cout << "shower shape variables filled"<<endl;
01036       }
01037       */
01038       
01039 
01040       photons.push_back(myPhoton);
01041 
01042     }
01043 
01044   //std::cout << "end of createPhotons"<<std::endl;
01045 }
01046 
01047 
01048 const reco::PFCandidate & PFPhotonTranslator::correspondingDaughterCandidate(const reco::PFCandidate & cand, const reco::PFBlockElement & pfbe) const
01049 {
01050   unsigned refindex=pfbe.index();
01051   //  std::cout << " N daughters " << cand.numberOfDaughters() << std::endl;
01052   reco::PFCandidate::const_iterator myDaughterCandidate=cand.begin();
01053   reco::PFCandidate::const_iterator itend=cand.end();
01054 
01055   for(;myDaughterCandidate!=itend;++myDaughterCandidate)
01056     {
01057       const reco::PFCandidate * myPFCandidate = (const reco::PFCandidate*)&*myDaughterCandidate;
01058       if(myPFCandidate->elementsInBlocks().size()!=1)
01059         {
01060           //      std::cout << " Daughter with " << myPFCandidate.elementsInBlocks().size()<< " element in block " << std::endl;
01061           return cand;
01062         }
01063       if(myPFCandidate->elementsInBlocks()[0].second==refindex) 
01064         {
01065           //      std::cout << " Found it " << cand << std::endl;
01066           return *myPFCandidate;
01067         }      
01068     }
01069   return cand;
01070 }
01071
/data/refman/pasoursint/CMSSW_5_3_3/src/RecoParticleFlow/PFProducer/plugins/PFPhotonTranslator.cc
