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