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():
43 // Default constructor
46 ////////////////////////////////////////////////////////////////////////
48 AliUA1JetFinder::~AliUA1JetFinder()
54 ////////////////////////////////////////////////////////////////////////
57 # define ua1_jet_finder ua1_jet_finder_
62 # define ua1_jet_finder UA1_JET_FINDER
64 # define type_of_call _stdcall
67 extern "C" void type_of_call
68 ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
69 Float_t etc[60000], Float_t etac[60000],
71 Float_t& min_move, Float_t& max_move, Int_t& mode,
72 Float_t& prec_bg, Int_t& ierror);
74 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
76 ////////////////////////////////////////////////////////////////////////
78 void AliUA1JetFinder::FindJets()
81 // initialize event, look for jets, download jet info
83 // initialization in 2 steps
84 // 1) transform input to pt,eta,phi plus lego
87 // 1) transform input to pt,eta,phi plus lego
88 TClonesArray *lvArray = fReader->GetMomentumArray();
89 Int_t nIn = lvArray->GetEntries();
92 // local arrays for input
93 Float_t* enT = new Float_t[nIn];
94 Float_t* ptT = new Float_t[nIn];
95 Float_t* etaT = new Float_t[nIn];
96 Float_t* phiT = new Float_t[nIn];
99 for (Int_t i = 0; i < nIn; i++){
100 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
101 enT[i] = lv->Energy();
104 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
105 if (fReader->GetCutFlag(i) == 1)
106 fLego->Fill(etaT[i], phiT[i], ptT[i]);
108 fJets->SetNinput(nIn);
111 // check enough space! *to be done*
112 Float_t etCell[60000]; //! Cell Energy
113 Float_t etaCell[60000]; //! Cell eta
114 Float_t phiCell[60000]; //! Cell phi
117 TAxis* xaxis = fLego->GetXaxis();
118 TAxis* yaxis = fLego->GetYaxis();
120 for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
121 for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
122 e = fLego->GetBinContent(i,j);
123 if (e > 0.0) e -= fHeader->GetMinCellEt();
125 Float_t eta = xaxis->GetBinCenter(i);
126 Float_t phi = yaxis->GetBinCenter(j);
128 etaCell[nCell] = eta;
129 phiCell[nCell] = phi;
134 // run the algo. Parameters from header
135 Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
136 Float_t minmove = fHeader->GetMinMove();
137 Float_t maxmove = fHeader->GetMaxMove();
138 Int_t mode = fHeader->GetMode();
139 Float_t precbg = fHeader->GetPrecBg();
141 ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
142 minmove, maxmove, mode, precbg, ierror);
145 Int_t * idx = new Int_t[UA1JETS.njet];
146 TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
148 // download jet info.
149 for(Int_t i = 0; i < UA1JETS.njet; i++) {
150 // reject events outside acceptable eta range
151 if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
152 || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
154 Float_t px, py,pz,en; // convert to 4-vector
155 px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
156 py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
157 pz = UA1JETS.etj[idx[i]] /
158 TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
159 en = TMath::Sqrt(px * px + py * py + pz * pz);
160 fJets->AddJet(px, py, pz, en);
163 // find multiplicities and relationship jet-particle
164 // find also percentage of pt from pythia
165 Int_t* injet = new Int_t[nIn];
166 Int_t* sflag = new Int_t[nIn];
167 for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
168 Int_t* mult = new Int_t[fJets->GetNJets()];
169 Int_t* ncell = new Int_t[fJets->GetNJets()];
170 Float_t* percentage = new Float_t[fJets->GetNJets()];
172 for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
173 Float_t pt_sig = 0.0;
175 ncell[i] = UA1JETS.ncellj[i];
176 for (Int_t j = 0; j < nIn; j++) {
177 Float_t deta = etaT[j] - fJets->GetEta(i);
178 Float_t dphi = phiT[j] - fJets->GetPhi(i);
179 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
180 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
181 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
182 if (dr < fHeader->GetRadius() && injet[j] == 0) {
184 if ((fReader->GetCutFlag(j)) == 1 &&
185 (etaT[j] < fHeader->GetLegoEtaMax()) &&
186 (etaT[j] > fHeader->GetLegoEtaMin())) {
189 if (fReader->GetSignalFlag(j) == 1) {
196 percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
197 ((Double_t) fJets->GetPt(i));
199 fJets->SetNCells(ncell);
200 fJets->SetPtFromSignal(percentage);
201 fJets->SetMultiplicities(mult);
202 fJets->SetInJet(injet);
203 fJets->SetEtaIn(etaT);
204 fJets->SetPhiIn(phiT);
206 fJets->SetEtAvg(UA1JETS.etavg);
218 ////////////////////////////////////////////////////////////////////////
220 void AliUA1JetFinder::Reset()
226 ////////////////////////////////////////////////////////////////////////
228 void AliUA1JetFinder::WriteJHeaderToFile()
234 ////////////////////////////////////////////////////////////////////////
236 void AliUA1JetFinder::Init()
238 // initializes some variables
240 dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
241 /((Float_t) fHeader->GetLegoNbinEta());
242 dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
243 /((Float_t) fHeader->GetLegoNbinPhi());
245 UA1CELL.etaCellSize = dEta;
246 UA1CELL.phiCellSize = dPhi;
249 UA1PARA.coneRad = fHeader->GetRadius();
250 UA1PARA.etSeed = fHeader->GetEtSeed();
251 UA1PARA.ejMin = fHeader->GetMinJetEt();
252 UA1PARA.etMin = fHeader->GetMinCellEt();
256 TH2F("legoH","eta-phi",
257 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
258 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
259 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());