Fix for Coverity defect 22028
[u/mrichter/AliRoot.git] / JETAN / AliFastJetBkg.cxx
CommitLineData
d89b8229 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
46using namespace std;
47
48ClassImp(AliFastJetBkg)
49
50////////////////////////////////////////////////////////////////////////
51
52AliFastJetBkg::AliFastJetBkg():
53 TObject(),
54 fHeader(0),
55 fInputFJ(0)
56{
57 // Default constructor
58}
59
60//______________________________________________________________________
61AliFastJetBkg::AliFastJetBkg(const AliFastJetBkg& input):
62 TObject(input),
63 fHeader(input.fHeader),
64 fInputFJ(input.fInputFJ)
65{
66 // copy constructor
67}
68
69//______________________________________________________________________
70AliFastJetBkg& 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//___________________________________________________________________
84void AliFastJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma,
85 Double_t& meanarea)
86{
87 // Bkg estimation
88
89 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
37f45f89 90 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//_________________________________________________________________
105void 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;
37f45f89 112 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//____________________________________________________________________
126void 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;
37f45f89 134 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//____________________________________________________________________
193void 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;
37f45f89 200 Int_t debug = header->GetDebug(); // debug option
d89b8229 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);
230if(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//___________________________________________________________________
261Float_t AliFastJetBkg::BkgFastJet()
262{
263 // Return background
264
265 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
37f45f89 266 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//___________________________________________________________________
281Float_t AliFastJetBkg::BkgChargedFastJet()
282{
283 // Background for charged jets
284
285 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
37f45f89 286 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//___________________________________________________________________
304Float_t AliFastJetBkg::BkgStat()
305{
306 // background subtraction using statistical method
307
308 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
37f45f89 309 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//___________________________________________________________________
322Float_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;
37f45f89 328 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//___________________________________________________________________
448Double_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;
37f45f89 453 Int_t debug = header->GetDebug(); // debug option
d89b8229 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//___________________________________________________________________
531Double_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//___________________________________________________________________
7aa72c24 539Bool_t AliFastJetBkg::EmcalAcceptance(Float_t eta, Float_t phi, Float_t radius) const
d89b8229 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}