From d631e5e72d8149fa9f221862675131332ff0e886 Mon Sep 17 00:00:00 2001 From: belikov Date: Tue, 1 Dec 2009 13:43:51 +0000 Subject: [PATCH] Now, the estimator can work with the tracklets too (R.Shahoyan) --- PWG1/AliAnalysisTaskIPInfo.cxx | 195 +++++++++++++++++++++++-------- PWG1/AliAnalysisTaskIPInfo.h | 23 ++-- PWG1/AliIntSpotEstimator.cxx | 197 ++++++++++++++++++++++++-------- PWG1/AliIntSpotEstimator.h | 44 ++++--- PWG1/macros/AddTaskIntSpotESD.C | 20 +++- PWG1/macros/RunIPTask.C | 3 + 6 files changed, 355 insertions(+), 127 deletions(-) diff --git a/PWG1/AliAnalysisTaskIPInfo.cxx b/PWG1/AliAnalysisTaskIPInfo.cxx index e1ed42a89b1..edde24063cc 100644 --- a/PWG1/AliAnalysisTaskIPInfo.cxx +++ b/PWG1/AliAnalysisTaskIPInfo.cxx @@ -30,6 +30,7 @@ #include #include #include +#include #include #include "AliAnalysisTask.h" @@ -39,18 +40,22 @@ #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) 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 (GetInputData(0)); if (!tree) { Printf("ERROR: Could not read chain from input slot 0"); } else { + tree->SetBranchAddress("ESDfriend.",&fESDfriend); AliESDInputHandler *esdH = dynamic_cast (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;i1) { + 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,45 +170,125 @@ 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;ieGetVertexer()->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;iGetPhi(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 *) { diff --git a/PWG1/AliAnalysisTaskIPInfo.h b/PWG1/AliAnalysisTaskIPInfo.h index d747f89dcf7..d3a4a201a28 100644 --- a/PWG1/AliAnalysisTaskIPInfo.h +++ b/PWG1/AliAnalysisTaskIPInfo.h @@ -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 #include +#include #include #include #include @@ -15,23 +16,23 @@ 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())GetNIndices())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; +} diff --git a/PWG1/AliIntSpotEstimator.h b/PWG1/AliIntSpotEstimator.h index 585b055aa90..eb00e024047 100644 --- a/PWG1/AliIntSpotEstimator.h +++ b/PWG1/AliIntSpotEstimator.h @@ -8,6 +8,7 @@ #include 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 diff --git a/PWG1/macros/AddTaskIntSpotESD.C b/PWG1/macros/AddTaskIntSpotESD.C index 1aa76ea39e5..657d3e2abd3 100644 --- a/PWG1/macros/AddTaskIntSpotESD.C +++ b/PWG1/macros/AddTaskIntSpotESD.C @@ -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 diff --git a/PWG1/macros/RunIPTask.C b/PWG1/macros/RunIPTask.C index a7fe57404a1..1285a10d130 100644 --- a/PWG1/macros/RunIPTask.C +++ b/PWG1/macros/RunIPTask.C @@ -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); -- 2.39.3