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