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