#include <CMSMidpointAlgorithm.h>

Public Types
typedef std::vector< ProtoJet * >	InternalCollection

Public Member Functions
	CMSMidpointAlgorithm ()
	Default constructor. More...

	CMSMidpointAlgorithm (double fSeedThreshold, double fConeRadius, double fConeAreaFraction, int fMaxPairSize, int fMaxIterations, double fOverlapThreshold, int fDebugLevel)

void	run (const JetReco::InputCollection &fInput, JetReco::OutputCollection *fOutput)

Private Member Functions
void	addClustersToPairs (const JetReco::InputCollection &fInput, std::vector< int > &testPair, std::vector< std::vector< int > > &pairs, std::vector< std::vector< bool > > &distanceOK, int maxClustersInPair)

void	findStableConesFromMidPoints (const JetReco::InputCollection &fInput, InternalCollection *fOutput)

void	findStableConesFromSeeds (const JetReco::InputCollection &fInput, InternalCollection *fOutput)

void	iterateCone (const JetReco::InputCollection &fInput, double startRapidity, double startPhi, double startPt, bool reduceConeSize, InternalCollection *fOutput)

void	splitAndMerge (const JetReco::InputCollection &fInput, InternalCollection fProtoJets, JetReco::OutputCollection fFinalJets)

Private Attributes
double	theConeAreaFraction

double	theConeRadius

int	theDebugLevel

int	theMaxIterations

int	theMaxPairSize

double	theOverlapThreshold

double	theSeedThreshold

Detailed Description

CMSMidpointAlgorithm is an algorithm for CMS jet reconstruction baded on the CDF midpoint algorithm.

The algorithm is documented in the proceedings of the Physics at RUN II: QCD and Weak Boson Physics Workshop: hep-ex/0005012.

The algorithm runs off of generic Candidates

Author: Robert M Harris, Fermilab

Version: 1st Version Feb. 4, 2005 Based on the CDF Midpoint Algorithm code by Matthias Toennesmann.; 2nd Version Apr. 6, 2005 Modifications toward integration in new EDM.; 3rd Version Oct. 19, 2005 Modified to work with real CaloTowers from Jeremy Mans; F.Ratnikov, Mar. 8, 2006. Work from Candidate

Id:: CMSMidpointAlgorithm.h,v 1.1 2009/08/24 14:35:59 srappocc Exp

Definition at line 27 of file CMSMidpointAlgorithm.h.

Member Typedef Documentation

typedef std::vector<ProtoJet*> CMSMidpointAlgorithm::InternalCollection

Definition at line 30 of file CMSMidpointAlgorithm.h.

Constructor & Destructor Documentation

CMSMidpointAlgorithm::CMSMidpointAlgorithm ( )

inline

Default constructor.

Definition at line 33 of file CMSMidpointAlgorithm.h.

                           :
     theSeedThreshold(3.0),
     theConeRadius(0.5),
     theConeAreaFraction(1.0),
     theMaxPairSize(2),
     theMaxIterations(100),
     theOverlapThreshold(0.75),
     theDebugLevel(0)
   { }

CMSMidpointAlgorithm::CMSMidpointAlgorithm	(	double	fSeedThreshold,
		double	fConeRadius,
		double	fConeAreaFraction,
		int	fMaxPairSize,
		int	fMaxIterations,
		double	fOverlapThreshold,
		int	fDebugLevel
	)

inline

Constructor takes as input all the values of the algorithm that the user can change.

Parameters

