]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
More cluster iso and shape variables. (M. Cosentino)
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Sep 2011 07:16:43 +0000 (07:16 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Sep 2011 07:16:43 +0000 (07:16 +0000)
PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.cxx
PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.h

index 87b6feccd739fcb6599989099a30d0ed8f178886..1bbc53e98772e2c05d2e7af4ad410c33f98c5012 100644 (file)
@@ -1034,6 +1034,7 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
     cl->fIdMax    = id;
     cl->fEmax     = emax;
     cl->fTmax    =  cells->GetCellTime(id);
+    cl->fE2max   =  GetSecondMaxCell(clus);
     if (clus->GetDistanceToBadChannel()<10000)
       cl->fDbc    = clus->GetDistanceToBadChannel();
     if (!TMath::IsNaN(clus->GetDispersion()))
@@ -1046,6 +1047,11 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
     GetSigma(clus,maxAxis,minAxis);
     clus->SetTOF(maxAxis);     // store sigma in TOF
     cl->fSig      = maxAxis;
+    Double_t sEtaEta = 0;
+    Double_t sPhiPhi = 0;
+    GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
+    cl->fSigEtaEta = sEtaEta;
+    cl->fSigPhiPhi = sPhiPhi;
     Double_t clusterEcc = 0;
     if (maxAxis > 0)
       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
@@ -1054,8 +1060,20 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
     cl->fTrIso    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
     cl->fTrIso1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
     cl->fTrIso2   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
+    cl->fTrIsoD1    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
+    cl->fTrIso1D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
+    cl->fTrIso2D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
+    cl->fTrIsoD3    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
+    cl->fTrIso1D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
+    cl->fTrIso2D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
+    cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
     cl->fCeCore   = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
     cl->fCeIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
+    cl->fCeIso1   = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
+    cl->fCeIso3  = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
+    cl->fCeIso4x4 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
+    cl->fCeIso5x5 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
+    cl->fCeIso3x22 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
 
     if (fAmpInTrigger) { // fill trigger info if present
       Double_t trigpen = 0;
@@ -1101,11 +1119,13 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
       }
 
       Float_t tmpR=-1, tmpZ=-1;
+      Double_t dedx = 0;
       if (!fDoTrMatGeom) {
         AliExternalTrackParam *tParam = 0;
         if (fEsdEv) {
           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
           tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
+         dedx = esdTrack->GetTPCsignal();
         } else 
           tParam = new AliExternalTrackParam(track);
 
@@ -1141,6 +1161,7 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
         cl->fTrDz   = tmpZ;
         cl->fTrDr   = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
         cl->fTrEp   = clus->E()/track->P();
+       cl->fTrDedx = dedx;
         cl->fIsTrackM = 1;
       }
     }
@@ -1978,6 +1999,35 @@ Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t c
   }
   return cellIsolation;
 }
+//________________________________________________________________________
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
+{
+  // Compute isolation based on cell content, in a NxM rectangle.
+
+  AliVCaloCells *cells = fEsdCells;
+  if (!cells)
+    cells = fAodCells; 
+  if (!cells)
+    return 0;
+
+  Double_t cellIsolation = 0;
+  Int_t ncells = cells->GetNumberOfCells();
+  for (Int_t i = 0; i<ncells; ++i) {
+    Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
+    Float_t eta=-1, phi=-1;
+    fGeom->EtaPhiFromIndex(absID,eta,phi);
+    Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
+    Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
+    if(TMath::Abs(etadiff)/0.014>N)
+      continue;
+    if(TMath::Abs(phidiff)/0.014>M)
+      continue;
+    Double_t cellE = cells->GetAmplitude(i);
+    cellIsolation += cellE;
+  }
+  return cellIsolation;
+}
+
 
 //________________________________________________________________________
 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
