/* History of cvs commits:
*
* $Log$
+ * 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
*
#include "TPrincipal.h"
#include "TFile.h"
#include "TSystem.h"
+#include "TVector3.h"
// --- AliRoot header files ---
//#include "AliLog.h"
-#include "AliGenerator.h"
#include "AliPHOS.h"
#include "AliPHOSPIDv1.h"
#include "AliPHOSGetter.h"
+#include "AliESD.h"
+#include "AliESDVertex.h"
ClassImp( AliPHOSPIDv1)
if(gime->TrackSegments() && //Skip events, where no track segments made
gime->TrackSegments()->GetEntriesFast()) {
+ GetVertex() ;
MakeRecParticles() ;
if(fBayesian)
}
+//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)
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) ;
+ 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)))
emc->GetGlobalPosition(dir, dummy) ;
- //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) ;
+ TVector3 local ;
+ emc->GetLocalPosition(local) ;
+
+ 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()/vInc.Y() ;
+ depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/vInc.Y() ;
}
else{
- xo=yo=zo=0.;
+ AliError("Cluster with zero energy \n");
+ }
+ //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()/vInc.Y() ;
+ depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/vInc.Y() ;
}
- TVector3 origin(xo,yo,zo);
- dir = dir - origin ;
+
+ dir.SetXYZ(dir.X()-(depthxOld-depthx)*TMath::Sin(dir.Phi()),
+ dir.Y()-(depthxOld-depthx)*TMath::Cos(dir.Phi()),
+ dir.Z()+depthzOld-depthz) ;
+
+ //Correct for the vertex position
+ dir = dir - fVtx ;
dir.SetMag(1.) ;
return dir ;
if(ts->GetEmcIndex()>=0)
emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
- AliPHOSCpvRecPoint * cpv = 0 ;
- if(ts->GetCpvIndex()>=0)
- cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
-
-// Int_t track = 0 ;
-// track = ts->GetTrackIndex() ; //TPC tracks ?
+// AliPHOSCpvRecPoint * cpv = 0 ;
+// if(ts->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 )) ;
// ###########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
sdp[AliPID::kPhoton][index] = 1. ;
//########## CPV-EMC Distance#######################
// Info("MakePID", "Distance");
- Float_t x = TMath::Abs(GetDistance(emc, cpv, "X")) ;
- Float_t z = GetDistance(emc, cpv, "Z") ;
+ Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
+ Float_t z = ts->GetCpvDistance("Z") ;
Double_t pcpv = 0 ;
Double_t pcpvneutral = 0. ;
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 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();
// 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 "<<effPur<<endl;
}
rp->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);
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() ;
}
TString message ;
message = "\nevent " ;
- message += gAlice->GetEvNumber() ;
+ message += gime->EventNumber();
message += " found " ;
message += recParticles->GetEntriesFast();
message += " RecParticles\n" ;
gime->WritePID("OVERWRITE");
}
}
-
-
+//____________________________________________________________________________
+void AliPHOSPIDv1::GetVertex(void)
+{ //extract vertex either using ESD or generator
+
+ //Try to extract vertex from data
+ if(fESD){
+ const AliESDVertex *esdVtx = fESD->GetVertex() ;
+ if(esdVtx){
+ fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
+ return ;
+ }
+ }
+
+ 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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff