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
  88 //___________________________________________________________________
  89   void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma,
  90 Double_t& meanarea){
  91
  92     AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  93     Bool_t debug  = header->GetDebug();     // debug option
  94     if(debug)cout<<"===============  AliJetBkg::BkgFastJetb()  =========== "<<endl;
  95   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  96
  97   //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
  98
  99
 100
 101
 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,
 114 Double_t& meanarea){
 115     AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 116     Bool_t debug  = header->GetDebug();     // debug option
 117     if(debug)cout<<"===============  AliJetBkg::BkgWoHardest()  =========== "<<endl;
 118   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
 119
 120   //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
 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&
 133                           sigma,Double_t&
 134                           meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
 135                           rParamBkg,TString method){
 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());
 170   comment+= Form("Method: %s",method.Data());
 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&
 216 sigma,Double_t&
 217 meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
 218 rParamBkg,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());
 253   comment+= Form("Method: %s",method.Data());
 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,
 283 meanareab, false);
 284   median=medianb;
 285   sigma=sigmab;
 286   meanarea=meanareab;
 287
 288
 289  }
 290
 291
 292 //___________________________________________________________________
 293 Float_t AliJetBkg::BkgFastJet(){
 294 // Return background
 295   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 296   Bool_t debug  = header->GetDebug();     // debug option
 297
 298   if(debug)cout<<"===============  AliJetBkg::BkgFastJet()  =========== "<<endl;
 299   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
 300
 301   if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
 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
 308   double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
 309   Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
 310   if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
 311   return rho;
 312
 313 }
 314 //___________________________________________________________________
 315 Float_t  AliJetBkg::BkgChargedFastJet(){
 316 // Background for charged jets
 317   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 318   Bool_t debug  = header->GetDebug();     // debug option
 319
 320   if(debug)cout<<"===============  AliJetBkg::BkgChargedFastJet()  =========== "<<endl;
 321
 322   vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
 323
 324   if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
 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   }
 330
 331   double rParam = header->GetRparam();
 332
 333   Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
 334
 335   if(debug)cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
 336   return rho;
 337 }
 338
 339
 340
 341 Float_t AliJetBkg::BkgStat()
 342 {
 343   //background subtraction using statistical method
 344
 345   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 346   Bool_t debug  = header->GetDebug();     // debug option
 347
 348   if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
 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 }
 358
 359 ////////////////////////////////////////////////////////////////////////
 360 Float_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.
 365
 366   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 367   Bool_t debug  = header->GetDebug();     // debug option
 368
 369   if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
 370
 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...
 377   if(debug)cout<<"nJets:  "<<nJ<<endl;
 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) {
 400     AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
 401     jeteta=jettmp->Eta();
 402     jetphi=jettmp->Phi();
 403     acc=EmcalAcceptance(jeteta,jetphi,rCone);
 404     if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
 405     if(debug)cout<<" acc  "<<acc<<endl;
 406   }
 407
 408
 409   if(nJ>=2) {
 410     AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
 411     AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1));
 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
 422     if(debug)cout<<" acc1="<<acc<<"  acc2="<<acc1<<"  restarea="<<restarea<<endl;
 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
 488   if(debug)cout<<"total area left "<<restarea<<endl;
 489   if(debug)cout<<"sumpt="<<sumpt<<endl;
 490   // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
 491   //else rhoback=ptallback;
 492
 493   rhoback= ptallback/restarea;
 494   if(debug)cout<<"rhoback    "<<rhoback<<"     "<<nJ<<"   "<<endl;
 495
 496   return rhoback;
 497
 498 }
 499
 500
 501 Double_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);
 526   if(debug)cout<<"rParamBkg="<<rParamBkg<<"  ghostEtamax="<<ghostEtamax<<"  ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
 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);
 546
 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
 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;
 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
 579   if(debug)cout<<"bkg in R="<<rParamBkg<<"  : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" --  phi="<<phiMin<<","<<phiMax<<endl;
 580   return rho;
 581 }
 582
 583 Float_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
 590 Double_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
 596 Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius) const{
 597 // Apply emcal acceptance cuts
 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 }