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