[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.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.                  *
 **************************************************************************/

#include <Riostream.h> 
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <TArrayF.h>
#include <TClonesArray.h>
#include <TF1.h>
#include <TString.h>

#include "AliJetHeader.h"
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliFastJetFinder.h"
#include "AliFastJetHeaderV1.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
#include "AliJetUnitArray.h"
#include "AliFastJetInput.h"
#include "AliESDEvent.h"


#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/JetDefinition.hh"
// get info on how fastjet was configured
#include "fastjet/config.h"

#ifdef ENABLE_PLUGIN_SISCONE
#include "fastjet/SISConePlugin.hh"
#endif

#include<sstream>  // needed for internal io
#include<vector> 
#include <cmath> 

using namespace std;
#include "AliJetBkg.h"

ClassImp(AliJetBkg)

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

AliJetBkg::AliJetBkg():
  TObject(),
  fReader(0),
  fHeader(0),
  fInputFJ(0)
{
  // Default constructor

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

}
//______________________________________________________________________
AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
    // Assignment operator. 
    this->~AliJetBkg();
    new(this) AliJetBkg(source);
    return *this;

}
//___________________________________________________________________
Float_t AliJetBkg::BkgFastJet(){
  
  cout<<"===============  AliJetBkg::BkgFastJet()  =========== "<<endl;
  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  
  cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
  
  for(UInt_t i=0;i<inputParticles.size();i++){
    //  cout<<"                "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
    
  }
   
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
  Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
  cout<<"-------- rho (from all part)="<<rho<<endl; 
  return rho;
 
}
//___________________________________________________________________
Float_t  AliJetBkg::BkgChargedFastJet(){

  cout<<"===============  AliJetBkg::BkgChargedFastJet()  =========== "<<endl;

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

  for(UInt_t i=0;i<inputParticlesCharged.size();i++){
    // cout<<"                "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;

  } 
  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  double rParam = header->GetRparam();

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

  cout<<"-------- rho (from CHARGED part)="<<rho<<endl; 
  return rho;
}


