[u/mrichter/AliRoot.git] / JETAN / AliFastJetBkg.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//---------------------------------------------------------------------
// Class to calculate the background per unit area
// manages the search for jets
// Authors: Elena Bruna elena.bruna@yale.edu
//          Sevil Salur ssalur@lbl.gov
//
// 2011 :
// renamed from AliJetBkg to AliFastJetBkg as this class uses only FASTJET based algos.
//---------------------------------------------------------------------

#include <Riostream.h> 
#include <TClonesArray.h>
#include <TF1.h>
#include <TString.h>

#include "AliJetHeader.h"
#include "AliFastJetHeaderV1.h"
#include "AliAODJet.h"
#include "AliFastJetInput.h"
#include "AliFastJetBkg.h"

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/JetDefinition.hh"

#include<vector> 

using namespace std;

ClassImp(AliFastJetBkg)

////////////////////////////////////////////////////////////////////////

AliFastJetBkg::AliFastJetBkg():
  TObject(),
  fHeader(0),
  fInputFJ(0)
{
  // Default constructor
}

//______________________________________________________________________
AliFastJetBkg::AliFastJetBkg(const AliFastJetBkg& input):
  TObject(input),
  fHeader(input.fHeader),
  fInputFJ(input.fInputFJ)
{
  // copy constructor
}

//______________________________________________________________________
AliFastJetBkg& AliFastJetBkg::operator=(const AliFastJetBkg& source)
{
  // Assignment operator. 
  if(this!=&source){ 
   TObject::operator=(source);
   fHeader = source.fHeader;
   fInputFJ = source.fInputFJ;
  }
  
  return *this;

}

//___________________________________________________________________
void AliFastJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma, 
				Double_t& meanarea)
{
  // Bkg estimation
   
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; 
  Int_t debug  = header->GetDebug();     // debug option
  if(debug>0) cout<<"===============  AliFastJetBkg::BkgFastJetb()  =========== "<<endl;
  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  
  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation

  Double_t medianb,sigmab,meanareab;
  CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
  rho=medianb;
  sigma=sigmab;
  meanarea=meanareab;
 
}

//_________________________________________________________________
void AliFastJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma, 
					Double_t& meanarea)
{

  // Bkg estimation without hardest jet

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; 
  Int_t debug  = header->GetDebug();     // debug option
  if(debug) cout<<"===============  AliFastJetBkg::BkgWoHardest()  =========== "<<endl;
  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  
  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation  
  Double_t medianb,sigmab,meanareab;
  CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
  rho=medianb;
  sigma=sigmab;
  meanarea=meanareab;

}

//____________________________________________________________________
void AliFastJetBkg::CalcRhob(Double_t& median,Double_t& 
			     sigma,Double_t& 
			     meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t 
			     rParamBkg,TString method)
{
  // calculate rho using the fastjet method

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  fastjet::Strategy strategy = header->GetStrategy();
  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
  fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm();
  fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);

  // create an object that specifies how we to define the area
  fastjet::AreaDefinition areaDef;
  double ghostEtamax = header->GetGhostEtaMax(); 
  double ghostArea   = header->GetGhostArea(); 
  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 

  // now create the object that holds info about ghosts

  fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
  // and from that get an area definition
  fastjet::AreaType areaType = header->GetAreaType();
  areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
  
  //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
  TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
  comment+= "Jet definition: ";
  comment+= TString(jetDef.description());
  // comment+= ". Area definition: ";
  // comment+= TString(areaDef.description());
  comment+= ". Strategy adopted by FastJet: ";
  comment+= TString(clust_seq.strategy_string());
  comment+= Form("Method: %s",method.Data());
  header->SetComment(comment);
  if(debug>0){
    cout << "--------------------------------------------------------" << endl;
    cout << comment << endl;
    cout << "--------------------------------------------------------" << endl;
  }

  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); 

  double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;

  phiMin = 0;
  phiMax = 2*TMath::Pi();
 
  rapMax = ghostEtamax - rParamBkg;
  rapMin = - ghostEtamax + rParamBkg;

  fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);

  double medianb, sigmab, meanareab;
  clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false);
  median=medianb;
  sigma=sigmab;
  meanarea=meanareab; 
 
}

//____________________________________________________________________
void AliFastJetBkg::CalcRhoWoHardest(Double_t& median,Double_t& 
				     sigma,Double_t& meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t 
				     rParamBkg,TString method)
{
  // calculate rho (without the hardest jet) using the fastjet method

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  fastjet::Strategy strategy = header->GetStrategy();
  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
  fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm();
  fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);

  // create an object that specifies how we to define the area
  fastjet::AreaDefinition areaDef;
  double ghostEtamax = header->GetGhostEtaMax(); 
  double ghostArea   = header->GetGhostArea(); 
  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 

  // now create the object that holds info about ghosts

  fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
  // and from that get an area definition
  fastjet::AreaType areaType = header->GetAreaType();
  areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
  //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
  fastjet::ClusterSequenceArea clust_seq(inputParticles,jetDef,areaDef);
  TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
  comment+= "Jet definition: ";
  comment+= TString(jetDef.description());
  // comment+= ". Area definition: ";
  // comment+= TString(areaDef.description());
  comment+= ". Strategy adopted by FastJet: ";
  comment+= TString(clust_seq.strategy_string());
  comment+= Form("Method: %s",method.Data());
  header->SetComment(comment);
