EMCAL/AliEMCALJetFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*-- Authors: Andreas Morsch   (CERN)
  19 //*            J.L. Klay        (LBL)
  20 //*            Aleksei Pavlinov (WSU)
  21 //--
  22 //--
  23 //
  24
  25 #include <stdio.h>
  26 // From root ...
  27 #include <TArrayF.h>
  28 #include <TAxis.h>
  29 #include <TBranchElement.h>
  30 #include <TCanvas.h>
  31 #include <TClonesArray.h>
  32 #include <TFile.h>
  33 #include <TH1.h>
  34 #include <TH2.h>
  35 #include <TList.h>
  36 #include <TPDGCode.h>
  37 #include <TPad.h>
  38 #include <TParticle.h>
  39 #include <TParticlePDG.h>
  40 #include <TPaveText.h>
  41 #include <TPythia6Calls.h>
  42 #include <TROOT.h>
  43 #include <TStyle.h>
  44 #include <TTree.h>
  45 #include <TBrowser.h>
  46
  47 // From AliRoot ...
  48 #include "AliEMCAL.h"
  49 #include "AliEMCALDigit.h"
  50 #include "AliEMCALDigitizer.h"
  51 #include "AliEMCALFast.h"
  52 #include "AliEMCALGeometry.h"
  53 #include "AliEMCALHadronCorrection.h"
  54 #include "AliEMCALHit.h"
  55 #include "AliEMCALJetFinder.h"
  56 #include "AliEMCALJetMicroDst.h"
  57 #include "AliHeader.h"
  58 #include "AliMagF.h"
  59 #include "AliMagFCM.h"
  60 #include "AliRun.h"
  61 #include "AliGenerator.h"
  62 #include "AliEMCALGetter.h"
  63 // Interface to FORTRAN
  64 #include "Ecommon.h"
  65 #include "AliMC.h"
  66
  67
  68 ClassImp(AliEMCALJetFinder)
  69
  70 //____________________________________________________________________________
  71 AliEMCALJetFinder::AliEMCALJetFinder()
  72 {
  73 // Default constructor
  74     fJets             = 0;
  75     fNjets            = 0;
  76     fLego             = 0;
  77     fLegoB            = 0;
  78
  79     fTrackList        = 0;
  80     fPtT              = 0;
  81     fEtaT             = 0;
  82     fPhiT             = 0;
  83     fPdgT             = 0;
  84
  85     fTrackListB       = 0;
  86     fPtB              = 0;
  87     fEtaB             = 0;
  88     fPhiB             = 0;
  89     fPdgB             = 0;
  90
  91     fHCorrection      = 0;
  92     fHadronCorrector  = 0;
  93
  94     fWrite            = 0;
  95     fOutFileName      = 0;
  96     fOutFile          = 0;
  97     fInFile           = 0;
  98     fEvent            = 0;
  99
 100     fRandomBg         = 0;
 101
 102     SetParametersForBgSubtraction();
 103 }
 104
 105 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
 106     : TTask(name, title)
 107 {
 108 // Constructor
 109 // Title is used in method GetFileNameForParameters();
 110 //
 111     fJets  = new TClonesArray("AliEMCALJet",10000);
 112     fNjets = 0;
 113     for (Int_t i = 0; i < 30000; i++)
 114     {
 115         fEtCell[i]  = 0.;
 116         fEtaCell[i] = 0.;
 117         fPhiCell[i] = 0.;
 118     }
 119     fLego       = 0;
 120     fLegoB      = 0;
 121
 122     fTrackList  = 0;
 123     fPtT        = 0;
 124     fEtaT       = 0;
 125     fPhiT       = 0;
 126     fPdgT       = 0;
 127
 128     fTrackListB       = 0;
 129     fPtB        = 0;
 130     fEtaB       = 0;
 131     fPhiB       = 0;
 132     fPdgB       = 0;
 133
 134     fHCorrection      = 0;
 135     fHadronCorrector  = 0;
 136     fBackground       = 0;
 137     fWrite            = 0;
 138     fOutFileName      = 0;
 139     fOutFile          = 0;
 140     fInFile           = 0;
 141     fEvent            = 0;
 142 //
 143     SetPtCut();
 144     SetMomentumSmearing();
 145     SetEfficiencySim();
 146     SetDebug();
 147     SetHadronCorrection();
 148     SetIncludeK0andN();
 149
 150     fRandomBg         = 0;
 151     SetParametersForBgSubtraction();
 152 }
 153
 154 void AliEMCALJetFinder::SetParametersForBgSubtraction
 155 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
 156 {
 157 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
 158 // at WSU Linux cluster - 11-feb-2002
 159 // These parameters must be tuned more carefull !!!
 160   SetMode(mode);
 161   SetMinMove(minMove);
 162   SetMaxMove(maxMove);
 163   SetPrecBg(precBg);
 164 }
 165
 166 //____________________________________________________________________________
 167 AliEMCALJetFinder::~AliEMCALJetFinder()
 168 {
 169 // Destructor
 170 //
 171     if (fJets){
 172         fJets->Delete();
 173         delete fJets;
 174     }
 175     delete fLego;
 176     delete fLegoB;
 177     delete fhLegoTracks;
 178     delete fhLegoEMCAL;
 179     delete fhLegoHadrCorr;
 180     delete fhEff;
 181     delete fhCellEt;
 182     delete fhCellEMCALEt;
 183     delete fhTrackPt;
 184     delete fhTrackPtBcut;
 185     delete fhChPartMultInTpc;
 186
 187     delete[] fTrackList;
 188     delete[] fPtT;
 189     delete[] fEtaT;
 190     delete[] fPhiT;
 191     delete[] fPdgT;
 192
 193     delete[] fTrackListB;
 194     delete[] fPtB;
 195     delete[] fEtaB;
 196     delete[] fPhiB;
 197     delete[] fPdgB;
 198 }
 199
 200 #ifndef WIN32
 201 # define jet_finder_ua1 jet_finder_ua1_
 202 # define hf1 hf1_
 203 # define type_of_call
 204
 205 #else
 206 # define jet_finder_ua1 JET_FINDER_UA1
 207 # define hf1 HF1
 208 # define type_of_call _stdcall
 209 #endif
 210
 211 extern "C" void type_of_call
 212 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
 213                Float_t  etc[30000],  Float_t etac[30000],
 214                Float_t  phic[30000],
 215                Float_t& min_move, Float_t& max_move, Int_t& mode,
 216                Float_t& prec_bg,  Int_t& ierror);
 217
 218 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 219
 220
 221 void AliEMCALJetFinder::Init()
 222 {
 223 //
 224 // Geometry and I/O initialization
 225 //
 226 //
 227 //
 228 //  Get geometry parameters from EMCAL
 229 //
 230 //
 231 //  Geometry
 232   //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 233   //  AliEMCALGeometry* geom =
 234   //  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 235   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 236   AliEMCALGeometry* geom = gime->EMCALGeometry() ;
 237
 238 //    SetSamplingFraction(geom->GetSampling());
 239
 240     fNbinEta = geom->GetNZ();
 241     fNbinPhi = geom->GetNPhi();
 242     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 243     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 244     fEtaMin  = geom->GetArm1EtaMin();
 245     fEtaMax  = geom->GetArm1EtaMax();
 246     fDphi    = (fPhiMax-fPhiMin)/fNbinPhi;
 247     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
 248     fNtot    = fNbinPhi*fNbinEta;
 249     fWeightingMethod = kFALSE;
 250
 251 //
 252     SetCellSize(fDeta, fDphi);
 253 //
 254 //  I/O
 255     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
 256 //
 257 //
 258 }
 259
 260 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncelltot, Float_t etc[30000],
 261                              Float_t etac[30000], Float_t phic[30000],
 262                              Float_t minmove, Float_t maxmove, Int_t mode,
 263                              Float_t precbg,  Int_t   ierror)
 264 {
 265 // Wrapper for fortran coded jet finder
 266 // Full argument list
 267     jet_finder_ua1(ncell, ncelltot, etc, etac, phic,
 268                    minmove, maxmove, mode, precbg, ierror);
 269     // Write result to output
 270     if(fWrite) WriteJets();
 271     fEvent++;
 272 }
 273
 274 void AliEMCALJetFinder::Find()
 275 {
 276 // Wrapper for fortran coded jet finder using member data for
 277 // argument list
 278
 279     Float_t minmove = fMinMove;
 280     Float_t maxmove = fMaxMove;
 281     Int_t   mode     = fMode;
 282     Float_t precbg  = fPrecBg;
 283     Int_t   ierror;
 284     ResetJets(); // 4-feb-2002 by PAI
 285
 286     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 287                    minmove, maxmove, mode, precbg, ierror);
 288     fError = ierror;
 289     // Write result to output
 290     Int_t njet = Njets();
 291
 292     for (Int_t nj=0; nj<njet; nj++)
 293     {
 294         if (fWeightingMethod)
 295         {
 296           fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
 297                                       JetPhiW(nj),
 298                                       JetEtaW(nj));
 299
 300         }else{
 301
 302         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 303                                     JetPhiW(nj),
 304                                     JetEtaW(nj));
 305         }
 306         fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod);
 307         fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) );
 308         fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj)  ));
 309
 310     }
 311
 312     FindTracksInJetCone();
 313     if(fWrite) WriteJets();
 314     fEvent++;
 315 }
 316
 317
 318 Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
 319 {
 320 // Calculate the energy in the cone
 321 Float_t newenergy = 0.0;
 322 Float_t bineta,binphi;
 323 TAxis *x = fhLegoEMCAL->GetXaxis();
 324 TAxis *y = fhLegoEMCAL->GetYaxis();
 325 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 326         {
 327    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 328       {
 329          bineta = x->GetBinCenter(i);
 330          binphi = y->GetBinCenter(j);
 331        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 332       {
 333            newenergy += fhLegoEMCAL->GetBinContent(i,j);
 334       }
 335     }
 336 }
 337 return newenergy;
 338 }
 339
 340 Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
 341 {
 342 // Calculate the track energy in the cone
 343 Float_t newenergy = 0.0;
 344 Float_t bineta,binphi;
 345 TAxis *x = fhLegoTracks->GetXaxis();
 346 TAxis *y = fhLegoTracks->GetYaxis();
 347 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 348    {
 349    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 350    {
 351     bineta = x->GetBinCenter(i);
 352     binphi = y->GetBinCenter(j);
 353     if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 354       {
 355           newenergy += fhLegoTracks->GetBinContent(i,j);
 356        }
 357     }
 358 }
 359 return newenergy;
 360 }
 361
 362
 363 Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
 364 {
 365 // Calculate the weighted jet energy
 366
 367 Float_t newenergy = 0.0;
 368 Float_t bineta,binphi;
 369 TAxis *x = fhLegoEMCAL->GetXaxis();
 370 TAxis *y = fhLegoEMCAL->GetYaxis();
 371
 372
 373 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 374 {
 375    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 376    {
 377       bineta = x->GetBinCenter(i);
 378       binphi = y->GetBinCenter(j);
 379       if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 380       {
 381           newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
 382       }
 383     }
 384 }
 385
 386 return newenergy;
 387
 388 }
 389
 390
 391 Int_t AliEMCALJetFinder::Njets() const
 392 {
 393 // Get number of reconstructed jets
 394     return EMCALJETS.njet;
 395 }
 396
 397 Float_t AliEMCALJetFinder::JetEnergy(Int_t i) const
 398 {
 399 // Get reconstructed jet energy
 400     return EMCALJETS.etj[i];
 401 }
 402
 403 Float_t AliEMCALJetFinder::JetPhiL(Int_t i) const
 404 {
 405 // Get reconstructed jet phi from leading particle
 406     return EMCALJETS.phij[0][i];
 407 }
 408
 409 Float_t AliEMCALJetFinder::JetPhiW(Int_t i) const
 410 {
 411 // Get reconstructed jet phi from weighting
 412     return EMCALJETS.phij[1][i];
 413 }
 414
 415 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i) const
 416 {
 417 // Get reconstructed jet eta from leading particles
 418     return EMCALJETS.etaj[0][i];
 419 }
 420
 421
 422 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)  const
 423 {
 424 // Get reconstructed jet eta from weighting
 425     return EMCALJETS.etaj[1][i];
 426 }
 427
 428 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 429 {
 430 // Set grid cell size
 431     EMCALCELLGEO.etaCellSize = eta;
 432     EMCALCELLGEO.phiCellSize = phi;
 433 }
 434
 435 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 436 {
 437 // Set jet cone radius
 438     EMCALJETPARAM.coneRad = par;
 439     fConeRadius = par;
 440 }
 441
 442 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 443 {
 444 // Set et cut for seeds
 445     EMCALJETPARAM.etSeed = par;
 446     fEtSeed = par;
 447 }
 448
 449 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 450 {
 451 // Set minimum jet et
 452     EMCALJETPARAM.ejMin = par;
 453     fMinJetEt = par;
 454 }
 455
 456 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 457 {
 458 // Set et cut per cell
 459     EMCALJETPARAM.etMin = par;
 460     fMinCellEt = par;
 461 }
 462
 463 void AliEMCALJetFinder::SetPtCut(Float_t par)
 464 {
 465 // Set pT cut on charged tracks
 466     fPtCut = par;
 467 }
 468
 469
 470 void AliEMCALJetFinder::Test()
 471 {
 472 //
 473 // Test the finder call
 474 //
 475     const Int_t knmax = 30000;
 476     Int_t ncell      = 10;
 477     Int_t ncelltot  = 100;
 478
 479     Float_t etc[knmax];
 480     Float_t etac[knmax];
 481     Float_t phic[knmax];
 482     Float_t minmove = 0;
 483     Float_t maxmove = 0;
 484     Int_t   mode = 0;
 485     Float_t precbg = 0;
 486     Int_t   ierror = 0;
 487
 488
 489     Find(ncell, ncelltot, etc, etac, phic,
 490          minmove, maxmove, mode, precbg, ierror);
 491
 492 }
 493
 494 //
 495 //  I/O
 496 //
 497
 498 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 499 {
 500     //
 501     // Add a jet
 502     //
 503     TClonesArray &lrawcl = *fJets;
 504     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 505 }
 506
 507 void AliEMCALJetFinder::ResetJets()
 508 {
 509     //
 510     // Reset Jet List
 511     //
 512     fJets->Clear();
 513     fNjets = 0;
 514 }
 515
 516 void AliEMCALJetFinder::WriteJets()
 517 {
 518 //
 519 // Add all jets to the list
 520 //
 521     const Int_t kBufferSize = 4000;
 522     const char* file = 0;
 523
 524     Int_t njet = Njets();
 525
 526     for (Int_t nj = 0; nj < njet; nj++)
 527     {
 528         AddJet(*fJetT[nj]);
 529         delete fJetT[nj];
 530     }
 531
 532 // I/O
 533     AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 534     if (!fOutFileName) {
 535 //
 536 // output written to input file
 537 //
 538       AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 539       TTree* pK = gAlice->TreeK();
 540       file = (pK->GetCurrentFile())->GetName();
 541       TBranch * jetBranch ;
 542         if (fDebug > 1)
 543             printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
 544         //if (fJets && gAlice->TreeR()) {
 545         if (fJets && gime->TreeR()) {
 546           // pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 547           jetBranch = gime->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ;
 548           //pEMCAL->MakeBranchInTree(gime->TreeR(),
 549           //                         "EMCALJets",
 550           //                         &fJets,
 551           //                         kBufferSize,
 552           //                         file);
 553
 554           //Int_t nev = gAlice->GetHeader()->GetEvent();
 555           //gAlice->TreeR()->Fill();
 556           jetBranch->Fill();
 557           //char hname[30];
 558           //sprintf(hname,"TreeR%d", nev);
 559           //gAlice->TreeR()->Write(hname);
 560           //gAlice->TreeR()->Reset();
 561           gime->WriteRecPoints("OVERWRITE");
 562         }
 563     } else {
 564 //
 565 // Output written to user specified output file
 566 //
 567       //TTree* pK = gAlice->TreeK();
 568       TTree* pK = gAlice->TreeK();
 569       fInFile  = pK->GetCurrentFile();
 570
 571       fOutFile->cd();
 572       char hname[30];
 573       sprintf(hname,"TreeR%d", fEvent);
 574       TTree* treeJ = new TTree(hname, "EMCALJets");
 575       treeJ->Branch("EMCALJets", &fJets, kBufferSize);
 576       treeJ->Fill();
 577       treeJ->Write(hname);
 578       fInFile->cd();
 579     }
 580     ResetJets();
 581 }
 582
 583 void AliEMCALJetFinder::BookLego()
 584 {
 585 //
 586 //  Book histo for discretization
 587 //
 588
 589 //
 590 //  Don't add histos to the current directory
 591     if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
 592
 593     //    TH2::AddDirectory(0); // hists wil be put to the list from gROOT
 594     //    TH1::AddDirectory(0);
 595     gROOT->cd();
 596 //
 597 //  Signal map
 598     fLego = new TH2F("legoH","eta-phi",
 599                            fNbinEta, fEtaMin, fEtaMax,
 600                            fNbinPhi, fPhiMin, fPhiMax);
 601 //
 602 //  Background map
 603     fLegoB = new TH2F("legoB","eta-phi for BG event",
 604                            fNbinEta, fEtaMin, fEtaMax,
 605                            fNbinPhi, fPhiMin, fPhiMax);
 606
 607 //  Tracks map
 608     fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
 609     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 610 //  EMCAL map
 611     fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
 612     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 613 //  Hadron correction map
 614     fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
 615     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 616 //  Hists. for tuning jet finder
 617     fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
 618
 619     TArrayF eTmp(1101);
 620     eTmp[0] = 0.0;
 621     for(Int_t i=1; i<=1000; i++)      eTmp[i] = 0.1*i; // step 100 mev
 622     for(Int_t i=1001; i<=1100; i++)   eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
 623
 624     fhCellEt  = new TH1F("hCellEt","Cell E_{T} from fLego",
 625     eTmp.GetSize()-1, eTmp.GetArray());
 626     fhCellEMCALEt  = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
 627     eTmp.GetSize()-1, eTmp.GetArray());
 628     fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
 629     eTmp.GetSize()-1, eTmp.GetArray());
 630     fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
 631     eTmp.GetSize()-1, eTmp.GetArray());
 632
 633     fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
 634     "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
 635
 636     fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
 637     TAxis *xax = fhSinTheta->GetXaxis();
 638     for(Int_t i=1; i<=fNbinEta; i++) {
 639       Double_t eta = xax->GetBinCenter(i);
 640       fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
 641     }
 642
 643     //! first canvas for drawing
 644     fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
 645 }
 646
 647 void AliEMCALJetFinder::DumpLego()
 648 {
 649 //
 650 // Dump lego histo into array
 651 //
 652     fNcell = 0;
 653     TAxis* xaxis = fLego->GetXaxis();
 654     TAxis* yaxis = fLego->GetYaxis();
 655     //    fhCellEt->Clear();
 656     Float_t e, eH;
 657     for (Int_t i = 1; i <= fNbinEta; i++) {
 658         for (Int_t j = 1; j <= fNbinPhi; j++) {
 659
 660             e = fLego->GetBinContent(i,j);
 661             if (fRandomBg) {
 662                 if (gRandom->Rndm() < 0.5) {
 663                     Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.));
 664                     e += ebg;
 665                 }
 666             }
 667             if (e > 0.0) e -= fMinCellEt;
 668             if (e < 0.0) e = 0.;
 669             Float_t eta  = xaxis->GetBinCenter(i);
 670             Float_t phi  = yaxis->GetBinCenter(j);
 671             fEtCell[fNcell]  = e;
 672             fEtaCell[fNcell] = eta;
 673             fPhiCell[fNcell] = phi;
 674             fNcell++;
 675             fhCellEt->Fill(e);
 676             if(fhLegoEMCAL) {
 677               eH = fhLegoEMCAL->GetBinContent(i,j);
 678               if(eH > 0.0) fhCellEMCALEt->Fill(eH);
 679             }
 680         } // phi
 681     } // eta
 682 //  today
 683 //    fNcell--;
 684 }
 685
 686 void AliEMCALJetFinder::ResetMap()
 687 {
 688 //
 689 // Reset eta-phi array
 690
 691     for (Int_t i=0; i<30000; i++)
 692     {
 693         fEtCell[i]  = 0.;
 694         fEtaCell[i] = 0.;
 695         fPhiCell[i] = 0.;
 696     }
 697 }
 698
 699
 700 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 701 {
 702 // Which generator
 703 //
 704     const char* name = gAlice->Generator()->GetName();
 705     enum {kPythia, kHijing, kHijingPara};
 706     Int_t genType = 0;
 707
 708     if (!strcmp(name, "Hijing")){
 709         genType = kHijing;
 710     } else if (!strcmp(name, "Pythia")) {
 711         genType = kPythia;
 712     } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) {
 713         genType = kHijingPara;
 714     }
 715     if (fDebug>=2)
 716         printf("Fill tracks generated by %s %d\n", name, genType);
 717
 718
 719 //
 720 // Fill Cells with track information
 721 //
 722     if (fDebug >= 2)
 723     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
 724     fNChTpc = 0;
 725
 726     ResetMap();
 727
 728     if (!fLego) BookLego();
 729 // Reset
 730     if (flag == 0) fLego->Reset();
 731 //
 732 // Access particle information
 733     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 734     Int_t ntr   = (gAlice->GetHeader())->GetNtrack();
 735     printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n",
 736     npart, ntr, fLego->Integral());
 737
 738 // Create track list
 739 //
 740 // 0: not selected
 741 // 1: selected for jet finding
 742 // 2: inside jet
 743 // ....
 744     if (fTrackList) delete[] fTrackList;
 745     if (fPtT)       delete[] fPtT;
 746     if (fEtaT)      delete[] fEtaT;
 747     if (fPhiT)      delete[] fPhiT;
 748     if (fPdgT)      delete[] fPdgT;
 749
 750     fTrackList = new Int_t  [npart];
 751     fPtT       = new Float_t[npart];
 752     fEtaT      = new Float_t[npart];
 753     fPhiT      = new Float_t[npart];
 754     fPdgT      = new Int_t[npart];
 755
 756     fNt   = npart;
 757     fNtS  = 0;
 758     Float_t chTmp=0.0; // charge of current particle
 759     //    Int_t idGeant;
 760
 761     // this is for Pythia ??
 762     for (Int_t part = 0; part < npart; part++) {
 763         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
 764         Int_t mpart   = mPart->GetPdgCode();
 765         Int_t child1  = mPart->GetFirstDaughter();
 766         Float_t pT    = mPart->Pt();
 767         Float_t p     = mPart->P();
 768         Float_t phi   = mPart->Phi();
 769         Float_t eta   = -100.;
 770         if(pT > 0.001) eta   = mPart->Eta();
 771         Float_t theta = mPart->Theta();
 772         if  (fDebug>=2 && mPart->GetStatusCode()==1) {
 773            printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
 774            part, mpart, child1, mPart->GetLastDaughter(), mPart->GetStatusCode());
 775         }
 776
 777         if (fDebug >= 2 && genType == kPythia) {
 778             if (part == 6 || part == 7)
 779             {
 780                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
 781                        part-5, pT, eta, phi);
 782             }
 783         }
 784
 785         fTrackList[part] = 0;
 786         fPtT[part]       = pT; // must be change after correction for resolution !!!
 787         fEtaT[part]      = eta;
 788         fPhiT[part]      = phi;
 789         fPdgT[part]      = mpart;
 790
 791 //      final state particles only
 792
 793         if (genType == kPythia) {
 794             if (mPart->GetStatusCode() != 1) continue;
 795         } else if (genType == kHijing) {
 796             if (child1 >= 0 && child1 < npart) continue;
 797         }
 798
 799
 800         TParticlePDG* pdgP = 0;
 801 // charged or neutral
 802         pdgP  = mPart->GetPDG();
 803         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
 804
 805         if (ich == 0) {
 806             if (chTmp == 0) {
 807                 if (!fK0N) {
 808                     continue;
 809                 } else {
 810                     if (mpart != kNeutron    &&
 811                         mpart != kNeutronBar &&
 812                         mpart != kK0Long) continue;
 813                 }
 814             }
 815         } else if (ich == 2) {
 816           if (mpart == kNeutron    ||
 817               mpart == kNeutronBar ||
 818               mpart == kK0Long) continue;
 819         }
 820
 821         if (TMath::Abs(eta)<=0.9) fNChTpc++;
 822 // final state only
 823         if (child1 >= 0 && child1 < npart) continue;
 824 // acceptance cut
 825         if (eta > fEtaMax || eta < fEtaMin)    continue;
 826         if (phi > fPhiMax || phi < fPhiMin )   continue;
 827 //
 828         if (fDebug >= 3)
 829         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
 830         part, mpart, child1, eta, phi, pT, chTmp);
 831 //
 832 //
 833 // Momentum smearing goes here ...
 834 //
 835         fhTrackPt->Fill(pT);
 836         Float_t pw;
 837         if (fSmear && TMath::Abs(chTmp)) {
 838             pw = AliEMCALFast::SmearMomentum(1,p);
 839             // p changed - take into account when calculate pT,
 840             // pz and so on ..  ?
 841             pT = (pw/p) * pT;
 842             if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
 843             p  = pw;
 844         }
 845 //
 846 // Tracking Efficiency and TPC acceptance goes here ...
 847         Float_t eff;
 848         if (fEffic && TMath::Abs(chTmp)) {
 849             eff = 0.9; //  AliEMCALFast::Efficiency(2,p);
 850             if(fhEff) fhEff->Fill(p, eff);
 851             if (AliEMCALFast::RandomReject(eff)) {
 852                 if(fDebug >= 5) printf(" reject due to unefficiency ");
 853                 continue;
 854             }
 855         }
 856 //
 857 // Correction of Hadronic Energy goes here ...
 858 //
 859 //
 860 // phi propagation for hadronic correction
 861
 862         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
 863         Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
 864         if(TMath::Abs(chTmp)) {
 865             // hadr. correction only for charge particle
 866             dphi  = PropagatePhi(pT, chTmp, curls);
 867             phiHC = phi + dphi;
 868             if (fDebug >= 6) {
 869                 printf("\n Delta phi %f pT %f ", dphi, pT);
 870                 if (curls) printf("\n !! Track is curling");
 871             }
 872             if(!curls) fhTrackPtBcut->Fill(pT);
 873
 874             if (fHCorrection && !curls) {
 875                 if (!fHadronCorrector)
 876                     Fatal("AliEMCALJetFinder",
 877                           "Hadronic energy correction required but not defined !");
 878
 879                 dpH    = fHadronCorrector->GetEnergy(p, eta, 7);
 880                 eTdpH  = dpH*TMath::Sin(theta);
 881
 882                 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
 883                                         phi, phiHC, -eTdpH); // correction is negative
 884                 Int_t xbin,ybin;
 885                 xbin = fLego->GetXaxis()->FindBin(eta);
 886                 ybin = fLego->GetYaxis()->FindBin(phiHC);
 887                 cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 888                 fLego->Fill(eta, phi, -fSamplingF*eTdpH );
 889                 cout <<"Hadron Correction affected bin - contents after  correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 890                 fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
 891             }
 892         }
 893 //
 894 //  More to do ? Just think about it !
 895 //
 896         if (phi > fPhiMax || phi < fPhiMin )   continue;
 897
 898         if(TMath::Abs(chTmp) ) { // charge particle
 899             if (pT > fPtCut && !curls) {
 900                 if (fDebug >= 8) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 901                                         eta , phi, pT, fNtS);
 902                 fLego->Fill(eta, phi, pT);
 903                 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
 904                 fTrackList[part] = 1;
 905                 fNtS++;
 906             }
 907         } else if(ich > 0 || fK0N) {
 908             // case of n, nbar and K0L
 909             if (fDebug >= 9) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 910                                     eta , phi, pT, fNtS);
 911             fLego->Fill(eta, phi, pT);
 912             fTrackList[part] = 1;
 913             fNtS++;
 914         }
 915     } // primary loop
 916     Float_t etsum = 0.;
 917     for(Int_t i=0; i<fLego->GetSize(); i++) {
 918         Float_t etc =  (*fLego)[i];
 919         if (etc > fMinCellEt) etsum += etc;
 920     }
 921
 922     printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
 923     printf("                Track selected(fNtS) %i \n", fNtS);
 924
 925     DumpLego();
 926 }
 927
 928 void AliEMCALJetFinder::FillFromHits(Int_t flag)
 929 {
 930 //
 931 // Fill Cells with hit information
 932 //
 933 //
 934     if (fDebug >= 2)
 935     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
 936
 937     ResetMap();
 938
 939     if (!fLego) BookLego();
 940 //  Reset eta-phi maps if needed
 941     if (flag == 0) { // default behavior
 942       fLego->Reset();
 943       fhLegoTracks->Reset();
 944       fhLegoEMCAL->Reset();
 945       fhLegoHadrCorr->Reset();
 946     }
 947 //  Initialize from background event if available
 948 //
 949 // Access hit information
 950     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 951     AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 952     TTree *treeH = gime->TreeH();
 953     Int_t ntracks = (Int_t) treeH->GetEntries();
 954 //
 955 //   Loop over tracks
 956 //
 957     Int_t nbytes = 0;
 958     //    Double_t etH = 0.0;
 959
 960     for (Int_t track=0; track<ntracks;track++) {
 961         gAlice->ResetHits();
 962         nbytes += treeH->GetEvent(track);
 963 //
 964 //  Loop over hits
 965 //
 966         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 967             mHit;
 968             mHit=(AliEMCALHit*) pEMCAL->NextHit())
 969         {
 970             Float_t x      =    mHit->X();         // x-pos of hit
 971             Float_t y      =    mHit->Y();         // y-pos
 972             Float_t z      =    mHit->Z();         // z-pos
 973             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 974
 975             Float_t r      =    TMath::Sqrt(x*x+y*y);
 976             Float_t theta  =    TMath::ATan2(r,z);
 977             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 978             Float_t phi    =    TMath::ATan2(y,x);
 979
 980             if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
 981
 982             //            etH = fSamplingF*eloss*TMath::Sin(theta);
 983             fLego->Fill(eta, phi, eloss);
 984         } // Hit Loop
 985     } // Track Loop
 986
 987     // Transition from deposit energy to eT (eT = de*SF*sin(theta))
 988     Double_t etsum = 0;
 989     for(Int_t i=1; i<=fLego->GetNbinsX(); i++){   // eta
 990        Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
 991        for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
 992           eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
 993           fLego->SetBinContent(i,j,eT);
 994     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
 995           fhLegoEMCAL->SetBinContent(i,j,eT);
 996           if (eT > fMinCellEt) etsum += eT;
 997        }
 998     }
 999
