[u/mrichter/AliRoot.git] / JETAN / jetan2004 / AliTkConeJet.cxx

//$Id$

#include <Riostream.h>
#include <TROOT.h>
#include <TClonesArray.h>
#include <TParticle.h>
#include <TMath.h>

#include "AliTkTowerV2.h"
#include "AliTkConeJet.h"

#define Ntowers__ 100

AliTkConeJet::AliTkConeJet() : TObject(), 
                         fEta(-999),fPhi(-999),fEt(-999),fType(0),
                         fNTowers(0),fTowers(0),fPtCut(0.),
			 fCEta(-999),fCPhi(-999),fCEt(-999),
			 fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
			 fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
{
  fTowers = new TClonesArray("AliTkTowerV2",Ntowers__);
}

AliTkConeJet::AliTkConeJet(Float_t pt, Float_t eta, Float_t phi,Int_t type) 
                        : TObject(), 
			  fEta(eta),fPhi(phi),fEt(pt),fType(type), 
			  fNTowers(0),fTowers(0),fPtCut(0.),
			  fCEta(eta),fCPhi(phi),fCEt(pt),
			  fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
			  fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
{
  fTowers = new TClonesArray("AliTkTowerV2",0);
}

AliTkConeJet::AliTkConeJet(const AliTkConeJet &j) 
                        : TObject(), 
			  fEta(-999), fPhi(-999), fEt(-999),fType(0), 
			  fNTowers(0),fTowers(0),fPtCut(0.),
			  fCEta(-999),fCPhi(-999),fCEt(-999),
			  fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
			  fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
{
  fEta = j.getEta();
  fPhi = j.getPhi();
  fEt = j.getEt();
  fCEta = j.getCEta();
  fCPhi = j.getCPhi();
  fCEt = j.getCEt();
  fNTowers = j.getNTowers();
  fPtCut = j.getPtCut();
  fCEta = j.getCEta();
  fCPhi = j.getCPhi();
  fCEt = j.getCEt();
  fPLength = j.getPLength();
  fXAxis = j.getXAxis();   
  fYAxis = j.getYAxis();
  fZAxis = j.getZAxis();
  fPtLead = j.getPtLead();
  fLeadPart = new TParticle(*j.getLeadPart());
  fNParts = j.getNParts();

  // need to build a copy of the old TClonesArray...
  fTowers = new TClonesArray("AliTkTowerV2",Ntowers__);
  TClonesArray *otherTowers = j.getTowerList();
  AliTkTowerV2 *myTower;
  TIterator *iter = otherTowers->MakeIterator();
  Int_t i =0;
  while ((myTower = (AliTkTowerV2 *) iter->Next()) != NULL) {
    new ((*fTowers)[i]) AliTkTowerV2(*myTower);
    i++;
  }
  if(i!=fNTowers)
    cerr << "AliTkConeJet: should not happen " << i << " " << fNTowers << endl;
  delete iter;

  // need to build a copy of the old TClonesArray...
  fParts = new TClonesArray("TParticle",j.getNParts());
  TClonesArray *otherParts = j.getParts();
  TParticle *myParticle;
  iter = otherParts->MakeIterator();
  i =0;
  while ((myParticle = (TParticle *) iter->Next()) != NULL) {
    new ((*fParts)[i]) TParticle(*myParticle);
    i++;
  }
  if(i!=fNParts)
    cerr << "AliTkConeJet: should not happen " << i << " " << fNParts << endl;
  delete iter;
}
 
AliTkConeJet::~AliTkConeJet() 
{
  delete fTowers;
  if(fLeadPart) delete fLeadPart;
  if(fParts) delete fParts;
}

void AliTkConeJet::addTower(AliTkTowerV2 *tower) 
{
  new ((*fTowers)[fNTowers]) AliTkTowerV2(*tower);
  fNTowers++; 
}

void AliTkConeJet::Clear(Option_t *option) 
{
  TObject::Clear(option);
  fTowers->Delete();
  fNTowers=0;
  if(fParts) delete fParts;
  fParts=0;
  if(fLeadPart) delete fLeadPart;
  fLeadPart=0;
  fEta=-999; 
  fPhi=-999;
  fEt=-999;
  fPtCut=0.;
  fCEta=-999; 
  fCPhi=-999;
  fCEt=-999;
  fPLength=0.;
  fXAxis=0.;
  fYAxis=0.;
  fZAxis=0.;
  fPtLead=0.;
  fLeadPart=0;
  fNParts=0;
  fParts=0;
}

