Transition to NewIO
[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 "AliPHOSDigit.h"
34 #include "AliPHOSGetter.h"
35 #include "AliPHOSJetFinder.h"
36 #include "AliPHOSDigitizer.h"
37
38 ClassImp(AliPHOSJetFinder)
39
40
41 //____________________________________________________________________________ 
42   AliPHOSJetFinder::AliPHOSJetFinder():TNamed("AliPHOSJetFinder","") 
43 {
44   //Initialize jet parameters
45   fNJets = 0 ; 
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 {
66   //dtor
67   if(fParticles){
68     delete fParticles ;
69     fParticles = 0 ;
70   }
71   if(fJets){
72     delete fJets ;
73     fJets = 0 ;
74   } 
75 }
76
77 //____________________________________________________________________________ 
78 void  AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist) 
79 {
80   //Find jets in the case without detector.
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 //____________________________________________________________________________ 
165 void AliPHOSJetFinder::FindJetsFromDigits(const TClonesArray * digits, TObjArray * jets){
166   //Find jets in the case witht detector at the level of digits.
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 ;
189   Float_t * energy = new Float_t[copyDigits->GetEntries()] ;
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.
199   Int_t * index = new Int_t[copyDigits->GetEntries()] ;
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   }
281   delete [] energy;
282   delete [] index;
283   copyDigits->Delete() ;
284   
285 }
286 //____________________________________________________________________________ 
287 Double_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
301     AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
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 //____________________________________________________________________________ 
319 void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
320   //Calculate direction of the jet
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 //____________________________________________________________________________ 
329 void AliPHOSJetFinder::Print(Option_t * option){        
330   //Print parameters of the found jet
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 }