[u/mrichter/AliRoot.git] / JETAN / AliJetFillUnitArrayTracks.cxx


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//======================================================================
// *** November 2006
// Author: magali.estienne@ires.in2p3.fr
// 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and  
// EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance
//     2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment)
//                              fGrid==1, work with a part of the TPC acceptance  
// 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them
// 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//======================================================================
// ***September 2006
// TTask : Fill Unit Array for the Tracks information 
// Called by ESD reader for jet analysis
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//======================================================================
// *** July 2006
// 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex
// and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the
// index or the eta, phi position of a grid.
// 2) Define a grid for TPC
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//======================================================================
// ***July 2006
// Fill Unit Array class 
// Class used by AliJetESDReader to fill a UnitArray from the information extracted 
// from the particle tracks
// Author: magali.estienne@ires.in2p3.fr
//======================================================================


// --- Standard library ---
#include <Riostream.h>

// --- ROOT system ---
#include <TSystem.h>
#include <TLorentzVector.h>
#include <TVector3.h>
//#include "Math/Vector3D.h"
//#include "Math/Vector3Dfwd.h"
#include "TTask.h"
#include <TGeoManager.h>
#include <TMatrixD.h>
#include <TArrayD.h>

// --- AliRoot header files ---
#include "AliJetFinder.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
#include "AliESD.h"
#include "AliEMCALGeometry.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayTracks.h"
#include "AliJetGrid.h"

ClassImp(AliJetFillUnitArrayTracks)

AliEMCALGeometry *AliJetFillUnitArrayTracks::fGeom=0; 

//_____________________________________________________________________________
AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
  : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information")
{
  // constructor
   fNIn = 0;
  fOpt = 0;
  fDebug = 0;
  fNphi = 0;
  fNeta = 0;
  fPhiMin = 0;
  fPhiMax = 0;
  fEtaMin = 0;
  fEtaMax = 0;
  fEtaBinInTPCAcc = 0;
  fPhiBinInTPCAcc = 0;
  fEtaBinInEMCalAcc = 0;
  fPhiBinInEMCalAcc = 0;
  fNbinPhi = 0;
}

//____________________________________________________________________________
void AliJetFillUnitArrayTracks::SetEMCALGeometry()
{
  // Set EMCAL geometry information
  fGeom = AliEMCALGeometry::GetInstance();
  if (fGeom == 0)
    fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
  if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
}

//____________________________________________________________________________
void AliJetFillUnitArrayTracks::InitParameters()
{
  fHCorrection    = 0;     // For hadron correction
  fHadCorr        = 0;     // For hadron correction
  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
  fDebug = fReaderHeader->GetDebug();

  fEtaMinCal = fGeom->GetArm1EtaMin();
  fEtaMaxCal = fGeom->GetArm1EtaMax();
  fPhiMinCal = fGeom->GetArm1PhiMin();
  fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !

  if(fDebug>30){
    cout << "fEtaMinCal : " << fEtaMinCal << endl;
    cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
    cout << "fPhiMinCal : " << fPhiMinCal << endl;
    cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
  }

  fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
			fPhiMax,fEtaMin,fEtaMax);
  fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, 
			fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);

  fEta   = fTPCGrid->GetArrayEta();
  fPhi   = fTPCGrid->GetArrayPhi();
  fIndex = fTPCGrid->GetIndexObject();

  if(fDebug>20){
    for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
    for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
    
    for(Int_t i=0; i<fNphi+1; i++)
      for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
	  (*fIndex)(i,j) << endl; }
  } 
  if(fDebug>1) printf("\n Parameters initiated ! \n");
}

//_____________________________________________________________________________
AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
{
  // destructor
}

