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