X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSPIDv1.cxx;h=5816f07b7aab3d59bf4bfa833f4c29774948c647;hb=21d4a8d15d73f222ad4aba7cffa316925e0c6e60;hp=0b8d387a3a4b1f728601e82e84e203a6cf48bb30;hpb=35adb63890872b28c22d07d1d528fd702e1b0d35;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSPIDv1.cxx b/PHOS/AliPHOSPIDv1.cxx index 0b8d387a3a4..5816f07b7aa 100644 --- a/PHOS/AliPHOSPIDv1.cxx +++ b/PHOS/AliPHOSPIDv1.cxx @@ -15,6 +15,50 @@ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * Revision 1.113 2007/08/07 14:12:03 kharlov + * Quality assurance added (Yves Schutz) + * + * Revision 1.112 2007/07/11 13:43:30 hristov + * New class AliESDEvent, backward compatibility with the old AliESD (Christian) + * + * Revision 1.111 2007/05/04 14:49:29 policheh + * AliPHOSRecPoint inheritance from AliCluster + * + * Revision 1.110 2007/04/24 10:08:03 kharlov + * Vertex extraction from GenHeader + * + * Revision 1.109 2007/04/18 09:34:05 kharlov + * Geometry bug fixes + * + * Revision 1.108 2007/04/16 09:03:37 kharlov + * Incedent angle correction fixed + * + * Revision 1.107 2007/04/02 15:00:16 cvetan + * No more calls to gAlice in the reconstruction + * + * Revision 1.106 2007/04/01 15:40:15 kharlov + * Correction for actual vertex position implemented + * + * Revision 1.105 2007/03/06 06:57:46 kharlov + * DP:calculation of distance to CPV done in TSM + * + * Revision 1.104 2006/12/15 10:46:26 hristov + * Using TMath::Abs instead of fabs + * + * Revision 1.103 2006/09/07 18:31:08 kharlov + * Effective c++ corrections (T.Pocheptsov) + * + * Revision 1.102 2006/01/23 17:51:48 hristov + * Using the recommended way of forward declarations for TVector and TMatrix (see v5-08-00 release notes). Additional clean-up + * + * Revision 1.101 2005/05/28 14:19:04 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ // Implementation version v1 of the PHOS particle identifier // Particle identification based on the @@ -79,39 +123,57 @@ // // // --- ROOT system --- -#include "TROOT.h" -#include "TTree.h" -#include "TFile.h" -#include "TF2.h" + + +// --- Standard library --- +#include #include "TFormula.h" -#include "TCanvas.h" -#include "TFolder.h" -#include "TSystem.h" #include "TBenchmark.h" -#include "TMatrixD.h" #include "TPrincipal.h" +#include "TFile.h" #include "TSystem.h" - -// --- Standard library --- - +#include "TVector3.h" // --- AliRoot header files --- - -#include "AliGenerator.h" + //#include "AliLog.h" #include "AliPHOS.h" #include "AliPHOSPIDv1.h" -#include "AliPHOSClusterizerv1.h" -#include "AliPHOSEmcRecPoint.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" #include "AliPHOSTrackSegment.h" -#include "AliPHOSTrackSegmentMakerv1.h" +#include "AliPHOSEmcRecPoint.h" #include "AliPHOSRecParticle.h" -#include "AliPHOSGeometry.h" -#include "AliPHOSGetter.h" ClassImp( AliPHOSPIDv1) //____________________________________________________________________________ -AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() +AliPHOSPIDv1::AliPHOSPIDv1() : + AliPHOSPID(), + fBayesian(kFALSE), + fDefaultInit(kFALSE), + fWrite(kFALSE), + fFileNamePrincipalPhoton(), + fFileNamePrincipalPi0(), + fFileNameParameters(), + fPrincipalPhoton(0), + fPrincipalPi0(0), + fX(0), + fPPhoton(0), + fPPi0(0), + fParameters(0), + fVtx(0.), + fTFphoton(0), + fTFpiong(0), + fTFkaong(0), + fTFkaonl(0), + fTFhhadrong(0), + fTFhhadronl(0), + fDFmuon(0), + fERecWeight(0), + fChargedNeutralThreshold(0.), + fTOFEnThreshold(0), + fDispEnThreshold(0), + fDispMultThreshold(0) { // default ctor @@ -120,21 +182,73 @@ AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() } //____________________________________________________________________________ -AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid) +AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) : + AliPHOSPID(pid), + fBayesian(kFALSE), + fDefaultInit(kFALSE), + fWrite(kFALSE), + fFileNamePrincipalPhoton(), + fFileNamePrincipalPi0(), + fFileNameParameters(), + fPrincipalPhoton(0), + fPrincipalPi0(0), + fX(0), + fPPhoton(0), + fPPi0(0), + fParameters(0), + fVtx(0.), + fTFphoton(0), + fTFpiong(0), + fTFkaong(0), + fTFkaonl(0), + fTFhhadrong(0), + fTFhhadronl(0), + fDFmuon(0), + fERecWeight(0), + fChargedNeutralThreshold(0.), + fTOFEnThreshold(0), + fDispEnThreshold(0), + fDispMultThreshold(0) + { // ctor InitParameters() ; - Init() ; } //____________________________________________________________________________ -AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName) +AliPHOSPIDv1::AliPHOSPIDv1(AliPHOSGeometry *geom): + AliPHOSPID(geom), + fBayesian(kFALSE), + fDefaultInit(kFALSE), + fWrite(kFALSE), + fFileNamePrincipalPhoton(), + fFileNamePrincipalPi0(), + fFileNameParameters(), + fPrincipalPhoton(0), + fPrincipalPi0(0), + fX(0), + fPPhoton(0), + fPPi0(0), + fParameters(0), + fVtx(0.), + fTFphoton(0), + fTFpiong(0), + fTFkaong(0), + fTFkaonl(0), + fTFhhadrong(0), + fTFhhadronl(0), + fDFmuon(0), + fERecWeight(0), + fChargedNeutralThreshold(0.), + fTOFEnThreshold(0), + fDispEnThreshold(0), + fDispMultThreshold(0) + { //ctor with the indication on where to look for the track segments InitParameters() ; - Init() ; fDefaultInit = kFALSE ; } @@ -142,94 +256,101 @@ AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFold AliPHOSPIDv1::~AliPHOSPIDv1() { // dtor + fPrincipalPhoton = 0; + fPrincipalPi0 = 0; delete [] fX ; // Principal input delete [] fPPhoton ; // Photon Principal components delete [] fPPi0 ; // Pi0 Principal components -} -//____________________________________________________________________________ -const TString AliPHOSPIDv1::BranchName() const -{ - return GetName() ; + delete fParameters; + delete fTFphoton; + delete fTFpiong; + delete fTFkaong; + delete fTFkaonl; + delete fTFhhadrong; + delete fTFhhadronl; + delete fDFmuon; } -//____________________________________________________________________________ -void AliPHOSPIDv1::Init() -{ - // Make all memory allocations that are not possible in default constructor - // Add the PID task to the list of PHOS tasks - - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - if(!gime) - gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; - - if ( !gime->PID() ) - gime->PostPID(this) ; -} - //____________________________________________________________________________ void AliPHOSPIDv1::InitParameters() { // Initialize PID parameters fWrite = kTRUE ; - fRecParticlesInRun = 0 ; - fNEvent = 0 ; - fRecParticlesInRun = 0 ; fBayesian = kTRUE ; SetParameters() ; // fill the parameters matrix from parameters file - SetEventRange(0,-1) ; // initialisation of response function parameters // Tof + +// // Photons +// fTphoton[0] = 0.218 ; +// fTphoton[1] = 1.55E-8 ; +// fTphoton[2] = 5.05E-10 ; +// fTFphoton = new TFormula("ToF response to photons" , "gaus") ; +// fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; + +// // Pions +// //Gaus (0 to max probability) +// fTpiong[0] = 0.0971 ; +// fTpiong[1] = 1.58E-8 ; +// fTpiong[2] = 5.69E-10 ; +// fTFpiong = new TFormula("ToF response to pions" , "gaus") ; +// fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; + +// // Kaons +// //Gaus (0 to max probability) +// fTkaong[0] = 0.0542 ; +// fTkaong[1] = 1.64E-8 ; +// fTkaong[2] = 6.07E-10 ; +// fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; +// fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; +// //Landau (max probability to inf) +// fTkaonl[0] = 0.264 ; +// fTkaonl[1] = 1.68E-8 ; +// fTkaonl[2] = 4.10E-10 ; +// fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; +// fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; + +// //Heavy Hadrons +// //Gaus (0 to max probability) +// fThhadrong[0] = 0.0302 ; +// fThhadrong[1] = 1.73E-8 ; +// fThhadrong[2] = 9.52E-10 ; +// fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; +// fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; +// //Landau (max probability to inf) +// fThhadronl[0] = 0.139 ; +// fThhadronl[1] = 1.745E-8 ; +// fThhadronl[2] = 1.00E-9 ; +// fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; +// fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; + // Photons - fTphoton[0] = 0.218 ; - //fTphoton[0] = 1. ; + fTphoton[0] = 7.83E8 ; fTphoton[1] = 1.55E-8 ; - fTphoton[2] = 5.05E-10 ; + fTphoton[2] = 5.09E-10 ; fTFphoton = new TFormula("ToF response to photons" , "gaus") ; fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; -// // Electrons -// fTelectron[0] = 0.2 ; -// fTelectron[1] = 1.55E-8 ; -// fTelectron[2] = 5.35E-10 ; -// fTFelectron = new TFormula("ToF response to electrons" , "gaus") ; -// fTFelectron->SetParameters( fTelectron[0], fTelectron[1], fTelectron[2]) ; -// // Muons -// fTmuon[0] = 0.2 ; -// fTmuon[1] = 1.55E-8 ; -// fTmuon[2] = 5.1E-10 ; -// fTFmuon = new TFormula("ToF response to muons" , "gaus") ; -// fTFmuon->SetParameters( fTmuon[0], fTmuon[1], fTmuon[2]) ; // Pions //Gaus (0 to max probability) - fTpiong[0] = 0.0971 ; - //fTpiong[0] = 1. ; + fTpiong[0] = 6.73E8 ; fTpiong[1] = 1.58E-8 ; - fTpiong[2] = 5.69E-10 ; + fTpiong[2] = 5.87E-10 ; fTFpiong = new TFormula("ToF response to pions" , "gaus") ; fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; - // Landau (max probability to inf) -// fTpionl[0] = 0.05 ; -// //fTpionl[0] = 5.53 ; -// fTpionl[1] = 1.68E-8 ; -// fTpionl[2] = 5.38E-10 ; -// fTFpionl = new TFormula("ToF response to pions" , "landau") ; -// fTFpionl->SetParameters( fTpionl[0], fTpionl[1], fTpionl[2]) ; - // Kaons //Gaus (0 to max probability) - fTkaong[0] = 0.0542 ; - //fTkaong[0] = 1. ; + fTkaong[0] = 3.93E8 ; fTkaong[1] = 1.64E-8 ; - fTkaong[2] = 6.07-10 ; + fTkaong[2] = 6.07E-10 ; fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; //Landau (max probability to inf) - fTkaonl[0] = 0.264 ; - //fTkaonl[0] = 5.53 ; + fTkaonl[0] = 2.0E9 ; fTkaonl[1] = 1.68E-8 ; fTkaonl[2] = 4.10E-10 ; fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; @@ -237,104 +358,132 @@ void AliPHOSPIDv1::InitParameters() //Heavy Hadrons //Gaus (0 to max probability) - fThhadrong[0] = 0.0302 ; - //fThhadrong[0] = 1. ; + fThhadrong[0] = 2.02E8 ; fThhadrong[1] = 1.73E-8 ; fThhadrong[2] = 9.52E-10 ; fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; //Landau (max probability to inf) - fThhadronl[0] = 0.139 ; - //fThhadronl[0] = 5.53 ; - fThhadronl[1] = 1.745E-8 ; - fThhadronl[2] = 1.00E-9 ; + fThhadronl[0] = 1.10E9 ; + fThhadronl[1] = 1.74E-8 ; + fThhadronl[2] = 1.00E-9 ; fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; -/// /gaussian parametrization for pions -// fTpion[0] = 3.93E-2 ; fTpion[1] = 0.130 ; fTpion[2] =-6.37E-2 ;//constant -// fTpion[3] = 1.65E-8 ; fTpion[4] =-1.40E-9 ; fTpion[5] = 5.96E-10;//mean -// fTpion[6] = 8.09E-10; fTpion[7] =-4.65E-10; fTpion[8] = 1.50E-10;//sigma - -// //landau parametrization for kaons -// fTkaon[0] = 0.107 ; fTkaon[1] = 0.166 ; fTkaon[2] = 0.243 ;//constant -// fTkaon[3] = 1.80E-8 ; fTkaon[4] =-2.96E-9 ; fTkaon[5] = 9.60E-10;//mean -// fTkaon[6] = 1.37E-9 ; fTkaon[7] =-1.80E-9 ; fTkaon[8] = 6.74E-10;//sigma - -// //landau parametrization for nucleons -// fThhadron[0] = 6.33E-2 ; fThhadron[1] = 2.52E-2 ; fThhadron[2] = 2.16E-2 ;//constant -// fThhadron[3] = 1.94E-8 ; fThhadron[4] =-7.06E-10; fThhadron[5] =-4.69E-10;//mean -// fThhadron[6] = 2.55E-9 ; fThhadron[7] =-1.90E-9 ; fThhadron[8] = 5.41E-10;//sigma // Shower shape: dispersion gaussian parameters // Photons - - fDphoton[0] = 0.1 ; fDphoton[1] = 0. ; fDphoton[2] = 0. ;//constant - //fDphoton[0] = 1.0 ; fDphoton[1] = 0. ; fDphoton[2] = 0. ;//constant - fDphoton[3] = 1.55 ; fDphoton[4] =-0.0863 ; fDphoton[5] = 0.287 ;//mean - fDphoton[6] = 0.0451 ; fDphoton[7] =-0.0803 ; fDphoton[8] = 0.314 ;//sigma - - fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant - //fDpi0[0] = 1.0 ; fDpi0[1] = 0.0 ; fDpi0[2] = 0. ;//constant - fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean - fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma - //landau -// fDhadron[0] = 0.007 ; fDhadron[1] = 0. ; fDhadron[2] = 0. ;//constant -// //fDhadron[0] = 5.53 ; fDhadron[1] = 0. ; fDhadron[2] = 0. ;//constant -// fDhadron[3] = 3.38 ; fDhadron[4] = 0.0833 ; fDhadron[5] =-0.845 ;//mean -// fDhadron[6] = 0.627 ; fDhadron[7] = 0.012 ; fDhadron[8] =-0.170 ;//sigma - - fDhadron[0] =-5.10E-3 ; fDhadron[1] =-5.35E-3 ; fDhadron[2] = 3.77E-2 ;//constant - fDhadron[3] = 4.03 ; fDhadron[4] = 0.292 ; fDhadron[5] =-1.50 ;//mean - fDhadron[6] = 0.958 ; fDhadron[7] = 0.117 ; fDhadron[8] =-0.598 ;//sigma - // Muons - fDmuon[0] = 0.0631 ; - fDmuon[1] = 1.4 ; + +// fDphoton[0] = 4.62e-2; fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant +// fDphoton[3] = 1.53 ; fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339 ;//mean +// fDphoton[6] = 6.89e-2; fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194 ;//sigma + +// fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant +// fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean +// fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma + +// fDhadron[0] = 1.61E-2 ; fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant +// fDhadron[3] = 3.81 ; fDhadron[4] = 0.232 ; fDhadron[5] =-1.25 ;//mean +// fDhadron[6] = 0.897 ; fDhadron[7] = 0.0987 ; fDhadron[8] =-0.534 ;//sigma + + fDphoton[0] = 1.5 ; fDphoton[1] = 0.49 ; fDphoton[2] =-1.7E-2 ;//constant + fDphoton[3] = 1.5 ; fDphoton[4] = 4.0E-2 ; fDphoton[5] = 0.21 ;//mean + fDphoton[6] = 4.8E-2 ; fDphoton[7] =-0.12 ; fDphoton[8] = 0.27 ;//sigma + fDphoton[9] = 16.; //for E> fDphoton[9] parameters calculated at fDphoton[9] + + fDpi0[0] = 0.25 ; fDpi0[1] = 3.3E-2 ; fDpi0[2] =-1.0e-5 ;//constant + fDpi0[3] = 1.50 ; fDpi0[4] = 398. ; fDpi0[5] = 12. ;//mean + fDpi0[6] =-7.0E-2 ; fDpi0[7] =-524. ; fDpi0[8] = 22. ;//sigma + fDpi0[9] = 110.; //for E> fDpi0[9] parameters calculated at fDpi0[9] + + fDhadron[0] = 6.5 ; fDhadron[1] =-5.3 ; fDhadron[2] = 1.5 ;//constant + fDhadron[3] = 3.8 ; fDhadron[4] = 0.23 ; fDhadron[5] =-1.2 ;//mean + fDhadron[6] = 0.88 ; fDhadron[7] = 9.3E-2 ; fDhadron[8] =-0.51 ;//sigma + fDhadron[9] = 2.; //for E> fDhadron[9] parameters calculated at fDhadron[9] + + fDmuon[0] = 0.0631 ; + fDmuon[1] = 1.4 ; fDmuon[2] = 0.0557 ; fDFmuon = new TFormula("Shower shape response to muons" , "landau") ; fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ; - // CPV-EMC distance gaussian parameters - fCPVelectron[0] = 0.0 ; fCPVelectron[1] = 0.0160 ; fCPVelectron[2] = 0. ;//constant - //fCPVelectron[0] = 1.0 ; fCPVelectron[1] = 0. ; fCPVelectron[2] = 0. ;//constant - fCPVelectron[3] = 0.0682 ; fCPVelectron[4] =-1.32 ; fCPVelectron[5] = 6.67 ;//mean - fCPVelectron[6] = 0.276 ; fCPVelectron[7] = 0.234 ; fCPVelectron[8] = 0.356 ;//sigma + // x(CPV-EMC) distance gaussian parameters + +// fXelectron[0] = 8.06e-2 ; fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant +// fXelectron[3] = 0.202 ; fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55 ;//mean +// fXelectron[6] = 0.334 ; fXelectron[7] = 0.186 ; fXelectron[8] = 4.32e-2;//sigma + +// //charged hadrons gaus +// fXcharged[0] = 6.43e-3 ; fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant +// fXcharged[3] = 2.75 ; fXcharged[4] =-0.40 ; fXcharged[5] = 1.68 ;//mean +// fXcharged[6] = 3.135 ; fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma + +// // z(CPV-EMC) distance gaussian parameters + +// fZelectron[0] = 8.22e-2 ; fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant +// fZelectron[3] = 3.09e-2 ; fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean +// fZelectron[6] = 0.263 ; fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma + +// //charged hadrons gaus + +// fZcharged[0] = 1.00e-2 ; fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant +// fZcharged[3] =-4.68e-2 ; fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean +// fZcharged[6] = 1.425 ; fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma - //all charged landau - // fCPVcharged[0] = 0.0 ; fCPVcharged[1] = 0.0464 ; fCPVcharged[2] = 0. ;//constant -// //fCPVcharged[0] = 5.53 ; fCPVcharged[1] = 0. ; fCPVcharged[2] = 0. ;//constant -// fCPVcharged[3] = 0.297 ; fCPVcharged[4] =-0.652 ; fCPVcharged[5] = 1.91 ;//mean -// fCPVcharged[6] = 0.0786 ; fCPVcharged[7] =-0.237 ; fCPVcharged[8] = 0.752 ;//sigma -// //charged hadrons landau -// fCPVchargedl[0] = 0.103 ; fCPVchargedl[1] = 8.84E-3 ; fCPVchargedl[2] =-2.40E-2 ;//constant -// fCPVchargedl[3] = 2.86 ; fCPVchargedl[4] =-0.214 ; fCPVchargedl[5] = 0.817 ;//mean -// fCPVchargedl[6] = 0.182 ; fCPVchargedl[7] =-0.0693 ; fCPVchargedl[8] = 0.319 ;//sigma -// //charged hadrons gaus -// fCPVchargedg[0] = 0.0135 ; fCPVchargedg[1] = 2.43E-5 ; fCPVchargedg[2] = 3.01E-3 ;//constant -// fCPVchargedg[3] = 2.37 ; fCPVchargedg[4] =-0.181 ; fCPVchargedg[5] = 0.726 ;//mean -// fCPVchargedg[6] = 0.908 ; fCPVchargedg[7] =-0.0558 ; fCPVchargedg[8] = 0.219 ;//sigma + fXelectron[0] =-1.6E-2 ; fXelectron[1] = 0.77 ; fXelectron[2] =-0.15 ;//constant + fXelectron[3] = 0.35 ; fXelectron[4] = 0.25 ; fXelectron[5] = 4.12 ;//mean + fXelectron[6] = 0.30 ; fXelectron[7] = 0.11 ; fXelectron[8] = 0.16 ;//sigma + fXelectron[9] = 3.; //for E> fXelectron[9] parameters calculated at fXelectron[9] + //charged hadrons gaus + fXcharged[0] = 0.14 ; fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0 ;//constant + fXcharged[3] = 1.4 ; fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4 ;//mean + fXcharged[6] = 5.7 ; fXcharged[7] = 0.27 ; fXcharged[8] =-1.8 ;//sigma + fXcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] - //charged hadrons landau - fCPVcharged[0] = 6.06E-2 ; fCPVcharged[1] = 3.80E-3 ; fCPVcharged[2] =-1.40E-2 ;//constant - fCPVcharged[3] = 1.15 ; fCPVcharged[4] =-0.563 ; fCPVcharged[5] = 2.63 ;//mean - fCPVcharged[6] = 0.915 ; fCPVcharged[7] =-0.0790 ; fCPVcharged[8] = 0.307 ;//sigma + // z(CPV-EMC) distance gaussian parameters + + fZelectron[0] = 0.49 ; fZelectron[1] = 0.53 ; fZelectron[2] =-9.8E-2 ;//constant + fZelectron[3] = 2.8E-2 ; fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;//mean + fZelectron[6] = 0.25 ; fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17 ;//sigma + fZelectron[9] = 3.; //for E> fZelectron[9] parameters calculated at fZelectron[9] - for (Int_t i =0; i< AliESDtrack::kSPECIESN ; i++) - fInitPID[i] = 1.; + //charged hadrons gaus + fZcharged[0] = 0.46 ; fZcharged[1] =-0.65 ; fZcharged[2] = 0.52 ;//constant + fZcharged[3] = 1.1E-2 ; fZcharged[4] = 0. ; fZcharged[5] = 0. ;//mean + fZcharged[6] = 0.60 ; fZcharged[7] =-8.2E-2 ; fZcharged[8] = 0.45 ;//sigma + fZcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] + + //Threshold to differentiate between charged and neutral + fChargedNeutralThreshold = 1e-5; + fTOFEnThreshold = 2; //Maximum energy to use TOF + fDispEnThreshold = 0.5; //Minimum energy to use shower shape + fDispMultThreshold = 3; //Minimum multiplicity to use shower shape + + //Weight to hadrons recontructed energy + + fERecWeightPar[0] = 0.32 ; + fERecWeightPar[1] = 3.8 ; + fERecWeightPar[2] = 5.4E-3 ; + fERecWeightPar[3] = 5.6E-2 ; + fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ; + fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ; + + + for (Int_t i =0; i< AliPID::kSPECIESN ; i++) + fInitPID[i] = 1.; + } //________________________________________________________________________ -void AliPHOSPIDv1::Exec(Option_t *option) +void AliPHOSPIDv1::TrackSegments2RecParticles(Option_t *option) { // Steering method to perform particle reconstruction and identification // for the event range from fFirstEvent to fLastEvent. - // This range is optionally set by SetEventRange(). - // if fLastEvent=-1 (by default), then process events until the end. if(strstr(option,"tim")) gBenchmark->Start("PHOSPID"); @@ -344,79 +493,74 @@ void AliPHOSPIDv1::Exec(Option_t *option) return ; } + if(fTrackSegments && //Skip events, where no track segments made + fTrackSegments->GetEntriesFast()) { - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - - if (fLastEvent == -1) - fLastEvent = gime->MaxEvent() - 1 ; - else - fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); - Int_t nEvents = fLastEvent - fFirstEvent + 1; - - Int_t ievent ; - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - gime->Event(ievent,"TR") ; - if(gime->TrackSegments() && //Skip events, where no track segments made - gime->TrackSegments()->GetEntriesFast()) { - MakeRecParticles() ; - - if(fBayesian) - MakePID() ; + GetVertex() ; + MakeRecParticles() ; + + if(fBayesian) + MakePID() ; - WriteRecParticles(); - if(strstr(option,"deb")) - PrintRecParticles(option) ; - //increment the total number of rec particles per run - fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; - } + if(strstr(option,"deb")) + PrintRecParticles(option) ; } + if(strstr(option,"deb")) PrintRecParticles(option); if(strstr(option,"tim")){ gBenchmark->Stop("PHOSPID"); - Info("Exec", "took %f seconds for PID %f seconds per event", - gBenchmark->GetCpuTime("PHOSPID"), - gBenchmark->GetCpuTime("PHOSPID")/nEvents) ; + AliInfo(Form("took %f seconds for PID", + gBenchmark->GetCpuTime("PHOSPID"))); } - if(fWrite) - Unload(); } //________________________________________________________________________ Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par) { - Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; + //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), + //this method returns a density probability of this parameter, given by a gaussian + //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b + //Float_t xorg = x; + if (x > par[9]) x = par[9]; + + //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; + Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; - //cout<<"c "<< cnt << " mean "< 30) +// cout<<"En_in = "<SetParameters(cnt,mean,sigma); - //cout<<"gaus value "<Eval(y)<Eval(y) ; - return arg; + if(TMath::Abs(sigma) > 1.e-10){ + return cnt*TMath::Gaus(y,mean,sigma); } else return 0.; - + } //________________________________________________________________________ Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par) { + //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), + //this method returns a density probability of this parameter, given by a gaussian + //function whose parameters depend with the energy like second order polinomial + Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; Double_t mean = par[3] + par[4] * x + par[5] * x * x ; Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; - if(mean != 0. && sigma/mean > 1.e-4 ){ - TF1 * f = new TF1("gaus","gaus",0,100); - f->SetParameters(cnt,mean,sigma); - Double_t arg = f->Eval(y) ; - return arg; + if(TMath::Abs(sigma) > 1.e-10){ + return cnt*TMath::Gaus(y,mean,sigma); } else return 0.; + + + } //____________________________________________________________________________ @@ -425,9 +569,13 @@ const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const //Get file name that contains the PCA for a particle ("photon or pi0") particle.ToLower(); TString name; - if (particle=="photon") name = fFileNamePrincipalPhoton ; - else if (particle=="pi0" ) name = fFileNamePrincipalPi0 ; - else Error("GetFileNamePrincipal","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); + if (particle=="photon") + name = fFileNamePrincipalPhoton ; + else if (particle=="pi0" ) + name = fFileNamePrincipalPi0 ; + else + AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", + particle.Data())); return name; } @@ -436,9 +584,9 @@ Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const { // Get the i-th parameter "Calibration" Float_t param = 0.; - if (i>2 || i<0) - Error("GetParameterCalibration","Invalid parameter number: %d",i); - else + if (i>2 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else param = (*fParameters)(0,i); return param; } @@ -464,13 +612,17 @@ Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const { // Get the i-th parameter "CPV-EMC distance" for the specified axis Float_t param = 0.; - if(i>2 || i<0) - Error("GetParameterCpv2Emc","Invalid parameter number: %d",i); - else { + if(i>2 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else { axis.ToLower(); - if (axis == "x") param = (*fParameters)(1,i); - else if (axis == "z") param = (*fParameters)(2,i); - else Error("GetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); + if (axis == "x") + param = (*fParameters)(1,i); + else if (axis == "z") + param = (*fParameters)(2,i); + else { + AliError(Form("Invalid axis name: %s",axis.Data())); + } } return param; } @@ -518,9 +670,9 @@ Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const // Get the parameter "i" to calculate the boundary on the moment M2x // for photons at high p_T Float_t param = 0; - if (i>3 || i<0) - Error("GetParameterPhotonBoundary","Wrong parameter number: %d\n",i); - else + if (i>3 || i<0) { + AliError(Form("Wrong parameter number: %d\n",i)); + } else param = (*fParameters)(14,i) ; return param; } @@ -531,9 +683,9 @@ Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const // Get the parameter "i" to calculate the boundary on the moment M2x // for pi0 at high p_T Float_t param = 0; - if (i>2 || i<0) - Error("GetParameterPi0Boundary","Wrong parameter number: %d\n",i); - else + if (i>2 || i<0) { + AliError(Form("Wrong parameter number: %d\n",i)); + } else param = (*fParameters)(15,i) ; return param; } @@ -544,9 +696,9 @@ Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const // Get TimeGate parameter depending on Purity-Efficiency i: // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity Float_t param = 0.; - if(i>2 || i<0) - Error("GetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); - else + if(i>2 || i<0) { + AliError(Form("Invalid Efficiency-Purity choice %d",i)); + } else param = (*fParameters)(3,i) ; return param; } @@ -560,10 +712,13 @@ Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString particle.ToLower(); param. ToLower(); Int_t offset = -1; - if (particle == "photon") offset=0; - else if (particle == "pi0") offset=5; + if (particle == "photon") + offset=0; + else if (particle == "pi0") + offset=5; else - Error("GetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); + AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", + particle.Data())); Int_t p= -1; Float_t par = 0; @@ -574,58 +729,65 @@ Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString else if(param.Contains("x0"))p=7+offset; else if(param.Contains("y0"))p=8+offset; - if (i>4 || i<0) - Error("GetParameterToCalculateEllipse", "No parameter with index", i) ; - else if (p==-1) - Error("GetParameterToCalculateEllipse", "No parameter with name %s", param.Data() ) ; - else + if (i>4 || i<0) { + AliError(Form("No parameter with index %d", i)) ; + } else if (p==-1) { + AliError(Form("No parameter with name %s", param.Data() )) ; + } else par = (*fParameters)(p,i) ; return par; } +//DP____________________________________________________________________________ +//Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const +//{ +// // Calculates the distance between the EMC RecPoint and the PPSD RecPoint +// +// const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; +// TVector3 vecEmc ; +// TVector3 vecCpv ; +// if(cpv){ +// emc->GetLocalPosition(vecEmc) ; +// cpv->GetLocalPosition(vecCpv) ; +// +// if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ +// // Correct to difference in CPV and EMC position due to different distance to center. +// // we assume, that particle moves from center +// Float_t dCPV = geom->GetIPtoOuterCoverDistance(); +// Float_t dEMC = geom->GetIPtoCrystalSurface() ; +// dEMC = dEMC / dCPV ; +// vecCpv = dEMC * vecCpv - vecEmc ; +// if (axis == "X") return vecCpv.X(); +// if (axis == "Y") return vecCpv.Y(); +// if (axis == "Z") return vecCpv.Z(); +// if (axis == "R") return vecCpv.Mag(); +// } +// return 100000000 ; +// } +// return 100000000 ; +//} //____________________________________________________________________________ -Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const -{ - // Calculates the distance between the EMC RecPoint and the PPSD RecPoint - - const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; - TVector3 vecEmc ; - TVector3 vecCpv ; - if(cpv){ - emc->GetLocalPosition(vecEmc) ; - cpv->GetLocalPosition(vecCpv) ; - - if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ - // Correct to difference in CPV and EMC position due to different distance to center. - // we assume, that particle moves from center - Float_t dCPV = geom->GetIPtoOuterCoverDistance(); - Float_t dEMC = geom->GetIPtoCrystalSurface() ; - dEMC = dEMC / dCPV ; - vecCpv = dEMC * vecCpv - vecEmc ; - if (axis == "X") return vecCpv.X(); - if (axis == "Y") return vecCpv.Y(); - if (axis == "Z") return vecCpv.Z(); - if (axis == "R") return vecCpv.Mag(); - } - return 100000000 ; - } - return 100000000 ; -} -//____________________________________________________________________________ -Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Int_t effPur, Float_t e) const +Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const { + //Calculates the pid bit for the CPV selection per each purity. if(effPur>2 || effPur<0) - Error("GetCPVBit","Invalid Efficiency-Purity choice %d",effPur); + AliError(Form("Invalid Efficiency-Purity choice %d",effPur)); + +//DP if(ts->GetCpvIndex()<0) +//DP return 1 ; //no CPV cluster Float_t sigX = GetCpv2EmcDistanceCut("X",e); Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); - Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); - Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); - //Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ ) ; - if((deltaX>sigX*(effPur+1))&&(deltaZ>sigZ*(effPur+1))) + Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X")); + Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z")); +// Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ; + + //if(deltaX>sigX*(effPur+1)) + //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1))) + if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1))) return 1;//Neutral else return 0;//Charged @@ -638,9 +800,9 @@ Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t particle.ToLower(); Int_t prinbit = 0 ; - Float_t a = GetEllipseParameter(particle,"a" , e); - Float_t b = GetEllipseParameter(particle,"b" , e); - Float_t c = GetEllipseParameter(particle,"c" , e); + Float_t a = GetEllipseParameter(particle,"a" , e); + Float_t b = GetEllipseParameter(particle,"b" , e); + Float_t c = GetEllipseParameter(particle,"c" , e); Float_t x0 = GetEllipseParameter(particle,"x0", e); Float_t y0 = GetEllipseParameter(particle,"y0", e); @@ -653,7 +815,7 @@ Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t if((effPur==0) && (r<9./2.)) prinbit= 1; if(r<0) - Error("GetPrincipalBit", "Negative square?") ; + AliError("Negative square?") ; return prinbit; @@ -671,7 +833,8 @@ Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ TMath::Power(GetParameterPhotonBoundary(2),2)) + GetParameterPhotonBoundary(3); - //Info("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); + AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f", + e,m2x,m2xBoundary)); if (m2x < m2xBoundary) return 1;// A hard photon else @@ -689,7 +852,7 @@ Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const Float_t m2x = emc->GetM2x(); Float_t m2xBoundary = GetParameterPi0Boundary(0) + e * GetParameterPi0Boundary(1); - //Info("GetHardPi0Bit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); + AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary)); if (m2x > m2xBoundary) return 1;// A hard pi0 else @@ -705,26 +868,48 @@ TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpv // However because of the poor position resolution of PPSD the direction is always taken as if we were // in case 1. - TVector3 dir(0,0,0) ; - - TVector3 emcglobalpos ; - TMatrix dummy ; - - emc->GetGlobalPosition(emcglobalpos, dummy) ; - - - dir = emcglobalpos ; - dir.SetMag(1.) ; + TVector3 local ; + emc->GetLocalPosition(local) ; - //account correction to the position of IP - Float_t xo,yo,zo ; //Coordinates of the origin - if(gAlice && gAlice->GetMCApp() && gAlice->Generator()) - gAlice->Generator()->GetOrigin(xo,yo,zo) ; + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + //Correct for the non-perpendicular incidence + // Correction for the depth of the shower starting point (TDR p 127) + Float_t para = 0.925 ; + Float_t parb = 6.52 ; + + //Remove Old correction (vertex at 0,0,0) + TVector3 vtxOld(0.,0.,0.) ; + TVector3 vInc ; + Float_t x=local.X() ; + Float_t z=local.Z() ; + phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ; + Float_t depthxOld = 0.; + Float_t depthzOld = 0.; + Float_t energy = emc->GetEnergy() ; + if (energy > 0 && vInc.Y()!=0.) { + depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; + depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; + } else{ - xo=yo=zo=0.; + AliError("Cluster with zero energy \n"); } - TVector3 origin(xo,yo,zo); - dir = dir - origin ; + //Apply Real vertex + phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ; + Float_t depthx = 0.; + Float_t depthz = 0.; + if (energy > 0 && vInc.Y()!=0.) { + depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; + depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; + } + + //Correct for the vertex position and shower depth + Double_t xd=x+(depthxOld-depthx) ; + Double_t zd=z+(depthzOld-depthz) ; + TVector3 dir(0,0,0) ; + phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ; + + dir-=fVtx ; + dir.SetMag(1.) ; return dir ; } @@ -732,19 +917,19 @@ TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpv //________________________________________________________________________ Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) { - Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; + //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), + //this method returns a density probability of this parameter, given by a landau + //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b + + if (x > par[9]) x = par[9]; + + //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; + Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; - // Double_t mean = par[3] + par[4] * x + par[5] * x * x ; -// Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; - - //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ; - //return cnt * TMath::Exp(arg) ; - if(mean != 0. && sigma/mean > 1.e-4 ){ - TF1 * f = new TF1("landau","landau",0.,100.); - f->SetParameters(cnt,mean,sigma); - Double_t arg = f->Eval(y) ; - return arg; + + if(TMath::Abs(sigma) > 1.e-10){ + return cnt*TMath::Landau(y,mean,sigma); } else return 0.; @@ -753,18 +938,22 @@ Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) //________________________________________________________________________ Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) { + + //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), + //this method returns a density probability of this parameter, given by a landau + //function whose parameters depend with the energy like second order polinomial + Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; - Double_t mean = par[4] * (x*x) + par[5] * x + par[3] ; - Double_t sigma = par[7] * (x*x) + par[8] * x + par[6] ; - - if(mean != 0. && sigma/mean > 1.e-4 ){ - TF1 * f = new TF1("landau","landau",0.,100.); - f->SetParameters(cnt,mean,sigma); - Double_t arg = f->Eval(y) ; - return arg; + Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ; + Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ; + + if(TMath::Abs(sigma) > 1.e-10){ + return cnt*TMath::Landau(y,mean,sigma); } else return 0.; + + } // //________________________________________________________________________ // Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl) @@ -802,255 +991,338 @@ Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) void AliPHOSPIDv1::MakePID() { // construct the PID weight from a Bayesian Method + + const Int_t kSPECIES = AliPID::kSPECIESN ; + + Int_t nparticles = fRecParticles->GetEntriesFast() ; - Int_t index ; - const Int_t kSPECIES = AliESDtrack::kSPECIESN ; - Int_t nparticles = AliPHOSGetter::Instance()->RecParticles()->GetEntriesFast() ; - -// const Int_t kMAXPARTICLES = 2000 ; -// if (nparticles >= kMAXPARTICLES) -// Error("MakePID", "Change size of MAXPARTICLES") ; -// Double_t stof[kSPECIES][kMAXPARTICLES] ; + if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) { + AliFatal("RecPoints or TrackSegments not found !") ; + } + TIter next(fTrackSegments) ; + AliPHOSTrackSegment * ts ; + Int_t index = 0 ; Double_t * stof[kSPECIES] ; Double_t * sdp [kSPECIES] ; Double_t * scpv[kSPECIES] ; - + Double_t * sw [kSPECIES] ; //Info("MakePID","Begin MakePID"); for (Int_t i =0; i< kSPECIES; i++){ stof[i] = new Double_t[nparticles] ; sdp [i] = new Double_t[nparticles] ; scpv[i] = new Double_t[nparticles] ; + sw [i] = new Double_t[nparticles] ; } - // make the normalized distribution of pid for this event - // w(pid) in the Bayesian formulation - for(index = 0 ; index < nparticles ; index ++) { - - AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; - AliPHOSEmcRecPoint * emc = AliPHOSGetter::Instance()->EmcRecPoint(index) ; - AliPHOSCpvRecPoint * cpv = AliPHOSGetter::Instance()->CpvRecPoint(index) ; + + while ( (ts = (AliPHOSTrackSegment *)next()) ) { - Float_t en = emc->GetEnergy(); + //cout<<">>>>>> Bayesian Index "<GetEmcIndex()>=0) + emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; - // Tof - Double_t time = recpar->ToF() ; - //cout<<">>>>>>>Energy "<GetCpvIndex()>=0) +// cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; +// +//// Int_t track = 0 ; +//// track = ts->GetTrackIndex() ; //TPC tracks ? + if (!emc) { + AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; + } + + + // ############Tof############################# + + // Info("MakePID", "TOF"); + Float_t en = emc->GetEnergy(); + Double_t time = emc->GetTime() ; + // cout<<">>>>>>>Energy "<Eval(time) ; - // stof[AliESDtrack::kElectron][index] = stof[AliESDtrack::kPhoton][index] ; - // if(time < fTpionl[1]) - // stof[AliESDtrack::kPion][index] = fTFpiong ->Eval(time) ; //gaus distribution - // else - // stof[AliESDtrack::kPion][index] = fTFpionl ->Eval(time) ; //landau distribution - // if(time < fTkaonl[1]) - // stof[AliESDtrack::kKaon][index] = fTFkaong ->Eval(time) ; //gaus distribution - // else - // stof[AliESDtrack::kKaon][index] = fTFkaonl ->Eval(time) ; //landau distribution - // if(time < fThhadronl[1]) - // stof[AliESDtrack::kProton][index] = fTFhhadrong ->Eval(time) ; //gaus distribution - // else - // stof[AliESDtrack::kProton][index] = fTFhhadronl ->Eval(time) ; //landau distribution - // stof[AliESDtrack::kNeutron][index] = stof[AliESDtrack::kProton][index] ; - // stof[AliESDtrack::kEleCon][index] = stof[AliESDtrack::kPhoton][index] ; - // // a conversion electron has the photon ToF - // stof[AliESDtrack::kKaon0][index] = stof[AliESDtrack::kKaon][index] ; - // stof[AliESDtrack::kMuon][index] = stof[AliESDtrack::kPhoton][index] ; - if(en < 2.) { - // cout<<"TOF: pi "<< GausPol2(en, time, fTpion)<Eval(time) ; - stof[AliESDtrack::kElectron][index] = stof[AliESDtrack::kPhoton][index] ; -// stof[AliESDtrack::kPion][index] = GausPol2(en, time, fTpion) ; //gaus distribution -// stof[AliESDtrack::kKaon][index] = LandauPol2(en, time, fTkaon) ; //gaus distribution -// stof[AliESDtrack::kProton][index] = LandauPol2(en, time, fThhadron); //gaus distribution - stof[AliESDtrack::kPion][index] = fTFpiong ->Eval(time) ; //landau distribution + stof[AliPID::kPhoton][index] = 1.; + stof[AliPID::kElectron][index] = 1.; + stof[AliPID::kEleCon][index] = 1.; + //We assing the same prob to charged hadrons, sum is 1 + stof[AliPID::kPion][index] = 1./3.; + stof[AliPID::kKaon][index] = 1./3.; + stof[AliPID::kProton][index] = 1./3.; + //We assing the same prob to neutral hadrons, sum is 1 + stof[AliPID::kNeutron][index] = 1./2.; + stof[AliPID::kKaon0][index] = 1./2.; + stof[AliPID::kMuon][index] = 1.; + + if(en < fTOFEnThreshold) { + + Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution + Double_t pTofKaon = 0; + if(time < fTkaonl[1]) - stof[AliESDtrack::kKaon][index] = fTFkaong ->Eval(time) ; //gaus distribution + pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution else - stof[AliESDtrack::kKaon][index] = fTFkaonl ->Eval(time) ; //landau distribution + pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution + + Double_t pTofNucleon = 0; + if(time < fThhadronl[1]) - stof[AliESDtrack::kProton][index] = fTFhhadrong ->Eval(time) ; //gaus distribution + pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution else - stof[AliESDtrack::kProton][index] = fTFhhadronl ->Eval(time) ; //landau distribution + pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution + //We assing the same prob to neutral hadrons, sum is the average prob + Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ; + //We assing the same prob to charged hadrons, sum is the average prob + Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ; + + stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ; + //gaus distribution + stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ; + //a conversion electron has the photon ToF + stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ; + + stof[AliPID::kElectron][index] = pTofPion ; + + stof[AliPID::kPion][index] = pTofChHadron ; + stof[AliPID::kKaon][index] = pTofChHadron ; + stof[AliPID::kProton][index] = pTofChHadron ; - stof[AliESDtrack::kNeutron][index] = stof[AliESDtrack::kProton][index] ; - stof[AliESDtrack::kEleCon][index] = stof[AliESDtrack::kPhoton][index] ; - // a conversion electron has the photon ToF - stof[AliESDtrack::kKaon0][index] = stof[AliESDtrack::kKaon][index] ; - stof[AliESDtrack::kMuon][index] = stof[AliESDtrack::kPhoton][index] ; + stof[AliPID::kKaon0][index] = pTofNeHadron ; + stof[AliPID::kNeutron][index] = pTofNeHadron ; } - else { - stof[AliESDtrack::kPhoton][index] = 1.; - stof[AliESDtrack::kElectron][index] = 1.; - stof[AliESDtrack::kPion][index] = 1.; - stof[AliESDtrack::kKaon][index] = 1.; - stof[AliESDtrack::kProton][index] = 1.; - stof[AliESDtrack::kNeutron][index] = 1.; - stof[AliESDtrack::kEleCon][index] = 1.; - stof[AliESDtrack::kKaon0][index] = 1.; - stof[AliESDtrack::kMuon][index] = 1.; - } - // Info("MakePID", "TOF passed"); - // Shower shape: Dispersion + // Info("MakePID", "Dispersion"); + + // ###########Shower shape: Dispersion#################### Float_t dispersion = emc->GetDispersion(); + //DP: Correct for non-perpendicular incidence + //DP: still to be done + //dispersion is not well defined if the cluster is only in few crystals - // Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), - // dispersion); - // Info("MakePID","ss: photon %f, hadron %f ", GausF (en , dispersion, fDphoton), - // LandauF(en , dispersion, fDhadron ) ); - if(emc->GetMultiplicity() > 4){ - sdp[AliESDtrack::kPhoton][index] = GausF (en , dispersion, fDphoton) ; - sdp[AliESDtrack::kElectron][index] = sdp[AliESDtrack::kPhoton][index] ; - sdp[AliESDtrack::kPion][index] = LandauF(en , dispersion, fDhadron ) ; - sdp[AliESDtrack::kKaon][index] = sdp[AliESDtrack::kPion][index] ; - sdp[AliESDtrack::kProton][index] = sdp[AliESDtrack::kPion][index] ; - sdp[AliESDtrack::kNeutron][index] = sdp[AliESDtrack::kPion][index] ; - sdp[AliESDtrack::kEleCon][index] = sdp[AliESDtrack::kPhoton][index]; - sdp[AliESDtrack::kKaon0][index] = sdp[AliESDtrack::kPion][index] ; - sdp[AliESDtrack::kMuon][index] = fDFmuon ->Eval(dispersion) ; //landau distribution - } - else{ - sdp[AliESDtrack::kPhoton][index] = 1. ; - sdp[AliESDtrack::kElectron][index] = 1. ; - sdp[AliESDtrack::kPion][index] = 1. ; - sdp[AliESDtrack::kKaon][index] = 1. ; - sdp[AliESDtrack::kProton][index] = 1. ; - sdp[AliESDtrack::kNeutron][index] = 1. ; - sdp[AliESDtrack::kEleCon][index] = 1. ; - sdp[AliESDtrack::kKaon0][index] = 1. ; - sdp[AliESDtrack::kMuon][index] = 1. ; - } + sdp[AliPID::kPhoton][index] = 1. ; + sdp[AliPID::kElectron][index] = 1. ; + sdp[AliPID::kPion][index] = 1. ; + sdp[AliPID::kKaon][index] = 1. ; + sdp[AliPID::kProton][index] = 1. ; + sdp[AliPID::kNeutron][index] = 1. ; + sdp[AliPID::kEleCon][index] = 1. ; + sdp[AliPID::kKaon0][index] = 1. ; + sdp[AliPID::kMuon][index] = 1. ; + if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){ + sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ; + sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ; + sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ; + sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ; + sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ; + sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ; + sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index]; + sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ; + sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ; + //landau distribution + } - // CPV-EMC Distance - Float_t distance = GetDistance(emc, cpv, "R") ; - // Info("MakePID", "Distance %f", distance); - Float_t pcpv = 0 ; - Float_t pcpvneutral = 0. ; - Float_t pcpvelectron = GausF (en , distance, fCPVelectron) ; - Float_t pcpvcharged = LandauF(en , distance, fCPVcharged) ; - //Float_t pcpvcharged = ChargedHadronDistProb(en , distance, fCPVchargedg, fCPVchargedl) ; - // Info("MakePID", "CPV: electron %f, hadron %f", pcpvelectron, pcpvcharged); +// Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion); +// Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]); +// cout<<">>>>>multiplicity "<GetMultiplicity()<<", dispersion "<< dispersion<GetCpvDistance("X")) ; + Float_t z = ts->GetCpvDistance("Z") ; + + Double_t pcpv = 0 ; + Double_t pcpvneutral = 0. ; + + Double_t elprobx = GausF(en , x, fXelectron) ; + Double_t elprobz = GausF(en , z, fZelectron) ; + Double_t chprobx = GausF(en , x, fXcharged) ; + Double_t chprobz = GausF(en , z, fZcharged) ; + Double_t pcpvelectron = elprobx * elprobz; + Double_t pcpvcharged = chprobx * chprobz; + +// cout<<">>>>energy "<>>>electron : x "<>>>hadron : x "<>>>electron : px*pz "<= pcpvcharged) pcpv = pcpvelectron ; else pcpv = pcpvcharged ; - if(pcpv < 1e-4) + if(pcpv < fChargedNeutralThreshold) { pcpvneutral = 1. ; pcpvcharged = 0. ; pcpvelectron = 0. ; } + // else + // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"< 30.){ // pi0 are detected via decay photon - stof[AliESDtrack::kPi0][index] = fTFphoton ->Eval(time) ; - scpv[AliESDtrack::kPi0][index] = pcpvneutral ; - if(emc->GetMultiplicity() > 4) - sdp [AliESDtrack::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; - else - sdp [AliESDtrack::kPi0][index] = 1. ; - } - else{ - stof[AliESDtrack::kPi0][index] = 0. ; - scpv[AliESDtrack::kPi0][index] = 0. ; - sdp [AliESDtrack::kPi0][index] = 0. ; - fInitPID[AliESDtrack::kPi0] = 0. ; + stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index]; + scpv[AliPID::kPi0][index] = pcpvneutral ; + if(emc->GetMultiplicity() > fDispMultThreshold) + sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ; + //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; +// cout<<"E = "<Eval(en) ; + + sw[AliPID::kPhoton][index] = 1. ; + sw[AliPID::kElectron][index] = 1. ; + sw[AliPID::kPion][index] = weight ; + sw[AliPID::kKaon][index] = weight ; + sw[AliPID::kProton][index] = weight ; + sw[AliPID::kNeutron][index] = weight ; + sw[AliPID::kEleCon][index] = 1. ; + sw[AliPID::kKaon0][index] = weight ; + sw[AliPID::kMuon][index] = weight ; + sw[AliPID::kPi0][index] = 1. ; + +// if(en > 0.5){ +// cout<<"######################################################"<>>>>multiplicity "<GetMultiplicity()<>>>electron : xprob "<>>>hadron : xprob "<>>>electron : px*pz "<(fRecParticles->At(index)); + + //Conversion electron? + + if(recpar->IsEleCon()){ + fInitPID[AliPID::kEleCon] = 1. ; + fInitPID[AliPID::kPhoton] = 0. ; + fInitPID[AliPID::kElectron] = 0. ; + } + else{ + fInitPID[AliPID::kEleCon] = 0. ; + fInitPID[AliPID::kPhoton] = 1. ; + fInitPID[AliPID::kElectron] = 1. ; + } + // fInitPID[AliPID::kEleCon] = 0. ; + + // calculates the Bayesian weight + Int_t jndex ; Double_t wn = 0.0 ; for (jndex = 0 ; jndex < kSPECIES ; jndex++) - //wn += stof[jndex][index] * pid[jndex] ; - wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * fInitPID[jndex] ; - //cout<<"*************wn "<RecParticle(index) ; + wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * + sw[jndex][index] * fInitPID[jndex] ; + + // cout<<"*************wn "<0) for (jndex = 0 ; jndex < kSPECIES ; jndex++) { //cout<<"jndex "<SetPID(jndex, stof[jndex][index] * sdp[jndex][index] * - scpv[jndex][index] * fInitPID[jndex] / wn) ; -// cout<<"final prob "<SetPID(jndex, stof[jndex][index] * fInitPID[jndex] / wn) ; - //cout<<"After SetPID"<Print(); + sw[jndex][index] * scpv[jndex][index] * + fInitPID[jndex] / wn) ; } } // Info("MakePID", "Delete"); - // for (Int_t i =0; i< kSPECIES; i++){ - // delete [] stof[i]; - // cout<EmcRecPoints() ; - TObjArray * cpvRecPoints = gime->CpvRecPoints() ; - TClonesArray * trackSegments = gime->TrackSegments() ; - if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { - Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; + if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) { + AliFatal("RecPoints or TrackSegments not found !") ; } - TClonesArray * recParticles = gime->RecParticles() ; - recParticles->Clear(); + fRecParticles->Clear(); - TIter next(trackSegments) ; + TIter next(fTrackSegments) ; AliPHOSTrackSegment * ts ; Int_t index = 0 ; AliPHOSRecParticle * rp ; while ( (ts = (AliPHOSTrackSegment *)next()) ) { - - new( (*recParticles)[index] ) AliPHOSRecParticle() ; - rp = (AliPHOSRecParticle *)recParticles->At(index) ; + // cout<<">>>>>>>>>>>>>>>PCA Index "<At(index) ; rp->SetTrackSegment(index) ; rp->SetIndexInList(index) ; AliPHOSEmcRecPoint * emc = 0 ; if(ts->GetEmcIndex()>=0) - emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; + emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; AliPHOSCpvRecPoint * cpv = 0 ; if(ts->GetCpvIndex()>=0) - cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; + cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ; Int_t track = 0 ; track = ts->GetTrackIndex() ; @@ -1096,32 +1363,32 @@ void AliPHOSPIDv1::MakeRecParticles() // Choose the cluster energy range if (!emc) { - Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ; + AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; } Float_t e = emc->GetEnergy() ; Float_t lambda[2] ; emc->GetElipsAxis(lambda) ; - + if((lambda[0]>0.01) && (lambda[1]>0.01)){ // Looking PCA. Define and calculate the data (X), // introduce in the function X2P that gives the components (P). - Float_t Spher = 0. ; - Float_t Emaxdtotal = 0. ; + Float_t spher = 0. ; + Float_t emaxdtotal = 0. ; if((lambda[0]+lambda[1])!=0) - Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); + spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); - Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); + emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); fX[0] = lambda[0] ; fX[1] = lambda[1] ; fX[2] = emc->GetDispersion() ; - fX[3] = Spher ; + fX[3] = spher ; fX[4] = emc->GetMultiplicity() ; - fX[5] = Emaxdtotal ; + fX[5] = emaxdtotal ; fX[6] = emc->GetCoreEnergy() ; fPrincipalPhoton->X2P(fX,fPPhoton); @@ -1145,7 +1412,7 @@ void AliPHOSPIDv1::MakeRecParticles() // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) // is set to 1 - if(GetCPVBit(emc, cpv, effPur,e) == 1 ){ + if(GetCPVBit(ts, effPur,e) == 1 ){ rp->SetPIDBit(effPur) ; //cout<<"CPV bit "<SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; rp->SetCalcMass(0); rp->Name(); //If photon sets the particle pdg name to gamma - rp->SetProductionVertex(0,0,0,0); + rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0); rp->SetFirstMother(-1); rp->SetLastMother(-1); rp->SetFirstDaughter(-1); rp->SetLastDaughter(-1); rp->SetPolarisation(0,0,0); //Set the position in global coordinate system from the RecPoint - AliPHOSGeometry * geom = gime->PHOSGeometry() ; - AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; - AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; + AliPHOSTrackSegment * ts = static_cast(fTrackSegments->At(rp->GetPHOSTSIndex())); + AliPHOSEmcRecPoint * erp = static_cast(fEMCRecPoints->At(ts->GetEmcIndex())); TVector3 pos ; - geom->GetGlobal(erp, pos) ; + fGeom->GetGlobalPHOS(erp, pos) ; rp->SetPos(pos); index++ ; } } //____________________________________________________________________________ -void AliPHOSPIDv1::Print() const +void AliPHOSPIDv1::Print(const Option_t *) const { // Print the parameters used for the particle type identification - Info("Print", "=============== AliPHOSPIDv1 ================") ; + AliInfo("=============== AliPHOSPIDv1 ================") ; printf("Making PID\n") ; printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; @@ -1212,7 +1478,7 @@ void AliPHOSPIDv1::Print() const printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; - Printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; + printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; fParameters->Print() ; } @@ -1223,23 +1489,17 @@ void AliPHOSPIDv1::PrintRecParticles(Option_t * option) { // Print table of reconstructed particles - AliPHOSGetter *gime = AliPHOSGetter::Instance() ; - - TClonesArray * recParticles = gime->RecParticles() ; - TString message ; - message = "\nevent " ; - message += gAlice->GetEvNumber() ; - message += " found " ; - message += recParticles->GetEntriesFast(); + message = " found " ; + message += fRecParticles->GetEntriesFast(); message += " RecParticles\n" ; if(strstr(option,"all")) { // printing found TS message += "\n PARTICLE Index \n" ; Int_t index ; - for (index = 0 ; index < recParticles->GetEntries() ; index++) { - AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; + for (index = 0 ; index < fRecParticles->GetEntries() ; index++) { + AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ; message += "\n" ; message += rp->Name().Data() ; message += " " ; @@ -1248,7 +1508,7 @@ void AliPHOSPIDv1::PrintRecParticles(Option_t * option) message += rp->GetType() ; } } - Info("Print", message.Data() ) ; + AliInfo(message.Data() ) ; } //____________________________________________________________________________ @@ -1287,26 +1547,26 @@ void AliPHOSPIDv1::SetParameters() // lines 14-15: parameters to calculate border for high-pt photons and pi0 fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); - fParameters = new TMatrix(16,4) ; - const Int_t maxLeng=255; - char string[maxLeng]; + fParameters = new TMatrixF(16,4) ; + const Int_t kMaxLeng=255; + char string[kMaxLeng]; // Open a text file with PID parameters FILE *fd = fopen(fFileNameParameters.Data(),"r"); if (!fd) - Fatal("SetParameter","File %s with a PID parameters cannot be opened\n", - fFileNameParameters.Data()); + AliFatal(Form("File %s with a PID parameters cannot be opened\n", + fFileNameParameters.Data())); Int_t i=0; // Read parameter file line-by-line and skip empty line and comments - while (fgets(string,maxLeng,fd) != NULL) { + while (fgets(string,kMaxLeng,fd) != NULL) { if (string[0] == '\n' ) continue; if (string[0] == '!' ) continue; sscanf(string, "%f %f %f %f", &(*fParameters)(i,0), &(*fParameters)(i,1), &(*fParameters)(i,2), &(*fParameters)(i,3)); i++; - //Info("SetParameters", "line %d: %s",i,string); + AliDebug(1, Form("SetParameters", "line %d: %s",i,string)); } fclose(fd); } @@ -1315,9 +1575,9 @@ void AliPHOSPIDv1::SetParameters() void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) { // Set parameter "Calibration" i to a value param - if(i>2 || i<0) - Error("SetParameterCalibration","Invalid parameter number: %d",i); - else + if(i>2 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else (*fParameters)(0,i) = param ; } @@ -1327,13 +1587,15 @@ void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on // Purity-Efficiency point i - if(i>2 || i<0) - Error("SetParameterCpv2Emc","Invalid parameter number: %d",i); - else { + if(i>2 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else { axis.ToLower(); if (axis == "x") (*fParameters)(1,i) = cut; else if (axis == "z") (*fParameters)(2,i) = cut; - else Error("SetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); + else { + AliError(Form("Invalid axis name: %s",axis.Data())); + } } } @@ -1341,9 +1603,9 @@ void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) { // Set parameter "Hard photon boundary" i to a value param - if(i>4 || i<0) - Error("SetParameterPhotonBoundary","Invalid parameter number: %d",i); - else + if(i>4 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else (*fParameters)(14,i) = param ; } @@ -1351,9 +1613,9 @@ void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) { // Set parameter "Hard pi0 boundary" i to a value param - if(i>1 || i<0) - Error("SetParameterPi0Boundary","Invalid parameter number: %d",i); - else + if(i>1 || i<0) { + AliError(Form("Invalid parameter number: %d",i)); + } else (*fParameters)(15,i) = param ; } @@ -1361,9 +1623,9 @@ void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) { // Set the parameter TimeGate depending on Purity-Efficiency point i - if (i>2 || i<0) - Error("SetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); - else + if (i>2 || i<0) { + AliError(Form("Invalid Efficiency-Purity choice %d",i)); + } else (*fParameters)(3,i)= gate ; } @@ -1381,61 +1643,48 @@ void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString par if (particle == "photon") offset=0; else if (particle == "pi0") offset=5; else - Error("SetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); + AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", + particle.Data())); if (param.Contains("a")) p=4+offset; else if(param.Contains("b")) p=5+offset; else if(param.Contains("c")) p=6+offset; else if(param.Contains("x0"))p=7+offset; else if(param.Contains("y0"))p=8+offset; - if((i>4)||(i<0)) - Error("SetEllipseParameter", "No parameter with index %d", i) ; - else if(p==-1) - Error("SetEllipseParameter", "No parameter with name %s", param.Data() ) ; - else + if((i>4)||(i<0)) { + AliError(Form("No parameter with index %d", i)) ; + } else if(p==-1) { + AliError(Form("No parameter with name %s", param.Data() )) ; + } else (*fParameters)(p,i) = par ; } //____________________________________________________________________________ -void AliPHOSPIDv1::Unload() -{ - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - gime->PhosLoader()->UnloadRecPoints() ; - gime->PhosLoader()->UnloadTracks() ; - gime->PhosLoader()->UnloadRecParticles() ; -} - -//____________________________________________________________________________ -void AliPHOSPIDv1::WriteRecParticles() -{ +void AliPHOSPIDv1::GetVertex(void) +{ //extract vertex either using ESD or generator - AliPHOSGetter *gime = AliPHOSGetter::Instance() ; - - TClonesArray * recParticles = gime->RecParticles() ; - recParticles->Expand(recParticles->GetEntriesFast() ) ; - if(fWrite){ - TTree * treeP = gime->TreeP(); - - //First rp - Int_t bufferSize = 32000 ; - TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); - rpBranch->SetTitle(BranchName()); - - rpBranch->Fill() ; - - gime->WriteRecParticles("OVERWRITE"); - gime->WritePID("OVERWRITE"); + //Try to extract vertex from data + if(fESD){ + const AliESDVertex *esdVtx = fESD->GetVertex() ; + if(esdVtx && esdVtx->GetChi2()!=0.){ + fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ; + return ; + } } -} - + // Use vertex diamond from CDB GRP folder if the one from ESD is missing + // PLEASE FIX IT + AliWarning("Can not read vertex from data, use fixed \n") ; + fVtx.SetXYZ(0.,0.,0.) ; + +} //_______________________________________________________________________ void AliPHOSPIDv1::SetInitPID(const Double_t *p) { // Sets values for the initial population of each particle type - for (Int_t i=0; i