New track selection and new track matching
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Apr 2011 13:10:35 +0000 (13:10 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Apr 2011 13:10:35 +0000 (13:10 +0000)
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.h

index c1fb152..06e0f61 100644 (file)
@@ -27,6 +27,7 @@
 #include "AliESDVertex.h"
 #include "AliESDtrack.h"
 #include "AliESDtrackCuts.h"
+#include "AliEventplane.h"
 #include "AliGeomManager.h"
 #include "AliInputEventHandler.h"
 #include "AliLog.h"
@@ -53,13 +54,12 @@ AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
     fMinErat(0),
     fMinEcc(0),
     fGeoName("EMCAL_FIRSTYEARV1"),
-    fMinNClustPerTrack(50),
-    fMinPtPerTrack(1.0), 
+    fMinNClusPerTr(50),
     fIsoDist(0.2),
     fTrClassNames(""),
     fTrCuts(0),
-    fDoTrackMatWithGeom(0),
-    fDoConstrain(0),
+    fPrimTrCuts(0),
+    fDoTrMatGeom(0),
     fIsGeoMatsSet(0),
     fNEvs(0),
     fGeom(0),
@@ -75,6 +75,7 @@ AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
     fAodCells(0),
     fPtRanges(0),
     fSelTracks(0),
+    fSelPrimTracks(0),
     fNtuple(0),
     fHeader(0),
     fPrimVert(0),
@@ -106,6 +107,9 @@ AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
     fHClustEnergySigma(0x0),
     fHClustSigmaSigma(0x0),
     fHClustNCellEnergyRatio(0x0),
+    fHMatchDr(0x0),
+    fHMatchDz(0x0),
+    fHMatchEp(0x0),
     fHPionEtaPhi(0x0),
     fHPionMggPt(0x0),
     fHPionMggAsym(0x0),
@@ -135,6 +139,7 @@ AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
   delete fReco; fReco = 0;
   delete fTrClassNamesArr;
   delete fSelTracks;
+  delete fSelPrimTracks;
   delete [] fHColuRow;
   delete [] fHColuRowE;
   delete [] fHCellMult;
@@ -154,6 +159,7 @@ void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
   fOutput = new TList();
   fOutput->SetOwner();
   fSelTracks = new TObjArray;
+  fSelPrimTracks = new TObjArray;
 
   if (fDoNtuple) {
     TFile *f = OpenFile(1);
@@ -179,8 +185,8 @@ void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
   }
 
   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
-  Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
-  Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
+  Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
+  Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
 
   // histograms
   TH1::SetDefaultSumw2(kTRUE);
@@ -306,6 +312,13 @@ void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
   fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
   fOutput->Add(fHClustNCellEnergyRatio);
 
+  // histograms for track matching
+  fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
+  fOutput->Add(fHMatchDr);
+  fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
+  fOutput->Add(fHMatchDz);
+  fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
+  fOutput->Add(fHMatchEp);
   // histograms for pion candidates
   fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
   fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
@@ -363,7 +376,7 @@ void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
     am->LoadBranch("header");
   }
 
-  if (fDoTrackMatWithGeom && !AliGeomManager::GetGeometry()) { // get geometry 
+  if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry 
     AliWarning("Accessing geometry from OCDB, this is not very efficient!");
     AliCDBManager *cdb = AliCDBManager::Instance();
     if (!cdb->IsDefaultStorageSet())
@@ -488,9 +501,9 @@ void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
   }
 
   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
-  fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
+  fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
   fEsdCells    = 0; // will be set if ESD input used
-  fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
+  fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
   fAodCells    = 0; // will be set if AOD input used
 
@@ -565,6 +578,7 @@ void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
     ClusterAfterburner();
 
   CalcTracks();
+  CalcPrimTracks();
   CalcClusterProps();
 
   FillCellHists();
@@ -593,7 +607,7 @@ void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 //________________________________________________________________________
 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
 {
-  // Calculate track properties.
+  // Calculate track properties (including secondaries).
 
   fSelTracks->Clear();
 
@@ -612,13 +626,74 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
   if (!tracks)
     return;
 
+  if (fEsdEv) {
+    const Int_t Ntracks = tracks->GetEntries();
+    for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
+      AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
+      if (!track) {
+        AliWarning(Form("Could not receive track %d\n", iTracks));
+        continue;
+      }
+      if (fTrCuts && !fTrCuts->IsSelected(track))
+        continue;
+      Double_t eta = track->Eta();
+      if (eta<-1||eta>1)
+        continue;
+      fSelTracks->Add(track);
+    }
+  } else {
+    Int_t ntracks = tracks->GetEntries();
+    for (Int_t i=0; i<ntracks; ++i) {
+      AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
+      if (!track)
+        continue;
+      Double_t eta = track->Eta();
+      if (eta<-1||eta>1)
+        continue;
+      if(track->GetTPCNcls()<fMinNClusPerTr)
+        continue;
+
+      if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL 
+        AliExternalTrackParam tParam(track);
+        if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
+          Double_t trkPos[3];
+          tParam.GetXYZ(trkPos);
+          track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
+        }
+      }
+      fSelTracks->Add(track);
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
+{
+  // Calculate track properties.
+
+  fSelPrimTracks->Clear();
+
+  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+  TClonesArray *tracks = 0;
+  if (fEsdEv) {
+    am->LoadBranch("Tracks");
+    TList *l = fEsdEv->GetList();
+    tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
+  } else {
+    am->LoadBranch("tracks");
+    TList *l = fAodEv->GetList();
+    tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
+  }
+
+  if (!tracks)
+    return;
+
   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
-  Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
-  Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
+  Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
+  Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
 
   if (fEsdEv) {
-    if (fDoConstrain)
-      fSelTracks->SetOwner(kTRUE);
+    fSelPrimTracks->SetOwner(kTRUE);
     am->LoadBranch("PrimaryVertex.");
     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
     am->LoadBranch("SPDVertex.");
@@ -638,13 +713,6 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
         continue;
       if (track->Phi()<phimin||track->Phi()>phimax)
         continue;
-      if(track->GetTPCNcls()<fMinNClustPerTrack)
-        continue;
-
-      if (!fDoConstrain) {
-        fSelTracks->Add(track);
-        continue;
-      }
 
       AliESDtrack copyt(*track);
       Double_t bfield[3];
@@ -685,7 +753,7 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
         aTrack->SetChi2perNDF(-1);
       aTrack->SetFlags(copyt.GetStatus());
       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
-      fSelTracks->Add(aTrack);
+      fSelPrimTracks->Add(aTrack);
     }
   } else {
     Int_t ntracks = tracks->GetEntries();
@@ -698,18 +766,10 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
         continue;
       if (track->Phi()<phimin||track->Phi()>phimax)
         continue;
-      if(track->GetTPCNcls()<fMinNClustPerTrack)
+      if(track->GetTPCNcls()<fMinNClusPerTr)
         continue;
-
-      if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL 
-        AliExternalTrackParam tParam(track);
-        if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
-          Double_t trkPos[3];
-          tParam.GetXYZ(trkPos);
-          track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
-        }
-      }
-      fSelTracks->Add(track);
+      //todo: Learn how to set/filter AODs for prim/sec tracks
+      fSelPrimTracks->Add(track);
     }
   }
 }
@@ -728,6 +788,9 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
   Int_t nclus = clusters->GetEntries();
   Int_t ntrks = fSelTracks->GetEntries();
 
+  Bool_t btracks[6][ntrks];
+  memset(btracks,0,sizeof(btracks));
+
   for(Int_t i = 0; i<nclus; ++i) {
     fClusProps[i].Reset();
 
@@ -739,7 +802,7 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
     if (clus->E()<fMinE)
       continue;
 
-    Float_t  clsPos[3] = {0};
+    Float_t clsPos[3] = {0};
     clus->GetPosition(clsPos);
     TVector3 clsVec(clsPos);
     Double_t vertex[3] = {0};
@@ -753,15 +816,18 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
       AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
       if (!track)
         continue;
-      Double_t pt  = track->Pt();
-      if (pt<fMinPtPerTrack) 
-        continue;
+
       if (TMath::Abs(clsEta-track->Eta())>0.5)
         continue;
 
+      TVector3 vec(clsPos);
+      Int_t index =  (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
+      if (btracks[index-4][j]) {
+        continue;
+      }
+
       Float_t tmpR=-1, tmpZ=-1;
-      if (!fDoTrackMatWithGeom) {
-        
+      if (!fDoTrMatGeom) {
         AliExternalTrackParam *tParam = 0;
         if (fEsdEv) {
           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
@@ -769,18 +835,18 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
         } else 
           tParam = new AliExternalTrackParam(track);
 
-        Double_t bfield[3];
+        Double_t bfield[3] = {0};
         track->GetBxByBz(bfield);
-        TVector3 vec(clsPos);
-        Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
-        vec.RotateZ(-alpha);  //Rotate the cluster to the local extrapolation coordinate system
+        Double_t alpha = (index+0.5)*20*TMath::DegToRad();
+        vec.RotateZ(-alpha);   //Rotate the cluster to the local extrapolation coordinate system
         tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
         Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
         if (!ret) {
+          btracks[index-4][j]=1;
           delete tParam;
           continue;
         }
-        Double_t trkPos[3];
+        Double_t trkPos[3] = {0};
         tParam->GetXYZ(trkPos); //Get the extrapolated global position
         tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
         tmpZ = clsPos[2]-trkPos[2];
@@ -801,16 +867,26 @@ void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
         fClusProps[i].fTrDr    = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
         fClusProps[i].fTrDist  = d2;
         fClusProps[i].fTrEp    = clus->E()/track->P();
+        fClusProps[i].fPhiInd  = index;
       }
     }
-
-    if (0 && (fClusProps[i].fTrIndex>=0)) {
-      cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
+    
+    if (fClusProps[i].fTrIndex>=0) {
+      Int_t tid = fClusProps[i].fTrIndex;
+      if ( (btracks[0][tid] && fClusProps[i].fPhiInd==0) ||
+           (btracks[1][tid] && fClusProps[i].fPhiInd==1) ) {
+        cout << i << " " << tid << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
+        cout << btracks[0][tid] << " " << btracks[1][tid] << endl;
+      }
+      fHMatchDr->Fill(fClusProps[i].fTrDr);
+      fHMatchDz->Fill(fClusProps[i].fTrDz);
+      fHMatchEp->Fill(fClusProps[i].fTrEp);
     }
 
-    fClusProps[i].fTrIso      = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
-    fClusProps[i].fTrLowPtIso = 0;
-    fClusProps[i].fCellIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
+    fClusProps[i].fTrIso   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
+    fClusProps[i].fTrIso1  = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
+    fClusProps[i].fTrIso2  = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
+    fClusProps[i].fCellIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
   }
 }
 
@@ -1037,7 +1113,18 @@ void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
   fHeader->fCqual     = cent->GetQuality();
+  
+  AliEventplane *ep = InputEvent()->GetEventplane();
+  if (ep) {
+    if (ep->GetQVector())
+      fHeader->fPsi     = ep->GetQVector()->Phi()/2. ;
+    fHeader->fPsi = -1;
+    if (ep->GetQsub1()&&ep->GetQsub2())
+      fHeader->fPsiRes  = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
+    else 
+      fHeader->fPsiRes = 0;
+  }
+
   Double_t val = 0;
   TString trgclasses(fHeader->fFiredTriggers);
   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
@@ -1090,7 +1177,7 @@ void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
     cl->fR        = vpos.Perp();
     cl->fEta      = vpos.Eta();
     cl->fPhi      = vpos.Phi();
-    cl->fN        =  clus->GetNCells();
+    cl->fN        = clus->GetNCells();
     cl->fN1       = GetNCells(clus,0.100);
     cl->fN3       = GetNCells(clus,0.300);
     Short_t id    = -1;
@@ -1107,6 +1194,8 @@ void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
     cl->fTrDr     = fClusProps[i].fTrDr;;
     cl->fTrEp     = fClusProps[i].fTrEp;;
     cl->fTrIso    = fClusProps[i].fTrIso;
+    cl->fTrIso1   = fClusProps[i].fTrIso1;
+    cl->fTrIso2   = fClusProps[i].fTrIso2;
     cl->fCeIso    = fClusProps[i].fCellIso;
   }
   fNtuple->Fill();
@@ -1343,17 +1432,19 @@ Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) cons
 }
 
 //________________________________________________________________________
-Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
 {
   // Compute isolation based on tracks.
   
   Double_t trkIsolation = 0;
   Double_t rad2 = radius*radius;
-  Int_t ntrks = fSelTracks->GetEntries();
+  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_0_2pi(phi-cPhi);
index 78ad716..23d4e0d 100644 (file)
@@ -36,30 +36,30 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   void         UserExec(Option_t *option);
   void         Terminate(Option_t *);
 
-  void         SetAsymMax(Double_t asymMax)                   { fAsymMax = asymMax;         }
-  void         SetCentrality(const char *name)                { fCentVar = name;            }
+  void         SetAsymMax(Double_t asymMax)                   { fAsymMax       = asymMax;   }
+  void         SetCentrality(const char *n)                   { fCentVar       = n;         }
   void         SetCentralityRange(Double_t from, Double_t to) { fCentFrom=from; fCentTo=to; }
-  void         SetClusName(const char *name)                  { fClusName = name;           }
+  void         SetClusName(const char *n)                     { fClusName      = n;         }
   void         SetDoAfterburner(Bool_t b)                     { fDoAfterburner = b;         }
-  void         SetDoTrackMatWithGeom(Bool_t b)                { fDoTrackMatWithGeom = b;    }
-  void         SetDoTrackVtxConstrain(Bool_t b)               { fDoConstrain = b;           }
-  void         SetFillNtuple(Bool_t b)                        { fDoNtuple = b;              }
-  void         SetGeoName(const char *n)                      { fGeoName = n;               }
-  void         SetIsoDist(Double_t d)                         { fIsoDist = d;               }
-  void         SetMinClusEnergy(Double_t e)                   { fMinE = e;                  }
-  void         SetMinEcc(Double_t ecc)                        { fMinEcc = ecc;              }
-  void         SetMinErat(Double_t erat)                      { fMinErat = erat;            }
-  void         SetMinNClustersPerTrack(Double_t mct)          { fMinNClustPerTrack = mct;   }
-  void         SetMinPtPerMatchedTrack(Double_t mpt)          { fMinPtPerTrack = mpt;       }
-  void         SetNminCells(Int_t n)                          { fNminCells = n;             }
-  void         SetTrClassNames(const char *n)                 { fTrClassNames = n;          }
-  void         SetTrackCuts(AliESDtrackCuts *c)               { fTrCuts = c;                }
-  void         SetUseQualFlag(Bool_t b)                       { fUseQualFlag = b;           }
+  void         SetDoTrackMatWithGeom(Bool_t b)                { fDoTrMatGeom   = b;         }
+  void         SetFillNtuple(Bool_t b)                        { fDoNtuple      = b;         }
+  void         SetGeoName(const char *n)                      { fGeoName       = n;         }
+  void         SetIsoDist(Double_t d)                         { fIsoDist       = d;         }
+ void          SetMinNClustersPerTrack(Double_t m)            { fMinNClusPerTr = m;         }
+  void         SetMinClusEnergy(Double_t e)                   { fMinE          = e;         }
+  void         SetMinEcc(Double_t ecc)                        { fMinEcc        = ecc;       }
+  void         SetMinErat(Double_t erat)                      { fMinErat       = erat;      }
+  void         SetNminCells(Int_t n)                          { fNminCells     = n;         }
+  void         SetTrClassNames(const char *n)                 { fTrClassNames  = n;         }
+  void         SetTrackCuts(AliESDtrackCuts *c)               { fTrCuts        = c;         }
+  void         SetPrimTrackCuts(AliESDtrackCuts *c)           { fPrimTrCuts    = c;         }
+  void         SetUseQualFlag(Bool_t b)                       { fUseQualFlag   = b;         }
   void         SetVertexRange(Double_t z1, Double_t z2)       { fVtxZMin=z1; fVtxZMax=z2;   }
 
  protected:
   virtual void CalcClusterProps();
   virtual void CalcTracks();
+  virtual void CalcPrimTracks();
   virtual void ClusterAfterburner();
   virtual void FillCellHists();
   virtual void FillClusHists();
@@ -69,26 +69,28 @@ 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     GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2)                  const;
   Double_t     GetMaxCellEnergy(AliVCluster *c) const { Short_t id=-1; return GetMaxCellEnergy(c,id); }
