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