a3ef8e47efe946e612b2f10f488c0da0292c5ad8
[u/mrichter/AliRoot.git] / PHOS / AliPHOSJet.cxx
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(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
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(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
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(dEta*dEta + dPhi*dPhi) < 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(dEta*dEta + dPhi*dPhi) < 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(dEta*dEta + dPhi*dPhi) ;
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 }