[u/mrichter/AliRoot.git] / JETAN / jetan2004 / AliTkConeJetFinderV2.h

// $Id$

// Version 2 of my Cone Jet Finder
// improved storage (see also version 1)

#ifndef ALITKCONJETFINDERV2_H
#define ALITKCONJETFINDERV2_H

#include <Riostream.h>
#include <TParticle.h>
#include <list>

// helper classes
class tower {

 public:
  tower();
  tower(const tower& t);
  tower(Float_t phimin, Float_t phimax, Float_t etamin, Float_t etamax);
  ~tower();

  tower& operator+=(const Float_t E);
  tower& operator+=(const TParticle *part);

  Float_t getEtaMin() const {return fEta_min;}
  Float_t getEtaMax() const {return fEta_max;}
  Float_t getEta()    const {return fEta_center;}
  Float_t getPhiMin() const {return fPhi_min;}
  Float_t getPhiMax() const {return fPhi_max;}
  Float_t getPhi() const {return fPhi_center;}
  Float_t getEt()  const {return fEt;}
  Float_t getEtaWidth() const {return (fEta_max-fEta_min);}
  Float_t getPhiWidth() const {return (fPhi_max-fPhi_min);}

  void setEtaMin(Float_t f) {fEta_min=f;}
  void setEtaMax(Float_t f) {fEta_max=f;}
  void setEta(Float_t f)    {fEta_center=f;}
  void setPhiMin(Float_t f) {fPhi_min=f;}
  void setPhiMax(Float_t f) {fPhi_max=f;}
  void setPhi(Float_t f)    {fPhi_center=f;}

  void addParticle(const TParticle* part) {fParticles.push_back(part);}
  void clearParticles() {fParticles.erase(fParticles.begin(),fParticles.end());}
  void clear() {fEt = 0; clearParticles();}

  list<const TParticle*> *getParticles();

 private:
  Float_t fEta_min;
  Float_t fEta_max;
  Float_t fEta_center;
  Float_t fPhi_min;
  Float_t fPhi_max;
  Float_t fPhi_center;
  Float_t fEt;

  list<const TParticle*> fParticles;
};

ostream& operator<<(ostream& s,const tower& t);

//-----------------------------------------------------

#include "AliTkEtaPhiVector.h"

class protojet {

 public:
  protojet();
  protojet(const protojet& p);
  protojet(const protojet *p);
  ~protojet(){eraseTowers();}

  Float_t Eta() const {return fCentroid.Eta();}
  Float_t Phi() const {return fCentroid.Phi();}
  Float_t getEt() const {
    if(fUpdate) cerr << "Eta update needed" << endl;
    return fEt;
  }  
  Float_t getEt() {
    if(fUpdate) update();
    return fEt;
  }

  Int_t getNTowers() const { return fTowers.size(); }

  void setCentroidPosition(AliTkEtaPhiVector center) {fCentroid = center;}

  AliTkEtaPhiVector getCentroidPosition()   const {
    return fCentroid; 
  }
  AliTkEtaPhiVector getEtWeightedPosition() const {
    if(fUpdate) cerr << "Weighted update needed" << endl;
    return fEtWeightedCentroid;
  }
  AliTkEtaPhiVector getEtWeightedPosition() {
    if(fUpdate) update();
    return fEtWeightedCentroid;
  }

  void addTower(tower *tower) {fTowers.push_back(tower);fUpdate=kTRUE;}
  void eraseTowers() {fTowers.erase(fTowers.begin(),fTowers.end()); }
  void eraseTower(tower *pTower) {fTowers.remove(pTower);}
  void clear(){eraseTowers();}

  list<tower*> getTowerList() const;
  list<tower*> getSharedTowerList(const protojet *other) const;

  bool hasTower(tower *pTower) const;
  bool shareTowers(const protojet *other) const;
  Float_t diffToCenter(tower *pTower) const;
  bool operator<(const protojet &p1) const {return (getEt() < p1.getEt());}
  bool operator==(const protojet &p1);
  void update();

