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