[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetsTM.cxx

#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TList.h"
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"

#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
#include "AliAnalysisTaskJetsTM.h"

ClassImp(AliAnalysisTaskJetsTM)
//
// Thrust Major (TM) analysis of reconstructed jets.
// TM is the thrust in the plane perpendicular to the jet axis
// The present amalysis performs the following steps:
// (a) Construct to orthogonal unit vectors (e1, e2) in the plane perpendicular to the jet axis
// (b) Calculate the components of all particles with jT > 1 GeV with respect to e1, e2
// (c) Construct the sphericity matrix
// (d) Find the two orthogonal eigenvectors of the spericity matrix
// (e) Caluclate the components of all particles with jT < 1 GeV in the reference frame spanned by the eigenvectors
// (f) Calculate the azimuthal angle in this frame
//
//
// Author: andreas.morsch@cern.ch


AliAnalysisTaskJetsTM::AliAnalysisTaskJetsTM(const char *name) 
    : AliAnalysisTaskSE(name) 
      ,fHists(0)
      ,fPtH(0)
      ,fPtTH(0)
      ,fPhiM(0)
      ,fPhiMPt(0)
      ,fPhiMPtJ(0)
      ,fPtSum(0)
{
  //
  // Constructor
  //
  DefineOutput(1,  TList::Class());
}


//________________________________________________________________________
void AliAnalysisTaskJetsTM::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
    fHists    = new TList();
    fPtH      = new TH1F("fPtH" ,     " p_{T}       distribution of jets"  ,            50, 0., 100.0);
    fPtSum    = new TH2F("fPtSum" ,   " p_{T} sum in cone"                 ,            50, 0., 100.0, 50, 0, 100.);
    fPtTH     = new TH1F("fPtTH",     " p_{T}       distribution of tracks",            50, 0., 100.0);
    fPhiM     = new TH1F("fPhiM",     " phi^{major} distribution of tracks",            64, 0.,   3.2);
    fPhiMPt   = new TH2F("fPhiMPt",   " phi^{major} distribution of tracks vs pt",      64, 0.,   3.2, 20, 0., 2.);
    fPhiMPtJ  = new TH2F("fPhiMPtJ",  " phi^{major} distribution of tracks vs pt Jet",  64, 0.,   3.2, 40, 0., 100.);

    fHists->SetOwner();
    fHists->Add(fPtH);
    fHists->Add(fPtSum);
    fHists->Add(fPtTH);
    fHists->Add(fPhiM);
    fHists->Add(fPhiMPt);
    fHists->Add(fPhiMPtJ);
}


//________________________________________________________________________
void AliAnalysisTaskJetsTM::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event

  if (!fInputEvent) {
    Printf("ERROR: AOD not available");
    return;
  }

  AliAODEvent* aodE  = dynamic_cast<AliAODEvent*>  (fInputEvent);
  TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
  Int_t nJ = jets->GetEntries();

  Float_t ptmax = 0.;
  Int_t   imax  = -1;
  
//
// Find highest pT jet with pt > 20 GeV
//
  for (Int_t i = 0; i < nJ; i++) {
      AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
      Float_t ptJ  = jet->Pt();
      Float_t etaJ = TMath::Abs(jet->Eta());
      if ((ptJ > 20.) && (ptJ  > ptmax) && etaJ < 0.5) {
	  ptmax = ptJ;
	  imax = i;
      }
  }

  if (imax == -1) return;
  

  AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
  Float_t ptJ  = jet->Pt();
  Float_t phiJ = jet->Phi();
  Float_t etaJ = jet->Eta();
  
  fPtH->Fill(ptJ);

//
// The transverse plane
//
// 2 normalized vectors perpendicular to the jet direction
//
      Float_t pxJ = jet->Px();
      Float_t pyJ = jet->Py();
      Float_t pzJ = jet->Pz();
      
      TVector3  ppJ1(pxJ, pyJ, pzJ);
      TVector3  ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
      ppJ3.SetMag(1.);
      TVector3  ppJ2(-pyJ, pxJ, 0);
      ppJ2.SetMag(1.);

