[u/mrichter/AliRoot.git] / JETAN / AliTkKtJetFinder.cxx

// $Id$

// includes 
#include <Riostream.h>
#include <list>
#include <queue>
#include <iterator>

#include <TROOT.h>
#include <TObject.h>
#include <TMath.h>
#include <TH2.h>
#include <TClonesArray.h>
#include <TParticle.h>
#include <TStopwatch.h>

#include "AliTkKtJetFinder.h"

void AliTkKtJetFinder::init() {
  DebugOutput("TkJetFinder::init() called");

  cout << "TkJetFinder::init() - finder initialized" << endl;
  cout << "precluster parameter:" << endl;
  cout << " - phi binning: " << phiBins << " bins from "
       << phiMin << " to " << phiMax << ", bin width "
       << (phiMax-phiMin)/(Float_t)phiBins << endl;
  cout << " - theta binning: " << thetaBins << " bins from "
       << thetaMin << " to " << thetaMax << ", bin width "
       << (thetaMax-thetaMin)/(Float_t)thetaBins << endl;
  cout << "jetfinder parameter:" << endl;
  cout << " - D:" << finder_D << endl
       << " - D_Cut:" << finder_DCut << endl;

  DebugOutput("TkJetFinder::init() finished");
}

void AliTkKtJetFinder::makeTParticles(TClonesArray *particles) {
  DebugOutput("TkJetFinder::makeTParticle() called");
  initEvent();
  preclusterTParticles(particles);
  DebugOutput("TkJetFinder::makeTParticle() preclustering done");
  findJets();
  DebugOutput("TkJetFinder::makeTParticle() jetfinding done");
  DebugOutput("TkJetFinder::makeTParticle() finished");
}

void AliTkKtJetFinder::clear() {
  DebugOutput("TkJetFinder::clear() called");

  DebugOutput("TkJetFinder::clear() finished");
}

void AliTkKtJetFinder::finish() {
  DebugOutput("TkJetFinder::finish() called");
  
  DebugOutput("TkJetFinder::finish() finished");
}

void AliTkKtJetFinder::initEvent() {
  DebugOutput("TkJetFinder::initEvent() called");

  // clean precluster list
  if (firstPreCluster) {
    TPreCluster *act = firstPreCluster;
    TPreCluster *next;
    while (act != NULL) {
      next = act->next;
      delete act;
      act = next;
    }
  }

  // delete precluster array
  delete[] preClusterArray;

  // reset list
  firstPreCluster = NULL;
  lastPreCluster = NULL;
  preClusterUID = 0;
  preClusterArray = NULL;

  // delete dlist;
  for (list<Tdlist *>::iterator i = myDList.begin(); i != myDList.end(); ++i) {
    Tdlist *dlist = *i;
    delete dlist;
    dlist = NULL;
  }
  myDList.clear();
  
  
  DebugOutput("TkJetFinder::initEvent() finished");
}

