[u/mrichter/AliRoot.git] / PWGJE / 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);

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

  TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
  if (!jets) {
    Printf("ERROR: Jet branch not available");
    return;
  }

  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));
      if (!jet) continue;
      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 = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
	  if(!track) AliFatal("Not a standard AOD");
//
//    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 = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
	  if(!track) AliFatal("Not a standard AOD");
//
//    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);
	    
	  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
0eba9d77	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);
139a69f1	84
	85	if (!aodE) {
	86	Printf("ERROR: AOD not available");
	87	return;
	88	}
	89
0eba9d77	90	TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
139a69f1	91	if (!jets) {
	92	Printf("ERROR: Jet branch not available");
	93	return;
	94	}
	95
0eba9d77	96	Int_t nJ = jets->GetEntries();
	97
	98	Float_t ptmax = 0.;
	99	Int_t imax = -1;
	100
	101	//
	102	// Find highest pT jet with pt > 20 GeV
	103	//
	104	for (Int_t i = 0; i < nJ; i++) {
	105	AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
781ef0e6	106	if (!jet) continue;
0eba9d77	107	Float_t ptJ = jet->Pt();
	108	Float_t etaJ = TMath::Abs(jet->Eta());
	109	if ((ptJ > 20.) && (ptJ > ptmax) && etaJ < 0.5) {
	110	ptmax = ptJ;
	111	imax = i;
	112	}
	113	}
	114
	115	if (imax == -1) return;
	116
	117
	118	AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
	119	Float_t ptJ = jet->Pt();
	120	Float_t phiJ = jet->Phi();
	121	Float_t etaJ = jet->Eta();
	122
	123	fPtH->Fill(ptJ);
	124
	125	//
	126	// The transverse plane
	127	//
	128	// 2 normalized vectors perpendicular to the jet direction
	129	//
	130	Float_t pxJ = jet->Px();
	131	Float_t pyJ = jet->Py();
	132	Float_t pzJ = jet->Pz();
	133
	134	TVector3 ppJ1(pxJ, pyJ, pzJ);
	135	TVector3 ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
	136	ppJ3.SetMag(1.);
	137	TVector3 ppJ2(-pyJ, pxJ, 0);
	138	ppJ2.SetMag(1.);
	139
	140	//
	141	// 1st track loop to determine the sphericity matrix
	142	//
	143	Int_t nT = aodE->GetNumberOfTracks();
	144	// Initialize Shericity matrix
	145	Float_t mxx = 0.;
	146	Float_t myy = 0.;
	147	Float_t mxy = 0.;
	148	Int_t nc = 0;
	149	Float_t sump2 = 0.;
	150	Float_t ptMax = 0.;
	151	Float_t etaMax = 0.;
	152	Float_t phiMax = 0.;
	153	Int_t iMax = -1;
	154	Float_t ptsum = 0.;
	155
	156	for (Int_t i = 0; i < nT; i++) {
f15c1f69	157	AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
f15c1f69	158	if(!track) AliFatal("Not a standard AOD");
0eba9d77	159	//
	160	// Track quality cuts
	161	if (!track->TestFilterBit(1<<4)) continue;
	162	fPtTH->Fill(track->Pt());
	163
	164	//
	165	//
	166	TVector3 pp(track->Px(), track->Py(), track->Pz());
	167
	168	Float_t phi = track->Phi();
	169	Float_t eta = track->Eta();
	170	Float_t pt = track->Pt();
	171	Float_t jT = pp.Perp(ppJ1);
	172
	173
	174	// pT > 2.
	175
	176
	177	Float_t deta = eta - etaJ;
	178	Float_t dphi = phi - phiJ;
	179
	180	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
	181	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
	182	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
	183
	184	if (r < 0.4) ptsum += pt;
	185	if (r < 0.5 && jT > 1.) {
	186	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	187	TVector3 pPerp = pp - pLong;
	188
	189	Float_t ppjX = pPerp.Dot(ppJ2);
	190	Float_t ppjY = pPerp.Dot(ppJ3);
	191	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	192
	193	mxx += (ppjX * ppjX / ppjT);
	194	myy += (ppjY * ppjY / ppjT);
	195	mxy += (ppjX * ppjY / ppjT);
	196	nc++;
	197	sump2 += ppjT;
	198	// max pt
	199	if (pt > ptMax) {
	200	ptMax = pt;
	201	iMax = i;
	202	etaMax = deta;
	203	phiMax = dphi;
	204	}
	205	} // R < 0.4
	206	} // 1st Track Loop
	207	fPtSum->Fill(ptJ, ptsum);
	208
	209	//
	210	// At this point we have mxx, myy, mxy
	211	if (nc == 0) return;
	212	// Shericity Matrix
	213	const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
	214	TMatrixDSym m0(2, ele);
	215	// Find eigenvectors
	216	TMatrixDSymEigen m(m0);
	217	TVectorD eval(2);
	218	TMatrixD evecm = m.GetEigenVectors();
	219	eval = m.GetEigenValues();
	220	// Largest eigenvector
	221	Int_t jev = 0;
	222	if (eval[0] < eval[1]) jev = 1;
223	TVectorD evec0(2);
224	// Principle axis
225	evec0 = TMatrixDColumn(evecm, jev);
226	TVector2 evec(evec0[0], evec0[1]);
227	// Principle axis from leading partice
228	Float_t phiM = TMath::ATan2(phiMax, etaMax);
229	TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
230	Float_t phistM = evecM.DeltaPhi(evec);
231	if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
232
233	//
234	// 2nd correlation with principal axis
235	//
236	for (Int_t i = 0; i < nT; i++) {
f15c1f69	237	AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
f15c1f69	238	if(!track) AliFatal("Not a standard AOD");
0eba9d77	239	//
	240	// Track quality cuts
	241	if (!track->TestFilterBit(1<<4)) continue;
	242	//
	243	//
	244	TVector3 pp(track->Px(), track->Py(), track->Pz());
	245	Float_t phi = track->Phi();
	246	Float_t eta = track->Eta();
	247	Float_t pt = track->Pt();
	248	Float_t jT = pp.Perp(ppJ1);
	249	// pT < 2.
	250	if (jT > 1.) continue;
	251	Float_t deta = eta - etaJ;
	252	Float_t dphi = phi - phiJ;
	253	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
	254	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
	255
	256	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
	257
	258	if (r > 0.7) continue;
	259
	260	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	261	TVector3 pPerp = pp - pLong;
	262
	263	Float_t ppjX = pPerp.Dot(ppJ2);
	264	Float_t ppjY = pPerp.Dot(ppJ3);
0eba9d77	265
	266	TVector2 vr(ppjX, ppjY) ;
	267	Float_t phistr = evec.DeltaPhi(vr);
	268	fPhiM->Fill(phistr);
	269	fPhiMPt ->Fill(phistr, pt);
	270	fPhiMPtJ->Fill(phistr, ptJ);
	271
	272	} // 2nd Track loop
	273
	274
	275	// Post output data.
	276	PostData(1, fHists);
	277	}
	278
	279	//________________________________________________________________________
	280	void AliAnalysisTaskJetsTM::Terminate(Option_t *)
	281	{
	282	// Terminate
	283	}