CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/RecoParticleFlow/PFProducer/src/PFPhotonAlgo.cc

Go to the documentation of this file.
00001 //
00002 // Original Authors: Fabian Stoeckli: fabian.stoeckli@cern.ch
00003 //                   Nicholas Wardle: nckw@cern.ch
00004 //                   Rishi Patel rpatel@cern.ch
00005 //
00006 
00007 #include "RecoParticleFlow/PFProducer/interface/PFPhotonAlgo.h"
00008 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
00009 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00010 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
00011 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementSuperCluster.h"
00012 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00013 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
00014 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00015 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00016 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00017 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
00018 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00019 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00020 #include "DataFormats/Math/interface/deltaPhi.h"
00021 #include "DataFormats/Math/interface/deltaR.h"
00022 #include <TFile.h>
00023 #include <iomanip>
00024 #include <algorithm>
00025 #include <TMath.h>
00026 using namespace std;
00027 using namespace reco;
00028 
00029 
00030 PFPhotonAlgo::PFPhotonAlgo(std::string mvaweightfile,  
00031                            double mvaConvCut, 
00032                            bool useReg,
00033                            std::string X0_Map,
00034                            const reco::Vertex& primary,
00035                            const boost::shared_ptr<PFEnergyCalibration>& thePFEnergyCalibration,
00036                            double sumPtTrackIsoForPhoton,
00037                            double sumPtTrackIsoSlopeForPhoton
00038                            ) : 
00039   isvalid_(false), 
00040   verbosityLevel_(Silent), 
00041   MVACUT(mvaConvCut),
00042   useReg_(useReg),
00043   thePFEnergyCalibration_(thePFEnergyCalibration),
00044   sumPtTrackIsoForPhoton_(sumPtTrackIsoForPhoton),
00045   sumPtTrackIsoSlopeForPhoton_(sumPtTrackIsoSlopeForPhoton)
00046 {  
00047     primaryVertex_=primary;  
00048     //Book MVA  
00049     tmvaReader_ = new TMVA::Reader("!Color:Silent");  
00050     tmvaReader_->AddVariable("del_phi",&del_phi);  
00051     tmvaReader_->AddVariable("nlayers", &nlayers);  
00052     tmvaReader_->AddVariable("chi2",&chi2);  
00053     tmvaReader_->AddVariable("EoverPt",&EoverPt);  
00054     tmvaReader_->AddVariable("HoverPt",&HoverPt);  
00055     tmvaReader_->AddVariable("track_pt", &track_pt);  
00056     tmvaReader_->AddVariable("STIP",&STIP);  
00057     tmvaReader_->AddVariable("nlost", &nlost);  
00058     tmvaReader_->BookMVA("BDT",mvaweightfile.c_str());  
00059 
00060     //Material Map
00061     TFile *XO_File = new TFile(X0_Map.c_str(),"READ");
00062     X0_sum=(TH2D*)XO_File->Get("TrackerSum");
00063     X0_inner = (TH2D*)XO_File->Get("Inner");
00064     X0_middle = (TH2D*)XO_File->Get("Middle");
00065     X0_outer = (TH2D*)XO_File->Get("Outer");
00066     
00067 }
00068 
00069 void PFPhotonAlgo::RunPFPhoton(const reco::PFBlockRef&  blockRef,
00070                                std::vector<bool>& active,
00071                                std::auto_ptr<PFCandidateCollection> &pfCandidates,
00072                                std::vector<reco::PFCandidatePhotonExtra>& pfPhotonExtraCandidates,
00073                                std::vector<reco::PFCandidate> 
00074                                &tempElectronCandidates
00075 ){
00076   
00077   //std::cout<<" calling RunPFPhoton "<<std::endl;
00078   
00079   /*      For now we construct the PhotonCandidate simply from 
00080           a) adding the CORRECTED energies of each participating ECAL cluster
00081           b) build the energy-weighted direction for the Photon
00082   */
00083 
00084 
00085   // define how much is printed out for debugging.
00086   // ... will be setable via CFG file parameter
00087   verbosityLevel_ = Chatty;          // Chatty mode.
00088   
00089 
00090   // loop over all elements in the Block
00091   const edm::OwnVector< reco::PFBlockElement >&          elements         = blockRef->elements();
00092   edm::OwnVector< reco::PFBlockElement >::const_iterator ele              = elements.begin();
00093   std::vector<bool>::const_iterator                      actIter          = active.begin();
00094   PFBlock::LinkData                                      linkData         = blockRef->linkData();
00095   bool                                                   isActive         = true;
00096 
00097 
00098   if(elements.size() != active.size()) {
00099     // throw excpetion...
00100     //std::cout<<" WARNING: Size of collection and active-vectro don't agree!"<<std::endl;
00101     return;
00102   }
00103   
00104   // local vecotr to keep track of the indices of the 'elements' for the Photon candidate
00105   // once we decide to keep the candidate, the 'active' entriesd for them must be set to false
00106   std::vector<unsigned int> elemsToLock;
00107   elemsToLock.resize(0);
00108   
00109   for( ; ele != elements.end(); ++ele, ++actIter ) {
00110 
00111     // if it's not a SuperCluster, go to the next element
00112     if( !( ele->type() == reco::PFBlockElement::SC ) ) continue;
00113     
00114     // Photon kienmatics, will be updated for each identified participating element
00115     float photonEnergy_        =   0.;
00116     float photonX_             =   0.;
00117     float photonY_             =   0.;
00118     float photonZ_             =   0.;
00119     float RawEcalEne           =   0.;
00120 
00121     // Total pre-shower energy
00122     float ps1TotEne      = 0.;
00123     float ps2TotEne      = 0.;
00124     
00125     bool hasConvTrack=false;  
00126     bool hasSingleleg=false;  
00127     std::vector<unsigned int> AddClusters(0);  
00128     std::vector<unsigned int> IsoTracks(0);  
00129     std::multimap<unsigned int, unsigned int>ClusterAddPS1;  
00130     std::multimap<unsigned int, unsigned int>ClusterAddPS2;
00131     std::vector<reco::TrackRef>singleLegRef;
00132     std::vector<float>MVA_values(0);
00133     std::vector<float>MVALCorr;
00134     std::vector<const CaloCluster*>PFClusters;
00135     reco::ConversionRefVector ConversionsRef_;
00136     isActive = *(actIter);
00137     //cout << " Found a SuperCluster.  Energy " ;
00138     const reco::PFBlockElementSuperCluster *sc = dynamic_cast<const reco::PFBlockElementSuperCluster*>(&(*ele));
00139     //std::cout << sc->superClusterRef()->energy () << " Track/Ecal/Hcal Iso " << sc->trackIso()<< " " << sc->ecalIso() ;
00140     //std::cout << " " << sc->hcalIso() <<std::endl;
00141     if (!(sc->fromPhoton()))continue;
00142 
00143     // check the status of the SC Element... 
00144     // ..... I understand it should *always* be active, since PFElectronAlgo does not touch this (yet?) RISHI: YES
00145     if( !isActive ) {
00146       //std::cout<<" SuperCluster is NOT active.... "<<std::endl;
00147       continue;
00148     }
00149     elemsToLock.push_back(ele-elements.begin()); //add SC to elements to lock
00150     // loop over its constituent ECAL cluster
00151     std::multimap<double, unsigned int> ecalAssoPFClusters;
00152     blockRef->associatedElements( ele-elements.begin(), 
00153                                   linkData,
00154                                   ecalAssoPFClusters,
00155                                   reco::PFBlockElement::ECAL,
00156                                   reco::PFBlock::LINKTEST_ALL );
00157     
00158     // loop over the ECAL clusters linked to the iEle 
00159     if( ! ecalAssoPFClusters.size() ) {
00160       // This SC element has NO ECAL elements asigned... *SHOULD NOT HAPPEN*
00161       //std::cout<<" Found SC element with no ECAL assigned "<<std::endl;
00162       continue;
00163     }
00164     
00165     // This is basically CASE 2
00166     // .... we loop over all ECAL cluster linked to each other by this SC
00167     for(std::multimap<double, unsigned int>::iterator itecal = ecalAssoPFClusters.begin(); 
00168         itecal != ecalAssoPFClusters.end(); ++itecal) { 
00169       
00170       // to get the reference to the PF clusters, this is needed.
00171       reco::PFClusterRef clusterRef = elements[itecal->second].clusterRef();    
00172       
00173       // from the clusterRef get the energy, direction, etc
00174       //      float ClustRawEnergy = clusterRef->energy();
00175       //      float ClustEta = clusterRef->position().eta();
00176       //      float ClustPhi = clusterRef->position().phi();
00177 
00178       // initialize the vectors for the PS energies
00179       vector<double> ps1Ene(0);
00180       vector<double> ps2Ene(0);
00181       double ps1=0;  
00182       double ps2=0;  
00183       hasSingleleg=false;  
00184       hasConvTrack=false;
00185 
00186       /*
00187       cout << " My cluster index " << itecal->second 
00188            << " energy " <<  ClustRawEnergy
00189            << " eta " << ClustEta
00190            << " phi " << ClustPhi << endl;
00191       */
00192       // check if this ECAL element is still active (could have been eaten by PFElectronAlgo)
00193       // ......for now we give the PFElectron Algo *ALWAYS* Shot-Gun on the ECAL elements to the PFElectronAlgo
00194       
00195       if( !( active[itecal->second] ) ) {
00196         //std::cout<< "  .... this ECAL element is NOT active anymore. Is skipped. "<<std::endl;
00197         continue;
00198       }
00199       
00200       // ------------------------------------------------------------------------------------------
00201       // TODO: do some tests on the ECAL cluster itself, deciding to use it or not for the Photons
00202       // ..... ??? Do we need this?
00203       if ( false ) {
00204         // Check if there are a large number tracks that do not pass pre-ID around this ECAL cluster
00205         bool useIt = true;
00206         int mva_reject=0;  
00207         bool isClosest=false;  
00208         std::multimap<double, unsigned int> Trackscheck;  
00209         blockRef->associatedElements( itecal->second,  
00210                                       linkData,  
00211                                       Trackscheck,  
00212                                       reco::PFBlockElement::TRACK,  
00213                                       reco::PFBlock::LINKTEST_ALL);  
00214         for(std::multimap<double, unsigned int>::iterator track = Trackscheck.begin();  
00215             track != Trackscheck.end(); ++track) {  
00216            
00217           // first check if is it's still active  
00218           if( ! (active[track->second]) ) continue;  
00219           hasSingleleg=EvaluateSingleLegMVA(blockRef,  primaryVertex_, track->second);  
00220           //check if it is the closest linked track  
00221           std::multimap<double, unsigned int> closecheck;  
00222           blockRef->associatedElements(track->second,  
00223                                        linkData,  
00224                                        closecheck,  
00225                                        reco::PFBlockElement::ECAL,  
00226                                        reco::PFBlock::LINKTEST_ALL);  
00227           if(closecheck.begin()->second ==itecal->second)isClosest=true;  
00228           if(!hasSingleleg)mva_reject++;  
00229         }  
00230          
00231         if(mva_reject>0 &&  isClosest)useIt=false;  
00232         //if(mva_reject==1 && isClosest)useIt=false;
00233         if( !useIt ) continue;    // Go to next ECAL cluster within SC
00234       }
00235       // ------------------------------------------------------------------------------------------
00236       
00237       // We decided to keep the ECAL cluster for this Photon Candidate ...
00238       elemsToLock.push_back(itecal->second);
00239       
00240       // look for PS in this Block linked to this ECAL cluster      
00241       std::multimap<double, unsigned int> PS1Elems;
00242       std::multimap<double, unsigned int> PS2Elems;
00243       //PS Layer 1 linked to ECAL cluster
00244       blockRef->associatedElements( itecal->second,
00245                                     linkData,
00246                                     PS1Elems,
00247                                     reco::PFBlockElement::PS1,
00248                                     reco::PFBlock::LINKTEST_ALL );
00249       //PS Layer 2 linked to the ECAL cluster
00250       blockRef->associatedElements( itecal->second,
00251                                     linkData,
00252                                     PS2Elems,
00253                                     reco::PFBlockElement::PS2,
00254                                     reco::PFBlock::LINKTEST_ALL );
00255       
00256       // loop over all PS1 and compute energy
00257       for(std::multimap<double, unsigned int>::iterator iteps = PS1Elems.begin();
00258           iteps != PS1Elems.end(); ++iteps) {
00259 
00260         // first chekc if it's still active
00261         if( !(active[iteps->second]) ) continue;
00262         
00263         //Check if this PS1 is not closer to another ECAL cluster in this Block          
00264         std::multimap<double, unsigned int> ECALPS1check;  
00265         blockRef->associatedElements( iteps->second,  
00266                                       linkData,  
00267                                       ECALPS1check,  
00268                                       reco::PFBlockElement::ECAL,  
00269                                       reco::PFBlock::LINKTEST_ALL );  
00270         if(itecal->second==ECALPS1check.begin()->second)//then it is closest linked  
00271           {
00272             reco::PFClusterRef ps1ClusterRef = elements[iteps->second].clusterRef();
00273             ps1Ene.push_back( ps1ClusterRef->energy() );
00274             ps1=ps1+ps1ClusterRef->energy(); //add to total PS1
00275             // incativate this PS1 Element
00276             elemsToLock.push_back(iteps->second);
00277           }
00278       }
00279       for(std::multimap<double, unsigned int>::iterator iteps = PS2Elems.begin();
00280           iteps != PS2Elems.end(); ++iteps) {
00281 
00282         // first chekc if it's still active
00283         if( !(active[iteps->second]) ) continue;
00284         
00285         // Check if this PS2 is not closer to another ECAL cluster in this Block:
00286         std::multimap<double, unsigned int> ECALPS2check;  
00287         blockRef->associatedElements( iteps->second,  
00288                                       linkData,  
00289                                       ECALPS2check,  
00290                                       reco::PFBlockElement::ECAL,  
00291                                       reco::PFBlock::LINKTEST_ALL );  
00292         if(itecal->second==ECALPS2check.begin()->second)//is closest linked  
00293           {
00294             reco::PFClusterRef ps2ClusterRef = elements[iteps->second].clusterRef();
00295             ps2Ene.push_back( ps2ClusterRef->energy() );
00296             ps2=ps2ClusterRef->energy()+ps2; //add to total PS2
00297             // incativate this PS2 Element
00298             elemsToLock.push_back(iteps->second);
00299           }
00300       }
00301             
00302       // loop over the HCAL Clusters linked to the ECAL cluster (CASE 6)
00303       std::multimap<double, unsigned int> hcalElems;
00304       blockRef->associatedElements( itecal->second,linkData,
00305                                     hcalElems,
00306                                     reco::PFBlockElement::HCAL,
00307                                     reco::PFBlock::LINKTEST_ALL );
00308 
00309       for(std::multimap<double, unsigned int>::iterator ithcal = hcalElems.begin();
00310           ithcal != hcalElems.end(); ++ithcal) {
00311 
00312         if ( ! (active[ithcal->second] ) ) continue; // HCAL Cluster already used....
00313         
00314         // TODO: Decide if this HCAL cluster is to be used
00315         // .... based on some Physics
00316         // .... To we need to check if it's closer to any other ECAL/TRACK?
00317 
00318         bool useHcal = false;
00319         if ( !useHcal ) continue;
00320         //not locked
00321         //elemsToLock.push_back(ithcal->second);
00322       }
00323 
00324       // This is entry point for CASE 3.
00325       // .... we loop over all Tracks linked to this ECAL and check if it's labeled as conversion
00326       // This is the part for looping over all 'Conversion' Tracks
00327       std::multimap<double, unsigned int> convTracks;
00328       blockRef->associatedElements( itecal->second,
00329                                     linkData,
00330                                     convTracks,
00331                                     reco::PFBlockElement::TRACK,
00332                                     reco::PFBlock::LINKTEST_ALL);
00333       for(std::multimap<double, unsigned int>::iterator track = convTracks.begin();
00334           track != convTracks.end(); ++track) {
00335 
00336         // first check if is it's still active
00337         if( ! (active[track->second]) ) continue;
00338         
00339         // check if it's a CONV track
00340         const reco::PFBlockElementTrack * trackRef = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[track->second]));        
00341         
00342         //Check if track is a Single leg from a Conversion  
00343         mvaValue=-999;  
00344         hasSingleleg=EvaluateSingleLegMVA(blockRef,  primaryVertex_, track->second);  
00345 
00346         // Daniele; example for mvaValues, do the same for single leg trackRef and convRef
00347         //          
00348         //      if(hasSingleleg)
00349         //        mvaValues.push_back(mvaValue);
00350 
00351         //If it is not then it will be used to check Track Isolation at the end  
00352         if(!hasSingleleg)  
00353           {  
00354             bool included=false;  
00355             //check if this track is already included in the vector so it is linked to an ECAL cluster that is already examined  
00356             for(unsigned int i=0; i<IsoTracks.size(); i++)  
00357               {if(IsoTracks[i]==track->second)included=true;}  
00358             if(!included)IsoTracks.push_back(track->second);  
00359           }  
00360         //For now only Pre-ID tracks that are not already identified as Conversions  
00361         if(hasSingleleg &&!(trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)))  
00362           {  
00363             elemsToLock.push_back(track->second);
00364             
00365             reco::TrackRef t_ref=elements[track->second].trackRef();
00366             bool matched=false;
00367             for(unsigned int ic=0; ic<singleLegRef.size(); ic++)
00368               if(singleLegRef[ic]==t_ref)matched=true;
00369             
00370             if(!matched){
00371               singleLegRef.push_back(t_ref);
00372               MVA_values.push_back(mvaValue);
00373             }
00374             //find all the clusters linked to this track  
00375             std::multimap<double, unsigned int> moreClusters;  
00376             blockRef->associatedElements( track->second,  
00377                                           linkData,  
00378                                           moreClusters,  
00379                                           reco::PFBlockElement::ECAL,  
00380                                           reco::PFBlock::LINKTEST_ALL);  
00381              
00382             float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x() * elements[track->second].trackRef()->innerMomentum().x() +  
00383                             elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+  
00384                             elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());  
00385             float linked_E=0;  
00386             for(std::multimap<double, unsigned int>::iterator clust = moreClusters.begin();  
00387                 clust != moreClusters.end(); ++clust)  
00388               {  
00389                 if(!active[clust->second])continue;  
00390                 //running sum of linked energy  
00391                 linked_E=linked_E+elements[clust->second].clusterRef()->energy();  
00392                 //prevent too much energy from being added  
00393                 if(linked_E/p_in>1.5)break;  
00394                 bool included=false;  
00395                 //check if these ecal clusters are already included with the supercluster  
00396                 for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();  
00397                     cluscheck != ecalAssoPFClusters.end(); ++cluscheck)  
00398                   {  
00399                     if(cluscheck->second==clust->second)included=true;  
00400                   }  
00401                 if(!included)AddClusters.push_back(clust->second);//Add to a container of clusters to be Added to the Photon candidate  
00402               }  
00403           }
00404 
00405         // Possibly need to be more smart about them (CASE 5)
00406         // .... for now we simply skip non id'ed tracks
00407         if( ! (trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) ) ) continue;  
00408         hasConvTrack=true;  
00409         elemsToLock.push_back(track->second);
00410         //again look at the clusters linked to this track  
00411         //if(elements[track->second].convRef().isNonnull())
00412         //{         
00413         //  ConversionsRef_.push_back(elements[track->second].convRef());
00414         //}
00415         std::multimap<double, unsigned int> moreClusters;  
00416         blockRef->associatedElements( track->second,  
00417                                       linkData,  
00418                                       moreClusters,  
00419                                       reco::PFBlockElement::ECAL,  
00420                                       reco::PFBlock::LINKTEST_ALL);
00421         
00422         float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x() * elements[track->second].trackRef()->innerMomentum().x() +  
00423                         elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+  
00424                         elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());  
00425         float linked_E=0;  
00426         for(std::multimap<double, unsigned int>::iterator clust = moreClusters.begin();  
00427             clust != moreClusters.end(); ++clust)  
00428           {  
00429             if(!active[clust->second])continue;  
00430             linked_E=linked_E+elements[clust->second].clusterRef()->energy();  
00431             if(linked_E/p_in>1.5)break;  
00432             bool included=false;  
00433             for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();  
00434                 cluscheck != ecalAssoPFClusters.end(); ++cluscheck)  
00435               {  
00436                 if(cluscheck->second==clust->second)included=true;  
00437               }  
00438             if(!included)AddClusters.push_back(clust->second);//again only add if it is not already included with the supercluster  
00439           }
00440         
00441         // we need to check for other TRACKS linked to this conversion track, that point possibly no an ECAL cluster not included in the SC
00442         // .... This is basically CASE 4.
00443         
00444         std::multimap<double, unsigned int> moreTracks;
00445         blockRef->associatedElements( track->second,
00446                                       linkData,
00447                                       moreTracks,
00448                                       reco::PFBlockElement::TRACK,
00449                                       reco::PFBlock::LINKTEST_ALL);
00450         
00451         for(std::multimap<double, unsigned int>::iterator track2 = moreTracks.begin();
00452             track2 != moreTracks.end(); ++track2) {
00453           
00454           // first check if is it's still active
00455           if( ! (active[track2->second]) ) continue;
00456           //skip over the 1st leg already found above  
00457           if(track->second==track2->second)continue;      
00458           // check if it's a CONV track
00459           const reco::PFBlockElementTrack * track2Ref = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[track2->second]));    
00460           if( ! (track2Ref->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) ) ) continue;  // Possibly need to be more smart about them (CASE 5)
00461           elemsToLock.push_back(track2->second);
00462           // so it's another active conversion track, that is in the Block and linked to the conversion track we already found
00463           // find the ECAL cluster linked to it...
00464           std::multimap<double, unsigned int> convEcal;
00465           blockRef->associatedElements( track2->second,
00466                                         linkData,
00467                                         convEcal,
00468                                         reco::PFBlockElement::ECAL,
00469                                         reco::PFBlock::LINKTEST_ALL);
00470           float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x()*elements[track->second].trackRef()->innerMomentum().x()+
00471                           elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+  
00472                           elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());  
00473           
00474           
00475           float linked_E=0;
00476           for(std::multimap<double, unsigned int>::iterator itConvEcal = convEcal.begin();
00477               itConvEcal != convEcal.end(); ++itConvEcal) {
00478             
00479             if( ! (active[itConvEcal->second]) ) continue;
00480             bool included=false;  
00481             for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();  
00482                 cluscheck != ecalAssoPFClusters.end(); ++cluscheck)  
00483               {  
00484                 if(cluscheck->second==itConvEcal->second)included=true;  
00485               }
00486             linked_E=linked_E+elements[itConvEcal->second].clusterRef()->energy();
00487             if(linked_E/p_in>1.5)break;
00488             if(!included){AddClusters.push_back(itConvEcal->second);
00489             }
00490             
00491             // it's still active, so we have to add it.
00492             // CAUTION: we don't care here if it's part of the SC or not, we include it anyways
00493             
00494             // loop over the HCAL Clusters linked to the ECAL cluster (CASE 6)
00495             std::multimap<double, unsigned int> hcalElems_conv;
00496             blockRef->associatedElements( itecal->second,linkData,
00497                                           hcalElems_conv,
00498                                           reco::PFBlockElement::HCAL,
00499                                           reco::PFBlock::LINKTEST_ALL );
00500             
00501             for(std::multimap<double, unsigned int>::iterator ithcal2 = hcalElems_conv.begin();
00502                 ithcal2 != hcalElems_conv.end(); ++ithcal2) {
00503               
00504               if ( ! (active[ithcal2->second] ) ) continue; // HCAL Cluster already used....
00505               
00506               // TODO: Decide if this HCAL cluster is to be used
00507               // .... based on some Physics
00508               // .... To we need to check if it's closer to any other ECAL/TRACK?
00509               
00510               bool useHcal = true;
00511               if ( !useHcal ) continue;
00512               
00513               //elemsToLock.push_back(ithcal2->second);
00514 
00515             } // end of loop over HCAL clusters linked to the ECAL cluster from second CONVERSION leg
00516             
00517           } // end of loop over ECALs linked to second T_FROM_GAMMACONV
00518           
00519         } // end of loop over SECOND conversion leg
00520 
00521         // TODO: Do we need to check separatly if there are HCAL cluster linked to the track?
00522         
00523       } // end of loop over tracks
00524       
00525             
00526       // Calibrate the Added ECAL energy
00527       float addedCalibEne=0;
00528       float addedRawEne=0;
00529       std::vector<double>AddedPS1(0);
00530       std::vector<double>AddedPS2(0);  
00531       double addedps1=0;  
00532       double addedps2=0;  
00533       for(unsigned int i=0; i<AddClusters.size(); i++)  
00534         {  
00535           std::multimap<double, unsigned int> PS1Elems_conv;  
00536           std::multimap<double, unsigned int> PS2Elems_conv;  
00537           blockRef->associatedElements(AddClusters[i],  
00538                                        linkData,  
00539                                        PS1Elems_conv,  
00540                                        reco::PFBlockElement::PS1,  
00541                                        reco::PFBlock::LINKTEST_ALL );  
00542           blockRef->associatedElements( AddClusters[i],  
00543                                         linkData,  
00544                                         PS2Elems_conv,  
00545                                         reco::PFBlockElement::PS2,  
00546                                         reco::PFBlock::LINKTEST_ALL );  
00547            
00548           for(std::multimap<double, unsigned int>::iterator iteps = PS1Elems_conv.begin();  
00549               iteps != PS1Elems_conv.end(); ++iteps)  
00550             {  
00551               if(!active[iteps->second])continue;  
00552               std::multimap<double, unsigned int> PS1Elems_check;  
00553               blockRef->associatedElements(iteps->second,  
00554                                            linkData,  
00555                                            PS1Elems_check,  
00556                                            reco::PFBlockElement::ECAL,  
00557                                            reco::PFBlock::LINKTEST_ALL );  
00558               if(PS1Elems_check.begin()->second==AddClusters[i])  
00559                 {  
00560                    
00561                   reco::PFClusterRef ps1ClusterRef = elements[iteps->second].clusterRef();  
00562                   AddedPS1.push_back(ps1ClusterRef->energy());  
00563                   addedps1=addedps1+ps1ClusterRef->energy();  
00564                   elemsToLock.push_back(iteps->second);  
00565                 }  
00566             }  
00567            
00568           for(std::multimap<double, unsigned int>::iterator iteps = PS2Elems_conv.begin();  
00569               iteps != PS2Elems_conv.end(); ++iteps) {  
00570             if(!active[iteps->second])continue;  
00571             std::multimap<double, unsigned int> PS2Elems_check;  
00572             blockRef->associatedElements(iteps->second,  
00573                                          linkData,  
00574                                          PS2Elems_check,  
00575                                          reco::PFBlockElement::ECAL,  
00576                                          reco::PFBlock::LINKTEST_ALL );  
00577              
00578             if(PS2Elems_check.begin()->second==AddClusters[i])  
00579               {  
00580                 reco::PFClusterRef ps2ClusterRef = elements[iteps->second].clusterRef();  
00581                 AddedPS2.push_back(ps2ClusterRef->energy());  
00582                 addedps2=addedps2+ps2ClusterRef->energy();  
00583                 elemsToLock.push_back(iteps->second);  
00584               }  
00585           }  
00586           reco::PFClusterRef AddclusterRef = elements[AddClusters[i]].clusterRef();  
00587           addedRawEne = AddclusterRef->energy()+addedRawEne;  
00588           addedCalibEne = thePFEnergyCalibration_->energyEm(*AddclusterRef,AddedPS1,AddedPS2,false)+addedCalibEne;  
00589           AddedPS2.clear(); 
00590           AddedPS1.clear();  
00591           elemsToLock.push_back(AddClusters[i]);  
00592         }  
00593       AddClusters.clear();
00594       float EE=thePFEnergyCalibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,false)+addedCalibEne;
00595       PFClusters.push_back(&*clusterRef);
00596       if(useReg_){
00597         if(clusterRef->layer()==PFLayer::ECAL_BARREL){
00598           float LocCorr=EvaluateLCorrMVA(clusterRef);
00599           EE=LocCorr*clusterRef->energy()+addedCalibEne;
00600         }
00601         else{
00602           EE = thePFEnergyCalibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,false)+addedCalibEne; 
00603           MVALCorr.push_back(thePFEnergyCalibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,false));
00604         }       
00605       }
00606       else{
00607         if(clusterRef->layer()==PFLayer::ECAL_BARREL){
00608           float LocCorr=EvaluateLCorrMVA(clusterRef);
00609           MVALCorr.push_back(LocCorr*clusterRef->energy());
00610             }
00611         else{
00612           
00613           MVALCorr.push_back(thePFEnergyCalibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,false));         
00614         }
00615       }
00616       
00617       //cout<<"Original Energy "<<EE<<"Added Energy "<<addedCalibEne<<endl;
00618       
00619       photonEnergy_ +=  EE;
00620       RawEcalEne    +=  clusterRef->energy()+addedRawEne;
00621       photonX_      +=  EE * clusterRef->position().X();
00622       photonY_      +=  EE * clusterRef->position().Y();
00623       photonZ_      +=  EE * clusterRef->position().Z();                
00624       ps1TotEne     +=  ps1+addedps1;
00625       ps2TotEne     +=  ps2+addedps2;
00626     } // end of loop over all ECAL cluster within this SC
00627     AddFromElectron_.clear();
00628     float Elec_energy=0;
00629     float Elec_rawEcal=0;
00630     float Elec_totPs1=0;
00631     float Elec_totPs2=0;
00632     float ElectronX=0;
00633     float ElectronY=0;
00634     float ElectronZ=0;  
00635     std::vector<double>AddedPS1(0);
00636     std::vector<double>AddedPS2(0);
00637     
00638     EarlyConversion(    
00639             tempElectronCandidates,
00640             sc
00641             );   
00642     
00643     if(AddFromElectron_.size()>0)
00644       { 
00645         //collect elements from early Conversions that are reconstructed as Electrons
00646         
00647         for(std::vector<unsigned int>::const_iterator it = 
00648               AddFromElectron_.begin();
00649             it != AddFromElectron_.end(); ++it)
00650           {
00651             
00652             if(elements[*it].type()== reco::PFBlockElement::ECAL)
00653               {
00654                 //cout<<"Cluster ind "<<*it<<endl;
00655                 AddedPS1.clear();
00656                 AddedPS2.clear();
00657                 unsigned int index=*it;
00658                 reco::PFClusterRef clusterRef = 
00659                 elements[index].clusterRef();
00660                 //match to PS1 and PS2 to this cluster for calibration
00661                 Elec_rawEcal=Elec_rawEcal+
00662                   elements[index].clusterRef()->energy();
00663                 std::multimap<double, unsigned int> PS1Elems;  
00664                 std::multimap<double, unsigned int> PS2Elems;  
00665                 
00666                 blockRef->associatedElements(index,  
00667                                              linkData,  
00668                                              PS1Elems,                                               reco::PFBlockElement::PS1,  
00669                                              reco::PFBlock::LINKTEST_ALL );  
00670                 blockRef->associatedElements( index,  
00671                                               linkData,  
00672                                               PS2Elems,  
00673                                               reco::PFBlockElement::PS2,  
00674                                               reco::PFBlock::LINKTEST_ALL );
00675                 
00676                 
00677                 for(std::multimap<double, unsigned int>::iterator iteps = 
00678                       PS1Elems.begin();  
00679                     iteps != PS1Elems.end(); ++iteps) 
00680                   {  
00681                     std::multimap<double, unsigned int> Clustcheck;                         blockRef->associatedElements( iteps->second,                                                                   linkData,  
00682                                                                                                                           Clustcheck,  
00683                                                                                                                           reco::PFBlockElement::ECAL,  
00684                                                                                                                           reco::PFBlock::LINKTEST_ALL );
00685                     if(Clustcheck.begin()->second==index)
00686                       {
00687                         AddedPS1.push_back(elements[iteps->second].clusterRef()->energy());
00688                         Elec_totPs1=Elec_totPs1+elements[iteps->second].clusterRef()->energy();
00689                       }
00690                   }
00691                 
00692                 for(std::multimap<double, unsigned int>::iterator iteps = 
00693                       PS2Elems.begin();  
00694                     iteps != PS2Elems.end(); ++iteps) 
00695                   {  
00696                     std::multimap<double, unsigned int> Clustcheck;                         blockRef->associatedElements( iteps->second,                                                                   linkData,  
00697                                                                                                                           Clustcheck,  
00698                                                                                                                           reco::PFBlockElement::ECAL,  
00699                                                                                                                           reco::PFBlock::LINKTEST_ALL );
00700                     if(Clustcheck.begin()->second==index)
00701                       {
00702                         AddedPS2.push_back(elements[iteps->second].clusterRef()->energy());
00703                         Elec_totPs2=Elec_totPs2+elements[iteps->second].clusterRef()->energy();
00704                       }
00705                   }
00706                 
00707               //energy calibration 
00708                 float EE=thePFEnergyCalibration_->
00709                   energyEm(*clusterRef,AddedPS1,AddedPS2,false);
00710                 PFClusters.push_back(&*clusterRef);
00711                 if(useReg_){
00712                   if(clusterRef->layer()==PFLayer::ECAL_BARREL){
00713                     float LocCorr=EvaluateLCorrMVA(clusterRef);
00714                     EE=LocCorr*clusterRef->energy();      
00715                     MVALCorr.push_back(LocCorr*clusterRef->energy());
00716                   }
00717                   else {
00718                     EE=thePFEnergyCalibration_->energyEm(*clusterRef,AddedPS1,AddedPS2,false);
00719                     MVALCorr.push_back(EE);
00720                   }
00721                 }
00722                 else{
00723                   if(clusterRef->layer()==PFLayer::ECAL_BARREL){
00724                     float LocCorr=EvaluateLCorrMVA(clusterRef);
00725                     MVALCorr.push_back(LocCorr*clusterRef->energy());
00726                   }
00727                   else{
00728                     
00729                     MVALCorr.push_back(thePFEnergyCalibration_->energyEm(*clusterRef, AddedPS1,AddedPS2,false));          
00730                   }
00731                 }
00732                 
00733                 Elec_energy    += EE;
00734                 ElectronX      +=  EE * clusterRef->position().X();
00735                 ElectronY      +=  EE * clusterRef->position().Y();
00736                 ElectronZ      +=  EE * clusterRef->position().Z();
00737                 
00738               }
00739           }
00740         
00741       }
00742     
00743     //std::cout<<"Added Energy to Photon "<<Elec_energy<<" to "<<photonEnergy_<<std::endl;   
00744     photonEnergy_ +=  Elec_energy;
00745     RawEcalEne    +=  Elec_rawEcal;
00746     photonX_      +=  ElectronX;
00747     photonY_      +=  ElectronY;
00748     photonZ_      +=  ElectronZ;                
00749     ps1TotEne     +=  Elec_totPs1;
00750     ps2TotEne     +=  Elec_totPs2;
00751     
00752     // we've looped over all ECAL clusters, ready to generate PhotonCandidate
00753     if( ! (photonEnergy_ > 0.) ) continue;    // This SC is not a Photon Candidate
00754     float sum_track_pt=0;
00755     //Now check if there are tracks failing isolation outside of the Jurassic isolation region  
00756     for(unsigned int i=0; i<IsoTracks.size(); i++)sum_track_pt=sum_track_pt+elements[IsoTracks[i]].trackRef()->pt();  
00757     
00758 
00759 
00760     math::XYZVector photonPosition(photonX_,
00761                                    photonY_,
00762                                    photonZ_);
00763 
00764     math::XYZVector photonDirection=photonPosition.Unit();
00765     
00766     math::XYZTLorentzVector photonMomentum(photonEnergy_* photonDirection.X(),
00767                                            photonEnergy_* photonDirection.Y(),
00768                                            photonEnergy_* photonDirection.Z(),
00769                                            photonEnergy_           );
00770 
00771     if(sum_track_pt>(sumPtTrackIsoForPhoton_ + sumPtTrackIsoSlopeForPhoton_ * photonMomentum.pt()))
00772       {
00773         //cout<<"Hit Continue "<<endl;
00774         match_ind.clear(); //release the matched Electron candidates
00775         continue;
00776       }
00777 
00778         //THIS SC is not a Photon it fails track Isolation
00779     //if(sum_track_pt>(2+ 0.001* photonMomentum.pt()))
00780     //continue;//THIS SC is not a Photon it fails track Isolation
00781 
00782     /*
00783     std::cout<<" Created Photon with energy = "<<photonEnergy_<<std::endl;
00784     std::cout<<"                         pT = "<<photonMomentum.pt()<<std::endl;
00785     std::cout<<"                     RawEne = "<<RawEcalEne<<std::endl;
00786     std::cout<<"                          E = "<<photonMomentum.e()<<std::endl;
00787     std::cout<<"                        eta = "<<photonMomentum.eta()<<std::endl;
00788     std::cout<<"             TrackIsolation = "<< sum_track_pt <<std::endl;
00789     */
00790 
00791     reco::PFCandidate photonCand(0,photonMomentum, reco::PFCandidate::gamma);
00792     
00793     photonCand.setPs1Energy(ps1TotEne);
00794     photonCand.setPs2Energy(ps2TotEne);
00795     photonCand.setEcalEnergy(RawEcalEne,photonEnergy_);
00796     photonCand.setHcalEnergy(0.,0.);
00797     photonCand.set_mva_nothing_gamma(1.);  
00798     photonCand.setSuperClusterRef(sc->superClusterRef());
00799     math::XYZPoint v(primaryVertex_.x(), primaryVertex_.y(), primaryVertex_.z());
00800     photonCand.setVertex( v  );
00801     if(hasConvTrack || hasSingleleg)photonCand.setFlag( reco::PFCandidate::GAMMA_TO_GAMMACONV, true);
00802     int matches=match_ind.size();
00803     int count=0;
00804     for ( std::vector<reco::PFCandidate>::const_iterator ec=tempElectronCandidates.begin();   ec != tempElectronCandidates.end(); ++ec ){
00805       for(int i=0; i<matches; i++)
00806         {
00807           if(count==match_ind[i])photonCand.addDaughter(*ec);
00808           count++;
00809         }
00810     }
00811     //photonCand.setPositionAtECALEntrance(math::XYZPointF(photonMom_.position()));
00812     // set isvalid_ to TRUE since we've found at least one photon candidate
00813     isvalid_ = true;
00814     // push back the candidate into the collection ...
00815     //Add Elements from Electron
00816         for(std::vector<unsigned int>::const_iterator it = 
00817               AddFromElectron_.begin();
00818             it != AddFromElectron_.end(); ++it)photonCand.addElementInBlock(blockRef,*it);
00819 
00820 
00821     // ... and lock all elemts used
00822     for(std::vector<unsigned int>::const_iterator it = elemsToLock.begin();
00823         it != elemsToLock.end(); ++it)
00824       {
00825         if(active[*it])
00826           {
00827             photonCand.addElementInBlock(blockRef,*it);
00828             if( elements[*it].type() == reco::PFBlockElement::TRACK  )
00829               {
00830                 if(elements[*it].convRef().isNonnull())
00831                   {
00832                     //make sure it is not stored already as the partner track
00833                     bool matched=false;
00834                     for(unsigned int ic = 0; ic < ConversionsRef_.size(); ic++)
00835                       {
00836                         if(ConversionsRef_[ic]==elements[*it].convRef())matched=true;
00837                       }
00838                     if(!matched)ConversionsRef_.push_back(elements[*it].convRef());
00839                   }
00840               }
00841           }
00842         active[*it] = false;    
00843       }
00844     //Do Global Corrections here:
00845    float GCorr=EvaluateGCorrMVA(photonCand);
00846     if(useReg_){
00847     math::XYZTLorentzVector photonCorrMomentum(GCorr*photonEnergy_* photonDirection.X(),
00848                                                GCorr*photonEnergy_* photonDirection.Y(),
00849                                                GCorr*photonEnergy_* photonDirection.Z(),
00850                                                GCorr * photonEnergy_           );
00851     photonCand.setP4(photonCorrMomentum);
00852     }
00853     // here add the extra information
00854     PFCandidatePhotonExtra myExtra(sc->superClusterRef());
00855     //Mustache ID variables
00856     int excl=0;
00857     float Mustache_et_out=0; 
00858     Mustache Must;
00859     Must.MustacheID(PFClusters, excl, Mustache_et_out);
00860     myExtra.setMustache_Et(Mustache_et_out);
00861     myExtra.setExcludedClust(excl);
00862     //Store regressed energy
00863     myExtra.setMVAGlobalCorrE(GCorr * photonEnergy_);
00864     float Res=EvaluateResMVA(photonCand);
00865     myExtra.SetPFPhotonRes(Res*1.253);
00866     for(unsigned int l=0; l<MVALCorr.size(); ++l)
00867       {
00868         myExtra.addLCorrClusEnergy(MVALCorr[l]);
00869       }
00870 
00871     //    Daniele example for mvaValues
00872     //    do the same for single leg trackRef and convRef
00873     for(unsigned int ic = 0; ic < MVA_values.size(); ic++)
00874       {
00875         myExtra.addSingleLegConvMva(MVA_values[ic]);
00876         myExtra.addSingleLegConvTrackRef(singleLegRef[ic]);
00877         //cout<<"Single Leg Tracks "<<singleLegRef[ic]->pt()<<" MVA "<<MVA_values[ic]<<endl;
00878       }
00879     for(unsigned int ic = 0; ic < ConversionsRef_.size(); ic++)
00880       {
00881         myExtra.addConversionRef(ConversionsRef_[ic]);
00882         //cout<<"Conversion Pairs "<<ConversionsRef_[ic]->pairMomentum()<<endl;
00883       }
00884     pfPhotonExtraCandidates.push_back(myExtra);
00885     pfCandidates->push_back(photonCand);
00886     // ... and reset the vector
00887     elemsToLock.resize(0);
00888     hasConvTrack=false;
00889     hasSingleleg=false;
00890   } // end of loops over all elements in block
00891   
00892   return;
00893 }
00894 
00895 float PFPhotonAlgo::EvaluateResMVA(reco::PFCandidate photon){
00896   float BDTG=1;
00897   PFPhoEta_=photon.eta();
00898   PFPhoPhi_=photon.phi();
00899   PFPhoEt_=photon.pt();
00900   //recalculate R9 from sum PFClusterEnergy and E3x3 from Highest PFCluster Energy
00901   SCPhiWidth_=photon.superClusterRef()->phiWidth();
00902   SCEtaWidth_=photon.superClusterRef()->etaWidth();  
00903   //get from track with min R;
00904   RConv_=130;
00905   float ClustSumEt=0;
00906   std::vector<float>Clust_E(0);
00907   std::vector<float>Clust_Et(0);
00908   std::vector<float>Clust_Eta(0);
00909   std::vector<float>Clust_Phi(0);
00910   std::vector<reco::PFClusterRef> PFClusRef(0);
00911   std::vector<const CaloCluster*> PFclusters;
00912   //Multimap to sort clusters by energy
00913   std::multimap<float, int>Clust;
00914   PFCandidate::ElementsInBlocks eleInBlocks = photon.elementsInBlocks();
00915   for(unsigned i=0; i<eleInBlocks.size(); i++)
00916     {
00917       PFBlockRef blockRef = eleInBlocks[i].first;
00918       unsigned indexInBlock = eleInBlocks[i].second;
00919       const edm::OwnVector< reco::PFBlockElement >&  elements=eleInBlocks[i].first->elements();
00920       const reco::PFBlockElement& element = elements[indexInBlock];
00921       if(element.type()==reco::PFBlockElement::TRACK){
00922         float R=sqrt(element.trackRef()->innerPosition().X()*element.trackRef()->innerPosition().X()+element.trackRef()->innerPosition().Y()*element.trackRef()->innerPosition().Y());
00923         if(RConv_>R)RConv_=R;
00924       }
00925       
00926       if(element.type()==reco::PFBlockElement::ECAL){
00927         reco::PFClusterRef ClusterRef = element.clusterRef();
00928         //PFClusRef.push_back(ClusterRef);
00929         //clusters.push_back(ClusterRef);
00930         Clust_E.push_back(ClusterRef->energy());
00931         Clust_Et.push_back(ClusterRef->pt());   
00932         ClustSumEt=ClustSumEt+ClusterRef->pt();
00933         Clust_Eta.push_back(ClusterRef->eta());
00934         Clust_Phi.push_back(ClusterRef->phi());
00935         PFclusters.push_back(&*ClusterRef);
00936       }
00937       
00938       if(element.type()==reco::PFBlockElement::GSF)
00939         {
00940           //elements[indexInBlock].GsftrackRef();
00941           //  RConv_=sqrt(element.trackRef()->innerPosition().X()*element.trackRef()->innerPosition().X() + element.trackRef()->innerPosition().Y()*element.trackRef()->innerPosition().Y());
00942         }
00943       
00944     }
00945   
00946   nPFClus_=Clust_Et.size();
00947   
00948   //std::vector<int>included(0);
00949   //included=getPFMustacheClus(nPFClus_, Clust_Et, Clust_Eta, Clust_Phi);
00950   //if(included.size()>0)excluded_=nPFClus_-included.size();
00951   //else excluded_=0;
00952   //cout<<"initial excluded "<<excluded_<<endl;
00953   float Mustache_Et_out=0;
00954   int PFClus=0;
00955   Mustache Must;
00956   Must.MustacheID(PFclusters, PFClus, Mustache_Et_out);
00957   excluded_=PFClus;
00958   //order the clusters by energy
00959   //  float Mustache_Et=0;
00960   float ClusSum=0;
00961   for(unsigned int i=0; i<Clust_E.size(); ++i)
00962     {
00963       
00964       Clust.insert(make_pair(Clust_E[i], i));
00965       //Mustache_Et=Mustache_Et+Clust_Et[i];
00966       ClusSum=ClusSum+Clust_E[i];
00967     }
00968   Mustache_EtRatio_=(Mustache_Et_out)/ClustSumEt;
00969   std::multimap<float, int>::reverse_iterator it;
00970   it=Clust.rbegin();
00971   int max_c=(*it).second;
00972   it=Clust.rend();
00973   int min_c=(*it).second;
00974   if(nPFClus_>1)LowClusE_=Clust_E[min_c];
00975   else LowClusE_=0;
00976   if(nPFClus_>1){
00977     dEta_=fabs(Clust_Eta[max_c]-Clust_Eta[min_c]);
00978     dPhi_=acos(cos(Clust_Phi[max_c]-Clust_Phi[min_c]));
00979   }
00980   else{
00981     dEta_=0;
00982     dPhi_=0;
00983   }
00984   
00985   float dRmin=999;
00986   float SCphi=photon.superClusterRef()->position().phi();
00987   float SCeta=photon.superClusterRef()->position().eta();
00988   for(unsigned i=0; i<eleInBlocks.size(); i++)
00989     {
00990       PFBlockRef blockRef = eleInBlocks[i].first;
00991       unsigned indexInBlock = eleInBlocks[i].second;
00992       const edm::OwnVector< reco::PFBlockElement >&  elements=eleInBlocks[i].first->elements();
00993       const reco::PFBlockElement& element = elements[indexInBlock];
00994       if(element.type()==reco::PFBlockElement::ECAL){
00995         reco::PFClusterRef ClusterRef = element.clusterRef();
00996         float eta=ClusterRef->position().eta();
00997         float phi=ClusterRef->position().phi();
00998         float dR=deltaR(SCeta, SCphi, eta, phi);
00999         if(dR<dRmin){
01000           dRmin=dR;
01001           fill5x5Map(ClusterRef);
01002           PFPhoR9_=e3x3_/ClusSum;
01003         }
01004       }
01005     } 
01006   //fill Material Map:
01007   int ix = X0_sum->GetXaxis()->FindBin(PFPhoEta_);
01008   int iy = X0_sum->GetYaxis()->FindBin(PFPhoPhi_);
01009   x0inner_= X0_inner->GetBinContent(ix,iy);
01010   x0middle_=X0_middle->GetBinContent(ix,iy);
01011   x0outer_=X0_outer->GetBinContent(ix,iy);
01012   
01013   float GC_Var[16];
01014   GC_Var[0]=PFPhoEta_;
01015   GC_Var[1]=PFPhoEt_;
01016   GC_Var[2]=PFPhoR9_;
01017   GC_Var[3]=PFPhoPhi_;
01018   GC_Var[4]=SCEtaWidth_;
01019   GC_Var[5]=SCPhiWidth_;
01020   GC_Var[6]=x0inner_;  
01021   GC_Var[7]=x0middle_;
01022   GC_Var[8]=x0outer_;
01023   GC_Var[9]=nPFClus_;
01024   GC_Var[10]=RConv_;
01025   GC_Var[11]=LowClusE_;
01026   GC_Var[12]=dEta_;
01027   GC_Var[13]=dPhi_;
01028   GC_Var[14]=excluded_;
01029   GC_Var[15]= Mustache_EtRatio_;
01030   
01031   BDTG=ReaderRes_->GetResponse(GC_Var);
01032   //  cout<<"Res "<<BDTG<<endl;
01033   
01034   //  cout<<"BDTG Parameters X0"<<x0inner_<<", "<<x0middle_<<", "<<x0outer_<<endl;
01035   // cout<<"Et, Eta, Phi "<<PFPhoEt_<<", "<<PFPhoEta_<<", "<<PFPhoPhi_<<endl;
01036   // cout<<"PFPhoR9 "<<PFPhoR9_<<endl;
01037   // cout<<"R "<<RConv_<<endl;
01038   
01039   return BDTG;
01040 
01041 }
01042 
01043 float PFPhotonAlgo::EvaluateGCorrMVA(reco::PFCandidate photon){
01044   float BDTG=1;
01045   PFPhoEta_=photon.eta();
01046   PFPhoPhi_=photon.phi();
01047   PFPhoEt_=photon.pt();
01048   //recalculate R9 from sum PFClusterEnergy and E3x3 from Highest PFCluster Energy
01049   SCPhiWidth_=photon.superClusterRef()->phiWidth();
01050   SCEtaWidth_=photon.superClusterRef()->etaWidth();  
01051   //get from track with min R;
01052   RConv_=130;
01053   float ClustSumEt=0;
01054   std::vector<float>Clust_E(0);
01055   std::vector<float>Clust_Et(0);
01056   std::vector<float>Clust_Eta(0);
01057   std::vector<float>Clust_Phi(0);
01058   std::vector<reco::PFClusterRef> PFClusRef(0);
01059   std::vector<const CaloCluster*> PFclusters;
01060   //Multimap to sort clusters by energy
01061   std::multimap<float, int>Clust;
01062   PFCandidate::ElementsInBlocks eleInBlocks = photon.elementsInBlocks();
01063   for(unsigned i=0; i<eleInBlocks.size(); i++)
01064     {
01065       PFBlockRef blockRef = eleInBlocks[i].first;
01066       unsigned indexInBlock = eleInBlocks[i].second;
01067       const edm::OwnVector< reco::PFBlockElement >&  elements=eleInBlocks[i].first->elements();
01068       const reco::PFBlockElement& element = elements[indexInBlock];
01069       if(element.type()==reco::PFBlockElement::TRACK){
01070         float R=sqrt(element.trackRef()->innerPosition().X()*element.trackRef()->innerPosition().X()+element.trackRef()->innerPosition().Y()*element.trackRef()->innerPosition().Y());
01071         if(RConv_>R)RConv_=R;
01072       }
01073       
01074       if(element.type()==reco::PFBlockElement::ECAL){
01075         reco::PFClusterRef ClusterRef = element.clusterRef();
01076         //PFClusRef.push_back(ClusterRef);
01077         //clusters.push_back(ClusterRef);
01078         Clust_E.push_back(ClusterRef->energy());
01079         Clust_Et.push_back(ClusterRef->pt());   
01080         ClustSumEt=ClustSumEt+ClusterRef->pt();
01081         Clust_Eta.push_back(ClusterRef->eta());
01082         Clust_Phi.push_back(ClusterRef->phi());
01083         PFclusters.push_back(&*ClusterRef);
01084       }
01085       
01086       if(element.type()==reco::PFBlockElement::GSF)
01087         {
01088           //elements[indexInBlock].GsftrackRef();
01089           //  RConv_=sqrt(element.trackRef()->innerPosition().X()*element.trackRef()->innerPosition().X() + element.trackRef()->innerPosition().Y()*element.trackRef()->innerPosition().Y());
01090         }
01091       
01092     }
01093   
01094   nPFClus_=Clust_Et.size();
01095   
01096   //std::vector<int>included(0);
01097   //included=getPFMustacheClus(nPFClus_, Clust_Et, Clust_Eta, Clust_Phi);
01098   //if(included.size()>0)excluded_=nPFClus_-included.size();
01099   //else excluded_=0;
01100   //cout<<"initial excluded "<<excluded_<<endl;
01101   float Mustache_Et_out=0;
01102   int PFClus=0;
01103   Mustache Must;
01104   Must.MustacheID(PFclusters, PFClus, Mustache_Et_out);
01105   excluded_=PFClus;
01106   //order the clusters by energy
01107   //  float Mustache_Et=0;
01108   float ClusSum=0;
01109   for(unsigned int i=0; i<Clust_E.size(); ++i)
01110     {
01111       
01112       Clust.insert(make_pair(Clust_E[i], i));
01113       //Mustache_Et=Mustache_Et+Clust_Et[i];
01114       ClusSum=ClusSum+Clust_E[i];
01115     }
01116   Mustache_EtRatio_=(Mustache_Et_out)/ClustSumEt;
01117   std::multimap<float, int>::reverse_iterator it;
01118   it=Clust.rbegin();
01119   int max_c=(*it).second;
01120   it=Clust.rend();
01121   int min_c=(*it).second;
01122   if(nPFClus_>1)LowClusE_=Clust_E[min_c];
01123   else LowClusE_=0;
01124   if(nPFClus_>1){
01125     dEta_=fabs(Clust_Eta[max_c]-Clust_Eta[min_c]);
01126     dPhi_=acos(cos(Clust_Phi[max_c]-Clust_Phi[min_c]));
01127   }
01128   else{
01129     dEta_=0;
01130     dPhi_=0;
01131   }
01132   float dRmin=999;
01133   float SCphi=photon.superClusterRef()->position().phi();
01134   float SCeta=photon.superClusterRef()->position().eta();
01135   for(unsigned i=0; i<eleInBlocks.size(); i++)
01136     {
01137       PFBlockRef blockRef = eleInBlocks[i].first;
01138       unsigned indexInBlock = eleInBlocks[i].second;
01139       const edm::OwnVector< reco::PFBlockElement >&  elements=eleInBlocks[i].first->elements();
01140       const reco::PFBlockElement& element = elements[indexInBlock];
01141       if(element.type()==reco::PFBlockElement::ECAL){
01142         reco::PFClusterRef ClusterRef = element.clusterRef();
01143         float eta=ClusterRef->position().eta();
01144         float phi=ClusterRef->position().phi();
01145         float dR=deltaR(SCeta, SCphi, eta, phi);
01146         if(dR<dRmin){
01147           dRmin=dR;
01148           fill5x5Map(ClusterRef);
01149           PFPhoR9_=e3x3_/ClusSum;
01150         }
01151       }
01152     } 
01153 
01154   //fill Material Map:
01155   int ix = X0_sum->GetXaxis()->FindBin(PFPhoEta_);
01156   int iy = X0_sum->GetYaxis()->FindBin(PFPhoPhi_);
01157   x0inner_= X0_inner->GetBinContent(ix,iy);
01158   x0middle_=X0_middle->GetBinContent(ix,iy);
01159   x0outer_=X0_outer->GetBinContent(ix,iy);
01160 
01161   float GC_Var[16];
01162   GC_Var[0]=PFPhoEta_;
01163   GC_Var[1]=PFPhoEt_;
01164   GC_Var[2]=PFPhoR9_;
01165   GC_Var[3]=PFPhoPhi_;
01166   GC_Var[4]=SCEtaWidth_;
01167   GC_Var[5]=SCPhiWidth_;
01168   GC_Var[6]=x0inner_;  
01169   GC_Var[7]=x0middle_;
01170   GC_Var[8]=x0outer_;
01171   GC_Var[9]=nPFClus_;
01172   GC_Var[10]=RConv_;
01173   GC_Var[11]=LowClusE_/photon.energy();
01174   GC_Var[12]=dEta_;
01175   GC_Var[13]=dPhi_;
01176   GC_Var[14]=excluded_;
01177   GC_Var[15]= Mustache_EtRatio_;
01178   
01179 
01180   BDTG=ReaderGC_->GetResponse(GC_Var);
01181   //  cout<<"GC "<<BDTG<<endl;  
01182   
01183   //  cout<<"BDTG Parameters X0"<<x0inner_<<", "<<x0middle_<<", "<<x0outer_<<endl;
01184   // cout<<"Et, Eta, Phi "<<PFPhoEt_<<", "<<PFPhoEta_<<", "<<PFPhoPhi_<<endl;
01185   // cout<<"PFPhoR9 "<<PFPhoR9_<<endl;
01186   // cout<<"R "<<RConv_<<endl;
01187   
01188   return BDTG;
01189 
01190 }
01191 
01192 float PFPhotonAlgo::EvaluateLCorrMVA(reco::PFClusterRef clusterRef ){
01193   float BDTG=1;
01194   
01195   GetCrysCoordinates(clusterRef);
01196   fill5x5Map(clusterRef);
01197   VtxZ_=primaryVertex_.z();
01198   ClusPhi_=clusterRef->position().phi(); 
01199   ClusEta_=fabs(clusterRef->position().eta());
01200   EB=fabs(clusterRef->position().eta())/clusterRef->position().eta();
01201   logPFClusE_=log(clusterRef->energy());
01202    float LC_Var[20];
01203    LC_Var[0]=VtxZ_;
01204    LC_Var[1]=EB;
01205    LC_Var[2]=ClusEta_;
01206    LC_Var[3]=ClusPhi_;
01207    LC_Var[4]=logPFClusE_;
01208    LC_Var[5]=eSeed_;
01209    LC_Var[6]=ClusR9_;
01210    LC_Var[7]=e1x3_;
01211    LC_Var[8]=e3x1_;
01212    LC_Var[9]=Clus5x5ratio_;
01213    LC_Var[10]=e1x5_;
01214    LC_Var[11]=e2x5Max_;
01215    LC_Var[12]=e2x5Top_;
01216    LC_Var[13]=e2x5Bottom_;
01217    LC_Var[14]=e2x5Left_;
01218    LC_Var[15]=e2x5Right_;
01219    LC_Var[16]=CrysEta_;
01220    LC_Var[17]=CrysPhi_;
01221    LC_Var[18]=PFCrysPhiCrack_;
01222    LC_Var[19]=PFCrysEtaCrack_;
01223    BDTG=ReaderLC_->GetResponse(LC_Var);   
01224    //   cout<<"LC "<<BDTG<<endl;  
01225    return BDTG;
01226   
01227 }
01228 
01229 
01230 void PFPhotonAlgo::GetCrysCoordinates(reco::PFClusterRef clusterRef){
01231   // Commented by Florian. Source of warnings as not used 
01232   //  float PFSeedEta=99;
01233   float PFSeedPhi=99;
01234   float PFSeedTheta=99;
01235   double PFSeedE=0;
01236   //  unsigned int SeedDetId=-1;
01237   // float seedPhi=0;
01238   //  float seedEta=0;
01239   DetId idseed;
01240   const std::vector< reco::PFRecHitFraction >& PFRecHits=
01241     clusterRef->recHitFractions();
01242   for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin();
01243         it != PFRecHits.end(); ++it){
01244     const PFRecHitRef& RefPFRecHit = it->recHitRef();
01245     unsigned index=it-PFRecHits.begin();
01246     float frac=clusterRef->hitsAndFractions()[index].second;
01247     float E= RefPFRecHit->energy()* frac;
01248     if(E>PFSeedE){
01249       //      SeedDetId=RefPFRecHit.index();
01250       PFSeedE=E;  
01251       //      PFSeedEta=RefPFRecHit->positionREP().eta(); 
01252       PFSeedPhi=RefPFRecHit->positionREP().phi();
01253       PFSeedTheta=RefPFRecHit->positionREP().theta();
01254       RefPFRecHit->positionREP().theta();
01255       idseed = RefPFRecHit->detId();
01256     }
01257   }
01258   EBDetId EBidSeed=EBDetId(idseed.rawId());
01259   CrysIEta_=EBidSeed.ieta();
01260   CrysIPhi_=EBidSeed.iphi();
01261   
01262   //Crystal Coordinates:
01263   double Pi=TMath::Pi();
01264   float Phi=clusterRef->position().phi(); 
01265   //  float Eta=clusterRef->position().eta();
01266   double Theta = -(clusterRef->position().theta())+0.5* Pi;
01267   double PhiCentr = TVector2::Phi_mpi_pi(PFSeedPhi);
01268   double PhiWidth = (Pi/180.);
01269   double PhiCry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
01270   double ThetaCentr = -PFSeedTheta+0.5*Pi;
01271   double ThetaWidth = (Pi/180.)*cos(ThetaCentr);
01272   
01273   //cout<<"Clust Theta "<<Theta<<" Crys Theta "<<ThetaCentr<<endl;
01274   //cout<<" Width "<< ThetaWidth<<endl;
01275   double EtaCry = (Theta-ThetaCentr)/ThetaWidth; 
01276   CrysEta_=EtaCry;
01277   CrysPhi_=PhiCry;
01278  
01279   //check Module and crack:
01280   int iphi=CrysIPhi_;
01281   int phimod=iphi%20;
01282   if(phimod>1)PFCrysPhiCrack_=2;
01283   else PFCrysPhiCrack_=phimod; //should be 0, 1
01284   
01285   if(abs(CrysIEta_)==1 || abs(CrysIEta_)==2 )
01286     PFCrysEtaCrack_=abs(CrysIEta_);
01287   if(abs(CrysIEta_)>2 && abs(CrysIEta_)<24)
01288     PFCrysEtaCrack_=3;
01289   if(abs(CrysIEta_)==24)
01290     PFCrysEtaCrack_=4;
01291   if(abs(CrysIEta_)==25)
01292         PFCrysEtaCrack_=5;
01293       if(abs(CrysIEta_)==26)
01294         PFCrysEtaCrack_=6;
01295       if(abs(CrysIEta_)==27)
01296                 PFCrysEtaCrack_=7;
01297       if(abs(CrysIEta_)>27 &&  abs(CrysIEta_)<44)
01298         PFCrysEtaCrack_=8;
01299       if(abs(CrysIEta_)==44)
01300                 PFCrysEtaCrack_=9;
01301       if(abs(CrysIEta_)==45)
01302         PFCrysEtaCrack_=10;
01303       if(abs(CrysIEta_)==46)
01304                 PFCrysEtaCrack_=11;
01305       if(abs(CrysIEta_)==47)
01306         PFCrysEtaCrack_=12;
01307       if(abs(CrysIEta_)>47 &&  abs(CrysIEta_)<64)
01308         PFCrysEtaCrack_=13;
01309       if(abs(CrysIEta_)==64)
01310         PFCrysEtaCrack_=14;
01311       if(abs(CrysIEta_)==65)
01312         PFCrysEtaCrack_=15;
01313       if(abs(CrysIEta_)==66)
01314         PFCrysEtaCrack_=16;
01315       if(abs(CrysIEta_)==67)
01316         PFCrysEtaCrack_=17;
01317       if(abs(CrysIEta_)>67 &&  abs(CrysIEta_)<84)
01318         PFCrysEtaCrack_=18;
01319       if(abs(CrysIEta_)==84)
01320         PFCrysEtaCrack_=19;
01321       if(abs(CrysIEta_)==85)
01322         PFCrysEtaCrack_=20;
01323       
01324 }
01325 
01326 void PFPhotonAlgo::fill5x5Map(reco::PFClusterRef clusterRef){
01327 
01328   // commented out by Florian. Not used in the code
01329   //  float PFSeedEta=99;
01330   //  float PFSeedPhi=99;
01331   double PFSeedE=0;
01332   //  unsigned int SeedDetId=-1;
01333   //  float seedPhi=0;
01334   //float seedEta=0;
01335   DetId idseed;
01336   const std::vector< reco::PFRecHitFraction >& PFRecHits=
01337     clusterRef->recHitFractions();
01338   for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin();
01339         it != PFRecHits.end(); ++it){
01340     const PFRecHitRef& RefPFRecHit = it->recHitRef();
01341     unsigned index=it-PFRecHits.begin();
01342     //DetId id=eclusterRef->hitsAndFractions()[index].first;
01343     float frac=clusterRef->hitsAndFractions()[index].second;
01344     float E= RefPFRecHit->energy()* frac;
01345     if(E>PFSeedE){
01346       //      SeedDetId=RefPFRecHit.index();
01347       PFSeedE=E;  
01348       //      PFSeedEta=RefPFRecHit->positionREP().eta(); 
01349       //      PFSeedPhi=RefPFRecHit->positionREP().phi(); 
01350       idseed = RefPFRecHit->detId();
01351     }
01352   }
01353   
01354   
01355   //initialize 5x5 map
01356   for(int i=0; i<5; ++i)
01357     for(int j=0; j<5; ++j)e5x5Map[i][j]=0;
01358   float E3x3=0;
01359   float E5x5=0;
01360   //  int count=0;
01361   for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin();
01362         it != PFRecHits.end(); ++it){
01363     unsigned index=it-PFRecHits.begin();
01364     const PFRecHitRef& RefPFRecHit = it->recHitRef();
01365     float frac=clusterRef->hitsAndFractions()[index].second;
01366     DetId id = RefPFRecHit->detId();
01367     if(idseed.subdetId()==EcalBarrel){
01368       int deta=EBDetId::distanceEta(id,idseed);
01369       int dphi=EBDetId::distancePhi(id,idseed); 
01370       EBDetId EBidSeed=EBDetId(idseed.rawId());
01371       if(abs(dphi)<=1 && abs(deta)<=1)
01372         {
01373           E3x3=E3x3+(RefPFRecHit->energy()*frac);
01374         }
01375       if(abs(dphi)<=2 && abs(deta)<=2)
01376         {
01377           //center the array on [2][2] to be the Seed
01378           //deta and deta are unsigned so you want to recompute 
01379           //to get top bottom left right 
01380           EBDetId EBid=EBDetId(id.rawId());
01381           //note this is actually opposite to make it consistent to the lazy tools left right which inverts them
01382           int i=EBidSeed.ieta()-EBid.ieta();
01383           
01384           int j=EBid.iphi()-EBidSeed.iphi();
01385           int iEta=i+2;
01386           int iPhi=j+2;
01387           
01388           e5x5Map[iEta][iPhi]=RefPFRecHit->energy()*frac;
01389           E5x5=E5x5+(RefPFRecHit->energy()*frac);
01390         }
01391     }
01392     if(idseed.subdetId()==EcalEndcap){
01393       //dx and dy are unsigned so you want to recompute 
01394       //to get top bottom left right 
01395       int dx=EEDetId::distanceX(id,idseed);
01396       int dy=EEDetId::distanceY(id,idseed);
01397       EEDetId EEidSeed=EEDetId(idseed.rawId());
01398       //dx and dy are unsigned so you want to recompute 
01399       //to get top bottom left right 
01400       if(abs(dx)<=1 && abs(dy)<=1)
01401         {
01402           E3x3=E3x3+(RefPFRecHit->energy()*frac);
01403         }
01404       if(abs(dx)<=2 && abs(dy)<=2)
01405         {
01406           EEDetId EEid=EEDetId(id.rawId());
01407           int i=EEid.ix()-EEidSeed.ix();
01408           int j=EEid.iy()-EEidSeed.iy();
01409           //center the array on [2][2] to be the Seed
01410           int ix=i+2;
01411           int iy=j+2;
01412           e5x5Map[ix][iy]=RefPFRecHit->energy()*frac;
01413           E5x5=E5x5+(RefPFRecHit->energy()*frac);
01414         }
01415     }
01416     
01417   }
01418   e3x3_=E3x3;
01419   ClusR9_=E3x3/clusterRef->energy();
01420   Clus5x5ratio_=E5x5/clusterRef->energy();
01421   eSeed_= e5x5Map[2][2]/clusterRef->energy();
01422   e1x3_=(e5x5Map[2][2]+e5x5Map[1][2]+e5x5Map[3][2])/clusterRef->energy();
01423   e3x1_=(e5x5Map[2][2]+e5x5Map[2][1]+e5x5Map[2][3])/clusterRef->energy();
01424   e1x5_=e5x5Map[2][2]+e5x5Map[2][0]+e5x5Map[2][1]+e5x5Map[2][3]+e5x5Map[2][4];
01425   e2x5Top_=(e5x5Map[0][4]+e5x5Map[1][4]+e5x5Map[2][4]
01426             +e5x5Map[3][4]+e5x5Map[4][4]
01427     +e5x5Map[0][3]+e5x5Map[1][3]+e5x5Map[2][3]
01428             +e5x5Map[3][3]+e5x5Map[4][3])/clusterRef->energy();
01429   
01430   e2x5Bottom_=(e5x5Map[0][0]+e5x5Map[1][0]+e5x5Map[2][0]
01431                +e5x5Map[3][0]+e5x5Map[4][0]
01432                +e5x5Map[0][1]+e5x5Map[1][1]
01433                +e5x5Map[2][1]+e5x5Map[3][1]+e5x5Map[4][1])/clusterRef->energy();
01434   e2x5Left_= ( e5x5Map[0][1]+e5x5Map[0][1]
01435                +e5x5Map[0][2]+e5x5Map[0][3]+e5x5Map[0][4]
01436                +e5x5Map[1][0]+e5x5Map[1][1]+e5x5Map[1][2]
01437                +e5x5Map[1][3]+e5x5Map[1][4])/clusterRef->energy();
01438   
01439   e2x5Right_ =(e5x5Map[4][0]+e5x5Map[4][1]
01440                +e5x5Map[4][2]+e5x5Map[4][3]+e5x5Map[4][4]
01441                +e5x5Map[3][0]+e5x5Map[3][1]+e5x5Map[3][2]
01442                +e5x5Map[3][3]+e5x5Map[3][4])/clusterRef->energy(); 
01443   float centerstrip=e5x5Map[2][2]+e5x5Map[2][0]
01444     +e5x5Map[2][1]+e5x5Map[2][3]+e5x5Map[2][4];
01445   float rightstrip=e5x5Map[3][2]+e5x5Map[3][0]
01446     +e5x5Map[3][1]+e5x5Map[3][3]+e5x5Map[3][4];
01447   float leftstrip=e5x5Map[1][2]+e5x5Map[1][0]+e5x5Map[1][1]
01448     +e5x5Map[1][3]+e5x5Map[1][4];
01449   if(rightstrip>leftstrip)e2x5Max_=rightstrip+centerstrip;
01450   else e2x5Max_=leftstrip+centerstrip;
01451   e2x5Max_=e2x5Max_/clusterRef->energy();
01452 }
01453 
01454 
01455 
01456 
01457 bool PFPhotonAlgo::EvaluateSingleLegMVA(const reco::PFBlockRef& blockref, const reco::Vertex& primaryvtx, unsigned int track_index)  
01458 {  
01459   bool convtkfound=false;  
01460   const reco::PFBlock& block = *blockref;  
01461   const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();  
01462   //use this to store linkdata in the associatedElements function below  
01463   PFBlock::LinkData linkData =  block.linkData();  
01464   //calculate MVA Variables  
01465   chi2=elements[track_index].trackRef()->chi2()/elements[track_index].trackRef()->ndof();  
01466   nlost=elements[track_index].trackRef()->trackerExpectedHitsInner().numberOfLostHits();  
01467   nlayers=elements[track_index].trackRef()->hitPattern().trackerLayersWithMeasurement();  
01468   track_pt=elements[track_index].trackRef()->pt();  
01469   STIP=elements[track_index].trackRefPF()->STIP();  
01470    
01471   float linked_e=0;  
01472   float linked_h=0;  
01473   std::multimap<double, unsigned int> ecalAssoTrack;  
01474   block.associatedElements( track_index,linkData,  
01475                             ecalAssoTrack,  
01476                             reco::PFBlockElement::ECAL,  
01477                             reco::PFBlock::LINKTEST_ALL );  
01478   std::multimap<double, unsigned int> hcalAssoTrack;  
01479   block.associatedElements( track_index,linkData,  
01480                             hcalAssoTrack,  
01481                             reco::PFBlockElement::HCAL,  
01482                             reco::PFBlock::LINKTEST_ALL );  
01483   if(ecalAssoTrack.size() > 0) {  
01484     for(std::multimap<double, unsigned int>::iterator itecal = ecalAssoTrack.begin();  
01485         itecal != ecalAssoTrack.end(); ++itecal) {  
01486       linked_e=linked_e+elements[itecal->second].clusterRef()->energy();  
01487     }  
01488   }  
01489   if(hcalAssoTrack.size() > 0) {  
01490     for(std::multimap<double, unsigned int>::iterator ithcal = hcalAssoTrack.begin();  
01491         ithcal != hcalAssoTrack.end(); ++ithcal) {  
01492       linked_h=linked_h+elements[ithcal->second].clusterRef()->energy();  
01493     }  
01494   }  
01495   EoverPt=linked_e/elements[track_index].trackRef()->pt();  
01496   HoverPt=linked_h/elements[track_index].trackRef()->pt();  
01497   GlobalVector rvtx(elements[track_index].trackRef()->innerPosition().X()-primaryvtx.x(),  
01498                     elements[track_index].trackRef()->innerPosition().Y()-primaryvtx.y(),  
01499                     elements[track_index].trackRef()->innerPosition().Z()-primaryvtx.z());  
01500   double vtx_phi=rvtx.phi();  
01501   //delta Phi between conversion vertex and track  
01502   del_phi=fabs(deltaPhi(vtx_phi, elements[track_index].trackRef()->innerMomentum().Phi()));  
01503   mvaValue = tmvaReader_->EvaluateMVA("BDT");  
01504   if(mvaValue > MVACUT)convtkfound=true;  
01505   return convtkfound;  
01506 }
01507 
01508 //Recover Early Conversions reconstructed as PFelectrons
01509 void PFPhotonAlgo::EarlyConversion(    
01510                                    //std::auto_ptr< reco::PFCandidateCollection > 
01511                                    //&pfElectronCandidates_,
01512                                    std::vector<reco::PFCandidate>& 
01513                                    tempElectronCandidates,
01514                                    const reco::PFBlockElementSuperCluster* sc
01515                                    ){
01516   //step 1 check temp electrons for clusters that match Photon Supercluster:
01517   // permElectronCandidates->clear();
01518   int count=0;
01519   for ( std::vector<reco::PFCandidate>::const_iterator ec=tempElectronCandidates.begin();   ec != tempElectronCandidates.end(); ++ec ) 
01520     {
01521       //      bool matched=false;
01522       int mh=ec->gsfTrackRef()->trackerExpectedHitsInner().numberOfLostHits();
01523       if(mh==0)continue;//Case where missing hits greater than zero
01524       
01525       reco::GsfTrackRef gsf=ec->gsfTrackRef();
01526       //some hoopla to get Electron SC ref
01527       
01528       if(gsf->extra().isAvailable() && gsf->extra()->seedRef().isAvailable()) 
01529         {
01530           reco::ElectronSeedRef seedRef=  gsf->extra()->seedRef().castTo<reco::ElectronSeedRef>();
01531           if(seedRef.isAvailable() && seedRef->isEcalDriven()) 
01532             {
01533               reco::SuperClusterRef ElecscRef = seedRef->caloCluster().castTo<reco::SuperClusterRef>();
01534               
01535               if(ElecscRef.isNonnull()){
01536                 //finally see if it matches:
01537                 reco::SuperClusterRef PhotscRef=sc->superClusterRef();
01538                 
01539                 if(PhotscRef==ElecscRef)
01540                   {
01541                     match_ind.push_back(count);
01542                     //  matched=true; 
01543                     //cout<<"Matched Electron with Index "<<count<<" This is the electron "<<*ec<<endl;
01544                     //find that they have the same SC footprint start to collect Clusters and tracks and these will be passed to PFPhoton
01545                     reco::PFCandidate::ElementsInBlocks eleInBlocks = ec->elementsInBlocks();
01546                     for(unsigned i=0; i<eleInBlocks.size(); i++) 
01547                       {
01548                         reco::PFBlockRef blockRef = eleInBlocks[i].first;
01549                         unsigned indexInBlock = eleInBlocks[i].second;   
01550                         //const edm::OwnVector< reco::PFBlockElement >&  elements=eleInBlocks[i].first->elements();
01551                         //const reco::PFBlockElement& element = elements[indexInBlock];                 
01552                         
01553                         AddFromElectron_.push_back(indexInBlock);               
01554                       }             
01555                   }             
01556               }
01557             }     
01558         }           
01559       count++;
01560     }
01561 }