-
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* 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
+// Class used by AliJetESDReader to fill a UnitArray from the information
+// extracted from the particle tracks
// Author: magali.estienne@ires.in2p3.fr
-//======================================================================
+//----------------------------------------------------------------------
// --- Standard library ---
// --- ROOT system ---
#include <TSystem.h>
#include <TLorentzVector.h>
+#include <TRefArray.h>
#include <TVector3.h>
-//#include "Math/Vector3D.h"
-//#include "Math/Vector3Dfwd.h"
#include "TTask.h"
#include <TGeoManager.h>
#include <TMatrixD.h>
#include <TArrayD.h>
+#include <TMath.h>
+#include <TClonesArray.h>
// --- AliRoot header files ---
#include "AliJetFinder.h"
#include "AliJetReader.h"
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
-#include "AliESD.h"
-#include "AliEMCALGeometry.h"
+//#include "AliESD.h"
+#include "AliESDEvent.h"
+#include "AliJetDummyGeo.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayTracks.h"
+#include "AliJetHadronCorrectionv1.h"
#include "AliJetGrid.h"
ClassImp(AliJetFillUnitArrayTracks)
-AliEMCALGeometry *AliJetFillUnitArrayTracks::fGeom=0;
-
//_____________________________________________________________________________
AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
- : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information")
+ : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
+ fNumUnits(0),
+ fEtaMinCal(0),
+ fEtaMaxCal(0),
+ fPhiMinCal(0),
+ fPhiMaxCal(0),
+ fHadCorr(0),
+ fHCorrection(0),
+ fNTracks(0),
+ fNTracksCut(0),
+ fOpt(0),
+ fDZ(0),
+ fDebug(0),
+ fReaderHeader(0x0),
+ fMomentumArray(0x0),
+ fUnitArray(0x0),
+ fRefArray(0x0),
+ fTPCGrid(0x0),
+ fEMCalGrid(0x0),
+ fGeom(0x0),
+ fESD(0x0),
+ fGrid0(0x0),
+ fGrid1(0x0),
+ fGrid2(0x0),
+ fGrid3(0x0),
+ fGrid4(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
- 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()
+//_____________________________________________________________________________
+AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
+ : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
+ fNumUnits(0),
+ fEtaMinCal(0),
+ fEtaMaxCal(0),
+ fPhiMinCal(0),
+ fPhiMaxCal(0),
+ fHadCorr(0),
+ fHCorrection(0),
+ fNTracks(0),
+ fNTracksCut(0),
+ fOpt(0),
+ fDZ(0),
+ fDebug(0),
+ fReaderHeader(0x0),
+ fMomentumArray(0x0),
+ fUnitArray(0x0),
+ fRefArray(0x0),
+ fTPCGrid(0x0),
+ fEMCalGrid(0x0),
+ fGeom(0x0),
+ fESD(0x0),
+ fGrid0(0x0),
+ fGrid1(0x0),
+ fGrid2(0x0),
+ fGrid3(0x0),
+ fGrid4(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)
{
- // 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");
+ // constructor
}
//____________________________________________________________________________
void AliJetFillUnitArrayTracks::InitParameters()
{
- fHCorrection = 0; // For hadron correction
+ // fHCorrection = 0; // For hadron correction
fHadCorr = 0; // For hadron correction
fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
fDebug = fReaderHeader->GetDebug();
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,
}
//_____________________________________________________________________________
-void AliJetFillUnitArrayTracks::Exec(Option_t* option)
+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();
+ // get number of tracks in event (for the loop)
+ Int_t goodTrack = 0;
+ Int_t nt = 0;
+ Float_t pt, eta,phi;
+ TVector3 p3;
+ nt = fESD->GetNumberOfTracks();
+ if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
+
+ // temporary storage of signal and pt cut flag
+ Int_t* sflag = new Int_t[nt];
+ Int_t* cflag = new Int_t[nt];
+
+ // get cuts set by user
+ Float_t ptMin = fReaderHeader->GetPtCut();
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
+ fOpt = fReaderHeader->GetDetector();
+ fDZ = fReaderHeader->GetDZ();
- if(fDebug>3){
- cout << "nElements : " << nElements << endl;
- cout << "nElements2 : " << nElements2 << endl;
- cout << "fNumUnits : " << fNumUnits << endl;
- cout << "sum : " << fNumUnits+nElements << endl;
- }
+ Int_t nTracksEmcal = 0;
+ Int_t nTracksEmcalDZ = 0;
+ Int_t nTracksTpc = 0;
+ Int_t nTracksTpcOnly = 0;
+ Int_t nTracksEmcalCut = 0;
+ Int_t nTracksEmcalDZCut = 0;
+ Int_t nTracksTpcCut = 0;
+ Int_t nTracksTpcOnlyCut = 0;
- // 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){
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "TPC id : " << idTPC << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "phiT[" << i << "] : " << phiT[i] << endl;
- cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl;
- cout << "fUnitArray[idTPC].GetUnitEnergy(): " <<
- fUnitArray[idTPC].GetUnitEnergy() << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
+ fGrid = fTPCGrid->GetGridType();
- // Fill energy in TPC acceptance
- fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]);
- if(fDebug > 10){
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "unitEnergy in TPC acceptance after : " <<
- fUnitArray[idTPC].GetUnitEnergy() << endl;
- }
-
- // Pt cut flag
- if(fUnitArray[idTPC].GetUnitEnergy()<ptMin)
- fUnitArray[idTPC].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[idTPC].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[idTPC].SetUnitDetectorFlag(kTpc);
+ //loop over tracks
+ for (Int_t it = 0; it < nt; it++) {
+ AliESDtrack *track = fESD->GetTrack(it);
+ UInt_t status = track->GetStatus();
+
+ Double_t mom[3];
+ track->GetPxPyPz(mom);
+ p3.SetXYZ(mom[0],mom[1],mom[2]);
+ pt = p3.Pt();
+ Float_t mass = 0.;
+ mass = track->GetMass();
+
+ if (((status & AliESDtrack::kITSrefit) == 0) ||
+ ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check
+ if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
+ && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
+ if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
+ && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
+ eta = p3.Eta();
+ phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
+
+ if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
+
+ // sflag -> not yet implemented !!!!
+
+ if(fGrid==0)
+ {
+ // Only TPC filled from its grid in its total acceptance
+
+ Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+ Bool_t ok = kFALSE;
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+ uArray->SetUnitTrackID(it);
+
+ Float_t unitEnergy = 0.;
+ unitEnergy = uArray->GetUnitEnergy();
+ if(unitEnergy==0.){
+ nTracksTpcOnly++;
+ ok = kTRUE;
+ }
+ // Fill energy in TPC acceptance
+ uArray->SetUnitEnergy(unitEnergy + pt);
+ uArray->SetUnitPxPyPz(mom);
+ uArray->SetUnitMass(mass);
+
+ // Pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin){
+ uArray->SetUnitCutFlag(kPtSmaller);
+ }
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok) nTracksTpcOnlyCut++;
}
- if(fGrid==1)
- {
- // Fill track information in EMCAL acceptance
- if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) &&
- (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// &&
- // GetCutFlag(i) == 1)
- // ptT[i] > ptMin)
- {
- trackInEmcalAcc += 1;
-
- if(fDebug>20){
- cout << "before : " << endl;
- cout << "etaT[i] : " << etaT[i] << endl;
- cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
- }
-
- // This function should be modified soon
-// Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete
- Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified
-// Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]); // Using an EMCal grid -> calculated
-// Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]); // Using an EMCal grid -> tabulated (faster)
-
- Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy();
-
- if(fDebug>20) {
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "after : " << endl;
- cout << "towerID : " << towerID << endl;
- cout << "etaT[i] : " << etaT[i] << endl;
- cout << "phiT[i](rad) : " << phiT[i] << endl;
- cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
- cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl;
- cout << "fHadCorr : " << fHadCorr << endl;
- cout << "fHCorrection : " << fHCorrection << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
-
- //OLD WAY: //Do Hadron Correction
- if (fHCorrection != 0)
- {
- Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
- unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]);
- if(fDebug>20){
- cout << "Inside loop for hadron correction ==========================" << endl;
- cout << "enT[" << i << "] : " << enT[i] << endl;
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "Energy correction : " << hCEnergy << endl;
- cout << "unitEnergy : " << unitEnergy << endl;
- cout << "fUnitArray[towerID].GetUnitEnergy() : " <<
- fUnitArray[towerID].GetUnitEnergy() << endl;
- }
- } //end Hadron Correction loop
-
- fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]);
-
- // Put a pt cut flag
- if(fUnitArray[towerID].GetUnitEnergy()<ptMin)
- fUnitArray[towerID].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[towerID].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[towerID].SetUnitDetectorFlag(kTpc);
-
- if(fDebug>10){
- cout << "After pT filled ===============================================" << endl;
- cout << "PtCut : " << ptMin << endl;
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "unitEnergy : " << unitEnergy << endl;
- cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl;
- cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl;
- cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl;
- cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl;
- }
- } // end loop on EMCal acceptance cut + tracks quality
- else{
- // Outside EMCal acceptance
- if(fDebug>2)
- cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
-
- trackInTpcAcc += 1;
-
- // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
- Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
-
- Float_t unitEnergy2 = 0.;
- unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
-
- if(unitEnergy2 > 0. && fDebug >10){
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "TPC id : " << idTPC << endl;
- cout << "Total id : " << fNumUnits-1+idTPC << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "phiT[" << i << "] : " << phiT[i] << endl;
- cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl;
- cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " <<
- fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
+ // Detector flag
+ if(unitEnergy>0) {
+ uArray->SetUnitDetectorFlag(kAll);
+ }
+ if(uArray->GetUnitEnergy()>0){
+ fRefArray->Add(uArray);
+ }
- // Fill energy outside emcal acceptance
- fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
- // Pt cut flag
- if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()<ptMin)
- fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[fNumUnits-1+idTPC].SetUnitDetectorFlag(kTpc);
- }
- } // end fGrid==1
- } // end loop on entries (tpc tracks)
-
- if(fDebug>3){
- printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc);
- printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc);
- }
-
- if(fGrid==0) fNIn = trackInTpcAccOnly;
- if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc;
-
- fOpt = fReaderHeader->GetDetector();
-
- if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt);
-
- if(fOpt==1 || option=="tpc") // if only TPC
- { //Set all unit flags, Eta, Phi
-
- if(fGrid==0)
- {
- if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
- for(Int_t i=0; i<nElements; i++)
- {
- Float_t eta = 0.;
- Float_t phi = 0.;
- if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
- fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
- fUnitArray[i].SetUnitEntries(nElements);
- fUnitArray[i].SetUnitID(i);
- fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta);
- fUnitArray[i].SetUnitEta(eta);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
+ }
+
+ if(fGrid==1)
+ {
+ Int_t nElements = fTPCGrid->GetNEntries();
+
+ // Fill track information in EMCAL acceptance
+ if((eta >= fEtaMin && eta <= fEtaMax) &&
+ (phi >= fPhiMin && phi <= fPhiMax))// &&
+ {
+
+ // Include dead-zones
+ if(fDZ)
+ {
+ Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
+ Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
+ fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
+ fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+ fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+ fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+ fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+ Int_t n0 = fGrid0->GetNEntries();
+ Int_t n1 = fGrid1->GetNEntries();
+ Int_t n2 = fGrid2->GetNEntries();
+ Int_t n3 = fGrid3->GetNEntries();
+
+ if(phi >= phimin0 && phi <= phimax0){
+ Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
+ AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+ uArray0->SetUnitTrackID(it);
+ Float_t uEnergy0 = uArray0->GetUnitEnergy();
+ Bool_t ok0 = kFALSE;
+ if(uEnergy0==0.){
+ nTracksEmcalDZ++;
+ ok0 = kTRUE;
+ }
+ uArray0->SetUnitEnergy(uEnergy0+pt);
+ if(uArray0->GetUnitEnergy()<ptMin)
+ uArray0->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray0->SetUnitCutFlag(kPtHigher);
+ if(ok0) nTracksEmcalDZCut++;
+ }
+ if(uArray0->GetUnitEnergy()>0)
+ fRefArray->Add(uArray0);
}
- }
-
- if(fGrid==1)
- {
- if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
- for(Int_t i=0; i<fNumUnits+nElements; i++)
- {
- Float_t eta = 0.;
- Float_t phi = 0.;
- if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
- //Set all units to be outside a jet initially
- fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
- fUnitArray[i].SetUnitEntries(fNumUnits+nElements);
- fUnitArray[i].SetUnitID(i);
-
- if(fUnitArray[i].GetUnitID()<fNumUnits)
- {
- // fGeom->EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63
- fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry
- phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
- // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid
- // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid
- fUnitArray[i].SetUnitEta(eta);
- //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi());
- }
+ if(phi >= phimin1 && phi <= phimax1){
+ Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
+ AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
+ uArray1->SetUnitTrackID(it);
+ Float_t uEnergy1 = uArray1->GetUnitEnergy();
+ Bool_t ok1 = kFALSE;
+ if(uEnergy1==0.){
+ nTracksEmcalDZ++;
+ ok1 = kTRUE;
+ }
+ uArray1->SetUnitEnergy(uEnergy1+pt);
+ if(uArray1->GetUnitEnergy()<ptMin)
+ uArray1->SetUnitCutFlag(kPtSmaller);
else {
- fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta);
- fUnitArray[i].SetUnitEta(eta);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
-
- if(fDebug>10)
- {
- if(fUnitArray[i].GetUnitEnergy()!=0.){
- cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " <<
- fUnitArray[i].GetUnitID()+1-fNumUnits << endl;
- cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " <<
- fUnitArray[i].GetUnitEnergy() << endl;
- cout << "(fUnitArray[" << i << "].GetUnitEta() : " <<
- fUnitArray[i].GetUnitEta() << endl;
- cout << "(fUnitArray[" << i << "].GetUnitPhi() : " <<
- fUnitArray[i].GetUnitPhi() << endl;
- }
- }
+ uArray1->SetUnitCutFlag(kPtHigher);
+ if(ok1) nTracksEmcalDZCut++;
}
- fUnitArray[i].SetUnitClusterID(0);
- }//end loop over all units in array (same as all towers in EMCAL)
- }
-
- if(fDebug>20)
- {
- for(Int_t i=0; i<fNumUnits+nElements; i++)
- {
- if(fUnitArray[i].GetUnitEnergy()!=0.){
- cout << "######################################################### " << endl;
- cout << "Final UnitArray filled with energy != 0" << i << endl;
- cout << "Pointeur UnitArray : " << fUnitArray << " ID : " <<
- fUnitArray[i].GetUnitID() << " Energy : " << fUnitArray[i].GetUnitEnergy() <<
- " Eta : " << fUnitArray[i].GetUnitEta() << " Phi : " << fUnitArray[i].GetUnitPhi() << endl;
+ if(uArray1->GetUnitEnergy()>0)
+ fRefArray->Add(uArray1);
+ }
+ if(phi >= phimin2 && phi <= phimax2){
+ Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
+ AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
+ uArray2->SetUnitTrackID(it);
+ Float_t uEnergy2 = uArray2->GetUnitEnergy();
+ Bool_t ok2 = kFALSE;
+ if(uEnergy2==0.){
+ nTracksEmcalDZ++;
+ ok2 = kTRUE;
+ }
+ uArray2->SetUnitEnergy(uEnergy2+pt);
+ if(uArray2->GetUnitEnergy()<ptMin)
+ uArray2->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray2->SetUnitCutFlag(kPtHigher);
+ if(ok2) nTracksEmcalDZCut++;
+ }
+ if(uArray2->GetUnitEnergy()>0)
+ fRefArray->Add(uArray2);
+ }
+ if(phi >= phimin3 && phi <= phimax3){
+ Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
+ AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
+ uArray3->SetUnitTrackID(it);
+ Float_t uEnergy3 = uArray3->GetUnitEnergy();
+ Bool_t ok3 = kFALSE;
+ if(uEnergy3==0.){
+ nTracksEmcalDZ++;
+ ok3 = kTRUE;
+ }
+ uArray3->SetUnitEnergy(uEnergy3+pt);
+ if(uArray3->GetUnitEnergy()<ptMin)
+ uArray3->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray3->SetUnitCutFlag(kPtHigher);
+ if(ok3) nTracksEmcalDZCut++;
+ }
+ if(uArray3->GetUnitEnergy()>0)
+ fRefArray->Add(uArray3);
+ }
+ if(phi >= phimin4 && phi <= phimax4){
+ Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
+ AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
+ uArray4->SetUnitTrackID(it);
+ Float_t uEnergy4 = uArray4->GetUnitEnergy();
+ Bool_t ok4 = kFALSE;
+ if(uEnergy4==0.){
+ nTracksEmcalDZ++;
+ ok4 = kTRUE;
+ }
+ uArray4->SetUnitEnergy(uEnergy4+pt);
+ if(uArray4->GetUnitEnergy()<ptMin)
+ uArray4->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray4->SetUnitCutFlag(kPtHigher);
+ if(ok4) nTracksEmcalDZCut++;
}
+ if(uArray4->GetUnitEnergy()>0)
+ fRefArray->Add(uArray4);
}
+ } // end fDZ
+
+ Int_t towerID = 0;
+ fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
+
+ if(towerID==-1) continue;
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
+ uArray->SetUnitTrackID(it);
+ Float_t unitEnergy = uArray->GetUnitEnergy();
+ Bool_t ok = kFALSE;
+ if(unitEnergy==0.){
+ nTracksEmcal++;
+ ok=kTRUE;
+ }
+
+ // Do Hadron Correction
+ // Parametrization to be added
+ if (fHCorrection != 0)
+ {
+ // Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
+ Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
+ unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
+
+ } //end Hadron Correction loop
+
+ uArray->SetUnitEnergy(unitEnergy + pt);
+
+ // Put a pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin)
+ uArray->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok) nTracksEmcalCut++;
}
+ // Detector flag
+ if(unitEnergy > 0)
+ uArray->SetUnitDetectorFlag(kAll);
+
+ if(uArray->GetUnitEnergy()>0)
+ fRefArray->Add(uArray);
+
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
+
+ } // end loop on EMCal acceptance cut + tracks quality
+ else{
+ // Outside EMCal acceptance
+
+ // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
+ Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
+ uArray->SetUnitTrackID(it);
+
+ Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
+ Bool_t ok2 = kFALSE;
+ if(unitEnergy2==0.){
+ nTracksTpc++;
+ ok2=kTRUE;
+ }
+ // Fill energy outside emcal acceptance
+ uArray->SetUnitEnergy(unitEnergy2 + pt);
+
+ // Pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin){
+ uArray->SetUnitCutFlag(kPtSmaller);
+ }
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok2) nTracksTpcCut++;
+ }
+ // Detector flag
+ if(unitEnergy2 > 0)
+ uArray->SetUnitDetectorFlag(kTpc);
+ if(uArray->GetUnitEnergy()>0)
+ fRefArray->Add(uArray);
+
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
+
+ }
+ } // end fGrid==1
+ } // end loop on entries (tpc tracks)
- } // end if(fOpt==1 || option=="tpc")
+// // set the signal flags
+// fSignalFlag.Set(goodTrack,sflag);
+// fCutFlag.Set(goodTrack,cflag);
- delete[] enT;
- delete[] ptT;
- delete[] etaT;
- delete[] phiT;
- delete[] thetaT;
+// delete sflag;
+// delete cflag;
+
+ // } // end loop on entries (tpc tracks)
+ if(fGrid==0) {
+ fNTracks = nTracksTpcOnly;
+ fNTracksCut = nTracksTpcOnlyCut;
+ if(fDebug>10){
+ cout << "fNTracks : " << fNTracks << endl;
+ cout << "fNTracksCut : " << fNTracksCut << endl;
+ }
+ }
+ if(fGrid==1) {
+ fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
+ fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
+ if(fDebug>10){
+ cout << "fNTracks : " << fNTracks << endl;
+ cout << "fNTracksCut : " << fNTracksCut << endl;
+ }
+ }
+
}
//__________________________________________________________
void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
{
+//
+// Obtain the index from eta and phi
+//
for(Int_t j=0; j<fNphi+1; j++) {
for(Int_t i=0; i<fNeta+1; i++) {
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff