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