#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
// 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;
+ //
}
//________________________________________________________________________
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;
fEstMax[estID] = estmax;
fPMin[estID] = pmn;
fPMax[estID] = pmx;
+ fFillNt[estID] = fillNt;
}
//________________________________________________________________________
// 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();
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;
}
//________________________________________________________________________
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 *)
{
/* 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;
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();
//
// 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
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:
//
#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();
//
delete fEstimVtx;
delete fEstimTrc;
delete fVertexer;
+ delete fNtuple;
}
//______________________________________________________________________________________
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;
//
}
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();
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
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);
//______________________________________________________________________________________
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
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();
}
//______________________________________________________________________________________
-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;}
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;}
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;
}
//______________________________________________________________________________________
{
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);
}
}
// 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();
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;
}
//
char buff[200];
//
- double sigIPtot = GetIPSigma()*1e4;
- //
cnv->Divide(2,2);
cnv->cd(1);
//
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();
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);
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);
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);
//
}
+//_____________________________________________________________________
+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;
+}
#include <TAxis.h>
class TH1;
+class TNtuple;
class TCanvas;
class AliESDEvent;
class AliESDtrack;
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();}
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;}
//
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
// 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
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);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff