[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

//Change the bool depending on what information you want to print
// when all FALSE, prints minimum cluster information.
Bool_t kPrintKine = kTRUE; //Do not use for raw data.
Bool_t kPrintCaloCells = kFALSE;
Bool_t kPrintTrackMatches = kFALSE;
Bool_t kPrintClusterCells = kFALSE;

void TestESD() {

  TGeoManager::Import("geometry.root");
  AliEMCALGeometry *geom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETE") ;  

  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);
	  
	//In case you want to play with MC data
	AliStack *stack = 0;
	if(kPrintKine){  
		AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),  "read");
		rl->LoadKinematics();  
		rl->GetEvent(iev);
		stack = rl->Stack();
	}  
  	  
	//get reconstructed vertex position
	Double_t vertex_position[3];
	esd->GetVertex()->GetXYZ(vertex_position);
	  
	//GetCellsClusters Array  
    AliESDCaloCells &cells= *(esd->GetEMCALCells());
	  // Uncomment to see the full list of digit amplitudes and times.
	if(kPrintCaloCells){  
		Int_t nTotalCells = cells.GetNumberOfCells() ;  
		Int_t type        = cells.GetType();
		for (Int_t icell=  0; icell <  nTotalCells; icell++) {
			cout<<"Cell   : "<<icell<<"/"<<nTotalCells<<" ID: "<<cells.GetCellNumber(icell)<<" Amplitude: "<<cells.GetAmplitude(icell)<<" Time: "<<cells.GetTime(icell)<<endl;	  
		}// cell loop
	}
	  
	//GetCaloClusters Array
    TRefArray* caloClusters = new TRefArray();
    esd->GetEMCALClusters(caloClusters);

    //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*TMath::RadToDeg() << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched 
	   << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: " 
	   << labelIndex << endl;

	if(kPrintTrackMatches && 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;
		}

	}
		
	//Get CaloCells of cluster and print
	if(kPrintClusterCells){	
		UShort_t * index = clus->GetCellsAbsId() ;
		Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
		for(Int_t i = 0; i < nCells ; i++){
			Int_t absId =   index[i]; // or clus->GetCellNumber(i) ;
			Double_t ampFract =  fraction[i];
			Float_t amp       = cells.GetCellAmplitude(absId) ;
			Float_t time      = cells.GetCellTime(absId);
			cout<<"         Cluster Cell: AbsID : "<< absId << "; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time<<endl;
			//Geometry methods  
			Int_t iSupMod =  0 ;
			Int_t iTower  =  0 ;
			Int_t iIphi   =  0 ;
			Int_t iIeta   =  0 ;
			Int_t iphi    =  0 ;
			Int_t ieta    =  0 ;
			if(geom){
			  geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
			  //Gives SuperModule and Tower numbers
			  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
							    iIphi, iIeta,iphi,ieta);
			  //Gives label of cell in eta-phi position per each supermodule
			  cout<< "SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
			      <<"; Phi "<<iIphi<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
			}
		}
	}
		
	//Print primary info
	if(!stack || !kPrintKine) continue;
	if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
		TParticle * particle = stack->Particle(labelIndex);
		//Print primary values
		cout<<"         More  contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<"; Eta "<<particle->Eta()<<"; Phi "<<particle->Phi()*TMath::RadToDeg()<<endl;   
		for(Int_t i = 1; i < nLabels; i++){
			particle = stack->Particle((clus->GetLabels())->At(i));
			cout<<"         Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
		}
	}
	else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
	else cout<<"Negative label!!!  : "<<labelIndex<<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
9480984d	31	//Change the bool depending on what information you want to print
	32	// when all FALSE, prints minimum cluster information.
	33	Bool_t kPrintKine = kTRUE; //Do not use for raw data.
	34	Bool_t kPrintCaloCells = kFALSE;
	35	Bool_t kPrintTrackMatches = kFALSE;
	36	Bool_t kPrintClusterCells = kFALSE;
	37
