EMCAL/AliEMCALJetFinderAlgoUA1Revised.cxx

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