]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliJetBkg.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug in storage manager making possible to delete permanent event fixed
[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.cxx
CommitLineData
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
36using namespace std;
37
50c7d9f7 38ClassImp(AliJetBkg)
39
40////////////////////////////////////////////////////////////////////////
41
42AliJetBkg::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 57AliJetBkg::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//----------------------------------------------------------------
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 }
294056ab 83
139cbd96 84}
674c4c04 85
139cbd96 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;
294056ab 93
139cbd96 94}
674c4c04 95
139cbd96 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;
294056ab 103
294056ab 104}
294056ab 105
139cbd96 106//----------------------------------------------------------------
28ab787b 107Float_t AliJetBkg::CalcJetAreaEtaCut(Float_t radius, Float_t etaJet)
139cbd96 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//----------------------------------------------------------------
28ab787b 128void AliJetBkg::CalcJetAndBckgAreaEtaCut(Bool_t calcOutsideArea, Float_t radius, Int_t nJets, const Float_t* etaJet, Float_t* &areaJet, Float_t &areaOut)
139cbd96 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//----------------------------------------------------------------
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 //
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//----------------------------------------------------------------
224void 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//----------------------------------------------------------------
298void 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//----------------------------------------------------------------
433void 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