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