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

//_________________________________________________________________________
// Class to calculate jet chararacteristics
//
//*-- Author :  D.Peressounko
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---
#include "TParticle.h"

// --- Standard library ---

// --- AliRoot header files ---
#include "AliPHOSJet.h"
//  #include "AliPHOSGetter.h"

ClassImp(AliPHOSJet)
  
//____________________________________________________________________________ 
AliPHOSJet::AliPHOSJet():TObject() {
  fNpart = 0 ;
  fList = 0 ;
  //  fMode = 0 ;
  fConeRad = 0 ;
  fMaxConeMove = 0 ;
  fMinConeMove = 0 ;
  
  fSumEnergy = 0 ;
  fSumEta = 0 ;
  fSumPhi = 0 ; 
  fEnergy = 0 ;
  fEta = 0 ;
  fPhi = 0 ; 
  fLEnergy = 0 ;
  fLEta = 0 ;
  fLPhi = 0 ; 
}

//____________________________________________________________________________ 
AliPHOSJet::~AliPHOSJet(){
  if(fList){
    delete fList ;
    fList = 0 ;
  }
  
}
//____________________________________________________________________________ 
void AliPHOSJet::AddParticle(const TParticle * p,const Int_t index){
  //adds particle to jet. Calculates change in jet direction, 
  //due to addition of this particle and if it is smaller, than fMaxDev, 
  //add particle, axcept new direction and return true.
  
  fSumEnergy+=p->Energy() ;
  if(p->Pt()/p->Energy() > 0.001)
    fSumEta+=p->Eta()*p->Energy() ;
  else
    fSumEta+=100.*p->Pz() ; 
  fSumPhi+=p->Phi()*p->Energy() ;    
  
  //check if this a leading particle?
  if(fLEnergy < p->Energy()){
    fLEnergy = p->Energy() ;
    if(p->Pt()/p->Energy() > 0.001)
      fLEta = p->Eta() ;
    else 
      fLEta = 100.*p->Pz()/p->Energy() ;
    fLPhi = p->Phi() ;
  }

  if(index >=0){ //add index to list of indexes
    if(!fList){
      fList = new TArrayI(100) ;
    }
    if(fList->GetSize()<=fNpart+1){
      TArrayI * tmp = new TArrayI(fNpart*2) ;
      tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
      fList=tmp ;
    }
    fList->AddAt(index,fNpart++) ;
  }
}
//____________________________________________________________________________ 
void AliPHOSJet::AddDigit(const Double_t e,const Double_t eta,const Double_t phi, const Int_t index){
  //adds particle to jet. Calculates change in jet direction, 
  //due to addition of this particle and if it is smaller, than fMaxDev, 
  //add particle, axcept new direction and return true.
  
  fSumEnergy+=e ;
  fSumEta+=eta*e ;
  fSumPhi+=phi*e ;    
  
  //check if this a leading particle?
  if(fLEnergy < e){
    fLEnergy = e;
    fLEta = eta ;
    fLPhi = phi ;
  }

  if(index >=0){ //add index to list of indexes
    if(!fList){
      fList = new TArrayI(100) ;
    }
    if(fList->GetSize()<=fNpart+1){
      TArrayI * tmp = new TArrayI(fNpart*2) ;
      tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
      fList=tmp ;
    }
    fList->AddAt(index,fNpart++) ;
  }
}
// //____________________________________________________________________________ 
// Bool_t AliPHOSJet::IsInJet(TParticle * p,Int_t mode)
// {
//   Double_t dEta ;
//   Double_t dPhi ;
//   Double_t energy ;
//   if(!fEnergy){ //Final values not calculated yet, use intermediate
//     if(fSumEnergy==0)
//       return kTRUE ; //First particle
// note p->Eta() causes fpe! 
//     dEta=(p->Eta() - fSumEta/fSumEnergy) ;
//     dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
//     energy = fSumEnergy ;
//   }
//   else{ //calculated with respect to final 
//     dEta=(p->Eta() - fEta) ;
//     dPhi=(p->Phi() - fPhi) ;
//     energy = fEnergy ;    
//   }
  
