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