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