void AliTkConeJet::calculateFromParticles(Float_t &et, Float_t &eta, Float_t &phi, Float_t ptcut)
{
  Float_t px=0.,py=0.,pz=0.;
  TIterator *titer = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)titer->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      px+=particle->Px();
      py+=particle->Py();
      pz+=particle->Pz();
    }
    delete part;
  }
  delete titer;

  et = TMath::Sqrt(px*px+py*py);
  Float_t p = TMath::Sqrt(px*px+py*py+pz*pz);
  Float_t theta = (pz==0)?TMath::PiOver2():TMath::ACos(pz/p);
  phi = TMath::Pi()+TMath::ATan2(-py,-px);
  eta = -TMath::Log(TMath::Tan(theta/2.));
}

/* or take e = p_t * cosh(eta) */
Float_t AliTkConeJet::getE(Float_t ptcut) const 
{
  Float_t e=0.;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      e+=particle->P();
    }
    delete part;
  }
  delete iter;
  return e;
}

/* or take e = p_t * cosh(eta) */
Float_t AliTkConeJet::getECharged(Float_t ptcut) const 
{
  Float_t e=0.;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getChargedParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      e+=particle->P();
    }
    particles->Delete();
    delete particles;
    delete part;
  }
  delete iter;
  return e;
}

Float_t AliTkConeJet::getEtCharged() const 
{
  Float_t EtCharged = 0;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    EtCharged += tower->getEtCharged();
  }
  delete iter;
  return EtCharged;
}

/* use to get fake rate if 
   in mixed case particles were marked */
Float_t AliTkConeJet::getEtChargedMarked(Float_t ptcut) const 
{
  Float_t et=0.;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getChargedParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->GetWeight()>-100) continue; //(-123)
      if(particle->Pt()<ptcut) continue;
      et+=particle->Pt();
    }
    particles->Delete();
    delete particles;
   delete part;
  }
  delete iter;
  return et;
}

/* use to get fake rate if 
   in mixed case particles were marked */
Float_t AliTkConeJet::getEtMarked(Float_t ptcut) const 
{
  Float_t et=0.;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->GetWeight()>-100) continue; //(-123)
      if(particle->Pt()<ptcut) continue;
      et+=particle->Pt();
    }
    delete part;
  }
  delete iter;
  return et;
}

Float_t AliTkConeJet::getEtMarkedFrac(Float_t ptcut) const 
{
  Float_t et=0.,etall=0.;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      etall+=particle->Pt();
      if(particle->GetWeight()>-100) continue; //(-123)
      et+=particle->Pt();
    }
    delete part;
  }
  delete iter;
  if(etall>0) return et/etall;
  return 1;
}

Float_t AliTkConeJet::getEtEM() const 
{
  Float_t EtEM = 0;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    EtEM += tower->getEtEM();
  }
  delete iter;
  return EtEM; 
}

TClonesArray *AliTkConeJet::getParticles(Float_t ptcut) const 
{
  TClonesArray *allParticles = new TClonesArray("TParticle",10000);
  Int_t nParticles = 0;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      new ((*allParticles)[nParticles]) TParticle(*particle);
      nParticles++;
    }
    delete part;
  }
  delete iter;
  return allParticles;
}

TClonesArray *AliTkConeJet::getChargedParticles(Float_t ptcut) const 
{
  TClonesArray *chargedParticles = new TClonesArray("TParticle",10000);
  Int_t nChargedParticles = 0;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getChargedParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      new ((*chargedParticles)[nChargedParticles]) TParticle(*particle);
      nChargedParticles++;
    }
    particles->Delete();
    delete particles;
    delete part;
  }
  delete iter;
  return chargedParticles;
}

TClonesArray *AliTkConeJet::getNeutralParticles(Float_t ptcut) const 
{
  TClonesArray *neutralParticles = new TClonesArray("TParticle",10000);
  Int_t nNeutralParticles = 0;
  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getNeutralParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      new ((*neutralParticles)[nNeutralParticles]) TParticle(*particle);
      nNeutralParticles++;
    }
    particles->Delete();
    delete particles;
    delete part;
  }
  delete iter;
  return neutralParticles;
}

Int_t AliTkConeJet::getNParticles() const 
{
  Int_t n = -1;
  TClonesArray *p = getParticles();
  n = p->GetEntries();
  p->Delete();
  delete p;
  return n;
}

Int_t AliTkConeJet::getNChargedParticles() const 
{
  Int_t n = -1;
  TClonesArray *p = getChargedParticles();
  n = p->GetEntries();
  p->Delete();
  delete p;
  return n;
}

Int_t AliTkConeJet::getNNeutralParticles() const 
{
  Int_t n = -1;
  TClonesArray *p = getNeutralParticles();
  n = p->GetEntries();
  p->Delete();
  delete p;
  return n;
}

Int_t AliTkConeJet::getNParticles(Float_t ptCut) const 
{
  Int_t n = 0;
  TClonesArray *particles = getParticles();
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  while ((particle = (TParticle *)iter->Next()) != NULL) {
    if (particle->Pt() > ptCut) {
      n++;
    }
  }
  delete iter;
  particles->Delete();
  delete particles;
  return n;
}

Int_t AliTkConeJet::getNChargedParticles(Float_t ptCut) const 
{
  Int_t n = 0;
  TClonesArray *particles = getChargedParticles();
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  while ((particle = (TParticle *)iter->Next()) != NULL) {
    if (particle->Pt() > ptCut) {
      n++;
    }
  }
  delete iter;
  particles->Delete();
  delete particles;
  return n;
}

Int_t AliTkConeJet::getNNeutralParticles(Float_t ptCut) const 
{
  Int_t n = 0;
  TClonesArray *particles = getNeutralParticles();
  TIterator *iter = particles->MakeIterator();
  TParticle *particle;
  while ((particle = (TParticle *)iter->Next()) != NULL) {
    if (particle->Pt() > ptCut) {
      n++;
    }
  }
  delete iter;
  particles->Delete();
  delete particles;
  return n;
}

void AliTkConeJet::getAxis(Float_t &x,Float_t &y,Float_t &z,Float_t ptcut) const
{
  x=0.;
  y=0.;
  z=0.;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      x+=particle->Px();
      y+=particle->Py();
      z+=particle->Pz();
    }
    delete part;
  }
  delete iter;
}

void AliTkConeJet::getChAxis(Float_t &x,Float_t &y,Float_t &z,Float_t ptcut) const
{
  x=0.;
  y=0.;
  z=0.;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getChargedParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      x+=particle->Px();
      y+=particle->Py();
      z+=particle->Pz();
    }
    particles->Delete();
    delete particles;
    delete part;
  }
  delete iter;
  Float_t length=TMath::Sqrt(x*x+y*y+z*z);
  if(length!=0){
    x/=length;
    y/=length;
    z/=length;
  }
}

TParticle* AliTkConeJet::getLeadingPart(Float_t ptcut) const
{
  TParticle* plead=new TParticle; 
  Float_t ptlead=0;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      if(particle->Pt()>ptlead){
	ptlead=particle->Pt();
	delete plead;
        plead=new TParticle(*particle);
      }
    }
    delete part;
  }
  delete iter;
  return plead;
}

TParticle* AliTkConeJet::getLeadingChPart(Float_t ptcut) const
{
  TParticle* plead=new TParticle; 
  Float_t ptlead=0;

  TIterator *iter = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
    TClonesArray *particles = tower->getChargedParticleList();
    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      if(particle->Pt()>ptlead){
	ptlead=particle->Pt();
	delete plead;
        plead=new TParticle(*particle);
      }
    }
    particles->Delete();
    delete particles;
    delete part;
  }
  delete iter;
  return plead;
}

void AliTkConeJet::calculateValues(Float_t ptcut)
{
  fPtCut=ptcut;
  if(fParts) delete fParts;
  fParts=new TClonesArray("TParticle",fNTowers*25);

  Float_t px=0.,py=0.,pz=0.;
  TParticle* plead=new TParticle; 
  Float_t ptlead=0;
  Int_t counter=0;
  TIterator *titer = fTowers->MakeIterator();
  AliTkTowerV2 *tower = NULL;
  while ((tower = (AliTkTowerV2 *)titer->Next()) != NULL) {
    TClonesArray *particles = tower->getParticleList();

    TParticle *particle = NULL;
    TIterator *part = particles->MakeIterator();
    while ((particle = (TParticle *)part->Next()) != NULL) {
      if(particle->Pt()<ptcut) continue;
      px+=particle->Px();
      py+=particle->Py();
      pz+=particle->Pz();
      new ((*fParts)[counter]) TParticle(*particle);
      counter++;
      if(particle->Pt()>ptlead){
	ptlead=particle->Pt();
	delete plead;
        plead=new TParticle(*particle);
      }
    }
    delete part;
  }
  delete titer;

  fCEt = TMath::Sqrt(px*px+py*py);
  Float_t p = TMath::Sqrt(px*px+py*py+pz*pz);
  fPLength = p;
  //Float_t theta = ((pz==0)||(p==0))?TMath::PiOver2():TMath::ACos(pz/p);
  fCPhi = TMath::Pi()+TMath::ATan2(-py,-px);
  //fCEta = (theta==0)?-100:TMath::Log(TMath::Tan(theta/2.));
  fCEta = (p==pz)?-100:0.5*TMath::Log((p+pz)/(p-pz));
  if(p){
    fXAxis= px/p;
    fYAxis= py/p;
    fZAxis= pz/p;
  }
  fPtLead=ptlead;
  if(fLeadPart) delete fLeadPart;
  fLeadPart=new TParticle(*plead);
  delete plead;
  fNParts=counter;
}

