EMCAL/AliEMCALJetFinderAlgoUA1.cxx

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