8c8a1f39fc3c7aafe93caf854a64eaae870e2a02
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinder.cxx
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
16 /* $Id$ */
17
18 //*-- Authors: Andreas Morsch   (CERN)
19 //*            J.L. Klay        (LBL)
20 //*            Aleksei Pavlinov (WSU) 
21
22 #include <stdio.h>
23 // From root ...
24 #include <TArrayF.h>
25 #include <TAxis.h>
26 #include <TBranchElement.h>
27 #include <TCanvas.h>
28 #include <TClonesArray.h>
29 #include <TFile.h>
30 #include <TH1.h>
31 #include <TH2.h>
32 #include <TList.h>
33 #include <TPDGCode.h>
34 #include <TPad.h>
35 #include <TParticle.h>
36 #include <TParticlePDG.h>
37 #include <TPaveText.h>
38 #include <TPythia6Calls.h>
39 #include <TROOT.h>
40 #include <TStyle.h>
41 #include <TTree.h>
42 #include <TBrowser.h>
43
44 // From AliRoot ...
45 #include "AliEMCAL.h"
46 #include "AliEMCALDigit.h"
47 #include "AliEMCALDigitizer.h"
48 #include "AliEMCALFast.h"
49 #include "AliEMCALGeometry.h"
50 #include "AliEMCALHadronCorrection.h"
51 #include "AliEMCALHit.h"
52 #include "AliEMCALJetFinder.h"
53 #include "AliEMCALJetMicroDst.h"
54 #include "AliHeader.h"
55 #include "AliMagF.h"
56 #include "AliMagFCM.h"
57 #include "AliRun.h"
58 #include "AliGenerator.h"
59 #include "AliEMCALGetter.h"
60 // Interface to FORTRAN
61 #include "Ecommon.h"
62 #include "AliMC.h"
63
64
65 ClassImp(AliEMCALJetFinder)
66
67 //____________________________________________________________________________
68 AliEMCALJetFinder::AliEMCALJetFinder()
69 {
70 // Default constructor
71     fJets             = 0;
72     fNjets            = 0;
73     fLego             = 0;
74     fLegoB            = 0;
75
76     fTrackList        = 0;
77     fPtT              = 0;
78     fEtaT             = 0;
79     fPhiT             = 0;
80     fPdgT             = 0;
81     
82     fTrackListB       = 0;
83     fPtB              = 0;
84     fEtaB             = 0;
85     fPhiB             = 0;
86     fPdgB             = 0;
87
88     fHCorrection      = 0;
89     fHadronCorrector  = 0;
90
91     fWrite            = 0;
92     fOutFileName      = 0;
93     fOutFile          = 0;
94     fInFile           = 0;
95     fEvent            = 0;
96
97     fRandomBg         = 0;
98     
99     SetParametersForBgSubtraction();
100 }
101
102 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
103     : TTask(name, title)
104 {
105 // Constructor 
106 // Title is used in method GetFileNameForParameters();
107 //
108     fJets  = new TClonesArray("AliEMCALJet",10000);
109     fNjets = 0;
110     for (Int_t i = 0; i < 30000; i++)
111     {
112         fEtCell[i]  = 0.;
113         fEtaCell[i] = 0.;
114         fPhiCell[i] = 0.;
115     }
116     fLego       = 0;
117     fLegoB      = 0;
118
119     fTrackList  = 0;
120     fPtT        = 0;
121     fEtaT       = 0;
122     fPhiT       = 0;
123     fPdgT       = 0;
124
125     fTrackListB       = 0;
126     fPtB        = 0;
127     fEtaB       = 0;
128     fPhiB       = 0;
129     fPdgB       = 0;
130
131     fHCorrection      = 0;
132     fHadronCorrector  = 0;
133     fBackground       = 0;
134     fWrite            = 0;
135     fOutFileName      = 0;
136     fOutFile          = 0;
137     fInFile           = 0;
138     fEvent            = 0;
139 //
140     SetPtCut();
141     SetMomentumSmearing();
142     SetEfficiencySim();
143     SetDebug();
144     SetHadronCorrection();
145     SetIncludeK0andN();
146
147     fRandomBg         = 0;
148     SetParametersForBgSubtraction();
149 }
150
151 void AliEMCALJetFinder::SetParametersForBgSubtraction
152 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
153 {
154 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
155 // at WSU Linux cluster - 11-feb-2002
156 // These parameters must be tuned more carefull !!!
157   SetMode(mode);
158   SetMinMove(minMove);
159   SetMaxMove(maxMove);
160   SetPrecBg(precBg);
161 }
162
163 //____________________________________________________________________________
164 AliEMCALJetFinder::~AliEMCALJetFinder()
165 {
166 // Destructor
167 //
168     if (fJets){
169         fJets->Delete();
170         delete fJets;
171     }
172     delete fLego;            
173     delete fLegoB;
174     delete fhLegoTracks;  
175     delete fhLegoEMCAL;   
176     delete fhLegoHadrCorr;
177     delete fhEff;         
178     delete fhCellEt;      
179     delete fhCellEMCALEt; 
180     delete fhTrackPt;     
181     delete fhTrackPtBcut; 
182     delete fhChPartMultInTpc;
183
184     delete[] fTrackList;
185     delete[] fPtT;      
186     delete[] fEtaT;     
187     delete[] fPhiT;     
188     delete[] fPdgT;     
189  
190     delete[] fTrackListB;
191     delete[] fPtB;       
192     delete[] fEtaB;      
193     delete[] fPhiB;      
194     delete[] fPdgB;      
195 }
196
197 #ifndef WIN32
198 # define jet_finder_ua1 jet_finder_ua1_
199 # define hf1 hf1_
200 # define type_of_call
201
202 #else
203 # define jet_finder_ua1 JET_FINDER_UA1
204 # define hf1 HF1
205 # define type_of_call _stdcall
206 #endif
207
208 extern "C" void type_of_call 
209 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
210                Float_t  etc[30000],  Float_t etac[30000],
211                Float_t  phic[30000], 
212                Float_t& min_move, Float_t& max_move, Int_t& mode, 
213                Float_t& prec_bg,  Int_t& ierror);
214
215 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
216
217
218 void AliEMCALJetFinder::Init()
219 {
220 //
221 // Geometry and I/O initialization
222 //
223 //
224 //
225 //  Get geometry parameters from EMCAL
226 //
227 //
228 //  Geometry 
229   //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
230   //  AliEMCALGeometry* geom = 
231   //  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
232   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
233   AliEMCALGeometry* geom = gime->EMCALGeometry() ;  
234
235 //    SetSamplingFraction(geom->GetSampling());
236
237     fNbinEta = geom->GetNZ();
238     fNbinPhi = geom->GetNPhi();
239     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
240     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
241     fEtaMin  = geom->GetArm1EtaMin();
242     fEtaMax  = geom->GetArm1EtaMax();
243     fDphi    = (fPhiMax-fPhiMin)/fNbinPhi;
244     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
245     fNtot    = fNbinPhi*fNbinEta;
246     fWeightingMethod = kFALSE;
247
248 //
249     SetCellSize(fDeta, fDphi);
250 //
251 //  I/O
252     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
253 //
254 //
255 }
256
257 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], 
258                              Float_t etac[30000], Float_t phic[30000],
259                              Float_t min_move, Float_t max_move, Int_t mode, 
260                              Float_t prec_bg,  Int_t   ierror)
261 {
262 // Wrapper for fortran coded jet finder
263 // Full argument list
264     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic, 
265                    min_move, max_move, mode, prec_bg, ierror);
266     // Write result to output
267     if(fWrite) WriteJets();
268     fEvent++;
269 }
270
271 void AliEMCALJetFinder::Find()
272 {
273 // Wrapper for fortran coded jet finder using member data for 
274 // argument list
275
276     Float_t min_move = fMinMove;
277     Float_t max_move = fMaxMove;
278     Int_t   mode     = fMode;
279     Float_t prec_bg  = fPrecBg;
280     Int_t   ierror;
281     ResetJets(); // 4-feb-2002 by PAI
282
283     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, 
284                    min_move, max_move, mode, prec_bg, ierror);
285     fError = ierror;
286     // Write result to output
287     Int_t njet = Njets();
288     
289     for (Int_t nj=0; nj<njet; nj++)
290     {
291         if (fWeightingMethod)
292         {
293           fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
294                                       JetPhiW(nj),
295                                       JetEtaW(nj));
296
297         }else{
298                 
299         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
300                                     JetPhiW(nj),
301                                     JetEtaW(nj));
302         }
303         fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod); 
304         fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) ); 
305         fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj)  )); 
306         fJetT[nj]->SetHCEnergy(HCConeEnergy( JetEtaW(nj),JetPhiW(nj)  )); 
307         
308     }
309
310     FindTracksInJetCone();
311     if(fWrite) WriteJets();
312     fEvent++;
313 }
314
315
316 Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
317 {
318 Float_t newenergy = 0.0;
319 Float_t bineta,binphi;
320 TAxis *x = fhLegoEMCAL->GetXaxis();
321 TAxis *y = fhLegoEMCAL->GetYaxis();
322 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
323         {
324    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
325       {
326          bineta = x->GetBinCenter(i);
327          binphi = y->GetBinCenter(j);
328        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
329       {
330            newenergy += fhLegoEMCAL->GetBinContent(i,j);
331       }
332     }
333 }
334 return newenergy;
335 }       
336
337 Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
338 {
339 Float_t newenergy = 0.0;
340 Float_t bineta,binphi;
341 TAxis *x = fhLegoTracks->GetXaxis();
342 TAxis *y = fhLegoTracks->GetYaxis();
343 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
344    {
345    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
346    {
347     bineta = x->GetBinCenter(i);
348     binphi = y->GetBinCenter(j);
349     if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
350       {
351           newenergy += fhLegoTracks->GetBinContent(i,j);
352        }
353     }
354 }
355 return newenergy;
356 }
357
358 Float_t AliEMCALJetFinder::HCConeEnergy(Float_t eta, Float_t phi)
359 {
360 //Float_t newenergy = 0.0;
361 //Float_t bineta,binphi;
362 //TAxis *x = fhLegoTracks->GetXaxis();
363 //TAxis *y = fhLegoTracks->GetYaxis();
364 //for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
365 //   {
366 //      for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
367 //      {
368 //        bineta = x->GetBinCenter(i);
369 //        binphi = y->GetBinCenter(j);
370 //        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
371 //        {
372 //            newenergy += fhLegoTracks->GetBinContent(i,j);
373 //         }
374 //    }
375 //}
376 //return newenergy;
377         
378 return eta*phi*0.0;     
379         
380 }
381
382
383
384 Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
385 {
386
387
388 Float_t newenergy = 0.0;
389 Float_t bineta,binphi;
390 TAxis *x = fhLegoEMCAL->GetXaxis();
391 TAxis *y = fhLegoEMCAL->GetYaxis();
392
393
394 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
395 {
396    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
397    {
398       bineta = x->GetBinCenter(i);
399       binphi = y->GetBinCenter(j);
400       if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
401       {
402           newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
403       }
404     }
405 }
406
407 return newenergy;
408
409 }
410
411         
412 Int_t AliEMCALJetFinder::Njets()
413 {
414 // Get number of reconstructed jets
415     return EMCALJETS.njet;
416 }
417
418 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
419 {
420 // Get reconstructed jet energy
421     return EMCALJETS.etj[i];
422 }
423
424 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
425 {
426 // Get reconstructed jet phi from leading particle
427     return EMCALJETS.phij[0][i];
428 }
429
430 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
431 {
432 // Get reconstructed jet phi from weighting
433     return EMCALJETS.phij[1][i];
434 }
435
436 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
437 {
438 // Get reconstructed jet eta from leading particles
439     return EMCALJETS.etaj[0][i];
440 }
441
442
443 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)  
444 {
445 // Get reconstructed jet eta from weighting
446     return EMCALJETS.etaj[1][i];
447 }
448
449 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
450 {
451 // Set grid cell size
452     EMCALCELLGEO.etaCellSize = eta;
453     EMCALCELLGEO.phiCellSize = phi;    
454 }
455
456 void AliEMCALJetFinder::SetConeRadius(Float_t par)
457 {
458 // Set jet cone radius
459     EMCALJETPARAM.coneRad = par;
460     fConeRadius = par;
461 }
462
463 void AliEMCALJetFinder::SetEtSeed(Float_t par)
464 {
465 // Set et cut for seeds
466     EMCALJETPARAM.etSeed = par;
467     fEtSeed = par;
468 }
469
470 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
471 {
472 // Set minimum jet et
473     EMCALJETPARAM.ejMin = par;
474     fMinJetEt = par;
475 }
476
477 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
478 {
479 // Set et cut per cell
480     EMCALJETPARAM.etMin = par;
481     fMinCellEt = par;
482 }
483
484 void AliEMCALJetFinder::SetPtCut(Float_t par)
485 {
486 // Set pT cut on charged tracks
487     fPtCut = par;
488 }
489
490
491 void AliEMCALJetFinder::Test()
492 {
493 //
494 // Test the finder call
495 //
496     const Int_t nmax = 30000;
497     Int_t ncell      = 10;
498     Int_t ncell_tot  = 100;
499
500     Float_t etc[nmax];
501     Float_t etac[nmax];
502     Float_t phic[nmax];
503     Float_t min_move = 0;
504     Float_t max_move = 0;
505     Int_t   mode = 0;
506     Float_t prec_bg = 0;
507     Int_t   ierror = 0;
508
509     
510     Find(ncell, ncell_tot, etc, etac, phic, 
511          min_move, max_move, mode, prec_bg, ierror);
512
513 }
514
515 //
516 //  I/O
517 //      
518
519 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
520 {
521     //
522     // Add a jet 
523     //
524     TClonesArray &lrawcl = *fJets;
525     new(lrawcl[fNjets++]) AliEMCALJet(jet);
526 }
527
528 void AliEMCALJetFinder::ResetJets()
529 {
530     //
531     // Reset Jet List 
532     //
533     fJets->Clear();
534     fNjets = 0;
535 }
536
537 void AliEMCALJetFinder::WriteJets()
538 {
539 //
540 // Add all jets to the list
541 //
542     const Int_t kBufferSize = 4000;
543     const char* file = 0;
544
545     Int_t njet = Njets();
546
547     for (Int_t nj = 0; nj < njet; nj++)
548     {
549         AddJet(*fJetT[nj]);
550         delete fJetT[nj];
551     }
552
553 // I/O
554     AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
555     if (!fOutFileName) {
556 //
557 // output written to input file
558 //
559       AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
560       TTree* pK = gAlice->TreeK();
561       file = (pK->GetCurrentFile())->GetName();
562       TBranch * jetBranch ;  
563         if (fDebug > 1)
564             printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
565         //if (fJets && gAlice->TreeR()) {
566         if (fJets && gime->TreeR()) {
567           // pEMCAL->MakeBranchInTree(gAlice->TreeR(), 
568           jetBranch = gime->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ; 
569           //pEMCAL->MakeBranchInTree(gime->TreeR(), 
570           //                         "EMCALJets", 
571           //                         &fJets, 
572           //                         kBufferSize, 
573           //                         file);
574           
575           //Int_t nev = gAlice->GetHeader()->GetEvent();
576           //gAlice->TreeR()->Fill();
577           jetBranch->Fill();
578           //char hname[30];
579           //sprintf(hname,"TreeR%d", nev);
580           //gAlice->TreeR()->Write(hname);
581           //gAlice->TreeR()->Reset();
582           gime->WriteRecPoints("OVERWRITE");
583         }
584     } else {
585 //
586 // Output written to user specified output file
587 //
588       //TTree* pK = gAlice->TreeK();
589       TTree* pK = gAlice->TreeK();
590       fInFile  = pK->GetCurrentFile();
591
592       fOutFile->cd();
593       char hname[30];
594       sprintf(hname,"TreeR%d", fEvent);
595       TTree* treeJ = new TTree(hname, "EMCALJets");
596       treeJ->Branch("EMCALJets", &fJets, kBufferSize);
597       treeJ->Fill();
598       treeJ->Write(hname);
599       fInFile->cd();
600     }
601     ResetJets();        
602 }
603
604 void AliEMCALJetFinder::BookLego()
605 {
606 //
607 //  Book histo for discretization
608 //
609
610 //
611 //  Don't add histos to the current directory
612     if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
613  
614     //    TH2::AddDirectory(0); // hists wil be put to the list from gROOT
615     //    TH1::AddDirectory(0);
616     gROOT->cd();
617 //    
618 //  Signal map
619     fLego = new TH2F("legoH","eta-phi",
620                            fNbinEta, fEtaMin, fEtaMax, 
621                            fNbinPhi, fPhiMin, fPhiMax);
622 //
623 //  Background map
624     fLegoB = new TH2F("legoB","eta-phi for BG event",
625                            fNbinEta, fEtaMin, fEtaMax, 
626                            fNbinPhi, fPhiMin, fPhiMax);
627
628 //  Tracks map
629     fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
630     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
631 //  EMCAL map
632     fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
633     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
634 //  Hadron correction map
635     fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
636     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
637 //  Hists. for tuning jet finder 
638     fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
639
640     TArrayF eTmp(1101);
641     eTmp[0] = 0.0;
642     for(Int_t i=1; i<=1000; i++)      eTmp[i] = 0.1*i; // step 100 mev
643     for(Int_t i=1001; i<=1100; i++)   eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
644
645     fhCellEt  = new TH1F("hCellEt","Cell E_{T} from fLego", 
646     eTmp.GetSize()-1, eTmp.GetArray()); 
647     fhCellEMCALEt  = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself", 
648     eTmp.GetSize()-1, eTmp.GetArray()); 
649     fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
650     eTmp.GetSize()-1, eTmp.GetArray()); 
651     fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
652     eTmp.GetSize()-1, eTmp.GetArray()); 
653
654     fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
655     "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
656
657     fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
658     TAxis *xax = fhSinTheta->GetXaxis();
659     for(Int_t i=1; i<=fNbinEta; i++) {
660       Double_t eta = xax->GetBinCenter(i);
661       fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
662     }
663     
664     //! first canvas for drawing
665     fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
666 }
667
668 void AliEMCALJetFinder::DumpLego()
669 {
670 //
671 // Dump lego histo into array
672 //
673     fNcell = 0;
674     TAxis* Xaxis = fLego->GetXaxis();
675     TAxis* Yaxis = fLego->GetYaxis();
676     //    fhCellEt->Clear();
677     Float_t e, eH;
678     for (Int_t i = 1; i <= fNbinEta; i++) {
679         for (Int_t j = 1; j <= fNbinPhi; j++) {
680             
681             e = fLego->GetBinContent(i,j);
682             if (fRandomBg) {
683                 if (gRandom->Rndm() < 0.5) {
684                     Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.));
685                     e += ebg;
686                 }
687             }
688             if (e > 0.0) e -= fMinCellEt;
689             if (e < 0.0) e = 0.;
690             Float_t eta  = Xaxis->GetBinCenter(i);
691             Float_t phi  = Yaxis->GetBinCenter(j);          
692             fEtCell[fNcell]  = e;
693             fEtaCell[fNcell] = eta;
694             fPhiCell[fNcell] = phi;
695             fNcell++;
696             fhCellEt->Fill(e);
697             if(fhLegoEMCAL) {
698               eH = fhLegoEMCAL->GetBinContent(i,j);
699               if(eH > 0.0) fhCellEMCALEt->Fill(eH);
700             }
701         } // phi
702     } // eta
703 //  today
704 //    fNcell--;
705 }
706
707 void AliEMCALJetFinder::ResetMap()
708 {
709 //
710 // Reset eta-phi array
711
712     for (Int_t i=0; i<30000; i++)
713     {
714         fEtCell[i]  = 0.;
715         fEtaCell[i] = 0.;
716         fPhiCell[i] = 0.;
717     }
718 }
719
720
721 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
722 {
723 // Which generator
724 //
725     const char* name = gAlice->Generator()->GetName();
726     enum {kPythia, kHijing, kHijingPara};
727     Int_t genType = 0;
728     
729     if (!strcmp(name, "Hijing")){
730         genType = kHijing;
731     } else if (!strcmp(name, "Pythia")) {
732         genType = kPythia;
733     } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) {
734         genType = kHijingPara;
735     }
736     if (fDebug>=2) 
737         printf("Fill tracks generated by %s %d\n", name, genType);
738         
739             
740 //
741 // Fill Cells with track information
742 //
743     if (fDebug >= 2)
744     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
745     fNChTpc = 0;
746
747     ResetMap();
748     
749     if (!fLego) BookLego();
750 // Reset
751     if (flag == 0) fLego->Reset();
752 //
753 // Access particle information    
754     Int_t npart = (gAlice->GetHeader())->GetNprimary();
755     Int_t ntr   = (gAlice->GetHeader())->GetNtrack();
756     printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n", 
757     npart, ntr, fLego->Integral());
758  
759 // Create track list
760 //
761 // 0: not selected
762 // 1: selected for jet finding
763 // 2: inside jet
764 // ....
765     if (fTrackList) delete[] fTrackList;
766     if (fPtT)       delete[] fPtT;
767     if (fEtaT)      delete[] fEtaT;
768     if (fPhiT)      delete[] fPhiT;
769     if (fPdgT)      delete[] fPdgT;
770     
771     fTrackList = new Int_t  [npart];
772     fPtT       = new Float_t[npart];
773     fEtaT      = new Float_t[npart];
774     fPhiT      = new Float_t[npart];
775     fPdgT      = new Int_t[npart];
776
777     fNt   = npart;
778     fNtS  = 0;
779     Float_t chTmp=0.0; // charge of current particle 
780     //    Int_t idGeant;
781
782     // this is for Pythia ??
783     for (Int_t part = 0; part < npart; part++) {
784         TParticle *MPart = gAlice->GetMCApp()->Particle(part);
785         Int_t mpart   = MPart->GetPdgCode();
786         Int_t child1  = MPart->GetFirstDaughter();
787         Float_t pT    = MPart->Pt();
788         Float_t p     = MPart->P();
789         Float_t phi   = MPart->Phi();
790         Float_t eta   = -100.;
791         if(pT > 0.001) eta   = MPart->Eta();
792         Float_t theta = MPart->Theta(); 
793         if  (fDebug>=2 && MPart->GetStatusCode()==1) { 
794            printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n", 
795            part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
796         }
797         
798         if (fDebug >= 2 && genType == kPythia) {
799             if (part == 6 || part == 7)
800             {
801                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n", 
802                        part-5, pT, eta, phi);
803             }
804         }
805         
806         fTrackList[part] = 0; 
807         fPtT[part]       = pT; // must be change after correction for resolution !!!
808         fEtaT[part]      = eta;
809         fPhiT[part]      = phi;
810         fPdgT[part]      = mpart;
811         
812 //      final state particles only
813
814         if (genType == kPythia) {
815             if (MPart->GetStatusCode() != 1) continue;
816         } else if (genType == kHijing) {
817             if (child1 >= 0 && child1 < npart) continue;
818         }
819         
820         
821         TParticlePDG* pdgP = 0;
822 // charged or neutral 
823         pdgP  = MPart->GetPDG();
824         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!  
825
826         if (ich == 0) {
827             if (chTmp == 0) {
828                 if (!fK0N) { 
829                     continue;
830                 } else {
831                     if (mpart != kNeutron    &&
832                         mpart != kNeutronBar &&
833                         mpart != kK0Long) continue;
834                 }
835             }
836         } else if (ich == 2) {
837           if (mpart == kNeutron    ||
838               mpart == kNeutronBar ||
839               mpart == kK0Long) continue;
840         }
841
842         if (TMath::Abs(eta)<=0.9) fNChTpc++;
843 // final state only
844         if (child1 >= 0 && child1 < npart) continue;
845 // acceptance cut
846         if (eta > fEtaMax || eta < fEtaMin)    continue;
847         if (phi > fPhiMax || phi < fPhiMin )   continue;
848 //
849         if (fDebug >= 3) 
850         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
851         part, mpart, child1, eta, phi, pT, chTmp);
852 //
853 //
854 // Momentum smearing goes here ...
855 //
856         fhTrackPt->Fill(pT);
857         Float_t pw;
858         if (fSmear && TMath::Abs(chTmp)) {
859             pw = AliEMCALFast::SmearMomentum(1,p);
860             // p changed - take into account when calculate pT,
861             // pz and so on ..  ?
862             pT = (pw/p) * pT;
863             if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
864             p  = pw;
865         }
866 //
867 // Tracking Efficiency and TPC acceptance goes here ...
868         Float_t eff;
869         if (fEffic && TMath::Abs(chTmp)) {
870             eff = 0.9; //  AliEMCALFast::Efficiency(2,p);
871             if(fhEff) fhEff->Fill(p, eff);
872             if (AliEMCALFast::RandomReject(eff)) {
873                 if(fDebug >= 5) printf(" reject due to unefficiency ");
874                 continue;
875             }
876         }
877 //
878 // Correction of Hadronic Energy goes here ...
879 //
880 //
881 // phi propagation for hadronic correction
882
883         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
884         Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
885         if(TMath::Abs(chTmp)) {
886             // hadr. correction only for charge particle 
887             dphi  = PropagatePhi(pT, chTmp, curls);
888             phiHC = phi + dphi;
889             if (fDebug >= 6) {
890                 printf("\n Delta phi %f pT %f ", dphi, pT);
891                 if (curls) printf("\n !! Track is curling");
892             }
893             if(!curls) fhTrackPtBcut->Fill(pT);
894             
895             if (fHCorrection && !curls) {
896                 if (!fHadronCorrector)
897                     Fatal("AliEMCALJetFinder",
898                           "Hadronic energy correction required but not defined !");
899                 
900                 dpH    = fHadronCorrector->GetEnergy(p, eta, 7);
901                 eTdpH  = dpH*TMath::Sin(theta);
902                 
903                 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n", 
904                                         phi, phiHC, -eTdpH); // correction is negative
905                 Int_t xbin,ybin;
906                 xbin = fLego->GetXaxis()->FindBin(eta);
907                 ybin = fLego->GetYaxis()->FindBin(phiHC);
908                 cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
909                 fLego->Fill(eta, phi, -fSamplingF*eTdpH );
910                 cout <<"Hadron Correction affected bin - contents after  correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
911                 fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
912             }
913         }
914 //
915 //  More to do ? Just think about it !
916 //
917         if (phi > fPhiMax || phi < fPhiMin )   continue;
918
919         if(TMath::Abs(chTmp) ) { // charge particle
920             if (pT > fPtCut && !curls) {
921                 if (fDebug >= 8) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
922                                         eta , phi, pT, fNtS); 
923                 fLego->Fill(eta, phi, pT);
924                 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
925                 fTrackList[part] = 1;
926                 fNtS++;
927             }
928         } else if(ich > 0 || fK0N) {
929             // case of n, nbar and K0L
930             if (fDebug >= 9) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
931                                     eta , phi, pT, fNtS); 
932             fLego->Fill(eta, phi, pT);
933             fTrackList[part] = 1;
934             fNtS++;
935         }
936     } // primary loop    
937     Float_t etsum = 0.;
938     for(Int_t i=0; i<fLego->GetSize(); i++) {
939         Float_t etc =  (*fLego)[i];
940         if (etc > fMinCellEt) etsum += etc;
941     }
942     
943     printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
944     printf("                Track selected(fNtS) %i \n", fNtS);
945
946     DumpLego();
947 }
948
949 void AliEMCALJetFinder::FillFromHits(Int_t flag)
950 {
951 //
952 // Fill Cells with hit information
953 //
954 //
955     if (fDebug >= 2)
956     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
957
958     ResetMap();
959     
960     if (!fLego) BookLego();
961 //  Reset eta-phi maps if needed
962     if (flag == 0) { // default behavior
963       fLego->Reset();
964       fhLegoTracks->Reset();
965       fhLegoEMCAL->Reset();
966       fhLegoHadrCorr->Reset();
967     }
968 //  Initialize from background event if available
969 //
970 // Access hit information    
971     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
972     AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
973     TTree *treeH = gime->TreeH();
974     Int_t ntracks = (Int_t) treeH->GetEntries();
975 //
976 //   Loop over tracks
977 //
978     Int_t nbytes = 0;
979     //    Double_t etH = 0.0;
980
981     for (Int_t track=0; track<ntracks;track++) {
982         gAlice->ResetHits();
983         nbytes += treeH->GetEvent(track);
984 //
985 //  Loop over hits
986 //
987         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); 
988             mHit;
989             mHit=(AliEMCALHit*) pEMCAL->NextHit()) 
990         {
991             Float_t x      =    mHit->X();         // x-pos of hit
992             Float_t y      =    mHit->Y();         // y-pos
993             Float_t z      =    mHit->Z();         // z-pos
994             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
995
996             Float_t r      =    TMath::Sqrt(x*x+y*y);
997             Float_t theta  =    TMath::ATan2(r,z);
998             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
999             Float_t phi    =    TMath::ATan2(y,x);
1000
1001             if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
1002             
1003             //            etH = fSamplingF*eloss*TMath::Sin(theta);
1004             fLego->Fill(eta, phi, eloss);
1005         } // Hit Loop
1006     } // Track Loop
1007
1008     // Transition from deposit energy to eT (eT = de*SF*sin(theta))
1009     Double_t etsum = 0;    
1010     for(Int_t i=1; i<=fLego->GetNbinsX(); i++){   // eta
1011        Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
1012        for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
1013           eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
1014           fLego->SetBinContent(i,j,eT);
1015     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)    
1016           fhLegoEMCAL->SetBinContent(i,j,eT);
1017           if (eT > fMinCellEt) etsum += eT;
1018        }
1019     }
1020
1021     //    for(Int_t i=0; i<fLego->GetSize(); i++) {
1022     //  (*fhLegoEMCAL)[i] = (*fLego)[i];
1023     //  Float_t etc =  (*fLego)[i];
1024     //  if (etc > fMinCellEt) etsum += etc;
1025     //    }
1026     
1027     printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
1028     //    DumpLego(); ??
1029 }
1030
1031 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
1032 {
1033 //
1034 // Fill Cells with digit information
1035 //
1036 //
1037     ResetMap();
1038     
1039     if (!fLego) BookLego();
1040     if (flag == 0) fLego->Reset();
1041     Int_t nbytes=0;
1042     
1043
1044 //
1045 //  Connect digits    
1046 //
1047     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
1048     TTree *treeD = gAlice->TreeD();
1049     TBranchElement* branchDg = (TBranchElement*)
1050         treeD->GetBranch("EMCAL");
1051
1052     if (!branchDg) Fatal("AliEMCALJetFinder", 
1053                          "Reading digits requested but no digits in file !");
1054     
1055     branchDg->SetAddress(&digs);
1056     Int_t nent = (Int_t) branchDg->GetEntries();
1057 //
1058 //  Connect digitizer
1059 //
1060     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
1061     TBranchElement* branchDr = (TBranchElement*) 
1062         treeD->GetBranch("AliEMCALDigitizer");
1063     branchDr->SetAddress(&digr);
1064 //
1065 //
1066     nbytes += branchDg->GetEntry(0);
1067     nbytes += branchDr->GetEntry(0);
1068 //
1069 //  Get digitizer parameters
1070     Float_t preADCped = digr->GetPREpedestal();
1071     Float_t preADCcha = digr->GetPREchannel();
1072     Float_t ecADCped  = digr->GetECApedestal();
1073     Float_t ecADCcha  = digr->GetECAchannel();
1074     Float_t hcADCped  = digr->GetHCApedestal();
1075     Float_t hcADCcha  = digr->GetHCAchannel();
1076
1077     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1078     AliEMCALGeometry* geom = 
1079         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
1080     
1081     if (fDebug) {
1082         Int_t   ndig = digs->GetEntries();
1083         Info("FillFromDigits","Number of Digits: %d %d\n Parameters: PRE : %f %f EC: %f %f HC: %f %f\n Geometry: %d %d", 
1084              ndig, nent, preADCped, preADCcha, ecADCped, ecADCcha, hcADCped, hcADCcha, geom->GetNEta(), geom->GetNPhi());
1085     }
1086     
1087 //
1088 //  Loop over digits
1089     AliEMCALDigit* sdg;
1090     TIter next(digs);
1091     while ((sdg = (AliEMCALDigit*)(next())))
1092     {
1093         Double_t pedestal = 0.;
1094         Double_t channel  = 0.;
1095         if (geom->IsInPRE(sdg->GetId())) {
1096           pedestal = preADCped;
1097           channel  = preADCcha; 
1098         } 
1099         else if (geom->IsInECA(sdg->GetId()))  {
1100           pedestal = ecADCped;
1101           channel  = ecADCcha; 
1102         }
1103         else if (geom->IsInHCA(sdg->GetId()))  {
1104           pedestal = hcADCped;
1105           channel  = hcADCcha; 
1106         }
1107         else 
1108           Fatal("FillFromDigits", "unexpected digit is number!") ; 
1109         
1110         Float_t eta = sdg->GetEta();
1111         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1112         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1113         
1114         if (fDebug > 1) 
1115           Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1116                            eta, phi, amp, sdg->GetAmp());
1117         
1118         fLego->Fill(eta, phi, fSamplingF*amp);
1119     } // digit loop
1120 // 
1121 //  Dump lego hist
1122     DumpLego();
1123 }
1124
1125
1126 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1127 {
1128 //
1129 // Fill Cells with hit information
1130 //
1131 //
1132     ResetMap();
1133     
1134     if (!fLego) BookLego();
1135 // Reset
1136     if (flag == 0) fLego->Reset();
1137
1138 // Flag charged tracks passing through TPC which 
1139 // are associated to EMCAL Hits
1140     BuildTrackFlagTable();
1141
1142 //
1143 // Access particle information    
1144     TTree *treeK = gAlice->TreeK();
1145     Int_t ntracks = (Int_t) treeK->GetEntries();
1146
1147     if (fPtT)       delete[] fPtT;   
1148     if (fEtaT)      delete[] fEtaT;    
1149     if (fPhiT)      delete[] fPhiT;   
1150     if (fPdgT)      delete[] fPdgT;   
1151    
1152     fPtT       = new Float_t[ntracks];
1153     fEtaT      = new Float_t[ntracks];
1154     fPhiT      = new Float_t[ntracks];
1155     fPdgT      = new Int_t[ntracks];
1156
1157     fNt   = ntracks;
1158     fNtS  = 0;
1159     
1160     for (Int_t track = 0; track < ntracks; track++) {
1161           TParticle *MPart = gAlice->GetMCApp()->Particle(track);
1162           Float_t pT    = MPart->Pt();
1163           Float_t phi   = MPart->Phi();
1164           Float_t eta   = MPart->Eta();
1165
1166         if(fTrackList[track]) {
1167           fPtT[track]       = pT;
1168           fEtaT[track]      = eta;
1169           fPhiT[track]      = phi;
1170           fPdgT[track]      = MPart->GetPdgCode();
1171           
1172           if (track < 2) continue;      //Colliding particles?
1173           if (pT == 0 || pT < fPtCut) continue;
1174           fNtS++;
1175           fLego->Fill(eta, phi, pT);
1176         }
1177       } // track loop
1178     DumpLego();
1179 }
1180
1181 void AliEMCALJetFinder::FillFromParticles()
1182 {
1183 // 26-feb-2002 PAI - for checking all chain
1184 // Work on particles level; accept all particle (not neutrino )
1185     
1186     Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law 
1187     fNChTpc = 0;
1188
1189     ResetMap();
1190     if (!fLego) BookLego();
1191     fLego->Reset();
1192 //
1193 // Access particles information    
1194     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1195     if (fDebug >= 2 || npart<=0) {
1196        printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1197        if(npart<=0) return;
1198     }
1199     fNt   = npart;
1200     fNtS  = 0;
1201     RearrangeParticlesMemory(npart);
1202  
1203 //  Go through the particles
1204     Int_t mpart, child1, child2, geantPdg;
1205     Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1206     TParticle *MPart=0;
1207     for (Int_t part = 0; part < npart; part++) {
1208
1209         fTrackList[part] = 0;
1210
1211         MPart   = gAlice->GetMCApp()->Particle(part);
1212         mpart   = MPart->GetPdgCode();
1213         child1  = MPart->GetFirstDaughter();
1214         child2  = MPart->GetLastDaughter();
1215         pT      = MPart->Pt();
1216         phi     = MPart->Phi();
1217         eta     = MPart->Eta();
1218
1219         px      = MPart->Px();
1220         py      = MPart->Py();
1221         pz      = MPart->Pz();
1222         e       = MPart->Energy();
1223
1224 // see pyedit in Pythia's text
1225         geantPdg = mpart;
1226         if (IsThisPartonsOrDiQuark(mpart)) continue;
1227         printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n", 
1228         part, mpart, geantPdg, px, py, pz, e, child1, child2);
1229         
1230 //  exclude partons (21 - gluon, 92 - string) 
1231         
1232
1233 // exclude neutrinous also ??
1234         if (fDebug >= 11 && pT>0.01) 
1235         printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1236         part, mpart, eta, phi, pT);
1237
1238         fPtT[part]       = pT;
1239         fEtaT[part]      = eta;
1240         fPhiT[part]      = phi;
1241         fPdgT[part]      = mpart;
1242         fNtS++;
1243         
1244 // final state only
1245         if (child1 >= 0 && child1 < npart) continue;
1246
1247         //        printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n", 
1248         //        part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1249         PX += px; 
1250         PY += py; 
1251         PZ += pz;
1252         E  += e; 
1253
1254         
1255         if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1256 // acceptance cut
1257         if (eta > fEtaMax || eta < fEtaMin)    continue;
1258         if (phi > fPhiMax || phi < fPhiMin )   continue;
1259 //
1260         if(fK0N==0 ) { // exclude neutral hadrons
1261           if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue; 
1262         }
1263         fTrackList[part] = 1;
1264         fLego->Fill(eta, phi, pT);
1265
1266     } // primary loop
1267     printf("\n                PX %8.2f  PY %8.2f  PZ %8.2f  E %8.2f \n", 
1268     PX, PY, PZ, E);
1269     DumpLego();
1270     if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1271 }
1272
1273 void AliEMCALJetFinder::FillFromPartons()
1274 {
1275 // function under construction - 13-feb-2002 PAI
1276     
1277     if (fDebug >= 2)
1278     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1279     // 
1280
1281     ResetMap();
1282     if (!fLego) BookLego();
1283     fLego->Reset();
1284 //
1285 // Access particle information    
1286     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1287     if (fDebug >= 2 || npart<=0)
1288     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1289     fNt   = 0; // for FindTracksInJetCone
1290     fNtS  = 0;
1291  
1292 //  Go through the partons
1293     Int_t statusCode=0;
1294     for (Int_t part = 8; part < npart; part++) {
1295         TParticle *MPart = gAlice->GetMCApp()->Particle(part);
1296         Int_t mpart   = MPart->GetPdgCode();
1297         //      Int_t child1  = MPart->GetFirstDaughter();
1298         Float_t pT    = MPart->Pt();
1299         //      Float_t p     = MPart->P();
1300         Float_t phi   = MPart->Phi();
1301         Float_t eta   = MPart->Eta();
1302         //      Float_t theta = MPart->Theta();
1303         statusCode    = MPart->GetStatusCode();
1304         
1305 // accept partons (21 - gluon, 92 - string) 
1306         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1307         if (fDebug > 1 && pT>0.01) 
1308         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
1309         part, mpart, statusCode, eta, phi, pT);
1310         //      if (fDebug >= 3) MPart->Print(); 
1311 // accept partons before fragmentation - p.57 in Pythia manual
1312 //        if(statusCode != 1) continue;
1313 // acceptance cut
1314         if (eta > fEtaMax || eta < fEtaMin)    continue;
1315         if (phi > fPhiMax || phi < fPhiMin )   continue;
1316 // final state only
1317 //      if (child1 >= 0 && child1 < npart) continue;
1318 //
1319 //
1320         fLego->Fill(eta, phi, pT);
1321
1322     } // primary loop
1323     DumpLego();
1324 }
1325
1326 void AliEMCALJetFinder::BuildTrackFlagTable() {
1327
1328 // Method to generate a lookup table for TreeK particles
1329 // which are linked to hits in the EMCAL
1330 //
1331 // --Author: J.L. Klay
1332 //
1333 // Access hit information    
1334     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1335     
1336     TTree *TK = gAlice->TreeK();                // Get the number of entries in the kine tree
1337     Int_t nKTrks = (Int_t) TK->GetEntries();    // (Number of particles created somewhere)
1338     
1339     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
1340     fTrackList = new Int_t[nKTrks];             //before generating a new one
1341     
1342     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
1343         fTrackList[i] = 0;
1344     }
1345     
1346     AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
1347     // TTree *treeH = gAlice->TreeH();
1348     TTree *treeH = gime->TreeH();
1349    Int_t ntracks = (Int_t) treeH->GetEntries();
1350 //
1351 //   Loop over tracks
1352 //
1353     Int_t nbytes = 0;
1354     
1355     for (Int_t track=0; track<ntracks;track++) {
1356         gAlice->ResetHits();
1357         nbytes += treeH->GetEvent(track);
1358         if (pEMCAL)  {
1359             
1360 //
1361 //  Loop over hits
1362 //
1363             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); 
1364                 mHit;
1365                 mHit=(AliEMCALHit*) pEMCAL->NextHit()) 
1366             {
1367                 Int_t   iTrk   =    mHit->Track();       // track number
1368                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1369                 
1370                 //Determine the origin point of this particle - it made a hit in the EMCAL
1371                 TParticle         *trkPart = gAlice->GetMCApp()->Particle(iTrk);
1372                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1373                 Int_t      trkCode = trkPart->GetPdgCode();
1374                 Double_t           trkChg;
1375                 if (trkCode < 10000) {          //Big Ions cause problems for 
1376                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1377                 } else {                                //likely to make it very far, set
1378                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1379                 }
1380                 Float_t            vxTrk   = trkPart->Vx();
1381                 Float_t            vyTrk   = trkPart->Vy();
1382                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk); 
1383                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1384                 
1385                 //Loop through the ancestry of the EMCAL entrance particles
1386                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1387                 while (ancestor != -1) {
1388                     TParticle     *ancPart = gAlice->GetMCApp()->Particle(ancestor);  //get particle info on ancestor
1389                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1390                     Int_t             ancCode = ancPart->GetPdgCode();
1391                     Double_t       ancChg;
1392                     if (ancCode < 10000) {
1393                         ancChg  = ancPdg->Charge();
1394                     } else {
1395                         ancChg  = 0.0;
1396                     }
1397                     Float_t           vxAnc   = ancPart->Vx();
1398                     Float_t           vyAnc   = ancPart->Vy();
1399                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1400                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1401                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1402                 }
1403                 
1404                 //Determine the origin point of the primary particle associated with the hit
1405                 TParticle     *primPart = gAlice->GetMCApp()->Particle(idprim);
1406                 TParticlePDG  *primPdg  = primPart->GetPDG();
1407                 Int_t      primCode = primPart->GetPdgCode();
1408                 Double_t           primChg;
1409                 if (primCode < 10000) {
1410                     primChg  = primPdg->Charge();
1411                 } else {
1412                     primChg  = 0.0;
1413                 }
1414                 Float_t    vxPrim = primPart->Vx();
1415                 Float_t    vyPrim = primPart->Vy();
1416                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1417                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1418             } // Hit Loop
1419         } //if (pEMCAL)
1420     } // Track Loop
1421 }
1422
1423 Int_t AliEMCALJetFinder
1424 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1425
1426     Int_t flag    = 0; 
1427     Int_t parton  = 0; 
1428     Int_t neutral = 0;
1429     
1430     if (charge == 0) neutral = 1;
1431     
1432     if (TMath::Abs(code) <= 6   ||      
1433         code == 21              ||   
1434         code == 92) parton = 1;
1435     
1436     //It's not a parton, it's charged and it went through the TPC
1437     if (!parton && !neutral && radius <= 84.0) flag = 1;
1438     
1439     return flag;
1440 }
1441
1442
1443
1444 void AliEMCALJetFinder::SaveBackgroundEvent(Char_t *name)
1445 {
1446 // Saves the eta-phi lego and the tracklist and name of file with BG events
1447 //
1448     if (fLegoB) {
1449        fLegoB->Reset();
1450        (*fLegoB) = (*fLegoB) + (*fLego); 
1451        //       if(fDebug) 
1452        printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n", 
1453        fLegoB->Integral(), fLego->Integral()); 
1454     }
1455     
1456     if (fPtB)        delete[] fPtB;   
1457     if (fEtaB)       delete[] fEtaB;    
1458     if (fPhiB)       delete[] fPhiB;   
1459     if (fPdgB)       delete[] fPdgB;   
1460     if (fTrackListB) delete[] fTrackListB;   
1461    
1462     fPtB          = new Float_t[fNtS];
1463     fEtaB         = new Float_t[fNtS];
1464     fPhiB         = new Float_t[fNtS];
1465     fPdgB         = new Int_t  [fNtS];
1466     fTrackListB   = new Int_t  [fNtS];
1467     
1468     fNtB = 0;
1469     
1470     for (Int_t i = 0; i < fNt; i++) {
1471         if (!fTrackList[i]) continue;
1472         fPtB [fNtB]       = fPtT [i];
1473         fEtaB[fNtB]       = fEtaT[i];
1474         fPhiB[fNtB]       = fPhiT[i];
1475         fPdgB[fNtB]       = fPdgT[i];
1476         fTrackListB[fNtB] = 1;
1477         fNtB++;
1478     }
1479     fBackground = 1;
1480     printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1481
1482     if(strlen(name) == 0) {
1483        TSeqCollection *li = gROOT->GetListOfFiles();
1484        TString nf;
1485        for(Int_t i=0; i<li->GetSize(); i++) {
1486           nf = ((TFile*)li->At(i))->GetName();
1487           if(nf.Contains("backgorund")) break;
1488        }
1489        fBGFileName = nf;
1490     } else {
1491        fBGFileName = name;
1492     }
1493     printf("BG file name is \n %s\n", fBGFileName.Data());
1494 }
1495
1496 void AliEMCALJetFinder::InitFromBackground()
1497 {
1498 //
1499 //    
1500     if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1501     
1502     if (fLego) {
1503         fLego->Reset(); 
1504         (*fLego) = (*fLego) + (*fLegoB);
1505         if(fDebug) 
1506             printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n", 
1507                    fLego->Integral(), fLegoB->Integral()); 
1508     } else {
1509         printf(" => fLego undefined \n");
1510     }
1511 }
1512
1513     
1514 void AliEMCALJetFinder::FindTracksInJetCone()
1515 {
1516 //
1517 //  Build list of tracks inside jet cone
1518 //
1519 //  Loop over jets
1520     Int_t njet = Njets();
1521     for (Int_t nj = 0; nj < njet; nj++)
1522     {
1523         Float_t etaj = JetEtaW(nj);
1524         Float_t phij = JetPhiW(nj);     
1525         Int_t   nT   = 0;           // number of associated tracks
1526         
1527 // Loop over particles in current event 
1528         for (Int_t part = 0; part < fNt; part++) {
1529             if (!fTrackList[part]) continue;
1530             Float_t phi      = fPhiT[part];
1531             Float_t eta      = fEtaT[part];
1532             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1533                                            (phij-phi)*(phij-phi));
1534             if (dr < fConeRadius) {
1535                 fTrackList[part] = nj+2;
1536                 nT++;
1537             } // < ConeRadius ?
1538         } // particle loop
1539         
1540 // Same for background event if available
1541         Int_t nTB = 0;
1542         if (fBackground) {
1543             for (Int_t part = 0; part < fNtB; part++) {
1544                 Float_t phi      = fPhiB[part];
1545                 Float_t eta      = fEtaB[part];
1546                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1547                                                (phij-phi)*(phij-phi));
1548                 fTrackListB[part] = 1;
1549
1550                 if (dr < fConeRadius) {
1551                     fTrackListB[part] = nj+2;
1552                     nTB++;
1553                 } // < ConeRadius ?
1554             } // particle loop
1555         } // Background available ?
1556         
1557         Int_t nT0 = nT + nTB;
1558         printf("Total number of tracks %d\n", nT0);
1559         
1560         if (nT0 > 1000) nT0 = 1000;
1561         
1562         Float_t* ptT  = new Float_t[nT0];
1563         Float_t* etaT = new Float_t[nT0];
1564         Float_t* phiT = new Float_t[nT0];
1565         Int_t*   pdgT = new Int_t[nT0];
1566
1567         Int_t iT = 0;
1568         Int_t j;
1569         
1570         for (Int_t part = 0; part < fNt; part++) {
1571             if (fTrackList[part] == nj+2) {
1572                 Int_t index = 0;
1573                 for (j=iT-1; j>=0; j--) {
1574                     if (fPtT[part] > ptT[j]) {
1575                         index = j+1;
1576                         break;
1577                     }
1578                 }
1579                 for (j=iT-1; j>=index; j--) {
1580                     ptT [j+1]  = ptT [j];
1581                     etaT[j+1]  = etaT[j];
1582                     phiT[j+1]  = phiT[j];
1583                     pdgT[j+1]  = pdgT[j];
1584                 }
1585                 ptT [index] = fPtT [part];
1586                 etaT[index] = fEtaT[part];
1587                 phiT[index] = fPhiT[part];
1588                 pdgT[index] = fPdgT[part];
1589                 iT++;
1590             } // particle associated
1591             if (iT > nT0) break;
1592         } // particle loop
1593         
1594         if (fBackground) {
1595             for (Int_t part = 0; part < fNtB; part++) {
1596                 if (fTrackListB[part] == nj+2) {
1597                     Int_t index = 0;
1598                     for (j=iT-1; j>=0; j--) {
1599                         if (fPtB[part] > ptT[j]) {
1600                             index = j+1;
1601
1602                             break;
1603                         }
1604                     }
1605                     for (j=iT-1; j>=index; j--) {
1606                         ptT [j+1]  = ptT [j];
1607                         etaT[j+1]  = etaT[j];
1608                         phiT[j+1]  = phiT[j];
1609                         pdgT[j+1]  = pdgT[j];
1610                     }
1611                     ptT [index] = fPtB [part];
1612                     etaT[index] = fEtaB[part];
1613                     phiT[index] = fPhiB[part];
1614                     pdgT[index] = fPdgB[part];
1615                     iT++;
1616                 } // particle associated
1617                 if (iT > nT0) break;
1618             } // particle loop
1619         } // Background available ?
1620
1621         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1622         delete[] ptT;
1623         delete[] etaT;
1624         delete[] phiT;
1625         delete[] pdgT;
1626         
1627     } // jet loop loop
1628 }
1629
1630 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1631 {
1632 // Propagates phi angle to EMCAL radius
1633 //
1634   static Float_t b = 0.0, rEMCAL = -1.0;
1635 // Get field in kGS
1636   b = gAlice->Field()->SolenoidField();
1637 // Get EMCAL radius in cm 
1638   rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1639   Float_t dPhi = 0.;
1640 //
1641 //
1642 // bending radies
1643 // pt [Gev]
1644 // B  [kG]
1645 //
1646     Float_t rB = 3335.6 * pt / b;  // [cm]  (case of |charge|=1)
1647 //
1648 // check if particle is curling below EMCAL
1649     if (2.*rB < rEMCAL) {
1650         curls = kTRUE;
1651         return dPhi;
1652     }
1653 //
1654 // if not calculate delta phi
1655     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1656     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1657     dPhi = -TMath::Sign(dPhi, charge);
1658 //    
1659     return dPhi;
1660 }
1661
1662 void hf1(Int_t& , Float_t& , Float_t& )
1663 {
1664 // dummy for hbook calls
1665     ;
1666 }
1667
1668 void AliEMCALJetFinder::DrawLego(Char_t *opt) 
1669 {if(fLego) fLego->Draw(opt);}
1670
1671 void  AliEMCALJetFinder::DrawLegoBackground(Char_t *opt) 
1672 {if(fLegoB) fLegoB->Draw(opt);}
1673
1674 void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt) 
1675 {if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);}
1676
1677 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1678
1679   static TPaveText *varLabel=0;
1680   if(!fC1) {
1681     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1682   }
1683   fC1->Clear();
1684   fC1->Divide(2,2);
1685   fC1->cd(1);
1686   gPad->SetLogy(1);
1687   fhCellEt->Draw();
1688
1689   fC1->cd(2);
1690   gPad->SetLogy(1);
1691   fhCellEMCALEt->Draw();
1692
1693   fC1->cd(3);
1694   gPad->SetLogy(1);
1695   fhTrackPt->Draw();
1696   fhTrackPtBcut->SetLineColor(2);
1697   fhTrackPtBcut->Draw("same");
1698  
1699   fC1->cd(4);
1700   if(!varLabel) {
1701     PrintParameters(1);
1702     varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1703     varLabel->SetTextAlign(12);
1704     varLabel->SetFillColor(19); // see TAttFill
1705     TString tmp(GetTitle());
1706     varLabel->ReadFile(GetFileNameForParameters());
1707   }
1708   varLabel->Draw();
1709   fC1->Update();
1710   if(mode) { // for saving picture to the file
1711     TString stmp(GetFileNameForParameters());
1712     stmp.ReplaceAll("_Par.txt",".ps");
1713     fC1->Print(stmp.Data());
1714   }
1715 }
1716
1717 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1718 {
1719   FILE *file=0;
1720   if(mode==0) file = stdout; // output to terminal
1721   else {
1722     file = fopen(GetFileNameForParameters(),"w");
1723     if(file==0) file = stdout; 
1724   }
1725   fprintf(file,"====   Filling lego   ==== \n");
1726   fprintf(file,"Sampling fraction %6.3f  ", fSamplingF);
1727   fprintf(file,"Smearing          %6i  ", fSmear);
1728   fprintf(file,"Efficiency        %6i\n", fEffic);
1729   fprintf(file,"Hadr.Correct.     %6i  ", fHCorrection);
1730   fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1731   fprintf(file,"====  Jet finding     ==== \n");
1732   fprintf(file,"Cone radius       %6.2f  ", fConeRadius);
1733   fprintf(file,"Seed E_{T}           %6.1f\n", fEtSeed);
1734   fprintf(file,"Min E_{T} of cell    %6.1f  ", fMinCellEt);
1735   fprintf(file,"Min E_{T} of jet     %6.1f\n", fMinJetEt);
1736   if(fMode) {
1737     fprintf(file,"====  Bg subtraction     ==== \n");
1738     fprintf(file,"BG subtraction  %6i  ", fMode);
1739     fprintf(file,"Min cone move   %6.2f\n", fMinMove);
1740     fprintf(file,"Max cone move   %6.2f  ", fMaxMove);
1741     fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1742   } else
1743   fprintf(file,"==== No Bg subtraction     ==== \n");
1744   // 
1745   printf("BG file name is %s \n", fBGFileName.Data());
1746   if(file != stdout) fclose(file); 
1747 }
1748
1749 void AliEMCALJetFinder::DrawLegos()
1750
1751   if(!fC1) {
1752     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1753   }
1754   fC1->Clear();
1755   fC1->Divide(2,2);
1756   gStyle->SetOptStat(111111);
1757
1758   Int_t nent1, nent2, nent3, nent4;
1759   Double_t int1, int2, int3, int4;
1760   nent1 = (Int_t)fLego->GetEntries();
1761   int1  = fLego->Integral();
1762   fC1->cd(1);
1763   if(int1) fLego->Draw("lego");
1764
1765   nent2 = (Int_t)fhLegoTracks->GetEntries();
1766   int2  = fhLegoTracks->Integral();
1767   fC1->cd(2);
1768   if(int2) fhLegoTracks->Draw("lego");
1769
1770   nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1771   int3  = fhLegoEMCAL->Integral();
1772   fC1->cd(3);
1773   if(int3) fhLegoEMCAL->Draw("lego");
1774
1775   nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1776   int4  = fhLegoHadrCorr->Integral();
1777   fC1->cd(4);
1778   if(int4) fhLegoHadrCorr->Draw("lego");
1779
1780   // just for checking 
1781   printf(" Integrals \n");
1782   printf("lego   %10.3f \ntrack  %10.3f \nhits   %10.3f \nHCorr   %10.3f\n--      %10.3f(must be 0)\n", 
1783   int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1784 }
1785
1786 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1787 {
1788   static TString tmp;
1789   if(strlen(dir)) tmp = dir;
1790   tmp += GetTitle();
1791   tmp += "_Par.txt";
1792   return tmp.Data();
1793 }
1794
1795 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1796 { // See FillFromTracks() - npart must be positive
1797     if (fTrackList) delete[] fTrackList;
1798     if (fPtT)       delete[] fPtT;
1799     if (fEtaT)      delete[] fEtaT;
1800     if (fPhiT)      delete[] fPhiT;
1801     if (fPdgT)      delete[] fPdgT;
1802     
1803     if(npart>0) { 
1804        fTrackList = new Int_t  [npart];
1805        fPtT       = new Float_t[npart];
1806        fEtaT      = new Float_t[npart];
1807        fPhiT      = new Float_t[npart];
1808        fPdgT      = new Int_t[npart];
1809     } else {
1810        printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1811     }
1812 }
1813
1814 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1815 {
1816   Int_t absPdg = TMath::Abs(pdg);
1817   if(absPdg<=6) return kTRUE; // quarks
1818   if(pdg == 21) return kTRUE; // gluon 
1819   if(pdg == 92) return kTRUE; // string 
1820
1821   // see p.51 of Pythia Manual
1822   // Not include diquarks with c and b quark - 4-mar-2002
1823   //                 ud_0,sd_0,su_0; dd_1,ud_1,uu_1;  sd_1,su_1,ss_1
1824   static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203,  3103,3203,3303};
1825   for(Int_t i=0; i<9; i++) if(absPdg == diquark[i])  return kTRUE; // diquarks
1826
1827   return kFALSE;
1828 }
1829
1830 void AliEMCALJetFinder::FindChargedJet()
1831 {
1832 //
1833 // Find jet from charged particle information only
1834 //
1835     
1836 //
1837 //  Look for seeds
1838 //
1839     Int_t njets = 0;
1840     Int_t part  = 0;
1841     Int_t nseed = 0;
1842   
1843 //
1844 //
1845     ResetJets();
1846     
1847 //  
1848     for (part = 0; part < fNt; part++) {
1849         if (!fTrackList[part]) continue;
1850         if (fPtT[part] > fEtSeed) nseed++;
1851     }
1852     printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1853     Int_t* iSeeds = new Int_t[nseed];
1854     nseed = 0;
1855     
1856     for (part = 0; part < fNt; part++) {
1857         if (!fTrackList[part]) continue;
1858         if (fPtT[part] > fEtSeed) iSeeds[nseed++] =  part;
1859     }
1860
1861 //
1862 // Loop over seeds
1863 //
1864     Int_t seed = 0;
1865     Float_t pt;
1866     
1867     while(1){
1868 //
1869 // Find seed with highest pt
1870 // 
1871         Float_t ptmax = -1.;
1872         Int_t   index = -1;
1873         Int_t   jndex = -1;
1874         for (seed = 0; seed < nseed; seed++) {
1875             if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1876                 ptmax = pt;
1877                 index = seed;
1878             } // ptmax ?
1879         } // seeds 
1880         if (ptmax < 0.) break;
1881         jndex = iSeeds[index];
1882 //
1883 // Remove from the list   
1884         iSeeds[index] = -1;
1885         printf("\n Next Seed %d %f", jndex, ptmax);
1886 //
1887 // Find tracks in cone around seed
1888 //
1889         Float_t phiSeed = fPhiT[jndex];
1890         Float_t etaSeed = fEtaT[jndex];
1891         Float_t eT = 0.;
1892         Float_t pxJ = 0.;
1893         Float_t pyJ = 0.;
1894         Float_t pzJ = 0.;
1895         
1896         for (part = 0; part < fNt; part++) {
1897             if (!fTrackList[part]) continue;
1898             Float_t deta = fEtaT[part] - etaSeed;
1899             Float_t dphi = fPhiT[part] - phiSeed;
1900             Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
1901             if (dR < fConeRadius) {
1902                 eT += fPtT[part];
1903                 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1904                 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1905                 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1906                 Float_t pz = fPtT[part] / TMath::Tan(theta);
1907                 pxJ += px;
1908                 pyJ += py;
1909                 pzJ += pz;
1910                 //
1911                 // if seed, remove it
1912                 //
1913                 for (seed = 0; seed < nseed; seed++) {
1914                     if (part == iSeeds[seed]) iSeeds[seed] = -1;
1915                 } // seed ?
1916             } // < cone radius
1917         } // particle loop
1918
1919 //
1920 //      Estimate of jet direction
1921         Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
1922         Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1923         Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
1924 //      Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1925         
1926 //
1927 //      Sum up all energy
1928 //
1929         Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1930         Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1931         Int_t dIphi = Int_t(fConeRadius / fDphi);
1932         Int_t dIeta = Int_t(fConeRadius / fDeta);
1933         Int_t iPhi, iEta;
1934         Float_t sumE = 0;
1935         for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1936             for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1937                 if (iPhi < 0 || iEta < 0) continue;
1938                 Float_t dPhi = fPhiMin + iPhi * fDphi;
1939                 Float_t dEta = fEtaMin + iEta * fDeta;
1940                 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1941                 sumE += fLego->GetBinContent(iEta, iPhi);
1942             } // eta
1943         } // phi
1944 //
1945 //
1946 //
1947         fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1948         FindTracksInJetCone();
1949         printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1950         printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1951         printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1952     } // while(1)
1953     EMCALJETS.njet = njets;
1954     if (fWrite) WriteJets();
1955     fEvent++;
1956 }
1957 // 16-jan-2003 - just for convenience
1958 void AliEMCALJetFinder::Browse(TBrowser* b)
1959 {
1960    if(fHistsList)  b->Add((TObject*)fHistsList);
1961 }
1962
1963 Bool_t AliEMCALJetFinder::IsFolder() const
1964 {
1965   if(fHistsList) return kTRUE;
1966   else           return kFALSE;
1967 }
1968
1969 const Char_t* AliEMCALJetFinder::GetNameOfVariant()
1970 {// generate the literal string with info about jet finder
1971   Char_t name[200];
1972   sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
1973           fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
1974   TString nt(name);
1975   nt.ReplaceAll(" ","");
1976   if(fBGFileName.Length()) {
1977      Int_t i1 = fBGFileName.Index("kBackground");
1978      Int_t i2 = fBGFileName.Index("/0000") - 1;
1979      if(i1>=0 && i2>=0) {
1980         TString bg(fBGFileName(i1,i2-i1+1));
1981         nt += bg;
1982      }
1983   }
1984   printf("<I> Name of variant %s \n", nt.Data());
1985   return nt.Data();
1986 }