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