JETAN/AliTkConeJetFinderV2.cxx

   1 // $Id$
   2
   3 #include <Riostream.h>
   4 #include <vector>
   5 #include <list>
   6 #include <map>
   7
   8 #include <TROOT.h>
   9 #include <TClonesArray.h>
  10 #include <TParticle.h>
  11 #include <TH2.h>
  12 #include <TCanvas.h>
  13 #include <TFile.h>
  14 #include <TTree.h>
  15 #include <TString.h>
  16
  17 #include "AliTkConeJet.h"
  18 #include "AliTkTowerV2.h"
  19 #include "AliTkConeJetEvent.h"
  20 #include "AliTkConeJetFinderV2.h"
  21
  22 #ifdef ALICEINTERFACE
  23 #include <AliJetParticle.h>
  24 #include <AliJetEvent.h>
  25 #include <AliJetEventParticles.h>
  26 #endif
  27
  28 //numerical precision
  29 #define __epsilon__ 1e-6
  30
  31 //==== Implementation of tower =====
  32 tower::tower() : fEta_min(-999),fEta_max(-999),fEta_center(-999),
  33                  fPhi_min(-999),fPhi_max(-999),fPhi_center(-999),
  34                  fEt(0), fParticles(0)
  35 {
  36 }
  37
  38 tower::tower(const tower& t)
  39 {
  40   fEta_min = t.getEtaMin();
  41   fEta_max = t.getEtaMax();
  42   fEta_center = t.getEta();
  43   fPhi_min = t.getPhiMin();
  44   fPhi_max = t.getPhiMax();
  45   fPhi_center = t.getPhi();
  46   fEt = t.getEt();
  47 }
  48
  49 tower::tower(Float_t phimin, Float_t phimax, Float_t etamin, Float_t etamax)
  50 {
  51   fEta_min = etamin;
  52   fEta_max = etamax;
  53   fEta_center = (etamax+etamin)/2.;
  54   fPhi_min = phimin;
  55   fPhi_max = phimax;
  56   fPhi_center = (phimax+phimin)/2.;
  57   fEt = 0;
  58 }
  59
  60 tower::~tower()
  61 {
  62   clearParticles();
  63 }
  64
  65 tower& tower::operator+=(const Float_t E)
  66 {
  67   fEt += E;
  68   return *this;
  69 }
  70
  71 tower& tower::operator+=(const TParticle *part)
  72 {
  73   fEt += part->Pt();   // Et for a massless particle...
  74   addParticle(part);
  75   return *this;
  76 }
  77
  78 ostream& operator<<(ostream& s,const tower& t)
  79 {
  80   return s << "tower info: etamin=" << t.getEtaMin()
  81            << " etamax=" << t.getEtaMax()
  82            << " phimin=" << t.getPhiMin()
  83            << " phimax=" << t.getPhiMax()
  84            << " Et=" << t.getEt();
  85 }
  86
  87 list<const TParticle*> *tower::getParticles()
  88 {
  89   list<const TParticle*> *newList = new list<const TParticle*>;
  90   list<const TParticle*> &myList = *newList;
  91   copy(fParticles.begin(),fParticles.end(),back_inserter(myList));
  92   return newList;
  93 }
  94
  95 //==== Implementation of protojet =====
  96 protojet:: protojet() : fCentroid(-999,-999),
  97                         fEtWeightedCentroid(-999,-999),
  98                         fEt(-999),fUpdate(kFALSE),fTowers(0)
  99 {
 100 }
 101
 102 protojet::protojet(const protojet& p) :
 103                         fCentroid(-999,-999),
 104                         fEtWeightedCentroid(-999,-999),
 105                         fEt(-999),fUpdate(kFALSE),fTowers(0)
 106 {
 107   setCentroidPosition(p.getCentroidPosition());
 108 }
 109
 110 protojet::protojet(const protojet *p) :
 111                         fCentroid(-999,-999),
 112                         fEtWeightedCentroid(-999,-999),
 113                         fEt(-999),fUpdate(kFALSE),fTowers(0)
 114 {
 115   setCentroidPosition(p->getCentroidPosition());
 116 }
 117
 118 bool protojet::operator==(const protojet &p1)
 119 {
 120   AliTkEtaPhiVector otherCenter = p1.getCentroidPosition();
 121   if (fCentroid.diffSq(otherCenter) > __epsilon__)
 122     return false;
 123
 124   if ((getEt() - p1.getEt())*(getEt() - p1.getEt()) > __epsilon__)
 125    return false;
 126
 127   return true;
 128 }
 129
 130 void protojet::update()
 131 {
 132   Float_t Et = 0;
 133   Float_t phi = 0;
 134   Float_t eta = 0;
 135
 136   for (list<tower*>::const_iterator iter = fTowers.begin();
 137        iter != fTowers.end(); ++iter) {
 138     Et += (*iter)->getEt();
 139     phi += (*iter)->getEt() * (*iter)->getPhi();
 140     eta += (*iter)->getEt() * (*iter)->getEta();
 141   }
 142   if (Et != 0) {
 143     phi /= Et;
 144     eta /= Et;
 145   }
 146   fEtWeightedCentroid.setVector(eta,phi);
 147   fEt=Et;
 148   fUpdate=kFALSE;
 149 }
 150
 151 list<tower*> protojet::getTowerList() const
 152 {
 153   list<tower*> newList;
 154   copy(fTowers.begin(),fTowers.end(),back_inserter(newList));
 155   return newList;
 156 }
 157
 158 list<tower*> protojet::getSharedTowerList(const protojet *other) const
 159 {
 160   list<tower*> newList;
 161   for (list<tower*>::const_iterator i = fTowers.begin();
 162        i != fTowers.end(); ++i) {
 163     if (other->hasTower((*i))) {
 164       newList.push_back((*i));
 165     }
 166   }
 167   return newList;
 168 }
 169
 170 bool protojet::hasTower(tower *pTower) const
 171 {
 172   bool bHasTower = false;
 173   for (list<tower*>::const_iterator iter = fTowers.begin();
 174        iter != fTowers.end(); ++iter) {
 175     if (pTower == (*iter)) {
 176       bHasTower = true;
 177       break;
 178     }
 179   }
 180   return bHasTower;
 181 }
 182
 183 bool protojet::shareTowers(const protojet *other) const
 184 {
 185   bool bShareTowers = false;
 186   for (list<tower*>::const_iterator iter = fTowers.begin();
 187        iter != fTowers.end(); ++iter) {
 188     if (other->hasTower(*iter)) {
 189       bShareTowers = true;
 190       break;
 191     }
 192   }
 193   return bShareTowers;
 194 }
 195
 196 Float_t protojet::diffToCenter(tower *pTower) const
 197 {
 198   AliTkEtaPhiVector v(pTower->getEta(),pTower->getPhi());
 199   AliTkEtaPhiVector center = getCentroidPosition();
 200
 201   return v.diff(center);
 202 }
 203
 204 ostream& operator<<(ostream& s,const protojet& p)
 205 {
 206   return s << "Protojet info: eta=" << p.Eta()
 207            << " phi=" << p.Phi() << " Et=" << p.getEt();
 208 }
 209
 210
 211 //================== implementation of AliTkConeJetFinderV2 ==============
 212 ClassImp(AliTkConeJetFinderV2)
 213
 214 AliTkConeJetFinderV2::AliTkConeJetFinderV2()
 215                   : TObject(),
 216                     fOutput(0),fNTowers(0),fEtaBins(0),fEtaMin(0),
 217                     fEtaMax(0),fEtaWidth(0),fPhiBins(0),fPhiMin(0),fPhiMax(0),
 218                     fPhiWidth(0),fEtCut(0),fEtMinJet(0),fPtCut(0),fRadius(0),fTowers(0),
 219                     fSeedPointsNew(),fProtojets(),fJets(0),
 220                     fEvoutfile(0),fEvout_name(0),fEvoutevent(0),fEvouttree(0)
 221 #ifdef ALICEINTERFACE
 222                     ,fAliParticles(0)
 223 #endif
 224
 225 #ifdef DOHISTOS
 226                     ,fHistos(),fEventHistos(),fHistFile(0)
 227 #endif
 228 {
 229   defaultSettings();
 230 }
 231
 232
 233 AliTkConeJetFinderV2::~AliTkConeJetFinderV2()
 234 {
 235   if(fEvoutevent) delete fEvoutevent;
 236   if(fEvoutfile) delete fEvoutfile;
 237   if(fEvout_name) delete[] fEvout_name;
 238 #ifdef ALICEINTERFACE
 239   if(fAliParticles) delete fAliParticles;
 240 #endif
 241
 242 #ifdef DOHISTOS
 243   if(fHistFile) delete fHistFile;
 244 #endif
 245 }
 246
 247 void  AliTkConeJetFinderV2::defaultSettings()
 248 {
 249   fOutput = kFALSE;
 250   fRadius = 0.7;
 251   // some more or less usefull default settings
 252   // only possibilty to control the jet finder without setters/getters
 253   fPhiBins = (Int_t) (2 * TMath::Pi() / 0.1);
 254   fPhiMin = 0;
 255   fPhiMax = 2 * TMath::Pi();
 256   fEtaBins = 20;
 257   fEtaMin = -1;
 258   fEtaMax = 1;
 259   fEtCut = 0;
 260   fPtCut = 0;
 261   fEtMinJet = 0;
 262   fNTowers = fPhiBins*fEtaBins;
 263 }
 264
 265 void AliTkConeJetFinderV2::setSettings(Int_t phibins,Int_t etabins)
 266 {
 267   fPhiBins = phibins;
 268   fPhiMin = 0;
 269   fPhiMax = 2 * TMath::Pi();
 270   fEtaBins = etabins;
 271   fEtaMin = -1;
 272   fEtaMax = 1;
 273   fNTowers = fPhiBins*fEtaBins;
 274 }
 275
 276 void AliTkConeJetFinderV2::init()
 277 {
 278   createTowers();
 279
 280   fEvoutevent = new AliTkConeJetEvent();
 281   fEvoutfile = 0;
 282   fEvouttree = 0;
 283
 284 #ifdef DOHISTOS
 285   createHistos();
 286   createEventHistos();
 287 #endif
 288
 289   if (getEvOutFilename()) {
 290     if(fOutput)cout << "Writing finder output to " << getEvOutFilename() << endl;
 291     fEvoutfile = new TFile(getEvOutFilename(),"RECREATE");
 292     fEvouttree = new TTree("jets","TKConeJetFinderV2 jets");
 293     fEvouttree->Branch("ConeFinder","AliTkConeJetEvent",&(this->fEvoutevent),32000,0);
 294   }
 295 }
 296
 297 #ifdef ALICEINTERFACE
 298 void AliTkConeJetFinderV2::initEvent(const AliJetEventParticles *p,TString desc)
 299 {
 300   if(fEvoutevent){
 301     fEvoutevent->Clear();
 302     fEvoutevent->setDesc(desc);
 303     //fEvoutevent->setJetParticles(p);
 304     fEvoutevent->setJetParticles(0);
 305   }
 306   const TClonesArray *parts=p->GetParticles();
 307   initEvent_(parts,2);
 308 }
 309 #endif
 310
 311 void AliTkConeJetFinderV2::initEvent(const TClonesArray *particles,Int_t type,TString desc)
 312 {
 313   if(fEvoutevent){
 314     fEvoutevent->Clear();
 315     fEvoutevent->setDesc(desc);
 316   }
 317   initEvent_(particles,type);
 318 }
 319
 320 void AliTkConeJetFinderV2::initEvent(const TClonesArray *particles,Int_t type)
 321 {
 322 #ifdef DOHISTOS
 323   // clear event histograms for new event
 324   clearEventHistos();
 325 #endif
 326
 327   if(fEvoutevent){
 328     fEvoutevent->Clear();
 329   }
 330
 331   initEvent_(particles,type);
 332 }
 333
 334
 335 void AliTkConeJetFinderV2::initEvent_(const TClonesArray *particles,Int_t type)
 336 {
 337   fEvoutevent->setRadius(fRadius);
 338   fEvoutevent->setEtCut(fEtCut);
 339   fEvoutevent->setPtCut(fPtCut);
 340
 341   // reset all towers
 342   for (vector<tower>::iterator i = fTowers->begin(); i != fTowers->end();++i) {
 343      (*i).clear();
 344   }
 345   // reset seed points
 346   fSeedPointsNew.erase(fSeedPointsNew.begin(),fSeedPointsNew.end());
 347   // reset protojet list
 348   for (list<protojet*>::iterator i = fProtojets.begin();i != fProtojets.end(); ++i) {
 349     protojet *p = (*i);
 350     if (p) {
 351       delete p;
 352     }
 353   }
 354   fProtojets.erase(fProtojets.begin(),fProtojets.end());
 355
 356   // reset jet list
 357   for (list<protojet*>::iterator i = fJets.begin();
 358        i != fJets.end(); ++i) {
 359     protojet *p = (*i);
 360     if (p) {
 361       delete p;
 362     }
 363   }
 364   fJets.erase(fJets.begin(),fJets.end());
 365
 366   // fill Et towers from particles
 367   switch (type) {
 368   case 1:
 369     fillTowersFromTParticles(particles);
 370     break;
 371 #ifdef ALICEINTERFACE
 372   case 2:
 373     fillTowersFromAliParticles(particles);
 374     break;
 375 #endif
 376   default:
 377     cerr << "AliTkConeJetFinderV2: don't know how to fill Et hist from TClonesArray with type " << type << endl;
 378   }
 379   createSeedPoints();
 380
 381 #ifdef DOHISTOS
 382   fillTowerHist();
 383 #endif
 384 }
 385
 386 Bool_t AliTkConeJetFinderV2::isJetEnergy(Float_t min, Float_t max)
 387 { // Checks if there was an jet with Et between min and max
 388   Float_t men=maxJetEnergy();
 389
 390   if((men>min)&&(men<max)) return kTRUE;
 391   return kFALSE;
 392 }
 393
 394 Float_t AliTkConeJetFinderV2::maxJetEnergy()
 395 {
 396   Float_t men=-1;
 397   for (list<protojet*>::const_iterator iter = fJets.begin();
 398        iter != fJets.end(); ++iter) {
 399
 400     Float_t et=(*iter)->getEt();
 401     if(men<et) men=et;
 402   }
 403   return men;
 404 }
 405
 406 void AliTkConeJetFinderV2::run()
 407 {
 408   findProtojets();
 409   findJets();
 410 }
 411
 412 void AliTkConeJetFinderV2::finishEvent()
 413 {
 414   if (fEvoutevent){
 415     AliTkConeJet *jet = new AliTkConeJet();
 416     AliTkTowerV2 *myTower = new AliTkTowerV2();
 417     for (list<protojet*>::const_iterator iter = fJets.begin();
 418          iter != fJets.end(); ++iter) {
 419
 420       Float_t jetet=(*iter)->getEt();
 421       if(jetet<fEtMinJet) continue;
 422       jet->Clear();
 423       jet->setEta((*iter)->Eta());
 424       jet->setPhi((*iter)->Phi());
 425       jet->setEt(jetet);
 426
 427       list<tower*> jettowers = (*iter)->getTowerList();
 428       for (list<tower*>::const_iterator twiter = jettowers.begin();
 429            twiter != jettowers.end(); ++twiter) {
 430         myTower->Clear();
 431         myTower->setEta((*twiter)->getEta());
 432         myTower->setPhi((*twiter)->getPhi());
 433         myTower->setEt((*twiter)->getEt());
 434         list<const TParticle*> &twparts = *((*twiter)->getParticles());
 435         for (list<const TParticle*>::const_iterator partiter = twparts.begin();
 436             partiter != twparts.end(); ++partiter) {
 437           myTower->addParticle(*partiter);
 438         }
 439         delete &twparts;
 440         jet->addTower(myTower);
 441       }
 442       jet->calculateValues();
 443       fEvoutevent->addJet(jet);
 444     }
 445     fEvoutevent->sortJets();
 446
 447     if (fEvouttree)
 448       fEvouttree->Fill();
 449     delete jet;
 450     delete myTower;
 451   }
 452 #ifdef DOHISTOS
 453   writeEventHistos();
 454 #endif
 455 }
 456
 457 void AliTkConeJetFinderV2::finish()
 458 {
 459   if ((fEvoutfile) && (fEvoutfile->IsOpen())) {
 460     fEvoutfile->cd();
 461     fEvouttree->Write();
 462     fEvoutfile->Close();
 463   }
 464
 465 #ifdef DOHISTOS
 466   writeHistos();
 467   fHistFile->Close();
 468 #endif
 469 }
 470
 471 void AliTkConeJetFinderV2::createTowers()
 472 {
 473   if(fOutput) cout << "Creating " << fNTowers << " tower" << endl;
 474   if (fPhiMax > fPhiMin) {
 475     fPhiWidth = (fPhiMax-fPhiMin)/(Float_t)fPhiBins;
 476   } else {
 477     fPhiWidth = (fPhiMin-fPhiMax)/(Float_t)fPhiBins;
 478   }
 479   if (fEtaMax > fEtaMin) {
 480     fEtaWidth = (fEtaMax-fEtaMin)/(Float_t)fEtaBins;
 481   } else {
 482     fEtaWidth = (fEtaMin-fEtaMax)/(Float_t)fEtaBins;
 483   }
 484   if(fOutput) cout << "Delta(eta)=" << fEtaWidth << " Delta(phi)=" << fPhiWidth << endl;
 485
 486   // let's create the container
 487   fTowers = new vector<tower>(fNTowers);
 488
 489   // set tower boundaries...
 490   Float_t eta = fEtaMin;
 491   Float_t phi = fPhiMin;
 492   for (Int_t e = 0; e < fEtaBins; e++) {
 493     for (Int_t p = 0; p < fPhiBins; p++) {
 494       Int_t tower = e*fPhiBins + p;
 495       (*fTowers)[tower].setPhiMin(phi);
 496       (*fTowers)[tower].setPhiMax(phi + fPhiWidth);
 497       (*fTowers)[tower].setPhi(phi + fPhiWidth/2.);
 498       (*fTowers)[tower].setEtaMin(eta);
 499       (*fTowers)[tower].setEtaMax(eta + fEtaWidth);
 500       (*fTowers)[tower].setEta(eta + fEtaWidth/2.);
 501
 502       phi += fPhiWidth;
 503     }
 504     eta += fEtaWidth;
 505     phi = fPhiMin;
 506   }
 507 }
 508
 509 void AliTkConeJetFinderV2::fillTowersFromTParticles(const TClonesArray *particles)
 510 {
 511   // input TClonesArray with TParticles, e.g. from PYTHIA
 512   // fills the Et grid from these particles
 513   if (!particles) {
 514     return;
 515   }
 516
 517   // loop over all particles...
 518   TParticle *particle = NULL;
 519   TIterator *iter = particles->MakeIterator();
 520   while ((particle = (TParticle *) iter->Next()) != NULL) {
 521     // check if particle is accepted
 522     if (1) { //isTParticleAccepted(particle)) {
 523       Float_t pt=particle->Pt();
 524       if(pt<fPtCut) continue;
 525       Int_t tower = findTower(particle->Phi(),particle->Eta());
 526       if (tower >= 0) {
 527         (*fTowers)[tower] += particle;
 528 #ifdef DOHISTOS
 529       // calculate Et for a massless particle...
 530       // idea for Et check: Et = E *sin(theta);
 531       Float_t Et = TMath::Sqrt(pt*pt);
 532         // just as check - fill particle Et in histogram
 533         ((TH2F *)fEventHistos.At(1))->Fill(particle->Eta(),particle->Phi(),Et);
 534 #endif
 535       }
 536     }
 537   }
 538   delete iter;
 539 }
 540
 541 #ifdef ALICEINTERFACE
 542 void AliTkConeJetFinderV2::fillTowersFromAliParticles(const TClonesArray *particles)
 543 {
 544   // input TClonesArray with AliParticles
 545   // fills the Et grid from these particles
 546   if (!particles) {
 547     return;
 548   }
 549
 550   if(fAliParticles) fAliParticles->Clear();
 551   else fAliParticles=new TClonesArray("TParticle",0);
 552   fAliParticles->Expand(particles->GetEntriesFast());
 553
 554   // loop over all particles...
 555   AliJetParticle *aliparticle = NULL;
 556   TIterator *iter = particles->MakeIterator();
 557
 558   Int_t i=0;
 559   //Float_t meanet=0,totet=0;
 560   while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
 561     Float_t pt=aliparticle->Pt();
 562     if(pt<fPtCut) continue;
 563     // particle is accepted through reader in JETAN
 564     TParticle *particle = new((*fAliParticles)[i]) TParticle(0,0,0,0,0,0,aliparticle->Px(),aliparticle->Py(),
 565                                                    aliparticle->Pz(),aliparticle->Energy(),0,0,0,0);
 566     //mark particle (-123 for Pythia)
 567     particle->SetWeight(aliparticle->GetType());
 568     Int_t tower = findTower(particle->Phi(),particle->Eta());
 569     if (tower >= 0) {
 570       (*fTowers)[tower] += particle;
 571       //meanet+=pt;totet+=pt;
 572 #ifdef DOHISTOS
 573       // calculate Et for a massless particle...
 574       // idea for Et check: Et = E *sin(theta);
 575       Float_t Et = TMath::Sqrt(particle->Pt()*particle->Pt());
 576       // just as check - fill particle Et in histogram
 577       ((TH2F *)fEventHistos.At(1))->Fill(particle->Eta(),particle->Phi(),Et);
 578 #endif
 579     }
 580     i++;
 581   }
 582   delete iter;
 583   //cout << "Mean particle " << meanet/i << " " << totet << endl;
 584 }
 585
 586 #endif
 587
 588 Int_t AliTkConeJetFinderV2::findTower(Float_t phi,Float_t eta)
 589 {
 590   if ((phi < fPhiMin) || (phi > fPhiMax) ||
 591       (eta < fEtaMin) || (eta > fEtaMax)) {
 592     return -1;
 593   }
 594   Int_t phibins = (Int_t) ((phi - fPhiMin) / fPhiWidth);
 595   Int_t etabins = (Int_t) ((eta - fEtaMin) / fEtaWidth);
 596   return (etabins * fPhiBins + phibins);
 597 }
 598
 599 void AliTkConeJetFinderV2::createSeedPoints()
 600 {
 601   // function should decide if it makes sense
 602   // to use fEtCut for seed towers
 603   // + midpoints or simply tower list
 604   // uses only tower list so far...
 605
 606 #if 0
 607   Float_t met=0.;Float_t set=0.;
 608   Int_t counter=0;
 609   for(vector<tower>::iterator iter = fTowers->begin();
 610       iter != fTowers->end(); ++iter) {
 611     met+=iter->getEt();
 612     set+=iter->getEt()*iter->getEt();
 613     counter++;
 614   }
 615   if(counter>1){
 616     met/=counter;
 617     set=set/counter-met*met;
 618     if(set>0) set=TMath::Sqrt(set)/(counter-1);
 619     else set=0;
 620     if(fOutput) cout << "Tower Mean Et: " << met << "+-" << set << " " << counter << endl;
 621   }
 622   Float_t meanEt=met-set; //store mean tower et
 623
 624   for(vector<tower>::iterator pos = fTowers->begin();
 625       pos != fTowers->end(); ++pos) {
 626     if(pos->getEt()<meanEt){
 627       pos->clear();
 628     }
 629   }
 630 #endif
 631
 632   for(vector<tower>::iterator iter = fTowers->begin();
 633       iter != fTowers->end(); ++iter) {
 634     if(iter->getEt()>fEtCut){
 635       AliTkEtaPhiVector seedPoint((*iter).getEta(),(*iter).getPhi());
 636       fSeedPointsNew.push_back(seedPoint);
 637     }
 638   }
 639   if(fOutput) cout << "created " << fSeedPointsNew.size() << " seed points" << endl;
 640 }
 641
 642 void AliTkConeJetFinderV2::findProtojets()
 643 {
 644   const Float_t radiusSq = fRadius*fRadius;
 645   const Int_t maxIterations = 100;
 646   const Bool_t geoCut = kTRUE;
 647   const Float_t geoCutRadiusSq = 2*(fPhiWidth*fPhiWidth + fEtaWidth*fEtaWidth);
 648
 649   // loop over all seedpoints
 650   for(list<AliTkEtaPhiVector>::iterator iter = fSeedPointsNew.begin();
 651       iter != fSeedPointsNew.end(); ++iter) {
 652     // create a new protojet at seedpoint position
 653     protojet *pj = new protojet();
 654     pj->setCentroidPosition(*iter);
 655     // loop over all towers and add all within "radius" to the protojet...
 656     for (vector<tower>::iterator tower = fTowers->begin();
 657          tower != fTowers->end(); ++tower) {
 658       if(tower->getEt()<=0) continue;
 659       AliTkEtaPhiVector TwCenter((*tower).getEta(),(*tower).getPhi());
 660       AliTkEtaPhiVector protojetCenter(pj->getCentroidPosition());
 661       if (TwCenter.diffSq(protojetCenter) < radiusSq) {
 662         pj->addTower(&(*tower));
 663       }
 664     }
 665     //update mean values
 666     pj->update();
 667     // have added all towers within "radius" to protojet...
 668     // iterate seedpoint until stable
 669     // lets get the pT-weigthed center of the protojet...
 670     AliTkEtaPhiVector centerGeo = pj->getCentroidPosition();
 671     AliTkEtaPhiVector centerEt = pj->getEtWeightedPosition();
 672     Int_t iteration = 0;
 673     Bool_t wasGeoCut = kFALSE;
 674     while ((centerGeo.diffSq(centerEt) > __epsilon__) &&
 675            (iteration < maxIterations)) {
 676       iteration++;
 677       pj->eraseTowers();
 678       pj->setCentroidPosition(centerEt);
 679       centerGeo = pj->getCentroidPosition();
 680
 681       // if geoCut == kTRUE, break if it leaves original bin...
 682       if ((geoCut == kTRUE) &&
 683           ((*iter).diffSq(centerGeo) > geoCutRadiusSq)) {
 684         wasGeoCut = kTRUE;
 685         break;
 686       }
 687       // loop over all towers and add all within "radius" to the protojet...
 688       for (vector<tower>::iterator tower = fTowers->begin();
 689            tower != fTowers->end(); ++tower) {
 690         if(tower->getEt()<=0) continue;
 691         AliTkEtaPhiVector TwCenter((*tower).getEta(),(*tower).getPhi());
 692         if (TwCenter.diffSq(centerEt) < radiusSq) {
 693           pj->addTower(&(*tower));
 694         }
 695       }
 696       // have added all towers within "radius" to protojet...
 697       centerEt = pj->getEtWeightedPosition();
 698     }
 699
 700     // we have a stable protojet (or run out of iterations...)
 701 #if 0
 702     if((fOutput) &&(iteration==maxIterations))
 703       cout << " FindProtojets warning: max iterations of "
 704       << maxIterations << " reached !!!" << endl;
 705 #endif
 706     // let's add it to the protojet list...
 707     if (!wasGeoCut) {
 708       addProtojet(pj);
 709     } else {
 710       delete pj;
 711     }
 712   }
 713   //if(fOutput) cout << "found " << fProtojets.size() << "proto jets" << endl;
 714 }
 715
 716 void AliTkConeJetFinderV2::addProtojet(protojet *pj)
 717 {
 718   if (!pj) {
 719     return;
 720   }
 721   for(list<protojet*>::const_iterator iter = fProtojets.begin();
 722       iter != fProtojets.end(); ++iter) {
 723     if ((*pj) == (*(*iter))) {
 724       delete pj;
 725       return;
 726     }
 727   }
 728   fProtojets.push_back(pj);
 729 }
 730
 731 void AliTkConeJetFinderV2::dumpProtojets(Float_t etmin)
 732 {
 733   for(list<protojet*>::const_iterator iter = fProtojets.begin();
 734       iter != fProtojets.end(); ++iter) {
 735     if ((*iter)->getEt() > etmin) {
 736       cout << (*(*iter)) << endl;
 737     }
 738   }
 739 }
 740
 741 void AliTkConeJetFinderV2::findJets()
 742 {
 743   // loop over all protojets until list is empty
 744   while(!fProtojets.empty()) {
 745     // find the protojet with maximum Et
 746     // should be easy to have a sorted list...
 747     list<protojet*>::iterator maxEtProtojet = fProtojets.begin();
 748     Float_t maxEt = 0;
 749     for (list<protojet*>::iterator iter = fProtojets.begin();
 750          iter != fProtojets.end(); ++iter) {
 751       if ((*iter)->getEt() > maxEt) {
 752         maxEt = (*iter)->getEt();
 753         maxEtProtojet = iter;
 754       }
 755     }
 756     // we've found the protojet with the highest Et - remove it from the list
 757     protojet *jet1 = *maxEtProtojet;
 758     fProtojets.erase(maxEtProtojet);
 759     // loop again over all protojets to find sharing jet with highest Et
 760     list<protojet*>::iterator maxEtNeighbor = fProtojets.begin();
 761     maxEt = 0;
 762     for (list<protojet*>::iterator iter = fProtojets.begin();
 763          iter != fProtojets.end(); ++iter) {
 764       if (((*iter)->getEt() > maxEt) &&
 765           (jet1->shareTowers(*iter))) {
 766         maxEt = (*iter)->getEt();
 767         maxEtNeighbor = iter;
 768       }
 769     }
 770     if (maxEt>0) {
 771       // jet's share towers
 772       // merging splitting step...
 773       protojet *jet2 = (*maxEtNeighbor);
 774       fProtojets.erase(maxEtNeighbor);
 775       splitMergeJets(jet1,jet2);
 776     } else {
 777       // protojet 1 doesn't share towers with other protojets, make it a jet...
 778       addJet(jet1);
 779     }
 780   }
 781   if(fOutput){
 782     cout << "found " << fJets.size() << " jets" << endl;
 783     dumpJets();
 784   }
 785 }
 786
 787 void AliTkConeJetFinderV2::splitMergeJets(protojet *jet1,protojet *jet2)
 788 {
 789   const Float_t EtRatioCut = 0.5;
 790
 791   // let's calcualte the shared energy...
 792   Float_t fEtShared = 0.0;
 793   list<tower*> sharedTowers = jet1->getSharedTowerList(jet2);
 794   for (list<tower*>::const_iterator iter = sharedTowers.begin();
 795        iter != sharedTowers.end(); ++iter) {
 796     fEtShared += (*iter)->getEt();
 797   }
 798   Float_t fEtJet2 = jet2->getEt();
 799   // calculate the ratio of the shared energy and decided if split or merge
 800   Float_t fEtRatio = EtRatioCut + 1.;
 801   if(fEtJet2) fEtRatio = fEtShared / fEtJet2;
 802   if (fEtRatio > EtRatioCut) {
 803     mergeJets(jet1,jet2);
 804   } else {
 805     splitJets(jet1,jet2);
 806   }
 807 }
 808
 809 void AliTkConeJetFinderV2::mergeJets(protojet *jet1,protojet *jet2) {
 810   // merge protojets...
 811   list<tower*> jet2Towers = jet2->getTowerList();
 812   for (list<tower*>::const_iterator iter = jet2Towers.begin();
 813        iter != jet2Towers.end(); ++iter) {
 814     // add towers from protojet2 to protojet1...
 815     // don't add shared towers...
 816     if (!jet1->hasTower((*iter))) {
 817       jet1->addTower((*iter));
 818     }
 819   }
 820   addProtojet(jet1);
 821   delete jet2;
 822 }
 823
 824 void AliTkConeJetFinderV2::splitJets(protojet *jet1,protojet *jet2)
 825 {
 826   // split protojets...
 827   list<tower*> sharedTowers = jet1->getSharedTowerList(jet2);
 828   for (list<tower*>::const_iterator iter = sharedTowers.begin();
 829        iter != sharedTowers.end(); ++iter) {
 830     // search nearest jet
 831     if (jet1->diffToCenter(*iter) < jet2->diffToCenter(*iter)) {
 832       // tower closer to jet1
 833       jet2->eraseTower(*iter);
 834     } else {
 835       // tower closer to jet2
 836       jet1->eraseTower(*iter);
 837     }
 838   }
 839   if (jet1->shareTowers(jet2)) {
 840     cerr << "!!! SplitJets: Something is wrong !!!" << endl;
 841     addProtojet(jet1);
 842     delete jet2;
 843     return;
 844   }
 845   addProtojet(jet1);
 846   addProtojet(jet2);
 847 }
 848
 849 void AliTkConeJetFinderV2::dumpJets()
 850 {
 851   cout << "----- found jets > " << fEtMinJet << " GeV -----" << endl;
 852   for (list<protojet*>::const_iterator iter = fJets.begin();
 853        iter != fJets.end(); ++iter) {
 854     if((*iter)->getEt()>fEtMinJet) cout << (*(*iter)) << endl;
 855   }
 856 }
 857
 858 #ifdef DOHISTOS
 859 void AliTkConeJetFinderV2::createHistos()
 860 {
 861 }
 862 #endif
 863
 864 #ifdef DOHISTOS
 865 void AliTkConeJetFinderV2::createEventHistos()
 866 {
 867   fHistFile = new TFile("$JF_DATADIR/ConeFinderV2.root","RECREATE");
 868   TH2F *h = new TH2F("Etbins","Etbins",
 869                      fEtaBins,fEtaMin,fEtaMax,
 870                      fPhiBins,fPhiMin,fPhiMax);
 871   h->GetXaxis()->SetTitle("#eta");
 872   h->GetYaxis()->SetTitle("#phi");
 873   h->GetZaxis()->SetTitle("E_{t} (GeV)");
 874   fEventHistos.Add(h);
 875
 876   h = new TH2F("Etbin_check","Etbin_check",
 877                fEtaBins,fEtaMin,fEtaMax,
 878                fPhiBins,fPhiMin,fPhiMax);
 879   h->GetXaxis()->SetTitle("#eta");
 880   h->GetYaxis()->SetTitle("#phi");
 881   h->GetZaxis()->SetTitle("E_{t} (GeV)");
 882   fEventHistos.Add(h);
 883
 884   h = new TH2F("seedpoints","seedpoints",
 885                fEtaBins,fEtaMin,fEtaMax,
 886                fPhiBins,fPhiMin,fPhiMax);
 887   h->GetXaxis()->SetTitle("#eta");
 888   h->GetYaxis()->SetTitle("#phi");
 889   fEventHistos.Add(h);
 890 }
 891 #endif
 892 #ifdef DOHISTOS
 893 void AliTkConeJetFinderV2::clearEventHistos()
 894 {
 895   ((TH2F *)fEventHistos.At(0))->Reset();
 896   ((TH2F *)fEventHistos.At(1))->Reset();
 897   ((TH2F *)fEventHistos.At(2))->Reset();
 898 }
 899 #endif
 900
 901 #ifdef DOHISTOS
 902 void AliTkConeJetFinderV2::writeEventHistos()
 903 {
 904   fHistFile->cd();
 905   ((TH2F *)fEventHistos.At(0))->Write();
 906   ((TH2F *)fEventHistos.At(1))->Write();
 907   ((TH2F *)fEventHistos.At(2))->Write();
 908 }
 909 #endif
 910
 911 #ifdef DOHISTOS
 912 void AliTkConeJetFinderV2::writeHistos()
 913 {
 914 }
 915 #endif
 916
 917 #ifdef DOHISTOS
 918 void AliTkConeJetFinderV2::fillTowerHist()
 919 {
 920   for (vector<tower>::iterator i=(*fTowers).begin(); i!=(*fTowers).end(); ++i) {
 921     tower t = *i;
 922     Float_t eta = (t.getEtaMax() + t.getEtaMin()) / 2.;
 923     Float_t phi = (t.getPhiMax() + t.getPhiMin()) / 2.;
 924     ((TH2F *)fEventHistos.At(0))->Fill(eta,phi,t.getEt());
 925   }
 926 }
 927 #endif
 928
 929 void AliTkConeJetFinderV2::setEvOutFilename(const Char_t *filename)
 930 {
 931   if (!fEvout_name) {
 932     fEvout_name = new Char_t[4096];
 933   }
 934
 935   strcpy(fEvout_name,filename);
 936 }
 937