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