if(debug>0){
    cout << "--------------------------------------------------------" << endl;
    cout << comment << endl;
    cout << "--------------------------------------------------------" << endl;
  }

  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); 
  vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets);
  if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1);
    
  double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;

  phiMin = 0;
  phiMax = 2*TMath::Pi();
 
  rapMax = ghostEtamax - rParamBkg;
  rapMin = - ghostEtamax + rParamBkg;

  fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);

  double medianb, sigmab, meanareab;
  clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab, 
				     meanareab, false);
  median=medianb;
  sigma=sigmab;
  meanarea=meanareab; 

}

//___________________________________________________________________
Float_t AliFastJetBkg::BkgFastJet()
{
  // Return background  
 
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  if(debug>0) cout<<"===============  AliFastJetBkg::BkgFastJet()  =========== "<<endl;
  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  
  if(debug>0) cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
  
  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
  Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
  if(debug) cout<<"-------- rho (from all part)="<<rho<<endl; 
  return rho;
 
}

//___________________________________________________________________
Float_t AliFastJetBkg::BkgChargedFastJet()
{
  // Background for charged jets
  
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  if(debug>0) cout<<"===============  AliFastJetBkg::BkgChargedFastJet()  =========== "<<endl;

  vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
  
  if(debug>0) cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;

  double rParam = header->GetRparam();

  Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");

  if(debug>0) cout<<"-------- rho (from CHARGED part)="<<rho<<endl; 
  return rho;

}

//___________________________________________________________________
Float_t AliFastJetBkg::BkgStat()
{
  // background subtraction using statistical method
 
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  if(debug>0) cout<<"==============AliFastJetBkg::BkgStat()============="<<endl;
  //TO BE IMPLEMENTED 
  Int_t nTracks= 0;
  TF1 fun("fun",BkgFunction,0,800,1);
  Double_t enTot=fun.Eval(nTracks);
  Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
  return enTot/accEMCal;

}

//___________________________________________________________________
Float_t AliFastJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
{
  // Cone background subtraction method applied on the fastjet: REmove the particles of the
  // two largest jets with the given R from the estimation of new rho. 

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  if(debug>0) cout<<"==============AliFastJetBkg::SubtractFastJetBackgCone()============="<<endl;

  Float_t rc= header->GetRparam();

  //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
  Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal

  Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets... 
  if(debug>0) cout<<"nJets:  "<<nJ<<endl;
  
  // Information extracted from fInputParticle
  // load input vectors and calculate total energy in array
  Float_t pt,eta,phi;
  Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
  Float_t rhoback=0.0;

  Float_t ptallback=0.0; //particles without the jet
  Float_t restarea=accEMCal; //initial area set 
  Bool_t acc=0;
  Bool_t acc1=0;
  Float_t rCone=0.4;
  
  if(nJ==1) { 
    AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
    jeteta=jettmp->Eta();
    jetphi=jettmp->Phi();
    acc=EmcalAcceptance(jeteta,jetphi,rCone);
    if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
    if(debug) cout<<" acc  "<<acc<<endl;
  }
  
  if(nJ>=2) { 
    AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
    AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1));
    jeteta=jettmp->Eta();
    jetphi=jettmp->Phi();
    jeteta1=jettmp1->Eta();
    jetphi1=jettmp1->Phi(); 
    acc=EmcalAcceptance(jeteta,jetphi,rCone);
    acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
    if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
    if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
    if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;

    if(debug) cout<<" acc1="<<acc<<"  acc2="<<acc1<<"  restarea="<<restarea<<endl;

  }
  
  // cout<<" nIn = "<<nIn<<endl;
  Float_t sumpt=0;
  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  for(UInt_t i=0; i<inputParticles.size(); i++)
    { // Loop over input list of particles
      pt    = inputParticles[i].perp();
      eta   = inputParticles[i].eta();
      phi   = inputParticles[i].phi();

      // To be updated
      //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
	
      Float_t deta=0.0, dphi=0.0, dr=100.0;
      Float_t deta1=0.0, dphi1=0.0, dr1=100.0;

      //cout<<i<<"  pt="<<pt<<"  eta="<<eta<<"  phi="<<phi<<endl;
      if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
	sumpt+=pt;
	//if(i<30)cout<<i<<" pt = "<<pt<<endl;

	if(nJ==1 && acc==1) { 
	  deta = eta - jeteta;
	  dphi = phi - jetphi;
	  if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
	  if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
	  dr = TMath::Sqrt(deta * deta + dphi * dphi);
	  if(dr<=rc)sumpt-=pt;
	}
	
	if(nJ>=2) { 
	  if(acc==1){
	    deta = eta - jeteta;
	    dphi = phi - jetphi;
	    if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
	    if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
	    dr = TMath::Sqrt(deta * deta + dphi * dphi);
	    if(dr<=rc)sumpt-=pt;
	  }
	  if(acc1==1){
	    deta1 = eta - jeteta1;
	    dphi1 = phi - jetphi1;
	    if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
	    if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
	    dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
	    if(dr1<=rc)sumpt-=pt;
	  }
	}

	if(dr >= rc && dr1 >=rc) { 
	  // particles outside both cones
	  if(debug>1) cout<<" out of the cone  "<<dr<<"    "<<deta<<"  deltaeta  "<<dphi<<"  dphi "<<i<<"  particle  "<<endl;
	  if(debug>1) cout<<" out of the cone  "<<dr1<<"    "<<deta1<<"  deltaeta1  "<<dphi1<<"  dphi1 "<<i<<"  particle  "<<endl;
	  ptallback+=pt;
	}
      }
    } // End loop on input list of particles 
  
  if(debug>0) cout<<"total area left "<<restarea<<endl;
  if(debug>0) cout<<"sumpt="<<sumpt<<endl;
  // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
  // else rhoback=ptallback;

  rhoback= ptallback/restarea;
  if(debug)cout<<"rhoback    "<<rhoback<<"     "<<nJ<<"   "<<endl;

  return rhoback;
   
}

