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.29 2002/11/21 17:01:40 alibrary
19 Removing AliMCProcess and AliMC
21 Revision 1.28 2002/11/20 14:13:16 morsch
22 - FindChargedJets() added.
23 - Destructor corrected.
24 - Geometry cuts taken from AliEMCALGeometry.
26 Revision 1.27 2002/11/15 17:39:10 morsch
27 GetPythiaParticleName removed.
29 Revision 1.26 2002/10/14 14:55:35 hristov
30 Merging the VirtualMC branch to the main development branch (HEAD)
32 Revision 1.20.4.3 2002/10/10 15:07:49 hristov
33 Updating VirtualMC to v3-09-02
35 Revision 1.25 2002/09/13 10:24:57 morsch
38 Revision 1.24 2002/09/13 10:21:13 morsch
41 Revision 1.23 2002/06/27 09:24:26 morsch
42 Uncomment the TH1::AddDirectory statement.
44 Revision 1.22 2002/05/22 13:48:43 morsch
45 Pdg code added to track list.
47 Revision 1.21 2002/04/27 07:43:08 morsch
48 Calculation of fDphi corrected (Renan Cabrera)
50 Revision 1.20 2002/03/12 01:06:23 pavlinov
51 Testin output from generator
53 Revision 1.19 2002/02/27 00:46:33 pavlinov
54 Added method FillFromParticles()
56 Revision 1.18 2002/02/21 08:48:59 morsch
57 Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
59 Revision 1.17 2002/02/14 08:52:53 morsch
60 Major updates by Aleksei Pavlinov:
61 FillFromPartons, FillFromTracks, jetfinder configuration.
63 Revision 1.16 2002/02/08 11:43:05 morsch
64 SetOutputFileName(..) allows to specify an output file to which the
65 reconstructed jets are written. If not set output goes to input file.
66 Attention call Init() before processing.
68 Revision 1.15 2002/02/02 08:37:18 morsch
69 Formula for rB corrected.
71 Revision 1.14 2002/02/01 08:55:30 morsch
72 Fill map with Et and pT.
74 Revision 1.13 2002/01/31 09:37:36 morsch
75 Geometry parameters in constructor and call of SetCellSize()
77 Revision 1.12 2002/01/23 13:40:23 morsch
78 Fastidious debug print statement removed.
80 Revision 1.11 2002/01/22 17:25:47 morsch
81 Some corrections for event mixing and bg event handling.
83 Revision 1.10 2002/01/22 10:31:44 morsch
84 Some correction for bg mixing.
86 Revision 1.9 2002/01/21 16:28:42 morsch
89 Revision 1.8 2002/01/21 16:05:31 morsch
90 - different phi-bin for hadron correction
91 - provisions for background mixing.
93 Revision 1.7 2002/01/21 15:47:18 morsch
94 Bending radius correctly in cm.
96 Revision 1.6 2002/01/21 12:53:50 morsch
99 Revision 1.5 2002/01/21 12:47:47 morsch
100 Possibility to include K0long and neutrons.
102 Revision 1.4 2002/01/21 11:03:21 morsch
103 Phi propagation introduced in FillFromTracks.
105 Revision 1.3 2002/01/18 05:07:56 morsch
106 - hadronic correction
108 - track selection upon EMCAL information
112 //*-- Authors: Andreas Morsch (CERN)
114 //* Aleksei Pavlinov (WSU)
119 #include <TClonesArray.h>
121 #include <TBranchElement.h>
129 #include <TPaveText.h>
132 #include <TParticle.h>
133 #include <TParticlePDG.h>
134 #include <TPythia6Calls.h>
137 #include "AliEMCALJetFinder.h"
138 #include "AliEMCALFast.h"
139 #include "AliEMCALGeometry.h"
140 #include "AliEMCALHit.h"
141 #include "AliEMCALDigit.h"
142 #include "AliEMCALDigitizer.h"
143 #include "AliEMCALHadronCorrection.h"
144 #include "AliEMCALJetMicroDst.h"
147 #include "AliMagFCM.h"
148 #include "AliEMCAL.h"
149 #include "AliHeader.h"
152 // Interface to FORTRAN
156 ClassImp(AliEMCALJetFinder)
158 //____________________________________________________________________________
159 AliEMCALJetFinder::AliEMCALJetFinder()
161 // Default constructor
180 fHadronCorrector = 0;
188 SetParametersForBgSubtraction();
191 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
195 // Title is used in method GetFileNameForParameters();
197 fJets = new TClonesArray("AliEMCALJet",10000);
199 for (Int_t i = 0; i < 30000; i++)
221 fHadronCorrector = 0;
230 SetMomentumSmearing();
233 SetHadronCorrection();
234 SetSamplingFraction();
237 SetParametersForBgSubtraction();
240 void AliEMCALJetFinder::SetParametersForBgSubtraction
241 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
243 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
244 // at WSU Linux cluster - 11-feb-2002
245 // These parameters must be tuned more carefull !!!
252 //____________________________________________________________________________
253 AliEMCALJetFinder::~AliEMCALJetFinder()
265 delete fhLegoHadrCorr;
268 delete fhCellEMCALEt;
270 delete fhTrackPtBcut;
271 delete fhChPartMultInTpc;
279 delete[] fTrackListB;
287 # define jet_finder_ua1 jet_finder_ua1_
289 # define type_of_call
292 # define jet_finder_ua1 JET_FINDER_UA1
294 # define type_of_call _stdcall
297 extern "C" void type_of_call
298 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
299 Float_t etc[30000], Float_t etac[30000],
301 Float_t& min_move, Float_t& max_move, Int_t& mode,
302 Float_t& prec_bg, Int_t& ierror);
304 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
307 void AliEMCALJetFinder::Init()
310 // Geometry and I/O initialization
314 // Get geometry parameters from EMCAL
318 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
319 AliEMCALGeometry* geom =
320 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
321 fNbinEta = geom->GetNZ();
322 fNbinPhi = geom->GetNPhi();
323 fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
324 fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
325 fEtaMin = geom->GetArm1EtaMin();
326 fEtaMax = geom->GetArm1EtaMax();
327 fDphi = (fPhiMax-fPhiMin)/fNbinPhi;
328 fDeta = (fEtaMax-fEtaMin)/fNbinEta;
329 fNtot = fNbinPhi*fNbinEta;
331 SetCellSize(fDeta, fDphi);
334 if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
337 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
338 Float_t etac[30000], Float_t phic[30000],
339 Float_t min_move, Float_t max_move, Int_t mode,
340 Float_t prec_bg, Int_t ierror)
342 // Wrapper for fortran coded jet finder
343 // Full argument list
344 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
345 min_move, max_move, mode, prec_bg, ierror);
346 // Write result to output
347 if(fWrite) WriteJets();
351 void AliEMCALJetFinder::Find()
353 // Wrapper for fortran coded jet finder using member data for
356 Float_t min_move = fMinMove;
357 Float_t max_move = fMaxMove;
359 Float_t prec_bg = fPrecBg;
362 ResetJets(); // 4-feb-2002 by PAI
364 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
365 min_move, max_move, mode, prec_bg, ierror);
367 // Write result to output
368 Int_t njet = Njets();
370 for (Int_t nj=0; nj<njet; nj++)
373 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
378 FindTracksInJetCone();
379 if(fWrite) WriteJets();
383 Int_t AliEMCALJetFinder::Njets()
385 // Get number of reconstructed jets
386 return EMCALJETS.njet;
389 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
391 // Get reconstructed jet energy
392 return EMCALJETS.etj[i];
395 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
397 // Get reconstructed jet phi from leading particle
398 return EMCALJETS.phij[0][i];
401 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
403 // Get reconstructed jet phi from weighting
404 return EMCALJETS.phij[1][i];
407 Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
409 // Get reconstructed jet eta from leading particles
410 return EMCALJETS.etaj[0][i];
414 Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
416 // Get reconstructed jet eta from weighting
417 return EMCALJETS.etaj[1][i];
420 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
422 // Set grid cell size
423 EMCALCELLGEO.etaCellSize = eta;
424 EMCALCELLGEO.phiCellSize = phi;
427 void AliEMCALJetFinder::SetConeRadius(Float_t par)
429 // Set jet cone radius
430 EMCALJETPARAM.coneRad = par;
434 void AliEMCALJetFinder::SetEtSeed(Float_t par)
436 // Set et cut for seeds
437 EMCALJETPARAM.etSeed = par;
441 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
443 // Set minimum jet et
444 EMCALJETPARAM.ejMin = par;
448 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
450 // Set et cut per cell
451 EMCALJETPARAM.etMin = par;
455 void AliEMCALJetFinder::SetPtCut(Float_t par)
457 // Set pT cut on charged tracks
462 void AliEMCALJetFinder::Test()
465 // Test the finder call
467 const Int_t nmax = 30000;
469 Int_t ncell_tot = 100;
474 Float_t min_move = 0;
475 Float_t max_move = 0;
481 Find(ncell, ncell_tot, etc, etac, phic,
482 min_move, max_move, mode, prec_bg, ierror);
490 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
495 TClonesArray &lrawcl = *fJets;
496 new(lrawcl[fNjets++]) AliEMCALJet(jet);
499 void AliEMCALJetFinder::ResetJets()
508 void AliEMCALJetFinder::WriteJets()
511 // Add all jets to the list
513 const Int_t kBufferSize = 4000;
514 const char* file = 0;
516 Int_t njet = Njets();
518 for (Int_t nj = 0; nj < njet; nj++)
527 // output written to input file
529 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
530 TTree* pK = gAlice->TreeK();
531 file = (pK->GetCurrentFile())->GetName();
533 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
534 if (fJets && gAlice->TreeR()) {
535 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
541 Int_t nev = gAlice->GetHeader()->GetEvent();
542 gAlice->TreeR()->Fill();
544 sprintf(hname,"TreeR%d", nev);
545 gAlice->TreeR()->Write(hname);
546 gAlice->TreeR()->Reset();
549 // Output written to user specified output file
551 TTree* pK = gAlice->TreeK();
552 fInFile = pK->GetCurrentFile();
556 sprintf(hname,"TreeR%d", fEvent);
557 TTree* treeJ = new TTree(hname, "EMCALJets");
558 treeJ->Branch("EMCALJets", &fJets, kBufferSize);
566 void AliEMCALJetFinder::BookLego()
569 // Book histo for discretisation
573 // Don't add histos to the current directory
574 if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
576 TH2::AddDirectory(0);
577 TH1::AddDirectory(0);
581 fLego = new TH2F("legoH","eta-phi",
582 fNbinEta, fEtaMin, fEtaMax,
583 fNbinPhi, fPhiMin, fPhiMax);
586 fLegoB = new TH2F("legoB","eta-phi for BG event",
587 fNbinEta, fEtaMin, fEtaMax,
588 fNbinPhi, fPhiMin, fPhiMax);
591 fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
592 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
594 fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
595 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
596 // Hadron correction map
597 fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
598 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
599 // Hists. for tuning jet finder
600 fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
604 for(Int_t i=1; i<=1000; i++) eTmp[i] = 0.1*i; // step 100 mev
605 for(Int_t i=1001; i<=1100; i++) eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
607 fhCellEt = new TH1F("hCellEt","Cell E_{T} from fLego",
608 eTmp.GetSize()-1, eTmp.GetArray());
609 fhCellEMCALEt = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
610 eTmp.GetSize()-1, eTmp.GetArray());
611 fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
612 eTmp.GetSize()-1, eTmp.GetArray());
613 fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
614 eTmp.GetSize()-1, eTmp.GetArray());
616 fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
617 "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
619 //! first canvas for drawing
620 fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kTRUE);
623 void AliEMCALJetFinder::DumpLego()
626 // Dump lego histo into array
629 TAxis* Xaxis = fLego->GetXaxis();
630 TAxis* Yaxis = fLego->GetYaxis();
631 // fhCellEt->Clear();
633 for (Int_t i = 1; i <= fNbinEta; i++) {
634 for (Int_t j = 1; j <= fNbinPhi; j++) {
635 e = fLego->GetBinContent(i,j);
637 Float_t eta = Xaxis->GetBinCenter(i);
638 Float_t phi = Yaxis->GetBinCenter(j);
640 fEtaCell[fNcell] = eta;
641 fPhiCell[fNcell] = phi;
646 eH = fhLegoEMCAL->GetBinContent(i,j);
647 if(eH > 0.0) fhCellEMCALEt->Fill(eH);
654 void AliEMCALJetFinder::ResetMap()
657 // Reset eta-phi array
659 for (Int_t i=0; i<30000; i++)
668 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
671 // Fill Cells with track information
674 printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
679 if (!fLego) BookLego();
681 if (flag == 0) fLego->Reset();
683 // Access particle information
684 Int_t npart = (gAlice->GetHeader())->GetNprimary();
685 Int_t ntr = (gAlice->GetHeader())->GetNtrack();
686 printf(" : #primary particles %i # tracks %i \n", npart, ntr);
691 // 1: selected for jet finding
694 if (fTrackList) delete[] fTrackList;
695 if (fPtT) delete[] fPtT;
696 if (fEtaT) delete[] fEtaT;
697 if (fPhiT) delete[] fPhiT;
698 if (fPdgT) delete[] fPdgT;
700 fTrackList = new Int_t [npart];
701 fPtT = new Float_t[npart];
702 fEtaT = new Float_t[npart];
703 fPhiT = new Float_t[npart];
704 fPdgT = new Int_t[npart];
708 Float_t chTmp=0.0; // charge of current particle
711 // this is for Pythia ??
712 for (Int_t part = 0; part < npart; part++) {
713 TParticle *MPart = gAlice->Particle(part);
714 Int_t mpart = MPart->GetPdgCode();
715 Int_t child1 = MPart->GetFirstDaughter();
716 Float_t pT = MPart->Pt();
717 Float_t p = MPart->P();
718 Float_t phi = MPart->Phi();
720 if(pT > 0.001) eta = MPart->Eta();
721 Float_t theta = MPart->Theta();
723 printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
724 part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
728 if (part == 6 || part == 7)
730 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
731 part-5, pT, eta, phi);
736 // printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
737 // part, mpart, pT, eta, phi);
740 fTrackList[part] = 0;
741 fPtT[part] = pT; // must be change after correction for resolution !!!
747 if (part < 2) continue;
749 // move to fLego->Fill because hadron correction must apply
750 // if particle hit to EMCAL - 11-feb-2002
751 // if (pT == 0 || pT < fPtCut) continue;
752 TParticlePDG* pdgP = 0;
753 // charged or neutral
754 pdgP = MPart->GetPDG();
755 chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
762 if (mpart != kNeutron &&
763 mpart != kNeutronBar &&
764 mpart != kK0Long) continue;
770 if (TMath::Abs(mpart) <= 6 ||
772 mpart == 92) continue;
774 if (TMath::Abs(eta)<=0.9) fNChTpc++;
776 if (child1 >= 0 && child1 < npart) continue;
778 if (eta > fEtaMax || eta < fEtaMin) continue;
779 if (phi > fPhiMax || phi < fPhiMin ) continue;
782 printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
783 part, mpart, child1, eta, phi, pT, chTmp);
786 // Momentum smearing goes here ...
790 if (fSmear && TMath::Abs(chTmp)) {
791 pw = AliEMCALFast::SmearMomentum(1,p);
792 // p changed - take into account when calculate pT,
795 if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
799 // Tracking Efficiency and TPC acceptance goes here ...
801 if (fEffic && TMath::Abs(chTmp)) {
802 // eff = AliEMCALFast::Efficiency(1,p);
803 eff = 0.95; // for testing 25-feb-2002
804 if(fhEff) fhEff->Fill(p, eff);
805 if (AliEMCALFast::RandomReject(eff)) {
806 if(fDebug >= 5) printf(" reject due to unefficiency ");
811 // Correction of Hadronic Energy goes here ...
814 // phi propagation for hadronic correction
816 Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
817 Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
818 if(TMath::Abs(chTmp)) {
819 // hadr. correction only for charge particle
820 dphi = PropagatePhi(pT, chTmp, curls);
823 printf("\n Delta phi %f pT %f ", dphi, pT);
824 if (curls) printf("\n !! Track is curling");
826 if(!curls) fhTrackPtBcut->Fill(pT);
828 if (fHCorrection && !curls) {
829 if (!fHadronCorrector)
830 Fatal("AliEMCALJetFinder",
831 "Hadronic energy correction required but not defined !");
833 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
834 eTdpH = dpH*TMath::Sin(theta);
836 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
837 phi, phiHC, -eTdpH); // correction is negative
838 fLego->Fill(eta, phiHC, -eTdpH);
839 fhLegoHadrCorr->Fill(eta, phiHC, eTdpH);
843 // More to do ? Just think about it !
845 if (phi > fPhiMax || phi < fPhiMin ) continue;
847 if(TMath::Abs(chTmp) ) { // charge particle
848 if (pT > fPtCut && !curls) {
849 if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
850 eta , phi, pT, fNtS);
851 fLego->Fill(eta, phi, pT);
852 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
853 fTrackList[part] = 1;
856 } else if(ich==0 && fK0N) {
857 // case of n, nbar and K0L
858 if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
859 eta , phi, pT, fNtS);
860 fLego->Fill(eta, phi, pT);
861 fTrackList[part] = 1;
869 void AliEMCALJetFinder::FillFromHits(Int_t flag)
872 // Fill Cells with hit information
876 printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
880 if (!fLego) BookLego();
881 // Reset eta-phi maps if needed
882 if (flag == 0) { // default behavior
884 fhLegoTracks->Reset();
885 fhLegoEMCAL->Reset();
886 fhLegoHadrCorr->Reset();
888 // Initialize from background event if available
890 // Access hit information
891 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
892 TTree *treeH = gAlice->TreeH();
893 Int_t ntracks = (Int_t) treeH->GetEntries();
900 for (Int_t track=0; track<ntracks;track++) {
902 nbytes += treeH->GetEvent(track);
906 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
908 mHit=(AliEMCALHit*) pEMCAL->NextHit())
910 Float_t x = mHit->X(); // x-pos of hit
911 Float_t y = mHit->Y(); // y-pos
912 Float_t z = mHit->Z(); // z-pos
913 Float_t eloss = mHit->GetEnergy(); // deposited energy
915 Float_t r = TMath::Sqrt(x*x+y*y);
916 Float_t theta = TMath::ATan2(r,z);
917 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
918 Float_t phi = TMath::ATan2(y,x);
920 if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
921 // printf("\n Hit %f %f %f %f", x, y, z, eloss);
923 etH = fSamplingF*eloss*TMath::Sin(theta);
924 fLego->Fill(eta, phi, etH);
925 // fhLegoEMCAL->Fill(eta, phi, etH);
928 // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
929 for(Int_t i=0; i<fLego->GetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i];
933 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
936 // Fill Cells with digit information
941 if (!fLego) BookLego();
942 if (flag == 0) fLego->Reset();
949 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
950 TTree *treeD = gAlice->TreeD();
951 TBranchElement* branchDg = (TBranchElement*)
952 treeD->GetBranch("EMCAL");
954 if (!branchDg) Fatal("AliEMCALJetFinder",
955 "Reading digits requested but no digits in file !");
957 branchDg->SetAddress(&digs);
958 Int_t nent = (Int_t) branchDg->GetEntries();
962 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
963 TBranchElement* branchDr = (TBranchElement*)
964 treeD->GetBranch("AliEMCALDigitizer");
965 branchDr->SetAddress(&digr);
968 nbytes += branchDg->GetEntry(0);
969 nbytes += branchDr->GetEntry(0);
971 // Get digitizer parameters
972 Float_t towerADCped = digr->GetTowerpedestal();
973 Float_t towerADCcha = digr->GetTowerchannel();
974 Float_t preshoADCped = digr->GetPreShopedestal();
975 Float_t preshoADCcha = digr->GetPreShochannel();
977 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
978 AliEMCALGeometry* geom =
979 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
982 Int_t ndig = digs->GetEntries();
983 printf("\n Number of Digits: %d %d\n", ndig, nent);
984 printf("\n Parameters: %f %f %f %f\n",
985 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
986 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
993 while ((sdg = (AliEMCALDigit*)(next())))
997 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
999 pedestal = preshoADCped;
1000 channel = preshoADCcha;
1002 pedestal = towerADCped;
1003 channel = towerADCcha;
1006 Float_t eta = sdg->GetEta();
1007 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1008 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1010 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1011 eta, phi, amp, sdg->GetAmp());
1013 fLego->Fill(eta, phi, fSamplingF*amp);
1021 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1024 // Fill Cells with hit information
1029 if (!fLego) BookLego();
1031 if (flag == 0) fLego->Reset();
1033 // Flag charged tracks passing through TPC which
1034 // are associated to EMCAL Hits
1035 BuildTrackFlagTable();
1038 // Access particle information
1039 TTree *treeK = gAlice->TreeK();
1040 Int_t ntracks = (Int_t) treeK->GetEntries();
1042 if (fPtT) delete[] fPtT;
1043 if (fEtaT) delete[] fEtaT;
1044 if (fPhiT) delete[] fPhiT;
1045 if (fPdgT) delete[] fPdgT;
1047 fPtT = new Float_t[ntracks];
1048 fEtaT = new Float_t[ntracks];
1049 fPhiT = new Float_t[ntracks];
1050 fPdgT = new Int_t[ntracks];
1055 for (Int_t track = 0; track < ntracks; track++) {
1056 TParticle *MPart = gAlice->Particle(track);
1057 Float_t pT = MPart->Pt();
1058 Float_t phi = MPart->Phi();
1059 Float_t eta = MPart->Eta();
1061 if(fTrackList[track]) {
1065 fPdgT[track] = MPart->GetPdgCode();
1067 if (track < 2) continue; //Colliding particles?
1068 if (pT == 0 || pT < fPtCut) continue;
1070 fLego->Fill(eta, phi, pT);
1076 void AliEMCALJetFinder::FillFromParticles()
1078 // 26-feb-2002 PAI - for checking all chain
1079 // Work on particles level; accept all particle (not neutrino )
1081 Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law
1085 if (!fLego) BookLego();
1088 // Access particles information
1089 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1090 if (fDebug >= 2 || npart<=0) {
1091 printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1092 if(npart<=0) return;
1096 RearrangeParticlesMemory(npart);
1098 // Go through the particles
1099 Int_t mpart, child1, child2, geantPdg;
1100 Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1102 for (Int_t part = 0; part < npart; part++) {
1104 fTrackList[part] = 0;
1106 MPart = gAlice->Particle(part);
1107 mpart = MPart->GetPdgCode();
1108 child1 = MPart->GetFirstDaughter();
1109 child2 = MPart->GetLastDaughter();
1117 e = MPart->Energy();
1119 // see pyedit in Pythia's text
1121 if (IsThisPartonsOrDiQuark(mpart)) continue;
1122 printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1123 part, mpart, geantPdg, px, py, pz, e, child1, child2);
1125 // exclude partons (21 - gluon, 92 - string)
1128 // exclude neutrinous also ??
1129 if (fDebug >= 11 && pT>0.01)
1130 printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1131 part, mpart, eta, phi, pT);
1136 fPdgT[part] = mpart;
1140 if (child1 >= 0 && child1 < npart) continue;
1142 // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1143 // part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1150 if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1152 if (eta > fEtaMax || eta < fEtaMin) continue;
1153 if (phi > fPhiMax || phi < fPhiMin ) continue;
1155 if(fK0N==0 ) { // exclude neutral hadrons
1156 if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1158 fTrackList[part] = 1;
1159 fLego->Fill(eta, phi, pT);
1162 printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n",
1165 if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1168 void AliEMCALJetFinder::FillFromPartons()
1170 // function under construction - 13-feb-2002 PAI
1173 printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1177 if (!fLego) BookLego();
1180 // Access particle information
1181 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1182 if (fDebug >= 2 || npart<=0)
1183 printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1184 fNt = 0; // for FindTracksInJetCone
1187 // Go through the partons
1189 for (Int_t part = 8; part < npart; part++) {
1190 TParticle *MPart = gAlice->Particle(part);
1191 Int_t mpart = MPart->GetPdgCode();
1192 // Int_t child1 = MPart->GetFirstDaughter();
1193 Float_t pT = MPart->Pt();
1194 // Float_t p = MPart->P();
1195 Float_t phi = MPart->Phi();
1196 Float_t eta = MPart->Eta();
1197 // Float_t theta = MPart->Theta();
1198 statusCode = MPart->GetStatusCode();
1200 // accept partons (21 - gluon, 92 - string)
1201 if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1202 if (fDebug > 1 && pT>0.01)
1203 printf("\n part:%5d mpart %5d status %5d eta %8.2f phi %8.2f pT %8.2f ",
1204 part, mpart, statusCode, eta, phi, pT);
1205 // if (fDebug >= 3) MPart->Print();
1206 // accept partons before fragmentation - p.57 in Pythia manual
1207 // if(statusCode != 1) continue;
1209 if (eta > fEtaMax || eta < fEtaMin) continue;
1210 if (phi > fPhiMax || phi < fPhiMin ) continue;
1212 // if (child1 >= 0 && child1 < npart) continue;
1215 fLego->Fill(eta, phi, pT);
1221 void AliEMCALJetFinder::BuildTrackFlagTable() {
1223 // Method to generate a lookup table for TreeK particles
1224 // which are linked to hits in the EMCAL
1226 // --Author: J.L. Klay
1228 // Access hit information
1229 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1231 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
1232 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
1234 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
1235 fTrackList = new Int_t[nKTrks]; //before generating a new one
1237 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
1241 TTree *treeH = gAlice->TreeH();
1242 Int_t ntracks = (Int_t) treeH->GetEntries();
1248 for (Int_t track=0; track<ntracks;track++) {
1249 gAlice->ResetHits();
1250 nbytes += treeH->GetEvent(track);
1256 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1258 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1260 Int_t iTrk = mHit->Track(); // track number
1261 Int_t idprim = mHit->GetPrimary(); // primary particle
1263 //Determine the origin point of this particle - it made a hit in the EMCAL
1264 TParticle *trkPart = gAlice->Particle(iTrk);
1265 TParticlePDG *trkPdg = trkPart->GetPDG();
1266 Int_t trkCode = trkPart->GetPdgCode();
1268 if (trkCode < 10000) { //Big Ions cause problems for
1269 trkChg = trkPdg->Charge(); //this function. Since they aren't
1270 } else { //likely to make it very far, set
1271 trkChg = 0.0; //their charge to 0 for the Flag test
1273 Float_t vxTrk = trkPart->Vx();
1274 Float_t vyTrk = trkPart->Vy();
1275 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1276 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1278 //Loop through the ancestry of the EMCAL entrance particles
1279 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
1280 while (ancestor != -1) {
1281 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
1282 TParticlePDG *ancPdg = ancPart->GetPDG();
1283 Int_t ancCode = ancPart->GetPdgCode();
1285 if (ancCode < 10000) {
1286 ancChg = ancPdg->Charge();
1290 Float_t vxAnc = ancPart->Vx();
1291 Float_t vyAnc = ancPart->Vy();
1292 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1293 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1294 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
1297 //Determine the origin point of the primary particle associated with the hit
1298 TParticle *primPart = gAlice->Particle(idprim);
1299 TParticlePDG *primPdg = primPart->GetPDG();
1300 Int_t primCode = primPart->GetPdgCode();
1302 if (primCode < 10000) {
1303 primChg = primPdg->Charge();
1307 Float_t vxPrim = primPart->Vx();
1308 Float_t vyPrim = primPart->Vy();
1309 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1310 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1316 Int_t AliEMCALJetFinder
1317 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1323 if (charge == 0) neutral = 1;
1325 if (TMath::Abs(code) <= 6 ||
1327 code == 92) parton = 1;
1329 //It's not a parton, it's charged and it went through the TPC
1330 if (!parton && !neutral && radius <= 84.0) flag = 1;
1337 void AliEMCALJetFinder::SaveBackgroundEvent()
1339 // Saves the eta-phi lego and the tracklist
1343 (*fLegoB) = (*fLegoB) + (*fLego);
1345 printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1346 fLegoB->Integral(), fLego->Integral());
1349 if (fPtB) delete[] fPtB;
1350 if (fEtaB) delete[] fEtaB;
1351 if (fPhiB) delete[] fPhiB;
1352 if (fPdgB) delete[] fPdgB;
1353 if (fTrackListB) delete[] fTrackListB;
1355 fPtB = new Float_t[fNtS];
1356 fEtaB = new Float_t[fNtS];
1357 fPhiB = new Float_t[fNtS];
1358 fPdgB = new Int_t [fNtS];
1359 fTrackListB = new Int_t [fNtS];
1363 for (Int_t i = 0; i < fNt; i++) {
1364 if (!fTrackList[i]) continue;
1365 fPtB [fNtB] = fPtT [i];
1366 fEtaB[fNtB] = fEtaT[i];
1367 fPhiB[fNtB] = fPhiT[i];
1368 fPdgB[fNtB] = fPdgT[i];
1369 fTrackListB[fNtB] = 1;
1373 printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1376 void AliEMCALJetFinder::InitFromBackground()
1380 if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1384 (*fLego) = (*fLego) + (*fLegoB);
1386 printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1387 fLego->Integral(), fLegoB->Integral());
1389 printf(" => fLego undefined \n");
1394 void AliEMCALJetFinder::FindTracksInJetCone()
1397 // Build list of tracks inside jet cone
1400 Int_t njet = Njets();
1401 for (Int_t nj = 0; nj < njet; nj++)
1403 Float_t etaj = JetEtaW(nj);
1404 Float_t phij = JetPhiW(nj);
1405 Int_t nT = 0; // number of associated tracks
1407 // Loop over particles in current event
1408 for (Int_t part = 0; part < fNt; part++) {
1409 if (!fTrackList[part]) continue;
1410 Float_t phi = fPhiT[part];
1411 Float_t eta = fEtaT[part];
1412 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1413 (phij-phi)*(phij-phi));
1414 if (dr < fConeRadius) {
1415 fTrackList[part] = nj+2;
1420 // Same for background event if available
1423 for (Int_t part = 0; part < fNtB; part++) {
1424 Float_t phi = fPhiB[part];
1425 Float_t eta = fEtaB[part];
1426 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1427 (phij-phi)*(phij-phi));
1428 fTrackListB[part] = 1;
1430 if (dr < fConeRadius) {
1431 fTrackListB[part] = nj+2;
1435 } // Background available ?
1437 Int_t nT0 = nT + nTB;
1438 printf("Total number of tracks %d\n", nT0);
1440 if (nT0 > 1000) nT0 = 1000;
1442 Float_t* ptT = new Float_t[nT0];
1443 Float_t* etaT = new Float_t[nT0];
1444 Float_t* phiT = new Float_t[nT0];
1445 Int_t* pdgT = new Int_t[nT0];
1450 for (Int_t part = 0; part < fNt; part++) {
1451 if (fTrackList[part] == nj+2) {
1453 for (j=iT-1; j>=0; j--) {
1454 if (fPtT[part] > ptT[j]) {
1459 for (j=iT-1; j>=index; j--) {
1460 ptT [j+1] = ptT [j];
1461 etaT[j+1] = etaT[j];
1462 phiT[j+1] = phiT[j];
1463 pdgT[j+1] = pdgT[j];
1465 ptT [index] = fPtT [part];
1466 etaT[index] = fEtaT[part];
1467 phiT[index] = fPhiT[part];
1468 pdgT[index] = fPdgT[part];
1470 } // particle associated
1471 if (iT > nT0) break;
1475 for (Int_t part = 0; part < fNtB; part++) {
1476 if (fTrackListB[part] == nj+2) {
1478 for (j=iT-1; j>=0; j--) {
1479 if (fPtB[part] > ptT[j]) {
1485 for (j=iT-1; j>=index; j--) {
1486 ptT [j+1] = ptT [j];
1487 etaT[j+1] = etaT[j];
1488 phiT[j+1] = phiT[j];
1489 pdgT[j+1] = pdgT[j];
1491 ptT [index] = fPtB [part];
1492 etaT[index] = fEtaB[part];
1493 phiT[index] = fPhiB[part];
1494 pdgT[index] = fPdgB[part];
1496 } // particle associated
1497 if (iT > nT0) break;
1499 } // Background available ?
1501 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1510 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1512 // Propagates phi angle to EMCAL radius
1514 static Float_t b = 0.0, rEMCAL = -1.0;
1517 b = gAlice->Field()->SolenoidField();
1518 // Get EMCAL radius in cm
1519 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1520 printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
1529 Float_t rB = 3335.6 * pt / b; // [cm] (case of |charge|=1)
1531 // check if particle is curling below EMCAL
1532 if (2.*rB < rEMCAL) {
1537 // if not calculate delta phi
1538 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1539 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1540 dPhi = -TMath::Sign(dPhi, charge);
1545 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1547 // dummy for hbook calls
1551 void AliEMCALJetFinder::DrawLego(Char_t *opt)
1554 void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
1555 {fhLegoEMCAL->Draw(opt);}
1557 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1559 static TPaveText *varLabel=0;
1561 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1571 fhCellEMCALEt->Draw();
1576 fhTrackPtBcut->SetLineColor(2);
1577 fhTrackPtBcut->Draw("same");
1582 varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1583 varLabel->SetTextAlign(12);
1584 varLabel->SetFillColor(19); // see TAttFill
1585 TString tmp(GetTitle());
1586 varLabel->ReadFile(GetFileNameForParameters());
1590 if(mode) { // for saving picture to the file
1591 TString stmp(GetFileNameForParameters());
1592 stmp.ReplaceAll("_Par.txt",".ps");
1593 fC1->Print(stmp.Data());
1597 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1600 if(mode==0) file = stdout; // output to terminal
1602 file = fopen(GetFileNameForParameters(),"w");
1603 if(file==0) file = stdout;
1605 fprintf(file,"==== Filling lego ==== \n");
1606 fprintf(file,"Smearing %6i ", fSmear);
1607 fprintf(file,"Efficiency %6i\n", fEffic);
1608 fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
1609 fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1610 fprintf(file,"==== Jet finding ==== \n");
1611 fprintf(file,"Cone radius %6.2f ", fConeRadius);
1612 fprintf(file,"Seed E_{T} %6.1f\n", fEtSeed);
1613 fprintf(file,"Min E_{T} of cell %6.1f ", fMinCellEt);
1614 fprintf(file,"Min E_{T} of jet %6.1f\n", fMinJetEt);
1616 fprintf(file,"==== Bg subtraction ==== \n");
1617 fprintf(file,"BG subtraction %6i ", fMode);
1618 fprintf(file,"Min cone move %6.2f\n", fMinMove);
1619 fprintf(file,"Max cone move %6.2f ", fMaxMove);
1620 fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1622 fprintf(file,"==== No Bg subtraction ==== \n");
1623 if(file != stdout) fclose(file);
1626 void AliEMCALJetFinder::DrawLegos()
1629 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1633 gStyle->SetOptStat(111111);
1635 Int_t nent1, nent2, nent3, nent4;
1636 Double_t int1, int2, int3, int4;
1637 nent1 = (Int_t)fLego->GetEntries();
1638 int1 = fLego->Integral();
1640 if(int1) fLego->Draw("lego");
1642 nent2 = (Int_t)fhLegoTracks->GetEntries();
1643 int2 = fhLegoTracks->Integral();
1645 if(int2) fhLegoTracks->Draw("lego");
1647 nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1648 int3 = fhLegoEMCAL->Integral();
1650 if(int3) fhLegoEMCAL->Draw("lego");
1652 nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1653 int4 = fhLegoHadrCorr->Integral();
1655 if(int4) fhLegoHadrCorr->Draw("lego");
1657 // just for checking
1658 printf(" Integrals \n");
1659 printf("lego %10.3f \ntrack %10.3f \nhits %10.3f \nHCorr %10.3f\n-- %10.3f(must be 0)\n",
1660 int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1663 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1666 if(strlen(dir)) tmp = dir;
1672 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1673 { // See FillFromTracks() - npart must be positive
1674 if (fTrackList) delete[] fTrackList;
1675 if (fPtT) delete[] fPtT;
1676 if (fEtaT) delete[] fEtaT;
1677 if (fPhiT) delete[] fPhiT;
1678 if (fPdgT) delete[] fPdgT;
1681 fTrackList = new Int_t [npart];
1682 fPtT = new Float_t[npart];
1683 fEtaT = new Float_t[npart];
1684 fPhiT = new Float_t[npart];
1685 fPdgT = new Int_t[npart];
1687 printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1691 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1693 Int_t absPdg = TMath::Abs(pdg);
1694 if(absPdg<=6) return kTRUE; // quarks
1695 if(pdg == 21) return kTRUE; // gluon
1696 if(pdg == 92) return kTRUE; // string
1698 // see p.51 of Pythia Manual
1699 // Not include diquarks with c and b quark - 4-mar-2002
1700 // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1
1701 static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303};
1702 for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks
1707 void AliEMCALJetFinder::FindChargedJet()
1710 // Find jet from charged particle information only
1725 for (part = 0; part < fNt; part++) {
1726 if (!fTrackList[part]) continue;
1727 if (fPtT[part] > fEtSeed) nseed++;
1729 printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1730 Int_t* iSeeds = new Int_t[nseed];
1733 for (part = 0; part < fNt; part++) {
1734 if (!fTrackList[part]) continue;
1735 if (fPtT[part] > fEtSeed) iSeeds[nseed++] = part;
1746 // Find seed with highest pt
1748 Float_t ptmax = -1.;
1751 for (seed = 0; seed < nseed; seed++) {
1752 if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1757 if (ptmax < 0.) break;
1758 jndex = iSeeds[index];
1760 // Remove from the list
1762 printf("\n Next Seed %d %f", jndex, ptmax);
1764 // Find tracks in cone around seed
1766 Float_t phiSeed = fPhiT[jndex];
1767 Float_t etaSeed = fEtaT[jndex];
1773 for (part = 0; part < fNt; part++) {
1774 if (!fTrackList[part]) continue;
1775 Float_t deta = fEtaT[part] - etaSeed;
1776 Float_t dphi = fPhiT[part] - phiSeed;
1777 Float_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1778 if (dR < fConeRadius) {
1780 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1781 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1782 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1783 Float_t pz = fPtT[part] / TMath::Tan(theta);
1788 // if seed, remove it
1790 for (seed = 0; seed < nseed; seed++) {
1791 if (part == iSeeds[seed]) iSeeds[seed] = -1;
1797 // Estimate of jet direction
1798 Float_t phiJ = TMath::ATan2(pyJ, pxJ);
1799 Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1800 Float_t etaJ = TMath::Log(TMath::Tan(thetaJ / 2.));
1801 // Float_t ptJ = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1804 // Sum up all energy
1806 Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1807 Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1808 Int_t dIphi = Int_t(fConeRadius / fDphi);
1809 Int_t dIeta = Int_t(fConeRadius / fDeta);
1812 for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1813 for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1814 if (iPhi < 0 || iEta < 0) continue;
1815 Float_t dPhi = fPhiMin + iPhi * fDphi;
1816 Float_t dEta = fEtaMin + iEta * fDeta;
1817 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1818 sumE += fLego->GetBinContent(iEta, iPhi);
1824 fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1825 FindTracksInJetCone();
1826 printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1827 printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1828 printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1830 EMCALJETS.njet = njets;
1831 if (fWrite) WriteJets();