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