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