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