Adding rho's depenence on trigger track (M. Verweij)
[u/mrichter/AliRoot.git] / 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 }