#include "AliUA1JetHeaderV1.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
-#include "AliJet.h"
+#include "AliJetHeader.h"
+
+
#include "AliAODJet.h"
#include "AliLog.h"
AliUA1JetFinderV1::AliUA1JetFinderV1() :
AliJetFinder(),
- fLego(0)
+ fLego(0),
+ fhEtBackg(0),
+ fhAreaBackg(0)
{
// Constructor
+ for(int i = 0;i < kMaxJets;i++){
+ fhAreaJet[i] = fhEtJet[i] = 0;
+ }
}
////////////////////////////////////////////////////////////////////////
// destructor
delete fLego;
fLego = 0;
+ if(fhEtBackg)delete fhEtBackg;
+ fhEtBackg = 0;
+ if( fhAreaBackg) delete fhAreaBackg;
+ fhAreaBackg = 0;
+ for(int i = 0;i < kMaxJets;i++){
+ if(fhAreaJet[i])delete fhAreaJet[i];
+ if(fhEtJet[i]) delete fhEtJet[i];
+ fhAreaJet[i] = fhEtJet[i] = 0;
+ }
+
}
////////////////////////////////////////////////////////////////////////
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
TClonesArray *lvArray = fReader->GetMomentumArray();
Int_t nIn = lvArray->GetEntries();
- if (nIn == 0) return;
+ if (nIn <= 0) return;
// local arrays for input
+ // ToDo: check memory fragmentation, maybe better to
+ // define them globally and resize as needed
+ // Fragementation should be worse for low mult...
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];
+ memset(ptT,0,sizeof(Float_t)*nIn);
+ memset(etaT,0,sizeof(Float_t)*nIn);
+ memset(phiT,0,sizeof(Float_t)*nIn);
+
+
+ // load input vectors and calculate total energy in array
+
//total energy in array
Float_t etbgTotal = 0.0;
- TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
+ Float_t npart = 0;
+ Float_t etbg2 = 0;
- // load input vectors and calculate total energy in array
for (Int_t i = 0; i < nIn; i++){
TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
ptT[i] = lv->Pt();
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]);
- hPtTotal->Fill(ptT[i]);
+ npart += 1;
etbgTotal+= ptT[i];
+ etbg2 += ptT[i]*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 meanpt = 0;
+ Double_t ptRMS = 0;
+ if(npart>0){
+ meanpt = etbgTotal/npart;
+ etbg2 = etbg2/npart;
+ if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
+ ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
+ }
+ }
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[kMaxJets];
+ Float_t phiJet[kMaxJets];
+ Float_t etJet[kMaxJets];
+ Float_t etsigJet[kMaxJets]; //signal et in jet
+ Float_t etallJet[kMaxJets]; // total et in jet (tmp variable)
+ Int_t ncellsJet[kMaxJets];
+ Int_t multJet[kMaxJets];
Int_t nJets; // to hold number of jets found by algorithm
Int_t nj; // number of jets accepted
Float_t prec = header->GetPrecBg();
Float_t bgprec = 1;
while(bgprec > prec){
//reset jet arrays in memory
- memset(etaJet,0,sizeof(Float_t)*30);
- memset(phiJet,0,sizeof(Float_t)*30);
- memset(etJet,0,sizeof(Float_t)*30);
- memset(etallJet,0,sizeof(Float_t)*30);
- memset(etsigJet,0,sizeof(Float_t)*30);
- memset(ncellsJet,0,sizeof(Int_t)*30);
- memset(multJet,0,sizeof(Int_t)*30);
+ memset(etaJet,0,sizeof(Float_t)*kMaxJets);
+ memset(phiJet,0,sizeof(Float_t)*kMaxJets);
+ memset(etJet,0,sizeof(Float_t)*kMaxJets);
+ memset(etallJet,0,sizeof(Float_t)*kMaxJets);
+ memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
+ memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
+ memset(multJet,0,sizeof(Int_t)*kMaxJets);
nJets = 0;
nj = 0;
// reset particles-jet array in memory
}
// add jets to list
- Int_t* idxjets = new Int_t[nj];
+ Int_t idxjets[kMaxJets];
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();
+ Float_t rc= header->GetRadius();
for(Int_t kj=0; kj<nj; kj++){
if ((etaJet[kj] > (header->GetJetEtaMax())) ||
(etaJet[kj] < (header->GetJetEtaMin())) ||
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){
//jet.Print("");
+ // calculate the area of the jet
+ Float_t detamax = etaJet[kj] + rc;
+ Float_t detamin = etaJet[kj] - rc;
+ Float_t accmax = 0.0; Float_t accmin = 0.0;
+ if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
+ Float_t h = header->GetLegoEtaMax() - etaJet[kj];
+ accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+ if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
+ Float_t h = header->GetLegoEtaMax() + etaJet[kj];
+ accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+ Float_t areaJet = rc*rc*TMath::Pi() - accmax - accmin;
+ // set both areas
+ jet.SetEffArea(areaJet,areaJet);
+
AddJet(jet);
idxjets[nselectj] = kj;
} //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];
+ Float_t percentage[kMaxJets];
+ Int_t ncells[kMaxJets];
+ Int_t mult[kMaxJets];
for(Int_t i = 0; i< nselectj; i++){
percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
ncells[i] = ncellsJet[idxjets[i]];
}
if(bflag == 0) injet[bj] = -1; // set as background particle
}
- 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*(header->GetLegoEtaMax())*TMath::Pi()));
-
//delete
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;
-
-
}
////////////////////////////////////////////////////////////////////////
// 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[kMaxJets];
+ Float_t phiAlgoJet[kMaxJets];
+ Float_t etAlgoJet[kMaxJets];
+ Int_t ncellsAlgoJet[kMaxJets];
//run algorithm//
// sort cells by et
- Int_t * index = new Int_t[nCell];
+ Int_t index[nBinsMax];
TMath::Sort(nCell, etCell, index);
// variable used in centroide loop
Float_t eta = 0.0;
}
//store tmp jet info !!!
if(etbmax < etcmin) {
- etaAlgoJet[nJets] = eta;
- phiAlgoJet[nJets] = phi;
- etAlgoJet[nJets] = etCone;
- ncellsAlgoJet[nJets] = nCellIn;
- nJets++;
- }
-
+ if(nJets<kMaxJets){
+ etaAlgoJet[nJets] = eta;
+ phiAlgoJet[nJets] = phi;
+ etAlgoJet[nJets] = etCone;
+ ncellsAlgoJet[nJets] = nCellIn;
+ nJets++;
+ }
+ else{
+ AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
+ break;
+ }
+ }
} // end of cells loop
//reorder jets by et in cone
//sort jets by energy
- Int_t * idx = new Int_t[nJets];
- TMath::Sort(nJets, etAlgoJet, idx);
+ Int_t idx[kMaxJets];
+ TMath::Sort(nJets, etAlgoJet, idx); // sort only the found jets
for(Int_t p = 0; p < nJets; p++){
etaJet[p] = etaAlgoJet[idx[p]];
phiJet[p] = phiAlgoJet[idx[p]];
ncellsJet[p] = ncellsAlgoJet[idx[p]];
}
-
- //delete
- delete[] index;
- delete[] idx;
-
}
////////////////////////////////////////////////////////////////////////
//calculate energy inside and outside cones
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
Float_t rc= header->GetRadius();
- Float_t etIn[30];
+ Float_t etIn[kMaxJets];
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
} //end particle loop
//estimate jets and background areas
- Float_t areaJet[30];
+ Float_t areaJet[kMaxJets];
Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
for(Int_t k=0; k<nJ; k++){
Float_t detamax = etaJet[k] + rc;
//calculate energy inside
Float_t rc= header->GetRadius();
- Float_t etIn[30];
+ Float_t etIn[kMaxJets];
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
//if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
} //end particle loop
//calc jets areas
- Float_t areaJet[30];
+ Float_t areaJet[kMaxJets];
Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
for(Int_t k=0; k<nJ; k++){
Float_t detamax = etaJet[k] + rc;
Int_t ndiv = 100;
// jet energy and area arrays
- TH1F* hEtJet[30];
- TH1F* hAreaJet[30];
+ char hEtname[256];char hAreaname[256];
+
for(Int_t mjet=0; mjet<nJ; mjet++){
- char hEtname[256]; char hAreaname[256];
- sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "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);
- }
+ if(!fhEtJet[mjet]){
+ sprintf(hEtname, "hEtJet%d", mjet);
+ fhEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
+ }
+ if(!fhAreaJet[mjet]){
+ sprintf(hAreaname, "hAreaJet%d", mjet);
+ fhAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
+ }
+ fhEtJet[mjet]->Reset();
+ fhAreaJet[mjet]->Reset();
+ }
// background energy and area
- TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
- TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+ if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
+ fhEtBackg->Reset();
+ if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+ fhAreaBackg->Reset();
//fill energies
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
injet[jpart] = ijet;
multJet[ijet]++;
if((fReader->GetCutFlag(jpart)) == 1){// pt cut
- hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
+ fhEtJet[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] == -1 && fReader->GetCutFlag(jpart) == 1)
- hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
//calc areas
if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
}
- hAreaJet[ijet]->Fill(etac,areaj);
+ fhAreaJet[ijet]->Fill(etac,areaj);
areabg = areabg - areaj;
} // end jets loop
- hAreaBackg->Fill(etac,areabg);
+ fhAreaBackg->Fill(etac,areabg);
eta0 = eta1;
eta1 = eta1 + etaw;
} // end loop for all eta bins
for(Int_t kjet=0; kjet<nJ; kjet++){
etJet[kjet] = 0.0; // first clear etJet for this jet
for(Int_t bin = 0; bin< ndiv; bin++){
- if(hAreaJet[kjet]->GetBinContent(bin)){
- Float_t areab = hAreaBackg->GetBinContent(bin);
- Float_t etb = hEtBackg->GetBinContent(bin);
- Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
- etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
+ if(fhAreaJet[kjet]->GetBinContent(bin)){
+ Float_t areab = fhAreaBackg->GetBinContent(bin);
+ Float_t etb = fhEtBackg->GetBinContent(bin);
+ Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
+ etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
}
}
}
// calc background total
- Double_t etOut = hEtBackg->Integral();
- Double_t areaOut = hAreaBackg->Integral();
+ Double_t etOut = fhEtBackg->Integral();
+ Double_t areaOut = fhAreaBackg->Integral();
Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
-
- //delete
- for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
- delete hEtJet[ljet];
- delete hAreaJet[ljet];
- }
-
- delete hEtBackg;
- delete hAreaBackg;
}
////////////////////////////////////////////////////////////////////////
Int_t ndiv = 100;
// jet energy and area arrays
- TH1F* hEtJet[30];
- TH1F* hAreaJet[30];
+ // jet energy and area arrays
+ char hEtname[256];char hAreaname[256];
+
for(Int_t mjet=0; mjet<nJ; mjet++){
- char hEtname[256]; char hAreaname[256];
- sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "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
- }
+ if(!fhEtJet[mjet]){
+ sprintf(hEtname, "hEtJet%d", mjet);
+ fhEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
+ }
+ if(!fhAreaJet[mjet]){
+ sprintf(hAreaname, "hAreaJet%d", mjet);
+ fhAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
+ }
+ fhEtJet[mjet]->Reset();
+ fhAreaJet[mjet]->Reset();
+ }
// background energy and area
- TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range
- TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range
+ if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
+ fhEtBackg->Reset();
+ if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+ fhAreaBackg->Reset();
//fill energies
for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
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
+ fhEtJet[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] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
//calc areas
if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
}
- hAreaJet[ijet]->Fill(etac,areaj);
+ fhAreaJet[ijet]->Fill(etac,areaj);
areabg = areabg - areaj;
} // end jets loop
- hAreaBackg->Fill(etac,areabg);
+ fhAreaBackg->Fill(etac,areabg);
eta0 = eta1;
eta1 = eta1 + etaw;
} // end loop for all eta bins
for(Int_t kjet=0; kjet<nJ; kjet++){
etJet[kjet] = 0.0; // first clear etJet for this jet
for(Int_t bin = 0; bin< ndiv; bin++){
- if(hAreaJet[kjet]->GetBinContent(bin)){
- Float_t areab = hAreaBackg->GetBinContent(bin);
- Float_t etb = hEtBackg->GetBinContent(bin);
- Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
- etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
+ if(fhAreaJet[kjet]->GetBinContent(bin)){
+ Float_t areab = fhAreaBackg->GetBinContent(bin);
+ Float_t etb = fhEtBackg->GetBinContent(bin);
+ Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
+ etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
}
}
}
// calc background total
- Double_t etOut = hEtBackg->Integral();
- Double_t areaOut = hAreaBackg->Integral();
+ Double_t etOut = fhEtBackg->Integral();
+ Double_t areaOut = fhAreaBackg->Integral();
Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
-
- //delete
- for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
- delete hEtJet[ljet];
- delete hAreaJet[ljet];
- }
-
- delete hEtBackg;
- delete hAreaBackg;
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinderV1::Reset()
{
fLego->Reset();
- fJets->ClearJets();
AliJetFinder::Reset();
}
////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV1::WriteJHeaderToFile()
+void AliUA1JetFinderV1::WriteJHeaderToFile() const
{
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
header->Write();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff