[u/mrichter/AliRoot.git] / JETAN / jetan2004 / AliTkJetTrigger.cxx

// $Id$

/****************************************************************************
 * AliTkJetTrigger.cxx                                                         *
 * Test version of jet-trigger, includes many analysis functions to tune    *
 * paramters                                                                *
 * Thorsten Kollegger <kollegge@ikf.physik.uni-frankfurt.de>                *
 ***************************************************************************/

#include <Riostream.h>
#include <TFile.h>
#include <TROOT.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TTree.h>
#include "AliTkChargedJetFinder.h"
#include "AliTkJetTriggerDecision.h"
#include "AliTkJetTriggerEvent.h"
#include "AliTkJetTrigger.h"

void AliTkJetTrigger::init() {
  // open file and create tree...
  if (getEvOutFilename() != NULL) {
    cout << "Writing events to " << getEvOutFilename() << endl;
    mEvOutFile = new TFile(getEvOutFilename(),"RECREATE");
    mEvOutEvent = new AliTkJetTriggerEvent();
    mEvOutTree = new TTree("jetTrigger","");
    mEvOutTree->Branch("event","AliTkJetTriggerEvent",&mEvOutEvent,32000,99);
  }
}

void AliTkJetTrigger::initEvent() {
  // can not be implemented, we need the particles...
}

void AliTkJetTrigger::initEvent(TClonesArray *particles, Int_t type) {
  setParticleType(type);
  initEvent(particles);
}

void AliTkJetTrigger::initEvent(TClonesArray *particles) {
  const UInt_t mParticlesDefaultSize = getExpectedParticleNumber();
  // make a copy of all accepted particles
  // just to be sure that they are not modified...

  clearParticles();
  if (!particles) {
    return;
  }
  mParticles = new TClonesArray("TParticle",mParticlesDefaultSize);
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  UInt_t pos = 0;
  while ((particle = (TParticle *) iter->Next()) != NULL) {
    if (isParticleAccepted(particle)) {
      if(particle->Pt()>fParticleMinPt)
	new ((*mParticles)[pos++]) TParticle(*particle);
    }
  }
  delete iter;
  createSeedPoints();
}

void AliTkJetTrigger::make() {
  // fill the matrix containing trigger decision for all different cone radia,
  // nParticles and PtThresholds...
  // loop over all seedpoints...
  cout << "generated " << seedPoints.size() << " seedpoints" << endl;
  for (list<AliTkEtaPhiVector>::const_iterator seedIter = seedPoints.begin();
       seedIter != seedPoints.end(); ++seedIter) {
    AliTkEtaPhiVector coneCenter(*seedIter);
    for (list<Float_t>::const_iterator radiusIter = mConeRadia.begin();
	 radiusIter != mConeRadia.end(); ++radiusIter) {
      // get particles in cone
      list<TParticle*> *myParticles = getParticlesInCone(coneCenter,
							 *radiusIter);
      // loop over all PtThresholds...
      for (list<Float_t>::const_iterator PtThreshold =mConePtThreshold.begin();
	   PtThreshold != mConePtThreshold.end(); ++PtThreshold) {
	// loop over all nParticles thresholds...
	for (list<UInt_t>::const_iterator nParticles = mConeNParticles.begin();
	     nParticles != mConeNParticles.end(); ++nParticles) {
	  // what's the decision???
	  if (decide(myParticles,*PtThreshold,*nParticles)) {
	    AliTkJetTriggerDecision *d = new AliTkJetTriggerDecision();
	    d->setConeRadius(*radiusIter);
	    d->setPtThreshold(*PtThreshold);
	    d->setNParticles(*nParticles);
	    mEvOutEvent->addDecision(d);
	    // new change, requires testing 
	    delete d;
	  }
	} // end of nParticle threshold loop
      } // end of PtThreshold loop
      myParticles->erase(myParticles->begin(),myParticles->end());
      delete myParticles;
    } // end of cone radia loop
  } // end of seed point loop
}

void AliTkJetTrigger::finishEvent() {
  // let's write out all triggered settings...
  // not finish yet...
  clearParticles();
  if (mEvOutFile) {
    mEvOutTree->Fill();
  }
  mEvOutEvent->clear();
  seedPoints.erase(seedPoints.begin(),seedPoints.end());
}

void AliTkJetTrigger::finish() {
  // write object containing settings
  // close file
  if (mEvOutFile) {
    if (!mEvOutFile->IsOpen()) { cerr << "WHAT?" << endl; }
    mEvOutFile->cd();
    mEvOutTree->Write();
    mEvOutFile->Close();
  }
}

void AliTkJetTrigger::setDefaults() {
  fParticleMinPt=0;
  fParticleMinPtSeed=1;
  setParticleType(3);
  addConeRadius(0.7);
  addConeRadius(0.5);
  addConeRadius(0.3);
  addConeRadius(0.2);
  addConeRadius(0.1);
  addConeNParticles(1);
  addConeNParticles(2);
  addConeNParticles(3);
  addConeNParticles(4);
  addConeNParticles(5);
  addConePtThreshold(1.);
  addConePtThreshold(2.);
  addConePtThreshold(3.);
  addConePtThreshold(4.);
  addConePtThreshold(5.);

  setEvOutFilename("$JF_DATADIR/trigger.evout.root");
  setExpectedParticleNumber(30000);
}

void AliTkJetTrigger::setParticleType(Int_t type) {
  mParticleType = type;
  // current used particles types:
  // 1 == TParticle from TPythia6
  // 2 == TParticle from THijing
  // 3 == TParticle from TPythia6 + ALICE acceptance approximation
  // 4 == TParticle from THijing  + ALICE acceptance approximation
}

void AliTkJetTrigger::addConeRadius(Float_t radius) {
  mConeRadia.push_back(radius);
}

void AliTkJetTrigger::addConeNParticles(UInt_t nParticles) {
  mConeNParticles.push_back(nParticles);
}

void AliTkJetTrigger::addConePtThreshold(Float_t ptThreshold) {
  mConePtThreshold.push_back(ptThreshold);
}

void AliTkJetTrigger::mixParticles(TClonesArray *particles, Int_t type) {
  // mixes particles to the existing event
  if (!mParticles) {
    return;
  }
  if (!particles) {
    return;
  }
  TParticle *particle = NULL;
  UInt_t pos = 0;
  pos = mParticles->GetLast()+1;
  TIterator *iter = particles->MakeIterator();
  while ((particle = (TParticle *) iter->Next()) != NULL) {
    if (isParticleAccepted(particle,type)) {
      new ((*mParticles)[pos++]) TParticle(*particle);
    }
  }
  delete iter;
  // we need to generate new seedPoints...
  seedPoints.erase(seedPoints.begin(),seedPoints.end());
  createSeedPoints();
}

void AliTkJetTrigger::createSeedPoints() {
  TParticle *particle;
  TIterator *iter = mParticles->MakeIterator();
  while ((particle = (TParticle *)iter->Next()) != NULL) {
    // particles are checked for acceptance during initEvent, mixParticles
    //    if (isParticleAccepted(particle)) {
      if (particle->Pt() > fParticleMinPtSeed) {
	AliTkEtaPhiVector v(particle->Eta(),particle->Phi());
	seedPoints.push_back(v);
      }
      //    }
  }
  delete iter;
}

Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle) {
  return isParticleAccepted(particle,getParticleType());
}

Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle, Int_t type) {
  switch (type) {
  case 1:
    return isParticleAcceptedPythia(particle);
    break;
  case 2:
    return isParticleAcceptedHijing(particle);
    break;
  case 3:
    return isParticleAcceptedPythiaAliceGeo(particle);
    break;
  case 4:
    return isParticleAcceptedHijingAliceGeo(particle);
    break;
  };
  return kFALSE;
}

