--- /dev/null
+
+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)
+
+
+
--- /dev/null
+#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("../");
+}
+
--- /dev/null
+// ------------------------------------------------------
+//
+// 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)
+};
+
+
--- /dev/null
+#include "dNdEtaCorrection.h"
+
+#include <TCanvas.h>
+
+//____________________________________________________________________
+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 (nBinsVtxCount<nBinsVtx)
+ for (Int_t bx=0; bx<=nBinsX; bx++) {
+ for (Int_t by=0; by<=nBinsY; by++) {
+ if ((hEtaVsVtx_corr->GetBinContent(bx+1,by)==0)||
+ (hEtaVsVtx_corr->GetBinContent(bx-1,by)==0))
+ tmp->SetBinContent(bx,by,1);
+
+ }
+ }
+ if (nBinsEtaCount<nBinsEta)
+ for (Int_t bx=0; bx<=nBinsX; bx++) {
+ for (Int_t by=0; by<=nBinsY; by++) {
+ if ((hEtaVsVtx_corr->GetBinContent(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");
+}
--- /dev/null
+// ------------------------------------------------------
+//
+// 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)
+};
+
+
--- /dev/null
+void drawCorrection()
+{
+ gSystem->Load("libdNdEta.so");
+
+ dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
+ dNdEtaMap->LoadCorrection("correction_map.root");
+
+ dNdEtaMap->DrawHistograms();
+}
--- /dev/null
+//
+// 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; r<nDirs; r++) {
+
+ TSystemFile* presentDir = (TSystemFile*)dirList->At(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; i<nEvents; i++) {
+
+ esdBranch->GetEntry(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_mc<nPrim; i_mc++) {
+
+ TParticle* part = particleStack->Particle(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; t<nTracks; t++) {
+ AliESDtrack* esdTrack = esd->GetTrack(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();
+
+}
--- /dev/null
+//
+// 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; i<nEvents; i++) {
+
+ esdBranch->GetEntry(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; t<nTracks; t++) {
+ AliESDtrack* esdTrack = esd->GetTrack(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();
+
+}
--- /dev/null
+#include "ESDtrackQualityCuts.h"
+
+#include <Riostream.h>
+
+//____________________________________________________________________
+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; i<fNCuts; i++) cuts[i]=kFALSE;
+
+ if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
+ cuts[0]=kTRUE;
+ if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
+ cuts[1]=kTRUE;
+ if (nClustersTPC<fCut_MinNClusterTPC)
+ cuts[2]=kTRUE;
+ if (nClustersITS<fCut_MinNClusterITS)
+ cuts[3]=kTRUE;
+ if (chi2PerClusterTPC>fCut_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; i<fNCuts; i++)
+ if (cuts[i]) cut = kTRUE;
+
+ //########################################################################
+ // filling histograms
+ if (fHistogramsOn) {
+ hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
+
+ if (cut)
+ hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
+
+ for (Int_t i=0; i<fNCuts; i++) {
+ if (cuts[i])
+ hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
+
+ for (Int_t j=i; j<fNCuts; j++) {
+ if (cuts[i] && cuts[j]) {
+ Float_t x = hCutCorrelation->GetXaxis()->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; i<fNCuts; i++) {
+ hCutStatistics->GetXaxis()->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("../");
+}
+
+
+
--- /dev/null
+//****************************************************************
+//
+// 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)
+};
+
--- /dev/null
+
+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)
+
+
+
--- /dev/null
+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; i<nEvents; i++) {
+
+ esdBranch->GetEntry(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; t<nTracks; t++) {
+ AliESDtrack* esdTrack = esd->GetTrack(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();
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff