[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 = 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);
	    
	  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++) {
	157	AliAODTrack* track = aodE->GetTrack(i);
	158	//
	159	// Track quality cuts
	160	if (!track->TestFilterBit(1<<4)) continue;
	161	fPtTH->Fill(track->Pt());
	162
	163	//
	164	//
	165	TVector3 pp(track->Px(), track->Py(), track->Pz());
	166
	167	Float_t phi = track->Phi();
	168	Float_t eta = track->Eta();
	169	Float_t pt = track->Pt();
	170	Float_t jT = pp.Perp(ppJ1);
171
172
173	// pT > 2.
174
175
176	Float_t deta = eta - etaJ;
177	Float_t dphi = phi - phiJ;
178
179	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
180	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
181	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
182
183	if (r < 0.4) ptsum += pt;
184	if (r < 0.5 && jT > 1.) {
185	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
186	TVector3 pPerp = pp - pLong;
187
188	Float_t ppjX = pPerp.Dot(ppJ2);
189	Float_t ppjY = pPerp.Dot(ppJ3);
190	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
191
192	mxx += (ppjX * ppjX / ppjT);
193	myy += (ppjY * ppjY / ppjT);
194	mxy += (ppjX * ppjY / ppjT);
195	nc++;
196	sump2 += ppjT;
197	// max pt
198	if (pt > ptMax) {
199	ptMax = pt;
200	iMax = i;
201	etaMax = deta;
202	phiMax = dphi;
203	}
204	} // R < 0.4
205	} // 1st Track Loop
206	fPtSum->Fill(ptJ, ptsum);
207
208	//
209	// At this point we have mxx, myy, mxy
210	if (nc == 0) return;
211	// Shericity Matrix
212	const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
213	TMatrixDSym m0(2, ele);
214	// Find eigenvectors
215	TMatrixDSymEigen m(m0);
216	TVectorD eval(2);
217	TMatrixD evecm = m.GetEigenVectors();
218	eval = m.GetEigenValues();
219	// Largest eigenvector
220	Int_t jev = 0;
221	if (eval[0] < eval[1]) jev = 1;
222	TVectorD evec0(2);
223	// Principle axis
224	evec0 = TMatrixDColumn(evecm, jev);
225	TVector2 evec(evec0[0], evec0[1]);
226	// Principle axis from leading partice
227	Float_t phiM = TMath::ATan2(phiMax, etaMax);
228	TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
229	Float_t phistM = evecM.DeltaPhi(evec);
230	if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
231
232	//
233	// 2nd correlation with principal axis
234	//
235	for (Int_t i = 0; i < nT; i++) {
236	AliAODTrack* track = aodE->GetTrack(i);
237	//
238	// Track quality cuts
239	if (!track->TestFilterBit(1<<4)) continue;
240	//
241	//
242	TVector3 pp(track->Px(), track->Py(), track->Pz());
243	Float_t phi = track->Phi();
244	Float_t eta = track->Eta();
245	Float_t pt = track->Pt();
246	Float_t jT = pp.Perp(ppJ1);
247	// pT < 2.
248	if (jT > 1.) continue;
249	Float_t deta = eta - etaJ;
250	Float_t dphi = phi - phiJ;
251	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
252	if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
253
254	Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
255
256	if (r > 0.7) continue;
257
258	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
259	TVector3 pPerp = pp - pLong;
260
261	Float_t ppjX = pPerp.Dot(ppJ2);
262	Float_t ppjY = pPerp.Dot(ppJ3);
0eba9d77	263
	264	TVector2 vr(ppjX, ppjY) ;
	265	Float_t phistr = evec.DeltaPhi(vr);
	266	fPhiM->Fill(phistr);
	267	fPhiMPt ->Fill(phistr, pt);
	268	fPhiMPtJ->Fill(phistr, ptJ);
	269
	270	} // 2nd Track loop
	271
	272
	273	// Post output data.
	274	PostData(1, fHists);
	275	}
	276
	277	//________________________________________________________________________
	278	void AliAnalysisTaskJetsTM::Terminate(Option_t *)
	279	{
	280	// Terminate
	281	}