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