1258825da39bcd673a14e8a58524227e020695ca
[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 //
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
37 ClassImp(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 //____________________________________________________________________________ 
75 void  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 //____________________________________________________________________________ 
162 void 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 = new Float_t[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 = new Int_t[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   delete [] energy;
279   delete [] index;
280   copyDigits->Delete() ;
281   
282 }
283 //____________________________________________________________________________ 
284 Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){ 
285 //   if(fPedestals || fGains ){  //use calibration data
286 //     if(!fPedestals || !fGains ){
287 //       Error("Calibrate","Either Pedestals of Gains not set!") ;
288 //       return 0 ;
289 //     }
290 //     Float_t en=(digit->GetAmp() - fPedestals->Data(digit->GetId)()))*fGains->Data(digit->GetId()) ;
291 //     if(en>0)
292 //       return en ;
293 //     else
294 //       return 0. ;
295 //   }
296 //   else{ //simulation
297   if(fSimGain==0){ //read simulation parameters
298     AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
299     if(!gime){
300       Error("Calibrate","Can not read Calibration parameters") ;
301       return 0 ;
302     }
303     const TTask * task = gime->Digitizer() ;
304     if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
305       const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
306       fSimGain     = dig->GetEMCchannel() ;
307       fSimPedestal = dig->GetEMCpedestal();
308     }
309   }
310   
311   return fSimPedestal + digit->GetAmp()*fSimGain ;        
312   
313   //  }  
314 }
315 //____________________________________________________________________________ 
316 void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
317   e=Calibrate(d) ;
318   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
319   TVector3 pos ;
320   geom->RelPosInAlice(d->GetId(), pos) ;
321   eta = pos.Eta() ;
322   phi = pos.Phi() ;
323 }
324 //____________________________________________________________________________ 
325 void AliPHOSJetFinder::Print(Option_t * option){        
326   printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
327   printf(" Jets found .........%d \n",fNJets) ;
328   printf(" Seed energy cut ....%f \n",fEtSeed) ;
329   printf(" Cone radius ........%f \n",fConeRad) ;
330   printf(" Minimal cone move ..%f \n",fMinConeMove) ;
331   printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
332   printf("------------------------------------------------- \n") ;
333 }