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