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

//---------------------------------------------------------------------
// FastJet v2.3.4 finder algorithm interface
//
// Author: swensy.jangal@ires.in2p3.fr 
//
// Last modification: Neutral cell energy included in the jet reconstruction
//
// Author: Magali.estienne@subatech.in2p3.fr 
//---------------------------------------------------------------------


#include <Riostream.h>
#include <TArrayF.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TLorentzVector.h>
#include <TRandom.h>
#include <TClonesArray.h>

#include "AliHeader.h"
#include "AliJet.h"
#include "AliJetKineReader.h"
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliJetUnitArray.h"
#include "AliSISConeJetFinder.h"
#include "AliSISConeJetHeader.h"

#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.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;


ClassImp(AliSISConeJetFinder)


//____________________________________________________________________________

AliSISConeJetFinder::AliSISConeJetFinder():
AliJetFinder()
{
  // Constructor
}

//____________________________________________________________________________

AliSISConeJetFinder::~AliSISConeJetFinder()
{
  // destructor
}

//______________________________________________________________________________
void AliSISConeJetFinder::FindJets()
{

  // Pick up siscone header
  AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
  Bool_t debug  = header->GetDebug();     // debug option
  Int_t fOpt    = fReader->GetReaderHeader()->GetDetector();

  // Check if we are reading AOD jets
  TRefArray *refs = 0;
  Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
  if (fromAod) { refs = fReader->GetReferences(); }
  

//******************************** SISCONE PLUGIN CONFIGURATION
// Here we look for SISCone parameters in the header and we define our plugin.  

  Double_t coneRadius                = header->GetConeRadius();                 // cone radius
  Double_t overlapThreshold          = header->GetOverlapThreshold();           // overlap parameter
  Int_t    nPassMax                  = header->GetNPassMax();                   // maximum number of passes
  Double_t ptProtoJetMin             = header->GetPtProtojetMin();              // pT min of protojets
  Double_t caching                   = header->GetCaching();                    // do we record found cones for this set of data?

//  if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
//  Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets

  fastjet::JetDefinition::Plugin * plugin;
  plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);

  vector<fastjet::PseudoJet> inputParticles; 
//********************************  READING OF INPUT PARTICLES
// Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles. 

    if(fOpt==0) 
      {
	TClonesArray *lvArray = fReader->GetMomentumArray();
	Int_t nIn = lvArray->GetEntries();
	
	// We check if lvArray is ok
	if(lvArray == 0)
	  {
	    cout << "Could not get the momentum array" << endl;
	    return;
	  }
	
	if(nIn == 0)// nIn = Number of particles in the event
	  {
	    if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
	    return;
	  }
	
	fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects
	Float_t px,py,pz,en;
	
	// Load input vectors
	for(Int_t i = 0; i < nIn; i++)
	  { 
	    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
	    px = lv->Px();
	    py = lv->Py();
	    pz = lv->Pz();
	    en = lv->Energy();
	    
	    fastjet::PseudoJet inputPart(px,py,pz,en); 
	    inputPart.set_user_index(i);
	    inputParticles.push_back(inputPart); 
	    
	  }
      }
    else {
      TClonesArray* fUnit = fReader->GetUnitArray();
      if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
      Int_t         nCandidate = fReader->GetNumCandidate();
      Int_t         nIn = fUnit->GetEntries();
      if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
      fJets->SetNinput(nCandidate); // number of input objects // ME
      // Information extracted from fUnitArray
      // load input vectors and calculate total energy in array
      Float_t pt,eta,phi,theta,px,py,pz,en;
      Int_t ipart = 0;
      for(Int_t i=0; i<nIn; i++) 
	{
	  AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
	  
	  if(uArray->GetUnitEnergy()>0.){
	    
	    // It is not necessary anymore to cut on particle pt
	    pt    = uArray->GetUnitEnergy();
	    eta   = uArray->GetUnitEta();
	    phi   = uArray->GetUnitPhi();
	    theta = EtaToTheta(eta);
	    en    = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
	    px    = TMath::Cos(phi)*pt;
	    py    = TMath::Sin(phi)*pt;
	    pz    = en*TMath::TanH(eta);
	    if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
	    
	    fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
	    inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
	    inputParticles.push_back(inputPart);  // back of the input_particles vector 
	    ipart++;
	  }
	} // End loop on UnitArray 
    }
  
//******************************** CHOICE OF JET AREA  
// Here we determine jets area for subtracting background later
// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez

  Double_t ghostEtamax        = header->GetGhostEtaMax();       // maximum eta in which a ghost can be generated
  Double_t ghostArea          = header->GetGhostArea();         // area of a ghost
  Int_t    activeAreaRepeats  = header->GetActiveAreaRepeats(); // do we repeat area calculation?
  Double_t gridScatter        = header->GetGridScatter();       // fractional random fluctuations of the position of the ghosts on the y-phi grid
  Double_t ktScatter          = header->GetKtScatter();         // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
  Double_t meanGhostKt        = header->GetMeanGhostKt();       // average transverse momentum of the ghosts.

  Double_t areaTypeNumber  = header->GetAreaTypeNumber();       // the number determines jet area type 
  fastjet::AreaType areaType = fastjet::active_area;
  if (areaTypeNumber == 1)  areaType = fastjet::active_area;                 
  if (areaTypeNumber == 2)  areaType = fastjet::active_area_explicit_ghosts; 
  if (areaTypeNumber == 3)  areaType = fastjet::one_ghost_passive_area;      
  if (areaTypeNumber == 4)  areaType = fastjet::passive_area;                
  if (areaTypeNumber == 5)  areaType = fastjet::voronoi_area;                

  fastjet::AreaDefinition areaDef;

  if (areaTypeNumber < 5) 
  {
    fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
    areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
  }

  if (areaTypeNumber == 5)
  {
    Double_t effectiveRFact = header->GetEffectiveRFact();
    fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
    areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
  }

//********************************
//********************************

  Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not

//********************************
//********************************

//********************************
//******************************** BG SUBTRACTION
  if (bgMode == 1)// BG subtraction
  {
    //******************************** JETS FINDING AND EXTRACTION
    fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
    // Here we extract inclusive jets with pt > ptmin, sorted by pt 
    Double_t ptMin = header->GetMinJetPt(); 
    vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
    vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
  
    //***************************** BACKGROUND SUBTRACTION

    // Set the rapidity-azimuth range within which to study background
    Double_t rapMin = header->GetRapMin();
    Double_t rapMax = header->GetRapMax();
    Double_t phiMin = header->GetPhiMin();
    Double_t phiMax = header->GetPhiMax();
    fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);

    // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
    Int_t algo = header->GetBGAlgorithm();
    fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
    if (algo == 0) algorithm = fastjet::kt_algorithm;
    if (algo == 1) algorithm = fastjet::cambridge_algorithm;
    fastjet::JetDefinition jetDefForRho(algorithm, 0.5);
    fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
    Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
    cout<<"rho = "<<rho<<endl;

    // Vector of corrected jets
    vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
  
    //***************************** JETS DISPLAY

    for (size_t j = 0; j < jets.size(); j++)
    {
      // If the jet is only made of ghosts, continue.
      if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;

      // If the correction is > jet energy px = py = pz = e = 0
      if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;

      cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
      cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
//     cout<<"px = "<<jets[j].px()<<endl;
//     cout<<"py = "<<jets[j].py()<<endl;
//     cout<<"pz = "<<jets[j].pz()<<endl;
      cout<<"e = "<<jets[j].E()<<endl;

      cout<<"********************************** Corrected jet(s)"<<endl;
      cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
//     cout<<"px = "<<corrJets[j].px()<<endl;
//     cout<<"py = "<<corrJets[j].py()<<endl;
//     cout<<"pz = "<<corrJets[j].pz()<<endl;
      cout<<"e = "<<corrJets[j].E()<<endl;
    
      // Go to write AOD info
      AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
      if(debug) aodjet.Print("");
      AddJet(aodjet);
    }
  }

//********************************
//******************************** NO BG SUBTRACTION

  if (bgMode == 0)// No BG subtraction
  {
    //******************************** JETS FINDING AND EXTRACTION
    fastjet::ClusterSequence clustSeq(inputParticles, plugin);
    // Here we extract inclusive jets with pt > ptmin, sorted by pt 
    Double_t ptMin = header->GetMinJetPt(); 
    vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
    vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
  
    //***************************** JETS DISPLAY

    for (size_t k = 0; k < jets.size(); k++)
    {
      cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
      cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
//     cout<<"px = "<<jets[k].px()<<endl;
//     cout<<"py = "<<jets[k].py()<<endl;
//     cout<<"pz = "<<jets[k].pz()<<endl;
      cout<<"e = "<<jets[k].E()<<endl;

      // Go to write AOD info
      AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
      if(debug) aodjet.Print("");
      AddJet(aodjet);
    }
  }
}
 
