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