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