* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
//---------------------------------------------------------------------
// UA1 Cone Algorithm Jet finder
// (version in c++)
//---------------------------------------------------------------------
-#include <TLorentzVector.h>
+#include <TArrayF.h>
+#include <TClonesArray.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
-#include <TArrayF.h>
+#include <TLorentzVector.h>
+
#include "AliUA1JetFinderV1.h"
#include "AliUA1JetHeaderV1.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
#include "AliJet.h"
+#include "AliAODJet.h"
ClassImp(AliUA1JetFinderV1)
-////////////////////////////////////////////////////////////////////////
-
-AliUA1JetFinderV1::AliUA1JetFinderV1()
+/////////////////////////////////////////////////////////////////////
+AliUA1JetFinderV1::AliUA1JetFinderV1() :
+ AliJetFinder(),
+ fLego(0)
{
// Constructor
- fHeader = 0x0;
- fLego = 0x0;
}
////////////////////////////////////////////////////////////////////////
if (nIn == 0) return;
// local arrays for input
- Float_t* ptT = new Float_t[nIn];
- Float_t* etaT = new Float_t[nIn];
- Float_t* phiT = new Float_t[nIn];
+ Float_t* ptT = new Float_t[nIn];
+ Float_t* etaT = new Float_t[nIn];
+ Float_t* phiT = new Float_t[nIn];
Int_t* injet = new Int_t[nIn];
//total energy in array
etaT[i] = lv->Eta();
phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
if (fReader->GetCutFlag(i) != 1) continue;
- fLego->Fill(etaT[i], phiT[i], ptT[i]);
+ fLego ->Fill(etaT[i], phiT[i], ptT[i]);
hPtTotal->Fill(ptT[i]);
etbgTotal+= ptT[i];
}
+
fJets->SetNinput(nIn);
// calculate total energy and fluctuation in map
Double_t meanpt = hPtTotal->GetMean();
- Double_t ptRMS = hPtTotal->GetRMS();
- Double_t npart = hPtTotal->GetEntries();
+ Double_t ptRMS = hPtTotal->GetRMS();
+ Double_t npart = hPtTotal->GetEntries();
Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
// arrays to hold jets
// add jets to list
Int_t* idxjets = new Int_t[nj];
Int_t nselectj = 0;
+// printf("Found %d jets \n", nj);
+
for(Int_t kj=0; kj<nj; kj++){
if ((etaJet[kj] > (header->GetJetEtaMax())) ||
(etaJet[kj] < (header->GetJetEtaMin())) ||
pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
en = TMath::Sqrt(px * px + py * py + pz * pz);
fJets->AddJet(px, py, pz, en);
+ AliAODJet jet(px, py, pz, en);
+ //jet.Print("");
+
+ AddJet(jet);
+
idxjets[nselectj] = kj;
nselectj++;
}
//delete
- delete ptT;
- delete etaT;
- delete phiT;
- delete injet;
+ delete [] ptT;
+ delete [] etaT;
+ delete [] phiT;
+ delete [] injet;
delete hPtTotal;
- delete etaJet;
- delete phiJet;
- delete etJet;
- delete etsigJet;
- delete etallJet;
- delete ncellsJet;
- delete multJet;
- delete idxjets;
- delete percentage;
- delete ncells;
- delete mult;
+ delete [] etaJet;
+ delete [] phiJet;
+ delete [] etJet;
+ delete [] etsigJet;
+ delete [] etallJet;
+ delete [] ncellsJet;
+ delete [] multJet;
+ delete [] idxjets;
+ delete [] percentage;
+ delete [] ncells;
+ delete [] mult;
}
//dump lego
// check enough space! *to be done*
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
- Float_t etCell[60000]; //! Cell Energy
- Float_t etaCell[60000]; //! Cell eta
- Float_t phiCell[60000]; //! Cell phi
- Int_t flagCell[60000]; //! Cell flag
+ const Int_t nBinsMax = 70000;
+
+ Float_t etCell[nBinsMax]; //! Cell Energy
+ Float_t etaCell[nBinsMax]; //! Cell eta
+ Float_t phiCell[nBinsMax]; //! Cell phi
+ Int_t flagCell[nBinsMax]; //! Cell flag
Int_t nCell = 0;
TAxis* xaxis = fLego->GetXaxis();
etCell[nCell] = e;
etaCell[nCell] = eta;
phiCell[nCell] = phi;
- flagCell[nCell] = 0; //default
+ flagCell[nCell] = 0; //default
nCell++;
}
}
Int_t * index = new Int_t[nCell];
TMath::Sort(nCell, etCell, index);
// variable used in centroide loop
- Float_t eta = 0.0;
- Float_t phi = 0.0;
- Float_t eta0 = 0.0;
- Float_t phi0 = 0.0;
- Float_t etab = 0.0;
- Float_t phib = 0.0;
- Float_t etas = 0.0;
- Float_t phis = 0.0;
- Float_t ets = 0.0;
- Float_t deta = 0.0;
- Float_t dphi = 0.0;
- Float_t dr = 0.0;
- Float_t etsb = 0.0;
+ Float_t eta = 0.0;
+ Float_t phi = 0.0;
+ Float_t eta0 = 0.0;
+ Float_t phi0 = 0.0;
+ Float_t etab = 0.0;
+ Float_t phib = 0.0;
+ Float_t etas = 0.0;
+ Float_t phis = 0.0;
+ Float_t ets = 0.0;
+ Float_t deta = 0.0;
+ Float_t dphi = 0.0;
+ Float_t dr = 0.0;
+ Float_t etsb = 0.0;
Float_t etasb = 0.0;
Float_t phisb = 0.0;
-
+ Float_t dphib = 0.0;
+
for(Int_t icell = 0; icell < nCell; icell++){
Int_t jcell = index[icell];
phisb = 0.0;
for(Int_t kcell =0; kcell < nCell; kcell++){
Int_t lcell = index[kcell];
- if(lcell == jcell) continue; // cell itself
- if(flagCell[lcell] != 0) continue; // cell used before
- if(etCell[lcell] > etCell[jcell]) continue;
+ if(lcell == jcell) continue; // cell itself
+ if(flagCell[lcell] != 0) continue; // cell used before
+ if(etCell[lcell] > etCell[jcell]) continue; // can this happen
//calculate dr
deta = etaCell[lcell] - eta;
- dphi = phiCell[lcell] - phi;
- if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
- if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
- dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ dphi = TMath::Abs(phiCell[lcell] - phi);
+ if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
if(dr <= rc){
// calculate offset from initiate cell
deta = etaCell[lcell] - eta0;
dphi = phiCell[lcell] - phi0;
- if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
- if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ if (dphi < - TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
+
etas = etas + etCell[lcell]*deta;
phis = phis + etCell[lcell]*dphi;
ets = ets + etCell[lcell];
eta = eta0 + etas/ets;
phi = phi0 + phis/ets;
// if cone does not move much, just go to next step
- dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib));
+ dphib = TMath::Abs(phi - phib);
+ if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
+ dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
if(dr <= minmove) break;
// cone should not move more than max_mov
dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
if(dr > maxmove){
- eta = etab;
- phi = phib;
- ets = etsb;
- etas = etasb;
- phis = phisb;
- }else{ // store this loop information
- etab=eta;
- phib=phi;
- etsb = ets;
- etasb = etas;
- phisb = phis;
+ eta = etab;
+ phi = phib;
+ ets = etsb;
+ etas = etasb;
+ phis = phisb;
+ } else { // store this loop information
+ etab = eta;
+ phib = phi;
+ etsb = ets;
+ etasb = etas;
+ phisb = phis;
}
- }
+ } // inside cone
}//end of cells loop looking centroide
//avoid cones overloap (to be implemented in the future)
//delete
- delete index;
- delete idx;
+ delete[] index;
+ delete[] idx;
}
////////////////////////////////////////////////////////////////////////
{
fLego->Reset();
fJets->ClearJets();
+ AliJetFinder::Reset();
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV1::WriteJHeaderToFile()
{
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
- fOut->cd();
header->Write();
}
header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
header->GetLegoPhiMin(), header->GetLegoPhiMax());
+
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff