X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSPIDv1.cxx;h=86b1598a3177754da67ba1c7311ad0342789800d;hb=9cef37509a0e12154a1b9ffd3938b5c4c2558ba9;hp=25079b84f357ef002466e509484694a605a49a3c;hpb=702ab87e845509d9379dca12e47c7e8a37eba7f6;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSPIDv1.cxx b/PHOS/AliPHOSPIDv1.cxx index 25079b84f35..86b1598a317 100644 --- a/PHOS/AliPHOSPIDv1.cxx +++ b/PHOS/AliPHOSPIDv1.cxx @@ -18,6 +18,45 @@ /* 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. + * */ //_________________________________________________________________________ @@ -87,6 +126,7 @@ // --- Standard library --- +#include #include "TFormula.h" #include "TBenchmark.h" #include "TPrincipal.h" @@ -95,15 +135,44 @@ // --- AliRoot header files --- //#include "AliLog.h" -#include "AliGenerator.h" #include "AliPHOS.h" #include "AliPHOSPIDv1.h" -#include "AliPHOSGetter.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSRecParticle.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 @@ -112,21 +181,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 ; } @@ -150,38 +271,14 @@ AliPHOSPIDv1::~AliPHOSPIDv1() delete fTFhhadronl; delete fDFmuon; } -//____________________________________________________________________________ -const TString AliPHOSPIDv1::BranchName() const -{ - - return GetName() ; -} -//____________________________________________________________________________ -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 @@ -382,12 +479,10 @@ void AliPHOSPIDv1::InitParameters() } //________________________________________________________________________ -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"); @@ -397,43 +492,26 @@ 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"); - AliInfo(Form("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(); } //________________________________________________________________________ @@ -512,22 +590,6 @@ Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const return param; } -//____________________________________________________________________________ -Float_t AliPHOSPIDv1::GetCalibratedEnergy(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. - - Float_t p[]={0.,0.,0.}; - for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); - Float_t enerec = p[0] + p[1]*e + p[2]*e*e; - return enerec ; - -} - //____________________________________________________________________________ Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const { @@ -659,49 +721,22 @@ Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString return par; } - - //____________________________________________________________________________ -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))) @@ -751,8 +786,7 @@ Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ TMath::Power(GetParameterPhotonBoundary(2),2)) + GetParameterPhotonBoundary(3); - AliDebug(1, Form("GetHardPhotonBit","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 photon else @@ -786,21 +820,47 @@ 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) ; - TMatrix dummy ; - - emc->GetGlobalPosition(dir, dummy) ; + 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"); + } + //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()) ; } - TVector3 origin(xo,yo,zo); - dir = dir - origin ; + + //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 ; @@ -886,17 +946,13 @@ void AliPHOSPIDv1::MakePID() const Int_t kSPECIES = AliPID::kSPECIESN ; - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; + Int_t nparticles = fRecParticles->GetEntriesFast() ; - Int_t nparticles = gime->RecParticles()->GetEntriesFast() ; - - TObjArray * emcRecPoints = gime->EmcRecPoints() ; - TObjArray * cpvRecPoints = gime->CpvRecPoints() ; - TClonesArray * trackSegments = gime->TrackSegments() ; - if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { + if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) { AliFatal("RecPoints or TrackSegments not found !") ; } - TIter next(trackSegments) ; + + TIter next(fTrackSegments) ; AliPHOSTrackSegment * ts ; Int_t index = 0 ; @@ -920,17 +976,17 @@ void AliPHOSPIDv1::MakePID() 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()) ; - -// 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 )) ; + AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ; } @@ -997,6 +1053,9 @@ void AliPHOSPIDv1::MakePID() // ###########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. ; @@ -1030,8 +1089,8 @@ void AliPHOSPIDv1::MakePID() //########## 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. ; @@ -1165,7 +1224,7 @@ void AliPHOSPIDv1::MakePID() for(index = 0 ; index < nparticles ; index ++) { - AliPHOSRecParticle * recpar = gime->RecParticle(index) ; + AliPHOSRecParticle * recpar = static_cast(fRecParticles->At(index)); //Conversion electron? @@ -1224,34 +1283,29 @@ void AliPHOSPIDv1::MakeRecParticles() { // Makes a RecParticle out of a TrackSegment - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - TObjArray * emcRecPoints = gime->EmcRecPoints() ; - TObjArray * cpvRecPoints = gime->CpvRecPoints() ; - TClonesArray * trackSegments = gime->TrackSegments() ; - if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { + 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()) ) { // cout<<">>>>>>>>>>>>>>>PCA Index "<At(index) ; + new( (*fRecParticles)[index] ) AliPHOSRecParticle() ; + rp = (AliPHOSRecParticle *)fRecParticles->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() ; @@ -1261,14 +1315,14 @@ void AliPHOSPIDv1::MakeRecParticles() // Choose the cluster energy range if (!emc) { - AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; + AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ; } Float_t e = emc->GetEnergy() ; - Float_t lambda[2] ; + Float_t lambda[2]={0.,0.} ; 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). @@ -1277,7 +1331,7 @@ void AliPHOSPIDv1::MakeRecParticles() 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(); @@ -1310,7 +1364,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 "<SetPIDBit(14) ; //Set momentum, energy and other parameters - Float_t encal = GetCalibratedEnergy(e); TVector3 dir = GetMomentumDirection(emc,cpv) ; - dir.SetMag(encal) ; - rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; + dir.SetMag(e) ; + rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ; 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 * ts1 = static_cast(fTrackSegments->At(rp->GetPHOSTSIndex())); + AliPHOSEmcRecPoint * erp = static_cast(fEMCRecPoints->At(ts1->GetEmcIndex())); TVector3 pos ; - geom->GetGlobal(erp, pos) ; + fGeom->GetGlobalPHOS(erp, pos) ; rp->SetPos(pos); index++ ; } @@ -1377,7 +1429,7 @@ void AliPHOSPIDv1::Print(const Option_t *) 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() ; } @@ -1388,23 +1440,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 += " " ; @@ -1452,7 +1498,7 @@ 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) ; + fParameters = new TMatrixF(16,4) ; const Int_t kMaxLeng=255; char string[kMaxLeng]; @@ -1471,7 +1517,7 @@ void AliPHOSPIDv1::SetParameters() &(*fParameters)(i,0), &(*fParameters)(i,1), &(*fParameters)(i,2), &(*fParameters)(i,3)); i++; - AliDebug(1, Form("SetParameters", "line %d: %s",i,string)); + AliDebug(1, Form("Line %d: %s",i,string)); } fclose(fd); } @@ -1565,40 +1611,24 @@ void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString par } //____________________________________________________________________________ -void AliPHOSPIDv1::Unload() -{ - //Unloads RecPoints, Tracks and RecParticles - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - gime->PhosLoader()->UnloadRecPoints() ; - gime->PhosLoader()->UnloadTracks() ; - gime->PhosLoader()->UnloadRecParticles() ; -} - -//____________________________________________________________________________ -void AliPHOSPIDv1::WriteRecParticles() -{ - //It writes reconstructed particles and pid to file - - 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"); +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 && 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