-  Double_t     GetMaxCellEnergy(AliVCluster *c, Short_t &id)                                      const;
-  Int_t        GetNCells(AliVCluster *c, Double_t emin=0.)                                        const;
-  void         GetSigma(AliVCluster *c, Double_t &sigmaMax, Double_t &sigmaMin)                   const;
-  Double_t     GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2)               const;
+  Double_t     GetMaxCellEnergy(AliVCluster *c, Short_t &id)                                        const;
+  Int_t        GetNCells(AliVCluster *c, Double_t emin=0.)                                          const;
+  void         GetSigma(AliVCluster *c, Double_t &sigmaMax, Double_t &sigmaMin)                     const;
+  Double_t     GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.) const;
 
   class ClusProps {
     public:
       ClusProps() : fTrIndex(-1), fTrDz(-1), fTrDr(-1), fTrDist(-1), fTrEp(0), 
-                    fTrIso(0), fTrLowPtIso(0), fCellIso(0) {}
-      void Reset() { fTrIndex=-1; fTrDz=-1; fTrDr=-1; fTrDist=-1; fTrEp=0; fTrIso=0; fTrLowPtIso=0; fCellIso=0; }
+                    fTrIso(0), fTrIso1(0), fTrIso2(0), fCellIso(0), fPhiInd(-1) {}
+      void Reset() { fTrIndex=-1; fTrDz=-1; fTrDr=-1; fTrDist=-1; fTrEp=0; fTrIso=0; fTrIso1=0; fTrIso2=0; fCellIso=0; fPhiInd=-1; }
       Int_t    fTrIndex;
       Double_t fTrDz;
       Double_t fTrDr;
       Double_t fTrDist;
       Double_t fTrEp;
       Double_t fTrIso;
-      Double_t fTrLowPtIso;
+      Double_t fTrIso1;
+      Double_t fTrIso2;
       Double_t fCellIso;
+      Short_t  fPhiInd;
   };
     // input members
   TString                fCentVar;                // variable for centrality determination
@@ -106,13 +108,12 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   Double_t               fMinErat;                // minimum emax/ec ratio (def=0)
   Double_t               fMinEcc;                 // minimum eccentricity (def=0)
   TString                fGeoName;                // geometry name (def = EMCAL_FIRSTYEARV1)
-  Double_t               fMinNClustPerTrack;      // minimum number of cluster per track (def=50)
-  Double_t               fMinPtPerTrack;          // minimum pT per track (def=0.25 GeV/c)
+  Double_t               fMinNClusPerTr;          // minimum number of cluster per track (def=50)
   Double_t               fIsoDist;                // isolation distance (def=0.2)
   TString                fTrClassNames;           // trigger class names
   AliESDtrackCuts       *fTrCuts;                 // track cuts
-  Bool_t                 fDoTrackMatWithGeom;     // track matching including geometry
-  Bool_t                 fDoConstrain;            // if true constrain tracks to vertex 
+  AliESDtrackCuts       *fPrimTrCuts;                 // track cuts
+  Bool_t                 fDoTrMatGeom;            // track matching including geometry
 
     // derived members (ie with ! after //)
   Bool_t                 fIsGeoMatsSet;           //!indicate that geo matrices are set 
@@ -130,6 +131,7 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   AliAODCaloCells       *fAodCells;               //!pointer to aod cells
   TAxis                 *fPtRanges;               //!pointer to pt ranges
   TObjArray             *fSelTracks;              //!pointer to selected tracks
+  TObjArray             *fSelPrimTracks;          //!pointer to selected tracks
   ClusProps              fClusProps[1000];        //!array of cluster properties
     // ntuple
   TTree                 *fNtuple;                 //!pointer to ntuple
@@ -167,6 +169,10 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   TH2                   *fHClustEnergySigma;      //!histo for cluster energy vs. variance over long axis 
   TH2                   *fHClustSigmaSigma;       //!histo for sigma vs. lambda_0 comparison
   TH2                   *fHClustNCellEnergyRatio; //!histo for cluster n tow vs. energy ratio
