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