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