Float_t AliJetBkg::BkgStat()
{
  //background subtraction using statistical method
 
  cout<<"==============AliJetBkg::BkgStat()============="<<endl;
  //TO BE IMPLEMENTED 
  // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
  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 AliJetBkg::BkgRemoveJetLeading(TClonesArray* fAODJets)
{
  // Remove the particles of the
  // two largest jets using the track references  stored in the AODJet from the estimation of new rho. 

  cout<<"==============AliJetBkg::BkgRemoveJetLeading()============="<<endl;

 
  // check if we are reading AOD jets
  TRefArray *refs = 0;
  Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetESDReader");
  if (fromAod) { refs = fReader->GetReferences(); }
  
  //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... 
  cout<<"nJets:  "<<nJ<<endl;
  

  //begin unit array 
  TClonesArray* fUnit = fReader->GetUnitArray();
  if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }

  Int_t nIn = fUnit->GetEntries();
  if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
  
  Float_t rhoback=0.0;
  Float_t jetarea1=0.0,jetarea2=0.0;
    
  Int_t particlejet1=-99;
  Int_t particlejet2=-99;
  TRefArray *refarray1 = 0;
  TRefArray *refarray2 = 0;
  Int_t nJettracks1 = 0, nJettracks2 = 0;
  Int_t acc=0,acc1=0;
  AliAODJet *jet1;
  AliAODJet *jet2;

  if(nJ==1){
    jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
    jetarea1=jet1->EffectiveAreaCharged();
    Float_t jetPhi=jet1->Phi();
    Float_t jetEta=jet1->Eta();
    if(jetPhi>1.396 && jetPhi<3.316 && jetEta>-0.7 && jetEta<0.7)acc=1;
    refarray1=jet1->GetRefTracks();
    nJettracks1=refarray1->GetEntries();
    cout<<"nJ = 1, acc="<<acc<<"  jetarea1="<<jetarea1<<endl;

  }

  if(nJ>=2){
    jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
    jetarea1=jet1->EffectiveAreaCharged();
    Float_t jetPhi1=jet1->Phi();
    Float_t jetEta1=jet1->Eta();
    if(jetPhi1>1.396 && jetPhi1<3.316 && jetEta1>-0.7 && jetEta1<0.7)acc=1;
    refarray1=jet1->GetRefTracks();
    nJettracks1=refarray1->GetEntries();
    cout<<"npart = "<<nJettracks1<<endl;
   
    jet2 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
    jetarea2=jet2->EffectiveAreaCharged();
    Float_t jetPhi2=jet2->Phi();
    Float_t jetEta2=jet2->Eta();
    if(jetPhi2>1.396 && jetPhi2<3.316 && jetEta2>-0.7 && jetEta2<0.7)acc1=1;
    refarray2=jet2->GetRefTracks();
    nJettracks2=refarray2->GetEntries();
    cout<<"nJ = "<<nJ<<", acc="<<acc<<"  acc1="<<acc1<<"  jetarea1="<<jetarea1<<"  jetarea2="<<jetarea2<<endl;
  }


  // cout<<" nIn = "<<nIn<<endl;
  Float_t sumPt=0;
  Float_t eta,phi,pt;
  Int_t ipart=0;

  for(Int_t i=0; i<nIn; i++) 
    { //Unit Array Loop
      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);

      if(uArray->GetUnitEnergy()>0.){
	eta   = uArray->GetUnitEta();
	phi   = uArray->GetUnitPhi();
	pt = uArray->GetUnitEnergy();
	//	cout<<"ipart = "<<ipart<<" eta="<<eta<<"  phi="<<phi<<endl;
 	if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
	  //cout<<sumPt<<endl;
 	    sumPt+=pt;

	    
	    if(nJ==1 && acc==1){
	      for(Int_t ii=0; ii<nJettracks1;ii++){    
		
		particlejet1  = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
	 
		if(ipart==particlejet1) {
		  sumPt-=pt;
		}
	      }
	    }
	  
	  
	    if(nJ>=2){
	      
	      //first jet
	      if(acc==1){
		for(Int_t ii=0; ii<nJettracks1;ii++){    
		  particlejet1  = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
		  
		//cout<<"uArr loop = "<<i<<"  ipart in uArr (1/2)="<<ipart<<"  part in jet="<<ii<<"  partID="<<particlejet1<<" sumPt="<<sumPt<<endl; 
		  if(ipart==particlejet1) {
		    sumPt-=pt;
		  }
		}
	      }
	      if(acc1==1){
		//second jet
		for(Int_t ii=0; ii<nJettracks2;ii++){ 
		  particlejet2  = ((AliJetUnitArray*)refarray2->At(ii))->GetUnitTrackID();
		  //cout<<"uArr loop = "<<i<<"  ipart in uArr (2/2)="<<ipart<<"  part in jet="<<ii<<"  partID="<<particlejet2<<" sumPt="<<sumPt<<endl; 
		  if(ipart==particlejet2) {
		    sumPt-=pt;
		  }
		}
	      }
	      
	      
	    }
	
	}//if phi,eta
	ipart++;
      }//end if energy
    }// end unit array loop	      
  

  Float_t areasum=areasum=accEMCal-acc*jetarea1-acc1*jetarea2;
  cout<<"pt sum   "<<sumPt<<" area  "<<areasum<<endl;
   
  if(nJ>0) rhoback=sumPt/areasum;
  else rhoback=0.;
  cout<<" rho from leading jet paricle array removed   "<<rhoback<<endl;
 
  return rhoback;
 
}