 private:
  AliTkEtaPhiVector fCentroid;
  AliTkEtaPhiVector fEtWeightedCentroid;
  Float_t fEt;
  Bool_t fUpdate;
  list<tower*> fTowers;
};

ostream& operator<<(ostream& s,const protojet& p);

//-----------------------------------------------------
// Real ConeJetFinderClass
//-----------------------------------------------------

#include <TObject.h>
#include <TClonesArray.h>
#include <TTree.h>
#include <vector>
#include "AliTkConeJet.h"
#include "AliTkConeJetEvent.h"
#ifdef ALICEINTERFACE
#include <AliJetEventParticles.h>
#endif

//if defined produce some debugging histos
//#define DOHISTOS

class AliTkConeJetFinderV2 : public TObject {
 public:
  // constructor
  AliTkConeJetFinderV2();
  virtual ~AliTkConeJetFinderV2();
  
  // run control
  void defaultSettings();
  void setSettings(Int_t phibins,Int_t etabins);
  void setEtMinJet(Float_t et){fEtMinJet=et;} //minimum jet energy required
  void setEtCut(Float_t et){fEtCut=et;}       //min et for seedpoints     
  void setPtCut(Float_t pt){fPtCut=pt;}       //set pt cut 
  void setOutput(Bool_t out){fOutput=out;}
  void setRadius(Float_t r=0.7) {fRadius=r;}

  void init();
  void initEvent(const TClonesArray *particles,Int_t type);
  void initEvent(const TClonesArray *particles,Int_t type,TString desc); 
  void run();
  void finishEvent();
  void finish();

  // real physics functions...
  void createTowers();
  void fillTowersFromTParticles(const TClonesArray *particles);
#ifdef ALICEINTERFACE
  void fillTowersFromAliParticles(const TClonesArray *particles);
  void initEvent(const AliJetEventParticles *p,TString desc);
#endif

  void createSeedPoints();
  void findProtojets();
  void findJets();

#ifdef DOHISTOS
  // analysis functions...
  void createHistos();
  void createEventHistos();
  void clearEventHistos();
  void writeEventHistos();
  void writeHistos();
#endif

  // evout function
  void setEvOutFilename(const Char_t *filename);
  Char_t *getEvOutFilename() {return fEvout_name;}

  void addProtojet(protojet *pj);
  void dumpProtojets(Float_t etmin = 5.0);

  void splitMergeJets(protojet *jet1, protojet *jet2);
  void splitJets(protojet *jet1,protojet *jet2);
  void mergeJets(protojet *jet1,protojet *jet2);
  void addJet(protojet *pj) {fJets.push_back(pj);}
  void dumpJets();

  Bool_t isJetEnergy(Float_t min,Float_t max);
  Float_t maxJetEnergy(); //of protojets (before calling finishEvent)
  AliTkConeJetEvent* getEvent() const {return fEvoutevent;}

  //getters
  Bool_t getOutput()    const {return fOutput;}
  Int_t getNTowers()    const {return fNTowers;}
  Int_t getNEtaBins()   const {return fEtaBins;}
  Float_t getEtaMin()   const {return fEtaMin;}
  Float_t getEtaMax()   const {return fEtaMax;}
  Float_t getEtaWidth() const {return fEtaWidth;}
  Int_t getPhiBins()    const {return fPhiBins;}
  Float_t getPhiMin()   const {return fPhiMin;}
  Float_t getPhiMax()   const {return fPhiMax;}
  Float_t getPhiWidth() const {return fPhiWidth;}
  Float_t getEtCut()    const {return fEtCut;};  
  Float_t getPtCut()    const {return fPtCut;};  
  Float_t getEtMinJet() const {return fEtMinJet;}
  Float_t getRadius()   const {return fRadius;}

 protected:
  void initEvent_(const TClonesArray *particles,Int_t type);
  void fillTowerHist();
  Int_t findTower(Float_t phi, Float_t eta);
  Bool_t isTParticleAccepted(TParticle *particle);

