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