////////////////////////////////////////////////////////////////////////
Float_t AliJetBkg::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. 
  cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;

  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  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... 
  cout<<"nJets:  "<<nJ<<endl;
  

  //begin unit array 
  TClonesArray* fUnit = fReader->GetUnitArray();
  if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }

  Int_t nIn = fUnit->GetEntries();
  if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
  
  // Information extracted from fUnitArray
  // 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 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
    jeteta=jettmp->Eta();
    jetphi=jettmp->Phi();
    acc=EmcalAcceptance(jeteta,jetphi,rCone);
    if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
    cout<<" acc  "<<acc<<endl;
  }

  
  if(nJ>=2) { 
    AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
    AliAODJet *jettmp1 = dynamic_cast<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;

    cout<<" acc1="<<acc<<"  acc2="<<acc1<<"  restarea="<<restarea<<endl;

  }
  
  // cout<<" nIn = "<<nIn<<endl;
  Float_t sumpt=0;
  for(Int_t i=0; i<nIn; i++) 
    { //Unit Array Loop
      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
      if(uArray->GetUnitEnergy()>0.){
	
	pt    = uArray->GetUnitEnergy();
       	eta   = uArray->GetUnitEta();
	phi   = uArray->GetUnitPhi();

	//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 unit = "<<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
	  
	    //cout<<" out of the cone  "<<dr<<"    "<<deta<<"  deltaeta  "<<dphi<<"  dphi "<<i<<"  particle  "<<endl;
	    //cout<<" out of the cone  "<<dr1<<"    "<<deta1<<"  deltaeta1  "<<dphi1<<"  dphi1 "<<i<<"  particle  "<<endl;
	    ptallback+=pt;
	  }
	}
	//cout<<"  ipart  "<<ipart<<"  rhointegral  "<<rhoback <<endl;
      }		
    } // End loop on UnitArray 
  
  cout<<"total area left "<<restarea<<endl;
  cout<<"sumpt="<<sumpt<<endl;
  // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
  //else rhoback=ptallback;

  rhoback= ptallback/restarea;
  cout<<"rhoback    "<<rhoback<<"     "<<nJ<<"   "<<endl;

  return rhoback;
   
}


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

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

  fastjet::Strategy strategy = header->GetStrategy();
  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
  fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
  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);
  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){
    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); 
  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 phiMax = header->GetPhiMax();
  //double phiMin = header->GetPhiMin();
  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;
    //phiMin = 0;
    //phiMax = 2*TMath::Pi();
  }
  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);


  cout<<"bkg in R="<<rParamBkg<<"  : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" --  phi="<<phiMin<<","<<phiMax<<endl;
  return rho;
}

Float_t  AliJetBkg::EtaToTheta(Float_t arg)
{
  //  return (180./TMath::Pi())*2.*atan(exp(-arg));
  return 2.*atan(exp(-arg));
}


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


Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
 
  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
50c7d9f7	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():
856618e7	60	TObject(),
	61	fReader(0),
	62	fHeader(0),
	63	fInputFJ(0)
50c7d9f7	64	{
	65	// Default constructor
	66
	67	}
	68	//______________________________________________________________________
	69	AliJetBkg::AliJetBkg(const AliJetBkg& input):
