/* $Id$ */
//---------------------------------------------------------------------
-// UA1 Cone Algorithm Jet finder
-// manages the search for jets
-// Author: Rafael.Diaz.Valdes@cern.ch
-// (version in c++)
-// Modified to include neutral particles (magali.estienne@ires.in2p3.fr)
+// UA1 Cone Algorithm Jet finder for charged + neutral jet studies
+// manages the search for jets using charged particle momentum and
+// neutral cell energy information
+// Based on UA1 V1 (from R. Diaz)
+// Author: magali.estienne@subatech.in2p3.fr
//---------------------------------------------------------------------
-#include <Riostream.h>
-
-#include <TArrayF.h>
#include <TClonesArray.h>
-#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <TRefArray.h>
+#include "TFile.h"
#include "AliUA1JetFinderV2.h"
#include "AliUA1JetHeaderV1.h"
#include "AliJetUnitArray.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
-#include "AliJet.h"
-#include "AliAODJet.h"
+#include "AliJetHeader.h"
+class TArrayF;
+class TFile;
+class AliJetReader;
+class AliAODJet;
ClassImp(AliUA1JetFinderV2)
-////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
AliUA1JetFinderV2::AliUA1JetFinderV2() :
- AliJetFinder(),
- fLego(0),
- fDebug(0),
- fOpt(0)
+ AliJetFinder(),
+ fLego(0),
+ fOpt(0)
{
+ //
// Constructor
+ //
}
////////////////////////////////////////////////////////////////////////
-
AliUA1JetFinderV2::~AliUA1JetFinderV2()
-
{
- // destructor
+ //
+ // Destructor
+ //
}
////////////////////////////////////////////////////////////////////////
+void AliUA1JetFinderV2::FindJetsC()
+{
+ //
+ // Used to find jets using charged particle momentum information
+ //
+ // 1) Fill cell map array
+ // 2) calculate total energy and fluctuation level
+ // 3) Run algorithm
+ // 3.1) look centroides in cell map
+ // 3.2) calculate total energy in cones
+ // 3.3) flag as a possible jet
+ // 3.4) reorder cones by energy
+ // 4) subtract backg in accepted jets
+ // 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];
+ 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];
+
+ //total energy in array
+ Float_t etbgTotal = 0.0;
+ TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.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();
+ etaT[i] = lv->Eta();
+ phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+ cFlagT[i] = fReader->GetCutFlag(i);
+ sFlagT[i] = fReader->GetSignalFlag(i);
+
+ if (fReader->GetCutFlag(i) != 1) continue;
+ fLego->Fill(etaT[i], phiT[i], ptT[i]);
+ hPtTotal->Fill(ptT[i]);
+ etbgTotal+= ptT[i];
+ }
+
+
+ // calculate total energy and fluctuation in map
+ Double_t meanpt = hPtTotal->GetMean();
+ 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
+ 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];
+ //--- 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];
+ //--------------------------
+ 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);
+ //--- Added for jet reordering at the end of the jet finding procedure
+ memset(etaJetOk,0,sizeof(Float_t)*30);
+ memset(phiJetOk,0,sizeof(Float_t)*30);
+ memset(etJetOk,0,sizeof(Float_t)*30);
+ memset(etallJetOk,0,sizeof(Float_t)*30);
+ memset(etsigJetOk,0,sizeof(Float_t)*30);
+ memset(ncellsJetOk,0,sizeof(Int_t)*30);
+ memset(multJetOk,0,sizeof(Int_t)*30);
+ //--------------------------
+ nJets = 0;
+ nj = 0;
+
+ // reset particles-jet array in memory
+ memset(injet,-1,sizeof(Int_t)*nIn);
+ //run cone algorithm finder
+ RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
+
+ //run background subtraction
+ 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(header->GetBackgMode() == 1) // standard
+ SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ if(header->GetBackgMode() == 2) //cone
+ SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ if(header->GetBackgMode() == 3) //ratio
+ SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ if(header->GetBackgMode() == 4) //statistic
+ SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ //calc precision
+ if(TMath::Abs(etbgTotalN) > 0.001)
+ bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
+ else
+ bgprec = 0;
+ etbgTotal = etbgTotalN; // update with new background estimation
+ } //end while
+
+ // add jets to list
+ Int_t* idxjets = new Int_t[nj];
+ Int_t nselectj = 0;
+ printf("Found %d jets \n", nj);
+ // Reorder jets by et in cone
+ Int_t * idx = new Int_t[nJets];
+ TMath::Sort(nJets, etJet, idx);
+ for(Int_t p = 0; p < nJets; p++){
+ etaJetOk[p] = etaJet[idx[p]];
+ phiJetOk[p] = phiJet[idx[p]];
+ etJetOk[p] = etJet[idx[p]];
+ etallJetOk[p] = etJet[idx[p]];
+ etsigJetOk[p] = etsigJet[idx[p]];
+ ncellsJetOk[p] = ncellsJet[idx[p]];
+ multJetOk[p] = multJet[idx[p]];
+ }
+
+ for(Int_t kj=0; kj<nj; kj++)
+ {
+ if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
+ (etaJetOk[kj] < (header->GetJetEtaMin())) ||
+ (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
+ Float_t px, py,pz,en; // convert to 4-vector
+ px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
+ py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
+ pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
+ en = TMath::Sqrt(px * px + py * py + pz * pz);
+
+ AliAODJet jet(px, py, pz, en);
+ jet.Print("");
+
+ AddJet(jet);
+
+ idxjets[nselectj] = kj;
+ nselectj++;
+ }
-void AliUA1JetFinderV2::FindJets()
+ //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] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
+ ncells[i] = ncellsJetOk[idxjets[i]];
+ mult[i] = multJetOk[idxjets[i]];
+ }
+
+ //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[] cFlagT;
+ delete[] sFlagT;
+ 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;
+ //--- Added for jet reordering
+ delete etaJetOk;
+ delete phiJetOk;
+ delete etJetOk;
+ delete etsigJetOk;
+ delete etallJetOk;
+ delete ncellsJetOk;
+ delete multJetOk;
+ //--------------------------
+
+}
+////////////////////////////////////////////////////////////////////////
+void AliUA1JetFinderV2::FindJets()
{
- //1) Fill cell map array
- //2) calculate total energy and fluctuation level
- //3) Run algorithm
- // 3.1) look centroides in cell map
- // 3.2) calculate total energy in cones
- // 3.3) flag as a possible jet
- // 3.4) reorder cones by energy
- //4) subtract backg in accepted jets
- //5) fill AliJet list
+ //
+ // Used to find jets using charged particle momentum information
+ // & neutral energy from calo cells
+ //
+ // 1) Fill cell map array
+ // 2) calculate total energy and fluctuation level
+ // 3) Run algorithm
+ // 3.1) look centroides in cell map
+ // 3.2) calculate total energy in cones
+ // 3.3) flag as a possible jet
+ // 3.4) reorder cones by energy
+ // 4) subtract backg in accepted jets
+ // 5) fill AliJet list
// transform input to pt,eta,phi plus lego
- AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
- TClonesArray* fUnit = fReader->GetUnitArray();
- Int_t nCandidate = fReader->GetNumCandidate();
- Int_t nIn = fUnit->GetEntries();
-
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ TClonesArray* fUnit = fReader->GetUnitArray();
+ Int_t nCand = fReader->GetNumCandidate();
+ Int_t nCandCut = fReader->GetNumCandidateCut();
+ 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 nCandidate = 0;
+
+ nCandidate = nCand;
+ nCandidateCut = nCandCut;
+
// local arrays for input No Cuts
// Both pt < ptMin and pt > ptMin
- Float_t* enT = new Float_t[nCandidate];
- Float_t* ptT = new Float_t[nCandidate];
- Float_t* etaT = new Float_t[nCandidate];
- Float_t* phiT = new Float_t[nCandidate];
- Float_t* detaT = new Float_t[nCandidate];
- Float_t* dphiT = new Float_t[nCandidate];
- Float_t* cFlagT = new Float_t[nCandidate];
- Float_t* cClusterT = new Float_t[nCandidate];
- Float_t* idT = new Float_t[nCandidate];
- Int_t loop1 = 0;
- Int_t* injet = new Int_t[nCandidate];
- Int_t* sflag = new Int_t[nCandidate];
-
+ Float_t* ptT = new Float_t[nCandidate];
+ Float_t* en2T = new Float_t[nCandidate];
+ Float_t* pt2T = new Float_t[nCandidate];
+ Int_t* detT = new Int_t[nCandidate];
+ Float_t* etaT = new Float_t[nCandidate];
+ Float_t* phiT = 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];
+ TRefArray* trackRef = new TRefArray();
//total energy in array
Float_t etbgTotal = 0.0;
- TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
+ TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
// Input cell info
- Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
- Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
- Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
- Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
-
+ Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
+ Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
+ Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
+ Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
+ Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
+ 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
+
// Information extracted from fUnitArray
- // load input vectors and calculate total energy in array
+ // Load input vectors and calculate total energy in array
for(Int_t i=0; i<nIn; i++)
{
+ // Recover particle information from UnitArray
+
AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
- if(uArray->GetUnitCutFlag()==1){
- etCell[i] = uArray->GetUnitEnergy();
- if (etCell[i] > 0.0) etCell[i] -= header->GetMinCellEt();
- if (etCell[i] < 0.0) etCell[i] = 0.;
- etaCell[i] = uArray->GetUnitEta();
- phiCell[i] = uArray->GetUnitPhi();
- flagCell[i] = 0; // default
- }
- else {
- etCell[i] = 0.;
- etaCell[i] = uArray->GetUnitEta();
- phiCell[i] = uArray->GetUnitPhi();
- flagCell[i] = 0;
- }
+ TRefArray* ref = uArray->GetUnitTrackRef();
+ Int_t nRef = ref->GetEntries();
if(uArray->GetUnitEnergy()>0.){
+
+ for(Int_t jpart=0; jpart<nRef;jpart++)
+ trackRef->Add((AliVTrack*)ref->At(jpart));
ptT[loop1] = uArray->GetUnitEnergy();
- enT[loop1] = uArray->GetUnitEnergy();
+ detT[loop1] = uArray->GetUnitDetectorFlag();
etaT[loop1] = uArray->GetUnitEta();
phiT[loop1] = uArray->GetUnitPhi();
- detaT[loop1] = uArray->GetUnitDeta();
- dphiT[loop1] = uArray->GetUnitDphi();
- cFlagT[loop1]= uArray->GetUnitCutFlag();
- idT[loop1] = uArray->GetUnitID();
- if(cFlagT[loop1] == 1) {
- hPtTotal->Fill(ptT[loop1]);
- // fLego->Fill(etaT[i], phiT[i], ptT[i]);
- etbgTotal+= ptT[loop1];
+ cFlagT[loop1]= uArray->GetUnitCutFlag(); // pt cut tpc
+ cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
+ sFlagT[loop1]= uArray->GetUnitSignalFlag();
+ vectT[loop1] = nRef;
+ if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
+ pt2T[loop1] = 0.;
+ en2T[loop1] = 0.;
+ if(detT[loop1]==1){
+ en2T[loop1] = ptT[loop1] - header->GetMinCellEt();
+ if(en2T[loop1] < 0) en2T[loop1]=0;
+ hPtTotal->Fill(en2T[loop1]);
+ etbgTotal += en2T[loop1];
+ }
+ if(detT[loop1]==0){ // TPC+ITS
+ Float_t pt = 0.;
+ for(Int_t j=0; j<nRef;j++){
+ Float_t x=0.; Float_t y=0.; Float_t z=0.;
+ x = ((AliVTrack*)ref->At(j))->Px();
+ y = ((AliVTrack*)ref->At(j))->Py();
+ z = ((AliVTrack*)ref->At(j))->Pz();
+ pt = TMath::Sqrt(x*x+y*y);
+ if(pt>ptMin) {
+ pt2T[loop1] += pt;
+ en2T[loop1] += pt;
+ hPtTotal->Fill(pt);
+ etbgTotal+= pt;
+ }
+ }
+ }
+ if(detT[loop1]==2) { // EMCal
+ Float_t ptCTot = 0.;
+ Float_t pt = 0.;
+ Float_t enC = 0.;
+ for(Int_t j=0; j<nRef;j++){
+ Float_t x=0.; Float_t y=0.; Float_t z=0.;
+ x = ((AliVTrack*)ref->At(j))->Px();
+ y = ((AliVTrack*)ref->At(j))->Py();
+ z = ((AliVTrack*)ref->At(j))->Pz();
+ pt = TMath::Sqrt(x*x+y*y);
+ if(pt>ptMin) {
+ pt2T[loop1]+=pt;
+ en2T[loop1]+=pt;
+ hPtTotal->Fill(pt);
+ etbgTotal+= pt;
+ }
+ ptCTot += pt;
+ }
+ enC = ptT[loop1] - ptCTot - header->GetMinCellEt();
+ if(enC < 0.) enC=0.;
+ en2T[loop1] += enC;
+ hPtTotal->Fill(enC);
+ etbgTotal+= enC;
+ }
}
loop1++;
}
- }
- // fJets->SetNinput(nIn);
- fJets->SetNinput(nCandidate);
+ if(uArray->GetUnitCutFlag()==1) {
+ if(uArray->GetUnitDetectorFlag()==1){ // EMCal case
+ etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt();
+ if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.;
+ etaCell[i] = uArray->GetUnitEta();
+ phiCell[i] = uArray->GetUnitPhi();
+ flagCell[i] = 0; // default
+ etCell2[i] = etCell[i];
+ etaCell2[i] = uArray->GetUnitEta();
+ phiCell2[i] = uArray->GetUnitPhi();
+ flagCell2[i] = 0; // default
+ }
+ if(uArray->GetUnitDetectorFlag()==0){ // TPC case
+ Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
+ for(Int_t j=0; j<nRef;j++)
+ {
+ Float_t x=0.; Float_t y=0.; Float_t z=0.;
+ x = ((AliVTrack*)ref->At(j))->Px();
+ y = ((AliVTrack*)ref->At(j))->Py();
+ z = ((AliVTrack*)ref->At(j))->Pz();
+ pt = TMath::Sqrt(x*x+y*y);
+ if(pt>ptMin) {
+ et1 += pt;
+ et2 += pt;
+ }
+ }
+ etCell[i] = et1;
+ etCell2[i] = et2;
+ if(et1 < 0.) etCell[i] = etCell2[i] = 0.;
+ etaCell[i] = uArray->GetUnitEta();
+ phiCell[i] = uArray->GetUnitPhi();
+ flagCell[i] = 0; // default
+ etaCell2[i] = uArray->GetUnitEta();
+ phiCell2[i] = uArray->GetUnitPhi();
+ flagCell2[i] = 0; // default
+ }
+ if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case
+ Float_t ptCTot = 0.;
+ Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
+ Float_t enC = 0.;
+ for(Int_t j=0; j<nRef;j++)
+ {
+ Float_t x=0.; Float_t y=0.; Float_t z=0.;
+ x = ((AliVTrack*)ref->At(j))->Px();
+ y = ((AliVTrack*)ref->At(j))->Py();
+ z = ((AliVTrack*)ref->At(j))->Pz();
+ pt = TMath::Sqrt(x*x+y*y);
+ if(pt>ptMin) {
+ et1 += pt;
+ et2 += pt;
+ }
+ ptCTot += pt;
+ }
+ enC = uArray->GetUnitEnergy() - ptCTot;
+ etCell[i] = et1 + enC - header->GetMinCellEt();
+ etCell2[i] = et2 + enC - header->GetMinCellEt();
+ if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.;
+ etaCell[i] = uArray->GetUnitEta();
+ phiCell[i] = uArray->GetUnitPhi();
+ flagCell[i] = 0; // default
+ etaCell2[i] = uArray->GetUnitEta();
+ phiCell2[i] = uArray->GetUnitPhi();
+ flagCell2[i] = 0; // default
+ }
+ }
+ else {
+ etCell[i] = 0.;
+ etaCell[i] = uArray->GetUnitEta();
+ phiCell[i] = uArray->GetUnitPhi();
+ flagCell[i] = 0;
+ etCell2[i] = 0.;
+ etaCell2[i] = uArray->GetUnitEta();
+ phiCell2[i] = uArray->GetUnitPhi();
+ flagCell2[i] = 0;
+ }
+ } // end loop on nCandidate
+
// calculate total energy and fluctuation in map
Double_t meanpt = hPtTotal->GetMean();
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* 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];
- 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;
+ 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];
+ //--- 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];
+ //--------------------------
+ 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);
- nJets = 0;
- nj = 0;
- // reset particles-jet array in memory
- memset(injet,-1,sizeof(Int_t)*nCandidate);
- //run cone algorithm finder
- RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
- //run background subtraction
- 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(header->GetBackgMode() == 1) // standar
- SubtractBackg(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(header->GetBackgMode() == 2) //cone
- SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(header->GetBackgMode() == 3) //ratio
- SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- if(header->GetBackgMode() == 4) //statistic
- SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
- //calc precision
- if(etbgTotalN != 0.0)
- bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
- else
- bgprec = 0;
- etbgTotal = etbgTotalN; // update with new background estimation
- } //end while
+ //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);
+ //--- Added by me for jet reordering at the end of the jet finding procedure
+ memset(etaJetOk,0,sizeof(Float_t)*30);
+ memset(phiJetOk,0,sizeof(Float_t)*30);
+ memset(etJetOk,0,sizeof(Float_t)*30);
+ memset(etallJetOk,0,sizeof(Float_t)*30);
+ memset(etsigJetOk,0,sizeof(Float_t)*30);
+ memset(ncellsJetOk,0,sizeof(Int_t)*30);
+ memset(multJetOk,0,sizeof(Int_t)*30);
+
+ nJets = 0;
+ nj = 0;
+
+ // reset particles-jet array in memory
+ memset(injet,-1,sizeof(Int_t)*nCandidate);
+ //run cone algorithm finder
+ RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2,
+ flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,
+ etallJet,ncellsJet);
+
+ //run background subtraction
+ 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(header->GetBackgMode() == 1) // standard
+ SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ // To be modified ------------------------
+ if(header->GetBackgMode() == 2) //cone
+ SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ if(header->GetBackgMode() == 3) //ratio
+ SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ if(header->GetBackgMode() == 4) //statistic
+ SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+ //----------------------------------------
+ //calc precision
+ if(etbgTotalN != 0.0)
+ bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
+ else
+ bgprec = 0;
+ etbgTotal = etbgTotalN; // update with new background estimation
+ } //end while
+
// 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())) ||
- (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
+ // Reorder jets by et in cone
+ // Sort jets by energy
+ Int_t * idx = new Int_t[nJets];
+ TMath::Sort(nJets, etJet, idx);
+ for(Int_t p = 0; p < nJets; p++)
+ {
+ etaJetOk[p] = etaJet[idx[p]];
+ phiJetOk[p] = phiJet[idx[p]];
+ etJetOk[p] = etJet[idx[p]];
+ etallJetOk[p] = etJet[idx[p]];
+ etsigJetOk[p] = etsigJet[idx[p]];
+ ncellsJetOk[p] = ncellsJet[idx[p]];
+ 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())) ||
+ (etaJetOk[kj] < (header->GetJetEtaMin())) ||
+ (etJetOk[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])));
+ px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
+ py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
+ pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
en = TMath::Sqrt(px * px + py * py + pz * pz);
- fJets->AddJet(px, py, pz, en);
+
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;
nselectj++;
- }
+ }
+
//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]];
+ for(Int_t i = 0; i< nselectj; i++)
+ {
+ percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
+ ncells[i] = ncellsJetOk[idxjets[i]];
+ mult[i] = multJetOk[idxjets[i]];
+ }
+
+ //add particle-injet relationship ///
+ for(Int_t bj = 0; bj < nCandidate; 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
}
- //add particle-injet relationship ///
- // for(Int_t bj = 0; bj < nIn; bj++){
- for(Int_t bj = 0; bj < nCandidate; 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
- }
- 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 enT;
delete ptT;
+ delete en2T;
+ delete pt2T;
delete etaT;
delete phiT;
- delete detaT;
- delete dphiT;
+ trackRef->Delete();
+ delete trackRef;
+ delete detT;
delete cFlagT;
+ delete cFlag2T;
+ delete sFlagT;
delete cClusterT;
- delete idT;
+ delete vectT;
delete injet;
delete sflag;
delete hPtTotal;
delete etaCell;
delete phiCell;
delete flagCell;
+ delete etCell2;
+ delete etaCell2;
+ delete phiCell2;
+ delete flagCell2;
delete etaJet;
delete phiJet;
delete etJet;
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 trackinjet;
delete idxjets;
delete percentage;
delete ncells;
delete mult;
-
}
////////////////////////////////////////////////////////////////////////
-
void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
- Int_t* flagCell, Float_t etbgTotal, Double_t dEtTotal,
- Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
- Float_t* etallJet, Int_t* ncellsJet)
+ 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)
{
+ //
+ // Main method for jet finding
+ // UA1 base cone finder
+ //
- Int_t nCell = nIn;
+ Int_t nCell = nIn;
+ fDebug = fReader->GetReaderHeader()->GetDebug();
+ // Dump lego
+ // Check enough space! *to be done*
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ for(Int_t i=0; i<nCell; i++){
+ etCell[i] = etCell2[i];
+ etaCell[i] = etaCell2[i];
+ phiCell[i] = phiCell2[i];
+ flagCell[i] = flagCell2[i];
+ }
// Parameters from header
Float_t minmove = header->GetMinMove();
Float_t rc = header->GetRadius();
Float_t etseed = header->GetEtSeed();
- // tmp array of jets form algoritm
+ // Tmp array of jets form algoritm
Float_t etaAlgoJet[30];
Float_t phiAlgoJet[30];
Float_t etAlgoJet[30];
Int_t ncellsAlgoJet[30];
- //run algorithm//
+ // Run algorithm//
- // sort cells by et
+ // Sort cells by et
Int_t * index = new Int_t[nCell];
TMath::Sort(nCell, etCell, index);
- // variable used in centroide loop
+ // Variable used in centroide loop
Float_t eta = 0.0;
Float_t phi = 0.0;
Float_t eta0 = 0.0;
Float_t phisb = 0.0;
Float_t dphib = 0.0;
-
- for(Int_t icell = 0; icell < nCell; icell++){
+ for(Int_t icell = 0; icell < nCell; icell++)
+ {
Int_t jcell = index[icell];
if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
if(flagCell[jcell] != 0) continue; // if cell was used before
+
eta = etaCell[jcell];
phi = phiCell[jcell];
eta0 = eta;
etsb = ets;
etasb = 0.0;
phisb = 0.0;
- for(Int_t kcell =0; kcell < nCell; kcell++){
+ 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(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();
+ dphi = TMath::Abs(phiCell[lcell] - phi);
if (dphi > TMath::Pi()) dphi = 2.0 * 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 = dphi - 2.0 * TMath::Pi();
- etas = etas + etCell[lcell]*deta;
- phis = phis + etCell[lcell]*dphi;
- ets = ets + etCell[lcell];
- //new weighted eta and phi including this cell
- eta = eta0 + etas/ets;
- phi = phi0 + phis/ets;
- // if cone does not move much, just go to next step
- 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;
- }
- }
+ // 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 = dphi - 2.0 * TMath::Pi();
+ etas = etas + etCell[lcell]*deta;
+ phis = phis + etCell[lcell]*dphi;
+ ets = ets + etCell[lcell];
+ //new weighted eta and phi including this cell
+ eta = eta0 + etas/ets;
+ phi = phi0 + phis/ets;
+ // if cone does not move much, just go to next step
+ 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;
+ }
+ } // 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 = header->GetRadius();
- for(Int_t ncell =0; ncell < nCell; ncell++){
- if(flagCell[ncell] != 0) continue; // cell used before
- //calculate dr
- deta = etaCell[ncell] - eta;
- dphi = phiCell[ncell] - 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);
- if(dr <= rc){ // cell in cone
- flagCell[ncell] = -1;
- etCone+=etCell[ncell];
- nCellIn++;
- }
+ Float_t etCone = 0.0;
+ Int_t nCellIn = 0;
+ Int_t nCellOut = 0;
+ rc = header->GetRadius();
+
+ for(Int_t ncell =0; ncell < nCell; ncell++)
+ {
+ if(flagCell[ncell] != 0) continue; // cell used before
+ //calculate dr
+ deta = etaCell[ncell] - eta;
+ // if(deta <= rc){ // Added to improve velocity -> to be tested
+ dphi = phiCell[ncell] - phi;
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ // if(dphi <= rc){ // Added to improve velocity -> to be tested
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if(dr <= rc){ // cell in cone
+ flagCell[ncell] = -1;
+ etCone+=etCell[ncell];
+ nCellIn++;
+ }
+ else nCellOut++;
+ // } // end deta <= rc
+ // } // end dphi <= rc
}
- // select jets with et > background
- // estimate max fluctuation of background in cone
- Double_t ncellin = (Double_t)nCellIn;
- Double_t ntcell = (Double_t)nCell;
- Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
- // min cone et
- Double_t etcmin = etCone ; // could be used etCone - etmin !!
- //desicions !! etbmax < etcmin
- for(Int_t mcell =0; mcell < nCell; mcell++){
- if(flagCell[mcell] == -1){
- if(etbmax < etcmin)
- flagCell[mcell] = 1; //flag cell as used
- else
- flagCell[mcell] = 0; // leave it free
- }
- }
- //store tmp jet info !!!
- if(etbmax < etcmin) {
- etaAlgoJet[nJets] = eta;
- phiAlgoJet[nJets] = phi;
- etAlgoJet[nJets] = etCone;
- ncellsAlgoJet[nJets] = nCellIn;
- nJets++;
+ // select jets with et > background
+ // estimate max fluctuation of background in cone
+ Double_t ncellin = (Double_t)nCellIn;
+ Double_t ntcell = (Double_t)nCell;
+ Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/(ntcell));
+ // min cone et
+ Double_t etcmin = etCone ; // could be used etCone - etmin !!
+ //decisions !! etbmax < etcmin
+
+ for(Int_t mcell =0; mcell < nCell; mcell++)
+ {
+ if(flagCell[mcell] == -1){
+ if(etbmax < etcmin)
+ flagCell[mcell] = 1; //flag cell as used
+ else
+ flagCell[mcell] = 0; // leave it free
+ }
}
+ //store tmp jet info !!!
+ if(etbmax < etcmin)
+ {
+ etaAlgoJet[nJets] = eta;
+ phiAlgoJet[nJets] = phi;
+ etAlgoJet[nJets] = etCone;
+ ncellsAlgoJet[nJets] = nCellIn;
+ nJets++;
+ }
+
+ } // end of cells loop
+
+ for(Int_t p = 0; p < nJets; p++)
+ {
+ etaJet[p] = etaAlgoJet[p];
+ phiJet[p] = phiAlgoJet[p];
+ etJet[p] = etAlgoJet[p];
+ etallJet[p] = etAlgoJet[p];
+ ncellsJet[p] = ncellsAlgoJet[p];
+ }
+
+ //delete
+ 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)
+{
+ // 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
+
+ Int_t nCell = 0;
+ TAxis* xaxis = fLego->GetXaxis();
+ TAxis* yaxis = fLego->GetYaxis();
+ Float_t e = 0.0;
+ for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++)
+ {
+ for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++)
+ {
+ e = fLego->GetBinContent(i,j);
+ if (e < 0.0) continue; // don't include this cells
+ Float_t eta = xaxis->GetBinCenter(i);
+ Float_t phi = yaxis->GetBinCenter(j);
+ etCell[nCell] = e;
+ etaCell[nCell] = eta;
+ phiCell[nCell] = phi;
+ flagCell[nCell] = 0; //default
+ nCell++;
+ }
+ }
+
+ // Parameters from header
+ Float_t minmove = header->GetMinMove();
+ Float_t maxmove = header->GetMaxMove();
+ Float_t rc = header->GetRadius();
+ Float_t etseed = header->GetEtSeed();
- } // end of cells loop
+ // Tmp array of jets form algoritm
+ Float_t etaAlgoJet[30];
+ Float_t phiAlgoJet[30];
+ Float_t etAlgoJet[30];
+ Int_t ncellsAlgoJet[30];
+
+ // Run algorithm//
+
+ // Sort cells by et
+ 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 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];
+ if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
+ if(flagCell[jcell] != 0) continue; // if cell was used before
+
+ eta = etaCell[jcell];
+ phi = phiCell[jcell];
+ eta0 = eta;
+ phi0 = phi;
+ etab = eta;
+ phib = phi;
+ ets = etCell[jcell];
+ etas = 0.0;
+ phis = 0.0;
+ etsb = ets;
+ etasb = 0.0;
+ 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; // can this happen
+ //calculate dr
+ deta = etaCell[lcell] - eta;
+ dphi = TMath::Abs(phiCell[lcell] - phi);
+ if (dphi > TMath::Pi()) dphi = 2.0 * 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 = dphi - 2.0 * TMath::Pi();
+ etas = etas + etCell[lcell]*deta;
+ phis = phis + etCell[lcell]*dphi;
+ ets = ets + etCell[lcell];
+ //new weighted eta and phi including this cell
+ eta = eta0 + etas/ets;
+ phi = phi0 + phis/ets;
+ // if cone does not move much, just go to next step
+ 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;
+ }
+ } // 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;
+ Int_t nCellOut = 0;
+ rc = header->GetRadius();
+ for(Int_t ncell =0; ncell < nCell; ncell++)
+ {
+ if(flagCell[ncell] != 0) continue; // cell used before
+ //calculate dr
+ deta = etaCell[ncell] - eta;
+ dphi = phiCell[ncell] - 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);
+ if(dr <= rc){ // cell in cone
+ flagCell[ncell] = -1;
+ etCone+=etCell[ncell];
+ nCellIn++;
+ }
+ else nCellOut++;
+ }
+
+ // Select jets with et > background
+ // estimate max fluctuation of background in cone
+ Double_t ncellin = (Double_t)nCellIn;
+ Double_t ntcell = (Double_t)nCell;
+ Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
+ // min cone et
+ Double_t etcmin = etCone ; // could be used etCone - etmin !!
+ //decisions !! etbmax < etcmin
+
+ for(Int_t mcell =0; mcell < nCell; mcell++){
+ if(flagCell[mcell] == -1){
+ if(etbmax < etcmin)
+ flagCell[mcell] = 1; //flag cell as used
+ else
+ flagCell[mcell] = 0; // leave it free
+ }
+ }
+ //store tmp jet info !!!
+
+ if(etbmax < etcmin) {
+ etaAlgoJet[nJets] = eta;
+ phiAlgoJet[nJets] = phi;
+ etAlgoJet[nJets] = etCone;
+ ncellsAlgoJet[nJets] = nCellIn;
+ nJets++;
+ }
+
+ } // 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);
- for(Int_t p = 0; p < nJets; p++){
- etaJet[p] = etaAlgoJet[idx[p]];
- phiJet[p] = phiAlgoJet[idx[p]];
- etJet[p] = etAlgoJet[idx[p]];
- etallJet[p] = etAlgoJet[idx[p]];
- ncellsJet[p] = ncellsAlgoJet[idx[p]];
- }
-
-
+ for(Int_t p = 0; p < nJets; p++)
+ {
+ etaJet[p] = etaAlgoJet[idx[p]];
+ phiJet[p] = phiAlgoJet[idx[p]];
+ etJet[p] = etAlgoJet[idx[p]];
+ etallJet[p] = etAlgoJet[idx[p]];
+ ncellsJet[p] = ncellsAlgoJet[idx[p]];
+ }
+
//delete
delete index;
delete idx;
}
-////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::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 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* 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
-
+ //
+ // Background subtraction using cone method but without correction in dE/deta distribution
+ // Cases to take into account the EMCal geometry are included
+ //
+
//calculate energy inside and outside cones
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ fOpt = fReader->GetReaderHeader()->GetDetector();
+ fDebug = fReader->GetReaderHeader()->GetDebug();
Float_t rc= header->GetRadius();
Float_t etIn[30];
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
- // 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= -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
- 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] == -1 && fReader->GetCutFlag(jpart) == 1)
- etOut += ptT[jpart]; // particle outside cones and 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= -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
+ multJet[ijet]+=vectT[jpart];
+ injet[jpart] = ijet;
+
+ if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut
+ etIn[ijet] += ptT[jpart];
+ if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart];
+ }
+ break;
+ }
+ }// end jets loop
+
+ if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[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*(header->GetLegoEtaMax())*TMath::Pi();
- for(Int_t k=0; k<nJ; k++){
+ // TPC case
+ if(fOpt == 0 || fOpt == 1){
+ Float_t areaJet[30];
+ 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 > 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);
+ Float_t h = header->GetLegoEtaMax() - etaJet[k];
+ 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[k];
- accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ 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;
areaOut = areaOut - areaJet[k];
+ }
+ //subtract background using area method
+ for(Int_t ljet=0; ljet<nJ; ljet++){
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
+ }
+
+ // estimate new total background
+ Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
+ etbgTotalN = etOut*areaT/areaOut;
+ }
+ else { // If EMCal included
+ Float_t areaJet[30];
+ Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
+ for(Int_t k=0; k<nJ; k++){
+ Float_t detamax = etaJet[k] + rc;
+ Float_t detamin = etaJet[k] - rc;
+ Float_t dphimax = phiJet[k] + rc;
+ Float_t dphimin = phiJet[k] - rc;
+ Float_t eMax = header->GetLegoEtaMax();
+ Float_t eMin = header->GetLegoEtaMin();
+ Float_t pMax = header->GetLegoPhiMax();
+ Float_t pMin = header->GetLegoPhiMin();
+ Float_t accetamax = 0.0; Float_t accetamin = 0.0;
+ Float_t accphimax = 0.0; Float_t accphimin = 0.0;
+ if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
+ (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
+ Float_t h = eMax - etaJet[k];
+ accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+ if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
+ (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
+ Float_t h = eMax + etaJet[k];
+ accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+ if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )||
+ (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){
+ Float_t h = pMax - phiJet[k];
+ accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+ if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )||
+ (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){
+ Float_t h = phiJet[k] - pMin;
+ accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ }
+
+ if(detamax > eMax && dphimax > pMax ){
+ Float_t he = eMax - etaJet[k];
+ Float_t hp = pMax - phiJet[k];
+ Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
+ Float_t alphae = TMath::ACos(he/rc);
+ Float_t alphap = TMath::ACos(hp/rc);
+ Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
+ if(rlim <= rc){
+ accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
+ }
+ if(rlim > rc){
+ accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
+ ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
+ rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
+ }
+ }
+
+ if(detamax > eMax && dphimin < pMin ){
+ Float_t he = eMax - etaJet[k];
+ Float_t hp = phiJet[k] - pMin;
+ Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
+ Float_t alphae = TMath::ACos(he/rc);
+ Float_t alphap = TMath::ACos(hp/rc);
+ Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
+ if(rlim <= rc){
+ accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
+ }
+ if(rlim > rc){
+ accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
+ ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
+ rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
+ }
+ }
+
+ if(detamin < eMin && dphimax > pMax ){
+ Float_t he = eMax + etaJet[k];
+ Float_t hp = pMax - phiJet[k];
+ Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
+ Float_t alphae = TMath::ACos(he/rc);
+ Float_t alphap = TMath::ACos(hp/rc);
+ Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
+ if(rlim <= rc){
+ accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
+ }
+ if(rlim > rc){
+ accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
+ ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
+ rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
+ }
+ }
+
+ if(detamin < eMin && dphimin < pMin ){
+ Float_t he = eMax + etaJet[k];
+ Float_t hp = phiJet[k] - pMin;
+ Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
+ Float_t alphae = TMath::ACos(he/rc);
+ Float_t alphap = TMath::ACos(hp/rc);
+ Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
+ if(rlim <= rc){
+ accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
+ }
+ if(rlim > rc){
+ accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
+ accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
+ ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
+ rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
+ }
+ }
+ areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin;
+ areaOut = areaOut - areaJet[k];
+ } // end loop on jets
+
+ //subtract background using area method
+ for(Int_t ljet=0; ljet<nJ; ljet++){
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
+ }
+
+ // estimate new total background
+ Float_t areaT = 2*(header->GetLegoEtaMax()*header->GetLegoPhiMax());
+ etbgTotalN = etOut*areaT/areaOut;
+ }
+
+}
+
+////////////////////////////////////////////////////////////////////////
+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)
+{
+ //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 etOut = 0;
+ 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++){
+ 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
+ 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] == -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*(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 > 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 < 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;
+ areaOut = areaOut - areaJet[k];
}
//subtract background using area method
for(Int_t ljet=0; ljet<nJ; ljet++){
- Float_t areaRatio = areaJet[ljet]/areaOut;
- etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
}
-
+
// estimate new total background
Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
etbgTotalN = etOut*areaT/areaOut;
-
-
+
}
-////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::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 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 Int_t* cFlagT,
+ const Int_t* sFlagT, 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
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
-
+
//calculate energy inside
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
- 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
- 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
- } //end particle loop
-
+
+ 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
+ multJet[ijet]++;
+ injet[jpart] = ijet;
+ if(cFlagT[jpart] == 1){ // pt cut
+ etIn[ijet]+= ptT[jpart];
+ if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
+ }
+ break;
+ }
+ }// end jets loop
+ } //end particle loop
+
//calc jets areas
Float_t areaJet[30];
Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
- for(Int_t k=0; k<nJ; k++){
+ 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 > 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);
+ Float_t h = header->GetLegoEtaMax() - etaJet[k];
+ 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[k];
- accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+ 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;
- }
+ }
//subtract background using area method
for(Int_t ljet=0; ljet<nJ; ljet++){
- Float_t areaRatio = areaJet[ljet]/areaOut;
- etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
}
-
+
etbgTotalN = etbgStat;
-
}
////////////////////////////////////////////////////////////////////////
-
-void AliUA1JetFinderV2::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 AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
+ Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT, const Int_t* sFlagT,
+ 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= header->GetRadius();
- Float_t etamax = header->GetLegoEtaMax();
- Float_t etamin = header->GetLegoEtaMin();
- Int_t ndiv = 100;
-
- // jet energy and area arrays
- TH1F* hEtJet[30];
- 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);
- hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
- hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
+ // Cone background subtraction method taking into acount dEt/deta distribution
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ //general
+ Float_t rc= header->GetRadius();
+ Float_t etamax = header->GetLegoEtaMax();
+ Float_t etamin = header->GetLegoEtaMin();
+ Int_t ndiv = 100;
+
+ // jet energy and area arrays
+ TH1F* hEtJet[30];
+ 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);
+ hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
+ hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
}
- // 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);
-
- //fill energies
- for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- 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= -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;
- 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] == -1 && fReader->GetCutFlag(jpart) == 1)
- hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ // 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);
+
+ //fill energies
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ 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= -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;
+ multJet[ijet]++;
+ if(cFlagT[jpart] == 1){// pt cut
+ hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
+ if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
+ }
+ break;
+ }
+ }// end jets loop
+
+ if(injet[jpart] == -1 && cFlagT[jpart] == 1)
+ hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
- //calc areas
- Float_t eta0 = etamin;
- Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
- Float_t eta1 = eta0 + etaw;
- for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
- Float_t etac = eta0 + etaw/2.0;
- Float_t areabg = etaw*2.0*TMath::Pi();
- for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
- Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
- Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
- Float_t acc0 = 0.0; Float_t acc1 = 0.0;
- Float_t areaj = 0.0;
- if(deta0 > rc && deta1 < rc){
- acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
- areaj = acc1;
- }
- if(deta0 < rc && deta1 > rc){
- acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
- areaj = acc0;
- }
- if(deta0 < rc && deta1 < rc){
- acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
- acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
- if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
- 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);
- areabg = areabg - areaj;
- } // end jets loop
- hAreaBackg->Fill(etac,areabg);
- eta0 = eta1;
- eta1 = eta1 + etaw;
- } // end loop for all eta bins
-
- //subtract background
- 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
- }
- }
- }
-
- // calc background total
- Double_t etOut = hEtBackg->Integral();
- Double_t areaOut = hAreaBackg->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;
-}
+ //calc areas
+ Float_t eta0 = etamin;
+ Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
+ Float_t eta1 = eta0 + etaw;
+ for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
+ Float_t etac = eta0 + etaw/2.0;
+ Float_t areabg = etaw*2.0*TMath::Pi();
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+ Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+ Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+ Float_t areaj = 0.0;
+ if(deta0 > rc && deta1 < rc){
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ areaj = acc1;
+ }
+ if(deta0 < rc && deta1 > rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ areaj = acc0;
+ }
+ if(deta0 < rc && deta1 < rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
+ 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);
+ areabg = areabg - areaj;
+ } // end jets loop
+ hAreaBackg->Fill(etac,areabg);
+ eta0 = eta1;
+ eta1 = eta1 + etaw;
+ } // end loop for all eta bins
+
+ //subtract background
+ 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
+ }
+ }
+ }
-////////////////////////////////////////////////////////////////////////
+ // calc background total
+ Double_t etOut = hEtBackg->Integral();
+ Double_t areaOut = hAreaBackg->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 AliUA1JetFinderV2::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 AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
+ Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Int_t* cFlagT, const Int_t* sFlagT,
+ 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 = header->GetBackgCutRatio();
-
-
- //general
- Float_t rc= header->GetRadius();
- Float_t etamax = header->GetLegoEtaMax();
- Float_t etamin = header->GetLegoEtaMin();
- Int_t ndiv = 100;
-
- // jet energy and area arrays
- TH1F* hEtJet[30];
- 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);
- 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
+ // Ratio background subtraction method taking into acount dEt/deta distribution
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ //factor F calc before
+ Float_t bgRatioCut = header->GetBackgCutRatio();
+
+ //general
+ Float_t rc= header->GetRadius();
+ Float_t etamax = header->GetLegoEtaMax();
+ Float_t etamin = header->GetLegoEtaMin();
+ Int_t ndiv = 100;
+
+ // jet energy and area arrays
+ TH1F* hEtJet[30];
+ 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);
+ 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
}
- // 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
-
- //fill energies
- for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
- //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= -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
- 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] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ // 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
+
+ //fill energies
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ 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= -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
+ multJet[ijet]++;
+ injet[jpart] = ijet;
+ if(cFlagT[jpart] == 1){ //pt cut
+ hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
+ if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
+ }
+ break;
+ }
+ }// end jets loop
+ if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
} //end particle loop
- //calc areas
- Float_t eta0 = etamin;
- Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
- Float_t eta1 = eta0 + etaw;
- for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
- Float_t etac = eta0 + etaw/2.0;
- Float_t areabg = etaw*2.0*TMath::Pi();
- for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
- Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
- Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
- Float_t acc0 = 0.0; Float_t acc1 = 0.0;
- Float_t areaj = 0.0;
- if(deta0 > rc && deta1 < rc){
- acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
- areaj = acc1;
- }
- if(deta0 < rc && deta1 > rc){
- acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
- areaj = acc0;
- }
- if(deta0 < rc && deta1 < rc){
- acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
- acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
- if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
- 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);
- areabg = areabg - areaj;
- } // end jets loop
- hAreaBackg->Fill(etac,areabg);
- eta0 = eta1;
- eta1 = eta1 + etaw;
- } // end loop for all eta bins
-
- //subtract background
- 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
- }
- }
- }
-
- // calc background total
- Double_t etOut = hEtBackg->Integral();
- Double_t areaOut = hAreaBackg->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;
+ //calc areas
+ Float_t eta0 = etamin;
+ Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
+ Float_t eta1 = eta0 + etaw;
+ for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
+ Float_t etac = eta0 + etaw/2.0;
+ Float_t areabg = etaw*2.0*TMath::Pi();
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+ Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+ Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+ Float_t areaj = 0.0;
+ if(deta0 > rc && deta1 < rc){
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ areaj = acc1;
+ }
+ if(deta0 < rc && deta1 > rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ areaj = acc0;
+ }
+ if(deta0 < rc && deta1 < rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
+ 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);
+ areabg = areabg - areaj;
+ } // end jets loop
+ hAreaBackg->Fill(etac,areabg);
+ eta0 = eta1;
+ eta1 = eta1 + etaw;
+ } // end loop for all eta bins
+
+ //subtract background
+ 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
+ }
+ }
+ }
+
+ // calc background total
+ Double_t etOut = hEtBackg->Integral();
+ Double_t areaOut = hAreaBackg->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 AliUA1JetFinderV2::Reset()
{
fLego->Reset();
- fJets->ClearJets();
AliJetFinder::Reset();
}
////////////////////////////////////////////////////////////////////////
-
-void AliUA1JetFinderV2::WriteJHeaderToFile()
+void AliUA1JetFinderV2::WriteJHeaderToFile() const
{
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
header->Write();
}
////////////////////////////////////////////////////////////////////////
-
-void AliUA1JetFinderV2::Init()
+void AliUA1JetFinderV2::InitTask(TChain* tree)
{
+
// initializes some variables
AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
- // book lego
- fLego = new
- TH2F("legoH","eta-phi",
- header->GetLegoNbinEta(), header->GetLegoEtaMin(),
- header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
- header->GetLegoPhiMin(), header->GetLegoPhiMax());
-
+ // book lego
+ fLego = new TH2F("legoH","eta-phi",
+ header->GetLegoNbinEta(), header->GetLegoEtaMin(),
+ header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
+ header->GetLegoPhiMin(), header->GetLegoPhiMax());
+
fDebug = fReader->GetReaderHeader()->GetDebug();
fOpt = fReader->GetReaderHeader()->GetDetector();
-
+
+ // Tasks initialization
if(fOpt>0)
- fReader->CreateTasks();
+ fReader->CreateTasks(tree);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff