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