]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliDAJetFinder.cxx
avoid compilation warnings when adding selection macros
[u/mrichter/AliRoot.git] / JETAN / AliDAJetFinder.cxx
index 8fad2a45c3d3eac4102fec0c700887970b83e95b..130dc076180f43077e0b290190df8e195ab5c90d 100644 (file)
 #include <TMath.h>
 #include <TRandom.h>
 #include <TClonesArray.h>
+#include "AliJetReaderHeader.h"
 #include "AliJetReader.h"
 #include "AliDAJetHeader.h"
 #include "AliDAJetFinder.h"
 
-
 ClassImp(AliDAJetFinder)
 
 
 //-----------------------------------------------------------------------------------
 AliDAJetFinder::AliDAJetFinder():
+       AliJetFinder(),
        fAlpha(1.01),
        fDelta(1e-8),
        fAvDist(1e-6),
@@ -42,12 +43,6 @@ AliDAJetFinder::AliDAJetFinder():
        fNloopMax(100),
        fBeta(0.1),
        fNclustMax(0),
-       fPyx(0x0),
-       fY(0x0),
-       fPx(0x0),
-       fPy(0x0),
-       fXEta(0x0),
-       fXPhi(0x0),
        fNin(0)
 {
        // Constructor
@@ -57,12 +52,6 @@ AliDAJetFinder::AliDAJetFinder():
 AliDAJetFinder::~AliDAJetFinder()
 {
        // Destructor
-       delete fPyx;
-       delete fY;
-       delete fPx;
-       delete fPy;
-       delete [] fXEta;
-       delete [] fXPhi;
 }
 
 //-----------------------------------------------------------------------------------
@@ -73,76 +62,92 @@ void AliDAJetFinder::FindJets()
        Float_t betaStop=100.;
 
        Double_t dEtSum=0;
-       InitDetAnn(dEtSum);
+       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) return;
 
-       Int_t nc=1,nk;
-       DoubleClusters(nc,nk);
+       Int_t nc=1, nk;
+       DoubleClusters(nc,nk,vPy,mY);
        do{                                     //loop over beta
                fBeta*=fAlpha;
-               Annealing(nk);
-               NumCl(nc,nk);
+               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[fNin];
-       EndDetAnn(nk,xx,dEtSum);
-       StoreJets(nk,xx);
+       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() 
-           : 
+           ((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);
+       TClonesArray *lvArray = fReader->GetMomentumArray();
+       Int_t nEntr = lvArray->GetEntries();
+       fNin=0;
+       for (Int_t iEn=0; iEn<nEntr; iEn++) if (fReader->GetCutFlag(iEn)==1) fNin++;
+       Double_t *xEta = new Double_t[fNin];
+       Double_t *xPhi = new Double_t[fNin];
+       xData[0]=xEta; xData[1]=xPhi;
+       vPx->ResizeTo(fNin);
+       Int_t iIn=0;
+       for (Int_t iEn=0; iEn<nEntr; iEn++){
+               if (fReader->GetCutFlag(iEn)==0) continue;
+               TLorentzVector *lv=(TLorentzVector*)lvArray->At(iEn);
+               xEta[iIn] = lv->Eta();
+               xPhi[iIn] = lv->Phi()<0 ? lv->Phi() + 2*TMath::Pi() : lv->Phi();
+               (*vPx)(iIn)=lv->Pt();
+               dEtSum+=(*vPx)(iIn);
+               iIn++;
        }
-       for (Int_t iIn=0; iIn<fNin; iIn++) (*fPx)(iIn)=(*fPx)(iIn)/dEtSum;
+       for (iIn=0; iIn<fNin; iIn++) (*vPx)(iIn)=(*vPx)(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;
+       mPyx->ResizeTo(fNin,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 (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]);
+       for (iIn=0; iIn<fNin; iIn++){
+               (*mY)(0,0)+=(*vPx)(iIn)*xEta[iIn];
+               ypos+=(*vPx)(iIn)*TMath::Sin(xPhi[iIn]);
+               xpos+=(*vPx)(iIn)*TMath::Cos(xPhi[iIn]);
        }
-       (*fY)(1,0)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
-       lvArray->Delete();
+       (*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);
+               (*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,  TVectorD *vPx,  TVectorD *vPy,  TMatrixD *mPyx,  TMatrixD *mY)
 {
 // Main part of the algorithm
        const Double_t pi=TMath::Pi();
@@ -151,17 +156,19 @@ void AliDAJetFinder::Annealing(Int_t 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)->GetEtaCut(),pi};
+       Double_t df[2]={fReader->GetReaderHeader()->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)=(*fY)(i,iClust);
-                       (*py)(iClust)=(*fPy)(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);
@@ -176,31 +183,31 @@ void AliDAJetFinder::Annealing(Int_t nk)
        //recalculate conditional probabilities
                        nloop++;
                        for (Int_t iIn=0; iIn<fNin; iIn++){
-                               vPart(0)=fXEta[iIn]; vPart(1)=fXPhi[iIn];
+                               vPart(0)=xEta[iIn]; vPart(1)=xPhi[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);
+                                       (*mPyx)(iIn,iClust)=-log((*py)(iClust))+fBeta*Dist(vPart,TMatrixDColumn(*y1,iClust));
+                                       m[iClust]=(*mPyx)(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);
+                                       (*mPyx)(iIn,iClust)-=minPyx;
+                                       (*mPyx)(iIn,iClust)=exp(-(*mPyx)(iIn,iClust));
+                                       pyxNorm+=(*mPyx)(iIn,iClust);
                                }
-                               for (Int_t iClust=0; iClust<nk; iClust++) (*fPyx)(iIn,iClust)/=pyxNorm;
+                               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<fNin; iIn++) (*p)(iClust)+=(*fPx)(iIn)*(*fPyx)(iIn,iClust);
+                               for (Int_t iIn=0; iIn<fNin; iIn++) (*p)(iClust)+=(*vPx)(iIn)*(*mPyx)(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];
+                                       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;
                        }
@@ -223,19 +230,19 @@ void AliDAJetFinder::Annealing(Int_t nk)
                }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);