//
// 1st track loop to determine the sphericity matrix
//   
      Int_t nT = aodE->GetNumberOfTracks();
      // Initialize Shericity matrix
      Float_t mxx    = 0.;
      Float_t myy    = 0.;
      Float_t mxy    = 0.;
      Int_t   nc     = 0;
      Float_t sump2  = 0.;
      Float_t ptMax  = 0.;
      Float_t etaMax = 0.;
      Float_t phiMax = 0.;
      Int_t   iMax   = -1;
      Float_t ptsum  = 0.;
      
      for (Int_t i = 0; i < nT; i++) {
	  AliAODTrack* track = aodE->GetTrack(i);
//
//    Track quality cuts
	  if (!track->TestFilterBit(1<<4)) continue;
	  fPtTH->Fill(track->Pt());

//
//
	  TVector3 pp(track->Px(), track->Py(), track->Pz());
	      
	  Float_t phi = track->Phi();
	  Float_t eta = track->Eta();
	  Float_t pt  = track->Pt();
	  Float_t jT  = pp.Perp(ppJ1);
	  
	  
	  // pT > 2.
	  
	    
	  Float_t deta = eta - etaJ;
	  Float_t dphi = phi - phiJ;

	  if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
	  if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
	  Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);

	  if (r < 0.4) ptsum += pt;
	  if (r < 0.5 && jT > 1.) {
	      TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	      TVector3 pPerp = pp - pLong;

	      Float_t ppjX = pPerp.Dot(ppJ2);
	      Float_t ppjY = pPerp.Dot(ppJ3);
	      Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);

	      mxx += (ppjX * ppjX / ppjT);
	      myy += (ppjY * ppjY / ppjT);
	      mxy += (ppjX * ppjY / ppjT);
	      nc++;
	      sump2 += ppjT;
	      // max pt
	      if (pt > ptMax) {
		  ptMax  = pt;
		  iMax   = i;
		  etaMax = deta;
		  phiMax = dphi;
	      }
	  } // R < 0.4
      } // 1st Track Loop
      fPtSum->Fill(ptJ, ptsum);
      
//
// At this point we have mxx, myy, mxy
      if (nc == 0) return;      
// Shericity Matrix	
      const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};	
      TMatrixDSym m0(2, ele);
// Find eigenvectors
      TMatrixDSymEigen m(m0);
      TVectorD eval(2);
      TMatrixD evecm = m.GetEigenVectors();
      eval  = m.GetEigenValues();
// Largest eigenvector
      Int_t jev = 0;
      if (eval[0] < eval[1]) jev = 1;
      TVectorD evec0(2);
// Principle axis
      evec0 = TMatrixDColumn(evecm, jev);
      TVector2 evec(evec0[0], evec0[1]); 
// Principle axis from leading partice
      Float_t phiM = TMath::ATan2(phiMax, etaMax);
      TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM)); 
      Float_t phistM = evecM.DeltaPhi(evec);
      if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);

//
// 2nd correlation with principal axis
//   
      for (Int_t i = 0; i < nT; i++) {
	  AliAODTrack* track = aodE->GetTrack(i);
//
//    Track quality cuts
	  if (!track->TestFilterBit(1<<4)) continue;
//
//
	  TVector3 pp(track->Px(), track->Py(), track->Pz());
	  Float_t phi = track->Phi();
	  Float_t eta = track->Eta();
	  Float_t pt  = track->Pt();
	  Float_t jT  = pp.Perp(ppJ1);
	  // pT < 2.
	  if (jT > 1.) continue;	    
	  Float_t deta = eta - etaJ;
	  Float_t dphi = phi - phiJ;
	  if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
	  if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;

	  Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);

	  if (r > 0.7) continue;

	  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	  TVector3 pPerp = pp - pLong;
	  
	  Float_t ppjX = pPerp.Dot(ppJ2);
	  Float_t ppjY = pPerp.Dot(ppJ3);
	  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	    
	  TVector2 vr(ppjX, ppjY) ;
	  Float_t phistr = evec.DeltaPhi(vr);
	  fPhiM->Fill(phistr);
	  fPhiMPt ->Fill(phistr, pt);
	  fPhiMPtJ->Fill(phistr, ptJ);
	  
      } // 2nd Track loop
      
  
  // Post output data.
  PostData(1, fHists);
}      

