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