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