/************************************************************************** * 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 Jet finder // manages the search for jets // Author: jgcn@mda.cinvestav.mx // (code adapted from EMCAL directory) //--------------------------------------------------------------------- #include #include #include #include #include "AliUA1JetFinder.h" #include "AliUA1JetHeader.h" #include "UA1Common.h" #include "AliJetReaderHeader.h" #include "AliJetReader.h" #include "AliJet.h" ClassImp(AliUA1JetFinder) //////////////////////////////////////////////////////////////////////// AliUA1JetFinder::AliUA1JetFinder() { // Constructor fHeader = 0x0; fLego = 0x0; } //////////////////////////////////////////////////////////////////////// AliUA1JetFinder::~AliUA1JetFinder() { // destructor } //////////////////////////////////////////////////////////////////////// #ifndef WIN32 # define ua1_jet_finder ua1_jet_finder_ # define hf1 hf1_ # define type_of_call #else # define ua1_jet_finder UA1_JET_FINDER # define hf1 HF1 # define type_of_call _stdcall #endif extern "C" void type_of_call ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot, Float_t etc[60000], Float_t etac[60000], Float_t phic[60000], Float_t& min_move, Float_t& max_move, Int_t& mode, Float_t& prec_bg, Int_t& ierror); extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt); //////////////////////////////////////////////////////////////////////// void AliUA1JetFinder::FindJets() { // initialize event, look for jets, download jet info // initialization in 2 steps // 1) transform input to pt,eta,phi plus lego // 2) dump lego // 1) 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* enT = new Float_t[nIn]; Float_t* ptT = new Float_t[nIn]; Float_t* etaT = new Float_t[nIn]; Float_t* phiT = new Float_t[nIn]; // load input vectors for (Int_t i = 0; i < nIn; i++){ TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); enT[i] = lv->Energy(); 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) fLego->Fill(etaT[i], phiT[i], ptT[i]); } fJets->SetNinput(nIn); // 2) 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 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) e -= fHeader->GetMinCellEt(); if (e < 0.0) e = 0.; Float_t eta = xaxis->GetBinCenter(i); Float_t phi = yaxis->GetBinCenter(j); etCell[nCell] = e; etaCell[nCell] = eta; phiCell[nCell] = phi; nCell++; } } // run the algo. Parameters from header Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi()); Float_t minmove = fHeader->GetMinMove(); Float_t maxmove = fHeader->GetMaxMove(); Int_t mode = fHeader->GetMode(); Float_t precbg = fHeader->GetPrecBg(); Int_t ierror; ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, minmove, maxmove, mode, precbg, ierror); // sort jets Int_t * idx = new Int_t[UA1JETS.njet]; TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx); // download jet info. for(Int_t i = 0; i < UA1JETS.njet; i++) { // reject events outside acceptable eta range if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax())) || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin()))) continue; Float_t px, py,pz,en; // convert to 4-vector px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]); py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]); pz = UA1JETS.etj[idx[i]] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]]))); en = TMath::Sqrt(px * px + py * py + pz * pz); fJets->AddJet(px, py, pz, en); } // find multiplicities and relationship jet-particle // find also percentage of pt from pythia Int_t* injet = new Int_t[nIn]; Int_t* sflag = new Int_t[nIn]; for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;} Int_t* mult = new Int_t[fJets->GetNJets()]; Int_t* ncell = new Int_t[fJets->GetNJets()]; Float_t* percentage = new Float_t[fJets->GetNJets()]; for(Int_t i = 0; i < (fJets->GetNJets()); i++) { Float_t pt_sig = 0.0; mult[i] = 0; ncell[i] = UA1JETS.ncellj[i]; for (Int_t j = 0; j < nIn; j++) { Float_t deta = etaT[j] - fJets->GetEta(i); Float_t dphi = phiT[j] - fJets->GetPhi(i); 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 < fHeader->GetRadius() && injet[j] == 0) { injet[j] = -(i+1); if ((fReader->GetCutFlag(j)) == 1 && (etaT[j] < fHeader->GetLegoEtaMax()) && (etaT[j] > fHeader->GetLegoEtaMin())) { injet[j] = i+1; mult[i]++; if (fReader->GetSignalFlag(j) == 1) { pt_sig+=ptT[j]; sflag[j]=1; } } } } percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/ ((Double_t) fJets->GetPt(i)); } fJets->SetNCells(ncell); fJets->SetPtFromSignal(percentage); fJets->SetMultiplicities(mult); fJets->SetInJet(injet); fJets->SetEtaIn(etaT); fJets->SetPhiIn(phiT); fJets->SetPtIn(ptT); fJets->SetEtAvg(UA1JETS.etavg); delete ncell; delete enT; delete ptT; delete etaT; delete phiT; delete injet; delete idx; delete mult; delete percentage; } //////////////////////////////////////////////////////////////////////// void AliUA1JetFinder::Reset() { fLego->Reset(); fJets->ClearJets(); } //////////////////////////////////////////////////////////////////////// void AliUA1JetFinder::WriteJHeaderToFile() { fOut->cd(); fHeader->Write(); } //////////////////////////////////////////////////////////////////////// void AliUA1JetFinder::Init() { // initializes some variables Float_t dEta, dPhi; dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin()) /((Float_t) fHeader->GetLegoNbinEta()); dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin()) /((Float_t) fHeader->GetLegoNbinPhi()); UA1CELL.etaCellSize = dEta; UA1CELL.phiCellSize = dPhi; // jet parameters UA1PARA.coneRad = fHeader->GetRadius(); UA1PARA.etSeed = fHeader->GetEtSeed(); UA1PARA.ejMin = fHeader->GetMinJetEt(); UA1PARA.etMin = fHeader->GetMinCellEt(); // book lego fLego = new TH2F("legoH","eta-phi", fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(), fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax()); }