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