Fix for Coverity defect 22028
[u/mrichter/AliRoot.git] / 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 }