void AliTkKtJetFinder::preclusterTParticles(TClonesArray *particles) {
  Char_t buffer[1024];
  Int_t acceptedParticles = 0;
  Int_t acceptedPreCluster = 0;

  DetailedOutput("TkJetFinder::preclusterTParticle() called");

  TStopwatch timer;
  timer.Reset();

  sprintf(buffer,"TkJetFinder::preclusterTParticle() - found %d particles",
	  particles->GetEntries());
  DetailedOutput(buffer);

  cout << "WARNING - TkJetFinder::preclusterTParticle()" << endl
       << "          no preclustering implemented - each particle accepted!"
       << endl;

  timer.Start();
  // start the particle loop
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  TPreCluster *precluster;
  while((particle = (TParticle *) iter->Next()) != NULL) {
    if (isTParticleAccepted(particle)) {
      acceptedParticles++;
      precluster = new TPreCluster;
      // to be kind of realistic use maseless particles 
      // - TPC doesn't measure energy and there is no hadron calorimeter...
      precluster->id = preClusterUID;
      precluster->E = TMath::Sqrt(particle->P()*particle->P());
      precluster->px = particle->Px();
      precluster->py = particle->Py();
      precluster->pz = particle->Pz();
      precluster->next = NULL;
      if (true) {
	// add precluster to list of preclusters
	addPreCluster(precluster);
	acceptedPreCluster++;
	preClusterUID++;
      } else {
	delete precluster;
      }
    }
  }

  // build an array of pointers for direct access on the precluster...
  preClusterArray = new TPreCluster*[preClusterUID];
  precluster = firstPreCluster;
  while (precluster) {
    preClusterArray[precluster->id] = precluster;
    precluster = precluster->next;
  }

  timer.Stop();

  // if you want to be really sure what the precluster algorithm does...
  //dumpPreClusters();
  //dumpPreClusterArray();

  // status information at the end...
  sprintf(buffer,"TkJetFinder::preclusterTParticle() - used %d particles",
	  acceptedParticles);
  DebugOutput(buffer);
  sprintf(buffer,
	  "TkJetFinder::preclusterTParticle() - created %d preclusters",
	  acceptedPreCluster);
  DebugOutput(buffer);
  
  sprintf(buffer,
	  "TkJetFinder::preclusterTParticle() - timing CPU: %f, real %f",
	  timer.CpuTime(),timer.RealTime());
  TimingOutput(buffer);
  
  cout << "some stupid tests..." << endl;
  Tdlist *p1 = new Tdlist;
  Tdlist *p2 = new Tdlist;
  Tdlist *p3 = new Tdlist;
  p1->d = 3;
  p2->d = 2;
  p3->d = 1;
  if (p1==p2) {
    cout << "p1==p2 (" << p1->d << "==" << p2->d << ")" << endl;
  } else {
    cout << "p1!=p2 (" << p1->d << "!=" << p2->d << ")" << endl;
  }
  if (p1 < p2) {
    cout << "p1 < p2 (" << p1->d << ">" << p2->d << ")" << endl;
  } else {
    cout << "p1 > p2 (" << p1->d << "<" << p2->d << ")" << endl;
  }
  priority_queue<Tdlist *> myHeap;
  myHeap.push(p2);
  myHeap.push(p1);
  myHeap.push(p3);

  cout << myHeap.size() << endl;
  while (!myHeap.empty()) {
    Tdlist *myEntry = myHeap.top();
    myHeap.pop();
    cout << "d = " << myEntry->d << endl;
    delete myEntry;
  }

  DetailedOutput("TkJetFinder::preclusterTParticle() finished");
}

void AliTkKtJetFinder::findJets() {
  Char_t buffer[1024];

  TStopwatch timer;
  timer.Reset();

  DetailedOutput("TkJetFinder::findJets() called");

  timer.Start();
  // let's start...
  buildNewDList();
  // find the precluster/pair with lowest dmin...
  //list<Tdlist *>::iterator i;
  list<Tdlist *>::iterator min;

  myDList.sort();
  while(!myDList.empty()) {
    min = myDList.begin();

    // we have the entry with the lowest min...
    //cout << (*min)->d << endl;
    delete *min;
    myDList.erase(min);
  }
    

  timer.Stop();

  sprintf(buffer,
	  "TkJetFinder::findJets() - timing CPU: %f, real %f",
	  timer.CpuTime(),timer.RealTime());
  TimingOutput(buffer);

  DetailedOutput("TkJetFinder::findJets() finish");
}

void AliTkKtJetFinder::setPhiBins(Int_t nPhiBins, 
			       Float_t fPhiMin, Float_t fPhiMax) {
  setNPhiBins(nPhiBins);
  setPhiMin(fPhiMin);
  setPhiMax(fPhiMax);
}