+    // histograms for track matching
+  TH1                   *fHMatchDr;               //!histo for dR track cluster matching
+  TH1                   *fHMatchDz;               //!histo for dZ track cluster matching
+  TH1                   *fHMatchEp;               //!histo for E/p track cluster matching
     // histograms for pion candidates
   TH2                   *fHPionEtaPhi;            //!histo for pion eta vs. phi
   TH2                   *fHPionMggPt;             //!histo for pion mass vs. pT
@@ -178,7 +184,7 @@ class AliAnalysisTaskEMCALPi0PbPb : public AliAnalysisTaskSE {
   AliAnalysisTaskEMCALPi0PbPb(const AliAnalysisTaskEMCALPi0PbPb&);            // not implemented
   AliAnalysisTaskEMCALPi0PbPb &operator=(const AliAnalysisTaskEMCALPi0PbPb&); // not implemented
 
-  ClassDef(AliAnalysisTaskEMCALPi0PbPb, 5); // Analysis task for neutral pions in Pb+Pb
+  ClassDef(AliAnalysisTaskEMCALPi0PbPb, 6) // Analysis task for neutral pions in Pb+Pb
 };
 #endif
 
@@ -189,13 +195,9 @@ class AliStaHeader
  public:
   AliStaHeader() : fRun(0), fOrbit(0), fPeriod(0), fBx(0), fL0(0), fL1(0), fL2(0),
                    fTrClassMask(0), fTrCluster(0), fOffTriggers(0), fFiredTriggers(),
-                   fTcls(0), fV0Cent(0), fCl1Cent(0), fTrCent(0), fCqual(-1) {;}
+                   fTcls(0), fV0Cent(0), fCl1Cent(0), fTrCent(0), fCqual(-1),
+                   fPsi(0), fPsiRes(0) {;}
   virtual ~AliStaHeader() {;}
-  ULong64_t     GetEventId() const {
-                  return (((ULong64_t)fPeriod << 36) |
-                          ((ULong64_t)fOrbit  << 12) |
-                          (ULong64_t)fBx); 
-                }
 
  public:
   Int_t         fRun;            //         run number
@@ -214,8 +216,10 @@ class AliStaHeader
   Double32_t    fCl1Cent;        //[0,0,16] cl1 cent
   Double32_t    fTrCent;         //[0,0,16] tr cent
   Int_t         fCqual;          //         centrality quality
+  Double32_t    fPsi;            //[0,0,16] event-plane angle
+  Double32_t    fPsiRes;         //[0,0,16] event-plane ange resolution
 
-  ClassDef(AliStaHeader,1) // Header class
+  ClassDef(AliStaHeader,2) // Header class
 };
 
 class AliStaVertex
@@ -243,11 +247,9 @@ class AliStaVertex
 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),  
-                      fDbc(0), fDisp(0), fM20(0), fM02(0), fEcc(0), fSig(0), fTrDz(0), fTrDr(-1), fTrEp(0), 
-                      fTrIso(0), fCeIso(0) {;}
-
-//  void          GetMom(TLorentzVector& p, Double_t *vertex=0);
+  AliStaCluster() : TObject(), fE(0), fR(0), fEta(0), fPhi(0), fN(0), fN1(0), fN3(0), fIdMax(0), fEmax(0),  
+                    fDbc(0), fDisp(0), fM20(0), fM02(0), fEcc(0), fSig(0), fTrDz(0), fTrDr(-1), fTrEp(0), 
+                    fTrIso(0), fTrIso1(0), fTrIso2(0), fCeIso(0) {;}
 
  public:
   Double32_t    fE;                //[0,0,16] energy
@@ -269,8 +271,10 @@ class AliStaCluster : public TObject
   Double32_t    fTrDr;             //[0,0,16] dR to nearest track (in x,y; if neg then no match)
   Double32_t    fTrEp;             //[0,0,16] E/P 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
 
-  ClassDef(AliStaCluster,1) // Cluster class
+  ClassDef(AliStaCluster,2) // Cluster class
 };
 #endif