--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// Class to calculate jet chararacteristics
+//
+//*-- Author : D.Peressounko
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+#include "TParticle.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+#include "AliPHOSJet.h"
+// #include "AliPHOSGetter.h"
+
+ClassImp(AliPHOSJet)
+
+//____________________________________________________________________________
+AliPHOSJet::AliPHOSJet():TObject() {
+ fNpart = 0 ;
+ fList = 0 ;
+ // fMode = 0 ;
+ fConeRad = 0 ;
+ fMaxConeMove = 0 ;
+ fMinConeMove = 0 ;
+
+ fSumEnergy = 0 ;
+ fSumEta = 0 ;
+ fSumPhi = 0 ;
+ fEnergy = 0 ;
+ fEta = 0 ;
+ fPhi = 0 ;
+ fLEnergy = 0 ;
+ fLEta = 0 ;
+ fLPhi = 0 ;
+}
+
+//____________________________________________________________________________
+AliPHOSJet::~AliPHOSJet(){
+ if(fList){
+ delete fList ;
+ fList = 0 ;
+ }
+
+}
+//____________________________________________________________________________
+void AliPHOSJet::AddParticle(const TParticle * p,const Int_t index){
+ //adds particle to jet. Calculates change in jet direction,
+ //due to addition of this particle and if it is smaller, than fMaxDev,
+ //add particle, axcept new direction and return true.
+
+ fSumEnergy+=p->Energy() ;
+ if(p->Pt()/p->Energy() > 0.001)
+ fSumEta+=p->Eta()*p->Energy() ;
+ else
+ fSumEta+=100.*p->Pz() ;
+ fSumPhi+=p->Phi()*p->Energy() ;
+
+ //check if this a leading particle?
+ if(fLEnergy < p->Energy()){
+ fLEnergy = p->Energy() ;
+ if(p->Pt()/p->Energy() > 0.001)
+ fLEta = p->Eta() ;
+ else
+ fLEta = 100.*p->Pz()/p->Energy() ;
+ fLPhi = p->Phi() ;
+ }
+
+ if(index >=0){ //add index to list of indexes
+ if(!fList){
+ fList = new TArrayI(100) ;
+ }
+ if(fList->GetSize()<=fNpart+1){
+ TArrayI * tmp = new TArrayI(fNpart*2) ;
+ tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
+ fList=tmp ;
+ }
+ fList->AddAt(index,fNpart++) ;
+ }
+}
+//____________________________________________________________________________
+void AliPHOSJet::AddDigit(const Double_t e,const Double_t eta,const Double_t phi, const Int_t index){
+ //adds particle to jet. Calculates change in jet direction,
+ //due to addition of this particle and if it is smaller, than fMaxDev,
+ //add particle, axcept new direction and return true.
+
+ fSumEnergy+=e ;
+ fSumEta+=eta*e ;
+ fSumPhi+=phi*e ;
+
+ //check if this a leading particle?
+ if(fLEnergy < e){
+ fLEnergy = e;
+ fLEta = eta ;
+ fLPhi = phi ;
+ }
+
+ if(index >=0){ //add index to list of indexes
+ if(!fList){
+ fList = new TArrayI(100) ;
+ }
+ if(fList->GetSize()<=fNpart+1){
+ TArrayI * tmp = new TArrayI(fNpart*2) ;
+ tmp->Adopt(fList->GetSize(),fList->GetArray()) ; //note, we should not delete old array!
+ fList=tmp ;
+ }
+ fList->AddAt(index,fNpart++) ;
+ }
+}
+// //____________________________________________________________________________
+// Bool_t AliPHOSJet::IsInJet(TParticle * p,Int_t mode)
+// {
+// Double_t dEta ;
+// Double_t dPhi ;
+// Double_t energy ;
+// if(!fEnergy){ //Final values not calculated yet, use intermediate
+// if(fSumEnergy==0)
+// return kTRUE ; //First particle
+// note p->Eta() causes fpe!
+// dEta=(p->Eta() - fSumEta/fSumEnergy) ;
+// dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
+// energy = fSumEnergy ;
+// }
+// else{ //calculated with respect to final
+// dEta=(p->Eta() - fEta) ;
+// dPhi=(p->Phi() - fPhi) ;
+// energy = fEnergy ;
+// }
+
+// switch (mode){
+// case 0:
+// return IsInCone(eta,phi) ; //pure geometrical comparison
+// case 1:
+// return AcceptConeDeviation(dEta,dPhi,p->Energy() );
+// default:
+// Error("IsInCone","Unknown mode of cone calculation %d \n",mode );
+// }
+// return kFALSE ;
+//}
+//____________________________________________________________________________
+Bool_t AliPHOSJet::AcceptConeDeviation(const TParticle * p)const
+{ //Calculate cone deviation in case of inclusion of the given
+ //particle to jet.
+
+ Double_t tmpEnergy = fSumEnergy + p->Energy() ;
+ Double_t tmpEta = fSumEta ;
+ if(p->Pt()/p->Energy() >0.001)
+ tmpEta+= p->Eta()*p->Energy() ;
+ else
+ tmpEta+= 100.*p->Pz() ;
+ tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) - 1 ;
+ Double_t tmpPhi = fSumPhi + p->Phi()*p->Energy() ;
+ tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) - 1 ;
+ Double_t dev = TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
+ if(dev < fMaxConeMove)
+ return kTRUE ;
+ else
+ return kFALSE ;
+}
+//____________________________________________________________________________
+Bool_t AliPHOSJet::AcceptConeDeviation(const Double_t e,const Double_t eta,const Double_t phi)const
+{ //Calculate cone deviation in case of inclusion of the given
+ //particle to jet.
+
+ Double_t tmpEnergy = fSumEnergy + e ;
+ Double_t tmpEta = fSumEta + eta*e ;
+ tmpEta = tmpEta/tmpEnergy / (fSumEta/fSumEnergy) ;
+ Double_t tmpPhi = fSumPhi + phi*e ;
+ tmpPhi = tmpPhi/tmpEnergy /(fSumPhi/fSumEnergy) ;
+ Double_t dev = TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi) ;
+ if(dev<fMaxConeMove && dev > fMinConeMove)
+ return kTRUE ;
+ else
+ return kFALSE ;
+}
+//____________________________________________________________________________
+Bool_t AliPHOSJet::IsInCone(const TParticle * p)const
+{//
+ Double_t dEta ;
+ Double_t dPhi ;
+ if(!fEnergy){ //Final values not calculated yet, use intermediate
+ if(fSumEnergy==0)
+ return kTRUE ; //First particle
+ if(p->Pt()/p->Energy() > 0.001)
+ dEta=(p->Eta() - fSumEta/fSumEnergy) ;
+ else
+ dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;
+ dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
+ }
+ else{ //calculated with respect to final
+ if(p->Pt()/p->Energy() > 0.001)
+ dEta=(p->Eta() - fEta) ;
+ else
+ dEta=(100.*p->Pz()/p->Energy() - fEta) ;
+ dPhi=(p->Phi() - fPhi) ;
+ }
+ if(TMath::Sqrt(dEta*dEta + dPhi*dPhi) < fConeRad)
+ return kTRUE ;
+ else
+ return kFALSE ;
+}
+//____________________________________________________________________________
+Bool_t AliPHOSJet::IsInCone(const Double_t e,const Double_t eta,const Double_t phi)const
+{//
+ Double_t dEta ;
+ Double_t dPhi ;
+ if(!fEnergy){ //Final values not calculated yet, use intermediate
+ if(fSumEnergy==0)
+ return kTRUE ; //First particle
+ dEta=(eta - fSumEta/fSumEnergy) ;
+ dPhi=(phi - fSumPhi/fSumEnergy) ;
+ }
+ else{ //calculated with respect to final
+ dEta=(eta - fEta) ;
+ dPhi=(phi - fPhi) ;
+ }
+ if(TMath::Sqrt(dEta*dEta + dPhi*dPhi) < fConeRad)
+ return kTRUE ;
+ else
+ return kFALSE ;
+}
+//____________________________________________________________________________
+Double_t AliPHOSJet::DistanceToJet(const TParticle *p)const{
+ //Calculate radius
+ Double_t dEta ;
+ Double_t dPhi ;
+ if(!fEnergy){ //Final values not calculated yet, use intermediate
+ if(fSumEnergy==0)
+ return kTRUE ; //First particle
+ if(p->Pt()/p->Energy() > 0.001)
+ dEta=(p->Eta() - fSumEta/fSumEnergy) ;
+ else
+ dEta=(100.*p->Pz()/p->Energy() - fSumEta/fSumEnergy) ;
+ dPhi=(p->Phi() - fSumPhi/fSumEnergy) ;
+ }
+ else{ //calculated with respect to final
+ if(p->Pt()/p->Energy() > 0.001)
+ dEta=(p->Eta() - fEta) ;
+ else
+ dEta=(100*p->Pz()/p->Energy() - fEta) ;
+ dPhi=(p->Phi() - fPhi) ;
+ }
+ return TMath::Sqrt(dEta*dEta + dPhi*dPhi) ;
+
+}
+//____________________________________________________________________________
+void AliPHOSJet::CalculateAll(void){
+
+ if(fSumEnergy==0)
+ return ; //Nothing to calculate
+
+ fEta = fSumEta/fSumEnergy ;
+ fPhi = fSumPhi/fSumEnergy ;
+ fEnergy = fSumEnergy ;
+
+ fSumEnergy = 0. ;
+ fSumEta = 0. ;
+ fSumPhi = 0. ;
+}
+//____________________________________________________________________________
+void AliPHOSJet::Print(Option_t * option){
+ printf("-------------- AliPHOSJet ------------\n") ;
+ printf(" Energy............. %f \n",fEnergy) ;
+ printf(" Eta................ %f \n",fEta ) ;
+ printf(" Phi................ %f \n",fPhi ) ;
+ printf(" Leading Energy..... %f \n",fLEnergy) ;
+ printf(" Leading Eta........ %f \n",fLEta) ;
+ printf(" Leading Phi........ %f \n",fLPhi) ;
+ printf(" N particles in jet %d \n",fNpart) ;
+ printf("----------------------------------\n") ;
+}
--- /dev/null
+#ifndef ALIJET_H
+#define ALIJET_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+// Base Class for Jets in ALICE
+//
+//*-- Author: D.Peressounko
+
+
+// --- ROOT system ---
+#include "TObject.h"
+class TParticle ;
+#include "TArrayI.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+class AliPHOSJet : public TObject {
+
+public:
+ AliPHOSJet() ; // ctor
+ AliPHOSJet(const AliPHOSJet & jet){} ; // ctor
+ virtual ~AliPHOSJet() ;
+
+ void AddDigit(const Double_t e,const Double_t eta,const Double_t phi,const Int_t index) ;
+ void AddParticle(const TParticle * p,const Int_t index) ;
+ //adds particle p to jet. index: index of p in list of all particles in event
+
+ Double_t DistanceToJet(const TParticle *p)const ;
+ //calculates distance to Jet in accordance with some scheam:
+ //geometrical, inv mass, etc
+
+ void CalculateAll(void) ;
+ //calculate final Energy, Eta & phi from intermediate ones.
+
+ const Int_t * Indexs(Int_t & nIndexs)const{nIndexs = fNpart; return fList->GetArray() ;}
+
+ Bool_t IsInCone(const TParticle * p)const ;
+ Bool_t IsInCone(const Double_t e,const Double_t eta,const Double_t phi)const ;
+ Bool_t AcceptConeDeviation(const TParticle *p)const ;
+ Bool_t AcceptConeDeviation(const Double_t e,const Double_t eta,const Double_t phi)const ;
+
+ void SetConeRadius(Double_t r){fConeRad = r ;} ;
+ void SetMaxConeMove(Double_t max = 0.15){fMaxConeMove = max ;} ;
+ void SetMinConeMove(Double_t min = 0.05){fMinConeMove = min ;} ;
+
+ Double_t Energy(void)const{if(fEnergy) return fEnergy ;
+ else return fSumEnergy ;}
+ Double_t Eta(void)const{if(fEta) return fEta;
+ else return fSumEta/fSumEnergy ;}
+ Double_t Phi(void)const{if(fPhi) return fPhi;
+ else return fSumPhi/fSumEnergy ;}
+ Int_t GetNJetParticles(void)const{return fNpart;}
+
+ void Print(Option_t * option="") ;
+
+private:
+ Int_t fNpart ; //Number of particles in jet
+ TArrayI * fList ; //Indexes of particles in list
+
+ Double_t fConeRad ;
+ Double_t fMaxConeMove ;
+ Double_t fMinConeMove ;
+ Double_t fSumEnergy ; //! Intermediate energy
+ Double_t fSumEta ; //! Intermediate eta
+ Double_t fSumPhi ; //! Intermediate phi
+ Double_t fEnergy ; //Energy of the jet
+ Double_t fEta ; //Eta directtion of the jet
+ Double_t fPhi ; //Phi direction of the jet
+ Double_t fLEnergy ; //Energy of Leading particle of jet
+ Double_t fLEta ; //Eta directtion of Leading particle of the jet
+ Double_t fLPhi ; //Phi direction of leading particles of the jet
+
+ ClassDef(AliPHOSJet,1) // description
+
+};
+
+#endif // ALIJET_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// C++ version of UA2 and/or Lund jet finding algorithm
+// UA1 jet algorithm from LUND JETSET (LUCELL)
+//
+//*-- Author : D.Peressounko after UA1 coll. etc
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+#include "TClonesArray.h"
+// #include "TIter.h"
+#include "TParticle.h"
+#include "TTask.h"
+// --- Standard library ---
+
+// --- AliRoot header files ---
+#include "AliPHOSJet.h"
+#include "AliPHOSDigit.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSJetFinder.h"
+#include "AliPHOSDigitizer.h"
+
+ClassImp(AliPHOSJetFinder)
+
+
+//____________________________________________________________________________
+ AliPHOSJetFinder::AliPHOSJetFinder():TNamed("AliPHOSJetFinder","")
+{
+ fNJets = 0 ;
+ fMode = 0 ; //No iterations
+ fStatusCode = -999 ; //no selection
+
+ fConeRad = 1.; //Radius of jet value?????????
+ fMaxConeMove = 0.15 ; //value???????
+ fMinConeMove = 0.05 ; //value???????
+ fEtSeed = 4. ; //Energy of seed particle value??????????
+ fEtMin = 5.; //minimal energy of jet value??????????
+ fPrecBg = 0.00035 ; //value????????
+ fSimGain = 0.;
+ fSimPedestal = 0.;
+
+ fParticles = 0;
+ fJets = 0 ;
+
+}
+
+//____________________________________________________________________________
+ AliPHOSJetFinder::~AliPHOSJetFinder()
+{
+ if(fParticles){
+ delete fParticles ;
+ fParticles = 0 ;
+ }
+ if(fJets){
+ delete fJets ;
+ fJets = 0 ;
+ }
+}
+
+//____________________________________________________________________________
+void AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist)
+{
+
+ TIter next(plist) ;
+
+ TIter nextJet(jetslist) ;
+
+ fNJets = 0 ;
+ TParticle * p ;
+ AliPHOSJet * jet ;
+ //In this cicle we find number of jets and define approx. their directions
+ //note, that we do not really add particles to jet (index =-1)
+ while((p=static_cast<TParticle*>(next()))){
+ if(fStatusCode==-999 || p->GetStatusCode()==fStatusCode){
+ if(p->Energy() >= fEtSeed){ //Energetic enough
+ //cout << "p " << p->Energy() << endl ;
+ //cout << "Status "<<fStatusCode<<" "<<p->GetName()<<" " << p->Energy() << " "<<p->Eta()<< " "<<p->Phi()<<endl ;
+ Bool_t startnew = kTRUE ;
+ //Do not start new jet if part of older jet
+ nextJet.Reset() ;
+ while((jet=static_cast<AliPHOSJet*>(nextJet()))){
+ //jet->Print() ;
+ if(jet->AcceptConeDeviation(p)){
+ startnew = kFALSE ;
+ //cout << "false" << endl ;
+ break ;
+ }
+ }
+ if(startnew){
+ //cout << "new " << endl ;
+ jet = new AliPHOSJet() ;
+ jetslist->Add(jet) ;
+ // jet = static_cast<AliPHOSJet*>(jetslist->Last()) ;
+ jet->SetConeRadius(fConeRad) ;
+ jet->SetMaxConeMove(fMaxConeMove) ;
+ // jet->SetMinConeMove(fMinConeMove) ;
+ jet->AddParticle(p,-1) ;
+ fNJets++;
+ }
+ }
+ while((jet=static_cast<AliPHOSJet*>(nextJet()))){
+ if(jet->AcceptConeDeviation(p))
+ jet->AddParticle(p,-1) ; //Just recalculate direction of jet
+ }
+ }
+ }
+
+ //now calculate directions of jets using collected information
+ nextJet.Reset() ;
+ while((jet=static_cast<AliPHOSJet*>(nextJet()))){
+ jet->CalculateAll() ;
+ if(jet->Energy() < fEtMin){
+ jetslist->Remove(jet) ;
+ delete jet ;
+ }
+ }
+
+ jetslist->Compress() ;
+ //And finally, really add particles to jets
+ for(Int_t iPart=0; iPart<plist->GetEntries();iPart++){
+ p=static_cast<TParticle*>(plist->At(iPart)) ;
+ if(fStatusCode == -999 || p->GetStatusCode()==fStatusCode){
+ Double_t dist = 999999. ; //big distance
+ Int_t iJet = -1 ;
+ for(Int_t i=0; i<jetslist->GetEntriesFast();i++){
+ jet=static_cast<AliPHOSJet*>(jetslist->At(i)) ;
+ if(jet->IsInCone(p)){
+ Double_t cdist = jet->DistanceToJet(p);
+ if(cdist < dist){
+ dist = cdist ;
+ iJet = i ;
+ }
+ }
+ }
+ if(iJet>-1)
+ (static_cast<AliPHOSJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
+ }
+ }
+
+ //Calculate jet parameters
+ nextJet.Reset() ;
+ while((jet=static_cast<AliPHOSJet*>(nextJet()))){
+ jet->CalculateAll() ;
+ }
+
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::FindJetsFromDigits(const TClonesArray * digits, TObjArray * jets){
+
+ if(digits->GetEntries()==0){
+ Error("JetsFromDigits","No entries in digits list \n") ;
+ return ;
+ }
+
+
+ TClonesArray * copyDigits = new TClonesArray(*digits) ;
+
+ //Remove CPV digits if any
+ AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
+ Int_t iDigit ;
+ AliPHOSDigit * digit ;
+ for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
+ if(!geom->IsInEMC(digit->GetId()))
+ copyDigits->RemoveAt(iDigit) ;
+ else
+ break ;
+ }
+ copyDigits->Compress() ;
+
+ Double_t totalEnergy = 0 ;
+ Float_t energy[copyDigits->GetEntries()] ;
+ //calculate average energy of digits
+ //fill array of energies
+ for(iDigit=0;iDigit<copyDigits->GetEntries();iDigit++){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
+ energy[iDigit] = Calibrate(digit) ;
+ totalEnergy+=energy[iDigit] ;
+ }
+
+ //Sort digits in decreasing energy.
+ Int_t index[copyDigits->GetEntries()] ;
+ TMath::Sort(copyDigits->GetEntries(),energy,index) ;
+
+ Double_t eAverage = totalEnergy/copyDigits->GetEntries() ;
+ //remove digits below average energy
+ for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(index[iDigit])) ;
+ if(energy[index[iDigit]] < eAverage)
+ copyDigits->RemoveAt(iDigit) ;
+ else
+ break ;
+ }
+
+
+ AliPHOSJet * jet ;
+ Int_t iIter = 0 ;
+ while(iIter < 10){//less than 10 iterations
+
+ //while digits above seed
+ for(Int_t ind=0;ind<copyDigits->GetEntriesFast();ind++){
+ digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
+ if(energy[index[ind]] > fEtSeed && digit){ //start new jet
+ jet = new AliPHOSJet() ;
+ Double_t e,eta,phi ;
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ jet->AddDigit(e,eta,phi,-1) ;
+ //loop over left digits
+ for(Int_t iDigit = 0 ; iDigit < copyDigits->GetEntries() ; iDigit++){
+ if(iDigit!= ind){ //first digit already in jet
+ digit = static_cast<AliPHOSDigit *>(copyDigits->At(iDigit));
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ if(jet->IsInCone(e,eta,phi) && //is cell in cone
+ jet->AcceptConeDeviation(e,eta,phi)){//if cone does not move too much
+ jet->AddDigit(e,eta,phi,-1) ; //accept new direction
+ }
+ }
+ }//end of loop over cells
+
+ //accept all anused cells incide cone
+ //note, that digits might be returned as anused later
+ for(Int_t icell = 0 ; icell < copyDigits->GetEntries() ; icell++){
+ digit = static_cast<AliPHOSDigit *>(copyDigits->At(icell));
+ if(jet->IsInCone(e,eta,phi)){ //is cell in cone
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ jet->AddDigit(e,eta,phi,digit->GetIndexInList()) ;
+ }
+ }
+
+ //Accept Jet with Et > Et_min and remove all belonging digits
+ if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
+ Int_t nIndxs ;
+ const Int_t * indxs = jet->Indexs(nIndxs) ;
+ for(Int_t i=0;i<nIndxs;i++){
+ copyDigits->RemoveAt(indxs[i]) ;
+ }
+ jet->CalculateAll() ;
+ jets->AddAt(jet,fNJets++);
+ }
+ else{ //remove jet and do not touch digits
+ delete jet ;
+ }
+ }
+ else{
+ if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
+ break ;
+ }
+ }
+
+ iIter++ ;
+ //calculate new energy of backrgound
+ Double_t oldTotalEnergy = totalEnergy ;
+ totalEnergy = 0 ;
+ for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
+ digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
+ if(digit)
+ totalEnergy+=energy[i] ;
+ }
+ if(!fMode || (oldTotalEnergy != 0) &&
+ (TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg))
+ break ;
+ }
+ }
+
+ copyDigits->Delete() ;
+
+}
+//____________________________________________________________________________
+Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
+// if(fPedestals || fGains ){ //use calibration data
+// if(!fPedestals || !fGains ){
+// Error("Calibrate","Either Pedestals of Gains not set!") ;
+// return 0 ;
+// }
+// Float_t en=(digit->GetAmp() - fPedestals->Data(digit->GetId)()))*fGains->Data(digit->GetId()) ;
+// if(en>0)
+// return en ;
+// else
+// return 0. ;
+// }
+// else{ //simulation
+ if(fSimGain==0){ //read simulation parameters
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ if(!gime){
+ Error("Calibrate","Can not read Calibration parameters") ;
+ return 0 ;
+ }
+ const TTask * task = gime->Digitizer() ;
+ if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
+ const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
+ fSimGain = dig->GetEMCchannel() ;
+ fSimPedestal = dig->GetEMCpedestal();
+ }
+ }
+
+ return fSimPedestal + digit->GetAmp()*fSimGain ;
+
+ // }
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
+ e=Calibrate(d) ;
+ AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
+ TVector3 pos ;
+ geom->RelPosInAlice(d->GetId(), pos) ;
+ eta = pos.Eta() ;
+ phi = pos.Phi() ;
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::Print(Option_t * option){
+ printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
+ printf(" Jets found .........%d \n",fNJets) ;
+ printf(" Seed energy cut ....%f \n",fEtSeed) ;
+ printf(" Cone radius ........%f \n",fConeRad) ;
+ printf(" Minimal cone move ..%f \n",fMinConeMove) ;
+ printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
+ printf("------------------------------------------------- \n") ;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// C++ version of UA2 and/or Lund jet finding algorithm
+// UA1 jet algorithm from LUND JETSET (LUCELL)
+//
+//*-- Author : D.Peressounko after UA1 coll. etc
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+#include "TClonesArray.h"
+// #include "TIter.h"
+#include "TParticle.h"
+#include "TTask.h"
+// --- Standard library ---
+
+// --- AliRoot header files ---
+#include "AliJet.h"
+#include "AliPHOSDigit.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSJetFinder.h"
+#include "AliPHOSDigitizer.h"
+
+ClassImp(AliPHOSJetFinder)
+
+
+//____________________________________________________________________________
+ AliPHOSJetFinder::AliPHOSJetFinder():TNamed("AliPHOSJetFinder","")
+{
+ fNJets = 0 ;
+ fMode = 0 ; //No iterations
+ fStatusCode = -999 ; //no selection
+
+ fConeRad = 1.; //Radius of jet value?????????
+ fMaxConeMove = 0.15 ; //value???????
+ fMinConeMove = 0.05 ; //value???????
+ fEtSeed = 4. ; //Energy of seed particle value??????????
+ fEtMin = 5.; //minimal energy of jet value??????????
+ fPrecBg = 0.00035 ; //value????????
+ fSimGain = 0.;
+ fSimPedestal = 0.;
+
+ fParticles = 0;
+ fJets = 0 ;
+
+}
+
+//____________________________________________________________________________
+ AliPHOSJetFinder::~AliPHOSJetFinder()
+{
+ if(fParticles){
+ delete fParticles ;
+ fParticles = 0 ;
+ }
+ if(fJets){
+ delete fJets ;
+ fJets = 0 ;
+ }
+}
+
+//____________________________________________________________________________
+void AliPHOSJetFinder::FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist)
+{
+
+ TIter next(plist) ;
+
+ TIter nextJet(jetslist) ;
+
+ fNJets = 0 ;
+ TParticle * p ;
+ AliJet * jet ;
+ //In this cicle we find number of jets and define approx. their directions
+ //note, that we do not really add particles to jet (index =-1)
+ while((p=static_cast<TParticle*>(next()))){
+ if(fStatusCode==-999 || p->GetStatusCode()==fStatusCode){
+ if(p->Energy() >= fEtSeed){ //Energetic enough
+ //cout << "p " << p->Energy() << endl ;
+ //cout << "Status "<<fStatusCode<<" "<<p->GetName()<<" " << p->Energy() << " "<<p->Eta()<< " "<<p->Phi()<<endl ;
+ Bool_t startnew = kTRUE ;
+ //Do not start new jet if part of older jet
+ nextJet.Reset() ;
+ while((jet=static_cast<AliJet*>(nextJet()))){
+ //jet->Print() ;
+ if(jet->AcceptConeDeviation(p)){
+ startnew = kFALSE ;
+ //cout << "false" << endl ;
+ break ;
+ }
+ }
+ if(startnew){
+ //cout << "new " << endl ;
+ jet = new AliJet() ;
+ jetslist->Add(jet) ;
+ // jet = static_cast<AliJet*>(jetslist->Last()) ;
+ jet->SetConeRadius(fConeRad) ;
+ jet->SetMaxConeMove(fMaxConeMove) ;
+ // jet->SetMinConeMove(fMinConeMove) ;
+ jet->AddParticle(p,-1) ;
+ fNJets++;
+ }
+ }
+ while((jet=static_cast<AliJet*>(nextJet()))){
+ if(jet->AcceptConeDeviation(p))
+ jet->AddParticle(p,-1) ; //Just recalculate direction of jet
+ }
+ }
+ }
+
+ //now calculate directions of jets using collected information
+ nextJet.Reset() ;
+ while((jet=static_cast<AliJet*>(nextJet()))){
+ jet->CalculateAll() ;
+ if(jet->Energy() < fEtMin){
+ jetslist->Remove(jet) ;
+ delete jet ;
+ }
+ }
+
+ jetslist->Compress() ;
+ //And finally, really add particles to jets
+ for(Int_t iPart=0; iPart<plist->GetEntries();iPart++){
+ p=static_cast<TParticle*>(plist->At(iPart)) ;
+ if(fStatusCode == -999 || p->GetStatusCode()==fStatusCode){
+ Double_t dist = 999999. ; //big distance
+ Int_t iJet = -1 ;
+ for(Int_t i=0; i<jetslist->GetEntriesFast();i++){
+ jet=static_cast<AliJet*>(jetslist->At(i)) ;
+ if(jet->IsInCone(p)){
+ Double_t cdist = jet->DistanceToJet(p);
+ if(cdist < dist){
+ dist = cdist ;
+ iJet = i ;
+ }
+ }
+ }
+ if(iJet>-1)
+ (static_cast<AliJet*>(jetslist->At(iJet)))->AddParticle(p,iPart); //assign particle to closest jet
+ }
+ }
+
+ //Calculate jet parameters
+ nextJet.Reset() ;
+ while((jet=static_cast<AliJet*>(nextJet()))){
+ jet->CalculateAll() ;
+ }
+
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::FindJetsFromDigits(const TClonesArray * digits, TObjArray * jets){
+
+ if(digits->GetEntries()==0){
+ Error("JetsFromDigits","No entries in digits list \n") ;
+ return ;
+ }
+
+
+ TClonesArray * copyDigits = new TClonesArray(*digits) ;
+
+ //Remove CPV digits if any
+ AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
+ Int_t iDigit ;
+ AliPHOSDigit * digit ;
+ for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
+ if(!geom->IsInEMC(digit->GetId()))
+ copyDigits->RemoveAt(iDigit) ;
+ else
+ break ;
+ }
+ copyDigits->Compress() ;
+
+ Double_t totalEnergy = 0 ;
+ Float_t energy[copyDigits->GetEntries()] ;
+ //calculate average energy of digits
+ //fill array of energies
+ for(iDigit=0;iDigit<copyDigits->GetEntries();iDigit++){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(iDigit)) ;
+ energy[iDigit] = Calibrate(digit) ;
+ totalEnergy+=energy[iDigit] ;
+ }
+
+ //Sort digits in decreasing energy.
+ Int_t index[copyDigits->GetEntries()] ;
+ TMath::Sort(copyDigits->GetEntries(),energy,index) ;
+
+ Double_t eAverage = totalEnergy/copyDigits->GetEntries() ;
+ //remove digits below average energy
+ for(iDigit=copyDigits->GetEntries()-1;iDigit>=0;iDigit--){
+ digit=static_cast<AliPHOSDigit *>(copyDigits->At(index[iDigit])) ;
+ if(energy[index[iDigit]] < eAverage)
+ copyDigits->RemoveAt(iDigit) ;
+ else
+ break ;
+ }
+
+
+ AliJet * jet ;
+ Int_t iIter = 0 ;
+ while(iIter < 10){//less than 10 iterations
+
+ //while digits above seed
+ for(Int_t ind=0;ind<copyDigits->GetEntriesFast();ind++){
+ digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
+ if(energy[index[ind]] > fEtSeed && digit){ //start new jet
+ jet = new AliJet() ;
+ Double_t e,eta,phi ;
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ jet->AddDigit(e,eta,phi,-1) ;
+ //loop over left digits
+ for(Int_t iDigit = 0 ; iDigit < copyDigits->GetEntries() ; iDigit++){
+ if(iDigit!= ind){ //first digit already in jet
+ digit = static_cast<AliPHOSDigit *>(copyDigits->At(iDigit));
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ if(jet->IsInCone(e,eta,phi) && //is cell in cone
+ jet->AcceptConeDeviation(e,eta,phi)){//if cone does not move too much
+ jet->AddDigit(e,eta,phi,-1) ; //accept new direction
+ }
+ }
+ }//end of loop over cells
+
+ //accept all anused cells incide cone
+ //note, that digits might be returned as anused later
+ for(Int_t icell = 0 ; icell < copyDigits->GetEntries() ; icell++){
+ digit = static_cast<AliPHOSDigit *>(copyDigits->At(icell));
+ if(jet->IsInCone(e,eta,phi)){ //is cell in cone
+ CalculateEEtaPhi(digit,e,eta,phi) ;
+ jet->AddDigit(e,eta,phi,digit->GetIndexInList()) ;
+ }
+ }
+
+ //Accept Jet with Et > Et_min and remove all belonging digits
+ if(jet->Energy()/TMath::CosH(jet->Eta()) > fEtMin){
+ Int_t nIndxs ;
+ const Int_t * indxs = jet->Indexs(nIndxs) ;
+ for(Int_t i=0;i<nIndxs;i++){
+ copyDigits->RemoveAt(indxs[i]) ;
+ }
+ jet->CalculateAll() ;
+ jets->AddAt(jet,fNJets++);
+ }
+ else{ //remove jet and do not touch digits
+ delete jet ;
+ }
+ }
+ else{
+ if(energy[index[ind]] < fEtSeed){ // no more digits above threshold left, return from loop
+ break ;
+ }
+ }
+
+ iIter++ ;
+ //calculate new energy of backrgound
+ Double_t oldTotalEnergy = totalEnergy ;
+ totalEnergy = 0 ;
+ for(Int_t i=0 ; i<copyDigits->GetEntriesFast() ; i++){
+ digit=static_cast<AliPHOSDigit*>(copyDigits->At(index[ind])) ;
+ if(digit)
+ totalEnergy+=energy[i] ;
+ }
+ if(!fMode || (oldTotalEnergy != 0) &&
+ (TMath::Abs(oldTotalEnergy - totalEnergy)/oldTotalEnergy < fPrecBg))
+ break ;
+ }
+ }
+
+ copyDigits->Delete() ;
+
+}
+//____________________________________________________________________________
+Double_t AliPHOSJetFinder::Calibrate(const AliPHOSDigit * digit){
+// if(fPedestals || fGains ){ //use calibration data
+// if(!fPedestals || !fGains ){
+// Error("Calibrate","Either Pedestals of Gains not set!") ;
+// return 0 ;
+// }
+// Float_t en=(digit->GetAmp() - fPedestals->Data(digit->GetId)()))*fGains->Data(digit->GetId()) ;
+// if(en>0)
+// return en ;
+// else
+// return 0. ;
+// }
+// else{ //simulation
+ if(fSimGain==0){ //read simulation parameters
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ if(!gime){
+ Error("Calibrate","Can not read Calibration parameters") ;
+ return 0 ;
+ }
+ const TTask * task = gime->Digitizer() ;
+ if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
+ const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
+ fSimGain = dig->GetEMCchannel() ;
+ fSimPedestal = dig->GetEMCpedestal();
+ }
+ }
+
+ return fSimPedestal + digit->GetAmp()*fSimGain ;
+
+ // }
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &eta, Double_t &phi){
+ e=Calibrate(d) ;
+ AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance("GPS2","") ;
+ TVector3 pos ;
+ geom->RelPosInAlice(d->GetId(), pos) ;
+ eta = pos.Eta() ;
+ phi = pos.Phi() ;
+}
+//____________________________________________________________________________
+void AliPHOSJetFinder::Print(Option_t * option){
+ printf("\n --------------- AliPHOSJetFinder --------------- \n") ;
+ printf(" Jets found .........%d \n",fNJets) ;
+ printf(" Seed energy cut ....%f \n",fEtSeed) ;
+ printf(" Cone radius ........%f \n",fConeRad) ;
+ printf(" Minimal cone move ..%f \n",fMinConeMove) ;
+ printf(" Maximal cone move ..%f \n",fMaxConeMove) ;
+ printf("------------------------------------------------- \n") ;
+}
--- /dev/null
+#ifndef ALIPHOSJETFINDER_H
+#define ALIPHOSJETFINDER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+//
+//*-- Author: D.Peressounko
+
+
+// --- ROOT system ---
+#include "TNamed.h"
+class TClonesArray ;
+class TObjArray ;
+class AliPHOSDigit ;
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+class AliPHOSJetFinder : public TNamed {
+
+public:
+ AliPHOSJetFinder() ; // ctor
+
+ virtual ~AliPHOSJetFinder() ; // dtor
+
+ void FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist) ; //Do the job
+ void FindJetsFromDigits(const TClonesArray * digits,TObjArray * jetslist) ; //Do the job
+
+ void Print(Option_t * option = "") ;
+
+ void SetEtSeed(Double_t etseed){fEtSeed = etseed ;} ;
+ void SetEtMin(Double_t etmin){fEtMin = etmin ;} ;
+ void SetConRad(Double_t r){fConeRad = r ;} ;
+ void SetMaxConeMove(Double_t move){fMaxConeMove=move ; } ;
+ void SetMinConeMove(Double_t move){fMinConeMove=move ; } ;
+ void SetStatusCode(Int_t stc = 1){fStatusCode=stc ;} ;
+
+private:
+ Double_t Calibrate(const AliPHOSDigit * digit) ;
+ void CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &Eta, Double_t &phi);
+
+private:
+ Int_t fNJets ;
+ Int_t fStatusCode ; //Status code of particles to handle
+ Int_t fMode ; //Mode for background calculation
+
+ Double_t fConeRad ; //Maximal radius of jet
+ Double_t fMaxConeMove ;
+ Double_t fMinConeMove ;
+ Double_t fEtSeed ;
+ Double_t fEtMin ;
+ Double_t fPrecBg ;
+ Double_t fSimGain ;
+ Double_t fSimPedestal ;
+
+
+ TClonesArray * fParticles ;
+ TObjArray * fJets ;
+
+ ClassDef(AliPHOSJetFinder,1) //Class to find Jets
+
+};
+
+#endif // AliPHOSJETFINDER_H
--- /dev/null
+#ifndef ALIPHOSJETFINDER_H
+#define ALIPHOSJETFINDER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+//
+//*-- Author: D.Peressounko
+
+
+// --- ROOT system ---
+#include "TNamed.h"
+class TClonesArray ;
+class TObjArray ;
+class AliPHOSDigit ;
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+class AliPHOSJetFinder : public TNamed {
+
+public:
+ AliPHOSJetFinder() ; // ctor
+
+ virtual ~AliPHOSJetFinder() ; // dtor
+
+ void FindJetsFromParticles(const TClonesArray * plist,TObjArray * jetslist) ; //Do the job
+ void FindJetsFromDigits(const TClonesArray * digits,TObjArray * jetslist) ; //Do the job
+
+ void Print(Option_t * option = "") ;
+
+ void SetEtSeed(Double_t etseed){fEtSeed = etseed ;} ;
+ void SetEtMin(Double_t etmin){fEtMin = etmin ;} ;
+ void SetConRad(Double_t r){fConeRad = r ;} ;
+ void SetMaxConeMove(Double_t move){fMaxConeMove=move ; } ;
+ void SetMinConeMove(Double_t move){fMinConeMove=move ; } ;
+ void SetStatusCode(Int_t stc = 1){fStatusCode=stc ;} ;
+
+private:
+ Double_t Calibrate(const AliPHOSDigit * digit) ;
+ void CalculateEEtaPhi(const AliPHOSDigit * d,Double_t &e, Double_t &Eta, Double_t &phi);
+
+private:
+ Int_t fNJets ;
+ Int_t fStatusCode ; //Status code of particles to handle
+ Int_t fMode ; //Mode for background calculation
+
+ Double_t fConeRad ; //Maximal radius of jet
+ Double_t fMaxConeMove ;
+ Double_t fMinConeMove ;
+ Double_t fEtSeed ;
+ Double_t fEtMin ;
+ Double_t fPrecBg ;
+ Double_t fSimGain ;
+ Double_t fSimPedestal ;
+
+
+ TClonesArray * fParticles ;
+ TObjArray * fJets ;
+
+ ClassDef(AliPHOSJetFinder,1) //Class to find Jets
+
+};
+
+#endif // AliPHOSJETFINDER_H
#pragma link C++ class AliPHOSConTableDB+;
#pragma link C++ class AliPHOSCalibrator+;
#pragma link C++ class AliPHOSCalibrationData+;
+#pragma link C++ class AliPHOSJetFinder+;
+#pragma link C++ class AliPHOSJet+;
+
#endif
AliPHOSRaw2Digits.cxx AliPHOSBeamTestEvent.cxx \
AliPHOSCalibrManager.cxx \
AliPHOSConTableDB.cxx AliPHOSCalibrator.cxx\
- AliPHOSCalibrationData.cxx
+ AliPHOSCalibrationData.cxx AliPHOSJet.cxx AliPHOSJetFinder.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff