[u/mrichter/AliRoot.git] / EMCAL / macros / TestESD.C

#if !defined(__CINT__) || defined(__MAKECINT__)

//Root include files 
#include <Riostream.h>
#include <TFile.h>
#include <TChain.h>
#include <TParticle.h>
#include <TNtuple.h>
#include <TCanvas.h>
#include <TObjArray.h>
#include <TSystem.h>
#include <TString.h>
#include <TH1F.h>
#include <TVector.h>
#include <TParticle.h>
#include <TRefArray.h>
#include <TArrayS.h>

//AliRoot include files 
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliPID.h"
#include "AliLog.h"

#endif

void TestESD() {

  TFile* f = new TFile("AliESDs.root");
  TTree* esdTree = (TTree*)f->Get("esdTree");
  
  AliESDEvent* esd = new AliESDEvent();
  esd->ReadFromTree(esdTree);

  Int_t nEvt = esdTree->GetEntries();
  Float_t pos[3] = {0.,0.,0.};

  for(Int_t iev = 0; iev < nEvt; iev++) {
    cout << "Event: " << iev+1 << "/" << nEvt << endl;
    esdTree->GetEvent(iev);

    TRefArray* caloClusters = new TRefArray();
    esd->GetEMCALClusters(caloClusters);

    //get reconstructed vertex position
    Double_t vertex_position[3];
    esd->GetVertex()->GetXYZ(vertex_position);

    //loop over clusters
    Int_t nclus = caloClusters->GetEntries();
    for (Int_t icl = 0; icl < nclus; icl++) {

      AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
      Float_t energy = clus->E();
      clus->GetPosition(pos);
      TVector3 vpos(pos[0],pos[1],pos[2]);
      TLorentzVector p;
      clus->GetMomentum(p,vertex_position);
      Double_t cphi = vpos.Phi();
      Double_t ceta = vpos.Eta();

      Int_t nMatched = clus->GetNTracksMatched();
      Int_t trackIndex = clus->GetTrackMatched();
      Int_t nLabels = clus->GetNLabels();
      Int_t labelIndex = clus->GetLabel();

      Int_t nCells = clus->GetNCells();

      //For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME

      cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: " 
	   << cphi << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched 
	   << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: " 
	   << labelIndex << endl;

      if(trackIndex >= 0) {
	AliESDtrack* track = esd->GetTrack(trackIndex);
	Double_t tphi = track->GetOuterParam()->Phi();
	Double_t teta = track->GetOuterParam()->Eta();
	Double_t tmom = track->GetOuterParam()->P();
	cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;

	Double_t deta = teta - ceta;
	Double_t dphi = tphi - cphi;
	if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
	if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
	Double_t dR = sqrt(dphi*dphi + deta*deta);

	Double_t pOverE = tmom/energy;

	if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
	  cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
	}

      }

    }


  }


}
Commit	Line	Data
cd259404	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2
	3	//Root include files
	4	#include <Riostream.h>
	5	#include <TFile.h>
	6	#include <TChain.h>
	7	#include <TParticle.h>
	8	#include <TNtuple.h>
	9	#include <TCanvas.h>
	10	#include <TObjArray.h>
	11	#include <TSystem.h>
	12	#include <TString.h>
	13	#include <TH1F.h>
	14	#include <TVector.h>
	15	#include <TParticle.h>
	16	#include <TRefArray.h>
	17	#include <TArrayS.h>
	18
	19	//AliRoot include files
	20	#include "AliRunLoader.h"
	21	#include "AliStack.h"
	22	#include "AliESDEvent.h"
	23	#include "AliESDVertex.h"
	24	#include "AliESDCaloCluster.h"
	25	#include "AliESDCaloCells.h"
	26	#include "AliPID.h"
	27	#include "AliLog.h"
	28
	29	#endif
	30
	31	void TestESD() {
	32
	33	TFile* f = new TFile("AliESDs.root");
	34	TTree* esdTree = (TTree*)f->Get("esdTree");
	35
	36	AliESDEvent* esd = new AliESDEvent();
	37	esd->ReadFromTree(esdTree);
	38
	39	Int_t nEvt = esdTree->GetEntries();
	40	Float_t pos[3] = {0.,0.,0.};
	41
	42	for(Int_t iev = 0; iev < nEvt; iev++) {
	43	cout << "Event: " << iev+1 << "/" << nEvt << endl;
	44	esdTree->GetEvent(iev);
	45
	46	TRefArray* caloClusters = new TRefArray();
	47	esd->GetEMCALClusters(caloClusters);
	48
	49	//get reconstructed vertex position
	50	Double_t vertex_position[3];
	51	esd->GetVertex()->GetXYZ(vertex_position);
	52
	53	//loop over clusters
	54	Int_t nclus = caloClusters->GetEntries();
	55	for (Int_t icl = 0; icl < nclus; icl++) {
	56
	57	AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
	58	Float_t energy = clus->E();
	59	clus->GetPosition(pos);
	60	TVector3 vpos(pos[0],pos[1],pos[2]);
	61	TLorentzVector p;
	62	clus->GetMomentum(p,vertex_position);
	63	Double_t cphi = vpos.Phi();
	64	Double_t ceta = vpos.Eta();
65
66	Int_t nMatched = clus->GetNTracksMatched();
67	Int_t trackIndex = clus->GetTrackMatched();
68	Int_t nLabels = clus->GetNLabels();
69	Int_t labelIndex = clus->GetLabel();
70
71	Int_t nCells = clus->GetNCells();
72
73	//For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME
74
75	cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: "
76	<< cphi << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched
77	<< " Index: " << trackIndex << " #Labels: " << nLabels << " Index: "
78	<< labelIndex << endl;
79
80	if(trackIndex >= 0) {
81	AliESDtrack* track = esd->GetTrack(trackIndex);
82	Double_t tphi = track->GetOuterParam()->Phi();
83	Double_t teta = track->GetOuterParam()->Eta();
84	Double_t tmom = track->GetOuterParam()->P();
85	cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
86
87	Double_t deta = teta - ceta;
88	Double_t dphi = tphi - cphi;
89	if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
90	if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
91	Double_t dR = sqrt(dphidphi + detadeta);
92
93	Double_t pOverE = tmom/energy;
94
95	if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
96	cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
97	}
98
99	}
100
101	}
102
103
104	}
105
106
107
108	}