1000     //    for(Int_t i=0; i<fLego->GetSize(); i++) {
1001     //  (*fhLegoEMCAL)[i] = (*fLego)[i];
1002     //  Float_t etc =  (*fLego)[i];
1003     //  if (etc > fMinCellEt) etsum += etc;
1004     //    }
1005
1006     printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
1007     //    DumpLego(); ??
1008 }
1009
1010 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
1011 {
1012 //
1013 // Fill Cells with digit information
1014 //
1015 //
1016     ResetMap();
1017
1018     if (!fLego) BookLego();
1019     if (flag == 0) fLego->Reset();
1020     Int_t nbytes=0;
1021
1022
1023 //
1024 //  Connect digits
1025 //
1026     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
1027     TTree *treeD = gAlice->TreeD();
1028     TBranchElement* branchDg = (TBranchElement*)
1029         treeD->GetBranch("EMCAL");
1030
1031     if (!branchDg) Fatal("AliEMCALJetFinder",
1032                          "Reading digits requested but no digits in file !");
1033
1034     branchDg->SetAddress(&digs);
1035     Int_t nent = (Int_t) branchDg->GetEntries();
1036 //
1037 //  Connect digitizer
1038 //
1039     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
1040     TBranchElement* branchDr = (TBranchElement*)
1041         treeD->GetBranch("AliEMCALDigitizer");
1042     branchDr->SetAddress(&digr);
1043 //
1044 //
1045     nbytes += branchDg->GetEntry(0);
1046     nbytes += branchDr->GetEntry(0);
1047 //
1048 //  Get digitizer parameters
1049     Float_t ecADCped  = digr->GetECApedestal();
1050     Float_t ecADCcha  = digr->GetECAchannel();
1051
1052     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1053     AliEMCALGeometry* geom =
1054         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
1055
1056     if (fDebug) {
1057         Int_t   ndig = digs->GetEntries();
1058         Info("FillFromDigits","Number of Digits: %d %d\n Parameters: EC: %f %f\n Geometry: %d %d",
1059              ndig, nent, ecADCped, ecADCcha, geom->GetNEta(), geom->GetNPhi());
1060     }
1061
1062 //
1063 //  Loop over digits
1064     AliEMCALDigit* sdg;
1065     TIter next(digs);
1066     while ((sdg = (AliEMCALDigit*)(next())))
1067     {
1068         Double_t pedestal = 0.;
1069         Double_t channel  = 0.;
1070         if (geom->IsInECA(sdg->GetId()))  {
1071           pedestal = ecADCped;
1072           channel  = ecADCcha;
1073         }
1074         else
1075           Fatal("FillFromDigits", "unexpected digit is number!") ;
1076
1077         Float_t eta = sdg->GetEta();
1078         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1079         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1080
1081         if (fDebug > 1)
1082           Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1083                            eta, phi, amp, sdg->GetAmp());
1084
1085         fLego->Fill(eta, phi, fSamplingF*amp);
1086     } // digit loop
1087 //
1088 //  Dump lego hist
1089     DumpLego();
1090 }
1091
1092
1093 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1094 {
1095 //
1096 // Fill Cells with hit information
1097 //
1098 //
1099     ResetMap();
1100
1101     if (!fLego) BookLego();
1102 // Reset
1103     if (flag == 0) fLego->Reset();
1104
1105 // Flag charged tracks passing through TPC which
1106 // are associated to EMCAL Hits
1107     BuildTrackFlagTable();
1108
1109 //
1110 // Access particle information
1111     TTree *treeK = gAlice->TreeK();
1112     Int_t ntracks = (Int_t) treeK->GetEntries();
1113
1114     if (fPtT)       delete[] fPtT;
1115     if (fEtaT)      delete[] fEtaT;
1116     if (fPhiT)      delete[] fPhiT;
1117     if (fPdgT)      delete[] fPdgT;
1118
1119     fPtT       = new Float_t[ntracks];
1120     fEtaT      = new Float_t[ntracks];
1121     fPhiT      = new Float_t[ntracks];
1122     fPdgT      = new Int_t[ntracks];
1123
1124     fNt   = ntracks;
1125     fNtS  = 0;
1126
1127     for (Int_t track = 0; track < ntracks; track++) {
1128           TParticle *mPart = gAlice->GetMCApp()->Particle(track);
1129           Float_t pT    = mPart->Pt();
1130           Float_t phi   = mPart->Phi();
1131           Float_t eta   = mPart->Eta();
1132
1133         if(fTrackList[track]) {
1134           fPtT[track]       = pT;
1135           fEtaT[track]      = eta;
1136           fPhiT[track]      = phi;
1137           fPdgT[track]      = mPart->GetPdgCode();
1138
1139           if (track < 2) continue;      //Colliding particles?
1140           if (pT == 0 || pT < fPtCut) continue;
1141           fNtS++;
1142           fLego->Fill(eta, phi, pT);
1143         }
1144       } // track loop
1145     DumpLego();
1146 }
1147
1148 void AliEMCALJetFinder::FillFromParticles()
1149 {
1150 // 26-feb-2002 PAI - for checking all chain
1151 // Work on particles level; accept all particle (not neutrino )
1152
1153     Double_t pX=0, pY=0, pZ=0, energy=0; // checking conservation law
1154     fNChTpc = 0;
1155
1156     ResetMap();
1157     if (!fLego) BookLego();
1158     fLego->Reset();
1159 //
1160 // Access particles information
1161     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1162     if (fDebug >= 2 || npart<=0) {
1163        printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1164        if(npart<=0) return;
1165     }
1166     fNt   = npart;
1167     fNtS  = 0;
1168     RearrangeParticlesMemory(npart);
1169
1170 //  Go through the particles
1171     Int_t mpart, child1, child2, geantPdg;
1172     Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1173     TParticle *mPart=0;
1174     for (Int_t part = 0; part < npart; part++) {
1175
1176         fTrackList[part] = 0;
1177
1178         mPart   = gAlice->GetMCApp()->Particle(part);
1179         mpart   = mPart->GetPdgCode();
1180         child1  = mPart->GetFirstDaughter();
1181         child2  = mPart->GetLastDaughter();
1182         pT      = mPart->Pt();
1183         phi     = mPart->Phi();
1184         eta     = mPart->Eta();
1185
1186         px      = mPart->Px();
1187         py      = mPart->Py();
1188         pz      = mPart->Pz();
1189         e       = mPart->Energy();
1190
1191 // see pyedit in Pythia's text
1192         geantPdg = mpart;
1193         if (IsThisPartonsOrDiQuark(mpart)) continue;
1194         printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1195         part, mpart, geantPdg, px, py, pz, e, child1, child2);
1196
1197 //  exclude partons (21 - gluon, 92 - string)
1198
1199
1200 // exclude neutrinous also ??
1201         if (fDebug >= 11 && pT>0.01)
1202         printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1203         part, mpart, eta, phi, pT);
1204
1205         fPtT[part]       = pT;
1206         fEtaT[part]      = eta;
1207         fPhiT[part]      = phi;
1208         fPdgT[part]      = mpart;
1209         fNtS++;
1210
1211 // final state only
1212         if (child1 >= 0 && child1 < npart) continue;
1213
1214         //        printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1215         //        part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1216         pX += px;
1217         pY += py;
1218         pZ += pz;
1219         energy  += e;
1220
1221
1222         if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1223 // acceptance cut
1224         if (eta > fEtaMax || eta < fEtaMin)    continue;
1225         if (phi > fPhiMax || phi < fPhiMin )   continue;
1226 //
1227         if(fK0N==0 ) { // exclude neutral hadrons
1228           if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1229         }
1230         fTrackList[part] = 1;
1231         fLego->Fill(eta, phi, pT);
1232
1233     } // primary loop
1234     printf("\n                PX %8.2f  PY %8.2f  PZ %8.2f  E %8.2f \n",
1235     pX, pY, pZ, energy);
1236     DumpLego();
1237     if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1238 }
1239
1240 void AliEMCALJetFinder::FillFromPartons()
1241 {
1242 // function under construction - 13-feb-2002 PAI
1243
1244     if (fDebug >= 2)
1245     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1246     //
1247
1248     ResetMap();
1249     if (!fLego) BookLego();
1250     fLego->Reset();
1251 //
1252 // Access particle information
1253     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1254     if (fDebug >= 2 || npart<=0)
1255     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1256     fNt   = 0; // for FindTracksInJetCone
1257     fNtS  = 0;
1258
1259 //  Go through the partons
1260     Int_t statusCode=0;
1261     for (Int_t part = 8; part < npart; part++) {
1262         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
1263         Int_t mpart   = mPart->GetPdgCode();
1264         //      Int_t child1  = MPart->GetFirstDaughter();
1265         Float_t pT    = mPart->Pt();
1266         //      Float_t p     = MPart->P();
1267         Float_t phi   = mPart->Phi();
1268         Float_t eta   = mPart->Eta();
1269         //      Float_t theta = MPart->Theta();
1270         statusCode    = mPart->GetStatusCode();
1271
1272 // accept partons (21 - gluon, 92 - string)
1273         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1274         if (fDebug > 1 && pT>0.01)
1275         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
1276         part, mpart, statusCode, eta, phi, pT);
1277         //      if (fDebug >= 3) MPart->Print();
1278 // accept partons before fragmentation - p.57 in Pythia manual
1279 //        if(statusCode != 1) continue;
1280 // acceptance cut
1281         if (eta > fEtaMax || eta < fEtaMin)    continue;
1282         if (phi > fPhiMax || phi < fPhiMin )   continue;
1283 // final state only
1284 //      if (child1 >= 0 && child1 < npart) continue;
1285 //
1286 //
1287         fLego->Fill(eta, phi, pT);
1288
1289     } // primary loop
1290     DumpLego();
1291 }
1292
1293 void AliEMCALJetFinder::BuildTrackFlagTable() {
1294
1295 // Method to generate a lookup table for TreeK particles
1296 // which are linked to hits in the EMCAL
1297 //
1298 // --Author: J.L. Klay
1299 //
1300 // Access hit information
1301     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1302
1303     TTree *treeK = gAlice->TreeK();             // Get the number of entries in the kine tree
1304     Int_t nKTrks = (Int_t) treeK->GetEntries(); // (Number of particles created somewhere)
1305
1306     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
1307     fTrackList = new Int_t[nKTrks];             //before generating a new one
1308
1309     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
1310         fTrackList[i] = 0;
1311     }
1312
1313     AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
1314     // TTree *treeH = gAlice->TreeH();
1315     TTree *treeH = gime->TreeH();
1316    Int_t ntracks = (Int_t) treeH->GetEntries();
1317 //
1318 //   Loop over tracks
1319 //
1320     Int_t nbytes = 0;
1321
1322     for (Int_t track=0; track<ntracks;track++) {
1323         gAlice->ResetHits();
1324         nbytes += treeH->GetEvent(track);
1325         if (pEMCAL)  {
1326
1327 //
1328 //  Loop over hits
1329 //
1330             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1331                 mHit;
1332                 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1333             {
1334                 Int_t   iTrk   =    mHit->Track();       // track number
1335                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1336
1337                 //Determine the origin point of this particle - it made a hit in the EMCAL
1338                 TParticle         *trkPart = gAlice->GetMCApp()->Particle(iTrk);
1339                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1340                 Int_t      trkCode = trkPart->GetPdgCode();
1341                 Double_t           trkChg;
1342                 if (trkCode < 10000) {          //Big Ions cause problems for
1343                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1344                 } else {                                //likely to make it very far, set
1345                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1346                 }
1347                 Float_t            vxTrk   = trkPart->Vx();
1348                 Float_t            vyTrk   = trkPart->Vy();
1349                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1350                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1351
1352                 //Loop through the ancestry of the EMCAL entrance particles
1353                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1354                 while (ancestor != -1) {
1355                     TParticle     *ancPart = gAlice->GetMCApp()->Particle(ancestor);  //get particle info on ancestor
1356                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1357                     Int_t             ancCode = ancPart->GetPdgCode();
1358                     Double_t       ancChg;
1359                     if (ancCode < 10000) {
1360                         ancChg  = ancPdg->Charge();
1361                     } else {
1362                         ancChg  = 0.0;
1363                     }
1364                     Float_t           vxAnc   = ancPart->Vx();
1365                     Float_t           vyAnc   = ancPart->Vy();
1366                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1367                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1368                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1369                 }
1370
1371                 //Determine the origin point of the primary particle associated with the hit
1372                 TParticle     *primPart = gAlice->GetMCApp()->Particle(idprim);
1373                 TParticlePDG  *primPdg  = primPart->GetPDG();
1374                 Int_t      primCode = primPart->GetPdgCode();
1375                 Double_t           primChg;
1376                 if (primCode < 10000) {
1377                     primChg  = primPdg->Charge();
1378                 } else {
1379                     primChg  = 0.0;
1380                 }
1381                 Float_t    vxPrim = primPart->Vx();
1382                 Float_t    vyPrim = primPart->Vy();
1383                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1384                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1385             } // Hit Loop
1386         } //if (pEMCAL)
1387     } // Track Loop
1388 }
1389
1390 Int_t AliEMCALJetFinder
1391 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1392 // Set the flag for the track
1393     Int_t flag    = 0;
1394     Int_t parton  = 0;
1395     Int_t neutral = 0;
1396
1397     if (charge == 0) neutral = 1;
1398
1399     if (TMath::Abs(code) <= 6   ||
1400         code == 21              ||
1401         code == 92) parton = 1;
1402
1403     //It's not a parton, it's charged and it went through the TPC
1404     if (!parton && !neutral && radius <= 84.0) flag = 1;
1405
1406     return flag;
1407 }
1408
1409
1410
1411 void AliEMCALJetFinder::SaveBackgroundEvent(Char_t *name)
1412 {
1413 // Saves the eta-phi lego and the tracklist and name of file with BG events
1414 //
1415     if (fLegoB) {
1416        fLegoB->Reset();
1417        (*fLegoB) = (*fLegoB) + (*fLego);
1418        //       if(fDebug)
1419        printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1420        fLegoB->Integral(), fLego->Integral());
1421     }
1422
1423     if (fPtB)        delete[] fPtB;
1424     if (fEtaB)       delete[] fEtaB;
1425     if (fPhiB)       delete[] fPhiB;
1426     if (fPdgB)       delete[] fPdgB;
1427     if (fTrackListB) delete[] fTrackListB;
1428
1429     fPtB          = new Float_t[fNtS];
1430     fEtaB         = new Float_t[fNtS];
1431     fPhiB         = new Float_t[fNtS];
1432     fPdgB         = new Int_t  [fNtS];
1433     fTrackListB   = new Int_t  [fNtS];
1434
1435     fNtB = 0;
1436
1437     for (Int_t i = 0; i < fNt; i++) {
1438         if (!fTrackList[i]) continue;
1439         fPtB [fNtB]       = fPtT [i];
1440         fEtaB[fNtB]       = fEtaT[i];
1441         fPhiB[fNtB]       = fPhiT[i];
1442         fPdgB[fNtB]       = fPdgT[i];
1443         fTrackListB[fNtB] = 1;
1444         fNtB++;
1445     }
1446     fBackground = 1;
1447     printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1448
1449     if(strlen(name) == 0) {
1450        TSeqCollection *li = gROOT->GetListOfFiles();
1451        TString nf;
1452        for(Int_t i=0; i<li->GetSize(); i++) {
1453           nf = ((TFile*)li->At(i))->GetName();
1454           if(nf.Contains("backgorund")) break;
1455        }
1456        fBGFileName = nf;
1457     } else {
1458        fBGFileName = name;
1459     }
1460     printf("BG file name is \n %s\n", fBGFileName.Data());
1461 }
1462
1463 void AliEMCALJetFinder::InitFromBackground()
1464 {
1465 //
1466 //
1467     if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1468
1469     if (fLego) {
1470         fLego->Reset();
1471         (*fLego) = (*fLego) + (*fLegoB);
1472         if(fDebug)
1473             printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1474                    fLego->Integral(), fLegoB->Integral());
1475     } else {
1476         printf(" => fLego undefined \n");
1477     }
1478 }
1479
1480
1481 void AliEMCALJetFinder::FindTracksInJetCone()
1482 {
1483 //
1484 //  Build list of tracks inside jet cone
1485 //
1486 //  Loop over jets
1487     Int_t njet = Njets();
1488     for (Int_t nj = 0; nj < njet; nj++)
1489     {
1490         Float_t etaj = JetEtaW(nj);
1491         Float_t phij = JetPhiW(nj);
1492         Int_t   nT   = 0;           // number of associated tracks
1493
1494 // Loop over particles in current event
1495         for (Int_t part = 0; part < fNt; part++) {
1496             if (!fTrackList[part]) continue;
1497             Float_t phi      = fPhiT[part];
1498             Float_t eta      = fEtaT[part];
1499             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1500                                            (phij-phi)*(phij-phi));
1501             if (dr < fConeRadius) {
1502                 fTrackList[part] = nj+2;
1503                 nT++;
1504             } // < ConeRadius ?
1505         } // particle loop
1506
1507 // Same for background event if available
1508         Int_t nTB = 0;
1509         if (fBackground) {
1510             for (Int_t part = 0; part < fNtB; part++) {
1511                 Float_t phi      = fPhiB[part];
1512                 Float_t eta      = fEtaB[part];
1513                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1514                                                (phij-phi)*(phij-phi));
1515                 fTrackListB[part] = 1;
1516
1517                 if (dr < fConeRadius) {
1518                     fTrackListB[part] = nj+2;
1519                     nTB++;
1520                 } // < ConeRadius ?
1521             } // particle loop
1522         } // Background available ?
1523
1524         Int_t nT0 = nT + nTB;
1525         printf("Total number of tracks %d\n", nT0);
1526
1527         if (nT0 > 1000) nT0 = 1000;
1528
1529         Float_t* ptT  = new Float_t[nT0];
1530         Float_t* etaT = new Float_t[nT0];
1531         Float_t* phiT = new Float_t[nT0];
1532         Int_t*   pdgT = new Int_t[nT0];
1533
1534         Int_t iT = 0;
1535         Int_t j;
1536
1537         for (Int_t part = 0; part < fNt; part++) {
1538             if (fTrackList[part] == nj+2) {
1539                 Int_t index = 0;
1540                 for (j=iT-1; j>=0; j--) {
1541                     if (fPtT[part] > ptT[j]) {
1542                         index = j+1;
1543                         break;
1544                     }
1545                 }
1546                 for (j=iT-1; j>=index; j--) {
1547                     ptT [j+1]  = ptT [j];
1548                     etaT[j+1]  = etaT[j];
1549                     phiT[j+1]  = phiT[j];
1550                     pdgT[j+1]  = pdgT[j];
1551                 }
1552                 ptT [index] = fPtT [part];
1553                 etaT[index] = fEtaT[part];
1554                 phiT[index] = fPhiT[part];
1555                 pdgT[index] = fPdgT[part];
1556                 iT++;
1557             } // particle associated
1558             if (iT > nT0) break;
1559         } // particle loop
1560
1561         if (fBackground) {
1562             for (Int_t part = 0; part < fNtB; part++) {
1563                 if (fTrackListB[part] == nj+2) {
1564                     Int_t index = 0;
1565                     for (j=iT-1; j>=0; j--) {
1566                         if (fPtB[part] > ptT[j]) {
1567                             index = j+1;
1568
1569                             break;
1570                         }
1571                     }
1572                     for (j=iT-1; j>=index; j--) {
1573                         ptT [j+1]  = ptT [j];
1574                         etaT[j+1]  = etaT[j];
1575                         phiT[j+1]  = phiT[j];
1576                         pdgT[j+1]  = pdgT[j];
1577                     }
1578                     ptT [index] = fPtB [part];
1579                     etaT[index] = fEtaB[part];
1580                     phiT[index] = fPhiB[part];
1581                     pdgT[index] = fPdgB[part];
1582                     iT++;
1583                 } // particle associated
1584                 if (iT > nT0) break;
1585             } // particle loop
1586         } // Background available ?
1587
1588         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1589         delete[] ptT;
1590         delete[] etaT;
1591         delete[] phiT;
1592         delete[] pdgT;
1593
1594     } // jet loop loop
1595 }
1596
1597 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1598 {
1599 // Propagates phi angle to EMCAL radius
1600 //
1601   static Float_t b = 0.0, rEMCAL = -1.0;
1602 // Get field in kGS
1603   b = gAlice->Field()->SolenoidField();
1604 // Get EMCAL radius in cm
1605   rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1606   Float_t dPhi = 0.;
1607 //
1608 //
1609 // bending radies
1610 // pt [Gev]
1611 // B  [kG]
1612 //
1613     Float_t rB = 3335.6 * pt / b;  // [cm]  (case of |charge|=1)
1614 //
1615 // check if particle is curling below EMCAL
1616     if (2.*rB < rEMCAL) {
1617         curls = kTRUE;
1618         return dPhi;
1619     }
1620 //
1621 // if not calculate delta phi
1622     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1623     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1624     dPhi = -TMath::Sign(dPhi, charge);
1625 //
1626     return dPhi;
1627 }
1628
1629 void hf1(Int_t& , Float_t& , Float_t& )
1630 {
1631 // dummy for hbook calls
1632     ;
1633 }
1634
1635 void AliEMCALJetFinder::DrawLego(Char_t *opt)
1636 {
1637 // Draw lego
1638         if(fLego) fLego->Draw(opt);
1639 }
1640
1641 void  AliEMCALJetFinder::DrawLegoBackground(Char_t *opt)
1642 {
1643 // Draw background lego
1644         if(fLegoB) fLegoB->Draw(opt);
1645 }
1646
1647 void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
1648 {
1649 // Draw EMCAL Lego
1650         if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);
1651 }
1652
1653 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1654 {
1655 // Draw all hists
1656   static TPaveText *varLabel=0;
1657   if(!fC1) {
1658     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1659   }
1660   fC1->Clear();
1661   fC1->Divide(2,2);
1662   fC1->cd(1);
1663   gPad->SetLogy(1);
1664   fhCellEt->Draw();
1665
1666   fC1->cd(2);
1667   gPad->SetLogy(1);
1668   fhCellEMCALEt->Draw();
1669
1670   fC1->cd(3);
1671   gPad->SetLogy(1);
1672   fhTrackPt->Draw();
1673   fhTrackPtBcut->SetLineColor(2);
1674   fhTrackPtBcut->Draw("same");
1675
1676   fC1->cd(4);
1677   if(!varLabel) {
1678     PrintParameters(1);
1679     varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1680     varLabel->SetTextAlign(12);
1681     varLabel->SetFillColor(19); // see TAttFill
1682     TString tmp(GetTitle());
1683     varLabel->ReadFile(GetFileNameForParameters());
1684   }
1685   varLabel->Draw();
1686   fC1->Update();
1687   if(mode) { // for saving picture to the file
1688     TString stmp(GetFileNameForParameters());
1689     stmp.ReplaceAll("_Par.txt",".ps");
1690     fC1->Print(stmp.Data());
1691   }
1692 }
1693
1694 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1695 {
1696 // Print all parameters out
1697
1698   FILE *file=0;
1699   if(mode==0) file = stdout; // output to terminal
1700   else {
1701     file = fopen(GetFileNameForParameters(),"w");
1702     if(file==0) file = stdout;
1703   }
1704   fprintf(file,"====   Filling lego   ==== \n");
1705   fprintf(file,"Sampling fraction %6.3f  ", fSamplingF);
1706   fprintf(file,"Smearing          %6i  ", fSmear);
1707   fprintf(file,"Efficiency        %6i\n", fEffic);
1708   fprintf(file,"Hadr.Correct.     %6i  ", fHCorrection);
1709   fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1710   fprintf(file,"====  Jet finding     ==== \n");
1711   fprintf(file,"Cone radius       %6.2f  ", fConeRadius);
1712   fprintf(file,"Seed E_{T}           %6.1f\n", fEtSeed);
1713   fprintf(file,"Min E_{T} of cell    %6.1f  ", fMinCellEt);
1714   fprintf(file,"Min E_{T} of jet     %6.1f\n", fMinJetEt);
1715   if(fMode) {
1716     fprintf(file,"====  Bg subtraction     ==== \n");
1717     fprintf(file,"BG subtraction  %6i  ", fMode);
1718     fprintf(file,"Min cone move   %6.2f\n", fMinMove);
1719     fprintf(file,"Max cone move   %6.2f  ", fMaxMove);
1720     fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1721   } else
1722   fprintf(file,"==== No Bg subtraction     ==== \n");
1723   //
1724   printf("BG file name is %s \n", fBGFileName.Data());
1725   if(file != stdout) fclose(file);
1726 }
1727
1728 void AliEMCALJetFinder::DrawLegos()
1729 {
1730 // Draw all legos
1731   if(!fC1) {
1732     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1733   }
1734   fC1->Clear();
1735   fC1->Divide(2,2);
1736   gStyle->SetOptStat(111111);
1737
1738   Int_t nent1, nent2, nent3, nent4;
1739   Double_t int1, int2, int3, int4;
1740   nent1 = (Int_t)fLego->GetEntries();
1741   int1  = fLego->Integral();
1742   fC1->cd(1);
1743   if(int1) fLego->Draw("lego");
1744
1745   nent2 = (Int_t)fhLegoTracks->GetEntries();
1746   int2  = fhLegoTracks->Integral();
1747   fC1->cd(2);
1748   if(int2) fhLegoTracks->Draw("lego");
1749
1750   nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1751   int3  = fhLegoEMCAL->Integral();
1752   fC1->cd(3);
1753   if(int3) fhLegoEMCAL->Draw("lego");
1754
1755   nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1756   int4  = fhLegoHadrCorr->Integral();
1757   fC1->cd(4);
1758   if(int4) fhLegoHadrCorr->Draw("lego");
1759
1760   // just for checking
1761   printf(" Integrals \n");
1762   printf("lego   %10.3f \ntrack  %10.3f \nhits   %10.3f \nHCorr   %10.3f\n--      %10.3f(must be 0)\n",
1763   int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1764 }
1765
1766 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1767 {
1768 // Get paramters from a file
1769   static TString tmp;
1770   if(strlen(dir)) tmp = dir;
1771   tmp += GetTitle();
1772   tmp += "_Par.txt";
1773   return tmp.Data();
1774 }
1775
1776 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1777 {
1778         // See FillFromTracks() - npart must be positive
1779     if (fTrackList) delete[] fTrackList;
1780     if (fPtT)       delete[] fPtT;
1781     if (fEtaT)      delete[] fEtaT;
1782     if (fPhiT)      delete[] fPhiT;
1783     if (fPdgT)      delete[] fPdgT;
1784
1785     if(npart>0) {
1786        fTrackList = new Int_t  [npart];
1787        fPtT       = new Float_t[npart];
1788        fEtaT      = new Float_t[npart];
1789        fPhiT      = new Float_t[npart];
1790        fPdgT      = new Int_t[npart];
1791     } else {
1792        printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1793     }
1794 }
1795
1796 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1797 {
1798 // Return quark info
1799   Int_t absPdg = TMath::Abs(pdg);
1800   if(absPdg<=6) return kTRUE; // quarks
1801   if(pdg == 21) return kTRUE; // gluon
1802   if(pdg == 92) return kTRUE; // string
1803
1804   // see p.51 of Pythia Manual
1805   // Not include diquarks with c and b quark - 4-mar-2002
1806   //                 ud_0,sd_0,su_0; dd_1,ud_1,uu_1;  sd_1,su_1,ss_1
1807   static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203,  3103,3203,3303};
1808   for(Int_t i=0; i<9; i++) if(absPdg == diquark[i])  return kTRUE; // diquarks
1809
1810   return kFALSE;
1811 }
1812
1813 void AliEMCALJetFinder::FindChargedJet()
1814 {
1815 //
1816 // Find jet from charged particle information only
1817 //
1818
1819 //
1820 //  Look for seeds
1821 //
1822     Int_t njets = 0;
1823     Int_t part  = 0;
1824     Int_t nseed = 0;
1825
1826 //
1827 //
1828     ResetJets();
1829
1830 //
1831     for (part = 0; part < fNt; part++) {
1832         if (!fTrackList[part]) continue;
1833         if (fPtT[part] > fEtSeed) nseed++;
1834     }
1835     printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1836     Int_t* iSeeds = new Int_t[nseed];
1837     nseed = 0;
1838
1839     for (part = 0; part < fNt; part++) {
1840         if (!fTrackList[part]) continue;
1841         if (fPtT[part] > fEtSeed) iSeeds[nseed++] =  part;
1842     }
1843
1844 //
1845 // Loop over seeds
1846 //
1847     Int_t seed = 0;
1848     Float_t pt;
1849
1850     while(1){
1851 //
1852 // Find seed with highest pt
1853 //
1854         Float_t ptmax = -1.;
1855         Int_t   index = -1;
1856         Int_t   jndex = -1;
1857         for (seed = 0; seed < nseed; seed++) {
1858             if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1859                 ptmax = pt;
1860                 index = seed;
1861             } // ptmax ?
1862         } // seeds
1863         if (ptmax < 0.) break;
1864         jndex = iSeeds[index];
1865 //
1866 // Remove from the list
1867         iSeeds[index] = -1;
1868         printf("\n Next Seed %d %f", jndex, ptmax);
1869 //
1870 // Find tracks in cone around seed
1871 //
1872         Float_t phiSeed = fPhiT[jndex];
1873         Float_t etaSeed = fEtaT[jndex];
1874         Float_t eT = 0.;
1875         Float_t pxJ = 0.;
1876         Float_t pyJ = 0.;
1877         Float_t pzJ = 0.;
1878
1879         for (part = 0; part < fNt; part++) {
1880             if (!fTrackList[part]) continue;
1881             Float_t deta = fEtaT[part] - etaSeed;
1882             Float_t dphi = fPhiT[part] - phiSeed;
1883             Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
1884             if (dR < fConeRadius) {
1885                 eT += fPtT[part];
1886                 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1887                 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1888                 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1889                 Float_t pz = fPtT[part] / TMath::Tan(theta);
1890                 pxJ += px;
1891                 pyJ += py;
1892                 pzJ += pz;
1893                 //
1894                 // if seed, remove it
1895                 //
1896                 for (seed = 0; seed < nseed; seed++) {
1897                     if (part == iSeeds[seed]) iSeeds[seed] = -1;
1898                 } // seed ?
1899             } // < cone radius
1900         } // particle loop
1901
1902 //
1903 //      Estimate of jet direction
1904         Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
1905         Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1906         Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
1907 //      Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1908
1909 //
1910 //      Sum up all energy
1911 //
1912         Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1913         Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1914         Int_t dIphi = Int_t(fConeRadius / fDphi);
1915         Int_t dIeta = Int_t(fConeRadius / fDeta);
1916         Int_t iPhi, iEta;
1917         Float_t sumE = 0;
1918         for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1919             for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1920                 if (iPhi < 0 || iEta < 0) continue;
1921                 Float_t dPhi = fPhiMin + iPhi * fDphi;
1922                 Float_t dEta = fEtaMin + iEta * fDeta;
1923                 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1924                 sumE += fLego->GetBinContent(iEta, iPhi);
1925             } // eta
1926         } // phi
1927 //
1928 //
1929 //
1930         fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1931         FindTracksInJetCone();
1932         printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1933         printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1934         printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1935     } // while(1)
1936     EMCALJETS.njet = njets;
1937     if (fWrite) WriteJets();
1938     fEvent++;
1939 }
1940 // 16-jan-2003 - just for convenience
1941 void AliEMCALJetFinder::Browse(TBrowser* b)
1942 {
1943 // Add to browser
1944    if(fHistsList)  b->Add((TObject*)fHistsList);
1945 }
1946
1947 Bool_t AliEMCALJetFinder::IsFolder() const
1948 {
1949 // Return folder status
1950   if(fHistsList) return kTRUE;
1951   else           return kFALSE;
1952 }
1953
1954 const Char_t* AliEMCALJetFinder::GetNameOfVariant()
1955 {
1956         // generate the literal string with info about jet finder
1957   Char_t name[200];
1958   sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
1959           fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
1960   TString nt(name);
1961   nt.ReplaceAll(" ","");
1962   if(fBGFileName.Length()) {
1963      Int_t i1 = fBGFileName.Index("kBackground");
1964      Int_t i2 = fBGFileName.Index("/0000") - 1;
1965      if(i1>=0 && i2>=0) {
1966         TString bg(fBGFileName(i1,i2-i1+1));
1967         nt += bg;
1968      }
1969   }
1970   printf("<I> Name of variant %s \n", nt.Data());
1971   return nt.Data();
1972 }