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