//________________________________________________________________________
void AliAnalysisTaskJetsTM::Terminate(Option_t *) 
{
// Terminate
}
Commit	Line	Data
	1	#include "TH1F.h"
	2	#include "TH2F.h"
	3	#include "TProfile.h"
	4	#include "TClonesArray.h"
	5	#include "TCanvas.h"
	6	#include "TList.h"
	7	#include "TMatrixD.h"
	8	#include "TMatrixDSym.h"
	9	#include "TMatrixDSymEigen.h"
	10
	11	#include "AliAODEvent.h"
	12	#include "AliAODTrack.h"
	13	#include "AliAODJet.h"
	14	#include "AliAnalysisTaskJetsTM.h"
	15
	16	ClassImp(AliAnalysisTaskJetsTM)
	17	//
	18	// Thrust Major (TM) analysis of reconstructed jets.
	19	// TM is the thrust in the plane perpendicular to the jet axis
	20	// The present amalysis performs the following steps:
	21	// (a) Construct to orthogonal unit vectors (e1, e2) in the plane perpendicular to the jet axis
	22	// (b) Calculate the components of all particles with jT > 1 GeV with respect to e1, e2
	23	// (c) Construct the sphericity matrix
	24	// (d) Find the two orthogonal eigenvectors of the spericity matrix
	25	// (e) Caluclate the components of all particles with jT < 1 GeV in the reference frame spanned by the eigenvectors
	26	// (f) Calculate the azimuthal angle in this frame
	27	//
	28	//
	29	// Author: andreas.morsch@cern.ch
	30
	31
	32	AliAnalysisTaskJetsTM::AliAnalysisTaskJetsTM(const char *name)
	33	: AliAnalysisTaskSE(name)
	34	,fHists(0)
	35	,fPtH(0)
	36	,fPtTH(0)
	37	,fPhiM(0)
	38	,fPhiMPt(0)
	39	,fPhiMPtJ(0)
	40	,fPtSum(0)
	41	{
	42	//
	43	// Constructor
	44	//
	45	DefineOutput(1, TList::Class());
	46	}
	47
	48
	49	//________________________________________________________________________
	50	void AliAnalysisTaskJetsTM::UserCreateOutputObjects()
	51	{
	52	// Create histograms
	53	// Called once
	54	fHists = new TList();
	55	fPtH = new TH1F("fPtH" , " p_{T} distribution of jets" , 50, 0., 100.0);
	56	fPtSum = new TH2F("fPtSum" , " p_{T} sum in cone" , 50, 0., 100.0, 50, 0, 100.);
	57	fPtTH = new TH1F("fPtTH", " p_{T} distribution of tracks", 50, 0., 100.0);
	58	fPhiM = new TH1F("fPhiM", " phi^{major} distribution of tracks", 64, 0., 3.2);
	59	fPhiMPt = new TH2F("fPhiMPt", " phi^{major} distribution of tracks vs pt", 64, 0., 3.2, 20, 0., 2.);
	60	fPhiMPtJ = new TH2F("fPhiMPtJ", " phi^{major} distribution of tracks vs pt Jet", 64, 0., 3.2, 40, 0., 100.);
	61
	62	fHists->SetOwner();
	63	fHists->Add(fPtH);
	64	fHists->Add(fPtSum);
	65	fHists->Add(fPtTH);
	66	fHists->Add(fPhiM);
	67	fHists->Add(fPhiMPt);
	68	fHists->Add(fPhiMPtJ);
	69	}
	70
	71
	72	//________________________________________________________________________
	73	void AliAnalysisTaskJetsTM::UserExec(Option_t *)
	74	{
	75	// Main loop
	76	// Called for each event
	77
	78	if (!fInputEvent) {
	79	Printf("ERROR: AOD not available");
	80	return;
	81	}
	82
	83	AliAODEvent* aodE = dynamic_cast<AliAODEvent*> (fInputEvent);
	84	TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
	85	Int_t nJ = jets->GetEntries();
	86
	87	Float_t ptmax = 0.;
	88	Int_t imax = -1;
	89
	90	//
	91	// Find highest pT jet with pt > 20 GeV
	92	//
	93	for (Int_t i = 0; i < nJ; i++) {
	94	AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
	95	Float_t ptJ = jet->Pt();
	96	Float_t etaJ = TMath::Abs(jet->Eta());
	97	if ((ptJ > 20.) && (ptJ > ptmax) && etaJ < 0.5) {
	98	ptmax = ptJ;
	99	imax = i;
	100	}
	101	}
	102
	103	if (imax == -1) return;
	104
	105
	106	AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
	107	Float_t ptJ = jet->Pt();
	108	Float_t phiJ = jet->Phi();
	109	Float_t etaJ = jet->Eta();
	110
	111	fPtH->Fill(ptJ);
	112
	113	//
	114	// The transverse plane
	115	//
	116	// 2 normalized vectors perpendicular to the jet direction
	117	//
	118	Float_t pxJ = jet->Px();
	119	Float_t pyJ = jet->Py();
	120	Float_t pzJ = jet->Pz();
	121
	122	TVector3 ppJ1(pxJ, pyJ, pzJ);
	123	TVector3 ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
	124	ppJ3.SetMag(1.);
	125	TVector3 ppJ2(-pyJ, pxJ, 0);
	126	ppJ2.SetMag(1.);
	127
	128	//
	129	// 1st track loop to determine the sphericity matrix
	130	//
	131	Int_t nT = aodE->GetNumberOfTracks();
	132	// Initialize Shericity matrix
	133	Float_t mxx = 0.;
	134	Float_t myy = 0.;
	135	Float_t mxy = 0.;
	136	Int_t nc = 0;
	137	Float_t sump2 = 0.;
	138	Float_t ptMax = 0.;
	139	Float_t etaMax = 0.;
	140	Float_t phiMax = 0.;
	141	Int_t iMax = -1;
	142	Float_t ptsum = 0.;
	143
	144	for (Int_t i = 0; i < nT; i++) {
	145	AliAODTrack* track = aodE->GetTrack(i);
	146	//
	147	// Track quality cuts
	148	if (!track->TestFilterBit(1<<4)) continue;
	149	fPtTH->Fill(track->Pt());
	150
	151	//
	152	//
	153	TVector3 pp(track->Px(), track->Py(), track->Pz());
	154
	155	Float_t phi = track->Phi();
	156	Float_t eta = track->Eta();
	157	Float_t pt = track->Pt();
	158	Float_t jT = pp.Perp(ppJ1);
	159
	160
	161	// pT > 2.
	162
	163
	164	Float_t deta = eta - etaJ;
	165	Float_t dphi = phi - phiJ;
	166
	167	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
	168	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
	169	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
	170
	171	if (r < 0.4) ptsum += pt;
	172	if (r < 0.5 && jT > 1.) {
	173	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	174	TVector3 pPerp = pp - pLong;
	175
	176	Float_t ppjX = pPerp.Dot(ppJ2);
	177	Float_t ppjY = pPerp.Dot(ppJ3);
	178	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	179
	180	mxx += (ppjX * ppjX / ppjT);
	181	myy += (ppjY * ppjY / ppjT);
	182	mxy += (ppjX * ppjY / ppjT);
	183	nc++;
	184	sump2 += ppjT;
	185	// max pt
	186	if (pt > ptMax) {
	187	ptMax = pt;
	188	iMax = i;
	189	etaMax = deta;
	190	phiMax = dphi;
	191	}
	192	} // R < 0.4
	193	} // 1st Track Loop
	194	fPtSum->Fill(ptJ, ptsum);
	195
	196	//
	197	// At this point we have mxx, myy, mxy
	198	if (nc == 0) return;
	199	// Shericity Matrix
	200	const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
	201	TMatrixDSym m0(2, ele);
	202	// Find eigenvectors
	203	TMatrixDSymEigen m(m0);
	204	TVectorD eval(2);
	205	TMatrixD evecm = m.GetEigenVectors();
	206	eval = m.GetEigenValues();
	207	// Largest eigenvector
	208	Int_t jev = 0;
	209	if (eval[0] < eval[1]) jev = 1;
	210	TVectorD evec0(2);
	211	// Principle axis
	212	evec0 = TMatrixDColumn(evecm, jev);
	213	TVector2 evec(evec0[0], evec0[1]);
	214	// Principle axis from leading partice
	215	Float_t phiM = TMath::ATan2(phiMax, etaMax);
	216	TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
	217	Float_t phistM = evecM.DeltaPhi(evec);
	218	if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
	219
	220	//
	221	// 2nd correlation with principal axis
	222	//
	223	for (Int_t i = 0; i < nT; i++) {
	224	AliAODTrack* track = aodE->GetTrack(i);
	225	//
	226	// Track quality cuts
	227	if (!track->TestFilterBit(1<<4)) continue;
	228	//
	229	//
	230	TVector3 pp(track->Px(), track->Py(), track->Pz());
	231	Float_t phi = track->Phi();
	232	Float_t eta = track->Eta();
	233	Float_t pt = track->Pt();
	234	Float_t jT = pp.Perp(ppJ1);
	235	// pT < 2.
	236	if (jT > 1.) continue;
	237	Float_t deta = eta - etaJ;
	238	Float_t dphi = phi - phiJ;
	239	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
	240	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
	241
	242	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
	243
	244	if (r > 0.7) continue;
	245
	246	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	247	TVector3 pPerp = pp - pLong;
	248
	249	Float_t ppjX = pPerp.Dot(ppJ2);
	250	Float_t ppjY = pPerp.Dot(ppJ3);
	251	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	252
	253	TVector2 vr(ppjX, ppjY) ;
	254	Float_t phistr = evec.DeltaPhi(vr);
	255	fPhiM->Fill(phistr);
	256	fPhiMPt ->Fill(phistr, pt);
	257	fPhiMPtJ->Fill(phistr, ptJ);
	258
	259	} // 2nd Track loop
	260
	261
	262	// Post output data.
	263	PostData(1, fHists);
	264	}
	265
	266	//________________________________________________________________________
	267	void AliAnalysisTaskJetsTM::Terminate(Option_t *)
	268	{
	269	// Terminate
	270	}