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