From 43c9dae15a722bf623d55249bf2a480e5a3a5709 Mon Sep 17 00:00:00 2001 From: belikov Date: Tue, 11 Mar 2008 10:43:46 +0000 Subject: [PATCH] Tuning of primary vertex reco with TPC stand-alone (A. Dainese) --- STEER/AliReconstruction.cxx | 29 ++++++++ STEER/AliReconstruction.h | 4 +- STEER/AliVertexerTracks.cxx | 26 ++++--- STEER/AliVertexerTracks.h | 3 +- STEER/AliVertexerTracksTest.C | 127 ++++++++++++++++++++++------------ 5 files changed, 133 insertions(+), 56 deletions(-) diff --git a/STEER/AliReconstruction.cxx b/STEER/AliReconstruction.cxx index 0246b988c52..9d4e9004df4 100644 --- a/STEER/AliReconstruction.cxx +++ b/STEER/AliReconstruction.cxx @@ -253,6 +253,7 @@ AliReconstruction::AliReconstruction(const char* gAliceFilename, fVertexer(NULL), fDiamondProfile(NULL), + fDiamondProfileTPC(NULL), fMeanVertexConstraint(kTRUE), fGRPList(NULL), @@ -327,6 +328,7 @@ AliReconstruction::AliReconstruction(const AliReconstruction& rec) : fVertexer(NULL), fDiamondProfile(NULL), + fDiamondProfileTPC(NULL), fMeanVertexConstraint(rec.fMeanVertexConstraint), fGRPList(NULL), @@ -773,6 +775,16 @@ Bool_t AliReconstruction::Run(const char* input, Bool_t IsOnline) AliError("No diamond profile found in OCDB!"); } + entry = 0; + entry = AliCDBManager::Instance() + ->Get("GRP/Calib/MeanVertexTPC"); + + if(entry) { + fDiamondProfileTPC = dynamic_cast (entry->GetObject()); + } else { + AliError("No diamond profile found in OCDB!"); + } + AliVertexerTracks tVertexer(AliTracker::GetBz()); if(fDiamondProfile && fMeanVertexConstraint) tVertexer.SetVtxStart(fDiamondProfile); @@ -973,6 +985,9 @@ Bool_t AliReconstruction::Run(const char* input, Bool_t IsOnline) if (tpcTrack) ok = AliTracker:: PropagateTrackTo(tpcTrack,kRadius,track->GetMass(),kMaxStep,kTRUE); + + + if (ok) { Int_t n=trkArray.GetEntriesFast(); selectedIdx[n]=track->GetID(); @@ -999,6 +1014,12 @@ Bool_t AliReconstruction::Run(const char* input, Bool_t IsOnline) } if (fRunVertexFinderTracks) { // TPC + ITS primary vertex + tVertexer.SetITSrefitRequired(); + if(fDiamondProfile && fMeanVertexConstraint) { + tVertexer.SetVtxStart(fDiamondProfile); + } else { + tVertexer.SetConstraintOff(); + } AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd); if (pvtx) { if (pvtx->GetStatus()) { @@ -1011,6 +1032,12 @@ Bool_t AliReconstruction::Run(const char* input, Bool_t IsOnline) } // TPC-only primary vertex + tVertexer.SetITSrefitNotRequired(); + if(fDiamondProfileTPC && fMeanVertexConstraint) { + tVertexer.SetVtxStart(fDiamondProfileTPC); + } else { + tVertexer.SetConstraintOff(); + } pvtx=tVertexer.FindPrimaryVertex(&trkArray,selectedIdx); if (pvtx) { if (pvtx->GetStatus()) { @@ -2076,6 +2103,8 @@ void AliReconstruction::CleanUp(TFile* file, TFile* fileOld) if(!(AliCDBManager::Instance()->GetCacheFlag())) { delete fDiamondProfile; fDiamondProfile = NULL; + delete fDiamondProfileTPC; + fDiamondProfileTPC = NULL; delete fGRPList; fGRPList = NULL; } diff --git a/STEER/AliReconstruction.h b/STEER/AliReconstruction.h index 585a6a79fda..85f12d87857 100644 --- a/STEER/AliReconstruction.h +++ b/STEER/AliReconstruction.h @@ -83,6 +83,7 @@ public: void SetWriteAOD(Bool_t flag=kTRUE){fWriteAOD=flag;} void SetFillTriggerESD(Bool_t flag=kTRUE){fFillTriggerESD=flag;} void SetDiamondProfile(AliESDVertex *dp) {fDiamondProfile=dp;} + void SetDiamondProfileTPC(AliESDVertex *dp) {fDiamondProfileTPC=dp;} void SetMeanVertexConstraint(Bool_t flag=kTRUE){fMeanVertexConstraint=flag;} void SetCleanESD(Bool_t flag=kTRUE){fCleanESD=flag;} @@ -213,6 +214,7 @@ private: AliVertexer* fVertexer; //! vertexer for ITS AliTracker* fTracker[fgkNDetectors]; //! trackers AliESDVertex* fDiamondProfile; // (x,y) diamond profile for AliVertexerTracks + AliESDVertex* fDiamondProfileTPC; // (x,y) diamond profile from TPC for AliVertexerTracks Bool_t fMeanVertexConstraint; // use fDiamondProfile in AliVertexerTracks TList* fGRPList; // TList from the GRP/GRP/Data CDB folder @@ -233,7 +235,7 @@ private: // Plane Efficiency Evaluation Bool_t fRunPlaneEff ; // Evaluate Plane Efficiency - ClassDef(AliReconstruction, 20) // class for running the reconstruction + ClassDef(AliReconstruction, 21) // class for running the reconstruction }; #endif diff --git a/STEER/AliVertexerTracks.cxx b/STEER/AliVertexerTracks.cxx index f6b8014464e..47cd1474355 100644 --- a/STEER/AliVertexerTracks.cxx +++ b/STEER/AliVertexerTracks.cxx @@ -199,6 +199,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig, // read tracks from ESD Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast(); + if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig); if(nTrksOrig<=0) { if(fDebug) printf("TooFewTracks\n"); TooFewTracks(); @@ -516,11 +517,11 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig, Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast(); Int_t nTrksSel = 0; Double_t maxd0rphi; - Double_t maxd0z0 = fMaxd0z0; // default is 5 mm + Double_t maxd0z0 = (fITSrefit ? fMaxd0z0 : 10.*fMaxd0z0); Double_t sigmaCurr[3]; Double_t normdistx,normdisty; Double_t d0z0[2],covd0z0[3]; - Double_t sigma; + Double_t sigmad0; Double_t field = GetFieldkG(); AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1); @@ -533,12 +534,13 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig, track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i)); // only TPC tracks in |eta|<0.9 if(!fITSrefit && TMath::Abs(track->GetTgl())>1.) { + if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl()); delete track; continue; } // propagate track to vertex - if(optImpParCut<=1 || fOnlyFitter) { // optImpParCut==1 or 0 + if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0 track->PropagateToDCA(initVertex,field,100.,d0z0,covd0z0); } else { // optImpParCut==2 fCurrentVertex->GetSigmaXYZ(sigmaCurr); @@ -552,21 +554,26 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig, } } - sigma = TMath::Sqrt(covd0z0[0]); - maxd0rphi = fNSigma*sigma; + sigmad0 = TMath::Sqrt(covd0z0[0]); + maxd0rphi = fNSigma*sigmad0; if(optImpParCut==1) maxd0rphi *= 5.; + //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement + //maxd0z0 = 10.*fNSigma*sigmad0z0; if(fDebug) printf("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),maxd0z0); - // during iterations 1 and 2, if fConstraint=kFALSE, + // if fConstraint=kFALSE, during iterations 1 and 2 || + // if fConstraint=kTRUE, during iteration 2, // select tracks with d0oplusz0 < maxd0z0 - if(optImpParCut>=1 && !fConstraint && nTrksOrig>=3 && - fVert.GetNContributors()>0) { - if(TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]) > maxd0z0) { + if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) || + ( fConstraint && optImpParCut==2)) { + if(nTrksOrig>=3 && + TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>maxd0z0) { if(fDebug) printf(" rejected\n"); delete track; continue; } } + // select tracks with d0rphi < maxd0rphi if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { @@ -1388,3 +1395,4 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray, return fCurrentVertex; } //-------------------------------------------------------------------------- + diff --git a/STEER/AliVertexerTracks.h b/STEER/AliVertexerTracks.h index 2b06f4a8337..d7d7a8b18c2 100644 --- a/STEER/AliVertexerTracks.h +++ b/STEER/AliVertexerTracks.h @@ -112,7 +112,7 @@ class AliVertexerTracks : public TObject { Double_t fDCAcut; // maximum DCA between 2 tracks used for vertex Int_t fAlgo; // option for vertex finding algorythm Double_t fNSigma; // number of sigmas for d0 cut in PrepareTracks() - Double_t fMaxd0z0; // value [mm] for sqrt(d0d0+z0z0) cut + Double_t fMaxd0z0; // value for sqrt(d0d0+z0z0) cut // in PrepareTracks(1) if fConstraint=kFALSE Bool_t fITSrefit; // if kTRUE (default), use only kITSrefit tracks // if kFALSE, use all tracks (also TPC only) @@ -140,3 +140,4 @@ class AliVertexerTracks : public TObject { }; #endif + diff --git a/STEER/AliVertexerTracksTest.C b/STEER/AliVertexerTracksTest.C index 41093ef0371..23b2aafc79b 100644 --- a/STEER/AliVertexerTracksTest.C +++ b/STEER/AliVertexerTracksTest.C @@ -1,17 +1,19 @@ -void AliVertexerTracksTest(Double_t nSigma=3., +void AliVertexerTracksTest(TString outname="AliVertexerTracksTest.root", + Bool_t tpconly=kTRUE, Bool_t useMeanVtx=kFALSE, - Int_t itsin=1, + Double_t nSigma=3., Int_t itsrefit=1, Double_t maxd0z0=0.5, Int_t minitscls=5, Int_t mintracks=1, - TString outname="AliVertexerTracksTest.root", Bool_t onlyfit=kFALSE) { //------------------------------------------------------------------------ // Run vertex reconstruction and store results in histos and ntuple //------------------------------------------------------------------------ // WITHOUT Kinematics.root // + AliGeomManager::LoadGeometry("geometry.root"); + if (gAlice) { delete gAlice->GetRunLoader(); delete gAlice; @@ -93,8 +95,7 @@ void AliVertexerTracksTest(Double_t nSigma=3., if(!useMeanVtx) vertexer->SetConstraintOff(); if(onlyfit) vertexer->SetOnlyFitter(); vertexer->SetNSigmad0(nSigma); - if(!itsrefit) vertexer->SetITSrefitNotRequired(); - if(!itsin) vertexer->SetITSNotRequired(); + if(!itsrefit || tpconly) vertexer->SetITSrefitNotRequired(); vertexer->SetMinTracks(mintracks); vertexer->SetMinITSClusters(minitscls); vertexer->SetMaxd0z0(maxd0z0); @@ -102,11 +103,11 @@ void AliVertexerTracksTest(Double_t nSigma=3., //---------------------------------------------------------- - if(gSystem->AccessPathName("AliESDs.root",kFileExists)) { - printf("AliESDs.root not found!\n"); + if(gSystem->AccessPathName("AliESDs_itstpc.root",kFileExists)) { + printf("AliESDs_itstpc.root not found!\n"); return; } - TFile *fin = TFile::Open("AliESDs.root"); + TFile *fin = TFile::Open("AliESDs_itstpc.root"); AliESDEvent *esd = new AliESDEvent(); TTree *esdTree = (TTree*)fin->Get("esdTree"); if(!esdTree) return; @@ -119,12 +120,12 @@ void AliVertexerTracksTest(Double_t nSigma=3., for(Int_t iev=0; ievFindPrimaryVertex(esd); + AliESDVertex *vertex = 0; + if(!tpconly) { + vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd); + } else { + TObjArray trkArrayOrig(ntracks); + UShort_t *idOrig = new UShort_t[ntracks]; + const Double_t kRadius = 2.8; //something less than the beam pipe radius + const Double_t kMaxStep = 5; //max step over the material + Bool_t ok; + Int_t nTrksOrig=0; + for(Int_t i=0; iGetTrack(i); + AliExternalTrackParam *tpcTrack = + (AliExternalTrackParam *)esdt->GetTPCInnerParam(); + ok = kFALSE; + if (tpcTrack) + ok = AliTracker:: + PropagateTrackTo(tpcTrack,kRadius,esdt->GetMass(),kRadius,kTRUE); + if (ok) { + trkArrayOrig.AddLast(tpcTrack); + idOrig[nTrksOrig]=(UShort_t)esdt->GetID(); + nTrksOrig++; + } + } + + vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(&trkArrayOrig,idOrig); + delete [] idOrig; + } + nc = (Int_t)vertex->GetNContributors(); if(nc>=2) chi2red = vertex->GetChi2toNDF(); @@ -434,19 +463,26 @@ Int_t GetBinZ(Float_t z) { return 14; } //---------------------------------------------------------------------------- -void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { +void PlotVtxRes(TString inname="AliVertexerTracksTest.root", + Bool_t tpconly=kFALSE) { //--------------------------------------------------------------------- // Plot efficiency, resolutions, pulls // [at least 10k pp events are needed] //--------------------------------------------------------------------- - Float_t zMCmin=0.; - Float_t zMCmax=20.; - Float_t ncminforzshift=1.; - gStyle->SetOptStat(0); gStyle->SetOptFit(10001); + Float_t zMCmin=0.; + Float_t zMCmax=15.; + Float_t ncminforzshift=1.; + Float_t maxdiffx = 1e3; + Float_t rangefactor=1.; + if(tpconly) { + maxdiffx *= 1e2; + rangefactor = 15.; + } + Int_t nEvVtx=0; Int_t nEvNoVtx=0; Int_t ev,nUsedTrks; @@ -458,15 +494,15 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { // // HISTOGRAMS // - TH1F *hdx = new TH1F("hdx","",50,-600,600); + TH1F *hdx = new TH1F("hdx","",50,-600*rangefactor,600*rangefactor); hdx->SetXTitle("#Delta x [#mu m]"); hdx->SetYTitle("events"); // - TH1F *hdy = new TH1F("hdy","",50,-600,600); + TH1F *hdy = new TH1F("hdy","",50,-600*rangefactor,600*rangefactor); hdy->SetXTitle("#Delta y [#mu m]"); hdy->SetYTitle("events"); // - TH1F *hdz = new TH1F("hdz","",200,-600,600); + TH1F *hdz = new TH1F("hdz","",200,-600*rangefactor,600*rangefactor); hdz->SetXTitle("#Delta z [#mu m]"); hdz->SetYTitle("events"); // @@ -558,24 +594,24 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { // // OTHER HISTOGRAMS // - TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500,500); - TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500,500); - TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500,500); - TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400,400); - TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300,300); - TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300,300); - TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500,500); - TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500,500); - TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500,500); - TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400,400); - TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300,300); - TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300,300); - TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000,1000); - TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800,800); - TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800,800); - TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600); - TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500); - TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500); + TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500*rangefactor,500*rangefactor); + TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500*rangefactor,500*rangefactor); + TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500*rangefactor,500*rangefactor); + TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400*rangefactor,400*rangefactor); + TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300*rangefactor,300*rangefactor); + TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300*rangefactor,300*rangefactor); + TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500*rangefactor,500*rangefactor); + TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500*rangefactor,500*rangefactor); + TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500*rangefactor,500*rangefactor); + TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400*rangefactor,400*rangefactor); + TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300*rangefactor,300*rangefactor); + TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300*rangefactor,300*rangefactor); + TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000*rangefactor,1000*rangefactor); + TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800*rangefactor,800*rangefactor); + TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800*rangefactor,800*rangefactor); + TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600*rangefactor,600*rangefactor); + TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500*rangefactor,500*rangefactor); + TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500*rangefactor,500*rangefactor); // TH1F *hdz_z = NULL; hdz_z = new TH1F[15]; for(Int_t l=0;l<15;l++) hdz_z[l] = *hdz; @@ -624,7 +660,7 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { for(Int_t i=0; iGetEvent(i); // only triggered events - if(triggered<0.5) continue; + // if(triggered<0.5) continue; // consider only events with zMCmin<|zMC|zMCmax) continue; // @@ -636,7 +672,7 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { hntrackletsVSzMC->Fill(zMC,ntracklets); hnitstracksVSzMC->Fill(zMC,nTrks); hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6); - if(TMath::Abs(diffX) < 1000.) { // vtx found - closer than 1 mm + if(TMath::Abs(diffX) < maxdiffx) { // vtx found - closer than maxdiffx micron nEvVtx++; mult2[bin] += ntracklets; vtx[bin]++; @@ -1097,7 +1133,7 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { // res VS ntracklets TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500); - TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240); + TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240*rangefactor); fr6->SetXTitle("SPD tracklets"); fr6->SetYTitle("#sigma [#mu m]"); fr6->Draw(); @@ -1157,7 +1193,7 @@ void PlotVtxRes(TString inname="AliVertexerTracksTest.root") { ezmc[l]=1.; } TCanvas *zzz = new TCanvas("zzz","zzz",0,0,600,500); - TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50,50); + TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50*rangefactor,50*rangefactor); fr9->SetXTitle("z_{true} [cm]"); fr9->SetYTitle(" [#mu m]"); fr9->Draw(); @@ -1874,3 +1910,4 @@ void AliComputeVtxMeanFromTree(TString file="AliESDVertices.root", return; } //--------------------------------------------------------------------------- + -- 2.39.3