Bug fix
[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
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 */
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
50ClassImp(AliPHOSJetFinder)
51
52
53//____________________________________________________________________________
3f7dbdb7 54AliPHOSJetFinder::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//____________________________________________________________________________
74AliPHOSJetFinder::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//____________________________________________________________________________
109void 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//____________________________________________________________________________
196void 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;
313 delete [] index;
5240525e 314 copyDigits->Delete() ;
315
316}
317//____________________________________________________________________________
318Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
319// if(fPedestals || fGains ){ //use calibration data
320// if(!fPedestals || !fGains ){
351dd634 321// AliError(Form("Either Pedestals of Gains not set!")) ;
5240525e 322// return 0 ;
323// }
324// Float_t en=(digit->GetAmp() - fPedestals->Data(digit->GetId)()))*fGains->Data(digit->GetId()) ;
325// if(en>0)
326// return en ;
327// else
328// return 0. ;
329// }
330// else{ //simulation
331 if(fSimGain==0){ //read simulation parameters
88cb7938 332 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
5240525e 333 if(!gime){
351dd634 334 AliError(Form("Can not read Calibration parameters")) ;
5240525e 335 return 0 ;
336 }
337 const TTask * task = gime->Digitizer() ;
338 if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
339 const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
340 fSimGain = dig->GetEMCchannel() ;
341 fSimPedestal = dig->GetEMCpedestal();
342 }
343 }
344
345 return fSimPedestal + digit->GetAmp()*fSimGain ;
346
347 // }
348}
349//____________________________________________________________________________
350void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
0bc3b8ed 351 //Calculate direction of the jet
5240525e 352 e=Calibrate(d) ;
353 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
354 TVector3 pos ;
355 geom->RelPosInAlice(d->GetId(), pos) ;
356 eta = pos.Eta() ;
357 phi = pos.Phi() ;
358}
359//____________________________________________________________________________
702ab87e 360void AliPHOSJetFinder::Print(const Option_t *) const {
0bc3b8ed 361 //Print parameters of the found jet
5240525e 362 printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
363 printf(" Jets found .........%d \n",fNJets) ;
364 printf(" Seed energy cut ....%f \n",fEtSeed) ;
365 printf(" Cone radius ........%f \n",fConeRad) ;
366 printf(" Minimal cone move ..%f \n",fMinConeMove) ;
367 printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
368 printf("------------------------------------------------- \n") ;
369}