JETAN/DEV/AliDAJetFinder.cxx

   1
   2 //  **************************************************************************
   3 //  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4 //  *                                                                        *
   5 //  * Author: The ALICE Off-line Project.                                    *
   6 //  * Contributors are mentioned in the code where appropriate.              *
   7 //  *                                                                        *
   8 //  * Permission to use, copy, modify and distribute this software and its   *
   9 //  * documentation strictly for non-commercial purposes is hereby granted   *
  10 //  * without fee, provided that the above copyright notice appears in all   *
  11 //  * copies and that both the copyright notice and this permission notice   *
  12 //  * appear in the supporting documentation. The authors make no claims     *
  13 //  * about the suitability of this software for any purpose. It is          *
  14 //  * provided "as is" without express or implied warranty.                  *
  15 //  **************************************************************************
  16
  17 /* $Id$ */
  18
  19 //-----------------------------------------------------------------------------------
  20 // Jet finder based on Deterministic Annealing
  21 // For further informations about the DA working features see:
  22 // Phys.Lett. B601 (2004) 56-63 (http://arxiv.org/abs/hep-ph/0407214)
  23 // Author: Davide Perrino (davide.perrino@ba.infn.it, davide.perrino@cern.ch)
  24 //         alexandre.shabetai@cern.ch (Modification of the input object (reader/finder splitting))
  25 // ** 2011
  26 // Modified accordingly to reader/finder splitting and new handling of neutral information
  27 //-----------------------------------------------------------------------------------
  28
  29 #include <TMath.h>
  30 #include <TRandom2.h>
  31
  32 #include "AliAODJet.h"
  33 #include "AliDAJetHeader.h"
  34 #include "AliDAJetFinder.h"
  35
  36 ClassImp(AliDAJetFinder)
  37
  38 ///////////////////////////////////////////////////////////////////////
  39
  40 AliDAJetFinder::AliDAJetFinder():
  41   AliJetFinder(),
  42   fAlpha(1.01),
  43   fDelta(1e-8),
  44   fAvDist(1e-6),
  45   fEps(1e-4),
  46   fEpsMax(1e-2),
  47   fNloopMax(100),
  48   fBeta(0.1),
  49   fNclustMax(0),
  50   fNin(0),
  51   fNeff(0)
  52 {
  53   // Constructor
  54 }
  55
  56 //-----------------------------------------------------------------------------------
  57 AliDAJetFinder::~AliDAJetFinder()
  58 {
  59   // Destructor
  60 }
  61
  62 //-----------------------------------------------------------------------------------
  63 void AliDAJetFinder::FindJets()
  64 {
  65   // Find the jets in current event
  66   //
  67   Float_t betaStop=100.;
  68   fDebug = fHeader->GetDebug();
  69
  70   Double_t dEtSum=0;
  71   Double_t *xData[2];
  72   TVectorD *vPx = new TVectorD();
  73   TVectorD *vPy = new TVectorD();
  74   TMatrixD *mPyx= new TMatrixD();
  75   TMatrixD *mY  = new TMatrixD();
  76   InitDetAnn(dEtSum,xData,vPx,vPy,mPyx,mY);
  77   if (fNin < fNclustMax){
  78     delete [] xData[0], delete [] xData[1];
  79     delete vPx;
  80     delete vPy;
  81     delete mPyx;
  82     delete mY;
  83     return;
  84   }
  85   Int_t nc=1, nk;
  86   DoubleClusters(nc,nk,vPy,mY);
  87   do{                                   //loop over beta
  88     fBeta*=fAlpha;
  89     Annealing(nk,xData,vPx,vPy,mPyx,mY);
  90     NumCl(nc,nk,vPy,mPyx,mY);
  91   }while((fBeta<betaStop || nc<4) && nc<fNclustMax);
  92
  93   Int_t *xx=new Int_t[fNeff];
  94   for (Int_t i = 0; i < fNeff; i++) xx[i] = 0;
  95
  96   EndDetAnn(nk,xData,xx,dEtSum,vPx,vPy,mPyx,mY);
  97   StoreJets(nk,xData,xx,mY);
  98   delete [] xx;
  99
 100   delete [] xData[0], delete [] xData[1];
 101   delete mPyx;
 102   delete mY;
 103   delete vPx;
 104   delete vPy;
 105
 106 }
 107
 108 //-----------------------------------------------------------------------------------
 109 void AliDAJetFinder::InitDetAnn(Double_t &dEtSum,Double_t **xData,TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
 110 {
 111   //Initialise the variables used by the algorithm
 112   fBeta=0.1;
 113   fNclustMax = ((AliDAJetHeader*)fHeader)->GetFixedCl() ?
 114     ((AliDAJetHeader*)fHeader)->GetNclustMax() :
 115     TMath::Max((Int_t)TMath::Sqrt(fNin),5);
 116   Float_t etaEff = ((AliDAJetHeader*)fHeader)->GetEtaEff();
 117
 118   fNin=0;
 119   for (Int_t iTr=0; iTr<GetCalTrkEvent()->GetNCalTrkTracks(); iTr++) if (GetCalTrkEvent()->GetCalTrkTrack(iTr)->GetCutFlag()==1) fNin++;
 120
 121   fNeff = ((AliDAJetHeader*)fHeader)->GetNeff();
 122   fNeff = TMath::Max(fNeff,fNin);
 123   Double_t *xEta = new Double_t[fNeff];
 124   Double_t *xPhi = new Double_t[fNeff];
 125   xData[0]=xEta; xData[1]=xPhi;
 126   vPx->ResizeTo(fNeff);
 127   Int_t iIn=0;
 128
 129   for (Int_t iTr=0; iTr<GetCalTrkEvent()->GetNCalTrkTracks(); iTr++){
 130     AliJetCalTrkTrack* ctT = GetCalTrkEvent()->GetCalTrkTrack(iTr);
 131     if (ctT->GetCutFlag()==0) continue;
 132     xEta[iIn] = ctT->GetEta();
 133     xPhi[iIn] = ctT->GetPhi()<0 ? ctT->GetPhi() + 2*TMath::Pi() : ctT->GetPhi();
 134     (*vPx)(iIn)=ctT->GetPt();
 135     dEtSum+=(*vPx)(iIn);
 136     iIn++;
 137   }
 138
 139   TRandom2 r;
 140   r.SetSeed(0);
 141   for (iIn=fNin; iIn<fNeff; iIn++){
 142     xEta[iIn]=r.Uniform(-1*etaEff,etaEff);
 143     xPhi[iIn]=r.Uniform(0.,2*TMath::Pi());
 144     (*vPx)(iIn)=r.Uniform(0.01,0.02);
 145     dEtSum+=(*vPx)(iIn);
 146   }
 147   for (iIn=0; iIn<fNeff; iIn++) (*vPx)(iIn)=(*vPx)(iIn)/dEtSum;
 148
 149   Int_t njdim=2*fNclustMax+1;
 150   mPyx->ResizeTo(fNeff,njdim);
 151   mY->ResizeTo(4,njdim);
 152   vPy->ResizeTo(njdim);
 153   mY->Zero();mPyx->Zero();vPy->Zero();
 154   (*vPy)(0)=1;
 155   TMatrixDColumn(*mPyx,0)=1;
 156   Double_t ypos=0,xpos=0;
 157   for (iIn=0; iIn<fNeff; iIn++){
 158     (*mY)(0,0)+=(*vPx)(iIn)*xEta[iIn];
 159     ypos+=(*vPx)(iIn)*TMath::Sin(xPhi[iIn]);
 160     xpos+=(*vPx)(iIn)*TMath::Cos(xPhi[iIn]);
 161   }
 162   (*mY)(1,0)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
 163
 164 }
 165
 166 //-----------------------------------------------------------------------------------
 167 void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk, TVectorD *vPy, TMatrixD *mY) const
 168 {
 169   // Return double clusters
 170   for(Int_t iClust=0; iClust<nc; iClust++){
 171     (*vPy)(iClust)=(*vPy)(iClust)/2;
 172     (*vPy)(nc+iClust)=(*vPy)(iClust);
 173     for(Int_t iComp=0; iComp<3; iComp++) (*mY)(iComp,nc+iClust)=(*mY)(iComp,iClust);
 174   }
 175   nk=2*nc;
 176
 177 }
 178
 179 //-----------------------------------------------------------------------------------
 180 void AliDAJetFinder::Annealing(Int_t nk,Double_t **xData,  const TVectorD *vPx,  TVectorD *vPy,  TMatrixD *mPyx,  TMatrixD *mY)
 181 {
 182   // Main part of the algorithm
 183   const Double_t pi=TMath::Pi();
 184   TVectorD *py = new TVectorD(nk);
 185   TVectorD *p  = new TVectorD(nk);
 186   TMatrixD *y  = new TMatrixD(4,nk);
 187   TMatrixD *y1 = new TMatrixD(4,nk);
 188   TMatrixD *ry = new TMatrixD(2,nk);
 189   Double_t *xEta = xData[0];
 190   Double_t *xPhi = xData[1];
 191   Double_t Dist(TVectorD,TVectorD);
 192
 193   Double_t df[2]={((AliDAJetHeader*)fHeader)->GetFiducialEtaMax(),pi};
 194   TVectorD vPart(2);
 195   Double_t *m = new Double_t[nk];
 196   Double_t chi,chi1;
 197   do{
 198     Int_t nloop=0;
 199     for (Int_t iClust=0; iClust<nk; iClust++){
 200       for (Int_t i=0; i<3; i++)(*y1)(i,iClust)=(*mY)(i,iClust);
 201       (*py)(iClust)=(*vPy)(iClust);
 202     }
 203     //perturbation of codevectors
 204     Double_t seed=1000000*gRandom->Rndm(24);
 205     ry->Randomize(-0.5,0.5,seed);
 206     for (Int_t i=0; i<2; i++){
 207       for (Int_t iClust=0; iClust<nk/2; iClust++)
 208         (*y1)(i,iClust)+=((*ry)(i,iClust)+TMath::Sign(0.5,(*ry)(i,iClust)))*fDelta*df[i];
 209       for (Int_t iClust=nk/2; iClust<nk; iClust++)
 210         (*y1)(i,iClust)-=((*ry)(i,iClust-nk/2)+TMath::Sign(0.5,(*ry)(i,iClust-nk/2)))*fDelta*df[i];
 211     }
 212     do{
 213       //recalculate conditional probabilities
 214       nloop++;
 215       for (Int_t iIn=0; iIn<fNeff; iIn++){
 216         vPart(0)=xEta[iIn]; vPart(1)=xPhi[iIn];
 217         for(Int_t iClust=0; iClust<nk; iClust++){
 218           (*mPyx)(iIn,iClust)=-log((*py)(iClust))+fBeta*Dist(vPart,TMatrixDColumn(*y1,iClust));
 219           m[iClust]=(*mPyx)(iIn,iClust);
 220         }
 221         Double_t pyxNorm=0;
 222         Double_t minPyx=TMath::MinElement(nk,m);
 223         for (Int_t iClust=0; iClust<nk; iClust++){
 224           (*mPyx)(iIn,iClust)-=minPyx;
 225           (*mPyx)(iIn,iClust)=exp(-(*mPyx)(iIn,iClust));
 226           pyxNorm+=(*mPyx)(iIn,iClust);
 227         }
 228         for (Int_t iClust=0; iClust<nk; iClust++) (*mPyx)(iIn,iClust)/=pyxNorm;
 229       }
 230       p->Zero();
 231       y->Zero();
 232       //recalculate codevectors
 233       for (Int_t iClust=0; iClust<nk; iClust++){
 234         Double_t xpos=0,ypos=0,pxy;
 235         for (Int_t iIn=0; iIn<fNeff; iIn++) (*p)(iClust)+=(*vPx)(iIn)*(*mPyx)(iIn,iClust);
 236         for (Int_t iIn=0; iIn<fNeff; iIn++){
 237           pxy=(*vPx)(iIn)*(*mPyx)(iIn,iClust)/(*p)(iClust);
 238           ypos+=pxy*TMath::Sin(xPhi[iIn]);
 239           xpos+=pxy*TMath::Cos(xPhi[iIn]);
 240           (*y)(0,iClust)+=pxy*xEta[iIn];
 241         }
 242         (*y)(1,iClust)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*pi;
 243       }
 244       //verify codevectors' stability
 245       chi=0;
 246       for (Int_t iClust=0; iClust<nk; iClust++){
 247         chi1=TMath::CosH((*y1)(0,iClust)-(*y)(0,iClust))-TMath::Cos((*y1)(1,iClust)-(*y)(1,iClust));
 248         chi1/=(2*TMath::CosH((*y1)(0,iClust))*TMath::CosH((*y)(0,iClust)));
 249         chi1*=chi1;
 250         if (chi1>chi) chi=chi1;
 251       }
 252       chi=TMath::Sqrt(chi);
 253       for (Int_t iClust=0; iClust<nk; iClust++){
 254         for (Int_t i=0; i<2; i++) (*y1)(i,iClust)=(*y)(i,iClust);
 255         (*py)(iClust)=(*p)(iClust);
 256       }
 257       if (nloop>fNloopMax){
 258         if (chi<fEpsMax || nloop>500) break;
 259       }
 260     }while (chi>fEps);
 261   }while (chi>fEpsMax);
 262   for (Int_t iClust=0; iClust<nk; iClust++){                            //set codevectors and probability equal to those calculated
 263     for (Int_t i=0; i<2; i++) (*mY)(i,iClust)=(*y)(i,iClust);
 264     (*vPy)(iClust)=(*p)(iClust);
 265   }
 266   delete py;
 267   delete p;
 268   delete y;
 269   delete y1;
 270   delete ry;
 271   delete [] m;
 272
 273 }
 274
 275 //-----------------------------------------------------------------------------------
 276 void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy, TMatrixD *mPyx,TMatrixD *mY)
 277 {
 278   // Number of clusters
 279   static Bool_t growcl=true;
 280
 281   if (nk==2) growcl=true;
 282   if (growcl){
 283     //verify if two codevectors are equal within fAvDist
 284     Int_t *nSame = new Int_t[nk];
 285     Int_t **iSame = new Int_t*[nk];
 286     Int_t **cont = new Int_t*[nk];
 287     for (Int_t iClust=0; iClust<nk; iClust++) {
 288       cont[iClust] =new Int_t[nk];
 289       iSame[iClust]=new Int_t[nk];
 290     }
 291
 292     for (Int_t iClust=0; iClust<nk; iClust++){
 293       iSame[iClust][iClust]=1;
 294       for (Int_t iClust1=iClust+1; iClust1<nk; iClust1++){
 295         Double_t eta  = (*mY)(0,iClust) ; Double_t phi  = (*mY)(1,iClust);
 296         Double_t eta1 = (*mY)(0,iClust1); Double_t phi1 = (*mY)(1,iClust1);
 297         Double_t distCl=(TMath::CosH(eta-eta1)-TMath::Cos(phi-phi1))/(2*TMath::CosH(eta)*TMath::CosH(eta1));
 298         if (distCl < fAvDist) iSame[iClust][iClust1]=iSame[iClust1][iClust]=1;
 299         else iSame[iClust][iClust1]=iSame[iClust1][iClust]=0;
 300       }
 301     }
 302     ReduceClusters(iSame,nk,nc,cont,nSame);
 303     if (nc >= fNclustMax) growcl=false;
 304     //recalculate the nc distinct codevectors
 305     TMatrixD *pyx = new TMatrixD(fNeff,2*nc);
 306     TVectorD *py = new TVectorD(nk);
 307     TMatrixD *y1  = new TMatrixD(3,nk);
 308     for (Int_t iClust=0; iClust<nc; iClust++){
 309       for(Int_t j=0; j<nSame[iClust]; j++){
 310         Int_t iClust1 = cont[iClust][j];
 311         for (Int_t iIn=0; iIn<fNeff; iIn++) (*pyx)(iIn,iClust)+=(*mPyx)(iIn,iClust1);
 312         (*py)(iClust)+=(*vPy)(iClust1);
 313         for (Int_t i=0; i<2; i++) (*y1)(i,iClust)+=(*mY)(i,iClust1);
 314       }
 315       for (Int_t i=0; i<2; i++) (*y1)(i,iClust)/=nSame[iClust];
 316     }
 317     for (Int_t iClust=0; iClust<nk; iClust++) delete [] cont[iClust], delete [] iSame[iClust];
 318     delete [] iSame;
 319     delete [] cont;
 320     delete [] nSame;
 321     if (nc > nk/2){
 322       for (Int_t iClust=0; iClust<nc; iClust++){
 323         for (Int_t iIn=0; iIn<fNeff; iIn++) (*mPyx)(iIn,iClust)=(*pyx)(iIn,iClust);
 324         for (Int_t iComp=0; iComp<2; iComp++) (*mY)(iComp,iClust)=(*y1)(iComp,iClust);
 325         (*vPy)(iClust)=(*py)(iClust);
 326       }
 327       nk=nc;
 328       if (growcl) DoubleClusters(nc,nk,vPy,mY);
 329     }
 330     delete pyx;
 331     delete py;
 332     delete y1;
 333   }
 334
 335 }
 336
 337 //-----------------------------------------------------------------------------------
 338 void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut) const
 339 {
 340   // Reduction step
 341   Int_t *nSame = new Int_t[nc];
 342   Int_t *iperm = new Int_t[nc];
 343   Int_t *go = new Int_t[nc];
 344   for (Int_t iCl=0; iCl<nc; iCl++){
 345     nSame[iCl]=0;
 346     for (Int_t jCl=0; jCl<nc; jCl++) nSame[iCl]+=iSame[iCl][jCl], cont[iCl][jCl]=0;
 347     iperm[iCl]=iCl;
 348     go[iCl]=1;
 349   }
 350   TMath::Sort(nc,nSame,iperm,true);
 351   Int_t l=0;
 352   for (Int_t iCl=0; iCl<nc; iCl++){
 353     Int_t k=iperm[iCl];
 354     if (go[k] == 1){
 355       Int_t m=0;
 356       for (Int_t jCl=0; jCl<nc; jCl++){
 357         if (iSame[k][jCl] == 1){
 358           cont[l][m]=jCl;
 359           go[jCl]=0;
 360           m++;
 361         }
 362       }
 363       nSameOut[l]=m;
 364       l++;
 365     }
 366   }
 367   ncout=l;
 368   delete [] nSame;
 369   delete [] iperm;
 370   delete [] go;
 371
 372 }
 373
 374 //-----------------------------------------------------------------------------------
 375 void AliDAJetFinder::EndDetAnn(Int_t &nk,Double_t **xData,Int_t *xx,Double_t etx,const TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
 376 {
 377   //now assign each particle to only one cluster
 378   Double_t *clusters=new Double_t[nk];
 379   for (Int_t iIn=0; iIn<fNeff; iIn++){
 380     for (Int_t iClust=0; iClust<nk; iClust++) clusters[iClust]=(*mPyx)(iIn,iClust);
 381     xx[iIn]=TMath::LocMax(nk,clusters);
 382   }
 383   delete [] clusters;
 384
 385   //recalculate codevectors, having all p(y|x)=0 or 1
 386   Double_t *xEta = xData[0];
 387   Double_t *xPhi = xData[1];
 388   mY->Zero();
 389   mPyx->Zero();
 390   vPy->Zero();
 391   for (Int_t iIn=0; iIn<fNin; iIn++){
 392     Int_t iClust=xx[iIn];
 393     (*mPyx)(iIn,iClust)=1;
 394     (*vPy)(iClust)+=(*vPx)(iIn);
 395     (*mY)(0,iClust)+=(*vPx)(iIn)*xEta[iIn];
 396     (*mY)(3,iClust)+=(*vPx)(iIn)*etx;
 397   }
 398   Int_t k=0;
 399   for (Int_t iClust=0; iClust<nk; iClust++){
 400     if ((*vPy)(iClust)>0){
 401       Double_t xpos=0,ypos=0,pxy;
 402       for (Int_t iIn=0; iIn<fNin; iIn++){
 403         pxy=(*vPx)(iIn)*(*mPyx)(iIn,iClust)/(*vPy)(iClust);
 404         ypos+=pxy*TMath::Sin(xPhi[iIn]);
 405         xpos+=pxy*TMath::Cos(xPhi[iIn]);
 406         if (xx[iIn]==iClust) xx[iIn]=k;
 407       }
 408       (*mY)(0,k)=(*mY)(0,iClust)/(*vPy)(iClust);
 409       (*mY)(1,k)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
 410       (*mY)(3,k)=(*mY)(3,iClust);
 411       k++;
 412     }
 413   }
 414   nk=k;
 415
 416 }
 417
 418 //-----------------------------------------------------------------------------------
 419 void AliDAJetFinder::StoreJets(Int_t nk, Double_t **xData, const Int_t *xx, const TMatrixD *mY)
 420 {
 421   //evaluate significant clusters properties
 422   const Double_t pi=TMath::Pi();
 423   Float_t dFidEtaMax = ((AliDAJetHeader*)fHeader)->GetFiducialEtaMax();
 424   Float_t dFidEtaMin = ((AliDAJetHeader*)fHeader)->GetFiducialEtaMin();
 425   Float_t dFiducialEta= dFidEtaMax - dFidEtaMin;
 426   Double_t *xEta = xData[0];
 427   Double_t *xPhi = xData[1];
 428   Int_t nEff = 0;
 429   for (Int_t i=0; i<fNeff; i++) if (xEta[i]<dFidEtaMax && xEta[i]>dFidEtaMin) nEff++;
 430   Double_t dMeanDist=0.;
 431   if (nEff > 0)
 432     dMeanDist=TMath::Sqrt(2*dFiducialEta*pi/nEff);
 433   Bool_t   *isJet = new Bool_t[nk];
 434   Double_t *etNoBg= new Double_t[nk];
 435   Double_t *dDeltaEta=new Double_t[nk];
 436   Double_t *dDeltaPhi=new Double_t[nk];
 437   Double_t *surf  = new Double_t[nk];
 438   Double_t *etDens= new Double_t[nk];
 439   for (Int_t iClust=0; iClust<nk; iClust++){
 440     isJet[iClust]=false;
 441     Double_t dEtaMin=10.,dEtaMax=-10.,dPhiMin=10.,dPhiMax=0.;
 442     for (Int_t iIn=0; iIn<fNeff; iIn++){
 443       if (xx[iIn]!=iClust || xEta[iIn]>dFidEtaMax || xEta[iIn]<dFidEtaMin) continue;
 444       if (xEta[iIn] < dEtaMin) dEtaMin=xEta[iIn];
 445       if (xEta[iIn] > dEtaMax) dEtaMax=xEta[iIn];
 446       Double_t dPhi=xPhi[iIn]-(*mY)(1,iClust);
 447       if      (dPhi > pi     ) dPhi-=2*pi;
 448       else if (dPhi < (-1)*pi) dPhi+=2*pi;
 449       if      (dPhi < dPhiMin) dPhiMin=dPhi;
 450       else if (dPhi > dPhiMax) dPhiMax=dPhi;
 451     }
 452     dDeltaEta[iClust]=dEtaMax-dEtaMin+dMeanDist;
 453     dDeltaPhi[iClust]=dPhiMax-dPhiMin+dMeanDist;
 454     surf[iClust]=0.25*pi*dDeltaEta[iClust]*dDeltaPhi[iClust];
 455     etDens[iClust]=(*mY)(3,iClust)/surf[iClust];
 456   }
 457
 458   if (((AliDAJetHeader*)fHeader)->GetSelJets()){
 459     for (Int_t iClust=0; iClust<nk; iClust++){
 460       if (!isJet[iClust] && (*mY)(0,iClust)<dFidEtaMax && (*mY)(0,iClust)>dFidEtaMin){
 461         Double_t etDensMed=0.;
 462         Double_t etDensSqr=0.;
 463         Int_t norm=0;
 464         for (Int_t iClust1=0; iClust1<nk; iClust1++){
 465           if(iClust1!=iClust && (*mY)(0,iClust)<dFidEtaMax && (*mY)(0,iClust)>dFidEtaMin){
 466             etDensMed+=etDens[iClust1];
 467             etDensSqr+=TMath::Power(etDens[iClust1],2);
 468             norm++;
 469           }
 470         }
 471         etDensMed/=TMath::Max(norm,1);
 472         etDensSqr/=TMath::Max(norm,1);
 473         Double_t deltaEtDens=TMath::Sqrt(etDensSqr-TMath::Power(etDensMed,2));
 474         if ((*mY)(3,iClust) > (etDensMed+deltaEtDens)*surf[iClust]) isJet[iClust]=kTRUE;
 475         etNoBg[iClust]=(*mY)(3,iClust)-etDensMed*surf[iClust];
 476       }
 477     }
 478     for (Int_t iClust=0; iClust<nk; iClust++){
 479       if (isJet[iClust]){
 480         Double_t etDensMed=0;
 481         Double_t extSurf=2*dFiducialEta*pi;
 482         for (Int_t iClust1=0; iClust1<nk; iClust1++){
 483           if (!isJet[iClust1]) etDensMed+=(*mY)(3,iClust1);
 484           else extSurf-=surf[iClust1];
 485         }
 486         etDensMed/=extSurf;
 487         etNoBg[iClust]=(*mY)(3,iClust)-etDensMed*surf[iClust];
 488         if (etNoBg[iClust]<((AliDAJetHeader*)fHeader)->GetEtMin()){
 489           isJet[iClust]=kFALSE;
 490           iClust=-1;
 491         }
 492       }
 493     }
 494   } else {
 495     for (Int_t iClust=0; iClust<nk; iClust++){
 496       isJet[iClust]=true;
 497       etNoBg[iClust]=(*mY)(3,iClust);
 498     }
 499   }
 500   delete [] etDens;
 501   delete [] surf;
 502
 503   //now add selected jets to the list
 504   Int_t *iSort = new Int_t[nk];
 505   TMath::Sort(nk,etNoBg,iSort,true);
 506   Int_t iCl = 0;
 507
 508   for (Int_t iClust=0; iClust<nk; iClust++){ //clusters loop
 509     iCl=iSort[iClust];
 510     if (isJet[iCl]){ //choose cluster
 511       Float_t px,py,pz,en;
 512       px = (*mY)(3,iCl)*TMath::Cos((*mY)(1,iCl));
 513       py = (*mY)(3,iCl)*TMath::Sin((*mY)(1,iCl));
 514       pz = (*mY)(3,iCl)/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-(*mY)(0,iCl))));
 515       en = TMath::Sqrt(px * px + py * py + pz * pz);
 516       AliAODJet jet(px, py, pz, en);
 517       Int_t iIn=0;
 518       Int_t nTr = GetCalTrkEvent()->GetNCalTrkTracks();
 519       for (Int_t iTr=0; iTr<nTr; iTr++){
 520         AliJetCalTrkTrack* ctT = GetCalTrkEvent()->GetCalTrkTrack(iTr);
 521         if (ctT->GetCutFlag()==0) continue;
 522         if (xx[iIn]==iCl) jet.AddTrack(ctT->GetTrackObject());
 523         iIn++;
 524       }
 525       AddJet(jet);
 526       if (fDebug > 0) printf("jet %d, Eta: %f, Phi: %f, Et: %f\n",iCl,jet.Eta(),jet.Phi(),jet.Pt());
 527     }
 528   }
 529   delete [] dDeltaEta; delete [] dDeltaPhi;
 530   delete [] etNoBg;
 531   delete [] isJet;
 532   delete [] iSort;
 533
 534 }
 535
 536 //-----------------------------------------------------------------------------------
 537 Double_t Dist(TVectorD x,TVectorD y)
 538 {
 539   // Squared distance
 540   const Double_t pi=TMath::Pi();
 541   Double_t dphi=TMath::Abs(x(1)-y(1));
 542   if (dphi > pi) dphi=2*pi-dphi;
 543   Double_t dist=pow(x(0)-y(0),2)+pow(dphi,2);
 544   return dist;
 545
 546 }