JETAN/AliFastJetBkg.cxx

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