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