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 / b; // [cm] (case of |charge|=1)
1669 // check if particle is curling below EMCAL
1670 if (2.*rB < rEMCAL) {
1675 // if not calculate delta phi
1676 Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1677 dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1678 dPhi = -TMath::Sign(dPhi, charge);
1683 void hf1(Int_t& , Float_t& , Float_t& )
1685 // dummy for hbook calls
1689 void AliEMCALJetFinder::DrawLego(const char *opt)
1692 if(fLego) fLego->Draw(opt);
1695 void AliEMCALJetFinder::DrawLegoBackground(const char *opt)
1697 // Draw background lego
1698 if(fLegoB) fLegoB->Draw(opt);
1701 void AliEMCALJetFinder::DrawLegoEMCAL(const char *opt)
1704 if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);
1707 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1710 static TPaveText *varLabel=0;
1712 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1722 fhCellEMCALEt->Draw();
1727 fhTrackPtBcut->SetLineColor(2);
1728 fhTrackPtBcut->Draw("same");
1733 varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1734 varLabel->SetTextAlign(12);
1735 varLabel->SetFillColor(19); // see TAttFill
1736 TString tmp(GetTitle());
1737 varLabel->ReadFile(GetFileNameForParameters());
1741 if(mode) { // for saving picture to the file
1742 TString stmp(GetFileNameForParameters());
1743 stmp.ReplaceAll("_Par.txt",".ps");
1744 fC1->Print(stmp.Data());
1748 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1750 // Print all parameters out
1753 if(mode==0) file = stdout; // output to terminal
1755 file = fopen(GetFileNameForParameters(),"w");
1756 if(file==0) file = stdout;
1758 fprintf(file,"==== Filling lego ==== \n");
1759 fprintf(file,"Sampling fraction %6.3f ", fSamplingF);
1760 fprintf(file,"Smearing %6i ", fSmear);
1761 fprintf(file,"Efficiency %6i\n", fEffic);
1762 fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
1763 fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1764 fprintf(file,"==== Jet finding ==== \n");
1765 fprintf(file,"Cone radius %6.2f ", fConeRadius);
1766 fprintf(file,"Seed E_{T} %6.1f\n", fEtSeed);
1767 fprintf(file,"Min E_{T} of cell %6.1f ", fMinCellEt);
1768 fprintf(file,"Min E_{T} of jet %6.1f\n", fMinJetEt);
1770 fprintf(file,"==== Bg subtraction ==== \n");
1771 fprintf(file,"BG subtraction %6i ", fMode);
1772 fprintf(file,"Min cone move %6.2f\n", fMinMove);
1773 fprintf(file,"Max cone move %6.2f ", fMaxMove);
1774 fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1776 fprintf(file,"==== No Bg subtraction ==== \n");
1778 printf("BG file name is %s \n", fBGFileName.Data());
1779 if(file != stdout) fclose(file);
1782 void AliEMCALJetFinder::DrawLegos()
1786 fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1790 gStyle->SetOptStat(111111);
1792 Int_t nent1, nent2, nent3, nent4;
1793 Double_t int1, int2, int3, int4;
1794 nent1 = (Int_t)fLego->GetEntries();
1795 int1 = fLego->Integral();
1797 if(int1) fLego->Draw("lego");
1799 nent2 = (Int_t)fhLegoTracks->GetEntries();
1800 int2 = fhLegoTracks->Integral();
1802 if(int2) fhLegoTracks->Draw("lego");
1804 nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1805 int3 = fhLegoEMCAL->Integral();
1807 if(int3) fhLegoEMCAL->Draw("lego");
1809 nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1810 int4 = fhLegoHadrCorr->Integral();
1812 if(int4) fhLegoHadrCorr->Draw("lego");
1814 // just for checking
1815 printf(" Integrals \n");
1816 printf("lego %10.3f \ntrack %10.3f \nhits %10.3f \nHCorr %10.3f\n-- %10.3f(must be 0)\n",
1817 int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1820 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(const char* dir)
1822 // Get paramters from a file
1824 if(strlen(dir)) tmp = dir;
1830 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1832 // See FillFromTracks() - npart must be positive
1833 if (fTrackList) delete[] fTrackList;
1834 if (fPtT) delete[] fPtT;
1835 if (fEtaT) delete[] fEtaT;
1836 if (fPhiT) delete[] fPhiT;
1837 if (fPdgT) delete[] fPdgT;
1840 fTrackList = new Int_t [npart];
1841 fPtT = new Float_t[npart];
1842 fEtaT = new Float_t[npart];
1843 fPhiT = new Float_t[npart];
1844 fPdgT = new Int_t[npart];
1846 printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1850 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1852 // Return quark info
1853 Int_t absPdg = TMath::Abs(pdg);
1854 if(absPdg<=6) return kTRUE; // quarks
1855 if(pdg == 21) return kTRUE; // gluon
1856 if(pdg == 92) return kTRUE; // string
1858 // see p.51 of Pythia Manual
1859 // Not include diquarks with c and b quark - 4-mar-2002
1860 // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1
1861 static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303};
1862 for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks
1867 void AliEMCALJetFinder::FindChargedJet()
1870 // Find jet from charged particle information only
1885 for (part = 0; part < fNt; part++) {
1886 if (!fTrackList[part]) continue;
1887 if (fPtT[part] > fEtSeed) nseed++;
1889 printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1890 Int_t* iSeeds = new Int_t[nseed];
1893 for (part = 0; part < fNt; part++) {
1894 if (!fTrackList[part]) continue;
1895 if (fPtT[part] > fEtSeed) iSeeds[nseed++] = part;
1906 // Find seed with highest pt
1908 Float_t ptmax = -1.;
1911 for (seed = 0; seed < nseed; seed++) {
1912 if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1917 if (ptmax < 0.) break;
1918 jndex = iSeeds[index];
1920 // Remove from the list
1922 printf("\n Next Seed %d %f", jndex, ptmax);
1924 // Find tracks in cone around seed
1926 Float_t phiSeed = fPhiT[jndex];
1927 Float_t etaSeed = fEtaT[jndex];
1933 for (part = 0; part < fNt; part++) {
1934 if (!fTrackList[part]) continue;
1935 Float_t deta = fEtaT[part] - etaSeed;
1936 Float_t dphi = fPhiT[part] - phiSeed;
1937 Float_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1938 if (dR < fConeRadius) {
1940 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1941 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1942 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1943 Float_t pz = fPtT[part] / TMath::Tan(theta);
1948 // if seed, remove it
1950 for (seed = 0; seed < nseed; seed++) {
1951 if (part == iSeeds[seed]) iSeeds[seed] = -1;
1957 // Estimate of jet direction
1958 Float_t phiJ = TMath::ATan2(pyJ, pxJ);
1959 Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1960 Float_t etaJ = TMath::Log(TMath::Tan(thetaJ / 2.));
1961 // Float_t ptJ = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1964 // Sum up all energy
1966 Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1967 Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1968 Int_t dIphi = Int_t(fConeRadius / fDphi);
1969 Int_t dIeta = Int_t(fConeRadius / fDeta);
1972 for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1973 for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1974 if (iPhi < 0 || iEta < 0) continue;
1975 Float_t dPhi = fPhiMin + iPhi * fDphi;
1976 Float_t dEta = fEtaMin + iEta * fDeta;
1977 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1978 sumE += fLego->GetBinContent(iEta, iPhi);
1984 fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1985 FindTracksInJetCone();
1986 printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1987 printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1988 printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1990 EMCALJETS.njet = njets;
1991 if (fWrite) WriteJets();
1994 // 16-jan-2003 - just for convenience
1995 void AliEMCALJetFinder::Browse(TBrowser* b)
1998 if(fHistsList) b->Add((TObject*)fHistsList);
2001 Bool_t AliEMCALJetFinder::IsFolder() const
2003 // Return folder status
2004 if(fHistsList) return kTRUE;
2008 const Char_t* AliEMCALJetFinder::GetNameOfVariant()
2010 // generate the literal string with info about jet finder
2012 sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
2013 fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
2015 nt.ReplaceAll(" ","");
2016 if(fBGFileName.Length()) {
2017 Int_t i1 = fBGFileName.Index("kBackground");
2018 Int_t i2 = fBGFileName.Index("/0000") - 1;
2019 if(i1>=0 && i2>=0) {
2020 TString bg(fBGFileName(i1,i2-i1+1));
2024 printf("<I> Name of variant %s \n", nt.Data());