1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //*-- Authors: Andreas Morsch (CERN)
20 //* Aleksei Pavlinov (WSU)
29 #include <TBranchElement.h>
31 #include <TClonesArray.h>
38 #include <TParticle.h>
39 #include <TParticlePDG.h>
40 #include <TPaveText.h>
41 #include <TPythia6Calls.h>
48 #include "AliEMCALJetFinder.h"
49 #include "AliHeader.h"
51 #include "AliMagFCM.h"
53 #include "AliGenerator.h"
54 #include "AliRunLoader.h"
56 #include "AliEMCALLoader.h"
57 #include "AliEMCALDigit.h"
58 #include "AliEMCALDigitizer.h"
59 #include "AliEMCALFast.h"
60 #include "AliEMCALGeometry.h"
61 #include "AliEMCALHadronCorrection.h"
62 #include "AliEMCALHit.h"
63 #include "AliEMCALJetMicroDst.h"
65 // Interface to FORTRAN
70 ClassImp(AliEMCALJetFinder)
72 //____________________________________________________________________________
73 AliEMCALJetFinder::AliEMCALJetFinder()
74 : fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE),
75 fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0),
76 fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0),
77 fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0),
78 fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0),
79 fDebug(0), fBackground(0), fConeRadius(0.3),
80 fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.),
81 fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0),
82 fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0),
83 fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.),
84 fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0),
85 fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0),
86 fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.),
87 fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0)
89 // Default constructor
108 fHadronCorrector = 0;
118 SetParametersForBgSubtraction();
121 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
122 : TTask(name, title),
123 fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE),
124 fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0),
125 fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0),
126 fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0),
127 fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0),
128 fDebug(0), fBackground(0), fConeRadius(0.3),
129 fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.),
130 fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0),
131 fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0),
132 fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.),
133 fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0),
134 fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0),
135 fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.),
136 fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0)
139 // Title is used in method GetFileNameForParameters();
141 fJets = new TClonesArray("AliEMCALJet",10000);
143 for (Int_t i = 0; i < 30000; i++)
165 fHadronCorrector = 0;
174 SetMomentumSmearing();
177 SetHadronCorrection();
181 SetParametersForBgSubtraction();
185 AliEMCALJetFinder::AliEMCALJetFinder(const AliEMCALJetFinder& jf)
186 : TTask(jf.GetName(), jf.GetTitle()),
187 fBGFileName(jf.fBGFileName), fEMCALWeight(jf.fEMCALWeight), fTrackWeight(jf.fTrackWeight),
188 fRandomBg(jf.fRandomBg),fWrite(jf.fWrite), fWeightingMethod(jf.fWeightingMethod), fJets(jf.fJets),
189 fLego(jf.fLego), fLegoB(jf.fLegoB), fhLegoTracks(jf.fhLegoTracks), fhLegoEMCAL(jf.fhLegoEMCAL),
190 fhLegoHadrCorr(jf.fhLegoHadrCorr), fhEff(jf.fhEff), fhCellEt(jf.fhCellEt),fhCellEMCALEt(jf.fhCellEMCALEt),
191 fhTrackPt(jf.fhTrackPt), fhTrackPtBcut(jf.fhTrackPtBcut), fhChPartMultInTpc(jf.fhChPartMultInTpc),
192 fhSinTheta(jf.fhSinTheta), fC1(jf.fC1), fHistsList(jf.fHistsList), fHadronCorrector(jf.fHadronCorrector),
193 fHCorrection(jf.fHCorrection),fDebug(jf.fDebug), fBackground(jf.fBackground), fConeRadius(jf.fConeRadius),
194 fPtCut(jf.fPtCut), fEtSeed(jf.fEtSeed), fMinJetEt(jf.fMinJetEt), fMinCellEt(jf.fMinCellEt),
195 fSamplingF(jf.fSamplingF), fSmear(jf.fSmear), fEffic(jf.fEffic), fK0N(jf.fK0N), fNjets(jf.fNjets),
196 fDeta(jf.fDeta), fDphi(jf.fDphi), fNcell(jf.fNcell), fNtot(jf.fNtot), fNbinEta(jf.fNbinEta),
197 fNbinPhi(jf.fNbinPhi), fEtaMin(jf.fEtaMin), fEtaMax(jf.fEtaMax), fPhiMin(jf.fPhiMin),
198 fPhiMax(jf.fPhiMax), fNt(jf.fNt), fNChTpc(jf.fNChTpc), fNtS(jf.fNtS), fTrackList(jf.fTrackList),
199 fPtT(jf.fPtT), fEtaT(jf.fEtaT), fPhiT(jf.fPhiT), fPdgT(jf.fPdgT), fNtB(jf.fNtB), fTrackListB(jf.fTrackListB),
200 fPtB(jf.fPtB), fEtaB(jf.fEtaB),fPhiB(jf.fPhiB), fPdgB(jf.fPdgB), fMode(jf.fMode), fMinMove(jf.fMinMove),
201 fMaxMove(jf.fMaxMove), fPrecBg(jf.fPrecBg), fError(jf.fError), fOutFileName(jf.fOutFileName),
202 fOutFile(jf.fOutFile), fInFile(jf.fInFile), fEvent(jf.fEvent)
207 void AliEMCALJetFinder::SetParametersForBgSubtraction
208 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
210 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
211 // at WSU Linux cluster - 11-feb-2002
212 // These parameters must be tuned more carefull !!!
219 //____________________________________________________________________________
220 AliEMCALJetFinder::~AliEMCALJetFinder()
232 delete fhLegoHadrCorr;
235 delete fhCellEMCALEt;
237 delete fhTrackPtBcut;
238 delete fhChPartMultInTpc;
246 delete[] fTrackListB;
254 # define jet_finder_ua1 jet_finder_ua1_
256 # define type_of_call
259 # define jet_finder_ua1 JET_FINDER_UA1
261 # define type_of_call _stdcall
264 extern "C" void type_of_call
265 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
266 Float_t etc[30000], Float_t etac[30000],
268 Float_t& min_move, Float_t& max_move, Int_t& mode,
269 Float_t& prec_bg, Int_t& ierror);
271 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
274 void AliEMCALJetFinder::Init()
277 // Geometry and I/O initialization
281 // Get geometry parameters from EMCAL
285 //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
286 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
288 // SetSamplingFraction(geom->GetSampling());
290 fNbinEta = geom->GetNZ();
291 fNbinPhi = geom->GetNPhi();
292 fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
293 fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
294 fEtaMin = geom->GetArm1EtaMin();
295 fEtaMax = geom->GetArm1EtaMax();
296 fDphi = (fPhiMax-fPhiMin)/fNbinPhi;
297 fDeta = (fEtaMax-fEtaMin)/fNbinEta;
298 fNtot = fNbinPhi*fNbinEta;
299 fWeightingMethod = kFALSE;
302 SetCellSize(fDeta, fDphi);
305 if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
310 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncelltot, Float_t etc[30000],
311 Float_t etac[30000], Float_t phic[30000],
312 Float_t minmove, Float_t maxmove, Int_t mode,
313 Float_t precbg, Int_t ierror)
315 // Wrapper for fortran coded jet finder
316 // Full argument list
317 jet_finder_ua1(ncell, ncelltot, etc, etac, phic,
318 minmove, maxmove, mode, precbg, ierror);
319 // Write result to output
320 if(fWrite) WriteJets();
324 void AliEMCALJetFinder::Find()
326 // Wrapper for fortran coded jet finder using member data for
329 Float_t minmove = fMinMove;
330 Float_t maxmove = fMaxMove;
332 Float_t precbg = fPrecBg;
334 ResetJets(); // 4-feb-2002 by PAI
336 jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
337 minmove, maxmove, mode, precbg, ierror);
339 // Write result to output
340 Int_t njet = Njets();
342 for (Int_t nj=0; nj<njet; nj++)
344 if (fWeightingMethod)
346 fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
352 fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
356 fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod);
357 fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) );
358 fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj) ));
362 FindTracksInJetCone();
363 if(fWrite) WriteJets();
368 Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
370 // Calculate the energy in the cone
371 Float_t newenergy = 0.0;
372 Float_t bineta,binphi;
373 TAxis *x = fhLegoEMCAL->GetXaxis();
374 TAxis *y = fhLegoEMCAL->GetYaxis();
375 for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
377 for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
379 bineta = x->GetBinCenter(i);
380 binphi = y->GetBinCenter(j);
381 if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
383 newenergy += fhLegoEMCAL->GetBinContent(i,j);
390 Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
392 // Calculate the track energy in the cone
393 Float_t newenergy = 0.0;
394 Float_t bineta,binphi;
395 TAxis *x = fhLegoTracks->GetXaxis();
396 TAxis *y = fhLegoTracks->GetYaxis();
397 for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
399 for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
401 bineta = x->GetBinCenter(i);
402 binphi = y->GetBinCenter(j);
403 if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
405 newenergy += fhLegoTracks->GetBinContent(i,j);
413 Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
415 // Calculate the weighted jet energy
417 Float_t newenergy = 0.0;
418 Float_t bineta,binphi;
419 TAxis *x = fhLegoEMCAL->GetXaxis();
420 TAxis *y = fhLegoEMCAL->GetYaxis();
423 for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
425 for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
427 bineta = x->GetBinCenter(i);
428 binphi = y->GetBinCenter(j);
429 if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
431 newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
441 Int_t AliEMCALJetFinder::Njets() const
443 // Get number of reconstructed jets
444 return EMCALJETS.njet;
447 Float_t AliEMCALJetFinder::JetEnergy(Int_t i) const
449 // Get reconstructed jet energy
450 return EMCALJETS.etj[i];
453 Float_t AliEMCALJetFinder::JetPhiL(Int_t i) const
455 // Get reconstructed jet phi from leading particle
456 return EMCALJETS.phij[0][i];
459 Float_t AliEMCALJetFinder::JetPhiW(Int_t i) const
461 // Get reconstructed jet phi from weighting
462 return EMCALJETS.phij[1][i];
465 Float_t AliEMCALJetFinder::JetEtaL(Int_t i) const
467 // Get reconstructed jet eta from leading particles
468 return EMCALJETS.etaj[0][i];
472 Float_t AliEMCALJetFinder::JetEtaW(Int_t i) const
474 // Get reconstructed jet eta from weighting
475 return EMCALJETS.etaj[1][i];
478 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
480 // Set grid cell size
481 EMCALCELLGEO.etaCellSize = eta;
482 EMCALCELLGEO.phiCellSize = phi;
485 void AliEMCALJetFinder::SetConeRadius(Float_t par)
487 // Set jet cone radius
488 EMCALJETPARAM.coneRad = par;
492 void AliEMCALJetFinder::SetEtSeed(Float_t par)
494 // Set et cut for seeds
495 EMCALJETPARAM.etSeed = par;
499 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
501 // Set minimum jet et
502 EMCALJETPARAM.ejMin = par;
506 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
508 // Set et cut per cell
509 EMCALJETPARAM.etMin = par;
513 void AliEMCALJetFinder::SetPtCut(Float_t par)
515 // Set pT cut on charged tracks
520 void AliEMCALJetFinder::Test()
523 // Test the finder call
525 const Int_t knmax = 30000;
527 Int_t ncelltot = 100;
539 Find(ncell, ncelltot, etc, etac, phic,
540 minmove, maxmove, mode, precbg, ierror);
548 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
553 TClonesArray &lrawcl = *fJets;
554 new(lrawcl[fNjets++]) AliEMCALJet(jet);
557 void AliEMCALJetFinder::ResetJets()
566 void AliEMCALJetFinder::WriteJets()
569 // Add all jets to the list
571 const Int_t kBufferSize = 4000;
572 const char* file = 0;
574 Int_t njet = Njets();
576 for (Int_t nj = 0; nj < njet; nj++)
583 AliRunLoader *rl = AliRunLoader::GetRunLoader();
584 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
588 // output written to input file
590 AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
591 TTree* pK = rl->TreeK();
592 file = (pK->GetCurrentFile())->GetName();
593 TBranch * jetBranch ;
595 printf("Make Branch - TreeR address %p %p\n",(void*)gAlice->TreeR(), (void*)pEMCAL);
596 //if (fJets && gAlice->TreeR()) {
597 if (fJets && emcalLoader->TreeR()) {
598 // pEMCAL->MakeBranchInTree(gAlice->TreeR(),
599 jetBranch = emcalLoader->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ;
600 //pEMCAL->MakeBranchInTree(gime->TreeR(),
606 //Int_t nev = gAlice->GetHeader()->GetEvent();
607 //gAlice->TreeR()->Fill();
610 //sprintf(hname,"TreeR%d", nev);
611 //gAlice->TreeR()->Write(hname);
612 //gAlice->TreeR()->Reset();
613 //gime->WriteRecPoints("OVERWRITE");
614 emcalLoader->WriteRecPoints("OVERWRITE");
618 // Output written to user specified output file
620 //TTree* pK = gAlice->TreeK();
621 TTree* pK = rl->TreeK();
622 fInFile = pK->GetCurrentFile();
626 sprintf(hname,"TreeR%d", fEvent);
627 TTree* treeJ = new TTree(hname, "EMCALJets");
628 treeJ->Branch("EMCALJets", &fJets, kBufferSize);
636 void AliEMCALJetFinder::BookLego()
639 // Book histo for discretization
643 // Don't add histos to the current directory
644 if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
646 // TH2::AddDirectory(0); // hists wil be put to the list from gROOT
647 // TH1::AddDirectory(0);
651 fLego = new TH2F("legoH","eta-phi",
652 fNbinEta, fEtaMin, fEtaMax,
653 fNbinPhi, fPhiMin, fPhiMax);
656 fLegoB = new TH2F("legoB","eta-phi for BG event",
657 fNbinEta, fEtaMin, fEtaMax,
658 fNbinPhi, fPhiMin, fPhiMax);
661 fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
662 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
664 fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
665 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
666 // Hadron correction map
667 fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
668 fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
669 // Hists. for tuning jet finder
670 fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
674 for(Int_t i=1; i<=1000; i++) eTmp[i] = 0.1*i; // step 100 mev
675 for(Int_t i=1001; i<=1100; i++) eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
677 fhCellEt = new TH1F("hCellEt","Cell E_{T} from fLego",
678 eTmp.GetSize()-1, eTmp.GetArray());
679 fhCellEMCALEt = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
680 eTmp.GetSize()-1, eTmp.GetArray());
681 fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
682 eTmp.GetSize()-1, eTmp.GetArray());
683 fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
684 eTmp.GetSize()-1, eTmp.GetArray());
686 fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
687 "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
689 fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
690 TAxis *xax = fhSinTheta->GetXaxis();
691 for(Int_t i=1; i<=fNbinEta; i++) {
692 Double_t eta = xax->GetBinCenter(i);
693 fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
696 //! first canvas for drawing
697 fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
700 void AliEMCALJetFinder::DumpLego()
703 // Dump lego histo into array
706 TAxis* xaxis = fLego->GetXaxis();
707 TAxis* yaxis = fLego->GetYaxis();
708 // fhCellEt->Clear();
710 for (Int_t i = 1; i <= fNbinEta; i++) {
711 for (Int_t j = 1; j <= fNbinPhi; j++) {
713 e = fLego->GetBinContent(i,j);
715 if (gRandom->Rndm() < 0.5) {
716 Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.));
720 if (e > 0.0) e -= fMinCellEt;
722 Float_t eta = xaxis->GetBinCenter(i);
723 Float_t phi = yaxis->GetBinCenter(j);
725 fEtaCell[fNcell] = eta;
726 fPhiCell[fNcell] = phi;
730 eH = fhLegoEMCAL->GetBinContent(i,j);
731 if(eH > 0.0) fhCellEMCALEt->Fill(eH);
739 void AliEMCALJetFinder::ResetMap()
742 // Reset eta-phi array
744 for (Int_t i=0; i<30000; i++)
753 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
757 const char* name = gAlice->Generator()->GetName();
758 enum {kPythia, kHijing, kHijingPara};
761 if (!strcmp(name, "Hijing")){
763 } else if (!strcmp(name, "Pythia")) {
765 } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) {
766 genType = kHijingPara;
769 printf("Fill tracks generated by %s %d\n", name, genType);
773 // Fill Cells with track information
776 printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
781 if (!fLego) BookLego();
783 if (flag == 0) fLego->Reset();
785 // Access particle information
786 Int_t npart = (gAlice->GetHeader())->GetNprimary();
787 Int_t ntr = (gAlice->GetHeader())->GetNtrack();
788 printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n",
789 npart, ntr, fLego->Integral());
794 // 1: selected for jet finding
797 if (fTrackList) delete[] fTrackList;
798 if (fPtT) delete[] fPtT;
799 if (fEtaT) delete[] fEtaT;
800 if (fPhiT) delete[] fPhiT;
801 if (fPdgT) delete[] fPdgT;
803 fTrackList = new Int_t [npart];
804 fPtT = new Float_t[npart];
805 fEtaT = new Float_t[npart];
806 fPhiT = new Float_t[npart];
807 fPdgT = new Int_t[npart];
811 Float_t chTmp=0.0; // charge of current particle
814 // this is for Pythia ??
815 for (Int_t part = 0; part < npart; part++) {
816 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
817 Int_t mpart = mPart->GetPdgCode();
818 Int_t child1 = mPart->GetFirstDaughter();
819 Float_t pT = mPart->Pt();
820 Float_t p = mPart->P();
821 Float_t phi = mPart->Phi();
823 if(pT > 0.001) eta = mPart->Eta();
824 Float_t theta = mPart->Theta();
825 if (fDebug>=2 && mPart->GetStatusCode()==1) {
826 printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
827 part, mpart, child1, mPart->GetLastDaughter(), mPart->GetStatusCode());
830 if (fDebug >= 2 && genType == kPythia) {
831 if (part == 6 || part == 7)
833 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
834 part-5, pT, eta, phi);
838 fTrackList[part] = 0;
839 fPtT[part] = pT; // must be change after correction for resolution !!!
844 // final state particles only
846 if (genType == kPythia) {
847 if (mPart->GetStatusCode() != 1) continue;
848 } else if (genType == kHijing) {
849 if (child1 >= 0 && child1 < npart) continue;
853 TParticlePDG* pdgP = 0;
854 // charged or neutral
855 pdgP = mPart->GetPDG();
856 chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
863 if (mpart != kNeutron &&
864 mpart != kNeutronBar &&
865 mpart != kK0Long) continue;
868 } else if (ich == 2) {
869 if (mpart == kNeutron ||
870 mpart == kNeutronBar ||
871 mpart == kK0Long) continue;
874 if (TMath::Abs(eta)<=0.9) fNChTpc++;
876 if (child1 >= 0 && child1 < npart) continue;
878 if (eta > fEtaMax || eta < fEtaMin) continue;
879 if (phi > fPhiMax || phi < fPhiMin ) continue;
882 printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
883 part, mpart, child1, eta, phi, pT, chTmp);
886 // Momentum smearing goes here ...
890 if (fSmear && TMath::Abs(chTmp)) {
891 pw = AliEMCALFast::SmearMomentum(1,p);
892 // p changed - take into account when calculate pT,
895 if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
899 // Tracking Efficiency and TPC acceptance goes here ...
901 if (fEffic && TMath::Abs(chTmp)) {
902 eff = 0.9; // AliEMCALFast::Efficiency(2,p);
903 if(fhEff) fhEff->Fill(p, eff);
904 if (AliEMCALFast::RandomReject(eff)) {
905 if(fDebug >= 5) printf(" reject due to unefficiency ");
910 // Correction of Hadronic Energy goes here ...
913 // phi propagation for hadronic correction
915 Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
916 Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
917 if(TMath::Abs(chTmp)) {
918 // hadr. correction only for charge particle
919 dphi = PropagatePhi(pT, chTmp, curls);
922 printf("\n Delta phi %f pT %f ", dphi, pT);
923 if (curls) printf("\n !! Track is curling");
925 if(!curls) fhTrackPtBcut->Fill(pT);
927 if (fHCorrection && !curls) {
928 if (!fHadronCorrector)
929 Fatal("AliEMCALJetFinder",
930 "Hadronic energy correction required but not defined !");
932 dpH = fHadronCorrector->GetEnergy(p, eta, 7);
933 eTdpH = dpH*TMath::Sin(theta);
935 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
936 phi, phiHC, -eTdpH); // correction is negative
938 xbin = fLego->GetXaxis()->FindBin(eta);
939 ybin = fLego->GetYaxis()->FindBin(phiHC);
940 cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
941 fLego->Fill(eta, phi, -fSamplingF*eTdpH );
942 cout <<"Hadron Correction affected bin - contents after correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
943 fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
947 // More to do ? Just think about it !
949 if (phi > fPhiMax || phi < fPhiMin ) continue;
951 if(TMath::Abs(chTmp) ) { // charge particle
952 if (pT > fPtCut && !curls) {
953 if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
954 eta , phi, pT, fNtS);
955 fLego->Fill(eta, phi, pT);
956 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
957 fTrackList[part] = 1;
960 } else if(ich > 0 || fK0N) {
961 // case of n, nbar and K0L
962 if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
963 eta , phi, pT, fNtS);
964 fLego->Fill(eta, phi, pT);
965 fTrackList[part] = 1;
970 for(Int_t i=0; i<fLego->GetSize(); i++) {
971 Float_t etc = (*fLego)[i];
972 if (etc > fMinCellEt) etsum += etc;
975 printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
976 printf(" Track selected(fNtS) %i \n", fNtS);
981 void AliEMCALJetFinder::FillFromHits(Int_t flag)
984 // Fill Cells with hit information
988 printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
992 if (!fLego) BookLego();
993 // Reset eta-phi maps if needed
994 if (flag == 0) { // default behavior
996 fhLegoTracks->Reset();
997 fhLegoEMCAL->Reset();
998 fhLegoHadrCorr->Reset();
1000 // Initialize from background event if available
1002 // Access hit information
1003 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1004 AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL");
1005 TTree *treeH = emcalLoader->TreeH();
1006 Int_t ntracks = (Int_t) treeH->GetEntries();
1011 // Double_t etH = 0.0;
1013 for (Int_t track=0; track<ntracks;track++) {
1014 gAlice->ResetHits();
1015 nbytes += treeH->GetEvent(track);
1019 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1021 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1023 Float_t x = mHit->X(); // x-pos of hit
1024 Float_t y = mHit->Y(); // y-pos
1025 Float_t z = mHit->Z(); // z-pos
1026 Float_t eloss = mHit->GetEnergy(); // deposited energy
1028 Float_t r = TMath::Sqrt(x*x+y*y);
1029 Float_t theta = TMath::ATan2(r,z);
1030 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
1031 Float_t phi = TMath::ATan2(y,x);
1033 if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
1035 // etH = fSamplingF*eloss*TMath::Sin(theta);
1036 fLego->Fill(eta, phi, eloss);
1040 // Transition from deposit energy to eT (eT = de*SF*sin(theta))
1042 for(Int_t i=1; i<=fLego->GetNbinsX(); i++){ // eta
1043 Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
1044 for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
1045 eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
1046 fLego->SetBinContent(i,j,eT);
1047 // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
1048 fhLegoEMCAL->SetBinContent(i,j,eT);
1049 if (eT > fMinCellEt) etsum += eT;
1053 // for(Int_t i=0; i<fLego->GetSize(); i++) {
1054 // (*fhLegoEMCAL)[i] = (*fLego)[i];
1055 // Float_t etc = (*fLego)[i];
1056 // if (etc > fMinCellEt) etsum += etc;
1059 printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
1063 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
1066 // Fill Cells with digit information
1071 if (!fLego) BookLego();
1072 if (flag == 0) fLego->Reset();
1079 TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
1080 TTree *treeD = gAlice->TreeD();
1081 TBranchElement* branchDg = (TBranchElement*)
1082 treeD->GetBranch("EMCAL");
1084 if (!branchDg) Fatal("AliEMCALJetFinder",
1085 "Reading digits requested but no digits in file !");
1087 branchDg->SetAddress(&digs);
1088 Int_t nent = (Int_t) branchDg->GetEntries();
1090 // Connect digitizer
1092 AliEMCALDigitizer* digr = new AliEMCALDigitizer();
1093 TBranchElement* branchDr = (TBranchElement*)
1094 treeD->GetBranch("AliEMCALDigitizer");
1095 branchDr->SetAddress(&digr);
1098 nbytes += branchDg->GetEntry(0);
1099 nbytes += branchDr->GetEntry(0);
1101 // Get digitizer parameters
1102 Float_t ecADCped = digr->GetECApedestal();
1103 Float_t ecADCcha = digr->GetECAchannel();
1105 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1106 AliEMCALGeometry* geom =
1107 AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
1110 Int_t ndig = digs->GetEntries();
1111 Info("FillFromDigits","Number of Digits: %d %d\n Parameters: EC: %f %f\n Geometry: %d %d",
1112 ndig, nent, ecADCped, ecADCcha, geom->GetNEta(), geom->GetNPhi());
1119 while ((sdg = (AliEMCALDigit*)(next())))
1121 Double_t pedestal = 0.;
1122 Double_t channel = 0.;
1123 if (geom->CheckAbsCellId(sdg->GetId())) { // May 31, 2006; IsInECA(Int_t) -> CheckAbsCellId
1124 pedestal = ecADCped;
1128 Fatal("FillFromDigits", "unexpected digit is number!") ;
1130 Float_t eta = sdg->GetEta();
1131 Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1132 Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1135 Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1136 eta, phi, amp, sdg->GetAmp());
1138 fLego->Fill(eta, phi, fSamplingF*amp);
1146 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1149 // Fill Cells with hit information
1154 if (!fLego) BookLego();
1156 if (flag == 0) fLego->Reset();
1158 // Flag charged tracks passing through TPC which
1159 // are associated to EMCAL Hits
1160 BuildTrackFlagTable();
1163 // Access particle information
1164 TTree *treeK = gAlice->TreeK();
1165 Int_t ntracks = (Int_t) treeK->GetEntries();
1167 if (fPtT) delete[] fPtT;
1168 if (fEtaT) delete[] fEtaT;
1169 if (fPhiT) delete[] fPhiT;
1170 if (fPdgT) delete[] fPdgT;
1172 fPtT = new Float_t[ntracks];
1173 fEtaT = new Float_t[ntracks];
1174 fPhiT = new Float_t[ntracks];
1175 fPdgT = new Int_t[ntracks];
1180 for (Int_t track = 0; track < ntracks; track++) {
1181 TParticle *mPart = gAlice->GetMCApp()->Particle(track);
1182 Float_t pT = mPart->Pt();
1183 Float_t phi = mPart->Phi();
1184 Float_t eta = mPart->Eta();
1186 if(fTrackList[track]) {
1190 fPdgT[track] = mPart->GetPdgCode();
1192 if (track < 2) continue; //Colliding particles?
1193 if (pT == 0 || pT < fPtCut) continue;
1195 fLego->Fill(eta, phi, pT);
1201 void AliEMCALJetFinder::FillFromParticles()
1203 // 26-feb-2002 PAI - for checking all chain
1204 // Work on particles level; accept all particle (not neutrino )
1206 Double_t pX=0, pY=0, pZ=0, energy=0; // checking conservation law
1210 if (!fLego) BookLego();
1213 // Access particles information
1214 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1215 if (fDebug >= 2 || npart<=0) {
1216 printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1217 if(npart<=0) return;
1221 RearrangeParticlesMemory(npart);
1223 // Go through the particles
1224 Int_t mpart, child1, child2, geantPdg;
1225 Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1227 for (Int_t part = 0; part < npart; part++) {
1229 fTrackList[part] = 0;
1231 mPart = gAlice->GetMCApp()->Particle(part);
1232 mpart = mPart->GetPdgCode();
1233 child1 = mPart->GetFirstDaughter();
1234 child2 = mPart->GetLastDaughter();
1242 e = mPart->Energy();
1244 // see pyedit in Pythia's text
1246 if (IsThisPartonsOrDiQuark(mpart)) continue;
1247 printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1248 part, mpart, geantPdg, px, py, pz, e, child1, child2);
1250 // exclude partons (21 - gluon, 92 - string)
1253 // exclude neutrinous also ??
1254 if (fDebug >= 11 && pT>0.01)
1255 printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1256 part, mpart, eta, phi, pT);
1261 fPdgT[part] = mpart;
1265 if (child1 >= 0 && child1 < npart) continue;
1267 // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1268 // part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1275 if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1277 if (eta > fEtaMax || eta < fEtaMin) continue;
1278 if (phi > fPhiMax || phi < fPhiMin ) continue;
1280 if(fK0N==0 ) { // exclude neutral hadrons
1281 if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1283 fTrackList[part] = 1;
1284 fLego->Fill(eta, phi, pT);
1287 printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n",
1288 pX, pY, pZ, energy);
1290 if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1293 void AliEMCALJetFinder::FillFromPartons()
1295 // function under construction - 13-feb-2002 PAI
1298 printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1302 if (!fLego) BookLego();
1305 // Access particle information
1306 Int_t npart = (gAlice->GetHeader())->GetNprimary();
1307 if (fDebug >= 2 || npart<=0)
1308 printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1309 fNt = 0; // for FindTracksInJetCone
1312 // Go through the partons
1314 for (Int_t part = 8; part < npart; part++) {
1315 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
1316 Int_t mpart = mPart->GetPdgCode();
1317 // Int_t child1 = MPart->GetFirstDaughter();
1318 Float_t pT = mPart->Pt();
1319 // Float_t p = MPart->P();
1320 Float_t phi = mPart->Phi();
1321 Float_t eta = mPart->Eta();
1322 // Float_t theta = MPart->Theta();
1323 statusCode = mPart->GetStatusCode();
1325 // accept partons (21 - gluon, 92 - string)
1326 if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1327 if (fDebug > 1 && pT>0.01)
1328 printf("\n part:%5d mpart %5d status %5d eta %8.2f phi %8.2f pT %8.2f ",
1329 part, mpart, statusCode, eta, phi, pT);
1330 // if (fDebug >= 3) MPart->Print();
1331 // accept partons before fragmentation - p.57 in Pythia manual
1332 // if(statusCode != 1) continue;
1334 if (eta > fEtaMax || eta < fEtaMin) continue;
1335 if (phi > fPhiMax || phi < fPhiMin ) continue;
1337 // if (child1 >= 0 && child1 < npart) continue;
1340 fLego->Fill(eta, phi, pT);
1346 void AliEMCALJetFinder::BuildTrackFlagTable() {
1348 // Method to generate a lookup table for TreeK particles
1349 // which are linked to hits in the EMCAL
1351 // --Author: J.L. Klay
1353 // Access hit information
1354 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1356 TTree *treeK = gAlice->TreeK(); // Get the number of entries in the kine tree
1357 Int_t nKTrks = (Int_t) treeK->GetEntries(); // (Number of particles created somewhere)
1359 if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
1360 fTrackList = new Int_t[nKTrks]; //before generating a new one
1362 for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
1366 AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL");
1367 //AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
1368 // TTree *treeH = gAlice->TreeH();
1369 TTree *treeH = emcalLoader->TreeH();
1370 Int_t ntracks = (Int_t) treeH->GetEntries();
1376 for (Int_t track=0; track<ntracks;track++) {
1377 gAlice->ResetHits();
1378 nbytes += treeH->GetEvent(track);
1384 for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1386 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1388 Int_t iTrk = mHit->Track(); // track number
1389 Int_t idprim = mHit->GetPrimary(); // primary particle
1391 //Determine the origin point of this particle - it made a hit in the EMCAL
1392 TParticle *trkPart = gAlice->GetMCApp()->Particle(iTrk);
1393 TParticlePDG *trkPdg = trkPart->GetPDG();
1394 Int_t trkCode = trkPart->GetPdgCode();
1396 if (trkCode < 10000) { //Big Ions cause problems for
1397 trkChg = trkPdg->Charge(); //this function. Since they aren't
1398 } else { //likely to make it very far, set
1399 trkChg = 0.0; //their charge to 0 for the Flag test
1401 Float_t vxTrk = trkPart->Vx();
1402 Float_t vyTrk = trkPart->Vy();
1403 Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1404 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1406 //Loop through the ancestry of the EMCAL entrance particles
1407 Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
1408 while (ancestor != -1) {
1409 TParticle *ancPart = gAlice->GetMCApp()->Particle(ancestor); //get particle info on ancestor
1410 TParticlePDG *ancPdg = ancPart->GetPDG();
1411 Int_t ancCode = ancPart->GetPdgCode();
1413 if (ancCode < 10000) {
1414 ancChg = ancPdg->Charge();
1418 Float_t vxAnc = ancPart->Vx();
1419 Float_t vyAnc = ancPart->Vy();
1420 Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1421 fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1422 ancestor = ancPart->GetFirstMother(); //Get the next ancestor
1425 //Determine the origin point of the primary particle associated with the hit
1426 TParticle *primPart = gAlice->GetMCApp()->Particle(idprim);
1427 TParticlePDG *primPdg = primPart->GetPDG();
1428 Int_t primCode = primPart->GetPdgCode();
1430 if (primCode < 10000) {
1431 primChg = primPdg->Charge();
1435 Float_t vxPrim = primPart->Vx();
1436 Float_t vyPrim = primPart->Vy();
1437 Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1438 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1444 Int_t AliEMCALJetFinder
1445 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1446 // Set the flag for the track
1451 if (charge == 0) neutral = 1;
1453 if (TMath::Abs(code) <= 6 ||
1455 code == 92) parton = 1;
1457 //It's not a parton, it's charged and it went through the TPC
1458 if (!parton && !neutral && radius <= 84.0) flag = 1;
1465 void AliEMCALJetFinder::SaveBackgroundEvent(const char *name)
1467 // Saves the eta-phi lego and the tracklist and name of file with BG events
1471 (*fLegoB) = (*fLegoB) + (*fLego);
1473 printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1474 fLegoB->Integral(), fLego->Integral());
1477 if (fPtB) delete[] fPtB;
1478 if (fEtaB) delete[] fEtaB;
1479 if (fPhiB) delete[] fPhiB;
1480 if (fPdgB) delete[] fPdgB;
1481 if (fTrackListB) delete[] fTrackListB;
1483 fPtB = new Float_t[fNtS];
1484 fEtaB = new Float_t[fNtS];
1485 fPhiB = new Float_t[fNtS];
1486 fPdgB = new Int_t [fNtS];
1487 fTrackListB = new Int_t [fNtS];
1491 for (Int_t i = 0; i < fNt; i++) {
1492 if (!fTrackList[i]) continue;
1493 fPtB [fNtB] = fPtT [i];
1494 fEtaB[fNtB] = fEtaT[i];
1495 fPhiB[fNtB] = fPhiT[i];
1496 fPdgB[fNtB] = fPdgT[i];
1497 fTrackListB[fNtB] = 1;
1501 printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1503 if(strlen(name) == 0) {
1504 TSeqCollection *li = gROOT->GetListOfFiles();
1506 for(Int_t i=0; i<li->GetSize(); i++) {
1507 nf = ((TFile*)li->At(i))->GetName();
1508 if(nf.Contains("backgorund")) break;
1514 printf("BG file name is \n %s\n", fBGFileName.Data());
1517 void AliEMCALJetFinder::InitFromBackground()
1521 if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1525 (*fLego) = (*fLego) + (*fLegoB);
1527 printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1528 fLego->Integral(), fLegoB->Integral());
1530 printf(" => fLego undefined \n");
1535 void AliEMCALJetFinder::FindTracksInJetCone()
1538 // Build list of tracks inside jet cone
1541 Int_t njet = Njets();
1542 for (Int_t nj = 0; nj < njet; nj++)
1544 Float_t etaj = JetEtaW(nj);
1545 Float_t phij = JetPhiW(nj);
1546 Int_t nT = 0; // number of associated tracks
1548 // Loop over particles in current event
1549 for (Int_t part = 0; part < fNt; part++) {
1550 if (!fTrackList[part]) continue;
1551 Float_t phi = fPhiT[part];
1552 Float_t eta = fEtaT[part];
1553 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1554 (phij-phi)*(phij-phi));
1555 if (dr < fConeRadius) {
1556 fTrackList[part] = nj+2;
1561 // Same for background event if available
1564 for (Int_t part = 0; part < fNtB; part++) {
1565 Float_t phi = fPhiB[part];
1566 Float_t eta = fEtaB[part];
1567 Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1568 (phij-phi)*(phij-phi));
1569 fTrackListB[part] = 1;
1571 if (dr < fConeRadius) {
1572 fTrackListB[part] = nj+2;
1576 } // Background available ?
1578 Int_t nT0 = nT + nTB;
1579 printf("Total number of tracks %d\n", nT0);
1581 if (nT0 > 1000) nT0 = 1000;
1583 Float_t* ptT = new Float_t[nT0];
1584 Float_t* etaT = new Float_t[nT0];
1585 Float_t* phiT = new Float_t[nT0];
1586 Int_t* pdgT = new Int_t[nT0];
1591 for (Int_t part = 0; part < fNt; part++) {
1592 if (fTrackList[part] == nj+2) {
1594 for (j=iT-1; j>=0; j--) {
1595 if (fPtT[part] > ptT[j]) {
1600 for (j=iT-1; j>=index; j--) {
1601 ptT [j+1] = ptT [j];
1602 etaT[j+1] = etaT[j];
1603 phiT[j+1] = phiT[j];
1604 pdgT[j+1] = pdgT[j];
1606 ptT [index] = fPtT [part];
1607 etaT[index] = fEtaT[part];
1608 phiT[index] = fPhiT[part];
1609 pdgT[index] = fPdgT[part];
1611 } // particle associated
1612 if (iT > nT0) break;
1616 for (Int_t part = 0; part < fNtB; part++) {
1617 if (fTrackListB[part] == nj+2) {
1619 for (j=iT-1; j>=0; j--) {
1620 if (fPtB[part] > ptT[j]) {
1626 for (j=iT-1; j>=index; j--) {
1627 ptT [j+1] = ptT [j];
1628 etaT[j+1] = etaT[j];
1629 phiT[j+1] = phiT[j];
1630 pdgT[j+1] = pdgT[j];
1632 ptT [index] = fPtB [part];
1633 etaT[index] = fEtaB[part];
1634 phiT[index] = fPhiB[part];
1635 pdgT[index] = fPdgB[part];
1637 } // particle associated
1638 if (iT > nT0) break;
1640 } // Background available ?
1642 fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1651 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1653 // Propagates phi angle to EMCAL radius
1655 static Float_t b = 0.0, rEMCAL = -1.0;
1657 b = gAlice->Field()->SolenoidField();
1658 // Get EMCAL radius in cm
1659 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1667 Float_t rB = 3335.6 * pt / TMath::Abs(b); // [cm] (case of |charge|=1)
1669 // check if particle is curling below EMCAL
1670 if (2.*rB < rEMCAL) {
1672 AliDebug(1, Form(" Low pt %f \n", pt));
1676 // if not calculate delta phi
1677 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1678 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1679 dPhi = -TMath::Sign(dPhi, charge);
1681 AliDebug(1, Form(" B %7.3f kG : rEMCAL %7.2f : dphi %7.5f(%7.5f)\n",
1682 b, rEMCAL, dPhi,dPhi*TMath::RadToDeg()));
1686 void hf1(Int_t& , Float_t& , Float_t& )
1688 // dummy for hbook calls
1692 void AliEMCALJetFinder::DrawLego(const char *opt)
1695 if(fLego) fLego->Draw(opt);
1698 void AliEMCALJetFinder::DrawLegoBackground(const char *opt)
1700 // Draw background lego
1701 if(fLegoB) fLegoB->Draw(opt);
1704 void AliEMCALJetFinder::DrawLegoEMCAL(const char *opt)
1707 if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);
1710 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1713 static TPaveText *varLabel=0;
1715 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1725 fhCellEMCALEt->Draw();
1730 fhTrackPtBcut->SetLineColor(2);
1731 fhTrackPtBcut->Draw("same");
1736 varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1737 varLabel->SetTextAlign(12);
1738 varLabel->SetFillColor(19); // see TAttFill
1739 TString tmp(GetTitle());
1740 varLabel->ReadFile(GetFileNameForParameters());
1744 if(mode) { // for saving picture to the file
1745 TString stmp(GetFileNameForParameters());
1746 stmp.ReplaceAll("_Par.txt",".ps");
1747 fC1->Print(stmp.Data());
1751 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1753 // Print all parameters out
1756 if(mode==0) file = stdout; // output to terminal
1758 file = fopen(GetFileNameForParameters(),"w");
1759 if(file==0) file = stdout;
1761 fprintf(file,"==== Filling lego ==== \n");
1762 fprintf(file,"Sampling fraction %6.3f ", fSamplingF);
1763 fprintf(file,"Smearing %6i ", fSmear);
1764 fprintf(file,"Efficiency %6i\n", fEffic);
1765 fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
1766 fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1767 fprintf(file,"==== Jet finding ==== \n");
1768 fprintf(file,"Cone radius %6.2f ", fConeRadius);
1769 fprintf(file,"Seed E_{T} %6.1f\n", fEtSeed);
1770 fprintf(file,"Min E_{T} of cell %6.1f ", fMinCellEt);
1771 fprintf(file,"Min E_{T} of jet %6.1f\n", fMinJetEt);
1773 fprintf(file,"==== Bg subtraction ==== \n");
1774 fprintf(file,"BG subtraction %6i ", fMode);
1775 fprintf(file,"Min cone move %6.2f\n", fMinMove);
1776 fprintf(file,"Max cone move %6.2f ", fMaxMove);
1777 fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1779 fprintf(file,"==== No Bg subtraction ==== \n");
1781 printf("BG file name is %s \n", fBGFileName.Data());
1782 if(file != stdout) fclose(file);
1785 void AliEMCALJetFinder::DrawLegos()
1789 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1793 gStyle->SetOptStat(111111);
1795 Int_t nent1, nent2, nent3, nent4;
1796 Double_t int1, int2, int3, int4;
1797 nent1 = (Int_t)fLego->GetEntries();
1798 int1 = fLego->Integral();
1800 if(int1) fLego->Draw("lego");
1802 nent2 = (Int_t)fhLegoTracks->GetEntries();
1803 int2 = fhLegoTracks->Integral();
1805 if(int2) fhLegoTracks->Draw("lego");
1807 nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1808 int3 = fhLegoEMCAL->Integral();
1810 if(int3) fhLegoEMCAL->Draw("lego");
1812 nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1813 int4 = fhLegoHadrCorr->Integral();
1815 if(int4) fhLegoHadrCorr->Draw("lego");
1817 // just for checking
1818 printf(" Integrals \n");
1819 printf("lego %10.3f \ntrack %10.3f \nhits %10.3f \nHCorr %10.3f\n-- %10.3f(must be 0)\n",
1820 int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1823 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(const char* dir)
1825 // Get paramters from a file
1827 if(strlen(dir)) tmp = dir;
1833 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1835 // See FillFromTracks() - npart must be positive
1836 if (fTrackList) delete[] fTrackList;
1837 if (fPtT) delete[] fPtT;
1838 if (fEtaT) delete[] fEtaT;
1839 if (fPhiT) delete[] fPhiT;
1840 if (fPdgT) delete[] fPdgT;
1843 fTrackList = new Int_t [npart];
1844 fPtT = new Float_t[npart];
1845 fEtaT = new Float_t[npart];
1846 fPhiT = new Float_t[npart];
1847 fPdgT = new Int_t[npart];
1849 printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1853 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1855 // Return quark info
1856 Int_t absPdg = TMath::Abs(pdg);
1857 if(absPdg<=6) return kTRUE; // quarks
1858 if(pdg == 21) return kTRUE; // gluon
1859 if(pdg == 92) return kTRUE; // string
1861 // see p.51 of Pythia Manual
1862 // Not include diquarks with c and b quark - 4-mar-2002
1863 // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1
1864 static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303};
1865 for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks
1870 void AliEMCALJetFinder::FindChargedJet()
1873 // Find jet from charged particle information only
1888 for (part = 0; part < fNt; part++) {
1889 if (!fTrackList[part]) continue;
1890 if (fPtT[part] > fEtSeed) nseed++;
1892 printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1893 Int_t* iSeeds = new Int_t[nseed];
1896 for (part = 0; part < fNt; part++) {
1897 if (!fTrackList[part]) continue;
1898 if (fPtT[part] > fEtSeed) iSeeds[nseed++] = part;
1909 // Find seed with highest pt
1911 Float_t ptmax = -1.;
1914 for (seed = 0; seed < nseed; seed++) {
1915 if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1920 if (ptmax < 0.) break;
1921 jndex = iSeeds[index];
1923 // Remove from the list
1925 printf("\n Next Seed %d %f", jndex, ptmax);
1927 // Find tracks in cone around seed
1929 Float_t phiSeed = fPhiT[jndex];
1930 Float_t etaSeed = fEtaT[jndex];
1936 for (part = 0; part < fNt; part++) {
1937 if (!fTrackList[part]) continue;
1938 Float_t deta = fEtaT[part] - etaSeed;
1939 Float_t dphi = fPhiT[part] - phiSeed;
1940 Float_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1941 if (dR < fConeRadius) {
1943 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1944 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1945 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1946 Float_t pz = fPtT[part] / TMath::Tan(theta);
1951 // if seed, remove it
1953 for (seed = 0; seed < nseed; seed++) {
1954 if (part == iSeeds[seed]) iSeeds[seed] = -1;
1960 // Estimate of jet direction
1961 Float_t phiJ = TMath::ATan2(pyJ, pxJ);
1962 Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1963 Float_t etaJ = TMath::Log(TMath::Tan(thetaJ / 2.));
1964 // Float_t ptJ = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1967 // Sum up all energy
1969 Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1970 Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1971 Int_t dIphi = Int_t(fConeRadius / fDphi);
1972 Int_t dIeta = Int_t(fConeRadius / fDeta);
1975 for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1976 for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1977 if (iPhi < 0 || iEta < 0) continue;
1978 Float_t dPhi = fPhiMin + iPhi * fDphi;
1979 Float_t dEta = fEtaMin + iEta * fDeta;
1980 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1981 sumE += fLego->GetBinContent(iEta, iPhi);
1987 fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1988 FindTracksInJetCone();
1989 printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1990 printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1991 printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1993 EMCALJETS.njet = njets;
1994 if (fWrite) WriteJets();
1997 // 16-jan-2003 - just for convenience
1998 void AliEMCALJetFinder::Browse(TBrowser* b)
2001 if(fHistsList) b->Add((TObject*)fHistsList);
2004 Bool_t AliEMCALJetFinder::IsFolder() const
2006 // Return folder status
2007 if(fHistsList) return kTRUE;
2011 const Char_t* AliEMCALJetFinder::GetNameOfVariant()
2013 // generate the literal string with info about jet finder
2015 sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
2016 fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
2018 nt.ReplaceAll(" ","");
2019 if(fBGFileName.Length()) {
2020 Int_t i1 = fBGFileName.Index("kBackground");
2021 Int_t i2 = fBGFileName.Index("/0000") - 1;
2022 if(i1>=0 && i2>=0) {
2023 TString bg(fBGFileName(i1,i2-i1+1));
2027 printf("<I> Name of variant %s \n", nt.Data());