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