Bool_t AliTkJetTrigger::isParticleAcceptedPythia(TParticle *particle) {
  // check for stable particle
  if (particle->GetStatusCode() != 1) {
    return kFALSE;
  }
  return kTRUE;
}

Bool_t AliTkJetTrigger::isParticleAcceptedPythiaAliceGeo(TParticle *particle) {
  // check if valid PYTHIA particle
  if (!isParticleAcceptedPythia(particle)) {
    return kFALSE;
  }
  // check if particle is in |eta| < 1
  if (TMath::Abs(particle->Eta()) > 1) {
    return kFALSE;
  }
  // check if charged particle
  TParticlePDG *partPDG = particle->GetPDG();
  if (partPDG->Charge() == 0) {
    return kFALSE;
  }
  return kTRUE;
}

Bool_t AliTkJetTrigger::isParticleAcceptedHijing(TParticle */*particle*/) {
  return kTRUE;
}

Bool_t AliTkJetTrigger::isParticleAcceptedHijingAliceGeo(TParticle *particle) {
  if (!isParticleAcceptedHijing(particle)) {
    return kFALSE;
  }
  // check if particle is in |eta| < 1
  if (TMath::Abs(particle->Eta()) > 1) {
    return kFALSE;
  }
  // check if charged particle
  TParticlePDG *partPDG = particle->GetPDG();
  if (partPDG->Charge() == 0) {
    return kFALSE;
  }
  return kTRUE;
}

void AliTkJetTrigger::clearParticles() {
  // delete my particle copy...
  if (mParticles) {
    mParticles->Delete();
    delete mParticles;
  }
  mParticles = NULL;
}

list<TParticle*> *AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
						   Float_t radius) {
  return getParticlesInCone(center,radius,mParticles);
}

list<TParticle*> *AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
						   Float_t radius,
						   TClonesArray *particles) {
  // returns a list with all particles which are in a cone with radius around
  // center
  // one can easily speed up this part by saving the last cone+parameters and
  // calculate the new one from this one if there is some overlap
  list<TParticle*> *myParticles = new list<TParticle*>;
  Float_t radiusSq = radius*radius;
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  while ((particle = (TParticle *)iter->Next()) != NULL) {
    AliTkEtaPhiVector pos(particle->Eta(),particle->Phi());
    if (center.diffSq(pos) < radiusSq) {
      myParticles->push_back(particle);
    }
  }
  delete iter;
  return myParticles;
}

Bool_t AliTkJetTrigger::decide(list<TParticle*> *particles, Float_t ptThreshold, 
			    UInt_t nParticles) {
  UInt_t n = 0;
  for (list<TParticle*>::const_iterator iter = particles->begin();
       iter != particles->end(); ++iter) {
    if ((*iter)->Pt() > ptThreshold) {
      n++;
    }
  }
  if (n >= nParticles) {
    return kTRUE;
  }
  return kFALSE;
}

void AliTkJetTrigger::setEvOutFilename(Char_t *filename) {
  if (!mEvOutName) {
    mEvOutName = new Char_t[4096];
  }
  strcpy(mEvOutName,filename);
}


