1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //--------------------------------------------------
19 // Method implementation for background studies and background subtraction with UA1 algorithms
21 // Author: magali.estienne@subatech.in2p3.fr
22 //-------------------------------------------------
24 #include <Riostream.h>
26 #include <TClonesArray.h>
30 #include "AliAODJetEventBackground.h"
31 #include "AliUA1JetHeaderV1.h"
32 #include "AliJetCalTrk.h"
33 #include "AliJetBkg.h"
40 ////////////////////////////////////////////////////////////////////////
42 AliJetBkg::AliJetBkg():
50 // Default constructor
51 for(int i = 0;i < kMaxJets;i++){
52 fhAreaJet[i] = fhEtJet[i] = 0;
56 //----------------------------------------------------------------
57 AliJetBkg::AliJetBkg(const AliJetBkg& input):
60 fHeader(input.fHeader),
62 fhEtBackg(input.fhEtBackg),
63 fhAreaBackg(input.fhAreaBackg)
66 for(int i = 0;i < kMaxJets;i++){
67 fhAreaJet[i] = input.fhAreaJet[i];
68 fhEtJet[i] = input.fhEtJet[i];
73 //----------------------------------------------------------------
74 AliJetBkg::~AliJetBkg()
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];
86 //----------------------------------------------------------------
87 Bool_t AliJetBkg::PtCutPass(Int_t id, Int_t nTracks)
89 // Check if track or cell passes the cut flag
90 if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetCutFlag() == 1)
96 //----------------------------------------------------------------
97 Bool_t AliJetBkg::SignalCutPass(Int_t id, Int_t nTracks)
99 // Check if track or cell passes the cut flag
100 if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetSignalFlag() == 1)
106 //----------------------------------------------------------------
107 Float_t AliJetBkg::CalcJetAreaEtaCut(Float_t radius, const Float_t etaJet)
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);
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);
123 return radius*radius*TMath::Pi() - accmax - accmin;
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)
130 // Calculate jet and bacground areas taking into account an acceptance cut in eta
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];
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)
149 // Background subtraction using cone method but without correction in dE/deta distribution
150 // Cases to take into account the EMCal geometry are included
153 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
154 //calculate energy inside and outside cones
155 fDebug = header->GetDebug();
156 Float_t rc = header->GetRadius();
158 // Get number of tracks from EventCalTrk
159 Int_t nTracks = fEvent->GetNCalTrkTracks();
161 Float_t etIn[kMaxJets] = {0};
162 Float_t areaOut = 0.;
164 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
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]++;
177 if(PtCutPass(jpart,nTracks)){ // pt cut
178 etIn[ijet] += ptT[jpart];
179 if(SignalCutPass(jpart,nTracks))
180 etsigJet[ijet]+= ptT[jpart];
186 if((injet[jpart] == -1) &&
187 (PtCutPass(jpart,nTracks))){
188 etOut += ptT[jpart]; // particle outside cones and pt cut
190 } //end particle loop
192 // Calculate jet and background areas
193 Bool_t calcAreaOut = kTRUE;
194 CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
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
202 // estimate new total background
204 areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
205 etbgTotalN = etOut*areaT/areaOut;
207 // estimate standard deviation of background
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 ?)
219 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
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)
231 //background subtraction using statistical method
232 // Cases to take into account the EMCal geometry are included
235 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
236 Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
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.;
245 for(Int_t jpart = 0; jpart < nIn; jpart++)
246 { // loop for all particles in array
248 for(Int_t ijet=0; ijet<nJ; ijet++)
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]++;
261 if(PtCutPass(jpart,nTracks)){ // pt cut
262 etIn[ijet] += ptT[jpart];
263 if(SignalCutPass(jpart,nTracks))
264 etsigJet[ijet]+= ptT[jpart];
269 } //end particle loop
271 // Calculate jet and background areas
272 Bool_t calcAreaOut = kFALSE;
273 CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
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
281 etbgTotalN = etbgStat;
283 // estimate standard deviation of background
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];
293 if(count>0)sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
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*/)
304 // Cone background subtraction method taking into acount dEt/deta distribution
305 // Cases to take into account the EMCal geometry are not included
308 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
310 Float_t rc= header->GetRadius();
311 Float_t etamax = header->GetLegoEtaMax();
312 Float_t etamin = header->GetLegoEtaMin();
314 // Get number of tracks from EventCalTrk
315 Int_t nTracks = fEvent->GetNCalTrkTracks();
317 // jet energy and area arrays
318 Bool_t oldStatus = TH1::AddDirectoryStatus();
319 TH1::AddDirectory(kFALSE);
321 for(Int_t mjet=0; mjet<nJ; mjet++){
323 fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
325 if(!fhAreaJet[mjet]){
326 fhAreaJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
328 fhEtJet[mjet]->Reset();
329 fhAreaJet[mjet]->Reset();
331 // background energy and area
332 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
334 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
335 fhAreaBackg->Reset();
336 TH1::AddDirectory(oldStatus);
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]++;
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];
361 if((injet[jpart] == -1) &&
362 (PtCutPass(jpart,nTracks) == 1))
363 fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
364 } //end particle loop
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;
378 if(deta0 > rc && deta1 < rc){
379 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
382 if(deta0 < rc && deta1 > rc){
383 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
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
393 fhAreaJet[ijet]->Fill(etac,areaj);
394 areabg = areabg - areaj;
396 fhAreaBackg->Fill(etac,areabg);
399 } // end loop for all eta bins
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
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;
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];
428 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
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*/)
439 // Ratio background subtraction method taking into acount dEt/deta distribution
440 // Cases to take into account the EMCal geometry are not included
442 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
443 //factor F calc before
444 Float_t bgRatioCut = header->GetBackgCutRatio();
447 Float_t rc= header->GetRadius();
448 Float_t etamax = header->GetLegoEtaMax();
449 Float_t etamin = header->GetLegoEtaMin();
451 // Get number of tracks from EventCalTrk
452 Int_t nTracks = fEvent->GetNCalTrkTracks();
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++){
459 fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
461 if(!fhAreaJet[mjet]){
462 fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
464 fhEtJet[mjet]->Reset();
465 fhAreaJet[mjet]->Reset();
467 // background energy and area
468 if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
470 if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
471 fhAreaBackg->Reset();
472 TH1::AddDirectory(oldStatus);
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]++;
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];
496 if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
497 } //end particle loop
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;
511 if(deta0 > rc && deta1 < rc){
512 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
515 if(deta0 < rc && deta1 > rc){
516 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
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
526 fhAreaJet[ijet]->Fill(etac,areaj);
527 areabg = areabg - areaj;
529 fhAreaBackg->Fill(etac,areabg);
532 } // end loop for all eta bins
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
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;
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];
562 sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);