//   switch (mode){
//   case 0: 
//     return IsInCone(eta,phi) ; //pure geometrical comparison
//   case 1: 
//     return AcceptConeDeviation(dEta,dPhi,p->Energy() );
//   default:
//     Error("IsInCone","Unknown mode of cone calculation %d \n",mode );
//   }
//   return kFALSE ;
//}
//____________________________________________________________________________ 
Bool_t AliPHOSJet::AcceptConeDeviation(const TParticle * p)const
{ //Calculate cone deviation in case of inclusion of the given
  //particle to jet. 

  Double_t tmpEnergy = fSumEnergy + p->Energy() ;  
  Double_t tmpEta = fSumEta ;
  if(p->Pt()/p->Energy() >0.001)
    tmpEta+= p->Eta()*p->Energy() ;
  else
    tmpEta+= 100.*p->Pz() ;
  tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) - 1 ;
  Double_t tmpPhi = fSumPhi + p->Phi()*p->Energy() ;
  tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) - 1 ;
  Double_t dev = TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
  if(dev < fMaxConeMove)
    return kTRUE ;
  else
    return kFALSE ;
}
//____________________________________________________________________________ 
Bool_t AliPHOSJet::AcceptConeDeviation(const Double_t e,const Double_t eta,const Double_t phi)const
{ //Calculate cone deviation in case of inclusion of the given
  //particle to jet. 

  Double_t tmpEnergy = fSumEnergy + e ;
  Double_t tmpEta = fSumEta + eta*e ;
  tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) ;
  Double_t tmpPhi = fSumPhi + phi*e ;
  tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) ;
  Double_t dev = TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
  if(dev<fMaxConeMove && dev > fMinConeMove)
    return kTRUE ;
  else
    return kFALSE ;
}
//____________________________________________________________________________ 
Bool_t AliPHOSJet::IsInCone(const TParticle * p)const
{//
  Double_t dEta ;
  Double_t dPhi ;
  if(!fEnergy){ //Final values not calculated yet, use intermediate
    if(fSumEnergy==0)
      return kTRUE ; //First particle    
    if(p->Pt()/p->Energy() > 0.001)
      dEta=(p->Eta() - fSumEta/fSumEnergy) ;
    else
      dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;
    dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
  }
  else{ //calculated with respect to final 
    if(p->Pt()/p->Energy() > 0.001)
      dEta=(p->Eta() - fEta) ;
    else
      dEta=(100.*p->Pz()/p->Energy() - fEta) ;
    dPhi=(p->Phi() - fPhi) ;
  }
  if(TMath::Sqrt(dEta*dEta + dPhi*dPhi) < fConeRad)
    return kTRUE ;
  else
    return kFALSE ;
}
//____________________________________________________________________________ 
Bool_t AliPHOSJet::IsInCone(const Double_t e,const Double_t eta,const Double_t phi)const
{//
  Double_t dEta ;
  Double_t dPhi ;
  if(!fEnergy){ //Final values not calculated yet, use intermediate
    if(fSumEnergy==0)
      return kTRUE ; //First particle    
    dEta=(eta - fSumEta/fSumEnergy) ;
    dPhi=(phi - fSumPhi/fSumEnergy) ;
  }
  else{ //calculated with respect to final 
    dEta=(eta - fEta) ;
    dPhi=(phi - fPhi) ;
  }
  if(TMath::Sqrt(dEta*dEta + dPhi*dPhi) < fConeRad)
    return kTRUE ;
  else
    return kFALSE ;
}
//____________________________________________________________________________ 
Double_t AliPHOSJet::DistanceToJet(const TParticle *p)const{
  //Calculate radius 
  Double_t dEta ;
  Double_t dPhi ;
  if(!fEnergy){ //Final values not calculated yet, use intermediate
    if(fSumEnergy==0)
      return kTRUE ; //First particle    
    if(p->Pt()/p->Energy() > 0.001)
      dEta=(p->Eta() - fSumEta/fSumEnergy) ;
    else
      dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;      
    dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
  }
  else{ //calculated with respect to final 
    if(p->Pt()/p->Energy() > 0.001)
      dEta=(p->Eta() - fEta) ;
    else
      dEta=(100*p->Pz()/p->Energy() - fEta) ;    
    dPhi=(p->Phi() - fPhi) ;
  }
  return TMath::Sqrt(dEta*dEta + dPhi*dPhi) ;

}
//____________________________________________________________________________ 
void AliPHOSJet::CalculateAll(void){

  if(fSumEnergy==0)
    return  ; //Nothing to calculate    
  
  fEta = fSumEta/fSumEnergy ;
  fPhi = fSumPhi/fSumEnergy ;
  fEnergy = fSumEnergy ;
  
  fSumEnergy = 0. ;
  fSumEta = 0. ;
  fSumPhi = 0. ;
}
//____________________________________________________________________________ 
void AliPHOSJet::Print(Option_t * option){
  printf("-------------- AliPHOSJet ------------\n") ;
  printf(" Energy............. %f \n",fEnergy) ;
  printf(" Eta................ %f \n",fEta ) ;
  printf(" Phi................ %f \n",fPhi ) ;
  printf(" Leading Energy..... %f \n",fLEnergy) ;
  printf(" Leading Eta........ %f \n",fLEta) ;
  printf(" Leading Phi........ %f \n",fLPhi) ;
  printf(" N particles in jet  %d \n",fNpart) ;
  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	// Class to calculate jet chararacteristics
	18	//
	19	//*-- Author : D.Peressounko
	20	//////////////////////////////////////////////////////////////////////////////
	21
	22	// --- ROOT system ---
	23	#include "TParticle.h"
	24
	25	// --- Standard library ---
	26
	27	// --- AliRoot header files ---
	28	#include "AliPHOSJet.h"
	29	// #include "AliPHOSGetter.h"
	30
	31	ClassImp(AliPHOSJet)
	32
	33	//____________________________________________________________________________
	34	AliPHOSJet::AliPHOSJet():TObject() {
	35	fNpart = 0 ;
	36	fList = 0 ;
	37	// fMode = 0 ;
	38	fConeRad = 0 ;
	39	fMaxConeMove = 0 ;
	40	fMinConeMove = 0 ;
	41
	42	fSumEnergy = 0 ;
	43	fSumEta = 0 ;
	44	fSumPhi = 0 ;
	45	fEnergy = 0 ;
	46	fEta = 0 ;
	47	fPhi = 0 ;
	48	fLEnergy = 0 ;
	49	fLEta = 0 ;
	50	fLPhi = 0 ;
	51	}
	52
	53	//____________________________________________________________________________
	54	AliPHOSJet::~AliPHOSJet(){
	55	if(fList){
	56	delete fList ;
	57	fList = 0 ;
	58	}
	59
	60	}
	61	//____________________________________________________________________________
	62	void AliPHOSJet::AddParticle(const TParticle * p,const Int_t index){
	63	//adds particle to jet. Calculates change in jet direction,
	64	//due to addition of this particle and if it is smaller, than fMaxDev,
65	//add particle, axcept new direction and return true.
66
67	fSumEnergy+=p->Energy() ;
68	if(p->Pt()/p->Energy() > 0.001)
69	fSumEta+=p->Eta()*p->Energy() ;
70	else
71	fSumEta+=100.*p->Pz() ;
72	fSumPhi+=p->Phi()*p->Energy() ;
73
74	//check if this a leading particle?
75	if(fLEnergy < p->Energy()){
76	fLEnergy = p->Energy() ;
77	if(p->Pt()/p->Energy() > 0.001)
78	fLEta = p->Eta() ;
79	else
80	fLEta = 100.*p->Pz()/p->Energy() ;
81	fLPhi = p->Phi() ;
82	}
83
84	if(index >=0){ //add index to list of indexes
85	if(!fList){
86	fList = new TArrayI(100) ;
87	}
88	if(fList->GetSize()<=fNpart+1){
89	TArrayI * tmp = new TArrayI(fNpart*2) ;
90	tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
91	fList=tmp ;
92	}
93	fList->AddAt(index,fNpart++) ;
94	}
95	}
96	//____________________________________________________________________________
97	void AliPHOSJet::AddDigit(const Double_t e,const Double_t eta,const Double_t phi, const Int_t index){
98	//adds particle to jet. Calculates change in jet direction,
99	//due to addition of this particle and if it is smaller, than fMaxDev,
100	//add particle, axcept new direction and return true.
101
102	fSumEnergy+=e ;
103	fSumEta+=eta*e ;
104	fSumPhi+=phi*e ;
105
106	//check if this a leading particle?
107	if(fLEnergy < e){
108	fLEnergy = e;
109	fLEta = eta ;
110	fLPhi = phi ;
111	}
112
113	if(index >=0){ //add index to list of indexes
114	if(!fList){
115	fList = new TArrayI(100) ;
116	}
117	if(fList->GetSize()<=fNpart+1){
118	TArrayI * tmp = new TArrayI(fNpart*2) ;
119	tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
120	fList=tmp ;
121	}
122	fList->AddAt(index,fNpart++) ;
123	}
124	}
125	// //____________________________________________________________________________
126	// Bool_t AliPHOSJet::IsInJet(TParticle * p,Int_t mode)
127	// {
128	// Double_t dEta ;
129	// Double_t dPhi ;
130	// Double_t energy ;
131	// if(!fEnergy){ //Final values not calculated yet, use intermediate
132	// if(fSumEnergy==0)
133	// return kTRUE ; //First particle
134	// note p->Eta() causes fpe!
135	// dEta=(p->Eta() - fSumEta/fSumEnergy) ;
136	// dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
137	// energy = fSumEnergy ;
138	// }
139	// else{ //calculated with respect to final
140	// dEta=(p->Eta() - fEta) ;
141	// dPhi=(p->Phi() - fPhi) ;
142	// energy = fEnergy ;
143	// }
144
145	// switch (mode){
146	// case 0:
147	// return IsInCone(eta,phi) ; //pure geometrical comparison
148	// case 1:
149	// return AcceptConeDeviation(dEta,dPhi,p->Energy() );
150	// default:
151	// Error("IsInCone","Unknown mode of cone calculation %d \n",mode );
152	// }
153	// return kFALSE ;
154	//}
155	//____________________________________________________________________________
156	Bool_t AliPHOSJet::AcceptConeDeviation(const TParticle * p)const
157	{ //Calculate cone deviation in case of inclusion of the given
158	//particle to jet.
159
160	Double_t tmpEnergy = fSumEnergy + p->Energy() ;
161	Double_t tmpEta = fSumEta ;
162	if(p->Pt()/p->Energy() >0.001)
163	tmpEta+= p->Eta()*p->Energy() ;
164	else
165	tmpEta+= 100.*p->Pz() ;
166	tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) - 1 ;
167	Double_t tmpPhi = fSumPhi + p->Phi()*p->Energy() ;
168	tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) - 1 ;
169	Double_t dev = TMath::Sqrt(tmpEtatmpEta + tmpPhitmpPhi) ;
170	if(dev < fMaxConeMove)
171	return kTRUE ;
172	else
173	return kFALSE ;
174	}
175	//____________________________________________________________________________
176	Bool_t AliPHOSJet::AcceptConeDeviation(const Double_t e,const Double_t eta,const Double_t phi)const
177	{ //Calculate cone deviation in case of inclusion of the given
178	//particle to jet.
179
180	Double_t tmpEnergy = fSumEnergy + e ;
181	Double_t tmpEta = fSumEta + eta*e ;
182	tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) ;
183	Double_t tmpPhi = fSumPhi + phi*e ;
184	tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) ;
185	Double_t dev = TMath::Sqrt(tmpEtatmpEta + tmpPhitmpPhi) ;
186	if(dev<fMaxConeMove && dev > fMinConeMove)
187	return kTRUE ;
188	else
189	return kFALSE ;
190	}
191	//____________________________________________________________________________
192	Bool_t AliPHOSJet::IsInCone(const TParticle * p)const
193	{//
194	Double_t dEta ;
195	Double_t dPhi ;
196	if(!fEnergy){ //Final values not calculated yet, use intermediate
197	if(fSumEnergy==0)
198	return kTRUE ; //First particle
199	if(p->Pt()/p->Energy() > 0.001)
200	dEta=(p->Eta() - fSumEta/fSumEnergy) ;
201	else
202	dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;
203	dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
204	}
205	else{ //calculated with respect to final
206	if(p->Pt()/p->Energy() > 0.001)
207	dEta=(p->Eta() - fEta) ;
208	else
209	dEta=(100.*p->Pz()/p->Energy() - fEta) ;
210	dPhi=(p->Phi() - fPhi) ;
211	}
212	if(TMath::Sqrt(dEtadEta + dPhidPhi) < fConeRad)
213	return kTRUE ;
214	else
215	return kFALSE ;
216	}
217	//____________________________________________________________________________
218	Bool_t AliPHOSJet::IsInCone(const Double_t e,const Double_t eta,const Double_t phi)const
219	{//
220	Double_t dEta ;
221	Double_t dPhi ;
222	if(!fEnergy){ //Final values not calculated yet, use intermediate
223	if(fSumEnergy==0)
224	return kTRUE ; //First particle
225	dEta=(eta - fSumEta/fSumEnergy) ;
226	dPhi=(phi - fSumPhi/fSumEnergy) ;
227	}
228	else{ //calculated with respect to final
229	dEta=(eta - fEta) ;
230	dPhi=(phi - fPhi) ;
231	}
232	if(TMath::Sqrt(dEtadEta + dPhidPhi) < fConeRad)
233	return kTRUE ;
234	else
235	return kFALSE ;
236	}
237	//____________________________________________________________________________
238	Double_t AliPHOSJet::DistanceToJet(const TParticle *p)const{
239	//Calculate radius
240	Double_t dEta ;
241	Double_t dPhi ;
242	if(!fEnergy){ //Final values not calculated yet, use intermediate
243	if(fSumEnergy==0)
244	return kTRUE ; //First particle
245	if(p->Pt()/p->Energy() > 0.001)
246	dEta=(p->Eta() - fSumEta/fSumEnergy) ;
247	else
248	dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;
249	dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
250	}
251	else{ //calculated with respect to final
252	if(p->Pt()/p->Energy() > 0.001)
253	dEta=(p->Eta() - fEta) ;
254	else
255	dEta=(100*p->Pz()/p->Energy() - fEta) ;
256	dPhi=(p->Phi() - fPhi) ;
257	}
258	return TMath::Sqrt(dEtadEta + dPhidPhi) ;
259
260	}
261	//____________________________________________________________________________
262	void AliPHOSJet::CalculateAll(void){
263
264	if(fSumEnergy==0)
265	return ; //Nothing to calculate
266
267	fEta = fSumEta/fSumEnergy ;
268	fPhi = fSumPhi/fSumEnergy ;
269	fEnergy = fSumEnergy ;
270
271	fSumEnergy = 0. ;
272	fSumEta = 0. ;
273	fSumPhi = 0. ;
274	}
275	//____________________________________________________________________________
276	void AliPHOSJet::Print(Option_t * option){
277	printf("-------------- AliPHOSJet ------------\n") ;
278	printf(" Energy............. %f \n",fEnergy) ;
279	printf(" Eta................ %f \n",fEta ) ;
280	printf(" Phi................ %f \n",fPhi ) ;
281	printf(" Leading Energy..... %f \n",fLEnergy) ;
282	printf(" Leading Eta........ %f \n",fLEta) ;
283	printf(" Leading Phi........ %f \n",fLPhi) ;
284	printf(" N particles in jet %d \n",fNpart) ;
285	printf("----------------------------------\n") ;
286	}