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