From be6e5811083a745dcdf28ea035aa5ded8961ea72 Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 30 Jul 2009 09:10:18 +0000 Subject: [PATCH] Updates provided by Magali Estienne 1 - remove the coding convention violations you sent to me few weeks ago 2 - look for jets using charged + neutral information from AODs --- JETAN/AliJet.cxx | 37 +- JETAN/AliJet.h | 17 +- JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx | 390 ++++++++++++ JETAN/AliJetAODFillUnitArrayEMCalDigits.h | 76 +++ JETAN/AliJetAODFillUnitArrayTracks.cxx | 616 +++++++++++++++++++ JETAN/AliJetAODFillUnitArrayTracks.h | 64 ++ JETAN/AliJetAODReader.cxx | 411 ++++++++++++- JETAN/AliJetAODReader.h | 58 +- JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx | 415 +++++++++++++ JETAN/AliJetESDFillUnitArrayEMCalDigits.h | 77 +++ JETAN/AliJetESDFillUnitArrayTracks.cxx | 624 ++++++++++++++++++++ JETAN/AliJetESDFillUnitArrayTracks.h | 64 ++ JETAN/AliJetESDReader.cxx | 92 ++- JETAN/AliJetESDReader.h | 20 +- JETAN/AliJetFillUnitArray.cxx | 143 +++++ JETAN/AliJetFillUnitArray.h | 122 ++++ JETAN/AliJetFinder.cxx | 11 +- JETAN/AliJetFinder.h | 14 +- JETAN/AliJetReader.cxx | 14 +- JETAN/AliJetReader.h | 25 +- JETAN/AliJetUnitArray.cxx | 171 +++--- JETAN/AliJetUnitArray.h | 75 ++- JETAN/AliUA1JetFinderV2.cxx | 202 +++---- JETAN/AliUA1JetFinderV2.h | 44 +- JETAN/JETANLinkDef.h | 7 +- JETAN/libJETAN.pkg | 4 +- 26 files changed, 3431 insertions(+), 362 deletions(-) create mode 100755 JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx create mode 100755 JETAN/AliJetAODFillUnitArrayEMCalDigits.h create mode 100644 JETAN/AliJetAODFillUnitArrayTracks.cxx create mode 100644 JETAN/AliJetAODFillUnitArrayTracks.h create mode 100755 JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx create mode 100755 JETAN/AliJetESDFillUnitArrayEMCalDigits.h create mode 100644 JETAN/AliJetESDFillUnitArrayTracks.cxx create mode 100644 JETAN/AliJetESDFillUnitArrayTracks.h create mode 100644 JETAN/AliJetFillUnitArray.cxx create mode 100644 JETAN/AliJetFillUnitArray.h diff --git a/JETAN/AliJet.cxx b/JETAN/AliJet.cxx index 23e4e68739e..4367b669758 100644 --- a/JETAN/AliJet.cxx +++ b/JETAN/AliJet.cxx @@ -23,7 +23,6 @@ #include -//#include #include #include @@ -44,12 +43,8 @@ AliJet::AliJet(): fPtIn(0), fPtChPtCutIn(0), fEnTotChPtCutIn(0), - fVectorSizeIn(0), fDetIn(0), - fVPx(0), - fVPy(0), - fVPz(0) - // fVectorIn(0) + fTrackRef(new TRefArray()) { // Default constructor fJets = new TClonesArray("TLorentzVector",1000); @@ -62,7 +57,6 @@ AliJet::AliJet(): fNCells = TArrayI(); fPtChPtCutIn = TArrayF(); fEnTotChPtCutIn = TArrayF(); - fVectorSizeIn = TArrayI(); fDetIn = TArrayI(); } @@ -280,34 +274,6 @@ void AliJet::SetDetectorFlagIn(Int_t* x) //////////////////////////////////////////////////////////////////////// -void AliJet::SetVectorSizeIn(Int_t* x) -{ - if (fNInput>0) fVectorSizeIn.Set(fNInput, x); -} - -//////////////////////////////////////////////////////////////////////// - -void AliJet::SetVectorPxIn(vector< vector > pxT) -{ - fVPx = pxT;; -} - -//////////////////////////////////////////////////////////////////////// - -void AliJet::SetVectorPyIn(vector< vector > pyT) -{ - fVPy = pyT;; -} - -//////////////////////////////////////////////////////////////////////// - -void AliJet::SetVectorPzIn(vector< vector > pzT) -{ - fVPz = pzT;; -} - -//////////////////////////////////////////////////////////////////////// - void AliJet::SetPtFromSignal(Float_t* p) { // set information of percentage of pt of jets @@ -348,7 +314,6 @@ void AliJet::ClearJets(Option_t *option) fNCells.Set(0); fPtChPtCutIn.Set(0); fEnTotChPtCutIn.Set(0); - fVectorSizeIn.Set(0); fDetIn.Set(0); } diff --git a/JETAN/AliJet.h b/JETAN/AliJet.h index 74689e95edf..2d05825d0ae 100644 --- a/JETAN/AliJet.h +++ b/JETAN/AliJet.h @@ -12,11 +12,11 @@ //--------------------------------------------------------------------- #include -#include #include #include #include +#include class TClonesArray; class TLorentzVector; @@ -41,10 +41,7 @@ class AliJet : public TObject TArrayF GetPtIn() const { return fPtIn; } TArrayF GetPtChargedPtCutIn() const { return fPtChPtCutIn; } TArrayF GetEnTotChargedPtCutIn() const {return fEnTotChPtCutIn; } - TArrayI GetVectorSizeIn() const { return fVectorSizeIn; } - vector< vector > GetVectorPxIn() const { return fVPx; } - vector< vector > GetVectorPyIn() const { return fVPy; } - vector< vector > GetVectorPzIn() const { return fVPz; } + TRefArray* GetTrackRef() const { return fTrackRef;} TArrayI GetDetectorFlagIn() const { return fDetIn; } Double_t GetEtAvg() const { return fEtAvg; } @@ -74,10 +71,7 @@ class AliJet : public TObject void SetInJet(Int_t* idx); void SetPtChargedPtCutIn(Float_t* pt2T); void SetEnTotChargedPtCutIn(Float_t* en2T); - void SetVectorSizeIn(Int_t* vectT); - void SetVectorPxIn(vector< vector > pxT); - void SetVectorPyIn(vector< vector > pyT); - void SetVectorPzIn(vector< vector > pzT); + void SetTrackReferences(TRefArray* ref) {fTrackRef = ref;} void SetDetectorFlagIn(Int_t* detT); void SetEtAvg(Double_t et) { fEtAvg = et; } // others @@ -105,11 +99,8 @@ class AliJet : public TObject TArrayF fPtIn; // Arrays of input particles kine:Pt TArrayF fPtChPtCutIn; // Arrays of input particles kin:Pt Charged with pt cut TArrayF fEnTotChPtCutIn; // Arrays of total energy with pt cut on charged + cut min on cell - TArrayI fVectorSizeIn; // Arrays of number of charged tracks in each unitArray TArrayI fDetIn; // Arrays of detector type of each UnitArray - vector< vector > fVPx; //|| - vector< vector > fVPy; //|| - vector< vector > fVPz; //|| + TRefArray* fTrackRef; //|| Reference to tracks which could belongs to the jet ClassDef(AliJet,2) }; diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx new file mode 100755 index 00000000000..2e9429d913b --- /dev/null +++ b/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx @@ -0,0 +1,390 @@ + +/************************************************************************** + * 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. * + **************************************************************************/ + +//====================================================================== +// Fill Unit Array +// Called by AOD reader for jet analysis +// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) +//====================================================================== + +class TSystem; +class TLorentzVector; +class TGeoManager; +class TArrayS; + +// --- AliRoot header files --- +#include "AliJetUnitArray.h" +#include "AliJetAODFillUnitArrayEMCalDigits.h" +class AliJetFinder; + +// #include "AliEMCALCalibData.h" +// #include "AliCDBManager.h" +// class AliCDBStorage; +// #include "AliCDBEntry.h" + + +ClassImp(AliJetAODFillUnitArrayEMCalDigits) + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits() + : AliJetFillUnitArray(), + fAOD(0), + fNIn(0), + fCluster(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fApplyElectronCorrection(kFALSE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod) + : AliJetFillUnitArray(), + fAOD(aod), + fNIn(0), + fCluster(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fApplyElectronCorrection(kFALSE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det) + : AliJetFillUnitArray(det), + fAOD(det.fAOD), + fNIn(det.fNIn), + fCluster(det.fCluster), + fNCEMCAL(det.fNCEMCAL), + fNCPHOS(det.fNCPHOS), + fNCCalo(det.fNCCalo), + fApplyElectronCorrection(det.fApplyElectronCorrection), + fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection), + fFractionHadronicCorrection(det.fFractionHadronicCorrection), + fClus(det.fClus), + fNDigitEmcal(det.fNDigitEmcal), + fNDigitEmcalCut(det.fNDigitEmcalCut) +{ + // Copy constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayEMCalDigits& AliJetAODFillUnitArrayEMCalDigits::operator=(const AliJetAODFillUnitArrayEMCalDigits& other) +{ + // Assignment + + fAOD = other.fAOD; + fNIn = other.fNIn; + fCluster = other.fCluster; + fNCEMCAL = other.fNCEMCAL; + fNCPHOS = other.fNCPHOS; + fNCCalo = other.fNCCalo; + fApplyElectronCorrection = other.fApplyElectronCorrection; + fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection; + fFractionHadronicCorrection = other.fFractionHadronicCorrection; + fClus = other.fClus; + fNDigitEmcal = other.fNDigitEmcal; + fNDigitEmcalCut = other.fNDigitEmcalCut; + + return (*this); + +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits() +{ + // destructor +} + + + +//_____________________________________________________________________________ +void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) +//void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* option) +{ + // + // Main method. + // Fill the unit array with the neutral particle information from the EMCal cells in AOD + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + fCluster = fReaderHeader->GetCluster(); + + // Init parameters + InitParameters(); + +//(not used ?) Int_t nDigitTot = 0; + Int_t goodDigit = 0; +//(not used ?) Int_t beg = 0; +//(not used ?) Int_t end = 0; + Int_t index = 0; +//(not used ?) Int_t count = 0; + +//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); + + if(!fCluster) { // Keep all digit information + // Loop over all cell information + //------------------------------------------------------------------ + AliAODCaloCells &cells= *(fAOD->GetEMCALCells()); + Int_t ncell = cells.GetNumberOfCells() ; +//(not used ?) Int_t type = cells.GetType(); + + for (Int_t icell= 0; icell < ncell; icell++) { + Int_t digitID = cells.GetCellNumber(icell); + Float_t digitAmp = cells.GetAmplitude(icell); + Float_t digitEn = digitAmp*0.0153; // Last correct + // Float_t digitEn = Calibrate(digitAmp,digitID); + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); + + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); + + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + uArray->SetUnitEnergy(unitEnergy+digitEt); + + uArray->SetUnitCutFlag(kPtHigher); + + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 12) printf("goodDigit : %d\n", goodDigit); + + } // End loop over cells + } // end if !fCluster + else { // Keep digit information from clusterization + + // Loop over calo clusters + //------------------------------------------------------------------ + + // select EMCAL clusters only + TRefArray * caloClusters = new TRefArray(); + fAOD->GetEMCALClusters(caloClusters); + + // Total number of EMCAL cluster + Int_t nclus = caloClusters->GetEntries() ; + Int_t beg = 0; + Float_t pos[3] ; + + // Get CaloCells + AliAODCaloCells &cells= *(fAOD->GetEMCALCells()); +//(not used ?) Int_t ncell = cells.GetNumberOfCells() ; +//(not used ?) Int_t type = cells.GetType(); + + for(Int_t j = beg; j < nclus; j++) { // loop over clusters + // Retrieve cluster from aod + AliAODCaloCluster *fClus = (AliAODCaloCluster *) caloClusters->At(j) ; + + // Get the cluster info +//(not used ?) Float_t energy = fClus->E() ; + + fClus->GetPosition(pos) ; + TVector3 vpos(pos[0],pos[1],pos[2]) ; + + Int_t digMult = fClus->GetNCells() ; + UShort_t *digID = fClus->GetCellsAbsId() ; +//(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ; + Int_t trackIndex = -1; + if(fClus->GetNTracksMatched()!=0) + trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID(); +//(not used ?) Int_t cLabel = fClus->GetLabel(0); + // Do double-counted electron correction + if (fApplyElectronCorrection != 0 && trackIndex !=-1 ) + { + // The electron correction go there + // Under construction !!!! + } // End of Electron correction + + // Get CaloCells of cluster and fill the unitArray + for(Int_t i = 0; i < digMult ; i++) { + Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ; +//(not used ?) Double_t digitAmpFrac = digAmpFrac[i]; + Float_t digitAmp = cells.GetCellAmplitude(digitID) ; + + // Calibration for an energy in GeV + Float_t digitEn = digitAmp*0.0153; + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); + + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + uArray->SetUnitEnergy(unitEnergy+digitEt); + + uArray->SetUnitCutFlag(kPtHigher); + + // Signal or background + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 12) printf("goodDigit : %d\n", goodDigit); + + }// End loop over cells + } // End loop over clusters + } // end else + + fNIn += goodDigit; + + if(fDebug>1) + { + printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"); + printf("goodDigit : %d\n", goodDigit); + } +} + +/* +//____________________________________________________________________________ +void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters() +{ + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, they are taken from digitizer. + // + // It is a user responsilibity to open CDB before reconstruction, + // for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet())) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + + if(!fCalibData) + printf("************* Calibration parameters not found in CDB! ****************"); +// AliFatal("Calibration parameters not found in CDB!"); + + +} + +//____________________________________________________________________________ +Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId) +{ + + // Convert digitized amplitude into energy. + // Calibration parameters are taken from calibration data base for raw data, + // or from digitizer parameters for simulated data. + + if(fCalibData){ + + if (fGeom==0) + printf("************* Did not get geometry from EMCALLoader ***************"); +// AliFatal("Did not get geometry from EMCALLoader") ; + + Int_t iSupMod = -1; + Int_t nModule = -1; + Int_t nIphi = -1; + Int_t nIeta = -1; + Int_t iphi = -1; + Int_t ieta = -1; + + Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ; + if(!bCell) { + // fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } + + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + + fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); + fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + + return -fADCpedestalECA + amp * fADCchannelECA ; + + } + else //Return energy with default parameters if calibration is not available + return -fADCpedestalECA + amp * fADCchannelECA ; + +} +*/ + + diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.h b/JETAN/AliJetAODFillUnitArrayEMCalDigits.h new file mode 100755 index 00000000000..eff0287165e --- /dev/null +++ b/JETAN/AliJetAODFillUnitArrayEMCalDigits.h @@ -0,0 +1,76 @@ +#ifndef ALIJETAODFILLUNITARRAYEMCALDIGITS_H +#define ALIJETAODFILLUNITARRAYEMCALDIGITS_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Jet Fill Unit Array +// Called by ESD Reader for jet analysis +// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr) +//--------------------------------------------------------------------- + +#include "AliJetFillUnitArray.h" + +class AliESDCaloCluster; +class AliJetReader; +class AliJetAODReader; + +class AliJetAODFillUnitArrayEMCalDigits : public AliJetFillUnitArray +{ + public: + AliJetAODFillUnitArrayEMCalDigits(); + AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *fAOD); + virtual ~AliJetAODFillUnitArrayEMCalDigits(); + + // Setter + void SetApplyElectronCorrection(Int_t flag = 1) {fApplyElectronCorrection = flag;} + void SetApplyFractionHadronicCorrection(Bool_t val) {fApplyFractionHadronicCorrection = val;} + void SetFractionHadronicCorrection(Double_t val) {fFractionHadronicCorrection = val;} + void SetAOD(AliAODEvent* const aod) {fAOD = aod;} + void SetInitMult(Int_t mult) {fNDigitEmcal = mult;} + void SetInitMultCut(Int_t multcut) {fNDigitEmcalCut = multcut;} + + // Getter + Int_t GetMult() const {return fNDigitEmcal;} + Int_t GetMultCut() const {return fNDigitEmcalCut;} + + // For calibration + // virtual Float_t Calibrate(Int_t amp, Int_t cellId) ; // Tranforms Amp to energy + + // Other + void Exec(Option_t* const option); + + protected: + AliAODEvent *fAOD; // ESD + Int_t fNIn; // Number of Array filled in UnitArray + Int_t fCluster; // Use all cells or cells in clusters for jet finding + Int_t fNCEMCAL; // Number of clusters in EMCAL + Int_t fNCPHOS; // Number of clusters in PHOS + Int_t fNCCalo; // Number of cluster in EMCAL + PHOS calorimeters + + Bool_t fApplyElectronCorrection; // Electron correction flag + Bool_t fApplyFractionHadronicCorrection; // Fraction hadronic correction flag + Bool_t fFractionHadronicCorrection; // Fraction hadronic correction + + + // geometry info + AliESDCaloCluster *fClus; //! + Int_t fNDigitEmcal; //! + Int_t fNDigitEmcalCut; //! + //Calibration parameters... to be replaced by database +// AliEMCALCalibData *fCalibData; //! Calibration database if aval +// Float_t fADCchannelECA; // width of one ADC channel for EC section (GeV) +// Float_t fADCpedestalECA;// pedestal of ADC for EC section (GeV) + + + private: + AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det); + AliJetAODFillUnitArrayEMCalDigits &operator=(const AliJetAODFillUnitArrayEMCalDigits &det); + +// void GetCalibrationParameters(void) ; + + ClassDef(AliJetAODFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information +}; + +#endif diff --git a/JETAN/AliJetAODFillUnitArrayTracks.cxx b/JETAN/AliJetAODFillUnitArrayTracks.cxx new file mode 100644 index 00000000000..59da52ce936 --- /dev/null +++ b/JETAN/AliJetAODFillUnitArrayTracks.cxx @@ -0,0 +1,616 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +//====================================================================== +// ***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 +//====================================================================== + + +// --- ROOT system --- +#include +#include + +#include "AliJetUnitArray.h" +#include "AliJetAODFillUnitArrayTracks.h" + +// --- ROOT system --- +class TSystem; +class TLorentzVector; +class TGeoManager; + +// --- AliRoot header files --- +class AliJetFinder; + +ClassImp(AliJetAODFillUnitArrayTracks) + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks() + : AliJetFillUnitArray(), + fNumUnits(0), + fHadCorr(0), + fApplyMIPCorrection(kTRUE), + fAOD(0x0), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod) + : AliJetFillUnitArray(), + fNumUnits(0), + fHadCorr(0), + fApplyMIPCorrection(kTRUE), + fAOD(aod), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det) + : AliJetFillUnitArray(det), + fNumUnits(det.fNumUnits), + fHadCorr(det.fHadCorr), + fApplyMIPCorrection(det.fApplyMIPCorrection), + fAOD(det.fAOD), + fGrid0(det.fGrid0), + fGrid1(det.fGrid1), + fGrid2(det.fGrid2), + fGrid3(det.fGrid3), + fGrid4(det.fGrid4) +{ + // Copy constructor +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other) +{ + // Assignment + + fNumUnits = other.fNumUnits; + fHadCorr = other.fHadCorr; + fApplyMIPCorrection = other.fApplyMIPCorrection; + fAOD = other.fAOD; + fGrid0 = other.fGrid0; + fGrid1 = other.fGrid1; + fGrid2 = other.fGrid2; + fGrid3 = other.fGrid3; + fGrid4 = other.fGrid4; + + return (*this); +} + +//____________________________________________________________________________ +void AliJetAODFillUnitArrayTracks::InitParameters() +{ + fHadCorr = 0; // For hadron correction + fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL + + fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, + fPhiMax,fEtaMin,fEtaMax); + fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, + fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi); + + fIndex = fTPCGrid->GetIndexObject(); + + if(fDebug>20){ + for(Int_t i=0; i1) printf("\n Parameters initiated ! \n"); +} + +//_____________________________________________________________________________ +AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/) +//void AliJetAODFillUnitArrayTracks::Exec(Option_t* option) +{ + // + // Main method. + // Fill the unit array with the charged particle information in AOD + // + + fDebug = fReaderHeader->GetDebug(); + + // Set parameters + InitParameters(); + + // get number of tracks in event (for the loop) + Int_t goodTrack = 0; + Int_t nt = 0; +//(not used ?) Int_t nt2 = 0; + Int_t nmax = 0; + Float_t pt,eta,phi; + // Int_t sflag = 0; + TVector3 p3; + + nt = fAOD->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(); + UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask(); + + 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; + +//(not used ?) Int_t nElements = fTPCGrid->GetNEntries(); +//(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries(); + fGrid = fTPCGrid->GetGridType(); + +//(not used ?) Int_t nEntries = fUnitArray->GetEntries(); + +//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); + + //loop over tracks + nmax = nt; + for (Int_t it = 0; it < nmax; it++) { + AliAODTrack *track; + track = fAOD->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(); + +//(not used ?) Float_t mass = 0.; + + eta = p3.Eta(); + phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi()); + if (status == 0) continue; + if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue; + + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + + if(fGrid==0) + { + // Only TPC filled from its grid in its total acceptance + + Int_t idTPC = fTPCGrid->GetIndex(phi,eta); + Bool_t ok = kFALSE; + Bool_t ref = kFALSE; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy = 0.; + unitEnergy = uArray->GetUnitEnergy(); + // nTracksTpcOnly is necessary to count the number of candidate cells + // but it doesn't give the real multiplicity -> it will be extracted + // in the jet finder through the reference to tracks + if(unitEnergy==0.){ + nTracksTpcOnly++; + ok = kTRUE; + ref = kTRUE; + } + + // Fill energy in TPC acceptance + uArray->SetUnitEnergy(unitEnergy + pt); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksTpcOnlyCut++; + } + + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId) { + fProcId = kTRUE; + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + } + fRefArray->Add(uArray); + } + } + + 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(); +//(not used ?) Int_t n4 = fGrid4->GetNEntries(); + + // cout << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl; + if(phi >= phimin0 && phi <= phimax0){ + Int_t id0 = fGrid0->GetIndex(phi,eta)-1; + AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements); + TRefArray *reference0 = uArray0->GetUnitTrackRef(); + reference0->Add(track); + uArray0->SetUnitTrackID(it); + + Float_t uEnergy0 = uArray0->GetUnitEnergy(); + Bool_t ok0 = kFALSE; + Bool_t ref0 = kFALSE; + if(uEnergy0==0.){ + nTracksEmcalDZ++; + ok0 = kTRUE; + ref0 = kTRUE; + } + uArray0->SetUnitEnergy(uEnergy0+pt); + + if(uArray0->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray0->SetUnitCutFlag(kPtHigher); + if(ok0) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray0->SetUnitSignalFlag(kGood); + uArray0->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray0->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray0->GetUnitEnergy()>0 && ref0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0)); + fProcId = kTRUE; + } + fRefArray->Add(uArray0); + } + } + + if(phi >= phimin1 && phi <= phimax1){ + Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0; + AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements); + TRefArray *reference1 = uArray1->GetUnitTrackRef(); + reference1->Add(track); + + Float_t uEnergy1 = uArray1->GetUnitEnergy(); + Bool_t ok1 = kFALSE; + Bool_t ref1 = kFALSE; + if(uEnergy1==0.){ + nTracksEmcalDZ++; + ok1 = kTRUE; + ref1 = kTRUE; + } + uArray1->SetUnitEnergy(uEnergy1+pt); + + if(uArray1->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray1->SetUnitCutFlag(kPtHigher); + if(ok1) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray1->SetUnitSignalFlag(kGood); + uArray1->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray1->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray1->GetUnitEnergy()>0 && ref1){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1)); + fProcId = kTRUE; + } + 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); + TRefArray *reference2 = uArray2->GetUnitTrackRef(); + reference2->Add(track); + + Float_t uEnergy2 = uArray2->GetUnitEnergy(); + Bool_t ok2 = kFALSE; + Bool_t ref2 = kFALSE; + if(uEnergy2==0.){ + nTracksEmcalDZ++; + ok2 = kTRUE; + ref2 = kTRUE; + } + uArray2->SetUnitEnergy(uEnergy2+pt); + + if(uArray2->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray2->SetUnitCutFlag(kPtHigher); + if(ok2) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray2->SetUnitSignalFlag(kGood); + uArray2->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray2->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray2->GetUnitEnergy()>0 && ref2){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2)); + fProcId = kTRUE; + } + 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); + TRefArray *reference3 = uArray3->GetUnitTrackRef(); + reference3->Add(track); + + Float_t uEnergy3 = uArray3->GetUnitEnergy(); + Bool_t ok3 = kFALSE; + Bool_t ref3 = kFALSE; + if(uEnergy3==0.){ + nTracksEmcalDZ++; + ok3 = kTRUE; + ref3 = kTRUE; + } + uArray3->SetUnitEnergy(uEnergy3+pt); + + if(uArray3->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray3->SetUnitCutFlag(kPtHigher); + if(ok3) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray3->SetUnitSignalFlag(kGood); + uArray3->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray3->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray3->GetUnitEnergy()>0 && ref3){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3)); + fProcId = kTRUE; + } + 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); + TRefArray *reference4 = uArray4->GetUnitTrackRef(); + reference4->Add(track); + + Float_t uEnergy4 = uArray4->GetUnitEnergy(); + Bool_t ok4 = kFALSE; + Bool_t ref4 = kFALSE; + if(uEnergy4==0.){ + nTracksEmcalDZ++; + ok4 = kTRUE; + ref4 = kTRUE; + } + uArray4->SetUnitEnergy(uEnergy4+pt); + + if(uArray4->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray4->SetUnitCutFlag(kPtHigher); + if(ok4) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray4->SetUnitSignalFlag(kGood); + uArray4->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray4->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray4->GetUnitEnergy()>0 && ref4){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4)); + fProcId = kTRUE; + } + fRefArray->Add(uArray4); + } + } + } // end fDZ + + Int_t towerID = 0; + fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID); + if(towerID==-1) continue; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy = uArray->GetUnitEnergy(); + Bool_t ok = kFALSE; + Bool_t ref = kFALSE; + if(unitEnergy==0.){ + nTracksEmcal++; + ok=kTRUE; + ref=kTRUE; + } + + // Do Hadron Correction + // This is under construction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + // Parametrization to be added + if (fApplyMIPCorrection != 0) + { +// Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); +// 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()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksEmcalCut++; + } + + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } + } // 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); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1 + Bool_t okout = kFALSE; + Bool_t refout = kFALSE; + if(unitEnergy2==0.){ + nTracksTpc++; + okout=kTRUE; + refout=kTRUE; + } + // Fill energy outside emcal acceptance + uArray->SetUnitEnergy(unitEnergy2 + pt); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(okout) nTracksTpcCut++; + } + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray->GetUnitEnergy()>0 && refout){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } + } + } // end fGrid==1 + + goodTrack++; + + } // end loop on entries (tpc tracks) + + + if(fDebug>0) + { + cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl; + cout << "goodTracks: " << goodTrack << endl; + } + + 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; + } + } + +} + + + + + + + + + + + + + diff --git a/JETAN/AliJetAODFillUnitArrayTracks.h b/JETAN/AliJetAODFillUnitArrayTracks.h new file mode 100644 index 00000000000..67dd8fa0eb2 --- /dev/null +++ b/JETAN/AliJetAODFillUnitArrayTracks.h @@ -0,0 +1,64 @@ +#ifndef ALIJETAODFILLUNITARRAYTRACKS_H +#define ALIJETAODFILLUNITARRAYTRACKS_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Jet Fill Unit Array +// Called by ESD Reader for jet analysis +// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr) +//--------------------------------------------------------------------- + +#include "AliJetFillUnitArray.h" + +class AliJetReader; +class AliJetAODReader; +class AliEMCALGeometry; + +class AliJetAODFillUnitArrayTracks : public AliJetFillUnitArray +{ + public: + AliJetAODFillUnitArrayTracks(); + AliJetAODFillUnitArrayTracks(AliAODEvent *fAOD); + virtual ~AliJetAODFillUnitArrayTracks(); + + // Setter + void SetHadCorrector(AliJetHadronCorrection* const corr) {fHadCorr = corr;} + void SetApplyMIPCorrection(Bool_t const val) {fApplyMIPCorrection = val;} + void SetAOD(AliAODEvent* const aod) {fAOD = aod;} + void SetGrid0(AliJetGrid* const grid0) {fGrid0 = grid0;} + void SetGrid1(AliJetGrid* const grid1) {fGrid1 = grid1;} + void SetGrid2(AliJetGrid* const grid2) {fGrid2 = grid2;} + void SetGrid3(AliJetGrid* const grid3) {fGrid3 = grid3;} + void SetGrid4(AliJetGrid* const grid4) {fGrid4 = grid4;} + + // Getter + Int_t GetHadCorrection() const {return fApplyMIPCorrection;} + Int_t GetMult() const {return fNTracks;} + Int_t GetMultCut() const {return fNTracksCut;} + + // Other + void Exec(Option_t* const option); + + protected: + Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) + AliJetHadronCorrection *fHadCorr; // Pointer to Hadron Correction Object + Bool_t fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection + + AliAODEvent *fAOD; // ESD + AliJetGrid *fGrid0; // Grid used for dead zones definition + AliJetGrid *fGrid1; // Grid used for dead zones definition + AliJetGrid *fGrid2; // Grid used for dead zones definition + AliJetGrid *fGrid3; // Grid used for dead zones definition + AliJetGrid *fGrid4; // Grid used for dead zones definition + + private: + AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det); + AliJetAODFillUnitArrayTracks &operator=(const AliJetAODFillUnitArrayTracks &det); + void InitParameters(); + + ClassDef(AliJetAODFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information +}; + +#endif diff --git a/JETAN/AliJetAODReader.cxx b/JETAN/AliJetAODReader.cxx index 89c7d622f57..b62ee7fc0d7 100644 --- a/JETAN/AliJetAODReader.cxx +++ b/JETAN/AliJetAODReader.cxx @@ -28,21 +28,57 @@ #include #include #include +#include +#include +#include #include "AliJetAODReader.h" #include "AliJetAODReaderHeader.h" #include "AliAODEvent.h" #include "AliAODTrack.h" +#include "AliJetDummyGeo.h" +#include "AliJetAODFillUnitArrayTracks.h" +#include "AliJetAODFillUnitArrayEMCalDigits.h" +#include "AliJetHadronCorrection.h" +#include "AliJetUnitArray.h" ClassImp(AliJetAODReader) AliJetAODReader::AliJetAODReader(): AliJetReader(), fChain(0x0), + fTree(0x0), fAOD(0x0), fRef(new TRefArray), fDebug(0), - fOpt(0) + fOpt(0), + fGeom(0), + fHadCorr(0x0), + fTpcGrid(0x0), + fEmcalGrid(0x0), + fGrid0(0), + fGrid1(0), + fGrid2(0), + fGrid3(0), + fGrid4(0), + fPtCut(0), + fHCorrection(0), + fApplyElectronCorrection(kFALSE), + fEFlag(kFALSE), + fApplyMIPCorrection(kTRUE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), + fNumUnits(0), + fMass(0), + fSign(0), + fNIn(0), + fDZ(0), + fNeta(0), + fNphi(0), + fArrayInitialised(0), + fRefArray(0x0), + fProcId(kFALSE) + { // Constructor } @@ -55,6 +91,17 @@ AliJetAODReader::~AliJetAODReader() delete fChain; delete fAOD; delete fRef; + delete fTpcGrid; + delete fEmcalGrid; + if(fDZ) + { + delete fGrid0; + delete fGrid1; + delete fGrid2; + delete fGrid3; + delete fGrid4; + } + } //____________________________________________________________________________ @@ -88,13 +135,12 @@ void AliJetAODReader::OpenInputFiles() gSystem->FreeDirectory(dir); - fAOD = 0; fChain->SetBranchAddress("AOD",&fAOD); int nMax = fChain->GetEntries(); - if(fDebug>1)printf("\n AliJetAODReader::OpenInputFiles Total number of events in chain= %d \n",nMax); + printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax); // set number of events in header if (fReaderHeader->GetLastEvent() == -1) @@ -113,7 +159,7 @@ void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) { fChain = (TChain*) tree; Int_t nMax = fChain->GetEntries(); - if(fDebug>1)printf("\n AliJetAODReader::ConnectTree Total number of events in chain= %5d \n", nMax); + printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax); // set number of events in header if (fReaderHeader->GetLastEvent() == -1) fReaderHeader->SetLastEvent(nMax); @@ -138,7 +184,7 @@ Bool_t AliJetAODReader::FillMomentumArray() // get number of tracks in event (for the loop) Int_t nt = fAOD->GetNTracks(); - if(fDebug>2)printf(" AliJetAODReader::FillMomentumArray AOD tracks: %5d \t", nt); + printf("AOD tracks: %5d \t", nt); // temporary storage of signal and pt cut flag Int_t* sflag = new Int_t[nt]; @@ -173,7 +219,7 @@ Bool_t AliJetAODReader::FillMomentumArray() aodTrack++; fRef->Add(track); } - if(fDebug>2)printf("Used AOD tracks: %5d \n", aodTrack); + printf("Used AOD tracks: %5d \n", aodTrack); // set the signal flags fSignalFlag.Set(aodTrack,sflag); fCutFlag.Set(aodTrack,cflag); @@ -183,3 +229,356 @@ Bool_t AliJetAODReader::FillMomentumArray() return kTRUE; } + +//__________________________________________________________ +void AliJetAODReader::SetApplyMIPCorrection(Bool_t val) +{ + // + // Set flag to apply MIP correction fApplyMIPCorrection + // - exclusive with fApplyFractionHadronicCorrection + // + + fApplyMIPCorrection = val; + if(fApplyMIPCorrection == kTRUE) + { + SetApplyFractionHadronicCorrection(kFALSE); + printf("Enabling MIP Correction \n"); + } + else + { + printf("Disabling MIP Correction \n"); + } +} + +//__________________________________________________________ +void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val) +{ + // + // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection + // - exclusive with fApplyMIPCorrection + // + + fApplyFractionHadronicCorrection = val; + if(fApplyFractionHadronicCorrection == kTRUE) + { + SetApplyMIPCorrection(kFALSE); + printf("Enabling Fraction Hadronic Correction \n"); + } + else + { + printf("Disabling Fraction Hadronic Correction \n"); + } +} + +//__________________________________________________________ +void AliJetAODReader::SetFractionHadronicCorrection(Double_t val) +{ + // + // Set value to fFractionHadronicCorrection (default is 0.3) + // apply EMC hadronic correction fApplyFractionHadronicCorrection + // - exclusive with fApplyMIPCorrection + // + + fFractionHadronicCorrection = val; + if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0) + { + SetApplyFractionHadronicCorrection(kTRUE); + printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection); + } + else + { + SetApplyFractionHadronicCorrection(kFALSE); + } +} + +//____________________________________________________________________________ +void AliJetAODReader::CreateTasks(TChain* tree) +{ + fDebug = fReaderHeader->GetDebug(); + fDZ = fReaderHeader->GetDZ(); + fTree = tree; + + // Init EMCAL geometry and create UnitArray object + SetEMCALGeometry(); + // cout << "In create task" << endl; + InitParameters(); + InitUnitArray(); + + fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); + fFillUAFromTracks = new AliJetAODFillUnitArrayTracks(); + fFillUAFromTracks->SetReaderHeader(fReaderHeader); + fFillUAFromTracks->SetGeom(fGeom); + fFillUAFromTracks->SetTPCGrid(fTpcGrid); + fFillUAFromTracks->SetEMCalGrid(fEmcalGrid); + + if(fDZ) + { + fFillUAFromTracks->SetGrid0(fGrid0); + fFillUAFromTracks->SetGrid1(fGrid1); + fFillUAFromTracks->SetGrid2(fGrid2); + fFillUAFromTracks->SetGrid3(fGrid3); + fFillUAFromTracks->SetGrid4(fGrid4); + } + fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection); + fFillUAFromTracks->SetHadCorrector(fHadCorr); + fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits(); + fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); + fFillUAFromEMCalDigits->SetGeom(fGeom); + fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); + fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); + fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection); + fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection); + fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection); + + fFillUnitArray->Add(fFillUAFromTracks); + fFillUnitArray->Add(fFillUAFromEMCalDigits); + fFillUAFromTracks->SetActive(kFALSE); + fFillUAFromEMCalDigits->SetActive(kFALSE); + + cout << "Tasks instantiated at that stage ! " << endl; + cout << "You can loop over events now ! " << endl; + +} + +//____________________________________________________________________________ +Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray) +{ + // + // Main function + // Fill the reader part + // + + fProcId = procid; + fRefArray = refArray; +//(not used ?) Int_t nEntRef = fRefArray->GetEntries(); +//(not used ?) Int_t nEntUnit = fUnitArray->GetEntries(); + + // clear momentum array + ClearArray(); + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + + if(!fAOD) { + return kFALSE; + } + + // TPC only or Digits+TPC or Clusters+TPC + if(fOpt%2==!0 && fOpt!=0){ + fFillUAFromTracks->SetAOD(fAOD); + fFillUAFromTracks->SetActive(kTRUE); + fFillUAFromTracks->SetUnitArray(fUnitArray); + fFillUAFromTracks->SetRefArray(fRefArray); + fFillUAFromTracks->SetProcId(fProcId); + // fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed + fFillUAFromTracks->Exec("tpc"); + if(fOpt==1){ + fNumCandidate = fFillUAFromTracks->GetMult(); + fNumCandidateCut = fFillUAFromTracks->GetMultCut(); + } + } + + // Digits only or Digits+TPC + if(fOpt>=2 && fOpt<=3){ + fFillUAFromEMCalDigits->SetAOD(fAOD); + fFillUAFromEMCalDigits->SetActive(kTRUE); + fFillUAFromEMCalDigits->SetUnitArray(fUnitArray); + fFillUAFromEMCalDigits->SetRefArray(fRefArray); + fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId()); + fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult()); + fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut()); + fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added + fNumCandidate = fFillUAFromEMCalDigits->GetMult(); + fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut(); + } + + // fFillUnitArray->ExecuteTask(); // => Temporarily commented + + return kTRUE; +} + +//____________________________________________________________________________ +Bool_t AliJetAODReader::SetEMCALGeometry() +{ + // + // Set the EMCal Geometry + // + + if (!fTree->GetFile()) + return kFALSE; + + TString geomFile(fTree->GetFile()->GetName()); + geomFile.ReplaceAll("AliESDs", "geometry"); + + // temporary workaround for PROOF bug #18505 + geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root"); + if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data()); + + // Define EMCAL geometry to be able to read ESDs + fGeom = AliJetDummyGeo::GetInstance(); + if (fGeom == 0) + fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL"); + + // To be setted to run some AliEMCALGeometry functions + TGeoManager::Import(geomFile); + fGeom->GetTransformationForSM(); + printf("\n EMCal Geometry set ! \n"); + + return kTRUE; + +} + +//____________________________________________________________________________ +void AliJetAODReader::InitParameters() +{ + // Initialise parameters + fOpt = fReaderHeader->GetDetector(); + // fHCorrection = 0; // For hadron correction + fHadCorr = 0; // For hadron correction + if(fEFlag==kFALSE){ + if(fOpt==0 || fOpt==1) + fECorrection = 0; // For electron correction + else fECorrection = 1; // For electron correction + } + fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL + if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); +} + +//____________________________________________________________________________ +void AliJetAODReader::InitUnitArray() +{ + //Initialises unit arrays + Int_t nElements = fTpcGrid->GetNEntries(); + Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.; + if(fArrayInitialised) fUnitArray->Delete(); + + if(fTpcGrid->GetGridType()==0) + { // Fill the following quantities : + // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi, + // detector flag, in/out jet, pt cut, mass, cluster ID) + for(Int_t nBin = 1; nBin < nElements+1; nBin++) + { + // fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi); + fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta); + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fTpcGrid->GetDeta(); + Dphi = fTpcGrid->GetDphi(); + new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + } + + if(fTpcGrid->GetGridType()==1) + { + Int_t nGaps = 0; + Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0; + + if(fDZ) + { + // Define a grid of cell for the gaps between SM + 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); + fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015 + fGrid0->SetGridType(0); + fGrid0->SetMatrixIndexes(); + fGrid0->SetIndexIJ(); + n0 = fGrid0->GetNEntries(); + fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1); + fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015 + fGrid1->SetGridType(0); + fGrid1->SetMatrixIndexes(); + fGrid1->SetIndexIJ(); + n1 = fGrid1->GetNEntries(); + fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2); + fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015 + fGrid2->SetGridType(0); + fGrid2->SetMatrixIndexes(); + fGrid2->SetIndexIJ(); + n2 = fGrid2->GetNEntries(); + fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3); + fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015 + fGrid3->SetGridType(0); + fGrid3->SetMatrixIndexes(); + fGrid3->SetIndexIJ(); + n3 = fGrid3->GetNEntries(); + fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4); + fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015 + fGrid4->SetGridType(0); + fGrid4->SetMatrixIndexes(); + fGrid4->SetIndexIJ(); + n4 = fGrid4->GetNEntries(); + + nGaps = n0+n1+n2+n3+n4; + + } + + for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) + { + if(nBinEtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry + // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values + Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else { + if(nBin>=fNumUnits && nBinGetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta); + phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); + Deta = fTpcGrid->GetDeta(); + Dphi = fTpcGrid->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else { + if(fDZ) { + if(nBin>=fNumUnits+nElements && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta); + Deta = fGrid0->GetDeta(); + Dphi = fGrid0->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta); + Deta = fGrid1->GetDeta(); + Dphi = fGrid1->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta); + Deta = fGrid2->GetDeta(); + Dphi = fGrid2->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta); + Deta = fGrid3->GetDeta(); + Dphi = fGrid3->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta); + Deta = fGrid4->GetDeta(); + Dphi = fGrid4->GetDphi(); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); + } + } + } // end if(fDZ) + } // end else 2 + } // end else 1 + } // end loop on nBin + } // end grid type == 1 + fArrayInitialised = 1; +} diff --git a/JETAN/AliJetAODReader.h b/JETAN/AliJetAODReader.h index 889f15a3c89..d11bbb2507b 100644 --- a/JETAN/AliJetAODReader.h +++ b/JETAN/AliJetAODReader.h @@ -11,6 +11,11 @@ //--------------------------------------------------------------------- #include "AliJetReader.h" +#include "AliJetUnitArray.h" +#include "AliJetGrid.h" +class AliJetUnitArray; +class AliJetDummyGeo; +class AliJetHadronCorrection; class AliJetAODReaderHeader; class AliJetReaderHeader; class AliAODEvent; @@ -27,17 +32,66 @@ class AliJetAODReader : public AliJetReader Bool_t FillMomentumArray(); void OpenInputFiles(); void ConnectTree(TTree* tree, TObject* data); - void SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* /*mc*/) {fAOD = (AliAODEvent*) aod;} + void InitUnitArray(); + void CreateTasks(TChain* tree); + Bool_t ExecTasks(Bool_t procid, TRefArray* refArray); + + void SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* /*mc*/) {fAOD = (AliAODEvent*) aod;} + void SetTPCGrid(AliJetGrid *grid) {fTpcGrid = grid;} + void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;} + // Correction of hadronic energy + void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} + void SetHadronCorrector(AliJetHadronCorrection* corr) {fHadCorr = corr;} + void SetApplyElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;} + void SetApplyMIPCorrection(Bool_t val); + void SetApplyFractionHadronicCorrection(Bool_t val); + void SetFractionHadronicCorrection(Double_t val); + + private: + Bool_t SetEMCALGeometry(); + void InitParameters(); AliJetAODReader(const AliJetAODReader &det); AliJetAODReader &operator=(const AliJetAODReader &det); private: TChain *fChain; //! chain for reconstructed tracks + TChain *fTree; //! tree to get the EMCal geometry AliAODEvent *fAOD; //! pointer to aod - TRefArray *fRef; //! pointer to array of references to tracks + TRefArray *fRef; // pointer to array of references to tracks Int_t fDebug; // Debug option Int_t fOpt; // Detector to be used for jet reconstruction + AliJetDummyGeo *fGeom; //! EMCAL Geometry + + AliJetHadronCorrection *fHadCorr; //! Pointer to Hadron Correction Object + AliJetGrid *fTpcGrid; //! Pointer to grid object + AliJetGrid *fEmcalGrid; //! Pointer to grid object + AliJetGrid *fGrid0; // Pointer to grid object + AliJetGrid *fGrid1; // Pointer to grid object + AliJetGrid *fGrid2; // Pointer to grid object + AliJetGrid *fGrid3; // Pointer to grid object + AliJetGrid *fGrid4; // Pointer to grid object + Float_t fPtCut; // Pt cut for tracks to minimise background contribution + Int_t fHCorrection; // Hadron correction flag + Int_t fApplyElectronCorrection; // Electron correction flag + Bool_t fEFlag; // if (fEFlag == kFALSE) => fECorrection automatically setted + Bool_t fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection + Bool_t fApplyFractionHadronicCorrection; // Another type of charged particle energy deposition in EMC + Double_t fFractionHadronicCorrection; // Fraction of momentum of the TPC track to be subtracted from EMC tower + Int_t fNumUnits; // Number of units in the unit object array + // (same as num towers in EMCAL) + Float_t fMass; // Particle mass + Int_t fSign; // Particle sign + Int_t fNIn; // Number of Array filled in UnitArray + Bool_t fDZ; // Use or not dead zones + Int_t fNeta; // Number of bins in eta of tpc grid + Int_t fNphi; // Number of bins in phi of tpc grid + Bool_t fArrayInitialised; // To check that array of units is initialised + + TRefArray *fRefArray; // array of digit position and energy + Bool_t fProcId; // Bool_t for TProcessID synchronization + + ClassDef(AliJetAODReader,1) }; diff --git a/JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx new file mode 100755 index 00000000000..199d64d20d6 --- /dev/null +++ b/JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx @@ -0,0 +1,415 @@ + +/************************************************************************** + * 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. * + **************************************************************************/ + +//--------------------------------------------------------------------- +// Fill Unit Array +// Called by ESD reader for jet analysis +// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) +//--------------------------------------------------------------------- + +// --- ROOT system --- +#include +#include + +// --- AliRoot header files --- +#include "AliESDCaloCluster.h" +#include "AliESDCaloCells.h" +#include "AliJetUnitArray.h" +#include "AliJetESDFillUnitArrayEMCalDigits.h" +// #include "AliEMCALCalibData.h" +// #include "AliCDBManager.h" +// // class AliCDBStorage; +// #include "AliCDBEntry.h" + +// --- ROOT system --- +class TSystem; +class TGeoManager; +class TArrayS; + +// --- AliRoot header files --- +class AliJetFinder; +class AliJetReader; +class AliJetESDReader; +// class AliEMCALCalibData; +// class AliCDBManager; +// class AliCDBStorage; +// class AliCDBEntry; + +ClassImp(AliJetESDFillUnitArrayEMCalDigits) + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits() + : AliJetFillUnitArray(), + fESD(0), + fNIn(0), + fOpt(0), + fCluster(0), + fDebug(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fApplyElectronCorrection(kFALSE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd) + : AliJetFillUnitArray(), + fESD(esd), + fNIn(0), + fOpt(0), + fCluster(0), + fDebug(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fApplyElectronCorrection(kFALSE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det) + : AliJetFillUnitArray(det), + fESD(det.fESD), + fNIn(det.fNIn), + fOpt(det.fOpt), + fCluster(det.fCluster), + fDebug(det.fDebug), + fNCEMCAL(det.fNCEMCAL), + fNCPHOS(det.fNCPHOS), + fNCCalo(det.fNCCalo), + fApplyElectronCorrection(det.fApplyElectronCorrection), + fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection), + fFractionHadronicCorrection(det.fFractionHadronicCorrection), + fClus(det.fClus), + fNDigitEmcal(det.fNDigitEmcal), + fNDigitEmcalCut(det.fNDigitEmcalCut) +{ + // Copy constructor +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other) +{ + // Assignment + + fESD = other.fESD; + fNIn = other.fNIn; + fOpt = other.fOpt; + fCluster = other.fCluster; + fDebug = other.fDebug; + fNCEMCAL = other.fNCEMCAL; + fNCPHOS = other.fNCPHOS; + fNCCalo = other.fNCCalo; + fApplyElectronCorrection = other.fApplyElectronCorrection; + fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection; + fFractionHadronicCorrection = other.fFractionHadronicCorrection; + fClus = other.fClus; + fNDigitEmcal = other.fNDigitEmcal; + fNDigitEmcalCut = other.fNDigitEmcalCut; + + return (*this); + +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/) +{ + // + // Main method. + // Fill the unit array with the neutral particle information from the EMCal cells in ESD + // + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + fCluster = fReaderHeader->GetCluster(); + + // Init parameters + InitParameters(); + +//(not used ?) Int_t nDigitTot = 0; + Int_t goodDigit = 0; +//(not used ?) Int_t beg = 0; +//(not used ?) Int_t end = 0; + Int_t index = 0; +//(not used ?) Int_t count = 0; + +//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); + + if(!fCluster) { // Keep all digit information + // Loop over all cell information + //------------------------------------------------------------------ + AliESDCaloCells &cells= *(fESD->GetEMCALCells()); + Int_t ncell = cells.GetNumberOfCells() ; +//(not used ?) Int_t type = cells.GetType(); + + for (Int_t icell= 0; icell < ncell; icell++) { + Int_t digitID = cells.GetCellNumber(icell); + Float_t digitAmp = cells.GetAmplitude(icell); +//(not used ?) Float_t digitTime = cells.GetTime(icell); + Float_t digitEn = digitAmp*0.0153; // Last correct + // Float_t digitEn = Calibrate(digitAmp,digitID); + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); + + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); + + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + uArray->SetUnitEnergy(unitEnergy+digitEt); + + uArray->SetUnitCutFlag(kPtHigher); + + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 1) printf("goodDigit : %d\n", goodDigit); + + } // End loop over cells + } // end if !fCluster + else { // Keep digit information from clusterization + + // Loop over calo clusters + //------------------------------------------------------------------ + + //select EMCAL clusters only + TRefArray * caloClusters = new TRefArray(); + fESD->GetEMCALClusters(caloClusters); + + // Total number of EMCAL cluster + Int_t nclus = caloClusters->GetEntries() ; + Int_t beg = 0; + Float_t pos[3] ; + + // Get reconstructed vertex position + Double_t vertexPosition[3] ; + fESD->GetVertex()->GetXYZ(vertexPosition) ; + + // Get CaloCells + AliESDCaloCells &cells= *(fESD->GetEMCALCells()); +//(not used ?) Int_t ncell = cells.GetNumberOfCells() ; +//(not used ?) Int_t type = cells.GetType(); + + for(Int_t j = beg; j < nclus; j++) { // loop over clusters + // Retrieve cluster from esd + AliESDCaloCluster *fClus = (AliESDCaloCluster *) caloClusters->At(j) ; + + // Get the cluster info +//(not used ?) Float_t energy = fClus->E() ; +//(not used ?) Int_t iprim = fClus->GetLabel(); + + fClus->GetPosition(pos) ; + TVector3 vpos(pos[0],pos[1],pos[2]) ; + TLorentzVector p ; + fClus->GetMomentum(p,vertexPosition); + + Int_t digMult = fClus->GetNCells() ; + UShort_t *digID = fClus->GetCellsAbsId() ; +//(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ; + Int_t trackIndex = fClus->GetTrackMatched(); + + // Do double-counted electron correction + if (fApplyElectronCorrection != 0 && trackIndex !=-1 ) + { + // The electron correction go there + // Under construction !!!! + } // End of Electron correction + + // Get CaloCells of cluster and fill the unitArray + for(Int_t i = 0; i < digMult ; i++){ + Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ; +//(not used ?) Double_t digitAmpFrac = digAmpFrac[i]; + Float_t digitAmp = cells.GetCellAmplitude(digitID) ; +//(not used ?) Float_t digitTime = cells.GetCellTime(digitID); + + // Calibration for an energy in GeV + Float_t digitEn = digitAmp*0.0153; + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); + + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + uArray->SetUnitEnergy(unitEnergy+digitEt); + + uArray->SetUnitCutFlag(kPtHigher); + + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 12) printf("goodDigit : %d\n", goodDigit); + + } // End loop over cells + } // End loop over clusters + } // end else + + fNIn += goodDigit; + + if(fDebug>1) + { + printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"); + printf("goodDigit : %d\n", goodDigit); + } +} + +/* +//____________________________________________________________________________ +void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters() +{ + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, they are taken from digitizer. + // + // It is a user responsilibity to open CDB before reconstruction, + // for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet())) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + + if(!fCalibData) + printf("************* Calibration parameters not found in CDB! ****************"); +// AliFatal("Calibration parameters not found in CDB!"); + + +} + +//____________________________________________________________________________ +Float_t AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId) +{ + + // Convert digitized amplitude into energy. + // Calibration parameters are taken from calibration data base for raw data, + // or from digitizer parameters for simulated data. + + if(fCalibData){ + + if (fGeom==0) + printf("************* Did not get geometry from EMCALLoader ***************"); +// AliFatal("Did not get geometry from EMCALLoader") ; + + Int_t iSupMod = -1; + Int_t nModule = -1; + Int_t nIphi = -1; + Int_t nIeta = -1; + Int_t iphi = -1; + Int_t ieta = -1; + + Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ; + if(!bCell) { + // fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } + + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + + fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); + fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + + return -fADCpedestalECA + amp * fADCchannelECA ; + + } + else //Return energy with default parameters if calibration is not available + return -fADCpedestalECA + amp * fADCchannelECA ; + +} +*/ + diff --git a/JETAN/AliJetESDFillUnitArrayEMCalDigits.h b/JETAN/AliJetESDFillUnitArrayEMCalDigits.h new file mode 100755 index 00000000000..230adba1d16 --- /dev/null +++ b/JETAN/AliJetESDFillUnitArrayEMCalDigits.h @@ -0,0 +1,77 @@ +#ifndef ALIJETESDFILLUNITARRAYEMCALDIGITS_H +#define ALIJETESDFILLUNITARRAYEMCALDIGITS_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Jet Fill Unit Array +// Called by ESD Reader for jet analysis +// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr) +//--------------------------------------------------------------------- + +#include "AliJetFillUnitArray.h" + +class AliESDCaloCluster; +class AliJetReader; +class AliJetESDReader; +class AliJetESDReaderHeader; +//class AliEMCALCalibData ; + +class AliJetESDFillUnitArrayEMCalDigits : public AliJetFillUnitArray +{ + public: + AliJetESDFillUnitArrayEMCalDigits(); + AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *fESD); + virtual ~AliJetESDFillUnitArrayEMCalDigits(); + + // Setter + void SetApplyElectronCorrection(Int_t flag = 1) {fApplyElectronCorrection = flag;} + void SetApplyFractionHadronicCorrection(Bool_t val) {fApplyFractionHadronicCorrection = val;} + void SetFractionHadronicCorrection(Double_t val) {fFractionHadronicCorrection = val;} + void SetESD(AliESDEvent* const esd) {fESD = esd;} + void SetInitMult(Int_t mult) {fNDigitEmcal = mult;} + void SetInitMultCut(Int_t multcut) {fNDigitEmcalCut = multcut;} + + // Getter + Int_t GetMult() const {return fNDigitEmcal;} + Int_t GetMultCut() const {return fNDigitEmcalCut;} + + // Other + void Exec(Option_t* const option); + // For calibration + // virtual Float_t Calibrate(Int_t amp, Int_t cellId) ; // Tranforms Amp to energy + + protected: + AliESDEvent *fESD; // ESD + Int_t fNIn; // Number of Array filled in UnitArray + Int_t fOpt; // Detector to be used for jet reconstruction + Int_t fCluster; // Use all cells or cells in clusters for jet finding + Int_t fDebug; // Debug option + Int_t fNCEMCAL; // Number of clusters in EMCAL + Int_t fNCPHOS; // Number of clusters in PHOS + Int_t fNCCalo; // Number of cluster in EMCAL + PHOS calorimeters + + Bool_t fApplyElectronCorrection; // Electron correction flag + Bool_t fApplyFractionHadronicCorrection; // Fraction hadronic correction flag + Bool_t fFractionHadronicCorrection; // Fraction hadronic correction + + AliESDCaloCluster *fClus; //! + Int_t fNDigitEmcal; //! + Int_t fNDigitEmcalCut; //! + // Calibration parameters... to be replaced by database +/* AliEMCALCalibData *fCalibData; //! Calibration database if aval */ +/* Float_t fADCchannelECA; // width of one ADC channel for EC section (GeV) */ +/* Float_t fADCpedestalECA; // pedestal of ADC for EC section (GeV) */ + + + private: + AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det); + AliJetESDFillUnitArrayEMCalDigits &operator=(const AliJetESDFillUnitArrayEMCalDigits &det); + +// void GetCalibrationParameters(void) ; + + ClassDef(AliJetESDFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information +}; + +#endif diff --git a/JETAN/AliJetESDFillUnitArrayTracks.cxx b/JETAN/AliJetESDFillUnitArrayTracks.cxx new file mode 100644 index 00000000000..837fc333d29 --- /dev/null +++ b/JETAN/AliJetESDFillUnitArrayTracks.cxx @@ -0,0 +1,624 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +//====================================================================== +// ***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 +//====================================================================== + + +// --- ROOT system --- +#include +#include + +// --- AliRoot header files --- +#include "AliJetUnitArray.h" +#include "AliJetESDFillUnitArrayTracks.h" + +// --- ROOT system --- +class TSystem; +class TLorentzVector; +class TGeoManager; + +// --- AliRoot header files --- +class AliJetFinder; +class AliJetReader; +class AliJetESDReader; + +ClassImp(AliJetESDFillUnitArrayTracks) + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks() + : AliJetFillUnitArray(), + fNumUnits(0), + fHadCorr(0), + fApplyMIPCorrection(kTRUE), + fESD(0x0), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(AliESDEvent* esd) + : AliJetFillUnitArray(), + fNumUnits(0), + fHadCorr(0), + fApplyMIPCorrection(kTRUE), + fESD(esd), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det) + : AliJetFillUnitArray(det), + fNumUnits(det.fNumUnits), + fHadCorr(det.fHadCorr), + fApplyMIPCorrection(det.fApplyMIPCorrection), + fESD(det.fESD), + fGrid0(det.fGrid0), + fGrid1(det.fGrid1), + fGrid2(det.fGrid2), + fGrid3(det.fGrid3), + fGrid4(det.fGrid4) +{ + // Copy constructor +} + +//------------------------------------------------------------------ +AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJetESDFillUnitArrayTracks& other) +{ + // Assignment + + fNumUnits = other.fNumUnits; + fHadCorr = other.fHadCorr; + fApplyMIPCorrection = other.fApplyMIPCorrection; + fESD = other.fESD; + fGrid0 = other.fGrid0; + fGrid1 = other.fGrid1; + fGrid2 = other.fGrid2; + fGrid3 = other.fGrid3; + fGrid4 = other.fGrid4; + + return (*this); +} + +//____________________________________________________________________________ +void AliJetESDFillUnitArrayTracks::InitParameters() +{ + fHadCorr = 0; // For hadron correction + fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL + cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl; + + fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, + fPhiMax,fEtaMin,fEtaMax); + fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, + fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi); + + fIndex = fTPCGrid->GetIndexObject(); + + if(fDebug>20){ + for(Int_t i=0; i1) printf("\n Parameters initiated ! \n"); +} + +//_____________________________________________________________________________ +AliJetESDFillUnitArrayTracks::~AliJetESDFillUnitArrayTracks() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/) +{ + // + // Main method. + // Fill the unit array with the charged particle information in ESD + // + + fDebug = fReaderHeader->GetDebug(); + + // Set parameters + InitParameters(); + + // get number of tracks in event (for the loop) + Int_t goodTrack = 0; + Int_t nt = 0; +//(not used ?) Int_t nt2 = 0; + Int_t nmax = 0; + Float_t pt,eta,phi; + // Int_t sflag = 0; + TVector3 p3; + + if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl; + + 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(); + + 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; + +//(not used ?) Int_t nElements = fTPCGrid->GetNEntries(); +//(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries(); + fGrid = fTPCGrid->GetGridType(); + +//(not used ?) Int_t nEntries = fUnitArray->GetEntries(); + +//(not used ?) Int_t nRefEnt = fRefArray->GetEntries(); + + //loop over tracks + nmax = nt; + for (Int_t it = 0; it < nmax; it++) { + AliESDtrack *track; + 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::kTPCrefit) == 0) continue; // quality check + if ((status & AliESDtrack::kITSrefit) == 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[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + + if(fGrid==0) + { + // Only TPC filled from its grid in its total acceptance + + Int_t idTPC = fTPCGrid->GetIndex(phi,eta); + Bool_t ok = kFALSE; + Bool_t ref = kFALSE; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy = 0.; + unitEnergy = uArray->GetUnitEnergy(); + // nTracksTpcOnly is necessary to count the number of candidate cells + // but it doesn't give the real multiplicity -> it will be extracted + // in the jet finder through the reference to tracks + if(unitEnergy==0.){ + nTracksTpcOnly++; + ok = kTRUE; + ref = kTRUE; + } + + // Fill energy in TPC acceptance + uArray->SetUnitEnergy(unitEnergy + pt); + uArray->SetUnitMass(mass); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksTpcOnlyCut++; + } + + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId) { + fProcId = kTRUE; + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + } + fRefArray->Add(uArray); + } + } + + 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(); +//(not used ?) Int_t n4 = fGrid4->GetNEntries(); + + if(phi >= phimin0 && phi <= phimax0){ + Int_t id0 = fGrid0->GetIndex(phi,eta)-1; + AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements); + TRefArray *reference0 = uArray0->GetUnitTrackRef(); + reference0->Add(track); + + Float_t uEnergy0 = uArray0->GetUnitEnergy(); + Bool_t ok0 = kFALSE; + Bool_t ref0 = kFALSE; + if(uEnergy0==0.){ + nTracksEmcalDZ++; + ok0 = kTRUE; + ref0 = kTRUE; + } + uArray0->SetUnitEnergy(uEnergy0+pt); + + if(uArray0->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray0->SetUnitCutFlag(kPtHigher); + if(ok0) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray0->SetUnitSignalFlag(kGood); + uArray0->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray0->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray0->GetUnitEnergy()>0 && ref0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0)); + fProcId = kTRUE; + } + fRefArray->Add(uArray0); + } + } + + if(phi >= phimin1 && phi <= phimax1){ + Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0; + AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements); + TRefArray *reference1 = uArray1->GetUnitTrackRef(); + reference1->Add(track); + + Float_t uEnergy1 = uArray1->GetUnitEnergy(); + Bool_t ok1 = kFALSE; + Bool_t ref1 = kFALSE; + if(uEnergy1==0.){ + nTracksEmcalDZ++; + ok1 = kTRUE; + ref1 = kTRUE; + } + uArray1->SetUnitEnergy(uEnergy1+pt); + + if(uArray1->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray1->SetUnitCutFlag(kPtHigher); + if(ok1) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray1->SetUnitSignalFlag(kGood); + uArray1->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray1->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray1->GetUnitEnergy()>0 && ref1){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1)); + fProcId = kTRUE; + } + 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); + TRefArray *reference2 = uArray2->GetUnitTrackRef(); + reference2->Add(track); + + Float_t uEnergy2 = uArray2->GetUnitEnergy(); + Bool_t ok2 = kFALSE; + Bool_t ref2 = kFALSE; + if(uEnergy2==0.){ + nTracksEmcalDZ++; + ok2 = kTRUE; + ref2 = kTRUE; + } + uArray2->SetUnitEnergy(uEnergy2+pt); + + if(uArray2->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray2->SetUnitCutFlag(kPtHigher); + if(ok2) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray2->SetUnitSignalFlag(kGood); + uArray2->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray2->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray2->GetUnitEnergy()>0 && ref2){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2)); + fProcId = kTRUE; + } + 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); + TRefArray *reference3 = uArray3->GetUnitTrackRef(); + reference3->Add(track); + + Float_t uEnergy3 = uArray3->GetUnitEnergy(); + Bool_t ok3 = kFALSE; + Bool_t ref3 = kFALSE; + if(uEnergy3==0.){ + nTracksEmcalDZ++; + ok3 = kTRUE; + ref3 = kTRUE; + } + uArray3->SetUnitEnergy(uEnergy3+pt); + + if(uArray3->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray3->SetUnitCutFlag(kPtHigher); + if(ok3) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray3->SetUnitSignalFlag(kGood); + uArray3->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray3->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray3->GetUnitEnergy()>0 && ref3){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3)); + fProcId = kTRUE; + } + 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); + TRefArray *reference4 = uArray4->GetUnitTrackRef(); + reference4->Add(track); + + Float_t uEnergy4 = uArray4->GetUnitEnergy(); + Bool_t ok4 = kFALSE; + Bool_t ref4 = kFALSE; + if(uEnergy4==0.){ + nTracksEmcalDZ++; + ok4 = kTRUE; + ref4 = kTRUE; + } + uArray4->SetUnitEnergy(uEnergy4+pt); + + if(uArray4->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray4->SetUnitCutFlag(kPtHigher); + if(ok4) nTracksEmcalDZCut++; + } + if(sflag[goodTrack] == 1) { + uArray4->SetUnitSignalFlag(kGood); + uArray4->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray4->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray4->GetUnitEnergy()>0 && ref4){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4)); + fProcId = kTRUE; + } + fRefArray->Add(uArray4); + } + } + } // end fDZ + + Int_t towerID = 0; + fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID); + if(towerID==-1) continue; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy = uArray->GetUnitEnergy(); + Bool_t ok = kFALSE; + Bool_t ref = kFALSE; + if(unitEnergy==0.){ + nTracksEmcal++; + ok=kTRUE; + ref=kTRUE; + } + + // Do Hadron Correction + // This is under construction !!!!!!!!!!!!!!!!!!!!!!! + // Parametrization to be added + if (fApplyMIPCorrection != 0) + { +// Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0); +// 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()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksEmcalCut++; + } + + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } + } // end loop on EMCal acceptance cut + tracks quality + else{ + // Outside EMCal acceptance + + Int_t idTPC = fTPCGrid->GetIndex(phi,eta); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC); + TRefArray *reference = uArray->GetUnitTrackRef(); + reference->Add(track); + + Float_t unitEnergy2 = uArray->GetUnitEnergy(); + Bool_t okout = kFALSE; + Bool_t refout = kFALSE; + if(unitEnergy2==0.){ + nTracksTpc++; + okout=kTRUE; + refout=kTRUE; + } + // Fill energy outside emcal acceptance + uArray->SetUnitEnergy(unitEnergy2 + pt); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(okout) nTracksTpcCut++; + } + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + + if(uArray->GetUnitEnergy()>0 && refout){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } + } + } // end fGrid==1 + + goodTrack++; + + } // end loop on entries (tpc tracks) + + if(fDebug>0) + { + cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl; + cout << "goodTracks: " << goodTrack << endl; + } + + delete[] sflag; + delete[] cflag; + + 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; + } + } + +} + + + + + + + + + + + + + + diff --git a/JETAN/AliJetESDFillUnitArrayTracks.h b/JETAN/AliJetESDFillUnitArrayTracks.h new file mode 100644 index 00000000000..200773bfe18 --- /dev/null +++ b/JETAN/AliJetESDFillUnitArrayTracks.h @@ -0,0 +1,64 @@ +#ifndef ALIJETESDFILLUNITARRAYTRACKS_H +#define ALIJETESDFILLUNITARRAYTRACKS_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Jet Fill Unit Array +// Called by ESD Reader for jet analysis +// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr) +//--------------------------------------------------------------------- + +#include "AliJetFillUnitArray.h" +#include "AliJetESDReaderHeader.h" + +class AliEMCALGeometry; +class AliJetReader; +class AliJetESDReader; + +class AliJetESDFillUnitArrayTracks : public AliJetFillUnitArray +{ + public: + AliJetESDFillUnitArrayTracks(); + AliJetESDFillUnitArrayTracks(AliESDEvent *fESD); + virtual ~AliJetESDFillUnitArrayTracks(); + + // Setter + void SetHadCorrector(AliJetHadronCorrection* const corr) {fHadCorr = corr;} + void SetApplyMIPCorrection(Bool_t const val) {fApplyMIPCorrection = val;} + void SetESD(AliESDEvent* const esd) {fESD = esd;} + void SetGrid0(AliJetGrid* const grid0) {fGrid0 = grid0;} + void SetGrid1(AliJetGrid* const grid1) {fGrid1 = grid1;} + void SetGrid2(AliJetGrid* const grid2) {fGrid2 = grid2;} + void SetGrid3(AliJetGrid* const grid3) {fGrid3 = grid3;} + void SetGrid4(AliJetGrid* const grid4) {fGrid4 = grid4;} + + // Getter + Int_t GetHadCorrection() const {return fApplyMIPCorrection;} + Int_t GetMult() const {return fNTracks;} + Int_t GetMultCut() const {return fNTracksCut;} + + // Other + void Exec(Option_t* const option); + + protected: + Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) + AliJetHadronCorrection *fHadCorr; // Pointer to Hadron Correction Object + Bool_t fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection + AliESDEvent *fESD; // ESD + AliJetGrid *fGrid0; // Grid used for dead zones definition + AliJetGrid *fGrid1; // Grid used for dead zones definition + AliJetGrid *fGrid2; // Grid used for dead zones definition + AliJetGrid *fGrid3; // Grid used for dead zones definition + AliJetGrid *fGrid4; // Grid used for dead zones definition + + private: + AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det); + AliJetESDFillUnitArrayTracks &operator=(const AliJetESDFillUnitArrayTracks &det); + void InitParameters(); + + ClassDef(AliJetESDFillUnitArrayTracks,1) +}; + +#endif diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 6675a5d5b48..6db1134a341 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -45,8 +45,8 @@ #include "AliESD.h" #include "AliESDtrack.h" #include "AliJetDummyGeo.h" -#include "AliJetFillUnitArrayTracks.h" -#include "AliJetFillUnitArrayEMCalDigits.h" +#include "AliJetESDFillUnitArrayTracks.h" +#include "AliJetESDFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" #include "AliAnalysisTask.h" @@ -67,9 +67,11 @@ AliJetESDReader::AliJetESDReader(): fGrid3(0), fGrid4(0), fPtCut(0), - fHCorrection(0), - fECorrection(0), + fApplyElectronCorrection(kFALSE), fEFlag(kFALSE), + fApplyMIPCorrection(kTRUE), + fApplyFractionHadronicCorrection(kFALSE), + fFractionHadronicCorrection(0.3), fNumUnits(0), fDebug(0), fMass(0), @@ -151,6 +153,68 @@ void AliJetESDReader::OpenInputFiles() } +//__________________________________________________________ +void AliJetESDReader::SetApplyMIPCorrection(Bool_t val) +{ + // + // Set flag to apply MIP correction fApplyMIPCorrection + // - exclusive with fApplyFractionHadronicCorrection + // + + fApplyMIPCorrection = val; + if(fApplyMIPCorrection == kTRUE) + { + SetApplyFractionHadronicCorrection(kFALSE); + printf("Enabling MIP Correction \n"); + } + else + { + printf("Disabling MIP Correction \n"); + } +} + +//__________________________________________________________ +void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val) +{ + // + // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection + // - exclusive with fApplyMIPCorrection + // + + fApplyFractionHadronicCorrection = val; + if(fApplyFractionHadronicCorrection == kTRUE) + { + SetApplyMIPCorrection(kFALSE); + printf("Enabling Fraction Hadronic Correction \n"); + } + else + { + printf("Disabling Fraction Hadronic Correction \n"); + } +} + +//__________________________________________________________ +void AliJetESDReader::SetFractionHadronicCorrection(Double_t val) +{ + // + // Set value to fFractionHadronicCorrection (default is 0.3) + // apply EMC hadronic correction fApplyFractionHadronicCorrection + // - exclusive with fApplyMIPCorrection + // + + fFractionHadronicCorrection = val; + if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0) + { + SetApplyFractionHadronicCorrection(kTRUE); + printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection); + } + else + { + SetApplyFractionHadronicCorrection(kFALSE); + } +} + + //____________________________________________________________________________ void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) { // Connect the tree @@ -247,7 +311,7 @@ void AliJetESDReader::CreateTasks(TChain* tree) // Create global reader task for analysis fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); // Create a task for to fill the charged particle information - fFillUAFromTracks = new AliJetFillUnitArrayTracks(); + fFillUAFromTracks = new AliJetESDFillUnitArrayTracks(); fFillUAFromTracks->SetReaderHeader(fReaderHeader); fFillUAFromTracks->SetGeom(fGeom); fFillUAFromTracks->SetTPCGrid(fTpcGrid); @@ -260,15 +324,17 @@ void AliJetESDReader::CreateTasks(TChain* tree) fFillUAFromTracks->SetGrid3(fGrid3); fFillUAFromTracks->SetGrid4(fGrid4); } - fFillUAFromTracks->SetHadCorrection(fHCorrection); + fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection); fFillUAFromTracks->SetHadCorrector(fHadCorr); // Create a task for to fill the neutral particle information - fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits(); + fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits(); fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); fFillUAFromEMCalDigits->SetGeom(fGeom); fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); - fFillUAFromEMCalDigits->SetEleCorrection(fECorrection); + fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection); + fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection); + fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection); // Add the task to global task fFillUnitArray->Add(fFillUAFromTracks); fFillUnitArray->Add(fFillUAFromEMCalDigits); @@ -289,7 +355,6 @@ Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray) fProcId = procid; fRefArray = refArray; - vector vtmp(3); // clear momentum array ClearArray(); @@ -446,15 +511,6 @@ void AliJetESDReader::InitUnitArray() fGrid4->SetIndexIJ(); n4 = fGrid4->GetNEntries(); - if(fDebug>1) - { - cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl; - cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl; - cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl; - cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl; - cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl; - } - nGaps = n0+n1+n2+n3+n4; } diff --git a/JETAN/AliJetESDReader.h b/JETAN/AliJetESDReader.h index d42595c8718..e4d6b54f8d4 100755 --- a/JETAN/AliJetESDReader.h +++ b/JETAN/AliJetESDReader.h @@ -14,8 +14,6 @@ // Author : magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- -#include - #include "AliJetReader.h" #include "AliJetUnitArray.h" #include "AliJetGrid.h" @@ -50,16 +48,18 @@ class AliJetESDReader : public AliJetReader void SetTPCGrid(AliJetGrid *grid) {fTpcGrid = grid;} void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;} // Correction of hadronic energy - void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} - void SetHadronCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} - void SetElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;} + void SetHadronCorrector(AliJetHadronCorrection* corr) {fHadCorr = corr;} + void SetApplyElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;} + void SetApplyMIPCorrection(Bool_t val); + void SetApplyFractionHadronicCorrection(Bool_t val); + void SetFractionHadronicCorrection(Double_t val); protected: AliJetDummyGeo *fGeom; //! EMCAL Geometry TChain *fChain; //! chain for reconstructed tracks - TChain *fTree; //! tree for reconstructed tracks + TChain *fTree; //! tree to get the EMCal geometry AliESDEvent *fESD; //! pointer to esd - AliJetHadronCorrectionv1 *fHadCorr; //! Pointer to Hadron Correction Object + AliJetHadronCorrection *fHadCorr; //! Pointer to Hadron Correction Object AliJetGrid *fTpcGrid; //! Pointer to grid object AliJetGrid *fEmcalGrid; //! Pointer to grid object AliJetGrid *fGrid0; // Pointer to grid object @@ -68,9 +68,11 @@ class AliJetESDReader : public AliJetReader AliJetGrid *fGrid3; // Pointer to grid object AliJetGrid *fGrid4; // Pointer to grid object Float_t fPtCut; // Pt cut for tracks to minimise background contribution - Int_t fHCorrection; // Hadron correction flag - Int_t fECorrection; // Electron correction flag + Int_t fApplyElectronCorrection; // Electron correction flag Bool_t fEFlag; // if (fEFlag == kFALSE) => fECorrection automatically setted + Bool_t fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection + Bool_t fApplyFractionHadronicCorrection; // Another type of charged particle energy deposition in EMC + Double_t fFractionHadronicCorrection; // Fraction of momentum of the TPC track to be subtracted from EMC tower Int_t fNumUnits; // Number of units in the unit object array // (same as num towers in EMCAL) Int_t fDebug; //! Debug option diff --git a/JETAN/AliJetFillUnitArray.cxx b/JETAN/AliJetFillUnitArray.cxx new file mode 100644 index 00000000000..e2a91abf9a4 --- /dev/null +++ b/JETAN/AliJetFillUnitArray.cxx @@ -0,0 +1,143 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +//====================================================================== +// ***July 2009 +// Fill Unit Array class +// Base class used by AliJetESDReader to fill a UnitArray from the information extracted +// from the particle tracks or emcal cells +// Author: magali.estienne@subatech.in2p3.fr +//====================================================================== + + +#include "AliJetFillUnitArray.h" + +// --- ROOT system --- +class TSystem; +class TLorentzVector; +class TVector3; +class TGeoManager; +class TProcessID; + +// --- AliRoot header files --- +class AliJetFinder; +class AliJetReader; +class AliJetESDReader; +class AliJetESDReaderHeader; +class AliJetUnitArray; + +ClassImp(AliJetFillUnitArray) + +//_____________________________________________________________________________ +AliJetFillUnitArray::AliJetFillUnitArray() + : TTask("AliJetFillUnitArray","Fill Unit Array with tpc/its and emcal information"), + fNTracks(0), + fNTracksCut(0), + fOpt(0), + fDZ(0), + fDebug(0), + fReaderHeader(0x0), + fMomentumArray(0x0), + fUnitArray(0x0), + fRefArray(0x0), + fProcId(kFALSE), + fTPCGrid(0x0), + fEMCalGrid(0x0), + fGeom(0x0), + fNphi(0), + fNeta(0), + fGrid(0), + fPhi2(0), + fEta2(0), + fIndex(0x0), + fParams(0x0), + fPhiMin(0), + fPhiMax(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinInTPCAcc(0), + fPhiBinInTPCAcc(0), + fEtaBinInEMCalAcc(0), + fPhiBinInEMCalAcc(0), + fNbinPhi(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetFillUnitArray::~AliJetFillUnitArray() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliJetFillUnitArray::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi) +{ + // Get the eta,phi position from the index + + for(Int_t j=0; jAt(i); + phi = fPhi2->At(j); + } + } + + // TPC-EMCAL grid + //------------------------------------- + Int_t ii = 0; + if(i==0) ii = 0; + if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; + if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && iAt(i); + phi = fPhi2->At(j); + } + + if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && + ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) { + if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);} + else eta = fEta2->At(i); + phi = fPhi2->At(j); + } + + if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) { + eta = fEta2->At(i); + phi = fPhi2->At(j); + } + } + } + } +} + +//_____________________________________________________________________________ +Float_t AliJetFillUnitArray::EtaToTheta(Float_t arg) +{ + return 2.*atan(exp(-arg)); +} + + + + + + + diff --git a/JETAN/AliJetFillUnitArray.h b/JETAN/AliJetFillUnitArray.h new file mode 100644 index 00000000000..1c190b57ed3 --- /dev/null +++ b/JETAN/AliJetFillUnitArray.h @@ -0,0 +1,122 @@ +#ifndef ALIJETFILLUNITARRAY_H +#define ALIJETFILLUNITARRAY_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// Base class used to fill Unit Array +// Called by ESD Reader for jet analysis +// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr) +//--------------------------------------------------------------------- + +//#include +#include +#include + +#include +#include +#include + +#include "AliJetReaderHeader.h" +#include "AliJetAODReaderHeader.h" +#include "AliJetDummyGeo.h" +#include "AliJetGrid.h" +#include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "AliJetHadronCorrection.h" + +class Riostream; + +class AliEMCALGeometry; +class AliJetReader; +class AliJetESDReader; + +class AliJetFillUnitArray : public TTask +{ + public: + AliJetFillUnitArray(); + virtual ~AliJetFillUnitArray(); + AliJetFillUnitArray(const AliJetFillUnitArray &det); + AliJetFillUnitArray &operator=(const AliJetFillUnitArray &det); + + // Setter + virtual void SetReaderHeader(AliJetReaderHeader* const readerHeader) {fReaderHeader = readerHeader;} + virtual void SetGeom(AliJetDummyGeo* const geom) {fGeom = geom;} + virtual void SetMomentumArray(TClonesArray* const momentumArray) {fMomentumArray = momentumArray;} + virtual void SetUnitArray(TClonesArray* const unitArray) {fUnitArray = unitArray;} + virtual void SetRefArray(TRefArray* const refArray) {fRefArray = refArray;} + virtual void SetTPCGrid(AliJetGrid* const grid) {fTPCGrid = grid;} + virtual void SetEMCalGrid(AliJetGrid* const grid) {fEMCalGrid = grid;} + virtual void SetProcId(Bool_t id) {fProcId = id;} + virtual void SetGrid0(AliJetGrid */*grid0*/) {;} + virtual void SetGrid1(AliJetGrid */*grid1*/) {;} + virtual void SetGrid2(AliJetGrid */*grid2*/) {;} + virtual void SetGrid3(AliJetGrid */*grid3*/) {;} + virtual void SetGrid4(AliJetGrid */*grid4*/) {;} + virtual void SetHadCorrector(AliJetHadronCorrection* /*corr*/) {;} + virtual void SetApplyMIPCorrection(Bool_t /*val*/) {;} + virtual void SetESD(AliESDEvent */*esd*/) {;} + virtual void SetAOD(AliAODEvent */*aod*/) {;} + virtual void SetApplyElectronCorrection(Int_t /*flag*/) {;} + virtual void SetApplyFractionHadronicCorrection(Bool_t /*val*/) {;} + virtual void SetFractionHadronicCorrection(Double_t /*val*/) {;} + virtual void SetInitMult(Int_t /*mult*/) {;} + virtual void SetInitMultCut(Int_t /*multcut*/) {;} + + // Getter + virtual TClonesArray* GetUnitArray() const {return fUnitArray;} + virtual TRefArray* GetRefArray() const {return fRefArray;} + virtual void GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi); + virtual Int_t GetNeta() const {return fNeta;} + virtual Int_t GetNphi() const {return fNphi;} + virtual Int_t GetMult() const {return 0;} // To be checked + virtual Int_t GetMultCut() const {return 0;} // To be checked + virtual Bool_t GetProcId() const {return fProcId;} + + // Other + virtual void Exec(Option_t* const /*option*/) {;} + virtual Float_t EtaToTheta(Float_t arg); + virtual void InitParameters() {;} + + protected: + Int_t fNTracks; // Number of tracks stored in UnitArray + Int_t fNTracksCut; // Number of tracks stored in UnitArray with a pt cut + Int_t fOpt; // Detector to be used for jet reconstruction + Bool_t fDZ; // Use or not dead zones + Int_t fDebug; // Debug option + + AliJetReaderHeader *fReaderHeader; // ReaderHeader + TClonesArray *fMomentumArray; // MomentumArray + TClonesArray *fUnitArray; // UnitArray + TRefArray *fRefArray; // UnitArray + Bool_t fProcId; // Bool_t for TProcessID synchronization + AliJetGrid *fTPCGrid; // Define filled grid + AliJetGrid *fEMCalGrid; // Define filled grid + AliJetDummyGeo *fGeom; // Define EMCal geometry + + Int_t fNphi; // number of points in the grid: phi + Int_t fNeta; // " eta + Int_t fGrid; // Select the grid acceptance you want to fill + // 0 = TPC acceptance, 1 = TPC-EMCal acceptance + TArrayD *fPhi2; // grid points in phi + TArrayD *fEta2; // grid points in eta + TMatrixD *fIndex; // grid points in (phi,eta) + TMatrixD *fParams; // matrix of parameters in the grid points + Float_t fPhiMin; // EMCal acceptance + Float_t fPhiMax; // EMCal acceptance + Float_t fEtaMin; // EMCal acceptance + Float_t fEtaMax; // EMCal acceptance + Int_t fEtaBinInTPCAcc; // Number of bins in Eta in TPC acceptance + Int_t fPhiBinInTPCAcc; // Number of bins in phi in TPC acceptance + Int_t fEtaBinInEMCalAcc;// Number of bins in Eta in EMCal acceptance + Int_t fPhiBinInEMCalAcc;// Number of bins in phi in EMCal acceptance + Int_t fNbinPhi; // Number of phi bins + + private: + + + ClassDef(AliJetFillUnitArray,1) // Fill Unit Array with tpc and/or emcal information +}; + +#endif diff --git a/JETAN/AliJetFinder.cxx b/JETAN/AliJetFinder.cxx index 25239894afb..bbdaa671b30 100644 --- a/JETAN/AliJetFinder.cxx +++ b/JETAN/AliJetFinder.cxx @@ -25,19 +25,18 @@ #include #include -#include -#include #include "AliJetFinder.h" #include "AliJet.h" #include "AliAODJet.h" -#include "AliJetReader.h" -#include "AliJetReaderHeader.h" #include "AliJetControlPlots.h" #include "AliLeading.h" #include "AliAODEvent.h" #include "AliJetUnitArray.h" +class TProcessID; +class TClonesArray; + ClassImp(AliJetFinder) AliJetFinder::AliJetFinder(): @@ -252,7 +251,6 @@ Bool_t AliJetFinder::ProcessEvent2() FindJets(); Int_t nEntRef = ref->GetEntries(); - vector vtmp; for(Int_t i=0; iAt(i))->SetUnitEnergy(0.); ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller); ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller); - ((AliJetUnitArray*)ref->At(i))->SetUnitPxPyPz(kTRUE,vtmp); ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad); + ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad); ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc); ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet); + ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef(); // Reset process ID AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i); diff --git a/JETAN/AliJetFinder.h b/JETAN/AliJetFinder.h index a782d7c8e3b..3d6a100be0d 100755 --- a/JETAN/AliJetFinder.h +++ b/JETAN/AliJetFinder.h @@ -12,14 +12,20 @@ // magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- +//#include +//#include + #include #include "AliAODJet.h" #include "AliJetReader.h" +#include "AliJetHeader.h" +#include "AliJet.h" +#include "AliJetReaderHeader.h" class TFile; class TChain; -class AliJet; -class AliJetHeader; +//class AliJet; +//class AliJetHeader; class AliJetControlPlots; class AliLeading; class AliAODJet; @@ -53,7 +59,7 @@ class AliJetFinder : public TObject virtual void Reset() {fNAODjets = 0;} virtual void FindJets() {} virtual void FindJetsC(){} - virtual void WriteJHeaderToFile() { } + virtual void WriteJHeaderToFile() {} // some methods to allow steering from the outside virtual Bool_t ProcessEvent(); virtual Bool_t ProcessEvent2(); @@ -63,7 +69,7 @@ class AliJetFinder : public TObject virtual void ConnectAODNonStd(AliAODEvent* aod,const char* bname); virtual TTree *MakeTreeJ(char* name); virtual void WriteHeaders(); - virtual void WriteJetsToFile() {;} + virtual void WriteJetsToFile() const {;} virtual void TestJet() {;} protected: diff --git a/JETAN/AliJetReader.cxx b/JETAN/AliJetReader.cxx index 2364b971627..88a3a4c50e0 100755 --- a/JETAN/AliJetReader.cxx +++ b/JETAN/AliJetReader.cxx @@ -29,8 +29,8 @@ #include "AliJetReaderHeader.h" #include "AliESDEvent.h" #include "AliHeader.h" -#include "AliJetFillUnitArrayTracks.h" -#include "AliJetFillUnitArrayEMCalDigits.h" +#include "AliJetESDFillUnitArrayTracks.h" +#include "AliJetESDFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" #include "AliJetHadronCorrectionv1.h" @@ -51,10 +51,9 @@ AliJetReader::AliJetReader(): fSignalFlag(0), fCutFlag(0), fUnitArray(new TClonesArray("AliJetUnitArray",60000)), - fUnitArrayNoCuts(new TClonesArray("AliJetUnitArray",60000)), fArrayInitialised(0), - fFillUAFromTracks(new AliJetFillUnitArrayTracks()), - fFillUAFromEMCalDigits(new AliJetFillUnitArrayEMCalDigits()), + fFillUAFromTracks(new AliJetESDFillUnitArrayTracks()), + fFillUAFromEMCalDigits(new AliJetESDFillUnitArrayEMCalDigits()), fNumCandidate(0), fNumCandidateCut(0), fHadronCorrector(0), @@ -82,11 +81,6 @@ AliJetReader::~AliJetReader() delete fUnitArray; } - if (fUnitArrayNoCuts) { - fUnitArrayNoCuts->Delete(); - delete fUnitArrayNoCuts; - } - if (fFillUnitArray) { delete fFillUnitArray; } diff --git a/JETAN/AliJetReader.h b/JETAN/AliJetReader.h index 457f3886a29..72aa1d41e69 100755 --- a/JETAN/AliJetReader.h +++ b/JETAN/AliJetReader.h @@ -22,10 +22,9 @@ class AliJetReaderHeader; class AliESDEvent; class AliHeader; class AliJetUnitArray; -class AliJetHadronCorrectionv1; +class AliJetHadronCorrection; class AliJet; -class AliJetFillUnitArrayTracks; -class AliJetFillUnitArrayEMCalDigits; +class AliJetFillUnitArray; class AliJetReader : public TObject { @@ -37,8 +36,6 @@ class AliJetReader : public TObject virtual TClonesArray* GetMomentumArray() const {return fMomentumArray;} virtual TRefArray* GetReferences() const {return 0;} virtual TClonesArray *GetUnitArray() const {return fUnitArray;} - virtual TClonesArray *GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;} - virtual AliJetReaderHeader* GetReaderHeader() const {return fReaderHeader;} virtual AliHeader *GetAliHeader() const {return fAliHeader;} virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];} @@ -50,17 +47,16 @@ class AliJetReader : public TObject // Setters virtual Bool_t FillMomentumArray() {return kTRUE;} virtual Bool_t ReadEventLoader(Int_t) {return kTRUE;} - virtual void FillUnitArrayFromTPCTracks(Int_t) {} // temporarily not used - virtual void FillUnitArrayFromEMCALHits() {} // temporarily not used - virtual void FillUnitArrayFromEMCALDigits(Int_t) {} // temporarily not used virtual void InitUnitArray() {} virtual void InitParameters() {} virtual void CreateTasks(TChain* /*tree*/) {} virtual Bool_t ExecTasks(Bool_t /*procid*/, TRefArray* /*refArray*/) {return kFALSE;} // Correction of hadronic energy - virtual void SetHadronCorrector(AliJetHadronCorrectionv1*) {;} - virtual void SetHadronCorrection(Int_t) {;} - virtual void SetElectronCorrection(Int_t) {;} + virtual void SetHadronCorrector(AliJetHadronCorrection*) {;} + virtual void SetApplyMIPCorrection(Bool_t /*val*/){;} + virtual void SetApplyFractionHadronicCorrection(Bool_t /*val*/){;} + virtual void SetFractionHadronicCorrection(Double_t /*val*/){;} + virtual void SetApplyElectronCorrection(Int_t /*flag*/) {;} virtual void SetReaderHeader(AliJetReaderHeader* header) {fReaderHeader = header;} virtual void SetESD(AliESDEvent* esd) { fESD = esd;} @@ -89,14 +85,13 @@ class AliJetReader : public TObject // from the underlying event TArrayI fCutFlag; // to flag if a particle passed the pt cut or not TClonesArray *fUnitArray; // array of digit position and energy - TClonesArray *fUnitArrayNoCuts; //|| array of digit position and energy Bool_t fArrayInitialised; // To check that array of units is initialised - AliJetFillUnitArrayTracks *fFillUAFromTracks; // For charged particle task - AliJetFillUnitArrayEMCalDigits *fFillUAFromEMCalDigits; // For neutral particle task + AliJetFillUnitArray *fFillUAFromTracks; // For charged particle task + AliJetFillUnitArray *fFillUAFromEMCalDigits; // For neutral particle task Int_t fNumCandidate; // Number of entries different from zero in unitarray Int_t fNumCandidateCut; // Number of entries different from zero in unitarray // which pass pt cut - AliJetHadronCorrectionv1 *fHadronCorrector; //! Pointer to hadronic correction + AliJetHadronCorrection *fHadronCorrector; //! Pointer to hadronic correction Int_t fHCorrection; // Hadron correction flag Int_t fECorrection; // Electron correction flag Bool_t fEFlag; // Electron correction flag diff --git a/JETAN/AliJetUnitArray.cxx b/JETAN/AliJetUnitArray.cxx index 5b55049f268..28ffcd65eee 100644 --- a/JETAN/AliJetUnitArray.cxx +++ b/JETAN/AliJetUnitArray.cxx @@ -22,16 +22,14 @@ //*-- Author: Sarah Blyth (LBL/UCT) // -- // Revised Version for JETAN -// -- Magali Estienne (IReS) - -//#include - -#include -#include -#include +// -- Magali Estienne (magali.estienne@subatech.in2p3.fr) #include "AliJetUnitArray.h" +class TVector3; +class TLorentzVector; +class TClonesArray; + ClassImp(AliJetUnitArray) AliJetUnitArray::AliJetUnitArray(): @@ -53,10 +51,11 @@ AliJetUnitArray::AliJetUnitArray(): fUnitPy(0.), fUnitPz(0.), fUnitMass(0.), - fV(0), fVc(0), fVn(0), - fVet(0) + fUnitTrackRef(new TRefArray), + fUnitCellRef(new TRefArray), + fUnitClusterRef(new TRefArray) { // Default constructor } @@ -74,16 +73,17 @@ AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t fUnitFlag(inout), fUnitCutFlag(cut), fUnitCutFlag2(cut2), - fUnitSignalFlag(signal), + fUnitSignalFlag(kBad), fUnitDetectorFlag(det), fUnitPx(0.), fUnitPy(0.), fUnitPz(0.), fUnitMass(mass), - fV(0), fVc(0), fVn(0), - fVet(0) + fUnitTrackRef(new TRefArray), + fUnitCellRef(new TRefArray), + fUnitClusterRef(new TRefArray) { //abs ID (in a eta,phi grid, track ID in ESD, eta, phi, energy, px, py, pz, Deta, Dphi, detector flag, in/out jet, mass @@ -109,18 +109,94 @@ AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t fUnitPy(py), fUnitPz(pz), fUnitMass(mass), - fV(0), fVc(0), fVn(0), - fVet(0) + fUnitTrackRef(new TRefArray), + fUnitCellRef(new TRefArray), + fUnitClusterRef(new TRefArray) { // Constructor 2 } +AliJetUnitArray::AliJetUnitArray(const AliJetUnitArray& rUnit): + TObject(rUnit), + fUnitEnergy(rUnit.fUnitEnergy), + fUnitEta(rUnit.fUnitEta), + fUnitPhi(rUnit.fUnitPhi), + fUnitDeta(rUnit.fUnitDeta), + fUnitDphi(rUnit.fUnitDphi), + fUnitID(rUnit.fUnitID), + fUnitTrackID(rUnit.fUnitTrackID), + fUnitNum(rUnit.fUnitNum), + fUnitClusterID(rUnit.fUnitClusterID), + fUnitFlag(rUnit.fUnitFlag), + fUnitCutFlag(rUnit.fUnitCutFlag), + fUnitCutFlag2(rUnit.fUnitCutFlag2), + fUnitSignalFlag(rUnit.fUnitSignalFlag), + fUnitDetectorFlag(rUnit.fUnitDetectorFlag), + fUnitPx(rUnit.fUnitPx), + fUnitPy(rUnit.fUnitPy), + fUnitPz(rUnit.fUnitPz), + fUnitMass(rUnit.fUnitMass), + fVc(rUnit.fVc), + fVn(rUnit.fVn), + fUnitTrackRef(rUnit.fUnitTrackRef), + fUnitCellRef(rUnit.fUnitCellRef), + fUnitClusterRef(rUnit.fUnitClusterRef) +{ + // Copy constructor +} + +AliJetUnitArray& AliJetUnitArray::operator=(const AliJetUnitArray& rhs) +{ + // Assignment + fUnitEnergy = rhs.fUnitEnergy; + fUnitEta = rhs.fUnitEta; + fUnitPhi = rhs.fUnitPhi; + fUnitDeta = rhs.fUnitDeta; + fUnitDphi = rhs.fUnitDphi; + fUnitID = rhs.fUnitID; + fUnitTrackID = rhs.fUnitTrackID; + fUnitNum = rhs.fUnitNum; + fUnitClusterID = rhs.fUnitClusterID; + fUnitFlag = rhs.fUnitFlag; + fUnitCutFlag = rhs.fUnitCutFlag; + fUnitCutFlag2 = rhs.fUnitCutFlag2; + fUnitSignalFlag = rhs.fUnitSignalFlag; + fUnitDetectorFlag = rhs.fUnitDetectorFlag; + fUnitPx = rhs.fUnitPx; + fUnitPy = rhs.fUnitPy; + fUnitPz = rhs.fUnitPz; + fUnitMass = rhs.fUnitMass; + fVc = rhs.fVc; + fVn = rhs.fVn; + fUnitTrackRef = rhs.fUnitTrackRef; + fUnitCellRef = rhs.fUnitCellRef; + fUnitClusterRef = rhs.fUnitClusterRef; + + return *this; + +} + + //------------------------------------------------------------------------ AliJetUnitArray::~AliJetUnitArray() { // Destructor + delete fUnitTrackRef; + delete fUnitCellRef; + delete fUnitClusterRef; + +} + +void AliJetUnitArray::ClearUnitTrackRef() +{ + fUnitTrackRef->Clear(); +} + +void AliJetUnitArray::ClearUnitCellRef() +{ + fUnitCellRef->Clear(); } //------------------------------------------------------------------------ @@ -145,32 +221,6 @@ void AliJetUnitArray::SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlag else fVn.push_back(flag); } -//------------------------------------------------------------------------ -void AliJetUnitArray::SetUnitEtN(Bool_t init, Float_t et) -{ - // Set transverse energy of the neutral cell - if(init){ - if(!fVet.empty()) - fVet.clear(); - } - else fVet.push_back(et); -} - - -//------------------------------------------------------------------------ -void AliJetUnitArray::SetUnitPxPyPz(Bool_t init, vector v3) -{ - // Set momentum components of the charged particle - if(init) - { - if(!fV.empty()){ - fV.clear(); - } - } - else{ - fV.push_back(v3); - } -} //------------------------------------------------------------------------ Bool_t AliJetUnitArray::GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc) @@ -196,47 +246,6 @@ Bool_t AliJetUnitArray::GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlag else return kFALSE; } -//------------------------------------------------------------------------ -Bool_t AliJetUnitArray::GetUnitEtN(Int_t ind, Float_t &et) -{ - // Get transverse energy of the neutral cell - if(ind <= (Int_t)fVet.size()) - { - et = (Float_t)fVet[ind]; - return kTRUE; - } - else return kFALSE; -} - -//------------------------------------------------------------------------ -Bool_t AliJetUnitArray::GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz) -{ - // Get momentum components of the charged particle - if(ind <= (Int_t)fV.size()) - { - px = (Float_t)fV[ind][0]; - py = (Float_t)fV[ind][1]; - pz = (Float_t)fV[ind][2]; - return kTRUE; - } - else return kFALSE; -} - -//------------------------------------------------------------------------ -Bool_t AliJetUnitArray::GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en) -{ -// Get 4-momentum components of the charged particle - if(ind <= (Int_t)fV.size()) - { - px = (Float_t)fV[ind][0]; - py = (Float_t)fV[ind][1]; - pz = (Float_t)fV[ind][2]; - en = TMath::Sqrt(px*px+py*py+pz*pz); - return kTRUE; - } - else return kFALSE; -} - //------------------------------------------------------------------------ Float_t AliJetUnitArray::EtaToTheta(Float_t arg) const { diff --git a/JETAN/AliJetUnitArray.h b/JETAN/AliJetUnitArray.h index 1673cda21be..5300bc2546b 100644 --- a/JETAN/AliJetUnitArray.h +++ b/JETAN/AliJetUnitArray.h @@ -15,14 +15,22 @@ // // -#include -#include - +#include #include -#include -#include +#include +#include +#include "AliESDtrack.h" +#include "AliAODTrack.h" +#include "AliVParticle.h" +#include "AliESDCaloCells.h" +#include "AliAODCaloCells.h" +#include "AliESDCaloCluster.h" +#include "AliAODCaloCluster.h" #include "AliJetFinderTypes.h" +class TList; +class TVector3; + class AliJetUnitArray : public TObject { public: @@ -68,36 +76,36 @@ class AliJetUnitArray : public TObject { fUnitDetectorFlag = detectorflag; } - void SetUnitPxPyPz(Bool_t init, vector v3); - void SetUnitEtN(Bool_t init, Float_t et); // Added for background studies void SetUnitSignalFlagC(Bool_t init, AliJetFinderUnitSignalFlagType_t flag); void SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlagType_t flag); void SetUnitMass(Float_t mass) {fUnitMass = mass;} + void SetUnitTrackRef(TRefArray* trackref) {fUnitTrackRef = trackref;} + void SetUnitCellRef(TRefArray* cellref) {fUnitCellRef = cellref;} + void SetUnitClusterRef(TRefArray* clusterref) {fUnitClusterRef = clusterref;} // Getter - Float_t GetUnitEnergy() const {return fUnitEnergy;} - Float_t GetUnitEta() const {return fUnitEta;} - Float_t GetUnitPhi() const {return fUnitPhi;} - Float_t GetUnitPx() const {return fUnitPx;} - Float_t GetUnitPy() const {return fUnitPy;} - Float_t GetUnitPz() const {return fUnitPz;} - Float_t GetUnitDeta() const {return fUnitDeta;} - Float_t GetUnitDphi() const {return fUnitDphi;} - Int_t GetUnitID() const {return fUnitID;} - Int_t GetUnitTrackID() const {return fUnitTrackID;} - Int_t GetUnitEntries() const {return fUnitNum;} - Int_t GetUnitClusterID() const {return fUnitClusterID;} - Float_t GetUnitMass() const {return fUnitMass;} - Bool_t GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz); - Bool_t GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en); - Bool_t GetUnitEtN(Int_t ind, Float_t &et); // Added for background studies - Bool_t GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc); - Bool_t GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn); - Int_t GetUnitVectorSize() {return fV.size();} - Int_t GetUnitVectorSizeC() {return fVc.size();} - Int_t GetUnitVectorSizeN() {return fVn.size();} - - Float_t EtaToTheta(Float_t arg) const; + Float_t GetUnitEnergy() const {return fUnitEnergy;} + Float_t GetUnitEta() const {return fUnitEta;} + Float_t GetUnitPhi() const {return fUnitPhi;} + Float_t GetUnitPx() const {return fUnitPx;} + Float_t GetUnitPy() const {return fUnitPy;} + Float_t GetUnitPz() const {return fUnitPz;} + Float_t GetUnitDeta() const {return fUnitDeta;} + Float_t GetUnitDphi() const {return fUnitDphi;} + Int_t GetUnitID() const {return fUnitID;} + Int_t GetUnitTrackID() const {return fUnitTrackID;} + Int_t GetUnitEntries() const {return fUnitNum;} + Int_t GetUnitClusterID() const {return fUnitClusterID;} + Float_t GetUnitMass() const {return fUnitMass;} + Bool_t GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc); + Bool_t GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn); + TRefArray* GetUnitTrackRef() const {return fUnitTrackRef;} + TRefArray* GetUnitCellRef() const {return fUnitCellRef;} + TRefArray* GetUnitClusterRef() const {return fUnitClusterRef;} + + Float_t EtaToTheta(Float_t arg) const; + void ClearUnitTrackRef(); + void ClearUnitCellRef(); AliJetFinderUnitFlagType_t GetUnitFlag() const { @@ -121,6 +129,8 @@ class AliJetUnitArray : public TObject } protected: + AliJetUnitArray(const AliJetUnitArray& rUnit); + AliJetUnitArray& operator = (const AliJetUnitArray& rhs); Bool_t operator> ( AliJetUnitArray unit1) const; Bool_t operator< ( AliJetUnitArray unit1) const; Bool_t operator== ( AliJetUnitArray unit1) const; @@ -143,10 +153,11 @@ class AliJetUnitArray : public TObject Float_t fUnitPy; // Py of charged track Float_t fUnitPz; // Pz of charged track Float_t fUnitMass; // Mass of particle - vector< vector< Float_t > > fV; //|| vector to store part information in each cell vector< AliJetFinderUnitSignalFlagType_t > fVc; //|| added for background studies vector< AliJetFinderUnitSignalFlagType_t > fVn; //|| added for background studies - vector< Float_t > fVet; //|| added for background studies + TRefArray* fUnitTrackRef; //! pointer to array of references to esd tracks + TRefArray* fUnitCellRef; //! pointer to array of references to esd cells + TRefArray* fUnitClusterRef; //! pointer to array of references to esd clusters ClassDef(AliJetUnitArray,1) diff --git a/JETAN/AliUA1JetFinderV2.cxx b/JETAN/AliUA1JetFinderV2.cxx index 9e11e175c5c..04decb9fc9b 100644 --- a/JETAN/AliUA1JetFinderV2.cxx +++ b/JETAN/AliUA1JetFinderV2.cxx @@ -23,26 +23,22 @@ // Author: magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- -#include -#include - -#include #include -#include #include #include #include #include #include +#include "TFile.h" #include "AliUA1JetFinderV2.h" #include "AliUA1JetHeaderV1.h" #include "AliJetUnitArray.h" -#include "AliJetReaderHeader.h" -#include "AliJetReader.h" -#include "AliJet.h" -#include "AliAODJet.h" +class TArrayF; +class TFile; +class AliJetReader; +class AliAODJet; ClassImp(AliUA1JetFinderV2) @@ -109,8 +105,8 @@ void AliUA1JetFinderV2::FindJetsC() ptT[i] = lv->Pt(); etaT[i] = lv->Eta(); phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi()); - cFlagT[i] = fReader->GetCutFlag(i); // Temporarily added - sFlagT[i] = fReader->GetSignalFlag(i); // Temporarily added peut-etre a mettre apres cut en pt !!! + cFlagT[i] = fReader->GetCutFlag(i); + sFlagT[i] = fReader->GetSignalFlag(i); if (fReader->GetCutFlag(i) != 1) continue; fLego->Fill(etaT[i], phiT[i], ptT[i]); @@ -127,10 +123,10 @@ void AliUA1JetFinderV2::FindJetsC() Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS); // arrays to hold jets - Float_t* etaJet = new Float_t[30]; - Float_t* phiJet = new Float_t[30]; - Float_t* etJet = new Float_t[30]; - Float_t* etsigJet = new Float_t[30]; //signal et in jet + Float_t* etaJet = new Float_t[30]; // eta jet + Float_t* phiJet = new Float_t[30]; // phi jet + Float_t* etJet = new Float_t[30]; // et jet + Float_t* etsigJet = new Float_t[30]; // signal et in jet Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable) Int_t* ncellsJet = new Int_t[30]; Int_t* multJet = new Int_t[30]; @@ -138,7 +134,7 @@ void AliUA1JetFinderV2::FindJetsC() Float_t* etaJetOk = new Float_t[30]; Float_t* phiJetOk = new Float_t[30]; Float_t* etJetOk = new Float_t[30]; - Float_t* etsigJetOk = new Float_t[30]; //signal et in jet + Float_t* etsigJetOk = new Float_t[30]; // signal et in jet Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable) Int_t* ncellsJetOk = new Int_t[30]; Int_t* multJetOk = new Int_t[30]; @@ -181,7 +177,6 @@ void AliUA1JetFinderV2::FindJetsC() //subtract background Float_t etbgTotalN = 0.0; //new background if(header->GetBackgMode() == 1) // standard - // SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); if(header->GetBackgMode() == 2) //cone SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); @@ -210,6 +205,7 @@ void AliUA1JetFinderV2::FindJetsC() phiJetOk[p] = phiJet[idx[p]]; etJetOk[p] = etJet[idx[p]]; etallJetOk[p] = etJet[idx[p]]; + etsigJetOk[p] = etsigJet[idx[p]]; ncellsJetOk[p] = ncellsJet[idx[p]]; multJetOk[p] = multJet[idx[p]]; } @@ -325,8 +321,9 @@ void AliUA1JetFinderV2::FindJets() Int_t nCand = fReader->GetNumCandidate(); Int_t nCandCut = fReader->GetNumCandidateCut(); Int_t nIn = fUnit->GetEntries(); - Float_t fPtMin = fReader->GetReaderHeader()->GetPtCut(); + Float_t ptMin = fReader->GetReaderHeader()->GetPtCut(); + fDebug = fReader->GetReaderHeader()->GetDebug(); if (nIn == 0) return; Int_t nCandidateCut = 0; @@ -337,36 +334,34 @@ void AliUA1JetFinderV2::FindJets() // local arrays for input No Cuts // Both pt < ptMin and pt > ptMin - Float_t* ptT = new Float_t[nCandidate]; - Float_t* en2T = new Float_t[nCandidate]; - Float_t* pt2T = new Float_t[nCandidate]; - Int_t* detT = new Int_t[nCandidate]; - Float_t* etaT = new Float_t[nCandidate]; - Float_t* phiT = new Float_t[nCandidate]; - Float_t* cFlagT = new Float_t[nCandidate]; - Float_t* cFlag2T = new Float_t[nCandidate]; - Float_t* sFlagT = new Float_t[nCandidate]; - Float_t* cClusterT = new Float_t[nCandidate]; - Int_t* vectT = new Int_t[nCandidate]; - Int_t loop1 = 0; - Int_t* injet = new Int_t[nCandidate]; - Int_t* sflag = new Int_t[nCandidate]; - vector< vector > pxT; - vector< vector > pyT; - vector< vector > pzT; + Float_t* ptT = new Float_t[nCandidate]; + Float_t* en2T = new Float_t[nCandidate]; + Float_t* pt2T = new Float_t[nCandidate]; + Int_t* detT = new Int_t[nCandidate]; + Float_t* etaT = new Float_t[nCandidate]; + Float_t* phiT = new Float_t[nCandidate]; + Float_t* cFlagT = new Float_t[nCandidate]; + Float_t* cFlag2T = new Float_t[nCandidate]; + Float_t* sFlagT = new Float_t[nCandidate]; + Float_t* cClusterT = new Float_t[nCandidate]; + Int_t* vectT = new Int_t[nCandidate]; + Int_t loop1 = 0; + Int_t* injet = new Int_t[nCandidate]; + Int_t* sflag = new Int_t[nCandidate]; + TRefArray* trackRef = new TRefArray(); //total energy in array Float_t etbgTotal = 0.0; - TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0); + TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0); // Input cell info - Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray - Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray - Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray - Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray - Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray - Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray - Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray + Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray + Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray + Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray + Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray + Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray + Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray + Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray Int_t *flagCell2 = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray // Information extracted from fUnitArray @@ -376,25 +371,13 @@ void AliUA1JetFinderV2::FindJets() // Recover particle information from UnitArray AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i); - + TRefArray* ref = uArray->GetUnitTrackRef(); + Int_t nRef = ref->GetEntries(); + if(uArray->GetUnitEnergy()>0.){ - vector vtmpx; - vector vtmpy; - vector vtmpz; - for(Int_t j=0; jGetUnitVectorSize();j++) - { - Float_t x=0.; Float_t y=0.; Float_t z=0.; - uArray->GetUnitPxPyPz(j,x,y,z); - vtmpx.push_back(x); - vtmpy.push_back(y); - vtmpz.push_back(z); - } - pxT.push_back(vtmpx); - pyT.push_back(vtmpy); - pzT.push_back(vtmpz); - vtmpx.clear(); - vtmpy.clear(); - vtmpz.clear(); + + for(Int_t jpart=0; jpartAdd((AliVTrack*)ref->At(jpart)); ptT[loop1] = uArray->GetUnitEnergy(); detT[loop1] = uArray->GetUnitDetectorFlag(); etaT[loop1] = uArray->GetUnitEta(); @@ -402,7 +385,7 @@ void AliUA1JetFinderV2::FindJets() cFlagT[loop1]= uArray->GetUnitCutFlag(); // pt cut tpc cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal sFlagT[loop1]= uArray->GetUnitSignalFlag(); - vectT[loop1] = uArray->GetUnitVectorSize(); + vectT[loop1] = nRef; if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) { pt2T[loop1] = 0.; en2T[loop1] = 0.; @@ -414,11 +397,13 @@ void AliUA1JetFinderV2::FindJets() } if(detT[loop1]==0){ // TPC+ITS Float_t pt = 0.; - for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + x = ((AliVTrack*)ref->At(j))->Px(); + y = ((AliVTrack*)ref->At(j))->Py(); + z = ((AliVTrack*)ref->At(j))->Pz(); pt = TMath::Sqrt(x*x+y*y); - if(pt>fPtMin) { + if(pt>ptMin) { pt2T[loop1] += pt; en2T[loop1] += pt; hPtTotal->Fill(pt); @@ -430,11 +415,13 @@ void AliUA1JetFinderV2::FindJets() Float_t ptCTot = 0.; Float_t pt = 0.; Float_t enC = 0.; - for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + x = ((AliVTrack*)ref->At(j))->Px(); + y = ((AliVTrack*)ref->At(j))->Py(); + z = ((AliVTrack*)ref->At(j))->Pz(); pt = TMath::Sqrt(x*x+y*y); - if(pt>fPtMin) { + if(pt>ptMin) { pt2T[loop1]+=pt; en2T[loop1]+=pt; hPtTotal->Fill(pt); @@ -466,12 +453,14 @@ void AliUA1JetFinderV2::FindJets() } if(uArray->GetUnitDetectorFlag()==0){ // TPC case Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.; - for(Int_t j=0; jGetUnitVectorSize();j++) + for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + x = ((AliVTrack*)ref->At(j))->Px(); + y = ((AliVTrack*)ref->At(j))->Py(); + z = ((AliVTrack*)ref->At(j))->Pz(); pt = TMath::Sqrt(x*x+y*y); - if(pt>fPtMin) { + if(pt>ptMin) { et1 += pt; et2 += pt; } @@ -490,12 +479,14 @@ void AliUA1JetFinderV2::FindJets() Float_t ptCTot = 0.; Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.; Float_t enC = 0.; - for(Int_t j=0; jGetUnitVectorSize();j++) + for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + x = ((AliVTrack*)ref->At(j))->Px(); + y = ((AliVTrack*)ref->At(j))->Py(); + z = ((AliVTrack*)ref->At(j))->Pz(); pt = TMath::Sqrt(x*x+y*y); - if(pt>fPtMin) { + if(pt>ptMin) { et1 += pt; et2 += pt; } @@ -625,6 +616,7 @@ void AliUA1JetFinderV2::FindJets() phiJetOk[p] = phiJet[idx[p]]; etJetOk[p] = etJet[idx[p]]; etallJetOk[p] = etJet[idx[p]]; + etsigJetOk[p] = etsigJet[idx[p]]; ncellsJetOk[p] = ncellsJet[idx[p]]; multJetOk[p] = multJet[idx[p]]; } @@ -682,10 +674,7 @@ void AliUA1JetFinderV2::FindJets() fJets->SetEtaIn(etaT); fJets->SetPhiIn(phiT); fJets->SetPtIn(ptT); - fJets->SetVectorSizeIn(vectT); - fJets->SetVectorPxIn(pxT); - fJets->SetVectorPyIn(pyT); - fJets->SetVectorPzIn(pzT); + fJets->SetTrackReferences(trackRef); fJets->SetDetectorFlagIn(detT); fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()))); @@ -695,9 +684,8 @@ void AliUA1JetFinderV2::FindJets() delete pt2T; delete etaT; delete phiT; - pxT.clear(); - pyT.clear(); - pzT.clear(); + trackRef->Delete(); + delete trackRef; delete detT; delete cFlagT; delete cFlag2T; @@ -740,13 +728,18 @@ void AliUA1JetFinderV2::FindJets() //////////////////////////////////////////////////////////////////////// void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, - Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, - Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, - Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, - Float_t* etallJet, Int_t* ncellsJet) + Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2, + const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal, + Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet, + Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet) { - + // + // Main method for jet finding + // UA1 base cone finder + // + Int_t nCell = nIn; + fDebug = fReader->GetReaderHeader()->GetDebug(); // Dump lego // Check enough space! *to be done* @@ -1137,10 +1130,10 @@ void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, Float_t* ptT, - Int_t*vectT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, - Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, - Float_t* etsigJet, Int_t* multJet, Int_t* injet) +void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT, + const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, + const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, + const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet) { // // Background subtraction using cone method but without correction in dE/deta distribution @@ -1149,7 +1142,8 @@ void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, //calculate energy inside and outside cones AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - fOpt = fReader->GetReaderHeader()->GetDetector(); + Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); + fDebug = fReader->GetReaderHeader()->GetDebug(); Float_t rc= header->GetRadius(); Float_t etIn[30]; Float_t etOut = 0; @@ -1342,9 +1336,9 @@ void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, +void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet) { //background subtraction using cone method but without correction in dE/deta distribution @@ -1408,10 +1402,10 @@ void AliUA1JetFinderV2::SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, - Int_t* multJet, Int_t* injet) +void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, + const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, + Float_t* etsigJet, Int_t* multJet, Int_t* injet) { //background subtraction using statistical method @@ -1473,9 +1467,9 @@ void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotal } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT, - Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, +void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT, + Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet) { // Cone background subtraction method taking into acount dEt/deta distribution @@ -1587,9 +1581,9 @@ void AliUA1JetFinderV2::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTota } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, +void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, + Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet) { // Ratio background subtraction method taking into acount dEt/deta distribution diff --git a/JETAN/AliUA1JetFinderV2.h b/JETAN/AliUA1JetFinderV2.h index eef1429af29..5c64d88fce9 100644 --- a/JETAN/AliUA1JetFinderV2.h +++ b/JETAN/AliUA1JetFinderV2.h @@ -12,9 +12,8 @@ // (version in c++) //--------------------------------------------------------------------- -#include - #include "AliJetFinder.h" + class AliUA1JetHeaderV1; class TH2F; class TChain; @@ -34,45 +33,44 @@ class AliUA1JetFinderV2 : public AliJetFinder Float_t* etallJet, Int_t* ncellsJet); void RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, - Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, - Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, + Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2, const Float_t* phiCell2, + const Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet); - void SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + void SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&EtbgTotalN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,Int_t* multJet, Int_t* injet); - void SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, Float_t* ptT, Int_t* vectT, - Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, - Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + void SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&EtbgTotalN, const Float_t* ptT, const Int_t* vectT, + const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* cFlag2T, + const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); - void SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + void SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& EtbgTotalN, + Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); - void SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + void SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& EtbgTotalN, + Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); - void SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + void SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&EtbgTotalN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, + const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); void Reset(); void InitTask(TChain* tree); void WriteJHeaderToFile(); - protected: AliUA1JetFinderV2(const AliUA1JetFinderV2& rJetF1); AliUA1JetFinderV2& operator = (const AliUA1JetFinderV2& rhsf); - TH2F * fLego; //Lego Histo - Int_t fDebug; - Int_t fOpt; + TH2F * fLego; // Lego Histo + Int_t fDebug; // Debug option + Int_t fOpt; // Detector option (charged only or charged+neutral) ClassDef(AliUA1JetFinderV2,1) }; diff --git a/JETAN/JETANLinkDef.h b/JETAN/JETANLinkDef.h index d870cca2556..5d2a01f29fc 100644 --- a/JETAN/JETANLinkDef.h +++ b/JETAN/JETANLinkDef.h @@ -28,8 +28,11 @@ #pragma link C++ class AliJetHadronCorrection+; #pragma link C++ class AliJetHadronCorrectionv0+; #pragma link C++ class AliJetHadronCorrectionv1+; -#pragma link C++ class AliJetFillUnitArrayTracks+; -#pragma link C++ class AliJetFillUnitArrayEMCalDigits+; +#pragma link C++ class AliJetFillUnitArray+; +#pragma link C++ class AliJetESDFillUnitArrayTracks+; +#pragma link C++ class AliJetESDFillUnitArrayEMCalDigits+; +#pragma link C++ class AliJetAODFillUnitArrayTracks+; +#pragma link C++ class AliJetAODFillUnitArrayEMCalDigits+; #pragma link C++ class AliJetDummyGeo+; #pragma link C++ class AliJetDummyShishKebabTrd1Module+; #pragma link C++ class AliAnalysisTaskJets+; diff --git a/JETAN/libJETAN.pkg b/JETAN/libJETAN.pkg index 8ce13b1f7b6..31220eee46c 100644 --- a/JETAN/libJETAN.pkg +++ b/JETAN/libJETAN.pkg @@ -11,7 +11,9 @@ SRCS = AliJet.cxx AliJetHeader.cxx \ AliCdfJetFinder.cxx AliCdfJetHeader.cxx \ AliJetGrid.cxx AliJetUnitArray.cxx AliJetHadronCorrection.cxx \ AliJetHadronCorrectionv0.cxx AliJetHadronCorrectionv1.cxx \ - AliJetFillUnitArrayTracks.cxx AliJetFillUnitArrayEMCalDigits.cxx\ + AliJetFillUnitArray.cxx AliJetESDFillUnitArrayTracks.cxx \ + AliJetESDFillUnitArrayEMCalDigits.cxx AliJetAODFillUnitArrayTracks.cxx \ + AliJetAODFillUnitArrayEMCalDigits.cxx \ AliJetDummyGeo.cxx AliJetDummyShishKebabTrd1Module.cxx\ AliJetFinderTypes.cxx \ AliAnalysisTaskJets.cxx \ -- 2.43.0