JETAN/jetan2004/AliTkChargedJetFinder.cxx

   1 // $Id$
   2 //--------------------------------------------------------------------------
   3 // AliTkChargedJetFinder.cxx
   4 // implementation of a simple jet finder for charged particles
   5 // based on CDF PRD 65, 092002 (2002)
   6 // T. Kollegger <kollegge@ikf.physik.uni-frankfurt.de> 10/12/2002
   7 //--------------------------------------------------------------------------
   8
   9 #include <Riostream.h>
  10 //-----------------------------------------------------------------------
  11 // STL includes
  12 #include <list>
  13 #include <map>
  14 //-----------------------------------------------------------------------
  15 // ROOT includes
  16 #include <TROOT.h>
  17 #include <TFile.h>
  18 #include <TTree.h>
  19 #include <TH1.h>
  20 #include <TVector2.h>
  21 #include <TClonesArray.h>
  22 #include <TParticle.h>
  23 //--------------------------------------------------------------------------
  24 #include "AliTkEtaPhiVector.h"
  25 #ifdef DOCHAERGED
  26 #include "TkChargedJet.h"
  27 #else
  28 #include "AliTkConeJet.h"
  29 #include "AliTkConeJetEvent.h"
  30 #include "AliTkTowerV2.h"
  31 #endif
  32 //--------------------------------------------------------------------------
  33 #ifdef ALICEINTERFACE
  34 #include <AliJetParticle.h>
  35 #include <AliJetEvent.h>
  36 #include <AliJetEventParticles.h>
  37 #endif
  38 //-------------------------------------------------------------------------
  39 #include "AliTkChargedJetFinder.h"
  40
  41 //------------------------------------------------------------------------
  42 // implementation of the jet helper class
  43 //------------------------------------------------------------------------
  44
  45 jet::jet(){
  46   fPt=-999;
  47   fEta=-999;
  48   fPhi=-999;
  49   fNParticles=0;
  50 }
  51
  52 AliTkEtaPhiVector jet::getCentroid() const
  53 {
  54   return AliTkEtaPhiVector(fEta,fPhi);
  55 }
  56
  57 void jet::addParticle(TParticle *particle)
  58 {
  59   // add a particle to the jet
  60   fParticles.push_back(particle);
  61   if(fNParticles==0) {
  62     fEta=particle->Eta();
  63     fPhi=particle->Phi();
  64     fPt=0;
  65   }
  66   fPt+=particle->Pt();
  67   fNParticles++;
  68 }
  69
  70 Double_t jet::getPt() const
  71 {
  72 #if 1
  73   return fPt;
  74 #else
  75   // returns the pt of the jet
  76   Double_t pt = 0.0;
  77
  78   // loop over all particles which are assigned to the jet...
  79   // and sum the pt of them
  80   for(list<TParticle *>::const_iterator i = particles.begin();
  81       i != particles.end(); ++i) {
  82     TParticle *particle = *i;
  83     pt += particle->Pt();
  84   }
  85
  86   return pt;
  87 #endif
  88 }
  89
  90 Int_t jet::getNParticles() const
  91 {
  92   return fNParticles;
  93 }
  94
  95 TParticle *jet::getParticle(Int_t i)  const
  96 {
  97   if ((i < 0) || ((UInt_t)i > fParticles.size()-1)) {
  98     cerr << "jet:: out of range" << i << endl;
  99     return NULL;
 100   }
 101   list<TParticle *>::const_iterator iter = fParticles.begin();
 102   for(Int_t pos = 0; pos < i+1; pos++) {
 103     ++iter;
 104   }
 105   return (*iter);
 106 }
 107
 108 TClonesArray *jet::getParticles() const
 109 {
 110   TClonesArray *parts = new TClonesArray("TParticle",this->getNParticles());
 111   Int_t i = 0;
 112   for(list<TParticle *>::const_iterator iter = fParticles.begin();
 113       iter != fParticles.end(); ++iter, i++) {
 114     TParticle *particle = *iter;
 115     if (particle == NULL) {
 116       cerr << "jet: What's wrong? NULL pointer in particles" << endl;
 117       continue;
 118     }
 119     new ((*parts)[i]) TParticle(*particle);
 120   }
 121   return parts;
 122 }
 123
 124 Double_t jet::getDiff(TParticle *particle) const
 125 {
 126   // calculate the difference between jet center and particle in eta-phi space
 127   return sqrt(getDiffSq(particle));
 128 }
 129
 130 Double_t jet::getDiffSq(TParticle *particle) const
 131 {
 132   // calculate the square of the difference between jet center
 133   // and particle in eta-phi space
 134   AliTkEtaPhiVector jetCenter = getCentroid();
 135   AliTkEtaPhiVector particlePos(particle->Eta(),particle->Phi());
 136   return jetCenter.diffSq(particlePos);
 137 }
 138
 139 ostream& operator<< (ostream& s,jet& j)
 140 {
 141   return s << "jet: position "  << j.getCentroid()
 142            << " #particles="    << j.getNParticles()
 143            << " pT="            << j.getPt();
 144 }
 145
 146
 147 //------------------------------------------------------------------------
 148 // implementation of the jet finder class
 149 //------------------------------------------------------------------------
 150
 151 ClassImp(AliTkChargedJetFinder)
 152
 153 AliTkChargedJetFinder::AliTkChargedJetFinder() : TObject()
 154 {
 155   fOutput=0;
 156   fR=0;
 157   fRSq=0;
 158   fEtaMin=0;
 159   fEtaMax=0;
 160   fPtCut=0;
 161   fPtSeed=0;
 162   fMinJetPt=0;
 163
 164 #ifdef DOCHARGED
 165   fMyTJet=0;
 166 #else
 167   fEvoutevent=0;
 168 #endif
 169
 170 #ifdef DOHISTOS
 171   // histograms...
 172   fOutput_name=0;
 173   fHistos=0;
 174   fHistoutfile=0;
 175 #endif
 176
 177   fEvout_name=0;
 178   fEvoutfile=0;
 179   fEvouttree=0;
 180
 181 #ifdef ALICEINTERFACE
 182   fAliParticles=0;
 183 #endif
 184
 185   defaultSettings();
 186 }
 187
 188 AliTkChargedJetFinder::~AliTkChargedJetFinder()
 189 {
 190   if(fEvoutevent) delete fEvoutevent;
 191   if(fEvoutfile) delete fEvoutfile;
 192   if(fEvout_name) delete[] fEvout_name;
 193 #ifdef ALICEINTERFACE
 194   if(fAliParticles) delete fAliParticles;
 195 #endif
 196 #ifdef DOHISTOS
 197   if(fHistos) delete fHistos;
 198   if(fHistoutFile) delete fHistoutFile;
 199   if(fOutput_name) delete[] fOutput_name;
 200 #endif
 201 }
 202
 203 void AliTkChargedJetFinder::defaultSettings()
 204 {
 205   // default settings for a first study...
 206   // jet finder radius
 207   setFinderR(0.7);
 208   // lower bound of eta range
 209   setEtaMin(-1);
 210   // high bound of eta range
 211   setEtaMax(1);
 212   // seed value
 213   setPtSeed(5);
 214   // minimum jet value
 215   setMinJetPt(5);
 216
 217   fEvout_name=0;
 218   setEvOutFilename("$JF_DATADIR/charged_jets.evout.root");
 219 #ifdef DOHISTOS
 220   fOutput_name=0;
 221   setHistFilename("$JF_DATADIR/charged_jets.hist.root");
 222 #endif
 223 }
 224
 225 void AliTkChargedJetFinder::init()
 226 {
 227 #ifdef DOHISTOS
 228   // initalization of the jet finder
 229   fHistos = new TObjArray(10);
 230   TH1F *hist = new TH1F("ncharged","charged tracks",31,-0.5,30.5);
 231   hist->GetXaxis()->SetTitle("n_{charged} tracks");
 232   hist->GetYaxis()->SetTitle("entries");
 233   fHistos->Add(hist);
 234   hist = new TH1F("jet_pt","jet p_{T}",100,0,100);
 235   hist->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 236   hist->GetYaxis()->SetTitle("entries");
 237   fHistos->Add(hist);
 238
 239   // declare all histograms
 240   if (fOutput_name) {
 241     if(fOutput) cout << "Writing summary output to " << getHistFilename() << endl;
 242     fHistoutfile = new TFile(getHistFilename(),"RECREATE");
 243   }
 244 #endif
 245
 246 #ifdef DOCHARGED
 247   if (fEvout_name) {
 248     if(fOutput) cout << "Writing event output to " << getEvOutFilename() << endl;
 249     fEvoutfile = new TFile(getEvOutFilename(),"RECREATE");
 250     fMyTJet = new TkChargedJet();
 251     fEvouttree = new TTree("jets","");
 252     fEvouttree->Branch("jet","TkChargedJet",&(this->fMyTJet));
 253   }
 254 #else
 255   fEvoutfile = 0;
 256   fEvouttree = 0;
 257   fEvoutevent = new AliTkConeJetEvent();
 258
 259   if (getEvOutFilename()) {
 260     if(fOutput) cout << "Writing finder output to " << getEvOutFilename() << endl;
 261     fEvoutfile = new TFile(getEvOutFilename(),"RECREATE");
 262     fEvouttree = new TTree("jets","TKConeJetFinderV2 jets");
 263     fEvouttree->Branch("ConeFinder","AliTkConeJetEvent",&(this->fEvoutevent),32000,0);
 264   }
 265 #endif
 266 }
 267
 268 void AliTkChargedJetFinder::initEvent(TClonesArray *newParticles,Int_t type)
 269 {
 270   // initalizes the jet finder for a new event
 271 #ifndef DOCHARGED
 272   if(fEvoutevent){
 273     fEvoutevent->Clear();
 274     fEvoutevent->setRadius(fR);
 275     fEvoutevent->setEtCut(fPtSeed);
 276     fEvoutevent->setPtCut(fPtCut);
 277     fEvoutevent->setDesc(*new TString("AliTkChargedJetFinder"));
 278   }
 279 #endif
 280
 281   // delete old particle list and old jets
 282   fParticles.erase(fParticles.begin(),fParticles.end());
 283   fJets.erase(fJets.begin(),fJets.end());
 284
 285   TParticle *particle;
 286   TIterator *iter = newParticles->MakeIterator();
 287 #ifdef testptr
 288   list<TParticle *>::const_iterator siter;
 289 #endif
 290
 291   // create new particle list
 292   switch (type) {
 293   case 1:
 294     // TParticles in TClonesArray
 295     while ((particle = (TParticle *) iter->Next()) != NULL) {
 296       if (isTParticleAccepted(particle)) {
 297         addParticle(particle);
 298       }
 299     }
 300     if(fOutput) cout << "found particles for jets: " << fParticles.size() << endl;
 301     break;
 302   case 2:
 303     cerr << "AliTkChargedJetFinder: don't know how to fill Et hist from TClonesArray with type " << type << endl;
 304     ;
 305     break;
 306   }
 307 #ifdef testptr
 308     for (siter=particles.begin(); siter != particles.end(); ++siter) {
 309       cout << "particle eta=" << (*siter)->Eta()
 310            << " phi=" << (*siter)->Phi()
 311            << " pT=" << (*siter)->Pt()
 312            << " pointer=" << (*siter)
 313            << endl;
 314     }
 315 #endif
 316 }
 317
 318 #ifdef ALICEINTERFACE
 319 void AliTkChargedJetFinder::initEvent(const AliJetEventParticles *p,TString desc)
 320 {
 321   if(fEvoutevent){
 322     fEvoutevent->Clear();
 323     fEvoutevent->setRadius(fR);
 324     fEvoutevent->setEtCut(fPtSeed);
 325     fEvoutevent->setPtCut(fPtCut);
 326     fEvoutevent->setDesc(desc);
 327     //fEvoutevent->setJetParticles(p);
 328   }
 329
 330   // delete old particle list and old jets
 331   fParticles.erase(fParticles.begin(),fParticles.end());
 332   fJets.erase(fJets.begin(),fJets.end());
 333
 334   const TClonesArray *particles=p->GetParticles();
 335   if(fAliParticles) fAliParticles->Clear();
 336   else fAliParticles=new TClonesArray("TParticle",0);
 337   fAliParticles->Expand(particles->GetEntriesFast());
 338
 339   // loop over all particles...
 340   Int_t i=0;
 341   AliJetParticle *aliparticle = NULL;
 342   TIterator *iter = particles->MakeIterator();
 343   while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
 344     Float_t pt=aliparticle->Pt();
 345     if(pt<fPtCut) continue;
 346     // particle is accepted through reader in JETAN
 347     TParticle *particle = new((*fAliParticles)[i]) TParticle(0,0,0,0,0,0,aliparticle->Px(),aliparticle->Py(),
 348                                                    aliparticle->Pz(),aliparticle->Energy(),0,0,0,0);
 349     //mark particle (-123 for Pythia)
 350     particle->SetWeight(aliparticle->GetType());
 351     //particle is accepted through reader
 352     addParticle(particle);
 353     i++;
 354   }
 355   delete iter;
 356 }
 357 #endif
 358
 359 void AliTkChargedJetFinder::run()
 360 {
 361   // loop over particles as long as there are
 362   // some which are not assigned to a jet
 363   if (fParticles.empty()) {
 364     return;
 365   }
 366
 367   list<TParticle *>::iterator iter;
 368   while (!fParticles.empty()) {
 369     // get particle with highest pT...
 370     iter = fParticles.begin();
 371     if((*iter)->Pt()<fPtSeed) {
 372       fParticles.erase(iter);
 373       continue;
 374     }
 375     // create a new jet...
 376     jet myJet;
 377     // and add the particle with highest pT...
 378     myJet.addParticle(*iter);
 379     // delete it from the list of unassigned particles...
 380     fParticles.erase(iter);
 381     // loop over all particles to see if they are within the new jet...
 382     for(iter = fParticles.begin(); iter != fParticles.end(); ++iter) {
 383       if (myJet.getDiffSq(*iter) < fRSq) {
 384         // add this particle to the jet...
 385         myJet.addParticle(*iter);
 386         // remove this particle from the list of unassigned ones...
 387         fParticles.erase(iter);
 388         // and start again from the beginnig of the list
 389         // -- the pT-weigthed center of the jet might have changed,
 390         //    so one might add now a particle which was before rejected
 391         iter = fParticles.begin();
 392         iter--;
 393       }
 394     }
 395     // we have our final jet - add it to the jet list...
 396     fJets.push_back(myJet);
 397   }
 398 }
 399
 400 void AliTkChargedJetFinder::finishEvent()
 401 {
 402   // analyse event...
 403   checkJets();
 404
 405 #ifdef DOCHARGED
 406   list<jet>::iterator highjet = findHighestJet();
 407   if (highjet == NULL) {
 408     cerr << "AliTkChargedJetFinder:: no jet found - no output! " << endl;
 409     return;
 410   }
 411   if(fOutput) {
 412     cout << "high jet ncharged=" << (*highjet).getNParticles()
 413          << " pt="               << (*highjet).getPt() << endl;
 414     cout << *highjet << endl;
 415   }
 416
 417 #ifdef DOHISTOS
 418   // write out event histograms
 419   if ((*highjet).getPt() > 30) {
 420     ((TH1F *)fHistos->At(0))->Fill((*highjet).getNParticles());
 421   }
 422   ((TH1F *)fHistos->At(1))->Fill((*highjet).getPt());
 423 #endif
 424
 425   TkChargedJet *jet = new TkChargedJet(*highjet);
 426   if ((fEvoutfile) && (fEvoutfile->IsOpen())) {
 427     delete fMyTJet;
 428     if (jet) {
 429       fMyTJet = jet;
 430     } else {
 431       fMyTJet = new TkChargedJet();
 432     }
 433     fEvoutfile->cd();
 434     fEvouttree->Fill();
 435   }
 436 #else /******************************************************/
 437   if (fEvoutevent){
 438     AliTkConeJet *conejet = new AliTkConeJet();
 439     AliTkTowerV2 *myTower = new AliTkTowerV2();
 440     list<jet>::iterator iter;
 441     for (iter = fJets.begin();
 442          iter != fJets.end(); ++iter) {
 443
 444       Float_t jetet=(*iter).getPt();
 445       if(jetet<fMinJetPt) continue;
 446       conejet->Clear();
 447       conejet->setEta((*iter).Eta());
 448       conejet->setPhi((*iter).Phi());
 449       conejet->setEt(jetet);
 450
 451       myTower->Clear();
 452       myTower->setEta((*iter).Eta());
 453       myTower->setPhi((*iter).Phi());
 454       myTower->setEt(jetet);
 455       myTower->setParticleList((*iter).getParticles());
 456       conejet->addTower(myTower);
 457       conejet->calculateValues();
 458       fEvoutevent->addJet(conejet);
 459     }
 460     fEvoutevent->sortJets();
 461     fEvoutevent->Print("");
 462     if (fEvouttree)
 463       fEvouttree->Fill();
 464     delete conejet;
 465     delete myTower;
 466   }
 467 #endif
 468 }
 469
 470 void AliTkChargedJetFinder::finish()
 471 {
 472 #ifdef DOHISTOS
 473   // write out histograms
 474   if ((fHistoutfile) && (fHistoutfile->IsOpen())) {
 475     fHistoutfile->cd();
 476     ((TH1F *)fHistos->At(0))->Write();
 477     ((TH1F *)fHistos->At(1))->Write();
 478     fHistoutfile->Write();
 479     fHistoutfile->Close();
 480   }
 481 #endif
 482   if ((fEvoutfile) && (fEvoutfile->IsOpen())) {
 483     fEvoutfile->cd();
 484     fEvouttree->Write();
 485 #ifdef DOHISTOS
 486     ((TH1F *)fHistos->At(0))->Write();
 487 #endif
 488     fEvoutfile->Write();
 489     fEvoutfile->Close();
 490   }
 491 }
 492
 493 void AliTkChargedJetFinder::setEvOutFilename(const Char_t *filename)
 494 {
 495   if (!fEvout_name) {
 496     fEvout_name = new Char_t[4096];
 497   }
 498   strcpy(fEvout_name,filename);
 499 }
 500
 501 #ifdef DOHISTOS
 502 void AliTkChargedJetFinder::setHistFilename(const Char_t *filename)
 503 {
 504   if (!fOutput_name) {
 505     fOutput_name = new Char_t[4096];
 506   }
 507   strcpy(fOutput_name,filename);
 508 }
 509 #endif
 510
 511 bool AliTkChargedJetFinder::isTParticleAccepted(TParticle *particle)
 512 {
 513   // particle cuts...
 514
 515   // check if particle is stable
 516   if (particle->GetStatusCode() != 1) return false;
 517
 518   // check for charge
 519   TParticlePDG *partPDG;
 520   partPDG = particle->GetPDG();
 521   if (partPDG->Charge() == 0) {
 522     return false;
 523   }
 524
 525   // check eta range
 526   if ((particle->Eta() < fEtaMin) || (particle->Eta() > fEtaMax)) {
 527     return false;
 528   }
 529
 530   return true;
 531 }
 532
 533 void AliTkChargedJetFinder::addParticle(TParticle *particle)
 534 {
 535   // add particle at right position
 536   // first point to speed up program, use a sorted tree for example...
 537
 538   // if first particle, store at beginning
 539   if (fParticles.empty()) {
 540     fParticles.push_front(particle);
 541     return;
 542   }
 543
 544   // search insert position
 545   list<TParticle *>::iterator iter;
 546   for (iter = fParticles.begin(); iter != fParticles.end(); ++iter) {
 547     if ((*iter)->Pt() < particle->Pt()) {
 548       fParticles.insert(iter,particle);
 549       break;
 550     }
 551   }
 552 }
 553
 554 list<jet>::iterator AliTkChargedJetFinder::findHighestJet()
 555 {
 556   list<jet>::iterator iter;
 557   list<jet>::iterator high;
 558
 559   if (fJets.empty()) {
 560     return fJets.begin();
 561   }
 562
 563   high = fJets.begin();
 564   for (iter = fJets.begin(); iter != fJets.end(); ++iter) {
 565     //cout << *iter << endl;
 566     if ((*iter).getPt() > (*high).getPt()) {
 567       high = iter;
 568     }
 569   }
 570   return high;
 571 }
 572
 573 void AliTkChargedJetFinder::checkJets()
 574 {
 575   list<jet>::iterator i;
 576   list<jet>::iterator j;
 577   for (i = fJets.begin(); i != fJets.end(); ++i) {
 578     for (j = i, ++j; j != fJets.end(); ++j) {
 579       AliTkEtaPhiVector v1 = (*i).getCentroid();
 580       AliTkEtaPhiVector v2 = (*j).getCentroid();
 581       Double_t diff = v1.diff(v2);
 582       if (diff < fR) {
 583         cerr << "AliTkChargedJetFinder: Something is wrong - jets to close..." << endl;
 584         cerr << v1 << " pt=" << (*i).getPt() <<  " -- "
 585              << v2 << " pt=" << (*j).getPt() << endl;
 586       }
 587     }
 588   }
 589 }