fSeedThreshold,:	minimum ET in GeV of an CaloTower that can seed a jet.
fTowerThreshold,:	minimum ET in GeV of an CaloTower that is included in a jet
fConeRadius,:	nominal radius of the jet in eta-phi space
fConeAreaFraction,:	multiplier to reduce the search cone area during the iteration phase. introduced by CDF in 2002 to avoid energy loss due to proto-jet migration. Actively being discussed in 2005. A value of 1.0 gives the original run II algorithm in hep-ex/0005012. CDF value of 0.25 gives new search cone of 0.5 theConeRadius during iteration, but keeps the final cone at theConeRadius for the last iteration.
fMaxPairSize,:	Maximum size of proto-jet pair, triplet, etc for defining midpoint. Both CDF and D0 use 2.
fMaxIterations,:	Maximum number of iterations before finding a stable cone.
fOverlapThreshold,:	When two proto-jets overlap, this is the merging threshold on the fraction of PT in the overlap region compared to the lower Pt Jet. If the overlapPt/lowerJetPt > theOverlapThreshold the 2 proto-jets will be merged into one final jet. if the overlapPt/lowerJetPt < theOverlapThreshold the towers in the two proto-jets will be seperated into two distinct sets of towers: two final jets. D0 has used 0.5, CDF has used both 0.5 and 0.75 in run 2, and 0.75 in run 1.
fDebugLevel,:	integer level of diagnostic printout: 0 = no printout, 1 = minimal printout, 2 = pages per event.

Definition at line 65 of file CMSMidpointAlgorithm.h.

                                                                                                 : 
     theSeedThreshold(fSeedThreshold),
     theConeRadius(fConeRadius),
     theConeAreaFraction(fConeAreaFraction),
     theMaxPairSize(fMaxPairSize),
     theMaxIterations(fMaxIterations),
     theOverlapThreshold(fOverlapThreshold),
     theDebugLevel(fDebugLevel)
   { }

Member Function Documentation

void CMSMidpointAlgorithm::addClustersToPairs	(	const JetReco::InputCollection &	fInput,
		std::vector< int > &	testPair,
		std::vector< std::vector< int > > &	pairs,
		std::vector< std::vector< bool > > &	distanceOK,
		int	maxClustersInPair
	)

private

Add proto-jets to pairs, triplets, etc, prior to finding their midpoints. Called by findStableConesFromMidPoints but public method to allow testing and studies.

Definition at line 278 of file CMSMidpointAlgorithm.cc.

 {
   // Recursively adds clusters to pairs, triplets, ... whose mid-points are then calculated.
   
   // Find StableCone number to start with (either 0 at the beginning or last element of testPair + 1).
   int nextClusterStart = 0;
   if(testPair.size())
     nextClusterStart = testPair.back() + 1;
   for(unsigned int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); ++nextCluster){
     // Is new SeedCone less than 2*_theConeRadius apart from all clusters in testPair?
     bool addCluster = true;
     for(unsigned int iCluster = 0; iCluster < testPair.size() && addCluster; ++iCluster)
       if(!distanceOK[nextCluster - 1][testPair[iCluster]])
     addCluster = false;
     if(addCluster){
       // Add it to the testPair.
       testPair.push_back(nextCluster);
       // If testPair is a pair, add it to pairs.
       if(testPair.size() > 1)
     pairs.push_back(testPair);
       // If not bigger than allowed, find more clusters within 2*theConeRadius.
       if(testPair.size() < unsigned(maxClustersInPair))
     addClustersToPairs(fInput, testPair,pairs,distanceOK,maxClustersInPair);
       // All combinations containing testPair found. Remove last element.
       testPair.pop_back();
     }
   }
 }

void CMSMidpointAlgorithm::findStableConesFromMidPoints	(	const JetReco::InputCollection &	fInput,
		InternalCollection *	fOutput
	)

private

Add to the list of proto-jets the list of midpoints. Called by run but public method to allow testing and studies.

Definition at line 236 of file CMSMidpointAlgorithm.cc.

