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