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