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