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