856618e7	70	TObject(input),
50c7d9f7	71	fReader(input.fReader),
50c7d9f7	72	fHeader(input.fHeader),
856618e7	73	fInputFJ(input.fInputFJ)
50c7d9f7	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	Float_t AliJetBkg::BkgFastJet(){
	88
	89	cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
	90	vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
	91
	92	cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
	93
	94	for(UInt_t i=0;i<inputParticles.size();i++){
	95	// cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
	96
	97	}
	98
	99	AliFastJetHeaderV1 header = (AliFastJetHeaderV1)fHeader;
	100	double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
	101	Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
	102	cout<<"-------- rho (from all part)="<<rho<<endl;
	103	return rho;
	104
	105	}
	106	//___________________________________________________________________
	107	Float_t AliJetBkg::BkgChargedFastJet(){
	108
	109	cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
	110
	111	vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
	112
	113	cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
	114
	115	for(UInt_t i=0;i<inputParticlesCharged.size();i++){
	116	// cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
	117
	118	}
	119	AliFastJetHeaderV1 header = (AliFastJetHeaderV1)fHeader;
	120	double rParam = header->GetRparam();
	121
	122	Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
	123
	124	cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
	125	return rho;
	126	}
	127
	128
	129
	130	Float_t AliJetBkg::BkgStat()
	131	{
	132	//background subtraction using statistical method
	133
	134	cout<<"==============AliJetBkg::BkgStat()============="<<endl;
	135	//TO BE IMPLEMENTED
	136	// Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
	137	Int_t nTracks= 0;
138	TF1 fun("fun",BkgFunction,0,800,1);
139	Double_t enTot=fun.Eval(nTracks);
140	Double_t accEMCal=20.7110./180*TMath::Pi();//2.68 area of EMCal
141	return enTot/accEMCal;
142
143	}
144	/////////////////////////////////
145	Float_t AliJetBkg::BkgRemoveJetLeading(TClonesArray* fAODJets)
146	{
147	// Remove the particles of the
148	// two largest jets using the track references stored in the AODJet from the estimation of new rho.
149
150	cout<<"==============AliJetBkg::BkgRemoveJetLeading()============="<<endl;
151
152
153	// check if we are reading AOD jets
154	TRefArray *refs = 0;
155	Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetESDReader");
156	if (fromAod) { refs = fReader->GetReferences(); }
157
158	//Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
159	Double_t accEMCal=20.7110./180*TMath::Pi();//2.68 area of EMCal
160
161	Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
162	cout<<"nJets: "<<nJ<<endl;
163
164
165	//begin unit array
166	TClonesArray* fUnit = fReader->GetUnitArray();
167	if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
168
169	Int_t nIn = fUnit->GetEntries();
170	if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
171
172	Float_t rhoback=0.0;
173	Float_t jetarea1=0.0,jetarea2=0.0;
174
175	Int_t particlejet1=-99;
176	Int_t particlejet2=-99;
177	TRefArray *refarray1 = 0;
178	TRefArray *refarray2 = 0;
179	Int_t nJettracks1 = 0, nJettracks2 = 0;
180	Int_t acc=0,acc1=0;
181	AliAODJet *jet1;
182	AliAODJet *jet2;
183
184	if(nJ==1){
185	jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
186	jetarea1=jet1->EffectiveAreaCharged();
187	Float_t jetPhi=jet1->Phi();
188	Float_t jetEta=jet1->Eta();
189	if(jetPhi>1.396 && jetPhi<3.316 && jetEta>-0.7 && jetEta<0.7)acc=1;
190	refarray1=jet1->GetRefTracks();
191	nJettracks1=refarray1->GetEntries();
192	cout<<"nJ = 1, acc="<<acc<<" jetarea1="<<jetarea1<<endl;
193
194	}
195
196	if(nJ>=2){
197	jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
198	jetarea1=jet1->EffectiveAreaCharged();
199	Float_t jetPhi1=jet1->Phi();
200	Float_t jetEta1=jet1->Eta();
201	if(jetPhi1>1.396 && jetPhi1<3.316 && jetEta1>-0.7 && jetEta1<0.7)acc=1;
202	refarray1=jet1->GetRefTracks();
203	nJettracks1=refarray1->GetEntries();
204	cout<<"npart = "<<nJettracks1<<endl;
205
206	jet2 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
207	jetarea2=jet2->EffectiveAreaCharged();
208	Float_t jetPhi2=jet2->Phi();
209	Float_t jetEta2=jet2->Eta();
210	if(jetPhi2>1.396 && jetPhi2<3.316 && jetEta2>-0.7 && jetEta2<0.7)acc1=1;
211	refarray2=jet2->GetRefTracks();
212	nJettracks2=refarray2->GetEntries();
213	cout<<"nJ = "<<nJ<<", acc="<<acc<<" acc1="<<acc1<<" jetarea1="<<jetarea1<<" jetarea2="<<jetarea2<<endl;
214	}
215
216
217
218	// cout<<" nIn = "<<nIn<<endl;
219	Float_t sumPt=0;
220	Float_t eta,phi,pt;
221	Int_t ipart=0;
222
223	for(Int_t i=0; i<nIn; i++)
224	{ //Unit Array Loop
225	AliJetUnitArray uArray = (AliJetUnitArray)fUnit->At(i);
226
227	if(uArray->GetUnitEnergy()>0.){
228	eta = uArray->GetUnitEta();
229	phi = uArray->GetUnitPhi();
230	pt = uArray->GetUnitEnergy();
231	// cout<<"ipart = "<<ipart<<" eta="<<eta<<" phi="<<phi<<endl;
232	if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
233	//cout<<sumPt<<endl;
234	sumPt+=pt;
235
236
237	if(nJ==1 && acc==1){
238	for(Int_t ii=0; ii<nJettracks1;ii++){
239
240	particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
241
242	if(ipart==particlejet1) {
243	sumPt-=pt;
244	}
245	}
246	}
247
248
249	if(nJ>=2){
250
251	//first jet
252	if(acc==1){
253	for(Int_t ii=0; ii<nJettracks1;ii++){
254	particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
255
256	//cout<<"uArr loop = "<<i<<" ipart in uArr (1/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet1<<" sumPt="<<sumPt<<endl;
257	if(ipart==particlejet1) {
258	sumPt-=pt;
259	}
260	}
261	}
262	if(acc1==1){
263	//second jet
264	for(Int_t ii=0; ii<nJettracks2;ii++){
265	particlejet2 = ((AliJetUnitArray*)refarray2->At(ii))->GetUnitTrackID();
266	//cout<<"uArr loop = "<<i<<" ipart in uArr (2/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet2<<" sumPt="<<sumPt<<endl;
267	if(ipart==particlejet2) {
268	sumPt-=pt;
269	}
270	}
271	}
272
273
274	}
275
276	}//if phi,eta
277	ipart++;
278	}//end if energy
279	}// end unit array loop
280
281
282	Float_t areasum=areasum=accEMCal-accjetarea1-acc1jetarea2;
283	cout<<"pt sum "<<sumPt<<" area "<<areasum<<endl;
284
285	if(nJ>0) rhoback=sumPt/areasum;
286	else rhoback=0.;
287	cout<<" rho from leading jet paricle array removed "<<rhoback<<endl;
288
289	return rhoback;
290
291	}
292
293
294
295
296	////////////////////////////////////////////////////////////////////////
297	Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
298	{
299
300	// Cone background subtraction method applied on the fastjet: REmove the particles of the
301	// two largest jets with the given R from the estimation of new rho.
302	cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
303
304	AliFastJetHeaderV1 header = (AliFastJetHeaderV1)fHeader;
305	Float_t rc= header->GetRparam();
306
307	//Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
308	Double_t accEMCal=20.7110./180*TMath::Pi();//2.68 area of EMCal
309
310	Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
311	cout<<"nJets: "<<nJ<<endl;
312
313
314	//begin unit array
315	TClonesArray* fUnit = fReader->GetUnitArray();
316	if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
317
318	Int_t nIn = fUnit->GetEntries();
319	if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
320
321	// Information extracted from fUnitArray
322	// load input vectors and calculate total energy in array
323	Float_t pt,eta,phi;
324	Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
325	Float_t rhoback=0.0;
326
327	Float_t ptallback=0.0; //particles without the jet
328	Float_t restarea=accEMCal; //initial area set
329	Bool_t acc=0;
330	Bool_t acc1=0;
331	Float_t rCone=0.4;
332
333	if(nJ==1) {
334	AliAODJet jettmp = dynamic_cast<AliAODJet>(fAODJets->At(0));
335	jeteta=jettmp->Eta();
336	jetphi=jettmp->Phi();
337	acc=EmcalAcceptance(jeteta,jetphi,rCone);
338	if(acc==1)restarea= accEMCal-TMath::Pi()rcrc;
339	cout<<" acc "<<acc<<endl;
340	}
341
342
343	if(nJ>=2) {
344	AliAODJet jettmp = dynamic_cast<AliAODJet>(fAODJets->At(0));
345	AliAODJet jettmp1 = dynamic_cast<AliAODJet>(fAODJets->At(1));
346	jeteta=jettmp->Eta();
347	jetphi=jettmp->Phi();
348	jeteta1=jettmp1->Eta();
349	jetphi1=jettmp1->Phi();
350	acc=EmcalAcceptance(jeteta,jetphi,rCone);
351	acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
352	if(acc1==1 && acc==1)restarea= accEMCal-2TMath::Pi()rc*rc;
353	if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()rcrc;
354	if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()rcrc;
355
356	cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
357
358	}
359
360	// cout<<" nIn = "<<nIn<<endl;
361	Float_t sumpt=0;
362	for(Int_t i=0; i<nIn; i++)
363	{ //Unit Array Loop
364	AliJetUnitArray uArray = (AliJetUnitArray)fUnit->At(i);
365	if(uArray->GetUnitEnergy()>0.){
366
367	pt = uArray->GetUnitEnergy();
368	eta = uArray->GetUnitEta();
369	phi = uArray->GetUnitPhi();
370
371	//cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
372
373	Float_t deta=0.0, dphi=0.0, dr=100.0;
374	Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
375
376	//cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
377	if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
378	sumpt+=pt;
379	//if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
380
381	if(nJ==1 && acc==1) {
382	deta = eta - jeteta;
383	dphi = phi - jetphi;
384	if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
385	if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
386	dr = TMath::Sqrt(deta * deta + dphi * dphi);
387	if(dr<=rc)sumpt-=pt;
388	}
389
390	if(nJ>=2) {
391	if(acc==1){
392	deta = eta - jeteta;
393	dphi = phi - jetphi;
394	if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
395	if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
396	dr = TMath::Sqrt(deta * deta + dphi * dphi);
397	if(dr<=rc)sumpt-=pt;
398	}
399	if(acc1==1){
400	deta1 = eta - jeteta1;
401	dphi1 = phi - jetphi1;
402	if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
403	if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
404	dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
405	if(dr1<=rc)sumpt-=pt;
406	}
407
408	}
409
410	if(dr >= rc && dr1 >=rc) {
411	// particles outside both cones
412
413	//cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl;
414	//cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl;
415	ptallback+=pt;
416	}
417	}
418	//cout<<" ipart "<<ipart<<" rhointegral "<<rhoback <<endl;
419	}
420	} // End loop on UnitArray
421
422	cout<<"total area left "<<restarea<<endl;
423	cout<<"sumpt="<<sumpt<<endl;
424	// if(acc==1 \|\| acc1==1) rhoback= ptallback/restarea;
425	//else rhoback=ptallback;
426
427	rhoback= ptallback/restarea;
428	cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
429
430	return rhoback;
431
432	}
433
434
435	Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
436	//calculate rho using the fastjet method
437
438	AliFastJetHeaderV1 header = (AliFastJetHeaderV1)fHeader;
439	Bool_t debug = header->GetDebug(); // debug option
440
441	fastjet::Strategy strategy = header->GetStrategy();
442	fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
443	fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
444	fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
445
446	// create an object that specifies how we to define the area
447	fastjet::AreaDefinition areaDef;
448	double ghostEtamax = header->GetGhostEtaMax();
449	double ghostArea = header->GetGhostArea();
450	int activeAreaRepeats = header->GetActiveAreaRepeats();
451
452	// now create the object that holds info about ghosts
453
454	if (method.Contains("Charg"))ghostEtamax=0.9;
455
456	fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
457	// and from that get an area definition
458	fastjet::AreaType areaType = header->GetAreaType();
459	areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
460	cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
461	//fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
462	fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
463	TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
464	comment+= "Jet definition: ";
465	comment+= TString(jetDef.description());
466	// comment+= ". Area definition: ";
467	//comment+= TString(areaDef.description());
468	comment+= ". Strategy adopted by FastJet: ";
469	comment+= TString(clust_seq.strategy_string());
470	header->SetComment(comment);
471	if(debug){
472	cout << "--------------------------------------------------------" << endl;
473	cout << comment << endl;
474	cout << "--------------------------------------------------------" << endl;
475	}
476
477	double ptmin = header->GetPtMin();
478	vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
479	vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
480	cout<<"# of BKG jets = "<<jets.size()<<endl;
481	for (size_t j = 0; j < jets.size(); j++) { // loop for jets
482
483	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]));
484	}
485
486	// double phiMax = header->GetPhiMax();
487	//double phiMin = header->GetPhiMin();
488	double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
489
490	if (method.Contains("All")){
491	phiMin = 80.*TMath::Pi()/180+rParamBkg;
492	phiMax = 190.*TMath::Pi()/180-rParamBkg;
493	//phiMin = 0;
494	//phiMax = 2*TMath::Pi();
495	}
496	if (method.Contains("Charg")){
497	phiMin = 0;
498	phiMax = 2*TMath::Pi();
499	}
500	rapMax = ghostEtamax - rParamBkg;
501	rapMin = - ghostEtamax + rParamBkg;
502
503
504	fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
505
506
507	Double_t rho=clust_seq.median_pt_per_unit_area(range);
508	// double median, sigma, meanArea;
509	//clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
510	//fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
511
512	// fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
513
514
515	cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
516	return rho;
517	}
518
519	Float_t AliJetBkg::EtaToTheta(Float_t arg)
520	{
521	// return (180./TMath::Pi())2.atan(exp(-arg));
522	return 2.*atan(exp(-arg));
523	}
524
525
526	Double_t AliJetBkg::BkgFunction(Double_t /x/,Double_t /par/){
527	//to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
528	return 1;
529	}
530
531
532	Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
533
534	Float_t meanPhi=190./180.TMath::Pi()-110./180.TMath::Pi()/2;
535	Float_t deltaphi=110./180.*TMath::Pi();
536	Float_t phicut=deltaphi/2.-radius;
537	Float_t etacut=0.7-radius;
538	//cout<<" eta "<<eta<<" phi "<<phi<<endl;
539	//cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl;
540	if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
541	else return 0;
542	}