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