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