[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinder.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */
//*-- Author: Andreas Morsch (CERN)
#include <TClonesArray.h>
#include <TTree.h>
#include <TFile.h>
#include <TH2.h>
#include <TAxis.h>
#include <TParticle.h>
#include <TParticlePDG.h>
#include "AliEMCALJetFinder.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALHit.h"
#include "Ecommon.h"
#include "AliRun.h"
#include "AliEMCAL.h"
#include "AliHeader.h"

ClassImp(AliEMCALJetFinder)

//____________________________________________________________________________
AliEMCALJetFinder::AliEMCALJetFinder()
{
// Default constructor
    fJets  = 0;
    fNjets = 0;
    fLego  = 0;
}

AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
    : TTask(name, title)
{
// Constructor 
    fJets  = new TClonesArray("AliEMCALJet",10000);
    fNjets = 0;
    for (Int_t i=0; i<30000; i++)
    {
	fEtCell[i]  = 0.;
	fEtaCell[i] = 0.;
        fPhiCell[i] = 0.;
    }
    fLego = 0;
}


//____________________________________________________________________________
AliEMCALJetFinder::~AliEMCALJetFinder()
{
// Destructor
//
    if (fJets){
	fJets->Delete();
	delete fJets;
    }
}


#ifndef WIN32
# define jet_finder_ua1 jet_finder_ua1_
# define hf1 hf1_
# define type_of_call

#else
# define jet_finder_ua1 J
# define hf1 HF1
# define type_of_call _stdcall
#endif

extern "C" void type_of_call 
jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
	       Float_t  etc[30000],  Float_t etac[30000],
	       Float_t  phic[30000], 
	       Float_t& min_move, Float_t& max_move, Int_t& mode, 
	       Float_t& prec_bg,  Int_t& ierror);

extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);


void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], 
			     Float_t etac[30000], Float_t phic[30000],
			     Float_t min_move, Float_t max_move, Int_t mode, 
			     Float_t prec_bg,  Int_t   ierror)
{
// Wrapper for fortran coded jet finder
// Full argument list
    jet_finder_ua1(ncell, ncell_tot, etc, etac, phic, 
		   min_move, max_move, mode, prec_bg, ierror);
    // Write result to output
    WriteJets();
}

void AliEMCALJetFinder::Find()
{
// Wrapper for fortran coded jet finder using member data for 
// argument list

    Float_t min_move = 0;
    Float_t max_move = 0;
    Int_t   mode     = 0;
    Float_t prec_bg  = 0.;
    Int_t   ierror   = 0;

    jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, 
		   min_move, max_move, mode, prec_bg, ierror);
    // Write result to output
    WriteJets();
}


Int_t AliEMCALJetFinder::Njets()
{
// Get number of reconstructed jets
    return EMCALJETS.njet;
}

Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
{
// Get reconstructed jet energy
    return EMCALJETS.etj[i];
}

Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
{
// Get reconstructed jet phi from leading particle
    return EMCALJETS.phij[0][i];
}

Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
{
// Get reconstructed jet phi from weighting
    return EMCALJETS.phij[1][i];
}

Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
{
// Get reconstructed jet eta from leading particles
    return EMCALJETS.etaj[0][i];
}


Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)  
{
// Get reconstructed jet eta from weighting
    return EMCALJETS.etaj[1][i];
}

void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
{
// Set grid cell size
    EMCALCELLGEO.etaCellSize = eta;
    EMCALCELLGEO.phiCellSize = phi;    
}

void AliEMCALJetFinder::SetConeRadius(Float_t par)
{
// Set jet cone radius
    EMCALJETPARAM.coneRad = par;
}

void AliEMCALJetFinder::SetEtSeed(Float_t par)
{
// Set et cut for seeds
    EMCALJETPARAM.etSeed = par;
}

void AliEMCALJetFinder::SetMinJetEt(Float_t par)
{
// Set minimum jet et
    EMCALJETPARAM.ejMin = par;
}

