From 70e58892af2e25391f6199e3d3d340f9fa1ecf57 Mon Sep 17 00:00:00 2001 From: hristov Date: Tue, 13 Jun 2006 09:43:53 +0000 Subject: [PATCH] New version of the jet finder (Rafael) --- JETAN/AliUA1JetFinderV1.cxx | 794 ++++++++++++++++++++++++++++++++++++ JETAN/AliUA1JetFinderV1.h | 67 +++ JETAN/AliUA1JetHeaderV1.cxx | 70 ++++ JETAN/AliUA1JetHeaderV1.h | 89 ++++ JETAN/JetAnalysisLinkDef.h | 2 + JETAN/libJETAN.pkg | 3 +- 6 files changed, 1024 insertions(+), 1 deletion(-) create mode 100644 JETAN/AliUA1JetFinderV1.cxx create mode 100644 JETAN/AliUA1JetFinderV1.h create mode 100644 JETAN/AliUA1JetHeaderV1.cxx create mode 100644 JETAN/AliUA1JetHeaderV1.h diff --git a/JETAN/AliUA1JetFinderV1.cxx b/JETAN/AliUA1JetFinderV1.cxx new file mode 100644 index 00000000000..2b1030b8d58 --- /dev/null +++ b/JETAN/AliUA1JetFinderV1.cxx @@ -0,0 +1,794 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//--------------------------------------------------------------------- +// UA1 Cone Algorithm Jet finder +// manages the search for jets +// Author: Rafael.Diaz.Valdes@cern.ch +// (version in c++) +//--------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include "AliUA1JetFinderV1.h" +#include "AliUA1JetHeaderV1.h" +#include "AliJetReaderHeader.h" +#include "AliJetReader.h" +#include "AliJet.h" + + +ClassImp(AliUA1JetFinderV1) + +//////////////////////////////////////////////////////////////////////// + +AliUA1JetFinderV1::AliUA1JetFinderV1() + +{ + // Constructor + fHeader = 0x0; + fLego = 0x0; +} + +//////////////////////////////////////////////////////////////////////// + +AliUA1JetFinderV1::~AliUA1JetFinderV1() + +{ + // destructor +} + +//////////////////////////////////////////////////////////////////////// + + +void AliUA1JetFinderV1::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 + + // transform input to pt,eta,phi plus lego + 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* 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()); + if (fReader->GetCutFlag(i) != 1) continue; + fLego->Fill(etaT[i], phiT[i], ptT[i]); + hPtTotal->Fill(ptT[i]); + etbgTotal+= ptT[i]; + } + fJets->SetNinput(nIn); + + // calculate total energy and fluctuation in map + Double_t meanpt = hPtTotal->GetMean(); + Double_t ptRMS = hPtTotal->GetRMS(); + Double_t npart = hPtTotal->GetEntries(); + Double_t 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 = fHeader->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,0,sizeof(Int_t)*nIn); + //run cone algorithm finder + RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet); + //run background subtraction + if(nJets > fHeader->GetNAcceptJets()) // limited number of accepted jets per event + nj = fHeader->GetNAcceptJets(); + else + nj = nJets; + //subtract background + Float_t etbgTotalN = 0.0; //new background + if(fHeader->GetBackgMode() == 1) // standar + SubtractBackg(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(fHeader->GetBackgMode() == 2) //cone + SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(fHeader->GetBackgMode() == 3) //ratio + SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(fHeader->GetBackgMode() == 4) //statistic + SubtractBackgStat(nIn,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 + + // add jets to list + Int_t* idxjets = new Int_t[nj]; + Int_t nselectj = 0; + for(Int_t kj=0; kj fHeader->GetMinJetEt()){ // check if et jets > et min + Float_t px, py,pz,en; // convert to 4-vector + px = etJet[kj] * TMath::Cos(phiJet[kj]); + py = etJet[kj] * TMath::Sin(phiJet[kj]); + pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj]))); + en = TMath::Sqrt(px * px + py * py + pz * pz); + fJets->AddJet(px, py, pz, en); + 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]]; + } + fJets->SetNCells(ncells); + fJets->SetPtFromSignal(percentage); + fJets->SetMultiplicities(mult); + fJets->SetInJet(injet); + fJets->SetEtaIn(etaT); + fJets->SetPhiIn(phiT); + fJets->SetPtIn(ptT); + fJets->SetEtAvg(etbgTotal/(4*(fHeader->GetLegoEtaMax())*TMath::Pi())); + + + //delete + delete ptT; + delete etaT; + delete phiT; + delete injet; + delete hPtTotal; + delete etaJet; + delete phiJet; + delete etJet; + delete etsigJet; + delete etallJet; + delete ncellsJet; + delete multJet; + delete idxjets; + delete percentage; + delete ncells; + delete mult; + + +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV1::RunAlgoritm(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* + 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 <= fHeader->GetLegoNbinEta(); i++) { + for (Int_t j = 1; j <= fHeader->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 = fHeader->GetMinMove(); + Float_t maxmove = fHeader->GetMaxMove(); + Float_t rc = fHeader->GetRadius(); + Float_t etseed = fHeader->GetEtSeed(); + //Float_t etmin = fHeader->GetMinJetEt(); + + + + // 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; + + + 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; + } + } + }//end of cells loop looking centroide + + //avoid cones overloap (to be implemented in the future) + + //flag cells in Rc, estimate total energy in cone + Float_t etCone = 0.0; + Int_t nCellIn = 0; + rc = fHeader->GetRadius(); + 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++; + } + } + + // 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 !!! + // reject tmp outside acceptable eta range + if ((eta < (fHeader->GetJetEtaMax())) && + (eta > (fHeader->GetJetEtaMin())) && + (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]]; + } + + + //delete + delete index; + delete idx; + +} +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV1::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, + Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, + Int_t* multJet, Int_t* injet) +{ + //background subtraction using cone method but without correction in dE/deta distribution + + //calculate energy inside and outside cones + Float_t rc= fHeader->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 + etIn[ijet] += ptT[jpart]; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart]; + multJet[ijet]++; + injet[jpart] = ijet; + break; + } + }// end jets loop + if(injet[jpart] == 0) etOut += ptT[jpart]; // particle outside cones + } //end particle loop + + //estimate jets and background areas + Float_t areaJet[30]; + Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi(); + for(Int_t k=0; k fHeader->GetLegoEtaMax()){ // sector outside etamax + Float_t h = fHeader->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin + Float_t h = fHeader->GetLegoEtaMax() + etaJet[k]; + 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 AliUA1JetFinderV1::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN, + Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, + Int_t* multJet, Int_t* injet) +{ + + //background subtraction using statistical method + + Float_t etbgStat = fHeader->GetBackgStat(); // pre-calculated background + + //calculate energy inside + Float_t rc= fHeader->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 + etIn[ijet]+= ptT[jpart]; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + multJet[ijet]++; + injet[jpart] = ijet; + break; + } + }// end jets loop + } //end particle loop + + //calc jets areas + Float_t areaJet[30]; + Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi(); + for(Int_t k=0; k fHeader->GetLegoEtaMax()){ // sector outside etamax + Float_t h = fHeader->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin + Float_t h = fHeader->GetLegoEtaMax() + etaJet[k]; + 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 = fHeader->GetLegoEtaMax(); + Float_t etamin = fHeader->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; // 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 + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + multJet[ijet]++; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + injet[jpart] = ijet; + break; + } + }// end jets loop + if(injet[jpart] == 0) 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*(fHeader->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetGetBackgCutRatio(); + + + //general + Float_t rc= fHeader->GetRadius(); + Float_t etamax = fHeader->GetLegoEtaMax(); + Float_t etamin = fHeader->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 + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + multJet[ijet]++; + injet[jpart] = ijet; + break; + } + }// end jets loop + if(injet[jpart] == 0) 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*(fHeader->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetReset(); + fJets->ClearJets(); +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV1::WriteJHeaderToFile() +{ + fOut->cd(); + fHeader->Write(); +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV1::Init() +{ + // initializes some variables + + // book lego + fLego = new + TH2F("legoH","eta-phi", + fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), + fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(), + fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax()); + +} diff --git a/JETAN/AliUA1JetFinderV1.h b/JETAN/AliUA1JetFinderV1.h new file mode 100644 index 00000000000..a71abde7262 --- /dev/null +++ b/JETAN/AliUA1JetFinderV1.h @@ -0,0 +1,67 @@ +#ifndef ALIUA1JETFINDERV1_H +#define ALIUA1JETFINDERV1_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + +//--------------------------------------------------------------------- +// UA1 Cone Algorithm Finder V1 +// manages the search for jets +// Author: Rafael.Diaz.Valdes@cern.ch +// (version in c++) +//--------------------------------------------------------------------- + +#include "AliJetFinder.h" +class AliUA1JetHeaderV1; +class TH2F; + +class AliUA1JetFinderV1 : public AliJetFinder +{ + public: + + AliUA1JetFinderV1(); + ~AliUA1JetFinderV1(); + + // getters + + // setters + void SetJetHeader(AliUA1JetHeaderV1* h) {fHeader= h;} + // others + void FindJets(); + void RunAlgoritm(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 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 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 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 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 Reset(); + void Init(); + void WriteJHeaderToFile(); + + protected: + + AliUA1JetHeaderV1* fHeader; // pointer to jet header + TH2F * fLego; //! Lego Histo + + ClassDef(AliUA1JetFinderV1,1) +}; + +#endif diff --git a/JETAN/AliUA1JetHeaderV1.cxx b/JETAN/AliUA1JetHeaderV1.cxx new file mode 100644 index 00000000000..cd25614f2d5 --- /dev/null +++ b/JETAN/AliUA1JetHeaderV1.cxx @@ -0,0 +1,70 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//--------------------------------------------------------------------- +// Jet Gen header class +// Stores parameters of particle algoritm +// Author: Rafael.Diaz.Valdes@cern.ch +//--------------------------------------------------------------------- + +#include +#include "AliUA1JetHeaderV1.h" +ClassImp(AliUA1JetHeaderV1) + +//////////////////////////////////////////////////////////////////////// + +AliUA1JetHeaderV1::AliUA1JetHeaderV1(): + AliJetHeader("AliUA1JetHeaderV1") +{ + // Constructor + fConeRadius = 0.3; + fEtSeed = 3.0; + fMinJetEt = 10.0; + fMinMove = 0.05; + fMaxMove = 0.15; + fBackgMode = 1; // subtract backg + fPrecBg = 0.035; //background prec + fBackgStat = 0.0; // pre-calculated background used in statistic subtraction method + fBackgCutRatio = 1.0; // pre-calculated pt-cut ratio used in ratio subtraction method + fNAcceptJets = 3; // number of accepted jets per events + fLegoNbinEta = 36; + fLegoNbinPhi = 124; + fLegoPhiMin = 0.; + fLegoPhiMax = 2. * TMath::Pi(); + fLegoEtaMin = -0.9; + fLegoEtaMax = 0.9; +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetHeaderV1::PrintParameters() const +{ + // prints out parameters of jet algorithm + + cout << " UA version1 jet algorithm " << endl; + cout << " * Jet parameters: " << endl; + cout << " Cone size: " << fConeRadius<< endl; + cout << " Minimum energy for a seed: " << fEtSeed << endl; + cout << " Minumum energy for a jet: " << fMinJetEt << endl; + cout << " Minimum allowed move: " << fMinMove << endl; + cout << " Maximum allowed move: " << fMaxMove << endl; + cout << " * Lego parameters: " << endl; + cout << " Number of bins in eta: " << fLegoNbinEta<< endl; + cout << " Number of bins in phi: " << fLegoNbinPhi<< endl; + cout << " Minimum azimuthal angle: " << fLegoPhiMin<< endl; + cout << " Maximum azimuthal angle: " << fLegoPhiMax<< endl; + cout << " Minimum rapidity angle: " << fLegoEtaMin<< endl; + cout << " Maximum rapidity angle: " << fLegoEtaMax<< endl; +} diff --git a/JETAN/AliUA1JetHeaderV1.h b/JETAN/AliUA1JetHeaderV1.h new file mode 100644 index 00000000000..4ed7e98cad3 --- /dev/null +++ b/JETAN/AliUA1JetHeaderV1.h @@ -0,0 +1,89 @@ +#ifndef ALIUA1JETHEADERV1_H +#define ALIUA1JETHEADERV1_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Jet Finder header class for algoritm using particles +// Stores the parameters of the particles algoritm +// Author: Rafael.Diaz.Valdes@cern.ch +//--------------------------------------------------------------------- + +#include "AliJetHeader.h" + + +class AliUA1JetHeaderV1 : public AliJetHeader +{ + public: + + AliUA1JetHeaderV1(); + virtual ~AliUA1JetHeaderV1() { } + + // Getters + Float_t GetRadius() const {return fConeRadius;} + Float_t GetMinMove() const {return fMinMove;} + Float_t GetMaxMove() const {return fMaxMove;} + Float_t GetEtSeed() const {return fEtSeed;} + Float_t GetMinJetEt() const {return fMinJetEt;} + Int_t GetLegoNbinEta() const {return fLegoNbinEta;} + Int_t GetLegoNbinPhi() const {return fLegoNbinPhi;} + Float_t GetLegoEtaMin() const {return fLegoEtaMin;} + Float_t GetLegoEtaMax() const {return fLegoEtaMax;} + Float_t GetLegoPhiMin() const {return fLegoPhiMin;} + Float_t GetLegoPhiMax() const {return fLegoPhiMax;} + Int_t GetBackgMode() const {return fBackgMode;} + Float_t GetPrecBg() const {return fPrecBg;} + Float_t GetBackgStat() const {return fBackgStat;} + Float_t GetBackgCutRatio() const {return fBackgCutRatio;} + Int_t GetNAcceptJets() const {return fNAcceptJets;} + + + // Setters + void SetRadius(Float_t f) {fConeRadius=f;} + void SetMinMove(Float_t f) {fMinMove=f;} + void SetMaxMove(Float_t f) {fMaxMove=f;} + void SetEtSeed(Float_t f) {fEtSeed=f;} + void SetMinJetEt(Float_t f) {fMinJetEt=f;} + void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;} + void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;} + void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;} + void SetLegoEtaMax(Float_t f) {fLegoEtaMax=f;} + void SetLegoPhiMin(Float_t f) {fLegoPhiMin=f;} + void SetLegoPhiMax(Float_t f) {fLegoPhiMax=f;} + void BackgMode(Int_t mod ) {fBackgMode = mod;} + void SetPrecBg(Float_t f) {fPrecBg=f;} + void SetBackgStat(Float_t f) {fBackgStat=f;} + void SetBackgCutRatio(Float_t f) {fBackgCutRatio=f;} + void SetNAcceptJets(Int_t ajets) {fNAcceptJets = ajets;} + + // others + void PrintParameters() const; + +protected: + + // parameters of algorithm + Float_t fConeRadius; // Cone radius + Float_t fEtSeed; // Min. Et for seed + Float_t fMinJetEt; // Min Et of jet + // parameters of backgound substraction + Float_t fMinMove; // min cone move + Float_t fMaxMove; // max cone move + Int_t fBackgMode; // background subtraction mode + Float_t fPrecBg; // max value of change for BG (in %) + Float_t fBackgStat; // pre-calculated background used in statistic subtraction method + Float_t fBackgCutRatio; // pre-calculated pt-cut ratio used in ratio subtraction method + Int_t fNAcceptJets; // number of accepted jets per events + + // parameters for legos + Int_t fLegoNbinEta; //! number of cells in eta + Int_t fLegoNbinPhi; //! number of cells in phi + Float_t fLegoEtaMin; //! minimum eta + Float_t fLegoEtaMax; //! maximum eta + Float_t fLegoPhiMin; //! minimun phi + Float_t fLegoPhiMax; //! maximum phi + + ClassDef(AliUA1JetHeaderV1,1) +}; + +#endif diff --git a/JETAN/JetAnalysisLinkDef.h b/JETAN/JetAnalysisLinkDef.h index d8abac04ec5..ad26b6c8df3 100644 --- a/JETAN/JetAnalysisLinkDef.h +++ b/JETAN/JetAnalysisLinkDef.h @@ -23,4 +23,6 @@ #pragma link C++ class AliJetProductionDataPDC2004+; #pragma link C++ class AliJetAnalysis+; #pragma link C++ class AliJetDistributions+; +#pragma link C++ class AliUA1JetFinderV1+; +#pragma link C++ class AliUA1JetHeaderV1+; #endif diff --git a/JETAN/libJETAN.pkg b/JETAN/libJETAN.pkg index 6d0d1a0f8e8..befd8f3c457 100644 --- a/JETAN/libJETAN.pkg +++ b/JETAN/libJETAN.pkg @@ -7,7 +7,8 @@ SRCS = AliJet.cxx AliJetHeader.cxx AliPxconeJetHeader.cxx \ AliLeading.cxx \ AliJetKineReader.cxx AliJetKineReaderHeader.cxx \ AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \ - AliJetAnalysis.cxx AliJetDistributions.cxx + AliJetAnalysis.cxx AliJetDistributions.cxx \ + AliUA1JetFinderV1.cxx AliUA1JetHeaderV1.cxx FSRCS = pxcone.F ua1_jet_finder.F -- 2.43.0