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