- Changed module names to be the same as the libraries names;
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.cxx
CommitLineData
99e5fe42 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
99e5fe42 15
16//---------------------------------------------------------------------
17// UA1 Jet finder
18// manages the search for jets
19// Author: jgcn@mda.cinvestav.mx
20// (code adapted from EMCAL directory)
21//---------------------------------------------------------------------
22
23#include <Riostream.h>
24#include <TLorentzVector.h>
83a444b1 25#include <TFile.h>
99e5fe42 26#include <TH2F.h>
27#include "AliUA1JetFinder.h"
28#include "AliUA1JetHeader.h"
29#include "UA1Common.h"
8011d399 30#include "AliJetReaderHeader.h"
99e5fe42 31#include "AliJetReader.h"
32#include "AliJet.h"
33
34
35ClassImp(AliUA1JetFinder)
36
37////////////////////////////////////////////////////////////////////////
38
1b7d5d7e 39AliUA1JetFinder::AliUA1JetFinder():
40 fHeader(0x0),
41 fLego(0x0)
99e5fe42 42{
1b7d5d7e 43 // Default constructor
99e5fe42 44}
45
46////////////////////////////////////////////////////////////////////////
47
48AliUA1JetFinder::~AliUA1JetFinder()
49
50{
99e5fe42 51 // destructor
99e5fe42 52}
53
54////////////////////////////////////////////////////////////////////////
55
56#ifndef WIN32
57# define ua1_jet_finder ua1_jet_finder_
58# define hf1 hf1_
59# define type_of_call
60
61#else
62# define ua1_jet_finder UA1_JET_FINDER
63# define hf1 HF1
64# define type_of_call _stdcall
65#endif
66
67extern "C" void type_of_call
68ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
69 Float_t etc[60000], Float_t etac[60000],
70 Float_t phic[60000],
71 Float_t& min_move, Float_t& max_move, Int_t& mode,
72 Float_t& prec_bg, Int_t& ierror);
73
74extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
75
76////////////////////////////////////////////////////////////////////////
77
b45b0c92 78void AliUA1JetFinder::FindJetsTPC()
99e5fe42 79
80{
81 // initialize event, look for jets, download jet info
82
99e5fe42 83 // initialization in 2 steps
84 // 1) transform input to pt,eta,phi plus lego
85 // 2) dump lego
86
87 // 1) transform input to pt,eta,phi plus lego
88 TClonesArray *lvArray = fReader->GetMomentumArray();
89 Int_t nIn = lvArray->GetEntries();
90 if (nIn == 0) return;
91
92 // local arrays for input
83a444b1 93 Float_t* enT = new Float_t[nIn];
99e5fe42 94 Float_t* ptT = new Float_t[nIn];
95 Float_t* etaT = new Float_t[nIn];
96 Float_t* phiT = new Float_t[nIn];
97
98 // load input vectors
99 for (Int_t i = 0; i < nIn; i++){
83a444b1 100 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
101 enT[i] = lv->Energy();
102 ptT[i] = lv->Pt();
103 etaT[i] = lv->Eta();
104 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
105 if (fReader->GetCutFlag(i) == 1)
99e5fe42 106 fLego->Fill(etaT[i], phiT[i], ptT[i]);
107 }
108 fJets->SetNinput(nIn);
109
110 // 2) dump lego
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
115
116 Int_t nCell = 0;
117 TAxis* xaxis = fLego->GetXaxis();
118 TAxis* yaxis = fLego->GetYaxis();
119 Float_t e = 0.0;
83a444b1 120 for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
121 for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
99e5fe42 122 e = fLego->GetBinContent(i,j);
123 if (e > 0.0) e -= fHeader->GetMinCellEt();
124 if (e < 0.0) e = 0.;
125 Float_t eta = xaxis->GetBinCenter(i);
126 Float_t phi = yaxis->GetBinCenter(j);
127 etCell[nCell] = e;
128 etaCell[nCell] = eta;
129 phiCell[nCell] = phi;
130 nCell++;
131 }
132 }
133
134 // run the algo. Parameters from header
83a444b1 135 Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
99e5fe42 136 Float_t minmove = fHeader->GetMinMove();
137 Float_t maxmove = fHeader->GetMaxMove();
138 Int_t mode = fHeader->GetMode();
139 Float_t precbg = fHeader->GetPrecBg();
140 Int_t ierror;
141 ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
142 minmove, maxmove, mode, precbg, ierror);
143
99e5fe42 144 // sort jets
145 Int_t * idx = new Int_t[UA1JETS.njet];
146 TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
8011d399 147
148 // download jet info.
8011d399 149 for(Int_t i = 0; i < UA1JETS.njet; i++) {
150 // reject events outside acceptable eta range
83a444b1 151 if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
152 || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
153 continue;
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);
8011d399 161 }
99e5fe42 162
8011d399 163 // find multiplicities and relationship jet-particle
164 // find also percentage of pt from pythia
99e5fe42 165 Int_t* injet = new Int_t[nIn];
83a444b1 166 Int_t* sflag = new Int_t[nIn];
167 for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
8011d399 168 Int_t* mult = new Int_t[fJets->GetNJets()];
83a444b1 169 Int_t* ncell = new Int_t[fJets->GetNJets()];
8011d399 170 Float_t* percentage = new Float_t[fJets->GetNJets()];
171
172 for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
83a444b1 173 Float_t pt_sig = 0.0;
174 mult[i] = 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) {
183 injet[j] = -(i+1);
184 if ((fReader->GetCutFlag(j)) == 1 &&
185 (etaT[j] < fHeader->GetLegoEtaMax()) &&
186 (etaT[j] > fHeader->GetLegoEtaMin())) {
187 injet[j] = i+1;
188 mult[i]++;
189 if (fReader->GetSignalFlag(j) == 1) {
190 pt_sig+=ptT[j];
191 sflag[j]=1;
192 }
193 }
194 }
99e5fe42 195 }
83a444b1 196 percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
197 ((Double_t) fJets->GetPt(i));
99e5fe42 198 }
83a444b1 199 fJets->SetNCells(ncell);
8011d399 200 fJets->SetPtFromSignal(percentage);
99e5fe42 201 fJets->SetMultiplicities(mult);
202 fJets->SetInJet(injet);
83a444b1 203 fJets->SetEtaIn(etaT);
204 fJets->SetPhiIn(phiT);
205 fJets->SetPtIn(ptT);
206 fJets->SetEtAvg(UA1JETS.etavg);
207 delete ncell;
208 delete enT;
209 delete ptT;
210 delete etaT;
211 delete phiT;
212 delete injet;
213 delete idx;
214 delete mult;
215 delete percentage;
99e5fe42 216}
217
b45b0c92 218void AliUA1JetFinder::FindJets()
219{ // Find jets with TPC or EMCAL or TPC + EMCAL information
220 // initialize event, look for jets, download jet info
221
222 // 1) transform input to pt,eta,phi plus lego
223 AliJetUnitArray *fUnit = fReader->GetUnitArray();
224 Int_t nIn = fUnit->GetUnitEntries();
225 Int_t fOption = fReader->GetReaderHeader()->GetDetector();
226 Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
227
228 if(fDebug>1){
229 printf("Inside FindJets function ! \n");
230 printf("fOption : %d \n",fOption);
231 }
232
233 if(fDebug>10){
234 cout << "fUnit : " << fUnit << endl;
235 printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
236 }
237
238 if(fDebug>10){
239 printf(" -----------------------------------------------------------------\n");
240 printf(" --- All inputs in fUnitArray ---\n");
241 for(Int_t i=0; i<nIn ; i++){
242 if(fUnit[i].GetUnitEnergy()!=0){
243 printf(" -----------------------------------------------------------------\n");
244 cout << " | ID : " << fUnit[i].GetUnitID() << " | Eta : " << fUnit[i].GetUnitEta() << " | Phi : " << fUnit[i].GetUnitPhi() << " | Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
245 }
246 }
247 }
248
249 if (nIn == 0) return;
250
251 Int_t nCandidateCut = 0;
252 Int_t nCandidateCut2 = 0;
253 Int_t nCandidate = 0;
254 for (Int_t i = 0; i < nIn; i++){
255 if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
256 // Number of good tracks/digits which have passed the pt cut
257 nCandidateCut += 1;
258 }
259 if(fUnit[i].GetUnitEnergy()>0.) {
260 // Number of good tracks/digits without pt cut
261 nCandidate += 1;
262 }
263 }
264
265 if(fDebug>=3){
266 cout << "nCandidate : " << nCandidate << endl;
267 cout << "nCandidateCut : " << nCandidateCut << endl;
268 cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
269 }
270
271 // local arrays for input No Cuts
272 // Both pt < ptMin and pt > ptMin
273 Float_t* enT = new Float_t[nCandidate];
274 Float_t* ptT = new Float_t[nCandidate];
275 Float_t* etaT = new Float_t[nCandidate];
276 Float_t* phiT = new Float_t[nCandidate];
277 Float_t* detaT = new Float_t[nCandidate];
278 Float_t* dphiT = new Float_t[nCandidate];
279 Float_t* cFlagT = new Float_t[nCandidate];
280 Float_t* cClusterT = new Float_t[nCandidate];
281 Float_t* idT = new Float_t[nCandidate];
282 Int_t loop1 = 0;
283
284 fJets->SetNinput(nCandidate);
285
286 if(fDebug>3){
287 cout << "nCandidate : " << nCandidate << endl;
288 // cout << "nMultCandidate : " << nMultCandidate << endl;
289 }
290
291 Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
292 Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
293 Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
294
295 // Information extracted from fUnitArray
296 for(Int_t i=0; i<nIn; i++)
297 {
298 if(fUnit[i].GetUnitCutFlag()==1){
299 etCell[i] = fUnit[i].GetUnitEnergy();
300 if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
301 if (etCell[i] < 0.0) etCell[i] = 0.;
302 etaCell[i] = fUnit[i].GetUnitEta();
303 phiCell[i] = fUnit[i].GetUnitPhi();
304 }
305 else {
306 etCell[i] = 0.;
307 etaCell[i] = fUnit[i].GetUnitEta();
308 phiCell[i] = fUnit[i].GetUnitPhi();
309 }
310
311 if(fUnit[i].GetUnitEnergy()>0.){
312 ptT[loop1] = fUnit[i].GetUnitEnergy();
313 enT[loop1] = fUnit[i].GetUnitEnergy();
314 etaT[loop1] = fUnit[i].GetUnitEta();
315 phiT[loop1] = fUnit[i].GetUnitPhi();
316 detaT[loop1] = fUnit[i].GetUnitDeta();
317 dphiT[loop1] = fUnit[i].GetUnitDphi();
318 cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
319 idT[loop1] = fUnit[i].GetUnitID();
320 loop1++;
321 }
322 }
323
324
325 if(fDebug > 40) // For comparison
326 {
327 for(Int_t j=0; j<nIn; j++) {
328 if(etCell[j]>0){
329 cout << "etCell[" << j << "] : " << etCell[j] << endl;
330 cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
331 cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
332 }
333 }
334 }
335
336 // Run the algo. Parameters from header
337 // Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
338 Int_t nTot = nIn;
339 Float_t minmove = fHeader->GetMinMove();
340 Float_t maxmove = fHeader->GetMaxMove();
341 Int_t mode = fHeader->GetMode();
342 Float_t precbg = fHeader->GetPrecBg();
343 Int_t ierror;
344
345 ua1_jet_finder(nIn, nTot, etCell, etaCell, phiCell,
346 minmove, maxmove, mode, precbg, ierror);
347
348 // sort jets
349 Int_t * idx = new Int_t[UA1JETS.njet];
350 TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
351
352 if(fDebug > 20)
353 {
354 for(Int_t i = 0; i < UA1JETS.njet; i++)
355 {
356 cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
357 cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
358 cout << "idx[" << i << "] : " << idx[i] << endl;
359 cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
360 cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
361 cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
362 }
363 }
364
365 // download jet info.
366 for(Int_t i = 0; i < UA1JETS.njet; i++) {
367 // reject events outside acceptable eta range
368 if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
369 || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
370 {
371 continue;
372 }
373
374 Float_t px, py,pz,en,pT; // convert to 4-vector
375 px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
376 py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
377 pz = UA1JETS.etj[idx[i]] /
378 TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
379 en = TMath::Sqrt(px * px + py * py + pz * pz);
380 pT = TMath::Sqrt(px * px + py * py);
381
382 fJets->AddJet(px, py, pz, en);
383
384 }
385
386 // find multiplicities and relationship jet-particle
387 // find also percentage of pt from pythia
388
389 Int_t* injet = new Int_t[nCandidate];
390 Int_t* sflag = new Int_t[nCandidate];
391 for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
392 Int_t* mult = new Int_t[fJets->GetNJets()];
393 Int_t* ncell = new Int_t[fJets->GetNJets()];
394 Float_t* percentage = new Float_t[fJets->GetNJets()];
395
396 // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
397 // With the particle ID, it is possible to to have access to its physical properties and one can,
398 // for example, set if the corresponding particle is inside or outside the jet with the flag
399 // kOutJet/kInJet, other possibilities...
400
401 for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
402 Float_t pt_sig = 0.0;
403 mult[i] = 0;
404 ncell[i] = UA1JETS.ncellj[i];
405 for (Int_t j = 0; j < nCandidate; j++) {
406 Float_t deta = etaT[j] - fJets->GetEta(i);
407 Float_t dphi = phiT[j] - fJets->GetPhi(i);
408 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
409 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
410 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
411 if (dr < fHeader->GetRadius() && injet[j] == 0) {
412 injet[j] = -(i+1);
413 if(cFlagT[j] == 1 &&
414 (etaT[j] < fHeader->GetLegoEtaMax()) &&
415 (etaT[j] > fHeader->GetLegoEtaMin())) {
416 injet[j] = i+1;
417 mult[i]++;
418 pt_sig+=enT[j];
419 sflag[j]=1;
420 }
421 }
422 if(fDebug>10){
423 cout << "mult[" << i << "] : " << mult[i] << endl;
424 cout << "ncell[" << i << "] : " << ncell[i] << endl;
425 }
426 }
427 percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
428 ((Double_t) fJets->GetPt(i));
429 }
430
431 fJets->SetNCells(ncell);
432 fJets->SetPtFromSignal(percentage);
433 fJets->SetMultiplicities(mult);
434 fJets->SetInJet(injet);
435 fJets->SetEtaIn(etaT);
436 if(fDebug>10){
437 for(Int_t i=0; i<nCandidate ; i++){
438 cout << "phiT[" << i << "] : " << phiT[i] << endl;
439 cout << "etaT[" << i << "] : " << etaT[i] << endl;
440 }
441 }
442 fJets->SetPhiIn(phiT);
443 fJets->SetPtIn(enT);
444 fJets->SetEtAvg(UA1JETS.etavg);
445 delete etCell;
446 delete etaCell;
447 delete phiCell;
448 delete ncell;
449 delete cFlagT;
450 delete cClusterT;
451 delete enT;
452 delete ptT;
453 delete etaT;
454 delete phiT;
455 delete detaT;
456 delete dphiT;
457 delete injet;
458 delete idx;
459 delete mult;
460 delete percentage;
461
462}
463
464
99e5fe42 465////////////////////////////////////////////////////////////////////////
466
467void AliUA1JetFinder::Reset()
99e5fe42 468{
469 fLego->Reset();
470 fJets->ClearJets();
471}
472
473////////////////////////////////////////////////////////////////////////
474
475void AliUA1JetFinder::WriteJHeaderToFile()
476{
477 fOut->cd();
478 fHeader->Write();
479}
480
99e5fe42 481////////////////////////////////////////////////////////////////////////
482
483void AliUA1JetFinder::Init()
484{
485 // initializes some variables
486 Float_t dEta, dPhi;
83a444b1 487 dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
488 /((Float_t) fHeader->GetLegoNbinEta());
489 dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
490 /((Float_t) fHeader->GetLegoNbinPhi());
99e5fe42 491
492 UA1CELL.etaCellSize = dEta;
493 UA1CELL.phiCellSize = dPhi;
494
495 // jet parameters
496 UA1PARA.coneRad = fHeader->GetRadius();
497 UA1PARA.etSeed = fHeader->GetEtSeed();
498 UA1PARA.ejMin = fHeader->GetMinJetEt();
499 UA1PARA.etMin = fHeader->GetMinCellEt();
500
501 // book lego
502 fLego = new
503 TH2F("legoH","eta-phi",
83a444b1 504 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
505 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
506 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
99e5fe42 507}