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