JETAN/AliTkConeJetFinder.cxx

   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
  12 ClassImp(AliTkConeJetFinder)
  13
  14 AliTkConeJetFinder::AliTkConeJetFinder() : TObject() {
  15   setEtaGrid(-999,999,-999);
  16   setPhiGrid(-999,999,-999);
  17   hEtHist = NULL;
  18 }
  19
  20 AliTkConeJetFinder::~AliTkConeJetFinder() {
  21   if (hEtHist) {
  22     delete hEtHist;
  23   }
  24 }
  25
  26 void 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
  34 void 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
 107 void 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
 127 void 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
 167 Bool_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
 187 void 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
 203 Int_t AliTkConeJetFinder::run() {
 204   Int_t status;
 205   status = FindProtoJets();
 206   cout << "Found " << protojets->GetEntries() << " stable protojets"
 207        << endl;
 208   return status;
 209 }
 210
 211 Int_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
 307 Int_t AliTkConeJetFinder::FindJets() {
 308   //let's look for the highest Et-Jet
 309   //AliTkProtoJet *pj;
 310
 311   return 0;
 312 }
 313
 314
 315 Int_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
 424 Float_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
 441 Bool_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
 449 Bool_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
 465 Bool_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
 476 Int_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...
 500 void AliTkConeJetFinder::setEtaNBins(Int_t nbins) {
 501   nEtaBins = nbins;
 502 }
 503
 504 Int_t AliTkConeJetFinder::getEtaNBins() {
 505   return nEtaBins;
 506 }
 507
 508 void AliTkConeJetFinder::setEtaRange(Float_t min, Float_t max) {
 509   nEtaMin = min;
 510   nEtaMax = max;
 511 }
 512
 513 Float_t AliTkConeJetFinder::getEtaRangeMin() {
 514   return nEtaMin;
 515 }
 516
 517 Float_t AliTkConeJetFinder::getEtaRangeMax() {
 518   return nEtaMax;
 519 }
 520
 521 void AliTkConeJetFinder::setEtaGrid(Int_t nbins, Float_t min, Float_t max) {
 522   setEtaNBins(nbins);
 523   setEtaRange(min,max);
 524 }
 525
 526 void AliTkConeJetFinder::setPhiNBins(Int_t nbins) {
 527   nPhiBins = nbins;
 528 }
 529
 530 Int_t AliTkConeJetFinder::getPhiNBins() {
 531   return nPhiBins;
 532 }
 533
 534 void AliTkConeJetFinder::setPhiRange(Float_t min, Float_t max) {
 535   nPhiMin = min;
 536   nPhiMax = max;
 537 }
 538
 539 Float_t AliTkConeJetFinder::getPhiRangeMin() {
 540   return nPhiMin;
 541 }
 542
 543 Float_t AliTkConeJetFinder::getPhiRangeMax() {
 544   return nPhiMax;
 545 }
 546
 547 void AliTkConeJetFinder::setPhiGrid(Int_t nbins, Float_t min, Float_t max) {
 548   setPhiNBins(nbins);
 549   setPhiRange(min,max);
 550 }
 551
 552 void AliTkConeJetFinder::setJetConeRadius(Float_t r) {
 553   jetRadius = r;
 554 }
 555
 556 Float_t AliTkConeJetFinder::getJetConeRadius() {
 557   return jetRadius;
 558 }
 559
 560 TObjArray *AliTkConeJetFinder::getProtoJetList() {
 561   return protojets;
 562 }
 563
 564 //==========================================================================
 565 ClassImp(AliTkTower)
 566
 567 //==========================================================================
 568 ClassImp(AliTkProtoJet)
 569
 570 Bool_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
 586 bool 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 }