Bug fixes. Log replaced by Id
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderAlgoOmni.cxx
1
2 //THIS Also includes summing ALL cells in the jetcone towards the jet energy NOT just those above threshold!!!!!
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 "AliEMCALJetFinderAlgoOmni.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(AliEMCALJetFinderAlgoOmni)
42
43   AliEMCALJetFinderAlgoOmni::AliEMCALJetFinderAlgoOmni()
44 {
45   //Default constructor
46 if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
47
48   fNumIter           = 0;
49   fNumUnits          = 13824;     //Number of towers in EMCAL
50   fESeed             = 5.0;       //Default value
51   fConeRad           = 0.3;       //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   fBGType            = kRatio;   //Set Ratio method as default BG subtraction method 
82   fBGPar             = -1.0;      //Set to 1 to start
83 }
84
85  AliEMCALJetFinderAlgoOmni::~AliEMCALJetFinderAlgoOmni()
86    {
87      //Destructor
88      if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Destructor");
89      delete[] fUnit;
90      delete[] fUnitNoCuts;
91    }
92
93  void AliEMCALJetFinderAlgoOmni::SetJetFindingParameters
94                                (Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin, 
95                                Float_t minMove, Float_t maxMove, Float_t bgMaxMove)
96    {
97      //Sets parameters for the JetFinding algorithm
98      if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");
99
100      SetNumUnits(numUnits);
101      SetJetESeed(eSeed);
102      SetConeRad(coneRad);
103      SetJetEMin(jetEMin);
104      SetEtMin(etMin);
105      SetMinMove(minMove);
106      SetMaxMove(maxMove);
107      SetBGMaxMove(bgMaxMove);
108    }
109
110  void AliEMCALJetFinderAlgoOmni::SetJetFindingParameters
111                                (Int_t numUnits, Float_t eSeed, Float_t coneRad, Float_t jetEMin, Float_t etMin)
112    {
113      //Sets fewer parameters for the JetFinding algorithm
114      if (fDebug>1) Info("SetJetFindingParameters","Setting parameters for JetFinding");
115
116      SetNumUnits(numUnits);
117      SetJetESeed(eSeed);
118      SetConeRad(coneRad);
119      SetJetEMin(jetEMin);
120      SetEtMin(etMin);
121      SetMinMove(fMinMove);
122      SetMaxMove(fMaxMove);
123      SetBGMaxMove(fBGMaxMove);
124    }
125
126  void AliEMCALJetFinderAlgoOmni::InitUnitArray()
127    {
128      //Initialises unit arrays
129      if(fArrayInitialised) delete[] fUnit;
130      fUnit = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits];
131      fUnitNoCuts = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits];
132      fArrayInitialised = 1;
133    }
134
135  void AliEMCALJetFinderAlgoOmni::FillUnitArray(AliEMCALJetFinderAlgoUA1FillUnitFlagType_t flag)
136    {
137      if (fDebug>1) Info("FillUnitArray","Beginning FillUnitArray");
138      //     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
139
140          //   if (pEMCAL){ 
141          //          AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
142          //     }else
143          //    {
144      //AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
145      AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry() ; 
146        //    }
147          
148      AliEMCALJetFinderAlgoUA1FillUnitFlagType_t option = flag;
149      Int_t         numTracks, numDigits;
150     
151      //Loops over all elements in the AliEMCALJetFinderAlgoUA1Unit array and 
152      //fills the objects with relevant values from the Data Input object
153      if (fDebug>10) Info("FillUnitArray","Filling array with Unit objects");
154      if (fDebug>15) Info("FillUnitArray","NTracks %i NDigits %i",fInputPointer->GetNTracks(),fInputPointer->GetNDigits());
155          numTracks = fInputPointer->GetNTracks();
156          numDigits = fInputPointer->GetNDigits();
157          TParticle         *myPart;
158          AliEMCALDigit     *myDigit;
159
160          //Fill units with Track info if appropriate
161          if(option==kFillTracksOnly || option ==kFillAll) 
162            {          
163             for(Int_t j=0; j<numTracks; j++)
164             {
165              myPart = fInputPointer->GetTrack(j);
166              Float_t eta = myPart->Eta();
167              Float_t  phi = myPart->Phi();
168              Int_t towerID = geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi);
169              Float_t  pT = myPart->Pt();
170              Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy(); 
171              Float_t unitEnergyNoCuts = fUnitNoCuts[towerID-1].GetUnitEnergy();
172
173              /*
174              //OLD WAY:              //Do Hadron Correction
175               if(fHadCorr != 0)
176                {
177                  Double_t   fullP = myPart->P();
178                  Double_t   hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
179                  unitEnergy -= hCEnergy*TMath::Sin(myPart->Theta());
180                  fUnit[towerID-1].SetUnitEnergy(unitEnergy);
181                } //end Hadron Correction loop 
182              */      
183             
184               //Do Hadron Correction with propagate phi for the track
185               if(fHadCorr != 0)
186                 {
187                   Bool_t curl = 1;
188                   Float_t deltaPhi;
189                   TParticlePDG *pdg = myPart->GetPDG();
190                   if(pdg->Charge() < 0)
191                     {
192                       deltaPhi = PropagatePhi(myPart->Pt(), -1.0, curl);
193                     }
194                   else{
195                     deltaPhi = PropagatePhi(myPart->Pt(), 1.0, curl);
196                   }
197                   phi += deltaPhi;
198                   //Get new tower id for cell that track would curve into
199                   Int_t towerID2;
200                   if(phi<(1.0/3.0)*TMath::Pi() || phi>TMath::Pi())
201                     {
202                       towerID2 = -1;
203                     }
204                   else{
205                       towerID2 = geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi);
206                     }
207                   
208                   if(towerID2 != -1)
209                     {
210                       //Find unit energy of new tower
211                       Float_t unitEnergy2 = fUnit[towerID2-1].GetUnitEnergy();
212                       Float_t unitEnergy2NoCuts = fUnitNoCuts[towerID2-1].GetUnitEnergy();
213                       Double_t   fullP = myPart->P();
214                       Double_t   hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
215                       unitEnergy2 -= hCEnergy*TMath::Sin(myPart->Theta());
216                       unitEnergy2NoCuts -= hCEnergy*TMath::Sin(myPart->Theta());
217                       fUnit[towerID2-1].SetUnitEnergy(unitEnergy2);
218                       fUnitNoCuts[towerID2-1].SetUnitEnergy(unitEnergy2NoCuts);
219                     }//end if for towerID2
220                 }//end Hadron Correction loop
221
222               fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts + pT);
223              //Do Pt cut on tracks
224              if(fPtCut != 0 && pT < fPtCut) continue;
225
226              fUnit[towerID-1].SetUnitEnergy(unitEnergy+pT);
227
228              }//end tracks loop
229            }//end Tracks condition
230
231
232          //Fill units with Digit info if appropriate
233          if(option ==kFillDigitsOnly || option ==kFillAll)
234            {
235             for(Int_t k=0; k<numDigits; k++)
236             {
237              myDigit = fInputPointer->GetDigit(k);
238              if (fDebug>10) Info("FillUnitArray","getting digits %i %i numdigits",k,numDigits );
239              Int_t towerID = myDigit->GetId();
240              Int_t amplitude = myDigit->GetAmp();     //Gets the integer valued amplitude of the digit
241              Float_t amp = (Float_t)amplitude;        //Need to typecast to Float_t before doing real energy conversion
242              Float_t digitEnergy = amp/10000000.0;    //Factor of 10 million needed to convert!
243              Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy() + digitEnergy;
244              Float_t unitEnergyNoCuts = fUnitNoCuts[towerID-1].GetUnitEnergy() + digitEnergy;
245              fUnit[towerID-1].SetUnitEnergy(unitEnergy);
246              fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts);
247             }//end digits loop
248            }//end digits condition
249
250          //Set all unit flags, Eta, Phi
251          for(Int_t i=0; i<fNumUnits; i++)
252            {
253              if (fDebug>10) Info("FillUnitArray","Setting all units outside jets");
254              fUnit[i].SetUnitFlag(kOutJet);           //Set all units to be outside a jet initially
255              fUnit[i].SetUnitID(i+1);
256              Float_t eta;
257              Float_t phi;
258              geom->EtaPhiFromIndex(fUnit[i].GetUnitID(), eta, phi);
259              fUnit[i].SetUnitEta(eta);
260              fUnit[i].SetUnitPhi(phi*TMath::Pi()/180.0);
261
262              fUnitNoCuts[i].SetUnitFlag(kOutJet);           //Set all units to be outside a jet initially
263              fUnitNoCuts[i].SetUnitID(i+1);
264              eta = 0.0;
265              phi = 0.0;
266              geom->EtaPhiFromIndex(fUnitNoCuts[i].GetUnitID(), eta, phi);
267              fUnitNoCuts[i].SetUnitEta(eta);
268              fUnitNoCuts[i].SetUnitPhi(phi*TMath::Pi()/180.0);
269              //      if(i>13000) cout<<"!!!!!!!!!!!!!!!!!For unit0, eta="<<eta<<" and phi="<<phi*TMath::Pi()/180.0<<" and ID="<<fUnit[i].GetUnitID()<<endl;
270              //  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;
271            }//end loop over all units in array (same as all towers in EMCAL)
272    }
273
274
275  void AliEMCALJetFinderAlgoOmni::Sort(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t length)
276  {
277    //Calls the recursive quicksort method to sort unit objects in decending order of Energy
278    if (fDebug>1) Info("Sort","Sorting Unit objects");
279    QS(unit, 0, length-1);
280  }
281   
282
283  void AliEMCALJetFinderAlgoOmni::QS(AliEMCALJetFinderAlgoUA1Unit *unit, Int_t left, Int_t right)
284  {
285   //Sorts the AliEMCALJetFinderAlgoUA1Unit objects in decending order of Energy
286    if (fDebug>111) Info("QS","QuickSorting Unit objects");   
287
288    Int_t    i;
289    Int_t    j;
290    AliEMCALJetFinderAlgoUA1Unit  unitFirst;
291    AliEMCALJetFinderAlgoUA1Unit  unitSecond;
292
293    i = left;
294    j = right;
295    unitFirst = unit[(left+right)/2];
296
297  do
298   {
299     while( (unit[i].GetUnitEnergy() > unitFirst.GetUnitEnergy()) && (i < right)) i++;
300     while( (unitFirst.GetUnitEnergy() > unit[j].GetUnitEnergy()) && (j > left)) j--;
301
302     if(i <= j)
303       {
304         unitSecond = unit[i];
305         unit[i] = unit[j];
306         unit[j] = unitSecond;
307         i++;
308         j--;
309       }//end if
310   }while(i <= j);
311
312  if(left < j) QS(unit, left, j);
313  if(i < right) QS(unit, i, right);
314  }
315
316
317  void AliEMCALJetFinderAlgoOmni::FindBG()
318  {
319    if(fBGType == kRatio) RatioBG();
320    else if(fBGType == kCone) ConeBG();
321    else if(fBGType == kConstant) ConstantBG();
322  }
323
324  void AliEMCALJetFinderAlgoOmni::RatioBG()
325    {
326      //Finds the background energy for the iteration
327      //using the Ratio method
328      if (fDebug>1) Info("FindBG","Finding Average Background"); 
329      //Store BGEperCell from previous iteration!
330      fEBGTotalOld = fEBGTotal;
331      fEBGTotal          = 0.0;
332      Int_t numCone      = 0;
333
334      //If user has not set fBGPar, set it to the default
335      //for TPC = 90% efficiency, PtCut = 2GeV/c, timecut = 30ns
336      if(fBGPar == -1) fBGPar = 0.4685;
337
338      //Loop over all unit objects in the Unit array and link to same
339      //unit ID in NoCuts Unit array
340      for(Int_t i=0; i<fNumUnits; i++)
341        {
342          if(fUnit[i].GetUnitFlag() != kInJet)
343            {
344              Int_t id = fUnit[i].GetUnitID();
345              fEBGTotal += fUnitNoCuts[id-1].GetUnitEnergy();  
346            }
347          else numCone++;
348        }//end for
349
350      fEBGTotal *= fBGPar;
351      fEBGAve = fEBGTotal / (fNumUnits - numCone);
352      if (fDebug>5) Info("FindBG","Average BG is %f: ",fEBGAve); 
353
354      for(Int_t count=0; count<fNumUnits;count++)
355        {
356          fUnit[count].SetUnitFlag(kOutJet);
357        }//end for
358    }
359
360  void AliEMCALJetFinderAlgoOmni::ConeBG()
361    {
362      //Finds the background energy for the iteration
363      //using all energy not contained inside a jet
364      if (fDebug>1) Info("FindBG","Finding Average Background"); 
365      //Store old value of BGEperCell!
366      fEBGTotalOld = fEBGTotal;
367      fEBGTotal          = 0.0;
368      Int_t numCone      = 0;
369
370      //Loop over all unit objects in the array and sum the energy of those not in a jet
371      for(Int_t i=0; i<fNumUnits; i++)
372        {
373          if(fUnit[i].GetUnitFlag() != kInJet)
374            fEBGTotal += fUnit[i].GetUnitEnergy();
375          else numCone++;
376        }//end for
377
378      fEBGAve = fEBGTotal / (fNumUnits - numCone);
379      if (fDebug>5) Info("FindBG","Average BG is %f: ",fEBGAve);      
380
381      for(Int_t count=0; count<fNumUnits;count++)
382        {
383          fUnit[count].SetUnitFlag(kOutJet);
384        }//end for
385    }
386
387  void AliEMCALJetFinderAlgoOmni::ConstantBG()
388    {     
389      //Finds the background energy for the iteration
390      //using all energy not contained inside a jet
391      if (fDebug>1) Info("FindBG","Finding Average Background"); 
392
393      //If user has not set fBGPar, set it to the default
394      //for TPC = 90% efficiency, PtCut = 2GeV/c, timecut = 30ns
395      if(fBGPar == -1) fBGPar = 0.03378;
396
397      fEBGAve = fBGPar;
398      if (fDebug>5) Info("FindBG","Average BG is %f: ",fEBGAve);      
399
400      fEBGTotal          = 0.0;
401      Int_t numCone      = 0;
402      for(Int_t count=0; count<fNumUnits;count++)
403        {
404          if(fUnit[count].GetUnitFlag() == kInJet)
405            {
406              numCone++;
407            }
408          fUnit[count].SetUnitFlag(kOutJet);
409        }//end for
410      fEBGTotal = fEBGAve * (fNumUnits-numCone);
411      fEBGTotalOld = fEBGTotal;
412    }
413
414  void AliEMCALJetFinderAlgoOmni::FindJetEtaPhi(Int_t counter)
415    {
416      //Finds the eta and phi of the jet axis
417      if (fDebug>10) Info("FindJetEtaPhi","Finding Jet Eta and Phi");
418
419      fDEta = fUnit[counter].GetUnitEta() - fEtaInit;
420      fDPhi = fUnit[counter].GetUnitPhi() - fPhiInit;
421
422      fEnergy = fUnit[counter].GetUnitEnergy() - fEBGAve;
423      fJetEtaSum += fEnergy * fDEta;
424      fJetPhiSum += fEnergy * fDPhi;
425      fJetESum += fEnergy;
426      fJetEta = fEtaInit + (fJetEtaSum / fJetESum);
427      fJetPhi = fPhiInit + (fJetPhiSum / fJetESum);
428    }
429
430
431  void AliEMCALJetFinderAlgoOmni::FindJetEnergy()
432    {
433      //Finds the energy of the jet after the final axis has been found
434      if (fDebug>1) Info("FindJetEnergy","Finding Jet Energy");
435
436      for(Int_t i=0; i<fNumUnits; i++)
437        {
438          //Loop over all unit objects in the array and find if within cone radius
439          Float_t dEta = fUnit[i].GetUnitEta() - fJetEta;
440          Float_t dPhi = fUnit[i].GetUnitPhi() - fJetPhi;
441          Float_t rad = TMath::Sqrt( (dEta*dEta) + (dPhi*dPhi) );
442
443          if(fUnit[i].GetUnitFlag()==kOutJet && rad<= fConeRad)
444            {
445              fUnit[i].SetUnitFlag(kInCurrentJet);
446              Float_t energy = fUnit[i].GetUnitEnergy() - fEBGAve;
447              fJetESum += energy;                             
448              fJetEtaSum += energy * dEta;
449              fJetPhiSum += energy * dPhi;
450              fNumInCone++;                     //Increment the number of cells in the jet cone
451            }//end if
452        }//end for
453    }
454
455
456  void AliEMCALJetFinderAlgoOmni::StoreJetInfo()
457    {
458      //Stores the resulting jet information in appropriate storage structure (TO BE DECIDED!!!!)
459      if (fDebug>1) Info("StoreJetInfo","Storing Jet Information");
460     
461      //Store:
462      //fJetESum is the final jet energy (background has been subtracted)
463      //fJetEta is the final jet Eta
464      //fJetPhi is the final jet Phi
465      //fNumInCone is the final number of cells included in the jet cone
466      //fEtaInit is the eta of the initiator cell
467      //fPhiInit is the phi of the initiator cell
468      fJet.SetEnergy(fJetESum);
469      fJet.SetEta(fJetEta);
470      fJet.SetPhi(fJetPhi);
471
472       cout<<"For iteration "<<fNumIter <<" and Jet number " <<fNumJets <<endl;
473       cout<<"The jet energy is: " <<fJetESum <<endl;
474       cout<<"The jet eta is ---->" <<fJetEta <<endl;
475       cout<<"The jet phi is ---->" <<fJetPhi <<endl;
476
477      Int_t             numberTracks = fInputPointer->GetNTracks();
478      TParticle         *myP;
479      Int_t             numTracksInCone = 0;
480
481      for(Int_t counter=0; counter<numberTracks; counter++)
482        {
483         myP = fInputPointer->GetTrack(counter);
484         Float_t eta = myP->Eta();
485         Float_t phi = myP->Phi(); 
486         Float_t deta = fJetEta-eta;
487         Float_t dphi = fJetPhi -phi;
488         Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi));
489         if(rad<=fConeRad) numTracksInCone++;
490        }//end for
491
492      Float_t    *pTArray = new Float_t[numTracksInCone];
493      Float_t    *etaArray = new Float_t[numTracksInCone];
494      Float_t    *phiArray = new Float_t[numTracksInCone];
495      Int_t      *pdgArray = new Int_t[numTracksInCone];
496      Int_t             index = 0;
497
498      for(Int_t counter2=0; counter2<numberTracks; counter2++)
499        {
500          myP = fInputPointer->GetTrack(counter2);
501          Float_t eta = myP->Eta();
502          Float_t phi = myP->Phi(); 
503          Float_t deta = fJetEta-eta;
504          Float_t dphi = fJetPhi -phi;
505          Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi));
506          if(rad<=fConeRad)
507            {
508              pTArray[index] = myP->Pt();
509              etaArray[index] = eta;
510              phiArray[index] = phi;
511              pdgArray[index] = myP->GetPdgCode();
512              index++;
513            }//end if
514        }//end for
515
516      fJet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray);
517      fOutputObject.AddJet(&fJet);
518      delete[] pTArray;
519      delete[] etaArray;
520      delete[] phiArray;
521      delete[] pdgArray;
522    }
523
524
525  void AliEMCALJetFinderAlgoOmni::FindJets()
526    {
527      //Runs the complete UA1 JetFinding algorithm to find jets!
528      if (fDebug>1) Info("FindJets","Starting Jet Finding!!!");
529
530      //If the array of JetFinderUnit objects has not been initialised then initialise with default settings
531      if(!fArrayInitialised) 
532       {
533        InitUnitArray();
534        FillUnitArray(kFillAll);
535       }//end if
536      if (fDebug>1) Info("FindJets","Unit array filled");
537
538      //Step 1. Sort the array in descending order of Energy
539      Sort(fUnit,fNumUnits);
540
541      //Step 2. Set the number of iterations and Number of jets found to zero to start
542      fNumIter = 0;
543      fNumJets = 0;
544
545      //Step 3. Begin the iteration loop to find jets
546      //Need to iterate the algorithm while number of iterations<2 OR number of iterations<10 AND 
547      //the value of the average background has changed more than specified amount
548      //Min iterations = 2, Max iterations = 10
549      //while(fNumIter<2 || (fNumIter <10 && ( (fEBGTotal-fEBGTotalOld)/fEBGTotal) > fBGMaxMove) )
550
551      while(fNumIter<2 || (fNumIter <10 && ( fEBGTotal-fEBGTotalOld) > fEBGTotal*fBGMaxMove) )
552        {
553         if (fDebug>1) Info("FindJets","Starting BIG iteration ---> %i",fNumIter);
554
555          //Step 4. Find the value of the average background energy
556          FindBG();
557          fOutputObject.Reset(kResetJets); //Reset output object to store info for new iteration
558          fNumJets=0;
559
560          //Loop over the array of unit objects and flag those with energy below MinCellEt
561          Int_t numbelow = 0;
562          for(Int_t j=0; j<fNumUnits; j++)
563            {
564              if( (fUnit[j].GetUnitEnergy()-fEBGAve) < fEtMin)
565                 {       
566                   //          fUnit[j].SetUnitFlag(kBelowMinEt);    TAKING OUT kBelow flag
567                   numbelow++;
568                 }//end if
569            }//end for
570          //cout<<"THERE WERE "<<numbelow<<" units with E <EtMin!!!!!!!!!!!!!!!"<<endl;
571
572          //Do quick check if there are no jets upfront
573          // if(fUnit[0].GetUnitFlag() == kBelowMinEt)
574          if( (fUnit[0].GetUnitEnergy()-fEBGAve) < fEtMin)
575            {
576             cout <<"There are no jets for this event!" <<endl;
577             break;
578            }//end if
579
580          //Step 5. Begin with the first jet candidate cell (JET SEED LOOP)
581          if (fDebug>5) Info("FindJets","Beginning JET SEED LOOP");
582          for(Int_t count=0; count<fNumUnits; count++)
583            {
584
585 //CHECK CONDITION HERE _ NOT SURE IF SHOULD MAYBE BE: GetUnitEnergy()-fEBGAve >fESeed?????????????????????????????
586              if(fUnit[count].GetUnitEnergy()>=fESeed && fUnit[count].GetUnitFlag()==kOutJet)
587                {
588                  fEnergy = fUnit[count].GetUnitEnergy() - fEBGAve;
589                  fJetEta = fUnit[count].GetUnitEta();
590                  fJetPhi = fUnit[count].GetUnitPhi();
591                  Int_t seedID = fUnit[count].GetUnitID();
592                  if (fDebug>5) Info("FindJets","Inside first candidate jet seed loop for time : %i", count);
593                  if (fDebug>5) Info("FindJets","Found candidate energy %f ",fEnergy);
594                  if (fDebug>5) Info("FindJets","Found candidate eta %f ", fJetEta);
595                  if (fDebug>5) Info("FindJets","Found candidate phi %f ", fJetPhi);
596                  if (fDebug>5) Info("FindJets","Found candidate ID %i", seedID);
597
598                  fEtaInit = fJetEta;
599                  fPhiInit = fJetPhi;
600                  fEtaB = fJetEta;
601                  fPhiB = fJetPhi;
602                  fJetESum = 0.0;
603                  fJetEtaSum = 0.0;
604                  fJetPhiSum = 0.0;
605        
606          //Step 6. Find Jet Eta and Phi
607                  //Loop over all units in the array to find the ones in the jet cone and determine contrib to Jet eta, phi
608                  do
609                    {
610                      for(Int_t count1=0; count1<fNumUnits; count1++)               
611                        {
612                          if(fUnit[count1].GetUnitID() == seedID) continue;   //skip unit if the jetseed to avoid doublecounting
613                          if(fUnit[count1].GetUnitFlag() == kOutJet)
614                            {
615                              fDEta = fUnit[count1].GetUnitEta() - fJetEta;
616                              fDPhi = fUnit[count1].GetUnitPhi() - fJetPhi;
617                              fRad = TMath::Sqrt( (fDEta*fDEta) + (fDPhi*fDPhi) );
618                              if(fRad <= fConeRad)
619                                {
620                                  FindJetEtaPhi(count1); 
621                                }//end if
622                            }//end if
623                        }//end for (Jet Eta, Phi LOOP)
624                              
625                       //Find the distance cone centre moved from previous cone centre
626                       if (fDebug>10) Info("FindJets","Checking if cone move small enough");
627                       fDistP = TMath::Sqrt( ((fJetEta-fEtaB)*(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)*(fJetPhi-fPhiB)) );
628                       //     if(fDistP <= fMinMove) break;
629                              
630
631                       //Find the distance cone centre is from initiator cell
632                       if (fDebug>10) Info("FindJets","Checking if cone move too large");
633                       fDistI = TMath::Sqrt( ((fJetEtaSum/fJetESum)*(fJetEtaSum/fJetESum)) + ((fJetPhiSum/fJetESum)*
634                                                                                                     (fJetPhiSum/fJetESum)));
635
636                       if(fDistP>fMinMove && fDistI<fMaxMove)
637                         {
638                           fEtaB = fJetEta;
639                           fPhiB = fJetPhi;
640                         }//end if
641                  
642                    }while(fDistP>fMinMove && fDistI<fMaxMove);
643                           
644                  fJetEta = fEtaB;
645                  fJetPhi = fPhiB;
646
647
648        //Step 7. Find the Jet Energy
649                  if (fDebug>1) Info("FindJets","Looking for Jet energy");
650                  fJetESum = 0.0;
651                  fJetEtaSum = 0.0;
652                  fJetPhiSum = 0.0;
653                  fNumInCone = 0;
654                  FindJetEnergy();
655
656                  //cout<<"Number of cells in jet cone is: "<<fNumInCone<<endl;
657
658        //Step 8. Check if the jet is a valid jet
659                  //Check if cluster energy is above Min allowed to be a jet
660 //DID NOT DO THE COSH COMPARISON HERE -> NEED TO CHECK WHICH COMPARISON IS BEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661                  if (fDebug>5) Info("FindJets","Checking cluster is valid jet");
662                  if(fJetESum < fJetEMin)
663                    {
664                      for(Int_t count2=0; count2<fNumUnits; count2++)
665                        {
666                          if(fUnit[count2].GetUnitFlag()==kInCurrentJet || fUnit[count2].GetUnitFlag()==kOutJet)
667                            fUnit[count2].SetUnitFlag(kOutJet);
668                        }//end for
669                    if (fDebug>10) Info("FindJets","NOT a valid jet cell");
670                   }else
671                     {
672                      for(Int_t count2=0; count2<fNumUnits; count2++)
673                        {
674                          if(fUnit[count2].GetUnitFlag()==kInCurrentJet)
675                            {
676                              //      cout<<"Setting unit #"<<count2 <<" to be officially in a jet!"<<endl;
677                            fUnit[count2].SetUnitFlag(kInJet);
678                            }
679                        }//end for                       
680
681  //NEED TO CHECK FINAL WEIRD ITERATION OF ETA AND PHI CHANGES!!!!!!!!!
682                      //  fJetPhi += fJetPhiSum/fJetESum;        //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
683                      //  fJetEta += fJetEtaSum/fJetESum;        //CHECK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
684
685                      fNumJets++;              //Incrementing number of jets found
686                      StoreJetInfo();          //Storing jet info
687
688                  }//end if (check cluster above Min Jet Energy)
689                }//end if (Jet Seed condition)
690            }//end (JET SEED LOOP)
691
692 if (fDebug>5) Info("FindJets","End of BIG iteration number %i",fNumIter);
693 // this->Dump();
694          fNumIter++;
695        }//end 10 iteration WHILE LOOP
696  }
697
698
699
700
701
702
703
704
705
706
707