JETAN/AliJetBkg.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 #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
  52 using namespace std;
  53 #include "AliJetBkg.h"
  54
  55 ClassImp(AliJetBkg)
  56
  57 ////////////////////////////////////////////////////////////////////////
  58
  59 AliJetBkg::AliJetBkg():
  60   TObject(),
  61   fReader(0),
  62   fHeader(0),
  63   fInputFJ(0)
  64 {
  65   // Default constructor
  66
  67 }
  68 //______________________________________________________________________
  69 AliJetBkg::AliJetBkg(const AliJetBkg& input):
  70   TObject(input),
  71     fReader(input.fReader),
  72     fHeader(input.fHeader),
  73     fInputFJ(input.fInputFJ)
  74 {
  75   // copy constructor
  76
  77 }
  78 //______________________________________________________________________
  79 AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
  80     // Assignment operator.
  81     this->~AliJetBkg();
  82     new(this) AliJetBkg(source);
  83     return *this;
  84
  85 }
  86 //___________________________________________________________________
  87 Float_t AliJetBkg::BkgFastJet(){
  88
  89   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  90   Bool_t debug  = header->GetDebug();     // debug option
  91
  92   if(debug)cout<<"===============  AliJetBkg::BkgFastJet()  =========== "<<endl;
  93   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  94
  95   if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
  96
  97   for(UInt_t i=0;i<inputParticles.size();i++){
  98     //  cout<<"                "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
  99
 100   }
 101
 102   double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
 103   Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
 104   if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
 105   return rho;
 106
 107 }
 108 //___________________________________________________________________
 109 Float_t  AliJetBkg::BkgChargedFastJet(){
 110
 111   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 112   Bool_t debug  = header->GetDebug();     // debug option
 113
 114   if(debug)cout<<"===============  AliJetBkg::BkgChargedFastJet()  =========== "<<endl;
 115
 116   vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
 117
 118   if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
 119
 120   for(UInt_t i=0;i<inputParticlesCharged.size();i++){
 121     // cout<<"                "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
 122
 123   }
 124
 125   double rParam = header->GetRparam();
 126
 127   Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
 128
 129   cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
 130   return rho;
 131 }
 132
 133
 134
 135 Float_t AliJetBkg::BkgStat()
 136 {
 137   //background subtraction using statistical method
 138
 139   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 140   Bool_t debug  = header->GetDebug();     // debug option
 141
 142   if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
 143   //TO BE IMPLEMENTED
 144   // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
 145   Int_t nTracks= 0;
 146   TF1 fun("fun",BkgFunction,0,800,1);
 147   Double_t enTot=fun.Eval(nTracks);
 148   Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
 149   return enTot/accEMCal;
 150
 151 }
 152 /////////////////////////////////
 153 Float_t AliJetBkg::BkgRemoveJetLeading(TClonesArray* fAODJets)
 154 {
 155   // Remove the particles of the
 156   // two largest jets using the track references  stored in the AODJet from the estimation of new rho.
 157
 158   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 159   Bool_t debug  = header->GetDebug();     // debug option
 160
 161   if(debug)cout<<"==============AliJetBkg::BkgRemoveJetLeading()============="<<endl;
 162
 163
 164   // check if we are reading AOD jets
 165   TRefArray *refs = 0;
 166   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetESDReader");
 167   if (fromAod) { refs = fReader->GetReferences(); }
 168
 169   //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
 170   Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
 171
 172   Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
 173   if(debug)cout<<"nJets:  "<<nJ<<endl;
 174
 175
 176   //begin unit array
 177   TClonesArray* fUnit = fReader->GetUnitArray();
 178   if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
 179
 180   Int_t nIn = fUnit->GetEntries();
 181   if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
 182
 183   Float_t rhoback=0.0;
 184   Float_t jetarea1=0.0,jetarea2=0.0;
 185
 186   Int_t particlejet1=-99;
 187   Int_t particlejet2=-99;
 188   TRefArray *refarray1 = 0;
 189   TRefArray *refarray2 = 0;
 190   Int_t nJettracks1 = 0, nJettracks2 = 0;
 191   Int_t acc=0,acc1=0;
 192   AliAODJet *jet1;
 193   AliAODJet *jet2;
 194
 195   if(nJ==1){
 196     jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
 197     jetarea1=jet1->EffectiveAreaCharged();
 198     Float_t jetPhi=jet1->Phi();
 199     Float_t jetEta=jet1->Eta();
 200     if(jetPhi>1.396 && jetPhi<3.316 && jetEta>-0.7 && jetEta<0.7)acc=1;
 201     refarray1=jet1->GetRefTracks();
 202     nJettracks1=refarray1->GetEntries();
 203     if(debug)cout<<"nJ = 1, acc="<<acc<<"  jetarea1="<<jetarea1<<endl;
 204
 205   }
 206
 207   if(nJ>=2){
 208     jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
 209     jetarea1=jet1->EffectiveAreaCharged();
 210     Float_t jetPhi1=jet1->Phi();
 211     Float_t jetEta1=jet1->Eta();
 212     if(jetPhi1>1.396 && jetPhi1<3.316 && jetEta1>-0.7 && jetEta1<0.7)acc=1;
 213     refarray1=jet1->GetRefTracks();
 214     nJettracks1=refarray1->GetEntries();
 215     if(debug)cout<<"npart = "<<nJettracks1<<endl;
 216
 217     jet2 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
 218     jetarea2=jet2->EffectiveAreaCharged();
 219     Float_t jetPhi2=jet2->Phi();
 220     Float_t jetEta2=jet2->Eta();
 221     if(jetPhi2>1.396 && jetPhi2<3.316 && jetEta2>-0.7 && jetEta2<0.7)acc1=1;
 222     refarray2=jet2->GetRefTracks();
 223     nJettracks2=refarray2->GetEntries();
 224     if(debug)cout<<"nJ = "<<nJ<<", acc="<<acc<<"  acc1="<<acc1<<"  jetarea1="<<jetarea1<<"  jetarea2="<<jetarea2<<endl;
 225   }
 226
 227
 228
 229   // cout<<" nIn = "<<nIn<<endl;
 230   Float_t sumPt=0;
 231   Float_t eta,phi,pt;
 232   Int_t ipart=0;
 233
 234   for(Int_t i=0; i<nIn; i++)
 235     { //Unit Array Loop
 236       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 237
 238       if(uArray->GetUnitEnergy()>0.){
 239         eta   = uArray->GetUnitEta();
 240         phi   = uArray->GetUnitPhi();
 241         pt = uArray->GetUnitEnergy();
 242         //      cout<<"ipart = "<<ipart<<" eta="<<eta<<"  phi="<<phi<<endl;
 243         if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
 244           //cout<<sumPt<<endl;
 245             sumPt+=pt;
 246
 247
 248             if(nJ==1 && acc==1){
 249               for(Int_t ii=0; ii<nJettracks1;ii++){
 250
 251                 particlejet1  = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
 252
 253                 if(ipart==particlejet1) {
 254                   sumPt-=pt;
 255                 }
 256               }
 257             }
 258
 259
 260             if(nJ>=2){
 261
 262               //first jet
 263               if(acc==1){
 264                 for(Int_t ii=0; ii<nJettracks1;ii++){
 265                   particlejet1  = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
 266
 267                 //cout<<"uArr loop = "<<i<<"  ipart in uArr (1/2)="<<ipart<<"  part in jet="<<ii<<"  partID="<<particlejet1<<" sumPt="<<sumPt<<endl;
 268                   if(ipart==particlejet1) {
 269                     sumPt-=pt;
 270                   }
 271                 }
 272               }
 273               if(acc1==1){
 274                 //second jet
 275                 for(Int_t ii=0; ii<nJettracks2;ii++){
 276                   particlejet2  = ((AliJetUnitArray*)refarray2->At(ii))->GetUnitTrackID();
 277                   //cout<<"uArr loop = "<<i<<"  ipart in uArr (2/2)="<<ipart<<"  part in jet="<<ii<<"  partID="<<particlejet2<<" sumPt="<<sumPt<<endl;
 278                   if(ipart==particlejet2) {
 279                     sumPt-=pt;
 280                   }
 281                 }
 282               }
 283
 284
 285             }
 286
 287         }//if phi,eta
 288         ipart++;
 289       }//end if energy
 290     }// end unit array loop
 291
 292
 293   Float_t areasum = accEMCal-acc*jetarea1-acc1*jetarea2;
 294   if(debug)cout<<"pt sum   "<<sumPt<<" area  "<<areasum<<endl;
 295
 296   if(nJ>0) rhoback=sumPt/areasum;
 297   else rhoback=0.;
 298   if(debug)cout<<" rho from leading jet paricle array removed   "<<rhoback<<endl;
 299
 300   return rhoback;
 301
 302 }
 303
 304
 305
 306
 307 ////////////////////////////////////////////////////////////////////////
 308 Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
 309 {
 310
 311   // Cone background subtraction method applied on the fastjet: REmove the particles of the
 312   // two largest jets with the given R from the estimation of new rho.
 313
 314   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 315   Bool_t debug  = header->GetDebug();     // debug option
 316
 317   if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
 318
 319   Float_t rc= header->GetRparam();
 320
 321   //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
 322   Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
 323
 324   Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
 325   if(debug)cout<<"nJets:  "<<nJ<<endl;
 326
 327
 328   //begin unit array
 329   TClonesArray* fUnit = fReader->GetUnitArray();
 330   if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
 331
 332   Int_t nIn = fUnit->GetEntries();
 333   if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
 334
 335   // Information extracted from fUnitArray
 336   // load input vectors and calculate total energy in array
 337   Float_t pt,eta,phi;
 338   Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
 339   Float_t rhoback=0.0;
 340
 341   Float_t ptallback=0.0; //particles without the jet
 342   Float_t restarea=accEMCal; //initial area set
 343   Bool_t acc=0;
 344   Bool_t acc1=0;
 345   Float_t rCone=0.4;
 346
 347   if(nJ==1) {
 348     AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
 349     jeteta=jettmp->Eta();
 350     jetphi=jettmp->Phi();
 351     acc=EmcalAcceptance(jeteta,jetphi,rCone);
 352     if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
 353     if(debug)cout<<" acc  "<<acc<<endl;
 354   }
 355
 356
 357   if(nJ>=2) {
 358     AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
 359     AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
 360     jeteta=jettmp->Eta();
 361     jetphi=jettmp->Phi();
 362     jeteta1=jettmp1->Eta();
 363     jetphi1=jettmp1->Phi();
 364     acc=EmcalAcceptance(jeteta,jetphi,rCone);
 365     acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
 366     if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
 367     if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
 368     if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
 369
 370     if(debug)cout<<" acc1="<<acc<<"  acc2="<<acc1<<"  restarea="<<restarea<<endl;
 371
 372   }
 373
 374   // cout<<" nIn = "<<nIn<<endl;
 375   Float_t sumpt=0;
 376   for(Int_t i=0; i<nIn; i++)
 377     { //Unit Array Loop
 378       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 379       if(uArray->GetUnitEnergy()>0.){
 380
 381         pt    = uArray->GetUnitEnergy();
 382         eta   = uArray->GetUnitEta();
 383         phi   = uArray->GetUnitPhi();
 384
 385         //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
 386
 387         Float_t deta=0.0, dphi=0.0, dr=100.0;
 388         Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
 389
 390         //cout<<i<<"  pt="<<pt<<"  eta="<<eta<<"  phi="<<phi<<endl;
 391         if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
 392           sumpt+=pt;
 393           //if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
 394
 395           if(nJ==1 && acc==1) {
 396             deta = eta - jeteta;
 397             dphi = phi - jetphi;
 398             if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 399             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 400             dr = TMath::Sqrt(deta * deta + dphi * dphi);
 401             if(dr<=rc)sumpt-=pt;
 402           }
 403
 404           if(nJ>=2) {
 405             if(acc==1){
 406               deta = eta - jeteta;
 407               dphi = phi - jetphi;
 408               if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 409               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 410               dr = TMath::Sqrt(deta * deta + dphi * dphi);
 411               if(dr<=rc)sumpt-=pt;
 412             }
 413             if(acc1==1){
 414               deta1 = eta - jeteta1;
 415               dphi1 = phi - jetphi1;
 416               if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
 417               if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
 418               dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
 419               if(dr1<=rc)sumpt-=pt;
 420             }
 421
 422           }
 423
 424           if(dr >= rc && dr1 >=rc) {
 425             // particles outside both cones
 426
 427             //cout<<" out of the cone  "<<dr<<"    "<<deta<<"  deltaeta  "<<dphi<<"  dphi "<<i<<"  particle  "<<endl;
 428             //cout<<" out of the cone  "<<dr1<<"    "<<deta1<<"  deltaeta1  "<<dphi1<<"  dphi1 "<<i<<"  particle  "<<endl;
 429             ptallback+=pt;
 430           }
 431         }
 432         //cout<<"  ipart  "<<ipart<<"  rhointegral  "<<rhoback <<endl;
 433       }
 434     } // End loop on UnitArray
 435
 436   if(debug)cout<<"total area left "<<restarea<<endl;
 437   if(debug)cout<<"sumpt="<<sumpt<<endl;
 438   // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
 439   //else rhoback=ptallback;
 440
 441   rhoback= ptallback/restarea;
 442   if(debug)cout<<"rhoback    "<<rhoback<<"     "<<nJ<<"   "<<endl;
 443
 444   return rhoback;
 445
 446 }
 447
 448
 449 Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
 450   //calculate rho using the fastjet method
 451
 452   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 453   Bool_t debug  = header->GetDebug();     // debug option
 454
 455   fastjet::Strategy strategy = header->GetStrategy();
 456   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
 457   fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
 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)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){
 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) {
 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
 516   fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 517
 518
 519   Double_t rho=clust_seq.median_pt_per_unit_area(range);
 520   //    double median, sigma, meanArea;
 521   //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
 522   //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
 523
 524   // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
 525
 526
 527   if(debug)cout<<"bkg in R="<<rParamBkg<<"  : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" --  phi="<<phiMin<<","<<phiMax<<endl;
 528   return rho;
 529 }
 530
 531 Float_t  AliJetBkg::EtaToTheta(Float_t arg)
 532 {
 533   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 534   return 2.*atan(exp(-arg));
 535 }
 536
 537
 538 Double_t  AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
 539   //to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
 540   return 1;
 541 }
 542
 543
 544 Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
 545
 546   Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
 547   Float_t deltaphi=110./180.*TMath::Pi();
 548   Float_t phicut=deltaphi/2.-radius;
 549   Float_t etacut=0.7-radius;
 550   //cout<<"  eta    "<<eta<<"  phi    "<<phi<<endl;
 551   //cout<<etacut<<"    "<<phicut<<"    "<<meanPhi<<"    "<<endl;
 552   if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
 553   else return 0;
 554 }