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