+               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;
+    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;
        
@@ -249,8 +256,8 @@ void AliDAJetFinder::NumCl(Int_t &nc,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 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;
@@ -265,9 +272,9 @@ void AliDAJetFinder::NumCl(Int_t &nc,Int_t &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 iIn=0; iIn<fNin; 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];
                }
@@ -277,12 +284,12 @@ void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk)
                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);
+                               for (Int_t iIn=0; iIn<fNin; 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);
+                       if (growcl) DoubleClusters(nc,nk,vPy,mY);
                }
                delete pyx;
                delete py;
@@ -292,7 +299,7 @@ void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk)
 }
 
 //-----------------------------------------------------------------------------------
-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];
@@ -327,40 +334,43 @@ void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **
 }
 
 //-----------------------------------------------------------------------------------
-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,
+                               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);
+               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();
+    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];
-               (*fPyx)(iIn,iClust)=1;
-               (*fPy)(iClust)+=(*fPx)(iIn);
-               (*fY)(0,iClust)+=(*fPx)(iIn)*fXEta[iIn];
-               (*fY)(3,iClust)+=(*fPx)(iIn)*etx;
+               (*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 ((*fPy)(iClust)>0){
+               if ((*vPy)(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]);
+                               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;
                        }
-                       (*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);
+                       (*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++;
                }
        }
@@ -368,11 +378,15 @@ void AliDAJetFinder::EndDetAnn(Int_t &nk,Int_t *xx,Double_t etx)
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
+void AliDAJetFinder::StoreJets(Int_t nk,Double_t **xData,  Int_t *xx,  TMatrixD *mY)
 {
 //evaluate significant clusters properties
        const Double_t pi=TMath::Pi();
-       Double_t dMeanDist=TMath::Sqrt(4*((AliDAJetHeader*)fHeader)->GetEtaCut()*pi/fNin);
+       AliJetReaderHeader *rHeader=fReader->GetReaderHeader();
+       Float_t dFiducialEta=rHeader->GetFiducialEtaMax()-rHeader->GetFiducialEtaMin();
+       Double_t dMeanDist=TMath::Sqrt(2*dFiducialEta*pi/fNin);
+    Double_t *xEta = xData[0];
+    Double_t *xPhi = xData[1];
        Bool_t   *isJet = new Bool_t[nk];
        Double_t *etNoBg= new Double_t[nk];
        Double_t *dDeltaEta=new Double_t[nk];
@@ -384,9 +398,9 @@ void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
                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 (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;
@@ -395,7 +409,7 @@ void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
                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];
+               etDens[iClust]=(*mY)(3,iClust)/surf[iClust];
        }
 
        if (((AliDAJetHeader*)fHeader)->GetSelJets()){
@@ -414,20 +428,20 @@ void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
                                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];
+                               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=4*((AliDAJetHeader*)fHeader)->GetEtaCut()*pi;
+                               Double_t extSurf=2*dFiducialEta*pi;
                                for (Int_t iClust1=0; iClust1<nk; iClust1++){
-                                       if (!isJet[iClust1]) etDensMed+=(*fY)(3,iClust1);
+                                       if (!isJet[iClust1]) etDensMed+=(*mY)(3,iClust1);
                                        else extSurf-=surf[iClust1];
                                }
                                etDensMed/=extSurf;
-                               etNoBg[iClust]=(*fY)(3,iClust)-etDensMed*surf[iClust];
+                               etNoBg[iClust]=(*mY)(3,iClust)-etDensMed*surf[iClust];
                                if (etNoBg[iClust]<((AliDAJetHeader*)fHeader)->GetEtMin()){
                                        isJet[iClust]=kFALSE;
                                        iClust=-1;
@@ -448,13 +462,20 @@ void AliDAJetFinder::StoreJets(Int_t nk,Int_t *xx)
        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))));
+                       px = (*mY)(3,iClust)*TMath::Cos((*mY)(1,iClust));
+                       py = (*mY)(3,iClust)*TMath::Sin((*mY)(1,iClust));
+                       pz = (*mY)(3,iClust)/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-(*mY)(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));
+                       if (fromAod){
+                               Int_t iIn=0;
+                               Int_t nEntr = fReader->GetMomentumArray()->GetEntries();
+                               for (Int_t iEn=0; iEn<nEntr; iEn++){
+                                       if (fReader->GetCutFlag(iEn)==0) continue;
+                                       if (xx[iIn]==iClust) jet.AddTrack(refs->At(iEn));
+                                       iIn++;
+                               }
+                       }
                        AddJet(jet);
                        printf("jet %d, Eta: %f, Phi: %f, Et: %f\n",iClust,jet.Eta(),jet.Phi(),jet.Pt());
                }