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