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