]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALJetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
A macro to test the coordinate transformations of AliITSgeom.
[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
d594b22e 16/*
17$Log$
e4717047 18Revision 1.11 2002/01/22 17:25:47 morsch
19Some corrections for event mixing and bg event handling.
20
610f27c9 21Revision 1.10 2002/01/22 10:31:44 morsch
22Some correction for bg mixing.
23
46955412 24Revision 1.9 2002/01/21 16:28:42 morsch
25Correct sign of dphi.
26
f8f69f71 27Revision 1.8 2002/01/21 16:05:31 morsch
28- different phi-bin for hadron correction
29- provisions for background mixing.
30
ff114de3 31Revision 1.7 2002/01/21 15:47:18 morsch
32Bending radius correctly in cm.
33
b3feface 34Revision 1.6 2002/01/21 12:53:50 morsch
35authors
36
0c207dbc 37Revision 1.5 2002/01/21 12:47:47 morsch
38Possibility to include K0long and neutrons.
39
d53b7b0b 40Revision 1.4 2002/01/21 11:03:21 morsch
41Phi propagation introduced in FillFromTracks.
42
26dadf3a 43Revision 1.3 2002/01/18 05:07:56 morsch
44- hadronic correction
45- filling of digits
46- track selection upon EMCAL information
47
d594b22e 48*/
49
0c207dbc 50//*-- Authors: Andreas Morsch (CERN)
51//* J.L. Klay (LBL)
52//* Aleksei Pavlinov (WSU)
53
d594b22e 54// From root ...
471f69dc 55#include <TClonesArray.h>
56#include <TTree.h>
d594b22e 57#include <TBranchElement.h>
471f69dc 58#include <TFile.h>
59#include <TH2.h>
60#include <TAxis.h>
61#include <TParticle.h>
62#include <TParticlePDG.h>
d594b22e 63
64// From AliRoot ...
471f69dc 65#include "AliEMCALJetFinder.h"
d594b22e 66#include "AliEMCALFast.h"
471f69dc 67#include "AliEMCALGeometry.h"
68#include "AliEMCALHit.h"
d594b22e 69#include "AliEMCALDigit.h"
70#include "AliEMCALDigitizer.h"
71#include "AliEMCALHadronCorrection.h"
471f69dc 72#include "AliRun.h"
26dadf3a 73#include "AliMagF.h"
74#include "AliMagFCM.h"
471f69dc 75#include "AliEMCAL.h"
76#include "AliHeader.h"
d53b7b0b 77#include "AliPDG.h"
471f69dc 78
ff114de3 79// Interface to FORTRAN
80#include "Ecommon.h"
81
82
471f69dc 83ClassImp(AliEMCALJetFinder)
84
85//____________________________________________________________________________
86AliEMCALJetFinder::AliEMCALJetFinder()
87{
88// Default constructor
d594b22e 89 fJets = 0;
90 fNjets = 0;
91 fLego = 0;
ff114de3 92 fLegoB = 0;
d594b22e 93 fTrackList = 0;
94 fPtT = 0;
95 fEtaT = 0;
96 fPhiT = 0;
ff114de3 97 fPtB = 0;
98 fEtaB = 0;
99 fPhiB = 0;
d594b22e 100 fHadronCorrector = 0;
471f69dc 101}
102
103AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
104 : TTask(name, title)
105{
106// Constructor
107 fJets = new TClonesArray("AliEMCALJet",10000);
108 fNjets = 0;
d594b22e 109 for (Int_t i = 0; i < 30000; i++)
471f69dc 110 {
111 fEtCell[i] = 0.;
112 fEtaCell[i] = 0.;
113 fPhiCell[i] = 0.;
114 }
115 fLego = 0;
ff114de3 116 fTrackList = 0;
117 fPtT = 0;
118 fEtaT = 0;
119 fPhiT = 0;
120 fPtB = 0;
121 fEtaB = 0;
122 fPhiB = 0;
d594b22e 123 fHadronCorrector = 0;
ff114de3 124 fBackground = 0;
0c207dbc 125//
d594b22e 126 SetPtCut();
127 SetMomentumSmearing();
128 SetEfficiencySim();
129 SetDebug();
130 SetHadronCorrection();
131 SetSamplingFraction();
d53b7b0b 132 SetIncludeK0andN();
471f69dc 133}
134
135
136
137
138//____________________________________________________________________________
139AliEMCALJetFinder::~AliEMCALJetFinder()
140{
141// Destructor
142//
143 if (fJets){
144 fJets->Delete();
145 delete fJets;
146 }
147}
148
149
150
3bc109a9 151
471f69dc 152#ifndef WIN32
153# define jet_finder_ua1 jet_finder_ua1_
154# define hf1 hf1_
155# define type_of_call
156
157#else
158# define jet_finder_ua1 J
159# define hf1 HF1
160# define type_of_call _stdcall
161#endif
162
163extern "C" void type_of_call
164jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
165 Float_t etc[30000], Float_t etac[30000],
166 Float_t phic[30000],
167 Float_t& min_move, Float_t& max_move, Int_t& mode,
168 Float_t& prec_bg, Int_t& ierror);
169
170extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
171
172
173
174
175
176void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
177 Float_t etac[30000], Float_t phic[30000],
178 Float_t min_move, Float_t max_move, Int_t mode,
179 Float_t prec_bg, Int_t ierror)
180{
181// Wrapper for fortran coded jet finder
182// Full argument list
183 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
184 min_move, max_move, mode, prec_bg, ierror);
185 // Write result to output
186 WriteJets();
187}
188
189void AliEMCALJetFinder::Find()
190{
191// Wrapper for fortran coded jet finder using member data for
192// argument list
193
194 Float_t min_move = 0;
195 Float_t max_move = 0;
196 Int_t mode = 0;
197 Float_t prec_bg = 0.;
198 Int_t ierror = 0;
199
200 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
201 min_move, max_move, mode, prec_bg, ierror);
202 // Write result to output
3bc109a9 203 Int_t njet = Njets();
d594b22e 204
3bc109a9 205 for (Int_t nj=0; nj<njet; nj++)
206 {
d594b22e 207
208
3bc109a9 209 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
210 JetPhiW(nj),
211 JetEtaW(nj));
212 }
213 FindTracksInJetCone();
471f69dc 214 WriteJets();
215}
216
217
218Int_t AliEMCALJetFinder::Njets()
219{
220// Get number of reconstructed jets
221 return EMCALJETS.njet;
222}
223
224Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
225{
226// Get reconstructed jet energy
227 return EMCALJETS.etj[i];
228}
229
230Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
231{
232// Get reconstructed jet phi from leading particle
233 return EMCALJETS.phij[0][i];
234}
235
236Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
237{
238// Get reconstructed jet phi from weighting
239 return EMCALJETS.phij[1][i];
240}
241
242Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
243{
244// Get reconstructed jet eta from leading particles
245 return EMCALJETS.etaj[0][i];
246}
247
248
249Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
250{
251// Get reconstructed jet eta from weighting
252 return EMCALJETS.etaj[1][i];
253}
254
255void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
256{
257// Set grid cell size
258 EMCALCELLGEO.etaCellSize = eta;
259 EMCALCELLGEO.phiCellSize = phi;
d594b22e 260 fDeta = eta;
261 fDphi = phi;
471f69dc 262}
263
264void AliEMCALJetFinder::SetConeRadius(Float_t par)
265{
266// Set jet cone radius
267 EMCALJETPARAM.coneRad = par;
3bc109a9 268 fConeRadius = par;
471f69dc 269}
270
271void AliEMCALJetFinder::SetEtSeed(Float_t par)
272{
273// Set et cut for seeds
274 EMCALJETPARAM.etSeed = par;
d594b22e 275 fEtSeed = par;
471f69dc 276}
277
278void AliEMCALJetFinder::SetMinJetEt(Float_t par)
279{
280// Set minimum jet et
281 EMCALJETPARAM.ejMin = par;
d594b22e 282 fMinJetEt = par;
471f69dc 283}
284
285void AliEMCALJetFinder::SetMinCellEt(Float_t par)
286{
287// Set et cut per cell
288 EMCALJETPARAM.etMin = par;
d594b22e 289 fMinCellEt = par;
471f69dc 290}
291
3bc109a9 292void AliEMCALJetFinder::SetPtCut(Float_t par)
293{
294// Set pT cut on charged tracks
295 fPtCut = par;
296}
297
471f69dc 298
299void AliEMCALJetFinder::Test()
300{
301//
302// Test the finder call
303//
304 const Int_t nmax = 30000;
305 Int_t ncell = 10;
306 Int_t ncell_tot = 100;
307
308 Float_t etc[nmax];
309 Float_t etac[nmax];
310 Float_t phic[nmax];
311 Float_t min_move = 0;
312 Float_t max_move = 0;
313 Int_t mode = 0;
314 Float_t prec_bg = 0;
315 Int_t ierror = 0;
316
317
318 Find(ncell, ncell_tot, etc, etac, phic,
319 min_move, max_move, mode, prec_bg, ierror);
320
321}
322
323//
324// I/O
325//
326
327void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
328{
329 //
330 // Add a jet
331 //
332 TClonesArray &lrawcl = *fJets;
333 new(lrawcl[fNjets++]) AliEMCALJet(jet);
334}
335
336void AliEMCALJetFinder::ResetJets()
337{
338 //
339 // Reset Jet List
340 //
341 fJets->Clear();
342 fNjets = 0;
343}
344
345void AliEMCALJetFinder::WriteJets()
346{
347//
348// Add all jets to the list
349//
350 const Int_t kBufferSize = 4000;
351 TTree *pK = gAlice->TreeK();
352 const char* file = (pK->GetCurrentFile())->GetName();
353// I/O
354 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
d594b22e 355
356 if (fDebug > 1)
471f69dc 357 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
d594b22e 358
471f69dc 359 if (fJets && gAlice->TreeR()) {
360 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
d594b22e 361 "EMCALJets",
471f69dc 362 &fJets,
363 kBufferSize,
364 file);
365 }
366 Int_t njet = Njets();
3bc109a9 367 for (Int_t nj = 0; nj < njet; nj++)
471f69dc 368 {
3bc109a9 369 AddJet(*fJetT[nj]);
370 delete fJetT[nj];
471f69dc 371 }
3bc109a9 372
471f69dc 373 Int_t nev = gAlice->GetHeader()->GetEvent();
374 gAlice->TreeR()->Fill();
375 char hname[30];
376 sprintf(hname,"TreeR%d", nev);
377 gAlice->TreeR()->Write(hname);
378 gAlice->TreeR()->Reset();
379 ResetJets();
380}
381
382void AliEMCALJetFinder::BookLego()
383{
384//
385// Book histo for discretisation
386//
387//
388// Get geometry parameters from
389 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
390 AliEMCALGeometry* geom =
391 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
392 fNbinEta = geom->GetNZ();
393 fNbinPhi = geom->GetNPhi();
394 const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
395 const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
610f27c9 396 const Float_t etaMin = geom->GetArm1EtaMin();
397 const Float_t etaMax = geom->GetArm1EtaMax();
471f69dc 398 fDphi = (phiMax-phiMin)/fNbinEta;
399 fDeta = 1.4/fNbinEta;
400 fNtot = fNbinPhi*fNbinEta;
610f27c9 401//
402// Don't add histos to the current directory
403 TH2::AddDirectory(0);
471f69dc 404//
46955412 405// Signal map
471f69dc 406 fLego = new TH2F("legoH","eta-phi",
46955412 407 fNbinEta, etaMin, etaMax,
408 fNbinPhi, phiMin, phiMax);
409//
410// Background map
411 fLegoB = new TH2F("legoB","eta-phi",
412 fNbinEta, etaMin, etaMax,
471f69dc 413 fNbinPhi, phiMin, phiMax);
610f27c9 414
471f69dc 415}
416
417void AliEMCALJetFinder::DumpLego()
418{
419//
420// Dump lego histo into array
421//
422 fNcell = 0;
423 for (Int_t i = 1; i < fNbinEta; i++) {
424 for (Int_t j = 1; j < fNbinPhi; j++) {
425 Float_t e = fLego->GetBinContent(i,j);
46955412 426 if (e <=0.) continue;
427
471f69dc 428 TAxis* Xaxis = fLego->GetXaxis();
429 TAxis* Yaxis = fLego->GetYaxis();
430 Float_t eta = Xaxis->GetBinCenter(i);
431 Float_t phi = Yaxis->GetBinCenter(j);
432 fEtCell[fNcell] = e;
433 fEtaCell[fNcell] = eta;
434 fPhiCell[fNcell] = phi;
435 fNcell++;
436 } // phi
437 } // eta
438 fNcell--;
439}
440
441void AliEMCALJetFinder::ResetMap()
442{
443//
444// Reset eta-phi array
445
446 for (Int_t i=0; i<30000; i++)
447 {
448 fEtCell[i] = 0.;
449 fEtaCell[i] = 0.;
450 fPhiCell[i] = 0.;
451 }
452}
453
454
455void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
456{
457//
d594b22e 458// Fill Cells with track information
471f69dc 459//
460//
461 ResetMap();
462
463 if (!fLego) BookLego();
464// Reset
465 if (flag == 0) fLego->Reset();
466//
467// Access particle information
468 Int_t npart = (gAlice->GetHeader())->GetNprimary();
3bc109a9 469// Create track list
470//
471// 0: not selected
472// 1: selected for jet finding
473// 2: inside jet
474// ....
475 if (fTrackList) delete[] fTrackList;
476 if (fPtT) delete[] fPtT;
477 if (fEtaT) delete[] fEtaT;
478 if (fPhiT) delete[] fPhiT;
479
480 fTrackList = new Int_t [npart];
481 fPtT = new Float_t[npart];
482 fEtaT = new Float_t[npart];
483 fPhiT = new Float_t[npart];
ff114de3 484
3bc109a9 485 fNt = npart;
ff114de3 486 fNtS = 0;
487
d594b22e 488 for (Int_t part = 2; part < npart; part++) {
471f69dc 489 TParticle *MPart = gAlice->Particle(part);
490 Int_t mpart = MPart->GetPdgCode();
491 Int_t child1 = MPart->GetFirstDaughter();
492 Float_t pT = MPart->Pt();
d594b22e 493 Float_t p = MPart->P();
471f69dc 494 Float_t phi = MPart->Phi();
3bc109a9 495 Float_t eta = MPart->Eta();
46955412 496
d594b22e 497 if (fDebug > 0) {
498 if (part == 6 || part == 7)
499 {
610f27c9 500 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
d594b22e 501 part-5, pT, eta, phi);
502 }
503
504// if (mpart == 21)
505
506// printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
507// part, mpart, pT, eta, phi);
508 }
509
3bc109a9 510
511 fTrackList[part] = 0;
512 fPtT[part] = pT;
513 fEtaT[part] = eta;
514 fPhiT[part] = phi;
515
516 if (part < 2) continue;
517 if (pT == 0 || pT < fPtCut) continue;
26dadf3a 518 TParticlePDG* pdgP = 0;
471f69dc 519// charged or neutral
520 if (ich == 0) {
26dadf3a 521 pdgP = MPart->GetPDG();
d53b7b0b 522 if (pdgP->Charge() == 0) {
523 if (!fK0N) {
524 continue;
525 } else {
526 if (mpart != kNeutron &&
527 mpart != kNeutronBar &&
528 mpart != kK0Long) continue;
529 }
530 }
471f69dc 531 }
532// skip partons
533 if (TMath::Abs(mpart) <= 6 ||
534 mpart == 21 ||
535 mpart == 92) continue;
536// acceptance cut
537 if (TMath::Abs(eta) > 0.7) continue;
538 if (phi*180./TMath::Pi() > 120.) continue;
539// final state only
540 if (child1 >= 0 && child1 < npart) continue;
d594b22e 541//
542 if (fDebug > 1)
543 printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
544 part, mpart, child1, eta, phi, pT);
545//
26dadf3a 546//
ff114de3 547// phi propagation for hadronic correction
26dadf3a 548
549 Bool_t curls = kFALSE;
550 Float_t dphi = PropagatePhi(pT, pdgP->Charge(), curls);
610f27c9 551 if (fDebug > 1) printf("\n Delta phi %f pT%f", dphi, pT);
552 if (curls && fDebug > 1) printf("\n Track is curling %f", pT);
ff114de3 553 Float_t phiHC = phi + dphi;
26dadf3a 554//
d594b22e 555// Momentum smearing goes here ...
556//
557 if (fSmear) {
558 p = AliEMCALFast::SmearMomentum(1,p);
559 }
560//
561// Tracking Efficiency and TPC acceptance goes here ...
562 if (fEffic) {
563 Float_t eff = AliEMCALFast::Efficiency(1,p);
564 if (AliEMCALFast::RandomReject(eff)) continue;
565 }
566//
567// Correction of Hadronic Energy goes here ...
568//
ff114de3 569 Float_t dpH = 0.;
570
d594b22e 571 if (fHCorrection) {
572 if (!fHadronCorrector)
573 Fatal("AliEMCALJetFinder",
574 "Hadronic energy correction required but not defined !");
ff114de3 575 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
d594b22e 576 }
577//
578// More to do ? Just think about it !
579//
3bc109a9 580 fTrackList[part] = 1;
ff114de3 581//
582 fNtS++;
583
d594b22e 584 fLego->Fill(eta, phi, p);
610f27c9 585 if (fHCorrection && !curls) fLego->Fill(eta, phiHC, -dpH);
471f69dc 586 } // primary loop
587 DumpLego();
588}
589
590void AliEMCALJetFinder::FillFromHits(Int_t flag)
591{
592//
d594b22e 593// Fill Cells with hit information
471f69dc 594//
595//
596 ResetMap();
597
598 if (!fLego) BookLego();
ff114de3 599// Reset eta-phi map if needed
600 if (flag == 0) fLego->Reset();
601// Initialize from background event if available
471f69dc 602//
603// Access hit information
604 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
471f69dc 605
606 TTree *treeH = gAlice->TreeH();
607 Int_t ntracks = (Int_t) treeH->GetEntries();
608//
609// Loop over tracks
610//
611 Int_t nbytes = 0;
612
613
614 for (Int_t track=0; track<ntracks;track++) {
615 gAlice->ResetHits();
616 nbytes += treeH->GetEvent(track);
617//
618// Loop over hits
619//
620 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
621 mHit;
622 mHit=(AliEMCALHit*) pEMCAL->NextHit())
623 {
624 Float_t x = mHit->X(); // x-pos of hit
625 Float_t y = mHit->Y(); // y-pos
626 Float_t z = mHit->Z(); // z-pos
627 Float_t eloss = mHit->GetEnergy(); // deposited energy
d594b22e 628
629 Float_t r = TMath::Sqrt(x*x+y*y);
630 Float_t theta = TMath::ATan2(r,z);
631 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
632 Float_t phi = TMath::ATan2(y,x);
610f27c9 633
634 if (fDebug > 1) printf("\n Hit %f %f %f %f", x, y, z, eloss);
635
d594b22e 636 fLego->Fill(eta, phi, fSamplingF*eloss);
471f69dc 637 } // Hit Loop
638 } // Track Loop
639 DumpLego();
640}
641
d594b22e 642void AliEMCALJetFinder::FillFromDigits(Int_t flag)
643{
644//
645// Fill Cells with digit information
646//
647//
648 ResetMap();
649
650 if (!fLego) BookLego();
651 if (flag == 0) fLego->Reset();
652 Int_t nbytes;
653
654
655//
656// Connect digits
657//
658 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
659 TTree *treeD = gAlice->TreeD();
660 TBranchElement* branchDg = (TBranchElement*)
661 treeD->GetBranch("EMCAL");
662
663 if (!branchDg) Fatal("AliEMCALJetFinder",
664 "Reading digits requested but no digits in file !");
665
666 branchDg->SetAddress(&digs);
667 Int_t nent = (Int_t) branchDg->GetEntries();
668//
669// Connect digitizer
670//
671 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
672 TBranchElement* branchDr = (TBranchElement*)
673 treeD->GetBranch("AliEMCALDigitizer");
674 branchDr->SetAddress(&digr);
675//
676//
677 nbytes += branchDg->GetEntry(0);
678 nbytes += branchDr->GetEntry(0);
679//
680// Get digitizer parameters
681 Float_t towerADCped = digr->GetTowerpedestal();
682 Float_t towerADCcha = digr->GetTowerchannel();
683 Float_t preshoADCped = digr->GetPreShopedestal();
684 Float_t preshoADCcha = digr->GetPreShochannel();
685
686 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
687 AliEMCALGeometry* geom =
688 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
689
690 if (fDebug) {
691 Int_t ndig = digs->GetEntries();
692 printf("\n Number of Digits: %d %d\n", ndig, nent);
693 printf("\n Parameters: %f %f %f %f\n",
694 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
695 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
696 }
697
698//
699// Loop over digits
700 AliEMCALDigit* sdg;
701 TIter next(digs);
702 while ((sdg = (AliEMCALDigit*)(next())))
703 {
704 Double_t pedestal;
705 Double_t channel;
706 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
707 {
708 pedestal = preshoADCped;
709 channel = preshoADCcha;
710 } else {
711 pedestal = towerADCped;
712 channel = towerADCcha;
713 }
714
715 Float_t eta = sdg->GetEta();
716 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
717 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
718
719 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
720 eta, phi, amp, sdg->GetAmp());
721
722 fLego->Fill(eta, phi, fSamplingF*amp);
723 } // digit loop
724//
725// Dump lego hist
726 DumpLego();
727}
728
729
730void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
731{
732//
733// Fill Cells with hit information
734//
735//
736 ResetMap();
737
738 if (!fLego) BookLego();
739// Reset
740 if (flag == 0) fLego->Reset();
741
ff114de3 742// Flag charged tracks passing through TPC which
743// are associated to EMCAL Hits
d594b22e 744 BuildTrackFlagTable();
745
746//
747// Access particle information
748 TTree *treeH = gAlice->TreeH();
749 Int_t ntracks = (Int_t) treeH->GetEntries();
750
ff114de3 751 if (fPtT) delete[] fPtT;
752 if (fEtaT) delete[] fEtaT;
753 if (fPhiT) delete[] fPhiT;
754
d594b22e 755 fPtT = new Float_t[ntracks];
756 fEtaT = new Float_t[ntracks];
757 fPhiT = new Float_t[ntracks];
d594b22e 758
ff114de3 759 fNt = ntracks;
760 fNtS = 0;
761
d594b22e 762 for (Int_t track = 0; track < ntracks; track++) {
763 TParticle *MPart = gAlice->Particle(track);
764 Float_t pT = MPart->Pt();
765 Float_t phi = MPart->Phi();
766 Float_t eta = MPart->Eta();
767
768 if(fTrackList[track]) {
769 fPtT[track] = pT;
770 fEtaT[track] = eta;
771 fPhiT[track] = phi;
772
773 if (track < 2) continue; //Colliding particles?
774 if (pT == 0 || pT < fPtCut) continue;
ff114de3 775 fNtS++;
d594b22e 776 fLego->Fill(eta, phi, pT);
777 }
778 } // track loop
779 DumpLego();
780}
781
782void AliEMCALJetFinder::BuildTrackFlagTable() {
783
784// Method to generate a lookup table for TreeK particles
785// which are linked to hits in the EMCAL
786//
787// --Author: J.L. Klay
788//
d594b22e 789// Access hit information
790 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
791
792 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
793 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
794
795 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
796 fTrackList = new Int_t[nKTrks]; //before generating a new one
797
798 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
799 fTrackList[i] = 0;
800 }
801
802 TTree *treeH = gAlice->TreeH();
803 Int_t ntracks = (Int_t) treeH->GetEntries();
804//
805// Loop over tracks
806//
807 Int_t nbytes = 0;
808
809 for (Int_t track=0; track<ntracks;track++) {
810 gAlice->ResetHits();
811 nbytes += treeH->GetEvent(track);
812 if (pEMCAL) {
813
814//
815// Loop over hits
816//
817 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
818 mHit;
819 mHit=(AliEMCALHit*) pEMCAL->NextHit())
820 {
821 Int_t iTrk = mHit->Track(); // track number
822 Int_t idprim = mHit->GetPrimary(); // primary particle
823
824 //Determine the origin point of this particle - it made a hit in the EMCAL
825 TParticle *trkPart = gAlice->Particle(iTrk);
826 TParticlePDG *trkPdg = trkPart->GetPDG();
827 Int_t trkCode = trkPart->GetPdgCode();
828 Double_t trkChg;
829 if (trkCode < 10000) { //Big Ions cause problems for
830 trkChg = trkPdg->Charge(); //this function. Since they aren't
831 } else { //likely to make it very far, set
832 trkChg = 0.0; //their charge to 0 for the Flag test
833 }
834 Float_t vxTrk = trkPart->Vx();
835 Float_t vyTrk = trkPart->Vy();
836 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
837 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
838
839 //Loop through the ancestry of the EMCAL entrance particles
840 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
841 while (ancestor != -1) {
842 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
843 TParticlePDG *ancPdg = ancPart->GetPDG();
844 Int_t ancCode = ancPart->GetPdgCode();
845 Double_t ancChg;
846 if (ancCode < 10000) {
847 ancChg = ancPdg->Charge();
848 } else {
849 ancChg = 0.0;
850 }
851 Float_t vxAnc = ancPart->Vx();
852 Float_t vyAnc = ancPart->Vy();
853 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
854 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
855 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
856 }
857
858 //Determine the origin point of the primary particle associated with the hit
859 TParticle *primPart = gAlice->Particle(idprim);
860 TParticlePDG *primPdg = primPart->GetPDG();
861 Int_t primCode = primPart->GetPdgCode();
862 Double_t primChg;
863 if (primCode < 10000) {
864 primChg = primPdg->Charge();
865 } else {
866 primChg = 0.0;
867 }
868 Float_t vxPrim = primPart->Vx();
869 Float_t vyPrim = primPart->Vy();
870 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
871 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
872 } // Hit Loop
873 } //if (pEMCAL)
874 } // Track Loop
875}
876
877Int_t AliEMCALJetFinder
878::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
879
ff114de3 880 Int_t flag = 0;
881 Int_t parton = 0;
d594b22e 882 Int_t neutral = 0;
883
884 if (charge == 0) neutral = 1;
885
886 if (TMath::Abs(code) <= 6 ||
887 code == 21 ||
888 code == 92) parton = 1;
889
890 //It's not a parton, it's charged and it went through the TPC
891 if (!parton && !neutral && radius <= 84.0) flag = 1;
892
893 return flag;
894}
895
ff114de3 896
897
898void AliEMCALJetFinder::SaveBackgroundEvent()
899{
900// Saves the eta-phi lego and the tracklist
901//
46955412 902 if (fLegoB) fLegoB->Reset();
ff114de3 903
46955412 904 fLego->Copy(*fLegoB);
ff114de3 905
906 if (fPtB) delete[] fPtB;
907 if (fEtaB) delete[] fEtaB;
908 if (fPhiB) delete[] fPhiB;
909 if (fTrackListB) delete[] fTrackListB;
910
911 fPtB = new Float_t[fNtS];
912 fEtaB = new Float_t[fNtS];
913 fPhiB = new Float_t[fNtS];
914 fTrackListB = new Int_t [fNtS];
46955412 915
916 fNtB = 0;
917
ff114de3 918 for (Int_t i = 0; i < fNt; i++) {
919 if (!fTrackList[i]) continue;
610f27c9 920 fPtB [fNtB] = fPtT [i];
921 fEtaB[fNtB] = fEtaT[i];
922 fPhiB[fNtB] = fPhiT[i];
46955412 923 fTrackListB[fNtB] = 1;
ff114de3 924 fNtB++;
925 }
610f27c9 926 fBackground = 1;
ff114de3 927}
928
929void AliEMCALJetFinder::InitFromBackground()
930{
931//
932//
46955412 933 if (fDebug) printf("\n Initializing from Background");
934
935 if (fLego) fLego->Reset();
936 fLegoB->Copy(*fLego);
ff114de3 937}
938
939
940
3bc109a9 941void AliEMCALJetFinder::FindTracksInJetCone()
942{
943//
944// Build list of tracks inside jet cone
945//
946// Loop over jets
947 Int_t njet = Njets();
948 for (Int_t nj = 0; nj < njet; nj++)
949 {
950 Float_t etaj = JetEtaW(nj);
951 Float_t phij = JetPhiW(nj);
952 Int_t nT = 0; // number of associated tracks
953
ff114de3 954// Loop over particles in current event
3bc109a9 955 for (Int_t part = 0; part < fNt; part++) {
956 if (!fTrackList[part]) continue;
957 Float_t phi = fPhiT[part];
958 Float_t eta = fEtaT[part];
959 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
960 (phij-phi)*(phij-phi));
961 if (dr < fConeRadius) {
962 fTrackList[part] = nj+2;
963 nT++;
964 } // < ConeRadius ?
965 } // particle loop
d594b22e 966
ff114de3 967// Same for background event if available
968 Int_t nTB = 0;
969 if (fBackground) {
970 for (Int_t part = 0; part < fNtB; part++) {
971 Float_t phi = fPhiB[part];
972 Float_t eta = fEtaB[part];
973 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
974 (phij-phi)*(phij-phi));
610f27c9 975 fTrackListB[part] = 1;
976
ff114de3 977 if (dr < fConeRadius) {
610f27c9 978 fTrackListB[part] = nj+2;
ff114de3 979 nTB++;
980 } // < ConeRadius ?
981 } // particle loop
982 } // Background available ?
983
984 Int_t nT0 = nT + nTB;
985
986 if (nT0 > 50) nT0 = 50;
d594b22e 987
610f27c9 988 Float_t* ptT = new Float_t[nT0];
989 Float_t* etaT = new Float_t[nT0];
990 Float_t* phiT = new Float_t[nT0];
d594b22e 991 Int_t iT = 0;
992 Int_t j;
993
3bc109a9 994 for (Int_t part = 0; part < fNt; part++) {
995 if (fTrackList[part] == nj+2) {
d594b22e 996 Int_t index = 0;
997 for (j=iT-1; j>=0; j--) {
998 if (fPtT[part] > ptT[j]) {
999 index = j+1;
1000 break;
1001 }
1002 }
1003 for (j=iT-1; j>=index; j--) {
1004 ptT [j+1] = ptT [j];
1005 etaT[j+1] = etaT[j];
1006 phiT[j+1] = phiT[j];
1007 }
1008 ptT [index] = fPtT [part];
1009 etaT[index] = fEtaT[part];
1010 phiT[index] = fPhiT[part];
ff114de3 1011 iT++;
3bc109a9 1012 } // particle associated
ff114de3 1013 if (iT > nT0) break;
3bc109a9 1014 } // particle loop
d594b22e 1015
ff114de3 1016 if (fBackground) {
1017 for (Int_t part = 0; part < fNtB; part++) {
610f27c9 1018 if (fTrackListB[part] == nj+2) {
ff114de3 1019 Int_t index = 0;
1020 for (j=iT-1; j>=0; j--) {
1021 if (fPtB[part] > ptT[j]) {
1022 index = j+1;
610f27c9 1023
ff114de3 1024 break;
1025 }
1026 }
1027 for (j=iT-1; j>=index; j--) {
1028 ptT [j+1] = ptT [j];
1029 etaT[j+1] = etaT[j];
1030 phiT[j+1] = phiT[j];
1031 }
1032 ptT [index] = fPtB [part];
1033 etaT[index] = fEtaB[part];
1034 phiT[index] = fPhiB[part];
1035 iT++;
1036 } // particle associated
1037 if (iT > nT0) break;
1038 } // particle loop
1039 } // Background available ?
1040
610f27c9 1041 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
3bc109a9 1042 delete[] ptT;
1043 delete[] etaT;
1044 delete[] phiT;
1045 } // jet loop loop
1046}
1047
26dadf3a 1048Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1049{
1050// Propagates phi angle to EMCAL radius
1051//
1052 Float_t dPhi = 0.;
1053// Get field
1054 Float_t b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
1055// Get EMCAL radius
1056 Float_t rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1057//
1058//
1059// bending radies
ff114de3 1060 Float_t rB = 333.56 * pt / b; // [cm]
26dadf3a 1061//
1062// check if particle is curling below EMCAL
1063 if (2.*rB < rEMCAL) {
1064 curls = kTRUE;
1065 return dPhi;
1066 }
1067//
1068// if not calculate delta phi
1069 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1070 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
f8f69f71 1071 dPhi = -TMath::Sign(dPhi, charge);
26dadf3a 1072//
1073 return dPhi;
1074
1075}
1076
1077
471f69dc 1078void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1079{
1080// dummy for hbook calls
1081 ;
1082}
1083
1084
d594b22e 1085
1086
1087
1088
1089
1090