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