[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 "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;
  AliJetReaderHeader *readerheader = (AliJetReaderHeader*)fReader->GetReaderHeader();

  Int_t debug = header->GetDebug();                        // debug option
  Int_t fOpt  = fReader->GetReaderHeader()->GetDetector();


//******************************** 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. 

  Double_t pttrackcut = readerheader->GetPtCut(); 
  if (debug)cout<<"pT track cut for SISCone jet finder = "<<pttrackcut<<" GeV"<<endl;

  if(fOpt==0) 
  {
    TClonesArray *lvArray    = fReader->GetMomentumArray();

    // We check if lvArray is ok
    if(lvArray == 0)
    {
      cout << "Could not get the momentum array" << endl;
      delete plugin;
      return;
    }
	
    Int_t nIn = lvArray->GetEntries();
	
    if(nIn == 0)// nIn = Number of particles in the event
    {
      if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
      delete plugin;
      return;
    }
	
    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();

      if (pttrackcut > TMath::Sqrt(px*px+py*py)) continue; //added by syssy
	    
      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; 	    delete plugin; return; }
    Int_t         nIn = fUnit->GetEntries();
    if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;  	    delete plugin; return; }
    // 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 (pttrackcut > TMath::Sqrt(px*px+py*py)) continue; 

	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);

    // As SISCone gives too small areas to estimate background,
    // we have to choose between kT and Cambridge/Aachen algorithms
    // to estimate mean background rho
    fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
    Int_t algo = header->GetBGAlgorithm();
    if (algo == 0) algorithm = fastjet::kt_algorithm;
    if (algo == 1) algorithm = fastjet::cambridge_algorithm;

    // Rho estimation :
    // Radius used to calculate rho, can be different from the one used to reconstruct jets
    Double_t RRho = header->GetRForRho(); 
    fastjet::JetDefinition jetDefinitionForRhoEstimation(algorithm, RRho);
    fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefinitionForRhoEstimation, areaDef);
    vector<fastjet::PseudoJet> inclusiveJetsForRhoEstimation = sorted_by_pt(csForRho.inclusive_jets());

    // Number of hard jets not to be used to estimate rho
    Int_t NHardJets = header->GetNumberOfJetsToErase();
    if (debug) cout<<"Number of jets before subtraction : "<<inclusiveJetsForRhoEstimation.size()<<endl;
    if (debug) cout<<"Number of jets not to count into rho estimation : "<<NHardJets<<endl;
    Bool_t rhoOnJet = 0;
    if (inclusiveJetsForRhoEstimation.size() == 1 && inclusiveJetsForRhoEstimation[0].perp() > 5) rhoOnJet = 1;//if there's only one jet (not bg kind), rho = 0
    if (inclusiveJetsForRhoEstimation.size() <= (UInt_t)NHardJets && inclusiveJetsForRhoEstimation.size() != 1  && inclusiveJetsForRhoEstimation.size() != 0)
    {
      //      inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(),inclusiveJetsForRhoEstimation.begin()+1);
      inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin());
      cout<<"Only the 1st (hardest) jet of the event hasn't been taken into account for rho estimation"<<endl;
    }
    if (inclusiveJetsForRhoEstimation.size() > (UInt_t)NHardJets) inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(), inclusiveJetsForRhoEstimation.begin()+NHardJets);
    
    // Estimation of rho and fluctuations sigma

    //1st method
    Double_t rho      = 0;
    Double_t sigma    = 0;
    Double_t meanarea = 0;
    // 3rd argument   : 1 = use area 4 vector rather than simple area
    // last argument  : 0 = in case of explicit ghosts use
    if (inclusiveJetsForRhoEstimation.size() != 0) csForRho.get_median_rho_and_sigma(inclusiveJetsForRhoEstimation, range, 1, rho, sigma, meanarea, 1); // this gives the fluctuations also

    if (rhoOnJet) rho = 0;

    if(debug) cout<<"rho = "<<rho<<endl;

    // Vector of corrected jets
    vector<fastjet::PseudoJet> corrJets = sorted_by_pt(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(corrJets[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;

      if(debug)
      {
	cout<<"********************************** Reconstructed jet(s) (not 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<<"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<<"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::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);

    //***************************** JETS DISPLAY

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

      // Go to write AOD info
      Double_t area = clustSeq.area(jets[k]);
      AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
      aodjet.SetEffArea(area,0);
      if(debug) aodjet.Print("");
      AddJet(aodjet);
    }
  }

  delete plugin;

}
 
//____________________________________________________________________________

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() const
{
  fHeader->Write();
}

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