JETAN/AliUA1JetFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //---------------------------------------------------------------------
  17 // UA1 Jet finder
  18 // manages the search for jets
  19 // Author: jgcn@mda.cinvestav.mx
  20 // (code adapted from EMCAL directory)
  21 //---------------------------------------------------------------------
  22
  23 #include <Riostream.h>
  24 #include <TLorentzVector.h>
  25 #include <TFile.h>
  26 #include <TH2F.h>
  27 #include "AliUA1JetFinder.h"
  28 #include "AliUA1JetHeader.h"
  29 #include "UA1Common.h"
  30 #include "AliJetReaderHeader.h"
  31 #include "AliJetReader.h"
  32 #include "AliJet.h"
  33
  34
  35 ClassImp(AliUA1JetFinder)
  36
  37 ////////////////////////////////////////////////////////////////////////
  38
  39 AliUA1JetFinder::AliUA1JetFinder()
  40
  41 {
  42   // Constructor
  43   fHeader = 0x0;
  44   fLego   = 0x0;
  45 }
  46
  47 ////////////////////////////////////////////////////////////////////////
  48
  49 AliUA1JetFinder::~AliUA1JetFinder()
  50
  51 {
  52   // destructor
  53 }
  54
  55 ////////////////////////////////////////////////////////////////////////
  56
  57 #ifndef WIN32
  58 # define ua1_jet_finder ua1_jet_finder_
  59 # define hf1 hf1_
  60 # define type_of_call
  61
  62 #else
  63 # define ua1_jet_finder UA1_JET_FINDER
  64 # define hf1 HF1
  65 # define type_of_call _stdcall
  66 #endif
  67
  68 extern "C" void type_of_call
  69 ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
  70                Float_t  etc[60000],  Float_t etac[60000],
  71                Float_t  phic[60000],
  72                Float_t& min_move, Float_t& max_move, Int_t& mode,
  73                Float_t& prec_bg,  Int_t& ierror);
  74
  75 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
  76
  77 ////////////////////////////////////////////////////////////////////////
  78
  79 void AliUA1JetFinder::FindJets()
  80
  81 {
  82   // initialize event, look for jets, download jet info
  83
  84   // initialization in 2 steps
  85   // 1) transform input to pt,eta,phi plus lego
  86   // 2) dump lego
  87
  88   // 1) transform input to pt,eta,phi plus lego
  89   TClonesArray *lvArray = fReader->GetMomentumArray();
  90   Int_t nIn =  lvArray->GetEntries();
  91   if (nIn == 0) return;
  92
  93   // local arrays for input
  94   Float_t* enT  = new Float_t[nIn];
  95   Float_t* ptT  = new Float_t[nIn];
  96   Float_t* etaT = new Float_t[nIn];
  97   Float_t* phiT = new Float_t[nIn];
  98
  99   // load input vectors
 100   for (Int_t i = 0; i < nIn; i++){
 101     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 102     enT[i]  = lv->Energy();
 103     ptT[i]  = lv->Pt();
 104     etaT[i] = lv->Eta();
 105     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
 106     if (fReader->GetCutFlag(i) == 1)
 107       fLego->Fill(etaT[i], phiT[i], ptT[i]);
 108   }
 109   fJets->SetNinput(nIn);
 110
 111   // 2) dump lego
 112   // check enough space! *to be done*
 113   Float_t etCell[60000];   //! Cell Energy
 114   Float_t etaCell[60000];  //! Cell eta
 115   Float_t phiCell[60000];  //! Cell phi
 116
 117   Int_t nCell = 0;
 118   TAxis* xaxis = fLego->GetXaxis();
 119   TAxis* yaxis = fLego->GetYaxis();
 120   Float_t e = 0.0;
 121   for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
 122       for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
 123           e = fLego->GetBinContent(i,j);
 124           if (e > 0.0) e -= fHeader->GetMinCellEt();
 125           if (e < 0.0) e = 0.;
 126           Float_t eta  = xaxis->GetBinCenter(i);
 127           Float_t phi  = yaxis->GetBinCenter(j);
 128           etCell[nCell]  = e;
 129           etaCell[nCell] = eta;
 130           phiCell[nCell] = phi;
 131           nCell++;
 132       }
 133   }
 134
 135   // run the algo. Parameters from header
 136   Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
 137   Float_t minmove = fHeader->GetMinMove();
 138   Float_t maxmove = fHeader->GetMaxMove();
 139   Int_t mode      = fHeader->GetMode();
 140   Float_t precbg  = fHeader->GetPrecBg();
 141   Int_t ierror;
 142   ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
 143                  minmove, maxmove, mode, precbg, ierror);
 144
 145   // sort jets
 146   Int_t * idx  = new Int_t[UA1JETS.njet];
 147   TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
 148
 149   // download jet info.
 150   for(Int_t i = 0; i < UA1JETS.njet; i++) {
 151     // reject events outside acceptable eta range
 152     if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
 153         || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
 154           continue;
 155       Float_t px, py,pz,en; // convert to 4-vector
 156       px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
 157       py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
 158       pz = UA1JETS.etj[idx[i]] /
 159           TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
 160       en = TMath::Sqrt(px * px + py * py + pz * pz);
 161       fJets->AddJet(px, py, pz, en);
 162   }
 163
 164   // find multiplicities and relationship jet-particle
 165   // find also percentage of pt from pythia
 166   Int_t* injet = new Int_t[nIn];
 167   Int_t* sflag = new Int_t[nIn];
 168   for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
 169   Int_t* mult  = new Int_t[fJets->GetNJets()];
 170   Int_t* ncell  = new Int_t[fJets->GetNJets()];
 171   Float_t* percentage  = new Float_t[fJets->GetNJets()];
 172
 173   for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
 174       Float_t pt_sig = 0.0;
 175       mult[i] = 0;
 176       ncell[i] = UA1JETS.ncellj[i];
 177       for (Int_t j = 0; j < nIn; j++) {
 178           Float_t deta = etaT[j] - fJets->GetEta(i);
 179           Float_t dphi = phiT[j] - fJets->GetPhi(i);
 180           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 181           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 182           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 183           if (dr < fHeader->GetRadius() && injet[j] == 0) {
 184               injet[j] = -(i+1);
 185               if ((fReader->GetCutFlag(j)) == 1 &&
 186                   (etaT[j] < fHeader->GetLegoEtaMax()) &&
 187                   (etaT[j] > fHeader->GetLegoEtaMin())) {
 188                   injet[j] = i+1;
 189                   mult[i]++;
 190                   if (fReader->GetSignalFlag(j) == 1) {
 191                     pt_sig+=ptT[j];
 192                     sflag[j]=1;
 193                   }
 194               }
 195           }
 196       }
 197       percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
 198         ((Double_t) fJets->GetPt(i));
 199   }
 200   fJets->SetNCells(ncell);
 201   fJets->SetPtFromSignal(percentage);
 202   fJets->SetMultiplicities(mult);
 203   fJets->SetInJet(injet);
 204   fJets->SetEtaIn(etaT);
 205   fJets->SetPhiIn(phiT);
 206   fJets->SetPtIn(ptT);
 207   fJets->SetEtAvg(UA1JETS.etavg);
 208   delete ncell;
 209   delete enT;
 210   delete ptT;
 211   delete etaT;
 212   delete phiT;
 213   delete injet;
 214   delete idx;
 215   delete mult;
 216   delete percentage;
 217 }
 218
 219 ////////////////////////////////////////////////////////////////////////
 220
 221 void AliUA1JetFinder::Reset()
 222 {
 223   fLego->Reset();
 224   fJets->ClearJets();
 225 }
 226
 227 ////////////////////////////////////////////////////////////////////////
 228
 229 void AliUA1JetFinder::WriteJHeaderToFile()
 230 {
 231   fOut->cd();
 232   fHeader->Write();
 233 }
 234
 235 ////////////////////////////////////////////////////////////////////////
 236
 237 void AliUA1JetFinder::Init()
 238 {
 239   // initializes some variables
 240   Float_t dEta, dPhi;
 241   dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
 242     /((Float_t) fHeader->GetLegoNbinEta());
 243   dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
 244     /((Float_t) fHeader->GetLegoNbinPhi());
 245
 246   UA1CELL.etaCellSize = dEta;
 247   UA1CELL.phiCellSize = dPhi;
 248
 249   // jet parameters
 250   UA1PARA.coneRad = fHeader->GetRadius();
 251   UA1PARA.etSeed  = fHeader->GetEtSeed();
 252   UA1PARA.ejMin   = fHeader->GetMinJetEt();
 253   UA1PARA.etMin   = fHeader->GetMinCellEt();
 254
 255   // book lego
 256   fLego = new
 257     TH2F("legoH","eta-phi",
 258          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
 259          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
 260          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
 261 }