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