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