X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG0%2FesdTrackCuts%2FAliESDtrackCuts.cxx;h=b0fa62775058f255cc55c1611926ed556e2281cb;hb=70d782ef4ffa19c13c7a87a9e517afbd3095bc41;hp=de124fd8209c17473b4975a7249d62ffb2669273;hpb=10ebe68d7e3b94ddd4c0617414019e675ee19721;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG0/esdTrackCuts/AliESDtrackCuts.cxx b/PWG0/esdTrackCuts/AliESDtrackCuts.cxx index de124fd8209..b0fa6277505 100644 --- a/PWG0/esdTrackCuts/AliESDtrackCuts.cxx +++ b/PWG0/esdTrackCuts/AliESDtrackCuts.cxx @@ -1,10 +1,31 @@ -#include "AliESDtrackCuts.h" +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ +#include "AliESDtrackCuts.h" #include #include +#include #include +#include +#include +#include + //____________________________________________________________________ ClassImp(AliESDtrackCuts) @@ -24,7 +45,6 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = { "trk-to-vtx", "trk-to-vtx failed", "kink daughters", - "p", "p_{T}", "p_{x}", @@ -35,7 +55,38 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = { }; //____________________________________________________________________ -AliESDtrackCuts::AliESDtrackCuts() +AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title), + fCutMinNClusterTPC(0), + fCutMinNClusterITS(0), + fCutMaxChi2PerClusterTPC(0), + fCutMaxChi2PerClusterITS(0), + fCutMaxC11(0), + fCutMaxC22(0), + fCutMaxC33(0), + fCutMaxC44(0), + fCutMaxC55(0), + fCutAcceptKinkDaughters(0), + fCutRequireTPCRefit(0), + fCutRequireITSRefit(0), + fCutNsigmaToVertex(0), + fCutSigmaToVertexRequired(0), + fPMin(0), + fPMax(0), + fPtMin(0), + fPtMax(0), + fPxMin(0), + fPxMax(0), + fPyMin(0), + fPyMax(0), + fPzMin(0), + fPzMax(0), + fEtaMin(0), + fEtaMax(0), + fRapMin(0), + fRapMax(0), + fHistogramsOn(0), + fhCutStatistics(0), + fhCutCorrelation(0) { // // constructor @@ -67,7 +118,39 @@ AliESDtrackCuts::AliESDtrackCuts() } //_____________________________________________________________________________ -AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : TObject(c) +AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c), + fCutMinNClusterTPC(0), + fCutMinNClusterITS(0), + fCutMaxChi2PerClusterTPC(0), + fCutMaxChi2PerClusterITS(0), + fCutMaxC11(0), + fCutMaxC22(0), + fCutMaxC33(0), + fCutMaxC44(0), + fCutMaxC55(0), + fCutAcceptKinkDaughters(0), + fCutRequireTPCRefit(0), + fCutRequireITSRefit(0), + fCutNsigmaToVertex(0), + fCutSigmaToVertexRequired(0), + fPMin(0), + fPMax(0), + fPtMin(0), + fPtMax(0), + fPxMin(0), + fPxMax(0), + fPyMin(0), + fPyMax(0), + fPzMin(0), + fPzMax(0), + fEtaMin(0), + fEtaMax(0), + fRapMin(0), + fRapMax(0), + fHistogramsOn(0), + ffDTheoretical(0), + fhCutStatistics(0), + fhCutCorrelation(0) { // // copy constructor @@ -82,7 +165,51 @@ AliESDtrackCuts::~AliESDtrackCuts() // destructor // - // ## TODO to be implemented + for (Int_t i=0; i<2; i++) { + + if (fhNClustersITS[i]) + delete fhNClustersITS[i]; + if (fhNClustersTPC[i]) + delete fhNClustersTPC[i]; + if (fhChi2PerClusterITS[i]) + delete fhChi2PerClusterITS[i]; + if (fhChi2PerClusterTPC[i]) + delete fhChi2PerClusterTPC[i]; + if (fhC11[i]) + delete fhC11[i]; + if (fhC22[i]) + delete fhC22[i]; + if (fhC33[i]) + delete fhC33[i]; + if (fhC44[i]) + delete fhC44[i]; + if (fhC55[i]) + delete fhC55[i]; + + if (fhDXY[i]) + delete fhDXY[i]; + if (fhDZ[i]) + delete fhDZ[i]; + if (fhDXYvsDZ[i]) + delete fhDXYvsDZ[i]; + + if (fhDXYNormalized[i]) + delete fhDXYNormalized[i]; + if (fhDZNormalized[i]) + delete fhDZNormalized[i]; + if (fhDXYvsDZNormalized[i]) + delete fhDXYvsDZNormalized[i]; + if (fhNSigmaToVertex[i]) + delete fhNSigmaToVertex[i]; + } + + if (ffDTheoretical) + delete ffDTheoretical; + + if (fhCutStatistics) + delete fhCutStatistics; + if (fhCutCorrelation) + delete fhCutCorrelation; } void AliESDtrackCuts::Init() @@ -148,7 +275,9 @@ void AliESDtrackCuts::Init() fhDXYNormalized[i] = 0; fhDZNormalized[i] = 0; fhDXYvsDZNormalized[i] = 0; + fhNSigmaToVertex[i] = 0; } + ffDTheoretical = 0; fhCutStatistics = 0; fhCutCorrelation = 0; @@ -233,56 +362,93 @@ void AliESDtrackCuts::Copy(TObject &c) const if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone(); if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone(); if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone(); + if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone(); } + if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone(); if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone(); if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone(); - TObject::Copy(c); + TNamed::Copy(c); } -//____________________________________________________________________ -Bool_t -AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { - // - // figure out if the tracks survives all the track cuts defined - // - // the different quality parameter and kinematic values are first - // retrieved from the track. then it is found out what cuts the - // track did not survive and finally the cuts are imposed. +//_____________________________________________________________________________ +Long64_t AliESDtrackCuts::Merge(TCollection* list) { + // Merge a list of AliESDtrackCuts objects with this (needed for PROOF) + // Returns the number of merged objects (including this) + if (!list) + return 0; + + if (list->IsEmpty()) + return 1; + if (!fHistogramsOn) + return 0; - UInt_t status = esdTrack->GetStatus(); + TIterator* iter = list->MakeIterator(); + TObject* obj; - // dummy array - Int_t fIdxInt[200]; - // getting quality parameters from the ESD track - Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt); - Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt); - + // collection of measured and generated histograms + Int_t count = 0; + while ((obj = iter->Next())) { + AliESDtrackCuts* entry = dynamic_cast(obj); + if (entry == 0) + continue; - 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); + if (!entry->fHistogramsOn) + continue; + + for (Int_t i=0; i<2; i++) { + + fhNClustersITS[i] ->Add(entry->fhNClustersITS[i] ); + fhNClustersTPC[i] ->Add(entry->fhNClustersTPC[i] ); + + fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); + fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); + + fhC11[i] ->Add(entry->fhC11[i] ); + fhC22[i] ->Add(entry->fhC22[i] ); + fhC33[i] ->Add(entry->fhC33[i] ); + fhC44[i] ->Add(entry->fhC44[i] ); + fhC55[i] ->Add(entry->fhC55[i] ); + + fhDXY[i] ->Add(entry->fhDXY[i] ); + fhDZ[i] ->Add(entry->fhDZ[i] ); + fhDXYvsDZ[i] ->Add(entry->fhDXYvsDZ[i] ); + + fhDXYNormalized[i] ->Add(entry->fhDXYNormalized[i] ); + fhDZNormalized[i] ->Add(entry->fhDZNormalized[i] ); + fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]); + fhNSigmaToVertex[i] ->Add(entry->fhNSigmaToVertex[i]); + + } + + fhCutStatistics ->Add(entry->fhCutStatistics); + fhCutCorrelation ->Add(entry->fhCutCorrelation); + + count++; + } - Double_t extCov[15]; - esdTrack->GetExternalCovariance(extCov); + return count+1; +} + + +//____________________________________________________________________ +Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack) +{ + // Calculates the number of sigma to the vertex. - // getting the track to vertex parameters Float_t b[2]; Float_t bRes[2]; Float_t bCov[3]; - esdTrack->GetImpactParameters(b,bCov); + esdTrack->GetImpactParameters(b,bCov); if (bCov[0]<=0 || bCov[2]<=0) { AliDebug(1, "Estimated b resolution lower or equal zero!"); bCov[0]=0; bCov[2]=0; - } + } bRes[0] = TMath::Sqrt(bCov[0]); bRes[1] = TMath::Sqrt(bCov[2]); @@ -296,22 +462,93 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { // // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2) // Can this be expressed in a different way? + + if (bRes[0] == 0 || bRes[1] ==0) + return -1; + + Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); + + // stupid rounding problem screws up everything: + // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :( + if (TMath::Exp(-d * d / 2) < 1e-10) + return 1000; + + d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); + return d; +} + +void AliESDtrackCuts::EnableNeededBranches(TTree* tree) +{ + // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack + + tree->SetBranchStatus("fTracks.fFlags", 1); + tree->SetBranchStatus("fTracks.fITSncls", 1); + tree->SetBranchStatus("fTracks.fTPCncls", 1); + tree->SetBranchStatus("fTracks.fITSchi2", 1); + tree->SetBranchStatus("fTracks.fTPCchi2", 1); + tree->SetBranchStatus("fTracks.fC*", 1); + tree->SetBranchStatus("fTracks.fD", 1); + tree->SetBranchStatus("fTracks.fZ", 1); + tree->SetBranchStatus("fTracks.fCdd", 1); + tree->SetBranchStatus("fTracks.fCdz", 1); + tree->SetBranchStatus("fTracks.fCzz", 1); + tree->SetBranchStatus("fTracks.fP*", 1); + tree->SetBranchStatus("fTracks.fR*", 1); + tree->SetBranchStatus("fTracks.fKinkIndexes*", 1); +} + +//____________________________________________________________________ +Bool_t +AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { + // + // figure out if the tracks survives all the track cuts defined // - // - // FIX: I don't think this is correct!!! Keeping d as n_sigma for now... + // the different quality parameter and kinematic values are first + // retrieved from the track. then it is found out what cuts the + // track did not survive and finally the cuts are imposed. - Float_t nSigmaToVertex = -1; - if (bRes[0]!=0 && bRes[1]!=0) { - Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); - nSigmaToVertex = d; + // this function needs the following branches: + // fTracks.fFlags + // fTracks.fITSncls + // fTracks.fTPCncls + // fTracks.fITSchi2 + // fTracks.fTPCchi2 + // fTracks.fC //GetExternalCovariance + // fTracks.fD //GetImpactParameters + // fTracks.fZ //GetImpactParameters + // fTracks.fCdd //GetImpactParameters + // fTracks.fCdz //GetImpactParameters + // fTracks.fCzz //GetImpactParameters + // fTracks.fP //GetPxPyPz + // fTracks.fR //GetMass + // fTracks.fP //GetMass + // fTracks.fKinkIndexes - // JF solution: - //nSigmaToVertex = TMath::ErfInverse(TMath::Sqrt(2.0/TMath::Pi()) * TMath::Erf(d / TMath::Sqrt(2))) * TMath::Sqrt(2); - // Claus solution: - //nSigmaToVertex = TMath::Sqrt(2)*(TMath::ErfInverse(1 - TMath::Exp(0.5*(-d*d)))); - } + UInt_t status = esdTrack->GetStatus(); - // getting the kinematic variables of the track + // dummy array + Int_t fIdxInt[200]; + + // getting quality parameters from the ESD track + 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); + + // getting the track to vertex parameters + Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack); + + // getting the kinematic variables of the track // (assuming the mass is known) Double_t p[3]; esdTrack->GetPxPyPz(p); @@ -340,7 +577,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { cuts[0]=kTRUE; if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0) cuts[1]=kTRUE; - if (nClustersTPC fCutMaxC55) cuts[10]=kTRUE; - if (nSigmaToVertex > fCutNsigmaToVertex) + if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired) cuts[11] = kTRUE; // if n sigma could not be calculated if (nSigmaToVertex<0 && fCutSigmaToVertexRequired) cuts[12]=kTRUE; - if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) + if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) cuts[13]=kTRUE; // track kinematics cut if((momentum < fPMin) || (momentum > fPMax)) @@ -374,7 +611,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { cuts[16] = kTRUE; if((p[1] < fPyMin) || (p[1] > fPyMax)) cuts[17] = kTRUE; - if((p[2] < fPzMin) || (p[2] > fPzMax)) + if((p[2] < fPzMin) || (p[2] > fPzMax)) cuts[18] = kTRUE; if((eta < fEtaMin) || (eta > fEtaMax)) cuts[19] = kTRUE; @@ -407,48 +644,71 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { } - fhNClustersITS[0]->Fill(nClustersITS); - fhNClustersTPC[0]->Fill(nClustersTPC); + fhNClustersITS[0]->Fill(nClustersITS); + fhNClustersTPC[0]->Fill(nClustersTPC); fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS); - fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC); - - fhC11[0]->Fill(extCov[0]); - fhC22[0]->Fill(extCov[2]); - fhC33[0]->Fill(extCov[5]); - fhC44[0]->Fill(extCov[9]); - fhC55[0]->Fill(extCov[14]); - - fhDZ[0]->Fill(b[1]); - fhDXY[0]->Fill(b[0]); + fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC); + + fhC11[0]->Fill(extCov[0]); + fhC22[0]->Fill(extCov[2]); + fhC33[0]->Fill(extCov[5]); + fhC44[0]->Fill(extCov[9]); + fhC55[0]->Fill(extCov[14]); + + 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 lower or equal zero!"); + bCov[0]=0; bCov[2]=0; + } + bRes[0] = TMath::Sqrt(bCov[0]); + bRes[1] = TMath::Sqrt(bCov[2]); + + fhDZ[0]->Fill(b[1]); + fhDXY[0]->Fill(b[0]); fhDXYvsDZ[0]->Fill(b[1],b[0]); if (bRes[0]!=0 && bRes[1]!=0) { fhDZNormalized[0]->Fill(b[1]/bRes[1]); - fhDXYNormalized[0]->Fill(b[0]/bRes[0]); + fhDXYNormalized[0]->Fill(b[0]/bRes[0]); fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]); + fhNSigmaToVertex[0]->Fill(nSigmaToVertex); } } - //######################################################################## + //######################################################################## // cut the track! if (cut) return kFALSE; - //######################################################################## + //######################################################################## // filling histograms after cut if (fHistogramsOn) { - fhNClustersITS[1]->Fill(nClustersITS); - fhNClustersTPC[1]->Fill(nClustersTPC); + fhNClustersITS[1]->Fill(nClustersITS); + fhNClustersTPC[1]->Fill(nClustersTPC); fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS); - fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC); - - fhC11[1]->Fill(extCov[0]); - fhC22[1]->Fill(extCov[2]); + fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC); + + fhC11[1]->Fill(extCov[0]); + fhC22[1]->Fill(extCov[2]); fhC33[1]->Fill(extCov[5]); - fhC44[1]->Fill(extCov[9]); - fhC55[1]->Fill(extCov[14]); - - fhDZ[1]->Fill(b[1]); - fhDXY[1]->Fill(b[0]); + fhC44[1]->Fill(extCov[9]); + fhC55[1]->Fill(extCov[14]); + + 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 lower or equal zero!"); + bCov[0]=0; bCov[2]=0; + } + bRes[0] = TMath::Sqrt(bCov[0]); + bRes[1] = TMath::Sqrt(bCov[2]); + + fhDZ[1]->Fill(b[1]); + fhDXY[1]->Fill(b[0]); fhDXYvsDZ[1]->Fill(b[1],b[0]); if (bRes[0]!=0 && bRes[1]!=0) @@ -456,15 +716,55 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) { fhDZNormalized[1]->Fill(b[1]/bRes[1]); fhDXYNormalized[1]->Fill(b[0]/bRes[0]); fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]); + fhNSigmaToVertex[1]->Fill(nSigmaToVertex); } } - + return kTRUE; } //____________________________________________________________________ -TObjArray* -AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) +TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) +{ + // + // returns an array of all tracks that pass the cuts + // + + TObjArray* acceptedTracks = new TObjArray(); + + // loop over esd tracks + for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) { + AliESDtrack* track = esd->GetTrack(iTrack); + + if (AcceptTrack(track)) + acceptedTracks->Add(track); + } + + return acceptedTracks; +} + +//____________________________________________________________________ +Int_t AliESDtrackCuts::CountAcceptedTracks(AliESD* esd) +{ + // + // returns an the number of tracks that pass the cuts + // + + Int_t count = 0; + + // loop over esd tracks + for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) { + AliESDtrack* track = esd->GetTrack(iTrack); + + if (AcceptTrack(track)) + count++; + } + + return count; +} + +//____________________________________________________________________ +TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd) { // // returns an array of all tracks that pass the cuts @@ -484,8 +784,7 @@ AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) } //____________________________________________________________________ -Int_t -AliESDtrackCuts::CountAcceptedTracks(AliESD* esd) +Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd) { // // returns an the number of tracks that pass the cuts @@ -511,7 +810,10 @@ AliESDtrackCuts::CountAcceptedTracks(AliESD* esd) // fHistogramsOn=kTRUE; - + + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + //################################################################################### // defining histograms @@ -538,46 +840,48 @@ AliESDtrackCuts::CountAcceptedTracks(AliESD* esd) if (i==0) sprintf(str," "); else sprintf(str,"_cut"); - fhNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5); - fhNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5); + fhNClustersITS[i] = new TH1F(Form("nClustersITS%s",str) ,"",8,-0.5,7.5); + fhNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str) ,"",165,-0.5,164.5); fhChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10); fhChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10); - fhC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5); - fhC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5); - fhC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5); + fhC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",2000,0,20); + fhC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",2000,0,20); + fhC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,1); fhC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5); fhC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5); - fhDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10); - fhDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10); + fhDXY[i] = new TH1F(Form("dXY%s",str) ,"",500,-10,10); + fhDZ[i] = new TH1F(Form("dZ%s",str) ,"",500,-10,10); fhDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10); - fhDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10); - fhDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10); + fhDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str) ,"",500,-10,10); + fhDZNormalized[i] = new TH1F(Form("dZNormalized%s",str) ,"",500,-10,10); fhDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10); + fhNSigmaToVertex[i] = new TH1F(Form("nSigmaToVertex%s",str),"",500,0,50); - fhNClustersITS[i]->SetXTitle("n ITS clusters"); - fhNClustersTPC[i]->SetXTitle("n TPC clusters"); - fhChi2PerClusterITS[i]->SetXTitle("#Chi^{2} per ITS cluster"); - fhChi2PerClusterTPC[i]->SetXTitle("#Chi^{2} per TPC cluster"); + fhNClustersITS[i]->SetTitle("n ITS clusters"); + fhNClustersTPC[i]->SetTitle("n TPC clusters"); + fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster"); + fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster"); - fhC11[i]->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]"); - fhC22[i]->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]"); - fhC33[i]->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}"); - fhC44[i]->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}"); - fhC55[i]->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]"); + fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]"); + fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]"); + fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}"); + fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}"); + fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]"); - fhDXY[i]->SetXTitle("transverse impact parameter"); - fhDZ[i]->SetXTitle("longitudinal impact parameter"); - fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter"); + fhDXY[i]->SetTitle("transverse impact parameter"); + fhDZ[i]->SetTitle("longitudinal impact parameter"); + fhDXYvsDZ[i]->SetTitle("longitudinal impact parameter"); fhDXYvsDZ[i]->SetYTitle("transverse impact parameter"); - fhDXYNormalized[i]->SetXTitle("normalized trans impact par"); - fhDZNormalized[i]->SetXTitle("normalized long impact par"); - fhDXYvsDZNormalized[i]->SetXTitle("normalized long impact par"); + fhDXYNormalized[i]->SetTitle("normalized trans impact par"); + fhDZNormalized[i]->SetTitle("normalized long impact par"); + fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par"); fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par"); + fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex"); fhNClustersITS[i]->SetLineColor(color); fhNClustersITS[i]->SetLineWidth(2); fhNClustersTPC[i]->SetLineColor(color); fhNClustersTPC[i]->SetLineWidth(2); @@ -594,39 +898,125 @@ AliESDtrackCuts::CountAcceptedTracks(AliESD* esd) fhDZ[i]->SetLineColor(color); fhDZ[i]->SetLineWidth(2); fhDXYNormalized[i]->SetLineColor(color); fhDXYNormalized[i]->SetLineWidth(2); - fhDZNormalized[i]->SetLineColor(color); fhDZNormalized[i]->SetLineWidth(2); + fhDZNormalized[i]->SetLineColor(color); fhDZNormalized[i]->SetLineWidth(2); + fhNSigmaToVertex[i]->SetLineColor(color); fhNSigmaToVertex[i]->SetLineWidth(2); } + + // The number of sigmas to the vertex is per definition gaussian + ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50); + ffDTheoretical->SetParameter(0,1); + + TH1::AddDirectory(oldStatus); } //____________________________________________________________________ -void -AliESDtrackCuts::Print(const Option_t*) const { +Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir) +{ // - // print method - still to be implemented + // loads the histograms from a file + // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) // - AliInfo("AliESDtrackCuts..."); -} + if (!dir) + dir = GetName(); + + if (!gDirectory->cd(dir)) + return kFALSE; + ffDTheoretical = dynamic_cast (gDirectory->Get("nSigmaToVertexTheory")); + + fhCutStatistics = dynamic_cast (gDirectory->Get("cut_statistics")); + fhCutCorrelation = dynamic_cast (gDirectory->Get("cut_correlation")); + + Char_t str[5]; + for (Int_t i=0; i<2; i++) { + if (i==0) + { + gDirectory->cd("before_cuts"); + str[0] = 0; + } + else + { + gDirectory->cd("after_cuts"); + sprintf(str,"_cut"); + } + + fhNClustersITS[i] = dynamic_cast (gDirectory->Get(Form("nClustersITS%s",str) )); + fhNClustersTPC[i] = dynamic_cast (gDirectory->Get(Form("nClustersTPC%s",str) )); + fhChi2PerClusterITS[i] = dynamic_cast (gDirectory->Get(Form("chi2PerClusterITS%s",str))); + fhChi2PerClusterTPC[i] = dynamic_cast (gDirectory->Get(Form("chi2PerClusterTPC%s",str))); + + fhC11[i] = dynamic_cast (gDirectory->Get(Form("covMatrixDiagonal11%s",str))); + fhC22[i] = dynamic_cast (gDirectory->Get(Form("covMatrixDiagonal22%s",str))); + fhC33[i] = dynamic_cast (gDirectory->Get(Form("covMatrixDiagonal33%s",str))); + fhC44[i] = dynamic_cast (gDirectory->Get(Form("covMatrixDiagonal44%s",str))); + fhC55[i] = dynamic_cast (gDirectory->Get(Form("covMatrixDiagonal55%s",str))); + + fhDXY[i] = dynamic_cast (gDirectory->Get(Form("dXY%s",str) )); + fhDZ[i] = dynamic_cast (gDirectory->Get(Form("dZ%s",str) )); + fhDXYvsDZ[i] = dynamic_cast (gDirectory->Get(Form("dXYvsDZ%s",str))); + + fhDXYNormalized[i] = dynamic_cast (gDirectory->Get(Form("dXYNormalized%s",str) )); + fhDZNormalized[i] = dynamic_cast (gDirectory->Get(Form("dZNormalized%s",str) )); + fhDXYvsDZNormalized[i] = dynamic_cast (gDirectory->Get(Form("dXYvsDZNormalized%s",str))); + + fhNSigmaToVertex[i] = dynamic_cast (gDirectory->Get(Form("nSigmaToVertex%s",str))); + + // TODO only temporary + /*fhNClustersITS[i]->SetTitle("n ITS clusters"); + fhNClustersTPC[i]->SetTitle("n TPC clusters"); + fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster"); + fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster"); + + fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]"); + fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]"); + fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}"); + fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}"); + fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]"); + + fhDXY[i]->SetTitle("transverse impact parameter"); + fhDZ[i]->SetTitle("longitudinal impact parameter"); + fhDXYvsDZ[i]->SetTitle("longitudinal impact parameter"); + fhDXYvsDZ[i]->SetYTitle("transverse impact parameter"); + + fhDXYNormalized[i]->SetTitle("normalized trans impact par"); + fhDZNormalized[i]->SetTitle("normalized long impact par"); + fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par"); + fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par"); + fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");*/ + + gDirectory->cd("../"); + } + + gDirectory->cd(".."); + + return kTRUE; +} //____________________________________________________________________ -void AliESDtrackCuts::SaveHistograms(Char_t* dir) { - // +void AliESDtrackCuts::SaveHistograms(const Char_t* dir) { + // // saves the histograms in a directory (dir) - // + // - if (!fHistogramsOn) { AliDebug(0, "Histograms not on - cannot save histograms!!!"); return; } + if (!dir) + dir = GetName(); + gDirectory->mkdir(dir); gDirectory->cd(dir); gDirectory->mkdir("before_cuts"); gDirectory->mkdir("after_cuts"); + // a factor of 2 is needed since n sigma is positive + ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width")); + ffDTheoretical->Write("nSigmaToVertexTheory"); + fhCutStatistics->Write(); fhCutCorrelation->Write(); @@ -635,12 +1025,12 @@ void AliESDtrackCuts::SaveHistograms(Char_t* dir) { gDirectory->cd("before_cuts"); else gDirectory->cd("after_cuts"); - + fhNClustersITS[i] ->Write(); fhNClustersTPC[i] ->Write(); fhChi2PerClusterITS[i] ->Write(); fhChi2PerClusterTPC[i] ->Write(); - + fhC11[i] ->Write(); fhC22[i] ->Write(); fhC33[i] ->Write(); @@ -650,10 +1040,11 @@ void AliESDtrackCuts::SaveHistograms(Char_t* dir) { fhDXY[i] ->Write(); fhDZ[i] ->Write(); fhDXYvsDZ[i] ->Write(); - + fhDXYNormalized[i] ->Write(); fhDZNormalized[i] ->Write(); fhDXYvsDZNormalized[i] ->Write(); + fhNSigmaToVertex[i] ->Write(); gDirectory->cd("../"); } @@ -661,5 +1052,132 @@ void AliESDtrackCuts::SaveHistograms(Char_t* dir) { gDirectory->cd("../"); } - +//____________________________________________________________________ +void AliESDtrackCuts::DrawHistograms() +{ + // draws some histograms + + TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800); + canvas1->Divide(2, 2); + + canvas1->cd(1); + fhNClustersTPC[0]->SetStats(kFALSE); + fhNClustersTPC[0]->Draw(); + + canvas1->cd(2); + fhChi2PerClusterTPC[0]->SetStats(kFALSE); + fhChi2PerClusterTPC[0]->Draw(); + + canvas1->cd(3); + fhNSigmaToVertex[0]->SetStats(kFALSE); + fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10); + fhNSigmaToVertex[0]->Draw(); + + canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName())); + + TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800); + canvas2->Divide(3, 2); + + canvas2->cd(1); + fhC11[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhC11[0]->Draw(); + + canvas2->cd(2); + fhC22[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhC22[0]->Draw(); + + canvas2->cd(3); + fhC33[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhC33[0]->Draw(); + + canvas2->cd(4); + fhC44[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhC44[0]->Draw(); + + canvas2->cd(5); + fhC55[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhC55[0]->Draw(); + + canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName())); + + TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800); + canvas3->Divide(3, 2); + + canvas3->cd(1); + fhDXY[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhDXY[0]->Draw(); + + canvas3->cd(2); + fhDZ[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhDZ[0]->Draw(); + + canvas3->cd(3); + fhDXYvsDZ[0]->SetStats(kFALSE); + gPad->SetLogz(); + gPad->SetRightMargin(0.15); + fhDXYvsDZ[0]->Draw("COLZ"); + + canvas3->cd(4); + fhDXYNormalized[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhDXYNormalized[0]->Draw(); + + canvas3->cd(5); + fhDZNormalized[0]->SetStats(kFALSE); + gPad->SetLogy(); + fhDZNormalized[0]->Draw(); + + canvas3->cd(6); + fhDXYvsDZNormalized[0]->SetStats(kFALSE); + gPad->SetLogz(); + gPad->SetRightMargin(0.15); + fhDXYvsDZNormalized[0]->Draw("COLZ"); + + canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName())); + + TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500); + canvas4->Divide(2, 1); + + canvas4->cd(1); + fhCutStatistics->SetStats(kFALSE); + fhCutStatistics->LabelsOption("v"); + gPad->SetBottomMargin(0.3); + fhCutStatistics->Draw(); + + canvas4->cd(2); + fhCutCorrelation->SetStats(kFALSE); + fhCutCorrelation->LabelsOption("v"); + gPad->SetBottomMargin(0.3); + gPad->SetLeftMargin(0.3); + fhCutCorrelation->Draw("COLZ"); + + canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName())); + + /*canvas->cd(1); + fhDXYvsDZNormalized[0]->SetStats(kFALSE); + fhDXYvsDZNormalized[0]->DrawCopy("COLZ"); + + canvas->cd(2); + fhNClustersTPC[0]->SetStats(kFALSE); + fhNClustersTPC[0]->DrawCopy(); + + canvas->cd(3); + fhChi2PerClusterITS[0]->SetStats(kFALSE); + fhChi2PerClusterITS[0]->DrawCopy(); + fhChi2PerClusterITS[1]->SetLineColor(2); + fhChi2PerClusterITS[1]->DrawCopy("SAME"); + + canvas->cd(4); + fhChi2PerClusterTPC[0]->SetStats(kFALSE); + fhChi2PerClusterTPC[0]->DrawCopy(); + fhChi2PerClusterTPC[1]->SetLineColor(2); + fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/ +}