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