X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliDAJetFinder.cxx;h=3bd12ee2422b970cca71807697edc4cedb3f68f2;hb=48f32a76d93bab329a856ef9883dde81fd415c38;hp=8fad2a45c3d3eac4102fec0c700887970b83e95b;hpb=852db00e23c8c716f89942536743989de41a5388;p=u%2Fmrichter%2FAliRoot.git

diff --git a/JETAN/AliDAJetFinder.cxx b/JETAN/AliDAJetFinder.cxx
index 8fad2a45c3d..3bd12ee2422 100644
--- a/JETAN/AliDAJetFinder.cxx
+++ b/JETAN/AliDAJetFinder.cxx
@@ -14,464 +14,533 @@
 //  * provided "as is" without express or implied warranty.                  *
 //  **************************************************************************
 
+/* $Id$ */
+
 //-----------------------------------------------------------------------------------
 // Jet finder based on Deterministic Annealing
 // For further informations about the DA working features see:
 // Phys.Lett. B601 (2004) 56-63 (http://arxiv.org/abs/hep-ph/0407214)
 // Author: Davide Perrino (davide.perrino@ba.infn.it, davide.perrino@cern.ch)
+//         alexandre.shabetai@cern.ch (Modification of the input object (reader/finder splitting)) 
+// ** 2011
+// Modified accordingly to reader/finder splitting and new handling of neutral information   
 //-----------------------------------------------------------------------------------
 
 #include <TMath.h>
-#include <TRandom.h>
-#include <TClonesArray.h>
-#include "AliJetReader.h"
+#include <TRandom2.h>
+
+#include "AliAODJet.h"
 #include "AliDAJetHeader.h"
 #include "AliDAJetFinder.h"
 
-
 ClassImp(AliDAJetFinder)
 
+///////////////////////////////////////////////////////////////////////
 
-//-----------------------------------------------------------------------------------
 AliDAJetFinder::AliDAJetFinder():