@@ -2080,6 +2130,63 @@ void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigma
   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin)); 
 }
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
+{
+  // Calculate the (E) weighted variance along the pseudorapidity.
+  
+  sEtaEta = 0;
+  sPhiPhi = 0;
+
+  Double_t Ec  = c->E(); // cluster energy
+  if(Ec<=0)
+    return;
+
+  const Int_t ncells = c->GetNCells();
+
+  Double_t EtaC    = 0;  // cluster first moment along eta
+  Double_t PhiC    = 0;  // cluster first moment along phi
+  Double_t Setaeta = 0;  // cluster second central moment along eta
+  Double_t Sphiphi = 0;  // cluster second central moment along phi
+  Double_t w[ncells];    // weight max(0,4.5*log(E_i/Ec))
+  Double_t sumw = 0;
+  Int_t id[ncells];
+
+  AliVCaloCells *cells = fEsdCells;
+  if (!cells)
+    cells = fAodCells;
+
+  if (!cells)
+    return;
+
+  if (ncells==1)
+    return;
+
+  for(Int_t j=0; j<ncells; ++j) {
+    id[j] = TMath::Abs(c->GetCellAbsId(j));
+    Double_t cellen = cells->GetCellAmplitude(id[j]);
+    w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
+    TVector3 pos;
+    fGeom->GetGlobal(id[j],pos);
+    EtaC += w[j]*pos.Eta();
+    PhiC += w[j]*pos.Phi();
+    sumw += w[j];
+  }
+  EtaC /= sumw;
+  PhiC /= sumw;
+  
+  for(Int_t j=0; j<ncells; ++j) {
+    TVector3 pos;
+    fGeom->GetGlobal(id[j],pos);
+    Setaeta =  w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
+    Sphiphi =  w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
+  }
+  Setaeta /= sumw;
+  sEtaEta = TMath::Sqrt(Setaeta);
+  Sphiphi /= sumw;
+  sPhiPhi = TMath::Sqrt(Sphiphi);
+}
+
 
 //________________________________________________________________________
 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
@@ -2127,6 +2234,30 @@ Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t
   } 
   return trkIsolation;
 }
+//________________________________________________________________________
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
+{
+  // Compute isolation based on tracks.
+  
+  Double_t trkIsolation = 0;
+  Int_t ntrks = fSelPrimTracks->GetEntries();
+  for(Int_t j = 0; j<ntrks; ++j) {
+    AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
+    if (!track)
+      continue;
+    if (track->Pt()<pt)
+      continue;
+    Float_t eta = track->Eta();
+    Float_t phi = track->Phi();
+    Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
+    Double_t etadiff = (eta-cEta);
+    if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
+      continue;
+    trkIsolation += track->Pt();
+  } 
+  return trkIsolation;
+}
+
 
 //________________________________________________________________________
 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
@@ -2323,3 +2454,31 @@ void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index
     ProcessDaughters(dmc,i,arr);
   }
 }
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCell(AliVCluster *clus)
+{
+  // Get second maximum cell.
+
+  AliVCaloCells *cells = fEsdCells;
+  if (!cells)
+    cells = fAodCells;
+  if (!cells)
+    return -1;
+  Double_t secondEmax=0, firstEmax=0;
+  Double_t cellen;
+  for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
+    Int_t absId = clus->GetCellAbsId(iCell);
+    cellen = cells->GetCellAmplitude(absId);
+    if(cellen > firstEmax)
+      firstEmax = cellen;
+  }
+  for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
+    Int_t absId = clus->GetCellAbsId(iCell);
+    cellen = cells->GetCellAmplitude(absId);
+    if(cellen < firstEmax && cellen > secondEmax)
+      secondEmax = cellen;
+  }
+  return secondEmax;
+}
index ec2398790c0713c2af8015ad0a05cf820cd10c61..525a2fdbd6f20b9f02c5ec3ccd15c8cdac21bcde 100644 (file)
@@ -84,12 +84,15 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   void         FillVertex(AliStaVertex *v, const AliESDVertex *esdv);
   void         FillVertex(AliStaVertex *v, const AliAODVertex *aodv);
   Double_t     GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2)                        const;
