From 75ec0f41ce2f1329911197321d7b68f3a8fc1862 Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Thu, 11 May 2006 11:08:01 +0000 Subject: [PATCH] first checkin of Claus' code with a few modifications --- PWG0/dNdEta/Makefile | 72 ++++ PWG0/dNdEta/dNdEtaAnalysis.cxx | 118 +++++++ PWG0/dNdEta/dNdEtaAnalysis.h | 51 +++ PWG0/dNdEta/dNdEtaCorrection.cxx | 176 ++++++++++ PWG0/dNdEta/dNdEtaCorrection.h | 60 ++++ PWG0/dNdEta/drawCorrection.C | 9 + PWG0/dNdEta/makeCorrection.C | 240 ++++++++++++++ PWG0/dNdEta/testAnalysis.C | 194 +++++++++++ PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx | 383 ++++++++++++++++++++++ PWG0/esdTrackCuts/ESDtrackQualityCuts.h | 116 +++++++ PWG0/esdTrackCuts/Makefile | 70 ++++ PWG0/esdTrackCuts/testESDtrackCuts.C | 141 ++++++++ 12 files changed, 1630 insertions(+) create mode 100644 PWG0/dNdEta/Makefile create mode 100644 PWG0/dNdEta/dNdEtaAnalysis.cxx create mode 100644 PWG0/dNdEta/dNdEtaAnalysis.h create mode 100644 PWG0/dNdEta/dNdEtaCorrection.cxx create mode 100644 PWG0/dNdEta/dNdEtaCorrection.h create mode 100644 PWG0/dNdEta/drawCorrection.C create mode 100644 PWG0/dNdEta/makeCorrection.C create mode 100644 PWG0/dNdEta/testAnalysis.C create mode 100644 PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx create mode 100644 PWG0/esdTrackCuts/ESDtrackQualityCuts.h create mode 100644 PWG0/esdTrackCuts/Makefile create mode 100644 PWG0/esdTrackCuts/testESDtrackCuts.C diff --git a/PWG0/dNdEta/Makefile b/PWG0/dNdEta/Makefile new file mode 100644 index 00000000000..edb37819c2b --- /dev/null +++ b/PWG0/dNdEta/Makefile @@ -0,0 +1,72 @@ + +ROOTCINT = rootcint +ROOTLIBS = $(shell root-config --glibs) -lEG +ROOTCFLAGS = $(shell root-config --cflags) + +# +# Compile variables +# +LIBS = $(ROOTLIBS) +CPPFLAGS = $(ROOTCFLAGS) -I./ + +CXX = g++ +CXXFLAGS = -c -g -Wall -fPIC +LD = g++ +LDFLAGS = -rdynamic -Wl,-rpath,./ $(LIBS) -o +SOFLAGS = -shared -Wl,-soname, +CP = cp + +LIB = libdNdEta.so +LIBO = dNdEtaAnalysis.o \ + dNdEtaAnalysisCint.o \ + dNdEtaCorrection.o \ + dNdEtaCorrectionCint.o + +# +# Rules + +%Cint.cxx:%.h + rootcint -f $@ -c $(ROOTCFLAGS) $(CPPFLAGS) -p -I./ $^ + +%.o:%.cxx + $(CXX) $(CPPFLAGS) $(CXXFLAGS) $< + +%.so: + $(LD) $(SOFLAGS)$@ $^ -o $@ + +%:%.o + $(LD) $(LDFLAGS) $(LIBS) $^ -o $@ + + + +# +# +#---------------------------------------------------------------------------- + + +LIBRARIES = $(LIB) + +# +# Targets +# +.PHONY:all module + +all : $(LIBRARIES) + +clean: + @echo "Cleaning up ..." + @$(RM) -f *~ core *Cint.* *.o + +realclean: clean + @echo "Being very tidy ... " + @$(RM) -f $(LIBRARIES)* + @echo "done" + + +# +# Dependencies +# +$(LIB):$(LIBO) + + + diff --git a/PWG0/dNdEta/dNdEtaAnalysis.cxx b/PWG0/dNdEta/dNdEtaAnalysis.cxx new file mode 100644 index 00000000000..3e63e58e466 --- /dev/null +++ b/PWG0/dNdEta/dNdEtaAnalysis.cxx @@ -0,0 +1,118 @@ +#include "dNdEtaAnalysis.h" + +//____________________________________________________________________ +ClassImp(dNdEtaAnalysis); + +//____________________________________________________________________ +dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name) { + + fName = TString(name); + + hEtaVsVtx = new TH2F("eta_vs_vtx","",80,-20,20,120,-6,6); + hVtx = (TH1F*)hEtaVsVtx->ProjectionX("vtx"); + hdNdEta = (TH1F*)hEtaVsVtx->ProjectionY("dNdEta"); + + hEtaVsVtx->SetXTitle("vtx z [cm]"); + hEtaVsVtx->SetYTitle("#eta"); + hVtx->SetXTitle("vtx z [cm]"); + hdNdEta->SetXTitle("#eta"); + hdNdEta->SetYTitle("dN/d#eta"); + + hEtaVsVtx->Sumw2(); + hVtx->Sumw2(); + +} + +//____________________________________________________________________ +void +dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) { + hEtaVsVtx->Fill(vtx, eta, weight); +} + +//____________________________________________________________________ +void +dNdEtaAnalysis::FillEvent(Float_t vtx) { + hVtx->Fill(vtx); +} + +//____________________________________________________________________ +void +dNdEtaAnalysis::Finish() { + + // first normalize with n events (per vtx) + for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) { + Float_t nEvents = hVtx->GetBinContent(i_vtx); + Float_t nEventsError = hVtx->GetBinError(i_vtx); + + if (nEvents==0) continue; + + for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) { + Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents; + if (value==0) continue; + Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents; + error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/ + hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) + + TMath::Power(nEventsError/nEvents,2)); + hEtaVsVtx->SetBinContent(i_vtx, i_eta, value); + hEtaVsVtx->SetBinError(i_vtx, i_eta, error); + } + } + + // then take the wieghted average for each eta + // is this the right way to do it??? + for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) { + Float_t sum = 0; + Float_t sumError2 = 0; + Int_t nMeasurements = 0; + + Float_t sum_xw = 0; + Float_t sum_w = 0; + + for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) { + if (hVtx->GetBinContent(i_vtx)==0) continue; + if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue; + + Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2); + sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w; + sum_w = sum_w + w; + + sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta); + sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2); + nMeasurements++; + } + Float_t dndeta = 0; + Float_t error = 0; + + if (nMeasurements!=0) { + dndeta = sum/Float_t(nMeasurements); + error = TMath::Sqrt(sumError2)/Float_t(nMeasurements); + + dndeta = sum_xw/sum_w; + error = 1/TMath::Sqrt(sum_w); + + dndeta = dndeta/hdNdEta->GetBinWidth(i_eta); + error = error/hdNdEta->GetBinWidth(i_eta); + + hdNdEta->SetBinContent(i_eta, dndeta); + hdNdEta->SetBinError(i_eta, error); + } + + } +} + + + +//____________________________________________________________________ +void +dNdEtaAnalysis::SaveHistograms() { + + gDirectory->mkdir(fName.Data()); + gDirectory->cd(fName.Data()); + + hEtaVsVtx ->Write(); + hVtx ->Write(); + hdNdEta ->Write(); + + gDirectory->cd("../"); +} + diff --git a/PWG0/dNdEta/dNdEtaAnalysis.h b/PWG0/dNdEta/dNdEtaAnalysis.h new file mode 100644 index 00000000000..295f80278db --- /dev/null +++ b/PWG0/dNdEta/dNdEtaAnalysis.h @@ -0,0 +1,51 @@ +// ------------------------------------------------------ +// +// Class for dn/deta analysis +// +// ------------------------------------------------------ +// +// TODO: +// - more documentation +// - add debug statements +// - add more histograms +// - add functionality to set the bin sizes +// - figure out correct way to treat the errors +// - add functionality to make dn/deta for different mult classes? + +#ifndef ROOT_TObject +#include "TObject.h" +#endif +#ifndef ROOT_TFile +#include "TFile.h" +#endif +#ifndef ROOT_TH2 +#include "TH2.h" +#endif +#ifndef ROOT_TMath +#include "TMath.h" +#endif + +class dNdEtaAnalysis : public TObject +{ +protected: + + TString fName; + + TH2F* hEtaVsVtx; + TH1F* hVtx; + TH1F* hdNdEta; + +public: + dNdEtaAnalysis(Char_t* name="dndeta_correction"); + + void FillTrack(Float_t vtx, Float_t eta, Float_t weight); + void FillEvent(Float_t vtx); + + void Finish(); + + void SaveHistograms(); + + ClassDef(dNdEtaAnalysis,0) +}; + + diff --git a/PWG0/dNdEta/dNdEtaCorrection.cxx b/PWG0/dNdEta/dNdEtaCorrection.cxx new file mode 100644 index 00000000000..0f3a1592116 --- /dev/null +++ b/PWG0/dNdEta/dNdEtaCorrection.cxx @@ -0,0 +1,176 @@ +#include "dNdEtaCorrection.h" + +#include + +//____________________________________________________________________ +ClassImp(dNdEtaCorrection); + +//____________________________________________________________________ +dNdEtaCorrection::dNdEtaCorrection(Char_t* name) { + + fName = TString(name); + + hEtaVsVtx_meas = new TH2F("etaVsVtx_meas", "etaVsVtx_meas" ,80,-20,20,120,-6,6); + hEtaVsVtx_gene = new TH2F("etaVsVtx_gene", "etaVsVtx_gene" ,80,-20,20,120,-6,6); + hEtaVsVtx_corr = new TH2F("etaVsVtx_corr", "etaVsVtx_corr", 80,-20,20,120,-6,6); + + hEtaVsVtx_ratio = new TH2F("etaVsVtx_ratio","etaVsVtx_ratio",80,-20,20,120,-6,6); + + hEtaVsVtx_meas ->SetXTitle("vtx z [cm]"); hEtaVsVtx_meas ->SetYTitle("#eta"); + hEtaVsVtx_gene ->SetXTitle("vtx z [cm]"); hEtaVsVtx_gene ->SetYTitle("#eta"); + hEtaVsVtx_corr ->SetXTitle("vtx z [cm]"); hEtaVsVtx_corr ->SetYTitle("#eta"); + hEtaVsVtx_ratio->SetXTitle("vtx z [cm]"); hEtaVsVtx_ratio->SetYTitle("#eta"); + +} + +//____________________________________________________________________ +void +dNdEtaCorrection::Finish() { + + hEtaVsVtx_ratio->Divide(hEtaVsVtx_meas, hEtaVsVtx_gene, 1,1,"B"); + hEtaVsVtx_corr->Divide(hEtaVsVtx_gene, hEtaVsVtx_meas, 1,1,"B"); + + Int_t nBinsVtx = hEtaVsVtx_corr->GetNbinsX(); + Int_t nBinsEta = hEtaVsVtx_corr->GetNbinsY(); + + TH2F* tmp = (TH2F*)hEtaVsVtx_corr->Clone("tmp"); + + // cut at 0.2 + for (Int_t bx=0; bx<=nBinsVtx; bx++) { + for (Int_t by=0; by<=nBinsEta; by++) { + if (tmp->GetBinContent(bx,by)<0.2) { + hEtaVsVtx_corr->SetBinContent(bx,by,0); + hEtaVsVtx_corr->SetBinError(bx,by,0); + + tmp->SetBinContent(bx,by,0); + } + else + tmp->SetBinContent(bx,by,1); + } + } + +} + +//____________________________________________________________________ +void +dNdEtaCorrection::RemoveEdges(Float_t cut, Int_t nBinsVtx, Int_t nBinsEta) { + + // remove edges of correction histogram by removing + // - bins with content less than cut + // - bins next to bins with zero bin content + + Int_t nBinsX = hEtaVsVtx_corr->GetNbinsX(); + Int_t nBinsY = hEtaVsVtx_corr->GetNbinsY(); + + // set bin content to zero for bins with content smaller cut + for (Int_t bx=0; bx<=nBinsX; bx++) { + for (Int_t by=0; by<=nBinsY; by++) { + if (hEtaVsVtx_corr->GetBinContent(bx,by)>cut) { + hEtaVsVtx_corr->SetBinContent(bx,by,0); + hEtaVsVtx_corr->SetBinError(bx,by,0); + } + } + } + + // set bin content to zero for bins next to bins with zero + TH2F* tmp = (TH2F*)hEtaVsVtx_corr->Clone("tmp"); + tmp->Reset(); + + Bool_t done = kFALSE; + Int_t nBinsVtxCount = 0; + Int_t nBinsEtaCount = 0; + while (!done) { + if (nBinsVtxCountGetBinContent(bx+1,by)==0)|| + (hEtaVsVtx_corr->GetBinContent(bx-1,by)==0)) + tmp->SetBinContent(bx,by,1); + + } + } + if (nBinsEtaCountGetBinContent(bx,by+1)==0)|| + (hEtaVsVtx_corr->GetBinContent(bx,by-1)==0)) + tmp->SetBinContent(bx,by,1); + } + } + for (Int_t bx=0; bx<=nBinsX; bx++) { + for (Int_t by=0; by<=nBinsY; by++) { + if (tmp->GetBinContent(bx,by)==1) { + hEtaVsVtx_corr->SetBinContent(bx,by,0); + hEtaVsVtx_corr->SetBinError(bx,by,0); + } + } + } + nBinsVtxCount++; + nBinsEtaCount++; + if ((nBinsVtxCount>=nBinsVtx)&&(nBinsEtaCount>=nBinsEta)) done=kTRUE; + } + tmp->Delete(); + +} + +//____________________________________________________________________ +Bool_t +dNdEtaCorrection::LoadHistograms(Char_t* fileName, Char_t* dir) { + + TFile* fin = TFile::Open(fileName); + + // add test of file + // return kFALSE + + hEtaVsVtx_meas = (TH2F*)fin->Get(Form("%s/etaVsVtx_meas",dir)); + hEtaVsVtx_gene = (TH2F*)fin->Get(Form("%s/etaVsVtx_gene",dir)); + hEtaVsVtx_corr = (TH2F*)fin->Get(Form("%s/etaVsVtx_corr",dir)); + + hEtaVsVtx_ratio = (TH2F*)fin->Get(Form("%s/etaVsVtx_ratio",dir)); + + return kTRUE; +} + + +//____________________________________________________________________ +void +dNdEtaCorrection::SaveHistograms() { + + gDirectory->mkdir(fName.Data()); + gDirectory->cd(fName.Data()); + + hEtaVsVtx_meas ->Write(); + hEtaVsVtx_gene ->Write(); + + if (hEtaVsVtx_corr) + hEtaVsVtx_corr->Write(); + + if (hEtaVsVtx_ratio) + hEtaVsVtx_ratio->Write(); + + + gDirectory->cd("../"); +} + +//____________________________________________________________________ +void dNdEtaCorrection::DrawHistograms() +{ + TCanvas* canvas = new TCanvas("dNdEtaCorrection", "dNdEtaCorrection", 800, 800); + canvas->Divide(2, 2); + + canvas->cd(1); + if (hEtaVsVtx_meas) + hEtaVsVtx_meas->Draw("COLZ"); + + canvas->cd(2); + if (hEtaVsVtx_gene) + hEtaVsVtx_gene->Draw("COLZ"); + + canvas->cd(3); + if (hEtaVsVtx_ratio) + hEtaVsVtx_ratio->Draw("COLZ"); + + canvas->cd(4); + if (hEtaVsVtx_corr) + hEtaVsVtx_corr->Draw("COLZ"); +} diff --git a/PWG0/dNdEta/dNdEtaCorrection.h b/PWG0/dNdEta/dNdEtaCorrection.h new file mode 100644 index 00000000000..73001d494f4 --- /dev/null +++ b/PWG0/dNdEta/dNdEtaCorrection.h @@ -0,0 +1,60 @@ +// ------------------------------------------------------ +// +// Class to handle corrections for dN/dEta measurements +// +// ------------------------------------------------------ +// +// TODO: +// - add documentation +// - add status: generate or use maps +// - add functionality to set the bin sizes +// - add histograms with errors (for error visualization) +// + +#ifndef ROOT_TObject +#include "TObject.h" +#endif +#ifndef ROOT_TFile +#include "TFile.h" +#endif +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + + +class dNdEtaCorrection : public TObject +{ +protected: + + TString fName; + + TH2F* hEtaVsVtx_meas; + TH2F* hEtaVsVtx_gene; + + TH2F* hEtaVsVtx_corr; + TH2F* hEtaVsVtx_ratio; + +public: + dNdEtaCorrection(Char_t* name="dndeta_correction"); + + void FillMeas(Float_t vtx, Float_t eta) {hEtaVsVtx_meas->Fill(vtx, eta);} + void FillGene(Float_t vtx, Float_t eta) {hEtaVsVtx_gene->Fill(vtx, eta);} + + void Finish(); + + void SaveHistograms(); + Bool_t LoadHistograms(Char_t* fileName, Char_t* dir = "dndeta_correction"); + Bool_t LoadCorrection(Char_t* fileName, Char_t* dir = "dndeta_correction") + {return LoadHistograms(fileName, dir);} + + void DrawHistograms(); + + void RemoveEdges(Float_t cut=2, Int_t nBinsVtx=0, Int_t nBinsEta=0); + + Float_t GetCorrection(Float_t vtx, Float_t eta) + {return hEtaVsVtx_corr->GetBinContent(hEtaVsVtx_corr->FindBin(vtx,eta));} + + ClassDef(dNdEtaCorrection,0) +}; + + diff --git a/PWG0/dNdEta/drawCorrection.C b/PWG0/dNdEta/drawCorrection.C new file mode 100644 index 00000000000..e151d020214 --- /dev/null +++ b/PWG0/dNdEta/drawCorrection.C @@ -0,0 +1,9 @@ +void drawCorrection() +{ + gSystem->Load("libdNdEta.so"); + + dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection(); + dNdEtaMap->LoadCorrection("correction_map.root"); + + dNdEtaMap->DrawHistograms(); +} diff --git a/PWG0/dNdEta/makeCorrection.C b/PWG0/dNdEta/makeCorrection.C new file mode 100644 index 00000000000..374e979aafe --- /dev/null +++ b/PWG0/dNdEta/makeCorrection.C @@ -0,0 +1,240 @@ +// +// Script to make correction maps for dndeta measurements using the +// dNdEtaCorrection class. +// + +makeCorrection(Char_t* dataDir, Int_t nRuns=20) { + + Char_t str[256]; + + gSystem->Load("../esdTrackCuts/libESDtrackQuality.so"); + gSystem->Load("libdNdEta.so"); + + // ######################################################## + // selection of esd tracks + ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts(); + esdTrackCuts->DefineHistograms(1); + + esdTrackCuts->SetMinNClustersTPC(50); + esdTrackCuts->SetMaxChi2PerClusterTPC(3.5); + esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + + esdTrackCuts->SetMinNsigmaToVertex(3); + esdTrackCuts->SetAcceptKingDaughters(kFALSE); + + AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1); + + // ######################################################## + // definition of dNdEta correction object + + dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection(); + + // ######################################################## + // get the data dir + Char_t execDir[256]; + sprintf(execDir,gSystem->pwd()); + TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir); + TList* dirList = baseDir->GetListOfFiles(); + Int_t nDirs = dirList->GetEntries(); + // go back to the dir where this script is executed + gSystem->cd(execDir); + + // ######################################################## + // definition of used pointers + TFile* esdFile; + TTree* esdTree; + TBranch* esdBranch; + + AliESD* esd =0; + + // ######################################################## + // loop over runs + Int_t nRunCounter = 0; + for (Int_t r=0; rAt(r); + if (!presentDir || !presentDir->IsDirectory()) + continue; + // check that the files are there + TString currentDataDir; + currentDataDir.Form("%s/%s",dataDir, presentDir->GetName()); + cout << "Processing directory " << currentDataDir.Data() << endl; + if ((!gSystem->Which(currentDataDir,"galice.root")) || + (!gSystem->Which(currentDataDir,"AliESDs.root"))) + continue; + + if (nRunCounter++ >= nRuns) + break; + + cout << "run #" << nRunCounter << endl; + + // ######################################################### + // setup galice and runloader + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + + sprintf(str,"%s/galice.root",currentDataDir.Data()); + AliRunLoader* runLoader = AliRunLoader::Open(str); + + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + runLoader->LoadKinematics(); + runLoader->LoadHeader(); + + // ######################################################### + // open esd file and get the tree + + sprintf(str,"%s/AliESDs.root",currentDataDir.Data()); + // close it first to avoid memory leak + if (esdFile) + if (esdFile->IsOpen()) + esdFile->Close(); + + esdFile = TFile::Open(str); + esdTree = (TTree*)esdFile->Get("esdTree"); + if (!esdTree) + continue; + esdBranch = esdTree->GetBranch("ESD"); + esdBranch->SetAddress(&esd); + if (!esdBranch) + continue; + + // ######################################################## + // Magnetic field + AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field + //AliKalmanTrack::SetConvConst(1000/0.299792458/5.); + + // ######################################################## + // getting number of events + Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents(); + Int_t nESDEvents = esdBranch->GetEntries(); + + if (nEvents!=nESDEvents) { + cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl; + return; + } + // ######################################################## + // loop over number of events + cout << " looping over events..." << endl; + for(Int_t i=0; iGetEntry(i); + runLoader->GetEvent(i); + + // ######################################################## + // get the EDS vertex + AliESDVertex* vtxESD = esd->GetVertex(); + + Double_t vtx[3]; + Double_t vtx_res[3]; + vtxESD->GetXYZ(vtx); + + vtx_res[0] = vtxESD->GetXRes(); + vtx_res[1] = vtxESD->GetYRes(); + vtx_res[2] = vtxESD->GetZRes(); + + // the vertex should be reconstructed + if (strcmp(vtxESD->GetName(),"default")==0) + continue; + + // the resolution should be reasonable??? + if (vtx_res[2]==0 || vtx_res[2]>0.1) + continue; + + + // ######################################################## + // get the MC vertex + AliGenPythiaEventHeader* genHeader = + (AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader(); + + TArrayF vtxMC(3); + genHeader->PrimaryVertex(vtxMC); + Double_t vtx[3]; + vtx[0] = vtxMC[0]; + vtx[1] = vtxMC[1]; + vtx[2] = vtxMC[2]; + + + // ######################################################## + // loop over mc particles + AliStack* particleStack = runLoader->Stack(); + Int_t nPrim = particleStack->GetNprimary(); + + for (Int_t i_mc=0; i_mcParticle(i_mc); + if (!part || strcmp(part->GetName(),"XXX")==0) + continue; + + TParticlePDG* pdgPart = part->GetPDG(); + + Bool_t prim = kFALSE; + // check if it's a primary particle - is there a better way ??? + if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) { + if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0) + prim = kTRUE; + } + if (!prim) + continue; + + if (pdgPart->Charge()==0) + continue; + + Float_t eta = part->Eta(); + + if (prim) + dNdEtaMap->FillGene(vtx[2], eta); + + }// end of mc particle + + // ######################################################## + // loop over esd tracks + Int_t nTracks = esd->GetNumberOfTracks(); + for (Int_t t=0; tGetTrack(t); + + // cut the esd track? + if (!esdTrackCuts->AcceptTrack(esdTrack)) + continue; + + AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack); + if (tpcTrack->GetAlpha()==0) { + cout << " Warning esd track alpha = 0" << endl; + continue; + } + + Float_t theta = tpcTrack->Theta(); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + + // using the eta of the mc particle + Int_t label = TMath::Abs(esdTrack->GetLabel()); + if (label<0) { + cout << " cannot find corresponding mc part !!! " << label << endl; + continue; + } + TParticle* mcPart = particleStack->Particle(label); + eta = mcPart->Eta(); + + dNdEtaMap->FillMeas(vtx[2], eta); + + } // end of track loop + } // end of event loop + } // end of run loop + + dNdEtaMap->Finish(); + + TFile* fout = new TFile("correction_map.root","RECREATE"); + + esdTrackCuts->SaveHistograms("esd_track_cuts"); + dNdEtaMap->SaveHistograms(); + + fout->Write(); + fout->Close(); + + dNdEtaMap->DrawHistograms(); + +} diff --git a/PWG0/dNdEta/testAnalysis.C b/PWG0/dNdEta/testAnalysis.C new file mode 100644 index 00000000000..c509646a9ce --- /dev/null +++ b/PWG0/dNdEta/testAnalysis.C @@ -0,0 +1,194 @@ +// +// Script to test the dN/dEta analysis using the dNdEtaAnalysis and +// dNdEtaCorrection classes. Note that there is a cut on the events, +// so the measurement will be biassed. +// + +testAnalysis(Char_t* dataDir, Int_t nRuns=20) { + + Char_t str[256]; + + gSystem->Load("../esdTrackCuts/libESDtrackQuality.so"); + gSystem->Load("libdNdEta.so"); + + // ######################################################## + // selection of esd tracks + ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts(); + esdTrackCuts->DefineHistograms(1); + + esdTrackCuts->SetMinNClustersTPC(50); + esdTrackCuts->SetMaxChi2PerClusterTPC(3.5); + esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + + esdTrackCuts->SetMinNsigmaToVertex(3); + esdTrackCuts->SetAcceptKingDaughters(kFALSE); + + AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1); + + // ######################################################## + // definition of dNdEta objects + + dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection(); + + dNdEtaMap->LoadHistograms("correction_map.root","dndeta_correction"); + dNdEtaMap->RemoveEdges(2,0,2); + + dNdEtaAnalysis* analyse = new dNdEtaAnalysis("dndeta"); + + // ######################################################## + // get the data dir + Char_t execDir[256]; + sprintf(execDir,gSystem->pwd()); + TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir); + TList* dirList = baseDir->GetListOfFiles(); + Int_t nDirs = dirList->GetEntries(); + // go back to the dir where this script is executed + gSystem->cd(execDir); + + + // ######################################################## + // definition of used pointers + TFile* esdFile; + TTree* esdTree; + TBranch* esdBranch; + + AliESD* esd =0; + + // ######################################################## + // loop over runs + Int_t nRunCounter = 0; + for (Int_t r=1; r<=nDirs; r++) { + + TSystemFile* presentDir = (TSystemFile*)dirList->At(r); + if (!presentDir->IsDirectory()) + continue; + // check that the files are there + sprintf(str,"%s/%s",dataDir, presentDir->GetName()); + if ((!gSystem->Which(str,"galice.root")) || + (!gSystem->Which(str,"AliESDs.root"))) + continue; + + if (nRunCounter++ >= nRuns) + break; + + cout << "run #" << nRunCounter << endl; + + // ######################################################### + // setup galice and runloader + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + + sprintf(str,"%s/run%d/galice.root",dataDir,r); + AliRunLoader* runLoader = AliRunLoader::Open(str); + + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + runLoader->LoadKinematics(); + runLoader->LoadHeader(); + + // ######################################################### + // open esd file and get the tree + + sprintf(str,"%s/run%d/AliESDs.root",dataDir,r); + // close it first to avoid memory leak + if (esdFile) + if (esdFile->IsOpen()) + esdFile->Close(); + + esdFile = TFile::Open(str); + esdTree = (TTree*)esdFile->Get("esdTree"); + if (!esdTree) + continue; + esdBranch = esdTree->GetBranch("ESD"); + esdBranch->SetAddress(&esd); + if (!esdBranch) + continue; + + // ######################################################## + // Magnetic field + AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field + //AliKalmanTrack::SetConvConst(1000/0.299792458/5.); + + // ######################################################## + // getting number of events + Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents(); + Int_t nESDEvents = esdBranch->GetEntries(); + + if (nEvents!=nESDEvents) { + cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl; + return; + } + // ######################################################## + // loop over number of events + cout << " looping over events..." << endl; + for(Int_t i=1; iGetEntry(i); + runLoader->GetEvent(i); + + + // ######################################################## + // get the EDS vertex + AliESDVertex* vtxESD = esd->GetVertex(); + + Double_t vtx[3]; + Double_t vtx_res[3]; + vtxESD->GetXYZ(vtx); + + vtx_res[0] = vtxESD->GetXRes(); + vtx_res[1] = vtxESD->GetYRes(); + vtx_res[2] = vtxESD->GetZRes(); + + // we need a good vertex + // => there will be a bias on the measurement, since this cuts away some events + if (strcmp(vtxESD->GetName(),"default")==0) + continue; + if (vtx_res[2]==0 || vtx_res[2]>0.1) + continue; + + // ######################################################## + // loop over esd tracks + Int_t nTracks = esd->GetNumberOfTracks(); + for (Int_t t=0; tGetTrack(t); + + if (!esdTrackCuts->AcceptTrack(esdTrack)) + continue; + + AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack); + + if (tpcTrack->GetAlpha()==0) { + cout << " Warning esd track alpha = 0" << endl; + continue; + } + + Float_t theta = tpcTrack->Theta(); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta); + + dNdEtaMap->FillMeas(vtx[2], eta); + + analyse ->FillTrack(vtx[2], eta, correction); + + } // end of track loop + analyse->FillEvent(vtx[2]); + + } // end of event loop + } // end of run loop + + analyse->Finish(); + + TFile* fout = new TFile("out.root","RECREATE"); + + esdTrackCuts->SaveHistograms("esd_tracks_cuts"); + dNdEtaMap->SaveHistograms(); + analyse ->SaveHistograms(); + + fout->Write(); + fout->Close(); + +} diff --git a/PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx b/PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx new file mode 100644 index 00000000000..2a3287e7d1f --- /dev/null +++ b/PWG0/esdTrackCuts/ESDtrackQualityCuts.cxx @@ -0,0 +1,383 @@ +#include "ESDtrackQualityCuts.h" + +#include + +//____________________________________________________________________ +ClassImp(ESDtrackQualityCuts); + +//____________________________________________________________________ +ESDtrackQualityCuts::ESDtrackQualityCuts() { + + //############################################################################## + // setting default cuts + SetMinNClustersTPC(); + SetMinNClustersITS(); + SetMaxChi2PerClusterTPC(); + SetMaxChi2PerClusterITS(); + SetMaxCovDiagonalElements(); + SetRequireTPCRefit(); + SetRequireITSRefit(); + SetAcceptKingDaughters(); + SetMinNsigmaToVertex(); + SetRequireSigmaToVertex(); + + SetHistogramsOn(); + + // set the cut names + fCutNames[0] = "require TPC refit"; + fCutNames[1] = "require ITS refit"; + fCutNames[2] = "n clusters TPC"; + fCutNames[3] = "n clusters ITS"; + fCutNames[4] = "#Chi^{2}/clusters TPC"; + fCutNames[5] = "#Chi^{2}/clusters ITS"; + fCutNames[6] = "cov 11"; + fCutNames[7] = "cov 22"; + fCutNames[8] = "cov 33"; + fCutNames[9] = "cov 44"; + fCutNames[10] = "cov 55"; + fCutNames[11] = "trk-to-vtx"; + fCutNames[12] = "trk-to-vtx failed"; + fCutNames[13] = "kink daughters"; +} + +//____________________________________________________________________ +Bool_t +ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) { + + // relate the track to the new vertex + esdTrack->RelateToVertex(esdVtx, field, 999999); + + return AcceptTrack(esdTrack); +} + +//____________________________________________________________________ +Bool_t +ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) { + + AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex"); + esdTrack->RelateToVertex(esdVtx, field, 1e99); + return AcceptTrack(esdTrack); +} + +//____________________________________________________________________ +Bool_t +ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack) { + + UInt_t status = esdTrack->GetStatus(); + + // getting variables from ESD + Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt); + Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt); + + Float_t chi2PerClusterITS = -1; + Float_t chi2PerClusterTPC = -1; + if (nClustersITS!=0) + chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS); + if (nClustersTPC!=0) + chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC); + + Double_t extCov[15]; + esdTrack->GetExternalCovariance(extCov); + + Float_t b[2]; + Float_t bRes[2]; + Float_t bCov[3]; + esdTrack->GetImpactParameters(b,bCov); + if (bCov[0]<=0 || bCov[2]<=0) { + AliDebug(1, "Estimated b resolution zero!"); + bCov[0]=0; bCov[1]=0; + } + bRes[0] = TMath::Sqrt(bCov[0]); + bRes[1] = TMath::Sqrt(bCov[2]); + + // + // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + // + // this is not correct - it will not give n sigma!!! + // + // + Float_t nSigmaToVertex = -1; + if (bRes[0]!=0 && bRes[1]!=0) + nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); + + + //######################################################################## + // cut the track? + + Bool_t cuts[fNCuts]; + for (Int_t i=0; ifCut_MaxChi2PerClusterTPC) + cuts[4]=kTRUE; + if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS) + cuts[5]=kTRUE; + if (extCov[0] > fCut_MaxC11) + cuts[6]=kTRUE; + if (extCov[2] > fCut_MaxC22) + cuts[7]=kTRUE; + if (extCov[5] > fCut_MaxC33) + cuts[8]=kTRUE; + if (extCov[9] > fCut_MaxC44) + cuts[9]=kTRUE; + if (extCov[14] > fCut_MaxC55) + cuts[10]=kTRUE; + if (nSigmaToVertex > fCut_NsigmaToVertex) + cuts[11] = kTRUE; + // if n sigma could not be calculated + if (nSigmaToVertex==-1 && fCut_SigmaToVertexRequired) + cuts[12]=kTRUE; + if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) + cuts[13]=kTRUE; + + Bool_t cut=kFALSE; + for (Int_t i=0; iFill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks"))); + + if (cut) + hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks"))); + + for (Int_t i=0; iFill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i]))); + + for (Int_t j=i; jGetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i])); + Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j])); + hCutCorrelation->Fill(x,y); + } + } + } + + + hNClustersITS[0]->Fill(nClustersITS); + hNClustersTPC[0]->Fill(nClustersTPC); + hChi2PerClusterITS[0]->Fill(chi2PerClusterITS); + hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC); + + hC11[0]->Fill(extCov[0]); + hC22[0]->Fill(extCov[2]); + hC33[0]->Fill(extCov[5]); + hC44[0]->Fill(extCov[9]); + hC55[0]->Fill(extCov[14]); + + hDZ[0]->Fill(b[1]); + hDXY[0]->Fill(b[0]); + hDXYvsDZ[0]->Fill(b[1],b[0]); + + if (bRes[0]!=0 && bRes[1]!=0) { + hDZNormalized[0]->Fill(b[1]/bRes[1]); + hDXYNormalized[0]->Fill(b[0]/bRes[0]); + hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]); + } + } + + //######################################################################## + // cut the track! + if (cut) return kFALSE; + + //######################################################################## + // filling histograms after cut + if (fHistogramsOn) { + hNClustersITS[1]->Fill(nClustersITS); + hNClustersTPC[1]->Fill(nClustersTPC); + hChi2PerClusterITS[1]->Fill(chi2PerClusterITS); + hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC); + + hC11[1]->Fill(extCov[0]); + hC22[1]->Fill(extCov[2]); + hC33[1]->Fill(extCov[5]); + hC44[1]->Fill(extCov[9]); + hC55[1]->Fill(extCov[14]); + + hDZ[1]->Fill(b[1]); + hDXY[1]->Fill(b[0]); + hDXYvsDZ[1]->Fill(b[1],b[0]); + + hDZNormalized[1]->Fill(b[1]/bRes[1]); + hDXYNormalized[1]->Fill(b[0]/bRes[0]); + hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]); + } + + return kTRUE; +} + +//____________________________________________________________________ +void +ESDtrackQualityCuts::DefineHistograms(Int_t color) { + + fHistogramsOn=kTRUE; + + //################################################################################### + // defining histograms + + hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5); + + hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks"); + hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks"); + + hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);; + + for (Int_t i=0; iGetXaxis()->SetBinLabel(i+4,fCutNames[i]); + hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]); + hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]); + } + + hCutStatistics ->SetLineColor(color); + hCutCorrelation ->SetLineColor(color); + hCutStatistics ->SetLineWidth(2); + hCutCorrelation ->SetLineWidth(2); + + + hNClustersITS = new TH1F*[2]; + hNClustersTPC = new TH1F*[2]; + hChi2PerClusterITS = new TH1F*[2]; + hChi2PerClusterTPC = new TH1F*[2]; + + hC11 = new TH1F*[2]; + hC22 = new TH1F*[2]; + hC33 = new TH1F*[2]; + hC44 = new TH1F*[2]; + hC55 = new TH1F*[2]; + + hDXY = new TH1F*[2]; + hDZ = new TH1F*[2]; + hDXYvsDZ = new TH2F*[2]; + + hDXYNormalized = new TH1F*[2]; + hDZNormalized = new TH1F*[2]; + hDXYvsDZNormalized = new TH2F*[2]; + + + + Char_t str[256]; + for (Int_t i=0; i<2; i++) { + if (i==0) sprintf(str,""); + else sprintf(str,"_cut"); + + hNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5); + hNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5); + hChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10); + hChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10); + + hC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5); + hC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5); + hC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5); + hC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5); + hC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5); + + hDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10); + hDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10); + hDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10); + + hDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10); + hDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10); + hDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10); + + + hNClustersITS[i] ->SetXTitle("n ITS clusters"); + hNClustersTPC[i] ->SetXTitle("n TPC clusters"); + hChi2PerClusterITS[i] ->SetXTitle("#Chi^{2} per ITS cluster"); + hChi2PerClusterTPC[i] ->SetXTitle("#Chi^{2} per TPC cluster"); + + hC11[i] ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]"); + hC22[i] ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]"); + hC33[i] ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}"); + hC44[i] ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}"); + hC55[i] ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]"); + + hDXY[i] ->SetXTitle("transverse impact parameter"); + hDZ[i] ->SetXTitle("longitudinal impact parameter"); + hDXYvsDZ[i] ->SetXTitle("longitudinal impact parameter"); + hDXYvsDZ[i] ->SetYTitle("transverse impact parameter"); + + hDXYNormalized[i] ->SetXTitle("normalized trans impact par"); + hDZNormalized[i] ->SetXTitle("normalized long impact par"); + hDXYvsDZNormalized[i] ->SetXTitle("normalized long impact par"); + hDXYvsDZNormalized[i] ->SetYTitle("normalized trans impact par"); + + hNClustersITS[i] ->SetLineColor(color); hNClustersITS[i] ->SetLineWidth(2); + hNClustersTPC[i] ->SetLineColor(color); hNClustersTPC[i] ->SetLineWidth(2); + hChi2PerClusterITS[i] ->SetLineColor(color); hChi2PerClusterITS[i] ->SetLineWidth(2); + hChi2PerClusterTPC[i] ->SetLineColor(color); hChi2PerClusterTPC[i] ->SetLineWidth(2); + + hC11[i] ->SetLineColor(color); hC11[i] ->SetLineWidth(2); + hC22[i] ->SetLineColor(color); hC22[i] ->SetLineWidth(2); + hC33[i] ->SetLineColor(color); hC33[i] ->SetLineWidth(2); + hC44[i] ->SetLineColor(color); hC44[i] ->SetLineWidth(2); + hC55[i] ->SetLineColor(color); hC55[i] ->SetLineWidth(2); + + hDXY[i] ->SetLineColor(color); hDXY[i] ->SetLineWidth(2); + hDZ[i] ->SetLineColor(color); hDZ[i] ->SetLineWidth(2); + + hDXYNormalized[i] ->SetLineColor(color); hDXYNormalized[i] ->SetLineWidth(2); + hDZNormalized[i] ->SetLineColor(color); hDZNormalized[i] ->SetLineWidth(2); + + } +} + +//____________________________________________________________________ +void +ESDtrackQualityCuts::SaveHistograms(Char_t* dir) { + + if (!fHistogramsOn) { + AliDebug(0, "Histograms not on - cannot save histograms!!!"); + return; + } + + gDirectory->mkdir(dir); + gDirectory->cd(dir); + + gDirectory->mkdir("before_cuts"); + gDirectory->mkdir("after_cuts"); + + hCutStatistics->Write(); + hCutCorrelation->Write(); + + for (Int_t i=0; i<2; i++) { + if (i==0) + gDirectory->cd("before_cuts"); + else + gDirectory->cd("after_cuts"); + + hNClustersITS[i] ->Write(); + hNClustersTPC[i] ->Write(); + hChi2PerClusterITS[i] ->Write(); + hChi2PerClusterTPC[i] ->Write(); + + hC11[i] ->Write(); + hC22[i] ->Write(); + hC33[i] ->Write(); + hC44[i] ->Write(); + hC55[i] ->Write(); + + hDXY[i] ->Write(); + hDZ[i] ->Write(); + hDXYvsDZ[i] ->Write(); + + hDXYNormalized[i] ->Write(); + hDZNormalized[i] ->Write(); + hDXYvsDZNormalized[i] ->Write(); + + gDirectory->cd("../"); + } + + gDirectory->cd("../"); +} + + + diff --git a/PWG0/esdTrackCuts/ESDtrackQualityCuts.h b/PWG0/esdTrackCuts/ESDtrackQualityCuts.h new file mode 100644 index 00000000000..f93895aa483 --- /dev/null +++ b/PWG0/esdTrackCuts/ESDtrackQualityCuts.h @@ -0,0 +1,116 @@ +//**************************************************************** +// +// Class for handling of ESD track quality cuts +// +// TODO: +// - add functionality to save and load cuts +// - fix the n sigma cut so it is really a n sigma cut +// + +#ifndef ROOT_TObject +#include "TObject.h" +#endif +#ifndef ROOT_TTree +#include "TTree.h" +#endif +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +#include "AliESDtrack.h" +#include "AliESDVertex.h" +#include "AliLog.h" + +class ESDtrackQualityCuts : public TObject +{ +protected: + + //###################################################### + // esd track quality cuts + static const Int_t fNCuts = 14; + Char_t* fCutNames[14]; + + Int_t fCut_MinNClusterTPC; // min number of tpc clusters + Int_t fCut_MinNClusterITS; // min number of its clusters + + Float_t fCut_MaxChi2PerClusterTPC; // max tpc fit chi2 per tpc cluster + Float_t fCut_MaxChi2PerClusterITS; // max its fit chi2 per its cluster + + Float_t fCut_MaxC11; // max resolutions of covariance matrix diag. elements + Float_t fCut_MaxC22; + Float_t fCut_MaxC33; + Float_t fCut_MaxC44; + Float_t fCut_MaxC55; + + Float_t fCut_NsigmaToVertex; // max number of estimated sigma from track-to-vertex + Bool_t fCut_SigmaToVertexRequired; // cut track if sigma from track-to-vertex could not be calculated + + Bool_t fCut_AcceptKinkDaughters; // accepting kink daughters? + Bool_t fCut_RequireTPCRefit; // require TPC refit + Bool_t fCut_RequireITSRefit; // require ITS refit + + //###################################################### + // diagnostics histograms + Bool_t fHistogramsOn; + + TH1F** hNClustersITS; + TH1F** hNClustersTPC; + + TH1F** hChi2PerClusterITS; + TH1F** hChi2PerClusterTPC; + + TH1F** hC11; + TH1F** hC22; + TH1F** hC33; + TH1F** hC44; + TH1F** hC55; + + TH1F** hDXY; + TH1F** hDZ; + TH2F** hDXYvsDZ; + + TH1F** hDXYNormalized; + TH1F** hDZNormalized; + TH2F** hDXYvsDZNormalized; + + TH1F* hCutStatistics; + TH2F* hCutCorrelation; + + + // dummy array + Int_t fIdxInt[200]; + +public: + ESDtrackQualityCuts(); + + Bool_t AcceptTrack(AliESDtrack* esdTrack); + Bool_t AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field); + Bool_t AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field); + Bool_t AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Float_t field) + {return AcceptTrack(esdTrack,esdVtx, Double_t(field));} + + //###################################################### + // track quality cut setters + void SetMinNClustersTPC(Int_t min=-1) {fCut_MinNClusterTPC=min;} + void SetMinNClustersITS(Int_t min=-1) {fCut_MinNClusterITS=min;} + void SetMaxChi2PerClusterTPC(Float_t max=1e99) {fCut_MaxChi2PerClusterTPC=max;} + void SetMaxChi2PerClusterITS(Float_t max=1e99) {fCut_MaxChi2PerClusterITS=max;} + void SetMinNsigmaToVertex(Float_t sigma=3) {fCut_NsigmaToVertex=sigma;} + void SetRequireSigmaToVertex(Bool_t b=kTRUE) { fCut_SigmaToVertexRequired = b; } + void SetRequireTPCRefit(Bool_t b=kFALSE) {fCut_RequireTPCRefit=b;} + void SetRequireITSRefit(Bool_t b=kFALSE) {fCut_RequireITSRefit=b;} + void SetAcceptKingDaughters(Bool_t b=kFALSE) {fCut_AcceptKinkDaughters=b;} + void SetMaxCovDiagonalElements(Float_t c1=1e99, Float_t c2=1e99, Float_t c3=1e99, Float_t c4=1e99, Float_t c5=1e99) + {fCut_MaxC11=c1; fCut_MaxC22=c2; fCut_MaxC33=c3; fCut_MaxC44=c4; fCut_MaxC55=c5;} + + //###################################################### + void SetHistogramsOn(Bool_t b=kFALSE) {fHistogramsOn = b;} + void DefineHistograms(Int_t color=1); + void SaveHistograms(Char_t* dir="track_selection"); + + // void SaveQualityCuts(Char_t* file) + // void LoadQualityCuts(Char_t* file) + + ClassDef(ESDtrackQualityCuts,0) +}; + diff --git a/PWG0/esdTrackCuts/Makefile b/PWG0/esdTrackCuts/Makefile new file mode 100644 index 00000000000..0c772a891fa --- /dev/null +++ b/PWG0/esdTrackCuts/Makefile @@ -0,0 +1,70 @@ + +ROOTCINT = rootcint +ROOTLIBS = $(shell root-config --glibs) -lEG +ROOTCFLAGS = $(shell root-config --cflags) + +# +# Compile variables +# +LIBS = $(ROOTLIBS) -L$(ALICE_ROOT)/lib/tgt_linux/ -lESD +CPPFLAGS = $(ROOTCFLAGS) -I$(ALICE_ROOT)/include -I$(ALICE_ROOT)/STEER/ -I./ + +CXX = g++ +CXXFLAGS = -c -g -Wall -fPIC +LD = g++ +LDFLAGS = -rdynamic -Wl,-rpath,./ $(LIBS) -o +SOFLAGS = -shared -Wl,-soname, +CP = cp + +LIB = libESDtrackQuality.so +LIBO = ESDtrackQualityCutsCint.o \ + ESDtrackQualityCuts.o + +# +# Rules + +%Cint.cxx:%.h + rootcint -f $@ -c $(ROOTCFLAGS) $(CPPFLAGS) -p -I./ $^ + +%.o:%.cxx + $(CXX) $(CPPFLAGS) $(CXXFLAGS) $< + +%.so: + $(LD) $(SOFLAGS)$@ $^ -o $@ + +%:%.o + $(LD) $(LDFLAGS) $(LIBS) $^ -o $@ + + + +# +# +#---------------------------------------------------------------------------- + + +LIBRARIES = $(LIB) + +# +# Targets +# +.PHONY:all module + +all : $(LIBRARIES) + +clean: + @echo "Cleaning up ..." + @$(RM) -f *~ core *Cint.* *.o + +realclean: clean + @echo "Being very tidy ... " + @$(RM) -f $(LIBRARIES)* + @echo "done" + + +# +# Dependencies +# +$(LIB):$(LIBO) + + + diff --git a/PWG0/esdTrackCuts/testESDtrackCuts.C b/PWG0/esdTrackCuts/testESDtrackCuts.C new file mode 100644 index 00000000000..2d75edae45c --- /dev/null +++ b/PWG0/esdTrackCuts/testESDtrackCuts.C @@ -0,0 +1,141 @@ +testESDtrackCuts(Char_t* dataDir=, Int_t nRuns=10) { + + Char_t str[256]; + + gSystem->Load("libESDtrackQuality.so"); + + // ######################################################## + // definition of ESD track cuts + + ESDtrackQualityCuts* trackCuts = new ESDtrackQualityCuts(); + trackCuts->DefineHistograms(4); + + trackCuts->SetMinNClustersTPC(50); + trackCuts->SetMaxChi2PerClusterTPC(3.5); + trackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + trackCuts->SetRequireTPCRefit(kTRUE); + + trackCuts->SetMinNsigmaToVertex(3); + trackCuts->SetAcceptKingDaughters(kFALSE); + + AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1); + + // ######################################################## + // definition of used pointers + TFile* esdFile; + TTree* esdTree; + TBranch* esdBranch; + + AliESD* esd = 0; + + // ######################################################## + // get the data dir + Char_t execDir[256]; + sprintf(execDir,gSystem->pwd()); + TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir); + TList* dirList = baseDir->GetListOfFiles(); + Int_t nDirs = dirList->GetEntries(); + // go back to the dir where this script is executed + gSystem->cd(execDir); + + // ######################################################## + // loop over runs + Int_t nRunCounter = 0; + for (Int_t r=1; r<=nDirs; r++) { + + TSystemFile* presentDir = (TSystemFile*)dirList->At(r); + if (!presentDir->IsDirectory()) + continue; + // first check that the files are there + sprintf(str,"%s/%s",dataDir, presentDir->GetName()); + if ((!gSystem->Which(str,"galice.root")) || + (!gSystem->Which(str,"AliESDs.root"))) + continue; + + if (nRunCounter++ >= nRuns) + break; + + cout << "run #" << nRunCounter << endl; + + // ######################################################### + // setup galice and runloader + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + + sprintf(str,"%s/run%d/galice.root",dataDir,r); + AliRunLoader* runLoader = AliRunLoader::Open(str); + + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + runLoader->LoadHeader(); + + // ######################################################### + // open esd file and get the tree + + sprintf(str,"%s/run%d/AliESDs.root",dataDir,r); + // close it first to avoid memory leak + if (esdFile) + if (esdFile->IsOpen()) + esdFile->Close(); + + esdFile = TFile::Open(str); + esdTree = (TTree*)esdFile->Get("esdTree"); + if (!esdTree) + continue; + esdBranch = esdTree->GetBranch("ESD"); + esdBranch->SetAddress(&esd); + if (!esdBranch) + continue; + + // ######################################################## + // Magnetic field + AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field + + // ######################################################## + // getting number of events + Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents(); + Int_t nESDEvents = esdBranch->GetEntries(); + + if (nEvents!=nESDEvents) + cout << " Warning: Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl; + + // ######################################################## + // loop over number of events + cout << " looping over events..." << endl; + for(Int_t i=1; iGetEntry(i); + runLoader->GetEvent(i); + + // ######################################################## + // get the EDS vertex + AliESDVertex* vtxESD = esd->GetVertex(); + + Double_t vtxSigma[3]; + vtxESD->GetSigmaXYZ(vtxSigma); + + // ######################################################## + // loop over esd tracks + Int_t nTracks = esd->GetNumberOfTracks(); + + for (Int_t t=0; tGetTrack(t); + + //trackCuts->AcceptTrack(esdTrack, vtxESD, esd->GetMagneticField()); + trackCuts->AcceptTrack(esdTrack); + + } // end of track loop + } // end of event loop + } // end of run loop + + TFile* fout = new TFile("out.root","RECREATE"); + + trackCuts->SaveHistograms("esd_track_cuts"); + + fout->Write(); + fout->Close(); + +} -- 2.39.3