[u/mrichter/AliRoot.git] / PHOS / AliPHOSJetFinder.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.                  *
 **************************************************************************/

//_________________________________________________________________________
// C++ version of UA2 and/or Lund jet finding algorithm
// UA1 jet algorithm from LUND JETSET (LUCELL)
// Find jets at the level of no detector and Digits.
// Needs modifications. 
//*-- Author : D.Peressounko after UA1 coll. etc
//////////////////////////////////////////////////////////////////////////////

/* $Id$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.8  2005/05/28 14:19:04  schutz
 * Compilation warnings fixed by T.P.
 *
 */

// --- ROOT system ---
#include "TClonesArray.h"
//      #include "TIter.h"
#include "TParticle.h"
#include "TTask.h"
// --- Standard library ---

// --- AliRoot header files ---
#include "AliLog.h"
#include "AliPHOSJet.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSDigit.h"
#include "AliPHOSGetter.h"
#include "AliPHOSJetFinder.h"
#include "AliPHOSDigitizer.h"

ClassImp(AliPHOSJetFinder)


//____________________________________________________________________________ 
AliPHOSJetFinder::AliPHOSJetFinder():
  TNamed("AliPHOSJetFinder",""),
  fNJets(0),
  fStatusCode(-999),
  fMode(0),
  fConeRad(1.),
  fMaxConeMove(0.15),
  fMinConeMove(0.05),
  fEtSeed(4.),
  fEtMin(5.),
  fPrecBg(0.00035),
  fSimGain(0.),
  fSimPedestal(0.),
  fParticles(0),
  fJets(0)
{
  //Initialize jet parameters
}

//____________________________________________________________________________ 
AliPHOSJetFinder::AliPHOSJetFinder(const AliPHOSJetFinder & jet) : 
  TNamed(jet),
  fNJets(0),
  fStatusCode(-999),
  fMode(0),
  fConeRad(1.),
  fMaxConeMove(0.15),
  fMinConeMove(0.05),
  fEtSeed(4.),
  fEtMin(5.),
  fPrecBg(0.00035),
  fSimGain(0.),
  fSimPedestal(0.),
  fParticles(0),
  fJets(0)
{
  // copy ctor: no implementation yet
  Fatal("cpy ctor", "not implemented");
}

//____________________________________________________________________________ 
  AliPHOSJetFinder::~AliPHOSJetFinder()
{
  //dtor
  if(fParticles){
    delete fParticles ;
    fParticles = 0 ;
  }
  if(fJets){
    delete fJets ;
    fJets = 0 ;
  } 
}

//____________________________________________________________________________ 
void  AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist) 
{
  //Find jets in the case without detector.
  TIter next(plist) ;

  TIter nextJet(jetslist) ;

  fNJets = 0 ;
  TParticle * p ;
  AliPHOSJet * jet ;
  //In this cicle we find number of jets and define approx. their directions
  //note, that we do not really add particles to jet (index =-1)
  while((p=static_cast<TParticle*>(next()))){
    if(fStatusCode==-999 || p->GetStatusCode()==fStatusCode){
    if(p->Energy() >= fEtSeed){ //Energetic enough
      //cout << "p " << p->Energy() << endl ;
      //cout << "Status "<<fStatusCode<<" "<<p->GetName()<<" " << p->Energy() << " "<<p->Eta()<< " "<<p->Phi()<<endl ;
      Bool_t startnew = kTRUE ;
      //Do not start new jet if part of older jet
      nextJet.Reset() ;
      while((jet=static_cast<AliPHOSJet*>(nextJet()))){
	//jet->Print() ;
	if(jet->AcceptConeDeviation(p)){
	  startnew = kFALSE ;
	  //cout << "false" << endl ;
	  break ;
	}
      }
      if(startnew){
	//cout << "new " << endl ;
	jet = new AliPHOSJet() ;
	jetslist->Add(jet) ;
	//	jet = static_cast<AliPHOSJet*>(jetslist->Last()) ;
	jet->SetConeRadius(fConeRad) ;
	jet->SetMaxConeMove(fMaxConeMove) ;
	//	jet->SetMinConeMove(fMinConeMove) ;
	jet->AddParticle(p,-1) ;
	fNJets++;
      }
    }
    while((jet=static_cast<AliPHOSJet*>(nextJet()))){
      if(jet->AcceptConeDeviation(p))
	jet->AddParticle(p,-1) ; //Just recalculate direction of jet
    }
    }
  }
    
  //now calculate directions of jets using collected information
  nextJet.Reset() ;
  while((jet=static_cast<AliPHOSJet*>(nextJet()))){
    jet->CalculateAll() ;
    if(jet->Energy() < fEtMin){
      jetslist->Remove(jet) ;
      delete jet ;
    }
  }

  jetslist->Compress() ;
  //And finally, really add particles to jets
  for(Int_t iPart=0; iPart<plist->GetEntries();iPart++){
    p=static_cast<TParticle*>(plist->At(iPart)) ;
    if(fStatusCode == -999 || p->GetStatusCode()==fStatusCode){
    Double_t dist = 999999. ; //big distance
    Int_t iJet = -1 ;
    for(Int_t i=0; i<jetslist->GetEntriesFast();i++){
      jet=static_cast<AliPHOSJet*>(jetslist->At(i)) ;
      if(jet->IsInCone(p)){
	Double_t cdist = jet->DistanceToJet(p); 
	if(cdist < dist){
	  dist = cdist ;
	  iJet = i ;
	}
      }
    }
    if(iJet>-1)
      (static_cast<AliPHOSJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
  }
  }
  
  //Calculate jet parameters 
  nextJet.Reset() ;
  while((jet=static_cast<AliPHOSJet*>(nextJet()))){
    jet->CalculateAll() ;
  }

}
//____________________________________________________________________________ 
void AliPHOSJetFinder::FindJetsFromDigits(const TClonesArray * digits, TObjArray * jets){
  //Find jets in the case witht detector at the level of digits.
  if(digits->GetEntries()==0){
    AliError(Form("No entries in digits list \n")) ;
    return ;
  }


  TClonesArray * copyDigits = new TClonesArray(*digits) ;

  //Remove CPV digits if any
  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
  Int_t iDigit ;
  AliPHOSDigit * digit ;
  for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
    digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
    if(!geom->IsInEMC(digit->GetId()))
      copyDigits->RemoveAt(iDigit) ;
    else
      break ;
  }
  copyDigits->Compress() ;

  Double_t totalEnergy = 0 ;
  Float_t * energy = new Float_t[copyDigits->GetEntries()] ;
  //calculate average energy of digits
  //fill array of energies
  for(iDigit=0;iDigit<copyDigits->GetEntries();iDigit++){
    digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
    energy[iDigit] = Calibrate(digit) ;
    totalEnergy+=energy[iDigit] ;
  }
  
  //Sort digits in decreasing energy.
  Int_t * index = new Int_t[copyDigits->GetEntries()] ;
  TMath::Sort(copyDigits->GetEntries(),energy,index) ;
  
  Double_t eAverage = totalEnergy/copyDigits->GetEntries()  ;
  //remove digits below average energy
  for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
    digit=static_cast<AliPHOSDigit *>(copyDigits->At(index[iDigit])) ;
    if(energy[index[iDigit]] < eAverage)
      copyDigits->RemoveAt(iDigit) ;
    else
      break ;
  }
  
  
  AliPHOSJet * jet ;
  Int_t iIter = 0 ;
  while(iIter < 10){//less than 10 iterations
    
    //while digits above seed
    for(Int_t ind=0;ind<copyDigits->GetEntriesFast();ind++){
      digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
      if(energy[index[ind]] > fEtSeed && digit){ //start new jet      
 	jet = new AliPHOSJet() ;
	Double_t e,eta,phi ;
	CalculateEEtaPhi(digit,e,eta,phi) ;
	jet->AddDigit(e,eta,phi,-1) ;  
	//loop over left digits
	for(Int_t iDigit = 0 ; iDigit < copyDigits->GetEntries() ; iDigit++){
	  if(iDigit!= ind){ //first digit already in jet
	    digit = static_cast<AliPHOSDigit *>(copyDigits->At(iDigit));
	    CalculateEEtaPhi(digit,e,eta,phi) ;	    
	    if(jet->IsInCone(eta,phi) && //is cell in cone
	       jet->AcceptConeDeviation(e,eta,phi)){//if cone does not move too much	      
	      jet->AddDigit(e,eta,phi,-1) ;  //accept new direction
	    }
	  }
	}//end of loop over cells
	
	//accept  all anused cells incide cone
	//note, that digits might be returned as anused later
	for(Int_t icell = 0 ; icell < copyDigits->GetEntries() ; icell++){
	  digit = static_cast<AliPHOSDigit *>(copyDigits->At(icell));
	  if(jet->IsInCone(eta,phi)){ //is cell in cone
	    CalculateEEtaPhi(digit,e,eta,phi) ;
	    jet->AddDigit(e,eta,phi,digit->GetIndexInList()) ; 
	  }
	}
	
	//Accept Jet with Et > Et_min and remove all belonging digits
	if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
	  Int_t nIndxs ;
	  const Int_t * indxs = jet->Indexs(nIndxs) ;
	  for(Int_t i=0;i<nIndxs;i++){
	    copyDigits->RemoveAt(indxs[i]) ;
	  }
	  jet->CalculateAll() ;
	  jets->AddAt(jet,fNJets++);
	}
	else{ //remove jet and do not touch digits
	  delete jet ;
	}
      }
      else{ 
	if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
	  break ;
	}
      }
      
      iIter++ ;
      //calculate new energy of backrgound
      Double_t oldTotalEnergy = totalEnergy ;
      totalEnergy = 0 ;
      for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
	digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
	if(digit)
	  totalEnergy+=energy[i] ;
      }
      if(!fMode || (oldTotalEnergy != 0) && 
	 (TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg))
	break ;	
    }
  }
  delete [] energy;
  delete [] index;
  copyDigits->Delete() ;
  
}
//____________________________________________________________________________ 
Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){ 
  //Both simulated and raw digits are already calibrated
  
  return digit->GetEnergy() ;        
  
  //  }  
}
//____________________________________________________________________________ 
void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
  //Calculate direction of the jet
  e=Calibrate(d) ;
  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
  TVector3 pos ;
  geom->RelPosInAlice(d->GetId(), pos) ;
  eta = pos.Eta() ;
  phi = pos.Phi() ;
}
//____________________________________________________________________________ 
void AliPHOSJetFinder::Print(const Option_t *) const {	
  //Print parameters of the found jet
  printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
  printf(" Jets found .........%d \n",fNJets) ;
  printf(" Seed energy cut ....%f \n",fEtSeed) ;
  printf(" Cone radius ........%f \n",fConeRad) ;
  printf(" Minimal cone move ..%f \n",fMinConeMove) ;
  printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
  printf("------------------------------------------------- \n") ;
}
Commit	Line	Data
5240525e	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	//_________________________________________________________________________
	17	// C++ version of UA2 and/or Lund jet finding algorithm
