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