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