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