#include "AliUA1JetFinderV2.h"
#include "AliUA1JetHeaderV1.h"
#include "AliJetUnitArray.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetReader.h"
+#include "AliJetHeader.h"
class TArrayF;
class TFile;
AliUA1JetFinderV2::AliUA1JetFinderV2() :
AliJetFinder(),
fLego(0),
- fDebug(0),
fOpt(0)
{
//
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
TClonesArray* lvArray = fReader->GetMomentumArray();
Int_t nIn = lvArray->GetEntries();
+ fDebug = fHeader->GetDebug();
if (nIn == 0) return;
Float_t* ptT = new Float_t[nIn];
Float_t* etaT = new Float_t[nIn];
Float_t* phiT = new Float_t[nIn];
- Float_t* cFlagT = new Float_t[nIn]; // Temporarily added
- Float_t* sFlagT = new Float_t[nIn]; // Temporarily added
+ Int_t* cFlagT = new Int_t[nIn]; // Temporarily added
+ Int_t* sFlagT = new Int_t[nIn]; // Temporarily added
Int_t* injet = new Int_t[nIn];
-
+
+ for (Int_t i = 0; i < nIn; i++) {
+ ptT[i] = 0.;
+ etaT[i] = 0.;
+ phiT[i] = 0.;
+ cFlagT[i] = 0;
+ sFlagT[i] = 0;
+ injet[i] = 0;
+ }
//total energy in array
Float_t etbgTotal = 0.0;
TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
// arrays to hold jets
- Float_t* etaJet = new Float_t[30]; // eta jet
- Float_t* phiJet = new Float_t[30]; // phi jet
- Float_t* etJet = new Float_t[30]; // et jet
- Float_t* etsigJet = new Float_t[30]; // signal et in jet
- Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
- Int_t* ncellsJet = new Int_t[30];
- Int_t* multJet = new Int_t[30];
+ Float_t etaJet[30]; // eta jet
+ Float_t phiJet[30]; // phi jet
+ Float_t etJet[30]; // et jet
+ Float_t etsigJet[30]; // signal et in jet
+ Float_t etallJet[30]; // total et in jet (tmp variable)
+ Int_t ncellsJet[30];
+ Int_t multJet[30];
//--- Added for jet reordering at the end of the jet finding procedure
- Float_t* etaJetOk = new Float_t[30];
- Float_t* phiJetOk = new Float_t[30];
- Float_t* etJetOk = new Float_t[30];
- Float_t* etsigJetOk = new Float_t[30]; // signal et in jet
- Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
- Int_t* ncellsJetOk = new Int_t[30];
- Int_t* multJetOk = new Int_t[30];
+ Float_t etaJetOk[30];
+ Float_t phiJetOk[30];
+ Float_t etJetOk[30];
+ Float_t etsigJetOk[30]; // signal et in jet
+ Float_t etallJetOk[30]; // total et in jet (tmp variable)
+ Int_t ncellsJetOk[30];
+ Int_t multJetOk[30];
//--------------------------
Int_t nJets; // to hold number of jets found by algorithm
Int_t nj; // number of jets accepted
if(header->GetBackgMode() == 4) //statistic
SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
//calc precision
- if(etbgTotalN != 0.0)
+ if(TMath::Abs(etbgTotalN) > 0.001)
bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
else
bgprec = 0;
delete[] cFlagT;
delete[] sFlagT;
delete[] injet;
- delete[] hPtTotal;
- delete[] etaJet;
- delete[] phiJet;
- delete[] etJet;
- delete[] etsigJet;
- delete[] etallJet;
- delete[] ncellsJet;
- delete[] multJet;
+ delete hPtTotal;
delete[] idxjets;
+ delete[] idx;
+
delete[] percentage;
delete[] ncells;
delete[] mult;
- //--- Added for jet reordering
- delete etaJetOk;
- delete phiJetOk;
- delete etJetOk;
- delete etsigJetOk;
- delete etallJetOk;
- delete ncellsJetOk;
- delete multJetOk;
//--------------------------
}
Int_t nIn = fUnit->GetEntries();
Float_t ptMin = fReader->GetReaderHeader()->GetPtCut();
- fDebug = fReader->GetReaderHeader()->GetDebug();
if (nIn == 0) return;
Int_t nCandidateCut = 0;
Int_t* detT = new Int_t[nCandidate];
Float_t* etaT = new Float_t[nCandidate];
Float_t* phiT = new Float_t[nCandidate];
- Float_t* cFlagT = new Float_t[nCandidate];
- Float_t* cFlag2T = new Float_t[nCandidate];
- Float_t* sFlagT = new Float_t[nCandidate];
+ Int_t* cFlagT = new Int_t[nCandidate];
+ Int_t* cFlag2T = new Int_t[nCandidate];
+ Int_t* sFlagT = new Int_t[nCandidate];
Float_t* cClusterT = new Float_t[nCandidate];
Int_t* vectT = new Int_t[nCandidate];
Int_t loop1 = 0;
Int_t* injet = new Int_t[nCandidate];
Int_t* sflag = new Int_t[nCandidate];
+
+ for(Int_t i = 0; i < nCandidate; i++) {
+ ptT[i] = 0;
+ en2T[i] = 0;
+ pt2T[i] = 0;
+ detT[i] = 0;
+ etaT[i] = 0;
+ phiT[i] = 0;
+ cFlagT[i] = 0;
+ cFlag2T[i] = 0;
+ sFlagT[i] = 0;
+ cClusterT[i] = 0;
+ vectT[i] = 0;
+ injet[i] = 0;
+ sflag[i] = 0;
+}
+
TRefArray* trackRef = new TRefArray();
//total energy in array
- Float_t etbgTotal = 0.0;
+ Float_t etbgTotal = 0.0;
TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
// Input cell info
Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
Int_t *flagCell2 = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
-
+ for(Int_t i = 0; i < nIn; i++) {
+ etCell[i] = 0.;
+ etaCell[i] = 0.;
+ phiCell[i] = 0.;
+ flagCell[i] = 0;
+ etCell2[i] = 0.;
+ etaCell2[i] = 0.;
+ phiCell2[i] = 0.;
+ flagCell2[i] = 0;
+ }
// Information extracted from fUnitArray
// Load input vectors and calculate total energy in array
for(Int_t i=0; i<nIn; i++)
Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
// arrays to hold jets
- Float_t* etaJet = new Float_t[30];
- Float_t* phiJet = new Float_t[30];
- Float_t* etJet = new Float_t[30];
- Float_t* etsigJet = new Float_t[30]; //signal et in jet
- Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
- Int_t* ncellsJet = new Int_t[30];
- Int_t* multJet = new Int_t[30];
+ Float_t etaJet[30];
+ Float_t phiJet[30];
+ Float_t etJet[30];
+ Float_t etsigJet[30]; //signal et in jet
+ Float_t etallJet[30]; // total et in jet (tmp variable)
+ Int_t ncellsJet[30];
+ Int_t multJet[30];
//--- Added by me for jet reordering at the end of the jet finding procedure
- Float_t* etaJetOk = new Float_t[30];
- Float_t* phiJetOk = new Float_t[30];
- Float_t* etJetOk = new Float_t[30];
- Float_t* etsigJetOk = new Float_t[30]; //signal et in jet
- Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
- Int_t* ncellsJetOk = new Int_t[30];
- Int_t* multJetOk = new Int_t[30];
+ Float_t etaJetOk[30];
+ Float_t phiJetOk[30];
+ Float_t etJetOk[30];
+ Float_t etsigJetOk[30]; //signal et in jet
+ Float_t etallJetOk[30]; // total et in jet (tmp variable)
+ Int_t ncellsJetOk[30];
+ Int_t multJetOk[30];
//--------------------------
Int_t nJets; // to hold number of jets found by algorithm
Int_t nj; // number of jets accepted
multJetOk[p] = multJet[idx[p]];
}
+ TRefArray *refs = 0;
+ Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+ if (fromAod) refs = fReader->GetReferences();
+ Int_t nTracks = 0;
+ if (fromAod) nTracks = ((TRefArray*)refs)->GetEntries();
+ Int_t* trackinjet = new Int_t[nTracks];
+ for(Int_t it=0; it<nTracks; it++) trackinjet[it]=-1;
+
for(Int_t kj=0; kj<nj; kj++)
{
if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
AliAODJet jet(px, py, pz, en);
jet.Print("");
+ if (fromAod){
+ for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
+ Float_t deta = ((AliAODTrack*)refs->At(jpart))->Eta() - etaJetOk[kj];
+ Float_t dphi = ((AliAODTrack*)refs->At(jpart))->Phi() - phiJetOk[kj];
+ 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 <= header->GetRadius() && fReader->GetCutFlag(jpart) == 1) {
+ // particles inside this cone
+ if(trackinjet[jpart]==-1) {
+ trackinjet[jpart] = kj;
+ } else if(fDebug>10) {
+ printf("The track already belongs to jet %d \n",trackinjet[jpart]);
+ }
+ }
+ if(trackinjet[jpart]==kj)
+ jet.AddTrack(refs->At(jpart)); // check if the particle belongs to the jet and add the ref
+ }
+ }
+
AddJet(jet);
idxjets[nselectj] = kj;
//delete
- delete ptT;
- delete en2T;
- delete pt2T;
- delete etaT;
- delete phiT;
+ delete [] ptT;
+ delete [] en2T;
+ delete [] pt2T;
+ delete [] etaT;
+ delete [] phiT;
+ delete [] detT;
+ delete [] cFlagT;
+ delete [] cFlag2T;
+ delete [] sFlagT;
+ delete [] cClusterT;
+ delete [] vectT;
+ delete [] injet;
+ delete [] sflag;
trackRef->Delete();
delete trackRef;
- delete detT;
- delete cFlagT;
- delete cFlag2T;
- delete sFlagT;
- delete cClusterT;
- delete vectT;
- delete injet;
- delete sflag;
+
delete hPtTotal;
- delete etCell;
- delete etaCell;
- delete phiCell;
- delete flagCell;
- delete etCell2;
- delete etaCell2;
- delete phiCell2;
- delete flagCell2;
- delete etaJet;
- delete phiJet;
- delete etJet;
- delete etsigJet;
- delete etallJet;
- delete ncellsJet;
- delete multJet;
- //--- Added for jet reordering
- delete etaJetOk;
- delete phiJetOk;
- delete etJetOk;
- delete etsigJetOk;
- delete etallJetOk;
- delete ncellsJetOk;
- delete multJetOk;
+ delete [] etCell;
+ delete [] etaCell;
+ delete [] phiCell;
+ delete [] flagCell;
+ delete [] etCell2;
+ delete [] etaCell2;
+ delete [] phiCell2;
+ delete [] flagCell2;
//--------------------------
- delete idxjets;
- delete percentage;
- delete ncells;
- delete mult;
+ delete [] idxjets;
+ delete [] idx;
+ delete [] trackinjet;
+
+ delete [] percentage;
+ delete [] ncells;
+ delete [] mult;
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
- Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2,
- const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal,
- Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet,
- Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
+void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* const etaCell, Float_t* phiCell,
+ Int_t* const flagCell, const Float_t* etCell2, const Float_t* etaCell2, const Float_t* phiCell2,
+ const Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal,
+ Int_t& nJets, Float_t* const etJet, Float_t* const etaJet, Float_t* const phiJet,
+ Float_t* const etallJet, Int_t* const ncellsJet)
{
//
// Main method for jet finding
// UA1 base cone finder
//
- Int_t nCell = nIn;
- fDebug = fReader->GetReaderHeader()->GetDebug();
+ Int_t nCell = nIn;
// Dump lego
// Check enough space! *to be done*
Float_t etseed = header->GetEtSeed();
// Tmp array of jets form algoritm
- Float_t etaAlgoJet[30];
- Float_t phiAlgoJet[30];
- Float_t etAlgoJet[30];
- Int_t ncellsAlgoJet[30];
+ Float_t etaAlgoJet[30] = {0.};
+ Float_t phiAlgoJet[30] = {0.};
+ Float_t etAlgoJet[30] = {0.};
+ Int_t ncellsAlgoJet[30] = {0};
// Run algorithm//
}
//delete
- delete index;
+ delete[] index;
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
- Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
- Float_t* etallJet, Int_t* ncellsJet)
+ Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
+ Float_t* const etallJet, Int_t* const ncellsJet)
{
// Dump lego
// Check enough space! *to be done*
Float_t etseed = header->GetEtSeed();
// Tmp array of jets form algoritm
- Float_t etaAlgoJet[30];
- Float_t phiAlgoJet[30];
- Float_t etAlgoJet[30];
- Int_t ncellsAlgoJet[30];
+ Float_t etaAlgoJet[30] = {0.};
+ Float_t phiAlgoJet[30] = {0.};
+ Float_t etAlgoJet[30] = {0.};
+ Int_t ncellsAlgoJet[30] = {0};
// Run algorithm//
}
//delete
- delete index;
- delete idx;
+ delete [] index;
+ delete [] idx;
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT,
- const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
- const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet,
- const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t& etbgTotalN, const Float_t* ptT, const Int_t* vectT,
+ const Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT, const Int_t* cFlag2T,
+ const Int_t* sFlagT, Float_t* const etJet, const Float_t* etaJet, const Float_t* phiJet,
+ Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet)
{
//
// Background subtraction using cone method but without correction in dE/deta distribution
//calculate energy inside and outside cones
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
- Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
- fDebug = fReader->GetReaderHeader()->GetDebug();
+ fOpt = fReader->GetReaderHeader()->GetDetector();
Float_t rc= header->GetRadius();
- Float_t etIn[30];
+ Float_t etIn[30] = {0.};
Float_t etOut = 0;
-
- for(Int_t j=0;j<30;j++){etIn[j]=0.;}
-
+
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
for(Int_t ijet=0; ijet<nJ; ijet++){
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV2::SubtractBackgC(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)
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+ Float_t* const etJet, const Float_t* etaJet, const Float_t* phiJet,
+ Float_t* const etsigJet,Int_t* const multJet, Int_t* const injet)
{
//background subtraction using cone method but without correction in dE/deta distribution
//calculate energy inside and outside cones
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
Float_t rc= header->GetRadius();
- Float_t etIn[30];
+ Float_t etIn[30] = {0.};
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
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
- const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
- const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
- Float_t* etsigJet, Int_t* multJet, Int_t* injet)
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT,
+ const Int_t* sFlagT, Float_t* const etJet, const Float_t* etaJet, const Float_t* phiJet,
+ Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet)
{
//background subtraction using statistical method
//calculate energy inside
Float_t rc= header->GetRadius();
- Float_t etIn[30];
+ Float_t etIn[30] = {0.};
for(Int_t jpart = 0; jpart < nIn; jpart++)
{ // loop for all particles in array
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
- Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
- Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT, const Int_t* sFlagT,
+ Float_t* const etJet, const Float_t* etaJet, const Float_t* phiJet,
+ Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet)
{
// Cone background subtraction method taking into acount dEt/deta distribution
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
TH1F* hAreaJet[30];
for(Int_t mjet=0; mjet<nJ; mjet++){
char hEtname[256]; char hAreaname[256];
- sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+ snprintf(hEtname, 256, "hEtJet%d", mjet); snprintf(hAreaname, 256, "hAreaJet%d", mjet);
hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
- Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
- Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
- Int_t* multJet, Int_t* injet)
+ const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT, const Int_t* sFlagT,
+ Float_t* const etJet, const Float_t* etaJet, const Float_t* phiJet,
+ Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet)
{
// Ratio background subtraction method taking into acount dEt/deta distribution
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
TH1F* hAreaJet[30];
for(Int_t mjet=0; mjet<nJ; mjet++){
char hEtname[256]; char hAreaname[256];
- sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+ snprintf(hEtname, 256, "hEtJet%d", mjet); snprintf(hAreaname, 256, "hAreaJet%d", mjet);
hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax); // change range
hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax); // change range
}
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::WriteJHeaderToFile()
+void AliUA1JetFinderV2::WriteJHeaderToFile() const
{
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
header->Write();
header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
header->GetLegoPhiMin(), header->GetLegoPhiMax());
- fDebug = fReader->GetReaderHeader()->GetDebug();
+ fDebug = fHeader->GetDebug();
fOpt = fReader->GetReaderHeader()->GetDetector();
// Tasks initialization