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