//____________________________________________________________________________

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


}

//____________________________________________________________________________

void AliSISConeJetFinder::InitTask(TChain *tree)
{

  printf("SISCone jet finder initialization ******************");
  fReader->CreateTasks(tree);

}


//____________________________________________________________________________

void AliSISConeJetFinder::WriteJHeaderToFile()
{
  fHeader->Write();
}

//____________________________________________________________________________
Commit	Line	Data
2551af9d	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
	17	//---------------------------------------------------------------------
8838ab7a	18	// FastJet v2.3.4 finder algorithm interface
2551af9d	19	//
2551af9d	20	// Author: swensy.jangal@ires.in2p3.fr
8838ab7a	21	//
	22	// Last modification: Neutral cell energy included in the jet reconstruction
	23	//
	24	// Author: Magali.estienne@subatech.in2p3.fr
2551af9d	25	//---------------------------------------------------------------------
2551af9d	26
8838ab7a	27
2551af9d	28	#include <Riostream.h>
2551af9d	29	#include <TArrayF.h>
2551af9d	30	#include <TFile.h>
	31	#include <TH1F.h>
	32	#include <TH2F.h>
	33	#include <TLorentzVector.h>
	34	#include <TRandom.h>
8838ab7a	35	#include <TClonesArray.h>
2551af9d	36
	37	#include "AliHeader.h"
	38	#include "AliJet.h"
	39	#include "AliJetKineReader.h"
	40	#include "AliJetReader.h"
	41	#include "AliJetReaderHeader.h"
8838ab7a	42	#include "AliJetUnitArray.h"
2551af9d	43	#include "AliSISConeJetFinder.h"
	44	#include "AliSISConeJetHeader.h"
	45
	46	#include "fastjet/AreaDefinition.hh"
	47	#include "fastjet/ClusterSequenceArea.hh"
	48	#include "fastjet/JetDefinition.hh"
	49	#include "fastjet/PseudoJet.hh"
	50
8838ab7a	51	// get info on how fastjet was configured
2551af9d	52	#include "fastjet/config.h"
	53
	54	#ifdef ENABLE_PLUGIN_SISCONE
	55	#include "fastjet/SISConePlugin.hh"
	56	#endif
	57
8838ab7a	58	#include<sstream> // needed for internal io
2551af9d	59	#include<vector>
	60	#include <cmath>
	61
	62	using namespace std;
	63