ClassImp(AliTkJetTrigger)
Commit	Line	Data
b9a6a391	1	// $Id$
	2
	3	/****************************************************************************
	4	* AliTkJetTrigger.cxx *
	5	* Test version of jet-trigger, includes many analysis functions to tune *
	6	* paramters *
	7	* Thorsten Kollegger <kollegge@ikf.physik.uni-frankfurt.de> *
	8	***************************************************************************/
	9
	10	#include <Riostream.h>
	11	#include <TFile.h>
	12	#include <TROOT.h>
	13	#include <TParticle.h>
	14	#include <TClonesArray.h>
	15	#include <TTree.h>
	16	#include "AliTkChargedJetFinder.h"
	17	#include "AliTkJetTriggerDecision.h"
	18	#include "AliTkJetTriggerEvent.h"
	19	#include "AliTkJetTrigger.h"
	20
	21	void AliTkJetTrigger::init() {
	22	// open file and create tree...
	23	if (getEvOutFilename() != NULL) {
	24	cout << "Writing events to " << getEvOutFilename() << endl;
	25	mEvOutFile = new TFile(getEvOutFilename(),"RECREATE");
	26	mEvOutEvent = new AliTkJetTriggerEvent();
	27	mEvOutTree = new TTree("jetTrigger","");
	28	mEvOutTree->Branch("event","AliTkJetTriggerEvent",&mEvOutEvent,32000,99);
	29	}
	30	}
	31
	32	void AliTkJetTrigger::initEvent() {
	33	// can not be implemented, we need the particles...
	34	}
	35
	36	void AliTkJetTrigger::initEvent(TClonesArray *particles, Int_t type) {
	37	setParticleType(type);
	38	initEvent(particles);
	39	}
	40
	41	void AliTkJetTrigger::initEvent(TClonesArray *particles) {
	42	const UInt_t mParticlesDefaultSize = getExpectedParticleNumber();
	43	// make a copy of all accepted particles
	44	// just to be sure that they are not modified...
	45
	46	clearParticles();
	47	if (!particles) {
	48	return;
	49	}
	50	mParticles = new TClonesArray("TParticle",mParticlesDefaultSize);
	51	TIterator *iter = particles->MakeIterator();
	52	TParticle *particle;
	53	UInt_t pos = 0;
	54	while ((particle = (TParticle *) iter->Next()) != NULL) {
	55	if (isParticleAccepted(particle)) {
	56	if(particle->Pt()>fParticleMinPt)
	57	new ((mParticles)[pos++]) TParticle(particle);
	58	}
	59	}
	60	delete iter;
	61	createSeedPoints();
	62	}
	63
	64	void AliTkJetTrigger::make() {
65	// fill the matrix containing trigger decision for all different cone radia,
66	// nParticles and PtThresholds...
67	// loop over all seedpoints...
68	cout << "generated " << seedPoints.size() << " seedpoints" << endl;
69	for (list<AliTkEtaPhiVector>::const_iterator seedIter = seedPoints.begin();
70	seedIter != seedPoints.end(); ++seedIter) {
71	AliTkEtaPhiVector coneCenter(*seedIter);
72	for (list<Float_t>::const_iterator radiusIter = mConeRadia.begin();
73	radiusIter != mConeRadia.end(); ++radiusIter) {
74	// get particles in cone
75	list<TParticle> myParticles = getParticlesInCone(coneCenter,
76	*radiusIter);
77	// loop over all PtThresholds...
78	for (list<Float_t>::const_iterator PtThreshold =mConePtThreshold.begin();
79	PtThreshold != mConePtThreshold.end(); ++PtThreshold) {
80	// loop over all nParticles thresholds...
81	for (list<UInt_t>::const_iterator nParticles = mConeNParticles.begin();
82	nParticles != mConeNParticles.end(); ++nParticles) {
83	// what's the decision???
84	if (decide(myParticles,PtThreshold,nParticles)) {
85	AliTkJetTriggerDecision *d = new AliTkJetTriggerDecision();
86	d->setConeRadius(*radiusIter);
87	d->setPtThreshold(*PtThreshold);
88	d->setNParticles(*nParticles);
89	mEvOutEvent->addDecision(d);
90	// new change, requires testing
91	delete d;
92	}
93	} // end of nParticle threshold loop
94	} // end of PtThreshold loop
95	myParticles->erase(myParticles->begin(),myParticles->end());
96	delete myParticles;
97	} // end of cone radia loop
98	} // end of seed point loop
99	}
100
101	void AliTkJetTrigger::finishEvent() {
102	// let's write out all triggered settings...
103	// not finish yet...
104	clearParticles();
105	if (mEvOutFile) {
106	mEvOutTree->Fill();
107	}
108	mEvOutEvent->clear();
109	seedPoints.erase(seedPoints.begin(),seedPoints.end());
110	}
111
112	void AliTkJetTrigger::finish() {
113	// write object containing settings
114	// close file
115	if (mEvOutFile) {
116	if (!mEvOutFile->IsOpen()) { cerr << "WHAT?" << endl; }
117	mEvOutFile->cd();
118	mEvOutTree->Write();
119	mEvOutFile->Close();
120	}
121	}
122
123	void AliTkJetTrigger::setDefaults() {
124	fParticleMinPt=0;
125	fParticleMinPtSeed=1;
126	setParticleType(3);
127	addConeRadius(0.7);
128	addConeRadius(0.5);
129	addConeRadius(0.3);
130	addConeRadius(0.2);
131	addConeRadius(0.1);
132	addConeNParticles(1);
133	addConeNParticles(2);
134	addConeNParticles(3);
135	addConeNParticles(4);
136	addConeNParticles(5);
137	addConePtThreshold(1.);
138	addConePtThreshold(2.);
139	addConePtThreshold(3.);
140	addConePtThreshold(4.);
141	addConePtThreshold(5.);
142
143	setEvOutFilename("$JF_DATADIR/trigger.evout.root");
144	setExpectedParticleNumber(30000);
145	}
146
147	void AliTkJetTrigger::setParticleType(Int_t type) {
148	mParticleType = type;
149	// current used particles types:
150	// 1 == TParticle from TPythia6
151	// 2 == TParticle from THijing
152	// 3 == TParticle from TPythia6 + ALICE acceptance approximation
153	// 4 == TParticle from THijing + ALICE acceptance approximation
154	}
155
156	void AliTkJetTrigger::addConeRadius(Float_t radius) {
157	mConeRadia.push_back(radius);
158	}
159
160	void AliTkJetTrigger::addConeNParticles(UInt_t nParticles) {
161	mConeNParticles.push_back(nParticles);
162	}
163
164	void AliTkJetTrigger::addConePtThreshold(Float_t ptThreshold) {
165	mConePtThreshold.push_back(ptThreshold);
166	}
167
168	void AliTkJetTrigger::mixParticles(TClonesArray *particles, Int_t type) {
169	// mixes particles to the existing event
170	if (!mParticles) {
171	return;
172	}
173	if (!particles) {
174	return;
175	}
176	TParticle *particle = NULL;
177	UInt_t pos = 0;
178	pos = mParticles->GetLast()+1;
179	TIterator *iter = particles->MakeIterator();
180	while ((particle = (TParticle *) iter->Next()) != NULL) {
181	if (isParticleAccepted(particle,type)) {
182	new ((mParticles)[pos++]) TParticle(particle);
183	}
184	}
185	delete iter;
186	// we need to generate new seedPoints...
187	seedPoints.erase(seedPoints.begin(),seedPoints.end());
188	createSeedPoints();
189	}
190
191	void AliTkJetTrigger::createSeedPoints() {
192	TParticle *particle;
193	TIterator *iter = mParticles->MakeIterator();
194	while ((particle = (TParticle *)iter->Next()) != NULL) {
195	// particles are checked for acceptance during initEvent, mixParticles
196	// if (isParticleAccepted(particle)) {
197	if (particle->Pt() > fParticleMinPtSeed) {
198	AliTkEtaPhiVector v(particle->Eta(),particle->Phi());
199	seedPoints.push_back(v);
200	}
201	// }
202	}
203	delete iter;
204	}
205
206	Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle) {
207	return isParticleAccepted(particle,getParticleType());
208	}
209
210	Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle, Int_t type) {
211	switch (type) {
212	case 1:
213	return isParticleAcceptedPythia(particle);
214	break;
215	case 2:
216	return isParticleAcceptedHijing(particle);
217	break;
218	case 3:
219	return isParticleAcceptedPythiaAliceGeo(particle);
220	break;
221	case 4:
222	return isParticleAcceptedHijingAliceGeo(particle);
223	break;
224	};
225	return kFALSE;
226	}
227
228	Bool_t AliTkJetTrigger::isParticleAcceptedPythia(TParticle *particle) {
229	// check for stable particle
230	if (particle->GetStatusCode() != 1) {
231	return kFALSE;
232	}
233	return kTRUE;
234	}
235
236	Bool_t AliTkJetTrigger::isParticleAcceptedPythiaAliceGeo(TParticle *particle) {
237	// check if valid PYTHIA particle
238	if (!isParticleAcceptedPythia(particle)) {
239	return kFALSE;
240	}
241	// check if particle is in \|eta\| < 1
242	if (TMath::Abs(particle->Eta()) > 1) {
243	return kFALSE;
244	}
245	// check if charged particle
246	TParticlePDG *partPDG = particle->GetPDG();
247	if (partPDG->Charge() == 0) {
248	return kFALSE;
249	}
250	return kTRUE;
251	}
252
253	Bool_t AliTkJetTrigger::isParticleAcceptedHijing(TParticle /particle*/) {
254	return kTRUE;
255	}
256
257	Bool_t AliTkJetTrigger::isParticleAcceptedHijingAliceGeo(TParticle *particle) {
258	if (!isParticleAcceptedHijing(particle)) {
259	return kFALSE;
260	}
261	// check if particle is in \|eta\| < 1
262	if (TMath::Abs(particle->Eta()) > 1) {
263	return kFALSE;
264	}
265	// check if charged particle
266	TParticlePDG *partPDG = particle->GetPDG();
267	if (partPDG->Charge() == 0) {
268	return kFALSE;
269	}
270	return kTRUE;
271	}
272
273	void AliTkJetTrigger::clearParticles() {
274	// delete my particle copy...
275	if (mParticles) {
276	mParticles->Delete();
277	delete mParticles;
278	}
279	mParticles = NULL;
280	}
281
282	list<TParticle> AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
283	Float_t radius) {
284	return getParticlesInCone(center,radius,mParticles);
285	}
286
287	list<TParticle> AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
288	Float_t radius,
289	TClonesArray *particles) {
290	// returns a list with all particles which are in a cone with radius around
291	// center
292	// one can easily speed up this part by saving the last cone+parameters and
293	// calculate the new one from this one if there is some overlap
294	list<TParticle> myParticles = new list<TParticle*>;
295	Float_t radiusSq = radius*radius;
296	TIterator *iter = particles->MakeIterator();
297	TParticle *particle;
298	while ((particle = (TParticle *)iter->Next()) != NULL) {
299	AliTkEtaPhiVector pos(particle->Eta(),particle->Phi());
300	if (center.diffSq(pos) < radiusSq) {
301	myParticles->push_back(particle);
302	}
303	}
304	delete iter;
305	return myParticles;
306	}
307
308	Bool_t AliTkJetTrigger::decide(list<TParticle> particles, Float_t ptThreshold,
309	UInt_t nParticles) {
310	UInt_t n = 0;
311	for (list<TParticle*>::const_iterator iter = particles->begin();
312	iter != particles->end(); ++iter) {
313	if ((*iter)->Pt() > ptThreshold) {
314	n++;
315	}
316	}
317	if (n >= nParticles) {
318	return kTRUE;
319	}
320	return kFALSE;
321	}
322
323	void AliTkJetTrigger::setEvOutFilename(Char_t *filename) {
324	if (!mEvOutName) {
325	mEvOutName = new Char_t[4096];
326	}
327	strcpy(mEvOutName,filename);
328	}
329
330
331	ClassImp(AliTkJetTrigger)