]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALJetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removed files for new Flugg structure
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinder.cxx
CommitLineData
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$
b9d0a01d 18Revision 1.20.4.3 2002/10/10 15:07:49 hristov
19Updating VirtualMC to v3-09-02
20
21Revision 1.25 2002/09/13 10:24:57 morsch
22idem
23
9b91a677 24Revision 1.24 2002/09/13 10:21:13 morsch
25No cast to AliMagFCM.
26
a42a0f1a 27Revision 1.23 2002/06/27 09:24:26 morsch
28Uncomment the TH1::AddDirectory statement.
29
396840bb 30Revision 1.22 2002/05/22 13:48:43 morsch
31Pdg code added to track list.
32
975127ed 33Revision 1.21 2002/04/27 07:43:08 morsch
34Calculation of fDphi corrected (Renan Cabrera)
35
6fe407f0 36Revision 1.20 2002/03/12 01:06:23 pavlinov
37Testin output from generator
38
b561574c 39Revision 1.19 2002/02/27 00:46:33 pavlinov
40Added method FillFromParticles()
41
c44ecd4c 42Revision 1.18 2002/02/21 08:48:59 morsch
43Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
44
ff8fdf97 45Revision 1.17 2002/02/14 08:52:53 morsch
46Major updates by Aleksei Pavlinov:
47FillFromPartons, FillFromTracks, jetfinder configuration.
48
b7a73953 49Revision 1.16 2002/02/08 11:43:05 morsch
50SetOutputFileName(..) allows to specify an output file to which the
51reconstructed jets are written. If not set output goes to input file.
52Attention call Init() before processing.
53
08a3fb06 54Revision 1.15 2002/02/02 08:37:18 morsch
55Formula for rB corrected.
56
db5026bf 57Revision 1.14 2002/02/01 08:55:30 morsch
58Fill map with Et and pT.
59
f719d9b7 60Revision 1.13 2002/01/31 09:37:36 morsch
61Geometry parameters in constructor and call of SetCellSize()
62
be9787fe 63Revision 1.12 2002/01/23 13:40:23 morsch
64Fastidious debug print statement removed.
65
e4717047 66Revision 1.11 2002/01/22 17:25:47 morsch
67Some corrections for event mixing and bg event handling.
68
610f27c9 69Revision 1.10 2002/01/22 10:31:44 morsch
70Some correction for bg mixing.
71
46955412 72Revision 1.9 2002/01/21 16:28:42 morsch
73Correct sign of dphi.
74
f8f69f71 75Revision 1.8 2002/01/21 16:05:31 morsch
76- different phi-bin for hadron correction
77- provisions for background mixing.
78
ff114de3 79Revision 1.7 2002/01/21 15:47:18 morsch
80Bending radius correctly in cm.
81
b3feface 82Revision 1.6 2002/01/21 12:53:50 morsch
83authors
84
0c207dbc 85Revision 1.5 2002/01/21 12:47:47 morsch
86Possibility to include K0long and neutrons.
87
d53b7b0b 88Revision 1.4 2002/01/21 11:03:21 morsch
89Phi propagation introduced in FillFromTracks.
90
26dadf3a 91Revision 1.3 2002/01/18 05:07:56 morsch
92- hadronic correction
93- filling of digits
94- track selection upon EMCAL information
95
d594b22e 96*/
97
0c207dbc 98//*-- Authors: Andreas Morsch (CERN)
99//* J.L. Klay (LBL)
100//* Aleksei Pavlinov (WSU)
101
c44ecd4c 102#include <stdio.h>
d594b22e 103// From root ...
b561574c 104#include <TROOT.h>
471f69dc 105#include <TClonesArray.h>
106#include <TTree.h>
d594b22e 107#include <TBranchElement.h>
471f69dc 108#include <TFile.h>
c44ecd4c 109#include <TH1.h>
471f69dc 110#include <TH2.h>
c44ecd4c 111#include <TArrayF.h>
112#include <TCanvas.h>
b561574c 113#include <TList.h>
c44ecd4c 114#include <TPad.h>
115#include <TPaveText.h>
471f69dc 116#include <TAxis.h>
c44ecd4c 117#include <TStyle.h>
471f69dc 118#include <TParticle.h>
119#include <TParticlePDG.h>
b561574c 120#include <TPythia6Calls.h>
d594b22e 121
122// From AliRoot ...
471f69dc 123#include "AliEMCALJetFinder.h"
d594b22e 124#include "AliEMCALFast.h"
471f69dc 125#include "AliEMCALGeometry.h"
126#include "AliEMCALHit.h"
d594b22e 127#include "AliEMCALDigit.h"
128#include "AliEMCALDigitizer.h"
129#include "AliEMCALHadronCorrection.h"
b561574c 130#include "AliEMCALJetMicroDst.h"
471f69dc 131#include "AliRun.h"
26dadf3a 132#include "AliMagF.h"
133#include "AliMagFCM.h"
471f69dc 134#include "AliEMCAL.h"
135#include "AliHeader.h"
d53b7b0b 136#include "AliPDG.h"
975127ed 137#include "AliMC.h"
471f69dc 138
ff114de3 139// Interface to FORTRAN
140#include "Ecommon.h"
141
142
471f69dc 143ClassImp(AliEMCALJetFinder)
144
145//____________________________________________________________________________
146AliEMCALJetFinder::AliEMCALJetFinder()
147{
148// Default constructor
d594b22e 149 fJets = 0;
150 fNjets = 0;
151 fLego = 0;
ff114de3 152 fLegoB = 0;
b561574c 153
d594b22e 154 fTrackList = 0;
155 fPtT = 0;
156 fEtaT = 0;
157 fPhiT = 0;
975127ed 158 fPdgT = 0;
159
b561574c 160 fTrackListB = 0;
ff114de3 161 fPtB = 0;
162 fEtaB = 0;
163 fPhiB = 0;
975127ed 164 fPdgB = 0;
b561574c 165
b7a73953 166 fHCorrection = 0;
d594b22e 167 fHadronCorrector = 0;
b7a73953 168
c44ecd4c 169 fWrite = 0;
08a3fb06 170 fOutFileName = 0;
171 fOutFile = 0;
172 fInFile = 0;
173 fEvent = 0;
b7a73953 174
175 SetParametersForBgSubtraction();
471f69dc 176}
177
178AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
179 : TTask(name, title)
180{
181// Constructor
c44ecd4c 182// Title is used in method GetFileNameForParameters();
183//
471f69dc 184 fJets = new TClonesArray("AliEMCALJet",10000);
185 fNjets = 0;
d594b22e 186 for (Int_t i = 0; i < 30000; i++)
471f69dc 187 {
188 fEtCell[i] = 0.;
189 fEtaCell[i] = 0.;
190 fPhiCell[i] = 0.;
191 }
b561574c 192 fLego = 0;
193 fLegoB = 0;
194
ff114de3 195 fTrackList = 0;
196 fPtT = 0;
197 fEtaT = 0;
198 fPhiT = 0;
975127ed 199 fPdgT = 0;
b561574c 200
201 fTrackListB = 0;
ff114de3 202 fPtB = 0;
203 fEtaB = 0;
204 fPhiB = 0;
975127ed 205 fPdgB = 0;
b561574c 206
b7a73953 207 fHCorrection = 0;
208 fHadronCorrector = 0;
209 fBackground = 0;
c44ecd4c 210 fWrite = 0;
08a3fb06 211 fOutFileName = 0;
212 fOutFile = 0;
213 fInFile = 0;
214 fEvent = 0;
0c207dbc 215//
d594b22e 216 SetPtCut();
217 SetMomentumSmearing();
218 SetEfficiencySim();
219 SetDebug();
220 SetHadronCorrection();
221 SetSamplingFraction();
d53b7b0b 222 SetIncludeK0andN();
471f69dc 223
b7a73953 224 SetParametersForBgSubtraction();
225}
471f69dc 226
b7a73953 227void AliEMCALJetFinder::SetParametersForBgSubtraction
228(Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
229{
230// see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
231// at WSU Linux cluster - 11-feb-2002
232// These parameters must be tuned more carefull !!!
233 SetMode(mode);
234 SetMinMove(minMove);
235 SetMaxMove(maxMove);
236 SetPrecBg(precBg);
237}
471f69dc 238
239//____________________________________________________________________________
240AliEMCALJetFinder::~AliEMCALJetFinder()
241{
242// Destructor
243//
244 if (fJets){
245 fJets->Delete();
246 delete fJets;
247 }
248}
249
471f69dc 250#ifndef WIN32
251# define jet_finder_ua1 jet_finder_ua1_
252# define hf1 hf1_
253# define type_of_call
254
255#else
975127ed 256# define jet_finder_ua1 JET_FINDER_UA1
471f69dc 257# define hf1 HF1
258# define type_of_call _stdcall
259#endif
260
261extern "C" void type_of_call
262jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
263 Float_t etc[30000], Float_t etac[30000],
264 Float_t phic[30000],
265 Float_t& min_move, Float_t& max_move, Int_t& mode,
266 Float_t& prec_bg, Int_t& ierror);
267
268extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
975127ed 269
471f69dc 270
08a3fb06 271void AliEMCALJetFinder::Init()
272{
273//
274// Geometry and I/O initialization
275//
276//
277//
278// Get geometry parameters from EMCAL
279//
280//
281// Geometry
282 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
283 AliEMCALGeometry* geom =
284 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
285 fNbinEta = geom->GetNZ();
286 fNbinPhi = geom->GetNPhi();
287 fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
288 fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
289 fEtaMin = geom->GetArm1EtaMin();
290 fEtaMax = geom->GetArm1EtaMax();
6fe407f0 291 fDphi = (fPhiMax-fPhiMin)/fNbinPhi;
08a3fb06 292 fDeta = (fEtaMax-fEtaMin)/fNbinEta;
293 fNtot = fNbinPhi*fNbinEta;
294//
295 SetCellSize(fDeta, fDphi);
296//
297// I/O
298 if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
299}
471f69dc 300
301void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
302 Float_t etac[30000], Float_t phic[30000],
303 Float_t min_move, Float_t max_move, Int_t mode,
304 Float_t prec_bg, Int_t ierror)
305{
306// Wrapper for fortran coded jet finder
307// Full argument list
308 jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
309 min_move, max_move, mode, prec_bg, ierror);
310 // Write result to output
b7a73953 311 if(fWrite) WriteJets();
312 fEvent++;
471f69dc 313}
314
315void AliEMCALJetFinder::Find()
316{
317// Wrapper for fortran coded jet finder using member data for
318// argument list
319
b7a73953 320 Float_t min_move = fMinMove;
321 Float_t max_move = fMaxMove;
322 Int_t mode = fMode;
323 Float_t prec_bg = fPrecBg;
324 Int_t ierror;
325
326 ResetJets(); // 4-feb-2002 by PAI
471f69dc 327
328 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
329 min_move, max_move, mode, prec_bg, ierror);
b7a73953 330 fError = ierror;
471f69dc 331 // Write result to output
3bc109a9 332 Int_t njet = Njets();
d594b22e 333
3bc109a9 334 for (Int_t nj=0; nj<njet; nj++)
335 {
d594b22e 336
3bc109a9 337 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
338 JetPhiW(nj),
339 JetEtaW(nj));
340 }
b7a73953 341
c44ecd4c 342 FindTracksInJetCone();
b7a73953 343 if(fWrite) WriteJets();
08a3fb06 344 fEvent++;
471f69dc 345}
346
471f69dc 347Int_t AliEMCALJetFinder::Njets()
348{
349// Get number of reconstructed jets
350 return EMCALJETS.njet;
351}
352
353Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
354{
355// Get reconstructed jet energy
356 return EMCALJETS.etj[i];
357}
358
359Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
360{
361// Get reconstructed jet phi from leading particle
362 return EMCALJETS.phij[0][i];
363}
364
365Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
366{
367// Get reconstructed jet phi from weighting
368 return EMCALJETS.phij[1][i];
369}
370
371Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
372{
373// Get reconstructed jet eta from leading particles
374 return EMCALJETS.etaj[0][i];
375}
376
377
378Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
379{
380// Get reconstructed jet eta from weighting
381 return EMCALJETS.etaj[1][i];
382}
383
384void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
385{
386// Set grid cell size
387 EMCALCELLGEO.etaCellSize = eta;
388 EMCALCELLGEO.phiCellSize = phi;
389}
390
391void AliEMCALJetFinder::SetConeRadius(Float_t par)
392{
393// Set jet cone radius
394 EMCALJETPARAM.coneRad = par;
3bc109a9 395 fConeRadius = par;
471f69dc 396}
397
398void AliEMCALJetFinder::SetEtSeed(Float_t par)
399{
400// Set et cut for seeds
401 EMCALJETPARAM.etSeed = par;
d594b22e 402 fEtSeed = par;
471f69dc 403}
404
405void AliEMCALJetFinder::SetMinJetEt(Float_t par)
406{
407// Set minimum jet et
408 EMCALJETPARAM.ejMin = par;
d594b22e 409 fMinJetEt = par;
471f69dc 410}
411
412void AliEMCALJetFinder::SetMinCellEt(Float_t par)
413{
414// Set et cut per cell
415 EMCALJETPARAM.etMin = par;
d594b22e 416 fMinCellEt = par;
471f69dc 417}
418
3bc109a9 419void AliEMCALJetFinder::SetPtCut(Float_t par)
420{
421// Set pT cut on charged tracks
422 fPtCut = par;
423}
424
471f69dc 425
426void AliEMCALJetFinder::Test()
427{
428//
429// Test the finder call
430//
431 const Int_t nmax = 30000;
432 Int_t ncell = 10;
433 Int_t ncell_tot = 100;
434
435 Float_t etc[nmax];
436 Float_t etac[nmax];
437 Float_t phic[nmax];
438 Float_t min_move = 0;
439 Float_t max_move = 0;
440 Int_t mode = 0;
441 Float_t prec_bg = 0;
442 Int_t ierror = 0;
443
444
445 Find(ncell, ncell_tot, etc, etac, phic,
446 min_move, max_move, mode, prec_bg, ierror);
447
448}
449
450//
451// I/O
452//
453
454void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
455{
456 //
457 // Add a jet
458 //
459 TClonesArray &lrawcl = *fJets;
460 new(lrawcl[fNjets++]) AliEMCALJet(jet);
461}
462
463void AliEMCALJetFinder::ResetJets()
464{
465 //
466 // Reset Jet List
467 //
468 fJets->Clear();
469 fNjets = 0;
470}
471
472void AliEMCALJetFinder::WriteJets()
473{
474//
475// Add all jets to the list
476//
477 const Int_t kBufferSize = 4000;
08a3fb06 478 const char* file = 0;
479
471f69dc 480 Int_t njet = Njets();
08a3fb06 481
3bc109a9 482 for (Int_t nj = 0; nj < njet; nj++)
471f69dc 483 {
3bc109a9 484 AddJet(*fJetT[nj]);
485 delete fJetT[nj];
471f69dc 486 }
3bc109a9 487
08a3fb06 488// I/O
489 if (!fOutFileName) {
490//
491// output written to input file
492//
493 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
494 TTree* pK = gAlice->TreeK();
495 file = (pK->GetCurrentFile())->GetName();
496 if (fDebug > 1)
497 printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
498 if (fJets && gAlice->TreeR()) {
499 pEMCAL->MakeBranchInTree(gAlice->TreeR(),
500 "EMCALJets",
501 &fJets,
502 kBufferSize,
503 file);
504 }
505 Int_t nev = gAlice->GetHeader()->GetEvent();
506 gAlice->TreeR()->Fill();
507 char hname[30];
508 sprintf(hname,"TreeR%d", nev);
509 gAlice->TreeR()->Write(hname);
510 gAlice->TreeR()->Reset();
511 } else {
512//
513// Output written to user specified output file
514//
515 TTree* pK = gAlice->TreeK();
516 fInFile = pK->GetCurrentFile();
517
518 fOutFile->cd();
519 char hname[30];
520 sprintf(hname,"TreeR%d", fEvent);
521 TTree* treeJ = new TTree(hname, "EMCALJets");
522 treeJ->Branch("EMCALJets", &fJets, kBufferSize);
523 treeJ->Fill();
524 treeJ->Write(hname);
525 fInFile->cd();
526 }
471f69dc 527 ResetJets();
528}
529
530void AliEMCALJetFinder::BookLego()
531{
532//
533// Book histo for discretisation
534//
be9787fe 535
610f27c9 536//
537// Don't add histos to the current directory
b561574c 538 if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
539
396840bb 540 TH2::AddDirectory(0);
541 TH1::AddDirectory(0);
b561574c 542 gROOT->cd();
471f69dc 543//
46955412 544// Signal map
471f69dc 545 fLego = new TH2F("legoH","eta-phi",
be9787fe 546 fNbinEta, fEtaMin, fEtaMax,
547 fNbinPhi, fPhiMin, fPhiMax);
46955412 548//
549// Background map
b561574c 550 fLegoB = new TH2F("legoB","eta-phi for BG event",
be9787fe 551 fNbinEta, fEtaMin, fEtaMax,
552 fNbinPhi, fPhiMin, fPhiMax);
c44ecd4c 553
554// Tracks map
555 fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
556 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
557// EMCAL map
558 fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
559 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
560// Hadron correction map
561 fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
562 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
563// Hists. for tuning jet finder
564 fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
565
566 TArrayF eTmp(1101);
567 eTmp[0] = 0.0;
568 for(Int_t i=1; i<=1000; i++) eTmp[i] = 0.1*i; // step 100 mev
569 for(Int_t i=1001; i<=1100; i++) eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
570
571 fhCellEt = new TH1F("hCellEt","Cell E_{T} from fLego",
572 eTmp.GetSize()-1, eTmp.GetArray());
573 fhCellEMCALEt = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
574 eTmp.GetSize()-1, eTmp.GetArray());
575 fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
576 eTmp.GetSize()-1, eTmp.GetArray());
577 fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
578 eTmp.GetSize()-1, eTmp.GetArray());
b561574c 579
580 fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
581 "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
582
583 //! first canvas for drawing
584 fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kTRUE);
471f69dc 585}
586
587void AliEMCALJetFinder::DumpLego()
588{
589//
590// Dump lego histo into array
591//
592 fNcell = 0;
b7a73953 593 TAxis* Xaxis = fLego->GetXaxis();
594 TAxis* Yaxis = fLego->GetYaxis();
c44ecd4c 595 // fhCellEt->Clear();
596 Float_t e, eH;
597 for (Int_t i = 1; i <= fNbinEta; i++) {
598 for (Int_t j = 1; j <= fNbinPhi; j++) {
599 e = fLego->GetBinContent(i,j);
600 if (e > 0.0) {
601 Float_t eta = Xaxis->GetBinCenter(i);
602 Float_t phi = Yaxis->GetBinCenter(j);
603 fEtCell[fNcell] = e;
604 fEtaCell[fNcell] = eta;
605 fPhiCell[fNcell] = phi;
606 fNcell++;
607 fhCellEt->Fill(e);
608 }
609 if(fhLegoEMCAL) {
610 eH = fhLegoEMCAL->GetBinContent(i,j);
611 if(eH > 0.0) fhCellEMCALEt->Fill(eH);
612 }
471f69dc 613 } // phi
614 } // eta
615 fNcell--;
616}
617
618void AliEMCALJetFinder::ResetMap()
619{
620//
621// Reset eta-phi array
622
623 for (Int_t i=0; i<30000; i++)
624 {
625 fEtCell[i] = 0.;
626 fEtaCell[i] = 0.;
627 fPhiCell[i] = 0.;
628 }
629}
630
631
632void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
633{
634//
d594b22e 635// Fill Cells with track information
471f69dc 636//
b7a73953 637 if (fDebug >= 2)
638 printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
b561574c 639 fNChTpc = 0;
b7a73953 640
471f69dc 641 ResetMap();
642
643 if (!fLego) BookLego();
644// Reset
645 if (flag == 0) fLego->Reset();
646//
647// Access particle information
648 Int_t npart = (gAlice->GetHeader())->GetNprimary();
b561574c 649 Int_t ntr = (gAlice->GetHeader())->GetNtrack();
650 printf(" : #primary particles %i # tracks %i \n", npart, ntr);
b7a73953 651
3bc109a9 652// Create track list
653//
654// 0: not selected
655// 1: selected for jet finding
656// 2: inside jet
657// ....
658 if (fTrackList) delete[] fTrackList;
659 if (fPtT) delete[] fPtT;
660 if (fEtaT) delete[] fEtaT;
661 if (fPhiT) delete[] fPhiT;
975127ed 662 if (fPdgT) delete[] fPdgT;
3bc109a9 663
664 fTrackList = new Int_t [npart];
665 fPtT = new Float_t[npart];
666 fEtaT = new Float_t[npart];
667 fPhiT = new Float_t[npart];
975127ed 668 fPdgT = new Int_t[npart];
ff114de3 669
b561574c 670 fNt = npart;
ff114de3 671 fNtS = 0;
b7a73953 672 Float_t chTmp=0.0; // charge of current particle
b561574c 673 // Int_t idGeant;
ff114de3 674
b561574c 675 // this is for Pythia ??
c44ecd4c 676 for (Int_t part = 0; part < npart; part++) {
471f69dc 677 TParticle *MPart = gAlice->Particle(part);
678 Int_t mpart = MPart->GetPdgCode();
679 Int_t child1 = MPart->GetFirstDaughter();
680 Float_t pT = MPart->Pt();
d594b22e 681 Float_t p = MPart->P();
471f69dc 682 Float_t phi = MPart->Phi();
c44ecd4c 683 Float_t eta = -100.;
684 if(pT > 0.001) eta = MPart->Eta();
b7a73953 685 Float_t theta = MPart->Theta();
b561574c 686 if (fDebug>=2) {
687 printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
688 part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
689 }
46955412 690
b561574c 691 if (fDebug >= 2) {
d594b22e 692 if (part == 6 || part == 7)
693 {
610f27c9 694 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
d594b22e 695 part-5, pT, eta, phi);
696 }
697
698// if (mpart == 21)
699
700// printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
701// part, mpart, pT, eta, phi);
702 }
703
b7a73953 704 fTrackList[part] = 0;
705 fPtT[part] = pT; // must be change after correction for resolution !!!
3bc109a9 706 fEtaT[part] = eta;
707 fPhiT[part] = phi;
975127ed 708 fPdgT[part] = mpart;
709
3bc109a9 710
711 if (part < 2) continue;
b7a73953 712
713 // move to fLego->Fill because hadron correction must apply
714 // if particle hit to EMCAL - 11-feb-2002
715 // if (pT == 0 || pT < fPtCut) continue;
26dadf3a 716 TParticlePDG* pdgP = 0;
471f69dc 717// charged or neutral
b7a73953 718 pdgP = MPart->GetPDG();
719 chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
471f69dc 720 if (ich == 0) {
b7a73953 721 if (chTmp == 0) {
722 if (!fK0N) {
d53b7b0b 723 continue;
724 } else {
725 if (mpart != kNeutron &&
726 mpart != kNeutronBar &&
727 mpart != kK0Long) continue;
728 }
729 }
b7a73953 730 }
471f69dc 731// skip partons
732 if (TMath::Abs(mpart) <= 6 ||
733 mpart == 21 ||
734 mpart == 92) continue;
b561574c 735
736 if (TMath::Abs(eta)<=0.9) fNChTpc++;
737// final state only
738 if (child1 >= 0 && child1 < npart) continue;
471f69dc 739// acceptance cut
740 if (TMath::Abs(eta) > 0.7) continue;
b7a73953 741// Initial phi may be out of acceptance but track may
742// hit two the acceptance - see variable curls
b561574c 743 if (phi*180./TMath::Pi() > 120.) continue;
d594b22e 744//
b561574c 745 if (fDebug >= 3)
b7a73953 746 printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
747 part, mpart, child1, eta, phi, pT, chTmp);
26dadf3a 748//
26dadf3a 749//
d594b22e 750// Momentum smearing goes here ...
751//
c44ecd4c 752 fhTrackPt->Fill(pT);
b7a73953 753 Float_t pw;
754 if (fSmear && TMath::Abs(chTmp)) {
755 pw = AliEMCALFast::SmearMomentum(1,p);
756 // p changed - take into account when calculate pT,
757 // pz and so on .. ?
758 pT = (pw/p) * pT;
b561574c 759 if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
b7a73953 760 p = pw;
d594b22e 761 }
762//
763// Tracking Efficiency and TPC acceptance goes here ...
c44ecd4c 764 Float_t eff;
b7a73953 765 if (fEffic && TMath::Abs(chTmp)) {
c44ecd4c 766 // eff = AliEMCALFast::Efficiency(1,p);
767 eff = 0.95; // for testing 25-feb-2002
768 if(fhEff) fhEff->Fill(p, eff);
b7a73953 769 if (AliEMCALFast::RandomReject(eff)) {
b561574c 770 if(fDebug >= 5) printf(" reject due to unefficiency ");
b7a73953 771 continue;
772 }
d594b22e 773 }
774//
775// Correction of Hadronic Energy goes here ...
776//
b7a73953 777//
778// phi propagation for hadronic correction
779
780 Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
c44ecd4c 781 Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
b7a73953 782 if(TMath::Abs(chTmp)) {
783 // hadr. correction only for charge particle
784 dphi = PropagatePhi(pT, chTmp, curls);
785 phiHC = phi + dphi;
b561574c 786 if (fDebug >= 6) {
b7a73953 787 printf("\n Delta phi %f pT %f ", dphi, pT);
788 if (curls) printf("\n !! Track is curling");
789 }
c44ecd4c 790 if(!curls) fhTrackPtBcut->Fill(pT);
ff114de3 791
b7a73953 792 if (fHCorrection && !curls) {
793 if (!fHadronCorrector)
794 Fatal("AliEMCALJetFinder",
d594b22e 795 "Hadronic energy correction required but not defined !");
b7a73953 796
c44ecd4c 797 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
798 eTdpH = dpH*TMath::Sin(theta);
b7a73953 799
b561574c 800 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
c44ecd4c 801 phi, phiHC, -eTdpH); // correction is negative
802 fLego->Fill(eta, phiHC, -eTdpH);
803 fhLegoHadrCorr->Fill(eta, phiHC, eTdpH);
b7a73953 804 }
805 }
d594b22e 806//
807// More to do ? Just think about it !
808//
b7a73953 809
810 if (phi*180./TMath::Pi() > 120.) continue;
811
812 if(TMath::Abs(chTmp) ) { // charge particle
813 if (pT > fPtCut && !curls) {
b561574c 814 if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
b7a73953 815 eta , phi, pT);
816 fLego->Fill(eta, phi, pT);
c44ecd4c 817 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
b7a73953 818 fTrackList[part] = 1;
819 fNtS++;
820 }
c44ecd4c 821 } else if(ich==0 && fK0N) {
b7a73953 822 // case of n, nbar and K0L
b561574c 823 if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
b7a73953 824 eta , phi, pT);
825 fLego->Fill(eta, phi, pT);
826 fTrackList[part] = 1;
827 fNtS++;
828 }
829
c44ecd4c 830 } // primary loop
471f69dc 831 DumpLego();
832}
833
834void AliEMCALJetFinder::FillFromHits(Int_t flag)
835{
836//
d594b22e 837// Fill Cells with hit information
471f69dc 838//
839//
b7a73953 840 if (fDebug >= 2)
841 printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
842
471f69dc 843 ResetMap();
844
845 if (!fLego) BookLego();
c44ecd4c 846// Reset eta-phi maps if needed
847 if (flag == 0) { // default behavior
848 fLego->Reset();
849 fhLegoTracks->Reset();
850 fhLegoEMCAL->Reset();
851 fhLegoHadrCorr->Reset();
852 }
ff114de3 853// Initialize from background event if available
471f69dc 854//
855// Access hit information
856 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
471f69dc 857
858 TTree *treeH = gAlice->TreeH();
859 Int_t ntracks = (Int_t) treeH->GetEntries();
860//
861// Loop over tracks
862//
863 Int_t nbytes = 0;
c44ecd4c 864 Double_t etH = 0.0;
471f69dc 865
471f69dc 866 for (Int_t track=0; track<ntracks;track++) {
867 gAlice->ResetHits();
868 nbytes += treeH->GetEvent(track);
869//
870// Loop over hits
871//
872 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
873 mHit;
874 mHit=(AliEMCALHit*) pEMCAL->NextHit())
875 {
876 Float_t x = mHit->X(); // x-pos of hit
877 Float_t y = mHit->Y(); // y-pos
878 Float_t z = mHit->Z(); // z-pos
879 Float_t eloss = mHit->GetEnergy(); // deposited energy
d594b22e 880
881 Float_t r = TMath::Sqrt(x*x+y*y);
882 Float_t theta = TMath::ATan2(r,z);
883 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
884 Float_t phi = TMath::ATan2(y,x);
610f27c9 885
b7a73953 886 if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
610f27c9 887
c44ecd4c 888 etH = fSamplingF*eloss*TMath::Sin(theta);
889 fLego->Fill(eta, phi, etH);
890 // fhLegoEMCAL->Fill(eta, phi, etH);
471f69dc 891 } // Hit Loop
892 } // Track Loop
c44ecd4c 893 // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
894 for(Int_t i=0; i<fLego->GetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i];
895 // DumpLego(); ??
471f69dc 896}
897
d594b22e 898void AliEMCALJetFinder::FillFromDigits(Int_t flag)
899{
900//
901// Fill Cells with digit information
902//
903//
904 ResetMap();
905
906 if (!fLego) BookLego();
907 if (flag == 0) fLego->Reset();
908 Int_t nbytes;
909
910
911//
912// Connect digits
913//
914 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
915 TTree *treeD = gAlice->TreeD();
916 TBranchElement* branchDg = (TBranchElement*)
917 treeD->GetBranch("EMCAL");
918
919 if (!branchDg) Fatal("AliEMCALJetFinder",
920 "Reading digits requested but no digits in file !");
921
922 branchDg->SetAddress(&digs);
923 Int_t nent = (Int_t) branchDg->GetEntries();
924//
925// Connect digitizer
926//
927 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
928 TBranchElement* branchDr = (TBranchElement*)
929 treeD->GetBranch("AliEMCALDigitizer");
930 branchDr->SetAddress(&digr);
931//
932//
933 nbytes += branchDg->GetEntry(0);
934 nbytes += branchDr->GetEntry(0);
935//
936// Get digitizer parameters
937 Float_t towerADCped = digr->GetTowerpedestal();
938 Float_t towerADCcha = digr->GetTowerchannel();
939 Float_t preshoADCped = digr->GetPreShopedestal();
940 Float_t preshoADCcha = digr->GetPreShochannel();
941
942 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
943 AliEMCALGeometry* geom =
944 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
945
946 if (fDebug) {
947 Int_t ndig = digs->GetEntries();
948 printf("\n Number of Digits: %d %d\n", ndig, nent);
949 printf("\n Parameters: %f %f %f %f\n",
950 towerADCped, towerADCcha, preshoADCped, preshoADCcha );
951 printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
952 }
953
954//
955// Loop over digits
956 AliEMCALDigit* sdg;
957 TIter next(digs);
958 while ((sdg = (AliEMCALDigit*)(next())))
959 {
960 Double_t pedestal;
961 Double_t channel;
962 if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
963 {
964 pedestal = preshoADCped;
965 channel = preshoADCcha;
966 } else {
967 pedestal = towerADCped;
968 channel = towerADCcha;
969 }
970
971 Float_t eta = sdg->GetEta();
972 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
973 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
974
975 if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
976 eta, phi, amp, sdg->GetAmp());
977
978 fLego->Fill(eta, phi, fSamplingF*amp);
979 } // digit loop
980//
981// Dump lego hist
982 DumpLego();
983}
984
985
986void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
987{
988//
989// Fill Cells with hit information
990//
991//
992 ResetMap();
993
994 if (!fLego) BookLego();
995// Reset
996 if (flag == 0) fLego->Reset();
997
ff114de3 998// Flag charged tracks passing through TPC which
999// are associated to EMCAL Hits
d594b22e 1000 BuildTrackFlagTable();
1001
1002//
1003// Access particle information
ff8fdf97 1004 TTree *treeK = gAlice->TreeK();
1005 Int_t ntracks = (Int_t) treeK->GetEntries();
d594b22e 1006
ff114de3 1007 if (fPtT) delete[] fPtT;
1008 if (fEtaT) delete[] fEtaT;
1009 if (fPhiT) delete[] fPhiT;
975127ed 1010 if (fPdgT) delete[] fPdgT;
ff114de3 1011
d594b22e 1012 fPtT = new Float_t[ntracks];
1013 fEtaT = new Float_t[ntracks];
1014 fPhiT = new Float_t[ntracks];
975127ed 1015 fPdgT = new Int_t[ntracks];
d594b22e 1016
ff114de3 1017 fNt = ntracks;
1018 fNtS = 0;
1019
d594b22e 1020 for (Int_t track = 0; track < ntracks; track++) {
1021 TParticle *MPart = gAlice->Particle(track);
1022 Float_t pT = MPart->Pt();
1023 Float_t phi = MPart->Phi();
1024 Float_t eta = MPart->Eta();
1025
1026 if(fTrackList[track]) {
1027 fPtT[track] = pT;
1028 fEtaT[track] = eta;
1029 fPhiT[track] = phi;
975127ed 1030 fPdgT[track] = MPart->GetPdgCode();
1031
d594b22e 1032 if (track < 2) continue; //Colliding particles?
1033 if (pT == 0 || pT < fPtCut) continue;
ff114de3 1034 fNtS++;
d594b22e 1035 fLego->Fill(eta, phi, pT);
1036 }
1037 } // track loop
1038 DumpLego();
1039}
1040
c44ecd4c 1041void AliEMCALJetFinder::FillFromParticles()
1042{
1043// 26-feb-2002 PAI - for checking all chain
1044// Work on particles level; accept all particle (not neutrino )
1045
b561574c 1046 Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law
1047 fNChTpc = 0;
c44ecd4c 1048
1049 ResetMap();
1050 if (!fLego) BookLego();
1051 fLego->Reset();
1052//
1053// Access particles information
1054 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1055 if (fDebug >= 2 || npart<=0) {
b561574c 1056 printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1057 if(npart<=0) return;
c44ecd4c 1058 }
1059 fNt = npart;
1060 fNtS = 0;
1061 RearrangeParticlesMemory(npart);
1062
1063// Go through the particles
b561574c 1064 Int_t mpart, child1, child2, geantPdg;
1065 Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
c44ecd4c 1066 TParticle *MPart=0;
1067 for (Int_t part = 0; part < npart; part++) {
1068
1069 fTrackList[part] = 0;
c44ecd4c 1070
1071 MPart = gAlice->Particle(part);
1072 mpart = MPart->GetPdgCode();
1073 child1 = MPart->GetFirstDaughter();
b561574c 1074 child2 = MPart->GetLastDaughter();
c44ecd4c 1075 pT = MPart->Pt();
1076 phi = MPart->Phi();
1077 eta = MPart->Eta();
b561574c 1078
1079 px = MPart->Px();
1080 py = MPart->Py();
1081 pz = MPart->Pz();
1082 e = MPart->Energy();
1083
1084// see pyedit in Pythia's text
975127ed 1085 geantPdg = mpart;
1086// Int_t kc = pycomp_(&mpart);
1087// TString name = GetPythiaParticleName(mpart);
b561574c 1088 // printf(" mpart %6.6i;kc %6.6i -> gid %3.3i",mpart,kc,geantPdg);
1089 //printf(" (%s)\n", name.Data());
975127ed 1090 if (IsThisPartonsOrDiQuark(mpart)) continue;
1091 printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1092 part, mpart, geantPdg, px, py, pz, e, child1, child2);
c44ecd4c 1093
1094// exclude partons (21 - gluon, 92 - string)
975127ed 1095
1096
c44ecd4c 1097// exclude neutrinous also ??
b561574c 1098 if (fDebug >= 11 && pT>0.01)
1099 printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
c44ecd4c 1100 part, mpart, eta, phi, pT);
1101
1102 fPtT[part] = pT;
1103 fEtaT[part] = eta;
1104 fPhiT[part] = phi;
975127ed 1105 fPdgT[part] = mpart;
1106
b561574c 1107// final state only
1108 if (child1 >= 0 && child1 < npart) continue;
1109
1110 // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1111 // part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1112 PX += px;
1113 PY += py;
1114 PZ += pz;
1115 E += e;
1116
1117
1118 if (TMath::Abs(eta) <= 0.9) fNChTpc++;
c44ecd4c 1119// acceptance cut
1120 if (TMath::Abs(eta) > 0.7) continue;
1121 if (phi*180./TMath::Pi() > 120.) continue;
c44ecd4c 1122//
1123 if(fK0N==0 ) { // exclude neutral hadrons
1124 if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1125 }
1126 fTrackList[part] = 1;
1127 fLego->Fill(eta, phi, pT);
1128
1129 } // primary loop
b561574c 1130 printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n",
1131 PX, PY, PZ, E);
c44ecd4c 1132 DumpLego();
b561574c 1133 if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
c44ecd4c 1134}
1135
b7a73953 1136void AliEMCALJetFinder::FillFromPartons()
1137{
1138// function under construction - 13-feb-2002 PAI
1139
1140 if (fDebug >= 2)
1141 printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1142 //
1143
1144 ResetMap();
1145 if (!fLego) BookLego();
1146 fLego->Reset();
1147//
1148// Access particle information
1149 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1150 if (fDebug >= 2 || npart<=0)
1151 printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1152 fNt = 0; // for FindTracksInJetCone
1153 fNtS = 0;
1154
1155// Go through the partons
1156 Int_t statusCode=0;
1157 for (Int_t part = 8; part < npart; part++) {
1158 TParticle *MPart = gAlice->Particle(part);
1159 Int_t mpart = MPart->GetPdgCode();
1160 // Int_t child1 = MPart->GetFirstDaughter();
1161 Float_t pT = MPart->Pt();
1162 // Float_t p = MPart->P();
1163 Float_t phi = MPart->Phi();
1164 Float_t eta = MPart->Eta();
1165 // Float_t theta = MPart->Theta();
1166 statusCode = MPart->GetStatusCode();
1167
1168// accept partons (21 - gluon, 92 - string)
1169 if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1170 if (fDebug > 1 && pT>0.01)
1171 printf("\n part:%5d mpart %5d status %5d eta %8.2f phi %8.2f pT %8.2f ",
1172 part, mpart, statusCode, eta, phi, pT);
1173 // if (fDebug >= 3) MPart->Print();
1174// accept partons before fragmentation - p.57 in Pythia manual
1175// if(statusCode != 1) continue;
1176// acceptance cut
1177 if (TMath::Abs(eta) > 0.7) continue;
1178 if (phi*180./TMath::Pi() > 120.) continue;
1179// final state only
1180// if (child1 >= 0 && child1 < npart) continue;
1181//
1182//
1183 fLego->Fill(eta, phi, pT);
1184
1185 } // primary loop
1186 DumpLego();
1187}
1188
d594b22e 1189void AliEMCALJetFinder::BuildTrackFlagTable() {
1190
1191// Method to generate a lookup table for TreeK particles
1192// which are linked to hits in the EMCAL
1193//
1194// --Author: J.L. Klay
1195//
d594b22e 1196// Access hit information
1197 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1198
1199 TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
1200 Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
1201
1202 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
1203 fTrackList = new Int_t[nKTrks]; //before generating a new one
1204
1205 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
1206 fTrackList[i] = 0;
1207 }
1208
1209 TTree *treeH = gAlice->TreeH();
1210 Int_t ntracks = (Int_t) treeH->GetEntries();
1211//
1212// Loop over tracks
1213//
1214 Int_t nbytes = 0;
1215
1216 for (Int_t track=0; track<ntracks;track++) {
1217 gAlice->ResetHits();
1218 nbytes += treeH->GetEvent(track);
1219 if (pEMCAL) {
1220
1221//
1222// Loop over hits
1223//
1224 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1225 mHit;
1226 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1227 {
1228 Int_t iTrk = mHit->Track(); // track number
1229 Int_t idprim = mHit->GetPrimary(); // primary particle
1230
1231 //Determine the origin point of this particle - it made a hit in the EMCAL
1232 TParticle *trkPart = gAlice->Particle(iTrk);
1233 TParticlePDG *trkPdg = trkPart->GetPDG();
1234 Int_t trkCode = trkPart->GetPdgCode();
1235 Double_t trkChg;
1236 if (trkCode < 10000) { //Big Ions cause problems for
1237 trkChg = trkPdg->Charge(); //this function. Since they aren't
1238 } else { //likely to make it very far, set
1239 trkChg = 0.0; //their charge to 0 for the Flag test
1240 }
1241 Float_t vxTrk = trkPart->Vx();
1242 Float_t vyTrk = trkPart->Vy();
1243 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1244 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1245
1246 //Loop through the ancestry of the EMCAL entrance particles
1247 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
1248 while (ancestor != -1) {
1249 TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
1250 TParticlePDG *ancPdg = ancPart->GetPDG();
1251 Int_t ancCode = ancPart->GetPdgCode();
1252 Double_t ancChg;
1253 if (ancCode < 10000) {
1254 ancChg = ancPdg->Charge();
1255 } else {
1256 ancChg = 0.0;
1257 }
1258 Float_t vxAnc = ancPart->Vx();
1259 Float_t vyAnc = ancPart->Vy();
1260 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1261 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1262 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
1263 }
1264
1265 //Determine the origin point of the primary particle associated with the hit
1266 TParticle *primPart = gAlice->Particle(idprim);
1267 TParticlePDG *primPdg = primPart->GetPDG();
1268 Int_t primCode = primPart->GetPdgCode();
1269 Double_t primChg;
1270 if (primCode < 10000) {
1271 primChg = primPdg->Charge();
1272 } else {
1273 primChg = 0.0;
1274 }
1275 Float_t vxPrim = primPart->Vx();
1276 Float_t vyPrim = primPart->Vy();
1277 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1278 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1279 } // Hit Loop
1280 } //if (pEMCAL)
1281 } // Track Loop
1282}
1283
1284Int_t AliEMCALJetFinder
1285::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1286
ff114de3 1287 Int_t flag = 0;
1288 Int_t parton = 0;
d594b22e 1289 Int_t neutral = 0;
1290
1291 if (charge == 0) neutral = 1;
1292
1293 if (TMath::Abs(code) <= 6 ||
1294 code == 21 ||
1295 code == 92) parton = 1;
1296
1297 //It's not a parton, it's charged and it went through the TPC
1298 if (!parton && !neutral && radius <= 84.0) flag = 1;
1299
1300 return flag;
1301}
1302
ff114de3 1303
1304
1305void AliEMCALJetFinder::SaveBackgroundEvent()
1306{
1307// Saves the eta-phi lego and the tracklist
1308//
b561574c 1309 if (fLegoB) {
1310 fLegoB->Reset();
1311 (*fLegoB) = (*fLegoB) + (*fLego);
1312 if(fDebug)
1313 printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1314 fLegoB->Integral(), fLego->Integral());
1315 }
ff114de3 1316
1317 if (fPtB) delete[] fPtB;
1318 if (fEtaB) delete[] fEtaB;
1319 if (fPhiB) delete[] fPhiB;
975127ed 1320 if (fPdgB) delete[] fPdgB;
ff114de3 1321 if (fTrackListB) delete[] fTrackListB;
1322
1323 fPtB = new Float_t[fNtS];
1324 fEtaB = new Float_t[fNtS];
1325 fPhiB = new Float_t[fNtS];
975127ed 1326 fPdgB = new Int_t [fNtS];
ff114de3 1327 fTrackListB = new Int_t [fNtS];
46955412 1328
1329 fNtB = 0;
1330
ff114de3 1331 for (Int_t i = 0; i < fNt; i++) {
1332 if (!fTrackList[i]) continue;
610f27c9 1333 fPtB [fNtB] = fPtT [i];
1334 fEtaB[fNtB] = fEtaT[i];
1335 fPhiB[fNtB] = fPhiT[i];
975127ed 1336 fPdgB[fNtB] = fPdgT[i];
1337
46955412 1338 fTrackListB[fNtB] = 1;
ff114de3 1339 fNtB++;
1340 }
610f27c9 1341 fBackground = 1;
b561574c 1342 printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
ff114de3 1343}
1344
1345void AliEMCALJetFinder::InitFromBackground()
1346{
1347//
1348//
b561574c 1349 if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
46955412 1350
b561574c 1351 if (fLego) {
396840bb 1352 fLego->Reset();
1353 (*fLego) = (*fLego) + (*fLegoB);
1354 if(fDebug)
1355 printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1356 fLego->Integral(), fLegoB->Integral());
b561574c 1357 } else {
396840bb 1358 printf(" => fLego undefined \n");
b561574c 1359 }
ff114de3 1360}
1361
ff114de3 1362
3bc109a9 1363void AliEMCALJetFinder::FindTracksInJetCone()
1364{
1365//
1366// Build list of tracks inside jet cone
1367//
1368// Loop over jets
1369 Int_t njet = Njets();
1370 for (Int_t nj = 0; nj < njet; nj++)
1371 {
1372 Float_t etaj = JetEtaW(nj);
1373 Float_t phij = JetPhiW(nj);
1374 Int_t nT = 0; // number of associated tracks
1375
ff114de3 1376// Loop over particles in current event
3bc109a9 1377 for (Int_t part = 0; part < fNt; part++) {
1378 if (!fTrackList[part]) continue;
1379 Float_t phi = fPhiT[part];
1380 Float_t eta = fEtaT[part];
1381 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1382 (phij-phi)*(phij-phi));
1383 if (dr < fConeRadius) {
1384 fTrackList[part] = nj+2;
1385 nT++;
1386 } // < ConeRadius ?
1387 } // particle loop
d594b22e 1388
ff114de3 1389// Same for background event if available
1390 Int_t nTB = 0;
1391 if (fBackground) {
1392 for (Int_t part = 0; part < fNtB; part++) {
1393 Float_t phi = fPhiB[part];
1394 Float_t eta = fEtaB[part];
1395 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1396 (phij-phi)*(phij-phi));
610f27c9 1397 fTrackListB[part] = 1;
1398
ff114de3 1399 if (dr < fConeRadius) {
610f27c9 1400 fTrackListB[part] = nj+2;
ff114de3 1401 nTB++;
1402 } // < ConeRadius ?
1403 } // particle loop
1404 } // Background available ?
1405
1406 Int_t nT0 = nT + nTB;
1407
1408 if (nT0 > 50) nT0 = 50;
d594b22e 1409
610f27c9 1410 Float_t* ptT = new Float_t[nT0];
1411 Float_t* etaT = new Float_t[nT0];
1412 Float_t* phiT = new Float_t[nT0];
975127ed 1413 Int_t* pdgT = new Int_t[nT0];
1414
d594b22e 1415 Int_t iT = 0;
1416 Int_t j;
1417
3bc109a9 1418 for (Int_t part = 0; part < fNt; part++) {
1419 if (fTrackList[part] == nj+2) {
d594b22e 1420 Int_t index = 0;
1421 for (j=iT-1; j>=0; j--) {
1422 if (fPtT[part] > ptT[j]) {
1423 index = j+1;
1424 break;
1425 }
1426 }
1427 for (j=iT-1; j>=index; j--) {
1428 ptT [j+1] = ptT [j];
1429 etaT[j+1] = etaT[j];
1430 phiT[j+1] = phiT[j];
975127ed 1431 pdgT[j+1] = pdgT[j];
d594b22e 1432 }
1433 ptT [index] = fPtT [part];
1434 etaT[index] = fEtaT[part];
1435 phiT[index] = fPhiT[part];
975127ed 1436 pdgT[index] = fPdgT[part];
ff114de3 1437 iT++;
3bc109a9 1438 } // particle associated
ff114de3 1439 if (iT > nT0) break;
3bc109a9 1440 } // particle loop
d594b22e 1441
ff114de3 1442 if (fBackground) {
1443 for (Int_t part = 0; part < fNtB; part++) {
610f27c9 1444 if (fTrackListB[part] == nj+2) {
ff114de3 1445 Int_t index = 0;
1446 for (j=iT-1; j>=0; j--) {
1447 if (fPtB[part] > ptT[j]) {
1448 index = j+1;
610f27c9 1449
ff114de3 1450 break;
1451 }
1452 }
1453 for (j=iT-1; j>=index; j--) {
1454 ptT [j+1] = ptT [j];
1455 etaT[j+1] = etaT[j];
1456 phiT[j+1] = phiT[j];
975127ed 1457 pdgT[j+1] = pdgT[j];
ff114de3 1458 }
1459 ptT [index] = fPtB [part];
1460 etaT[index] = fEtaB[part];
1461 phiT[index] = fPhiB[part];
975127ed 1462 pdgT[index] = fPdgB[part];
ff114de3 1463 iT++;
1464 } // particle associated
1465 if (iT > nT0) break;
1466 } // particle loop
1467 } // Background available ?
1468
975127ed 1469 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
3bc109a9 1470 delete[] ptT;
1471 delete[] etaT;
1472 delete[] phiT;
975127ed 1473 delete[] pdgT;
1474
3bc109a9 1475 } // jet loop loop
1476}
1477
26dadf3a 1478Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1479{
1480// Propagates phi angle to EMCAL radius
1481//
b7a73953 1482 static Float_t b = 0.0, rEMCAL = -1.0;
1483 if(rEMCAL<0) {
1484// Get field in kGS
9b91a677 1485 b = gAlice->Field()->SolenoidField();
b7a73953 1486// Get EMCAL radius in cm
1487 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1488 printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
1489 }
26dadf3a 1490 Float_t dPhi = 0.;
26dadf3a 1491//
1492//
1493// bending radies
db5026bf 1494// pt [Gev]
1495// B [kG]
1496//
b7a73953 1497 Float_t rB = 3335.6 * pt / b; // [cm] (case of |charge|=1)
26dadf3a 1498//
1499// check if particle is curling below EMCAL
1500 if (2.*rB < rEMCAL) {
1501 curls = kTRUE;
1502 return dPhi;
1503 }
1504//
1505// if not calculate delta phi
1506 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1507 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
f8f69f71 1508 dPhi = -TMath::Sign(dPhi, charge);
26dadf3a 1509//
1510 return dPhi;
26dadf3a 1511}
1512
471f69dc 1513void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1514{
1515// dummy for hbook calls
1516 ;
1517}
1518
c44ecd4c 1519void AliEMCALJetFinder::DrawLego(Char_t *opt)
1520{fLego->Draw(opt);}
471f69dc 1521
c44ecd4c 1522void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
1523{fhLegoEMCAL->Draw(opt);}
d594b22e 1524
c44ecd4c 1525void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1526{
1527 static TPaveText *varLabel=0;
1528 if(!fC1) {
1529 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1530 }
1531 fC1->Clear();
1532 fC1->Divide(2,2);
1533 fC1->cd(1);
1534 gPad->SetLogy(1);
1535 fhCellEt->Draw();
1536
1537 fC1->cd(2);
1538 gPad->SetLogy(1);
1539 fhCellEMCALEt->Draw();
1540
1541 fC1->cd(3);
1542 gPad->SetLogy(1);
1543 fhTrackPt->Draw();
1544 fhTrackPtBcut->SetLineColor(2);
1545 fhTrackPtBcut->Draw("same");
1546
1547 fC1->cd(4);
1548 if(!varLabel) {
1549 PrintParameters(1);
1550 varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1551 varLabel->SetTextAlign(12);
1552 varLabel->SetFillColor(19); // see TAttFill
1553 TString tmp(GetTitle());
1554 varLabel->ReadFile(GetFileNameForParameters());
1555 }
1556 varLabel->Draw();
1557 fC1->Update();
1558 if(mode) { // for saving picture to the file
1559 TString stmp(GetFileNameForParameters());
1560 stmp.ReplaceAll("_Par.txt",".ps");
1561 fC1->Print(stmp.Data());
1562 }
1563}
d594b22e 1564
c44ecd4c 1565void AliEMCALJetFinder::PrintParameters(Int_t mode)
1566{
1567 FILE *file=0;
1568 if(mode==0) file = stdout; // output to terminal
1569 else {
1570 file = fopen(GetFileNameForParameters(),"w");
1571 if(file==0) file = stdout;
1572 }
1573 fprintf(file,"==== Filling lego ==== \n");
1574 fprintf(file,"Smearing %6i ", fSmear);
1575 fprintf(file,"Efficiency %6i\n", fEffic);
1576 fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
1577 fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1578 fprintf(file,"==== Jet finding ==== \n");
1579 fprintf(file,"Cone radius %6.2f ", fConeRadius);
1580 fprintf(file,"Seed E_{T} %6.1f\n", fEtSeed);
1581 fprintf(file,"Min E_{T} of cell %6.1f ", fMinCellEt);
1582 fprintf(file,"Min E_{T} of jet %6.1f\n", fMinJetEt);
1583 if(fMode) {
1584 fprintf(file,"==== Bg subtraction ==== \n");
1585 fprintf(file,"BG subtraction %6i ", fMode);
1586 fprintf(file,"Min cone move %6.2f\n", fMinMove);
1587 fprintf(file,"Max cone move %6.2f ", fMaxMove);
1588 fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1589 } else
1590 fprintf(file,"==== No Bg subtraction ==== \n");
1591 if(file != stdout) fclose(file);
1592}
d594b22e 1593
c44ecd4c 1594void AliEMCALJetFinder::DrawLegos()
1595{
1596 if(!fC1) {
1597 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1598 }
1599 fC1->Clear();
1600 fC1->Divide(2,2);
1601 gStyle->SetOptStat(111111);
1602
1603 Int_t nent1, nent2, nent3, nent4;
1604 Double_t int1, int2, int3, int4;
1605 nent1 = (Int_t)fLego->GetEntries();
1606 int1 = fLego->Integral();
1607 fC1->cd(1);
1608 if(int1) fLego->Draw("lego");
1609
1610 nent2 = (Int_t)fhLegoTracks->GetEntries();
1611 int2 = fhLegoTracks->Integral();
1612 fC1->cd(2);
1613 if(int2) fhLegoTracks->Draw("lego");
1614
1615 nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1616 int3 = fhLegoEMCAL->Integral();
1617 fC1->cd(3);
1618 if(int3) fhLegoEMCAL->Draw("lego");
1619
1620 nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1621 int4 = fhLegoHadrCorr->Integral();
1622 fC1->cd(4);
1623 if(int4) fhLegoHadrCorr->Draw("lego");
1624
1625 // just for checking
1626 printf(" Integrals \n");
1627 printf("lego %10.3f \ntrack %10.3f \nhits %10.3f \nHCorr %10.3f\n-- %10.3f(must be 0)\n",
1628 int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1629}
d594b22e 1630
c44ecd4c 1631const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1632{
1633 static TString tmp;
1634 if(strlen(dir)) tmp = dir;
1635 tmp += GetTitle();
1636 tmp += "_Par.txt";
1637 return tmp.Data();
1638}
d594b22e 1639
c44ecd4c 1640void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1641{ // See FillFromTracks() - npart must be positive
1642 if (fTrackList) delete[] fTrackList;
1643 if (fPtT) delete[] fPtT;
1644 if (fEtaT) delete[] fEtaT;
1645 if (fPhiT) delete[] fPhiT;
975127ed 1646 if (fPdgT) delete[] fPdgT;
c44ecd4c 1647
1648 if(npart>0) {
1649 fTrackList = new Int_t [npart];
1650 fPtT = new Float_t[npart];
1651 fEtaT = new Float_t[npart];
1652 fPhiT = new Float_t[npart];
975127ed 1653 fPdgT = new Int_t[npart];
c44ecd4c 1654 } else {
1655 printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1656 }
1657}
b561574c 1658
1659Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1660{
1661 Int_t absPdg = TMath::Abs(pdg);
1662 if(absPdg<=6) return kTRUE; // quarks
1663 if(pdg == 21) return kTRUE; // gluon
1664 if(pdg == 92) return kTRUE; // string
1665
1666 // see p.51 of Pythia Manual
1667 // Not include diquarks with c and b quark - 4-mar-2002
1668 // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1
1669 static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303};
1670 for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks
1671
1672 return kFALSE;
1673}
1674
1675TString &AliEMCALJetFinder::GetPythiaParticleName(Int_t kf)
1676{// see subroutine PYNAME in PYTHIA
1677 static TString sname;
1678 char name[16];
1679 pyname_(&kf, name, 16);
1680 for(Int_t i=0; i<16; i++){
1681 if(name[i] == ' ') {
1682 name[i] = '\0';
1683 break;
1684 }
1685 }
1686 sname = name;
1687 return sname;
1688}