a5c3b95f6035524082a23d287823806417c19831
[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.16  2002/02/08 11:43:05  morsch
19 SetOutputFileName(..) allows to specify an output file to which the
20 reconstructed jets are written. If not set output goes to input file.
21 Attention call Init() before processing.
22
23 Revision 1.15  2002/02/02 08:37:18  morsch
24 Formula for rB corrected.
25
26 Revision 1.14  2002/02/01 08:55:30  morsch
27 Fill map with Et and pT.
28
29 Revision 1.13  2002/01/31 09:37:36  morsch
30 Geometry parameters in constructor and call of SetCellSize()
31
32 Revision 1.12  2002/01/23 13:40:23  morsch
33 Fastidious debug print statement removed.
34
35 Revision 1.11  2002/01/22 17:25:47  morsch
36 Some corrections for event mixing and bg event handling.
37
38 Revision 1.10  2002/01/22 10:31:44  morsch
39 Some correction for bg mixing.
40
41 Revision 1.9  2002/01/21 16:28:42  morsch
42 Correct sign of dphi.
43
44 Revision 1.8  2002/01/21 16:05:31  morsch
45 - different phi-bin for hadron correction
46 - provisions for background mixing.
47
48 Revision 1.7  2002/01/21 15:47:18  morsch
49 Bending radius correctly in cm.
50
51 Revision 1.6  2002/01/21 12:53:50  morsch
52 authors
53
54 Revision 1.5  2002/01/21 12:47:47  morsch
55 Possibility to include K0long and neutrons.
56
57 Revision 1.4  2002/01/21 11:03:21  morsch
58 Phi propagation introduced in FillFromTracks.
59
60 Revision 1.3  2002/01/18 05:07:56  morsch
61 - hadronic correction
62 - filling of digits
63 - track selection upon EMCAL information
64
65 */
66
67 //*-- Authors: Andreas Morsch   (CERN)
68 //*            J.L. Klay        (LBL)
69 //*            Aleksei Pavlinov (WSU) 
70
71 // From root ...
72 #include <TClonesArray.h>
73 #include <TTree.h>
74 #include <TBranchElement.h>
75 #include <TFile.h>
76 #include <TH2.h>
77 #include <TAxis.h>
78 #include <TParticle.h>
79 #include <TParticlePDG.h>
80
81 // From AliRoot ...
82 #include "AliEMCALJetFinder.h"
83 #include "AliEMCALFast.h"
84 #include "AliEMCALGeometry.h"
85 #include "AliEMCALHit.h"
86 #include "AliEMCALDigit.h"
87 #include "AliEMCALDigitizer.h"
88 #include "AliEMCALHadronCorrection.h"
89 #include "AliRun.h"
90 #include "AliMagF.h"
91 #include "AliMagFCM.h"
92 #include "AliEMCAL.h"
93 #include "AliHeader.h"
94 #include "AliPDG.h"
95
96 // Interface to FORTRAN
97 #include "Ecommon.h"
98
99
100 ClassImp(AliEMCALJetFinder)
101
102 //____________________________________________________________________________
103 AliEMCALJetFinder::AliEMCALJetFinder()
104 {
105 // Default constructor
106     fJets             = 0;
107     fNjets            = 0;
108     fLego             = 0;
109     fLegoB            = 0;
110     fTrackList        = 0;
111     fPtT              = 0;
112     fEtaT             = 0;
113     fPhiT             = 0;
114     fPtB              = 0;
115     fEtaB             = 0;
116     fPhiB             = 0;
117     fHCorrection      = 0;
118     fHadronCorrector  = 0;
119
120     fOutFileName      = 0;
121     fOutFile          = 0;
122     fInFile           = 0;
123     fEvent            = 0;
124
125     SetParametersForBgSubtraction();
126 }
127
128 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
129     : TTask(name, title)
130 {
131 // Constructor 
132     fJets  = new TClonesArray("AliEMCALJet",10000);
133     fNjets = 0;
134     for (Int_t i = 0; i < 30000; i++)
135     {
136         fEtCell[i]  = 0.;
137         fEtaCell[i] = 0.;
138         fPhiCell[i] = 0.;
139     }
140     fLego = 0;
141     fTrackList  = 0;
142     fPtT        = 0;
143     fEtaT       = 0;
144     fPhiT       = 0;
145     fPtB        = 0;
146     fEtaB       = 0;
147     fPhiB       = 0;
148     fHCorrection      = 0;
149     fHadronCorrector  = 0;
150     fBackground       = 0;
151     fOutFileName      = 0;
152     fOutFile          = 0;
153     fInFile           = 0;
154     fEvent            = 0;
155 //
156     SetPtCut();
157     SetMomentumSmearing();
158     SetEfficiencySim();
159     SetDebug();
160     SetHadronCorrection();
161     SetSamplingFraction();
162     SetIncludeK0andN();
163
164     SetParametersForBgSubtraction();
165 }
166
167 void AliEMCALJetFinder::SetParametersForBgSubtraction
168 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
169 {
170 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
171 // at WSU Linux cluster - 11-feb-2002
172 // These parameters must be tuned more carefull !!!
173   SetMode(mode);
174   SetMinMove(minMove);
175   SetMaxMove(maxMove);
176   SetPrecBg(precBg);
177 }
178
179 //____________________________________________________________________________
180 AliEMCALJetFinder::~AliEMCALJetFinder()
181 {
182 // Destructor
183 //
184     if (fJets){
185         fJets->Delete();
186         delete fJets;
187     }
188 }
189
190
191
192
193 #ifndef WIN32
194 # define jet_finder_ua1 jet_finder_ua1_
195 # define hf1 hf1_
196 # define type_of_call
197
198 #else
199 # define jet_finder_ua1 J
200 # define hf1 HF1
201 # define type_of_call _stdcall
202 #endif
203
204 extern "C" void type_of_call 
205 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
206                Float_t  etc[30000],  Float_t etac[30000],
207                Float_t  phic[30000], 
208                Float_t& min_move, Float_t& max_move, Int_t& mode, 
209                Float_t& prec_bg,  Int_t& ierror);
210
211 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
212
213
214
215
216 void AliEMCALJetFinder::Init()
217 {
218 //
219 // Geometry and I/O initialization
220 //
221 //
222 //
223 //  Get geometry parameters from EMCAL
224 //
225 //
226 //  Geometry 
227     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
228     AliEMCALGeometry* geom = 
229         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
230     fNbinEta = geom->GetNZ();
231     fNbinPhi = geom->GetNPhi();
232     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
233     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
234     fEtaMin  = geom->GetArm1EtaMin();
235     fEtaMax  = geom->GetArm1EtaMax();
236     fDphi    = (fPhiMax-fPhiMin)/fNbinEta;
237     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
238     fNtot    = fNbinPhi*fNbinEta;
239 //
240     SetCellSize(fDeta, fDphi);
241 //
242 //  I/O
243     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
244 }
245
246 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], 
247                              Float_t etac[30000], Float_t phic[30000],
248                              Float_t min_move, Float_t max_move, Int_t mode, 
249                              Float_t prec_bg,  Int_t   ierror)
250 {
251 // Wrapper for fortran coded jet finder
252 // Full argument list
253     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic, 
254                    min_move, max_move, mode, prec_bg, ierror);
255     // Write result to output
256     if(fWrite) WriteJets();
257     fEvent++;
258 }
259
260 void AliEMCALJetFinder::Find()
261 {
262 // Wrapper for fortran coded jet finder using member data for 
263 // argument list
264
265     Float_t min_move = fMinMove;
266     Float_t max_move = fMaxMove;
267     Int_t   mode     = fMode;
268     Float_t prec_bg  = fPrecBg;
269     Int_t   ierror;
270
271     ResetJets(); // 4-feb-2002 by PAI
272
273     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, 
274                    min_move, max_move, mode, prec_bg, ierror);
275     fError = ierror;
276     // Write result to output
277     Int_t njet = Njets();
278     
279     for (Int_t nj=0; nj<njet; nj++)
280     {
281         
282         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
283                                     JetPhiW(nj),
284                                     JetEtaW(nj));
285     }
286
287     if(fNt)    FindTracksInJetCone();
288     if(fWrite) WriteJets();
289     fEvent++;
290 }
291
292 Int_t AliEMCALJetFinder::Njets()
293 {
294 // Get number of reconstructed jets
295     return EMCALJETS.njet;
296 }
297
298 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
299 {
300 // Get reconstructed jet energy
301     return EMCALJETS.etj[i];
302 }
303
304 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
305 {
306 // Get reconstructed jet phi from leading particle
307     return EMCALJETS.phij[0][i];
308 }
309
310 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
311 {
312 // Get reconstructed jet phi from weighting
313     return EMCALJETS.phij[1][i];
314 }
315
316 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
317 {
318 // Get reconstructed jet eta from leading particles
319     return EMCALJETS.etaj[0][i];
320 }
321
322
323 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)  
324 {
325 // Get reconstructed jet eta from weighting
326     return EMCALJETS.etaj[1][i];
327 }
328
329 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
330 {
331 // Set grid cell size
332     EMCALCELLGEO.etaCellSize = eta;
333     EMCALCELLGEO.phiCellSize = phi;    
334 }
335
336 void AliEMCALJetFinder::SetConeRadius(Float_t par)
337 {
338 // Set jet cone radius
339     EMCALJETPARAM.coneRad = par;
340     fConeRadius = par;
341 }
342
343 void AliEMCALJetFinder::SetEtSeed(Float_t par)
344 {
345 // Set et cut for seeds
346     EMCALJETPARAM.etSeed = par;
347     fEtSeed = par;
348 }
349
350 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
351 {
352 // Set minimum jet et
353     EMCALJETPARAM.ejMin = par;
354     fMinJetEt = par;
355 }
356
357 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
358 {
359 // Set et cut per cell
360     EMCALJETPARAM.etMin = par;
361     fMinCellEt = par;
362 }
363
364 void AliEMCALJetFinder::SetPtCut(Float_t par)
365 {
366 // Set pT cut on charged tracks
367     fPtCut = par;
368 }
369
370
371 void AliEMCALJetFinder::Test()
372 {
373 //
374 // Test the finder call
375 //
376     const Int_t nmax = 30000;
377     Int_t ncell      = 10;
378     Int_t ncell_tot  = 100;
379
380     Float_t etc[nmax];
381     Float_t etac[nmax];
382     Float_t phic[nmax];
383     Float_t min_move = 0;
384     Float_t max_move = 0;
385     Int_t   mode = 0;
386     Float_t prec_bg = 0;
387     Int_t   ierror = 0;
388
389     
390     Find(ncell, ncell_tot, etc, etac, phic, 
391          min_move, max_move, mode, prec_bg, ierror);
392
393 }
394
395 //
396 //  I/O
397 //      
398
399 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
400 {
401     //
402     // Add a jet 
403     //
404     TClonesArray &lrawcl = *fJets;
405     new(lrawcl[fNjets++]) AliEMCALJet(jet);
406 }
407
408 void AliEMCALJetFinder::ResetJets()
409 {
410     //
411     // Reset Jet List 
412     //
413     fJets->Clear();
414     fNjets = 0;
415 }
416
417 void AliEMCALJetFinder::WriteJets()
418 {
419 //
420 // Add all jets to the list
421 //
422     const Int_t kBufferSize = 4000;
423     const char* file = 0;
424
425     Int_t njet = Njets();
426
427     for (Int_t nj = 0; nj < njet; nj++)
428     {
429         AddJet(*fJetT[nj]);
430         delete fJetT[nj];
431     }
432
433 // I/O
434     if (!fOutFileName) {
435 //
436 // output written to input file
437 //
438         AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
439         TTree* pK = gAlice->TreeK();
440         file = (pK->GetCurrentFile())->GetName();
441         if (fDebug > 1)
442             printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
443         if (fJets && gAlice->TreeR()) {
444             pEMCAL->MakeBranchInTree(gAlice->TreeR(), 
445                                      "EMCALJets", 
446                                      &fJets, 
447                                      kBufferSize, 
448                                      file);
449         }
450         Int_t nev = gAlice->GetHeader()->GetEvent();
451         gAlice->TreeR()->Fill();
452         char hname[30];
453         sprintf(hname,"TreeR%d", nev);
454         gAlice->TreeR()->Write(hname);
455         gAlice->TreeR()->Reset();
456     } else {
457 //
458 // Output written to user specified output file
459 //
460         TTree* pK = gAlice->TreeK();
461         fInFile  = pK->GetCurrentFile();
462
463         fOutFile->cd();
464         char hname[30];
465         sprintf(hname,"TreeR%d", fEvent);
466         TTree* treeJ = new TTree(hname, "EMCALJets");
467         treeJ->Branch("EMCALJets", &fJets, kBufferSize);
468         treeJ->Fill();
469         treeJ->Write(hname);
470         fInFile->cd();
471     }
472     ResetJets();        
473 }
474
475 void AliEMCALJetFinder::BookLego()
476 {
477 //
478 //  Book histo for discretisation
479 //
480
481 //
482 //  Don't add histos to the current directory
483     TH2::AddDirectory(0);
484 //    
485 //  Signal map
486     fLego = new TH2F("legoH","eta-phi",
487                            fNbinEta, fEtaMin, fEtaMax, 
488                            fNbinPhi, fPhiMin, fPhiMax);
489 //
490 //  Background map
491     fLegoB = new TH2F("legoB","eta-phi",
492                            fNbinEta, fEtaMin, fEtaMax, 
493                            fNbinPhi, fPhiMin, fPhiMax);
494     
495 }
496
497 void AliEMCALJetFinder::DumpLego()
498 {
499 //
500 // Dump lego histo into array
501 //
502     fNcell = 0;
503     TAxis* Xaxis = fLego->GetXaxis();
504     TAxis* Yaxis = fLego->GetYaxis();
505     for (Int_t i = 1; i < fNbinEta; i++) {
506         for (Int_t j = 1; j < fNbinPhi; j++) {
507             Float_t e    = fLego->GetBinContent(i,j);
508             if (e <=0.) continue;
509             
510             Float_t eta  = Xaxis->GetBinCenter(i);
511             Float_t phi  = Yaxis->GetBinCenter(j);          
512             fEtCell[fNcell]  = e;
513             fEtaCell[fNcell] = eta;
514             fPhiCell[fNcell] = phi;
515             fNcell++;
516         } // phi
517     } // eta
518     fNcell--;
519 }
520
521 void AliEMCALJetFinder::ResetMap()
522 {
523 //
524 // Reset eta-phi array
525
526     for (Int_t i=0; i<30000; i++)
527     {
528         fEtCell[i]  = 0.;
529         fEtaCell[i] = 0.;
530         fPhiCell[i] = 0.;
531     }
532 }
533
534
535 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
536 {
537 //
538 // Fill Cells with track information
539 //
540     if (fDebug >= 2)
541     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
542
543     ResetMap();
544     
545     if (!fLego) BookLego();
546 // Reset
547     if (flag == 0) fLego->Reset();
548 //
549 // Access particle information    
550     Int_t npart = (gAlice->GetHeader())->GetNprimary();
551     if (fDebug >= 2 || npart<=0) printf(" : npart %i\n", npart);
552  
553 // Create track list
554 //
555 // 0: not selected
556 // 1: selected for jet finding
557 // 2: inside jet
558 // ....
559     if (fTrackList) delete[] fTrackList;
560     if (fPtT)       delete[] fPtT;
561     if (fEtaT)      delete[] fEtaT;
562     if (fPhiT)      delete[] fPhiT;
563     
564     fTrackList = new Int_t  [npart];
565     fPtT       = new Float_t[npart];
566     fEtaT      = new Float_t[npart];
567     fPhiT      = new Float_t[npart];
568
569     fNt        = npart;
570     fNtS  = 0;
571     Float_t chTmp=0.0; // charge of current particle 
572
573     for (Int_t part = 2; part < npart; part++) {
574         TParticle *MPart = gAlice->Particle(part);
575         Int_t mpart   = MPart->GetPdgCode();
576         Int_t child1  = MPart->GetFirstDaughter();
577         Float_t pT    = MPart->Pt();
578         Float_t p     = MPart->P();
579         Float_t phi   = MPart->Phi();
580         Float_t eta   = MPart->Eta();
581         Float_t theta = MPart->Theta(); 
582         
583         if (fDebug > 1) {
584             if (part == 6 || part == 7)
585             {
586                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n", 
587                        part-5, pT, eta, phi);
588             }
589
590 //          if (mpart == 21)
591                 
592 //              printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
593 //                     part, mpart, pT, eta, phi); 
594         }
595         
596         
597         fTrackList[part] = 0; 
598         fPtT[part]       = pT; // must be change after correction for resolution !!!
599         fEtaT[part]      = eta;
600         fPhiT[part]      = phi;
601
602         if (part < 2) continue;
603
604         // move to fLego->Fill because hadron correction must apply 
605         // if particle hit to EMCAL - 11-feb-2002
606         //      if (pT == 0 || pT < fPtCut) continue;
607         TParticlePDG* pdgP = 0;
608 // charged or neutral 
609         pdgP  = MPart->GetPDG();
610         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!  
611         if (ich == 0) {
612             if (chTmp == 0) {
613                 if (!fK0N) { 
614                     continue;
615                 } else {
616                     if (mpart != kNeutron    &&
617                         mpart != kNeutronBar &&
618                         mpart != kK0Long) continue;
619                 }
620             }
621         }
622 // skip partons
623         if (TMath::Abs(mpart) <= 6         ||
624             mpart == 21                    ||
625             mpart == 92) continue;
626 // acceptance cut
627         if (TMath::Abs(eta) > 0.7)         continue;
628 //   Initial phi may be out of acceptance but track may 
629 //   hit two the acceptance  - see variable curls
630 //      if (phi*180./TMath::Pi() > 120.)   continue;
631 // final state only
632         if (child1 >= 0 && child1 < npart) continue;
633 //
634         if (fDebug > 1) 
635         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
636         part, mpart, child1, eta, phi, pT, chTmp);
637 //
638 //
639 // Momentum smearing goes here ...
640 //
641         Float_t pw;
642         if (fSmear && TMath::Abs(chTmp)) {
643             pw = AliEMCALFast::SmearMomentum(1,p);
644         // p changed - take into account when calculate pT,
645         // pz and so on ..  ?
646             pT = (pw/p) * pT;
647             if(fDebug > 1) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
648             p  = pw;
649         }
650 //
651 // Tracking Efficiency and TPC acceptance goes here ...
652         if (fEffic && TMath::Abs(chTmp)) {
653             Float_t eff =  AliEMCALFast::Efficiency(1,p);
654             if (AliEMCALFast::RandomReject(eff)) {
655               if(fDebug > 1) printf(" reject due to unefficiency ");
656               continue;
657             }
658         }
659 //
660 // Correction of Hadronic Energy goes here ...
661 //
662 //
663 // phi propagation for hadronic correction
664
665         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
666         Float_t phiHC=0.0, dpH=0.0, dphi=0.0;
667         if(TMath::Abs(chTmp)) {
668         // hadr. correction only for charge particle 
669           dphi  = PropagatePhi(pT, chTmp, curls);
670           phiHC = phi + dphi;
671           if (fDebug >= 2) {
672              printf("\n Delta phi %f pT %f ", dphi, pT);
673              if (curls) printf("\n !! Track is curling");
674           }
675         
676           if (fHCorrection && !curls) {
677               if (!fHadronCorrector)
678                  Fatal("AliEMCALJetFinder",
679                     "Hadronic energy correction required but not defined !");
680
681               dpH = fHadronCorrector->GetEnergy(p, eta, 7);
682
683               if (fDebug >= 2) printf(" phi %f phiHC %f eTcorr %f\n", 
684               phi, phiHC, -dpH*TMath::Sin(theta));   
685               fLego->Fill(eta, phiHC, -dpH*TMath::Sin(theta));
686           }
687         }
688 //
689 //  More to do ? Just think about it !
690 //
691         
692         if (phi*180./TMath::Pi() > 120.)   continue;
693
694         if(TMath::Abs(chTmp) ) { // charge particle
695           if (pT > fPtCut && !curls) {
696              if (fDebug >= 2) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
697                                       eta , phi, pT); 
698              fLego->Fill(eta, phi, pT);
699              fTrackList[part] = 1;
700              fNtS++;
701           }
702         } else if(ich == 0) {
703           // case of n, nbar and K0L
704              if (fDebug >= 2) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
705                                       eta , phi, pT); 
706             fLego->Fill(eta, phi, pT);
707             fTrackList[part] = 1;
708             fNtS++;
709         }
710
711     } // primary loop
712     DumpLego();
713 }
714
715 void AliEMCALJetFinder::FillFromHits(Int_t flag)
716 {
717 //
718 // Fill Cells with hit information
719 //
720 //
721     if (fDebug >= 2)
722     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
723
724     ResetMap();
725     
726     if (!fLego) BookLego();
727 //  Reset eta-phi map if needed
728     if (flag == 0)    fLego->Reset();
729 //  Initialize from background event if available
730 //
731 // Access hit information    
732     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
733     
734     TTree *treeH = gAlice->TreeH();
735     Int_t ntracks = (Int_t) treeH->GetEntries();
736 //
737 //   Loop over tracks
738 //
739     Int_t nbytes = 0;
740
741     
742     for (Int_t track=0; track<ntracks;track++) {
743         gAlice->ResetHits();
744         nbytes += treeH->GetEvent(track);
745 //
746 //  Loop over hits
747 //
748         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); 
749             mHit;
750             mHit=(AliEMCALHit*) pEMCAL->NextHit()) 
751         {
752             Float_t x      =    mHit->X();         // x-pos of hit
753             Float_t y      =    mHit->Y();         // y-pos
754             Float_t z      =    mHit->Z();         // z-pos
755             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
756
757             Float_t r      =    TMath::Sqrt(x*x+y*y);
758             Float_t theta  =    TMath::ATan2(r,z);
759             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
760             Float_t phi    =    TMath::ATan2(y,x);
761
762             if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
763             
764             fLego->Fill(eta, phi, fSamplingF*eloss*TMath::Sin(theta));
765         } // Hit Loop
766     } // Track Loop
767     DumpLego();
768 }
769
770 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
771 {
772 //
773 // Fill Cells with digit information
774 //
775 //
776     ResetMap();
777     
778     if (!fLego) BookLego();
779     if (flag == 0) fLego->Reset();
780     Int_t nbytes;
781     
782
783 //
784 //  Connect digits    
785 //
786     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
787     TTree *treeD = gAlice->TreeD();
788     TBranchElement* branchDg = (TBranchElement*)
789         treeD->GetBranch("EMCAL");
790
791     if (!branchDg) Fatal("AliEMCALJetFinder", 
792                          "Reading digits requested but no digits in file !");
793     
794     branchDg->SetAddress(&digs);
795     Int_t nent = (Int_t) branchDg->GetEntries();
796 //
797 //  Connect digitizer
798 //
799     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
800     TBranchElement* branchDr = (TBranchElement*) 
801         treeD->GetBranch("AliEMCALDigitizer");
802     branchDr->SetAddress(&digr);
803 //
804 //
805     nbytes += branchDg->GetEntry(0);
806     nbytes += branchDr->GetEntry(0);
807 //
808 //  Get digitizer parameters
809     Float_t towerADCped  = digr->GetTowerpedestal();
810     Float_t towerADCcha  = digr->GetTowerchannel();
811     Float_t preshoADCped = digr->GetPreShopedestal();
812     Float_t preshoADCcha = digr->GetPreShochannel();
813
814     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
815     AliEMCALGeometry* geom = 
816         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
817     
818     if (fDebug) {
819         Int_t   ndig = digs->GetEntries();
820         printf("\n Number of Digits: %d %d\n", ndig, nent);
821         printf("\n Parameters: %f %f %f %f\n", 
822                towerADCped, towerADCcha, preshoADCped, preshoADCcha );
823         printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
824     }
825     
826 //
827 //  Loop over digits
828     AliEMCALDigit* sdg;
829     TIter next(digs);
830     while ((sdg = (AliEMCALDigit*)(next())))
831     {
832         Double_t pedestal;
833         Double_t channel;
834         if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi())) 
835         {
836             pedestal = preshoADCped;
837             channel  = preshoADCcha; 
838         } else {
839             pedestal = towerADCped;
840             channel  = towerADCcha; 
841         }
842         
843         Float_t eta = sdg->GetEta();
844         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
845         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
846         
847         if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
848                            eta, phi, amp, sdg->GetAmp());
849         
850         fLego->Fill(eta, phi, fSamplingF*amp);
851     } // digit loop
852 // 
853 //  Dump lego hist
854     DumpLego();
855 }
856
857
858 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
859 {
860 //
861 // Fill Cells with hit information
862 //
863 //
864     ResetMap();
865     
866     if (!fLego) BookLego();
867 // Reset
868     if (flag == 0) fLego->Reset();
869
870 // Flag charged tracks passing through TPC which 
871 // are associated to EMCAL Hits
872     BuildTrackFlagTable();
873
874 //
875 // Access particle information    
876     TTree *treeH = gAlice->TreeH();
877     Int_t ntracks = (Int_t) treeH->GetEntries();
878
879     if (fPtT)       delete[] fPtT;   
880     if (fEtaT)      delete[] fEtaT;    
881     if (fPhiT)      delete[] fPhiT;   
882    
883     fPtT       = new Float_t[ntracks];
884     fEtaT      = new Float_t[ntracks];
885     fPhiT      = new Float_t[ntracks];
886
887     fNt   = ntracks;
888     fNtS  = 0;
889     
890     for (Int_t track = 0; track < ntracks; track++) {
891           TParticle *MPart = gAlice->Particle(track);
892           Float_t pT    = MPart->Pt();
893           Float_t phi   = MPart->Phi();
894           Float_t eta   = MPart->Eta();
895
896         if(fTrackList[track]) {
897           fPtT[track]       = pT;
898           fEtaT[track]      = eta;
899           fPhiT[track]      = phi;
900
901           if (track < 2) continue;      //Colliding particles?
902           if (pT == 0 || pT < fPtCut) continue;
903           fNtS++;
904           fLego->Fill(eta, phi, pT);
905         }
906       } // track loop
907     DumpLego();
908 }
909
910 void AliEMCALJetFinder::FillFromPartons()
911 {
912 // function under construction - 13-feb-2002 PAI
913     
914     if (fDebug >= 2)
915     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
916     // 
917
918     ResetMap();
919     if (!fLego) BookLego();
920     fLego->Reset();
921 //
922 // Access particle information    
923     Int_t npart = (gAlice->GetHeader())->GetNprimary();
924     if (fDebug >= 2 || npart<=0)
925     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
926     fNt   = 0; // for FindTracksInJetCone
927     fNtS  = 0;
928  
929 //  Go through the partons
930     Int_t statusCode=0;
931     for (Int_t part = 8; part < npart; part++) {
932         TParticle *MPart = gAlice->Particle(part);
933         Int_t mpart   = MPart->GetPdgCode();
934         //      Int_t child1  = MPart->GetFirstDaughter();
935         Float_t pT    = MPart->Pt();
936         //      Float_t p     = MPart->P();
937         Float_t phi   = MPart->Phi();
938         Float_t eta   = MPart->Eta();
939         //      Float_t theta = MPart->Theta();
940         statusCode    = MPart->GetStatusCode();
941         
942 // accept partons (21 - gluon, 92 - string) 
943         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
944         if (fDebug > 1 && pT>0.01) 
945         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
946         part, mpart, statusCode, eta, phi, pT);
947         //      if (fDebug >= 3) MPart->Print(); 
948 // accept partons before fragmentation - p.57 in Pythia manual
949 //        if(statusCode != 1) continue;
950 // acceptance cut
951         if (TMath::Abs(eta) > 0.7)         continue;
952         if (phi*180./TMath::Pi() > 120.)   continue;
953 // final state only
954 //      if (child1 >= 0 && child1 < npart) continue;
955 //
956 //
957         fLego->Fill(eta, phi, pT);
958
959     } // primary loop
960     DumpLego();
961 }
962
963 void AliEMCALJetFinder::BuildTrackFlagTable() {
964
965 // Method to generate a lookup table for TreeK particles
966 // which are linked to hits in the EMCAL
967 //
968 // --Author: J.L. Klay
969 //
970 // Access hit information    
971     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
972     
973     TTree *TK = gAlice->TreeK();                // Get the number of entries in the kine tree
974     Int_t nKTrks = (Int_t) TK->GetEntries();    // (Number of particles created somewhere)
975     
976     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
977     fTrackList = new Int_t[nKTrks];             //before generating a new one
978     
979     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
980         fTrackList[i] = 0;
981     }
982     
983     TTree *treeH = gAlice->TreeH();
984     Int_t ntracks = (Int_t) treeH->GetEntries();
985 //
986 //   Loop over tracks
987 //
988     Int_t nbytes = 0;
989     
990     for (Int_t track=0; track<ntracks;track++) {
991         gAlice->ResetHits();
992         nbytes += treeH->GetEvent(track);
993         if (pEMCAL)  {
994             
995 //
996 //  Loop over hits
997 //
998             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); 
999                 mHit;
1000                 mHit=(AliEMCALHit*) pEMCAL->NextHit()) 
1001             {
1002                 Int_t   iTrk   =    mHit->Track();       // track number
1003                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1004                 
1005                 //Determine the origin point of this particle - it made a hit in the EMCAL
1006                 TParticle         *trkPart = gAlice->Particle(iTrk);
1007                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1008                 Int_t      trkCode = trkPart->GetPdgCode();
1009                 Double_t           trkChg;
1010                 if (trkCode < 10000) {          //Big Ions cause problems for 
1011                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1012                 } else {                                //likely to make it very far, set
1013                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1014                 }
1015                 Float_t            vxTrk   = trkPart->Vx();
1016                 Float_t            vyTrk   = trkPart->Vy();
1017                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk); 
1018                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1019                 
1020                 //Loop through the ancestry of the EMCAL entrance particles
1021                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1022                 while (ancestor != -1) {
1023                     TParticle     *ancPart = gAlice->Particle(ancestor);  //get particle info on ancestor
1024                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1025                     Int_t             ancCode = ancPart->GetPdgCode();
1026                     Double_t       ancChg;
1027                     if (ancCode < 10000) {
1028                         ancChg  = ancPdg->Charge();
1029                     } else {
1030                         ancChg  = 0.0;
1031                     }
1032                     Float_t           vxAnc   = ancPart->Vx();
1033                     Float_t           vyAnc   = ancPart->Vy();
1034                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1035                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1036                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1037                 }
1038                 
1039                 //Determine the origin point of the primary particle associated with the hit
1040                 TParticle     *primPart = gAlice->Particle(idprim);
1041                 TParticlePDG  *primPdg  = primPart->GetPDG();
1042                 Int_t      primCode = primPart->GetPdgCode();
1043                 Double_t           primChg;
1044                 if (primCode < 10000) {
1045                     primChg  = primPdg->Charge();
1046                 } else {
1047                     primChg  = 0.0;
1048                 }
1049                 Float_t    vxPrim = primPart->Vx();
1050                 Float_t    vyPrim = primPart->Vy();
1051                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1052                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1053             } // Hit Loop
1054         } //if (pEMCAL)
1055     } // Track Loop
1056 }
1057
1058 Int_t AliEMCALJetFinder
1059 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1060
1061     Int_t flag    = 0; 
1062     Int_t parton  = 0; 
1063     Int_t neutral = 0;
1064     
1065     if (charge == 0) neutral = 1;
1066     
1067     if (TMath::Abs(code) <= 6   ||      
1068         code == 21              ||   
1069         code == 92) parton = 1;
1070     
1071     //It's not a parton, it's charged and it went through the TPC
1072     if (!parton && !neutral && radius <= 84.0) flag = 1;
1073     
1074     return flag;
1075 }
1076
1077
1078
1079 void AliEMCALJetFinder::SaveBackgroundEvent()
1080 {
1081 // Saves the eta-phi lego and the tracklist
1082 //
1083     if (fLegoB) fLegoB->Reset();
1084     
1085     fLego->Copy(*fLegoB);
1086     
1087     if (fPtB)        delete[] fPtB;   
1088     if (fEtaB)       delete[] fEtaB;    
1089     if (fPhiB)       delete[] fPhiB;   
1090     if (fTrackListB) delete[] fTrackListB;   
1091    
1092     fPtB          = new Float_t[fNtS];
1093     fEtaB         = new Float_t[fNtS];
1094     fPhiB         = new Float_t[fNtS];
1095     fTrackListB   = new Int_t  [fNtS];
1096     
1097     fNtB = 0;
1098     
1099     for (Int_t i = 0; i < fNt; i++) {
1100         if (!fTrackList[i]) continue;
1101         fPtB [fNtB]       = fPtT [i];
1102         fEtaB[fNtB]       = fEtaT[i];
1103         fPhiB[fNtB]       = fPhiT[i];
1104         fTrackListB[fNtB] = 1;
1105         fNtB++;
1106     }
1107     fBackground = 1;
1108 }
1109
1110 void AliEMCALJetFinder::InitFromBackground()
1111 {
1112 //
1113 //    
1114     if (fDebug) printf("\n Initializing from Background");
1115     
1116     if (fLego) fLego->Reset(); 
1117     fLegoB->Copy(*fLego);
1118 }
1119
1120     
1121 void AliEMCALJetFinder::FindTracksInJetCone()
1122 {
1123 //
1124 //  Build list of tracks inside jet cone
1125 //
1126 //  Loop over jets
1127     Int_t njet = Njets();
1128     for (Int_t nj = 0; nj < njet; nj++)
1129     {
1130         Float_t etaj = JetEtaW(nj);
1131         Float_t phij = JetPhiW(nj);     
1132         Int_t   nT   = 0;           // number of associated tracks
1133         
1134 // Loop over particles in current event 
1135         for (Int_t part = 0; part < fNt; part++) {
1136             if (!fTrackList[part]) continue;
1137             Float_t phi      = fPhiT[part];
1138             Float_t eta      = fEtaT[part];
1139             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1140                                            (phij-phi)*(phij-phi));
1141             if (dr < fConeRadius) {
1142                 fTrackList[part] = nj+2;
1143                 nT++;
1144             } // < ConeRadius ?
1145         } // particle loop
1146         
1147 // Same for background event if available
1148         Int_t nTB = 0;
1149         if (fBackground) {
1150             for (Int_t part = 0; part < fNtB; part++) {
1151                 Float_t phi      = fPhiB[part];
1152                 Float_t eta      = fEtaB[part];
1153                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1154                                                (phij-phi)*(phij-phi));
1155                 fTrackListB[part] = 1;
1156
1157                 if (dr < fConeRadius) {
1158                     fTrackListB[part] = nj+2;
1159                     nTB++;
1160                 } // < ConeRadius ?
1161             } // particle loop
1162         } // Background available ?
1163         
1164         Int_t nT0 = nT + nTB;
1165         
1166         if (nT0 > 50) nT0 = 50;
1167         
1168         Float_t* ptT  = new Float_t[nT0];
1169         Float_t* etaT = new Float_t[nT0];
1170         Float_t* phiT = new Float_t[nT0];
1171         Int_t iT = 0;
1172         Int_t j;
1173         
1174         for (Int_t part = 0; part < fNt; part++) {
1175             if (fTrackList[part] == nj+2) {
1176                 Int_t index = 0;
1177                 for (j=iT-1; j>=0; j--) {
1178                     if (fPtT[part] > ptT[j]) {
1179                         index = j+1;
1180                         break;
1181                     }
1182                 }
1183                 for (j=iT-1; j>=index; j--) {
1184                     ptT [j+1]  = ptT [j];
1185                     etaT[j+1]  = etaT[j];
1186                     phiT[j+1]  = phiT[j];
1187                 }
1188                 ptT [index] = fPtT [part];
1189                 etaT[index] = fEtaT[part];
1190                 phiT[index] = fPhiT[part];
1191                 iT++;
1192             } // particle associated
1193             if (iT > nT0) break;
1194         } // particle loop
1195         
1196         if (fBackground) {
1197             for (Int_t part = 0; part < fNtB; part++) {
1198                 if (fTrackListB[part] == nj+2) {
1199                     Int_t index = 0;
1200                     for (j=iT-1; j>=0; j--) {
1201                         if (fPtB[part] > ptT[j]) {
1202                             index = j+1;
1203
1204                             break;
1205                         }
1206                     }
1207                     for (j=iT-1; j>=index; j--) {
1208                         ptT [j+1]  = ptT [j];
1209                         etaT[j+1]  = etaT[j];
1210                         phiT[j+1]  = phiT[j];
1211                     }
1212                     ptT [index] = fPtB [part];
1213                     etaT[index] = fEtaB[part];
1214                     phiT[index] = fPhiB[part];
1215                     iT++;
1216                 } // particle associated
1217                 if (iT > nT0) break;
1218             } // particle loop
1219         } // Background available ?
1220
1221         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
1222         delete[] ptT;
1223         delete[] etaT;
1224         delete[] phiT;
1225     } // jet loop loop
1226 }
1227
1228 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1229 {
1230 // Propagates phi angle to EMCAL radius
1231 //
1232   static Float_t b = 0.0, rEMCAL = -1.0;
1233   if(rEMCAL<0) {
1234 // Get field in kGS
1235     b =  ((AliMagFCM*) gAlice->Field())->SolenoidField();
1236 // Get EMCAL radius in cm 
1237     rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1238     printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
1239   }
1240     Float_t dPhi = 0.;
1241 //
1242 //
1243 // bending radies
1244 // pt [Gev]
1245 // B  [kG]
1246 //
1247     Float_t rB = 3335.6 * pt / b;  // [cm]  (case of |charge|=1)
1248 //
1249 // check if particle is curling below EMCAL
1250     if (2.*rB < rEMCAL) {
1251         curls = kTRUE;
1252         return dPhi;
1253     }
1254 //
1255 // if not calculate delta phi
1256     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1257     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1258     dPhi = -TMath::Sign(dPhi, charge);
1259 //    
1260     return dPhi;
1261 }
1262
1263
1264 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1265 {
1266 // dummy for hbook calls
1267     ;
1268 }
1269
1270
1271
1272
1273
1274
1275
1276