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