* about the suitability of this software for any purpose. It is *
* 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 "AliJetHeader.h"
+
+
+#include "AliAODJet.h"
+#include "AliLog.h"
ClassImp(AliUA1JetFinderV1)
-////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////
-AliUA1JetFinderV1::AliUA1JetFinderV1():
- fHeader(0x0),
- fLego(0x0)
+AliUA1JetFinderV1::AliUA1JetFinderV1() :
+ AliJetFinder(),
+ fLego(0)
{
- // Default constructor
-
+ // Constructor
}
////////////////////////////////////////////////////////////////////////
{
// destructor
+ delete fLego;
+ fLego = 0;
}
////////////////////////////////////////////////////////////////////////
//5) fill AliJet list
// transform input to pt,eta,phi plus lego
+
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
TClonesArray *lvArray = fReader->GetMomentumArray();
Int_t nIn = lvArray->GetEntries();
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
Int_t* multJet = new Int_t[30];
Int_t nJets; // to hold number of jets found by algorithm
Int_t nj; // number of jets accepted
- Float_t prec = fHeader->GetPrecBg();
+ Float_t prec = header->GetPrecBg();
Float_t bgprec = 1;
while(bgprec > prec){
//reset jet arrays in memory
nJets = 0;
nj = 0;
// reset particles-jet array in memory
- memset(injet,0,sizeof(Int_t)*nIn);
+ memset(injet,-1,sizeof(Int_t)*nIn);
//run cone algorithm finder
RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
//run background subtraction
- if(nJets > fHeader->GetNAcceptJets()) // limited number of accepted jets per event
- nj = fHeader->GetNAcceptJets();
+ if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
+ nj = header->GetNAcceptJets();
else
nj = nJets;
//subtract background
Float_t etbgTotalN = 0.0; //new background
- if(fHeader->GetBackgMode() == 1) // standar
+ if(header->GetBackgMode() == 1) // standar
SubtractBackg(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(fHeader->GetBackgMode() == 2) //cone
+ if(header->GetBackgMode() == 2) //cone
SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(fHeader->GetBackgMode() == 3) //ratio
+ if(header->GetBackgMode() == 3) //ratio
SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(fHeader->GetBackgMode() == 4) //statistic
+ if(header->GetBackgMode() == 4) //statistic
SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
//calc precision
if(etbgTotalN != 0.0)
etbgTotal = etbgTotalN; // update with new background estimation
} //end while
+ // add tracks to the jet if it wasn't yet done
+ if (header->GetBackgMode() == 0){
+ Float_t rc= header->GetRadius();
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ for(Int_t ijet=0; ijet<nj; ijet++){
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if(dr <= rc){ // particles inside this cone
+ injet[jpart] = ijet;
+ break;
+ }
+ }// end jets loop
+ } //end particle loop
+ }
+
// add jets to list
Int_t* idxjets = new Int_t[nj];
Int_t nselectj = 0;
+// printf("Found %d jets \n", nj);
+
+ TRefArray *refs = 0;
+ Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+ if (fromAod) refs = fReader->GetReferences();
for(Int_t kj=0; kj<nj; kj++){
- if(etJet[kj] > fHeader->GetMinJetEt()){ // check if et jets > et min
+ if ((etaJet[kj] > (header->GetJetEtaMax())) ||
+ (etaJet[kj] < (header->GetJetEtaMin())) ||
+ (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
Float_t px, py,pz,en; // convert to 4-vector
px = etJet[kj] * TMath::Cos(phiJet[kj]);
py = etJet[kj] * TMath::Sin(phiJet[kj]);
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);
+
+ if (fromAod){
+ for(Int_t jpart = 0; jpart < nIn; jpart++) // loop for all particles in array
+ if (injet[jpart] == kj && fReader->GetCutFlag(jpart) == 1)
+ jet.AddTrack(refs->At(jpart)); // check if the particle belongs to the jet and add the ref
+ }
+
+ //jet.Print("");
+
+ AddJet(jet);
+
idxjets[nselectj] = kj;
nselectj++;
- }
- }
+ } //end particle loop
+
//add signal percentage and total signal in AliJets for analysis tool
Float_t* percentage = new Float_t[nselectj];
Int_t* ncells = new Int_t[nselectj];
Int_t* mult = new Int_t[nselectj];
-
for(Int_t i = 0; i< nselectj; i++){
percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
ncells[i] = ncellsJet[idxjets[i]];
mult[i] = multJet[idxjets[i]];
}
- fJets->SetNCells(ncells);
- fJets->SetPtFromSignal(percentage);
- fJets->SetMultiplicities(mult);
- fJets->SetInJet(injet);
- fJets->SetEtaIn(etaT);
- fJets->SetPhiIn(phiT);
- fJets->SetPtIn(ptT);
- fJets->SetEtAvg(etbgTotal/(4*(fHeader->GetLegoEtaMax())*TMath::Pi()));
-
+ //add particle-injet relationship ///
+ for(Int_t bj = 0; bj < nIn; bj++){
+ if(injet[bj] == -1) continue; //background particle
+ Int_t bflag = 0;
+ for(Int_t ci = 0; ci< nselectj; ci++){
+ if(injet[bj] == idxjets[ci]){
+ injet[bj]= ci;
+ bflag++;
+ break;
+ }
+ }
+ if(bflag == 0) injet[bj] = -1; // set as background particle
+ }
//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*
- Float_t etCell[60000]; //! Cell Energy
- Float_t etaCell[60000]; //! Cell eta
- Float_t phiCell[60000]; //! Cell phi
- Int_t flagCell[60000]; //! Cell flag
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
+
+ const Int_t nBinEta = header->GetLegoNbinEta();
+ const Int_t nBinPhi = header->GetLegoNbinPhi();
+ if((nBinPhi*nBinEta)>nBinsMax){
+ AliError("Too many bins of the ETA-PHI histogram");
+ }
+
+ Float_t etCell[nBinsMax]; //! Cell Energy
+ Float_t etaCell[nBinsMax]; //! Cell eta
+ Float_t phiCell[nBinsMax]; //! Cell phi
+ Short_t flagCell[nBinsMax]; //! Cell flag
Int_t nCell = 0;
TAxis* xaxis = fLego->GetXaxis();
TAxis* yaxis = fLego->GetYaxis();
Float_t e = 0.0;
- for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
- for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
+ for (Int_t i = 1; i <= nBinEta; i++) {
+ for (Int_t j = 1; j <= nBinPhi; j++) {
e = fLego->GetBinContent(i,j);
if (e < 0.0) continue; // don't include this cells
Float_t eta = xaxis->GetBinCenter(i);
etCell[nCell] = e;
etaCell[nCell] = eta;
phiCell[nCell] = phi;
- flagCell[nCell] = 0; //default
+ flagCell[nCell] = 0; //default
nCell++;
}
}
// Parameters from header
- Float_t minmove = fHeader->GetMinMove();
- Float_t maxmove = fHeader->GetMaxMove();
- Float_t rc = fHeader->GetRadius();
- Float_t etseed = fHeader->GetEtSeed();
- //Float_t etmin = fHeader->GetMinJetEt();
+ Float_t minmove = header->GetMinMove();
+ Float_t maxmove = header->GetMaxMove();
+ Float_t rc = header->GetRadius();
+ Float_t etseed = header->GetEtSeed();
+ //Float_t etmin = header->GetMinJetEt();
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)
//flag cells in Rc, estimate total energy in cone
Float_t etCone = 0.0;
Int_t nCellIn = 0;
- rc = fHeader->GetRadius();
+ rc = header->GetRadius();
for(Int_t ncell =0; ncell < nCell; ncell++){
if(flagCell[ncell] != 0) continue; // cell used before
//calculate dr
}
}
//store tmp jet info !!!
- // reject tmp outside acceptable eta range
- if ((eta < (fHeader->GetJetEtaMax())) &&
- (eta > (fHeader->GetJetEtaMin())) &&
- (etbmax < etcmin) ){
+ if(etbmax < etcmin) {
etaAlgoJet[nJets] = eta;
phiAlgoJet[nJets] = phi;
etAlgoJet[nJets] = etCone;
//delete
- delete index;
- delete idx;
+ delete[] index;
+ delete[] idx;
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV1::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
- Float_t* ptT, Float_t* etaT, Float_t* phiT,
- Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV1::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+ Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
+ Int_t* multJet, Int_t* injet)
{
//background subtraction using cone method but without correction in dE/deta distribution
//calculate energy inside and outside cones
- Float_t rc= fHeader->GetRadius();
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ Float_t rc= header->GetRadius();
Float_t etIn[30];
Float_t etOut = 0;
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
+ // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
for(Int_t ijet=0; ijet<nJ; ijet++){
Float_t deta = etaT[jpart] - etaJet[ijet];
Float_t dphi = phiT[jpart] - phiJet[ijet];
if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
if(dr <= rc){ // particles inside this cone
- etIn[ijet] += ptT[jpart];
- if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
multJet[ijet]++;
injet[jpart] = ijet;
+ if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
+ etIn[ijet] += ptT[jpart];
+ if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
+ }
break;
}
}// end jets loop
- if(injet[jpart] == 0) etOut += ptT[jpart]; // particle outside cones
+ if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
+ etOut += ptT[jpart]; // particle outside cones and pt cut
} //end particle loop
//estimate jets and background areas
Float_t areaJet[30];
- Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+ Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
for(Int_t k=0; k<nJ; k++){
Float_t detamax = etaJet[k] + rc;
Float_t detamin = etaJet[k] - rc;
Float_t accmax = 0.0; Float_t accmin = 0.0;
- if(detamax > fHeader->GetLegoEtaMax()){ // sector outside etamax
- Float_t h = fHeader->GetLegoEtaMax() - etaJet[k];
+ if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
+ Float_t h = header->GetLegoEtaMax() - etaJet[k];
accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
}
- if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin
- Float_t h = fHeader->GetLegoEtaMax() + etaJet[k];
+ if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
+ Float_t h = header->GetLegoEtaMax() + etaJet[k];
accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
}
areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
}
// estimate new total background
- Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+ Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV1::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN,
- Float_t* ptT, Float_t* etaT, Float_t* phiT,
- Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV1::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+ Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
+ Int_t* multJet, Int_t* injet)
{
//background subtraction using statistical method
-
- Float_t etbgStat = fHeader->GetBackgStat(); // pre-calculated background
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
//calculate energy inside
- Float_t rc= fHeader->GetRadius();
+ Float_t rc= header->GetRadius();
Float_t etIn[30];
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
+ //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
for(Int_t ijet=0; ijet<nJ; ijet++){
Float_t deta = etaT[jpart] - etaJet[ijet];
Float_t dphi = phiT[jpart] - phiJet[ijet];
if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
if(dr <= rc){ // particles inside this cone
- etIn[ijet]+= ptT[jpart];
- if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
multJet[ijet]++;
injet[jpart] = ijet;
+ if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
+ etIn[ijet]+= ptT[jpart];
+ if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+ }
break;
}
}// end jets loop
//calc jets areas
Float_t areaJet[30];
- Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+ Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
for(Int_t k=0; k<nJ; k++){
Float_t detamax = etaJet[k] + rc;
Float_t detamin = etaJet[k] - rc;
Float_t accmax = 0.0; Float_t accmin = 0.0;
- if(detamax > fHeader->GetLegoEtaMax()){ // sector outside etamax
- Float_t h = fHeader->GetLegoEtaMax() - etaJet[k];
+ if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
+ Float_t h = header->GetLegoEtaMax() - etaJet[k];
accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
}
- if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin
- Float_t h = fHeader->GetLegoEtaMax() + etaJet[k];
+ if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
+ Float_t h = header->GetLegoEtaMax() + etaJet[k];
accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
}
areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV1::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
- Float_t* ptT, Float_t* etaT, Float_t* phiT,
- Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV1::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+ Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
+ Int_t* multJet, Int_t* injet)
{
// Cone background subtraction method taking into acount dEt/deta distribution
-
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
//general
- Float_t rc= fHeader->GetRadius();
- Float_t etamax = fHeader->GetLegoEtaMax();
- Float_t etamin = fHeader->GetLegoEtaMin();
+ Float_t rc= header->GetRadius();
+ Float_t etamax = header->GetLegoEtaMax();
+ Float_t etamin = header->GetLegoEtaMin();
Int_t ndiv = 100;
// jet energy and area arrays
//fill energies
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
Float_t deta = etaT[jpart] - etaJet[ijet];
- Float_t dphi = phiT[jpart] - phiJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
- if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
- Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
if(dr <= rc){ // particles inside this cone
- hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
- multJet[ijet]++;
- if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
injet[jpart] = ijet;
+ multJet[ijet]++;
+ if((fReader->GetCutFlag(jpart)) == 1){// pt cut
+ hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
+ if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+ }
break;
}
}// end jets loop
- if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
+ hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
//calc areas
// calc background total
Double_t etOut = hEtBackg->Integral();
Double_t areaOut = hAreaBackg->Integral();
- Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+ Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
//delete
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV1::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
- Float_t* ptT, Float_t* etaT, Float_t* phiT,
- Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV1::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ, Float_t& etbgTotalN,
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+ Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
+ Int_t* multJet, Int_t* injet)
{
// Ratio background subtraction method taking into acount dEt/deta distribution
-
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
//factor F calc before
- Float_t bgRatioCut = fHeader->GetBackgCutRatio();
+ Float_t bgRatioCut = header->GetBackgCutRatio();
//general
- Float_t rc= fHeader->GetRadius();
- Float_t etamax = fHeader->GetLegoEtaMax();
- Float_t etamin = fHeader->GetLegoEtaMin();
+ Float_t rc= header->GetRadius();
+ Float_t etamax = header->GetLegoEtaMax();
+ Float_t etamin = header->GetLegoEtaMin();
Int_t ndiv = 100;
// jet energy and area arrays
//fill energies
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- if((fReader->GetCutFlag(jpart)) != 1) continue;
+ //if((fReader->GetCutFlag(jpart)) != 1) continue;
for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
Float_t deta = etaT[jpart] - etaJet[ijet];
Float_t dphi = phiT[jpart] - phiJet[ijet];
if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
if(dr <= rc){ // particles inside this cone
- hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
- if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
multJet[ijet]++;
injet[jpart] = ijet;
+ if((fReader->GetCutFlag(jpart)) == 1){ //pt cut
+ hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
+ if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+ }
break;
}
}// end jets loop
- if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
//calc areas
// calc background total
Double_t etOut = hEtBackg->Integral();
Double_t areaOut = hAreaBackg->Integral();
- Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+ Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
//delete
void AliUA1JetFinderV1::Reset()
{
fLego->Reset();
- fJets->ClearJets();
+ AliJetFinder::Reset();
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV1::WriteJHeaderToFile()
{
- fOut->cd();
- fHeader->Write();
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ header->Write();
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV1::Init()
{
// initializes some variables
-
- // book lego
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
fLego = new
TH2F("legoH","eta-phi",
- fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
- fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
- fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
+ header->GetLegoNbinEta(), header->GetLegoEtaMin(),
+ header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
+ header->GetLegoPhiMin(), header->GetLegoPhiMax());
+ // Do not store in current dir
+ fLego->SetDirectory(0);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff