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