+  Double_t     GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M)                              const;
   Double_t     GetCellEnergy(const AliVCluster *c)    const;
   Double_t     GetMaxCellEnergy(const AliVCluster *c) const { Short_t id=-1; return GetMaxCellEnergy(c,id); }
   Double_t     GetMaxCellEnergy(const AliVCluster *c, Short_t &id)                                        const;
   Int_t        GetNCells(const AliVCluster *c, Double_t emin=0.)                                          const;
   void         GetSigma(const AliVCluster *c, Double_t &sigmaMax, Double_t &sigmaMin)                     const;
+  void         GetSigmaEtaEta(const AliVCluster *c, Double_t &sigmaEtaEta, Double_t &sigmaPhiPhi)         const;
   Double_t     GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.)       const;
+  Double_t     GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta=0.015, Double_t dPhi=0.3, Double_t pt=0.)       const;
   Double_t     GetTrigEnergy(const AliVCluster *c)                                                        const;
   Bool_t       IsShared(const AliVCluster *c)                                                             const;
   void         PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level=0)                 const;
@@ -97,6 +100,7 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   void         PrintTrackRefs(AliMCParticle *p)                                                           const;
   void         ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr);
   void         ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr);
+  Double_t     GetSecondMaxCell(AliVCluster *clus);
 
     // input members
   TString                fCentVar;                // variable for centrality determination
