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