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