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