fParticleBuffer points to object which doesn't belong to AliStack, do not delete...
[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$
b7a73953 18Revision 1.16 2002/02/08 11:43:05 morsch
19SetOutputFileName(..) allows to specify an output file to which the
20reconstructed jets are written. If not set output goes to input file.
21Attention call Init() before processing.
22
08a3fb06 23Revision 1.15 2002/02/02 08:37:18 morsch
24Formula for rB corrected.
25
db5026bf 26Revision 1.14 2002/02/01 08:55:30 morsch
27Fill map with Et and pT.
28
f719d9b7 29Revision 1.13 2002/01/31 09:37:36 morsch
30Geometry parameters in constructor and call of SetCellSize()
31
be9787fe 32Revision 1.12 2002/01/23 13:40:23 morsch
33Fastidious debug print statement removed.
34
e4717047 35Revision 1.11 2002/01/22 17:25:47 morsch
36Some corrections for event mixing and bg event handling.
37
610f27c9 38Revision 1.10 2002/01/22 10:31:44 morsch
39Some correction for bg mixing.
40
46955412 41Revision 1.9 2002/01/21 16:28:42 morsch
42Correct sign of dphi.
43
f8f69f71 44Revision 1.8 2002/01/21 16:05:31 morsch
45- different phi-bin for hadron correction
46- provisions for background mixing.
47
ff114de3 48Revision 1.7 2002/01/21 15:47:18 morsch
49Bending radius correctly in cm.
50
b3feface 51Revision 1.6 2002/01/21 12:53:50 morsch
52authors
53
0c207dbc 54Revision 1.5 2002/01/21 12:47:47 morsch
55Possibility to include K0long and neutrons.
56
d53b7b0b 57Revision 1.4 2002/01/21 11:03:21 morsch
58Phi propagation introduced in FillFromTracks.
59
26dadf3a 60Revision 1.3 2002/01/18 05:07:56 morsch
61- hadronic correction
62- filling of digits
63- track selection upon EMCAL information
64
d594b22e 65*/
66
0c207dbc 67//*-- Authors: Andreas Morsch (CERN)
68//* J.L. Klay (LBL)
69//* Aleksei Pavlinov (WSU)
70
d594b22e 71// From root ...
471f69dc 72#include <TClonesArray.h>
73#include <TTree.h>
d594b22e 74#include <TBranchElement.h>
471f69dc 75#include <TFile.h>
76#include <TH2.h>
77#include <TAxis.h>
78#include <TParticle.h>
79#include <TParticlePDG.h>
d594b22e 80
81// From AliRoot ...
471f69dc 82#include "AliEMCALJetFinder.h"
d594b22e 83#include "AliEMCALFast.h"
471f69dc 84#include "AliEMCALGeometry.h"
85#include "AliEMCALHit.h"
d594b22e 86#include "AliEMCALDigit.h"
87#include "AliEMCALDigitizer.h"
88#include "AliEMCALHadronCorrection.h"
471f69dc 89#include "AliRun.h"
26dadf3a 90#include "AliMagF.h"
91#include "AliMagFCM.h"
471f69dc 92#include "AliEMCAL.h"
93#include "AliHeader.h"
d53b7b0b 94#include "AliPDG.h"
471f69dc 95
ff114de3 96// Interface to FORTRAN
97#include "Ecommon.h"
98
99
471f69dc 100ClassImp(AliEMCALJetFinder)
101
102//____________________________________________________________________________
103AliEMCALJetFinder::AliEMCALJetFinder()
104{
105// Default constructor
d594b22e 106 fJets = 0;
107 fNjets = 0;
108 fLego = 0;
ff114de3 109 fLegoB = 0;
d594b22e 110 fTrackList = 0;
111 fPtT = 0;
112 fEtaT = 0;
113 fPhiT = 0;
ff114de3 114 fPtB = 0;
115 fEtaB = 0;
116 fPhiB = 0;
b7a73953 117 fHCorrection = 0;
d594b22e 118 fHadronCorrector = 0;
b7a73953 119
08a3fb06 120 fOutFileName = 0;
121 fOutFile = 0;
122 fInFile = 0;
123 fEvent = 0;
b7a73953 124
125 SetParametersForBgSubtraction();
471f69dc 126}
127
128AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
129 : TTask(name, title)
130{
131// Constructor
132 fJets = new TClonesArray("AliEMCALJet",10000);
133 fNjets = 0;
d594b22e 134 for (Int_t i = 0; i < 30000; i++)
471f69dc 135 {
136 fEtCell[i] = 0.;
137 fEtaCell[i] = 0.;
138 fPhiCell[i] = 0.;
139 }
140 fLego = 0;
ff114de3 141 fTrackList = 0;
142 fPtT = 0;
143 fEtaT = 0;
144 fPhiT = 0;
145 fPtB = 0;
146 fEtaB = 0;
147 fPhiB = 0;
b7a73953 148 fHCorrection = 0;
149 fHadronCorrector = 0;
150 fBackground = 0;
08a3fb06 151 fOutFileName = 0;
152 fOutFile = 0;
153 fInFile = 0;
154 fEvent = 0;
0c207dbc 155//
d594b22e 156 SetPtCut();
157 SetMomentumSmearing();
158 SetEfficiencySim();
159 SetDebug();
160 SetHadronCorrection();
161 SetSamplingFraction();
d53b7b0b 162 SetIncludeK0andN();
471f69dc 163
b7a73953 164 SetParametersForBgSubtraction();
165}
471f69dc 166
b7a73953 167void AliEMCALJetFinder::SetParametersForBgSubtraction
168(Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
169{
170// see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
171// at WSU Linux cluster - 11-feb-2002
172// These parameters must be tuned more carefull !!!
173 SetMode(mode);
174 SetMinMove(minMove);
175 SetMaxMove(maxMove);
176 SetPrecBg(precBg);
177}
471f69dc 178
179//____________________________________________________________________________
180AliEMCALJetFinder::~AliEMCALJetFinder()
181{
182// Destructor
183//
184 if (fJets){
185 fJets->Delete();
186 delete fJets;
187 }
188}
189
190
191
3bc109a9 192
471f69dc 193#ifndef WIN32
194# define jet_finder_ua1 jet_finder_ua1_
195# define hf1 hf1_
196# define type_of_call
197
198#else
199# define jet_finder_ua1 J
200# define hf1 HF1
201# define type_of_call _stdcall
202#endif
203
204extern "C" void type_of_call
205jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
206 Float_t etc[30000], Float_t etac[30000],
207 Float_t phic[30000],
208 Float_t& min_move, Float_t& max_move, Int_t& mode,
209 Float_t& prec_bg, Int_t& ierror);
210
211extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
212
213
214
215
08a3fb06 216void AliEMCALJetFinder::Init()
217{
218//
219// Geometry and I/O initialization
220//
221//
222//
223// Get geometry parameters from EMCAL
224//
225//
226// Geometry
227 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
228 AliEMCALGeometry* geom =
229 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
230 fNbinEta = geom->GetNZ();
231 fNbinPhi = geom->GetNPhi();
232 fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
233 fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
234 fEtaMin = geom->GetArm1EtaMin();
235 fEtaMax = geom->GetArm1EtaMax();
236 fDphi = (fPhiMax-fPhiMin)/fNbinEta;
237 fDeta = (fEtaMax-fEtaMin)/fNbinEta;
238 fNtot = fNbinPhi*fNbinEta;
239//
240 SetCellSize(fDeta, fDphi);
241//
242// I/O
243 if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
244}
471f69dc 245
246void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
247 Float_t etac[30000], Float_t phic[30000],
248 Float_t min_move, Float_t max_move, Int_t mode,
249 Float_t prec_bg, Int_t ierror)
250{
251// Wrapper for fortran coded jet finder
252// Full argument list
253 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
254 min_move, max_move, mode, prec_bg, ierror);
255 // Write result to output
b7a73953 256 if(fWrite) WriteJets();
257 fEvent++;
471f69dc 258}
259
260void AliEMCALJetFinder::Find()
261{
262// Wrapper for fortran coded jet finder using member data for
263// argument list
264
b7a73953 265 Float_t min_move = fMinMove;
266 Float_t max_move = fMaxMove;
267 Int_t mode = fMode;
268 Float_t prec_bg = fPrecBg;
269 Int_t ierror;
270
271 ResetJets(); // 4-feb-2002 by PAI
471f69dc 272
273 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
274 min_move, max_move, mode, prec_bg, ierror);
b7a73953 275 fError = ierror;
471f69dc 276 // Write result to output
3bc109a9 277 Int_t njet = Njets();
d594b22e 278
3bc109a9 279 for (Int_t nj=0; nj<njet; nj++)
280 {
d594b22e 281
3bc109a9 282 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
283 JetPhiW(nj),
284 JetEtaW(nj));
285 }
b7a73953 286
287 if(fNt) FindTracksInJetCone();
288 if(fWrite) WriteJets();
08a3fb06 289 fEvent++;
471f69dc 290}
291
471f69dc 292Int_t AliEMCALJetFinder::Njets()
293{
294// Get number of reconstructed jets
295 return EMCALJETS.njet;
296}
297
298Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
299{
300// Get reconstructed jet energy
301 return EMCALJETS.etj[i];
302}
303
304Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
305{
306// Get reconstructed jet phi from leading particle
307 return EMCALJETS.phij[0][i];
308}
309
310Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
311{
312// Get reconstructed jet phi from weighting
313 return EMCALJETS.phij[1][i];
314}
315
316Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
317{
318// Get reconstructed jet eta from leading particles
319 return EMCALJETS.etaj[0][i];
320}
321
322
323Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
324{
325// Get reconstructed jet eta from weighting
326 return EMCALJETS.etaj[1][i];
327}
328
329void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
330{
331// Set grid cell size
332 EMCALCELLGEO.etaCellSize = eta;
333 EMCALCELLGEO.phiCellSize = phi;
334}
335
336void AliEMCALJetFinder::SetConeRadius(Float_t par)
337{
338// Set jet cone radius
339 EMCALJETPARAM.coneRad = par;
3bc109a9 340 fConeRadius = par;
471f69dc 341}
342
343void AliEMCALJetFinder::SetEtSeed(Float_t par)
344{
345// Set et cut for seeds
346 EMCALJETPARAM.etSeed = par;
d594b22e 347 fEtSeed = par;
471f69dc 348}
349
350void AliEMCALJetFinder::SetMinJetEt(Float_t par)
351{
352// Set minimum jet et
353 EMCALJETPARAM.ejMin = par;
d594b22e 354 fMinJetEt = par;
471f69dc 355}
356
357void AliEMCALJetFinder::SetMinCellEt(Float_t par)
358{
359// Set et cut per cell
360 EMCALJETPARAM.etMin = par;
d594b22e 361 fMinCellEt = par;
471f69dc 362}
363
3bc109a9 364void AliEMCALJetFinder::SetPtCut(Float_t par)
365{
366// Set pT cut on charged tracks
367 fPtCut = par;
368}
369
471f69dc 370
371void AliEMCALJetFinder::Test()
372{
373//
374// Test the finder call
375//
376 const Int_t nmax = 30000;
377 Int_t ncell = 10;
378 Int_t ncell_tot = 100;
379
380 Float_t etc[nmax];
381 Float_t etac[nmax];
382 Float_t phic[nmax];
383 Float_t min_move = 0;
384 Float_t max_move = 0;
385 Int_t mode = 0;
386 Float_t prec_bg = 0;
387 Int_t ierror = 0;
388
389
390 Find(ncell, ncell_tot, etc, etac, phic,
391 min_move, max_move, mode, prec_bg, ierror);
392
393}
394
395//
396// I/O
397//
398
399void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
400{
401 //
402 // Add a jet
403 //
404 TClonesArray &lrawcl = *fJets;
405 new(lrawcl[fNjets++]) AliEMCALJet(jet);
406}
407
408void AliEMCALJetFinder::ResetJets()
409{
410 //
411 // Reset Jet List
412 //
413 fJets->Clear();
414 fNjets = 0;
415}
416
417void AliEMCALJetFinder::WriteJets()
418{
419//
420// Add all jets to the list
421//
422 const Int_t kBufferSize = 4000;
08a3fb06 423 const char* file = 0;
424
471f69dc 425 Int_t njet = Njets();
08a3fb06 426
3bc109a9 427 for (Int_t nj = 0; nj < njet; nj++)
471f69dc 428 {
3bc109a9 429 AddJet(*fJetT[nj]);
430 delete fJetT[nj];
471f69dc 431 }
3bc109a9 432
08a3fb06 433// I/O
434 if (!fOutFileName) {
435//
436// output written to input file
437//
438 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
439 TTree* pK = gAlice->TreeK();
440 file = (pK->GetCurrentFile())->GetName();
441 if (fDebug > 1)
442 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
443 if (fJets && gAlice->TreeR()) {
444 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
445 "EMCALJets",
446 &fJets,
447 kBufferSize,
448 file);
449 }
450 Int_t nev = gAlice->GetHeader()->GetEvent();
451 gAlice->TreeR()->Fill();
452 char hname[30];
453 sprintf(hname,"TreeR%d", nev);
454 gAlice->TreeR()->Write(hname);
455 gAlice->TreeR()->Reset();
456 } else {
457//
458// Output written to user specified output file
459//
460 TTree* pK = gAlice->TreeK();
461 fInFile = pK->GetCurrentFile();
462
463 fOutFile->cd();
464 char hname[30];
465 sprintf(hname,"TreeR%d", fEvent);
466 TTree* treeJ = new TTree(hname, "EMCALJets");
467 treeJ->Branch("EMCALJets", &fJets, kBufferSize);
468 treeJ->Fill();
469 treeJ->Write(hname);
470 fInFile->cd();
471 }
471f69dc 472 ResetJets();
473}
474
475void AliEMCALJetFinder::BookLego()
476{
477//
478// Book histo for discretisation
479//
be9787fe 480
610f27c9 481//
482// Don't add histos to the current directory
483 TH2::AddDirectory(0);
471f69dc 484//
46955412 485// Signal map
471f69dc 486 fLego = new TH2F("legoH","eta-phi",
be9787fe 487 fNbinEta, fEtaMin, fEtaMax,
488 fNbinPhi, fPhiMin, fPhiMax);
46955412 489//
490// Background map
491 fLegoB = new TH2F("legoB","eta-phi",
be9787fe 492 fNbinEta, fEtaMin, fEtaMax,
493 fNbinPhi, fPhiMin, fPhiMax);
610f27c9 494
471f69dc 495}
496
497void AliEMCALJetFinder::DumpLego()
498{
499//
500// Dump lego histo into array
501//
502 fNcell = 0;
b7a73953 503 TAxis* Xaxis = fLego->GetXaxis();
504 TAxis* Yaxis = fLego->GetYaxis();
471f69dc 505 for (Int_t i = 1; i < fNbinEta; i++) {
506 for (Int_t j = 1; j < fNbinPhi; j++) {
507 Float_t e = fLego->GetBinContent(i,j);
46955412 508 if (e <=0.) continue;
509
471f69dc 510 Float_t eta = Xaxis->GetBinCenter(i);
511 Float_t phi = Yaxis->GetBinCenter(j);
512 fEtCell[fNcell] = e;
513 fEtaCell[fNcell] = eta;
514 fPhiCell[fNcell] = phi;
515 fNcell++;
516 } // phi
517 } // eta
518 fNcell--;
519}
520
521void AliEMCALJetFinder::ResetMap()
522{
523//
524// Reset eta-phi array
525
526 for (Int_t i=0; i<30000; i++)
527 {
528 fEtCell[i] = 0.;
529 fEtaCell[i] = 0.;
530 fPhiCell[i] = 0.;
531 }
532}
533
534
535void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
536{
537//
d594b22e 538// Fill Cells with track information
471f69dc 539//
b7a73953 540 if (fDebug >= 2)
541 printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
542
471f69dc 543 ResetMap();
544
545 if (!fLego) BookLego();
546// Reset
547 if (flag == 0) fLego->Reset();
548//
549// Access particle information
550 Int_t npart = (gAlice->GetHeader())->GetNprimary();
b7a73953 551 if (fDebug >= 2 || npart<=0) printf(" : npart %i\n", npart);
552
3bc109a9 553// Create track list
554//
555// 0: not selected
556// 1: selected for jet finding
557// 2: inside jet
558// ....
559 if (fTrackList) delete[] fTrackList;
560 if (fPtT) delete[] fPtT;
561 if (fEtaT) delete[] fEtaT;
562 if (fPhiT) delete[] fPhiT;
563
564 fTrackList = new Int_t [npart];
565 fPtT = new Float_t[npart];
566 fEtaT = new Float_t[npart];
567 fPhiT = new Float_t[npart];
ff114de3 568
3bc109a9 569 fNt = npart;
ff114de3 570 fNtS = 0;
b7a73953 571 Float_t chTmp=0.0; // charge of current particle
ff114de3 572
d594b22e 573 for (Int_t part = 2; part < npart; part++) {
471f69dc 574 TParticle *MPart = gAlice->Particle(part);
575 Int_t mpart = MPart->GetPdgCode();
576 Int_t child1 = MPart->GetFirstDaughter();
577 Float_t pT = MPart->Pt();
d594b22e 578 Float_t p = MPart->P();
471f69dc 579 Float_t phi = MPart->Phi();
3bc109a9 580 Float_t eta = MPart->Eta();
b7a73953 581 Float_t theta = MPart->Theta();
46955412 582
b7a73953 583 if (fDebug > 1) {
d594b22e 584 if (part == 6 || part == 7)
585 {
610f27c9 586 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
d594b22e 587 part-5, pT, eta, phi);
588 }
589
590// if (mpart == 21)
591
592// printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
593// part, mpart, pT, eta, phi);
594 }
595
3bc109a9 596
b7a73953 597 fTrackList[part] = 0;
598 fPtT[part] = pT; // must be change after correction for resolution !!!
3bc109a9 599 fEtaT[part] = eta;
600 fPhiT[part] = phi;
601
602 if (part < 2) continue;
b7a73953 603
604 // move to fLego->Fill because hadron correction must apply
605 // if particle hit to EMCAL - 11-feb-2002
606 // if (pT == 0 || pT < fPtCut) continue;
26dadf3a 607 TParticlePDG* pdgP = 0;
471f69dc 608// charged or neutral
b7a73953 609 pdgP = MPart->GetPDG();
610 chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
471f69dc 611 if (ich == 0) {
b7a73953 612 if (chTmp == 0) {
613 if (!fK0N) {
d53b7b0b 614 continue;
615 } else {
616 if (mpart != kNeutron &&
617 mpart != kNeutronBar &&
618 mpart != kK0Long) continue;
619 }
620 }
b7a73953 621 }
471f69dc 622// skip partons
623 if (TMath::Abs(mpart) <= 6 ||
624 mpart == 21 ||
625 mpart == 92) continue;
626// acceptance cut
627 if (TMath::Abs(eta) > 0.7) continue;
b7a73953 628// Initial phi may be out of acceptance but track may
629// hit two the acceptance - see variable curls
630// if (phi*180./TMath::Pi() > 120.) continue;
471f69dc 631// final state only
632 if (child1 >= 0 && child1 < npart) continue;
d594b22e 633//
634 if (fDebug > 1)
b7a73953 635 printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
636 part, mpart, child1, eta, phi, pT, chTmp);
26dadf3a 637//
26dadf3a 638//
d594b22e 639// Momentum smearing goes here ...
640//
b7a73953 641 Float_t pw;
642 if (fSmear && TMath::Abs(chTmp)) {
643 pw = AliEMCALFast::SmearMomentum(1,p);
644 // p changed - take into account when calculate pT,
645 // pz and so on .. ?
646 pT = (pw/p) * pT;
647 if(fDebug > 1) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
648 p = pw;
d594b22e 649 }
650//
651// Tracking Efficiency and TPC acceptance goes here ...
b7a73953 652 if (fEffic && TMath::Abs(chTmp)) {
d594b22e 653 Float_t eff = AliEMCALFast::Efficiency(1,p);
b7a73953 654 if (AliEMCALFast::RandomReject(eff)) {
655 if(fDebug > 1) printf(" reject due to unefficiency ");
656 continue;
657 }
d594b22e 658 }
659//
660// Correction of Hadronic Energy goes here ...
661//
b7a73953 662//
663// phi propagation for hadronic correction
664
665 Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
666 Float_t phiHC=0.0, dpH=0.0, dphi=0.0;
667 if(TMath::Abs(chTmp)) {
668 // hadr. correction only for charge particle
669 dphi = PropagatePhi(pT, chTmp, curls);
670 phiHC = phi + dphi;
671 if (fDebug >= 2) {
672 printf("\n Delta phi %f pT %f ", dphi, pT);
673 if (curls) printf("\n !! Track is curling");
674 }
ff114de3 675
b7a73953 676 if (fHCorrection && !curls) {
677 if (!fHadronCorrector)
678 Fatal("AliEMCALJetFinder",
d594b22e 679 "Hadronic energy correction required but not defined !");
b7a73953 680
681 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
682
683 if (fDebug >= 2) printf(" phi %f phiHC %f eTcorr %f\n",
684 phi, phiHC, -dpH*TMath::Sin(theta));
685 fLego->Fill(eta, phiHC, -dpH*TMath::Sin(theta));
686 }
687 }
d594b22e 688//
689// More to do ? Just think about it !
690//
b7a73953 691
692 if (phi*180./TMath::Pi() > 120.) continue;
693
694 if(TMath::Abs(chTmp) ) { // charge particle
695 if (pT > fPtCut && !curls) {
696 if (fDebug >= 2) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
697 eta , phi, pT);
698 fLego->Fill(eta, phi, pT);
699 fTrackList[part] = 1;
700 fNtS++;
701 }
702 } else if(ich == 0) {
703 // case of n, nbar and K0L
704 if (fDebug >= 2) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
705 eta , phi, pT);
706 fLego->Fill(eta, phi, pT);
707 fTrackList[part] = 1;
708 fNtS++;
709 }
710
471f69dc 711 } // primary loop
712 DumpLego();
713}
714
715void AliEMCALJetFinder::FillFromHits(Int_t flag)
716{
717//
d594b22e 718// Fill Cells with hit information
471f69dc 719//
720//
b7a73953 721 if (fDebug >= 2)
722 printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
723
471f69dc 724 ResetMap();
725
726 if (!fLego) BookLego();
ff114de3 727// Reset eta-phi map if needed
728 if (flag == 0) fLego->Reset();
729// Initialize from background event if available
471f69dc 730//
731// Access hit information
732 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
471f69dc 733
734 TTree *treeH = gAlice->TreeH();
735 Int_t ntracks = (Int_t) treeH->GetEntries();
736//
737// Loop over tracks
738//
739 Int_t nbytes = 0;
740
741
742 for (Int_t track=0; track<ntracks;track++) {
743 gAlice->ResetHits();
744 nbytes += treeH->GetEvent(track);
745//
746// Loop over hits
747//
748 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
749 mHit;
750 mHit=(AliEMCALHit*) pEMCAL->NextHit())
751 {
752 Float_t x = mHit->X(); // x-pos of hit
753 Float_t y = mHit->Y(); // y-pos
754 Float_t z = mHit->Z(); // z-pos
755 Float_t eloss = mHit->GetEnergy(); // deposited energy
d594b22e 756
757 Float_t r = TMath::Sqrt(x*x+y*y);
758 Float_t theta = TMath::ATan2(r,z);
759 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
760 Float_t phi = TMath::ATan2(y,x);
610f27c9 761
b7a73953 762 if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
610f27c9 763
f719d9b7 764 fLego->Fill(eta, phi, fSamplingF*eloss*TMath::Sin(theta));
471f69dc 765 } // Hit Loop
766 } // Track Loop
767 DumpLego();
768}
769
d594b22e 770void AliEMCALJetFinder::FillFromDigits(Int_t flag)
771{
772//
773// Fill Cells with digit information
774//
775//
776 ResetMap();
777
778 if (!fLego) BookLego();
779 if (flag == 0) fLego->Reset();
780 Int_t nbytes;
781
782
783//
784// Connect digits
785//
786 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
787 TTree *treeD = gAlice->TreeD();
788 TBranchElement* branchDg = (TBranchElement*)
789 treeD->GetBranch("EMCAL");
790
791 if (!branchDg) Fatal("AliEMCALJetFinder",
792 "Reading digits requested but no digits in file !");
793
794 branchDg->SetAddress(&digs);
795 Int_t nent = (Int_t) branchDg->GetEntries();
796//
797// Connect digitizer
798//
799 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
800 TBranchElement* branchDr = (TBranchElement*)
801 treeD->GetBranch("AliEMCALDigitizer");
802 branchDr->SetAddress(&digr);
803//
804//
805 nbytes += branchDg->GetEntry(0);
806 nbytes += branchDr->GetEntry(0);
807//
808// Get digitizer parameters
809 Float_t towerADCped = digr->GetTowerpedestal();
810 Float_t towerADCcha = digr->GetTowerchannel();
811 Float_t preshoADCped = digr->GetPreShopedestal();
812 Float_t preshoADCcha = digr->GetPreShochannel();
813
814 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
815 AliEMCALGeometry* geom =
816 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
817
818 if (fDebug) {
819 Int_t ndig = digs->GetEntries();
820 printf("\n Number of Digits: %d %d\n", ndig, nent);
821 printf("\n Parameters: %f %f %f %f\n",
822 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
823 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
824 }
825
826//
827// Loop over digits
828 AliEMCALDigit* sdg;
829 TIter next(digs);
830 while ((sdg = (AliEMCALDigit*)(next())))
831 {
832 Double_t pedestal;
833 Double_t channel;
834 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
835 {
836 pedestal = preshoADCped;
837 channel = preshoADCcha;
838 } else {
839 pedestal = towerADCped;
840 channel = towerADCcha;
841 }
842
843 Float_t eta = sdg->GetEta();
844 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
845 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
846
847 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
848 eta, phi, amp, sdg->GetAmp());
849
850 fLego->Fill(eta, phi, fSamplingF*amp);
851 } // digit loop
852//
853// Dump lego hist
854 DumpLego();
855}
856
857
858void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
859{
860//
861// Fill Cells with hit information
862//
863//
864 ResetMap();
865
866 if (!fLego) BookLego();
867// Reset
868 if (flag == 0) fLego->Reset();
869
ff114de3 870// Flag charged tracks passing through TPC which
871// are associated to EMCAL Hits
d594b22e 872 BuildTrackFlagTable();
873
874//
875// Access particle information
876 TTree *treeH = gAlice->TreeH();
877 Int_t ntracks = (Int_t) treeH->GetEntries();
878
ff114de3 879 if (fPtT) delete[] fPtT;
880 if (fEtaT) delete[] fEtaT;
881 if (fPhiT) delete[] fPhiT;
882
d594b22e 883 fPtT = new Float_t[ntracks];
884 fEtaT = new Float_t[ntracks];
885 fPhiT = new Float_t[ntracks];
d594b22e 886
ff114de3 887 fNt = ntracks;
888 fNtS = 0;
889
d594b22e 890 for (Int_t track = 0; track < ntracks; track++) {
891 TParticle *MPart = gAlice->Particle(track);
892 Float_t pT = MPart->Pt();
893 Float_t phi = MPart->Phi();
894 Float_t eta = MPart->Eta();
895
896 if(fTrackList[track]) {
897 fPtT[track] = pT;
898 fEtaT[track] = eta;
899 fPhiT[track] = phi;
900
901 if (track < 2) continue; //Colliding particles?
902 if (pT == 0 || pT < fPtCut) continue;
ff114de3 903 fNtS++;
d594b22e 904 fLego->Fill(eta, phi, pT);
905 }
906 } // track loop
907 DumpLego();
908}
909
b7a73953 910void AliEMCALJetFinder::FillFromPartons()
911{
912// function under construction - 13-feb-2002 PAI
913
914 if (fDebug >= 2)
915 printf("\n AliEMCALJetFinder::FillFromPartons()\n");
916 //
917
918 ResetMap();
919 if (!fLego) BookLego();
920 fLego->Reset();
921//
922// Access particle information
923 Int_t npart = (gAlice->GetHeader())->GetNprimary();
924 if (fDebug >= 2 || npart<=0)
925 printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
926 fNt = 0; // for FindTracksInJetCone
927 fNtS = 0;
928
929// Go through the partons
930 Int_t statusCode=0;
931 for (Int_t part = 8; part < npart; part++) {
932 TParticle *MPart = gAlice->Particle(part);
933 Int_t mpart = MPart->GetPdgCode();
934 // Int_t child1 = MPart->GetFirstDaughter();
935 Float_t pT = MPart->Pt();
936 // Float_t p = MPart->P();
937 Float_t phi = MPart->Phi();
938 Float_t eta = MPart->Eta();
939 // Float_t theta = MPart->Theta();
940 statusCode = MPart->GetStatusCode();
941
942// accept partons (21 - gluon, 92 - string)
943 if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
944 if (fDebug > 1 && pT>0.01)
945 printf("\n part:%5d mpart %5d status %5d eta %8.2f phi %8.2f pT %8.2f ",
946 part, mpart, statusCode, eta, phi, pT);
947 // if (fDebug >= 3) MPart->Print();
948// accept partons before fragmentation - p.57 in Pythia manual
949// if(statusCode != 1) continue;
950// acceptance cut
951 if (TMath::Abs(eta) > 0.7) continue;
952 if (phi*180./TMath::Pi() > 120.) continue;
953// final state only
954// if (child1 >= 0 && child1 < npart) continue;
955//
956//
957 fLego->Fill(eta, phi, pT);
958
959 } // primary loop
960 DumpLego();
961}
962
d594b22e 963void AliEMCALJetFinder::BuildTrackFlagTable() {
964
965// Method to generate a lookup table for TreeK particles
966// which are linked to hits in the EMCAL
967//
968// --Author: J.L. Klay
969//
d594b22e 970// Access hit information
971 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
972
973 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
974 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
975
976 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
977 fTrackList = new Int_t[nKTrks]; //before generating a new one
978
979 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
980 fTrackList[i] = 0;
981 }
982
983 TTree *treeH = gAlice->TreeH();
984 Int_t ntracks = (Int_t) treeH->GetEntries();
985//
986// Loop over tracks
987//
988 Int_t nbytes = 0;
989
990 for (Int_t track=0; track<ntracks;track++) {
991 gAlice->ResetHits();
992 nbytes += treeH->GetEvent(track);
993 if (pEMCAL) {
994
995//
996// Loop over hits
997//
998 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
999 mHit;
1000 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1001 {
1002 Int_t iTrk = mHit->Track(); // track number
1003 Int_t idprim = mHit->GetPrimary(); // primary particle
1004
1005 //Determine the origin point of this particle - it made a hit in the EMCAL
1006 TParticle *trkPart = gAlice->Particle(iTrk);
1007 TParticlePDG *trkPdg = trkPart->GetPDG();
1008 Int_t trkCode = trkPart->GetPdgCode();
1009 Double_t trkChg;
1010 if (trkCode < 10000) { //Big Ions cause problems for
1011 trkChg = trkPdg->Charge(); //this function. Since they aren't
1012 } else { //likely to make it very far, set
1013 trkChg = 0.0; //their charge to 0 for the Flag test
1014 }
1015 Float_t vxTrk = trkPart->Vx();
1016 Float_t vyTrk = trkPart->Vy();
1017 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1018 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1019
1020 //Loop through the ancestry of the EMCAL entrance particles
1021 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
1022 while (ancestor != -1) {
1023 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
1024 TParticlePDG *ancPdg = ancPart->GetPDG();
1025 Int_t ancCode = ancPart->GetPdgCode();
1026 Double_t ancChg;
1027 if (ancCode < 10000) {
1028 ancChg = ancPdg->Charge();
1029 } else {
1030 ancChg = 0.0;
1031 }
1032 Float_t vxAnc = ancPart->Vx();
1033 Float_t vyAnc = ancPart->Vy();
1034 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1035 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1036 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
1037 }
1038
1039 //Determine the origin point of the primary particle associated with the hit
1040 TParticle *primPart = gAlice->Particle(idprim);
1041 TParticlePDG *primPdg = primPart->GetPDG();
1042 Int_t primCode = primPart->GetPdgCode();
1043 Double_t primChg;
1044 if (primCode < 10000) {
1045 primChg = primPdg->Charge();
1046 } else {
1047 primChg = 0.0;
1048 }
1049 Float_t vxPrim = primPart->Vx();
1050 Float_t vyPrim = primPart->Vy();
1051 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1052 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1053 } // Hit Loop
1054 } //if (pEMCAL)
1055 } // Track Loop
1056}
1057
1058Int_t AliEMCALJetFinder
1059::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1060
ff114de3 1061 Int_t flag = 0;
1062 Int_t parton = 0;
d594b22e 1063 Int_t neutral = 0;
1064
1065 if (charge == 0) neutral = 1;
1066
1067 if (TMath::Abs(code) <= 6 ||
1068 code == 21 ||
1069 code == 92) parton = 1;
1070
1071 //It's not a parton, it's charged and it went through the TPC
1072 if (!parton && !neutral && radius <= 84.0) flag = 1;
1073
1074 return flag;
1075}
1076
ff114de3 1077
1078
1079void AliEMCALJetFinder::SaveBackgroundEvent()
1080{
1081// Saves the eta-phi lego and the tracklist
1082//
46955412 1083 if (fLegoB) fLegoB->Reset();
ff114de3 1084
46955412 1085 fLego->Copy(*fLegoB);
ff114de3 1086
1087 if (fPtB) delete[] fPtB;
1088 if (fEtaB) delete[] fEtaB;
1089 if (fPhiB) delete[] fPhiB;
1090 if (fTrackListB) delete[] fTrackListB;
1091
1092 fPtB = new Float_t[fNtS];
1093 fEtaB = new Float_t[fNtS];
1094 fPhiB = new Float_t[fNtS];
1095 fTrackListB = new Int_t [fNtS];
46955412 1096
1097 fNtB = 0;
1098
ff114de3 1099 for (Int_t i = 0; i < fNt; i++) {
1100 if (!fTrackList[i]) continue;
610f27c9 1101 fPtB [fNtB] = fPtT [i];
1102 fEtaB[fNtB] = fEtaT[i];
1103 fPhiB[fNtB] = fPhiT[i];
46955412 1104 fTrackListB[fNtB] = 1;
ff114de3 1105 fNtB++;
1106 }
610f27c9 1107 fBackground = 1;
ff114de3 1108}
1109
1110void AliEMCALJetFinder::InitFromBackground()
1111{
1112//
1113//
46955412 1114 if (fDebug) printf("\n Initializing from Background");
1115
1116 if (fLego) fLego->Reset();
1117 fLegoB->Copy(*fLego);
ff114de3 1118}
1119
ff114de3 1120
3bc109a9 1121void AliEMCALJetFinder::FindTracksInJetCone()
1122{
1123//
1124// Build list of tracks inside jet cone
1125//
1126// Loop over jets
1127 Int_t njet = Njets();
1128 for (Int_t nj = 0; nj < njet; nj++)
1129 {
1130 Float_t etaj = JetEtaW(nj);
1131 Float_t phij = JetPhiW(nj);
1132 Int_t nT = 0; // number of associated tracks
1133
ff114de3 1134// Loop over particles in current event
3bc109a9 1135 for (Int_t part = 0; part < fNt; part++) {
1136 if (!fTrackList[part]) continue;
1137 Float_t phi = fPhiT[part];
1138 Float_t eta = fEtaT[part];
1139 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1140 (phij-phi)*(phij-phi));
1141 if (dr < fConeRadius) {
1142 fTrackList[part] = nj+2;
1143 nT++;
1144 } // < ConeRadius ?
1145 } // particle loop
d594b22e 1146
ff114de3 1147// Same for background event if available
1148 Int_t nTB = 0;
1149 if (fBackground) {
1150 for (Int_t part = 0; part < fNtB; part++) {
1151 Float_t phi = fPhiB[part];
1152 Float_t eta = fEtaB[part];
1153 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1154 (phij-phi)*(phij-phi));
610f27c9 1155 fTrackListB[part] = 1;
1156
ff114de3 1157 if (dr < fConeRadius) {
610f27c9 1158 fTrackListB[part] = nj+2;
ff114de3 1159 nTB++;
1160 } // < ConeRadius ?
1161 } // particle loop
1162 } // Background available ?
1163
1164 Int_t nT0 = nT + nTB;
1165
1166 if (nT0 > 50) nT0 = 50;
d594b22e 1167
610f27c9 1168 Float_t* ptT = new Float_t[nT0];
1169 Float_t* etaT = new Float_t[nT0];
1170 Float_t* phiT = new Float_t[nT0];
d594b22e 1171 Int_t iT = 0;
1172 Int_t j;
1173
3bc109a9 1174 for (Int_t part = 0; part < fNt; part++) {
1175 if (fTrackList[part] == nj+2) {
d594b22e 1176 Int_t index = 0;
1177 for (j=iT-1; j>=0; j--) {
1178 if (fPtT[part] > ptT[j]) {
1179 index = j+1;
1180 break;
1181 }
1182 }
1183 for (j=iT-1; j>=index; j--) {
1184 ptT [j+1] = ptT [j];
1185 etaT[j+1] = etaT[j];
1186 phiT[j+1] = phiT[j];
1187 }
1188 ptT [index] = fPtT [part];
1189 etaT[index] = fEtaT[part];
1190 phiT[index] = fPhiT[part];
ff114de3 1191 iT++;
3bc109a9 1192 } // particle associated
ff114de3 1193 if (iT > nT0) break;
3bc109a9 1194 } // particle loop
d594b22e 1195
ff114de3 1196 if (fBackground) {
1197 for (Int_t part = 0; part < fNtB; part++) {
610f27c9 1198 if (fTrackListB[part] == nj+2) {
ff114de3 1199 Int_t index = 0;
1200 for (j=iT-1; j>=0; j--) {
1201 if (fPtB[part] > ptT[j]) {
1202 index = j+1;
610f27c9 1203
ff114de3 1204 break;
1205 }
1206 }
1207 for (j=iT-1; j>=index; j--) {
1208 ptT [j+1] = ptT [j];
1209 etaT[j+1] = etaT[j];
1210 phiT[j+1] = phiT[j];
1211 }
1212 ptT [index] = fPtB [part];
1213 etaT[index] = fEtaB[part];
1214 phiT[index] = fPhiB[part];
1215 iT++;
1216 } // particle associated
1217 if (iT > nT0) break;
1218 } // particle loop
1219 } // Background available ?
1220
610f27c9 1221 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
3bc109a9 1222 delete[] ptT;
1223 delete[] etaT;
1224 delete[] phiT;
1225 } // jet loop loop
1226}
1227
26dadf3a 1228Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1229{
1230// Propagates phi angle to EMCAL radius
1231//
b7a73953 1232 static Float_t b = 0.0, rEMCAL = -1.0;
1233 if(rEMCAL<0) {
1234// Get field in kGS
1235 b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
1236// Get EMCAL radius in cm
1237 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1238 printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
1239 }
26dadf3a 1240 Float_t dPhi = 0.;
26dadf3a 1241//
1242//
1243// bending radies
db5026bf 1244// pt [Gev]
1245// B [kG]
1246//
b7a73953 1247 Float_t rB = 3335.6 * pt / b; // [cm] (case of |charge|=1)
26dadf3a 1248//
1249// check if particle is curling below EMCAL
1250 if (2.*rB < rEMCAL) {
1251 curls = kTRUE;
1252 return dPhi;
1253 }
1254//
1255// if not calculate delta phi
1256 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1257 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
f8f69f71 1258 dPhi = -TMath::Sign(dPhi, charge);
26dadf3a 1259//
1260 return dPhi;
26dadf3a 1261}
1262
1263
471f69dc 1264void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1265{
1266// dummy for hbook calls
1267 ;
1268}
1269
1270
d594b22e 1271
1272
1273
1274
1275
1276