0bc3b8ed	18	// UA1 jet algorithm from LUND JETSET (LUCELL)
	19	// Find jets at the level of no detector and Digits.
	20	// Needs modifications.
5240525e	21	//*-- Author : D.Peressounko after UA1 coll. etc
	22	//////////////////////////////////////////////////////////////////////////////
	23
702ab87e	24	/* $Id$ */
	25
	26	/* History of cvs commits:
	27	*
	28	* $Log$
3f7dbdb7	29	* Revision 1.8 2005/05/28 14:19:04 schutz
	30	* Compilation warnings fixed by T.P.
	31	*
702ab87e	32	*/
702ab87e	33
5240525e	34	// --- ROOT system ---
	35	#include "TClonesArray.h"
	36	// #include "TIter.h"
	37	#include "TParticle.h"
	38	#include "TTask.h"
	39	// --- Standard library ---
	40
	41	// --- AliRoot header files ---
351dd634	42	#include "AliLog.h"
5240525e	43	#include "AliPHOSJet.h"
e957fea8	44	#include "AliPHOSGeometry.h"
5240525e	45	#include "AliPHOSDigit.h"
	46	#include "AliPHOSGetter.h"
	47	#include "AliPHOSJetFinder.h"
	48	#include "AliPHOSDigitizer.h"
	49
	50	ClassImp(AliPHOSJetFinder)
	51
	52
	53	//____________________________________________________________________________
