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