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 **************************************************************************/
16 //---------------------------------------------------------------------
17 // UA1 Cone Algorithm Jet finder
18 // manages the search for jets
19 // Author: Rafael.Diaz.Valdes@cern.ch
21 //---------------------------------------------------------------------
23 #include <TLorentzVector.h>
28 #include "AliUA1JetFinderV1.h"
29 #include "AliUA1JetHeaderV1.h"
30 #include "AliJetReaderHeader.h"
31 #include "AliJetReader.h"
35 ClassImp(AliUA1JetFinderV1)
37 ////////////////////////////////////////////////////////////////////////
39 AliUA1JetFinderV1::AliUA1JetFinderV1()
47 ////////////////////////////////////////////////////////////////////////
49 AliUA1JetFinderV1::~AliUA1JetFinderV1()
55 ////////////////////////////////////////////////////////////////////////
58 void AliUA1JetFinderV1::FindJets()
61 //1) Fill cell map array
62 //2) calculate total energy and fluctuation level
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
71 // transform input to pt,eta,phi plus lego
72 TClonesArray *lvArray = fReader->GetMomentumArray();
73 Int_t nIn = lvArray->GetEntries();
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];
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);
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);
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]);
97 fJets->SetNinput(nIn);
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);
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();
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);
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();
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);
148 if(etbgTotalN != 0.0)
149 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
152 etbgTotal = etbgTotalN; // update with new background estimation
156 Int_t* idxjets = new Int_t[nj];
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;
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];
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]];
180 fJets->SetNCells(ncells);
181 fJets->SetPtFromSignal(percentage);
182 fJets->SetMultiplicities(mult);
183 fJets->SetInJet(injet);
184 fJets->SetEtaIn(etaT);
185 fJets->SetPhiIn(phiT);
187 fJets->SetEtAvg(etbgTotal/(4*(fHeader->GetLegoEtaMax())*TMath::Pi()));
211 ////////////////////////////////////////////////////////////////////////
213 void 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)
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
226 TAxis* xaxis = fLego->GetXaxis();
227 TAxis* yaxis = fLego->GetYaxis();
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);
236 etaCell[nCell] = eta;
237 phiCell[nCell] = phi;
238 flagCell[nCell] = 0; //default
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();
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];
261 Int_t * index = new Int_t[nCell];
262 TMath::Sort(nCell, etCell, index);
263 // variable used in centroide loop
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];
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;
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);
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));
331 }else{ // store this loop information
339 }//end of cells loop looking centroide
341 //avoid cones overloap (to be implemented in the future)
343 //flag cells in Rc, estimate total energy in cone
344 Float_t etCone = 0.0;
346 rc = fHeader->GetRadius();
347 for(Int_t ncell =0; ncell < nCell; ncell++){
348 if(flagCell[ncell] != 0) continue; // cell used before
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];
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);
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){
373 flagCell[mcell] = 1; //flag cell as used
375 flagCell[mcell] = 0; // leave it free
378 //store tmp jet info !!!
379 // reject tmp outside acceptable eta range
380 if ((eta < (fHeader->GetJetEtaMax())) &&
381 (eta > (fHeader->GetJetEtaMin())) &&
383 etaAlgoJet[nJets] = eta;
384 phiAlgoJet[nJets] = phi;
385 etAlgoJet[nJets] = etCone;
386 ncellsAlgoJet[nJets] = nCellIn;
390 } // end of cells loop
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]];
410 ////////////////////////////////////////////////////////////////////////
412 void 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)
417 //background subtraction using cone method but without correction in dE/deta distribution
419 //calculate energy inside and outside cones
420 Float_t rc= fHeader->GetRadius();
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];
439 if(injet[jpart] == 0) etOut += ptT[jpart]; // particle outside cones
440 } //end particle loop
442 //estimate jets and background areas
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);
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);
457 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
458 areaOut = areaOut - areaJet[k];
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
466 // estimate new total background
467 Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
468 etbgTotalN = etOut*areaT/areaOut;
473 ////////////////////////////////////////////////////////////////////////
475 void 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)
481 //background subtraction using statistical method
483 Float_t etbgStat = fHeader->GetBackgStat(); // pre-calculated background
485 //calculate energy inside
486 Float_t rc= fHeader->GetRadius();
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];
505 } //end particle loop
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);
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);
522 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
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
531 etbgTotalN = etbgStat;
535 ////////////////////////////////////////////////////////////////////////
537 void 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)
542 // Cone background subtraction method taking into acount dEt/deta distribution
545 Float_t rc= fHeader->GetRadius();
546 Float_t etamax = fHeader->GetLegoEtaMax();
547 Float_t etamin = fHeader->GetLegoEtaMin();
550 // jet energy and area arrays
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);
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);
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
575 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
580 if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
581 } //end particle loop
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;
595 if(deta0 > rc && deta1 < rc){
596 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
599 if(deta0 < rc && deta1 > rc){
600 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
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
610 hAreaJet[ijet]->Fill(etac,areaj);
611 areabg = areabg - areaj;
613 hAreaBackg->Fill(etac,areabg);
616 } // end loop for all eta bins
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
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;
638 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
640 delete hAreaJet[ljet];
647 ////////////////////////////////////////////////////////////////////////
650 void 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)
655 // Ratio background subtraction method taking into acount dEt/deta distribution
657 //factor F calc before
658 Float_t bgRatioCut = fHeader->GetBackgCutRatio();
662 Float_t rc= fHeader->GetRadius();
663 Float_t etamax = fHeader->GetLegoEtaMax();
664 Float_t etamin = fHeader->GetLegoEtaMin();
667 // jet energy and area arrays
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
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
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];
697 if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
698 } //end particle loop
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;
712 if(deta0 > rc && deta1 < rc){
713 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
716 if(deta0 < rc && deta1 > rc){
717 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
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
727 hAreaJet[ijet]->Fill(etac,areaj);
728 areabg = areabg - areaj;
730 hAreaBackg->Fill(etac,areabg);
733 } // end loop for all eta bins
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
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;
755 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
757 delete hAreaJet[ljet];
764 ////////////////////////////////////////////////////////////////////////
767 void AliUA1JetFinderV1::Reset()
773 ////////////////////////////////////////////////////////////////////////
775 void AliUA1JetFinderV1::WriteJHeaderToFile()
781 ////////////////////////////////////////////////////////////////////////
783 void AliUA1JetFinderV1::Init()
785 // initializes some variables
789 TH2F("legoH","eta-phi",
790 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
791 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
792 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());