]>
Commit | Line | Data |
---|---|---|
50c7d9f7 | 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 | ||
139cbd96 | 16 | /* $Id$ */ |
17 | ||
18 | //-------------------------------------------------- | |
19 | // Method implementation for background studies and background subtraction with UA1 algorithms | |
20 | // | |
21 | // Author: magali.estienne@subatech.in2p3.fr | |
22 | //------------------------------------------------- | |
23 | ||
50c7d9f7 | 24 | #include <Riostream.h> |
50c7d9f7 | 25 | #include <TList.h> |
50c7d9f7 | 26 | #include <TClonesArray.h> |
139cbd96 | 27 | #include <TH1F.h> |
28 | #include <TH2F.h> | |
50c7d9f7 | 29 | |
139cbd96 | 30 | #include "AliAODJetEventBackground.h" |
31 | #include "AliUA1JetHeaderV1.h" | |
32 | #include "AliJetCalTrk.h" | |
50c7d9f7 | 33 | #include "AliJetBkg.h" |
34 | ||
139cbd96 | 35 | |
36 | using namespace std; | |
37 | ||
50c7d9f7 | 38 | ClassImp(AliJetBkg) |
39 | ||
40 | //////////////////////////////////////////////////////////////////////// | |
41 | ||
42 | AliJetBkg::AliJetBkg(): | |
856618e7 | 43 | TObject(), |
139cbd96 | 44 | fEvent(0x0), |
45 | fHeader(0x0), | |
46 | fDebug(0), | |
47 | fhEtBackg(0x0), | |
48 | fhAreaBackg(0x0) | |
50c7d9f7 | 49 | { |
50 | // Default constructor | |
139cbd96 | 51 | for(int i = 0;i < kMaxJets;i++){ |
52 | fhAreaJet[i] = fhEtJet[i] = 0; | |
53 | } | |
50c7d9f7 | 54 | } |
139cbd96 | 55 | |
56 | //---------------------------------------------------------------- | |
50c7d9f7 | 57 | AliJetBkg::AliJetBkg(const AliJetBkg& input): |
856618e7 | 58 | TObject(input), |
139cbd96 | 59 | fEvent(input.fEvent), |
60 | fHeader(input.fHeader), | |
61 | fDebug(input.fDebug), | |
62 | fhEtBackg(input.fhEtBackg), | |
63 | fhAreaBackg(input.fhAreaBackg) | |
50c7d9f7 | 64 | { |
65 | // copy constructor | |
139cbd96 | 66 | for(int i = 0;i < kMaxJets;i++){ |
67 | fhAreaJet[i] = input.fhAreaJet[i]; | |
68 | fhEtJet[i] = input.fhEtJet[i]; | |
69 | } | |
50c7d9f7 | 70 | |
71 | } | |
294056ab | 72 | |
139cbd96 | 73 | //---------------------------------------------------------------- |
74 | AliJetBkg::~AliJetBkg() | |
75 | { | |
76 | // Destructor | |
77 | if(fhEtBackg) delete fhEtBackg; | |
78 | if(fhAreaBackg) delete fhAreaBackg; | |
79 | for(int i = 0;i < kMaxJets;i++){ | |
80 | if(fhAreaJet[i]) delete fhAreaJet[i]; | |
81 | if(fhEtJet[i]) delete fhEtJet[i]; | |
82 | } | |
294056ab | 83 | |
139cbd96 | 84 | } |
674c4c04 | 85 | |
139cbd96 | 86 | //---------------------------------------------------------------- |
87 | Bool_t AliJetBkg::PtCutPass(Int_t id, Int_t nTracks) | |
88 | { | |
89 | // Check if track or cell passes the cut flag | |
90 | if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetCutFlag() == 1) | |
91 | return kTRUE; | |
92 | else return kFALSE; | |
294056ab | 93 | |
139cbd96 | 94 | } |
674c4c04 | 95 | |
139cbd96 | 96 | //---------------------------------------------------------------- |
97 | Bool_t AliJetBkg::SignalCutPass(Int_t id, Int_t nTracks) | |
98 | { | |
99 | // Check if track or cell passes the cut flag | |
100 | if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetSignalFlag() == 1) | |
101 | return kTRUE; | |
102 | else return kFALSE; | |
294056ab | 103 | |
294056ab | 104 | } |
294056ab | 105 | |
139cbd96 | 106 | //---------------------------------------------------------------- |
107 | Float_t AliJetBkg::CalcJetAreaEtaCut(Float_t radius, const Float_t etaJet) | |
108 | { | |
109 | // Calculate jet area taking into account an acceptance cut in eta | |
110 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; | |
111 | Float_t detamax = etaJet + radius; | |
112 | Float_t detamin = etaJet - radius; | |
113 | Float_t accmax = 0.0; Float_t accmin = 0.0; | |
114 | if(detamax > header->GetLegoEtaMax()){ // sector outside etamax | |
115 | Float_t h = header->GetLegoEtaMax() - etaJet; | |
116 | accmax = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h); | |
294056ab | 117 | } |
139cbd96 | 118 | if(detamin < header->GetLegoEtaMin()){ // sector outside etamin |
119 | Float_t h = header->GetLegoEtaMax() + etaJet; | |
120 | accmin = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h); | |
294056ab | 121 | } |
294056ab | 122 | |
139cbd96 | 123 | return radius*radius*TMath::Pi() - accmax - accmin; |
294056ab | 124 | |
139cbd96 | 125 | } |
294056ab | 126 | |
139cbd96 | 127 | //---------------------------------------------------------------- |
128 | void AliJetBkg::CalcJetAndBckgAreaEtaCut(Bool_t calcOutsideArea, Float_t radius, const Int_t nJets, const Float_t* etaJet, Float_t* &areaJet, Float_t &areaOut) | |
129 | { | |
130 | // Calculate jet and bacground areas taking into account an acceptance cut in eta | |
294056ab | 131 | |
139cbd96 | 132 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; |
133 | areaOut = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax() - header->GetLegoPhiMin()); | |
134 | for(Int_t k=0; k<nJets; k++){ | |
135 | areaJet[k] = CalcJetAreaEtaCut(radius, etaJet[k]); | |
136 | if(calcOutsideArea) areaOut = areaOut - areaJet[k]; | |
294056ab | 137 | } |
138 | ||
50c7d9f7 | 139 | } |
50c7d9f7 | 140 | |
139cbd96 | 141 | //---------------------------------------------------------------- |
142 | void AliJetBkg::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, Float_t&sigmaN, | |
143 | const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, | |
144 | Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, | |
145 | Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet, | |
146 | Int_t* injet, Float_t* &areaJet) | |
147 | { | |
148 | // | |
149 | // Background subtraction using cone method but without correction in dE/deta distribution | |
150 | // Cases to take into account the EMCal geometry are included | |
151 | // | |
50c7d9f7 | 152 | |
139cbd96 | 153 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; |
154 | //calculate energy inside and outside cones | |
155 | fDebug = header->GetDebug(); | |
156 | Float_t rc = header->GetRadius(); | |
157 | Float_t etOut = 0; | |
158 | // Get number of tracks from EventCalTrk | |
159 | Int_t nTracks = fEvent->GetNCalTrkTracks(); | |
160 | ||
161 | Float_t etIn[kMaxJets] = {0}; | |
162 | Float_t areaOut = 0.; | |
163 | ||
164 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
165 | ||
166 | for(Int_t ijet=0; ijet<nJ; ijet++){ | |
167 | Float_t deta = etaT[jpart] - etaJet[ijet]; | |
168 | Float_t dphi = phiT[jpart] - phiJet[ijet]; | |
169 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
170 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
171 | Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
172 | if(dr <= rc){ // particles inside this cone | |
173 | if(jpart < nTracks) multJetT[ijet]++; | |
174 | else multJetC[ijet]++; | |
175 | multJet[ijet]++; | |
176 | injet[jpart] = ijet; | |
177 | if(PtCutPass(jpart,nTracks)){ // pt cut | |
178 | etIn[ijet] += ptT[jpart]; | |
179 | if(SignalCutPass(jpart,nTracks)) | |
180 | etsigJet[ijet]+= ptT[jpart]; | |
181 | } | |
182 | break; | |
183 | } | |
184 | }// end jets loop | |
185 | ||
186 | if((injet[jpart] == -1) && | |
187 | (PtCutPass(jpart,nTracks))){ | |
188 | etOut += ptT[jpart]; // particle outside cones and pt cut | |
189 | } | |
190 | } //end particle loop | |
191 | ||
192 | // Calculate jet and background areas | |
193 | Bool_t calcAreaOut = kTRUE; | |
194 | CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut); | |
195 | ||
196 | //subtract background using area method | |
197 | for(Int_t ljet=0; ljet<nJ; ljet++){ | |
198 | Float_t areaRatio = areaJet[ljet]/areaOut; | |
199 | etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction | |
200 | } | |
50c7d9f7 | 201 | |
139cbd96 | 202 | // estimate new total background |
203 | Float_t areaT = 0; | |
204 | areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); | |
205 | etbgTotalN = etOut*areaT/areaOut; | |
206 | ||
207 | // estimate standard deviation of background | |
208 | Int_t count = 0; | |
209 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
210 | if((injet[jpart] == -1) && | |
211 | (PtCutPass(jpart,nTracks))){ | |
212 | sigmaN += etbgTotalN/areaT - ptT[jpart]; | |
213 | // To be checked (Division by jet area to obtain standard deviation of rho ?) | |
214 | ||
215 | count=count+1; | |
216 | } | |
217 | } | |
218 | if (count>0) | |
219 | sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); | |
50c7d9f7 | 220 | |
50c7d9f7 | 221 | } |
222 | ||
139cbd96 | 223 | //---------------------------------------------------------------- |
224 | void AliJetBkg::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN, Float_t&sigmaN, | |
225 | const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, | |
226 | Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, | |
227 | Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet, | |
228 | Int_t* injet, Float_t* &areaJet) | |
50c7d9f7 | 229 | { |
139cbd96 | 230 | // |
50c7d9f7 | 231 | //background subtraction using statistical method |
139cbd96 | 232 | // Cases to take into account the EMCal geometry are included |
233 | // | |
50c7d9f7 | 234 | |
139cbd96 | 235 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; |
236 | Float_t etbgStat = header->GetBackgStat(); // pre-calculated background | |
50c7d9f7 | 237 | |
139cbd96 | 238 | //calculate energy inside |
239 | Float_t rc= header->GetRadius(); | |
240 | Float_t etIn[kMaxJets] = {0.0}; | |
241 | // Get number of tracks from EventCalTrk | |
242 | Int_t nTracks = fEvent->GetNCalTrkTracks(); | |
243 | Float_t areaOut = 0.; | |
244 | ||
245 | for(Int_t jpart = 0; jpart < nIn; jpart++) | |
246 | { // loop for all particles in array | |
247 | ||
248 | for(Int_t ijet=0; ijet<nJ; ijet++) | |
249 | { | |
250 | Float_t deta = etaT[jpart] - etaJet[ijet]; | |
251 | Float_t dphi = phiT[jpart] - phiJet[ijet]; | |
252 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
253 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
254 | Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
255 | if(dr <= rc){ // particles inside this cone | |
256 | if(jpart < nTracks) multJetT[ijet]++; | |
257 | else multJetC[ijet]++; | |
258 | multJet[ijet]++; | |
259 | injet[jpart] = ijet; | |
260 | ||
261 | if(PtCutPass(jpart,nTracks)){ // pt cut | |
262 | etIn[ijet] += ptT[jpart]; | |
263 | if(SignalCutPass(jpart,nTracks)) | |
264 | etsigJet[ijet]+= ptT[jpart]; | |
265 | } | |
266 | break; | |
267 | } | |
268 | }// end jets loop | |
269 | } //end particle loop | |
50c7d9f7 | 270 | |
139cbd96 | 271 | // Calculate jet and background areas |
272 | Bool_t calcAreaOut = kFALSE; | |
273 | CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut); | |
274 | ||
275 | //subtract background using area method | |
276 | for(Int_t ljet=0; ljet<nJ; ljet++){ | |
277 | Float_t areaRatio = areaJet[ljet]/areaOut; | |
278 | etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction | |
50c7d9f7 | 279 | } |
139cbd96 | 280 | Int_t count=0; |
281 | etbgTotalN = etbgStat; | |
282 | ||
283 | // estimate standard deviation of background | |
284 | Float_t areaT = 0; | |
285 | areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); | |
286 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
287 | if((injet[jpart] == -1) && | |
288 | (PtCutPass(jpart,nTracks))){ | |
289 | sigmaN += etbgTotalN/areaT - ptT[jpart]; | |
290 | count=count+1; | |
291 | } | |
50c7d9f7 | 292 | } |
139cbd96 | 293 | if(count>0)sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); |
50c7d9f7 | 294 | |
50c7d9f7 | 295 | } |
296 | ||
139cbd96 | 297 | //---------------------------------------------------------------- |
298 | void AliJetBkg::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN, | |
299 | const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, Float_t* etJet, | |
300 | const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, | |
301 | Int_t* multJetT, Int_t* multJetC, Int_t* multJet, Int_t* injet, Float_t* &/*areaJet*/) | |
302 | { | |
303 | // | |
304 | // Cone background subtraction method taking into acount dEt/deta distribution | |
305 | // Cases to take into account the EMCal geometry are not included | |
306 | // | |
307 | ||
308 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; | |
309 | //general | |
310 | Float_t rc= header->GetRadius(); | |
311 | Float_t etamax = header->GetLegoEtaMax(); | |
312 | Float_t etamin = header->GetLegoEtaMin(); | |
313 | Int_t ndiv = 100; | |
314 | // Get number of tracks from EventCalTrk | |
315 | Int_t nTracks = fEvent->GetNCalTrkTracks(); | |
316 | ||
317 | // jet energy and area arrays | |
318 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
319 | TH1::AddDirectory(kFALSE); | |
320 | ||
321 | for(Int_t mjet=0; mjet<nJ; mjet++){ | |
322 | if(!fhEtJet[mjet]){ | |
323 | fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax); | |
324 | } | |
325 | if(!fhAreaJet[mjet]){ | |
326 | fhAreaJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax); | |
327 | } | |
328 | fhEtJet[mjet]->Reset(); | |
329 | fhAreaJet[mjet]->Reset(); | |
50c7d9f7 | 330 | } |
139cbd96 | 331 | // background energy and area |
332 | if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); | |
333 | fhEtBackg->Reset(); | |
334 | if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); | |
335 | fhAreaBackg->Reset(); | |
336 | TH1::AddDirectory(oldStatus); | |
337 | ||
338 | //fill energies | |
339 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
340 | for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets | |
341 | Float_t deta = etaT[jpart] - etaJet[ijet]; | |
342 | Float_t dphi = phiT[jpart] - phiJet[ijet]; | |
343 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
344 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
345 | Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
346 | if(dr <= rc){ // particles inside this cone | |
347 | if(jpart < nTracks) multJetT[ijet]++; | |
348 | else multJetC[ijet]++; | |
349 | multJet[ijet]++; | |
350 | injet[jpart] = ijet; | |
351 | ||
352 | if(PtCutPass(jpart,nTracks)){ // pt cut | |
353 | fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone | |
354 | if(SignalCutPass(jpart,nTracks)) | |
355 | etsigJet[ijet]+= ptT[jpart]; | |
356 | } | |
357 | break; | |
358 | } | |
359 | }// end jets loop | |
360 | ||
361 | if((injet[jpart] == -1) && | |
362 | (PtCutPass(jpart,nTracks) == 1)) | |
363 | fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones | |
364 | } //end particle loop | |
365 | ||
366 | //calc areas | |
367 | Float_t eta0 = etamin; | |
368 | Float_t etaw = (etamax - etamin)/((Float_t)ndiv); | |
369 | Float_t eta1 = eta0 + etaw; | |
370 | for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins | |
371 | Float_t etac = eta0 + etaw/2.0; | |
372 | Float_t areabg = etaw*2.0*TMath::Pi(); | |
373 | for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets | |
374 | Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]); | |
375 | Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]); | |
376 | Float_t acc0 = 0.0; Float_t acc1 = 0.0; | |
377 | Float_t areaj = 0.0; | |
378 | if(deta0 > rc && deta1 < rc){ | |
379 | acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); | |
380 | areaj = acc1; | |
381 | } | |
382 | if(deta0 < rc && deta1 > rc){ | |
383 | acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); | |
384 | areaj = acc0; | |
385 | } | |
386 | if(deta0 < rc && deta1 < rc){ | |
387 | acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); | |
388 | acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); | |
389 | if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1 | |
390 | if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2 | |
391 | if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3 | |
392 | } | |
393 | fhAreaJet[ijet]->Fill(etac,areaj); | |
394 | areabg = areabg - areaj; | |
395 | } // end jets loop | |
396 | fhAreaBackg->Fill(etac,areabg); | |
397 | eta0 = eta1; | |
398 | eta1 = eta1 + etaw; | |
399 | } // end loop for all eta bins | |
400 | ||
401 | //subtract background | |
402 | for(Int_t kjet=0; kjet<nJ; kjet++){ | |
403 | etJet[kjet] = 0.0; // first clear etJet for this jet | |
404 | for(Int_t bin = 0; bin< ndiv; bin++){ | |
405 | if(fhAreaJet[kjet]->GetBinContent(bin)){ | |
406 | Float_t areab = fhAreaBackg->GetBinContent(bin); | |
407 | Float_t etb = fhEtBackg->GetBinContent(bin); | |
408 | Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab; | |
409 | etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction | |
db2272d5 | 410 | } |
139cbd96 | 411 | } |
db2272d5 | 412 | } |
50c7d9f7 | 413 | |
139cbd96 | 414 | // calc background total |
415 | Double_t etOut = fhEtBackg->Integral(); | |
416 | Double_t areaOut = fhAreaBackg->Integral(); | |
417 | Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); | |
418 | etbgTotalN = etOut*areaT/areaOut; | |
419 | ||
420 | Int_t count=0; | |
421 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
422 | if((injet[jpart] == -1) && | |
423 | (PtCutPass(jpart,nTracks))){ | |
424 | sigmaN += etbgTotalN/areaT - ptT[jpart]; | |
425 | count=count+1; | |
426 | } | |
50c7d9f7 | 427 | } |
139cbd96 | 428 | sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); |
429 | ||
50c7d9f7 | 430 | } |
431 | ||
139cbd96 | 432 | //---------------------------------------------------------------- |
433 | void AliJetBkg::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN, | |
434 | const Float_t* ptT,const Float_t* etaT, const Float_t* phiT, | |
435 | Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, | |
436 | Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet, | |
437 | Int_t* injet, Float_t* &/*areaJet*/) | |
50c7d9f7 | 438 | { |
139cbd96 | 439 | // Ratio background subtraction method taking into acount dEt/deta distribution |
440 | // Cases to take into account the EMCal geometry are not included | |
50c7d9f7 | 441 | |
139cbd96 | 442 | AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; |
443 | //factor F calc before | |
444 | Float_t bgRatioCut = header->GetBackgCutRatio(); | |
445 | ||
446 | //general | |
447 | Float_t rc= header->GetRadius(); | |
448 | Float_t etamax = header->GetLegoEtaMax(); | |
449 | Float_t etamin = header->GetLegoEtaMin(); | |
450 | Int_t ndiv = 100; | |
451 | // Get number of tracks from EventCalTrk | |
452 | Int_t nTracks = fEvent->GetNCalTrkTracks(); | |
453 | ||
454 | // jet energy and area arrays | |
455 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
456 | TH1::AddDirectory(kFALSE); | |
457 | for(Int_t mjet=0; mjet<nJ; mjet++){ | |
458 | if(!fhEtJet[mjet]){ | |
459 | fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax); | |
460 | } | |
461 | if(!fhAreaJet[mjet]){ | |
462 | fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax); | |
463 | } | |
464 | fhEtJet[mjet]->Reset(); | |
465 | fhAreaJet[mjet]->Reset(); | |
466 | } | |
467 | // background energy and area | |
468 | if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); | |
469 | fhEtBackg->Reset(); | |
470 | if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); | |
471 | fhAreaBackg->Reset(); | |
472 | TH1::AddDirectory(oldStatus); | |
473 | ||
474 | //fill energies | |
475 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
476 | for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets | |
477 | Float_t deta = etaT[jpart] - etaJet[ijet]; | |
478 | Float_t dphi = phiT[jpart] - phiJet[ijet]; | |
479 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
480 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
481 | Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
482 | if(dr <= rc){ // particles inside this cone | |
483 | if(jpart < nTracks) multJetT[ijet]++; | |
484 | else multJetC[ijet]++; | |
485 | multJet[ijet]++; | |
486 | injet[jpart] = ijet; | |
487 | ||
488 | if(PtCutPass(jpart,nTracks)){ // pt cut | |
489 | fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone | |
490 | if(SignalCutPass(jpart,nTracks)) | |
491 | etsigJet[ijet]+= ptT[jpart]; | |
492 | } | |
493 | break; | |
494 | } | |
495 | }// end jets loop | |
496 | if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones | |
497 | } //end particle loop | |
498 | ||
499 | //calc areas | |
500 | Float_t eta0 = etamin; | |
501 | Float_t etaw = (etamax - etamin)/((Float_t)ndiv); | |
502 | Float_t eta1 = eta0 + etaw; | |
503 | for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins | |
504 | Float_t etac = eta0 + etaw/2.0; | |
505 | Float_t areabg = etaw*2.0*TMath::Pi(); | |
506 | for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets | |
507 | Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]); | |
508 | Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]); | |
509 | Float_t acc0 = 0.0; Float_t acc1 = 0.0; | |
510 | Float_t areaj = 0.0; | |
511 | if(deta0 > rc && deta1 < rc){ | |
512 | acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); | |
513 | areaj = acc1; | |
514 | } | |
515 | if(deta0 < rc && deta1 > rc){ | |
516 | acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); | |
517 | areaj = acc0; | |
518 | } | |
519 | if(deta0 < rc && deta1 < rc){ | |
520 | acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); | |
521 | acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); | |
522 | if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1 | |
523 | if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2 | |
524 | if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3 | |
525 | } | |
526 | fhAreaJet[ijet]->Fill(etac,areaj); | |
527 | areabg = areabg - areaj; | |
528 | } // end jets loop | |
529 | fhAreaBackg->Fill(etac,areabg); | |
530 | eta0 = eta1; | |
531 | eta1 = eta1 + etaw; | |
532 | } // end loop for all eta bins | |
533 | ||
534 | //subtract background | |
535 | for(Int_t kjet=0; kjet<nJ; kjet++){ | |
536 | etJet[kjet] = 0.0; // first clear etJet for this jet | |
537 | for(Int_t bin = 0; bin< ndiv; bin++){ | |
538 | if(fhAreaJet[kjet]->GetBinContent(bin)){ | |
539 | Float_t areab = fhAreaBackg->GetBinContent(bin); | |
540 | Float_t etb = fhEtBackg->GetBinContent(bin); | |
541 | Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab; | |
542 | etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction | |
543 | } | |
544 | } | |
545 | } | |
50c7d9f7 | 546 | |
139cbd96 | 547 | // calc background total |
548 | Double_t etOut = fhEtBackg->Integral(); | |
549 | Double_t areaOut = fhAreaBackg->Integral(); | |
550 | Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); | |
551 | etbgTotalN = etOut*areaT/areaOut; | |
552 | ||
553 | Int_t count=0; | |
554 | ||
555 | for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array | |
556 | if((injet[jpart] == -1) && | |
557 | (PtCutPass(jpart,nTracks))){ | |
558 | sigmaN += etbgTotalN/areaT - ptT[jpart]; | |
559 | count=count+1; | |
560 | } | |
561 | } | |
562 | sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); | |
50c7d9f7 | 563 | |
50c7d9f7 | 564 | } |
139cbd96 | 565 |