//___________________________________________________________________
Double_t AliFastJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method)
{
  // calculate rho using the fastjet method

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  Int_t debug  = header->GetDebug();     // debug option

  fastjet::Strategy strategy = header->GetStrategy();
  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
  fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm();
  fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);

  // create an object that specifies how we to define the area
  fastjet::AreaDefinition areaDef;
  double ghostEtamax       = header->GetGhostEtaMax(); 
  double ghostArea         = header->GetGhostArea(); 
  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 
  
  // now create the object that holds info about ghosts

  if (method.Contains("Charg"))ghostEtamax=0.9;

  fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
  // and from that get an area definition
  fastjet::AreaType areaType = header->GetAreaType();
  areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
  if(debug>0) cout<<"rParamBkg="<<rParamBkg<<"  ghostEtamax="<<ghostEtamax<<"  ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<< endl;
  //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
  TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
  comment+= "Jet definition: ";
  comment+= TString(jetDef.description());
  // comment+= ". Area definition: ";
  // comment+= TString(areaDef.description());
  comment+= ". Strategy adopted by FastJet: ";
  comment+= TString(clust_seq.strategy_string());
  header->SetComment(comment);
  if(debug>0){
    cout << "--------------------------------------------------------" << endl;
    cout << comment << endl;
    cout << "--------------------------------------------------------" << endl;
  }

  double ptmin = header->GetPtMin(); 
  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); 

  if (debug>0) {
    cout<<"# of BKG jets = "<<jets.size()<<endl;
    for (size_t j = 0; j < jets.size(); j++) { // loop for jets   
      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]));
    }
  }
  
  double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;

  if (method.Contains("All")){
    phiMin = 80.*TMath::Pi()/180+rParamBkg;
    phiMax = 190.*TMath::Pi()/180-rParamBkg;
  }
  if (method.Contains("Charg")){
    phiMin = 0;
    phiMax = 2*TMath::Pi();
  }
  rapMax = ghostEtamax - rParamBkg;
  rapMin = - ghostEtamax + rParamBkg;

  fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);

  Double_t rho=clust_seq.median_pt_per_unit_area(range);
  // double median, sigma, meanArea;
  // clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
  // fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);

  // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);

  if(debug>0) cout<<"bkg in R="<<rParamBkg<<"  : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" --  phi="<<phiMin<<","<<phiMax<<endl;

  return rho;

}

//___________________________________________________________________
Double_t  AliFastJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/)
{
  // to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
  return 1;

}

//___________________________________________________________________
Bool_t AliFastJetBkg::EmcalAcceptance(Float_t eta, Float_t phi, Float_t radius) const
{
  // Apply emcal acceptance cuts
  // To be updated

  Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
  Float_t deltaphi=110./180.*TMath::Pi();
  Float_t phicut=deltaphi/2.-radius;
  Float_t etacut=0.7-radius;
  //cout<<"  eta    "<<eta<<"  phi    "<<phi<<endl;
  //cout<<etacut<<"    "<<phicut<<"    "<<meanPhi<<"    "<<endl;
  if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
  else return 0; 

}
Commit	Line	Data
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
	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;
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	//_________________________________________________________________
	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;
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	//____________________________________________________________________
	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;
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	//____________________________________________________________________
	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;
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);
	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;
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	//___________________________________________________________________
	281	Float_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	//___________________________________________________________________
	304	Float_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=20.7110./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;
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=20.7110./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()rcrc;
	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-2TMath::Pi()rc*rc;
	371	if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()rcrc;
	372	if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()rcrc;
	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;
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	//___________________________________________________________________
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	//___________________________________________________________________
7aa72c24	539	Bool_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	}