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