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