* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.18 2002/02/21 08:48:59 morsch
-Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
-
-Revision 1.17 2002/02/14 08:52:53 morsch
-Major updates by Aleksei Pavlinov:
-FillFromPartons, FillFromTracks, jetfinder configuration.
-
-Revision 1.16 2002/02/08 11:43:05 morsch
-SetOutputFileName(..) allows to specify an output file to which the
-reconstructed jets are written. If not set output goes to input file.
-Attention call Init() before processing.
-
-Revision 1.15 2002/02/02 08:37:18 morsch
-Formula for rB corrected.
-
-Revision 1.14 2002/02/01 08:55:30 morsch
-Fill map with Et and pT.
-
-Revision 1.13 2002/01/31 09:37:36 morsch
-Geometry parameters in constructor and call of SetCellSize()
-
-Revision 1.12 2002/01/23 13:40:23 morsch
-Fastidious debug print statement removed.
-
-Revision 1.11 2002/01/22 17:25:47 morsch
-Some corrections for event mixing and bg event handling.
-
-Revision 1.10 2002/01/22 10:31:44 morsch
-Some correction for bg mixing.
-
-Revision 1.9 2002/01/21 16:28:42 morsch
-Correct sign of dphi.
-
-Revision 1.8 2002/01/21 16:05:31 morsch
-- different phi-bin for hadron correction
-- provisions for background mixing.
-
-Revision 1.7 2002/01/21 15:47:18 morsch
-Bending radius correctly in cm.
-
-Revision 1.6 2002/01/21 12:53:50 morsch
-authors
-
-Revision 1.5 2002/01/21 12:47:47 morsch
-Possibility to include K0long and neutrons.
-
-Revision 1.4 2002/01/21 11:03:21 morsch
-Phi propagation introduced in FillFromTracks.
-
-Revision 1.3 2002/01/18 05:07:56 morsch
-- hadronic correction
-- filling of digits
-- track selection upon EMCAL information
-
-*/
+/* $Id$ */
//*-- Authors: Andreas Morsch (CERN)
//* J.L. Klay (LBL)
#include <stdio.h>
// From root ...
-#include <TClonesArray.h>
-#include <TTree.h>
+#include <TArrayF.h>
+#include <TAxis.h>
#include <TBranchElement.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
-#include <TArrayF.h>
-#include <TCanvas.h>
+#include <TList.h>
+#include <TPDGCode.h>
#include <TPad.h>
-#include <TPaveText.h>
-#include <TAxis.h>
-#include <TStyle.h>
#include <TParticle.h>
#include <TParticlePDG.h>
+#include <TPaveText.h>
+#include <TPythia6Calls.h>
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TTree.h>
+#include <TBrowser.h>
// From AliRoot ...
-#include "AliEMCALJetFinder.h"
-#include "AliEMCALFast.h"
-#include "AliEMCALGeometry.h"
-#include "AliEMCALHit.h"
+#include "AliEMCAL.h"
#include "AliEMCALDigit.h"
#include "AliEMCALDigitizer.h"
+#include "AliEMCALFast.h"
+#include "AliEMCALGeometry.h"
#include "AliEMCALHadronCorrection.h"
-#include "AliRun.h"
+#include "AliEMCALHit.h"
+#include "AliEMCALJetFinder.h"
+#include "AliEMCALJetMicroDst.h"
+#include "AliHeader.h"
#include "AliMagF.h"
#include "AliMagFCM.h"
-#include "AliEMCAL.h"
-#include "AliHeader.h"
-#include "AliPDG.h"
-
+#include "AliRun.h"
+#include "AliGenerator.h"
+#include "AliEMCALGetter.h"
// Interface to FORTRAN
#include "Ecommon.h"
fNjets = 0;
fLego = 0;
fLegoB = 0;
+
fTrackList = 0;
fPtT = 0;
fEtaT = 0;
fPhiT = 0;
+ fPdgT = 0;
+
+ fTrackListB = 0;
fPtB = 0;
fEtaB = 0;
fPhiB = 0;
+ fPdgB = 0;
+
fHCorrection = 0;
fHadronCorrector = 0;
fInFile = 0;
fEvent = 0;
+ fRandomBg = 0;
+
SetParametersForBgSubtraction();
}
fEtaCell[i] = 0.;
fPhiCell[i] = 0.;
}
- fLego = 0;
+ fLego = 0;
+ fLegoB = 0;
+
fTrackList = 0;
fPtT = 0;
fEtaT = 0;
fPhiT = 0;
+ fPdgT = 0;
+
+ fTrackListB = 0;
fPtB = 0;
fEtaB = 0;
fPhiB = 0;
+ fPdgB = 0;
+
fHCorrection = 0;
fHadronCorrector = 0;
fBackground = 0;
SetEfficiencySim();
SetDebug();
SetHadronCorrection();
- SetSamplingFraction();
SetIncludeK0andN();
+ fRandomBg = 0;
SetParametersForBgSubtraction();
}
fJets->Delete();
delete fJets;
}
+ delete fLego;
+ delete fLegoB;
+ delete fhLegoTracks;
+ delete fhLegoEMCAL;
+ delete fhLegoHadrCorr;
+ delete fhEff;
+ delete fhCellEt;
+ delete fhCellEMCALEt;
+ delete fhTrackPt;
+ delete fhTrackPtBcut;
+ delete fhChPartMultInTpc;
+
+ delete[] fTrackList;
+ delete[] fPtT;
+ delete[] fEtaT;
+ delete[] fPhiT;
+ delete[] fPdgT;
+
+ delete[] fTrackListB;
+ delete[] fPtB;
+ delete[] fEtaB;
+ delete[] fPhiB;
+ delete[] fPdgB;
}
-
-
-
#ifndef WIN32
# define jet_finder_ua1 jet_finder_ua1_
# define hf1 hf1_
# define type_of_call
#else
-# define jet_finder_ua1 J
+# define jet_finder_ua1 JET_FINDER_UA1
# define hf1 HF1
# define type_of_call _stdcall
#endif
extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
-
-
void AliEMCALJetFinder::Init()
{
//
//
//
// Geometry
- AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
- AliEMCALGeometry* geom =
- AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+ //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+ // AliEMCALGeometry* geom =
+ // AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ AliEMCALGeometry* geom = gime->EMCALGeometry() ;
+
+// SetSamplingFraction(geom->GetSampling());
+
fNbinEta = geom->GetNZ();
fNbinPhi = geom->GetNPhi();
fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
fEtaMin = geom->GetArm1EtaMin();
fEtaMax = geom->GetArm1EtaMax();
- fDphi = (fPhiMax-fPhiMin)/fNbinEta;
+ fDphi = (fPhiMax-fPhiMin)/fNbinPhi;
fDeta = (fEtaMax-fEtaMin)/fNbinEta;
fNtot = fNbinPhi*fNbinEta;
+ fWeightingMethod = kFALSE;
+
//
SetCellSize(fDeta, fDphi);
//
// I/O
if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
+//
+//
}
void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
Int_t mode = fMode;
Float_t prec_bg = fPrecBg;
Int_t ierror;
-
ResetJets(); // 4-feb-2002 by PAI
jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
for (Int_t nj=0; nj<njet; nj++)
{
-
+ if (fWeightingMethod)
+ {
+ fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
+ JetPhiW(nj),
+ JetEtaW(nj));
+
+ }else{
+
fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
JetPhiW(nj),
JetEtaW(nj));
+ }
+ fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod);
+ fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) );
+ fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj) ));
+ fJetT[nj]->SetHCEnergy(HCConeEnergy( JetEtaW(nj),JetPhiW(nj) ));
+
}
FindTracksInJetCone();
fEvent++;
}
+
+Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
+{
+Float_t newenergy = 0.0;
+Float_t bineta,binphi;
+TAxis *x = fhLegoEMCAL->GetXaxis();
+TAxis *y = fhLegoEMCAL->GetYaxis();
+for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
+ {
+ for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
+ {
+ bineta = x->GetBinCenter(i);
+ binphi = y->GetBinCenter(j);
+ if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
+ {
+ newenergy += fhLegoEMCAL->GetBinContent(i,j);
+ }
+ }
+}
+return newenergy;
+}
+
+Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
+{
+Float_t newenergy = 0.0;
+Float_t bineta,binphi;
+TAxis *x = fhLegoTracks->GetXaxis();
+TAxis *y = fhLegoTracks->GetYaxis();
+for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
+ {
+ for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
+ {
+ bineta = x->GetBinCenter(i);
+ binphi = y->GetBinCenter(j);
+ if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
+ {
+ newenergy += fhLegoTracks->GetBinContent(i,j);
+ }
+ }
+}
+return newenergy;
+}
+
+Float_t AliEMCALJetFinder::HCConeEnergy(Float_t eta, Float_t phi)
+{
+//Float_t newenergy = 0.0;
+//Float_t bineta,binphi;
+//TAxis *x = fhLegoTracks->GetXaxis();
+//TAxis *y = fhLegoTracks->GetYaxis();
+//for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
+// {
+// for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
+// {
+// bineta = x->GetBinCenter(i);
+// binphi = y->GetBinCenter(j);
+// if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
+// {
+// newenergy += fhLegoTracks->GetBinContent(i,j);
+// }
+// }
+//}
+//return newenergy;
+
+return 0.0;
+
+}
+
+
+
+Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
+{
+
+
+Float_t newenergy = 0.0;
+Float_t bineta,binphi;
+TAxis *x = fhLegoEMCAL->GetXaxis();
+TAxis *y = fhLegoEMCAL->GetYaxis();
+
+
+for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord
+{
+ for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord
+ {
+ bineta = x->GetBinCenter(i);
+ binphi = y->GetBinCenter(j);
+ if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
+ {
+ newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
+ }
+ }
+}
+
+return newenergy;
+
+}
+
+
Int_t AliEMCALJetFinder::Njets()
{
// Get number of reconstructed jets
}
// I/O
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
if (!fOutFileName) {
//
// output written to input file
//
- AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
- TTree* pK = gAlice->TreeK();
- file = (pK->GetCurrentFile())->GetName();
+ AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
+ TTree* pK = gAlice->TreeK();
+ file = (pK->GetCurrentFile())->GetName();
+ TBranch * jetBranch ;
if (fDebug > 1)
printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
- if (fJets && gAlice->TreeR()) {
- pEMCAL->MakeBranchInTree(gAlice->TreeR(),
- "EMCALJets",
- &fJets,
- kBufferSize,
- file);
+ //if (fJets && gAlice->TreeR()) {
+ if (fJets && gime->TreeR()) {
+ // pEMCAL->MakeBranchInTree(gAlice->TreeR(),
+ jetBranch = gime->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ;
+ //pEMCAL->MakeBranchInTree(gime->TreeR(),
+ // "EMCALJets",
+ // &fJets,
+ // kBufferSize,
+ // file);
+
+ //Int_t nev = gAlice->GetHeader()->GetEvent();
+ //gAlice->TreeR()->Fill();
+ jetBranch->Fill();
+ //char hname[30];
+ //sprintf(hname,"TreeR%d", nev);
+ //gAlice->TreeR()->Write(hname);
+ //gAlice->TreeR()->Reset();
+ gime->WriteRecPoints("OVERWRITE");
}
- Int_t nev = gAlice->GetHeader()->GetEvent();
- gAlice->TreeR()->Fill();
- char hname[30];
- sprintf(hname,"TreeR%d", nev);
- gAlice->TreeR()->Write(hname);
- gAlice->TreeR()->Reset();
} else {
//
// Output written to user specified output file
//
- TTree* pK = gAlice->TreeK();
- fInFile = pK->GetCurrentFile();
-
- fOutFile->cd();
- char hname[30];
- sprintf(hname,"TreeR%d", fEvent);
- TTree* treeJ = new TTree(hname, "EMCALJets");
- treeJ->Branch("EMCALJets", &fJets, kBufferSize);
- treeJ->Fill();
- treeJ->Write(hname);
- fInFile->cd();
+ //TTree* pK = gAlice->TreeK();
+ TTree* pK = gAlice->TreeK();
+ fInFile = pK->GetCurrentFile();
+
+ fOutFile->cd();
+ char hname[30];
+ sprintf(hname,"TreeR%d", fEvent);
+ TTree* treeJ = new TTree(hname, "EMCALJets");
+ treeJ->Branch("EMCALJets", &fJets, kBufferSize);
+ treeJ->Fill();
+ treeJ->Write(hname);
+ fInFile->cd();
}
ResetJets();
}
void AliEMCALJetFinder::BookLego()
{
//
-// Book histo for discretisation
+// Book histo for discretization
//
//
// Don't add histos to the current directory
- TH2::AddDirectory(0);
- TH1::AddDirectory(0);
+ if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
+
+ // TH2::AddDirectory(0); // hists wil be put to the list from gROOT
+ // TH1::AddDirectory(0);
+ gROOT->cd();
//
// Signal map
fLego = new TH2F("legoH","eta-phi",
fNbinPhi, fPhiMin, fPhiMax);
//
// Background map
- fLegoB = new TH2F("legoB","eta-phi",
+ fLegoB = new TH2F("legoB","eta-phi for BG event",
fNbinEta, fEtaMin, fEtaMax,
fNbinPhi, fPhiMin, fPhiMax);
eTmp.GetSize()-1, eTmp.GetArray());
fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
eTmp.GetSize()-1, eTmp.GetArray());
+
+ fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
+ "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
+
+ fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
+ TAxis *xax = fhSinTheta->GetXaxis();
+ for(Int_t i=1; i<=fNbinEta; i++) {
+ Double_t eta = xax->GetBinCenter(i);
+ fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
+ }
+
+ //! first canvas for drawing
+ fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
}
void AliEMCALJetFinder::DumpLego()
Float_t e, eH;
for (Int_t i = 1; i <= fNbinEta; i++) {
for (Int_t j = 1; j <= fNbinPhi; j++) {
+
e = fLego->GetBinContent(i,j);
- if (e > 0.0) {
- Float_t eta = Xaxis->GetBinCenter(i);
- Float_t phi = Yaxis->GetBinCenter(j);
- fEtCell[fNcell] = e;
- fEtaCell[fNcell] = eta;
- fPhiCell[fNcell] = phi;
- fNcell++;
- fhCellEt->Fill(e);
- }
+ if (fRandomBg) {
+ if (gRandom->Rndm() < 0.5) {
+ Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.));
+ e += ebg;
+ }
+ }
+ if (e > 0.0) e -= fMinCellEt;
+ if (e < 0.0) e = 0.;
+ Float_t eta = Xaxis->GetBinCenter(i);
+ Float_t phi = Yaxis->GetBinCenter(j);
+ fEtCell[fNcell] = e;
+ fEtaCell[fNcell] = eta;
+ fPhiCell[fNcell] = phi;
+ fNcell++;
+ fhCellEt->Fill(e);
if(fhLegoEMCAL) {
eH = fhLegoEMCAL->GetBinContent(i,j);
if(eH > 0.0) fhCellEMCALEt->Fill(eH);
}
} // phi
} // eta
- fNcell--;
+// today
+// fNcell--;
}
void AliEMCALJetFinder::ResetMap()
void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
{
+// Which generator
+//
+ const char* name = gAlice->Generator()->GetName();
+ enum {kPythia, kHijing, kHijingPara};
+ Int_t genType = 0;
+
+ if (!strcmp(name, "Hijing")){
+ genType = kHijing;
+ } else if (!strcmp(name, "Pythia")) {
+ genType = kPythia;
+ } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) {
+ genType = kHijingPara;
+ }
+ if (fDebug>=2)
+ printf("Fill tracks generated by %s %d\n", name, genType);
+
+
//
// Fill Cells with track information
//
if (fDebug >= 2)
printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
+ fNChTpc = 0;
ResetMap();
//
// Access particle information
Int_t npart = (gAlice->GetHeader())->GetNprimary();
- if (fDebug >= 2 || npart<=0) printf(" : npart %i\n", npart);
+ Int_t ntr = (gAlice->GetHeader())->GetNtrack();
+ printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n",
+ npart, ntr, fLego->Integral());
// Create track list
//
if (fPtT) delete[] fPtT;
if (fEtaT) delete[] fEtaT;
if (fPhiT) delete[] fPhiT;
+ if (fPdgT) delete[] fPdgT;
fTrackList = new Int_t [npart];
fPtT = new Float_t[npart];
fEtaT = new Float_t[npart];
fPhiT = new Float_t[npart];
+ fPdgT = new Int_t[npart];
- fNt = npart;
+ fNt = npart;
fNtS = 0;
Float_t chTmp=0.0; // charge of current particle
+ // Int_t idGeant;
+ // this is for Pythia ??
for (Int_t part = 0; part < npart; part++) {
TParticle *MPart = gAlice->Particle(part);
Int_t mpart = MPart->GetPdgCode();
Float_t eta = -100.;
if(pT > 0.001) eta = MPart->Eta();
Float_t theta = MPart->Theta();
+ if (fDebug>=2 && MPart->GetStatusCode()==1) {
+ printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
+ part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
+ }
- if (fDebug > 1) {
+ if (fDebug >= 2 && genType == kPythia) {
if (part == 6 || part == 7)
{
printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
part-5, pT, eta, phi);
}
-
-// if (mpart == 21)
-
-// printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
-// part, mpart, pT, eta, phi);
}
fTrackList[part] = 0;
fPtT[part] = pT; // must be change after correction for resolution !!!
fEtaT[part] = eta;
fPhiT[part] = phi;
+ fPdgT[part] = mpart;
+
+// final state particles only
- if (part < 2) continue;
-
- // move to fLego->Fill because hadron correction must apply
- // if particle hit to EMCAL - 11-feb-2002
- // if (pT == 0 || pT < fPtCut) continue;
+ if (genType == kPythia) {
+ if (MPart->GetStatusCode() != 1) continue;
+ } else if (genType == kHijing) {
+ if (child1 >= 0 && child1 < npart) continue;
+ }
+
+
TParticlePDG* pdgP = 0;
// charged or neutral
pdgP = MPart->GetPDG();
chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
+
if (ich == 0) {
if (chTmp == 0) {
if (!fK0N) {
mpart != kK0Long) continue;
}
}
+ } else if (ich == 2) {
+ if (mpart == kNeutron ||
+ mpart == kNeutronBar ||
+ mpart == kK0Long) continue;
}
-// skip partons
- if (TMath::Abs(mpart) <= 6 ||
- mpart == 21 ||
- mpart == 92) continue;
-// acceptance cut
- if (TMath::Abs(eta) > 0.7) continue;
-// Initial phi may be out of acceptance but track may
-// hit two the acceptance - see variable curls
-// if (phi*180./TMath::Pi() > 120.) continue;
+
+ if (TMath::Abs(eta)<=0.9) fNChTpc++;
// final state only
if (child1 >= 0 && child1 < npart) continue;
- // printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
- //part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
+// acceptance cut
+ if (eta > fEtaMax || eta < fEtaMin) continue;
+ if (phi > fPhiMax || phi < fPhiMin ) continue;
//
- if (fDebug > 1)
+ if (fDebug >= 3)
printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
part, mpart, child1, eta, phi, pT, chTmp);
//
Float_t pw;
if (fSmear && TMath::Abs(chTmp)) {
pw = AliEMCALFast::SmearMomentum(1,p);
- // p changed - take into account when calculate pT,
- // pz and so on .. ?
+ // p changed - take into account when calculate pT,
+ // pz and so on .. ?
pT = (pw/p) * pT;
- if(fDebug > 1) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
+ if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
p = pw;
}
//
// Tracking Efficiency and TPC acceptance goes here ...
Float_t eff;
if (fEffic && TMath::Abs(chTmp)) {
- // eff = AliEMCALFast::Efficiency(1,p);
- eff = 0.95; // for testing 25-feb-2002
+ eff = 0.9; // AliEMCALFast::Efficiency(2,p);
if(fhEff) fhEff->Fill(p, eff);
if (AliEMCALFast::RandomReject(eff)) {
- if(fDebug > 1) printf(" reject due to unefficiency ");
- continue;
+ if(fDebug >= 5) printf(" reject due to unefficiency ");
+ continue;
}
}
//
Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
if(TMath::Abs(chTmp)) {
- // hadr. correction only for charge particle
- dphi = PropagatePhi(pT, chTmp, curls);
- phiHC = phi + dphi;
- if (fDebug >= 2) {
- printf("\n Delta phi %f pT %f ", dphi, pT);
- if (curls) printf("\n !! Track is curling");
- }
- if(!curls) fhTrackPtBcut->Fill(pT);
-
- if (fHCorrection && !curls) {
- if (!fHadronCorrector)
- Fatal("AliEMCALJetFinder",
- "Hadronic energy correction required but not defined !");
-
- dpH = fHadronCorrector->GetEnergy(p, eta, 7);
- eTdpH = dpH*TMath::Sin(theta);
-
- if (fDebug >= 2) printf(" phi %f phiHC %f eTcorr %f\n",
- phi, phiHC, -eTdpH); // correction is negative
- fLego->Fill(eta, phiHC, -eTdpH);
- fhLegoHadrCorr->Fill(eta, phiHC, eTdpH);
- }
+ // hadr. correction only for charge particle
+ dphi = PropagatePhi(pT, chTmp, curls);
+ phiHC = phi + dphi;
+ if (fDebug >= 6) {
+ printf("\n Delta phi %f pT %f ", dphi, pT);
+ if (curls) printf("\n !! Track is curling");
+ }
+ if(!curls) fhTrackPtBcut->Fill(pT);
+
+ if (fHCorrection && !curls) {
+ if (!fHadronCorrector)
+ Fatal("AliEMCALJetFinder",
+ "Hadronic energy correction required but not defined !");
+
+ dpH = fHadronCorrector->GetEnergy(p, eta, 7);
+ eTdpH = dpH*TMath::Sin(theta);
+
+ if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
+ phi, phiHC, -eTdpH); // correction is negative
+ Int_t xbin,ybin;
+ xbin = fLego->GetXaxis()->FindBin(eta);
+ ybin = fLego->GetYaxis()->FindBin(phiHC);
+ cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
+ fLego->Fill(eta, phi, -fSamplingF*eTdpH );
+ cout <<"Hadron Correction affected bin - contents after correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
+ fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
+ }
}
//
// More to do ? Just think about it !
//
-
- if (phi*180./TMath::Pi() > 120.) continue;
+ if (phi > fPhiMax || phi < fPhiMin ) continue;
if(TMath::Abs(chTmp) ) { // charge particle
- if (pT > fPtCut && !curls) {
- if (fDebug >= 2) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
- eta , phi, pT);
- fLego->Fill(eta, phi, pT);
- fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
- fTrackList[part] = 1;
- fNtS++;
- }
- } else if(ich==0 && fK0N) {
- // case of n, nbar and K0L
- if (fDebug >= 2) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f)\n",
- eta , phi, pT);
+ if (pT > fPtCut && !curls) {
+ if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
+ eta , phi, pT, fNtS);
+ fLego->Fill(eta, phi, pT);
+ fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
+ fTrackList[part] = 1;
+ fNtS++;
+ }
+ } else if(ich > 0 || fK0N) {
+ // case of n, nbar and K0L
+ if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
+ eta , phi, pT, fNtS);
fLego->Fill(eta, phi, pT);
fTrackList[part] = 1;
fNtS++;
}
-
} // primary loop
+ Float_t etsum = 0.;
+ for(Int_t i=0; i<fLego->GetSize(); i++) {
+ Float_t etc = (*fLego)[i];
+ if (etc > fMinCellEt) etsum += etc;
+ }
+
+ printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
+ printf(" Track selected(fNtS) %i \n", fNtS);
+
DumpLego();
}
//
// Access hit information
AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
-
- TTree *treeH = gAlice->TreeH();
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ TTree *treeH = gime->TreeH();
Int_t ntracks = (Int_t) treeH->GetEntries();
//
// Loop over tracks
//
Int_t nbytes = 0;
- Double_t etH = 0.0;
+ // Double_t etH = 0.0;
for (Int_t track=0; track<ntracks;track++) {
gAlice->ResetHits();
Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
Float_t phi = TMath::ATan2(y,x);
- if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
+ if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
- etH = fSamplingF*eloss*TMath::Sin(theta);
- fLego->Fill(eta, phi, etH);
- // fhLegoEMCAL->Fill(eta, phi, etH);
+ // etH = fSamplingF*eloss*TMath::Sin(theta);
+ fLego->Fill(eta, phi, eloss);
} // Hit Loop
} // Track Loop
- // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
- for(Int_t i=0; i<fLego->GetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i];
+
+ // Transition from deposit energy to eT (eT = de*SF*sin(theta))
+ Double_t etsum = 0;
+ for(Int_t i=1; i<=fLego->GetNbinsX(); i++){ // eta
+ Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
+ for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
+ eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
+ fLego->SetBinContent(i,j,eT);
+ // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
+ fhLegoEMCAL->SetBinContent(i,j,eT);
+ if (eT > fMinCellEt) etsum += eT;
+ }
+ }
+
+ // for(Int_t i=0; i<fLego->GetSize(); i++) {
+ // (*fhLegoEMCAL)[i] = (*fLego)[i];
+ // Float_t etc = (*fLego)[i];
+ // if (etc > fMinCellEt) etsum += etc;
+ // }
+
+ printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
// DumpLego(); ??
}
if (!fLego) BookLego();
if (flag == 0) fLego->Reset();
- Int_t nbytes;
+ Int_t nbytes=0;
//
nbytes += branchDr->GetEntry(0);
//
// Get digitizer parameters
- Float_t towerADCped = digr->GetTowerpedestal();
- Float_t towerADCcha = digr->GetTowerchannel();
- Float_t preshoADCped = digr->GetPreShopedestal();
- Float_t preshoADCcha = digr->GetPreShochannel();
+ Float_t preADCped = digr->GetPREpedestal();
+ Float_t preADCcha = digr->GetPREchannel();
+ Float_t ecADCped = digr->GetECApedestal();
+ Float_t ecADCcha = digr->GetECAchannel();
+ Float_t hcADCped = digr->GetHCApedestal();
+ Float_t hcADCcha = digr->GetHCAchannel();
AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
AliEMCALGeometry* geom =
if (fDebug) {
Int_t ndig = digs->GetEntries();
- printf("\n Number of Digits: %d %d\n", ndig, nent);
- printf("\n Parameters: %f %f %f %f\n",
- towerADCped, towerADCcha, preshoADCped, preshoADCcha );
- printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
+ Info("FillFromDigits","Number of Digits: %d %d\n Parameters: PRE : %f %f EC: %f %f HC: %f %f\n Geometry: %d %d",
+ ndig, nent, preADCped, preADCcha, ecADCped, ecADCcha, hcADCped, hcADCcha, geom->GetNEta(), geom->GetNPhi());
}
//
TIter next(digs);
while ((sdg = (AliEMCALDigit*)(next())))
{
- Double_t pedestal;
- Double_t channel;
- if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
- {
- pedestal = preshoADCped;
- channel = preshoADCcha;
- } else {
- pedestal = towerADCped;
- channel = towerADCcha;
+ Double_t pedestal = 0.;
+ Double_t channel = 0.;
+ if (geom->IsInPRE(sdg->GetId())) {
+ pedestal = preADCped;
+ channel = preADCcha;
+ }
+ else if (geom->IsInECA(sdg->GetId())) {
+ pedestal = ecADCped;
+ channel = ecADCcha;
+ }
+ else if (geom->IsInHCA(sdg->GetId())) {
+ pedestal = hcADCped;
+ channel = hcADCcha;
}
+ else
+ Fatal("FillFromDigits", "unexpected digit is number!") ;
Float_t eta = sdg->GetEta();
Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
- if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
+ if (fDebug > 1)
+ Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
eta, phi, amp, sdg->GetAmp());
fLego->Fill(eta, phi, fSamplingF*amp);
if (fPtT) delete[] fPtT;
if (fEtaT) delete[] fEtaT;
if (fPhiT) delete[] fPhiT;
+ if (fPdgT) delete[] fPdgT;
fPtT = new Float_t[ntracks];
fEtaT = new Float_t[ntracks];
fPhiT = new Float_t[ntracks];
+ fPdgT = new Int_t[ntracks];
fNt = ntracks;
fNtS = 0;
fPtT[track] = pT;
fEtaT[track] = eta;
fPhiT[track] = phi;
-
+ fPdgT[track] = MPart->GetPdgCode();
+
if (track < 2) continue; //Colliding particles?
if (pT == 0 || pT < fPtCut) continue;
fNtS++;
// 26-feb-2002 PAI - for checking all chain
// Work on particles level; accept all particle (not neutrino )
- if (fDebug >= 2)
- printf("\n AliEMCALJetFinder::FillFromParticles()\n");
+ Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law
+ fNChTpc = 0;
ResetMap();
if (!fLego) BookLego();
// Access particles information
Int_t npart = (gAlice->GetHeader())->GetNprimary();
if (fDebug >= 2 || npart<=0) {
- printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
- return;
+ printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
+ if(npart<=0) return;
}
fNt = npart;
fNtS = 0;
RearrangeParticlesMemory(npart);
// Go through the particles
- Int_t mpart, child1;
- Float_t pT, phi, eta;
+ Int_t mpart, child1, child2, geantPdg;
+ Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
TParticle *MPart=0;
for (Int_t part = 0; part < npart; part++) {
fTrackList[part] = 0;
- if(part<8) continue; // skip initial parton configuration
MPart = gAlice->Particle(part);
mpart = MPart->GetPdgCode();
child1 = MPart->GetFirstDaughter();
+ child2 = MPart->GetLastDaughter();
pT = MPart->Pt();
phi = MPart->Phi();
eta = MPart->Eta();
+
+ px = MPart->Px();
+ py = MPart->Py();
+ pz = MPart->Pz();
+ e = MPart->Energy();
+
+// see pyedit in Pythia's text
+ geantPdg = mpart;
+ if (IsThisPartonsOrDiQuark(mpart)) continue;
+ printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
+ part, mpart, geantPdg, px, py, pz, e, child1, child2);
// exclude partons (21 - gluon, 92 - string)
- if ((TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
+
+
// exclude neutrinous also ??
- if (fDebug > 1 && pT>0.01)
- printf("\n part:%5d mpart %5d eta %8.2f phi %8.2f pT %8.2f ",
+ if (fDebug >= 11 && pT>0.01)
+ printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
part, mpart, eta, phi, pT);
fPtT[part] = pT;
fEtaT[part] = eta;
fPhiT[part] = phi;
-// acceptance cut
- if (TMath::Abs(eta) > 0.7) continue;
- if (phi*180./TMath::Pi() > 120.) continue;
+ fPdgT[part] = mpart;
+ fNtS++;
+
// final state only
if (child1 >= 0 && child1 < npart) continue;
+
+ // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
+ // part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
+ PX += px;
+ PY += py;
+ PZ += pz;
+ E += e;
+
+
+ if (TMath::Abs(eta) <= 0.9) fNChTpc++;
+// acceptance cut
+ if (eta > fEtaMax || eta < fEtaMin) continue;
+ if (phi > fPhiMax || phi < fPhiMin ) continue;
//
if(fK0N==0 ) { // exclude neutral hadrons
if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
fLego->Fill(eta, phi, pT);
} // primary loop
+ printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n",
+ PX, PY, PZ, E);
DumpLego();
+ if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
}
void AliEMCALJetFinder::FillFromPartons()
// accept partons before fragmentation - p.57 in Pythia manual
// if(statusCode != 1) continue;
// acceptance cut
- if (TMath::Abs(eta) > 0.7) continue;
- if (phi*180./TMath::Pi() > 120.) continue;
+ if (eta > fEtaMax || eta < fEtaMin) continue;
+ if (phi > fPhiMax || phi < fPhiMin ) continue;
// final state only
// if (child1 >= 0 && child1 < npart) continue;
//
fTrackList[i] = 0;
}
- TTree *treeH = gAlice->TreeH();
- Int_t ntracks = (Int_t) treeH->GetEntries();
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ // TTree *treeH = gAlice->TreeH();
+ TTree *treeH = gime->TreeH();
+ Int_t ntracks = (Int_t) treeH->GetEntries();
//
// Loop over tracks
//
-void AliEMCALJetFinder::SaveBackgroundEvent()
+void AliEMCALJetFinder::SaveBackgroundEvent(Char_t *name)
{
-// Saves the eta-phi lego and the tracklist
-//
- if (fLegoB) fLegoB->Reset();
-
- fLego->Copy(*fLegoB);
+// Saves the eta-phi lego and the tracklist and name of file with BG events
+//
+ if (fLegoB) {
+ fLegoB->Reset();
+ (*fLegoB) = (*fLegoB) + (*fLego);
+ // if(fDebug)
+ printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
+ fLegoB->Integral(), fLego->Integral());
+ }
if (fPtB) delete[] fPtB;
if (fEtaB) delete[] fEtaB;
if (fPhiB) delete[] fPhiB;
+ if (fPdgB) delete[] fPdgB;
if (fTrackListB) delete[] fTrackListB;
fPtB = new Float_t[fNtS];
fEtaB = new Float_t[fNtS];
fPhiB = new Float_t[fNtS];
+ fPdgB = new Int_t [fNtS];
fTrackListB = new Int_t [fNtS];
fNtB = 0;
fPtB [fNtB] = fPtT [i];
fEtaB[fNtB] = fEtaT[i];
fPhiB[fNtB] = fPhiT[i];
+ fPdgB[fNtB] = fPdgT[i];
fTrackListB[fNtB] = 1;
fNtB++;
}
fBackground = 1;
+ printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
+
+ if(strlen(name) == 0) {
+ TSeqCollection *li = gROOT->GetListOfFiles();
+ TString nf;
+ for(Int_t i=0; i<li->GetSize(); i++) {
+ nf = ((TFile*)li->At(i))->GetName();
+ if(nf.Contains("backgorund")) break;
+ }
+ fBGFileName = nf;
+ } else {
+ fBGFileName = name;
+ }
+ printf("BG file name is \n %s\n", fBGFileName.Data());
}
void AliEMCALJetFinder::InitFromBackground()
{
//
//
- if (fDebug) printf("\n Initializing from Background");
+ if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
- if (fLego) fLego->Reset();
- fLegoB->Copy(*fLego);
+ if (fLego) {
+ fLego->Reset();
+ (*fLego) = (*fLego) + (*fLegoB);
+ if(fDebug)
+ printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
+ fLego->Integral(), fLegoB->Integral());
+ } else {
+ printf(" => fLego undefined \n");
+ }
}
} // Background available ?
Int_t nT0 = nT + nTB;
+ printf("Total number of tracks %d\n", nT0);
- if (nT0 > 50) nT0 = 50;
+ if (nT0 > 1000) nT0 = 1000;
Float_t* ptT = new Float_t[nT0];
Float_t* etaT = new Float_t[nT0];
Float_t* phiT = new Float_t[nT0];
+ Int_t* pdgT = new Int_t[nT0];
+
Int_t iT = 0;
Int_t j;
ptT [j+1] = ptT [j];
etaT[j+1] = etaT[j];
phiT[j+1] = phiT[j];
+ pdgT[j+1] = pdgT[j];
}
ptT [index] = fPtT [part];
etaT[index] = fEtaT[part];
phiT[index] = fPhiT[part];
+ pdgT[index] = fPdgT[part];
iT++;
} // particle associated
if (iT > nT0) break;
ptT [j+1] = ptT [j];
etaT[j+1] = etaT[j];
phiT[j+1] = phiT[j];
+ pdgT[j+1] = pdgT[j];
}
ptT [index] = fPtB [part];
etaT[index] = fEtaB[part];
phiT[index] = fPhiB[part];
+ pdgT[index] = fPdgB[part];
iT++;
} // particle associated
if (iT > nT0) break;
} // particle loop
} // Background available ?
- fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT);
+ fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
delete[] ptT;
delete[] etaT;
delete[] phiT;
+ delete[] pdgT;
+
} // jet loop loop
}
// Propagates phi angle to EMCAL radius
//
static Float_t b = 0.0, rEMCAL = -1.0;
- if(rEMCAL<0) {
// Get field in kGS
- b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
+ b = gAlice->Field()->SolenoidField();
// Get EMCAL radius in cm
- rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
- printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
- }
- Float_t dPhi = 0.;
+ rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
+ Float_t dPhi = 0.;
//
//
// bending radies
}
void AliEMCALJetFinder::DrawLego(Char_t *opt)
-{fLego->Draw(opt);}
+{if(fLego) fLego->Draw(opt);}
+
+void AliEMCALJetFinder::DrawLegoBackground(Char_t *opt)
+{if(fLegoB) fLegoB->Draw(opt);}
void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
-{fhLegoEMCAL->Draw(opt);}
+{if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);}
void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
{
if(file==0) file = stdout;
}
fprintf(file,"==== Filling lego ==== \n");
+ fprintf(file,"Sampling fraction %6.3f ", fSamplingF);
fprintf(file,"Smearing %6i ", fSmear);
fprintf(file,"Efficiency %6i\n", fEffic);
fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
} else
fprintf(file,"==== No Bg subtraction ==== \n");
+ //
+ printf("BG file name is %s \n", fBGFileName.Data());
if(file != stdout) fclose(file);
}
if (fPtT) delete[] fPtT;
if (fEtaT) delete[] fEtaT;
if (fPhiT) delete[] fPhiT;
+ if (fPdgT) delete[] fPdgT;
if(npart>0) {
fTrackList = new Int_t [npart];
fPtT = new Float_t[npart];
fEtaT = new Float_t[npart];
fPhiT = new Float_t[npart];
+ fPdgT = new Int_t[npart];
} else {
printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
}
}
+
+Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
+{
+ Int_t absPdg = TMath::Abs(pdg);
+ if(absPdg<=6) return kTRUE; // quarks
+ if(pdg == 21) return kTRUE; // gluon
+ if(pdg == 92) return kTRUE; // string
+
+ // see p.51 of Pythia Manual
+ // Not include diquarks with c and b quark - 4-mar-2002
+ // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1
+ static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303};
+ for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks
+
+ return kFALSE;
+}
+
+void AliEMCALJetFinder::FindChargedJet()
+{
+//
+// Find jet from charged particle information only
+//
+
+//
+// Look for seeds
+//
+ Int_t njets = 0;
+ Int_t part = 0;
+ Int_t nseed = 0;
+
+//
+//
+ ResetJets();
+
+//
+ for (part = 0; part < fNt; part++) {
+ if (!fTrackList[part]) continue;
+ if (fPtT[part] > fEtSeed) nseed++;
+ }
+ printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
+ Int_t* iSeeds = new Int_t[nseed];
+ nseed = 0;
+
+ for (part = 0; part < fNt; part++) {
+ if (!fTrackList[part]) continue;
+ if (fPtT[part] > fEtSeed) iSeeds[nseed++] = part;
+ }
+
+//
+// Loop over seeds
+//
+ Int_t seed = 0;
+ Float_t pt;
+
+ while(1){
+//
+// Find seed with highest pt
+//
+ Float_t ptmax = -1.;
+ Int_t index = -1;
+ Int_t jndex = -1;
+ for (seed = 0; seed < nseed; seed++) {
+ if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
+ ptmax = pt;
+ index = seed;
+ } // ptmax ?
+ } // seeds
+ if (ptmax < 0.) break;
+ jndex = iSeeds[index];
+//
+// Remove from the list
+ iSeeds[index] = -1;
+ printf("\n Next Seed %d %f", jndex, ptmax);
+//
+// Find tracks in cone around seed
+//
+ Float_t phiSeed = fPhiT[jndex];
+ Float_t etaSeed = fEtaT[jndex];
+ Float_t eT = 0.;
+ Float_t pxJ = 0.;
+ Float_t pyJ = 0.;
+ Float_t pzJ = 0.;
+
+ for (part = 0; part < fNt; part++) {
+ if (!fTrackList[part]) continue;
+ Float_t deta = fEtaT[part] - etaSeed;
+ Float_t dphi = fPhiT[part] - phiSeed;
+ Float_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
+ if (dR < fConeRadius) {
+ eT += fPtT[part];
+ Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
+ Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
+ Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
+ Float_t pz = fPtT[part] / TMath::Tan(theta);
+ pxJ += px;
+ pyJ += py;
+ pzJ += pz;
+ //
+ // if seed, remove it
+ //
+ for (seed = 0; seed < nseed; seed++) {
+ if (part == iSeeds[seed]) iSeeds[seed] = -1;
+ } // seed ?
+ } // < cone radius
+ } // particle loop
+
+//
+// Estimate of jet direction
+ Float_t phiJ = TMath::ATan2(pyJ, pxJ);
+ Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
+ Float_t etaJ = TMath::Log(TMath::Tan(thetaJ / 2.));
+// Float_t ptJ = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
+
+//
+// Sum up all energy
+//
+ Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
+ Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
+ Int_t dIphi = Int_t(fConeRadius / fDphi);
+ Int_t dIeta = Int_t(fConeRadius / fDeta);
+ Int_t iPhi, iEta;
+ Float_t sumE = 0;
+ for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
+ for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
+ if (iPhi < 0 || iEta < 0) continue;
+ Float_t dPhi = fPhiMin + iPhi * fDphi;
+ Float_t dEta = fEtaMin + iEta * fDeta;
+ if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
+ sumE += fLego->GetBinContent(iEta, iPhi);
+ } // eta
+ } // phi
+//
+//
+//
+ fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
+ FindTracksInJetCone();
+ printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
+ printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
+ printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
+ } // while(1)
+ EMCALJETS.njet = njets;
+ if (fWrite) WriteJets();
+ fEvent++;
+}
+// 16-jan-2003 - just for convenience
+void AliEMCALJetFinder::Browse(TBrowser* b)
+{
+ if(fHistsList) b->Add((TObject*)fHistsList);
+}
+
+Bool_t AliEMCALJetFinder::IsFolder() const
+{
+ if(fHistsList) return kTRUE;
+ else return kFALSE;
+}
+
+const Char_t* AliEMCALJetFinder::GetNameOfVariant()
+{// generate the literal string with info about jet finder
+ Char_t name[200];
+ sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
+ fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
+ TString nt(name);
+ nt.ReplaceAll(" ","");
+ if(fBGFileName.Length()) {
+ Int_t i1 = fBGFileName.Index("kBackground");
+ Int_t i2 = fBGFileName.Index("/0000") - 1;
+ if(i1>=0 && i2>=0) {
+ TString bg(fBGFileName(i1,i2-i1+1));
+ nt += bg;
+ }
+ }
+ printf("<I> Name of variant %s \n", nt.Data());
+ return nt.Data();
+}