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 **************************************************************************/
16 //---------------------------------------------------------------------
18 // manages the search for jets
19 // Author: jgcn@mda.cinvestav.mx
20 // (code adapted from EMCAL directory)
21 //---------------------------------------------------------------------
23 #include <Riostream.h>
24 #include <TLorentzVector.h>
27 #include "AliUA1JetFinder.h"
28 #include "AliUA1JetHeader.h"
29 #include "UA1Common.h"
30 #include "AliJetReaderHeader.h"
31 #include "AliJetReader.h"
35 ClassImp(AliUA1JetFinder)
37 ////////////////////////////////////////////////////////////////////////
39 AliUA1JetFinder::AliUA1JetFinder()
47 ////////////////////////////////////////////////////////////////////////
49 AliUA1JetFinder::~AliUA1JetFinder()
55 ////////////////////////////////////////////////////////////////////////
58 # define ua1_jet_finder ua1_jet_finder_
63 # define ua1_jet_finder UA1_JET_FINDER
65 # define type_of_call _stdcall
68 extern "C" void type_of_call
69 ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
70 Float_t etc[60000], Float_t etac[60000],
72 Float_t& min_move, Float_t& max_move, Int_t& mode,
73 Float_t& prec_bg, Int_t& ierror);
75 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
77 ////////////////////////////////////////////////////////////////////////
79 void AliUA1JetFinder::FindJets()
82 // initialize event, look for jets, download jet info
84 // initialization in 2 steps
85 // 1) transform input to pt,eta,phi plus lego
88 // 1) transform input to pt,eta,phi plus lego
89 TClonesArray *lvArray = fReader->GetMomentumArray();
90 Int_t nIn = lvArray->GetEntries();
93 // local arrays for input
94 Float_t* enT = new Float_t[nIn];
95 Float_t* ptT = new Float_t[nIn];
96 Float_t* etaT = new Float_t[nIn];
97 Float_t* phiT = new Float_t[nIn];
100 for (Int_t i = 0; i < nIn; i++){
101 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
102 enT[i] = lv->Energy();
105 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
106 if (fReader->GetCutFlag(i) == 1)
107 fLego->Fill(etaT[i], phiT[i], ptT[i]);
109 fJets->SetNinput(nIn);
112 // check enough space! *to be done*
113 Float_t etCell[60000]; //! Cell Energy
114 Float_t etaCell[60000]; //! Cell eta
115 Float_t phiCell[60000]; //! Cell phi
118 TAxis* xaxis = fLego->GetXaxis();
119 TAxis* yaxis = fLego->GetYaxis();
121 for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
122 for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
123 e = fLego->GetBinContent(i,j);
124 if (e > 0.0) e -= fHeader->GetMinCellEt();
126 Float_t eta = xaxis->GetBinCenter(i);
127 Float_t phi = yaxis->GetBinCenter(j);
129 etaCell[nCell] = eta;
130 phiCell[nCell] = phi;
135 // run the algo. Parameters from header
136 Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
137 Float_t minmove = fHeader->GetMinMove();
138 Float_t maxmove = fHeader->GetMaxMove();
139 Int_t mode = fHeader->GetMode();
140 Float_t precbg = fHeader->GetPrecBg();
142 ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
143 minmove, maxmove, mode, precbg, ierror);
146 Int_t * idx = new Int_t[UA1JETS.njet];
147 TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
149 // download jet info.
150 for(Int_t i = 0; i < UA1JETS.njet; i++) {
151 // reject events outside acceptable eta range
152 if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
153 || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
155 Float_t px, py,pz,en; // convert to 4-vector
156 px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
157 py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
158 pz = UA1JETS.etj[idx[i]] /
159 TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
160 en = TMath::Sqrt(px * px + py * py + pz * pz);
161 fJets->AddJet(px, py, pz, en);
164 // find multiplicities and relationship jet-particle
165 // find also percentage of pt from pythia
166 Int_t* injet = new Int_t[nIn];
167 Int_t* sflag = new Int_t[nIn];
168 for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
169 Int_t* mult = new Int_t[fJets->GetNJets()];
170 Int_t* ncell = new Int_t[fJets->GetNJets()];
171 Float_t* percentage = new Float_t[fJets->GetNJets()];
173 for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
174 Float_t pt_sig = 0.0;
176 ncell[i] = UA1JETS.ncellj[i];
177 for (Int_t j = 0; j < nIn; j++) {
178 Float_t deta = etaT[j] - fJets->GetEta(i);
179 Float_t dphi = phiT[j] - fJets->GetPhi(i);
180 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
181 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
182 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
183 if (dr < fHeader->GetRadius() && injet[j] == 0) {
185 if ((fReader->GetCutFlag(j)) == 1 &&
186 (etaT[j] < fHeader->GetLegoEtaMax()) &&
187 (etaT[j] > fHeader->GetLegoEtaMin())) {
190 if (fReader->GetSignalFlag(j) == 1) {
197 percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
198 ((Double_t) fJets->GetPt(i));
200 fJets->SetNCells(ncell);
201 fJets->SetPtFromSignal(percentage);
202 fJets->SetMultiplicities(mult);
203 fJets->SetInJet(injet);
204 fJets->SetEtaIn(etaT);
205 fJets->SetPhiIn(phiT);
207 fJets->SetEtAvg(UA1JETS.etavg);
219 ////////////////////////////////////////////////////////////////////////
221 void AliUA1JetFinder::Reset()
227 ////////////////////////////////////////////////////////////////////////
229 void AliUA1JetFinder::WriteJHeaderToFile()
235 ////////////////////////////////////////////////////////////////////////
237 void AliUA1JetFinder::Init()
239 // initializes some variables
241 dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
242 /((Float_t) fHeader->GetLegoNbinEta());
243 dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
244 /((Float_t) fHeader->GetLegoNbinPhi());
246 UA1CELL.etaCellSize = dEta;
247 UA1CELL.phiCellSize = dPhi;
250 UA1PARA.coneRad = fHeader->GetRadius();
251 UA1PARA.etSeed = fHeader->GetEtSeed();
252 UA1PARA.ejMin = fHeader->GetMinJetEt();
253 UA1PARA.etMin = fHeader->GetMinCellEt();
257 TH2F("legoH","eta-phi",
258 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
259 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
260 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());