void AliTkKtJetFinder::setThetaBins(Int_t nThetaBins,
				 Float_t fThetaMin, Float_t fThetaMax) {
  setNThetaBins(nThetaBins);
  setThetaMin(fThetaMin);
  setThetaMax(fThetaMax);
}


void AliTkKtJetFinder::setDefaultOptions() {
  // precluster options
  setPhiBins((Int_t)(TMath::Pi()*2./0.1), 0. ,2.*TMath::Pi());
  setThetaBins(10,(TMath::Pi()/2. - 0.5),(TMath::Pi()/2. + 0.5));

  // finder options
  setD(1.);
  setDCut(-1.);

  // debug options
  setDebugLevel(0);
  setDebugFilename("KtJetFinderDebug.root");

  cout << "TkJetFinder::setDefaultOptions() - DEFAULT options set" << endl;
}

Bool_t AliTkKtJetFinder::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
  if (particle->GetStatusCode() != 1) {
    return kFALSE;
  }
  // and now just for fun
  TParticlePDG *partPDG;
  partPDG = particle->GetPDG();
//   if (partPDG->Charge() == 0) {
//     cout << partPDG->GetName() << endl;
//     return kFALSE;
//   }

  // default case: accept
  return kTRUE;
}

void AliTkKtJetFinder::addPreCluster(TPreCluster *precluster) {
  if (firstPreCluster == NULL) {
    firstPreCluster = precluster;
    lastPreCluster = precluster;
    return;
  }

  lastPreCluster->next = precluster;
  precluster->next = NULL;
  lastPreCluster = precluster;
}

void AliTkKtJetFinder::deletePreCluster(Int_t UID) {
  TPreCluster *p1 = NULL;
  TPreCluster *p2 = firstPreCluster;
  while (p2!=NULL) {
    if (p2->id == UID) {
      p1->next = p2->next;
      delete p2;
      break;
    }
    p1 = p2;
    p2 = p2->next;
  }
}

void AliTkKtJetFinder::dumpPreClusters() {
  TPreCluster *pre = firstPreCluster;
  while (pre) {
    cout << "UID: " << pre->id
	 << " E: " << pre->E
	 << " px: " << pre->px
	 << " py: " << pre->py
	 << " pz: " << pre->pz
	 << " next: " << pre->next
	 << endl;
    pre = pre->next;
  }
}

void AliTkKtJetFinder::dumpPreClusterArray() {
  for (Int_t i = 0; i < preClusterUID; i++) {
    cout << "UID: " << preClusterArray[i]->id
	 << " E: " << preClusterArray[i]->E
	 << " px: " << preClusterArray[i]->px
	 << " py: " << preClusterArray[i]->py
	 << " pz: " << preClusterArray[i]->pz
	 << " next: " << preClusterArray[i]->next
	 << endl;
  }
}

void AliTkKtJetFinder::addD(Tdlist *newD) {
    myDList.push_back(newD);
}

void AliTkKtJetFinder::buildNewDList() {
  Char_t buffer[1024];

  // step 1 - add all preclusters to dlist
  struct Tdlist *dlist;
  for (Int_t i = 0; i < preClusterUID; i++) {
    dlist = new Tdlist;
    dlist->d = calcD((TPreCluster *)preClusterArray[i]);
    dlist->id1 = preClusterArray[i]->id;
    dlist->id2 = -1;
    dlist->prev = NULL;
    dlist->next = NULL;
    addD(dlist);
  }

  sprintf(buffer,
	  "AliTkKtJetFinder::buildNewDList - added %d preclusters to new list",
	  (Int_t)myDList.size());
  DetailedOutput(buffer);

  // and now for all combinations...
  for (Int_t i = 0; i < preClusterUID; i++) {
    for (Int_t j = i + 1; j < preClusterUID; j++) {
      if (i == j) {
	// makes no sense...
	continue;
      }
      dlist = new Tdlist;
      dlist->d = calcD((TPreCluster *)preClusterArray[i],
		       (TPreCluster *)preClusterArray[j]);
      dlist->id1 = preClusterArray[i]->id;
      dlist->id2 = preClusterArray[j]->id;
      dlist->prev = NULL;
      dlist->next = NULL;
      addD(dlist);
    }
  }

  sprintf(buffer,
	  "AliTkKtJetFinder::buildNewDList - created list with %d entries",
	  (Int_t)myDList.size());
  DetailedOutput(buffer);
}

