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