Int_t AliTkConeJet::Compare(const TObject *obj) const
{
  Double_t val=((AliTkConeJet*)obj)->getEt();

  if(fEt>val) return -1; //qsort is ascending
  else if (fEt<val) return 1;
  else return 0;
}

ostream& operator<<(ostream& s,const AliTkConeJet &p) 
{
  return s << "AliTkConeJet info: eta=" << p.getEta()
	   << " phi=" << p.getPhi() << " Et=" << p.getEt();
}

Float_t AliTkConeJet::Diff(const AliTkConeJet *jet1, const AliTkConeJet *jet2, Float_t &etdiff, Float_t &phidiff, Float_t &etadiff)
{
  Float_t ret=0;
  phidiff=TMath::Abs(jet1->getPhi()-jet2->getPhi());
  if(phidiff>TMath::Pi()) ret=TMath::TwoPi()-phidiff;
  etadiff=TMath::Abs(jet1->getEta()-jet2->getEta());
  ret=TMath::Sqrt(phidiff*phidiff+etadiff*etadiff);
  etdiff=jet1->getEt()-jet2->getEt();
  return ret;
}

ClassImp(AliTkConeJet)
Commit	Line	Data
b9a6a391	1	//$Id$
	2
	3	#include <Riostream.h>
	4	#include <TROOT.h>
	5	#include <TClonesArray.h>
	6	#include <TParticle.h>
	7	#include <TMath.h>
	8
	9	#include "AliTkTowerV2.h"
	10	#include "AliTkConeJet.h"
	11
	12	#define Ntowers__ 100
	13
	14	AliTkConeJet::AliTkConeJet() : TObject(),
	15	fEta(-999),fPhi(-999),fEt(-999),fType(0),
	16	fNTowers(0),fTowers(0),fPtCut(0.),
	17	fCEta(-999),fCPhi(-999),fCEt(-999),
	18	fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
	19	fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
	20	{
	21	fTowers = new TClonesArray("AliTkTowerV2",Ntowers__);
	22	}
	23
	24	AliTkConeJet::AliTkConeJet(Float_t pt, Float_t eta, Float_t phi,Int_t type)
	25	: TObject(),
	26	fEta(eta),fPhi(phi),fEt(pt),fType(type),
	27	fNTowers(0),fTowers(0),fPtCut(0.),
	28	fCEta(eta),fCPhi(phi),fCEt(pt),
	29	fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
	30	fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
	31	{
	32	fTowers = new TClonesArray("AliTkTowerV2",0);
	33	}
	34
	35	AliTkConeJet::AliTkConeJet(const AliTkConeJet &j)
	36	: TObject(),
	37	fEta(-999), fPhi(-999), fEt(-999),fType(0),
	38	fNTowers(0),fTowers(0),fPtCut(0.),
	39	fCEta(-999),fCPhi(-999),fCEt(-999),
	40	fPLength(0.),fXAxis(0.),fYAxis(0.),fZAxis(0.),
	41	fPtLead(0.),fLeadPart(0),fNParts(0),fParts(0)
	42	{
	43	fEta = j.getEta();
	44	fPhi = j.getPhi();
	45	fEt = j.getEt();
	46	fCEta = j.getCEta();
	47	fCPhi = j.getCPhi();
	48	fCEt = j.getCEt();
	49	fNTowers = j.getNTowers();
	50	fPtCut = j.getPtCut();
	51	fCEta = j.getCEta();
	52	fCPhi = j.getCPhi();
	53	fCEt = j.getCEt();
	54	fPLength = j.getPLength();
	55	fXAxis = j.getXAxis();
	56	fYAxis = j.getYAxis();
	57	fZAxis = j.getZAxis();
	58	fPtLead = j.getPtLead();
	59	fLeadPart = new TParticle(*j.getLeadPart());
	60	fNParts = j.getNParts();
	61
	62	// need to build a copy of the old TClonesArray...
	63	fTowers = new TClonesArray("AliTkTowerV2",Ntowers__);
	64	TClonesArray *otherTowers = j.getTowerList();
65	AliTkTowerV2 *myTower;
66	TIterator *iter = otherTowers->MakeIterator();
67	Int_t i =0;
68	while ((myTower = (AliTkTowerV2 *) iter->Next()) != NULL) {
69	new ((fTowers)[i]) AliTkTowerV2(myTower);
70	i++;
71	}
72	if(i!=fNTowers)
73	cerr << "AliTkConeJet: should not happen " << i << " " << fNTowers << endl;
74	delete iter;
75
76	// need to build a copy of the old TClonesArray...
77	fParts = new TClonesArray("TParticle",j.getNParts());
78	TClonesArray *otherParts = j.getParts();
79	TParticle *myParticle;
80	iter = otherParts->MakeIterator();
81	i =0;
82	while ((myParticle = (TParticle *) iter->Next()) != NULL) {
83	new ((fParts)[i]) TParticle(myParticle);
84	i++;
85	}
86	if(i!=fNParts)
87	cerr << "AliTkConeJet: should not happen " << i << " " << fNParts << endl;
88	delete iter;
89	}
90
91	AliTkConeJet::~AliTkConeJet()
92	{
93	delete fTowers;
94	if(fLeadPart) delete fLeadPart;
95	if(fParts) delete fParts;
96	}
97
98	void AliTkConeJet::addTower(AliTkTowerV2 *tower)
99	{
100	new ((fTowers)[fNTowers]) AliTkTowerV2(tower);
101	fNTowers++;
102	}
103
104	void AliTkConeJet::Clear(Option_t *option)
105	{
106	TObject::Clear(option);
107	fTowers->Delete();
108	fNTowers=0;
109	if(fParts) delete fParts;
110	fParts=0;
111	if(fLeadPart) delete fLeadPart;
112	fLeadPart=0;
113	fEta=-999;
114	fPhi=-999;
115	fEt=-999;
116	fPtCut=0.;
117	fCEta=-999;
118	fCPhi=-999;
119	fCEt=-999;
120	fPLength=0.;
121	fXAxis=0.;
122	fYAxis=0.;
123	fZAxis=0.;
124	fPtLead=0.;
125	fLeadPart=0;
126	fNParts=0;
127	fParts=0;
128	}
129
130	void AliTkConeJet::calculateFromParticles(Float_t &et, Float_t &eta, Float_t &phi, Float_t ptcut)
131	{
132	Float_t px=0.,py=0.,pz=0.;
133	TIterator *titer = fTowers->MakeIterator();
134	AliTkTowerV2 *tower = NULL;
135	while ((tower = (AliTkTowerV2 *)titer->Next()) != NULL) {
136	TClonesArray *particles = tower->getParticleList();
137	TParticle *particle = NULL;
138	TIterator *part = particles->MakeIterator();
139	while ((particle = (TParticle *)part->Next()) != NULL) {
140	if(particle->Pt()<ptcut) continue;
141	px+=particle->Px();
142	py+=particle->Py();
143	pz+=particle->Pz();
144	}
145	delete part;
146	}
147	delete titer;
148
149	et = TMath::Sqrt(pxpx+pypy);
150	Float_t p = TMath::Sqrt(pxpx+pypy+pz*pz);
151	Float_t theta = (pz==0)?TMath::PiOver2():TMath::ACos(pz/p);
152	phi = TMath::Pi()+TMath::ATan2(-py,-px);
153	eta = -TMath::Log(TMath::Tan(theta/2.));
154	}
155
156	/* or take e = p_t * cosh(eta) */
157	Float_t AliTkConeJet::getE(Float_t ptcut) const
158	{
159	Float_t e=0.;
160
161	TIterator *iter = fTowers->MakeIterator();
162	AliTkTowerV2 *tower = NULL;
163	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
164	TClonesArray *particles = tower->getParticleList();
165	TParticle *particle = NULL;
166	TIterator *part = particles->MakeIterator();
167	while ((particle = (TParticle *)part->Next()) != NULL) {
168	if(particle->Pt()<ptcut) continue;
169	e+=particle->P();
170	}
171	delete part;
172	}
173	delete iter;
174	return e;
175	}
176
177	/* or take e = p_t * cosh(eta) */
178	Float_t AliTkConeJet::getECharged(Float_t ptcut) const
179	{
180	Float_t e=0.;
181
182	TIterator *iter = fTowers->MakeIterator();
183	AliTkTowerV2 *tower = NULL;
184	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
185	TClonesArray *particles = tower->getChargedParticleList();
186	TParticle *particle = NULL;
187	TIterator *part = particles->MakeIterator();
188	while ((particle = (TParticle *)part->Next()) != NULL) {
189	if(particle->Pt()<ptcut) continue;
190	e+=particle->P();
191	}
192	particles->Delete();
193	delete particles;
194	delete part;
195	}
196	delete iter;
197	return e;
198	}
199
200	Float_t AliTkConeJet::getEtCharged() const
201	{
202	Float_t EtCharged = 0;
203	TIterator *iter = fTowers->MakeIterator();
204	AliTkTowerV2 *tower = NULL;
205	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
206	EtCharged += tower->getEtCharged();
207	}
208	delete iter;
209	return EtCharged;
210	}
211
212	/* use to get fake rate if
213	in mixed case particles were marked */
214	Float_t AliTkConeJet::getEtChargedMarked(Float_t ptcut) const
215	{
216	Float_t et=0.;
217
218	TIterator *iter = fTowers->MakeIterator();
219	AliTkTowerV2 *tower = NULL;
220	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
221	TClonesArray *particles = tower->getChargedParticleList();
222	TParticle *particle = NULL;
223	TIterator *part = particles->MakeIterator();
224	while ((particle = (TParticle *)part->Next()) != NULL) {
225	if(particle->GetWeight()>-100) continue; //(-123)
226	if(particle->Pt()<ptcut) continue;
227	et+=particle->Pt();
228	}
229	particles->Delete();
230	delete particles;
231	delete part;
232	}
233	delete iter;
234	return et;
235	}
236
237	/* use to get fake rate if
238	in mixed case particles were marked */
239	Float_t AliTkConeJet::getEtMarked(Float_t ptcut) const
240	{
241	Float_t et=0.;
242	TIterator *iter = fTowers->MakeIterator();
243	AliTkTowerV2 *tower = NULL;
244	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
245	TClonesArray *particles = tower->getParticleList();
246	TParticle *particle = NULL;
247	TIterator *part = particles->MakeIterator();
248	while ((particle = (TParticle *)part->Next()) != NULL) {
249	if(particle->GetWeight()>-100) continue; //(-123)
250	if(particle->Pt()<ptcut) continue;
251	et+=particle->Pt();
252	}
253	delete part;
254	}
255	delete iter;
256	return et;
257	}
258
259	Float_t AliTkConeJet::getEtMarkedFrac(Float_t ptcut) const
260	{
261	Float_t et=0.,etall=0.;
262	TIterator *iter = fTowers->MakeIterator();
263	AliTkTowerV2 *tower = NULL;
264	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
265	TClonesArray *particles = tower->getParticleList();
266	TParticle *particle = NULL;
267	TIterator *part = particles->MakeIterator();
268	while ((particle = (TParticle *)part->Next()) != NULL) {
269	if(particle->Pt()<ptcut) continue;
270	etall+=particle->Pt();
271	if(particle->GetWeight()>-100) continue; //(-123)
272	et+=particle->Pt();
273	}
274	delete part;
275	}
276	delete iter;
277	if(etall>0) return et/etall;
278	return 1;
279	}
280
281	Float_t AliTkConeJet::getEtEM() const
282	{
283	Float_t EtEM = 0;
284	TIterator *iter = fTowers->MakeIterator();
285	AliTkTowerV2 *tower = NULL;
286	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
287	EtEM += tower->getEtEM();
288	}
289	delete iter;
290	return EtEM;
291	}
292
293	TClonesArray *AliTkConeJet::getParticles(Float_t ptcut) const
294	{
295	TClonesArray *allParticles = new TClonesArray("TParticle",10000);
296	Int_t nParticles = 0;
297	TIterator *iter = fTowers->MakeIterator();
298	AliTkTowerV2 *tower = NULL;
299	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
300	TClonesArray *particles = tower->getParticleList();
301	TParticle *particle = NULL;
302	TIterator *part = particles->MakeIterator();
303	while ((particle = (TParticle *)part->Next()) != NULL) {
304	if(particle->Pt()<ptcut) continue;
305	new ((allParticles)[nParticles]) TParticle(particle);
306	nParticles++;
307	}
308	delete part;
309	}
310	delete iter;
311	return allParticles;
312	}
313
314	TClonesArray *AliTkConeJet::getChargedParticles(Float_t ptcut) const
315	{
316	TClonesArray *chargedParticles = new TClonesArray("TParticle",10000);
317	Int_t nChargedParticles = 0;
318	TIterator *iter = fTowers->MakeIterator();
319	AliTkTowerV2 *tower = NULL;
320	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
321	TClonesArray *particles = tower->getChargedParticleList();
322	TParticle *particle = NULL;
323	TIterator *part = particles->MakeIterator();
324	while ((particle = (TParticle *)part->Next()) != NULL) {
325	if(particle->Pt()<ptcut) continue;
326	new ((chargedParticles)[nChargedParticles]) TParticle(particle);
327	nChargedParticles++;
328	}
329	particles->Delete();
330	delete particles;
331	delete part;
332	}
333	delete iter;
334	return chargedParticles;
335	}
336
337	TClonesArray *AliTkConeJet::getNeutralParticles(Float_t ptcut) const
338	{
339	TClonesArray *neutralParticles = new TClonesArray("TParticle",10000);
340	Int_t nNeutralParticles = 0;
341	TIterator *iter = fTowers->MakeIterator();
342	AliTkTowerV2 *tower = NULL;
343	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
344	TClonesArray *particles = tower->getNeutralParticleList();
345	TParticle *particle = NULL;
346	TIterator *part = particles->MakeIterator();
347	while ((particle = (TParticle *)part->Next()) != NULL) {
348	if(particle->Pt()<ptcut) continue;
349	new ((neutralParticles)[nNeutralParticles]) TParticle(particle);
350	nNeutralParticles++;
351	}
352	particles->Delete();
353	delete particles;
354	delete part;
355	}
356	delete iter;
357	return neutralParticles;
358	}
359
360	Int_t AliTkConeJet::getNParticles() const
361	{
362	Int_t n = -1;
363	TClonesArray *p = getParticles();
364	n = p->GetEntries();
365	p->Delete();
366	delete p;
367	return n;
368	}
369
370	Int_t AliTkConeJet::getNChargedParticles() const
371	{
372	Int_t n = -1;
373	TClonesArray *p = getChargedParticles();
374	n = p->GetEntries();
375	p->Delete();
376	delete p;
377	return n;
378	}
379
380	Int_t AliTkConeJet::getNNeutralParticles() const
381	{
382	Int_t n = -1;
383	TClonesArray *p = getNeutralParticles();
384	n = p->GetEntries();
385	p->Delete();
386	delete p;
387	return n;
388	}
389
390	Int_t AliTkConeJet::getNParticles(Float_t ptCut) const
391	{
392	Int_t n = 0;
393	TClonesArray *particles = getParticles();
394	TIterator *iter = particles->MakeIterator();
395	TParticle *particle;
396	while ((particle = (TParticle *)iter->Next()) != NULL) {
397	if (particle->Pt() > ptCut) {
398	n++;
399	}
400	}
401	delete iter;
402	particles->Delete();
403	delete particles;
404	return n;
405	}
406
407	Int_t AliTkConeJet::getNChargedParticles(Float_t ptCut) const
408	{
409	Int_t n = 0;
410	TClonesArray *particles = getChargedParticles();
411	TIterator *iter = particles->MakeIterator();
412	TParticle *particle;
413	while ((particle = (TParticle *)iter->Next()) != NULL) {
414	if (particle->Pt() > ptCut) {
415	n++;
416	}
417	}
418	delete iter;
419	particles->Delete();
420	delete particles;
421	return n;
422	}
423
424	Int_t AliTkConeJet::getNNeutralParticles(Float_t ptCut) const
425	{
426	Int_t n = 0;
427	TClonesArray *particles = getNeutralParticles();
428	TIterator *iter = particles->MakeIterator();
429	TParticle *particle;
430	while ((particle = (TParticle *)iter->Next()) != NULL) {
431	if (particle->Pt() > ptCut) {
432	n++;
433	}
434	}
435	delete iter;
436	particles->Delete();
437	delete particles;
438	return n;
439	}
440
441	void AliTkConeJet::getAxis(Float_t &x,Float_t &y,Float_t &z,Float_t ptcut) const
442	{
443	x=0.;
444	y=0.;
445	z=0.;
446
447	TIterator *iter = fTowers->MakeIterator();
448	AliTkTowerV2 *tower = NULL;
449	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
450	TClonesArray *particles = tower->getParticleList();
451	TParticle *particle = NULL;
452	TIterator *part = particles->MakeIterator();
453	while ((particle = (TParticle *)part->Next()) != NULL) {
454	if(particle->Pt()<ptcut) continue;
455	x+=particle->Px();
456	y+=particle->Py();
457	z+=particle->Pz();
458	}
459	delete part;
460	}
461	delete iter;
462	}
463
464	void AliTkConeJet::getChAxis(Float_t &x,Float_t &y,Float_t &z,Float_t ptcut) const
465	{
466	x=0.;
467	y=0.;
468	z=0.;
469
470	TIterator *iter = fTowers->MakeIterator();
471	AliTkTowerV2 *tower = NULL;
472	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
473	TClonesArray *particles = tower->getChargedParticleList();
474	TParticle *particle = NULL;
475	TIterator *part = particles->MakeIterator();
476	while ((particle = (TParticle *)part->Next()) != NULL) {
477	if(particle->Pt()<ptcut) continue;
478	x+=particle->Px();
479	y+=particle->Py();
480	z+=particle->Pz();
481	}
482	particles->Delete();
483	delete particles;
484	delete part;
485	}
486	delete iter;
487	Float_t length=TMath::Sqrt(xx+yy+z*z);
488	if(length!=0){
489	x/=length;
490	y/=length;
491	z/=length;
492	}
493	}
494
495	TParticle* AliTkConeJet::getLeadingPart(Float_t ptcut) const
496	{
497	TParticle* plead=new TParticle;
498	Float_t ptlead=0;
499
500	TIterator *iter = fTowers->MakeIterator();
501	AliTkTowerV2 *tower = NULL;
502	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
503	TClonesArray *particles = tower->getParticleList();
504	TParticle *particle = NULL;
505	TIterator *part = particles->MakeIterator();
506	while ((particle = (TParticle *)part->Next()) != NULL) {
507	if(particle->Pt()<ptcut) continue;
508	if(particle->Pt()>ptlead){
509	ptlead=particle->Pt();
510	delete plead;
511	plead=new TParticle(*particle);
512	}
513	}
514	delete part;
515	}
516	delete iter;
517	return plead;
518	}
519
520	TParticle* AliTkConeJet::getLeadingChPart(Float_t ptcut) const
521	{
522	TParticle* plead=new TParticle;
523	Float_t ptlead=0;
524
525	TIterator *iter = fTowers->MakeIterator();
526	AliTkTowerV2 *tower = NULL;
527	while ((tower = (AliTkTowerV2 *)iter->Next()) != NULL) {
528	TClonesArray *particles = tower->getChargedParticleList();
529	TParticle *particle = NULL;
530	TIterator *part = particles->MakeIterator();
531	while ((particle = (TParticle *)part->Next()) != NULL) {
532	if(particle->Pt()<ptcut) continue;
533	if(particle->Pt()>ptlead){
534	ptlead=particle->Pt();
535	delete plead;
536	plead=new TParticle(*particle);
537	}
538	}
539	particles->Delete();
540	delete particles;
541	delete part;
542	}
543	delete iter;
544	return plead;
545	}
546
547	void AliTkConeJet::calculateValues(Float_t ptcut)
548	{
549	fPtCut=ptcut;
550	if(fParts) delete fParts;
551	fParts=new TClonesArray("TParticle",fNTowers*25);
552
553	Float_t px=0.,py=0.,pz=0.;
554	TParticle* plead=new TParticle;
555	Float_t ptlead=0;
556	Int_t counter=0;
557	TIterator *titer = fTowers->MakeIterator();
558	AliTkTowerV2 *tower = NULL;
559	while ((tower = (AliTkTowerV2 *)titer->Next()) != NULL) {
560	TClonesArray *particles = tower->getParticleList();
561
562	TParticle *particle = NULL;
563	TIterator *part = particles->MakeIterator();
564	while ((particle = (TParticle *)part->Next()) != NULL) {
565	if(particle->Pt()<ptcut) continue;
566	px+=particle->Px();
567	py+=particle->Py();
568	pz+=particle->Pz();
569	new ((fParts)[counter]) TParticle(particle);
570	counter++;
571	if(particle->Pt()>ptlead){
572	ptlead=particle->Pt();
573	delete plead;
574	plead=new TParticle(*particle);
575	}
576	}
577	delete part;
578	}
579	delete titer;
580
581	fCEt = TMath::Sqrt(pxpx+pypy);
582	Float_t p = TMath::Sqrt(pxpx+pypy+pz*pz);
583	fPLength = p;
584	//Float_t theta = ((pz==0)\|\|(p==0))?TMath::PiOver2():TMath::ACos(pz/p);
585	fCPhi = TMath::Pi()+TMath::ATan2(-py,-px);
586	//fCEta = (theta==0)?-100:TMath::Log(TMath::Tan(theta/2.));
587	fCEta = (p==pz)?-100:0.5*TMath::Log((p+pz)/(p-pz));
588	if(p){
589	fXAxis= px/p;
590	fYAxis= py/p;
591	fZAxis= pz/p;
592	}
593	fPtLead=ptlead;
594	if(fLeadPart) delete fLeadPart;
595	fLeadPart=new TParticle(*plead);
596	delete plead;
597	fNParts=counter;
598	}
599
600	Int_t AliTkConeJet::Compare(const TObject *obj) const
601	{
602	Double_t val=((AliTkConeJet*)obj)->getEt();
603
604	if(fEt>val) return -1; //qsort is ascending
605	else if (fEt<val) return 1;
606	else return 0;
607	}
608
609	ostream& operator<<(ostream& s,const AliTkConeJet &p)
610	{
611	return s << "AliTkConeJet info: eta=" << p.getEta()
612	<< " phi=" << p.getPhi() << " Et=" << p.getEt();
613	}
614
615	Float_t AliTkConeJet::Diff(const AliTkConeJet jet1, const AliTkConeJet jet2, Float_t &etdiff, Float_t &phidiff, Float_t &etadiff)
616	{
617	Float_t ret=0;
618	phidiff=TMath::Abs(jet1->getPhi()-jet2->getPhi());
619	if(phidiff>TMath::Pi()) ret=TMath::TwoPi()-phidiff;
620	etadiff=TMath::Abs(jet1->getEta()-jet2->getEta());
621	ret=TMath::Sqrt(phidiffphidiff+etadiffetadiff);
622	etdiff=jet1->getEt()-jet2->getEt();
623	return ret;
624	}
625
626	ClassImp(AliTkConeJet)