3f7dbdb7	54	AliPHOSJetFinder::AliPHOSJetFinder():
	55	TNamed("AliPHOSJetFinder",""),
	56	fNJets(0),
	57	fStatusCode(-999),
	58	fMode(0),
	59	fConeRad(1.),
	60	fMaxConeMove(0.15),
	61	fMinConeMove(0.05),
	62	fEtSeed(4.),
	63	fEtMin(5.),
	64	fPrecBg(0.00035),
	65	fSimGain(0.),
	66	fSimPedestal(0.),
	67	fParticles(0),
	68	fJets(0)
5240525e	69	{
0bc3b8ed	70	//Initialize jet parameters
3f7dbdb7	71	}
5240525e	72
3f7dbdb7	73	//____________________________________________________________________________
	74	AliPHOSJetFinder::AliPHOSJetFinder(const AliPHOSJetFinder & jet) :
	75	TNamed(jet),
	76	fNJets(0),
	77	fStatusCode(-999),
	78	fMode(0),
	79	fConeRad(1.),
	80	fMaxConeMove(0.15),
	81	fMinConeMove(0.05),
	82	fEtSeed(4.),
	83	fEtMin(5.),
	84	fPrecBg(0.00035),
	85	fSimGain(0.),
	86	fSimPedestal(0.),
	87	fParticles(0),
	88	fJets(0)
	89	{
	90	// copy ctor: no implementation yet
	91	Fatal("cpy ctor", "not implemented");
5240525e	92	}
	93
	94	//____________________________________________________________________________
	95	AliPHOSJetFinder::~AliPHOSJetFinder()
	96	{
0bc3b8ed	97	//dtor
5240525e	98	if(fParticles){
	99	delete fParticles ;
	100	fParticles = 0 ;
	101	}
	102	if(fJets){
	103	delete fJets ;
	104	fJets = 0 ;
	105	}
	106	}
	107
	108	//____________________________________________________________________________
	109	void AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist)
	110	{
0bc3b8ed	111	//Find jets in the case without detector.
5240525e	112	TIter next(plist) ;
	113
	114	TIter nextJet(jetslist) ;
	115
	116	fNJets = 0 ;
	117	TParticle * p ;
	118	AliPHOSJet * jet ;
	119	//In this cicle we find number of jets and define approx. their directions
	120	//note, that we do not really add particles to jet (index =-1)
	121	while((p=static_cast<TParticle*>(next()))){
	122	if(fStatusCode==-999 \|\| p->GetStatusCode()==fStatusCode){
	123	if(p->Energy() >= fEtSeed){ //Energetic enough
	124	//cout << "p " << p->Energy() << endl ;
	125	//cout << "Status "<<fStatusCode<<" "<<p->GetName()<<" " << p->Energy() << " "<<p->Eta()<< " "<<p->Phi()<<endl ;
	126	Bool_t startnew = kTRUE ;
	127	//Do not start new jet if part of older jet
	128	nextJet.Reset() ;
	129	while((jet=static_cast<AliPHOSJet*>(nextJet()))){
	130	//jet->Print() ;
	131	if(jet->AcceptConeDeviation(p)){
	132	startnew = kFALSE ;
	133	//cout << "false" << endl ;
	134	break ;
	135	}
	136	}
	137	if(startnew){
	138	//cout << "new " << endl ;
	139	jet = new AliPHOSJet() ;
	140	jetslist->Add(jet) ;
	141	// jet = static_cast<AliPHOSJet*>(jetslist->Last()) ;
	142	jet->SetConeRadius(fConeRad) ;
	143	jet->SetMaxConeMove(fMaxConeMove) ;
	144	// jet->SetMinConeMove(fMinConeMove) ;
	145	jet->AddParticle(p,-1) ;
	146	fNJets++;
	147	}
	148	}
	149	while((jet=static_cast<AliPHOSJet*>(nextJet()))){
	150	if(jet->AcceptConeDeviation(p))
	151	jet->AddParticle(p,-1) ; //Just recalculate direction of jet
	152	}
	153	}
	154	}
	155
	156	//now calculate directions of jets using collected information
	157	nextJet.Reset() ;
	158	while((jet=static_cast<AliPHOSJet*>(nextJet()))){
	159	jet->CalculateAll() ;
	160	if(jet->Energy() < fEtMin){
	161	jetslist->Remove(jet) ;
	162	delete jet ;
	163	}
	164	}
	165
	166	jetslist->Compress() ;
	167	//And finally, really add particles to jets
	168	for(Int_t iPart=0; iPart<plist->GetEntries();iPart++){
	169	p=static_cast<TParticle*>(plist->At(iPart)) ;
	170	if(fStatusCode == -999 \|\| p->GetStatusCode()==fStatusCode){
	171	Double_t dist = 999999. ; //big distance
	172	Int_t iJet = -1 ;
	173	for(Int_t i=0; i<jetslist->GetEntriesFast();i++){
	174	jet=static_cast<AliPHOSJet*>(jetslist->At(i)) ;
	175	if(jet->IsInCone(p)){
176	Double_t cdist = jet->DistanceToJet(p);
177	if(cdist < dist){
178	dist = cdist ;
179	iJet = i ;
180	}
181	}
182	}
183	if(iJet>-1)
184	(static_cast<AliPHOSJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
185	}
186	}
187
188	//Calculate jet parameters
189	nextJet.Reset() ;
190	while((jet=static_cast<AliPHOSJet*>(nextJet()))){
191	jet->CalculateAll() ;
192	}
193
194	}
195	//____________________________________________________________________________
196	void AliPHOSJetFinder::FindJetsFromDigits(const TClonesArray * digits, TObjArray * jets){
0bc3b8ed	197	//Find jets in the case witht detector at the level of digits.
5240525e	198	if(digits->GetEntries()==0){
351dd634	199	AliError(Form("No entries in digits list \n")) ;
5240525e	200	return ;
	201	}
	202
	203
	204	TClonesArray * copyDigits = new TClonesArray(*digits) ;
	205
	206	//Remove CPV digits if any
	207	AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
	208	Int_t iDigit ;
	209	AliPHOSDigit * digit ;
	210	for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
	211	digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
	212	if(!geom->IsInEMC(digit->GetId()))
	213	copyDigits->RemoveAt(iDigit) ;
	214	else
	215	break ;
	216	}
	217	copyDigits->Compress() ;
	218
	219	Double_t totalEnergy = 0 ;
2fa01244	220	Float_t * energy = new Float_t[copyDigits->GetEntries()] ;
5240525e	221	//calculate average energy of digits
	222	//fill array of energies
	223	for(iDigit=0;iDigit<copyDigits->GetEntries();iDigit++){
	224	digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
	225	energy[iDigit] = Calibrate(digit) ;
	226	totalEnergy+=energy[iDigit] ;
	227	}
	228
	229	//Sort digits in decreasing energy.
2fa01244	230	Int_t * index = new Int_t[copyDigits->GetEntries()] ;
5240525e	231	TMath::Sort(copyDigits->GetEntries(),energy,index) ;
	232
	233	Double_t eAverage = totalEnergy/copyDigits->GetEntries() ;
	234	//remove digits below average energy
	235	for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
	236	digit=static_cast<AliPHOSDigit *>(copyDigits->At(index[iDigit])) ;
	237	if(energy[index[iDigit]] < eAverage)
	238	copyDigits->RemoveAt(iDigit) ;
	239	else
	240	break ;
	241	}
	242
	243
	244	AliPHOSJet * jet ;
	245	Int_t iIter = 0 ;
	246	while(iIter < 10){//less than 10 iterations
	247
	248	//while digits above seed
	249	for(Int_t ind=0;ind<copyDigits->GetEntriesFast();ind++){
	250	digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
	251	if(energy[index[ind]] > fEtSeed && digit){ //start new jet
	252	jet = new AliPHOSJet() ;
	253	Double_t e,eta,phi ;
	254	CalculateEEtaPhi(digit,e,eta,phi) ;
	255	jet->AddDigit(e,eta,phi,-1) ;
	256	//loop over left digits
	257	for(Int_t iDigit = 0 ; iDigit < copyDigits->GetEntries() ; iDigit++){
	258	if(iDigit!= ind){ //first digit already in jet
	259	digit = static_cast<AliPHOSDigit *>(copyDigits->At(iDigit));
	260	CalculateEEtaPhi(digit,e,eta,phi) ;
90cceaf6	261	if(jet->IsInCone(eta,phi) && //is cell in cone
5240525e	262	jet->AcceptConeDeviation(e,eta,phi)){//if cone does not move too much
	263	jet->AddDigit(e,eta,phi,-1) ; //accept new direction
	264	}
	265	}
	266	}//end of loop over cells
	267
	268	//accept all anused cells incide cone
	269	//note, that digits might be returned as anused later
	270	for(Int_t icell = 0 ; icell < copyDigits->GetEntries() ; icell++){
	271	digit = static_cast<AliPHOSDigit *>(copyDigits->At(icell));
90cceaf6	272	if(jet->IsInCone(eta,phi)){ //is cell in cone
5240525e	273	CalculateEEtaPhi(digit,e,eta,phi) ;
	274	jet->AddDigit(e,eta,phi,digit->GetIndexInList()) ;
	275	}
	276	}
	277
	278	//Accept Jet with Et > Et_min and remove all belonging digits
	279	if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
	280	Int_t nIndxs ;
	281	const Int_t * indxs = jet->Indexs(nIndxs) ;
	282	for(Int_t i=0;i<nIndxs;i++){
	283	copyDigits->RemoveAt(indxs[i]) ;
	284	}
	285	jet->CalculateAll() ;
	286	jets->AddAt(jet,fNJets++);
	287	}
	288	else{ //remove jet and do not touch digits
	289	delete jet ;
	290	}
	291	}
	292	else{
	293	if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
	294	break ;
	295	}
	296	}
	297
	298	iIter++ ;
	299	//calculate new energy of backrgound
	300	Double_t oldTotalEnergy = totalEnergy ;
	301	totalEnergy = 0 ;
	302	for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
	303	digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
	304	if(digit)
	305	totalEnergy+=energy[i] ;
	306	}
	307	if(!fMode \|\| (oldTotalEnergy != 0) &&
	308	(TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg))
	309	break ;
	310	}
	311	}
