CMSSW Reference Manual

00001 #include "RecoTauTag/RecoTau/interface/CaloRecoTauAlgorithm.h"
00002 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00003 
00004 CaloRecoTauAlgorithm::CaloRecoTauAlgorithm() : TransientTrackBuilder_(0),MagneticField_(0),chargedpi_mass_(0.13957018){}  
00005 CaloRecoTauAlgorithm::CaloRecoTauAlgorithm(const ParameterSet& iConfig) : TransientTrackBuilder_(0),MagneticField_(0),chargedpi_mass_(0.13957018){
00006   LeadTrack_minPt_                    = iConfig.getParameter<double>("LeadTrack_minPt");
00007   Track_minPt_                        = iConfig.getParameter<double>("Track_minPt");
00008   UseTrackLeadTrackDZconstraint_      = iConfig.getParameter<bool>("UseTrackLeadTrackDZconstraint");
00009   TrackLeadTrack_maxDZ_               = iConfig.getParameter<double>("TrackLeadTrack_maxDZ");
00010   ECALRecHit_minEt_                   = iConfig.getParameter<double>("ECALRecHit_minEt");       
00011   
00012   MatchingConeMetric_                 = iConfig.getParameter<string>("MatchingConeMetric");
00013   MatchingConeSizeFormula_            = iConfig.getParameter<string>("MatchingConeSizeFormula");
00014   MatchingConeSize_min_               = iConfig.getParameter<double>("MatchingConeSize_min");
00015   MatchingConeSize_max_               = iConfig.getParameter<double>("MatchingConeSize_max");
00016   TrackerSignalConeMetric_            = iConfig.getParameter<string>("TrackerSignalConeMetric");
00017   TrackerSignalConeSizeFormula_       = iConfig.getParameter<string>("TrackerSignalConeSizeFormula");
00018   TrackerSignalConeSize_min_          = iConfig.getParameter<double>("TrackerSignalConeSize_min");
00019   TrackerSignalConeSize_max_          = iConfig.getParameter<double>("TrackerSignalConeSize_max");
00020   TrackerIsolConeMetric_              = iConfig.getParameter<string>("TrackerIsolConeMetric"); 
00021   TrackerIsolConeSizeFormula_         = iConfig.getParameter<string>("TrackerIsolConeSizeFormula"); 
00022   TrackerIsolConeSize_min_            = iConfig.getParameter<double>("TrackerIsolConeSize_min");
00023   TrackerIsolConeSize_max_            = iConfig.getParameter<double>("TrackerIsolConeSize_max");
00024   ECALSignalConeMetric_               = iConfig.getParameter<string>("ECALSignalConeMetric");
00025   ECALSignalConeSizeFormula_          = iConfig.getParameter<string>("ECALSignalConeSizeFormula");    
00026   ECALSignalConeSize_min_             = iConfig.getParameter<double>("ECALSignalConeSize_min");
00027   ECALSignalConeSize_max_             = iConfig.getParameter<double>("ECALSignalConeSize_max");
00028   ECALIsolConeMetric_                 = iConfig.getParameter<string>("ECALIsolConeMetric");
00029   ECALIsolConeSizeFormula_            = iConfig.getParameter<string>("ECALIsolConeSizeFormula");      
00030   ECALIsolConeSize_min_               = iConfig.getParameter<double>("ECALIsolConeSize_min");
00031   ECALIsolConeSize_max_               = iConfig.getParameter<double>("ECALIsolConeSize_max");
00032   
00033   AreaMetric_recoElements_maxabsEta_  = iConfig.getParameter<double>("AreaMetric_recoElements_maxabsEta");
00034 
00035   //Computing the TFormula
00036   myTrackerSignalConeSizeTFormula=TauTagTools::computeConeSizeTFormula(TrackerSignalConeSizeFormula_,"Tracker signal cone size");
00037   myTrackerIsolConeSizeTFormula=TauTagTools::computeConeSizeTFormula(TrackerIsolConeSizeFormula_,"Tracker isolation cone size");
00038   myECALSignalConeSizeTFormula=TauTagTools::computeConeSizeTFormula(ECALSignalConeSizeFormula_,"ECAL signal cone size");
00039   myECALIsolConeSizeTFormula=TauTagTools::computeConeSizeTFormula(ECALIsolConeSizeFormula_,"ECAL isolation cone size");
00040   myMatchingConeSizeTFormula=TauTagTools::computeConeSizeTFormula(MatchingConeSizeFormula_,"Matching cone size");
00041 
00042  
00043 }
00044 void CaloRecoTauAlgorithm::setTransientTrackBuilder(const TransientTrackBuilder* x){TransientTrackBuilder_=x;}
00045 void CaloRecoTauAlgorithm::setMagneticField(const MagneticField* x){MagneticField_=x;} 
00046 
00047 CaloTau CaloRecoTauAlgorithm::buildCaloTau(Event& iEvent,const EventSetup& iSetup,const CaloTauTagInfoRef& myCaloTauTagInfoRef,const Vertex& myPV){
00048   CaloJetRef myCaloJet=(*myCaloTauTagInfoRef).calojetRef(); // catch a ref to the initial CaloJet  
00049   const vector<CaloTowerPtr> myCaloTowers=(*myCaloJet).getCaloConstituents();
00050   CaloTau myCaloTau(numeric_limits<int>::quiet_NaN(),myCaloJet->p4()); // create the CaloTau with the initial CaloJet Lorentz-vector
00051   
00052   myCaloTau.setcaloTauTagInfoRef(myCaloTauTagInfoRef);
00053   
00054   TrackRefVector myTks=(*myCaloTauTagInfoRef).Tracks();
00055   
00056   CaloTauElementsOperators myCaloTauElementsOperators(myCaloTau);
00057   double myMatchingConeSize=myCaloTauElementsOperators.computeConeSize(myMatchingConeSizeTFormula,MatchingConeSize_min_,MatchingConeSize_max_);
00058   TrackRef myleadTk=myCaloTauElementsOperators.leadTk(MatchingConeMetric_,myMatchingConeSize,LeadTrack_minPt_);
00059   double myCaloTau_refInnerPosition_x=0.;
00060   double myCaloTau_refInnerPosition_y=0.;
00061   double myCaloTau_refInnerPosition_z=0.;
00062   if(myleadTk.isNonnull()){
00063     myCaloTau.setleadTrack(myleadTk);
00064     double myleadTkDZ=(*myleadTk).dz();
00065     if(TransientTrackBuilder_!=0){ 
00066       const TransientTrack myleadTransientTk=TransientTrackBuilder_->build(&(*myleadTk));
00067       GlobalVector myCaloJetdir((*myCaloJet).px(),(*myCaloJet).py(),(*myCaloJet).pz());
00068       if(IPTools::signedTransverseImpactParameter(myleadTransientTk,myCaloJetdir,myPV).first)
00069         myCaloTau.setleadTracksignedSipt(IPTools::signedTransverseImpactParameter(myleadTransientTk,myCaloJetdir,myPV).second.significance());
00070     }
00071     if((*myleadTk).innerOk()){
00072       myCaloTau_refInnerPosition_x=(*myleadTk).innerPosition().x(); 
00073       myCaloTau_refInnerPosition_y=(*myleadTk).innerPosition().y(); 
00074       myCaloTau_refInnerPosition_z=(*myleadTk).innerPosition().z(); 
00075     }
00076     
00077     if(MagneticField_!=0){ 
00078       math::XYZPoint mypropagleadTrackECALSurfContactPoint=TauTagTools::propagTrackECALSurfContactPoint(MagneticField_,myleadTk);
00079       if(mypropagleadTrackECALSurfContactPoint.R()!=0.){
00080         double myleadTrackHCAL3x3hottesthitDEta=0.;
00081         double myleadTrackHCAL3x3hottesthitEt=0.;
00082         double myleadTrackHCAL3x3hitsEtSum=0.;
00083         ESHandle<CaloGeometry> myCaloGeometry;
00084         iSetup.get<CaloGeometryRecord>().get(myCaloGeometry);
00085         const CaloSubdetectorGeometry* myCaloSubdetectorGeometry=(*myCaloGeometry).getSubdetectorGeometry(DetId::Calo,CaloTowerDetId::SubdetId);
00086         CaloTowerDetId mypropagleadTrack_closestCaloTowerId((*myCaloSubdetectorGeometry).getClosestCell(GlobalPoint(mypropagleadTrackECALSurfContactPoint.x(),
00087                                                                                                                     mypropagleadTrackECALSurfContactPoint.y(),
00088                                                                                                                     mypropagleadTrackECALSurfContactPoint.z())));
00089         vector<CaloTowerDetId> mypropagleadTrack_closestCaloTowerNeighbourIds=getCaloTowerneighbourDetIds(myCaloSubdetectorGeometry,mypropagleadTrack_closestCaloTowerId);
00090         for(vector<CaloTowerPtr>::const_iterator iCaloTower=myCaloTowers.begin();iCaloTower!=myCaloTowers.end();iCaloTower++){
00091           CaloTowerDetId iCaloTowerId((**iCaloTower).id());
00092           bool CaloTower_inside3x3matrix=false;
00093           if (iCaloTowerId==mypropagleadTrack_closestCaloTowerId) CaloTower_inside3x3matrix=true;
00094           if (!CaloTower_inside3x3matrix){
00095             for(vector<CaloTowerDetId>::const_iterator iCaloTowerDetId=mypropagleadTrack_closestCaloTowerNeighbourIds.begin();iCaloTowerDetId!=mypropagleadTrack_closestCaloTowerNeighbourIds.end();iCaloTowerDetId++){
00096               if (iCaloTowerId==(*iCaloTowerDetId)){ 
00097                 CaloTower_inside3x3matrix=true;
00098                 break;
00099               }
00100             }
00101           }
00102           if (!CaloTower_inside3x3matrix) continue;       
00103           myleadTrackHCAL3x3hitsEtSum+=(**iCaloTower).hadEt();
00104           if((**iCaloTower).hadEt()>=myleadTrackHCAL3x3hottesthitEt ){
00105             if ((**iCaloTower).hadEt()!=myleadTrackHCAL3x3hottesthitEt || 
00106                 ((**iCaloTower).hadEt()==myleadTrackHCAL3x3hottesthitEt && fabs((**iCaloTower).eta()-mypropagleadTrackECALSurfContactPoint.Eta())<myleadTrackHCAL3x3hottesthitDEta)) myleadTrackHCAL3x3hottesthitDEta = fabs((**iCaloTower).eta()-mypropagleadTrackECALSurfContactPoint.Eta());
00107             myleadTrackHCAL3x3hottesthitEt=(**iCaloTower).hadEt();
00108           }     
00109         }       
00110         myCaloTau.setleadTrackHCAL3x3hitsEtSum(myleadTrackHCAL3x3hitsEtSum);
00111         if (myleadTrackHCAL3x3hottesthitEt!=0.) myCaloTau.setleadTrackHCAL3x3hottesthitDEta(myleadTrackHCAL3x3hottesthitDEta);
00112       }
00113     }
00114     
00115     if (UseTrackLeadTrackDZconstraint_){
00116       TrackRefVector myTksbis;
00117       for (TrackRefVector::const_iterator iTrack=myTks.begin();iTrack!=myTks.end();++iTrack) {
00118         if (fabs((**iTrack).dz()-myleadTkDZ)<=TrackLeadTrack_maxDZ_) myTksbis.push_back(*iTrack);
00119       }
00120       myTks=myTksbis;
00121     }
00122 
00123 
00124     double myTrackerSignalConeSize=myCaloTauElementsOperators.computeConeSize(myTrackerSignalConeSizeTFormula,TrackerSignalConeSize_min_,TrackerSignalConeSize_max_);
00125         double myTrackerIsolConeSize=myCaloTauElementsOperators.computeConeSize(myTrackerIsolConeSizeTFormula,TrackerIsolConeSize_min_,TrackerIsolConeSize_max_);       
00126     double myECALSignalConeSize=myCaloTauElementsOperators.computeConeSize(myECALSignalConeSizeTFormula,ECALSignalConeSize_min_,ECALSignalConeSize_max_);
00127     double myECALIsolConeSize=myCaloTauElementsOperators.computeConeSize(myECALIsolConeSizeTFormula,ECALIsolConeSize_min_,ECALIsolConeSize_max_);       
00128 
00129     TrackRefVector mySignalTks;
00130     if (UseTrackLeadTrackDZconstraint_) mySignalTks=myCaloTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,Track_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ);
00131     else mySignalTks=myCaloTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,Track_minPt_);
00132     myCaloTau.setsignalTracks(mySignalTks);
00133     
00134     // setting invariant mass of the signal Tracks system
00135     math::XYZTLorentzVector mySignalTksInvariantMass(0.,0.,0.,0.);
00136     if((int)(mySignalTks.size())!=0){
00137       int mySignalTks_qsum=0;       
00138       for(int i=0;i<(int)mySignalTks.size();i++){
00139         mySignalTks_qsum+=mySignalTks[i]->charge();
00140         math::XYZTLorentzVector mychargedpicand_fromTk_LorentzVect(mySignalTks[i]->momentum().x(),mySignalTks[i]->momentum().y(),mySignalTks[i]->momentum().z(),sqrt(pow((double)mySignalTks[i]->momentum().r(),2)+pow(chargedpi_mass_,2)));
00141         mySignalTksInvariantMass+=mychargedpicand_fromTk_LorentzVect;
00142       }
00143       myCaloTau.setCharge(mySignalTks_qsum);    
00144     }
00145     myCaloTau.setsignalTracksInvariantMass(mySignalTksInvariantMass.mass());
00146     
00147     TrackRefVector myIsolTks;
00148     if (UseTrackLeadTrackDZconstraint_) myIsolTks=myCaloTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,Track_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ);
00149     else myIsolTks=myCaloTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,Track_minPt_);
00150     myCaloTau.setisolationTracks(myIsolTks);
00151     
00152     // setting sum of Pt of the isolation annulus Tracks
00153     float myIsolTks_Ptsum=0.;
00154     for(int i=0;i<(int)myIsolTks.size();i++) myIsolTks_Ptsum+=myIsolTks[i]->pt();
00155     myCaloTau.setisolationTracksPtSum(myIsolTks_Ptsum);
00156     
00157     // setting sum of Et of the isolation annulus ECAL RecHits
00158     float myIsolEcalRecHits_EtSum=0.;
00159     vector<pair<math::XYZPoint,float> > myIsolPositionAndEnergyEcalRecHits=myCaloTauElementsOperators.EcalRecHitsInAnnulus((*myleadTk).momentum(),ECALSignalConeMetric_,myECALSignalConeSize,ECALIsolConeMetric_,myECALIsolConeSize,ECALRecHit_minEt_);
00160     for(vector<pair<math::XYZPoint,float> >::const_iterator iEcalRecHit=myIsolPositionAndEnergyEcalRecHits.begin();iEcalRecHit!=myIsolPositionAndEnergyEcalRecHits.end();iEcalRecHit++){
00161       myIsolEcalRecHits_EtSum+=(*iEcalRecHit).second*fabs(sin((*iEcalRecHit).first.theta()));
00162     }
00163     myCaloTau.setisolationECALhitsEtSum(myIsolEcalRecHits_EtSum);    
00164   }
00165 
00166   math::XYZTLorentzVector myTks_XYZTLorentzVect(0.,0.,0.,0.);
00167   math::XYZTLorentzVector alternatLorentzVect(0.,0.,0.,0.);
00168   for(TrackRefVector::iterator iTrack=myTks.begin();iTrack!=myTks.end();iTrack++) {
00169     // build a charged pion candidate Lorentz vector from a Track
00170     math::XYZTLorentzVector iChargedPionCand_XYZTLorentzVect((**iTrack).momentum().x(),(**iTrack).momentum().y(),(**iTrack).momentum().z(),sqrt(pow((double)(**iTrack).momentum().r(),2)+pow(chargedpi_mass_,2)));
00171     myTks_XYZTLorentzVect+=iChargedPionCand_XYZTLorentzVect;
00172     alternatLorentzVect+=iChargedPionCand_XYZTLorentzVect;
00173   }
00174   myCaloTau.setTracksInvariantMass(myTks_XYZTLorentzVect.mass());
00175   vector<BasicClusterRef> myneutralECALBasicClusters=(*myCaloTauTagInfoRef).neutralECALBasicClusters();
00176   for(vector<BasicClusterRef>::const_iterator iBasicCluster=myneutralECALBasicClusters.begin();iBasicCluster!=myneutralECALBasicClusters.end();iBasicCluster++) {
00177     // build a gamma candidate Lorentz vector from a neutral ECAL BasicCluster
00178     double iGammaCand_px=(**iBasicCluster).energy()*sin((**iBasicCluster).position().theta())*cos((**iBasicCluster).position().phi());
00179     double iGammaCand_py=(**iBasicCluster).energy()*sin((**iBasicCluster).position().theta())*sin((**iBasicCluster).position().phi());
00180     double iGammaCand_pz=(**iBasicCluster).energy()*cos((**iBasicCluster).position().theta());
00181     math::XYZTLorentzVector iGammaCand_XYZTLorentzVect(iGammaCand_px,iGammaCand_py,iGammaCand_pz,(**iBasicCluster).energy());
00182     alternatLorentzVect+=iGammaCand_XYZTLorentzVect;
00183   }
00184   myCaloTau.setalternatLorentzVect(alternatLorentzVect);
00185   
00186   
00187   myCaloTau.setVertex(math::XYZPoint(myCaloTau_refInnerPosition_x,myCaloTau_refInnerPosition_y,myCaloTau_refInnerPosition_z));
00188     
00189   // setting Et of the highest Et HCAL CaloTower
00190   double mymaxEtHCALtower_Et=0.; 
00191   for(unsigned int iTower=0;iTower<myCaloTowers.size();iTower++){
00192     if((*myCaloTowers[iTower]).hadEt()>=mymaxEtHCALtower_Et) mymaxEtHCALtower_Et=(*myCaloTowers[iTower]).hadEt();
00193   }
00194   myCaloTau.setmaximumHCALhitEt(mymaxEtHCALtower_Et);
00195 
00196   return myCaloTau;  
00197 }
00198 
00199 vector<CaloTowerDetId> CaloRecoTauAlgorithm::getCaloTowerneighbourDetIds(const CaloSubdetectorGeometry* myCaloSubdetectorGeometry,CaloTowerDetId myCaloTowerDetId){
00200   CaloTowerTopology myCaloTowerTopology;
00201   vector<CaloTowerDetId> myCaloTowerneighbourDetIds;
00202   vector<DetId> northDetIds=myCaloTowerTopology.north(myCaloTowerDetId);
00203   vector<DetId> westDetIds=myCaloTowerTopology.west(myCaloTowerDetId);
00204   vector<DetId> northwestDetIds,southwestDetIds;
00205   if (westDetIds.size()>0){
00206     northwestDetIds=myCaloTowerTopology.north(westDetIds[0]);
00207     southwestDetIds=myCaloTowerTopology.south(westDetIds[(int)westDetIds.size()-1]);
00208   }
00209   vector<DetId> southDetIds=myCaloTowerTopology.south(myCaloTowerDetId);
00210   vector<DetId> eastDetIds=myCaloTowerTopology.east(myCaloTowerDetId);
00211   vector<DetId> northeastDetIds,southeastDetIds;
00212   if (eastDetIds.size()>0){
00213     northeastDetIds=myCaloTowerTopology.north(eastDetIds[0]);
00214     southeastDetIds=myCaloTowerTopology.south(eastDetIds[(int)eastDetIds.size()-1]);
00215   }
00216   vector<DetId> myneighbourDetIds=northDetIds;
00217   myneighbourDetIds.insert(myneighbourDetIds.end(),westDetIds.begin(),westDetIds.end());
00218   myneighbourDetIds.insert(myneighbourDetIds.end(),northwestDetIds.begin(),northwestDetIds.end());
00219   myneighbourDetIds.insert(myneighbourDetIds.end(),southwestDetIds.begin(),southwestDetIds.end());
00220   myneighbourDetIds.insert(myneighbourDetIds.end(),southDetIds.begin(),southDetIds.end());
00221   myneighbourDetIds.insert(myneighbourDetIds.end(),eastDetIds.begin(),eastDetIds.end());
00222   myneighbourDetIds.insert(myneighbourDetIds.end(),northeastDetIds.begin(),northeastDetIds.end());
00223   myneighbourDetIds.insert(myneighbourDetIds.end(),southeastDetIds.begin(),southeastDetIds.end());
00224   for(vector<DetId>::const_iterator iDetId=myneighbourDetIds.begin();iDetId!=myneighbourDetIds.end();iDetId++){
00225     CaloTowerDetId iCaloTowerId(*iDetId);
00226     myCaloTowerneighbourDetIds.push_back(iCaloTowerId);
00227   }
00228   return myCaloTowerneighbourDetIds;
00229 }
00230
CaloRecoTauAlgorithm.cc
