EMCAL/AliEMCALJetFinder.cxx

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