List of tracks added.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinder.cxx
CommitLineData
471f69dc 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 **************************************************************************/
15
16/* $Id$ */
17//*-- Author: Andreas Morsch (CERN)
18#include <TClonesArray.h>
19#include <TTree.h>
20#include <TFile.h>
21#include <TH2.h>
22#include <TAxis.h>
23#include <TParticle.h>
24#include <TParticlePDG.h>
25#include "AliEMCALJetFinder.h"
26#include "AliEMCALGeometry.h"
27#include "AliEMCALHit.h"
28#include "Ecommon.h"
29#include "AliRun.h"
30#include "AliEMCAL.h"
31#include "AliHeader.h"
32
33ClassImp(AliEMCALJetFinder)
34
35//____________________________________________________________________________
36AliEMCALJetFinder::AliEMCALJetFinder()
37{
38// Default constructor
39 fJets = 0;
40 fNjets = 0;
41 fLego = 0;
42}
43
44AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
45 : TTask(name, title)
46{
47// Constructor
48 fJets = new TClonesArray("AliEMCALJet",10000);
49 fNjets = 0;
50 for (Int_t i=0; i<30000; i++)
51 {
52 fEtCell[i] = 0.;
53 fEtaCell[i] = 0.;
54 fPhiCell[i] = 0.;
55 }
56 fLego = 0;
57}
58
59
60
61
62//____________________________________________________________________________
63AliEMCALJetFinder::~AliEMCALJetFinder()
64{
65// Destructor
66//
67 if (fJets){
68 fJets->Delete();
69 delete fJets;
70 }
71}
72
73
74
75#ifndef WIN32
76# define jet_finder_ua1 jet_finder_ua1_
77# define hf1 hf1_
78# define type_of_call
79
80#else
81# define jet_finder_ua1 J
82# define hf1 HF1
83# define type_of_call _stdcall
84#endif
85
86extern "C" void type_of_call
87jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
88 Float_t etc[30000], Float_t etac[30000],
89 Float_t phic[30000],
90 Float_t& min_move, Float_t& max_move, Int_t& mode,
91 Float_t& prec_bg, Int_t& ierror);
92
93extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
94
95
96
97
98
99void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
100 Float_t etac[30000], Float_t phic[30000],
101 Float_t min_move, Float_t max_move, Int_t mode,
102 Float_t prec_bg, Int_t ierror)
103{
104// Wrapper for fortran coded jet finder
105// Full argument list
106 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
107 min_move, max_move, mode, prec_bg, ierror);
108 // Write result to output
109 WriteJets();
110}
111
112void AliEMCALJetFinder::Find()
113{
114// Wrapper for fortran coded jet finder using member data for
115// argument list
116
117 Float_t min_move = 0;
118 Float_t max_move = 0;
119 Int_t mode = 0;
120 Float_t prec_bg = 0.;
121 Int_t ierror = 0;
122
123 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
124 min_move, max_move, mode, prec_bg, ierror);
125 // Write result to output
126 WriteJets();
127}
128
129
130Int_t AliEMCALJetFinder::Njets()
131{
132// Get number of reconstructed jets
133 return EMCALJETS.njet;
134}
135
136Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
137{
138// Get reconstructed jet energy
139 return EMCALJETS.etj[i];
140}
141
142Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
143{
144// Get reconstructed jet phi from leading particle
145 return EMCALJETS.phij[0][i];
146}
147
148Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
149{
150// Get reconstructed jet phi from weighting
151 return EMCALJETS.phij[1][i];
152}
153
154Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
155{
156// Get reconstructed jet eta from leading particles
157 return EMCALJETS.etaj[0][i];
158}
159
160
161Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
162{
163// Get reconstructed jet eta from weighting
164 return EMCALJETS.etaj[1][i];
165}
166
167void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
168{
169// Set grid cell size
170 EMCALCELLGEO.etaCellSize = eta;
171 EMCALCELLGEO.phiCellSize = phi;
172}
173
174void AliEMCALJetFinder::SetConeRadius(Float_t par)
175{
176// Set jet cone radius
177 EMCALJETPARAM.coneRad = par;
178}
179
180void AliEMCALJetFinder::SetEtSeed(Float_t par)
181{
182// Set et cut for seeds
183 EMCALJETPARAM.etSeed = par;
184}
185
186void AliEMCALJetFinder::SetMinJetEt(Float_t par)
187{
188// Set minimum jet et
189 EMCALJETPARAM.ejMin = par;
190}
191
192void AliEMCALJetFinder::SetMinCellEt(Float_t par)
193{
194// Set et cut per cell
195 EMCALJETPARAM.etMin = par;
196}
197
198
199void AliEMCALJetFinder::Test()
200{
201//
202// Test the finder call
203//
204 const Int_t nmax = 30000;
205 Int_t ncell = 10;
206 Int_t ncell_tot = 100;
207
208 Float_t etc[nmax];
209 Float_t etac[nmax];
210 Float_t phic[nmax];
211 Float_t min_move = 0;
212 Float_t max_move = 0;
213 Int_t mode = 0;
214 Float_t prec_bg = 0;
215 Int_t ierror = 0;
216
217
218 Find(ncell, ncell_tot, etc, etac, phic,
219 min_move, max_move, mode, prec_bg, ierror);
220
221}
222
223//
224// I/O
225//
226
227void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
228{
229 //
230 // Add a jet
231 //
232 TClonesArray &lrawcl = *fJets;
233 new(lrawcl[fNjets++]) AliEMCALJet(jet);
234}
235
236void AliEMCALJetFinder::ResetJets()
237{
238 //
239 // Reset Jet List
240 //
241 fJets->Clear();
242 fNjets = 0;
243}
244
245void AliEMCALJetFinder::WriteJets()
246{
247//
248// Add all jets to the list
249//
250 const Int_t kBufferSize = 4000;
251 TTree *pK = gAlice->TreeK();
252 const char* file = (pK->GetCurrentFile())->GetName();
253// I/O
254 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
255 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
256 if (fJets && gAlice->TreeR()) {
257 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
258 "Jets",
259 &fJets,
260 kBufferSize,
261 file);
262 }
263 Int_t njet = Njets();
264 for (Int_t nj=0; nj<njet; nj++)
265 {
266 AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
267 JetPhiW(nj),
268 JetEtaW(nj));
269 AddJet(*jet);
270 delete jet;
271 }
272 Int_t nev = gAlice->GetHeader()->GetEvent();
273 gAlice->TreeR()->Fill();
274 char hname[30];
275 sprintf(hname,"TreeR%d", nev);
276 gAlice->TreeR()->Write(hname);
277 gAlice->TreeR()->Reset();
278 ResetJets();
279}
280
281void AliEMCALJetFinder::BookLego()
282{
283//
284// Book histo for discretisation
285//
286//
287// Get geometry parameters from
288 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
289 AliEMCALGeometry* geom =
290 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
291 fNbinEta = geom->GetNZ();
292 fNbinPhi = geom->GetNPhi();
293 const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
294 const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
295 fDphi = (phiMax-phiMin)/fNbinEta;
296 fDeta = 1.4/fNbinEta;
297 fNtot = fNbinPhi*fNbinEta;
298
299//
300 fLego = new TH2F("legoH","eta-phi",
301 fNbinEta, -0.7, 0.7,
302 fNbinPhi, phiMin, phiMax);
303}
304
305void AliEMCALJetFinder::DumpLego()
306{
307//
308// Dump lego histo into array
309//
310 fNcell = 0;
311 for (Int_t i = 1; i < fNbinEta; i++) {
312 for (Int_t j = 1; j < fNbinPhi; j++) {
313 Float_t e = fLego->GetBinContent(i,j);
314 TAxis* Xaxis = fLego->GetXaxis();
315 TAxis* Yaxis = fLego->GetYaxis();
316 Float_t eta = Xaxis->GetBinCenter(i);
317 Float_t phi = Yaxis->GetBinCenter(j);
318 fEtCell[fNcell] = e;
319 fEtaCell[fNcell] = eta;
320 fPhiCell[fNcell] = phi;
321 fNcell++;
322 } // phi
323 } // eta
324 fNcell--;
325}
326
327void AliEMCALJetFinder::ResetMap()
328{
329//
330// Reset eta-phi array
331
332 for (Int_t i=0; i<30000; i++)
333 {
334 fEtCell[i] = 0.;
335 fEtaCell[i] = 0.;
336 fPhiCell[i] = 0.;
337 }
338}
339
340
341void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
342{
343//
344// Fill Cells with hit information
345//
346//
347 ResetMap();
348
349 if (!fLego) BookLego();
350// Reset
351 if (flag == 0) fLego->Reset();
352//
353// Access particle information
354 Int_t npart = (gAlice->GetHeader())->GetNprimary();
355 for (Int_t part=2; part<npart; part++) {
356 TParticle *MPart = gAlice->Particle(part);
357 Int_t mpart = MPart->GetPdgCode();
358 Int_t child1 = MPart->GetFirstDaughter();
359 Float_t pT = MPart->Pt();
360 Float_t phi = MPart->Phi();
361 Float_t theta = MPart->Theta();
362 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
363// if (part == 6 || part == 7)
364// {
365// printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
366// part-5, pT, eta, phi);
367// }
368
369 if (pT == 0.) continue;
370// charged or neutral
371 if (ich == 0) {
372 TParticlePDG* pdgP = MPart->GetPDG();
373 if (pdgP->Charge() == 0) continue;
374 }
375// skip partons
376 if (TMath::Abs(mpart) <= 6 ||
377 mpart == 21 ||
378 mpart == 92) continue;
379// acceptance cut
380 if (TMath::Abs(eta) > 0.7) continue;
381 if (phi*180./TMath::Pi() > 120.) continue;
382// final state only
383 if (child1 >= 0 && child1 < npart) continue;
384// printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
385// part, mpart, child1, eta, phi, pT);
386 fLego->Fill(eta, phi, pT);
387 } // primary loop
388 DumpLego();
389}
390
391void AliEMCALJetFinder::FillFromHits(Int_t flag)
392{
393//
394// Fill Cells with track information
395//
396//
397 ResetMap();
398
399 if (!fLego) BookLego();
400 if (flag == 0) fLego->Reset();
401
402//
403// Access hit information
404 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
405// AliEMCALGeometry* geom =
406// AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
407
408 TTree *treeH = gAlice->TreeH();
409 Int_t ntracks = (Int_t) treeH->GetEntries();
410//
411// Loop over tracks
412//
413 Int_t nbytes = 0;
414
415
416 for (Int_t track=0; track<ntracks;track++) {
417 gAlice->ResetHits();
418 nbytes += treeH->GetEvent(track);
419//
420// Loop over hits
421//
422 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
423 mHit;
424 mHit=(AliEMCALHit*) pEMCAL->NextHit())
425 {
426 Float_t x = mHit->X(); // x-pos of hit
427 Float_t y = mHit->Y(); // y-pos
428 Float_t z = mHit->Z(); // z-pos
429 Float_t eloss = mHit->GetEnergy(); // deposited energy
430// Int_t index = mHit->GetId(); // cell index
431// Float_t eta, phi;
432// geom->EtaPhiFromIndex(index, eta, phi);
433 Float_t r = TMath::Sqrt(x*x+y*y);
434 Float_t theta = TMath::ATan2(r,z);
435 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
436 Float_t phi = TMath::ATan2(y,x);
437 fLego->Fill(eta, phi, eloss);
438// if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f",
439// r, z, eta, phi, eloss);
440// printf("\n Max %f", fLego->GetMaximum());
441 } // Hit Loop
442 } // Track Loop
443 DumpLego();
444}
445
446
447void hf1(Int_t& id, Float_t& x, Float_t& wgt)
448{
449// dummy for hbook calls
450 ;
451}
452
453