// ================================================
 // JetMaker class Implementation
 // ================================================
 //
 // THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
 //
 // Namespace Atlfast::
 //
 #include "AtlfastAlgs/JetMaker.h"
 #include "AtlfastAlgs/JetSmearer.h"
 #include "AtlfastAlgs/MissingMomentum.h"
 #include "AtlfastAlgs/GlobalEventData.h"
 
 #include "AtlfastEvent/ContainerAssocsDispatcher.h"
 #include "AtlfastEvent/Jet.h"
 #include "AtlfastEvent/ICell.h"
 #include "AtlfastEvent/CollectionDefs.h"
 #include "AtlfastEvent/ReconstructedParticle.h"
 #include "AtlfastEvent/ParticleCodes.h"
 #include "AtlfastEvent/MsgStreamDefs.h"
 #include "AtlfastEvent/ICluster.h"
 #include "AtlfastEvent/MCparticleCollection.h"
 #include "AtlfastEvent/SimpleKinematic.h"
 #include "AtlfastEvent/Phi.h"
 
 #include "AtlfastUtils/TesIO.h"
 #include "AtlfastUtils/HeaderPrinter.h"
 #include "AtlfastUtils/FunctionObjects.h"
 #include "AtlfastUtils/IsAssociated.h"
 #include "AtlfastUtils/SmearCell.h"
 
 #include "TruthHelper/GenIMCselector.h"
 #include "TruthHelper/IsGenType.h"
 #include "TruthHelper/IsGenStable.h"
 #include "TruthHelper/NCutter.h"
 
 #include <cmath> 
 #include <algorithm>
 #include <numeric>
 #include <assert.h> 
 
 // Gaudi includes
 #include "GaudiKernel/DataSvc.h"
 #include "GaudiKernel/ISvcLocator.h"
 #include "GaudiKernel/MsgStream.h"
 
 // CLHEP,HepMC
 #include "GeneratorObjects/McEventCollection.h"
 #include "HepMC/GenEvent.h"
 #include "HepMC/GenVertex.h"
 
 #include "CLHEP/Units/SystemOfUnits.h"
 
   /** @brief Makes Atlfast::Jets from acceptable Atlfast::Clusters
    *
    * Acceptability is based on ET and eta cuts. Once accepted as a jet,
    * the energy is smeared. Initial labelling is performed by searching
    * in the truth block for proximity to b, c and hadronic taus.
    *
    * @image html JetMaker.jpg "JetMaker sequence"
    */
 
 namespace Atlfast{
 
 //--------------------------------
 // Constructors and destructors
 //--------------------------------
   using std::abs;
   using std::sort;
   using std::partition;
   
   JetMaker::JetMaker 
   ( const std::string& name, ISvcLocator* pSvcLocator ) 
     : Algorithm( name, pSvcLocator ),
       m_bSelector(0),
       m_cSelector(0),
       m_tauSelector(0),
       m_smearer(0),
       m_cellSmearer(0),
       m_tesIO(0)
   {
     
     // Setting the parameter defaults.
     m_minPT            = 10*GeV;     // Min jet trans momentum
     m_maxEta           = 5.0;          //Max jet eta 
     m_doSmearing       = true;       //smearing is done
     m_rconeb           = 0.400;
     m_rconef           = 0.400;
     
     m_bPtMin           = 5.0*GeV;    //min pt for bquark to label jet
     m_bMaxDeltaR       = 0.3;        //max dR for bquark to label jet
     
     m_cPtMin           = 5.0*GeV;    //min pt for cquark to label jet
     m_cMaxDeltaR       = 0.3;
     
     m_tauPtMin         = 10.0*GeV;
     m_tauMaxDeltaR     = 0.3;
     
     m_etaTagMax        = 2.5;
 
     m_tauJetPtRatio    = 0;
     
     m_adjustMissETForIsolation = true;
 
     // Default paths for entities in the TES
     m_inputLocation             = "/Event/AtlfastClusters";
     m_outputLocation            = "/Event/AtlfastJets";
     m_muonLocation              = "/Event/AtlfastNonIsolatedMuons";
     m_unusedCellLocation        = "/Event/AtlfastUnusedCells";
     m_missingMomentumLocation   = "/Event/AtlfastMissingMomentum";
     m_isolatedElectronLocation  = "/Event/AtlfastIsolatedElectrons";
     m_isolatedPhotonLocation    = "/Event/AtlfastIsolatedPhotons";
     
     // This is how you declare the paramemters to Gaudi so that
     // they can be over-written via the job options file
     declareProperty( "MinimumPT",                m_minPT ) ;
     declareProperty( "MaximumEta",               m_maxEta ) ;
     declareProperty( "DoSmearing",               m_doSmearing ); 
     declareProperty( "RconeB",                   m_rconeb );
     declareProperty( "RconeF",                   m_rconef );
     declareProperty( "bPtMin",                   m_bPtMin );
     declareProperty( "bMaxDeltaR",               m_bMaxDeltaR );
     declareProperty( "cPtMin",                   m_cPtMin );
     declareProperty( "cMaxDeltaR",               m_cMaxDeltaR );
     declareProperty( "tauPtMin",                 m_tauPtMin );
     declareProperty( "tauMaxDeltaR",             m_tauMaxDeltaR );
     declareProperty( "etaTagMax",                m_etaTagMax );
     declareProperty( "tauJetPtRatio",            m_tauJetPtRatio);
     declareProperty( "adjustMissETForIsolation", m_adjustMissETForIsolation);
 
     declareProperty( "InputLocation",            m_inputLocation ) ;
     declareProperty( "OutputLocation",           m_outputLocation ) ;
     declareProperty( "UnusedCellLocation",       m_unusedCellLocation ) ;
     declareProperty( "MuonLocation",             m_muonLocation );
     declareProperty( "MissingMomentumLocation",  m_missingMomentumLocation );
     declareProperty( "IsolatedElectronLocation", m_isolatedElectronLocation );
     declareProperty( "IsolatedPhotonLocation",   m_isolatedPhotonLocation );
 
   }
   
   JetMaker::~JetMaker() {
     
     MsgStream log( messageService(), name() ) ;
     log << MSG::INFO << "Destructor Called" << endreq;
     
     if (m_smearer) {
       log << MSG::INFO << "Deleting jet Smearer" << endreq;
       delete m_smearer;
     }
     if (m_cellSmearer) {
       log << MSG::INFO << "Deleting Cell Smearer" << endreq;
       delete m_cellSmearer;
     }
     if (m_tesIO) {
       delete m_tesIO;
     }
   } 
   
   //---------------------------------
   // initialise() 
   //---------------------------------
   
   StatusCode JetMaker::initialize(){
     MsgStream log( messageService(), name() ) ;
     
     log << MSG::DEBUG << "instantiating a JetSmearer" << endreq;
     
     
     //get the Global Event Data using singleton pattern
     GlobalEventData* ged = GlobalEventData::Instance();
     int lumi                  = ged->lumi();
     int randSeed              = ged->randSeed() ;
     m_barrelForwardEta        = ged->barrelForwardEta();
     m_mcLocation              = ged->mcLocation();
 
     // Pass this into GlobalEventData for EventHeaderMaker to read later
     ged->set_adjustMissEtForIsolation(m_adjustMissETForIsolation);
 
     //    m_tesIO = new TesIO(eventDataService());
     m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
     
     if(m_doSmearing){
       m_smearer                 = new JetSmearer(randSeed, 
                                                  lumi,
                                                  m_rconeb, 
                                                  m_rconef, 
                                                  m_barrelForwardEta);    
       m_cellSmearer             = new JetSmearer(randSeed, 
                                                  lumi,
                                                  0., 
                                                  0., 
                                                  m_barrelForwardEta);    
     } else {
       m_smearer = 0;
       m_cellSmearer = 0;
     }
 
     const int BQUARK(ParticleCodes::BQUARK);
     const int CQUARK(ParticleCodes::CQUARK);
 
     //.........................................................
     // Construct objects to select b, c and tau quarks with kine conditions
     
     TruthHelper::GenIMCselector* selector;
     
     std::vector<TruthHelper::GenIMCselector*> bSelectors;
     selector = new HepMC_helper::SelectJetTag(BQUARK, m_bPtMin, m_etaTagMax);
     bSelectors.push_back( selector );
       m_bSelector = new TruthHelper::NCutter(bSelectors);
 
     std::vector<TruthHelper::GenIMCselector*> cSelectors;
     selector = new HepMC_helper::SelectJetTag(CQUARK, m_cPtMin, m_etaTagMax);
     cSelectors.push_back( selector );
      m_cSelector = new TruthHelper::NCutter(cSelectors);
 
     std::vector<TruthHelper::GenIMCselector*> tauSelectors;
     selector = new HepMC_helper::SelectTauTag(m_tauPtMin, m_etaTagMax);
     tauSelectors.push_back( selector );
     m_tauSelector = new TruthHelper::NCutter(tauSelectors);
 
     //.........................................................
     // Tidy up selectors
     std::vector<TruthHelper::GenIMCselector*>::iterator iter;
     for(iter=bSelectors.begin();   iter!=bSelectors.end();   delete *iter, ++iter);
     for(iter=cSelectors.begin();   iter!=cSelectors.end();   delete *iter, ++iter);
     for(iter=tauSelectors.begin(); iter!=tauSelectors.end(); delete *iter, ++iter);
     
     
 
     HeaderPrinter hp("Atlfast JetMaker:", log);
     
     hp.add("TES Locations:              ");
     hp.add(" Muons from                 ", m_muonLocation); 
     hp.add(" Clusters from              ", m_inputLocation); 
     hp.add(" MC Input Location          ", m_mcLocation);
     hp.add(" Jets to                    ", m_outputLocation);
     hp.add(" unused cell,cluster sum to ", m_missingMomentumLocation);
     hp.add("Smearing on                 ", m_doSmearing);
     hp.add("Cluster min Pt              ", m_minPT);
     hp.add("Cluster max Eta             ", m_maxEta);           
     hp.add("Luminosity                  ", lumi);
     hp.add("Initial random seed         ", randSeed);      
     hp.add("End cap Cone size           ", m_rconef);
     hp.add("Barrel  Cone size           ", m_rconeb);
     hp.add("Eta of start of end cap     ", m_barrelForwardEta);
     hp.add("Parameters for labels       ", "       ");
     hp.add(" Max eta                    ", m_etaTagMax);
     hp.add(" b-quark min pt             ", m_bPtMin);
     hp.add(" b-quark max DeltaR         ", m_bMaxDeltaR);
     hp.add(" c-quark min pt             ", m_cPtMin);
     hp.add(" c-quark max DeltaR         ", m_cMaxDeltaR);
     hp.add(" tau lep min pt             ", m_tauPtMin);
     hp.add(" tau lep max DeltaR         ", m_tauMaxDeltaR);
     hp.add(" tau/Jet Pt Ratio           ", m_tauJetPtRatio);
     hp.print();
     
 
 
     return StatusCode::SUCCESS ;
   }
   
   //---------------------------------
   // finalise() 
   //---------------------------------
   
   StatusCode JetMaker::finalize(){
 
     MsgStream log( messageService(), name() ) ;
     log << MSG::INFO << "finalizing" << endreq;    
 
     if (m_bSelector)    delete m_bSelector;
     if (m_cSelector)    delete m_cSelector;
     if (m_tauSelector)  delete m_tauSelector;
 
     return StatusCode::SUCCESS ;
   }
   
   
   //----------------------------------------------
   // execute() method called once per event
   //----------------------------------------------
   //
   //  This algorithm creates smeared Jets passing cuts. 
   //  It scans the list of Clusters from the Monte Carlo truth. 
   //  If a Cluster of the right flavour is found it creates 
   //  a candidate Jet object by smearing the four vector
   //  of the Cluster and adding non-isolated muons. 
   // 
   //  It then applies kinematic criteria on the properties of the Jet
   //  Those which pass the cuts are kept.
   //  Finally all successful Jets are added to the event store.
   //
   //    //
   
   StatusCode JetMaker::execute( ){
     
     MsgStream log( messageService(), name() ) ;    
 
     std::string mess;
     // make a message logging stream
     
     //.............................
     // Make some locals stores which be used to keep pointers to
     // all of the entities from the event store which are needed by
     // this algorithm. These are all local and are typedefed in this class
     
     localClusterCollection      my_Clusters ;
     MuonCollection              my_Muons ;
     //................................
     
     //    if( ! m_tesIO->copy( my_Clusters, m_inputLocation ) ) {
     if( ! m_tesIO->copy<IClusterCollection> ( my_Clusters, m_inputLocation ) ) {
       log << MSG::INFO << "No Clusters in TES " << endreq;
     } else {
       log << MSG::DEBUG << my_Clusters.size()<<" Clusters from TES " << endreq;
     }
     // get non-isolated muons for jet construction
     
     if( ! m_tesIO->copy<ReconstructedParticleCollection>
         ( my_Muons, m_muonLocation ) ) {
       log << MSG::DEBUG << "No muons in TES " << endreq;
     }  
 
     // ......................
     // Make a container object which will store pointers to all isolated 
     // Jetss which are successfully created.  Since it is going to go 
     // into the event store, it has to be a special Gaudi type of container. 
     // This is typedefined in the entity include file
     // As it is going into the TES it needs to be created on the heap with new.
     
     JetCollection* myJets = new JetCollection ;
     
     //.........................................................
     // From each Cluster create a reconstructed Jet. 
     // If all requirements are satisfied add to the collection
     
     Jet* candidate ;
     localClusterIterator src ;
     m_sumET = 0.;
     HepLorentzVector missingMomentum(0.,0.,0.,0.);
     
     MCparticleCollection  mcParticles_b ;
     MCparticleCollection  mcParticles_c ;
     MCparticleCollection  mcParticles_tau ;
 
     //.........................................................
     // Get quarks used in tagging from TES using tesio
 
     TesIoStat stat = m_tesIO->getMC( mcParticles_b, m_bSelector ) ;
     if(!stat){
       log << MSG::ERROR <<"Problem getting b-quarks" << endreq; 
       return stat;
     }
     stat = m_tesIO->getMC( mcParticles_c, m_cSelector ) ;
     if(!stat){
       log << MSG::ERROR <<"Problem getting c-quarks" << endreq; 
       return stat;
     }
     stat = m_tesIO->getMC( mcParticles_tau, m_tauSelector ) ;
     if(!stat){
       log << MSG::ERROR <<"Problem getting taus" << endreq; 
       return stat;
     }
  
     for( src = my_Clusters.begin() ; src != my_Clusters.end() ; ++src ) { 
       //
       //      ReconstructedParticle rp;
       //      std::vector<IAOO*> v;
       //      v.push_back(*src);
       //      AssocTypeRecoverer atr(v.begin(), v.end());
       //      if ((atr.typeVector(rp))->size() == 0){
       if ( !IsAssociated<ReconstructedParticle>( *src ) ){
         // Make a new jet
         candidate = this->create( *src );
         
         //Make a temp HepLorentzVector from candidates 4-mom
         //before muon is added in
         HepLorentzVector temp4Vec = candidate->momentum();
         
         // add isolated muons
         this->addMuons(candidate, my_Muons);
         
         
         // add to missing momentum vector if not used in a jet.
         // remember, these are smeared momenta if smearing is enabled
         if( this->isAcceptable( candidate ) ) {    
           
           // if the jet is acceptable label the jet
           this->label(candidate,mcParticles_b,mcParticles_c,mcParticles_tau);
           
           log << MSG::DEBUG << "Pushing back a Jet" << endreq;
           myJets->push_back( candidate ) ; 
           log<<MSG::DEBUG<<"jet "<<*candidate<<endreq;
         }else {
           log << MSG::DEBUG << "Jet fails cuts" << endreq;
           //add in temp4Vec to avoid muon double counting
           missingMomentum += temp4Vec;
           m_sumET += temp4Vec.perp();
           delete candidate;
         }
       }else{
         log << MSG::DEBUG << "Cluster is associated, so ignored" << endreq;
       }
     }
     
     log<<MSG::DEBUG<<"Missing momentum: clusters "<<missingMomentum<<endreq;
     missingMomentum += addCells(log);
     log << MSG::DEBUG << "sumET after adding cells: " << m_sumET << endreq;
     log<<MSG::DEBUG<<"Missing momentum: cells+clstrs "
        <<missingMomentum<<endreq;
 
     //......................................
     // Register any Jets which have been successfilly created in the 
     // transient event store. Allow modification further along the line
     // for tagging routines.
     log<<MSG::DEBUG<<"Storing "<< myJets->size() <<" Jets"<<endreq;
     stat = m_tesIO->store( myJets, m_outputLocation, true) ;
     if(!stat){
       log<<MSG::ERROR<<"Could not store jets"<<endreq;
       return stat;
     }
 
 
     MissingMomentum* mm = new MissingMomentum(missingMomentum,m_sumET);
     stat = m_tesIO->store(mm,m_missingMomentumLocation ) ;
     std::string message = stat? "Stored missing mom":"Failed to store missing mom";
     log<<MSG::DEBUG<<message<<endreq;
   
 
     return StatusCode::SUCCESS;
   }
   
   
   
   //----------------------------------
   // create()
   //----------------------------------
   
   // Takes a single Cluster and creates a Jet
   // according to the desired criteria
   //
   // Note that we must NEW these jets so they can go to the TES
   // and if we decide that we do not want them, we must DELETE them
   // Once they are put to the TES, they are no longer our responsibility to delete
   
   Jet* JetMaker::create ( ICluster* src )
   {
     Jet* jet;
     MsgStream log( messageService(), name() ) ;
     HepLorentzVector vec = src->momentum();
     
     if (m_doSmearing) 
       {
         // Ask our Smearer make a smeared 4-vector
         vec = m_smearer->smear(vec);      
       }
     
     jet = new Jet(vec, *src);
     log << MSG::DEBUG << "Unsmeared Cluster: " << *src << endreq;
     log << MSG::DEBUG << "Smeared Jet      : " << jet << endreq ;
     
     //IAOO* js = jet;
     //assert(!js->unAssociated());
     // return Jet
     return jet;
     
     
   }
   //-------------------------------------------
   // addmuons
   //-------------------------------------------
   
   // adds non-isolated muons to the jets
   
   void JetMaker::addMuons(Jet* jet, MuonCollection& muons )
   {
     
     MsgStream log( messageService(), name() ) ;
     
     if (muons.size() == 0) return;
     
     MuonIterator first = muons.begin() ;
     MuonIterator last  = muons.end() ;
     
     // Partition muons to get those that satisfy the dR condition in either 
     // eta region
     MuonIterator partitionLast = 
       partition(first,last,PartitionCondition::BelowThresholdDeltaR
       (jet,m_rconeb,m_rconef,m_barrelForwardEta));
     
     if (partitionLast == first){
       log << MSG::DEBUG << "No muons in R cone" << endreq;
       return;
     }
 
     log << MSG::DEBUG << "Adding " << partitionLast-first << " muons to Jet" << endreq;
 
     // Add muon momenta to jet
     for (MuonIterator muItr = first; muItr != partitionLast;++muItr) jet->addMuon(*muItr);
     
     // Remove muons from list
     muons.erase(first,partitionLast);
     
   }
   
   
   //-------------------------------------------
   // isAcceptable
   //------------------------------------------
   
   // Decides whether a candidate is acceptable to this Maker, i.e. whether
   // it wants to proceed and add it to its output product collection
   
   bool JetMaker::isAcceptable ( Jet* candidate)
   {
     
     MsgStream log( messageService(), name() ) ;
     
     // Apply kinematic criteria to the candidate jet
     
     if( candidate->pT() < m_minPT ) {        
       log << MSG::DEBUG << "Candidate failed pt cut: pt " << candidate->pT() 
           << " min was " << m_minPT << endreq;
       return false ;
     }
     
     if( abs(candidate->eta()) > m_maxEta ) {
       log << MSG::DEBUG << "Candidate failed max eta cut: eta " 
           << candidate->eta() 
           << " max was " << m_maxEta << endreq;
       return false ;
     }
     
     return true ;
   }
   
   //-------------------------------------------
   // label
   //------------------------------------------
   //
   
   void JetMaker::label(Jet* jet, 
                        MCparticleCollection& mcParticles_b, 
                        MCparticleCollection& mcParticles_c, 
                        MCparticleCollection& mcParticles_tau)const {
     
     MsgStream log( messageService(), name() ) ;
     
     double dR = 0.;
 
     // The getJetLabel and getTauLabel functions loop over the input particle collections
     // and return the PDGID of a particle if within maxdeltaR of the jet (0 otherwise). 
     // dR is also set to the lowest found particle-jet delta R value
 
     int bLabel = getJetLabel(mcParticles_b, jet, m_bMaxDeltaR, dR);
     jet->setdRbquark(dR);
     int cLabel = getJetLabel(mcParticles_c, jet, m_cMaxDeltaR, dR);
     jet->setdRcquark(dR);
     int tauLabel = getTauLabel( mcParticles_tau, jet, m_tauMaxDeltaR, dR);
     jet->setdRhadtau(dR);
 
     // Set jet label (order dependent)
     if(bLabel){
       jet->setPdg_id(bLabel);
       return;
     }
     if(cLabel){
       jet->setPdg_id(cLabel);
       return;
     }
     if(tauLabel){
       jet->setPdg_id(tauLabel);
       return;
     }
 
     jet->setPdg_id(0);
     return; 
 
   }
   //-------------------------------------------
   // getJetLabel
   //------------------------------------------
   int JetMaker::getJetLabel(MCparticleCollection&  mcParticles, 
                             Jet* jet,
                             double maxDeltaR,
                             double &deltaR) const {
     
     MsgStream log( messageService(), name() ) ;
 
     int pdgId    = 0;
     double drMin = 999.;
     
     //Construct an object to select IKinematics within DeltaR of the current jet
     PartitionCondition::BelowThresholdDeltaR labelSelector(jet, maxDeltaR); 
     MCparticleCollectionCIter src ;
     for( src = mcParticles.begin(); src != mcParticles.end(); ++src ){
       SimpleKinematic sk(*src);
       if(labelSelector(sk)){
         pdgId = (*src)->pdg_id();
       }      
       Phi dPhi( sk.phi() - jet->phi() );
       double dR = sqrt((dPhi*dPhi) + 
                        (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
       if (dR < drMin) drMin = dR;
     }
     deltaR = drMin;
     return pdgId;
   }
   //-------------------------------------------
   // getTauLabel
   //------------------------------------------
   int JetMaker::getTauLabel(MCparticleCollection&  mcParticles,
                             Jet* jet,
                             double maxDeltaR,
                             double &deltaR) const {
     
     MsgStream log( messageService(), name() ) ;
     
     int pdgId    = 0;
     double drMin = 999.;
 
     //Construct an object to select IKinematics within DeltaR of the current jet
     PartitionCondition::BelowThresholdDeltaR labelSelector(jet, maxDeltaR);
     MCparticleCollectionCIter src ;
     for( src = mcParticles.begin() ; src != mcParticles.end() ; ++src )
       {
         HepLorentzVector tauHlv((*src)->momentum().px(),(*src)->momentum().py()
         ,(*src)->momentum().pz(),(*src)->momentum().e());
         HepLorentzVector neuHlv;
 
         // Find Tau children
         HepMC::GenVertex*  endVertex = (*src)->end_vertex();
         if(!endVertex) continue;
 
         HepMC::GenVertex::particle_iterator firstChild = 
           endVertex->particles_begin(HepMC::children);
         
         HepMC::GenVertex::particle_iterator endChild = 
           endVertex->particles_end(HepMC::children);
         
         HepMC::GenVertex::particle_iterator thisChild = firstChild; 
         
         // Find Pt of Taus tau-neutrino
         for(; thisChild!=endChild; ++thisChild){
           if(  abs((*thisChild)->pdg_id()) == 16){
             neuHlv.setPx((*thisChild)->momentum().px());
             neuHlv.setPy((*thisChild)->momentum().py());
             neuHlv.setPz((*thisChild)->momentum().pz());
             neuHlv.setE((*thisChild)->momentum().e());
           }
         }
         double sigma;
         double sqrtene = sqrt(jet->momentum().e()/GeV);
         if ( fabs(jet->eta()) < m_barrelForwardEta )
           sigma = 0.5/sqrtene;
         else
           sigma = 1.0/sqrtene;
 
         double cutQuantity = (m_tauJetPtRatio) ? m_tauJetPtRatio : 1.0 - 2*sigma;
 
         log << MSG::DEBUG << "m_tauJetPtRatio = " << m_tauJetPtRatio
             << ", cutQuantity = " << cutQuantity << endreq;
 
         HepMC::GenParticle trueHadSystem;
         trueHadSystem.set_momentum(tauHlv-neuHlv); 
         SimpleKinematic sk(trueHadSystem);
         if (labelSelector(sk) &&
             (tauHlv-neuHlv).perp()/jet->pT() >= cutQuantity){
           Phi dPhi( sk.phi() - jet->phi() );
           deltaR = sqrt((dPhi*dPhi) +
                         (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
           pdgId = (*src)->pdg_id();
         }
         Phi dPhi( sk.phi() - jet->phi() );
         double dR = sqrt((dPhi*dPhi) + 
                          (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
         if (dR < drMin) drMin = dR;
       }
     
     deltaR = drMin;
     return pdgId;
   }   
   
   HepLorentzVector JetMaker::addCells(MsgStream& log){
     
     std::vector<const ITwoCptCell*> cellsToSmear(0);
     std::vector<const ITwoCptCell*> newCellsVector(0);
     
     getUnusedCells(cellsToSmear,log);
     
     if (m_adjustMissETForIsolation){
       getIsolationCorrectionCells(cellsToSmear,newCellsVector,m_isolatedElectronLocation,log);
       getIsolationCorrectionCells(cellsToSmear,newCellsVector,m_isolatedPhotonLocation,log);
     }
     
     log << MSG::DEBUG << "Smearing a total of " << cellsToSmear.size() << " cells" << endreq;
     log << MSG::DEBUG << "sumET before adding cells: " << m_sumET << endreq;
     
     // accumulate used to compute vectorial sum of Cell HepLorentzVectors
     // m_sumET also passed into SmearCell to compute sum of smeared cell ETs
     HepLorentzVector temp(0., 0., 0., 0.);
     HepLorentzVector sum = std::accumulate(cellsToSmear.begin(),
                                            cellsToSmear.end(),
                                            temp,
                                            SmearCell(m_cellSmearer,&m_sumET));
     //if (eHitCellPtr) delete eHitCellPtr;
 
     std::vector<const ITwoCptCell*>::iterator itccItr   = newCellsVector.begin();
     std::vector<const ITwoCptCell*>::iterator itccItrE  = newCellsVector.end();
 
     for (;itccItr!=itccItrE;++itccItr) delete *itccItr;
     newCellsVector.clear();
 
     return sum;
   }
   
   void JetMaker::getUnusedCells(std::vector<const ITwoCptCell*> &cellsToSmear, MsgStream& log){
     
     std::vector<const TwoCptCell*> unusedCells(0);
     if( !m_tesIO->copy<std::vector<const TwoCptCell*> >(unusedCells,
                                                         m_unusedCellLocation)){
       log << MSG::INFO << "No unused cells  in TES " << endreq;
       return;
     }
     
     log << MSG::DEBUG << unusedCells.size() << " unused cells in TES" << endreq;
     
     std::copy(unusedCells.begin(),unusedCells.end(),std::back_inserter(cellsToSmear));
     
     log << MSG::DEBUG  << "Added " << unusedCells.size() << " unused cells to cellToSmear" << endreq;
     
     return;
     
   }
   
   void JetMaker::getIsolationCorrectionCells(std::vector<const ITwoCptCell*> &cellsToSmear, 
                                              std::vector<const ITwoCptCell*> &newcellsVector,
                                              std::string tesAddress,
                                              MsgStream& log){
     
     
     std::vector<const ITwoCptCell*> isolationCorrectionCells(0);
     
     const ReconstructedParticleCollection* dh(0);
     if( !m_tesIO->getDH(dh, tesAddress) ) throw "Error in numberInList";
     
     ReconstructedParticleCollection::const_iterator iter = dh->begin();
     
     for (; iter != dh->end(); ++iter){
       
       const HepMC::GenParticle* gp = (*iter)->truth();
 
       // Get associations of ReconstructedParticle associations 
       const IAOO* iaooRP = (*iter);
       TypeVisitor tv =
         ContainerAssocsDispatcher( iaooRP->begin(), iaooRP->end(), TypeVisitor() );
       
       // Select the TwoCptCells
       std::vector<const TwoCptCell*> tccv = tv.typeVector(TwoCptCell());
       log << MSG::DEBUG << "TypeVisitor found " << tccv.size() 
           << " TwoCptCells in this cluster" << endreq;
       
       std::vector<const TwoCptCell*>::iterator tccvi;
 
       for ( tccvi = tccv.begin(); tccvi != tccv.end();){
         const ICell* ic = (*tccvi);
 
         // See if this TwoCptCell contains the GenParticle that made the
         // original ReconstructedParticle
         std::vector<const GenParticle*> gpv = ic->particles();
         std::vector<const GenParticle*>::iterator gpi;
         gpi = std::find(gpv.begin(),gpv.end(),gp);
 
         if ( gpi == gpv.end() ) ++tccvi;     // If not, go onto next TwoCptCell
         else{                                // It's in this one
           log << MSG::DEBUG << "GenParticle found in this cell, one of " 
               << gpv.size() << " GenParticles" << endreq;
           log << MSG::DEBUG << "GenParticle momentum = " << (*gpi)->momentum() << endreq;
           log << MSG::DEBUG << "Cell momentum = " << ic->momentum() << endreq;
 
           if ( gpv.size() > 1){
             // Make a new cell with pT reduced by that of the originating GenParticle
             ITwoCptCell* newcell_itcc = (*tccvi)->cloneITCC();
             newcell_itcc->resetCell();
             ICell* newcell = newcell_itcc;
             double newpT = ic->momentum().perp() - (*gpi)->momentum().perp();
 //          newpT = ( fabs(newpT) < 1e-10 ) ? 0 : newpT; // problem with precision
 //          newcell->newHit( newpT );
 
 
             // The above commented out lines have been updated with the lines below.
             // This is because the calculation of the extra energy (not belonging to the GenParticle)
             // namely:
             //         double newpT = ic->momentum().perp() - (*gpi)->momentum().perp();
             // crashes and is not valid if a detector effect like scaling energy is applied to a cell
             // since the cell energy (ic->momentum().perp()) will be scaled (possible lower that original)
             // but the GenParticle energy ((*gpi)->momentum().perp()) will not be, resulting in a new pt
             // being set that is <0 for this now correction cell. This causes a crash as might be expected.
             // Solution here is to have the subtraction performed by the cell. This way the scaling is done
             // on the fly by the cell which knows its scale factor (which incidently we don't know here and
             // so couldn't scale the GenParticle properly)
             if( fabs(newpT) < 1e-10 ){ newcell->newHit( 0. ); }// there seems to have bee a problem with precision
             else{                                              // so keep this check just in case.
               newcell->setPtSum( ic->momentum().perp() ); //add back the original momentum not with newHit as that would scale it at the cell energy is already scaled... thats the problem
               newcell->newHit( -1*((*gpi)->momentum().perp()) ); //subtract the GenParticle momentum
             }
             
             // Add the new, adjusted hit cell
             isolationCorrectionCells.push_back(newcell_itcc);
             newcellsVector.push_back(newcell_itcc);
           }
           
           // Remove the pointer to the hit cell
           tccvi = tccv.erase(tccvi);
           
         }
       }
       // After adjustment, copy all cells into isolationCorrectionCells
       std::copy(tccv.begin(),tccv.end(),std::back_inserter(isolationCorrectionCells));
     }
     
     log << MSG::DEBUG << "Added " << isolationCorrectionCells.size() 
         << " isolation correction cells to cellsToSmear" << endreq;
    
     // Add isolationCorrectionCells cells to vector of cells for smearing
     std::copy(isolationCorrectionCells.begin(),
               isolationCorrectionCells.end(),
               std::back_inserter(cellsToSmear));
     
     return;
     
   }
   
 } // end namespace bracket
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.