References gather_cfg::cout, Geom::deltaR2(), ProtoJet::phi(), python.multivaluedict::sort(), and ProtoJet::y().

                                                                                                                      {
   // We take the previous list of stable protojets from seeds as input and add to it
   // those from the midpoints between proto-jet pairs, triplets, etc.
   // distanceOK[i-1][j] = Is distance between stableCones i and j (i>j) less than 2*_theConeRadius?
   vector< vector<bool> > distanceOK;         // A vector of vectors
   distanceOK.resize(stableCones->size() - 1); // Set the outer vector size, num protojets - 1
   for(unsigned int nCluster1 = 1; nCluster1 < stableCones->size(); ++nCluster1){  // Loop over the protojets
     distanceOK[nCluster1 - 1].resize(nCluster1);               //Set inner vector size: monotonically increasing.
     const ProtoJet* cluster1 = (*stableCones)[nCluster1];
     for(unsigned int nCluster2 = 0; nCluster2 < nCluster1; ++nCluster2){         // Loop over the other proto-jets
       const ProtoJet* cluster2 = (*stableCones)[nCluster2];
       double dR2 = deltaR2 (cluster1->y(), cluster1->phi(), cluster2->y(), cluster2->phi());
       distanceOK[nCluster1 - 1][nCluster2] = dR2 < 4*theConeRadius*theConeRadius;
     }
   }
 
   // Find all pairs (triplets, ...) of stableCones which are less than 2*theConeRadius apart from each other.
   vector< vector<int> > pairs(0);  
   vector<int> testPair(0);
   int maxClustersInPair = theMaxPairSize;  // Set maximum proto-jets to a pair or a triplet, etc.
   if(!maxClustersInPair)maxClustersInPair = stableCones->size();  // If zero, then skys the limit!
   addClustersToPairs(fInput, testPair,pairs,distanceOK,maxClustersInPair);  // Make the pairs, triplets, etc.
   
   // Loop over all combinations. Calculate MidPoint. Make midPointClusters.
   bool reduceConeSize = false;  // Note that here we keep the iteration cone size fixed.
   for(vector<vector<int> >::const_iterator iPair = pairs.begin(); iPair != pairs.end(); ++iPair) {
     const vector<int> & Pair = *iPair;
     // Calculate rapidity, phi and energy of MidPoint.
     reco::Particle::LorentzVector midPoint (0,0,0,0);
     for(vector<int>::const_iterator iPairMember = Pair.begin(); iPairMember != Pair.end(); ++iPairMember) {
       midPoint += (*stableCones)[*iPairMember]->p4();
     }
     if(theDebugLevel>=2)cout << endl << "[CMSMidpointAlgorithm] midpoint " << iPair-pairs.begin() << ": y = " << midPoint.Rapidity() << ", phi=" << midPoint.Phi() <<
               ", size=" << Pair.size() << endl;
     iterateCone(fInput, midPoint.Rapidity(),midPoint.Phi(),midPoint.e(),reduceConeSize,stableCones);
   }
   GreaterByPtPtr<ProtoJet> compJets;
   sort (stableCones->begin(), stableCones->end(), compJets);
 }

void CMSMidpointAlgorithm::findStableConesFromSeeds	(	const JetReco::InputCollection &	fInput,
		InternalCollection *	fOutput
	)

private

Find the list of proto-jets from the seeds. Called by run, but public method to allow testing and studies.

Definition at line 130 of file CMSMidpointAlgorithm.cc.

References gather_cfg::cout, and i.

                                                                                                               {
   // This dictates that the cone size will be reduced in the iterations procedure,
   // to prevent excessive cone movement, and then will be enlarged at the end of
   // the iteration procedure (ala CDF).
   bool reduceConeSize = true;  
   
   // Get the Seed Towers sorted by Et.  
   vector<InputItem> seedTowers = etOrderedCaloTowers(fInput, theSeedThreshold);  //This gets towers
 
   // Loop over all Seeds
   for(vector<InputItem>::const_iterator i = seedTowers.begin(); i != seedTowers.end(); ++i) {
     double seedEta = (*i)->eta ();
     double seedPhi = (*i)->phi (); 
 
     // Find stable cone from seed.  
     // This iterates the cone centroid, makes a proto-jet, and adds it to the list.
     if(theDebugLevel>=2) cout << endl << "[CMSMidpointAlgorithm] seed " << i-seedTowers.begin() << ": eta=" << seedEta << ", phi=" << seedPhi << endl;
     iterateCone(fInput, seedEta, seedPhi, 0, reduceConeSize, fOutput);
     
   }
 
 }

