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