2fa01244	312	delete [] energy;
2fa01244	313	delete [] index;
5240525e	314	copyDigits->Delete() ;
	315
	316	}
	317	//____________________________________________________________________________
	318	Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
f898e0f3	319	//Both simulated and raw digits are already calibrated
5240525e	320
f898e0f3	321	return digit->GetEnergy() ;
5240525e	322
	323	// }
	324	}
	325	//____________________________________________________________________________
	326	void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
0bc3b8ed	327	//Calculate direction of the jet
5240525e	328	e=Calibrate(d) ;
	329	AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
	330	TVector3 pos ;
	331	geom->RelPosInAlice(d->GetId(), pos) ;
	332	eta = pos.Eta() ;
	333	phi = pos.Phi() ;
	334	}
	335	//____________________________________________________________________________
702ab87e	336	void AliPHOSJetFinder::Print(const Option_t *) const {
0bc3b8ed	337	//Print parameters of the found jet
5240525e	338	printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
	339	printf(" Jets found .........%d \n",fNJets) ;
	340	printf(" Seed energy cut ....%f \n",fEtSeed) ;
	341	printf(" Cone radius ........%f \n",fConeRad) ;
	342	printf(" Minimal cone move ..%f \n",fMinConeMove) ;
	343	printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
	344	printf("------------------------------------------------- \n") ;
	345	}