cd259404	38	void TestESD() {
cd259404	39
a6bf15ef	40	TGeoManager::Import("geometry.root");
	41	AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETE") ;
	42
cd259404	43	TFile* f = new TFile("AliESDs.root");
	44	TTree* esdTree = (TTree*)f->Get("esdTree");
	45
	46	AliESDEvent* esd = new AliESDEvent();
	47	esd->ReadFromTree(esdTree);
	48
	49	Int_t nEvt = esdTree->GetEntries();
	50	Float_t pos[3] = {0.,0.,0.};
	51
	52	for(Int_t iev = 0; iev < nEvt; iev++) {
9480984d	53	cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
cd259404	54	esdTree->GetEvent(iev);
9480984d	55
	56	//In case you want to play with MC data
	57	AliStack *stack = 0;
	58	if(kPrintKine){
	59	AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(), "read");
	60	rl->LoadKinematics();
	61	rl->GetEvent(iev);
	62	stack = rl->Stack();
	63	}
	64
	65	//get reconstructed vertex position
	66	Double_t vertex_position[3];
	67	esd->GetVertex()->GetXYZ(vertex_position);
	68
	69	//GetCellsClusters Array
	70	AliESDCaloCells &cells= *(esd->GetEMCALCells());
	71	// Uncomment to see the full list of digit amplitudes and times.
	72	if(kPrintCaloCells){
	73	Int_t nTotalCells = cells.GetNumberOfCells() ;
	74	Int_t type = cells.GetType();
	75	for (Int_t icell= 0; icell < nTotalCells; icell++) {
	76	cout<<"Cell : "<<icell<<"/"<<nTotalCells<<" ID: "<<cells.GetCellNumber(icell)<<" Amplitude: "<<cells.GetAmplitude(icell)<<" Time: "<<cells.GetTime(icell)<<endl;
	77	}// cell loop
	78	}
	79
	80	//GetCaloClusters Array
cd259404	81	TRefArray* caloClusters = new TRefArray();
	82	esd->GetEMCALClusters(caloClusters);
	83
cd259404	84	//loop over clusters
	85	Int_t nclus = caloClusters->GetEntries();
	86	for (Int_t icl = 0; icl < nclus; icl++) {
	87
	88	AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
	89	Float_t energy = clus->E();
	90	clus->GetPosition(pos);
	91	TVector3 vpos(pos[0],pos[1],pos[2]);
	92	TLorentzVector p;
	93	clus->GetMomentum(p,vertex_position);
	94	Double_t cphi = vpos.Phi();
	95	Double_t ceta = vpos.Eta();
	96
	97	Int_t nMatched = clus->GetNTracksMatched();
	98	Int_t trackIndex = clus->GetTrackMatched();
	99	Int_t nLabels = clus->GetNLabels();
	100	Int_t labelIndex = clus->GetLabel();
	101
	102	Int_t nCells = clus->GetNCells();
	103
	104	//For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME
	105
	106	cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: "
a6bf15ef	107	<< cphi*TMath::RadToDeg() << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched
cd259404	108	<< " Index: " << trackIndex << " #Labels: " << nLabels << " Index: "
	109	<< labelIndex << endl;
	110
9480984d	111	if(kPrintTrackMatches && trackIndex >= 0) {
	112	AliESDtrack* track = esd->GetTrack(trackIndex);
	113	Double_t tphi = track->GetOuterParam()->Phi();
	114	Double_t teta = track->GetOuterParam()->Eta();
	115	Double_t tmom = track->GetOuterParam()->P();
	116	cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
cd259404	117
9480984d	118	Double_t deta = teta - ceta;
	119	Double_t dphi = tphi - cphi;
	120	if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
	121	if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
	122	Double_t dR = sqrt(dphidphi + detadeta);
cd259404	123
9480984d	124	Double_t pOverE = tmom/energy;
cd259404	125
9480984d	126	if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
	127	cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
	128	}
cd259404	129
9480984d	130	}
	131
	132	//Get CaloCells of cluster and print
	133	if(kPrintClusterCells){
	134	UShort_t * index = clus->GetCellsAbsId() ;
	135	Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
	136	for(Int_t i = 0; i < nCells ; i++){
	137	Int_t absId = index[i]; // or clus->GetCellNumber(i) ;
	138	Double_t ampFract = fraction[i];
	139	Float_t amp = cells.GetCellAmplitude(absId) ;
	140	Float_t time = cells.GetCellTime(absId);
	141	cout<<" Cluster Cell: AbsID : "<< absId << "; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time<<endl;
a6bf15ef	142	//Geometry methods
	143	Int_t iSupMod = 0 ;
	144	Int_t iTower = 0 ;
	145	Int_t iIphi = 0 ;
	146	Int_t iIeta = 0 ;
	147	Int_t iphi = 0 ;
	148	Int_t ieta = 0 ;
	149	if(geom){
	150	geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
	151	//Gives SuperModule and Tower numbers
	152	geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
	153	iIphi, iIeta,iphi,ieta);
	154	//Gives label of cell in eta-phi position per each supermodule
	155	cout<< "SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
	156	<<"; Phi "<<iIphi<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
	157	}
9480984d	158	}
	159	}
	160
	161	//Print primary info
	162	if(!stack \|\| !kPrintKine) continue;
	163	if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
	164	TParticle * particle = stack->Particle(labelIndex);
	165	//Print primary values
a6bf15ef	166	cout<<" More contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<"; Eta "<<particle->Eta()<<"; Phi "<<particle->Phi()*TMath::RadToDeg()<<endl;
9480984d	167	for(Int_t i = 1; i < nLabels; i++){
a6bf15ef	168	particle = stack->Particle((clus->GetLabels())->At(i));
9480984d	169	cout<<" Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
	170	}
	171	}
	172	else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
	173	else cout<<"Negative label!!! : "<<labelIndex<<endl;
	174
	175	}
cd259404	176
	177
	178	}
	179
	180
	181
	182	}