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 "AliPHOSJetFinder.h"
  47 #include "AliPHOSDigitizer.h"
  48
  49 ClassImp(AliPHOSJetFinder)
  50
  51
  52 //____________________________________________________________________________
  53 AliPHOSJetFinder::AliPHOSJetFinder():
  54   TNamed("AliPHOSJetFinder",""),
  55   fNJets(0),
  56   fStatusCode(-999),
  57   fMode(0),
  58   fConeRad(1.),
  59   fMaxConeMove(0.15),
  60   fMinConeMove(0.05),
  61   fEtSeed(4.),
  62   fEtMin(5.),
  63   fPrecBg(0.00035),
  64   fSimGain(0.),
  65   fSimPedestal(0.),
  66   fParticles(0),
  67   fJets(0)
  68 {
  69   //Initialize jet parameters
  70 }
  71
  72 //____________________________________________________________________________
  73 AliPHOSJetFinder::AliPHOSJetFinder(const AliPHOSJetFinder & jet) :
  74   TNamed(jet),
  75   fNJets(0),
  76   fStatusCode(-999),
  77   fMode(0),
  78   fConeRad(1.),
  79   fMaxConeMove(0.15),
  80   fMinConeMove(0.05),
  81   fEtSeed(4.),
  82   fEtMin(5.),
  83   fPrecBg(0.00035),
  84   fSimGain(0.),
  85   fSimPedestal(0.),
  86   fParticles(0),
  87   fJets(0)
  88 {
  89   // copy ctor: no implementation yet
  90   Fatal("cpy ctor", "not implemented");
  91 }
  92
  93 //____________________________________________________________________________
  94   AliPHOSJetFinder::~AliPHOSJetFinder()
  95 {
  96   //dtor
  97   if(fParticles){
  98     delete fParticles ;
  99     fParticles = 0 ;
 100   }
 101   if(fJets){
 102     delete fJets ;
 103     fJets = 0 ;
 104   }
 105 }
 106
 107 //____________________________________________________________________________
 108 void  AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist)
 109 {
 110   //Find jets in the case without detector.
 111   TIter next(plist) ;
 112
 113   TIter nextJet(jetslist) ;
 114
 115   fNJets = 0 ;
 116   TParticle * p ;
 117   AliPHOSJet * jet ;
 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
 127       nextJet.Reset() ;
 128       while((jet=static_cast<AliPHOSJet*>(nextJet()))){
 129         //jet->Print() ;
 130         if(jet->AcceptConeDeviation(p)){
 131           startnew = kFALSE ;
 132           //cout << "false" << endl ;
 133           break ;
 134         }
 135       }
 136       if(startnew){
 137         //cout << "new " << endl ;
 138         jet = new AliPHOSJet() ;
 139         jetslist->Add(jet) ;
 140         //      jet = static_cast<AliPHOSJet*>(jetslist->Last()) ;
 141         jet->SetConeRadius(fConeRad) ;
 142         jet->SetMaxConeMove(fMaxConeMove) ;
 143         //      jet->SetMinConeMove(fMinConeMove) ;
 144         jet->AddParticle(p,-1) ;
 145         fNJets++;
 146       }
 147     }
 148     while((jet=static_cast<AliPHOSJet*>(nextJet()))){
 149       if(jet->AcceptConeDeviation(p))
 150         jet->AddParticle(p,-1) ; //Just recalculate direction of jet
 151     }
 152     }
 153   }
 154
 155   //now calculate directions of jets using collected information
 156   nextJet.Reset() ;
 157   while((jet=static_cast<AliPHOSJet*>(nextJet()))){
 158     jet->CalculateAll() ;
 159     if(jet->Energy() < fEtMin){
 160       jetslist->Remove(jet) ;
 161       delete jet ;
 162     }
 163   }
 164
 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
 171     Int_t iJet = -1 ;
 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);
 176         if(cdist < dist){
 177           dist = cdist ;
 178           iJet = i ;
 179         }
 180       }
 181     }
 182     if(iJet>-1)
 183       (static_cast<AliPHOSJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
 184   }
 185   }
 186
 187   //Calculate jet parameters
 188   nextJet.Reset() ;
 189   while((jet=static_cast<AliPHOSJet*>(nextJet()))){
 190     jet->CalculateAll() ;
 191   }
 192
 193 }
 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")) ;
 199     return ;
 200   }
 201
 202
 203   TClonesArray * copyDigits = new TClonesArray(*digits) ;
 204
 205   //Remove CPV digits if any
 206   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
 207   Int_t iDigit ;
 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) ;
 213     else
 214       break ;
 215   }
 216   copyDigits->Compress() ;
 217
 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] ;
 226   }
 227
 228   //Sort digits in decreasing energy.
 229   Int_t * index = new Int_t[copyDigits->GetEntries()] ;
 230   TMath::Sort(copyDigits->GetEntries(),energy,index) ;
 231
 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) ;
 238     else
 239       break ;
 240   }
 241
 242
 243   AliPHOSJet * jet ;
 244   Int_t iIter = 0 ;
 245   while(iIter < 10){//less than 10 iterations
 246
 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() ;
 252         Double_t e,eta,phi ;
 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
 263             }
 264           }
 265         }//end of loop over cells
 266
 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()) ;
 274           }
 275         }
 276
 277         //Accept Jet with Et > Et_min and remove all belonging digits
 278         if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
 279           Int_t nIndxs ;
 280           const Int_t * indxs = jet->Indexs(nIndxs) ;
 281           for(Int_t i=0;i<nIndxs;i++){
 282             copyDigits->RemoveAt(indxs[i]) ;
 283           }
 284           jet->CalculateAll() ;
 285           jets->AddAt(jet,fNJets++);
 286         }
 287         else{ //remove jet and do not touch digits
 288           delete jet ;
 289         }
 290       }
 291       else{
 292         if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
 293           break ;
 294         }
 295       }
 296
 297       iIter++ ;
 298       //calculate new energy of backrgound
 299       Double_t oldTotalEnergy = totalEnergy ;
 300       totalEnergy = 0 ;
 301       for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
 302         digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
 303         if(digit)
 304           totalEnergy+=energy[i] ;
 305       }
 306       if(!fMode || ((oldTotalEnergy != 0) &&
 307                     (TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg)))
 308         break ;
 309     }
 310   }
 311   delete [] energy;
 312   delete [] index;
 313   copyDigits->Delete() ;
 314
 315 }
 316 //____________________________________________________________________________
 317 Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
 318   //Both simulated and raw digits are already calibrated
 319
 320   return digit->GetEnergy() ;
 321
 322   //  }
 323 }
 324 //____________________________________________________________________________
 325 void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
 326   //Calculate direction of the jet
 327   e=Calibrate(d) ;
 328   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
 329   TVector3 pos ;
 330   geom->RelPosInAlice(d->GetId(), pos) ;
 331   eta = pos.Eta() ;
 332   phi = pos.Phi() ;
 333 }
 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") ;
 344 }