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