Reference Manual

Go to the documentation of this file.
00001 #include "RecoTauTag/RecoTau/interface/PFRecoTauAlgorithm.h"
00002 
00003 // Turn off filtering by pt.  Min pt is set to zero,
00004 //  as this functionality is implemented in the underlying
00005 //  PFTauTagInfo production.  Additional pt filters are applied
00006 //  the discriminators.
00007 
00008 #define PFTauAlgo_NeutrHadrCand_minPt_   (0.0)
00009 #define PFTauAlgo_GammaCand_minPt_       (0.0)
00010 #define PFTauAlgo_PFCand_minPt_          (0.0)
00011 #define PFTauAlgo_Track_minPt_           (0.0)
00012 #define PFTauAlgo_ChargedHadrCand_minPt_ (0.0)
00013 
00014 using namespace reco;
00015 
00016 PFRecoTauAlgorithm::PFRecoTauAlgorithm() : PFRecoTauAlgorithmBase(){}  
00017 PFRecoTauAlgorithm::PFRecoTauAlgorithm(const edm::ParameterSet& iConfig):PFRecoTauAlgorithmBase(iConfig){
00018    LeadPFCand_minPt_                   = iConfig.getParameter<double>("LeadPFCand_minPt"); 
00019 
00020    UseChargedHadrCandLeadChargedHadrCand_tksDZconstraint_ 
00021       = iConfig.getParameter<bool>("UseChargedHadrCandLeadChargedHadrCand_tksDZconstraint");
00022 
00023    ChargedHadrCandLeadChargedHadrCand_tksmaxDZ_ 
00024       = iConfig.getParameter<double>("ChargedHadrCandLeadChargedHadrCand_tksmaxDZ");
00025 
00026    LeadTrack_minPt_                    = iConfig.getParameter<double>("LeadTrack_minPt");
00027    UseTrackLeadTrackDZconstraint_      = iConfig.getParameter<bool>("UseTrackLeadTrackDZconstraint");
00028    TrackLeadTrack_maxDZ_               = iConfig.getParameter<double>("TrackLeadTrack_maxDZ");
00029 
00030    MatchingConeMetric_                 = iConfig.getParameter<std::string>("MatchingConeMetric");
00031    MatchingConeSizeFormula_            = iConfig.getParameter<std::string>("MatchingConeSizeFormula");
00032    MatchingConeSize_min_               = iConfig.getParameter<double>("MatchingConeSize_min");
00033    MatchingConeSize_max_               = iConfig.getParameter<double>("MatchingConeSize_max");
00034    TrackerSignalConeMetric_            = iConfig.getParameter<std::string>("TrackerSignalConeMetric");
00035    TrackerSignalConeSizeFormula_       = iConfig.getParameter<std::string>("TrackerSignalConeSizeFormula");
00036    TrackerSignalConeSize_min_          = iConfig.getParameter<double>("TrackerSignalConeSize_min");
00037    TrackerSignalConeSize_max_          = iConfig.getParameter<double>("TrackerSignalConeSize_max");
00038    TrackerIsolConeMetric_              = iConfig.getParameter<std::string>("TrackerIsolConeMetric"); 
00039    TrackerIsolConeSizeFormula_         = iConfig.getParameter<std::string>("TrackerIsolConeSizeFormula"); 
00040    TrackerIsolConeSize_min_            = iConfig.getParameter<double>("TrackerIsolConeSize_min");
00041    TrackerIsolConeSize_max_            = iConfig.getParameter<double>("TrackerIsolConeSize_max");
00042    ECALSignalConeMetric_               = iConfig.getParameter<std::string>("ECALSignalConeMetric");
00043    ECALSignalConeSizeFormula_          = iConfig.getParameter<std::string>("ECALSignalConeSizeFormula");    
00044    ECALSignalConeSize_min_             = iConfig.getParameter<double>("ECALSignalConeSize_min");
00045    ECALSignalConeSize_max_             = iConfig.getParameter<double>("ECALSignalConeSize_max");
00046    ECALIsolConeMetric_                 = iConfig.getParameter<std::string>("ECALIsolConeMetric");
00047    ECALIsolConeSizeFormula_            = iConfig.getParameter<std::string>("ECALIsolConeSizeFormula");      
00048    ECALIsolConeSize_min_               = iConfig.getParameter<double>("ECALIsolConeSize_min");
00049    ECALIsolConeSize_max_               = iConfig.getParameter<double>("ECALIsolConeSize_max");
00050    HCALSignalConeMetric_               = iConfig.getParameter<std::string>("HCALSignalConeMetric");
00051    HCALSignalConeSizeFormula_          = iConfig.getParameter<std::string>("HCALSignalConeSizeFormula");    
00052    HCALSignalConeSize_min_             = iConfig.getParameter<double>("HCALSignalConeSize_min");
00053    HCALSignalConeSize_max_             = iConfig.getParameter<double>("HCALSignalConeSize_max");
00054    HCALIsolConeMetric_                 = iConfig.getParameter<std::string>("HCALIsolConeMetric");
00055    HCALIsolConeSizeFormula_            = iConfig.getParameter<std::string>("HCALIsolConeSizeFormula");      
00056    HCALIsolConeSize_min_               = iConfig.getParameter<double>("HCALIsolConeSize_min");
00057    HCALIsolConeSize_max_               = iConfig.getParameter<double>("HCALIsolConeSize_max");
00058 
00059    // get paramaeters for ellipse EELL
00060    Rphi_                                      = iConfig.getParameter<double>("Rphi");
00061    MaxEtInEllipse_                    = iConfig.getParameter<double>("MaxEtInEllipse");
00062    AddEllipseGammas_                  = iConfig.getParameter<bool>("AddEllipseGammas");
00063    // EELL
00064 
00065    AreaMetric_recoElements_maxabsEta_    = iConfig.getParameter<double>("AreaMetric_recoElements_maxabsEta");
00066    ChargedHadrCand_IsolAnnulus_minNhits_ = iConfig.getParameter<uint32_t>("ChargedHadrCand_IsolAnnulus_minNhits");
00067    Track_IsolAnnulus_minNhits_           = iConfig.getParameter<uint32_t>("Track_IsolAnnulus_minNhits");
00068 
00069    ElecPreIDLeadTkMatch_maxDR_           = iConfig.getParameter<double>("ElecPreIDLeadTkMatch_maxDR");
00070    EcalStripSumE_minClusEnergy_          = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
00071    EcalStripSumE_deltaEta_               = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
00072    EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
00073    EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
00074    maximumForElectrionPreIDOutput_       = iConfig.getParameter<double>("maximumForElectrionPreIDOutput");
00075 
00076    DataType_ = iConfig.getParameter<std::string>("DataType");
00077 
00078    //TFormula computation
00079    myMatchingConeSizeTFormula      = TauTagTools::computeConeSizeTFormula(MatchingConeSizeFormula_,"Matching cone size");
00080    //Charged particles cones
00081    myTrackerSignalConeSizeTFormula = TauTagTools::computeConeSizeTFormula(TrackerSignalConeSizeFormula_,"Tracker signal cone size");
00082    myTrackerIsolConeSizeTFormula   = TauTagTools::computeConeSizeTFormula(TrackerIsolConeSizeFormula_,"Tracker isolation cone size");
00083    //Gamma candidates cones
00084    myECALSignalConeSizeTFormula    = TauTagTools::computeConeSizeTFormula(ECALSignalConeSizeFormula_,"ECAL signal cone size");
00085    myECALIsolConeSizeTFormula      = TauTagTools::computeConeSizeTFormula(ECALIsolConeSizeFormula_,"ECAL isolation cone size");
00086    //Neutral hadrons cones
00087    myHCALSignalConeSizeTFormula    = TauTagTools::computeConeSizeTFormula(HCALSignalConeSizeFormula_,"HCAL signal cone size");
00088    myHCALIsolConeSizeTFormula      = TauTagTools::computeConeSizeTFormula(HCALIsolConeSizeFormula_,"HCAL isolation cone size");
00089 }
00090 
00091 
00092 PFTau PFRecoTauAlgorithm::buildPFTau(const PFTauTagInfoRef& myPFTauTagInfoRef, const Vertex& myPV)
00093 {
00094    PFJetRef myPFJet=(*myPFTauTagInfoRef).pfjetRef();  // catch a ref to the initial PFJet  
00095    PFTau myPFTau(std::numeric_limits<int>::quiet_NaN(),myPFJet->p4());   // create the PFTau
00096 
00097    myPFTau.setpfTauTagInfoRef(myPFTauTagInfoRef);
00098 
00099    PFCandidateRefVector myPFCands=(*myPFTauTagInfoRef).PFCands();
00100 
00101    PFTauElementsOperators myPFTauElementsOperators(myPFTau);
00102    double myMatchingConeSize=myPFTauElementsOperators.computeConeSize(myMatchingConeSizeTFormula,MatchingConeSize_min_,MatchingConeSize_max_);
00103 
00104    PFCandidateRef myleadPFChargedCand=myPFTauElementsOperators.leadPFChargedHadrCand(MatchingConeMetric_,myMatchingConeSize,PFTauAlgo_PFCand_minPt_);
00105 
00106    // These two quantities always taken from the signal cone
00107    PFCandidateRef myleadPFNeutralCand;
00108    PFCandidateRef myleadPFCand;
00109 
00110    bool myleadPFCand_rectkavailable = false;
00111    double myleadPFCand_rectkDZ      = 0.;
00112 
00113    // Determine the SIPT of the lead track
00114    if(myleadPFChargedCand.isNonnull()) {
00115       myPFTau.setleadPFChargedHadrCand(myleadPFChargedCand);
00116       TrackRef myleadPFCand_rectk=(*myleadPFChargedCand).trackRef();
00117       if(myleadPFCand_rectk.isNonnull()) {
00118          myleadPFCand_rectkavailable=true;
00119          myleadPFCand_rectkDZ=(*myleadPFCand_rectk).dz(myPV.position());
00120          if(TransientTrackBuilder_!=0) { 
00121             const TransientTrack myleadPFCand_rectransienttk=TransientTrackBuilder_->build(&(*myleadPFCand_rectk));
00122             GlobalVector myPFJetdir((*myPFJet).px(),(*myPFJet).py(),(*myPFJet).pz());
00123             if(IPTools::signedTransverseImpactParameter(myleadPFCand_rectransienttk,myPFJetdir,myPV).first)
00124                myPFTau.setleadPFChargedHadrCandsignedSipt(
00125                      IPTools::signedTransverseImpactParameter(myleadPFCand_rectransienttk,myPFJetdir,myPV).second.significance());
00126          }
00127       }
00128    }
00129 
00130    //Building PF Components
00131    if (myleadPFChargedCand.isNonnull())
00132    {
00133       math::XYZVector tauAxis = myleadPFChargedCand->momentum();
00134       // Compute energy of the PFTau considering only inner constituents 
00135       // (inner == pfcandidates inside a cone which is equal to the maximum value of the signal cone)
00136       // The axis is built about the lead charged hadron
00137       PFCandidateRefVector myTmpPFCandsInSignalCone = 
00138          myPFTauElementsOperators.PFCandsInCone(tauAxis,TrackerSignalConeMetric_,TrackerSignalConeSize_max_,0.5);
00139       math::XYZTLorentzVector tmpLorentzVect(0.,0.,0.,0.);
00140 
00141       double jetOpeningAngle = 0.0;
00142       for (PFCandidateRefVector::const_iterator iCand = myTmpPFCandsInSignalCone.begin(); 
00143             iCand != myTmpPFCandsInSignalCone.end(); iCand++) 
00144       {
00145          //find the maximum opening angle of the jet (now a parameter in available TFormulas)
00146          double deltaRToSeed = TauTagTools::computeDeltaR(tauAxis, (**iCand).momentum());
00147          if (deltaRToSeed > jetOpeningAngle)
00148             jetOpeningAngle = deltaRToSeed;
00149 
00150          tmpLorentzVect+=(**iCand).p4();
00151       }
00152 
00153       //Setting the myPFTau four momentum as the one made from the signal cone constituents.
00154       double energy = tmpLorentzVect.energy();
00155       double transverseEnergy = tmpLorentzVect.pt();
00156       myPFTau.setP4(tmpLorentzVect);
00157 
00158       // Compute the cone sizes
00159       double myTrackerSignalConeSize = myPFTauElementsOperators.computeConeSize(
00160             myTrackerSignalConeSizeTFormula, TrackerSignalConeSize_min_, TrackerSignalConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00161       double myTrackerIsolConeSize = myPFTauElementsOperators.computeConeSize(
00162             myTrackerIsolConeSizeTFormula, TrackerIsolConeSize_min_, TrackerIsolConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00163       double myECALSignalConeSize = myPFTauElementsOperators.computeConeSize(
00164             myECALSignalConeSizeTFormula, ECALSignalConeSize_min_, ECALSignalConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00165       double myECALIsolConeSize = myPFTauElementsOperators.computeConeSize(
00166             myECALIsolConeSizeTFormula, ECALIsolConeSize_min_, ECALIsolConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00167       double myHCALSignalConeSize = myPFTauElementsOperators.computeConeSize(
00168             myHCALSignalConeSizeTFormula, HCALSignalConeSize_min_, HCALSignalConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00169       double myHCALIsolConeSize = myPFTauElementsOperators.computeConeSize(
00170             myHCALIsolConeSizeTFormula, HCALIsolConeSize_min_, HCALIsolConeSize_max_, transverseEnergy, energy, jetOpeningAngle);
00171 
00172       // Signal cone collections
00173       PFCandidateRefVector mySignalPFChargedHadrCands, mySignalPFNeutrHadrCands, mySignalPFGammaCands, mySignalPFCands;
00174 
00175       if (UseChargedHadrCandLeadChargedHadrCand_tksDZconstraint_ && myleadPFCand_rectkavailable) {
00176          mySignalPFChargedHadrCands=myPFTauElementsOperators.PFChargedHadrCandsInCone(tauAxis,
00177                TrackerSignalConeMetric_, myTrackerSignalConeSize, PFTauAlgo_ChargedHadrCand_minPt_,
00178                ChargedHadrCandLeadChargedHadrCand_tksmaxDZ_, myleadPFCand_rectkDZ, myPV);
00179       }
00180       else {
00181          mySignalPFChargedHadrCands=myPFTauElementsOperators.PFChargedHadrCandsInCone(tauAxis,
00182                TrackerSignalConeMetric_, myTrackerSignalConeSize, PFTauAlgo_ChargedHadrCand_minPt_);
00183       }
00184 
00185       // Set the Charged hadronics that live in the signal cones
00186       myPFTau.setsignalPFChargedHadrCands(mySignalPFChargedHadrCands);
00187 
00188       // Set the neurtral hadrons that live in the signal cone
00189       mySignalPFNeutrHadrCands=myPFTauElementsOperators.PFNeutrHadrCandsInCone(tauAxis,
00190             HCALSignalConeMetric_, myHCALSignalConeSize, PFTauAlgo_NeutrHadrCand_minPt_);
00191 
00192       myPFTau.setsignalPFNeutrHadrCands(mySignalPFNeutrHadrCands);
00193 
00194       // Compute the gammas that live in the signal cone
00195       mySignalPFGammaCands=myPFTauElementsOperators.PFGammaCandsInCone(tauAxis,
00196             ECALSignalConeMetric_,myECALSignalConeSize,PFTauAlgo_GammaCand_minPt_);
00197 
00198       myPFTau.setsignalPFGammaCands(mySignalPFGammaCands);
00199 
00200       // Add charged objects to signal cone, and calculate charge
00201       if(mySignalPFChargedHadrCands.size() != 0) {
00202          int mySignalPFChargedHadrCands_qsum=0;       
00203          for(size_t i = 0; i < mySignalPFChargedHadrCands.size(); i++) {
00204             mySignalPFChargedHadrCands_qsum += mySignalPFChargedHadrCands[i]->charge();
00205             mySignalPFCands.push_back(mySignalPFChargedHadrCands[i]);
00206          }
00207          myPFTau.setCharge(mySignalPFChargedHadrCands_qsum);    
00208       }
00209 
00210       //Add neutral objects to signal cone
00211       for(size_t i = 0; i < mySignalPFNeutrHadrCands.size(); i++) {
00212          mySignalPFCands.push_back(mySignalPFNeutrHadrCands[i]);
00213       }
00214 
00215       // For the signal gammas, keep track of the highest pt object 
00216       double maxSignalGammaPt = 0.;
00217       for(size_t i = 0; i < mySignalPFGammaCands.size(); i++) {
00218          if(mySignalPFGammaCands[i]->pt() > maxSignalGammaPt) {
00219             myleadPFNeutralCand = mySignalPFGammaCands[i];
00220             maxSignalGammaPt = mySignalPFGammaCands[i]->pt();
00221          }
00222          mySignalPFCands.push_back(mySignalPFGammaCands[i]);
00223       }
00224       myPFTau.setsignalPFCands(mySignalPFCands);
00225       // Set leading gamma
00226       myPFTau.setleadPFNeutralCand(myleadPFNeutralCand); 
00227 
00228       // Logic to determine lead PFCand.  If the lead charged object
00229       // is above the threshold, take that.  If the lead charged object is less
00230       // than the threshold (but exists), AND there exists a gamma above the threshold
00231       // take the gamma as the leadPFCand.  Otherwise it is null.
00232 
00233       if(myleadPFChargedCand->pt() > LeadPFCand_minPt_) {
00234          myPFTau.setleadPFCand(myleadPFChargedCand);
00235       } else if (maxSignalGammaPt > LeadPFCand_minPt_) {
00236          myPFTau.setleadPFCand(myleadPFNeutralCand);
00237       }
00238 
00239       // Declare isolation collections
00240       PFCandidateRefVector myUnfilteredIsolPFChargedHadrCands, myIsolPFNeutrHadrCands, myIsolPFGammaCands, myIsolPFCands;
00241 
00242       // Build unfiltered isolation collection
00243       if(UseChargedHadrCandLeadChargedHadrCand_tksDZconstraint_ && myleadPFCand_rectkavailable) {
00244          myUnfilteredIsolPFChargedHadrCands=myPFTauElementsOperators.PFChargedHadrCandsInAnnulus(
00245                tauAxis,TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,
00246                PFTauAlgo_ChargedHadrCand_minPt_,ChargedHadrCandLeadChargedHadrCand_tksmaxDZ_,myleadPFCand_rectkDZ, myPV);
00247       } else {
00248          myUnfilteredIsolPFChargedHadrCands=myPFTauElementsOperators.PFChargedHadrCandsInAnnulus(
00249                tauAxis,TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,
00250                PFTauAlgo_ChargedHadrCand_minPt_);
00251       }
00252 
00253       // Filter isolation annulus charge dhadrons with additional nHits quality cut 
00254       // (note that other cuts [pt, chi2, are already cut on])
00255       PFCandidateRefVector myIsolPFChargedHadrCands;
00256       myIsolPFChargedHadrCands = TauTagTools::filteredPFChargedHadrCandsByNumTrkHits(
00257             myUnfilteredIsolPFChargedHadrCands, ChargedHadrCand_IsolAnnulus_minNhits_); 
00258 
00259       myPFTau.setisolationPFChargedHadrCands(myIsolPFChargedHadrCands);
00260 
00261       // Fill neutral hadrons
00262       myIsolPFNeutrHadrCands = myPFTauElementsOperators.PFNeutrHadrCandsInAnnulus(
00263             tauAxis, HCALSignalConeMetric_, myHCALSignalConeSize, HCALIsolConeMetric_, 
00264             myHCALIsolConeSize, PFTauAlgo_NeutrHadrCand_minPt_);
00265       myPFTau.setisolationPFNeutrHadrCands(myIsolPFNeutrHadrCands);
00266 
00267       // Fill gamma candidates
00268       myIsolPFGammaCands = myPFTauElementsOperators.PFGammaCandsInAnnulus(
00269             tauAxis, ECALSignalConeMetric_, myECALSignalConeSize, ECALIsolConeMetric_, 
00270             myECALIsolConeSize, PFTauAlgo_GammaCand_minPt_);  
00271       myPFTau.setisolationPFGammaCands(myIsolPFGammaCands);
00272 
00273       //Incorporate converted gammas from isolation ellipse into signal  ... ELLL
00274       //Get pair with in/out elements using the isoPFGammaCandidates set by default
00275       if(AddEllipseGammas_) {
00276          double rPhi;
00277          if(Rphi_ >= 1.) 
00278             rPhi = Rphi_*myECALSignalConeSize;
00279          else 
00280             rPhi = Rphi_;
00281 
00282          std::pair<PFCandidateRefVector,PFCandidateRefVector> elementsInOutEllipse = 
00283             myPFTauElementsOperators.PFGammaCandsInOutEllipse(myIsolPFGammaCands, *myleadPFCand, rPhi, myECALSignalConeSize, MaxEtInEllipse_);
00284 
00285          PFCandidateRefVector elementsInEllipse = elementsInOutEllipse.first;
00286          PFCandidateRefVector elementsOutEllipse = elementsInOutEllipse.second;
00287          //add the inside elements to signal PFCandidates and reset signal PFCands
00288          for(PFCandidateRefVector::const_iterator inEllipseIt = elementsInEllipse.begin(); inEllipseIt != elementsInEllipse.end(); inEllipseIt++){
00289             mySignalPFCands.push_back(*inEllipseIt);
00290             mySignalPFGammaCands.push_back(*inEllipseIt);
00291          }
00292          myPFTau.setsignalPFCands(mySignalPFCands);
00293          //redefine isoPFGammaCandidates to be the outside elements
00294          myIsolPFGammaCands=elementsOutEllipse;
00295          myPFTau.setisolationPFGammaCands(myIsolPFGammaCands);
00296       }
00297 
00298 
00299       // Fill isolation collections, and calculate pt sum in isolation cone
00300       float myIsolPFChargedHadrCands_Ptsum = 0.;
00301       float myIsolPFGammaCands_Etsum       = 0.;
00302       for(size_t i = 0; i < myIsolPFChargedHadrCands.size(); i++) {
00303          myIsolPFChargedHadrCands_Ptsum += myIsolPFChargedHadrCands[i]->pt();
00304          myIsolPFCands.push_back(myIsolPFChargedHadrCands[i]);
00305       }
00306       myPFTau.setisolationPFChargedHadrCandsPtSum(myIsolPFChargedHadrCands_Ptsum);
00307 
00308       // Put neutral hadrons into collection
00309       for(size_t i = 0; i < myIsolPFNeutrHadrCands.size(); i++) {
00310          myIsolPFCands.push_back(myIsolPFNeutrHadrCands[i]);
00311       }
00312 
00313       for(size_t i = 0; i < myIsolPFGammaCands.size(); i++) {
00314          myIsolPFGammaCands_Etsum += myIsolPFGammaCands[i]->et();
00315          myIsolPFCands.push_back(myIsolPFGammaCands[i]);
00316       } 
00317       myPFTau.setisolationPFGammaCandsEtSum(myIsolPFGammaCands_Etsum);
00318       myPFTau.setisolationPFCands(myIsolPFCands);
00319 
00320       //Making the alternateLorentzVector, i.e. direction with only signal components
00321       math::XYZTLorentzVector alternatLorentzVect(0.,0.,0.,0.);
00322       for (PFCandidateRefVector::const_iterator iGammaCand = mySignalPFGammaCands.begin(); 
00323             iGammaCand != mySignalPFGammaCands.end(); iGammaCand++) {
00324          alternatLorentzVect+=(**iGammaCand).p4();
00325       }
00326 
00327       for (PFCandidateRefVector::const_iterator iChargedHadrCand = mySignalPFChargedHadrCands.begin();
00328             iChargedHadrCand != mySignalPFChargedHadrCands.end(); iChargedHadrCand++) {
00329          alternatLorentzVect+=(**iChargedHadrCand).p4();  
00330       }
00331       myPFTau.setalternatLorentzVect(alternatLorentzVect);
00332       myPFTau.setP4(alternatLorentzVect);
00333 
00334       // Set tau vertex as PV vertex
00335       myPFTau.setVertex(math::XYZPoint(myPV.x(), myPV.y(), myPV.z()));
00336    }  
00337 
00338    // set the leading, signal cone and isolation annulus Tracks (the initial list of Tracks was catched through a JetTracksAssociation 
00339    // object, not through the charged hadr. PFCandidates inside the PFJet ; 
00340    // the motivation for considering these objects is the need for checking that a selection by the 
00341    // charged hadr. PFCandidates is equivalent to a selection by the rec. Tracks.)
00342    TrackRef myleadTk=myPFTauElementsOperators.leadTk(MatchingConeMetric_,myMatchingConeSize,LeadTrack_minPt_);
00343    myPFTau.setleadTrack(myleadTk);
00344    if(myleadTk.isNonnull()){
00345       double myleadTkDZ = (*myleadTk).dz(myPV.position());
00346       double myTrackerSignalConeSize=myPFTauElementsOperators.computeConeSize(myTrackerSignalConeSizeTFormula,TrackerSignalConeSize_min_,TrackerSignalConeSize_max_);
00347       double myTrackerIsolConeSize=myPFTauElementsOperators.computeConeSize(myTrackerIsolConeSizeTFormula,TrackerIsolConeSize_min_,TrackerIsolConeSize_max_);     
00348       if (UseTrackLeadTrackDZconstraint_){
00349          myPFTau.setsignalTracks(myPFTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,PFTauAlgo_Track_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ, myPV));
00350 
00351          TrackRefVector myUnfilteredTracks = myPFTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,PFTauAlgo_Track_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ, myPV);
00352          TrackRefVector myFilteredTracks = TauTagTools::filteredTracksByNumTrkHits(myUnfilteredTracks, Track_IsolAnnulus_minNhits_);
00353          myPFTau.setisolationTracks(myFilteredTracks);
00354 
00355       }else{
00356          myPFTau.setsignalTracks(myPFTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,PFTauAlgo_Track_minPt_));
00357 
00358          TrackRefVector myUnfilteredTracks = myPFTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,PFTauAlgo_Track_minPt_);
00359          TrackRefVector myFilteredTracks = TauTagTools::filteredTracksByNumTrkHits(myUnfilteredTracks, Track_IsolAnnulus_minNhits_);
00360          myPFTau.setisolationTracks(myFilteredTracks);
00361       }
00362    }
00363 
00364 
00365    /* For elecron rejection */
00366    double myECALenergy             =  0.;
00367    double myHCALenergy             =  0.;
00368    double myHCALenergy3x3          =  0.;
00369    double myMaximumHCALPFClusterE  =  0.;
00370    double myMaximumHCALPFClusterEt =  0.;
00371    double myStripClusterE          =  0.;
00372    double myEmfrac                 = -1.;
00373    double myElectronPreIDOutput    = -1111.;
00374    bool   myElecPreid              =  false;
00375    reco::TrackRef myElecTrk;
00376 
00377    typedef std::pair<reco::PFBlockRef, unsigned> ElementInBlock;
00378    typedef std::vector< ElementInBlock > ElementsInBlocks; 
00379 
00380    //Use the electron rejection only in case there is a charged leading pion
00381    if(myleadPFChargedCand.isNonnull()){
00382       myElectronPreIDOutput = myleadPFChargedCand->mva_e_pi();
00383 
00384       math::XYZPointF myElecTrkEcalPos = myleadPFChargedCand->positionAtECALEntrance();
00385       myElecTrk = myleadPFChargedCand->trackRef();//Electron candidate
00386 
00387       if(myElecTrk.isNonnull()) {
00388          //FROM AOD
00389          if(DataType_ == "AOD"){
00390             // Corrected Cluster energies
00391             for(int i=0;i<(int)myPFCands.size();i++){
00392                myHCALenergy += myPFCands[i]->hcalEnergy();
00393                myECALenergy += myPFCands[i]->ecalEnergy();
00394 
00395                math::XYZPointF candPos;
00396                if (myPFCands[i]->particleId()==1 || myPFCands[i]->particleId()==2)//if charged hadron or electron
00397                   candPos = myPFCands[i]->positionAtECALEntrance();
00398                else
00399                   candPos = math::XYZPointF(myPFCands[i]->px(),myPFCands[i]->py(),myPFCands[i]->pz());
00400 
00401                double deltaR   = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos,candPos);
00402                double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos,candPos);
00403                double deltaEta = std::abs(myElecTrkEcalPos.eta()-candPos.eta());
00404                double deltaPhiOverQ = deltaPhi/(double)myElecTrk->charge();
00405 
00406                if (myPFCands[i]->ecalEnergy() >= EcalStripSumE_minClusEnergy_ && deltaEta < EcalStripSumE_deltaEta_ &&
00407                      deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_  && deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) {
00408                   myStripClusterE += myPFCands[i]->ecalEnergy();
00409                }
00410                if (deltaR<0.184) {
00411                   myHCALenergy3x3 += myPFCands[i]->hcalEnergy();
00412                }
00413                if (myPFCands[i]->hcalEnergy()>myMaximumHCALPFClusterE) {
00414                   myMaximumHCALPFClusterE = myPFCands[i]->hcalEnergy();
00415                }
00416                if ((myPFCands[i]->hcalEnergy()*fabs(sin(candPos.Theta())))>myMaximumHCALPFClusterEt) {
00417                   myMaximumHCALPFClusterEt = (myPFCands[i]->hcalEnergy()*fabs(sin(candPos.Theta())));
00418                }
00419             }
00420 
00421          } else if(DataType_ == "RECO"){ //From RECO
00422             // Against double counting of clusters
00423             std::vector<math::XYZPoint> hcalPosV; hcalPosV.clear();
00424             std::vector<math::XYZPoint> ecalPosV; ecalPosV.clear();
00425             for(int i=0;i<(int)myPFCands.size();i++){
00426                const ElementsInBlocks& elts = myPFCands[i]->elementsInBlocks();
00427                for(ElementsInBlocks::const_iterator it=elts.begin(); it!=elts.end(); ++it) {
00428                   const reco::PFBlock& block = *(it->first);
00429                   unsigned indexOfElementInBlock = it->second;
00430                   const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
00431                   assert(indexOfElementInBlock<elements.size());
00432 
00433                   const reco::PFBlockElement& element = elements[indexOfElementInBlock];
00434 
00435                   if(element.type()==reco::PFBlockElement::HCAL) {
00436                      math::XYZPoint clusPos = element.clusterRef()->position();
00437                      double en = (double)element.clusterRef()->energy();
00438                      double et = (double)element.clusterRef()->energy()*fabs(sin(clusPos.Theta()));
00439                      if (en>myMaximumHCALPFClusterE) {
00440                         myMaximumHCALPFClusterE = en;
00441                      }
00442                      if (et>myMaximumHCALPFClusterEt) {
00443                         myMaximumHCALPFClusterEt = et;
00444                      }
00445                      if (!checkPos(hcalPosV,clusPos)) {
00446                         hcalPosV.push_back(clusPos);
00447                         myHCALenergy += en;
00448                         double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos,clusPos);
00449                         if (deltaR<0.184) {
00450                            myHCALenergy3x3 += en;
00451                         }
00452                      }
00453                   } else if(element.type()==reco::PFBlockElement::ECAL) {
00454                      double en = (double)element.clusterRef()->energy();
00455                      math::XYZPoint clusPos = element.clusterRef()->position();
00456                      if (!checkPos(ecalPosV,clusPos)) {
00457                         ecalPosV.push_back(clusPos);
00458                         myECALenergy += en;
00459                         double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos,clusPos);
00460                         double deltaEta = std::abs(myElecTrkEcalPos.eta()-clusPos.eta());
00461                         double deltaPhiOverQ = deltaPhi/(double)myElecTrk->charge();
00462                         if (en >= EcalStripSumE_minClusEnergy_ && deltaEta<EcalStripSumE_deltaEta_ && deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ && deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) { 
00463                            myStripClusterE += en;
00464                         }
00465                      }    
00466                   }
00467                } //end elements in blocks
00468             } //end loop over PFcands
00469          } //end RECO case       
00470       } // end check for null electrk  
00471    } // end check for null pfChargedHadrCand
00472 
00473    if ((myHCALenergy+myECALenergy)>0.)
00474       myEmfrac = myECALenergy/(myHCALenergy+myECALenergy);
00475    myPFTau.setemFraction((float)myEmfrac);
00476 
00477    // scale the appropriate quantities by the momentum of the electron if it exists
00478    if (myElecTrk.isNonnull())
00479    {
00480       float myElectronMomentum = (float)myElecTrk->p();
00481       if (myElectronMomentum > 0.)
00482       {
00483          myHCALenergy            /= myElectronMomentum;
00484          myMaximumHCALPFClusterE /= myElectronMomentum;
00485          myHCALenergy3x3         /= myElectronMomentum;
00486          myStripClusterE         /= myElectronMomentum;
00487       }
00488    }
00489    myPFTau.sethcalTotOverPLead((float)myHCALenergy);
00490    myPFTau.sethcalMaxOverPLead((float)myMaximumHCALPFClusterE);
00491    myPFTau.sethcal3x3OverPLead((float)myHCALenergy3x3);
00492    myPFTau.setecalStripSumEOverPLead((float)myStripClusterE);
00493    myPFTau.setmaximumHCALPFClusterEt(myMaximumHCALPFClusterEt);
00494    myPFTau.setelectronPreIDOutput(myElectronPreIDOutput);
00495    if (myElecTrk.isNonnull()) 
00496       myPFTau.setelectronPreIDTrack(myElecTrk);
00497    if (myElectronPreIDOutput > maximumForElectrionPreIDOutput_)
00498       myElecPreid = true;
00499    myPFTau.setelectronPreIDDecision(myElecPreid);
00500 
00501    // These need to be filled!
00502    //myPFTau.setbremsRecoveryEOverPLead(my...);
00503 
00504    /* End elecron rejection */
00505 
00506    return myPFTau;  
00507 }
00508 
00509 bool
00510 PFRecoTauAlgorithm::checkPos(std::vector<math::XYZPoint> CalPos,math::XYZPoint CandPos) const{
00511    bool flag = false;
00512    for (unsigned int i=0;i<CalPos.size();i++) {
00513       if (CalPos[i] == CandPos) {
00514          flag = true;
00515          break;
00516       }
00517    }
00518    return flag;
00519    //return false;
00520 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoTauTag/RecoTau/src/PFRecoTauAlgorithm.cc
