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