]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliUA1JetFinderV2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding violations corrected.
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinderV2.cxx
CommitLineData
ee7de0dd 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//---------------------------------------------------------------------
8838ab7a 19// UA1 Cone Algorithm Jet finder for charged + neutral jet studies
20// manages the search for jets using charged particle momentum and
21// neutral cell energy information
22// Based on UA1 V1 (from R. Diaz)
23// Author: magali.estienne@subatech.in2p3.fr
ee7de0dd 24//---------------------------------------------------------------------
25
ee7de0dd 26#include <TClonesArray.h>
ee7de0dd 27#include <TH1F.h>
28#include <TH2F.h>
29#include <TLorentzVector.h>
30#include <TMath.h>
31#include <TRefArray.h>
be6e5811 32#include "TFile.h"
ee7de0dd 33
34#include "AliUA1JetFinderV2.h"
35#include "AliUA1JetHeaderV1.h"
36#include "AliJetUnitArray.h"
ee7de0dd 37
be6e5811 38class TArrayF;
39class TFile;
40class AliJetReader;
41class AliAODJet;
ee7de0dd 42
43ClassImp(AliUA1JetFinderV2)
44
ee7de0dd 45
8838ab7a 46////////////////////////////////////////////////////////////////////////
9e4cc50d 47AliUA1JetFinderV2::AliUA1JetFinderV2() :
8838ab7a 48 AliJetFinder(),
49 fLego(0),
50 fDebug(0),
51 fOpt(0)
ee7de0dd 52{
8838ab7a 53 //
ee7de0dd 54 // Constructor
8838ab7a 55 //
ee7de0dd 56}
57
58////////////////////////////////////////////////////////////////////////
ee7de0dd 59AliUA1JetFinderV2::~AliUA1JetFinderV2()
ee7de0dd 60{
8838ab7a 61 //
62 // Destructor
63 //
ee7de0dd 64}
65
66////////////////////////////////////////////////////////////////////////
8838ab7a 67void AliUA1JetFinderV2::FindJetsC()
68{
69 //
70 // Used to find jets using charged particle momentum information
71 //
72 // 1) Fill cell map array
73 // 2) calculate total energy and fluctuation level
74 // 3) Run algorithm
75 // 3.1) look centroides in cell map
76 // 3.2) calculate total energy in cones
77 // 3.3) flag as a possible jet
78 // 3.4) reorder cones by energy
79 // 4) subtract backg in accepted jets
80 // 5) fill AliJet list
81
82 // Transform input to pt,eta,phi plus lego
83
84 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
85 TClonesArray* lvArray = fReader->GetMomentumArray();
86 Int_t nIn = lvArray->GetEntries();
87
88 if (nIn == 0) return;
89
90 // local arrays for input
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 Float_t* cFlagT = new Float_t[nIn]; // Temporarily added
95 Float_t* sFlagT = new Float_t[nIn]; // Temporarily added
96 Int_t* injet = new Int_t[nIn];
97
98 //total energy in array
99 Float_t etbgTotal = 0.0;
100 TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
101
102 // load input vectors and calculate total energy in array
103 for (Int_t i = 0; i < nIn; i++){
104 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
105 ptT[i] = lv->Pt();
106 etaT[i] = lv->Eta();
107 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
be6e5811 108 cFlagT[i] = fReader->GetCutFlag(i);
109 sFlagT[i] = fReader->GetSignalFlag(i);
8838ab7a 110
111 if (fReader->GetCutFlag(i) != 1) continue;
112 fLego->Fill(etaT[i], phiT[i], ptT[i]);
113 hPtTotal->Fill(ptT[i]);
114 etbgTotal+= ptT[i];
115 }
116
117 fJets->SetNinput(nIn);
118
119 // calculate total energy and fluctuation in map
120 Double_t meanpt = hPtTotal->GetMean();
121 Double_t ptRMS = hPtTotal->GetRMS();
122 Double_t npart = hPtTotal->GetEntries();
123 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
124
125 // arrays to hold jets
be6e5811 126 Float_t* etaJet = new Float_t[30]; // eta jet
127 Float_t* phiJet = new Float_t[30]; // phi jet
128 Float_t* etJet = new Float_t[30]; // et jet
129 Float_t* etsigJet = new Float_t[30]; // signal et in jet
8838ab7a 130 Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
131 Int_t* ncellsJet = new Int_t[30];
132 Int_t* multJet = new Int_t[30];
133 //--- Added for jet reordering at the end of the jet finding procedure
134 Float_t* etaJetOk = new Float_t[30];
135 Float_t* phiJetOk = new Float_t[30];
136 Float_t* etJetOk = new Float_t[30];
be6e5811 137 Float_t* etsigJetOk = new Float_t[30]; // signal et in jet
8838ab7a 138 Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
139 Int_t* ncellsJetOk = new Int_t[30];
140 Int_t* multJetOk = new Int_t[30];
141 //--------------------------
142 Int_t nJets; // to hold number of jets found by algorithm
143 Int_t nj; // number of jets accepted
144 Float_t prec = header->GetPrecBg();
145 Float_t bgprec = 1;
146 while(bgprec > prec){
147 //reset jet arrays in memory
148 memset(etaJet,0,sizeof(Float_t)*30);
149 memset(phiJet,0,sizeof(Float_t)*30);
150 memset(etJet,0,sizeof(Float_t)*30);
151 memset(etallJet,0,sizeof(Float_t)*30);
152 memset(etsigJet,0,sizeof(Float_t)*30);
153 memset(ncellsJet,0,sizeof(Int_t)*30);
154 memset(multJet,0,sizeof(Int_t)*30);
155 //--- Added for jet reordering at the end of the jet finding procedure
156 memset(etaJetOk,0,sizeof(Float_t)*30);
157 memset(phiJetOk,0,sizeof(Float_t)*30);
158 memset(etJetOk,0,sizeof(Float_t)*30);
159 memset(etallJetOk,0,sizeof(Float_t)*30);
160 memset(etsigJetOk,0,sizeof(Float_t)*30);
161 memset(ncellsJetOk,0,sizeof(Int_t)*30);
162 memset(multJetOk,0,sizeof(Int_t)*30);
163 //--------------------------
164 nJets = 0;
165 nj = 0;
166
167 // reset particles-jet array in memory
168 memset(injet,-1,sizeof(Int_t)*nIn);
169 //run cone algorithm finder
170 RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
171
172 //run background subtraction
173 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
174 nj = header->GetNAcceptJets();
175 else
176 nj = nJets;
177 //subtract background
178 Float_t etbgTotalN = 0.0; //new background
179 if(header->GetBackgMode() == 1) // standard
8838ab7a 180 SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
181 if(header->GetBackgMode() == 2) //cone
182 SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
183 if(header->GetBackgMode() == 3) //ratio
184 SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
185 if(header->GetBackgMode() == 4) //statistic
186 SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
187 //calc precision
188 if(etbgTotalN != 0.0)
189 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
190 else
191 bgprec = 0;
192 etbgTotal = etbgTotalN; // update with new background estimation
193 } //end while
194
195 // add jets to list
196 Int_t* idxjets = new Int_t[nj];
197 Int_t nselectj = 0;
198 printf("Found %d jets \n", nj);
ee7de0dd 199
8838ab7a 200 // Reorder jets by et in cone
201 Int_t * idx = new Int_t[nJets];
202 TMath::Sort(nJets, etJet, idx);
203 for(Int_t p = 0; p < nJets; p++){
204 etaJetOk[p] = etaJet[idx[p]];
205 phiJetOk[p] = phiJet[idx[p]];
206 etJetOk[p] = etJet[idx[p]];
207 etallJetOk[p] = etJet[idx[p]];
be6e5811 208 etsigJetOk[p] = etsigJet[idx[p]];
8838ab7a 209 ncellsJetOk[p] = ncellsJet[idx[p]];
210 multJetOk[p] = multJet[idx[p]];
211 }
212
213 for(Int_t kj=0; kj<nj; kj++)
214 {
215 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
216 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
217 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
218 Float_t px, py,pz,en; // convert to 4-vector
219 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
220 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
221 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
222 en = TMath::Sqrt(px * px + py * py + pz * pz);
223 fJets->AddJet(px, py, pz, en);
224
225 AliAODJet jet(px, py, pz, en);
226 jet.Print("");
227
228 AddJet(jet);
229
230 idxjets[nselectj] = kj;
231 nselectj++;
232 }
ee7de0dd 233
8838ab7a 234 //add signal percentage and total signal in AliJets for analysis tool
235 Float_t* percentage = new Float_t[nselectj];
236 Int_t* ncells = new Int_t[nselectj];
237 Int_t* mult = new Int_t[nselectj];
238 for(Int_t i = 0; i< nselectj; i++)
239 {
240 percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
241 ncells[i] = ncellsJetOk[idxjets[i]];
242 mult[i] = multJetOk[idxjets[i]];
243 }
244
245 //add particle-injet relationship ///
246 for(Int_t bj = 0; bj < nIn; bj++)
247 {
248 if(injet[bj] == -1) continue; //background particle
249 Int_t bflag = 0;
250 for(Int_t ci = 0; ci< nselectj; ci++){
251 if(injet[bj] == idxjets[ci]){
252 injet[bj]= ci;
253 bflag++;
254 break;
255 }
256 }
257 if(bflag == 0) injet[bj] = -1; // set as background particle
258 }
259
260 fJets->SetNCells(ncells);
261 fJets->SetPtFromSignal(percentage);
262 fJets->SetMultiplicities(mult);
263 fJets->SetInJet(injet);
264 fJets->SetEtaIn(etaT);
265 fJets->SetPhiIn(phiT);
266 fJets->SetPtIn(ptT);
267 fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi()));
268
269 //delete
270 delete[] ptT;
271 delete[] etaT;
272 delete[] phiT;
273 delete[] cFlagT;
274 delete[] sFlagT;
275 delete[] injet;
276 delete[] hPtTotal;
277 delete[] etaJet;
278 delete[] phiJet;
279 delete[] etJet;
280 delete[] etsigJet;
281 delete[] etallJet;
282 delete[] ncellsJet;
283 delete[] multJet;
284 delete[] idxjets;
285 delete[] percentage;
286 delete[] ncells;
287 delete[] mult;
288 //--- Added for jet reordering
289 delete etaJetOk;
290 delete phiJetOk;
291 delete etJetOk;
292 delete etsigJetOk;
293 delete etallJetOk;
294 delete ncellsJetOk;
295 delete multJetOk;
296 //--------------------------
297
298}
ee7de0dd 299
8838ab7a 300////////////////////////////////////////////////////////////////////////
301void AliUA1JetFinderV2::FindJets()
ee7de0dd 302{
8838ab7a 303 //
304 // Used to find jets using charged particle momentum information
305 // & neutral energy from calo cells
306 //
307 // 1) Fill cell map array
308 // 2) calculate total energy and fluctuation level
309 // 3) Run algorithm
310 // 3.1) look centroides in cell map
311 // 3.2) calculate total energy in cones
312 // 3.3) flag as a possible jet
313 // 3.4) reorder cones by energy
314 // 4) subtract backg in accepted jets
315 // 5) fill AliJet list
ee7de0dd 316
317 // transform input to pt,eta,phi plus lego
318
8838ab7a 319 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
320 TClonesArray* fUnit = fReader->GetUnitArray();
321 Int_t nCand = fReader->GetNumCandidate();
322 Int_t nCandCut = fReader->GetNumCandidateCut();
323 Int_t nIn = fUnit->GetEntries();
be6e5811 324 Float_t ptMin = fReader->GetReaderHeader()->GetPtCut();
ee7de0dd 325
be6e5811 326 fDebug = fReader->GetReaderHeader()->GetDebug();
ee7de0dd 327 if (nIn == 0) return;
328
8838ab7a 329 Int_t nCandidateCut = 0;
330 Int_t nCandidate = 0;
331
332 nCandidate = nCand;
333 nCandidateCut = nCandCut;
334
ee7de0dd 335 // local arrays for input No Cuts
336 // Both pt < ptMin and pt > ptMin
be6e5811 337 Float_t* ptT = new Float_t[nCandidate];
338 Float_t* en2T = new Float_t[nCandidate];
339 Float_t* pt2T = new Float_t[nCandidate];
340 Int_t* detT = new Int_t[nCandidate];
341 Float_t* etaT = new Float_t[nCandidate];
342 Float_t* phiT = new Float_t[nCandidate];
343 Float_t* cFlagT = new Float_t[nCandidate];
344 Float_t* cFlag2T = new Float_t[nCandidate];
345 Float_t* sFlagT = new Float_t[nCandidate];
346 Float_t* cClusterT = new Float_t[nCandidate];
347 Int_t* vectT = new Int_t[nCandidate];
348 Int_t loop1 = 0;
349 Int_t* injet = new Int_t[nCandidate];
350 Int_t* sflag = new Int_t[nCandidate];
351 TRefArray* trackRef = new TRefArray();
ee7de0dd 352
353 //total energy in array
354 Float_t etbgTotal = 0.0;
be6e5811 355 TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
ee7de0dd 356
357 // Input cell info
be6e5811 358 Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
359 Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
360 Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
361 Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
362 Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
363 Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
364 Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
8838ab7a 365 Int_t *flagCell2 = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
366
ee7de0dd 367 // Information extracted from fUnitArray
8838ab7a 368 // Load input vectors and calculate total energy in array
ee7de0dd 369 for(Int_t i=0; i<nIn; i++)
370 {
8838ab7a 371 // Recover particle information from UnitArray
372
ee7de0dd 373 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
be6e5811 374 TRefArray* ref = uArray->GetUnitTrackRef();
375 Int_t nRef = ref->GetEntries();
376
ee7de0dd 377 if(uArray->GetUnitEnergy()>0.){
be6e5811 378
379 for(Int_t jpart=0; jpart<nRef;jpart++)
380 trackRef->Add((AliVTrack*)ref->At(jpart));
ee7de0dd 381 ptT[loop1] = uArray->GetUnitEnergy();
8838ab7a 382 detT[loop1] = uArray->GetUnitDetectorFlag();
ee7de0dd 383 etaT[loop1] = uArray->GetUnitEta();
384 phiT[loop1] = uArray->GetUnitPhi();
8838ab7a 385 cFlagT[loop1]= uArray->GetUnitCutFlag(); // pt cut tpc
386 cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
387 sFlagT[loop1]= uArray->GetUnitSignalFlag();
be6e5811 388 vectT[loop1] = nRef;
8838ab7a 389 if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
390 pt2T[loop1] = 0.;
391 en2T[loop1] = 0.;
392 if(detT[loop1]==1){
393 en2T[loop1] = ptT[loop1] - header->GetMinCellEt();
394 if(en2T[loop1] < 0) en2T[loop1]=0;
395 hPtTotal->Fill(en2T[loop1]);
396 etbgTotal += en2T[loop1];
397 }
398 if(detT[loop1]==0){ // TPC+ITS
399 Float_t pt = 0.;
be6e5811 400 for(Int_t j=0; j<nRef;j++){
8838ab7a 401 Float_t x=0.; Float_t y=0.; Float_t z=0.;
be6e5811 402 x = ((AliVTrack*)ref->At(j))->Px();
403 y = ((AliVTrack*)ref->At(j))->Py();
404 z = ((AliVTrack*)ref->At(j))->Pz();
8838ab7a 405 pt = TMath::Sqrt(x*x+y*y);
be6e5811 406 if(pt>ptMin) {
8838ab7a 407 pt2T[loop1] += pt;
408 en2T[loop1] += pt;
409 hPtTotal->Fill(pt);
410 etbgTotal+= pt;
411 }
412 }
413 }
414 if(detT[loop1]==2) { // EMCal
415 Float_t ptCTot = 0.;
416 Float_t pt = 0.;
417 Float_t enC = 0.;
be6e5811 418 for(Int_t j=0; j<nRef;j++){
8838ab7a 419 Float_t x=0.; Float_t y=0.; Float_t z=0.;
be6e5811 420 x = ((AliVTrack*)ref->At(j))->Px();
421 y = ((AliVTrack*)ref->At(j))->Py();
422 z = ((AliVTrack*)ref->At(j))->Pz();
8838ab7a 423 pt = TMath::Sqrt(x*x+y*y);
be6e5811 424 if(pt>ptMin) {
8838ab7a 425 pt2T[loop1]+=pt;
426 en2T[loop1]+=pt;
427 hPtTotal->Fill(pt);
428 etbgTotal+= pt;
429 }
430 ptCTot += pt;
431 }
432 enC = ptT[loop1] - ptCTot - header->GetMinCellEt();
433 if(enC < 0.) enC=0.;
434 en2T[loop1] += enC;
435 hPtTotal->Fill(enC);
436 etbgTotal+= enC;
437 }
ee7de0dd 438 }
439 loop1++;
440 }
ee7de0dd 441
8838ab7a 442 if(uArray->GetUnitCutFlag()==1) {
443 if(uArray->GetUnitDetectorFlag()==1){ // EMCal case
444 etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt();
445 if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.;
446 etaCell[i] = uArray->GetUnitEta();
447 phiCell[i] = uArray->GetUnitPhi();
448 flagCell[i] = 0; // default
449 etCell2[i] = etCell[i];
450 etaCell2[i] = uArray->GetUnitEta();
451 phiCell2[i] = uArray->GetUnitPhi();
452 flagCell2[i] = 0; // default
453 }
454 if(uArray->GetUnitDetectorFlag()==0){ // TPC case
455 Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
be6e5811 456 for(Int_t j=0; j<nRef;j++)
8838ab7a 457 {
458 Float_t x=0.; Float_t y=0.; Float_t z=0.;
be6e5811 459 x = ((AliVTrack*)ref->At(j))->Px();
460 y = ((AliVTrack*)ref->At(j))->Py();
461 z = ((AliVTrack*)ref->At(j))->Pz();
8838ab7a 462 pt = TMath::Sqrt(x*x+y*y);
be6e5811 463 if(pt>ptMin) {
8838ab7a 464 et1 += pt;
465 et2 += pt;
466 }
467 }
468 etCell[i] = et1;
469 etCell2[i] = et2;
470 if(et1 < 0.) etCell[i] = etCell2[i] = 0.;
471 etaCell[i] = uArray->GetUnitEta();
472 phiCell[i] = uArray->GetUnitPhi();
473 flagCell[i] = 0; // default
474 etaCell2[i] = uArray->GetUnitEta();
475 phiCell2[i] = uArray->GetUnitPhi();
476 flagCell2[i] = 0; // default
477 }
478 if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case
479 Float_t ptCTot = 0.;
480 Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
481 Float_t enC = 0.;
be6e5811 482 for(Int_t j=0; j<nRef;j++)
8838ab7a 483 {
484 Float_t x=0.; Float_t y=0.; Float_t z=0.;
be6e5811 485 x = ((AliVTrack*)ref->At(j))->Px();
486 y = ((AliVTrack*)ref->At(j))->Py();
487 z = ((AliVTrack*)ref->At(j))->Pz();
8838ab7a 488 pt = TMath::Sqrt(x*x+y*y);
be6e5811 489 if(pt>ptMin) {
8838ab7a 490 et1 += pt;
491 et2 += pt;
492 }
493 ptCTot += pt;
494 }
495 enC = uArray->GetUnitEnergy() - ptCTot;
496 etCell[i] = et1 + enC - header->GetMinCellEt();
497 etCell2[i] = et2 + enC - header->GetMinCellEt();
498 if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.;
499 etaCell[i] = uArray->GetUnitEta();
500 phiCell[i] = uArray->GetUnitPhi();
501 flagCell[i] = 0; // default
502 etaCell2[i] = uArray->GetUnitEta();
503 phiCell2[i] = uArray->GetUnitPhi();
504 flagCell2[i] = 0; // default
505 }
506 }
507 else {
508 etCell[i] = 0.;
509 etaCell[i] = uArray->GetUnitEta();
510 phiCell[i] = uArray->GetUnitPhi();
511 flagCell[i] = 0;
512 etCell2[i] = 0.;
513 etaCell2[i] = uArray->GetUnitEta();
514 phiCell2[i] = uArray->GetUnitPhi();
515 flagCell2[i] = 0;
516 }
517 } // end loop on nCandidate
518
ee7de0dd 519 fJets->SetNinput(nCandidate);
520
521 // calculate total energy and fluctuation in map
522 Double_t meanpt = hPtTotal->GetMean();
523 Double_t ptRMS = hPtTotal->GetRMS();
524 Double_t npart = hPtTotal->GetEntries();
525 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
526
527 // arrays to hold jets
8838ab7a 528 Float_t* etaJet = new Float_t[30];
529 Float_t* phiJet = new Float_t[30];
530 Float_t* etJet = new Float_t[30];
ee7de0dd 531 Float_t* etsigJet = new Float_t[30]; //signal et in jet
8838ab7a 532 Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
533 Int_t* ncellsJet = new Int_t[30];
534 Int_t* multJet = new Int_t[30];
535 //--- Added by me for jet reordering at the end of the jet finding procedure
536 Float_t* etaJetOk = new Float_t[30];
537 Float_t* phiJetOk = new Float_t[30];
538 Float_t* etJetOk = new Float_t[30];
539 Float_t* etsigJetOk = new Float_t[30]; //signal et in jet
540 Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
541 Int_t* ncellsJetOk = new Int_t[30];
542 Int_t* multJetOk = new Int_t[30];
543 //--------------------------
544 Int_t nJets; // to hold number of jets found by algorithm
545 Int_t nj; // number of jets accepted
546 Float_t prec = header->GetPrecBg();
547 Float_t bgprec = 1;
548
ee7de0dd 549 while(bgprec > prec){
ee7de0dd 550
8838ab7a 551 //reset jet arrays in memory
552 memset(etaJet,0,sizeof(Float_t)*30);
553 memset(phiJet,0,sizeof(Float_t)*30);
554 memset(etJet,0,sizeof(Float_t)*30);
555 memset(etallJet,0,sizeof(Float_t)*30);
556 memset(etsigJet,0,sizeof(Float_t)*30);
557 memset(ncellsJet,0,sizeof(Int_t)*30);
558 memset(multJet,0,sizeof(Int_t)*30);
559 //--- Added by me for jet reordering at the end of the jet finding procedure
560 memset(etaJetOk,0,sizeof(Float_t)*30);
561 memset(phiJetOk,0,sizeof(Float_t)*30);
562 memset(etJetOk,0,sizeof(Float_t)*30);
563 memset(etallJetOk,0,sizeof(Float_t)*30);
564 memset(etsigJetOk,0,sizeof(Float_t)*30);
565 memset(ncellsJetOk,0,sizeof(Int_t)*30);
566 memset(multJetOk,0,sizeof(Int_t)*30);
567
568 nJets = 0;
569 nj = 0;
570
571 // reset particles-jet array in memory
572 memset(injet,-1,sizeof(Int_t)*nCandidate);
573 //run cone algorithm finder
574 RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2,
575 flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,
576 etallJet,ncellsJet);
577
578 //run background subtraction
579 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
580 nj = header->GetNAcceptJets();
581 else
582 nj = nJets;
583
584 //subtract background
585 Float_t etbgTotalN = 0.0; //new background
586 if(header->GetBackgMode() == 1) // standard
587 SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
588 // To be modified ------------------------
589 if(header->GetBackgMode() == 2) //cone
590 SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
591 if(header->GetBackgMode() == 3) //ratio
592 SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
593 if(header->GetBackgMode() == 4) //statistic
594 SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
595 //----------------------------------------
596 //calc precision
597 if(etbgTotalN != 0.0)
598 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
599 else
600 bgprec = 0;
601 etbgTotal = etbgTotalN; // update with new background estimation
602 } //end while
603
ee7de0dd 604 // add jets to list
605 Int_t* idxjets = new Int_t[nj];
606 Int_t nselectj = 0;
607 printf("Found %d jets \n", nj);
608
8838ab7a 609 // Reorder jets by et in cone
610 // Sort jets by energy
611 Int_t * idx = new Int_t[nJets];
612 TMath::Sort(nJets, etJet, idx);
613 for(Int_t p = 0; p < nJets; p++)
614 {
615 etaJetOk[p] = etaJet[idx[p]];
616 phiJetOk[p] = phiJet[idx[p]];
617 etJetOk[p] = etJet[idx[p]];
618 etallJetOk[p] = etJet[idx[p]];
be6e5811 619 etsigJetOk[p] = etsigJet[idx[p]];
8838ab7a 620 ncellsJetOk[p] = ncellsJet[idx[p]];
621 multJetOk[p] = multJet[idx[p]];
622 }
623
624 for(Int_t kj=0; kj<nj; kj++)
625 {
626 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
627 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
628 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
ee7de0dd 629 Float_t px, py,pz,en; // convert to 4-vector
8838ab7a 630 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
631 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
632 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
ee7de0dd 633 en = TMath::Sqrt(px * px + py * py + pz * pz);
634 fJets->AddJet(px, py, pz, en);
635 AliAODJet jet(px, py, pz, en);
636 jet.Print("");
637
638 AddJet(jet);
639
640 idxjets[nselectj] = kj;
641 nselectj++;
8838ab7a 642 }
643
ee7de0dd 644 //add signal percentage and total signal in AliJets for analysis tool
645 Float_t* percentage = new Float_t[nselectj];
646 Int_t* ncells = new Int_t[nselectj];
647 Int_t* mult = new Int_t[nselectj];
8838ab7a 648 for(Int_t i = 0; i< nselectj; i++)
649 {
650 percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
651 ncells[i] = ncellsJetOk[idxjets[i]];
652 mult[i] = multJetOk[idxjets[i]];
653 }
654
655 //add particle-injet relationship ///
656 for(Int_t bj = 0; bj < nCandidate; bj++)
657 {
658 if(injet[bj] == -1) continue; //background particle
659 Int_t bflag = 0;
660 for(Int_t ci = 0; ci< nselectj; ci++){
661 if(injet[bj] == idxjets[ci]){
662 injet[bj]= ci;
663 bflag++;
664 break;
665 }
666 }
667 if(bflag == 0) injet[bj] = -1; // set as background particle
ee7de0dd 668 }
8838ab7a 669
ee7de0dd 670 fJets->SetNCells(ncells);
671 fJets->SetPtFromSignal(percentage);
672 fJets->SetMultiplicities(mult);
673 fJets->SetInJet(injet);
674 fJets->SetEtaIn(etaT);
675 fJets->SetPhiIn(phiT);
676 fJets->SetPtIn(ptT);
be6e5811 677 fJets->SetTrackReferences(trackRef);
8838ab7a 678 fJets->SetDetectorFlagIn(detT);
679 fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin())));
ee7de0dd 680
681 //delete
ee7de0dd 682 delete ptT;
8838ab7a 683 delete en2T;
684 delete pt2T;
ee7de0dd 685 delete etaT;
686 delete phiT;
be6e5811 687 trackRef->Delete();
688 delete trackRef;
8838ab7a 689 delete detT;
ee7de0dd 690 delete cFlagT;
8838ab7a 691 delete cFlag2T;
692 delete sFlagT;
ee7de0dd 693 delete cClusterT;
8838ab7a 694 delete vectT;
ee7de0dd 695 delete injet;
696 delete sflag;
697 delete hPtTotal;
698 delete etCell;
699 delete etaCell;
700 delete phiCell;
701 delete flagCell;
8838ab7a 702 delete etCell2;
703 delete etaCell2;
704 delete phiCell2;
705 delete flagCell2;
ee7de0dd 706 delete etaJet;
707 delete phiJet;
708 delete etJet;
709 delete etsigJet;
710 delete etallJet;
711 delete ncellsJet;
712 delete multJet;
8838ab7a 713 //--- Added for jet reordering
714 delete etaJetOk;
715 delete phiJetOk;
716 delete etJetOk;
717 delete etsigJetOk;
718 delete etallJetOk;
719 delete ncellsJetOk;
720 delete multJetOk;
721 //--------------------------
ee7de0dd 722 delete idxjets;
723 delete percentage;
724 delete ncells;
725 delete mult;
726
ee7de0dd 727}
728
729////////////////////////////////////////////////////////////////////////
ee7de0dd 730void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
be6e5811 731 Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2,
732 const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal,
733 Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet,
734 Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
ee7de0dd 735{
be6e5811 736 //
737 // Main method for jet finding
738 // UA1 base cone finder
739 //
740
8838ab7a 741 Int_t nCell = nIn;
be6e5811 742 fDebug = fReader->GetReaderHeader()->GetDebug();
ee7de0dd 743
8838ab7a 744 // Dump lego
745 // Check enough space! *to be done*
ee7de0dd 746 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
8838ab7a 747 for(Int_t i=0; i<nCell; i++){
748 etCell[i] = etCell2[i];
749 etaCell[i] = etaCell2[i];
750 phiCell[i] = phiCell2[i];
751 flagCell[i] = flagCell2[i];
752 }
ee7de0dd 753
754 // Parameters from header
755 Float_t minmove = header->GetMinMove();
756 Float_t maxmove = header->GetMaxMove();
757 Float_t rc = header->GetRadius();
758 Float_t etseed = header->GetEtSeed();
759
8838ab7a 760 // Tmp array of jets form algoritm
ee7de0dd 761 Float_t etaAlgoJet[30];
762 Float_t phiAlgoJet[30];
763 Float_t etAlgoJet[30];
764 Int_t ncellsAlgoJet[30];
765
8838ab7a 766 // Run algorithm//
ee7de0dd 767
8838ab7a 768 // Sort cells by et
ee7de0dd 769 Int_t * index = new Int_t[nCell];
770 TMath::Sort(nCell, etCell, index);
771
8838ab7a 772 // Variable used in centroide loop
9dda5307 773 Float_t eta = 0.0;
774 Float_t phi = 0.0;
775 Float_t eta0 = 0.0;
776 Float_t phi0 = 0.0;
777 Float_t etab = 0.0;
778 Float_t phib = 0.0;
779 Float_t etas = 0.0;
780 Float_t phis = 0.0;
781 Float_t ets = 0.0;
782 Float_t deta = 0.0;
783 Float_t dphi = 0.0;
784 Float_t dr = 0.0;
785 Float_t etsb = 0.0;
ee7de0dd 786 Float_t etasb = 0.0;
787 Float_t phisb = 0.0;
9dda5307 788 Float_t dphib = 0.0;
ee7de0dd 789
8838ab7a 790 for(Int_t icell = 0; icell < nCell; icell++)
791 {
9dda5307 792 Int_t jcell = index[icell];
793 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
794 if(flagCell[jcell] != 0) continue; // if cell was used before
8838ab7a 795
9dda5307 796 eta = etaCell[jcell];
797 phi = phiCell[jcell];
798 eta0 = eta;
799 phi0 = phi;
800 etab = eta;
801 phib = phi;
802 ets = etCell[jcell];
803 etas = 0.0;
804 phis = 0.0;
805 etsb = ets;
806 etasb = 0.0;
807 phisb = 0.0;
8838ab7a 808 for(Int_t kcell =0; kcell < nCell; kcell++)
809 {
9dda5307 810 Int_t lcell = index[kcell];
811 if(lcell == jcell) continue; // cell itself
812 if(flagCell[lcell] != 0) continue; // cell used before
8838ab7a 813 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
9dda5307 814 //calculate dr
815 deta = etaCell[lcell] - eta;
8838ab7a 816 dphi = TMath::Abs(phiCell[lcell] - phi);
9dda5307 817 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
818 dr = TMath::Sqrt(deta * deta + dphi * dphi);
819 if(dr <= rc){
8838ab7a 820 // calculate offset from initiate cell
821 deta = etaCell[lcell] - eta0;
822 dphi = phiCell[lcell] - phi0;
823 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
824 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
825 etas = etas + etCell[lcell]*deta;
826 phis = phis + etCell[lcell]*dphi;
827 ets = ets + etCell[lcell];
828 //new weighted eta and phi including this cell
829 eta = eta0 + etas/ets;
830 phi = phi0 + phis/ets;
831 // if cone does not move much, just go to next step
832 dphib = TMath::Abs(phi - phib);
833 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
834 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
835 if(dr <= minmove) break;
836 // cone should not move more than max_mov
837 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
838 if(dr > maxmove){
839 eta = etab;
840 phi = phib;
841 ets = etsb;
842 etas = etasb;
843 phis = phisb;
844 } else { // store this loop information
845 etab = eta;
846 phib = phi;
847 etsb = ets;
848 etasb = etas;
849 phisb = phis;
850 }
851 } // inside cone
ee7de0dd 852 }//end of cells loop looking centroide
853
854 //avoid cones overloap (to be implemented in the future)
855
856 //flag cells in Rc, estimate total energy in cone
8838ab7a 857 Float_t etCone = 0.0;
858 Int_t nCellIn = 0;
859 Int_t nCellOut = 0;
860 rc = header->GetRadius();
861
862 for(Int_t ncell =0; ncell < nCell; ncell++)
863 {
864 if(flagCell[ncell] != 0) continue; // cell used before
865 //calculate dr
866 deta = etaCell[ncell] - eta;
867 // if(deta <= rc){ // Added to improve velocity -> to be tested
868 dphi = phiCell[ncell] - phi;
869 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
870 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
871 // if(dphi <= rc){ // Added to improve velocity -> to be tested
872 dr = TMath::Sqrt(deta * deta + dphi * dphi);
873 if(dr <= rc){ // cell in cone
874 flagCell[ncell] = -1;
875 etCone+=etCell[ncell];
876 nCellIn++;
877 }
878 else nCellOut++;
879 // } // end deta <= rc
880 // } // end dphi <= rc
ee7de0dd 881 }
882
8838ab7a 883 // select jets with et > background
884 // estimate max fluctuation of background in cone
885 Double_t ncellin = (Double_t)nCellIn;
886 Double_t ntcell = (Double_t)nCell;
887 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/(ntcell));
888 // min cone et
889 Double_t etcmin = etCone ; // could be used etCone - etmin !!
890 //decisions !! etbmax < etcmin
891
892 for(Int_t mcell =0; mcell < nCell; mcell++)
893 {
894 if(flagCell[mcell] == -1){
895 if(etbmax < etcmin)
896 flagCell[mcell] = 1; //flag cell as used
897 else
898 flagCell[mcell] = 0; // leave it free
899 }
ee7de0dd 900 }
8838ab7a 901 //store tmp jet info !!!
902 if(etbmax < etcmin)
903 {
904 etaAlgoJet[nJets] = eta;
905 phiAlgoJet[nJets] = phi;
906 etAlgoJet[nJets] = etCone;
907 ncellsAlgoJet[nJets] = nCellIn;
908 nJets++;
909 }
910
911 } // end of cells loop
912
913 for(Int_t p = 0; p < nJets; p++)
914 {
915 etaJet[p] = etaAlgoJet[p];
916 phiJet[p] = phiAlgoJet[p];
917 etJet[p] = etAlgoJet[p];
918 etallJet[p] = etAlgoJet[p];
919 ncellsJet[p] = ncellsAlgoJet[p];
920 }
921
922 //delete
923 delete index;
924
925}
926
927////////////////////////////////////////////////////////////////////////
928void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
929 Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
930 Float_t* etallJet, Int_t* ncellsJet)
931{
932 // Dump lego
933 // Check enough space! *to be done*
934 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
935 Float_t etCell[60000]; //! Cell Energy
936 Float_t etaCell[60000]; //! Cell eta
937 Float_t phiCell[60000]; //! Cell phi
938 Int_t flagCell[60000]; //! Cell flag
939
940 Int_t nCell = 0;
941 TAxis* xaxis = fLego->GetXaxis();
942 TAxis* yaxis = fLego->GetYaxis();
943 Float_t e = 0.0;
944 for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++)
945 {
946 for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++)
947 {
948 e = fLego->GetBinContent(i,j);
949 if (e < 0.0) continue; // don't include this cells
950 Float_t eta = xaxis->GetBinCenter(i);
951 Float_t phi = yaxis->GetBinCenter(j);
952 etCell[nCell] = e;
953 etaCell[nCell] = eta;
954 phiCell[nCell] = phi;
955 flagCell[nCell] = 0; //default
956 nCell++;
957 }
958 }
959
960 // Parameters from header
961 Float_t minmove = header->GetMinMove();
962 Float_t maxmove = header->GetMaxMove();
963 Float_t rc = header->GetRadius();
964 Float_t etseed = header->GetEtSeed();
965
966 // Tmp array of jets form algoritm
967 Float_t etaAlgoJet[30];
968 Float_t phiAlgoJet[30];
969 Float_t etAlgoJet[30];
970 Int_t ncellsAlgoJet[30];
971
972 // Run algorithm//
973
974 // Sort cells by et
975 Int_t * index = new Int_t[nCell];
976 TMath::Sort(nCell, etCell, index);
977 // variable used in centroide loop
978 Float_t eta = 0.0;
979 Float_t phi = 0.0;
980 Float_t eta0 = 0.0;
981 Float_t phi0 = 0.0;
982 Float_t etab = 0.0;
983 Float_t phib = 0.0;
984 Float_t etas = 0.0;
985 Float_t phis = 0.0;
986 Float_t ets = 0.0;
987 Float_t deta = 0.0;
988 Float_t dphi = 0.0;
989 Float_t dr = 0.0;
990 Float_t etsb = 0.0;
991 Float_t etasb = 0.0;
992 Float_t phisb = 0.0;
993 Float_t dphib = 0.0;
ee7de0dd 994
8838ab7a 995 for(Int_t icell = 0; icell < nCell; icell++)
996 {
997 Int_t jcell = index[icell];
998 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
999 if(flagCell[jcell] != 0) continue; // if cell was used before
1000
1001 eta = etaCell[jcell];
1002 phi = phiCell[jcell];
1003 eta0 = eta;
1004 phi0 = phi;
1005 etab = eta;
1006 phib = phi;
1007 ets = etCell[jcell];
1008 etas = 0.0;
1009 phis = 0.0;
1010 etsb = ets;
1011 etasb = 0.0;
1012 phisb = 0.0;
1013 for(Int_t kcell =0; kcell < nCell; kcell++)
1014 {
1015 Int_t lcell = index[kcell];
1016 if(lcell == jcell) continue; // cell itself
1017 if(flagCell[lcell] != 0) continue; // cell used before
1018 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
1019 //calculate dr
1020 deta = etaCell[lcell] - eta;
1021 dphi = TMath::Abs(phiCell[lcell] - phi);
1022 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1023 dr = TMath::Sqrt(deta * deta + dphi * dphi);
1024 if(dr <= rc)
1025 {
1026 // calculate offset from initiate cell
1027 deta = etaCell[lcell] - eta0;
1028 dphi = phiCell[lcell] - phi0;
1029 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
1030 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
1031 etas = etas + etCell[lcell]*deta;
1032 phis = phis + etCell[lcell]*dphi;
1033 ets = ets + etCell[lcell];
1034 //new weighted eta and phi including this cell
1035 eta = eta0 + etas/ets;
1036 phi = phi0 + phis/ets;
1037 // if cone does not move much, just go to next step
1038 dphib = TMath::Abs(phi - phib);
1039 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
1040 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
1041 if(dr <= minmove) break;
1042 // cone should not move more than max_mov
1043 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
1044 if(dr > maxmove){
1045 eta = etab;
1046 phi = phib;
1047 ets = etsb;
1048 etas = etasb;
1049 phis = phisb;
1050 } else { // store this loop information
1051 etab=eta;
1052 phib=phi;
1053 etsb = ets;
1054 etasb = etas;
1055 phisb = phis;
1056 }
1057 } // inside cone
1058 }//end of cells loop looking centroide
1059
1060 // Avoid cones overloap (to be implemented in the future)
1061
1062 // Flag cells in Rc, estimate total energy in cone
1063 Float_t etCone = 0.0;
1064 Int_t nCellIn = 0;
1065 Int_t nCellOut = 0;
1066 rc = header->GetRadius();
1067 for(Int_t ncell =0; ncell < nCell; ncell++)
1068 {
1069 if(flagCell[ncell] != 0) continue; // cell used before
1070 //calculate dr
1071 deta = etaCell[ncell] - eta;
1072 dphi = phiCell[ncell] - phi;
1073 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1074 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1075 dr = TMath::Sqrt(deta * deta + dphi * dphi);
1076 if(dr <= rc){ // cell in cone
1077 flagCell[ncell] = -1;
1078 etCone+=etCell[ncell];
1079 nCellIn++;
1080 }
1081 else nCellOut++;
1082 }
1083
1084 // Select jets with et > background
1085 // estimate max fluctuation of background in cone
1086 Double_t ncellin = (Double_t)nCellIn;
1087 Double_t ntcell = (Double_t)nCell;
1088 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
1089 // min cone et
1090 Double_t etcmin = etCone ; // could be used etCone - etmin !!
1091 //decisions !! etbmax < etcmin
1092
1093 for(Int_t mcell =0; mcell < nCell; mcell++){
1094 if(flagCell[mcell] == -1){
1095 if(etbmax < etcmin)
1096 flagCell[mcell] = 1; //flag cell as used
1097 else
1098 flagCell[mcell] = 0; // leave it free
1099 }
1100 }
1101 //store tmp jet info !!!
1102
1103 if(etbmax < etcmin) {
1104 etaAlgoJet[nJets] = eta;
1105 phiAlgoJet[nJets] = phi;
1106 etAlgoJet[nJets] = etCone;
1107 ncellsAlgoJet[nJets] = nCellIn;
1108 nJets++;
1109 }
1110
1111 } // end of cells loop
ee7de0dd 1112
1113 //reorder jets by et in cone
1114 //sort jets by energy
1115 Int_t * idx = new Int_t[nJets];
1116 TMath::Sort(nJets, etAlgoJet, idx);
8838ab7a 1117 for(Int_t p = 0; p < nJets; p++)
1118 {
1119 etaJet[p] = etaAlgoJet[idx[p]];
1120 phiJet[p] = phiAlgoJet[idx[p]];
1121 etJet[p] = etAlgoJet[idx[p]];
1122 etallJet[p] = etAlgoJet[idx[p]];
1123 ncellsJet[p] = ncellsAlgoJet[idx[p]];
1124 }
1125
ee7de0dd 1126 //delete
1127 delete index;
1128 delete idx;
1129
1130}
ee7de0dd 1131
8838ab7a 1132////////////////////////////////////////////////////////////////////////
be6e5811 1133void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT,
1134 const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1135 const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet,
1136 const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet)
ee7de0dd 1137{
8838ab7a 1138 //
1139 // Background subtraction using cone method but without correction in dE/deta distribution
1140 // Cases to take into account the EMCal geometry are included
1141 //
1142
ee7de0dd 1143 //calculate energy inside and outside cones
1144 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
be6e5811 1145 Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
1146 fDebug = fReader->GetReaderHeader()->GetDebug();
ee7de0dd 1147 Float_t rc= header->GetRadius();
1148 Float_t etIn[30];
1149 Float_t etOut = 0;
8838ab7a 1150
1151 for(Int_t j=0;j<30;j++){etIn[j]=0.;}
1152
ee7de0dd 1153 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
8838ab7a 1154
1155 for(Int_t ijet=0; ijet<nJ; ijet++){
1156
1157 Float_t deta = etaT[jpart] - etaJet[ijet];
1158 Float_t dphi = phiT[jpart] - phiJet[ijet];
1159 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1160 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1161
1162 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1163 if(dr <= rc){ // particles inside this cone
1164 multJet[ijet]+=vectT[jpart];
1165 injet[jpart] = ijet;
1166
1167 if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut
1168 etIn[ijet] += ptT[jpart];
1169 if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart];
1170 }
1171 break;
1172 }
1173 }// end jets loop
1174
1175 if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[jpart] == 1)){
1176 etOut += ptT[jpart]; // particle outside cones and pt cut
1177 }
ee7de0dd 1178 } //end particle loop
1179
1180 //estimate jets and background areas
8838ab7a 1181 // TPC case
1182 if(fOpt == 0 || fOpt == 1){
1183 Float_t areaJet[30];
1184 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1185
1186 for(Int_t k=0; k<nJ; k++){
ee7de0dd 1187 Float_t detamax = etaJet[k] + rc;
1188 Float_t detamin = etaJet[k] - rc;
1189 Float_t accmax = 0.0; Float_t accmin = 0.0;
1190 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
8838ab7a 1191 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1192 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
ee7de0dd 1193 }
1194 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
8838ab7a 1195 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1196 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
ee7de0dd 1197 }
1198 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1199 areaOut = areaOut - areaJet[k];
8838ab7a 1200 }
1201 //subtract background using area method
1202 for(Int_t ljet=0; ljet<nJ; ljet++){
1203 Float_t areaRatio = areaJet[ljet]/areaOut;
1204 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1205 }
1206
1207 // estimate new total background
1208 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1209 etbgTotalN = etOut*areaT/areaOut;
1210 }
1211 else { // If EMCal included
1212 Float_t areaJet[30];
1213 Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
1214 for(Int_t k=0; k<nJ; k++){
1215 Float_t detamax = etaJet[k] + rc;
1216 Float_t detamin = etaJet[k] - rc;
1217 Float_t dphimax = phiJet[k] + rc;
1218 Float_t dphimin = phiJet[k] - rc;
1219 Float_t eMax = header->GetLegoEtaMax();
1220 Float_t eMin = header->GetLegoEtaMin();
1221 Float_t pMax = header->GetLegoPhiMax();
1222 Float_t pMin = header->GetLegoPhiMin();
1223 Float_t accetamax = 0.0; Float_t accetamin = 0.0;
1224 Float_t accphimax = 0.0; Float_t accphimin = 0.0;
1225 if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1226 (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1227 Float_t h = eMax - etaJet[k];
1228 accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1229 }
1230 if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1231 (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1232 Float_t h = eMax + etaJet[k];
1233 accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1234 }
1235 if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )||
1236 (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){
1237 Float_t h = pMax - phiJet[k];
1238 accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1239 }
1240 if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )||
1241 (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){
1242 Float_t h = phiJet[k] - pMin;
1243 accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1244 }
1245
1246 if(detamax > eMax && dphimax > pMax ){
1247 Float_t he = eMax - etaJet[k];
1248 Float_t hp = pMax - phiJet[k];
1249 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1250 Float_t alphae = TMath::ACos(he/rc);
1251 Float_t alphap = TMath::ACos(hp/rc);
1252 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1253 if(rlim <= rc){
1254 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1255 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1256 }
1257 if(rlim > rc){
1258 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1259 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1260 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1261 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1262 }
1263 }
1264
1265 if(detamax > eMax && dphimin < pMin ){
1266 Float_t he = eMax - etaJet[k];
1267 Float_t hp = phiJet[k] - pMin;
1268 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1269 Float_t alphae = TMath::ACos(he/rc);
1270 Float_t alphap = TMath::ACos(hp/rc);
1271 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1272 if(rlim <= rc){
1273 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1274 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1275 }
1276 if(rlim > rc){
1277 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1278 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1279 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1280 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1281 }
1282 }
1283
1284 if(detamin < eMin && dphimax > pMax ){
1285 Float_t he = eMax + etaJet[k];
1286 Float_t hp = pMax - phiJet[k];
1287 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1288 Float_t alphae = TMath::ACos(he/rc);
1289 Float_t alphap = TMath::ACos(hp/rc);
1290 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1291 if(rlim <= rc){
1292 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1293 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1294 }
1295 if(rlim > rc){
1296 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1297 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1298 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1299 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1300 }
1301 }
1302
1303 if(detamin < eMin && dphimin < pMin ){
1304 Float_t he = eMax + etaJet[k];
1305 Float_t hp = phiJet[k] - pMin;
1306 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1307 Float_t alphae = TMath::ACos(he/rc);
1308 Float_t alphap = TMath::ACos(hp/rc);
1309 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1310 if(rlim <= rc){
1311 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1312 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1313 }
1314 if(rlim > rc){
1315 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1316 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1317 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1318 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1319 }
1320 }
1321 areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin;
1322 areaOut = areaOut - areaJet[k];
1323 } // end loop on jets
1324
1325 //subtract background using area method
1326 for(Int_t ljet=0; ljet<nJ; ljet++){
1327 Float_t areaRatio = areaJet[ljet]/areaOut;
1328 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1329 }
1330
1331 // estimate new total background
1332 Float_t areaT = 2*(header->GetLegoEtaMax()*header->GetLegoPhiMax());
1333 etbgTotalN = etOut*areaT/areaOut;
1334 }
1335
1336}
1337
1338////////////////////////////////////////////////////////////////////////
be6e5811 1339void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
1340 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
1341 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
8838ab7a 1342 Int_t* multJet, Int_t* injet)
1343{
1344 //background subtraction using cone method but without correction in dE/deta distribution
1345
1346 //calculate energy inside and outside cones
1347 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1348 Float_t rc= header->GetRadius();
1349 Float_t etIn[30];
1350 Float_t etOut = 0;
1351 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1352 // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
1353 for(Int_t ijet=0; ijet<nJ; ijet++){
1354 Float_t deta = etaT[jpart] - etaJet[ijet];
1355 Float_t dphi = phiT[jpart] - phiJet[ijet];
1356 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1357 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1358 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1359 if(dr <= rc){ // particles inside this cone
1360 multJet[ijet]++;
1361 injet[jpart] = ijet;
1362 if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
1363 etIn[ijet] += ptT[jpart];
1364 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
1365 }
1366 break;
1367 }
1368 }// end jets loop
1369 if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
1370 etOut += ptT[jpart]; // particle outside cones and pt cut
1371 } //end particle loop
1372
1373 //estimate jets and background areas
1374 Float_t areaJet[30];
1375 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1376 for(Int_t k=0; k<nJ; k++){
1377 Float_t detamax = etaJet[k] + rc;
1378 Float_t detamin = etaJet[k] - rc;
1379 Float_t accmax = 0.0; Float_t accmin = 0.0;
1380 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1381 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1382 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1383 }
1384 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1385 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1386 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1387 }
1388 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1389 areaOut = areaOut - areaJet[k];
ee7de0dd 1390 }
1391 //subtract background using area method
1392 for(Int_t ljet=0; ljet<nJ; ljet++){
8838ab7a 1393 Float_t areaRatio = areaJet[ljet]/areaOut;
1394 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
ee7de0dd 1395 }
8838ab7a 1396
ee7de0dd 1397 // estimate new total background
1398 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1399 etbgTotalN = etOut*areaT/areaOut;
8838ab7a 1400
ee7de0dd 1401}
1402
ee7de0dd 1403
8838ab7a 1404////////////////////////////////////////////////////////////////////////
be6e5811 1405void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
1406 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1407 const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
1408 Float_t* etsigJet, Int_t* multJet, Int_t* injet)
ee7de0dd 1409{
1410
1411 //background subtraction using statistical method
1412 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1413 Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
8838ab7a 1414
ee7de0dd 1415 //calculate energy inside
1416 Float_t rc= header->GetRadius();
1417 Float_t etIn[30];
8838ab7a 1418
1419 for(Int_t jpart = 0; jpart < nIn; jpart++)
1420 { // loop for all particles in array
1421
1422 for(Int_t ijet=0; ijet<nJ; ijet++)
1423 {
1424 Float_t deta = etaT[jpart] - etaJet[ijet];
1425 Float_t dphi = phiT[jpart] - phiJet[ijet];
1426 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1427 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1428 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1429 if(dr <= rc){ // particles inside this cone
1430 multJet[ijet]++;
1431 injet[jpart] = ijet;
1432 if(cFlagT[jpart] == 1){ // pt cut
1433 etIn[ijet]+= ptT[jpart];
1434 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1435 }
1436 break;
1437 }
1438 }// end jets loop
1439 } //end particle loop
1440
ee7de0dd 1441 //calc jets areas
1442 Float_t areaJet[30];
1443 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
8838ab7a 1444 for(Int_t k=0; k<nJ; k++)
1445 {
ee7de0dd 1446 Float_t detamax = etaJet[k] + rc;
1447 Float_t detamin = etaJet[k] - rc;
1448 Float_t accmax = 0.0; Float_t accmin = 0.0;
1449 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
8838ab7a 1450 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1451 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
ee7de0dd 1452 }
1453 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
8838ab7a 1454 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1455 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
ee7de0dd 1456 }
1457 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
8838ab7a 1458 }
ee7de0dd 1459
1460 //subtract background using area method
1461 for(Int_t ljet=0; ljet<nJ; ljet++){
8838ab7a 1462 Float_t areaRatio = areaJet[ljet]/areaOut;
1463 etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
ee7de0dd 1464 }
8838ab7a 1465
ee7de0dd 1466 etbgTotalN = etbgStat;
ee7de0dd 1467}
1468
1469////////////////////////////////////////////////////////////////////////
be6e5811 1470void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
1471 Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1472 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
8838ab7a 1473 Int_t* multJet, Int_t* injet)
ee7de0dd 1474{
8838ab7a 1475 // Cone background subtraction method taking into acount dEt/deta distribution
1476 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1477 //general
1478 Float_t rc= header->GetRadius();
1479 Float_t etamax = header->GetLegoEtaMax();
1480 Float_t etamin = header->GetLegoEtaMin();
1481 Int_t ndiv = 100;
1482
1483 // jet energy and area arrays
1484 TH1F* hEtJet[30];
1485 TH1F* hAreaJet[30];
1486 for(Int_t mjet=0; mjet<nJ; mjet++){
1487 char hEtname[256]; char hAreaname[256];
1488 sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1489 hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
1490 hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
ee7de0dd 1491 }
8838ab7a 1492 // background energy and area
1493 TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
1494 TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
1495
1496 //fill energies
1497 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1498 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1499 Float_t deta = etaT[jpart] - etaJet[ijet];
1500 Float_t dphi = phiT[jpart] - phiJet[ijet];
1501 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1502 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1503 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1504 if(dr <= rc){ // particles inside this cone
1505 injet[jpart] = ijet;
1506 multJet[ijet]++;
1507 if(cFlagT[jpart] == 1){// pt cut
1508 hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
1509 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1510 }
1511 break;
1512 }
1513 }// end jets loop
1514
1515 if(injet[jpart] == -1 && cFlagT[jpart] == 1)
1516 hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
ee7de0dd 1517 } //end particle loop
1518
8838ab7a 1519 //calc areas
1520 Float_t eta0 = etamin;
1521 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1522 Float_t eta1 = eta0 + etaw;
1523 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1524 Float_t etac = eta0 + etaw/2.0;
1525 Float_t areabg = etaw*2.0*TMath::Pi();
1526 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1527 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1528 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1529 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1530 Float_t areaj = 0.0;
1531 if(deta0 > rc && deta1 < rc){
1532 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1533 areaj = acc1;
1534 }
1535 if(deta0 < rc && deta1 > rc){
1536 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1537 areaj = acc0;
1538 }
1539 if(deta0 < rc && deta1 < rc){
1540 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1541 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1542 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
1543 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1544 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1545 }
1546 hAreaJet[ijet]->Fill(etac,areaj);
1547 areabg = areabg - areaj;
1548 } // end jets loop
1549 hAreaBackg->Fill(etac,areabg);
1550 eta0 = eta1;
1551 eta1 = eta1 + etaw;
1552 } // end loop for all eta bins
1553
1554 //subtract background
1555 for(Int_t kjet=0; kjet<nJ; kjet++){
1556 etJet[kjet] = 0.0; // first clear etJet for this jet
1557 for(Int_t bin = 0; bin< ndiv; bin++){
1558 if(hAreaJet[kjet]->GetBinContent(bin)){
1559 Float_t areab = hAreaBackg->GetBinContent(bin);
1560 Float_t etb = hEtBackg->GetBinContent(bin);
1561 Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1562 etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
1563 }
1564 }
1565 }
ee7de0dd 1566
8838ab7a 1567 // calc background total
1568 Double_t etOut = hEtBackg->Integral();
1569 Double_t areaOut = hAreaBackg->Integral();
1570 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1571 etbgTotalN = etOut*areaT/areaOut;
1572
1573 //delete
1574 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
1575 delete hEtJet[ljet];
1576 delete hAreaJet[ljet];
1577 }
ee7de0dd 1578
8838ab7a 1579 delete hEtBackg;
1580 delete hAreaBackg;
1581}
ee7de0dd 1582
8838ab7a 1583////////////////////////////////////////////////////////////////////////
be6e5811 1584void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
1585 Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1586 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
8838ab7a 1587 Int_t* multJet, Int_t* injet)
ee7de0dd 1588{
8838ab7a 1589 // Ratio background subtraction method taking into acount dEt/deta distribution
1590 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1591 //factor F calc before
1592 Float_t bgRatioCut = header->GetBackgCutRatio();
1593
1594 //general
1595 Float_t rc= header->GetRadius();
1596 Float_t etamax = header->GetLegoEtaMax();
1597 Float_t etamin = header->GetLegoEtaMin();
1598 Int_t ndiv = 100;
1599
1600 // jet energy and area arrays
1601 TH1F* hEtJet[30];
1602 TH1F* hAreaJet[30];
1603 for(Int_t mjet=0; mjet<nJ; mjet++){
1604 char hEtname[256]; char hAreaname[256];
1605 sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1606 hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax); // change range
1607 hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax); // change range
ee7de0dd 1608 }
8838ab7a 1609 // background energy and area
1610 TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range
1611 TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range
1612
1613 //fill energies
1614 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1615 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1616 Float_t deta = etaT[jpart] - etaJet[ijet];
1617 Float_t dphi = phiT[jpart] - phiJet[ijet];
1618 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1619 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1620 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1621 if(dr <= rc){ // particles inside this cone
1622 multJet[ijet]++;
1623 injet[jpart] = ijet;
1624 if(cFlagT[jpart] == 1){ //pt cut
1625 hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
1626 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1627 }
1628 break;
1629 }
1630 }// end jets loop
1631 if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
ee7de0dd 1632 } //end particle loop
1633
8838ab7a 1634 //calc areas
1635 Float_t eta0 = etamin;
1636 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1637 Float_t eta1 = eta0 + etaw;
1638 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1639 Float_t etac = eta0 + etaw/2.0;
1640 Float_t areabg = etaw*2.0*TMath::Pi();
1641 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1642 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1643 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1644 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1645 Float_t areaj = 0.0;
1646 if(deta0 > rc && deta1 < rc){
1647 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1648 areaj = acc1;
1649 }
1650 if(deta0 < rc && deta1 > rc){
1651 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1652 areaj = acc0;
1653 }
1654 if(deta0 < rc && deta1 < rc){
1655 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1656 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1657 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
1658 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1659 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1660 }
1661 hAreaJet[ijet]->Fill(etac,areaj);
1662 areabg = areabg - areaj;
1663 } // end jets loop
1664 hAreaBackg->Fill(etac,areabg);
1665 eta0 = eta1;
1666 eta1 = eta1 + etaw;
1667 } // end loop for all eta bins
1668
1669 //subtract background
1670 for(Int_t kjet=0; kjet<nJ; kjet++){
1671 etJet[kjet] = 0.0; // first clear etJet for this jet
1672 for(Int_t bin = 0; bin< ndiv; bin++){
1673 if(hAreaJet[kjet]->GetBinContent(bin)){
1674 Float_t areab = hAreaBackg->GetBinContent(bin);
1675 Float_t etb = hEtBackg->GetBinContent(bin);
1676 Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1677 etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
1678 }
1679 }
1680 }
1681
1682 // calc background total
1683 Double_t etOut = hEtBackg->Integral();
1684 Double_t areaOut = hAreaBackg->Integral();
1685 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1686 etbgTotalN = etOut*areaT/areaOut;
1687
1688 //delete
1689 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
1690 delete hEtJet[ljet];
1691 delete hAreaJet[ljet];
1692 }
1693
1694 delete hEtBackg;
1695 delete hAreaBackg;
ee7de0dd 1696}
1697
1698////////////////////////////////////////////////////////////////////////
ee7de0dd 1699void AliUA1JetFinderV2::Reset()
1700{
1701 fLego->Reset();
1702 fJets->ClearJets();
1703 AliJetFinder::Reset();
1704}
1705
1706////////////////////////////////////////////////////////////////////////
ee7de0dd 1707void AliUA1JetFinderV2::WriteJHeaderToFile()
1708{
1709 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1710 header->Write();
1711}
1712
1713////////////////////////////////////////////////////////////////////////
8838ab7a 1714void AliUA1JetFinderV2::InitTask(TChain* tree)
ee7de0dd 1715{
8838ab7a 1716
ee7de0dd 1717 // initializes some variables
1718 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
8838ab7a 1719 // book lego
1720 fLego = new TH2F("legoH","eta-phi",
1721 header->GetLegoNbinEta(), header->GetLegoEtaMin(),
1722 header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
1723 header->GetLegoPhiMin(), header->GetLegoPhiMax());
1724
ee7de0dd 1725 fDebug = fReader->GetReaderHeader()->GetDebug();
1726 fOpt = fReader->GetReaderHeader()->GetDetector();
8838ab7a 1727
1728 // Tasks initialization
ee7de0dd 1729 if(fOpt>0)
8838ab7a 1730 fReader->CreateTasks(tree);
ee7de0dd 1731
1732}