From ee7de0dd4cd0478eb8d46cf642ff78304ecf8f51 Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 13 Sep 2007 08:12:57 +0000 Subject: [PATCH] - Use UnitArray to store track and emcal information. - Hadron correction. (Magali Estienne) --- JETAN/AliJetAODReader.cxx | 2 + JETAN/AliJetControlPlots.h | 29 +- JETAN/AliJetDummyGeo.cxx | 571 ++++++++++++++- JETAN/AliJetDummyGeo.h | 120 +++- JETAN/AliJetESDReader.cxx | 381 ++++++++-- JETAN/AliJetESDReader.h | 19 +- JETAN/AliJetFillUnitArrayEMCalDigits.cxx | 248 +++++++ JETAN/AliJetFillUnitArrayEMCalDigits.h | 95 +++ JETAN/AliJetFillUnitArrayTracks.cxx | 831 ++++++++++----------- JETAN/AliJetFillUnitArrayTracks.h | 63 +- JETAN/AliJetFinder.cxx | 2 - JETAN/AliJetHadronCorrectionv1.cxx | 1 - JETAN/AliJetHadronCorrectionv1.h | 1 - JETAN/AliJetMCReader.cxx | 1 + JETAN/AliJetReader.cxx | 20 +- JETAN/AliJetReader.h | 53 +- JETAN/AliJetReaderHeader.cxx | 8 + JETAN/AliJetReaderHeader.h | 13 +- JETAN/AliJetUnitArray.cxx | 65 +- JETAN/AliJetUnitArray.h | 58 +- JETAN/AliUA1JetFinderV2.cxx | 878 +++++++++++++++++++++++ JETAN/AliUA1JetFinderV2.h | 69 ++ JETAN/AliUA1JetHeaderV1.cxx | 1 + JETAN/AliUA1JetHeaderV1.h | 3 + JETAN/JETANLinkDef.h | 2 + JETAN/libJETAN.pkg | 4 +- 26 files changed, 2946 insertions(+), 592 deletions(-) create mode 100644 JETAN/AliJetFillUnitArrayEMCalDigits.cxx create mode 100644 JETAN/AliJetFillUnitArrayEMCalDigits.h create mode 100644 JETAN/AliUA1JetFinderV2.cxx create mode 100644 JETAN/AliUA1JetFinderV2.h diff --git a/JETAN/AliJetAODReader.cxx b/JETAN/AliJetAODReader.cxx index 592ac61f9e9..bb286315b29 100644 --- a/JETAN/AliJetAODReader.cxx +++ b/JETAN/AliJetAODReader.cxx @@ -24,6 +24,8 @@ #include #include #include +#include + #include "AliJetAODReader.h" #include "AliJetAODReaderHeader.h" #include "AliAODEvent.h" diff --git a/JETAN/AliJetControlPlots.h b/JETAN/AliJetControlPlots.h index d1e305ff84b..48c6bdad6e4 100755 --- a/JETAN/AliJetControlPlots.h +++ b/JETAN/AliJetControlPlots.h @@ -18,6 +18,9 @@ class TFile; class TClonesArray; class TH1I; class TH1D; +class TH1F; +class TH1; + class AliJetReader; class AliJet; @@ -29,19 +32,19 @@ class AliJetControlPlots : public TObject // setter // getters - TH1I *GetNJetsH() {return fNJetsH;} - TH1I *GetMultH() {return fMultH;} - TH1D *GetPhiH() {return fPhiH;} - TH1D *GetFractionInJetH() {return fInJetH;} - TH1D *GetEneH() {return fEneH;} - TH1D *GetPtH() {return fPtH;} - TH1D *GetEtaH() {return fEtaH;} - TH1D *GetFragH() {return fFragH;} - TH1D *GetFragLnH() {return fFragLnH;} - TH1D *GetFragrH() {return fFragrH;} - TH1D *GetFragLnrH() {return fFragLnrH;} - TH1D *GetShapeH() {return fShapeH;} - TH1D *GetShaperH() {return fShaperH;} + TH1I *GetNJetsH() const {return fNJetsH;} + TH1I *GetMultH() const {return fMultH;} + TH1D *GetPhiH() const {return fPhiH;} + TH1D *GetFractionInJetH() const {return fInJetH;} + TH1D *GetEneH() const {return fEneH;} + TH1D *GetPtH() const {return fPtH;} + TH1D *GetEtaH() const {return fEtaH;} + TH1D *GetFragH() const {return fFragH;} + TH1D *GetFragLnH() const {return fFragLnH;} + TH1D *GetFragrH() const {return fFragrH;} + TH1D *GetFragLnrH() const {return fFragLnrH;} + TH1D *GetShapeH() const {return fShapeH;} + TH1D *GetShaperH() const {return fShaperH;} // others void FillHistos(AliJet *j); diff --git a/JETAN/AliJetDummyGeo.cxx b/JETAN/AliJetDummyGeo.cxx index b2a1700d77b..7e5010e3ef8 100644 --- a/JETAN/AliJetDummyGeo.cxx +++ b/JETAN/AliJetDummyGeo.cxx @@ -12,7 +12,576 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - + +// +// Temporarily added to define part of the EMCal geometry +// necessary for the jet finder +// + +#include + +// --- Root header files --- +#include +#include +#include +#include + +// --- AliRoot header files --- #include "AliJetDummyGeo.h" ClassImp(AliJetDummyGeo) + +AliJetDummyGeo::AliJetDummyGeo(): + TObject(), + fArm1PhiMin(80.0), + fArm1PhiMax(200.0), + fArm1EtaMin(-0.7), + fArm1EtaMax(+0.7), + fNumberOfSuperModules(12), + fSteelFrontThick(0.0), + fIPDistance(428.0), + fZLength(), + fPhiGapForSM(2.), + fNPhi(12), + fNZ(24), + fPhiModuleSize(12.26 - fPhiGapForSM / Float_t(fNPhi)), + fEtaModuleSize(fPhiModuleSize), + fNPHIdiv(2), + fNETAdiv(2), + fNECLayers(77), + fECScintThick(0.16), + fECPbRadThickness(0.16), + fSampling(12.327), + fTrd1Angle(1.5), + fNCellsInModule(fNPHIdiv*fNETAdiv), + fNCellsInSupMod(fNCellsInModule*fNPhi*fNZ), + fNCells(fNCellsInSupMod*fNumberOfSuperModules-fNCellsInSupMod), + fLongModuleSize(fNECLayers*(fECScintThick + fECPbRadThickness)), + f2Trd1Dx2(fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.)), + fShellThickness(TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2)+fSteelFrontThick), + fEtaMaxOfTRD1(0.67064) // Value extracted from ShishKebab +{ + // Constructor + // Local coordinates + fParSM[0] = GetShellThickness()/2.; + fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.; + fParSM[2] = 350./2.; + + fZLength = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage + fEnvelop[0] = fIPDistance; // mother volume inner radius + fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r. + fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume. + + // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 + fPhiBoundariesOfSM.Set(fNumberOfSuperModules); + fPhiCentersOfSM.Set(fNumberOfSuperModules/2); + fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules) + fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance); + fPhiCentersOfSM[0] = TMath::PiOver2(); + for(int i=1; i<=4; i++) { // from 2th ro 9th + fPhiBoundariesOfSM[2*i] = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i; + fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i; + fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i; + } + fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad(); + fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance); + fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; + + // for(int i=0; i=360.) phiy -= 360.; + TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0); + fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind), + xpos,ypos,-zpos, geoRot1); + } // for + +} + +AliJetDummyGeo::AliJetDummyGeo(const AliJetDummyGeo& geom): + TObject(), + fArm1PhiMin(geom.fArm1PhiMin), + fArm1PhiMax(geom.fArm1PhiMax), + fArm1EtaMin(geom.fArm1EtaMin), + fArm1EtaMax(geom.fArm1EtaMax), + fNumberOfSuperModules(geom.fNumberOfSuperModules), + fSteelFrontThick(geom.fSteelFrontThick), + fIPDistance(geom.fIPDistance), + fPhiGapForSM(geom.fPhiGapForSM), + fNPhi(geom.fNPhi), + fNZ(geom.fNZ), + fPhiModuleSize(geom.fPhiModuleSize), + fEtaModuleSize(geom.fEtaModuleSize), + fNPHIdiv(geom.fNPHIdiv), + fNETAdiv(geom.fNETAdiv), + fNECLayers(geom.fNECLayers), + fECScintThick(geom.fECScintThick), + fECPbRadThickness(geom.fECPbRadThickness), + fSampling(geom.fSampling), + fTrd1Angle(geom.fTrd1Angle), + fNCellsInModule(geom.fNCellsInModule), + fNCellsInSupMod(geom.fNCellsInSupMod), + fNCells(geom.fNCells), + fLongModuleSize(geom.fLongModuleSize), + f2Trd1Dx2(geom.f2Trd1Dx2), + fShellThickness(geom.fShellThickness), + fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1) +{ + // Constructor + // Local coordinates + fParSM[0] = GetShellThickness()/2.; + fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.; + fParSM[2] = 350./2.; + + // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 + fPhiBoundariesOfSM.Set(fNumberOfSuperModules); + fPhiCentersOfSM.Set(fNumberOfSuperModules/2); + fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules) + fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance); + fPhiCentersOfSM[0] = TMath::PiOver2(); + for(int i=1; i<=4; i++) { // from 2th ro 9th + fPhiBoundariesOfSM[2*i] = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i; + fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i; + fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i; + } + fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad(); + fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance); + fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; + + // for(int i=0; i=360.) phiy -= 360.; + TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0); + fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind), + xpos,ypos,-zpos, geoRot1); + + delete geoRot0; + delete geoRot1; + + } // for + +} + +AliJetDummyGeo::~AliJetDummyGeo() +{ + // Destructor + delete [] fMatrixOfSM; +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::EtaPhiFromIndex(Int_t absId, Float_t& eta, Float_t& phi) +{ + // Nov 16, 2006- float to double + // version for TRD1 only + static TVector3 vglob; + GetGlobal(absId, vglob); + eta = vglob.Eta(); + phi = vglob.Phi(); + +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const +{ + // Figure out the global numbering of a given supermodule from the + // local numbering + // Alice numbering - Jun 03,2006 + // if(fMatrixOfSM[0] == 0) GetTransformationForSM(); + + if(ind>=0 && ind < GetNumberOfSuperModules()) { + fMatrixOfSM[ind]->LocalToMaster(loc, glob); + } +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetGlobal(Int_t absId , double glob[3]) const +{ + // Alice numbering scheme - Jun 03, 2006 + static Int_t nSupMod, nModule, nIphi, nIeta; + static double loc[3]; + + glob[0]=glob[1]=glob[2]=0.0; // bad case + if(RelPosCellInSModule(absId, loc)) { + GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta); + fMatrixOfSM[nSupMod]->LocalToMaster(loc, glob); + } +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetGlobal(Int_t absId , TVector3 &vglob) const +{ + // Alice numbering scheme - Jun 03, 2006 + static Double_t glob[3]; + + GetGlobal(absId, glob); + vglob.SetXYZ(glob[0], glob[1], glob[2]); + +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const +{ + // Alice numbering scheme - Jun 03, 2006 + loc[0] = loc[1] = loc[2]=0.0; + if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) { + return kTRUE; + } + return kFALSE; +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const +{ + // Look to see what the relative position inside a given cell is + // for a recpoint. + // Alice numbering scheme - Jun 08, 2006 + // In: + // absId - cell is as in Geant, 0<= absId < fNCells; + // OUT: + // xr,yr,zr - x,y,z coordinates of cell with absId inside SM + + // Shift index taking into account the difference between standard SM + // and SM of half size in phi direction + const Int_t phiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2 + static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta; + if(!CheckAbsCellId(absId)) return kFALSE; + + GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta); + GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); + + xr = fCentersOfCellsXDir.At(ieta); + zr = fCentersOfCellsEtaDir.At(ieta); + + if(nSupMod<10) { + yr = fCentersOfCellsPhiDir.At(iphi); + } else { + yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift); + } + + return kTRUE; +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::CheckAbsCellId(Int_t absId) const +{ + // May 31, 2006; only trd1 now + if(absId<0 || absId >= fNCells) return kFALSE; + else return kTRUE; +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const +{ + // 21-sep-04; 19-oct-05; + // May 31, 2006; ALICE numbering scheme: + // + // In: + // absId - cell is as in Geant, 0<= absId < fNCells; + // Out: + // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules; + // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th); + // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; + // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; + // + static Int_t tmp=0, sm10=0; + if(!CheckAbsCellId(absId)) return kFALSE; + + sm10 = fNCellsInSupMod*10; + if(absId >= sm10) { // 110 degree case; last two supermodules + nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10; + tmp = (absId-sm10) % (fNCellsInSupMod/2); + } else { + nSupMod = absId / fNCellsInSupMod; + tmp = absId % fNCellsInSupMod; + } + + nModule = tmp / fNCellsInModule; + tmp = tmp % fNCellsInModule; + nIphi = tmp / fNPHIdiv; + nIeta = tmp % fNPHIdiv; + + return kTRUE; +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, int &iphi, int &ieta) const +{ + // + // Added nSupMod; Nov 25, 05 + // Alice numbering scheme - Jun 01,2006 + // IN: + // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules; + // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th); + // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; + // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; + // + // OUT: + // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM + // ieta - have to change from 0 to (fNZ*fNETAdiv-1) + // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1) + // + static Int_t iphim, ietam; + + GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam); + // ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) + ieta = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM) + iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM) + +} + + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, int &iphim, int &ietam) const +{ + // added nSupMod; - 19-oct-05 ! + // Alice numbering scheme - Jun 01,2006 + // ietam, iphi - indexes of module in two dimensional grid of SM + // ietam - have to change from 0 to fNZ-1 + // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1) + static Int_t nphi; + + if(nSupMod>=10) nphi = fNPhi/2; + else nphi = fNPhi; + + ietam = nModule/nphi; + iphim = nModule%nphi; +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const +{ + // Nov 17,2006 + // stay here - phi problem as usual + static Int_t nSupMod, i, ieta, iphi, etaShift, nphi; + static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc; + absId = nSupMod = - 1; + if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) { + // phi index first + phi = TVector2::Phi_0_2pi(phi); + phiLoc = phi - fPhiCentersOfSM[nSupMod/2]; + nphi = fPhiCentersOfCells.GetSize(); + if(nSupMod>=10) { + phiLoc = phi - 190.*TMath::DegToRad(); + nphi /= 2; + } + + dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc); + iphi = 0; + for(i=1; i fEtaMaxOfTRD1) return kFALSE; + + phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries + for(i=0; i<6; i++) { + if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) { + nSupMod = 2*i; + if(eta < 0.0) nSupMod++; + return kTRUE; + } + } + return kFALSE; +} + +//------------------------------------------------------------------------------------ +Int_t AliJetDummyGeo::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const +{ + // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId + static Int_t ietam, iphim, nModule; + static Int_t nIeta, nIphi; // cell indexes in module + + GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule); + + nIeta = ieta%fNETAdiv; + nIeta = fNETAdiv - 1 - nIeta; + nIphi = iphi%fNPHIdiv; + + return GetAbsCellId(nSupMod, nModule, nIphi, nIeta); +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, + Int_t &iphim, Int_t &ietam, Int_t &nModule) const +{ + // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule) + static Int_t nphi; + nphi = GetNumberOfModuleInPhiDirection(nSupMod); + + ietam = ieta/fNETAdiv; + iphim = iphi/fNPHIdiv; + nModule = ietam * nphi + iphim; +} + +//------------------------------------------------------------------------------------ +Int_t AliJetDummyGeo::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const +{ + // 27-aug-04; + // corr. 21-sep-04; + // 13-oct-05; 110 degree case + // May 31, 2006; ALICE numbering scheme: + // 0 <= nSupMod < fNumberOfSuperModules + // 0 <= nModule < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1) + // 0 <= nIphi < fNPHIdiv + // 0 <= nIeta < fNETAdiv + // 0 <= absid < fNCells + static Int_t id=0; // have to change from 0 to fNCells-1 + if(nSupMod >= 10) { // 110 degree case; last two supermodules + id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10); + } else { + id = fNCellsInSupMod*nSupMod; + } + id += fNCellsInModule *nModule; + id += fNPHIdiv *nIphi; + id += nIeta; + if(id<0 || id >= fNCells) { + id = -TMath::Abs(id); // if negative something wrong + } + return id; +} + +//------------------------------------------------------------------------------------ +Bool_t AliJetDummyGeo::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const +{ + // 0<= nPhiSec <=4; phi in rad + // 0; gap boundaries between 0th&2th | 1th&3th SM + // 1; gap boundaries between 2th&4th | 3th&5th SM + // 2; gap boundaries between 4th&6th | 5th&7th SM + // 3; gap boundaries between 6th&8th | 7th&9th SM + // 4; gap boundaries between 8th&10th | 9th&11th SM + if(nPhiSec<0 || nPhiSec >4) return kFALSE; + phiMin = fPhiBoundariesOfSM[2*nPhiSec+1]; + phiMax = fPhiBoundariesOfSM[2*nPhiSec+2]; + return kTRUE; +} + +//------------------------------------------------------------------------------------ +void AliJetDummyGeo::GetTransformationForSM() +{ + // Uses the geometry manager to load the transformation matrix + // for the supermodules + // Unused after 19 Jan, 2007 - keep for compatibility; + + return; + static Bool_t transInit=kFALSE; + if(transInit) return; + + int i=0; + if(gGeoManager == 0) { + Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()"); + assert(0); + } + + TGeoNode *tn = gGeoManager->GetTopNode(); + TGeoNode *node=0, *xen1 = 0; + for(i=0; iGetNdaughters(); i++) { + node = tn->GetDaughter(i); + TString ns(node->GetName()); + if(ns.Contains(GetNameOfEMCALEnvelope())) { + xen1 = node; + break; + } + } + + if(!xen1) { + Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s", + GetNameOfEMCALEnvelope()); + assert(0); + } + printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters()); + for(i=0; iGetNdaughters(); i++) { + TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i); + fMatrixOfSM[i] = sm->GetMatrix(); + } + printf("transInit %d: ", transInit); + transInit = kTRUE; +} + diff --git a/JETAN/AliJetDummyGeo.h b/JETAN/AliJetDummyGeo.h index 475d5ad9fbc..2d29bddf7c2 100644 --- a/JETAN/AliJetDummyGeo.h +++ b/JETAN/AliJetDummyGeo.h @@ -4,27 +4,115 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ +// +// Temporarily added to define part of the EMCal geometry +// necessary for the jet finder +// + #include +#include +#include +#include + +class TGeoMatrix; class AliJetDummyGeo : public TObject { public: - AliJetDummyGeo(){;} - virtual ~AliJetDummyGeo(){;} - static AliJetDummyGeo* GetInstance() {return new AliJetDummyGeo();} - static AliJetDummyGeo* GetInstance(char* /*name*/, char* /*title*/) - {return new AliJetDummyGeo();} - Int_t GetNCells(){ return 0;} - Float_t GetArm1EtaMin() {return 0.;} - Float_t GetArm1EtaMax() {return 0.;} - Float_t GetArm1PhiMin() {return 0.;} - Float_t GetArm1PhiMax() {return 0.;} - void EtaPhiFromIndex(Int_t /*id*/, Float_t& /*eta*/, Float_t& /*phi*/) - {;} - Int_t TowerIndexFromEtaPhi2(Float_t /*eta*/, Float_t /*phi*/) {return 0;} - void GetTransformationForSM(){;} - Float_t GetSampling() {return 0.;} - ClassDef(AliJetDummyGeo,1) + AliJetDummyGeo(); + AliJetDummyGeo(const AliJetDummyGeo& geom); + virtual ~AliJetDummyGeo(); + static AliJetDummyGeo* GetInstance() {return new AliJetDummyGeo();} + static AliJetDummyGeo* GetInstance(char* /*name*/, char* /*title*/) + {return new AliJetDummyGeo();} + Char_t* GetNameOfEMCALEnvelope() const {return "XEN1";} + Float_t GetEnvelop(Int_t index) const { return fEnvelop[index];} + Float_t AngleFromEta(Float_t eta) const { + // returns theta in radians for a given pseudorapidity + return 2.0*TMath::ATan(TMath::Exp(-eta)); + } + Float_t ZFromEtaR(Float_t r,Float_t eta) const { + // returns z in for a given + // pseudorapidity and r=sqrt(x*x+y*y). + return r/TMath::Tan(AngleFromEta(eta)); + } + Int_t GetNCells() {return fNCells;} + Float_t GetPhiModuleSize() const {return fPhiModuleSize;} + Float_t GetEtaModuleSize() const {return fEtaModuleSize;} + Float_t GetShellThickness() const {return fShellThickness;} + Int_t GetNPhi() const {return fNPhi;} + Int_t GetNumberOfSuperModules() const {return fNumberOfSuperModules;} + Float_t GetArm1EtaMin() {return fArm1EtaMin;} + Float_t GetArm1EtaMax() {return fArm1EtaMax;} + Float_t GetArm1PhiMin() {return fArm1PhiMin;} + Float_t GetArm1PhiMax() {return fArm1PhiMax;} + void EtaPhiFromIndex(Int_t /*id*/, Float_t& /*eta*/, Float_t& /*phi*/); + void GetGlobal(const Double_t *loc, Double_t *glob, int ind) const; + void GetGlobal(Int_t absId, Double_t glob[3]) const; + void GetGlobal(Int_t absId, TVector3 &vglob) const; + Bool_t RelPosCellInSModule(Int_t absId, Double_t loc[3]) const; + Bool_t RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const; + Bool_t CheckAbsCellId(Int_t absId) const; + Bool_t GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const; + void GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, int &iphi, int &ieta) const; + void GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, int &iphim, int &ietam) const; + Bool_t GetAbsCellIdFromEtaPhi(Double_t eta,Double_t phi, Int_t &absId) const; + Bool_t SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const; + Int_t GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const; + void GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, + Int_t &iphim, Int_t &ietam, Int_t &nModule) const; + Int_t GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const; + Int_t GetNumberOfModuleInPhiDirection(Int_t nSupMod) const + { + // inline function + if(nSupMod>=10) return fNPhi/2; + else return fNPhi; + } + Bool_t GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const; + void GetTransformationForSM(); + Float_t GetSampling() {return fSampling;} + + protected: + Float_t fArm1PhiMin; // Minimum angular position of EMCAL in Phi (degrees) + Float_t fArm1PhiMax; // Maximum angular position of EMCAL in Phi (degrees) + Float_t fArm1EtaMin; // Minimum pseudorapidity position of EMCAL in Eta + Float_t fArm1EtaMax; // Maximum pseudorapidity position of EMCAL in Eta + Int_t fNumberOfSuperModules; + Float_t fSteelFrontThick; // Thickness of the front stell face of the support box - 9-sep-04 + Float_t fEnvelop[3]; // the GEANT TUB for the detector + Float_t fIPDistance; // Radial Distance of the inner surface of the EMCAL + Float_t fZLength; // Total length in z direction + Float_t fPhiGapForSM; // Gap betweeen supermodules in phi direction + Int_t fNPhi; // Number of Towers in the PHI direction + Int_t fNZ; // Number of Towers in the Z direction + Float_t fPhiModuleSize; // Phi -> X + Float_t fEtaModuleSize; // Eta -> Y + Int_t fNPHIdiv; // number phi divizion of module + Int_t fNETAdiv; // number eta divizion of module + Int_t fNECLayers; // number of scintillator layers + Float_t fECScintThick; // cm, Thickness of the scintillators + Float_t fECPbRadThickness; // cm, Thickness of the Pb radiators + Float_t fSampling; // Sampling factor + Float_t fTrd1Angle; // angle in x-z plane (in degree) + Int_t fNCellsInModule; // number cell in module + Int_t fNCellsInSupMod; // number cell in super module + Int_t fNCells; // number of cells in calo + Float_t fLongModuleSize; // Size of long module + Float_t f2Trd1Dx2; // 2*dx2 for TRD1 + Float_t fShellThickness; // Total thickness in (x,y) direction + Float_t fEtaMaxOfTRD1; // max eta in case of TRD1 geometry (see AliEMCALShishKebabTrd1Module) + Float_t fParSM[3]; // SM sizes as in GEANT (TRD1) + TArrayD fPhiBoundariesOfSM; // phi boundaries of SM in rad; size is fNumberOfSuperModules; + TArrayD fPhiCentersOfSM; // phi of centers of SMl size is fNumberOfSuperModules/2 + TGeoMatrix* fMatrixOfSM[12]; //![fNumberOfSuperModules]; get from gGeoManager; + TArrayD fCentersOfCellsEtaDir; // size fNEta*fNETAdiv (for TRD1 only) (eta or z in SM, in cm) + TArrayD fCentersOfCellsXDir; // size fNEta*fNETAdiv (for TRD1 only) ( x in SM, in cm) + TArrayD fCentersOfCellsPhiDir; // size fNPhi*fNPHIdiv (for TRD1 only) (phi or y in SM, in cm) + TArrayD fEtaCentersOfCells; // [fNEta*fNETAdiv*fNPhi*fNPHIdiv], positive direction (eta>0); + // eta depend from phi position; + TArrayD fPhiCentersOfCells; // [fNPhi*fNPHIdiv] from center of SM (-10. < phi < +10.) + + ClassDef(AliJetDummyGeo,1) }; #endif diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 3975f25e876..733b26ac207 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -20,48 +20,65 @@ // Magali Estienne //------------------------------------------------------------------------- - +// --- Standard library --- #include + +// --- ROOT system --- #include +#include #include #include +#include #include +#include +#include +#include +#include + +// --- AliRoot header files --- #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" #include "AliESDEvent.h" #include "AliESD.h" #include "AliESDtrack.h" -//#include "AliEMCALGeometry.h" #include "AliJetDummyGeo.h" #include "AliJetFillUnitArrayTracks.h" +#include "AliJetFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" ClassImp(AliJetESDReader) AliJetESDReader::AliJetESDReader(): - AliJetReader(), - fGeom(0), - fChain(0x0), - fESD(0x0), - fHadCorr(0x0), - fTpcGrid(0x0), - fEmcalGrid(0x0), - fPtCut(0), - fHCorrection(0), - fNumUnits(0), - fDebug(0), - fNIn(0), - fOpt(0), - fNeta(0), - fNphi(0), - fArrayInitialised(0) + AliJetReader(), + fGeom(0), + fChain(0x0), + fESD(0x0), + fHadCorr(0x0), + fTpcGrid(0x0), + fEmcalGrid(0x0), + fGrid0(0), + fGrid1(0), + fGrid2(0), + fGrid3(0), + fGrid4(0), + fPtCut(0), + fHCorrection(0), + fNumUnits(0), + fDebug(0), + fMass(0), + fSign(0), + fNIn(0), + fOpt(0), + fDZ(0), + fNeta(0), + fNphi(0), + fArrayInitialised(0) { - // Constructor + // Constructor } //____________________________________________________________________________ - AliJetESDReader::~AliJetESDReader() { // Destructor @@ -69,10 +86,17 @@ AliJetESDReader::~AliJetESDReader() delete fESD; delete fTpcGrid; delete fEmcalGrid; + if(fDZ) + { + delete fGrid0; + delete fGrid1; + delete fGrid2; + delete fGrid3; + delete fGrid4; + } } //____________________________________________________________________________ - void AliJetESDReader::OpenInputFiles() { // chain for the ESDs @@ -82,8 +106,7 @@ void AliJetESDReader::OpenInputFiles() const char* dirName=fReaderHeader->GetDirectory(); const char* pattern=fReaderHeader->GetPattern(); -// // Add files matching patters to the chain - + // Add files matching patters to the chain void *dir = gSystem->OpenDirectory(dirName); const char *name = 0x0; int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd(); @@ -92,19 +115,20 @@ void AliJetESDReader::OpenInputFiles() if (a>=nesd) continue; if (strstr(name,pattern)){ - char path[256]; - sprintf(path,"%s/%s/AliESDs.root",dirName,name); - fChain->AddFile(path); - a++; + printf("Adding %s\n",name); + char path[256]; + // sprintf(path,"%s/%s/AliESDs.root",dirName,name); + sprintf(path,"%s/%s/AliESDs.root",dirName,name); + fChain->AddFile(path); + a++; } } gSystem->FreeDirectory(dir); + fESD = new AliESDEvent(); + fESD->ReadFromTree(fChain); - fESD = 0; - fChain->SetBranchAddress("ESD", &fESD); - int nMax = fChain->GetEntries(); printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax); @@ -118,13 +142,13 @@ void AliJetESDReader::OpenInputFiles() } } +//____________________________________________________________________________ void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) { // Connect the tree fESD = (AliESDEvent*) esd; } //____________________________________________________________________________ - Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/) { // Fill momentum array @@ -137,11 +161,11 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/) // clear momentum array ClearArray(); fDebug = fReaderHeader->GetDebug(); - + fOpt = fReaderHeader->GetDetector(); + if (!fESD) { return kFALSE; } - // get number of tracks in event (for the loop) nt = fESD->GetNumberOfTracks(); @@ -184,50 +208,138 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/) // set the signal flags fSignalFlag.Set(goodTrack,sflag); fCutFlag.Set(goodTrack,cflag); + return kTRUE; -// -// - if (fTpcGrid || fEmcalGrid) { - SetEMCALGeometry(); - InitParameters(); - AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks(); - fillUAFromTracks->SetReaderHeader(fReaderHeader); - fillUAFromTracks->SetMomentumArray(fMomentumArray); - fillUAFromTracks->SetTPCGrid(fTpcGrid); - fillUAFromTracks->SetEMCalGrid(fEmcalGrid); - fillUAFromTracks->SetHadCorrection(fHCorrection); - fillUAFromTracks->SetHadCorrector(fHadCorr); - fNeta = fillUAFromTracks->GetNeta(); - fNphi = fillUAFromTracks->GetNphi(); - fillUAFromTracks->SetActive(kFALSE); - // TPC only or Digits+TPC or Clusters+TPC - if(fOpt%2==!0 && fOpt!=0) { - fillUAFromTracks->SetActive(kTRUE); - fillUAFromTracks->SetUnitArray(fUnitArray); - fillUAFromTracks->ExecuteTask("tpc"); - } +} + +//____________________________________________________________________________ +void AliJetESDReader::CreateTasks() +{ + fDebug = fReaderHeader->GetDebug(); + fDZ = fReaderHeader->GetDZ(); + + // Init EMCAL geometry and create UnitArray object + SetEMCALGeometry(); + InitParameters(); + InitUnitArray(); + + fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); + fFillUAFromTracks = new AliJetFillUnitArrayTracks(); + 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->SetHadCorrection(fHCorrection); + fFillUAFromTracks->SetHadCorrector(fHadCorr); + fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits(); + fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); + fFillUAFromEMCalDigits->SetGeom(fGeom); + fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); + fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); + // 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; + +} + +//____________________________________________________________________________ +//void AliJetESDReader::ExecTasks(Int_t event) +Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/) +{ + // clear momentum array + Int_t nEntRef = fRefArray->GetEntries(); + + for(Int_t i=0; iAt(i))->SetUnitTrackID(0); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc); + ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet); + } + + ClearArray(); + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + // InitParameters(); + + if(!fESD) { + return kFALSE; + } - delete fillUAFromTracks; + /* + // get event from chain + // For TSelectors + // fChain->GetTree()->GetEntry(event); + // For interactive process + // fChain->GetEntry(event); + fChain->GetEvent(event); + */ + + // TPC only or Digits+TPC or Clusters+TPC + if(fOpt%2==!0 && fOpt!=0){ + fFillUAFromTracks->SetESD(fESD); + fFillUAFromTracks->SetActive(kTRUE); + fFillUAFromTracks->SetUnitArray(fUnitArray); + fFillUAFromTracks->SetRefArray(fRefArray); + // fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed !!! + // Incompatibility with Andrei's analysis framework + fFillUAFromTracks->Exec("tpc"); + if(fOpt==1){ + fNumCandidate = fFillUAFromTracks->GetMult(); + fNumCandidateCut = fFillUAFromTracks->GetMultCut(); + } + } + + // Digits only or Digits+TPC + if(fOpt>=2 && fOpt<=3){ + fFillUAFromEMCalDigits->SetESD(fESD); + fFillUAFromEMCalDigits->SetActive(kTRUE); + fFillUAFromEMCalDigits->SetUnitArray(fUnitArray); + fFillUAFromEMCalDigits->SetRefArray(fRefArray); + fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult()); + fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut()); + fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!! + fNumCandidate = fFillUAFromEMCalDigits->GetMult(); + fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut(); } + // fFillUnitArray->ExecuteTask(); // => Temporarily commented + return kTRUE; } - +//____________________________________________________________________________ void AliJetESDReader::SetEMCALGeometry() { // Define EMCAL geometry to be able to read ESDs fGeom = AliJetDummyGeo::GetInstance(); if (fGeom == 0) fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL"); - + // To be setted to run some AliEMCALGeometry functions TGeoManager::Import("geometry.root"); - fGeom->GetTransformationForSM(); + // fGeom->GetTransformationForSM(); printf("\n EMCal Geometry set ! \n"); } - + + +//____________________________________________________________________________ void AliJetESDReader::InitParameters() { // Initialise parameters @@ -237,16 +349,153 @@ void AliJetESDReader::InitParameters() if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); } +//____________________________________________________________________________ void AliJetESDReader::InitUnitArray() { //Initialises unit arrays - Int_t nElements = fTpcGrid->GetNEntries(); - if(fArrayInitialised) delete [] fUnitArray; - if(fTpcGrid->GetGridType()==0) - fUnitArray = new AliJetUnitArray[nElements]; - if(fTpcGrid->GetGridType()==1) - fUnitArray = new AliJetUnitArray[fNumUnits+nElements]; - fArrayInitialised = 1; + 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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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(); + + 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; + + } + + 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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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/AliJetESDReader.h b/JETAN/AliJetESDReader.h index d5cdf59f3cd..55d41cb87fd 100755 --- a/JETAN/AliJetESDReader.h +++ b/JETAN/AliJetESDReader.h @@ -18,7 +18,6 @@ #include "AliJetUnitArray.h" #include "AliJetGrid.h" class AliJetESDReaderHeader; -class AliEMCALGeometry; class AliJetDummyGeo; class AliJetHadronCorrection; class AliJetUnitArray; @@ -30,10 +29,18 @@ class AliJetESDReader : public AliJetReader public: AliJetESDReader(); virtual ~AliJetESDReader(); + + // Getters + Float_t GetTrackMass() const {return fMass;} // returns mass of the track + Int_t GetTrackSign() const {return fSign;} // returns sign of the track + // Setters Bool_t FillMomentumArray(Int_t event); void OpenInputFiles(); void InitUnitArray(); + void CreateTasks(); + // void ExecTasks(Int_t event); + Bool_t ExecTasks(Int_t event); void SetInputEvent(TObject* esd, TObject* aod, TObject* mc); virtual void SetTPCGrid(AliJetGrid *grid) {fTpcGrid = grid;} virtual void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;} @@ -50,19 +57,25 @@ class AliJetESDReader : public AliJetReader AliJetHadronCorrectionv1 *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 fNumUnits; // Number of units in the unit object array // (same as num towers in EMCAL) Int_t fDebug; // Debug option + Float_t fMass; // Particle mass + Int_t fSign; // Particle sign Int_t fNIn; // Number of Array filled in UnitArray Int_t fOpt; // Detector to be used for jet reconstruction + 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 - - ClassDef(AliJetESDReader,1) }; diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetFillUnitArrayEMCalDigits.cxx new file mode 100644 index 00000000000..a07836489ef --- /dev/null +++ b/JETAN/AliJetFillUnitArrayEMCalDigits.cxx @@ -0,0 +1,248 @@ + +/************************************************************************** + * 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 with EMCal information +// Called by ESD reader for jet analysis +// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) +//------------------------------------------------------------- + +// --- Standard library --- +#include + +// --- ROOT system --- +#include +#include +#include +#include "TTask.h" +#include +#include +#include +#include + +// --- AliRoot header files --- +#include "AliJetFinder.h" +#include "AliJetReaderHeader.h" +#include "AliJetReader.h" +#include "AliJetESDReader.h" +#include "AliJetESDReaderHeader.h" +//#include "AliESD.h" +#include "AliESDEvent.h" +#include "AliJetDummyGeo.h" +#include "AliESDCaloCluster.h" +#include "AliJetUnitArray.h" +#include "AliJetFillUnitArrayEMCalDigits.h" +#include "AliCDBManager.h" +#include "AliCDBStorage.h" +#include "AliCDBEntry.h" + +ClassImp(AliJetFillUnitArrayEMCalDigits) + +//_____________________________________________________________________________ +AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits() + : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"), + fESD(0), + fNumUnits(0), + fEtaMinCal(0.), + fEtaMaxCal(0.), + fPhiMinCal(0.), + fPhiMaxCal(0.), + fNIn(0), + fOpt(0), + fDebug(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fTPCGrid(0x0), + fEMCalGrid(0x0), + fReaderHeader(0x0), + fMomentumArray(0x0), + fUnitArray(0x0), + fRefArray(0x0), + fGeom(0x0), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + +//_____________________________________________________________________________ +AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/) + : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"), + fESD(0), + fNumUnits(0), + fEtaMinCal(0.), + fEtaMaxCal(0.), + fPhiMinCal(0.), + fPhiMaxCal(0.), + fNIn(0), + fOpt(0), + fDebug(0), + fNCEMCAL(0), + fNCPHOS(0), + fNCCalo(0), + fTPCGrid(0x0), + fEMCalGrid(0x0), + fReaderHeader(0x0), + fMomentumArray(0x0), + fUnitArray(0x0), + fRefArray(0x0), + fGeom(0x0), + fClus(0x0), + fNDigitEmcal(0), + fNDigitEmcalCut(0) +{ + // constructor +} + + +//____________________________________________________________________________ +void AliJetFillUnitArrayEMCalDigits::InitParameters() +{ + fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL + + fEtaMinCal = fGeom->GetArm1EtaMin(); + fEtaMaxCal = fGeom->GetArm1EtaMax(); + fPhiMinCal = fGeom->GetArm1PhiMin(); + fPhiMaxCal = fGeom->GetArm1PhiMax(); + fClus = 0; + + if(fDebug>1) printf("\n EMCAL parameters initiated ! \n"); + +} + +//_____________________________________________________________________________ +AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/) +{ + // + // Main method. + // Explain + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + + // Init parameters + InitParameters(); + + // Get number of clusters from EMCAL + + Int_t nDigitTot = 0; + Int_t goodDigit = 0; + Int_t beg = 0; + Int_t end = 0; + Float_t ptMin = fReaderHeader->GetPtCut(); + + // Loop over calo clusters + //------------------------------------------------------------------ + Int_t type = 0; + Int_t index = 0; + + // Total number of EMCAL cluster + end = fESD->GetNumberOfCaloClusters(); + + for(Int_t j = beg; j < end; j++) { + fClus = fESD->GetCaloCluster(j); + if(!fClus->IsEMCAL()) continue; + + type = fClus->GetClusterType(); + index = fClus->GetID(); + nDigitTot = fClus->GetNumberOfDigits(); + + // Keep clusters or pseudo clusters + if (type != AliESDCaloCluster::kClusterv1) continue; + // if (type != AliESDCaloCluster::kPseudoCluster) continue; + + // Get the digit index and the digit information + //============================================================ + + // Get number of digits in a cluster + Int_t nD = fClus->GetNumberOfDigits(); + + TArrayS *digID = fClus->GetDigitIndex(); + TArrayS *digEnergy = fClus->GetDigitAmplitude(); + Float_t *digitEnergy = new Float_t[nD]; + // Float_t digitEn = 0.; + + // Loop over digits + for(Int_t k=0; kAt(k); + // Calibration for an energy in GeV + digitEnergy[k] = (Float_t)digEnergy->At(k)/500.; + + // Second method to extract eta, phi positions of a digit + //================================================================= + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(idF,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD); + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idF); + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + if(unitEnergy==0){ + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + uArray->SetUnitEnergy(unitEnergy+etDigit); + // Put a cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) fNDigitEmcalCut++; + } + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 12) printf("goodDigit : %d\n", goodDigit); + + } // End loop over digits + + } // End loop over clusters + + fNIn += goodDigit; + +} + +//_____________________________________________________________________________ +Float_t AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg) +{ + // return (180./TMath::Pi())*2.*atan(exp(-arg)); + return 2.*atan(exp(-arg)); + + +} diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.h b/JETAN/AliJetFillUnitArrayEMCalDigits.h new file mode 100644 index 00000000000..e59332c3a8a --- /dev/null +++ b/JETAN/AliJetFillUnitArrayEMCalDigits.h @@ -0,0 +1,95 @@ +#ifndef ALIJETFILLUNITARRAYEMCALDIGITS_H +#define ALIJETFILLUNITARRAYEMCALDIGITS_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) +//--------------------------------------------------------------------- + +#ifndef ROOT_TTask +#include "TTask.h" +#endif + +class AliJetDummyGeo; +class AliESDCaloCluster; +class AliEMCALCalibData; +class AliJetReader; +class AliJetESDReader; +class TClonesArray; +class TRefArray; +class AliJetUnitArray; +//class AliESD; +class AliESDEvent; +class AliJetGrid; + +class AliJetFillUnitArrayEMCalDigits : public TTask +{ + public: + AliJetFillUnitArrayEMCalDigits(); + AliJetFillUnitArrayEMCalDigits(Int_t event); + AliJetFillUnitArrayEMCalDigits(AliESD *fESD); + AliJetFillUnitArrayEMCalDigits(AliESDEvent *fESD); + virtual ~AliJetFillUnitArrayEMCalDigits(); + + // Setter + void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;} + void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;} + void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;} + void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;} + void SetRefArray(TRefArray *refArray) {fRefArray = refArray;} + void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;} + void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;} + // void SetESD(AliESD *esd) {fESD = esd;} + void SetESD(AliESDEvent *esd) {fESD = esd;} + void SetInitMult(Int_t mult) {fNDigitEmcal = mult;} + void SetInitMultCut(Int_t multcut) {fNDigitEmcalCut = multcut;} + + // Getter + TClonesArray* GetUnitArray() {return fUnitArray;} + TRefArray* GetRefArray() {return fRefArray;} + Int_t GetMult() {return fNDigitEmcal;} + Int_t GetMultCut() {return fNDigitEmcalCut;} + + // Other + void Exec(Option_t*); + Float_t EtaToTheta(Float_t arg); + private: + void InitParameters(); + + protected: + AliESDEvent *fESD; // ESD + Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) + Float_t fEtaMinCal; // Define EMCAL acceptance in Eta + Float_t fEtaMaxCal; // Define EMCAL acceptance in Eta + Float_t fPhiMinCal; // Define EMCAL acceptance in Phi + Float_t fPhiMaxCal; // Define EMCAL acceptance in Phi + Int_t fNIn; // Number of Array filled in UnitArray + Int_t fOpt; // Detector to be used for jet reconstruction + 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 + + AliJetGrid *fTPCGrid; // Define filled grid + AliJetGrid *fEMCalGrid; // Define filled grid + + AliJetReaderHeader *fReaderHeader; // ReaderHeader + TClonesArray *fMomentumArray; // MomentumArray + TClonesArray *fUnitArray; // UnitArray + TRefArray *fRefArray; // UnitArray + AliJetDummyGeo *fGeom; // Set EMCal geometry + + AliESDCaloCluster *fClus; //! + Int_t fNDigitEmcal; //! + Int_t fNDigitEmcalCut; //! + + + + ClassDef(AliJetFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information +}; + +#endif diff --git a/JETAN/AliJetFillUnitArrayTracks.cxx b/JETAN/AliJetFillUnitArrayTracks.cxx index 1cf9d0e2b9b..e5cdc2efe0e 100644 --- a/JETAN/AliJetFillUnitArrayTracks.cxx +++ b/JETAN/AliJetFillUnitArrayTracks.cxx @@ -1,4 +1,3 @@ - /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -14,35 +13,13 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//====================================================================== -// *** November 2006 -// Author: magali.estienne@ires.in2p3.fr -// 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and -// EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance -// 2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment) -// fGrid==1, work with a part of the TPC acceptance -// 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them -// 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry -// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) -//====================================================================== -// ***September 2006 -// TTask : Fill Unit Array for the Tracks information -// Called by ESD reader for jet analysis -// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) -//====================================================================== -// *** July 2006 -// 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex -// and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the -// index or the eta, phi position of a grid. -// 2) Define a grid for TPC -// Author: Magali Estienne (magali.estienne@ires.in2p3.fr) -//====================================================================== -// ***July 2006 + +//---------------------------------------------------------------------- // Fill Unit Array class -// Class used by AliJetESDReader to fill a UnitArray from the information extracted -// from the particle tracks +// Class used by AliJetESDReader to fill a UnitArray from the information +// extracted from the particle tracks // Author: magali.estienne@ires.in2p3.fr -//====================================================================== +//---------------------------------------------------------------------- // --- Standard library --- @@ -51,13 +28,14 @@ // --- ROOT system --- #include #include +#include #include -//#include "Math/Vector3D.h" -//#include "Math/Vector3Dfwd.h" #include "TTask.h" #include #include #include +#include +#include // --- AliRoot header files --- #include "AliJetFinder.h" @@ -65,8 +43,8 @@ #include "AliJetReader.h" #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" -#include "AliESD.h" -//#include "AliEMCALGeometry.h" +//#include "AliESD.h" +#include "AliESDEvent.h" #include "AliJetDummyGeo.h" #include "AliJetUnitArray.h" #include "AliJetFillUnitArrayTracks.h" @@ -77,59 +55,108 @@ ClassImp(AliJetFillUnitArrayTracks) //_____________________________________________________________________________ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks() - : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"), - fNumUnits(0), - fEtaMinCal(0.), - fEtaMaxCal(0.), - fPhiMinCal(0.), - fPhiMaxCal(0.), - fHadCorr(0x0), - fHCorrection(0), - fNIn(0), - fOpt(0), - fDebug(0), - fReaderHeader(0x0), - fMomentumArray(0x0), - fUnitArray(0x0), - fTPCGrid(0x0), - fEMCalGrid(0x0), - fGeom(0x0), - fNphi(0), - fNeta(0), - fPhi2(0x0), - fEta2(0x0), - fPhi(0x0), - fEta(0x0), - fIndex(0x0), - fParams(0x0), - fGrid(0), - fPhiMin(0.), - fPhiMax(0.), - fEtaMin(0.), - fEtaMax(0.), - fEtaBinInTPCAcc(0), - fPhiBinInTPCAcc(0), - fEtaBinInEMCalAcc(0), - fPhiBinInEMCalAcc(0), - fNbinPhi(0) + : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"), + fNumUnits(0), + fEtaMinCal(0), + fEtaMaxCal(0), + fPhiMinCal(0), + fPhiMaxCal(0), + fHadCorr(0), + fHCorrection(0), + fNTracks(0), + fNTracksCut(0), + fOpt(0), + fDZ(0), + fDebug(0), + fReaderHeader(0x0), + fMomentumArray(0x0), + fUnitArray(0x0), + fRefArray(0x0), + fTPCGrid(0x0), + fEMCalGrid(0x0), + fGeom(0x0), + fESD(0x0), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0), + fNphi(0), + fNeta(0), + fPhi2(0x0), + fEta2(0x0), + fPhi(0x0), + fEta(0x0), + fIndex(0x0), + fParams(0x0), + fGrid(0), + fPhiMin(0), + fPhiMax(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinInTPCAcc(0), + fPhiBinInTPCAcc(0), + fEtaBinInEMCalAcc(0), + fPhiBinInEMCalAcc(0), + fNbinPhi(0) { // constructor } -//____________________________________________________________________________ -void AliJetFillUnitArrayTracks::SetEMCALGeometry() +//_____________________________________________________________________________ +AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/) + : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"), + fNumUnits(0), + fEtaMinCal(0), + fEtaMaxCal(0), + fPhiMinCal(0), + fPhiMaxCal(0), + fHadCorr(0), + fHCorrection(0), + fNTracks(0), + fNTracksCut(0), + fOpt(0), + fDZ(0), + fDebug(0), + fReaderHeader(0x0), + fMomentumArray(0x0), + fUnitArray(0x0), + fRefArray(0x0), + fTPCGrid(0x0), + fEMCalGrid(0x0), + fGeom(0x0), + fESD(0x0), + fGrid0(0x0), + fGrid1(0x0), + fGrid2(0x0), + fGrid3(0x0), + fGrid4(0x0), + fNphi(0), + fNeta(0), + fPhi2(0x0), + fEta2(0x0), + fPhi(0x0), + fEta(0x0), + fIndex(0x0), + fParams(0x0), + fGrid(0), + fPhiMin(0), + fPhiMax(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinInTPCAcc(0), + fPhiBinInTPCAcc(0), + fEtaBinInEMCalAcc(0), + fPhiBinInEMCalAcc(0), + fNbinPhi(0) { - // Set EMCAL geometry information - fGeom = AliJetDummyGeo::GetInstance(); - if (fGeom == 0) - fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL"); - if(fDebug>1) printf("\n EMCAL Geometry setted ! \n"); + // constructor } //____________________________________________________________________________ void AliJetFillUnitArrayTracks::InitParameters() { - fHCorrection = 0; // For hadron correction + // fHCorrection = 0; // For hadron correction fHadCorr = 0; // For hadron correction fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL fDebug = fReaderHeader->GetDebug(); @@ -139,13 +166,6 @@ void AliJetFillUnitArrayTracks::InitParameters() fPhiMinCal = fGeom->GetArm1PhiMin(); fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur ! - if(fDebug>30){ - cout << "fEtaMinCal : " << fEtaMinCal << endl; - cout << "fEtaMaxCal : " << fEtaMaxCal << endl; - cout << "fPhiMinCal : " << fPhiMinCal << endl; - cout << "fPhiMaxCal : " << fPhiMaxCal << endl; - } - fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, fPhiMax,fEtaMin,fEtaMax); fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, @@ -173,349 +193,360 @@ AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks() } //_____________________________________________________________________________ -void AliJetFillUnitArrayTracks::Exec(Option_t* option) +void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) { // // Main method. - // Explain + // fDebug = fReaderHeader->GetDebug(); - if(fDebug>1) printf("In AliJetFillUnitArrayTracks !"); - if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells()); - // Set EMCal Geometry - SetEMCALGeometry(); + // Set parameters InitParameters(); - TClonesArray *lvArray = fMomentumArray; // Correct checked ! - Int_t nInT = lvArray->GetEntries(); // Correct checked ! - Float_t ptMin = fReaderHeader->GetPtCut(); // Correct checked ! - - // sflag -> not yet implemented !!!! - - if(fDebug>3) cout << "nInT : " << nInT << endl; - - if (nInT == 0) return; - - // local arrays for input - Float_t* enT = new Float_t[nInT]; - Float_t* ptT = new Float_t[nInT]; - Float_t* etaT = new Float_t[nInT]; - Float_t* phiT = new Float_t[nInT]; - Float_t* thetaT = new Float_t[nInT]; - - Int_t trackInEmcalAcc = 0; - Int_t trackInTpcAcc = 0; - Int_t trackInTpcAccOnly = 0; - - Int_t nElements = fTPCGrid->GetNEntries(); - Int_t nElements2 = fEMCalGrid->GetNEntries(); - fGrid = fTPCGrid->GetGridType(); + // get number of tracks in event (for the loop) + Int_t goodTrack = 0; + Int_t nt = 0; + Float_t pt, eta,phi; + TVector3 p3; + nt = fESD->GetNumberOfTracks(); + if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl; + + // temporary storage of signal and pt cut flag + Int_t* sflag = new Int_t[nt]; + Int_t* cflag = new Int_t[nt]; + + // get cuts set by user + Float_t ptMin = fReaderHeader->GetPtCut(); + Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); + Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); + fOpt = fReaderHeader->GetDetector(); + fDZ = fReaderHeader->GetDZ(); - if(fDebug>3){ - cout << "nElements : " << nElements << endl; - cout << "nElements2 : " << nElements2 << endl; - cout << "fNumUnits : " << fNumUnits << endl; - cout << "sum : " << fNumUnits+nElements << endl; - } + Int_t nTracksEmcal = 0; + Int_t nTracksEmcalDZ = 0; + Int_t nTracksTpc = 0; + Int_t nTracksTpcOnly = 0; + Int_t nTracksEmcalCut = 0; + Int_t nTracksEmcalDZCut = 0; + Int_t nTracksTpcCut = 0; + Int_t nTracksTpcOnlyCut = 0; - // Set energy exactly to zero - if(fGrid==0) - for(Int_t k=0; kAt(i); - enT[i] = lv->Energy(); - ptT[i] = lv->Pt(); - etaT[i] = lv->Eta(); - phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi()); - thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i])); - - if(fDebug>20){ - cout << "enT[" << i << "] : " << enT[i] << endl; - cout << "ptT[" << i << "] : " << ptT[i] << endl; - cout << "etaT[" << i << "] : " << etaT[i] << endl; - cout << "phiT[" << i << "] : " << phiT[i] << endl; - cout << "thetaT[" << i << "] : " << thetaT[i] << endl; - cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl; - cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl; - cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl; - cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl; - } - - if(fGrid==0) - { - // For the moment, only TPC filled from its grid in its total acceptance - if(fDebug>2) - cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl; - - trackInTpcAccOnly += 1; - - Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]); - - Float_t unitEnergy = 0.; - unitEnergy = fUnitArray[idTPC].GetUnitEnergy(); - - if(unitEnergy > 0. && fDebug >10){ - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "TPC id : " << idTPC << endl; - cout << "etaT[" << i << "] : " << etaT[i] << endl; - cout << "phiT[" << i << "] : " << phiT[i] << endl; - cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl; - cout << "fUnitArray[idTPC].GetUnitEnergy(): " << - fUnitArray[idTPC].GetUnitEnergy() << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - } + fGrid = fTPCGrid->GetGridType(); - // Fill energy in TPC acceptance - fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]); - if(fDebug > 10){ - cout << "ptT[" << i << "] : " << ptT[i] << endl; - cout << "unitEnergy in TPC acceptance after : " << - fUnitArray[idTPC].GetUnitEnergy() << endl; - } - - // Pt cut flag - if(fUnitArray[idTPC].GetUnitEnergy()GetTrack(it); + UInt_t status = track->GetStatus(); + + Double_t mom[3]; + track->GetPxPyPz(mom); + p3.SetXYZ(mom[0],mom[1],mom[2]); + pt = p3.Pt(); + Float_t mass = 0.; + mass = track->GetMass(); + + if (((status & AliESDtrack::kITSrefit) == 0) || + ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check + if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() + && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check + if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() + && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check + eta = p3.Eta(); + phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi()); + + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + + // sflag -> not yet implemented !!!! + + if(fGrid==0) + { + // Only TPC filled from its grid in its total acceptance + + Int_t idTPC = fTPCGrid->GetIndex(phi,eta); + Bool_t ok = kFALSE; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1); + uArray->SetUnitTrackID(it); + + Float_t unitEnergy = 0.; + unitEnergy = uArray->GetUnitEnergy(); + if(unitEnergy==0.){ + nTracksTpcOnly++; + ok = kTRUE; + } + // Fill energy in TPC acceptance + uArray->SetUnitEnergy(unitEnergy + pt); + uArray->SetUnitPxPyPz(mom); + uArray->SetUnitMass(mass); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksTpcOnlyCut++; } - if(fGrid==1) - { - // Fill track information in EMCAL acceptance - if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) && - (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// && - // GetCutFlag(i) == 1) - // ptT[i] > ptMin) - { - trackInEmcalAcc += 1; - - if(fDebug>20){ - cout << "before : " << endl; - cout << "etaT[i] : " << etaT[i] << endl; - cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl; - } - - // This function should be modified soon -// Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete - Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified -// Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]); // Using an EMCal grid -> calculated -// Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]); // Using an EMCal grid -> tabulated (faster) - - Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy(); - - if(fDebug>20) { - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "after : " << endl; - cout << "towerID : " << towerID << endl; - cout << "etaT[i] : " << etaT[i] << endl; - cout << "phiT[i](rad) : " << phiT[i] << endl; - cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl; - cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl; - cout << "fHadCorr : " << fHadCorr << endl; - cout << "fHCorrection : " << fHCorrection << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - } - - //OLD WAY: //Do Hadron Correction - if (fHCorrection != 0) - { - Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]); - unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]); - if(fDebug>20){ - cout << "Inside loop for hadron correction ==========================" << endl; - cout << "enT[" << i << "] : " << enT[i] << endl; - cout << "ptT[" << i << "] : " << ptT[i] << endl; - cout << "etaT[" << i << "] : " << etaT[i] << endl; - cout << "Energy correction : " << hCEnergy << endl; - cout << "unitEnergy : " << unitEnergy << endl; - cout << "fUnitArray[towerID].GetUnitEnergy() : " << - fUnitArray[towerID].GetUnitEnergy() << endl; - } - } //end Hadron Correction loop - - fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]); - - // Put a pt cut flag - if(fUnitArray[towerID].GetUnitEnergy()10){ - cout << "After pT filled ===============================================" << endl; - cout << "PtCut : " << ptMin << endl; - cout << "ptT[" << i << "] : " << ptT[i] << endl; - cout << "unitEnergy : " << unitEnergy << endl; - cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl; - cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl; - cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl; - cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl; - } - } // end loop on EMCal acceptance cut + tracks quality - else{ - // Outside EMCal acceptance - if(fDebug>2) - cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl; - - trackInTpcAcc += 1; - - // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]); - Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]); - - Float_t unitEnergy2 = 0.; - unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1 - - if(unitEnergy2 > 0. && fDebug >10){ - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "TPC id : " << idTPC << endl; - cout << "Total id : " << fNumUnits-1+idTPC << endl; - cout << "etaT[" << i << "] : " << etaT[i] << endl; - cout << "phiT[" << i << "] : " << phiT[i] << endl; - cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl; - cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " << - fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl; - } + // Detector flag + if(unitEnergy>0) { + uArray->SetUnitDetectorFlag(kAll); + } + if(uArray->GetUnitEnergy()>0){ + fRefArray->Add(uArray); + } - // Fill energy outside emcal acceptance - fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]); + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + goodTrack++; - // Pt cut flag - if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()3){ - printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc); - printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc); - } - - if(fGrid==0) fNIn = trackInTpcAccOnly; - if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc; - - fOpt = fReaderHeader->GetDetector(); - - if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt); - - if(fOpt==1 || option=="tpc") // if only TPC - { //Set all unit flags, Eta, Phi - - if(fGrid==0) - { - if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n"); - for(Int_t i=0; i10) Info("FillUnitArray","Setting all units outside jets \n"); - fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2... - fUnitArray[i].SetUnitEntries(nElements); - fUnitArray[i].SetUnitID(i); - fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); - fUnitArray[i].SetUnitEta(eta); - fUnitArray[i].SetUnitPhi(phi); - fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta()); - fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi()); + } + + if(fGrid==1) + { + Int_t nElements = fTPCGrid->GetNEntries(); + + // Fill track information in EMCAL acceptance + if((eta >= fEtaMin && eta <= fEtaMax) && + (phi >= fPhiMin && phi <= fPhiMax))// && + { + + // Include dead-zones + if(fDZ) + { + Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.; + Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.; + fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0); + fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1); + fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2); + fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3); + fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4); + Int_t n0 = fGrid0->GetNEntries(); + Int_t n1 = fGrid1->GetNEntries(); + Int_t n2 = fGrid2->GetNEntries(); + Int_t n3 = fGrid3->GetNEntries(); + + if(phi >= phimin0 && phi <= phimax0){ + Int_t id0 = fGrid0->GetIndex(phi,eta)-1; + AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements); + uArray0->SetUnitTrackID(it); + Float_t uEnergy0 = uArray0->GetUnitEnergy(); + Bool_t ok0 = kFALSE; + if(uEnergy0==0.){ + nTracksEmcalDZ++; + ok0 = kTRUE; + } + uArray0->SetUnitEnergy(uEnergy0+pt); + if(uArray0->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray0->SetUnitCutFlag(kPtHigher); + if(ok0) nTracksEmcalDZCut++; + } + if(uArray0->GetUnitEnergy()>0) + fRefArray->Add(uArray0); } - } - - if(fGrid==1) - { - if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n"); - for(Int_t i=0; i10) Info("FillUnitArray","Setting all units outside jets \n"); - //Set all units to be outside a jet initially - fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2... - fUnitArray[i].SetUnitEntries(fNumUnits+nElements); - fUnitArray[i].SetUnitID(i); - - if(fUnitArray[i].GetUnitID()EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63 - fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry - phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); - // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid - // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid - fUnitArray[i].SetUnitEta(eta); - //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0); - fUnitArray[i].SetUnitPhi(phi); - fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta()); - fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi()); - } + if(phi >= phimin1 && phi <= phimax1){ + Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0; + AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements); + uArray1->SetUnitTrackID(it); + Float_t uEnergy1 = uArray1->GetUnitEnergy(); + Bool_t ok1 = kFALSE; + if(uEnergy1==0.){ + nTracksEmcalDZ++; + ok1 = kTRUE; + } + uArray1->SetUnitEnergy(uEnergy1+pt); + if(uArray1->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { - fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta); - fUnitArray[i].SetUnitEta(eta); - fUnitArray[i].SetUnitPhi(phi); - fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta()); - fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi()); - - if(fDebug>10) - { - if(fUnitArray[i].GetUnitEnergy()!=0.){ - cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " << - fUnitArray[i].GetUnitID()+1-fNumUnits << endl; - cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " << - fUnitArray[i].GetUnitEnergy() << endl; - cout << "(fUnitArray[" << i << "].GetUnitEta() : " << - fUnitArray[i].GetUnitEta() << endl; - cout << "(fUnitArray[" << i << "].GetUnitPhi() : " << - fUnitArray[i].GetUnitPhi() << endl; - } - } + uArray1->SetUnitCutFlag(kPtHigher); + if(ok1) nTracksEmcalDZCut++; } - fUnitArray[i].SetUnitClusterID(0); - }//end loop over all units in array (same as all towers in EMCAL) - } - - if(fDebug>20) - { - for(Int_t i=0; iGetUnitEnergy()>0) + fRefArray->Add(uArray1); + } + if(phi >= phimin2 && phi <= phimax2){ + Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1; + AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements); + uArray2->SetUnitTrackID(it); + Float_t uEnergy2 = uArray2->GetUnitEnergy(); + Bool_t ok2 = kFALSE; + if(uEnergy2==0.){ + nTracksEmcalDZ++; + ok2 = kTRUE; + } + uArray2->SetUnitEnergy(uEnergy2+pt); + if(uArray2->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray2->SetUnitCutFlag(kPtHigher); + if(ok2) nTracksEmcalDZCut++; + } + if(uArray2->GetUnitEnergy()>0) + fRefArray->Add(uArray2); + } + if(phi >= phimin3 && phi <= phimax3){ + Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2; + AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements); + uArray3->SetUnitTrackID(it); + Float_t uEnergy3 = uArray3->GetUnitEnergy(); + Bool_t ok3 = kFALSE; + if(uEnergy3==0.){ + nTracksEmcalDZ++; + ok3 = kTRUE; + } + uArray3->SetUnitEnergy(uEnergy3+pt); + if(uArray3->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray3->SetUnitCutFlag(kPtHigher); + if(ok3) nTracksEmcalDZCut++; + } + if(uArray3->GetUnitEnergy()>0) + fRefArray->Add(uArray3); + } + if(phi >= phimin4 && phi <= phimax4){ + Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3; + AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements); + uArray4->SetUnitTrackID(it); + Float_t uEnergy4 = uArray4->GetUnitEnergy(); + Bool_t ok4 = kFALSE; + if(uEnergy4==0.){ + nTracksEmcalDZ++; + ok4 = kTRUE; + } + uArray4->SetUnitEnergy(uEnergy4+pt); + if(uArray4->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray4->SetUnitCutFlag(kPtHigher); + if(ok4) nTracksEmcalDZCut++; } + if(uArray4->GetUnitEnergy()>0) + fRefArray->Add(uArray4); } + } // end fDZ + + Int_t towerID = 0; + fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID); + + if(towerID==-1) continue; + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID); + uArray->SetUnitTrackID(it); + Float_t unitEnergy = uArray->GetUnitEnergy(); + Bool_t ok = kFALSE; + if(unitEnergy==0.){ + nTracksEmcal++; + ok=kTRUE; + } + + // Do Hadron Correction + // Parametrization to be added + if (fHCorrection != 0) + { + // Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]); + Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta); + unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); + + } //end Hadron Correction loop + + uArray->SetUnitEnergy(unitEnergy + pt); + + // Put a pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok) nTracksEmcalCut++; } + // Detector flag + if(unitEnergy > 0) + uArray->SetUnitDetectorFlag(kAll); + + if(uArray->GetUnitEnergy()>0) + fRefArray->Add(uArray); + + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + goodTrack++; + + } // end loop on EMCal acceptance cut + tracks quality + else{ + // Outside EMCal acceptance + + // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]); + Int_t idTPC = fTPCGrid->GetIndex(phi,eta); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC); + uArray->SetUnitTrackID(it); + + Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1 + Bool_t ok2 = kFALSE; + if(unitEnergy2==0.){ + nTracksTpc++; + ok2=kTRUE; + } + // Fill energy outside emcal acceptance + uArray->SetUnitEnergy(unitEnergy2 + pt); + + // Pt cut flag + if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } + else { + uArray->SetUnitCutFlag(kPtHigher); + if(ok2) nTracksTpcCut++; + } + // Detector flag + if(unitEnergy2 > 0) + uArray->SetUnitDetectorFlag(kTpc); + if(uArray->GetUnitEnergy()>0) + fRefArray->Add(uArray); + + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + goodTrack++; + + } + } // end fGrid==1 + } // end loop on entries (tpc tracks) - } // end if(fOpt==1 || option=="tpc") +// // set the signal flags +// fSignalFlag.Set(goodTrack,sflag); +// fCutFlag.Set(goodTrack,cflag); - delete[] enT; - delete[] ptT; - delete[] etaT; - delete[] phiT; - delete[] thetaT; +// delete sflag; +// delete cflag; + + // } // end loop on entries (tpc tracks) + if(fGrid==0) { + fNTracks = nTracksTpcOnly; + fNTracksCut = nTracksTpcOnlyCut; + if(fDebug>10){ + cout << "fNTracks : " << fNTracks << endl; + cout << "fNTracksCut : " << fNTracksCut << endl; + } + } + if(fGrid==1) { + fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc; + fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut; + if(fDebug>10){ + cout << "fNTracks : " << fNTracks << endl; + cout << "fNTracksCut : " << fNTracksCut << endl; + } + } + } //__________________________________________________________ diff --git a/JETAN/AliJetFillUnitArrayTracks.h b/JETAN/AliJetFillUnitArrayTracks.h index 08e6356b9a0..42fc6ead4e1 100644 --- a/JETAN/AliJetFillUnitArrayTracks.h +++ b/JETAN/AliJetFillUnitArrayTracks.h @@ -17,41 +17,56 @@ #include #include -class AliEMCALGeometry; -class AliJetDummyGeo; class AliJetHadronCorrection; class AliJetReader; class AliJetESDReader; class TClonesArray; -class AliJetUnitArray; +class TRefArray; +//class AliJetUnitArray; //class AliJetReaderHeader; -class AliJetReader; class AliJetGrid; +class AliJetDummyGeo; +class AliESD; +class AliESDEvent; class AliJetFillUnitArrayTracks : public TTask { public: AliJetFillUnitArrayTracks(); + AliJetFillUnitArrayTracks(AliESD *fESD); + AliJetFillUnitArrayTracks(AliESDEvent *fESD); virtual ~AliJetFillUnitArrayTracks(); // Setter void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;} + void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;} void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;} - void SetUnitArray(AliJetUnitArray *unitArray) {fUnitArray = unitArray;} + void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;} + void SetRefArray(TRefArray *refArray) {fRefArray = refArray;} void SetHadCorrection(Int_t flag = 1) {fHCorrection = flag;} void SetHadCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;} void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;} void SetGrid(Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax); + // void SetESD(AliESD *esd) {fESD = esd;} + void SetESD(AliESDEvent *esd) {fESD = esd;} + void SetGrid0(AliJetGrid *grid0){fGrid0 = grid0;} + void SetGrid1(AliJetGrid *grid1){fGrid1 = grid1;} + void SetGrid2(AliJetGrid *grid2){fGrid2 = grid2;} + void SetGrid3(AliJetGrid *grid3){fGrid3 = grid3;} + void SetGrid4(AliJetGrid *grid4){fGrid4 = grid4;} // Getter - AliJetUnitArray* GetUnitArray() {return fUnitArray;} - // Int_t GetIndexFromEtaPhi(Double_t eta,Double_t phi) const; - void GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi); - Int_t GetNeta() {return fNeta;} - Int_t GetNphi() {return fNphi;} - - void Exec(Option_t*); + TClonesArray* GetUnitArray() {return fUnitArray;} + TRefArray* GetRefArray() {return fRefArray;} + // Int_t GetIndexFromEtaPhi(Double_t eta,Double_t phi) const; + void GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi); + Int_t GetNeta() {return fNeta;} + Int_t GetNphi() {return fNphi;} + Int_t GetHadCorrection() {return fHCorrection;} + Int_t GetMult() {return fNTracks;} + Int_t GetMultCut() {return fNTracksCut;} + void Exec(Option_t*); protected: Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) @@ -60,19 +75,27 @@ class AliJetFillUnitArrayTracks : public TTask Float_t fPhiMinCal; // Define EMCal acceptance in Phi Float_t fPhiMaxCal; // Define EMCal acceptance in Phi AliJetHadronCorrectionv1 *fHadCorr; // Pointer to Hadron Correction Object - Int_t fHCorrection; // Hadron correction flag - Int_t fNIn; // Number of Array filled in UnitArray + Int_t fHCorrection; // Hadron correction flag + 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 - AliJetUnitArray *fUnitArray; // UnitArray + TClonesArray *fUnitArray; // UnitArray + TRefArray *fRefArray; // UnitArray AliJetGrid *fTPCGrid; // Define filled grid AliJetGrid *fEMCalGrid; // Define filled grid - - // geometry info - AliJetDummyGeo *fGeom; //! + AliJetDummyGeo *fGeom; // Define EMCal geometry + // AliESD *fESD; // ESD + AliESDEvent *fESD; // ESD + AliJetGrid *fGrid0; + AliJetGrid *fGrid1; + AliJetGrid *fGrid2; + AliJetGrid *fGrid3; + AliJetGrid *fGrid4; Int_t fNphi; // number of points in the grid: phi Int_t fNeta; // " eta @@ -95,11 +118,9 @@ class AliJetFillUnitArrayTracks : public TTask Int_t fNbinPhi; private: - void SetEMCALGeometry(); + // void SetEMCALGeometry(); void InitParameters(); - // Int_t fEvent; - ClassDef(AliJetFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information }; diff --git a/JETAN/AliJetFinder.cxx b/JETAN/AliJetFinder.cxx index 58d19fb697e..02416cc4b6c 100644 --- a/JETAN/AliJetFinder.cxx +++ b/JETAN/AliJetFinder.cxx @@ -253,7 +253,5 @@ void AliJetFinder::AddJet(AliAODJet p) void AliJetFinder::ConnectAOD(AliAODEvent* aod) { // Connect to the AOD - printf("Connect AOD \n"); fAODjets = aod->GetJets(); - printf("Connect AOD %p \n", fAODjets); } diff --git a/JETAN/AliJetHadronCorrectionv1.cxx b/JETAN/AliJetHadronCorrectionv1.cxx index 61503d00bc0..5ac309a2af3 100644 --- a/JETAN/AliJetHadronCorrectionv1.cxx +++ b/JETAN/AliJetHadronCorrectionv1.cxx @@ -26,7 +26,6 @@ // --- AliRoot header files --- //#include "../EMCAL/AliJetGeometry.h" #include "AliJetDummyGeo.h" -//#include "AliEMCALGeometry.h" #include "AliJetHadronCorrectionv1.h" ClassImp(AliJetHadronCorrectionv1) diff --git a/JETAN/AliJetHadronCorrectionv1.h b/JETAN/AliJetHadronCorrectionv1.h index b8da7fc7920..36dc6631be9 100644 --- a/JETAN/AliJetHadronCorrectionv1.h +++ b/JETAN/AliJetHadronCorrectionv1.h @@ -11,7 +11,6 @@ #define HCPARAMETERS 6 #define HCPARAMETERSETS 2 -class AliEMCALGeometry; class AliJetDummyGeo; diff --git a/JETAN/AliJetMCReader.cxx b/JETAN/AliJetMCReader.cxx index 60078620d90..7ebb8ac7ff0 100644 --- a/JETAN/AliJetMCReader.cxx +++ b/JETAN/AliJetMCReader.cxx @@ -25,6 +25,7 @@ #include #include #include +#include // From AliRoot ... #include "AliJetMCReader.h" #include "AliJetMCReaderHeader.h" diff --git a/JETAN/AliJetReader.cxx b/JETAN/AliJetReader.cxx index 2f6dd92fa48..db62f0f3c38 100755 --- a/JETAN/AliJetReader.cxx +++ b/JETAN/AliJetReader.cxx @@ -22,9 +22,15 @@ // root #include +#include +#include "TTask.h" //AliRoot #include "AliJetReader.h" #include "AliJetReaderHeader.h" +#include "AliESDEvent.h" +#include "AliHeader.h" +#include "AliJetFillUnitArrayTracks.h" +#include "AliJetFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" #include "AliJetHadronCorrectionv1.h" @@ -33,15 +39,23 @@ ClassImp(AliJetReader) //////////////////////////////////////////////////////////////////////// AliJetReader::AliJetReader(): + // Constructor + fChain(0), fMomentumArray(new TClonesArray("TLorentzVector",2000)), fArrayMC(0), fFillUnitArray(new TTask("fillUnitArray","Fill unit array jet finder")), + fESD(0), fReaderHeader(0), fSignalFlag(0), fCutFlag(0), - fUnitArray(new AliJetUnitArray[60000]), - fUnitArrayNoCuts(new AliJetUnitArray[60000]), - fArrayInitialised(0) + fUnitArray(new TClonesArray("AliJetUnitArray",60000)), + fRefArray(new TRefArray()), + fUnitArrayNoCuts(new TClonesArray("AliJetUnitArray",60000)), + fArrayInitialised(0), + fFillUAFromTracks(new AliJetFillUnitArrayTracks()), + fFillUAFromEMCalDigits(new AliJetFillUnitArrayEMCalDigits()) + // fHadronCorrector(0), + // fHCorrection(0) { // Default constructor fSignalFlag = TArrayI(); diff --git a/JETAN/AliJetReader.h b/JETAN/AliJetReader.h index 3b63c3282cd..a104c345c7c 100755 --- a/JETAN/AliJetReader.h +++ b/JETAN/AliJetReader.h @@ -12,18 +12,20 @@ #include #include #include -#ifndef ROOT_TTask -#include "TTask.h" -#endif class TTree; class TTask; +class TChain; class TClonesArray; class TRefArray; class AliJetReaderHeader; +class AliESDEvent; +class AliHeader; class AliJetUnitArray; class AliJetHadronCorrectionv1; class AliJet; +class AliJetFillUnitArrayTracks; +class AliJetFillUnitArrayEMCalDigits; class AliJetReader : public TObject { @@ -32,15 +34,18 @@ class AliJetReader : public TObject virtual ~AliJetReader(); // Getters - virtual TClonesArray* GetMomentumArray() const {return fMomentumArray;} - virtual TRefArray* GetReferences() const {return 0;} - virtual AliJetUnitArray* GetUnitArray() const {return fUnitArray;} - virtual AliJetUnitArray* GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;} - - virtual AliJetReaderHeader* GetReaderHeader() const {return fReaderHeader;} - virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];} - virtual Int_t GetCutFlag(Int_t i) const {return fCutFlag[i];} - virtual Int_t GetArrayInitialised() const {return fArrayInitialised;} + virtual TClonesArray* GetMomentumArray() const {return fMomentumArray;} + virtual TRefArray* GetReferences() const {return 0;} + virtual TClonesArray *GetUnitArray() const {return fUnitArray;} + virtual TRefArray *GetRefArray() const {return fRefArray;} + virtual TClonesArray *GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;} + + virtual AliJetReaderHeader* GetReaderHeader() const {return fReaderHeader;} + virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];} + virtual Int_t GetCutFlag(Int_t i) const {return fCutFlag[i];} + virtual Int_t GetArrayInitialised() const {return fArrayInitialised;} + virtual Int_t GetNumCandidate() const {return fNumCandidate;} + virtual Int_t GetNumCandidateCut() const {return fNumCandidateCut;} // Setters virtual Bool_t FillMomentumArray(Int_t) {return kTRUE;} @@ -49,9 +54,20 @@ class AliJetReader : public TObject virtual void FillUnitArrayFromEMCALDigits(Int_t) {} // temporarily not used virtual void FillUnitArrayFromEMCALClusters(Int_t) {} // temporarily not used virtual void InitUnitArray() {} + virtual void InitParameters() {} + virtual void CreateTasks() {} + // virtual void ExecTasks(Int_t) {} + virtual Bool_t ExecTasks(Int_t) {return kTRUE;} + /* // Correction of hadronic energy + virtual void SetHadronCorrector(AliEMCALHadronCorrectionv1* corr) {fHadronCorrector = corr;} + virtual void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} */ virtual void SetReaderHeader(AliJetReaderHeader* header) {fReaderHeader = header;} - + virtual void SetESD(AliESDEvent* esd) { fESD = esd;} + // virtual Int_t SetNumCandidate(Int_t cand) {fNumCandidate = cand;} + // virtual Int_t SetNumCandidateCut(Int_t candcut) {fNumCandidateCut = candcut;} + + // Others virtual void OpenInputFiles() {} virtual void SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* /*mc*/) {;} @@ -63,16 +79,23 @@ class AliJetReader : public TObject protected: AliJetReader(const AliJetReader& rJetReader); AliJetReader& operator = (const AliJetReader& rhsr); + TChain *fChain; // chain for reconstructed tracks TClonesArray *fMomentumArray; // array of particle momenta TClonesArray *fArrayMC; //! array of mc particles TTask *fFillUnitArray; //! task list for filling the UnitArray + AliESDEvent *fESD; // pointer to esd AliJetReaderHeader *fReaderHeader; // pointer to header TArrayI fSignalFlag; // to flag if a particle comes from pythia or // from the underlying event TArrayI fCutFlag; // to flag if a particle passed the pt cut or not - AliJetUnitArray *fUnitArray; //! array of digit position and energy - AliJetUnitArray *fUnitArrayNoCuts; //! array of digit position and energy + TClonesArray *fUnitArray; // array of digit position and energy + TRefArray *fRefArray; // 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; + AliJetFillUnitArrayEMCalDigits *fFillUAFromEMCalDigits; + Int_t fNumCandidate; + Int_t fNumCandidateCut; ClassDef(AliJetReader,1) }; diff --git a/JETAN/AliJetReaderHeader.cxx b/JETAN/AliJetReaderHeader.cxx index f172d845d83..7642a8dee19 100755 --- a/JETAN/AliJetReaderHeader.cxx +++ b/JETAN/AliJetReaderHeader.cxx @@ -19,6 +19,8 @@ // Author: jgcn@mda.cinvestav.mx //------------------------------------------------------------------------- +#include + #include "AliJetReaderHeader.h" ClassImp(AliJetReaderHeader) @@ -30,9 +32,12 @@ AliJetReaderHeader::AliJetReaderHeader(): fFirst(0), fLast(-1), fOption(0), + fDZ(0), fSignalPerBg(0), fFiducialEtaMin(-0.9), fFiducialEtaMax(0.9), + fFiducialPhiMin(0.), + fFiducialPhiMax(2*TMath::Pi()), fPtCut(2.0), fComment("No comment"), fDir(""), @@ -49,9 +54,12 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name): fFirst(0), fLast(-1), fOption(0), + fDebug(0), fSignalPerBg(0), fFiducialEtaMin(-0.9), fFiducialEtaMax(0.9), + fFiducialPhiMin(0.), + fFiducialPhiMax(2*TMath::Pi()), fPtCut(2.0), fComment("No comment"), fDir(""), diff --git a/JETAN/AliJetReaderHeader.h b/JETAN/AliJetReaderHeader.h index 2fbbf8426e2..cbf8208d1ed 100755 --- a/JETAN/AliJetReaderHeader.h +++ b/JETAN/AliJetReaderHeader.h @@ -27,12 +27,15 @@ class AliJetReaderHeader : public TNamed virtual const char* GetBgDirectory(){return fBgDir.Data();} virtual const char* GetPattern() {return fPattern.Data();} virtual Float_t GetFiducialEtaMin() const {return fFiducialEtaMin;} - virtual Float_t GetFiducialEtaMax() const {return fFiducialEtaMax;} + virtual Float_t GetFiducialEtaMax() const {return fFiducialEtaMax;} + virtual Float_t GetFiducialPhiMin() const {return fFiducialPhiMin;} + virtual Float_t GetFiducialPhiMax() const {return fFiducialPhiMax;} virtual Float_t GetPtCut() const {return fPtCut;} Int_t GetNEvents() const {return fLast-fFirst;} Int_t GetFirstEvent() const {return fFirst;} Int_t GetLastEvent() const {return fLast;} Int_t GetDetector() const {return fOption;} + Bool_t GetDZ() const {return fDZ;} Int_t GetDebug() const {return fDebug;} Int_t GetSignalPerBg() const {return fSignalPerBg;} @@ -46,18 +49,24 @@ class AliJetReaderHeader : public TNamed virtual void SetLastEvent(Int_t i=-1) {fLast=i;} virtual void SetFiducialEta(Float_t etamin, Float_t etamax) { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;} + virtual void SetFiducialPhi(Float_t phimin, Float_t phimax) + { fFiducialPhiMin = phimin; fFiducialPhiMax = phimax;} virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;} - + virtual void SetDZ(Bool_t deadzone = 0) {fDZ = deadzone;} virtual void SetDetector(Int_t option = 0) {fOption = option;} virtual void SetDebug(Int_t debug = 0) {fDebug = debug;} + protected: Int_t fFirst; // First and last events analyzed Int_t fLast; // in current set of files Int_t fOption; // detector used for jet reconstruction Int_t fDebug; // debug option + Bool_t fDZ; // include dead zones or not Int_t fSignalPerBg; Float_t fFiducialEtaMin; // Fiducial minimum eta Float_t fFiducialEtaMax; // Fiducial maximum eta + Float_t fFiducialPhiMin; // Fiducial minimum phi + Float_t fFiducialPhiMax; // Fiducial maximum phi Float_t fPtCut; // pt cut TString fComment; // a comment TString fDir; // directory with input files for signal diff --git a/JETAN/AliJetUnitArray.cxx b/JETAN/AliJetUnitArray.cxx index 353a7e6bbcb..1c443d4b43b 100644 --- a/JETAN/AliJetUnitArray.cxx +++ b/JETAN/AliJetUnitArray.cxx @@ -13,40 +13,61 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ - -/* $Id$ */ - -//_________________________________________________________________________ +//------------------------------------------------------------------ // Unit used by UA1 algorithm -// -- -//*-- Author: Sarah Blyth (LBL/UCT) -// -- -// Revised Version for JETAN - 30/03/2006 -// -- Magali Estienne (IReS) +// Authors: Sarah Blyth (LBL/UCT) +// Magali Estienne (IReS) (new version for JETAN) +//------------------------------------------------------------------ #include "AliJetUnitArray.h" ClassImp(AliJetUnitArray) - AliJetUnitArray::AliJetUnitArray(): - fUnitEnergy(0.), - fUnitEta(0.), - fUnitPhi(0.), - fUnitDeta(0.), - fUnitDphi(0.), - fUnitID(0), - fUnitNum(0), - fUnitClusterID(0), - fUnitFlag(kOutJet), - fUnitCutFlag(kPtSmaller), - fUnitSignalFlag(kGood), - fUnitDetectorFlag(kAll) + fUnitEnergy(0.0), + fUnitEta(0.0), + fUnitPhi(0.0), + fUnitDeta(0.), + fUnitDphi(0.), + fUnitID(0), + fUnitTrackID(0), + fUnitNum(0), + fUnitClusterID(0), + fUnitFlag(kOutJet), + fUnitCutFlag(kPtSmaller), + fUnitSignalFlag(kBad), + fUnitDetectorFlag(kTpc), + fUnitPx(0.), + fUnitPy(0.), + fUnitPz(0.), + fUnitMass(0.) { // Default constructor } +AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId): + fUnitEnergy(en), + fUnitEta(eta), + fUnitPhi(phi), + fUnitDeta(Deta), + fUnitDphi(Dphi), + fUnitID(absId), + fUnitTrackID(esdId), + fUnitNum(0), + fUnitClusterID(clusId), + fUnitFlag(inout), + fUnitCutFlag(cut), + fUnitSignalFlag(kBad), + fUnitDetectorFlag(det), + fUnitPx(px), + fUnitPy(py), + fUnitPz(pz), + fUnitMass(mass) +{ + // Constructor 2 +} + AliJetUnitArray::~AliJetUnitArray() { // Destructor diff --git a/JETAN/AliJetUnitArray.h b/JETAN/AliJetUnitArray.h index 79d5fa512b1..20e810c2f45 100644 --- a/JETAN/AliJetUnitArray.h +++ b/JETAN/AliJetUnitArray.h @@ -1,19 +1,15 @@ -#ifndef ALIUA1JETFINDERUNIT_H -#define ALIUA1JETFINDERUNIT_H +#ifndef ALIJETUNITARRAY_H +#define ALIJETUNITARRAY_H /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * See cxx source for full Copyright notice */ -/* $Id$ */ - -//====================================================================== -// Revised Version for JETAN - 08/11/2006 -// Functions added -// Author: magali.estienne@ires.in2p3.fr -//====================================================================== +//------------------------------------------------------------------ // Unit used by UA1 algorithm -// -//*-- Author: Sarah Blyth (LBL/UCT) +// Authors: Sarah Blyth (LBL/UCT) +// Magali Estienne (IReS) (new version for JETAN) +//------------------------------------------------------------------ + #include #include "AliJetFinderTypes.h" @@ -22,6 +18,7 @@ class AliJetUnitArray : public TObject { public: AliJetUnitArray(); + AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId); ~AliJetUnitArray(); // Setter @@ -31,6 +28,7 @@ class AliJetUnitArray : public TObject void SetUnitDeta(Float_t deta) {fUnitDeta = deta;} void SetUnitDphi(Float_t dphi) {fUnitDphi = dphi;} void SetUnitID(Int_t id) {fUnitID = id;} + void SetUnitTrackID(Int_t esdid) {fUnitTrackID = esdid;} void SetUnitEntries(Int_t num) {fUnitNum = num;} void SetUnitClusterID(Int_t id) {fUnitClusterID = id;} void SetUnitFlag(AliJetFinderUnitFlagType_t flag) @@ -49,6 +47,8 @@ class AliJetUnitArray : public TObject { fUnitDetectorFlag = detectorflag; } + void SetUnitPxPyPz(Double_t *pxyz) {fUnitPx = pxyz[0]; fUnitPy = pxyz[1]; fUnitPz = pxyz[2];} + void SetUnitMass(Float_t mass) {fUnitMass = mass;} // Getter Float_t GetUnitEnergy() const {return fUnitEnergy;} @@ -57,8 +57,12 @@ class AliJetUnitArray : public TObject 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(Double_t* p) const {p[0]=fUnitPx; p[1]=fUnitPy; p[2]=fUnitPz; return kTRUE;} + AliJetFinderUnitFlagType_t GetUnitFlag() const { return fUnitFlag; @@ -76,25 +80,31 @@ class AliJetUnitArray : public TObject return fUnitDetectorFlag; } + Bool_t operator> ( AliJetUnitArray unit1) const; Bool_t operator< ( AliJetUnitArray unit1) const; Bool_t operator== ( AliJetUnitArray unit1) const; protected: - Float_t fUnitEnergy; // Energy of the unit - Float_t fUnitEta; // Eta of the unit - Float_t fUnitPhi; // Phi of the unit - Float_t fUnitDeta; // Delta Eta of the unit - Float_t fUnitDphi; // Delta Phi of the unit - Int_t fUnitID; // ID of the unit - Int_t fUnitNum; // number of units - Int_t fUnitClusterID; // ID of the unit - AliJetFinderUnitFlagType_t fUnitFlag; // Flag of the unit - AliJetFinderUnitCutFlagType_t fUnitCutFlag; // Flag of the unit - AliJetFinderUnitSignalFlagType_t fUnitSignalFlag; // Flag of the unit - AliJetFinderUnitDetectorFlagType_t fUnitDetectorFlag; // Detector flag of the unit + Float_t fUnitEnergy; // Energy of the unit + Float_t fUnitEta; // Eta of the unit + Float_t fUnitPhi; // Phi of the unit + Float_t fUnitDeta; // Delta Eta of the unit + Float_t fUnitDphi; // Delta Phi of the unit + Int_t fUnitID; // ID of the unit + Int_t fUnitTrackID; // ID of the unit + Int_t fUnitNum; // number of units + Int_t fUnitClusterID; // ID of the unit + AliJetFinderUnitFlagType_t fUnitFlag; // Flag of the unit + AliJetFinderUnitCutFlagType_t fUnitCutFlag; // Flag of the unit + AliJetFinderUnitSignalFlagType_t fUnitSignalFlag; // Flag of the unit + AliJetFinderUnitDetectorFlagType_t fUnitDetectorFlag; // Detector flag of the unit + Float_t fUnitPx; // Px of charged track + Float_t fUnitPy; // Py of charged track + Float_t fUnitPz; // Pz of charged track + Float_t fUnitMass; // Mass of charged particle - ClassDef(AliJetUnitArray,3) + ClassDef(AliJetUnitArray,1) }; #endif diff --git a/JETAN/AliUA1JetFinderV2.cxx b/JETAN/AliUA1JetFinderV2.cxx new file mode 100644 index 00000000000..81da6b2d735 --- /dev/null +++ b/JETAN/AliUA1JetFinderV2.cxx @@ -0,0 +1,878 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* $Id$ */ + +//--------------------------------------------------------------------- +// UA1 Cone Algorithm Jet finder +// manages the search for jets +// Author: Rafael.Diaz.Valdes@cern.ch +// (version in c++) +// Modified to include neutral particles (magali.estienne@ires.in2p3.fr) +//--------------------------------------------------------------------- + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliUA1JetFinderV2.h" +#include "AliUA1JetHeaderV1.h" +#include "AliJetUnitArray.h" +#include "AliJetReaderHeader.h" +#include "AliJetReader.h" +#include "AliJet.h" +#include "AliAODJet.h" + + +ClassImp(AliUA1JetFinderV2) + +//////////////////////////////////////////////////////////////////////// + + AliUA1JetFinderV2::AliUA1JetFinderV2(): + fDebug(0), + fOpt(0) +{ + // Constructor + fHeader = 0x0; + fLego = 0x0; +} + +//////////////////////////////////////////////////////////////////////// + +AliUA1JetFinderV2::~AliUA1JetFinderV2() + +{ + // destructor +} + +//////////////////////////////////////////////////////////////////////// + + +void AliUA1JetFinderV2::FindJets() + +{ + //1) Fill cell map array + //2) calculate total energy and fluctuation level + //3) Run algorithm + // 3.1) look centroides in cell map + // 3.2) calculate total energy in cones + // 3.3) flag as a possible jet + // 3.4) reorder cones by energy + //4) subtract backg in accepted jets + //5) fill AliJet list + + // transform input to pt,eta,phi plus lego + + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + TClonesArray* fUnit = fReader->GetUnitArray(); + Int_t nCandidate = fReader->GetNumCandidate(); + Int_t nIn = fUnit->GetEntries(); + + if (nIn == 0) return; + + // local arrays for input No Cuts + // Both pt < ptMin and pt > ptMin + Float_t* enT = new Float_t[nCandidate]; + Float_t* ptT = new Float_t[nCandidate]; + Float_t* etaT = new Float_t[nCandidate]; + Float_t* phiT = new Float_t[nCandidate]; + Float_t* detaT = new Float_t[nCandidate]; + Float_t* dphiT = new Float_t[nCandidate]; + Float_t* cFlagT = new Float_t[nCandidate]; + Float_t* cClusterT = new Float_t[nCandidate]; + Float_t* idT = new Float_t[nCandidate]; + Int_t loop1 = 0; + Int_t* injet = new Int_t[nCandidate]; + Int_t* sflag = new Int_t[nCandidate]; + + + //total energy in array + Float_t etbgTotal = 0.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 + + // Information extracted from fUnitArray + // load input vectors and calculate total energy in array + for(Int_t i=0; iAt(i); + if(uArray->GetUnitCutFlag()==1){ + etCell[i] = uArray->GetUnitEnergy(); + if (etCell[i] > 0.0) etCell[i] -= header->GetMinCellEt(); + if (etCell[i] < 0.0) etCell[i] = 0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; // default + } + else { + etCell[i] = 0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; + } + + if(uArray->GetUnitEnergy()>0.){ + ptT[loop1] = uArray->GetUnitEnergy(); + enT[loop1] = uArray->GetUnitEnergy(); + etaT[loop1] = uArray->GetUnitEta(); + phiT[loop1] = uArray->GetUnitPhi(); + detaT[loop1] = uArray->GetUnitDeta(); + dphiT[loop1] = uArray->GetUnitDphi(); + cFlagT[loop1]= uArray->GetUnitCutFlag(); + idT[loop1] = uArray->GetUnitID(); + if(cFlagT[loop1] == 1) { + hPtTotal->Fill(ptT[loop1]); + // fLego->Fill(etaT[i], phiT[i], ptT[i]); + etbgTotal+= ptT[loop1]; + } + loop1++; + } + } + + // fJets->SetNinput(nIn); + fJets->SetNinput(nCandidate); + + // calculate total energy and fluctuation in map + Double_t meanpt = hPtTotal->GetMean(); + Double_t ptRMS = hPtTotal->GetRMS(); + Double_t npart = hPtTotal->GetEntries(); + 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* 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]; + Int_t nJets; // to hold number of jets found by algorithm + Int_t nj; // number of jets accepted + Float_t prec = header->GetPrecBg(); + Float_t bgprec = 1; + while(bgprec > prec){ + //reset jet arrays in memory + memset(etaJet,0,sizeof(Float_t)*30); + memset(phiJet,0,sizeof(Float_t)*30); + memset(etJet,0,sizeof(Float_t)*30); + memset(etallJet,0,sizeof(Float_t)*30); + memset(etsigJet,0,sizeof(Float_t)*30); + memset(ncellsJet,0,sizeof(Int_t)*30); + memset(multJet,0,sizeof(Int_t)*30); + nJets = 0; + nj = 0; + // reset particles-jet array in memory + memset(injet,-1,sizeof(Int_t)*nCandidate); + //run cone algorithm finder + RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet); + //run background subtraction + if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event + nj = header->GetNAcceptJets(); + else + nj = nJets; + //subtract background + Float_t etbgTotalN = 0.0; //new background + if(header->GetBackgMode() == 1) // standar + SubtractBackg(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 2) //cone + SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 3) //ratio + SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 4) //statistic + SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + //calc precision + if(etbgTotalN != 0.0) + bgprec = (etbgTotal - etbgTotalN)/etbgTotalN; + else + bgprec = 0; + etbgTotal = etbgTotalN; // update with new background estimation + } //end while + + // add jets to list + Int_t* idxjets = new Int_t[nj]; + Int_t nselectj = 0; + printf("Found %d jets \n", nj); + + for(Int_t kj=0; kj (header->GetJetEtaMax())) || + (etaJet[kj] < (header->GetJetEtaMin())) || + (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin + Float_t px, py,pz,en; // convert to 4-vector + px = etJet[kj] * TMath::Cos(phiJet[kj]); + py = etJet[kj] * TMath::Sin(phiJet[kj]); + pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj]))); + en = TMath::Sqrt(px * px + py * py + pz * pz); + fJets->AddJet(px, py, pz, en); + AliAODJet jet(px, py, pz, en); + jet.Print(""); + + AddJet(jet); + + idxjets[nselectj] = kj; + nselectj++; + } + //add signal percentage and total signal in AliJets for analysis tool + Float_t* percentage = new Float_t[nselectj]; + Int_t* ncells = new Int_t[nselectj]; + Int_t* mult = new Int_t[nselectj]; + for(Int_t i = 0; i< nselectj; i++){ + percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]]; + ncells[i] = ncellsJet[idxjets[i]]; + mult[i] = multJet[idxjets[i]]; + } + //add particle-injet relationship /// + // for(Int_t bj = 0; bj < nIn; bj++){ + for(Int_t bj = 0; bj < nCandidate; bj++){ + if(injet[bj] == -1) continue; //background particle + Int_t bflag = 0; + for(Int_t ci = 0; ci< nselectj; ci++){ + if(injet[bj] == idxjets[ci]){ + injet[bj]= ci; + bflag++; + break; + } + } + if(bflag == 0) injet[bj] = -1; // set as background particle + } + fJets->SetNCells(ncells); + fJets->SetPtFromSignal(percentage); + fJets->SetMultiplicities(mult); + fJets->SetInJet(injet); + fJets->SetEtaIn(etaT); + fJets->SetPhiIn(phiT); + fJets->SetPtIn(ptT); + fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi())); + + + //delete + delete enT; + delete ptT; + delete etaT; + delete phiT; + delete detaT; + delete dphiT; + delete cFlagT; + delete cClusterT; + delete idT; + delete injet; + delete sflag; + delete hPtTotal; + delete etCell; + delete etaCell; + delete phiCell; + delete flagCell; + delete etaJet; + delete phiJet; + delete etJet; + delete etsigJet; + delete etallJet; + delete ncellsJet; + delete multJet; + delete idxjets; + delete percentage; + delete ncells; + delete mult; + + +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, + Int_t* flagCell, 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 nCell = nIn; + + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + + // Parameters from header + Float_t minmove = header->GetMinMove(); + Float_t maxmove = header->GetMaxMove(); + Float_t rc = header->GetRadius(); + Float_t etseed = header->GetEtSeed(); + + // tmp array of jets form algoritm + Float_t etaAlgoJet[30]; + Float_t phiAlgoJet[30]; + Float_t etAlgoJet[30]; + Int_t ncellsAlgoJet[30]; + + //run algorithm// + + // sort cells by et + Int_t * index = new Int_t[nCell]; + TMath::Sort(nCell, etCell, index); + + // variable used in centroide loop + Float_t eta = 0.0; + Float_t phi = 0.0; + Float_t eta0 = 0.0; + Float_t phi0 = 0.0; + Float_t etab = 0.0; + Float_t phib = 0.0; + Float_t etas = 0.0; + Float_t phis = 0.0; + Float_t ets = 0.0; + Float_t deta = 0.0; + Float_t dphi = 0.0; + Float_t dr = 0.0; + Float_t etsb = 0.0; + Float_t etasb = 0.0; + Float_t phisb = 0.0; + + + for(Int_t icell = 0; icell < nCell; icell++){ + Int_t jcell = index[icell]; + if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed + if(flagCell[jcell] != 0) continue; // if cell was used before + eta = etaCell[jcell]; + phi = phiCell[jcell]; + eta0 = eta; + phi0 = phi; + etab = eta; + phib = phi; + ets = etCell[jcell]; + etas = 0.0; + phis = 0.0; + etsb = ets; + etasb = 0.0; + phisb = 0.0; + for(Int_t kcell =0; kcell < nCell; kcell++){ + Int_t lcell = index[kcell]; + if(lcell == jcell) continue; // cell itself + if(flagCell[lcell] != 0) continue; // cell used before + if(etCell[lcell] > etCell[jcell]) continue; + //calculate dr + deta = etaCell[lcell] - eta; + dphi = phiCell[lcell] - phi; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ + // calculate offset from initiate cell + deta = etaCell[lcell] - eta0; + dphi = phiCell[lcell] - phi0; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + etas = etas + etCell[lcell]*deta; + phis = phis + etCell[lcell]*dphi; + ets = ets + etCell[lcell]; + //new weighted eta and phi including this cell + eta = eta0 + etas/ets; + phi = phi0 + phis/ets; + // if cone does not move much, just go to next step + dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib)); + if(dr <= minmove) break; + // cone should not move more than max_mov + dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets)); + if(dr > maxmove){ + eta = etab; + phi = phib; + ets = etsb; + etas = etasb; + phis = phisb; + }else{ // store this loop information + etab=eta; + phib=phi; + etsb = ets; + etasb = etas; + phisb = phis; + } + } + }//end of cells loop looking centroide + + //avoid cones overloap (to be implemented in the future) + + //flag cells in Rc, estimate total energy in cone + Float_t etCone = 0.0; + Int_t nCellIn = 0; + rc = header->GetRadius(); + for(Int_t ncell =0; ncell < nCell; ncell++){ + if(flagCell[ncell] != 0) continue; // cell used before + //calculate dr + deta = etaCell[ncell] - eta; + dphi = phiCell[ncell] - phi; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // cell in cone + flagCell[ncell] = -1; + etCone+=etCell[ncell]; + nCellIn++; + } + } + + // select jets with et > background + // estimate max fluctuation of background in cone + Double_t ncellin = (Double_t)nCellIn; + Double_t ntcell = (Double_t)nCell; + Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell); + // min cone et + Double_t etcmin = etCone ; // could be used etCone - etmin !! + //desicions !! etbmax < etcmin + for(Int_t mcell =0; mcell < nCell; mcell++){ + if(flagCell[mcell] == -1){ + if(etbmax < etcmin) + flagCell[mcell] = 1; //flag cell as used + else + flagCell[mcell] = 0; // leave it free + } + } + //store tmp jet info !!! + if(etbmax < etcmin) { + etaAlgoJet[nJets] = eta; + phiAlgoJet[nJets] = phi; + etAlgoJet[nJets] = etCone; + ncellsAlgoJet[nJets] = nCellIn; + nJets++; + } + + } // end of cells loop + + //reorder jets by et in cone + //sort jets by energy + Int_t * idx = new Int_t[nJets]; + TMath::Sort(nJets, etAlgoJet, idx); + for(Int_t p = 0; p < nJets; p++){ + etaJet[p] = etaAlgoJet[idx[p]]; + phiJet[p] = phiAlgoJet[idx[p]]; + etJet[p] = etAlgoJet[idx[p]]; + etallJet[p] = etAlgoJet[idx[p]]; + ncellsJet[p] = ncellsAlgoJet[idx[p]]; + } + + + //delete + delete index; + delete idx; + +} +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV2::SubtractBackg(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, + Int_t* multJet, Int_t* injet) +{ + //background subtraction using cone method but without correction in dE/deta distribution + + //calculate energy inside and outside cones + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Float_t rc= header->GetRadius(); + Float_t etIn[30]; + Float_t etOut = 0; + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if((fReader->GetCutFlag(jpart)) == 1){ // pt cut + etIn[ijet] += ptT[jpart]; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) + etOut += ptT[jpart]; // particle outside cones and pt cut + } //end particle loop + + //estimate jets and background areas + Float_t areaJet[30]; + Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); + for(Int_t k=0; k header->GetLegoEtaMax()){ // sector outside etamax + Float_t h = header->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < header->GetLegoEtaMin()){ // sector outside etamin + Float_t h = header->GetLegoEtaMax() + etaJet[k]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin; + areaOut = areaOut - areaJet[k]; + } + //subtract background using area method + for(Int_t ljet=0; ljetGetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV2::SubtractBackgStat(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, + Int_t* multJet, Int_t* injet) +{ + + //background subtraction using statistical method + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Float_t etbgStat = header->GetBackgStat(); // pre-calculated background + + //calculate energy inside + Float_t rc= header->GetRadius(); + Float_t etIn[30]; + + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if((fReader->GetCutFlag(jpart)) == 1){ // pt cut + etIn[ijet]+= ptT[jpart]; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + } //end particle loop + + //calc jets areas + Float_t areaJet[30]; + Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); + for(Int_t k=0; k header->GetLegoEtaMax()){ // sector outside etamax + Float_t h = header->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < header->GetLegoEtaMin()){ // sector outside etamin + Float_t h = header->GetLegoEtaMax() + etaJet[k]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin; + } + + //subtract background using area method + for(Int_t ljet=0; ljetGetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + + // jet energy and area arrays + TH1F* hEtJet[30]; + TH1F* hAreaJet[30]; + for(Int_t mjet=0; mjet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + injet[jpart] = ijet; + multJet[ijet]++; + if((fReader->GetCutFlag(jpart)) == 1){// pt cut + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) + hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + } //end particle loop + + //calc areas + Float_t eta0 = etamin; + Float_t etaw = (etamax - etamin)/((Float_t)ndiv); + Float_t eta1 = eta0 + etaw; + for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins + Float_t etac = eta0 + etaw/2.0; + Float_t areabg = etaw*2.0*TMath::Pi(); + for(Int_t ijet=0; ijet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + hAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = hAreaBackg->GetBinContent(bin); + Float_t etb = hEtBackg->GetBinContent(bin); + Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction + } + } + } + + // calc background total + Double_t etOut = hEtBackg->Integral(); + Double_t areaOut = hAreaBackg->Integral(); + Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetGetBackgCutRatio(); + + + //general + Float_t rc= header->GetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + + // jet energy and area arrays + TH1F* hEtJet[30]; + TH1F* hAreaJet[30]; + for(Int_t mjet=0; mjetGetCutFlag(jpart)) != 1) continue; + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if((fReader->GetCutFlag(jpart)) == 1){ //pt cut + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + } //end particle loop + + //calc areas + Float_t eta0 = etamin; + Float_t etaw = (etamax - etamin)/((Float_t)ndiv); + Float_t eta1 = eta0 + etaw; + for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins + Float_t etac = eta0 + etaw/2.0; + Float_t areabg = etaw*2.0*TMath::Pi(); + for(Int_t ijet=0; ijet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + hAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = hAreaBackg->GetBinContent(bin); + Float_t etb = hEtBackg->GetBinContent(bin); + Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction + } + } + } + + // calc background total + Double_t etOut = hEtBackg->Integral(); + Double_t areaOut = hAreaBackg->Integral(); + Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetReset(); + fJets->ClearJets(); + AliJetFinder::Reset(); +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV2::WriteJHeaderToFile() +{ + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + header->Write(); +} + +//////////////////////////////////////////////////////////////////////// + +void AliUA1JetFinderV2::Init() +{ + // initializes some variables + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + // book lego + fLego = new + TH2F("legoH","eta-phi", + header->GetLegoNbinEta(), header->GetLegoEtaMin(), + header->GetLegoEtaMax(), header->GetLegoNbinPhi(), + header->GetLegoPhiMin(), header->GetLegoPhiMax()); + + fDebug = fReader->GetReaderHeader()->GetDebug(); + fOpt = fReader->GetReaderHeader()->GetDetector(); + + if(fOpt>0) + fReader->CreateTasks(); + +} diff --git a/JETAN/AliUA1JetFinderV2.h b/JETAN/AliUA1JetFinderV2.h new file mode 100644 index 00000000000..b3823458681 --- /dev/null +++ b/JETAN/AliUA1JetFinderV2.h @@ -0,0 +1,69 @@ +#ifndef ALIUA1JETFINDERV2_H +#define ALIUA1JETFINDERV2_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + +//--------------------------------------------------------------------- +// UA1 Cone Algorithm Finder V1 +// manages the search for jets +// Author: Rafael.Diaz.Valdes@cern.ch +// (version in c++) +// Modified to include neutral particles (magali.estienne@ires.in2p3.fr) +//--------------------------------------------------------------------- + +#include "AliJetFinder.h" +class AliUA1JetHeaderV1; +class TH2F; + +class AliUA1JetFinderV2 : public AliJetFinder +{ + public: + + AliUA1JetFinderV2(); + ~AliUA1JetFinderV2(); + + // others + void FindJets(); + + void RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, + Int_t* flagCell, 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 SubtractBackg(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,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* etJet,Float_t* etaJet, 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* etJet,Float_t* etaJet, 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* etJet,Float_t* etaJet, Float_t* phiJet, + Float_t* etsigJet, Int_t* multJet, Int_t* injet); + void Reset(); + void Init(); + void WriteJHeaderToFile(); + + protected: + AliUA1JetFinderV2(const AliUA1JetFinderV2& rJetF1); + AliUA1JetFinderV2& operator = (const AliUA1JetFinderV2& rhsf); + TH2F * fLego; //Lego Histo + Int_t fDebug; + Int_t fOpt; + + ClassDef(AliUA1JetFinderV2,1) +}; + +#endif diff --git a/JETAN/AliUA1JetHeaderV1.cxx b/JETAN/AliUA1JetHeaderV1.cxx index 3933bfca88f..73b1f156826 100644 --- a/JETAN/AliUA1JetHeaderV1.cxx +++ b/JETAN/AliUA1JetHeaderV1.cxx @@ -35,6 +35,7 @@ AliUA1JetHeaderV1::AliUA1JetHeaderV1(): fConeRadius(0.3), fEtSeed(3.0), fMinJetEt(10.), + fMinCellEt(0.0), // Temporarily added -> To be removed if not necessary fMinMove(0.05), fMaxMove(0.15), fBackgMode(1), diff --git a/JETAN/AliUA1JetHeaderV1.h b/JETAN/AliUA1JetHeaderV1.h index 69f2adffc0e..f54eba98d99 100644 --- a/JETAN/AliUA1JetHeaderV1.h +++ b/JETAN/AliUA1JetHeaderV1.h @@ -26,6 +26,7 @@ class AliUA1JetHeaderV1 : public AliJetHeader Float_t GetMaxMove() const {return fMaxMove;} Float_t GetEtSeed() const {return fEtSeed;} Float_t GetMinJetEt() const {return fMinJetEt;} + Float_t GetMinCellEt() const {return fMinCellEt;} // Added temporarily !!! To be removed if not necessary Int_t GetLegoNbinEta() const {return fLegoNbinEta;} Int_t GetLegoNbinPhi() const {return fLegoNbinPhi;} Float_t GetLegoEtaMin() const {return fLegoEtaMin;} @@ -45,6 +46,7 @@ class AliUA1JetHeaderV1 : public AliJetHeader void SetMaxMove(Float_t f) {fMaxMove=f;} void SetEtSeed(Float_t f) {fEtSeed=f;} void SetMinJetEt(Float_t f) {fMinJetEt=f;} + void SetMinCellEt(Float_t f) {fMinCellEt=f;} void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;} void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;} void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;} @@ -66,6 +68,7 @@ protected: Float_t fConeRadius; // Cone radius Float_t fEtSeed; // Min. Et for seed Float_t fMinJetEt; // Min Et of jet + Float_t fMinCellEt; // Min Et in one cell // parameters of backgound substraction Float_t fMinMove; // min cone move Float_t fMaxMove; // max cone move diff --git a/JETAN/JETANLinkDef.h b/JETAN/JETANLinkDef.h index e4668bb4b38..0731d4f3ab2 100644 --- a/JETAN/JETANLinkDef.h +++ b/JETAN/JETANLinkDef.h @@ -19,6 +19,7 @@ #pragma link C++ class AliJetAnalysis+; #pragma link C++ class AliJetDistributions+; #pragma link C++ class AliUA1JetFinderV1+; +#pragma link C++ class AliUA1JetFinderV2+; #pragma link C++ class AliUA1JetHeaderV1+; #pragma link C++ class AliJetGrid+; #pragma link C++ class AliJetUnitArray+; @@ -26,6 +27,7 @@ #pragma link C++ class AliJetHadronCorrectionv0+; #pragma link C++ class AliJetHadronCorrectionv1+; #pragma link C++ class AliJetFillUnitArrayTracks+; +#pragma link C++ class AliJetFillUnitArrayEMCalDigits+; #pragma link C++ class AliJetDummyGeo+; #pragma link C++ class AliAnalysisTaskJets+; #pragma link C++ class AliDAJetFinder+; diff --git a/JETAN/libJETAN.pkg b/JETAN/libJETAN.pkg index 1fe5d253d32..98100120f25 100644 --- a/JETAN/libJETAN.pkg +++ b/JETAN/libJETAN.pkg @@ -7,10 +7,10 @@ SRCS = AliJet.cxx AliJetHeader.cxx \ AliLeading.cxx \ AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \ AliJetAnalysis.cxx AliJetDistributions.cxx \ - AliUA1JetFinderV1.cxx AliUA1JetHeaderV1.cxx \ + AliUA1JetFinderV1.cxx AliUA1JetFinderV2.cxx AliUA1JetHeaderV1.cxx \ AliJetGrid.cxx AliJetUnitArray.cxx AliJetHadronCorrection.cxx \ AliJetHadronCorrectionv0.cxx AliJetHadronCorrectionv1.cxx \ - AliJetFillUnitArrayTracks.cxx \ + AliJetFillUnitArrayTracks.cxx AliJetFillUnitArrayEMCalDigits.cxx\ AliJetDummyGeo.cxx \ AliJetFinderTypes.cxx \ AliAnalysisTaskJets.cxx \ -- 2.43.0