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