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.10 2002/01/22 10:31:44 morsch
19 Some correction for bg mixing.
21 Revision 1.9 2002/01/21 16:28:42 morsch
24 Revision 1.8 2002/01/21 16:05:31 morsch
25 - different phi-bin for hadron correction
26 - provisions for background mixing.
28 Revision 1.7 2002/01/21 15:47:18 morsch
29 Bending radius correctly in cm.
31 Revision 1.6 2002/01/21 12:53:50 morsch
34 Revision 1.5 2002/01/21 12:47:47 morsch
35 Possibility to include K0long and neutrons.
37 Revision 1.4 2002/01/21 11:03:21 morsch
38 Phi propagation introduced in FillFromTracks.
40 Revision 1.3 2002/01/18 05:07:56 morsch
43 - track selection upon EMCAL information
47 //*-- Authors: Andreas Morsch (CERN)
49 //* Aleksei Pavlinov (WSU)
52 #include <TClonesArray.h>
54 #include <TBranchElement.h>
58 #include <TParticle.h>
59 #include <TParticlePDG.h>
62 #include "AliEMCALJetFinder.h"
63 #include "AliEMCALFast.h"
64 #include "AliEMCALGeometry.h"
65 #include "AliEMCALHit.h"
66 #include "AliEMCALDigit.h"
67 #include "AliEMCALDigitizer.h"
68 #include "AliEMCALHadronCorrection.h"
71 #include "AliMagFCM.h"
73 #include "AliHeader.h"
76 // Interface to FORTRAN
80 ClassImp(AliEMCALJetFinder)
82 //____________________________________________________________________________
83 AliEMCALJetFinder::AliEMCALJetFinder()
85 // Default constructor
100 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
104 fJets = new TClonesArray("AliEMCALJet",10000);
106 for (Int_t i = 0; i < 30000; i++)
120 fHadronCorrector = 0;
124 SetMomentumSmearing();
127 SetHadronCorrection();
128 SetSamplingFraction();
135 //____________________________________________________________________________
136 AliEMCALJetFinder::~AliEMCALJetFinder()
150 # define jet_finder_ua1 jet_finder_ua1_
152 # define type_of_call
155 # define jet_finder_ua1 J
157 # define type_of_call _stdcall
160 extern "C" void type_of_call
161 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
162 Float_t etc[30000], Float_t etac[30000],
164 Float_t& min_move, Float_t& max_move, Int_t& mode,
165 Float_t& prec_bg, Int_t& ierror);
167 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
173 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
174 Float_t etac[30000], Float_t phic[30000],
175 Float_t min_move, Float_t max_move, Int_t mode,
176 Float_t prec_bg, Int_t ierror)
178 // Wrapper for fortran coded jet finder
179 // Full argument list
180 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
181 min_move, max_move, mode, prec_bg, ierror);
182 // Write result to output
186 void AliEMCALJetFinder::Find()
188 // Wrapper for fortran coded jet finder using member data for
191 Float_t min_move = 0;
192 Float_t max_move = 0;
194 Float_t prec_bg = 0.;
197 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
198 min_move, max_move, mode, prec_bg, ierror);
199 // Write result to output
200 Int_t njet = Njets();
202 for (Int_t nj=0; nj<njet; nj++)
206 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
210 FindTracksInJetCone();
215 Int_t AliEMCALJetFinder::Njets()
217 // Get number of reconstructed jets
218 return EMCALJETS.njet;
221 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
223 // Get reconstructed jet energy
224 return EMCALJETS.etj[i];
227 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
229 // Get reconstructed jet phi from leading particle
230 return EMCALJETS.phij[0][i];
233 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
235 // Get reconstructed jet phi from weighting
236 return EMCALJETS.phij[1][i];
239 Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
241 // Get reconstructed jet eta from leading particles
242 return EMCALJETS.etaj[0][i];
246 Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
248 // Get reconstructed jet eta from weighting
249 return EMCALJETS.etaj[1][i];
252 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
254 // Set grid cell size
255 EMCALCELLGEO.etaCellSize = eta;
256 EMCALCELLGEO.phiCellSize = phi;
261 void AliEMCALJetFinder::SetConeRadius(Float_t par)
263 // Set jet cone radius
264 EMCALJETPARAM.coneRad = par;
268 void AliEMCALJetFinder::SetEtSeed(Float_t par)
270 // Set et cut for seeds
271 EMCALJETPARAM.etSeed = par;
275 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
277 // Set minimum jet et
278 EMCALJETPARAM.ejMin = par;
282 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
284 // Set et cut per cell
285 EMCALJETPARAM.etMin = par;
289 void AliEMCALJetFinder::SetPtCut(Float_t par)
291 // Set pT cut on charged tracks
296 void AliEMCALJetFinder::Test()
299 // Test the finder call
301 const Int_t nmax = 30000;
303 Int_t ncell_tot = 100;
308 Float_t min_move = 0;
309 Float_t max_move = 0;
315 Find(ncell, ncell_tot, etc, etac, phic,
316 min_move, max_move, mode, prec_bg, ierror);
324 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
329 TClonesArray &lrawcl = *fJets;
330 new(lrawcl[fNjets++]) AliEMCALJet(jet);
333 void AliEMCALJetFinder::ResetJets()
342 void AliEMCALJetFinder::WriteJets()
345 // Add all jets to the list
347 const Int_t kBufferSize = 4000;
348 TTree *pK = gAlice->TreeK();
349 const char* file = (pK->GetCurrentFile())->GetName();
351 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
354 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
356 if (fJets && gAlice->TreeR()) {
357 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
363 Int_t njet = Njets();
364 for (Int_t nj = 0; nj < njet; nj++)
370 Int_t nev = gAlice->GetHeader()->GetEvent();
371 gAlice->TreeR()->Fill();
373 sprintf(hname,"TreeR%d", nev);
374 gAlice->TreeR()->Write(hname);
375 gAlice->TreeR()->Reset();
379 void AliEMCALJetFinder::BookLego()
382 // Book histo for discretisation
385 // Get geometry parameters from
386 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
387 AliEMCALGeometry* geom =
388 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
389 fNbinEta = geom->GetNZ();
390 fNbinPhi = geom->GetNPhi();
391 const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
392 const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
393 const Float_t etaMin = geom->GetArm1EtaMin();
394 const Float_t etaMax = geom->GetArm1EtaMax();
395 fDphi = (phiMax-phiMin)/fNbinEta;
396 fDeta = 1.4/fNbinEta;
397 fNtot = fNbinPhi*fNbinEta;
399 // Don't add histos to the current directory
400 TH2::AddDirectory(0);
403 fLego = new TH2F("legoH","eta-phi",
404 fNbinEta, etaMin, etaMax,
405 fNbinPhi, phiMin, phiMax);
408 fLegoB = new TH2F("legoB","eta-phi",
409 fNbinEta, etaMin, etaMax,
410 fNbinPhi, phiMin, phiMax);
414 void AliEMCALJetFinder::DumpLego()
417 // Dump lego histo into array
420 for (Int_t i = 1; i < fNbinEta; i++) {
421 for (Int_t j = 1; j < fNbinPhi; j++) {
422 Float_t e = fLego->GetBinContent(i,j);
423 if (e <=0.) continue;
425 TAxis* Xaxis = fLego->GetXaxis();
426 TAxis* Yaxis = fLego->GetYaxis();
427 Float_t eta = Xaxis->GetBinCenter(i);
428 Float_t phi = Yaxis->GetBinCenter(j);
430 fEtaCell[fNcell] = eta;
431 fPhiCell[fNcell] = phi;
438 void AliEMCALJetFinder::ResetMap()
441 // Reset eta-phi array
443 for (Int_t i=0; i<30000; i++)
452 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
455 // Fill Cells with track information
460 if (!fLego) BookLego();
462 if (flag == 0) fLego->Reset();
464 // Access particle information
465 Int_t npart = (gAlice->GetHeader())->GetNprimary();
469 // 1: selected for jet finding
472 if (fTrackList) delete[] fTrackList;
473 if (fPtT) delete[] fPtT;
474 if (fEtaT) delete[] fEtaT;
475 if (fPhiT) delete[] fPhiT;
477 fTrackList = new Int_t [npart];
478 fPtT = new Float_t[npart];
479 fEtaT = new Float_t[npart];
480 fPhiT = new Float_t[npart];
485 for (Int_t part = 2; part < npart; part++) {
486 TParticle *MPart = gAlice->Particle(part);
487 Int_t mpart = MPart->GetPdgCode();
488 Int_t child1 = MPart->GetFirstDaughter();
489 Float_t pT = MPart->Pt();
490 Float_t p = MPart->P();
491 Float_t phi = MPart->Phi();
492 Float_t eta = MPart->Eta();
495 if (part == 6 || part == 7)
497 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
498 part-5, pT, eta, phi);
503 // printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
504 // part, mpart, pT, eta, phi);
508 fTrackList[part] = 0;
513 if (part < 2) continue;
514 if (pT == 0 || pT < fPtCut) continue;
515 TParticlePDG* pdgP = 0;
516 // charged or neutral
518 pdgP = MPart->GetPDG();
519 if (pdgP->Charge() == 0) {
523 if (mpart != kNeutron &&
524 mpart != kNeutronBar &&
525 mpart != kK0Long) continue;
530 if (TMath::Abs(mpart) <= 6 ||
532 mpart == 92) continue;
534 if (TMath::Abs(eta) > 0.7) continue;
535 if (phi*180./TMath::Pi() > 120.) continue;
537 if (child1 >= 0 && child1 < npart) continue;
540 printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
541 part, mpart, child1, eta, phi, pT);
544 // phi propagation for hadronic correction
546 Bool_t curls = kFALSE;
547 Float_t dphi = PropagatePhi(pT, pdgP->Charge(), curls);
548 if (fDebug > 1) printf("\n Delta phi %f pT%f", dphi, pT);
549 if (curls && fDebug > 1) printf("\n Track is curling %f", pT);
550 Float_t phiHC = phi + dphi;
552 // Momentum smearing goes here ...
555 p = AliEMCALFast::SmearMomentum(1,p);
558 // Tracking Efficiency and TPC acceptance goes here ...
560 Float_t eff = AliEMCALFast::Efficiency(1,p);
561 if (AliEMCALFast::RandomReject(eff)) continue;
564 // Correction of Hadronic Energy goes here ...
569 if (!fHadronCorrector)
570 Fatal("AliEMCALJetFinder",
571 "Hadronic energy correction required but not defined !");
572 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
575 // More to do ? Just think about it !
577 fTrackList[part] = 1;
581 fLego->Fill(eta, phi, p);
582 if (fHCorrection && !curls) fLego->Fill(eta, phiHC, -dpH);
587 void AliEMCALJetFinder::FillFromHits(Int_t flag)
590 // Fill Cells with hit information
595 if (!fLego) BookLego();
596 // Reset eta-phi map if needed
597 if (flag == 0) fLego->Reset();
598 // Initialize from background event if available
600 // Access hit information
601 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
603 TTree *treeH = gAlice->TreeH();
604 Int_t ntracks = (Int_t) treeH->GetEntries();
611 for (Int_t track=0; track<ntracks;track++) {
613 nbytes += treeH->GetEvent(track);
617 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
619 mHit=(AliEMCALHit*) pEMCAL->NextHit())
621 Float_t x = mHit->X(); // x-pos of hit
622 Float_t y = mHit->Y(); // y-pos
623 Float_t z = mHit->Z(); // z-pos
624 Float_t eloss = mHit->GetEnergy(); // deposited energy
626 Float_t r = TMath::Sqrt(x*x+y*y);
627 Float_t theta = TMath::ATan2(r,z);
628 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
629 Float_t phi = TMath::ATan2(y,x);
631 if (fDebug > 1) printf("\n Hit %f %f %f %f", x, y, z, eloss);
633 fLego->Fill(eta, phi, fSamplingF*eloss);
639 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
642 // Fill Cells with digit information
647 if (!fLego) BookLego();
648 if (flag == 0) fLego->Reset();
655 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
656 TTree *treeD = gAlice->TreeD();
657 TBranchElement* branchDg = (TBranchElement*)
658 treeD->GetBranch("EMCAL");
660 if (!branchDg) Fatal("AliEMCALJetFinder",
661 "Reading digits requested but no digits in file !");
663 branchDg->SetAddress(&digs);
664 Int_t nent = (Int_t) branchDg->GetEntries();
668 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
669 TBranchElement* branchDr = (TBranchElement*)
670 treeD->GetBranch("AliEMCALDigitizer");
671 branchDr->SetAddress(&digr);
674 nbytes += branchDg->GetEntry(0);
675 nbytes += branchDr->GetEntry(0);
677 // Get digitizer parameters
678 Float_t towerADCped = digr->GetTowerpedestal();
679 Float_t towerADCcha = digr->GetTowerchannel();
680 Float_t preshoADCped = digr->GetPreShopedestal();
681 Float_t preshoADCcha = digr->GetPreShochannel();
683 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
684 AliEMCALGeometry* geom =
685 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
688 Int_t ndig = digs->GetEntries();
689 printf("\n Number of Digits: %d %d\n", ndig, nent);
690 printf("\n Parameters: %f %f %f %f\n",
691 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
692 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
699 while ((sdg = (AliEMCALDigit*)(next())))
703 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
705 pedestal = preshoADCped;
706 channel = preshoADCcha;
708 pedestal = towerADCped;
709 channel = towerADCcha;
712 Float_t eta = sdg->GetEta();
713 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
714 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
716 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
717 eta, phi, amp, sdg->GetAmp());
719 fLego->Fill(eta, phi, fSamplingF*amp);
727 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
730 // Fill Cells with hit information
735 if (!fLego) BookLego();
737 if (flag == 0) fLego->Reset();
739 // Flag charged tracks passing through TPC which
740 // are associated to EMCAL Hits
741 BuildTrackFlagTable();
744 // Access particle information
745 TTree *treeH = gAlice->TreeH();
746 Int_t ntracks = (Int_t) treeH->GetEntries();
748 if (fPtT) delete[] fPtT;
749 if (fEtaT) delete[] fEtaT;
750 if (fPhiT) delete[] fPhiT;
752 fPtT = new Float_t[ntracks];
753 fEtaT = new Float_t[ntracks];
754 fPhiT = new Float_t[ntracks];
759 for (Int_t track = 0; track < ntracks; track++) {
760 TParticle *MPart = gAlice->Particle(track);
761 Float_t pT = MPart->Pt();
762 Float_t phi = MPart->Phi();
763 Float_t eta = MPart->Eta();
765 if(fTrackList[track]) {
770 if (track < 2) continue; //Colliding particles?
771 if (pT == 0 || pT < fPtCut) continue;
773 fLego->Fill(eta, phi, pT);
779 void AliEMCALJetFinder::BuildTrackFlagTable() {
781 // Method to generate a lookup table for TreeK particles
782 // which are linked to hits in the EMCAL
784 // --Author: J.L. Klay
786 // Access hit information
787 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
789 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
790 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
792 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
793 fTrackList = new Int_t[nKTrks]; //before generating a new one
795 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
799 TTree *treeH = gAlice->TreeH();
800 Int_t ntracks = (Int_t) treeH->GetEntries();
806 for (Int_t track=0; track<ntracks;track++) {
808 nbytes += treeH->GetEvent(track);
814 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
816 mHit=(AliEMCALHit*) pEMCAL->NextHit())
818 Int_t iTrk = mHit->Track(); // track number
819 Int_t idprim = mHit->GetPrimary(); // primary particle
821 //Determine the origin point of this particle - it made a hit in the EMCAL
822 TParticle *trkPart = gAlice->Particle(iTrk);
823 TParticlePDG *trkPdg = trkPart->GetPDG();
824 Int_t trkCode = trkPart->GetPdgCode();
826 if (trkCode < 10000) { //Big Ions cause problems for
827 trkChg = trkPdg->Charge(); //this function. Since they aren't
828 } else { //likely to make it very far, set
829 trkChg = 0.0; //their charge to 0 for the Flag test
831 Float_t vxTrk = trkPart->Vx();
832 Float_t vyTrk = trkPart->Vy();
833 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
834 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
836 //Loop through the ancestry of the EMCAL entrance particles
837 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
838 while (ancestor != -1) {
839 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
840 TParticlePDG *ancPdg = ancPart->GetPDG();
841 Int_t ancCode = ancPart->GetPdgCode();
843 if (ancCode < 10000) {
844 ancChg = ancPdg->Charge();
848 Float_t vxAnc = ancPart->Vx();
849 Float_t vyAnc = ancPart->Vy();
850 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
851 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
852 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
855 //Determine the origin point of the primary particle associated with the hit
856 TParticle *primPart = gAlice->Particle(idprim);
857 TParticlePDG *primPdg = primPart->GetPDG();
858 Int_t primCode = primPart->GetPdgCode();
860 if (primCode < 10000) {
861 primChg = primPdg->Charge();
865 Float_t vxPrim = primPart->Vx();
866 Float_t vyPrim = primPart->Vy();
867 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
868 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
874 Int_t AliEMCALJetFinder
875 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
881 if (charge == 0) neutral = 1;
883 if (TMath::Abs(code) <= 6 ||
885 code == 92) parton = 1;
887 //It's not a parton, it's charged and it went through the TPC
888 if (!parton && !neutral && radius <= 84.0) flag = 1;
895 void AliEMCALJetFinder::SaveBackgroundEvent()
897 // Saves the eta-phi lego and the tracklist
899 if (fLegoB) fLegoB->Reset();
901 fLego->Copy(*fLegoB);
903 if (fPtB) delete[] fPtB;
904 if (fEtaB) delete[] fEtaB;
905 if (fPhiB) delete[] fPhiB;
906 if (fTrackListB) delete[] fTrackListB;
908 fPtB = new Float_t[fNtS];
909 fEtaB = new Float_t[fNtS];
910 fPhiB = new Float_t[fNtS];
911 fTrackListB = new Int_t [fNtS];
915 for (Int_t i = 0; i < fNt; i++) {
916 if (!fTrackList[i]) continue;
917 fPtB [fNtB] = fPtT [i];
918 fEtaB[fNtB] = fEtaT[i];
919 fPhiB[fNtB] = fPhiT[i];
920 fTrackListB[fNtB] = 1;
926 void AliEMCALJetFinder::InitFromBackground()
930 if (fDebug) printf("\n Initializing from Background");
932 if (fLego) fLego->Reset();
933 fLegoB->Copy(*fLego);
938 void AliEMCALJetFinder::FindTracksInJetCone()
941 // Build list of tracks inside jet cone
944 Int_t njet = Njets();
945 for (Int_t nj = 0; nj < njet; nj++)
947 Float_t etaj = JetEtaW(nj);
948 Float_t phij = JetPhiW(nj);
949 Int_t nT = 0; // number of associated tracks
951 // Loop over particles in current event
952 for (Int_t part = 0; part < fNt; part++) {
953 if (!fTrackList[part]) continue;
954 Float_t phi = fPhiT[part];
955 Float_t eta = fEtaT[part];
956 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
957 (phij-phi)*(phij-phi));
958 if (dr < fConeRadius) {
959 fTrackList[part] = nj+2;
964 // Same for background event if available
967 for (Int_t part = 0; part < fNtB; part++) {
968 Float_t phi = fPhiB[part];
969 Float_t eta = fEtaB[part];
970 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
971 (phij-phi)*(phij-phi));
972 fTrackListB[part] = 1;
974 if (dr < fConeRadius) {
975 printf("\n B sel %d %d %d", part, nj+2, nTB);
976 fTrackListB[part] = nj+2;
980 } // Background available ?
982 Int_t nT0 = nT + nTB;
984 if (nT0 > 50) nT0 = 50;
986 Float_t* ptT = new Float_t[nT0];
987 Float_t* etaT = new Float_t[nT0];
988 Float_t* phiT = new Float_t[nT0];
992 for (Int_t part = 0; part < fNt; part++) {
993 if (fTrackList[part] == nj+2) {
995 for (j=iT-1; j>=0; j--) {
996 if (fPtT[part] > ptT[j]) {
1001 for (j=iT-1; j>=index; j--) {
1002 ptT [j+1] = ptT [j];
1003 etaT[j+1] = etaT[j];
1004 phiT[j+1] = phiT[j];
1006 ptT [index] = fPtT [part];
1007 etaT[index] = fEtaT[part];
1008 phiT[index] = fPhiT[part];
1010 } // particle associated
1011 if (iT > nT0) break;
1015 for (Int_t part = 0; part < fNtB; part++) {
1016 if (fTrackListB[part] == nj+2) {
1018 for (j=iT-1; j>=0; j--) {
1019 if (fPtB[part] > ptT[j]) {
1025 for (j=iT-1; j>=index; j--) {
1026 ptT [j+1] = ptT [j];
1027 etaT[j+1] = etaT[j];
1028 phiT[j+1] = phiT[j];
1030 ptT [index] = fPtB [part];
1031 etaT[index] = fEtaB[part];
1032 phiT[index] = fPhiB[part];
1034 } // particle associated
1035 if (iT > nT0) break;
1037 } // Background available ?
1039 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
1046 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1048 // Propagates phi angle to EMCAL radius
1052 Float_t b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
1054 Float_t rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1058 Float_t rB = 333.56 * pt / b; // [cm]
1060 // check if particle is curling below EMCAL
1061 if (2.*rB < rEMCAL) {
1066 // if not calculate delta phi
1067 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1068 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1069 dPhi = -TMath::Sign(dPhi, charge);
1076 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1078 // dummy for hbook calls