[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderAlgoUA1Hybrid.cxx

//SARAH'S HYBRID PERSONAL COPY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

//THIS IS SARAH'S REVISED UA1 CODE WITH CHANGES FOR ETA/PHI ITERATION INCLUDED!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//THIS Also includes summing ALL cells in the jetcone towards the jet energy NOT just those above threshold!!!!!


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
 
$Log$
Revision 1.3  2003/06/18 17:00:25  horner
Cleaned up more cout statements - Mark

Revision 1.2  2003/06/18 00:20:07  horner
Removed cout statement - Mark


*/


//*--Author: Sarah Blyth (LBL)
//*--Based on UA1 jet algorithm from LUND JETSET called from EMC-erj

#include "TTask.h"
#include "AliEMCALJetFinderInput.h"
#include "AliEMCALJetFinderOutput.h"
#include "AliEMCALJetFinderAlgo.h"
#include "AliEMCALJetFinderAlgoUA1Hybrid.h"
#include "AliEMCALJetFinderAlgoUA1Unit.h"
#include "AliEMCALGeometry.h"
#include "AliEMCAL.h"
#include "AliEMCALGetter.h"
#include "AliEMCALDigit.h"
#include "TParticle.h"
#include "AliRun.h"
#include "AliEMCALJet.h"
#include "TMath.h"


ClassImp(AliEMCALJetFinderAlgoUA1Hybrid)

  AliEMCALJetFinderAlgoUA1Hybrid::AliEMCALJetFinderAlgoUA1Hybrid()
{
  //Default constructor
if (fDebug>0) Info("AliEMCALJetFinderAlgoUA1Hybrid","Beginning Default Constructor");

  fNumIter           = 0;
  fNumUnits          = 13824;     //Number of towers in EMCAL
  fESeed             = 5.0;       //Default value
  fConeRad           = 0.5;       //Default value
  fJetEMin           = 10.0;      //Default value
  fEtMin             = 0.28;      //Default value
  fMinMove           = 0.05;      //From original UA1 JetFinder
  fMaxMove           = 0.15;      //From original UA1 JetFinder
  fBGMaxMove         = 0.035;     //From original UA1 JetFinder
  fPtCut             = 0;         
  fHadCorr           = 0;       
  fEBGTotal          = 1.0;       //Set to 1 so that no div by zero in first FindJets() loop
  fEBGTotalOld       = 0.0;
  fEBGAve            = 0.0;
  fEnergy            = 0.0;
  fJetEta            = 0.0;
  fJetPhi            = 0.0;
  fEtaInit           = 0.0;
  fPhiInit           = 0.0;
  fEtaB              = 0.0;
  fPhiB              = 0.0;
  fJetESum           = 0.0;
  fJetEtaSum         = 0.0;
  fJetPhiSum         = 0.0;
  fDEta              = 0.0;
  fDPhi              = 0.0;
  fDistP             = 0.0;
  fDistI             = 0.0;
  fTempE             = 0.0;
  fRad               = 2.0;      //Set to 2 to start 
  fNumInCone         = 0;
  fNumJets           = 0;
  fArrayInitialised  = 0;        //Set to FALSE to start
}

 AliEMCALJetFinderAlgoUA1Hybrid::~AliEMCALJetFinderAlgoUA1Hybrid()
   {
     //Destructor
     if (fDebug>0) Info("AliEMCALJetFinderAlgoUA1Revised","Beginning Destructor");
     delete[] fUnit;
   }

 void AliEMCALJetFinderAlgoUA1Hybrid::SetJetFindingParameters
                               (Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin, 
                               Float_t minMove, Float_t maxMove, Float_t bgMaxMove)
   {
     //Sets parameters for the JetFinding algorithm
     if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");

     SetNumUnits(numUnits);
     SetJetESeed(eSeed);
     SetConeRad(coneRad);
     SetJetEMin(jetEMin);
     SetEtMin(etMin);
     SetMinMove(minMove);
     SetMaxMove(maxMove);
     SetBGMaxMove(bgMaxMove);
   }

 void AliEMCALJetFinderAlgoUA1Hybrid::SetJetFindingParameters
                               (Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin)
   {
     //Sets fewer parameters for the JetFinding algorithm
     if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");

     SetNumUnits(numUnits);
     SetJetESeed(eSeed);
     SetConeRad(coneRad);
     SetJetEMin(jetEMin);
     SetEtMin(etMin);
     SetMinMove(fMinMove);
     SetMaxMove(fMaxMove);
     SetBGMaxMove(fBGMaxMove);
   }

 void AliEMCALJetFinderAlgoUA1Hybrid::InitUnitArray()
   {
     //Initialises unit array
     if(fArrayInitialised) delete[] fUnit;
     fUnit = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits];
     fArrayInitialised = 1;
   }

 void AliEMCALJetFinderAlgoUA1Hybrid::FillUnitArray(AliEMCALJetFinderAlgoUA1FillUnitFlagType_t flag)
   {
     if (fDebug>1) Info("FillUnitArray","Beginning FillUnitArray");
     //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");

         //   if (pEMCAL){ 
         //	     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
         //     }else
         //    {
     //AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
     AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry() ; 
	  
     //    }
         
     AliEMCALJetFinderAlgoUA1FillUnitFlagType_t option = flag;
     Int_t         numTracks, numDigits;
    
     //Loops over all elements in the AliEMCALJetFinderAlgoUA1Unit array and 
     //fills the objects with relevant values from the Data Input object
     if (fDebug>1) Info("FillUnitArray","Filling array with Unit objects");
     if (fDebug>1) Info("FillUnitArray","NTracks %i NDigits %i",fInputPointer->GetNTracks(),fInputPointer->GetNDigits());
	 numTracks = fInputPointer->GetNTracks();
	 numDigits = fInputPointer->GetNDigits();
	 TParticle         *myPart;
	 AliEMCALDigit     *myDigit;

	 //Fill units with Track info if appropriate
	 if(option==kFillTracksOnly || option ==kFillAll) 
	   {          
	    for(Int_t j=0; j<numTracks; j++)
	    {
	     myPart = fInputPointer->GetTrack(j);
	     Float_t eta = myPart->Eta();
	     Float_t  phi = myPart->Phi();
	     Int_t towerID = geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi);
	     Float_t  pT = myPart->Pt();
	     Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy(); 

   	     //Do Hadron Correction
              if(fHadCorr != 0)
	       {
		 Double_t   fullP = myPart->P();
		 Double_t   hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
		 unitEnergy -= hCEnergy*TMath::Sin(myPart->Theta());
		 fUnit[towerID-1].SetUnitEnergy(unitEnergy);
	       } //end Hadron Correction loop 
	     
	     //Do Pt cut on tracks
	     if(fPtCut != 0 && pT < fPtCut) continue;

	     fUnit[towerID-1].SetUnitEnergy(unitEnergy+pT);

             }//end tracks loop
	   }//end Tracks condition


	 //Fill units with Digit info if appropriate
	 if(option ==kFillDigitsOnly || option ==kFillAll)
	   {
            for(Int_t k=0; k<numDigits; k++)
	    {
	     myDigit = fInputPointer->GetDigit(k);
             if (fDebug>1) Info("FillUnitArray","getting digits %i %i numdigits",k,numDigits );
	     Int_t towerID = myDigit->GetId();
	     Int_t amplitude = myDigit->GetAmp();     //Gets the integer valued amplitude of the digit
	     Float_t amp = (Float_t)amplitude;        //Need to typecast to Float_t before doing real energy conversion
	     Float_t digitEnergy = amp/10000000.0;    //Factor of 10 million needed to convert!
	     Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy() + digitEnergy;
	     fUnit[towerID-1].SetUnitEnergy(unitEnergy);

	    }//end digits loop
	   }//end digits condition

	 //Set all unit flags, Eta, Phi
	 for(Int_t i=0; i<fNumUnits; i++)
	   {
             if (fDebug>1) Info("FillUnitArray","Setting all units outside jets");
	     fUnit[i].SetUnitFlag(kOutJet);           //Set all units to be outside a jet initially
	     fUnit[i].SetUnitID(i+1);
	     Float_t eta;
	     Float_t phi;
	     geom->EtaPhiFromIndex(fUnit[i].GetUnitID(), eta, phi);
	     fUnit[i].SetUnitEta(eta);
	     fUnit[i].SetUnitPhi(phi*TMath::Pi()/180.0);
	     //	     if(i>13000) cout<<"!!!!!!!!!!!!!!!!!For unit0, eta="<<eta<<" and phi="<<phi*TMath::Pi()/180.0<<" and ID="<<fUnit[i].GetUnitID()<<endl;
	     //  if(fUnit[i].GetUnitEnergy()>0) cout<<"Unit ID "<<fUnit[i].GetUnitID() <<"with eta="<<eta<<" and phi="<<phi*TMath::Pi()/180.0<<" has energy="<<fUnit[i].GetUnitEnergy()<<endl;
	   }//end loop over all units in array (same as all towers in EMCAL)
   }


 void AliEMCALJetFinderAlgoUA1Hybrid::Sort(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t length)
 {
   //Calls the recursive quicksort method to sort unit objects in decending order of Energy
   if (fDebug>1) Info("Sort","Sorting Unit objects");
   QS(unit, 0, length-1);
 }
  

 void AliEMCALJetFinderAlgoUA1Hybrid::QS(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t left, Int_t right)
 {
  //Sorts the AliEMCALJetFinderAlgoUA1Unit objects in decending order of Energy
   if (fDebug>111) Info("QS","QuickSorting Unit objects");   

   Int_t    i;
   Int_t    j;
   AliEMCALJetFinderAlgoUA1Unit  unitFirst;
   AliEMCALJetFinderAlgoUA1Unit  unitSecond;

   i = left;
   j = right;
   unitFirst = unit[(left+right)/2];

 do
  {
    while( (unit[i].GetUnitEnergy() > unitFirst.GetUnitEnergy()) && (i < right)) i++;
    while( (unitFirst.GetUnitEnergy() > unit[j].GetUnitEnergy()) && (j > left)) j--;

    if(i <= j)
      {
	unitSecond = unit[i];
	unit[i] = unit[j];
	unit[j] = unitSecond;
	i++;
	j--;
      }//end if
  }while(i <= j);

 if(left < j) QS(unit, left, j);
 if(i < right) QS(unit, i, right);
 }


 void AliEMCALJetFinderAlgoUA1Hybrid::FindBG()
   {
     //Finds the background energy for the iteration
     if (fDebug>1) Info("FindBG","Finding Average Background"); 

     fEBGTotal          = 0.0;
     Int_t numCone      = 0;

     //Loop over all unit objects in the array and sum the energy of those not in a jet
     for(Int_t i=0; i<fNumUnits; i++)
       {
	 if(fUnit[i].GetUnitFlag() != kInJet)
	   fEBGTotal += fUnit[i].GetUnitEnergy();
	 else numCone++;
       }//end for

     fEBGTotalOld = fEBGTotal;
     fEBGAve = fEBGTotal / (fNumUnits - numCone);
     if (fDebug>5) Info("FindBG","Average BG is %f: ",fEBGAve);      

     for(Int_t count=0; count<fNumUnits;count++)
       {
	 fUnit[count].SetUnitFlag(kOutJet);
       }//end for
   }


 void AliEMCALJetFinderAlgoUA1Hybrid::FindJetEtaPhi(Int_t counter)
   {
     //Finds the eta and phi of the jet axis
     if (fDebug>1) Info("FindJetEtaPhi","Finding Jet Eta and Phi");

     fDEta = fUnit[counter].GetUnitEta() - fEtaInit;
     fDPhi = fUnit[counter].GetUnitPhi() - fPhiInit;

     fEnergy = fUnit[counter].GetUnitEnergy() - fEBGAve;
     fJetEtaSum += fEnergy * fDEta;
     fJetPhiSum += fEnergy * fDPhi;
     fJetESum += fEnergy;
     fJetEta = fEtaInit + (fJetEtaSum / fJetESum);
     fJetPhi = fPhiInit + (fJetPhiSum / fJetESum);
   }


 void AliEMCALJetFinderAlgoUA1Hybrid::FindJetEnergy()
   {
     //Finds the energy of the jet after the final axis has been found
     if (fDebug>1) Info("FindJetEnergy","Finding Jet Energy");

     for(Int_t i=0; i<fNumUnits; i++)
       {
	 //Loop over all unit objects in the array and find if within cone radius
	 Float_t dEta = fUnit[i].GetUnitEta() - fJetEta;
	 Float_t dPhi = fUnit[i].GetUnitPhi() - fJetPhi;
	 Float_t rad = TMath::Sqrt( (dEta*dEta) + (dPhi*dPhi) );

	 if(fUnit[i].GetUnitFlag()==kOutJet && rad<= fConeRad)
	   {
	     fUnit[i].SetUnitFlag(kInCurrentJet);
	     Float_t energy = fUnit[i].GetUnitEnergy() - fEBGAve;
	     fJetESum += energy;                             
	     fJetEtaSum += energy * dEta;
	     fJetPhiSum += energy * dPhi;
	     fNumInCone++;                     //Increment the number of cells in the jet cone
	   }//end if
       }//end for
   }


 void AliEMCALJetFinderAlgoUA1Hybrid::StoreJetInfo()
   {
     //Stores the resulting jet information in appropriate storage structure (TO BE DECIDED!!!!)
     if (fDebug>1) Info("StoreJetInfo","Storing Jet Information");
    
     //Store:
     //fJetESum is the final jet energy (background has been subtracted)
     //fJetEta is the final jet Eta
     //fJetPhi is the final jet Phi
     //fNumInCone is the final number of cells included in the jet cone
     //fEtaInit is the eta of the initiator cell
     //fPhiInit is the phi of the initiator cell
     fJet.SetEnergy(fJetESum);
     fJet.SetEta(fJetEta);
     fJet.SetPhi(fJetPhi);

      cout<<"For iteration "<<fNumIter <<" and Jet number " <<fNumJets <<endl;
      cout<<"The jet energy is: " <<fJetESum <<endl;
      cout<<"The jet eta is ---->" <<fJetEta <<endl;
      cout<<"The jet phi is ---->" <<fJetPhi <<endl;

     Int_t             numberTracks = fInputPointer->GetNTracks();
     TParticle         *myP;
     Int_t             numTracksInCone = 0;

     for(Int_t counter=0; counter<numberTracks; counter++)
       {
	myP = fInputPointer->GetTrack(counter);
	Float_t eta = myP->Eta();
	Float_t phi = myP->Phi(); 
	Float_t deta = fJetEta-eta;
	Float_t dphi = fJetPhi -phi;
	Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi));
	if(rad<=fConeRad) numTracksInCone++;
       }//end for

     Float_t    *pTArray = new Float_t[numTracksInCone];
     Float_t    *etaArray = new Float_t[numTracksInCone];
     Float_t    *phiArray = new Float_t[numTracksInCone];
     Int_t      *pdgArray = new Int_t[numTracksInCone];
     Int_t             index = 0;

     for(Int_t counter2=0; counter2<numberTracks; counter2++)
       {
	 myP = fInputPointer->GetTrack(counter2);
 	 Float_t eta = myP->Eta();
         Float_t phi = myP->Phi(); 
	 Float_t deta = fJetEta-eta;
	 Float_t dphi = fJetPhi -phi;
	 Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi));
	 if(rad<=fConeRad)
	   {
	     pTArray[index] = myP->Pt();
	     etaArray[index] = eta;
	     phiArray[index] = phi;
	     pdgArray[index] = myP->GetPdgCode();
	     index++;
	   }//end if
       }//end for

     fJet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray);
     fOutputObject.AddJet(&fJet);
     delete[] pTArray;
     delete[] etaArray;
     delete[] phiArray;
     delete[] pdgArray;
   }


 void AliEMCALJetFinderAlgoUA1Hybrid::FindJets()
   {
     //Runs the complete UA1 JetFinding algorithm to find jets!
     if (fDebug>1) Info("FindJets","Starting Jet Finding!!!");

     //If the array of JetFinderUnit objects has not been initialised then initialise with default settings
     if(!fArrayInitialised) 
      {
       InitUnitArray();
       FillUnitArray(kFillAll);
      }//end if
     if (fDebug>1) Info("FindJets","Unit array filled");

     //Step 1. Sort the array in descending order of Energy
     Sort(fUnit,fNumUnits);

     //Step 2. Set the number of iterations and Number of jets found to zero to start
     fNumIter = 0;
     fNumJets = 0;

     //Step 3. Begin the iteration loop to find jets
     //Need to iterate the algorithm while number of iterations<2 OR number of iterations<10 AND 
     //the value of the average background has changed more than specified amount
     //Min iterations = 2, Max iterations = 10
     //while(fNumIter<2 || (fNumIter <10 && ( (fEBGTotal-fEBGTotalOld)/fEBGTotal) > fBGMaxMove) )

     while(fNumIter<2 || (fNumIter <10 && ( fEBGTotal-fEBGTotalOld) > fEBGTotal*fBGMaxMove) )
       {
        if (fDebug>1) Info("FindJets","Starting BIG iteration ---> %i",fNumIter);

         //Step 4. Find the value of the average background energy
	 FindBG();
	 fOutputObject.Reset(kResetJets); //Reset output object to store info for new iteration
	 fNumJets=0;

	 //Loop over the array of unit objects and flag those with energy below MinCellEt
         Int_t numbelow = 0;
	 for(Int_t j=0; j<fNumUnits; j++)
	   {
	     if( (fUnit[j].GetUnitEnergy()-fEBGAve) < fEtMin)
        	{       
		  //          fUnit[j].SetUnitFlag(kBelowMinEt);    TAKING OUT kBelow flag
                  numbelow++;
                }//end if
	   }//end for
	 //cout<<"THERE WERE "<<numbelow<<" units with E <EtMin!!!!!!!!!!!!!!!"<<endl;

	 //Do quick check if there are no jets upfront
	 // if(fUnit[0].GetUnitFlag() == kBelowMinEt)
	 if( (fUnit[0].GetUnitEnergy()-fEBGAve) < fEtMin)
	   {
            cout <<"There are no jets for this event!" <<endl;
	    break;
	   }//end if

         //Step 5. Begin with the first jet candidate cell (JET SEED LOOP)
         if (fDebug>5) Info("FindJets","Beginning JET SEED LOOP");
	 for(Int_t count=0; count<fNumUnits; count++)
	   {

//CHECK CONDITION HERE _ NOT SURE IF SHOULD MAYBE BE: GetUnitEnergy()-fEBGAve >fESeed?????????????????????????????
	     if(fUnit[count].GetUnitEnergy()>=fESeed && fUnit[count].GetUnitFlag()==kOutJet)
	       {
		 fEnergy = fUnit[count].GetUnitEnergy() - fEBGAve;
		 fJetEta = fUnit[count].GetUnitEta();
		 fJetPhi = fUnit[count].GetUnitPhi();
		 Int_t seedID = fUnit[count].GetUnitID();
                 if (fDebug>5) Info("FindJets","Inside first candidate jet seed loop for time : %i", count);
                 if (fDebug>5) Info("FindJets","Found candidate energy %f ",fEnergy);
                 if (fDebug>5) Info("FindJets","Found candidate eta %f ", fJetEta);
                 if (fDebug>5) Info("FindJets","Found candidate phi %f ", fJetPhi);
		 if (fDebug>5) Info("FindJets","Found candidate ID %i", seedID);

		 fEtaInit = fJetEta;
		 fPhiInit = fJetPhi;
		 fEtaB = fJetEta;
		 fPhiB = fJetPhi;
		 fJetESum = 0.0;
		 fJetEtaSum = 0.0;
		 fJetPhiSum = 0.0;
       
         //Step 6. Find Jet Eta and Phi
		 //Loop over all units in the array to find the ones in the jet cone and determine contrib to Jet eta, phi
		 do
		   {
		     for(Int_t count1=0; count1<fNumUnits; count1++)		   
		       {
			 if(fUnit[count1].GetUnitID() == seedID) continue;   //skip unit if the jetseed to avoid doublecounting
			 if(fUnit[count1].GetUnitFlag() == kOutJet)
			   {
			     fDEta = fUnit[count1].GetUnitEta() - fJetEta;
			     fDPhi = fUnit[count1].GetUnitPhi() - fJetPhi;
			     fRad = TMath::Sqrt( (fDEta*fDEta) + (fDPhi*fDPhi) );
			     if(fRad <= fConeRad)
			       {
				 FindJetEtaPhi(count1); 
			       }//end if
			   }//end if
		       }//end for (Jet Eta, Phi LOOP)
			     
		      //Find the distance cone centre moved from previous cone centre
		      if (fDebug>10) Info("FindJets","Checking if cone move small enough");
		      fDistP = TMath::Sqrt( ((fJetEta-fEtaB)*(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)*(fJetPhi-fPhiB)) );
		      //     if(fDistP <= fMinMove) break;
			     

		      //Find the distance cone centre is from initiator cell
		      if (fDebug>10) Info("FindJets","Checking if cone move too large");
		      fDistI = TMath::Sqrt( ((fJetEtaSum/fJetESum)*(fJetEtaSum/fJetESum)) + ((fJetPhiSum/fJetESum)*
												    (fJetPhiSum/fJetESum)));

		      if(fDistP>fMinMove && fDistI<fMaxMove)
			{
			  fEtaB = fJetEta;
			  fPhiB = fJetPhi;
			}//end if
		 
		   }while(fDistP>fMinMove && fDistI<fMaxMove);
			  
		 fJetEta = fEtaB;
		 fJetPhi = fPhiB;


       //Step 7. Find the Jet Energy
                 if (fDebug>1) Info("FindJets","Looking for Jet energy");
		 fJetESum = 0.0;
		 fJetEtaSum = 0.0;
		 fJetPhiSum = 0.0;
		 fNumInCone = 0;
		 FindJetEnergy();

		 //cout<<"Number of cells in jet cone is: "<<fNumInCone<<endl;

       //Step 8. Check if the jet is a valid jet
		 //Check if cluster energy is above Min allowed to be a jet
//DID NOT DO THE COSH COMPARISON HERE -> NEED TO CHECK WHICH COMPARISON IS BEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                 if (fDebug>5) Info("FindJets","Checking cluster is valid jet");
		 if(fJetESum < fJetEMin)
		   {
		     for(Int_t count2=0; count2<fNumUnits; count2++)
		       {
			 if(fUnit[count2].GetUnitFlag()==kInCurrentJet || fUnit[count2].GetUnitFlag()==kOutJet)
			   fUnit[count2].SetUnitFlag(kOutJet);
		       }//end for
                   if (fDebug>10) Info("FindJets","NOT a valid jet cell");
		  }else
		    {
		     for(Int_t count2=0; count2<fNumUnits; count2++)
		       {
			 if(fUnit[count2].GetUnitFlag()==kInCurrentJet)
			   {
			     //	     cout<<"Setting unit #"<<count2 <<" to be officially in a jet!"<<endl;
                           fUnit[count2].SetUnitFlag(kInJet);
			   }
		       }//end for			

 //NEED TO CHECK FINAL WEIRD ITERATION OF ETA AND PHI CHANGES!!!!!!!!!
		     //	 fJetPhi += fJetPhiSum/fJetESum;        //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		     //  fJetEta += fJetEtaSum/fJetESum;        //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		     fNumJets++;              //Incrementing number of jets found
		     StoreJetInfo();          //Storing jet info

		 }//end if (check cluster above Min Jet Energy)
	       }//end if (Jet Seed condition)
	   }//end (JET SEED LOOP)

if (fDebug>5) Info("FindJets","End of BIG iteration number %i",fNumIter);
// this->Dump();
	 fNumIter++;
       }//end 10 iteration WHILE LOOP
 }
Commit	Line	Data
f7d5860b	1	//SARAH'S HYBRID PERSONAL COPY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2
	3	//THIS IS SARAH'S REVISED UA1 CODE WITH CHANGES FOR ETA/PHI ITERATION INCLUDED!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4	//THIS Also includes summing ALL cells in the jetcone towards the jet energy NOT just those above threshold!!!!!
	5
	6
	7	/**************************************************************************
	8	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	9	* *
	10	* Author: The ALICE Off-line Project. *
	11	* Contributors are mentioned in the code where appropriate. *
	12	* *
	13	* Permission to use, copy, modify and distribute this software and its *
	14	* documentation strictly for non-commercial purposes is hereby granted *
	15	* without fee, provided that the above copyright notice appears in all *
	16	* copies and that both the copyright notice and this permission notice *
	17	* appear in the supporting documentation. The authors make no claims *
	18	* about the suitability of this software for any purpose. It is *
	19	* provided "as is" without express or implied warranty. *
	20	**************************************************************************/
	21
	22	/*
	23
	24	$Log$
	25	Revision 1.3 2003/06/18 17:00:25 horner
	26	Cleaned up more cout statements - Mark
	27
	28	Revision 1.2 2003/06/18 00:20:07 horner
	29	Removed cout statement - Mark
	30
	31
	32	*/
	33
	34
	35	//*--Author: Sarah Blyth (LBL)
	36	//*--Based on UA1 jet algorithm from LUND JETSET called from EMC-erj
	37
	38	#include "TTask.h"
	39	#include "AliEMCALJetFinderInput.h"
	40	#include "AliEMCALJetFinderOutput.h"
	41	#include "AliEMCALJetFinderAlgo.h"
	42	#include "AliEMCALJetFinderAlgoUA1Hybrid.h"
	43	#include "AliEMCALJetFinderAlgoUA1Unit.h"
	44	#include "AliEMCALGeometry.h"
	45	#include "AliEMCAL.h"
	46	#include "AliEMCALGetter.h"
	47	#include "AliEMCALDigit.h"
	48	#include "TParticle.h"
	49	#include "AliRun.h"
	50	#include "AliEMCALJet.h"
	51	#include "TMath.h"
	52
	53
	54	ClassImp(AliEMCALJetFinderAlgoUA1Hybrid)
	55
	56	AliEMCALJetFinderAlgoUA1Hybrid::AliEMCALJetFinderAlgoUA1Hybrid()
	57	{
	58	//Default constructor
	59	if (fDebug>0) Info("AliEMCALJetFinderAlgoUA1Hybrid","Beginning Default Constructor");
	60
	61	fNumIter = 0;
	62	fNumUnits = 13824; //Number of towers in EMCAL
	63	fESeed = 5.0; //Default value
	64	fConeRad = 0.5; //Default value
65	fJetEMin = 10.0; //Default value
66	fEtMin = 0.28; //Default value
67	fMinMove = 0.05; //From original UA1 JetFinder
68	fMaxMove = 0.15; //From original UA1 JetFinder
69	fBGMaxMove = 0.035; //From original UA1 JetFinder
70	fPtCut = 0;
71	fHadCorr = 0;
72	fEBGTotal = 1.0; //Set to 1 so that no div by zero in first FindJets() loop
73	fEBGTotalOld = 0.0;
74	fEBGAve = 0.0;
75	fEnergy = 0.0;
76	fJetEta = 0.0;
77	fJetPhi = 0.0;
78	fEtaInit = 0.0;
79	fPhiInit = 0.0;
80	fEtaB = 0.0;
81	fPhiB = 0.0;
82	fJetESum = 0.0;
83	fJetEtaSum = 0.0;
84	fJetPhiSum = 0.0;
85	fDEta = 0.0;
86	fDPhi = 0.0;
87	fDistP = 0.0;
88	fDistI = 0.0;
89	fTempE = 0.0;
90	fRad = 2.0; //Set to 2 to start
91	fNumInCone = 0;
92	fNumJets = 0;
93	fArrayInitialised = 0; //Set to FALSE to start
94	}
95
96	AliEMCALJetFinderAlgoUA1Hybrid::~AliEMCALJetFinderAlgoUA1Hybrid()
97	{
98	//Destructor
99	if (fDebug>0) Info("AliEMCALJetFinderAlgoUA1Revised","Beginning Destructor");
100	delete[] fUnit;
101	}
102
103	void AliEMCALJetFinderAlgoUA1Hybrid::SetJetFindingParameters
104	(Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin,
105	Float_t minMove, Float_t maxMove, Float_t bgMaxMove)
106	{
107	//Sets parameters for the JetFinding algorithm
108	if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");
109
110	SetNumUnits(numUnits);
111	SetJetESeed(eSeed);
112	SetConeRad(coneRad);
113	SetJetEMin(jetEMin);
114	SetEtMin(etMin);
115	SetMinMove(minMove);
116	SetMaxMove(maxMove);
117	SetBGMaxMove(bgMaxMove);
118	}
119
120	void AliEMCALJetFinderAlgoUA1Hybrid::SetJetFindingParameters
121	(Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin)
122	{
123	//Sets fewer parameters for the JetFinding algorithm
124	if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");
125
126	SetNumUnits(numUnits);
127	SetJetESeed(eSeed);
128	SetConeRad(coneRad);
129	SetJetEMin(jetEMin);
130	SetEtMin(etMin);
131	SetMinMove(fMinMove);
132	SetMaxMove(fMaxMove);
133	SetBGMaxMove(fBGMaxMove);
134	}
135
136	void AliEMCALJetFinderAlgoUA1Hybrid::InitUnitArray()
137	{
138	//Initialises unit array
139	if(fArrayInitialised) delete[] fUnit;
140	fUnit = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits];
141	fArrayInitialised = 1;
142	}
143
144	void AliEMCALJetFinderAlgoUA1Hybrid::FillUnitArray(AliEMCALJetFinderAlgoUA1FillUnitFlagType_t flag)
145	{
146	if (fDebug>1) Info("FillUnitArray","Beginning FillUnitArray");
147	//AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
148
149	// if (pEMCAL){
150	// AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
151	// }else
152	// {
153	//AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
154	AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
155
156	// }
157
158	AliEMCALJetFinderAlgoUA1FillUnitFlagType_t option = flag;
159	Int_t numTracks, numDigits;
160
161	//Loops over all elements in the AliEMCALJetFinderAlgoUA1Unit array and
162	//fills the objects with relevant values from the Data Input object
163	if (fDebug>1) Info("FillUnitArray","Filling array with Unit objects");
164	if (fDebug>1) Info("FillUnitArray","NTracks %i NDigits %i",fInputPointer->GetNTracks(),fInputPointer->GetNDigits());
165	numTracks = fInputPointer->GetNTracks();
166	numDigits = fInputPointer->GetNDigits();
167	TParticle *myPart;
168	AliEMCALDigit *myDigit;
169
170	//Fill units with Track info if appropriate
171	if(option==kFillTracksOnly \|\| option ==kFillAll)
172	{
173	for(Int_t j=0; j<numTracks; j++)
174	{
175	myPart = fInputPointer->GetTrack(j);
176	Float_t eta = myPart->Eta();
177	Float_t phi = myPart->Phi();
178	Int_t towerID = geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi);
179	Float_t pT = myPart->Pt();
180	Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy();
181
182	//Do Hadron Correction
183	if(fHadCorr != 0)
184	{
185	Double_t fullP = myPart->P();
186	Double_t hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
187	unitEnergy -= hCEnergy*TMath::Sin(myPart->Theta());
188	fUnit[towerID-1].SetUnitEnergy(unitEnergy);
189	} //end Hadron Correction loop
190
191	//Do Pt cut on tracks
192	if(fPtCut != 0 && pT < fPtCut) continue;
193
194	fUnit[towerID-1].SetUnitEnergy(unitEnergy+pT);
195
196	}//end tracks loop
197	}//end Tracks condition
198
199
200	//Fill units with Digit info if appropriate
201	if(option ==kFillDigitsOnly \|\| option ==kFillAll)
202	{
203	for(Int_t k=0; k<numDigits; k++)
204	{
205	myDigit = fInputPointer->GetDigit(k);
206	if (fDebug>1) Info("FillUnitArray","getting digits %i %i numdigits",k,numDigits );
207	Int_t towerID = myDigit->GetId();
208	Int_t amplitude = myDigit->GetAmp(); //Gets the integer valued amplitude of the digit
209	Float_t amp = (Float_t)amplitude; //Need to typecast to Float_t before doing real energy conversion
210	Float_t digitEnergy = amp/10000000.0; //Factor of 10 million needed to convert!
211	Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy() + digitEnergy;
212	fUnit[towerID-1].SetUnitEnergy(unitEnergy);
213
214	}//end digits loop
215	}//end digits condition
216
217	//Set all unit flags, Eta, Phi
218	for(Int_t i=0; i<fNumUnits; i++)
219	{
220	if (fDebug>1) Info("FillUnitArray","Setting all units outside jets");
221	fUnit[i].SetUnitFlag(kOutJet); //Set all units to be outside a jet initially
222	fUnit[i].SetUnitID(i+1);
223	Float_t eta;
224	Float_t phi;
225	geom->EtaPhiFromIndex(fUnit[i].GetUnitID(), eta, phi);
226	fUnit[i].SetUnitEta(eta);
227	fUnit[i].SetUnitPhi(phi*TMath::Pi()/180.0);
228	// if(i>13000) cout<<"!!!!!!!!!!!!!!!!!For unit0, eta="<<eta<<" and phi="<<phi*TMath::Pi()/180.0<<" and ID="<<fUnit[i].GetUnitID()<<endl;
229	// if(fUnit[i].GetUnitEnergy()>0) cout<<"Unit ID "<<fUnit[i].GetUnitID() <<"with eta="<<eta<<" and phi="<<phi*TMath::Pi()/180.0<<" has energy="<<fUnit[i].GetUnitEnergy()<<endl;
230	}//end loop over all units in array (same as all towers in EMCAL)
231	}
232
233
234	void AliEMCALJetFinderAlgoUA1Hybrid::Sort(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t length)
235	{
236	//Calls the recursive quicksort method to sort unit objects in decending order of Energy
237	if (fDebug>1) Info("Sort","Sorting Unit objects");
238	QS(unit, 0, length-1);
239	}
240
241
242	void AliEMCALJetFinderAlgoUA1Hybrid::QS(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t left, Int_t right)
243	{
244	//Sorts the AliEMCALJetFinderAlgoUA1Unit objects in decending order of Energy
245	if (fDebug>111) Info("QS","QuickSorting Unit objects");
246
247	Int_t i;
248	Int_t j;
249	AliEMCALJetFinderAlgoUA1Unit unitFirst;
250	AliEMCALJetFinderAlgoUA1Unit unitSecond;
251
252	i = left;
253	j = right;
254	unitFirst = unit[(left+right)/2];
255
256	do
257	{
258	while( (unit[i].GetUnitEnergy() > unitFirst.GetUnitEnergy()) && (i < right)) i++;
259	while( (unitFirst.GetUnitEnergy() > unit[j].GetUnitEnergy()) && (j > left)) j--;
260
261	if(i <= j)
262	{
263	unitSecond = unit[i];
264	unit[i] = unit[j];
265	unit[j] = unitSecond;
266	i++;
267	j--;
268	}//end if
269	}while(i <= j);
270
271	if(left < j) QS(unit, left, j);
272	if(i < right) QS(unit, i, right);
273	}
274
275
276	void AliEMCALJetFinderAlgoUA1Hybrid::FindBG()
277	{
278	//Finds the background energy for the iteration
279	if (fDebug>1) Info("FindBG","Finding Average Background");
280
281	fEBGTotal = 0.0;
282	Int_t numCone = 0;
283
284	//Loop over all unit objects in the array and sum the energy of those not in a jet
285	for(Int_t i=0; i<fNumUnits; i++)
286	{
287	if(fUnit[i].GetUnitFlag() != kInJet)
288	fEBGTotal += fUnit[i].GetUnitEnergy();
289	else numCone++;
290	}//end for
291
292	fEBGTotalOld = fEBGTotal;
293	fEBGAve = fEBGTotal / (fNumUnits - numCone);
294	if (fDebug>5) Info("FindBG","Average BG is %f: ",fEBGAve);
295
296	for(Int_t count=0; count<fNumUnits;count++)
297	{
298	fUnit[count].SetUnitFlag(kOutJet);
299	}//end for
300	}
301
302
303	void AliEMCALJetFinderAlgoUA1Hybrid::FindJetEtaPhi(Int_t counter)
304	{
305	//Finds the eta and phi of the jet axis
306	if (fDebug>1) Info("FindJetEtaPhi","Finding Jet Eta and Phi");
307
308	fDEta = fUnit[counter].GetUnitEta() - fEtaInit;
309	fDPhi = fUnit[counter].GetUnitPhi() - fPhiInit;
310
311	fEnergy = fUnit[counter].GetUnitEnergy() - fEBGAve;
312	fJetEtaSum += fEnergy * fDEta;
313	fJetPhiSum += fEnergy * fDPhi;
314	fJetESum += fEnergy;
315	fJetEta = fEtaInit + (fJetEtaSum / fJetESum);
316	fJetPhi = fPhiInit + (fJetPhiSum / fJetESum);
317	}
318
319
320	void AliEMCALJetFinderAlgoUA1Hybrid::FindJetEnergy()
321	{
322	//Finds the energy of the jet after the final axis has been found
323	if (fDebug>1) Info("FindJetEnergy","Finding Jet Energy");
324
325	for(Int_t i=0; i<fNumUnits; i++)
326	{
327	//Loop over all unit objects in the array and find if within cone radius
328	Float_t dEta = fUnit[i].GetUnitEta() - fJetEta;
329	Float_t dPhi = fUnit[i].GetUnitPhi() - fJetPhi;
330	Float_t rad = TMath::Sqrt( (dEtadEta) + (dPhidPhi) );
331
332	if(fUnit[i].GetUnitFlag()==kOutJet && rad<= fConeRad)
333	{
334	fUnit[i].SetUnitFlag(kInCurrentJet);
335	Float_t energy = fUnit[i].GetUnitEnergy() - fEBGAve;
336	fJetESum += energy;
337	fJetEtaSum += energy * dEta;
338	fJetPhiSum += energy * dPhi;
339	fNumInCone++; //Increment the number of cells in the jet cone
340	}//end if
341	}//end for
342	}
343
344
345	void AliEMCALJetFinderAlgoUA1Hybrid::StoreJetInfo()
346	{
347	//Stores the resulting jet information in appropriate storage structure (TO BE DECIDED!!!!)
348	if (fDebug>1) Info("StoreJetInfo","Storing Jet Information");
349
350	//Store:
351	//fJetESum is the final jet energy (background has been subtracted)
352	//fJetEta is the final jet Eta
353	//fJetPhi is the final jet Phi
354	//fNumInCone is the final number of cells included in the jet cone
355	//fEtaInit is the eta of the initiator cell
356	//fPhiInit is the phi of the initiator cell
357	fJet.SetEnergy(fJetESum);
358	fJet.SetEta(fJetEta);
359	fJet.SetPhi(fJetPhi);
360
361	cout<<"For iteration "<<fNumIter <<" and Jet number " <<fNumJets <<endl;
362	cout<<"The jet energy is: " <<fJetESum <<endl;
363	cout<<"The jet eta is ---->" <<fJetEta <<endl;
364	cout<<"The jet phi is ---->" <<fJetPhi <<endl;
365
366	Int_t numberTracks = fInputPointer->GetNTracks();
367	TParticle *myP;
368	Int_t numTracksInCone = 0;
369
370	for(Int_t counter=0; counter<numberTracks; counter++)
371	{
372	myP = fInputPointer->GetTrack(counter);
373	Float_t eta = myP->Eta();
374	Float_t phi = myP->Phi();
375	Float_t deta = fJetEta-eta;
376	Float_t dphi = fJetPhi -phi;
377	Float_t rad = TMath::Sqrt( (detadeta) + (dphidphi));
378	if(rad<=fConeRad) numTracksInCone++;
379	}//end for
380
381	Float_t *pTArray = new Float_t[numTracksInCone];
382	Float_t *etaArray = new Float_t[numTracksInCone];
383	Float_t *phiArray = new Float_t[numTracksInCone];
384	Int_t *pdgArray = new Int_t[numTracksInCone];
385	Int_t index = 0;
386
387	for(Int_t counter2=0; counter2<numberTracks; counter2++)
388	{
389	myP = fInputPointer->GetTrack(counter2);
390	Float_t eta = myP->Eta();
391	Float_t phi = myP->Phi();
392	Float_t deta = fJetEta-eta;
393	Float_t dphi = fJetPhi -phi;
394	Float_t rad = TMath::Sqrt( (detadeta) + (dphidphi));
395	if(rad<=fConeRad)
396	{
397	pTArray[index] = myP->Pt();
398	etaArray[index] = eta;
399	phiArray[index] = phi;
400	pdgArray[index] = myP->GetPdgCode();
401	index++;
402	}//end if
403	}//end for
404
405	fJet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray);
406	fOutputObject.AddJet(&fJet);
407	delete[] pTArray;
408	delete[] etaArray;
409	delete[] phiArray;
410	delete[] pdgArray;
411	}
412
413
414	void AliEMCALJetFinderAlgoUA1Hybrid::FindJets()
415	{
416	//Runs the complete UA1 JetFinding algorithm to find jets!
417	if (fDebug>1) Info("FindJets","Starting Jet Finding!!!");
418
419	//If the array of JetFinderUnit objects has not been initialised then initialise with default settings
420	if(!fArrayInitialised)
421	{
422	InitUnitArray();
423	FillUnitArray(kFillAll);
424	}//end if
425	if (fDebug>1) Info("FindJets","Unit array filled");
426
427	//Step 1. Sort the array in descending order of Energy
428	Sort(fUnit,fNumUnits);
429
430	//Step 2. Set the number of iterations and Number of jets found to zero to start
431	fNumIter = 0;
432	fNumJets = 0;
433
434	//Step 3. Begin the iteration loop to find jets
435	//Need to iterate the algorithm while number of iterations<2 OR number of iterations<10 AND
436	//the value of the average background has changed more than specified amount
437	//Min iterations = 2, Max iterations = 10
438	//while(fNumIter<2 \|\| (fNumIter <10 && ( (fEBGTotal-fEBGTotalOld)/fEBGTotal) > fBGMaxMove) )
439
440	while(fNumIter<2 \|\| (fNumIter <10 && ( fEBGTotal-fEBGTotalOld) > fEBGTotal*fBGMaxMove) )
441	{
442	if (fDebug>1) Info("FindJets","Starting BIG iteration ---> %i",fNumIter);
443
444	//Step 4. Find the value of the average background energy
445	FindBG();
446	fOutputObject.Reset(kResetJets); //Reset output object to store info for new iteration
447	fNumJets=0;
448
449	//Loop over the array of unit objects and flag those with energy below MinCellEt
450	Int_t numbelow = 0;
451	for(Int_t j=0; j<fNumUnits; j++)
452	{
453	if( (fUnit[j].GetUnitEnergy()-fEBGAve) < fEtMin)
454	{
455	// fUnit[j].SetUnitFlag(kBelowMinEt); TAKING OUT kBelow flag
456	numbelow++;
457	}//end if
458	}//end for
459	//cout<<"THERE WERE "<<numbelow<<" units with E <EtMin!!!!!!!!!!!!!!!"<<endl;
460
461	//Do quick check if there are no jets upfront
462	// if(fUnit[0].GetUnitFlag() == kBelowMinEt)
463	if( (fUnit[0].GetUnitEnergy()-fEBGAve) < fEtMin)
464	{
465	cout <<"There are no jets for this event!" <<endl;
466	break;
467	}//end if
468
469	//Step 5. Begin with the first jet candidate cell (JET SEED LOOP)
470	if (fDebug>5) Info("FindJets","Beginning JET SEED LOOP");
471	for(Int_t count=0; count<fNumUnits; count++)
472	{
473
474	//CHECK CONDITION HERE _ NOT SURE IF SHOULD MAYBE BE: GetUnitEnergy()-fEBGAve >fESeed?????????????????????????????
475	if(fUnit[count].GetUnitEnergy()>=fESeed && fUnit[count].GetUnitFlag()==kOutJet)
476	{
477	fEnergy = fUnit[count].GetUnitEnergy() - fEBGAve;
478	fJetEta = fUnit[count].GetUnitEta();
479	fJetPhi = fUnit[count].GetUnitPhi();
480	Int_t seedID = fUnit[count].GetUnitID();
481	if (fDebug>5) Info("FindJets","Inside first candidate jet seed loop for time : %i", count);
482	if (fDebug>5) Info("FindJets","Found candidate energy %f ",fEnergy);
483	if (fDebug>5) Info("FindJets","Found candidate eta %f ", fJetEta);
484	if (fDebug>5) Info("FindJets","Found candidate phi %f ", fJetPhi);
485	if (fDebug>5) Info("FindJets","Found candidate ID %i", seedID);
486
487	fEtaInit = fJetEta;
488	fPhiInit = fJetPhi;
489	fEtaB = fJetEta;
490	fPhiB = fJetPhi;
491	fJetESum = 0.0;
492	fJetEtaSum = 0.0;
493	fJetPhiSum = 0.0;
494
495	//Step 6. Find Jet Eta and Phi
496	//Loop over all units in the array to find the ones in the jet cone and determine contrib to Jet eta, phi
497	do
498	{
499	for(Int_t count1=0; count1<fNumUnits; count1++)
500	{
501	if(fUnit[count1].GetUnitID() == seedID) continue; //skip unit if the jetseed to avoid doublecounting
502	if(fUnit[count1].GetUnitFlag() == kOutJet)
503	{
504	fDEta = fUnit[count1].GetUnitEta() - fJetEta;
505	fDPhi = fUnit[count1].GetUnitPhi() - fJetPhi;
506	fRad = TMath::Sqrt( (fDEtafDEta) + (fDPhifDPhi) );
507	if(fRad <= fConeRad)
508	{
509	FindJetEtaPhi(count1);
510	}//end if
511	}//end if
512	}//end for (Jet Eta, Phi LOOP)
513
514	//Find the distance cone centre moved from previous cone centre
515	if (fDebug>10) Info("FindJets","Checking if cone move small enough");
516	fDistP = TMath::Sqrt( ((fJetEta-fEtaB)(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)(fJetPhi-fPhiB)) );
517	// if(fDistP <= fMinMove) break;
518
519
520	//Find the distance cone centre is from initiator cell
521	if (fDebug>10) Info("FindJets","Checking if cone move too large");
522	fDistI = TMath::Sqrt( ((fJetEtaSum/fJetESum)(fJetEtaSum/fJetESum)) + ((fJetPhiSum/fJetESum)
523	(fJetPhiSum/fJetESum)));
524
525	if(fDistP>fMinMove && fDistI<fMaxMove)
526	{
527	fEtaB = fJetEta;
528	fPhiB = fJetPhi;
529	}//end if
530
531	}while(fDistP>fMinMove && fDistI<fMaxMove);
532
533	fJetEta = fEtaB;
534	fJetPhi = fPhiB;
535
536
537	//Step 7. Find the Jet Energy
538	if (fDebug>1) Info("FindJets","Looking for Jet energy");
539	fJetESum = 0.0;
540	fJetEtaSum = 0.0;
541	fJetPhiSum = 0.0;
542	fNumInCone = 0;
543	FindJetEnergy();
544
545	//cout<<"Number of cells in jet cone is: "<<fNumInCone<<endl;
546
547	//Step 8. Check if the jet is a valid jet
548	//Check if cluster energy is above Min allowed to be a jet
549	//DID NOT DO THE COSH COMPARISON HERE -> NEED TO CHECK WHICH COMPARISON IS BEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
550	if (fDebug>5) Info("FindJets","Checking cluster is valid jet");
551	if(fJetESum < fJetEMin)
552	{
553	for(Int_t count2=0; count2<fNumUnits; count2++)
554	{
555	if(fUnit[count2].GetUnitFlag()==kInCurrentJet \|\| fUnit[count2].GetUnitFlag()==kOutJet)
556	fUnit[count2].SetUnitFlag(kOutJet);
557	}//end for
558	if (fDebug>10) Info("FindJets","NOT a valid jet cell");
559	}else
560	{
561	for(Int_t count2=0; count2<fNumUnits; count2++)
562	{
563	if(fUnit[count2].GetUnitFlag()==kInCurrentJet)
564	{
565	// cout<<"Setting unit #"<<count2 <<" to be officially in a jet!"<<endl;
566	fUnit[count2].SetUnitFlag(kInJet);
567	}
568	}//end for
569
570	//NEED TO CHECK FINAL WEIRD ITERATION OF ETA AND PHI CHANGES!!!!!!!!!
571	// fJetPhi += fJetPhiSum/fJetESum; //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
572	// fJetEta += fJetEtaSum/fJetESum; //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
573
574	fNumJets++; //Incrementing number of jets found
575	StoreJetInfo(); //Storing jet info
576
577	}//end if (check cluster above Min Jet Energy)
578	}//end if (Jet Seed condition)
579	}//end (JET SEED LOOP)
580
581	if (fDebug>5) Info("FindJets","End of BIG iteration number %i",fNumIter);
582	// this->Dump();
583	fNumIter++;
584	}//end 10 iteration WHILE LOOP
585	}
586
587
588
589
590
591
592
593
594
595
596