]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliUA1JetFinderV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
undo previous changes
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinderV1.cxx
CommitLineData
70e58892 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 **************************************************************************/
98e98c1c 15
7ca4655f 16/* $Id$ */
70e58892 17
18//---------------------------------------------------------------------
19// UA1 Cone Algorithm Jet finder
20// manages the search for jets
21// Author: Rafael.Diaz.Valdes@cern.ch
22// (version in c++)
23//---------------------------------------------------------------------
24
7ca4655f 25#include <TArrayF.h>
26#include <TClonesArray.h>
70e58892 27#include <TFile.h>
28#include <TH1F.h>
29#include <TH2F.h>
7ca4655f 30#include <TLorentzVector.h>
31
70e58892 32#include "AliUA1JetFinderV1.h"
33#include "AliUA1JetHeaderV1.h"
34#include "AliJetReaderHeader.h"
35#include "AliJetReader.h"
952b368c 36#include "AliJetHeader.h"
37
38
1d27ecd2 39#include "AliAODJet.h"
d9eba389 40#include "AliLog.h"
70e58892 41
42
43ClassImp(AliUA1JetFinderV1)
44
9e4cc50d 45/////////////////////////////////////////////////////////////////////
1b7d5d7e 46
9e4cc50d 47AliUA1JetFinderV1::AliUA1JetFinderV1() :
48 AliJetFinder(),
5e380e7e 49 fLego(0),
50 fhEtBackg(0),
51 fhAreaBackg(0)
98e98c1c 52{
53 // Constructor
5e380e7e 54 for(int i = 0;i < kMaxJets;i++){
55 fhAreaJet[i] = fhEtJet[i] = 0;
56 }
70e58892 57}
58
59////////////////////////////////////////////////////////////////////////
60
61AliUA1JetFinderV1::~AliUA1JetFinderV1()
62
63{
64 // destructor
c1f55a27 65 delete fLego;
66 fLego = 0;
5e380e7e 67 if(fhEtBackg)delete fhEtBackg;
68 fhEtBackg = 0;
69 if( fhAreaBackg) delete fhAreaBackg;
70 fhAreaBackg = 0;
71 for(int i = 0;i < kMaxJets;i++){
72 if(fhAreaJet[i])delete fhAreaJet[i];
73 if(fhEtJet[i]) delete fhEtJet[i];
74 fhAreaJet[i] = fhEtJet[i] = 0;
75 }
76
70e58892 77}
78
79////////////////////////////////////////////////////////////////////////
80
81
82void AliUA1JetFinderV1::FindJets()
83
84{
85 //1) Fill cell map array
86 //2) calculate total energy and fluctuation level
87 //3) Run algorithm
88 // 3.1) look centroides in cell map
89 // 3.2) calculate total energy in cones
90 // 3.3) flag as a possible jet
91 // 3.4) reorder cones by energy
92 //4) subtract backg in accepted jets
93 //5) fill AliJet list
94
95 // transform input to pt,eta,phi plus lego
7d0f353c 96
97 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
70e58892 98 TClonesArray *lvArray = fReader->GetMomentumArray();
99 Int_t nIn = lvArray->GetEntries();
100 if (nIn == 0) return;
101
102 // local arrays for input
5e380e7e 103 // ToDo: check memory fragmentation, maybe better to
104 // define them globally and resize as needed
105 // Fragementation should be worse for low mult...
1d27ecd2 106 Float_t* ptT = new Float_t[nIn];
107 Float_t* etaT = new Float_t[nIn];
108 Float_t* phiT = new Float_t[nIn];
70e58892 109 Int_t* injet = new Int_t[nIn];
110
5e380e7e 111
112
113 // load input vectors and calculate total energy in array
114
70e58892 115 //total energy in array
116 Float_t etbgTotal = 0.0;
5e380e7e 117 Float_t npart = 0;
118 Float_t etbg2 = 0;
70e58892 119
70e58892 120 for (Int_t i = 0; i < nIn; i++){
121 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
122 ptT[i] = lv->Pt();
123 etaT[i] = lv->Eta();
124 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
125 if (fReader->GetCutFlag(i) != 1) continue;
9dda5307 126 fLego ->Fill(etaT[i], phiT[i], ptT[i]);
5e380e7e 127 npart += 1;
70e58892 128 etbgTotal+= ptT[i];
5e380e7e 129 etbg2 += ptT[i]*ptT[i];
70e58892 130 }
9dda5307 131
70e58892 132 // calculate total energy and fluctuation in map
5e380e7e 133 Double_t meanpt = 0;
134 Double_t ptRMS = 0;
135 if(npart>0){
136 meanpt = etbgTotal/npart;
137 etbg2 = etbg2/npart;
138 if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
139 ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
140 }
141 }
70e58892 142 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
143
144 // arrays to hold jets
5382c085 145 Float_t etaJet[kMaxJets];
146 Float_t phiJet[kMaxJets];
147 Float_t etJet[kMaxJets];
148 Float_t etsigJet[kMaxJets]; //signal et in jet
149 Float_t etallJet[kMaxJets]; // total et in jet (tmp variable)
150 Int_t ncellsJet[kMaxJets];
151 Int_t multJet[kMaxJets];
70e58892 152 Int_t nJets; // to hold number of jets found by algorithm
153 Int_t nj; // number of jets accepted
7d0f353c 154 Float_t prec = header->GetPrecBg();
70e58892 155 Float_t bgprec = 1;
156 while(bgprec > prec){
157 //reset jet arrays in memory
5382c085 158 memset(etaJet,0,sizeof(Float_t)*kMaxJets);
159 memset(phiJet,0,sizeof(Float_t)*kMaxJets);
160 memset(etJet,0,sizeof(Float_t)*kMaxJets);
161 memset(etallJet,0,sizeof(Float_t)*kMaxJets);
162 memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
163 memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
164 memset(multJet,0,sizeof(Int_t)*kMaxJets);
70e58892 165 nJets = 0;
166 nj = 0;
167 // reset particles-jet array in memory
98e98c1c 168 memset(injet,-1,sizeof(Int_t)*nIn);
70e58892 169 //run cone algorithm finder
170 RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
171 //run background subtraction
7d0f353c 172 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
173 nj = header->GetNAcceptJets();
70e58892 174 else
175 nj = nJets;
176 //subtract background
177 Float_t etbgTotalN = 0.0; //new background
7d0f353c 178 if(header->GetBackgMode() == 1) // standar
70e58892 179 SubtractBackg(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
7d0f353c 180 if(header->GetBackgMode() == 2) //cone
70e58892 181 SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
7d0f353c 182 if(header->GetBackgMode() == 3) //ratio
70e58892 183 SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
7d0f353c 184 if(header->GetBackgMode() == 4) //statistic
70e58892 185 SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
186 //calc precision
187 if(etbgTotalN != 0.0)
188 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
189 else
190 bgprec = 0;
191 etbgTotal = etbgTotalN; // update with new background estimation
192 } //end while
193
f1744ef9 194 // add tracks to the jet if it wasn't yet done
195 if (header->GetBackgMode() == 0){
196 Float_t rc= header->GetRadius();
197 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
198 for(Int_t ijet=0; ijet<nj; ijet++){
199 Float_t deta = etaT[jpart] - etaJet[ijet];
200 Float_t dphi = phiT[jpart] - phiJet[ijet];
201 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
202 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
203 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
204 if(dr <= rc){ // particles inside this cone
205 injet[jpart] = ijet;
206 break;
207 }
208 }// end jets loop
209 } //end particle loop
210 }
211
70e58892 212 // add jets to list
5e380e7e 213 Int_t idxjets[kMaxJets];
70e58892 214 Int_t nselectj = 0;
5e380e7e 215
c02b94ea 216 TRefArray *refs = 0;
217 Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
218 if (fromAod) refs = fReader->GetReferences();
70e58892 219 for(Int_t kj=0; kj<nj; kj++){
7d0f353c 220 if ((etaJet[kj] > (header->GetJetEtaMax())) ||
221 (etaJet[kj] < (header->GetJetEtaMin())) ||
222 (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
70e58892 223 Float_t px, py,pz,en; // convert to 4-vector
224 px = etJet[kj] * TMath::Cos(phiJet[kj]);
225 py = etJet[kj] * TMath::Sin(phiJet[kj]);
226 pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
227 en = TMath::Sqrt(px * px + py * py + pz * pz);
42b0ac89 228
1d27ecd2 229 AliAODJet jet(px, py, pz, en);
1eef1eb4 230
c02b94ea 231 if (fromAod){
232 for(Int_t jpart = 0; jpart < nIn; jpart++) // loop for all particles in array
233 if (injet[jpart] == kj && fReader->GetCutFlag(jpart) == 1)
234 jet.AddTrack(refs->At(jpart)); // check if the particle belongs to the jet and add the ref
235 }
1eef1eb4 236
99513b18 237 //jet.Print("");
1d27ecd2 238
239 AddJet(jet);
240
70e58892 241 idxjets[nselectj] = kj;
242 nselectj++;
1eef1eb4 243 } //end particle loop
244
70e58892 245 //add signal percentage and total signal in AliJets for analysis tool
5e380e7e 246 Float_t percentage[kMaxJets];
247 Int_t ncells[kMaxJets];
248 Int_t mult[kMaxJets];
70e58892 249 for(Int_t i = 0; i< nselectj; i++){
250 percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
251 ncells[i] = ncellsJet[idxjets[i]];
252 mult[i] = multJet[idxjets[i]];
253 }
98e98c1c 254 //add particle-injet relationship ///
255 for(Int_t bj = 0; bj < nIn; bj++){
256 if(injet[bj] == -1) continue; //background particle
257 Int_t bflag = 0;
258 for(Int_t ci = 0; ci< nselectj; ci++){
259 if(injet[bj] == idxjets[ci]){
260 injet[bj]= ci;
261 bflag++;
262 break;
263 }
264 }
265 if(bflag == 0) injet[bj] = -1; // set as background particle
266 }
70e58892 267
268 //delete
0ab34e91 269 delete [] ptT;
270 delete [] etaT;
271 delete [] phiT;
272 delete [] injet;
70e58892 273}
274
275////////////////////////////////////////////////////////////////////////
276
277void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
278 Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
279 Float_t* etallJet, Int_t* ncellsJet)
280{
281
282 //dump lego
7d0f353c 283 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
d9eba389 284 const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
9dda5307 285
d9eba389 286 const Int_t nBinEta = header->GetLegoNbinEta();
287 const Int_t nBinPhi = header->GetLegoNbinPhi();
288 if((nBinPhi*nBinEta)>nBinsMax){
289 AliError("Too many bins of the ETA-PHI histogram");
290 }
291
0ab34e91 292 Float_t etCell[nBinsMax]; //! Cell Energy
293 Float_t etaCell[nBinsMax]; //! Cell eta
294 Float_t phiCell[nBinsMax]; //! Cell phi
d9eba389 295 Short_t flagCell[nBinsMax]; //! Cell flag
70e58892 296
297 Int_t nCell = 0;
298 TAxis* xaxis = fLego->GetXaxis();
299 TAxis* yaxis = fLego->GetYaxis();
300 Float_t e = 0.0;
d9eba389 301 for (Int_t i = 1; i <= nBinEta; i++) {
302 for (Int_t j = 1; j <= nBinPhi; j++) {
70e58892 303 e = fLego->GetBinContent(i,j);
304 if (e < 0.0) continue; // don't include this cells
305 Float_t eta = xaxis->GetBinCenter(i);
306 Float_t phi = yaxis->GetBinCenter(j);
307 etCell[nCell] = e;
308 etaCell[nCell] = eta;
309 phiCell[nCell] = phi;
0ab34e91 310 flagCell[nCell] = 0; //default
70e58892 311 nCell++;
312 }
313 }
314
315 // Parameters from header
7d0f353c 316 Float_t minmove = header->GetMinMove();
317 Float_t maxmove = header->GetMaxMove();
318 Float_t rc = header->GetRadius();
319 Float_t etseed = header->GetEtSeed();
320 //Float_t etmin = header->GetMinJetEt();
70e58892 321
322
323
324 // tmp array of jets form algoritm
5e380e7e 325 Float_t etaAlgoJet[kMaxJets];
326 Float_t phiAlgoJet[kMaxJets];
327 Float_t etAlgoJet[kMaxJets];
328 Int_t ncellsAlgoJet[kMaxJets];
70e58892 329
330 //run algorithm//
331
332 // sort cells by et
5e380e7e 333 Int_t index[nBinsMax];
70e58892 334 TMath::Sort(nCell, etCell, index);
335 // variable used in centroide loop
9dda5307 336 Float_t eta = 0.0;
337 Float_t phi = 0.0;
338 Float_t eta0 = 0.0;
339 Float_t phi0 = 0.0;
340 Float_t etab = 0.0;
341 Float_t phib = 0.0;
342 Float_t etas = 0.0;
343 Float_t phis = 0.0;
344 Float_t ets = 0.0;
345 Float_t deta = 0.0;
346 Float_t dphi = 0.0;
347 Float_t dr = 0.0;
348 Float_t etsb = 0.0;
70e58892 349 Float_t etasb = 0.0;
350 Float_t phisb = 0.0;
9dda5307 351 Float_t dphib = 0.0;
352
70e58892 353
354 for(Int_t icell = 0; icell < nCell; icell++){
355 Int_t jcell = index[icell];
356 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
357 if(flagCell[jcell] != 0) continue; // if cell was used before
358 eta = etaCell[jcell];
359 phi = phiCell[jcell];
360 eta0 = eta;
361 phi0 = phi;
362 etab = eta;
363 phib = phi;
364 ets = etCell[jcell];
365 etas = 0.0;
366 phis = 0.0;
367 etsb = ets;
368 etasb = 0.0;
369 phisb = 0.0;
370 for(Int_t kcell =0; kcell < nCell; kcell++){
371 Int_t lcell = index[kcell];
9dda5307 372 if(lcell == jcell) continue; // cell itself
373 if(flagCell[lcell] != 0) continue; // cell used before
374 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
70e58892 375 //calculate dr
376 deta = etaCell[lcell] - eta;
9dda5307 377 dphi = TMath::Abs(phiCell[lcell] - phi);
378 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
379 dr = TMath::Sqrt(deta * deta + dphi * dphi);
70e58892 380 if(dr <= rc){
381 // calculate offset from initiate cell
382 deta = etaCell[lcell] - eta0;
383 dphi = phiCell[lcell] - phi0;
9dda5307 384 if (dphi < - TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
385 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
386
70e58892 387 etas = etas + etCell[lcell]*deta;
388 phis = phis + etCell[lcell]*dphi;
389 ets = ets + etCell[lcell];
390 //new weighted eta and phi including this cell
391 eta = eta0 + etas/ets;
392 phi = phi0 + phis/ets;
393 // if cone does not move much, just go to next step
9dda5307 394 dphib = TMath::Abs(phi - phib);
395 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
396 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
70e58892 397 if(dr <= minmove) break;
398 // cone should not move more than max_mov
399 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
400 if(dr > maxmove){
9dda5307 401 eta = etab;
402 phi = phib;
403 ets = etsb;
404 etas = etasb;
405 phis = phisb;
406 } else { // store this loop information
407 etab = eta;
408 phib = phi;
409 etsb = ets;
410 etasb = etas;
411 phisb = phis;
70e58892 412 }
9dda5307 413 } // inside cone
70e58892 414 }//end of cells loop looking centroide
415
416 //avoid cones overloap (to be implemented in the future)
417
418 //flag cells in Rc, estimate total energy in cone
419 Float_t etCone = 0.0;
420 Int_t nCellIn = 0;
7d0f353c 421 rc = header->GetRadius();
70e58892 422 for(Int_t ncell =0; ncell < nCell; ncell++){
423 if(flagCell[ncell] != 0) continue; // cell used before
424 //calculate dr
425 deta = etaCell[ncell] - eta;
426 dphi = phiCell[ncell] - phi;
427 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
428 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
429 dr = TMath::Sqrt(deta * deta + dphi * dphi);
430 if(dr <= rc){ // cell in cone
431 flagCell[ncell] = -1;
432 etCone+=etCell[ncell];
433 nCellIn++;
434 }
435 }
436
437 // select jets with et > background
438 // estimate max fluctuation of background in cone
439 Double_t ncellin = (Double_t)nCellIn;
440 Double_t ntcell = (Double_t)nCell;
441 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
442 // min cone et
443 Double_t etcmin = etCone ; // could be used etCone - etmin !!
444 //desicions !! etbmax < etcmin
445 for(Int_t mcell =0; mcell < nCell; mcell++){
446 if(flagCell[mcell] == -1){
447 if(etbmax < etcmin)
448 flagCell[mcell] = 1; //flag cell as used
449 else
450 flagCell[mcell] = 0; // leave it free
451 }
452 }
453 //store tmp jet info !!!
98e98c1c 454 if(etbmax < etcmin) {
5e380e7e 455 if(nJets<kMaxJets){
456 etaAlgoJet[nJets] = eta;
457 phiAlgoJet[nJets] = phi;
458 etAlgoJet[nJets] = etCone;
459 ncellsAlgoJet[nJets] = nCellIn;
460 nJets++;
461 }
462 else{
463 AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
464 break;
465 }
466 }
70e58892 467 } // end of cells loop
468
469 //reorder jets by et in cone
470 //sort jets by energy
5e380e7e 471 Int_t idx[kMaxJets];
472 TMath::Sort(nJets, etAlgoJet, idx); // sort only the found jets
70e58892 473 for(Int_t p = 0; p < nJets; p++){
474 etaJet[p] = etaAlgoJet[idx[p]];
475 phiJet[p] = phiAlgoJet[idx[p]];
476 etJet[p] = etAlgoJet[idx[p]];
477 etallJet[p] = etAlgoJet[idx[p]];
478 ncellsJet[p] = ncellsAlgoJet[idx[p]];
479 }
480
70e58892 481}
482////////////////////////////////////////////////////////////////////////
483
63cdee97 484void AliUA1JetFinderV1::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
485 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
486 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
487 Int_t* multJet, Int_t* injet)
70e58892 488{
489 //background subtraction using cone method but without correction in dE/deta distribution
490
491 //calculate energy inside and outside cones
7d0f353c 492 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
493 Float_t rc= header->GetRadius();
5e380e7e 494 Float_t etIn[kMaxJets];
70e58892 495 Float_t etOut = 0;
496 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
98e98c1c 497 // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
70e58892 498 for(Int_t ijet=0; ijet<nJ; ijet++){
499 Float_t deta = etaT[jpart] - etaJet[ijet];
500 Float_t dphi = phiT[jpart] - phiJet[ijet];
501 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
502 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
503 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
504 if(dr <= rc){ // particles inside this cone
70e58892 505 multJet[ijet]++;
506 injet[jpart] = ijet;
98e98c1c 507 if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
508 etIn[ijet] += ptT[jpart];
509 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
510 }
70e58892 511 break;
512 }
513 }// end jets loop
970a3bbc 514 if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
98e98c1c 515 etOut += ptT[jpart]; // particle outside cones and pt cut
70e58892 516 } //end particle loop
517
518 //estimate jets and background areas
5e380e7e 519 Float_t areaJet[kMaxJets];
7d0f353c 520 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
70e58892 521 for(Int_t k=0; k<nJ; k++){
522 Float_t detamax = etaJet[k] + rc;
523 Float_t detamin = etaJet[k] - rc;
524 Float_t accmax = 0.0; Float_t accmin = 0.0;
7d0f353c 525 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
526 Float_t h = header->GetLegoEtaMax() - etaJet[k];
70e58892 527 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
528 }
7d0f353c 529 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
530 Float_t h = header->GetLegoEtaMax() + etaJet[k];
70e58892 531 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
532 }
533 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
534 areaOut = areaOut - areaJet[k];
535 }
536 //subtract background using area method
537 for(Int_t ljet=0; ljet<nJ; ljet++){
538 Float_t areaRatio = areaJet[ljet]/areaOut;
539 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
540 }
541
542 // estimate new total background
7d0f353c 543 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
70e58892 544 etbgTotalN = etOut*areaT/areaOut;
545
546
547}
548
549////////////////////////////////////////////////////////////////////////
550
63cdee97 551void AliUA1JetFinderV1::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
552 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
553 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
554 Int_t* multJet, Int_t* injet)
70e58892 555{
556
557 //background subtraction using statistical method
7d0f353c 558 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
559 Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
70e58892 560
561 //calculate energy inside
7d0f353c 562 Float_t rc= header->GetRadius();
5e380e7e 563 Float_t etIn[kMaxJets];
70e58892 564
565 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
98e98c1c 566 //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
70e58892 567 for(Int_t ijet=0; ijet<nJ; ijet++){
568 Float_t deta = etaT[jpart] - etaJet[ijet];
569 Float_t dphi = phiT[jpart] - phiJet[ijet];
570 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
571 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
572 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
573 if(dr <= rc){ // particles inside this cone
70e58892 574 multJet[ijet]++;
575 injet[jpart] = ijet;
98e98c1c 576 if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
577 etIn[ijet]+= ptT[jpart];
578 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
579 }
70e58892 580 break;
581 }
582 }// end jets loop
583 } //end particle loop
584
585 //calc jets areas
5e380e7e 586 Float_t areaJet[kMaxJets];
7d0f353c 587 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
70e58892 588 for(Int_t k=0; k<nJ; k++){
589 Float_t detamax = etaJet[k] + rc;
590 Float_t detamin = etaJet[k] - rc;
591 Float_t accmax = 0.0; Float_t accmin = 0.0;
7d0f353c 592 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
593 Float_t h = header->GetLegoEtaMax() - etaJet[k];
70e58892 594 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
595 }
7d0f353c 596 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
597 Float_t h = header->GetLegoEtaMax() + etaJet[k];
70e58892 598 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
599 }
600 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
601 }
602
603 //subtract background using area method
604 for(Int_t ljet=0; ljet<nJ; ljet++){
605 Float_t areaRatio = areaJet[ljet]/areaOut;
606 etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
607 }
608
609 etbgTotalN = etbgStat;
610
611}
612
613////////////////////////////////////////////////////////////////////////
614
63cdee97 615void AliUA1JetFinderV1::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
616 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
617 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
618 Int_t* multJet, Int_t* injet)
70e58892 619{
620 // Cone background subtraction method taking into acount dEt/deta distribution
7d0f353c 621 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
70e58892 622 //general
7d0f353c 623 Float_t rc= header->GetRadius();
624 Float_t etamax = header->GetLegoEtaMax();
625 Float_t etamin = header->GetLegoEtaMin();
70e58892 626 Int_t ndiv = 100;
627
628 // jet energy and area arrays
5e380e7e 629 char hEtname[256];char hAreaname[256];
630
70e58892 631 for(Int_t mjet=0; mjet<nJ; mjet++){
5e380e7e 632 if(!fhEtJet[mjet]){
633 sprintf(hEtname, "hEtJet%d", mjet);
634 fhEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
635 }
636 if(!fhAreaJet[mjet]){
637 sprintf(hAreaname, "hAreaJet%d", mjet);
638 fhAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
639 }
640 fhEtJet[mjet]->Reset();
641 fhAreaJet[mjet]->Reset();
642 }
70e58892 643 // background energy and area
5e380e7e 644 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
645 fhEtBackg->Reset();
646 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
647 fhAreaBackg->Reset();
70e58892 648
649 //fill energies
650 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
70e58892 651 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
652 Float_t deta = etaT[jpart] - etaJet[ijet];
98e98c1c 653 Float_t dphi = phiT[jpart] - phiJet[ijet];
70e58892 654 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
98e98c1c 655 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
656 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
70e58892 657 if(dr <= rc){ // particles inside this cone
70e58892 658 injet[jpart] = ijet;
98e98c1c 659 multJet[ijet]++;
660 if((fReader->GetCutFlag(jpart)) == 1){// pt cut
5e380e7e 661 fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
98e98c1c 662 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
663 }
70e58892 664 break;
665 }
666 }// end jets loop
970a3bbc 667 if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
5e380e7e 668 fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
70e58892 669 } //end particle loop
670
671 //calc areas
672 Float_t eta0 = etamin;
673 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
674 Float_t eta1 = eta0 + etaw;
675 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
676 Float_t etac = eta0 + etaw/2.0;
677 Float_t areabg = etaw*2.0*TMath::Pi();
678 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
679 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
680 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
681 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
682 Float_t areaj = 0.0;
683 if(deta0 > rc && deta1 < rc){
684 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
685 areaj = acc1;
686 }
687 if(deta0 < rc && deta1 > rc){
688 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
689 areaj = acc0;
690 }
691 if(deta0 < rc && deta1 < rc){
692 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
693 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
694 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
695 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
696 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
697 }
5e380e7e 698 fhAreaJet[ijet]->Fill(etac,areaj);
70e58892 699 areabg = areabg - areaj;
700 } // end jets loop
5e380e7e 701 fhAreaBackg->Fill(etac,areabg);
70e58892 702 eta0 = eta1;
703 eta1 = eta1 + etaw;
704 } // end loop for all eta bins
705
706 //subtract background
707 for(Int_t kjet=0; kjet<nJ; kjet++){
708 etJet[kjet] = 0.0; // first clear etJet for this jet
709 for(Int_t bin = 0; bin< ndiv; bin++){
5e380e7e 710 if(fhAreaJet[kjet]->GetBinContent(bin)){
711 Float_t areab = fhAreaBackg->GetBinContent(bin);
712 Float_t etb = fhEtBackg->GetBinContent(bin);
713 Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
714 etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
70e58892 715 }
716 }
717 }
718
719 // calc background total
5e380e7e 720 Double_t etOut = fhEtBackg->Integral();
721 Double_t areaOut = fhAreaBackg->Integral();
7d0f353c 722 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
70e58892 723 etbgTotalN = etOut*areaT/areaOut;
70e58892 724}
725
726////////////////////////////////////////////////////////////////////////
727
728
63cdee97 729void AliUA1JetFinderV1::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ, Float_t& etbgTotalN,
730 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
731 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
732 Int_t* multJet, Int_t* injet)
70e58892 733{
734 // Ratio background subtraction method taking into acount dEt/deta distribution
7d0f353c 735 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
70e58892 736 //factor F calc before
7d0f353c 737 Float_t bgRatioCut = header->GetBackgCutRatio();
70e58892 738
739
740 //general
7d0f353c 741 Float_t rc= header->GetRadius();
742 Float_t etamax = header->GetLegoEtaMax();
743 Float_t etamin = header->GetLegoEtaMin();
70e58892 744 Int_t ndiv = 100;
745
746 // jet energy and area arrays
5e380e7e 747 // jet energy and area arrays
748 char hEtname[256];char hAreaname[256];
749
70e58892 750 for(Int_t mjet=0; mjet<nJ; mjet++){
5e380e7e 751 if(!fhEtJet[mjet]){
752 sprintf(hEtname, "hEtJet%d", mjet);
753 fhEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
754 }
755 if(!fhAreaJet[mjet]){
756 sprintf(hAreaname, "hAreaJet%d", mjet);
757 fhAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
758 }
759 fhEtJet[mjet]->Reset();
760 fhAreaJet[mjet]->Reset();
761 }
70e58892 762 // background energy and area
5e380e7e 763 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
764 fhEtBackg->Reset();
765 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
766 fhAreaBackg->Reset();
70e58892 767
768 //fill energies
769 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
98e98c1c 770 //if((fReader->GetCutFlag(jpart)) != 1) continue;
70e58892 771 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
772 Float_t deta = etaT[jpart] - etaJet[ijet];
773 Float_t dphi = phiT[jpart] - phiJet[ijet];
774 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
775 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
776 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
777 if(dr <= rc){ // particles inside this cone
70e58892 778 multJet[ijet]++;
779 injet[jpart] = ijet;
98e98c1c 780 if((fReader->GetCutFlag(jpart)) == 1){ //pt cut
5e380e7e 781 fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
98e98c1c 782 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
783 }
70e58892 784 break;
785 }
786 }// end jets loop
5e380e7e 787 if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
70e58892 788 } //end particle loop
789
790 //calc areas
791 Float_t eta0 = etamin;
792 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
793 Float_t eta1 = eta0 + etaw;
794 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
795 Float_t etac = eta0 + etaw/2.0;
796 Float_t areabg = etaw*2.0*TMath::Pi();
797 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
798 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
799 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
800 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
801 Float_t areaj = 0.0;
802 if(deta0 > rc && deta1 < rc){
803 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
804 areaj = acc1;
805 }
806 if(deta0 < rc && deta1 > rc){
807 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
808 areaj = acc0;
809 }
810 if(deta0 < rc && deta1 < rc){
811 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
812 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
813 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
814 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
815 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
816 }
5e380e7e 817 fhAreaJet[ijet]->Fill(etac,areaj);
70e58892 818 areabg = areabg - areaj;
819 } // end jets loop
5e380e7e 820 fhAreaBackg->Fill(etac,areabg);
70e58892 821 eta0 = eta1;
822 eta1 = eta1 + etaw;
823 } // end loop for all eta bins
824
825 //subtract background
826 for(Int_t kjet=0; kjet<nJ; kjet++){
827 etJet[kjet] = 0.0; // first clear etJet for this jet
828 for(Int_t bin = 0; bin< ndiv; bin++){
5e380e7e 829 if(fhAreaJet[kjet]->GetBinContent(bin)){
830 Float_t areab = fhAreaBackg->GetBinContent(bin);
831 Float_t etb = fhEtBackg->GetBinContent(bin);
832 Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
833 etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
70e58892 834 }
835 }
836 }
837
838 // calc background total
5e380e7e 839 Double_t etOut = fhEtBackg->Integral();
840 Double_t areaOut = fhAreaBackg->Integral();
7d0f353c 841 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
70e58892 842 etbgTotalN = etOut*areaT/areaOut;
70e58892 843}
844
845////////////////////////////////////////////////////////////////////////
846
847
848void AliUA1JetFinderV1::Reset()
849{
850 fLego->Reset();
ad8ada26 851 AliJetFinder::Reset();
70e58892 852}
853
854////////////////////////////////////////////////////////////////////////
855
446dbc09 856void AliUA1JetFinderV1::WriteJHeaderToFile() const
70e58892 857{
7d0f353c 858 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
7d0f353c 859 header->Write();
70e58892 860}
861
862////////////////////////////////////////////////////////////////////////
863
864void AliUA1JetFinderV1::Init()
865{
866 // initializes some variables
7d0f353c 867 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
70e58892 868 fLego = new
869 TH2F("legoH","eta-phi",
7d0f353c 870 header->GetLegoNbinEta(), header->GetLegoEtaMin(),
871 header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
872 header->GetLegoPhiMin(), header->GetLegoPhiMax());
c1f55a27 873 // Do not store in current dir
874 fLego->SetDirectory(0);
9dda5307 875
70e58892 876}