Float_t AliTkKtJetFinder::calcD(TPreCluster *p1, TPreCluster *p2) {
  if (!p1) {
    // something is wrong - we need a precluster
    return -9999;
  }
  if (!p2) {
    // called with one particle only
    // return pT^2
    return (p1->px*p1->px + p1->py*p1->py);
  } 
  // we have two particles... 
  // calculate the minimum pt...
  Float_t pt1Sq = p1->px*p1->px + p1->py*p1->py;
  Float_t pt2Sq = p2->px*p2->px + p2->py*p2->py;
  Float_t ptMinSq;
  if (pt1Sq < pt2Sq) {
    ptMinSq = pt1Sq;
  } else {
    ptMinSq = pt2Sq;
  }
  // calculate rapidities - !move this into precluster structure!!!
  Float_t y1 = 0.5 * TMath::Log((p1->E+p1->pz)/(p1->E-p1->pz));
  Float_t y2 = 0.5 * TMath::Log((p2->E+p2->pz)/(p2->E-p2->pz));
  // calculate phi - !move this into precluster structure!!!
  Float_t phi1 = 0;
  Float_t phi2 = 0;
  // calculate ((y1-y2)^2 + (phi1-phi2)^2)/(D^2)
  // !save D^2 instead of D to save time!!!
  Float_t R = ((y1-y2)*(y1-y2) + (phi1-phi2)*(phi1-phi2))/(finder_D*finder_D);
  
  return ptMinSq * R;
}

void AliTkKtJetFinder::setDebugLevel(Int_t nDebugLevel) {
  debugLevel = nDebugLevel;

  // description of debug levels
  // levels are additive!!!
  //  0 -> no debug output - only normal output (default)
  //  1 -> normal debug output
  //       major function calls, summary of major results
  //  2 -> timing output
  //       timing information for major steps
  //  4 -> detailed results
  //       output of much more results
  //  8 -> result histos
  //       histograms of event results for internal consitency checks
  // 16 -> timing histos
  //       histograms with timing info - not cleared per event!
  // 32 -> write debug histos to file
}

void AliTkKtJetFinder::DebugOutput(Char_t *output) {
  if ((debugLevel & 1) == 1) {
    cout << output << endl;
  }
}

void AliTkKtJetFinder::TimingOutput(Char_t *output) {
  if ((debugLevel & 2) == 2) {
    cout << output << endl;
  }
}

void AliTkKtJetFinder::DetailedOutput(Char_t *output) {
  if ((debugLevel & 4) == 4) {
    cout << output << endl;
  }
}


ClassImp(AliTkKtJetFinder)
Commit	Line	Data
b9a6a391	1	// $Id$
	2
	3	// includes
	4	#include <Riostream.h>
	5	#include <list>
	6	#include <queue>
	7	#include <iterator>
	8
	9	#include <TROOT.h>
	10	#include <TObject.h>
	11	#include <TMath.h>
	12	#include <TH2.h>
	13	#include <TClonesArray.h>
	14	#include <TParticle.h>
	15	#include <TStopwatch.h>
	16
	17	#include "AliTkKtJetFinder.h"
	18
	19	void AliTkKtJetFinder::init() {
	20	DebugOutput("TkJetFinder::init() called");
	21
	22	cout << "TkJetFinder::init() - finder initialized" << endl;
	23	cout << "precluster parameter:" << endl;
	24	cout << " - phi binning: " << phiBins << " bins from "
	25	<< phiMin << " to " << phiMax << ", bin width "
	26	<< (phiMax-phiMin)/(Float_t)phiBins << endl;
	27	cout << " - theta binning: " << thetaBins << " bins from "
	28	<< thetaMin << " to " << thetaMax << ", bin width "
	29	<< (thetaMax-thetaMin)/(Float_t)thetaBins << endl;
	30	cout << "jetfinder parameter:" << endl;
	31	cout << " - D:" << finder_D << endl
	32	<< " - D_Cut:" << finder_DCut << endl;
	33
	34	DebugOutput("TkJetFinder::init() finished");
	35	}
	36
	37	void AliTkKtJetFinder::makeTParticles(TClonesArray *particles) {
	38	DebugOutput("TkJetFinder::makeTParticle() called");
	39	initEvent();
	40	preclusterTParticles(particles);
	41	DebugOutput("TkJetFinder::makeTParticle() preclustering done");
	42	findJets();
	43	DebugOutput("TkJetFinder::makeTParticle() jetfinding done");
	44	DebugOutput("TkJetFinder::makeTParticle() finished");
	45	}
	46
	47	void AliTkKtJetFinder::clear() {
	48	DebugOutput("TkJetFinder::clear() called");
	49
	50	DebugOutput("TkJetFinder::clear() finished");
	51	}
	52
	53	void AliTkKtJetFinder::finish() {
	54	DebugOutput("TkJetFinder::finish() called");
	55
	56	DebugOutput("TkJetFinder::finish() finished");
	57	}
	58
	59	void AliTkKtJetFinder::initEvent() {
	60	DebugOutput("TkJetFinder::initEvent() called");
	61
	62	// clean precluster list
	63	if (firstPreCluster) {
	64	TPreCluster *act = firstPreCluster;
65	TPreCluster *next;
66	while (act != NULL) {
67	next = act->next;
68	delete act;
69	act = next;
70	}
71	}
72
73	// delete precluster array
74	delete[] preClusterArray;
75
76	// reset list
77	firstPreCluster = NULL;
78	lastPreCluster = NULL;
79	preClusterUID = 0;
80	preClusterArray = NULL;
81
82	// delete dlist;
83	for (list<Tdlist *>::iterator i = myDList.begin(); i != myDList.end(); ++i) {
84	Tdlist dlist = i;
85	delete dlist;
86	dlist = NULL;
87	}
88	myDList.clear();
89
90
91	DebugOutput("TkJetFinder::initEvent() finished");
92	}
93
94	void AliTkKtJetFinder::preclusterTParticles(TClonesArray *particles) {
95	Char_t buffer[1024];
96	Int_t acceptedParticles = 0;
97	Int_t acceptedPreCluster = 0;
98
99	DetailedOutput("TkJetFinder::preclusterTParticle() called");
100
101	TStopwatch timer;
102	timer.Reset();
103
104	sprintf(buffer,"TkJetFinder::preclusterTParticle() - found %d particles",
105	particles->GetEntries());
106	DetailedOutput(buffer);
107
108	cout << "WARNING - TkJetFinder::preclusterTParticle()" << endl
109	<< " no preclustering implemented - each particle accepted!"
110	<< endl;
111
112	timer.Start();
113	// start the particle loop
114	TIterator *iter = particles->MakeIterator();
115	TParticle *particle;
116	TPreCluster *precluster;
117	while((particle = (TParticle *) iter->Next()) != NULL) {
118	if (isTParticleAccepted(particle)) {
119	acceptedParticles++;
120	precluster = new TPreCluster;
121	// to be kind of realistic use maseless particles
122	// - TPC doesn't measure energy and there is no hadron calorimeter...
123	precluster->id = preClusterUID;
124	precluster->E = TMath::Sqrt(particle->P()*particle->P());
125	precluster->px = particle->Px();
126	precluster->py = particle->Py();
127	precluster->pz = particle->Pz();
128	precluster->next = NULL;
129	if (true) {
130	// add precluster to list of preclusters
131	addPreCluster(precluster);
132	acceptedPreCluster++;
133	preClusterUID++;
134	} else {
135	delete precluster;
136	}
137	}
138	}
139
140	// build an array of pointers for direct access on the precluster...
141	preClusterArray = new TPreCluster*[preClusterUID];
142	precluster = firstPreCluster;
143	while (precluster) {
144	preClusterArray[precluster->id] = precluster;
145	precluster = precluster->next;
146	}
147
148	timer.Stop();
149
150	// if you want to be really sure what the precluster algorithm does...
151	//dumpPreClusters();
152	//dumpPreClusterArray();
153
154	// status information at the end...
155	sprintf(buffer,"TkJetFinder::preclusterTParticle() - used %d particles",
156	acceptedParticles);
157	DebugOutput(buffer);
158	sprintf(buffer,
159	"TkJetFinder::preclusterTParticle() - created %d preclusters",
160	acceptedPreCluster);
161	DebugOutput(buffer);
162
163	sprintf(buffer,
164	"TkJetFinder::preclusterTParticle() - timing CPU: %f, real %f",
165	timer.CpuTime(),timer.RealTime());
166	TimingOutput(buffer);
167
168	cout << "some stupid tests..." << endl;
169	Tdlist *p1 = new Tdlist;
170	Tdlist *p2 = new Tdlist;
171	Tdlist *p3 = new Tdlist;
172	p1->d = 3;
173	p2->d = 2;
174	p3->d = 1;
175	if (p1==p2) {
176	cout << "p1==p2 (" << p1->d << "==" << p2->d << ")" << endl;
177	} else {
178	cout << "p1!=p2 (" << p1->d << "!=" << p2->d << ")" << endl;
179	}
180	if (p1 < p2) {
181	cout << "p1 < p2 (" << p1->d << ">" << p2->d << ")" << endl;
182	} else {
183	cout << "p1 > p2 (" << p1->d << "<" << p2->d << ")" << endl;
184	}
185	priority_queue<Tdlist *> myHeap;
186	myHeap.push(p2);
187	myHeap.push(p1);
188	myHeap.push(p3);
189
190	cout << myHeap.size() << endl;
191	while (!myHeap.empty()) {
192	Tdlist *myEntry = myHeap.top();
193	myHeap.pop();
194	cout << "d = " << myEntry->d << endl;
195	delete myEntry;
196	}
197
198	DetailedOutput("TkJetFinder::preclusterTParticle() finished");
199	}
200
201	void AliTkKtJetFinder::findJets() {
202	Char_t buffer[1024];
203
204	TStopwatch timer;
205	timer.Reset();
206
207	DetailedOutput("TkJetFinder::findJets() called");
208
209	timer.Start();
210	// let's start...
211	buildNewDList();
212	// find the precluster/pair with lowest dmin...
213	//list<Tdlist *>::iterator i;
214	list<Tdlist *>::iterator min;
215
216	myDList.sort();
217	while(!myDList.empty()) {
218	min = myDList.begin();
219
220	// we have the entry with the lowest min...
221	//cout << (*min)->d << endl;
222	delete *min;
223	myDList.erase(min);
224	}
225
226
227	timer.Stop();
228
229	sprintf(buffer,
230	"TkJetFinder::findJets() - timing CPU: %f, real %f",
231	timer.CpuTime(),timer.RealTime());
232	TimingOutput(buffer);
233
234	DetailedOutput("TkJetFinder::findJets() finish");
235	}
236
237	void AliTkKtJetFinder::setPhiBins(Int_t nPhiBins,
238	Float_t fPhiMin, Float_t fPhiMax) {
239	setNPhiBins(nPhiBins);
240	setPhiMin(fPhiMin);
241	setPhiMax(fPhiMax);
242	}
243
244	void AliTkKtJetFinder::setThetaBins(Int_t nThetaBins,
245	Float_t fThetaMin, Float_t fThetaMax) {
246	setNThetaBins(nThetaBins);
247	setThetaMin(fThetaMin);
248	setThetaMax(fThetaMax);
249	}
250
251
252	void AliTkKtJetFinder::setDefaultOptions() {
253	// precluster options
254	setPhiBins((Int_t)(TMath::Pi()2./0.1), 0. ,2.TMath::Pi());
255	setThetaBins(10,(TMath::Pi()/2. - 0.5),(TMath::Pi()/2. + 0.5));
256
257	// finder options
258	setD(1.);
259	setDCut(-1.);
260
261	// debug options
262	setDebugLevel(0);
263	setDebugFilename("KtJetFinderDebug.root");
264
265	cout << "TkJetFinder::setDefaultOptions() - DEFAULT options set" << endl;
266	}
267
268	Bool_t AliTkKtJetFinder::isTParticleAccepted(TParticle *particle) {
269	// check if particle is accepted
270	// makes sense to write this into a own class, but now I'm lazy
271
272	// check if particle is stable -> a detectable particle
273	if (particle->GetStatusCode() != 1) {
274	return kFALSE;
275	}
276	// and now just for fun
277	TParticlePDG *partPDG;
278	partPDG = particle->GetPDG();
279	// if (partPDG->Charge() == 0) {
280	// cout << partPDG->GetName() << endl;
281	// return kFALSE;
282	// }
283
284	// default case: accept
285	return kTRUE;
286	}
287
288	void AliTkKtJetFinder::addPreCluster(TPreCluster *precluster) {
289	if (firstPreCluster == NULL) {
290	firstPreCluster = precluster;
291	lastPreCluster = precluster;
292	return;
293	}
294
295	lastPreCluster->next = precluster;
296	precluster->next = NULL;
297	lastPreCluster = precluster;
298	}
299
300	void AliTkKtJetFinder::deletePreCluster(Int_t UID) {
301	TPreCluster *p1 = NULL;
302	TPreCluster *p2 = firstPreCluster;
303	while (p2!=NULL) {
304	if (p2->id == UID) {
305	p1->next = p2->next;
306	delete p2;
307	break;
308	}
309	p1 = p2;
310	p2 = p2->next;
311	}
312	}
313
314	void AliTkKtJetFinder::dumpPreClusters() {
315	TPreCluster *pre = firstPreCluster;
316	while (pre) {
317	cout << "UID: " << pre->id
318	<< " E: " << pre->E
319	<< " px: " << pre->px
320	<< " py: " << pre->py
321	<< " pz: " << pre->pz
322	<< " next: " << pre->next
323	<< endl;
324	pre = pre->next;
325	}
326	}
327
328	void AliTkKtJetFinder::dumpPreClusterArray() {
329	for (Int_t i = 0; i < preClusterUID; i++) {
330	cout << "UID: " << preClusterArray[i]->id
331	<< " E: " << preClusterArray[i]->E
332	<< " px: " << preClusterArray[i]->px
333	<< " py: " << preClusterArray[i]->py
334	<< " pz: " << preClusterArray[i]->pz
335	<< " next: " << preClusterArray[i]->next
336	<< endl;
337	}
338	}
339
340	void AliTkKtJetFinder::addD(Tdlist *newD) {
341	myDList.push_back(newD);
342	}
343
344	void AliTkKtJetFinder::buildNewDList() {
345	Char_t buffer[1024];
346
347	// step 1 - add all preclusters to dlist
348	struct Tdlist *dlist;
349	for (Int_t i = 0; i < preClusterUID; i++) {
350	dlist = new Tdlist;
351	dlist->d = calcD((TPreCluster *)preClusterArray[i]);
352	dlist->id1 = preClusterArray[i]->id;
353	dlist->id2 = -1;
354	dlist->prev = NULL;
355	dlist->next = NULL;
356	addD(dlist);
357	}
358
359	sprintf(buffer,
360	"AliTkKtJetFinder::buildNewDList - added %d preclusters to new list",
361	(Int_t)myDList.size());
362	DetailedOutput(buffer);
363
364	// and now for all combinations...
365	for (Int_t i = 0; i < preClusterUID; i++) {
366	for (Int_t j = i + 1; j < preClusterUID; j++) {
367	if (i == j) {
368	// makes no sense...
369	continue;
370	}
371	dlist = new Tdlist;
372	dlist->d = calcD((TPreCluster *)preClusterArray[i],
373	(TPreCluster *)preClusterArray[j]);
374	dlist->id1 = preClusterArray[i]->id;
375	dlist->id2 = preClusterArray[j]->id;
376	dlist->prev = NULL;
377	dlist->next = NULL;
378	addD(dlist);
379	}
380	}
381
382	sprintf(buffer,
383	"AliTkKtJetFinder::buildNewDList - created list with %d entries",
384	(Int_t)myDList.size());
385	DetailedOutput(buffer);
386	}
387
388	Float_t AliTkKtJetFinder::calcD(TPreCluster p1, TPreCluster p2) {
389	if (!p1) {
390	// something is wrong - we need a precluster
391	return -9999;
392	}
393	if (!p2) {
394	// called with one particle only
395	// return pT^2
396	return (p1->pxp1->px + p1->pyp1->py);
397	}
398	// we have two particles...
399	// calculate the minimum pt...
400	Float_t pt1Sq = p1->pxp1->px + p1->pyp1->py;
401	Float_t pt2Sq = p2->pxp2->px + p2->pyp2->py;
402	Float_t ptMinSq;
403	if (pt1Sq < pt2Sq) {
404	ptMinSq = pt1Sq;
405	} else {
406	ptMinSq = pt2Sq;
407	}
408	// calculate rapidities - !move this into precluster structure!!!
409	Float_t y1 = 0.5 * TMath::Log((p1->E+p1->pz)/(p1->E-p1->pz));
410	Float_t y2 = 0.5 * TMath::Log((p2->E+p2->pz)/(p2->E-p2->pz));
411	// calculate phi - !move this into precluster structure!!!
412	Float_t phi1 = 0;
413	Float_t phi2 = 0;
414	// calculate ((y1-y2)^2 + (phi1-phi2)^2)/(D^2)
415	// !save D^2 instead of D to save time!!!
416	Float_t R = ((y1-y2)(y1-y2) + (phi1-phi2)(phi1-phi2))/(finder_D*finder_D);
417
418	return ptMinSq * R;
419	}
420
421	void AliTkKtJetFinder::setDebugLevel(Int_t nDebugLevel) {
422	debugLevel = nDebugLevel;
423
424	// description of debug levels
425	// levels are additive!!!
426	// 0 -> no debug output - only normal output (default)
427	// 1 -> normal debug output
428	// major function calls, summary of major results
429	// 2 -> timing output
430	// timing information for major steps
431	// 4 -> detailed results
432	// output of much more results
433	// 8 -> result histos
434	// histograms of event results for internal consitency checks
435	// 16 -> timing histos
436	// histograms with timing info - not cleared per event!
437	// 32 -> write debug histos to file
438	}
439
440	void AliTkKtJetFinder::DebugOutput(Char_t *output) {
441	if ((debugLevel & 1) == 1) {
442	cout << output << endl;
443	}
444	}
445
446	void AliTkKtJetFinder::TimingOutput(Char_t *output) {
447	if ((debugLevel & 2) == 2) {
448	cout << output << endl;
449	}
450	}
451
452	void AliTkKtJetFinder::DetailedOutput(Char_t *output) {
453	if ((debugLevel & 4) == 4) {
454	cout << output << endl;
455	}
456	}
457
458
459	ClassImp(AliTkKtJetFinder)