@@ -300,9 +304,11 @@ class AliStaCluster : public TObject
 {
  public:
     AliStaCluster() : TObject(), fE(0), fR(0), fEta(0), fPhi(0), fN(0), fN1(0), fN3(0), fIdMax(0), fEmax(0), fTmax(0), 
-                      fDbc(-1), fDisp(-1), fM20(0), fM02(0), fEcc(0), fSig(0), fIsTrackM(0), fTrDz(0), fTrDr(-1), 
-                      fTrEp(0), fTrIso(0), fTrIso1(0), fTrIso2(0), fCeIso(0), fCeCore(0), fIsTrigM(0), fTrigE(-1), 
-                      fTrigMaskE(-1), fIsShared(0), fMcLabel(-1), fEmbE(0) {;}
+                      fE2max(0),fDbc(-1), fDisp(-1), fM20(0), fM02(0), fEcc(0), fSig(0), fSigEtaEta(0), fSigPhiPhi(0),
+                      fIsTrackM(0), fTrDz(0), fTrDr(-1), fTrEp(0), fTrDedx(0), fTrIso(0), fTrIso1(0), fTrIso2(0),  
+                      fTrIsoD1(0), fTrIso1D1(0), fTrIso2D1(0), fTrIsoD3(0), fTrIso1D3(0), fTrIso2D3(0),fTrIsoStrip(0),
+                      fCeIso(0), fCeIso1(0), fCeIso3(0), fCeIso4x4(0), fCeIso5x5(0), fCeCore(0), fCeIso3x22(0), 
+                      fIsTrigM(0), fTrigE(-1), fTrigMaskE(-1), fIsShared(0), fMcLabel(-1), fEmbE(0) {;}
 
  public:
   Double32_t    fE;                //[0,0,16] energy
@@ -315,21 +321,37 @@ class AliStaCluster : public TObject
   UShort_t      fIdMax;            //         id maximum cell
   Double32_t    fEmax;             //[0,0,16] energy of maximum cell
   Double32_t    fTmax;             //[0,0,16] time of maximum cell
+  Double32_t    fE2max;            //[0,0,16] energy of second maximum cell
   Double32_t    fDbc;              //[0,0,16] distance to nearest bad channel
   Double32_t    fDisp;             //[0,0,16] cluster dispersion, for shape analysis
   Double32_t    fM20;              //[0,0,16] 2-nd moment along the main eigen axis
   Double32_t    fM02;              //[0,0,16] 2-nd moment along the second eigen axis
   Double32_t    fEcc;              //[0,0,16] eccentricity
   Double32_t    fSig;              //[0,0,16] sigma
+  Double32_t    fSigEtaEta;        //[0,0,16] sigma eta-eta
+  Double32_t    fSigPhiPhi;        //[0,0,16] sigma phi-phi
   Bool_t        fIsTrackM;         //         if true then track values are set
   Double32_t    fTrDz;             //[0,0,16] dZ to nearest track
   Double32_t    fTrDr;             //[0,0,16] dR to nearest track (in x,y)
   Double32_t    fTrEp;             //[0,0,16] E/P to nearest track 
+  Double32_t    fTrDedx;           //[0,0,16] dE/dx (TPC signal) to nearest track 
   Double32_t    fTrIso;            //[0,0,16] track isolation
   Double32_t    fTrIso1;           //[0,0,16] track isolation (pt>1GeV/c)
   Double32_t    fTrIso2;           //[0,0,16] track isolation (pt>2GeV/c)
-  Double32_t    fCeIso;            //[0,0,16] cell isolation
-  Double32_t    fCeCore;           //[0,0,16] cell content in R=0.025
+  Double32_t    fTrIsoD1;          //[0,0,16] track isolation, iso dist 0.25
+  Double32_t    fTrIso1D1;         //[0,0,16] track isolation (pt>1GeV/c), iso dist 0.1
+  Double32_t    fTrIso2D1;         //[0,0,16] track isolation (pt>2GeV/c), iso dist 0.1
+  Double32_t    fTrIsoD3;          //[0,0,16] track isolation, iso dist 0.3
+  Double32_t    fTrIso1D3;         //[0,0,16] track isolation (pt>1GeV/c), iso dist 0.3
+  Double32_t    fTrIso2D3;         //[0,0,16] track isolation (pt>2GeV/c), iso dist 0.3
+  Double32_t    fTrIsoStrip;       //[0,0,16] track isolation strip, dEtaXdPhi=0.015x0.3
+  Double32_t    fCeIso;            //[0,0,16] cell isolation in R=0.20
+  Double32_t    fCeIso1;           //[0,0,16] cell isolation in R=0.10
+  Double32_t    fCeIso3;          //[0,0,16] cell isolation in  R=0.30
+  Double32_t    fCeIso4x4;         //[0,0,16] cell isolation in  4x4 cells
+  Double32_t    fCeIso5x5;         //[0,0,16] cell isolation in  5x5 cells
+  Double32_t    fCeCore;           //[0,0,16] cell content in R=0.05 
+  Double32_t    fCeIso3x22;        //[0,0,16] cell isolation in rectangular strip of dEtaXdPhi=0.042x0.308
   Bool_t        fIsTrigM;          //         if true then trigger values are set
   Double32_t    fTrigE;            //[0,0,16] trigger tower energy
   Double32_t    fTrigMaskE;        //[0,0,16] masked trigger tower energy
@@ -360,7 +382,7 @@ class AliStaPart : public TObject
 {
  public:
   AliStaPart() : TObject(), fPt(0), fEta(0), fPhi(0), fVR(0), fVEta(0), fVPhi(0), fPid(0), fMo(-1), fDet(-2), 
-                 fLab(-1), fNs(0) { memset(fDs,-1,sizeof(Short_t)*9); }
+                 fLab(-1), fNs(0) { memset(fDs,-1,sizeof(Short_t)*99); }
 
   Int_t         OnEmcal() const { return (fDet==8);  }
   Int_t         IsSim()   const { return (fDet!=-2); }
@@ -378,7 +400,7 @@ class AliStaPart : public TObject
     // the following must be filled before first usage
   Short_t       fLab;              //!        label (index in array)
   Short_t       fNs;               //!        number of daughters
-  Short_t       fDs[9];            //!        daughters
+  Short_t       fDs[99];           //!        daughters
 
   ClassDef(AliStaPart,1) // Particle class
 };