]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/jetan2004/AliTkConeJetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
SetInput call corrected for generated jets.
[u/mrichter/AliRoot.git] / JETAN / jetan2004 / AliTkConeJetFinder.cxx
CommitLineData
b9a6a391 1//$Id$
2/**************************************************************************
3 * AliTkConeJetFinder.cxx *
4 * Thorsten Kollegger <kollegge@ikf.physik.uni-frankfurt.de *
5 * Jet finder based on a cone algorithm with seeds and addition of *
6 * midpoints, for a description of the algorithm see hep-ex/0005012 *
7 *************************************************************************/
8
9// all includes in AliTkConeJetFinder.h
10#include "AliTkConeJetFinder.h"
11
12ClassImp(AliTkConeJetFinder)
13
14AliTkConeJetFinder::AliTkConeJetFinder() : TObject() {
15 setEtaGrid(-999,999,-999);
16 setPhiGrid(-999,999,-999);
17 hEtHist = NULL;
18}
19
20AliTkConeJetFinder::~AliTkConeJetFinder() {
21 if (hEtHist) {
22 delete hEtHist;
23 }
24}
25
26void AliTkConeJetFinder::setDefaultSettings() {
27 setEtaGrid(100,-5.,5.);
28 setPhiGrid((Int_t) (2*TMath::Pi()/0.1),
29 0,
30 2*TMath::Pi());
31 setJetConeRadius(0.7);
32}
33
34void AliTkConeJetFinder::Init() {
35 // create an array for the tower info
36 // and at this stage we know how many towers we have...
37 towers = NULL;
38 nTower = nEtaBins*nPhiBins;
39 towers = new AliTkTower[nTower];
40 if (!towers) {
41 cerr << "Couldn't allocate space for tower info!!!" << endl;
42 }
43 // let's create the towers;
44 Float_t etaBinSize = (nEtaMax-nEtaMin)/(Float_t)nEtaBins;
45 cout << "Init: etaBinSize=" << etaBinSize << endl;
46 Float_t phiBinSize = (nPhiMax-nPhiMin)/(Float_t)nPhiBins;
47 cout << "Init: phiBinSize=" << phiBinSize << endl;
48 //AliTkTower *tower;
49 Float_t eta = nEtaMin;
50 Float_t phi = nPhiMin;
51 for (Int_t uid = 0; uid < nTower; uid++) {
52 towers[uid].uid = uid;
53 towers[uid].etaMin = eta;
54 towers[uid].phiMin = phi;
55 towers[uid].etaMax = eta + etaBinSize;
56 towers[uid].phiMax = phi + phiBinSize;
57 towers[uid].Et = 0;
58 phi += phiBinSize;
59 if (phi > nPhiMax) {
60 phi = nPhiMin;
61 eta += etaBinSize;
62 }
63 }
64
65 // create an array for the seed list...
66 SeedTower = new Int_t[nTower];
67
68 // create the array to store the protojets
69 const Int_t expectedNumberOfProtoJets = 1000;
70 protojets = new TObjArray(expectedNumberOfProtoJets);
71 protojets->SetOwner(kTRUE);
72
73 // called once to initalize the finder...
74 // create the histogram which stores the Et information
75 hEtHist = new TH2D("Tower_EtHist","Tower_EtHist",
76 nEtaBins,nEtaMin,nEtaMax,
77 nPhiBins,nPhiMin,nPhiMax);
78 hEtHist->GetXaxis()->SetTitle("#eta");
79 hEtHist->GetYaxis()->SetTitle("#phi (rad)");
80 hEtHist->GetZaxis()->SetTitle("E_{t} in bins (GeV)");
81
82 // create histogram to hold Et cone centeres around tower center
83 hEtConeHist = new TH2D("Cone_EtHist","Cone_EtHist",
84 nEtaBins,nEtaMin,nEtaMax,
85 nPhiBins,nPhiMin,nPhiMax);
86 hEtConeHist->GetXaxis()->SetTitle("#eta");
87 hEtConeHist->GetYaxis()->SetTitle("#phi (rad)");
88 hEtConeHist->GetZaxis()->SetTitle("E_{t} in cones (GeV)");
89
90 // create histogram to show the stable protojets before merging
91 hStableProtoJetHist = new TH2D("StableProtoJetHist","StableProtoJetHist",
92 nEtaBins,nEtaMin,nEtaMax,
93 nPhiBins,nPhiMin,nPhiMax);
94 hStableProtoJetHist->GetXaxis()->SetTitle("#eta");
95 hStableProtoJetHist->GetYaxis()->SetTitle("#phi (rad)");
96 hStableProtoJetHist->GetZaxis()->SetTitle("E_{t} in stable cones (GeV)");
97
98 hTowerEt = new TH2D("tower_EtCheck","tower EtCheck",
99 nEtaBins,nEtaMin,nEtaMax,
100 nPhiBins,nPhiMin,nPhiMax);
101 hTowerEt->GetXaxis()->SetTitle("#eta");
102 hTowerEt->GetYaxis()->SetTitle("#phi (rad)");
103 hTowerEt->GetZaxis()->SetTitle("E_{t} in bins (GeV)");
104
105}
106
107void AliTkConeJetFinder::InitEvent() {
108 // called at each events
109 // clears the jet finder for the next event
110
111 // clear tower array - only Et
112 for (Int_t i = 0; i < nTower; i++) {
113 towers[i].Et = 0.;
114 }
115
116 // clear protojet list
117 protojets->Clear();
118
119 // clear control histograms
120 hEtHist->Reset();
121 hEtConeHist->Reset();
122 hTowerEt->Reset();
123 hStableProtoJetHist->Reset();
124
125}
126
127void AliTkConeJetFinder::FillEtHistFromTParticles(TClonesArray *particles) {
128 // input TClonesArray with TParticles, e.g. from PYTHIA
129 // fills the Et grid from these particles
130 if (!particles) {
131 return;
132 }
133 TParticle *particle = NULL;
134 Float_t Et = 0;
135
136 // loop over all particles...
137 TIterator *iter = particles->MakeIterator();
138 while ((particle = (TParticle *) iter->Next()) != NULL) {
139 // check if particle is accepted
140 if (isTParticleAccepted(particle)) {
141 // calculate Et for a massless particle...
142 Et = TMath::Sqrt(particle->Pt()*particle->Pt());
143 // idea for Et check: Et = E *sin(theta);
144 Int_t tower = findTower(particle->Eta(),particle->Phi());
145 if (tower >= 0) {
146 towers[tower].Et += Et;
147 } else {
148// cerr << "Couldn't find tower!!! eta="
149// << particle->Eta() << " phi=" << particle->Phi()
150// << endl;
151 }
152 // and to the histogram for control reasons
153 AddEtHist(particle->Eta(),particle->Phi(),Et);
154 }
155 }
156
157 // at this point we have filled all towers...
158 // let's make sure tower info is right
159 for (Int_t i = 0; i < nTower; i++) {
160 Float_t etacenter = (towers[i].etaMin + towers[i].etaMax)/2.;
161 Float_t phicenter = (towers[i].phiMin + towers[i].phiMax)/2.;
162 hTowerEt->Fill(etacenter,phicenter,towers[i].Et);
163 }
164
165}
166
167Bool_t AliTkConeJetFinder::isTParticleAccepted(TParticle *particle) {
168 // check if particle is accepted
169 // makes sense to write this into a own class, but now I'm lazy
170
171 // check if particle is stable -> a detectable particle
172 if (particle->GetStatusCode() != 1) {
173 return kFALSE;
174 }
175 // and now just for fun
176 TParticlePDG *partPDG;
177 partPDG = particle->GetPDG();
178// if (partPDG->Charge() == 0) {
179// cout << partPDG->GetName() << endl;
180// return kFALSE;
181// }
182
183 // default case: accept
184 return kTRUE;
185}
186
187void AliTkConeJetFinder::AddEtHist(Float_t eta, Float_t phi, Double_t Et) {
188 // add particle/tower withEt at position eta,phi to the current Et grid...
189 if (hEtHist) {
190 Double_t oldEt;
191 // get the Et already stored in the bin eta,phi
192 oldEt = hEtHist->GetBinContent(hEtHist->GetXaxis()->FindFixBin(eta),
193 hEtHist->GetYaxis()->FindFixBin(phi));
194 // add the new one
195 Et += oldEt;
196 // and set the bin to the sum of old+new Et.
197 hEtHist->SetBinContent(hEtHist->GetXaxis()->FindFixBin(eta),
198 hEtHist->GetYaxis()->FindFixBin(phi),
199 Et);
200 }
201}
202
203Int_t AliTkConeJetFinder::run() {
204 Int_t status;
205 status = FindProtoJets();
206 cout << "Found " << protojets->GetEntries() << " stable protojets"
207 << endl;
208 return status;
209}
210
211Int_t AliTkConeJetFinder::FindProtoJets() {
212 Float_t etaTower;
213 Float_t phiTower;
214
215 Float_t eta;
216 Float_t phi;
217 Float_t Et;
218
219 Float_t etaDiff;
220 Float_t phiDiff;
221
222 AliTkProtoJet protojet;
223
224 // let's take each tower as seed
225 // this means seedless search...
226 // introduce later step with seed+midpoints...
227 for (Int_t i=0; i < nTower; i++) {
228 // loop over all towers
229 etaTower = (towers[i].etaMax+towers[i].etaMin)/2.;
230 phiTower = (towers[i].phiMax+towers[i].phiMin)/2.;
231
232 eta = etaTower;
233 phi = phiTower;
234
235 // let's calculate a protojet with centered around tower center
236 if ((CalculateConeCentroid(&eta,&phi,&Et) > 0)) {
237 // ok, we found a valid centroid around tower pos
238 // let's fill the initial seed into the EtConeHist
239 hEtConeHist->Fill(etaTower,phiTower,Et);
240 // let's calculate the vector to the weighted cone center
241 etaDiff = (eta - etaTower)*(eta - etaTower);
242 phiDiff = calcPhiDiff(phi,phiTower);
243
244 // let's iterate the proto-jet while it's in the tower boundaries
245 // or it is stable
246 // or the maximum number of iterations is reached
247 Bool_t isStable = kFALSE;
248 Bool_t isInTower = isProtoJetInTower(i,eta,phi);
249 //Bool_t isInTower = isProtoJetInTower(etaBin,phiBin,eta,phi);
250 Int_t nIterations = 0;
251 Float_t etaOld;
252 Float_t phiOld;
253
254 // make const variables private members and set in default...
255 const Int_t nMaxIterations = 100;
256 while ((nIterations < nMaxIterations) &&
257 (!isStable) &&
258 (isInTower)) {
259 etaOld = eta;
260 phiOld = phi;
261 CalculateConeCentroid(&eta,&phi,&Et);
262 etaDiff = (eta - etaOld)*(eta - etaOld);
263 phiDiff = calcPhiDiff(phi,phiOld)*calcPhiDiff(phi,phiOld);
264 // and check protojet again...
265 isStable = isProtoJetStable(etaDiff,phiDiff);
266 // we allow the jet to flow away...
267 //isInTower = kTRUE;
268 // we force protojets to stay in tower
269 isInTower = isProtoJetInTower(i,eta,phi);
270 nIterations++;
271 }
272
273 if ((isStable) || (nIterations == nMaxIterations)) {
274// cout << "Found a stable protojet at eta=" << eta
275// << " phi=" << phi << " with Et=" << Et << endl;
276 if (nIterations == nMaxIterations) {
277 cout << "To many iterations on protojet" << endl;
278 }
279 // create a new protojet object
280 AliTkProtoJet *protojet = new AliTkProtoJet();
281 protojet->eta = eta;
282 protojet->phi = phi;
283 protojet->Et = Et;
284
285 // let's check if we allready found this protojet
286 TIter *iter = new TIter(protojets);
287 AliTkProtoJet *pj;
288 Bool_t wasFound = kFALSE;
289 while ((pj = (AliTkProtoJet *) iter->Next()) && !wasFound) {
290 wasFound = protojet->IsEqual(pj);
291 }
292
293 // and add it to the list of protojets...
294 if (!wasFound) {
295 protojets->Add(protojet);
296 //cout << protojet->Et << endl;
297 hStableProtoJetHist->Fill(eta,phi,Et);
298 } else {
299 delete protojet;
300 }
301 }
302 }
303 }
304 return 0;
305}
306
307Int_t AliTkConeJetFinder::FindJets() {
308 //let's look for the highest Et-Jet
309 //AliTkProtoJet *pj;
310
311 return 0;
312}
313
314
315Int_t AliTkConeJetFinder::CalculateConeCentroid(Float_t *eta, Float_t *phi,
316 Float_t *Et) {
317
318 //====================================================================
319 // WORKS STILL ON HISTOGRAM
320 // MUST WORK ON TOWERS
321 //====================================================================
322
323
324 // let's copy some stuff
325 Float_t etaCenter = *eta;
326 Float_t phiCenter = *phi;
327
328 /*
329 cout << "CalculateConeCentroid: Initial Cone Center eta="
330 << etaCenter
331 << " phi="
332 << phiCenter
333 << endl;
334 */
335
336 Float_t etaTower = 0;
337 Float_t phiTower = 0;
338 Float_t EtTower = 0;
339
340 Float_t etaDiff = 0;
341 Float_t phiDiff = 0;
342
343 Float_t etaJet = 0;
344 Float_t phiJet = 0;
345 Float_t EtJet = 0;
346
347 // lets calculate boundaries from jet cone
348 Float_t etaMin = etaCenter - jetRadius;
349 Float_t etaMax = etaCenter + jetRadius;
350 Float_t phiMin = phiCenter - jetRadius;
351
352 if (phiMin < 0) {
353 phiMin = 2*TMath::Pi() + phiMin;
354 }
355 Float_t phiMax = phiCenter + jetRadius;
356 if (phiMax > 2*TMath::Pi()) {
357 phiMax = phiMax - 2*TMath::Pi();
358 }
359
360 Int_t etaBinMin = hEtHist->GetXaxis()->FindFixBin(etaMin);
361 if (etaBinMin == 0) {
362 return -1;
363 }
364 Int_t etaBinMax = hEtHist->GetXaxis()->FindFixBin(etaMax);
365 if (etaBinMax == hEtHist->GetXaxis()->GetNbins()+1) {
366 return -1;
367 }
368 Int_t phiBinMin = hEtHist->GetYaxis()->FindFixBin(phiMin);
369 Int_t phiBinMax = hEtHist->GetYaxis()->FindFixBin(phiMax);
370
371 // let's loop over all bins/towers...
372 for (Int_t etaBin = etaBinMin; etaBin != etaBinMax+1; etaBin++) {
373 for (Int_t phiBin = phiBinMin; phiBin != phiBinMax+1; phiBin++) {
374 // roll over in phi...
375 if (phiBin == hEtHist->GetYaxis()->GetNbins()+1) {
376 phiBin = 0;
377 continue;
378 }
379 // let's get the tower
380 etaTower = hEtHist->GetXaxis()->GetBinCenter(etaBin);
381 phiTower = hEtHist->GetYaxis()->GetBinCenter(phiBin);
382 EtTower = hEtHist->GetCellContent(etaBin,phiBin);
383
384 // let's check if the bin is in the cone
385 etaDiff = (etaTower - etaCenter)*(etaTower - etaCenter);
386 phiDiff = calcPhiDiff(phiTower,phiCenter)*
387 calcPhiDiff(phiTower,phiCenter);
388 if (TMath::Sqrt(etaDiff+phiDiff) < jetRadius) {
389 // tower is in cone - add it to the cone
390 EtJet += EtTower;
391 if (EtTower != 0) {
392 etaJet += EtTower * etaTower;
393 phiJet += EtTower * phiTower;
394 }
395 }
396 }
397 }
398
399 // normalize eta and phi
400 if (EtJet != 0) {
401 etaJet /= EtJet;
402 phiJet /= EtJet;
403 }
404
405 /*
406 cout << "CalculateConeCentroid: Final Cone Center eta="
407 << etaJet
408 << " phi="
409 << phiJet
410 << " Et="
411 << EtJet
412 << endl;
413 */
414
415 // return the jet cone parameters
416 *eta = etaJet;
417 *phi = phiJet;
418 *Et = EtJet;
419
420
421 return 1;
422}
423
424Float_t AliTkConeJetFinder::calcPhiDiff(Float_t phi1, Float_t phi2) {
425 Float_t phidiff1 = TMath::Sqrt((phi1-phi2)*(phi1-phi2));
426 Float_t phidiff2 = TMath::Sqrt((phi1-phi2-2*TMath::Pi())*
427 (phi1-phi2-2*TMath::Pi()));
428 Float_t phidiff3 = TMath::Sqrt((phi1-phi2+2*TMath::Pi())*
429 (phi1-phi2+2*TMath::Pi()));
430 // return the minmum;
431 Float_t phiret = phidiff1;
432 if (phiret > phidiff2) {
433 phiret = phidiff2;
434 }
435 if (phiret > phidiff3) {
436 phiret = phidiff3;
437 }
438 return phiret;
439}
440
441Bool_t AliTkConeJetFinder::isProtoJetStable(Float_t etaDiff, Float_t phiDiff) {
442 const Float_t epsilon = 0.00001;
443 if (TMath::Sqrt(etaDiff*phiDiff) < epsilon) {
444 return kTRUE;
445 }
446 return kFALSE;
447}
448
449Bool_t AliTkConeJetFinder::isProtoJetInTower(Int_t etaBin, Int_t phiBin,
450 Float_t eta, Float_t phi) {
451 Float_t etaBinMax = hEtHist->GetXaxis()->GetBinUpEdge(etaBin);
452 Float_t etaBinMin = hEtHist->GetXaxis()->GetBinLowEdge(etaBin);
453 Float_t phiBinMax = hEtHist->GetYaxis()->GetBinUpEdge(phiBin);
454 Float_t phiBinMin = hEtHist->GetYaxis()->GetBinLowEdge(phiBin);
455 if ((eta < etaBinMax) &&
456 (eta > etaBinMin) &&
457 (phi < phiBinMax) &&
458 (phi > phiBinMin)) {
459 return kTRUE;
460 }
461
462 return kFALSE;
463}
464
465Bool_t AliTkConeJetFinder::isProtoJetInTower(Int_t tower, Float_t eta,
466 Float_t phi) {
467 if ((eta > towers[tower].etaMin) &&
468 (eta < towers[tower].etaMax) &&
469 (phi > towers[tower].phiMin) &&
470 (phi < towers[tower].phiMax)) {
471 return kTRUE;
472 }
473 return kFALSE;
474}
475
476Int_t AliTkConeJetFinder::findTower(Float_t eta, Float_t phi) {
477 Float_t etaBinSize = (nEtaMax - nEtaMin)/(Float_t)nEtaBins;
478 Float_t phiBinSize = (nPhiMax - nPhiMin)/(Float_t)nPhiBins;
479
480 if ((eta < nEtaMin) || (eta > nEtaMax)) {
481 return -1;
482 }
483 if ((phi < nPhiMin) || (phi > nPhiMax)) {
484 return -1;
485 }
486
487 Int_t etaBin = (Int_t) ((eta - nEtaMin)/etaBinSize);
488 Int_t phiBin = (Int_t) ((phi - nPhiMin)/phiBinSize);
489
490 Int_t bin = etaBin * nPhiBins + phiBin;
491
492 if (bin >= nTower) {
493 return -2;
494 }
495
496 return bin;
497}
498
499// lot of setters/getters... nothing special, no comments...
500void AliTkConeJetFinder::setEtaNBins(Int_t nbins) {
501 nEtaBins = nbins;
502}
503
504Int_t AliTkConeJetFinder::getEtaNBins() {
505 return nEtaBins;
506}
507
508void AliTkConeJetFinder::setEtaRange(Float_t min, Float_t max) {
509 nEtaMin = min;
510 nEtaMax = max;
511}
512
513Float_t AliTkConeJetFinder::getEtaRangeMin() {
514 return nEtaMin;
515}
516
517Float_t AliTkConeJetFinder::getEtaRangeMax() {
518 return nEtaMax;
519}
520
521void AliTkConeJetFinder::setEtaGrid(Int_t nbins, Float_t min, Float_t max) {
522 setEtaNBins(nbins);
523 setEtaRange(min,max);
524}
525
526void AliTkConeJetFinder::setPhiNBins(Int_t nbins) {
527 nPhiBins = nbins;
528}
529
530Int_t AliTkConeJetFinder::getPhiNBins() {
531 return nPhiBins;
532}
533
534void AliTkConeJetFinder::setPhiRange(Float_t min, Float_t max) {
535 nPhiMin = min;
536 nPhiMax = max;
537}
538
539Float_t AliTkConeJetFinder::getPhiRangeMin() {
540 return nPhiMin;
541}
542
543Float_t AliTkConeJetFinder::getPhiRangeMax() {
544 return nPhiMax;
545}
546
547void AliTkConeJetFinder::setPhiGrid(Int_t nbins, Float_t min, Float_t max) {
548 setPhiNBins(nbins);
549 setPhiRange(min,max);
550}
551
552void AliTkConeJetFinder::setJetConeRadius(Float_t r) {
553 jetRadius = r;
554}
555
556Float_t AliTkConeJetFinder::getJetConeRadius() {
557 return jetRadius;
558}
559
560TObjArray *AliTkConeJetFinder::getProtoJetList() {
561 return protojets;
562}
563
564//==========================================================================
565ClassImp(AliTkTower)
566
567//==========================================================================
568ClassImp(AliTkProtoJet)
569
570Bool_t AliTkProtoJet::IsEqual(AliTkProtoJet *other) {
571 const Float_t epsilon = 0.0000001;
572
573 Float_t etaDiff = (other->eta - eta)*(other->eta - eta);
574 Float_t phiDiff = (other->phi - phi)*(other->phi - phi);
575 Float_t EtDiff = (other->Et - Et)*(other->Et - Et);
576
577 Bool_t isEqual = kFALSE;
578 if ((etaDiff < epsilon) &&
579 (phiDiff < epsilon) &&
580 (EtDiff < epsilon)) {
581 isEqual = kTRUE;
582 }
583 return isEqual;
584}
585
586bool SProtoJet::operator==(const SProtoJet &s1) {
587 Float_t epsilon = 0.0000001;
588 bool b;
589 if (eta > s1.eta) {
590 b = ((eta - s1.eta) < epsilon);
591 } else {
592 b = ((s1.eta - eta) < epsilon);
593 }
594 if (phi > s1.phi) {
595 b = b && ((phi - s1.phi) < epsilon);
596 } else {
597 b = b && ((s1.phi - phi) < epsilon);
598 }
599 if (Et > s1.Et) {
600 b = b && ((Et - s1.Et) < epsilon);
601 } else {
602 b = b && ((s1.Et - Et) < epsilon);
603 }
604 return b;
605}