PHOS/AliPHOSJetFinder.cxx.old

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