void AliEMCALJetFinder::SetMinCellEt(Float_t par)
{
// Set et cut per cell
    EMCALJETPARAM.etMin = par;
}


void AliEMCALJetFinder::Test()
{
//
// Test the finder call
//
    const Int_t nmax = 30000;
    Int_t ncell      = 10;
    Int_t ncell_tot  = 100;

    Float_t etc[nmax];
    Float_t etac[nmax];
    Float_t phic[nmax];
    Float_t min_move = 0;
    Float_t max_move = 0;
    Int_t   mode = 0;
    Float_t prec_bg = 0;
    Int_t   ierror = 0;

    
    Find(ncell, ncell_tot, etc, etac, phic, 
	 min_move, max_move, mode, prec_bg, ierror);

}

//
//  I/O
//	

void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
{
    //
    // Add a jet 
    //
    TClonesArray &lrawcl = *fJets;
    new(lrawcl[fNjets++]) AliEMCALJet(jet);
}

void AliEMCALJetFinder::ResetJets()
{
    //
    // Reset Jet List 
    //
    fJets->Clear();
    fNjets = 0;
}

void AliEMCALJetFinder::WriteJets()
{
//
// Add all jets to the list
//
    const Int_t kBufferSize = 4000;
    TTree *pK = gAlice->TreeK();
    const char* file = (pK->GetCurrentFile())->GetName();
// I/O
    AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
    printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
    if (fJets && gAlice->TreeR()) {
	pEMCAL->MakeBranchInTree(gAlice->TreeR(), 
				 "Jets", 
				 &fJets, 
				 kBufferSize, 
				 file);
    }
    Int_t njet = Njets();
    for (Int_t nj=0; nj<njet; nj++)
    {
	AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
					   JetPhiW(nj),
					   JetEtaW(nj));
	AddJet(*jet);
	delete jet;
    }
    Int_t nev = gAlice->GetHeader()->GetEvent();
    gAlice->TreeR()->Fill();
    char hname[30];
    sprintf(hname,"TreeR%d", nev);
    gAlice->TreeR()->Write(hname);
    gAlice->TreeR()->Reset();
    ResetJets();        
}

void AliEMCALJetFinder::BookLego()
{
//
//  Book histo for discretisation
//
//
//  Get geometry parameters from 
    AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
    AliEMCALGeometry* geom = 
	AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
    fNbinEta = geom->GetNZ();
    fNbinPhi = geom->GetNPhi();
    const Float_t  phiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
    const Float_t  phiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
    fDphi   = (phiMax-phiMin)/fNbinEta;
    fDeta   = 1.4/fNbinEta;
    fNtot   = fNbinPhi*fNbinEta;
    
//    
    fLego = new TH2F("legoH","eta-phi",
			   fNbinEta, -0.7,  0.7, 
			   fNbinPhi, phiMin, phiMax);
}

void AliEMCALJetFinder::DumpLego()
{
//
// Dump lego histo into array
//
    fNcell = 0;
    for (Int_t i = 1; i < fNbinEta; i++) {
	for (Int_t j = 1; j < fNbinPhi; j++) {
	    Float_t e    = fLego->GetBinContent(i,j);
	    TAxis* Xaxis = fLego->GetXaxis();
	    TAxis* Yaxis = fLego->GetYaxis();
	    Float_t eta  = Xaxis->GetBinCenter(i);
	    Float_t phi  = Yaxis->GetBinCenter(j);	    
	    fEtCell[fNcell]  = e;
	    fEtaCell[fNcell] = eta;
	    fPhiCell[fNcell] = phi;
	    fNcell++;
	} // phi
    } // eta
    fNcell--;
}

void AliEMCALJetFinder::ResetMap()
{
//
// Reset eta-phi array

    for (Int_t i=0; i<30000; i++)
    {
	fEtCell[i]  = 0.;
	fEtaCell[i] = 0.;
	fPhiCell[i] = 0.;
    }
}


