]>
Commit | Line | Data |
---|---|---|
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 | ||
50 | ClassImp(AliPHOSJetFinder) | |
51 | ||
52 | ||
53 | //____________________________________________________________________________ | |
3f7dbdb7 | 54 | AliPHOSJetFinder::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 | //____________________________________________________________________________ |
74 | AliPHOSJetFinder::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 | //____________________________________________________________________________ | |
109 | void 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 | //____________________________________________________________________________ | |
196 | void 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 | //____________________________________________________________________________ | |
318 | Double_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 | //____________________________________________________________________________ | |
350 | void 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 | 360 | void 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 | } |