| 45a58699 |
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 | /* $Id$ */ |
| 17 | |
| 18 | //*-- Authors: Andreas Morsch (CERN) |
| 19 | //* J.L. Klay (LBL) |
| 20 | //* Aleksei Pavlinov (WSU) |
| 21 | //-- |
| 22 | //-- |
| 23 | // |
| 24 | |
| 25 | #include <stdio.h> |
| 26 | // From root ... |
| 27 | #include <TArrayF.h> |
| 28 | #include <TAxis.h> |
| 29 | #include <TBranchElement.h> |
| 30 | #include <TCanvas.h> |
| 31 | #include <TClonesArray.h> |
| 32 | #include <TFile.h> |
| 33 | #include <TH1.h> |
| 34 | #include <TH2.h> |
| 35 | #include <TList.h> |
| 36 | #include <TPDGCode.h> |
| 37 | #include <TPad.h> |
| 38 | #include <TParticle.h> |
| 39 | #include <TParticlePDG.h> |
| 40 | #include <TPaveText.h> |
| 41 | #include <TPythia6Calls.h> |
| 42 | #include <TROOT.h> |
| 43 | #include <TStyle.h> |
| 44 | #include <TTree.h> |
| 45 | #include <TBrowser.h> |
| 46 | |
| 47 | // From AliRoot ... |
| 48 | #include "AliEMCALJetFinder.h" |
| 49 | #include "AliHeader.h" |
| 50 | #include "AliMagF.h" |
| 45a58699 |
51 | #include "AliRun.h" |
| 52 | #include "AliGenerator.h" |
| 53 | #include "AliRunLoader.h" |
| 54 | #include "AliEMCAL.h" |
| 55 | #include "AliEMCALLoader.h" |
| 56 | #include "AliEMCALDigit.h" |
| 57 | #include "AliEMCALDigitizer.h" |
| 58 | #include "AliEMCALFast.h" |
| 59 | #include "AliEMCALGeometry.h" |
| 60 | #include "AliEMCALHadronCorrection.h" |
| 61 | #include "AliEMCALHit.h" |
| 62 | #include "AliEMCALJetMicroDst.h" |
| 63 | |
| 64 | // Interface to FORTRAN |
| 65 | #include "Ecommon.h" |
| 66 | #include "AliMC.h" |
| 67 | |
| 68 | |
| 69 | ClassImp(AliEMCALJetFinder) |
| 70 | |
| 71 | //____________________________________________________________________________ |
| 72 | AliEMCALJetFinder::AliEMCALJetFinder() |
| 18a21c7c |
73 | : fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE), |
| 74 | fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0), |
| 75 | fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0), |
| 76 | fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0), |
| 77 | fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0), |
| 78 | fDebug(0), fBackground(0), fConeRadius(0.3), |
| 79 | fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.), |
| 80 | fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0), |
| 81 | fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0), |
| 82 | fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.), |
| 83 | fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0), |
| 84 | fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0), |
| 85 | fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.), |
| 86 | fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0) |
| 45a58699 |
87 | { |
| 88 | // Default constructor |
| 89 | fJets = 0; |
| 90 | fNjets = 0; |
| 91 | fLego = 0; |
| 92 | fLegoB = 0; |
| 93 | |
| 94 | fTrackList = 0; |
| 95 | fPtT = 0; |
| 96 | fEtaT = 0; |
| 97 | fPhiT = 0; |
| 98 | fPdgT = 0; |
| 99 | |
| 100 | fTrackListB = 0; |
| 101 | fPtB = 0; |
| 102 | fEtaB = 0; |
| 103 | fPhiB = 0; |
| 104 | fPdgB = 0; |
| 105 | |
| 106 | fHCorrection = 0; |
| 107 | fHadronCorrector = 0; |
| 108 | |
| 109 | fWrite = 0; |
| 110 | fOutFileName = 0; |
| 111 | fOutFile = 0; |
| 112 | fInFile = 0; |
| 113 | fEvent = 0; |
| 114 | |
| 115 | fRandomBg = 0; |
| 116 | |
| 117 | SetParametersForBgSubtraction(); |
| 118 | } |
| 119 | |
| 120 | AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title) |
| 18a21c7c |
121 | : TTask(name, title), |
| 122 | fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE), |
| 123 | fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0), |
| 124 | fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0), |
| 125 | fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0), |
| 126 | fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0), |
| 127 | fDebug(0), fBackground(0), fConeRadius(0.3), |
| 128 | fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.), |
| 129 | fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0), |
| 130 | fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0), |
| 131 | fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.), |
| 132 | fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0), |
| 133 | fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0), |
| 134 | fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.), |
| 135 | fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0) |
| 45a58699 |
136 | { |
| 137 | // Constructor |
| 138 | // Title is used in method GetFileNameForParameters(); |
| 139 | // |
| 140 | fJets = new TClonesArray("AliEMCALJet",10000); |
| 141 | fNjets = 0; |
| 142 | for (Int_t i = 0; i < 30000; i++) |
| 143 | { |
| 144 | fEtCell[i] = 0.; |
| 145 | fEtaCell[i] = 0.; |
| 146 | fPhiCell[i] = 0.; |
| 147 | } |
| 148 | fLego = 0; |
| 149 | fLegoB = 0; |
| 150 | |
| 151 | fTrackList = 0; |
| 152 | fPtT = 0; |
| 153 | fEtaT = 0; |
| 154 | fPhiT = 0; |
| 155 | fPdgT = 0; |
| 156 | |
| 157 | fTrackListB = 0; |
| 158 | fPtB = 0; |
| 159 | fEtaB = 0; |
| 160 | fPhiB = 0; |
| 161 | fPdgB = 0; |
| 162 | |
| 163 | fHCorrection = 0; |
| 164 | fHadronCorrector = 0; |
| 165 | fBackground = 0; |
| 166 | fWrite = 0; |
| 167 | fOutFileName = 0; |
| 168 | fOutFile = 0; |
| 169 | fInFile = 0; |
| 170 | fEvent = 0; |
| 171 | // |
| 172 | SetPtCut(); |
| 173 | SetMomentumSmearing(); |
| 174 | SetEfficiencySim(); |
| 175 | SetDebug(); |
| 176 | SetHadronCorrection(); |
| 177 | SetIncludeK0andN(); |
| 178 | |
| 179 | fRandomBg = 0; |
| 180 | SetParametersForBgSubtraction(); |
| 181 | } |
| 182 | |
| 18a21c7c |
183 | |
| 184 | AliEMCALJetFinder::AliEMCALJetFinder(const AliEMCALJetFinder& jf) |
| 185 | : TTask(jf.GetName(), jf.GetTitle()), |
| 186 | fBGFileName(jf.fBGFileName), fEMCALWeight(jf.fEMCALWeight), fTrackWeight(jf.fTrackWeight), |
| 187 | fRandomBg(jf.fRandomBg),fWrite(jf.fWrite), fWeightingMethod(jf.fWeightingMethod), fJets(jf.fJets), |
| 188 | fLego(jf.fLego), fLegoB(jf.fLegoB), fhLegoTracks(jf.fhLegoTracks), fhLegoEMCAL(jf.fhLegoEMCAL), |
| 189 | fhLegoHadrCorr(jf.fhLegoHadrCorr), fhEff(jf.fhEff), fhCellEt(jf.fhCellEt),fhCellEMCALEt(jf.fhCellEMCALEt), |
| 190 | fhTrackPt(jf.fhTrackPt), fhTrackPtBcut(jf.fhTrackPtBcut), fhChPartMultInTpc(jf.fhChPartMultInTpc), |
| 191 | fhSinTheta(jf.fhSinTheta), fC1(jf.fC1), fHistsList(jf.fHistsList), fHadronCorrector(jf.fHadronCorrector), |
| 192 | fHCorrection(jf.fHCorrection),fDebug(jf.fDebug), fBackground(jf.fBackground), fConeRadius(jf.fConeRadius), |
| 193 | fPtCut(jf.fPtCut), fEtSeed(jf.fEtSeed), fMinJetEt(jf.fMinJetEt), fMinCellEt(jf.fMinCellEt), |
| 194 | fSamplingF(jf.fSamplingF), fSmear(jf.fSmear), fEffic(jf.fEffic), fK0N(jf.fK0N), fNjets(jf.fNjets), |
| 195 | fDeta(jf.fDeta), fDphi(jf.fDphi), fNcell(jf.fNcell), fNtot(jf.fNtot), fNbinEta(jf.fNbinEta), |
| 196 | fNbinPhi(jf.fNbinPhi), fEtaMin(jf.fEtaMin), fEtaMax(jf.fEtaMax), fPhiMin(jf.fPhiMin), |
| 197 | fPhiMax(jf.fPhiMax), fNt(jf.fNt), fNChTpc(jf.fNChTpc), fNtS(jf.fNtS), fTrackList(jf.fTrackList), |
| 198 | fPtT(jf.fPtT), fEtaT(jf.fEtaT), fPhiT(jf.fPhiT), fPdgT(jf.fPdgT), fNtB(jf.fNtB), fTrackListB(jf.fTrackListB), |
| 199 | fPtB(jf.fPtB), fEtaB(jf.fEtaB),fPhiB(jf.fPhiB), fPdgB(jf.fPdgB), fMode(jf.fMode), fMinMove(jf.fMinMove), |
| 200 | fMaxMove(jf.fMaxMove), fPrecBg(jf.fPrecBg), fError(jf.fError), fOutFileName(jf.fOutFileName), |
| 201 | fOutFile(jf.fOutFile), fInFile(jf.fInFile), fEvent(jf.fEvent) |
| 202 | { |
| 203 | // Copy Constructor |
| 204 | } |
| 205 | |
| 45a58699 |
206 | void AliEMCALJetFinder::SetParametersForBgSubtraction |
| 207 | (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg) |
| 208 | { |
| 209 | // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc |
| 210 | // at WSU Linux cluster - 11-feb-2002 |
| 211 | // These parameters must be tuned more carefull !!! |
| 212 | SetMode(mode); |
| 213 | SetMinMove(minMove); |
| 214 | SetMaxMove(maxMove); |
| 215 | SetPrecBg(precBg); |
| 216 | } |
| 217 | |
| 218 | //____________________________________________________________________________ |
| 219 | AliEMCALJetFinder::~AliEMCALJetFinder() |
| 220 | { |
| 221 | // Destructor |
| 222 | // |
| 223 | if (fJets){ |
| 224 | fJets->Delete(); |
| 225 | delete fJets; |
| 226 | } |
| 227 | delete fLego; |
| 228 | delete fLegoB; |
| 229 | delete fhLegoTracks; |
| 230 | delete fhLegoEMCAL; |
| 231 | delete fhLegoHadrCorr; |
| 232 | delete fhEff; |
| 233 | delete fhCellEt; |
| 234 | delete fhCellEMCALEt; |
| 235 | delete fhTrackPt; |
| 236 | delete fhTrackPtBcut; |
| 237 | delete fhChPartMultInTpc; |
| 238 | |
| 239 | delete[] fTrackList; |
| 240 | delete[] fPtT; |
| 241 | delete[] fEtaT; |
| 242 | delete[] fPhiT; |
| 243 | delete[] fPdgT; |
| 244 | |
| 245 | delete[] fTrackListB; |
| 246 | delete[] fPtB; |
| 247 | delete[] fEtaB; |
| 248 | delete[] fPhiB; |
| 249 | delete[] fPdgB; |
| 250 | } |
| 251 | |
| 252 | #ifndef WIN32 |
| 253 | # define jet_finder_ua1 jet_finder_ua1_ |
| 254 | # define hf1 hf1_ |
| 255 | # define type_of_call |
| 256 | |
| 257 | #else |
| 258 | # define jet_finder_ua1 JET_FINDER_UA1 |
| 259 | # define hf1 HF1 |
| 260 | # define type_of_call _stdcall |
| 261 | #endif |
| 262 | |
| 263 | extern "C" void type_of_call |
| 264 | jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot, |
| 265 | Float_t etc[30000], Float_t etac[30000], |
| 266 | Float_t phic[30000], |
| 267 | Float_t& min_move, Float_t& max_move, Int_t& mode, |
| 268 | Float_t& prec_bg, Int_t& ierror); |
| 269 | |
| 270 | extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt); |
| 271 | |
| 272 | |
| 273 | void AliEMCALJetFinder::Init() |
| 274 | { |
| 275 | // |
| 276 | // Geometry and I/O initialization |
| 277 | // |
| 278 | // |
| 279 | // |
| 280 | // Get geometry parameters from EMCAL |
| 281 | // |
| 282 | // |
| 283 | // Geometry |
| 284 | //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); |
| 285 | AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(); |
| 286 | |
| 287 | // SetSamplingFraction(geom->GetSampling()); |
| 288 | |
| 289 | fNbinEta = geom->GetNZ(); |
| 290 | fNbinPhi = geom->GetNPhi(); |
| 291 | fPhiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.; |
| 292 | fPhiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.; |
| 293 | fEtaMin = geom->GetArm1EtaMin(); |
| 294 | fEtaMax = geom->GetArm1EtaMax(); |
| 295 | fDphi = (fPhiMax-fPhiMin)/fNbinPhi; |
| 296 | fDeta = (fEtaMax-fEtaMin)/fNbinEta; |
| 297 | fNtot = fNbinPhi*fNbinEta; |
| 298 | fWeightingMethod = kFALSE; |
| 299 | |
| 300 | // |
| 301 | SetCellSize(fDeta, fDphi); |
| 302 | // |
| 303 | // I/O |
| 304 | if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate"); |
| 305 | // |
| 306 | // |
| 307 | } |
| 308 | |
| 309 | void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncelltot, Float_t etc[30000], |
| 310 | Float_t etac[30000], Float_t phic[30000], |
| 311 | Float_t minmove, Float_t maxmove, Int_t mode, |
| 312 | Float_t precbg, Int_t ierror) |
| 313 | { |
| 314 | // Wrapper for fortran coded jet finder |
| 315 | // Full argument list |
| 316 | jet_finder_ua1(ncell, ncelltot, etc, etac, phic, |
| 317 | minmove, maxmove, mode, precbg, ierror); |
| 318 | // Write result to output |
| 319 | if(fWrite) WriteJets(); |
| 320 | fEvent++; |
| 321 | } |
| 322 | |
| 323 | void AliEMCALJetFinder::Find() |
| 324 | { |
| 325 | // Wrapper for fortran coded jet finder using member data for |
| 326 | // argument list |
| 327 | |
| 328 | Float_t minmove = fMinMove; |
| 329 | Float_t maxmove = fMaxMove; |
| 330 | Int_t mode = fMode; |
| 331 | Float_t precbg = fPrecBg; |
| 332 | Int_t ierror; |
| 333 | ResetJets(); // 4-feb-2002 by PAI |
| 334 | |
| 335 | jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, |
| 336 | minmove, maxmove, mode, precbg, ierror); |
| 337 | fError = ierror; |
| 338 | // Write result to output |
| 339 | Int_t njet = Njets(); |
| 340 | |
| 341 | for (Int_t nj=0; nj<njet; nj++) |
| 342 | { |
| 343 | if (fWeightingMethod) |
| 344 | { |
| 345 | fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ), |
| 346 | JetPhiW(nj), |
| 347 | JetEtaW(nj)); |
| 348 | |
| 349 | }else{ |
| 350 | |
| 351 | fJetT[nj] = new AliEMCALJet(JetEnergy(nj), |
| 352 | JetPhiW(nj), |
| 353 | JetEtaW(nj)); |
| 354 | } |
| 355 | fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod); |
| 356 | fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) ); |
| 357 | fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj) )); |
| 358 | |
| 359 | } |
| 360 | |
| 361 | FindTracksInJetCone(); |
| 362 | if(fWrite) WriteJets(); |
| 363 | fEvent++; |
| 364 | } |
| 365 | |
| 366 | |
| 367 | Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi) |
| 368 | { |
| 369 | // Calculate the energy in the cone |
| 370 | Float_t newenergy = 0.0; |
| 371 | Float_t bineta,binphi; |
| 372 | TAxis *x = fhLegoEMCAL->GetXaxis(); |
| 373 | TAxis *y = fhLegoEMCAL->GetYaxis(); |
| 374 | for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord |
| 375 | { |
| 376 | for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord |
| 377 | { |
| 378 | bineta = x->GetBinCenter(i); |
| 379 | binphi = y->GetBinCenter(j); |
| 380 | if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) |
| 381 | { |
| 382 | newenergy += fhLegoEMCAL->GetBinContent(i,j); |
| 383 | } |
| 384 | } |
| 385 | } |
| 386 | return newenergy; |
| 387 | } |
| 388 | |
| 389 | Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi) |
| 390 | { |
| 391 | // Calculate the track energy in the cone |
| 392 | Float_t newenergy = 0.0; |
| 393 | Float_t bineta,binphi; |
| 394 | TAxis *x = fhLegoTracks->GetXaxis(); |
| 395 | TAxis *y = fhLegoTracks->GetYaxis(); |
| 396 | for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord |
| 397 | { |
| 398 | for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord |
| 399 | { |
| 400 | bineta = x->GetBinCenter(i); |
| 401 | binphi = y->GetBinCenter(j); |
| 402 | if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) |
| 403 | { |
| 404 | newenergy += fhLegoTracks->GetBinContent(i,j); |
| 405 | } |
| 406 | } |
| 407 | } |
| 408 | return newenergy; |
| 409 | } |
| 410 | |
| 411 | |
| 412 | Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi) |
| 413 | { |
| 414 | // Calculate the weighted jet energy |
| 415 | |
| 416 | Float_t newenergy = 0.0; |
| 417 | Float_t bineta,binphi; |
| 418 | TAxis *x = fhLegoEMCAL->GetXaxis(); |
| 419 | TAxis *y = fhLegoEMCAL->GetYaxis(); |
| 420 | |
| 421 | |
| 422 | for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord |
| 423 | { |
| 424 | for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord |
| 425 | { |
| 426 | bineta = x->GetBinCenter(i); |
| 427 | binphi = y->GetBinCenter(j); |
| 428 | if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) |
| 429 | { |
| 430 | newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j); |
| 431 | } |
| 432 | } |
| 433 | } |
| 434 | |
| 435 | return newenergy; |
| 436 | |
| 437 | } |
| 438 | |
| 439 | |
| 440 | Int_t AliEMCALJetFinder::Njets() const |
| 441 | { |
| 442 | // Get number of reconstructed jets |
| 443 | return EMCALJETS.njet; |
| 444 | } |
| 445 | |
| 446 | Float_t AliEMCALJetFinder::JetEnergy(Int_t i) const |
| 447 | { |
| 448 | // Get reconstructed jet energy |
| 449 | return EMCALJETS.etj[i]; |
| 450 | } |
| 451 | |
| 452 | Float_t AliEMCALJetFinder::JetPhiL(Int_t i) const |
| 453 | { |
| 454 | // Get reconstructed jet phi from leading particle |
| 455 | return EMCALJETS.phij[0][i]; |
| 456 | } |
| 457 | |
| 458 | Float_t AliEMCALJetFinder::JetPhiW(Int_t i) const |
| 459 | { |
| 460 | // Get reconstructed jet phi from weighting |
| 461 | return EMCALJETS.phij[1][i]; |
| 462 | } |
| 463 | |
| 464 | Float_t AliEMCALJetFinder::JetEtaL(Int_t i) const |
| 465 | { |
| 466 | // Get reconstructed jet eta from leading particles |
| 467 | return EMCALJETS.etaj[0][i]; |
| 468 | } |
| 469 | |
| 470 | |
| 471 | Float_t AliEMCALJetFinder::JetEtaW(Int_t i) const |
| 472 | { |
| 473 | // Get reconstructed jet eta from weighting |
| 474 | return EMCALJETS.etaj[1][i]; |
| 475 | } |
| 476 | |
| 477 | void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi) |
| 478 | { |
| 479 | // Set grid cell size |
| 480 | EMCALCELLGEO.etaCellSize = eta; |
| 481 | EMCALCELLGEO.phiCellSize = phi; |
| 482 | } |
| 483 | |
| 484 | void AliEMCALJetFinder::SetConeRadius(Float_t par) |
| 485 | { |
| 486 | // Set jet cone radius |
| 487 | EMCALJETPARAM.coneRad = par; |
| 488 | fConeRadius = par; |
| 489 | } |
| 490 | |
| 491 | void AliEMCALJetFinder::SetEtSeed(Float_t par) |
| 492 | { |
| 493 | // Set et cut for seeds |
| 494 | EMCALJETPARAM.etSeed = par; |
| 495 | fEtSeed = par; |
| 496 | } |
| 497 | |
| 498 | void AliEMCALJetFinder::SetMinJetEt(Float_t par) |
| 499 | { |
| 500 | // Set minimum jet et |
| 501 | EMCALJETPARAM.ejMin = par; |
| 502 | fMinJetEt = par; |
| 503 | } |
| 504 | |
| 505 | void AliEMCALJetFinder::SetMinCellEt(Float_t par) |
| 506 | { |
| 507 | // Set et cut per cell |
| 508 | EMCALJETPARAM.etMin = par; |
| 509 | fMinCellEt = par; |
| 510 | } |
| 511 | |
| 512 | void AliEMCALJetFinder::SetPtCut(Float_t par) |
| 513 | { |
| 514 | // Set pT cut on charged tracks |
| 515 | fPtCut = par; |
| 516 | } |
| 517 | |
| 518 | |
| 519 | void AliEMCALJetFinder::Test() |
| 520 | { |
| 521 | // |
| 522 | // Test the finder call |
| 523 | // |
| 524 | const Int_t knmax = 30000; |
| 525 | Int_t ncell = 10; |
| 526 | Int_t ncelltot = 100; |
| 527 | |
| 528 | Float_t etc[knmax]; |
| 529 | Float_t etac[knmax]; |
| 530 | Float_t phic[knmax]; |
| 531 | Float_t minmove = 0; |
| 532 | Float_t maxmove = 0; |
| 533 | Int_t mode = 0; |
| 534 | Float_t precbg = 0; |
| 535 | Int_t ierror = 0; |
| 536 | |
| 537 | |
| 538 | Find(ncell, ncelltot, etc, etac, phic, |
| 539 | minmove, maxmove, mode, precbg, ierror); |
| 540 | |
| 541 | } |
| 542 | |
| 543 | // |
| 544 | // I/O |
| 545 | // |
| 546 | |
| 547 | void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet) |
| 548 | { |
| 549 | // |
| 550 | // Add a jet |
| 551 | // |
| 552 | TClonesArray &lrawcl = *fJets; |
| 553 | new(lrawcl[fNjets++]) AliEMCALJet(jet); |
| 554 | } |
| 555 | |
| 556 | void AliEMCALJetFinder::ResetJets() |
| 557 | { |
| 558 | // |
| 559 | // Reset Jet List |
| 560 | // |
| 561 | fJets->Clear(); |
| 562 | fNjets = 0; |
| 563 | } |
| 564 | |
| 565 | void AliEMCALJetFinder::WriteJets() |
| 566 | { |
| 567 | // |
| 568 | // Add all jets to the list |
| 569 | // |
| 570 | const Int_t kBufferSize = 4000; |
| 571 | const char* file = 0; |
| 572 | |
| 573 | Int_t njet = Njets(); |
| 574 | |
| 575 | for (Int_t nj = 0; nj < njet; nj++) |
| 576 | { |
| 577 | AddJet(*fJetT[nj]); |
| 578 | delete fJetT[nj]; |
| 579 | } |
| 580 | |
| 581 | // I/O |
| 582 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
| 583 | AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
| 584 | |
| 585 | if (!fOutFileName) { |
| 586 | // |
| 587 | // output written to input file |
| 588 | // |
| 589 | AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL"); |
| 590 | TTree* pK = rl->TreeK(); |
| 591 | file = (pK->GetCurrentFile())->GetName(); |
| 592 | TBranch * jetBranch ; |
| 593 | if (fDebug > 1) |
| 594 | printf("Make Branch - TreeR address %p %p\n",(void*)gAlice->TreeR(), (void*)pEMCAL); |
| 595 | //if (fJets && gAlice->TreeR()) { |
| 596 | if (fJets && emcalLoader->TreeR()) { |
| 597 | // pEMCAL->MakeBranchInTree(gAlice->TreeR(), |
| 598 | jetBranch = emcalLoader->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ; |
| 599 | //pEMCAL->MakeBranchInTree(gime->TreeR(), |
| 600 | // "EMCALJets", |
| 601 | // &fJets, |
| 602 | // kBufferSize, |
| 603 | // file); |
| 604 | |
| 605 | //Int_t nev = gAlice->GetHeader()->GetEvent(); |
| 606 | //gAlice->TreeR()->Fill(); |
| 607 | jetBranch->Fill(); |
| 608 | //char hname[30]; |
| 609 | //sprintf(hname,"TreeR%d", nev); |
| 610 | //gAlice->TreeR()->Write(hname); |
| 611 | //gAlice->TreeR()->Reset(); |
| 612 | //gime->WriteRecPoints("OVERWRITE"); |
| 613 | emcalLoader->WriteRecPoints("OVERWRITE"); |
| 614 | } |
| 615 | } else { |
| 616 | // |
| 617 | // Output written to user specified output file |
| 618 | // |
| 619 | //TTree* pK = gAlice->TreeK(); |
| 620 | TTree* pK = rl->TreeK(); |
| 621 | fInFile = pK->GetCurrentFile(); |
| 622 | |
| 623 | fOutFile->cd(); |
| 624 | char hname[30]; |
| 625 | sprintf(hname,"TreeR%d", fEvent); |
| 626 | TTree* treeJ = new TTree(hname, "EMCALJets"); |
| 627 | treeJ->Branch("EMCALJets", &fJets, kBufferSize); |
| 628 | treeJ->Fill(); |
| 629 | treeJ->Write(hname); |
| 630 | fInFile->cd(); |
| 631 | } |
| 632 | ResetJets(); |
| 633 | } |
| 634 | |
| 635 | void AliEMCALJetFinder::BookLego() |
| 636 | { |
| 637 | // |
| 638 | // Book histo for discretization |
| 639 | // |
| 640 | |
| 641 | // |
| 642 | // Don't add histos to the current directory |
| 643 | if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n"); |
| 644 | |
| 645 | // TH2::AddDirectory(0); // hists wil be put to the list from gROOT |
| 646 | // TH1::AddDirectory(0); |
| 647 | gROOT->cd(); |
| 648 | // |
| 649 | // Signal map |
| 650 | fLego = new TH2F("legoH","eta-phi", |
| 651 | fNbinEta, fEtaMin, fEtaMax, |
| 652 | fNbinPhi, fPhiMin, fPhiMax); |
| 653 | // |
| 654 | // Background map |
| 655 | fLegoB = new TH2F("legoB","eta-phi for BG event", |
| 656 | fNbinEta, fEtaMin, fEtaMax, |
| 657 | fNbinPhi, fPhiMin, fPhiMax); |
| 658 | |
| 659 | // Tracks map |
| 660 | fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks", |
| 661 | fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax); |
| 662 | // EMCAL map |
| 663 | fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL", |
| 664 | fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax); |
| 665 | // Hadron correction map |
| 666 | fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction", |
| 667 | fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax); |
| 668 | // Hists. for tuning jet finder |
| 669 | fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.); |
| 670 | |
| 671 | TArrayF eTmp(1101); |
| 672 | eTmp[0] = 0.0; |
| 673 | for(Int_t i=1; i<=1000; i++) eTmp[i] = 0.1*i; // step 100 mev |
| 674 | for(Int_t i=1001; i<=1100; i++) eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV |
| 675 | |
| 676 | fhCellEt = new TH1F("hCellEt","Cell E_{T} from fLego", |
| 677 | eTmp.GetSize()-1, eTmp.GetArray()); |
| 678 | fhCellEMCALEt = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself", |
| 679 | eTmp.GetSize()-1, eTmp.GetArray()); |
| 680 | fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ", |
| 681 | eTmp.GetSize()-1, eTmp.GetArray()); |
| 682 | fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut", |
| 683 | eTmp.GetSize()-1, eTmp.GetArray()); |
| 684 | |
| 685 | fhChPartMultInTpc = new TH1F("hChPartMultInTpc", |
| 686 | "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000); |
| 687 | |
| 688 | fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax); |
| 689 | TAxis *xax = fhSinTheta->GetXaxis(); |
| 690 | for(Int_t i=1; i<=fNbinEta; i++) { |
| 691 | Double_t eta = xax->GetBinCenter(i); |
| 692 | fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta) |
| 693 | } |
| 694 | |
| 695 | //! first canvas for drawing |
| 696 | fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE); |
| 697 | } |
| 698 | |
| 699 | void AliEMCALJetFinder::DumpLego() |
| 700 | { |
| 701 | // |
| 702 | // Dump lego histo into array |
| 703 | // |
| 704 | fNcell = 0; |
| 705 | TAxis* xaxis = fLego->GetXaxis(); |
| 706 | TAxis* yaxis = fLego->GetYaxis(); |
| 707 | // fhCellEt->Clear(); |
| 708 | Float_t e, eH; |
| 709 | for (Int_t i = 1; i <= fNbinEta; i++) { |
| 710 | for (Int_t j = 1; j <= fNbinPhi; j++) { |
| 711 | |
| 712 | e = fLego->GetBinContent(i,j); |
| 713 | if (fRandomBg) { |
| 714 | if (gRandom->Rndm() < 0.5) { |
| 715 | Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.)); |
| 716 | e += ebg; |
| 717 | } |
| 718 | } |
| 719 | if (e > 0.0) e -= fMinCellEt; |
| 720 | if (e < 0.0) e = 0.; |
| 721 | Float_t eta = xaxis->GetBinCenter(i); |
| 722 | Float_t phi = yaxis->GetBinCenter(j); |
| 723 | fEtCell[fNcell] = e; |
| 724 | fEtaCell[fNcell] = eta; |
| 725 | fPhiCell[fNcell] = phi; |
| 726 | fNcell++; |
| 727 | fhCellEt->Fill(e); |
| 728 | if(fhLegoEMCAL) { |
| 729 | eH = fhLegoEMCAL->GetBinContent(i,j); |
| 730 | if(eH > 0.0) fhCellEMCALEt->Fill(eH); |
| 731 | } |
| 732 | } // phi |
| 733 | } // eta |
| 734 | // today |
| 735 | // fNcell--; |
| 736 | } |
| 737 | |
| 738 | void AliEMCALJetFinder::ResetMap() |
| 739 | { |
| 740 | // |
| 741 | // Reset eta-phi array |
| 742 | |
| 743 | for (Int_t i=0; i<30000; i++) |
| 744 | { |
| 745 | fEtCell[i] = 0.; |
| 746 | fEtaCell[i] = 0.; |
| 747 | fPhiCell[i] = 0.; |
| 748 | } |
| 749 | } |
| 750 | |
| 751 | |
| 752 | void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) |
| 753 | { |
| 754 | // Which generator |
| 755 | // |
| 756 | const char* name = gAlice->Generator()->GetName(); |
| 757 | enum {kPythia, kHijing, kHijingPara}; |
| 758 | Int_t genType = 0; |
| 759 | |
| 760 | if (!strcmp(name, "Hijing")){ |
| 761 | genType = kHijing; |
| 762 | } else if (!strcmp(name, "Pythia")) { |
| 763 | genType = kPythia; |
| 764 | } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) { |
| 765 | genType = kHijingPara; |
| 766 | } |
| 767 | if (fDebug>=2) |
| 768 | printf("Fill tracks generated by %s %d\n", name, genType); |
| 769 | |
| 770 | |
| 771 | // |
| 772 | // Fill Cells with track information |
| 773 | // |
| 774 | if (fDebug >= 2) |
| 775 | printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich); |
| 776 | fNChTpc = 0; |
| 777 | |
| 778 | ResetMap(); |
| 779 | |
| 780 | if (!fLego) BookLego(); |
| 781 | // Reset |
| 782 | if (flag == 0) fLego->Reset(); |
| 783 | // |
| 784 | // Access particle information |
| 785 | Int_t npart = (gAlice->GetHeader())->GetNprimary(); |
| 786 | Int_t ntr = (gAlice->GetHeader())->GetNtrack(); |
| 787 | printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n", |
| 788 | npart, ntr, fLego->Integral()); |
| 789 | |
| 790 | // Create track list |
| 791 | // |
| 792 | // 0: not selected |
| 793 | // 1: selected for jet finding |
| 794 | // 2: inside jet |
| 795 | // .... |
| 796 | if (fTrackList) delete[] fTrackList; |
| 797 | if (fPtT) delete[] fPtT; |
| 798 | if (fEtaT) delete[] fEtaT; |
| 799 | if (fPhiT) delete[] fPhiT; |
| 800 | if (fPdgT) delete[] fPdgT; |
| 801 | |
| 802 | fTrackList = new Int_t [npart]; |
| 803 | fPtT = new Float_t[npart]; |
| 804 | fEtaT = new Float_t[npart]; |
| 805 | fPhiT = new Float_t[npart]; |
| 806 | fPdgT = new Int_t[npart]; |
| 807 | |
| 808 | fNt = npart; |
| 809 | fNtS = 0; |
| 810 | Float_t chTmp=0.0; // charge of current particle |
| 811 | // Int_t idGeant; |
| 812 | |
| 813 | // this is for Pythia ?? |
| 814 | for (Int_t part = 0; part < npart; part++) { |
| 815 | TParticle *mPart = gAlice->GetMCApp()->Particle(part); |
| 816 | Int_t mpart = mPart->GetPdgCode(); |
| 817 | Int_t child1 = mPart->GetFirstDaughter(); |
| 818 | Float_t pT = mPart->Pt(); |
| 819 | Float_t p = mPart->P(); |
| 820 | Float_t phi = mPart->Phi(); |
| 821 | Float_t eta = -100.; |
| 822 | if(pT > 0.001) eta = mPart->Eta(); |
| 823 | Float_t theta = mPart->Theta(); |
| 824 | if (fDebug>=2 && mPart->GetStatusCode()==1) { |
| 825 | printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n", |
| 826 | part, mpart, child1, mPart->GetLastDaughter(), mPart->GetStatusCode()); |
| 827 | } |
| 828 | |
| 829 | if (fDebug >= 2 && genType == kPythia) { |
| 830 | if (part == 6 || part == 7) |
| 831 | { |
| 832 | printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n", |
| 833 | part-5, pT, eta, phi); |
| 834 | } |
| 835 | } |
| 836 | |
| 837 | fTrackList[part] = 0; |
| 838 | fPtT[part] = pT; // must be change after correction for resolution !!! |
| 839 | fEtaT[part] = eta; |
| 840 | fPhiT[part] = phi; |
| 841 | fPdgT[part] = mpart; |
| 842 | |
| 843 | // final state particles only |
| 844 | |
| 845 | if (genType == kPythia) { |
| 846 | if (mPart->GetStatusCode() != 1) continue; |
| 847 | } else if (genType == kHijing) { |
| 848 | if (child1 >= 0 && child1 < npart) continue; |
| 849 | } |
| 850 | |
| 851 | |
| 852 | TParticlePDG* pdgP = 0; |
| 853 | // charged or neutral |
| 854 | pdgP = mPart->GetPDG(); |
| 855 | chTmp = pdgP->Charge() / 3.; // 13-feb-2001!! |
| 856 | |
| 857 | if (ich == 0) { |
| 858 | if (chTmp == 0) { |
| 859 | if (!fK0N) { |
| 860 | continue; |
| 861 | } else { |
| 862 | if (mpart != kNeutron && |
| 863 | mpart != kNeutronBar && |
| 864 | mpart != kK0Long) continue; |
| 865 | } |
| 866 | } |
| 867 | } else if (ich == 2) { |
| 868 | if (mpart == kNeutron || |
| 869 | mpart == kNeutronBar || |
| 870 | mpart == kK0Long) continue; |
| 871 | } |
| 872 | |
| 873 | if (TMath::Abs(eta)<=0.9) fNChTpc++; |
| 874 | // final state only |
| 875 | if (child1 >= 0 && child1 < npart) continue; |
| 876 | // acceptance cut |
| 877 | if (eta > fEtaMax || eta < fEtaMin) continue; |
| 878 | if (phi > fPhiMax || phi < fPhiMin ) continue; |
| 879 | // |
| 880 | if (fDebug >= 3) |
| 881 | printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ", |
| 882 | part, mpart, child1, eta, phi, pT, chTmp); |
| 883 | // |
| 884 | // |
| 885 | // Momentum smearing goes here ... |
| 886 | // |
| 887 | fhTrackPt->Fill(pT); |
| 888 | Float_t pw; |
| 889 | if (fSmear && TMath::Abs(chTmp)) { |
| 890 | pw = AliEMCALFast::SmearMomentum(1,p); |
| 891 | // p changed - take into account when calculate pT, |
| 892 | // pz and so on .. ? |
| 893 | pT = (pw/p) * pT; |
| 894 | if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw); |
| 895 | p = pw; |
| 896 | } |
| 897 | // |
| 898 | // Tracking Efficiency and TPC acceptance goes here ... |
| 899 | Float_t eff; |
| 900 | if (fEffic && TMath::Abs(chTmp)) { |
| 901 | eff = 0.9; // AliEMCALFast::Efficiency(2,p); |
| 902 | if(fhEff) fhEff->Fill(p, eff); |
| 903 | if (AliEMCALFast::RandomReject(eff)) { |
| 904 | if(fDebug >= 5) printf(" reject due to unefficiency "); |
| 905 | continue; |
| 906 | } |
| 907 | } |
| 908 | // |
| 909 | // Correction of Hadronic Energy goes here ... |
| 910 | // |
| 911 | // |
| 912 | // phi propagation for hadronic correction |
| 913 | |
| 914 | Bool_t curls = kFALSE; // hit two the EMCAL (no curl) |
| 915 | Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0; |
| 916 | if(TMath::Abs(chTmp)) { |
| 917 | // hadr. correction only for charge particle |
| 918 | dphi = PropagatePhi(pT, chTmp, curls); |
| 919 | phiHC = phi + dphi; |
| 920 | if (fDebug >= 6) { |
| 921 | printf("\n Delta phi %f pT %f ", dphi, pT); |
| 922 | if (curls) printf("\n !! Track is curling"); |
| 923 | } |
| 924 | if(!curls) fhTrackPtBcut->Fill(pT); |
| 925 | |
| 926 | if (fHCorrection && !curls) { |
| 927 | if (!fHadronCorrector) |
| 928 | Fatal("AliEMCALJetFinder", |
| 929 | "Hadronic energy correction required but not defined !"); |
| 930 | |
| 931 | dpH = fHadronCorrector->GetEnergy(p, eta, 7); |
| 932 | eTdpH = dpH*TMath::Sin(theta); |
| 933 | |
| 934 | if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n", |
| 935 | phi, phiHC, -eTdpH); // correction is negative |
| 936 | Int_t xbin,ybin; |
| 937 | xbin = fLego->GetXaxis()->FindBin(eta); |
| 938 | ybin = fLego->GetYaxis()->FindBin(phiHC); |
| 939 | cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl; |
| 940 | fLego->Fill(eta, phi, -fSamplingF*eTdpH ); |
| 941 | cout <<"Hadron Correction affected bin - contents after correction : "<<fLego->GetBinContent(xbin,ybin)<<endl; |
| 942 | fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH); |
| 943 | } |
| 944 | } |
| 945 | // |
| 946 | // More to do ? Just think about it ! |
| 947 | // |
| 948 | if (phi > fPhiMax || phi < fPhiMin ) continue; |
| 949 | |
| 950 | if(TMath::Abs(chTmp) ) { // charge particle |
| 951 | if (pT > fPtCut && !curls) { |
| 952 | if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", |
| 953 | eta , phi, pT, fNtS); |
| 954 | fLego->Fill(eta, phi, pT); |
| 955 | fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking |
| 956 | fTrackList[part] = 1; |
| 957 | fNtS++; |
| 958 | } |
| 959 | } else if(ich > 0 || fK0N) { |
| 960 | // case of n, nbar and K0L |
| 961 | if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", |
| 962 | eta , phi, pT, fNtS); |
| 963 | fLego->Fill(eta, phi, pT); |
| 964 | fTrackList[part] = 1; |
| 965 | fNtS++; |
| 966 | } |
| 967 | } // primary loop |
| 968 | Float_t etsum = 0.; |
| 969 | for(Int_t i=0; i<fLego->GetSize(); i++) { |
| 970 | Float_t etc = (*fLego)[i]; |
| 971 | if (etc > fMinCellEt) etsum += etc; |
| 972 | } |
| 973 | |
| 974 | printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral()); |
| 975 | printf(" Track selected(fNtS) %i \n", fNtS); |
| 976 | |
| 977 | DumpLego(); |
| 978 | } |
| 979 | |
| 980 | void AliEMCALJetFinder::FillFromHits(Int_t flag) |
| 981 | { |
| 982 | // |
| 983 | // Fill Cells with hit information |
| 984 | // |
| 985 | // |
| 986 | if (fDebug >= 2) |
| 987 | printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag); |
| 988 | |
| 989 | ResetMap(); |
| 990 | |
| 991 | if (!fLego) BookLego(); |
| 992 | // Reset eta-phi maps if needed |
| 993 | if (flag == 0) { // default behavior |
| 994 | fLego->Reset(); |
| 995 | fhLegoTracks->Reset(); |
| 996 | fhLegoEMCAL->Reset(); |
| 997 | fhLegoHadrCorr->Reset(); |
| 998 | } |
| 999 | // Initialize from background event if available |
| 1000 | // |
| 1001 | // Access hit information |
| 1002 | AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); |
| 1003 | AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"); |
| 1004 | TTree *treeH = emcalLoader->TreeH(); |
| 1005 | Int_t ntracks = (Int_t) treeH->GetEntries(); |
| 1006 | // |
| 1007 | // Loop over tracks |
| 1008 | // |
| 1009 | Int_t nbytes = 0; |
| 1010 | // Double_t etH = 0.0; |
| 1011 | |
| 1012 | for (Int_t track=0; track<ntracks;track++) { |
| 1013 | gAlice->ResetHits(); |
| 1014 | nbytes += treeH->GetEvent(track); |
| 1015 | // |
| 1016 | // Loop over hits |
| 1017 | // |
| 1018 | for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); |
| 1019 | mHit; |
| 1020 | mHit=(AliEMCALHit*) pEMCAL->NextHit()) |
| 1021 | { |
| 1022 | Float_t x = mHit->X(); // x-pos of hit |
| 1023 | Float_t y = mHit->Y(); // y-pos |
| 1024 | Float_t z = mHit->Z(); // z-pos |
| 1025 | Float_t eloss = mHit->GetEnergy(); // deposited energy |
| 1026 | |
| 1027 | Float_t r = TMath::Sqrt(x*x+y*y); |
| 1028 | Float_t theta = TMath::ATan2(r,z); |
| 1029 | Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); |
| 1030 | Float_t phi = TMath::ATan2(y,x); |
| 1031 | |
| 1032 | if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF); |
| 1033 | |
| 1034 | // etH = fSamplingF*eloss*TMath::Sin(theta); |
| 1035 | fLego->Fill(eta, phi, eloss); |
| 1036 | } // Hit Loop |
| 1037 | } // Track Loop |
| 1038 | |
| 1039 | // Transition from deposit energy to eT (eT = de*SF*sin(theta)) |
| 1040 | Double_t etsum = 0; |
| 1041 | for(Int_t i=1; i<=fLego->GetNbinsX(); i++){ // eta |
| 1042 | Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0; |
| 1043 | for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi |
| 1044 | eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta; |
| 1045 | fLego->SetBinContent(i,j,eT); |
| 1046 | // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical) |
| 1047 | fhLegoEMCAL->SetBinContent(i,j,eT); |
| 1048 | if (eT > fMinCellEt) etsum += eT; |
| 1049 | } |
| 1050 | } |
| 1051 | |
| 1052 | // for(Int_t i=0; i<fLego->GetSize(); i++) { |
| 1053 | // (*fhLegoEMCAL)[i] = (*fLego)[i]; |
| 1054 | // Float_t etc = (*fLego)[i]; |
| 1055 | // if (etc > fMinCellEt) etsum += etc; |
| 1056 | // } |
| 1057 | |
| 1058 | printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum); |
| 1059 | // DumpLego(); ?? |
| 1060 | } |
| 1061 | |
| 1062 | void AliEMCALJetFinder::FillFromDigits(Int_t flag) |
| 1063 | { |
| 1064 | // |
| 1065 | // Fill Cells with digit information |
| 1066 | // |
| 1067 | // |
| 1068 | ResetMap(); |
| 1069 | |
| 1070 | if (!fLego) BookLego(); |
| 1071 | if (flag == 0) fLego->Reset(); |
| 1072 | Int_t nbytes=0; |
| 1073 | |
| 1074 | |
| 1075 | // |
| 1076 | // Connect digits |
| 1077 | // |
| 1078 | TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000); |
| 1079 | TTree *treeD = gAlice->TreeD(); |
| 1080 | TBranchElement* branchDg = (TBranchElement*) |
| 1081 | treeD->GetBranch("EMCAL"); |
| 1082 | |
| 1083 | if (!branchDg) Fatal("AliEMCALJetFinder", |
| 1084 | "Reading digits requested but no digits in file !"); |
| 1085 | |
| 1086 | branchDg->SetAddress(&digs); |
| 1087 | Int_t nent = (Int_t) branchDg->GetEntries(); |
| 1088 | // |
| 1089 | // Connect digitizer |
| 1090 | // |
| 1091 | AliEMCALDigitizer* digr = new AliEMCALDigitizer(); |
| 1092 | TBranchElement* branchDr = (TBranchElement*) |
| 1093 | treeD->GetBranch("AliEMCALDigitizer"); |
| 1094 | branchDr->SetAddress(&digr); |
| 1095 | // |
| 1096 | // |
| 1097 | nbytes += branchDg->GetEntry(0); |
| 1098 | nbytes += branchDr->GetEntry(0); |
| 1099 | // |
| 1100 | // Get digitizer parameters |
| 1101 | Float_t ecADCped = digr->GetECApedestal(); |
| 1102 | Float_t ecADCcha = digr->GetECAchannel(); |
| 1103 | |
| 1104 | AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); |
| 1105 | AliEMCALGeometry* geom = |
| 1106 | AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), ""); |
| 1107 | |
| 1108 | if (fDebug) { |
| 1109 | Int_t ndig = digs->GetEntries(); |
| 1110 | Info("FillFromDigits","Number of Digits: %d %d\n Parameters: EC: %f %f\n Geometry: %d %d", |
| 1111 | ndig, nent, ecADCped, ecADCcha, geom->GetNEta(), geom->GetNPhi()); |
| 1112 | } |
| 1113 | |
| 1114 | // |
| 1115 | // Loop over digits |
| 1116 | AliEMCALDigit* sdg; |
| 1117 | TIter next(digs); |
| 1118 | while ((sdg = (AliEMCALDigit*)(next()))) |
| 1119 | { |
| 1120 | Double_t pedestal = 0.; |
| 1121 | Double_t channel = 0.; |
| dc7da436 |
1122 | if (geom->CheckAbsCellId(sdg->GetId())) { // May 31, 2006; IsInECA(Int_t) -> CheckAbsCellId |
| 45a58699 |
1123 | pedestal = ecADCped; |
| 1124 | channel = ecADCcha; |
| 1125 | } |
| 1126 | else |
| 1127 | Fatal("FillFromDigits", "unexpected digit is number!") ; |
| 1128 | |
| 1129 | Float_t eta = sdg->GetEta(); |
| 1130 | Float_t phi = sdg->GetPhi() * TMath::Pi()/180.; |
| 1131 | Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal); |
| 1132 | |
| 1133 | if (fDebug > 1) |
| 1134 | Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d", |
| 1135 | eta, phi, amp, sdg->GetAmp()); |
| 1136 | |
| 1137 | fLego->Fill(eta, phi, fSamplingF*amp); |
| 1138 | } // digit loop |
| 1139 | // |
| 1140 | // Dump lego hist |
| 1141 | DumpLego(); |
| 1142 | } |
| 1143 | |
| 1144 | |
| 1145 | void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag) |
| 1146 | { |
| 1147 | // |
| 1148 | // Fill Cells with hit information |
| 1149 | // |
| 1150 | // |
| 1151 | ResetMap(); |
| 1152 | |
| 1153 | if (!fLego) BookLego(); |
| 1154 | // Reset |
| 1155 | if (flag == 0) fLego->Reset(); |
| 1156 | |
| 1157 | // Flag charged tracks passing through TPC which |
| 1158 | // are associated to EMCAL Hits |
| 1159 | BuildTrackFlagTable(); |
| 1160 | |
| 1161 | // |
| 1162 | // Access particle information |
| 1163 | TTree *treeK = gAlice->TreeK(); |
| 1164 | Int_t ntracks = (Int_t) treeK->GetEntries(); |
| 1165 | |
| 1166 | if (fPtT) delete[] fPtT; |
| 1167 | if (fEtaT) delete[] fEtaT; |
| 1168 | if (fPhiT) delete[] fPhiT; |
| 1169 | if (fPdgT) delete[] fPdgT; |
| 1170 | |
| 1171 | fPtT = new Float_t[ntracks]; |
| 1172 | fEtaT = new Float_t[ntracks]; |
| 1173 | fPhiT = new Float_t[ntracks]; |
| 1174 | fPdgT = new Int_t[ntracks]; |
| 1175 | |
| 1176 | fNt = ntracks; |
| 1177 | fNtS = 0; |
| 1178 | |
| 1179 | for (Int_t track = 0; track < ntracks; track++) { |
| 1180 | TParticle *mPart = gAlice->GetMCApp()->Particle(track); |
| 1181 | Float_t pT = mPart->Pt(); |
| 1182 | Float_t phi = mPart->Phi(); |
| 1183 | Float_t eta = mPart->Eta(); |
| 1184 | |
| 1185 | if(fTrackList[track]) { |
| 1186 | fPtT[track] = pT; |
| 1187 | fEtaT[track] = eta; |
| 1188 | fPhiT[track] = phi; |
| 1189 | fPdgT[track] = mPart->GetPdgCode(); |
| 1190 | |
| 1191 | if (track < 2) continue; //Colliding particles? |
| 1192 | if (pT == 0 || pT < fPtCut) continue; |
| 1193 | fNtS++; |
| 1194 | fLego->Fill(eta, phi, pT); |
| 1195 | } |
| 1196 | } // track loop |
| 1197 | DumpLego(); |
| 1198 | } |
| 1199 | |
| 1200 | void AliEMCALJetFinder::FillFromParticles() |
| 1201 | { |
| 1202 | // 26-feb-2002 PAI - for checking all chain |
| 1203 | // Work on particles level; accept all particle (not neutrino ) |
| 1204 | |
| 1205 | Double_t pX=0, pY=0, pZ=0, energy=0; // checking conservation law |
| 1206 | fNChTpc = 0; |
| 1207 | |
| 1208 | ResetMap(); |
| 1209 | if (!fLego) BookLego(); |
| 1210 | fLego->Reset(); |
| 1211 | // |
| 1212 | // Access particles information |
| 1213 | Int_t npart = (gAlice->GetHeader())->GetNprimary(); |
| 1214 | if (fDebug >= 2 || npart<=0) { |
| 1215 | printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart); |
| 1216 | if(npart<=0) return; |
| 1217 | } |
| 1218 | fNt = npart; |
| 1219 | fNtS = 0; |
| 1220 | RearrangeParticlesMemory(npart); |
| 1221 | |
| 1222 | // Go through the particles |
| 1223 | Int_t mpart, child1, child2, geantPdg; |
| 1224 | Float_t pT, phi, eta, e=0, px=0, py=0, pz=0; |
| 1225 | TParticle *mPart=0; |
| 1226 | for (Int_t part = 0; part < npart; part++) { |
| 1227 | |
| 1228 | fTrackList[part] = 0; |
| 1229 | |
| 1230 | mPart = gAlice->GetMCApp()->Particle(part); |
| 1231 | mpart = mPart->GetPdgCode(); |
| 1232 | child1 = mPart->GetFirstDaughter(); |
| 1233 | child2 = mPart->GetLastDaughter(); |
| 1234 | pT = mPart->Pt(); |
| 1235 | phi = mPart->Phi(); |
| 1236 | eta = mPart->Eta(); |
| 1237 | |
| 1238 | px = mPart->Px(); |
| 1239 | py = mPart->Py(); |
| 1240 | pz = mPart->Pz(); |
| 1241 | e = mPart->Energy(); |
| 1242 | |
| 1243 | // see pyedit in Pythia's text |
| 1244 | geantPdg = mpart; |
| 1245 | if (IsThisPartonsOrDiQuark(mpart)) continue; |
| 1246 | printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n", |
| 1247 | part, mpart, geantPdg, px, py, pz, e, child1, child2); |
| 1248 | |
| 1249 | // exclude partons (21 - gluon, 92 - string) |
| 1250 | |
| 1251 | |
| 1252 | // exclude neutrinous also ?? |
| 1253 | if (fDebug >= 11 && pT>0.01) |
| 1254 | printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ", |
| 1255 | part, mpart, eta, phi, pT); |
| 1256 | |
| 1257 | fPtT[part] = pT; |
| 1258 | fEtaT[part] = eta; |
| 1259 | fPhiT[part] = phi; |
| 1260 | fPdgT[part] = mpart; |
| 1261 | fNtS++; |
| 1262 | |
| 1263 | // final state only |
| 1264 | if (child1 >= 0 && child1 < npart) continue; |
| 1265 | |
| 1266 | // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n", |
| 1267 | // part, mpart, geantPdg, px, py, pz, e, child1, name.Data()); |
| 1268 | pX += px; |
| 1269 | pY += py; |
| 1270 | pZ += pz; |
| 1271 | energy += e; |
| 1272 | |
| 1273 | |
| 1274 | if (TMath::Abs(eta) <= 0.9) fNChTpc++; |
| 1275 | // acceptance cut |
| 1276 | if (eta > fEtaMax || eta < fEtaMin) continue; |
| 1277 | if (phi > fPhiMax || phi < fPhiMin ) continue; |
| 1278 | // |
| 1279 | if(fK0N==0 ) { // exclude neutral hadrons |
| 1280 | if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue; |
| 1281 | } |
| 1282 | fTrackList[part] = 1; |
| 1283 | fLego->Fill(eta, phi, pT); |
| 1284 | |
| 1285 | } // primary loop |
| 1286 | printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n", |
| 1287 | pX, pY, pZ, energy); |
| 1288 | DumpLego(); |
| 1289 | if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc); |
| 1290 | } |
| 1291 | |
| 1292 | void AliEMCALJetFinder::FillFromPartons() |
| 1293 | { |
| 1294 | // function under construction - 13-feb-2002 PAI |
| 1295 | |
| 1296 | if (fDebug >= 2) |
| 1297 | printf("\n AliEMCALJetFinder::FillFromPartons()\n"); |
| 1298 | // |
| 1299 | |
| 1300 | ResetMap(); |
| 1301 | if (!fLego) BookLego(); |
| 1302 | fLego->Reset(); |
| 1303 | // |
| 1304 | // Access particle information |
| 1305 | Int_t npart = (gAlice->GetHeader())->GetNprimary(); |
| 1306 | if (fDebug >= 2 || npart<=0) |
| 1307 | printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart); |
| 1308 | fNt = 0; // for FindTracksInJetCone |
| 1309 | fNtS = 0; |
| 1310 | |
| 1311 | // Go through the partons |
| 1312 | Int_t statusCode=0; |
| 1313 | for (Int_t part = 8; part < npart; part++) { |
| 1314 | TParticle *mPart = gAlice->GetMCApp()->Particle(part); |
| 1315 | Int_t mpart = mPart->GetPdgCode(); |
| 1316 | // Int_t child1 = MPart->GetFirstDaughter(); |
| 1317 | Float_t pT = mPart->Pt(); |
| 1318 | // Float_t p = MPart->P(); |
| 1319 | Float_t phi = mPart->Phi(); |
| 1320 | Float_t eta = mPart->Eta(); |
| 1321 | // Float_t theta = MPart->Theta(); |
| 1322 | statusCode = mPart->GetStatusCode(); |
| 1323 | |
| 1324 | // accept partons (21 - gluon, 92 - string) |
| 1325 | if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue; |
| 1326 | if (fDebug > 1 && pT>0.01) |
| 1327 | printf("\n part:%5d mpart %5d status %5d eta %8.2f phi %8.2f pT %8.2f ", |
| 1328 | part, mpart, statusCode, eta, phi, pT); |
| 1329 | // if (fDebug >= 3) MPart->Print(); |
| 1330 | // accept partons before fragmentation - p.57 in Pythia manual |
| 1331 | // if(statusCode != 1) continue; |
| 1332 | // acceptance cut |
| 1333 | if (eta > fEtaMax || eta < fEtaMin) continue; |
| 1334 | if (phi > fPhiMax || phi < fPhiMin ) continue; |
| 1335 | // final state only |
| 1336 | // if (child1 >= 0 && child1 < npart) continue; |
| 1337 | // |
| 1338 | // |
| 1339 | fLego->Fill(eta, phi, pT); |
| 1340 | |
| 1341 | } // primary loop |
| 1342 | DumpLego(); |
| 1343 | } |
| 1344 | |
| 1345 | void AliEMCALJetFinder::BuildTrackFlagTable() { |
| 1346 | |
| 1347 | // Method to generate a lookup table for TreeK particles |
| 1348 | // which are linked to hits in the EMCAL |
| 1349 | // |
| 1350 | // --Author: J.L. Klay |
| 1351 | // |
| 1352 | // Access hit information |
| 1353 | AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); |
| 1354 | |
| 1355 | TTree *treeK = gAlice->TreeK(); // Get the number of entries in the kine tree |
| 1356 | Int_t nKTrks = (Int_t) treeK->GetEntries(); // (Number of particles created somewhere) |
| 1357 | |
| 1358 | if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one |
| 1359 | fTrackList = new Int_t[nKTrks]; //before generating a new one |
| 1360 | |
| 1361 | for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0 |
| 1362 | fTrackList[i] = 0; |
| 1363 | } |
| 1364 | |
| 1365 | AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"); |
| 1366 | //AliEMCALGetter * gime = AliEMCALGetter::Instance() ; |
| 1367 | // TTree *treeH = gAlice->TreeH(); |
| 1368 | TTree *treeH = emcalLoader->TreeH(); |
| 1369 | Int_t ntracks = (Int_t) treeH->GetEntries(); |
| 1370 | // |
| 1371 | // Loop over tracks |
| 1372 | // |
| 1373 | Int_t nbytes = 0; |
| 1374 | |
| 1375 | for (Int_t track=0; track<ntracks;track++) { |
| 1376 | gAlice->ResetHits(); |
| 1377 | nbytes += treeH->GetEvent(track); |
| 1378 | if (pEMCAL) { |
| 1379 | |
| 1380 | // |
| 1381 | // Loop over hits |
| 1382 | // |
| 1383 | for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1); |
| 1384 | mHit; |
| 1385 | mHit=(AliEMCALHit*) pEMCAL->NextHit()) |
| 1386 | { |
| 1387 | Int_t iTrk = mHit->Track(); // track number |
| 1388 | Int_t idprim = mHit->GetPrimary(); // primary particle |
| 1389 | |
| 1390 | //Determine the origin point of this particle - it made a hit in the EMCAL |
| 1391 | TParticle *trkPart = gAlice->GetMCApp()->Particle(iTrk); |
| 1392 | TParticlePDG *trkPdg = trkPart->GetPDG(); |
| 1393 | Int_t trkCode = trkPart->GetPdgCode(); |
| 1394 | Double_t trkChg; |
| 1395 | if (trkCode < 10000) { //Big Ions cause problems for |
| 1396 | trkChg = trkPdg->Charge(); //this function. Since they aren't |
| 1397 | } else { //likely to make it very far, set |
| 1398 | trkChg = 0.0; //their charge to 0 for the Flag test |
| 1399 | } |
| 1400 | Float_t vxTrk = trkPart->Vx(); |
| 1401 | Float_t vyTrk = trkPart->Vy(); |
| 1402 | Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk); |
| 1403 | fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg); |
| 1404 | |
| 1405 | //Loop through the ancestry of the EMCAL entrance particles |
| 1406 | Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother |
| 1407 | while (ancestor != -1) { |
| 1408 | TParticle *ancPart = gAlice->GetMCApp()->Particle(ancestor); //get particle info on ancestor |
| 1409 | TParticlePDG *ancPdg = ancPart->GetPDG(); |
| 1410 | Int_t ancCode = ancPart->GetPdgCode(); |
| 1411 | Double_t ancChg; |
| 1412 | if (ancCode < 10000) { |
| 1413 | ancChg = ancPdg->Charge(); |
| 1414 | } else { |
| 1415 | ancChg = 0.0; |
| 1416 | } |
| 1417 | Float_t vxAnc = ancPart->Vx(); |
| 1418 | Float_t vyAnc = ancPart->Vy(); |
| 1419 | Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc); |
| 1420 | fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg); |
| 1421 | ancestor = ancPart->GetFirstMother(); //Get the next ancestor |
| 1422 | } |
| 1423 | |
| 1424 | //Determine the origin point of the primary particle associated with the hit |
| 1425 | TParticle *primPart = gAlice->GetMCApp()->Particle(idprim); |
| 1426 | TParticlePDG *primPdg = primPart->GetPDG(); |
| 1427 | Int_t primCode = primPart->GetPdgCode(); |
| 1428 | Double_t primChg; |
| 1429 | if (primCode < 10000) { |
| 1430 | primChg = primPdg->Charge(); |
| 1431 | } else { |
| 1432 | primChg = 0.0; |
| 1433 | } |
| 1434 | Float_t vxPrim = primPart->Vx(); |
| 1435 | Float_t vyPrim = primPart->Vy(); |
| 1436 | Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim); |
| 1437 | fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg); |
| 1438 | } // Hit Loop |
| 1439 | } //if (pEMCAL) |
| 1440 | } // Track Loop |
| 1441 | } |
| 1442 | |
| 1443 | Int_t AliEMCALJetFinder |
| 1444 | ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) { |
| 1445 | // Set the flag for the track |
| 1446 | Int_t flag = 0; |
| 1447 | Int_t parton = 0; |
| 1448 | Int_t neutral = 0; |
| 1449 | |
| 1450 | if (charge == 0) neutral = 1; |
| 1451 | |
| 1452 | if (TMath::Abs(code) <= 6 || |
| 1453 | code == 21 || |
| 1454 | code == 92) parton = 1; |
| 1455 | |
| 1456 | //It's not a parton, it's charged and it went through the TPC |
| 1457 | if (!parton && !neutral && radius <= 84.0) flag = 1; |
| 1458 | |
| 1459 | return flag; |
| 1460 | } |
| 1461 | |
| 1462 | |
| 1463 | |
| 1464 | void AliEMCALJetFinder::SaveBackgroundEvent(const char *name) |
| 1465 | { |
| 1466 | // Saves the eta-phi lego and the tracklist and name of file with BG events |
| 1467 | // |
| 1468 | if (fLegoB) { |
| 1469 | fLegoB->Reset(); |
| 1470 | (*fLegoB) = (*fLegoB) + (*fLego); |
| 1471 | // if(fDebug) |
| 1472 | printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n", |
| 1473 | fLegoB->Integral(), fLego->Integral()); |
| 1474 | } |
| 1475 | |
| 1476 | if (fPtB) delete[] fPtB; |
| 1477 | if (fEtaB) delete[] fEtaB; |
| 1478 | if (fPhiB) delete[] fPhiB; |
| 1479 | if (fPdgB) delete[] fPdgB; |
| 1480 | if (fTrackListB) delete[] fTrackListB; |
| 1481 | |
| 1482 | fPtB = new Float_t[fNtS]; |
| 1483 | fEtaB = new Float_t[fNtS]; |
| 1484 | fPhiB = new Float_t[fNtS]; |
| 1485 | fPdgB = new Int_t [fNtS]; |
| 1486 | fTrackListB = new Int_t [fNtS]; |
| 1487 | |
| 1488 | fNtB = 0; |
| 1489 | |
| 1490 | for (Int_t i = 0; i < fNt; i++) { |
| 1491 | if (!fTrackList[i]) continue; |
| 1492 | fPtB [fNtB] = fPtT [i]; |
| 1493 | fEtaB[fNtB] = fEtaT[i]; |
| 1494 | fPhiB[fNtB] = fPhiT[i]; |
| 1495 | fPdgB[fNtB] = fPdgT[i]; |
| 1496 | fTrackListB[fNtB] = 1; |
| 1497 | fNtB++; |
| 1498 | } |
| 1499 | fBackground = 1; |
| 1500 | printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt); |
| 1501 | |
| 1502 | if(strlen(name) == 0) { |
| 1503 | TSeqCollection *li = gROOT->GetListOfFiles(); |
| 1504 | TString nf; |
| 1505 | for(Int_t i=0; i<li->GetSize(); i++) { |
| 1506 | nf = ((TFile*)li->At(i))->GetName(); |
| 1507 | if(nf.Contains("backgorund")) break; |
| 1508 | } |
| 1509 | fBGFileName = nf; |
| 1510 | } else { |
| 1511 | fBGFileName = name; |
| 1512 | } |
| 1513 | printf("BG file name is \n %s\n", fBGFileName.Data()); |
| 1514 | } |
| 1515 | |
| 1516 | void AliEMCALJetFinder::InitFromBackground() |
| 1517 | { |
| 1518 | // |
| 1519 | // |
| 1520 | if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() "); |
| 1521 | |
| 1522 | if (fLego) { |
| 1523 | fLego->Reset(); |
| 1524 | (*fLego) = (*fLego) + (*fLegoB); |
| 1525 | if(fDebug) |
| 1526 | printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n", |
| 1527 | fLego->Integral(), fLegoB->Integral()); |
| 1528 | } else { |
| 1529 | printf(" => fLego undefined \n"); |
| 1530 | } |
| 1531 | } |
| 1532 | |
| 1533 | |
| 1534 | void AliEMCALJetFinder::FindTracksInJetCone() |
| 1535 | { |
| 1536 | // |
| 1537 | // Build list of tracks inside jet cone |
| 1538 | // |
| 1539 | // Loop over jets |
| 1540 | Int_t njet = Njets(); |
| 1541 | for (Int_t nj = 0; nj < njet; nj++) |
| 1542 | { |
| 1543 | Float_t etaj = JetEtaW(nj); |
| 1544 | Float_t phij = JetPhiW(nj); |
| 1545 | Int_t nT = 0; // number of associated tracks |
| 1546 | |
| 1547 | // Loop over particles in current event |
| 1548 | for (Int_t part = 0; part < fNt; part++) { |
| 1549 | if (!fTrackList[part]) continue; |
| 1550 | Float_t phi = fPhiT[part]; |
| 1551 | Float_t eta = fEtaT[part]; |
| 1552 | Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) + |
| 1553 | (phij-phi)*(phij-phi)); |
| 1554 | if (dr < fConeRadius) { |
| 1555 | fTrackList[part] = nj+2; |
| 1556 | nT++; |
| 1557 | } // < ConeRadius ? |
| 1558 | } // particle loop |
| 1559 | |
| 1560 | // Same for background event if available |
| 1561 | Int_t nTB = 0; |
| 1562 | if (fBackground) { |
| 1563 | for (Int_t part = 0; part < fNtB; part++) { |
| 1564 | Float_t phi = fPhiB[part]; |
| 1565 | Float_t eta = fEtaB[part]; |
| 1566 | Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) + |
| 1567 | (phij-phi)*(phij-phi)); |
| 1568 | fTrackListB[part] = 1; |
| 1569 | |
| 1570 | if (dr < fConeRadius) { |
| 1571 | fTrackListB[part] = nj+2; |
| 1572 | nTB++; |
| 1573 | } // < ConeRadius ? |
| 1574 | } // particle loop |
| 1575 | } // Background available ? |
| 1576 | |
| 1577 | Int_t nT0 = nT + nTB; |
| 1578 | printf("Total number of tracks %d\n", nT0); |
| 1579 | |
| 1580 | if (nT0 > 1000) nT0 = 1000; |
| 1581 | |
| 1582 | Float_t* ptT = new Float_t[nT0]; |
| 1583 | Float_t* etaT = new Float_t[nT0]; |
| 1584 | Float_t* phiT = new Float_t[nT0]; |
| 1585 | Int_t* pdgT = new Int_t[nT0]; |
| 1586 | |
| 1587 | Int_t iT = 0; |
| 1588 | Int_t j; |
| 1589 | |
| 1590 | for (Int_t part = 0; part < fNt; part++) { |
| 1591 | if (fTrackList[part] == nj+2) { |
| 1592 | Int_t index = 0; |
| 1593 | for (j=iT-1; j>=0; j--) { |
| 1594 | if (fPtT[part] > ptT[j]) { |
| 1595 | index = j+1; |
| 1596 | break; |
| 1597 | } |
| 1598 | } |
| 1599 | for (j=iT-1; j>=index; j--) { |
| 1600 | ptT [j+1] = ptT [j]; |
| 1601 | etaT[j+1] = etaT[j]; |
| 1602 | phiT[j+1] = phiT[j]; |
| 1603 | pdgT[j+1] = pdgT[j]; |
| 1604 | } |
| 1605 | ptT [index] = fPtT [part]; |
| 1606 | etaT[index] = fEtaT[part]; |
| 1607 | phiT[index] = fPhiT[part]; |
| 1608 | pdgT[index] = fPdgT[part]; |
| 1609 | iT++; |
| 1610 | } // particle associated |
| 1611 | if (iT > nT0) break; |
| 1612 | } // particle loop |
| 1613 | |
| 1614 | if (fBackground) { |
| 1615 | for (Int_t part = 0; part < fNtB; part++) { |
| 1616 | if (fTrackListB[part] == nj+2) { |
| 1617 | Int_t index = 0; |
| 1618 | for (j=iT-1; j>=0; j--) { |
| 1619 | if (fPtB[part] > ptT[j]) { |
| 1620 | index = j+1; |
| 1621 | |
| 1622 | break; |
| 1623 | } |
| 1624 | } |
| 1625 | for (j=iT-1; j>=index; j--) { |
| 1626 | ptT [j+1] = ptT [j]; |
| 1627 | etaT[j+1] = etaT[j]; |
| 1628 | phiT[j+1] = phiT[j]; |
| 1629 | pdgT[j+1] = pdgT[j]; |
| 1630 | } |
| 1631 | ptT [index] = fPtB [part]; |
| 1632 | etaT[index] = fEtaB[part]; |
| 1633 | phiT[index] = fPhiB[part]; |
| 1634 | pdgT[index] = fPdgB[part]; |
| 1635 | iT++; |
| 1636 | } // particle associated |
| 1637 | if (iT > nT0) break; |
| 1638 | } // particle loop |
| 1639 | } // Background available ? |
| 1640 | |
| 1641 | fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT); |
| 1642 | delete[] ptT; |
| 1643 | delete[] etaT; |
| 1644 | delete[] phiT; |
| 1645 | delete[] pdgT; |
| 1646 | |
| 1647 | } // jet loop loop |
| 1648 | } |
| 1649 | |
| 1650 | Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls) |
| 1651 | { |
| 1652 | // Propagates phi angle to EMCAL radius |
| 1653 | // |
| 1654 | static Float_t b = 0.0, rEMCAL = -1.0; |
| 1655 | // Get field in kGS |
| f7a1cc68 |
1656 | b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField(); |
| 45a58699 |
1657 | // Get EMCAL radius in cm |
| 1658 | rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance(); |
| 1659 | Float_t dPhi = 0.; |
| 1660 | // |
| 1661 | // |
| 1662 | // bending radies |
| 1663 | // pt [Gev] |
| 1664 | // B [kG] |
| 1665 | // |
| 4bba84bd |
1666 | Float_t rB = 3335.6 * pt / TMath::Abs(b); // [cm] (case of |charge|=1) |
| 45a58699 |
1667 | // |
| 1668 | // check if particle is curling below EMCAL |
| 1669 | if (2.*rB < rEMCAL) { |
| 1670 | curls = kTRUE; |
| 4bba84bd |
1671 | AliDebug(1, Form(" Low pt %f \n", pt)); |
| 45a58699 |
1672 | return dPhi; |
| 1673 | } |
| 1674 | // |
| 1675 | // if not calculate delta phi |
| 1676 | Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB)); |
| 1677 | dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi)); |
| 1678 | dPhi = -TMath::Sign(dPhi, charge); |
| 1679 | // |
| 4bba84bd |
1680 | AliDebug(1, Form(" B %7.3f kG : rEMCAL %7.2f : dphi %7.5f(%7.5f)\n", |
| 1681 | b, rEMCAL, dPhi,dPhi*TMath::RadToDeg())); |
| 45a58699 |
1682 | return dPhi; |
| 1683 | } |
| 1684 | |
| 1685 | void hf1(Int_t& , Float_t& , Float_t& ) |
| 1686 | { |
| 1687 | // dummy for hbook calls |
| 1688 | ; |
| 1689 | } |
| 1690 | |
| 1691 | void AliEMCALJetFinder::DrawLego(const char *opt) |
| 1692 | { |
| 1693 | // Draw lego |
| 1694 | if(fLego) fLego->Draw(opt); |
| 1695 | } |
| 1696 | |
| 1697 | void AliEMCALJetFinder::DrawLegoBackground(const char *opt) |
| 1698 | { |
| 1699 | // Draw background lego |
| 1700 | if(fLegoB) fLegoB->Draw(opt); |
| 1701 | } |
| 1702 | |
| 1703 | void AliEMCALJetFinder::DrawLegoEMCAL(const char *opt) |
| 1704 | { |
| 1705 | // Draw EMCAL Lego |
| 1706 | if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt); |
| 1707 | } |
| 1708 | |
| 1709 | void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode) |
| 1710 | { |
| 1711 | // Draw all hists |
| 1712 | static TPaveText *varLabel=0; |
| 1713 | if(!fC1) { |
| 1714 | fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800); |
| 1715 | } |
| 1716 | fC1->Clear(); |
| 1717 | fC1->Divide(2,2); |
| 1718 | fC1->cd(1); |
| 1719 | gPad->SetLogy(1); |
| 1720 | fhCellEt->Draw(); |
| 1721 | |
| 1722 | fC1->cd(2); |
| 1723 | gPad->SetLogy(1); |
| 1724 | fhCellEMCALEt->Draw(); |
| 1725 | |
| 1726 | fC1->cd(3); |
| 1727 | gPad->SetLogy(1); |
| 1728 | fhTrackPt->Draw(); |
| 1729 | fhTrackPtBcut->SetLineColor(2); |
| 1730 | fhTrackPtBcut->Draw("same"); |
| 1731 | |
| 1732 | fC1->cd(4); |
| 1733 | if(!varLabel) { |
| 1734 | PrintParameters(1); |
| 1735 | varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC"); |
| 1736 | varLabel->SetTextAlign(12); |
| 1737 | varLabel->SetFillColor(19); // see TAttFill |
| 1738 | TString tmp(GetTitle()); |
| 1739 | varLabel->ReadFile(GetFileNameForParameters()); |
| 1740 | } |
| 1741 | varLabel->Draw(); |
| 1742 | fC1->Update(); |
| 1743 | if(mode) { // for saving picture to the file |
| 1744 | TString stmp(GetFileNameForParameters()); |
| 1745 | stmp.ReplaceAll("_Par.txt",".ps"); |
| 1746 | fC1->Print(stmp.Data()); |
| 1747 | } |
| 1748 | } |
| 1749 | |
| 1750 | void AliEMCALJetFinder::PrintParameters(Int_t mode) |
| 1751 | { |
| 1752 | // Print all parameters out |
| 1753 | |
| 1754 | FILE *file=0; |
| 1755 | if(mode==0) file = stdout; // output to terminal |
| 1756 | else { |
| 1757 | file = fopen(GetFileNameForParameters(),"w"); |
| 1758 | if(file==0) file = stdout; |
| 1759 | } |
| 1760 | fprintf(file,"==== Filling lego ==== \n"); |
| 1761 | fprintf(file,"Sampling fraction %6.3f ", fSamplingF); |
| 1762 | fprintf(file,"Smearing %6i ", fSmear); |
| 1763 | fprintf(file,"Efficiency %6i\n", fEffic); |
| 1764 | fprintf(file,"Hadr.Correct. %6i ", fHCorrection); |
| 1765 | fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut); |
| 1766 | fprintf(file,"==== Jet finding ==== \n"); |
| 1767 | fprintf(file,"Cone radius %6.2f ", fConeRadius); |
| 1768 | fprintf(file,"Seed E_{T} %6.1f\n", fEtSeed); |
| 1769 | fprintf(file,"Min E_{T} of cell %6.1f ", fMinCellEt); |
| 1770 | fprintf(file,"Min E_{T} of jet %6.1f\n", fMinJetEt); |
| 1771 | if(fMode) { |
| 1772 | fprintf(file,"==== Bg subtraction ==== \n"); |
| 1773 | fprintf(file,"BG subtraction %6i ", fMode); |
| 1774 | fprintf(file,"Min cone move %6.2f\n", fMinMove); |
| 1775 | fprintf(file,"Max cone move %6.2f ", fMaxMove); |
| 1776 | fprintf(file,"%% change for BG %6.4f\n", fPrecBg); |
| 1777 | } else |
| 1778 | fprintf(file,"==== No Bg subtraction ==== \n"); |
| 1779 | // |
| 1780 | printf("BG file name is %s \n", fBGFileName.Data()); |
| 1781 | if(file != stdout) fclose(file); |
| 1782 | } |
| 1783 | |
| 1784 | void AliEMCALJetFinder::DrawLegos() |
| 1785 | { |
| 1786 | // Draw all legos |
| 1787 | if(!fC1) { |
| 1788 | fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800); |
| 1789 | } |
| 1790 | fC1->Clear(); |
| 1791 | fC1->Divide(2,2); |
| 1792 | gStyle->SetOptStat(111111); |
| 1793 | |
| 1794 | Int_t nent1, nent2, nent3, nent4; |
| 1795 | Double_t int1, int2, int3, int4; |
| 1796 | nent1 = (Int_t)fLego->GetEntries(); |
| 1797 | int1 = fLego->Integral(); |
| 1798 | fC1->cd(1); |
| 1799 | if(int1) fLego->Draw("lego"); |
| 1800 | |
| 1801 | nent2 = (Int_t)fhLegoTracks->GetEntries(); |
| 1802 | int2 = fhLegoTracks->Integral(); |
| 1803 | fC1->cd(2); |
| 1804 | if(int2) fhLegoTracks->Draw("lego"); |
| 1805 | |
| 1806 | nent3 = (Int_t)fhLegoEMCAL->GetEntries(); |
| 1807 | int3 = fhLegoEMCAL->Integral(); |
| 1808 | fC1->cd(3); |
| 1809 | if(int3) fhLegoEMCAL->Draw("lego"); |
| 1810 | |
| 1811 | nent4 = (Int_t)fhLegoHadrCorr->GetEntries(); |
| 1812 | int4 = fhLegoHadrCorr->Integral(); |
| 1813 | fC1->cd(4); |
| 1814 | if(int4) fhLegoHadrCorr->Draw("lego"); |
| 1815 | |
| 1816 | // just for checking |
| 1817 | printf(" Integrals \n"); |
| 1818 | printf("lego %10.3f \ntrack %10.3f \nhits %10.3f \nHCorr %10.3f\n-- %10.3f(must be 0)\n", |
| 1819 | int1, int2, int3, int4, int1 - (int2 + int3 - int4)); |
| 1820 | } |
| 1821 | |
| 1822 | const Char_t* AliEMCALJetFinder::GetFileNameForParameters(const char* dir) |
| 1823 | { |
| 1824 | // Get paramters from a file |
| 1825 | static TString tmp; |
| 1826 | if(strlen(dir)) tmp = dir; |
| 1827 | tmp += GetTitle(); |
| 1828 | tmp += "_Par.txt"; |
| 1829 | return tmp.Data(); |
| 1830 | } |
| 1831 | |
| 1832 | void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart) |
| 1833 | { |
| 1834 | // See FillFromTracks() - npart must be positive |
| 1835 | if (fTrackList) delete[] fTrackList; |
| 1836 | if (fPtT) delete[] fPtT; |
| 1837 | if (fEtaT) delete[] fEtaT; |
| 1838 | if (fPhiT) delete[] fPhiT; |
| 1839 | if (fPdgT) delete[] fPdgT; |
| 1840 | |
| 1841 | if(npart>0) { |
| 1842 | fTrackList = new Int_t [npart]; |
| 1843 | fPtT = new Float_t[npart]; |
| 1844 | fEtaT = new Float_t[npart]; |
| 1845 | fPhiT = new Float_t[npart]; |
| 1846 | fPdgT = new Int_t[npart]; |
| 1847 | } else { |
| 1848 | printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart); |
| 1849 | } |
| 1850 | } |
| 1851 | |
| 1852 | Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg) |
| 1853 | { |
| 1854 | // Return quark info |
| 1855 | Int_t absPdg = TMath::Abs(pdg); |
| 1856 | if(absPdg<=6) return kTRUE; // quarks |
| 1857 | if(pdg == 21) return kTRUE; // gluon |
| 1858 | if(pdg == 92) return kTRUE; // string |
| 1859 | |
| 1860 | // see p.51 of Pythia Manual |
| 1861 | // Not include diquarks with c and b quark - 4-mar-2002 |
| 1862 | // ud_0,sd_0,su_0; dd_1,ud_1,uu_1; sd_1,su_1,ss_1 |
| 1863 | static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203, 3103,3203,3303}; |
| 1864 | for(Int_t i=0; i<9; i++) if(absPdg == diquark[i]) return kTRUE; // diquarks |
| 1865 | |
| 1866 | return kFALSE; |
| 1867 | } |
| 1868 | |
| 1869 | void AliEMCALJetFinder::FindChargedJet() |
| 1870 | { |
| 1871 | // |
| 1872 | // Find jet from charged particle information only |
| 1873 | // |
| 1874 | |
| 1875 | // |
| 1876 | // Look for seeds |
| 1877 | // |
| 1878 | Int_t njets = 0; |
| 1879 | Int_t part = 0; |
| 1880 | Int_t nseed = 0; |
| 1881 | |
| 1882 | // |
| 1883 | // |
| 1884 | ResetJets(); |
| 1885 | |
| 1886 | // |
| 1887 | for (part = 0; part < fNt; part++) { |
| 1888 | if (!fTrackList[part]) continue; |
| 1889 | if (fPtT[part] > fEtSeed) nseed++; |
| 1890 | } |
| 1891 | printf("\nThere are %d seeds (%d)\n", nseed, fNtS); |
| 1892 | Int_t* iSeeds = new Int_t[nseed]; |
| 1893 | nseed = 0; |
| 1894 | |
| 1895 | for (part = 0; part < fNt; part++) { |
| 1896 | if (!fTrackList[part]) continue; |
| 1897 | if (fPtT[part] > fEtSeed) iSeeds[nseed++] = part; |
| 1898 | } |
| 1899 | |
| 1900 | // |
| 1901 | // Loop over seeds |
| 1902 | // |
| 1903 | Int_t seed = 0; |
| 1904 | Float_t pt; |
| 1905 | |
| 1906 | while(1){ |
| 1907 | // |
| 1908 | // Find seed with highest pt |
| 1909 | // |
| 1910 | Float_t ptmax = -1.; |
| 1911 | Int_t index = -1; |
| 1912 | Int_t jndex = -1; |
| 1913 | for (seed = 0; seed < nseed; seed++) { |
| 1914 | if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) { |
| 1915 | ptmax = pt; |
| 1916 | index = seed; |
| 1917 | } // ptmax ? |
| 1918 | } // seeds |
| 1919 | if (ptmax < 0.) break; |
| 1920 | jndex = iSeeds[index]; |
| 1921 | // |
| 1922 | // Remove from the list |
| 1923 | iSeeds[index] = -1; |
| 1924 | printf("\n Next Seed %d %f", jndex, ptmax); |
| 1925 | // |
| 1926 | // Find tracks in cone around seed |
| 1927 | // |
| 1928 | Float_t phiSeed = fPhiT[jndex]; |
| 1929 | Float_t etaSeed = fEtaT[jndex]; |
| 1930 | Float_t eT = 0.; |
| 1931 | Float_t pxJ = 0.; |
| 1932 | Float_t pyJ = 0.; |
| 1933 | Float_t pzJ = 0.; |
| 1934 | |
| 1935 | for (part = 0; part < fNt; part++) { |
| 1936 | if (!fTrackList[part]) continue; |
| 1937 | Float_t deta = fEtaT[part] - etaSeed; |
| 1938 | Float_t dphi = fPhiT[part] - phiSeed; |
| 1939 | Float_t dR = TMath::Sqrt(deta * deta + dphi * dphi); |
| 1940 | if (dR < fConeRadius) { |
| 1941 | eT += fPtT[part]; |
| 1942 | Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part])); |
| 1943 | Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]); |
| 1944 | Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]); |
| 1945 | Float_t pz = fPtT[part] / TMath::Tan(theta); |
| 1946 | pxJ += px; |
| 1947 | pyJ += py; |
| 1948 | pzJ += pz; |
| 1949 | // |
| 1950 | // if seed, remove it |
| 1951 | // |
| 1952 | for (seed = 0; seed < nseed; seed++) { |
| 1953 | if (part == iSeeds[seed]) iSeeds[seed] = -1; |
| 1954 | } // seed ? |
| 1955 | } // < cone radius |
| 1956 | } // particle loop |
| 1957 | |
| 1958 | // |
| 1959 | // Estimate of jet direction |
| 1960 | Float_t phiJ = TMath::ATan2(pyJ, pxJ); |
| 1961 | Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ); |
| 1962 | Float_t etaJ = TMath::Log(TMath::Tan(thetaJ / 2.)); |
| 1963 | // Float_t ptJ = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ); |
| 1964 | |
| 1965 | // |
| 1966 | // Sum up all energy |
| 1967 | // |
| 1968 | Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi); |
| 1969 | Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta); |
| 1970 | Int_t dIphi = Int_t(fConeRadius / fDphi); |
| 1971 | Int_t dIeta = Int_t(fConeRadius / fDeta); |
| 1972 | Int_t iPhi, iEta; |
| 1973 | Float_t sumE = 0; |
| 1974 | for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) { |
| 1975 | for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) { |
| 1976 | if (iPhi < 0 || iEta < 0) continue; |
| 1977 | Float_t dPhi = fPhiMin + iPhi * fDphi; |
| 1978 | Float_t dEta = fEtaMin + iEta * fDeta; |
| 1979 | if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue; |
| 1980 | sumE += fLego->GetBinContent(iEta, iPhi); |
| 1981 | } // eta |
| 1982 | } // phi |
| 1983 | // |
| 1984 | // |
| 1985 | // |
| 1986 | fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ); |
| 1987 | FindTracksInJetCone(); |
| 1988 | printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets); |
| 1989 | printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]); |
| 1990 | printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]); |
| 1991 | } // while(1) |
| 1992 | EMCALJETS.njet = njets; |
| 1993 | if (fWrite) WriteJets(); |
| 1994 | fEvent++; |
| 1995 | } |
| 1996 | // 16-jan-2003 - just for convenience |
| 1997 | void AliEMCALJetFinder::Browse(TBrowser* b) |
| 1998 | { |
| 1999 | // Add to browser |
| 2000 | if(fHistsList) b->Add((TObject*)fHistsList); |
| 2001 | } |
| 2002 | |
| 2003 | Bool_t AliEMCALJetFinder::IsFolder() const |
| 2004 | { |
| 2005 | // Return folder status |
| 2006 | if(fHistsList) return kTRUE; |
| 2007 | else return kFALSE; |
| 2008 | } |
| 2009 | |
| 2010 | const Char_t* AliEMCALJetFinder::GetNameOfVariant() |
| 2011 | { |
| 2012 | // generate the literal string with info about jet finder |
| 2013 | Char_t name[200]; |
| 2014 | sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f", |
| 2015 | fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF); |
| 2016 | TString nt(name); |
| 2017 | nt.ReplaceAll(" ",""); |
| 2018 | if(fBGFileName.Length()) { |
| 2019 | Int_t i1 = fBGFileName.Index("kBackground"); |
| 2020 | Int_t i2 = fBGFileName.Index("/0000") - 1; |
| 2021 | if(i1>=0 && i2>=0) { |
| 2022 | TString bg(fBGFileName(i1,i2-i1+1)); |
| 2023 | nt += bg; |
| 2024 | } |
| 2025 | } |
| 2026 | printf("<I> Name of variant %s \n", nt.Data()); |
| 2027 | return nt.Data(); |
| 2028 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame