Now, the estimator can work with the tracklets too (R.Shahoyan)
authorbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 1 Dec 2009 13:43:51 +0000 (13:43 +0000)
committerbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 1 Dec 2009 13:43:51 +0000 (13:43 +0000)
PWG1/AliAnalysisTaskIPInfo.cxx
PWG1/AliAnalysisTaskIPInfo.h
PWG1/AliIntSpotEstimator.cxx
PWG1/AliIntSpotEstimator.h
PWG1/macros/AddTaskIntSpotESD.C
PWG1/macros/RunIPTask.C

index e1ed42a..edde240 100644 (file)
@@ -30,6 +30,7 @@
 #include <TObjArray.h>
 #include <TH1F.h>
 #include <TH2F.h>  
+#include <TNtuple.h>  
 #include <TCanvas.h>
 
 #include "AliAnalysisTask.h"
 #include "AliExternalTrackParam.h"
 #include "AliESDVertex.h"
 #include "AliESDEvent.h"
+#include "AliESDfriend.h"
 #include "AliVertexerTracks.h"
 #include "AliESDInputHandler.h"
 #include "AliAnalysisTaskIPInfo.h"
 #include "AliIntSpotEstimator.h"
+#include "AliMultiplicity.h"
 
 
 ClassImp(AliAnalysisTaskIPInfo)
 
+const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};
+
 //________________________________________________________________________
 AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name) 
 : AliAnalysisTask(name, "IP analysis"),
-  fESD(0),fOutput(0),fTracks(50)
+  fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
 {
 
   // Define input and output slots here
@@ -59,10 +64,8 @@ AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name)
   // Output slot #0 writes into a TList container
   DefineOutput(0, TList::Class());  //My private output
   //
-  for (int i=0;i<kNEst;i++) { // default settings
-    fIPEst[i] = 0;
-    SetOptions(i);
-  }
+  for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
+  //
 }
 
 //________________________________________________________________________
@@ -74,17 +77,27 @@ AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
     fOutput = 0;
   }
 }
+//________________________________________________________________________
+void AliAnalysisTaskIPInfo::SetIPCenIni(Int_t estID, Double_t x,Double_t y,Double_t z)
+{
+  // set initial estimate of the IP center
+  if (estID<0 || estID>= kNEst) return;
+  fIPCenIni[estID][0] = x;
+  fIPCenIni[estID][1] = y;
+  fIPCenIni[estID][2] = z;
+}
 
 //________________________________________________________________________
 void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
-                                      Double_t outcut,Int_t nPhiBins,Int_t nestb,
+                                      Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
                                       Double_t estmin,Double_t estmax,
                                       Int_t ntrBins,Int_t ntMn,Int_t ntMx,
-                                      Int_t nPBins,Double_t pmn,Double_t pmx)
+                                      Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
 {
   // set options for estimators
   if (estID<0 || estID>= kNEst) return;
-  fRecoVtx[estID] = recoVtx;
+  fNTrMinIP[estID] = ntrIP;
+  fRecoVtx[estID]  = recoVtx;
   fNPhiBins[estID] = nPhiBins;
   fNEstb[estID]    = nestb;
   fNTrBins[estID]  = ntrBins;
@@ -96,6 +109,7 @@ void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
   fEstMax[estID]   = estmax;
   fPMin[estID]     = pmn;
   fPMax[estID]     = pmx;
+  fFillNt[estID]   = fillNt;
 }
 
 //________________________________________________________________________
@@ -104,11 +118,13 @@ void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *)
   // Connect ESD or AOD here
   // Called once
   //
+  AliInfo("HERE");
   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
   if (!tree) {
     Printf("ERROR: Could not read chain from input slot 0");
   } 
   else {
+    tree->SetBranchAddress("ESDfriend.",&fESDfriend);
     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
     if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
     else fESD = esdH->GetEvent();
@@ -125,27 +141,28 @@ void AliAnalysisTaskIPInfo::CreateOutputObjects()
   fOutput = new TList;
   fOutput->SetOwner();
   //
-  TString nm = GetName();
-  nm += "_TPC";
-  fIPEst[kTPC] = new AliIntSpotEstimator(nm.Data(),fOutCut[kTPC],fNPhiBins[kTPC],
-                                        fNEstb[kTPC],fEstMin[kTPC],fEstMax[kTPC],
-                                        fNTrBins[kTPC],fNTrMin[kTPC],fNTrMax[kTPC],
-                                        fNPBins[kTPC],fPMin[kTPC],fPMax[kTPC]);
-  fIPEst[kTPC]->GetVertexer()->SetTPCMode();
-  fIPEst[kTPC]->GetVertexer()->SetConstraintOff();
-  //
-  nm = GetName();
-  nm += "_ITSTPC";
-  fIPEst[kITS] = new AliIntSpotEstimator(nm.Data(),fOutCut[kITS],fNPhiBins[kITS],
-                                        fNEstb[kITS],fEstMin[kITS],fEstMax[kITS],
-                                        fNTrBins[kITS],fNTrMin[kITS],fNTrMax[kITS],
-                                        fNPBins[kITS],fPMin[kITS],fPMax[kITS]);
-  //
-  fIPEst[kITS]->GetVertexer()->SetITSMode();
-  fIPEst[kITS]->GetVertexer()->SetConstraintOff();
+  Double_t err[3]={0.0400,0.0400,7.5};
   //
-  fOutput->Add(fIPEst[kTPC]);
-  fOutput->Add(fIPEst[kITS]);
+  for (int i=0;i<kNEst;i++) if (fNEstb[i]>1) {
+      TString nm   = GetName();
+      TString nmes = fEstNames[i];
+      nm += nmes;
+      fIPEst[i] = new AliIntSpotEstimator(nm.Data(),fOutCut[i],fNTrMinIP[i],fNPhiBins[i],
+                                         fNEstb[i],fEstMin[i],fEstMax[i],
+                                        fNTrBins[i],fNTrMin[i],fNTrMax[i],
+                                        fNPBins[i],fPMin[i],fPMax[i],fFillNt[i]);
+      AliESDVertex *initVertex = new AliESDVertex(fIPCenIni[i],err);
+      fIPEst[i]->GetVertexer()->SetVtxStart(initVertex);
+      delete initVertex;
+      fIPEst[i]->GetVertexer()->SetConstraintOff();
+      fIPEst[i]->GetVertexer()->SetMinClusters(2);
+      fIPEst[i]->SetIPCenIni(fIPCenIni[i]);
+      if (nmes == "TPC") fIPEst[i]->GetVertexer()->SetTPCMode();
+      else               fIPEst[i]->GetVertexer()->SetITSMode();
+      //
+      fOutput->Add(fIPEst[i]);
+      if (fIPEst[i]->GetNtuple()) fOutput->Add(fIPEst[i]->GetNtuple());
+    }
   //
   return;
 }
@@ -153,46 +170,126 @@ void AliAnalysisTaskIPInfo::CreateOutputObjects()
 //________________________________________________________________________
 void AliAnalysisTaskIPInfo::Exec(Option_t *) 
 {
+  static TClonesArray tracks("AliExternalTrackParam",50);
+  //
   // Main loop
   // Called for each event
-  
   if (!fESD) {
     Printf("ERROR: fESD not available");
     return;
   }
+  fESD->SetESDfriend(fESDfriend);
   //
-  fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
-  fIPEst[kITS]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
   const AliESDVertex *vtx;
   UShort_t *trackID;
   Int_t ntracks;
   //
-  // Process vertex made of TPC tracks only
-  vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
-  if (vtx) {
-    ntracks = vtx->GetNIndices();
-    trackID = (UShort_t*)vtx->GetIndices();
-    for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
-    fIPEst[kTPC]->ProcessEvent(&fTracks);
-    fTracks.Clear();
+  for (int ie=0;ie<kNEst;ie++) {
+    //
+    if (!fIPEst[ie]) continue;
+    if (ie==kTPC) {
+      fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
+      vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
+      if (vtx) {
+       ntracks = vtx->GetNIndices();
+       trackID = (UShort_t*)vtx->GetIndices();
+       for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
+       fIPEst[kTPC]->ProcessEvent(&fTracks);
+       fTracks.Clear();
+      }
+    }
+    else if (ie==kITSTPC) {
+      fIPEst[kITSTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
+      vtx = fRecoVtx[kITSTPC] ? fIPEst[kITSTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
+      if (vtx) {
+       ntracks = vtx->GetNIndices();
+       trackID = (UShort_t*)vtx->GetIndices();
+       for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
+       fIPEst[kITSTPC]->ProcessEvent(&fTracks);
+       fTracks.Clear();
+      }
+    }
+    else if (ie==kSPD) {
+      fIPEst[kSPD]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
+      ntracks = CreateSPDTracklets(tracks);
+      for (int i=ntracks;i--;) fTracks.Add((TObject*)tracks[i]);
+      fIPEst[kSPD]->ProcessEvent(&fTracks);
+      fTracks.Clear();
+    }
   }
   //
-  // Process vertex made of TPC+ITS tracks
-  vtx = fRecoVtx[kITS] ? fIPEst[kITS]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
-  if (vtx) {
-    ntracks = vtx->GetNIndices();
-    trackID = (UShort_t*)vtx->GetIndices();
-    for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
-    fIPEst[kITS]->ProcessEvent(&fTracks);
-    fTracks.Clear();
-  }
-  //  
   PostData(0, fOutput);
   //
   return;
 }      
 
 //________________________________________________________________________
+Int_t AliAnalysisTaskIPInfo::CreateSPDTracklets(TClonesArray& tracks)
+{
+  // create traclets from multiplicity class
+  double cv[21] = {
+    25e-4,
+    0   ,  25e-4,
+    0   ,     0,  40e-2,
+    0   ,     0,      0,   1e-2,
+    0   ,     0,      0,      0,   1e-2,
+    0   ,     0,      0,      0,      0,   1e-2
+  };
+  //
+  double xyzt[3],pxyz[3];
+  int nSPDtracks = 0;
+  tracks.Delete();
+  const AliMultiplicity *mult = fESD->GetMultiplicity();
+  const AliESDVertex *spdVtx  = fESD->GetPrimaryVertexSPD();
+  int nTracklets = 0;
+  if (mult && spdVtx && (nTracklets=mult->GetNumberOfTracklets())>2 ) {
+    const Double_t kRLay1=3.9, kRLay2=7.6;
+    double xv = spdVtx->GetXv();
+    double yv = spdVtx->GetYv();
+    double zv = spdVtx->GetZv();
+    for (int i=0;i<nTracklets;i++) { // get cluster coordinates from tracklet
+      double phi1 = mult->GetPhi(i);
+      double tht1 = mult->GetTheta(i);
+      double phi2 = phi1 - mult->GetDeltaPhi(i);
+      double tht2 = tht1 - mult->GetDeltaTheta(i);
+      double cs = TMath::Cos(phi1);
+      double sn = TMath::Sin(tht1);
+      double det = xv*sn+yv*sn;
+      det = det*det + kRLay1*kRLay1;
+      double t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
+      double x1 = cs*t;
+      double y1 = sn*t;
+      double z1 = zv + TMath::Sqrt(x1*x1+y1*y1)/TMath::Tan(tht1);
+      x1 += xv;
+      y1 += yv;
+      //
+      cs = TMath::Cos(phi2);
+      sn = TMath::Sin(tht2);
+      det = xv*sn+yv*sn;
+      det = det*det + kRLay2*kRLay2;
+      t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
+      double dx = cs*t;
+      double dy = sn*t;
+      double dz = zv + TMath::Sqrt(dx*dx+dy*dy)/TMath::Tan(tht2);
+      dx += xv-x1;
+      dy += yv-y1;
+      dz += -z1;
+      double dr = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
+      pxyz[0] = dx/dr; // direction cosines
+      pxyz[1] = dy/dr;
+      pxyz[2] = dz/dr;
+      t = (xv-x1)*pxyz[0] + (yv-y1)*pxyz[1] + (zv-z1)*pxyz[2];
+      xyzt[0] = x1 + t*pxyz[0];  // PCA to vertex
+      xyzt[1] = y1 + t*pxyz[1];
+      xyzt[2] = z1 + t*pxyz[2];
+      //
+      new(tracks[nSPDtracks++]) AliExternalTrackParam(xyzt,pxyz,cv,0);
+    }
+  }
+  return nSPDtracks;
+}
+
+//________________________________________________________________________
 void AliAnalysisTaskIPInfo::Terminate(Option_t *) 
 {
   // Draw result to the screen
index d747f89..d3a4a20 100644 (file)
@@ -4,6 +4,7 @@
 /* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
+class AliESDfriend;
 class AliESDEvent;
 class AliESDVertex;
 class AliIntSpotEstimator;
@@ -13,20 +14,21 @@ class AliIntSpotEstimator;
 class AliAnalysisTaskIPInfo : public AliAnalysisTask 
 {
  public:
-  enum {kTPC,kITS,kNEst};
+  enum {kITSTPC,kTPC,kSPD,kNEst};
   //
   AliAnalysisTaskIPInfo(const char *name = "IPInfo");
                        
   virtual ~AliAnalysisTaskIPInfo(); 
   //
-  AliIntSpotEstimator* GetEstTPC()                          const {return fIPEst[kTPC];}
-  AliIntSpotEstimator* GetEstITS()                          const {return fIPEst[kITS];}
-
-  void SetOptions(Int_t estID, Bool_t recoVtx=kFALSE,
-                 Double_t outcut=1e-4, Int_t nPhiBins=12,Int_t nestb=1000,
+  AliIntSpotEstimator* GetEstimator(Int_t i)                const {return i>=0&&i<kNEst ? fIPEst[i] : 0;}
+  //
+  void SetOptions(Int_t estID, Bool_t recoVtx=kFALSE, 
+                 Double_t outcut=1e-4, Int_t ntrIP=2,Int_t nPhiBins=12,Int_t nestb=1000,
                  Double_t estmin=-4e-2,Double_t estmax=6e-2,
                  Int_t ntrBins=10,Int_t ntMn=2,Int_t ntMx=32,
-                 Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3.);
+                 Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3.,Bool_t fillNt=kFALSE);
+  void SetIPCenIni(Int_t esdID, Double_t x=0,Double_t y=0,Double_t z=0);
+  Int_t CreateSPDTracklets(TClonesArray& tracks);
   //
   virtual void   ConnectInputData(Option_t *);
   virtual void   CreateOutputObjects();
@@ -37,6 +39,7 @@ class AliAnalysisTaskIPInfo : public AliAnalysisTask
   //
   // options for estimators creation
   Bool_t    fRecoVtx[kNEst];               //! request to refit the vertex for given estimator
+  Int_t     fNTrMinIP[kNEst];              //! min tracks for IP estimator
   Int_t     fNPhiBins[kNEst];              //! n bins in phi for IP
   Int_t     fNEstb[kNEst];                 //! n of estimator bins
   Int_t     fNTrBins[kNEst];               //! n of vtx.mult. bins
@@ -48,11 +51,15 @@ class AliAnalysisTaskIPInfo : public AliAnalysisTask
   Double_t  fEstMax[kNEst];                //! upper estimator boundary
   Double_t  fPMin[kNEst];                  //! lower P cut
   Double_t  fPMax[kNEst];                  //! upper P cut
+  Double_t  fIPCenIni[kNEst][3];           //! initial estimate of IP Center
+  Bool_t    fFillNt[kNEst];                //! request to fill ntuple
   //
   AliIntSpotEstimator* fIPEst[kNEst];     //! estimators
-  AliESDEvent *fESD;                      //! ESD object
+  AliESDEvent  *fESD;                     //! ESD object
+  AliESDfriend *fESDfriend;               //! ESD friend object
   TList       *fOutput;                   //! list send on output slot 0
   TObjArray fTracks;                      //! temporary storage for extracted tracks
+  static const Char_t* fEstNames[kNEst];  //  estimator names
   //
  private:    
   //
index 364f0e2..785b695 100644 (file)
@@ -6,6 +6,7 @@
 #include "AliLog.h"
 #include <TH1.h>
 #include <TH1F.h>
+#include <TNtuple.h>
 #include <TObjArray.h>
 #include <TCanvas.h>
 #include <TPaveText.h>
 ClassImp(AliIntSpotEstimator)
 
 //______________________________________________________________________________________
-AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,
+AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,Int_t ntrIP,
                                         Int_t nPhiBins,Int_t nestb,
                                         Double_t estmin,Double_t estmax,
                                         Int_t ntrBins,Int_t ntMn,Int_t ntMx,
-                                        Int_t nPBins,Double_t pmn,Double_t pmx) 
+                                        Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple) 
 : TNamed(name,""),
-  fEvProc(0),fIPCenterStat(0),fOutlierCut(outcut),fEstimIP(0),fEstimVtx(0),fEstimTrc(0),
-  fVertexer(0),fTracks(0)
+  fEvProc(0),fIPCenterStat(0),fMinTracksForIP(ntrIP>2?ntrIP:2),fOutlierCut(outcut),fEstimIP(0),
+  fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
 {
-  InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx);
+  InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
 }
 
 //______________________________________________________________________________________                                    
 AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef) 
   : TNamed("IPEstimator",""),
-    fEvProc(0),fIPCenterStat(0),fOutlierCut(1e-4),
-    fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fVertexer(0),fTracks(0)
+    fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
+    fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
 {
   if (initDef) InitEstimators();
   //
@@ -44,6 +45,7 @@ AliIntSpotEstimator::~AliIntSpotEstimator()
   delete fEstimVtx;
   delete fEstimTrc;
   delete fVertexer;
+  delete fNtuple;
 }
 
 //______________________________________________________________________________________
@@ -54,9 +56,13 @@ Bool_t AliIntSpotEstimator::ProcessIPCenter(const AliESDVertex* vtx)
   double xyz[3];
   vtx->GetXYZ(xyz);
   double r = xyz[0]*xyz[0] + xyz[1]*xyz[1];
-  if (r>1.) return kFALSE;
-  for (int i=3;i--;) fIPCenter[i] += xyz[i];
+  if (r>2.) return kFALSE;
+  for (int i=3;i--;) {
+    fIPCenter[i] += xyz[i];
+    fIPCen2[i] += xyz[i]*xyz[i];
+  }
   fIPCenterStat++;
+  fHVtxXY->Fill(xyz[0],xyz[1]);
   return kTRUE;
   //
 }
@@ -116,14 +122,13 @@ Bool_t AliIntSpotEstimator::ProcessTracks()
   UShort_t selTrackID,movTrackID=0;
   //
   AliESDVertex* recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
-  if (!recNewVtx || 
-      ((nTracks=recNewVtx->GetNIndices())<GetMinTracks()) ||  // in case some tracks were rejected
-      !ProcessIPCenter(recNewVtx)) { 
+  if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
     if (recNewVtx) delete recNewVtx;
     fTracks->Clear();
     delete[] trackID;
     return kFALSE;
   }
+  if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
   //
   double pmn = GetTrackMinP();
   double pmx = GetTrackMaxP();
@@ -152,10 +157,22 @@ Bool_t AliIntSpotEstimator::ProcessTracks()
       double phiTrack = selTrack->Phi();
       double cs = TMath::Cos(phiTrack);
       double sn = TMath::Sin(phiTrack);
-      double trDCA = xyzDCA[0]         *sn - xyzDCA[1]         *cs;  // track signed DCA to origin
-      double vtDCA = recNewVtx->GetXv()*sn - recNewVtx->GetYv()*cs;  // vertex signed DCA to origin
+      double trDCA = (xyzDCA[0]-fIPCenIni[0])         *sn - (xyzDCA[1]-fIPCenIni[1])         *cs;  // track signed DCA to origin
+      double vtDCA = (recNewVtx->GetXv()-fIPCenIni[0])*sn - (recNewVtx->GetYv()-fIPCenIni[1])*cs;  // vertex signed DCA to origin
       UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
       selTrack->PropagateTo(told,fieldVal);    // restore the track
+      if (fNtuple) {
+       static float ntf[8];
+       ntf[0] = float(nTracks1);
+       ntf[1] = recNewVtx->GetXv();
+       ntf[2] = recNewVtx->GetYv();
+       ntf[3] = recNewVtx->GetZv();
+       ntf[4] = xyzDCA[0];
+       ntf[5] = xyzDCA[1];
+       ntf[6] = phiTrack;
+       ntf[7] = pTrack;
+       fNtuple->Fill(ntf);
+      }
     }
     delete recNewVtx;
     // restore the track indices
@@ -181,7 +198,7 @@ void AliIntSpotEstimator::UpdateEstimators(double rvD, double rtD, double nTrack
   double estVtx = rvD*(rvD - rtD);
   double estTrc = rtD*(rtD - rvD);
   //
-  fEstimIP->Fill(phiTrack, estIP);
+  if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
   fEstimVtx->Fill(nTracks, estVtx);
   if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
   fEstimTrc->Fill(1./pTrack,estTrc);
@@ -191,7 +208,7 @@ void AliIntSpotEstimator::UpdateEstimators(double rvD, double rtD, double nTrack
 //______________________________________________________________________________________
 void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
                                         Int_t ntrBins,Int_t ntMn,Int_t ntMx,
-                                        Int_t nPBins,Double_t pmn,Double_t pmx)
+                                        Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
 {
   Clear();
   // regularize binning/limits for DCA resolution
@@ -227,6 +244,16 @@ void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t est
   nm = GetName();  
   nm += "dcaEst";   
   fEstimTrc = new TH2F(nm.Data(),nm.Data(),nPBins,1./pmx,1./pmn,      nestb,estmin,estmax);
+  // 
+  nm = GetName();  
+  nm += "VtxXY";   
+  fHVtxXY = new TH2F(nm.Data(),nm.Data(),200, -1,1, 200,-1,1);
+  //
+  if (ntuple) {
+    nm = GetName();  
+    nm += "ntuple";
+    fNtuple = new TNtuple(nm.Data(),nm.Data(),"ntrack:xv:yv:zv:xt:yt:phi:p");
+  }
   //
   fVertexer = new AliVertexerTracks(); // the field is set dynamically
   fVertexer->SetConstraintOff();
@@ -236,20 +263,35 @@ void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t est
 }
 
 //______________________________________________________________________________________                                                                                             
-Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin) const
+Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
 {                                                                                                                                                                                    
   // get estimate for the IP sigma
   if (!IsValid()) {AliError("Not initialized yet"); return -999;}
-  double cx  = GetIPCenter(0);
-  double cy  = GetIPCenter(1);
+  double cxe,cye;
+  double cx  = GetIPCenter(0,&cxe) - GetIPCenIni(0);
+  double cy  = GetIPCenter(1,&cye) - GetIPCenIni(1);
   TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
-  double est = CalcMean(proj, fOutlierCut) - (cx*cx + cy*cy)/2.;
+  double merr = 0;
+  double est = CalcMean(proj, fOutlierCut, &merr) - (cx*cx + cy*cy)/2.;
+  if (est>0) {
+    est = TMath::Sqrt(est);
+    if (err) {
+      *err = 0;
+      *err = merr*merr;
+      *err += cx*cx*cxe*cxe + cy*cy*cye*cye;
+      *err = TMath::Sqrt(*err)/est/2.;
+    }
+  }
+  else {
+    est = 0;
+    if (err) *err = 0;
+  }
   delete proj;
-  return est>0 ? TMath::Sqrt(est) : 0;
+  return est;
 }
 
 //______________________________________________________________________________________                                                                                             
-Double_t AliIntSpotEstimator::GetVtxSigma(int ntr) const
+Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
 {                                                                                                                                                                                    
   // get estimate for the IP sigma^2
   if (!IsValid()) {AliError("Not initialized yet"); return -999;}
@@ -260,14 +302,22 @@ Double_t AliIntSpotEstimator::GetVtxSigma(int ntr) const
     return -1;
   }
   TH1* proj = fEstimVtx->ProjectionY("vrProj",bin,bin,"e");
-  double est = CalcMean(proj, fOutlierCut);
+  double est = CalcMean(proj, fOutlierCut, err);
   delete proj;
-  return est>0 ? TMath::Sqrt(est) : 0;
+  if (est>0) {
+    est = TMath::Sqrt(est);
+    if (err) *err /= 2*est;
+  }
+  else {
+    est = 0;
+    if (err) *err = 0;
+  }
+  return est;
   //
 }
 
 //______________________________________________________________________________________                                                                                             
-Double_t AliIntSpotEstimator::GetDCASigma(double pt) const
+Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
 {                                                                                                                                                                                    
   // get estimate for the IP sigma^2
   if (!IsValid()) {AliError("Not initialized yet"); return -999;}
@@ -279,30 +329,46 @@ Double_t AliIntSpotEstimator::GetDCASigma(double pt) const
     return -1;
   }
   TH1* proj = fEstimTrc->ProjectionY("trProj",bin,bin,"e");
-  double est = CalcMean(proj, fOutlierCut);
+  double est = CalcMean(proj, fOutlierCut, err);
   delete proj;
-  return est>0 ? TMath::Sqrt(est) : 0;
+   if (est>0) {
+    est = TMath::Sqrt(est);
+    if (err) *err /= 2*est;
+  }
+  else {
+    est = 0;
+    if (err) *err = 0;
+  }
+  return est;
   //
 }
 
 //______________________________________________________________________________________                                                                                             
-Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact)
+Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
 {
   // calculate mean applying the cut on the outliers
   //
   double max = histo->GetMaximum();
-  double cut = (ctfact>0&&ctfact<1.) ? TMath::Max(1.0,max*ctfact) : 0;
+  double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
   int nb = histo->GetNbinsX();
-  double mean = 0.,cumul = 0;
+  double mean = 0.,cumul = 0, rms = 0;
   for (int i=1;i<=nb;i++) {
     double vl = histo->GetBinContent(i) - cut;
     if (vl<1e-10) continue;
-    mean  += vl*histo->GetBinCenter(i);
+    double x = histo->GetBinCenter(i);
+    mean  += vl*x;
+    rms   += vl*x*x;
     cumul += vl;
   }
   //
-  return cumul>0 ? mean/cumul : 0;
+  mean = cumul>0 ? mean/cumul : 0;
+  rms -= mean*mean*cumul;
+  if (err) {
+    *err = cumul > 1 ? rms/(cumul-1) : 0;
+    if (*err>0) *err = TMath::Sqrt(*err/cumul);
+  }
   //
+  return mean;
 }
 
 //______________________________________________________________________________________
@@ -310,29 +376,35 @@ void AliIntSpotEstimator::Print(Option_t *) const
 {   
   if (!IsValid()) {printf("Not initialized yet\n"); return;}
   //
+  double cx,cy,cz,cxe,cye,cze;
+  cx = GetIPCenter(0,&cxe);
+  cy = GetIPCenter(1,&cye);
+  cz = GetIPCenter(2,&cze);
   printf("Processed %d events\n",fEvProc);
-  printf("Estimator for IP center: %+.4e %+.4e %+.4e\n",
-        GetIPCenter(0),GetIPCenter(1),GetIPCenter(2));
-  double sgIP = GetIPSigma();
-  printf("Estimator for IP sigma : %.4e\n",sgIP);
+  printf("Estimator for IP center: %+.4e+-%.3e | %+.4e+-%.3e | %+.4e+-%.3e\n",
+        cx,cxe,cy,cye,cz,cze);
+  printf("Initial IP center was  : %+.4e %+.4e %+.4e\n",
+        GetIPCenIni(0),GetIPCenIni(1),GetIPCenIni(2));
+  double sgIP = GetIPSigma(0,&cxe);
+  printf("Estimator for IP sigma : %.4e+-%.3e\n",sgIP,cxe);
   //
   printf("Estimators for vertex resolution vs Ntracks:\n");
   
   for (int i=1;i<=GetNTrackBins();i++) {
-    double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i) );
+    double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
     if (IsZero(sig)) continue;
     int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
     int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
-    printf("%3d-%3d : %.4e\n",tmin,tmax,sig);
+    printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
   }
   //
   printf("Estimators for track DCA resolution vs P:\n");
   for (int i=1;i<=GetNPBins();i++) {
-    double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i) );
+    double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
     if (IsZero(sig)) continue;
     double pmax = 1./fEstimTrc->GetXaxis()->GetBinLowEdge(i);
     double pmin = 1./(fEstimTrc->GetXaxis()->GetBinLowEdge(i)+fEstimTrc->GetXaxis()->GetBinWidth(i));
-    printf("%.2f-%.2f : %.4e\n",pmin,pmax,sig);
+    printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
   }
 }
 
@@ -342,7 +414,7 @@ void AliIntSpotEstimator::Clear(Option_t *)
   // reset all
   fEvProc = 0;
   fIPCenterStat = 0;
-  fIPCenter[0]  = fIPCenter[1]  = fIPCenter[2]  = 0;
+  for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
   if (fEstimIP)  fEstimIP->Reset();
   if (fEstimVtx) fEstimVtx->Reset();
   if (fEstimTrc) fEstimTrc->Reset();
@@ -357,6 +429,7 @@ AliIntSpotEstimator &AliIntSpotEstimator::operator += (const AliIntSpotEstimator
   if (fEstimIP  && src.fEstimIP ) fEstimIP->Add(src.fEstimIP);
   if (fEstimVtx && src.fEstimVtx) fEstimVtx->Add(src.fEstimVtx);
   if (fEstimTrc && src.fEstimTrc) fEstimTrc->Add(src.fEstimTrc);
+  if (fHVtxXY   && src.fHVtxXY)   fHVtxXY->Add(src.fHVtxXY);
   //
   return *this;
 }
@@ -374,8 +447,6 @@ TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
   //
   char buff[200];
   //
-  double sigIPtot = GetIPSigma()*1e4;
-  //  
   cnv->Divide(2,2);
   cnv->cd(1);
   //
@@ -389,10 +460,17 @@ TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
   sprintf(buff,"Accepted  Events\n%d",fIPCenterStat);
   pt->AddText(buff);
   //
-  sprintf(buff,"Int.Spot X:Y:Z (cm)\n%+.4e : %+.4e : %+.4e",GetIPCenter(0),GetIPCenter(1),GetIPCenter(2));
+  double cx,cy,cz,cxe,cye,cze;
+  cx = GetIPCenter(0,&cxe);
+  cy = GetIPCenter(1,&cye);
+  cz = GetIPCenter(2,&cze);
+  //  
+  sprintf(buff,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
+         int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
   pt->AddText(buff);
   //
-  sprintf(buff,"Int.Spot #sigma (#mum):\n%d",int(sigIPtot));
+  cx = GetIPSigma(0,&cxe);
+  sprintf(buff,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
   pt->AddText(buff);
   pt->Draw();
   gPad->Modified();
@@ -402,7 +480,11 @@ TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
   TH1* iph = fEstimIP->ProjectionX();
   iph->Reset();
   iph->SetTitle("Int.Spot size vs #phi");
-  for (int i=1;i<=iph->GetNbinsX();i++) iph->SetBinContent(i,GetIPSigma(i)*1e4);
+  for (int i=1;i<=iph->GetNbinsX();i++) {
+    cx = GetIPSigma(i,&cxe);
+    iph->SetBinContent(i,cx*1e4);
+    iph->SetBinError(i,cxe*1e4);
+  }
   iph->GetXaxis()->SetTitle("#phi");
   iph->GetYaxis()->SetTitle("#sigma_{IP} [#mum]");
   iph->SetMarkerStyle(20);
@@ -415,7 +497,11 @@ TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
   TH1* vrh = fEstimVtx->ProjectionX();
   vrh->Reset();
   vrh->SetTitle("Vertex resolution vs N tracks");
-  for (int i=1;i<=vrh->GetNbinsX();i++) vrh->SetBinContent(i,GetVtxSigma(TMath::Nint(vrh->GetBinCenter(i)))*1e4);
+  for (int i=1;i<=vrh->GetNbinsX();i++) {
+    cx = GetVtxSigma( TMath::Nint(vrh->GetBinCenter(i)), &cxe);
+    vrh->SetBinContent(i,cx*1e4);
+    vrh->SetBinError(i,cxe*1e4);
+  }
   vrh->GetXaxis()->SetTitle("n tracks");
   vrh->GetYaxis()->SetTitle("#sigma_{VTX} [#mum]");
   vrh->SetMarkerStyle(20);
@@ -428,7 +514,11 @@ TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
   TH1* trh = fEstimTrc->ProjectionX();
   trh->Reset();
   trh->SetTitle("Track DCA resolution vs 1/P");
-  for (int i=1;i<=trh->GetNbinsX();i++) trh->SetBinContent(i,GetDCASigma(1./trh->GetBinCenter(i))*1e4);
+  for (int i=1;i<=trh->GetNbinsX();i++) {
+    cx = GetDCASigma(1./trh->GetBinCenter(i), &cxe);
+    trh->SetBinContent(i,cx*1e4);
+    trh->SetBinError(i,cxe*1e4);
+  }
   trh->GetXaxis()->SetTitle("1/p [GeV]");
   trh->GetYaxis()->SetTitle("#sigma_{DCA} [#mum]");
   trh->SetMarkerStyle(20);
@@ -463,3 +553,14 @@ Long64_t AliIntSpotEstimator::Merge(TCollection *coll)
   //
 }
 
+//_____________________________________________________________________
+Double_t AliIntSpotEstimator::GetIPCenter(Int_t id,Double_t *err) const
+{
+  // calculate IP center in axis id and error
+  double cen = fIPCenterStat>0 ? fIPCenter[id]/fIPCenterStat:0;
+  if (err) {
+    *err = fIPCenterStat>1 ? (fIPCen2[id] - cen*cen*fIPCenterStat)/(fIPCenterStat-1) : 0;
+    *err = *err > 0 ? TMath::Sqrt(*err/fIPCenterStat) : 0;
+  }
+  return cen;
+}
index 585b055..eb00e02 100644 (file)
@@ -8,6 +8,7 @@
 #include <TAxis.h>
 
 class TH1;
+class TNtuple;
 class TCanvas;
 class AliESDEvent;
 class AliESDtrack;
@@ -19,36 +20,42 @@ class AliIntSpotEstimator : public TNamed {
  public:
   //
   AliIntSpotEstimator(Bool_t initDef=kFALSE);
-  AliIntSpotEstimator(const char* name, Double_t outcut=1e-4,
+  AliIntSpotEstimator(const char* name, Double_t outcut=1e-4,Int_t ntrIP=2,
                      Int_t nPhiBins=12,Int_t nestb=500,
                      Double_t estmin=-2e-2,Double_t estmax=6e-2,
                      Int_t ntrBins=10,Int_t ntMn=2,Int_t ntMx=32,
-                     Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3.);
+                     Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3., Bool_t ntuple=kFALSE);
   ~AliIntSpotEstimator();
   AliIntSpotEstimator &operator += (const AliIntSpotEstimator &src);  
   //
   void          InitEstimators(Int_t nPhiBins=12,Int_t nestb=500,
                               Double_t estmin=-2e-2,Double_t estmax=6e-2,
                               Int_t ntrBins=10,Int_t ntMn=2,Int_t ntMx=32,
-                              Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3.);
+                              Int_t nPBins=14,Double_t pmn=0.2,Double_t pmx=3.,Bool_t ntuple=kFALSE);
   //
   Bool_t        ProcessEvent(const AliESDEvent* esd, const AliESDVertex* vtx=0);
   Bool_t        ProcessEvent(const TObjArray* tracks);
   //
-  Double_t      GetIPCenter(Int_t id)                 const  {return fIPCenterStat>0 ? fIPCenter[id]/fIPCenterStat:0;}
-  Double_t      GetIPSigma(Int_t phibin=0)            const;
-  Double_t      GetVtxSigma(int ntr)                  const;
-  Double_t      GetDCASigma(double p)                 const;
+  Double_t      GetIPCenter(Int_t id,Double_t *err=0)              const;
+  Double_t      GetIPCenIni(Int_t id)                              const  {return fIPCenIni[id];}
+  Double_t      GetIPSigma(Int_t phibin=0,Double_t *err=0)         const;
+  Double_t      GetVtxSigma(int ntr, Double_t *err=0)              const;
+  Double_t      GetDCASigma(double p, Double_t *err=0)             const;
   //
-  Int_t         GetEventsAccepted()                   const  {return fIPCenterStat;} 
-  Int_t         GetEventsProcessed()                  const  {return fEvProc;} 
+  Int_t         GetEventsAccepted()                                const  {return fIPCenterStat;} 
+  Int_t         GetEventsProcessed()                               const  {return fEvProc;} 
+  Int_t         GetMinTracksForIP()                                const  {return fMinTracksForIP;}
   //
-  void          SetOutlierCut(Double_t v=1e-4)               {fOutlierCut = v;}
+  void          SetOutlierCut(Double_t v=1e-4)                            {fOutlierCut = v;}
+  void          SetIPCenIni(Double_t *xyz)                                {for (int i=3;i--;) fIPCenIni[i] = xyz[i];}
+  void          SetMinTracksForIP(Int_t ntr=2)                            {fMinTracksForIP = ntr>2 ? ntr : 2;}
   //
-  TH2F*         GetHistoIP()                          const  {return fEstimIP;}
-  TH2F*         GetHistoVtx()                         const  {return fEstimVtx;}
-  TH2F*         GetHistoTrc()                         const  {return fEstimTrc;}
-  AliVertexerTracks* GetVertexer()                    const  {return fVertexer;}
+  TH2F*         GetHistoIP()                                       const  {return fEstimIP;}
+  TH2F*         GetHistoVtx()                                      const  {return fEstimVtx;}
+  TH2F*         GetHistoTrc()                                      const  {return fEstimTrc;}
+  TH2F*         GetHistoVtxXY()                                    const  {return fHVtxXY;}
+  TNtuple*      GetNtuple()                                        const  {return (TNtuple*)fNtuple;}
+  AliVertexerTracks* GetVertexer()                                 const  {return fVertexer;}
   //
   Int_t         GetNPhiBins()    const {return !IsValid() ? 0:fEstimIP->GetXaxis()->GetNbins();}
   Int_t         GetNTrackBins()  const {return !IsValid() ? 0:fEstimVtx->GetXaxis()->GetNbins();}
@@ -62,7 +69,7 @@ class AliIntSpotEstimator : public TNamed {
   virtual void  Print(Option_t *opt="")               const;
   virtual void  Clear(Option_t *opt="");
   virtual Long64_t Merge(TCollection *coll);
-  static Double_t CalcMean(TH1* histo, Double_t ctfact);
+  static Double_t CalcMean(TH1* histo, Double_t ctfact,Double_t *err=0);
   //
  protected:
   Bool_t        IsValid()       const {return fEstimIP!=0;}
@@ -80,12 +87,17 @@ class AliIntSpotEstimator : public TNamed {
   //
   Int_t         fEvProc;                         // number of events processed
   Int_t         fIPCenterStat;                   // number of events used for fIPCenter
+  Int_t         fMinTracksForIP;                 // account IP estimator only for vertices with >= tracks
   Double_t      fOutlierCut;                     // cut on outliers
-  Double_t      fIPCenter[3];                    // mean IP position XYZ
+  Double_t      fIPCenIni[3];                    // mean IP position XYZ (initial)
+  Double_t      fIPCenter[3];                    // IP position XYZ
+  Double_t      fIPCen2[3];                      // IP position XYZ^2
   //
   TH2F         *fEstimIP;                        // distribution of IP estimator
   TH2F         *fEstimVtx;                       // distribution of VTRes estimator
   TH2F         *fEstimTrc;                       // distribution of DCA res estimator
+  TH2F         *fHVtxXY;                         // XY histo
+  TNtuple      *fNtuple;                         //! optional ntuple with dca's
   //
   AliVertexerTracks* fVertexer;                  //! vertex fitter
   TObjArray    *fTracks;                         //! storage for processed tracks
index 1aa76ea..657d3e2 100644 (file)
@@ -14,12 +14,20 @@ AliAnalysisTaskIPInfo* AddTaskIntSpotESD()
 
   // Create the task
   AliAnalysisTaskIPInfo *taskIP = new AliAnalysisTaskIPInfo("IPInfo");
-  taskIP->SetOptions(AliAnalysisTaskIPInfo::kTPC, kFALSE, 1e-4,        12,1000,
-                    -5,5,  10,2,32,  14,0.2,3.);
-  taskIP->SetOptions(AliAnalysisTaskIPInfo::kITS, kFALSE, 1e-4,        12,1000,
-                    -4e-2,8e-2,  10,2,32,  14,0.2,3.);
-
-
+  taskIP->SetOptions(AliAnalysisTaskIPInfo::kITSTPC, kFALSE, 1e-4, 2, 12,1000,
+                    -4e-2,8e-2,  10,2,32,  14,0.2,3., kFALSE);
+  taskIP->SetIPCenIni(AliAnalysisTaskIPInfo::kITSTPC, -0.0764,0.2481,0);
+  //
+  //
+  taskIP->SetOptions(AliAnalysisTaskIPInfo::kTPC, kFALSE, 1e-4,        2, 12,1000,
+                    -4e-2,8e-2,  10,2,32,  14,0.2,3., kFALSE);
+  taskIP->SetIPCenIni(AliAnalysisTaskIPInfo::kTPC, -0.0764,0.2481,0);
+  //
+  //
+  taskIP->SetOptions(AliAnalysisTaskIPInfo::kSPD, kFALSE, 1e-4,        2, 12,1000,
+                    -4e-2,8e-2,  10,2,32,  14,0.2,3., kTRUE);
+  taskIP->SetIPCenIni(AliAnalysisTaskIPInfo::kSPD, -0.0764,0.2481,0);
+  //
   mgr->AddTask(taskIP);
 
   // Create containers for input/output
index a7fe574..1285a10 100644 (file)
@@ -15,13 +15,16 @@ void RunIPTask(const char* mode)
   TChain *chainESD = 0;
   gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
   TChain* chain = CreateESDChain(mode,-1);
+  //chain->SetBranchStatus("*ESDfriend*",1);
   //
   AliAnalysisManager *mgr  = new AliAnalysisManager("My Manager","My Manager");
   AliESDInputHandler *esdH = new AliESDInputHandler();
+  esdH->SetActiveBranches("ESDfriend");
   //
   mgr->SetInputEventHandler(esdH);
   gROOT->LoadMacro("$ALICE_ROOT/PWG1/macros/AddTaskIntSpotESD.C");
   AliAnalysisTaskIPInfo* iptask = AddTaskIntSpotESD();
+
   if(!mgr->InitAnalysis()) return;
   //
   mgr->StartAnalysis("local",chain);