From 50739f1510907d8669c9ee981a9083f4a746dd2d Mon Sep 17 00:00:00 2001 From: schutz Date: Mon, 21 Oct 2002 16:41:31 +0000 Subject: [PATCH] Particle identification improvment by Gustavo --- PHOS/AliPHOSPIDv1.cxx | 640 +++++++++++++++++------------------------- PHOS/AliPHOSPIDv1.h | 103 +++---- PHOS/Parameters.dat | 39 ++- 3 files changed, 314 insertions(+), 468 deletions(-) diff --git a/PHOS/AliPHOSPIDv1.cxx b/PHOS/AliPHOSPIDv1.cxx index b752e0031fd..2e4704bfed7 100644 --- a/PHOS/AliPHOSPIDv1.cxx +++ b/PHOS/AliPHOSPIDv1.cxx @@ -28,42 +28,56 @@ // -Bit 3 to 5: bit set if TOF < fTimeGate (each bit corresponds // to a different efficiency-purity point of the photon identification) // -Bit 6 to 9: bit set if Principal Components are -// inside an ellipse defined by fX_center, fY_center, fA, fB, fAngle +// inside an ellipse defined by the parameters a, b, c, x0 and y0. // (each bit corresponds to a different efficiency-purity point of the -// photon identification) +// photon identification) +// The PCA (Principal components analysis) needs a file that contains +// a previous analysis of the correlations between the particles. This +// file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis don for +// energies between 0.5 and 100 GeV. // A calibrated energy is calculated. The energy of the reconstructed -// cluster is corrected with the formula A + B * E + C * E^2, whose parameters -// where obtained thourgh the study of the reconstructed energy -// distribution of monoenergetic photons. -// +// cluster is corrected with the formula A + B * E + C * E^2, whose +// parameters where obtained thourgh the study of the reconstructed +// energy distribution of monoenergetic photons. // +// All the parameters (RCPV(6 rows-3 columns),TOF(6r-3c),PCA(5r-4c) +// and calibration(1r-3c))are stored in a file called +// $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is +// initialized, this parameters are copied to a Matrix (18,4), a +// TMatrixD object. // // use case: -// root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root","v1") +// root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") // Warning in : object already instantiated -// // reading headers from file galice1.root and create RecParticles with title v1 - // TrackSegments and RecPoints with title "v1" are used -// // set file name for the branch RecParticles +// // reading headers from file galice1.root and create RecParticles + // TrackSegments and RecPoints are used +// // set file name for the branch RecParticles // root [1] p->ExecuteTask("deb all time") -// // available options -// // "deb" - prints # of reconstructed particles -// // "deb all" - prints # and list of RecParticles -// // "time" - prints benchmarking results +// // available options +// // "deb" - prints # of reconstructed particles +// // "deb all" - prints # and list of RecParticles +// // "time" - prints benchmarking results // -// root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1","v0") +// root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) // Warning in : object already instantiated -// // reading headers from file galice1.root and create RecParticles with title v1 - // RecPoints and TrackSegments with title "v0" are used +// //Split mode. // root [3] p2->ExecuteTask() // -// There are two possible principal files available to do the analysis. -// One for energy ranges from 0.5 to 5 GeV, and another -// one from 5 to 100 GeV. This files are automatically called in function -// of the cluster energy. + //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & // Gustavo Conesa April 2002 - +// PCA redesigned by Gustavo Conesa October 2002: +// The way of using the PCA has changed. Instead of 2 +// files with the PCA, each one with different energy ranges +// of application, we use the wide one (0.5-100 GeV), and instead +// of fixing 3 elipses for different ranges of energy, it has been +// studied the dependency of the ellipses parameters with the +// energy, and they are implemented in the code as a funtion +// of the energy. +// +// +// // --- ROOT system --- #include "TROOT.h" #include "TTree.h" @@ -131,8 +145,7 @@ AliPHOSPIDv1::~AliPHOSPIDv1() delete [] fX ; // Principal input delete [] fP ; // Principal components // delete fParameters ; // Matrix of Parameters -// delete fParameters5 ; // Matrix of Parameters -// delete fParameters100 ; // Matrix of Parameters + if (!fDefaultInit) { @@ -240,7 +253,7 @@ void AliPHOSPIDv1::InitParameters() } //____________________________________________________________________________ -Double_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur) +const Double_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur) const { // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and // Purity-Efficiency point (possible options "HIGH EFFICIENCY" @@ -248,16 +261,15 @@ Double_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const T // EFFICIENCY by PURITY) Int_t eff_pur = GetEffPurOption(Eff_Pur); - - GetAnalysisParameters(Cluster_En) ; - if((fClusterrcpv!= -1)&&(eff_pur != -1)) - return (*fParameters)(fClusterrcpv,eff_pur) ; + Int_t cluster = GetClusterOption(Cluster_En) ; + if((cluster!= -1)&&(eff_pur != -1)) + return (*fParameters)(cluster,eff_pur) ; else return 0.0; } //____________________________________________________________________________ -Double_t AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) +const Double_t AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) const { // Get TimeGate parameter depending on the cluster energy and // Purity-Efficiency point (possible options "HIGH EFFICIENCY" @@ -265,16 +277,15 @@ Double_t AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_ // EFFICIENCY by PURITY) Int_t eff_pur = GetEffPurOption(Eff_Pur); - GetAnalysisParameters(Cluster_En) ; - - if((fCluster!= -1)&&(eff_pur != -1)) - return (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur) ; + Int_t cluster = GetClusterOption(Cluster_En) ; + if((cluster!= -1)&&(eff_pur != -1)) + return (*fParameters)(cluster+6,eff_pur) ; else return 0.0; } //_____________________________________________________________________________ -Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const +const Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const { // Calculates the distance between the EMC RecPoint and the PPSD RecPoint @@ -303,159 +314,46 @@ Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cp } //____________________________________________________________________________ -Double_t AliPHOSPIDv1::CalibratedEnergy(Float_t e){ - //It calibrates Energy depending on the recpoint energy. -// The energy of the reconstructed -// cluster is corrected with the formula A + B* E + C* E^2, whose parameters -// where obtained through the study of the reconstructed energy -// distribution of monoenergetic photons. - Double_t enerec; - enerec = fACalParameter + fBCalParameter * e+ fCCalParameter * e * e; +const Double_t AliPHOSPIDv1::GetCalibratedEnergy(const Float_t e) const +{ +// It calibrates Energy depending on the recpoint energy. +// The energy of the reconstructed cluster is corrected with +// the formula A + B* E + C* E^2, whose parameters where obtained +// through the study of the reconstructed energy distribution of +// monoenergetic photons. + + Double_t p[]={0.,0.,0.}; + Int_t i; + for(i=0;i<3;i++) p[i]= (*fParameters)(17,i); + Double_t enerec = p[0] + p[1]* e+ p[2] * e * e; return enerec ; } //____________________________________________________________________________ -Int_t AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur)const +const Int_t AliPHOSPIDv1::GetPrincipalSign(const Double_t* P,const Int_t eff_pur, const Float_t E)const { - //This method gives if the PCA of the particle are inside a defined ellipse - // Get the parameters that define the ellipse stored in the - // fParameters matrix. - Double_t X_center = (*fParameters)(cluster+6,eff_pur) ; - Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ; - Double_t A = (*fParameters)(cluster+12,eff_pur) ; - Double_t B = (*fParameters)(cluster+15,eff_pur) ; - Double_t Angle = (*fParameters)(cluster+18,eff_pur) ; - - Int_t prinsign; - Double_t Dx = 0. ; - Double_t Delta = 0. ; - Double_t Y = 0. ; - Double_t Y_1 = 0. ; - Double_t Y_2 = 0. ; - Double_t Pi = TMath::Pi() ; - Double_t Cos_Theta = TMath::Cos(Pi*Angle/180.) ; - Double_t Sin_Theta = TMath::Sin(Pi*Angle/180.) ; - - Dx = P[0] - X_center ; - Delta = 4.*A*A*B*B* (A*A*Cos_Theta*Cos_Theta - + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ; - if (Delta < 0.) - {prinsign=0;} - - else if (Delta == 0.) - { - Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx / - (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ; - Y += Y_center ; - if(P[1]==Y ) - {prinsign=1;} - else - {prinsign=0;} - } - else - { - Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx + - TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta + - B*B*Sin_Theta*Sin_Theta) ; - Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx - - TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta - + B*B*Sin_Theta*Sin_Theta) ; - Y_1 += Y_center ; - Y_2 += Y_center ; - if ((P[1]<=Y_1) && (P[1]>=Y_2)) - {prinsign=1;} - else - {prinsign=0;} - } + //Is the particle inside de PCA ellipse? + + Int_t prinsign= 0 ; + Double_t A = GetEllipseParameter("a", E); + Double_t B = GetEllipseParameter("b", E); + Double_t C = GetEllipseParameter("c", E); + Double_t X_center = GetEllipseParameter("x0", E); + Double_t Y_center = GetEllipseParameter("y0", E); + + Double_t R = TMath::Power((P[0] - X_center)/A,2) + + TMath::Power((P[1] - Y_center)/B,2) + + C*(P[0] - X_center)*(P[1] - Y_center)/(A*B) ; + //3 different ellipses defined + if((eff_pur==2)&&(R <1./2.)) prinsign= 1; + if((eff_pur==1)&&(R <2. )) prinsign= 1; + if((eff_pur==0)&&(R <9./2.)) prinsign= 1; + + if(R<0)cout<<"Error: Negative square?"<ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_5.dat"); - fFileName100 = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; - fFileNamePar100 = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_100.dat"); - - //SetPrincipalFileOptions(); - //fOptFileName); - TFile f5( fFileName5.Data(), "read" ) ; - fPrincipal5 = dynamic_cast (f5.Get("principal")) ; - f5.Close() ; - TFile f100( fFileName100.Data(), "read" ) ; - fPrincipal100 = dynamic_cast (f100.Get("principal")) ; - f100.Close() ; - TFile f( fFileName100.Data(), "read" ) ; + // Open principal and parameters files to be used + + fFileName = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; + fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); + TFile f( fFileName.Data(), "read" ) ; fPrincipal = dynamic_cast (f.Get("principal")) ; f.Close() ; - // Initialization of the Parameters matrix. In the File ParametersXX.dat - // are all the parameters. These are introduced in a matrix of 21x3 or 22x3 - // elements (depending on the principal file 21 rows for 0.5-5 GeV and 22 - // rows for 5-100). - // All the parameters defined in this file are, in order of row (there are - // 3 rows per parameter): CpvtoEmcDistanceCut(if the principal file is 5-100 - // GeV then 4 rows), TimeGate and the ellipse parameters, X_center, Y_center, - // a, b, angle. Each row of a given parameter depends on the cluster energy range - // (wich depends on the chosen principal file) - // Each column designs the parameters for a point in the Efficiency-Purity - // of the photon identification P1(96%,63%), P2(87%,0.88%) and P3(68%,94%) - // for the principal file from 0.5-5 GeV and for the other one P1(95%,79%), - // P2(89%,90%) and P3(72%,96%) - - fEnergyAnalysisCut = 5.; // Energy cut to change PCA - - fParameters5 = new TMatrixD(21,3) ; - fParameters100 = new TMatrixD(22,3) ; - fParameters = new TMatrixD(22,3) ; - - ifstream paramFile5(fFileNamePar5) ; - - Int_t i,j ; - for(i = 0; i< 21; i++){ - for(j = 0; j< 3; j++){ - paramFile5 >> (*fParameters5)(i,j) ; - } - } - paramFile5.close(); - - ifstream paramFile100(fFileNamePar100) ; + // Initialization of the Parameters matrix. In the File Parameters.dat + // are all the parameters. These are introduced in a matrix of 18x4 + // + // All the parameters defined in this file are, in order of row: + // CpvtoEmcDistanceCut (6 rows, each one depends on the energy range of the + // particle, and 3 columns, each one depending on the efficiency-purity + // point), TimeGate (the same) and the parameters of the functions that + // calculate the ellipse parameters, x0,y0,a, b, c. These 5 parameters + // (5 rows) depend on 4 parameters (columns). Finally there is a row with + // the energy calibration parameters, 3 parameters. - Int_t l,k ; + fParameters = new TMatrixD(18,4) ; - for(l = 0; l< 22; l++){ - for(k = 0; k< 3; k++){ - paramFile100 >> (*fParameters100)(l,k) ; - } - } - paramFile100.close(); - - ifstream paramFile(fFileNamePar100) ; + ifstream paramFile(fFileNamePar) ; Int_t h,n; - for(h = 0; h< 22; h++){ - for(n = 0; n< 3; n++){ + for(h = 0; h< 18; h++){ + for(n = 0; n< 4; n++){ paramFile >> (*fParameters)(h,n) ; } } paramFile.close(); - - fCluster = -1; - fClusterrcpv = -1; - fMatrixExtraRow = 0; - - //Calibration parameters Encal = C * E^2 + B * E + A (E is the energy from cluster) - fACalParameter = 0.0241 ; - fBCalParameter = 1.0504 ; - fCCalParameter = 0.000249 ; - - // fParameters->Print(); } //_____________________________________________________________________________ -void AliPHOSPIDv1::GetAnalysisParameters(Float_t Cluster_En) +const Int_t AliPHOSPIDv1::GetClusterOption(const Float_t Cluster_En) const { - if(Cluster_En <= fEnergyAnalysisCut){ - fPrincipal = fPrincipal5; - fParameters = fParameters5; - fMatrixExtraRow = 0; - GetClusterOption(Cluster_En,kFALSE) ; - } - else{ - fPrincipal = fPrincipal100; - fParameters = fParameters100; - fMatrixExtraRow = 1; - GetClusterOption(Cluster_En,kTRUE) ; - } -} - -//_____________________________________________________________________________ -void AliPHOSPIDv1::GetClusterOption(const Float_t Cluster_En, const Bool_t range) -{ - - // Gives the cluster energy range. - // range = kFALSE Default analysis range from 0.5 to 5 GeV - // range = kTRUE analysis range from 0.5 to 100 GeV - - - //Int_t cluster = -1 ; - - if((range == kFALSE)){ - if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)){ - fCluster = 0 ; - fClusterrcpv = 0 ; - } - if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)){ - fCluster = 1 ; - fClusterrcpv = 1 ; - } - if( Cluster_En > 2.0){ - fCluster = 2 ; - fClusterrcpv = 2 ; - } - } - else if(range == kTRUE){ - if((Cluster_En > 0.5 )&&(Cluster_En <= 20.0)) fCluster = 0 ; - if((Cluster_En > 20.0)&&(Cluster_En <= 50.0)) fCluster = 1 ; - if( Cluster_En > 50.0) fCluster = 2 ; - if((Cluster_En > 5.0 )&&(Cluster_En <= 10.0)) fClusterrcpv = 0 ; - if((Cluster_En > 10.0)&&(Cluster_En <= 20.0)) fClusterrcpv = 1 ; - if((Cluster_En > 20.0)&&(Cluster_En <= 30.0)) fClusterrcpv = 2 ; - if( Cluster_En > 30.0) fClusterrcpv = 3 ; - } - else { - fCluster = -1 ; - fClusterrcpv = -1; - cout<<"Invalid Energy option"< 0.0 )&&(Cluster_En <= 2.0 )) cluster = 0 ; + if((Cluster_En > 2.0 )&&(Cluster_En <= 5.0 )) cluster = 1 ; + if((Cluster_En > 5.0 )&&(Cluster_En <= 10.0)) cluster = 2 ; + if((Cluster_En > 10.0)&&(Cluster_En <= 20.0)) cluster = 3 ; + if((Cluster_En > 20.0)&&(Cluster_En <= 30.0)) cluster = 4 ; + if( Cluster_En > 30.0) cluster = 5 ; + + return cluster; } //____________________________________________________________________________ -Int_t AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const +const Int_t AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const { // Looks for the Purity-Efficiency point (possible options "HIGH EFFICIENCY" @@ -659,6 +462,92 @@ Int_t AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const return eff_pur; } +//________________________________________________________________________ +void AliPHOSPIDv1::SetEllipseParameter(TString Param, Int_t i, Double_t par) +{ + // Set the parameter "i" that is needed to calculate the ellipse + // parameter "Param". + + Int_t p= -1; + + if(Param.Contains("a"))p=12; + if(Param.Contains("b"))p=13; + if(Param.Contains("c"))p=14; + if(Param.Contains("x0"))p=15; + if(Param.Contains("y0"))p=16; + if((i>4)||(i<0)) + cout<<"Error:: No parameter with index "<4)||(i<0)) + cout<<"Error:: No parameter with index "<70.)E=70.; + } + + if(Param.Contains("b")){ + for(i=0;i<4;i++)p[i]=(*fParameters)(13,i); + if(E>70.)E=70.; + } + + if(Param.Contains("c")) + for(i=0;i<4;i++)p[i]=(*fParameters)(14,i); + + if(Param.Contains("x0")){ + for(i=0;i<4;i++)p[i]=(*fParameters)(15,i); + if(E<1.)E=1.1; + } + if(Param.Contains("y0")) + for(i=0;i<4;i++)p[i]=(*fParameters)(16,i); + + value = p[0]/TMath::Sqrt(E)+p[1]*E+p[2]*E*E+p[3]; + return value; +} //____________________________________________________________________________ void AliPHOSPIDv1::Exec(Option_t * option) @@ -787,61 +676,64 @@ void AliPHOSPIDv1::MakeRecParticles(){ abort(); } Float_t e = emc->GetEnergy() ; - - GetAnalysisParameters(e);// Gives value to fCluster, fClusterrcpv, fMatrixExtraRow, and to fPrincipal and fParameters depending on the energy. - - if((fCluster== -1)||(fClusterrcpv == -1)) continue ; - + Int_t cluster = GetClusterOption(e) ;// Gives value to cluster that defines the energy range parameter to be used in de RCPV, TOF and used in the PCA. + if(cluster== -1) continue ; + 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. ; + + if((lambda[0]+lambda[1])!=0) Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); + + Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); + + fX[0] = lambda[0] ; + fX[1] = lambda[1] ; + fX[2] = emc->GetDispersion() ; + fX[3] = Spher ; + fX[4] = emc->GetMultiplicity() ; + fX[5] = Emaxdtotal ; + fX[6] = emc->GetCoreEnergy() ; + + fPrincipal->X2P(fX,fP); + } + else{ + fP[0]=-100.0; //We do not accept clusters with + fP[1]=-100.0; //one cell as a photon-like + } + Float_t time =emc->GetTime() ; - if((lambda[0]>0.01) && (lambda[1]>0.01) && time > 0.){ + // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken + // into account to set the particle identification) + for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){ + + // 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(GetDistance(emc, cpv, "R") > (*fParameters)(cluster,eff_pur) ) + rp->SetPIDBit(eff_pur) ; - // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken - // into account to set the particle identification) - for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){ - - // 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(GetDistance(emc, cpv, "R") > (*fParameters)(fClusterrcpv,eff_pur) ) - rp->SetPIDBit(eff_pur) ; - - // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th - // bit (depending on the efficiency-purity point )is set to 1 - if(time< (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur)) - rp->SetPIDBit(eff_pur+3) ; - - // 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. ; - - if((lambda[0]+lambda[1])!=0) Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); - - Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); - - fX[0] = lambda[0] ; - fX[1] = lambda[1] ; - fX[2] = emc->GetDispersion() ; - fX[3] = Spher ; - fX[4] = emc->GetMultiplicity() ; - fX[5] = Emaxdtotal ; - fX[6] = emc->GetCoreEnergy() ; - - fPrincipal->X2P(fX,fP); - - //If we are inside the ellipse, 7th, 8th or 9th - // bit (depending on the efficiency-purity point )is set to 1 - if(GetPrincipalSign(fP,fCluster+fMatrixExtraRow,eff_pur) == 1) - rp->SetPIDBit(eff_pur+6) ; - + // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th + // bit (depending on the efficiency-purity point )is set to 1 + if(time< (*fParameters)(cluster+6,eff_pur)) { + rp->SetPIDBit(eff_pur+3) ; } + + //If we are inside the ellipse, 7th, 8th or 9th + // bit (depending on the efficiency-purity point )is set to 1 + if(GetPrincipalSign(fP,eff_pur,e) == 1) + rp->SetPIDBit(eff_pur+6) ; } - + //Set momentum, energy and other parameters - Float_t encal = CalibratedEnergy(e); + Float_t encal = GetCalibratedEnergy(e); TVector3 dir = GetMomentumDirection(emc,cpv) ; dir.SetMag(encal) ; rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; @@ -868,26 +760,14 @@ void AliPHOSPIDv1:: Print() // cout << " RecPoints branch title: " << fRecPointsTitle.Data() << endl ; // cout << " TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ; // cout << " RecParticles Branch title " << fRecParticlesTitle.Data() << endl; - - cout << " Pricipal analysis file from 0.5 to 5 " << fFileName5.Data() << endl; - cout << " Name of parameters file "<Print() ; - - cout << " Pricipal analysis file from 5 to 100 " << fFileName100.Data() << endl; - cout << " Name of parameters file "<Print() ; - - cout << " Energy Calibration Parameters A + B* E + C * E^2"<Print() ; cout << "============================================" << endl ; } diff --git a/PHOS/AliPHOSPIDv1.h b/PHOS/AliPHOSPIDv1.h index 2e26e29c6e8..abfa6679d34 100644 --- a/PHOS/AliPHOSPIDv1.h +++ b/PHOS/AliPHOSPIDv1.h @@ -27,92 +27,75 @@ class AliPHOSRecPoint ; #include "AliPHOSPID.h" class AliPHOSPIDv1 : public AliPHOSPID { - -public: - + + public: + AliPHOSPIDv1() ; // ctor AliPHOSPIDv1(const char* headerFile, const char * tsBranch = "Default", const Bool_t toSplit=kFALSE) ; - + virtual ~AliPHOSPIDv1() ; // dtor - + virtual void Exec(Option_t * option) ; // virtual char * GetRecParticlesBranch()const {return (char*) fRecParticlesTitle.Data() ;} // virtual char * GetTrackSegmentsBranch()const{return (char*) fTrackSegmentsTitle.Data(); } virtual const Int_t GetRecParticlesInRun() const {return fRecParticlesInRun ;} - + virtual void Print(Option_t * option) const {} void Print() ; - // Get CpvtoEmcDistanceCut and TimeGate parameters depending on the custer energy and - // Purity-Efficiency point (possible options "HIGH EFFICIENCY" "MEDIUM EFFICIENCY" "LOW - // EFFICIENCY" and 3 more options changing EFFICIENCY by PURITY) - Double_t GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur) ; - Double_t GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) ; - + //Get files that contain the PCA - const TString GetPrincipalFile5( )const {return fFileName5 ;} - const TString GetPrincipalFilePar5( )const {return fFileNamePar5 ;} - const TString GetPrincipalFile100( )const {return fFileName100 ;} - const TString GetPrincipalFilePar100( )const {return fFileNamePar100 ;} - - // Set all parameters necessary in the PID depending on the custer energy and - // Purity-Efficiency point (possible options "HIGH EFFICIENCY" "MEDIUM EFFICIENCY" "LOW - // EFFICIENCY" and 3 more options changing EFFICIENCY by PURITY) + const TString GetPrincipalFile( )const {return fFileName ;} + const TString GetPrincipalFilePar( )const {return fFileNamePar ;} + + // Set and Get all parameters necessary in the PID depending on the + // custer energy and Purity-Efficiency point (possible options "HIGH + // EFFICIENCY" "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options + // changing EFFICIENCY by PURITY) void SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut) ; void SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate) ; - void SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x) ; - void SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y) ; - void SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a) ; - void SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b) ; - void SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle) ; - void SetEllipseParameters(Float_t Cluster_En, TString Eff_Pur, Float_t x, Float_t y,Float_t a, Float_t b,Float_t angle) ; + const Double_t GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur) const ; + const Double_t GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) const; + + void SetEllipseParameter(TString Param, Int_t i, Double_t par) ; + const Double_t GetParameterToCalculateEllipse(const TString Param, const Int_t i) const ; + const Double_t GetEllipseParameter(const TString Param, Float_t E) const; //Get and Set energy calibration parameters - Float_t GetACalParameter() {return fACalParameter ;} - Float_t GetBCalParameter() {return fBCalParameter ;} - Float_t GetCCalParameter() {return fCCalParameter ;} - void SetACalParameter(Float_t a) { fACalParameter = a ;} - void SetBCalParameter(Float_t b) { fBCalParameter = b ;} - void SetCCalParameter(Float_t c) { fCCalParameter = c ;} - - - Float_t GetEnergyAnalysisCut() {return fEnergyAnalysisCut ;} - void SetEnergyAnalysisCut(Float_t e) { fEnergyAnalysisCut = e ;} - + + void SetCalibrationParameter(Int_t Param,Double_t param); + const Double_t GetCalibrationParameter(const Int_t i) const; + const Double_t GetCalibratedEnergy(const Float_t e) const; //Calibrates energy. + // virtual void SetTrackSegmentsBranch(const char* title) { fTrackSegmentsTitle = title;} // virtual void SetRecParticlesBranch (const char* title) { fRecParticlesTitle = title;} virtual const char * Version() const { return "pid-v1" ; } private: - + const TString BranchName() const ; virtual void Init() ; virtual void InitParameters() ; void MakeRecParticles(void ) ; // Relative Distance CPV-EMC - Float_t GetDistance(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv, Option_t * Axis)const ; - Int_t GetPrincipalSign(Double_t* P, Int_t ell, Int_t eff_pur)const ; //Principal cut + const Float_t GetDistance(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv, Option_t * Axis)const ; + const Int_t GetPrincipalSign(const Double_t* P, const Int_t eff_pur, const Float_t E)const ; //Principal cut TVector3 GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const ; void PrintRecParticles(Option_t * option) ; // Gives in wich cluster energy range is the event - void GetClusterOption(const Float_t Cluster_En,const Bool_t range) ; + const Int_t GetClusterOption(const Float_t Cluster_En) const; // Gives the Efficiency-Purity point. - Int_t GetEffPurOption(const TString Eff_Pur)const ; + const Int_t GetEffPurOption(const TString Eff_Pur)const ; virtual void WriteRecParticles(Int_t event) ; void SetParameters() ; //Fills the matrix of parameters - //Selects principal and parameters file in function of energy range. - void GetAnalysisParameters(Float_t ClusterEn) ; - Double_t CalibratedEnergy(Float_t e) ; //Calibrates energy. + private: Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) - TString fFileName5 ; // File that contains the Principal file for analysis from 0.5 to 5 GeV - TString fFileName100 ; // File that contains the Principal file for analysis from 0.5 to 100 GeV - TString fFileNamePar5 ; // File that contains the parameters for analysis from 0.5 to 5 GeV - TString fFileNamePar100 ;// File that contains the parameters for analysis from 0.5 to 100 GeV + TString fFileName ; // File that contains the Principal file for analysis + TString fFileNamePar ;// File that contains the parameters for analysis - // TString fFrom ; // name of Recpoints and TrackSegments // TString fHeaderFileName ; // file name with event header // TString fTrackSegmentsTitle; // branch name with track segments @@ -123,27 +106,13 @@ public: // AliPHOSClusterizer * fClusterizer ; //! clusterizer // AliPHOSTrackSegmentMaker * fTSMaker ; //! track segment maker - TPrincipal * fPrincipal5 ; //! TPrincipal from fFileName5 - TPrincipal * fPrincipal100 ; //! TPrincipal from fFileName100 - TPrincipal * fPrincipal ; //! TPrincipal copy + TPrincipal * fPrincipal ; //! TPrincipal from pca file Double_t * fX ; //! Principal data Double_t * fP ; //! Principal eigenvalues - Int_t fRecParticlesInRun ; //! Total number of recparticles in one run - - TMatrixD * fParameters5 ; //! Matrix of identification Parameters 0.5 to 5 GeV - TMatrixD * fParameters100 ; //! Matrix of identification Parameters 5-100 GeV - TMatrixD * fParameters; //! Matrix copy of identification Parameters - Float_t fEnergyAnalysisCut; // Energy to change from one PCA to the other. - Int_t fCluster; // Cluster energy range to choose parameters - Int_t fClusterrcpv; // Cluster energy range to choos rcpv parameters - Int_t fMatrixExtraRow; // Different size of the parameters file. Depends on range - - Float_t fACalParameter ;// A parameter energy calibration Encal=A+B*En+C*En^2 - Float_t fBCalParameter ;// B parameter energy calibration Encal=A+B*En+C*En^2 - Float_t fCCalParameter ;// B parameter energy calibration Encal=A+B*En+C*En^2 + TMatrixD * fParameters; //! Matrix of identification Parameters - ClassDef( AliPHOSPIDv1,5) // Particle identifier implementation version 1 + ClassDef( AliPHOSPIDv1,6) // Particle identifier implementation version 1 }; diff --git a/PHOS/Parameters.dat b/PHOS/Parameters.dat index dca881dbe47..10d906c2493 100644 --- a/PHOS/Parameters.dat +++ b/PHOS/Parameters.dat @@ -1,21 +1,18 @@ -1.0 6.0 6.0 -1.0 6.0 6.0 -1.0 4.0 4.0 -1.7e-8 1.7e-8 1.7e-8 -1.7e-8 1.7e-8 1.7e-8 -1.7e-8 1.7e-8 1.7e-8 -1.6 1.6 2.15 -2.1 1.9 2.5 -2.35 2.35 2.85 -1.1 1.1 1.5 -0.8 0.5 0.5 --0.1 -0.2 -0.2 -2.35 1.65 1.65 -2.15 1.9 1.9 -2.9 2.8 2.8 -1.6 1.05 1.05 -1.3 1.15 1.15 -1.75 1.5 1.5 -50. 50. 105. -50. 50. 70. -80. 80. 80. +1.0 6.0 6.0 0.0 +1.0 4.0 4.0 0.0 +1.0 3.5 3.5 0.0 +1.0 3.0 3.0 0.0 +1.0 2.5 2.5 0.0 +1.0 2.0 2.0 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +1.7e-8 1.7e-8 1.65e-8 0.0 +0.520 8.4e-3 -6.0e-5 0.235 +0.469 8.3e-3 -5.6e-5 0.522 +0.36 1.4e-2 -8.0e-5 -0.99 +-0.747 1.91e-2 4.0e-6 1.49 +0.91 -3.03e-2 1.08e-4 1.03 +0.0241 1.0504 0.000249 0.0 -- 2.31.1