]> 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
propagate changes to DEV
[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
16#include <Riostream.h>
17#include <TChain.h>
18#include <TFile.h>
19#include <TList.h>
20#include <TArrayF.h>
21#include <TClonesArray.h>
22#include <TF1.h>
23#include <TString.h>
24
25#include "AliJetHeader.h"
26#include "AliJetReader.h"
27#include "AliJetReaderHeader.h"
28#include "AliFastJetFinder.h"
29#include "AliFastJetHeaderV1.h"
30#include "AliJetReaderHeader.h"
31#include "AliJetReader.h"
32#include "AliJetUnitArray.h"
33#include "AliFastJetInput.h"
34#include "AliESDEvent.h"
35
36
37#include "fastjet/PseudoJet.hh"
38#include "fastjet/ClusterSequenceArea.hh"
39#include "fastjet/AreaDefinition.hh"
40#include "fastjet/JetDefinition.hh"
41// get info on how fastjet was configured
42#include "fastjet/config.h"
43
44#ifdef ENABLE_PLUGIN_SISCONE
45#include "fastjet/SISConePlugin.hh"
46#endif
47
48#include<sstream> // needed for internal io
49#include<vector>
50#include <cmath>
51
52using namespace std;
53#include "AliJetBkg.h"
54
55ClassImp(AliJetBkg)
56
57////////////////////////////////////////////////////////////////////////
58
59AliJetBkg::AliJetBkg():
856618e7 60 TObject(),
61 fReader(0),
62 fHeader(0),
63 fInputFJ(0)
50c7d9f7 64{
65 // Default constructor
66
67}
68//______________________________________________________________________
69AliJetBkg::AliJetBkg(const AliJetBkg& input):
856618e7 70 TObject(input),
50c7d9f7 71 fReader(input.fReader),
72 fHeader(input.fHeader),
856618e7 73 fInputFJ(input.fInputFJ)
50c7d9f7 74{
75 // copy constructor
76
77}
78//______________________________________________________________________
79AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
80 // Assignment operator.
81 this->~AliJetBkg();
82 new(this) AliJetBkg(source);
83 return *this;
84
85}
294056ab 86
87
88//___________________________________________________________________
89 void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma,
90Double_t& meanarea){
674c4c04 91
92 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
93 Bool_t debug = header->GetDebug(); // debug option
94 if(debug)cout<<"=============== AliJetBkg::BkgFastJetb() =========== "<<endl;
294056ab 95 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
96
97 //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
98
99
100
674c4c04 101
294056ab 102 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
103
104 Double_t medianb,sigmab,meanareab;
105 CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
106 rho=medianb;
107 sigma=sigmab;
108 meanarea=meanareab;
109
110}
111//_________________________________________________________________
112
113 void AliJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma,
114Double_t& meanarea){
674c4c04 115 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
116 Bool_t debug = header->GetDebug(); // debug option
117 if(debug)cout<<"=============== AliJetBkg::BkgWoHardest() =========== "<<endl;
294056ab 118 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
119
674c4c04 120 //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
294056ab 121 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
122 Double_t medianb,sigmab,meanareab;
123 CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
124 rho=medianb;
125 sigma=sigmab;
126 meanarea=meanareab;
127
128 }
129
130//____________________________________________________________________
131
132 void AliJetBkg::CalcRhob(Double_t& median,Double_t&
bffdedda 133 sigma,Double_t&
134 meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
135 rParamBkg,TString method){
294056ab 136 //calculate rho using the fastjet method
137
138 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
139 Bool_t debug = header->GetDebug(); // debug option
140
141 fastjet::Strategy strategy = header->GetStrategy();
142 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
143 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
144 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
145
146 // create an object that specifies how we to define the area
147 fastjet::AreaDefinition areaDef;
148 double ghostEtamax = header->GetGhostEtaMax();
149 double ghostArea = header->GetGhostArea();
150 int activeAreaRepeats = header->GetActiveAreaRepeats();
151
152 // now create the object that holds info about ghosts
153
154
155
156 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
157 // and from that get an area definition
158 fastjet::AreaType areaType = header->GetAreaType();
159 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
160
161 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
162 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
163 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
164 comment+= "Jet definition: ";
165 comment+= TString(jetDef.description());
166 // comment+= ". Area definition: ";
167 //comment+= TString(areaDef.description());
168 comment+= ". Strategy adopted by FastJet: ";
169 comment+= TString(clust_seq.strategy_string());
bffdedda 170 comment+= Form("Method: %s",method.Data());
294056ab 171 header->SetComment(comment);
172 if(debug){
173 cout << "--------------------------------------------------------" << endl;
174 cout << comment << endl;
175 cout << "--------------------------------------------------------" << endl;
176 }
177
178 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
179 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
180
181
182 //cout<<"# of BKG jets = "<<jets.size()<<endl;
183 // for (size_t j = 0; j < jets.size(); j++) { // loop for jets
184 //printf("BKG Jet found %5d %9.5f %8.5f %10.3f %4.4f
185 //\n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(),
186 //clust_seq.area(jets[j]));
187 // }
188
189
190 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
191
192
193 phiMin = 0;
194 phiMax = 2*TMath::Pi();
195
196 rapMax = ghostEtamax - rParamBkg;
197 rapMin = - ghostEtamax + rParamBkg;
198
199
200 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
201
202
203 double medianb, sigmab, meanareab;
204 clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false);
205 median=medianb;
206 sigma=sigmab;
207 meanarea=meanareab;
208
209
210 }
211
212//____________________________________________________________________
213
214
215 void AliJetBkg::CalcRhoWoHardest(Double_t& median,Double_t&
216sigma,Double_t&
217meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
218rParamBkg,TString method){
219 //calculate rho using the fastjet method
220
221 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
222 Bool_t debug = header->GetDebug(); // debug option
223
224 fastjet::Strategy strategy = header->GetStrategy();
225 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
226 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
227 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
228
229 // create an object that specifies how we to define the area
230 fastjet::AreaDefinition areaDef;
231 double ghostEtamax = header->GetGhostEtaMax();
232 double ghostArea = header->GetGhostArea();
233 int activeAreaRepeats = header->GetActiveAreaRepeats();
234
235 // now create the object that holds info about ghosts
236
237
238
239 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
240 // and from that get an area definition
241 fastjet::AreaType areaType = header->GetAreaType();
242 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
243 //cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
244 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
245 fastjet::ClusterSequenceArea clust_seq(inputParticles,jetDef,areaDef);
246 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
247 comment+= "Jet definition: ";
248 comment+= TString(jetDef.description());
249 // comment+= ". Area definition: ";
250 //comment+= TString(areaDef.description());
251 comment+= ". Strategy adopted by FastJet: ";
252 comment+= TString(clust_seq.strategy_string());
bffdedda 253 comment+= Form("Method: %s",method.Data());
294056ab 254 header->SetComment(comment);
255 if(debug){
256 cout << "--------------------------------------------------------" << endl;
257 cout << comment << endl;
258 cout << "--------------------------------------------------------" << endl;
259 }
260
261 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
262 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
263 vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets);
264 if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1);
265
266
267
268 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
269
270
271 phiMin = 0;
272 phiMax = 2*TMath::Pi();
273
274 rapMax = ghostEtamax - rParamBkg;
275 rapMin = - ghostEtamax + rParamBkg;
276
277
278 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
279
280
281 double medianb, sigmab, meanareab;
282 clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab,
283meanareab, false);
284 median=medianb;
285 sigma=sigmab;
286 meanarea=meanareab;
287
288
289 }
290
291
50c7d9f7 292//___________________________________________________________________
293Float_t AliJetBkg::BkgFastJet(){
1240edf5 294// Return background
967c266c 295 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
296 Bool_t debug = header->GetDebug(); // debug option
297
298 if(debug)cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
50c7d9f7 299 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
300
967c266c 301 if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
50c7d9f7 302
303 for(UInt_t i=0;i<inputParticles.size();i++){
304 // cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
305
306 }
307
50c7d9f7 308 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
309 Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
967c266c 310 if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
50c7d9f7 311 return rho;
312
313}
314//___________________________________________________________________
315Float_t AliJetBkg::BkgChargedFastJet(){
1240edf5 316// Background for charged jets
967c266c 317 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
318 Bool_t debug = header->GetDebug(); // debug option
319
320 if(debug)cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
50c7d9f7 321
322 vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
323
967c266c 324 if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
50c7d9f7 325
326 for(UInt_t i=0;i<inputParticlesCharged.size();i++){
327 // cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
328
329 }
967c266c 330
50c7d9f7 331 double rParam = header->GetRparam();
332
333 Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
334
674c4c04 335 if(debug)cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
50c7d9f7 336 return rho;
337}
338
339
340
341Float_t AliJetBkg::BkgStat()
342{
343 //background subtraction using statistical method
344
967c266c 345 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
346 Bool_t debug = header->GetDebug(); // debug option
347
348 if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
50c7d9f7 349 //TO BE IMPLEMENTED
350 // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
351 Int_t nTracks= 0;
352 TF1 fun("fun",BkgFunction,0,800,1);
353 Double_t enTot=fun.Eval(nTracks);
354 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
355 return enTot/accEMCal;
356
357}
50c7d9f7 358
359////////////////////////////////////////////////////////////////////////
360Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
361{
362
363 // Cone background subtraction method applied on the fastjet: REmove the particles of the
364 // two largest jets with the given R from the estimation of new rho.
50c7d9f7 365
366 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
967c266c 367 Bool_t debug = header->GetDebug(); // debug option
368
369 if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
370
50c7d9f7 371 Float_t rc= header->GetRparam();
372
373 //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
374 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
375
376 Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
967c266c 377 if(debug)cout<<"nJets: "<<nJ<<endl;
50c7d9f7 378
379
380 //begin unit array
381 TClonesArray* fUnit = fReader->GetUnitArray();
382 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
383
384 Int_t nIn = fUnit->GetEntries();
385 if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
386
387 // Information extracted from fUnitArray
388 // load input vectors and calculate total energy in array
389 Float_t pt,eta,phi;
390 Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
391 Float_t rhoback=0.0;
392
393 Float_t ptallback=0.0; //particles without the jet
394 Float_t restarea=accEMCal; //initial area set
395 Bool_t acc=0;
396 Bool_t acc1=0;
397 Float_t rCone=0.4;
398
399 if(nJ==1) {
225f0094 400 AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
50c7d9f7 401 jeteta=jettmp->Eta();
402 jetphi=jettmp->Phi();
403 acc=EmcalAcceptance(jeteta,jetphi,rCone);
404 if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
967c266c 405 if(debug)cout<<" acc "<<acc<<endl;
50c7d9f7 406 }
407
408
409 if(nJ>=2) {
225f0094 410 AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
411 AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1));
50c7d9f7 412 jeteta=jettmp->Eta();
413 jetphi=jettmp->Phi();
414 jeteta1=jettmp1->Eta();
415 jetphi1=jettmp1->Phi();
416 acc=EmcalAcceptance(jeteta,jetphi,rCone);
417 acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
418 if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
419 if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
420 if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
421
967c266c 422 if(debug)cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
50c7d9f7 423
424 }
425
426 // cout<<" nIn = "<<nIn<<endl;
427 Float_t sumpt=0;
428 for(Int_t i=0; i<nIn; i++)
429 { //Unit Array Loop
430 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
431 if(uArray->GetUnitEnergy()>0.){
432
433 pt = uArray->GetUnitEnergy();
434 eta = uArray->GetUnitEta();
435 phi = uArray->GetUnitPhi();
436
437 //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
438
439 Float_t deta=0.0, dphi=0.0, dr=100.0;
440 Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
441
442 //cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
443 if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
444 sumpt+=pt;
445 //if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
446
447 if(nJ==1 && acc==1) {
448 deta = eta - jeteta;
449 dphi = phi - jetphi;
450 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
451 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
452 dr = TMath::Sqrt(deta * deta + dphi * dphi);
453 if(dr<=rc)sumpt-=pt;
454 }
455
456 if(nJ>=2) {
457 if(acc==1){
458 deta = eta - jeteta;
459 dphi = phi - jetphi;
460 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
461 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
462 dr = TMath::Sqrt(deta * deta + dphi * dphi);
463 if(dr<=rc)sumpt-=pt;
464 }
465 if(acc1==1){
466 deta1 = eta - jeteta1;
467 dphi1 = phi - jetphi1;
468 if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
469 if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
470 dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
471 if(dr1<=rc)sumpt-=pt;
472 }
473
474 }
475
476 if(dr >= rc && dr1 >=rc) {
477 // particles outside both cones
478
479 //cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl;
480 //cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl;
481 ptallback+=pt;
482 }
483 }
484 //cout<<" ipart "<<ipart<<" rhointegral "<<rhoback <<endl;
485 }
486 } // End loop on UnitArray
487
967c266c 488 if(debug)cout<<"total area left "<<restarea<<endl;
489 if(debug)cout<<"sumpt="<<sumpt<<endl;
50c7d9f7 490 // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
491 //else rhoback=ptallback;
492
493 rhoback= ptallback/restarea;
967c266c 494 if(debug)cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
50c7d9f7 495
496 return rhoback;
497
498}
499
500
501Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
502 //calculate rho using the fastjet method
503
504 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
505 Bool_t debug = header->GetDebug(); // debug option
506
507 fastjet::Strategy strategy = header->GetStrategy();
508 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
509 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
510 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
511
512 // create an object that specifies how we to define the area
513 fastjet::AreaDefinition areaDef;
514 double ghostEtamax = header->GetGhostEtaMax();
515 double ghostArea = header->GetGhostArea();
516 int activeAreaRepeats = header->GetActiveAreaRepeats();
517
518 // now create the object that holds info about ghosts
519
520 if (method.Contains("Charg"))ghostEtamax=0.9;
521
522 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
523 // and from that get an area definition
524 fastjet::AreaType areaType = header->GetAreaType();
525 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
967c266c 526 if(debug)cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
50c7d9f7 527 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
528 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
529 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
530 comment+= "Jet definition: ";
531 comment+= TString(jetDef.description());
532 // comment+= ". Area definition: ";
533 //comment+= TString(areaDef.description());
534 comment+= ". Strategy adopted by FastJet: ";
535 comment+= TString(clust_seq.strategy_string());
536 header->SetComment(comment);
537 if(debug){
538 cout << "--------------------------------------------------------" << endl;
539 cout << comment << endl;
540 cout << "--------------------------------------------------------" << endl;
541 }
542
543 double ptmin = header->GetPtMin();
544 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
545 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
50c7d9f7 546
db2272d5 547 if (debug) {
548 cout<<"# of BKG jets = "<<jets.size()<<endl;
549 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
550 printf("BKG Jet found %5d %9.5f %8.5f %10.3f %4.4f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(),clust_seq.area(jets[j]));
551 }
552 }
553
50c7d9f7 554 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
555
556 if (method.Contains("All")){
557 phiMin = 80.*TMath::Pi()/180+rParamBkg;
558 phiMax = 190.*TMath::Pi()/180-rParamBkg;
50c7d9f7 559 }
560 if (method.Contains("Charg")){
561 phiMin = 0;
562 phiMax = 2*TMath::Pi();
563 }
564 rapMax = ghostEtamax - rParamBkg;
565 rapMin = - ghostEtamax + rParamBkg;
566
567
568 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
569
570
571 Double_t rho=clust_seq.median_pt_per_unit_area(range);
572 // double median, sigma, meanArea;
573 //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
574 //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
575
576 // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
577
578
967c266c 579 if(debug)cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
50c7d9f7 580 return rho;
581}
582
583Float_t AliJetBkg::EtaToTheta(Float_t arg)
584{
585 // return (180./TMath::Pi())*2.*atan(exp(-arg));
586 return 2.*atan(exp(-arg));
587}
588
589
590Double_t AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
591 //to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
592 return 1;
593}
594
595
1240edf5 596Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius) const{
597// Apply emcal acceptance cuts
50c7d9f7 598 Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
599 Float_t deltaphi=110./180.*TMath::Pi();
600 Float_t phicut=deltaphi/2.-radius;
601 Float_t etacut=0.7-radius;
602 //cout<<" eta "<<eta<<" phi "<<phi<<endl;
603 //cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl;
604 if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
605 else return 0;
606}