X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliUA1JetFinderV2.cxx;h=76a4743728ac3f4d0358e131e7a64e27b19ebec7;hb=b94f89e41d0ca20ed820faf513317da39cdcdf7f;hp=81da6b2d7351f625d1a4fb68990db94d448cb8cb;hpb=ee7de0dd4cd0478eb8d46cf642ff78304ecf8f51;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliUA1JetFinderV2.cxx b/JETAN/AliUA1JetFinderV2.cxx index 81da6b2d735..76a4743728a 100644 --- a/JETAN/AliUA1JetFinderV2.cxx +++ b/JETAN/AliUA1JetFinderV2.cxx @@ -16,145 +16,485 @@ /* $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 - -#include #include -#include #include #include #include #include #include +#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(): - fDebug(0), - fOpt(0) +//////////////////////////////////////////////////////////////////////// +AliUA1JetFinderV2::AliUA1JetFinderV2() : + AliJetFinder(), + fLego(0), + fOpt(0) { + // // Constructor - fHeader = 0x0; - fLego = 0x0; + // } //////////////////////////////////////////////////////////////////////// - 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(); + fDebug = fHeader->GetDebug(); + + 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[30]; // eta jet + Float_t phiJet[30]; // phi jet + Float_t etJet[30]; // et jet + Float_t etsigJet[30]; // signal et in jet + Float_t etallJet[30]; // total et in jet (tmp variable) + Int_t ncellsJet[30]; + Int_t multJet[30]; + //--- Added for jet reordering at the end of the jet finding procedure + Float_t etaJetOk[30]; + Float_t phiJetOk[30]; + Float_t etJetOk[30]; + Float_t etsigJetOk[30]; // signal et in jet + Float_t etallJetOk[30]; // total et in jet (tmp variable) + Int_t ncellsJetOk[30]; + Int_t multJetOk[30]; + //-------------------------- + Int_t nJets; // to hold number of jets found by algorithm + Int_t nj; // number of jets accepted + 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 (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[] idxjets; + delete[] idx; + + delete[] percentage; + delete[] ncells; + delete[] mult; + //-------------------------- + +} + +//////////////////////////////////////////////////////////////////////// +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(); 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; iAt(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; jpartAdd((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; jAt(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; jAt(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; jAt(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; jAt(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(); @@ -163,154 +503,238 @@ void AliUA1JetFinderV2::FindJets() 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]; - 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 etaJet[30]; + Float_t phiJet[30]; + Float_t etJet[30]; + Float_t etsigJet[30]; //signal et in jet + Float_t etallJet[30]; // total et in jet (tmp variable) + Int_t ncellsJet[30]; + Int_t multJet[30]; + //--- Added by me for jet reordering at the end of the jet finding procedure + Float_t etaJetOk[30]; + Float_t phiJetOk[30]; + Float_t etJetOk[30]; + Float_t etsigJetOk[30]; //signal et in jet + Float_t etallJetOk[30]; // total et in jet (tmp variable) + Int_t ncellsJetOk[30]; + Int_t multJetOk[30]; + //-------------------------- + Int_t nJets; // to hold number of jets found by algorithm + Int_t nj; // number of jets accepted + 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 (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 (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 etaT; - delete phiT; - delete detaT; - delete dphiT; - delete cFlagT; - delete cClusterT; - delete idT; - delete injet; - delete sflag; - delete hPtTotal; - delete etCell; - delete etaCell; - delete phiCell; - delete flagCell; - delete etaJet; - delete phiJet; - delete etJet; - delete etsigJet; - delete etallJet; - delete ncellsJet; - delete multJet; - delete idxjets; - delete percentage; - delete ncells; - delete mult; + delete [] ptT; + delete [] en2T; + delete [] pt2T; + delete [] etaT; + delete [] phiT; + delete [] detT; + delete [] cFlagT; + delete [] cFlag2T; + delete [] sFlagT; + delete [] cClusterT; + delete [] vectT; + delete [] injet; + delete [] sflag; + trackRef->Delete(); + delete trackRef; + delete hPtTotal; + delete [] etCell; + delete [] etaCell; + delete [] phiCell; + delete [] flagCell; + delete [] etCell2; + delete [] etaCell2; + delete [] phiCell2; + delete [] flagCell2; + //-------------------------- + delete [] idxjets; + delete [] idx; + delete [] trackinjet; + + delete [] percentage; + delete [] ncells; + delete [] mult; } //////////////////////////////////////////////////////////////////////// - -void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, - Int_t* flagCell, Float_t etbgTotal, Double_t dEtTotal, - Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, - Float_t* etallJet, Int_t* ncellsJet) +void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* const etaCell, Float_t* phiCell, + Int_t* const flagCell, const Float_t* etCell2, const Float_t* etaCell2, const Float_t* phiCell2, + const Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, + Int_t& nJets, Float_t* const etJet, Float_t* const etaJet, Float_t* const phiJet, + Float_t* const etallJet, Int_t* const ncellsJet) { + // + // Main method for jet finding + // UA1 base cone finder + // - Int_t nCell = nIn; + Int_t nCell = nIn; + // Dump lego + // Check enough space! *to be done* AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + for(Int_t i=0; iGetMinMove(); @@ -318,561 +742,974 @@ void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell 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 - 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; + // 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; - //calculate dr - deta = etaCell[lcell] - eta; - dphi = phiCell[lcell] - phi; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ - // calculate offset from initiate cell - deta = etaCell[lcell] - eta0; - dphi = phiCell[lcell] - phi0; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - 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 - dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib)); - 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; - } - } + 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; - 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* const etJet,Float_t* const etaJet, Float_t* const phiJet, + Float_t* const etallJet, Int_t* const ncellsJet) +{ + // Dump lego + // Check enough space! *to be done* + 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(); + + // 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; - } // end of cells loop + 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; + 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* const etJet, const Float_t* etaJet, const Float_t* phiJet, + Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet) { - //background subtraction using cone method but without correction in dE/deta distribution - + // + // 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(); 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 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 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; kGetLegoEtaMax())*TMath::Pi(); + + for(Int_t k=0; k 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; ljetGetLegoEtaMax())*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; kGetLegoEtaMax(); + 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; ljetGetLegoEtaMax()*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* const etJet, const Float_t* etaJet, const Float_t* phiJet, + Float_t* const etsigJet,Int_t* const multJet, Int_t* const injet) +{ + //background subtraction using cone method but without correction in dE/deta distribution + + //calculate energy inside and outside cones + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Float_t rc= header->GetRadius(); + Float_t etIn[30]; + Float_t 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 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 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; ljetGetLegoEtaMax())*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* const etJet, const Float_t* etaJet, const Float_t* phiJet, + Float_t* const etsigJet, Int_t* const multJet, Int_t* const injet) { //background subtraction using statistical method 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 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 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 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; ljetGetRadius(); - 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; mjetGetRadius(); + 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 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 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 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(eta1Fill(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; kjetGetBinContent(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 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(eta1Fill(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; kjetGetBinContent(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; ljetGetBackgCutRatio(); - - - //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; mjetGetBackgCutRatio(); + + //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; mjetGetCutFlag(jpart)) != 1) continue; - for(Int_t ijet=0; ijet 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 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 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(eta1Fill(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; kjetGetBinContent(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 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(eta1Fill(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; kjetGetBinContent(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; ljetReset(); - 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()); - - fDebug = fReader->GetReaderHeader()->GetDebug(); + // book lego + fLego = new TH2F("legoH","eta-phi", + header->GetLegoNbinEta(), header->GetLegoEtaMin(), + header->GetLegoEtaMax(), header->GetLegoNbinPhi(), + header->GetLegoPhiMin(), header->GetLegoPhiMax()); + + fDebug = fHeader->GetDebug(); fOpt = fReader->GetReaderHeader()->GetDetector(); - + + // Tasks initialization if(fOpt>0) - fReader->CreateTasks(); + fReader->CreateTasks(tree); }