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