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