]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/DEV/AliJetBkg.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding JETAN and FASTJETAN development libs for new i/o of tracks/particles for the...
[u/mrichter/AliRoot.git] / JETAN / DEV / AliJetBkg.cxx
CommitLineData
d89b8229 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//--------------------------------------------------
19// Method implementation for background studies and background subtraction with UA1 algorithms
20//
21// Author: magali.estienne@subatech.in2p3.fr
22//-------------------------------------------------
23
24#include <Riostream.h>
25#include <TList.h>
26#include <TClonesArray.h>
27#include <TH1F.h>
28#include <TH2F.h>
29
30#include "AliAODJetEventBackground.h"
31#include "AliUA1JetHeaderV1.h"
32#include "AliJetCalTrk.h"
33#include "AliJetBkg.h"
34
35
36using namespace std;
37
38ClassImp(AliJetBkg)
39
40////////////////////////////////////////////////////////////////////////
41
42AliJetBkg::AliJetBkg():
43 TObject(),
44 fEvent(0x0),
45 fHeader(0x0),
46 fDebug(0),
47 fhEtBackg(0x0),
48 fhAreaBackg(0x0)
49{
50 // Default constructor
51 for(int i = 0;i < kMaxJets;i++){
52 fhAreaJet[i] = fhEtJet[i] = 0;
53 }
54}
55
56//----------------------------------------------------------------
57AliJetBkg::AliJetBkg(const AliJetBkg& input):
58 TObject(input),
59 fEvent(input.fEvent),
60 fHeader(input.fHeader),
61 fDebug(input.fDebug),
62 fhEtBackg(input.fhEtBackg),
63 fhAreaBackg(input.fhAreaBackg)
64{
65 // copy constructor
66 for(int i = 0;i < kMaxJets;i++){
67 fhAreaJet[i] = input.fhAreaJet[i];
68 fhEtJet[i] = input.fhEtJet[i];
69 }
70
71}
72
73//----------------------------------------------------------------
74AliJetBkg::~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 }
83
84}
85
86//----------------------------------------------------------------
87Bool_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;
93
94}
95
96//----------------------------------------------------------------
97Bool_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;
103
104}
105
106//----------------------------------------------------------------
107Float_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);
117 }
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);
121 }
122
123 return radius*radius*TMath::Pi() - accmax - accmin;
124
125}
126
127//----------------------------------------------------------------
128void 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
131
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];
137 }
138
139}
140
141//----------------------------------------------------------------
142void 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 //
152
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 }
201
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 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
219
220}
221
222//----------------------------------------------------------------
223void AliJetBkg::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN, Float_t&sigmaN,
224 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
225 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
226 Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
227 Int_t* injet, Float_t* &areaJet)
228{
229 //
230 //background subtraction using statistical method
231 // Cases to take into account the EMCal geometry are included
232 //
233
234 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
235 Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
236
237 //calculate energy inside
238 Float_t rc= header->GetRadius();
239 Float_t etIn[kMaxJets] = {0.0};
240 // Get number of tracks from EventCalTrk
241 Int_t nTracks = fEvent->GetNCalTrkTracks();
242 Float_t areaOut = 0.;
243
244 for(Int_t jpart = 0; jpart < nIn; jpart++)
245 { // loop for all particles in array
246
247 for(Int_t ijet=0; ijet<nJ; ijet++)
248 {
249 Float_t deta = etaT[jpart] - etaJet[ijet];
250 Float_t dphi = phiT[jpart] - phiJet[ijet];
251 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
252 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
253 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
254 if(dr <= rc){ // particles inside this cone
255 if(jpart < nTracks) multJetT[ijet]++;
256 else multJetC[ijet]++;
257 multJet[ijet]++;
258 injet[jpart] = ijet;
259
260 if(PtCutPass(jpart,nTracks)){ // pt cut
261 etIn[ijet] += ptT[jpart];
262 if(SignalCutPass(jpart,nTracks))
263 etsigJet[ijet]+= ptT[jpart];
264 }
265 break;
266 }
267 }// end jets loop
268 } //end particle loop
269
270 // Calculate jet and background areas
271 Bool_t calcAreaOut = kFALSE;
272 CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
273
274 //subtract background using area method
275 for(Int_t ljet=0; ljet<nJ; ljet++){
276 Float_t areaRatio = areaJet[ljet]/areaOut;
277 etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
278 }
279 Int_t count=0;
280 etbgTotalN = etbgStat;
281
282 // estimate standard deviation of background
283 Float_t areaT = 0;
284 areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
285 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
286 if((injet[jpart] == -1) &&
287 (PtCutPass(jpart,nTracks))){
288 sigmaN += etbgTotalN/areaT - ptT[jpart];
289 count=count+1;
290 }
291 }
292 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
293
294}
295
296//----------------------------------------------------------------
297void AliJetBkg::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
298 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, Float_t* etJet,
299 const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
300 Int_t* multJetT, Int_t* multJetC, Int_t* multJet, Int_t* injet, Float_t* &/*areaJet*/)
301{
302 //
303 // Cone background subtraction method taking into acount dEt/deta distribution
304 // Cases to take into account the EMCal geometry are not included
305 //
306
307 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
308 //general
309 Float_t rc= header->GetRadius();
310 Float_t etamax = header->GetLegoEtaMax();
311 Float_t etamin = header->GetLegoEtaMin();
312 Int_t ndiv = 100;
313 // Get number of tracks from EventCalTrk
314 Int_t nTracks = fEvent->GetNCalTrkTracks();
315
316 // jet energy and area arrays
317 Bool_t oldStatus = TH1::AddDirectoryStatus();
318 TH1::AddDirectory(kFALSE);
319
320 for(Int_t mjet=0; mjet<nJ; mjet++){
321 if(!fhEtJet[mjet]){
322 fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
323 }
324 if(!fhAreaJet[mjet]){
325 fhAreaJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
326 }
327 fhEtJet[mjet]->Reset();
328 fhAreaJet[mjet]->Reset();
329 }
330 // background energy and area
331 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
332 fhEtBackg->Reset();
333 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
334 fhAreaBackg->Reset();
335 TH1::AddDirectory(oldStatus);
336
337 //fill energies
338 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
339 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
340 Float_t deta = etaT[jpart] - etaJet[ijet];
341 Float_t dphi = phiT[jpart] - phiJet[ijet];
342 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
343 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
344 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
345 if(dr <= rc){ // particles inside this cone
346 if(jpart < nTracks) multJetT[ijet]++;
347 else multJetC[ijet]++;
348 multJet[ijet]++;
349 injet[jpart] = ijet;
350
351 if(PtCutPass(jpart,nTracks)){ // pt cut
352 fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
353 if(SignalCutPass(jpart,nTracks))
354 etsigJet[ijet]+= ptT[jpart];
355 }
356 break;
357 }
358 }// end jets loop
359
360 if((injet[jpart] == -1) &&
361 (PtCutPass(jpart,nTracks) == 1))
362 fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
363 } //end particle loop
364
365 //calc areas
366 Float_t eta0 = etamin;
367 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
368 Float_t eta1 = eta0 + etaw;
369 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
370 Float_t etac = eta0 + etaw/2.0;
371 Float_t areabg = etaw*2.0*TMath::Pi();
372 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
373 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
374 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
375 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
376 Float_t areaj = 0.0;
377 if(deta0 > rc && deta1 < rc){
378 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
379 areaj = acc1;
380 }
381 if(deta0 < rc && deta1 > rc){
382 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
383 areaj = acc0;
384 }
385 if(deta0 < rc && deta1 < rc){
386 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
387 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
388 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
389 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
390 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
391 }
392 fhAreaJet[ijet]->Fill(etac,areaj);
393 areabg = areabg - areaj;
394 } // end jets loop
395 fhAreaBackg->Fill(etac,areabg);
396 eta0 = eta1;
397 eta1 = eta1 + etaw;
398 } // end loop for all eta bins
399
400 //subtract background
401 for(Int_t kjet=0; kjet<nJ; kjet++){
402 etJet[kjet] = 0.0; // first clear etJet for this jet
403 for(Int_t bin = 0; bin< ndiv; bin++){
404 if(fhAreaJet[kjet]->GetBinContent(bin)){
405 Float_t areab = fhAreaBackg->GetBinContent(bin);
406 Float_t etb = fhEtBackg->GetBinContent(bin);
407 Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
408 etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
409 }
410 }
411 }
412
413 // calc background total
414 Double_t etOut = fhEtBackg->Integral();
415 Double_t areaOut = fhAreaBackg->Integral();
416 Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
417 etbgTotalN = etOut*areaT/areaOut;
418
419 Int_t count=0;
420 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
421 if((injet[jpart] == -1) &&
422 (PtCutPass(jpart,nTracks))){
423 sigmaN += etbgTotalN/areaT - ptT[jpart];
424 count=count+1;
425 }
426 }
427 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
428
429}
430
431//----------------------------------------------------------------
432void AliJetBkg::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
433 const Float_t* ptT,const Float_t* etaT, const Float_t* phiT,
434 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
435 Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
436 Int_t* injet, Float_t* &/*areaJet*/)
437{
438 // Ratio background subtraction method taking into acount dEt/deta distribution
439 // Cases to take into account the EMCal geometry are not included
440
441 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
442 //factor F calc before
443 Float_t bgRatioCut = header->GetBackgCutRatio();
444
445 //general
446 Float_t rc= header->GetRadius();
447 Float_t etamax = header->GetLegoEtaMax();
448 Float_t etamin = header->GetLegoEtaMin();
449 Int_t ndiv = 100;
450 // Get number of tracks from EventCalTrk
451 Int_t nTracks = fEvent->GetNCalTrkTracks();
452
453 // jet energy and area arrays
454 Bool_t oldStatus = TH1::AddDirectoryStatus();
455 TH1::AddDirectory(kFALSE);
456 for(Int_t mjet=0; mjet<nJ; mjet++){
457 if(!fhEtJet[mjet]){
458 fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
459 }
460 if(!fhAreaJet[mjet]){
461 fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
462 }
463 fhEtJet[mjet]->Reset();
464 fhAreaJet[mjet]->Reset();
465 }
466 // background energy and area
467 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
468 fhEtBackg->Reset();
469 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
470 fhAreaBackg->Reset();
471 TH1::AddDirectory(oldStatus);
472
473 //fill energies
474 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
475 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
476 Float_t deta = etaT[jpart] - etaJet[ijet];
477 Float_t dphi = phiT[jpart] - phiJet[ijet];
478 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
479 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
480 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
481 if(dr <= rc){ // particles inside this cone
482 if(jpart < nTracks) multJetT[ijet]++;
483 else multJetC[ijet]++;
484 multJet[ijet]++;
485 injet[jpart] = ijet;
486
487 if(PtCutPass(jpart,nTracks)){ // pt cut
488 fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
489 if(SignalCutPass(jpart,nTracks))
490 etsigJet[ijet]+= ptT[jpart];
491 }
492 break;
493 }
494 }// end jets loop
495 if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
496 } //end particle loop
497
498 //calc areas
499 Float_t eta0 = etamin;
500 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
501 Float_t eta1 = eta0 + etaw;
502 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
503 Float_t etac = eta0 + etaw/2.0;
504 Float_t areabg = etaw*2.0*TMath::Pi();
505 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
506 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
507 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
508 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
509 Float_t areaj = 0.0;
510 if(deta0 > rc && deta1 < rc){
511 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
512 areaj = acc1;
513 }
514 if(deta0 < rc && deta1 > rc){
515 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
516 areaj = acc0;
517 }
518 if(deta0 < rc && deta1 < rc){
519 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
520 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
521 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
522 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
523 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
524 }
525 fhAreaJet[ijet]->Fill(etac,areaj);
526 areabg = areabg - areaj;
527 } // end jets loop
528 fhAreaBackg->Fill(etac,areabg);
529 eta0 = eta1;
530 eta1 = eta1 + etaw;
531 } // end loop for all eta bins
532
533 //subtract background
534 for(Int_t kjet=0; kjet<nJ; kjet++){
535 etJet[kjet] = 0.0; // first clear etJet for this jet
536 for(Int_t bin = 0; bin< ndiv; bin++){
537 if(fhAreaJet[kjet]->GetBinContent(bin)){
538 Float_t areab = fhAreaBackg->GetBinContent(bin);
539 Float_t etb = fhEtBackg->GetBinContent(bin);
540 Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
541 etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
542 }
543 }
544 }
545
546 // calc background total
547 Double_t etOut = fhEtBackg->Integral();
548 Double_t areaOut = fhAreaBackg->Integral();
549 Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
550 etbgTotalN = etOut*areaT/areaOut;
551
552 Int_t count=0;
553
554 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
555 if((injet[jpart] == -1) &&
556 (PtCutPass(jpart,nTracks))){
557 sigmaN += etbgTotalN/areaT - ptT[jpart];
558 count=count+1;
559 }
560 }
561 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
562
563}
564