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