//_____________________________________________________________________________
void AliJetFillUnitArrayTracks::Exec(Option_t* option)
{
  //
  // Main method.
  // Explain

  fDebug = fReaderHeader->GetDebug();
  if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");  
  if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
  // Set EMCal Geometry
  SetEMCALGeometry();
  // Set parameters
  InitParameters();

  TClonesArray *lvArray = fMomentumArray;     // Correct checked !
  Int_t nInT =  lvArray->GetEntries();        // Correct checked !
  Float_t ptMin = fReaderHeader->GetPtCut();  // Correct checked !
  
  // sflag -> not yet implemented !!!!
  
  if(fDebug>3) cout << "nInT : " << nInT << endl;

  if (nInT == 0) return;
  
  // local arrays for input
  Float_t* enT  = new Float_t[nInT];
  Float_t* ptT  = new Float_t[nInT];
  Float_t* etaT = new Float_t[nInT];
  Float_t* phiT = new Float_t[nInT];
  Float_t* thetaT = new Float_t[nInT];
  
  Int_t trackInEmcalAcc = 0;
  Int_t trackInTpcAcc = 0;
  Int_t trackInTpcAccOnly = 0;

  Int_t nElements = fTPCGrid->GetNEntries();
  Int_t nElements2 = fEMCalGrid->GetNEntries();
  fGrid = fTPCGrid->GetGridType();

  if(fDebug>3){
    cout << "nElements : " << nElements << endl;
    cout << "nElements2 : " << nElements2 << endl;
    cout << "fNumUnits : " << fNumUnits << endl;
    cout << "sum : " << fNumUnits+nElements << endl;
  }

  // Set energy exactly to zero
  if(fGrid==0)
    for(Int_t k=0; k<nElements; k++)
      fUnitArray[k].SetUnitEnergy(0.); 
  
  if(fGrid==1)
    for(Int_t k=0; k<fNumUnits+nElements; k++)
      fUnitArray[k].SetUnitEnergy(0.);  

  // load input vectors
  for (Int_t i = 0; i < nInT; i++)  
    {
      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
      enT[i]  = lv->Energy();
      ptT[i]  = lv->Pt();
      etaT[i] = lv->Eta();
      phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
      thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
      
      if(fDebug>20){    
	cout << "enT[" << i << "] : " <<  enT[i] << endl;
	cout << "ptT[" << i << "] : " <<  ptT[i] << endl;
	cout << "etaT[" << i << "] : " <<  etaT[i] << endl;
	cout << "phiT[" << i << "] : " <<  phiT[i] << endl;
	cout << "thetaT[" << i << "] : " <<  thetaT[i] << endl;
	cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
	cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
	cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
	cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
      }
      
      if(fGrid==0)
	{
	  // For the moment, only TPC filled from its grid in its total acceptance
	  if(fDebug>2)
	    cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
	    
	    trackInTpcAccOnly += 1;
	    
	    Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
	    
	    Float_t unitEnergy = 0.;
	    unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
	    
	    if(unitEnergy > 0. && fDebug >10){
Commit	Line	Data
4a01bb2c	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	//======================================================================
	18	// *** November 2006
	19	// Author: magali.estienne@ires.in2p3.fr
	20	// 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and
	21	// EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance
	22	// 2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment)
	23	// fGrid==1, work with a part of the TPC acceptance
	24	// 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them
	25	// 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry
	26	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
	27	//======================================================================
	28	// ***September 2006
	29	// TTask : Fill Unit Array for the Tracks information
	30	// Called by ESD reader for jet analysis
	31	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
	32	//======================================================================
	33	// *** July 2006
	34	// 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex
	35	// and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the
	36	// index or the eta, phi position of a grid.
	37	// 2) Define a grid for TPC
	38	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
	39	//======================================================================
	40	// ***July 2006
	41	// Fill Unit Array class
	42	// Class used by AliJetESDReader to fill a UnitArray from the information extracted
	43	// from the particle tracks
	44	// Author: magali.estienne@ires.in2p3.fr
	45	//======================================================================
	46
	47
	48	// --- Standard library ---
	49	#include <Riostream.h>
	50
	51	// --- ROOT system ---
	52	#include <TSystem.h>
	53	#include <TLorentzVector.h>
	54	#include <TVector3.h>
	55	//#include "Math/Vector3D.h"
	56	//#include "Math/Vector3Dfwd.h"
	57	#include "TTask.h"
	58	#include <TGeoManager.h>
	59	#include <TMatrixD.h>
	60	#include <TArrayD.h>
	61
	62	// --- AliRoot header files ---
	63	#include "AliJetFinder.h"
	64	#include "AliJetReaderHeader.h"
65	#include "AliJetReader.h"
66	#include "AliJetESDReader.h"
67	#include "AliJetESDReaderHeader.h"
68	#include "AliESD.h"
69	#include "AliEMCALGeometry.h"
70	#include "AliJetUnitArray.h"
71	#include "AliJetFillUnitArrayTracks.h"
72	#include "AliJetGrid.h"
73
74	ClassImp(AliJetFillUnitArrayTracks)
75
76	AliEMCALGeometry *AliJetFillUnitArrayTracks::fGeom=0;
77
78	//_____________________________________________________________________________
79	AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
80	: TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information")
81	{
82	// constructor
83	fNIn = 0;
84	fOpt = 0;
85	fDebug = 0;
86	fNphi = 0;
87	fNeta = 0;
88	fPhiMin = 0;
89	fPhiMax = 0;
90	fEtaMin = 0;
91	fEtaMax = 0;
92	fEtaBinInTPCAcc = 0;
93	fPhiBinInTPCAcc = 0;
94	fEtaBinInEMCalAcc = 0;
95	fPhiBinInEMCalAcc = 0;
96	fNbinPhi = 0;
97	}
98
99	//____________________________________________________________________________
100	void AliJetFillUnitArrayTracks::SetEMCALGeometry()
101	{
102	// Set EMCAL geometry information
103	fGeom = AliEMCALGeometry::GetInstance();
104	if (fGeom == 0)
105	fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
106	if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
107	}
108
109	//____________________________________________________________________________
110	void AliJetFillUnitArrayTracks::InitParameters()
111	{
112	fHCorrection = 0; // For hadron correction
113	fHadCorr = 0; // For hadron correction
114	fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
115	fDebug = fReaderHeader->GetDebug();
116
117	fEtaMinCal = fGeom->GetArm1EtaMin();
118	fEtaMaxCal = fGeom->GetArm1EtaMax();
119	fPhiMinCal = fGeom->GetArm1PhiMin();
120	fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
121
122	if(fDebug>30){
123	cout << "fEtaMinCal : " << fEtaMinCal << endl;
124	cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
125	cout << "fPhiMinCal : " << fPhiMinCal << endl;
126	cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
127	}
128
129	fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
130	fPhiMax,fEtaMin,fEtaMax);
131	fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
132	fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
133
134	fEta = fTPCGrid->GetArrayEta();
135	fPhi = fTPCGrid->GetArrayPhi();
136	fIndex = fTPCGrid->GetIndexObject();
137
138	if(fDebug>20){
139	for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
140	for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
141
142	for(Int_t i=0; i<fNphi+1; i++)
143	for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
144	(*fIndex)(i,j) << endl; }
145	}
146	if(fDebug>1) printf("\n Parameters initiated ! \n");
147	}
148
149	//_____________________________________________________________________________
150	AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
151	{
152	// destructor
153	}
154
155	//_____________________________________________________________________________
156	void AliJetFillUnitArrayTracks::Exec(Option_t* option)
157	{
158	//
159	// Main method.
160	// Explain
161
162	fDebug = fReaderHeader->GetDebug();
163	if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");
164	if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
165	// Set EMCal Geometry
166	SetEMCALGeometry();
167	// Set parameters
168	InitParameters();
169
170	TClonesArray *lvArray = fMomentumArray; // Correct checked !
171	Int_t nInT = lvArray->GetEntries(); // Correct checked !
172	Float_t ptMin = fReaderHeader->GetPtCut(); // Correct checked !
173
174	// sflag -> not yet implemented !!!!
175
176	if(fDebug>3) cout << "nInT : " << nInT << endl;
177
178	if (nInT == 0) return;
179
180	// local arrays for input
181	Float_t* enT = new Float_t[nInT];
182	Float_t* ptT = new Float_t[nInT];
183	Float_t* etaT = new Float_t[nInT];
184	Float_t* phiT = new Float_t[nInT];
185	Float_t* thetaT = new Float_t[nInT];
186
187	Int_t trackInEmcalAcc = 0;
188	Int_t trackInTpcAcc = 0;
189	Int_t trackInTpcAccOnly = 0;
190
191	Int_t nElements = fTPCGrid->GetNEntries();
192	Int_t nElements2 = fEMCalGrid->GetNEntries();
193	fGrid = fTPCGrid->GetGridType();
194
195	if(fDebug>3){
196	cout << "nElements : " << nElements << endl;
197	cout << "nElements2 : " << nElements2 << endl;
198	cout << "fNumUnits : " << fNumUnits << endl;
199	cout << "sum : " << fNumUnits+nElements << endl;
200	}
201
202	// Set energy exactly to zero
203	if(fGrid==0)
204	for(Int_t k=0; k<nElements; k++)
205	fUnitArray[k].SetUnitEnergy(0.);
206
207	if(fGrid==1)
208	for(Int_t k=0; k<fNumUnits+nElements; k++)
209	fUnitArray[k].SetUnitEnergy(0.);
210
211	// load input vectors
212	for (Int_t i = 0; i < nInT; i++)
213	{
214	TLorentzVector lv = (TLorentzVector) lvArray->At(i);
215	enT[i] = lv->Energy();
216	ptT[i] = lv->Pt();
217	etaT[i] = lv->Eta();
218	phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
219	thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
220
221	if(fDebug>20){
222	cout << "enT[" << i << "] : " << enT[i] << endl;
223	cout << "ptT[" << i << "] : " << ptT[i] << endl;
224	cout << "etaT[" << i << "] : " << etaT[i] << endl;
225	cout << "phiT[" << i << "] : " << phiT[i] << endl;
226	cout << "thetaT[" << i << "] : " << thetaT[i] << endl;
227	cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
228	cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
229	cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
230	cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
231	}
232
233	if(fGrid==0)
234	{
235	// For the moment, only TPC filled from its grid in its total acceptance
236	if(fDebug>2)
237	cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
238
239	trackInTpcAccOnly += 1;
240
241	Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
242
243	Float_t unitEnergy = 0.;
244	unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
245
246	if(unitEnergy > 0. && fDebug >10){
247