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