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