PFProducer/plugins/PFPhotonTranslator.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h"
0006 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0007 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
0008 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0009 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
0010 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0011 #include "DataFormats/VertexReco/interface/Vertex.h"
0012 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
0013 #include "DataFormats/Math/interface/Vector3D.h"
0014 #include "DataFormats/Math/interface/LorentzVector.h"
0015 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
0016 #include "FWCore/Framework/interface/stream/EDProducer.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "DataFormats/Common/interface/ValueMap.h"
0019 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidatePhotonExtra.h"
0023 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0024 #include "CondFormats/EcalObjects/interface/EcalMustacheSCParameters.h"
0025 #include "CondFormats/DataRecord/interface/EcalMustacheSCParametersRcd.h"
0026
0027 class CaloGeometry;
0028 class CaloTopology;
0029 class DetId;
0030 namespace edm {
0031   class EventSetup;
0032 }  // namespace edm
0033
0034 class PFPhotonTranslator : public edm::stream::EDProducer<> {
0035 public:
0036   explicit PFPhotonTranslator(const edm::ParameterSet &);
0037   ~PFPhotonTranslator() override;
0038
0039   void produce(edm::Event &, const edm::EventSetup &) override;
0040
0041   typedef std::vector<edm::Handle<edm::ValueMap<double> > > IsolationValueMaps;
0042
0043 private:
0044   // to retrieve the collection from the event
0045   bool fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection> &c,
0046                                 const edm::InputTag &tag,
0047                                 const edm::Event &iEvent) const;
0048
0049   // makes a basic cluster from PFBlockElement and add it to the collection ; the corrected energy is taken
0050   // from the PFCandidate
0051   void createBasicCluster(const reco::PFBlockElement &,
0052                           reco::BasicClusterCollection &basicClusters,
0053                           std::vector<const reco::PFCluster *> &,
0054                           const reco::PFCandidate &coCandidate) const;
0055   // makes a preshower cluster from of PFBlockElement and add it to the collection
0056   void createPreshowerCluster(const reco::PFBlockElement &PFBE,
0057                               reco::PreshowerClusterCollection &preshowerClusters,
0058                               unsigned plane) const;
0059
0060   // create the basic cluster Ptr
0061   void createBasicClusterPtrs(const edm::OrphanHandle<reco::BasicClusterCollection> &basicClustersHandle);
0062
0063   // create the preshower cluster Refs
0064   void createPreshowerClusterPtrs(const edm::OrphanHandle<reco::PreshowerClusterCollection> &preshowerClustersHandle);
0065
0066   // make a super cluster from its ingredients and add it to the collection
0067   void createSuperClusters(const reco::PFCandidateCollection &, reco::SuperClusterCollection &superClusters) const;
0068
0069   void createOneLegConversions(const edm::OrphanHandle<reco::SuperClusterCollection> &superClustersHandle,
0070                                reco::ConversionCollection &oneLegConversions);
0071
0072   //create photon cores
0073   void createPhotonCores(const edm::OrphanHandle<reco::SuperClusterCollection> &superClustersHandle,
0074                          const edm::OrphanHandle<reco::ConversionCollection> &oneLegConversionHandle,
0075                          reco::PhotonCoreCollection &photonCores);
0076
0077   void createPhotons(reco::VertexCollection &vertexCollection,
0078                      edm::Handle<reco::PhotonCollection> &egPhotons,
0079                      const edm::OrphanHandle<reco::PhotonCoreCollection> &photonCoresHandle,
0080                      const IsolationValueMaps &isolationValues,
0081                      reco::PhotonCollection &photons);
0082
0083   const reco::PFCandidate &correspondingDaughterCandidate(const reco::PFCandidate &cand,
0084                                                           const reco::PFBlockElement &pfbe) const;
0085
0086   edm::InputTag inputTagPFCandidates_;
0087   std::vector<edm::InputTag> inputTagIsoVals_;
0088   std::string PFBasicClusterCollection_;
0089   std::string PFPreshowerClusterCollection_;
0090   std::string PFSuperClusterCollection_;
0091   std::string PFPhotonCoreCollection_;
0092   std::string PFPhotonCollection_;
0093   std::string PFConversionCollection_;
0094   std::string EGPhotonCollection_;
0095   std::string vertexProducer_;
0096   edm::InputTag barrelEcalHits_;
0097   edm::InputTag endcapEcalHits_;
0098   edm::InputTag hcalTowers_;
0099   double hOverEConeSize_;
0100
0101   // the collection of basic clusters associated to a photon
0102   std::vector<reco::BasicClusterCollection> basicClusters_;
0103   // the correcsponding PFCluster ref
0104   std::vector<std::vector<const reco::PFCluster *> > pfClusters_;
0105   // the collection of preshower clusters associated to a photon
0106   std::vector<reco::PreshowerClusterCollection> preshowerClusters_;
0107   // the super cluster collection (actually only one) associated to a photon
0108   std::vector<reco::SuperClusterCollection> superClusters_;
0109   // the references to the basic clusters associated to a photon
0110   std::vector<reco::CaloClusterPtrVector> basicClusterPtr_;
0111   // the references to the basic clusters associated to a photon
0112   std::vector<reco::CaloClusterPtrVector> preshowerClusterPtr_;
0113   // keep track of the index of the PF Candidate
0114   std::vector<int> photPFCandidateIndex_;
0115   // the list of candidatePtr
0116   std::vector<reco::CandidatePtr> CandidatePtr_;
0117   // the e/g SC associated
0118   std::vector<reco::SuperClusterRef> egSCRef_;
0119   // the e/g photon associated
0120   std::vector<reco::PhotonRef> egPhotonRef_;
0121   // the PF MVA and regression
0122   std::vector<float> pfPhotonMva_;
0123   std::vector<float> energyRegression_;
0124   std::vector<float> energyRegressionError_;
0125
0126   //Vector of vector of Conversions Refs
0127   std::vector<reco::ConversionRefVector> pfConv_;
0128   std::vector<std::vector<reco::TrackRef> > pfSingleLegConv_;
0129   std::vector<std::vector<float> > pfSingleLegConvMva_;
0130
0131   std::vector<int> conv1legPFCandidateIndex_;
0132   std::vector<int> conv2legPFCandidateIndex_;
0133
0134   edm::ESHandle<CaloTopology> theCaloTopo_;
0135   edm::ESHandle<CaloGeometry> theCaloGeom_;
0136
0137   // Mustache SC parameters
0138   edm::ESGetToken<EcalMustacheSCParameters, EcalMustacheSCParametersRcd> ecalMustacheSCParametersToken_;
0139   const EcalMustacheSCParameters *mustacheSCParams_;
0140
0141   bool emptyIsOk_;
0142 };
0143
0144 DEFINE_FWK_MODULE(PFPhotonTranslator);
0145
0146 using namespace edm;
0147 using namespace std;
0148 using namespace reco;
0149
0150 typedef math::XYZTLorentzVector LorentzVector;
0151 typedef math::XYZPoint Point;
0152 typedef math::XYZVector Vector;
0153
0154 PFPhotonTranslator::PFPhotonTranslator(const edm::ParameterSet &iConfig) {
0155   //std::cout << "PFPhotonTranslator" << std::endl;
0156
0157   inputTagPFCandidates_ = iConfig.getParameter<edm::InputTag>("PFCandidate");
0158
0159   edm::ParameterSet isoVals = iConfig.getParameter<edm::ParameterSet>("isolationValues");
0160   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfChargedHadrons"));
0161   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfPhotons"));
0162   inputTagIsoVals_.push_back(isoVals.getParameter<edm::InputTag>("pfNeutralHadrons"));
0163
0164   PFBasicClusterCollection_ = iConfig.getParameter<std::string>("PFBasicClusters");
0165   PFPreshowerClusterCollection_ = iConfig.getParameter<std::string>("PFPreshowerClusters");
0166   PFSuperClusterCollection_ = iConfig.getParameter<std::string>("PFSuperClusters");
0167   PFConversionCollection_ = iConfig.getParameter<std::string>("PFConversionCollection");
0168   PFPhotonCoreCollection_ = iConfig.getParameter<std::string>("PFPhotonCores");
0169   PFPhotonCollection_ = iConfig.getParameter<std::string>("PFPhotons");
0170
0171   EGPhotonCollection_ = iConfig.getParameter<std::string>("EGPhotons");
0172
0173   vertexProducer_ = iConfig.getParameter<std::string>("primaryVertexProducer");
0174
0175   barrelEcalHits_ = iConfig.getParameter<edm::InputTag>("barrelEcalHits");
0176   endcapEcalHits_ = iConfig.getParameter<edm::InputTag>("endcapEcalHits");
0177
0178   hcalTowers_ = iConfig.getParameter<edm::InputTag>("hcalTowers");
0179   hOverEConeSize_ = iConfig.getParameter<double>("hOverEConeSize");
0180
0181   if (iConfig.exists("emptyIsOk"))
0182     emptyIsOk_ = iConfig.getParameter<bool>("emptyIsOk");
0183   else
0184     emptyIsOk_ = false;
0185
0186   ecalMustacheSCParametersToken_ = esConsumes<EcalMustacheSCParameters, EcalMustacheSCParametersRcd>();
0187
0188   produces<reco::BasicClusterCollection>(PFBasicClusterCollection_);
0189   produces<reco::PreshowerClusterCollection>(PFPreshowerClusterCollection_);
0190   produces<reco::SuperClusterCollection>(PFSuperClusterCollection_);
0191   produces<reco::PhotonCoreCollection>(PFPhotonCoreCollection_);
0192   produces<reco::PhotonCollection>(PFPhotonCollection_);
0193   produces<reco::ConversionCollection>(PFConversionCollection_);
0194 }
0195
0196 PFPhotonTranslator::~PFPhotonTranslator() {}
0197
0198 void PFPhotonTranslator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0199   //cout << "NEW EVENT"<<endl;
0200   mustacheSCParams_ = &iSetup.getData(ecalMustacheSCParametersToken_);
0201
0202   auto basicClusters_p = std::make_unique<reco::BasicClusterCollection>();
0203
0204   auto psClusters_p = std::make_unique<reco::PreshowerClusterCollection>();
0205
0206   /*
0207   auto SingleLeg_p = std::make_unique<reco::ConversionCollection>();
0208   */
0209
0210   reco::SuperClusterCollection outputSuperClusterCollection;
0211   reco::ConversionCollection outputOneLegConversionCollection;
0212   reco::PhotonCoreCollection outputPhotonCoreCollection;
0213   reco::PhotonCollection outputPhotonCollection;
0214
0215   outputSuperClusterCollection.clear();
0216   outputOneLegConversionCollection.clear();
0217   outputPhotonCoreCollection.clear();
0218   outputPhotonCollection.clear();
0219
0220   edm::Handle<reco::PFCandidateCollection> pfCandidates;
0221   bool status = fetchCandidateCollection(pfCandidates, inputTagPFCandidates_, iEvent);
0222
0223   edm::Handle<reco::PhotonCollection> egPhotons;
0224   iEvent.getByLabel(EGPhotonCollection_, egPhotons);
0225
0226   Handle<reco::VertexCollection> vertexHandle;
0227
0228   IsolationValueMaps isolationValues(inputTagIsoVals_.size());
0229   for (size_t j = 0; j < inputTagIsoVals_.size(); ++j) {
0230     iEvent.getByLabel(inputTagIsoVals_[j], isolationValues[j]);
0231   }
0232
0233   // clear the vectors
0234   photPFCandidateIndex_.clear();
0235   basicClusters_.clear();
0236   pfClusters_.clear();
0237   preshowerClusters_.clear();
0238   superClusters_.clear();
0239   basicClusterPtr_.clear();
0240   preshowerClusterPtr_.clear();
0241   CandidatePtr_.clear();
0242   egSCRef_.clear();
0243   egPhotonRef_.clear();
0244   pfPhotonMva_.clear();
0245   energyRegression_.clear();
0246   energyRegressionError_.clear();
0247   pfConv_.clear();
0248   pfSingleLegConv_.clear();
0249   pfSingleLegConvMva_.clear();
0250   conv1legPFCandidateIndex_.clear();
0251   conv2legPFCandidateIndex_.clear();
0252
0253   // loop on the candidates
0254   //CC@@
0255   // we need first to create AND put the SuperCluster,
0256   // basic clusters and presh clusters collection
0257   // in order to get a working Handle
0258   unsigned ncand = (status) ? pfCandidates->size() : 0;
0259
0260   unsigned iphot = 0;
0261   unsigned iconv1leg = 0;
0262   unsigned iconv2leg = 0;
0263
0264   for (unsigned i = 0; i < ncand; ++i) {
0265     const reco::PFCandidate &cand = (*pfCandidates)[i];
0266     if (cand.particleId() != reco::PFCandidate::gamma)
0267       continue;
0268     //cout << "cand.mva_nothing_gamma()="<<cand. mva_nothing_gamma()<<endl;
0269     if (cand.mva_nothing_gamma() > 0.001)  //Found PFPhoton with PFPhoton Extras saved
0270     {
0271       //cout << "NEW PHOTON" << endl;
0272
0273       //std::cout << "nDoubleLegConv="<<cand.photonExtraRef()->conversionRef().size()<<std::endl;
0274
0275       if (!cand.photonExtraRef()->conversionRef().empty()) {
0276         pfConv_.push_back(reco::ConversionRefVector());
0277
0278         const reco::ConversionRefVector &doubleLegConvColl = cand.photonExtraRef()->conversionRef();
0279         for (unsigned int iconv = 0; iconv < doubleLegConvColl.size(); iconv++) {
0280           pfConv_[iconv2leg].push_back(doubleLegConvColl[iconv]);
0281         }
0282
0283         conv2legPFCandidateIndex_.push_back(iconv2leg);
0284         iconv2leg++;
0285       } else
0286         conv2legPFCandidateIndex_.push_back(-1);
0287
0288       const std::vector<reco::TrackRef> &singleLegConvColl = cand.photonExtraRef()->singleLegConvTrackRef();
0289       const std::vector<float> &singleLegConvCollMva = cand.photonExtraRef()->singleLegConvMva();
0290
0291       //std::cout << "nSingleLegConv=" <<singleLegConvColl.size() << std::endl;
0292
0293       if (!singleLegConvColl.empty()) {
0294         pfSingleLegConv_.push_back(std::vector<reco::TrackRef>());
0295         pfSingleLegConvMva_.push_back(std::vector<float>());
0296
0297         //cout << "nTracks="<< singleLegConvColl.size()<<endl;
0298         for (unsigned int itk = 0; itk < singleLegConvColl.size(); itk++) {
0299           //cout << "Track pt="<< singleLegConvColl[itk]->pt() <<endl;
0300
0301           pfSingleLegConv_[iconv1leg].push_back(singleLegConvColl[itk]);
0302           pfSingleLegConvMva_[iconv1leg].push_back(singleLegConvCollMva[itk]);
0303         }
0304
0305         conv1legPFCandidateIndex_.push_back(iconv1leg);
0306
0307         iconv1leg++;
0308       } else
0309         conv1legPFCandidateIndex_.push_back(-1);
0310     }
0311
0312     photPFCandidateIndex_.push_back(i);
0313     pfPhotonMva_.push_back(cand.mva_nothing_gamma());
0314     energyRegression_.push_back(cand.photonExtraRef()->MVAGlobalCorrE());
0315     energyRegressionError_.push_back(cand.photonExtraRef()->MVAGlobalCorrEError());
0316     basicClusters_.push_back(reco::BasicClusterCollection());
0317     pfClusters_.push_back(std::vector<const reco::PFCluster *>());
0318     preshowerClusters_.push_back(reco::PreshowerClusterCollection());
0319     superClusters_.push_back(reco::SuperClusterCollection());
0320
0321     reco::PFCandidatePtr ptrToPFPhoton(pfCandidates, i);
0322     CandidatePtr_.push_back(ptrToPFPhoton);
0323     egSCRef_.push_back(cand.superClusterRef());
0324     //std::cout << "PFPhoton cand " << iphot << std::endl;
0325
0326     int iegphot = 0;
0327     for (reco::PhotonCollection::const_iterator gamIter = egPhotons->begin(); gamIter != egPhotons->end(); ++gamIter) {
0328       if (cand.superClusterRef() == gamIter->superCluster()) {
0329         reco::PhotonRef PhotRef(reco::PhotonRef(egPhotons, iegphot));
0330         egPhotonRef_.push_back(PhotRef);
0331       }
0332       iegphot++;
0333     }
0334
0335     //std::cout << "Cand elements in blocks : " << cand.elementsInBlocks().size() << std::endl;
0336
0337     for (unsigned iele = 0; iele < cand.elementsInBlocks().size(); ++iele) {
0338       // first get the block
0339       reco::PFBlockRef blockRef = cand.elementsInBlocks()[iele].first;
0340       //
0341       unsigned elementIndex = cand.elementsInBlocks()[iele].second;
0342       // check it actually exists
0343       if (blockRef.isNull())
0344         continue;
0345
0346       // then get the elements of the block
0347       const edm::OwnVector<reco::PFBlockElement> &elements = (*blockRef).elements();
0348
0349       const reco::PFBlockElement &pfbe(elements[elementIndex]);
0350       // The first ECAL element should be the cluster associated to the GSF; defined as the seed
0351       if (pfbe.type() == reco::PFBlockElement::ECAL) {
0352         //std::cout << "BlockElement ECAL" << std::endl;
0353         // the Brem photons are saved as daughter PFCandidate; this
0354         // is convenient to access the corrected energy
0355         //    std::cout << " Found candidate "  << correspondingDaughterCandidate(coCandidate,pfbe) << " " << coCandidate << std::endl;
0356         createBasicCluster(pfbe, basicClusters_[iphot], pfClusters_[iphot], correspondingDaughterCandidate(cand, pfbe));
0357       }
0358       if (pfbe.type() == reco::PFBlockElement::PS1) {
0359         //std::cout << "BlockElement PS1" << std::endl;
0360         createPreshowerCluster(pfbe, preshowerClusters_[iphot], 1);
0361       }
0362       if (pfbe.type() == reco::PFBlockElement::PS2) {
0363         //std::cout << "BlockElement PS2" << std::endl;
0364         createPreshowerCluster(pfbe, preshowerClusters_[iphot], 2);
0365       }
0366
0367     }  // loop on the elements
0368
0369     // save the basic clusters
0370     basicClusters_p->insert(basicClusters_p->end(), basicClusters_[iphot].begin(), basicClusters_[iphot].end());
0371     // save the preshower clusters
0372     psClusters_p->insert(psClusters_p->end(), preshowerClusters_[iphot].begin(), preshowerClusters_[iphot].end());
0373
0374     ++iphot;
0375
0376   }  // loop on PFCandidates
0377
0378   //Save the basic clusters and get an handle as to be able to create valid Refs (thanks to Claude)
0379   //  std::cout << " Number of basic clusters " << basicClusters_p->size() << std::endl;
0380   const edm::OrphanHandle<reco::BasicClusterCollection> bcRefProd =
0381       iEvent.put(std::move(basicClusters_p), PFBasicClusterCollection_);
0382
0383   //preshower clusters
0384   const edm::OrphanHandle<reco::PreshowerClusterCollection> psRefProd =
0385       iEvent.put(std::move(psClusters_p), PFPreshowerClusterCollection_);
0386
0387   // now that the Basic clusters are in the event, the Ref can be created
0388   createBasicClusterPtrs(bcRefProd);
0389   // now that the preshower clusters are in the event, the Ref can be created
0390   createPreshowerClusterPtrs(psRefProd);
0391
0392   // and now the Super cluster can be created with valid references
0393   //if(status) createSuperClusters(*pfCandidates,*superClusters_p);
0394   if (status)
0395     createSuperClusters(*pfCandidates, outputSuperClusterCollection);
0396
0397   //std::cout << "nb superclusters in collection : "<<outputSuperClusterCollection.size()<<std::endl;
0398
0399   // Let's put the super clusters in the event
0400   auto superClusters_p = std::make_unique<reco::SuperClusterCollection>(outputSuperClusterCollection);
0401   const edm::OrphanHandle<reco::SuperClusterCollection> scRefProd =
0402       iEvent.put(std::move(superClusters_p), PFSuperClusterCollection_);
0403
0404   /*
0405   int ipho=0;
0406   for (reco::SuperClusterCollection::const_iterator gamIter = scRefProd->begin(); gamIter != scRefProd->end(); ++gamIter){
0407     std::cout << "SC i="<<ipho<<" energy="<<gamIter->energy()<<std::endl;
0408     ipho++;
0409   }
0410   */
0411
0412   //1-leg conversions
0413
0414   if (status)
0415     createOneLegConversions(scRefProd, outputOneLegConversionCollection);
0416
0417   auto SingleLeg_p = std::make_unique<reco::ConversionCollection>(outputOneLegConversionCollection);
0418   const edm::OrphanHandle<reco::ConversionCollection> ConvRefProd =
0419       iEvent.put(std::move(SingleLeg_p), PFConversionCollection_);
0420   /*
0421   int iconv = 0;
0422   for (reco::ConversionCollection::const_iterator convIter = ConvRefProd->begin(); convIter != ConvRefProd->end(); ++convIter){
0423
0424     std::cout << "OneLegConv i="<<iconv<<" nTracks="<<convIter->nTracks()<<" EoverP="<<convIter->EoverP() <<std::endl;
0425     std::vector<edm::RefToBase<reco::Track> > convtracks = convIter->tracks();
0426     for (unsigned int itk=0; itk<convtracks.size(); itk++){
0427       std::cout << "Track pt="<< convtracks[itk]->pt() << std::endl;
0428     }
0429
0430     iconv++;
0431   }
0432   */
0433
0434   //create photon cores
0435   //if(status) createPhotonCores(pfCandidates, scRefProd, *photonCores_p);
0436   if (status)
0437     createPhotonCores(scRefProd, ConvRefProd, outputPhotonCoreCollection);
0438
0439   //std::cout << "nb photoncores in collection : "<<outputPhotonCoreCollection.size()<<std::endl;
0440
0441   // Put the photon cores in the event
0442   auto photonCores_p = std::make_unique<reco::PhotonCoreCollection>(outputPhotonCoreCollection);
0443   //std::cout << "photon core collection put in unique_ptr"<<std::endl;
0444   const edm::OrphanHandle<reco::PhotonCoreCollection> pcRefProd =
0445       iEvent.put(std::move(photonCores_p), PFPhotonCoreCollection_);
0446
0447   //std::cout << "photon core have been put in the event"<<std::endl;
0448   /*
0449   int ipho=0;
0450   for (reco::PhotonCoreCollection::const_iterator gamIter = pcRefProd->begin(); gamIter != pcRefProd->end(); ++gamIter){
0451     std::cout << "PhotonCore i="<<ipho<<" energy="<<gamIter->parentSuperCluster()->energy()<<std::endl;
0452     //for (unsigned int i=0; i<)
0453
0454     std::cout << "PhotonCore i="<<ipho<<" nconv2leg="<<gamIter->conversions().size()<<" nconv1leg="<<gamIter->conversionsOneLeg().size()<<std::endl;
0455
0456     const reco::ConversionRefVector & conv = gamIter->conversions();
0457     for (unsigned int iconv=0; iconv<conv.size(); iconv++){
0458       cout << "2-leg iconv="<<iconv<<endl;
0459       cout << "2-leg nTracks="<<conv[iconv]->nTracks()<<endl;
0460       cout << "2-leg EoverP="<<conv[iconv]->EoverP()<<endl;
0461       cout << "2-leg ConvAlgorithm="<<conv[iconv]->algo()<<endl;
0462     }
0463
0464     const reco::ConversionRefVector & convbis = gamIter->conversionsOneLeg();
0465     for (unsigned int iconv=0; iconv<convbis.size(); iconv++){
0466       cout << "1-leg iconv="<<iconv<<endl;
0467       cout << "1-leg nTracks="<<convbis[iconv]->nTracks()<<endl;
0468       cout << "1-leg EoverP="<<convbis[iconv]->EoverP()<<endl;
0469       cout << "1-leg ConvAlgorithm="<<convbis[iconv]->algo()<<endl;
0470     }
0471
0472     ipho++;
0473   }
0474   */
0475
0476   //load vertices
0477   reco::VertexCollection vertexCollection;
0478   bool validVertex = true;
0479   iEvent.getByLabel(vertexProducer_, vertexHandle);
0480   if (!vertexHandle.isValid()) {
0481     edm::LogWarning("PhotonProducer") << "Error! Can't get the product primary Vertex Collection "
0482                                       << "\n";
0483     validVertex = false;
0484   }
0485   if (validVertex)
0486     vertexCollection = *(vertexHandle.product());
0487
0488   /*
0489   //load Ecal rechits
0490   bool validEcalRecHits=true;
0491   Handle<EcalRecHitCollection> barrelHitHandle;
0492   EcalRecHitCollection barrelRecHits;
0493   iEvent.getByLabel(barrelEcalHits_, barrelHitHandle);
0494   if (!barrelHitHandle.isValid()) {
0495     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
0496     validEcalRecHits=false;
0497   }
0498   if (  validEcalRecHits)  barrelRecHits = *(barrelHitHandle.product());
0499
0500   Handle<EcalRecHitCollection> endcapHitHandle;
0501   iEvent.getByLabel(endcapEcalHits_, endcapHitHandle);
0502   EcalRecHitCollection endcapRecHits;
0503   if (!endcapHitHandle.isValid()) {
0504     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
0505     validEcalRecHits=false;
0506   }
0507   if( validEcalRecHits) endcapRecHits = *(endcapHitHandle.product());
0508
0509   //load detector topology & geometry
0510   iSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
0511
0512   edm::ESHandle<CaloTopology> pTopology;
0513   iSetup.get<CaloTopologyRecord>().get(theCaloTopo_);
0514   const CaloTopology *topology = theCaloTopo_.product();
0515
0516   // get Hcal rechits collection
0517   Handle<CaloTowerCollection> hcalTowersHandle;
0518   iEvent.getByLabel(hcalTowers_, hcalTowersHandle);
0519   */
0520
0521   //create photon collection
0522   //if(status) createPhotons(vertexCollection, pcRefProd, topology, &barrelRecHits, &endcapRecHits, hcalRecHitsHandle, isolationValues, outputPhotonCollection);
0523   if (status)
0524     createPhotons(vertexCollection, egPhotons, pcRefProd, isolationValues, outputPhotonCollection);
0525
0526   // Put the photons in the event
0527   auto photons_p = std::make_unique<reco::PhotonCollection>(outputPhotonCollection);
0528   //std::cout << "photon collection put in unique_ptr"<<std::endl;
0529   const edm::OrphanHandle<reco::PhotonCollection> photonRefProd = iEvent.put(std::move(photons_p), PFPhotonCollection_);
0530   //std::cout << "photons have been put in the event"<<std::endl;
0531
0532   /*
0533   ipho=0;
0534   for (reco::PhotonCollection::const_iterator gamIter = photonRefProd->begin(); gamIter != photonRefProd->end(); ++gamIter){
0535     std::cout << "Photon i="<<ipho<<" pfEnergy="<<gamIter->parentSuperCluster()->energy()<<std::endl;
0536
0537     const reco::ConversionRefVector & conv = gamIter->conversions();
0538     cout << "conversions obtained : conv.size()="<< conv.size()<<endl;
0539     for (unsigned int iconv=0; iconv<conv.size(); iconv++){
0540       cout << "iconv="<<iconv<<endl;
0541       cout << "nTracks="<<conv[iconv]->nTracks()<<endl;
0542       cout << "EoverP="<<conv[iconv]->EoverP()<<endl;
0543
0544       cout << "Vtx x="<<conv[iconv]->conversionVertex().x() << " y="<< conv[iconv]->conversionVertex().y()<<" z="<<conv[iconv]->conversionVertex().z()<< endl;
0545       cout << "VtxError x=" << conv[iconv]->conversionVertex().xError() << endl;
0546
0547       std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
0548       //cout << "nTracks="<< convtracks.size()<<endl;
0549       for (unsigned int itk=0; itk<convtracks.size(); itk++){
0550     double convtrackpt = convtracks[itk]->pt();
0551     const edm::RefToBase<reco::Track> & mytrack = convtracks[itk];
0552     cout << "Track pt="<< convtrackpt <<endl;
0553     cout << "Track origin="<<gamIter->conversionTrackOrigin(mytrack)<<endl;
0554       }
0555     }
0556
0557     //1-leg
0558     const reco::ConversionRefVector & convbis = gamIter->conversionsOneLeg();
0559     cout << "conversions obtained : conv.size()="<< convbis.size()<<endl;
0560     for (unsigned int iconv=0; iconv<convbis.size(); iconv++){
0561       cout << "iconv="<<iconv<<endl;
0562       cout << "nTracks="<<convbis[iconv]->nTracks()<<endl;
0563       cout << "EoverP="<<convbis[iconv]->EoverP()<<endl;
0564
0565       cout << "Vtx x="<<convbis[iconv]->conversionVertex().x() << " y="<< convbis[iconv]->conversionVertex().y()<<" z="<<convbis[iconv]->conversionVertex().z()<< endl;
0566       cout << "VtxError x=" << convbis[iconv]->conversionVertex().xError() << endl;
0567
0568       std::vector<edm::RefToBase<reco::Track> > convtracks = convbis[iconv]->tracks();
0569       //cout << "nTracks="<< convtracks.size()<<endl;
0570       for (unsigned int itk=0; itk<convtracks.size(); itk++){
0571     double convtrackpt = convtracks[itk]->pt();
0572     cout << "Track pt="<< convtrackpt <<endl;
0573     cout << "Track origin="<<gamIter->conversionTrackOrigin((convtracks[itk]))<<endl;
0574       }
0575     }
0576     ipho++;
0577   }
0578   */
0579 }
0580
0581 bool PFPhotonTranslator::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection> &c,
0582                                                   const edm::InputTag &tag,
0583                                                   const edm::Event &iEvent) const {
0584   bool found = iEvent.getByLabel(tag, c);
0585
0586   if (!found && !emptyIsOk_) {
0587     std::ostringstream err;
0588     err << " cannot get PFCandidates: " << tag << std::endl;
0589     edm::LogError("PFPhotonTranslator") << err.str();
0590   }
0591   return found;
0592 }
0593
0594 // The basic cluster is a copy of the PFCluster -> the energy is not corrected
0595 // It should be possible to get the corrected energy (including the associated PS energy)
0596 // from the PFCandidate daugthers ; Needs some work
0597 void PFPhotonTranslator::createBasicCluster(const reco::PFBlockElement &PFBE,
0598                                             reco::BasicClusterCollection &basicClusters,
0599                                             std::vector<const reco::PFCluster *> &pfClusters,
0600                                             const reco::PFCandidate &coCandidate) const {
0601   const reco::PFClusterRef &myPFClusterRef = PFBE.clusterRef();
0602   if (myPFClusterRef.isNull())
0603     return;
0604
0605   const reco::PFCluster &myPFCluster(*myPFClusterRef);
0606   pfClusters.push_back(&myPFCluster);
0607   //std::cout << " Creating BC " << myPFCluster.energy() << " " << coCandidate.ecalEnergy() <<" "<<  coCandidate.rawEcalEnergy() <<std::endl;
0608   //std::cout << " # hits " << myPFCluster.hitsAndFractions().size() << std::endl;
0609
0610   //  basicClusters.push_back(reco::CaloCluster(myPFCluster.energy(),
0611   basicClusters.push_back(reco::CaloCluster(  //coCandidate.rawEcalEnergy(),
0612       myPFCluster.energy(),
0613       myPFCluster.position(),
0614       myPFCluster.caloID(),
0615       myPFCluster.hitsAndFractions(),
0616       myPFCluster.algo(),
0617       myPFCluster.seed()));
0618 }
0619
0620 void PFPhotonTranslator::createPreshowerCluster(const reco::PFBlockElement &PFBE,
0621                                                 reco::PreshowerClusterCollection &preshowerClusters,
0622                                                 unsigned plane) const {
0623   const reco::PFClusterRef &myPFClusterRef = PFBE.clusterRef();
0624   preshowerClusters.push_back(reco::PreshowerCluster(
0625       myPFClusterRef->energy(), myPFClusterRef->position(), myPFClusterRef->hitsAndFractions(), plane));
0626 }
0627
0628 void PFPhotonTranslator::createBasicClusterPtrs(
0629     const edm::OrphanHandle<reco::BasicClusterCollection> &basicClustersHandle) {
0630   unsigned size = photPFCandidateIndex_.size();
0631   unsigned basicClusterCounter = 0;
0632   basicClusterPtr_.resize(size);
0633
0634   for (unsigned iphot = 0; iphot < size; ++iphot)  // loop on tracks
0635   {
0636     unsigned nbc = basicClusters_[iphot].size();
0637     for (unsigned ibc = 0; ibc < nbc; ++ibc)  // loop on basic clusters
0638     {
0639       //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
0640       reco::CaloClusterPtr bcPtr(basicClustersHandle, basicClusterCounter);
0641       basicClusterPtr_[iphot].push_back(bcPtr);
0642       ++basicClusterCounter;
0643     }
0644   }
0645 }
0646
0647 void PFPhotonTranslator::createPreshowerClusterPtrs(
0648     const edm::OrphanHandle<reco::PreshowerClusterCollection> &preshowerClustersHandle) {
0649   unsigned size = photPFCandidateIndex_.size();
0650   unsigned psClusterCounter = 0;
0651   preshowerClusterPtr_.resize(size);
0652
0653   for (unsigned iphot = 0; iphot < size; ++iphot)  // loop on tracks
0654   {
0655     unsigned nbc = preshowerClusters_[iphot].size();
0656     for (unsigned ibc = 0; ibc < nbc; ++ibc)  // loop on basic clusters
0657     {
0658       //      std::cout <<  "Track "<< iGSF << " ref " << basicClusterCounter << std::endl;
0659       reco::CaloClusterPtr psPtr(preshowerClustersHandle, psClusterCounter);
0660       preshowerClusterPtr_[iphot].push_back(psPtr);
0661       ++psClusterCounter;
0662     }
0663   }
0664 }
0665
0666 void PFPhotonTranslator::createSuperClusters(const reco::PFCandidateCollection &pfCand,
0667                                              reco::SuperClusterCollection &superClusters) const {
0668   unsigned nphot = photPFCandidateIndex_.size();
0669   for (unsigned iphot = 0; iphot < nphot; ++iphot) {
0670     //cout << "SC iphot=" << iphot << endl;
0671
0672     // Computes energy position a la e/gamma
0673     double sclusterE = 0;
0674     double posX = 0.;
0675     double posY = 0.;
0676     double posZ = 0.;
0677
0678     unsigned nbasics = basicClusters_[iphot].size();
0679     for (unsigned ibc = 0; ibc < nbasics; ++ibc) {
0680       //cout << "BC in SC : iphot="<<iphot<<endl;
0681
0682       double e = basicClusters_[iphot][ibc].energy();
0683       sclusterE += e;
0684       posX += e * basicClusters_[iphot][ibc].position().X();
0685       posY += e * basicClusters_[iphot][ibc].position().Y();
0686       posZ += e * basicClusters_[iphot][ibc].position().Z();
0687     }
0688     posX /= sclusterE;
0689     posY /= sclusterE;
0690     posZ /= sclusterE;
0691
0692     /*
0693       if(pfCand[gsfPFCandidateIndex_[iphot]].gsfTrackRef()!=GsfTrackRef_[iphot])
0694     {
0695       edm::LogError("PFElectronTranslator") << " Major problem in PFElectron Translator" << std::endl;
0696     }
0697       */
0698
0699     // compute the width
0700     PFClusterWidthAlgo pfwidth(pfClusters_[iphot]);
0701
0702     double correctedEnergy = pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
0703     reco::SuperCluster mySuperCluster(correctedEnergy, math::XYZPoint(posX, posY, posZ));
0704     // protection against empty basic cluster collection ; the value is -2 in this case
0705     if (nbasics) {
0706       //      std::cout << "SuperCluster creation; energy " << pfCand[gsfPFCandidateIndex_[iphot]].ecalEnergy();
0707       //      std::cout << " " <<   pfCand[gsfPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
0708       //      std::cout << "Seed energy from basic " << basicClusters_[iphot][0].energy() << std::endl;
0709       mySuperCluster.setSeed(basicClusterPtr_[iphot][0]);
0710     } else {
0711       //      std::cout << "SuperCluster creation ; seed energy " << 0 << std::endl;
0712       //std::cout << "SuperCluster creation ; energy " << pfCand[photPFCandidateIndex_[iphot]].ecalEnergy();
0713       //std::cout << " " <<   pfCand[photPFCandidateIndex_[iphot]].rawEcalEnergy() << std::endl;
0714       //      std::cout << " No seed found " << 0 << std::endl;
0715       //      std::cout << " MVA " << pfCand[gsfPFCandidateIndex_[iphot]].mva_e_pi() << std::endl;
0716       mySuperCluster.setSeed(reco::CaloClusterPtr());
0717     }
0718     // the seed should be the first basic cluster
0719
0720     for (unsigned ibc = 0; ibc < nbasics; ++ibc) {
0721       mySuperCluster.addCluster(basicClusterPtr_[iphot][ibc]);
0722       //      std::cout <<"Adding Ref to SC " << basicClusterPtr_[iphot][ibc].index() << std::endl;
0723       const std::vector<std::pair<DetId, float> > &v1 = basicClusters_[iphot][ibc].hitsAndFractions();
0724       //      std::cout << " Number of cells " << v1.size() << std::endl;
0725       for (std::vector<std::pair<DetId, float> >::const_iterator diIt = v1.begin(); diIt != v1.end(); ++diIt) {
0726         //      std::cout << " Adding DetId " << (diIt->first).rawId() << " " << diIt->second << std::endl;
0727         mySuperCluster.addHitAndFraction(diIt->first, diIt->second);
0728       }  // loop over rechits
0729     }
0730
0731     unsigned nps = preshowerClusterPtr_[iphot].size();
0732     for (unsigned ips = 0; ips < nps; ++ips) {
0733       mySuperCluster.addPreshowerCluster(preshowerClusterPtr_[iphot][ips]);
0734     }
0735
0736     // Set the preshower energy
0737     mySuperCluster.setPreshowerEnergy(pfCand[photPFCandidateIndex_[iphot]].pS1Energy() +
0738                                       pfCand[photPFCandidateIndex_[iphot]].pS2Energy());
0739
0740     // Set the cluster width
0741     mySuperCluster.setEtaWidth(pfwidth.pflowEtaWidth());
0742     mySuperCluster.setPhiWidth(pfwidth.pflowPhiWidth());
0743     // Force the computation of rawEnergy_ of the reco::SuperCluster
0744     mySuperCluster.rawEnergy();
0745
0746     //cout << "SC energy="<< mySuperCluster.energy()<<endl;
0747
0748     superClusters.push_back(mySuperCluster);
0749     //std::cout << "nb super clusters in collection : "<<superClusters.size()<<std::endl;
0750   }
0751 }
0752
0753 void PFPhotonTranslator::createOneLegConversions(
0754     const edm::OrphanHandle<reco::SuperClusterCollection> &superClustersHandle,
0755     reco::ConversionCollection &oneLegConversions) {
0756   //std::cout << "createOneLegConversions" << std::endl;
0757
0758   unsigned nphot = photPFCandidateIndex_.size();
0759   for (unsigned iphot = 0; iphot < nphot; ++iphot) {
0760     //if (conv1legPFCandidateIndex_[iphot]==-1) cout << "No OneLegConversions to add"<<endl;
0761     //else std::cout << "Phot "<<iphot<< " nOneLegConversions to add : "<<pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size()<<endl;
0762
0763     if (conv1legPFCandidateIndex_[iphot] > -1) {
0764       for (unsigned iConv = 0; iConv < pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++) {
0765         reco::CaloClusterPtrVector scPtrVec;
0766         std::vector<reco::CaloClusterPtr> matchingBC;
0767         math::Error<3>::type error;
0768         const reco::Vertex *convVtx =
0769             new reco::Vertex(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition(), error);
0770
0771         //cout << "Vtx x="<<convVtx->x() << " y="<< convVtx->y()<<" z="<<convVtx->z()<< endl;
0772         //cout << "VtxError x=" << convVtx->xError() << endl;
0773
0774         std::vector<reco::TrackRef> OneLegConvVector;
0775         OneLegConvVector.push_back(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]);
0776
0777         std::vector<reco::TrackRef> tr = OneLegConvVector;
0778         std::vector<math::XYZPointF> trackPositionAtEcalVec;
0779         std::vector<math::XYZPointF> innPointVec;
0780         std::vector<math::XYZVectorF> trackPinVec;
0781         std::vector<math::XYZVectorF> trackPoutVec;
0782         math::XYZPointF trackPositionAtEcal(
0783             pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerPosition().X(),
0784             pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerPosition().Y(),
0785             pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerPosition().Z());
0786         math::XYZPointF innPoint(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition().X(),
0787                                  pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition().Y(),
0788                                  pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerPosition().Z());
0789         math::XYZVectorF trackPin(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerMomentum().X(),
0790                                   pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerMomentum().Y(),
0791                                   pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->innerMomentum().Z());
0792         math::XYZVectorF trackPout(pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerMomentum().X(),
0793                                    pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerMomentum().Y(),
0794                                    pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->outerMomentum().Z());
0795         float DCA = pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]][iConv]->d0();
0796         trackPositionAtEcalVec.push_back(trackPositionAtEcal);
0797         innPointVec.push_back(innPoint);
0798         trackPinVec.push_back(trackPin);
0799         trackPoutVec.push_back(trackPout);
0800         std::vector<float> OneLegMvaVector;
0801         reco::Conversion myOneLegConversion(scPtrVec,
0802                                             OneLegConvVector,
0803                                             trackPositionAtEcalVec,
0804                                             *convVtx,
0805                                             matchingBC,
0806                                             DCA,
0807                                             innPointVec,
0808                                             trackPinVec,
0809                                             trackPoutVec,
0810                                             pfSingleLegConvMva_[conv1legPFCandidateIndex_[iphot]][iConv],
0811                                             reco::Conversion::pflow);
0812         OneLegMvaVector.push_back(pfSingleLegConvMva_[conv1legPFCandidateIndex_[iphot]][iConv]);
0813         myOneLegConversion.setOneLegMVA(OneLegMvaVector);
0814         //reco::Conversion myOneLegConversion(scPtrVec,
0815         //OneLegConvVector, *convVtx, reco::Conversion::pflow);
0816
0817         /*
0818         std::cout << "One leg conversion created" << endl;
0819         std::vector<edm::RefToBase<reco::Track> > convtracks = myOneLegConversion.tracks();
0820         const std::vector<float> mvalist = myOneLegConversion.oneLegMVA();
0821
0822         cout << "nTracks="<< convtracks.size()<<endl;
0823         for (unsigned int itk=0; itk<convtracks.size(); itk++){
0824           //double convtrackpt = convtracks[itk]->pt();
0825           std::cout << "Track pt="<< convtracks[itk]->pt() << std::endl;
0826           std::cout << "Track mva="<< mvalist[itk] << std::endl;
0827         }
0828         */
0829         oneLegConversions.push_back(myOneLegConversion);
0830
0831         //cout << "OneLegConv added"<<endl;
0832       }
0833     }
0834   }
0835 }
0836
0837 void PFPhotonTranslator::createPhotonCores(const edm::OrphanHandle<reco::SuperClusterCollection> &superClustersHandle,
0838                                            const edm::OrphanHandle<reco::ConversionCollection> &oneLegConversionHandle,
0839                                            reco::PhotonCoreCollection &photonCores) {
0840   //std::cout << "createPhotonCores" << std::endl;
0841
0842   unsigned nphot = photPFCandidateIndex_.size();
0843
0844   unsigned i1legtot = 0;
0845
0846   for (unsigned iphot = 0; iphot < nphot; ++iphot) {
0847     //std::cout << "iphot="<<iphot<<std::endl;
0848
0849     reco::PhotonCore myPhotonCore;
0850
0851     reco::SuperClusterRef SCref(reco::SuperClusterRef(superClustersHandle, iphot));
0852
0853     myPhotonCore.setPFlowPhoton(true);
0854     myPhotonCore.setStandardPhoton(false);
0855     myPhotonCore.setParentSuperCluster(SCref);
0856     myPhotonCore.setSuperCluster(egSCRef_[iphot]);
0857
0858     reco::ElectronSeedRefVector pixelSeeds = egPhotonRef_[iphot]->electronPixelSeeds();
0859     for (unsigned iseed = 0; iseed < pixelSeeds.size(); iseed++) {
0860       myPhotonCore.addElectronPixelSeed(pixelSeeds[iseed]);
0861     }
0862
0863     //cout << "PhotonCores : SC OK" << endl;
0864
0865     //cout << "conv1legPFCandidateIndex_[iphot]="<<conv1legPFCandidateIndex_[iphot]<<endl;
0866     //cout << "conv2legPFCandidateIndex_[iphot]="<<conv2legPFCandidateIndex_[iphot]<<endl;
0867
0868     if (conv1legPFCandidateIndex_[iphot] > -1) {
0869       for (unsigned int iConv = 0; iConv < pfSingleLegConv_[conv1legPFCandidateIndex_[iphot]].size(); iConv++) {
0870         const reco::ConversionRef &OneLegRef(reco::ConversionRef(oneLegConversionHandle, i1legtot));
0871         myPhotonCore.addOneLegConversion(OneLegRef);
0872
0873         //cout << "PhotonCores : 1-leg OK" << endl;
0874         /*
0875       cout << "Testing 1-leg :"<<endl;
0876       const reco::ConversionRefVector & conv = myPhotonCore.conversionsOneLeg();
0877       for (unsigned int iconv=0; iconv<conv.size(); iconv++){
0878         cout << "Testing 1-leg : iconv="<<iconv<<endl;
0879         cout << "Testing 1-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
0880         cout << "Testing 1-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
0881         std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
0882         for (unsigned int itk=0; itk<convtracks.size(); itk++){
0883           //double convtrackpt = convtracks[itk]->pt();
0884           std::cout << "Testing 1-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
0885           // std::cout << "Track mva="<< mvalist[itk] << std::endl;
0886         }
0887       }
0888       */
0889
0890         i1legtot++;
0891       }
0892     }
0893
0894     if (conv2legPFCandidateIndex_[iphot] > -1) {
0895       for (unsigned int iConv = 0; iConv < pfConv_[conv2legPFCandidateIndex_[iphot]].size(); iConv++) {
0896         const reco::ConversionRef &TwoLegRef(pfConv_[conv2legPFCandidateIndex_[iphot]][iConv]);
0897         myPhotonCore.addConversion(TwoLegRef);
0898       }
0899       //cout << "PhotonCores : 2-leg OK" << endl;
0900
0901       /*
0902     cout << "Testing 2-leg :"<<endl;
0903     const reco::ConversionRefVector & conv = myPhotonCore.conversions();
0904     for (unsigned int iconv=0; iconv<conv.size(); iconv++){
0905       cout << "Testing 2-leg : iconv="<<iconv<<endl;
0906       cout << "Testing 2-leg : nTracks="<<conv[iconv]->nTracks()<<endl;
0907       cout << "Testing 2-leg : EoverP="<<conv[iconv]->EoverP()<<endl;
0908       std::vector<edm::RefToBase<reco::Track> > convtracks = conv[iconv]->tracks();
0909       for (unsigned int itk=0; itk<convtracks.size(); itk++){
0910         //double convtrackpt = convtracks[itk]->pt();
0911         std::cout << "Testing 2-leg : Track pt="<< convtracks[itk]->pt() << std::endl;
0912         // std::cout << "Track mva="<< mvalist[itk] << std::endl;
0913       }
0914     }
0915     */
0916     }
0917
0918     photonCores.push_back(myPhotonCore);
0919   }
0920
0921   //std::cout << "end of createPhotonCores"<<std::endl;
0922 }
0923
0924 void PFPhotonTranslator::createPhotons(reco::VertexCollection &vertexCollection,
0925                                        edm::Handle<reco::PhotonCollection> &egPhotons,
0926                                        const edm::OrphanHandle<reco::PhotonCoreCollection> &photonCoresHandle,
0927                                        const IsolationValueMaps &isolationValues,
0928                                        reco::PhotonCollection &photons) {
0929   //cout << "createPhotons" << endl;
0930
0931   unsigned nphot = photPFCandidateIndex_.size();
0932
0933   for (unsigned iphot = 0; iphot < nphot; ++iphot) {
0934     //std::cout << "iphot="<<iphot<<std::endl;
0935
0936     reco::PhotonCoreRef PCref(reco::PhotonCoreRef(photonCoresHandle, iphot));
0937
0938     math::XYZPoint vtx(0., 0., 0.);
0939     if (!vertexCollection.empty())
0940       vtx = vertexCollection.begin()->position();
0941     //std::cout << "vtx made" << std::endl;
0942
0943     math::XYZVector direction = PCref->parentSuperCluster()->position() - vtx;
0944
0945     //It could be that pfSC energy gives not the best resolution : use smaller agregates for some cases ?
0946     math::XYZVector P3 = direction.unit() * PCref->parentSuperCluster()->energy();
0947     LorentzVector P4(P3.x(), P3.y(), P3.z(), PCref->parentSuperCluster()->energy());
0948
0949     reco::Photon myPhoton(P4, PCref->parentSuperCluster()->position(), PCref, vtx);
0950     //cout << "photon created"<<endl;
0951
0952     reco::Photon::ShowerShape showerShape;
0953     reco::Photon::SaturationInfo saturationInfo;
0954     reco::Photon::FiducialFlags fiducialFlags;
0955     reco::Photon::IsolationVariables isolationVariables03;
0956     reco::Photon::IsolationVariables isolationVariables04;
0957
0958     showerShape.e1x5 = egPhotonRef_[iphot]->e1x5();
0959     showerShape.e2x5 = egPhotonRef_[iphot]->e2x5();
0960     showerShape.e3x3 = egPhotonRef_[iphot]->e3x3();
0961     showerShape.e5x5 = egPhotonRef_[iphot]->e5x5();
0962     showerShape.maxEnergyXtal = egPhotonRef_[iphot]->maxEnergyXtal();
0963     showerShape.sigmaEtaEta = egPhotonRef_[iphot]->sigmaEtaEta();
0964     showerShape.sigmaIetaIeta = egPhotonRef_[iphot]->sigmaIetaIeta();
0965     for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) {
0966       showerShape.hcalOverEcal[id] = egPhotonRef_[iphot]->hcalOverEcal(id + 1);
0967     }
0968     myPhoton.setShowerShapeVariables(showerShape);
0969
0970     saturationInfo.nSaturatedXtals = egPhotonRef_[iphot]->nSaturatedXtals();
0971     saturationInfo.isSeedSaturated = egPhotonRef_[iphot]->isSeedSaturated();
0972     myPhoton.setSaturationInfo(saturationInfo);
0973
0974     fiducialFlags.isEB = egPhotonRef_[iphot]->isEB();
0975     fiducialFlags.isEE = egPhotonRef_[iphot]->isEE();
0976     fiducialFlags.isEBEtaGap = egPhotonRef_[iphot]->isEBEtaGap();
0977     fiducialFlags.isEBPhiGap = egPhotonRef_[iphot]->isEBPhiGap();
0978     fiducialFlags.isEERingGap = egPhotonRef_[iphot]->isEERingGap();
0979     fiducialFlags.isEEDeeGap = egPhotonRef_[iphot]->isEEDeeGap();
0980     fiducialFlags.isEBEEGap = egPhotonRef_[iphot]->isEBEEGap();
0981     myPhoton.setFiducialVolumeFlags(fiducialFlags);
0982
0983     isolationVariables03.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR03();
0984     for (uint id = 0; id < isolationVariables03.hcalRecHitSumEt.size(); ++id) {
0985       isolationVariables03.hcalRecHitSumEt[id] = egPhotonRef_[iphot]->hcalTowerSumEtConeDR03(id + 1);
0986     }
0987     isolationVariables03.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR03();
0988     isolationVariables03.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR03();
0989     isolationVariables03.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR03();
0990     isolationVariables03.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR03();
0991     isolationVariables04.ecalRecHitSumEt = egPhotonRef_[iphot]->ecalRecHitSumEtConeDR04();
0992     for (uint id = 0; id < isolationVariables04.hcalRecHitSumEt.size(); ++id) {
0993       isolationVariables04.hcalRecHitSumEt[id] = egPhotonRef_[iphot]->hcalTowerSumEtConeDR04(id + 1);
0994     }
0995     isolationVariables04.trkSumPtSolidCone = egPhotonRef_[iphot]->trkSumPtSolidConeDR04();
0996     isolationVariables04.trkSumPtHollowCone = egPhotonRef_[iphot]->trkSumPtHollowConeDR04();
0997     isolationVariables04.nTrkSolidCone = egPhotonRef_[iphot]->nTrkSolidConeDR04();
0998     isolationVariables04.nTrkHollowCone = egPhotonRef_[iphot]->nTrkHollowConeDR04();
0999     myPhoton.setIsolationVariables(isolationVariables04, isolationVariables03);
1000
1001     reco::Photon::PflowIsolationVariables myPFIso;
1002     myPFIso.chargedHadronIso = (*isolationValues[0])[CandidatePtr_[iphot]];
1003     myPFIso.photonIso = (*isolationValues[1])[CandidatePtr_[iphot]];
1004     myPFIso.neutralHadronIso = (*isolationValues[2])[CandidatePtr_[iphot]];
1005     myPhoton.setPflowIsolationVariables(myPFIso);
1006
1007     reco::Photon::PflowIDVariables myPFVariables;
1008
1009     reco::Mustache myMustache(mustacheSCParams_);
1010     myMustache.MustacheID(
1011         *(myPhoton.parentSuperCluster()), myPFVariables.nClusterOutsideMustache, myPFVariables.etOutsideMustache);
1012     myPFVariables.mva = pfPhotonMva_[iphot];
1013     myPhoton.setPflowIDVariables(myPFVariables);
1014
1015     //cout << "chargedHadronIso="<<myPhoton.chargedHadronIso()<<" photonIso="<<myPhoton.photonIso()<<" neutralHadronIso="<<myPhoton.neutralHadronIso()<<endl;
1016
1017     // set PF-regression energy
1018     myPhoton.setCorrectedEnergy(
1019         reco::Photon::regression2, energyRegression_[iphot], energyRegressionError_[iphot], false);
1020
1021     /*
1022       if (basicClusters_[iphot].size()>0){
1023       // Cluster shape variables
1024       //Algorithms from EcalClusterTools could be adapted to PF photons ? (not based on 5x5 BC)
1025       //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)
1026       const EcalRecHitCollection* hits = 0 ;
1027       int subdet = PCref->parentSuperCluster()->seed()->hitsAndFractions()[0].first.subdetId();
1028       if (subdet==EcalBarrel) hits = barrelRecHits;
1029       else if  (subdet==EcalEndcap) hits = endcapRecHits;
1030       const CaloGeometry* geometry = theCaloGeom_.product();
1031
1032       float maxXtal =   EcalClusterTools::eMax( *(PCref->parentSuperCluster()->seed()), &(*hits) );
1033       //cout << "maxXtal="<<maxXtal<<endl;
1034       float e1x5    =   EcalClusterTools::e1x5(  *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology));
1035       //cout << "e1x5="<<e1x5<<endl;
1036       float e2x5    =   EcalClusterTools::e2x5Max(  *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology));
1037       //cout << "e2x5="<<e2x5<<endl;
1038       float e3x3    =   EcalClusterTools::e3x3(  *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology));
1039       //cout << "e3x3="<<e3x3<<endl;
1040       float e5x5    =   EcalClusterTools::e5x5( *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology));
1041       //cout << "e5x5="<<e5x5<<endl;
1042       std::vector<float> cov =  EcalClusterTools::covariances( *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology), geometry);
1043       float sigmaEtaEta = sqrt(cov[0]);
1044       //cout << "sigmaEtaEta="<<sigmaEtaEta<<endl;
1045       std::vector<float> locCov =  EcalClusterTools::localCovariances( *(PCref->parentSuperCluster()->seed()), &(*hits), &(*topology));
1046       float sigmaIetaIeta = sqrt(locCov[0]);
1047       //cout << "sigmaIetaIeta="<<sigmaIetaIeta<<endl;
1048       //float r9 =e3x3/(PCref->parentSuperCluster()->rawEnergy());
1049
1050
1051       // calculate HoE
1052       const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
1053       EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;
1054       EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;
1055       double HoE1=towerIso1.getTowerESum(&(*PCref->parentSuperCluster()))/PCref->pfSuperCluster()->energy();
1056       double HoE2=towerIso2.getTowerESum(&(*PCref->parentSuperCluster()))/PCref->pfSuperCluster()->energy();
1057       //cout << "HoE1="<<HoE1<<endl;
1058       //cout << "HoE2="<<HoE2<<endl;
1059
1060       reco::Photon::ShowerShape  showerShape;
1061       showerShape.e1x5= e1x5;
1062       showerShape.e2x5= e2x5;
1063       showerShape.e3x3= e3x3;
1064       showerShape.e5x5= e5x5;
1065       showerShape.maxEnergyXtal =  maxXtal;
1066       showerShape.sigmaEtaEta =    sigmaEtaEta;
1067       showerShape.sigmaIetaIeta =  sigmaIetaIeta;
1068       showerShape.hcalDepth1OverEcal = HoE1;
1069       showerShape.hcalDepth2OverEcal = HoE2;
1070       myPhoton.setShowerShapeVariables ( showerShape );
1071       //cout << "shower shape variables filled"<<endl;
1072       }
1073       */
1074
1075     photons.push_back(myPhoton);
1076   }
1077
1078   //std::cout << "end of createPhotons"<<std::endl;
1079 }
1080
1081 const reco::PFCandidate &PFPhotonTranslator::correspondingDaughterCandidate(const reco::PFCandidate &cand,
1082                                                                             const reco::PFBlockElement &pfbe) const {
1083   unsigned refindex = pfbe.index();
1084   //  std::cout << " N daughters " << cand.numberOfDaughters() << std::endl;
1085   reco::PFCandidate::const_iterator myDaughterCandidate = cand.begin();
1086   reco::PFCandidate::const_iterator itend = cand.end();
1087
1088   for (; myDaughterCandidate != itend; ++myDaughterCandidate) {
1089     const reco::PFCandidate *myPFCandidate = (const reco::PFCandidate *)&*myDaughterCandidate;
1090     if (myPFCandidate->elementsInBlocks().size() != 1) {
1091       //      std::cout << " Daughter with " << myPFCandidate.elementsInBlocks().size()<< " element in block " << std::endl;
1092       return cand;
1093     }
1094     if (myPFCandidate->elementsInBlocks()[0].second == refindex) {
1095       //      std::cout << " Found it " << cand << std::endl;
1096       return *myPFCandidate;
1097     }
1098   }
1099   return cand;
1100 }