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