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