//-----------------------------------------------------------------------------------
#include <TMath.h>
-#include <TRandom.h>
+#include <TRandom2.h>
#include <TClonesArray.h>
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
fNloopMax(100),
fBeta(0.1),
fNclustMax(0),
- fNin(0)
+ fNin(0),
+ fNeff(0)
{
// Constructor
}
// Find the jets in current event
//
Float_t betaStop=100.;
- fDebug = fHeader->GetDebug();
+ 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();
+ 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) return;
-
+ 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
NumCl(nc,nk,vPy,mPyx,mY);
}while((fBeta<betaStop || nc<4) && nc<fNclustMax);
- Int_t *xx=new Int_t[fNin];
+ 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;
{
//Initialise the variables used by the algorithm
fBeta=0.1;
- fNclustMax= ((AliDAJetHeader*)fHeader)->GetFixedCl() ?
+ fNclustMax = ((AliDAJetHeader*)fHeader)->GetFixedCl() ?
((AliDAJetHeader*)fHeader)->GetNclustMax() :
TMath::Max((Int_t)TMath::Sqrt(fNin),5);
+ Float_t etaEff = ((AliDAJetHeader*)fHeader)->GetEtaEff();
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];
+
+ 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(fNin);
+ vPx->ResizeTo(fNeff);
Int_t iIn=0;
for (Int_t iEn=0; iEn<nEntr; iEn++){
if (fReader->GetCutFlag(iEn)==0) continue;
dEtSum+=(*vPx)(iIn);
iIn++;
}
- for (iIn=0; iIn<fNin; iIn++) (*vPx)(iIn)=(*vPx)(iIn)/dEtSum;
+ 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(fNin,njdim);
+ 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<fNin; iIn++){
+ 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]);
}
//-----------------------------------------------------------------------------------
-void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk, TVectorD *vPy, TMatrixD *mY) const
+void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk, TVectorD *vPy, TMatrixD *mY) const
{
+// Return double clusters
for(Int_t iClust=0; iClust<nc; iClust++){
(*vPy)(iClust)=(*vPy)(iClust)/2;
(*vPy)(nc+iClust)=(*vPy)(iClust);
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 *xEta = xData[0];
+ Double_t *xPhi = xData[1];
Double_t Dist(TVectorD,TVectorD);
Double_t df[2]={fReader->GetReaderHeader()->GetFiducialEtaMax(),pi};
do{
//recalculate conditional probabilities
nloop++;
- for (Int_t iIn=0; iIn<fNin; iIn++){
+ 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));
//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)+=(*vPx)(iIn)*(*mPyx)(iIn,iClust);
- for (Int_t iIn=0; iIn<fNin; iIn++){
+ 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]);
//-----------------------------------------------------------------------------------
void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy, TMatrixD *mPyx,TMatrixD *mY)
{
+ // Number of clusters
static Bool_t growcl=true;
if (nk==2) growcl=true;
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++) {
+ 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++){
ReduceClusters(iSame,nk,nc,cont,nSame);
if (nc >= fNclustMax) growcl=false;
//recalculate the nc distinct codevectors
- TMatrixD *pyx = new TMatrixD(fNin,2*nc);
+ 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<fNin; iIn++) (*pyx)(iIn,iClust)+=(*mPyx)(iIn,iClust1);
+ 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);
}
delete [] nSame;
if (nc > nk/2){
for (Int_t iClust=0; iClust<nc; iClust++){
- for (Int_t iIn=0; iIn<fNin; iIn++) (*mPyx)(iIn,iClust)=(*pyx)(iIn,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);
}
//-----------------------------------------------------------------------------------
void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut) const
{
+// Reduction step
Int_t *nSame = new Int_t[nc];
Int_t *iperm = new Int_t[nc];
Int_t *go = new Int_t[nc];
}
//-----------------------------------------------------------------------------------
-void AliDAJetFinder::EndDetAnn(Int_t &nk,Double_t **xData,Int_t *xx,Double_t etx,
- TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
+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 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;
+ delete [] clusters;
//recalculate codevectors, having all p(y|x)=0 or 1
- Double_t *xEta = xData[0];
- Double_t *xPhi = xData[1];
+ Double_t *xEta = xData[0];
+ Double_t *xPhi = xData[1];
mY->Zero();
mPyx->Zero();
vPy->Zero();
}
//-----------------------------------------------------------------------------------
-void AliDAJetFinder::StoreJets(Int_t nk,Double_t **xData, Int_t *xx, TMatrixD *mY)
+void AliDAJetFinder::StoreJets(Int_t nk, Double_t **xData, Int_t *xx, TMatrixD *mY)
{
//evaluate significant clusters properties
const Double_t pi=TMath::Pi();
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];
+ Float_t dFidEtaMax = rHeader->GetFiducialEtaMax();
+ Float_t dFidEtaMin = rHeader->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++){ //clusters loop
+ 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<fNin; iIn++){
- if (xx[iIn]!=iClust) continue;
+ 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 (((AliDAJetHeader*)fHeader)->GetSelJets()){
for (Int_t iClust=0; iClust<nk; iClust++){
- if (!isJet[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){
+ if(iClust1!=iClust && (*mY)(0,iClust)<dFidEtaMax && (*mY)(0,iClust)>dFidEtaMin){
etDensMed+=etDens[iClust1];
etDensSqr+=TMath::Power(etDens[iClust1],2);
norm++;
}
}
} else {
- for (Int_t iClust=0; iClust<nk; iClust++) isJet[iClust]=true;
+ 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;
+ TMath::Sort(nk,etNoBg,iSort,true);
+ Int_t iCl = 0;
TRefArray *refs = 0;
Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
if (fromAod) refs = fReader->GetReferences();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff