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