  Bool_t fOutput; //outpug some information if true
  Int_t fNTowers; //# of towers
  Int_t fEtaBins; //eta-phi grid
  Float_t fEtaMin;
  Float_t fEtaMax;
  Float_t fEtaWidth;
  Int_t fPhiBins;
  Float_t fPhiMin;
  Float_t fPhiMax;
  Float_t fPhiWidth;
  Float_t fEtCut;    // cut for seedpoints
  Float_t fEtMinJet; // cut for jets
  Float_t fPtCut; // cut for particles
  Float_t fRadius;

  //towers, seeds, protojets and jets
  vector<tower> *fTowers; //!
  list<AliTkEtaPhiVector> fSeedPointsNew; //!
  list<protojet*> fProtojets; //!
  list<protojet*> fJets; //!

  //store results
  TFile *fEvoutfile; //!
  Char_t *fEvout_name; //!
  AliTkConeJetEvent *fEvoutevent; //!
  TTree *fEvouttree; //!

#ifdef ALICEINTERFACE
  TClonesArray *fAliParticles; //!
#endif

#ifdef DOHISTOS
  //histograms
  TObjArray fHistos;
  TObjArray fEventHistos;
  TFile *fHistFile; //!
#endif

  ClassDef(AliTkConeJetFinderV2,1)
};

inline Bool_t AliTkConeJetFinderV2::isTParticleAccepted(TParticle * /* particle */) 
{
  // check if particle is accepted
  // makes sense to write this into a own class, but now I'm lazy

  // check if particle is stable -> a detectable particle
#ifndef ALICEINTERFACE 
  if (particle->GetStatusCode()%100 != 1) {
    return kFALSE;
  }
#endif

  // default case: accept
  return kTRUE;
}

#endif
Commit	Line	Data
b9a6a391	1	// $Id$
	2
	3	// Version 2 of my Cone Jet Finder
	4	// improved storage (see also version 1)
	5
	6	#ifndef ALITKCONJETFINDERV2_H
	7	#define ALITKCONJETFINDERV2_H
	8
	9	#include <Riostream.h>
	10	#include <TParticle.h>
	11	#include <list>
	12
	13	// helper classes
	14	class tower {
	15
	16	public:
	17	tower();
	18	tower(const tower& t);
	19	tower(Float_t phimin, Float_t phimax, Float_t etamin, Float_t etamax);
	20	~tower();
	21
	22	tower& operator+=(const Float_t E);
	23	tower& operator+=(const TParticle *part);
	24
	25	Float_t getEtaMin() const {return fEta_min;}
	26	Float_t getEtaMax() const {return fEta_max;}
	27	Float_t getEta() const {return fEta_center;}
	28	Float_t getPhiMin() const {return fPhi_min;}
	29	Float_t getPhiMax() const {return fPhi_max;}
	30	Float_t getPhi() const {return fPhi_center;}
	31	Float_t getEt() const {return fEt;}
	32	Float_t getEtaWidth() const {return (fEta_max-fEta_min);}
	33	Float_t getPhiWidth() const {return (fPhi_max-fPhi_min);}
	34
	35	void setEtaMin(Float_t f) {fEta_min=f;}
	36	void setEtaMax(Float_t f) {fEta_max=f;}
	37	void setEta(Float_t f) {fEta_center=f;}
	38	void setPhiMin(Float_t f) {fPhi_min=f;}
	39	void setPhiMax(Float_t f) {fPhi_max=f;}
	40	void setPhi(Float_t f) {fPhi_center=f;}
	41
	42	void addParticle(const TParticle* part) {fParticles.push_back(part);}
	43	void clearParticles() {fParticles.erase(fParticles.begin(),fParticles.end());}
	44	void clear() {fEt = 0; clearParticles();}
	45
	46	list<const TParticle> getParticles();
	47
	48	private:
	49	Float_t fEta_min;
	50	Float_t fEta_max;
	51	Float_t fEta_center;
	52	Float_t fPhi_min;
	53	Float_t fPhi_max;
	54	Float_t fPhi_center;
	55	Float_t fEt;
	56
	57	list<const TParticle*> fParticles;
	58	};
	59
	60	ostream& operator<<(ostream& s,const tower& t);
	61
	62	//-----------------------------------------------------
	63
	64	#include "AliTkEtaPhiVector.h"
