1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.11 2002/01/22 17:25:47 morsch
19 Some corrections for event mixing and bg event handling.
21 Revision 1.10 2002/01/22 10:31:44 morsch
22 Some correction for bg mixing.
24 Revision 1.9 2002/01/21 16:28:42 morsch
27 Revision 1.8 2002/01/21 16:05:31 morsch
28 - different phi-bin for hadron correction
29 - provisions for background mixing.
31 Revision 1.7 2002/01/21 15:47:18 morsch
32 Bending radius correctly in cm.
34 Revision 1.6 2002/01/21 12:53:50 morsch
37 Revision 1.5 2002/01/21 12:47:47 morsch
38 Possibility to include K0long and neutrons.
40 Revision 1.4 2002/01/21 11:03:21 morsch
41 Phi propagation introduced in FillFromTracks.
43 Revision 1.3 2002/01/18 05:07:56 morsch
46 - track selection upon EMCAL information
50 //*-- Authors: Andreas Morsch (CERN)
52 //* Aleksei Pavlinov (WSU)
55 #include <TClonesArray.h>
57 #include <TBranchElement.h>
61 #include <TParticle.h>
62 #include <TParticlePDG.h>
65 #include "AliEMCALJetFinder.h"
66 #include "AliEMCALFast.h"
67 #include "AliEMCALGeometry.h"
68 #include "AliEMCALHit.h"
69 #include "AliEMCALDigit.h"
70 #include "AliEMCALDigitizer.h"
71 #include "AliEMCALHadronCorrection.h"
74 #include "AliMagFCM.h"
76 #include "AliHeader.h"
79 // Interface to FORTRAN
83 ClassImp(AliEMCALJetFinder)
85 //____________________________________________________________________________
86 AliEMCALJetFinder::AliEMCALJetFinder()
88 // Default constructor
100 fHadronCorrector = 0;
103 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
107 fJets = new TClonesArray("AliEMCALJet",10000);
109 for (Int_t i = 0; i < 30000; i++)
123 fHadronCorrector = 0;
127 SetMomentumSmearing();
130 SetHadronCorrection();
131 SetSamplingFraction();
138 //____________________________________________________________________________
139 AliEMCALJetFinder::~AliEMCALJetFinder()
153 # define jet_finder_ua1 jet_finder_ua1_
155 # define type_of_call
158 # define jet_finder_ua1 J
160 # define type_of_call _stdcall
163 extern "C" void type_of_call
164 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
165 Float_t etc[30000], Float_t etac[30000],
167 Float_t& min_move, Float_t& max_move, Int_t& mode,
168 Float_t& prec_bg, Int_t& ierror);
170 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
176 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
177 Float_t etac[30000], Float_t phic[30000],
178 Float_t min_move, Float_t max_move, Int_t mode,
179 Float_t prec_bg, Int_t ierror)
181 // Wrapper for fortran coded jet finder
182 // Full argument list
183 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
184 min_move, max_move, mode, prec_bg, ierror);
185 // Write result to output
189 void AliEMCALJetFinder::Find()
191 // Wrapper for fortran coded jet finder using member data for
194 Float_t min_move = 0;
195 Float_t max_move = 0;
197 Float_t prec_bg = 0.;
200 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
201 min_move, max_move, mode, prec_bg, ierror);
202 // Write result to output
203 Int_t njet = Njets();
205 for (Int_t nj=0; nj<njet; nj++)
209 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
213 FindTracksInJetCone();
218 Int_t AliEMCALJetFinder::Njets()
220 // Get number of reconstructed jets
221 return EMCALJETS.njet;
224 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
226 // Get reconstructed jet energy
227 return EMCALJETS.etj[i];
230 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
232 // Get reconstructed jet phi from leading particle
233 return EMCALJETS.phij[0][i];
236 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
238 // Get reconstructed jet phi from weighting
239 return EMCALJETS.phij[1][i];
242 Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
244 // Get reconstructed jet eta from leading particles
245 return EMCALJETS.etaj[0][i];
249 Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
251 // Get reconstructed jet eta from weighting
252 return EMCALJETS.etaj[1][i];
255 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
257 // Set grid cell size
258 EMCALCELLGEO.etaCellSize = eta;
259 EMCALCELLGEO.phiCellSize = phi;
264 void AliEMCALJetFinder::SetConeRadius(Float_t par)
266 // Set jet cone radius
267 EMCALJETPARAM.coneRad = par;
271 void AliEMCALJetFinder::SetEtSeed(Float_t par)
273 // Set et cut for seeds
274 EMCALJETPARAM.etSeed = par;
278 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
280 // Set minimum jet et
281 EMCALJETPARAM.ejMin = par;
285 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
287 // Set et cut per cell
288 EMCALJETPARAM.etMin = par;
292 void AliEMCALJetFinder::SetPtCut(Float_t par)
294 // Set pT cut on charged tracks
299 void AliEMCALJetFinder::Test()
302 // Test the finder call
304 const Int_t nmax = 30000;
306 Int_t ncell_tot = 100;
311 Float_t min_move = 0;
312 Float_t max_move = 0;
318 Find(ncell, ncell_tot, etc, etac, phic,
319 min_move, max_move, mode, prec_bg, ierror);
327 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
332 TClonesArray &lrawcl = *fJets;
333 new(lrawcl[fNjets++]) AliEMCALJet(jet);
336 void AliEMCALJetFinder::ResetJets()
345 void AliEMCALJetFinder::WriteJets()
348 // Add all jets to the list
350 const Int_t kBufferSize = 4000;
351 TTree *pK = gAlice->TreeK();
352 const char* file = (pK->GetCurrentFile())->GetName();
354 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
357 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
359 if (fJets && gAlice->TreeR()) {
360 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
366 Int_t njet = Njets();
367 for (Int_t nj = 0; nj < njet; nj++)
373 Int_t nev = gAlice->GetHeader()->GetEvent();
374 gAlice->TreeR()->Fill();
376 sprintf(hname,"TreeR%d", nev);
377 gAlice->TreeR()->Write(hname);
378 gAlice->TreeR()->Reset();
382 void AliEMCALJetFinder::BookLego()
385 // Book histo for discretisation
388 // Get geometry parameters from
389 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
390 AliEMCALGeometry* geom =
391 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
392 fNbinEta = geom->GetNZ();
393 fNbinPhi = geom->GetNPhi();
394 const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
395 const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
396 const Float_t etaMin = geom->GetArm1EtaMin();
397 const Float_t etaMax = geom->GetArm1EtaMax();
398 fDphi = (phiMax-phiMin)/fNbinEta;
399 fDeta = 1.4/fNbinEta;
400 fNtot = fNbinPhi*fNbinEta;
402 // Don't add histos to the current directory
403 TH2::AddDirectory(0);
406 fLego = new TH2F("legoH","eta-phi",
407 fNbinEta, etaMin, etaMax,
408 fNbinPhi, phiMin, phiMax);
411 fLegoB = new TH2F("legoB","eta-phi",
412 fNbinEta, etaMin, etaMax,
413 fNbinPhi, phiMin, phiMax);
417 void AliEMCALJetFinder::DumpLego()
420 // Dump lego histo into array
423 for (Int_t i = 1; i < fNbinEta; i++) {
424 for (Int_t j = 1; j < fNbinPhi; j++) {
425 Float_t e = fLego->GetBinContent(i,j);
426 if (e <=0.) continue;
428 TAxis* Xaxis = fLego->GetXaxis();
429 TAxis* Yaxis = fLego->GetYaxis();
430 Float_t eta = Xaxis->GetBinCenter(i);
431 Float_t phi = Yaxis->GetBinCenter(j);
433 fEtaCell[fNcell] = eta;
434 fPhiCell[fNcell] = phi;
441 void AliEMCALJetFinder::ResetMap()
444 // Reset eta-phi array
446 for (Int_t i=0; i<30000; i++)
455 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
458 // Fill Cells with track information
463 if (!fLego) BookLego();
465 if (flag == 0) fLego->Reset();
467 // Access particle information
468 Int_t npart = (gAlice->GetHeader())->GetNprimary();
472 // 1: selected for jet finding
475 if (fTrackList) delete[] fTrackList;
476 if (fPtT) delete[] fPtT;
477 if (fEtaT) delete[] fEtaT;
478 if (fPhiT) delete[] fPhiT;
480 fTrackList = new Int_t [npart];
481 fPtT = new Float_t[npart];
482 fEtaT = new Float_t[npart];
483 fPhiT = new Float_t[npart];
488 for (Int_t part = 2; part < npart; part++) {
489 TParticle *MPart = gAlice->Particle(part);
490 Int_t mpart = MPart->GetPdgCode();
491 Int_t child1 = MPart->GetFirstDaughter();
492 Float_t pT = MPart->Pt();
493 Float_t p = MPart->P();
494 Float_t phi = MPart->Phi();
495 Float_t eta = MPart->Eta();
498 if (part == 6 || part == 7)
500 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
501 part-5, pT, eta, phi);
506 // printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
507 // part, mpart, pT, eta, phi);
511 fTrackList[part] = 0;
516 if (part < 2) continue;
517 if (pT == 0 || pT < fPtCut) continue;
518 TParticlePDG* pdgP = 0;
519 // charged or neutral
521 pdgP = MPart->GetPDG();
522 if (pdgP->Charge() == 0) {
526 if (mpart != kNeutron &&
527 mpart != kNeutronBar &&
528 mpart != kK0Long) continue;
533 if (TMath::Abs(mpart) <= 6 ||
535 mpart == 92) continue;
537 if (TMath::Abs(eta) > 0.7) continue;
538 if (phi*180./TMath::Pi() > 120.) continue;
540 if (child1 >= 0 && child1 < npart) continue;
543 printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
544 part, mpart, child1, eta, phi, pT);
547 // phi propagation for hadronic correction
549 Bool_t curls = kFALSE;
550 Float_t dphi = PropagatePhi(pT, pdgP->Charge(), curls);
551 if (fDebug > 1) printf("\n Delta phi %f pT%f", dphi, pT);
552 if (curls && fDebug > 1) printf("\n Track is curling %f", pT);
553 Float_t phiHC = phi + dphi;
555 // Momentum smearing goes here ...
558 p = AliEMCALFast::SmearMomentum(1,p);
561 // Tracking Efficiency and TPC acceptance goes here ...
563 Float_t eff = AliEMCALFast::Efficiency(1,p);
564 if (AliEMCALFast::RandomReject(eff)) continue;
567 // Correction of Hadronic Energy goes here ...
572 if (!fHadronCorrector)
573 Fatal("AliEMCALJetFinder",
574 "Hadronic energy correction required but not defined !");
575 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
578 // More to do ? Just think about it !
580 fTrackList[part] = 1;
584 fLego->Fill(eta, phi, p);
585 if (fHCorrection && !curls) fLego->Fill(eta, phiHC, -dpH);
590 void AliEMCALJetFinder::FillFromHits(Int_t flag)
593 // Fill Cells with hit information
598 if (!fLego) BookLego();
599 // Reset eta-phi map if needed
600 if (flag == 0) fLego->Reset();
601 // Initialize from background event if available
603 // Access hit information
604 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
606 TTree *treeH = gAlice->TreeH();
607 Int_t ntracks = (Int_t) treeH->GetEntries();
614 for (Int_t track=0; track<ntracks;track++) {
616 nbytes += treeH->GetEvent(track);
620 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
622 mHit=(AliEMCALHit*) pEMCAL->NextHit())
624 Float_t x = mHit->X(); // x-pos of hit
625 Float_t y = mHit->Y(); // y-pos
626 Float_t z = mHit->Z(); // z-pos
627 Float_t eloss = mHit->GetEnergy(); // deposited energy
629 Float_t r = TMath::Sqrt(x*x+y*y);
630 Float_t theta = TMath::ATan2(r,z);
631 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
632 Float_t phi = TMath::ATan2(y,x);
634 if (fDebug > 1) printf("\n Hit %f %f %f %f", x, y, z, eloss);
636 fLego->Fill(eta, phi, fSamplingF*eloss);
642 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
645 // Fill Cells with digit information
650 if (!fLego) BookLego();
651 if (flag == 0) fLego->Reset();
658 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
659 TTree *treeD = gAlice->TreeD();
660 TBranchElement* branchDg = (TBranchElement*)
661 treeD->GetBranch("EMCAL");
663 if (!branchDg) Fatal("AliEMCALJetFinder",
664 "Reading digits requested but no digits in file !");
666 branchDg->SetAddress(&digs);
667 Int_t nent = (Int_t) branchDg->GetEntries();
671 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
672 TBranchElement* branchDr = (TBranchElement*)
673 treeD->GetBranch("AliEMCALDigitizer");
674 branchDr->SetAddress(&digr);
677 nbytes += branchDg->GetEntry(0);
678 nbytes += branchDr->GetEntry(0);
680 // Get digitizer parameters
681 Float_t towerADCped = digr->GetTowerpedestal();
682 Float_t towerADCcha = digr->GetTowerchannel();
683 Float_t preshoADCped = digr->GetPreShopedestal();
684 Float_t preshoADCcha = digr->GetPreShochannel();
686 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
687 AliEMCALGeometry* geom =
688 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
691 Int_t ndig = digs->GetEntries();
692 printf("\n Number of Digits: %d %d\n", ndig, nent);
693 printf("\n Parameters: %f %f %f %f\n",
694 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
695 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
702 while ((sdg = (AliEMCALDigit*)(next())))
706 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
708 pedestal = preshoADCped;
709 channel = preshoADCcha;
711 pedestal = towerADCped;
712 channel = towerADCcha;
715 Float_t eta = sdg->GetEta();
716 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
717 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
719 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
720 eta, phi, amp, sdg->GetAmp());
722 fLego->Fill(eta, phi, fSamplingF*amp);
730 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
733 // Fill Cells with hit information
738 if (!fLego) BookLego();
740 if (flag == 0) fLego->Reset();
742 // Flag charged tracks passing through TPC which
743 // are associated to EMCAL Hits
744 BuildTrackFlagTable();
747 // Access particle information
748 TTree *treeH = gAlice->TreeH();
749 Int_t ntracks = (Int_t) treeH->GetEntries();
751 if (fPtT) delete[] fPtT;
752 if (fEtaT) delete[] fEtaT;
753 if (fPhiT) delete[] fPhiT;
755 fPtT = new Float_t[ntracks];
756 fEtaT = new Float_t[ntracks];
757 fPhiT = new Float_t[ntracks];
762 for (Int_t track = 0; track < ntracks; track++) {
763 TParticle *MPart = gAlice->Particle(track);
764 Float_t pT = MPart->Pt();
765 Float_t phi = MPart->Phi();
766 Float_t eta = MPart->Eta();
768 if(fTrackList[track]) {
773 if (track < 2) continue; //Colliding particles?
774 if (pT == 0 || pT < fPtCut) continue;
776 fLego->Fill(eta, phi, pT);
782 void AliEMCALJetFinder::BuildTrackFlagTable() {
784 // Method to generate a lookup table for TreeK particles
785 // which are linked to hits in the EMCAL
787 // --Author: J.L. Klay
789 // Access hit information
790 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
792 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
793 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
795 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
796 fTrackList = new Int_t[nKTrks]; //before generating a new one
798 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
802 TTree *treeH = gAlice->TreeH();
803 Int_t ntracks = (Int_t) treeH->GetEntries();
809 for (Int_t track=0; track<ntracks;track++) {
811 nbytes += treeH->GetEvent(track);
817 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
819 mHit=(AliEMCALHit*) pEMCAL->NextHit())
821 Int_t iTrk = mHit->Track(); // track number
822 Int_t idprim = mHit->GetPrimary(); // primary particle
824 //Determine the origin point of this particle - it made a hit in the EMCAL
825 TParticle *trkPart = gAlice->Particle(iTrk);
826 TParticlePDG *trkPdg = trkPart->GetPDG();
827 Int_t trkCode = trkPart->GetPdgCode();
829 if (trkCode < 10000) { //Big Ions cause problems for
830 trkChg = trkPdg->Charge(); //this function. Since they aren't
831 } else { //likely to make it very far, set
832 trkChg = 0.0; //their charge to 0 for the Flag test
834 Float_t vxTrk = trkPart->Vx();
835 Float_t vyTrk = trkPart->Vy();
836 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
837 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
839 //Loop through the ancestry of the EMCAL entrance particles
840 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
841 while (ancestor != -1) {
842 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
843 TParticlePDG *ancPdg = ancPart->GetPDG();
844 Int_t ancCode = ancPart->GetPdgCode();
846 if (ancCode < 10000) {
847 ancChg = ancPdg->Charge();
851 Float_t vxAnc = ancPart->Vx();
852 Float_t vyAnc = ancPart->Vy();
853 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
854 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
855 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
858 //Determine the origin point of the primary particle associated with the hit
859 TParticle *primPart = gAlice->Particle(idprim);
860 TParticlePDG *primPdg = primPart->GetPDG();
861 Int_t primCode = primPart->GetPdgCode();
863 if (primCode < 10000) {
864 primChg = primPdg->Charge();
868 Float_t vxPrim = primPart->Vx();
869 Float_t vyPrim = primPart->Vy();
870 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
871 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
877 Int_t AliEMCALJetFinder
878 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
884 if (charge == 0) neutral = 1;
886 if (TMath::Abs(code) <= 6 ||
888 code == 92) parton = 1;
890 //It's not a parton, it's charged and it went through the TPC
891 if (!parton && !neutral && radius <= 84.0) flag = 1;
898 void AliEMCALJetFinder::SaveBackgroundEvent()
900 // Saves the eta-phi lego and the tracklist
902 if (fLegoB) fLegoB->Reset();
904 fLego->Copy(*fLegoB);
906 if (fPtB) delete[] fPtB;
907 if (fEtaB) delete[] fEtaB;
908 if (fPhiB) delete[] fPhiB;
909 if (fTrackListB) delete[] fTrackListB;
911 fPtB = new Float_t[fNtS];
912 fEtaB = new Float_t[fNtS];
913 fPhiB = new Float_t[fNtS];
914 fTrackListB = new Int_t [fNtS];
918 for (Int_t i = 0; i < fNt; i++) {
919 if (!fTrackList[i]) continue;
920 fPtB [fNtB] = fPtT [i];
921 fEtaB[fNtB] = fEtaT[i];
922 fPhiB[fNtB] = fPhiT[i];
923 fTrackListB[fNtB] = 1;
929 void AliEMCALJetFinder::InitFromBackground()
933 if (fDebug) printf("\n Initializing from Background");
935 if (fLego) fLego->Reset();
936 fLegoB->Copy(*fLego);
941 void AliEMCALJetFinder::FindTracksInJetCone()
944 // Build list of tracks inside jet cone
947 Int_t njet = Njets();
948 for (Int_t nj = 0; nj < njet; nj++)
950 Float_t etaj = JetEtaW(nj);
951 Float_t phij = JetPhiW(nj);
952 Int_t nT = 0; // number of associated tracks
954 // Loop over particles in current event
955 for (Int_t part = 0; part < fNt; part++) {
956 if (!fTrackList[part]) continue;
957 Float_t phi = fPhiT[part];
958 Float_t eta = fEtaT[part];
959 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
960 (phij-phi)*(phij-phi));
961 if (dr < fConeRadius) {
962 fTrackList[part] = nj+2;
967 // Same for background event if available
970 for (Int_t part = 0; part < fNtB; part++) {
971 Float_t phi = fPhiB[part];
972 Float_t eta = fEtaB[part];
973 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
974 (phij-phi)*(phij-phi));
975 fTrackListB[part] = 1;
977 if (dr < fConeRadius) {
978 fTrackListB[part] = nj+2;
982 } // Background available ?
984 Int_t nT0 = nT + nTB;
986 if (nT0 > 50) nT0 = 50;
988 Float_t* ptT = new Float_t[nT0];
989 Float_t* etaT = new Float_t[nT0];
990 Float_t* phiT = new Float_t[nT0];
994 for (Int_t part = 0; part < fNt; part++) {
995 if (fTrackList[part] == nj+2) {
997 for (j=iT-1; j>=0; j--) {
998 if (fPtT[part] > ptT[j]) {
1003 for (j=iT-1; j>=index; j--) {
1004 ptT [j+1] = ptT [j];
1005 etaT[j+1] = etaT[j];
1006 phiT[j+1] = phiT[j];
1008 ptT [index] = fPtT [part];
1009 etaT[index] = fEtaT[part];
1010 phiT[index] = fPhiT[part];
1012 } // particle associated
1013 if (iT > nT0) break;
1017 for (Int_t part = 0; part < fNtB; part++) {
1018 if (fTrackListB[part] == nj+2) {
1020 for (j=iT-1; j>=0; j--) {
1021 if (fPtB[part] > ptT[j]) {
1027 for (j=iT-1; j>=index; j--) {
1028 ptT [j+1] = ptT [j];
1029 etaT[j+1] = etaT[j];
1030 phiT[j+1] = phiT[j];
1032 ptT [index] = fPtB [part];
1033 etaT[index] = fEtaB[part];
1034 phiT[index] = fPhiB[part];
1036 } // particle associated
1037 if (iT > nT0) break;
1039 } // Background available ?
1041 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
1048 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1050 // Propagates phi angle to EMCAL radius
1054 Float_t b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
1056 Float_t rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1060 Float_t rB = 333.56 * pt / b; // [cm]
1062 // check if particle is curling below EMCAL
1063 if (2.*rB < rEMCAL) {
1068 // if not calculate delta phi
1069 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1070 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1071 dPhi = -TMath::Sign(dPhi, charge);
1078 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1080 // dummy for hbook calls