- Hadron correction.
(Magali Estienne)
#include <TSystem.h>
#include <TLorentzVector.h>
#include <TVector3.h>
+#include <TChain.h>
+
#include "AliJetAODReader.h"
#include "AliJetAODReaderHeader.h"
#include "AliAODEvent.h"
class TClonesArray;
class TH1I;
class TH1D;
+class TH1F;
+class TH1;
+
class AliJetReader;
class AliJet;
// 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);
* 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 <assert.h>
+
+// --- Root header files ---
+#include <TVector3.h>
+#include <TGeoManager.h>
+#include <TGeoMatrix.h>
+#include <TGeoNode.h>
+
+// --- 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<fNumberOfSuperModules; i++) fMatrixOfSM[i] = 0;
+
+ fCentersOfCellsEtaDir.Set(fNZ*fNETAdiv);
+ fCentersOfCellsXDir.Set(fNZ*fNETAdiv);
+ fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
+ fEtaCentersOfCells.Set(fNZ*fNETAdiv*fNPhi*fNPHIdiv);
+ fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
+
+ int nphism = GetNumberOfSuperModules()/2;
+ double dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
+ double rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
+ double phi, phiRad, xpos, ypos, zpos;
+ for(int i=0; i<nphism; i++){
+ phi = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
+ phiRad = phi*TMath::Pi()/180.;
+ xpos = rpos * TMath::Cos(phiRad);
+ ypos = rpos * TMath::Sin(phiRad);
+ zpos = fParSM[2];
+ if(i==5) {
+ xpos += (fParSM[1]/2. * TMath::Sin(phiRad));
+ ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
+ }
+ // pozitive z
+ int ind = 2*i;
+ TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
+ fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+ xpos,ypos, zpos, geoRot0);
+ // negaive z
+ ind++;
+ double phiy = 90. + phi + 180.;
+ if(phiy>=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<fNumberOfSuperModules; i++) fMatrixOfSM[i] = 0;
+
+ fCentersOfCellsEtaDir.Set(fNZ*fNETAdiv);
+ fCentersOfCellsXDir.Set(fNZ*fNETAdiv);
+ fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
+ fEtaCentersOfCells.Set(fNZ*fNETAdiv*fNPhi*fNPHIdiv);
+ fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
+
+ int nphism = GetNumberOfSuperModules()/2;
+ double dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
+ double rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
+ double phi, phiRad, xpos, ypos, zpos;
+ for(int i=0; i<nphism; i++){
+ phi = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
+ phiRad = phi*TMath::Pi()/180.;
+ xpos = rpos * TMath::Cos(phiRad);
+ ypos = rpos * TMath::Sin(phiRad);
+ zpos = fParSM[2];
+ if(i==5) {
+ xpos += (fParSM[1]/2. * TMath::Sin(phiRad));
+ ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
+ }
+ // pozitive z
+ int ind = 2*i;
+ TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
+ fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+ xpos,ypos, zpos, geoRot0);
+ // negaive z
+ ind++;
+ double phiy = 90. + phi + 180.;
+ if(phiy>=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<nphi; i++) {
+ d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
+ if(d < dmin) {
+ dmin = d;
+ iphi = i;
+ }
+ }
+ // odd SM are turned with respect of even SM - reverse indexes
+
+ // eta index
+ absEta = TMath::Abs(eta);
+ etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
+ dmin = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
+ ieta = 0;
+ for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
+ d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
+ if(d < dmin) {
+ dmin = d;
+ ieta = i;
+ }
+ }
+
+ if(eta<0) iphi = (nphi-1) - iphi;
+ absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
+
+ return kTRUE;
+ }
+ return kFALSE;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
+{
+ // Return false if phi belongs a phi cracks between SM
+
+ static Int_t i;
+
+ if(TMath::Abs(eta) > 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; i<tn->GetNdaughters(); 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; i<xen1->GetNdaughters(); i++) {
+ TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
+ fMatrixOfSM[i] = sm->GetMatrix();
+ }
+ printf("transInit %d: ", transInit);
+ transInit = kTRUE;
+}
+
/* 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 <TObject.h>
+#include <TArrayD.h>
+#include <TMath.h>
+#include <TVector3.h>
+
+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
// Magali Estienne <magali.estienne@IReS.in2p3.fr>
//-------------------------------------------------------------------------
-
+// --- Standard library ---
#include <Riostream.h>
+
+// --- ROOT system ---
#include <TSystem.h>
+#include <TStopwatch.h>
#include <TLorentzVector.h>
#include <TVector3.h>
+#include <TTask.h>
#include <TGeoManager.h>
+#include <assert.h>
+#include <TRefArray.h>
+#include <TMath.h>
+#include <TChain.h>
+
+// --- 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
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
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();
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);
}
}
+//____________________________________________________________________________
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
// 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();
// 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; i<nEntRef; i++)
+ {
+ ((AliJetUnitArray*)fRefArray->At(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
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(nBin<fNumUnits)
+ {
+ fGeom->EtaPhiFromIndex(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 && nBin<fNumUnits+nElements){
+ fTpcGrid->GetEtaPhiFromIndex2(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 && nBin<fNumUnits+nElements+nGaps){
+ if(nBin<fNumUnits+nElements+n0)
+ {
+ Float_t phi = eta = 0.;
+ fGrid0->GetEtaPhiFromIndex2(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 && nBin<fNumUnits+nElements+n0+n1)
+ {
+ Float_t phi = eta = 0.;
+ fGrid1->GetEtaPhiFromIndex2(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 && nBin<fNumUnits+nElements+n0+n1+n2)
+ {
+ Float_t phi = eta = 0.;
+ fGrid2->GetEtaPhiFromIndex2(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 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
+ {
+ Float_t phi = eta = 0.;
+ fGrid3->GetEtaPhiFromIndex2(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 && nBin<fNumUnits+nElements+nGaps)
+ {
+ Float_t phi = eta = 0.;
+ fGrid4->GetEtaPhiFromIndex2(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;
}
+
#include "AliJetUnitArray.h"
#include "AliJetGrid.h"
class AliJetESDReaderHeader;
-class AliEMCALGeometry;
class AliJetDummyGeo;
class AliJetHadronCorrection;
class AliJetUnitArray;
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;}
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)
};
--- /dev/null
+
+/**************************************************************************
+ * 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 <Riostream.h>
+
+// --- ROOT system ---
+#include <TSystem.h>
+#include <TLorentzVector.h>
+#include <TRefArray.h>
+#include "TTask.h"
+#include <TGeoManager.h>
+#include <TMath.h>
+#include <TClonesArray.h>
+#include <TArrayS.h>
+
+// --- 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; k<nD; k++) {
+
+ // Convert energy in GeV
+ Int_t idF = (Int_t)digID->At(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()<ptMin)
+ uArray->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));
+
+
+}
--- /dev/null
+#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
-
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//======================================================================
-// *** November 2006
-// Author: magali.estienne@ires.in2p3.fr
-// 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and
-// EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance
-// 2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment)
-// fGrid==1, work with a part of the TPC acceptance
-// 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them
-// 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// ***September 2006
-// TTask : Fill Unit Array for the Tracks information
-// Called by ESD reader for jet analysis
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// *** July 2006
-// 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex
-// and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the
-// index or the eta, phi position of a grid.
-// 2) Define a grid for TPC
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// ***July 2006
+
+//----------------------------------------------------------------------
// Fill Unit Array class
-// Class used by AliJetESDReader to fill a UnitArray from the information extracted
-// from the particle tracks
+// Class used by AliJetESDReader to fill a UnitArray from the information
+// extracted from the particle tracks
// Author: magali.estienne@ires.in2p3.fr
-//======================================================================
+//----------------------------------------------------------------------
// --- Standard library ---
// --- ROOT system ---
#include <TSystem.h>
#include <TLorentzVector.h>
+#include <TRefArray.h>
#include <TVector3.h>
-//#include "Math/Vector3D.h"
-//#include "Math/Vector3Dfwd.h"
#include "TTask.h"
#include <TGeoManager.h>
#include <TMatrixD.h>
#include <TArrayD.h>
+#include <TMath.h>
+#include <TClonesArray.h>
// --- AliRoot header files ---
#include "AliJetFinder.h"
#include "AliJetReader.h"
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
-#include "AliESD.h"
-//#include "AliEMCALGeometry.h"
+//#include "AliESD.h"
+#include "AliESDEvent.h"
#include "AliJetDummyGeo.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayTracks.h"
//_____________________________________________________________________________
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();
fPhiMinCal = fGeom->GetArm1PhiMin();
fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
- if(fDebug>30){
- cout << "fEtaMinCal : " << fEtaMinCal << endl;
- cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
- cout << "fPhiMinCal : " << fPhiMinCal << endl;
- cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
- }
-
fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
fPhiMax,fEtaMin,fEtaMax);
fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
}
//_____________________________________________________________________________
-void AliJetFillUnitArrayTracks::Exec(Option_t* option)
+void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
{
//
// Main method.
- // Explain
+ //
fDebug = fReaderHeader->GetDebug();
- if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");
- if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
- // Set EMCal Geometry
- SetEMCALGeometry();
+
// Set parameters
InitParameters();
- TClonesArray *lvArray = fMomentumArray; // Correct checked !
- Int_t nInT = lvArray->GetEntries(); // Correct checked !
- Float_t ptMin = fReaderHeader->GetPtCut(); // Correct checked !
-
- // sflag -> not yet implemented !!!!
-
- if(fDebug>3) cout << "nInT : " << nInT << endl;
-
- if (nInT == 0) return;
-
- // local arrays for input
- Float_t* enT = new Float_t[nInT];
- Float_t* ptT = new Float_t[nInT];
- Float_t* etaT = new Float_t[nInT];
- Float_t* phiT = new Float_t[nInT];
- Float_t* thetaT = new Float_t[nInT];
-
- Int_t trackInEmcalAcc = 0;
- Int_t trackInTpcAcc = 0;
- Int_t trackInTpcAccOnly = 0;
-
- Int_t nElements = fTPCGrid->GetNEntries();
- Int_t nElements2 = fEMCalGrid->GetNEntries();
- fGrid = fTPCGrid->GetGridType();
+ // get number of tracks in event (for the loop)
+ Int_t goodTrack = 0;
+ Int_t nt = 0;
+ Float_t pt, eta,phi;
+ TVector3 p3;
+ nt = fESD->GetNumberOfTracks();
+ if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
+
+ // temporary storage of signal and pt cut flag
+ Int_t* sflag = new Int_t[nt];
+ Int_t* cflag = new Int_t[nt];
+
+ // get cuts set by user
+ Float_t ptMin = fReaderHeader->GetPtCut();
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
+ fOpt = fReaderHeader->GetDetector();
+ fDZ = fReaderHeader->GetDZ();
- if(fDebug>3){
- cout << "nElements : " << nElements << endl;
- cout << "nElements2 : " << nElements2 << endl;
- cout << "fNumUnits : " << fNumUnits << endl;
- cout << "sum : " << fNumUnits+nElements << endl;
- }
+ Int_t nTracksEmcal = 0;
+ Int_t nTracksEmcalDZ = 0;
+ Int_t nTracksTpc = 0;
+ Int_t nTracksTpcOnly = 0;
+ Int_t nTracksEmcalCut = 0;
+ Int_t nTracksEmcalDZCut = 0;
+ Int_t nTracksTpcCut = 0;
+ Int_t nTracksTpcOnlyCut = 0;
- // Set energy exactly to zero
- if(fGrid==0)
- for(Int_t k=0; k<nElements; k++)
- fUnitArray[k].SetUnitEnergy(0.);
-
- if(fGrid==1)
- for(Int_t k=0; k<fNumUnits+nElements; k++)
- fUnitArray[k].SetUnitEnergy(0.);
-
- // load input vectors
- for (Int_t i = 0; i < nInT; i++)
- {
- TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
- enT[i] = lv->Energy();
- ptT[i] = lv->Pt();
- etaT[i] = lv->Eta();
- phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
- thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
-
- if(fDebug>20){
- cout << "enT[" << i << "] : " << enT[i] << endl;
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "phiT[" << i << "] : " << phiT[i] << endl;
- cout << "thetaT[" << i << "] : " << thetaT[i] << endl;
- cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
- cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
- cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
- cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
- }
-
- if(fGrid==0)
- {
- // For the moment, only TPC filled from its grid in its total acceptance
- if(fDebug>2)
- cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
-
- trackInTpcAccOnly += 1;
-
- Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
-
- Float_t unitEnergy = 0.;
- unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
-
- if(unitEnergy > 0. && fDebug >10){
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "TPC id : " << idTPC << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "phiT[" << i << "] : " << phiT[i] << endl;
- cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl;
- cout << "fUnitArray[idTPC].GetUnitEnergy(): " <<
- fUnitArray[idTPC].GetUnitEnergy() << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
+ fGrid = fTPCGrid->GetGridType();
- // Fill energy in TPC acceptance
- fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]);
- if(fDebug > 10){
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "unitEnergy in TPC acceptance after : " <<
- fUnitArray[idTPC].GetUnitEnergy() << endl;
- }
-
- // Pt cut flag
- if(fUnitArray[idTPC].GetUnitEnergy()<ptMin)
- fUnitArray[idTPC].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[idTPC].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[idTPC].SetUnitDetectorFlag(kTpc);
+ //loop over tracks
+ for (Int_t it = 0; it < nt; it++) {
+ AliESDtrack *track = fESD->GetTrack(it);
+ UInt_t status = track->GetStatus();
+
+ Double_t mom[3];
+ track->GetPxPyPz(mom);
+ p3.SetXYZ(mom[0],mom[1],mom[2]);
+ pt = p3.Pt();
+ Float_t mass = 0.;
+ mass = track->GetMass();
+
+ if (((status & AliESDtrack::kITSrefit) == 0) ||
+ ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check
+ if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
+ && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
+ if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
+ && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
+ eta = p3.Eta();
+ phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
+
+ if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
+
+ // sflag -> not yet implemented !!!!
+
+ if(fGrid==0)
+ {
+ // Only TPC filled from its grid in its total acceptance
+
+ Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+ Bool_t ok = kFALSE;
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+ uArray->SetUnitTrackID(it);
+
+ Float_t unitEnergy = 0.;
+ unitEnergy = uArray->GetUnitEnergy();
+ if(unitEnergy==0.){
+ nTracksTpcOnly++;
+ ok = kTRUE;
+ }
+ // Fill energy in TPC acceptance
+ uArray->SetUnitEnergy(unitEnergy + pt);
+ uArray->SetUnitPxPyPz(mom);
+ uArray->SetUnitMass(mass);
+
+ // Pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin){
+ uArray->SetUnitCutFlag(kPtSmaller);
+ }
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok) nTracksTpcOnlyCut++;
}
- if(fGrid==1)
- {
- // Fill track information in EMCAL acceptance
- if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) &&
- (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// &&
- // GetCutFlag(i) == 1)
- // ptT[i] > ptMin)
- {
- trackInEmcalAcc += 1;
-
- if(fDebug>20){
- cout << "before : " << endl;
- cout << "etaT[i] : " << etaT[i] << endl;
- cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
- }
-
- // This function should be modified soon
-// Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete
- Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified
-// Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]); // Using an EMCal grid -> calculated
-// Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]); // Using an EMCal grid -> tabulated (faster)
-
- Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy();
-
- if(fDebug>20) {
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "after : " << endl;
- cout << "towerID : " << towerID << endl;
- cout << "etaT[i] : " << etaT[i] << endl;
- cout << "phiT[i](rad) : " << phiT[i] << endl;
- cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
- cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl;
- cout << "fHadCorr : " << fHadCorr << endl;
- cout << "fHCorrection : " << fHCorrection << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
-
- //OLD WAY: //Do Hadron Correction
- if (fHCorrection != 0)
- {
- Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
- unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]);
- if(fDebug>20){
- cout << "Inside loop for hadron correction ==========================" << endl;
- cout << "enT[" << i << "] : " << enT[i] << endl;
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "Energy correction : " << hCEnergy << endl;
- cout << "unitEnergy : " << unitEnergy << endl;
- cout << "fUnitArray[towerID].GetUnitEnergy() : " <<
- fUnitArray[towerID].GetUnitEnergy() << endl;
- }
- } //end Hadron Correction loop
-
- fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]);
-
- // Put a pt cut flag
- if(fUnitArray[towerID].GetUnitEnergy()<ptMin)
- fUnitArray[towerID].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[towerID].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[towerID].SetUnitDetectorFlag(kTpc);
-
- if(fDebug>10){
- cout << "After pT filled ===============================================" << endl;
- cout << "PtCut : " << ptMin << endl;
- cout << "ptT[" << i << "] : " << ptT[i] << endl;
- cout << "unitEnergy : " << unitEnergy << endl;
- cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl;
- cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl;
- cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl;
- cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl;
- }
- } // end loop on EMCal acceptance cut + tracks quality
- else{
- // Outside EMCal acceptance
- if(fDebug>2)
- cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
-
- trackInTpcAcc += 1;
-
- // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
- Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
-
- Float_t unitEnergy2 = 0.;
- unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
-
- if(unitEnergy2 > 0. && fDebug >10){
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "TPC id : " << idTPC << endl;
- cout << "Total id : " << fNumUnits-1+idTPC << endl;
- cout << "etaT[" << i << "] : " << etaT[i] << endl;
- cout << "phiT[" << i << "] : " << phiT[i] << endl;
- cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl;
- cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " <<
- fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
- }
+ // Detector flag
+ if(unitEnergy>0) {
+ uArray->SetUnitDetectorFlag(kAll);
+ }
+ if(uArray->GetUnitEnergy()>0){
+ fRefArray->Add(uArray);
+ }
- // Fill energy outside emcal acceptance
- fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
- // Pt cut flag
- if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()<ptMin)
- fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtSmaller);
- else fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtHigher);
- // Detector flag
- fUnitArray[fNumUnits-1+idTPC].SetUnitDetectorFlag(kTpc);
- }
- } // end fGrid==1
- } // end loop on entries (tpc tracks)
-
- if(fDebug>3){
- printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc);
- printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc);
- }
-
- if(fGrid==0) fNIn = trackInTpcAccOnly;
- if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc;
-
- fOpt = fReaderHeader->GetDetector();
-
- if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt);
-
- if(fOpt==1 || option=="tpc") // if only TPC
- { //Set all unit flags, Eta, Phi
-
- if(fGrid==0)
- {
- if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
- for(Int_t i=0; i<nElements; i++)
- {
- Float_t eta = 0.;
- Float_t phi = 0.;
- if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
- fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
- fUnitArray[i].SetUnitEntries(nElements);
- fUnitArray[i].SetUnitID(i);
- fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta);
- fUnitArray[i].SetUnitEta(eta);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
+ }
+
+ if(fGrid==1)
+ {
+ Int_t nElements = fTPCGrid->GetNEntries();
+
+ // Fill track information in EMCAL acceptance
+ if((eta >= fEtaMin && eta <= fEtaMax) &&
+ (phi >= fPhiMin && phi <= fPhiMax))// &&
+ {
+
+ // Include dead-zones
+ if(fDZ)
+ {
+ Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
+ Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
+ fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
+ fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+ fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+ fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+ fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+ Int_t n0 = fGrid0->GetNEntries();
+ Int_t n1 = fGrid1->GetNEntries();
+ Int_t n2 = fGrid2->GetNEntries();
+ Int_t n3 = fGrid3->GetNEntries();
+
+ if(phi >= phimin0 && phi <= phimax0){
+ Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
+ AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+ uArray0->SetUnitTrackID(it);
+ Float_t uEnergy0 = uArray0->GetUnitEnergy();
+ Bool_t ok0 = kFALSE;
+ if(uEnergy0==0.){
+ nTracksEmcalDZ++;
+ ok0 = kTRUE;
+ }
+ uArray0->SetUnitEnergy(uEnergy0+pt);
+ if(uArray0->GetUnitEnergy()<ptMin)
+ uArray0->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray0->SetUnitCutFlag(kPtHigher);
+ if(ok0) nTracksEmcalDZCut++;
+ }
+ if(uArray0->GetUnitEnergy()>0)
+ fRefArray->Add(uArray0);
}
- }
-
- if(fGrid==1)
- {
- if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
- for(Int_t i=0; i<fNumUnits+nElements; i++)
- {
- Float_t eta = 0.;
- Float_t phi = 0.;
- if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
- //Set all units to be outside a jet initially
- fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
- fUnitArray[i].SetUnitEntries(fNumUnits+nElements);
- fUnitArray[i].SetUnitID(i);
-
- if(fUnitArray[i].GetUnitID()<fNumUnits)
- {
- // fGeom->EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63
- fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry
- phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
- // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid
- // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid
- fUnitArray[i].SetUnitEta(eta);
- //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi());
- }
+ if(phi >= phimin1 && phi <= phimax1){
+ Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
+ AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
+ uArray1->SetUnitTrackID(it);
+ Float_t uEnergy1 = uArray1->GetUnitEnergy();
+ Bool_t ok1 = kFALSE;
+ if(uEnergy1==0.){
+ nTracksEmcalDZ++;
+ ok1 = kTRUE;
+ }
+ uArray1->SetUnitEnergy(uEnergy1+pt);
+ if(uArray1->GetUnitEnergy()<ptMin)
+ uArray1->SetUnitCutFlag(kPtSmaller);
else {
- fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta);
- fUnitArray[i].SetUnitEta(eta);
- fUnitArray[i].SetUnitPhi(phi);
- fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
- fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
-
- if(fDebug>10)
- {
- if(fUnitArray[i].GetUnitEnergy()!=0.){
- cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " <<
- fUnitArray[i].GetUnitID()+1-fNumUnits << endl;
- cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " <<
- fUnitArray[i].GetUnitEnergy() << endl;
- cout << "(fUnitArray[" << i << "].GetUnitEta() : " <<
- fUnitArray[i].GetUnitEta() << endl;
- cout << "(fUnitArray[" << i << "].GetUnitPhi() : " <<
- fUnitArray[i].GetUnitPhi() << endl;
- }
- }
+ uArray1->SetUnitCutFlag(kPtHigher);
+ if(ok1) nTracksEmcalDZCut++;
}
- fUnitArray[i].SetUnitClusterID(0);
- }//end loop over all units in array (same as all towers in EMCAL)
- }
-
- if(fDebug>20)
- {
- for(Int_t i=0; i<fNumUnits+nElements; i++)
- {
- if(fUnitArray[i].GetUnitEnergy()!=0.){
- cout << "######################################################### " << endl;
- cout << "Final UnitArray filled with energy != 0" << i << endl;
- cout << "Pointeur UnitArray : " << fUnitArray << " ID : " <<
- fUnitArray[i].GetUnitID() << " Energy : " << fUnitArray[i].GetUnitEnergy() <<
- " Eta : " << fUnitArray[i].GetUnitEta() << " Phi : " << fUnitArray[i].GetUnitPhi() << endl;
+ if(uArray1->GetUnitEnergy()>0)
+ fRefArray->Add(uArray1);
+ }
+ if(phi >= phimin2 && phi <= phimax2){
+ Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
+ AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
+ uArray2->SetUnitTrackID(it);
+ Float_t uEnergy2 = uArray2->GetUnitEnergy();
+ Bool_t ok2 = kFALSE;
+ if(uEnergy2==0.){
+ nTracksEmcalDZ++;
+ ok2 = kTRUE;
+ }
+ uArray2->SetUnitEnergy(uEnergy2+pt);
+ if(uArray2->GetUnitEnergy()<ptMin)
+ uArray2->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray2->SetUnitCutFlag(kPtHigher);
+ if(ok2) nTracksEmcalDZCut++;
+ }
+ if(uArray2->GetUnitEnergy()>0)
+ fRefArray->Add(uArray2);
+ }
+ if(phi >= phimin3 && phi <= phimax3){
+ Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
+ AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
+ uArray3->SetUnitTrackID(it);
+ Float_t uEnergy3 = uArray3->GetUnitEnergy();
+ Bool_t ok3 = kFALSE;
+ if(uEnergy3==0.){
+ nTracksEmcalDZ++;
+ ok3 = kTRUE;
+ }
+ uArray3->SetUnitEnergy(uEnergy3+pt);
+ if(uArray3->GetUnitEnergy()<ptMin)
+ uArray3->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray3->SetUnitCutFlag(kPtHigher);
+ if(ok3) nTracksEmcalDZCut++;
+ }
+ if(uArray3->GetUnitEnergy()>0)
+ fRefArray->Add(uArray3);
+ }
+ if(phi >= phimin4 && phi <= phimax4){
+ Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
+ AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
+ uArray4->SetUnitTrackID(it);
+ Float_t uEnergy4 = uArray4->GetUnitEnergy();
+ Bool_t ok4 = kFALSE;
+ if(uEnergy4==0.){
+ nTracksEmcalDZ++;
+ ok4 = kTRUE;
+ }
+ uArray4->SetUnitEnergy(uEnergy4+pt);
+ if(uArray4->GetUnitEnergy()<ptMin)
+ uArray4->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray4->SetUnitCutFlag(kPtHigher);
+ if(ok4) nTracksEmcalDZCut++;
}
+ if(uArray4->GetUnitEnergy()>0)
+ fRefArray->Add(uArray4);
}
+ } // end fDZ
+
+ Int_t towerID = 0;
+ fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
+
+ if(towerID==-1) continue;
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
+ uArray->SetUnitTrackID(it);
+ Float_t unitEnergy = uArray->GetUnitEnergy();
+ Bool_t ok = kFALSE;
+ if(unitEnergy==0.){
+ nTracksEmcal++;
+ ok=kTRUE;
+ }
+
+ // Do Hadron Correction
+ // Parametrization to be added
+ if (fHCorrection != 0)
+ {
+ // Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
+ Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
+ unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
+
+ } //end Hadron Correction loop
+
+ uArray->SetUnitEnergy(unitEnergy + pt);
+
+ // Put a pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin)
+ uArray->SetUnitCutFlag(kPtSmaller);
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok) nTracksEmcalCut++;
}
+ // Detector flag
+ if(unitEnergy > 0)
+ uArray->SetUnitDetectorFlag(kAll);
+
+ if(uArray->GetUnitEnergy()>0)
+ fRefArray->Add(uArray);
+
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
+
+ } // end loop on EMCal acceptance cut + tracks quality
+ else{
+ // Outside EMCal acceptance
+
+ // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
+ Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
+ uArray->SetUnitTrackID(it);
+
+ Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
+ Bool_t ok2 = kFALSE;
+ if(unitEnergy2==0.){
+ nTracksTpc++;
+ ok2=kTRUE;
+ }
+ // Fill energy outside emcal acceptance
+ uArray->SetUnitEnergy(unitEnergy2 + pt);
+
+ // Pt cut flag
+ if(uArray->GetUnitEnergy()<ptMin){
+ uArray->SetUnitCutFlag(kPtSmaller);
+ }
+ else {
+ uArray->SetUnitCutFlag(kPtHigher);
+ if(ok2) nTracksTpcCut++;
+ }
+ // Detector flag
+ if(unitEnergy2 > 0)
+ uArray->SetUnitDetectorFlag(kTpc);
+ if(uArray->GetUnitEnergy()>0)
+ fRefArray->Add(uArray);
+
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
+ goodTrack++;
+
+ }
+ } // end fGrid==1
+ } // end loop on entries (tpc tracks)
- } // end if(fOpt==1 || option=="tpc")
+// // set the signal flags
+// fSignalFlag.Set(goodTrack,sflag);
+// fCutFlag.Set(goodTrack,cflag);
- delete[] enT;
- delete[] ptT;
- delete[] etaT;
- delete[] phiT;
- delete[] thetaT;
+// delete sflag;
+// delete cflag;
+
+ // } // end loop on entries (tpc tracks)
+ if(fGrid==0) {
+ fNTracks = nTracksTpcOnly;
+ fNTracksCut = nTracksTpcOnlyCut;
+ if(fDebug>10){
+ cout << "fNTracks : " << fNTracks << endl;
+ cout << "fNTracksCut : " << fNTracksCut << endl;
+ }
+ }
+ if(fGrid==1) {
+ fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
+ fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
+ if(fDebug>10){
+ cout << "fNTracks : " << fNTracks << endl;
+ cout << "fNTracksCut : " << fNTracksCut << endl;
+ }
+ }
+
}
//__________________________________________________________
#include <TMatrixD.h>
#include <TArrayD.h>
-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)
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
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
};
void AliJetFinder::ConnectAOD(AliAODEvent* aod)
{
// Connect to the AOD
- printf("Connect AOD \n");
fAODjets = aod->GetJets();
- printf("Connect AOD %p \n", fAODjets);
}
// --- AliRoot header files ---
//#include "../EMCAL/AliJetGeometry.h"
#include "AliJetDummyGeo.h"
-//#include "AliEMCALGeometry.h"
#include "AliJetHadronCorrectionv1.h"
ClassImp(AliJetHadronCorrectionv1)
#define HCPARAMETERS 6
#define HCPARAMETERSETS 2
-class AliEMCALGeometry;
class AliJetDummyGeo;
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TSystem.h>
+#include <TChain.h>
// From AliRoot ...
#include "AliJetMCReader.h"
#include "AliJetMCReaderHeader.h"
// root
#include <TClonesArray.h>
+#include <TRefArray.h>
+#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"
////////////////////////////////////////////////////////////////////////
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();
#include <TObject.h>
#include <TChain.h>
#include <TArrayI.h>
-#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
{
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;}
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*/) {;}
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)
};
// Author: jgcn@mda.cinvestav.mx
//-------------------------------------------------------------------------
+#include <TMath.h>
+
#include "AliJetReaderHeader.h"
ClassImp(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(""),
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(""),
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;}
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
* 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
-#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 <TObject.h>
#include "AliJetFinderTypes.h"
{
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
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)
{
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;}
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;
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
--- /dev/null
+/**************************************************************************
+ * 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 <Riostream.h>
+
+#include <TArrayF.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TRefArray.h>
+
+#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; i<nIn; i++)
+ {
+ AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(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<nj; kj++){
+ if ((etaJet[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<nJ; ijet++){
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > 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<nJ; k++){
+ Float_t detamax = etaJet[k] + rc;
+ Float_t detamin = etaJet[k] - rc;
+ Float_t accmax = 0.0; Float_t accmin = 0.0;
+ if(detamax > 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; ljet<nJ; ljet++){
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
+ }
+
+ // estimate new total background
+ Float_t areaT = 4*(header->GetLegoEtaMax())*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<nJ; ijet++){
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > 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<nJ; k++){
+ Float_t detamax = etaJet[k] + rc;
+ Float_t detamin = etaJet[k] - rc;
+ Float_t accmax = 0.0; Float_t accmin = 0.0;
+ if(detamax > 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; ljet<nJ; ljet++){
+ Float_t areaRatio = areaJet[ljet]/areaOut;
+ etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
+ }
+
+ etbgTotalN = etbgStat;
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::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)
+{
+ // Cone background subtraction method taking into acount dEt/deta distribution
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ //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; mjet<nJ; mjet++){
+ char hEtname[256]; char hAreaname[256];
+ sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+ hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
+ hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
+ }
+ // background energy and area
+ TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
+ TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+
+ //fill energies
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > 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<nJ; ijet++){ // loop for all jets
+ Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+ Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+ Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+ Float_t areaj = 0.0;
+ if(deta0 > 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(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
+ if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
+ if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
+ }
+ hAreaJet[ijet]->Fill(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; kjet<nJ; kjet++){
+ etJet[kjet] = 0.0; // first clear etJet for this jet
+ for(Int_t bin = 0; bin< ndiv; bin++){
+ if(hAreaJet[kjet]->GetBinContent(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; ljet<nJ; ljet++){ // loop for all jets
+ delete hEtJet[ljet];
+ delete hAreaJet[ljet];
+ }
+
+ delete hEtBackg;
+ delete hAreaBackg;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliUA1JetFinderV2::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)
+{
+ // Ratio background subtraction method taking into acount dEt/deta distribution
+ AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+ //factor F calc before
+ Float_t bgRatioCut = header->GetBackgCutRatio();
+
+
+ //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; mjet<nJ; mjet++){
+ char hEtname[256]; char hAreaname[256];
+ sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+ hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax); // change range
+ hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax); // change range
+ }
+ // background energy and area
+ TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range
+ TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range
+
+ //fill energies
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ //if((fReader->GetCutFlag(jpart)) != 1) continue;
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > 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<nJ; ijet++){ // loop for all jets
+ Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+ Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+ Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+ Float_t areaj = 0.0;
+ if(deta0 > 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(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
+ if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
+ if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
+ }
+ hAreaJet[ijet]->Fill(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; kjet<nJ; kjet++){
+ etJet[kjet] = 0.0; // first clear etJet for this jet
+ for(Int_t bin = 0; bin< ndiv; bin++){
+ if(hAreaJet[kjet]->GetBinContent(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; ljet<nJ; ljet++){ // loop for all jets
+ delete hEtJet[ljet];
+ delete hAreaJet[ljet];
+ }
+
+ delete hEtBackg;
+ delete hAreaBackg;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliUA1JetFinderV2::Reset()
+{
+ fLego->Reset();
+ 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();
+
+}
--- /dev/null
+#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
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),
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;}
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;}
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
#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+;
#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+;
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 \
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff