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