Reference Manual

Go to the documentation of this file.
00001 #include "RecoTauTag/RecoTau/interface/CaloRecoTauAlgorithm.h"
00002 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00003 
00004 #include "DataFormats/JetReco/interface/JetCollection.h"
00005 
00006 using namespace reco;
00007 
00008 CaloRecoTauAlgorithm::CaloRecoTauAlgorithm() : TransientTrackBuilder_(0),MagneticField_(0),chargedpi_mass_(0.13957018){}  
00009 CaloRecoTauAlgorithm::CaloRecoTauAlgorithm(const edm::ParameterSet& iConfig) : TransientTrackBuilder_(0),MagneticField_(0),chargedpi_mass_(0.13957018){
00010   LeadTrack_minPt_                    = iConfig.getParameter<double>("LeadTrack_minPt");
00011   Track_minPt_                        = iConfig.getParameter<double>("Track_minPt");
00012   IsolationTrack_minPt_               = iConfig.getParameter<double>("IsolationTrack_minPt");
00013   IsolationTrack_minHits_             = iConfig.getParameter<unsigned int>("IsolationTrack_minHits");
00014   UseTrackLeadTrackDZconstraint_      = iConfig.getParameter<bool>("UseTrackLeadTrackDZconstraint");
00015   TrackLeadTrack_maxDZ_               = iConfig.getParameter<double>("TrackLeadTrack_maxDZ");
00016   ECALRecHit_minEt_                   = iConfig.getParameter<double>("ECALRecHit_minEt");       
00017   
00018   MatchingConeMetric_                 = iConfig.getParameter<std::string>("MatchingConeMetric");
00019   MatchingConeSizeFormula_            = iConfig.getParameter<std::string>("MatchingConeSizeFormula");
00020   MatchingConeSize_min_               = iConfig.getParameter<double>("MatchingConeSize_min");
00021   MatchingConeSize_max_               = iConfig.getParameter<double>("MatchingConeSize_max");
00022   TrackerSignalConeMetric_            = iConfig.getParameter<std::string>("TrackerSignalConeMetric");
00023   TrackerSignalConeSizeFormula_       = iConfig.getParameter<std::string>("TrackerSignalConeSizeFormula");
00024   TrackerSignalConeSize_min_          = iConfig.getParameter<double>("TrackerSignalConeSize_min");
00025   TrackerSignalConeSize_max_          = iConfig.getParameter<double>("TrackerSignalConeSize_max");
00026   TrackerIsolConeMetric_              = iConfig.getParameter<std::string>("TrackerIsolConeMetric"); 
00027   TrackerIsolConeSizeFormula_         = iConfig.getParameter<std::string>("TrackerIsolConeSizeFormula"); 
00028   TrackerIsolConeSize_min_            = iConfig.getParameter<double>("TrackerIsolConeSize_min");
00029   TrackerIsolConeSize_max_            = iConfig.getParameter<double>("TrackerIsolConeSize_max");
00030   ECALSignalConeMetric_               = iConfig.getParameter<std::string>("ECALSignalConeMetric");
00031   ECALSignalConeSizeFormula_          = iConfig.getParameter<std::string>("ECALSignalConeSizeFormula");    
00032   ECALSignalConeSize_min_             = iConfig.getParameter<double>("ECALSignalConeSize_min");
00033   ECALSignalConeSize_max_             = iConfig.getParameter<double>("ECALSignalConeSize_max");
00034   ECALIsolConeMetric_                 = iConfig.getParameter<std::string>("ECALIsolConeMetric");
00035   ECALIsolConeSizeFormula_            = iConfig.getParameter<std::string>("ECALIsolConeSizeFormula");      
00036   ECALIsolConeSize_min_               = iConfig.getParameter<double>("ECALIsolConeSize_min");
00037   ECALIsolConeSize_max_               = iConfig.getParameter<double>("ECALIsolConeSize_max");
00038   
00039   EBRecHitsLabel_                     = iConfig.getParameter<edm::InputTag>("EBRecHitsSource"); 
00040   EERecHitsLabel_                     = iConfig.getParameter<edm::InputTag>("EERecHitsSource"); 
00041   ESRecHitsLabel_                     = iConfig.getParameter<edm::InputTag>("ESRecHitsSource"); 
00042 
00043 
00044 
00045   AreaMetric_recoElements_maxabsEta_  = iConfig.getParameter<double>("AreaMetric_recoElements_maxabsEta");
00046 
00047   //Computing the TFormula
00048   myTrackerSignalConeSizeTFormula=TauTagTools::computeConeSizeTFormula(TrackerSignalConeSizeFormula_,"Tracker signal cone size");
00049   myTrackerIsolConeSizeTFormula=TauTagTools::computeConeSizeTFormula(TrackerIsolConeSizeFormula_,"Tracker isolation cone size");
00050   myECALSignalConeSizeTFormula=TauTagTools::computeConeSizeTFormula(ECALSignalConeSizeFormula_,"ECAL signal cone size");
00051   myECALIsolConeSizeTFormula=TauTagTools::computeConeSizeTFormula(ECALIsolConeSizeFormula_,"ECAL isolation cone size");
00052   myMatchingConeSizeTFormula=TauTagTools::computeConeSizeTFormula(MatchingConeSizeFormula_,"Matching cone size");
00053 
00054   mySelectedDetId_.clear();
00055 }
00056 void CaloRecoTauAlgorithm::setTransientTrackBuilder(const TransientTrackBuilder* x){TransientTrackBuilder_=x;}
00057 void CaloRecoTauAlgorithm::setMagneticField(const MagneticField* x){MagneticField_=x;} 
00058 
00059 CaloTau CaloRecoTauAlgorithm::buildCaloTau(edm::Event& iEvent,const edm::EventSetup& iSetup,const CaloTauTagInfoRef& myCaloTauTagInfoRef,const Vertex& myPV){
00060   CaloJetRef myCaloJet=(*myCaloTauTagInfoRef).calojetRef(); // catch a ref to the initial CaloJet  
00061   const std::vector<CaloTowerPtr> myCaloTowers=(*myCaloJet).getCaloConstituents();
00062   JetBaseRef jetRef = myCaloTauTagInfoRef->jetRef();
00063   CaloTau myCaloTau(std::numeric_limits<int>::quiet_NaN(),jetRef->p4()); // create the CaloTau with the corrected Lorentz-vector
00064   
00065   myCaloTau.setcaloTauTagInfoRef(myCaloTauTagInfoRef);
00066   
00067   TrackRefVector myTks=(*myCaloTauTagInfoRef).Tracks();
00068   
00069   CaloTauElementsOperators myCaloTauElementsOperators(myCaloTau);
00070   double myMatchingConeSize=myCaloTauElementsOperators.computeConeSize(myMatchingConeSizeTFormula,MatchingConeSize_min_,MatchingConeSize_max_);
00071   TrackRef myleadTk=myCaloTauElementsOperators.leadTk(MatchingConeMetric_,myMatchingConeSize,LeadTrack_minPt_);
00072   double myCaloTau_refInnerPosition_x=0.;
00073   double myCaloTau_refInnerPosition_y=0.;
00074   double myCaloTau_refInnerPosition_z=0.;
00075   if(myleadTk.isNonnull()){
00076     myCaloTau.setleadTrack(myleadTk);
00077     double myleadTkDZ=(*myleadTk).dz(myPV.position());
00078     if(TransientTrackBuilder_!=0)
00079     { 
00080       const TransientTrack myleadTransientTk=TransientTrackBuilder_->build(&(*myleadTk));
00081       GlobalVector myCaloJetdir((*myCaloJet).px(),(*myCaloJet).py(),(*myCaloJet).pz());
00082       if(IPTools::signedTransverseImpactParameter(myleadTransientTk,myCaloJetdir,myPV).first)
00083         myCaloTau.setleadTracksignedSipt(IPTools::signedTransverseImpactParameter(myleadTransientTk,myCaloJetdir,myPV).second.significance());
00084     }
00085     if((*myleadTk).innerOk()){
00086       myCaloTau_refInnerPosition_x=(*myleadTk).innerPosition().x(); 
00087       myCaloTau_refInnerPosition_y=(*myleadTk).innerPosition().y(); 
00088       myCaloTau_refInnerPosition_z=(*myleadTk).innerPosition().z(); 
00089     }
00090     
00091     if(MagneticField_!=0){ 
00092       math::XYZPoint mypropagleadTrackECALSurfContactPoint=TauTagTools::propagTrackECALSurfContactPoint(MagneticField_,myleadTk);
00093       if(mypropagleadTrackECALSurfContactPoint.R()!=0.){
00094         double myleadTrackHCAL3x3hottesthitDEta=0.;
00095         double myleadTrackHCAL3x3hottesthitEt=0.;
00096         double myleadTrackHCAL3x3hitsEtSum=0.;
00097         edm::ESHandle<CaloGeometry> myCaloGeometry;
00098         iSetup.get<CaloGeometryRecord>().get(myCaloGeometry);
00099         const CaloSubdetectorGeometry* myCaloSubdetectorGeometry=(*myCaloGeometry).getSubdetectorGeometry(DetId::Calo,CaloTowerDetId::SubdetId);
00100         CaloTowerDetId mypropagleadTrack_closestCaloTowerId((*myCaloSubdetectorGeometry).getClosestCell(GlobalPoint(mypropagleadTrackECALSurfContactPoint.x(),
00101                                                                                                                     mypropagleadTrackECALSurfContactPoint.y(),
00102                                                                                                                     mypropagleadTrackECALSurfContactPoint.z())));
00103         std::vector<CaloTowerDetId> mypropagleadTrack_closestCaloTowerNeighbourIds=getCaloTowerneighbourDetIds(myCaloSubdetectorGeometry,mypropagleadTrack_closestCaloTowerId);
00104         for(std::vector<CaloTowerPtr>::const_iterator iCaloTower=myCaloTowers.begin();iCaloTower!=myCaloTowers.end();iCaloTower++){
00105           CaloTowerDetId iCaloTowerId((**iCaloTower).id());
00106           bool CaloTower_inside3x3matrix=false;
00107           if (iCaloTowerId==mypropagleadTrack_closestCaloTowerId) CaloTower_inside3x3matrix=true;
00108           if (!CaloTower_inside3x3matrix){
00109             for(std::vector<CaloTowerDetId>::const_iterator iCaloTowerDetId=mypropagleadTrack_closestCaloTowerNeighbourIds.begin();iCaloTowerDetId!=mypropagleadTrack_closestCaloTowerNeighbourIds.end();iCaloTowerDetId++){
00110               if (iCaloTowerId==(*iCaloTowerDetId)){ 
00111                 CaloTower_inside3x3matrix=true;
00112                 break;
00113               }
00114             }
00115           }
00116           if (!CaloTower_inside3x3matrix) continue;       
00117           myleadTrackHCAL3x3hitsEtSum+=(**iCaloTower).hadEt();
00118           if((**iCaloTower).hadEt()>=myleadTrackHCAL3x3hottesthitEt ){
00119             if ((**iCaloTower).hadEt()!=myleadTrackHCAL3x3hottesthitEt || 
00120                 ((**iCaloTower).hadEt()==myleadTrackHCAL3x3hottesthitEt && fabs((**iCaloTower).eta()-mypropagleadTrackECALSurfContactPoint.Eta())<myleadTrackHCAL3x3hottesthitDEta)) myleadTrackHCAL3x3hottesthitDEta = fabs((**iCaloTower).eta()-mypropagleadTrackECALSurfContactPoint.Eta());
00121             myleadTrackHCAL3x3hottesthitEt=(**iCaloTower).hadEt();
00122           }     
00123         }       
00124         myCaloTau.setleadTrackHCAL3x3hitsEtSum(myleadTrackHCAL3x3hitsEtSum);
00125         if (myleadTrackHCAL3x3hottesthitEt!=0.) myCaloTau.setleadTrackHCAL3x3hottesthitDEta(myleadTrackHCAL3x3hottesthitDEta);
00126       }
00127     }
00128     
00129     if (UseTrackLeadTrackDZconstraint_){
00130       TrackRefVector myTksbis;
00131       for (TrackRefVector::const_iterator iTrack=myTks.begin();iTrack!=myTks.end();++iTrack) {
00132         if (fabs((**iTrack).dz(myPV.position())-myleadTkDZ)<=TrackLeadTrack_maxDZ_) myTksbis.push_back(*iTrack);
00133       }
00134       myTks=myTksbis;
00135     }
00136 
00137 
00138     double myTrackerSignalConeSize=myCaloTauElementsOperators.computeConeSize(myTrackerSignalConeSizeTFormula,TrackerSignalConeSize_min_,TrackerSignalConeSize_max_);
00139     double myTrackerIsolConeSize=myCaloTauElementsOperators.computeConeSize(myTrackerIsolConeSizeTFormula,TrackerIsolConeSize_min_,TrackerIsolConeSize_max_);           
00140     double myECALSignalConeSize=myCaloTauElementsOperators.computeConeSize(myECALSignalConeSizeTFormula,ECALSignalConeSize_min_,ECALSignalConeSize_max_);
00141     double myECALIsolConeSize=myCaloTauElementsOperators.computeConeSize(myECALIsolConeSizeTFormula,ECALIsolConeSize_min_,ECALIsolConeSize_max_);       
00142 
00143     TrackRefVector mySignalTks;
00144     if (UseTrackLeadTrackDZconstraint_) mySignalTks=myCaloTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,Track_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ, myPV);
00145     else mySignalTks=myCaloTauElementsOperators.tracksInCone((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,Track_minPt_);
00146     myCaloTau.setsignalTracks(mySignalTks);
00147     
00148     // setting invariant mass of the signal Tracks system
00149     math::XYZTLorentzVector mySignalTksInvariantMass(0.,0.,0.,0.);
00150     if((int)(mySignalTks.size())!=0){
00151       int mySignalTks_qsum=0;       
00152       for(int i=0;i<(int)mySignalTks.size();i++){
00153         mySignalTks_qsum+=mySignalTks[i]->charge();
00154         math::XYZTLorentzVector mychargedpicand_fromTk_LorentzVect(mySignalTks[i]->momentum().x(),mySignalTks[i]->momentum().y(),mySignalTks[i]->momentum().z(),sqrt(std::pow((double)mySignalTks[i]->momentum().r(),2)+std::pow(chargedpi_mass_,2)));
00155         mySignalTksInvariantMass+=mychargedpicand_fromTk_LorentzVect;
00156       }
00157       myCaloTau.setCharge(mySignalTks_qsum);    
00158     }
00159     myCaloTau.setsignalTracksInvariantMass(mySignalTksInvariantMass.mass());
00160     
00161     TrackRefVector myIsolTks;
00162     if (UseTrackLeadTrackDZconstraint_) myIsolTks=myCaloTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,IsolationTrack_minPt_,TrackLeadTrack_maxDZ_,myleadTkDZ, myPV);
00163     else myIsolTks=myCaloTauElementsOperators.tracksInAnnulus((*myleadTk).momentum(),TrackerSignalConeMetric_,myTrackerSignalConeSize,TrackerIsolConeMetric_,myTrackerIsolConeSize,IsolationTrack_minPt_);
00164     myIsolTks = TauTagTools::filteredTracksByNumTrkHits(myIsolTks,IsolationTrack_minHits_);
00165     myCaloTau.setisolationTracks(myIsolTks);
00166     
00167     // setting sum of Pt of the isolation annulus Tracks
00168     float myIsolTks_Ptsum=0.;
00169     for(int i=0;i<(int)myIsolTks.size();i++) myIsolTks_Ptsum+=myIsolTks[i]->pt();
00170     myCaloTau.setisolationTracksPtSum(myIsolTks_Ptsum);
00171 
00172 
00173     //getting the EcalRecHits. Just take them all
00174   std::vector<std::pair<math::XYZPoint,float> > thePositionAndEnergyEcalRecHits;
00175   mySelectedDetId_.clear();
00176   //  std::vector<CaloTowerPtr> theCaloTowers=myCaloJet->getCaloConstituents();
00177   edm::ESHandle<CaloGeometry> theCaloGeometry;
00178   iSetup.get<CaloGeometryRecord>().get(theCaloGeometry);
00179   const CaloSubdetectorGeometry* theCaloSubdetectorGeometry;  
00180   edm::Handle<EBRecHitCollection> EBRecHits;
00181   edm::Handle<EERecHitCollection> EERecHits; 
00182   edm::Handle<ESRecHitCollection> ESRecHits; 
00183   iEvent.getByLabel(EBRecHitsLabel_,EBRecHits);
00184   iEvent.getByLabel(EERecHitsLabel_,EERecHits);
00185   iEvent.getByLabel(ESRecHitsLabel_,ESRecHits);
00186   double maxDeltaR = 0.8;
00187     math::XYZPoint myCaloJetdir((*myCaloJet).px(),(*myCaloJet).py(),(*myCaloJet).pz());
00188     
00189   for(EBRecHitCollection::const_iterator theRecHit = EBRecHits->begin();theRecHit != EBRecHits->end(); theRecHit++){
00190     theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00191     const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(theRecHit->id());  
00192     math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00193     if(ROOT::Math::VectorUtil::DeltaR(myCaloJetdir,theRecHitCell_XYZPoint) < maxDeltaR){
00194       std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00195       thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00196       mySelectedDetId_.push_back(theRecHit->id());
00197     }
00198   }
00199 
00200 for(EERecHitCollection::const_iterator theRecHit = EERecHits->begin();theRecHit != EERecHits->end(); theRecHit++){
00201     theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00202     const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(theRecHit->id());  
00203     math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00204     if(ROOT::Math::VectorUtil::DeltaR(myCaloJetdir,theRecHitCell_XYZPoint) < maxDeltaR){
00205       std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00206       thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00207       mySelectedDetId_.push_back(theRecHit->id());
00208     }
00209 }
00210  for(ESRecHitCollection::const_iterator theRecHit = ESRecHits->begin();theRecHit != ESRecHits->end(); theRecHit++){
00211   theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalPreshower);
00212     const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(theRecHit->id());  
00213     math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00214     if(ROOT::Math::VectorUtil::DeltaR(myCaloJetdir,theRecHitCell_XYZPoint) < maxDeltaR){
00215       std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00216       thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00217       mySelectedDetId_.push_back(theRecHit->id());
00218     }
00219  }
00220 
00221   /*
00222   for(std::vector<CaloTowerPtr>::const_iterator i_Tower=theCaloTowers.begin();i_Tower!=theCaloTowers.end();i_Tower++){
00223     size_t numRecHits = (**i_Tower).constituentsSize();
00224     for(size_t j=0;j<numRecHits;j++) {
00225       DetId RecHitDetID=(**i_Tower).constituent(j);      
00226 
00227 
00228       DetId::Detector DetNum=RecHitDetID.det();     
00229       if(DetNum==DetId::Ecal){
00230         if((EcalSubdetector)RecHitDetID.subdetId()==EcalBarrel){
00231           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00232           EBDetId EcalID=RecHitDetID;
00233           EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);
00234           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00235           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00236           std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00237           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00238         }else if((EcalSubdetector)RecHitDetID.subdetId()==EcalEndcap){
00239           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00240           EEDetId EcalID = RecHitDetID;
00241           EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);     
00242           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00243           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00244           std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00245           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00246         }else if((EcalSubdetector)RecHitDetID.subdetId()==EcalPreshower){
00247           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalPreshower);
00248           ESDetId EcalID = RecHitDetID;
00249           ESRecHitCollection::const_iterator theRecHit=ESRecHits->find(EcalID);     
00250           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00251         
00252           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00253           std::pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00254           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00255         }        
00256       } 
00257     }
00258   }
00259   */
00260     
00261     // setting sum of Et of the isolation annulus ECAL RecHits
00262     float myIsolEcalRecHits_EtSum=0.;
00263 
00264     std::vector< std::pair<math::XYZPoint,float> > myIsolPositionAndEnergyEcalRecHits=myCaloTauElementsOperators.EcalRecHitsInAnnulus((*myleadTk).momentum(),ECALSignalConeMetric_,myECALSignalConeSize,ECALIsolConeMetric_,myECALIsolConeSize,ECALRecHit_minEt_,thePositionAndEnergyEcalRecHits);
00265     for(std::vector< std::pair<math::XYZPoint,float> >::const_iterator iEcalRecHit=myIsolPositionAndEnergyEcalRecHits.begin();iEcalRecHit!=myIsolPositionAndEnergyEcalRecHits.end();iEcalRecHit++){
00266       myIsolEcalRecHits_EtSum+=(*iEcalRecHit).second*fabs(sin((*iEcalRecHit).first.theta()));
00267     }
00268     myCaloTau.setisolationECALhitsEtSum(myIsolEcalRecHits_EtSum);    
00269   }
00270 
00271   math::XYZTLorentzVector myTks_XYZTLorentzVect(0.,0.,0.,0.);
00272   math::XYZTLorentzVector alternatLorentzVect(0.,0.,0.,0.);
00273   for(TrackRefVector::iterator iTrack=myTks.begin();iTrack!=myTks.end();iTrack++) {
00274     // build a charged pion candidate Lorentz vector from a Track
00275     math::XYZTLorentzVector iChargedPionCand_XYZTLorentzVect((**iTrack).momentum().x(),(**iTrack).momentum().y(),(**iTrack).momentum().z(),sqrt(std::pow((double)(**iTrack).momentum().r(),2)+std::pow(chargedpi_mass_,2)));
00276     myTks_XYZTLorentzVect+=iChargedPionCand_XYZTLorentzVect;
00277     alternatLorentzVect+=iChargedPionCand_XYZTLorentzVect;
00278   }
00279   myCaloTau.setTracksInvariantMass(myTks_XYZTLorentzVect.mass());
00280 
00281   std::vector<BasicClusterRef> myneutralECALBasicClusters=(*myCaloTauTagInfoRef).neutralECALBasicClusters();
00282   for(std::vector<BasicClusterRef>::const_iterator iBasicCluster=myneutralECALBasicClusters.begin();iBasicCluster!=myneutralECALBasicClusters.end();iBasicCluster++) {
00283     // build a gamma candidate Lorentz vector from a neutral ECAL BasicCluster
00284     double iGammaCand_px=(**iBasicCluster).energy()*sin((**iBasicCluster).position().theta())*cos((**iBasicCluster).position().phi());
00285     double iGammaCand_py=(**iBasicCluster).energy()*sin((**iBasicCluster).position().theta())*sin((**iBasicCluster).position().phi());
00286     double iGammaCand_pz=(**iBasicCluster).energy()*cos((**iBasicCluster).position().theta());
00287     math::XYZTLorentzVector iGammaCand_XYZTLorentzVect(iGammaCand_px,iGammaCand_py,iGammaCand_pz,(**iBasicCluster).energy());
00288     alternatLorentzVect+=iGammaCand_XYZTLorentzVect;
00289   }
00290   myCaloTau.setalternatLorentzVect(alternatLorentzVect);
00291   
00292   
00293   myCaloTau.setVertex(math::XYZPoint(myCaloTau_refInnerPosition_x,myCaloTau_refInnerPosition_y,myCaloTau_refInnerPosition_z));
00294     
00295   // setting Et of the highest Et HCAL CaloTower
00296   double mymaxEtHCALtower_Et=0.; 
00297   for(unsigned int iTower=0;iTower<myCaloTowers.size();iTower++){
00298     if((*myCaloTowers[iTower]).hadEt()>=mymaxEtHCALtower_Et) mymaxEtHCALtower_Et=(*myCaloTowers[iTower]).hadEt();
00299   }
00300   myCaloTau.setmaximumHCALhitEt(mymaxEtHCALtower_Et);
00301 
00302   return myCaloTau;  
00303 }
00304 
00305 std::vector<CaloTowerDetId> CaloRecoTauAlgorithm::getCaloTowerneighbourDetIds(const CaloSubdetectorGeometry* myCaloSubdetectorGeometry,CaloTowerDetId myCaloTowerDetId){
00306   CaloTowerTopology myCaloTowerTopology;
00307   std::vector<CaloTowerDetId> myCaloTowerneighbourDetIds;
00308   std::vector<DetId> northDetIds=myCaloTowerTopology.north(myCaloTowerDetId);
00309   std::vector<DetId> westDetIds=myCaloTowerTopology.west(myCaloTowerDetId);
00310   std::vector<DetId> northwestDetIds,southwestDetIds;
00311   if (westDetIds.size()>0){
00312     northwestDetIds=myCaloTowerTopology.north(westDetIds[0]);
00313     southwestDetIds=myCaloTowerTopology.south(westDetIds[(int)westDetIds.size()-1]);
00314   }
00315   std::vector<DetId> southDetIds=myCaloTowerTopology.south(myCaloTowerDetId);
00316   std::vector<DetId> eastDetIds=myCaloTowerTopology.east(myCaloTowerDetId);
00317   std::vector<DetId> northeastDetIds,southeastDetIds;
00318   if (eastDetIds.size()>0){
00319     northeastDetIds=myCaloTowerTopology.north(eastDetIds[0]);
00320     southeastDetIds=myCaloTowerTopology.south(eastDetIds[(int)eastDetIds.size()-1]);
00321   }
00322   std::vector<DetId> myneighbourDetIds=northDetIds;
00323   myneighbourDetIds.insert(myneighbourDetIds.end(),westDetIds.begin(),westDetIds.end());
00324   myneighbourDetIds.insert(myneighbourDetIds.end(),northwestDetIds.begin(),northwestDetIds.end());
00325   myneighbourDetIds.insert(myneighbourDetIds.end(),southwestDetIds.begin(),southwestDetIds.end());
00326   myneighbourDetIds.insert(myneighbourDetIds.end(),southDetIds.begin(),southDetIds.end());
00327   myneighbourDetIds.insert(myneighbourDetIds.end(),eastDetIds.begin(),eastDetIds.end());
00328   myneighbourDetIds.insert(myneighbourDetIds.end(),northeastDetIds.begin(),northeastDetIds.end());
00329   myneighbourDetIds.insert(myneighbourDetIds.end(),southeastDetIds.begin(),southeastDetIds.end());
00330   for(std::vector<DetId>::const_iterator iDetId=myneighbourDetIds.begin();iDetId!=myneighbourDetIds.end();iDetId++){
00331     CaloTowerDetId iCaloTowerId(*iDetId);
00332     myCaloTowerneighbourDetIds.push_back(iCaloTowerId);
00333   }
00334   return myCaloTowerneighbourDetIds;
00335 }
00336
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoTauTag/RecoTau/src/CaloRecoTauAlgorithm.cc