-	fAlpha(1.01),
-	fDelta(1e-8),
-	fAvDist(1e-6),
-	fEps(1e-4),
-	fEpsMax(1e-2),
-	fNloopMax(100),
-	fBeta(0.1),
-	fNclustMax(0),
-	fPyx(0x0),
-	fY(0x0),
-	fPx(0x0),
-	fPy(0x0),
-	fXEta(0x0),
-	fXPhi(0x0),
-	fNin(0)
+  AliJetFinder(),
+  fAlpha(1.01),
+  fDelta(1e-8),
+  fAvDist(1e-6),
+  fEps(1e-4),
+  fEpsMax(1e-2),
+  fNloopMax(100),
+  fBeta(0.1),
+  fNclustMax(0),
+  fNin(0),
+  fNeff(0)
 {
-	// Constructor
+  // Constructor
 }
 
 //-----------------------------------------------------------------------------------
 AliDAJetFinder::~AliDAJetFinder()
 {
-	// Destructor
-	delete fPyx;
-	delete fY;
-	delete fPx;
-	delete fPy;
-	delete [] fXEta;
-	delete [] fXPhi;
+  // Destructor
 }
 
 //-----------------------------------------------------------------------------------
 void AliDAJetFinder::FindJets()  
 {
-// Find the jets in current event
-// 
-	Float_t betaStop=100.;
-
-	Double_t dEtSum=0;
-	InitDetAnn(dEtSum);
-	if (!fNin) return;
-
-	Int_t nc=1,nk;
-	DoubleClusters(nc,nk);
-	do{					//loop over beta
-		fBeta*=fAlpha;
-		Annealing(nk);
-		NumCl(nc,nk);
-	}while((fBeta<betaStop || nc<4) && nc<fNclustMax);
-
-	Int_t *xx=new Int_t[fNin];
-	EndDetAnn(nk,xx,dEtSum);
-	StoreJets(nk,xx);
-	delete [] xx;
+  // Find the jets in current event
+  //
+  Float_t betaStop=100.;
+  fDebug = fHeader->GetDebug();
+
+  Double_t dEtSum=0;
+  Double_t *xData[2];
+  TVectorD *vPx = new TVectorD();
+  TVectorD *vPy = new TVectorD();
+  TMatrixD *mPyx= new TMatrixD();
+  TMatrixD *mY  = new TMatrixD();
+  InitDetAnn(dEtSum,xData,vPx,vPy,mPyx,mY);
+  if (fNin < fNclustMax){
+    delete [] xData[0], delete [] xData[1];
+    delete vPx;
+    delete vPy;
+    delete mPyx;
+    delete mY;
+    return;
+  }
+  Int_t nc=1, nk;
+  DoubleClusters(nc,nk,vPy,mY);
+  do{					//loop over beta
+    fBeta*=fAlpha;
+    Annealing(nk,xData,vPx,vPy,mPyx,mY);
+    NumCl(nc,nk,vPy,mPyx,mY);
+  }while((fBeta<betaStop || nc<4) && nc<fNclustMax);
+
+  Int_t *xx=new Int_t[fNeff];
+  for (Int_t i = 0; i < fNeff; i++) xx[i] = 0;
+
+  EndDetAnn(nk,xData,xx,dEtSum,vPx,vPy,mPyx,mY);
+  StoreJets(nk,xData,xx,mY);
+  delete [] xx;
+
+  delete [] xData[0], delete [] xData[1];
+  delete mPyx;
+  delete mY;
+  delete vPx;
+  delete vPy;
 
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::InitDetAnn(Double_t &dEtSum)
+void AliDAJetFinder::InitDetAnn(Double_t &dEtSum,Double_t **xData,TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
 {
-//Initialise the variables used by the algorithm
-	fBeta=0.1;
-	TClonesArray *lvArray = fReader->GetMomentumArray();
-	fNin = lvArray->GetEntries();
-	fNclustMax= ((AliDAJetHeader*)fHeader)->GetFixedCl() ? 
-	    ((AliDAJetHeader*)fHeader)->GetNclustMax() 
-	    : 
-	    TMath::Max((Int_t)TMath::Sqrt(fNin),5);
-	fXEta=new Double_t[fNin]; fXPhi=new Double_t[fNin];
-	fPx = new TVectorD(fNin);
-	for (Int_t iIn=0; iIn<fNin; iIn++){
-		TLorentzVector *lv=(TLorentzVector*)lvArray->At(iIn);
-		fXEta[iIn] = lv->Eta();
-		fXPhi[iIn] = lv->Phi()<0 ? lv->Phi() + 2*TMath::Pi() : lv->Phi();
-		(*fPx)(iIn)=lv->Pt();
-		dEtSum+=(*fPx)(iIn);
-	}
-	for (Int_t iIn=0; iIn<fNin; iIn++) (*fPx)(iIn)=(*fPx)(iIn)/dEtSum;
-
-	Int_t njdim=2*fNclustMax+1;
-	fPyx = new TMatrixD(fNin,njdim);
-	fY = new TMatrixD(4,njdim);
-	fPy= new TVectorD(njdim);
-	fY->Zero();fPyx->Zero();fPy->Zero();
-	(*fPy)(0)=1;
-	TMatrixDColumn(*fPyx,0)=1;
-	Double_t ypos=0,xpos=0;
-	for (Int_t iIn=0; iIn<fNin; iIn++){
-		(*fY)(0,0)+=(*fPx)(iIn)*fXEta[iIn];
-		ypos+=(*fPx)(iIn)*TMath::Sin(fXPhi[iIn]);
-		xpos+=(*fPx)(iIn)*TMath::Cos(fXPhi[iIn]);
-	}
-	(*fY)(1,0)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
-	lvArray->Delete();
+  //Initialise the variables used by the algorithm
+  fBeta=0.1;
+  fNclustMax = ((AliDAJetHeader*)fHeader)->GetFixedCl() ?
+    ((AliDAJetHeader*)fHeader)->GetNclustMax() :
+    TMath::Max((Int_t)TMath::Sqrt(fNin),5);
+  Float_t etaEff = ((AliDAJetHeader*)fHeader)->GetEtaEff();
+
+  fNin=0;
+  for (Int_t iTr=0; iTr<GetCalTrkEvent()->GetNCalTrkTracks(); iTr++) if (GetCalTrkEvent()->GetCalTrkTrack(iTr)->GetCutFlag()==1) fNin++;
+
+  fNeff = ((AliDAJetHeader*)fHeader)->GetNeff();
+  fNeff = TMath::Max(fNeff,fNin);
+  Double_t *xEta = new Double_t[fNeff];
+  Double_t *xPhi = new Double_t[fNeff];
+  xData[0]=xEta; xData[1]=xPhi;
+  vPx->ResizeTo(fNeff);
+  Int_t iIn=0;
+
+  for (Int_t iTr=0; iTr<GetCalTrkEvent()->GetNCalTrkTracks(); iTr++){
+    AliJetCalTrkTrack* ctT = GetCalTrkEvent()->GetCalTrkTrack(iTr);
+    if (ctT->GetCutFlag()==0) continue;
+    xEta[iIn] = ctT->GetEta();
+    xPhi[iIn] = ctT->GetPhi()<0 ? ctT->GetPhi() + 2*TMath::Pi() : ctT->GetPhi();
+    (*vPx)(iIn)=ctT->GetPt();
+    dEtSum+=(*vPx)(iIn);
+    iIn++;
+  }
+
+  TRandom2 r;
+  r.SetSeed(0);
+  for (iIn=fNin; iIn<fNeff; iIn++){
+    xEta[iIn]=r.Uniform(-1*etaEff,etaEff);
+    xPhi[iIn]=r.Uniform(0.,2*TMath::Pi());
+    (*vPx)(iIn)=r.Uniform(0.01,0.02);
+    dEtSum+=(*vPx)(iIn);
+  }
+  for (iIn=0; iIn<fNeff; iIn++) (*vPx)(iIn)=(*vPx)(iIn)/dEtSum;
+
+  Int_t njdim=2*fNclustMax+1;
+  mPyx->ResizeTo(fNeff,njdim);
+  mY->ResizeTo(4,njdim);
+  vPy->ResizeTo(njdim);
+  mY->Zero();mPyx->Zero();vPy->Zero();
+  (*vPy)(0)=1;
+  TMatrixDColumn(*mPyx,0)=1;
+  Double_t ypos=0,xpos=0;
+  for (iIn=0; iIn<fNeff; iIn++){
+    (*mY)(0,0)+=(*vPx)(iIn)*xEta[iIn];
+    ypos+=(*vPx)(iIn)*TMath::Sin(xPhi[iIn]);
+    xpos+=(*vPx)(iIn)*TMath::Cos(xPhi[iIn]);
+  }
+  (*mY)(1,0)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
+
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk)
+void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk, TVectorD *vPy, TMatrixD *mY) const
 {
-	for(Int_t iClust=0; iClust<nc; iClust++){
-		(*fPy)(iClust)=(*fPy)(iClust)/2;
-		(*fPy)(nc+iClust)=(*fPy)(iClust);
-		for(Int_t iComp=0; iComp<3; iComp++) (*fY)(iComp,nc+iClust)=(*fY)(iComp,iClust);
-	}
-	nk=2*nc;
+  // Return double clusters
+  for(Int_t iClust=0; iClust<nc; iClust++){
+    (*vPy)(iClust)=(*vPy)(iClust)/2;
+    (*vPy)(nc+iClust)=(*vPy)(iClust);
+    for(Int_t iComp=0; iComp<3; iComp++) (*mY)(iComp,nc+iClust)=(*mY)(iComp,iClust);
+  }
+  nk=2*nc;
+
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::Annealing(Int_t nk)
+void AliDAJetFinder::Annealing(Int_t nk,Double_t **xData,  const TVectorD *vPx,  TVectorD *vPy,  TMatrixD *mPyx,  TMatrixD *mY)
 {
-// Main part of the algorithm
-	const Double_t pi=TMath::Pi();
-	TVectorD *py = new TVectorD(nk);
-	TVectorD *p  = new TVectorD(nk);
-	TMatrixD *y  = new TMatrixD(4,nk);
-	TMatrixD *y1 = new TMatrixD(4,nk);
-	TMatrixD *ry = new TMatrixD(2,nk);
-	Double_t Dist(TVectorD,TVectorD);
-
-	Double_t df[2]={((AliDAJetHeader*)fHeader)->GetEtaCut(),pi};
-	TVectorD vPart(2);
-	Double_t *m = new Double_t[nk];
-	Double_t chi,chi1;
-	do{
-		Int_t nloop=0;
-		for (Int_t iClust=0; iClust<nk; iClust++){
-			for (Int_t i=0; i<3; i++)(*y1)(i,iClust)=(*fY)(i,iClust);
-			(*py)(iClust)=(*fPy)(iClust);
-		}
-	//perturbation of codevectors
-		Double_t seed=1000000*gRandom->Rndm(24);
-		ry->Randomize(-0.5,0.5,seed);
-		for (Int_t i=0; i<2; i++){
-			for (Int_t iClust=0; iClust<nk/2; iClust++)
-				(*y1)(i,iClust)+=((*ry)(i,iClust)+TMath::Sign(0.5,(*ry)(i,iClust)))*fDelta*df[i];
-			for (Int_t iClust=nk/2; iClust<nk; iClust++)
-				(*y1)(i,iClust)-=((*ry)(i,iClust-nk/2)+TMath::Sign(0.5,(*ry)(i,iClust-nk/2)))*fDelta*df[i];
-		}
-		do{
-	//recalculate conditional probabilities
-			nloop++;
-			for (Int_t iIn=0; iIn<fNin; iIn++){
-				vPart(0)=fXEta[iIn]; vPart(1)=fXPhi[iIn];
-				for(Int_t iClust=0; iClust<nk; iClust++){
-					(*fPyx)(iIn,iClust)=-log((*py)(iClust))+fBeta*Dist(vPart,TMatrixDColumn(*y1,iClust));
-					m[iClust]=(*fPyx)(iIn,iClust);
-				}
-				Double_t pyxNorm=0;
-				Double_t minPyx=TMath::MinElement(nk,m);
-				for (Int_t iClust=0; iClust<nk; iClust++){
-					(*fPyx)(iIn,iClust)-=minPyx;
-					(*fPyx)(iIn,iClust)=exp(-(*fPyx)(iIn,iClust));
-					pyxNorm+=(*fPyx)(iIn,iClust);
-				}
-				for (Int_t iClust=0; iClust<nk; iClust++) (*fPyx)(iIn,iClust)/=pyxNorm;
-			}
-			p->Zero();
-			y->Zero();
-	//recalculate codevectors
-			for (Int_t iClust=0; iClust<nk; iClust++){
-				Double_t xpos=0,ypos=0,pxy;
-				for (Int_t iIn=0; iIn<fNin; iIn++) (*p)(iClust)+=(*fPx)(iIn)*(*fPyx)(iIn,iClust);
-				for (Int_t iIn=0; iIn<fNin; iIn++){
-					pxy=(*fPx)(iIn)*(*fPyx)(iIn,iClust)/(*p)(iClust);
-					ypos+=pxy*TMath::Sin(fXPhi[iIn]);
-					xpos+=pxy*TMath::Cos(fXPhi[iIn]);
-					(*y)(0,iClust)+=pxy*fXEta[iIn];
-				}
-				(*y)(1,iClust)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*pi;
-			}
-	//verify codevectors' stability
-			chi=0;
-			for (Int_t iClust=0; iClust<nk; iClust++){
-				chi1=TMath::CosH((*y1)(0,iClust)-(*y)(0,iClust))-TMath::Cos((*y1)(1,iClust)-(*y)(1,iClust));
-				chi1/=(2*TMath::CosH((*y1)(0,iClust))*TMath::CosH((*y)(0,iClust)));
-				chi1*=chi1;
-				if (chi1>chi) chi=chi1;
-			}
-			chi=TMath::Sqrt(chi);
-			for (Int_t iClust=0; iClust<nk; iClust++){
-				for (Int_t i=0; i<2; i++) (*y1)(i,iClust)=(*y)(i,iClust);
-				(*py)(iClust)=(*p)(iClust);
-			}
-			if (nloop>fNloopMax){
-				if (chi<fEpsMax || nloop>500) break;
-			}
-		}while (chi>fEps);
-	}while (chi>fEpsMax);
-	for (Int_t iClust=0; iClust<nk; iClust++){				//set codevectors and probability equal to those calculated
-		for (Int_t i=0; i<2; i++) (*fY)(i,iClust)=(*y)(i,iClust);
-		(*fPy)(iClust)=(*p)(iClust);
+  // Main part of the algorithm
+  const Double_t pi=TMath::Pi();
+  TVectorD *py = new TVectorD(nk);
+  TVectorD *p  = new TVectorD(nk);
+  TMatrixD *y  = new TMatrixD(4,nk);
+  TMatrixD *y1 = new TMatrixD(4,nk);
+  TMatrixD *ry = new TMatrixD(2,nk);
+  Double_t *xEta = xData[0];
+  Double_t *xPhi = xData[1];
+  Double_t Dist(TVectorD,TVectorD);
+
+  Double_t df[2]={((AliDAJetHeader*)fHeader)->GetFiducialEtaMax(),pi};
+  TVectorD vPart(2);
+  Double_t *m = new Double_t[nk];
+  Double_t chi,chi1;
+  do{
+    Int_t nloop=0;
+    for (Int_t iClust=0; iClust<nk; iClust++){
+      for (Int_t i=0; i<3; i++)(*y1)(i,iClust)=(*mY)(i,iClust);
+      (*py)(iClust)=(*vPy)(iClust);
+    }
+    //perturbation of codevectors
+    Double_t seed=1000000*gRandom->Rndm(24);
+    ry->Randomize(-0.5,0.5,seed);
+    for (Int_t i=0; i<2; i++){
+      for (Int_t iClust=0; iClust<nk/2; iClust++)
+	(*y1)(i,iClust)+=((*ry)(i,iClust)+TMath::Sign(0.5,(*ry)(i,iClust)))*fDelta*df[i];
+      for (Int_t iClust=nk/2; iClust<nk; iClust++)
+	(*y1)(i,iClust)-=((*ry)(i,iClust-nk/2)+TMath::Sign(0.5,(*ry)(i,iClust-nk/2)))*fDelta*df[i];
+    }
+    do{
+      //recalculate conditional probabilities
+      nloop++;
+      for (Int_t iIn=0; iIn<fNeff; iIn++){
+	vPart(0)=xEta[iIn]; vPart(1)=xPhi[iIn];
+	for(Int_t iClust=0; iClust<nk; iClust++){
+	  (*mPyx)(iIn,iClust)=-log((*py)(iClust))+fBeta*Dist(vPart,TMatrixDColumn(*y1,iClust));
+	  m[iClust]=(*mPyx)(iIn,iClust);
 	}
-    delete py;
-    delete p;
-    delete y;
-    delete y1;
-    delete ry;
-	delete [] m;
+	Double_t pyxNorm=0;
+	Double_t minPyx=TMath::MinElement(nk,m);
+	for (Int_t iClust=0; iClust<nk; iClust++){
+	  (*mPyx)(iIn,iClust)-=minPyx;
+	  (*mPyx)(iIn,iClust)=exp(-(*mPyx)(iIn,iClust));
+	  pyxNorm+=(*mPyx)(iIn,iClust);
+	}
+	for (Int_t iClust=0; iClust<nk; iClust++) (*mPyx)(iIn,iClust)/=pyxNorm;
+      }
+      p->Zero();
+      y->Zero();
+      //recalculate codevectors
+      for (Int_t iClust=0; iClust<nk; iClust++){
+	Double_t xpos=0,ypos=0,pxy;
+	for (Int_t iIn=0; iIn<fNeff; iIn++) (*p)(iClust)+=(*vPx)(iIn)*(*mPyx)(iIn,iClust);
+	for (Int_t iIn=0; iIn<fNeff; iIn++){
+	  pxy=(*vPx)(iIn)*(*mPyx)(iIn,iClust)/(*p)(iClust);
+	  ypos+=pxy*TMath::Sin(xPhi[iIn]);
+	  xpos+=pxy*TMath::Cos(xPhi[iIn]);
+	  (*y)(0,iClust)+=pxy*xEta[iIn];
+	}
+	(*y)(1,iClust)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*pi;
+      }
+      //verify codevectors' stability
+      chi=0;
+      for (Int_t iClust=0; iClust<nk; iClust++){
+	chi1=TMath::CosH((*y1)(0,iClust)-(*y)(0,iClust))-TMath::Cos((*y1)(1,iClust)-(*y)(1,iClust));
+	chi1/=(2*TMath::CosH((*y1)(0,iClust))*TMath::CosH((*y)(0,iClust)));
+	chi1*=chi1;
+	if (chi1>chi) chi=chi1;
+      }
+      chi=TMath::Sqrt(chi);
+      for (Int_t iClust=0; iClust<nk; iClust++){
+	for (Int_t i=0; i<2; i++) (*y1)(i,iClust)=(*y)(i,iClust);
+	(*py)(iClust)=(*p)(iClust);
+      }
+      if (nloop>fNloopMax){
+	if (chi<fEpsMax || nloop>500) break;
+      }
+    }while (chi>fEps);
+  }while (chi>fEpsMax);
+  for (Int_t iClust=0; iClust<nk; iClust++){				//set codevectors and probability equal to those calculated
+    for (Int_t i=0; i<2; i++) (*mY)(i,iClust)=(*y)(i,iClust);
+    (*vPy)(iClust)=(*p)(iClust);
+  }
+  delete py;
+  delete p;
+  delete y;
+  delete y1;
+  delete ry;
+  delete [] m;
+
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk)
+void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy, TMatrixD *mPyx,TMatrixD *mY)
 {
-	static Bool_t growcl=true;
+  // Number of clusters
+  static Bool_t growcl=true;
 	
-	if (nk==2) growcl=true;
-	if (growcl){
-//verify if two codevectors are equal within fAvDist
-		Int_t *nSame = new Int_t[nk];
-		Int_t **iSame = new Int_t*[nk];
-		Int_t **cont = new Int_t*[nk];
-		for (Int_t iClust=0; iClust<nk; iClust++) cont[iClust]=new Int_t[nk],iSame[iClust]=new Int_t[nk];
-		for (Int_t iClust=0; iClust<nk; iClust++){
-			iSame[iClust][iClust]=1;
-			for (Int_t iClust1=iClust+1; iClust1<nk; iClust1++){
-				Double_t eta  = (*fY)(0,iClust) ; Double_t phi  = (*fY)(1,iClust);
-				Double_t eta1 = (*fY)(0,iClust1); Double_t phi1 = (*fY)(1,iClust1);
-				Double_t distCl=(TMath::CosH(eta-eta1)-TMath::Cos(phi-phi1))/(2*TMath::CosH(eta)*TMath::CosH(eta1));
-				if (distCl < fAvDist) iSame[iClust][iClust1]=iSame[iClust1][iClust]=1;
-				else iSame[iClust][iClust1]=iSame[iClust1][iClust]=0;
-			}
-		}
-		ReduceClusters(iSame,nk,nc,cont,nSame);
-		if (nc >= fNclustMax) growcl=false;
-//recalculate the nc distinct codevectors
-		TMatrixD *pyx = new TMatrixD(fNin,2*nc);
-		TVectorD *py = new TVectorD(nk);
-		TMatrixD *y1  = new TMatrixD(3,nk);
-		for (Int_t iClust=0; iClust<nc; iClust++){
-			for(Int_t j=0; j<nSame[iClust]; j++){
-				Int_t iClust1 = cont[iClust][j];
-				for (Int_t iIn=0; iIn<fNin; iIn++) (*pyx)(iIn,iClust)+=(*fPyx)(iIn,iClust1);
-				(*py)(iClust)+=(*fPy)(iClust1);
-				for (Int_t i=0; i<2; i++) (*y1)(i,iClust)+=(*fY)(i,iClust1);
-			}
-			for (Int_t i=0; i<2; i++) (*y1)(i,iClust)/=nSame[iClust];
-		}
-		for (Int_t iClust=0; iClust<nk; iClust++) delete [] cont[iClust], delete [] iSame[iClust];
-		delete [] iSame;
-		delete [] cont;
-		delete [] nSame;
-		if (nc > nk/2){
-			for (Int_t iClust=0; iClust<nc; iClust++){
-				for (Int_t iIn=0; iIn<fNin; iIn++) (*fPyx)(iIn,iClust)=(*pyx)(iIn,iClust);
-				for (Int_t iComp=0; iComp<2; iComp++) (*fY)(iComp,iClust)=(*y1)(iComp,iClust);
-				(*fPy)(iClust)=(*py)(iClust);
-			}
-			nk=nc;
-			if (growcl) DoubleClusters(nc,nk);
-		}
-		delete pyx;
-		delete py;
-		delete y1;
-	}
+  if (nk==2) growcl=true;
+  if (growcl){
+    //verify if two codevectors are equal within fAvDist
+    Int_t *nSame = new Int_t[nk];
+    Int_t **iSame = new Int_t*[nk];
+    Int_t **cont = new Int_t*[nk];
+    for (Int_t iClust=0; iClust<nk; iClust++) {
+      cont[iClust] =new Int_t[nk];
+      iSame[iClust]=new Int_t[nk];         
+    }
+
+    for (Int_t iClust=0; iClust<nk; iClust++){
+      iSame[iClust][iClust]=1;
+      for (Int_t iClust1=iClust+1; iClust1<nk; iClust1++){
+	Double_t eta  = (*mY)(0,iClust) ; Double_t phi  = (*mY)(1,iClust);
+	Double_t eta1 = (*mY)(0,iClust1); Double_t phi1 = (*mY)(1,iClust1);
+	Double_t distCl=(TMath::CosH(eta-eta1)-TMath::Cos(phi-phi1))/(2*TMath::CosH(eta)*TMath::CosH(eta1));
+	if (distCl < fAvDist) iSame[iClust][iClust1]=iSame[iClust1][iClust]=1;
+	else iSame[iClust][iClust1]=iSame[iClust1][iClust]=0;
+      }
+    }
+    ReduceClusters(iSame,nk,nc,cont,nSame);
+    if (nc >= fNclustMax) growcl=false;
+    //recalculate the nc distinct codevectors
+    TMatrixD *pyx = new TMatrixD(fNeff,2*nc);
+    TVectorD *py = new TVectorD(nk);
+    TMatrixD *y1  = new TMatrixD(3,nk);
+    for (Int_t iClust=0; iClust<nc; iClust++){
+      for(Int_t j=0; j<nSame[iClust]; j++){
+	Int_t iClust1 = cont[iClust][j];
+	for (Int_t iIn=0; iIn<fNeff; iIn++) (*pyx)(iIn,iClust)+=(*mPyx)(iIn,iClust1);
+	(*py)(iClust)+=(*vPy)(iClust1);
+	for (Int_t i=0; i<2; i++) (*y1)(i,iClust)+=(*mY)(i,iClust1);
+      }
+      for (Int_t i=0; i<2; i++) (*y1)(i,iClust)/=nSame[iClust];
+    }
+    for (Int_t iClust=0; iClust<nk; iClust++) delete [] cont[iClust], delete [] iSame[iClust];
+    delete [] iSame;
+    delete [] cont;
+    delete [] nSame;
+    if (nc > nk/2){
+      for (Int_t iClust=0; iClust<nc; iClust++){
+	for (Int_t iIn=0; iIn<fNeff; iIn++) (*mPyx)(iIn,iClust)=(*pyx)(iIn,iClust);
+	for (Int_t iComp=0; iComp<2; iComp++) (*mY)(iComp,iClust)=(*y1)(iComp,iClust);
+	(*vPy)(iClust)=(*py)(iClust);
+      }
+      nk=nc;
+      if (growcl) DoubleClusters(nc,nk,vPy,mY);
+    }
+    delete pyx;
+    delete py;
+    delete y1;
+  }
 
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut)
+void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut) const
 {
-	Int_t *nSame = new Int_t[nc];
-	Int_t *iperm = new Int_t[nc];
-	Int_t *go = new Int_t[nc];
-	for (Int_t iCl=0; iCl<nc; iCl++){
-		nSame[iCl]=0;
-		for (Int_t jCl=0; jCl<nc; jCl++) nSame[iCl]+=iSame[iCl][jCl], cont[iCl][jCl]=0;
-		iperm[iCl]=iCl;
-		go[iCl]=1;
-	}
-	TMath::Sort(nc,nSame,iperm,true);
-	Int_t l=0;
-	for (Int_t iCl=0; iCl<nc; iCl++){
-		Int_t k=iperm[iCl];
-		if (go[k] == 1){
-			Int_t m=0;
-			for (Int_t jCl=0; jCl<nc; jCl++){
-				if (iSame[k][jCl] == 1){
-					cont[l][m]=jCl;
-					go[jCl]=0;
-					m++;
-				}
-			}
-			nSameOut[l]=m;
-			l++;
-		}
+  // Reduction step
+  Int_t *nSame = new Int_t[nc];
+  Int_t *iperm = new Int_t[nc];
+  Int_t *go = new Int_t[nc];
+  for (Int_t iCl=0; iCl<nc; iCl++){
+    nSame[iCl]=0;
+    for (Int_t jCl=0; jCl<nc; jCl++) nSame[iCl]+=iSame[iCl][jCl], cont[iCl][jCl]=0;
+    iperm[iCl]=iCl;
+    go[iCl]=1;
+  }
+  TMath::Sort(nc,nSame,iperm,true);
+  Int_t l=0;
+  for (Int_t iCl=0; iCl<nc; iCl++){
+    Int_t k=iperm[iCl];
+    if (go[k] == 1){
+      Int_t m=0;
+      for (Int_t jCl=0; jCl<nc; jCl++){
+	if (iSame[k][jCl] == 1){
+	  cont[l][m]=jCl;
+	  go[jCl]=0;
+	  m++;
 	}
-	ncout=l;
-	delete [] nSame;
-	delete [] iperm;
-	delete [] go;
+      }
+      nSameOut[l]=m;
+      l++;
+    }
+  }
+  ncout=l;
+  delete [] nSame;
+  delete [] iperm;
+  delete [] go;
+
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::EndDetAnn(Int_t &nk,Int_t *xx,Double_t etx)
+void AliDAJetFinder::EndDetAnn(Int_t &nk,Double_t **xData,Int_t *xx,Double_t etx,const TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
 {
-//now assign each particle to only one cluster
-	Double_t *clusters=new Double_t[nk];
-	for (Int_t iIn=0; iIn<fNin; iIn++){
-		for (Int_t iClust=0; iClust<nk; iClust++) clusters[iClust]=(*fPyx)(iIn,iClust);
-		xx[iIn]=TMath::LocMax(nk,clusters);
-	}
-    delete [] clusters;
+  //now assign each particle to only one cluster
+  Double_t *clusters=new Double_t[nk];
+  for (Int_t iIn=0; iIn<fNeff; iIn++){
+    for (Int_t iClust=0; iClust<nk; iClust++) clusters[iClust]=(*mPyx)(iIn,iClust);
+    xx[iIn]=TMath::LocMax(nk,clusters);
+  }
+  delete [] clusters;
 	
-//recalculate codevectors, having all p(y|x)=0 or 1
-	fY->Zero();
-	fPyx->Zero();
-	fPy->Zero();
-	for (Int_t iIn=0; iIn<fNin; iIn++){
-		Int_t iClust=xx[iIn];
-		(*fPyx)(iIn,iClust)=1;
-		(*fPy)(iClust)+=(*fPx)(iIn);
-		(*fY)(0,iClust)+=(*fPx)(iIn)*fXEta[iIn];
-		(*fY)(3,iClust)+=(*fPx)(iIn)*etx;
-	}
-	Int_t k=0;
-	for (Int_t iClust=0; iClust<nk; iClust++){
-		if ((*fPy)(iClust)>0){
-			Double_t xpos=0,ypos=0,pxy;
-			for (Int_t iIn=0; iIn<fNin; iIn++){
-				pxy=(*fPx)(iIn)*(*fPyx)(iIn,iClust)/(*fPy)(iClust);
-				ypos+=pxy*TMath::Sin(fXPhi[iIn]);
-				xpos+=pxy*TMath::Cos(fXPhi[iIn]);
-				if (xx[iIn]==iClust) xx[iIn]=k;
-			}
-			(*fY)(0,k)=(*fY)(0,iClust)/(*fPy)(iClust);
-			(*fY)(1,k)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
-			(*fY)(3,k)=(*fY)(3,iClust);
-			k++;
-		}
-	}
-	nk=k;
+  //recalculate codevectors, having all p(y|x)=0 or 1
+  Double_t *xEta = xData[0];
+  Double_t *xPhi = xData[1];
+  mY->Zero();
+  mPyx->Zero();
+  vPy->Zero();
+  for (Int_t iIn=0; iIn<fNin; iIn++){
+    Int_t iClust=xx[iIn];
+    (*mPyx)(iIn,iClust)=1;
+    (*vPy)(iClust)+=(*vPx)(iIn);
+    (*mY)(0,iClust)+=(*vPx)(iIn)*xEta[iIn];
+    (*mY)(3,iClust)+=(*vPx)(iIn)*etx;
+  }
+  Int_t k=0;
+  for (Int_t iClust=0; iClust<nk; iClust++){
+    if ((*vPy)(iClust)>0){
+      Double_t xpos=0,ypos=0,pxy;
+      for (Int_t iIn=0; iIn<fNin; iIn++){
+	pxy=(*vPx)(iIn)*(*mPyx)(iIn,iClust)/(*vPy)(iClust);
+	ypos+=pxy*TMath::Sin(xPhi[iIn]);
+	xpos+=pxy*TMath::Cos(xPhi[iIn]);
+	if (xx[iIn]==iClust) xx[iIn]=k;
+      }
+      (*mY)(0,k)=(*mY)(0,iClust)/(*vPy)(iClust);
+      (*mY)(1,k)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
+      (*mY)(3,k)=(*mY)(3,iClust);
+      k++;
+    }
+  }
+  nk=k;
+
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
+void AliDAJetFinder::StoreJets(Int_t nk, Double_t **xData, const Int_t *xx, const TMatrixD *mY)
 {
-//evaluate significant clusters properties
-	const Double_t pi=TMath::Pi();
-	Double_t dMeanDist=TMath::Sqrt(4*((AliDAJetHeader*)fHeader)->GetEtaCut()*pi/fNin);
-	Bool_t   *isJet = new Bool_t[nk];
-	Double_t *etNoBg= new Double_t[nk];
-	Double_t *dDeltaEta=new Double_t[nk];
-	Double_t *dDeltaPhi=new Double_t[nk];
-	Double_t *surf  = new Double_t[nk];
-	Double_t *etDens= new Double_t[nk];
-	for (Int_t iClust=0; iClust<nk; iClust++){									//clusters loop
-		isJet[iClust]=false;
-		Double_t dEtaMin=10.,dEtaMax=-10.,dPhiMin=10.,dPhiMax=0.;
-		for (Int_t iIn=0; iIn<fNin; iIn++){
-			if (xx[iIn]!=iClust) continue;
-			if (fXEta[iIn] < dEtaMin) dEtaMin=fXEta[iIn];
-			if (fXEta[iIn] > dEtaMax) dEtaMax=fXEta[iIn];
-			Double_t dPhi=fXPhi[iIn]-(*fY)(1,iClust);
-			if      (dPhi > pi     ) dPhi-=2*pi;
-			else if (dPhi < (-1)*pi) dPhi+=2*pi;
-			if      (dPhi < dPhiMin) dPhiMin=dPhi;
-			else if (dPhi > dPhiMax) dPhiMax=dPhi;
-		}
-		dDeltaEta[iClust]=dEtaMax-dEtaMin+dMeanDist;
-		dDeltaPhi[iClust]=dPhiMax-dPhiMin+dMeanDist;
-		surf[iClust]=0.25*pi*dDeltaEta[iClust]*dDeltaPhi[iClust];
-		etDens[iClust]=(*fY)(3,iClust)/surf[iClust];
-	}
+  //evaluate significant clusters properties
+  const Double_t pi=TMath::Pi();
+  Float_t dFidEtaMax = ((AliDAJetHeader*)fHeader)->GetFiducialEtaMax();
+  Float_t dFidEtaMin = ((AliDAJetHeader*)fHeader)->GetFiducialEtaMin();
+  Float_t dFiducialEta= dFidEtaMax - dFidEtaMin;
+  Double_t *xEta = xData[0];
+  Double_t *xPhi = xData[1];
+  Int_t nEff = 0;
+  for (Int_t i=0; i<fNeff; i++) if (xEta[i]<dFidEtaMax && xEta[i]>dFidEtaMin) nEff++;
+  Double_t dMeanDist=0.;
+  if (nEff > 0)
+    dMeanDist=TMath::Sqrt(2*dFiducialEta*pi/nEff);
+  Bool_t   *isJet = new Bool_t[nk];
+  Double_t *etNoBg= new Double_t[nk];
+  Double_t *dDeltaEta=new Double_t[nk];
+  Double_t *dDeltaPhi=new Double_t[nk];
+  Double_t *surf  = new Double_t[nk];
+  Double_t *etDens= new Double_t[nk];
+  for (Int_t iClust=0; iClust<nk; iClust++){
+    isJet[iClust]=false;
+    Double_t dEtaMin=10.,dEtaMax=-10.,dPhiMin=10.,dPhiMax=0.;
+    for (Int_t iIn=0; iIn<fNeff; iIn++){
+      if (xx[iIn]!=iClust || xEta[iIn]>dFidEtaMax || xEta[iIn]<dFidEtaMin) continue;
+      if (xEta[iIn] < dEtaMin) dEtaMin=xEta[iIn];
+      if (xEta[iIn] > dEtaMax) dEtaMax=xEta[iIn];
+      Double_t dPhi=xPhi[iIn]-(*mY)(1,iClust);
+      if      (dPhi > pi     ) dPhi-=2*pi;
+      else if (dPhi < (-1)*pi) dPhi+=2*pi;
+      if      (dPhi < dPhiMin) dPhiMin=dPhi;
+      else if (dPhi > dPhiMax) dPhiMax=dPhi;
+    }
+    dDeltaEta[iClust]=dEtaMax-dEtaMin+dMeanDist;
+    dDeltaPhi[iClust]=dPhiMax-dPhiMin+dMeanDist;
+    surf[iClust]=0.25*pi*dDeltaEta[iClust]*dDeltaPhi[iClust];
+    etDens[iClust]=(*mY)(3,iClust)/surf[iClust];
+  }
 
-	if (((AliDAJetHeader*)fHeader)->GetSelJets()){
-		for (Int_t iClust=0; iClust<nk; iClust++){
-			if (!isJet[iClust]){
-				Double_t etDensMed=0.;
-				Double_t etDensSqr=0.;
-				Int_t norm=0;
-				for (Int_t iClust1=0; iClust1<nk; iClust1++){
-					if(iClust1!=iClust){
-						etDensMed+=etDens[iClust1];
-						etDensSqr+=TMath::Power(etDens[iClust1],2);
-						norm++;
-					}
-				}
-				etDensMed/=TMath::Max(norm,1);
-				etDensSqr/=TMath::Max(norm,1);
-				Double_t deltaEtDens=TMath::Sqrt(etDensSqr-TMath::Power(etDensMed,2));
-				if ((*fY)(3,iClust) > (etDensMed+deltaEtDens)*surf[iClust]) isJet[iClust]=kTRUE;
-				etNoBg[iClust]=(*fY)(3,iClust)-etDensMed*surf[iClust];
-			}
-		}
-		for (Int_t iClust=0; iClust<nk; iClust++){
-			if (isJet[iClust]){
-				Double_t etDensMed=0;
-				Double_t extSurf=4*((AliDAJetHeader*)fHeader)->GetEtaCut()*pi;
-				for (Int_t iClust1=0; iClust1<nk; iClust1++){
-					if (!isJet[iClust1]) etDensMed+=(*fY)(3,iClust1);
-					else extSurf-=surf[iClust1];
-				}
-				etDensMed/=extSurf;
-				etNoBg[iClust]=(*fY)(3,iClust)-etDensMed*surf[iClust];
-				if (etNoBg[iClust]<((AliDAJetHeader*)fHeader)->GetEtMin()){
-					isJet[iClust]=kFALSE;
-					iClust=-1;
-				}
-			}
-		}
-	} else {
-		for (Int_t iClust=0; iClust<nk; iClust++) isJet[iClust]=true;
+  if (((AliDAJetHeader*)fHeader)->GetSelJets()){
+    for (Int_t iClust=0; iClust<nk; iClust++){
+      if (!isJet[iClust] && (*mY)(0,iClust)<dFidEtaMax && (*mY)(0,iClust)>dFidEtaMin){
+	Double_t etDensMed=0.;
+	Double_t etDensSqr=0.;
+	Int_t norm=0;
+	for (Int_t iClust1=0; iClust1<nk; iClust1++){
+	  if(iClust1!=iClust && (*mY)(0,iClust)<dFidEtaMax && (*mY)(0,iClust)>dFidEtaMin){
+	    etDensMed+=etDens[iClust1];
+	    etDensSqr+=TMath::Power(etDens[iClust1],2);
+	    norm++;
+	  }
 	}
-	delete [] etDens;
-	delete [] surf;
-	
-//now add selected jets to the list
-	Int_t *inJet = new Int_t[fNin];
-	TRefArray *refs = 0;
-	Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
-	if (fromAod) refs = fReader->GetReferences();
-	for (Int_t iClust=0; iClust<nk; iClust++){									//clusters loop
-		if (isJet[iClust]){														//choose cluster
-			Float_t px,py,pz,en;
-			px = (*fY)(3,iClust)*TMath::Cos((*fY)(1,iClust));
-			py = (*fY)(3,iClust)*TMath::Sin((*fY)(1,iClust));
-			pz = (*fY)(3,iClust)/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-(*fY)(0,iClust))));
-			en = TMath::Sqrt(px * px + py * py + pz * pz);
-			AliAODJet jet(px, py, pz, en);
-			if (fromAod) 
-			    for (Int_t iIn=0; iIn<fNin; iIn++) if (xx[iIn]==iClust) jet.AddTrack(refs->At(iIn));
-			AddJet(jet);
-			printf("jet %d, Eta: %f, Phi: %f, Et: %f\n",iClust,jet.Eta(),jet.Phi(),jet.Pt());
-		}
+	etDensMed/=TMath::Max(norm,1);
+	etDensSqr/=TMath::Max(norm,1);
+	Double_t deltaEtDens=TMath::Sqrt(etDensSqr-TMath::Power(etDensMed,2));
+	if ((*mY)(3,iClust) > (etDensMed+deltaEtDens)*surf[iClust]) isJet[iClust]=kTRUE;
+	etNoBg[iClust]=(*mY)(3,iClust)-etDensMed*surf[iClust];
+      }
+    }
+    for (Int_t iClust=0; iClust<nk; iClust++){
+      if (isJet[iClust]){
+	Double_t etDensMed=0;
+	Double_t extSurf=2*dFiducialEta*pi;
+	for (Int_t iClust1=0; iClust1<nk; iClust1++){
+	  if (!isJet[iClust1]) etDensMed+=(*mY)(3,iClust1);
+	  else extSurf-=surf[iClust1];
 	}
-	delete [] dDeltaEta; delete [] dDeltaPhi;
-	delete [] etNoBg;
-	delete [] isJet;
-	delete [] inJet;
+	etDensMed/=extSurf;
+	etNoBg[iClust]=(*mY)(3,iClust)-etDensMed*surf[iClust];
+	if (etNoBg[iClust]<((AliDAJetHeader*)fHeader)->GetEtMin()){
+	  isJet[iClust]=kFALSE;
+	  iClust=-1;
+	}
+      }
+    }
+  } else {
+    for (Int_t iClust=0; iClust<nk; iClust++){
+      isJet[iClust]=true;
+      etNoBg[iClust]=(*mY)(3,iClust);
+    }
+  }
+  delete [] etDens;
+  delete [] surf;
+	
+  //now add selected jets to the list
+  Int_t *iSort = new Int_t[nk];
+  TMath::Sort(nk,etNoBg,iSort,true);
+  Int_t iCl = 0;
+
+  for (Int_t iClust=0; iClust<nk; iClust++){ //clusters loop
+    iCl=iSort[iClust];
+    if (isJet[iCl]){ //choose cluster
+      Float_t px,py,pz,en;
+      px = (*mY)(3,iCl)*TMath::Cos((*mY)(1,iCl));
+      py = (*mY)(3,iCl)*TMath::Sin((*mY)(1,iCl));
+      pz = (*mY)(3,iCl)/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-(*mY)(0,iCl))));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      AliAODJet jet(px, py, pz, en);
+      Int_t iIn=0;
+      Int_t nTr = GetCalTrkEvent()->GetNCalTrkTracks();
+      for (Int_t iTr=0; iTr<nTr; iTr++){
+	AliJetCalTrkTrack* ctT = GetCalTrkEvent()->GetCalTrkTrack(iTr);
+	if (ctT->GetCutFlag()==0) continue;
+	if (xx[iIn]==iCl) jet.AddTrack(ctT->GetTrackObject());
+	iIn++;
+      }
+      AddJet(jet);
+      if (fDebug > 0) printf("jet %d, Eta: %f, Phi: %f, Et: %f\n",iCl,jet.Eta(),jet.Phi(),jet.Pt());
+    }
+  }
+  delete [] dDeltaEta; delete [] dDeltaPhi;
+  delete [] etNoBg;
+  delete [] isJet;
+  delete [] iSort;
+
 }
 
 //-----------------------------------------------------------------------------------
 Double_t Dist(TVectorD x,TVectorD y)
 {
-// Squared distance
-	const Double_t pi=TMath::Pi();
-	Double_t dphi=TMath::Abs(x(1)-y(1));
-	if (dphi > pi) dphi=2*pi-dphi;
-	Double_t dist=pow(x(0)-y(0),2)+pow(dphi,2);
-	return dist;
+  // Squared distance
+  const Double_t pi=TMath::Pi();
+  Double_t dphi=TMath::Abs(x(1)-y(1));
+  if (dphi > pi) dphi=2*pi-dphi;
+  Double_t dist=pow(x(0)-y(0),2)+pow(dphi,2);
+  return dist;
+
 }