void CMSMidpointAlgorithm::iterateCone	(	const JetReco::InputCollection &	fInput,
		double	startRapidity,
		double	startPhi,
		double	startPt,
		bool	reduceConeSize,
		InternalCollection *	fOutput
	)

private

Iterate the proto-jet center until it is stable. Called by findStableConesFromSeeds and findStableConesFromMidPoints but public method to allow testing and studies.

Definition at line 155 of file CMSMidpointAlgorithm.cc.

References abs, gather_cfg::cout, ProtoJet::e(), reco::Candidate::et(), reco::Candidate::eta(), ProtoJet::getTowerList(), i, ProtoJet::p4(), ProtoJet::phi(), reco::Candidate::phi(), ProtoJet::pt(), ProtoJet::putTowers(), mathSSE::sqrt(), lumiQTWidget::t, and ProtoJet::y().

                                                         {
   //  The workhorse of the algorithm.
   //  Throws a cone around a position and iterates the cone until the centroid and energy of the clsuter is stable.
   //  Uses a reduced cone size to prevent excessive cluster drift, ala CDF.
   //  Adds unique clusters to the list of proto-jets.
   //
   int nIterations = 0;
   bool keepJet = true;
   ProtoJet* trialCone = new ProtoJet;
   double iterationtheConeRadius = theConeRadius;
   if(reduceConeSize)iterationtheConeRadius *= sqrt(theConeAreaFraction);
   while(++nIterations <= theMaxIterations + 1 && keepJet){
     
     // Last iteration uses the full cone size.  Others use reduced cone size.
     if(nIterations == theMaxIterations + 1)iterationtheConeRadius = theConeRadius;
     
     //Add all towers in cone and over threshold to the cluster
     vector<InputItem> towersInSeedCluster 
       = towersWithinCone(fInput, startRapidity, startPhi, iterationtheConeRadius);
     if(theDebugLevel>=2)cout << "[CMSMidpointAlgorithm] iter=" << nIterations << ", towers=" <<towersInSeedCluster.size();    
     
     if(towersInSeedCluster.size()<1) {  // Just in case there is an empty cone.
       keepJet = false;
       if(theDebugLevel>=2)cout << endl;
     } else{
       //Put seed cluster into trial cone 
       trialCone->putTowers(towersInSeedCluster);
       double endRapidity = trialCone->y();
       double endPhi      = trialCone->phi();
       double endPt       = trialCone->pt();
       double endE        = trialCone->e();
       if(theDebugLevel>=2)cout << ", y=" << endRapidity << ", phi=" << endPhi << ", PT=" << endPt << ", constituents=" << trialCone->getTowerList().size () << endl;
       if(nIterations <= theMaxIterations){
     // Do we have a stable cone?
     if(std::abs(endRapidity-startRapidity)<.001 && 
        std::abs(endPhi-startPhi)<.001 && 
        std::abs(endE - startE) < .001) {
       nIterations = theMaxIterations;       // If cone size is reduced, then still one more iteration.
       if(!reduceConeSize) ++nIterations; // Otherwise, this is the last iteration.
     }
     else{
       // Another iteration.
       startRapidity = endRapidity;
       startPhi      = endPhi;
       startE        = endE;
     }
       } 
       if(nIterations==theMaxIterations+1) {
     for(vector<InputItem>::const_iterator i = towersInSeedCluster.begin(); i != towersInSeedCluster.end(); ++i) {
       InputItem t = *i;
       if(theDebugLevel>=2) cout << "[CMSMidpointAlgorithm] Tower " <<
                  i-towersInSeedCluster.begin() << ": eta=" << t->eta() << 
                  ", phi=" << t->phi() << ", ET="  << t->et() << endl;
     }
       }         
     }
   }
   
   if(keepJet){  // Our trial cone is now stable
     bool identical = false;
     // Loop over proto-jets and check that our trial cone is a unique proto-jet
     for (unsigned icone = 0; icone < stableCones->size(); icone++) {
       if (trialCone->p4() == (*stableCones) [icone]->p4()) identical = true;  // This proto-jet is not unique.
     }
     if(!identical){
       stableCones->push_back(trialCone);  // Save the unique proto-jets
       trialCone = 0;
       if(theDebugLevel>=2)cout << "[CMSMidpointAlgorithm] Unique Proto-Jet Saved" << endl;
     }    
   }
   delete trialCone;
 }

