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