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