[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 "AliJetDummyGeo.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayTracks.h"
#include "AliJetHadronCorrectionv1.h"
#include "AliJetGrid.h"

ClassImp(AliJetFillUnitArrayTracks)

//_____________________________________________________________________________
AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
    : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
      fNumUnits(0),
      fEtaMinCal(0.),
      fEtaMaxCal(0.),
      fPhiMinCal(0.),
      fPhiMaxCal(0.),
      fHadCorr(0x0),
      fHCorrection(0),
      fNIn(0),
      fOpt(0),
      fDebug(0),
      fReaderHeader(0x0),
      fMomentumArray(0x0),
      fUnitArray(0x0),
      fTPCGrid(0x0),
      fEMCalGrid(0x0),
      fGeom(0x0),
      fNphi(0),
      fNeta(0),
      fPhi2(0x0),
      fEta2(0x0),
      fPhi(0x0),
      fEta(0x0),
      fIndex(0x0),
      fParams(0x0),
      fGrid(0),
      fPhiMin(0.),
      fPhiMax(0.),
      fEtaMin(0.),
      fEtaMax(0.),
      fEtaBinInTPCAcc(0),
      fPhiBinInTPCAcc(0),
      fEtaBinInEMCalAcc(0),
      fPhiBinInEMCalAcc(0),
      fNbinPhi(0)
{
  // constructor
}

//____________________________________________________________________________
void AliJetFillUnitArrayTracks::SetEMCALGeometry()
{
    // Set EMCAL geometry information
    fGeom = AliJetDummyGeo::GetInstance();
    if (fGeom == 0)
	fGeom = AliJetDummyGeo::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"
b92e2ccf	69	//#include "AliEMCALGeometry.h"
b92e2ccf	70	#include "AliJetDummyGeo.h"
4a01bb2c	71	#include "AliJetUnitArray.h"
4a01bb2c	72	#include "AliJetFillUnitArrayTracks.h"
7aec0427	73	#include "AliJetHadronCorrectionv1.h"
4a01bb2c	74	#include "AliJetGrid.h"
	75
	76	ClassImp(AliJetFillUnitArrayTracks)
	77
4a01bb2c	78	//_____________________________________________________________________________
4a01bb2c	79	AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
b92e2ccf	80	: TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
	81	fNumUnits(0),
	82	fEtaMinCal(0.),
	83	fEtaMaxCal(0.),
	84	fPhiMinCal(0.),
	85	fPhiMaxCal(0.),
	86	fHadCorr(0x0),
	87	fHCorrection(0),
	88	fNIn(0),
	89	fOpt(0),
	90	fDebug(0),
	91	fReaderHeader(0x0),
	92	fMomentumArray(0x0),
	93	fUnitArray(0x0),
	94	fTPCGrid(0x0),
	95	fEMCalGrid(0x0),
	96	fGeom(0x0),
	97	fNphi(0),
	98	fNeta(0),
	99	fPhi2(0x0),
	100	fEta2(0x0),
	101	fPhi(0x0),
	102	fEta(0x0),
	103	fIndex(0x0),
	104	fParams(0x0),
	105	fGrid(0),
	106	fPhiMin(0.),
	107	fPhiMax(0.),
	108	fEtaMin(0.),
	109	fEtaMax(0.),
	110	fEtaBinInTPCAcc(0),
	111	fPhiBinInTPCAcc(0),
	112	fEtaBinInEMCalAcc(0),
	113	fPhiBinInEMCalAcc(0),
	114	fNbinPhi(0)
4a01bb2c	115	{
4a01bb2c	116	// constructor
4a01bb2c	117	}
	118
	119	//____________________________________________________________________________
	120	void AliJetFillUnitArrayTracks::SetEMCALGeometry()
	121	{
b45b0c92	122	// Set EMCAL geometry information
b92e2ccf	123	fGeom = AliJetDummyGeo::GetInstance();
b45b0c92	124	if (fGeom == 0)
b92e2ccf	125	fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
b45b0c92	126	if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
4a01bb2c	127	}
	128
	129	//____________________________________________________________________________
	130	void AliJetFillUnitArrayTracks::InitParameters()
	131	{
	132	fHCorrection = 0; // For hadron correction
	133	fHadCorr = 0; // For hadron correction
	134	fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
	135	fDebug = fReaderHeader->GetDebug();
	136
	137	fEtaMinCal = fGeom->GetArm1EtaMin();
	138	fEtaMaxCal = fGeom->GetArm1EtaMax();
	139	fPhiMinCal = fGeom->GetArm1PhiMin();
	140	fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
	141
	142	if(fDebug>30){
	143	cout << "fEtaMinCal : " << fEtaMinCal << endl;
	144	cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
	145	cout << "fPhiMinCal : " << fPhiMinCal << endl;
	146	cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
	147	}
	148
	149	fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
	150	fPhiMax,fEtaMin,fEtaMax);
	151	fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
	152	fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
	153
	154	fEta = fTPCGrid->GetArrayEta();
	155	fPhi = fTPCGrid->GetArrayPhi();
	156	fIndex = fTPCGrid->GetIndexObject();
	157
	158	if(fDebug>20){
	159	for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
	160	for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
	161
	162	for(Int_t i=0; i<fNphi+1; i++)
	163	for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
	164	(*fIndex)(i,j) << endl; }
	165	}
	166	if(fDebug>1) printf("\n Parameters initiated ! \n");
	167	}
	168
	169	//_____________________________________________________________________________
	170	AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
	171	{
	172	// destructor
	173	}
	174
	175	//_____________________________________________________________________________
	176	void AliJetFillUnitArrayTracks::Exec(Option_t* option)
	177	{
	178	//
	179	// Main method.
	180	// Explain
	181
	182	fDebug = fReaderHeader->GetDebug();
	183	if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");
	184	if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
	185	// Set EMCal Geometry
	186	SetEMCALGeometry();
	187	// Set parameters
	188	InitParameters();
	189
	190	TClonesArray *lvArray = fMomentumArray; // Correct checked !
191	Int_t nInT = lvArray->GetEntries(); // Correct checked !
192	Float_t ptMin = fReaderHeader->GetPtCut(); // Correct checked !
193
194	// sflag -> not yet implemented !!!!
195
196	if(fDebug>3) cout << "nInT : " << nInT << endl;
197
198	if (nInT == 0) return;
199
200	// local arrays for input
201	Float_t* enT = new Float_t[nInT];
202	Float_t* ptT = new Float_t[nInT];
203	Float_t* etaT = new Float_t[nInT];
204	Float_t* phiT = new Float_t[nInT];
205	Float_t* thetaT = new Float_t[nInT];
206
207	Int_t trackInEmcalAcc = 0;
208	Int_t trackInTpcAcc = 0;
209	Int_t trackInTpcAccOnly = 0;
210
211	Int_t nElements = fTPCGrid->GetNEntries();
212	Int_t nElements2 = fEMCalGrid->GetNEntries();
213	fGrid = fTPCGrid->GetGridType();
214
215	if(fDebug>3){
216	cout << "nElements : " << nElements << endl;
217	cout << "nElements2 : " << nElements2 << endl;
218	cout << "fNumUnits : " << fNumUnits << endl;
219	cout << "sum : " << fNumUnits+nElements << endl;
220	}
221
222	// Set energy exactly to zero
223	if(fGrid==0)
224	for(Int_t k=0; k<nElements; k++)
225	fUnitArray[k].SetUnitEnergy(0.);
226
227	if(fGrid==1)
228	for(Int_t k=0; k<fNumUnits+nElements; k++)
229	fUnitArray[k].SetUnitEnergy(0.);
230
231	// load input vectors
232	for (Int_t i = 0; i < nInT; i++)
233	{
234	TLorentzVector lv = (TLorentzVector) lvArray->At(i);
235	enT[i] = lv->Energy();
236	ptT[i] = lv->Pt();
237	etaT[i] = lv->Eta();
238	phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
239	thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
240
241	if(fDebug>20){
242	cout << "enT[" << i << "] : " << enT[i] << endl;
243	cout << "ptT[" << i << "] : " << ptT[i] << endl;
244	cout << "etaT[" << i << "] : " << etaT[i] << endl;
245	cout << "phiT[" << i << "] : " << phiT[i] << endl;
246	cout << "thetaT[" << i << "] : " << thetaT[i] << endl;
247	cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
248	cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
249	cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
250	cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
251	}
252
253	if(fGrid==0)
254	{
255	// For the moment, only TPC filled from its grid in its total acceptance
256	if(fDebug>2)
257	cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
258
259	trackInTpcAccOnly += 1;
260
261	Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
262
263	Float_t unitEnergy = 0.;
264	unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
265
266	if(unitEnergy > 0. && fDebug >10){
267