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