1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //---------------------------------------------------------------------
19 // UA1 Cone Algorithm Jet finder for jet studies
20 // manages the search for jets using charged particle momentum and
21 // neutral cell energy information
22 // Authors: Rafael.Diaz.Valdes@cern.ch
23 // magali.estienne@subatech.in2p3.fr
24 // alexandre.shabetai@cern.ch
26 // Modified accordingly to reader/finder splitting and new handling of neutral information
27 // Versions V1 and V2 merged
28 //---------------------------------------------------------------------
33 #include "AliUA1JetFinder.h"
34 #include "AliUA1JetHeaderV1.h"
35 #include "AliJetCalTrk.h"
36 #include "AliJetBkg.h"
37 #include "AliAODJetEventBackground.h"
38 #include "AliAODJet.h"
40 ClassImp(AliUA1JetFinder)
42 ////////////////////////////////////////////////////////////////////////
44 AliUA1JetFinder::AliUA1JetFinder():
47 fJetBkg(new AliJetBkg())
49 // Default constructor
52 //-----------------------------------------------------------------------
53 AliUA1JetFinder::~AliUA1JetFinder()
61 //-----------------------------------------------------------------------
62 void AliUA1JetFinder::FindJets()
64 // Used to find jets using charged particle momentum information
65 // & neutral energy from calo cells
67 // 1) Fill cell map array
68 // 2) calculate total energy and fluctuation level
70 // 3.1) look centroides in cell map
71 // 3.2) calculate total energy in cones
72 // 3.3) flag as a possible jet
73 // 3.4) reorder cones by energy
74 // 4) subtract backg in accepted jets
75 // 5) fill AliJet list
77 // transform input to pt,eta,phi plus lego
79 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
80 Int_t nIn = fCalTrkEvent->GetNCalTrkTracks();
81 fDebug = fHeader->GetDebug();
84 fJetBkg->SetHeader(fHeader);
85 fJetBkg->SetCalTrkEvent(GetCalTrkEvent());
86 fJetBkg->SetDebug(fDebug);
87 // local arrays for input
88 // ToDo: check memory fragmentation, maybe better to
89 // define them globally and resize as needed
90 // Fragmentation should be worse for low mult...
91 Float_t* ptT = new Float_t[nIn];
92 Float_t* etaT = new Float_t[nIn];
93 Float_t* phiT = new Float_t[nIn];
94 Int_t* injet = new Int_t[nIn];
95 Int_t* injetOk = new Int_t[nIn];
97 memset(ptT,0,sizeof(Float_t)*nIn);
98 memset(etaT,0,sizeof(Float_t)*nIn);
99 memset(phiT,0,sizeof(Float_t)*nIn);
100 memset(injet,0,sizeof(Int_t)*nIn);
101 memset(injetOk,-1,sizeof(Int_t)*nIn);
103 // load input vectors and calculate total energy in array
105 // total energy in array
106 Float_t etbgTotal = 0.;
110 for (Int_t i = 0; i < fCalTrkEvent->GetNCalTrkTracks(); i++){
111 ptT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
112 etaT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
113 phiT[i] = ((fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() < 0) ? (fCalTrkEvent->GetCalTrkTrack(i)->GetPhi()) + 2 * TMath::Pi() : fCalTrkEvent->GetCalTrkTrack(i)->GetPhi());
114 //fCalTrkEvent->GetCalTrkTrack(i)->Print(Form("%d",i));
115 if (fCalTrkEvent->GetCalTrkTrack(i)->GetCutFlag() != 1) continue;
116 fLego->Fill(etaT[i], phiT[i], ptT[i]);
119 etbg2 += ptT[i]*ptT[i];
122 // calculate total energy and fluctuation in map
123 Double_t meanpt = 0.;
126 meanpt = etbgTotal/npart;
128 if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
129 ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
132 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
134 // arrays to hold jets
135 Float_t etaJet[kMaxJets];
136 Float_t phiJet[kMaxJets];
137 Float_t etJet[kMaxJets];
138 Float_t etsigJet[kMaxJets]; //signal et in jet
139 Float_t etallJet[kMaxJets]; // total et in jet (tmp variable)
140 Int_t ncellsJet[kMaxJets];
141 Int_t multJetT[kMaxJets];
142 Int_t multJetC[kMaxJets];
143 Int_t multJet[kMaxJets];
144 Float_t *areaJet = new Float_t[kMaxJets];
145 // Used for jet reordering at the end of the jet finding procedure
146 Float_t etaJetOk[kMaxJets];
147 Float_t phiJetOk[kMaxJets];
148 Float_t etJetOk[kMaxJets];
149 Float_t etsigJetOk[kMaxJets]; //signal et in jet
150 Float_t etallJetOk[kMaxJets]; // total et in jet (tmp variable)
151 Int_t ncellsJetOk[kMaxJets];
152 Int_t multJetOk[kMaxJets];
153 Float_t *areaJetOk = new Float_t[kMaxJets];
154 Int_t nJets; // to hold number of jets found by algorithm
155 Int_t nj; // number of jets accepted
156 Float_t prec = header->GetPrecBg();
159 while(bgprec > prec){
160 //reset jet arrays in memory
161 memset(etaJet,0,sizeof(Float_t)*kMaxJets);
162 memset(phiJet,0,sizeof(Float_t)*kMaxJets);
163 memset(etJet,0,sizeof(Float_t)*kMaxJets);
164 memset(etallJet,0,sizeof(Float_t)*kMaxJets);
165 memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
166 memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
167 memset(multJetT,0,sizeof(Int_t)*kMaxJets);
168 memset(multJetC,0,sizeof(Int_t)*kMaxJets);
169 memset(multJet,0,sizeof(Int_t)*kMaxJets);
170 memset(areaJet,0,sizeof(Float_t)*kMaxJets);
171 memset(etaJetOk,0,sizeof(Float_t)*kMaxJets);
172 memset(phiJetOk,0,sizeof(Float_t)*kMaxJets);
173 memset(etJetOk,0,sizeof(Float_t)*kMaxJets);
174 memset(etallJetOk,0,sizeof(Float_t)*kMaxJets);
175 memset(etsigJetOk,0,sizeof(Float_t)*kMaxJets);
176 memset(ncellsJetOk,0,sizeof(Int_t)*kMaxJets);
177 memset(multJetOk,0,sizeof(Int_t)*kMaxJets);
178 memset(areaJetOk,0,sizeof(Float_t)*kMaxJets);
181 // reset particles-jet array in memory
182 memset(injet,-1,sizeof(Int_t)*nIn);
183 //run cone algorithm finder
184 RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
185 //run background subtraction
186 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
187 nj = header->GetNAcceptJets();
191 //subtract background
192 Float_t etbgTotalN = 0.0; //new background
193 Float_t sigmaN = 0.0; //new background
194 if(header->GetBackgMode() == 1) {// standard
195 fJetBkg->SubtractBackg(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);}
196 if(header->GetBackgMode() == 2) //cone
197 fJetBkg->SubtractBackgCone(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
198 if(header->GetBackgMode() == 3) //ratio
199 fJetBkg->SubtractBackgRatio(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
200 if(header->GetBackgMode() == 4) //statistic
201 fJetBkg->SubtractBackgStat(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
204 if(etbgTotalN != 0.0)
205 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
208 etbgTotal = etbgTotalN; // update with new background estimation
212 // add tracks to the jet if it wasn't yet done
213 if (header->GetBackgMode() == 0){
214 Float_t rc= header->GetRadius();
215 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
216 for(Int_t ijet=0; ijet<nJets; ijet++){
217 Float_t deta = etaT[jpart] - etaJet[ijet];
218 Float_t dphi = phiT[jpart] - phiJet[ijet];
219 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
220 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
221 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
222 if(dr <= rc){ // particles inside this cone
227 } //end particle loop
231 if (fDebug>1) printf("Found %d jets \n", nj);
233 // Reorder jets by et in cone
234 // Sort jets by energy
236 TMath::Sort(nJets, etJet, idx);
237 for(Int_t p = 0; p < nJets; p++)
239 etaJetOk[p] = etaJet[idx[p]];
240 phiJetOk[p] = phiJet[idx[p]];
241 etJetOk[p] = etJet[idx[p]];
242 etallJetOk[p] = etJet[idx[p]];
243 etsigJetOk[p] = etsigJet[idx[p]];
244 ncellsJetOk[p] = ncellsJet[idx[p]];
245 multJetOk[p] = multJet[idx[p]];
246 areaJetOk[p] = areaJet[idx[p]];
249 //////////////////////////
251 Int_t nTracks = fCalTrkEvent->GetNCalTrkTracks();
253 for(Int_t kj=0; kj<nj; kj++)
255 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
256 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
257 (phiJetOk[kj] > (header->GetJetPhiMax())) ||
258 (phiJetOk[kj] < (header->GetJetPhiMin())) ||
259 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
260 Float_t px=-999, py=-999 ,pz=-999 ,en=-999; // convert to 4-vector
261 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
262 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
263 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
264 en = TMath::Sqrt(px * px + py * py + pz * pz);
265 AliAODJet jet(px, py, pz, en);
267 // Calc jet area if it wasn't yet done
268 if (header->GetBackgMode() == 0){
269 // calculate the area of the jet
270 Float_t rc= header->GetRadius();
271 areaJetOk[kj] = fJetBkg->CalcJetAreaEtaCut(rc,etaJetOk[kj]);
274 jet.SetEffArea(areaJetOk[kj],0.,0.,0.);
275 jet.SetBgEnergy(etbgTotal,0.);
276 if (fDebug>1) jet.Print(Form("%d",kj));
278 for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
279 // Track to jet reordering
280 if(injet[jpart] == idx[kj]){
283 // Check if the particle belongs to the jet and add the ref
284 if(injetOk[jpart] == kj && fCalTrkEvent->GetCalTrkTrack(jpart)->GetCutFlag() == 1) {
285 jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jpart)->GetTrackObject());
304 //-----------------------------------------------------------------------
305 void AliUA1JetFinder::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
306 Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
307 Float_t* const etallJet, Int_t* const ncellsJet)
310 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
311 const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
313 const Int_t nBinEta = header->GetLegoNbinEta();
314 const Int_t nBinPhi = header->GetLegoNbinPhi();
315 if((nBinPhi*nBinEta)>nBinsMax){
316 AliError("Too many bins of the ETA-PHI histogram");
319 Float_t etCell[nBinsMax] = {0.}; //! Cell Energy
320 Float_t etaCell[nBinsMax] = {0.}; //! Cell eta
321 Float_t phiCell[nBinsMax] = {0.}; //! Cell phi
322 Short_t flagCell[nBinsMax] = {0}; //! Cell flag
325 TAxis* xaxis = fLego->GetXaxis();
326 TAxis* yaxis = fLego->GetYaxis();
328 for (Int_t i = 1; i <= nBinEta; i++) {
329 for (Int_t j = 1; j <= nBinPhi; j++) {
330 e = fLego->GetBinContent(i,j);
331 if (e < 0.0) continue; // don't include this cells
332 Float_t eta = xaxis->GetBinCenter(i);
333 Float_t phi = yaxis->GetBinCenter(j);
335 etaCell[nCell] = eta;
336 phiCell[nCell] = phi;
337 flagCell[nCell] = 0; //default
341 // Parameters from header
342 Float_t minmove = header->GetMinMove();
343 Float_t maxmove = header->GetMaxMove();
344 Float_t rc = header->GetRadius();
345 Float_t etseed = header->GetEtSeed();
346 // Tmp array of jets form algoritm
347 Float_t etaAlgoJet[kMaxJets] = {0.0};
348 Float_t phiAlgoJet[kMaxJets] = {0.0};
349 Float_t etAlgoJet[kMaxJets] = {0.0};
350 Int_t ncellsAlgoJet[kMaxJets] = {0};
355 Int_t index[nBinsMax];
356 TMath::Sort(nCell, etCell, index);
357 // variable used in centroide loop
375 for(Int_t icell = 0; icell < nCell; icell++)
377 Int_t jcell = index[icell];
378 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
379 if(flagCell[jcell] != 0) continue; // if cell was used before
381 eta = etaCell[jcell];
382 phi = phiCell[jcell];
393 for(Int_t kcell =0; kcell < nCell; kcell++)
395 Int_t lcell = index[kcell];
396 if(lcell == jcell) continue; // cell itself
397 if(flagCell[lcell] != 0) continue; // cell used before
398 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
400 deta = etaCell[lcell] - eta;
401 dphi = TMath::Abs(phiCell[lcell] - phi);
402 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
403 dr = TMath::Sqrt(deta * deta + dphi * dphi);
406 // calculate offset from initiate cell
407 deta = etaCell[lcell] - eta0;
408 dphi = phiCell[lcell] - phi0;
409 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
410 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
411 etas = etas + etCell[lcell]*deta;
412 phis = phis + etCell[lcell]*dphi;
413 ets = ets + etCell[lcell];
414 //new weighted eta and phi including this cell
415 eta = eta0 + etas/ets;
416 phi = phi0 + phis/ets;
417 // if cone does not move much, just go to next step
418 dphib = TMath::Abs(phi - phib);
419 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
420 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
421 if(dr <= minmove) break;
422 // cone should not move more than max_mov
423 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
430 } else { // store this loop information
438 }//end of cells loop looking centroide
440 // Avoid cones overloap (to be implemented in the future)
442 // Flag cells in Rc, estimate total energy in cone
443 Float_t etCone = 0.0;
445 rc = header->GetRadius();
447 for(Int_t ncell =0; ncell < nCell; ncell++)
449 if(flagCell[ncell] != 0) continue; // cell used before
451 deta = etaCell[ncell] - eta;
452 dphi = phiCell[ncell] - phi;
453 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
454 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
455 dr = TMath::Sqrt(deta * deta + dphi * dphi);
456 if(dr <= rc){ // cell in cone
457 flagCell[ncell] = -1;
458 etCone+=etCell[ncell];
463 // Select jets with et > background
464 // estimate max fluctuation of background in cone
465 Double_t ncellin = (Double_t)nCellIn;
466 Double_t ntcell = (Double_t)nCell;
467 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
469 Double_t etcmin = etCone ; // could be used etCone - etmin !!
470 //decisions !! etbmax < etcmin
472 for(Int_t mcell =0; mcell < nCell; mcell++){
473 if(flagCell[mcell] == -1){
475 flagCell[mcell] = 1; //flag cell as used
477 flagCell[mcell] = 0; // leave it free
480 //store tmp jet info !!!
481 if(etbmax < etcmin) {
483 etaAlgoJet[nJets] = eta;
484 phiAlgoJet[nJets] = phi;
485 etAlgoJet[nJets] = etCone;
486 ncellsAlgoJet[nJets] = nCellIn;
490 AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
494 } // end of cells loop
496 //reorder jets by et in cone
497 //sort jets by energy
498 for(Int_t p = 0; p < nJets; p++)
500 etaJet[p] = etaAlgoJet[p];
501 phiJet[p] = phiAlgoJet[p];
502 etJet[p] = etAlgoJet[p];
503 etallJet[p] = etAlgoJet[p];
504 ncellsJet[p] = ncellsAlgoJet[p];
509 //-----------------------------------------------------------------------
510 void AliUA1JetFinder::Reset()
513 AliJetFinder::Reset();
517 //-----------------------------------------------------------------------
518 void AliUA1JetFinder::WriteJHeaderToFile() const
520 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
525 //-----------------------------------------------------------------------
526 void AliUA1JetFinder::Init()
529 // initializes some variables
530 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
532 fLego = new TH2F("legoH","eta-phi",
533 header->GetLegoNbinEta(), header->GetLegoEtaMin(),
534 header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
535 header->GetLegoPhiMin(), header->GetLegoPhiMax());