void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
{
//
// Fill Cells with hit information
//
//
    ResetMap();
    
    if (!fLego) BookLego();
// Reset
    if (flag == 0) fLego->Reset();
//
// Access particle information    
    Int_t npart = (gAlice->GetHeader())->GetNprimary();
    for (Int_t part=2; part<npart; part++) {
	TParticle *MPart = gAlice->Particle(part);
	Int_t mpart   = MPart->GetPdgCode();
	Int_t child1  = MPart->GetFirstDaughter();
	Float_t pT    = MPart->Pt();
	Float_t phi   = MPart->Phi();
	Float_t theta = MPart->Theta();
	Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
//	if (part == 6 || part == 7)
//	{
//	    printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f", 
//		   part-5, pT, eta, phi);
//	}
	    
	if (pT == 0.) continue;
// charged or neutral 
	if (ich == 0) {
	    TParticlePDG* pdgP = MPart->GetPDG();
	    if (pdgP->Charge() == 0) continue;
	} 
// skip partons
	if (TMath::Abs(mpart) <= 6         ||
	    mpart == 21                    ||
	    mpart == 92) continue;
// acceptance cut
	if (TMath::Abs(eta) > 0.7)         continue;
	if (phi*180./TMath::Pi() > 120.)   continue;
// final state only
	if (child1 >= 0 && child1 < npart) continue;
//	printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f", 
//	part, mpart, child1, eta, phi, pT);
	fLego->Fill(eta, phi, pT);
    } // primary loop
    DumpLego();
}

void AliEMCALJetFinder::FillFromHits(Int_t flag)
{
//
// Fill Cells with track information
//
//
    ResetMap();
    
    if (!fLego) BookLego();
    if (flag == 0) fLego->Reset();

//
// Access hit information    
    AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
//    AliEMCALGeometry* geom = 
//	AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
    
    TTree *treeH = gAlice->TreeH();
    Int_t ntracks = (Int_t) treeH->GetEntries();
//
//   Loop over tracks
//
    Int_t nbytes = 0;

    
    for (Int_t track=0; track<ntracks;track++) {
	gAlice->ResetHits();
	nbytes += treeH->GetEvent(track);
//
//  Loop over hits
//
	for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); 
	    mHit;
	    mHit=(AliEMCALHit*) pEMCAL->NextHit()) 
	{
	    Float_t x      =    mHit->X();         // x-pos of hit
	    Float_t y      =    mHit->Y();         // y-pos
	    Float_t z      =    mHit->Z();         // z-pos
	    Float_t eloss  =    mHit->GetEnergy(); // deposited energy
//	    Int_t   index  =    mHit->GetId();     // cell index
//	    Float_t eta, phi;
//	    geom->EtaPhiFromIndex(index,  eta, phi);
	    Float_t r     = TMath::Sqrt(x*x+y*y);
	    Float_t theta = TMath::ATan2(r,z);
	    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
	    Float_t phi   = TMath::ATan2(y,x);
	    fLego->Fill(eta, phi, eloss);
//	    if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f", 
//	    r, z, eta, phi, eloss);
//	    printf("\n Max %f", fLego->GetMaximum());
	} // Hit Loop
    } // Track Loop
    DumpLego();
}


