Bug fix
[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 "AliPHOSGetter.h"
47 #include "AliPHOSJetFinder.h"
48 #include "AliPHOSDigitizer.h"
49
50 ClassImp(AliPHOSJetFinder)
51
52
53 //____________________________________________________________________________ 
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)
69 {
70   //Initialize jet parameters
71 }
72
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");
92 }
93
94 //____________________________________________________________________________ 
95   AliPHOSJetFinder::~AliPHOSJetFinder()
96 {
97   //dtor
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 {
111   //Find jets in the case without detector.
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){
197   //Find jets in the case witht detector at the level of digits.
198   if(digits->GetEntries()==0){
199     AliError(Form("No entries in digits list \n")) ;
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 ;
220   Float_t * energy = new Float_t[copyDigits->GetEntries()] ;
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.
230   Int_t * index = new Int_t[copyDigits->GetEntries()] ;
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) ;     
261             if(jet->IsInCone(eta,phi) && //is cell in cone
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));
272           if(jet->IsInCone(eta,phi)){ //is cell in cone
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   }
312   delete [] energy;
313   delete [] index;
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 ){
321 //       AliError(Form("Either Pedestals of Gains not set!")) ;
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
332     AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
333     if(!gime){
334       AliError(Form("Can not read Calibration parameters")) ;
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){
351   //Calculate direction of the jet
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 //____________________________________________________________________________ 
360 void AliPHOSJetFinder::Print(const Option_t *) const {  
361   //Print parameters of the found jet
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 }