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