void CMSMidpointAlgorithm::run	(	const JetReco::InputCollection &	fInput,
		JetReco::OutputCollection *	fOutput
	)

Runs the algorithm and returns a list of caloJets. The user declares the vector and calls this method.

Definition at line 92 of file CMSMidpointAlgorithm.cc.

References dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, and getHLTprescales::index.

 {
   if (!fOutput) {
     std::cerr << "CMSMidpointAlgorithm::run-> ERROR: no output collection" << std::endl;
   }
   if(theDebugLevel>=1)cout << "[CMSMidpointAlgorithm] Num of input constituents = " << fInput.size() << endl;
   if(theDebugLevel>=2) {
     unsigned index = 0;
     for (; index < fInput.size (); index++) {
       cout << index << " consituent p/eta/phi: " 
        << fInput[index]->p() << '/'
        << fInput[index]->eta() << '/'
        << fInput[index]->phi() << std::endl;
     }
   }
   // Find proto-jets from the seeds.
   InternalCollection protoJets;         // Initialize working container
   // vector<ProtoJet> finalJets;         // Final proto-jets container
   findStableConesFromSeeds(fInput, &protoJets);   //Find the proto-jets from the seeds
   if(theDebugLevel>=1)cout << "[CMSMidpointAlgorithm] Num proto-jets from Seeds = " << protoJets.size() << endl;
 
   // Find proto-jets from the midpoints, and add them to the list
   if(protoJets.size()>0)findStableConesFromMidPoints(fInput, &protoJets);// Add midpoints
   if(theDebugLevel>=1)cout << "[CMSMidpointAlgorithm] Num proto-jets from Seeds and Midpoints = " << protoJets.size() << endl;
 
   // Split and merge the proto-jets, assigning each tower in the protojets to one and only one final jet.
   if(protoJets.size()>0)splitAndMerge(fInput, &protoJets, fOutput);    // Split and merge 
   if(theDebugLevel>=1)cout << "[CMSMidpointAlgorithm] Num final jets = " << fOutput->size() << endl;
 
   // Make the CaloJets from the final protojets
   //  MakeCaloJet(*theCtcp, finalJets, caloJets);
 
 }

void CMSMidpointAlgorithm::splitAndMerge	(	const JetReco::InputCollection &	fInput,
		InternalCollection *	fProtoJets,
		JetReco::OutputCollection *	fFinalJets
	)

private

Split and merge the proto-jets, assigning each tower in the protojets to one and only one final jet. Called by run but public method to allow testing and studies.

Definition at line 311 of file CMSMidpointAlgorithm.cc.

References dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, Geom::deltaR2(), spr::find(), ProtoJet::getTowerList(), i, j, muon::overlap(), ProtoJet::p4(), ProtoJet::phi(), ProtoJet::pt(), ProtoJet::putTowers(), python.multivaluedict::sort(), and ProtoJet::y().

                                                                                     {
   //
   // This can use quite a bit of optimization and simplification.
   //
   
   // Debugging
   if(theDebugLevel>=2){
     // Take a look at the proto-jets we were given
     int numProtojets = stableCones->size();
     for(int i = 0; i < numProtojets ; ++i){
       const ProtoJet* icone = (*stableCones)[i];
       int numTowers = icone->getTowerList().size();
       cout << endl << "[CMSMidpointAlgorithm] ProtoJet " << i << ": PT=" << (*stableCones)[i]->pt()
        << ", y="<< icone->y()
        << ", phi="<< icone->phi()  
        << ", ntow="<< numTowers << endl;
       ProtoJet::Constituents protojetTowers = icone->getTowerList(); 
       for(int j = 0; j < numTowers; ++j){
     cout << "[CMSMidpointAlgorithm] Tower " << j << ": ET=" << protojetTowers[j]->et() 
          << ", eta="<< protojetTowers[j]->eta()
          << ", phi="<<   protojetTowers[j]->phi() << endl;
       }
     }      
   }
   
   // Start of split and merge algorithm
   bool mergingNotFinished = true;
   while(mergingNotFinished){
     
     // Sort the stable cones (highest pt first).
     GreaterByPtPtr<ProtoJet> compJets;
     sort(stableCones->begin(),stableCones->end(),compJets);
     // clean removed clusters
     //InternalCollection::const_iterator i = find (stableCones->begin(), stableCones->end(), (ProtoJet*)0);
     //    for (; i != stableCones->end(); ++i) std::cout << "CMSMidpointAlgorithm::splitAndMerge-> removing pointer " << *i << std::endl; 
     stableCones->erase (find (stableCones->begin(), stableCones->end(), (ProtoJet*)0), stableCones->end()); 
     
     // Start with the highest pt cone left in list. List changes in loop,
     // getting smaller with each iteration.
     InternalCollection::iterator stableConeIter1 = stableCones->begin();
     
     if(stableConeIter1 == stableCones->end()) {  // Stable cone list empty?
       mergingNotFinished = false;
     }
     else {
       ProtoJet* stableCone1 = *stableConeIter1;
       if (!stableCone1) {
     std::cerr << "CMSMidpointAlgorithm::splitAndMerge-> Error: stableCone1 should never be 0" << std::endl;
     continue;
       }
       bool coneNotModified = true;
       
       // Determine whether highest pt cone has an overlap with other stable cones.
       InternalCollection::iterator stableConeIter2 = stableConeIter1; 
       ++stableConeIter2;   // Iterator for 2nd highest pt cone, just past 1st cone.
       
       while(coneNotModified && stableConeIter2 != stableCones->end()){
     ProtoJet* stableCone2 = *stableConeIter2;
     if (!stableCone2) {
       std::cerr << "CMSMidpointAlgorithm::splitAndMerge-> Error: stableCone2 should never be 0" << std::endl;
       continue;
     }
     
     // Calculate overlap of the two cones.
     bool overlap = false;
     vector<InputItem> overlapTowers;  // Make a list to hold the overlap towers
     //cout << "1st cone num towers=" << stableCone1->getTowerList().size() << endl;
     //int numTowers1=0;
     
     //Loop over towers in higher Pt cone
     for(ProtoJet::Constituents::const_iterator towerIter1 = stableCone1->getTowerList().begin();
         towerIter1 != stableCone1->getTowerList().end();
         ++towerIter1){
       //cout << "1st cone tower " << numTowers1 << endl;
       //++numTowers1;
       bool isInCone2 = false;
       //cout << "2nd cone num towers=" << stableCone2->getTowerList().size() << endl;
       //int numTowers2=0;
       
       // Loop over towers in lower Pt cone
       for(ProtoJet::Constituents::const_iterator towerIter2 = stableCone2->getTowerList().begin();
           towerIter2 != stableCone2->getTowerList().end();
           ++towerIter2) {
         //cout << "2st cone tower " << numTowers2 << endl;
         //++numTowers2;
         
         // Check if towers are the same by checking for unique eta, phi and energy values.
         // Will want to replace this with checking the tower index when available.
         if(sameTower(*towerIter1, *towerIter2)) {
           isInCone2 = true;   //The tower is in both cones
           //cout << "Merging found overlap tower: eta=" << (*towerIter1)->eta << ", phi=" << (*towerIter1)->phi << " in 1st protojet "  << 
           //                             " and eta=" << (*towerIter2)->eta << ", phi=" << (*towerIter2)->phi << " in 2nd protojet "  << endl;
         }
       }
       if(isInCone2){
         overlapTowers.push_back(*towerIter1);  //Add tower in both cones to the overlap list
         overlap = true;
       }          
     }
     if(overlap){   
       // non-empty overlap.  Decide on splitting or merging.
       
       // Make a proto-jet with the overlap towers so we can calculate things for the overlap
       ProtoJet overlap;
       overlap.putTowers(overlapTowers);
       coneNotModified = false;
       
       // Compare the overlap pt with the overlap fractcion threshold times the lower jet pt.
       if(overlap.pt() >= theOverlapThreshold*stableCone2->pt()){
         
         // Merge the two cones.
         // Get a copy of the list of towers in higher Pt proto-jet 
         ProtoJet::Constituents stableCone1Towers = stableCone1->getTowerList(); 
         
         //Loop over the list of towers lower Pt jet
         for(ProtoJet::Constituents::const_iterator towerIter2 = stableCone2->getTowerList().begin();
         towerIter2 != stableCone2->getTowerList().end();
         ++towerIter2){
           bool isInOverlap = false;
           
           //Check if that tower is in the overlap region
           for(vector<InputItem>::iterator overlapTowerIter = overlapTowers.begin();
           overlapTowerIter != overlapTowers.end();
           ++overlapTowerIter){
         // Check if towers are the same by checking for unique eta, phi and energy values.
         // Will want to replace this with checking the tower index when available.
         if(sameTower(*overlapTowerIter, *towerIter2))   isInOverlap = true;
           }
           //Add  non-overlap tower from proto-jet 2 into the proto-jet 1 tower list.
           if(!isInOverlap)stableCone1Towers.push_back(*towerIter2);  
         }
         
         if(theDebugLevel>=2)cout << endl << "[CMSMidpointAlgorithm] Merging: 1st Proto-jet grows: "
                   " y=" << stableCone1->y() << 
                   ", phi=" << stableCone1->phi() << 
                   " increases from " << stableCone1->getTowerList().size() << 
                   " to "  << stableCone1Towers.size() << " towers." << endl;
         
         // Put the new expanded list of towers into the first proto-jet
         stableCone1->putTowers(stableCone1Towers);
         
         if(theDebugLevel>=2)cout << "[CMSMidpointAlgorithm] Merging: 1st protojet now at y=" << stableCone1->y() <<
                   ", phi=" << stableCone1->phi() << endl;
         
         if(theDebugLevel>=2)cout << "[CMSMidpointAlgorithm] Merging: 2nd Proto-jet removed:" 
                   " y=" << stableCone2->y() << 
                   ", phi=" << stableCone2->phi() << endl;
         
         // Remove the second proto-jet.
         delete *stableConeIter2;
         *stableConeIter2 = 0;
         
       }
       else{
         // Split the two proto-jets.
         
         // Create lists of towers to remove from each proto-jet
         vector<InputItem> removeFromCone1,removeFromCone2;
         
         // Which tower goes where?
         // Loop over the overlap towers
         for(vector<InputItem>::iterator towerIter = overlapTowers.begin();
         towerIter != overlapTowers.end();
         ++towerIter){
           double dR2Jet1 = deltaR2 ((*towerIter)->p4().Rapidity(), (*towerIter)->phi(), 
                       stableCone1->y(), stableCone1->phi()); 
           // Calculate distance from proto-jet 2.
           double dR2Jet2 = deltaR2 ((*towerIter)->p4().Rapidity(), (*towerIter)->phi(),
                       stableCone2->y(), stableCone2->phi()); 
           
           if(dR2Jet1 < dR2Jet2){
         // Tower is closer to proto-jet 1. To be removed from proto-jet 2.
         removeFromCone2.push_back(*towerIter);
           }
           else {
         // Tower is closer to proto-jet 2. To be removed from proto-jet 1.
         removeFromCone1.push_back(*towerIter);
           }
         }
         // Remove towers in the overlap region from the cones to which they have the larger distance.
         
         // Remove towers from proto-jet 1.
         vector<InputItem> towerList1 (stableCone1->getTowerList().begin(), stableCone1->getTowerList().end()); 
 
         // Loop over towers in remove list
         for(vector<InputItem>::iterator towerIter = removeFromCone1.begin();
         towerIter != removeFromCone1.end();
         ++towerIter) {
           // Loop over towers in protojet
           for(vector<InputItem>::iterator towerIter1 = towerList1.begin(); towerIter1 != towerList1.end(); ++towerIter1) {
         // Check if they are equal
         if(sameTower(*towerIter, *towerIter1)) {
           // Remove the tower
           towerList1.erase(towerIter1);
           break;
         }
           }
         }
         
         if(theDebugLevel>=2)cout << endl << "[CMSMidpointAlgorithm] Splitting: 1st Proto-jet  shrinks: y=" << 
                   stableCone1->y() << 
                   ", phi=" << stableCone1->phi() << 
                   " decreases from" << stableCone1->getTowerList().size() << 
                   " to "  << towerList1.size() << " towers." << endl;
         
         //Put the new reduced list of towers into proto-jet 1.
         stableCone1->putTowers(towerList1); 
         // Remove towers from cone 2.
         vector<InputItem> towerList2 (stableCone2->getTowerList().begin(), stableCone2->getTowerList().end()); 
         
         // Loop over towers in remove list
         for(vector<InputItem>::iterator towerIter = removeFromCone2.begin();
         towerIter != removeFromCone2.end();
         ++towerIter) {
           // Loop over towers in protojet
           for(vector<InputItem>::iterator towerIter2 = towerList2.begin(); towerIter2 != towerList2.end(); ++towerIter2){
         // Check if they are equal
         if(sameTower(*towerIter, *towerIter2)) {
           // Remove the tower
           towerList2.erase(towerIter2);
           break;
         }
           }
         }
         
         if(theDebugLevel>=2)cout << "[CMSMidpointAlgorithm] Splitting: 2nd Proto-jet shrinks: y=" <<
                   stableCone2->y() << 
                   ", phi=" << stableCone2->phi() << 
                   " decreases from" << stableCone2->getTowerList().size() << 
                   " to "  << towerList2.size() << " towers." << endl;
         
         //Put the new reduced list of towers into proto-jet 2.
         stableCone2->putTowers(towerList2); 
       }
     }
     else {
       if(theDebugLevel>=2)cout << endl << 
                 "[CMSMidpointAlgorithm] no overlap between 1st protojet at  y=" << stableCone1->y() << 
                 ", phi=" << stableCone1->phi() << 
                 " and 2nd protojet at  y=" << stableCone2->y() << 
                 ", phi=" << stableCone2->phi() <<endl;
     }
     
     ++stableConeIter2;  //Increment iterator to the next highest Pt protojet
       }
       if(coneNotModified){
     
     if(theDebugLevel>=2)cout << 
                   "[CMSMidpointAlgorithm] Saving: Proto-jet  at y=" << stableCone1->y() << 
                   ", phi=" << stableCone1->phi() <<  " has no overlap" << endl;
     
     fFinalJets->push_back(ProtoJet (stableCone1->p4(), stableCone1->getTowerList()));
     delete *stableConeIter1;
     *stableConeIter1 = 0;
       }
     }
   }
   
   GreaterByPt<ProtoJet> compJets;
   sort(fFinalJets->begin(),fFinalJets->end(),compJets);
 }

Member Data Documentation

double CMSMidpointAlgorithm::theConeAreaFraction

private

Definition at line 109 of file CMSMidpointAlgorithm.h.

double CMSMidpointAlgorithm::theConeRadius

private

Definition at line 108 of file CMSMidpointAlgorithm.h.

int CMSMidpointAlgorithm::theDebugLevel

private

Definition at line 113 of file CMSMidpointAlgorithm.h.

int CMSMidpointAlgorithm::theMaxIterations

private

Definition at line 111 of file CMSMidpointAlgorithm.h.

int CMSMidpointAlgorithm::theMaxPairSize

private

Definition at line 110 of file CMSMidpointAlgorithm.h.

double CMSMidpointAlgorithm::theOverlapThreshold

private

Definition at line 112 of file CMSMidpointAlgorithm.h.

double CMSMidpointAlgorithm::theSeedThreshold

private

Definition at line 107 of file CMSMidpointAlgorithm.h.

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation