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