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