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