void hf1(Int_t& id, Float_t& x, Float_t& wgt)
{
// dummy for hbook calls
    ;
}
Commit	Line	Data
471f69dc	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17	//*-- Author: Andreas Morsch (CERN)
	18	#include <TClonesArray.h>
	19	#include <TTree.h>
	20	#include <TFile.h>
	21	#include <TH2.h>
	22	#include <TAxis.h>
	23	#include <TParticle.h>
	24	#include <TParticlePDG.h>
	25	#include "AliEMCALJetFinder.h"
	26	#include "AliEMCALGeometry.h"
	27	#include "AliEMCALHit.h"
	28	#include "Ecommon.h"
	29	#include "AliRun.h"
	30	#include "AliEMCAL.h"
	31	#include "AliHeader.h"
	32
	33	ClassImp(AliEMCALJetFinder)
	34
	35	//____________________________________________________________________________
	36	AliEMCALJetFinder::AliEMCALJetFinder()
	37	{
	38	// Default constructor
	39	fJets = 0;
	40	fNjets = 0;
	41	fLego = 0;
	42	}
	43
	44	AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
	45	: TTask(name, title)
	46	{
	47	// Constructor
	48	fJets = new TClonesArray("AliEMCALJet",10000);
	49	fNjets = 0;
	50	for (Int_t i=0; i<30000; i++)
	51	{
	52	fEtCell[i] = 0.;
	53	fEtaCell[i] = 0.;
	54	fPhiCell[i] = 0.;
	55	}
	56	fLego = 0;
	57	}
	58
	59
	60
	61
	62	//____________________________________________________________________________
	63	AliEMCALJetFinder::~AliEMCALJetFinder()
	64	{
65	// Destructor
66	//
67	if (fJets){
68	fJets->Delete();
69	delete fJets;
70	}
71	}
72
73
74
75	#ifndef WIN32
76	# define jet_finder_ua1 jet_finder_ua1_
77	# define hf1 hf1_
78	# define type_of_call
79
80	#else
81	# define jet_finder_ua1 J
82	# define hf1 HF1
83	# define type_of_call _stdcall
84	#endif
85
86	extern "C" void type_of_call
87	jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
88	Float_t etc[30000], Float_t etac[30000],
89	Float_t phic[30000],
90	Float_t& min_move, Float_t& max_move, Int_t& mode,
91	Float_t& prec_bg, Int_t& ierror);
92
93	extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
94
95
96
97
98
99	void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
100	Float_t etac[30000], Float_t phic[30000],
101	Float_t min_move, Float_t max_move, Int_t mode,
102	Float_t prec_bg, Int_t ierror)
103	{
104	// Wrapper for fortran coded jet finder
105	// Full argument list
106	jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
107	min_move, max_move, mode, prec_bg, ierror);
108	// Write result to output
109	WriteJets();
110	}
111
112	void AliEMCALJetFinder::Find()
113	{
114	// Wrapper for fortran coded jet finder using member data for
115	// argument list
116
117	Float_t min_move = 0;
118	Float_t max_move = 0;
119	Int_t mode = 0;
120	Float_t prec_bg = 0.;
121	Int_t ierror = 0;
122
123	jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
124	min_move, max_move, mode, prec_bg, ierror);
125	// Write result to output
126	WriteJets();
127	}
128
129
130	Int_t AliEMCALJetFinder::Njets()
131	{
132	// Get number of reconstructed jets
133	return EMCALJETS.njet;
134	}
135
136	Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
137	{
138	// Get reconstructed jet energy
139	return EMCALJETS.etj[i];
140	}
141
142	Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
143	{
144	// Get reconstructed jet phi from leading particle
145	return EMCALJETS.phij[0][i];
146	}
147
148	Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
149	{
150	// Get reconstructed jet phi from weighting
151	return EMCALJETS.phij[1][i];
152	}
153
154	Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
155	{
156	// Get reconstructed jet eta from leading particles
157	return EMCALJETS.etaj[0][i];
158	}
159
160
161	Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
162	{
163	// Get reconstructed jet eta from weighting
164	return EMCALJETS.etaj[1][i];
165	}
166
167	void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
168	{
169	// Set grid cell size
170	EMCALCELLGEO.etaCellSize = eta;
171	EMCALCELLGEO.phiCellSize = phi;
172	}
173
174	void AliEMCALJetFinder::SetConeRadius(Float_t par)
175	{
176	// Set jet cone radius
177	EMCALJETPARAM.coneRad = par;
178	}
179
180	void AliEMCALJetFinder::SetEtSeed(Float_t par)
181	{
182	// Set et cut for seeds
183	EMCALJETPARAM.etSeed = par;
184	}
185
186	void AliEMCALJetFinder::SetMinJetEt(Float_t par)
187	{
188	// Set minimum jet et
189	EMCALJETPARAM.ejMin = par;
190	}
191
192	void AliEMCALJetFinder::SetMinCellEt(Float_t par)
193	{
194	// Set et cut per cell
195	EMCALJETPARAM.etMin = par;
196	}
197
198
199	void AliEMCALJetFinder::Test()
200	{
201	//
202	// Test the finder call
203	//
204	const Int_t nmax = 30000;
205	Int_t ncell = 10;
206	Int_t ncell_tot = 100;
207
208	Float_t etc[nmax];
209	Float_t etac[nmax];
210	Float_t phic[nmax];
211	Float_t min_move = 0;
212	Float_t max_move = 0;
213	Int_t mode = 0;
214	Float_t prec_bg = 0;
215	Int_t ierror = 0;
216
217
218	Find(ncell, ncell_tot, etc, etac, phic,
219	min_move, max_move, mode, prec_bg, ierror);
220
221	}
222
223	//
224	// I/O
225	//
226
227	void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
228	{
229	//
230	// Add a jet
231	//
232	TClonesArray &lrawcl = *fJets;
233	new(lrawcl[fNjets++]) AliEMCALJet(jet);
234	}
235
236	void AliEMCALJetFinder::ResetJets()
237	{
238	//
239	// Reset Jet List
240	//
241	fJets->Clear();
242	fNjets = 0;
243	}
244
245	void AliEMCALJetFinder::WriteJets()
246	{
247	//
248	// Add all jets to the list
249	//
250	const Int_t kBufferSize = 4000;
251	TTree *pK = gAlice->TreeK();
252	const char* file = (pK->GetCurrentFile())->GetName();
253	// I/O
254	AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
255	printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
256	if (fJets && gAlice->TreeR()) {
257	pEMCAL->MakeBranchInTree(gAlice->TreeR(),
258	"Jets",
259	&fJets,
260	kBufferSize,
261	file);
262	}
263	Int_t njet = Njets();
264	for (Int_t nj=0; nj<njet; nj++)
265	{
266	AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
267	JetPhiW(nj),
268	JetEtaW(nj));
269	AddJet(*jet);
270	delete jet;
271	}
272	Int_t nev = gAlice->GetHeader()->GetEvent();
273	gAlice->TreeR()->Fill();
274	char hname[30];
275	sprintf(hname,"TreeR%d", nev);
276	gAlice->TreeR()->Write(hname);
277	gAlice->TreeR()->Reset();
278	ResetJets();
279	}
280
281	void AliEMCALJetFinder::BookLego()
282	{
283	//
284	// Book histo for discretisation
285	//
286	//
287	// Get geometry parameters from
288	AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
289	AliEMCALGeometry* geom =
290	AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
291	fNbinEta = geom->GetNZ();
292	fNbinPhi = geom->GetNPhi();
293	const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
294	const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
295	fDphi = (phiMax-phiMin)/fNbinEta;
296	fDeta = 1.4/fNbinEta;
297	fNtot = fNbinPhi*fNbinEta;
298
299	//
300	fLego = new TH2F("legoH","eta-phi",
301	fNbinEta, -0.7, 0.7,
302	fNbinPhi, phiMin, phiMax);
303	}
304
305	void AliEMCALJetFinder::DumpLego()
306	{
307	//
308	// Dump lego histo into array
309	//
310	fNcell = 0;
311	for (Int_t i = 1; i < fNbinEta; i++) {
312	for (Int_t j = 1; j < fNbinPhi; j++) {
313	Float_t e = fLego->GetBinContent(i,j);
314	TAxis* Xaxis = fLego->GetXaxis();
315	TAxis* Yaxis = fLego->GetYaxis();
316	Float_t eta = Xaxis->GetBinCenter(i);
317	Float_t phi = Yaxis->GetBinCenter(j);
318	fEtCell[fNcell] = e;
319	fEtaCell[fNcell] = eta;
320	fPhiCell[fNcell] = phi;
321	fNcell++;
322	} // phi
323	} // eta
324	fNcell--;
325	}
326
327	void AliEMCALJetFinder::ResetMap()
328	{
329	//
330	// Reset eta-phi array
331
332	for (Int_t i=0; i<30000; i++)
333	{
334	fEtCell[i] = 0.;
335	fEtaCell[i] = 0.;
336	fPhiCell[i] = 0.;
337	}
338	}
339
340
341	void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
342	{
343	//
344	// Fill Cells with hit information
345	//
346	//
347	ResetMap();
348
349	if (!fLego) BookLego();
350	// Reset
351	if (flag == 0) fLego->Reset();
352	//
353	// Access particle information
354	Int_t npart = (gAlice->GetHeader())->GetNprimary();
355	for (Int_t part=2; part<npart; part++) {
356	TParticle *MPart = gAlice->Particle(part);
357	Int_t mpart = MPart->GetPdgCode();
358	Int_t child1 = MPart->GetFirstDaughter();
359	Float_t pT = MPart->Pt();
360	Float_t phi = MPart->Phi();
361	Float_t theta = MPart->Theta();
362	Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
363	// if (part == 6 \|\| part == 7)
364	// {
365	// printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
366	// part-5, pT, eta, phi);
367	// }
368
369	if (pT == 0.) continue;
370	// charged or neutral
371	if (ich == 0) {
372	TParticlePDG* pdgP = MPart->GetPDG();
373	if (pdgP->Charge() == 0) continue;
374	}
375	// skip partons
376	if (TMath::Abs(mpart) <= 6 \|\|
377	mpart == 21 \|\|
378	mpart == 92) continue;
379	// acceptance cut
380	if (TMath::Abs(eta) > 0.7) continue;
381	if (phi*180./TMath::Pi() > 120.) continue;
382	// final state only
383	if (child1 >= 0 && child1 < npart) continue;
384	// printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
385	// part, mpart, child1, eta, phi, pT);
386	fLego->Fill(eta, phi, pT);
387	} // primary loop
388	DumpLego();
389	}
390
391	void AliEMCALJetFinder::FillFromHits(Int_t flag)
392	{
393	//
394	// Fill Cells with track information
395	//
396	//
397	ResetMap();
398
399	if (!fLego) BookLego();
400	if (flag == 0) fLego->Reset();
401
402	//
403	// Access hit information
404	AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
405	// AliEMCALGeometry* geom =
406	// AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
407
408	TTree *treeH = gAlice->TreeH();
409	Int_t ntracks = (Int_t) treeH->GetEntries();
410	//
411	// Loop over tracks
412	//
413	Int_t nbytes = 0;
414
415
416	for (Int_t track=0; track<ntracks;track++) {
417	gAlice->ResetHits();
418	nbytes += treeH->GetEvent(track);
419	//
420	// Loop over hits
421	//
422	for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
423	mHit;
424	mHit=(AliEMCALHit*) pEMCAL->NextHit())
425	{
426	Float_t x = mHit->X(); // x-pos of hit
427	Float_t y = mHit->Y(); // y-pos
428	Float_t z = mHit->Z(); // z-pos
429	Float_t eloss = mHit->GetEnergy(); // deposited energy
430	// Int_t index = mHit->GetId(); // cell index
431	// Float_t eta, phi;
432	// geom->EtaPhiFromIndex(index, eta, phi);
433	Float_t r = TMath::Sqrt(xx+yy);
434	Float_t theta = TMath::ATan2(r,z);
435	Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
436	Float_t phi = TMath::ATan2(y,x);
437	fLego->Fill(eta, phi, eloss);
438	// if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f",
439	// r, z, eta, phi, eloss);
440	// printf("\n Max %f", fLego->GetMaximum());
441	} // Hit Loop
442	} // Track Loop
443	DumpLego();
444	}
445
446
447	void hf1(Int_t& id, Float_t& x, Float_t& wgt)
448	{
449	// dummy for hbook calls
450	;
451	}
452
453