1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //////////////////////////////////////////////////////////////////////////////
26 /* History of cvs commits:
29 * Revision 1.8 2005/05/28 14:19:04 schutz
30 * Compilation warnings fixed by T.P.
34 // --- ROOT system ---
35 #include "TClonesArray.h"
37 #include "TParticle.h"
39 // --- Standard library ---
41 // --- AliRoot header files ---
43 #include "AliPHOSJet.h"
44 #include "AliPHOSGeometry.h"
45 #include "AliPHOSDigit.h"
46 #include "AliPHOSJetFinder.h"
47 #include "AliPHOSDigitizer.h"
49 ClassImp(AliPHOSJetFinder)
52 //____________________________________________________________________________
53 AliPHOSJetFinder::AliPHOSJetFinder():
54 TNamed("AliPHOSJetFinder",""),
69 //Initialize jet parameters
72 //____________________________________________________________________________
73 AliPHOSJetFinder::AliPHOSJetFinder(const AliPHOSJetFinder & jet) :
89 // copy ctor: no implementation yet
90 Fatal("cpy ctor", "not implemented");
93 //____________________________________________________________________________
94 AliPHOSJetFinder::~AliPHOSJetFinder()
107 //____________________________________________________________________________
108 void AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist)
110 //Find jets in the case without detector.
113 TIter nextJet(jetslist) ;
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
128 while((jet=static_cast<AliPHOSJet*>(nextJet()))){
130 if(jet->AcceptConeDeviation(p)){
132 //cout << "false" << endl ;
137 //cout << "new " << endl ;
138 jet = new AliPHOSJet() ;
140 // jet = static_cast<AliPHOSJet*>(jetslist->Last()) ;
141 jet->SetConeRadius(fConeRad) ;
142 jet->SetMaxConeMove(fMaxConeMove) ;
143 // jet->SetMinConeMove(fMinConeMove) ;
144 jet->AddParticle(p,-1) ;
148 while((jet=static_cast<AliPHOSJet*>(nextJet()))){
149 if(jet->AcceptConeDeviation(p))
150 jet->AddParticle(p,-1) ; //Just recalculate direction of jet
155 //now calculate directions of jets using collected information
157 while((jet=static_cast<AliPHOSJet*>(nextJet()))){
158 jet->CalculateAll() ;
159 if(jet->Energy() < fEtMin){
160 jetslist->Remove(jet) ;
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
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);
183 (static_cast<AliPHOSJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
187 //Calculate jet parameters
189 while((jet=static_cast<AliPHOSJet*>(nextJet()))){
190 jet->CalculateAll() ;
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")) ;
203 TClonesArray * copyDigits = new TClonesArray(*digits) ;
205 //Remove CPV digits if any
206 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
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) ;
216 copyDigits->Compress() ;
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] ;
228 //Sort digits in decreasing energy.
229 Int_t * index = new Int_t[copyDigits->GetEntries()] ;
230 TMath::Sort(copyDigits->GetEntries(),energy,index) ;
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) ;
245 while(iIter < 10){//less than 10 iterations
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() ;
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
265 }//end of loop over cells
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()) ;
277 //Accept Jet with Et > Et_min and remove all belonging digits
278 if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
280 const Int_t * indxs = jet->Indexs(nIndxs) ;
281 for(Int_t i=0;i<nIndxs;i++){
282 copyDigits->RemoveAt(indxs[i]) ;
284 jet->CalculateAll() ;
285 jets->AddAt(jet,fNJets++);
287 else{ //remove jet and do not touch digits
292 if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
298 //calculate new energy of backrgound
299 Double_t oldTotalEnergy = totalEnergy ;
301 for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
302 digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
304 totalEnergy+=energy[i] ;
306 if(!fMode || (oldTotalEnergy != 0) &&
307 (TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg))
313 copyDigits->Delete() ;
316 //____________________________________________________________________________
317 Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
318 //Both simulated and raw digits are already calibrated
320 return digit->GetEnergy() ;
324 //____________________________________________________________________________
325 void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
326 //Calculate direction of the jet
328 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
330 geom->RelPosInAlice(d->GetId(), pos) ;
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") ;