8838ab7a	64
2551af9d	65	ClassImp(AliSISConeJetFinder)
2551af9d	66
8838ab7a	67
2551af9d	68	//____________________________________________________________________________
	69
	70	AliSISConeJetFinder::AliSISConeJetFinder():
	71	AliJetFinder()
	72	{
	73	// Constructor
	74	}
	75
	76	//____________________________________________________________________________
	77
	78	AliSISConeJetFinder::~AliSISConeJetFinder()
	79	{
	80	// destructor
	81	}
	82
	83	//______________________________________________________________________________
	84	void AliSISConeJetFinder::FindJets()
	85	{
	86
2551af9d	87	// Pick up siscone header
2551af9d	88	AliSISConeJetHeader header = (AliSISConeJetHeader)fHeader;
8838ab7a	89	Bool_t debug = header->GetDebug(); // debug option
8838ab7a	90	Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
2551af9d	91
	92	// Check if we are reading AOD jets
	93	TRefArray *refs = 0;
	94	Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
	95	if (fromAod) { refs = fReader->GetReferences(); }
	96
	97
	98	//******************************** SISCONE PLUGIN CONFIGURATION
	99	// Here we look for SISCone parameters in the header and we define our plugin.
	100
	101	Double_t coneRadius = header->GetConeRadius(); // cone radius
	102	Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter
	103	Int_t nPassMax = header->GetNPassMax(); // maximum number of passes
	104	Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets
	105	Double_t caching = header->GetCaching(); // do we record found cones for this set of data?
	106
8838ab7a	107	// if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
8838ab7a	108	// Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
2551af9d	109
	110	fastjet::JetDefinition::Plugin * plugin;
	111	plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
	112
8838ab7a	113	vector<fastjet::PseudoJet> inputParticles;
2551af9d	114	//******************************** READING OF INPUT PARTICLES
	115	// Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles.
	116
8838ab7a	117	if(fOpt==0)
	118	{
	119	TClonesArray *lvArray = fReader->GetMomentumArray();
	120	Int_t nIn = lvArray->GetEntries();
	121
	122	// We check if lvArray is ok
	123	if(lvArray == 0)
	124	{
	125	cout << "Could not get the momentum array" << endl;
	126	return;
	127	}
	128
	129	if(nIn == 0)// nIn = Number of particles in the event
	130	{
	131	if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
	132	return;
	133	}
	134
	135	fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects
	136	Float_t px,py,pz,en;
	137
	138	// Load input vectors
	139	for(Int_t i = 0; i < nIn; i++)
	140	{
	141	TLorentzVector lv = (TLorentzVector) lvArray->At(i);
	142	px = lv->Px();
	143	py = lv->Py();
	144	pz = lv->Pz();
	145	en = lv->Energy();
	146
	147	fastjet::PseudoJet inputPart(px,py,pz,en);
	148	inputPart.set_user_index(i);
	149	inputParticles.push_back(inputPart);
	150
	151	}
	152	}
	153	else {
	154	TClonesArray* fUnit = fReader->GetUnitArray();
	155	if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
	156	Int_t nCandidate = fReader->GetNumCandidate();
	157	Int_t nIn = fUnit->GetEntries();
	158	if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
	159	fJets->SetNinput(nCandidate); // number of input objects // ME
	160	// Information extracted from fUnitArray
	161	// load input vectors and calculate total energy in array
	162	Float_t pt,eta,phi,theta,px,py,pz,en;
	163	Int_t ipart = 0;
	164	for(Int_t i=0; i<nIn; i++)
	165	{
	166	AliJetUnitArray uArray = (AliJetUnitArray)fUnit->At(i);
	167
	168	if(uArray->GetUnitEnergy()>0.){
	169
	170	// It is not necessary anymore to cut on particle pt
	171	pt = uArray->GetUnitEnergy();
	172	eta = uArray->GetUnitEta();
	173	phi = uArray->GetUnitPhi();
	174	theta = EtaToTheta(eta);
	175	en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
	176	px = TMath::Cos(phi)*pt;
	177	py = TMath::Sin(phi)*pt;
	178	pz = en*TMath::TanH(eta);
	179	if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
	180
181	fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
182	inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
183	inputParticles.push_back(inputPart); // back of the input_particles vector
184	ipart++;
185	}
186	} // End loop on UnitArray
187	}
188
189	//******************************** CHOICE OF JET AREA
190	// Here we determine jets area for subtracting background later
191	// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
192
193	Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated
194	Double_t ghostArea = header->GetGhostArea(); // area of a ghost
195	Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation?
196	Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid
197	Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
198	Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts.
199
200	Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type
201	fastjet::AreaType areaType = fastjet::active_area;
202	if (areaTypeNumber == 1) areaType = fastjet::active_area;
203	if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts;
204	if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area;
205	if (areaTypeNumber == 4) areaType = fastjet::passive_area;
206	if (areaTypeNumber == 5) areaType = fastjet::voronoi_area;
207
208	fastjet::AreaDefinition areaDef;
209
210	if (areaTypeNumber < 5)
2551af9d	211	{
8838ab7a	212	fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
8838ab7a	213	areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d	214	}
2551af9d	215
8838ab7a	216	if (areaTypeNumber == 5)
2551af9d	217	{
8838ab7a	218	Double_t effectiveRFact = header->GetEffectiveRFact();
	219	fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
	220	areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d	221	}
2551af9d	222
8838ab7a	223	//********************************
8838ab7a	224	//********************************
2551af9d	225
8838ab7a	226	Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
2551af9d	227
8838ab7a	228	//********************************
8838ab7a	229	//********************************
2551af9d	230
8838ab7a	231	//********************************
	232	//******************************** BG SUBTRACTION
	233	if (bgMode == 1)// BG subtraction
	234	{
	235	//******************************** JETS FINDING AND EXTRACTION
	236	fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
	237	// Here we extract inclusive jets with pt > ptmin, sorted by pt
	238	Double_t ptMin = header->GetMinJetPt();
	239	vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
	240	vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
	241
	242	//***************************** BACKGROUND SUBTRACTION
	243
	244	// Set the rapidity-azimuth range within which to study background
	245	Double_t rapMin = header->GetRapMin();
	246	Double_t rapMax = header->GetRapMax();
	247	Double_t phiMin = header->GetPhiMin();
	248	Double_t phiMax = header->GetPhiMax();
	249	fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
	250
	251	// rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
	252	Int_t algo = header->GetBGAlgorithm();
	253	fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
	254	if (algo == 0) algorithm = fastjet::kt_algorithm;
	255	if (algo == 1) algorithm = fastjet::cambridge_algorithm;
	256	fastjet::JetDefinition jetDefForRho(algorithm, 0.5);
	257	fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
	258	Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
	259	cout<<"rho = "<<rho<<endl;
	260
	261	// Vector of corrected jets
	262	vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
	263
	264	//***************************** JETS DISPLAY
	265
	266	for (size_t j = 0; j < jets.size(); j++)
	267	{
	268	// If the jet is only made of ghosts, continue.
	269	if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;
	270
	271	// If the correction is > jet energy px = py = pz = e = 0
	272	if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
	273
	274	cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
	275	cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
	276	// cout<<"px = "<<jets[j].px()<<endl;
	277	// cout<<"py = "<<jets[j].py()<<endl;
	278	// cout<<"pz = "<<jets[j].pz()<<endl;
	279	cout<<"e = "<<jets[j].E()<<endl;
	280
	281	cout<<"********************************** Corrected jet(s)"<<endl;
	282	cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
	283	// cout<<"px = "<<corrJets[j].px()<<endl;
	284	// cout<<"py = "<<corrJets[j].py()<<endl;
	285	// cout<<"pz = "<<corrJets[j].pz()<<endl;
	286	cout<<"e = "<<corrJets[j].E()<<endl;
	287
	288	// Go to write AOD info
	289	AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
	290	if(debug) aodjet.Print("");
	291	AddJet(aodjet);
	292	}
	293	}
2551af9d	294
8838ab7a	295	//********************************
8838ab7a	296	//******************************** NO BG SUBTRACTION
2551af9d	297
8838ab7a	298	if (bgMode == 0)// No BG subtraction
2551af9d	299	{
8838ab7a	300	//******************************** JETS FINDING AND EXTRACTION
	301	fastjet::ClusterSequence clustSeq(inputParticles, plugin);
	302	// Here we extract inclusive jets with pt > ptmin, sorted by pt
	303	Double_t ptMin = header->GetMinJetPt();
	304	vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
	305	vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
	306
	307	//***************************** JETS DISPLAY
	308
	309	for (size_t k = 0; k < jets.size(); k++)
	310	{
	311	cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
	312	cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
	313	// cout<<"px = "<<jets[k].px()<<endl;
	314	// cout<<"py = "<<jets[k].py()<<endl;
	315	// cout<<"pz = "<<jets[k].pz()<<endl;
	316	cout<<"e = "<<jets[k].E()<<endl;
	317
	318	// Go to write AOD info
	319	AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
	320	if(debug) aodjet.Print("");
	321	AddJet(aodjet);
	322	}
2551af9d	323	}
	324	}
	325
8838ab7a	326	//____________________________________________________________________________
	327
	328	Float_t AliSISConeJetFinder::EtaToTheta(Float_t arg)
	329	{
	330	// return (180./TMath::Pi())2.atan(exp(-arg));
	331	return 2.*atan(exp(-arg));
	332
	333
	334	}
	335
	336	//____________________________________________________________________________
	337
	338	void AliSISConeJetFinder::InitTask(TChain *tree)
	339	{
	340
	341	printf("SISCone jet finder initialization ******************");
	342	fReader->CreateTasks(tree);
	343
	344	}
	345
	346
2551af9d	347	//____________________________________________________________________________
	348
	349	void AliSISConeJetFinder::WriteJHeaderToFile()
	350	{
	351	fHeader->Write();
	352	}
8838ab7a	353
	354	//____________________________________________________________________________
	355