65
66	class protojet {
67
68	public:
69	protojet();
70	protojet(const protojet& p);
71	protojet(const protojet *p);
72	~protojet(){eraseTowers();}
73
74	Float_t Eta() const {return fCentroid.Eta();}
75	Float_t Phi() const {return fCentroid.Phi();}
76	Float_t getEt() const {
77	if(fUpdate) cerr << "Eta update needed" << endl;
78	return fEt;
79	}
80	Float_t getEt() {
81	if(fUpdate) update();
82	return fEt;
83	}
84
85	Int_t getNTowers() const { return fTowers.size(); }
86
87	void setCentroidPosition(AliTkEtaPhiVector center) {fCentroid = center;}
88
89	AliTkEtaPhiVector getCentroidPosition() const {
90	return fCentroid;
91	}
92	AliTkEtaPhiVector getEtWeightedPosition() const {
93	if(fUpdate) cerr << "Weighted update needed" << endl;
94	return fEtWeightedCentroid;
95	}
96	AliTkEtaPhiVector getEtWeightedPosition() {
97	if(fUpdate) update();
98	return fEtWeightedCentroid;
99	}
100
101	void addTower(tower *tower) {fTowers.push_back(tower);fUpdate=kTRUE;}
102	void eraseTowers() {fTowers.erase(fTowers.begin(),fTowers.end()); }
103	void eraseTower(tower *pTower) {fTowers.remove(pTower);}
104	void clear(){eraseTowers();}
105
106	list<tower*> getTowerList() const;
107	list<tower> getSharedTowerList(const protojet other) const;
108
109	bool hasTower(tower *pTower) const;
110	bool shareTowers(const protojet *other) const;
111	Float_t diffToCenter(tower *pTower) const;
112	bool operator<(const protojet &p1) const {return (getEt() < p1.getEt());}
113	bool operator==(const protojet &p1);
114	void update();
115
116	private:
117	AliTkEtaPhiVector fCentroid;
118	AliTkEtaPhiVector fEtWeightedCentroid;
119	Float_t fEt;
120	Bool_t fUpdate;
121	list<tower*> fTowers;
122	};
123
124	ostream& operator<<(ostream& s,const protojet& p);
125
126	//-----------------------------------------------------
127	// Real ConeJetFinderClass
128	//-----------------------------------------------------
129
130	#include <TObject.h>
131	#include <TClonesArray.h>
132	#include <TTree.h>
133	#include <vector>
134	#include "AliTkConeJet.h"
135	#include "AliTkConeJetEvent.h"
136	#ifdef ALICEINTERFACE
137	#include <AliJetEventParticles.h>
138	#endif
139
140	//if defined produce some debugging histos
141	//#define DOHISTOS
142
143	class AliTkConeJetFinderV2 : public TObject {
144	public:
145	// constructor
146	AliTkConeJetFinderV2();
147	virtual ~AliTkConeJetFinderV2();
148
149	// run control
150	void defaultSettings();
151	void setSettings(Int_t phibins,Int_t etabins);
152	void setEtMinJet(Float_t et){fEtMinJet=et;} //minimum jet energy required
153	void setEtCut(Float_t et){fEtCut=et;} //min et for seedpoints
154	void setPtCut(Float_t pt){fPtCut=pt;} //set pt cut
155	void setOutput(Bool_t out){fOutput=out;}
156	void setRadius(Float_t r=0.7) {fRadius=r;}
157
158	void init();
159	void initEvent(const TClonesArray *particles,Int_t type);
160	void initEvent(const TClonesArray *particles,Int_t type,TString desc);
161	void run();
162	void finishEvent();
163	void finish();
164
165	// real physics functions...
166	void createTowers();
167	void fillTowersFromTParticles(const TClonesArray *particles);
168	#ifdef ALICEINTERFACE
169	void fillTowersFromAliParticles(const TClonesArray *particles);
170	void initEvent(const AliJetEventParticles *p,TString desc);
171	#endif
172
173	void createSeedPoints();
174	void findProtojets();
175	void findJets();
176
177	#ifdef DOHISTOS
178	// analysis functions...
179	void createHistos();
180	void createEventHistos();
181	void clearEventHistos();
182	void writeEventHistos();
183	void writeHistos();
184	#endif
185
186	// evout function
187	void setEvOutFilename(const Char_t *filename);
188	Char_t *getEvOutFilename() {return fEvout_name;}
189
190	void addProtojet(protojet *pj);
191	void dumpProtojets(Float_t etmin = 5.0);
192
193	void splitMergeJets(protojet jet1, protojet jet2);
194	void splitJets(protojet jet1,protojet jet2);
195	void mergeJets(protojet jet1,protojet jet2);
196	void addJet(protojet *pj) {fJets.push_back(pj);}
197	void dumpJets();
198
199	Bool_t isJetEnergy(Float_t min,Float_t max);
200	Float_t maxJetEnergy(); //of protojets (before calling finishEvent)
201	AliTkConeJetEvent* getEvent() const {return fEvoutevent;}
202
203	//getters
204	Bool_t getOutput() const {return fOutput;}
205	Int_t getNTowers() const {return fNTowers;}
206	Int_t getNEtaBins() const {return fEtaBins;}
207	Float_t getEtaMin() const {return fEtaMin;}
208	Float_t getEtaMax() const {return fEtaMax;}
209	Float_t getEtaWidth() const {return fEtaWidth;}
210	Int_t getPhiBins() const {return fPhiBins;}
211	Float_t getPhiMin() const {return fPhiMin;}
212	Float_t getPhiMax() const {return fPhiMax;}
213	Float_t getPhiWidth() const {return fPhiWidth;}
214	Float_t getEtCut() const {return fEtCut;};
215	Float_t getPtCut() const {return fPtCut;};
216	Float_t getEtMinJet() const {return fEtMinJet;}
217	Float_t getRadius() const {return fRadius;}
218
219	protected:
220	void initEvent_(const TClonesArray *particles,Int_t type);
221	void fillTowerHist();
222	Int_t findTower(Float_t phi, Float_t eta);
223	Bool_t isTParticleAccepted(TParticle *particle);
224
225	Bool_t fOutput; //outpug some information if true
226	Int_t fNTowers; //# of towers
227	Int_t fEtaBins; //eta-phi grid
228	Float_t fEtaMin;
229	Float_t fEtaMax;
230	Float_t fEtaWidth;
231	Int_t fPhiBins;
232	Float_t fPhiMin;
233	Float_t fPhiMax;
234	Float_t fPhiWidth;
235	Float_t fEtCut; // cut for seedpoints
236	Float_t fEtMinJet; // cut for jets
237	Float_t fPtCut; // cut for particles
238	Float_t fRadius;
239
240	//towers, seeds, protojets and jets
241	vector<tower> *fTowers; //!
242	list<AliTkEtaPhiVector> fSeedPointsNew; //!
243	list<protojet*> fProtojets; //!
244	list<protojet*> fJets; //!
245
246	//store results
247	TFile *fEvoutfile; //!
248	Char_t *fEvout_name; //!
249	AliTkConeJetEvent *fEvoutevent; //!
250	TTree *fEvouttree; //!
251
252	#ifdef ALICEINTERFACE
253	TClonesArray *fAliParticles; //!
254	#endif
255
256	#ifdef DOHISTOS
257	//histograms
258	TObjArray fHistos;
259	TObjArray fEventHistos;
260	TFile *fHistFile; //!
261	#endif
262
263	ClassDef(AliTkConeJetFinderV2,1)
264	};
265
f114c6bf	266	inline Bool_t AliTkConeJetFinderV2::isTParticleAccepted(TParticle * /* particle */)
b9a6a391	267	{
	268	// check if particle is accepted
	269	// makes sense to write this into a own class, but now I'm lazy
	270
	271	// check if particle is stable -> a detectable particle
	272	#ifndef ALICEINTERFACE
	273	if (particle->GetStatusCode()%100 != 1) {
	274	return kFALSE;
	275	}
	276	#endif
	277
	278	// default case: accept
	279	return kTRUE;
	280	}
	281
	282	#endif